Dynamic analysis of the tensegrity structure of a rotary airborne wind energy machine

The dynamic behavior of the tensegrity structure (helix) of a rotary airborne wind energy (RAWE) machine was investigated by combining experimental and numerical techniques. Taking advantage of the slenderness of the helix, a dynamic model for the evolution of its center line and the torsional deformation was developed by using Cosserat theory. The constitutive relations for the axial, bending, and torsional stiffness, which are a fundamental component of the model, were obtained experimentally by carrying out laboratory tests. Three scenarios of increasing complexity were then studied with the numerical tool. Firstly, a stationary solution of the model, i.e., with fixed ends and no rotation, was found numerically and used to verify the correct implementation of a numerical code based on finite elements. The stability analysis of this solution, which corresponds to the state of the structure just after deployment but before operation, showed that the natural periods of longitudinal, lateral, and torsional modes of the RAWE structure under consideration are around 0.03, 0.2, and 0.4 s, respectively. Secondly, the dynamics in nominal operation was investigated by keeping both end tips fixed and implementing a controller that adjusts the torque at the ground to reach a target angular velocity of 120 rpm. Key characteristic variables like the tension and the response times of the helix were obtained. Thirdly, the dynamics of the helix when the lower end is fixed and the upper end is driven in a circular motion of frequency f₁ was studied experimentally and numerically. The tension of the helix in the experiment increased for f₁ above a certain threshold, and the structure collapsed at f₁≈5 Hz. Simulation analysis revealed a resonance of the structure at a higher frequency (around 13 Hz).

1 Introduction

The increase in wind power density with altitude (Archer and Caldeira, 2009) constitutes an important driver for wind technologies. Conventional wind turbines have notably increased their size during the last few decades to reach higher winds. However, although control strategies have been developed to meet efficiency and reliability requirements (Njiri and Soffker, 2016), structural considerations may set an upper limit. Airborne wind energy systems (AWESs), mainly based on soft kites, rigid wings, or rotors linked to the ground by tethers, can harvest energy at high altitudes by using the tether tension or onboard wind turbines (Schmehl, 2018; Malz et al., 2022). In ground-gen AWE concepts the conversion from mechanical to electrical energy happens on the ground, whereas in fly-gen systems such conversion is done on the aircraft (Cherubini et al., 2015). Some of these technologies are in a precommercial state (European Commission, 2018; Blanch et al., 2022). Although it is still unclear which of them will dominate the market among the plethora of existing architectures (van de Kaa and Kamp, 2021), system reliability, operational robustness, and safety are crucially important aspects of system development (Salma et al., 2020) together with their performances in realistic vertical wind velocity profiles (Sommerfeld et al., 2019; Schelbergen et al., 2020; Sommerfeld et al., 2023).

An interesting subfamily of AWE systems involves concepts based on rotating kites or rotors that use the well-known phenomenon of auto-rotation for harvesting wind energy (Rimkus and Das, 2013). An example is a tethered autogyro system with four rotors in a quadrotor configuration that harvests energy and transmits it to the ground via the tether (Mackertich and Das, 2016). Rotating reel parotors, which combine rotary ring kites with ground-based rotating reel conversion systems (Benhaiem and Schmehl, 2018) and gyrocopter-type airborne wind turbines (Rancourt et al., 2016; Roberts et al., 2007), have been proposed. Prototypes of AWE systems aimed at the direct transmission of the mechanical torque produced by a set of kites or a flying rotor to a generator on the ground have also been manufactured and tested. The Daisy stack by Windswept and Interesting Ltd (Read, 2018) and the rotary AWE (RAWE) machine of SomeAwe (Beaupoil, 2017) belong to this category. These two concepts have in common the use of a light tensegrity structure (Motro, 2003) to transmit the torque to the generator and an auxiliary kite to provide extra lift and fix the elevation angle of the machine.

Flight testing activities with the RAWE machine showed that the tensegrity structure exhibits a rich dynamics that involves longitudinal, lateral, and torsional waves. The tensegrity structure has a helix-like shape that acquires a high torsional and bending stiffness when a traction load is applied. Finding the operational limits of the helix is crucial for the reliability of the RAWE machine. According to experimental tests, the structure can collapse if the torque is above a certain threshold (Beaupoil, 2022). For this reason, characterizing the axial, bending, and torsional stiffness of the helix, as well as having access to numerical tools capturing longitudinal, lateral, and torsional dynamics, is important. Numerical simulations can help to predict the operational boundaries and avoid hardware failure. Nevertheless, the modeling of the structure represents an important challenge that should balance fidelity and computational cost. Two dynamic models of different complexity have been developed in a previous study for the tensegrity structure of the Daisy stack (Tulloch et al., 2020, 2023): a simple spring-disk model and a multi-spring/multi-punctual mass model. The latter involves a higher number of degrees of freedom and can capture the variation in axial tension along the length of the structure. Interestingly, two stable-equilibrium states were predicted by the numerical simulations when the rotor dynamics was coupled with the spring-disk model (Tulloch et al., 2020). A steady-state model that uses aerodynamics and structural modules to analyze the performance of RAWE machines has been proposed (Wacker et al., 2023). The control problem was also studied by using a model with a small number of degrees of freedom (De Schutter et al., 2018), and the aerodynamics of the rotor was studied with blade-element theory and vortex computations (Pfister and Blondel, 2020). To the best of our knowledge there is a lack of knowledge about the axial, bending, and torsional stiffness of existing RAWE machines.

An important property of the tensegrity structures of RAWE machines is that they are slender; i.e., their unstretched length (L₀) and characteristic radius (R_s) satisfy L₀≫R_s. This work exploits this key feature to model the structure by using Cosserat theory that, due to its versatility, has been used in a very broad range of applications (Ramézani et al., 2011; Riahi and Curran, 2009). In this case, we use Cosserat rod theory to find a set of partial differential equations that govern the dynamics of the central line and its torsional deformation (see Sect. 2). This set of equations is appropriate for investigating in detail the rich dynamics observed in RAWE experiments because it captures the stretching, bending, and twisting of the structure. A finite-element method is proposed in Sect. 3 to approximate the partial differential equations by a set of ordinary differential equations. The axial, bending, and torsional stiffness, which naturally appear in the dynamic equations, have been determined experimentally in Sect. 4 for the RAWE machine of SomeAwe. Three scenarios of increasing complexity are studied numerically and experimentally in Sect. 5. Section 6 summarizes the conclusions. The code presented in this work is part of the open-source software LAKSA (Sánchez-Arriaga et al., 2017, 2019).

2 A model for rotary AWE machines

2.1 Kinematics considerations

Figure 1 shows a RAWE machine made of a tensegrity structure (helix) of unstretched length L₀ that transmits the torque generated by a set of blades to a generator on the ground. The slenderness of the tensegrity structure, i.e., with its length being much larger than the characteristic length of its cross section, suggests using some approximations that help to simplify the mathematical analysis. In particular, we use Cosserat theory and only compute the evolution of the position vector of the central line of the helix and an additional angle that orients the cross section, which is assumed to be unshearable. These two variables are called r(s,t) and α(s,t), where $\mathrm{0}\le s\le L_{\mathrm{0}}$ and t are the unstretched arc length and time, respectively. Taking the helix as a one-dimensional Cosserat rod with an unshearable cross section has the advantage of extraordinarily simplifying the model while allowing us to capture the axial deformation, the bending of the helix, and the twisting of the helix. An inertial frame of reference 𝒮_E with its origin at the generator, plane x_E−y_E spanning the ground, and the z_E axis pointing downwards are introduced. We call i_E, j_E, and k_E to the unit vectors along the three axes of frame 𝒮_E. Therefore, $\mathit{r}(\mathrm{0},t)=\mathrm{0}$ is the origin of 𝒮_E, and $\mathit{r}(L_{\mathrm{0}},t)={\mathit{r}}_T$ is the position vector of the hub of the blades (point T in Fig. 1). In this work, we use bold-italic characters to denote vectors, and bold-roman characters to denote matrices.

Figure 1 Sketch and frames of reference of the rotary AWE machine.

Following Buckham (2003), we define at every point of the center line a Frenet frame (𝒮_F) with axes spanned by the unit vectors ($\mathit{t},\mathit{n},\mathit{b}$). For small axial deformations of the helix

$$\begin{array}{}\text{(1)}& \left|{\displaystyle \frac{\mathrm{d}\mathit{r}}{\mathrm{d}s}}\right|-\mathrm{1}\ll \mathrm{1},\end{array}$$

the vectors of the Frenet basis read

$$\begin{array}{}\text{(2)}& \mathit{t}={\mathit{r}}^{\prime },\mathit{n}={\displaystyle \frac{{\mathit{r}}^{\prime \prime }}{\mathit{\kappa }}},\mathit{b}=\mathit{t}\times \mathit{n}.\end{array}$$

Although the same character is used in this work, the tangent vector t should not be confused with time t, which is a scalar and always written with a plain letter. The curvature and torsion of the center line are given by

$$\begin{array}{}\text{(3)}& \mathit{\kappa }={\left({\mathit{r}}^{\prime \prime }\cdot {\mathit{r}}^{\prime \prime }\right)}^{\mathrm{1}/\mathrm{2}},\mathit{\gamma }={\displaystyle \frac{{\mathit{r}}^{\prime }\cdot \left({\mathit{r}}^{\prime \prime }\times {\mathit{r}}^{\prime \prime \prime }\right)}{{\mathit{r}}^{\prime \prime }\cdot {\mathit{r}}^{\prime \prime }}},\end{array}$$

where we denoted by the superscript ^′ the derivative with respect to the unstretched arc length s. The derivatives of this basis with respect to the arc length are

$$\begin{array}{}\text{(4)}& {\mathit{t}}^{\prime }=\mathbf{\Omega }\times \mathit{t},{\mathit{n}}^{\prime }=\mathbf{\Omega }\times \mathit{n},{\mathit{b}}^{\prime }=\mathbf{\Omega }\times \mathit{b},\end{array}$$

where $\mathbf{\Omega }=\mathit{\kappa }\mathit{b}+\mathit{\gamma }\mathit{t}$ is the Darboux vector.

At every point of the center line we also define a local frame 𝒮_D with director vectors $({\mathit{d}}_{\mathrm{1}},{\mathit{d}}_{\mathrm{2}},{\mathit{d}}_{\mathrm{3}})$ and origin at the particular point of the center line. The director vector d₃ points along the tangent vector t. Vectors d₁ and d₂ span the cross section of the helix and are directed along its principal axes of inertia. The director basis is related to the Frenet frame by a rotation of an angle α about the tangent vector t

$$\begin{array}{}\text{(5)}& \left(\begin{array}{c}{\mathit{d}}_{\mathrm{1}}\\ {\mathit{d}}_{\mathrm{2}}\end{array}\right)=\left(\begin{array}{cc}\cos \mathit{\alpha }& -\sin \mathit{\alpha }\\ \sin \mathit{\alpha }& \cos \mathit{\alpha }\end{array}\right)\left(\begin{array}{c}\mathit{n}\\ \mathit{b}\end{array}\right).\end{array}$$

The derivative of the director basis along the center line is

