A reduced-order cutFEM approach to model complex moving sound sources

Sjoerd van Ophem 1 and Wim Desmet KU Leuven, Department of Mechanical Engineering Celestijnenlaan 300B, B-3001, Heverlee, Belgium Flanders Make- DMMS corelab, Belgium

ABSTRACT While traditional finite element based vibro-acoustic analysis of structures is done on stationary structures, in many occasions the structure under investigation is actually moving. To perform accurate finite element analysis on the vibro-acoustic performance of such structures requires the evaluation of a time-varying system, leading to several challenges, such as the selection of a strategy to update the mesh to the new configuration and an adequate time integration scheme that does not introduce large numerical artifacts. Additionally, the time-varying nature leads to an increase in computational complexity. In this paper a modeling strategy based on cut finite elements (cutFEM) is introduced that aims to tackle these challenges. The cutFEM approach can be seen as a fictitious domain approach in which the source is moving over a static mesh, cutting through elements where the actual source is located, using a level set description. The key advantage of this approach is that no re-meshing and / or mesh morphing is required. It is shown how both implicit and explicit time integration schemes can be used to propagate the solution while the source is moving. Additionally, several model order reduction strategies are compared to speed up the online evaluation of the reduced order model.

1. INTRODUCTION

A common assumption that is done for numerical vibro-acoustic analysis, e.g. with the finite element method (FEM), is that the source position is stationary. However, in many cases the acoustic source is actually moving. For example, when vehicle pass-by noise is evaluated or in machines and vehicles with many moving parts. While a far-field moving sound source can often be modeled with a simplified source model, such as a moving point source, the accurate modeling of near-field sound propagation requires a more detailed source model. A moving sound source brings many challenges from a modeling perspective. The traditional approach with conformal finite elements would mean that at every time instance a re-meshing step is required. This is not only a costly operation, but also requires manual intervention to assure the mesh quality is su ffi cient. To circumvent or mitigate this problem, several alternative solutions have been proposed in the literature that can broadly be divided into mesh morphing and static mesh approaches.

1 sjoerd.vanophem@kuleuven.be

a slaty. inter.noise 21-24 AUGUST SCOTTISH EVENT CAMPUS O ¥, ? GLASGOW

A popular approach that is based on mesh morphing is to reformulate the problem to an Arbitrary Lagrangian-Eulerian (ALE) formulation [1]. The advantage of the ALE form is that the degrees of freedom (DOFs) do not change when the source is moving, but morph to a new location. However, when the elements become too distorted, which is possible when the source moves large distances as compared to the acoustic wavelengths of interest, a re-meshing step is required, which is highly undesirable. An alternative approach is to use a static mesh, such that the problem does not require a re-meshing step during simulation, and use techniques to enrich the mesh to include the proper source description. Several techniques to do this are known, such as the extended finite element method (XFEM) and the cut finite element method (cutFEM) [2–4]. The disadvantage of these methods with respect to the ALE formulation is that the DOFs will change during the simulation, since the DOFs that are covered by the source fall outside of the physical domain and will not contribute to the solution. In this paper, it is chosen to investigate a cutFEM approach, due to the convenience of not having to apply mesh morphing and / or re-meshing. In Sec. 2, A cutFEM approach is derived that can model moving acoustic sources, while minimizing numerical artifacts. Several time-integration schemes are compared to this end in Sec. 3, comparing both the time response directly and the frequency content of the signal. Besides accurate modeling, also a preliminary investigation is done on the creation of a reduced order model in Sec. 4.The computational cost of such simulation is especially high, since the system matrices have to be recalculated at each time step. A reduced order model would speed up calculation time and potentially enable (near-)real time vibro-acoustic virtual sensing, as has been shown in previous work for stationary sources [5,6].

2. CUTFEM FORMULATION

In the following section the cutFEM formulation will be elaborated. The main idea behind a cutFEM formulation is to make the mesh as independent as possible from the geometry of the problem under investigation, while keeping the accuracy and stability of conformal finite elements. The boundaries / interfaces of the geometry are described by cuts through elements in the static mesh and numerical integration of the cut elements is done until the cut. A general overview of the cutFEM formulation is given by Burman et al. [2].

2.1. Discretization of the wave equation with cutFEM To obtain the discretized acoustic cutFEM problem, we follow the approach given by Sticko and Kreiss [7]. Consider a 2-D domain in which the acoustic propagation is defined by the wave equation

