Calculation of impedances for axial symmetric foundation embedded in half-space medium using solutions for one layer stratum

The paper is devoted to develop a new scheme to calculate impedance matrices for axial-symmetric foundations embedded in halfspace medium. The half-space medium can be approximated by increasing thickness of one layer stratum on rigid bedrock due to the effect of material damping. However, as thickness increases, the numerical problem will arise. This problem is caused by the numerical contamination by some negligible reflection waves from rigid bedrock. In reality, the effect of these reflection waves on impedance is getting small as the thickness of the one layer stratum increases. Therefore, the scheme will employ the solutions for one layer stratum with suppressing these reflection waves to generate the impedance for the case of half-space. The numerical results by the presented scheme are compared with the results by other scheme in order to show that the new numerical scheme is effective and the solutions in layered medium can be extended to obtain the results for the cases of layered half-space medium. Some numerical results of torsional, vertical, horizontal, coupling and rocking impedances with different embedded depths will be presented and comments on the numerical scheme will be given.


INTRODUCTION
For many decades, the vibration of foundation embedded in half-space or layered media has attracted some attentions.Therefore, there are many approaches to deal with this problem.In the early time, the finite element method has been popular to determine the impedance functions for foundation embedded in an elastic stratum.In the model, the soil medium can be divided into nearfield and far-field.For the near-field, regular meshing of finite element method is employed.For the far-field, infinite element is used to model the semi-infinite domain [1,2].In order to improve the deficiency of infinite element model, Gupta et al [3] and Tzong and Penzien [4] employed analytical approach to simulate the outgoing waves in the far-fields of half-space and layered media respectively.
Regarding boundary element method (BEM), Green's function in frequency domain for layered half-space medium developed by Luco and Apsel [5] is employed.The impedance functions for foundation embedded in layered half-space using the Green's function was generated [6,7].However, using boundary element method to compute impedance functions, one has to deal with singularity problem by avoiding coincidence between source and observation points.Also, cone model was developed to calculate the dynamic response of a disk on surface of a soil layer resting on flexible rock subjected to harmonic excitations [8].Furthermore, the concept of cone model was extended to calculate the dynamic stiffness for foundation embedded in a multiple-layered half-space [9].
Regarding the method using analytical solution, Liou has developed a technique to decompose boundary conditions in the way to match the boundary values of general solution of wave equations in cylindrical coordinates for the cases of layered media [10].The technique has been successfully applied to find the impedance functions for foundations on layered half-space medium [11,12] and axial symmetric foundation embedded in layered medium [13].For the cases of foundation embedded in half-space medium, one can approximate the cases by increasing the thickness of a layered medium, if material damping exists in the half-space medium.
However, some numerical problems arise in the process of using the analytical solution technique as the thickness of layer increases.This is due to magnitude difference in z-direction between the exponential terms representing variation of upward reflecting waves and the exponential terms representing the variation of downward waves.The magnitude difference between upward reflecting waves and downward waves becomes enormously huge for the modes with large real parts of wave numbers in z-direction.This will make the numerical results unstable in calculation process.In the real situation, the contribution to vibrations near the embedded foundation from downward propagating wave should be much more important than that from upward reflecting wave, since the material damping will cause the reflecting waves decay a lot as they reach the foundation.
To remedy the numerical problem stated above, the upward propagating wave for the modes with very large real parts of wave numbers in z-direction should be suppressed or neglected.The paper is devoted to deal with this problem.The procedure developed by Liou and Chung [13] for layered medium will be modified to simulate the case of layered half-space medium.
The numerical results show that considering only downward propagating waves for the modes with very large real parts of wave numbers in z-direction is satisfactory.Also, the thickness of one layer stratum needed to simulate half-space medium will be investigated for different damping ratios.

