# A Time-Discontinuous Galerkin Finite Element Method for the Solution of Impact Problem of Gas-Saturated Coal

Authors: Jingfei Zhang; Deyong Guo; Wenhua Wu; Pan Guo
Journal: Shock and Vibration
 (2020)

---

## Abstract

Based on the general Biot theory of saturated porous media, a modified time-discontinuous Galerkin finite element method (MDGFEM) is presented to simulate the structural dynamics and wave propagation problems of gas-saturated coal subjected to impact loading. Numerical results of one dimension and two dimensions show that the present MDGFEM possesses better abilities and provides much more accurate solutions than the traditional Newmark method and previous DGFEM for the impact problem. It can effectively capture the discontinuities of the wave and filter out the effects of spurious numerical oscillation induced by high-frequency impulsive load. The results can provide a technological basis for the research of the prevention of coal and gas dynamic disasters under deep mining. And the method could be useful for the further numerical research of coal-rock-gas coupling problems and coal-gas-heat coupling problems subjected to impact loading.


---

## Body

## 1. Introduction

With the increase of coal mining depth, coal and gas outburst disasters are becoming more and more severe [1]. The mechanism of outburst accidents, which has been studied by many experts and scholars, demonstrates that the coupling scheme of gas pressure and coal structure plays a dominant role in the coal and gas outburst disasters, and the outburst accidents are always caused by different impact loading [1, 2]. However, understanding of the disaster caused by the impact load under complicated environment is far from enough and the risk of outburst accidents cannot be completely eliminated [1–4]. And one motivation of our research is to investigate the dynamic and wave propagation problem of gas-saturated coal subjected to impact loads.In the past, efforts have been made to find a porous media model, such as Biot’s theory, to discover the reflection and transmission of waves in a porous medium [5, 6]. Biot’s theory derives the equations of motion of each phase based on energy considerations which include the inertial, potential, and viscous coupling between the two phases. The equations governing the interaction of the solid and fluid media for dynamics phenomena were first established in 1956 by Biot [5]. Owing to its simplicity, many authors have followed Biot’s approach to research the problem of wave propagation in gas-saturated coal [3, 4].Developing a robust and efficient numerical method for the coupled problems subjected to the high-frequency impact loading is challenging, and the simulation of the problems is another motivation of our research. This is particularly the case for the problems of wave propagation in gas-saturated coal where discontinuities exist both in the time domain and the space domain. Despite the popularity, the traditional time integration method (such as Newmark method) [7–9] struggles to provide more accurate solutions and eliminate spurious numerical oscillations in numerical simulation of wave propagation involving sharp gradients and discontinuities in the time domain.Li et al. presented a novel time-discontinuous Galerkin finite element method (DGFEM) for solving dynamics and wave propagation in nonlinear solids and saturated porous media and heat wave propagation problem subjected to impact loading [10–12]. The remarkable characteristics of the DG method is that the basic unknown variables together with their temporal derivatives are assumed to be discontinuous at a defined time point and are assumed to be interpolated by third-order Hermite function and linear function in a time step, respectively. We then successfully applied the MDGFEM, based on an additional artificial damping scheme, to simulate heat wave propagation and generalized thermoelastic problems subjected to thermal shock [13, 14]. Numerical results demonstrate that the MDGFEM illustrates better performance in numerical simulation of wave propagation in eliminating spurious numerical oscillations and in providing more accurate solutions than that of the traditional time integration method and the DGFEM before being modified.For the high-speed motion and high pressure of the gas-saturated coal problem under high-frequency loading, the acceleration of porous fluid should be considered [6]. In this article, an additional artificial damping discontinuous Galerkin numerical algorithm for gas-saturated coal problems is developed on the basis of elastic-plastic saturated porous medium model [6] and previous research studies [12]. An additional stiffness proportional damping scheme is brought into the final finite element discretized form to reduce the numerical oscillations appearing in the wave-front stage of the elastic-plastic stress wave and the pressure wave. Numerical results show clearly that the MDGFEM method is valid for filtering out spurious oscillations in both waves and for providing much more accurate solutions than those obtained using the previous DGFEM and the traditional time integration methods (e.g., the Newmark method). The present modeling method provides accurate numerical method for the prediction and early warning of coal and gas outburst and could also be extended to the study of coal-rock-gas coupling problems and coal-gas-heat coupling problems subjected to impact loading.


