Issue 
Mechanics & Industry
Volume 24, 2023



Article Number  14  
Number of page(s)  18  
DOI  https://doi.org/10.1051/meca/2023003  
Published online  03 May 2023 
Research Article
Static reinforcement and vibration reduction of structures using topology optimization
Laboratoire de Mécanique des Structures et des Systèmes Couplés (LMSSC), Structural Mechanics and Coupled System Laboratory, Conservatoire National des Arts et Métiers (Cnam), 292 rue SaintMartin, 75003 Paris, France
^{*} email: antoine.legay@lecnam.net
Received:
23
August
2022
Accepted:
13
January
2023
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 preexisting 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 preexisting 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 CubeSatlike 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.
Key words: Static Reinforcement / Vibration Reduction / Topology optimization / Viscoelasticity
© S. Burri and A. Legay, Published by EDP Sciences, 2023
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1 Introduction 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, nonexhaustively: 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–11].
The present article uses an algorithm based on the well known SIMPmethod [12–14] in order to focus on these objectives and tries to address the situation when one wants to add to a preexisting given structure Ω_{1}, a second domain Ω_{2} with the aim of either reinforcing the preexisting 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–18], the idea of considering separated substructures is newly addressed here. This means that mass (for dynamic 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–22]. Among the various rheological existing models of viscoelastic behavior, the fourparameter Zener fractionalderivative model is used in this work. It allows to efficiently represent the frequencydependence 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 posttreatment (plotting).
2 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 . 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 . 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.
Fig. 1 Geometry and finite element mesh of the problem composed of a given structure Ω_{1} and a design domain Ω_{2}. 
2.2 Viscoelastic model
In this work, the viscoelastic behavior of the dissipative material in Ω_{2} is modeled using a fourparameter Zener fractionalderivative 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:
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 subpart Ω_{2}.
2.3 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 wellknown 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 SIMPlaw 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}.
2.4 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 subdomain Ω_{2} is introduced in order to use a topology optimization procedure. Using a subdomain decomposition between domain Ω_{1} and Ω_{2}, the global discretized system equation (6) is written:
where 1, 2 and I denote respectively subdomain Ω_{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 subdomain Ω_{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
2.5 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.
2.6 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 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 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.
3 Topology optimization problem formulation
3.1 Minimization problem
In this work, the aim of the topology optimization is to find the material density set of design variables minimizing the objective function f (ω, x) for a given angular frequency ω according to constraints:
such that
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 modifiedSIMP 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 nmax in this work).
3.2 Compliance based objective function
3.2.1 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 denotes the transposeconjugate 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
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).
3.2.2 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, 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).
3.3 Displacement based objective function
3.3.1 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.
The derivative of d(ω, x) with respect to x_{e} is
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.
3.3.2 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:
4 Applications
4.1 Static reinforcement
4.1.1 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 8node 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 by and ; 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 solution , 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 solution , 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 nonconverged elements (material density between 0.01 and 0.99) is about 1% in both cases.
Fig. 2 3D cantilever beam (Ω_{1} ) reinforced by a design domain (Ω_{2}). 
Aluminum material properties.
Fig. 3 Optimal converged final shapes for the 3D cantilever beam: elements in Ω_{2} such that x_{e} > 0.99 are in red, elements in Ω_{1} are in blue. 
Fig. 4 Evolution of the objective functions c_{s} (x_{i})/c_{s} (x_{1}) and d_{s} (x_{i})/d_{s} (x_{1}) over the optimization iterations for the 3D cantilever beam where x_{1} is the design at the end of the first iteration while x_{i} is the design at iteration i. 
4.1.2 Application to a 3D CubeSatlike structure
The application addressed in this section concerns a CubeSatlike structure which is a standard format created by PuigSuari 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 8node 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 by and ; 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 solution , 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 solution , 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 nonconverged elements (material density between 0.01 and 0.99) is less than 1% in both cases.
Fig. 5 CubeSatlike structure (Ω_{1}) reinforced by a design domain (Ω_{2}). Note that half of Ω_{2} volume is represented with a translation from Ω_{1} in order to see its shape. 
Fig. 6 Optimal converged final shapes for the 3D CubeSatlike structure: elements in Ω_{2} such that x_{e} > 0.99 are in red, elements in Ω_{1} are in blue. 
Fig. 7 Evolution of the objective functions c_{s} (x_{i})/c_{s}(x_{1}) and d_{s}(x_{i})/d_{s}(x_{1}) over the optimization iterations for the 3D CubeSatlike structure where x_{1} is the design at the end of the first iteration while x_{i} is the design at iteration i. 
4.2 Vibration reduction
4.2.1 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 4node 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 by , and . 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 , as well as 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 solution can be better than an other in terms of minimizing the compliance. For instance, the material density solution given by the minimization of the compliance at 500 Hz exhibits the minimum of compliance at 500 Hz compared to the other solutions, even the solution computed with a full layer of viscoelastic material. The same conclusions can be done with solution at the frequency 1000 Hz; and with solution at the frequency 2000 Hz.
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 by , and . 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 , as well as 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.
Fig. 8 Optimization of an internal viscoelastic layer (Ω_{2}) of a 2D sandwich beam composed of 2 external aluminum layers (Ω_{1}). 
Fig. 9 Optimal solutions of the viscoelastic layer of the 2D sandwich beam (black: x_{e} > 0.99) for the compliance criteria, from bottom to top: , . 
Fig. 10 Dynamic compliance c(x) evolution over the optimization iterations for the 2D sandwich beam and for the three different harmonic excitations: c(500 Hz,x_{i}), c(1000 Hz,x_{i}) and c(2000 Hz,x_{i}). 
Viscoelastic material properties.
Number of elements in Ω_{2} in terms of their final material densities for the 2D sandwich beam.
Fig. 11 Frequency response functions of the 2D sandwich beam in terms of dynamic compliance: , and . 
Fig. 12 Optimal solutions of the viscoelastic layer of the 2D sandwich beam (black: x_{e} > 0.99) for the displacement criteria, from bottom to top: , , . 
Fig. 13 Displacement objective function d(x) evolution over the optimization iterations for the 2D sandwich beam and for the three different harmonic excitations: d(500 Hz, x_{i}), d(1000 Hz, x_{i}) and d(2000 Hz, x_{i}). 
Fig. 14 Frequency response functions of the 2D sandwich beam in terms of displacement objective function: , and 
4.2.2 Optimization of a viscoelastic layer for a 2D CubeSatlike 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 CubeSatlike 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 is imposed on the upper left corner at the angular velocity .
The structure is discretized with a total of 3222 quadrangular 4node 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 by , and . 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). These oscillations are more important for , nevertheless if only elements with x_{e} > 0.99 are considered the solution is acceptable. Moreover, for this solution, there are 369 elements with a material density less than 0.1; and so only 128 elements between 0.1 and 0.99. The frequency response functions of the structure are then computed using these three material densities. The compliance functions c(ω, ), c(ω, ) as well as c(ω, ) 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 solution 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 solution at the frequency 100 Hz; and with solution 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 by , and . 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 : 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 the solution, the last three iterations are plotted on Figure 22. Nevertheless if only elements with x_{e} > 0.99 are considered the solution is acceptable. The frequency response functions of the structure are then computed using these three material densities. The displacement functions d(ω, ), d(ω, ), as well as d(ω, ), are shown 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.
Fig. 15 Optimization of an internal viscoelastic layer (ω_{2}) for a 2D CubeSatlike structure (ω_{1}). 
Fig. 16 Deformed structure for the three target frequencies and for the compliance objective function (scale factor of the deformed mesh is 5000; black correspond to no displacement, yellow to maximum). 
Fig. 17 Optimal solutions of the viscoelastic layer of the 2D CubeSatlike structure (black: x_{e} > 0.99) for the compliance criteria, from left to right: , . 
Fig. 18 Dynamic compliance c(x) evolution over the optimization iterations for the 2D CubeSatlike structure and for the three different harmonic excitations: c(50 Hz, x_{i}), c(100 Hz, x_{i}) and c(200 Hz, x_{i}). 
Number of elements in Ω_{2} in terms of their final material densities for the 2D CubeSatlike structure.
Fig. 19 Frequency response functions of the 2D CubeSatlike structure in terms of dynamic compliance: c(ω, ), c(ω, ) and c(ω, ). 
Fig. 20 Optimal solutions of the viscoelastic layer of the 2D CubeSatlike structure (black: x_{e} > 0.99) for the displacement criteria, from left to right:, , . 
Fig. 21 Displacement objective function evolution over the optimization iterations for the 2D CubeSatlike structure and for the three different harmonic excitations: d(50 Hz, x_{i}), d(100 Hz, x_{i}) and d(200 Hz, x_{i}). 
Fig. 22 Last 3 iterations for the solution of the 2D CubeSatlike structure, from left to right: iterations 98, 99, and 100 (red: x_{e} = 1, blue: x_{e} = 0). 
Fig. 23 Frequency response functions of the 2D CubeSatlike structure in terms of displacement amplitude: d(ω, ), d(ω, ) and d(ω, ). 
5 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 fourparameter Zener fractionalderivative 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.
Declaration of conflicting interests
The authors declare no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.
References
 F. Rigaud, M. Charlotte, C. Kerdreux, P. Marechal, Multiobjective optimization of rotarywing aircrafts at the predesign stage, Mech. Ind. 15, 267–277 (2014) [Google Scholar]
 K. Suzuki, N. Kikuchi, A homogenization method for shape and topology optimization, Comput. Methods Appl. Mech. Eng. 93, 291–318 (1991) [CrossRef] [Google Scholar]
 O. Sigmund, A 99 line topology optimization code written in Matlab, Struct. Multidiscipl. Optim. 21, 120–127 (2001) [Google Scholar]
 M. Stolpe, K. Svanberg, An alternative interpolation scheme for minimum compliance topology optimization, Struct. Multidiscipl. Optim. 22, 116–124 (2001) [Google Scholar]
 D. Tchierniak, Topology optimization of resonating structures using SIMP method, Int. J. Numer. Methods Eng. 54, 1605–1622 (2002) [CrossRef] [Google Scholar]
 A. Takezawa, M. Daifuku, Y. Nakano, K. Nakagawa, T. Yamamoto, M. Kitamuraa, Topology optimization of damping material for reducing resonance response based on complex dynamic compliance, J. Sound Vibr. 365, 230–243 (2016) [Google Scholar]
 R.B. Haber, C.S. Jog, M.P. Bendsøe, A new approach to variabletopology shape design using a constraint on perimeter, Struct. Optim. 11, 1–12 (1996) [Google Scholar]
 L. Shu, M.Y. Wang, Z. Fang, Z. Ma, P. Wei, Level set based structural topology optimization for minimizing frequency response, J. Sound Vibr. 330, 5820–5834 (2011) [Google Scholar]
 J. Cao, K. Cai, P. Fei Wang, D. Yan, J. Shi, Multiple materials layout optimization in a layered structure, Mech. Ind. 17, 404 (2016) [Google Scholar]
 O.M. Silva, M.M. Neves, A. Lenzi, A critical analysis of using the dynamic compliance as objective function in topology optimization of onematerial structures considering steadystate forced vibration problems, J. Sound Vibr. 444, 1–20 (2019) [Google Scholar]
 S. Burri, Contributions a l’optimisation topologique de liaisons amortissantes pour des applications spatiales. These de doctorat. CNAM, HESAM Université (2020). https://tel.archivesouvertes.fr/tel03179845 [Google Scholar]
 M.P. Bendsøe, Optimal shape design as a material distribution problem, Struct. Optim. 1, 193–202 (1989) [Google Scholar]
 N. Zhou, G.I.N. Rozvany, The COC algorithm. Part II: Topological, geometrical and generalized shape optimization. Comput. Methods Appl. Mech. Eng. 89, 309–336 (1991) [CrossRef] [Google Scholar]
 M.P. Bendsøe, O. Sigmund, Topology Optimization – Theory, Methods and Applications (Springer, Berlin, Heidelberg, 2004) [Google Scholar]
 A. Diaz, N. Kikuchi, Solutions to shape and topology eigenvalue optimization problems using a homogenization method. Int. J. Numer. Methods Eng. 35, 1487–1502 (1992) [CrossRef] [Google Scholar]
 Z.D. Ma, N. Kikuchi, I. Hagiwara, Structural topology and shape optimization for a frequency response problem. Comput. Mech. 13, 157–174 (1993) [CrossRef] [MathSciNet] [Google Scholar]
 Z.D. Ma, H.S. Cheng, N. Kikuchi, Structural design for obtaining desired eigenfrequencies by using the shape and topology optimization method. Comput. Syst. Eng. 5, 77–89 (1994) [CrossRef] [Google Scholar]
 M. Bruggi, A. Talierci, Topology optimization of the fiberreinforcement retrofitting existing structures, Int. J. Solids Struct. 50, 121–136 (2013) [CrossRef] [Google Scholar]
 J.P. Kruth, M.C. Leu, T. Nakagawa, Progress in additive manufacturing and rapid prototyping, CIRP Ann. 47, 525–540 (1998) [CrossRef] [Google Scholar]
 R.L Bagley, P. Torvik, Fractional calculus – a different approach to the analysis of viscoelastically damped structures, AIAA J. 21, 741–748 (1983) [CrossRef] [Google Scholar]
 A. Germant, XLV. On fractional differentials, Philos. Mag. Ser. 1 25, 540–549 (1938) [Google Scholar]
 B. Morin, A. Legay, J.F. Deü, Reduced order models for dynamic behavior of elastomer damping devices, Finite Elem. Anal. Des. 143, 66–75 (2018) [CrossRef] [Google Scholar]
 L. Rouleau, A. Legay, J.F. Deü, Modelisation vibroacoustique de structures sandwich munies de matériaux viscoélastiques. These de Doctorat CNAM (2013) [Google Scholar]
 C. Geuzaine, J.F. Remacle, Gmsh: a 3D finite element mesh generator with builtin pre and postprocessing facilities. Int. J. Numer. Methods Eng. 79, 1309–1331 (2019) [Google Scholar]
 A.C. Galucio, J.F. Deü, R. Ohayon, Finite element formulation of viscoelastic sandwich beams using fractional derivative operators. Comput. Mech. 33, 282–291 (2004) [CrossRef] [Google Scholar]
 M.P. Bendsøe, N. Kikuchi, Generating optimal topologies in structural design using a homogenization method, Comput. Methods Appl. Mech. Eng. 71, 197–224 (1988) [Google Scholar]
 O. Sigmund, Morphologybased black and white filters for topology optimization. Struct. Multidiscipl. Optim. 33, 401–424 (2007) [Google Scholar]
 G.I.N. Rozvany, Structural Design via Optimality Criteria: The Prager Approach to Structural Optimization. (Kluwer Academic Publishers, Springer Dordrecht, 1989) [Google Scholar]
 K. Svanberg, The method of moving asymptotes – a new method for structural optimization. Int. J. Numer. Methods Eng. 24, 359–373 (1987) [CrossRef] [Google Scholar]
 Z.D. Ma, N. Kikuchi, I. Hagiwara, Structural topology and shape optimization for a frequency response problem, Comput. Mech. 13, 157–174 (1993) [CrossRef] [MathSciNet] [Google Scholar]
 G.H. Yoon, Structural topology optimization for frequency response problem using model reduction schemes, Comput. Methods Appl. Mech. Eng. 199, 1744–1763 (2010) [CrossRef] [Google Scholar]
 K.S. Yun, S.K. Youn, Topology optimization of viscoelastic damping layers for attenuating transient response of shell structures, Finite Element Anal. Des. 141, 154–165 (2018) [CrossRef] [Google Scholar]
 C.S. Jog, Topology design of structures subjected to periodic loading, J. Sound Vibr. 253, 687–709 (2002) [Google Scholar]
 K. Liu, A. Tovar, An efficient 3D topology optimization code written in Matlab, Struct. Multidiscipl. Optim. 50, 1175–1196 (2014) [Google Scholar]
 H. Heidt, J. PuigSuari, A.S. Moore, S. Nakasuka, R.J. Twiggs, CubeSat: a new generation of picosatellite for education and industry lowcost space experimentation, in 14th Annual AIAA/USU Conference on Small Satellites (2000) [Google Scholar]
 R. Hevner, W. Holemans, J. PuigSuari, R.J. Twiggs, An advanced standard for CubeSats, in 25th Annual AIAA/USU Conference on Small Satellites (2011) [Google Scholar]
Cite this article as: S. Burri and A. Legay, Static Reinforcement and Vibration Reduction of Structures using Topology Optimization, Mechanics & Industry 24, 14 (2023)
All Tables
Number of elements in Ω_{2} in terms of their final material densities for the 2D sandwich beam.
Number of elements in Ω_{2} in terms of their final material densities for the 2D CubeSatlike structure.
All Figures
Fig. 1 Geometry and finite element mesh of the problem composed of a given structure Ω_{1} and a design domain Ω_{2}. 

In the text 
Fig. 2 3D cantilever beam (Ω_{1} ) reinforced by a design domain (Ω_{2}). 

In the text 
Fig. 3 Optimal converged final shapes for the 3D cantilever beam: elements in Ω_{2} such that x_{e} > 0.99 are in red, elements in Ω_{1} are in blue. 

In the text 
Fig. 4 Evolution of the objective functions c_{s} (x_{i})/c_{s} (x_{1}) and d_{s} (x_{i})/d_{s} (x_{1}) over the optimization iterations for the 3D cantilever beam where x_{1} is the design at the end of the first iteration while x_{i} is the design at iteration i. 

In the text 
Fig. 5 CubeSatlike structure (Ω_{1}) reinforced by a design domain (Ω_{2}). Note that half of Ω_{2} volume is represented with a translation from Ω_{1} in order to see its shape. 

