Noise Pollution Reduction through a Novel Optimization Procedure in Passive Control Methods

Haojie Lian,Leilei Chen,Xiao Lin,Wenchang Zhao,Stephane P.A.Bordas and Mingdong Zhou

1Key Laboratory of In-Situ Property-Improving Mining of Ministry of Education,Taiyuan University of Technology,Taiyuan,030024,China

2Henan International Joint Laboratory of Structural Mechanics and Computational Simulation,Huanghuai University,Zhumadian,463000,China

3College of Architecture and Civil Engineering,Xinyang Normal University,Xinyang,464000,China

4The York Management School,University of York,York,YO10 5DD,UK

5CAS Key Laboratory of Mechanical Behavior and Design of Materials, University of Science and Technology of China,Hefei,230027,China

6Institute for Computational Engineering, Faculty of Science, Technology and Communication, University of Luxembourg,Luxembourg,SA2 8PP,UK

7School of Engineering,CardiffUniversity,Cardiff,CF24 3AA,UK

8Shanghai Key Laboratory of Digital Manufacture for Thin-Walled Structures,School of Mechanical Engineering,Shanghai Jiao Tong University,Shanghai,200240,China

ABSTRACT This paper proposes a novel optimization framework in passive control techniques to reduce noise pollution.The geometries of the structures are represented by Catmull-Clark subdivision surfaces,which are able to build gap-free Computer-Aided Design models and meanwhile tackle the extraordinary points that are commonly encountered in geometric modelling.The acoustic fields are simulated using the isogeometric boundary element method,and a density-based topology optimization is conducted to optimize distribution of sound-absorbing materials adhered to structural surfaces.The approach enables one to perform acoustic optimization from Computer-Aided Design models directly without needing meshing and volume parameterization,thereby avoiding the geometric errors and time-consuming preprocessing steps in conventional simulation and optimization methods.The effectiveness of the present method is demonstrated by three dimensional numerical examples.

KEYWORDS Noise control; topology optimization; Catmull-Clark subdivision surfaces; isogeometric analysis; boundary element methods

1 Introduction

Installing sound-absorbing materials on structures is considered as an effective passive control technique to reduce noise level [1-3].Considering the increase in weight and cost, the full coverage of absorption materials is impractical.Under the volume constraint, the materials should be distributed reasonably to attain a satisfactory soundproofing effect.As a pace-setting technique,topology optimization serves as an effective tool to solve this engineering problem [4-6].Topology optimization is a mathematical method that optimizes material distribution for particular objective functions and constraints.In topology optimization, the design is improved iteratively until it converges to an optimal solution.The first application of topology optimization in acoustics is traced back to the research of Dhring et al.[6], in which outdoor sound barriers were optimized and remarkable noise reduction was achieved.

As a versatile numerical technique, the finite element method (FEM)is widely used in topology optimization.For acoustic problems, the boundary element method (BEM)is also commonly used for its advantages in dealing with unbounded domain problems [7-12].However, FEM or BEM relies on a preprocessing step that converts the structural geometries into polygonal meshes, which is time-consuming and results in geometric errors.The isogeometric boundary element method (IGABEM)has recently emerged as a competitive alternative to the conventional numerical methods.The core concept of the IGABEM is to employ basis functions used in Computer-Aided Design (CAD)for geometric modeling to solve boundary integral equation that are transformed from the partial differential equations.As such, IGABEM enables numerical simulation to be conducted from CAD directly which avoids the cumbersome meshing procedure and retains geometric accuracy [13-15].Unlike isogeometric finite element analysis, a volume parameterization is not needed in IGABEM, because both BEM and CAD are boundaryrepresented.IGABEM has been successfully applied in linear elasticity [16,17], acoustic analysis[18,19], structural optimization [20-23], etc.

Subdivision surface modeling is an important CAD technique to represent the complex surface geometries of structures [24-26].Compared to Non-Uniform Rational B-splines (NURBS),the main merit of subdivision surfaces lies in its ability of constructing gap-free models which are a prerequisite for the numerical analysis.Apart from that, subdivision surfaces are able to tackle extraordinary points at which the curvatures are not bounded.Although many different schemes have been proposed in the category of subdivision surface modeling since its inception [27-30], its first version, the Catmull-Clark subdivision surface [31-33], is still gaining popularity in practice and is being applied to a large variety of engineering problems [34-37].This may be attributed to its simplicity, generality, and well-established ecosystem.Although the Catmull-Clark subdivision surfaces have been successfully incorporated into IGABEM for numerical simulation [15,38], to the best of our knowledge, no study has been reported on topology optimization thus far.

