# Topology Optimization of Free Damping Treatments on Plates Using Level Set Method

Authors: Jian Pang; Weiguang Zheng; Liang Yang; Yuping Wan; Qibai Huang; Rongjiang Tang
Journal: Shock and Vibration
 (2020)

---

## Abstract

Application of level set method to optimize the topology of free damping treatments on plates is investigated. The objective function is defined as a combination of several desired modal loss factors solved by the finite element-modal strain energy method. The finite element model for the composite plate is described as combining the level set function. A clamped rectangle composite plate is numerically and experimentally analyzed. The optimized results for a single modal show that the proposed method has the possibility of nucleation of new holes inside the material domain, and the final design is insensitive to initial designs. The damping treatments are guided towards the areas with high modal strain energy. For the multimodal case, the optimized result matches the normalized modal strain energy of the base plate, which would provide a simple implementation way for industrial application. Experimental results show good agreements with the proposed method. The experimental results are in good agreement with the optimization results. It is very promising to see that the optimized result for each modal has almost the same damping effect as that of the full coverage case, and the result for multimodal gets moderate damping at each modal.


---

## Body

## 1. Introduction

In the automotive industry, the door, roof, dash, floor, and cab back panels of automobiles are always treated with damping materials to reduce the structure-borne noise. The effectiveness of damping treatments depends upon design parameters such as damping material types, locations, and size of the treatment. As the weight often plays a key role in the performance and cost in the industry, it is essential to get the optimal distribution of the damping treatments with a limited coverage rate.Traditionally, experimental techniques conduct laser vibro-meter type tests on the BIW structure or full-vehicle prototypes excited by a shaker at each suspension to generate velocity contours as optimal damping treatment configuration, which is prohibitively time consuming and expensive. Conventional CAE methodology employs the modal strain energy contour based on finite element modeling to determine the distribution of the damping treatment. The method is easy to implement and cost-effective, while the optimal results are not very accurate as the effectiveness of the damping treatment was not taken into account in the FE modeling, and the optimal configuration is determined by the engineer’s selection by trial-and-error based on the strain energy contour [1].During the recent years, extensive efforts have been exerted to optimally design damping treatments for vibration and noise reduction, resulting in a large number of studies in the field. Lall et al. [2] carried out optimum design studies for the parameters like the densities, thicknesses, and temperature of viscoelastic materials for a fully passive constrained layer damping (PCLD) covered plate with objective functions such as maximizing the modal loss factor and minimizing the displacement response. Marcelin et al. [3] used the finite element method (FEM) and genetic algorithm (GA) to maximize the damping factor for the partially covered beam. The design variables were the dimensions and locations of the patches. Zheng et al. [4] studied the optimal layout of PCLD patches on beams for minimization vibration energy based on the analytical model using GA with penalty functions for constraints.All of these studies mentioned above which performed shape or size optimization are at the level of macroscopic design, using a macroscopic definition of geometry given by, for example, dimensions or boundaries. One of the limitations of these conventional shape or size optimizations by boundary variations is that the topology of the structure is fixed during the iterative design process. A very interesting idea was proposed by Alvelid [5] who suggested a shape optimization method by adding the PCLD to the base structure piece by piece at permissible positions in a mosaique manner with the objective function as the square of the surface velocity. Though the shape is not limited in the optimization procedure, the optimized results may not be the best in a global perspective as the adding manner of the PCLD element is limited by predetermined distribution. Koruk and Sanliturk [6] proposed an optimization procedure based on the big bang-big crunch and modal strain energy (MSE) methods, which can be applied to general structures for locations and thicknesses of damping treatments.Since there is a limited type for the PCLD products on the market, the main problem is to determine the shape, location, and sizing of the PCLD patch with certain specification of materials and dimensions. Its very nature is to formulate this optimization problem as a topological design optimization. The idea of using topology optimization to design continuum structures was first introduced about twenty years ago by Bendsøe and Kikuchi [7]. Mohammed [8] dealt with the optimal distribution for viscoelastic materials on sandwiched damping beams or plates using the inverse homogenization method to maximize the modal loss factors. Ling et al. [9] extended Mohanmmed’s work to optimize the viscoelastic materials distribution in CLD using the solid isotropic material with penalization (SIMP) model and the method of moving asymptote (MMA) approach. Lots of optimization approaches for vibration and noise control using damping treatments [10, 11] are carried out based on the SIMP method due to its easy implementation.Unfortunately, the above topology optimization methods tend to suffer from numerical instability problems such as mesh dependency, checkerboard patterns, and gray scales. A different approach is used in level set-based structural optimization methods that have been proposed as a new type of structural optimization. Sethian and Wiegmann [12] introduced the level set method for the first time in the field of structural topology optimization. They used the function of stress as the velocity field to promote the structure boundary when studying the algorithm problem of the elastic structure boundary design. The main idea of the level set method is to introduce the motion interface in the N-dimensional space as the zero contour or the zero isosurface of the level set function in the N + 1-dimensional space by introducing the level set function, which is positive in the presence of the material, zero in the material boundary, and negative in the material-free area. The method describes the geometric boundary of the structure by the evolution of the level set function, which can simultaneously optimize the shape and topology. The level set method has the advantages of obtaining the smooth boundary, not relying on the meshing density.Ansari et al. [13] first optimized the modal loss factor as the optimization objective function and used the variational derivative to obtain the velocity field to optimize the PCLD layout on the plate structure. The classical level set optimization method was adopted, which has the drawbacks that the optimization result is dependent on the initial design and has difficulty in solving Hamilton–Jacobi partial differential equation (H-J PDE) numerical problems.The parameterized level set can well compensate for the shortcomings of numerical calculation difficulties [14]. The main idea of the parameterized level set method is to interpolate the level set function by a linear combination of radial basis function (RBF) and coefficients and participate in topology optimization instead of the original level set function. Since RBF is only related to space coordinates, the evolution of the level set function is transformed into an iterative update of the interpolation coefficient, which directly avoids the solution of the original H-J PDE, and effectively solves the problem of the numerical solution of the traditional level set method. Parametric modeling using RBF can maintain a relatively smooth level set function without having to reinitialize the level set function as often as the traditional method. Since all interpolation points are involved in the calculation of the velocity field, the velocity field is naturally extended to the entire design domain, which gives the parametric level set method the ability to autonomously open holes in the structure. Wang et al. [15] interpolated the level set function using the globally supported radial basis function (GSRBF) to avoid solving the complex H-J PDEs directly. The nucleation of new holes inside the material domain becomes possible, and thus the final design is fairly insensitive to initial designs. Furthermore, since the time step size and the number of RBF knots are determined by the requirement of optimization convergence rather than the need to solving the H-J equation, the constraint from the CFL condition on the temporal and spatial discretization can be relaxed. Wei et al. [16] proposed an 88-line MATLAB code for a parameterized level set method using multiquadric (MQ) RBFs with an approximate reinitialization operation.The main goal of this paper is to investigate the use of topology optimization based on a level set-based parameterization method to optimize free damping treatments with partial coverage to improve the damping characteristics of vibrating plates. The material of the paper is organized as follows. In Section2, we describe the topology optimization problem of the free damping treatment plate for the maximum modal loss factor. An introduction of FE modeling of the composite plate combining level set function is present in Section 3. The level set method of the proposed problem is deduced in Section 4. In Sections 5 and 6, taking a simple clamped rectangle plate as an example, numerical and experimental results are analyzed. Section 7 concludes the paper.


