© Q. Xia and T. Shi, Published by EDP Sciences, 2016

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 quasi-static 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 load-displacement 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 load-displacement curve of geometrically nonlinear structure.

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

Compliance or the so-called 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 multi-objective 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 multi-objective optimization. Therefore, in the present study, we propose to extend the inspiring previous studies by considering load-normalized strain energy increments between consecutive load steps. The increments are aggregated through the Kreisselmeier-Steinhauser (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 load-normalized 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 ⁰Ω and ⁰Γ_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 Piola-Kirchhoff stress tensor; is the virtual Green-Lagrange 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.Venant-Kirchhoff 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 Green-Lagrange strain ε_kl(u) is given by(6)the virtual Green-Lagrange strain is given by(7)

Because is nonlinear in the argument u, the Newton-Raphson method is used to solve the nonlinear equilibrium equation, i.e., equation (2), in an iterative manner. The displacement ^m+1u of the (m + 1)-th iteration is computed as the sum of ^mu 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 quasi-static condition).

The true work of external loads under quasi-static 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 s-th 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 Newton-Raphson method for the s-th 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 multi-objective optimization.

In the present study, load-normalized strain energy increments between consecutive load steps are aggregated through the Kreisselmeier-Steinhauser (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 load-normalized 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 y-axis is applied to a geometrically nonlinear structure. Then, in equation (22), the components of fictitious load are set to l₁ = 0 and l₂ = −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 load-normalized 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 perimeter-based 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 Newton-Raphson algorithm arrives at a converged solution u.

Ignore the shape derivatives defined on boundary ⁰Γ _N since ⁰Γ _N is not subject to optimization, and note that only ⁰Γ _H is optimized, we have the partial derivative of with respect to ⁰Ω 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 s-th 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 ⁰Γ _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 ⁰Γ_H of a structure.

Propagation of ⁰Γ_H is modeled by the Hamilton-Jacobi 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; 3-node 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⁻⁴ because the order of magnitude of KS function is 10⁻¹. 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⁴–10⁵. 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].

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 multi-objective 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 load-normalized 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 load-displacement 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¹ of the red curve is smaller than that of the blue curve.

Figure 8. Load-displacement 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¹ 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¹ 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¹ should not dominate the optimization (recall that Δu¹ 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¹ = 0.1345, Δu² = 0.1033, Δu³ = 0.0907, Δu⁴ = 0.0800, Δu⁵ = 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⁻⁴.

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¹ = 0.0108, Δu² = 0.0109, Δu³ = 0.0111, Δu⁴ = 0.0115, Δu⁵ = 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 well-known 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⁻⁴.

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¹ = 0.0157, Δu² = 0.0160, Δu³ = 0.0164, Δu⁴ = 0.0171, Δu⁵ = 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 load-normalized 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

All Tables

All Figures

	Figure 1. A load-displacement 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. Load-displacement 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