|Computer Modeling in Engineering & Sciences|
Noise Pollution Reduction through a Novel Optimization Procedure in Passive Control Methods
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, Cardiff University, 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
*Corresponding Authors: Wenchang Zhao. Email: firstname.lastname@example.org; Mingdong Zhou. Email: email@example.com
Received: 10 October 2021; Accepted: 02 November 2021
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
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. , 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 boundary-represented. 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  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 level k to the level k + 1 follows the following rule in Fig. 2.
• E-vertex: Let the two vertices of the inner edge be and , and the vertices of the two faces sharing this edge be , , and , respectively. Subsequently, the new edge point xk+1 generated by this inner edge is
• F-vertex: If the vertices on each face are , , and , then the new face point xk+1 generated by the face is
• V-vertex: For an internal vertex xk, we denote its one-neighbor vertices by , in which the vertices with odd subscripts are directly connected to xk, and the ones with even subscripts are on the diagonal line of the quadrilateral element. Correspondingly, the updated position of vertex xk+1 is
where , , .
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:
where ni is the number of basis functions, is the parametric element, and denotes the basis functions that are expressed by
where the symbols and / represent the remainder and division, respectively. The function represents the cubic uniform B-spline basis functions, i.e., .
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.
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 :
where is the Laplace operator, is the sound pressure, is the wave number with being the circular frequency and c the wave velocity. By transforming and infinitely approaching point in the domain towards the boundary and using Green’s second theorem, the boundary integral equation at point is obtained
where and denote the source point and field point, respectively, is the outer normal vector at the source point , denotes the structural surface, and is acoustic flux. By denoting as the incident acoustic pressure and as scattered acoustic pressure, the total acoustic pressure is written as . The Green’s function and its normalized derivative can be expressed as
where is the imaginary unit and represents the distance between the source and field points. and are 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 denotes the normalized acoustic admittance at field point .
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 for k > 1 and 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 functions Bi are used to discretize the acoustic field around the surface,
where pe and qe denote the acoustic pressure and its normal flux at a field point in a subdivision surface element, and and are the control point parameters for discretizing the sound pressure and normal flux on the boundary element , 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
where Ne is the number of elements. The above equations can be collected as
where and are the coefficient matrices of the IGABEM with the Catmull-Clark subdivision scheme, is a column vector collecting sound pressures at the collocation points, collects the incident acoustic pressure at the collocation points, and 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, collects the sound pressure at one or several inspection points, denotes the conjugate transpose of the vector, and ve are the density and volume of element , respectively, and fv is the constraint of the volume fraction. The value of is obtained using the discretized integral equation:
where matrices and as well as the vector are 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 . 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 and are vectors of adjoint variables whose values can be arbitrarily selected. By differentiating Eq. (17) with respect to , we obtain the following equation:
After rearrangements, Eq. (18) can be reformulated as
where represents the operation of determining the real part of the complex number. Finally, the derivative of the objective function with respect to design variable is expressed by
As mentioned above, the adjoint vectors and can be arbitrarily selected if they satisfy the following equation:
Eq. (21) is called adjoint equation. After solving it, the adjoint variables and can be obtained. Upon substitution of and into 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 , the normalized impedance of sound absorption materials reads
where is the flow resistance rate of the material , f is 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 .
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 stands for the objective function value at the i-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 . 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 the x direction, 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.
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.
Similarly, when f = 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.
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.
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.
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 the x direction, 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).
Herein, the volume ratio constraint is 0.5, that is, at most 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.
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.
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.
|This work is licensed under a Creative Commons Attribution 4.0 International License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.|