## 2. Problem Statement

Figure1 shows a base plate covered by the viscoelastic layer with an initial shape Hϕ0 to suppress flexural vibrations. Hϕ is the Heaviside function defined as follows to represent the existence of the damping layer at a point of the level set function value ϕ:(1)Hϕ=0,ϕ<0,1,ϕ≥0.Figure 1
Damping treatments with a shape ofHϕ0.The objective of the paper is to find the optimal shape of the viscoelastic layerHϕoptimum for the maximum vibration suppression. The coverage rate of the viscoelastic layer is constrained by a given value and the specification of the viscoelastic layer is predetermined.A linear and viscoelastic, frequency-independent, complex constant modulus is supposed to describe the properties of the viscoelastic materials as follows:(2)Ev∗=Ev1+jηv,where Ev and ηv are the elastic modulus and the material loss factor, respectively. The dynamic equations for free vibration of the structure with viscoelastic materials has the finite element form as(3)MX¨+KR+jKIX=0,where M is the global mass matrix, KR is the real part of the global stiffness matrix, KI is the imaginary part of the global stiffness matrix, and X is the nodal displacement vector. Then, the rth modal loss factor, which represents the vibration energy dissipation ratio of the modal, can be derived using the FE-MSE method as follows [11]:(4)ηr=φrTKIφrφrTKRφr,where φr is the rth real mode shape vector of the associated undamped system. In order to get an averaged effect over broadband frequencies, we define the objective function J as(5)J=∑rarηr,where ar is the associated weighting factor and ∑rar=1.In addition, the constraint should be considered to limit the consumption of damping material treatments, and the volume fraction is limited here. Then, the mathematical formulation based on the level set method of the optimization problem is defined as(6)MaxmizeϕJϕ=∫Ω∑rarφrTϕKIϕφrϕφrTϕKRϕφrϕdΩs.t.∫ΩHϕdΩΩ=βKRϕ−λrϕMϕφrϕ=0,where Ω represents the design domain and β is the fraction of the damping layer coverage.


## 3. Finite Element Modeling

WhenHϕ=0, the neutral plane of the composite is the midplane of the base plate as no damping layer is covered. When Hϕ=1, as shown in Figure 2, both the base plate and the damping layer contribute to the elasticity. The complex position of the neutral axis hN∗ is represented as follows [17]:(7)hN∗=Eh2+Ev∗h12h+h12Eh+Ev∗h1,where E is the elastic modulus of the base plate and h and h1 are the thickness of the base layer and the viscoelastic layer, respectively.Figure 2
The neutral axis of the composite layer.The composite layer can be meshed with the rectangular nonconforming plate-bending element. Based on the classical plate and finite element theory, the displacement field of the element can be written as follows [11]:(8)u=TNue,where(9)T=100−z00100−z00100,N=N1N2N3N4,Ni=ψi00000ψi00000gi1gi2gi300∂gi1/∂x∂gi2/∂x∂gi3/∂x00∂gi1/∂y∂gi2/∂y∂gi3/∂y,in which ψii=1,2,3,4 and gijj=1,2,3 are the linear and nonconforming Hermite cubic interpolation functions, respectively, which can be expressed as follows:(10)ψi=141+ξξi1+ηηi,gi1=181+ξξi1+ηηi2+ξξi+ηηi−ξ2−η2,gi2=18ξi1+ξξi2ξξi−11+ηηia,gi3=18ηi1+ηηi2ηηi−11+ξξib,where a and b are the half lengths of the rectangular element along the x and y directions, respectively, and ue is the nodal displacement vector in the neutral plane and each node has three translations and two rotation degrees of freedom.Then, the mass matrix and the stiffness matrix of the composite element can be derived from the variational principle as follows [11, 17]:(11)Me=∫−11∫−11ρNT∫−h/2h/2TTTdzNJdξdη+Hϕ∫−11∫−11ρ1NT∫−h1/2h1/2TTTdzNJacdξdη,Ke=1−Hϕ∫−11∫−11∫−h/2h/2BTQBJacdξdηdz+Hϕ∫−11∫−11∫−h/2−hN∗h/2−hN∗BTQBJacdξdηdz+∫−11∫−11∫h/2−hN∗h/2−hN∗+h1BTQ1BJacdξdηdz,where(12)B=B1B2B3B4=A−zC,Bi=Ai−zCi,Ai=∂ψi/∂x00000∂ψi/∂y000∂ψi/∂y∂ψi/∂x000,Ci=00∂2gi1/∂x2∂2gi2/∂x2∂2gi3/∂x200∂2gi1/∂y2∂2gi2/∂y2∂2gi3/∂y200∂2gi1/∂x∂y∂2gi2/∂x∂y∂2gi3/∂x∂y,Jac is the Jacobian matrix, ρ and ρ1 are the density of the base layer and the viscoelastic layer, respectively, and Q and Q1 are the elastic matrixes of the base layer and the viscoelastic layer, respectively.Then, the global mass matrix and the stiffness matrix of the composite structure in equation (3) can be represented as follows:(13)M=∑j=1nMje,KR=∑j=1nKR,je,KI=∑j=1nKI,je,where KRe and KIe are the real and imaginary parts of Ke and n is the number of the elements.


