Virtual elements for sound propagation in complex poroelastic media

Abstract

We develop a novel Virtual Element Method (VEM) to resolve the mixed Biot displacement-pressure formulation governing wave propagation in porous media. Within this setting, the weak form of the governing equations is discretized using implicitly defined canonical basis functions and the resulting integral forms are computed using appropriate polynomial projections. The projection operator accounting for the solid, fluid, and coupling phases of the problem are presented. Different boundary, interface and excitation conditions are accounted for. The convergence behaviour, accuracy, and efficiency of the method is examined through a set of illustrative examples. A node insertion strategy is proposed to resolve non-conforming interfaces that naturally arise in multilayered systems. Finally the power of the VEM is exploited to examine the acoustic response of composite materials with periodic and non-periodic inclusions of complex geometries.

Appendices

Appendix A monomial spaces

The contents of the monomial spaces $[{\mathbb {M}}_k({\mathcal {K}})]^2$ and $[{\mathbb {M}}_k({\mathcal {K}})]$ are iteratively defined in Table 11.

Table 11 Generalized scalar and vector valued monomials for ${\mathbb {M}}_k({\mathcal {K}})$ and $[{\mathbb {M}}_k({\mathcal {K}})]^2$, respectively

In Table 11, $\xi = \frac{x-x_{{\mathcal {K}}}}{h_{{\mathcal {K}}}}$ and $\eta = \frac{y-y_{{\mathcal {K}}}}{h_{{\mathcal {K}}}}$ denote scaled monomials in each parametric direction. The number of terms in $[{\mathbb {M}}_k({\mathcal {K}})]^2$ and $[{\mathbb {M}}_k({\mathcal {K}})]$ are $n^u_k = (k+1)(k+2)$ and $n^p_k = \frac{(k+1)(k+2)}{2}$, respectively. The operator-specific kernels are provided in Table 12. The contents of these kernels can be derived using kinematical decomposition relations mentioned in [73].

Table 12 Definition of operator kernels

The zero-energy modes contained in ${\mathbb {K}}^{\varepsilon }({\mathcal {K}})$ can be understood as rigid body motions, i.e., two translations and one rotation in 2-D physical space.

Appendix B computing B matrices

Using quadratures, Eq. (66) results in

$$\begin{aligned} {\mathbf {B}}^{\nabla u}_j={\mathbf {B}}_{bj}^{\nabla u}+{\mathbf {B}}_{dj}^{\nabla u}, \end{aligned}$$

(128)

where

$$\begin{aligned} {\mathbf {B}}^{\nabla u}_{bj}=\sum \limits _{e \in \partial {\mathcal {K}}}\int \limits _{e} {\mathbf {u}}_h \cdot (\delta ^{\star } \nabla {\mathbf {m}}_{j+2}\, {\mathbf {n}}^{\nabla }_e) \, \text {d}e, \end{aligned}$$

(129)

and

$$\begin{aligned} {\mathbf {B}}^{\nabla u}_{dj}=-\int \limits _{{\mathcal {K}}\setminus \partial {\mathcal {K}}} {\mathbf {u}}_h \cdot (\delta ^{\star }\varDelta {\mathbf {m}}_{j+2}) \, \text {d}{\mathcal {K}} \end{aligned}$$

(130)

respectively.

Similarly, Eq. (69) becomes

$$\begin{aligned} {\mathbf {B}}^{\nabla p}_j= {\mathbf {B}}_{bj}^{\nabla p}+{\mathbf {B}}_{dj}^{\nabla p}, \end{aligned}$$

(131)

where

$$\begin{aligned} {\mathbf {B}}^{\nabla p}_{bj}=\sum \limits _{e \in \partial {\mathcal {K}}}\int \limits _{e} \text {p}_h \cdot \Big ( \frac{1}{{\tilde{\rho }}_{\text {eq}}} \nabla \text {m}_{j+1}\cdot {\mathbf {n}}^{\nabla }\left( e\right) \Big ) \, \text {d}e, \end{aligned}$$

(132)

and

$$\begin{aligned} {\mathbf {B}}^{\nabla p}_{dj}=-\int \limits _{{\mathcal {K}}\setminus \partial {\mathcal {K}}} \text {p}_h \cdot \Big (\frac{1}{{\tilde{\rho }}_{\text {eq}}}\varDelta \text {m}_{j+1} \Big ) \, \text {d}{\mathcal {K}} \end{aligned}$$

(133)

respectively, where the array ${\mathbf {n}}^{\nabla }\left( e\right)$ is defined here as

