A Shakedown Strength Based Parametric Optimization Technique and Its Application on an Airtight Module

In aerospace engineering, design and optimization of mechanical structures are usually performed with respect to elastic limit. Besides causing insufficient use of the material, such design concept fails to meet the ever growing needs of the light weight design. To remedy this problem, in the present study, a shakedown theory based numerical approach for performing parametric optimization is presented. Within this approach, strength of the structure is measured by its shakedown limit calculated from the direct method. The numerical method developed for the structural optimization consists of nested loops: the inner loop adopts the interior point method to solve shakedown problems pertained to fixed design parameters, while the outer loop employs the genetic algorithm to find optimal design parameters leading to the greatest shakedown limit. The method established is first verified by the classic plate-with-a-circular-hole example, and after that it is applied to an airtight module for determining few key design parameters. By carefully analyzing results generated during the optimization process, it is convinced that the approach can become a viable means for designing similar aerospace structures.


Introduction
In aerospace engineering, the goals of reducing the launch costs and allowing multiple landers to be packaged on a single launch vehicle become urgent [1], therefore, the demand for lightweight structure is rising rapidly. In the design process of spacecraft, the traditional elastic limit rule and its affiliated relax factors are frequently used. Due to the intrinsic conflict between weight and loading capacity, this approach often leads to an unacceptable overall weight [2][3][4][5]. Although reusable spacecraft can significantly reduce the cost for space travel which makes it a clear trend for future development, the corresponding evaluation criteria of the structural reliability are still somehow missing [6]. Therefore, it is particularly important to establish a new design and evaluation approach customized for lightweight structures appeared in reusable spacecrafts. To this end, the shakedown theory could be a viable means since it is well recognized and has been implemented in numerous codes related to strength and safety assessment. According to Melan's renowned shakedown theorem, the maximum load bearing capacity of a structure can be obtained without following the entire load history. For this reason, the method is also named as direct method (DM) [7][8][9]. In DM, the plastic limit and shakedown limit, where the former one can be regarded as a particular case of the latter one, correspond to thresholds of plastic failure of a structure subjected to monotonic and cyclic loading, respectively.
When designing structures used in spacecrafts, complicated part geometries always lead to an inhomogeneous stress distribution. As a consequence, if in the whole structure the stress level is required to be well below the