$$\begin{array}{}\text{(6)}& {\mathit{d}}_{\mathrm{1}}^{\prime }=\mathit{K}\times {\mathit{d}}_{\mathrm{1}},{\mathit{d}}_{\mathrm{2}}^{\prime }=\mathit{K}\times {\mathit{d}}_{\mathrm{2}},{\mathit{d}}_{\mathrm{3}}^{\prime }=\mathit{K}\times {\mathit{d}}_{\mathrm{3}},\end{array}$$

where

$$\begin{array}{}\text{(7)}& \mathit{K}=\mathbf{\Omega }+{\displaystyle \frac{\partial \mathit{\alpha }}{\partial s}}{\mathit{d}}_{\mathrm{3}}=\mathit{\kappa }\mathit{b}+\mathit{\tau }\mathit{t}\end{array}$$

is the twist vector, and it involves the curvature κ and torsion γ of the center line. The scalar

$$\begin{array}{}\text{(8)}& \mathit{\tau }\equiv \mathit{\gamma }+{\displaystyle \frac{\partial \mathit{\alpha }}{\partial s}}\end{array}$$

is the total twist of the principal axes of inertia of the helix about the tangent vector t. The time evolution of the director vectors is given by an equation analogous to Eq. (6):

$$\begin{array}{}\text{(9)}& {\dot{\mathit{d}}}_{\mathrm{1}}={\mathit{\omega }}_{DE}\times {\mathit{d}}_{\mathrm{1}},{\dot{\mathit{d}}}_{\mathrm{2}}={\mathit{\omega }}_{DE}\times {\mathit{d}}_{\mathrm{2}},{\dot{\mathit{d}}}_{\mathrm{3}}={\mathit{\omega }}_{DE}\times {\mathit{d}}_{\mathrm{3}},\end{array}$$

where the dot denotes the time derivative, and ω_DEis the angular velocity of the local frame (𝒮_D) with respect to the inertial frame (𝒮_E). Vectors ω_DE and K are linked by the constraint

$$\begin{array}{}\text{(10)}& {\displaystyle \frac{\partial {\mathit{\omega }}_{DE}}{\partial s}}-{\displaystyle \frac{\partial \mathit{K}}{\partial t}}={\mathit{\omega }}_{DE}\times \mathit{K},\end{array}$$

which is found by differentiating Eq. (6) with respect to time and Eq. (9) with respect to the arc length. The angular velocities of frames 𝒮_D and 𝒮_F with respect to 𝒮_E are related by ${\mathit{\omega }}_{DE}={\mathit{\omega }}_{FE}-\dot{\mathit{\alpha }}\mathit{t}$.

2.2 Dynamic model

Momentum and angular momentum equations of the structure with total mass M read (see for instance Villaggio, 2005)

$$\begin{array}{}\text{(11)}& {\displaystyle }\mathit{\rho }A\ddot{\mathit{r}}={\mathit{f}}^{\prime }+\mathit{\rho }Ag{\mathit{k}}_E+{\mathit{F}}_{\mathrm{A}},\text{(12)}& {\displaystyle }\mathit{\rho }\left(I_{\mathrm{2}}{\ddot{\mathit{d}}}_{\mathrm{1}}\times {\mathit{d}}_{\mathrm{1}}+I_{\mathrm{1}}{\ddot{\mathit{d}}}_{\mathrm{2}}\times {\mathit{d}}_{\mathrm{2}}\right)={\mathit{m}}^{\prime }+{\mathit{r}}^{\prime }\times \mathit{f},\end{array}$$

where F_A is the aerodynamic force per unit length, and f and m are the internal force and torque, respectively. We also introduced the density $\mathit{\rho }=M/L_{\mathrm{0}}A$ and the two principal moments of the cross section of area A:

$$\begin{array}{}\text{(13)}& I_{\mathrm{1}}=\int \underset{A}{\int }{x}_{\mathrm{2}}^{\mathrm{2}}\mathrm{d}{x}_{\mathrm{1}}\mathrm{d}{x}_{\mathrm{2}},I_{\mathrm{2}}=\int \underset{A}{\int }{x}_{\mathrm{1}}^{\mathrm{2}}\mathrm{d}{x}_{\mathrm{1}}\mathrm{d}{x}_{\mathrm{2}}.\end{array}$$

Regarding F_A, we only consider the aerodynamic drag and write

$$\begin{array}{}\text{(14)}& {\mathit{F}}_{\mathrm{A}}=-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\mathit{\rho }}_{\mathrm{0}}{l}_a{C}_{\mathrm{D}}{v}_{\mathrm{A}}{\mathit{v}}_{\mathrm{A}},\end{array}$$

where C_D and ρ₀ are the drag coefficient and the air density. We also defined the aerodynamic velocity (${\mathit{v}}_{\mathrm{A}}=\mathit{v}-{\mathit{v}}_{\mathrm{w}}$); $\mathit{v}\equiv \mathrm{d}\mathit{r}/\mathrm{d}t$; the wind velocity (v_w); and the characteristic transversal length (l_a), which can be approximated by the sum of the diameters of the bars of the helix. In our analysis we consider a constant wind velocity vector ${\mathit{v}}_{\mathrm{w}}=-v_{\mathrm{w}\mathrm{0}}{\mathit{i}}_E$.

The model is completed by the constituent relations that depend on the elastic properties of the material and on the shape and the dimension of the cross section of the structure. Writing the twist vector as $\mathit{K}=K_{\mathrm{1}}{\mathit{d}}_{\mathrm{1}}+K_{\mathrm{2}}{\mathit{d}}_{\mathrm{2}}+K_{\mathrm{3}}{\mathit{d}}_{\mathrm{3}}$, the internal moment is (Love, 1944)

$$\begin{array}{}\text{(15)}& \mathit{m}=E{I}_{\mathrm{1}}{K}_{\mathrm{1}}{\mathit{d}}_{\mathrm{1}}+E{I}_{\mathrm{2}}{K}_{\mathrm{2}}{\mathit{d}}_{\mathrm{2}}+GJ{K}_{\mathrm{3}}{\mathit{d}}_{\mathrm{3}},\end{array}$$

where E is Young's modulus, G is the rigidity modulus, and J is the torsion constant. Hereafter, we model the structure as a hollow cylinder of radius R_s and thickness h_s. Therefore, we have from Eq. (13)

$$\begin{array}{}\text{(16)}& I_{\mathrm{1}}=I_{\mathrm{2}}=I\equiv \mathit{\pi }{R}_{\mathrm{s}}^{\mathrm{3}}{h}_{\mathrm{s}},J=\mathrm{2}\mathit{\pi }{R}_{\mathrm{s}}^{\mathrm{3}}{h}_t,\end{array}$$

and, using Eq. (7), Eq. (15) becomes

$$\begin{array}{}\text{(17)}& \mathit{m}=EI\mathit{\kappa }\mathit{b}+GJ\mathit{\tau }\mathit{t}.\end{array}$$

The set of equations is simplified if some considerations based on the geometry of the structure and the forces and torques that are expected to act on it are taken into account. We start by assuming that $\mid \dot{\mathit{\alpha }}\mid \gg \mid {\mathit{\omega }}_{FE}\mid$ and writing ${\mathit{\omega }}_{DE}\approx -\dot{\mathit{\alpha }}\mathit{t}$. The left-hand side of Eq. (12) then becomes

$$\begin{array}{}\text{(18)}& \mathit{\rho }\left(I_{\mathrm{2}}{\ddot{\mathit{d}}}_{\mathrm{1}}\times {\mathit{d}}_{\mathrm{1}}+I_{\mathrm{1}}{\ddot{\mathit{d}}}_{\mathrm{2}}\times {\mathit{d}}_{\mathrm{2}}\right)=\mathrm{2}\mathit{\rho }I\ddot{\mathit{\alpha }}\mathit{t},\end{array}$$

where we used Eq. (9). Regarding the right-hand side, we use Eq. (4) and write $\mathit{f}=T\mathit{t}+{\mathit{f}}_{⟂}$ with ${\mathit{f}}_{⟂}\cdot \mathit{t}=\mathrm{0}$ to find

$$\begin{array}{}\text{(19)}& {\mathit{m}}^{\prime }+{\mathit{r}}^{\prime }\times \mathit{f}=GJ{\mathit{\tau }}^{\prime }\mathit{t}+\mathit{t}\times \left[{\mathit{f}}_{⟂}-GJ\mathit{\kappa }\mathit{\tau }\mathit{b}+EI{\left(\mathit{\kappa }\mathit{n}\right)}^{\prime }\right].\end{array}$$

Equation (12) then becomes

$$\begin{array}{}\text{(20)}& \mathrm{2}\mathit{\rho }I\ddot{\mathit{\alpha }}\mathit{t}=GJ{\mathit{\tau }}^{\prime }\mathit{t}+\mathit{t}\times \left[{\mathit{f}}_{⟂}-GJ\mathit{\kappa }\mathit{\tau }\mathit{b}+EI{\left(\mathit{\kappa }\mathit{n}\right)}^{\prime }\right].\end{array}$$

The components normal to t of Eq. (20) give f_⟂, and the total internal force then becomes

$$\begin{array}{}\text{(21)}& \begin{array}{rl}\mathit{f}& =\left(T-EI{\mathit{\kappa }}^{\mathrm{2}}\right)\mathit{t}+GJ\mathit{\kappa }\mathit{\tau }\mathit{b}-EI{\left(\mathit{\kappa }\mathit{n}\right)}^{\prime }\\ & =\left(T-EI{\mathit{\kappa }}^{\mathrm{2}}\right){\mathit{r}}^{\prime }+GJ\mathit{\tau }{\mathit{r}}^{\prime }\times {\mathit{r}}^{\prime \prime }-EI{\mathit{r}}^{\prime \prime \prime },\end{array}\end{array}$$

where we used Eq. (2). For the tension we assume Hooke's law:

$$\begin{array}{}\text{(22)}& T=EA\left(\left|{\mathit{r}}^{\prime }\right|-\mathrm{1}\right).\end{array}$$

This law is consistent with our assumption $\left(\left|{\mathit{r}}^{\prime }\right|-\mathrm{1}\right)\ll \mathrm{1}$, which was also used in Eq. (2). For large deformations, the model should be revisited.

After substituting these results in Eqs. (11) and (20), one finds

$$\begin{array}{}\text{(23)}& {\displaystyle }\begin{array}{rl}\mathit{\rho }A\ddot{\mathit{r}}& ={\left\{\left(T-EI{\mathit{\kappa }}^{\mathrm{2}}\right){\mathit{r}}^{\prime }+GJ\mathit{\tau }\left({\mathit{r}}^{\prime }\times {\mathit{r}}^{\prime \prime }\right)\right\}}^{\prime }\\ & -EI{\mathit{r}}^{\prime \prime \prime \prime }+\mathit{\rho }Ag{\mathit{k}}_E+{\mathit{F}}_{\mathrm{A}},\end{array}\text{(24)}& {\displaystyle }\mathrm{2}\mathit{\rho }I\ddot{\mathit{\alpha }}=GJ\left({\mathit{\alpha }}^{\prime \prime }+{\mathit{\gamma }}^{\prime }\right).\end{array}$$

The above set of equations describes the dynamics of the center line of the helix and its torsional deformation. Given appropriate initial and boundary conditions, they can be integrated numerically to find their evolution.

For convenience, we introduce the normalized variables

