André Jacomel Torii

Rafael Holdorf Lopez

Fabio Biondini

André ; Torii

André Jacomel Torii

An approach

Keywords: RBDO, truss structures, geometry optimization, topology optimization reliability based optimization, truss structures, geometry optimization, topology optimization

come

Introduction

Methods for the optimization of trusses, where the members areas are taken as design variable, are well established and there is a rich literature on this subject (Hemp, 1973;[START_REF] Pedersen | On the minimum mass layout of trusses[END_REF][START_REF] Pedersen | Topology optimization of three-dimensional trusses[END_REF]. The problem of truss shape optimization, where the nodes' positions are taken as design variables, is also addressed in literature [START_REF] Achtziger | Simultaneous optimization of truss topology and geometry, revisited[END_REF][START_REF] Achtziger | On simultaneous optimization of truss geometry and topology[END_REF][START_REF] Kocvara | How mathematics can help in design of mechanical structures[END_REF], but not as extensively.

One of the main reasons why sizing optimization was studied more frequently than shape optimization, also called here geometry optimization, is that the sizing optimization problem can, in some cases, be stated as a linear programming problem.

For such cases, very efficient methods are available and it is possible to guarantee certain important mathematical properties of the solution, such as the existence of a global optimum [START_REF] Nocedal | Numerical optimization[END_REF][START_REF] Rao | Computational methods in optimization considering uncertainties -an overview[END_REF][START_REF] Arora | Introduction to optimum design[END_REF]. The truss geometry optimization problem is, instead, nonlinear by its nature, and therefore it needs to be solved by nonlinear optimization methods, which are in general more complex and computationally demanding than the linear programming techniques. Also, the geometry optimization of a truss may lead to a non-convex optimization problem, and several local minima may exist [START_REF] Achtziger | Simultaneous optimization of truss topology and geometry, revisited[END_REF][START_REF] Achtziger | On simultaneous optimization of truss geometry and topology[END_REF]).

In the last decades the increase in computational resources allowed an increase in the use of nonlinear optimization methods, and thus the problem of geometry optimization was studied more extensively. Some strategies have been developed to deal with such optimization problem, and with respect to [START_REF] Achtziger | On simultaneous optimization of truss geometry and topology[END_REF] these strategies can be grouped as: simultaneous optimization of truss topology and geometry, alternating optimization and implicit programming optimization.

In deterministic optimization, however, the uncertainties involved in the design problem, such as those affecting material properties and loads, among others, are not considered. Robust optimization or reliability based design optimization (RBDO) methods are usually employed to take such uncertainties into account (see, for instance, [START_REF] Beyer | Robust optimization -a comprehensive review[END_REF], Schuëller and Jensen (2009)). The former has as main goal the minimization of the variability of some parameters related to system response due to its uncertainties. For example, [START_REF] Calafiore | Optimization under uncertainty with applications to design of truss structures[END_REF] applied this concept in the field of design of truss structures.

