Transient Acoustic Wave Propagation Problems in Multilayered Pavement Using a Time Discontinuous Galerkin Finite Element Method

Abstract

This paper presents a modified time discontinuous Galerkin finite element method (MDGFEM) for transient acoustic wave propagation problem of multilayered pavement. The pavement consists of cement concrete pavement, semi-rigid base, and natural soil. The multilayered pavement is modeled as poroelastic mediums and assumed to be water-saturated. The well-known generalized Biot’s theory is employed to describe the wave propagation problem. The present MDGFEM, based on the artificial damping scheme, employs the Hermite (P3) functions and the linear (P1) shape functions to interpolate the global nodal vector and its temporal gradient respectively in a time step. Numerical results of 1D and 2D problems show that the MDGFEM can filter out spurious numerical oscillations before and after waves, boundaries of the hole, and the interface between the layers more effectively for the propagation of acoustic waves in multilayered pavement. Compared with widely used time-continuous methods such as the Newmark method, the method proposed in this paper presents better capabilities in the fluid–structure interaction behavior of multilayer pavements and provides a more accurate solution, which contributes to the further development of non-destructive testing of multilayer pavement structures.

1. Introduction

Cement concrete pavement is one of the main structural forms of highway pavement in the world and ultrasound systems are widely used in the detection of internal diseases in multilayer pavement [1]. The pavement structure is always composed of multiple layers of inhomogeneous porous materials. Due to the existence of internal defects such as cracks and voids, the propagation of transient acoustic waves in the vicinity of the defects changes, and ultrasound systems can identify the defects through these changes in the non-destructive testing of pavement [2]. On the other hand, the interaction between ultrasound and pavement is very complex [3], which motivates our interest in the study of wave propagation laws in multilayer pavement structures.

In the past few years, there has been a lot of effort to study the effect of pore water pressure in multilayer pavements. Cai et al. [4,5] described the mechanical response of a fully saturated poroelastic half-space pavement under different velocity loads based on Biot’s dynamic poroelasticity theory [4,5,6,7]. The results show that the influence of the pore water becomes more pronounced as the truck speed increases. Due to its simplicity and practicality, Biot’s method is widely used in the problem of wave propagation in multilayer saturated medium [8,9,10].

Developing a robust and effective numerical method for hyperbolic initial boundary value problems has been a research hotspot in computational mechanics. Effective representation of acoustic propagation in multilayer pavement with irregular interfaces is challenging, and simulation of this problem is also a motivation for our study. The time domain and space domain of sound wave propagation in multilayer pavement are discontinuous, which makes the traditional time continuous method no longer applicable.

In the traditional time-continuous Galerkin method represented by the Newmark method, time discontinuity is not allowed. In its general implementation, the discontinuity exists in the time domain, the space domain is not described accurately, and spurious numerical oscillations always appear in the wave-after stage. In the time discontinuous Galerkin finite element scheme, the time interval is further subdivided into independent time plates and time discontinuities between plates are allowed [11]. It has been shown that time discontinuous Galerkin finite element schemes have superior stability, robustness, and higher-order accuracy, and the schemes have evolved considerably in the last decades [12]. However, high-order discontinuous Galerkin methods always have high order polynomials inside each mesh element and high cost per iteration. Li et al. presented a novel time-discontinuous Galerkin finite element method [13,14], which has lower order polynomials and lower cost, for solving wave propagation problems with weak numerical oscillations generated at the interfaces between different layers. Based on the above studies, we propose an improved time-discontinuous Galerkin finite element method for simulating fluid–structure coupled problems, using an artificial damping scheme to reduce the numerical oscillations occurring at the wavefront and interlayer interfaces, with a particular focus on the solution of the impact problem of multi-layer saturated pavements [15,16].

In this work, we describe the acoustic wave propagation problem of multilayered pavement of water saturated by using the generalized Biot’s theory and extending the new version numerical algorithm to solve the equations. The results of three examples show that the present MDGFEM can effectively filter out the spurious numerical oscillations generated at the wave-front stage and interfaces between pavement layers. MDGFEM provides a more accurate solution than CGFEM (such as the Newmark method) under the premise of using the same time step size for the acoustic wave propagation problem of multilayered pavement of water saturated and subjected to impulsive loading.