In the text 
Fig. 6 Optimal converged final shapes for the 3D CubeSatlike structure: elements in Ω_{2} such that x_{e} > 0.99 are in red, elements in Ω_{1} are in blue. 

In the text 
Fig. 7 Evolution of the objective functions c_{s} (x_{i})/c_{s}(x_{1}) and d_{s}(x_{i})/d_{s}(x_{1}) over the optimization iterations for the 3D CubeSatlike structure where x_{1} is the design at the end of the first iteration while x_{i} is the design at iteration i. 

In the text 
Fig. 8 Optimization of an internal viscoelastic layer (Ω_{2}) of a 2D sandwich beam composed of 2 external aluminum layers (Ω_{1}). 

In the text 
Fig. 9 Optimal solutions of the viscoelastic layer of the 2D sandwich beam (black: x_{e} > 0.99) for the compliance criteria, from bottom to top: , . 

In the text 
Fig. 10 Dynamic compliance c(x) evolution over the optimization iterations for the 2D sandwich beam and for the three different harmonic excitations: c(500 Hz,x_{i}), c(1000 Hz,x_{i}) and c(2000 Hz,x_{i}). 

In the text 
Fig. 11 Frequency response functions of the 2D sandwich beam in terms of dynamic compliance: , and . 

In the text 
Fig. 12 Optimal solutions of the viscoelastic layer of the 2D sandwich beam (black: x_{e} > 0.99) for the displacement criteria, from bottom to top: , , . 