The main goal of the RBDO is to optimize a structure and guarantee that its probability of failure is lower than a certain level chosen a priori by the designer. Nakib (1991), Thampan and Krishnamoorty (2001), [START_REF] Togan | Reliability and reliability-based design optimization[END_REF] dealt with the RBDO of truss structures, regarding only the size optimization of such However, only a few articles have dealt with the reliability based shape and topology optimization of truss structures (Morutsu and Shao, 1990;[START_REF] Stocki | Study on discrete optimization techniques in reliability-based optimization of truss structures[END_REF]. Thus, this article proposes a probabilistic approach for solving the simultaneous optimization of geometry and topology for statically indeterminate trusses, taking into account the uncertainties on the applied forces as well as the yielding stresses. In this approach the applied forces and the yielding stresses are modeled as random variables, and the failure constraints are expressed in probabilistic terms.

This study has three main contributions: (i) It is shown how to efficiently pursue the sensitivity analysis required by the gradient based method used in the optimization process; (ii) Based on the assumptions of linear structural behavior and independent and normally distributed random variables, the RBDO problem is posed in such a way that the reliability of the structure is accessed directly, without using iterative methods such as a First Order Reliability Method (FORM) or Monte Carlo Simulation (MCS) [START_REF] Haldar | Reliability assessment using stochastic finite element analysis[END_REF].; (iii) It is shown that, when uncertainties are considered, the resulting optimum structure may not be a simply scaled version of the deterministic solution (e.g. higher cross-section areas), but there may be a change in the structural geometry as well.

Deterministic optimization

Deterministic optimization of geometry and topology of trusses is presented in the following, with emphasis on how the sensitivities may be efficiently evaluated.

Formulation of the optimization problem

The optimization problem is posed as the minimization of the volume of the structure subject to stress constraints by taking the nodal coordinates and cross-section areas as design variables:

Find: x and A that gives

) ( . ) , ( min x L A A x T V = , (1)
subject to

,...,m) , (j g t j j 2 1 0 = ≤ - + = σ σ , (2) ,...,m) , (j g c j m j 2 1 0 = ≤ + - = + σ σ . ( 3
)
where V is the volume of the structure, x is the vector of nodal coordinates, A is the vector of member areas, L is the vector of member lengths, g j are stress constraints, σ j is the stress on member j, σ t is the yielding stress in tension, σ c is the yielding stress in compression and m is the number of members subjected to stress constraints. In this article, buckling constraints are not introduced.

For convenience, the design variables A and x can be grouped into a single design vector X, and the constraints from Eq. ( 2) and Eq. ( 3) can be grouped into a single vector of constraints g. In this way, the previous problem is rewritten as follows:

Find: X that gives

L A X . ) ( min T V = , (4)
subject to where g is a vector with 2m components since there are two constraints defined for each bar of the structure.

0 g ≤ , (5)

Sensitivity analysis

Most nonlinear optimization methods need to evaluate the gradients of both the objective function and constraints with respect to the design variables. The partial derivative of the volume V with respect to each member area is given by the length of each bar. The partial derivative of V with respect to the nodal coordinates is also easy to obtain, and expressions similar to the ones presented by [START_REF] Martínez | Growth method for size, topology, and geometry optimization of truss structures[END_REF] and [START_REF] Torii | A simple geometry optimization method for statically indeterminate trusses[END_REF] can be used for this purpose. Moreover, all these derivatives can be efficiently evaluated by finite differences, if necessary, since the evaluation of V is not computationally demanding.

The major difficulty arises when evaluating the gradients of the constraints from Eq. ( 2) and Eq. (3), since these constraints involve stresses that depend implicitly on nodal displacements. From Eq. (2) and Eq. (3) it can be seen that the partial derivatives from these expressions will be the same with opposite sign and, consequently, only one of them must be evaluated numerically.

Working with an arbitrary chosen constraint g j , as defined in Eq. ( 2), its partial derivative with respect to a design variable X i can be written as

i T j i j i j X dX d X g ∂ ∂ ∂ ∂ + = ∂ ∂ u u . σ σ , ( 6
)
where the symbol d stands for an ordinary derivative and the symbol ∂ stands for a partial derivative. It is worth noting that the first term in Eq. ( 6), the derivative dσ j /dX i , is a derivative related to the direct dependence of σ j on X i . The second term, instead, is an implicit derivative related to the dependence of σ j on X i through the nodal displacements u.

The partial derivative ∂σ j /∂u can be evaluated analytically, or found by finite differences if necessary, since the evaluation of σ j when only displacements are changed is straightforward. However, the partial derivative ∂u/∂X i is difficult to evaluate, and its computation should be avoided. This can be accomplished by using the adjoint method, as described by [START_REF] Haftka | Elements of structural optimization[END_REF].

Adjoint method

The following system of linear equations relating the applied forces F and the displacements u is considered (Bathe 1996):

F K.u = , ( 7
)
where K is the stiffness matrix of the structure. The partial derivative of Eq. ( 7) with respect to a design variable is

0 F u K. .u K = ∂ ∂ - ∂ ∂ + ∂ ∂ i i i X X X , (8)
which has been rearranged in order to give an expression equal to a zero vector.

Since Eq. ( 8) is equal to a zero vector, Eq. ( 6) can be rewritten as where τ j is an arbitrary vector. Eq. ( 10) shows that the evaluation of ∂u/∂X i can be avoided if

        ∂ ∂ - ∂ ∂ + ∂ ∂ - ∂ ∂ ∂ ∂ + = ∂ ∂ i i i T j i T j i j i j X X X X dX d X g F u K. .u K τ u u . . σ σ (9) or         ∂ ∂ - ∂ ∂ - ∂ ∂         - ∂ ∂ + = ∂ ∂ i i T j i T j T j i j i j X X X dX d X g F .u K τ u K τ u . . . σ σ , (10)
0 K τ u = - ∂ ∂ . T j T j σ . ( 11
)
Thus, the vector τ j is chosen as the solution of the system of linear equations

u τ K ∂ ∂ = j j σ . , (12)
where the symmetry of the stiffness matrix has been used (Bathe 1996). Equation ( 10)

then becomes         ∂ ∂ - ∂ ∂ - = ∂ ∂ i i T j i j i j X X dX d X g F .u K τ . σ , (13)
that can be easily evaluated once the vector τ j is found from Eq. ( 12). In fact, the derivatives ∂K/∂X i , and ∂F/∂X i can be evaluated analytically [START_REF] Torii | A simple geometry optimization method for statically indeterminate trusses[END_REF], or numerically by finite differences.

It is worth noting that the evaluation of the partial derivative ∂u/∂X i could be also avoided by using Eq. ( 8) rewritten as follows:

        ∂ ∂ - ∂ ∂ = ∂ ∂ - .u K F K u i i i X X X . 1 . ( 14
)
However, in this approach the number of systems of linear equations to be solved is equal to the number of design variables. In the adjoint method, instead, the number of systems of linear equations is equal to the number of bars in the ground structure. Since the number of design variables is the number of members (cross section areas) plus the number of nodal coordinates taken as design variables, the adjoint method will always lead to a lower number of systems of linear equations to be solved.

Finally, it is outlined that if the constraints and the objective function of the optimization problem are both incorporated into a single objective function, by means of the augmented Lagrangean method or some other penalty methods [START_REF] Rao | Computational methods in optimization considering uncertainties -an overview[END_REF], the number of systems of linear equations to be solved during sensitivity analysis can be drastically reduced, as demonstrated by [START_REF] Pereira | Topology Optimization of continuum structures with material failure constraints[END_REF].

Alternative loading conditions

When a set of s alternative loading conditions { } 4) and Eq. ( 5) becomes Find: X that gives