2. Biot’s Theory of Pavement and Spatial Discretization

This section summarizes the fluid–structure coupling governing equations of the linearized model of saturated porous pavement subjected to impulsive loading. The equation of motion for the solid skeleton is expressed by [10,13]:

$${\sigma }_{ij,j}+\rho b_i-\rho {\ddot{u}}_i+{\rho }_{\mathit{f}}\left({\ddot{w}}_i+{\dot{w}}_k{\dot{w}}_{i,k}\right)=0$$

(1)

The equation of motion of pore fluid can be written as

$$-p_{,i}+{\rho }_{\mathit{f}}{b}_i={\rho }_{\mathit{f}}\left({\ddot{u}}_i+\frac{1}{n}\left({\ddot{w}}_i+{\dot{w}}_k{\dot{w}}_{i,k}\right)\right)+\frac{R_i}{n}=0$$

(2)

In which, $\frac{R_i}{n}=k_{w,ij}^{-1}{\dot{w}}_j$; ${\sigma }_{ij}$ is the total Cauchy stress; $\rho$ is the assembly density; ${\rho }_{\mathit{f}}$ is the fluid density; $b_i$ is the acceleration of body force; $w_k$ is the fluid displacement; $u_i$ is the solid displacement; and $n$ is the porosity.

The equation of solid skeleton equilibrium can be obtained by inserting $U_i=u_i+\frac{1}{n}{w}_i$ into Equations (1) and (2) and subtracting Equation (2) $\times n$ from Equation (1).

$${\sigma }_{ij,j}^{″}-\left(\alpha -n\right)p_{,i}+\left(1-n\right){\rho }_s{b}_i-\left(1-n\right){\rho }_s{\ddot{u}}_i-R_i=0$$

(3)

In which, ${\sigma }_{ij}^{″}={\sigma }_{ij}-\alpha p$ is the effective stress tensor in generalized Biot theory; $R_i=n^2{k}_{w,ij}^{-1}\left({\dot{U}}_j-{\dot{u}}_j\right)$. Since both solid and water are isotropic materials, the Biot parameter [10] $\alpha =1-\frac{K_T}{K_s}$, where $K_T=\frac{E}{3\left(1-2\nu \right)}$. Therefore, Equation (2) can be rewritten as

$$-n{p}_{,i}+n{\rho }_f{b}_i-n{\rho }_f{\ddot{U}}_i-R_i=0$$

(4)

Then the mass conservation equation of fluid phase flow is expressed as

$$-n{\dot{U}}_{i,i}=\left(\alpha -n\right){\dot{\epsilon }}_{ii}+\frac{1}{Q}\dot{p}$$

(5)

From Equation (5), the equations of motion for fluids and porous solids can be written as

$${\sigma }_{ij,j}^{″}-\left(\alpha -n\right)Q{\left[\left(\alpha -n\right)u_{k,k}+n{U}_{k,k}\right]}_{,i}-\left(1-n\right){\rho }_s{\ddot{u}}_i-n^2{k}_{w,ij}^{-1}\left({\dot{U}}_j-{\dot{u}}_j\right)=0$$

(6)

$$-nQ{\left[\left(\alpha -n\right)u_{k,k}+n{U}_{k,k}\right]}_{,i}-n{\rho }_f{\ddot{U}}_i-n^2{k}_{w,ij}^{-1}\left({\dot{U}}_j-{\dot{u}}_j\right)=0$$

(7)

In which, $\frac{1}{Q}=\frac{n}{K_f}+\frac{\alpha -n}{K_s}$; where $K_f$ and $K_s$ are the bulk module of fluid and solid, respectively.

After semi-discretization of Equations (6) and (7), the following equation can be obtained.

