Issue 
Int. J. Simul. Multisci. Des. Optim.
Volume 7, 2016



Article Number  A3  
Number of page(s)  13  
DOI  https://doi.org/10.1051/smdo/2016002  
Published online  26 February 2016 
Research Article
Stiffness optimization of geometrically nonlinear structures and the level set based solution
The State Key Laboratory of Digital Manufacturing Equipment and Technology, Huazhong University of Science and Technology, Wuhan, P.R. China
^{*} email: qxia@mail.hust.edu.cn
Received:
11
December
2015
Accepted:
8
January
2016
Loadnormalized strain energy increments between consecutive load steps are aggregated through the KreisselmeierSteinhauser (KS) function, and the KS function is proposed as a stiffness criterion of geometrically nonlinear structures. A topology optimization problem is defined to minimize the KS function together with the perimeter of structure and a volume constraint. The finite element analysis is done by remeshing, and artificial weak material is not used. The topology optimization problem is solved by using the level set method. Several numerical examples in two dimensions are provided. Other criteria of stiffness, i.e., the end compliance and the complementary work, are compared.
Key words: Topology optimization / Geometrically nonlinear / Level set method / Stiffness / KreisselmeierSteinhauser function
© Q. Xia and T. Shi, Published by EDP Sciences, 2016
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1 Introduction
Maximizing stiffness of structures is a major topic of structural optimization. However, most of the studies on stiffness maximization consider linear elastic structures, and the studies on geometrically nonlinear structures are limited.
As the magnitude of external load increases slowly from zero to the target load level (i.e., under quasistatic condition), the configuration of a geometrically nonlinear structure changes, as a result the stiffness of structure changes as well. For instance, the distance from the external load to the support of a structure may become shorter when the structure deforms, thus the stiffness of structure increases. One can see that the stiffness of a geometrically nonlinear structure depends on its configuration. Therefore, evaluation of the stiffness of a geometrically nonlinear structure is ambiguous. It may refer to an instantaneous stiffness of a specific configuration in the whole course of deformation, or it may have a more comprehensive meaning that involves with many configurations.
Stiffness of a structure is often represented by the compliance (denoted as W), the strain energy (denoted as U), or the complementary work (denoted as U_{C}) [1–4]. For a linear elastic structure, minimizing any one of the three functions leads to the same result because W = 2U = 2U_{C}. For a geometrically nonlinear structure, however, different results may be obtained if a different function is used in optimization because W ≠ 2U ≠ 2U_{C} [2]. The compliance, the strain energy, the complementary work, along with a typical loaddisplacement curve of a geometrically nonlinear structure are illustrated in Figure 1. Area of the rectangle is equal to the compliance; area of the region under or above the curve is respectively equal to the strain energy or the complementary work.
Figure 1. A loaddisplacement curve of geometrically nonlinear structure. 
Figure 2. (a) Design problem and (b) initial design. 
Compliance or the socalled end compliance is used in most of the studies that aim to maximize the stiffness of geometrically nonlinear structures. It characterizes the stiffness of the final equilibrium configuration of a geometrically nonlinear structure, and the result of compliance minimization depends on load level [2]. In other words, the shape and topology of the optimal structure will change if the magnitude of design load is changed. More interesting, it was observed that the optimal structure obtained with a big design load is not stiff for smaller design load [2]. Therefore, it was suggested that the end compliance may not be the best criterion of stiffness for geometrically nonlinear structures [2, 5].
In order to overcome such a problem, Buhl et al. [2] proposed to minimize the complementary work and numerically approximated it through a weighted combination of displacements of many deformed configurations including the final equilibrium configuration and those went through before the final equilibrium. The strain energy was also used in the optimization of geometrically nonlinear structures [6]. Again, it was numerically approximated through a weighted combination of displacements of many deformed configurations [6]. Minimizing the complementary work or the strain energy can be regarded as a multiobjective optimization whose objectives are the displacements of many deformed configurations. To conclude, the end compliance describes only the instantaneous stiffness of the equilibrium configuration, but the complementary work and the strain energy give more comprehensive descriptions of stiffness since they involve with many configurations.
Although minimization of the complementary work or the strain energy offers a control over the displacements of many deformed configurations, the fixed weights of displacements in the weighted combinations are not perfect when they are viewed in the perspective of multiobjective optimization. Therefore, in the present study, we propose to extend the inspiring previous studies by considering loadnormalized strain energy increments between consecutive load steps. The increments are aggregated through the KreisselmeierSteinhauser (KS) function [7], and the KS function is proposed as a stiffness criterion of geometrically nonlinear structures. Then, a topology optimization problem is defined. The different criteria of stiffness are investigated and compared in the present paper.
Another important issue in the topology optimization of geometrically nonlinear structures is concerned with the numerical instability that may arise in the finite element analysis. In most topology optimizations, a method employing a fixed mesh and ersatz material is used for the finite element analysis. Such a method uses an artificial weak material to represent void. However, in the topology optimization of geometrically nonlinear structures, if an element filled with the artificial weak material undergoes excessively large deformation, its tangent stiffness matrix may become negative, and this will cause numerical instability. In the literature, several methods have been proposed to deal with such numerical instability. Buhl et al. [2] observed that the numerical instability arised at the nodes surrounded by “void” elements and proposed to eliminate such nodes from the convergence criterion. Bruns and Tortorelli [8] presented a systematic method for removing and reintroducing low density elements. Yoon and Kim [9] proposed a method that uses connectivity between elements as the design variables of topology optimization. In addition, Lahuerta et al. [10] and Wang et al. [11] proposed to tailor the material model of the artificial weak material such that the numerical instability can be eliminated. Meshless methods were also proposed by Cho and Kwak [12] and He et al. [13] to overcome the numerical instability. In a level set based solution to the topology optimization of geometrically nonlinear structures, remesh is proposed by Ha and Cho [14] to obtain a mesh that is conformal to the geometry of the structure, thus eliminating the artificial weak material as well as the numerical instability coming with it. In the present study, being similar to the method proposed by Ha and Cho [14], we do the finite element analysis by modifying a fixed background mesh and do not use the artificial weak material, as proposed in our previous study [15].
The paper is organized as follows. In Section 2 some basics of geometrically nonlinear structure are outlined. In Section 3 the end compliance, the strain energy, the complementary work, and the proposed maximum of loadnormalized strain energy increments are discussed; formulation of the topology optimization problem is described. In Section 4 sensitivity analysis for the optimization problem is described. In Section 5 the level set method is briefly described. Section 6 describes the finite element analysis. Section 7 gives numerical examples and discussions. Section 8 concludes this paper.
2 Analysis of geometrically nonlinear structure
The open bounded set occupied by a structure is denoted as (d = 2 or 3). It is required that during optimization all admissible structures stay in a fixed reference domain , i.e., . The boundary of Ω consists of three disjoint segments, i.e.,(1)where an external load is applied on the Neumann boundary Γ_{N}; a displacement constraint is enforced on the Dirichlet boundary Γ_{D}; Γ_{H} denotes the free boundary. In the present study, only Γ_{H} is subjected to optimization, while both Γ_{N} and Γ_{D} are fixed during the optimization.
In total Lagrangian formulation, weak form of the equilibrium equation of a geometrically nonlinear structure is given by(2)where(3) (4)where ^{0}Ω and ^{0}Γ_{N} are the structure domain and the Neumann boundary at initial undeformed configuration (a left superscript on Ω and Γ_{N} denotes the configuration of structure); u is displacement; is virtual displacement; is the space of kinematically admissible displacement fields, i.e., ; S_{ij}(u) is the second PiolaKirchhoff stress tensor; is the virtual GreenLagrange strain tensor; t_{i} is surface traction. In the present study, surface traction are assumed to be independent of deformation (e.g., a point load that does not change its direction during deformation of the structure), and body force is not considered. Details of the analysis are referred to [16, 17].
In the present study, although a structure undergoes large displacement, the strain is assumed to be small. In such case, the St.VenantKirchhoff material model is used, and the constitutive relation is given by(5)where C_{ijkl} is the elasticity tensor, and in the present study the material is assumed to be isotropic. The GreenLagrange strain ε_{kl}(u) is given by(6)the virtual GreenLagrange strain is given by(7)
Because is nonlinear in the argument u, the NewtonRaphson method is used to solve the nonlinear equilibrium equation, i.e., equation (2), in an iterative manner. The displacement ^{m+1}u of the (m + 1)th iteration is computed as the sum of ^{m}u and an increment of displacement Δu (8)
The increment of displacement Δu is obtained by solving a linearized equilibrium equation given by(9)where(10) (11)and(12)
The residual in equation (11) defines the difference between and , i.e.,(13)
The equilibrium of a structure is found when R becomes a zero vector.
3 Criteria of stiffness and optimization problem
In order to quantify the stiffness of a geometrically nonlinear structure, several criteria were used and discussed in the literature. They include the end compliance, the strain energy, and the complementary work [1–4]. The relation between the three criteria is given by(14)where W is the end compliance; U is the strain energy; U_{C} is the complementary work.
The end compliance is used in most of the studies on stiffness optimization of geometrically nonlinear structures. It is defined by(15)
The last equality is obtained by substituting u for in equation (2). The form of end compliance given by equation (15) is the same as that of the virtual work of external loads (i.e., Eq. (4)), however, the end compliance is not the true work of external loads (under quasistatic condition).
The true work of external loads under quasistatic condition is equal to the strain energy stored in a structure, and it is given by(16)
It can be approximated by using a trapezoidal scheme as [1, 6](17)where t ^{s} and u ^{s} denote respectively the surface traction and displacement of the sth load step (here we use a right superscript to denote the index of load step in the numerical solution of equilibrium equation); the displacement or u ^{s} is the converged solution obtained through the NewtonRaphson method for the sth load step ; ΔU ^{s} denotes the increment of strain energy from the load step s − 1 to s.
Assuming that the increments of external load between consecutive load steps are the same, we use Δt_{i} to denote the increment. In the present study, because only dead load is considered, such an assumption is valid, but it will not be valid for follower load. Based on this assumption, the strain energy can be rewritten as(18)
The complementary work U ^{C} of external load is obtained by substituting equations (15) and (16) into equation (14) as [3](19)
It can also be approximated by using a trapezoidal scheme as [2](20)where denotes the increment of complementary work from the load step s−1 to s. From equation (19), one can see that the complementary work U_{C} of a geometrically nonlinear structure is not equal to the strain energy U.
As can be seen in equations (18) and (20), the strain energy and the complementary work can be regarded as a weighted combination of displacements. For instance, Δt_{i}/2 or Δt_{i} are the weights in the complementary work. Minimization of U or U_{C} offers a control over the displacements of intermediate load steps or the stiffness of intermediate configurations. Therefore, as a stiffness criterion of geometrically nonlinear structures, the strain energy and the complementary work give a more comprehensive description of stiffness than the end compliance. However, the fixed weights of displacements in the weighted combinations are not perfect when they are viewed in the perspective of multiobjective optimization.
In the present study, loadnormalized strain energy increments between consecutive load steps are aggregated through the KreisselmeierSteinhauser (KS) function [7], and the KS function is proposed as a stiffness criterion of geometrically nonlinear structures. Here, the KS function is defined as(21)where ρ > 0 is a parameter, and Δu ^{s} is defined by(22)where l_{i}, i = 1, … d are the weighting parameters, and they are regarded as the components of an unit fictitious load whose direction is the same as that of the external load t. Because only dead load (magnitude and direction of the load are independent of displacement) is considered in the present study, the weighting parameters l_{i} are fixed during the optimization.
Comparing equation (22) with equation (17), one can see that Δu^{s} can be regarded as a loadnormalized strain energy increment, i.e., the increment of strain energy ΔU^{s} is normalized by the average load . In addition, suppose that a point load along the negative direction of yaxis is applied to a geometrically nonlinear structure. Then, in equation (22), the components of fictitious load are set to l_{1} = 0 and l_{2} = −1. In this case, Δu^{s} can also be regarded as the increment of displacement at the traction boundary.
The KS function has been used in many optimization problems to aggregate a family of functions into a single function [18–21]. It was proved that the KS function is always bigger than the maximum of the aggregated functions and that the limit of KS function will be equal to the maximum when the parameter ρ approaches positive infinity [22], i.e.,(23)and(24)
Therefore, when the parameter ρ is very big, the KS function can be regarded as the maximum of loadnormalized strain energy increments.
Now, the optimization problem is defined as(25)where V(Ω) = ∫_{Ω} dx is the volume of a structure; α in the objective function is a fixed positive parameter; P(Ω) is the perimeter of a structure that is used to regularize the optimization problem [23–26](26)
The detailed discussions about the perimeterbased regularization in a level set based structure optimization are referred to our previous studies [25, 26] and the references therein.
4 Sensitivity analysis
4.1 Sensitivity of KS (Δu)
We consider first the sensitivity of a generic function of displacement given by(27)where w _{ i } is the weight of displacement; u _{ i } is the displacement of an equilibrium configuration (equilibrium configuration of an intermediate load step or the final equilibrium configuration). For the end compliance, the weight is the external load t _{ i }. For the strain energy and the complementary work, the weights can be inferred from equation (18) to equation (20). For the Δu ^{ s } in equation (22) the weights are the fictitious load l _{ i }.
The Lagrangian concerned with the generic function ψ(u) is defined as(28)where v is a Lagrange multiplier for equation (2), and it is regarded as the adjoint displacement.
The partial derivative of with respect to u in the direction is(29)where(30) (31)
The stationary condition of equation (29) leads to the adjoint equation(32)
Equation (32) is solved for the adjoint displacement v, and the discretized version of is the tangent stiffness matrix constructed when the NewtonRaphson algorithm arrives at a converged solution u.
Ignore the shape derivatives defined on boundary ^{0}Γ_{ N } since ^{0}Γ_{ N } is not subject to optimization, and note that only ^{0}Γ_{ H } is optimized, we have the partial derivative of with respect to ^{0}Ω in the direction θ (33)and shape derivative of the function ψ(u) is(34)
Through simple derivation, the shape derivative of strain energy can be obtained as(35)where(36)and u ^{ s } and v ^{ s } are the primary and adjoint displacements of the sth load step. The shape derivative of complementary work has the same form as equation (35), except that the adjoint displacements v ^{ s } are the solutions to different adjoint equations.
The shape derivative of the KS function equation (21) is(37)where(38) (39)where the adjoint displacements v ^{ s } in G ^{ s } and v ^{ s−1} in G ^{ s−1} are the solutions of the following two adjoint equations(40) (41)
The final result of shape derivative is rewritten as(42)
4.2 Sensitivity of the perimeter term and volume constraint
Now, we consider the objective function J augmented with the volume constraint(43)where λ is the Lagrange multiplier for the volume constraint, and it is updated during the optimization according to the the augmented Lagrange multiplier method [27].
Note that only ^{0}Γ_{ H } is optimized, the shape derivative of is(44)where(45)where κ is the curvature of boundary.
According to equation (45), we can readily obtain the steepest descent direction θ by setting(46)
It can be seen that such θ will yield which implies the descent of .
The properties of the optimization problem and the optimization algorithm in the present work do not guarantee global optima. The optimization problem in the present study is not convex, and there are many solutions such as one global and many local minima. The optimization algorithm, i.e., the steepest descent, is based on the gradient at the current point in the design space, thus it sees no global information about the objective function and constraints, and all improvements of the structure is made based on local information at the current point in the design space. Therefore, it cannot perform a global search [28].
5 The level set method
In the level set method [29, 30], the boundary of a shape Ω is represented implicitly through a level set function Φ(x) as its zero level set, i.e., , and Φ(x) is used to define the inside and outside regions with respect to the boundary as follows(47)
In the present study, Φ is constructed as the signed distance function to the free boundary ^{0}Γ_{H} of a structure.
Propagation of ^{0}Γ_{H} is modeled by the HamiltonJacobi equation(48)where(49)
A structured grid and the finite difference method (the first order upwind spatial differencing and forward Euler time differencing) are used to solve the equation, and a reinitialization procedure is periodically performed. A PDE based method is used for velocity extension [31]. In the velocity extension, we first compute a and b on Γ_{H}, then we extend these quantities to the entire reference domain . More details of the numerical computations of level set can be found in [29, 30].
A move limit strategy [32] is applied to the extended velocity to improve the stability and convergence rate of the level set based optimization. This is achieved by applying the following equation to the extended velocity(50)where is the original extended velocity; is the extended velocity after applying the move limit; h is the grid size of level set computation, γ > 0 is a parameter that controls the size of move limit. The term Φ(x)/h stands for the normalization of Φ(x). The normalization is introduced here so that the value of γ is independent of the grid size h. When the value of γ decreases, the size of move limit decreases. In our present work, γ is set to 0.1.
6 Finite element analysis
Conventionally, an Eulerian method employing a fixed mesh and ersatz material is used for the finite element analysis in the level set based method of topology optimization. Instead of solving the equilibrium equation on Ω, it solves the equation on the entire reference domain with the void being mimicked by an artificial weak material. However, in the present study of topology optimization of geometrically nonlinear structures, if the elements filled by the artificial weak material undergo excessively large deformation, their tangent stiffness matrix may become negative, and this will cause numerical instability.
In the present work, being similar to the method proposed by Ha and Cho [14], we do the finite element analysis by modifying a fixed background mesh and do not use the artificial weak material, as proposed in our previous study [15] and used in [25, 33]. A similar method was also proposed by Allaire et al. [34]. Before the optimization, we set up a background triangle mesh for the reference domain . In each iteration of the optimization, modification of the fixed background mesh is performed to obtain a mesh that is conformal to the geometry of the structure. The detail of the modification is referred to [15].
7 Numerical examples
In this section, several examples in two dimensions are given, although there exists no conceptual difficulty to do the numerical examples in three dimensions. In these examples, it is assumed that the solid material has a Young’s modulus E = 3 GPa and Poisson’s ratio . For all the examples, the plane stress state is assumed; thickness is 0.1 m; 3node triangle finite element is used.
7.1 Example 1
The optimal design problem of the first example is shown in Figure 1. The reference domain is a rectangle of size 1 m × 0.25 m. The left side is fixed. The upper bound of volume is 50% of the reference domain. There are 12,318 triangle elements in the background mesh for the reference domain . A 500 × 125 rectilinear grid is used for level set computations. The initial design is shown in Figure 1.
We set load to 100, 200, 300, and 500 kN respectively in five optimizations, and in each optimization the load is evenly divided into five load steps. The parameter in the KS function is set to ρ = 20, and the parameter for the perimeter term in the objective function is set to α = 1 × 10^{−4} because the order of magnitude of KS function is 10^{−1}. The optimal structures are shown in Figure 3. The history of convergence with the load in optimization being 100 kN is provided as an example in Figure 4.
Figure 3. Optimal structures with the KS function being used as the stiffness criterion. (a) t = 100 kN, (b) t = 200 kN, (c) t = 300 kN, (d) t = 500 kN. 
Figure 4. History of convergence. 
In order to compare the proposed KS function with the end compliance and the complementary work, optimizations are also done by replacing the KS function in the objective function with W or U _{ C }. In these optimizations, the parameter for the perimeter term in the objective function is set to α = 1 because the order of magnitude of W and U _{ C } is 10^{4}–10^{5}. The optimal structures obtained when the end compliance is used as the stiffness criterion are shown in Figure 5, and those of the complementary work are shown in Figure 6. The deformed shape of the optimal structures are shown in Figure 7. Figures 7a, d, g, j are the deformed optimal structures obtained when the KS function is chosen as the stiffness criterion; Figure 7b, e, h, k are of the end compliance; Figure 7c, f, i, l are of the complementary work.
Figure 5. Optimal structures with W being used as the stiffness criterion. (a) t = 100 kN, (b) t = 200 kN, (c) t = 300 kN, (d) t = 500 kN. 
Figure 6. Optimal structures with UC being used as the stiffness criterion. (a) t = 100 kN, (b) t = 200 kN, (c) t = 300 kN, (d) t = 500 kN. 
Figure 7. Deformed optimal structures of different stiffness criteria. (a) KS (t = 100 kN), (b) U (t = 100 kN), (c) U _{ C } (t = 100 kN), (d) KS (t = 200 kN), (e) U (t = 200 kN), (f) U _{ C } (t = 200 kN), (g) KS (t = 300 kN), (h) U (t = 300 kN), (i) U _{ C } (t = 300 kN), (j) KS (t = 500 kN), (k) U (t = 500 kN), (l) U _{ C } (t = 500 kN). 
From these results one can see that the shape and topology of all the optimal structures depend on the load level. In contrast, for a linear elastic structure, the results will be independent of load level due to the principle of superposition [5].
When the load is 500 kN, all the optimal structures have a slender beam at the right end except that their lengths are different. Such a slender beam is not stiff, but it is able to effectively reduce the length of load arm in the final equilibrium configuration, as shown in Figure 7, thus effectively increasing the stiffness of the final equilibrium configuration. Such a phenomenon is most prominent in the optimal structures obtained when the end compliance is chosen as the stiffness criterion, since the end compliance concerns only with the stiffness of the final equilibrium configuration. However, the slender beam at the right end, for instance the one shown in Figure 5d, is not efficient for smaller loads, for instance an 100 kN load. One can imagine that even if the slender beam in Figure 5d is changed to connect with the remaining structure through a hinge, the end compliance will not increase significantly.
Detailed quantitative comparisons are given in Table 1. As one can see that when the end compliance is chosen as the stiffness criterion in optimization, the corresponding optimal structures have the smallest end compliance and the smallest strain energy (except when load is 200 kN), and that such a trend is more prominent when the load level is high. But as explained above, this phenomenon only indicates that the end compliance is a good stiffness criterion for the final equilibrium configuration. Actually, it was suggested that the end compliance may not be the best criterion of stiffness for geometrically nonlinear structures [2, 5].
Quantitative results of optimization with different criteria of stiffness. (Although there is no parameter ρ in the optimization of end compliance and complementary work, ρ is set to 20 for the comparison. The volume of all the optimal structures are 50% of the reference domain. The unit of W, U, U_{C} is J.)
In order to overcome such a problem, Buhl et al. [2] proposed to minimize the complementary work. Minimizing the complementary work can be regarded as a multiobjective optimization whose objectives are the displacements of many deformed configurations. As one can see in the results shown in Figure 6, when the load is 200 kN or 300 kN there is no slender beam at the right end, and when the load is 500 kN the length of the slender beam is much shorter than that shown in Figure 5d. In addition, as can be seen in Table 1, when U_{C} is chosen as the stiffness criterion, the optimal structures have the smallest U_{C} (except when load is 200 kN).
The inspiring previous studies are extended in our present work by considering the aggregated loadnormalized strain energy increments between consecutive load steps. From the optimal structure shown in Figure 3d, one can see that when the load is 500 kN the length of slender beam is much shorter than that shown in Figure 5d and is comparable to that shown in Figure 6d. More importantly, the quantitative results given in Table 1 show that when the KS function is chosen as the stiffness criterion in optimization, the optimal structures have smaller end compliance and smaller strain energy as compared to the optimal structures obtained when the complimentary work is chosen as the stiffness criterion (except when load is 100 kN).
Comparison of loaddisplacement curves of the optimal structures obtained when load is 500 kN in the optimization are shown in Figure 8. As can be seen in this figure, when the load level is below 300 kN, the displacement of the optimal structure obtained when the end compliance is chosen as the stiffness criterion is the largest, which indicates that the slender beam at the right end shown in Figure 5d is not stiff. The red curve is almost the same as the blue curve, nevertheless Δu^{1} of the red curve is smaller than that of the blue curve.
Figure 8. Loaddisplacement curves of the optimal structures (load is 500 kN). 
Figure 9. (a) Design problem and (b) initial design of the second example 
As mentioned in Section 3, when the parameter ρ in the KS function approaches positive infinity, the limit of the KS function will be equal to the maximum of Δu^{s}. In this example, we observed that Δu^{1} is always the maximal one of Δu^{s} because when the structure deforms the length of load arm decreases and hence stiffness increases. Therefore, a bigger value of ρ in this example will let Δu^{1} get a more significant role in the optimization. It seems that a higher value of ρ can eliminate the slender beam at the right end. However, we also observed that when the optimal structure obtained with an 100 kN load is subjected to a 500 kN load, the structure buckles and is not stable. Therefore, in the optimization when load is 500 kN, Δu^{1} should not dominate the optimization (recall that Δu^{1} is the load normalized strain energy increment when load is 100 kN), and the value of ρ should not be too big. In the present example, we see that ρ = 20 is a good choice. For the structure shown in Figure 3d, the Δu^{s} are: Δu^{1} = 0.1345, Δu^{2} = 0.1033, Δu^{3} = 0.0907, Δu^{4} = 0.0800, Δu^{5} = 0.0709, and the KS function is 0.1829. In fact, the parameter ρ in the KS function offers to designer a way to adjust the weights of each Δu^{s}, but the weights in the complementary work are fixed.
7.2 Example 2
The optimal design problem of the second example is shown in Figure 8. The reference domain is a rectangle of size 2 m × 1 m. The upper bound of volume is 30% of the reference domain. There are 10,934 triangle elements in the background mesh for the reference domain . A 600 × 300 rectilinear grid is used for level set computations. The initial design is shown in Figure 8. The load is 1200 kN, and it is evenly divided into five load steps. The parameter in the KS function is set to ρ = 400. The value of ρ is much bigger than that of the first example because the Δu^{s} in this example are closer to zero, and we observed that when the values of the aggregated functions are closer to zero, ρ should be bigger to ensure the KS function reasonably approximate the maximum. The parameter for the perimeter term in the objective function is set to α = 1 × 10^{−4}.
The optimal structure is shown in Figure 10. The history of convergence is shown in Figure 11, from which one can see that the convergence is smooth. Four intermediate designs of the optimization are shown in Figure 12, from which one can see that the shape and topology were changed during the optimization. The end compliance is W = 67,823 J; the strain energy is U = 34,673 J; the complementary work is U_{C} = 33,150 J; the KS function is KS = 0.015. The Δu^{s} are: Δu^{1} = 0.0108, Δu^{2} = 0.0109, Δu^{3} = 0.0111, Δu^{4} = 0.0115, Δu^{5} = 0.0121; the volume ratio is 30%. From the values of Δu^{s}, one can see that the optimal structure is softening as the load increases, which is different from the hardening behavior in the first example.
Figure 10. The optimal structure of the second example. 
Figure 11. History of convergence of the second example. 
Figure 12. Intermediate designs of the optimization of the second example. (a) Step 55, (b) step 90, (c) step 110, (d) step 200. 
Figure 13. (a) Design problem and (b) initial design of the third example. 
If the structure is linear elastic, a wellknown Michell’s solution to the current design problem is available in the literature [35–37]. The shape and topology of the optimal geometrically nonlinear structure is almost the same as the well known solution.
7.3 Example 3
The optimal design problem of the third example is shown in Figure 12. The upper bound of volume is 40% of the reference domain. There are 10,934 triangle elements in the background mesh for the reference domain . A 600 × 300 rectilinear grid is used for level set computations. The initial design is shown in Figure 12. The load is 1500 kN, and it is evenly divided into five load steps. The parameter in the KS function is set to ρ = 400. The parameter for the perimeter term in the objective function is set to α = 1 × 10^{−4}.
The optimal structure is shown in Figure 14. The history of convergence is shown in Figure 15, from which one can see that the convergence is smooth. The end compliance is W = 124,948 J; the strain energy is U = 64,229 J; the complementary work is U_{C} = 60,719 J; the KS function is KS = 0.021. The Δu^{s} are: Δu^{1} = 0.0157, Δu^{2} = 0.0160, Δu^{3} = 0.0164, Δu^{4} = 0.0171, Δu^{5} = 0.0181; the volume ratio is 40%. From the values of Δu^{s}, one can see that the optimal structure is softening as the load increases.
Figure 14. The optimal structure of the third example. 
Figure 15. History of convergence of the third example. 
8 Conclusion
In the present paper, the aggregated loadnormalized strain energy increments between consecutive load steps is proposed as a stiffness criterion of geometrically nonlinear structures. A topology optimization problem is then defined to minimize the KS function together with the perimeter of structure and a volume constraint. The finite element analysis is done by remeshing, and artificial weak material is not used. Therefore, the numerical instability related to the artificial weak material is eliminated. The shape derivative of free boundary is derived by using a Lagrangian function and the adjoint method. The topology optimization problem is solved by using the level set method. Several numerical examples in two dimensions are provided. Other criteria of stiffness, i.e., the end compliance and the complementary work, are compared. The results shown that the proposed KS function gives a comprehensive description of the stiffness of geometrically nonlinear structures.
Acknowledgments
This research work is supported by the National Natural Science Foundation of China (Grant No. 51575203).
References
 Yuge K, Iwai N, Kikuchi N. 1999. Optimization of 2D structures subjected to nonlinear deformation using the homogenization method. Struct. Optim., 17, 286–299. [CrossRef] [Google Scholar]
 Buhl T, Pedersen CBW, Sigmund O. 2000. Stiffness design of geometrically nonlinear structures using topology optimization. Struct. Multidiscip. Optim., 19, 93–104. [CrossRef] [Google Scholar]
 Gea HC, Luo J. 2001. Topology optimization of structures with geometrical nonlinearities. Comput. Struct., 79, 1977–1985. [CrossRef] [Google Scholar]
 Kemmler R, Lipka A, Ramm E. 2005. Large deformations and stability in topology optimization. Struct. Multidiscip. Optim., 30, 459–476. [CrossRef] [MathSciNet] [Google Scholar]
 Bruns TE, Tortorelli DA. 2001. Topology optimization of nonlinear elastic structures and compliant mechanisms. Comput. Methods Appl. Mech. Eng., 190(26–27), 3443–3459. [CrossRef] [Google Scholar]
 Huang X, Xie XM. 2008. Topology optimization of nonlinear structures under displacement loading. Eng. Struct., 30, 2057–2068. [CrossRef] [Google Scholar]
 Kreisselmeier G, Steinhauser R. 1979. Systematic control design by optimizing a vector performance index, International Federation of Active Controls Symposium on ComputerAided Design of Control Systems, Zurich, Switzerland. [Google Scholar]
 Bruns TE, Tortorelli DA. 2003. An element removal and reintroduction strategy for the topology optimization of structures and compliant mechanisms. Int. J. Numer. Meth. Eng., 57, 1413–1430. [CrossRef] [Google Scholar]
 Yoon GH, Kim YY. 2005. Element connectivity parameterization for topology optimization of geometrically nonlinear structures. Int. J. Solids Struct., 42, 1983–2009. [CrossRef] [Google Scholar]
 Lahuerta RD, Oes ETS, Campello EMB, Pimenta PM, Silva ECN. 2013. Towards the stabilization of the low density elements in topology optimization with large deformation. Comput. Mech., 52, 779–797. [CrossRef] [MathSciNet] [Google Scholar]
 Wang F, Lazarov BS, Sigmund O, Jensen JS. 2014. Interpolation scheme for fictitious domain techniques and topology optimization of finite strain elastic problems. Comput. Methods Appl. Mech. Eng., 276, 453–472. [CrossRef] [MathSciNet] [Google Scholar]
 Cho S, Kwak J. 2006. Topology design optimization of geometrically nonlinear structures using meshfree method. Comput. Methods Appl. Mech. Eng., 195, 5909–5925. [CrossRef] [MathSciNet] [Google Scholar]
 He Q, Kang Z, Wang YQ. 2014. A topology optimization method for geometrically nonlinear structures with meshless analysis and independent density field interpolation. Comput. Mech., 54, 629–644. [CrossRef] [MathSciNet] [Google Scholar]
 Ha SH, Cho S. 2008. Level set based topological shape optimization of geometrically nonlinear structures using unstructured mesh. Comput. Struct., 86, 1447–1455. [CrossRef] [Google Scholar]
 Xia Q, Shi T, Liu S, Wang MY. 2012. A level set solution to the stressbased structural shape and topology optimization. Comput. Struct., 90–91, 55–64. [CrossRef] [Google Scholar]
 Bathe KJ. 2004. Finite element procedures. Prentice Hall: Upper Saddle River, New Jersey. [Google Scholar]
 Reddy JN. 2004. An introduction to nonlinear finite element analysis. Oxford University Press: Oxford. [CrossRef] [Google Scholar]
 Yang RJ, Chahande AI. 1995. Automotive applications of topology optimization. Struct. Optim., 9, 245–249. [CrossRef] [Google Scholar]
 Akgun MA, Haftka RT, Wu KC, Walsh JL, Garcelon JH. 2001. Efficient structural optimization for multiple load cases using adjoint sensitivities, AIAA J., 39, 511–516. [CrossRef] [Google Scholar]
 Poon NMK, Martins JRRA. 2007. An adaptive approach to constraint aggregation using adjoint sensitivity analysis. Struct. Multidiscip. Optim., 34, 61–73. [CrossRef] [Google Scholar]
 Luo YJ, Wang MY, Kang Z. 2013. An enhanced aggregation method for topology optimization with local stress constraints. Comput. Methods Appl. Mech. Eng., 254, 31–41. [CrossRef] [MathSciNet] [Google Scholar]
 Raspanti CG, Bandoni JA, Biegler LT. 2000. New strategies for flexibility analysis and design under uncertainty. Comput. Chem. Eng., 24, 2193–2209. [CrossRef] [Google Scholar]
 Allaire G, Jouve F, Toader AM. 2004. Structural optimization using sensitivity analysis and a levelset method. J. Comput. Phys., 194, 363–393. [NASA ADS] [CrossRef] [EDP Sciences] [MathSciNet] [PubMed] [Google Scholar]
 Wang MY, Wang XM. 2004. PDEdriven level sets, shape sensitivity and curvature flow for structural topology optimization. CMESComput. Model. Eng. Sci., 6, 373–395. [Google Scholar]
 Xia Q, Wang MY, Shi TL. 2014. A level set method for shape and topology optimization of both structure and support of continuum structures. Comput. Methods Appl. Mech. Eng., 272, 340–353. [CrossRef] [MathSciNet] [Google Scholar]
 Xia Q, Wang MY, Shi TL. 2015. Topology optimization with pressure load through a level set method. Comput. Methods Appl. Mech. Eng., 283, 177–195. [CrossRef] [MathSciNet] [Google Scholar]
 Nocedal J, Wright SJ. 1999. Numerical optimization. Springer: New York. [CrossRef] [MathSciNet] [Google Scholar]
 Rozvany GIN. 2001. Aims, scope, methods, history and unified terminology of computeraided topology optimization in structural mechanics. Struct. Multidiscip. Optim., 21, 90–108. [CrossRef] [Google Scholar]
 Sethian JA. 1999. Level set methods and fast marching methods: evolving interfaces in computational geometry, fluid mechanics, computer vision, and materials science, 2nd edn. Cambridge Monographs on Applied and Computational Mathematics, Cambridge University Press: Cambridge, UK. [Google Scholar]
 Osher S, Fedkiw R. 2002. Level set methods and dynamic implicit surfaces. SpringerVerlag: New York. [Google Scholar]
 Peng DP, Merriman B, Osher S, Zhao HK, Kang M. 1999. A PDEbased fast local level set method. J. Comput. Phys., 155, 410–438. [CrossRef] [MathSciNet] [Google Scholar]
 Xia Q, Wang MY, Shi TL. 2013. A move limit strategy for level set based structural optimization. Eng. Optim., 45, 1061–1072. [CrossRef] [MathSciNet] [Google Scholar]
 Xia Q, Shi T, Liu S, Wang MY. 2013. Shape and topology optimization for tailoring stress in a local region to enhance performance of piezoresistive sensors. Comput. Struct., 114–115, 98–105. [CrossRef] [Google Scholar]
 Allaire G, Dapogny C, Frey P. 2014. Shape optimization with a level set based mesh evolution method. Comput. Methods Appl. Mech. Eng., 282, 22–53. [CrossRef] [MathSciNet] [Google Scholar]
 Michell AGM. 1904. The limits of economy of material in framestructures. Phil. Mag., 8, 589–597. [CrossRef] [Google Scholar]
 Hemp WS. 1973. Michell’s structural continua, in Optimum Structures. Clarendon Press: Oxford, Ch. 4. [Google Scholar]
 Xie YM, Steven GP. 1993. A simple evolutionary procedure for structural optimization. Comput. Struct., 49, 885–896. [CrossRef] [Google Scholar]