L A X . ) ( min T V = , ( 15
) subject to ) ,..., 1 ( s k k = ≤ 0 g , ( 16
)
where the vectors g k from Eq. ( 16) are defined for each loading condition in Eq. ( 2)

and Eq. ( 3).

Working again with constraints related to the yielding stress in compression, since constraints related to tension will just have opposite sign, the derivative of constraint g jk , related to the stress in the bar j when the loading condition k is applied, is

F o r P e e r R e v i e w O n l y i k k T jk i jk i jk X dX d X g ∂ ∂ ∂ ∂ + = ∂ ∂ u u . σ σ . ( 17
)
The adjoint vector is now given by the solution of the system of linear

equations k jk jk u τ K ∂ ∂ = σ . ( 18
)
and Eq. ( 17) can be rewritten as

        ∂ ∂ - ∂ ∂ - = ∂ ∂ i k k i T jk i jk i jk X X dX d X g F .u K τ . σ . ( 19
)
This formulation leads to an increase in the number of systems of linear equations that must be solved during the sensitivity analysis, i.e. s times the number of systems of linear equations needed for one loading condition. However, it is noted that the stiffness matrix remains the same. Consequently, decomposition techniques can be employed efficiently for the sensitivity analysis, since it may be assumed that the main computational effort lies in the decomposition procedure itself (Bathe 1996).

Bounds on the design variables

There are two different approaches for defining bounds on the nodal coordinates. In the first approach, bounds can be defined locally for each design variable, as shown in Fig. 1a. In this case, there are different bounds for each design variable, and this can be accomplished by defining a rectangular feasible region around each node. The second approach is that of defining bounds for all the design variables at once, as shown in Fig. 1b. In this case, the bounds are the same for all design variables. These two approaches may lead to different results, since the feasible domain defined in the first approach is smaller. However, the first approach may prevent problems related to node superposition, if the bounds are defined properly. Therefore, this approach may be recommended when there are many nodes for which the coordinates are taken as design variables in the optimization procedure.

For members areas it is necessary to allow only positive values, defining in this way a lower bound for these design variables. Upper bounds on members areas are not strictly necessary, since the algorithm will seek a structure with minimum volume. However, defining upper bounds on all design variables may be recommended, since the optimization problem becomes better posed.

Probabilistic optimization

In the following, the deterministic optimization problem presented previously is formulated in probabilistic terms.

Formulation of the problem and reliability assessment

Consider the applied forces F and the yielding stresses

L A X . ) ( min T V = , ( 20
) subject to ( ) ( ) ) ,..., 1 ( 0 0 0 m j g P P G j j = ≤ ≤ - = Ξ , (21)
where P(•) is the probability of the constraint to be fulfilled and P 0 is a minimum probability level of the constraint or its reliability. Note that now the constraints of the optimization problem are affected by the random variable vector Ξ , becoming themselves random variables. The constraint from Eq. ( 21) states that the probability of ( )

Ξ j g
being respected must be higher than a minimum probability P 0 .

The random variable ( )

Ξ j g
is constructed based on a linear structural behavior, as occurs in most structural optimization procedures (Hemp, 1973;[START_REF] Pedersen | On the minimum mass layout of trusses[END_REF]Pedersen, 1973;[START_REF] Achtziger | Simultaneous optimization of truss topology and geometry, revisited[END_REF][START_REF] Achtziger | On simultaneous optimization of truss geometry and topology[END_REF][START_REF] Martínez | Growth method for size, topology, and geometry optimization of truss structures[END_REF][START_REF] Torii | A simple geometry optimization method for statically indeterminate trusses[END_REF][START_REF] Pereira | Topology Optimization of continuum structures with material failure constraints[END_REF], and by assuming F and σ independent normal random variables.

If for a given applied force vector F 0 , the resulting stress in a given member is σ 0 , then for an arbitrary applied force vector F, obtained by the multiplication of F 0 by a scalar, the stress in that member is (by the principle of superposition from structural mechanics)

F F . . 0 0 0 0 σ σ σ = = F F , ( 22
)
where ⋅ denotes the norm of a vector, and F and F 0 are the norms of F and F 0 , respectively. Equation ( 22) can be substituted into Eq. (2), giving ( )

0 . 0 0 ≤ - = t F F g σ σ Ξ , ( 23
)
where the index j has been dropped for convenience. Since both the applied forces and yielding stresses are considered as independent normal random variables, a linear combination of them is also a normal random variable. Consequently, by denoting µ 1 and s 1 the mean value and standard deviation of the applied force, respectively, and µ 2 and s 2 the mean value and standard deviation of the yielding stress, respectively, ( )

Ξ g
as given by Eq. ( 23) is a normal random variable with mean

( ) 2 1 0 0 . µ µ σ µ - = Ξ F g (24)
and standard deviation

( ) 2 2 2 1 2 0 0 2 . s s F s g +         = Ξ σ . ( 25
)
The normalized value of ( )

Ξ g is then: ( ) ( ) ( ) Ξ Ξ - = g g s g µ β Ξ , ( 26
)
where β is the reliability or also called reliability index of the constraint (Lemaire et al., 2005). Thus, the reliability index is related to the probability of the constraint being feasible by:

( ) ( ) [ ] ( ) ( ) [ ] β β Φ = ≤ ≤ Φ = - 0 or 0 1 Ξ Ξ j j g P g P , ( 27
)
where Φ = Φ[•] is the standard normal cumulative probability function. The required minimum reliability level of the structure P 0 is related to the so-called target reliability Substituting Eq. ( 27) and ( 28), the constraint given by Eq. ( 21) becomes

) ,..., 1 ( 0 m j G j t j = ≤ - = β β . ( 29
)
The index j emphasizes that a reliability index is computed for each bar.

Sensitivity analysis

The derivative of Eq. ( 29) with respect to a design variable is given by

i T j i j i j X dX d X G ∂ ∂ ∂ ∂ - - = ∂ ∂ u u . β β , ( 30
)
where the derivative of the constraint with respect to a given design variable is composed of an explicit plus an implicit derivative, as occurs for the deterministic case. The ordinary derivative dβ j /dX i can be evaluated by finite differences, since only perturbations to the design variables are necessary. The second part of Eq. ( 30) can be solved efficiently by applying the adjoint method again, in order to avoid the evaluation of the derivative ∂u/∂X i . Following the same procedure described for the deterministic approach:

        ∂ ∂ - ∂ ∂ + - = ∂ ∂ i i T j i j i j X X dX d X G F .u K τ . β , ( 31
)
where the adjoint vector τ j is found by solving

u τ K ∂ ∂ - = j j β . . ( 32
)
The derivative ∂β j /∂u can be efficiently evaluated by finite differences. In order to access the reliability of the structure, it is necessary to study the effect of each force, 1 F and 2 F , separately. Thus, considering each force separately (again by the superposition principle)

Loading conditions given by several forces with different standard deviations

1 1 0 1 0 1 1 0 1 0 1 . . F F σ σ σ = = F F (33) and 2 2 0 2 0 2 2 0 2 0 2 . . F F σ σ σ = = F F , (34)
where the same notations as of Eq. ( 22) hold, but here the superscripts 1 and 2 represent the quantities for forces 1 F and 2 F . That is, one structural analysis is made for each force 1 0 F and 2 0 F , which gives the stresses 1 0 σ and 2 0 σ , respectively (usually, these forces are taken as unit forces, for convenience). The stress caused by the

application of 1 F + 2 F is then 2 2 0 2 0 1 1 0 1 0 2 1 . . F F F F σ σ σ σ σ + = + = , ( 35
)
which can be rewritten, for convenience, as The important aspect here is that the effect inside the bar can be written as a linear combination of the separate effects of each force, as this is a linear structural analysis. This is precisely the superposition principle from structural mechanics.

2 2 1 1 . . F k F k + = σ . ( 36
Since the applied forces are Gaussian random variables, the mean value of the stress in the bar and its standard deviation due to the combination of forces 1 F and 2 F (Eq. ( 36)) are, respectively,

2 1 . . 2 1 F F k k µ µ µ σ + = (37)
and

( ) ( ) ( ) 2 2 2 1 2 2 1 . . F F s k s k s + = σ , ( 38
)
where 1 It is possible to evaluate Eq. ( 39) and Eq. ( 40) with Eq. ( 37) and Eq. ( 38), and consequently, to obtain the reliability index as described previously. Also, note that Eq. ( 39) and Eq. ( 40) are analogous to Eq. ( 24) and Eq. ( 25), respectively.

The previous results can also be extended to more than two forces acting at the same time. However, this leads to a significant increase in the number of structural analysis that must be made. Also, the sensitivity analysis is also changed, since now there will be one adjoint vector associated with each force.

Numerical examples

In this section, several numerical examples are solved in order to discuss the main aspects of the approach presented. Note that for all the following examples each figure has its own scale, and the reader can compare the different solutions by the volume of material V presented for each example. A lower bound for the areas is defined as 1E-10m 2 (i.e. 0.1mm 2 ), in order to avoid singularity of the stiffness matrix.

Finally, the optimization algorithm used here is the Sequential Quadratic Programming (SQP) [START_REF] Nocedal | Numerical optimization[END_REF][START_REF] Rao | Computational methods in optimization considering uncertainties -an overview[END_REF].

Example 1 -two bar structure

The first example discussed here is the optimization of a two bar structure. From the results presented in Fig. 3 it can be seen that changing only the standard deviation of the yielding stresses leads to changes in the topology and geometry of the optimum structure. The structure in Fig. 3c is symmetric, since the yielding stresses both in tension and compression have the same mean value, in magnitude, and the same standard deviation. The solution of the deterministic problem, from Fig. 3b, is also symmetric. However, the structure from Fig. 3d is not symmetric because of the different values of the standard deviation of the yielding stresses. Finally, note that increasing the standard deviation of the yielding stress leads to an increase in the volume of the structure, as expected.

Example 2 -Effect of the reliability index

The ground structure from Fig. 4a For the three different cases the optimum topology and the geometry remain the same. The only difference is that all bar areas are scaled, and the volume of material increases for increasing reliability indexes. From this example, one may conclude that increasing the reliability index will always lead to a scaled solution, but this is not true for all cases, as demonstrated by the next example.

Example 3 -Effect of the reliability index

In order to show that increasing the reliability index may lead to geometry and topology changes, return to the example from Fig.

Example 4 -Effect of limiting the maximum area of the bars

The ground structure of this example is shown in Fig. 6, and this example demonstrates the effect of defining a maximum value for the bar areas. The structure lengths are L x = 2m and L y = 1m. The material properties are the same as the ones defined for the example 1. The structure is subjected to one loading condition with mean value F = 1,000N and standard deviation σ F = 100N. The nodes are taken as fixed during the optimization procedure, and thus only sizing optimization is needed for this case. Additionally, for this problem the bar areas are not allowed to be bigger than 1E-8m 2 (i.e. 10mm 2 ). Three cases are studied, with β equal to 0 (deterministic problem), 3.1 and 4.75.

The results from Fig. 6 show that defining a maximum bar area may lead to different topologies when only the reliability index is changed. This happens since the maximum bar area may enforce the optimization algorithm to "use" more bars in order to increase the reliability of the structure, leading to different topologies.

Besides, it is reasonable to assume that different topologies will lead to different geometries if geometry optimization is pursued.

Example 5 -Loading conditions composed of several forces

Figure 7 presents a ground structure that is subjected to three loading conditions.

However, vertical forces and horizontal forces have different mean values and standard deviations. This example illustrates the application of the method when a loading condition is composed by forces with different mean values and standard deviations. The results for this example are presented in Fig. 8. As a consequence of the increase to the lateral load, for the same vertical load, there was an increase of the volume of the structure and a change of its geometry and topology. The structures designed for higher lateral loads are clearly more fitted to resisting the increased bending moment that develops in these cases.

For the same ground structure from Fig. 7, another interesting example can be conceived. Taking the vertical force F 1 as defined previously and the mean of the horizontal force F 2 = 250N, the same problem is now solved for three different standard deviations: σ F2 = 25N, σ F2 = 50N and σ F2 = 75N. The results are presented in Fig. 9, where it can be noted that increasing the standard deviation of the lateral load leads to an increase in the volume of the structure. Also, the geometry and the topology change in order to resist the increased bending moments that appear in this case. The example from Fig. 10 compares results from the simultaneous optimization of geometry and sizing with results given by sizing optimization only. The ground structure has a total height of 1m and a total length of 8m. The material properties are as defined for the example 1. The target reliability index of the structure is β =4.75, and all the nodes from the upper chord are allowed to be moved up and down, by the optimization algorithm, to positions as far as 0.5m from its original position. The applied force is F = 1,000N with a standard deviation σ F = 100N.

Example 6 -Comparison of the sizing and simultaneous sizing-geometry optimization

The results from Fig. 10 demonstrate that changes in the geometry of the structure can lead to significant improvements of the solution, even if the topology (the bars that comprise the structure) remains basically the same. This justifies the effort in including the nodal coordinates as design variables of the optimization problem.

Example 7 -Local minima

The last example is used to illustrate the existence of local minima that can occur for problems of truss optimization. For this task, the ground structure from the example of Fig. 6a is used, with β equal to 4.75, but the maximum cross section area is taken as 5.0E-7m 2 (i.e. 500mm 2 ). In this case the upper bound for the cross section area does not affect the results. Also, the optimization algorithm is now allowed to move the nodes in the upper chord (except the nodes of support) up and down to positions as far as 0.5m from its original positions. This optimization problem is solved several times with random initial solutions. Two different local optima found for this problem are presented in Fig. 11. Note that the existence of local optima is also common for the deterministic case [START_REF] Achtziger | Simultaneous optimization of truss topology and geometry, revisited[END_REF][START_REF] Achtziger | On simultaneous optimization of truss geometry and topology[END_REF][START_REF] Kocvara | How mathematics can help in design of mechanical structures[END_REF][START_REF] Arora | Introduction to optimum design[END_REF]. For the examples presented earlier, the initial bar areas were all taken equal to 1E-9m 2 (i.e. 1mm 2 ), since it was observed that starting with solutions for which the constraints were not respected gave, in general, better results. This can be explained by the following reasoning. If the initial solution is over-designed, then the algorithm will try to reduce the member lengths in order to reduce the volume of the structure, and will not "care" about stress constraints since the current solution is over-designed.

This may bring the solution to a geometry composed of shorter members, but which may not be the best geometry possible. However, when starting at solutions for which the constraints are violated, the algorithm will try to rearrange the geometry of the structure in order to reduce the stresses in the bars. This will probably bring the solution to a "good" geometry that was obtained with the intent of reducing the stresses inside the structure, instead of just trying to obtain shorter bars.

Further Comments

Note that the optimization problems presented in this article have several reliability constraints, and consequently several reliability analysis problems must be solved in order to evaluate a given design point of the optimization process. Here, all these reliability analysis problems are solved analytically, and consequently the main computational effort lies in finding the stresses inside the structures for each applied load, i.e. one needs to solve one system of linear equations for each applied load.

If the reliability analysis problems were solved using a traditional approach (e.g.

FORM, Monte Carlo Simulation), an entire reliability analysis problem would need to be solved for each constraint. It is well known that solving a reliability analysis problem using these approaches leads to the solution of several systems of linear equations (since each constraint must be evaluated several times to obtain each reliability index), which leads to an increased computational effort. Take for example the structure from Fig. 7. In this case, there are 21 bars and 3 loading conditions.

Since there are 2 constraints for each bar for each loading condition, the optimization problem has a total of 21x2x3 = 126 reliability constraints. In a traditional approach (e.g. FORM), one would need to solve 126 reliability analysis problems in order to evaluate a given design point, and each reliability analysis problem would need to solve several systems of linear equations.

However, if the approach presented in this article is employed, one only needs to find the stresses inside the structure for each of the 10 loads of this problem (2 from loading condition 1, 4 from loading condition 2 and 4 for loading condition 3), i.e.

one needs to solve 10 systems of linear equations in order to evaluate a given design point and the remaining information (the reliability index of each constraint) is evaluated analytically. Thus, the computational effort involved in the present approach is orders of magnitude smaller than if a FORM algorithm were used.

Regarding stability, the extension of the proposed approach to consider

Eulerian buckling is straightforward. The Eulerian buckling limit stress depends only on the moment of inertia, the length and the Young's Modulus of each bar and consequently does not depend on the random variables considered in this article (yielding stress and applied forces). That is, in a Eulerian buckling limit state function, only the buckling stress on the bar would be a random variable, different from the yielding stress case (Eq.( 3)) where both the stress on the bar and the yielding stress are random variables. Thus, the inclusion of a buckling constraint would be made by the same procedure as for a standard compression stress constraint, but by taking the allowable buckling stress as a deterministic value.

Conclusions

This article presented a formulation for the simultaneous optimization of topology and geometry of truss structures. General aspects such as constraints and several loading conditions were also discussed, together with a detailed discussion on sensitivity analysis.

The deterministic optimization scheme was then extended to the case when the yielding stresses and the applied forces are Gaussian random variables. The important aspect here is that considering these variables as Gaussian allows one to access the reliability of the structure directly, without using iterative methods such as a First

Order Reliability Method (FORM) or Monte Carlo Simulation (MCS). Also, the sensitivity analysis is also simplified by the use of the adjoint method. Thus, the probabilistic problem as defined here has a similar computational cost per optimization step as the deterministic problem, and that is an important aspect.

The main drawbacks of the proposed approach are: (i) not considering the bar areas and the Young Modulus as random variables; (ii) modeling the applied forces and the yielding stresses as Gaussian random variables. However, this decision is based on three main aspects.

First, the resulting optimization problem is simplified and the computational effort greatly reduced. Reliability Based Optimization problems often lead to a significant increase in computational effort, and may be even unaffordable from the computational point of view. Consequently, to reduce the computation effort is a main objective of this article.

Second, the cross-section areas are, in most cases, to be chosen from an array of discrete values. Consequently, the bar areas are expected to be rounded up to the nearest standard area available in the market. In this context, considering the bar areas as random variables would lose all sense, since the final areas would probably present significant modifications, in comparison to the standard deviation of this variable. In this context, it may be more reasonable to consider the areas as random variables in a subsequent design stage.

Third, taking the Young Modulus as a random variable is not expected to lead to significant changes to the problem, since the influence of this parameter is not so direct to the forces inside the structure. Changes to the Young Modulus may lead to a redistribution of forces inside the structure, but these changes are not expected to be significant, for most cases, unless the Young Modulus is drastically changed. The influence of the yielding stresses and applied forces, instead, is direct.

The numerical examples presented allowed three main conclusions about the reliability based optimization problem. First, changes to probabilistic parameters, such as standard deviations of yielding stresses, are expected to lead to changes to both optimum topologies and geometries. Second, some optimum structures are just scaled when the reliability index is changed, but this is not true for all cases. Finally, the problem of simultaneous optimization of geometry and topology can have more than one local minimum. However, this is not an exclusive feature of the probabilistic problem, since this is also observed in deterministic problems.

  et al. (2007) pursued the reliability based robust optimization of trusses grouping at the same time the goals of the robust and RBDO.

  problem from Eq. (

  distribution. For convenience of notation, they are grouped into the random variable vector Ξ . The optimization problem still searches for the minimum volume structure, but now subject to a minimum reliability level of the structure. A component level reliability constraint is considered, instead of dealing with the probabilistic failure constraint at the system level. In other words, a minimum reliability level is enforced for each bar. For a single loading condition this leads to the following problem:

2 F

 2 In order to describe how forces with different standard deviations can be taken into account into the same loading condition, consider the example from Fig.2. In this case, the structure is subjected to a single loading condition. That is, the forces 1 F and are applied at the same time, but have different mean values and standard deviations.

F µ and 1 Fs

 1 are the mean value and standard deviation of the force F 1 , respectively (the same notation is also used for the force F 2 ). Moreover, σ µ and σ s are the mean value and the standard deviation of the stress inside the bar, respectively. Note that the stress is a Gaussian random variable, as it is the linear combination of two Gaussian random variables, namely the two forces.Thus far it has been shown how to obtain the mean value and the standard deviation of the stress inside the bar, given two forces with magnitude F 1 and F 2 acting simultaneously. Consider now the stress constraint from Eq. (2). Since it is a linear combination of Gaussian random variables (the stress and the maximum allowable stress), the mean and the standard deviation of the maximum allowable stress in tension.

  The ground structure from Fig.3ais subjected to a single loading condition with mean F = 10,000N and standard deviation σ F = 1,000N. The lengths of the structure are L x = 2.5m and L y = 5m and the Young Modulus of the material is E = 200GPa. The minimum reliability index is taken as β = 3.1. The yielding stress in tension has mean σ t = +250MPa and the yielding stress in compression has mean σ c = -250MPa. For the case of Fig.3cthe yielding stresses in tension and in compression have standard deviation σ σ = 10MPa; while for the case of Fig.3dthe standard deviation of the yielding stress in tension is raised to σ σt = 50MPa. The node of the applied forces is allowed to be moved up and down by the optimization algorithm.

  has lengths L x = 9m and L y = 1m. The material properties are the same as the ones defined for the example 1. The structure is subjected to 8 loading conditions with mean F = 1,000N and standard deviation σ F = 100N. Three cases are studied, the deterministic problem, the problem for β equal to 3.1 and for β equal to 4.75. The nodes of the upper chord are allowed to be move up the optimization algorithm, to positions as far as 0.9m from the original position.

  3. The loads and dimensions of the ground structure remain unchanged as well as most of the material properties. The only modification is that the standard deviation of the yielding stress in compression is now equal to σ σc = 10MPa and of the yielding stress in tension equal to σ σt = 20MPa. The problem is then solved for reliability indexes equal to 0 (deterministic problem), 2, 4, 6, 8 and 10. The results are presented in Fig.5. It can be seen that increasing the reliability index leads, in this case, to changes in both topology and geometry. This example contradicts the assumption of only scaling the cross-section areas when different reliability indexes are used. The reason why the example from Fig.4gave scaled solutions for different reliability indexes lie in the fact that all the properties of that ground structure are symmetric. The mean values and standard deviation of the yielding stresses are equal in magnitude, together with the applied forces. The case from Fig.5, instead, has different standard deviations for the yielding stresses, and this leads to changes in topology and geometry when the reliability index is changed.

  The ground structure from Fig.7has a total height of 4m and a total width of 2m. The material properties are the same as from example 1. The target reliability index is 3.1, and all nodes (except the nodes of the supports) are allowed to be moved left and right by the optimization algorithm, to positions as far as 0.8m from its original position. The vertical force mean value is F 1 = 1,000N and its standard deviation is σ F1 = 100N. For the horizontal forces F 2 , three cases are studied: mean equal F 2 = 50N and standard deviation σ F2 = 5N; mean F 2 = 150N and standard deviation σ F2 = 15N; and mean F 2 = 250N and standard deviation σ F2 = 25N. Finally, symmetry of the geometry is enforced.

Figure 1 :Figure 2 :

 12 Figure 1: Bounds on nodal coordinates defined a) locally for each node, and b) globally for all nodes.

Figure 4

 4 Figure 4: a) Ground structure and optimum solutions b) deterministic problem, c) β equal to 3.1 and d) β equal to 4.75.