$$\left[\begin{array}{cc}{\mathit{M}}_s& 0\\ 0& {\mathit{M}}_f\end{array}\right]\left[\begin{array}{c}\ddot{\overline{\mathit{u}}}\left(t\right)\\ \ddot{\overline{\mathit{U}}}\left(t\right)\end{array}\right]+\left[\begin{array}{cc}{\mathit{C}}_{11}& {\mathit{C}}_{12}\\ {\mathit{C}}_{21}& {\mathit{C}}_{22}\end{array}\right]\left[\begin{array}{c}\dot{\overline{\mathit{u}}}\left(t\right)\\ \dot{\overline{\mathit{U}}}\left(t\right)\end{array}\right]+\left[\begin{array}{cc}{\mathit{K}}_{11}& {\mathit{K}}_{12}\\ {\mathit{K}}_{21}& {\mathit{K}}_{22}\end{array}\right]\left[\begin{array}{c}\overline{\mathit{u}}\left(t\right)\\ \overline{\mathit{U}}\left(t\right)\end{array}\right]=\left[\begin{array}{c}{\mathit{f}}_u\\ {\mathit{f}}_U\end{array}\right], t\in I=\left(0,T\right)$$

(8)

In which,

$$\begin{gathered} {\mathit{M}}_s={\int }_{\mathsf{\Omega }}{N}_K^u\left(1-n\right){\rho }_s{\delta }_{ij}{N}_L^{u}d\mathsf{\Omega }; {\mathit{M}}_f={\int }_{\mathsf{\Omega }}{N}_K^U\left(1-n\right){\rho }_f{\delta }_{ij}{N}_L^{U}d\mathsf{\Omega } \\ {\mathit{C}}_{11}={\int }_{\mathsf{\Omega }}{N}_K^u{n}^2{\delta }_{ij}{k}_{w,ij}^{-1}{N}_L^{u}d\mathsf{\Omega }; {\mathit{C}}_{12}={\int }_{\mathsf{\Omega }}{N}_K^u{n}^2{\delta }_{ij}{k}_{w,ij}^{-1}{N}_L^{U}d\mathsf{\Omega } \\ {\mathit{C}}_{22}={\int }_{\mathsf{\Omega }}{N}_K^U{n}^2{\delta }_{ij}{k}_{w,ij}^{-1}{N}_L^{U}d\mathsf{\Omega }; {\mathit{K}}_{11}={\mathit{K}}_{\mathrm{t}}+{\mathit{K}}_{11}^{*} \\ {\mathit{K}}_t={\int }_{\mathsf{\Omega }}{N}_{\mathit{K},k}^u{\left(\alpha -n\right)}^{2}Q{N}_{L,j}^{u}d\mathsf{\Omega }; {\mathit{K}}_{11}^{*}={\int }_{\mathsf{\Omega }}{N}_{\mathit{K},i}^u{D}_{ikjl}{N}_{L,l}^{u}d\mathsf{\Omega } \\ {\mathit{K}}_{12}={\int }_{\mathsf{\Omega }}{N}_{\mathit{K},i}^{u}n\left(\alpha -n\right)Q{N}_{L,j}^{U}d\mathsf{\Omega }; {\mathit{K}}_{22}={\int }_{\mathsf{\Omega }}{N}_{\mathit{K},i}^U{n}^{2}Q{N}_{L,j}^{U}d\mathsf{\Omega } \end{gathered}$$

(9)

In which, $\overline{\mathit{u}}\left(t\right)$ is the nodal displacement of the solid skeleton; and $\overline{\mathit{U}}\left(t\right)$ is the nodal intrinsic displacement of the pore fluid. $D_{ikjl}$ is the fourth-order tensor used to define the constitutive law of soil skeleton. $N_{K,i}^u$ and $N_{K,i}^U$ are the spatial derivative of the kth node shape function concerning the i-axis of the spatial finite element discreteness.

It is worth noting that we need to set appropriate boundary conditions, which need to satisfy Equations (6) and (7). At the same time, the displacement of the points on two surfaces ${\Gamma }_u$ and ${\Gamma }_U$ needs to be satisfied:

$$\mathit{u}={\mathit{u}}_0\left(t\right) \mathrm{on} {\Gamma }_u$$

(10)

$$\mathit{U}={\mathit{U}}_0\left(t\right) \mathrm{on} {\Gamma }_U$$

(11)