## 2. The Basic Dynamics Governing Equations and DGFEM Formulation for Coal and Rock Porous Media

This section summarizes the basic governing equations of gas-saturated coal porous media subjected to impact loading. As the porous fluid is compressible, the equations of motion for the coal skeleton can be expressed as [6](1)σij,j+ρbi−ρu¨i+ρfw¨i+w˙kw˙i,k=0.For the pore fluid, the equation of motion can be written as(2)−p,i+ρfbi=ρfu¨i+1nw¨i+w˙kw˙i,k+Rin=0,where Ri/n=kw,ij−1w˙j; σij is the total Cauchy stress; bi is the body force acceleration; ρ, ρf, and ρs are the densities of the assembly, the fluid and the solid phases, respectively; wk is the displacement of the fluid phase; ui is the displacement of the solid phase; and n is the porosity.Insertion ofUi=ui+1/nwi into equations (1) and (2) and subtraction of equation (2) ×n from equation (1) leads to the equation of solid skeleton equilibrium.(3)σij,j″−α−np,i+1−nρsbi−1−nρsu¨i−Ri=0,where σij″=σij−αp is the generalized Biot effective stress tensor and Ri=n2kw,ij−1U˙j−u˙j. For an isotropic material, the Biot parameter α=1−KT/Ks, where KT=E/31−2v.Equation (2) can be written as(4)−np,i+nρfbi−nρfU¨i−Ri=0.Then, the mass conservation equation of the fluid flow can be written as(5)−nU˙i,i=α−nε˙ii+1Qp˙.Using the definition of equation (5), equations (1) and (2) can rewritten as(6)σij,j″−α−nQα−nuk,k+nUk,ki−1−nρsu¨−n2kw,ij−1U˙j−u˙j=0,(7)−nQα−nuk,k+nUk,k,i−nρfU¨i−n2kw,ij−1U˙j−u˙j=0,where 1/Q=n/Kf+α−n/Ks; Kf is the bulk modulus of the fluid; and Ks is the bulk modulus of the solid. The linearized form of the semidiscretized systems of the motion equations (6) and (7) can be written as(8)Ms00Mfu¯¨tU¯¨t+C11C12C21C22u¯˙tU¯˙t+K11K12K21K22u¯tU¯t=fufU,t∈0,T,in which(9)Ms=∫ΩNKu1−nρsδijNNudΩ,Mf=∫ΩNKU1−nρfδijNLUdΩ,C11=∫ΩNKun2δijkw,ij−1NLudΩ,C12=∫ΩNKun2δijkw,ij−1NLUdΩ,C22=∫ΩNKUn2δijkw,ij−1NLUdΩ,K11=Kt+K11∗,Kt=∫ΩNK,kuα−n2QNL,judΩ,K11∗=∫ΩNK,kuDijklQNL,iudΩ,K12=∫ΩNK,kunα−nQNL,jUdΩ,K22=∫ΩNK,kUn2QNL,jUdΩ,where u¯t and U¯t are the nodal displacements of the solid skeleton and the nodal intrinsic displacements of the pore fluid, respectively, Dijkl is a fourth-order tensor defining a constitutive law for the solid skeleton, and NK,iu and NK,iU are the spatial derivatives of the shape functions of the Kth node with respect to the i coordinate axis for the spatial finite element discretization.It should be noted that the appropriate boundary conditions associated with the governing equations (6) and (7) must be adopted. When the displacements are prescribed on the surfaces Γu and ΓU, respectively, one has(10)u=u0tonΓu,(11)U=U0tonΓU,where U0t and u0t are the prescribed values of fluid and solid displacement vector.On the other hand, if surface loadings are applied to the corresponding surfacesΓu and ΓU, the following boundary conditions should be satisfied:(12)fu=fu0tonΓu,(13)fU=fU0tonΓU,where fu0t and fU0t are the given surface loadings, respectively.In the present paper, surface impulse loading is applied with the form as(14)fu0t=fU0t=P,t∈0,t0,0,t>t0,in which P and t0 are the amplitude of loading and the constant time.


