Aeroelastic Modeling of Full-Scale Airborne Wind Turbine Based on FSI Approach

1. Introduction

Since the last decade, new trends in the wind turbine industry have steadily grown with the objective of novel harvesting technologies in favor of renewable energies [1]. Especially the airborne wind energy systems (AWESs) [2] have many undeniable advantages over traditional wind turbines to capture untapped energy [3] referred to as airborne wind energy (AWE) [4] of atmosphere layers. On the other hand, tower-based wind turbines have limited potential to harness the consistent winds blowing at elevated heights [5] due to the high-rise buildings, tower constraints, uneven terrain [6] and meteorological conditions. Subsequently, a technique of buoyant airborne turbine (BAT) is introduced as a unique representative of the AWT. Technically, the hardware configuration of the BAT energy system consists of a small size of rotor placed inside a helium-filled buoyant shell [7] while the entire assembly is lifted to operate at 400 m height. Despite beneficial yield, the AWT often confronts the challenges of increased wind speeds that might result in aeroelastic instabilities [8] of the wind converter which can be devastating to the rotor blades. Hence, accurate FSI modeling is crucial to investigate the BAT energy system.

The aerodynamic behavior of the BAT can be characterized by the ducted wind turbine (DWT) [9,10] or diffuser augmented wind turbine (DAWT) [11,12]. In contrast to the solid duct in DWT, the BAT is equipped with a flexible buoyant shell. Consequently, the enclosed helium gas generates the aerostatic lift to raise the BAT system at desired altitudes. From a numerical point of view [13], many researchers opted to investigate the aerodynamic performance of the ducted or shell-based turbines using a numerical approach. Different studies were carried out to evaluate the DWTs by employing the blade element momentum (BEM) and computational fluid dynamics (CFD) numerical approaches [14]. The Reynolds-averaged Navier–Stokes (RANS) based CFD simulations were conducted while enclosing the rotor inside an airfoil-based shroud [15], duct [16] and stepped duct [17]. Likewise, in the BAT system, many authors incorporated the RANS-CFD approach for the evaluation of the aerodynamic performance of 2-bladed rotor [18,19] and 3-bladed rotor [20] in a buoyant shell configuration. However, the majority of these analyses focus on the aerodynamic aspect of the DWTs in terms of performance analysis. For this reason, in this present study, the CFD analysis together with finite element method (FEM) analysis [21] was performed to analyze the structural integrity of AWT through FSI analysis.

The applicability of FSI analysis that combines the fluid and structural modules to compute the power output, stresses, and deflection has been broadly discussed in the literature. From the FSI modeling perspective of the rotor [22], most researchers extensively applied the CFD and FEM techniques to measure the aeroelastic responses on horizontal axis wind turbines (HAWTs) [23] and NREL phase VI wind turbines [24]. Similarly, the elasto-flexible membrane blade [25], vertical axis wind turbines (VAWTs) [26] and floating offshore wind turbines [27] were also put forward for numerical analysis. A CFD study was conducted using 20 kW NREL wind turbine by implementing morphing features [28]. The reliability of 15MW blades rotor of HAWT was assessed through FSI Simulations by combining Kriging model and Monte Carlo simulation. The promising results of our study can be expand into the creation of digital twins for wind turbines [29]. The optimization of the structural topology of H-Rotor blade based on a one-way FSI approach by coupling URANS with density based topology optimization method. The blade mass can be easily reduced by around 60% with the devised approach [30]. The slender blades of a full-scale floating 15 MW wind turbine were analyzed by establishing CFD-FEA approach under wave-induced motions on the wind field [31]. A general FSI framework combining CFD, FEA and multi-body dynamics (MBD) was proposed to evaluate the detailed stress distributions on the composite structures of wind turbine blade [32].

The aeroelastic response of the NREL 5 MW wind turbine under diverse inflow conditions by integrating an improved actuator line model ALM with an equivalent beam model [33]. The dynamic response of the H-type floating VAWT platform was explored using the coupled CFD-FEM method [34]. At present, artificial intelligence (AI) [35] playing an important role in the design and optimization of offshore floating wind turbines [36]. A fully coupled aero-elastic-hydro-mooring framework by CFD and FEM, NREL 5 MW HAWT mounted on a semi-submersible platform, was investigated under combined wind and wave loading [36]. The ANSYS-based one-way FSI study was conducted to assess the blade lifespan of two small-scale wind turbines under different wind speeds [37]. Recently, a deep learning based framework that integrates CNN-BEM-FEM models to perform FSI analysis for a composite blade of a tidal turbine [38]. Two-way FSI for a 4.5 MW wind turbine was carried out by integrating the Large eddy simulation (LES) and structural dynamics to investigate the aeroelastic responses considering the flexibility of the blade and tower [39]. The simulated time series via BEM and CFD based modeling approaches can be incorporated with the data-driven approach of high fidelity aeroelastic analysis [40].

Most of these aeroelastic analyses generally provide reasonable results by taking into account the geometric nonlinearities to estimate the large blade deflections and blade-tower interaction. Despite the undertaken research, very limited literature is available regarding the FSI modeling of the ducted or shell turbine. There only a few attempts have been made to understand the FSI behavior of a buoyant shell [41] and ducted turbine [42]. Moreover, the available literature is incapable of providing detailed information regarding the aerodynamic and nonlinear structural design [43] of the rotor in ducted/shell configuration. Therefore, it is imperative to account for these research gaps and significantly contribute to a more detailed study to gain insight into the three-dimensional (3D) aeroelastic modeling of the AWT. Hence, a one-way FSI analysis that couples [44] the higher resolution methods of CFD and FEM is performed by investigating the structural behaviors [45] together with the safe functionality of the rotor. Most importantly, the characteristics of the analyzed problem are deeply examined with respect to the complex loading of the shell around the rotor blades.

