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.
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]
For the pore fluid, the equation of motion can be written aswhere ; is the total Cauchy stress; is the body force acceleration; , , and are the densities of the assembly, the fluid and the solid phases, respectively; is the displacement of the fluid phase; is the displacement of the solid phase; and is the porosity.
Insertion of into equations (1) and (2) and subtraction of equation (2) from equation (1) leads to the equation of solid skeleton equilibrium.where is the generalized Biot effective stress tensor and . For an isotropic material, the Biot parameter , where .
Equation (2) can be written as
Then, the mass conservation equation of the fluid flow can be written as
Using the definition of equation (5), equations (1) and (2) can rewritten aswhere ; is the bulk modulus of the fluid; and is the bulk modulus of the solid. The linearized form of the semidiscretized systems of the motion equations (6) and (7) can be written asin whichwhere and are the nodal displacements of the solid skeleton and the nodal intrinsic displacements of the pore fluid, respectively, is a fourth-order tensor defining a constitutive law for the solid skeleton, and and are the spatial derivatives of the shape functions of the node with respect to the 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 and , respectively, one haswhere and are the prescribed values of fluid and solid displacement vector.
On the other hand, if surface loadings are applied to the corresponding surfaces and , the following boundary conditions should be satisfied:where and are the given surface loadings, respectively.
In the present paper, surface impulse loading is applied with the form asin which and 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:in which
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, . The temporal jump of the discontinuous function at can be defined asin which
The global displacement vectors of solid and fluid at arbitrary time in a time step are interpolated by using the third-order Hermite shape functions asin which and stand for the global nodal value of displacement vectors at times and and and 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 as follows:
The global velocity vector of solid and fluid at arbitrary time in a time step is interpolated as an independent variable by using the linear shape functions as follows:
The weak forms of the semidiscretized equation (15) and the constraint condition along with the discontinuity conditions and 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 and .
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:in which is the stiffness proportional matrix, is the damping constant, is the damping ratio estimated from the global behavior of the system, and 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 asin which
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.
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
The 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 and for the example using Newmark method, DGFEM, and the MDGFEM with time step size . 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.
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 rule for the cohesion are applied to simulate the elastic-plastic response of the gas-saturated coal rock. Here, , , and are the initial cohesion, strain hardening parameter, and the equivalent plastic strain, respectively. In the present example , , the internal frictional angle , and the plastic potential angle 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 and for the example using the Newmark method (blue lines) and the MDGFEM (red lines) with time step size . 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.
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 Figure 6. The left face of the domain is symmetry boundary. The 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: , , , and . The impulse force can be expressed as
Figures 7 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.
(a)
(b)
(c)
(d)
(e)
(f)
(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.
Data Availability
The data used to support the findings of this study are included within the article.
Conflicts of Interest
The authors declare that they have no conflicts of interest.