The time-explicit nodal discontinuous Galerkin method applied toacoustic-structure interaction problems
Kirill SHAPOSHNIKOV(1) , Mads J. HERRING JENSEN(2) , Elin SVENSSON(3)(1) COMSOL A/S, Denmark, kirill@comsol.dk(2) COMSOL A/S, Denmark, mads@comsol.dk(3) COMSOL AB, Sweden, elins@comsol.seAbstract
In this work, we apply the time-explicit nodal discontinuous Galerkin (dG) method to study wave propagationin coupled linear elastic and acoustic media. The velocity-strain formulation is used for the structural part of thesystem, and the velocity-pressure formulation, for the acoustic part. The interface boundary conditions between
the structural and acoustic domains are treated by imposing the proper upwind fluxes. The method is validatedby comparison to the results obtained for the second order time-implicit finite element (FEM) formulation.
Keywords: Discontinuous Galerkin, Time-explicit, Acoustic-structure interaction

1 INTRODUCTION
The time-explicit nodal discontinuous Galerkin method has gained a good reputation as a numerical approach
to solving first-order hyperbolic problems for its high accuracy and low memory consumption, great scalability
on clusters, and ability to operate on unstructured as well as nonconforming meshes. The method has been
proven to be efficient for simulation of transient acoustic phenomena in acoustically large computation domains.
Application areas include ultrasound transducers, nondestructive testing, geophysics, etc.
We apply the velocity-strain and the velocity-pressure formulation to the elastic and acoustic problems, respectively. In the three-dimensional space, this results in nine and four coupled hyperbolic equations for thestructural and the acoustic part. Even though the number of unknowns is higher compared to that used in thesecond-order time-implicit FEM formulation, the time-explicit dG method becomes more effective and in somecases the only applicable approach to large transient coupled acoustic-structure problems. The consistent andaccurate multiphysics coupling between the acoustic and structure is ensured through the proper upwind flux
across the interfaces between the elastic/elastic, elastic/acoustic, and acoustic/acoustic domains [1, 2].
The proposed approach is validated by comparison with results obtained for the second-order time-implicit FEM
formulation of the acoustic-structure systems. The effectiveness of the time-explicit dG scheme is demonstrated
by its application to large-scale wave propagation problems in coupled acoustic-elastic media.

2 ACOUSTIC WAVES IN SOLIDS AND FLUIDS
2.1 Governing equations
Let Ωe ⊂ R3 be a domain occupied by a linear elastic material with the density ρ and the fourth order tensor
of elasticity [C]. Assuming infinitesimal deformations, the mechanical displacement field, u m , the strain tensor,[E], and the Cauchy stress tensor, [S], are related to each other as follows1i1h∇uum (xx,t) + (∇uum (xx,t))T ,2[S (xx,t)] = [C (xx)] : [E (xx,t)].

[E (xx,t)] =1 Further on, the arguments are only used if necessary.

5450

(1)
(2)

Eqs. (1), (2) can be rewritten in the vector-matrix form according to Voigt notation
E = [C]Buu
S = [C]E(3)

with
S = sxxsyyszzsyz

E = exx∂eyyezz

2eyz

2exz

0∂∂z∂
∂z∂∂y0∂x
B=000∂∂y000

∂
∂zsxz

∂∂xsxyT,2exy
T∂∂y∂ ∂x 

T,,0and [C] is the 6 × 6 symmetric stiffness tensor.
The equations of linear elasticity written as the first order hyperbolic system in Ωe readρ∂v− BT S = f ,∂t∂E− Bvv = 0 ,
∂t

(4)
(5)

where v = ∂ u m /∂t is the structural velocity and f is the external volume force density.
Let Ωa ⊂ R3 be a domain occupied by an acoustic fluid with the density ρa and the speed of sound c. Here, v
becomes the acoustic particle velocity, u , and the stress tensor reduces to the scalar acoustic pressure, p, which
results in the following stress-strain relation equivalent to Eq. (2)− p = ρa c2

tr ([E]).
3(6)Substituting Eq. (6) into Eqs. (4), (5) yields the first order hyperbolic system of acoustic equations in Ωa∂u+ ∇p = f ,∂t
1 ∂p+ ∇ · u = 0.
ρa c2 ∂t
ρa(7)(8)