$$\begin{aligned} {\mathbf {n}}^{\nabla }\left( e\right) = \begin{bmatrix} n_x\left( e\right)&n_y\left( e\right) \end{bmatrix}^T. \end{aligned}$$

(134)

The boundary integrals in Eqs. (128) and (131) are evaluated using Gauss-Lobatto quadratures in a similar manner to Eq. (74).

To evaluate the corresponding domain integrals, the terms $\varDelta {\mathbf {m}}_{j+2}$ in Eq. (130) is expanded over the $[{\mathbb {M}}_{k-2}({\mathcal {K}})]^2$ basis

$$\begin{aligned} \delta ^{\star } \varDelta {\mathbf {m}}_{j+2} = \sum \limits _{\beta =1}^{n^u_{k-2}} \text {d}_{j \beta }^{\nabla u} {\mathbf {m}}_{\beta } , \quad \forall \, {\mathbf {m}}_{\beta } \in [{\mathbb {M}}_{k-2}({\mathcal {K}})]^2. \end{aligned}$$

(135)

Inserting Eq. (135) in Eq. (130) the following expression is derived

$$\begin{aligned} {\mathbf {B}}^{\nabla u}_{dj} = - |{\mathcal {K}}| \sum \limits _{\beta = 1}^{n^u_{k-2}} \text {d}_{j \beta }^{\nabla u} \text {dof}_{kN_v+\beta }({\mathbf {u}}_h) = - |{\mathcal {K}}| {\mathbf {d}}^{\nabla u}, \end{aligned}$$

(136)

where

(137)

Similarly, expanding $\varDelta \text {m}_{j+1}$ over the basis ${\mathbb {M}}_{k-2}({\mathcal {K}})$

$$\begin{aligned} \frac{1}{{\tilde{\rho }}_{\text {eq}}}\varDelta \text {m}_{j+1} = \sum \limits _{\beta =1}^{n^p_{k-2}} \text {d}_{j \beta }^{\nabla p} \text {m}_{\beta } , \quad \forall \, \text {m}_{\beta } \in {\mathbb {M}}_{k-2}({\mathcal {K}}), \end{aligned}$$

(138)

where the coefficients $\text {d}_{j \beta }^{\nabla p}$ are also obtained through inspection and substituting in Eq. (133)

$$\begin{aligned} {\mathbf {B}}^{\nabla p}_{dj} = - |{\mathcal {K}}| \sum \limits _{\beta = 1}^{n^p_{k-2}} \text {d}_{j \beta }^{\nabla p} \text {dof}_{kN_v+\beta }(\text {p}_h) = - |{\mathcal {K}}| {\mathbf {d}}^{\nabla p}, \end{aligned}$$

(139)

where

(140)

Appendix C computing D matrices

The matrices ${\mathbf {D}}^{{\varvec{\varepsilon }}}$, ${\mathbf {D}}^{0 u}$, ${\mathbf {D}}^{\nabla p}$, and ${\mathbf {D}}^{0 p}$ in Eqs. (105)–(108) assume the following form

$$\begin{aligned}&\begin{aligned}{\mathbf {D}}^{{\varvec{\varepsilon }}} = \begin{bmatrix} \text {dof}_1({\mathbf {m}}_1) &{} \ldots &{} \text {dof}_1({\mathbf {m}}_{n_k^u-3}) \\ \vdots &{} \ddots &{} \vdots \\ \text {dof}_{n^u_{\text {dof}}}({\mathbf {m}}_1) &{} \ldots &{} \text {dof}_{n^u_{\text {dof}}}({\mathbf {m}}_{n_k^u-3}) \end{bmatrix}, \\ \forall \, {\mathbf {m}} \in [{\mathbb {M}} _{\text {k}}({\mathcal {K}})]^2 \setminus {\mathbb {K}}^{{\varvec{\varepsilon }}} ({\mathcal {K}}) \end{aligned} \end{aligned}$$

(141a)

$$\begin{aligned}&\begin{aligned}{\mathbf {D}}^{0 u} = \begin{bmatrix} \text {dof}_1({\mathbf {m}}_1) &{} \ldots &{} \text {dof}_1({\mathbf {m}}_{n_k^{u}}) \\ \vdots &{} \ddots &{} \vdots \\ \text {dof}_{n^u_{\text {dof}}}({\mathbf {m}}_1) &{} \ldots &{} \text {dof}_{n^u_{\text {dof}}}({\mathbf {m}}_{n_k^{u}}) \end{bmatrix}, \\ \forall \, {\mathbf {m}} \in [{\mathbb {M}}_{\text {k}}({\mathcal {K}})]^2 \end{aligned}\end{aligned}$$

(141b)

