Jie HOU, Xiojun GU, Jihong ZHU,b,c, Jie WANG, Weihong ZHANG

a State IJR Center of Aerospace Design and Additive Manufacturing, Northwestern Polytechnical University, 710072 Xi’an, China

b Institute of Intelligence Material and Structure, Unmanned System Technologies, Northwestern Polytechnical University,710072 Xi’an, China

c MIIT Lab of Metal Additive Manufacturing and Innovative Design, Northwestern Polytechnical University, 710072 Xi’an, China

KEYWORDS

Abstract This paper presents an extended topology optimization approach considering joint load constraints with geo-metrical nonlinearity in design of assembled structures. The geometrical nonlinearity is firstly included to reflect the structural response and the joint load distribution under large deformation. To avoid a failure of fastener joints, topology optimization is then carried out to minimize the structural end compliance in the equilibrium state while controlling joint loads intensities over fasteners. During nonlinear analysis and optimization, a novel implementation of adjoint vector sensitivity analysis along with super element condensation is introduced to address numerical instability issues. The degrees of freedom of weak regions are condensed so that their influences are excluded from the iterative Newton-Raphson (NR) solution. Numerical examples are presented to validate the efficiency and robustness of the proposed method.The effects of joint load constraints and geometrical nonlinearity are highlighted by comparing numerical optimization results.

1. Introduction

Topology optimization is an efficient way to achieve an idealized conceptual design satisfying various performance requirements. Since the pioneering work in 1980s by Bendsoe and Kikuchi1,it has been an everlasting hotspot in both academia and industry2–4. Till now, the application scope of topology optimization is widely expanded to aerospace,automobile, civil engineering, etc. However, the majority focuses on linear elastic problems, which assume that structures undergo infinitesimal strain. When structural deformation exceeds a certain range, the small deformation assumption will be invalid.

In engineering practice, a nonlinear effect is urgently integrated with design methodologies. Modern aircraft, especially high-altitude long-endurance (HALE) aircraft, usually use a high aspect ratio and flexible wings. Therefore, elastic deformation will increase and sometimes exceed 25% of the wingspan5. Meanwhile, the stress does not exceed the plastic yielding limit, and the constitutive relation remains linear. In such cases, geometrical nonlinearity should be considered in aircraft structure design6.

With regard to topology optimization considering nonlinearity, Jog7focused on the geometrical nonlinearity in topology optimization, and used an extended perimeter method to solve nonlinear thermoelasticity. Sigmund8applied topology optimization to the design of a compliant mechanism undergoing a large displacement. Later, he developed a topology optimization method for metaphysics actuators and electrothermo-mechanical systems considering nonlinear analysis9.Bruns and Tortoreli10studied the topology optimization of a snap-through structure.Jung and Gea11derived the sensitivity of a general displacement function, and carried out optimization with geometrical and material nonlinearity simultaneously. Hung and Xie12solved problems with geometrical and material nonlinearity in topology optimization using the Bidirectional Evolutionary Structural Optimization (BESO)method alternatively. Cho and Kwak13used a mesh-free method, specifically, reproducing the kernel method to optimize a structure with geometrical nonlinearity.He et al.14proposed an analysis-independent density variable approach for topology optimization of geometrically nonlinear structures based on the element-free Galerkin (EFG) method. Gea and Luo15used a microstructure-based design domain and sequential convex approximation method to solve the stiffness optimization problem with geometrical nonlinearity.

The main obstacles for topology optimization with nonlinearity are the non-convergence issue and numerical instability of the nonlinear finite element analysis brought by a densitybased optimization method. To solve these problems, Buhl et al.16modified the convergence criterion of nonlinear analysis to circumvent the influence of a weak element. Convergence problems were seldom encountered. Bruns and Tortorelli10suggested using a hyper-elastic constitutive model to make void elements stiffer for large deformations and thereby to prevent the elements from unreasonable deformation. They also proposed an element removal and reintroducing strategy to solve the computational problem caused by a low-density element17.Similarly,Kang and Luo18fulfilled the regeneration of elements in conjunction with the sensitivity filter technique.Luo and Tong19extended the moving iso-surface threshold method to exclude a void element in the nonlinear finite element analysis, and the reappearance of elements was implemented via interpolating, filtering, and intersecting of the constructed response function. Wang et al.20proposed an energy interpolation model. Based on an analysis in the solid region on the nonlinear stored energy and an analysis in the void regions on the linear stored energy, the intermediate region was interpolated via a density-related factor with a reduced geometrical effect.

