# Seabed Dynamic Response of Offshore Wind Turbine Foundation under Vertical Harmonic Loading: An Analytic Solution

Authors: Dimitrios G. Pavlou; Yuzhu Li
Journal: Mathematical Problems in Engineering
 (2018)

---

## Abstract

Wave-induced dynamic loads yield displacement and stress propagation on the seabed beneath the offshore wind turbine foundation. Existing attempts to calculate the displacement/stress propagation based on commercial software packages or numerical solvers usually ignore the inertial effects and provide quasi-static solution for the dynamic problem or ignore the damping effect of the soil. Since real time displacement measurements on the seabed are difficult to be carried out, the existing results are not able to be validated and the conclusions are rather arbitrary. Analytic solutions can provide accurate results and are less time-consuming and more computationally efficient than the numerical approaches but difficult to be achieved. In the present study an exact solution of the wave propagation in the seabed under harmonic excitation of the wind turbine foundation is achieved. The solution is based on integral transforms and stress functions. Verification of the solution is provided and a representative example is solved and commended.


---

## Body

## 1. Introduction

Existing attempts to calculate the displacement and stress propagation on the seabed beneath the offshore wind turbine foundation are mainly based on commercial software packages or numerical solvers [1–4]. In the most existing works, the solutions to the wave, structure, and seabed interaction problems are based on a ‘one-way coupling’ assumption (e.g., [1, 3, 4]). In [3–5], the Biot’s equation [6, 7] is adopted to model the soil behaviour. The inertial effects are ignored, yielding quasi-static solution for the dynamic displacement propagation in the soil. Moreover, the assumption of the elastic behaviour of the soil does not take into account the damping effect. Since real-time displacement measurements on the seabed beneath the foundation are difficult to be carried out, the existing numerical results and the conclusions are not validated. Also, the usage of the numerical tools can have a high computational cost and is time-consuming for solving the coupling problem in order to achieve a low residual.Analytic solutions to such problems are more accurate and more efficient compared to the numerical ways but are difficult to achieve. A fundamental contact mechanics solution based on Green’s functions and boundary integral equation solution has been achieved in [8] for the stress analysis of the soil-foundation interaction under static torsion. Analytical solutions for dissimilar or homogeneous elastic medium (space and half-space) containing interfacial cracks, under axial, radial, and torsional loading sources are provided in [9–16]. Using double integral transforms, the dual integral equations for wave propagation of a cracked medium under transient dynamic loading are originally solved in [17]. Analytical solutions for the dynamic response of foundation on elastic and viscoelastic soil are achieved in [18–20]. Interaction of foundation piles embedded in transversely isotropic soil under dynamic loading is analysed in [21, 22]. Structural analysis of wind turbine tower under environmental loads is performed in [23], and analytical solution for evaluating the dynamic response of plates under underwater transient loads is presented in [24].In the present study, an exact solution of the wave propagation in the seabed under vertical harmonic excitation of the wind turbine foundation is achieved. The solution is based on integral transforms and stress functions. A representative example is solved and commended in the present work. Verification of the solution is provided by solving a specific case of a sudden point load on the solid surface and the solution is verified against the existing solution in [25]. The proposed solution has two benefits: (a) it takes into account the inertial effects in the displacement propagation, and (b) the solution is analytic, yielding better accuracy and negligible computation cost. The limitation of the present version is that it does not include the effect of porous material and damping yet. The analytic solution in the present work can be extended for solving porous media with time-varying dynamic loads taking into account the material damping effect.


## 2. Formulation of the Problem