In which, ${\mathit{U}}_0$ and ${\mathit{u}}_0$ represent the specified values of the displacement vectors for the two phases—fluid and solid, respectively.

If surface loads act on the above two surfaces, on the other hand, the boundary conditions should satisfy the following equations.

$${\mathit{f}}_u={\mathit{f}}_{u0}\left(t\right) \mathrm{on} {\Gamma }_u$$

(12)

$${\mathit{f}}_U={\mathit{f}}_{U0}\left(t\right) \mathrm{on} {\Gamma }_U$$

(13)

In which, ${\mathit{f}}_{u0}\left(t\right)$ and ${\mathit{f}}_{U0}\left(t\right)$ are the given surface loads.

The surface impact load used in this paper is expressed as a piecewise function.

$${\mathit{f}}_{u0}\left(t\right)={\mathit{f}}_{U0}\left(t\right)=\{\begin{array}{c}P tϵ\left[0,t_0\right]\\ 0 t>t_0]\end{array}$$

(14)

In which, $P$ is the amplitude of load, and $t_0$ is the constant time.

3. Temporal Discretization and Time Discontinuous Galerkin Finite Element Method

3.1. Time Discontinuous Galerkin Finite Element Method

The final matrix form of the fluid–solid interaction problem can be rewritten as

$$\mathit{M}\ddot{\mathit{d}}\left(t\right)+\mathit{C}\dot{\mathit{d}}\left(t\right)+\mathit{K}\mathit{d}\left(t\right)=\mathit{f}\left(t\right), t\in I=\left(0,T\right)$$

(15)

In which,

$$\begin{gathered} \mathit{M}=\left[\begin{array}{cc}{\mathit{M}}_s& 0\\ 0& {\mathit{M}}_f\end{array}\right]; \mathit{C}=\left[\begin{array}{cc}{\mathit{C}}_{11}& {\mathit{C}}_{12}\\ {\mathit{C}}_{21}& {\mathit{C}}_{22}\end{array}\right]; \mathit{K}=\left[\begin{array}{cc}{\mathit{K}}_{11}& {\mathit{K}}_{12}\\ {\mathit{K}}_{21}& {\mathit{K}}_{22}\end{array}\right]; \mathit{f}\left(t\right)=\left[\begin{array}{c}{\mathit{f}}_u\\ {\mathit{f}}_U\end{array}\right]; \\ \ddot{\mathit{d}}\left(t\right)=\left[\begin{array}{c}\ddot{\overline{\mathit{u}}}\left(t\right)\\ \ddot{\overline{\mathit{U}}}\left(t\right)\end{array}\right]; \dot{\mathit{d}}\left(t\right)=\left[\begin{array}{c}\dot{\overline{\mathit{u}}}\left(t\right)\\ \dot{\overline{\mathit{U}}}\left(t\right)\end{array}\right]; \mathit{d}\left(t\right)=\left[\begin{array}{c}\overline{\mathit{u}}\left(t\right)\\ \overline{\mathit{U}}\left(t\right)\end{array}\right]; \end{gathered}$$

(16)

The difference between DGFEM and continuous Galerkin method is that DGFEM allows discontinuities of functions in discrete time series, $0<t_n<t_N.$ The temporal jump of the discontinuous functions at $t_n$ is defined as

$$⟦{\mathit{d}}_n\left(t\right)⟧ = \mathit{d}\left(t_n^{+}\right)-{\mathit{d}}_n\left(t_n^{-}\right)$$

(17)

In which, $\mathit{d}\left(t_n^{+}\right)$ and ${\mathit{d}}_n\left(t_n^{-}\right)$ are discontinuities of functions at $t_n$, where $w\left(t_n^{\pm }\right)=\underset{\epsilon \to 0^{\pm }}{\lim }w\left(t_n+\epsilon \right)$.

For each time sub-domain $I_n=\left(t_n^{+},t_{n+1}^{-}\right)$, we define a time step length $\mathsf{\Delta }t=t_{n+1}-t_n$, for $n=0,1\cdots N-1$. For each time sub-domain $I_n$, we use the following third-order Hermite (P3) time shape function to interpolate to obtain the global nodal temperature displacement vector.

