A combined application of boundary-element and Runge-Kutta methods in three-dimensional elasticity and poroelasticity

The report presents the development of the time-boundary element methodology and a description of the related software based on a stepped method of numerical inversion of the integral Laplace transform in combination with a family of Runge-Kutta methods for analyzing 3-D mixed initial boundary-value problems of the dynamics of inhomogeneous elastic and poro-elastic bodies. The results of the numerical investigation are presented. The investigation methodology is based on directapproach boundary integral equations of 3-D isotropic linear theories of elasticity and poroelasticity in Laplace transforms. Poroelastic media are described using Biot models with four and five base functions. With the help of the boundary-element method, solutions in time are obtained, using the stepped method of numerically inverting Laplace transform on the nodes of Runge-Kutta methods. The boundary-element method is used in combination with the collocation method, local element-byelement approximation based on the matched interpolation model. The results of analyzing wave problems of the effect of a non-stationary force on elastic and poroelastic finite bodies, a poroelastic half-space (also with a fictitious boundary) and a layered half-space weakened by a cavity, and a half-space with a trench are presented. Excitation of a slow wave in a poroelastic medium is studied, using the stepped BEM-scheme on the nodes of Runge-Kutta methods. 1. Problem formulation The governing equations for the elasticity problem in Laplace domain are as follows L ik ūk(x, s) + F̄i = ρ s 2 ūi (x, s), (1) ū(x, s) = ũ, x ∈ , t̄n(x, s) = t̃n, x ∈ σ , (2) where L0 ik = G + (K + G/3)grad div, u denotes Dirichlet boundary and σ – Neumann boundary, G, K are elastic moduli, F̄i is bulk body force, ρ is material density, s is the Laplace transform parameter. The boundary-value problem for full Biot’s model of linear saturated poroelastic continuum in Laplace domain concerning 4 basic functions – skeleton displacements ūi and pore pressure p̄ – takes the following form [1]: Gūi, j j + ( K + 3 ) ū j,i j − (α − β) p̄,i = s2(ρ − βρ f )ūi − F̄ (3) β sρ f p̄,i i − φ 2s R p̄ − (α − β)sūi,i = −ā, x ∈ , (4) ū′(x, s) = ũ′, x ∈ , ū′ = (ū1, ū2, ū3, p̄) , (5) t̄ ′ n(x, s) = t̃ ′ n, x ∈ σ , t̄ ′ = (t̄1, t̄2, t̄3, q̄) (6) where φ is porosity, F̄i , ā are bulk body forces, β = κρ f φ 2s φ2 + sκ(ρe + φρ f ) , α = 1 − K Ks , (7) a Corresponding author: Igumnov@mech.unn.ru R = φ2 K f K 2 s K f (Ks − K ) + φKs(Ks − K f ) (8) – constants describe the interaction between the skeleton and filler, κ is permeability, ρ, ρe, ρ f are material density, apparent mass density and filler density respectively, Ks, K f are elastic bulk moduli of the skeleton and filler respectively. The governing equations of partially saturated poroelasticity in the Laplace domain with five unknowns – solid displacements ui , the pore wetting fluid pressure p, and the pore non-wetting fluid pressure pa – given by [2] Gūi, j j + (K + 3 )ū j,i j − (ρ − βSwρw − γ Saρa)sūi = = (α − β)Sw p̄ ,i − (α − γ )Sa p̄a ,i − F̄i , (9) −(α − β)Swsūi,i − (ζ − Saa Sw + Su)s p̄a + βSw ρws p̄ ,i i