c 2 0 ∇ 2 p − ¨ p = 0 . (1)

In this equation p represents the acoustic pressure, ¨ p is the second derivate of p with respect to time and c 0 the speed of sound in the fluid. Additionally, we assume that the source is described with a non-zero Neumann boundary condition (BC) at boundary Γ N

∂ p ∂ n | Γ N = g N . (2)

Within the scope of this paper, we do not consider any other type of BC or source description in this work. Now, suppose that a su ffi ciently refined background mesh T b exists of linear elements that covers a fictitious domain Ω T B that includes the physical domain Ω . In particular, we would like to solve our problem at the DOFs that correspond to the smallest mesh T that covers the physical domain Ω , see Fig. 1. Several elements of the fictitious domain will be intersected by the source boundary Γ , for which we assume that each intersected element, is intersected only once. Furthermore, the intersection is approximated as a straight line between the intersection points. By applying a test function v to Eq.

Figure 1: The background mesh T B is shown on the left. Assuming that the acoustic source is given by the circle, the problem is solved on mesh T that covers the physical domain Ω (in gray on the right) and the fictitious domain Ω T (the domain covered by all the elements on the right).

(1), performing integration by parts and the weak enforcement of the Neumann BCs, the following weak form results Z

Z

( ¨ p h v + c 2 0 ∇ p h . ∇ v ) d Ω= c 2 0

Γ N g N v d Γ , (3)

where p h indicates the finite element approximation of p .

2.2. Numerical integration and ghost penalty stabilization While the numerical integration of the non-cut elements is done the same as for the standard finite element method (using Galerkin’s method, and applying Gauss integration), numerical integration of the cut elements is done only for the part of the element that is inside the physical domain, as is described by Burman et al. [2]. To make integration easier, this area is subdivided into triangular areas, on which the required Gauss integration points will be defined, see Fig. 2. Note that, although the integration takes place until the cut (only integrating Ω ), the resulting element level mass and sti ff ness matrices will be fully populated, thus will give non-zero contributions for the non-physical DOFs as well (in case of Fig. 2, the DOFs corresponding to node 3 and 4). When only a small part of the cut element is inside the physical domain, this will lead to ill- conditioned element matrices, which can cause numerical artifacts and instability in the solution. Therefore, a stabilizing penalty term is added to the weak form that is applied across all of the interior faces F of the elements cut by the boundary, known as ghost penalty stabilization [8]. This stabilization term is incorporated in the weak form as follows

Z

( ¨ p h v + c 2 0 ∇ p h . ∇ v ) d Ω+ γ M h 2 Z

Z

Z

F h [ˆ n ·∇ p h ] . [ˆ n ·∇ v h ]dF = c 2 0

F h [ˆ n ·∇ ¨ p h ] . [ˆ n ·∇ v h ]dF + γ K

Γ N g N v d Γ ,

(4) in which [ˆ n · ∇ p h ] = ˆ n · p h | + − ˆ n · p h | − , is the jump in the normal derivative over each face, with normal ˆ n . Additionally, γ M and γ K are positive scaling factors that influence the level of stabilization. Thus, numerical integration of these penalty terms leads to cross-element stabilization terms to the mass and sti ff ness matrix. After numerical integration the weak form leads to a mass matrix M , a sti ff ness matrix K and a force vector F .

3. MODELING A MOVING SOURCE

In this paper we will only consider subsonic source movement, which is a reasonable assumption for most practical applications. The movement of the source implies that at every time instance the cut elements will di ff er, thus resulting in time-varying mass and sti ff ness matrices M ( t ) ∈ R n × n , and

Tp =~

3 4

1 2

Figure 2: Example of the placement of Gauss points (indicated with crosses) on a cut element. The gray area is split into triangles to perform numerical integration, summing up the individual contributions to the element level mass and sti ff ness matrices.

y v=10 m/s

x

Figure 3: An acoustic source (Neumann BC) is placed in the middle of the domain ( with sound hard walls) and moves in the positive x-direction during simulation.