$$\mathit{d}\left(t\right)={\mathit{d}}_n^{+}{N}_{t1}\left(t\right)+{\mathit{d}}_{n+1}^{-}{N}_{t2}\left(t\right)+{\mathit{v}}_n^{+}{M}_{t1}\left(t\right)+{\mathit{v}}_{n+1}^{-}{M}_{t2}\left(t\right)$$

(18)

In which, ${\mathit{d}}_n^{+}$ and ${\mathit{d}}_{n+1}^{-}$ are the global nodal values of displacements vectors at times $t_n \mathrm{and} t_{n+1}$; ${\mathit{v}}_n^{+}$ and ${\mathit{v}}_{n+1}^{-}$ are the time-derivatives, respectively. The global nodal values of displacements vector and its time-derivative, i.e., ${\mathit{d}}_n^{-}$ and ${\mathit{v}}_n^{-}$ at time $t_n^{-}$ can be obtained at the end of the previous time step.

Dropping the superscripts of the vectors ${\mathit{d}}_n^{+}$, ${\mathit{d}}_{n+1}^{-}$, ${\mathit{v}}_n^{+}$, ${\mathit{v}}_{n+1}^{-}$, and the time variable $t$, Equation (10) can be rewritten as

$$\mathit{d}\left(\mathit{t}\right)={\mathit{d}}_n{N}_{t1}\left(t\right)+{\mathit{d}}_{n+1}{N}_{t2}\left(t\right)+{\mathit{v}}_n{M}_{t1}\left(t\right)+{\mathit{v}}_{n+1}{M}_{t2}\left(t\right)$$

(19)

$$\mathit{v}\left(t\right)={\mathit{v}}_n^{+}{\mathsf{\lambda }}_1\left(t\right)+{\mathit{v}}_{n+1}^{-}{\mathsf{\lambda }}_2\left(t\right)$$

(20)

In which,

$$\begin{gathered} N_{t1}\left(t\right)={\lambda }_1^2\left({\lambda }_1+3{\lambda }_2\right), N_{t2}\left(t\right)={\lambda }_2^2\left({\lambda }_2+3{\lambda }_1\right), M_{t1}\left(t\right)={\lambda }_1^2{\lambda }_2\mathsf{\Delta }t, M_{t2}\left(t\right)=-{\lambda }_2^2{\lambda }_1\mathsf{\Delta }t \\ {\lambda }_1=\frac{t_{n+1}-t_n}{\mathsf{\Delta }t}, {\lambda }_2=\frac{t_n-t_{n+1}}{\mathsf{\Delta }t} \end{gathered}$$

(21)

The global nodal value ${\mathit{v}}_n$ of the time derivative of the displacement vector at any time $(t_n,t_{n+1})$ is transformed into the independent variable after interpolation by a linear (P1) time shape function as Equation (20).

The weak form of the semi-discrete equation and the constraints of the system $\dot{\mathit{d}}-\mathit{v}=\mathbf{0}$, as well as the discontinuity conditions $\mathit{d},\mathit{v}$ for the typical time subdomain can be expressed as

$$\underset{I_n}{\overset{}{{\displaystyle \int }}}\mathbf{\delta }{\mathit{v}}^{\mathbf{T}}\left(\mathit{M}\dot{\mathit{v}}+\mathit{C}\mathit{v}+\mathit{K}\mathit{d}-\mathit{Q}\right)\mathit{d}\mathit{t}+\underset{I_n}{\overset{}{{\displaystyle \int }}}\mathbf{\delta }{\mathit{d}}^T\mathit{K}\left(\dot{\mathit{d}}-\mathit{v}\right)\mathit{d}\mathit{t}+\mathbf{\delta }{\mathit{d}}_{\mathit{n}}^{\mathit{T}} \mathit{K}⟦{\mathit{d}}_{\mathit{n}}⟧+\mathbf{\delta }{\mathit{v}}_{\mathit{n}}^{\mathit{T}} \mathit{M}⟦{\mathit{v}}_{\mathit{n}}⟧=\mathbf{0}$$

(22)