$$\begin{array}{}\text{(25)}& \stackrel{\mathrm{̃}}{t}\equiv \sqrt{{\displaystyle \frac{g}{L_{\mathrm{0}}}}}t,\stackrel{\mathrm{̃}}{s}\equiv {\displaystyle \frac{s}{L_{\mathrm{0}}}}\stackrel{\mathrm{̃}}{\mathit{r}}(\stackrel{\mathrm{̃}}{s},\mathit{\tau })\equiv {\displaystyle \frac{\mathit{r}}{L_{\mathrm{0}}}},{\stackrel{\mathrm{̃}}{\mathit{v}}}_{\mathrm{A}}\equiv {\displaystyle \frac{{\mathit{v}}_{\mathrm{A}}}{\sqrt{g{L}_{\mathrm{0}}}}},\end{array}$$

parameters

$$\begin{array}{}\text{(26)}& {\displaystyle }\mathit{\mu }\equiv {\displaystyle \frac{EI}{Mg{L}_{\mathrm{0}}^{\mathrm{2}}}},\mathit{\beta }\equiv {\displaystyle \frac{GJ}{Mg{L}_{\mathrm{0}}^{\mathrm{2}}}},\mathit{\nu }\equiv {\displaystyle \frac{{\mathit{\rho }}_{\mathrm{0}}{l}_a{L}_{\mathrm{0}}^{\mathrm{2}}}{\mathrm{2}M}}{C}_{\mathrm{D}},\text{(27)}& {\displaystyle }\mathit{\sigma }\equiv {\displaystyle \frac{EA}{Mg}},\mathit{\delta }\equiv {\displaystyle \frac{GJ}{\mathrm{2}\mathit{\rho }Ig{L}_{\mathrm{0}}}}={\displaystyle \frac{\mathit{\beta }\mathit{\sigma }}{\mathrm{2}\mathit{\mu }}},{\stackrel{\mathrm{̃}}{v}}_{\mathrm{w}\mathrm{0}}={\displaystyle \frac{v_{\mathrm{w}\mathrm{0}}}{\sqrt{g{L}_{\mathrm{0}}}}},\end{array}$$

and normalized forces

$$\begin{array}{}\text{(28)}& {\displaystyle }{\stackrel{\mathrm{̃}}{\mathit{F}}}_{\mathrm{A}}\left(s\right)\equiv {\displaystyle \frac{L_{\mathrm{0}}{\mathit{F}}_{\mathrm{A}}}{Mg}}=-\mathit{\nu }{\stackrel{\mathrm{̃}}{v}}_{\mathrm{A}}{\stackrel{\mathrm{̃}}{\mathit{v}}}_{\mathrm{A}},\text{(29)}& {\displaystyle }\stackrel{\mathrm{̃}}{T}\left(s\right)\equiv {\displaystyle \frac{T}{Mg}}\equiv \mathit{\sigma }\left(\left|{\stackrel{\mathrm{̃}}{\mathit{r}}}^{\prime }\right|-\mathrm{1}\right).\end{array}$$

We also changed the notation, and the primes now denote the derivative with respect to the normalized arc length $\stackrel{\mathrm{̃}}{s}$. The equations of motion become.

$$\begin{array}{}\text{(30)}& {\displaystyle \frac{{\partial }^{\mathrm{2}}\stackrel{\mathrm{̃}}{\mathit{r}}}{\partial {\stackrel{\mathrm{̃}}{t}}^{\mathrm{2}}}}={\left[{\stackrel{\mathrm{̃}}{\mathit{F}}}_{\mathrm{NL}}-{\stackrel{\mathrm{̃}}{\mathit{F}}}_{\mathrm{L}}\right]}^{\prime }+{\mathit{k}}_E+{\stackrel{\mathrm{̃}}{\mathit{F}}}_{\mathrm{A}},\text{(31)}& {\displaystyle \frac{{\partial }^{\mathrm{2}}\mathit{\alpha }}{\partial {\stackrel{\mathrm{̃}}{t}}^{\mathrm{2}}}}=\mathit{\delta }{\stackrel{\mathrm{̃}}{\mathit{\tau }}}^{\prime },\end{array}$$

with

$$\begin{array}{}\text{(32)}& {\displaystyle }{\stackrel{\mathrm{̃}}{\mathit{F}}}_{\mathrm{NL}}=\left(\stackrel{\mathrm{̃}}{T}-\mathit{\mu }{\stackrel{\mathrm{̃}}{\mathit{\kappa }}}^{\mathrm{2}}\right){\stackrel{\mathrm{̃}}{\mathit{r}}}^{\prime }+\mathit{\beta }\stackrel{\mathrm{̃}}{\mathit{\tau }}\left({\stackrel{\mathrm{̃}}{\mathit{r}}}^{\prime }\times {\stackrel{\mathrm{̃}}{\mathit{r}}}^{\prime \prime }\right),\text{(33)}& {\displaystyle }{\stackrel{\mathrm{̃}}{\mathit{F}}}_{\mathrm{L}}=\mathit{\mu }{\stackrel{\mathrm{̃}}{\mathit{r}}}^{\prime \prime \prime },\end{array}$$

and

$$\begin{array}{}\text{(34)}& \stackrel{\mathrm{̃}}{\mathit{\kappa }}={\left({\stackrel{\mathrm{̃}}{\mathit{r}}}^{\prime \prime }\cdot {\stackrel{\mathrm{̃}}{\mathit{r}}}^{\prime \prime }\right)}^{\mathrm{1}/\mathrm{2}},\stackrel{\mathrm{̃}}{\mathit{\gamma }}={\displaystyle \frac{{\stackrel{\mathrm{̃}}{\mathit{r}}}^{\prime }\cdot \left({\stackrel{\mathrm{̃}}{\mathit{r}}}^{\prime \prime }\times {\stackrel{\mathrm{̃}}{\mathit{r}}}^{\prime \prime \prime }\right)}{{\stackrel{\mathrm{̃}}{\mathit{r}}}^{\prime \prime }\cdot {\stackrel{\mathrm{̃}}{\mathit{r}}}^{\prime \prime }}},\stackrel{\mathrm{̃}}{\mathit{\tau }}\equiv \stackrel{\mathrm{̃}}{\mathit{\gamma }}+{\displaystyle \frac{\partial \mathit{\alpha }}{\partial \stackrel{\mathrm{̃}}{s}}}.\end{array}$$

Regarding the boundary conditions, we do not impose the exact conditions of a real rotary machine, i.e., coupling the dynamics of the structure and the rotor. Such a complete analysis is beyond the scope of this work that is focused on the helix. We assume that the position vector of the central line at $\stackrel{\mathrm{̃}}{s}=\mathrm{0}$ coincides with the origin of the inertial system 𝒮_E. At $\stackrel{\mathrm{̃}}{s}=\mathrm{1}$, we also impose a prescribed position. We then have

$$\begin{array}{}\text{(35)}& {\displaystyle }\stackrel{\mathrm{̃}}{\mathit{r}}(\mathrm{0},\stackrel{\mathrm{̃}}{t})=\mathrm{0},\text{(36)}& {\displaystyle }\stackrel{\mathrm{̃}}{\mathit{r}}(\mathrm{1},\stackrel{\mathrm{̃}}{t})={\stackrel{\mathrm{̃}}{\mathit{r}}}_N\left(\stackrel{\mathrm{̃}}{t}\right).\end{array}$$

Additionally, we impose zero curvature at both tips,

$$\begin{array}{}\text{(37)}& {\stackrel{\mathrm{̃}}{\mathit{r}}}^{\prime \prime }(\mathrm{0},\stackrel{\mathrm{̃}}{t})=\mathrm{0},{\stackrel{\mathrm{̃}}{\mathit{r}}}^{\prime \prime }(\mathrm{1},\stackrel{\mathrm{̃}}{t})=\mathrm{0},\end{array}$$

and known torsional torques at $\stackrel{\mathrm{̃}}{s}=\mathrm{0}$ and $\stackrel{\mathrm{̃}}{s}=\mathrm{1}$:

$$\begin{array}{}\text{(38)}& {\displaystyle }\stackrel{\mathrm{̃}}{\mathit{\tau }}(\mathrm{0},\stackrel{\mathrm{̃}}{t})=& {\displaystyle \frac{\partial \mathit{\alpha }}{\partial \stackrel{\mathrm{̃}}{s}}}{\mid }_{\stackrel{\mathrm{̃}}{s}=\mathrm{0},\stackrel{\mathrm{̃}}{t}}=m_{\mathrm{0}}\left(\stackrel{\mathrm{̃}}{t}\right)\equiv {\displaystyle \frac{L{M}_{\mathrm{0}}\left(t\right)}{GJ}},\text{(39)}& {\displaystyle }\stackrel{\mathrm{̃}}{\mathit{\tau }}(\mathrm{1},\stackrel{\mathrm{̃}}{t})=& {\displaystyle \frac{\partial \mathit{\alpha }}{\partial \stackrel{\mathrm{̃}}{s}}}{\mid }_{\stackrel{\mathrm{̃}}{s}=\mathrm{1},\stackrel{\mathrm{̃}}{t}}=m_N\left(\stackrel{\mathrm{̃}}{t}\right)\equiv {\displaystyle \frac{L{M}_N\left(t\right)}{GJ}},\end{array}$$

where M₀(t) and M_N(t) are the external torques imposed at the ends of the structure. Equations (38)–(39) involve the use of Eq. (3) and Eq. (37) in Eq. (8). As explained below, the initial conditions depend on the type of analysis.

Since our work does not couple the helix with a rotor, the torque at $\stackrel{\mathrm{̃}}{s}=\mathrm{1}$ is a given function $m_N\left(\stackrel{\mathrm{̃}}{t}\right)$. Regarding the torque at $\stackrel{\mathrm{̃}}{s}=\mathrm{0}$, we implement the controller

$$\begin{array}{}\text{(40)}& {\dot{m}}_{\mathrm{0}}\left(\stackrel{\mathrm{̃}}{t}\right)={\dot{m}}_N\left(\stackrel{\mathrm{̃}}{t}\right)+k_{\mathrm{1}}\left[{\dot{\mathit{\alpha }}}_{\mathrm{1}}\left(\stackrel{\mathrm{̃}}{t}\right)-{\dot{\mathit{\alpha }}}^{*}\right]+k_{\mathrm{2}}{\ddot{\mathit{\alpha }}}_{\mathrm{1}}\left(\stackrel{\mathrm{̃}}{t}\right),\end{array}$$

where k₁ and k₂ are constant, and ${\ddot{\mathit{\alpha }}}_{\mathrm{1}}$ is given by the first equation in Eq. (56). Such a proportional-derivative controller adjusts the torque exerted at the ground to reach a (constant) target angular velocity ${\dot{\mathit{\alpha }}}^{*}$ at the first node.

3 A finite-element method