K ( t ) ∈ R n × n , and a changing Neumann BC, F ( t ) ∈ R n × 1 . The determination of the locations of the cuts is done with the aid of a level set function ϕ , that is defined to be zero at the boundary, positive in the physical domain, and negative inside the source region. Also, the DOFs under investigation are time-dependent and updated with respect to the source location. For the DOFs that are moving outside of domain T the acoustic pressure is explicitly set to zero, while the DOFs that are moving back into domain T are initialized through the partially integrated cut elements.

3.1. Moving all the time-dependent matrices to the right-hand side Several time integration schemes were assessed to compare the performance. A straightforward way to write down the equations of motion is

M ( t ) ¨ p + K ( t ) p ( t ) = F ( t ) . (5)

However, the main issue with Eq. (5) is that the DOFs of interest will change when the source is moving, thus resulting in a change of dimensions during simulation. This can be circumvented by reformulating the time dependent problem to a form in which all of the time-dependent matrices are moved to the right-hand side, thus acting as time-dependent inputs. To this end, we first note that the sti ff ness matrix of a cut element K cut ∈ R n e × n e , where n e is the number of DOFs per element, can be written as follows

K cut ( t ) = K full − K cut , out ( t ) , (6)

in which K full ∈ R n e × n e is the sti ff ness matrix of the non-cut element (thus performing Gauss integration as usual) and K cut , out ∈ R n e × n e is the sti ff ness matrix of the element only considering the part of the element outside of the physical domain. If we look at the sti ff ness matrix at system level K ∈ R n × n , this means it can be rewritten as follows

K ( t ) = K b − K out ( t ) − K cut , out ( t ) + K GP ( t ) , (7)

where K b ∈ R n × n is the assembled sti ff ness matrix of the background mesh (all elements integrated fully with Gauss integration), K out ∈ R n × n is the assembled matrix of K cut , out ∈ R n × n , K out ∈ R n × n is the assembly of all the elements that fall totally outside of T and K GP ∈ R n × n is the part of sti ff ness matrix resulting from the ghost penalty description. Also the mass matrix can be rewritten in the same form. By bringing all the time dependent matrices to the right hand side, this gives

M b ¨ p ( t ) + K b p ( t ) = F ( t ) + ( K out ( t ) + K cut , out ( t ) − K GP ( t )) p ( t ) + ( M out ( t ) + M cut , out ( t ) − M GP ( t )) ¨ p ( t ) . (8)

As stated, an advantage of this form is that the DOFs stay the same when the source is moving. Furthermore, the absence of time-dependent matrices on the left hand side, in combination with a correctly chosen time-stepper, means that the required matrix inversion can be solved through a single LU-factorization. Also, it might be advantageous for an e ffi cient model order reduction procedure, since the terms on the left-hand side of the equation only have to be reduced once if the reduced basis is chosen to be time-independent.

3.2. Numerical time integration of the non-reduced system To assure the system is stable before reduction and to get a reference solution, several numerical time integration schemes are assessed. Besides the often used Newmark- β scheme, which has poor high frequency dissipation, the generalized- α scheme is investigated, since it has a tunable high frequency dissipation, which however might influence the accuracy of the solution if too much damping is applied [9]. Also several variants of these two integrators exist, where the integrator is combined with a Newton-Raphson procedure to iteratively improve the solution [10].

3.3. Time-dependent problem formulation

Volume acceleration [m 3 /s 2 ]

0.2

0.02

Pressure [Pa]

0.1

0

0

-0.02

-0.1

0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 Time [s]

0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 Time [s]

0

0

Magnitude [dB]

Magnitude [dB]

-20

-50

-40

Newmark- Generalized- implicit/explicit Gen-

-100

-60

0 100 200 300 400 500 600 700 800 900 1000 Frequency [Hz]

0 200 400 600 800 1000 Frequency [Hz]

Figure 5: Pressure at x = 4 m , y = 0 . 5 m calculated with Newmark- β and Generalized- α .

Figure 4: Bandlimited impulsive input signal, applied as normal Neumann BC on the boundary of the moving source.

A 2D problem is considered of an acoustic space that is 5 by 2 m, where all of the edges are sound hard walls. An acoustic source in the shape of a circle is placed in the middle and is moving with a constant velocity of v = 10 m / s in the positive x-direction, see Fig. 3, while emitting a pass- band filtered impulsive input signal with the cut-o ff frequencies at 50 and 500 Hz, see Fig. 4. The sampling frequency is set at f s = 5000 Hz. The domain is meshed with rectangular elements of 0.05 by 0.05 m with linear shape functions. The output is assessed at a node placed at x = 4 m , y = 0 . 5 m. At every time instance the cut elements are recalculated with the updated level-set function. Besides the proposed cutFEM formulation, the same problem was also recreated in the commercial