Figure 10 Figure 11 :

 1011 Figure 10: a) Ground structure, b) optimum topology for the deterministic problem, c) optimum geometry for the deterministic problem, d) optimum topology for the probabilistic problem, e) optimum geometry for the probabilistic problem.

URL: http:/mc.manuscriptcentral.com/geno Email: A.B.Templeman@liverpool.ac.uk Engineering Optimization


Franck Celestini

Thomas Frisch

Yves Pomeau

Take-off of small Leidenfrost droplets

Keywords: numbers: 47, 55D-,68, 03, -g

We put in evidence the unexpected behaviour of Leidenfrost droplets at the later stage of their evaporation. We predict and observe that, below a critical size R l , the droplets spontaneously takeoff due to the breakdown of the lubrication regime. We establish the theoretical relation between the droplet radius and its elevation. We predict that the vapour layer thickness increases when the droplets become smaller. A satisfactory agreement is found between the model and the experimental results performed on droplets of water and of ethanol.

The name of J.G. Leidenfrost (1715Leidenfrost ( -1794) ) is still remembered because he was the first to publish the observation that a puddle of water dropped on a very hot surface divides into droplets, which stand each above the surface and slowly evaporate. At the end of the process, he noticed [START_REF] Leidenfrost | De Aquae Communis Nonnullis Qualitatibus Tractatus[END_REF] : "... before the whole drop disappears. Which at last exceedingly diminished so that it can hardly any more be seen, with an audible crack, which with the ears one easily hears, it finishes its existence, and in the spoon [making the hot surface] it leaves a small particle of earth ..."