Equations (30)–(31) were solved numerically by using a finite-element method. We introduce N elements and N+1 nodes, ${\stackrel{\mathrm{̃}}{\mathit{r}}}_j\left(\stackrel{\mathrm{̃}}{t}\right)\equiv \stackrel{\mathrm{̃}}{\mathit{r}}({\stackrel{\mathrm{̃}}{s}}_j,\stackrel{\mathrm{̃}}{t})$ and ${\mathit{\alpha }}_j\left(\stackrel{\mathrm{̃}}{t}\right)=\mathit{\alpha }({\stackrel{\mathrm{̃}}{s}}_j,\stackrel{\mathrm{̃}}{t})$ with ${\stackrel{\mathrm{̃}}{s}}_j=j/N$, and $j=\mathrm{0},\mathrm{1},\mathrm{\dots },N$. For each element between nodes ${\stackrel{\mathrm{̃}}{\mathit{r}}}_j$ and ${\stackrel{\mathrm{̃}}{\mathit{r}}}_{j+\mathrm{1}}$, the variables $\stackrel{\mathrm{̃}}{\mathit{r}}(\stackrel{\mathrm{̃}}{s},\stackrel{\mathrm{̃}}{t})$ and $\mathit{\alpha }(\stackrel{\mathrm{̃}}{s},\stackrel{\mathrm{̃}}{t})$ are written as (Buckham, 2003)

$$\begin{array}{}\text{(41)}& {\displaystyle }\stackrel{\mathrm{̃}}{\mathit{r}}(\stackrel{\mathrm{̃}}{s},\stackrel{\mathrm{̃}}{t})=& {\displaystyle }\sum _{k=\mathrm{0}}^{\mathrm{1}}{\stackrel{\mathrm{̃}}{\mathit{r}}}_{j+k}\left(\stackrel{\mathrm{̃}}{t}\right)N_k^j\left(\stackrel{\mathrm{̃}}{s}\right)+{\stackrel{\mathrm{̃}}{\mathit{r}}}_{j+k}^{\prime \prime }\left(\stackrel{\mathrm{̃}}{t}\right)H_k^j\left(\stackrel{\mathrm{̃}}{s}\right),\text{(42)}& {\displaystyle }\mathit{\alpha }(\stackrel{\mathrm{̃}}{s},\stackrel{\mathrm{̃}}{t})=& {\displaystyle }\sum _{k=\mathrm{0}}^{\mathrm{1}}{\mathit{\alpha }}_{j+k}\left(\stackrel{\mathrm{̃}}{t}\right)N_k^j\left(\stackrel{\mathrm{̃}}{s}\right)+{\mathit{\alpha }}_{j+k}^{\prime \prime }\left(\stackrel{\mathrm{̃}}{t}\right)H_k^j\left(\stackrel{\mathrm{̃}}{s}\right),\end{array}$$

with ${\stackrel{\mathrm{̃}}{s}}_j\le \stackrel{\mathrm{̃}}{s}\le {\stackrel{\mathrm{̃}}{s}}_{j+\mathrm{1}}$. The shape functions are the twisted-spline polynomials:

$$\begin{array}{}\text{(43)}& {\displaystyle }{N}_{\mathrm{0}}^j\left(\stackrel{\mathrm{̃}}{s}\right)=\mathrm{1}-{\mathit{\xi }}_j,H_{\mathrm{0}}^j\left(\stackrel{\mathrm{̃}}{s}\right)=-{\displaystyle \frac{\mathrm{1}}{\mathrm{6}}}\left({\mathit{\xi }}_j^{\mathrm{3}}-\mathrm{3}{\mathit{\xi }}_j^{\mathrm{2}}+\mathrm{2}{\mathit{\xi }}_j\right),\text{(44)}& {\displaystyle }{N}_{\mathrm{1}}^j\left(\stackrel{\mathrm{̃}}{s}\right)={\mathit{\xi }}_j,H_{\mathrm{1}}^j\left(\stackrel{\mathrm{̃}}{s}\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{6}}}\left({\mathit{\xi }}_j^{\mathrm{3}}-{\mathit{\xi }}_j\right),\end{array}$$

where

$$\begin{array}{}\text{(45)}& {\mathit{\xi }}_j={\displaystyle \frac{\stackrel{\mathrm{̃}}{s}-{\stackrel{\mathrm{̃}}{s}}_j}{{\stackrel{\mathrm{̃}}{s}}_{j+\mathrm{1}}-{\stackrel{\mathrm{̃}}{s}}_j}}\end{array}$$

is the local coordinate. One readily verifies that the coefficients ${\stackrel{\mathrm{̃}}{\mathit{r}}}_j\left(\stackrel{\mathrm{̃}}{t}\right)$ and ${\mathit{\alpha }}_j\left(\stackrel{\mathrm{̃}}{t}\right)$ in Eqs. (41) and (42) coincide with the displacements and the rotation angles at the nodes and ${\stackrel{\mathrm{̃}}{\mathit{r}}}_j^{\prime \prime }\left(\stackrel{\mathrm{̃}}{t}\right)$ and ${\mathit{\alpha }}_j^{\prime \prime }\left(\stackrel{\mathrm{̃}}{t}\right)$ with their second derivatives. This choice for the shape function thus gives continuity for $\stackrel{\mathrm{̃}}{\mathit{r}}(\stackrel{\mathrm{̃}}{s},\stackrel{\mathrm{̃}}{t})$, ${\stackrel{\mathrm{̃}}{\mathit{r}}}^{\prime \prime }(\stackrel{\mathrm{̃}}{s},\stackrel{\mathrm{̃}}{t})$, $\mathit{\alpha }(\stackrel{\mathrm{̃}}{s},\stackrel{\mathrm{̃}}{t})$, and ${\mathit{\alpha }}^{\prime \prime }(\stackrel{\mathrm{̃}}{s},\stackrel{\mathrm{̃}}{t})$. As explained in Appendix A, by assuming that the first derivative is continuous across the elements, we find relations between ${\stackrel{\mathrm{̃}}{\mathit{r}}}_j^{\prime \prime }\left(\stackrel{\mathrm{̃}}{t}\right)$ and ${\stackrel{\mathrm{̃}}{\mathit{r}}}_j\left(\stackrel{\mathrm{̃}}{t}\right)$ and between ${\mathit{\alpha }}_j^{\prime \prime }\left(\stackrel{\mathrm{̃}}{t}\right)$ and ${\mathit{\alpha }}_j\left(\stackrel{\mathrm{̃}}{t}\right)$. Equations (41)–(42) for all the elements, together with relation Eqs. (A1)–(A4) and the boundary condition Eqs. (35)–(36), allow us to write the variables as

$$\begin{array}{}\text{(46)}& {\displaystyle }\stackrel{\mathrm{̃}}{\mathit{r}}(\stackrel{\mathrm{̃}}{s},\stackrel{\mathrm{̃}}{t})=\sum _{j=\mathrm{1}}^{N-\mathrm{1}}{\stackrel{\mathrm{̃}}{\mathit{r}}}_j\left(\stackrel{\mathrm{̃}}{t}\right){\mathit{\varphi }}_j\left(\stackrel{\mathrm{̃}}{s}\right)+{\stackrel{\mathrm{̃}}{\mathit{r}}}_N\left(\stackrel{\mathrm{̃}}{t}\right){\mathit{\varphi }}_N\left(\stackrel{\mathrm{̃}}{s}\right),\text{(47)}& {\displaystyle }\mathit{\alpha }(\stackrel{\mathrm{̃}}{s},\stackrel{\mathrm{̃}}{t})=\sum _{j=\mathrm{0}}^N{\mathit{\alpha }}_j\left(\stackrel{\mathrm{̃}}{t}\right){\mathit{\phi }}_j\left(\stackrel{\mathrm{̃}}{s}\right)+m(\stackrel{\mathrm{̃}}{s},\stackrel{\mathrm{̃}}{t}),\end{array}$$

where ϕ_i and φ_i are 𝒞² continuous functions, and the function $m(\stackrel{\mathrm{̃}}{s},\stackrel{\mathrm{̃}}{t})$ just involves the boundary conditions (m₀ and m_N) and the shape functions. One readily checks that ϕ_i(0)=0 for i=1…N, ϕ_N(1)=1, and ϕ_i(1)=0 for $i=\mathrm{0}\mathrm{\dots }N-\mathrm{1}$.

Following a Galerkin method, we write Eqs. (30)–(31) in the weak form. We multiply Eq. (30) by ϕ_i with $i=\mathrm{1}\mathrm{\dots }N-\mathrm{1}$ and Eq. (31) by φ_i with i=0…N and integrate between $\stackrel{\mathrm{̃}}{s}=\mathrm{0}$ and $\stackrel{\mathrm{̃}}{s}=\mathrm{1}$ to find

$$\begin{array}{}\text{(48)}& {\displaystyle }\begin{array}{rl}\sum _{j=\mathrm{1}}^{N-\mathrm{1}}{P}_{ij}^{\left(\mathrm{0}\right)}{\ddot{\stackrel{\mathrm{̃}}{\mathit{r}}}}_j& ={\mathit{a}}_{iN}+{\mathit{p}}_i^{\left(\mathrm{0}\right)}\\ & +\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}\left[{\left({\stackrel{\mathrm{̃}}{\mathit{F}}}_{\mathrm{NL}}-{\stackrel{\mathrm{̃}}{\mathit{F}}}_L\right)}^{\prime }+{\stackrel{\mathrm{̃}}{\mathit{F}}}_{\mathrm{A}}\right]{\mathit{\varphi }}_i\mathrm{d}\stackrel{\mathrm{̃}}{s},\\ & i=\mathrm{1}\mathrm{\dots }N-\mathrm{1},\end{array}\text{(49)}& {\displaystyle }\sum _{j=\mathrm{0}}^N{Q}_{ij}^{\left(\mathrm{0}\right)}{\ddot{\mathit{\alpha }}}_j={\ddot{m}}_i^{\mathrm{0}}+\mathit{\delta }\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}{\stackrel{\mathrm{̃}}{\mathit{\tau }}}^{\prime }{\mathit{\phi }}_i\left(\stackrel{\mathrm{̃}}{s}\right)\mathrm{d}\stackrel{\mathrm{̃}}{s},i=\mathrm{0}\mathrm{\dots }N,\end{array}$$

where we introduce the matrices and the vectors

$$\begin{array}{}\text{(50)}& {\displaystyle }{P}_{ij}^{\left(\mathrm{0}\right)}=\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}{\mathit{\varphi }}_j\left(\stackrel{\mathrm{̃}}{s}\right){\mathit{\varphi }}_i\left(\stackrel{\mathrm{̃}}{s}\right)\mathrm{d}\stackrel{\mathrm{̃}}{s},\text{(51)}& {\displaystyle }{\mathit{a}}_{iN}=-{\ddot{\stackrel{\mathrm{̃}}{\mathit{r}}}}_N\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}{\mathit{\varphi }}_i\left(\stackrel{\mathrm{̃}}{s}\right){\mathit{\varphi }}_N\left(\stackrel{\mathrm{̃}}{s}\right)\mathrm{d}s,\text{(52)}& {\displaystyle }{\mathit{p}}_i^{\left(\mathrm{0}\right)}={\mathit{k}}_E\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}{\mathit{\varphi }}_i\left(\stackrel{\mathrm{̃}}{s}\right)\mathrm{d}\stackrel{\mathrm{̃}}{s},\text{(53)}& {\displaystyle }{Q}_{ij}^{\left(\mathrm{0}\right)}=\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}{\mathit{\phi }}_j\left(\stackrel{\mathrm{̃}}{s}\right){\mathit{\phi }}_i\left(\stackrel{\mathrm{̃}}{s}\right)\mathrm{d}\stackrel{\mathrm{̃}}{s},\text{(54)}& {\displaystyle }{\ddot{m}}_i^{\mathrm{0}}=-\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}\ddot{m}(\stackrel{\mathrm{̃}}{s},\stackrel{\mathrm{̃}}{t}){\mathit{\phi }}_i\mathrm{d}\stackrel{\mathrm{̃}}{s}.\end{array}$$

