Static reinforcement and vibration reduction of structures using topology optimization

. This paper presents a topology optimization formulation based on the Solid Isotropic Material with Penalization (SIMP) method and solved by the Modified Optimality Criteria (MOC) algorithm. It addresses mechanical design problems such as structural reinforcement adding elastic material or vibration reduction using viscoelastic layers. The aim is thus to attach on a pre-existing given structure a design domain in order to improve the behavior of this elastic structure, according to an objective function. This can be useful when one wants to use, for example, additive manufacturing to reinforce a pre-existing structure or to maximize structural damping. Two objective functions are tested in both linear statics and dynamics: a compliance based objective function and a displacement based one. In the dynamic case, written in the frequency domain, the two proposed objective functions include the viscoelastic material model (a Zener fractional derivative one) used to fill the design domain. The displacement criteria is developed using a general formula able to take into account as many degree of freedom as necessary. Finally, some applications based on beams and CubeSat-like structures are shown in this article. The proposed examples show that in both statics and dynamics, the optimization of a restrained design domain attached to an existing structure can improve its behavior: stiffness improvement or vibration reduction.


Introduction and motivation
The structural mass is often a main issue which has to be cautiously considered during the design process of a mechanical structure, for instance for a vehicle or an equipment in the aerospace industry. This is a structural optimization problem where one wants to reduce the mass regarding to another criterion such as, non-exhaustively: stiffness maximization usually expressed as a compliance minimization problem, displacement lowering or vibration reduction [1]. Such objectives can be effectively addressed by using a topology optimization approach [2][3][4][5][6][7][8][9][10][11].
The present article uses an algorithm based on the well known SIMP-method [12][13][14] in order to focus on these objectives and tries to address the situation when one wants to add to a pre-existing given structure Ω 1 , a second domain Ω 2 with the aim of either reinforcing the pre-existing structure, lowering the displacement of a chosen area, or reducing vibrations of Ω 1 using the design domain Ω 2 . Although few papers exist on the subject of structural reinforcement using topology optimization [15][16][17][18], the idea of considering separated substructures is newly addressed here. This means that mass (for dynamic * e-mail: antoine.legay@lecnam.net formulation) and stiffness matrices of each substructure are considered separately and assembled afterwards in the global system. During the optimization process, the discretized operators of Ω 1 remain constant (stiffness and mass matrices) while they vary for Ω 2 according to its material density distribution. In order to perform this topology optimization, the design variables are the elemental material density distribution in Ω 2 .
One can therefore imagine to reinforce Ω 1 using an efficient manufacturing process such as additive manufacturing [19] to build Ω 2 .
Concerning the purpose of vibration reduction, it is performed in the frequency domain and damping is taken into account thanks to viscoelastic materials [20][21][22]. Among the various rheological existing models of viscoelastic behavior, the four-parameter Zener fractional-derivative model is used in this work. It allows to efficiently represent the frequency-dependence of the damping properties with only four parameters: correlations between numerical and experimental tests show a good accuracy [23].
The next part of the present paper is dedicated to the building of the discretized numerical model of the physical problem, including a sensitivity analysis of the solution with respect to the design variables, namely the elemental material densities in the design domain. The third section details the topology optimization problem formulation, introducing the two different used objective functions: compliance and displacement based objective functions. The fourth section presents a selection of applications corresponding to a static reinforcement and to a vibration reduction.
The numerical implementation is done using Python and Fortran languages. Gmsh [24] is used for pretreatment (meshing) and post-treatment (plotting).

Modeling of the problem 2.1 Geometry
The studied problem is composed of a purely elastic given domain Ω 1 and a design domain Ω 2 which can be either purely elastic or viscoelastic (Fig. 1). The whole problem is assumed to be linear. The interface between the two domains is denoted by I. In the static case, the aim is to reinforce the given structure Ω 1 subjected to the static load F . In the dynamic case, considered in the frequency domain, the aim is to dissipate energy in order to protect the given structure Ω 1 from spurious vibrations due to the harmonic excitation F . The first case can be degenerated from the second by taking the angular frequency ω to zero.
A design domain Ω 2 is then attached to Ω 1 , in which each element e has a material density variable x e . The aim is to find the best material distribution in Ω 2 according to an objective function associated to constrains functions. This optimization problem is solved using a topology optimization procedure described in Section 3.1.