A circular foundation based on an elastic half-space (e.g., rock seabed) is subjected to vertical harmonic load. The coordinate system of the foundation is shown in Figure1. A pressure distribution(1)σzzr,0,t=ZHr-Hr-asin⁡ωtis assumed to act on the surface of the half-space. The symbol Hr denotes the Heaviside function. The term Hr-Hr-a represents the uniform pressure distribution in the area 0≤r≤α (Figure 1(b)), the term sin⁡ωt is the harmonic function of time, and Z is the amplitude of the load. The harmonic term sin⁡ωt causes wave propagation in the domain of the elastic seabed.Figure 1
Circular foundation on elastic seabed: (a) physical model and (b) mechanical model.
(a) (b)Since no body-forces exist, the governing equation in terms of displacements (Navier’s equation) can be written [26] as(2)λ+μ∇∇∙u→+μ∇2u→=ρ∂2u→∂t2where u→ is the vector of the displacement components ur,uz,uϑ in the directions r,z,ϑ respectively, λ,μ are the Lame constants, and ρ is the material density. Taking into account the axisymmetric geometry and load, the displacement component uϑ is zero. Therefore, the displacement vector can be expressed by the following equation:(3)u→r,z,t=ure→r+uze→zwhere e→r,e→z are the unit vectors in the directions r,z respectively. In terms of scalar and vector potentials φ,h→, (3) can be written [25] as(4)u→=∇φ+∇×hϑe→ϑWith the aid of (4), (2) yields(5)∇2φ=1c12∂2φ∂t2(6)∇2hϑ-1r2hϑ=1c22∂2hϑ∂t2In above equations the parameters c1,c2 are given in [25] by the formulae(7)c1=λ+2μρ(8)c2=μρand express the propagation velocity of the dilatational and rotational waves, respectively. The operator ∇2 in polar coordinates is given by the following equation:(9)∇2=1r∂r∂/∂r∂r+∂2∂z2If we use [25] a function ψ, where(10)hϑ=-∂ψ∂r (6) can further be simplified:(11)∇2ψ=1c22∂2ψ∂t2With the aid of (4), (5), (10), and (11), the stresses (in polar coordinates) of Appendix A and the displacements can now be written in terms of the scalar functions φ,ψ.(a) Displacements (12) u r = ∂ φ ∂ r + ∂ 2 ψ ∂ r ∂ z (13) u ϑ = 0 (14) u z = ∂ φ ∂ z + ∂ 2 ψ ∂ z 2 - 1 c 2 2 ∂ 2 ψ ∂ t 2(b) Stresses (15) σ r r = λ ∇ 2 φ + 2 μ ∂ 2 φ ∂ r 2 + ∂ 3 ψ ∂ r 2 ∂ z (16) σ ϑ ϑ = λ ∇ 2 φ + 2 μ r ∂ φ ∂ r + ∂ 2 ψ ∂ r ∂ z (17) σ z z = λ ∇ 2 φ + 2 μ ∂ ∂ z ∂ φ ∂ z + ∂ 2 ψ ∂ z 2 - 1 c 2 2 ∂ 2 ψ ∂ t 2 (18) σ r z = μ ∂ ∂ r 2 ∂ φ ∂ z + 2 ∂ 2 ψ ∂ z 2 - 1 c 2 2 ∂ 2 ψ ∂ t 2 (19) σ ϑ z = 0 (20) σ r ϑ = 0


## 3. Solution of the Dynamic Problem

For the solution of the problem, Laplace and Hankel integral transforms will be applied on the governing equations and boundary conditions. The definition and some important properties of the above transforms are given below.


### 3.1. Laplace Transform of Any Functionf(t)

(21) L f t ; t → p = F p = p ∫ 0 ∞ f t e - p t d tThe inverse form of Laplace transform is given by the following formula:(22)L-1Fp;p→t=ft=12πi∫a-i∞a+i∞Fppetpdp


### 3.2. Hankel Transform of Any Functiong(r)

(23) H n g r ; r → ξ = g - ξ = ∫ 0 ∞ g r r J n r ξ d rThe inverse form of Hankel transform is given by the following formula:(24)Hn-1g-ξ;ξ→r=gr=∫0∞g-ξξJnrξdξTaking the Laplace transform of the displacement and stress components given by (12)-(20), formulae (B.1)-(B.7) in Appendix B will be obtained.


### 3.3. Boundary Conditions