Figure 6: Comparison of full pressure field with Newmark- β (left) and generalized- α (right) at t = 0 . 08s. The pressure at the edges of the cut elements should be ignored, since they lie outside of the physical domain.

software package COMSOL TM , to compare the cutFEM solution to an ALE solution, which we use as a benchmark in Sec. 4. During simulation, COMSOL TM automatically remeshes the problem 6 times to prevent a too distorted mesh, while the cutFEM problem uses the same mesh throughout. Using a standard Newmark- β time stepper to solve Eq. (5) leads to noise amplification that becomes worse when time progresses, see Fig. 5. The linear generalized- α scheme directly applied to Eq. (7) performs much better, using a numerical damping factor of ρ ∞ = 0 . 5. However, the implicit / explicit generalized- α scheme that can solve Eq. (8) is much more e ffi cient, since a single LU-decomposition can be performed using the time-independent mass and sti ff ness matrix. In Fig. 5 the results are compared, where it can be seen that also the noise reduces, especially at low frequencies. Also, several full field solutions are compared, which confirm that generalized- α leads to less noisy full field responses, shown at 0.08 s (Fig. 6).

a

4. MODEL ORDER REDUCTION

0.02

Pressure [Pa]

0

-0.02

0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 Time [s]

0

ROM state-space ALE implicit/explicit Gen-

Magnitude [dB]

-20

-40

-60

0 100 200 300 400 500 600 700 800 900 1000 Frequency [Hz]

Figure 7: (Left) Pressure at x = 4 m , y = 0 . 5 m calculated with the reduced order state-space form and compared to the ALE results and implicit / explicit generalized- α . (Right) Full field response of state-space model (top) and of ALE form (bottom) at t = 0 . 08s.

While the solvers in Sec. 2 are implicit, also explicit solvers are of interest in combination with model order reduction (MOR) for real-time implementation. Specifically, a state-space form discretized with

exponential integration has shown its merit, as an e ffi cient and stable first order time integrator, as is shown in recent research [5]. To perform MOR, as a first step the reduced basis V ∈ R n × k with k << n has to be constructed. This basis can then be used to obtain the reduced order matrices, as follows

M r = V ⊤ MV , K r = V ⊤ KV , F r = V ⊤ F , p ≈ V p r . (9)

Thus giving reduced matrices M r ∈ R k × k , K r ∈ R k × k , and reduced vectors F r ∈ R k × 1 , x r ∈ R k × 1 . One of the choices that has to be made is if the reduced bases will be local (dependent on time) or global (a single reduced basis for all time steps). In this work we will mainly look into Krylov subspace reduction using the Second Order Arnoldi Recursion (SOAR), due to their previously proven performance for time-invariant vibro-acoustic systems [6]. However, this means that whenever the input changes, the optimal reduced basis changes and has to be recalculated. One way to obtain a global reduced basis would be to concatenate all of the optimal reduced bases for each time step and perform a singular value decomposition (SVD) to remove linearly dependent vectors, as is commonly done for parametric model order reduction [11]. However, the resulting SVD would be very expensive to calculate, unless a smart sampling scheme would be applied. Instead, for this paper it was chosen to use a time dependent reduced basis V ( t ) that is recalculated at each time step. The main disadvantages of this approach are that the reduced basis calculation is quite expensive and that at every time step the system matrices have to be projected to the reduced order form, which is an expensive operation that can take even more time than solving the unreduced problem. In follow-up research this problem will be addressed in more detail, since this is a common issue for the reduction of non-linear systems [12]. In this paper, we do not focus on performance, but we will only show that reduction and transformation to a state-space form is possible and leads to stable and accurate results. The required equations for the conversion of the reduced order model to a state-space form and the following time discretization step are not included in this paper for the sake of brevity, but can be found in [6]. As mentioned, the reduced matrices will change as function of time, thus a conversion to state-space form and time discretization has to be repeated at every time step. The results are shown in Fig. 7, where the results are directly compared to the results from the ALE simulation in COMSOL TM and to the results obtained with the generalized- α scheme. It has to be noted here that the discretized state space model does not contain any mechanism to dissipate high frequency noise. As can be seen, the resulting response has a similar response to the ALE results, while having lower noise outside of these bands than both the generalized- α scheme and the results from the ALE form (calculated with a BDF time stepper). Thus, while the current MOR strategy still has some remaining issues concerning computational e ffi ciency, this example shows that it is possible to achieve stability and a high accuracy through reduced order modeling.

