Dynamics of a rotor-structure-soil system: transient response by iterative coupling

Abstract

In this paper, the transient dynamics of a rotor-foundation-structure-soil system is studied to obtain unbalance response of the rotor considering the influence of the structure and the unbounded soil. Transient responses are obtained through iterative coupling between the rotor subsystem and the frame-soil subsystem. The non-linear rotor subsystem is represented by a Laval rotor with rigid bearings and external and internal damping. These equations of motion are solved with the fourth-order Runge-Kutta method. The frame is modeled by the Finite Element Method, while the homogeneous half-space is modeled by the 3D Direct Boundary Element Method in the frequency domain. To derive time-domain equivalent equations of motion, a methodology is proposed based on extracting modal parameters from Frequency Response Functions of the coupled soil-foundation system, using the Rational Fraction Polynomial Method. The methodology renders transient response for the rotor and structure with small time steps, allowing an accurate simulation of the rotor runup phase and the dynamics of the system going through resonance frequencies.

Keywords:
Dynamics of soil-structure Interaction; rotor dynamics; modal analysis; transient response

Graphical Abstract

1 INTRODUCTION

Rotors play a significant role in the generation or transformation of energy, as exemplified by gas turbines, airplane jet engines and more recently as part of wind turbines. Understanding the dynamics of turbines has been very important for the development of energy systems. This is still valid today with the increasing importance of wind turbines (onshore or offshore) on the generation of a more sustainable energy matrix.

Systems composed of rotors, structures, foundations and the soil are complex and difficult to model numerically, mainly due to radiation or geometric damping presented by the soil dynamics. Kellezi and Hansen (2002) conducted a study on the dynamics of windmills in the time domain, employing a three-dimensional axisymmetric transient finite element model (FEM) to analyze an offshore windmill on mono-pile foundation. Their approach included radiation damping via absorbing boundary conditions and considered the nonlinear behavior of the soil to determine deformations, velocities, and accelerations at the pile's top. Liu et al. (2019) developed a mixed frequency/time domain methodology to estimate the transient response of offshore wind turbines (OWT), utilizing initial conditions other than zero by discretizing external loads by their eigenvalues and corresponding complex coefficients. A recent review by Jahani et al. (2022) reveals the enormous amount of research that has been undertaken to analyze the dynamics of offshore wind turbines (OWT). The article by Filho et al. (2021) presents a transient dynamic analysis of two eolic systems composed of structure, rotor and blades under seismic and wind excitation. In-plane and out-of-plane vibrations are considered. The influence of the soil has also been incorporated into the analysis by modelling the soil behavior as a rotational spring. The results presented by the authors point to significant influence of the soil response on the dynamics of the coupled system.

Dynamic coupled-field problems in mechanics have traditionally been solved by partitioning the governing equations into subsystems, which are handled by subsystem analyzers. Subsystems are selected based on considerations such as weak-coupling or differing time response characteristics. Various methods exist for iterative coupling between subsystems, as described in Felippa and Park (1980). The partitioned method performs temporal integration separately on each component of the system, using either sequential or parallel execution of single-field analyzers. A specific type of partitioned solution, the staggered coupling solution, which is the focus of the work of Felippa and Park (1980), organizes the procedure through sequential execution of single-field analyzers. This latter approach offers two potentially important advantages, being an improvement in program modularity and greater computational efficiency.

The Boundary Element Method (BEM) has been extensively used to describe dynamic soil-structure interaction (DSSI) problems, especially in the frequency domain. The BEM can consider in a natural way the geometric damping that results from outgoing waves generated at the soil-foundation interface and that are not reflected and, consequently, withdraw energy from the system (Manolis and Beskos, 1988; Dominguez, 1992, Carrion at al., 2007). A methodology to describe the transient behavior of soil-foundation-structures systems was presented by Ferraz et al. (2023). The present article extends the work of Ferraz et al. (2023) to include the transient analysis of a soil-foundation-structure-rotor systems. The interaction between rotor and soil has previously been investigated in the frequency domain by Gasch et al. (1984), who demonstrated that soil–foundation coupling can substantially reduce rotor vibration amplitudes and improve system stability.

In this paper, the rotor is a Laval rotor resting of a base supported by a frame structure. The frame structure interacts with the soil profile. The soil is treated as a 3D visco-elastic half space. The complete system is subdivided into two subsystems. The first subsystem is the structural frame and the soil. The solution of this coupled soil-structure system is obtained in the frequency domain. The dynamics of this subsystem is characterized by a series of Frequency Response Functions (FRFs). From these FRFs, the modal parameters of an equivalent time-domain dynamic system, containing the dynamics of the couple soil-structure system, is obtained. This set of orthogonal time-domain equations may be integrated by standard numerical methods to render the transient behavior of the soil-structure system. The second sub-system is the rotor and its rigid base. The rotor presents elastic stiffness as well as internal and external damping mechanism and is excited by an external torque. The equations of motion of this sub-system are non-linear and their integration will render the transient response of the rotor and the rigid base. The coupling of the two sub-systems follows the staggered coupling procedure as presented by Felippa and Park (1980).

This strategy is applied to obtain the transient behavior of the coupled soil-structure-foundation-rotor system and used to analyze the rotor response in the runup phase. During the runup phase of the rotor, the system passes through various resonances, in which all displacement amplitudes are amplified and, consequently, the damping mechanisms dissipate an increased amount of energy. If the energy being dissipated at a given resonance is larger than the energy the input torque is adding to the system, the rotor will not be able to pass through the resonances creating potential damage to the system. The analysis performed in this article results in the transient responses of this complex system with very short time steps. To highlight the influence of the soil on the rotor response, three models are considered. The first model is a rotor on rigid bearings. The second model consists of a rotor with a base supported by a frame structure. The structure is supported by a rigid soil. In the third model the frame structure is supported by a flexible soil. The analysis performed in this article helps to understand the influence of the frame structure stiffness and damping mechanisms on the unbalance rotor response. The inclusion of the soil helps to understand the role of geometric damping on the rotor response. Emphasis is placed on the torque requirements for the rotor to overcome the distinct resonances during the runup phase.

3 METHODOLOGY

Figure 1 presents the rotor-foundation-frame-soil system analyzed in the proposed methodology. The system comprises a Laval rotor mounted on a rigid foundation, which is in turn supported by a frame resting on a soil medium. For modeling purposes, the system is partitioned at the interface between the rotor foundation and the frame. The equations of motion for each subsystem are formulated and solved directly in the time domain. Coupling between the subsystems is achieved through an iterative coupling approach. System excitation is introduced via an external torque applied to the rotor.

Figure 1
Coupling of rotor-foundation and soil-structure systems.