Integrations by parts of Eqs. (48) and (49) yield

$$\begin{array}{}\text{(55)}& {\displaystyle }\begin{array}{rl}& \sum _{j=\mathrm{1}}^{N-\mathrm{1}}{P}_{ij}^{\left(\mathrm{0}\right)}{\ddot{\stackrel{\mathrm{̃}}{\mathit{r}}}}_j={\mathit{a}}_{iN}+{\mathit{p}}_i^{\left(\mathrm{0}\right)}+\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}{\stackrel{\mathrm{̃}}{\mathit{F}}}_{\mathrm{A}}{\mathit{\varphi }}_i\mathrm{d}\stackrel{\mathrm{̃}}{s}\\ & -\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}\left({\stackrel{\mathrm{̃}}{\mathit{F}}}_{\mathrm{NL}}-{\stackrel{\mathrm{̃}}{\mathit{F}}}_L\right){\mathit{\varphi }}_i^{\prime }\mathrm{d}\stackrel{\mathrm{̃}}{s}\\ & ={\mathit{a}}_{iN}+{\mathit{p}}_i^{\left(\mathrm{0}\right)}-\mathit{\mu }\sum _{j=\mathrm{1}}^N{P}_{ij}^{\left(\mathrm{4}\right)}{\stackrel{\mathrm{̃}}{\mathit{r}}}_j\\ & +\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}\left({\stackrel{\mathrm{̃}}{\mathit{F}}}_{\mathrm{A}}{\mathit{\varphi }}_i-{\stackrel{\mathrm{̃}}{\mathit{F}}}_{\mathrm{NL}}{\mathit{\varphi }}_i^{\prime }\right)\mathrm{d}\stackrel{\mathrm{̃}}{s},i=\mathrm{1}\mathrm{\dots }N-\mathrm{1},\end{array}\text{(56)}& {\displaystyle }\begin{array}{rl}& \sum _{j=\mathrm{0}}^N{Q}_{ij}^{\left(\mathrm{0}\right)}{\ddot{\mathit{\alpha }}}_j={\ddot{m}}_i^{\mathrm{0}}+\mathit{\delta }\left\{{\left[\stackrel{\mathrm{̃}}{\mathit{\tau }}{\mathit{\phi }}_i\right]}_{\mathrm{0}}^{\mathrm{1}}-\sum _{j=\mathrm{0}}^N{Q}_{ij}^{\mathrm{1}}{\mathit{\alpha }}_j\right.\\ & \left.-\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}\left(\stackrel{\mathrm{̃}}{\mathit{\gamma }}+{\displaystyle \frac{\partial m}{\partial \stackrel{\mathrm{̃}}{s}}}\right){\mathit{\phi }}_i^{\prime }\left(\stackrel{\mathrm{̃}}{s}\right)\mathrm{d}\stackrel{\mathrm{̃}}{s}\right\},i=\mathrm{0}\mathrm{\dots }N,\end{array}\end{array}$$

where we use ${\mathit{\varphi }}_i\left(\mathrm{0}\right)={\mathit{\varphi }}_i\left(\mathrm{1}\right)=\mathrm{0}$ for $i=\mathrm{1}\mathrm{\dots }N-\mathrm{1}$, Eq. (37), and define

$$\begin{array}{}\text{(57)}& {\displaystyle }{Q}_{ij}^{\left(\mathrm{1}\right)}=\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}{\mathit{\phi }}_j^{\prime }\left(\stackrel{\mathrm{̃}}{s}\right){\mathit{\phi }}_i^{\prime }\left(\stackrel{\mathrm{̃}}{s}\right)\mathrm{d}\stackrel{\mathrm{̃}}{s},\text{(58)}& {\displaystyle }{P}_{ij}^{\left(\mathrm{4}\right)}=\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}{\mathit{\varphi }}_j^{\prime \prime }\left(\stackrel{\mathrm{̃}}{s}\right){\mathit{\varphi }}_i^{\prime \prime }\left(\stackrel{\mathrm{̃}}{s}\right)\mathrm{d}\stackrel{\mathrm{̃}}{s}.\end{array}$$

Since i runs from 1 to N−1 in Eq. (55) and from 0 to N in Eq. (56), there are a total of $\mathrm{3}(N-\mathrm{1})+N+\mathrm{1}=\mathrm{4}N-\mathrm{2}$ equations for the 4N−2 unknowns, which are ${\mathit{r}}_j\left(\stackrel{\mathrm{̃}}{t}\right)$ with $j=\mathrm{1}\mathrm{\dots }N-\mathrm{1}$ and ${\mathit{\alpha }}_j\left(\stackrel{\mathrm{̃}}{t}\right)$ with j=0…N.

After constructing the vector

$$\begin{array}{}\text{(59)}& {\mathit{q}}_{\mathrm{s}}=\left[{\stackrel{\mathrm{̃}}{\mathit{r}}}_{\mathrm{1}}{\stackrel{\mathrm{̃}}{\mathit{r}}}_{\mathrm{2}}\mathrm{\cdots }{\stackrel{\mathrm{̃}}{\mathit{r}}}_{N-\mathrm{1}}{\mathit{\alpha }}_{\mathrm{0}}\mathrm{\cdots }{\mathit{\alpha }}_N\right]\end{array}$$

and the state vector

$$\begin{array}{}\text{(60)}& {\mathit{x}}_{\mathrm{s}}=\left[{\mathit{q}}_{\mathrm{s}}{\dot{\mathit{q}}}_{\mathrm{s}}{m}_{\mathrm{0}}\right],\end{array}$$

we find the following set of first-order ordinary differential equations:

$$\begin{array}{}\text{(61)}& \stackrel{\mathrm{‾}}{\mathbf{M}}{\dot{\mathit{x}}}_s=\mathit{g}\left({\mathit{x}}_{\mathrm{s}}\right),\mathrm{or},{\dot{\mathit{x}}}_s=\mathit{f}\left({\mathit{x}}_{\mathrm{s}}\right)\equiv {\stackrel{\mathrm{‾}}{\mathbf{M}}}^{-\mathrm{1}}\mathit{g}\left({\mathit{x}}_{\mathrm{s}}\right).\end{array}$$

The explicit form of matrix $\stackrel{\mathrm{‾}}{\mathbf{M}}$ and the vector field g can be found from Eqs. (55), (56), and (40).

Several tests were carried out to verify the correct implementation of Eq. (61). For instance, to test the correct implementation of functions ϕ_j and φ_j and matrices B_r, C_r, etc., some analytical functions satisfying the boundary conditions were proposed for $\stackrel{\mathrm{̃}}{\mathit{r}}$ and $\stackrel{\mathrm{̃}}{\mathit{\alpha }}$. The derivatives $\partial \stackrel{\mathrm{̃}}{\mathit{r}}/\partial \stackrel{\mathrm{̃}}{s}$, ${\partial }^{\mathrm{2}}\stackrel{\mathrm{̃}}{\mathit{r}}/\partial {\stackrel{\mathrm{̃}}{s}}^{\mathrm{2}}$, ${\partial }^{\mathrm{3}}\stackrel{\mathrm{̃}}{\mathit{r}}/\partial {\stackrel{\mathrm{̃}}{s}}^{\mathrm{3}}$, $\partial \mathit{\alpha }/\partial s$, and ${\partial }^{\mathrm{2}}\mathit{\alpha }/\partial s^{\mathrm{2}}$ were computed analytically and also by using the finite-element approach. It was verified that both results converge as the number of elements increases. The analysis of Sect. 5.1, where a stationary solution is found with an alternative numerical method, constitutes another test of the correct implementation of the finite-element method.

4 Experimental characterization of the structure

Three laboratory tests were carried out to find the values of EA, EI, and GJ for the tensegrity structure of SomeAWE Labs. Detailed information, drawings, and videos are in Beaupoil (2024). As shown in the three panels of Fig. 2, a segment of the helix with length equal to L₀=6.76 m was clamped on a wall. As explained below, different types of loads were applied to determine the axial, bending, and torsional stiffness. For all the tests, a tension force T_e was applied to the structure, and its values were taken of the order of around 200 N, which is the nominal tension during the regular operation of the RAWE machine. In cases (a) and (c) in Fig. 2, i.e., the axial and torsion tests, such axial force was applied by using a tension belt. In the bending experiment (case b), where there is a lateral displacement, a pulley and a ballast mass were used, and the position of the pulley was adjusted to apply the force normal to the displacement. Figure 2d displays the helix configuration during the axial test and the tension belt.

Figure 2 Panels (a–c) show the sketches of the axial, bending, and torsion tests. Panels (d) and (e) display the helix in the axial test and the motorized eccentric arm explained in Sect. 5.3.

As shown in the top panel of Fig. 2, the EA product was found by applying an axial force T_e and measuring the resulting length of the structure L. The axial force was varied by just increasing the mass of the ballast. For a linear material, the relation between the axial force and the length of the structure is

$$\begin{array}{}\text{(62)}& T_e=EA\times {\displaystyle \frac{L-L_{\mathrm{0}}}{L_{\mathrm{0}}}}.\end{array}$$

As expected, the experimental values of the strain are proportional to the applied force (see Fig. 3). The slope obtained with a linear regression is $EA=\mathrm{3.64}\times {\mathrm{10}}^{\mathrm{4}}\mathrm{N}$.

Figure 3 Experimental results of the axial test.

A classical bending test was carried out to determine the EI product of the helix. As shown in the middle panel of Fig. 2, a normal force of value F_N was applied at the free end of the structure, and the vertical displacement y_e was measured. Since the theoretical relation between both magnitudes is

$$\begin{array}{}\text{(63)}& F_N=EI\times {\displaystyle \frac{\mathrm{3}{y}_e}{L_{\mathrm{0}}^{\mathrm{3}}}},\end{array}$$

measuring the pairs (F_N,y_e) and making a linear regression to the experimental data give the EI product. However, the bending stiffness of the tensegrity structure depends on the internal tension force. For this reason, the experiment was carried out for three different ballast masses, which yielded the tension forces T_e= 90, 181, and 259 N. As shown in Fig. 4, for each value of the tension, the helix has a different bending stiffness (slope of the linear regression in the F_N versus $\mathrm{3}{y}_e/L_{\mathrm{0}}^{\mathrm{3}}$ plane). The inset in Fig. 4 shows that the relation between the EI product and the tension is also linear and follows the law

$$\begin{array}{}\text{(64)}& EI\left(N{\mathrm{m}}^{\mathrm{2}}\right)=\mathrm{52.4}\times T_e\left(N\right)+\mathrm{2290}.\end{array}$$

As expected, the higher the tension, the higher the bending stiffness of the tensegrity structure.

Figure 4 Experimental results of the bending test.