DERIVATIONS OF MODAL SHAPES FUNCTIONS FOR CONSIDERING DOWN-WARD WAVES ONLY
Impedance matrices for circular foundation embedded in layered medium was formulated by Liou and Chung [13].In formulating impedances, all the solutions (modal shapes) used at the boundaries of interior domain and exterior domain are summarized in Fig. 1.In the figure, vector P 2 (which is also defined in Fig. 2 of Reference 13) represents the nodal intensities of piecewise linear model for interaction tractions at surface S 2 , matrices N and G represent all modal shapes of displacements and tractions respectively at interfaces S 2 , S 3 and S 2 for interior or exterior domains, superscripts (i) and (e) of N and G indicate the functions defined in interior and exterior domains respective- ly, subscript h of N and G indicates homogeneous solution which satisfies the boundary condi- tions of free surface and rigid bedrock, subscript p indicates particular solutions which satisfies properly decomposed interactions tractions of P 2 on surface S 2 and rigid bedrock, vector α (i ) and α (e) represent the unknown participation factors of all the modes for interior and exterior domains respectively, a 0 is radius of cylindrical cavity, d is depth of cavity, and L is thickness of layer.By intuition, if material damping exists, it is possible to approximate impedance matrices for circular foundation embedded in half-space medium by the method of generating impedance matrices for the same foundation embedded in a very thick layer over rigid bedrock.Therefore, to use the method to find the impedance matrices with good precision, the thickness L of the layer must be large enough in order to make the upward reflection waves be damped out.However, as L increas- es, the magnitude difference between e ′ v L and e − ′ v L will be enormously huge for the modes with large real part of ′ v L , and the truncation error will become a big problem and contaminates the numerical results.This contamination will cause loss of significant figures in the process of calculation.
Since material damping exists and driving force is applied at the location of foundation, the upward propagating waves which reflect from bedrock and vary in vertical z -direction by the forms of e ′ v Z and e vZ attenuate exponentially.So, these upward propagating waves with large real part of ′ v L and vL can be suppressed in order to simulate the case of half-space medium.Therefore, the ex- pressions of displacements and tractions at the interfaces between interior and exterior domains, between interior domain and foundation, and between exterior domain and foundation for these modes must be re-derived by neglecting upward propagating waves.The followings will briefly give the revisions for the modes with large real part of ′ v L .Also, in the following equations, time har- monic variation e iωt has been cancelled at both sides of equations for convenience.
Sezawa have solved wave equations in cylindrical coordinates for a homogeneous half-space [14].In the following derivations, only downward propagating waves are taken into account.After some mathematical manipulation, the displacement field of each mode of n th Fourier component for exterior domain can be expressed as (1) where T is unknown coefficients and will be determined by boundary conditions, k is horizontal wave number of the mode, , C S is shear wave velocity and C P is compressional wave ve- locity, ω is frequency, H n (kr) is second kind of Hankel function of order n , and dH n (kr) dr , subscript n is n th Fourier component with respect to azimuth θ .For the interior domain, the expressions are the same as that for exterior domain above, expect Hankel functions H n (kr) in Eq.( 2) are replaced with Bessel functions J n (kr) for interior domain.
For the modes of particular solution for interior domain, the traction field on horizontal plane can be written as : (4) where G is shear modulus of the layer, and matrix J is similar to H of Eq. ( 2) except Hankel functions H n (kr) are replaced with Bessel functions J n (kr) .
The modes for exterior domain contain only homogeneous solutions which satisfy the homogeneous boundaries at free surface z = 0 .The unknown coefficients in vector A in Eq. ( 1) can be expressed in term of displacement vector at free surface (z = 0) u 0 (e) as follows: ) where u 0 (e) is equal to u(r, z) z=0 in Eq.( 1).Substituting Eq. ( 6) into (1) and letting r = a 0 , one obtains ( where According Reference 13, H −1 u 0 (e) can be expressed as 1 ξ i 0 ( ) Rayleigh modes and 0, 0 1 ( ) for Love modes.For the above expressions, α i (e) is unknown modal participation factor and ξ i is scale factor of the mode which is defined in Reference 13.Also, after some mathematical manipulation of using Eq. ( 7), the tractions at depth z on the vertical interface ( S 1 + S 3 in Fig. 1) for the exterior domain with only considering downward propagat- ing waves can be expressed as the follows: where For the interior domain shown in Fig. 1, the solution is the combination of particular solutions and homogeneous solutions.For the modes of homogeneous solution, the expression for displacement For the modes of particular solution, one can refer to Reference 13.The transformation matrix Q n (Eq.34 in Reference 13), which transform piecewise linear traction at z = d to displacement at z = d , can be revised and written explicitly as follows: The traction at the interface ( S 2 in Fig. 1) between interior and exterior domains for the modes of particular solution can be obtained by using Eq. ( 9) with replacement of Hankel functions with respective Bessel functions, e − ′ v z and e − vz with e − ′ v (z−d ) and e − v(z−d ) respectively, and u 0 (e) with u 0 (i) .The u 0 (i) can be expressed in terms of vector P 2 of nodal intensities as shown in Reference 13.
Therefore, as shown in Fig. 1, by summing all modes of homogeneous and particular solutions, the total displacements and tractions on surface S 2 of interior domain can be written as ∑ P 2 , and the total displacements and tractions on surface S 3 of interior domain u Similarly, the total displacement and traction on surface S 3 + S 1 of exterior domain can be expressed as u ∑ (z)α (e) and t The above expressions of Eqs.( 1) ~ ( 14) are only employed for the modes with large real part of ′ v L .For the modes with smaller real part of ′ v L , the expressions presented by Liou and Chung are employed [13].Then, compatibility condition, variation principle and reciprocal theorem are employed to generate the impedance functions as Reference 13 suggested.

NUMERICAL INVESTIGATIONS
By using desktop computer to do the floating point computation, quadruple precision is defined in computer program.After extensive study, the expressions (Eqs.1-14)derived in previous section will be employed if the real part of ′ v L for the modes is greater than 25.One should also note that real part of vL is always greater than that of ′ v L .For those modes with real part of ′ v L smaller than 25, the expressions in Reference 13 are employed.
In calculating the numerical results of impedance function, Poisson ratio of soil medium is selected to be 0.33 and damping ratios is assigned to be 0.02, 0.05 and 0.1.In Figs.2-8, all the impedance functions shown are non-dimensionalized by shear modulus G and foundation radius a 0 , and the frequencies are non-dimensionalized by a 0 and real part of complex shear wave velocity Re(C S ) , K TT is torsional impedance, K VV is vertical impedance, K HH is horizontal impedance, K RH = K HR are coupling impedances, K RR is rocking impedance, and ω is frequency.In order to find the ap- propriate thickness of one layer stratum for approximating half-space, the impedances for foundation on one layer stratum are employed to compare with that for foundation on half-space.Figs.2-3 show the comparisons.In these figures, damping ratio ξ = 0.05 is selected, the thickness L of layer is gradually increased from a 0 L = 0.5 to a 0 L = 0.1 , and the numerical results are compared to that of Liou's previous work [11,15].In general, one can observed that the results of impedance functions for the case of one layer stratum is approaching that for the case of half-space medium, as a 0 L is getting smaller.Also, the fluctuation of the result is getting less dramatic as the layer is getting thicker.The reason for the phenomenon is that the traveling path of reflection waves become longer as L increases.This longer path will make the energy loss of reflection waves from rigid bedrock (upward propagating waves) greater.Furthermore, if one compares Figs. 2 with Figs. 3, one can observe that the fluctuation in Figs. 2. is less severe and the torsional impedance for the case of one layer stratum approaches that for half-space medium more quickly as a 0 L is getting smaller.This is because only shear waves are involved in calculating torsional impedance and shear wave length is shorter than compressional wave length which dominates in vertical vibration of foundation.This means the reflection shear waves need shorter distance to damp out energy.Therefore, to simulate the case of half-space medium, thinner layer can be used for generating torsional impedance.Similar observations can also be found from the results of horizontal, coupling and rocking impedances which are not shown in the paper.From these two figures and other convergence study, one can conclude that a 0 L = 0.1 is small enough for one layer stratum to simulate half-space medium, if damping ratio is 0.05 and non-dimensionalized frequency ω a 0 2π Re(C S ) = 0 ~1.For the cases of damping ratios ξ = 0.02 and ξ = 0.1, a thorough investigation has been done just like the study for ξ = 0.05 .The behaviors of the cases with ξ = 0.02 and 0.1 are similar to those presented above except the thicknesses of one layer medium needed to simulate half-space are different for different damping ratio.In general, smaller damping ratio in the medium means thicker layer needed to simulate half-space medium.For examples, for the case with ξ = 0.02 , a 0 L must be not greater than 0.08 and a 0 L must be not greater than 0.15 for the case with ξ = 0.1.
Figs. 4~8 show the simulated numerical results of non-dimensionalized torsional, vertical, horizontal, coupling and rocking impedances respectively for circular foundation embedded in half-space medium with embedded depths , and 1 .In the figures, the results for ξ = 0.02 are calculated by using a 0 L = 0.08 , the results for ξ = 0.05 are calculated by using a 0 L = 0.1 and the results for ξ = 0.1 are calculated by using a 0 L = 0.15 .From these figures, one can observe that the impedance, in general, tends to be greater as the embedded depth increases.Also, higher damping will make the impedances a little higher.

CONCLUDING REMARKS
The procedure to calculate impedance functions for foundation embedded in layered medium can be extended to calculate that for the case of half-space medium by using large thickness of layered medium.Although increasing the thickness of layered medium will give rise to the numerical problem in the process of computation.One just needs to suppress the upward propagating waves for the modes with large real part of ′ v L .In the study, calculation with quadruple precision (about 32 significant figures) is defined in computer program.This leads to that the upward propagating waves for the modes with Re( ′ v L) ≥ 25 should be suppressed.Also, to simulate the case of half- space medium, the thickness of layer could be decreased, if the damping ratio is greater.For examples, a 0 L = 0.15 for ξ = 0.1, a 0 L = 0.1 for ξ = 0.05 , and a 0 L = 0.08 for ξ = 0.02 .

Figure 1
Figure 1 Solutions at interfaces for layered medium with cylindrical cavity.

Figure 2 0 = 0
Figure 2 Comparison of non-dimensionalized torsional impedance with Liou's results for d a 0 and traction are similar to those shown in Eqs.(1)~(13) except the exponential functions e − ′ v z and e − vz are replaced with e − ′ v (z−d ) and e − v(z−d ) , Hankel matrices H,H 1 ,H 2 are replaced with Bessel matrices J, J 1 , J 2 in which Hankel functions are replaced with Bessel functions, and u 0 (e) is replaced with u 0 (i ) which is displacement vector at z = d .

Figure 3 0 = 0 Figure 4
Figure 3 Comparison of non-dimensionalized vertical impedance with Liou's results for d a 0

Figure 5 Figure 6 Figure 7 Figure 8
Figure 5 Non-dimensionalized vertical impedance with different embedded depths