The prime aim of this paper is to present a one-way FSI analysis for an airborne wind turbine at full scale. The task of a reliable aeroelastic model is accomplished by coupling the CFD and FEM models of the blade as well as the shell geometry. To maintain the proper stiffness, the blade is reinforced geometrically through a shear web. Similarly, the structural aspects of the shell’s deformable envelope are supposed to be a semi-rigid structure. Furthermore, the blade and shell were simplified by orthotropic and anisotropic material properties, respectively. The pressure loads in terms of fully converged steady-state are computed using CFD and then directly transferred to the structural surface. To determine the structural responses, nonlinear static structural analysis is performed using the developed FEM model. The main contribution is being made to enhance the clarity of the current study by explaining the rotor and shell models in parallel from CFD and FEM viewpoints throughout the paper. The results obtained using the one-way FSI approach are based on six operational conditions while the tip deformations and maximum stresses are found to be within structural limits. In addition, the numerical framework of the presented FSI model can be applied to many similar applications, such as ducted wind devices, a ducted propeller for urban air mobility and tidal turbines.

2. Bat Energy System

A 3-bladed upwind rotor consisting of the twisted and tapered shape of blades [46] was employed in a ring-shaped buoyant shell. The outer blade shape was composed of three different airfoils and was originally designed to operate at an airborne altitude of 400 m height. The radius of the rotor blade $\left(R_b=\sqrt{P_r/0.5{\rho }_h\pi {V_r}^3{C}_{\mathrm{p}}{\eta }_{\mathrm{g}}}\right)$ has span of 2.51 m that is a representative of a power output of 30kW under rated tip speed ratio (TSR) ${(\lambda }_r=\omega R_b/V_r)$of 6.5. On the other hand, the outer aerodynamic shape of the shell [20] was modeled using a Delft University airfoil (DU 93-W-210). This airfoil was selected because of its characteristics of a moderate cambered of 2.85% and a maximum relative thickness of 21%. The main aim of the installation of the shell around the rotor blades is to provide beneficial velocity speed-up at the rotor plane. The major advantage of aerostatic lift comes from the buoyant force ${(F}_b={\rho }_h{S}_v\mathrm{g})$. This upward buoyancy force is generated from the helium gas enclosed inside the internal volume of the shell. The two main aerodynamic components of the BAT energy system i.e., rotor and shell are shown in Figure 1. The details of the geometric models will be discussed in subsection 4.1 whereas the specifications of the rotor and shell are outlined in Table 1.

3. One-Way FSI Analysis

The modeling framework of the fluid-structure interaction (FSI) analysis [47] is composed of the coupling of the fluid model (CFD model) and structural model (FEM model). There are two main approaches of coupling to solve the model for the FSI analysis, namely two-way coupling and one-way coupling. The two-way FSI is a direct coupling approach, the fluid model and structural model are solved simultaneously but separately to acquire and exchange the resulting data of pressure loads and displacement, respectively. On the other hand, in a one-way coupling model, the fully converged aerodynamic loads in the steady-state are obtained from CFD and are transferred to the FEM model to predict the response behavior of a single blade. The procedure is executed in such a way that the resulting displacement data is not mapped back to the fluid model. Furthermore, the one-way coupling is particularly suited for initial modeling purposes within a reasonable computing cost. For this reason, the one-way FSI coupling method is selected for aeroelastic modeling in this study as explained in Figure 2.

3.1. Fluid Model

The fluid model was developed using the ANSYS-CFX^® software [48] to evaluate the CFD simulations. The CFD solver employs the finite volume method (FVM) approach for spatial discretization. The steady and incompressible RANS equations are incorporated to calculate the averaging flow quantities as wind loads on the resulting model. Mathematically, the governing equations of the mass and momentum are given as [49,50]:

$${\displaystyle \frac{\partial u_i}{\partial x_i}}=0$$

(1)

$$\rho u_j{\displaystyle \frac{\partial \left(u_i\right)}{\partial x_j}}=-{\displaystyle \frac{\partial p}{\partial x_i}}+\mu {\displaystyle \frac{{\partial }^2{u}_i}{\partial x_j\partial x_j}}+{\displaystyle \frac{\partial }{\partial x_j}}\left(-\rho \overline{{\acute{u}}_i{\acute{u}}_j}\right)+{\rho f}_{b,i}$$

(2)

Here the terms $\rho ,p,\mu ,$ and $f_{b,i}$denote the fluid density, pressure, dynamic viscosity and body force, respectively. The unknown term $\left(-\rho \overline{{\acute{u}}_i{\acute{u}}_j}\right)$represents the Reynolds stress. To simulate turbulent flow, the shear stress transport (SST) model [51] was enabled at preprocessing level. The SST model was chosen because of its well-known capability in handling the external flows over the surface of wind turbines and wings [52]. The SST turbulence model is a two-equations model that expresses the two additional transport equations in the following form of Eq. (3) and Eq. (4):

$${\displaystyle \frac{\partial }{\partial x_i}}\left(u_i\rho k\right)={\displaystyle \frac{\partial }{\partial x_i}}\left({\mu }_k{\displaystyle \frac{\partial }{\partial x_i}}k\right)+{\tilde{P}}_k-{\beta }^{*}\rho \omega k$$

(3)

$${\displaystyle \frac{\partial }{\partial x_i}}\left(u_i\rho \omega \right)={\displaystyle \frac{\partial }{\partial x_i}}\left({\mu }_{\omega }{\displaystyle \frac{\partial }{\partial x_i}}\omega \right)+\alpha \rho S^2-\beta \rho {\omega }^2+2\rho \left(1-F_1\right){\sigma }_{\omega 2}{\displaystyle \frac{1}{\omega }}{\displaystyle \frac{\partial }{\partial x_i}}k{\displaystyle \frac{\partial }{\partial x_i}}\omega$$

(4)

