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

School of Emergency Management and Safety Engineering, China University of Mine and Technology, Beijing 100083, China School of Mechanics and Safety Engineering, Zhengzhou University, Zhouzhou 450001, China 'e State Key Laboratory for Structural Analysis of Industrial Equipment, Department of Engineering Mechanics, Faculty of Vehicle Engineering and Mechanics, Dalian University of Technology, Dalian 116024, China


Introduction
With the increase of coal mining depth, coal and gas outburst disasters are becoming more and more severe [1]. e 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][2][3][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. e 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.
is 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][8][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][11][12]. e 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 gassaturated 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 wavefront 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). e 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.

The Basic Dynamics Governing Equations and DGFEM Formulation for Coal and Rock Porous Media
is 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] For the pore fluid, the equation of motion can be written as where (R i /n) � k − 1 w,ij _ w j ; σ ij is the total Cauchy stress; b i is the body force acceleration; ρ, ρ f , and ρ s are the densities of the assembly, the fluid and the solid phases, respectively; w k is the displacement of the fluid phase; u i is the displacement of the solid phase; and n is the porosity.
Insertion of U i � u i + (1/n)w i into equations (1) and (2) and subtraction of equation (2) × n from equation (1) leads to the equation of solid skeleton equilibrium.
Equation (2) can be written as en, the mass conservation equation of the fluid flow can be written as Using the definition of equation (5), equations (1) and (2) can rewritten as where (1/Q) � (n/K f ) + (α − n/K s ); K f is the bulk modulus of the fluid; and K s is the bulk modulus of the solid. e linearized form of the semidiscretized systems of the motion equations (6) and (7) can be written as where u(t) and U(t) are the nodal displacements of the solid skeleton and the nodal intrinsic displacements of the pore fluid, respectively, D ijkl is a fourth-order tensor defining a constitutive law for the solid skeleton, and N u K,i and N U K,i are the spatial derivatives of the shape functions of the K th 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 where U 0 (t) and u 0 (t) 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: where f u0 (t) and f U0 (t) are the given surface loadings, respectively.
In the present paper, surface impulse loading is applied with the form as in which P and t 0 are the amplitude of loading and the constant time.

Time-Discontinuous Galerkin Finite Element
Method. e standard Galerkin discretized equation (8) can be rewritten as follows: in which e 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 · · · < t n < · · · t N . e temporal jump of the discontinuous function d(t n ) at t n can be defined as e global displacement vectors of solid and fluid at arbitrary time t in a time step I n � (t + n , t − n+1 ) are interpolated by using the third-order Hermite shape functions as Shock and Vibration in which d + n and d − n+1 stand for the global nodal value of displacement vectors at times t n and t n+1 and v + n and v − n+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: e global velocity vector of solid and fluid at arbitrary time t in a time step I n � (t + n , t − n+1 ) is interpolated as an independent variable by using the linear shape functions as follows: e 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 Substituting equations (21) and (22) into (15), we obtain the following matrix equation of DGFEM for solving from independent variations: in which f e 1 � (Δt/3)f e n + (Δt/6)f e n+1 and f e 2 � (Δt/6)f e n + (Δt/3)f e n+1 . 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][11][12][13][14][15]. To demonstrate the advantage of the present DGFEM, we have to solve several problems in our previous work [10][11][12][13][14][15]. For the elastic-plastic dynamic problems, the Newton-Raphson process is used [10,14].

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][14][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. e equations of the damping scheme can be written as follows: 4 Shock and Vibration 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][14][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 in which Equations (24)

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.

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. e length of the column is 200 m, and the area of the cross section is 1 m 2 . e column, as shown in Figure 1, with one end fixed and one end loaded by an axial force impulse can be expressed by e column is meshed in elements of 0.4 m along the column axis. Figures 2 and 3 give the resulting stress wave and pore pressure wave propagation solutions at different time levels of 0.041 s and 0.08 s for the example using Newmark method, DGFEM, and the MDGFEM with time step size Δt � 0.001 s. 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.

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.
e Drucker-Prager criterion and the linear strain hardening rule c � c 0 + h c p ε p for the cohesion are applied to simulate the elastic-plastic response of the gas-saturated coal rock. Here, c 0 , h c p , and ε p are the initial cohesion, strain hardening parameter, and the equivalent plastic strain, respectively. In the present example c 0 � 5.0MPa, h c p � 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.041 s and 0.08 s for the example using the Newmark method (blue lines) and the MDGFEM (red lines) with time step size Δt � 0.001 s. 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. is 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.

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 Figure 6. e left face of the domain is symmetry boundary. e domain, which is meshed with elements, is analyzed as a plane strain problem.
In the present example, the material property data are defined as follows: c 0 � 5.0MPa, h c p � 1MPa, ϕ � 22°, and ψ � 0. e impulse force can be expressed as Shock and Vibration  Shock and Vibration 7

Conclusions
e 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. e 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-gasheat coupling problems subjected to impact loading.

Data Availability
e data used to support the findings of this study are included within the article.

Conflicts of Interest
e authors declare that they have no conflicts of interest.