The boundary-value problem for full Biot's model of linear saturated poroelastic continuum in Laplace domain concerning 4 basic functions -skeleton displacements ūi and pore pressure p -takes the following form [1]: ū (x, s) = ũ , x ∈ u , ū = ( ū1 , ū2 , ū3 , p) , (5) t n (x, s) = t n , x ∈ σ , t = ( t1 , t2 , t3 , q) (6) where φ is porosity, Fi , ā are bulk body forces, a Corresponding author: Igumnov@mech.unn.ru -constants describe the interaction between the skeleton and filler, κ is permeability, ρ, ρ e , ρ f are material density, apparent mass density and filler density respectively, K s , K f are elastic bulk moduli of the skeleton and filler respectively.The governing equations of partially saturated poroelasticity in the Laplace domain with five unknowns -solid displacements u i , the pore wetting fluid pressure p w , and the pore non-wetting fluid pressure p a -given by [2] G ūi, EPJ Web of Conferences Where K w and K a are bulk moduli of the fluid, φ is porosity, Fi , Ī w , Ī a 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 V w or V a to the void volume, i.e. it holds The following abbreviations S aa = S a + ϑ(S w − S r w ) ( 16) are introduced, where S r w is the residual wetting fluid saturation and S ra is the non-wetting fluid entry saturation.The symbol p d 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 κ w = K r w k/η w and κ a = K ra k/η a respectively.K r w and K ra 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 ) S e denotes the effective wetting fluid saturation degree given by For the extreme case S w = 0, the equations turn out to be simple a elasticity problem, for the extreme case S w = 1saturated poroelacticity problem.
To approximate the boundary consider its decomposition to a set of quadrangular and triangular 8-node biquadratic elements, where triangular elements are treated as singular quadrangular.Every element is mapped to a reference one (square . Interpolation nodes for boundary unknowns are a subset of geometrical boundaryelement grid nodes.Local approximation follows the Goldshteyn's generalized displacement-stress matched model [4]: generalized boundary displacements are approximated by bilinear elements whereas generalized tractions are approximated by constant.For BIE discretization the collocation method is used; the set of collocation nodes coincides with the set of approximation nodes for the boundary functions.Integrals in discretized BIE's are calculated using Gaussian quadrature in combination with singularity decreasing and eliminating algorithms [4].

Laplace transform inversion
CQM and time-step method for numerical Laplace transform inversion are quite close to each other but whereas CQM is based on the convolution theorem and intended for convolution integral calculation, time-step method for numerical Laplace transform inversion is based on the integration theorem and dedicated to the calculation of the original function integral.
In order to employ time-step method for obtaining the original function f (t) such as f (0) = 0 from known image f (s) we need to apply the integration theorem to its derivative: Following the theorem, the integral can be found as (basing on [2][3][4]): Here R is the radius of the analyticity region for f (γ (z)/ t) and γ (z) is the characteristic function for the linear multistep method applied to the Cauchy problem arising within the integral evaluation.Backward 04026-p.2differentiation (BDF) based methods of order ≤ 6 are applicable in scope of this solution scheme.For BDF-2 we have γ (z) = 3/2 − 2z + z 2 /2.For Runge-Kutta method applied instead of linear multistep method we obtain (based on [5][6][7]): (26) Runge-Kutta method is referred to as Butcher tableau Formula for coefficients ω n ( t) retains its written form with the only difference that now it's a matrix formula.The characteristic function writes as

Poroelastic column
Consider a problem of a poroelastic column of length l with one clamped end and a Heaviside-shaped in time load f (t) = S 0 H (t) to another end [8].Assuming the values S 0 = 1 N/m 2 and l = 10 m, the pore pressures p w , p a at the clamped end and displacement u 3 at the another  The problem is solved with the help of BDF-2 scheme with Ł=4N=1000; a mesh of 896 boundary elements is chosen.Figures 1-3 show the displacements and pore pressures.
For the nearly saturated case S w = 0.99999, the results are in good agreement with those of the saturated case.
One of the research aspects was modeling slow compression wave in terms of the half-space.For this purpose point E at a depth 3 m under the excitation area is considered.To calculate generalized displacement and  The study validates and improves the results of [9].As compared to [9], here all basic functions including displacements, tractions, pore pressure and flux are drawn; a mesh of 1176 boundary elements is used whereas [9] uses 342 triangular elements.

Layered half-space
A problem of a 2-layer poroelastic half-space subjected to a Heaviside-type load is considered (Fig. 4).All the parameters are the same as for the problem of homogeneous half-space except that here the layers are of different material properties: a layer of sandstone is located over the rock half-space.The sandstone properties are: 7, 8 represent vertical displacements at point A 10 m from the excitation area for sandstone layer of thickness h = 5 m and h = 10 m respectively.
The Rayleigh wave appearance (t ≈ 0.043 sec) can be observed after the fast compression wave arrival (t ≈ 0.01 sec).Then (t > 0.05 sec), some perturbation can as well be seen, that argues for the arrival of compression waves, reflected by the sandstone/rock boundary.The chosen sandstone layer thickness still allows distinguishing the effect of Rayleigh wave thanks to the fact that reflected compression waves arrive later.The results here as compared to ones in [9] rely on the matched displacement-stress model whose usage is better motivated than of isoparametric one.

Half-space with a cavity
A problem of vertical load to a surface element of poroelastic rock half-space with spherical or cubic cavity is considered.The load P(t) = P 0 f (t), P 0 = 1000 N/m 2 acts on the area 1 m 2 ; the cube side/sphere diameter is 10 m; the center of cavity is located at a depth of h = 7.5 m under the excitation area.The half-space surface is impermeable: q = 0.The problem is treated with the help of Radau-based scheme with N = 250 and Ł= 1000 on [0, π].The boundary element mesh takes into account 2 symmetry planes; a quarter of the mesh contains 04026-p.4420 elements on the surface and 324 elements on the cavity boundary.Responses of displacements and pore pressure at point 15 m from the excitation area are shown in Fig. 9 in comparison with the case of poroelastic half-space without cavity.Figure 10 presents the results obtained using the drained and undrained models.It's seen from the diagrams that the shape of the cavity affects the responses form qualitatively. Displacements for the case without cavity are closer to the spherical cavity-free case than to the cubic cavity case.However there's still certain difference.For the cavity case vertical displacements gain one more peak (t ≈ 0.011 sec) following the arrival of Rayleigh wave in the central part of the diagram (t ≈ 0.0095 sec); in the horizontal displacements diagram the peak right before the arrival of Rayleigh wave becomes much stronger (t ≈ 0.009 sec).

Half-space with a trench
A problem of vertical load applied to a surface element of poroelastic half-space with open trench is considered [8].The load P(t) = P 0 f (t), where P 0 = 1000 N/m 2 and f (t) is Heaviside function, acts on the area 0.25 m 2 and remaining surface is traction free.The pore pressure are assumed to be zero all over the surface.The trench with parameters L = 1 m, W = 0.5 m, H = 1 m is located at a distance D = 1 m 2 from the loading area (Fig. 11   The Radau-based time-step scheme of the BEM gives results that are close to those based on the traditional stepped scheme but is much more economical: it requires only 30 time-steps as compared with 100 for the Eulerbased scheme, and 121 double nodes for argument ϕ as compared with 415 single ones for the Euler-based scheme.

Concluding remarks
A boundary element approach utilizing Runge-Kutta based time-step method for numerical Laplace transform inversion is developed and applied to 3d poroelastodynamic boundary value problems.In combination with varied integration step and highly oscillatory quadrature employed here it helps to reduce computational costs for solution.It's seen from figures that the scheme gives results close to ones from the analogous scheme relying on traditional Euler, and better.As compared to BDF2, Runge-Kutta gives no delay of the final solution and better keeps its form in time.