Difficulties were found for the experimental determination of the torsional stiffness. As shown in the bottom panel of Fig. 2, the helical structure was tightened with the ballast mass to induce a tension force T_e0 for zero torsion angle (α=0). A torsion torque M_t was then applied at the free end of the structure until the torsion angle took a certain value α. The tension force T_e was then measured with the hanging scale, and it was found that it was higher, T_e>T_e0. Measurements of M_t, α, and T_e were taken by first increasing M_t sequentially and with an increment of the torsion angle of around Δα=10°. After reaching a maximum torsion angle of about 50°, M_t was lowered to produce decrements in the torsion angle of $\mathrm{\Delta }\mathit{\alpha }=-\mathrm{10}\mathit{°}$ until α=0 was reached. Four complete cycles raising and lowering the torsional torque were completed, and measurements for the three variables (M_t, α, and T_e) were taken. The maximum value of M_t in the experiment was around 25 N m, which is consistent with the torque measured during the real operation of the RAWE machine.

In principle, a simple model for the relation between the torsion torque and angle is

$$\begin{array}{}\text{(65)}& M_{\mathrm{t}}=GJ\times {\displaystyle \frac{\mathit{\alpha }}{L_{\mathrm{0}}}}.\end{array}$$

However, the experiment revealed that the behavior of the helix is not so simple. As seen in Fig. 5, the ratio $M_t{L}_{\mathrm{0}}/\mathit{\alpha }$, which may coincide with GJ according to Eq. (65), versus the tension force T_e present hysteresis. For the three tests covering a tension range between 180 and 240 N, the repeatability of the measurements is satisfactory, but it is clear that the ratio $M_t{L}_{\mathrm{0}}/\mathit{\alpha }$ cannot be parametrized by the tension force T_e alone. In a fourth test, the torsional behavior of the structure was explored at an even larger tension value, and the amplitude of the hysteresis cycle was similar. These results indicate that the value of the torsional rigidity depends not only on the instantaneous value of the torsion angle and internal tension but also on the history of the helix. Such a feature may be due to the peculiar construction of the tensegrity helix, which involves bars under compression, tethers, and knots. The unavoidable friction in the experimental setup may also play a role. For this reason, we decided in later analysis to set a nominal value for GJ of 140 N m², which is consistent with the experimental results of Fig. 5 for T_e=200 N.

Figure 5 Experimental results of the torsion test.

5 Dynamics of the helix: numerical and experimental results

This section studies three different dynamic scenarios of the helix. For all of them we used the parameters shown in Table 1.

Table 1 Common parameters used in the analysis.

5.1 Stationary solutions with fixed ends

We consider a helix with the properties of Tables 1 and 2, which correspond to the characteristics of the structure of the RAWE machine of SomeAwe. Since its nominal tension during regular operation is around 200 N, we set the elongation of the initial condition as ${\mathit{ϵ}}_{\mathrm{0}}=\mathrm{200}N/EA$ and compute EI by using Eq. (64) with T_e=200 N.

Table 2 Parameters used in the analysis of solutions with fixed ends.

The analysis focuses on stationary solutions of Eq. (61); i.e., all the variables are time independent, and the time derivatives of $\stackrel{\mathrm{̃}}{r}$ and $\stackrel{\mathrm{̃}}{\mathit{\alpha }}$ vanish. The boundary conditions are given by Eqs. (35)–(36), and

$$\begin{array}{}\text{(66)}& {\stackrel{\mathrm{̃}}{\mathit{r}}}_N\left(\stackrel{\mathrm{̃}}{t}\right)\equiv -(\mathrm{1}+{\mathit{ϵ}}_{\mathrm{0}})\left(\cos {\mathrm{\Gamma }}_{\mathrm{0}}{\mathit{i}}_E+\sin {\mathrm{\Gamma }}_{\mathrm{0}}{\mathit{k}}_E\right),\end{array}$$

where Γ₀ is the elevation angle of the helix, ϵ₀ is a parameter that determines its elongation, and i_E and k_E are unit vectors of frame 𝒮_E. We restrict the analysis to solutions of the type $\stackrel{\mathrm{̃}}{\mathit{r}}\left(\stackrel{\mathrm{̃}}{s}\right)=x\left(\stackrel{\mathrm{̃}}{s}\right){\mathit{i}}_E+z\left(\stackrel{\mathrm{̃}}{s}\right){\mathit{k}}_E$ and $\mathit{\alpha }(\stackrel{\mathrm{̃}}{s},\stackrel{\mathrm{̃}}{t})=\mathrm{0}$ with $k_{\mathrm{1}}=k_{\mathrm{2}}=m_{\mathrm{0}}=m_N=\mathrm{0}$. Since the proposed solution is contained in the x_E−z_E plane and no torque is applied at the ends of the helix, the torsion vanishes ($\mathit{\tau }\left(\stackrel{\mathrm{̃}}{s}\right)=\mathrm{0}$), and Eq. (31) is automatically fulfilled. Regarding Eq. (30), and after using ($\ddot{\stackrel{\mathrm{̃}}{r}}=\dot{\stackrel{\mathrm{̃}}{r}}=\mathrm{0}$), it becomes the following set of ordinary differential equations:

$$\begin{array}{}\text{(67)}& {\displaystyle \frac{d}{\mathrm{d}\stackrel{\mathrm{̃}}{s}}}\left(\begin{array}{c}x\\ z\\ x_s\\ z_s\\ x_{ss}\\ z_{ss}\\ x_{sss}\\ z_{sss}\end{array}\right)=\left(\begin{array}{c}{x}_s\\ z_s\\ x_{ss}\\ z_{ss}\\ x_{sss}\\ z_{sss}\\ \left[f_s{x}_s+\left(\stackrel{\mathrm{̃}}{T}-\mathit{\mu }{\mathit{\kappa }}^{\mathrm{2}}\right)x_{ss}-\mathit{\nu }{\stackrel{\mathrm{̃}}{v}}_{\mathrm{w}\mathrm{0}}^{\mathrm{2}}\right]/\mathit{\mu }\\ \left[f_s{z}_s+\left(\stackrel{\mathrm{̃}}{T}-\mathit{\mu }{\mathit{\kappa }}^{\mathrm{2}}\right)z_{ss}+\mathrm{1}\right]/\mathit{\mu }\end{array}\right),\end{array}$$

with ${\mathit{\kappa }}^{\mathrm{2}}=x_{ss}^{\mathrm{2}}+z_{ss}^{\mathrm{2}}$, $\stackrel{\mathrm{̃}}{T}=\mathit{\sigma }\left(\sqrt{x_s^{\mathrm{2}}+z_s^{\mathrm{2}}}-\mathrm{1}\right)$, and

$$\begin{array}{}\text{(68)}& f_s={\displaystyle \frac{\mathit{\sigma }}{\sqrt{x_s^{\mathrm{2}}+z_s^{\mathrm{2}}}}}\left(x_s{x}_{ss}+z_s{z}_{ss}\right)-\mathrm{2}\mathit{\mu }\left(x_{ss}{x}_{sss}+z_{ss}{z}_{sss}\right).\end{array}$$

The solution of System (67) compatible with the boundary conditions (35)–(36) was found by using a shooting method. The equations were integrated numerically with a Runge–Kutta method, and the state at s=0 was given by $\left(\mathrm{0},\mathrm{0},x_{s\mathrm{0}},z_{s\mathrm{0}},\mathrm{0},\mathrm{0},x_{sss\mathrm{0}},z_{sss\mathrm{0}}\right)$. Constants $x_{s\mathrm{0}},z_{s\mathrm{0}},x_{sss\mathrm{0}}$ and z_sss0 were varied with a shooting method until the conditions $x\left(\mathrm{1}\right)=-\left(\mathrm{1}+{\mathit{ϵ}}_{\mathrm{0}}\right)\cos {\mathrm{\Gamma }}_{\mathrm{0}}$, $z\left(\mathrm{1}\right)=-\left(\mathrm{1}+{\mathit{ϵ}}_{\mathrm{0}}\right)\sin {\mathrm{\Gamma }}_{\mathrm{0}}$, and $x_{ss}\left(\mathrm{1}\right)=z_{ss}\left(\mathrm{1}\right)=\mathrm{0}$ were satisfied. Once the solution was found, we constructed from it the state vector x_s of Eq. (60). Its substitution on the right-hand side of Eq. (61) revealed that it is indeed a stationary solution because we found $\mid \mathit{g}\left({\mathit{x}}_{\mathrm{s}}\right)\mid \approx \mathrm{0}$. The larger the number of finite elements, the lower the value of $\mid \mathit{g}\left({\mathit{x}}_{\mathrm{s}}^{*}\right)\mid$. This result constitutes an important test about the correct implementation of the finite-element method because x_s was found with an alternative numerical algorithm (an integration of the equations with a Runge–Kutta method).

The linear stability of the stationary solution was then investigated in the framework of Eq. (61) and setting N=2. Using the vector x_s obtained with the shooting algorithm as an initial guess, we found the state vector satisfying the condition $\mathit{f}\left({\mathit{x}}_{\mathrm{s}}^{*}\right)=\mathrm{0}$ in Eq. (61) with a Newton–Raphson method. Afterwards, the Jacobian matrix of flow f at ${\mathit{x}}_{\mathrm{s}}^{*}$ and its eigenvalues were computed. As shown in Table 3, the stationary solution is stable because the eigenvalues are pure imaginary numbers or complex numbers with a negative real part. There are two longitudinal modes where the dominant components of the eigenvectors are ${\dot{x}}_{\mathrm{1}}$ and ${\dot{z}}_{\mathrm{1}}$. The first of them has a short period (fast mode), which can be approximated by $T_{\mathrm{Long}}\approx \mathrm{2}{L}_{\mathrm{0}}/\sqrt{E/\mathit{\rho }}$ that is valid for thin and elastic rods of Young's modulus E and density ρ. In our case it reads $T_{\mathrm{Long}}\approx \mathrm{2}\sqrt{M{L}_{\mathrm{0}}/EA}\approx \mathrm{0.03}$ s and matches with the first mode in Table 3. There is also a lateral mode dominated by displacements contained in the x_E−z_E plane. The period of this mode is well described by the one of a string of linear density ρ_L, length L₀, tension T, and fixed tips. The nth harmonic vibration has wavelength ${\mathit{\lambda }}_n=\mathrm{2}{L}_{\mathrm{0}}/n$ and velocity $\sqrt{T/{\mathit{\rho }}_L}$. The period for n=2 in our case reads $T_{\mathrm{Lat}}\approx \sqrt{M{L}_{\mathrm{0}}/EA{\mathit{ϵ}}_{\mathrm{0}}}\approx \mathrm{0.20}$ s, where the tension is of the order of EAϵ₀. Regarding the torsional wave, its velocity is $\sim \sqrt{\mathit{\delta }g{L}_{\mathrm{0}}}$ and the wavelength ∼2L₀. The period then reads $T_{\mathrm{Tor}}\approx \mathrm{2}\sqrt{L_{\mathrm{0}}/\mathit{\delta }g}\approx \mathrm{0.4}$ s. The last two modes in Table 3 have 0 eigenvalues and correspond to free rotations as a rigid body of the full structure around its mean line.

Table 3 Stability results of the stationary solution. N=2.

n/a: not applicable

5.2 Nominal operation of a RAWE machine with fixed ends

We now investigate the dynamic response of the helix when the torque at the upper end varies as