Over the years, the Leidenfrost phenomenon has attracted the attention of many investigators, see for instance [3][4][5][6][7][8][9], but seemingly little interest has been paid to the observation by Leidenfrost of the final stage of the droplet evaporation, in particular what happens when the droplets become very small before finally disappearing. We investigate in this letter theoretically and experimentally the final regime of the Leidenfrost droplets. We show that when the droplets become smaller than a well defined radius, they suddenly take-off due to the breakdown of the classical Leidenfrost lubrication regime. They thus reach an elevation which is much higher than their radius.

Let us consider first an almost spherical drop of volatile liquid with radius R and mass M standing slightly above a hot flat surface with a vapour film of thickness h in between. Let l be the horizontal width of this film, this length will be related to the other physical parameters of the problem.

Our estimate will rely on the fact that, in a range of parameters, the horizontal length scale l is much larger than h. Both l and h are less than R, the radius of the droplet. Furthermore, we shall consider a range of parameters where the vapour pressure in the film is negligible compared to Laplace's capillary pressure. Therefore the droplet remains almost exactly spherical in this regime. This will yield the condition (3) derived below, valid under the constraint h ≪ R. This condition hap- pens to be impossible to satisfy if the droplet becomes too small, that is if its weight is too small to balance the upward pressure in the gap between the hot plate and the droplet. Therefore, it is natural to guess that, when this happens, the droplet lifts off the plate. This defines a second regime where h ≫ R. The scaling laws for this second regime are given in this letter and we show that the droplet reaches an higher elevation as it gets smaller and smaller by evaporating. Energy conservation during the evaporation process of the droplet (namely Stefan's boundary condition on the liquid/vapour interface), which we assume to be just below the boiling temperature, and Fourier law for the heat transfer in the film yield the order of magnitude for the vertical velocity of vapour w near the surface of the droplet:

w = λδT hLρ v . (1)
Here δT = T 1 -T 0 is the temperature difference between the droplet and the hot plate, ρ v is the mass density of the vapour, λ the heat conductivity in the vapour and L the latent heat of evaporation per unit mass. Lubrication theory, valid if l ≫ h, yields the magnitude of the horizontal velocity u in term of the vertical velocity : u = lw h , where l is the horizontal extent of the vapour film. From Poiseuille relation such a flow is driven by a horizontal pressure gradient scaling as δP/l ∼ ηu/h 2 . Therefore the typical gradient is of order δP/l ∼ wlη/h 3 . Using now equation (1) we find (replacing sign ∼ by = for readability) :

δP = ηl 2 λδT h 4 Lρ v , (2)
The liquid drop is at mechanical equilibrium when M g = δP l 2 namely when its weight is equal to the vertical component of the force due to viscous pressure in the gap between the droplet and the hot plate. Note that the contribution of the viscous stress to this vertical force is of the same order of magnitude as the one of hydrostatic pressure. Therefore we find using the mechanical equilibrium relation given above that :

l 4 h 4 = Lρ v M g ηδT λ . ( 3
)
As stressed above the lubrication approximation becomes invalid when h ∼ l. From the M dependence of the right-hand side of equation ( 3) we find that the lubrication approximation becomes invalid if the radius of the droplet is less than a critical drop radius R l defined by the condition l = h. The breakdown of the lubrication approximation will then occur for :