The left-hand side of Eqs. (3) and (4) represent the change and transport of the turbulent kinetic energy $\left(k\right)$ and turbulent dissipation rate $\left(\omega \right)$, respectively. Likewise, on the right-hand side of Eqs. (3) and (4) express the rate of production, turbulent diffusion, and rate of diffusion of k and ω, respectively. The term $F_1$in Eq. (4) defines the blending function which acts as a switch for certain values between standard k–ε and k–ω turbulence models.

3.2. Structural Model

The FEM approach of the ANSYS static structural module [53] allows breaking the spatial domain of the model into quadrilaterals and triangle elements. The governing equation in static structural analysis (${\displaystyle \frac{d^{2}u}{{dt}^2}}=0$) is expressed in differential form as conservation of momentum and is known as Cauchy’s equation of motion [54].

$$\nabla .{\sigma }_s+{\rho f}_{b,i}=0$$

(5)

Here $\rho$ is the density, $f_{b,i}$stands for the body force per unit mass in i-th direction, $\nabla$ refers to the del, or nabla operator and ${\sigma }_s$ is the Cauchy stress tensor. Hence, the partial differential expression of Cauchy’s equation that is governing the static equilibrium of the body takes the following forms:

$${\displaystyle \frac{{\partial \sigma }_{xx}}{\partial x}}+{\displaystyle \frac{{\partial \sigma }_{xy}}{\partial y}}+{\displaystyle \frac{{\partial \sigma }_{xz}}{\partial z}}+f_x=0$$

(6)

$${\displaystyle \frac{{\partial \sigma }_{yx}}{\partial x}}+{\displaystyle \frac{{\partial \sigma }_{yy}}{\partial y}}+{\displaystyle \frac{{\partial \sigma }_{yz}}{\partial z}}+f_y=0$$

(7)

$${\displaystyle \frac{{\partial \sigma }_{zx}}{\partial x}}+{\displaystyle \frac{{\partial \sigma }_{zy}}{\partial y}}+{\displaystyle \frac{{\partial \sigma }_{zz}}{\partial z}}+f_z=0$$

(8)

Here the $f_x={\rho f}_{b,x}$ represents the body force per unit volume.With the assumption of large deformation, the strain-displacement relations reflect the following form of the compatibility equations [55]:

$${\epsilon }_{xx}={\displaystyle \frac{{\partial u}_x}{\partial x}}+{\displaystyle \frac{1}{2}}\left[{\left({\displaystyle \frac{{\partial u}_x}{\partial x}}\right)}^2+{\left({\displaystyle \frac{{\partial u}_y}{\partial x}}\right)}^2+{\left({\displaystyle \frac{{\partial u}_z}{\partial x}}\right)}^2\right]$$

(9)

$${\epsilon }_{yy}={\displaystyle \frac{{\partial u}_x}{\partial y}}+{\displaystyle \frac{1}{2}}\left[{\left({\displaystyle \frac{{\partial u}_x}{\partial y}}\right)}^2+{\left({\displaystyle \frac{{\partial u}_y}{\partial y}}\right)}^2+{\left({\displaystyle \frac{{\partial u}_z}{\partial y}}\right)}^2\right]$$

(10)

$${\epsilon }_{zz}={\displaystyle \frac{{\partial u}_z}{\partial z}}+{\displaystyle \frac{1}{2}}\left[{\left({\displaystyle \frac{{\partial u}_x}{\partial z}}\right)}^2+{\left({\displaystyle \frac{{\partial u}_y}{\partial z}}\right)}^2+{\left({\displaystyle \frac{{\partial u}_z}{\partial z}}\right)}^2\right]$$

(11)

$${\epsilon }_{xy}={\displaystyle \frac{1}{2}}\left({\displaystyle \frac{{\partial u}_x}{\partial y}}+{\displaystyle \frac{{\partial u}_y}{\partial x}}\right)+{\displaystyle \frac{1}{2}}\left[{\displaystyle \frac{{\partial u}_x}{\partial x}}{\displaystyle \frac{{\partial u}_x}{\partial y}}+{\displaystyle \frac{{\partial u}_y}{\partial x}}{\displaystyle \frac{{\partial u}_y}{\partial y}}+{\displaystyle \frac{{\partial u}_z}{\partial x}}{\displaystyle \frac{{\partial u}_z}{\partial y}}\right]$$

(12)

$${\epsilon }_{xz}={\displaystyle \frac{1}{2}}\left({\displaystyle \frac{{\partial u}_x}{\partial z}}+{\displaystyle \frac{{\partial u}_z}{\partial x}}\right)+{\displaystyle \frac{1}{2}}\left[{\displaystyle \frac{{\partial u}_x}{\partial x}}{\displaystyle \frac{{\partial u}_x}{\partial z}}+{\displaystyle \frac{{\partial u}_y}{\partial x}}{\displaystyle \frac{{\partial u}_y}{\partial z}}+{\displaystyle \frac{{\partial u}_z}{\partial x}}{\displaystyle \frac{{\partial u}_z}{\partial z}}\right]$$

(13)

$${\epsilon }_{yz}={\displaystyle \frac{1}{2}}\left({\displaystyle \frac{{\partial u}_y}{\partial z}}+{\displaystyle \frac{{\partial u}_z}{\partial y}}\right)+{\displaystyle \frac{1}{2}}\left[{\displaystyle \frac{{\partial u}_x}{\partial y}}{\displaystyle \frac{{\partial u}_x}{\partial z}}+{\displaystyle \frac{{\partial u}_y}{\partial y}}{\displaystyle \frac{{\partial u}_y}{\partial z}}+{\displaystyle \frac{{\partial u}_z}{\partial y}}{\displaystyle \frac{{\partial u}_z}{\partial z}}\right]$$

(14)

