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 Saint-Martin, 75003 Paris, France
* e-mail: 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 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.
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, 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–11].
The present article uses an algorithm based on the well known SIMP-method [12–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–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 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).
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 xe. 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 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 i2 = −1. Both moduli can be expressed as a function of the static stiffness E0, 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 sub-part Ω2.
2.3 Implementation of design variables
For the purpose of topology optimization, the design domain Ω2 contains a set of design variables x = [x1, x2, …, xn]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 ≤ xe ≤ 1 where 0 is the lower bound of xe 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 xe such that:
where E*(ω) is given by equation (3), Emin = ϵE0 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 xe 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 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 K2(ω, 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 S1(ω) does not depend on x and represents the contribution of Ω1, while S2(ω) depends on x and represents the contribution of the design domain Ω2.
The global dynamic equation (Eq. (6)) then becomes
By denoting S(ω, x) = S1(ω) + S2(ω, 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 xe. Since only matrix S2(ω, x) is depending on design variables, the derivative of equation (11) with respect to the design variable xe leads to the following equation:
The sensitivity of the solution with respect to the design variable xe is then given by:
where U (x) is the solution of equation (11). In this last equation, the derivative of the S2(ω, x) matrix needs to be computed.
2.6 Derivatives of the stiffness and mass matrices of the design domain
The derivative of S2(ω, x) with respect to xe involves the derivative of the stiffness and mass matrices of Ω2, denoted respectively by K2(ω, x) and M2(x):
The stiffness matrix K2(ω, x) can be expressed as the assembly of the elemental stiffness matrices kj(ω, xj) over the n elements of Ω2-domain mesh:
The elemental stiffness matrix kj(ω, xj) can be written as
where is the stiffness matrix of element j computed with a unit Young modulus. Thus, the derivative of K2(ω, x) with respect to xe 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 M2(x) is expressed as the assembly of the elemental mass matrices mj(xj) over the n elements of Ω2-domain mesh:
The elemental mass matrix mj(xj) can be written as
where is the mass matrix of element j computed with a unit volumetric mass density. Thus, the derivative of M2(x) with respect to xe 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, v0 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 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 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 xe is then given by:
This last expression involves the sensitivity of g(ω, x) with respect to xe, 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 xe. Once the sensitivity of the g(x) function with respect to xe 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 cs(x)) can be derived from the previous section by taking ω = 0:
where
In this case, since S2(0, x) is a real matrix. The derivative of g(0, x) with respect to xe becomes:
and the static compliance sensitivity is then given by
where E*(0) = E0 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 bi 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 xe 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 ds(x)) can be derived from the previous section by taking ω = 0:
where U (x) is the solution of:
The sensitivity of ds(x) with respect to xe 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 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 cs(x) and displacement based function ds(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 nmax of iterations is set to 100.
The optimal final shapes obtained by minimizing cs(x) and ds(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 xe > 0.99
6420 elements such that xe < 0.01
80 elements such that 0.01 < xe < 0.99
In the final shape for solution , there are:
2140 elements such that xe > 0.99
6400 elements such that xe < 0.01
100 elements such that 0.01 < xe < 0.99
The quantity of non-converged 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 xe > 0.99 are in red, elements in Ω1 are in blue. |
Fig. 4 Evolution of the objective functions cs (xi)/cs (x1) and ds (xi)/ds (x1) over the optimization iterations for the 3D cantilever beam where x1 is the design at the end of the first iteration while xi is the design at iteration i. |
4.1.2 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 cs(x) and displacement based objective functions ds(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 nmax is set to 200.
The optimal final shapes obtained by minimizing cs(x) and ds(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 nmax 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 xe > 0.99
7064 elements such that xe < 0.01
80 elements such that 0.01 < xe < 0.99
In the final shape for solution , there are:
2368 elements such that xe > 0.99
7072 elements such that xe < 0.01
64 elements such that 0.01 < xe < 0.99
The quantity of non-converged elements (material density between 0.01 and 0.99) is less than 1% in both cases.
Fig. 5 CubeSat-like 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 CubeSat-like structure: elements in Ω2 such that xe > 0.99 are in red, elements in Ω1 are in blue. |
Fig. 7 Evolution of the objective functions cs (xi)/cs(x1) and ds(xi)/ds(x1) over the optimization iterations for the 3D CubeSat-like structure where x1 is the design at the end of the first iteration while xi 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 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 nmax 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 nmax=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 (xe = 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 (xe = 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 visco-elastic layer of the 2D sandwich beam (black: xe > 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,xi), c(1000 Hz,xi) and c(2000 Hz,xi). |
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 visco-elastic layer of the 2D sandwich beam (black: xe > 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, xi), d(1000 Hz, xi) and d(2000 Hz, xi). |
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 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 nmax 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 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 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 xe > 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 (xe = 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 xe < 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 xe > 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 (xe = 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 CubeSat-like 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 CubeSat-like structure (black: xe > 0.99) for the compliance criteria, from left to right: , . |
Fig. 18 Dynamic compliance c(x) evolution over the optimization iterations for the 2D CubeSat-like structure and for the three different harmonic excitations: c(50 Hz, xi), c(100 Hz, xi) and c(200 Hz, xi). |
Number of elements in Ω2 in terms of their final material densities for the 2D CubeSat-like structure.
Fig. 19 Frequency response functions of the 2D CubeSat-like structure in terms of dynamic compliance: c(ω, ), c(ω, ) and c(ω, ). |
Fig. 20 Optimal solutions of the viscoelastic layer of the 2D CubeSat-like structure (black: xe > 0.99) for the displacement criteria, from left to right:, , . |
Fig. 21 Displacement objective function evolution over the optimization iterations for the 2D CubeSat-like structure and for the three different harmonic excitations: d(50 Hz, xi), d(100 Hz, xi) and d(200 Hz, xi). |
Fig. 22 Last 3 iterations for the solution of the 2D CubeSat-like structure, from left to right: iterations 98, 99, and 100 (red: xe = 1, blue: xe = 0). |
Fig. 23 Frequency response functions of the 2D CubeSat-like 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 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.
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 rotary-wing aircrafts at the predesign stage, Mech. Ind. 15, 267–277 (2014) [CrossRef] [EDP Sciences] [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) [CrossRef] [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) [CrossRef] [Google Scholar]
- R.B. Haber, C.S. Jog, M.P. Bendsøe, A new approach to variable-topology shape design using a constraint on perimeter, Struct. Optim. 11, 1–12 (1996) [CrossRef] [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) [CrossRef] [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) [CrossRef] [EDP Sciences] [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 one-material structures considering steady-state forced vibration problems, J. Sound Vibr. 444, 1–20 (2019) [CrossRef] [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.archives-ouvertes.fr/tel-03179845 [Google Scholar]
- M.P. Bendsøe, Optimal shape design as a material distribution problem, Struct. Optim. 1, 193–202 (1989) [CrossRef] [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) [CrossRef] [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 vibro-acoustique de structures sandwich munies de matériaux viscoélastiques. These de Doctorat CNAM (2013) [Google Scholar]
- C. Geuzaine, J.-F. Remacle, Gmsh: a 3-D finite element mesh generator with built-in pre- and post-processing 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, Morphology-based black and white filters for topology optimization. Struct. Multidiscipl. Optim. 33, 401–424 (2007) [CrossRef] [Google Scholar]
- G.I.N. Rozvany, Structural Design via Optimality Criteria: The Prager Approach to Structural Optimization. (Kluwer Academic Publishers, Springer Dordrecht, 1989) [CrossRef] [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) [CrossRef] [Google Scholar]
- K. Liu, A. Tovar, An efficient 3D topology optimization code written in Matlab, Struct. Multidiscipl. Optim. 50, 1175–1196 (2014) [CrossRef] [Google Scholar]
- H. Heidt, J. Puig-Suari, A.S. Moore, S. Nakasuka, R.J. Twiggs, CubeSat: a new generation of picosatellite for education and industry low-cost space experimentation, in 14th Annual AIAA/USU Conference on Small Satellites (2000) [Google Scholar]
- R. Hevner, W. Holemans, J. Puig-Suari, 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 CubeSat-like 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 xe > 0.99 are in red, elements in Ω1 are in blue. |
|
In the text |
Fig. 4 Evolution of the objective functions cs (xi)/cs (x1) and ds (xi)/ds (x1) over the optimization iterations for the 3D cantilever beam where x1 is the design at the end of the first iteration while xi is the design at iteration i. |
|
In the text |
Fig. 5 CubeSat-like 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 CubeSat-like structure: elements in Ω2 such that xe > 0.99 are in red, elements in Ω1 are in blue. |
|
In the text |
Fig. 7 Evolution of the objective functions cs (xi)/cs(x1) and ds(xi)/ds(x1) over the optimization iterations for the 3D CubeSat-like structure where x1 is the design at the end of the first iteration while xi 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 visco-elastic layer of the 2D sandwich beam (black: xe > 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,xi), c(1000 Hz,xi) and c(2000 Hz,xi). |
|
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 visco-elastic layer of the 2D sandwich beam (black: xe > 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, xi), d(1000 Hz, xi) and d(2000 Hz, xi). |
|
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 CubeSat-like 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 CubeSat-like structure (black: xe > 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 CubeSat-like structure and for the three different harmonic excitations: c(50 Hz, xi), c(100 Hz, xi) and c(200 Hz, xi). |
|
In the text |
Fig. 19 Frequency response functions of the 2D CubeSat-like structure in terms of dynamic compliance: c(ω, ), c(ω, ) and c(ω, ). |
|
In the text |
Fig. 20 Optimal solutions of the viscoelastic layer of the 2D CubeSat-like structure (black: xe > 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 CubeSat-like structure and for the three different harmonic excitations: d(50 Hz, xi), d(100 Hz, xi) and d(200 Hz, xi). |
|
In the text |
Fig. 22 Last 3 iterations for the solution of the 2D CubeSat-like structure, from left to right: iterations 98, 99, and 100 (red: xe = 1, blue: xe = 0). |
|
In the text |
Fig. 23 Frequency response functions of the 2D CubeSat-like structure in terms of displacement amplitude: d(ω, ), d(ω, ) and d(ω, ). |
|
In the text |
Current usage metrics show cumulative count of Article Views (full-text 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 48-96 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.