5. CONCLUSIONS

To model the sound radiating from a moving (vibro-)acoustic source with the FEM, several challenges appear that are not present for static sources. One factor of major importance is how to perform meshing when the source is moving, while achieving high accuracy with minimal numerical artifacts. Additionally, it is desired to apply model order reduction to speed-up online calculations to integrate such a model in a physics-based virtual sensing framework. In this paper a cutFEM method is proposed that solves the time-varying problem on a static background mesh, thus not requiring any remeshing during simulation. Numerical artifacts are minimized by using a ghost penalty description in combination with a generalized- α scheme as time integrator. Furthermore, it is shown that stability preserving model order reduction is possible using a Krylov subspace reduction technique. The resulting time-dependent mass and sti ff ness matrices are converted to a state-space form, showing that this integration scheme leads to desirable high- and low-frequency dissipation, while retaining accuracy in the frequency bands of interest. Since the

reduction step currently takes more time than actually solving the full order model, future research will be focused on speeding up the reduction step, i.e. through hyper-reduction techniques.

ACKNOWLEDGEMENTS

The research of S. van Ophem (fellowship no. 1277021N) is funded by a grant from the Research Foundation – Flanders (FWO). The Flanders Innovation & Entrepreneurship Agency within the BE- QUIET project is gratefully acknowledged for its support. Internal Funds KU Leuven are gratefully acknowledged for their support.

REFERENCES

[1] O. Guasch, M. Arnela, R. Codina, and H. Espinoza. Stabilized finite element formulation for the mixed convected wave equation in domains with driven flexible boundaries. In Proceedings of NOVEM 2015 , April 2015. [2] E. Burman, S. Claus, P. Hansbo, M. G. Larson, and A. Massing. CutFEM: Discretizing geometry and partial di ff erential equations. International Journal for Numerical Methods in Engineering , 104(7):472–501, 2015. [3] C. B. Dilgen and N. Aage. Generalized shape optimization of transient vibroacoustic problems using cut elements. International Journal for Numerical Methods in Engineering , 122(6):1578– 1601, 2021. [4] A. Legay. An extended finite element method approach for structural-acoustic problems involving immersed structures at arbitrary positions. International Journal for Numerical Methods in Engineering , 93(4):376–399, 2013. [5] S. van Ophem and W. Desmet. Physics-based sound radiation estimation from multiple speakers by combined lumped parameter and reduced-order finite element modeling. Mechanical Systems and Signal Processing , 167:108585, March 2022. [6] S. van Ophem, O. Atak, E Deckers, and W Desmet. Stable model order reduction for time- domain exterior vibro-acoustic finite element simulations. Computer Methods in Applied Mechanics and Engineering , 325:240–264, October 2017. [7] S. Sticko and G. Kreiss. A stabilized Nitsche cut element method for the wave equation. Computer Methods in Applied Mechanics and Engineering , 309:364–387, September 2016. [8] E. Burman. Ghost penalty. Comptes Rendus Mathematique , 348(21):1217–1220, November 2010. [9] J. Chung and G. M. Hulbert. A Time Integration Algorithm for Structural Dynamics With Improved Numerical Dissipation: The Generalized- α Method. Journal of Applied Mechanics , 60(2):371–375, June 1993. [10] W. J. T. Daniel. Explicit / implicit partitioning and a new explicit form of the generalized alpha method. Communications in Numerical Methods in Engineering , 19(11):909–920, 2003. [11] P. Benner, S. Gugercin, and K. Willcox. A Survey of Projection-Based Model Reduction Methods for Parametric Dynamical Systems. SIAM Review , 57(4):483–531, January 2015. [12] C. Farhat, T. Chapman, and P. Avery. Structure-preserving, stability, and accuracy properties of the energy-conserving sampling and weighting method for the hyper reduction of nonlinear finite element dynamic models. International Journal for Numerical Methods in Engineering , 102(5):1077–1110, 2015.