The stress and strain relationship $\left(\left[\sigma \right]=\left[E\right]\left[\epsilon \right]\right)$ of orthotropic material [56] can be described as a constitutive equation with respect to an orthonormal coordinate system. The symmetry transformations of orthotropic materials about all three orthogonal planes are as follows [57]:

$$\left[\begin{array}{c}{\sigma }_{xx}\\ {\sigma }_{yy}\\ {\sigma }_{zz}\\ {\sigma }_{yz}\\ {\sigma }_{zx}\\ {\sigma }_{xy}\end{array}\right]=\left[\begin{array}{cccccc}{\displaystyle \frac{E_x(1-{\nu }_{yz}{\nu }_{zy})}{\xi }}& {\displaystyle \frac{E_x({\nu }_{yz}+{\nu }_{zx}{\nu }_{yz})}{\xi }}& {\displaystyle \frac{E_x({\nu }_{zx}+{\nu }_{yx}{\nu }_{zy})}{\xi }}& 0& 0& 0\\ {\displaystyle \frac{E_y({\nu }_{xy}+{\nu }_{xz}{\nu }_{zy})}{\xi }}& {\displaystyle \frac{E_y(1-{\nu }_{zx}{\nu }_{xz})}{\xi }}& {\displaystyle \frac{E_y({\nu }_{zy}+{\nu }_{zx}{\nu }_{xy})}{\xi }}& 0& 0& 0\\ {\displaystyle \frac{E_z({\nu }_{xz}+{\nu }_{xy}{\nu }_{yz})}{\xi }}& {\displaystyle \frac{E_z({\nu }_{yz}+{\nu }_{xz}{\nu }_{yx})}{\xi }}& {\displaystyle \frac{E_z(1-{\nu }_{xy}{\nu }_{yx})}{\xi }}& 0& 0& 0\\ 0& 0& 0& G_{yz}& 0& 0\\ 0& 0& 0& 0& G_{zx}& 0\\ 0& 0& 0& 0& 0& G_{xy}\end{array}\right]\left[\begin{array}{c}{\epsilon }_{xx}\\ {\epsilon }_{yy}\\ {\epsilon }_{zz}\\ {2\epsilon }_{yz}\\ {2\epsilon }_{zx}\\ 2{\epsilon }_{xy}\end{array}\right]$$

(15)

$$\xi =1-{\nu }_{xy}{\nu }_{yx}-{\nu }_{yz}{\nu }_{zy}-{\nu }_{zx}{\nu }_{xz}-{2\nu }_{xy}{\nu }_{yz}{\nu }_{zx}$$

(16)

Here $E_i$ is Young’s modulus in the i-th direction, $G_{ij}$ is the shear modulus in ij plane, ${\nu }_{ij}$ is the Poisson ratios.

4. Method of Solution

4.1. Geometric Models

The geometry of the rotor blade was modeled using the S-series of thick-airfoils family i.e., S-820, S-819 and S-821[58]. Similarly, the DU 93-W-210 airfoil and NACA-0012 airfoil were used to constitute the torus body of the shell and the projected wing, respectively. The geometric configuration, dimensions and shape of incorporated airfoils of the single blade and shell are shown in Figure 3. In view of the structural analysis, geometries were further segmented in multiple surfaces and a shear web was additionally introduced as illustrated in Figure 4 (a) and (b).

4.2. Mesh Models

For the fluid model in CFD analysis, a hexahedron mesh of high resolution was applied to constructed geometries based on 3D blocking topologies in ANSYS ICEM. First, the mesh generation of a single rotor was accomplished by adopting a C-type blocking topology. Second, the shell geometry was meshed using a series of O-grid topologies. The generated mesh of rotor and shell geometry is shown in Figure 5 (a) and (b). For the structural model in FEM analysis, a surface mesh of the rotor and shell was generated by a mean of 4-node shell-181 elements based on 2D blocking topologies illustrated in Figure 5 (c) and (d). To build a body-fitted mesh, a combination of quadrilateral and triangle cell shapes was incorporated within the ANSYS ICEM.

4.3. Boundary Conditions (BCs) and Domains

For the CFD model, the stationary domain$\left({\Omega }_{out}^s\right)$ was established to act as the outer fluid domain. The shell was modeled inside ${\Omega }_{out}^s$while the rotating domain$\left({\mathsf{\Omega }}_{in}^r\right)$ carrying the mesh model of the rotor blades was placed inside the gap into ${\Omega }_{out}^s$. In steady-state analysis, a rotor blade with a faster simulation approach of Frozen Rotor was applied to impart the associated momentum terms into flow unlike the actuator disc approach where blade geometries need not be modeled. The non-conformal mesh connections were defined at overlapped interfaces. The domain sizes and boundary conditions (BCs) of the inlet and out are depicted in Figure 6 (a). The pressure loads obtained from the CFD analysis were then input into the static structural model to act as pressure BC. The other BCs input for the FEM software was the rotor’s rpm, fixed support at the root region and remote displacement (rigid support) for the shell geometry as shown in Figure 6 (b) and (c).

4.4. Shell Elements and Material Distribution

For model simplification, orthotropic and Kevlar-49 material properties based on the linear elastic model were selected to examine the rotor blade and shell in static structural analysis, respectively, as given in Table 2. Generally, the thin blades of wind turbines are produced with composite layups. The orthotropic material as a subset of composite material is primarily intended to reduce the material modeling efforts in structural analysis.

The selection of the shell-181 element is the most suitable choice when the thickness of the structure is significantly smaller than the other dimensions of the structure. The shell thickness of the blade is defined as sectional property and is calculated by the following relation [59].

$$t_b=b_1-{\left(b_1^3-3I/2a\right)}^{1/3}$$

(17)