## 4. Parameterized Level Set Method

### 4.1. Velocity Field

In this paper, the Lagrangian formulation is applied by means of the Lagrangian multiplierλ to combine the volume constraint and the objective in the optimization problem, as expressed as follows:(14)Jϕ,λ=∫Ω∑r−arφrTKIφrφrTKRφr−λHϕ−βΩdΩ.When the derivative of equation (14) is zero, the objective function gets the minimum value. The variation of the Heaviside function can be expressed as(15)δHϕ=∂Hϕ∂ϕδϕ=δϕδϕ=δϕϕ=0,where δϕ is the unit pulse function and δϕ is the variation of ϕ. Since only the normal velocity changes the boundary shape, the tangential velocity is usually ignored and equation (15) can be written as(16)δHϕ=δϕϕ=0=∇ϕδl,where δl is an infinitesimal amount along the normal direction which is expressed as n=∇ϕ/∇ϕ. The variation of δl for ϕ can be expressed as(17)δϕJ=∫Ω∑r−arφrTKIφrφrTKRφr−λδϕ∇ϕδldΩ.The sensitivity of the stiffness matrix aboutϕ can be obtained from the finite element model. According to the variational principle, at the extreme point, the following equation can be derived from equation (17):(18)∑r−arφrTKIφrφrTKRφr−λδϕ∇ϕ=0.The numerical solution form of the level set equation can be written as(19)∂ϕ∂t−∑r−arφrTKIφrφrTKRφr−λδϕ︸Vn∇ϕ=0.Therefore, the sensitivity of the objective function about shape, that is, the velocity field of the element can be expressed as(20)Vne=∑r−arφrTKIeφrφrTKRφr−λeδϕ.Please note that(21)ηre=φrTKIeφrφrTKRφr,can be called the element modal loss factor.


### 4.2. MQ RBFs

The MQ function is proposed and used as a common form of radial basis function, which had been performing well in function approximation, curve fitting, partial differential equation solving, etc. MQ functions can be expressed as(22)gix=x−xi2+c2,xi∈D,where xixi,yi is the coordinate of the ith interpolation point and c is the shape parameter, a constant with a small value. The MQ function and its partial derivatives for x and y are shown in Figure 3.Figure 3
MQ function and its partial derivatives aboutx and y: (a) gix, (b) ∂gix/∂x, and (c) ∂gix/∂y.
(a)
(b)
(c)Interpolating the level set function withn fixed points, the parameterized level set function can be expressed as(23)ϕx,t=∑i=1nαitgit+px,t,where αit is the interpolation coefficient of the MQ function at the ith interpolation node at t time and px,t is a linear polynomial used to express the linear and constant parts of ϕx,t to ensure the correctness of the solution. For two-dimensional problems, px,t can be expressed as(24)px,t=p0t+p1tx+p2ty,where p0t, p1t, and p2t are the polynomial coefficients.In order to guarantee the uniqueness of the solution of the level set function, the interpolation coefficients must satisfy the following constraints:(25)∑i=1nαit=0,∑i=1nαitxi=0,∑i=1nαityi=0.Equations (23)–(25) can be expressed in the matrix form as follows:(26)Gαt=ϕt,where(27)G=APPT0,A=g1x1g2x1…gnx1g1x2g2x2…gnx2⋮⋮⋱⋮g1xng2xn…gnxn,P=1x1y1⋮⋮⋮1xnyn,αt=α1t…αntp0tp1tp2tT,ϕt=ϕx1,t…ϕxn,t000T.Since the matrixG is theoretically invertible, the interpolation coefficient can be expressed as(28)αt=G−1ϕt.Then, the level set function (23) can be written as(29)ϕx,t=gxTαt,where(30)gx=g1x…gnx1xyT.


### 4.3. Level Set Function Reinitialization

The fast marching method and the partial differential method are commonly used in the reinitialization of the level set function. The purpose is to ensure that the level set function is a symbol distance function at least near the boundary, that is,∇ϕ=1. Both methods have a large computational burden. The approximate reinitialization method chosen in this paper controlled the value of the level set function near the boundary to be within a stable range [16]. For the entire design domain, reinitialization can be expressed as(31)ϕu=ϕmean∇ϕ10,∇ϕ20,…,∇ϕr0,where ϕr0 is the level set function value of the rth point near the level set zero contour and ∇ϕr0 is the second order norm of the gradient of ϕr0. mean∇ϕ10,∇ϕ20,…,∇ϕr0 is the average value of ∇ϕ10,∇ϕ20,…,∇ϕr0. Reinitialization is shown in Figure 4.Figure 4
Level set function reinitialization.There is a linear relationship between the level set functionϕ and the interpolation coefficient α, so after the level set function has been reinitialized, the interpolation coefficient appears as(32)αu=αmean∇ϕ10,∇ϕ20,…,∇ϕr0.Compared with the common reinitialization method, the method used in this paper was simple and would not hinder the ability of autonomous opening and can fully retain the boundary information and only change the relative size of the value.