Cite this article as: Xia Q & Shi T: Stiffness optimization of geometrically nonlinear structures and the level set based solution. Int. J. Simul. Multisci. Des. Optim., 2016, 7, A3.
All Tables
Quantitative results of optimization with different criteria of stiffness. (Although there is no parameter ρ in the optimization of end compliance and complementary work, ρ is set to 20 for the comparison. The volume of all the optimal structures are 50% of the reference domain. The unit of W, U, U_{C} is J.)
All Figures
Figure 1. A loaddisplacement curve of geometrically nonlinear structure. 

In the text 
Figure 2. (a) Design problem and (b) initial design. 

In the text 
Figure 3. Optimal structures with the KS function being used as the stiffness criterion. (a) t = 100 kN, (b) t = 200 kN, (c) t = 300 kN, (d) t = 500 kN. 

In the text 
Figure 4. History of convergence. 

In the text 
Figure 5. Optimal structures with W being used as the stiffness criterion. (a) t = 100 kN, (b) t = 200 kN, (c) t = 300 kN, (d) t = 500 kN. 

In the text 
Figure 6. Optimal structures with UC being used as the stiffness criterion. (a) t = 100 kN, (b) t = 200 kN, (c) t = 300 kN, (d) t = 500 kN. 

In the text 
Figure 7. Deformed optimal structures of different stiffness criteria. (a) KS (t = 100 kN), (b) U (t = 100 kN), (c) U _{ C } (t = 100 kN), (d) KS (t = 200 kN), (e) U (t = 200 kN), (f) U _{ C } (t = 200 kN), (g) KS (t = 300 kN), (h) U (t = 300 kN), (i) U _{ C } (t = 300 kN), (j) KS (t = 500 kN), (k) U (t = 500 kN), (l) U _{ C } (t = 500 kN). 

In the text 
Figure 8. Loaddisplacement curves of the optimal structures (load is 500 kN). 

In the text 
Figure 9. (a) Design problem and (b) initial design of the second example 

In the text 
Figure 10. The optimal structure of the second example. 

In the text 
Figure 11. History of convergence of the second example. 

In the text 
Figure 12. Intermediate designs of the optimization of the second example. (a) Step 55, (b) step 90, (c) step 110, (d) step 200. 

In the text 
Figure 13. (a) Design problem and (b) initial design of the third example. 

In the text 
Figure 14. The optimal structure of the third example. 

In the text 
Figure 15. History of convergence of the third example. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.