In general, the direction of the x-axis of shell-181 is oriented by the nodal connectivity running from I-node to J-node (see Figure 7 (a)). Therefore, the element orientation was determined by using an orientation definition of the element coordinate system corresponding to the surface geometry of the blade as shown in Figure 7 (b). A varying thickness ranging from 5 mm (root) to 2 mm (tip) was used for the blade surface while a uniform thickness of 1 mm was used for the shell surfaces. Referring to Figure 8 (a), a familiar distribution of chord length $\left(c_l\right)$ and twist angle $\left(\beta \right)$ is gradually changing from the root to the tip region of the blade. A similar trend just like geometric entities can be observed in the case of material thickness $\left(t_b\right)$ and mass distribution $\left(m_d\right)$ as shown in Figure 8 (b).

4.5. Mesh Sensitivity Analysis

An assessment of the sensitivity of the mesh was performed for both the CFD and FEM models well before the FSI analysis. The CFD models of the blade and shell were analyzed at rated conditions using four sets of mesh models and performance parameters were monitored in each simulation. The gradual refinement of the mesh model from mesh 1 to mesh 4 permits the selection of the appropriate mesh density in terms of performance indicator and iteration time. Finally, mesh 3 was selected based on the percentage difference (between mesh 3 and mesh 4) of <1% and average time saving of >7% as described in Table 3.

For the mesh density of the FEM blade model, a modal analysis was performed on four different sets of meshes. The modal analysis is the fundamental dynamic analysis that provides the natural frequencies and mode shapes of the blade. In each set of meshes from mesh 1 to mesh 4, the FEM model was constrained at the root region of the blade and six modal frequencies (Hz) were calculated and compared as shown in Table 4. Based on the percentage difference, the FEM mesh model with 74,404 mesh elements was finally selected for further in-depth analysis of the non-linear static structural analysis of the blade. Likewise, several FEM meshes of different mesh densities for the buoyant shell were analyzed and a mesh with an element size of 1,46,130 was selected for further analysis. However, for brevity, the mesh sensitivity table of the shell is not included here.

5. Results and Discussion

5.1. Model Validation

The CFD model of the shell rotor must be verified to correctly determine the pressure load on the surface of rotor blades and buoyant shell. For model validation, the CFD model of the shell rotor was simulated at the rated wind speed ($V_r=18\mathrm{m}/\mathrm{s}$) against varying TSRs ranging from $\lambda =3.5-9.5$. All the cases were run in steady-state flow for a minimum of 1200 iterations. The convergence criteria were set to $1\times {10}^{-4}$ to achieve a stable CFD solution. Next, the CFD model of the rotor in shell configuration was analyzed by processing the numerical data. The mechanical torque $\left(T_m\right)$ and thrust force $\left(F_T\right)$ acquired from the CFD model in shell configuration were used for the estimation of the total power output ${(P}_t=\omega T_m)$ produced by the rotor. Similarly, the power coefficient of the shell rotor ${(C}_{\mathrm{p},s})$, the rotor thrust coefficient ${(C}_{T,rot})$, rotor torque coefficient ${(C}_{Tr})$, can be calculated from the following relations:

$$C_{\mathrm{p},s}=P_t/0.5{\rho }_h{V}_r^3\pi R_b^2$$

(18)

$$C_{T,rot}=F_T/0.5{\rho }_h{V}_r^2\pi R_b^2$$

(19)

$$C_{Tr}=T_m/0.5{\rho }_h{V}_r^2\pi R_b^3$$

(20)

Here the ${\rho }_h$, $V_r$ and $R_b$denotes the density of air at 400 m height, rated wind velocity and the blade radius, respectively. This extracted CFD data of various simulations were employed to generate $C_{\mathrm{p},s}-C_{T,rot}$ curve. As the next step, $C_{\mathrm{p},s}-C_{T,rot}$ curve of the fully resolved CFD model was compared with the ducted turbines available in the literature. Among the variety of ducted configurations, the chosen geometries of NACA-5415 (D1) [60], NACA-5415 (D2) [61] and DonQi [16] resembled more closely the shell’s shape of the current study.

Figure 9 displays a higher peak value of $C_{\mathrm{p},\mathrm{s}}$in the present study (DU-93-W-210) is 8% higher with respect to NACA-5415 (D1). Contrariwise, $C_{\mathrm{p},\mathrm{s}}$value is approximately 14% and 18 % less than those of NACA-5415 (D2) and DonQi ducted models. The discrepancies in model results are thought to be inherently associated with multiple factors that include the tip clearance, function of airfoil thickness and camber line, bending of the duct cross-section and numerical techniques. In fact, a scenario of different numerical techniques seems intuitive herein to explain the inconsistencies in the comparative analysis of the results because the performance of the rotor inside the DonQi duct was computed by the actuator disc model (ADM). Similarly, the NACA-5415 (D1) and NACA-5415 (D2) calculations were carried out via a semi-analytical method. In this regard, the ADM technique mimics the rotatory analysis of the rotor inside the duct through surface forces instead of actual blade geometry. The other notable limitations are evenly distributed rotor forces and the disregarded kinematics of the rotor’s wake.

The FEM model of the blade was validated with the model results calculated from the Numerical Manufacturing And Design (NuMAD) [62] tool. Basically, NuMAD is a quite user-friendly, stand-alone, MATLAB-based pre-processor for the FEM solver. First of all, the original blade model was replicated in the NuMAD graphics user environment. The material properties along with the desired design thickness were assigned to each blade section and shear web accordingly. As a further step, another program PreComp [63] was called to compute the stiffness and inertial properties of the blade at user-specified stations. PreComp calculates the detailed structural properties of the blade using an improved laminate theory.

For the calculation of the natural frequencies and mode shapes, the main output files of the PreComp were utilized as input files for another program BModes [64]. BModes takes the inputs in terms of the blade geometry, rotor speed and material distribution across the span of the blade. A salient feature of BModes is that it calculated the desired solution data by employing a cantilevered boundary condition at the blade root. However, the available version of BModes can only handle calculating the isotropic material properties of the structure instead. Hence, the modal analysis of the actual blade FEM model has been separately performed again with isotropic material properties to achieve comparable results in terms of modal parameters. For validation purposes, the natural frequencies $\left({\omega }_i\right)$ of mode 1 and mode 2 obtained from both the solution solver are compared and found to be in close agreement as shown in Figure 10. The equation of FEM modal analysis is as follows;