R l = ηδT λ gLρ v ρ l 1/3 . (4)
Therefore one expects that, as R gets smaller than R l , the droplet takes off from the hot plate. The range of values of h ∼ R ∼ R l is beyond reach of simple order of magnitude estimates and its analysis requires a solution of the full set of equations for the multi-phase Stokes flow and the heat transfer model. Let us consider the asymptotic regime where R is much less than R l , where one expects the droplet to lift well above the hot plate (h ≫ R). In this regime, the equilibrium elevation of the droplet is reached when the gravity force is balanced by the dipolar component of the velocity field generated by the evaporation from the droplet surface and by its image with respect to the horizontal plane. This dipolar approximation is correct for h ≫ R. Let T 1 the temperature at infinity and T 0 the temperature of the droplet. The temperature field satisfying Laplace's equation around an isolated droplet reads : of order T im ∼ (T 1 -T 0 )R/2h near the physical droplet, with a nearly uniform vertical gradient. As justified below, in the limit of large h the net vertical heat flux is of order of the derivative of T im with respect to h, namely (T 1 -T 0 )R/h 2 . Using Stefan's law we find the order of magnitude of the net vertical velocity in the vapour near the droplet reads :

T (r) = T 1 + (T 0 -T 1 )R/
λ Lρ v (T 1 -T 0 )R h 2 . ( 5
)
This relation is similar to equation ( 1), but with δT /h replaced by (T 1 -T 0 )R/h 2 . Multiplying this velocity by ηR, as in Stokes drag law, we obtain the vertical force on the droplet. Balancing this force with gravity we find that, in the limit h ≫ R (namely for R ≪ R l ), the dimensionless elevation of the droplet at equilibrium h * = h/R should scale as :

h * = (R l /R) 3/2 . ( 6
)
Before going further we have to analyse the constraint on the pressure in the gap: it must be negligible compared to the Laplace's pressure in the droplet to ensure that it remains close to its spherical shape everywhere. It can be shown [10] that l ∼ √ Rh is a fair estimate of the horizontal extent of the layer for which the viscous pressure contributes significantly to the vertical force. Using eq.( 2) and this relation between l, R and h, the pressure induced by the Poiseuille flow in the gap is found to be of order δp = gρ l RR 3 l h 3 . It crosses the order of magnitude of Laplace's pressure inside the droplet when δP = σ R . This happens when the radius becomes of order

R i = σ gρ l 2/7 R 3/7
l . This relation is obtained using the above geometrical relation for l and equation (3) reformulated as h/l = (R l /R) 3/4 . This assumes that the capillary length σ gρ l 1/2 is much larger than R l , since otherwise the cross over would happen in a range of values of R smaller than R l : this is impossible because R i was derived under the assumption that the lubrication approximation holds in the gap between the sphere and the hot plate. The condition R i ≫ R l is fulfilled in the experiments described below, but there could be other situations at stakes for which this condition may not be fulfilled. This gives us also the opportunity to address another feature of our experiments, namely the difference of behaviour between droplets of water and of ethanol. Since the quantity R i /R l is independent on the radius of the droplet, but depends on the properties of the liquid and its vapour, the behaviour of Leidenfrost droplets is not unique up to rescalings, because it depends on a dimensionless ratio with different value depending on the liquid used. Let us also notice that for radii larger than R i , the shape of the droplet becomes rather complex, because it depends on the (large) deformation of the droplet surface due to the Poiseuille pressure inside the vapour gap [11,12]. This happens in a range of parameters outside of this study. Therefore we focus on the range of parameters such that R l ≪ R ≪ R i . In this range, using equation (3) and the geometrical relation for l, the dimensionless height satisfies:

h * = (R l /R) 3/2 . ( 7
)
For usual liquids the value of R i is of a few hundred of microns. It is therefore worth mentioning that this prediction is qualitatively different and even opposed to the one commonly used for which the vapour thickness is expected to decrease when the radius of the droplet decreases.

The previous relations rely on the assumption that the the temperature field between the drop and the substrate (the plane located at z = 0) is a solution of Laplace's equation. This approximation is valid in the limit in which the convective heat flux is negligible compared to the diffusive flux (small Peclet number limit) and it holds in our experiments. From the image solution of potential problems, the temperature field between the drop and the hot plate is the same than the one between two spheres at temperatures T 0 and 2T 1 -T 0 and separated by a distance 2h. Laplace's equation for two spheres was solved one hundred years ago by Jeffery [13] thanks to bi-spherical coordinates. We have numerically recovered the expected dependence of the temperature gradient in the limits of h/R large and h/R small. As given above by scaling arguments, they respectively read ∇T ∼ ∆T R/h 2 and ∇T ∼ ∆T /h.

Considering a water drop on a hot plate at 400 0 C the numerical estimate of the critical radius is R l ≃ 19µm. To our knowledge no quantitative experiment has been done at such small sizes. We therefore decided to perform experiments on ultra-distilled water and ethanol drops to verify the scaling laws derived just above. The experimental set-up is depicted in Fig. 1. A cloud of submillimetric droplets is rained on a silicon substrate kept at a controlled temperature T 1 . The two upper blocks are necessary to reduce the number of drops falling in the field of the camera. They are also useful for the thermal stability of the system. Experiments are performed under a laminar air flow protection in order to prevent dust contamination. A high-speed camera, at a frame rate of 2000 frame per second, was used to record several evaporation take-off processes. An image analysis allowed us to measure the radius and elevation of droplets as a function of time. The typical error made on these measurements is 2µm. On the top of Fig. 2, we show six pictures taken during the take-off of a water droplet. On the bottom of Fig. 2 we plot the height of the water droplet deposited on a silicon substrate at T = 375 • C and room pressure. It first bounces and lands on the silicon substrate. It then enters in the lubrication regimes in which both radius and vapour thickness decrease with time. At certain time the droplet takes-off from the substrate (see the movie M1 in the supplementary materials).

We recorded several evaporation processes and we represent in Fig. 3 the dimensionless elevation of the droplet h * as a function of its radius R. The measurements are in qualitative agreement with the pictures of Fig. 2 (top) : almost all droplets take-off from the substrate when their radii becomes of order R l . The theoretical prediction given above is tested against experiments. The full line is a best fit with a single parameter a multiplying the predicted critical radius R l . We can see that the theoretical prediction fits rather well the experimental data with a value a = 2.2 close to unity. One can see that there is an important dispersion in sizes at which the droplets takeoff. It seems that, for the system of water droplets under study, the lift-off is analogous to a sub-critical transition (in other terms for the same radius and temperature difference, there are two equilibrium elevation in a certain range of values of R/R l ). The same behaviour has been observed for water droplets on substrates at different temperatures. Decreasing the temperature tends to decrease the dispersion of radius's at take-off. Nevertheless the discontinuous character of the escape from the lubrication regime remains.

We investigated also the behaviour of ethanol Leidenfrost droplet. The same experimental procedure is used and the silicon substrate is kept at similar temperatures than for the study of water Leidenfrost droplets. Contrarily to water droplets almost all ethanol droplets do not rebound and land on the substrate (see the movie M2 in the supplementary materials). This is due to the lower latent heat of ethanol and therefore to its larger evaporation rate. As a consequence the ethanol droplets enter directly in the regime described above for which h >> l and we do not observe a transition from the lubrication regime. The experimental results are presented in Fig. 4 in which we plot, as for water droplets, the dimensionless elevation of the droplet h * as a function of its radius R. The data are presented for 10 different ethanol droplets on a substrate kept at T = 400 • C. As stressed just above the escape from the lubrication regime is no longer present and the height h of the droplets increases in a continuous manner. In the insert of Fig. 4 we represent the same experimental data but plotted in a log-log representation. The black line corresponds to the scaling law predicted in equation (6). A rather good agreement is found between the theoretical prediction and the experiments.

In this letter we have put in evidence the unexpected behaviour of a Leidenfrost droplet in the latter stage of its evaporation. Below a critical size R l the droplets escape from the lubrication regime and take-off from the substrate. In an intermediary regime (R l ≪ R ≪ R i ) the thickness of the vapour film increases as the radius of the droplet decreases. This latter prediction, based upon a detailed analysis of the various physical phenomena involved, is qualitatively different, and even opposed to the one commonly used: as a result of the standard approach to the Leidenfrost phenomenon the vapour thickness is predicted to decrease as the radius of the droplet decreases. The different scaling laws presented in this paper are in semi-quantitative agreement with the experimental data obtained on water and ethanol Leidenfrost droplets. In particular, as observed and predicted, the thickness of the vapour film increases as the droplet become smaller. Besides its fundamental interest this study should find implications in many domains such as for example in diesel combustion engines or in heat transfer us-ing cooling spray [14]. We hope that this study will motivate numerical approaches complementary to the present study.

where p0 is a constant so that the pressure decreases to zero at infinity. Using the expression given just above for h(r) it can be shown that the pressure decreases rapidly to 0 as r -6 on a scale length √ Rh0. Therefore √ Rh0 is a fair estimate of the horizontal extent of the layer where the viscous pressure due to the evaporation contributes significantly to the lifting force.

FIG. 1 :

 1 FIG.1: Experimental set-up. A cloud of sub-millimetric droplets is rained toward a silicon wafer inserted in a copper bloc kept at a controlled temperature T1. The droplet of radius R is at temperature T0 and at a distance h of the substrate.

FIG. 2 :

 2 FIG. 2: Top : Pictures of a Leidenfrost droplet and its reflected image on the silicon wafer. Bottom : Height of the droplet as a function of time. After being dropped, it first bounces and lands on the substrate. The droplet then stands on the substrate in the lubrication regime. It finally takes-off from the substrate.

FIG. 3 : 2 FIG. 4 :

 324 FIG.3: Dimensionless height of the droplets as a function of their radii. The red line is a best fit to the equation (6). Note that a minority of droplets fall back on the substrate. This is certainly due to dust contamination. The same data are represented in the insert but in log-log. The dotted line indicates the value of R l for water Leidenfrost droplets on a substrate kept at T = 400 • C.