Open Access
Chinese Journal of Mechanical Engineering *Correspondence: gengchen@bjtu.edu.cn yield limit, geometrical parameters have to be selected very carefully. Despite its simplicity in practice, such design concept is neither necessary nor efficient. This is because, on the one hand, local stress exceeding the yield stress does not mean that the structure completely loses its loading capacity, while on the other hand, adopting such a design concept would entail an unacceptable structural weight. For these reasons, one viable way to attain a desired strength-to-weight ratio is to change design criterion from elastic limit to shakedown limit. Calculating shakedown limit by DM requires solving large scale nonlinear convex optimization problems which are computationally expensive to be solved [10][11][12]. This prevents the method from being prevalently adopted as a standard technique for assessing the structural strength in engineering practices. In early years, optimization problem entailed by DM have to be first converted to linear ones and then solved by linear programming techniques such as Simplex [13]. It is only until early 90s when the interior-point method (IPM) [14,15] becomes matured that these problems can be solved by nonlinear programming algorithms [16]. One major advantage of the IPM is that, compared to active set methods such as sequential quadratic programming, it avoids the problem of selecting active constraints. Meanwhile, generally IPM requires only few iterations to converge [17]. Due to these merits, IPM is employed in many commercial optimization software such as CPLEX, LOQO, Mosek, IPOPT, and Gurobi, and some of these software have been successfully adopted to solve DM problems [9,11,[18][19][20].
Limited by the computational effort to solve the large-scale optimization problem, shakedown analysis is rarely used in conjunction with structural optimization algorithms. To the best knowledge of the authors, most existing researches focus on analyzing parametric sensitivities. For instance, Yang [21] proposed a revised sequential quadratic programming algorithm with an active-constraint set strategy to solve shakedown problems, and by studying various straight pipelines and pipe elbows with this method, the failure mode and the effect of shapes and sizes of typical part-through slots on plastic and shakedown limits were clarified. Sun et al. [22] introduced a numerical method to evaluate the ultimate bearing capacity of a semisubmersible brace considering cyclic wave loads and clarified how load angle, shell thickness, stiffener thickness and stiffener spacing influence the shakedown limit. Chen et al. [23] investigated the impact of few parameters, such as the bend radius, the pipe's mean radius and wall thickness on plastic and shakedown limits of a 90° pipe under two load cases. Zheng and Xuan [24] utilized a direct cycle method to study the shakedown limit of perforated thick-walled cylinders under cyclic thermo-mechanical loads and revealed that the shakedown limits decrease when the hole radius becomes larger while the wall becomes thinner, especially when the structure is in an axial compressive state. Cho and Chen [25] assessed the plastic, shakedown, and ratchet limits of a 90° back-to-back pipe bend structure and performed sensitivity analysis of few design parameters using the so-called Linear Matching Method (LMM). Peng et al. [26] carried out the limit and shakedown analysis by means of the Stress Compensation Method (SCM) on 45° piping elbows with various geometries subjected to a steady internal pressure and a cyclic in-plane closing, opening and reverse bending moments. Based on many calculations it was pointed out that the elbows' design parameters need to be carefully designed such that a sufficient margin between two values can be attained which prevents unexpected plastic collapse during cyclic loads from happening. In the spirit of the direct method, Kammoun and Smaoui [27,28] developed a numerical method for solving optimization problems where the objective function is to minimize the overall weight, and the constraints require stress field to satisfy all static conditions within the shakedown theorem.
Although all these works attempted to employ shakedown limit as an evaluation criterion and based on that have successfully established qualitative relationships between design parameters and admissible shakedown load domains, an optimization tool which allows computer to autonomously determine optimal design parameters is still missing. For this reason, the goal of the present work is to create such a tool to perform parametric optimization with respect to the shakedown limit. The algorithm for the tool consists of nested loops: the inner loop adopts the IPM to solve shakedown problems pertained to fixed design parameters, while the outer loop seeks for optimal design parameters that lead to the greatest shakedown limit. Both in above mentioned works and in our own studies [20,29,30] it is noticed that the sensitivity of the shakedown limit to a certain design parameter is neither monotonic nor linear, especially when structures own complex geometries. For this reason, if the outer loop employs a gradient-based optimizer, it is highly likely that the solver would converge very slowly or only to a local minimum. To solve this problem, genetic algorithm (GA) is chosen as the optimizer of the outer loop. As a derivative-free algorithm, GA evaluates only the objective function and does not require any auxiliary information. In this study, as we are only interested in the application of GA, a well-established GA tool was used and more detailed information on theory of GA can be found in Refs. [31][32][33][34].
The remaining part of this paper is organized as follows: followed by the introduction, the shakedown theory and its numerical implementation are briefly elaborated. Next to that, a shakedown based parametric optimization framework is discussed together with strategies to improve both robustness and efficiency of the solving process. The validity of the proposed method is first examined by the classic plate-with-a-circular-hole benchmark example, then it is applied to a manned airtight module to study how optimal design parameters changes when the objective function switches between elastic, plastic and shakedown limit. Finally, based on results of these numerical studies, the advantages of employing DM over elastic criteria for optimizing design parameters and the great potential of the method for designing similar aerospace structures would be clarified.

Methodology
Before shakedown based parametric optimization method is elucidated, the lower bound problem for numerical simulation is briefly introduced here.

Lower Bound Limit and Shakedown Analysis
In the present study, the structure is idealized as an elastic-perfectly plastic material which obeys the von Mises yield condition and omits both geometrical nonlinearity and ductile damage. Considering a body B without body forces and subjected to multiple arbitrary varying loads, the total stress field σ(x, t) can be separated into two parts: a purely elastic stress σ E (x, t) and a residual stress where x indicates the position and t indicates the time.
After the decomposition, equilibrium conditions and boundary conditions have the following form: Here, is the material body, Ŵ 1 is the static boundary, n is the surface normal, f A is the surface loading. For the given load condition, there exists a load level under which although the plastic deformation would take place during first few load cycles, within the subsequent process the structure behaves purely elastically. This phenomenon is the so-called elastic shakedown. According to Melan's theorem, elastic shakedown happens if a timeindependent residual stress field ρ(x) can be found and its superposition with σ E multiplied by α ( α>1) satisfies the yield condition at any time t and any point x . Mathematical formulation for Melan' s theorem can be written as: where F stands for the von-Mises yield condition and σ Y is the yield strength, respectively.
When the structure B is simultaneously submitted to NL linearly independent varying loads P n (x, t) , then H(x, t) can be defined as a superposition of these loads: Here, P 0n (x) is a time independent loading basis, µ − n and µ + n are the lower and upper bounds of load multiplier µ n , respectively. A convex load domain with NV = 2 NL load vertices P k (k = 1, 2,...,NV) can be constructed. Here NL stands for the number of loads and NV the number of load vertices. König [35] has proven that for materials with a convex yield surface, if the shakedown condition is satisfied at all NV load vertices, then for any specific load history H within the domain the shakedown state will be attained. For this reason, the bearing capacity of a structure to multiple time varying loads can be obtained by solving following load-path independent optimization problem Here, α is the load factor and σ E (µ k ) the elastic stress at load vertex µ k . When two loads are acting on the structure (NL=2), NV=1 refers to the case when both loads vary monotonically and in this circumstance α corresponds to the plastic limit. Moreover, when NV=2 or 4, then two loads are time varying. In the former case, two loads are dependent and vary proportionally while in the latter case, they are independent and vary disproportionally. For both cases, α calculated is the shakedown limit. To solve the shakedown problem Eq. (5) more efficiently, a reformulation proposed by Akoa et al. [11] is adopted. Details of this reformulation can be found in our previous researches [29,30] or in the Appendix of this paper.

Shakedown Based Parametric Optimization
When the load bearing capacity of a structure to multiple time-varying loads is considered to be the objective for a design. Then in the simplest sense, the structural design can be expressed as the following parametric optimization problem where r stands for a vector of design parameters, f (r) is the objective function. In the present study, the design objective considered is either the load factor α calculated from Eq. (5) or the strength-to-weight ratio η = α SD m where m is the overall weight of the structure. For a parameterized model characterized by r , elastic stress for load vertex µ k and the constant residual stress can be expressed as σ E (µ k , r) and ρ(r) , respectively. Because shakedown limit α for design parameters r is calculated from Eq. (5), the major difficulty for solving Eq. (6) lies in evaluating the objective function: on the one hand, as has been figured out in the introduction, α is a nonlinear and non-convex function of r , therefore a global optimum can be only obtained using non-gradient algorithm.
On the other hand, because the computational effort for calculating an α from a given r is enormous, one should subject to r min ≤ r ≤ r max , avoid performing unnecessary shakedown calculations with respect to similar r ′ s . To this end, a pragmatic database-based approach was adopted. In this approach, a database is constructed when the optimization process starts, and once a trial r is generated, its Euclidean distance to all existing records will be evaluated. Only when this distance is greater than a threshold value, which means no similar record can be found, then shakedown analysis Eq. (5) will be performed and the newly derived results will be added to the database as a new record.

Implementation of Shakedown Based Parametric Optimization
The flowchart of the proposed nested genetic-gradient algorithm for solving the optimization problem Eq. (6) is shown in Figure 1, and the procedure can be outlined as following steps.
Step 1. Initialization: Create a random initial population consisted of 50 r samples lying within predefined ranges; calculate α corresponds to each r and write them to a database. This information is necessary for a warm start.
Step 2. Fitness calculation: In GA, the fitness function is a combination of objective function and constraints. Because all constraints considered in the present work are bounds of design variables and the GA optimization is performed in the MATLAB GA Toolbox, as a consequence, the objective function and bounds of variables are required to be provided separately. Based on this information, MATLAB will automatically generate a fitness function and employ it to evaluate the fitness of each individual. This way, within each generation, all samples are evaluated by the following sub-steps. will be written to the database and sent back to the main algorithm process simultaneously.
Step 3. Examine the stop conditions and create the next generation: The algorithm stops when the average relative change in the fitness function value is less than the tolerance fixed to 1 × 10 −6 . If the termination conditions are not met, elite individuals are automatically survived to the next generation. Besides, the crossover children and mutation children are created based on an established GA Toolbox of MATLAB.
Step 4. Population re-evaluation: Re-evaluate the fitness of the new population. Go to Step 2 and Step 3 repetitively until termination conditions are met.

Test Example of the Proposed Work Flow
A plate with a circular central hole subjected to surface tractions is a classic example for shakedown analysis. By comparing results we obtained from this example to literatures, validity of our in-house program for shakedown analysis is confirmed in our previous studies [29,30]. Here, this example is used again to test the newly established parametric optimization process.   provided in Table 1. Due to symmetry, only 1/4 of the geometry is considered, and the structure is only prescribed with P x in order to be simple. The settings of GA used in this example can be summarized as follows: the maximal number of generation and stall generations are set to 200 and 50, respectively. The elite count is 5% of the population size. The crossover fraction is 80% and the uniform mutation rate is 10%. Initial samples are randomly distributed within the feasible range. When the model is prescribed with P x , distribution of the purely elastic stress σ E under P x are similar to the one shown in Figure 2b. The iterative optimization process is summarized as Figure 3. According to this figure, after evolved for 135 generations and made 6750 function evaluations, the optimality condition is eventually met and the average change in the fitness value drops below 10 −6 . Optimal design parameters obtained are L opt 1 = 50.00 mm and L opt 2 = 43.50 mm which corresponds to the shakedown limit of P x =178.30 MPa. One can notice from Figure 3 that, in the first 70 generations, mean fitness values increase steadily; while after that this value experiences only a slight change and the dominant tendency is that the interval between minimal and maximal fitness becomes smaller. This phenomenon indicates that, if the accuracy of the solution is not particularly emphasized, the optimality condition can be relaxed to save many iterations and function evaluations.
To examine whether or not the result obtained corresponds to the global minimum, the interval for each variable r i is evenly meshed to ten grids, and by calculating α at all 10 × 10 =100 samples, a response surface is constructed (Figure 4). Because in this benchmark model, both the geometry and load conditions are quite simple, the surface appears to be smooth, monotonic and convex. Taking this advantage the parametric optimization process is performed again using the gradient-based  optimizer IPM for the outer loop, and the optimality condition remains the same as GA. Due to the fact that giving the same return value to proximity will result in a zero derivative, which prevents the IPM from continuing the calculation, the IPM here does not utilize the above mentioned pragmatic database based approach. Figure 4 compares trial solutions evaluated by GA and IPM during the entire optimization process, while in Table 2 more detailed information for this comparative study can be found.
Analyzing results illustrated in Figure 4 and Table 2, following conclusions can be drawn: first, variety of samples evaluated by GA is significantly greater than IPM. This is because GA emphasizes on exploring the entire design domain and it prevents the iterations from quickly converging to a local minimum. IPM, in contrast, generates an iterative step only based on local information, and since α is a smooth convex function of r , the algorithms require only few iterations to converge. Second, although optimal solutions r opt derived from GA and IPM are considerably different, the corresponding shakedown limits α(r) are alike and they differ only slightly to the theoretical optimum where [ L 1 L 2 ]=[50.0 50.0] and the shakedown limit is 179.0 MPa. This means in this example, both algorithms are capable of finding optimal design parameters. Third, although the solution obtained by GA and IPM are equivalently good in terms of shakedown limit, GA requires a much greater number of function evaluations (6750/1094 function evaluations without/with pragmatic database-based approach) compared to IPM (290 function evaluations). This is still because the response surface for the present case study is quite particular, and generally such surface will not be so smooth and demonstrates a sound convexity. For those cases, one can only rely on non-gradient algorithms like GA to search for optimal design parameters. To summarize, using the classic plate-with-a-hole model as a benchmark, the established numerical tool for solving shakedown-based optimization problems is critically examined, and it is justified to conclude that the method can be used to conduct parametric optimizations with respect to the shakedown limit.

Description of the Problem
In Figure 5a, a manned airtight module belonging to an international space station is illustrated. The module generally consists of wall panels, connected frames, scuttles and other components where the wall panels and connected frames contribute most of the weight. The connected frame structure, which can be seen from Figure 5b, is not only used to support the entire structure but also designed for realizing some additional functions such as installing other equipment. As a consequence, the structural geometries might cause inhomogeneous stress distribution which leads to critical stress concentrations. Our previous research [20] clarified for this manned airtight module how design parameters influence both plastic and shakedown limits. In the present study, it will be further investigated how to optimize design parameters so as to achieve a maximal strength-to-weight ratio. The structure is exposed to an internal pressure 0.15 MPa, and after performing elastic analysis it is noticed that the critical region lies in the lower corner highlighted in Figure 5b. From this region, two dimensions r 1 and r 2 as depicted in Figure 5b are identified as design parameters, while all other dimensions are fixed to constant values. In the original design, r 1 =7 mm and r 2 =8 mm, and for the optimization problem, range of these parameters are configured as 5 ≤ r 1 ≤ 11 and 6 ≤ r 2 ≤ 12, respectively. Material of the manned airtight module is the aluminum alloy with properties provided in Table 3. To perform the parametric optimization, the airtight module is created as a parameter driven model discretized to quadratic axisymmetric elements CAX4. When parameters change within given intervals, number of elements vary between 45353-46456 while the nodes between 51224-52319.

Results and Discussion
Same as in the plate-with-a-hole example, interval for each design parameter is meshed evenly to 13 grids, and using elastic, plastic and shakedown load factors calculated from 13 × 13=169 grids three response surfaces as shown in Figure 6 are constructed. Following are some noteworthy phenomena one can notice from these surfaces: First, both the elastic and the shakedown surfaces are quite uneven. They contain numerous local peaks and valleys, which means gradient-based algorithms such as IPM are unsuitable. The complexity of these surfaces also discourages the use of surrogate models because a large amount of details would be lost. Second, surface corresponding to elastic limit is far below the plastic and shakedown limits, which means under this circumstance the load bearing capacity of the material is not exhaustively used. Third, when r 1 < 8, r 2 < 7.5, increasing r 1 and r 2 will lead to a greater structural strength. However, once r 1 > 8 and r 2 > 7.5, although the plastic limit increases and saturates, the shakedown limit may even decrease and the surface become uneven.
To better understand the optimization process of GA, samples for several selected generations are presented in the Figure 7. There is a clear trend that the distribution gradually becomes narrower. All 50 samples almost overlap with each other after evolved for 20 generations, and the sample with best fitness values survive to the final generation. Detailed information of optimization runs are exhibited in Table 4. In this table, besides various strengths, one can also find results to the optimization problem employing the strength-to-weight ratio η as the design objective. The convergence process of η can be observed from Figure 8. One can notice from this figure that it requires 52 generations in total for η to achieve a converged solution. The gap between minimal and maximal fitness narrowed rapidly before generation 12, and after 34th generation it becomes almost zero. The individual corresponding to the global optimum first appears at the generation 18, and after generation 35 the distribution of samples are almost unvaried-since then until the final 52nd generation only slight changes in design variables were made.
As can be seen from the Table 4, structure under the original design parameters can already meet the Table 3 Material properties used in manned airtight module of Section 5 ρ (g/cm 3  However, whether to achieve a greater elastic limit or a greater shakedown limit, the two parameters would be adjusted to increase parameter a and decrease parameter b simultaneously compared to the original design. On the other hand, the optimal results under plastic limit design almost reaching the upper bound of the constraint range, it seems that the more material is, the larger the plastic limit load airtight module structure can reach. Unsurprisingly, this trend can also be seen from Figure 6. However, the relationship between strength and weight does not seem to be so straightforward in terms of the shakedown strength, which means there exists a desired balance between shakedown strength and the structural weight. In addition, to attain a more advantageous strength-to-weight performance, the combination of a = 10.22 and b = 7.22 mm should be taken into consideration. Although this will sacrifice a small amount of the installation space in the aerospace vehicle, increasing the wall thickness of the key area is still a reasonable solution.
Since it leads to a decrease of the peak stress presented in the connected structure under the frame as shown in Figure 9. However, what is noteworthy is that there is no clear relationship between von Mises stress and shakedown load result. As the optimal elastic limit design result in the lowest maximum von Mises stress, but the bearing capacity is not the optimal solution.

Conclusions
(1) The current research presents a method to perform design optimization of structures with respect to their shakedown limits. Numerically, the numerical method developed consisted of nested optimization loops. After its validity is intensively tested from the classic plate-with-a-hole model, the method is applied to optimize design parameters pertained to a manned airtight module. Results of this study support the idea that shakedown analysis can be used in placement of the elastic one when designing components in spacecrafts, since it more objectively reflects the actual load capacity of the structure. (2) In addition, this study confirmed that GA is a viable means for determining optimal parameters when shakedown limits are considered as strength criteria. (3) Finally, the present work considers only two design parameters and one objective value, in our future researches the algorithm elaborated will be extend to design problems containing more design parameters and multi-objectives, such that the method can be used to solve problems emerged from real engineering practices.

Formulation of the Lower Bound Problem for Numerical Simulation
We consider a structure discretized into NE elements with a total of NK nodes and NG Gaussian points. Each element has NGK nodes and NGE Gaussian points. Then, the following Eq. (7) can be derived by noticing that the virtual work contribution does not contain self-equilibrated residual stress ρ related part, i.e., ρ is not doing any virtual work.
Here, w ij is referred to as the weight factor of Gaussian point j in element i, J is the Jacobian matrix, B is dubbed the strain-displacement matrix, δu is a variation of virtual displacement and C is a self-equilibrium matrix. Additionally, following Akoa's strategy [11] of tailoring to von Mises criterion, one can let where Then we may rephrase (5) as the following form for two load cases, i.e., NL=2.
It is worth noting that such mathematical form can be extended to the three or more load cases which is usually called high dimension shakedown form. �·� 2 here is the Euclidean norm of a vector, γ r represents a column vector in R 5 whose five components corresponding to the first five rows of S σ e r1 − σ e r2 . Observe the optimization problem (5) and (14), one may draw a conclusion that the rephrased optimization problem is regular, nonlinear and convex. In addition, one may view u r1 and y r1 as a variable substitution of σ r1 : maximize u r1 ,u r2 ,y 1 ,α −α, subject to NG r=1 M r (u r1 ) + N y 1 − αw 1 = 0, u r1 − u r2 = αγ r , �u r1 � 2 ≤ 1, �u r2 � 2 ≤ 1, r= 1,...,NG.
Such variable substitution also serves the tailoring strategy which will lead to a better computational efficiency. Eq. (14) is a typical second order cone programming (SOCP) problem which can be solved quickly using the interior-point methods (IPM). The problem is then submitted to Gurobi which is a SOCP powerful solver based on IPM to calculate the maximum load factor α and its associated self-equilibrated residual stress ρ.