The surfacez=0 of the half-space is loaded by the stress σzz(r,0,t) given by (1). Since the surface is free of shear stresses, the following condition is valid:(25)σzrr,0,t=0Taking into account the definition of the Hankel and inverse Hankel transform given by (23) and (24), (1) can be written in the form(26)σzzr,0,t=ZH0H0-1σzzr,0,t;r→ξ;ξ→ror(27)σzzr,0,t=Zsin⁡ωt∫0∞ξJ0ξrdξ∫0∞Hr-Hr-axJ0xξdxWith the aid of the properties [27](28)H0Ha-x;x→ξ=∫0∞Ha-xxJ0xξdx=ξ-1aJ1aξ(29)H0H-x;x→ξ=0 (27) yields(30)σzzr,0,t=Zasin⁡ωt∫0∞J0ξrJ1aξdξTaking the Laplace transform of above equation, the following formula can be obtained:(31)Szzr,0,p=Lσzzr,0,t;t→p=Za∫0∞ωp2+ω2J0ξrJ1aξdξAs a consequence of the definition of the Heaviside function, the following equations can be obtained:(32)Hx=1-H-x(33)Hx-a=1-Ha-xCombination of (28) and (33) yields(34)∫0∞Hx-axJ0xξdx=1-ξ-1aJ1aξWith the aid of (29), (32), (33), and (34), (27) can be written as(35)σzzr,0,t=Zasin⁡ωt∫0∞ξJ0ξr1ξJ1aξdξUsing the notation(36)g-ξ=1ξJ1aξ (35) can be written:(37)σzzr,0,t=Zasin⁡ωt∫0∞g-ξξJ0ξrdξAccording to (24), above equation yields(38)σzzr,0,t=Zasin⁡ωtH0-1g-ξ;ξ→rTaking the Hankel transform of above equation, the following formula can be obtained:(39)σ-zzξ,0,t=Zasin⁡ωtg-ξor(40)σ-zzξ,0,t=Zasin⁡ωt1ξJ1aξThen, the Laplace transform of (40) yields(41)S-zzξ,0,p=Zaωp2+ω21ξJ1aξTaking the Laplace and Hankel transform to the boundary condition given by (25), the following result can be obtained:(42)T-zrξ,0,p=0Applying the Laplace transform in (17) and (18),(43)Szzr,z,p=Lσzzr,z,t;t→p=λ∇2Φ+2μ∂∂z∂Φ∂z+∂2Ψ∂z2-p2c22Ψ(44)Trzr,z,p=Lτrzr,z,t;t→p=μ∂∂r2∂Φ∂z+2∂2Ψ∂z2-p2c22ΨTaking the Hankel transform H0∙;r→ξ to above equations, the following formulae can be obtained:(45)S-zzξ,z,p=λh2Φ-+2μ∂2Φ-∂z2+∂3Ψ-∂z3-κ2∂Ψ-∂z(46)T-rzξ,z,p=-μξ2∂Φ-∂z+2∂2Ψ-∂z2-κ2Ψ-where(47)h=pc1(48)κ=pc2Combining (41) with (45) and (42) with (46), the following system in terms of Φ-,Ψ- can be obtained:(49)λh2Φ-+2μ∂2Φ-∂z2+∂3Ψ-∂z3-κ2∂Ψ-∂z=Zaωp2+ω21ξJ1aξ,z=0(50)-μξ2∂Φ-∂z+2∂2Ψ-∂z2-κ2Ψ-=0,z=0


### 3.4. Determination of the Functionsφ,ψ

