Application of perfectly matched layer to soil-foundation interaction analysis

Mohammad Davoodi,Abbas Pourdeilami,Hosein Jahankhah,Mohammad Kazem Jafari

International Institute of Earthquake Engineering and Seismology,Tehran,Iran

Keywords:Perfectly matched layer(PML)Finite element method(FEM)Wave propagation Impedance/compliance

ABSTRACT Despite of the limitation in modeling in finite space,the finite element method(FEM)is one of the most used tools to numerically study the geotechnical problems regarding the capacity of simulating different geometries,conditions and material behaviors.A kind of absorbing layer named perfectly matched layer(PML)has been applied to modeling the radiation damping using FEM,which makes the dynamic analysis of soil-structure interaction more accurate.The PML is capable of absorbing incident waves under any angle and frequency,ensuring them to pass through the model boundaries without reflection.In this context,a new FEM program has been written and the PML formula has been implemented by rewriting the dynamic equation of motion and deriving new properties for the quadrilateral elements.The analysis of soil-foundation interaction by applying the PML is validated by the evaluation of impedance/compliance functions for different ground conditions.The results obtained from the PML model match the extended mesh results,even though the domain is small enough that other types of absorbing boundaries can reflect waves back to the foundation.The mechanism of the wave propagation in the region shows that the forced vibrations can be fully absorbed and damped by the boundaries surrounded by PMLs which is the role of radiation damping in FEM modeling.

1.Introduction

Simulation of wave propagation in geomechanics is an interesting issue in dynamic analysis.Particularly,the analysis of soilfoundation interaction is one of the most challenging issues in geotechnical engineering.The near- field part of the foundation is often modeled by finite element method(FEM)which is a powerful tool for analyzing complex geometries and nonlinear behaviors.In order to model the far- field part of the foundation,previous studies such as Lysmer boundary conditions(Lysmer and Kuhlemeyer,1969),in finite element(Kim and Yun,2000),rational boundary conditions(Feltrin,1997),Dirichlet to Neuman mappings(Givoli,1999),boundary element method(Yazdchi et al.,1999),scaled boundary element method(Song and Wolf,2000;Wolf and Song,2000),discontinuous Galerkin methods(Park and Tassoulas,2002;Park and Antin,2004)and high-order non-reflecting boundaryconditions(Givoli,2004)have beenpresented taking into account the wave propagation towards in finity in the analysis.However,it is needed to find a method that can be used to simulate the radiation damping efficiently.In the present study,a kind of absorbing layer named perfectly matched layer(PML)has been incorporated into the FEM to simulate the radiation condition in transient dynamic analysis of soil-foundation rock systems.

The PML is able to absorb incident waves at any angle and frequency,ensuring them to pass through the model boundaries without reflection.The PML procedure designed by Berenger(1994)was applied to an electromagnetic problem.The method can be interpreted as a complex transformation of space coordinates and a special mesh system is introduced to replace the in finitepart.The propagating and evanescent waves attenuate in an exponential trend after entering the PML.Parameters of the PML are chosen so as to ensure suf ficiently fast attenuation within the layer.Thus,the reflection caused by simple boundary conditions at the outer PML boundary will be negligibly small.

One of the requirements for near-perfect absorption characteristics is to utilize a PML grid whose properties and stretching functions vary slowly from element to element,which makes it possible to assume constant properties within any element.In the absence of discretization,even a step function selected as stretching function in the PML would be reflectionless.Discretization breaks this property,thus one falls back on the ‘a(chǎn)diabatic theorem’,i.e.any material change in a wave equation,if it is made gradual enough,will be reflectionless(Oskooi et al.,2008).In our case,the same situation holds true:all that matters are using the parameters to change gradually(Kausel and Barbosa,2012).

Considering advantages of the PML,researchers have applied the idea to different physical problems such as the solution of Maxwell equation(Jiao and Jin,2001),acoustic wave propagation(Nataf,2006),and poroelastic media(Martin et al.,2008).

Chew and Liu(1996)employed complex coordinates to define PMLs and showed that the resultant medium could absorb propagating waves.Chew et al.(1997)changed the variables to transform Maxwell’s equations in PML media into a complex coordinate system.By introducing complex coordinates,Liu(1999)developed PMLs in cylindrical and spherical coordinates in time domain.Teixeira and Chew(2000)discussed the interpretation of the PML absorbing boundary condition as an analytic continuation of the coordinate space to a spatial domain with complex variables(complex space).They reviewed the generalization of the PML to curvilinear coordinates and to general linear media using this rationale.Harari et al.(2000)presented a new finite element formulation to use PML for the time-harmonic analysis of acoustic waves.Collino and Tsogka(2001)established a PML model using the split- field approach for a general hyperbolic system and implemented it in the linear elastodynamic problem for an anisotropic medium.Using a finite difference method and the approach of complex coordinates,Zeng et al.(2001)extended the PML approach to truncate unbounded poroelastic media.Zheng and Huang(2002)developed anisotropic PMLs for elastic waves in Cartesian,cylindrical and spherical coordinates.Bécache et al.(2003)theoretically investigated well-posedness and stability using PMLs for anisotropic elastic waves.Festa and Nielsen(2003)evaluated analytically and numerically the reflection of body and Rayleigh waves caused by the discrete properties of the PML,under variable angles of incidence wavelengths.They showed that the thickness of PML can be kept minimal for studies involving relativelylow frequencies,and no rescaling with model size is required.In addition,they showed through numerical examples that a major advantage of PMLs is their efficiency in absorbing Rayleigh waves at the free surface,a point where other more classical methods perform rather poorly.

Basu and Chopra(2003)defined the PMLs by employing complex coordinates to solve time-harmonic elastodynamic equations by finite element implementation.Furthermore,they transformed the frequency domain equations into time domain equations and presented an approach to solve the resultant equations(Basu and Chopra,2004).They indicated that the finite element implementation of the anti-plane PML is symmetric,whereas that of the plane-strain PML is not.Moreover,Basu(2009)presented a threedimensional explicit finite element implementation of PML.

Appel? and Kreiss(2006)utilized a modal PML into the equations of linear elasticity,resulting in better stability properties than previous split- field models.Harari and Albocher(2006)suggested simplified finite element implementation for a PML in an elastic medium which utilizes standard shape functions.They also presented some guidelines for proper selection of the parameters of PML.By explicit FEM and one-point integration scheme,Ma and Liu(2006)presented an easy implementation of PML.Qin et al.(2009)introduced auxiliary variables to divide the PML wave equation in the frequency domain into two parts of normal and attenuated terms.By this,theyavoided convolution integrals in equations after the transformation to time domain.They adopted the finite difference method to propose a novel numerical implementation approach for PML absorbing boundary conditions.Liu et al.(2009)utilized the Crank-Nicolson scheme together with several algorithms to calculate the SH wave equations.Kang and Kallivokas(2010)used the PML with a hybrid finite element formula(displacement-stress).Lancioni(2012)compared the performance of the PML and high-order non-reflecting boundary conditions in a one-dimensional dispersive wave equation.They considered linear,quadratic and cubic polynomial stretching functions for timedependent wave problems.Comparing the performance of the PMLs with the other types of absorbing boundary conditions(ABCs),they pointed out the merits and drawbacks of the two methods.Kim and Pasciak(2010)developed a Cartesian PML for solving Helmholtz equation on a two-dimensional(2D)unbounded medium.

Khazaee and Lot fi(2014)employed the PMLs in the dynamic analysisofdam-reservoirsystems.They introduced proper boundary conditions to the formulation of the PML area in the reservoir.Their results show that the PML approach is a very ef ficient method for the time-harmonic and transient analysis of damreservoir systems if boundary conditions of the PML domain are included.Farzanian et al.(2016)developed a displacement-based finite element model for unbounded heterogeneous domains with radiation damping produced by the PML.They validated the heterogeneous model using the closed-form solution of a free rod with two-part modulus subjected to a specified time history.

Literature review shows that the application of the PML to solving the dynamic and seismic problems,particularly regarding wave propagation phenomena,is a challenging issue to be investigated.In fact,there are extensive situations in the field of elastodynamics and soil-structure interaction that would be possible to represent precise responses and they can be marked as exact solutions by utilizing PML.

As stated earlier,wave propagation modeling for in finite medium is not an easy task,because the wave reflection due to incident of the progressive motion to the artificial boundaries will disturb the responses.Moreover,in the FEM,wave propagation is dependent on domain discretization and time-integration techniques.The key to have a proper modeling in the FEM is implementing radiation damping condition.This kind of damping decreases the amplitude because of energy distribution on a larger volume of medium.

2.Problem definition

Fig.1.Input motion(normal probability distribution)at point A to determine the wave propagation velocity.

Fig.2.FEM model for measuring the velocities of compressive and Rayleigh waves(unit:m).

By using the FEM,a new program has been generated based on sparse matrix coding in the present research.This study focuses on the calibration of the transient dynamic analysis of soil-foundation by FE-PML formula.In fact,The PML is combined with a finite element scheme,which is used for the simulation of 2D(planestrain)wave propagation.The near- field part of the foundation is modeled by FEM surrounded by a PML which can model the far field by absorbing waves propagating towards in finity.

In the present study,we focused on the applicability of PML in radiation damping simulation.To do this,the vibration of rigid massless shallow foundation has been analyzed and the results are presented in the form of compliance/impedance functions.It is worth mentioning that,when using the PML,the design of machine foundation is one of the applications of the compliance/impedance functions.In fact,with the vibration of the foundation,different surface and body waves would propagate in the medium;as approaching to the inappropriate model boundaries,the waves may be reflected into the medium,leading to disturbed responses.Therefore,utilizing PML in the FEM could result in proper wave propagation in the domain of the solutions and the radiation damping can be simulated properly.

PML is capable of reducing the size of the modeling and time duration of the analysis with high accuracy so that it can be applied in problems regarding surface dynamic loading.The efficiency of the PML in the surface dynamic loading has been investigated by solving classic soil-structure interaction problems and comparing the results with the pertinent literature,which includes generally semi-analytical approaches.In addition,to validate the 2D planestrain approach,the results are compared with the experimental data of a physical modeling.

3.Governing equations

The formulation introduced by Basu and Chopra(2004)is implemented using a standard displacement-based FEM.The governing equations of a PML are most naturally defined in thefrequency domain using frequency-dependent,complex-valued coordinatestretching.Byapplying an inverse Fourier transform,the relations will be generated in time domain.To create a perfectly matched medium,the ximust be replaced with a coordinatestretch():

Table 1 Material properties of the region for measuring the velocities of compressive and Rayleigh waves.

Applying the relations in terms of,the equations of PMLs in the foundation will be defined as(Basu and Chopra,2004):

Because multiplication or division by the factor iωin the frequency domain corresponds to a derivative or an integral,respectively,in the time domain,time-harmonic equations are easily transformed into corresponding equations for transient motion if the frequency-dependence of the former is only a simple dependence on this factor.Therefore,the stretching functions can be written in the form of(Basu and Chopra,2004):

Fig.3.Vertical displacements for two points(a)with 10 m spacing in the depth of the model,and(b)with 10.18 m spacing on the surface of the model.

Fig.4.FEM model for shear wave velocity measurement(unit:m).

Fig.5.Horizontal displacements for three points of the model with a distance of 15 m in depth.

where csindicates the shear wave velocity;b indicates the foundation half-width;,,Fe,and Fpindicate the attenuation matrices.Eq.(3)is pre-multiplied by iωΛ-Tand postmultiplied by Λ-1,Eqs.(7)and(8)are substituted into Eq.(3),and then byinverse Fourier transform,the time domain equations for PML are obtained:

where fm,fcand fkare defined as follows:

and Σ and E are as follows:

Fig.6.Input force vibration of Ricker type.

Fig.7.Geometry of half-plane model(unit:m).

Eqs.(9a)-(9c)is implemented utilizing the standard displacement-based FEM.PML element matrices are derived using the principle of virtual work to obtain the weak form of the governing equations over the entire domainΩas follows:

whereΓ = ?Ω is the boundary ofΩ,and n is the unit normal.

Interpolating u as the displacement field and w as the weighting function in terms of shape functions N,and integrating over the element lead to the stiffness and mass matrices expressed in terms of nodal submatrices as follows:

The element-level internal force is defined as

In the above equation,σn+1and Σn+1are the vectorial parameters and matricesandare formed as

Table 2 Material properties of the half-plane.

Fig.8.(a)Vertical displacement of foundation over half-plane withν=0.4 andξ=0.05;(b)Horizontal displacement of foundation over half-plane withν=0.4 andξ=0.05;(c)Vertical displacement of foundation over half-plane withν=0.33 andξ=0,and(d)Horizontal displacement of foundation over half-plane withν=0.33 andξ=0.

4.Program verification

In the present research,a 2D finite element code has been developed based on the above PML formulation by using MATLAB mathematics language that is able to reduce the data storage and programming quickness through vectorization and sparsification.Utilizing four-noded quadrilateral elements enhanced with Newmark implicit time integration scheme,the code is specialized for the dynamic linear analysis in time domain.In order to obtain an unconditionally stable solution,Newmarkαandβparameters have to satisfy the following conditions:β≥0.5 andα≥0.25(0.5+β)2.For an average acceleration scheme,αandβcan be determined as 0.25 and 0.5,respectively(Bathe,1996).

However,selecting the logical time step for any loading analysis and noting the minimum element size to pass the loading frequencies through the model are mandatory so as to obtain proper results.Generally,the code is able to model different types of geometries,multiple layers,dynamic and seismic analysis and significantly,utilizing PML as radiation damping simulator;this is an advantage of the current code in comparison with the similar software.

For accurate representation of wave transmission through a model,the mesh size must be smaller than approximately onetenth to one-eighth of the wavelength associated with the highest frequency component of the input motion(Lysmer and Kuhlemeyer,1969).In the present code,Rayleigh damping which is proportional to the mass and stiffness of the system is used to provide material damping that is almost frequency-independent over a restricted range of frequencies.The PML formulation has been implemented in the code by rewriting the dynamic equation of motion and deriving new properties for the quadrilateral elements.

In order to verify the dynamic performance of the program,an example of a plane-strain model has been established under the vertical incident wave of a normal probability distribution type on the surface centerline.In this model,the mechanism of wave propagation,the velocities of compressive wave(P-wave)and Rayleigh wave have been measured.In addition,a base excitation has been applied to measuring the velocity of shear wave(SV-wave).The time history and frequency content of the input motion are shown in Fig.1.

The model discretizes a linear elastic region of 28 m×20 m using a mesh system with the maximum mesh size of 1.25 m(see Fig.2).The properties of the model are presented in Table 1.

The centerline of the model surface(point A)is subjected to vertical incident Ricker wavelet.Using the elasticity theory,the velocity of P-wave is calculated as

Fig.9.Contours of(a)horizontal and(b)vertical displacements for the half-plane(unit:m).

Fig.10.Geometry and discretization of half-space model without PML(unit:m).

Fig.11.(a)Horizontal and(b)vertical displacements of foundation on half-space with and without PML(ν=0.33 and ξ=0).

Fig.12.Impedance functions for foundation over half-plane withν=0.33 andξ=0.

Also,the velocity of Rayleigh wave can be calculated theoretically by an equation given by Knopoff(1952):

Specifying two nodes in the depth and two nodes on the surface of the model,time histories of the variables such as displacement and velocity are recorded.Having the distance between the nodes and marking the correspondence time of variable peaks,the wave velocity could be calculated.In Fig.3,the procedure is demonstrated through bodyand surface waves.Therefore,the velocities of compressive and Rayleigh waves have been obtained,equal to 45.5 m/s and 22.1 m/s,respectively.These results are nearly the same as that of elasticity theory.

The shear wave velocity has been calculated by applying a base excitation of normal probability distribution type,assuming that half-space is truncated by a 5 m length PML at the depth of 30 m(Fig.4).Using the elasticity theory,the shear wave velocity is calculated as where G indicates the shear modulus.

The outputof the base excitation is illustratedin Fig.5.The shear wave velocity has been obtained equal to 25.4 m/s,which matches the elasticity theory.

5.Comparison of FE-PML solution with others

For comparison,the time domain response of the massless rigid foundation is presented for the two models of extended mesh and PML.In addition,impedance/compliance functions are calculated and compared with analytical or semi-analytical works of the previous researchers.Next,a brief introduction to the impedance/compliance functions is presented and then the main studies from the results are declared.

5.1.Dynamic impedance functions

An important step in current methods used for dynamic analysis of rigid massive machine foundations is the determination(using analytical or numerical methods)of the dynamic impedance functions,K(ω),of an ’associated’rigid but massless foundation,as a function of the excitation frequency,ω(Gazetas,1983).

Fig.13.Compliance functions for foundation over half-plane withν=0.4 andξ=0.05.

Fig.13.(continued).

Fig.14.Geometry of foundation on a layer over bedrock(unit:m).

For each particular harmonic excitation with frequencyω,the dynamic impedance is defined as the ratio between the steadystate force(or moment)and the resulting displacement(or rotation)at the base of the massless foundation.For example,the vertical impedance of a foundationwhose plan is center-symmetric can be defined as

where Rv(t)=Rveiωtis the harmonic vertical force applied at the base of the disk,and V(t)=Veiωtis the uniform harmonic settlement of the soil-foundation interface.It is evident that Rvis the total soft reaction against the foundation;it is made up of the normal stresses against the basement plus,in case of embedded foundations,the shear stresses along the vertical side walls.

Similarly,one may define the horizontal impedances,Kh,from the horizontal forces and displacements along the principal axes of the base.Referring to Eq.(25),it is interesting to note that the dynamic force and displacement are generally out of phase.In fact,any dynamic displacement can be resolved into two components:one in phase and one 90°out of phase with the imposed harmonic load.It is convenient then to introduce complex notation to represent forces and displacements.Consequently,the impedances may also be written in the form of

The real and imaginary components are both functions of the vibrational frequency.The real component reflects the stiffness and inertia of the supporting soil;its dependence on frequency is attributed solely to the in fluence that frequency has on inertia,since soil properties are essentially frequency-independent.The imaginary component reflects the radiation and material damping of the system.The former,being the result of energy dissipated by waves propagating away from the foundation,is also frequencydependent;the latter,arising chie fly from the hysteretic cyclic behavior of soil,is practically frequency-independent.

5.2.Dynamic compliance functions

Also given the names of dynamic ’displacement’functions and dynamic ’ flexibility’functions,theyare essentiallythe ratiosbetween dynamic displacements(or rotations)and the dynamic reactive forces(or moments)at the base of a foundation.It is convenient to express compliance using complex notation(Gazetas,1983):

The real and imaginary parts represent the displacement components which are in-phase and 90°out of phase with the reactive force,respectively,and they both are functions of frequency,as discussed previously.

Table 3 Material properties of layer on bedrock.

Fig.15.(a)Vertical and(b)horizontal displacements of foundation on a layer over bedrock withν=0.4 andξ=0.05.

5.3.Application of the PML model

A force vibration of Ricker type has been applied in the centerline of the surface of the foundation(Fig.6).The analysis of foundation vibration by application of the PML is validated by evaluation of impedance/compliance for different soil conditions such as half-plane,layer over bedrock and layer over half-plane.Please note that the following results could be obtained for any region with any shear wave velocity,given that a small value for Vshas no effect on the response of the problem.

The purpose of the study is the application of PML in vibrating shallow foundations,mostly involving machine foundation analysis.Since the existing analysis and formulations in the literature are conducted for a rigid massless foundation,the foundation mass in the present modeling is considered equal to zero.In addition,the nodes placed on the foundation location are constrained to each other rigidly so that the displacements of the mentioned nodes are the same and they experience a rigid body motion.The foundation is assumed as a strip footing in 2D plane-strain modeling.

5.3.1.Rigid massless foundation on half-plane

The geometry of the rigid strip foundation over half-plane numerically investigated is presented in Fig.7.The width of the massless(m=0)foundation is B=2 m and the PML one is the halfwidth of the footing,LP=1 m.In the corresponding PML,the stretching functions(λi)are in the form of Eq.(6)with linear attenuation functions(slope of 10/LP).

The properties of the half-plane region are defined in Table 2.The maximum element size is 0.25 m,which is placed within 1/10 of the longest wavelength to provide accurate wave transmission,and the time step is equal to 0.004 s.For comparison purpose,the vertical and horizontal displacements of the rigid massless foundation for both extended mesh region and PML media are demonstrated in Fig.8.These deformations are produced by vertical and horizontal loading forces,respectively.

The results obtained from the PML model follow the extended mesh results closely,even though the domain is small enough for other kinds of absorbing boundaries to reflect waves back to the footing,as manifested in the higher response amplitudes.

In order to make a complete illustration of the PML efficiency,horizontal and vertical displacement contours in thewhole area are presented in Fig.9 for the moments of 0.2 s,0.28 s and 0.4 s.The mechanism of the wave propagation in the region shows that the forced vibrations produced by both of horizontal or vertical loadings can be fully absorbed and damped by the boundaries surrounded by PMLs.In other words,after 0.4 s of foundation vibration,all of the scattered waves through the medium are absorbed by PML and would damp without any reflection.The operation,in fact,is the role of radiation damping in FEM modeling.

Fig.16.Compliance functions of foundation on a layer over bedrock withν=0.4 andξ=0.05.

