Robust high-cycle fatigue stress threshold optimization under uncertain loadings

This paper proposes a strategy to achieve robust optimization of structures against high-cycle fatigue when a potentially large number of uncertain load cases are considered. The strategy is heavily based on a convexity property of some of the most commonly used high-cycle design criteria. The convexity property is rigorously proven for the Crossland fatigue criterion. The proof uses a perturbation technique and involves the principal stress components and analytical expressions for the applicable fatigue criteria. The multiplicity of load cases is treated using load ratios which are bounded but are otherwise free to vary within certain limits. The strategy is applied to a notched plate subject to traditional normal and shear loadings that possess uncertain or unspecified components.


INTRODUCTION
Traditionally, structural design does not take into account fatigue life issues.Engineers conceive projects based on parameters such as stress limits, maximum displacements, fundamental frequency and buckling.Nevertheless, structures break down even when admissible loads, determined by the design criteria above, are applied.Commonly, it occurs after a long time the structure becomes operational and the cause of failure is the end of the component fatigue life.
In the aerospace industry several aircraft structural components deserve full attention from the fatigue point of view.In the fuselage upper panels, windshield, door frames, pressure domes, wing and empennage connections and engine pylon are sensitive components.In the wing lower panels, spars, landing gear fixtures, flap and aileron connections, and engine pylon are all prone to fatigue.Moreover, vertical and horizontal tails should be carefully designed for fatigue.
In order to ensure infinite fatigue life for metallic structures loaded by multiaxial timevarying loadings, designers must include a multiaxial high-cycle fatigue (MHCF) criterion in the structural analysis [11].However, it is very difficult to evaluate the optimal structure configuration in the case of uncertain loadings, because it is required to integrate an optimization method with a robustness strategy.
Mrzyg lód and Zieliński [9,10], without using a robust technique, implemented a structural optimization analysis adapting the time load history in a car front suspension arm regarding eight geometric decision variables, using the Dang Van fatigue criterion as a constraint and minimizing the objective function computed by the arm mass.The optimization method applied was a probabilistic search based on evolutionary algorithms, with an optimum result of 5.6% decrease of arm mass at the 29th generation.Certainly, a large number of finite element simulations were performed to achieve this result, spending a large calculation time.
Steenackers et al. [13] integrated an optimization method with a robustness strategy applied to the slat track (airplane component) finite element model taking into account the uncertainty of design parameters.They concluded that in order to obtain accurate results, an extensive computation time is necessary to run a large number of Monte Carlo simulations in combination with FE models.
Jung and Lee [5] proposed a method for robust optimization problems which conjugates the probability in design process and sensitivity analysis.Moreover they used the advanced first-order second moment method in order to reduce the calculation time.This method has the advantages of gradient-based optimization methods, since the second-order sensitivity information is not required, but the first-order sensitivity information is still necessary.d'Ippolito et al. [3] introduce a methodology based on reliability analysis and design optimization in order to ensure robustness for their structural design considering fatigue life of its components.A probabilistic optimization method to fatigue life based on strain-life approach was implemented.The authors used a hybrid meta-model/FE strategy to save computation time.Besides a quadratic response surface model based on design of experiments (DOE) results was necessary in order to reduce the number of FE analysis.
This work presents a new robust approach in metallic structural design optimization problem considering (multiaxial) high-cycle fatigue criteria.For the purpose of saving computational time, removing sensitivity analysis of design variables and performing a minimum number of finite element analyses as a possible, the optimization strategy is based on extreme modeling which adds robustness to the optimal design.It is necessary to combine this technique with convexity of fatigue criteria.A numerical example was solved applying the proposed robust optimization technique and the results were discussed.