Eqs. (4), (5) and Eqs. (7), (8) are respectively supplemented by the initial conditionsv (xx, 0) = v 0 (xx) , E (xx, 0) = E 0 (xx) ,u (xx, 0) = u0 (xx) , p (xx, 0) = p0 (xx) .
Moreover, the interface boundary conditions apply the continuity constraints upon the boundaries between theelastic and acoustic subdomains as follows+−+
S−n = Sn , v = v+
on Γ = Ω−e ∩ Ωe ,(9)+p n = p+ n, u−n = un
S n = −pnn, vn = un+on Γ = Ω−a ∩ Ωa ,(10)on Γ = Ωe ∩ Ωa .

(11)

−

2.2 Weak form and discontinuous Galerkin method
Both Eqs. (4), (5) and Eqs. (7), (8) correspond to a hyperbolic system of conservation law that can be written
as follows∂qQ+ ∇ · Γ (qq) = f .
(12)∂t5451

Let Ω ' Ωh = k Dk , and let q 0 (xx) be a locally defined test function on Dk . Then the locally defined weak andstrong formulations of Eq. (12) read [3]Z ZZ∂q 0· q − Γ (qq) · ∇qq0 dxx +(nn · Γ ∗ (qq)) · q 0 dxx =f · q 0 dxx(13)Q∂tDk∂ DkDkSandZZ∂q 00Γ∗ (qq) − Γ (qq))) · q 0 dxx =Q· q + (∇ · Γ (qq)) · q dxx +(nn · (Γf · q 0 dxx,∂tDk∂ DkDkZ(14)respectively, with Γ ∗ being the numerical flux across ∂ Dk .
The choice of the numerical flux depends on the problem at hand. In general, the numerical flux along thenormal n can be expressed asΓ∗ (qq)} + F q − , q + ,n · Γ ∗ (qq) = n · {Γwhere {ξ } = 1 /2 (ξ − + ξ + ) is the average operator and F defines the flux type [3]. The central flux is obtainedchoosing F = 0 and the Lax-Friedrichs flux is given by

C − −C+ +−1 ∂ Γ
Q q − Q q , C = max Λ Q.
(15)
F = [Qqq] =q
22∂q
The definition of C given by Eq. (15) yields that it corresponds to the (quasi) pressure wave velocity.
Both the central and Lax-Friedrichs fluxes are suitable for the system of acoustic equations. For the equationsof linear elasticity, Ye et al. [4] suggest a hybrid flux that acts as the central flux along the normal direction andas Lax-Friedrichs along the tangential direction. Thus the hybrid flux treats the normal and tangential conditions
across the interface separately ensuring that no restrictions are applied to the shear stress.
2.3 Upwind numerical flux
Both the central and Lax-Friedrichs fluxes become inapplicable if the material parameters jump across the interface, because the proper upwind numerical flux is not obtained. However, it is obtained from the solution of
the Riemann problem across the interface [5].
The linearized form of Eq. (12) readsQ∂q∂q∂q∂q+ Ax+ Ay+ Az= f.
∂t∂x∂y∂z(16)

Let the material properties in the inward and the outward directions of the normal vector n be denoted by thesymbols “−” and “+”, respectively. Then the Riemann problem across the interface is to find a solution toEq. (16) satisfying the initial condition [1](q − (xx) , n · x < 0,
q 0 (xx) =q + (xx) , n · x > 0.
Let A = nx Ax + ny Ay + nz Az , then the numerical flux in the normal direction reads Aqq. The characteristics ofhyperbolic system are given by the spectral properties of the matrix Q−1 A, that is, by its spectral decomposition
ΛR−1 .
Q−1 A = RΛ(17)where the columns of R denote the eigenvectors corresponding to the eigenvalues, which in case of a generalanisotropic linear elastic material readΛ = diag (−c p , −cs,1 , −cs,2 , 0, 0, 0, c p , cs,1 , cs,2 ) .

5452

Figure 1. To the definition of unknown states for the Riemann problem
Note that cs,1 = cs,2 = cs for isotropic linear elastic material, and the only nonzero eigenvalues for the acoustic
fluid are ±c.
In case of anisotropic linear elastic materials on either side of the interface, the system has six unknown states,
as shown in Figure 1. The Rankine-Hugoniot jump conditions [5] read


−−
c−
q− − q−q− − q−pQ
a +A
a = 0,
−−−−c−q−q−
a − qb + A
a − q b = 0,
s,1 Q


−
∗−∗
q−
q−= 0,
c−s,2 Qb −q +Ab −q
A− q ∗ − A+ q ∗∗ = 0,++−c+q ∗∗ − q +q ∗∗ − q +s,2 Q
b +Ab = 0,++++q+−c+q+s,1 Q
b − q a = 0,b − qa + A
+++−c+q++ A+ q += 0.
pQa −qa −q
Taking Eq. (17) into account, this results in− −− −q ∗ − q − = α1− r −1 + α2 r 2 + α3 r 3 ,+ ++ +q + − q ∗∗ = α1+ r +1 + α2 r 2 + α3 r 3 ,where {rr ±i } are nonzero eigenvectors on either side of the interface. The proper upwind numerical flux is thenobtained as
1 3+ + + +n · Γ ∗ (qq) = {Aqq} − ∑ αi− λi− Q− r −(18)i + αi λi Q r i .
2 i=1
The coefficients {αi± } are then defined from the interface boundary conditions Eqs. (9)–(11). A closed form ofEq. (18) is derived in [1] on elastic/elastic, acoustic/elastic, and acoustic/acoustic interfaces in case of isotropic
linear elastic material. For a general anisotropic linear elastic material, an exact Riemann solver is developedin [2]. The solver yields the solution to the Riemann problem in the vector-matrix form.

3 IMPLEMENTATION AND NUMERICAL TESTSThe developed numerical algorithm is implemented in COMSOL Multiphysicsr software [6] using the Wave
Form PDE physics interface. The interface is dedicated to solving hyperbolic partial differential equations ofthe form Eq. (12) using nodal discontinuous Lagrange shape functions and time explicit integration.
For the problem at hand the implementation results in 9 and 4 unknowns for the elastic and acoustic parts,5453

respectively. We use the fourth order approximation of the unknowns that yields the largest mesh element size
cmincmin,.
(19)
hmax ∈
2 fmax 1.5 fmax
Eq. (19) provides an unambiguous expression of the mesh element size for the system of acoustic equations,because of the acoustic waves propagate with the same speed of sound, c. This is not the case for the problemof linear elasticity, where cmin can be different depending on the problem formulation. This is shown in theexample below, where the elastic waves propagate in the lower half-plane with a “rock” on its top that issimilar to the one discussed in [7]. The domain is occupied by an isotropic linear elastic material with theproperties ρ = 2200 kg/m3 , c p = 3.2 km/s, and cs = 1.8475 km/s. The traction free boundary condition applied
to the top surface results in the propagation of a Rayleigh wave that travels in the shallow layer below thesurface [8]. Its velocity, cr , can be approximately calculated ascr ≈0.86 + 1.12ν
cs1+ν

(20)

in case of an isotropic linear elastic material with ν being the Poison’s ration. In our example, Eq. (20) yields
cr ≈ 1.7 km/s.
The computational domain used for the numerical tests is the rectangle of the size 40 × 20 km. The “rock”
on its top has the shape of a Gaussian of the heigh 0.2 km and the standard deviation 0.3 km. Its centre islocated 6 km to the left of the sagittal line. The simulation is run for 12 s and the forth-order time explicit
Runge-Kutta integration scheme is used. The source is given at the center of the top side of the rectangle asthe point force in the y-direction
F (t) = Fm 2 (π f (t − t0 ))2 − 1 exp − (π f (t − t0 ))2 ,where Fm = 10 MN, f = 2 Hz, and t0 = 1[s]. The Absorbing Layers that have the same functionality as in
the Acoustics Module of COMSOL Multiphysics [9] ensure the propagation in unbounded domains with low
reflections. The layers are adjacent to the left, right, and bottom sides of the rectangle and have the widthof 2 km. The impedance-like boundary condition, Low-Reflecting Boundary [10], is imposed on the outerboundaries of the absorbing layer domains.

(a)

(b)

(c)

Figure 2. Structural velocity magnitude profile at 5 s (a), 8 s (b), and 11 s (c)Figure 2 shows the resulting velocity magnitude profiles at three different time steps. One can see a goodresolution of the waves travelling through the media: longitudinal, shear, von Schmidt, and Rayleigh waves. Itis also seen that there are no distinguishable reflections from the outer boundary. Figure 3 depicts the velocity
profiles on the domain top boundary. The results are obtained by the time-explicit dG and the second-ordertime-implicit FEM [9, 10] approaches (red and blue curves, respectively) and are normalized to 10 km/s. The
comparison shows a good agreement between the two numerical techniques. At the same time, the former is5454

(a)

(b)

Figure 3. Structural velocity profiles at t = 8 s: x-component (a), y component (b)more effective in the considered example. Either method takes about 2.5 hours to run the simulation on a workstation with Intelr CoreTM i7 CPU (3.60 GHz) and 32 GB of RAM. However, the implemented dG approachrequires only 2.4 GB of RAM, which is much less that 15.5 GB used by FEM. This memory efficiency isachieved despite the larger number of degrees of freedom solved for (3.5 · 106 for dG vs. 1.7 · 106 for FEM)
and the smaller time step (∆texpl ≈ ∆timpl /4). This suggests that the time-explicit dG approach will be the only
possible option for solving such problems in larger computational domains, especially in the three-dimensional
space, and time intervals. Moreover, the time-explicit dG method is better suited for the parallelization onclusters than the second-order time-implicit FEM.

(a)

(b)

Figure 4. The pressure distribution in the cylinder together with the scattered (a) and the total (b) acousticpressure around it at t = 4/ f0 s
The second example shows the scattering of an acoustic wave off a solid aluminium cylinder. That is, theacoustic-structure interaction is taken into account. The background acoustic field is given by a Gaussian modu-5455

lated plane wave of the center frequency f0 = 2 kHz propagating in the air. The radius of the scatterer is 0.1 m.
The plots of the scattered and the total acoustic pressure fields are depicted in Figure 4 (a) and (b), respectively.
One can see that the total acoustic pressure field and the pressure filed in the elastic scatterer are continuous
across the acoustic-structure interface. That is, the condition given by Eq. (11) is fulfilled by means of theupwind flux.

4 CONCLUSIONS
The described time-explicit dG method to solving large transient acoustic-structure multiphysics problems gives
good results when compared to the second-order time-implicit FEM approach. Its advantages are high accuracy,
low memory consumption, and scalability on clusters. The flexibility of COMSOL Multiphysics makes it possible to manually implement the governing equations as well as the upwind fluxes, ensuring correct coupling
between the structure and acoustics.
Further work will include extension of the proposed approach to the cases where the acoustic field is convectedby a steady-state background flow.

REFERENCES[1] Wilcox, L. C.; et al. A high-order discontinuous Galerkin method for wave propagation through coupledelastic-acoustic media. J. Comp. Phys., 229, 2010, pp. 9373–9396.
[2] Zhan, Q.; et al. An exact Riemann solver for wave propagation in arbitrary anisotropic elastic media withfluid coupling. Comput. Methods. Appl. Mech. Engng., 329, 2018, pp. 24–39.
[3] Hesthaven, J. S.; Warburton, T. Nodal discontinuous Galerkin methods: Algorithms, analysis, and applications. Springer Verlag, New York, 2008.
[4] Ye, R.; et al. A discontinuous Galerkin method for the anisotropic elastic wave equation. Proceedings ofthe Project Review, Geo-Mathematical Imaging Group (Purdue University, West Lafayette IN), Vol. 1, 2013,
pp. 103–130.
[5] LeVeque, R. J. Finite-volume methods for hyperbolic problems. Cambridge University Press, Cambridge,
2004.
[6] COMSOL Multiphysicsr v. 5.4, www.comsol.com, COMSOL AB, Stockholm, Sweden.
[7] Appelö, D.; Petersson, N. A. A stable finite difference method for the elastic wave equation on complex
geometries with free surfaces. Commun. Comput. Phys., 5(1), 2009, pp. 84–107.
[8] Viktorov, I. A. Rayleigh and Lamb waves: Physical theory and applications. Springer Verlag, New York,2013.
[9] Acoustics Module, COMSOL Multiphysicsr v. 5.4, www.comsol.com, COMSOL AB, Stockholm, Sweden.
[10] Structural Mechanics Module, COMSOL Multiphysicsr v. 5.4, www.comsol.com, COMSOL AB, Stockholm, Sweden.

5456