A quadratic boundary element for 3D elastodynamics

Abstract This article presents novel non-singular influence functions for homogeneous media. These solutions are displacement and stress fields of a three-dimensional, isotropic full-space under time-harmonic vertical and horizontal loads, which can be used within the framework of boundary element methods to solve elastodynamics problems in engineering practice. In order to account for sharply-varying contact tractions that may occur in such problems, the solutions in this article consider a biquadratic distribution of the loads within the loaded surface. In the present derivation, sets of Fourier transforms are used to uncouple the medium's equation of motion and enable the incorporation of boundary conditions directly as traction discontinuities. The article brings selected numerical results for various geometric and constitutive parameters. Graphical Abstract


INTRODUCTION
Boundary element methods are typical choices of methods to model unbounded media such as the soil.This success is partly due to their ability to comply with radiation conditions of unbounded media (Sommerfeld, 1949) and because they do not require discretization of the unbounded domain (Kuzuoglu and Mittra, 1997).The Direct Boundary Element Method (DBEM) is one example of such methods.In DBEM formulations, a discretized boundary integral equation is used to approximate the solution of the problem.In elastodynamic problems, for example, this equation relates displacement and stress Green's functions, which comprise the displacement fields and stress tensors of the medium of interest in response to point loads.
A comprehensive literature review of Green's functions for elasticity has been presented by Pan (2019).A wellknown difficulty in DBEM schemes is dealing with singularities resulting from collocation of Green's functions at the boundary of the domain (Dumont, 1994).One way to avoid this difficulty is by using non-singular influence functions, which can be collocated at the boundary without resulting in singularity problems.These functions can be used in Indirect-BEM (IBEM) schemes, in which displacement and stress influence functions are related through sets of fictitious loads, rather than by a boundary integral equation.Some of these non-singular influence functions have been derived by the authors of this article and their work group in the past decades.These solutions include isotropic (Mesquita et al., 2012b) and anisotropic media (Barros and Mesquita, 1999), harmonic and transient excitations (Mesquita et al., 2003), concentrated (Mesquita et al., 2009a) and distributed loads, for half-spaces (Mesquita et al., 2012a) and full-spaces (Labaki et al., 2019), and with a variety of different methods of derivation (Adolph et al., 2007;Romanini et al., 2019).Extensive investigations into the numerical evaluation of such solutions (Labaki et al., 2012) and strategies to deal with their high computational cost (Mesquita et al., 2009b) have been presented as well.In view of soil media applications, other groups have also invested significant effort to model influence functions that represent the soil's transversely isotropic, layered, and poroelastic characteristics.Selected examples of these functions, together with their applications within IBEM frameworks, are models of the dynamic response of rectangular plates embedded in layered, transversely isotropic soils (Fu et al., 2017;Fu et al., 2019), models of the response of strip foundations on transversely isotropic, layered, elastic (Ba et al., 2018a) and poroelastic (Ba et al., 2018b) soils, and models of the seismic response of alluvial basins (Ba et al., 2020).
Solutions for uniformly-distributed loading cases such as those presented by Romanini et al. (2019), however, face a difficult challenge when used to model discontinuous contact problems.Such problems are common in multi-media and multi-body interaction applications, and the challenge is to represent sharply-varying contact tractions that occur at the edges and interfaces between different bodies and media (Barros and Mesquita, 2000).Figure 1 illustrates this problem.It shows the horizontal and vertical contact tractions t X and t Z at the interface between a rigid strip footing and a half-space, resulting from the application of an external rocking moment M Y on the footing.These results were computed by Barros and Mesquita (2000) for the normalized frequency of excitation a 0 =ω/c s =1, in which c s is the shear wave speed in the half-space.Modeling quantities like these using uniformly-distributed loading solutions requires large numbers of elements and results in high computational cost, and in some problems is altogether unattainable (Barros and Mesquita, 2009).A contribution toward improved representation of sharply-varying contact quantities has been recently presented by Romanini et al. (2021), in which bilinearly-varying loadings were considered.The solution in this article is a significant improvement in that regard.This article presents displacement and stress solutions for a three-dimensional, isotropic full-space under timeharmonic vertical and horizontal excitation.External loads can have arbitrary biquadratic distribution over a rectangular patch within the full-space.Within DBEM and IBEM contexts, these non-singular influence functions correspond to quadratic boundary elements, which can be used in elastostatic and elastodynamic analyses with improved accuracy, convergence, and computational cost.The article presents original numerical results for various geometric and constitutive parameters, and shows that the presented implementation can be used within formulations of the boundary element method.

PROBLEM DEFINITION
Consider a three-dimensional, isotropic full-space, under horizontal (x-and y-directions) and vertical (z-direction) timeharmonic loads of circular frequency ω.In the absence of body forces, the equation of motion of this medium is given by in which μ and λ are Lamé's constants and ρ is the mass density of the medium, and u=u(x) is the displacement vector of point x=(x,y,z).The corresponding stress field σ=σ(x) can be obtained from u=u(x) through the constitutive relation in which δ ij is the Kroenecker Delta.
In this work, we derive solutions for both u and σ for the loading case illustrated in Figure 2: the loaded surface is a square patch defined by |x|<A; |y|<B; z=0, over which arbitrary, biquadratically-varying loads are applied.

SOLUTION PROCEDURE
The strategy to derive solutions for u and σ in this paper is based on the Helmholtz decomposition (Helmholtz, 1858).Each term u i of the displacement field u=u i ê i (i=x,y,z) is expressed as the linear superposition in which Δ and Ω are independent fields, and k L 2 =ω 2 ρ/(λ+2μ) and k S 2 =ω 2 ρ/μ are primary and secondary wave numbers of the full-space.Equation 2 can be rewritten in terms of the scalar dilatation field Δ and the vector rotation field Ω as Latin American Journal of Solids and Structures, 2023, 20(6), e508 4/26 in which n 2 =k L 2 /k S 2 .Trial solutions for Δ and Ω can be written in the wave number domain (β, γ) as in which the super-indices m=1,2 indicate m=1 (−∞<z≤0) and m=2 (0≤z<+∞) halves of the full-space.The solutions expressed in Eqs. 5 and 6 satisfy the condition that u must vanish for x→±∞ (Sommerfeld, 1949).
In view of the properties of Δ and Ω, Eqs. 5 and 6 yield The passage from Eqs. 5 and 6 to Eqs. 7 and 8 is detailed in Appendix A. Substituting Eqs. 5 and 6 into 3 and 4 yields Latin American Journal of Solids and Structures, 2023, 20(6), e508 5/26 in which A (m) and B n (m) (m=1,2; n=1,2,3) are arbitrary functions that depend on specific boundary conditions.

BOUNDARY CONDITIONS
Solutions for the full-space (−∞<z<+∞) are obtained from Eqs. 9 to 17 by imposing continuity and equilibrium conditions at the interface between media m=1 and m=2.The former is imposed for all displacement components u j (j=x,y,z) as , , 0 , , 0 , while the latter is imposed for all stress components σ Zj as , , 0 , , 0 , , in which P j (x,y)=p j e i(βx+γy) are external loads of amplitude p j in the j-direction (j=x,y,z) applied at the interface between media m=1 and m=2 (Fig. 2).The solution of all six algebraic equations involved in Eqs.18 and 19 results in A (m) and B n (m)  for the full-space boundary-value problem:

Biquadratic external loads
The external load distribution shown in Fig. 2, which is defined by its values at the nine points t i (i=1:9) can be described by direct biquadratic interpolation over p j (Q i )=t i : for|x|<A and |y|<B, and p j (x,y,z=0)=0 otherwise, in which

T t 
In order to allow its incorporation into Eq.19, Eq. 23 must be written in the transformed (β, γ) space, which is obtained through its direct Fourier transform with respect to (x, y): Latin American Journal of Solids and Structures, 2023, 20(6), e508 7/26 Substituting Eqs.20 to 22 into 9 to 17, together with Eq. 24, results in u(k β ,k γ ) and σ(k β ,k γ ).Applying the inverse Fourier transform into u(k β ,k γ ) and σ(k β ,k γ ) yields u(x) and σ(x), the displacement and stress fields in the physical domain.In order to improve the readability of this article, the final expressions for u(x) and σ(x) are listed in Appendix B. Two selected components are shown here in order to illustrate the general aspect of these solutions: the horizontal (x-direction) displacement due to vertical (z-direction) loads, and the shear stress in the y-z plane due to horizontal (y-direction) loads.The parameters involved in Eqs. 25 and 26 are listed in Appendix B as well.
Notice that these solutions are expressed in the [0,+∞) interval, rather than in the original Fourier transform (−∞,∞) interval.This can be obtained after careful consideration of whether the integrand of each term is an odd or even function.Moreover, the term z/|z| is incorporated into some components so that a single expression for each Latin American Journal of Solids and Structures, 2023, 20(6), e508 8/26 component can represent the entire full-space, rather than separate solutions for each media m=1,2.Both these modifications entail arduous mathematical manipulations and are aimed at facilitating the numerical evaluation of these expressions.

NUMERICAL RESULTS
A detailed description of the numerical scheme used in this work to evaluate Eqs.48 to 71 (Appendix B) is presented by the authors in Labaki et al. (2012).In summary, a combination of two routines is used: an adaptive quadrature scheme is used to integrate a finite region of the integrand containing singularities, while an extrapolation-based scheme is used to integrate the oscillatory-decaying remainder portion of the integrand (Piessens et al., 2012).Additionally, a small damping factor η=0.001 is incorporated into the constitutive parameters according to λ*=λ(1+iη) and μ*=μ(1+iη) (Christensen, 2010), which then replace μ and λ in Eq. 1.This causes the integration path to evade the real-axis singularities via a small detour through the complex plane (Michalski and Mosig, 2016).

Verification
Figure 3 shows a comparison of displacement components u YY and u ZX with the dynamic Kelvin solution (Kitahara, 2014).For this comparison, Kelvin's unit point-load solution was numerically integrated over a square 2A×2B area, and this case can be simulated with the present solution by making t 1:9 =1.These results were computed at point x=y=z=1.5, in terms of the normalized frequency a 0 =ω/c S =1, in which c S =√μ/ρ is the shear wave speed in the half-space.The rootmean-square errors (RMSE) between the present and the reference Kelvin solution are 2.43⋅10 −5 and 4.65⋅10 −6 for the real and imaginary parts of u YY , respectively (Fig. 3a), and 3.32⋅10 −5 and 2.97⋅10 −5 for the real and imaginary parts of u ZX , respectively (Fig. 3b).Since both solutions are synthetized numerically, they both have intrinsic numerical errors.These quantitative comparisons serve merely to verify that the solutions produce comparable results, rather than to validate them against some exact solution, since no such solution is available for this problem.As for the stress components, Figure 4a shows a comparison of selected components with the classical static, pointload Kelvin solution (Kane, 1994).For this comparison, we considered A=B=0.01,y=0.3, and z=0.5.Additionally, Fig. 4b shows a comparison with Barros and Mesquita's (1999) 2D plane-strain problem for various frequencies.In this Latin American Journal of Solids and Structures, 2023, 20(6), e508 9/26 comparison, t 1:9 =1, B/A=50, x=0.8, z=0.5, and y=0.The RMSE between the present and the reference solution is 8.47⋅10 −5 , 2.09⋅10 −4 and 1.34⋅10 −4 for σ YYZ , σ XXX and σ XXY , respectively (Fig. 4a), and 3.11⋅10 −4 and 3.24⋅10 −4 for the real and imaginary parts of σ XXZ , respectively (Fig. 4b).This shows that the present solution yields results that are comparable to the ones obtained by a variety of different methods.
Figure 5 shows a verification of the compliance of the stress components with the boundary conditions prescribed in Eq. 19.These results show that as the line in which these solutions are measured approaches the loaded surface (z=0), the stress components tend to the traction discontinuity prescribed at that surface.This is the static case, and results are measured on the x-z plane (y=0).Figures 5a and b consider respectively t i =i, and t 1:8 =0; t 9 =1.Finally, an additional verification of the displacement components is obtained by comparing stress components that were directly evaluated from the solutions in this article with stress components that were synthesized numerically from displacement components.These numerically-synthesized equivalents, indicated by the superscript D, are obtained by computing strain components through the numerical derivation of displacement components, and subsequent computation of stress components through Hooke's Law.This scheme is shown in detail in Romanini et al. (2021).Figure 6 shows selected comparisons for a 0 =2 and t i =i.Figures 6a and b consider respectively stresses along the line x=y=0.5;2≤z≤8 and y=z=0.5;1≤x≤5.The RMSE between the stress components evaluated directly from Eqs. 57 and 71 and the numerically-synthesized stress components is 7.72⋅10 −5 and 1.77 −4 for the real and imaginary parts of σ XYZ , respectively (Fig. 6a), and 5.61⋅10 −2 and 4.60⋅10 −2 for the real and imaginary parts of σ YZY , respectively (Fig. 6b).The considerably larger discrepancy between the two results in the case of σ YZY could be initially thought to be due to the increased difficulty in evaluating stress influence function near the loaded surface (x→A).However, we have shown in Fig. 5 that these solutions can be accurately computed even very close to the loaded surface.This indicates that the discrepancy in these results must come from numerical difficulties that are intrinsic to the finite difference scheme used to evaluate the numerically-synthesized solutions (Romanini et al., 2021).

Comparison between load distributions
Figures 8 to 11 show displacement and stress fields due to various load distributions.Cases A, B, and C consider t i =1/(4AB) (i=1:9) (constant load distribution), t i =3i/(80AB), and t 1:8 =0; t 9 =9/(16AB) (Figure 7).These were selected to have unit net magnitude, that is, Latin American Journal of Solids and Structures, 2023, 20(6), e508 10/26 so that the effect of the load distribution itself could be studied separately.All cases consider a 0 =1.These results show that displacement and stress fields are strongly dependent on the load distribution.Nonuniformly distributed loads result in much sharper variation of displacements than the uniformly-distributed load case.Stress states are more strongly dependent on the load distribution under the loading area (x,y<A), outside which the effect of load distribution quickly becomes negligible.As expected, the overall magnitude of the direct displacement and stress components u ZZ and σ XXX is larger than that of the cross displacement and stress components u XZ and σ XYZ , regardless of the load distribution, which is physically consistent.It is also physically consistent that the shear stress due to vertical load σ XYZ is nearly zero in the uniformly-distributed load case, due to the symmetry of the load and to the proximity of the measuring line (0<x<5; y=z=0.05) to the x-z plane.Figure 11 Effect of different load distributions on σ XYZ .

Quality of the implementation
Formulations of boundary element methods, the main applications of influence functions, often require such influence functions to be evaluated for relatively high frequencies and for field (measuring) points relatively far from the source (loaded) points.Figures 12 and 13 show selected displacement and stress components evaluated at distant field points, while Figure 14 and 15 show selected components evaluated at high frequencies.All results in this section consider t i =i.
Even though, due to the lack of comparable solutions in the literature, the results in this section cannot be verified in terms of accuracy, they are physically consistent.Some aspects showing this physical consistency are the overall decrease in the magnitude of the solutions for increasing distances from the source point and for increasing frequencies of excitation, and direct stress and displacement components posessing larger magnitude than their cross counterparts.The results are well-behaved, smooth curves, even for large values of x and a 0 .These results show that the strategy described in this paper for the computational implementation of these functions is able to handle a wide range of input data and parameters without failing, and that small changes in the input data do not lead to large changes in the output.

CONCLUSIONS
This article introduced original non-singular influence functions for three-dimensional, homogeneous, isotropic fullspaces.These solutions are aimed at facilitating the representation of sharply-varying quantities, such as contact tractions between soil and foundations, which routinely occur in soil-structure interaction problems.In order to represent these quantities, the solution considered biquadratically-distributed loads applied to the full-space, since loads that vary according to higher-order polynomials are better suited to represent sharply-varying quantities than uniformly-distributed loads.These time-harmonic external loads were imposed in the Fourier transformed domain, in which uncoupled stress and displacement components are accessible separately.Selected numerical results showed that the solution yields physically consistent results, and that the proposed numerical evaluation scheme yields viable solutions.This improved representation comes at the cost of lengthy solutions to be implemented.These influence functions can be used as quadratic boundary elements, with improved representation of sharply-varying quantities, for elastodynamic analyses.

APPENDIX A
This appendix expands the mathematical manipulation that yields Eqs.7 and 8 from Eqs. 5 and 6.The rotational Ω of the displacement field u is a vector given by for i,j,k=x,y,z.On the other hand, the dilatation Δ of the vector field u is a scalar given by Consider the following identity about the vector field Ω: since Ω 3,21 =Ω 3,12 , Ω 2,31 =Ω 2,13 , and Ω 1,32 =Ω 1,23 , which can be easily verified from Eq. 27.In view of Eq. 29, computing the gradient of the displacement field u (Eq. 3) yields Equation 5 expresses an Ansatz solution for Δ.Substituting 5 into 31 yields 1 e e e 0, , and the superindices (1,2) have been omitted from A (1,2) for conciseness.Since e q ≠0 and A (1,2) =0 is the trivial solution, Eq. 32 yields the portion of Eq. 7 referring to α L : Additionally, in view of Eq.
since u m,ijk =u m,ikj =u m,jik =u m,jki =u m,kij =u m,kji (m=1,2,3).In view of Eq. 34, computing the rotational of the displacement field u (Eq. 3) yields, for each of its components i, The left-hand side of Eq. 36 is (Graff, 2012) since e ijk u k is simply the rotational of the displacement field written in index notation (Eq.27).Omitting the term 1/k 2 S for conciseness, the rotational terms in the right-hand side of Eq. 36 can be expanded into e e e î ĵ k.
Consider initially the term of Eq. 38 in the î-direction.The following expansion can be made: In view of the terms Ω 1 , Ω 2 and Ω 3 of Ω expressed in Eq. 27, then Ω 1,11 =u 3,211 −u 2,311 , Ω 2,12 =u 1,312 −u 3,112 , and Ω 3,13 =u 2,113 −u 1,213 .Therefore, The same expansions can be made for the terms of Eq. 38 in the other directions.In view of these expansions, Eq. 38 yields .

APPENDIX B
This appendix lists the final expressions for the displacement and stress fields in the physical domain.The displacement components are given by and Latin American Journal of Solids and Structures, 2023, 20(6), e508 19/26 due to loads in the z-direction, due to loads in the y-direction, and due to loads in the x-direction, with u ZX =u XZ , u ZY =u YZ , and u YX =u XY , in which Latin American Journal of Solids and Structures, 2023, 20(6), e508 20/26 , , e , The corresponding stress fields in the physical domain are: Latin American Journal of Solids and Structures, 2023, 20(6), e508 21/26 A quadratic boundary element for 3D elastodynamics Edivaldo Romanini et al.
Latin American Journal of Solids and Structures, 2023, 20(6), e508 22/26 and due to loads in the z-direction, Latin American Journal of Solids and Structures, 2023, 20(6), e508 23/26 and due to loads in the x-direction, and Latin American Journal of Solids and Structures, 2023, 20(6), e508 25/26 due to loads in the y-direction, in which

Figure 1
Figure1Sharply-varying contact traction distribution at the interface between a half-space and rigid strip footing under rocking moment.

Figure 2
Figure 2 Biquadratically-varying load distribution within the full-space.

Figure 3
Figure 3 Comparison of selected displacement components with the Kelvin solution.

Figure 4
Figure 4 Comparison of selected stress components with (a) Kelvin problem solution and (b) Barros and Mesquita (1999) 2D problem.

Figure 5
Figure 5 Compliance of the stress solution with prescribed boundary conditions.

Figure 6
Figure 6 Comparison between stress components and displacement-based stress components.

Figure 7
Figure 7 Load distributions considered in this section.

Figure 8
Figure 8 Effect of different load distributions on u XZ .

Figure 9 26 Figure 10
Figure 9 Effect of different load distributions on u ZZ .

Figure 12 26 Figure 13
Figure 12Selected results from the evaluation of far-field displacement components.

Figure 14
Figure 14Selected results from the evaluation of high-frequency displacement components.

Figure 15
Figure 15Selected results from the evaluation of high-frequency stress components.
28, one can show that