Modified stepped scheme for modelling the dynamic behaviour of 3 D poroviscoelastic solids

The problem of the dynamic response of a soil medium under different kinds of loads is of significant importance in various areas of engineering, especially in connection with structures. The present paper is dedicated to the modification of the numerical approach for modelling the dynamic behaviour of three dimensional poroviscoelastic solids. The basic equations for fluid-saturated porous media proposed by Biot are modified by replacing the classical linear elastic model of the solid skeleton with the viscoelastic model. Classical models of viscoelasticity are employed, such as Kelvin-Voight model, standard linear solid model and model with weakly singular kernel. Boundary integral equations method is applied to solving three-dimensional boundary-value problems. Stepped schemes modifications based on the linear and quadratic approximation of function are employed. A numerical example of poroviscoelastic rod under Heaviside type load is provided. A problem of a poroviscoelastic cube with a cavity subjected to a normal internal pressure is considered. The comparison of dynamic responses when poroviscoelastic material is described by different viscoelastic models is presented.


Introduction
Wave motion in porous media is an important issue in several branches of engineering science.Mechanics of such advanced materials is relevant to many disciplines, such as geophysics, geo-and biomechanics, seismology, constricting, petroleum engineering.Study of wave propagation processes in saturated porous media originates from Y. I. Frenkel [1] and M. Biot [2,3].The implementation of the solid viscoelastic effects in the theory of poroelasticity was first introduced by Biot [4].After that waves in saturated porous continua have been studied by many authors.Common state of the art can be found in works [5][6][7].Important, that there is still absent an analytical solution for 3d poroelastodynamics and large part of transient response problems can only be solved via numerical methods.
Classical formulations for boundary integral equation (BIE) method with their discretized realization and traditional boundary element method (BEM) are successful approaches for solving three-dimensional problem.Nowadays there are two fundamentally different approaches for constructing such numerical scheme: direct solving in time domain [8,9], or solving in Laplace or Fourier domain [10,11].First approach with traditional stepping schemes has some disadvantages, e.g.requires fundamental solutions in time, significant computing costs and low stability rate.New Convolution Quadrature Method (CQM) which is widely applied in construction of time-step boundary element schemes on the basis of fundamental solutions in Laplace domain was introduced in [12,13].CQM helps to get rid of fundamental solutions in time requirement and shows better stability but is still costly.So the Laplace transform is the main method for solving three-dimensional transient problems in porous media.Shanz and Cheng [14] presented an analytical solution in the Laplace domain and analytical time-domain solution without considering viscous coupling effect, and then governing equations was extended for saturated poroviscoelastic media by using the Kelvin-Voigt model and obtained an analytical solution in the Laplace domain for the 1D problem [15].In present work we modify basic equations for fluidsaturated porous media proposed by Biot by replacing the classical linear elastic model of the solid skeleton with the viscoelastic model.Thus, the new theory can take into account the viscoelastic effect of the solid skeleton.A stepped boundary element scheme, based on the time-step method for numerical inversion of Laplace transform is considered.A modification of the scheme with variable integration step based on highly oscillatory quadrature approach is performed.Also modifications based on the linear and quadratic approximation of function are employed.

Mathematical model
We consider a homogeneous solid  in three-dimensional space 3  R , and is the boundary of  .Solid  is taken as isotropic poroviscoelastic.Basic poroelastic material is a two-phase material consisting of an elastic skeleton and compressible fluid or gas filler.Porous material of a volume V can be constructed as follows: V is the summary pore volume and s V is the volume of the skeleton.It is assumed that filler can openly seep through the pores and all closed pores are assumed as a part of the skeleton.Then a correspondence principle is applied to the skeleton, so we extend poroelastic formulation to poroviscoelasticity.
Considering a boundary-value problem for Biot's model of fully saturated poroelastic continuum in Laplace domain in terms of 4 unknowns (displacements i u and pore pressure p ) the set of differential equations take the following form [5]: where u  and   denotes boundaries for boundary conditions of 1 st and 2 nd kind respectively,

