Development of boundary-element time-step scheme in solving 3D poroelastodynamics problems

The development of time-step boundary-element scheme for the three dimensional boundaryvalue problems of poroelastodynamics is presented. The poroelastic continuum is described using Biot’s mathematical model. Poroelastic material is assumed to consist of a solid phase constituting an elastic formdefining skeleton and carrying most of the loading, and two fluid phases filling the pores. Dynamic equations of the poroelastic medium are written for unknown functions of displacement of the elastic skeleton and pore pressures of the filling materials. Green’s matrices and, based on it, boundary integral equations are written in Laplace domain. Discrete analogue are obtained by applying the collocation method to a regularized boundary integral equation. Boundary element scheme is based on time-step method of numerical inversion of Laplace transform. A modification of the time-step scheme on the nodes of RungeKutta methods is considered. The Runge-Kutta scheme is exemplified with 2and 3-stage Radau schemes. The results of comparing the two schemes in analyzing a numerical example are presented.


Introduction
In BEM-modeling of dynamic processes, two main approaches can be conventionally discerned: solving in time, using a time-step scheme [1], and solving in Laplace or Fourier transforms with the subsequent inversion of the transforms [2]. The potential of the traditional time-step schemes using spline approximations is limited due to the absence of fundamental solutions in time. In many cases, matrices of fundamental solutions can only be constructed in Fourier and Laplace images. That is why the pioneering formulations of the boundary-element method for the dynamics of poroelastic media of the Biot model were published in Laplace transforms [3,4]. In [5][6][7][8][9][10], to construct a time-step boundary-element scheme based on fundamental solutions in Laplace transforms, the quadrature convolution method introduced by Lyubich [11,12] is used. In [13], using the convolution quadrature method in combination with Runge-Kutta methods was considered. The combined use of Runge-Kutta methods and a BEM time-step scheme for increasing the accuracy of the result with reducing computational costs was further developed in [14][15][16]. Work [17] reviews a wide range of approaches to using a boundary-element scheme in combination with the convolution quadrature method, based on both Euler scheme and other schemes of the Runge-Kutta family. In [18][19][20][21], dynamic wave processes in elastic media with coupled fields are automatically modeled using an analogous boundary-element scheme based on the timestep method of numerically inverting Laplace transforms. The present study describes a modification of such a time-step method on the nodes of Runge-Kutta schemes and a time-step BEM scheme constructed on its basis.

Problem formulation
A set of fully coupled governing differential equations of a porous medium saturated by two compressible fluids (water and air) subjected to dynamic loadings is considered. In this formulation the solid skeleton displacements i u , water pressure w p and air pressure a p are presumed to be independent variables [22]. The final differential equations in Laplace domain yield where w K and a K are bulk moduli of the fluid, ϕ is porosity, are bulk body forces. The bulk density is denoted by where s ρ is the density of the solid, w ρ is wetting fluid density, a ρ is non-wetting fluid density. The saturation degrees are defined as the ratios of the volume occupied by the fluid w V or a V to the void volume, i.e. it holds The following abbreviations ) ( (17) are introduced, where rw S is the residual wetting fluid saturation and ra S is the non-wetting fluid entry saturation. The symbol d p is non-wetting fluid entry pressure, ϑ is the pore distribution index while the value of ϑ lies between 0.2 and 3, normally. The symbols β and γ are Laplace parameter dependent variables and expressed as where w κ and a κ the phase permeability of the wetting and the non-wetting fluid are given by respectively. rw K and ra K denotes the relative fluid phase permeability, k denotes the intrinsic fluid permeability, w η and a η are viscosity of the fluid. To evaluate relative phase permeability following equations are used

BEM application
The boundary-element technique is based on the use of a regularized BIE direct approach [23]: ( ) where ( , , ) s U x y and ( , , ) s T x y are matrices of fundamental and singular solutions, respectively [24], ( , ) 0 T x y contains isolated singularities, x is integration point, y is observation point, u is generalized displacement vector, t is generalized force vector.
To solve equation (21), the boundary surface is divided into generalized eight-node quadrangular elements; the coordinates of the points on the k element are determined from the relation x are global coordinates of nodes [20].
Generalized boundary functions of the first kind are approximated bilinearly, and generalized boundary functions of the second kind are assumed to be constant over the element: integration on an element that does not contain a collocation point, in addition to the Gauss integration formulas, a hierarchical integration algorithm is applied, wherein the element is subdivided until the specified accuracy is achieved.

Laplace transform inversion
Consider the Runge-Kutta method written down using Butcher's table: A correct formulation of a time-step scheme requires the following conditions to be met [15]: •The Runge-Kutta method must be A-stable; , the method is automatically L-stable.
Using the Runge-Kutta method instead of the linear multi-step method [25] makes it possible to obtain, taking into account [15], formulas for a modified time- In formulas (27)-(29), weighing coefficients ( ) n t ω ∆ and characteristic function ( ) z γ are m-order matrices. In the present study, 2-and 3-stage Radau schemes (Radau IIA) were chosen as a particular example of the Runge-Kutta schemes meeting the formulated conditions.

Numerical example
The problem of the effect of end force 2 ( ) 1 / F t N m = on a 3 m-long 1 m-wide prismatic poroelastic column with a rigidly fixed end is considered. The boundary-value problem is presented in Fig. 1. The parameters of the partially saturated porous material correspond to those of sandstone [25]. A boundary-element discretization consisting of 504 quadrangular elements (Fig. 1) is used in the computations. An analytical solution of the corresponding one-dimensional problem is used for the boundary element model verification [22].     . Depending on the number of stages in the scheme, in each case the computations required m L × nodes of the time-step scheme. It can be seen in Fig. 2 that use of the both schemes makes it possible to obtain good results, and that the boundary-element solutions are close to the analytical solution. In spite of the twice-longer time step and the lower number of the nodes, the solution at the nodes of the 3-stage scheme is smooth, and the oscillation amplitude at the wave front points does not exceed that of the solution at the nodes of the 2-stage scheme (Fig. 3). The effect of the time step length on the solution is demonstrated in Figs. 4-5. It is evident that choosing a shorter time step does not make it possible to improve the solution in the case of the 2-stage scheme, but results in increasing the oscillation frequency. This effect could not be suppressed by increasing the number of nodes (Fig. 4). A similar result was obtained when using the 3-stage scheme (Fig. 5). Thus, the initially chosen time step length for each of the schemes proves to be optimal for the given boundary-element grid. Besides, the optimal time step length is bigger in the case of the 3-stage scheme, which corroborates the results of the study [15].

Conclusion
A time-step boundary-element scheme on the nodes of a family of Runge-Kutta methods for analyzing initial boundary-value problems of the dynamics of partially saturated poroelastic bodies has been considered. The results of a numerical experiment obtained using 2-and 3-stage Radau schemes are presented. The boundaryelement solutions are in good agreement with the analytical solution. The results of the comparative analysis of using the two schemes agree with the results published by other authors.