Substituting Equations (11) and (13) into Equation (14) yields the following DGFEM matrix equation for solving from the independent changes.

$${\mathit{d}}_n^{+}={\mathit{d}}_n^{-}$$

(23)

$$\left[\begin{array}{cc}\mathit{M}+\frac{\mathsf{\Delta }t}{6}\mathit{C}-\frac{\mathsf{\Delta }{t}^2}{12}\mathit{K}& -\frac{\mathsf{\Delta }t}{6}\mathit{C}-\frac{\mathsf{\Delta }{t}^2}{12}\mathit{K}\\ \frac{\mathsf{\Delta }t}{2}\mathit{C}+\frac{\mathsf{\Delta }{t}^2}{3}\mathit{K}& \mathit{M}+\frac{\mathsf{\Delta }t}{2}\mathit{C}+\frac{\mathsf{\Delta }{t}^2}{6}\mathit{K}\end{array}\right]\left\{\begin{array}{c}{\mathit{v}}_n\\ {\mathit{v}}_{n+1}\end{array}\right\}=\left\{\begin{array}{c}{\mathit{f}}_1^e-{\mathit{f}}_2^e+\mathit{M}{\mathit{v}}_n^{-}\\ {\mathit{f}}_1^e+{\mathit{f}}_2^e+\mathit{M}{\mathit{v}}_n^{-}-\Delta t\mathit{K}{\mathit{d}}_n^{-}\end{array}\right\}$$

(24)

$${\mathit{d}}_{n+1}^{-}={\mathit{d}}_n^{-}+\frac{1}{2}\mathsf{\Delta }\mathrm{t}\left({\mathit{v}}_n+{\mathit{v}}_{n+1}\right)$$

(25)

In which, ${\mathit{f}}_1^e=\frac{\mathsf{\Delta }t}{3}{\mathit{f}}_n^e+\frac{\mathsf{\Delta }t}{6}{\mathit{f}}_{n+1}^e$, ${\mathit{f}}_2^e=\frac{\mathsf{\Delta }t}{6}{\mathit{f}}_n^e+\frac{\mathsf{\Delta }t}{3}{\mathit{f}}_{n+1}^e$.

The above formula shows that the continuity of the node displacement vector $\left\{\mathit{d}\right\}$ can be automatically satisfied at any time. While the nodal derivative vector $\left\{\mathit{v}\right\}$ at discrete time level is discontinuous.

3.2. Artificial Damping Scheme for DGFEM

As mentioned in articles [14,15,16], if high-frequency and low-frequency oscillations are to be reduced, stiffness proportional damping coefficient and mass proportional damping coefficient are required respectively. Therefore, we introduce an artificial stiffness proportional Rayleigh-type damping scheme in DGFEM to facilitate filtering out the wavefront oscillations caused by high-frequency shock loads, which takes the form:

$$\left[{\mathit{C}}_{\lambda }\right]={\beta }^{\left(i\right)}\mathit{K}$$

(26)

$${\beta }^{\left(i\right)}=2{\xi }^{\left(i\right)}/{\omega }_n$$

(27)

In which, $\left[{\mathit{C}}_{\lambda }\right]$ is the introduced stiffness proportional damping matrix, ${\beta }^{\left(i\right)}$ is the damping constant, ${\xi }^{\left(i\right)}$ is the damping ratio estimated from the overall behavior of the pavement structure, and ${\omega }_n$ is the fundamental frequency of the multilayer pavement system. Summarizing the results of several trial runs, it can also be determined that the value of the effective stiffness proportional damping constant ${\beta }^{\left(i\right)}$ depends on the step size, that is, ${\beta }^{\left(i\right)}$ cannot be greater than the time step size.

After substituting the stiffness proportional artificial damping matrix Equation (26) into the DGFEM matrix Equation (24), we can get:

$$\left[\begin{array}{cc}\mathit{M}+\frac{\mathsf{\Delta }t}{6}{\mathit{C}}_{\mathit{w}}-\frac{\mathsf{\Delta }{t}^2}{12}\mathit{K}& -\frac{\mathsf{\Delta }t}{6}{\mathit{C}}_{\mathit{w}}-\frac{\mathsf{\Delta }{t}^2}{12}\mathit{K}\\ \frac{\mathsf{\Delta }t}{2}{\mathit{C}}_{\mathit{w}}+\frac{\mathsf{\Delta }{t}^2}{3}\mathit{K}& \mathit{M}+\frac{\mathsf{\Delta }t}{2}{\mathit{C}}_{\mathit{w}}+\frac{\mathsf{\Delta }{t}^2}{6}\mathit{K}\end{array}\right]\left\{\begin{array}{c}{\mathit{v}}_n\\ {\mathit{v}}_{n+1}\end{array}\right\}=\left\{\begin{array}{c}{\mathit{f}}_1^e-{\mathit{f}}_2^e+\mathit{M}{\mathit{v}}_n^{-}\\ {\mathit{f}}_1^e+{\mathit{f}}_2^e+\mathit{M}{\mathit{v}}_n^{-}-\mathsf{\Delta }t\mathit{K}{\mathit{d}}_n^{-}\end{array}\right\}$$

(28)

In which,

$$\left[{\mathit{C}}_{\mathit{w}}\right]=\left[\mathit{C}\right]+\left[{\mathit{C}}_{\lambda }\right]$$

(29)

Equations (23), (25) and (28) form the final expressions for the present MDGFEM.

4. Numerical Results and Discussion

To demonstrate the superiority of MDGFEM over traditional time-continuous Galerkin methods, we investigate numerical examples of a multi-layer saturated pavement subjected to impact loads in this section. One-dimensional and two-dimensional numerical results show that the proposed MDGFEM formulation can be more accurate for such problems. The basic parameters of the pavement used in each calculation example are shown in

Table 1. Values of relevant parameters for calculations.

Parameters	Cement Concrete Pavement	Semi-Rigid Base	Natural Soil
Darcy permeability coefficient $k$ (m/s)	1.0 × 10−6	1.0 × 10−5	1.0 × 10−6
Biot parameter $\alpha$	1.0	1.0	1.0
Porosity $n$	0.05	0.04	0.4
Young’s module (GPa)	25	2.1	1.0
shear module of the solid (GPa)	10.42	0.84	0.38
$\mathrm{bulk} \mathrm{module} \mathrm{of} \mathrm{the} \mathrm{fluid} K_f$ (GPa)	0.02	0.02	0.02
$\mathrm{bulk} \mathrm{module} \mathrm{of} \mathrm{the} \mathrm{solid} K_s$ (GPa)	13.89	1.40	0.83
$\mathrm{densities} \mathrm{of} \mathrm{the} \mathrm{fluid} {\rho }_f$ (kg/m3)	1000	1000	1000
$\mathrm{densities} \mathrm{of} \mathrm{the} \mathrm{solid} {\rho }_s$ (kg/m3)	2400	2100	1000
Poisson’s ratios	0.20	0.25	0.30

.

4.1. 1-D Elastic Wave Propagation Problem

The effectiveness of the DGFEM to solve the problem of generalized Biot’s theory was verified in articles [13,16]. In order to prove the good performance of MDGFEM in cement concrete medium, we first consider the propagation of one-dimensional stress waves in a cement-concrete column with a length of 0.9 m and a cross-sectional area of 1 m². As shown in Figure 1, one end of the column is fixed, and the other end is subjected to the axial force impulse F(t), the magnitude of which is