### 4.4. Numerical Solution

After the level set function was interpolated with the MQ radial basis function, the numerical calculation was performed. Substituting formula (30) into the level set equation (20), the “level set equation” for the radial basis function interpolation coefficient can be expressed as(33)gxdαtdt−Vn∇gxαt=0,where(34)∇gxαt=∂gx∂xαt2+∂gx∂yαt2.For this time-varying interpolation problem, in order to ensure that the interpolation coefficients can be solved without being limited to the positive definition of the condition of the MQ function, a boundary constraint needs to be introduced.(35)∑i=1nα˙it=0,∑i=1nα˙itxi=0,∑i=1nα˙ityi=0.So far, the time-dependent H-J PDEs were discrete into a set of coupled ordinary differential equations. Equation (33) could also be written as(36)G∂α∂t−Bα,t=0,where(37)Bα,t=Vnx1,t∇gx1αt⋮Vnxn,t∇gxnαt000.The set of coupled nonlinear ordinary differential equations can be numerically solved by the first-order forward Euler method, and the approximate solution can be expressed as(38)αti+1=αti+ΔtG−1Bαti,ti,where Δt is the size of time step and ti is the ith time step.At the same time, theδ function is used to control the fast boundary growth rate during the solution process, which can be expressed as(39)δϕ=0,ϕ>Δ,τΔ1−ϕ2Δ2,ϕ≤Δ,where Δ is a parameter used to control the upper and lower limits of the level set function value. When the level set function value is higher than Δ or lower than −Δ, its normal velocity will be ignored, so as to avoid the level set function value tending to infinity. Generally, Δ is selected to be more than 5 times the size of the mesh. τ is usually set to 0.75 [16] and also can be adjusted to avoid instability during the optimization (see discussions in Section 5.1).Combined with the reinitialization method, the iterative equation (38) of the level set function can be expressed as(40)αti+1=αuti+ΔtG−1B^αuti,ti.Among them,δϕ can replace ∇ϕ, then there was(41)B^αuti,ti=Vnx1,tiδϕx1,αuti⋮Vnxn,tiδϕxn,αuti000.In addition, since the time step and the number of interpolation nodes of the radial basis function are determined by the requirements of the optimization convergence rather than the need to solve the H-J equation, the time and space discretization requirements of the CFL condition can be relaxed. Using equation (40), the interpolation coefficient is updated until the convergence condition is satisfied, then the level set function is interpolated, and the current optimized structure topology can be obtained by taking zero contour.


### 4.5. Lagrange Multiplier

In the calculation of equation (40), the iteration of the Lagrange multiplier λ can be expressed as(42)λk+1=μGk,k≤nR,λk+γkG,k>nR,where μ and γ are the parameters controlling the Lagrange multiplier, k represents the kth iteration, and nR is the minimum number of iterations to be optimized, and the iteration of γ is expressed as(43)γk+1=minγk+Δγ,γmax,k>nR.When the volume fraction of the initial design is not equal to the volume fraction of the constraint, the volume constraint can be relaxed during the initial iteration, expressed as follows:(44)Gk=∫ΩHϕdΩ−V0−V0−αΩknR,k>nR,where V0 is the initial volume.


### 4.6. Flow Chart of the Algorithm

The flow chart of the damping layer parameterization level set optimization algorithm based on the modal loss factor is shown in Figure5 and is described as follows:(1)
Given the initial design of the structure and meshing, initialize the level set function and radial basis function.(2)
Run a finite element analysis to calculate the modal loss factor of the current structure.(3)
Determine whether the number of iterations exceeds the specified minimum number of iterations, determine whether the current volume fraction reaches the constraint condition, and determine whether the difference between the current modal loss factor and the modal loss factor of the first five calculations is sufficiently small. When all three are satisfied, the loop ends, otherwise, the following steps are continued.(4)
Update the Lagrange multiplier according to the current volume fraction and related parameters, and calculate the evolution velocity field and update the level set function.(5)
Update the radial basis function interpolation coefficients and return to step (2).Figure 5
Flow chart of the parameterized level set method for optimizing the modal loss factor.


## 4.1. Velocity Field

In this paper, the Lagrangian formulation is applied by means of the Lagrangian multiplierλ to combine the volume constraint and the objective in the optimization problem, as expressed as follows:(14)Jϕ,λ=∫Ω∑r−arφrTKIφrφrTKRφr−λHϕ−βΩdΩ.When the derivative of equation (14) is zero, the objective function gets the minimum value. The variation of the Heaviside function can be expressed as(15)δHϕ=∂Hϕ∂ϕδϕ=δϕδϕ=δϕϕ=0,where δϕ is the unit pulse function and δϕ is the variation of ϕ. Since only the normal velocity changes the boundary shape, the tangential velocity is usually ignored and equation (15) can be written as(16)δHϕ=δϕϕ=0=∇ϕδl,where δl is an infinitesimal amount along the normal direction which is expressed as n=∇ϕ/∇ϕ. The variation of δl for ϕ can be expressed as(17)δϕJ=∫Ω∑r−arφrTKIφrφrTKRφr−λδϕ∇ϕδldΩ.The sensitivity of the stiffness matrix aboutϕ can be obtained from the finite element model. According to the variational principle, at the extreme point, the following equation can be derived from equation (17):(18)∑r−arφrTKIφrφrTKRφr−λδϕ∇ϕ=0.The numerical solution form of the level set equation can be written as(19)∂ϕ∂t−∑r−arφrTKIφrφrTKRφr−λδϕ︸Vn∇ϕ=0.Therefore, the sensitivity of the objective function about shape, that is, the velocity field of the element can be expressed as(20)Vne=∑r−arφrTKIeφrφrTKRφr−λeδϕ.Please note that(21)ηre=φrTKIeφrφrTKRφr,can be called the element modal loss factor.