## 3. Temporal Discretization and Time-Discontinuous Galerkin Finite Element Method

### 3.1. Time-Discontinuous Galerkin Finite Element Method

The standard Galerkin discretized equation (8) can be rewritten as follows:(15)Md¨t+Cd˙t+Kdt=ft,t∈0,T,in which(16)M=Ms00Mf;C=C11C12C21C22;K=K11K12K21K22;ft=fufU;d˙t=u¯¨tU¯¨t;d˙t=u¯˙tU¯˙t;dt=u¯tU¯t.The main features of the discontinuous Galerkin integration method in time domain have been described in our previous articles [13, 15]. Different from the continuous Galerkin method, the DGFEM permits the discontinuity of functions at discrete-time sequence, 0⋯<tn<⋯tN. The temporal jump of the discontinuous function dtn at tn can be defined as(17)dtn=dtn+−dtn−,in which(18)dtn±=limε⟶0±dtn+ε.The global displacement vectors of solid and fluid at arbitrary timet in a time step In=tn+,tn+1− are interpolated by using the third-order Hermite shape functions as(19)dt=dn+N1t+dn+1−N2t+vn+M1t+vn+1−M2t,(20)N1t=λ12λ1+3λ2,N2t=λ22λ2+3λ1,M1t=λ12λ2Δt,M2t=−λ22λ1Δt,λ1=tn+1−tnΔt,λ2=tn−tn+1Δt,in which dn+ and dn+1− stand for the global nodal value of displacement vectors at times tn and tn+1 and vn+ and vn+1− stand for the global nodal value of velocity vectors, respectively. We can rewrite equation (19) by dropping the superscripts of the displacement and velocity vectors and the time variable t as follows:(21)dt=dnN1t+dn+1N2t+vnM1t+vn+1M2t.The global velocity vector of solid and fluid at arbitrary timet in a time step In=tn+,tn+1− is interpolated as an independent variable by using the linear shape functions as follows:(22)vt=vn+λ1t+vn+1−λ2t.The weak forms of the semidiscretized equation (15) and the constraint condition d˙t−v=0 along with the discontinuity conditions d and v on a typical time subdomain can be expressed as(23)∫InδvTMd¨t+Cd˙t+Kdt−ftdt+∫InδdTKd˙t−vdt+δdnTKdn+δvnTMvn=0.Substituting equations (21) and (22) into (15), we obtain the following matrix equation of DGFEM for solving from independent variations:(24)dn+=dn−,(25)M+Δt6C−Δt212K−Δt6C−Δt212KΔt2C+Δt23KM+Δt2C+Δt26Kvnvn+1=f1e−f2e+Mvn−f1e+f2e+Mvn−−ΔtKdn−,(26)dn+1=dn−+12Δtvn+vn+1,in which f1e=Δt/3fne+Δt/6fn+1e and f2e=Δt/6fne+Δt/3fn+1e.It should be noted that the global nodal displacement vector of solid and fluid is continuous at any time level and the global nodal velocity vector is still discontinuous at any time level [10–15]. To demonstrate the advantage of the present DGFEM, we have to solve several problems in our previous work [10–15]. For the elastic-plastic dynamic problems, the Newton–Raphson process is used [10, 14].


### 3.2. Modified Time-Discontinuous Galerkin Finite Element Method

When simulating problems of structure dynamics and wave propagation under impulse load, the modified time-continuous Galerkin finite element method shows better ability to filter out the effects of spurious numerical oscillations. Having demonstrated the applicability and advantage of the MDGFEM as indicated in our articles [13–15], we apply and discuss the method to solve the problem of gas-saturated coal based on Biot’s model subjected to impulse load.As we know that the selections of a stiffness proportional and mass proportional damping coefficient are effective for high-frequency oscillations and low-frequency oscillations, respectively, we modify the present DGFEM by using an artificial stiffness proportional Rayleigh-type damping scheme. The equations of the damping scheme can be written as follows:(27)Cλ=βK,(28)β=2ξωn,in which Cλ is the stiffness proportional matrix, β is the damping constant, ξ is the damping ratio estimated from the global behavior of the system, and ωn is the basic frequency of the considered structure. From equation (27), we can see that the value of stiffness proportional damping constant β is complex to calculate. As discussed in our articles [13–15], the value should be less than or equal to the time step.Using the stiffness proportional artificial damping matrix equation (27), (25) can be written in matrix form as(29)M+Δt6Cw−Δt212K−Δt6Cw−Δt212KΔt2Cw+Δt23KM+Δt2Cw+Δt26Kvnvn+1=f1e−f2e+Mvn−f1e+f2e+Mvn−−ΔtKdn−,in which(30)Cw=C+Cλ.Equations (24), (26), and (29) form the final expressions of the MDGFEM to solve the problem of gas-saturated coal based on Biot's model.