LOADING REPRESENTATION AND PROBLEM FORMULATION
Consider an arbitrary structure with boundary conditions and L applied loadings shown in Fig. 1.A 2D sketch is drawn to facilitate visualization but 3D structures can also be resolved under the approach to be developed.The regions of application of the loadings are assumed to be fixed.This assumption rules out the class of moving loads [2] but is generally accepted in the context of fatigue investigations.Moreover, moving load problems induce relevant dynamic effects which are disregarded in high-cycle fatigue studies.
The stress state of an arbitrary point x is given by the Cauchy stress tensor σ resulting Latin American Journal of Solids and Structures 9(2012) 615 -631 Applied loadings from all the loadings applied.Admitting that the loadings vary with time, the stress tensor σ also varies with time.Additionally, since linear elasticity is assumed, the resulting stress state is the summation of the stresses produced by individual loadings: where σ l relates to the loading p l .The principle of superposition embedded in Eq. ( 1) is valid for linear elasticity.In nonlinear problems Eq.( 1) is no longer applicable.However, since the high-cycle fatigue stress thresholds are well below the yield stress, linearity is a very reasonable assumption in the present analysis.
The load ratios r 1 , ..., r L represent the relative contribution of the loadings and are bounded by 0 ≤ r l ≤ 1 for all l = 1, ..., L. Usually the loadings are not applied simultaneously and the load ratios may obey a given relationship to reflect it.A useful relation is the linear convex combination r 1 + ... + r L = 1 which assumes that the relative contribution of the loadings add up to 100%.A natural question that arises refers to the most dangerous loading combination.As an example, consider only two loadings: p 1 and p 2 .What combination leads to the worst situation in terms of fatigue: (i) r 1 = 1, r 2 = 0, (ii) r 1 = 0, r 2 = 1, (iii) r 1 = 0.5, r 2 = 0.5 or (iv) another combination where r 1 + r 2 = 1?This question will be answered in this paper.
Any loading p l is written in terms of two components as in Eq. ( 2) where the superscripts m and a were adopted to comply with usual terminology [6] in fatigue investigations: m relates to the mean (average) of the loading component whereas a reflects amplitude and time variability.Parameter r a l (t) varies with time and is usually harmonic (e.g.sin(ωt + ϕ)) in traditional fatigue investigations.p m l and p a l are nominal or reference loadings that may be known from the loading envelope of the structure.
When r a l (t) are specified functions of time and the load ratios r l introduced in Eq. ( 1) are fixed (or known) a traditional optimization procedure can be proposed considering the loading trajectory in time along with a suitable multiaxial fatigue criteria [11,12].On the other hand, whenever uncertainties are present in either r l or r a l (t), more elaborate strategies are required including, for instance, a statistical treatment of the applicable loadings [8].
The optimization strategy proposed in this work to handle the lack of specification in both r l and r a l (t) is based on extreme modeling.The extreme modeling adds robustness to the optimization effectively eliminating sensitivity of optimal designs against variations in the uncertain loadings.However, in order to be economical, the extreme modeling technique should work in conjunction with convexity properties.In essence the extreme modeling technique obtains optimal designs against not only one particular load case, but against an entire space of admissible loadings.
Uncertainty is introduced in the problem admitting that r l and r a l are not known beforehand but vary within bounds for all instants of time t.These box constraints, mathematically defined in Eq. ( 3), guarantee a bounded problem.In practice this means that a reasonable estimative of the range of variability of r l and r a l must be available.
Notice that Eq. ( 3) does not assume beforehand how the load ratios may vary with time.This is an important aspect of the loading representation since nonharmonic and asynchronous loadings can be considered.Therefore, time variability of r a l (t) becomes unimportant since what essentially matters are the side constraints (r a l ) MIN and (r a l ) MAX .Statistically the interpretation for time independence is that r a l is assumed to have uniform probability of attaining any value contained within the box defined in Eq. (3c).
The traditional optimization problem may be posed as volume minimization under the constraint that a given multiaxial fatigue criterion is satisfied, e.g.Dang Van [14], Crossland [1], McDiarmid [7], Findley [4], Mamiya [6], among others: where θ is the vector of design variables such as thickness distribution over a plate or cross section dimensions of stiffeners in a reinforced panel, V (θ) is the total structural volume, σ c is the critical fatigue stress which depends on the stress tensor σ at x, b is the fatigue stress threshold, and r(t Notice that r is a function of time and so is σ c . Equation (4) shows that the most demanding procedure in the solution of the optimization problem is the computation of the constraints which involves double maximization of σ c over x and r(t).Additionally, the constraint in Eq. ( 4) is nondifferentiable with respect to θ since the worst r certainly varies with θ.The next section shows how to significantly simplify the double maximization problem required to evaluate the constraints present in Eq. ( 4).

CONVEXITY OF FATIGUE CRITERIA
The critical fatigue stress σ c is calculated based on the stress tensor σ at a particular point x.The stress tensor, in turn, depends on the load ratios r.The approach to solve the double Latin American Journal of Solids and Structures 9(2012) 615 -631 maximization problem is twofold.Firstly, consider a fixed point x 1 and vary r.Once the most dangerous load ratios r 1 are obtained for point x 1 another point x 2 is to be assessed and its most dangerous load ratios r 2 obtained.Observe that, since time independence is assumed, usually r 1 ≠ r 2 .
Computation of the worst r is an important part of the constraint evaluation which can be sped up if σ c is more carefully examined.Since most multiaxial stress fatigue criteria depend on the principal stresses of σ it is worth investigating how the principal stresses vary as σ is perturbed.Assume that vector r produces a given stress tensor at point x which can be computed as in where σ m l and σ a l relate, respectively, to the loadings p m l and p a l defined in Eq. ( 2).The newly defined parameter r l (t) is simply the product r l × r a l (t) and, according to Eqs. (3a) and (3c), obey the relationship Equation ( 5) is linear in r 1 , r 1 , ..., r L , r L .Therefore, a first order perturbation in these parameters lead to a first order perturbation in σ of the type σ + δσ, not containing higher order terms δ 2 σ, δ 3 σ, etc.
The principal stresses are computed from the eigenvalue problem stated in Eq. ( 7) where I is the identity tensor, λ is the principal stress and n is the principal direction.The perturbed eigenvalue problem is Equation ( 8) can be split into higher order problems such that the first order problem is and the second order problem is Premultiplication of Eq. ( 9) by n T , use of Eq. ( 7), symmetry of σ and the fact that Equation (11) shows how the principal stresses vary when r is perturbed.However, not much useful information can be retrieved from it.Premultiplication of Eq. ( 10) by n T , use of Eqs. ( 7) and ( 9), and the fact that n T n = 1 yields Equation ( 12) gives important information regarding the sign of the second derivative of λ (the principal stresses).Notice that σ can be written as σ = NΛN T where N is the matrix of eigenvectors (stored columnwise) and Λ is the diagonal matrix whose trace is composed of the principal stresses σ I , σ II , σ III with σ I ≥ σ II ≥ σ III .It is clear that N T N = NN T = I holds.Hence, Eq. ( 12) can be rewritten as There are two important situations: (i) λ = σ I and (ii) λ = σ III .When λ = σ I the diagonal matrix in Eq. ( 13) is negative semi-definite and when λ = σ III it is positive semi-definite.In one case δ 2 λ = δ 2 σ I ≥ 0 and in the other The conclusions just drawn can be immediately employed to prove convexity with respect to r 1 , r 1 , ..., r L , r L of one of the most commonly used high-cycle fatigue criterion in finite element analyses: Crossland [1].
The Crossland (σ CR ) stresses can be expressed as where κ CR is a parameter that can be expressed from data of two fatigue tests (reversed bending and reversed torsion), σ H = (σ I + σ II + σ III )/3 is the hydrostatic stress (three times the trace of the stress tensor) and J 2 is the second invariant of the stress deviator tensor given by The square root of J 2 appearing in Eq. ( 14) is an amplitude and not an instantaneous value.The calculation of this amplitude necessitates the definition of a complete load cycle.From a complete load cycle one can calculate the amplitude of the square root of J 2 and the maximum hydrostatic stress σ H occurring within this cycle to apply the Crossland criterion.Notice, however, that time dependancy is not an issue anymore since extreme modeling is being applied (see Eqs. (3a) and ( 6)).Since trace (3σ H ) is a linear function of σ I , σ II , σ III its second variation (or derivative) is zero, i.e., δ 2 σ H = 0. On the other hand, J 2 is a quadratic function of σ I , σ II , σ III meaning that its second variation must be nonnegative, what also implies that the second variation of √ J 2 is also nonnegative, i.e., δ 2 √ J 2 ≥ 0. It follows that Latin American Journal of Solids and Structures 9(2012) 615 -631 As an example, consider two reference stress states given by and assume two load ratios r 1 = r and r 2 = 1 − r, 0 ≤ r ≤ 1, are associated with σ m 1 and σ m 2 , respectively, such that the resulting stress state is σ = r 1 σ m 1 + r 2 σ m 2 .Taking κ CR = 0.5 it is possible to draw the stress curves shown in Fig. 2. As expected, it is clear that σ CR is a convex function of r.Moreover, notice that trace is a linear function or r, what is also expected since δ 2 σ H = 0. Nothing can be stated about the convexity of the second principal stress σ II .Actually, it is seen that the sign of δ 2 σ II is either positive or negative.

ROBUST OPTIMIZATION STRATEGY
Equation ( 4) states the robust optimization problem that must be solved in order to find the optimal design θ.In the previous section it was proven that the critical fatigue stress σ c is convex with respect to variations in the applied loadings r, i.e., δ 2 σ c ≥ 0, at least for the Crossland criterion.
Latin American Journal of Solids and Structures 9(2012) 615 -631 Convexity of σ c with respect to variations in the load ratios can be schematically seen in Fig. 3a.The critical stress surface is constructed by computing σ c for different load ratios represented generically by r i and r j .It must be clear that only two axes of load ratios are drawn to permit visualization in 3D.However, many more may be present as explained in section 2. The load ratios are represented by dots in the r i r j plane.Notice that the loading envelope described in Eqs.(3a) and ( 6) is highlighted in Fig. 3a and termed convex hull.This terminology is justified by the fact that any other set of load ratios, represented by dots within the convex hull, can be written as a linear convex combination of the vertices V 1 − V 4 of the convex hull.Figure 3b shows the geometrical representation of Eq. (3b) where the plane r 1 + r 2 + ... + r L = 1 can be seen.
With regard to the optimization strategy, constraint evaluation in Eq. ( 4) requires solution of the inner loop that maximizes σ c with respect to r.This inner maximization problem can be solved by checking the projections P 1 −P 4 of the convex hull vertices V 1 −V 4 onto the critical stress surface.The worst possible load ratios is the one whose point P k is associated with the greatest σ c .Notice that, in a situation where there are several loadings applied, even using the convexity of σ c entails computation of a very large number of possibilities.For a structure with L loadings the number of possibilities is L × 2 L , where L comes from the vertices defined by Eq. (3b) and 2 L comes from the box constraints defined in Eq. ( 6).Since linearity of the mechanical response is assumed, 2L load cases can be computed beforehand and the result stored for evaluation of the worst σ c , being L loadings p m l and L loadings p a l .The next section presents an example to clarify this point.The outer loop necessary to compute the constraints in Eq. ( 4) requires a search over the complete structure domain since x is the position of an arbitrary point.Mathematically, searching over the structure domain is the most precise procedure.In practice, however, finite element models are usually used to model complex structures.Therefore, a simpler alternative for the outer maximization loop consists in evaluating σ c at the nodes.

NUMERICAL EXAMPLE
The rectangular notched plate shown in Fig. 4 is used as an example.It has length 1.0 m, width 0.5 m and circular notch with radius 0.1 m located at the plate center.Notice that the plate has a nonuniform thickness distribution and the loadings may have mean and amplitude components as expressed by Eq. ( 2).A finite element code was developed to obtain the critical fatigue stress σ c based on the Crossland criterion.This code employs the 4-noded quadrilateral element with the usual bilinear shape functions The element stiffness matrix K e can be computed by where and the displacements u and v can be interpolated using the respective nodal values Notice that the integral in Eq. ( 18) involves the thickness h(x, y) which varies over the element domain.The thickness distribution of the plate is assumed to be adequately represented by nodal thicknesses such that, similarly to Eq. ( 20), it can be interpolated as and integration in Eq. ( 18) is numerically carried out using 2 × 2 Gaussian quadrature.The mesh shown in Fig. 5 with bilinear elements was used to model the plate.Refinement around the notch boundary was ensured because this is the region where stress concentrations shall occur.The optimization strategy selected to solve Eq. ( 4) is topometric, i.e., plate volume minimization is achieved through thickness variation.The procedure implemented is simple: Figure 6 presents the flowchart of the optimization procedure described here and in the previous section.The stresses evaluated at the nodal positions are computed on average, i.e., if two or more elements share node I, then the stress is given by the average of the stresses computed using interpolation functions and nodal displacements of these two or more elements.
Closer examination of the 24 vertices reveals that only 4 are necessarily the most dangerous.
Latin American Journal of Solids and Structures 9(2012) 615 -631 Tables 2-4 give clear account about how this conclusion is drawn.Consider for instance Tab. 2 where r 1 = 1 and r 2 = r 3 = 0. r 1 relates to the normal load N xx , r 2 relates to the normal load N yy and r 3 relates to the shear load N xy .The first three columns of Tab. 2 correspond to all possible permutations of r a 1 , r a 2 and r a 3 which can assume only two values (-1 or +1) according to Tab. 1. Columns 4-6 of Tab. 2 are the values of r 1 , r 2 and r 3 computed by the product r l = r l r a l (see Eq. ( 5)).Finally, columns 7-9 of Tab. 2 are computed using Tab. 1 and Eq.(5).For example, the fifth row is obtained computing Examining Tabs.2-4 it can be concluded that most of the loading vertices correspond to null loadings such that only four of them must be evaluated: (i) . This simplification is not always as straightforward since the existence of nonzero loading vertices and complex geometries of structural components commonly lead to situations where obvious undangerous loadings cannot be discarded beforehand.

Table 2 Determination of loading vertices for r
Latin American Journal of Solids and Structures 9(2012) 615 -631 Table 4 Determination of loading vertices for In order to illustrate the weaknesses of the conventional optimization against single loadings five optimization cases are considered as presented in Tab. 5. Cases I-IV correspond to single loads whereas case V corresponds to the robust optimization case.The '±' in case V indicates that either positive or negative shear may be applied as just discussed.It must be clear that, as implied by Eq. ( 5), loadings in case V are not applied simultaneously; it is actually a varying load case where the four load cases I-IV of Tab. 5 are applied as a convex combination σ = ∑ 4 l=1 r l σ l , with ∑ 4 l=1 r l = 1 and 0 ≤ r l ≤ 1.As explained in section 4 the loading ratios r 1 , r 1 , ..., r L , r L are free to vary but must satisfy Eqs. ( 3) and (6).For optimization purposes, it was shown in section 4 that it is sufficient to check the vertices of the loading space corresponding to N xx = 800 kN/m, N yy = 200 kN/m, N xy = +400 kN/m and N xy = −400 kN/m, one at a time.However, the optimal design thereby obtained is guaranteed to sustain any linear convex combination of these four loadings.6 presents the volume of the optimal designs obtained and the critical fatigue stresses calculated at the loading vertices.The highest volume is associated with optimal design V since this is the robust design, i.e., the only one that can equally well sustain all load cases (σ c ≤ 160 MPa).The vulnerability of the traditional optimal designs, cases I-IV, is evident.For instance, optimal design I is able to sustain load case I but it fails for all the other load cases II-IV (stresses much higher than the threshold b = 160 MPa).The last column of Tab. 6 is simply the maximum value of the critical fatigue stresses.A visualization of the optimization procedure is shown in Fig. 7. Figures without scales are drawn since the intention is to give a qualitative view of the stress and thickness distributions.Rows 1 to 5 correspond to the five effective loadings described in Tab. 5.In the first column are the critical fatigue stress results for the uniform plate.In the second column are the thickness distributions of the optimal designs.In the third and last column are the the critical fatigue stress results for the optimal plates.All optimization processes began with a uniform plate 1.0 mm thick, which corresponds to an initial volume of 4.6858 × 10 −4 m 3 .As shown in Fig. 7 (first column) stress concentrations occur around the notch.Therefore, although the initial structural volume is small (4.6858 × 10 −4 m 3 ) compared to those obtained in Tab.6, the initial uniform plate suffers from fatigue under all loadings.
It is clear that stress concentration regions around the notch are observed in all cases of uniform thickness distribution (red regions), whereas stresses are much more spread in the optimal designs (last column in Fig. 7).Cases III and IV are equivalent since they correspond simply to a change in sign of the shear loads.The thickness distribution of the robust optimal design (last row) is such that it reinforces the region around the notch where stress concentrations seen in the first column arise.Figure 7 shows optimal designs with unusual shapes (second column).It is obvious that the stress is concentrated around the notch (with a maximum value of 3 times the applied stress in single tension) and increasing the thickness can decrease this maximum value but, the question whether these shapes are acceptable is dependant on the application.In practice, some kind of smoothing may be necessary prior to actual production of optimal designs, similarly to what is observed in topology optimization.
The critical stress surface can be visualized for the plate example simulated.for the robust optimal design, 0 < σ * ≤ 1.The highest values of σ * happen in the region around the edges of the reference plane.
For the sake of comparison against Fig. 8, Fig. 9 depicts the critical stress surface obtained for the optimal design under load case I solely.Though still concave, the critical stress surface is now severely distorted.The normalized maximum critical fatigue stress σ * assumes now values of up to 15.7 which is way beyond the fatigue threshold b.

CONCLUSIONS
Robust optimal designs may be heavier than designs optimized against one individual load case but they are definitely safer.However, if multiple load cases are to be handled by a traditional optimization strategy, the resulting optimal design would be even heavier than the one obtained through the robust optimization technique proposed.It must be clear that the robust optimal designs obtained through the optimization procedure proposed can safely sustain a potentially large number of individual loadings, identified as vertices of the convex hull shown in Fig. 3, as well as any convex combination of these vertices.This property ensures that eventual nonharmonic or asynchronous, not predicted in the loading envelope, are not as harmful.The optimal results are robust against fatigue.However, from a practical view point, that is, if one must design a real structure, even the robust optimal design would hardly be manufactured.As a suggestion to obtain a viable optimal design which can be used in practice, geometric constraints could be included on the volume minimization.The simplest possibility is to consider that the plate is composed of a base plate whose thickness is uniform and fixed and external layers of variable thickness exist.A more elaborate proposal is to add constraints on the maximum thickness gradient such that the maximum difference in thickness between adjacent elements must not exceed a certain amount.
The topometric numerical optimization adopted may not result in strictly optimal designs in the global sense.However, it is a fast procedure that does not require evaluation of gradients.Notice that, although gradient of volume is easily computed, the max x max r σ c in Eq. ( 4) is a nondifferentiable constraint.Hence, the use of gradient-based optimization procedures would require considerable effort and ingenuity to overcome the lack of necessary smoothness in the constraint.

Figure 2
Figure 2 Stress curves for varying load ratios

Figure 3
Figure 3 Convexity of critical fatigue stress

Figure 4
Figure 4 Notched plate with loadings with nonuniform thickness distribution

Figure 5
Figure 5 Meshed notched plate i.For a given thickness distribution θ i find the worst σ c at all nodes ii.For each node I, update thickness h I,i+1 = h I,i * σ cI /b iii.Stop when all h I do not change considerably Latin American Journal of Solids and Structures 9(2012) 615 -631

select initial thickness distribution i = i + 1 solve L × 2 Figure 6
Figure 6 Optimization flowchart

3 relate to r 4 .Figure 7
Figure 7 Optimal design results

3 *Figure 9
Figure 9 Maximum critical stress surface for optimal design of case I

Table 1
Definition of loadings and load ratios

Table 5
Effective loadings for optimization

Table 6
Critical fatigue stress for optimal designs * Figure 8 Maximum critical stress surface for robust optimal design *