## 4.2. MQ RBFs

The MQ function is proposed and used as a common form of radial basis function, which had been performing well in function approximation, curve fitting, partial differential equation solving, etc. MQ functions can be expressed as(22)gix=x−xi2+c2,xi∈D,where xixi,yi is the coordinate of the ith interpolation point and c is the shape parameter, a constant with a small value. The MQ function and its partial derivatives for x and y are shown in Figure 3.Figure 3
MQ function and its partial derivatives aboutx and y: (a) gix, (b) ∂gix/∂x, and (c) ∂gix/∂y.
(a)
(b)
(c)Interpolating the level set function withn fixed points, the parameterized level set function can be expressed as(23)ϕx,t=∑i=1nαitgit+px,t,where αit is the interpolation coefficient of the MQ function at the ith interpolation node at t time and px,t is a linear polynomial used to express the linear and constant parts of ϕx,t to ensure the correctness of the solution. For two-dimensional problems, px,t can be expressed as(24)px,t=p0t+p1tx+p2ty,where p0t, p1t, and p2t are the polynomial coefficients.In order to guarantee the uniqueness of the solution of the level set function, the interpolation coefficients must satisfy the following constraints:(25)∑i=1nαit=0,∑i=1nαitxi=0,∑i=1nαityi=0.Equations (23)–(25) can be expressed in the matrix form as follows:(26)Gαt=ϕt,where(27)G=APPT0,A=g1x1g2x1…gnx1g1x2g2x2…gnx2⋮⋮⋱⋮g1xng2xn…gnxn,P=1x1y1⋮⋮⋮1xnyn,αt=α1t…αntp0tp1tp2tT,ϕt=ϕx1,t…ϕxn,t000T.Since the matrixG is theoretically invertible, the interpolation coefficient can be expressed as(28)αt=G−1ϕt.Then, the level set function (23) can be written as(29)ϕx,t=gxTαt,where(30)gx=g1x…gnx1xyT.


## 4.3. Level Set Function Reinitialization

The fast marching method and the partial differential method are commonly used in the reinitialization of the level set function. The purpose is to ensure that the level set function is a symbol distance function at least near the boundary, that is,∇ϕ=1. Both methods have a large computational burden. The approximate reinitialization method chosen in this paper controlled the value of the level set function near the boundary to be within a stable range [16]. For the entire design domain, reinitialization can be expressed as(31)ϕu=ϕmean∇ϕ10,∇ϕ20,…,∇ϕr0,where ϕr0 is the level set function value of the rth point near the level set zero contour and ∇ϕr0 is the second order norm of the gradient of ϕr0. mean∇ϕ10,∇ϕ20,…,∇ϕr0 is the average value of ∇ϕ10,∇ϕ20,…,∇ϕr0. Reinitialization is shown in Figure 4.Figure 4
Level set function reinitialization.There is a linear relationship between the level set functionϕ and the interpolation coefficient α, so after the level set function has been reinitialized, the interpolation coefficient appears as(32)αu=αmean∇ϕ10,∇ϕ20,…,∇ϕr0.Compared with the common reinitialization method, the method used in this paper was simple and would not hinder the ability of autonomous opening and can fully retain the boundary information and only change the relative size of the value.


## 4.4. Numerical Solution

After the level set function was interpolated with the MQ radial basis function, the numerical calculation was performed. Substituting formula (30) into the level set equation (20), the “level set equation” for the radial basis function interpolation coefficient can be expressed as(33)gxdαtdt−Vn∇gxαt=0,where(34)∇gxαt=∂gx∂xαt2+∂gx∂yαt2.For this time-varying interpolation problem, in order to ensure that the interpolation coefficients can be solved without being limited to the positive definition of the condition of the MQ function, a boundary constraint needs to be introduced.(35)∑i=1nα˙it=0,∑i=1nα˙itxi=0,∑i=1nα˙ityi=0.So far, the time-dependent H-J PDEs were discrete into a set of coupled ordinary differential equations. Equation (33) could also be written as(36)G∂α∂t−Bα,t=0,where(37)Bα,t=Vnx1,t∇gx1αt⋮Vnxn,t∇gxnαt000.The set of coupled nonlinear ordinary differential equations can be numerically solved by the first-order forward Euler method, and the approximate solution can be expressed as(38)αti+1=αti+ΔtG−1Bαti,ti,where Δt is the size of time step and ti is the ith time step.At the same time, theδ function is used to control the fast boundary growth rate during the solution process, which can be expressed as(39)δϕ=0,ϕ>Δ,τΔ1−ϕ2Δ2,ϕ≤Δ,where Δ is a parameter used to control the upper and lower limits of the level set function value. When the level set function value is higher than Δ or lower than −Δ, its normal velocity will be ignored, so as to avoid the level set function value tending to infinity. Generally, Δ is selected to be more than 5 times the size of the mesh. τ is usually set to 0.75 [16] and also can be adjusted to avoid instability during the optimization (see discussions in Section 5.1).Combined with the reinitialization method, the iterative equation (38) of the level set function can be expressed as(40)αti+1=αuti+ΔtG−1B^αuti,ti.Among them,δϕ can replace ∇ϕ, then there was(41)B^αuti,ti=Vnx1,tiδϕx1,αuti⋮Vnxn,tiδϕxn,αuti000.In addition, since the time step and the number of interpolation nodes of the radial basis function are determined by the requirements of the optimization convergence rather than the need to solve the H-J equation, the time and space discretization requirements of the CFL condition can be relaxed. Using equation (40), the interpolation coefficient is updated until the convergence condition is satisfied, then the level set function is interpolated, and the current optimized structure topology can be obtained by taking zero contour.


## 4.5. Lagrange Multiplier