Viscoelastic model
In this work, the viscoelastic behavior of the dissipative material in Ω 2 is modeled using a four-parameter Zener fractional-derivative viscoelastic model, introduced by Bagley and Torvik [20]. For that purpose, the Young modulus is considered as a complex one such as: where E (ω) is the storage modulus, E (ω) the loss modulus and i 2 = −1. Both moduli can be expressed as a function of the static stiffness E 0 , the dynamic stiffness E ∞ , the relaxation time τ and the fractional order of the derivative α [25], which represent the four parameters of the model: (2) Finally, the Young modulus can be expressed using the four above parameters: This last expression is the one that is implemented to describe the material contained in the sub-part Ω 2 .

Implementation of design variables
For the purpose of topology optimization, the design domain Ω 2 contains a set of design variables x = [x 1 , x 2 , . . . , x n ] T which is the set of elemental material densities (n is the number of elements in Ω 2 -domain's mesh). These densities are such that 0 ≤ x e ≤ 1 where 0 is the lower bound of x e corresponding to an absence of material in the element and 1 is the upper bound meaning a presence of material [26] (Fig. 1). An efficient way to use this is the well-known penalization algorithm SIMP from Bendsøe [12] and Zhou & Rozvany [13]. However, a derivative from this law established by Sigmund [27] and called the Modified SIMP-law is used here. This formulation has several advantages, including the fact that it avoids stiffness (or mass) matrix to become singular.
The effective Young modulus in the element e of Ω 2 is expressed as a function of the elemental material densities x e such that: where E (ω) is given by equation (3), E min = E 0 is the imposed lower bound of the effective Young modulus ( < 1) and p is a penalization factor. In the same way, the effective volumetric mass density in the element e of Ω 2 is expressed as a function of the elemental material densities x e such that: where ρ 0 is the volumetric mass density of the material used for the design domain Ω 2 , ρ min = ρ 0 is the imposed lower bound of the effective volumetric mass density and m is a penalization factor. In this work, the penalization factor p related to the Young modulus is set to 3 and the penalization factor m related to the volumetric mass is set to 1 [5]; the coefficient is chosen to be 10 −3 .

Discretized system
Using the finite element method, the dynamic discretized system of the whole problem is written in the frequency domain as: where K is the complex stiffness matrix (taking into account viscoelastic terms), M is the mass matrix, ω is the angular frequency and F represents the external nodal forces. The design variable vector x, composed of elemental material densities in design sub-domain Ω 2 is introduced in order to use a topology optimization procedure. Using a sub-domain decomposition between domain Ω 1 and Ω 2 , the global discretized system equation (6) is written: where 1, 2 and I denote respectively sub-domain Ω 1 , subdomain Ω 2 and the interface I between Ω 1 and Ω 2 . The stiffness matrix K 2 (ω, x) is complex and frequency dependent due to the viscoelastic material of sub-domain Ω 2 (Eq. (3)). Moreover, this last matrix depends on the design variable vector x.
In order to simplify the notations, the following matrices are introduced: where S 1 (ω) does not depend on x and represents the contribution of Ω 1 , while S 2 (ω) depends on x and represents the contribution of the design domain Ω 2 . The global dynamic equation (Eq. (6)) then becomes By denoting S(ω, x) = S 1 (ω) + S 2 (ω, x), the equation can be written as

Sensitivity analysis of the solution with respect to the design variables
In the following work, a topology optimization procedure is used (Sect. 3.1). This algorithm needs to compute the sensibility of the solution according to the design variables x e . Since only matrix S 2 (ω, x) is depending on design variables, the derivative of equation (11) with respect to the design variable x e leads to the following equation: The sensitivity of the solution with respect to the design variable x e is then given by: where U (x) is the solution of equation (11). In this last equation, the derivative of the S 2 (ω, x) matrix needs to be computed.

Derivatives of the stiffness and mass matrices of the design domain
The derivative of S 2 (ω, x) with respect to x e involves the derivative of the stiffness and mass matrices of Ω 2 , denoted respectively by K 2 (ω, x) and M 2 (x): The stiffness matrix K 2 (ω, x) can be expressed as the assembly of the elemental stiffness matrices k j (ω, x j ) over the n elements of Ω 2 -domain mesh: The elemental stiffness matrix k j (ω, x j ) can be written as where k 0 j is the stiffness matrix of element j computed with a unit Young modulus. Thus, the derivative of K 2 (ω, x) with respect to x e is given by with where δ ej is the Kronecker symbol (no implicit summation). The same development is applied to the mass matrix. The mass matrix M 2 (x) is expressed as the assembly of the elemental mass matrices m j (x j ) over the n elements of Ω 2 -domain mesh: The elemental mass matrix m j (x j ) can be written as where m 0 j is the mass matrix of element j computed with a unit volumetric mass density. Thus, the derivative of M 2 (x) with respect to x e is given by with Since m = 1 in this work, the simplification gives The elemental stiffness and mass matrices obtained with a unit Young modulus and a unit volumetric mass density are computed once at the beginning of the computation and stored in order to be used later on during the topology optimization process.

Minimization problem
In this work, the aim of the topology optimization is to find the material density set of design variablesx ω minimizing the objective function f (ω, x) for a given angular frequency ω according to constraints: where v(x) is the volume of material in Ω 2 associated to the distribution of material density x, v 0 is the total volume of Ω 2 and γ is the final target of material volume ratio in Ω 2 , chosen by the user. Two different objective functions are considered in this work. The first function is based on the compliance (denoted by c(ω, x)) while the second one is based on the displacement (denoted by d(ω, x)). Both can be used under a static environment (ω = 0) as well as a dynamic one (ω > 0). A density filtering strategy [3] is used to avoid checkerboard patterns. This minimization problem is solved using a modified-SIMP algorithm (see Sect. 2.3) associated to a solver called Modified Optimality Criteria (MOC) method [16] which is an extension of the Optimality Criteria (OC) method [14,28]. This algorithm has been chosen thanks to its capability to be adjustable depending on the optimization problem to solve. Also, it is much more suitable than OC due the modified parameter introduced for dealing with dynamic criteria, and relatively easier to implement and analyze compare to other complex algorithms such as MMA [29].
The iterative process is stopped until either a convergence of the objective function is reached, meaning when the difference of the elemental densities between 2 consecutive iterations is less than a fixed relative criterion (noted as ζ in this work); or when the maximum number of iterations is reached (noted as n max in this work).

Dynamic case
The compliance is given in terms of the complex displacement field U (x), and the stiffness and mass matrices of the problem as [30,31] where U T (x) denotes the transpose-conjugate of the U (x) vector and |.| denotes the norm of a complex number. By introducing the following g(ω, x) complex function the compliance can be expressed as a function of g(ω, x): The sensitivity of c(ω, x) with respect to x e is then given by: This last expression involves the sensitivity of g(ω, x) with respect to x e , which can be computed from equation (29): Using equation (13), this equation becomes (33) It can be rewritten as: where X is the solution of the following system: The system of equation (35) has to be solved for each iteration of the topology optimization loop, but not for each design variable x e . Once the sensitivity of the g(x) function with respect to x e is obtained, the compliance sensitivity is given by equation (31).

Static case
The compliance based objective function in the static case (denoted by c s (x)) can be derived from the previous section by taking ω = 0: where In this case, S(0, x) = S(0, x) since S 2 (0, x) is a real matrix. The derivative of g(0, x) with respect to x e becomes: and the static compliance sensitivity is then given by where E (0) = E 0 in the expression of equation (3).

Dynamic case
For some applications, it can be relevant to use a local criterion (e.g. when focusing on a targeted component embedded in a the structure Ω 1 ), such as the minimization of the displacement of a local point (or set of points). The following displacement based function given in terms of the displacement field U (x) [32,33] is then introduced: where β is a diagonal localization matrix which can be defined by (using the Kronecker symbol δ ij with no implicit summation): The coefficient b i is prescribed by the user (0 or 1) depending on whether one wants to consider the degree of freedom (dof) i or not.
Using equation (13), this last equation becomes This expression can be rewritten as: where Y is the solution of the following system: This equation has to be solved for each topology optimization iterative step, while equation (45) gives the sensitivity of the displacement based objective function: where is the real part of the imaginary number.

Static case
The displacement based function in the static case (denoted by d s (x)) can be derived from the previous section by taking ω = 0: where U (x) is the solution of: The sensitivity of d s (x) with respect to x e becomes: where Y is the solution of the following system:

Reinforcement of a 3D cantilever beam
This example is adapted from the cantilever beam from [34] and is composed of an elastic plate (Ω 1 ) surrounded by two layers which compose Ω 2 , as shown in Figure 2.
The beam is fully clamped on its left side while a vertical load of 100 N is applied on the right side of Ω 1 . The structure is discretized using 8-node hexahedral elements: 2160 in Ω 1 and 8640 in Ω 2 , leading to a total of 10,800 elements and 37,719 degrees of freedom. The chosen point for the local displacement minimization is point A (Fig. 2), so matrix β only has 3 unit terms corresponding to the 3 dofs associated to the node A. The optimum solution for both objective functions, static compliance c s (x) and displacement based function d s (x), are compared. The material for both domains Ω 1 and Ω 2 is aluminum whose properties are given in Table 1.
In this case, the target for the final volume ratio is set to γ = 25% of the total volume of Ω 2 , the stop criterion ζ is empirically set to 2.5% and the maximum number n max of iterations is set to 100.
The optimal final shapes obtained by minimizing c s (x) and d s (x) are denoted respectively byx c andx d ; they are given in Figure 3. One can see that the final shapes are quite similar, with only few differences on a local aspect. In both cases, computations have converged towards a physical shape. On a numerical aspect, the evolution of both objective functions (Fig. 4) show that a convergence is reached around 80 iterations for both criteria. In the final shape for solutionx c , there are: 2140 elements such that x e > 0.99 6420 elements such that x e < 0.01 80 elements such that 0.01 < x e < 0.99 In the final shape for solutionx d , there are: 2140 elements such that x e > 0.99 6400 elements such that x e < 0.01 100 elements such that 0.01 < x e < 0.99 The quantity of non-converged elements (material density between 0.01 and 0.99) is about 1% in both cases.

Application to a 3D CubeSat-like structure
The application addressed in this section concerns a CubeSat-like structure which is a standard format created by Puig-Suari and Twiggs [35,36] for the sake of university projects, in order to send small satellites into space. The specifications for the main configuration is called "1U" and is a 10 × 10 × 10 cm cube whose weight can not exceed 1.33 kg. Therefore, it seems appropriate to use topology optimization on an initial CubeSat skeleton in order to optimize the stiffness of the overall structure while controlling the mass of the added material. The considered initial structure is shown in Figure 5 and is composed of a skeleton (Ω 1 ) standing on its 4 clamped feet. A vertical surface load is applied on the upper surface. Design domain Ω 2 is chosen to be the volume resulting from the 6 faces of the cube times the thickness of the skeleton in the third dimension. The structure is discretized using 8-node hexahedral elements: 2808 in Ω 1 and 9504 in Ω 2 , leading to a total number of 12,312 elements and 49,560 degrees of freedom. The material for both domains Ω 1 and Ω 2 is aluminum whose properties are given in Table 1. The two criteria, compliance c s (x) and displacement based objective functions d s (x), are used. The chosen point for the local displacement minimization is one of the upper corner (point A), so matrix β only has 3 unit terms corresponding to the 3 dofs associated to this node A. The target for the final volume ratio is set to γ = 25% of the total volume of Ω 2 , the stop criterion ζ is empirically set to 1% and the maximum number of iterations n max is set to 200.
The optimal final shapes obtained by minimizing c s (x) and d s (x) are denoted respectively byx c andx d ; they are given in Figure 6. One can observe two clearly defined reinforcement shapes depending on the objective function. The evolution of the objective functions are given in Figure 7. It shows that the criteria ζ is reached at around 100 iterations for the displacement objective function while there are spurious oscillations in the material densities for the compliance objective function and the maximum number of iteration n max is reached. Nevertheless it does not affect the final shape of the design domain since it concerns only a few elements. Indeed, in the final shape for solutionx c , there are: 2360 elements such that x e > 0.99 7064 elements such that x e < 0.01 80 elements such that 0.01 < x e < 0.99     In the final shape for solutionx d , there are: 2368 elements such that x e > 0.99 7072 elements such that x e < 0.01 64 elements such that 0.01 < x e < 0.99 The quantity of non-converged elements (material density between 0.01 and 0.99) is less than 1% in both cases.

Optimization of a visoelastic layer for a 2D sandwich beam
The problem is described in Figure 8. It is assumed a plane stress state. The sandwich beam is composed of 2 external aluminum layers (domain Ω 1 ) and a viscoelastic layer (domain Ω 2 ) in between them. The material properties are given in Table 1 (for aluminum) and 2 (for the viscoelastic layer). The left side of the beam is clamped while an harmonic force (±10 N) is imposed at the right side at the angular velocity ω. The structure is discretized using quadrangular 4-node elements. There are 6 elements in the thickness direction of the aluminum layers and 7 elements in the thickness direction of the viscoelastic layer. Along the length direction of the beam, there are 198 elements. There is a total of 1386 elements in Ω 2 . The target final volume ratio is γ = 50% of the total volume of Ω 2 , the stop criterion ζ is empirically set to 1% and the maximum number of iterations n max is set to 100. Firstly, the viscoelastic layer is optimized using the compliance objective function c(ω, x) at three different harmonic excitations ω (500 Hz, 1000 Hz and 2000 Hz). The solutions in terms of material densities are denoted respectively byx c 500 Hz ,x c 1000 Hz andx c 2000 Hz . These optimal material densities are plotted in Figure 9. The convergence of the compliance over the iterations is plotted in Figure 10. It can be seen that for the three harmonic excitations, there are oscillations of the objective function and the optimization process is stopped at a fixed number of iteration (here n max =100 iterations). However, these spurious oscillations concern a few elements and do not affect the final solution (Tab. 3). The frequency response functions of the beam are then computed using these three material densities. The compliance functions c(ω,x c 500 Hz ), c(ω,x c 1000 Hz ) as well as c(ω,x c 2000 Hz ) are shown in Figure 11. The black curve is obtained using a fully filled layer of viscoelastic material (x e = 1, ∀e). It can be seen that for a given angular frequency, a     Secondly, the viscoelastic layer is optimized using the displacement objective function d(ω, x) at the same three different harmonic excitations ω (500 Hz, 1000 Hz and 2000 Hz). The solutions in terms of material densities are denoted respectively byx d 500 Hz ,x d 1000 Hz andx d 2000 Hz . These optimal material densities are plotted in Figure 12. The convergence of the displacement objective function over the iterations is plotted in Figure 13. As for the compliance criteria, there are spurious oscillations for two harmonic excitations (1000 Hz and 2000 Hz) ( Table 3). The frequency response functions of the beam are then computed using these three material densities. The displacement functions d(ω,x d 500 Hz ), d(ω,x d 1000 Hz ) as well as d(ω,x d 2000 Hz ) are shown in Figure 14. The black curve is obtained using a fully filled layer of viscoelastic material (x e = 1, ∀e). As for the compliance based criteria,      the same conclusions can be made: the density optimized solution obtained for a given frequency gives the minimum of the displacement based function at this frequency. One can notice that in this case, for 500 Hz, it is not completely true; the topology optimization process can not give a better solution than the one given by a fully filled layer of viscoelastic material. Nevertheless, the solution given by the optimization is lighter.

Optimization of a viscoelastic layer for a 2D
CubeSat-like structure The problem is described in Figure 15. It is assumed a plane stress state. The structure is composed of external aluminum layers (domain Ω 1 ) and internal viscoelastic layers (domain Ω 2 ). Its shape is the same as the 3D CubeSat-like structure studied in the previous section. The material properties are given in Table 1 (for aluminum) and 2 (for the viscoelastic layers). The target final volume ratio is γ = 50% of the total volume of Ω 2 , the stop criterion ζ is empirically set to 1% and the maximum number of iterations n max is set to 100. The bottom side of the structure is clamped while an harmonic force F is imposed on the upper left corner at the angular velocity ω ( F = 1 N). The structure is discretized with a total of 3222 quadrangular 4-node elements, the number of elements in Ω 2 is 966. The mesh can be seen in Figure 16.
Firstly, the viscoelastic layer is optimized using the compliance objective function c(ω, x) at three different harmonic excitations ω (50 Hz, 100 Hz and 200 Hz). The solutions in terms of material densities are denoted respectively byx c 50 Hz ,x c 100 Hz andx c 200 Hz . These optimal material densities are plotted in Figure 17. The convergence of the compliance over the iterations is plotted in Figure 18. As for the 2D sandwich beam, it can be seen that for the three harmonic excitations, there are oscillations of the objective function and the optimization process is stopped at a fixed number of iteration (here 100 iterations). However, these spurious oscillations concern a few elements and do not affect the final solution (Table 4)      are then computed using these three material densities. The compliance functions c(ω,x c 50 Hz ), c(ω,x c 100 Hz ) as well as c(ω,x c 200 Hz ) are shown in Figure 19. The black curve is obtained using a fully filled layer of viscoelastic material (x e = 1, ∀e). As for the 2D sandwich beam, it can be seen that for a given angular frequency, a solution can be better than another one in terms of minimizing the compliance. For instance, the material density solutionx c 50 Hz given by the minimization of the compliance at 50 Hz exhibits the minimum of compliance at 50 Hz compared to the other solutions, even the solution computed with a full layer of viscoelastic material. The same conclusions can be done with solutionx c 100 Hz at the frequency 100 Hz; and with solutionx c 200 Hz at the frequency 200 Hz.
Secondly, the viscoelastic layer is optimized using the displacement objective function d(ω, x) at the same three different harmonic excitations ω (50 Hz, 100 Hz and 200 Hz). The solutions in terms of material densities are denoted respectively byx d 50 Hz ,x d 100 Hz andx d 200 Hz . These optimal material densities are plotted in Figure 20. The convergence of the displacement criteria over the iterations is plotted in Figure 21. As previously, two of the three optimizations have spurious oscillations in terms of objective function over the iteration process (Table 4). These oscillations are more important for x d 50 Hz : for this case, there are no elements with x e < 0.01 but there are 362 elements with a material density less than 0.1; and so only 154 elements between 0.1 and 0.99. In order to illustrate the oscillations in   Figure 23. The black curve is obtained using a fully filled layer of viscoelastic material (x e = 1, ∀e). As for the compliance based criteria, the same conclusions can be made: the density optimized solution obtained for a given frequency gives the minimum of the displacement based function at this frequency.
However, sometimes, for this application, the topology optimization process does not improve so much the solution compared to the one given by a fully filled layer of viscoelastic material. Nevertheless, the solution given by the optimization is lighter.

Conclusion
The aim of this work is to use a topology optimization process to design a domain Ω 2 interacting through an interface with a fixed given structure Ω 1 in order to improve the behavior of the initial structure according to a given criteria. The problem can be static or dynamic. In the dynamic case, the problem is treated in the frequency domain and the used material in the design domain has viscoelastic properties able to reduce the vibration amplitude. This viscoelastic material is modeled using a four-parameter Zener fractional-derivative viscoelastic model.
Two objective functions are introduced in this paper: a global compliance based one and a local displacement based one. They are both defined for the dynamic as well as for the static case. Results show that the optimum solutions obtained for the two objective function are very often quite similar.
The formulation of the problems allows to have a complex shape of the additional design, even a design domain separated into several volumes, as it is the case in a few presented applications. Moreover, the chosen material for the design domain can be different than the one used for the initial structure. The shape and the material of the design domain have simply to be coherent with functional surfaces of the system and with the manufacturing process.
The applications in the static case show that it is possible to improve a given structure by adding material around it using a topology optimization process. The strategy may be applied to a damaged structure needed to be repaired. Additive manufacturing may be used to manufacture the domain Ω 2 in order to reinforce Ω 1 .
The examples in the dynamic case show that the proposed strategy enables to design additive viscoelastic layers to damp the vibrations at a given frequency.
The proposed examples deal with sandwich structures for which the viscoelastic layer in between two layers of metal since it is known to be more efficient than on the external surface of the structure. The final shape of the viscoelastic layer drastically depends on the target frequency.
Nevertheless, especially in the context of dynamics and viscoelasticity, the topology optimization process often reaches the maximum number of iterations, exhibiting spurious oscillations of the material densities between few elements. However, these oscillations concern a few elements, and do not affect the quality of the final solution.
Conclusions show that it is possible to deal with two different substructures with the aim of controlling the first one using the second one. The way it is implemented here allows for a great variety of configurations for future users.

Funding information
ArianeGroup is gratefully acknowledged for its financial support.