$$(\left[K\right]-{{\omega }^2}_i\left[M\right]){\left\{\varnothing \right\}}_i=\left\{0\right\}$$

(21)

Here $\left[K\right]$ and $\left[M\right]$ are the stiffness and mass matrices, respectively. Here, ${\omega }_i$ is called the angular frequency of a structure and are related to natural frequency $(f_i={\omega }_i/2\pi )$ in cycle per second (Hz) and ${\left\{\varnothing \right\}}_i$ represents the mode shape corresponds to certain $f_i$.

5.2. Mode Shapes of Blade

Based on the modal characteristics [65], the free-vibration FEM procedure permits to evaluate the various mode shapes. A single FEM blade model rightly addresses the local solutions with non-rotating and fixed boundary conditions to the blade root. The first six vibration modes of the blade in flapwise and edgewise directions extracted from the FEM model are depicted in Figure 11 and Figure 12, respectively. As can be seen from mode shapes, the flapwise mode is dominating in the 1st, 3rd, 4th, 5th and 6th modes of free vibration. On the other hand, the edgewise mode is dominated by the 2nd mode and a slight glimpse appears in the 1st and 6th modes as well in Figure 12. These mode shape patterns are in line with the previous literature [23] which confirms the validity of the present FEM model.

5.3. Aerodynamic Forces and Augmentation

This section briefly describes the torque coefficient curves and the possible reasons behind the augmentation. Additionally, the aerodynamic forces on blade and shell structures are also discussed. An important curve of $C_{Tr}-V_0$ are replicated [20] here to illustrate the influence of wind speeds at a constant ${\lambda }_r=6.5$on the bare rotor and shell rotor as shown in Figure 13 (a). The predictions of rotor torque coefficient obtained from CFD indicates that the $C_{Tr}$ consistently raises with an increase in wind speed. As expected, the noticeable gain of ~66% higher torque coefficient is achieved with the same swept size of the rotor installed inside a shell under optimal conditions. The difference between the two examined configurations is quite profound at every given wind speed $(V_0=6-26\mathrm{m}/\mathrm{s})$. This $C_{Tr}$ the increment can be further associated with mass swallowed by the rotor plane enclosed in the shell configuration. Figure 13 (b) shows the mass flow rate $(\dot{m)}$ with varying $V_0$ values.

Despite the velocity disturbance at the rotor plane, using the ANSYS CFX built-in function accurately calculates the mass flow rate while accounting for these flow variations. This automatically integrates the local conditions over the specified surface, regardless of any disturbances caused by the rotor. One can notice a higher $\dot{m}$ in case of shell rotor than that of the bare rotor under rated conditions. It can be easily seen that $\dot{m}$ increases quite significantly with an increase in the given wind speed at constant TSR. This resulting increase in $\dot{m}$is mainly associated with shell geometry around rotor blades. Meanwhile, $\dot{m}$also holds a directly proportional relation $\dot{(m}={\rho }_h\pi R_b^2{V}_0)$ with the wind speed $V_0$.

An in-depth explanation requires the local velocity and pressure distributions. Figure 14 (a) presents a contour map of the normalized velocity in the yz-plane on a cross-section of the shell boundary. Under optimal conditions, the acceleration of airflow is less visible than the rotor tip. However, the relative magnitude of velocity appears to be comparatively high in the throat area and diminishes at the trailing edge. The disruptive action of the expansion ratio ${(A}_{exit}/A_{throat})$ in the divergent section slows down the axial component further downstream of the rotor. (b) demonstrates the significant effect of the normalized pressure on the shell. In this pressure contour, the magnitude of the normalized pressure is considerably high at stagnation points such as the leading edge of the shell’s profile and stabilizing wing. High pressure is evident in the entire frontal portion of the shell. Herein, it is reasonable to summarize that the converged section largely causes a faster movement of the flow. Meanwhile, abrupt changes that occur in a narrow and diverged section indicate a larger diffusion of momentum, eddying flows and pressure recovery.

and corresponding mass flow rate $\dot{(m})$.

To quantify the impact of tangential and normal forces with respect to the rotor blade, the rotor torque coefficient ${(C}_{Tr})$and rotor thrust coefficient ${(C}_{T,rot})$, are drawn across the blade span as shown in Figure 15 (a) and (b). The distribution of both $C_{Tr}$ and $C_{T,rot}$of a single blade exhibits a quite similar behavior that tends to rise at the inboard region $r/R_b=0.25$ and then continues to increase gradually between $r/R_b=0.4-0.8$. After that the maximum value of$C_{Tr}$ and $C_{T,rot}$ approximately happens to appear in the outboard region $r/R_b=0.90$ and finally diminishes afterward in both configurations. Overall, higher aerodynamic forces in terms of $C_{Tr}$ and $C_{T,rot}$ are achieved in the shell rotor compared to the bare rotor comes. Aerodynamically, the axial force ${(F}_a)$ experienced by the shell meaningfully contributes to the generate a higher thrust force in the case of the shell rotor. For this reason, a glimpse of the aerodynamic loadings on the shell is presented in terms of axial ${(F}_a)$, tangential ${(F}_t)$ and normal forces ${(F}_n)$ at the three different spans (08%, 50%, and 85%) of the shell as shown in Figure 16. It can be seen that notable transition occurs in all loadings at an angular position of 90°, 180°, 270° and 0°/360° because of the shadow effects of the four stabilizing wings at these locations.

and thrust coefficients ${(C}_{T,rot})$ in bare and shell rotors.