In order to solve the system of (49) and (50), Laplace and Hankel transform will be applied to the governing equations (5) and (11). The Laplace transformation of (5) and (11) yields(51)∇2Φ-p2c12Φ=0(52)∇2Ψ-p2c22Ψ=0Taking the Hankel transform H0∙;r→ξ to above equations, the following ordinary differential equations can be obtained:(53)d2Φ-dz2-ξ2+h2Φ-=0(54)d2Ψ-dz2-ξ2+κ2Ψ-=0The solution of above equations is(55)Φ-=A1e-αz+B1eαz(56)Ψ-=A2e-βz+B2eβzwhere(57)α=ξ2+h2(58)β=ξ2+κ2For z→∞ the functions Φ-,Ψ- will be unbounded because of the terms B1eαz and B2eβz, respectively. Therefore, the following conditions should be valid:(59)B1=0(60)B2=0Then,(61)Φ-=A1e-αz(62)Ψ-=A2e-βzWith the aid of above two equations, the system of (49) and (50), for z=0, yields(63)λh2+2μα2A1-2μββ2-κ2A2=Zaωp2+ω21ξJ1aξ(64)-2αA1+2β2-κ2A2=0Combining (47) with (57) and (48) with (58), the following formulae can be obtained:(65)λh2+2μα2=μκ2+2ξ2(66)β2-κ2=ξ2(67)2β2-κ2=2ξ2+κ2With the aid of (57)-(67), the system of algebraic equations (63) and (64) can be written as(68)2ξ2+κ2-2βξ2-2α2ξ2+κ2A1A2=Zaμωp2+ω21ξJ1aξ0yielding(69)A1A2=Zaμωp2+ω2J1aξD1ξ2ξ2+κ22αwhere(70)D=2ξ2+κ2-4αβξ2With the aid of (69) and (70), the functions given by (61), (62) can now be obtained:(71)Φ-Ψ-=Zaμωp2+ω2J1aξD1ξ2ξ2+κ2e-αz2αe-βzSince all the displacement and stress components given by (12)-(20) are expressed in terms of the functions φ,ψ the inverse Hankel transform H0-1∙;ξ→r of above equation yields the solution of the problem.


### 3.5. Inversion of the Hankel Transform

Taking the inverse Hankel transformH0-1∙;ξ→r in (71), the following result can be obtained:(72)ΦΨ=H0-1Φ-Ψ-;ξ→r=Zaμωp2+ω2∫0∞J1aξD2ξ2+κ2e-αz2αe-βzJ0ξrdξTaking into account (72), the Laplace transformed displacements given in Appendix B can now be written as(73)Urr,z,p=-Zaμωp2+ω2∫0∞ξ2ξ2+κ2e-αz-2αβe-βzDJ1ξrJ1aξdξ(74)Uzr,z,p=Zaμωp2+ω2∫0∞-α2ξ2+κ2e-αz+2αβ2-κ2e-βzDJ0ξrJ1aξdξTaking the inverse Laplace transform of above equations, the final solution for the displacements ur(r,z,t),uz(r,z,t) can be obtained:(75)urr,z,t=L-1Urr,z,p;p→t(76)uzr,z,t=L-1Uzr,z,p;p→tAbove inverse Laplace transforms can be carried out numerically with the aid of the Stehfest [28] method.


## 3.1. Laplace Transform of Any Functionf(t)

(21) L f t ; t → p = F p = p ∫ 0 ∞ f t e - p t d tThe inverse form of Laplace transform is given by the following formula:(22)L-1Fp;p→t=ft=12πi∫a-i∞a+i∞Fppetpdp


## 3.2. Hankel Transform of Any Functiong(r)

(23) H n g r ; r → ξ = g - ξ = ∫ 0 ∞ g r r J n r ξ d rThe inverse form of Hankel transform is given by the following formula:(24)Hn-1g-ξ;ξ→r=gr=∫0∞g-ξξJnrξdξTaking the Laplace transform of the displacement and stress components given by (12)-(20), formulae (B.1)-(B.7) in Appendix B will be obtained.


## 3.3. Boundary Conditions