$$\begin{array}{}\text{(69)}& m_N\left(\stackrel{\mathrm{̃}}{t}\right)=m_{N\mathrm{0}}\left[\mathrm{1}-\exp \left(-\stackrel{\mathrm{̃}}{t}/{\stackrel{\mathrm{̃}}{t}}_{\mathrm{0}}\right)\right],\end{array}$$

with m_N0=2.88 and ${\stackrel{\mathrm{̃}}{t}}_{\mathrm{0}}=\mathrm{2}$, which correspond to 25 N m and 2.6 s. For the controller in Eq. (40), we used $k_{\mathrm{1}}=k_{\mathrm{2}}=\mathrm{0.5}$ and a target angular velocity of 120 rpm, which is the nominal value for the operation of the RAWE machine. Both ends of the structure are fixed according to the boundary conditions (35)–(36) and Eq. (66). The initial condition of the simulation is the stationary solution explained in Sect. 5.1. System (61) was integrated numerically with a Runge–Kutta algorithm and for N=2. Since the torques at the ends of the helix depend on time, the term ${\stackrel{\mathrm{̃}}{m}}_i^{\mathrm{0}}$ in Eq. (56) involves the third derivative of α at s=0. For its evaluation at every time step during the numerical integration, we used backwards finite-difference approximations.

Figure 6a shows the evolution of $\dot{\mathit{\alpha }}$ for the three nodes of the structure (dotted, dashed, and solid lines) in the simulation. The controller successfully adjusted the torque on the ground (panel c) to make the structure rotate at the target angular velocity, which is 120 rpm. As shown in panel (b), the tension is almost constant throughout the helix and equal to around 200 N. The velocity of the nodes (not shown) tends to 0, and there is no lateral displacement. Therefore, the long-term state of the helix has its central line inside the x_E−z_E plane in a stationary condition, and the angles α for all the nodes rotate at a constant rate. The conclusions of this paragraph have been corroborated by running simulations with a higher number of elements (N=3 and N=5) and faster raising time of the torque at the upper end (${\stackrel{\mathrm{̃}}{t}}_{\mathrm{0}}=\mathrm{0.5}$ and ${\stackrel{\mathrm{̃}}{t}}_{\mathrm{0}}=\mathrm{1}$).

Figure 6 Panels (a), (b), and (c) show the evolution of $\dot{\mathit{\alpha }}$ at the nodes and the target angular velocity (thick solid line), the torque at the ground, and the tensions at the nodes, respectively. The ends of the structure are fixed.

5.3 Dynamics of the RAWE machine with a mobile end

The real operation of a RAWE machine is more complex than the scenario described in Sect. 5.2 because the upper-end point is mobile. Although an auxiliary kite provides additional lift and helps to anchor it, the motion of the rotor and wind fluctuations naturally drive a motion for the upper end of the structure. This mechanism injects energy into the helix, and, as shown below, its dynamics depend heavily on the characteristic frequency of the forcing. Insight into the behavior of the RAWE machine can be obtained by using the model of this work and still avoiding the coupling with a rotor by considering the following boundary condition for the upper end:

$$\begin{array}{}\text{(70)}& \begin{array}{rl}\stackrel{\mathrm{̃}}{\mathit{r}}(\mathrm{1},\stackrel{\mathrm{̃}}{t})=& -\left[(\mathrm{1}+{\mathit{ϵ}}_{\mathrm{0}})\cos {\mathrm{\Gamma }}_{\mathrm{0}}+r_{\mathrm{1}}\sin {\mathrm{\Gamma }}_{\mathrm{0}}\sin {\mathit{\omega }}_{\mathrm{1}}\stackrel{\mathrm{̃}}{t}\right]{\mathit{i}}_E\\ & +r_{\mathrm{1}}\cos {\mathit{\omega }}_{\mathrm{1}}\stackrel{\mathrm{̃}}{t}{\mathit{j}}_E\\ & -\left[(\mathrm{1}+{\mathit{ϵ}}_{\mathrm{0}})\sin {\mathrm{\Gamma }}_{\mathrm{0}}-r_{\mathrm{1}}\cos {\mathrm{\Gamma }}_{\mathrm{0}}\sin {\mathit{\omega }}_{\mathrm{1}}\stackrel{\mathrm{̃}}{t}\right]{\mathit{k}}_E.\end{array}\end{array}$$

It corresponds to a circular motion of the upper point of the structure with a radius and a frequency of R₁=r₁L₀ and $f_{\mathrm{1}}=\left({\mathit{\omega }}_{\mathrm{1}}/\mathrm{2}\mathit{\pi }\right)\sqrt{g/L_{\mathrm{0}}}$, respectively. Its effect on the dynamics of the structure has been studied experimentally and numerically.

The experiment, which was carried out inside the laboratory, used a structure with the characteristics in Tables 1 and 4. Although similar to the one considered in Sect. 5.2, we used a structure with a shorter length and lower mass due to space constraints in the laboratory. The nominal tension was also lowered to 80 N, and, according to Fig. 4 and Eq. (64), the bending stiffness is also lower (EI=6483 N m²). The nominal tension was adjusted by imposing an adequate distance between the tips of the structure. For this test, the torque was applied into the helix at the ground by a motor that was operated at a steady angular velocity. To impose the circular periodic motion given by Eq. (70), the helix’s upper-end point was attached to an eccentric arm that was powered with a second motor (see Fig. 2e). A radius of R₁=0.1 m was used, and the frequency f₁ was varied. An electromagnetic brake countered the torque on the eccentric arm. This is opposite to how the system operates in the generation mode of a RAWE machine, where the torque is generated by the rotor and countered by the generator on the ground. Nonetheless, the resulting torque in the helix during the experiment was the same as in a normal operation, and this experimental configuration avoided the need for having another motor with its power supply on the arm. Regarding diagnostics, the axial tension on the helix was measured with load cells on the attachment points at the ground.

Table 4 Parameters used in the analysis of solutions with a mobile end.

As shown in Fig. 7, the frequency f₁ of the eccentric arm was increased in multiple steps, and the total tension on the structure was measured with load cells. The lateral displacement of the helix was also monitored by using multiple cameras. The experimental data revealed that the tension oscillated at a frequency close to the driving frequency. The oscillations are regular for low frequencies and irregular for a driving frequency close to 5 Hz. It was also observed in the experiment that the amplitude of the oscillations, i.e., the lateral displacement of the central point of the structure, increased with the driving frequency. When the driving frequency reached 5 Hz, the structure collapsed because multiple tethers ripped, and the experiment was stopped. Interestingly, and as shown in Fig. 7, changing the forcing frequency from 0.7 to 2.1 Hz did not result in an enhancement of the tension. Only when the frequency was increased with up to 3.5 Hz and beyond were higher values measured by the load cells. According to these experimental results and the simulation results presented below, we conclude that the periodic forcing of the upper end induced a resonance in the structure. The irregular oscillations observed at high frequencies are a signature of the nonlinearities excited by the high values of the tension. They reached values equal to 3 times the nominal tension before the helix collapsed. We also mention that the frequency that produced the collapse (5 Hz) is of the order of the natural lateral and torsional frequencies of the helix. For the helix in Table 3, which has slightly different physical properties (compare Tables 2 and 4), the natural frequencies of the lateral and torsional modes are 1/0.2 s = 5 Hz and 1/0.4 s = 2.5 Hz. This result highlights the importance of having access to numerical tools for predicting the natural frequencies of the helix.

Figure 7 Evolution of the tension (left) and the forcing frequency f₁ (right) in the experiment with a mobile upper end.

The same scenario was also studied by using the simulator. Besides the parameters in Tables 1 and 4, we imposed the boundary condition given in Eq. (70) and set N=2. We started the analysis at a very low forcing frequency, and, for the initial condition, we used the procedure explained in Sect. 5.1. System (61) was then integrated numerically until the helix exhibited regular oscillations. The output tension at s=0 was then post-processed to identify its maxima (T_max) in the long term (once the transient phase died out). Afterwards, the value of f₁ was increased, and the full process was repeated but now using as initial condition the state of the helix at the end of the previous simulation.

The result of the analysis is the bifurcation diagram with T_max versus f₁ shown in Fig. 8. Interestingly, it brings out the strengths and weakness of the numerical tool. Regarding the former, the simulator predicts a resonance of the structure similar to the one observed in the experimental data. However, the maximum tension in the simulator rises for driving frequencies above around 9 Hz, whereas in the experiment an enhancement of the tension was detected for lower frequencies (above 3.5 Hz). The inset of Fig. 8 shows the evolution of the tension at s=0 for f₁=10 Hz in the long term (after the transient has died out). Similarly to the experimental results obtained for low frequencies, the tension oscillates regularly and at a frequency close to f₁. However, unlike the experimental data, the amplitude of the oscillations in the tension are small (a few newtons), whereas in the experiment they were of the order of tens of newtons. We also mention that the actual version of the simulator assumes that every cross section of the helix is unshearable. Therefore, the simulator cannot be used to predict the collapse of the helix. In fact, taking into account that the maximum tension in the experiment was about 200 N, it is likely that in the experiment we only observed the initial rise of the resonance curve shown in Fig. 8. For a tension of about 200 N, the driving frequency is about 11 Hz in Fig. 8.

Therefore, we conclude that the numerical tool is able to capture some interesting phenomena correctly, and it can be used to predict them qualitatively. However, if the interest is in getting quantitative values, then it should be used carefully. This conclusion is reasonable because the simulator relies on a set of simplifying hypotheses and also needs to be fed with inputs that have a certain level of uncertainty. As shown by the results of the experimental characterization of the structure in Sect. 4, the bending and torsional stiffness of the tensegrity structure is difficult to model. The stiffness depends on the tension and, in the case of torsional stiffness, it presents hysteresis. Since such hysteresis was found in the experiments, the results of this work open the question on how to incorporate such a phenomenon into the structural model. Theoretical models based on additional ordinary differential equations and artificial neural networks to fit complex hysteretic behaviors are two potential approaches. Interestingly, the former has recently been used in AWE to model the dynamic stall (Castro-Fernandez et al., 2024), which is another hysteretic phenomenon but in this case on unsteady aerodynamics. Another important topic is the modeling of the damping in the system, which could be improved in future work. The aerodynamic drag, already considered in the model, dissipates energy, but there are other mechanisms that may also be included. An example is the friction between elements of the tensegrity structure, which can affect the response of the system shown in Fig. 8. Consequently, the predictive capability of simulations is limited. Nonetheless, each simulation for a given f₁ has a relatively low computational cost (around 1 min). This opens the possibility to perform analysis varying the physical parameters that are difficult to estimate and tune them in the simulator to reproduce the experimental results quantitatively.

Figure 8 Maximum values of the tension at s=0 versus the forcing frequency in the simulations. The inset shows the evolution of the tension at s=0 for the forcing frequency f₁= 10 Hz.

6 Conclusions