5.4. Aeroelastic Structural Responses

The blade total deformations under three operational conditions of bare configuration and three operational conditions of shell configuration are depicted in Figure 17 (a)–(f). One can clearly observe that all the maximum deformations occur on the blade tip region for all operational conditions regardless of the type of configuration. This is also evidently visible that the blade-tip deflection consistently increases as the wind speed increases from $V_0=12\mathrm{m}/\mathrm{s}$to $V_0=22\mathrm{m}/\mathrm{s}.$ However, the tip deflection is considerably higher in order of magnitude under the shell configuration than that of the bare configuration. There are obvious results that stem from the fact that the rotor in the shell configuration experiences higher wind loads due to the flow augmentation effect of the shell structure around the rotor plane. It is therefore intuitive that the accelerated flow owing to the shell configuration is inherently associated with the increasing deflection on the rotor blade. Hence, relatively increased loads on the blade surface are supported by the higher power-producing cycle of the rotor in the shell configuration. For more clarity, the results can be described as percentages to further support the discussion related to tip deflection.

It is evident that when wind speed increases from 14–18 m/s the percentage increase in deflection is around 45% in the bare rotor configuration. This further increases and approaches 51% when wind speed shifts from 18 m/s to 22 m/s. This means that the tip deflection increased on an average of 47% corresponding to each time increase of 4 m/s in the free stream wind velocity. A similar trend is quite dominating in the shell configuration where 61% and 56% increase in deflection occur in 14–18 m/s and 18–22 m/s, respectively. On the other hand, when comparing the deflection of both the rotor configurations against each other few more aspects of the structural simulation can be uncovered. In the operational condition of 14 m/s, the deflection remains just below 6% in the shell configuration. However, the rated operational condition of 18 m/s approximately exhibits an 18% increase in tip deflection. Similarly, an increase of 21% is observed on the blade tip at higher wind loads of 22 m/s of the proposed airborne wind energy system. The observed increase in deformation with rising aerodynamic forces, despite the absence of a direct relationship, can be attributed to the complex, nonlinear response of rotor blades to wind loads. This complex interaction between aerodynamic loading, material response, rotational inertia, and geometric irregularities introduce fluctuating loads that further intensify these deformation effects.

Figure 18 shows the blade stress analysis of both the bare and shell configurations. The magnitude of the maximum stress was calculated normal to all three x-axis, y-axis and z-axis directions corresponding to $V_0=14,18\mathrm{m}/\mathrm{s}$ and $22\mathrm{m}/\mathrm{s}$. The stress value is almost similar in both configurations in the x-axis at 14 m/s. One can visually see the stress magnitude is considerably high (>150%) in the z-axis direction in all cases of wind velocities. On the other hand, the equivalent stress exhibits a value very close to the z-axis stress magnitude since the z-axis direction contributed a much higher percentage than those of the other two axes. The profound effect of the stress concentration on the z-axis (flapwise deflection) is owing to the axial velocity in the z-direction during all the aerodynamic simulations. By considering the individual configuration, the stress increment is identified on an average of >45% and >55% corresponds to the bare and shell configuration.

Overall, the equivalent stress distribution on the rotor blade in shell configuration increases in all three cases of wind velocity 14 m/s, 18 m/s and 22 m/s up to 5%, 17% and 20%, respectively. In the worst case with a wind speed of 40 m/s, the maximum stress would reach by doubling the maximum stress values at 22 m/s on approximately 325 MPa based on linear relation. On the other hand, the tensile strength of most of the reinforced polymer composites ranges from 1000–2400 MPa. Generally, the material safety factor for wind turbine composite blades is around 2.2. In this case, the material safety factor is about 3.07, which is higher than the traditional benchmark of 2.2. This indicates the blade is unlikely to experience material failure under the given operational conditions.

As mentioned earlier, the goal of introducing a single shear web was to enhance the representative loads carrying capability of the blade. Figure 19 displays the stresses and strain induced in the shear web under wind loads of $V_0=18\mathrm{m}/\mathrm{s}$. The normalized values of stress and strain are drawn to determine how internal response could be mapped spanwise in the shell and bare rotor. The equivalent strain in terms of normal and shear stress components is estimated using the following relations;