## 3.1. Time-Discontinuous Galerkin Finite Element Method

The standard Galerkin discretized equation (8) can be rewritten as follows:(15)Md¨t+Cd˙t+Kdt=ft,t∈0,T,in which(16)M=Ms00Mf;C=C11C12C21C22;K=K11K12K21K22;ft=fufU;d˙t=u¯¨tU¯¨t;d˙t=u¯˙tU¯˙t;dt=u¯tU¯t.The main features of the discontinuous Galerkin integration method in time domain have been described in our previous articles [13, 15]. Different from the continuous Galerkin method, the DGFEM permits the discontinuity of functions at discrete-time sequence, 0⋯<tn<⋯tN. The temporal jump of the discontinuous function dtn at tn can be defined as(17)dtn=dtn+−dtn−,in which(18)dtn±=limε⟶0±dtn+ε.The global displacement vectors of solid and fluid at arbitrary timet in a time step In=tn+,tn+1− are interpolated by using the third-order Hermite shape functions as(19)dt=dn+N1t+dn+1−N2t+vn+M1t+vn+1−M2t,(20)N1t=λ12λ1+3λ2,N2t=λ22λ2+3λ1,M1t=λ12λ2Δt,M2t=−λ22λ1Δt,λ1=tn+1−tnΔt,λ2=tn−tn+1Δt,in which dn+ and dn+1− stand for the global nodal value of displacement vectors at times tn and tn+1 and vn+ and vn+1− stand for the global nodal value of velocity vectors, respectively. We can rewrite equation (19) by dropping the superscripts of the displacement and velocity vectors and the time variable t as follows:(21)dt=dnN1t+dn+1N2t+vnM1t+vn+1M2t.The global velocity vector of solid and fluid at arbitrary timet in a time step In=tn+,tn+1− is interpolated as an independent variable by using the linear shape functions as follows:(22)vt=vn+λ1t+vn+1−λ2t.The weak forms of the semidiscretized equation (15) and the constraint condition d˙t−v=0 along with the discontinuity conditions d and v on a typical time subdomain can be expressed as(23)∫InδvTMd¨t+Cd˙t+Kdt−ftdt+∫InδdTKd˙t−vdt+δdnTKdn+δvnTMvn=0.Substituting equations (21) and (22) into (15), we obtain the following matrix equation of DGFEM for solving from independent variations:(24)dn+=dn−,(25)M+Δt6C−Δt212K−Δt6C−Δt212KΔt2C+Δt23KM+Δt2C+Δt26Kvnvn+1=f1e−f2e+Mvn−f1e+f2e+Mvn−−ΔtKdn−,(26)dn+1=dn−+12Δtvn+vn+1,in which f1e=Δt/3fne+Δt/6fn+1e and f2e=Δt/6fne+Δt/3fn+1e.It should be noted that the global nodal displacement vector of solid and fluid is continuous at any time level and the global nodal velocity vector is still discontinuous at any time level [10–15]. To demonstrate the advantage of the present DGFEM, we have to solve several problems in our previous work [10–15]. For the elastic-plastic dynamic problems, the Newton–Raphson process is used [10, 14].


## 3.2. Modified Time-Discontinuous Galerkin Finite Element Method