$$\begin{aligned}&\begin{aligned}{\mathbf {D}}^{\nabla p} = \begin{bmatrix} \text {dof}_1(\text {m}_1) &{} \ldots &{} \text {dof}_1(\text {m}_{n_k^{p}-1}) \\ \vdots &{} \ddots &{} \vdots \\ \text {dof}_{n^p_{\text {dof}}}(\text {m}_1) &{} \ldots &{} \text {dof}_{n^p_{\text {dof}}}(\text {m}_{n_k^{p}-1}) \end{bmatrix}, \\ \forall \,\text {m} \in {\mathbb {M}}_{\text {k}}({\mathcal {K}}) \setminus {\mathbb {K}}^{\nabla }({\mathcal {K}}) \end{aligned}\end{aligned}$$

(141c)

$$\begin{aligned}&\begin{aligned}{\mathbf {D}}^{0 p} = \begin{bmatrix} \text {dof}_1(\text {m}_1) &{} \ldots &{} \text {dof}_1(\text {m}_{n_k^{p}}) \\ \vdots &{} \ddots &{} \vdots \\ \text {dof}_{n^p_{\text {dof}}}(\text {m}_1) &{} \ldots &{} \text {dof}_{n^p_{\text {dof}}}(\text {m}_{n_k^{p}}) \end{bmatrix}, \\ \forall \, \text {m} \in {\mathbb {M}}_{\text {k}}({\mathcal {K}}). \end{aligned} \end{aligned}$$

(141d)

The quantities $\text {dof}_i({\mathbf {m}}_j)$ and $\text {dof}_i(\text {m}_j)$ are evaluated according to the following expressions (see Table 4), for the solid