The rotor subsystem is characterized by three degrees of freedom, namely the vertical (u_rz) and horizontal (u_ry) lateral displacements of the rotor shaft, and the angular position of the rotor’s center of mass (D) relative to its geometric center (θ_D): {u_r}={{u_rz, u_ry, θ_D}^T. Similarly, the rigid foundation supporting the rotor exhibits three degrees of freedom, denoted as {u_f}={u_fz, u_fy, θ_f}^T. The frame structure presents three degrees of freedom for the i-th node {u_ti}={u_tzi, u_tyi, θ_txi}^T.

3.1 Subsystem I: rotor-rigid foundation

Figure 2 illustrates the rotor subsystem, modeled as a Laval (Jeffcott) rotor with rigid bearings. The system includes a disk with an unbalanced mass m_D positioned at the center of the elastic shaft (c_gr), with both internal c_Iy, c_Iz and external damping c_Ey, c_Ez coefficients, respectively, in the y and z directions, considered. The displacement of the rotor’s center of mass (c_mr)=(D), denoted as u_rz and u_ry, is expressed as the superposition of the displacement of the rotor’s geometric center (c_gr), represented by u_z and u_y, and the relative displacement between the (c_gr) and (c_mr) denoted by v_z and v_y, or u_rz= u_z+ v_z and u_ry= u_y+ v_y. The elastic stiffness of the rotor shaft in the y and z directions is represented by k_ry and k_rz. The rotor is excited by an external torque T.

Figure 2
Scheme of Laval rotor.

The kinematic scheme of the rotor-foundation subsystem is depicted in Figure 3. The model accounts for lateral vertical u_fz and horizontal u_fy displacements, as well as rotational motion θ_f of the foundation. The inertial coordinate system is defined at point (F), which corresponds to the foundation center of mass (c_m)=(F). The parameter h_o denotes the vertical distance between the bearing location (O) and the point (F), while h_F represents the distance from the bottom of the foundation (Q) to its mass center (F). The eccentricity of the rotor, denoted by e, is defined as the distance between the geometric center of the rotor (c_gr) its mass center designated as point (c_mr)=(D).

Figure 3
Rotor and foundation kinematics.

To derive the equations of motion for the rotor and foundation, it is necessary to define the forces acting on each component, as illustrated in the diagram of Figure 4. Following the approach of Mesquita et al. (2006), the governing equations are defined. For the rotor, the forces acting at the geometric center of the disk (c_gr) include the elastic force F_Er, the external damping force F_Ar, the internal damping force F_Ir, and the weight P_D, which acts at the rotor’s center of mass (D)=(c_mr). The rotor weight causes a static displacement u_est. The reaction forces from the rotor, {R_R}={F_Er+ F_Ar + F_Ir + P_D }^T are transmitted to the rotor bearing at point O′. The rigid foundation mass is m_f and the weight (P_F) acts at mass center (F)=(c_m) The forces, acting at the bottom of the rotor rigid base, point (Q), that is, at the interface between the rigid foundation base and the frame structure, are: {F_q}={F_qz, F_qy ,M_qx}^T.

Figure 4
Scheme of rotor subsystem.

Considering g as the acceleration due to gravity and the unit vectors in the y and z directions represented as n_y and n_z, the forces acting on the rotor {R_R} can be expressed as:

$$F_{Er}=(-k_{ry}{u}_{ry}+k_{ry}{u}_{fy}+k_{ry}{h}_O{\theta }_f)n_y+(-k_{rz}{u}_{rz}+k_{rz}{u}_{fz}-k_{rz}{u}_{est})n_z$$ (1)

$$F_{Ar}=(-c_{Ey}{\dot{u}}_{ry}+c_{Ey}{\dot{u}}_{fy}+c_{Ey}{h}_O{\dot{\theta }}_f)n_y+(-c_{Ez}{\dot{u}}_{rz}+c_{Ez}{\dot{u}}_{fz})n_z$$ (2)

$$\begin{array}{l}{F}_{Ir}=(-c_{Iy}{\dot{u}}_{ry}+c_{Iy}{\dot{u}}_{fy}+c_{Iy}{h}_O{\dot{\theta }}_f-c_{Iy}{u}_{rz}{\dot{\theta }}_D+c_{Iy}{u}_{fz}{\dot{\theta }}_D-c_{Iy}{u}_{est}{\dot{\theta }}_D)n_y+\\ +(-c_{Iz}{\dot{u}}_{rz}+c_{Iz}{\dot{u}}_{fz}+c_{Iz}{u}_{ry}{\dot{\theta }}_D-c_{Iz}{u}_{fy}{\dot{\theta }}_D-c_{Iz}{h}_O{\theta }_f{\dot{\theta }}_D)n_z\end{array}$$ (3)

$$P_D=m_{D}g{n}_z$$ (4)

The forces acting on the foundation are the weight P_F, the reaction forces or the rotor {R_R}={F_Er+F_Ar+F_Ir+P_D}^T, acting on point (O) and the forces acting at point (Q) at the bottom of the foundation: {F_q}={F_qz, F_qy, M_qx}^T:

$$P_F=m_{F}g{n}_z$$ (5)

$$R_R=F_{Er}+F_{Ar}+F_{Ir}+P_D$$ (6)

The moments produced by the equations (1) through (4) with respect to the rotor center of mass (D) will lead to the equations of motion (7) through (9):

$$M^{\raisebox{1ex}{$F_{Er}$}\!\left/ \!\raisebox{-1ex}{$D$}\right.}=\left[\begin{array}{l}-k_{ry}esen{\theta }_D{u}_{ry}+k_{ry}esen{\theta }_D{u}_{fy}+k_{ry}{h}_{O}esen{\theta }_D{\theta }_f+\\ +k_{rz}e\cos {\theta }_D{u}_{rz}-k_{rz}e\cos {\theta }_D{u}_{fz}+k_{rz}e\cos {\theta }_D{u}_{est}\end{array}\right]n_x$$ (7)

$$M^{\raisebox{1ex}{$F_{Ar}$}\!\left/ \!\raisebox{-1ex}{$D$}\right.}=\left[-c_{Ey}esen{\theta }_D{\dot{u}}_{ry}+c_{Ey}esen{\theta }_D{\dot{u}}_{fy}+c_{Ey}{h}_{O}esen{\theta }_D{\dot{\theta }}_f+c_{Ez}e\cos {\theta }_D{\dot{u}}_{rz}-c_{Ez}e\cos {\theta }_D{\dot{u}}_{fz}\right]n_x$$ (8)

$$M^{\raisebox{1ex}{$F_{Ir}$}\!\left/ \!\raisebox{-1ex}{$D$}\right.}=\left[\begin{array}{l}-c_{Iy}esen{\theta }_D{\dot{u}}_{ry}+c_{Iy}esen{\theta }_D{\dot{u}}_{fy}+c_{Iy}{h}_{O}esen{\theta }_D{\dot{\theta }}_f-c_{Iy}{u}_{rz}esen{\theta }_D{\dot{\theta }}_D+\\ +c_{Iy}{u}_{fz}esen{\theta }_D{\dot{\theta }}_D-c_{Iy}{u}_{est}esen{\theta }_D{\dot{\theta }}_D+c_{Iz}e\cos {\theta }_D{\dot{u}}_{rz}-c_{Iz}e\cos {\theta }_D{\dot{u}}_{fz}-\\ +c_{Iz}{u}_{ry}e\cos {\theta }_D{\dot{\theta }}_D+c_{Iz}{u}_{fy}e\cos {\theta }_D{\dot{\theta }}_D+c_{Iz}{h}_O{\theta }_{f}e\cos {\theta }_D{\dot{\theta }}_D\end{array}\right]n_x$$ (9)

with

$$T_R=T{n}_x$$ (10)

The moments produced by the equations (5) and (6) and the vector {F_q} with respect to the foundation center of mass (F) will lead to the equations of motion (11) and (12):

$$M^{\raisebox{1ex}{$R_R$}\!\left/ \!\raisebox{-1ex}{$F\text{'}$}\right.}=\left[\begin{array}{l}(c_{Iy}{h}_O+c_{Ey}{h}_O){\dot{u}}_{ry}-(c_{Iy}{h}_O+c_{Ey}{h}_O){\dot{u}}_{fy}-(c_{Iy}{h}_O^2+c_{Ey}{h}_O^2){\dot{\theta }}_f+\\ +(c_{Iy}{\dot{\theta }}_D)h_O{u}_{rz}+k_{ry}{h}_O{u}_{ry}-(c_{Iy}{\dot{\theta }}_D)h_O{u}_{fz}-k_{ry}{h}_O{u}_{fy}-(k_{ry}{h}_O^2){\theta }_f\end{array}\right]n_x$$ (11)

$$M^{\raisebox{1ex}{$F_q$}\!\left/ \!\raisebox{-1ex}{$F\text{'}$}\right.}=(-F_{qz}{h}_F{\theta }_f-F_{qy}{h}_F-M_{qx})n_x$$ (12)

Denoting the acceleration of the rotor center of mass (D) with respect to the inertial frame I as ^Ia^D’/O and/or the foundation center of mass (F) by ^Ia^F’/F the conservation of linear momentum will lead to:

$$m_C{}^{I}a{}^{\raisebox{1ex}{$C$}\!\left/ \!\raisebox{-1ex}{$C\text{'}$}\right.}={\displaystyle \sum F_c}$$ (13)

Representing the resultant of the moments applied to the rotor center of mass (D) and the foundation center of mass (F) as M^C/C’’ and the quantity of angular momentum of the rotor or foundation as ^IH^C/C’, the conservation of angular momentum will lead to:

$${}^I\dot{H}{}^{\raisebox{1ex}{$C$}\!\left/ \!\raisebox{-1ex}{$C\text{'}$}\right.}={\displaystyle \sum M^{\raisebox{1ex}{$C$}\!\left/ \!\raisebox{-1ex}{$C\text{'}$}\right.}}$$ (14)

Where the for the rotor the indexes are given by c=D and c’=O. Analogously for the foundation c=F and c’=F’. Where (^IḢ^C/C’) represents the time derivative of the angular momentum of the rotor (C=D) or foundation (C=F). By substituting equations (1) to (4) and (7) to (10) into equations (13) and (14), respectively, the equations of motion for the rotor are obtained. Similarly, substituting equations (5) and (6), along with (11) and (12), into equations (13) and (14), respectively, yields the equations of motion for the foundation.

Since the coupling between the foundation and the frame occurs at the base of the foundation, point (Q), it is necessary to transform the degrees of freedom originally defined at the foundation's geometric center (F), designated by {u_f}={u_fy, u_fz, θ_f}^T to the degrees of freedom at the base (Q), {u_q}={u_qz, u_qY, θ_q}^T. This transformation is carried out using the relationships provided in equations (15).

$$\left\{\begin{array}{l}{u}_{fz}=u_{qz}\hfill \\ u_{fy}=u_{qy}+h_F{q}_q\hfill \\ q_f=q_q\hfill \end{array}\right.$$ (15)

The resulting equations relating the complete vector of rotor and foundation displacements {u}={u_ry, u_rz, θ_D, u_qy, u_qz, θ_q}^T to the torque excitation T applied at the rotor and the excitation vector {F_q}={F_qz, F_qy, M_qx}^T acting at the foundation base point (Q) are given in equation (16):

$$\left[M_I\right]\{\ddot{u}(t)\}+\left[A_I(u)\right]\{\dot{u}(t)\}+\left[B_I(u,\dot{u})\right]\{u(t)\}=\{F_I(t,u,\dot{u},\ddot{u})\}$$ (16)

For the case in which the shaft elastic stiffness and the external and internal damping are isotropic or, k_r=k_ry=k_rz, c_E=c_Ey=c_Ez and c_I=c_Iy=c_Iz, the detailed expressions for equation (16) are given in the Appendix. These time-domain equations are non-linear and describe the dynamics of the rotor-foundation subsystem I. These are the equations that will be coupled to the soil-structure equations of subsystem II by the staggered coupling procedure.

3.2 Soil-structure subsystem

Figures 5 show the subsystem II. It is composed of a frame structure, Figure 5b (index t) attached to a soil (index s) through a rigid and massless foundation shown in Figure 5c. The structure in this subsystem is modeled using the finite element method (FEM), discretized as frame elements. Each node of the frame presents three degrees of freedom, {u_ti}={u_tzi, u_tyi, θ_txi}^T. The bottom node of the frame structure {u_t1} is attached to the displacement degrees of freedom of the rigid and massless foundation resting on the surface of the half-space, {u_s}={u_sz, u_sy, θ_sx}^T. The displacement of the massless rigid foundation at the soil-frame interface is obtained via the Direct version of the Boundary Element Method (DBEM) in the frequency domain (Carrion et al, 2007). The vector {F_tn}={F_tzn, F_tyn, M_txn}^T describes the excitation acting on the top of the frame structure (node n, point Q). The vector of forces acting at the bottom and first node of the structure is {F_t1}={F_tz1, F_ty1, M_tx1}^T. The vector {F_s}={F_sz, F_sy, M_sx}^T represents the forces acting on the degrees of freedom of the rigid and massless foundation, {u_s}={u_sz, u_sy, θ_sx}^T, (point S). The coupling between the structure and the soil is established by prescribing cinematic continuity {u_t1}={u_s} and equilibrium conditions {F_t1}+{F_s}=0 at the interface, point S (Pacheco et al., 2024).

Figure 5
Scheme of soil-structure subsystem.

The frame structure. The mass and stiffness matrices associated with the frame elements are defined in Equation (17), where, ρ denotes the density of the structure material, E is Young’s modulus, A is the cross-sectional area of the frame, L is the length of the frame element, I_xx represents its moment of inertia and α is the rotational inertia coefficient (Cook et al., 2005).

$$\left[m_{el}\right]=\rho AL\left[\begin{array}{cc}\begin{array}{ccc}0.5& 0& 0\\ 0& 0.5& 0\\ 0& 0& \alpha L^2\end{array}& \begin{array}{ccc}0\text{}& 0& 0\\ 0\text{}& 0& 0\\ 0\text{}& 0& 0\end{array}\\ \begin{array}{ccc}0\text{}& 0& 0\\ 0\text{}& 0& 0\\ 0\text{}& 0& 0\end{array}& \begin{array}{ccc}0.5& 0& 0\\ 0& 0.5& 0\\ 0& 0& \alpha L^2\end{array}\end{array}\right];\left[k_{el}\right]=\left[\begin{array}{cc}\begin{array}{ccc}\raisebox{1ex}{$EA$}\!\left/ \!\raisebox{-1ex}{$L$}\right.& 0& 0\\ 0& \raisebox{1ex}{$12E{I}_{xx}$}\!\left/ \!\raisebox{-1ex}{$L^3$}\right.& \raisebox{1ex}{$6E{I}_{xx}$}\!\left/ \!\raisebox{-1ex}{$L^2$}\right.\\ 0& \raisebox{1ex}{$6E{I}_{xx}$}\!\left/ \!\raisebox{-1ex}{$L^2$}\right.& \raisebox{1ex}{$4E{I}_{xx}$}\!\left/ \!\raisebox{-1ex}{$L$}\right.\end{array}& \begin{array}{ccc}\raisebox{1ex}{$-EA$}\!\left/ \!\raisebox{-1ex}{$L$}\right.& 0& 0\\ 0& \raisebox{1ex}{$-12E{I}_{xx}$}\!\left/ \!\raisebox{-1ex}{$L^3$}\right.& \raisebox{1ex}{$6E{I}_{xx}$}\!\left/ \!\raisebox{-1ex}{$L^2$}\right.\\ 0& \raisebox{1ex}{$-6E{I}_{xx}$}\!\left/ \!\raisebox{-1ex}{$L^2$}\right.& \raisebox{1ex}{$2E{I}_{xx}$}\!\left/ \!\raisebox{-1ex}{$L$}\right.\end{array}\\ \begin{array}{ccc}\raisebox{1ex}{$-EA$}\!\left/ \!\raisebox{-1ex}{$L$}\right.& 0& 0\\ 0& \raisebox{1ex}{$-12E{I}_{xx}$}\!\left/ \!\raisebox{-1ex}{$L^3$}\right.& \raisebox{1ex}{$-6E{I}_{xx}$}\!\left/ \!\raisebox{-1ex}{$L^2$}\right.\\ 0& \raisebox{1ex}{$6E{I}_{xx}$}\!\left/ \!\raisebox{-1ex}{$L^2$}\right.& \raisebox{1ex}{$2E{I}_{xx}$}\!\left/ \!\raisebox{-1ex}{$L$}\right.\end{array}& \begin{array}{ccc}\raisebox{1ex}{$EA$}\!\left/ \!\raisebox{-1ex}{$L$}\right.& 0& 0\\ 0& \raisebox{1ex}{$12E{I}_{xx}$}\!\left/ \!\raisebox{-1ex}{$L^3$}\right.& \raisebox{1ex}{$-6E{I}_{xx}$}\!\left/ \!\raisebox{-1ex}{$L^2$}\right.\\ 0& \raisebox{1ex}{$-6E{I}_{xx}$}\!\left/ \!\raisebox{-1ex}{$L^2$}\right.& \raisebox{1ex}{$4E{I}_{xx}$}\!\left/ \!\raisebox{-1ex}{$L$}\right.\end{array}\end{array}\right]$$ (17)

To define the damping matrix [C_t], Rayleigh damping is adopted in this study, which assumes damping proportional to both the mass and stiffness matrices of the system, as defined in equation (18). The coefficients μ and β correspond to the mass and stiffness proportional damping factors, respectively (Caughey, 1960; Chopra, 2012).

$$\left[C_t\right]=\mu \left[M_t\right]+\beta \left[K_t\right]$$ (18)

The complete equations of motion for the frame structure are:

$${\left[M_t\right]}_{3n\times 3n}{\{{\ddot{u}}_t(t)\}}_{3n\times 1}+{\left[C_t\right]}_{3n\times 3n}{\{{\dot{u}}_t(t)\}}_{3n\times 1}+{\left[K_t\right]}_{3n\times 3n}{\{u_t(t)\}}_{3n\times 1}={\{F_t(t)\}}_{3n\times 1}$$ (19)

The soil dynamic response. The soil dynamic response at the surface of the rigid and massless foundation is described by equation (20), where a denotes half the width of the foundation, G is the shear modulus of the soil and [N_S(ω)] denotes the frequency-dependent matrix containing the flexibility functions of the rigid and massless surface foundation, relating the vector of the displacement at point S, {U_s}={U_sz, U_sy, θ_sx}^T and the forces action on the foundation {F_s}={F_sz, F_sy, M_sx}^T:

$$\underset{{\{U_S(\omega )\}}_{3\times 1}}{\underbrace{\left\{\begin{array}{c}{U}_{sz}(\omega )\\ U_{sy}(\omega )\\ {\theta }_{sx}(\omega )\end{array}\right\}}}=\underset{{\left[N_S(\omega )\right]}_{3\times 3}}{\underbrace{\frac{1}{Ga}\left[\begin{array}{ccc}{N}_{Szz}(\omega )& 0& 0\\ 0& N_{Syy}(\omega )& N_{Syx}(\omega )\\ 0& N_{Sxy}(\omega )& N_{Sxx}(\omega )\end{array}\right]}}\underset{{\{F_S(\omega )\}}_{3\times 1}}{\underbrace{\left\{\begin{array}{c}{F}_{sz}(\omega )\\ F_{sy}(\omega )\\ M_{sx}(\omega )\end{array}\right\}}}$$ (20)

In this study, the cross-coupling terms N_Sxy(ω) and N_Syx(ω) are assumed to be zero. Figure 6 illustrates the frequency-dependent behavior of the main diagonal components of the flexibility matrix, [N_S(ω)].

Figure 6
Flexibilities of rigid massless foundation on the surface of the half space: (a) vertical; (b) horizontal; (c) rocking

Coupling of frame structure and soil dynamic response. To enable the coupling of equations (19) and (20), the former one must be transformed into the frequency domain as shown in equation (21) (Cheng, 1972):

$$\underset{{\left[H_t(\omega )\right]}_{3n\times 3n}}{\underbrace{(-{\omega }^2{[M_t]}_{3n\times 3n}+i\omega {[C_t]}_{3n\times 3n}+{[K_t]}_{3n\times 3n})}}{\{U_t(\omega )\}}_{3n\times 1}={\{F_t(\omega )\}}_{3n\times 1}$$ (21)

In equation (21) the frequency dependent variables are written in capital letters. This equation can be partitioned into submatrices by isolating the first node {U_t1}={U_tz1i, U_ty1, θ_tx1i}^T, which will be coupled with the rigid foundation node S, {U_s}={U_sz, U_sy, θ_sx}:

$$\left[\begin{array}{cc}{\left[H_{t11}(\omega )\right]}_{3\times 3}& {\left[H_{t1m}(\omega )\right]}_{3\times 3m}\\ {\left[H_{tm1}(\omega )\right]}_{3m\times 3}& {\left[H{t}_{mm}(\omega )\right]}_{3m\times 3m}\end{array}\right]\left\{\begin{array}{c}{\{U_{t1}(\omega )\}}_{3\times 1}\\ {\{U_{tm}(\omega )\}}_{3m\times 1}\end{array}\right\}=\left\{\begin{array}{c}{\{F_{t1}(\omega )\}}_{3\times 1}\\ {\{F_{tm}(\omega )\}}_{3m\times 1}\end{array}\right\}$$ (22)

In equation (22) the index m=n-1. By imposing equilibrium at the interface node S, {F_t1}+{F_s}=0 and rearranging equations (20) and (22) will lead to:

$$\left[\begin{array}{cc}{\left[H_{11}(\omega )\right]}_{3\times 3}& {\left[H_{1m}(\omega )\right]}_{3\times 3m}\\ {\left[H_{m1}(\omega )\right]}_{3m\times 3}& {\left[H_{mm}(\omega )\right]}_{3m\times 3m}\end{array}\right]\left\{\begin{array}{c}{\{U_{t1}(\omega )\}}_{3\times 1}\\ {\{U_{tm}^{*}(\omega )\}}_{3m\times 1}\end{array}\right\}=\left\{\begin{array}{c}{\left[N_S(\omega )\right]}_{3\times 3}^{-1}{\{U_S(\omega )\}}_{3\times 1}\\ {\{F_{tm}(\omega )\}}_{3m\times 1}\end{array}\right\}$$ (23)

By reformulating the first expression of Equation (23), the degrees of freedom of the frame structure that are not in contact with the soil {U_tm} can be expressed as a function of the degrees of freedom at the interface {U_t1}, yielding the following relationship:

$${\{U_{tm}^{*}(\omega )\}}_{3m\times 1}=\underset{{\left[H_{\mathrm{mod}}(\omega )\right]}_{3m\times 3}}{\underbrace{{[H_{1m}(\omega )]}_{3\times 3m}^{-1}({[N_S(\omega )]}_{3\times 3}^{-1}-{[H_{11}(\omega )]}_{3\times 3})}}{\{U_{t1}(\omega )\}}_{3\times 1}$$ (24)

Equation (24) relates the displacement vector at the first structural node {U_t1}={U_s} to the displacements of the remaining structural notes (m=n-1). But now the influence of the soil is present through the dynamic soil flexibility matrix [N_s(ω)]. The original vector {U_tm} has been renamed {U^*_tm} to account for this new influence of the soil. The link between the displacements {U_t1} and {U^*_tm} is given by the matrix [H_mod(ω)].

By substituting Equation (24) into the second expression of equation (23), the response of the structural degrees of freedom {U^*_tm} can be expressed as a function of the force vector {F_tm}:

$${\{U_{tm}^{*}(\omega )\}}_{3m\times 1}=\underset{{\left[S_{str}(\omega )\right]}_{3m\times 3m}}{\underbrace{{({[H_{m1}(\omega )]}_{3m\times 3}{[H_{\mathrm{mod}}(\omega )]}_{3m\times 3}^{-1}+{[H_{mm}(\omega )]}_{3m\times 3m})}^{-1}}}{\{F_{tm}(\omega )\}}_{3m\times 1}$$ (25)

In equation (25) the index m stands for all the structure degrees of freedom that are not in contact with the soil. Equation (25) relates the original loading vector acting on the m structural DOFs {F_tm} to the vector of the m structural displacements {U^*_tm}, modified by the inclusion of the soil response. By substituting equation (25) into equation (24), the response at Node 1, {U_t1}={U_s}, can be obtained as a function of the load vector {F_tm}:

$${\{U_{t1}^{*}(\omega )\}}_{3\times 1}=\underset{{\left[S_{soil}(\omega )\right]}_{3\times 3m}}{\underbrace{{\left[H_{\mathrm{mod}}(\omega )\right]}_{3m\times 3}^{-1}{\left[S_{str}(\omega )\right]}_{3m\times 3m}}}{\{F_{tm}(\omega )\}}_{3m\times 1}$$ (26)

The matrices [S_str(ω)] and [S_soil(ω)], presented in equations (25) and (26), respectively, contain the components representing the structure's frequency response functions (FRFs) under the influence of soil-structure interaction. If the structure is excited by the load vector {F_tm}, the frequency response of all structural DOFs {u_t1}+{u_tm} can be calculated. Figure 7 illustrates the typical results obtained in such an analysis. The shown 5-DOF frame is analyzed. In the first case the structure is fixed to a rigid ground. In the second case the structure is supported by a half-space. Typical FRFs for the fixed base and the soil supported structure are shown in Figures 7. The main difference induced by the soil in this example is an increase in the damping systems.

Figure 7
FRFs for structure on fixed base and supported by the soil.

Once the modified FRFs for the structures interacting with the soil are obtained, using the Rational Fraction Polynomial Method (RFPM) (Ewins, 2000), modal parameters can be extracted to create a modal basis that accounts for the influence of the soil. From the extracted set of modal parameters, it is possible to define an equivalent set of orthogonal equations of motion in modal coordinates in the time domain (Ferraz et al., 2023):

$$\left[I\right]\{{\ddot{q}}_t(t)\}+\left[2{\xi }_n{\omega }_n\right]\{{\dot{q}}_t(t)\}+\left[{\omega }_n^2\right]\{q_t(t)\}={\left[\Phi \right]}^T\{F_{tn}^{}(t)\}$$ (27)

The response in physical coordinates can the retrieved by transforming the nodal coordinates {q} to the physical ones {u^*_t} by the know relation {u^*_t}=[Φ]{q_t}. The resulting equations of motion in time domain, represent the solution of the subsystem II (Figure 5) and can be integrated yielding the response of the soil-structure system.

3.3 Staggered coupling solution

Once the equations of motion for the two subsystems are defined, the next step is to couple their transient responses. Following the approach proposed by Felippa and Park (1980), and considering that the coupling occurs between structural subsystems, the staggered coupling solution is adopted. This method enhances computational efficiency by eliminating the need for iterative procedures within each time step to achieve convergence.

Figure 8 presents a step-by-step flowchart of the coupling procedure. The process begins with the definition of a predictor term—specifically, the interface forces and bending moment between the two subsystems (F_qz​, F_qy​ and M_qx​​) as defined in equation (16)). At the initial time step, the system's initial conditions of displacement and velocity are used to calculate the predictor term. For subsequent steps, the predictor term is calculated using data from the preceding time step (n). Once the predictor term is computed, it serves as the excitation input for the equations governing the soil-structure subsystem (Subsystem II), as defined in equation (27), enabling the calculation of the linear and angular displacements of the frame at time step n+1. These displacements at the rotor interface are then used to integrate the equations of motion for the rotor system (Subsystem I), also at time step n+1, as described by equation (16). The procedure then advances to the next time step by calculating the new predictor term.

Figure 8
Flowchart of coupling procedure.

The predictor term is derived from the force balance in Subsystem I, as formulated in Equation (16) and detailed in equations (28) to (30).

$$\begin{array}{l}{F}_{qz}{(}^t=-m_F{\ddot{u}}_{sz}+(c_E+c_I){\dot{u}}_{rz}-(c_E+c_I){\dot{u}}_{qz}+k_r{u}_{rz}-(c_I{\dot{\theta }}_D)u_{ry}-k_r{u}_{qz}+\\ +(c_I{\dot{\theta }}_D)u_{qy}+c_I{\dot{\theta }}_D(h_O+h_F){\theta }_q\end{array}$$ (28)

$$\begin{array}{l}{F}_{qy}{(}^t=-m_F({\ddot{u}}_{qy}+h_F{\ddot{\theta }}_q)+(c_E+c_I){\dot{u}}_{ry}-(c_E+c_I){\dot{u}}_{qy}-(c_E+c_I)(h_O+h_F){\dot{\theta }}_q+\\ +(c_I{\dot{\theta }}_D)u_{rz}+k_r{u}_{ry}-(c_I{\dot{\theta }}_D)u_{qz}-k_r{u}_{qy}-k_r(h_O+h_F){\theta }_q\end{array}$$ (29)

$$M_{qx}{(t)}^{(n)}=-I_{xx}^{\raisebox{1ex}{$R$}\!\left/ \!\raisebox{-1ex}{$F$}\right.}{\ddot{q}}_q-m_F{h}_F{\ddot{u}}_{qy}-m_F{h}_F^2{\ddot{q}}_q+(h_O+h_F)\left[\begin{array}{l}(c_I+c_E){\dot{u}}_{ry}-(c_I+c_E){\dot{u}}_{qy}-(c_I+c_E)(h_O+h_F){\dot{q}}_q+\\ +(c_I{\dot{q}}_D)u_{rz}+k_r{u}_{ry}-(c_I{\dot{q}}_D)u_{qz}-k_r{u}_{qy}-k_r(h_O+h_F)q_q\end{array}\right]$$ (30)

4 NUMERICAL RESULTS

Cases considered. The above presented methodology will be used to analyze the dynamic behavior of a rotor-foundation-structure soil system. Three cases will be considered and are illustrated in Figure 9. In the first case (a) the rotor is supported by rigid bearings. For the second case (b) the rotor, the foundation and the structure are supported by a rigid base. In the third case (c) the soil is included as the supporting medium. In all cases, the system is excited by the torque T applied directly to the rotor. In the examples that follow, the frame has only one element with 6 DOFs.

Figure 9
Cases for rotor-foundation-frame-soil systems.

System parameters. The mass, stiffness, and damping matrices of the frame structure, are presented in Table 1. The value for parameter α =1/50 (Cook et al., 2005).

Table 1
Frame parameters.

Young’s Module	Density	Poisson’s Ratio	Length	Area (Cross section)	Moment of inertia
E = 200 GPa	ρf = 7850 kg/m3	υf = 0.3	L = 2.5 m	A = 0.005 m2	Izz = 2.5E-5 m4

Similarly to the frame, the parameters of the Laval rotor, summarized in Table 2. These remain constant across all 3 cases. For the structural frame steel is assumed. The rotor has a mass of m_r = 50 kg, a stiffness of k_r = 904,778.7 N/m, and a damping coefficient of c_r = 190.9 N·s/m. The internal damping coefficient c_i=0. The natural frequency of the rotor on rigid bearings is calculated to be ω_r = 134.4 rad/s.

Table 2
Rotor parameters.

Young’s Module- Shaft	Density	Poisson’s Ratio	Mass Imbalance
E = 200 GPa	ρr = 7850 kg/m3	υr = 0.3	me = 0.048 kg.m

For the homogeneous half-space, material properties approximating those of real soil were adopted, as listed in Table 3. The corresponding flexibility profiles for vertical, horizontal, and rocking displacements are shown in Figure 6.

Table 3
Soil properties.

Young’s Module	Density	Poisson’s Ratio	Shear Modulus	Damping
E = 234 MPa	ρS = 2500 kg/m3	υS = 0.3	G = 90 MPa	η = 0.01

To nondimensionalize the analyses, mass ratios are defined between the rotor and its foundation (f_mRF), between the rotor foundation and the frame (f_mFP), and between the rotor and the frame (f_mRP). Likewise, stiffness ratios are defined between the rotor and the frame (f_kRP) and between the frame and the soil (f_kPS), as shown in equation (31).

$$\frac{m_R}{m_F}=f_{m_{RF}};\frac{m_F}{m_P}=f_{m_{FP}};\frac{m_R}{m_P}=f_{m_{RF}}{f}_{m_{FP}}=f_{m_{RP}};\frac{k_R}{k_P}=f_{k_{RP}};\frac{k_P}{k_S}=f_{m_{PS}}$$ (31)

The parameter values used in the analyses for masses and stiffnesses are summarized in Table 4.

Table 4
Adimensional factors.

Mass factor	Stiffness factor
fmRF = 0.05	fkRP = 0.032 fkPS = 107
fmFP = 0.21
fmRP = 0.01

Dynamics of the structure and soil-structure. The frequency response functions (FRFs) of the frame subsystem with a fixed base (case b) and those of the frame subsystem on a homogeneous half-space (case c) are presented in Figures 10 and 11, respectively. For the case b, the frame has only 3DOFs, shown in Figure 10 for the vertical, horizontal and rocking DOFs. System c has 6DOFs, corresponding to the two nodes of the frame. The resulting FRFs are shown in Figures 11.

Figure 10
FRFs of the frame attached to a fixed base (case b).

Figure 11
FRFs of the frame attached to a half space (case c).

Extraction of the modal parameters. The modal parameters extracted from the FRFs of the fixed-base frame (Figures 10) and those obtained from the FRFs of the frame supported by the half-space (Figure 11) are summarized in Table 5.

Table 5
Modal parameters.

Frame over fixed base
Natural frequency (rad/s)	ω1 = 15.89		ω2 = 100.10		ω3 = 338.93
Damping factor	ξ1 = 0.0200		ξ2 = 0.0087		ξ3 = 0.0200
Frame attached to the half space
Natural frequency (rad/s)	ω1 = 11.96	ω2 = 38.54	ω3 = 49.46	ω4 = 91.76	ω5 = 146.20	ω6 = 480.42
Damping factor	ξ1 = 0.0382	ξ2 = 0.1309	ξ3 = 0.1242	ξ4 = 0.0735	ξ5 = 0.1554	ξ6 = 0.0362

The extracted modal forms of frame subsystem of case b, Figure 10, are shown in Figure 12.

Figure 12
Modal forms – Frame attached to rigid base (case b).

The modal forms extracted from the FRFs of the frame over the half-space Figure 11, case c, are shown in Figure 13.

Figure 13
Modal forms – Frame attached to the half space (case c).

The calculations of the natural frequencies of rigid bodies are performed as defined in Gazetas (1983), in which the static stiffness for each degree of freedom of the soil is defined in equation (32). The natural frequencies of rigid body are calculated following the methodology proposed by Gazetas (1983) where the static stiffness associated with each degree of freedom of the soil is defined by equation (32).

$$k_z=\frac{4GR}{(1-\nu )};k_y=\frac{8GR}{(2-\nu )};k_{\theta }=\frac{8G{R}^3}{3(1-\nu )}$$ (32)

where, k_z denotes the static stiffness associated with vertical displacement, k_y corresponds to horizontal displacement stiffness, and k_θ represents the stiffness related to rocking motion. Using these static stiffness values along with the frame’s mass matrix, the natural frequencies for the frame rigid body were calculated as ω_nz = 38.75 rad/s, ω_ny = 35.17 rad/s, and ω_nθ = 70.34 rad/s.

4.1 Response of Rotor on Rigid base – case a

One important issue that will be investigated using the proposed methodology is the dynamics of the rotor during the runup phase. In this process an external torque T is applied to the rotor and it starts accelerating. During this acceleration process the unbalance excitation frequency increases. In the runup process the rotor passes through the resonances present in the system. Close to the resonances the displacement amplitude of the rotor, of the foundation, of the structure and soil DOF does increase. The energy dissipated by the vibration of the rotor through the viscous damping mechanism/coefficient c_E will increase. If the energy dissipated by the damping mechanisms is equal or larger to the energy that the torque I inputs in the system, the rotor cannot pass through the resonances. It ‘hangs’ on one resonance frequency with large vibration amplitudes. This is a very undesirable situation. The minimum torque required for the rotor to pass through a resonance is called in this article ‘minimum torque’, T_min. In this article the ‘minimum torque’ required for the rotor to pass through the resonances will be studied for the three cases described in the previous section.

Based on the parameters adopted in this study, the torque limit for case a, rotor on rigid bearings, is T = 17.8 Nm. This value corresponds to the threshold, below which the system ‘hangs’ at an angular velocity of approximately 132.91 rad/s. Figure 14, which presents the rotor displacements u_zr and u_yr as a function of time (blue curve). The angular velocity ${\dot{\theta }}_D$ as functions of time is given by the red curve. As it can be seen after a period of angular acceleration, the rotor hangs at a vibrating frequency close to its natural frequency of 134.4 rad/s. For this case, the maximum displacement rotor amplitude is 0.027 m.

Figure 14
Case a: Rotor on Rigid Bearings – T = 17.8 Nm.

When the applied torque is increased to T=18 Nm, the system no longer ‘hangs’ at an angular velocity near the natural frequency, as can be seen in Figure 15. The rotor is able to pass through the resonance. It should be noticed, however, that there is a considerable reduction in the rotor angular acceleration, when the rotor passes through the resonance. It experiences large vibration amplitudes in for many cycles in both directions, u_rz and u_ry.

Figure 15
Rotor DOFs u_rz and u_ry. Case a: Rotor on Rigid Bearings – T = 18 Nm.

Figure 16 presents the rotor responses for a higher applied torque of T=28 Nm. In this case, the reduction in rotor acceleration near the resonance region becomes nearly imperceptible; however, the response amplitudes in this region remain clearly increased. The rotor passes more quickly through the resonance.

Figure 16
Rotor DOFs u_rz and u_ry. Case a: Rotor on Rigid Bearings – T = 28 Nm.

4.2 Rotor-foundation-structure on a rigid base – case b

Rotor response. For case b, rotor-foundation on a frame structure, the limiting torque is found to be T=16.9 Nm. As shown in Figure 17, the system stagnates at an angular velocity of approximately 134.06 rad/s, with an average maximum displacement amplitude of 0.025 m. The inclusion of the frame mounted on a rigid base under the rotor-foundation system results in a reduction of both the response amplitude and the limiting torque. It should be considered that the frame structure has in internal damping coefficient c_t and thus helps to dissipate the energy in the system. In an initial analysis it could be considered that this extra withdrawal of energy through the structural damping c_t would require an increased torque to overcame the resonance. But it should be reminded that viscous dampers dissipate energy directly proportional to the velocity of the displacements. In the case of a more damped system, the vibration amplitudes near the resonance are considerably reduced, as is the velocity, reducing considerably the dissipated energy. With less energy dissipated, a smaller torque is required to pass the resonance.

Figure 17
Rotor DOFs u_rz and u_ry. Case b: rotor-foundation-structure on rigid base – T = 16.9 Nm.

When the applied torque is increased to T = 19 Nm, the system no longer stagnates at an angular velocity near the natural frequency, as shown in Figure 18, similar to the behavior observed in case a. A significant increase in response amplitude in both directions is observed in the vicinity of the resonance region, although the reduction in rotor acceleration is less pronounced.

Figure 18
Rotor DOFs u_rz and u_ry. Case b: rotor-foundation-structure on rigid base – T = 19 Nm.

Figure 19 presents the vertical u_tz2 displacement response of the frame at node 2 for the applied torque of T = 19 Nm. Although the amplitudes are significantly lower than those observed in the rotor, an increase in response amplitude is evident in the rotor’s resonance region. Additionally, a second increase in amplitude is observed when the angular velocity approaches the natural frequency of the frame over rigid base, approximately 300 rad/s, as indicated in Table 5.

Figure 19
Structure vertical DOF u_tz2. Case b: rotor-foundation-structure on rigid base – T = 19 Nm.

A similar behavior is observed in the horizontal u_ty2 and rotation θ_tx2 displacements of the foundation node 2, as shown in Figure 20. However, in this case, the second increase in amplitude occurs prior to the one associated with the rotor's resonance region. This earlier amplification is observed near the second natural frequency of the frame over rigid base, approximately 110 rad/s, as listed in Table 5.

Figure 20
Structure DOF u_ty2 and θ_tx2. Case b: rotor-foundation-structure on rigid base – T = 19 Nm.

4.3 Rotor-foundation-structure on soil – case c

For case c, the rotor, the foundation and the structure supported by the soil, the torque limit is T=16.5 Nm, the lowest among the three cases analyzed. As shown in Figure 21, the system ‘hangs’ at an angular velocity of approximately 134.6 rad/s, with an average maximum displacement amplitude of 0.024 m. The inclusion of the soil beneath the frame introduces more damping through the geometric damping mechanism and leads to a reduction in both the response amplitude and the limiting torque in comparison with case b.

Figure 21
Rotor DOFs u_rz and u_ry. Case c: rotor-foundation-structure on soil – T = 16.5 Nm.

As in the previous cases, increasing the applied torque to T = 19 Nm results in a similar behavior: the rotor exhibits nearly constant acceleration, except in the vicinity of its resonance region, where a significant increase in response amplitude is observed, Figure 22.

Figure 22
Rotor DOFs u_rz and u_ry. Case c: rotor-foundation-structure on soil – T = 19 Nm.

Figure 23 shows the displacement response u_tz2 of the frame at node 2 under an applied torque of T=19 Nm. In this case, in addition to the amplitude increase near the rotor’s resonance region, two additional peaks of smaller amplitude are observed as the angular velocity approaches the second and sixth natural frequencies of the soil-structure subsystem, approximately at 47 rad/s and 460 rad/s, respectively, as indicated in Table 5. The rotor passes through three distinct resonances.

Figure 23
Structure DOF u_tz2. Case c: rotor-foundation-structure on soil – T = 19 Nm.

A similar behavior is observed for the frame displacements u_ty2 and θ_tx2, as illustrated in Figure 24. However, in this case, two amplitude peaks appear prior to the rotor’s resonance region. These amplifications occur near the third and fourth natural frequencies of the soil–structure subsystem, approximately at 59 rad/s and 98 rad/s, respectively, as listed in Table 5. It is worth noting that the peak near 59 rad/s appears only in the u_ty2 response.

Figure 24
Structure DOFs u_ty2 and θ_tx2. Case c: rotor-foundation-structure on soil – T = 19 Nm.

What should be stressed is not the important and interesting results of the dynamic response of these complex systems per se, but the synthesis of a methodology that allows to perform these transient analyses in systems that contain non-linearities (rotor) and present geometric damping (soil).

5 CONCLUSIONS

The present paper introduced a methodology to analyze the transient behavior of complex systems constituted of a rotor, a foundation, a structure and a soil. The system is divided into two subsystems. The time domain non-linear equations of a Laval rotor supported by a rigid foundation are derived and form the subsystem I. These non-linear equations of motion are integrated using a 4th order Runge Kutta method. A frame structure resting on a viscoelastic half-space constitutes the subsystem II. The soil which is considered an unbounded domain and presents geometric damping is modelled by a 3D version of the Boundary Element Method. The BEM can naturally describe the dynamics of complex soil profiles and foundation systems. A new methodology to construct the time domain equations of a system that is dynamically equivalent to the original soil-structure system is presented. The methodology is based on extraction of modal parameters from the Frequency Response Functions of the coupled structure-soil problem. This equivalent system in modal coordinates can be numerically integrated by standard procedures. The response of the complete system is obtained by a staggered coupling procedure, with the rotor and the soil-structure systems being integrated by distinct strategies.

This methodology enables an accurate solution of a non-linear system coupled with another system presenting geometric or radiation damping, related to the dynamics of unbounded domains. The methodology renders transient response with very small time steps, allowing it to be applied to problems like the transient response of an unbalanced rotor interacting with a foundation-structure-soil system during the runup phase. The examples given in the article show that the method is able to capture the transient behavior of the rotor near the resonances. The proposed method allows to determine the minimum torque required to be applied to the rotor in order to keep it from ‘hanging’ at resonance frequencies. The methodology also is accurate enough to enable the study of the influence of the damping mechanisms existing on the rotor, structure and soil on the unbalance response of the former. In short, the presented methodology opens the possibility to numerically investigate, with great accuracy, the transient behavior of complex coupled systems as exemplified by onshore eolic systems.

Appendix

Time domain non-linear equations of the rotor-foundation system.

Displacement vector:

$$\left\{u\right\}={\left\{\begin{array}{cc}\begin{array}{ccc}{u}_{rz}& u_{ry}& {\theta }_D\end{array}& \begin{array}{ccc}{u}_{qz}& u_{qy}& {\theta }_q\end{array}\end{array}\right\}}^T$$ (A1)

Mass matrix [M_I] equation (16):

$$\left[M_I\right]=\left[\begin{array}{cc}\begin{array}{ccc}{m}_D& 0& 0\\ 0& m_D& 0\\ 0& 0& I_{xx}^{\raisebox{1ex}{$R$}\!\left/ \!\raisebox{-1ex}{$D$}\right.}\end{array}& \begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\\ \begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}& \begin{array}{ccc}{m}_F& 0& 0\\ 0& m_F& m_F{h}_F\\ 0& m_F{h}_F& m_F{h}_F^2+I_{xx}^{\raisebox{1ex}{$R$}\!\left/ \!\raisebox{-1ex}{$F$}\right.}\end{array}\end{array}\right]$$ (A2)

Matrix [A_I] and [B_I] equation (16):

$$\left[A_I\right]=\left[\begin{array}{cc}\begin{array}{ccc}\begin{array}{c}\begin{array}{c}(c_E+c_I)\\ 0\\ -(c_E+c_I)e\cos {\theta }_D\end{array}\\ \begin{array}{c}-(c_E+c_I)\\ 0\\ 0\end{array}\end{array}& \begin{array}{c}\begin{array}{c}0\\ (c_E+c_I)\\ (c_E+c_I)esen{\theta }_D\end{array}\\ \begin{array}{c}0\\ -(c_E+c_I)\\ -(c_E+c_I)(h_O+h_F)\end{array}\end{array}& \begin{array}{c}\begin{array}{c}0\\ 0\\ c_I{u}_{est}esen{\theta }_D\end{array}\\ \begin{array}{c}0\\ 0\\ 0\end{array}\end{array}\end{array}& \begin{array}{ccc}\begin{array}{c}\begin{array}{c}-(c_E+c_I)\\ 0\\ (c_E+c_I)e\cos {\theta }_D\end{array}\\ \begin{array}{c}(c_E+c_I)\\ 0\\ 0\end{array}\end{array}& \begin{array}{c}\begin{array}{c}0\\ -(c_E+c_I)\\ -(c_E+c_I)esen{\theta }_D\end{array}\\ \begin{array}{c}0\\ (c_E+c_I)\\ (c_E+c_I)(h_O+h_F)\end{array}\end{array}& \begin{array}{c}\begin{array}{c}0\\ (c_E+c_I)(h_O+h_F)\\ -(c_E+c_I)(h_O+h_F)esen{\theta }_D\end{array}\\ \begin{array}{c}0\\ (c_E+c_I)(h_O+h_F)\\ (c_E+c_I){(h_O+h_F)}^2\end{array}\end{array}\end{array}\end{array}\right]$$ (A3)

$$\left[B_I\right]=\left[\begin{array}{cc}\begin{array}{ccc}\begin{array}{c}\begin{array}{c}{k}_r\\ (c_I{\dot{\theta }}_D)\\ -k_{r}e\cos {\theta }_D\end{array}\\ \begin{array}{c}-k_r\\ -(c_I{\dot{\theta }}_D)\\ -(c_I{\dot{\theta }}_D)(h_O+h_F)\end{array}\end{array}& \begin{array}{c}\begin{array}{c}-(c_I{\dot{\theta }}_D)\\ k_r\\ k_{r}esen{\theta }_D\end{array}\\ \begin{array}{c}(c_I{\dot{\theta }}_D)\\ -k_r\\ -k_r(h_O+h_F)\end{array}\end{array}& \begin{array}{c}\begin{array}{c}0\\ 0\\ 0\end{array}\\ \begin{array}{c}0\\ 0\\ 0\end{array}\end{array}\end{array}& \begin{array}{ccc}\begin{array}{c}\begin{array}{c}-k_r\\ -(c_I{\dot{\theta }}_D)\\ k_{r}e\cos {\theta }_D\end{array}\\ \begin{array}{c}{k}_r\\ (c_I{\dot{\theta }}_D)\\ (c_I{\dot{\theta }}_D)(h_O+h_F)\end{array}\end{array}& \begin{array}{c}\begin{array}{c}(c_I{\dot{\theta }}_D)\\ -k_r\\ -k_{r}esen{\theta }_D\end{array}\\ \begin{array}{c}-(c_I{\dot{\theta }}_D)\\ k_r\\ k_r(h_O+h_F)\end{array}\end{array}& \begin{array}{c}\begin{array}{c}(c_I{\dot{\theta }}_D)(h_O+h_F)\\ -k_r(h_O+h_F)\\ -k_r(h_O+h_F)esen{\theta }_D\end{array}\\ \begin{array}{c}-(c_I{\dot{\theta }}_D)(h_O+h_F)\\ k_r(h_O+h_F)\\ k_r{(h_O+h_F)}^2\end{array}\end{array}\end{array}\end{array}\right]$$ (A4)

Vector {F_I} equation (16):

$$\left\{F_I\right\}=\left\{\begin{array}{c}\begin{array}{c}{m}_D(-{\ddot{\theta }}_{D}e\cos {\theta }_D+{\dot{\theta }}_D^{2}esen{\theta }_D)\\ m_D({\ddot{\theta }}_{D}esen{\theta }_D+{\dot{\theta }}_D^{2}e\cos {\theta }_D)\\ c_{I}esen{\theta }_D{\dot{\theta }}_D(-u_{rz}+u_{sz})+c_{I}e\cos {\theta }_D{\dot{\theta }}_D(-u_{ry}+u_{sy}+(h_O+h_F){\theta }_q)+k_{r}e\cos {\theta }_D{u}_{est}+T\end{array}\\ \begin{array}{c}-F_{qz}\\ -F_{qy}\\ -M_{qx}\end{array}\end{array}\right\}$$ (A5)

6 ACKNOWLEDGEMENTS

The research leading to this work was funded by the São Paulo Research Foundation (Fapesp) grant FAPESP CEPID Process 2013/08293-7 and Coordination for the Improvement of Higher Education Personnel - Brazil (CAPES) - Finance Code 001. The support of CNPq and Faepex/Unicamp is also gratefully acknowledged.

• https://doi.org/10.1590/1679-7825/e8784
• Data availability statement:
  Research data is only available upon request.

References

• Carrion, R., Sousa, A.D.O. and Mesquita, E, 2007. “The influence of soil damping mechanisms and geo-profiles on the stationary response of 3d rigid block foundations. In Proceeding of COBEM 2007, 19th International Congress of Mechanical Engineering, November 5-9, 2007, Brasília, DF, pp. 01-10.
• Caughey, T. K. Classical normal modes in damped linear dynamic systems. Journal of Applied Mechanics (ASME), v. 27, p. 269–271, 1960.
• Cheng, D. K. Analysis of Linear Systems - 1st ed. [S.l.]: Addison-Wesley, 1972. 431 p.
• Chopra, A. K. Dynamics of Structures: Theory and Applications to Earthquake Engineering - 4th ed. [S.l.]: Bookman Editora, 2012. 944 p.
• Cook, R. D., Malkus, D. S., Plesha, M. E., and Witt, R. J., 2005, Concepts and Applications of Finite Element Analysis, 4th ed., John Wiley and Sons, New York.
• Dominguez J. “Boundary Element Methods in Dynamics”. Wit Press, 1992.
• Ewins, D.J. “Modal testing: theory, practice and application”. Wiley, New York 2000.
• Felippa C.A. and Park K.C. "Staggered transient analysis procedures for coupled mechanical systems: formulation". Computer Methods in Applied Mechanics and Engineering, Vo 24, Issue 1, pp. 61-111, 1980.
• Ferraz, A. C.; Pacheco, L. A.; Carrion, R. and Mesquita, E. “Coupling Modal Analysis with the BEM for the Transient Response of Bar Structures Interacting with Three-Dimensional Soil Profiles”. Latin American Journal of Solids and Structures, Vol. 20 No. 6, 2023.
• Filho, H.A.N.C., Avila, S.M. & de Brito, J.L.V. “Dynamic analysis of onshore wind turbines including soil–structure interaction”. J Braz. Soc. Mech. Sci. Eng. 43, 143 (2021). https://doi.org/10.1007/s40430-021-02837-5
  » https://doi.org/10.1007/s40430-021-02837-5
• Gasch, R.; Maurer, J. and Sarfeld, W. “Soil influence on unbalance response and stability of a simple rotor-foundation system”. Journal of Sound and Vibration 93(4), pp. 549–566, 1984.
• Gazetas, G. Analysis of machine foundation vibrations: state of the art. Soil Dynamics and Earthquake Engineering, v. 2, p. 02-41, 1983.
• Jahani, K., Langlois, R.G, Afagh, F.F.: Structural dynamics of offshore Wind Turbines: A review. Ocean Engineering, Volume 251, 1 May 2022, paper 111136
• Kellezi, L.; Hansen, P. B. “Dynamic transient analysis of a mono-pile windmill foundation”. Proceedings of EURODYN, 4th Int. Conference on Structural Dynamics. Munich, Germany. pp. 1315-1320, 2002.
• Liu, F.; Li, X.; Tian, Z.; Zhang, J.; Wang, B. Transient Response Estimation of an Offshore Wind Turbine Support System. Energies 2019, 12, 891.
• Manolis G. D. and Beskos D. E. “Boundary Element Methods in Elastodynamics”. Unwin Hyman Ltd., London, 1988.
• Mesquita, E.; Cavalca, K. L.; Ramalho D.A.; Thomazo, L. H. On the Influence of Soil Stratification Profiles on Rotor-Foundation Unbalance Response. Proceedings 7th Iftomm International Conference on Rotor Dynamics, 2006. v. 1. p. 01-10.
• Pacheco, L.A., Ferraz, A. C., Mesquita, E., Lima, L.F.V. “Transient response of bar structures interacting with layered soil profiles”. In Proceedings of the 9th International Symposium on Solid Mechanics, MECSOL2024, P.T.R.de Mendonça, E. A. Fancello and C. R. M. Roesler (Editors), Florianópolis SC, Brazil, October 21st to 23th, 2024.

Edited by

• Editor:
  Eduardo Alberto Fancello and Paulo de Tarso Mendonça

Data availability

Research data is only available upon request.

Publication Dates

• Publication in this collection
  01 Dec 2025
• Date of issue
  2026

History

• Received
  31 July 2025
• Reviewed
  17 Sept 2025
• Accepted
  01 Oct 2025