Here, we consider a largely deformed assembled structure where different components are jointed as an integrity by bolts, rivets, and the like. In addition to the usual strength design in topology optimization problems with stress and fatigue constraints21,22, multi-fastener joint loads are considered.Existing literature studied the problem in the category of linear elasticity. For example, Ekh and Schon23evaluated the effects of different parameters on the load distribution in a multifastener joint. The joint layout was optimized to minimize bearing stresses. Chickermance et al.24performed a joint location optimization to control their bearing loads without a decrease of stiffness. Zhu et al.25and Hou et al.26proposed a topology optimization to control joint loads intensities by designing connected structures.Zhang et al.27studied the concentrated force diffusion problem in the connection section of a launch vehicle using topology optimization. The studies above have shown the strong influence of joint loads on structural performance design,but the limitation of linear elasticity will prevent the application of ideas in largely deformed engineering structures.

In this paper, geometrical nonlinearity is therefore taken into consideration to have an accurate structural deformation.To favor the application of commercial finite element codes,which are convenient in solving largely deformed assembled structures with various element types but are limited in achieving analytical design sensitivities and modified convergence criterion,enhanced strategies involving the super element method and the continuation method are used to circumvent the inherent numerical instability caused by a density-based method.A novel implementation of adjoint vector sensitivity analysis is used to derive the sensitivities of end compliance and joint loads. Based on these treatments, the standard topology optimization model is here extended with considerations of geometrical nonlinearity and joint load constraint to improve the load path with a control of the joint load magnitude.

2. Basic theory and consideration of optimization

2.1. Nonlinear finite element analysis

Geometrical nonlinearity generally exists in aircraft and aerospace thin-walled, slender structural members. In these cases,structures are assumed to undergo large displacements, but mechanical strains are relatively small. For large displacement gradients, the second-order terms in the strain measure are no longer negligible, but the constitutive relationship of linear elasticity remains valid. The Green-Lagrange strain measure ηijis used and defined as

whereuis the nodal displacement.

Correspondingly,the conjugated stress measure,the second Piola-Kirchhoff stresssij, is defined as

The constitutive law can thus be expressed as

whereCis the elastic matrix. The residual force, which is the error between external and internal forces, is defined as

The external loadfextis considered as design-independent,while the internal nodal forcefinis defined as

whereBis the geometrical matrix that transfers the displacement into strain, andVis the structural volume. Adopting the total Lagrange formulation, all the quantities are referred to the original undeformed configuration, andBreads

whereBLis identical to the linear geometrical matrix, whileBNLis the function of the displacement.Essentially the nonlinear effect is brought in viaBNL. The equilibrium state is achieved when the residual force is lower than a given tolerance. The Newton-Raphson method and the incremental method are used to solve the nonlinear equation iteratively,which requires the tangent stiffness matrix, i.e.,

whereUis the displacement vector for the structure in its equilibrium state. The tangent stiffness matrix in the equilibrium state is important to the derivation of sensitivity. Detailed derivation of the tangent stiffness matrix can be found in many finite element analysis textbooks28.

2.2. Optimization formulation

Once the nonlinearity is taken into consideration, the equilibrium state is achieved iteratively.The formulation of optimization should be rephrased considering the solving process.Topology optimization considering the joint load constraint and geometrical nonlinearity can be expressed as

where ρ is the vector of pseudo-density design variables describing the material distribution in the design domain.g（U） is the objective function. In early research, end compliance (compliance in the equilibrium configuration), complementary work, and weighted end compliance are chosen as the objective function16.

Additionally, material interpolation maps the design variable ρ to the physical properties of the structure. A popular material interpolation named SIMP (Solid Isotropic Material with Penalization) can be written as wherepis the penalty factor and typically assigned larger than 1.E0is the Young’s modulus of the selected material. Meanwhile, the artificial intermediate state generates an oscillation of the structural displacement or even the singularity of the tangent stiffness matrix, which will lead to non-convergence of nonlinear finite element analysis.This issue will be discussed later in this paper.