This study aims to fill this research gap.We propose a topology optimization procedure to design the distribution of sound-absorbing materials, in which the Catmull-Clark subdivision surface is adopted to construct geometric models and its basis functions are used to simulate acoustic fields.The objective is to minimize the noise level subject to the volume constraints of absorption materials.As an extension of our previous work [39] where the Loop subdivision surface was used in IGABEM for topology optimization, the incorporation of Catmull-Clark subdivision surfaces into the acoustic topology optimization procedure enables us to use quadrilateral meshes and thus enhances our capability of geometric modeling.The main strength of the present method arises from its capability of seamlessly integrating numerical analysis and CAD models for complex geometries, which enhances the efficiency and accuracy of topology optimization.

The remainder of this paper is organized as follows.In Section 2, the fundamentals of Catmull-Clark subdivision surface modelling are briefly reviewed for completeness.Section 3 introduces the IGABEM formulation in acoustic analysis that takes into account impedance boundary conditions related to the acoustical effects of absorbing materials.Section 4 formulates the density-based acoustic topology optimization method in which the interpolation scheme of absorption material impedance and the sensitivity analysis with adjoint variable method are highlighted.Section 5 provides numerical examples to verify the correctness and efficiency of the proposed approach, followed by the conclusions in Section 6.

2 Catmull-Clark Subdivision Surfaces

The first step for a designer to build geometries with Catmull-Clark subdivision surfaces is to construct a control grid, which is a polygon mesh with quadrilateral elements whose vertices are called control points.Subsequently, the control grid is subdivided, new control points are introduced, and the positions of the existing control points are updated in Fig.1.The control grid can be subdivided repeatedly, and the subdivision process from the subdivision levelkto the levelk+1 follows the following rule in Fig.2.

Figure 1:Catmull-Clark subdivision algorithm for surfaces.The black points indicate the control points before the subdivision, and the red points indicate the control points after the subdivision.(a)Initial control mesh (b)One-level refined mesh

Figure 2:Topology scheme of a Catmull-Clark subdivision surface.The green circle represents the E-vertex inserted at the middle of an edge, the purple circle is the F-vertex inserted at the middle of an element, and the red circle represents the V-vertex formed by the transfer of a vertex in a parent element

• E-vertex:Let the two vertices of the inner edge bexk1andxk4, and the vertices of the two faces sharing this edge bexk2,xk3,xk5andxk6, respectively.Subsequently, the new edge pointxk+1generated by this inner edge is

• F-vertex:If the vertices on each face arexk1,xk2,xk3andxk4, then the new face pointxk+1generated by the face is

• V-vertex:For an internal vertexxk, we denote its one-neighbor vertices byxki(i=1,2,...,2v), in which the vertices with odd subscripts are directly connected toxk, and the ones with even subscripts are on the diagonal line of the quadrilateral element.Correspondingly, the updated position of vertexxk+1is

After the control points of the subdivided control grid are determined, a smooth surface can be constructed through a linear combination of control points and basis functions:

whereniis the number of basis functions,(ξ,η)∈[0,1]2is the parametric element, andB(ξ,η)denotes the basis functions that are expressed by

where the symbols % and/represent the remainder and division, respectively.The functionNi(ξ)represents the cubic uniform B-spline basis functions, i.e.,N1(ξ)=(1 −3ξ+3ξ2+ξ3)/6,N2(ξ)=(4 −6ξ2+3ξ3)/6,N3(ξ)=(1+3ξ+3ξ2−3ξ3)/6,N4(ξ)=(ξ3)/6.

Fig.3 shows an example of a werewolf generated by Catmull-Clark subdivision surfaces.The second row of the figure represents the head, hands, and feet of the werewolf model after it is subdivided once.The example shows that the subdivision surface method is highly adaptable to complex geometry and capable of constructing smooth surfaces.It is noted that the the final smooth geometries are the same regardless of the subdivision times of the control grid.

Figure 3:A werewolf model generated from Catmull-Clark subdivision surfaces.The number of elements in the initial model was 3444, which increases to 13776 after subdivision once

3 Acoustic Simulation Using IGABEM

Numerical simulation using the IGABEM can be used to compute the objective function and evaluate the design performance at each iterative step in the topology optimization process.The acoustic field is governed by the Helmholtz equation in domainΩ∈R3:

where ∇2is the Laplace operator,p(x)is the sound pressure,k=ω/cis the wave number withωbeing the circular frequency andcthe wave velocity.By transforming and infinitely approaching point x in the domain towards the boundary and using Green’s second theorem, the boundary integral equation at point x is obtained

where x and y denote the source point and field point, respectively,n(y)is the outer normal vector at the source point y,Γdenotes the structural surface, andq(y)=∂p(y)/∂n(y)is acoustic flux.By denotingpincas the incident acoustic pressure andpscas scattered acoustic pressure, the total acoustic pressurep(x)is written asp(x)=pinc+psc.The Green’s functionG(x,y)and its normalized derivative can be expressed as

where i is the imaginary unit andr=|x −y|represents the distance between the source and field points.G(x,y)andare also called kernel functions or fundamental solutions.To characterize the sound-absorption properties of adhesive sound-absorption materials on the structural surface, the impedance boundary condition is introduced

whereβ(y)denotes the normalized acoustic admittance at field point y.

The conventional boundary condition may yield fictitious eigen-frequencies.This problem can be resolved by Burton-Miller method [40-45], which is a linear combination of conventional boundary integral equation (Eq.(7))and its normal derivative.Taking into account the impedance boundary condition, the Burton-Miller method can be formulated as

whereαis the coupling parameter which is equal to i/kfork>1 and i otherwise.The above equations hold for both exterior and interior domain problems with the only difference in the direction of surface normal vector.

For a structural surface that is modeled using Catmull-Clark subdivision surfaces, whose basis functionsBiare used to discretize the acoustic field around the surface,

wherepeandqedenote the acoustic pressure and its normal flux at a field point(ξ,η)in a subdivision surface element, and peiand qeiare the control point parameters for discretizing the sound pressure and normal flux on the boundary elementΓe, respectively.

To transform the boundary integral equation to linear algebraic equations, we placed the source points at a set of discrete collocation points on the boundary and enforce the governing equations to be satisfied on them.In conventional BEM, the nodes of the mesh are normally chosen as collocation points.In the present work, the collocation points are selected as the points on the smooth surface that are mapped from the nodes of the control grid.With the collocation scheme, the system equation can be rewritten as

whereNeis the number of elements.The above equations can be collected as

where H and G are the coefficient matrices of the IGABEM with the Catmull-Clark subdivision scheme, p is a column vector collecting sound pressures at the collocation points, pinccollects the incident acoustic pressure at the collocation points, and C is the admittance matrix which is written as

Because the kernel functions are singular at the source points, the singular integral arises in Eq.(12)which has to be addressed carefully [46,47].To overcome this problem, the singularity subtraction technique (SST)is adopted in this work.SST subtracts the singular part from the integrand and integrate it analytically.The remaining part of the integrand is regular which can be integrated numerically with Gaussian quadrature.

4 Topology Optimization

4.1 Optimization Model

We consider an optimization problem for the absorbing-material distribution to minimize sound pressures subject to material volume constraints.The topology optimization problem can be formulated as follows:

where the objective functionΠis the quadratic sum of sound-pressure amplitudes at the reference points, pfcollects the sound pressure at one or several inspection points,(·)Hdenotes the conjugate transpose of the vector,ρeandveare the density and volume of elementΓe, respectively,andfvis the constraint of the volume fraction.The value of pfis obtained using the discretized integral equation:

where matrices Hfand Gfas well as the vectorare defined in a similar way to those in Eq.(13)except that the source point is outside the domain.

In the present work, the density-based method is used for topology optimization.As a gradient-based optimization algorithm, it is more mathematically sound and converges to the optimized solution more rapidly than the gradient-less optimization methods like genetic algorithms[48].To tackle a large number of design variables, the adjoint variable method is adopted to compute the sensitivities that play a vital role in density-based topology optimization procedure.In viewing of (13)and (16), we rewrite the objective function by adding zero functions as

whereλ1andλ2are vectors of adjoint variables whose values can be arbitrarily selected.By differentiating Eq.(17)with respect toρe, we obtain the following equation:

After rearrangements, Eq.(18)can be reformulated as

where R represents the operation of determining the real part of the complex number.Finally,the derivative of the objective function with respect to design variableρeis expressed by

As mentioned above, the adjoint vectorsλT1andλT2can be arbitrarily selected if they satisfy the following equation:

Eq.(21)is called adjoint equation.After solving it, the adjoint variablesλ1andλ2can be obtained.Upon substitution ofλ1andλ2into Eq.(20), the sensitivities of the objective function can be evaluated.

4.2 Interpolation Scheme of Acoustic Absorbing Materials

According to the Delany-Bazley-Miki model [49], the normalized impedance of sound absorption materials reads

whereσis the flow resistance rate of the material(N·s/m4),fis the frequency (Hz).Correspondingly the normalized admittance value is expressed as

To make the material distribution close to a 0-1 design, a material interpolation model, the so-called solid isotropic material with penalization (SIMP)method, is used to interpolate the element admittance:

whereηis the penalty parameter that has a function in converging the intermediate density to 0 or 1 and generally assumes the value of 3.The larger the penalty parameter, the closer the design variable value is to 0-1, but the increase in the value ofηconverges the optimization result to the local minimum [4].

After the sensitivities are evaluated, a gradient-based optimizer can be used to update the design variables until it converges.In the present work, the method of moving asymptotes is used as the optimizer and the convergence criterion is set as

whereΠistands for the objective function value at thei-th iteration.Finally, the Heaviside function can be used as an additional filter to remove intermediate densities to achieve a 0-1 design.

5 Results and Discussion

A numerical example is presented in this section to investigate the effectiveness of the topology optimization approach using the IGABEM with Catmull-Clark subdivision surfaces.The iterative convergence criterion is set to 1.0×10−4.The geometries in the following examples have C1 continuity.

5.1 Werewolf Model

Now, we consider the model mentioned in Section 2, which is used to demonstrate the applicability of the proposed optimization method in addressing complex geometries.Herein, the plane wave spreads along thexdirection, and the coordinates of the test point used for calculation of objective function are (12, −5, 0).The optimized distribution of sound-absorbing materials at different frequencies is analyzed, and four different frequencies (f=50, 100, 150 and 200 Hz)are considered in Fig.4.The color red indicates that the area is covered with sound-absorbing materials, the color blue indicates that the area is rigid.From this figure, we can observe that the final optimized distribution is symmetric with respect to the XOY plane.When the excitation frequency is different, the optimized distribution obtained by calculation is different.This indicates that the optimized distribution of sound-absorbing material has a frequency dependence.

Figure 4:Optimized distribution of sound-absorbing materials of werewolf model at different frequencies

Fig.5 shows the real part, imaginary part, and amplitude of the sound pressure and its flux on the optimized structural surface at 100 Hz, respectively.A good symmetry along the XOY plane can be observed, and it verifies the correctness of the algorithm developed in this study.

Figure 5:Contour of the acoustic pressure and flux on the werewolf surface with an optimized distribution of sound-absorbing materials

Similarly, whenf=50, 100, 150 and 200 Hz, the distribution of the sound-pressure level after optimization is shown in Fig.6.We can observe that the higher the frequency, the more intense the change in the distribution of the sound-pressure level.In addition, the value of the sound-pressure level on many structural areas increases with the increase in frequency.

Figure 6:Contour of the acoustic pressure levels (SPL)on the werewolf surface with an optimized distribution of sound-absorbing materials

The effect of different parameters on the objective function and material volume constraints is investigated.First, the effect of penalty factor on optimization results is considered.Theoretically,the larger the value ofη, the closer the design variable value is to 0 and 1.In each iteration step,when the design variable is updated, the greater the change in the value of the objective function,and the more apparent the oscillation phenomenon of the error curve.In addition, the increase inηconverges the optimization result to the local minimum.When the value ofηdecreases, although these problems are avoided, it may be difficult to eliminate the intermediate variable elements,resulting in the optimized distribution containing many green elements.

Figs.7 and 8 show the objective function and volume fraction function in terms of the iteration step with different values ofη.The figures indicate that at the initial several iteration steps,the objective function values corresponding to different penalty factor values are considerably different.When the convergence approaches, the difference in objective function values decreases.Although the value of the penalty factor has a slight effect on the final objective function value,this does not imply that it has a slight effect on the optimized distribution result.On the contrary,the setting of the penalty factor causes a more apparent change in the optimized distribution map.Therefore, an appropriate penalty factor value must be selected in the topology optimization process.If the penalty factor is excessively small, it causes the intermediate density to accumulate,and if the penalty factor is excessively large, the penalty is excessive.A satisfactory result is to set the parameter value to 3 or 5.

Figure 7:Iteration history of the objective function at different values of η

Figure 8:Iteration history of the material volume fraction at different values of η

We consider the effects of the initial values of design variables on the optimization results.Figs.9 and 10 show the objective and material volume fraction functions in terms of iteration steps with different initial values, respectively.The volume constraint is still set as 0.5.We set the initial values as (0.2, 0.4, 0.6, 0.8, 1.0)for comparative analysis.Fig.9 shows that in the initial iteration step, the initial value of the design variables has a significant effect on the objective function value, but as the iteration steps increase, the effect decreases rapidly until it converges to the same level value.The results indicate that the initial value of design variables has a slight effect on the topology optimization of material distribution.

Figure 9:Iteration history of the objective function with different initial values

Figure 10:Iteration history of the material volume fraction with different initial values

5.2 Car Model

In this example, a car body shell model is constructed using Catmull-Clark subdivision surfaces with an initial control grid with 2288 elements and 2290 vertices in Fig.11.We test the proposed algorithm by solving the acoustic scattering problems and optimizing the layout of sound-absorbing materials attached to the car surfaces.The incident plane acoustic wave propagates forward along thexdirection, and the excitation frequency is set as 100 Hz.The optimization aim is to minimize the acoustic pressure at a sample Point A with coordinates(20, 5, 0).

Figure 11:Initial control grid of a car model (left); smooth Catmull-Clark surface of the car model (right)

Herein, the volume ratio constraint is 0.5, that is, at most 50% of the structure surface is covered by acoustic material.The initial value of the design variable is set as 1.The material volume fraction varies from 0 to 1, where 1 indicates that all the structural surfaces are covered by sound-absorbing materials and 0 represents rigid structural surfaces.

In the iteration process of the optimization, the objective and material volume fraction functions vary with the number of iterative steps, as shown in Fig.12.The objective function initially increases abruptly, decreases rapidly, and finally converges stably because of the full coverage of the initial design.In addition, the material volume ratio decreases rapidly from 1 to 0.35 in the initial several iteration steps, and then remains constant at approximately 0.5 after several steps.Note that although the change in volume ratio is significantly small after the 15-th iteration step, the corresponding material distribution may have an apparent change until convergence is achieved.

Figure 12:Iteration history of the car model after the first subdivision during the optimization process

As the number of iterative steps increases, the optimized distribution of sound-absorbing materials is shown in Fig.13.We can observe that the optimized distribution is symmetric around the XOY plane, and this is because the original optimization problem is symmetric around the XOY plane.Therefore, for such axisymmetric problems, the symmetric part of the structure can be selected as the optimized object, which can both guarantee the complete symmetry of the final topology design and further reduce the computation.The sound-pressure levels of the car model after the optimization of distribution of sound-absorbing materials is shown in Fig.14.The acoustic pressure level of the area covered by the acoustic absorbing material is significantly lower than that of the area not covered.Similarly, a noticeable symmetry can be observed.The results verify the correctness of the algorithm and the effectiveness of the optimization analysis of complex models.

Figure 13:Sound-absorbing material distribution at different iteration steps

Figure 14:Acoustic pressure (SPL)around the car surface after the distribution of soundabsorbing materials is optimized

6 Conclusions and Future Work

In this study, we propose a method to optimize the distribution of sound-absorbing materials in noise pass control techniques based on isogeometric boundary element methods.The complex geometries of components or structures are modeled using Catmull-Clark subdivision surfaces,whose basis functions are also employed for acoustic analysis.The proposed method eliminates geometric errors and cumbersome preprocessing steps in traditional topology optimization procedures.The numerical results indicate that sound pressure levels can be significantly reduced subject to the given constraint of the volume of absorbing materials.However, the present paper only considers topology optimization.Combining shape and topology optimization with subdivision surfaces will deliver a more flexible design [50,51].In addition, a main limitation of the present method is that the structural-acoustic interaction effect is not taken into account, which will be studied by coupling FEM and BEM in an isogeometric analysis framework in the future work.We will also extend the present work to geometric uncertainly qualification [52-54] because IGABEM enables rapid sampling in stochastic analysis by allowing for numerical simulation directly from CAD.

Funding Statement:We acknowledge the support of the National Natural Science Foundation of China (NSFC)under Grant Nos.51904202 and 11702238.Stephane Bordas thanks the financial support of Intuitive modeling and SIMulation platform (IntuiSIM)(PoC17/12253887)grant by Luxembourg National Research Fund.

Conflicts of Interest:The authors declare that they have no conflicts of interest to report regarding the present study.