The surfacez=0 of the half-space is loaded by the stress σzz(r,0,t) given by (1). Since the surface is free of shear stresses, the following condition is valid:(25)σzrr,0,t=0Taking into account the definition of the Hankel and inverse Hankel transform given by (23) and (24), (1) can be written in the form(26)σzzr,0,t=ZH0H0-1σzzr,0,t;r→ξ;ξ→ror(27)σzzr,0,t=Zsin⁡ωt∫0∞ξJ0ξrdξ∫0∞Hr-Hr-axJ0xξdxWith the aid of the properties [27](28)H0Ha-x;x→ξ=∫0∞Ha-xxJ0xξdx=ξ-1aJ1aξ(29)H0H-x;x→ξ=0 (27) yields(30)σzzr,0,t=Zasin⁡ωt∫0∞J0ξrJ1aξdξTaking the Laplace transform of above equation, the following formula can be obtained:(31)Szzr,0,p=Lσzzr,0,t;t→p=Za∫0∞ωp2+ω2J0ξrJ1aξdξAs a consequence of the definition of the Heaviside function, the following equations can be obtained:(32)Hx=1-H-x(33)Hx-a=1-Ha-xCombination of (28) and (33) yields(34)∫0∞Hx-axJ0xξdx=1-ξ-1aJ1aξWith the aid of (29), (32), (33), and (34), (27) can be written as(35)σzzr,0,t=Zasin⁡ωt∫0∞ξJ0ξr1ξJ1aξdξUsing the notation(36)g-ξ=1ξJ1aξ (35) can be written:(37)σzzr,0,t=Zasin⁡ωt∫0∞g-ξξJ0ξrdξAccording to (24), above equation yields(38)σzzr,0,t=Zasin⁡ωtH0-1g-ξ;ξ→rTaking the Hankel transform of above equation, the following formula can be obtained:(39)σ-zzξ,0,t=Zasin⁡ωtg-ξor(40)σ-zzξ,0,t=Zasin⁡ωt1ξJ1aξThen, the Laplace transform of (40) yields(41)S-zzξ,0,p=Zaωp2+ω21ξJ1aξTaking the Laplace and Hankel transform to the boundary condition given by (25), the following result can be obtained:(42)T-zrξ,0,p=0Applying the Laplace transform in (17) and (18),(43)Szzr,z,p=Lσzzr,z,t;t→p=λ∇2Φ+2μ∂∂z∂Φ∂z+∂2Ψ∂z2-p2c22Ψ(44)Trzr,z,p=Lτrzr,z,t;t→p=μ∂∂r2∂Φ∂z+2∂2Ψ∂z2-p2c22ΨTaking the Hankel transform H0∙;r→ξ to above equations, the following formulae can be obtained:(45)S-zzξ,z,p=λh2Φ-+2μ∂2Φ-∂z2+∂3Ψ-∂z3-κ2∂Ψ-∂z(46)T-rzξ,z,p=-μξ2∂Φ-∂z+2∂2Ψ-∂z2-κ2Ψ-where(47)h=pc1(48)κ=pc2Combining (41) with (45) and (42) with (46), the following system in terms of Φ-,Ψ- can be obtained:(49)λh2Φ-+2μ∂2Φ-∂z2+∂3Ψ-∂z3-κ2∂Ψ-∂z=Zaωp2+ω21ξJ1aξ,z=0(50)-μξ2∂Φ-∂z+2∂2Ψ-∂z2-κ2Ψ-=0,z=0


## 3.4. Determination of the Functionsφ,ψ