$$\begin{aligned} {\left\{ \begin{array}{ll} \text {dof}_i({\mathbf {m}}_j) = {\mathbf {m}}_j({\mathbf {x}}_i), \quad \forall \, i \le 2kN_v \\ \text {dof}_i({\mathbf {m}}_j) = \frac{1}{|{\mathcal {K}}|} \int \limits _{{\mathcal {K}}} {\mathbf {m}}_j \cdot {\mathbf {m}}_{\beta } \, \text {d}{\mathcal {K}}, \\ \forall \, {\mathbf {m}}_{\beta } \in [{\mathbb {M}}_{k-2}({\mathcal {K}})]^2, \quad i >2kN_v, \end{array}\right. } \end{aligned}$$

(142)

and the fluid phase

$$\begin{aligned} {\left\{ \begin{array}{ll} \text {dof}_i(\text {m}_j) = \text {m}_j({\mathbf {x}}_i), \quad \forall \, i \le kN_v \\ \text {dof}_i(\text {m}_j) = \frac{1}{|{\mathcal {K}}|} \int \limits _{{\mathcal {K}}} \text {m}_j \cdot \text {m}_{\beta } \, \text {d}{\mathcal {K}}, \\ \forall \, \text {m}_{\beta } \in {\mathbb {M}}_{k-2}({\mathcal {K}}), \quad i > kN_v \end{array}\right. }, \end{aligned}$$

(143)

respectively.

Appendix D computing sound absorption and transmission loss coefficients

The complex valued surface impedance at the incident face normalized with respect to the impedance of air $z(\theta )$ is computed

$$\begin{aligned} Z_{\text {sn}\, {\mathcal {K}}}(\omega ,\theta ) = \Big (\frac{\text {p}^{\text {in}}_{{\mathcal {K}}}}{\text {v}^{\text {in}}_{n{\mathcal {K}}}} \Big )_{\varGamma ^{{\mathcal {K}}}_{hI}} \big / z(\theta ) \end{aligned}$$

(144)

where $\text {p}^{\text {in}}_{{\mathcal {K}}}$ and $\text {v}^{\text {in}}_{n{\mathcal {K}}}$ denote inlet pressures and normal component of fluid velocities over the elementary incident face ${\varGamma ^{{\mathcal {K}}}_{hI}}$. This quantity is now used to obtain the elementary coefficient of reflection:

$$\begin{aligned} R_{{\mathcal {K}}}(\omega ,\theta ) = \frac{ Z_{\text {sn}\, {\mathcal {K}}}(\omega ,\theta )-1}{ Z_{\text {sn}\, {\mathcal {K}}}(\omega ,\theta )+1} \end{aligned}$$

(145)

As the VEM computes only resultant quantities (the net incident and reflected waves), one requires $R_{{\mathcal {K}}}(\omega ,\theta )$ to obtain purely incident pressures and normal fluid velocities over ${\varGamma ^{{\mathcal {K}}}_{hI}}$.

$$\begin{aligned} \text {p}^{\text {inc}}_{{\mathcal {K}}} = \bigg | \frac{\text {p}^{\text {in}}_{{\mathcal {K}}}}{1+ R_{{\mathcal {K}}}(\omega ,\theta )}\bigg |_{\varGamma ^{{\mathcal {K}}}_{hI}}, \quad \text {v}^{\text {inc}}_{n {\mathcal {K}}} = \frac{\text {p}^{\text {inc}}_{{\mathcal {K}}}}{z(\theta )} \end{aligned}$$

(146)

The time averaged powers are evaluated according to Eqs. (147) below

$$\begin{aligned}&\begin{aligned} {\mathbb {W}}^{\text {in}}(\omega ,\theta ) = \frac{1}{2} \mathfrak {R}\Bigg (\int \limits _{\varGamma _{hI}} \text {p}^{\text {in}} \cdot \text {v}_n^{\text {in} \, *}\, d \varGamma \Bigg )= \\ \frac{1}{2} \mathfrak {R}\Bigg (\sum \limits _{i} \int \limits _{\varGamma _{hI}^{{\mathcal {K}}(i)}} \text {p}^{\text {in}}_{{\mathcal {K}}(i)} \cdot \text {v}_{n{\mathcal {K}}(i)}^{\text {in} \, *}\, d \varGamma \Bigg ), \end{aligned} \end{aligned}$$

(147a)

$$\begin{aligned}&\begin{aligned} {\mathbb {W}}^{\text {inc}}(\omega ,\theta ) = \frac{1}{2} \int \limits _{\varGamma _{hI}} \text {p}^{\text {inc}} \cdot \text {v}^{\text {inc}}_n\, d \varGamma = \\ \frac{1}{2} \sum \limits _{i} \int \limits _{\varGamma _{hI}^{{\mathcal {K}}(i)}} \text {p}^{\text {inc}}_{{\mathcal {K}}(i)} \cdot \text {v}^{\text {inc}}_{n{\mathcal {K}}(i)}\, d \varGamma , \end{aligned} \end{aligned}$$

(147b)

$$\begin{aligned}&{\mathbb {W}}^{\text {ref}}(\omega ,\theta ) = {\mathbb {W}}^{\text {inc}}(\omega ,\theta )-{\mathbb {W}}^{\text {in}}(\omega ,\theta ), \end{aligned}$$

(147c)

$$\begin{aligned}&\begin{aligned} {\mathbb {W}}^{\text {trans}}(\omega ,\theta ) = \frac{1}{2} \mathfrak {R}\Bigg (\int \limits _{\varGamma _{hO}} \text {p}^{\text {out}} \cdot \text {v}_n^{\text {out} \, *}\, d \varGamma \Bigg )= \\ \frac{1}{2} \mathfrak {R}\Bigg (\sum \limits _{i} \int \limits _{\varGamma _{hO}^{{\mathcal {K}}(i)}} \text {p}^{\text {out}}_{{\mathcal {K}}(i)} \cdot \text {v}_{n{\mathcal {K}}(i)}^{\text {out} \, *}\, d \varGamma \Bigg ) , \end{aligned} \end{aligned}$$

(147d)

where ${\mathbb {W}}^{\text {in}}$, ${\mathbb {W}}^{\text {inc}}$, ${\mathbb {W}}^{\text {ref}}$ and ${\mathbb {W}}^{\text {trans}}$ represent inlet, incident, reflected and transmitted powers, respectively. The operator $\mathfrak {R}(\cdot )$ extracts real valued data, Complex conjugation is denoted by $(*)$. Outlet fluid pressures and normal components of fluid velocity $\text {p}^{\text {out}}_{{\mathcal {K}}(i)}$ and $\text {v}_{n{\mathcal {K}}(i)}^{\text {out}}$ are evaluated over an elementary outlet face ${\varGamma _{hO}^{{\mathcal {K}}(i)}}$. The SAC and STL are finally derived for a plane wave incident at an angle $\theta$ with a driving angular frequency of $\omega$:

$$\begin{aligned} \alpha (\omega ,\theta ) \!=\! 1 - \frac{{\mathbb {W}}^{\text {ref}}(\omega ,\theta )}{{\mathbb {W}}^{\text {inc}}(\omega ,\theta )}, \quad {\mathcal {T}}(\omega ,\theta ) \!=\! 10 \text {log}\frac{{\mathbb {W}}^{\text {inc}}(\omega ,\theta )}{{\mathbb {W}}^{\text {trans}}(\omega ,\theta )}. \end{aligned}$$

(148)

For a detailed report investigating the post-processing procedures involved in structural and porous vibro-acoustics, see, e.g., [89, 97].