In the text 
Fig. 13 Displacement objective function d(x) evolution over the optimization iterations for the 2D sandwich beam and for the three different harmonic excitations: d(500 Hz, x_{i}), d(1000 Hz, x_{i}) and d(2000 Hz, x_{i}). 

In the text 
Fig. 14 Frequency response functions of the 2D sandwich beam in terms of displacement objective function: , and 

In the text 
Fig. 15 Optimization of an internal viscoelastic layer (ω_{2}) for a 2D CubeSatlike structure (ω_{1}). 

In the text 
Fig. 16 Deformed structure for the three target frequencies and for the compliance objective function (scale factor of the deformed mesh is 5000; black correspond to no displacement, yellow to maximum). 

In the text 
Fig. 17 Optimal solutions of the viscoelastic layer of the 2D CubeSatlike structure (black: x_{e} > 0.99) for the compliance criteria, from left to right: , . 

In the text 
Fig. 18 Dynamic compliance c(x) evolution over the optimization iterations for the 2D CubeSatlike structure and for the three different harmonic excitations: c(50 Hz, x_{i}), c(100 Hz, x_{i}) and c(200 Hz, x_{i}). 

In the text 
Fig. 19 Frequency response functions of the 2D CubeSatlike structure in terms of dynamic compliance: c(ω, ), c(ω, ) and c(ω, ). 

In the text 
Fig. 20 Optimal solutions of the viscoelastic layer of the 2D CubeSatlike structure (black: x_{e} > 0.99) for the displacement criteria, from left to right:, , . 

In the text 
Fig. 21 Displacement objective function evolution over the optimization iterations for the 2D CubeSatlike structure and for the three different harmonic excitations: d(50 Hz, x_{i}), d(100 Hz, x_{i}) and d(200 Hz, x_{i}). 

In the text 
Fig. 22 Last 3 iterations for the solution of the 2D CubeSatlike structure, from left to right: iterations 98, 99, and 100 (red: x_{e} = 1, blue: x_{e} = 0). 

In the text 
Fig. 23 Frequency response functions of the 2D CubeSatlike structure in terms of displacement amplitude: d(ω, ), d(ω, ) and d(ω, ). 

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