Boundary element approach
Fundamental and singular solutions are considered in term of singularity isolation.Numerical scheme is based on the Green-Betti-Somigliana formula.Boundary-value problem ( 1)-( 2) can be reduced to the BIE system as follows [10,16]: where -fundamental and singular solutions, -is an arbitrary point.Coefficient   equals 1 in case of finite domain and -1 in case of infinite domain.Boundary surface of our homogeneous solid is discretized by quadrangular and triangular elements and triangular elements are assumed as singular quadrangular elements.We use reference elements: square , and each boundary element is mapped to a reference one by the following formula: , ) , ( where l -local node number in element k, ) , ( l k  -global node number, ) ( l N -shape functions.Goldshteyn's displacement-stress mixed model is performed.To discretize the boundary surface eight-node biquadratic quadrilateral elements are used, generalized displacements and tractions -are approximated by linear and constant shape functions, respectively.
Subsequent applying of collocation method leads to the system of linear equations.As the collocation nodes we take the approximation nodes of boundary functions.Gaussian quadrature is used to calculate integrals on regular elements.But if an element contains a singularity, algorithm of singularity avoiding or order reducing is applied.When singularity is excluded we use an adaptive integration algorithm.An appropriate order of Gaussian quadrature is chosen from primarily known necessary precision, if it is impossible, the element is subdivided to smaller elements recursively.Solving the system of linear equations leads to the solution of the initial boundary-value problem in Laplace domain.

The stepped scheme of numerical inversion of Laplace transform
Consider a method based on the theorem of the integration of the original -the stepped method of numerical inversion of Laplace transform.Consider the following integral: Integral (5) gives rise to Cauchy problem for an ordinary differential equation: Integral ( 5) is substituted for by a quadrature sum, weighting factors of which are determined using Laplace representation f and the linear multi-step method [5,[12][13].
Further derivation is based on the results of those works.The traditional stepped method of integrating the original consists in that integral (5) is calculated using the following relation: where where t ; R is parameter of the method.
The traditional method uses a constant-step trapezoid method for calculating such an integral.Consider the following formula of constructing n  based on a variable step: For the cases where  is a strongly oscillating function, a combined formula is introduced, which makes use of the specific character of the integration of such functions.
The modification based on the linear approximation of function   where 2 The potential of the constructed modifications of the convolution quadrature method is illustrated using numerical examples.Fig. 1 presents the originals of the sought functions for the following parameters of the scheme: 500 . Number 1 indicates the curve constructed using the traditional formula of the convolution quadrature method, numbers 2 and 3 indicate the curves constructed based on the results of using the modified method with linear and quadratic interpolations of function   . Use of the variable-step convolution quadrature method with linearly interpolating the integrand does not make it possible to solve the problem of stability of the numerical construction of the sought function over a chosen time interval without refining the computational grid.Fig. 2  , respectively [17].Usage of combined formulas ( 10), ( 11) makes it possible to obtain the sought result on the same grid (the smooth curve in Fig. 2).The curves for formulas (10) and (11)  The constructed modifications of the convolution quadrature method based on combined formulas (10), (11) made it possible to overcome such disadvantages of the traditional method as having to choose number of time N coinciding with number of nodes L along angle  .Use of the modifications makes it possible to reduce the number of discretization points required for providing an assigned accuracy.In the above example the number of such points was reduced by half, and by using smaller time steps the number of points is further reduced, as the informative part of function f with a decreasing time step will tend to points 0 and  2 (an example of a three-fold reduction of the number of points is given).
5 Numerical results

Method verification (rod under Heaviside type load)
The problem of poroviscoelastic rod under Heaviside type load Material constants are: The problem has analytical solution in Laplace domain [5] which just transforms back to time domain with the help of the stepped method for Laplace transform.On Fig. 3 the comparison of numerical and analytical solutions is shown.Maximum relative error is less than 2% at considered time interval.Our numerical approach demonstrates a satisfactory accuracy, so we can apply it to 3D problems.Captions should be typed in 9point Times.They should be centred above the tables and flush left beneath the figures.

Conclusions
Boundary integral equations method and boundary element method are applied in order to solve three dimensional boundary-value problems.Numerical approach based on BEM is verified by comparison with analytics.The poroviscoelastic m edia modelling is based on Biot's theory of porous material in combination with the elastic-viscoelastic corresponding principle.Viscous properties are described by Kelvin-Voight model, standard linear solid model and model with weakly singular kernel.Some poroviscoelastic solutions were earlier obtained with Durbin's method of numerical inversion usage in [18].In this work time-step method, based on CQM is developed.
A problem of Heaviside-type load acting on a poroviscoelastic clamped rod and a problem of poroviscoelastic cube with a cavity subjected to a normal internal pressure are solved.An influence of material viscosity on transient responses of displacements and pore pressure is presented.Obtained results for the transient response of the displacements show the effect of the transition from instantaneous to equilibrium viscoelastic moduli in case of the standard linear solid model.This work is supported by the Russian Science Foundation under grant 16-19-10237.
Biot to describe dynamic interaction between fluid and skeleton.C is a factor depending on the pores geometry and excitation frequency.Poroviscoelastic solution is obtained from poroelastic solution by means of the elasticviscoelastic correspondence principle, applied to skeleton's constants K and G in Laplace domain.Forms of functions ) chosen viscoelastic model.
The modification based on a quadratic approximation of function  

Fig. 3 .
Fig. 3. Numerical and analytical poroviscoelastic solutions comparison in case of weakly singular kernel model, h =1 and 7 .0  

5. 2
Cube with a cavity subjected to a normal internal pressure A poroviscoelastic cube, with side lengths m, 2  L and containing a spherical cavity of radius m 5 .0  R is considered (see Fig. 4).The cube is clamped at the face m 1 3   x and cavity is subjected to a normal internal pressure:

Fig. 5 -
Fig. 5-10 shows the transient response of the displacements ) ( 3 t u and pore pressure p in the center of upper face at the nodal point with the coordinates (0, 0, 1) m.For poroviscoelastic analysis, different viscoelastic models were considered.Viscosity parameter influence on dynamic responses is studied.
depicts numerical oscillation arising for the following parameters of the scheme: