Implementation and verification of direct infinite elements in axisymmetric problems

: In the field of dynamic soil-structure interaction, the simulation of geotechnical massifs under non-vanishing time-averaged loads is often challenging, depending on factors such as load variations, the 3D model representation and the chosen model boundaries. This study aims to evaluate the effectiveness of the method developed recently by the authors under non-vanishing time-averaged loads. The evaluation involves comparing the results of a homogeneous axisymmetric model at two control points and the simulation of an axisymmetric bilayer. The aim is to demonstrate the robustness of the method in contrast to widely used dissipative devices in scenarios involving sudden loading. Using an implicit integration scheme in the time domain, the method is shown to be effective and correlations with reliable solutions and comparisons with common dissipative devices for sudden-load models are presented. This analysis leads to the conclusion that the proposed procedure adeptly addresses the challenges associated with the simulation of geotechnical masses under non-vanishing time-averaged loads


Introduction
When performing dynamic simulations of soil-structure interaction, it is crucial to carefully consider the choice of model to ensure accurate representation of the physical parameters.The use of a 3D model is often associated with higher computational costs, while the choice of a simplified model can lead to significant discrepancies in response and problematic solutions.However, the representation should also take into account the loading scenario, the characteristics of the medium and the ability of the model to dissipate additional energy at the boundaries.If adequate boundaries are not established to manage excessive motion, waves may be reflected back into the structure rather than propagating out of the model.Consequently, a differential equation based solution may result in inaccurate field measurements (Donida et al, 1988).
The boundary element method is classified as a rigorous limit method that does not necessitate a high level of domain discretisation.Nevertheless, the coupled stiffness matrix is fully populated and asymmetric.In the lower part, the coefficient matrix no longer retains its band structure (von Estorff & Prabucki 1990).On the other hand, the application of the standard viscous boundary through the finite element method is straightforward.
However, this boundary depends on the angle of incidence of the wave and does not take into account the static behaviour of the far field, which leads to asymptotic divergences.
The use of infinite elements for this model has largely been in the frequency domain, cloning techniques that require the calculation of the inverse fast Fourier transform (Bettess & Bettess 1991, Yun et al 2007, Seo et al 2007, and Kazakov 2010, 2012).However, the existence of various types of waves complicates this method, prompting some research in the time domain.Haggblad and Nordgreen (1987), Su and Wang (2013), and Edip et al. (2013) suggested employing an absorbing layer combined with an inverse formulation of the infinite element to enhance its performance.Bakhtaoui andChelghoum (2018 &2020) have further expanded this approach to a direct formulation utilizing 4 and 6 node elements (Bakhtaoui 2024).However, although the validity has been demonstrated in defined dynamic applications, it remains useful to explore this new concept in different dynamic situations.
The aim of this paper is to present and test a method, developed by Bakhtaoui and Chelghoum (2020), for modelling the far field on two axisymmetric geotechnical cases under different mesh conditions.The main objective is to evaluate the performance of this method in dynamic scenarios, such as sudden actions on a foundation, such as the snapping of a transmission tower cable or the effect of pile driving.
The paper revisits the formulation of the method, applied to a direct infinite 4-node element, and emphasises tests conducted to evaluate its performance through an implicit integration scheme in the time domain.These tests include non-vanishing time-averaged loads, reproducing reference dynamic scenarios (Simon & Randolph 1986, Shridhar & Chandrasekaran 1995).In addition, the paper extends previous research (Bakhtaoui and Chelghoum 2020) by comparing results at two control points for a homogeneous axisymmetric model.It further extends this analysis to simulate an axisymmetric bilayer, as presented in the works of Simon & Randolph (1986) or Shridhar & Chandrasekaran (1995).The results of these comparisons confirm the robustness of the proposed method with respect to several dissipative devices commonly used in models subjected to sudden loading.

Formulation of direct infinite elements
As defined by Bettess (1977Bettess ( , 1980)), the direct infinite element can be obtained from interpolation functions obtained from a standard 3-nodal finite element based on Lagrange polynomials defined for i=1 to n-1 and specific decay formulations.Therefore, according to the conventions of Figure ( 1), these include the following interpolation expressions N i (,): Which, when combined with decay functions fi() in exponential form, lead to the expressions of the displacement functions H i (,) of the direct infinite element Bettess 1980Bettess , 1984) ) : The decay function are defined in this case by: Local coordinates η and ξ are defined such as , while L determine the decrease ratio.

Stiffness and mass matrices for the direct infinite elements
The stiffness matrix [K e ] for an element is given by the virtual work theorem in the same way as for the finite elements by: [B]: deformation matrix.
t: is the thickness of the element Following the approach proposed by Hinton et al. (1976), the evaluation of the mass to be assigned to each node of the infinite element is based on the concept of the consistent mass matrix [Mc], while the numerical integrations of the expressions are performed by the quadrature Gauss-Laguerre (Pissanetzky 1983).

Implementation of the viscous boundary in the proposed method
Developed and implemented by Lysmer and Kuhlemeyer (1969), viscous boundaries are widely used to simulate soil-structure interaction to avoid reflection by wave induction.Based on the formulation of normal and shear stresses, they are considered in this model for use in axisymmetric cases (Figure 2).Their expressions at the nodes of the boundary are given before plan generation by  = −..  .̇ (5)  = −..  .̇ ρ is the density of the medium Written in following nodal form: In these expressions, [D*] is the matrix of nodal damping and [ ̇] is the normal and tangential vector of velocities.
Vp and Vs are the velocities of P (compression) and S (Shear) waves, respectively, while a and b are coefficients suggested by White et al (1977), for taking into account of the direction of the incident wave.To maximize the absorption rate of P and S wave of the viscous boundary we take a = b=1.
Based on the work of Bakhtaoui and Chelghoum (2020), the scheme proposes an absorbing layer in the form of damping to ensure the dissipation of the wave, while the elastic recovery is performed by the direct infinite element.At the absorbing interface (Figure 3), the wave is Where P is the projection matrix relating normal and tangential velocities to the global Cartesian velocities components, while [H] is the matrix of displacement interpolation functions of interface finite element.The damping elements are calculated using a 2 points Gauss-Legendre numerical integration on each concerned side.Using diagonal damping matrix form, similar method basing on Hinton scheme (Hinton et al 1976) can be used.Hence, each type of dissipative system can be implicitly combined in the formulation of the global equations.

Applications
The above procedure is used to model a mass of soil subjected to a non-vanishing time-averaged load.This is carried out by original software developed in the FORTRAN language.In this test, the region of interest is represented by finite elements, while the far field domain is discretised using the direct-absorbing infinite elements described above.The masses and rigidities are evaluated using Gaussian numerical integration, while the resolution is performed using a step-by-step implicit integration algorithm.
The present applications consist in evaluating the response of geological masses with flexible bases, subjected to vertical loads with a non-zero time average.These applications are directly derived from the work of Simons and Randolph (1986), Wolf (1988), Shridhar and Chandrasekaran (1995).The use of these applications makes it possible to test the method by comparing the results obtained with those given in these references.This type of loading can be encountered in some dynamic problems, such as the foundation for a transmission line tower when the cables snap on one side, or pile driving analysis.For both selected applications, the load is applied suddenly at time t=0 and then held static.According to the Fourier transform, the load (Heavyside) gives low frequency components.The corresponding wavelengths are large, in which case the element models are suitable for solving the problem.It is worthy to note, in these examples, the parameters were given without units, in order to be consistent with references examples.

First Application: homogeneous elastic half-space case
A sudden vertical load is applied over a circular area of unit radius on the surface of a homogeneous linear elastic halfspace and then is maintained constant thereafter.This pulse could, representing the effect of a sudden load on an elastic infinite medium, can be evaluated by the analytic displacement solution given by Timoshenko and Goodier (1951): Where q is the unit load applied on the circular area.r is the unit radius of this disk, while E and υ are respectively the young and Poisson modulus.The values of these parameters are given on figure 4. The solution can also be checked by using a reference model, which can be obtained by extending the mesh elements with rigid boundary conditions.So that the wave after hitting the wall of reflection will not have time to disrupt the solution.The extended model used in this first application is represented by a 30x30 mesh (Figure 4 This gives a constant static displacement δ=0.75

Results
From the vertical displacement at point A (Figure 5) given by the extended model, the response reaches a maximum value and then stabilises at the static displacement value, which is considered to be the real solution to the problem.Compared to the extended model, the discrepancy shown when using the rigid model (on 9X9 mesh) is clearly visible from the 12th second and the solution becomes distorted thereafter.This is also the case for the response given by ordinary infinite elements, which oscillates around the static solution with significant peak differences from the 22nd second.This highlights the negative behaviour of static boundaries in dynamic problems.This trend can also be seen in Figure 6, where we can observe the evolution of the divergences through the responses at point B located at x=4m, and y=4m.The different curves obtained confirm the stability of the solution by the proposed procedure with respect to the extended model, while the response of the viscous boundary diverges asymptotically.
Using Simon and Randolf (1986) and Shridhar and Chandrasekaran (1995) results, the new solution based on a refined mesh (in order to approach the 8 noded elements used by authors), shows also good performance compared to those obtained by viscous boundary and extrapolation algorithm (Figure 7).

Second Application: Response of a two-layered medium
Used by Shridhar & Chandrasekaran (1995) for validation testing's on the extrapolation algorithm used on the boundaries, this application considers a two-layered medium subjected to the same sudden load applied by a circular area of unit-radius as previously defined.In order to approach the 8 noded elements used by authors, the finite element size is taken as Δx=Δz=0.5, that gives a 60x60 meshing for the reference model.The testing models used and geotechnical characteristics are given in Figure 8.

Results
Highlighted by results of the previous application, solutions obtained by rigid boundaries and infinite elements without absorbing layers present divergences more important than which obtained with other presented procedures (Figure 9).It is interesting to note that, in the absence of

Conclusions
This paper presents, implements and simulates a specific boundary for the simulation of axisymmetric geological massifs subjected to dynamic non-vanishing time-averaged loads in time domain.The models, represented by finite elements in the near field, employ direct infinite elements coupled with absorbing layers in the far field.In other words, from the viscous formulation, the impedances resulting from the passage of the wavefront are added to the direct infinite elements in order to ensure dissipation, while the latter is responsible for elastic recovery.
Simulations tested on the two cases of wave propagation, lead to the following conclusions: 1) The procedure is easy to implement using the classical finite element method.
2) The procedure, allowing a reduction in the number of elements and nodes, and can be applied to the solution of axisymmetric 3D problems according to a study window that can be determined by cutting off the oscillation axis.

Figure 1 .
Figure 1.Local and global representation of an infinite element

Figure 2 .Figure 3 .
Figure 2. planar viscous boundary .a), while the model used for testing consists of a 9x9 mesh (Figure 4.b) on which various boundary conditions will be tested, such as rigid condition, viscous boundaries, infinite boundaries and the proposed procedure (Figure 4.c).The finite elements or infinite elements size, composing the mesh, Δx=Δz are taken equal to unity.With a propagation velocity Vp = 1.732 m/s and Vs = 1m/s.The time integration is carried out with Δt=0.025 s.The complementary data used are υ = 0.25, E = 2.5, Density = 1, and Thickness = 1.

Figure 4 .Figure 5 .
Figure 4. (a) The large model, (b) The used model, (c) model of method proposed

Figure 6 .Figure 7 .
Figure 6.Comparison of verticals responses at point B