2.3. Sensitivity analysis with adjoint vector sensitivity analysis

Sensitivities are fundamental to gradient-based optimization algorithms. Theoretically, the sensitivities of nonlinear topology optimization problems can be analytically expressed as presented in existing literature such as the works of Buhl16,Kang and Luo18,and Wang et al.20,where the tangent stiffness matrix was directly obtained from self-programmed finite element codes. However, when we use commercial codes such as ANSYS,ABAQUS,etc.,obtaining the tangent stiffness matrix is unfortunately limited. Therefore, we propose a novel implementation of adjoint vector sensitivity analysis. Straightforward, the end compliance designated as the objective function can be expressed as

Using the adjoint method,we introduce a Lagrangian multiplier ψ. Assuming that the equilibrium has been found, the term ψTRis equal to zero and can be added to the objective function without any changes.The modified objective function reads

Here, we considerfextas design-independent, i.e., the derivatives offextwith respect to any design variables equal to zero.Based on the chain rule,the sensitivity of the modified function with respect to the pseudo-densities reads

Noticing thatR=0, the Lagrange multiplier vector ψ can be chosen freely. To eliminate the unknown term ?R/?U， ψ can be chosen so that

which corresponds to solving the linear equations as follows:

In the equilibrium state under the external loadfext,the tangent stiffness matrix is used as the stiffness matrix of the adjoint structure. ψ can be regarded as the linear elastic displacement of the adjoint structure. In most commercial finite element software, the tangent stiffness matrix is embedded in the program kernel,and there is no access to invoke an adjoint analysis. Notice that the equilibrium is solved using the Newton-Raphson method incrementally. Assuming that the load is divided into$N$uniform steps,equilibrium is achieved for every load increment.The load increment ΔPcan be calculated as

Take two equilibriums, namelyUNfor the final load step andUN-1for the step before the final. When the step number is large enough and near convergence, the tangent stiffness matrices of two equilibriums can be considered as approximately equal, i.e.,

In this fashion,the adjoint structure can be rewritten by the following linear equation:

Combining Eqs. (14)–(17), ψ can be given by

Therefore, the sensitivity of the end compliance can be expressed as

Thus,the only unknown term is the sensitivity of the residual with respect to design variables, which can be obtained by differentiation of Eq. (5). By assuming that the geometric matrix is independent of design variables, one obtains the derivative of the residual with design variables as

wherefiis the element nodal force.

The structural response of the adjoint structure is obtained by one finite difference using the equilibrium states of the last load increments. The accuracy of the approximated structural response is dependent on the perturbation, i.e., the magnitude of the load step. A proper choice of the load step is problemdependent and shall be decided by trial and error.

2.4. Formulation of the super element method

To eliminate the inherent numerical instability of a densitybased method in nonlinear problems, many strategies are proposed, among which removing void elements (low-density elements) is the most straightforward approach. However, an absence of the void elements leads to potential discontinuity of the structural mesh and singularity of the finite element analysis.In addition,topology optimization requires an evaluation of the sensitivities in the void domain for a reappearance of the void elements.

To solve this issue, the super elements method (SEM) is used.The idea of the SEM is originally brought up to promote computing efficiency of large-scale structures by condensing a group of elements into a single macro-element29.It is used here to replace void elements using macro-elements,which are confined to linear elastic deformation, and the numerical instability caused by the void elements in the void area is improved.

To integrate the SEM, the design domain is divided into a void part and a solid part by an element density threshold ξ.Naturally, the nodal degrees of freedoms are categorized into 3 types,the DoFs of the solid elementsus,the DoFs of the void elementsuv,and those on the interface shared by the two partsub.

In the void region, the structural displacement responses are derived by the following linear equilibrium equation:

where the stiffness matrix is linear elastic. Thus,uvcan be obtained as the following:

The nodal force of the void part on the common boundary can be expressed as

where

In the solid region, the equilibrium is calculated by the Newton-Raphson method. In the （t+1）th iteration of theNth load increment, the relationship between the force and the displacement yields

Adding the shared DoFs in the void region, the structural nonlinear equilibrium equation reads

Substituting Eq. (29) into Eq. (19), the sensitivity can be obtained.