Fig.17.Geometry of foundation on a layer over half-plane(unit:m).

With no PML in previous cases,the vibration would be permanent in the medium and the usual static boundary conditions lead to wave reflection into the solution domain.A sample of model without PML is shown in Fig.10;in addition,the outputs are presented in Fig.11,comparing the foundation displacements in cases of applying PML and without it.

The results in the form of impedance functions for vertical(kv,cv)and horizontal(kh,ch)oscillations and in form of compliance functions for vertical(Fv1,Fv2)and horizontal(Fh1,Fh2)oscillations are presented in Figs.12 and 13,respectively.It should be noted that the curves shown in Fig.13 are the responses of the foundation assuming a material damping of 5%for the medium of the wave propagation.

The impedance functions have been compared with the analytical results of Gazetas and Roesset(1986)and the compliance functions have been compared with the analytical results of Gazetas(1983).Using this small domain,the results obtained from the PML model are in good agreement with the previous studies;however,the existence of material damping in the main domain and in the PML makes more accurate responses.It is worth mentioning that the PML outputs are obtained at a low computational cost.In general,the trends of responses are the same as analytical studies and the little differences come from the method of solutions.

5.3.2.Rigid massless foundation on layer over bedrock

The second condition is a foundation on a layer over bedrock of which the geometry is illustrated in Fig.14.The dimensions of the massless foundation are B=2 m and the PML of the half-width of the footing,i.e.LP=1 m.The thickness of the layerover bedrock has been selected as H/B=1 to be comparable with the previous investigations.Moreover,the bedrock is modeled as a rigid base.In the corresponding PML,the stretching functions(λi)are in the form of Eq.(6)with f1(x1)=0 and linear f2(x2)(slope of 20/LP)in the PMLs.The properties of the medium are defined in Table 3.The maximumelement size is0.25 m,whichis placed within 1/10 of the longest wavelength to provide accurate wave transmission,and the time step is equal to 0.004 s.

Forcomparison,the vertical and horizontal displacements of the rigid massless foundation for both extended mesh region and the media with PML are demonstrated in Fig.15.These deformations are produced by vertical and horizontal forces,respectively.

The results obtained from the PML model follow the extended mesh results closely,even though the domain is small enough for other kinds of absorbing boundaries to reflect waves back to the footing.The results in the form of compliance functions for vertical(Fv1,Fv2)and horizontal(Fh1,Fh2)oscillations are presented in Fig.16.

The compliance functions have been compared with the analytical results of Gazetas(1983)and Huh and Schmid(1984).Using this small domain,the results obtained from the PML model match fairly well the previous studies.It is noted that the PML outputs are obtained at a low computational cost.In general,the trends of responses are the same as analytical studies and all of the peaks and valleys are exactly captured by the introduced model.

5.3.3.Rigid massless foundation on layer over half-plane

In the last condition,a foundation on a layer over half-plane is evaluated.The geometry is illustrated in Fig.17.The dimensions of the massless foundation are B=2 and the PML of the half-width of the footing,LP=1 m.The thickness of the layer over half-plane has been selected as H/B=1 to be comparable with the previous investigations.In the corresponding PML,the stretching functions(λi)are in the form of Eq.(6)with linear attenuation functions(slope of 10/LP)in the PMLs.The PMLs utilized for the layer and the half-plane have different moduli,corresponding to the moduli for the elastic media as defined in Table 4.The maximum element size is 0.25 m,which is placed within 1/10 of the longest wavelength to provide accurate wave transmission,and the time step is equal to 0.004 s.

The results show that displacements of the PML and the extended mesh model are similar,as shown in Fig.18.It is worth mentioning again that the domain is small enough for other kinds of absorbing boundaries to reflect waves back into the footing.Without PML inprevious case,the vibrationwould be permanent in the medium and the usual static boundary conditions lead to wavereflection into the solution domain.A sample of model without PML is shown in Fig.19.The outputs are presented in Fig.20,comparing the foundation displacements in cases of applying PML and without it.

Table 4 Material properties of layer on half-plane.

Fig.18.(a)Vertical and(b)horizontal displacements of foundation on a layer over half-plane withν=0.4 andξ=0.05.

The compliance functions for vertical and horizontal oscillations are presented in Fig.21,which have been compared with the analytical results of Gazetas(1983).Using this small domain,the results obtained from the PML model match fairly well the previous studies.It should be noted that the PML outputs are obtained at a low computational cost.In general,the trends of responses are the same as analytical studies and the slight differences come from the approach of calculating material damping which are of two different types of Rayleigh and viscous for the present study and the literature,respectively.