In order to solve the system of (49) and (50), Laplace and Hankel transform will be applied to the governing equations (5) and (11). The Laplace transformation of (5) and (11) yields(51)∇2Φ-p2c12Φ=0(52)∇2Ψ-p2c22Ψ=0Taking the Hankel transform H0∙;r→ξ to above equations, the following ordinary differential equations can be obtained:(53)d2Φ-dz2-ξ2+h2Φ-=0(54)d2Ψ-dz2-ξ2+κ2Ψ-=0The solution of above equations is(55)Φ-=A1e-αz+B1eαz(56)Ψ-=A2e-βz+B2eβzwhere(57)α=ξ2+h2(58)β=ξ2+κ2For z→∞ the functions Φ-,Ψ- will be unbounded because of the terms B1eαz and B2eβz, respectively. Therefore, the following conditions should be valid:(59)B1=0(60)B2=0Then,(61)Φ-=A1e-αz(62)Ψ-=A2e-βzWith the aid of above two equations, the system of (49) and (50), for z=0, yields(63)λh2+2μα2A1-2μββ2-κ2A2=Zaωp2+ω21ξJ1aξ(64)-2αA1+2β2-κ2A2=0Combining (47) with (57) and (48) with (58), the following formulae can be obtained:(65)λh2+2μα2=μκ2+2ξ2(66)β2-κ2=ξ2(67)2β2-κ2=2ξ2+κ2With the aid of (57)-(67), the system of algebraic equations (63) and (64) can be written as(68)2ξ2+κ2-2βξ2-2α2ξ2+κ2A1A2=Zaμωp2+ω21ξJ1aξ0yielding(69)A1A2=Zaμωp2+ω2J1aξD1ξ2ξ2+κ22αwhere(70)D=2ξ2+κ2-4αβξ2With the aid of (69) and (70), the functions given by (61), (62) can now be obtained:(71)Φ-Ψ-=Zaμωp2+ω2J1aξD1ξ2ξ2+κ2e-αz2αe-βzSince all the displacement and stress components given by (12)-(20) are expressed in terms of the functions φ,ψ the inverse Hankel transform H0-1∙;ξ→r of above equation yields the solution of the problem.


## 3.5. Inversion of the Hankel Transform

Taking the inverse Hankel transformH0-1∙;ξ→r in (71), the following result can be obtained:(72)ΦΨ=H0-1Φ-Ψ-;ξ→r=Zaμωp2+ω2∫0∞J1aξD2ξ2+κ2e-αz2αe-βzJ0ξrdξTaking into account (72), the Laplace transformed displacements given in Appendix B can now be written as(73)Urr,z,p=-Zaμωp2+ω2∫0∞ξ2ξ2+κ2e-αz-2αβe-βzDJ1ξrJ1aξdξ(74)Uzr,z,p=Zaμωp2+ω2∫0∞-α2ξ2+κ2e-αz+2αβ2-κ2e-βzDJ0ξrJ1aξdξTaking the inverse Laplace transform of above equations, the final solution for the displacements ur(r,z,t),uz(r,z,t) can be obtained:(75)urr,z,t=L-1Urr,z,p;p→t(76)uzr,z,t=L-1Uzr,z,p;p→tAbove inverse Laplace transforms can be carried out numerically with the aid of the Stehfest [28] method.


## 4. Numerical Example and Verification of the Solution

### 4.1. Numerical Example

The derived solution given by (75) and (76) can now be implemented to a representative example. The dynamic load is a uniform pressure acting on a circular surface with radius a. The load can be modeled by a harmonic function and can be represented in terms of the variables r,t by (1). The selected data for the proposed loading function are Z=1.0N,a=8.2m,ω=0.94rad/s, and the material properties for the soil are c1=355m/s,c2=309m/s, μ=2.137×108N/m2. For the inverse Laplace transforms of (75) and (76), the Stehfest [28] method is adopted. According to this method, the functions ur(r,z,t),uz(r,z,t) can be obtained by the following formula:(77)urr,z,tuzr,z,t≈ln⁡2t∑n=1NcnUrr,z,nln⁡2tUzr,z,nln⁡2twhere(78)cn=-1n+N/2∑k=n+1/2minn,N/2kN/22k!N/2-k!k!k-1!n-k!2k-n!In above equation, N is a free even number representing the number of terms used in the summation. Values within the range 10≤N≤24 yield sufficient accuracy (e.g., [29]). In the present study we adopted the value N=24.Forr=0, the Bessel function J1(ξr) is null. Therefore, (73) yields Ur=0. With the aid of the standard software “Mathematica” [30] the displacement distribution uz versus z for r=0 and t=0.3,0.5,0.7s is demonstrated in Figures 2(a), 2(b), and 2(c).Figure 2
Displacement distributionuz versus z for r=0 and t=0.3,0.5,0.7s.
(a) (b) (c)According to these findings, the displacement tends to reach the z-axis asymptotically due to the termse-az and e-βz in (74).


