Dynamic or moving loads on the soil surface or inside a tunnel lead to disturbing vibrations in structures, which can reduce the serviceability of structures. Hence, the prediction of wave propagation in soil and structures due to rail or road traffic is necessary. Furthermore, there is a need for a realistic model of the complete Emission-Transmission-Immission-System and therefore a model for the dynamic Soil-Structure-Interaction. This model an then be used for numerical testing of the effectiveness of vibration reduction measures or modeling of moving loads and their effects on soil vibrations.
The general procedure is the following: Emission- and immission-systems have to be modeled appropriately (typically with the Finite Element Method (FEM), as they usually consist of different materials arranged in complex geometries). The soil as transmission medium with infinite extension but rather simple geometry and material behavior can - for the linear case - be described by the Integral Transform Method (ITM). This method provides analytical solutions for a homogenous or layered halfspace with a cylindrical or spherical cavity. At the interaction surface the displacements of both subsystems are coupled using the substructure technique. Thus, the advantages of both methods can be combined to model the total system as accurate as possible.
Fundamentals in Elastodynamics
In linear Elastodynamics, closed analytical solutions are available for the fundamental systems halfspace, fullspace with cylindrical cavity and fullspace with spherical cavity. By superposition of these fundamental systems it is possible to generate a solution for the halfspace with a spherical cavity/indentation or a tunnel structure.
Fundamental Solution Halfspace
Every displacement field inside a linear elastic, homogeneous, isotropic continuum has to fulfill the following Lamé differential equation:
| \mu u^i|_j^j + \left( \lambda + \mu \right) u^j|_j^i - \rho \ddot u^i = 0 |
In order to solve this coupled system of equations for the unknown displacements, a decoupling of the equations is advantageous as a direct integration is difficult. The decoupling is done using Helmholtz’ principle. This principle states that every possible vector field can be replaced by the sum of an irrotational field and a solenoidal field. As the gradient of a scalar field \Phi is irrotational and the rotation of a vector field \Psi is solenoidal, the components of the displacement field u^{i} can be replaced by
| u^i=\Phi |^i+\Psi_l |_k \, \in^{ikl} |
This leads to the following system of decoupled, partial differential equations
| \begin{alignat}{2} & \left[ \frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2} + \frac{\partial^2}{\partial z^2} - \frac{1}{{c_p}^2} \frac{\partial^2}{\partial t^2} \right] \Phi \left( x , y , z , t \right) && = 0\\[1mm] & \left[ \frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2} + \frac{\partial^2}{\partial z^2} - \frac{1}{{c_s}^2} \frac{\partial^2}{\partial t^2} \right] \Psi_x \left( x , y , z , t \right) && = 0\\[1mm] & \left[ \frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2} + \frac{\partial^2}{\partial z^2} - \frac{1}{{c_s}^2} \frac{\partial^2}{\partial t^2} \right] \Psi_y \left( x , y , z , t \right) && = 0 \end{alignat} |
A threefold Fourier Transformation x \circ - \bullet k_x, y \; \circ - \bullet \; k_y, t \; \circ - \bullet \;\omega implying \frac{\partial}{\partial x^2} \,\, \bullet - \circ \,\, (i kx)^2, \frac{\partial}{\partial y^2} \,\, \bullet - \circ \,\, (i ky)^2, \frac{\partial}{\partial t^2} \,\, \bullet - \circ \,\, (i \omega)^2 leads to a system of ordinary DEs
| \begin{alignat}{2} & \left[ - {k_x}^2 - {k_y}^2 + {k_p}^2 + \frac{\partial^2}{\partial z^2} \right] \hat{\Phi} \left(k_x, k_y, z, \omega \right) && = 0 \label{eq:1_34} \\[1mm] & \left[ - {k_x}^2 - {k_y}^2 + {k_s}^2 + \frac{\partial^2}{\partial z^2} \right] \hat{\Psi}_i \left(k_x, k_y, z, \omega \right) && = 0 %\label{eq:1_35} \\[3mm] %& \left[ - {k_x}^2 - {k_y}^2 + {k_s}^2 + \frac{\partial^2}{\partial z^2} \right] \hat{\Psi}_y \left(k_x, k_y, z, \omega \right) && = 0 \label{eq:1_36} \end{alignat} |
with the wavenumbers of the P- and S-waves for a given frequency: k_p=\omega/c_p and k_s=\omega/c_s.
A solution is possible using an exponential approach for the potentials \hat{\Phi} and \hat{\Psi}_i
| \begin{alignat}{2} & \hat{\Phi} && = A_1 \, e^{\lambda_1 z} + A_2 \, e^{- \lambda_1 z} \\[1mm] & \hat{\Psi}_i && = B_{i1} \, e^{\lambda_2 z} + B_{i2} \, e^{- \lambda_2 z} \end{alignat} |
with \lambda_{1}^2=k_x^2+k_{y}^2-k_{p}^2 and \lambda_{2}^2=k_{x}^2+k_{y}^2-k_{s}^2.
According to the Sommerfeld radiation condition for negative frequencies A_1, B_{i1} can be set to zero. The results for positive frequencies are the conjugate complex ones of those calculated for the negative frequencies. A calculation of the unknown coefficients using the local boundary conditions on the halfspace surface is given by
| \begin{pmatrix} \hat{\sigma}_{zz} \\ \hat{\sigma}_{zy} \\ \hat{\sigma}_{zx} \end{pmatrix} = \mu \underbrace{\begin{bmatrix} % 2 {\lambda_1}^2 - \frac{\lambda}{\mu} {k_p}^2 & 2 i k_y \lambda_2 & - 2 i k_x \lambda_2 \\ % - 2 i k_y \lambda_1 & {\lambda_2}^2 + {k_y}^2 & - k_x k_y \\ % - 2 i k_x \lambda_1 & k_x k_y & {\lambda_2}^2 - {k_x}^2 \end{bmatrix}}_{\left[\hat{K}_{hs}\right]} \begin{pmatrix} A_2 \\ B_{x2} \\ B_{y2} \end{pmatrix} = \begin{pmatrix} - \hat{p}_{z} \\ - \hat{p}_{y} \\ - \hat{p}_{x} \end{pmatrix} |
Wave propagation characteristics (exemplarily for \lambda_1. For \lambda_2 analogously.)
Body waves: k_x^2+k_y^2<k_p^2 \;\text{resp.}\; k_s^2 \rightarrow \lambda_1 imaginary
Spatially propagating waves \rightarrow far field in the medium. Depending on sign of \lambda_1, the waves are propagating in the negative or positive z-direction.
Surface waves: k_x^2+k_y^2>k_p^2 \;\text{resp.}\; k_s^2 \rightarrow \lambda_1 real
Near field, that is exponentially increasing or decaying with the depth z in dependency on the sign of \lambda_1
All wavetypes are in the solution of the Lamé DE included
Fundamental Solution Fullspace with Cylindrical Cavity
The cavity is first expressed in cylindrical coordinates x, r, \phi, t (compare figure Fundamental Systems). Using a Fourier Transformation yields x \circ - \bullet k_x, t \; \circ - \bullet \;\omega respectively a Fourier series in direction of the diameter \phi \; \circ - \bullet \;n. This leads to the following system of ordinary DEs after Helmholtz Decomposition:
| \begin{alignat}{2} & \left[ - {k_x}^2 + \frac{\partial^2}{\partial r^2} + \frac{1}{r} \frac{\partial}{\partial r} - \frac{n^2}{r^2} + {k_p}^2 \right] \hat{\Phi} \left( k_x , r , n , \omega \right) && = 0 \label{eq:1_89} \\[3mm] & \left[ - {k_x}^2 + \frac{\partial^2}{\partial r^2} + \frac{1}{r} \frac{\partial}{\partial r} - \frac{n^2}{r^2} + {k_s}^2 \right] \hat{\psi} \left( k_x , r , n , \omega \right) && = 0 \label{eq:1_90} \\[3mm] & \left[ - {k_x}^2 + \frac{\partial^2}{\partial r^2} + \frac{1}{r} \frac{\partial}{\partial r} - \frac{n^2}{r^2} + {k_s}^2 \right] \hat{\chi} \left( k_x , r , n , \omega \right) && = 0 \label{eq:1_91} \end{alignat} |
The solution is obtained with Hankel functions in dependency on the coordinate r and {k_1}^2 = {k_p}^2 - {k_x}^2 respectively {k_2}^2 = {k_s}^2 - {k_x}^2
| \begin{alignat}{2} & \hat{\Phi} \left( k_x , r , n , \omega \right) && = C_{1n} \; H_n^{\left( 1 \right)} \left( k_1 r \right) + C_{4n} \; H_n^{\left( 2 \right)} \left( k_1 r \right) \label{eq:1_95} \\[1mm] & \hat{\psi} \left( k_x , r , n , \omega \right) && = C_{2n} \; H_n^{\left( 1 \right)} \left( k_2 r \right) + C_{5n} \; H_n^{\left( 2 \right)} \left( k_2 r \right) \label{eq:1_96} \\[1mm] & \hat{\chi} \left( k_x , r , n , \omega \right) && = C_{3n} \; H_n^{\left( 1 \right)} \left( k_2 r \right) + C_{6n} \; H_n^{\left( 2 \right)} \left( k_2 r \right) \label{eq:1_97} \end{alignat} |
The unknowns are calculated using the local boundary conditions and the Sommerfeld condition of radiation.
The Hankel functions of first and second kind are calculated out of the Bessel functions of first kind J_n(z) and second kind Y_n(z)
| \begin{align} H_n^1(z)&=J_n(z)+i \, Y_n(z)\\ H_n^2(z)&=J_n(z)-i \, Y_n(z) \end{align} |
ITM Substructure
The full solution is given by a superposition of the stresses at the boundaries of the halfspace surface and the cylinder surface. Applying a stress on the boundary \Lambda induces stress at the boundary \delta \Gamma, which is shown in the left summand in the figure. Vice versa, applying a stress at \Gamma leads to stresses at \delta \Lambda (right summand). For equilibrium, the sum of the stresses at \Lambda and \delta \Lambda must equal the external load \hat{p}_{iz}(s) at this boundary. The same holds for the stresses at \Gamma respectively \delta \Gamma and the load \hat{p}_{jr}(n). Note, that the stresses are expressed in different coordinates at the two boundaries.
| \begin{align} C_{iz}(s) \cdot \hat{\sigma}_{iz}(s)+\sum\limits_{N}\sum\limits_{J} C_{jr}(n)\hat{\sigma}_{iz}^{(jr,n)}(s) &= \hat{p}_{iz}(s)\\[3pt] \sum\limits_{S}\sum\limits_{I} C_{iz}(s) \hat{\sigma}_{jr}^{(iz,s)}(n) + C_{jr}(n)\hat{\sigma}_{jr}(n) &= \hat{p}_{jr}(n) \end{align} |
| \begin{align} \Rightarrow\; \begin{bmatrix} \hat{S}_{ITM} \end{bmatrix} \begin{pmatrix} C \end{pmatrix} = \begin{bmatrix} -\hat{P}_{ITM} \end{bmatrix}\;\Rightarrow\; \begin{pmatrix} C \end{pmatrix} = \begin{bmatrix} \hat{S}_{ITM} \end{bmatrix}^{-1} \begin{bmatrix} -\hat{P}_{ITM} \end{bmatrix} \end{align} |
Computation of the displacements at the halfspace surface and the cylinder surface
| \begin{align} \begin{pmatrix} \hat{u}_{ITM} \end{pmatrix} = \begin{bmatrix} \hat{V}_{ITM} \end{bmatrix} \begin{pmatrix} C \end{pmatrix} \end{align} |
Combining stresses and displacements
| \begin{align} \begin{pmatrix} \hat{u}_{ITM} \end{pmatrix} = \begin{bmatrix} \hat{V}_{ITM} \end{bmatrix} \begin{bmatrix} \hat{S}_{ITM} \end{bmatrix}^{-1} \begin{pmatrix} -\hat{P}_{ITM} \end{pmatrix} \end{align} |
With the flexibility resp. stiffness matrix
| \begin{align} \begin{bmatrix} \hat{N}_{ITM} \end{bmatrix} = \begin{bmatrix} \hat{K}_{ITM} \end{bmatrix}^{-1} = \begin{bmatrix} \hat{V}_{ITM} \end{bmatrix} \begin{bmatrix} \hat{S}_{ITM} \end{bmatrix}^{-1} \end{align} |
| \begin{align} \begin{bmatrix} \left[\hat{K}_{{\Lambda\Lambda}_{ITM}}\right] & \left[\hat{K}_{{\Lambda\Gamma}_{ITM}}\right] \\[2mm] \left[\hat{K}_{{\Gamma\Lambda}_{ITM}}\right] & \left[\hat{K}_{{\Gamma\Gamma}_{ITM}}\right] \end{bmatrix} \begin{pmatrix} \hat{u}_{{\Lambda}_{ITM}} \\[2mm] \hat{u}_{{\Gamma}_{ITM}} \end{pmatrix} = \begin{pmatrix} \hat{P}_{\Lambda_{ITM}} \\[2mm] \hat{P}_{\,\Gamma_{ITM}} \end{pmatrix} \end{align} |
Note, that \hat{u}_{{\Lambda}_{ITM}} and \hat{u}_{{\Gamma}_{ITM}} are expressed in different coordinates.
FEM Substructure
To model complex structures and inhomogeneous soil conditions in the surrounding area the Finite-Element-Method is used. The stiffness matrices for structures with infinite extension in one dimension (2,5D) are calculated using a Fourier Transformed Finite-Element-Method formulation. This means, that the third dimension in x-direction is transformed into x \circ - \bullet k_x. The structure in the yz-plane is then modelled using shell elements. The stiffness matrix is obtained using the principle of virtual work:
| \begin{alignat}{2} &\text{Virtual Work} && \delta W_i = - \frac{1}{2 \pi}\int\limits_{-\infty}^\infty \int\limits_{(A)} \delta \tilde{\pmb{\varepsilon}}^* (k_x,\eta,\zeta) \; \tilde {\pmb{\sigma}} (k_x,\eta,\zeta) \; dA \; dk_x\\[3pt] &\text{Stiffness matrix} && \left[\tilde K \right] = \sum\limits_{k=1}^{n_{GP}} \left[\tilde B (k_x, \eta_k, \zeta_k)\right]^H \; \left[D'\right] \;\left[\tilde B (k_x, \eta_k, \zeta_k)\right] \det (J) \; w_k\\[3pt] &\text{} && \end{alignat} |
With the \left[B \right]-matrix in the twofold Fourier-transformed domain with bilinear form functions
| \begin{equation*} \left[\tilde{B}\right] = \begin{bmatrix} i k_x N_1 \left(\eta, \zeta \right) & 0 & 0 & i k_x N_2 \left(\eta, \zeta \right) & 0 & \cdots \\ 0 & \frac{\partial N_1 \left(\eta, \zeta \right)}{\partial y} & 0 & 0 & \frac{\partial N_2 \left(\eta, \zeta \right)}{\partial y} & \cdots \\ 0 & 0 & \frac{\partial N_1 \left(\eta, \zeta \right)}{\partial z} & 0 & 0 & \cdots \\ \frac{\partial N_1 \left(\eta, \zeta \right)}{\partial y} & i k_x N_1 \left(\eta, \zeta \right) & 0 & \frac{\partial N_2 \left(\eta, \zeta \right)}{\partial y} & i k_x N_2 \left(\eta, \zeta \right) & \cdots \\ 0 & \frac{\partial N_1 \left(\eta, \zeta \right)}{\partial z} & \frac{\partial N_1 \left(\eta, \zeta \right)}{\partial y} & 0 & \frac{\partial N_2 \left(\eta, \zeta \right)}{\partial z} & \cdots \\ \frac{\partial N_1 \left(\eta, \zeta \right)}{\partial z} & 0 & i k_x N_1 \left(\eta, \zeta \right) & \frac{\partial N_2 \left(\eta, \zeta \right)}{\partial z} & 0 & \cdots \end{bmatrix} \end{equation*} |
| \begin{align*} N_1 \left(\eta, \zeta \right) = \frac{1}{4} \left( 1 - \eta \right) \left( 1 - \zeta \right) \\[1mm] N_2 \left(\eta, \zeta \right) = \frac{1}{4} \left( 1 - \eta \right) \left( 1 + \zeta \right) \\[1mm] N_3 \left(\eta, \zeta \right) = \frac{1}{4} \left( 1 + \eta \right) \left( 1 + \zeta \right) \\[1mm] N_4 \left(\eta, \zeta \right) = \frac{1}{4} \left( 1 + \eta \right) \left( 1 - \zeta \right) \end{align*} |
Coupling of substructures
Recall, that the individual solutions are all expressed in different coordinates, such that the coordinate system for each part is most suitable for this application
- surface \Lambda : k_x, k_y, z, \omega
- cylindrical cavity: k_x, r, n, \omega
- FEM substructure: k_x, y, z, t
Hence, a transformation in a mutual reference system is required on the coupling surface.
| \begin{equation} \begin{bmatrix} \left[\hat{K}_{{\Lambda\Lambda}_{ITM}}\right] & \left[\hat{K}_{{\Lambda\Gamma}_{ITM}}\right] \\[2mm] \left[\hat{K}_{{\Gamma\Lambda}_{ITM}}\right] & \left[\hat{K}_{{\Gamma\Gamma}_{ITM}}\right] \end{bmatrix} \begin{pmatrix} \hat{u}_{{\Lambda}_{ITM}} \\[2mm] \hat{u}_{{\Gamma}_{ITM}} \end{pmatrix} = \begin{pmatrix} \hat{P}_{\Lambda_{ITM}} \\[2mm] \hat{P}_{\Gamma_{ITM}} \end{pmatrix} \label{eq:5_30} \end{equation} |
| \begin{equation} \begin{bmatrix} \left[\tilde {K}_{{\Gamma\Gamma}_{FE}}\right] & \left[\tilde {K}_{{\Gamma_s\Omega}_{FE}}\right] \\[2mm] \left[\tilde {K}_{{\Omega\Gamma}_{FE}}\right] & \left[\tilde {K}_{{\Omega\Omega}_{FE}}\right] \end{bmatrix} \begin{pmatrix} \hat {u}_{{\Gamma}_{FE}} \\[2mm] \tilde {u}_{{\Omega}_{FE}} \end{pmatrix} = \begin{pmatrix} \hat {P}_{\Gamma_{FE}} \\[2mm] \tilde {P}_{\Omega_{FE}} \end{pmatrix} \label{eq:5_31} \end{equation} |
| \begin{equation} \begin{bmatrix} \left[\hat{K}_{{\Lambda\Lambda}_{ITM}}\right] & \left[\hat{K}_{{\Lambda\Gamma}_{ITM}}\right] & 0 \\[2mm] \left[\hat{K}_{{\Gamma\Lambda}_{ITM}}\right] & \left[\hat{K}_{{\Gamma\Gamma}_{ITM}}\right] + \frac{1}{ds} \left[T\right]^{-1} \left[\tilde{K}_{{\Gamma\Gamma}_{FE}}\right] \left[T\right] & \frac{1}{ds} \left[T\right]^{-1} \left[\tilde{K}_{{\Gamma\Omega}_{FE}}\right] \\[2mm] 0 & \left[\tilde{K}_{{\Omega\Gamma}_{FE}}\right] \left[T\right] & \left[\tilde{K}_{{\Omega\Omega}_{FE}}\right] \end{bmatrix} \begin{pmatrix} \hat{u}_{{\Lambda}_{ITM}} \\[2mm] \hat{u}_\Gamma \\[2mm] \tilde{u}_{\Omega_{FE}} \end{pmatrix} = \begin{pmatrix} \hat{P}_{\Lambda_{ITM}} \\[2mm] \hat{P}_\Gamma \\[2mm] \tilde{P}_{\Omega_{FE}} \end{pmatrix} \end{equation} |
Moving Load Effects
A harmonic moving load in transformed domain can be expressed by a frequency shift
| \begin{align} p_{mov}(x,y,z,t)=p_0(x-vt,y,z)\cdot f(t) \qquad \circ - \bullet \qquad \hat{p}_{mov}(k_x,k_y,z,\omega)=\hat p_0(k_x,k_y,z)\cdot \hat f(\omega+v k_x) \end{align} |
For a given function
| \begin{align} \hat{f}(k_x,\omega)=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty} f(x,t) e^{-i k_x x} e^{-i\omega t} dx dt \end{align} |
yields, when introducing a moving coordinate system \bar{x}=x-vt and therefore x=\bar{x}+vt
| \begin{align} \hat{f}(k_x,\omega)=\int\limits_{-\infty}^{\infty}\int\limits_{-\infty}^{\infty} f(\bar{x},t) e^{-i k_x (\bar{x}+vt)} e^{-i\omega t} d\bar{x} dt\\ \hat{f}(k_x,\omega)=\int\limits_{-\infty}^{\infty}\int\limits_{-\infty}^{\infty} f(\bar{x},t) e^{-i k_x \bar{x}} e^{-i(\omega +k_x v)t} d\bar{x} dt \end{align} |
The displacements in the threefold Fourier transformed domain are calculated using the fundamental solution \hat{g}(k_x, k_y, z, \omega):
| \begin{align} \hat{u}(k_x, k_y, z, \omega)=\hat{g}(k_x, k_y, z, \omega) \hat{p}_{mov}(k_x,k_y,z,\omega) \end{align} |
Thus, after a threefold inverse transformation and the introduction of the moving coordinate frame \bar{x}=x-vt the displacement results as:
| \begin{align} \rightarrow \;\bar{u}(\bar{x},y,z,t)=\frac{1}{4\pi^2} \operatorname{Re} \left(e^{-i\Omega t} \int\limits^{\infty}_{-\infty} \int\limits^{\infty}_{-\infty} \frac{\hat{p}_0(k_x,k_y)}{\hat{K}(k_x,k_y,z,-\omega-v k_x)} e^{i(\bar{x}k_x+y k_y)} dk_y dk_x \right) \end{align} |
with \hat{K}(k_x,k_y,z,z,\omega)=\frac{1}{\hat{g}(k_x, k_y, z, \omega)}
Numerical Examples
One possible application for the coupled ITM-FEM approach is the numerical prediction of the effectiveness of vibration reduction measures. The method facilitates to introduce mitigation measures at the source, the receiver or in the transmission path. Furthermore, the prediction of the environmental vibrations caused by machines or the reaction of structures due to an excitation in the soil is possible. Also moving load effects in the supersonic regime like the formation of a shock wave (Mach Cone) or the Doppler-Effect can be modeled with this approach.
Vibration Reduction
To reduce the train induced vibrations in tunnels an elastic layers can be inserted between track and tunnel thus creating a mass-spring system. The effectiveness of the mass-spring system can be predicted by a numerical model. The figure on the right shows the insertion Loss R=20 \, \log(\frac{w_{wo,el}}{w_{el}}) in dB due to a harmonic load within the tunnel evaluated at the center of the halfspace surface.
Moving Load Effects
The left image shows a shock wave formation due to supersonic load speeds c_i. On the right image, the vertical displacement on halfspace surface due to an oscillating load with f=64 \,\text{Hz}, moving with v=256 \, \frac{\text{m}}{\text{s}} corr. to M_s = c/c_s=1.058 can be seen.
This figure shows the Vertical displacement on the surface of a halfspace due to a harmonic, moving load in the tunnel for v=32\, \frac{\text{m}}{\text{s}}\rightarrow M_s=0.14 \,(a), v=64\, \frac{\text{m}}{\text{s}}\rightarrow M_s=0.28 \,(b), v=96\, \frac{\text{m}}{\text{s}}\rightarrow M_s=0.42 \,(c), and v=160\, \frac{\text{m}}{\text{s}}\rightarrow M_s=0.71 \,(d) cp.
References
Dinkel, Jens. Ein semi-analytisches Modell zur dynamischen Berechnung des gekoppelten Systems Fahrzeug-Fahrweg-Untergrund für das Oberbausystem Feste Fahrbahn. Diss. Technische Universität München, 2000.
Eringen, A. Cemal, Erdogan S. Suhubi, and C. C. Chao. "Elastodynamics, Vol. II, Linear Theory." (1978)
Freisinger, Julian, and Gerhard Müller. "Modellierung eines Halbraums mit sphärischem oder zylinderförmigem Hohlraum für dreidimensionale Boden-Bauwerk-Interaktion." 6. VDI-Fachtagung Baudynamik 2018. 2018.
Frühe, Georg. Überlagerung von Grundlösungen in der Elastodynamik zur Behandlung der dynamischen Tunnel-Boden-Bauwerk-Interaktion. Diss. Technische Universität München, 2011.
Freisinger, J., M. Hackenberg, and G. Müller. "A coupled Integral Transform Method-Finite Element Method approach to model the Soil Structure Interaction of finite (3D) and length invariant (2.5 D) systems." Journal of Sound and Vibration (2020): 115443.
Hackenberg, M., M. Dengler, and G. Müller. "Implementation of the finite element method in the Fourier-transformed domain and coupling with analytical solutions." Eurodyn. 2014.
Müller, Gerhard. "Ein Verfahren zur Erfassung der Fundament-Boden Wechselwirkung unter Einwirkung periodischer Lasten." (1989).
Studer, Jost A., Martin G. Koller, and Jan Laue. Bodendynamik. Springer Berlin Heidelberg, 2008.
Wagner, L. . 2.5D ITM-FEM-Ansatz zur Modellierung eines Halbraums mit zylinderförmigem Graben mit variablen Einschlüssen und Aufbauten. Master thesis, Technische Universität München, München, 2019.