In the calculation of equation (40), the iteration of the Lagrange multiplier λ can be expressed as(42)λk+1=μGk,k≤nR,λk+γkG,k>nR,where μ and γ are the parameters controlling the Lagrange multiplier, k represents the kth iteration, and nR is the minimum number of iterations to be optimized, and the iteration of γ is expressed as(43)γk+1=minγk+Δγ,γmax,k>nR.When the volume fraction of the initial design is not equal to the volume fraction of the constraint, the volume constraint can be relaxed during the initial iteration, expressed as follows:(44)Gk=∫ΩHϕdΩ−V0−V0−αΩknR,k>nR,where V0 is the initial volume.


## 4.6. Flow Chart of the Algorithm

The flow chart of the damping layer parameterization level set optimization algorithm based on the modal loss factor is shown in Figure5 and is described as follows:(1)
Given the initial design of the structure and meshing, initialize the level set function and radial basis function.(2)
Run a finite element analysis to calculate the modal loss factor of the current structure.(3)
Determine whether the number of iterations exceeds the specified minimum number of iterations, determine whether the current volume fraction reaches the constraint condition, and determine whether the difference between the current modal loss factor and the modal loss factor of the first five calculations is sufficiently small. When all three are satisfied, the loop ends, otherwise, the following steps are continued.(4)
Update the Lagrange multiplier according to the current volume fraction and related parameters, and calculate the evolution velocity field and update the level set function.(5)
Update the radial basis function interpolation coefficients and return to step (2).Figure 5
Flow chart of the parameterized level set method for optimizing the modal loss factor.


## 5. Numerical Analysis

The dimensions of the rectangle base plate is0.4m×0.3m, h = 0.001 m, E = 210 GPa, Poisson’s ratio v = 0.3, and ρ=7900Kg/m3. The 0.003 m thick viscoelastic layer has the constant properties of Ev = 100 MPa, Poisson’s ratio vv = 0.5, ρ1=1400Kg/m3, and ηv=0.25. The boundary conditions are taken as clamped at all four edges, and 40×30 elements are used in the FE model. The volume fraction is 0.5 as the constraint condition.


### 5.1. Single Modal

The evolution of the damping layer for the 5th modal with four holes as the initial design is shown in Figure6. We can see that the proposed method has the capability in nucleation of new holes inside the material domain. The corresponding level set surfaces are shown in Figure 7. Figure 8 shows the development of the objective function and volume constrain along the optimization process which converge smoothly and rapidly. Here, τ=0.6 in equation (39).Figure 6
Evolution history of the optimized design of the damping layer for the 5th modal with four holes as the initial design. (a) Step 1. (b) Step 15. (c) Step 25. (d) Step 81.
(a)
(b)
(c)
(d)Figure 7
Evolution history of the level set surface of the damping layer for the 5th modal with four holes as the initial design. (a) Step 1. (b) Step 15. (c) Step 25. (d) Step 81.
(a)
(b)
(c)
(d)Figure 8
Convergence curves of the modal loss factor (a) and volume fraction (b) of the damping layer for the 5th modal with four holes as the initial design (τ=0.6).
(a)
(b)Figure9 shows the evolution of the damping layer for the 5th modal with only one hole as the initial design. The optimal layout as shown in Figure 9(d) is quite the same as that shown in Figure 6(d). It can be concluded that the final design is fairly insensitive to initial designs. The corresponding level set surfaces are shown in Figure 10.Figure 9
Evolution history of the optimized design of the damping layer for the 5th modal with one hole as the initial design. (a) Step 1. (b) Step 14. (c) Step 24. (d) Step 76.
(a)
(b)
(c)
(d)Figure 10
Evolution history of the level set surface of the damping layer for the 5th modal with one hole as the initial design. (a) Step 1. (b) Step 12. (c) Step 25. (d) Step 75.
(a)
(b)
(c)
(d)The objective function and volume constrain have certain oscillations from step 10 to step 30 during the optimization procedure as shown in Figure11 when τ=0.6. Figure 12 shows the convergence curves when τ=0.8. The optimization result is almost consistent with Figure 9(d). Figure 13 shows the convergence curves when τ=0.5. The optimization result is shown in Figure 14. Comparing with these figures, we can see that higher τ leads to fast boundary growth, and a global optimal solution would be obtained when the optimization is convergence, while it also leads to numerical instability during the optimization procedure, or even nonconvergence, when τ=0.9 in this case. Lower τ leads to slow boundary growth, and the optimization procedure becomes more stable, while the optimization result may not the global optimal solution, see Figure 14. Thus, a moderate value of τ can control the boundary growth rate to guarantee accuracy of optimization results and stability of optimization process.Figure 11
Convergence curves of the modal loss factor (a) and volume fraction (b) of the damping layer for the 5th modal with one hole as the initial design (τ=0.6).
(a)
(b)Figure 12
Convergence curves of the modal loss factor (a) and volume fraction (b) of the damping layer for the 5th modal with one hole as the initial design (τ=0.8).
(a)
(b)Figure 13
Convergence curves of the modal loss factor (a) and volume fraction (b) of the damping layer for the 5th modal with one hole as the initial design (τ=0.5).
(a)
(b)Figure 14
Optimized design of the damping layer for the 5th modal with one hole as the initial design (τ=0.5).Figure15 shows the optimization results for the first 6 modals, respectively. The results are in close agreement with the modal strain energy contours of the base plate alone as shown in Figure 16. It shows that the conventional CAE methodology for location and size optimization based on the strain energy contours is a cost-effective approach, and the proposed level set method provides a more precise way to optimize surface damping treatments of panels at a single structural modal.Figure 15
Optimized results for the first 6 modals: (a) 1st modal, (b) 2nd modal, (c) 3rd modal, (d) 4th modal, (e) 5th modal, and (f) 6th modal.
(a)
(b)
(c)
(d)
(e)
(f)Figure 16
Strain energy contours for the first 6 modals: (a) 1st modal, (b) 2nd modal, (c) 3rd modal, (d) 4th modal, (e) 5th modal, and (f) 6th modal.
(a)
(b)
(c)
(d)
(e)
(f)