$${\epsilon }_{eq}={\displaystyle \frac{1}{\sqrt{2}\left(1+{\nu }^{\text{'}}\right)}}{\left[{\left({\epsilon }_{xx}-{\epsilon }_{yy}\right)}^2+{\left({\epsilon }_{yy}-{\epsilon }_{zz}\right)}^2{+\left({\epsilon }_{zz}-{\epsilon }_{xx}\right)}^2+6\left({\epsilon }_{xy}^2+{\epsilon }_{yz}^2+{\epsilon }_{xz}^2\right)\right]}^{{\displaystyle \frac{1}{2}}}$$

(22)

Likewise, the von Mises stress or equivalent stress is written as:

$${\sigma }_{eq}={\displaystyle \frac{1}{\sqrt{2}}}{\left[{\left({\sigma }_{xx}-{\sigma }_{yy}\right)}^2+{\left({\sigma }_{yy}-{\sigma }_{zz}\right)}^2{+\left({\sigma z}_z-{\sigma }_{xx}\right)}^2+6\left({\sigma }_{xy}^2+{\sigma }_{yz}^2+{\sigma }_{xz}^2\right)\right]}^{{\displaystyle \frac{1}{2}}}$$

(23)

The stress and strain data obtained from Eqs. (22) and (23) are then normalized by dividing the maximum value of the stress and strain data and plotted against the spanwise length of the shear web. The curves of stress and strain in each configuration are visually identical but differ in terms of the observed peak. The predicted values of stress and strain on the normalized numerical scale exhibit an increment of 17% and 20% in the shell rotor as compared to the bare rotor, respectively. The maximum value of stress and strain occurs about 0.4 span of the shear web. Considering the importance of region $(r/R_b=0.25-0.4)$ of this peak value, the outer shape of the blade has already been designed by employing an S-821 airfoil having a thickness of $(t/c_l=24\mathrm{\%})$. The mounting evidence of structural analysis clearly justifies the use of a thick family of airfoils in modeling the blade. Additionally, the positioning of the S-821 airfoil was originally designed to facilitate the structural strength at the root region of the blade.

5.5. Operational Stability

Figure 20 shows the Campbell diagram of the blade model plotted from the initial modal analysis. A typical Campbell diagram provides an overall view of the vibration excitation that can occur on a rotating component during routine operation. The rotor blade was fixed at the root region while rotational velocity was applied clockwise in tabular form with an interval of 100 rpm having zero initial velocity. The first ten-mode shapes were subjected to generate the diagram in free vibration analysis. The natural frequencies of the system are plotted along the y-axis and the rotational velocity of the rotor on the x-axis. To meet structural and aerodynamic criteria, this frequency band indicates that operating the turbine above 611 rpm (red triangle) should be prohibited in all cases. It is pertinent to mention here that the design or rated operational speed of the rotor is 446 rpm and it is way below the computed threshold of the prohibited speed. In the worst-case scenario, the one possible situation the rotor blade might encounter is when it tends to achieve the TSR of 9.5 corresponding to the rotational velocity of 651 rpm at a rated wind speed of 18 m/s. However, it is already known from the aerodynamic analysis that the rotor must rotate near to rated rpm to optimize the power-producing cycle which restricts to avoid blade failure. Furthermore, the working operation of the rotor is reasonably stable in all the given rotational speeds other than that of the higher TSR close to the critical speed.

5.6. Deformation of Shell

To predict the deformation on the surface of the shell, the pressure loads were imported from the CFD solver at various load conditions. The pressure data from the CFD was perfectly mapped on the shell main as well as the four projected wings. These FEM simulations for the non-linear static structural analysis of the shell were simulated separately other than the rotor blades. For structural analysis of the shell, the sufficient deformation capacity of lower regions of the shell corresponding to the tethers joint must be ensured. Efforts were being made to model the tether’s effects through rigid support links to represent a physical consistent behavior of the shell structure. The simulations were carried out at three wind speeds of $V_0=14\mathrm{m}/\mathrm{s}$, $18\mathrm{m}/\mathrm{s}$ and $22\mathrm{m}/\mathrm{s},$ while the results of the total deformation can be seen in Figure 21 (a), (b) and (c). A relatively less deformation appears at the bottom half of the shell where the shell is constrained at four points using remote displacement boundary conditions. It is also visible that most of the deformation appears on the upper half of the shell.

The increment of wind velocity of $4\mathrm{m}/\mathrm{s}$ from $V_0=14\mathrm{m}/\mathrm{s}$ to $V_0=18\mathrm{m}/\mathrm{s}$ shows a percentage change of 37% in deformation. On the other hand, when velocity changes from $V_0=18\mathrm{m}/\mathrm{s}$ to $V_0=22\mathrm{m}/\mathrm{s}$ a percentage change of about 20.4% is observed. Hence, the structural elements by carrying both the membrane stiffness and bending stiffness, the shell structure can be considered equivalent to the semi-rigid structure. Keeping in view the deformation of the shell, the rotor was purposefully placed with a tip clearance of around 590 mm (0.59 m) from the main body of the shell as shown in Figure 21 (d). In this perspective, the calculated deformation is quite reasonable that further ensures the smooth and safe working of the rotor inside the shell. From a static structural analysis viewpoint, the results are quite realistic with the underlying assumptions of the rigid tethers and semi-rigid structure of the shell.

6. Conclusions

The current study has been carried out to analyze the aeroelastic response of the airborne wind turbine (AWT). Two vital components of the AWT namely the rotor and shell have been studied through an approach based on the one-way fluid-structure interaction (FSI) analysis. The fluid model was first established to apply a range of wind loads which was undertaken for each simulation case using a fully-resolved CFD model. These wind loads in terms of CFD pressure load are transferred to the FEM model to predict the response behavior of the AWT for airborne applications. The blade and shell under augmented loads were further investigated for non-linearities that could indicate structural deformations and stresses. The important conclusions are listed below.

• The augmentation effects of the shell around the rotor yielded a substantial gain of 66% in generated torque as compared to the bare rotor. Conversely, the impact of the shell exerted a higher aerodynamic force on the rotor blades in the shell configuration. Thus, a trade-off between performance enhancement and the structural integrity of AWT is mandatory while operating at a high altitude.
• The blade structure reinforced by a shear web showed larger tip deflections (18% increment) in shell configuration than that of the bare rotor under rated conditions. This implies that the maximum tip deflection (~ 32 mm) at a rated wind speed of $V_0=18\mathrm{m}/\mathrm{s}$, indicating that the blade is least vulnerable to experience a strike with the surrounding structure. The results further revealed that the tip deflection increased by an average of 47% corresponding to an increment of 4 m/s in the wind speed.
• The maximum stresses built in the orthotropic blade were examined to be 160 MPa, which are well below the material strength limits of the composite. Additionally, there was no discernible material failure observed during non-linear static structural analysis. The rotational speed of the rotor was found to be quite stable and reasonably safe at a rated condition of 446 rpm.
• The predicted result of the FEM of the shell also ensured the functional operability of the rotor within the design requirement by indicating a limited deformation (~ 22 mm) of the shell surface. Overall the one-way FSI analysis is the preferable choice to evaluate the aerodynamic loads and non-linear structural responses using a dedicated CFD module coupled with the FEM module.

The present study fully focused on the aeroelastic responses for the design solution of an airborne wind energy harvester to operate at higher altitudes. An interesting line of future research would be to employ more realistic features on the rotor and shell by taking into account composite layups and membrane elements using a two-way FSI approach.