5.4.Comparing with physical modeling

Dobry et al.(1986)have presented the results of the studies by Stokoeand Erden (1974)in their paper as a case of comparison.Stokoe and Erden(1974)investigated the dynamic damping of the rectangular shallow foundations by physical modeling.In order to simulate absorbing boundary conditions,they utilized wet sawdust of which the properties have been obtained by resonant column test.The idealized vertical section of the apparatus along with the material properties is shown in Fig.22(βis the material damping ratio).The properties of the sand and sawdust differ drastically:the damping ratio(β)for the sawdust is 30%,which is much higher than that of the sand(2.5%),while VSandρvalues for the sawdust are much smaller than that of the sand.A part of the results obtained by these researchers are displayed in Fig.23.

In the present study,strip footing has been analyzed using the prepared code with PML;the results are comparable with the physical modeling of the rectangular foundation with L/B≥10.The material properties in the numerical modeling are the same as the physical modeling except for the geometry dimensions which are altered due to the use of the PML instead of sawdust layer.The specifications of the PML are assigned based on the former problems.The model geometry and discretization are shown in Fig.24.The output in the form of the dynamic damping(the equations described in the Dobry et al.(1986))is presented in Fig.25.

Fig.19.Geometry and discretization of layered model without PML(unit:m).

Fig.20.(a)Vertical and(b)horizontal displacements of foundation on layered model with and without PML.

Fig.21.Compliance functions of foundation on a layer over half-plane withν=0.4 andξ=0.05.

Fig.22.Cross-section of Stokoe and Erden test facility(Dobry et al.,1986).

Fig.23.Vertical radiation damping(Dobry et al.,1986).

As can be seen,the results of the finite element analysis using PML are in a good agreement with the physical modeling.Therefore,the considered dimensions and PML characteristics that have been used to model the vibrating shallow foundation are reliable and the proposed method can be applied in similar analyses.

Fig.24.Numerical model for simulating the physical experiment(unit:m).

Fig.25.Vertical damping calculated by numerical and physical modeling.

6.Conclusions

The key to properly model wave propagation using the FEM is implementing radiation damping condition.In this context,a new FEM program has been written which is able to analyze different geometries with any number of layers under dynamic or seismic loading in the time domain.The program is featured by implementing PML formulation to simulate the radiation damping.The PML formulation has been implemented in the code by rewriting the dynamic equation of motion and deriving new properties for the quadrilateral elements.This study focuses on the calibration of the transient dynamic analysis of soil-foundation by FE-PML formulation.The near- field part of the foundation is modeled by FEM surrounded by a PML which models the far- field by absorbing waves propagating towards in finity.

In order to verify the dynamic performance of the program,an example of a plane-strain model has been established to measure the compressive,Rayleigh and shear wavevelocities.The analysis of foundation vibration using the PML is validated by the evaluation of impedance/compliance for different soil conditions such as halfplane,layer over bedrock and layer over half-plane.In the same way,vibration of strip footing has been analyzed using PML and the results have been compared with the physical modeling.

The results obtained from the PML model follow the extended mesh results closely,even though the domain is small enough for other kinds of absorbing boundaries to reflect waves back to the footing.The mechanism of the wave propagation in the region shows that the forced vibrations produced by either horizontal or vertical loading can be fully absorbed and damped by the boundaries surrounded by PMLs.The operation,in fact,is the role of radiation damping in FEM modeling.

Using this small domain,the results obtained from the PML model match the previous studies.In general,the trends of responses are the same as analytical studies and the minor differences come from the method of solutions and the approach of calculating material damping which are of two different types of Rayleigh and viscous for the present study and the literature,respectively.

Compared with others,the adjusted model yields more accurate and stable solution in the case of wave propagation,independent of the frequency or angle of incidence.Moreover,a case of the PML application is presented,where the foundation response is evaluated for different soil/rock conditions.Finally,it is concluded that the PML offers a reliable and efficient approach for modeling radiation damping mechanism.The formulation is applicable to elastodynamic problems involving wave propagation into unbounded media.

Conflict of interest

The authors wish to confirm that there are no known conflicts of interest associated with this publication and there has been no significant financial support for this work that could have in fluenced its outcome.