By constructing super elements in the void region,the void elements are excluded from nonlinear analysis, but computed in the linear elastic regime. The computation procedures of the void elements given in Eqs.(22)–(25)are conducted following the classic routine of the SEM.In the solid region,equilibrium equations in Eqs. (26) and (27) are solved by finite element software. Thus, extra operations of the tangent stiffness matrix are avoided. The sensitivities are derived in the two regions, respectively. Hence, an introduction of the SEM not only holds a higher computational efficiency,but also prevents an appearance of element distortion. It is notable that the condensation will inevitably introduce numerical errors.Hence,the values of ξ should be selected carefully.If the value of ξ is too high, the difference between the condensed and uncondensed elements is big. In this situation, the condensation will increase the stiffness of an element compared with that of the nonlinear formulation.

3. Topology optimization-based joint load control

3.1. Measure of the joint load with p-norm

The implementation of joint load control in topology optimization requires a precise measurement of the joint load.To avoid extra complexities and computational difficulties in simulating the behaviors of multi-fastener joints in topology optimization, many simplified modeling techniques and optimization schemes have been proposed by Chickermance et al.24,Ekh and Schon23,and Zhu et al.25.In this paper,beam elements are used to model multi-fastener joints to obtain the joint loads in terms of design variables.

For a beam element with a shear effect sketched in Fig.1,Fzis the elemental axial force,whileFxandFyare the components of shear forces which are critical to structural safety and mostly concerned in joint load optimization. According to the formulation of a beam element with a shear effect, the joint load of the element can be expressed in terms of beam parameters and nodal displacement as

Fig.1 Free body diagram of a beam element.

whereEandGare the elastic and shear moduli of the beam,respectively. The cross section is assumed to be circular, andI,S, andLdenote the moment of inertia, the cross-section area, and the beam length, respectively. φ is the shear coefficient, and ω is the shear factor of the cross section and equals to 10/9 for a circular beam.uAxand θAxare the corresponding nodal displacement and rotation angle,respectively,of node A in thexdirection.Similar definitions are used foruAy,θAy,uBx,θBx,uBy, and θBy. For brevity of expression, extract the coefficient vectors α, β, and γ of a beam element from Eq. (30) so that

The shear loadFsand the tensile loadFtcan be expressed as

The joint load optimization problem will thus become a min-max problem. Naturally, the joint load is evaluated and constrained in each element.In this way,the overall stress level is restricted below the prescribed limit, but the computational cost increases dramatically. On the contrary, the approach of only considering the maximum load will fail when the location of the maximum joint load changes. To overcome the issue,aggregation methods combining all local evaluations into one or several Kreisselmeier-Steinhauser (KS) or p-norm global functions are generally used30,31.

Here,the standard p-norm measure is utilized,which reads

whereFPNis the p-norm of the measured load, andFican be either the shear load or the tensile load according to the design requirements. pn is the aggregation parameter. Hence, the optimization constraint can be written as

Obviously, with a sufficiently large value of pn, the global measure can exactly match the maximum joint load.However,a large value of pn may lead to oscillation and ill-posed problems32.Therefore,to balance the aggregation quality and convergence efficiency,the value of pn should be chosen properly.

3.2. Sensitivity analysis of the joint load with p-norm

Generally, joint load constraints can be expressed in terms of density design variables. With the aggregation of p-norm, the sensitivity can be expressed as

where the sensitivity can be regarded as the summation of the derivatives of the joint loads with coefficients. The derivatives of the joint loads with respect to the design variables can be derived from the differentiation of Eq. (32) and written as

Combining Eqs. (32), (36), and (37), the sensitivities of p-norm read

whereLsandLtare defined as

Following the derivation of the end compliance sensitivity,terms ?u/?ρican be regarded as the adjoint vector.The sensitivities of p-norm can be converted to the sensitivity ofu*,where

Here, the adjoint method is used by introducing the Lagrangian multiplier λ as follows:

Again,for the equilibrium state,Requals to zero.Then the sensitivity of p-norm can be rewritten as

where ?R/?Uequals to the minus tangent stiffness matrixKT,and the derivation of ?R/?ρihas been discussed previously.To eliminate the unknown term dU/dρi, the Lagrangian multiplier λ is chosen so that

which corresponds to solving the linear equations as

After the computation of λ, the sensitivity of the p-norm can be expressed as

4. Numerical examples

In this section, the optimization algorithm GCMMA (Globally Convergent Method of Moving Asymptotes33) implemented in the general-purpose design platform Boss-Quattro is used as an optimizer.To eliminate the checkerboard pattern,a sensitivity filter is used with a filter radius of three times of the element size.

4.1. SEM verification

To verify the super element method, a C-shaped problem is presented. The problem was firstly introduced by Yoon and Kim34to study the non-convergence issue of a nonlinear finite element analysis and then further tested with topology optimization. The dimension of the initial configuration is sketched in Fig.2 and meshed with 100×100 4-node plane stress elements. The thickness is 1 m. The Young’s modulus is 3 GPa, and Poisson’s ratio is 0.4. Design variables are related to the Young’s modulus by an SIMP model. The square area is artificially divided into two parts: the C-shape area is occupied with solid material whose design variables are initially assigned as 1, and the design variables of the rest area are assigned as a lower bound of 0.001.As a consequence,the Young’s modulus is 3×10-6GPa in the void region.

Referring to the result of the linear elastic finite element analysis shown in Fig.3, the maximum deformation in the structure is 772.693 mm, and the total elastic compliance is 314.176 J. In the deformed configuration, highly inverted elements are observed in the void region.Those inverted elements are formulated on the undeformed configuration, and their stiffness matrix is positive definite. It is notable that a nonlinear finite element analysis will not converge using a regular method because of the highly inverted elements. To overcome this issue, the SEM is used to replace the void elements in the nonlinear finite element analysis, and the simulated result is shown in Fig.4. For comparison, a nonlinear analysis with completely removed void elements is carried out, shown in Fig.4(b). The simulated results of the solid region are almost identical. The end compliance of the void elements is 0.017%of the total energy, which means that it has little influence on the simulated results. The removal of void elements leads to a converged analysis result, but possibly introduces discontinuity in topology optimization. Using the SEM can address the numerical instability while retaining the presence of void elements in finite element analysis.

Fig.2 Initial domain of the C-shaped problem and load condition.

Fig.3 Displacement contour of the C-shaped problem using a linear finite element analysis.

For further verification of the SEM,topology optimization is carried out on the C-shape problem to minimize the end compliance in the equilibrium configuration with a volume fraction of 28%. Two initial configurations are considered here.One is sketched in Fig.2,referred as initial configuration 1.The other holds homogeneous design variables of 1,referred as initial configuration 2.The threshold of condensed elements ξ is originally assigned as 0.25 with a decrease of 0.01 every 10th iteration.The target value of ξ is 0.1.Figs.5and 6 present the iteration process of the two initial configurations. Convergence of nonlinear finite element analysis can be achieved in both cases. The optimization problem of initial configuration 1 converges after 103 iterations, while the other only needs 83 iterations to converge. At the beginning of optimization,elements except those in the C-shaped area of configuration 1 are condensed, and none is condensed in the second configuration. The final designs of the two initial configurations are almost identical,which means that the deviation introduced by the SEM does not affect the optimization results.The iteration curves of the end compliance and the condensed element number are plotted in Fig.7,and the iteration curve of the volume is plotted in Fig.8. For the case of configuration 1, the void region is condensed at the beginning of iteration. The number of condensed elements changes continuously. Meanwhile,there initially exists no condensed region in configuration 2.In the 14th iteration,low-density elements appear which cause a sudden change of the condensed element number. The end compliance increases accordingly. Finally, all the curves converge stably.The end compliances of the two optimized results are 8.05 J and 7.99 J, respectively. The numbers of condensed elements are 6599 and 6608, respectively, which are also quite close. It can be seen that using the SEM in an optimization problem considering geometrical nonlinear effects can solve the nonlinear non-convergence problem caused by void elements and ensure the robustness of the optimization process.

Fig.4 Comparison of the displacement contours using nonlinear finite element analysis between with and without void elements.

Fig.5 Design iterations with the C-shaped initial configuration.

Fig.6 Design iterations with the homogeneous initial configuration.

Fig.7 Iteration histories of the end compliance and the number of condensed elements.

Fig.8 Iteration history of the volume.

4.2. I-beam with connecting sheets