### 5.2. Multimodals

For multimodal application using the proposed level set method, the 1st and 6th modals are under consideration in this section. Here,a1 = a6 = 0.5. Figure 17(a) shows the final optimized design, and Figure 17(b) is the corresponding level set surface.Figure 17
(a) Optimized result for the 1st and 6th modals, (b) corresponding level set surface of Figure11(a), (c) MSEtotal contour for the 1st and 6th modals, and (d) MSEtotalnormalized contour for the 1st and 6th modals.
(a)
(b)
(c)
(d)For the conventional CAE method [1], the optimum damping treatment should be determined by the composite strain energy contours, which can be derived for each element as(45)MSEr,j=φr,jTKjeφr,j.For multimodal,(46)MSEtotal=∑r∑j=1nφr,jTKjeφr,j=∑rφrTKRφr.Figure17(c) shows the composite strain energy contour of the 1st and 6th modals. It is shown that the contour is similar to that of the 6th mode, as shown in Figure 16(f), and the 1st mode has little or no effect on the total strain energy, which fails to match the optimized layout.In the conventional FE modal, orthonormal modes are often used. That is,(47)φrTMφr=1.Then,(48)φrTKRφr=ωr2,which means the total composite strain energy of each modal equals ωr2, and then low order modes have less effect to the total composite strain energy when it was calculated using equation (46) for multimodal case, which does not match with the actual generally. To address this issue, here, we proposed an improved CAE methodology for location and size optimization of damping treatment, which is simple for industrial application. The normalized MSE can be expressed as(49)MSEtotalnormalized=∑rMSErnormalized=∑r∑j=1narφr,jTKjeφr,jφrTKRφr,which means the total MSE in each modal is normalized to 1. Figure 17(d) shows the normalized MSE contour of the 1st and 6th modals, which has almost the same distribution as the level set surface shown in Figure 17(b). Then, the locations of the damping layer should be determined at the place where normalized MSE is high, which would be similar with Figure 17(a).


## 5.1. Single Modal

The evolution of the damping layer for the 5th modal with four holes as the initial design is shown in Figure6. We can see that the proposed method has the capability in nucleation of new holes inside the material domain. The corresponding level set surfaces are shown in Figure 7. Figure 8 shows the development of the objective function and volume constrain along the optimization process which converge smoothly and rapidly. Here, τ=0.6 in equation (39).Figure 6
Evolution history of the optimized design of the damping layer for the 5th modal with four holes as the initial design. (a) Step 1. (b) Step 15. (c) Step 25. (d) Step 81.
(a)
(b)
(c)
(d)Figure 7
Evolution history of the level set surface of the damping layer for the 5th modal with four holes as the initial design. (a) Step 1. (b) Step 15. (c) Step 25. (d) Step 81.
(a)
(b)
(c)
(d)Figure 8
Convergence curves of the modal loss factor (a) and volume fraction (b) of the damping layer for the 5th modal with four holes as the initial design (τ=0.6).
(a)
(b)Figure9 shows the evolution of the damping layer for the 5th modal with only one hole as the initial design. The optimal layout as shown in Figure 9(d) is quite the same as that shown in Figure 6(d). It can be concluded that the final design is fairly insensitive to initial designs. The corresponding level set surfaces are shown in Figure 10.Figure 9
Evolution history of the optimized design of the damping layer for the 5th modal with one hole as the initial design. (a) Step 1. (b) Step 14. (c) Step 24. (d) Step 76.
(a)
(b)
(c)
(d)Figure 10
Evolution history of the level set surface of the damping layer for the 5th modal with one hole as the initial design. (a) Step 1. (b) Step 12. (c) Step 25. (d) Step 75.
(a)
(b)
(c)
(d)The objective function and volume constrain have certain oscillations from step 10 to step 30 during the optimization procedure as shown in Figure11 when τ=0.6. Figure 12 shows the convergence curves when τ=0.8. The optimization result is almost consistent with Figure 9(d). Figure 13 shows the convergence curves when τ=0.5. The optimization result is shown in Figure 14. Comparing with these figures, we can see that higher τ leads to fast boundary growth, and a global optimal solution would be obtained when the optimization is convergence, while it also leads to numerical instability during the optimization procedure, or even nonconvergence, when τ=0.9 in this case. Lower τ leads to slow boundary growth, and the optimization procedure becomes more stable, while the optimization result may not the global optimal solution, see Figure 14. Thus, a moderate value of τ can control the boundary growth rate to guarantee accuracy of optimization results and stability of optimization process.Figure 11
Convergence curves of the modal loss factor (a) and volume fraction (b) of the damping layer for the 5th modal with one hole as the initial design (τ=0.6).
(a)
(b)Figure 12
Convergence curves of the modal loss factor (a) and volume fraction (b) of the damping layer for the 5th modal with one hole as the initial design (τ=0.8).
(a)
(b)Figure 13
Convergence curves of the modal loss factor (a) and volume fraction (b) of the damping layer for the 5th modal with one hole as the initial design (τ=0.5).
(a)
(b)Figure 14
Optimized design of the damping layer for the 5th modal with one hole as the initial design (τ=0.5).Figure15 shows the optimization results for the first 6 modals, respectively. The results are in close agreement with the modal strain energy contours of the base plate alone as shown in Figure 16. It shows that the conventional CAE methodology for location and size optimization based on the strain energy contours is a cost-effective approach, and the proposed level set method provides a more precise way to optimize surface damping treatments of panels at a single structural modal.Figure 15
Optimized results for the first 6 modals: (a) 1st modal, (b) 2nd modal, (c) 3rd modal, (d) 4th modal, (e) 5th modal, and (f) 6th modal.
(a)
(b)
(c)
(d)
(e)
(f)Figure 16
Strain energy contours for the first 6 modals: (a) 1st modal, (b) 2nd modal, (c) 3rd modal, (d) 4th modal, (e) 5th modal, and (f) 6th modal.
(a)
(b)
(c)
(d)
(e)
(f)