The results of this work show that the proposed model based on Cosserat theory is appropriate for getting insight into the dynamics of the helical structure of a RAWE machine. It offers some advantages as compared to previous models based on springs and disks or a large set of masses linked with springs and dampers. For instance, unlike previous models, our model captures the three interesting motions of the helix, thus improving the fidelity. It provides a couple of nonlinear partial differential equations that can be solved with well-known numerical methods like the finite-element method used in this work. Such a compact formulation, which includes stretching, bending, and torsional effects, highlights the role of the different terms in the dynamics of the structure and their coupling. It has also been used here to estimate the characteristic times for the propagation of longitudinal, lateral, and torsional waves in the helix. Moreover, the airborne rotor and the ground generator enter the model as boundary conditions that include the position and velocity of the end sections of the structures and the external torque applied to them. The model was completed with the axial, bending, and torsional stiffness of a real AWE machine obtained by conducting three dedicated experiments in the laboratory. Since axial, bending, and torsional stiffness is essential information for RAWE machine simulators, such experimental work fills an important gap in the field. The relation between the bending and the torsional stiffness and the tension was measured for the tensegrity structure. It was also shown that the torsional stiffness depends on both the tension without torsion and the history of the structure as the tension is increased. A hysteric behavior was measured, which poses a challenge for the modeling of the helix.

The analysis of the stationary solution with no torsion and fixed ends, which was found independently by using a shooting algorithm, allowed us to verify the correct implementation of the finite-element code. A linear stability analysis was performed to identify longitudinal, lateral, and torsional modes and their natural periods. The quickest mode of the RAWE machine of SomeAwe corresponds to a longitudinal mode with a period of 0.03 s, whereas lateral and torsional modes exhibit natural oscillations with 0.2 and 0.4 s. Therefore, the natural frequencies of the helix predicted by the simulator are 33.3, 5, and 2.5 Hz. Simple estimations based on classical results for beams are in agreement with these results and constitute a second test for the correct implementation of the code.

The simulation tool provided interesting information about the nominal operation of the RAWE machine like the transient of the angular velocities of the cross sections of the helix, the tension, and the torques. A numerical analysis was carried out by keeping the two end points of the structure fixed. It assumed that the rotor imposes a time-dependent external torque that approaches 25 N m, which is a typical value for the RAWE machine under consideration. It was shown that a simple proportional-derivative controller for the torque at the ground generator can stabilize the angular spinning velocity of the structure at the target value (120 rpm). For this configuration, the tension is almost constant in time and throughout the structure. The evolution of the variables and their values are well aligned with the data collected in field testing by SomeAWE Labs.

In a second numerical analysis, the lower end of the helix was kept fixed, and the upper end was mobile. Although an auxiliary kite is used to anchor the rotor, in a real RAWE machine the upper-end point is not fixed due to wind velocity fluctuations and the dynamic coupling of the structure and the rotor. This scenario was mimicked by imposing a circular periodic motion to the upper end of the structure with a forcing frequency f₁ and small amplitude. The experimental and simulation work revealed that there is a resonance that, in the case of the experiment, resulted in the collapse of the structure for a frequency of 5 Hz and for a tension level of about 200 N. The simulator captured essential features of the experiment, like the resonance and the order of magnitude of the critical frequency. The simulator cannot reproduce the collapse of the structure because the model was constructed by assuming that the cross sections of the helix are unshearable. For this reason, it was able to explore higher frequencies and tension values with the simulator, reaching a maximum above 900 N.

Appendix A: Auxiliary calculations

By assuming that the first derivative is continuous across the elements, we find the following relation between ${\stackrel{\mathrm{̃}}{\mathit{r}}}_j^{\prime \prime }\left(\stackrel{\mathrm{̃}}{t}\right)$ and ${\stackrel{\mathrm{̃}}{\mathit{r}}}_j\left(\stackrel{\mathrm{̃}}{t}\right)$ (Press et al., 1992):

$$\begin{array}{}\text{(A1)}& {\mathbf{B}}_r\left(\begin{array}{c}{\stackrel{\mathrm{̃}}{\mathit{r}}}_{\mathrm{1}}^{\prime \prime }\\ {\stackrel{\mathrm{̃}}{\mathit{r}}}_{\mathrm{2}}^{\prime \prime }\\ {\stackrel{\mathrm{̃}}{\mathit{r}}}_{\mathrm{3}}^{\prime \prime }\\ \mathrm{\dots }\\ {\stackrel{\mathrm{̃}}{\mathit{r}}}_{N-\mathrm{1}}^{\prime \prime }\end{array}\right)={\mathbf{C}}_r\left(\begin{array}{c}{\stackrel{\mathrm{̃}}{\mathit{r}}}_{\mathrm{0}}\\ {\stackrel{\mathrm{̃}}{\mathit{r}}}_{\mathrm{1}}\\ {\stackrel{\mathrm{̃}}{\mathit{r}}}_{\mathrm{2}}\\ \mathrm{\dots }\\ {\stackrel{\mathrm{̃}}{\mathit{r}}}_N\end{array}\right),\end{array}$$

with

$$\begin{array}{}\text{(A2)}& {\displaystyle }\begin{array}{rl}& {\mathbf{B}}_r=\\ & {\left(\begin{array}{cccccc}\mathrm{2}\mathit{I}/\mathrm{3}& \mathit{I}/\mathrm{6}& \mathrm{0}& \mathrm{0}& \mathrm{\dots }& \mathrm{0}\\ \mathit{I}/\mathrm{6}& \mathrm{2}\mathit{I}/\mathrm{3}& \mathit{I}/\mathrm{6}& \mathrm{0}& \mathrm{\dots }& \mathrm{0}\\ \mathrm{0}& \mathit{I}/\mathrm{6}& \mathrm{2}\mathit{I}/\mathrm{3}& \mathit{I}/\mathrm{6}& \mathrm{\dots }& \mathrm{0}\\ \mathrm{0}& \mathrm{\dots }& \mathrm{\dots }& \mathrm{\dots }& \mathrm{\dots }& \mathrm{\dots }\\ \mathrm{0}& \mathrm{\dots }& \mathrm{\dots }& \mathrm{\dots }& \mathit{I}/\mathrm{6}& \mathrm{2}\mathit{I}/\mathrm{3}\end{array}\right)}_{\mathrm{3}(N-\mathrm{1})\times \mathrm{3}(N-\mathrm{1}),}\end{array}\text{(A3)}& {\displaystyle }\begin{array}{rl}& {\mathbf{C}}_r=\\ & {\left(\begin{array}{cccccc}\mathit{I}& -\mathrm{2}\mathit{I}& \mathit{I}& \mathrm{0}& \mathrm{\dots }& \mathrm{0}\\ \mathrm{0}& \mathit{I}& -\mathrm{2}\mathit{I}& \mathit{I}& \mathrm{\dots }& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathit{I}& -\mathrm{2}\mathit{I}& \mathit{I}& \mathrm{0}\\ \mathrm{0}& \mathrm{\dots }& \mathrm{\dots }& \mathrm{\dots }& \mathrm{\dots }& \mathrm{\dots }\\ \mathrm{0}& \mathrm{\dots }& \mathrm{\dots }& \mathrm{\dots }& \mathrm{\dots }& \\ \mathrm{0}& \mathrm{0}& \mathrm{\dots }& \mathit{I}& -\mathrm{2}\mathit{I}& \mathit{I}\end{array}\right)}_{\mathrm{3}(N-\mathrm{1})\times \mathrm{3}(N+\mathrm{1}),}\end{array}\end{array}$$

and we used the boundary conditions ${\mathit{r}}^{\prime \prime }(\mathrm{0},\stackrel{\mathrm{̃}}{t})={\mathit{r}}^{\prime \prime }(\mathrm{1},\stackrel{\mathrm{̃}}{t})=\mathrm{0}$. For α, the continuity condition for α^′ across the elements gives

$$\begin{array}{}\text{(A4)}& {\mathbf{B}}_{\mathit{\alpha }}\left(\begin{array}{c}{\mathit{\alpha }}_{\mathrm{0}}^{\prime \prime }\\ {\mathit{\alpha }}_{\mathrm{1}}^{\prime \prime }\\ \mathrm{\dots }\\ {\mathit{\alpha }}_{N-\mathrm{1}}^{\prime \prime }{\mathit{\alpha }}_N^{\prime \prime }\end{array}\right)={\mathbf{C}}_{\mathit{\alpha }}\left(\begin{array}{c}{\mathit{\alpha }}_{\mathrm{0}}\\ {\mathit{\alpha }}_{\mathrm{1}}\\ \mathrm{\dots }\\ {\mathit{\alpha }}_{N-\mathrm{1}}\mathrm{\dots }{\mathit{\alpha }}_N\end{array}\right)+\left(\begin{array}{c}-m_{\mathrm{0}}\\ \mathrm{0}\\ \mathrm{\dots }\\ \mathrm{0}\\ m_N\end{array}\right),\end{array}$$

with

$$\begin{array}{}\text{(A5)}& {\displaystyle }\begin{array}{rl}& {\mathbf{B}}_{\mathit{\alpha }}=\\ & {\left(\begin{array}{cccccc}\mathrm{1}/\mathrm{3}& \mathrm{1}/\mathrm{6}& \mathrm{0}& \mathrm{0}& \mathrm{\dots }& \mathrm{0}\\ \mathrm{1}/\mathrm{6}& \mathrm{2}/\mathrm{3}& \mathrm{1}/\mathrm{6}& \mathrm{0}& \mathrm{\dots }& \mathrm{0}\\ \mathrm{0}& \mathrm{1}/\mathrm{6}& \mathrm{2}/\mathrm{3}& \mathrm{1}/\mathrm{6}& \mathrm{\dots }& \mathrm{0}\\ \mathrm{0}& \mathrm{\dots }& \mathrm{\dots }& \mathrm{\dots }& \mathrm{\dots }& \mathrm{\dots }\\ \mathrm{0}& \mathrm{\dots }& \mathrm{\dots }& \mathrm{1}/\mathrm{6}& \mathrm{1}/\mathrm{3}& \mathrm{1}/\mathrm{6}\\ \mathrm{0}& \mathrm{\dots }& \mathrm{\dots }& \mathrm{\dots }& \mathrm{1}/\mathrm{6}& \mathrm{1}/\mathrm{3}\end{array}\right)}_{(N+\mathrm{1})\times (N+\mathrm{1}),}\end{array}\text{(A6)}& {\displaystyle }\begin{array}{rl}& {\mathbf{C}}_{\mathit{\alpha }}\\ & ={\left(\begin{array}{cccccc}-\mathrm{1}& \mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{\dots }& \mathrm{0}\\ \mathrm{1}& -\mathrm{2}& \mathrm{1}& \mathrm{0}& \mathrm{\dots }& \mathrm{0}\\ \mathrm{0}& \mathrm{1}& -\mathrm{2}& \mathrm{1}& \mathrm{0}& \mathrm{0}\\ \mathrm{\dots }& \mathrm{\dots }& \mathrm{\dots }& \mathrm{\dots }& \mathrm{\dots }& \mathrm{\dots }\\ \mathrm{0}& \mathrm{\dots }& \mathrm{\dots }& \mathrm{1}& -\mathrm{2}& \mathrm{1}\\ \mathrm{0}& \mathrm{\dots }& \mathrm{\dots }& \mathrm{0}& \mathrm{1}& -\mathrm{1}\end{array}\right)}_{(N+\mathrm{1})\times (N+\mathrm{1}),}\end{array}\end{array}$$

where we used the boundary conditions ${\mathit{\alpha }}^{\prime }(\mathrm{0},\stackrel{\mathrm{̃}}{t})=m_{\mathrm{0}}$ and ${\mathit{\alpha }}^{\prime }(\mathrm{1},\stackrel{\mathrm{̃}}{t})=m_N$.