A slender I-beam fastened to two thin sheets sketched in Fig.9 is optimized with the joint load constraints here. The I-beam consists of two 10 mm-thick flanges and a 5 mm-thick web.The two sheets are 10-mm thick. These components are connected by 4 rows of parallel-arranged fasteners, namely rows of A, B, C, and D. The structure is globally meshed with 10 mm×10 mm quadrangular shell elements. Fasteners are modeled by circular-section beam elements with a radius of 5 mm. The I-beam is left-side clamped. A point load of 10 kN is applied. The Young’s modulus is 210 GPa, and Poisson’s ratio is 0.3.

Fig.9 An I-beam structure with a distribution of fasteners.

Fig.10 Optimized topologies of the I-beam web.

The optimization problem is set up to minimize the end compliance with 50% of the material in the initial design domain.The threshold of ξ is set as 0.1.Additionally,the joint loads on 92 fasteners are constrained. Here, p-norm is used to aggregate the 92 constraints with pn assigned as 6. The upper bound of the p-norm function is set as 3600. Besides, the geometrically nonlinear effect is not negligible due to the high slender ratio (the ratio of the length to the width is 7.5). For comparison, four different optimizations are carried out. The first case is simulated using a linear elastic model and optimized without joint load constraints, as shown in Fig.10(a).The second case is simulated using a linear elastic model and optimized with joint load constraints, as shown in Fig.10(b).The third case is simulated using a geometric nonlinear model and optimized without joint load constraints, as shown in Fig.10(c).The last case is simulated using a geometric nonlinear model and optimized with joint load constraints,as shown in Fig.10(d).

The optimized configurations in Fig.10(a)and(b)are symmetric along the neutral line of the structure. These two cases are simulated using a linear elastic model.The maximum joint load arises near the clamped end in both results.With the joint load constraints in Fig.10(b), the maximum joint load decreases from 3741.4 N to 3589.9 N. The end compliances of the optimized structures in Fig.10(a) and (b) are 19.41 J and 20.67 J, respectively. The loss of stiffness in case (b) can be considered as a tradeoff for load distribution amelioration.When geometrical nonlinearity is taken into consideration,the optimized results are significantly different,as shown in Fig.10(c) and (d).

Fig.11 Principal stress of the I-beam.

The asymmetry along the neutral line of the structure disappears due to the nonlinear effect.It can be explained that nonlinear finite element analysis takes the difference between the tensile and compressive stress fields into account. The joint loads of the upper and lower rows in the same position are different in both configurations. In case (c), the maximum joint load 4026.6 N arises in row B, and it is 3583.7 N in row A in case (d). The end compliances of the optimized structures in Fig.11(c) and (d) are 17.51 J and 18.05 J, respectively.

For an in-depth understanding of the mechanism of joint load constraints, the principal stresses of the optimized structure are presented in Fig.11.The most severe joint load always appears near the clamped end.Hence,more material is moved towards the clamped end with the joint load constraints. This phenomenon holds for both linear and nonlinear finite element analysis.For linear analysis,the structure in Fig.11(b)retains its symmetry, while a more inclined push is provided by the structural branch to offset the shear effect in the optimized topology. With regard to the optimized results obtained by nonlinear analysis, the maximum joint load shifts from row B to row A.The branches near the clamped end provide extra support to reduce the load intensity.

5. Conclusions

In this paper,a novel semi-analytic sensitivity derivation along with the super element method is presented to solve the geometrical nonlinear problem in topology optimization. Topology optimization is extended with geometrical nonlinearity to reflect the structural response in a more realistic way, and sensitivity analysis is provided in detail. The DoFs in the void region are condensed to eliminate the numerical instability using the super element method. The extended topology optimization with joint load constraints empowers the design ability with an accurate load path design. Validated numerical results show the efficiency and robustness to minimize the nonlinear end compliance with joint load constraints. The semianalytic method not only is restricted to geometrical nonlinearity, but also can further solve other nonlinearity caused by material and contact.

Acknowledgements

This study was co-supported by National Key Research and Development Program (No. 2017YFB1102800), NSFC for Excellent Young Scholars (No. 11722219), and Key Project of NSFC (Nos. 51790171, 5171101743, 51735005,11620101002, and 11432011).