## 5.2. Multimodals

For multimodal application using the proposed level set method, the 1st and 6th modals are under consideration in this section. Here,a1 = a6 = 0.5. Figure 17(a) shows the final optimized design, and Figure 17(b) is the corresponding level set surface.Figure 17
(a) Optimized result for the 1st and 6th modals, (b) corresponding level set surface of Figure11(a), (c) MSEtotal contour for the 1st and 6th modals, and (d) MSEtotalnormalized contour for the 1st and 6th modals.
(a)
(b)
(c)
(d)For the conventional CAE method [1], the optimum damping treatment should be determined by the composite strain energy contours, which can be derived for each element as(45)MSEr,j=φr,jTKjeφr,j.For multimodal,(46)MSEtotal=∑r∑j=1nφr,jTKjeφr,j=∑rφrTKRφr.Figure17(c) shows the composite strain energy contour of the 1st and 6th modals. It is shown that the contour is similar to that of the 6th mode, as shown in Figure 16(f), and the 1st mode has little or no effect on the total strain energy, which fails to match the optimized layout.In the conventional FE modal, orthonormal modes are often used. That is,(47)φrTMφr=1.Then,(48)φrTKRφr=ωr2,which means the total composite strain energy of each modal equals ωr2, and then low order modes have less effect to the total composite strain energy when it was calculated using equation (46) for multimodal case, which does not match with the actual generally. To address this issue, here, we proposed an improved CAE methodology for location and size optimization of damping treatment, which is simple for industrial application. The normalized MSE can be expressed as(49)MSEtotalnormalized=∑rMSErnormalized=∑r∑j=1narφr,jTKjeφr,jφrTKRφr,which means the total MSE in each modal is normalized to 1. Figure 17(d) shows the normalized MSE contour of the 1st and 6th modals, which has almost the same distribution as the level set surface shown in Figure 17(b). Then, the locations of the damping layer should be determined at the place where normalized MSE is high, which would be similar with Figure 17(a).


## 6. Experimental Validation

The images of the plates with different damping distributions under test are shown in Figure18. Shape (a), shape (b), and shape (c) are the optimized results for the 1st modal shown in Figure 15(a), the 6th modal shown in Figure 15(f), and multimodal for the 1st and 6th modals shown in Figure 17(a), respectively. Shape (d) is fully covered with damping materials. The experimental setup is shown in Figure 19. The impact modal test is carried out. Exciting force is applied by using a force hammer, and response signal is tested by using an acceleration sensor located at point (0.12 m, 0.08 m). The modal parameters are analyzed by using the PolyMAX Modal Analysis module of LMS test.lab.Figure 18
The images of the plates with different damping distributions.
(a)
(b)
(c)
(d)Figure 19
Experimental setup.The modal frequency and modal damping under test are shown in Table1. We can see that the resonances shift to high frequencies when damping treatments are applied to the plate. Treatments with shape (a) and shape (b) achieve max damping ratio at the 1st modal and 6th modal, respectively, for 50% coverage, and treatment with shape (c) gets moderate damping at each modal but total max damping ratio is achieved. The experimental results are in good agreement with the optimization results. The differences between full-covered case and optimized result for each modal are very small while 50% damping material is saved.Table 1
Modal frequency and damping ratio of test plates.
1st modal
6th modal
Frequency (Hz)
Damping (%)
Frequency (Hz)
Damping (%)
Bared
71.75
0.17
325.31
0.68
Shape (a)
86.87
3.81
364.23
1.28
Shape (b)
77.45
2.90
359.36
1.5
Shape (c)
83.47
3.77
362.83
1.39
Shape (d)
87.02
3.84
372.79
1.58Frequency response functions between the exciting point (0.2 m, 0.14 m) and response point (0.12 m, 0.08 m) are shown in Figure20. The parameters of the composite plates of FEM model are almost the same with that described in Section 5, except for the material damping ratio of aluminium, and damping distributions is set to 0.01 and 0.3 to simulate the actual composite, respectively. The analysis of numerical simulation and experiment has similar results and provides additional evidence to the conclusion of the modal test result from Table 1. The reason for difference between numerical simulation and experiment is that the damping ratio and elastic modulus of damping treatments actually vary with temperature and frequency, while in the FEM model we suppose it is constant. The optimized results are not sensitive with the variation of these parameters in a certain range in fact.Figure 20
Frequency response functions between the exciting point (0.2 m, 0.14 m) and response point (0.12 m, 0.08 m): (a) Finite element method and (b) experimental result.
(a)
(b)


## 7. Conclusions

Topology optimization of free damping treatments on plates has been presented in the paper using the parameterized level set method. The objective function is defined as a combination of several desired modal loss factors solved by the FE-MSE method. A clamped rectangle plate has been applied to demonstrate the validation of the proposed approach. The optimized results for a single modal show that the proposed method has the possibility of nucleation of new holes inside the material domain, and the final design is insensitive to initial designs. The damping treatments are guided towards the areas with high modal strain energy. For multimodal case, the optimized result matches the normalized modal strain energy of the base plate, which would provide a simple implementation way for industrial application. Experimental results show good agreements with the proposed method. The experimental results show good agreement with the optimization results. It is very promising to see that the optimized result for each modal has almost the same damping effect as that of the full-coverage case, and the result for multimodal gets moderate damping at each modal for suppressing the plate vibration with 50% of weight saved.



---

*Source: 5084167-2020-01-23.xml*