### 4.2. Verification of the Solution

The derived solution can be adjusted to provide the surface displacement distribution for the specific case of a suddenly applied point load given by the following equation:(79)qr,0,t=Pδr2πrHtwhere δ(r) is the Dirac Delta function of the variable r and Ht is the Heaviside step function of time t. Solution for this specific case is available in [25]:(80)Urr,z,p=P2πμp∫0∞ξ22ξ2+κ2e-αz-2αβe-βzDJ1ξrdξ(81)Uzr,z,p=P2πμp∫0∞-αξ2ξ2+κ2e-αz+2αξ3e-βzDJ0ξrdξFollowing the procedure of [14], the point load P is given by the following equation:(82)P=Zπa2Therefore,(83)Z=Pπa2The above term for the uniform distributed load Z should be placed in (73) and (74). On the other hand, the Laplace transform(84)ωω2+p2=Lsin⁡ωt;t→pshould be replaced in (73) and (74) by the Laplace transform:(85)1p=LHt;t→pFor the case of point load P, the radius a of the circular surface of the equivalent uniform load Z given by (83) is a→0. Therefore, (73) and (74) yield(86)Urr,z,p=lima→0-P/πa2aμ1p∫0∞ξ2ξ2+κ2e-αz-2αβe-βzDJ1ξrJ1aξdξ(87)Uzr,z,p=lima→0P/πa2aμ1p∫0∞-α2ξ2+κ2e-αz+2αβ2-κ2e-βzDJ0ξrJ1aξdξSince(88)lima→0J1aξa=ξ2 (86) and (87) yield the known [25] (80) and (81). With the aid of the Laplace inversion technique proposed in [25], (87) yields the wave propagation uz on the soil surface for the pulse point load of (79):(89)uzr,0,t=0forc2tr<13uzr,0,t=P32μπr6-3c2t/r2-0.25-33+53+3/4-c2t/r2+33-5-3+3/4+c2t/r2for13<c2tr<1uzr,0,t=P32μπr6-33+53+3/4-c2t/r2for1<c2tr<3+32uzr,0,t=3P8μπrforc2tr>3+32For a unit load P=1, the displacement distributions at time instants t=0.0005s, 0.001s, 0.0015s, and 0.002s are demonstrated in Figures 3(a), 3(b), 3(c), and 3(d).Figure 3
3D displacement distribution ofuz due to transient point load for various time instants. Implementation of the proposed solution to transient point load shows full agreement with the existing solution in [25].
(a) (b) (c) (d)


## 4.1. Numerical Example

The derived solution given by (75) and (76) can now be implemented to a representative example. The dynamic load is a uniform pressure acting on a circular surface with radius a. The load can be modeled by a harmonic function and can be represented in terms of the variables r,t by (1). The selected data for the proposed loading function are Z=1.0N,a=8.2m,ω=0.94rad/s, and the material properties for the soil are c1=355m/s,c2=309m/s, μ=2.137×108N/m2. For the inverse Laplace transforms of (75) and (76), the Stehfest [28] method is adopted. According to this method, the functions ur(r,z,t),uz(r,z,t) can be obtained by the following formula:(77)urr,z,tuzr,z,t≈ln⁡2t∑n=1NcnUrr,z,nln⁡2tUzr,z,nln⁡2twhere(78)cn=-1n+N/2∑k=n+1/2minn,N/2kN/22k!N/2-k!k!k-1!n-k!2k-n!In above equation, N is a free even number representing the number of terms used in the summation. Values within the range 10≤N≤24 yield sufficient accuracy (e.g., [29]). In the present study we adopted the value N=24.Forr=0, the Bessel function J1(ξr) is null. Therefore, (73) yields Ur=0. With the aid of the standard software “Mathematica” [30] the displacement distribution uz versus z for r=0 and t=0.3,0.5,0.7s is demonstrated in Figures 2(a), 2(b), and 2(c).Figure 2
Displacement distributionuz versus z for r=0 and t=0.3,0.5,0.7s.
(a) (b) (c)According to these findings, the displacement tends to reach the z-axis asymptotically due to the termse-az and e-βz in (74).