When simulating problems of structure dynamics and wave propagation under impulse load, the modified time-continuous Galerkin finite element method shows better ability to filter out the effects of spurious numerical oscillations. Having demonstrated the applicability and advantage of the MDGFEM as indicated in our articles [13–15], we apply and discuss the method to solve the problem of gas-saturated coal based on Biot’s model subjected to impulse load.As we know that the selections of a stiffness proportional and mass proportional damping coefficient are effective for high-frequency oscillations and low-frequency oscillations, respectively, we modify the present DGFEM by using an artificial stiffness proportional Rayleigh-type damping scheme. The equations of the damping scheme can be written as follows:(27)Cλ=βK,(28)β=2ξωn,in which Cλ is the stiffness proportional matrix, β is the damping constant, ξ is the damping ratio estimated from the global behavior of the system, and ωn is the basic frequency of the considered structure. From equation (27), we can see that the value of stiffness proportional damping constant β is complex to calculate. As discussed in our articles [13–15], the value should be less than or equal to the time step.Using the stiffness proportional artificial damping matrix equation (27), (25) can be written in matrix form as(29)M+Δt6Cw−Δt212K−Δt6Cw−Δt212KΔt2Cw+Δt23KM+Δt2Cw+Δt26Kvnvn+1=f1e−f2e+Mvn−f1e+f2e+Mvn−−ΔtKdn−,in which(30)Cw=C+Cλ.Equations (24), (26), and (29) form the final expressions of the MDGFEM to solve the problem of gas-saturated coal based on Biot's model.


## 4. Numerical Results and Discussion