$$F\left(t\right)=\{\begin{array}{c}1.0\times {10}^{3}N \left(0\le t\le 45 \mathsf{\mu }\mathrm{s}\right)\\ 0 (t>45 \mathsf{\mu }\mathrm{s})\end{array}$$

(30)

The cement concrete column is divided into 90 elements along the length, that is, each element is 0.01 m. Figure 2 presents the stress solutions at difference time levels of $t_1=60 \mathsf{\mu }\mathrm{s}$, $t_2=180 \mathsf{\mu }\mathrm{s}$ produced by Newmark’s method, and MDGFEM with time step size $\Delta t=1.0\times {10}^{-6} \mathrm{s}$. Compared with the polyline obtained by the Newmark method, the smooth curve obtained by MDGFEM shows the superiority of the algorithm in filtering out severe numerical oscillations in wave after wave.

4.2. 1-D Multi-Layer Acoustic Wave Propagation Problem in Multilayered Pavement

Secondly, we build a one-dimensional model of the multi-layer pavement structure to study the wave propagation problem. The lengths of layers are 0.25 m (Cement concrete layer), 0.2 m (Cement stabilized subgrade), 0.45 m (Natural soil). The cross-sectional area of the multi-layer pavement structure is still 1 m². As shown in Figure 3, one end of the column is fixed, and the other end is subjected to the axial force impulse F(t), the magnitude of which is

$$F\left(t\right)=\{\begin{array}{c}2.0\times {10}^{4}N \left(0\le t\le 40 \mathsf{\mu }\mathrm{s}\right)\\ 0 (t>40 \mathsf{\mu }\mathrm{s})\end{array}$$

(31)

The cement concrete column is divided into 45 elements along the length, that is, each element is 0.02 m. Figure 4 and Figure 5 give the wave propagation solutions of compression and pressure generation using Newmark method and MDGFEM with a time step size $\Delta t=4.0\times {10}^{-6} \mathrm{s}$ by taking two different time levels of the compression and pressure, resulting in wave propagation solutions at different time levels of $t_1=80 \mathsf{\mu }\mathrm{s}$ and $t_2=200 \mathsf{\mu }\mathrm{s}$ for the example. It is remarked that the stress and pressure time-history curves obtained by the Newmark method have obvious spurious numerical oscillations in the pavement interlayer and wave-front regions, while the smooth and continuous stress-pressure curve of MDGFEM proves that the artificial damping added provides a more accurate solution to the wave propagation problem in multilayer pavement.

4.3. 2-D Multi-Layer Elastic Wave Propagation Problem in Multilayered Pavement

As the third example, we model a 2D saturated porous pavement structure to study the propagation of elastic waves. As shown in Figure 6, a square area (0.90 m × 0.90 m) with holes (0.020 m × 0.20 m) is subjected to the impulsed inclined load at the upper left corner. The domain consists of 46 × 46 nodes and is analyzed as a plane strain problem. The thicknesses of layers are 0.25 m (Cement concrete layer), 0.2 m (Cement stabilized subgrade), and 0.45 m (Natural soil). The force impulse $F\left(t\right)$ is expressed as

$$F\left(t\right)=\{\begin{array}{c}2.0\times {10}^{4}N \left(0\le t\le 40 \mathsf{\mu }\mathrm{s}\right)\\ 0 (t>40 \mathsf{\mu }\mathrm{s})\end{array}$$

(32)

Figure 7 and Figure 8 shows the cloud diagram of fluid pressure and solid stress distribution in the pores of the multilayer pavement structure at time $t_1=0.15 \mathsf{\mu }\mathrm{s}$, $t_2=0.3 \mathsf{\mu }\mathrm{s}$, and $t_3=0.4 \mathsf{\mu }\mathrm{s}$ obtained using the MDGFEM method and the Newmark method with the time step $\Delta t=4.0\times {10}^{-6} \mathrm{s}$. After an impact load is applied in the upper left corner of this domain, the wave propagates forward as time increases. Pressure and stress nephograms at different times demonstrate the obvious advantage of MDGFEM in reproducing the silent properties in the area around the hole and at the layer interface compared to the Newmark method, which is helpful for the further development and improvement of nondestructive testing technology of cracks and bottom void diseases in cement concrete pavement.

5. Conclusions

In this article, we present the generalized Biot’s theory formulation and a new numerical treatment version of the generalized fluid–structure coupled theory of water-saturated multilayered pavement subjected to high-frequency impact loads. The results presented in this article are demonstrated by a number of numerical examples in both one and two dimensions. The advantages of the MDGFEM proposed in this paper are indicated by a comparison with traditional time-continuous methods such as Newmark method. The present method is valid to filter out the spurious numerical oscillations in wave-front and wave-after stage, boundaries of hole, and the interface between the layers successfully for the propagation of acoustic waves in multilayered pavement.