## 4.2. Verification of the Solution

The derived solution can be adjusted to provide the surface displacement distribution for the specific case of a suddenly applied point load given by the following equation:(79)qr,0,t=Pδr2πrHtwhere δ(r) is the Dirac Delta function of the variable r and Ht is the Heaviside step function of time t. Solution for this specific case is available in [25]:(80)Urr,z,p=P2πμp∫0∞ξ22ξ2+κ2e-αz-2αβe-βzDJ1ξrdξ(81)Uzr,z,p=P2πμp∫0∞-αξ2ξ2+κ2e-αz+2αξ3e-βzDJ0ξrdξFollowing the procedure of [14], the point load P is given by the following equation:(82)P=Zπa2Therefore,(83)Z=Pπa2The above term for the uniform distributed load Z should be placed in (73) and (74). On the other hand, the Laplace transform(84)ωω2+p2=Lsin⁡ωt;t→pshould be replaced in (73) and (74) by the Laplace transform:(85)1p=LHt;t→pFor the case of point load P, the radius a of the circular surface of the equivalent uniform load Z given by (83) is a→0. Therefore, (73) and (74) yield(86)Urr,z,p=lima→0-P/πa2aμ1p∫0∞ξ2ξ2+κ2e-αz-2αβe-βzDJ1ξrJ1aξdξ(87)Uzr,z,p=lima→0P/πa2aμ1p∫0∞-α2ξ2+κ2e-αz+2αβ2-κ2e-βzDJ0ξrJ1aξdξSince(88)lima→0J1aξa=ξ2 (86) and (87) yield the known [25] (80) and (81). With the aid of the Laplace inversion technique proposed in [25], (87) yields the wave propagation uz on the soil surface for the pulse point load of (79):(89)uzr,0,t=0forc2tr<13uzr,0,t=P32μπr6-3c2t/r2-0.25-33+53+3/4-c2t/r2+33-5-3+3/4+c2t/r2for13<c2tr<1uzr,0,t=P32μπr6-33+53+3/4-c2t/r2for1<c2tr<3+32uzr,0,t=3P8μπrforc2tr>3+32For a unit load P=1, the displacement distributions at time instants t=0.0005s, 0.001s, 0.0015s, and 0.002s are demonstrated in Figures 3(a), 3(b), 3(c), and 3(d).Figure 3
3D displacement distribution ofuz due to transient point load for various time instants. Implementation of the proposed solution to transient point load shows full agreement with the existing solution in [25].
(a) (b) (c) (d)


## 5. Conclusions and Discussion

In the present study, an exact solution of the vertical harmonic loading induced dynamic response in the elastic soil has been achieved. The present analytic solution can be utilized for solving the seabed response under vertical harmonic excitation of the wind turbine foundation in waves. The solution is based on integral transforms and stress functions. Verification to the solution has been conducted and a case study has been provided. The achieved solution for harmonic (with respect to time) pressure has been implemented for suddenly applied point load. The derived displacement functions for this transient point loading source are in full agreement with the existing [25] solution. For t=0.0005 s (Figure 3(a)) high values for uz displacement take place in the vicinity of the point r=0 (loading source location). The displacement peak is surrounded by a valley ring which takes place in the area 0.987<r<0.222 m. For t=0.0015 s (Figure 3(b)) the peak is flattened and expanded in the area 0.0<r<0.208 m, and the valley ring is shallow. For t=0.0015 s (Figure 3(c)) the propagating peak is more flat and wide and is expanded in the area 0.0<r<0.388 m. For t=0.025 (Figure 3(d)) the plotted area 0.0<r<0.5 m is covered by the flat peak and no valley (in this area) is visible. The proposed method takes into account the inertial effects and does not require much computation power. Since the solution is analytic, the results are more accurate than numerical ones. The present work can be extended for solving viscoelastic porous media response under time-varying dynamic loads.



---

*Source: 6250158-2018-10-25.xml*