In this section, three numerical examples of gas-saturated coal based on Biot’s model subjected to impulse load are investigated. Various results in 1-D and 2-D are presented to demonstrate that the MDGFEM formulations are capable of producing reliable results in the problem [14]. Table 1 gives the basic parameters of gas-saturated coal.Table 1
Values of relevant parameters for calculations.
ParametersCortical boneHydraulic permeabilitykw1.0×10−6Biot parameterα0.78Porosityn0.08Young’s modulusE (GPa)3.0Bulk modulus of the fluidKf (GPa)1.01×10−4Bulk modulus of the solidKs (GPa)3.33Densities of the fluidρf(kgm−30.716Densities of the solidρs(kgm−31500Poisson’s ratio0.25


### 4.1. 1-D Elastic Wave Propagation Problem of Gas-Saturated Coal Rock

We first consider a one-dimensional stress wave propagation problem in saturated gas-saturated coal. The length of the column is 200 m, and the area of the cross section is 1 m2. The column, as shown in Figure 1, with one end fixed and one end loaded by an axial force impulse can be expressed by(31)Ft=2000N,0≤t≤0.041s,0N,0.041s<t.Figure 1
Model of example 1.The column is meshed in elements of 0.4 m along the column axis. Figures2 and 3 give the resulting stress wave and pore pressure wave propagation solutions at different time levels of 0.041s and 0.08s for the example using Newmark method, DGFEM, and the MDGFEM with time step size Δt=0.001s. It is noted that the results obtained by different methods are different from each other. As compared with the Newmark method, both DGFEMs filter out the serious numerical oscillation in the wave after. And the modified DGFEM shows better ability in filtering out the serious numerical oscillation in wave than DGFEM.Figure 2
Stress profiles of different methods att = 0.041 s and t = 0.08 s.Figure 3
Gas pressure profiles of different methods att = 0.041 s and t = 0.08 s.


### 4.2. 1-D Elastic-Plastic Stress Wave Propagation Problem of Gas-Saturated Coal Rock

Secondly, we consider a column of gas-saturated coal rock with the same geometry, boundary condition, load, and finite element mesh as the first example. The Drucker–Prager criterion and the linear strain hardening rulec=c0+hpcε¯p for the cohesion are applied to simulate the elastic-plastic response of the gas-saturated coal rock. Here, c0, hpc, and ε¯p are the initial cohesion, strain hardening parameter, and the equivalent plastic strain, respectively. In the present example c0=5.0MPa, hpc=1MPa, the internal frictional angle ϕ=22°, and the plastic potential angle ψ=5° are used to account for the problem. Figures 4 and 5 give the resulting elastic-plastic effective stress wave and pore pressure wave propagation solutions at different time levels of 0.041s and 0.08s for the example using the Newmark method (blue lines) and the MDGFEM (red lines) with time step size Δt=0.001s. It is noted that the results obtained by different methods are different from each other. Serious spurious numerical oscillations could be observed both in effective stress and pressure profiles at different time levels by the Newmark method. By contrast, the smooth and continuous stress and pressure distributions are exhibited by using the MDGFEM. This indicates that artificial damping constant and discontinuous characteristic of basic unknown variables play significant roles in filtering out the spurious numerical oscillations and in providing much more accurate solutions for elastic-plastic wave propagation problem of gas-saturated coal rock subjected to a shock load.Figure 4
Elastic-plastic stress profiles of different methods att = 0.041 s and t = 0.08 s.Figure 5
Pressure profiles of different methods att = 0.041 s and t = 0.08 s.


### 4.3. 2-D Elastic-Plastic Wave Propagation Problem of Gas-Saturated Coal Rock

As the third example, we consider elastic-plastic wave propagation problem of gas-saturated coal rock in two dimensions. A square domain, which is 10 m deep, 10 m wide, and of infinite length in the horizontal direction, is subjected to an impulse inclined compression load at top left corner (0–1.4 m) as depicted in Figure6. The left face of the domain is symmetry boundary. The domain, which is meshed with elements, is analyzed as a plane strain problem.Figure 6
Modeling and mesh discretization for example 3.In the present example, the material property data are defined as follows:c0=5.0MPa, hpc=1MPa, ϕ=22°, and ψ=0. The impulse force can be expressed as(32)Ft=50KN,0≤t≤0.041s,0N,0.041s<t.Figures7 and 8 illustrate the numerical results for the mean stress distributions and the pore fluid pressure with the coal at time t = 0.001 s, 0.003 s, and 0.005 s, as obtained by the MDGFEM method and the Newmark method using the time step. It should be remarked that due to incapability of filtering out the effect of the high modes, the spurious numerical oscillation occurring at the free end after thewave tail passes in the numerical solution by the Newmark method. Figures 7 and 8 demonstrate the good performance of the MDGFEM in reducing the numerical oscillations in the wave-front stage, while the Newmark method fails to do so.Figure 7
Comparison of the stress distributions between the Newmark method and the DGFEM. (a) Newmark method 0.001 s, (b) MDGFEM 0.001 s, (c) Newmark method 0.003 s, (d) MDGFEM 0.003 s, (e) Newmark method 0.005 s, and (f) MDGFEM 0.005 s.
(a)(b)(c)(d)(e)(f)Figure 8
Comparison of the pressure distributions between the Newmark method and the MDGFEM. (a) Newmark method 0.001 s, (b) MDGFEM 0.001 s, (c) Newmark method 0.003 s, (d) MDGFEM 0.003 s, (e) Newmark method 0.005 s, and (f) MDGFEM 0.005 s.
(a)(b)(c)(d)(e)(f)


## 4.1. 1-D Elastic Wave Propagation Problem of Gas-Saturated Coal Rock

We first consider a one-dimensional stress wave propagation problem in saturated gas-saturated coal. The length of the column is 200 m, and the area of the cross section is 1 m2. The column, as shown in Figure 1, with one end fixed and one end loaded by an axial force impulse can be expressed by(31)Ft=2000N,0≤t≤0.041s,0N,0.041s<t.Figure 1
Model of example 1.The column is meshed in elements of 0.4 m along the column axis. Figures2 and 3 give the resulting stress wave and pore pressure wave propagation solutions at different time levels of 0.041s and 0.08s for the example using Newmark method, DGFEM, and the MDGFEM with time step size Δt=0.001s. It is noted that the results obtained by different methods are different from each other. As compared with the Newmark method, both DGFEMs filter out the serious numerical oscillation in the wave after. And the modified DGFEM shows better ability in filtering out the serious numerical oscillation in wave than DGFEM.Figure 2
Stress profiles of different methods att = 0.041 s and t = 0.08 s.Figure 3
Gas pressure profiles of different methods att = 0.041 s and t = 0.08 s.


## 4.2. 1-D Elastic-Plastic Stress Wave Propagation Problem of Gas-Saturated Coal Rock

Secondly, we consider a column of gas-saturated coal rock with the same geometry, boundary condition, load, and finite element mesh as the first example. The Drucker–Prager criterion and the linear strain hardening rulec=c0+hpcε¯p for the cohesion are applied to simulate the elastic-plastic response of the gas-saturated coal rock. Here, c0, hpc, and ε¯p are the initial cohesion, strain hardening parameter, and the equivalent plastic strain, respectively. In the present example c0=5.0MPa, hpc=1MPa, the internal frictional angle ϕ=22°, and the plastic potential angle ψ=5° are used to account for the problem. Figures 4 and 5 give the resulting elastic-plastic effective stress wave and pore pressure wave propagation solutions at different time levels of 0.041s and 0.08s for the example using the Newmark method (blue lines) and the MDGFEM (red lines) with time step size Δt=0.001s. It is noted that the results obtained by different methods are different from each other. Serious spurious numerical oscillations could be observed both in effective stress and pressure profiles at different time levels by the Newmark method. By contrast, the smooth and continuous stress and pressure distributions are exhibited by using the MDGFEM. This indicates that artificial damping constant and discontinuous characteristic of basic unknown variables play significant roles in filtering out the spurious numerical oscillations and in providing much more accurate solutions for elastic-plastic wave propagation problem of gas-saturated coal rock subjected to a shock load.Figure 4
Elastic-plastic stress profiles of different methods att = 0.041 s and t = 0.08 s.Figure 5
Pressure profiles of different methods att = 0.041 s and t = 0.08 s.


## 4.3. 2-D Elastic-Plastic Wave Propagation Problem of Gas-Saturated Coal Rock

As the third example, we consider elastic-plastic wave propagation problem of gas-saturated coal rock in two dimensions. A square domain, which is 10 m deep, 10 m wide, and of infinite length in the horizontal direction, is subjected to an impulse inclined compression load at top left corner (0–1.4 m) as depicted in Figure6. The left face of the domain is symmetry boundary. The domain, which is meshed with elements, is analyzed as a plane strain problem.Figure 6
Modeling and mesh discretization for example 3.In the present example, the material property data are defined as follows:c0=5.0MPa, hpc=1MPa, ϕ=22°, and ψ=0. The impulse force can be expressed as(32)Ft=50KN,0≤t≤0.041s,0N,0.041s<t.Figures7 and 8 illustrate the numerical results for the mean stress distributions and the pore fluid pressure with the coal at time t = 0.001 s, 0.003 s, and 0.005 s, as obtained by the MDGFEM method and the Newmark method using the time step. It should be remarked that due to incapability of filtering out the effect of the high modes, the spurious numerical oscillation occurring at the free end after thewave tail passes in the numerical solution by the Newmark method. Figures 7 and 8 demonstrate the good performance of the MDGFEM in reducing the numerical oscillations in the wave-front stage, while the Newmark method fails to do so.Figure 7
Comparison of the stress distributions between the Newmark method and the DGFEM. (a) Newmark method 0.001 s, (b) MDGFEM 0.001 s, (c) Newmark method 0.003 s, (d) MDGFEM 0.003 s, (e) Newmark method 0.005 s, and (f) MDGFEM 0.005 s.
(a)(b)(c)(d)(e)(f)Figure 8
Comparison of the pressure distributions between the Newmark method and the MDGFEM. (a) Newmark method 0.001 s, (b) MDGFEM 0.001 s, (c) Newmark method 0.003 s, (d) MDGFEM 0.003 s, (e) Newmark method 0.005 s, and (f) MDGFEM 0.005 s.
(a)(b)(c)(d)(e)(f)


## 5. Conclusions

The traditional Galerkin finite element method such as the Newmark method fails to capture the discontinuities or sharp gradients of the stress wave in solving the impact problem. An additional artificial damping discontinuous Galerkin numerical algorithm was incorporated into the final finite element discretised form to reduce the numerical oscillations in the wave-front stage for the impact problem of gas-saturated coal. Based on the MDGFEM, series of numerical examples of dynamic problem of gas-saturated coal under impact loading show that the modified discontinuous Galerkin finite element method can effectively filter out the spurious numerical oscillations in the wave front of the elastic-plastic stress wave and the gas pressure wave. The study of the method in this paper demonstrates that it may also serve as a viable method for the coal-rock-gas coupling problem and coal-gas-heat coupling problems subjected to impact loading.



---

*Source: 8845056-2020-08-29.xml*