**1. Introduction**

With the world's high demand of energy, and the limitation of the amount of fossil fuels, renewable energies have become a field of interest for many researchers. Wind energy is one of the most growing renewable sources of energy, in terms of usage and research topics. The global cumulative installed wind capacity has increased 2200% from the year 2001 to 2017 [1]. At present, horizontal axis wind turbines remain the dominant wind energy conversion technology. In the past decades, the trend was to increase the diameter of the rotor, since the power produced is proportional to the rotor area. Currently, the world's largest wind turbine has a 12 MW capacity, and a 107-meter blade length [2]. However, with this large size comes big challenges, such as the huge transportation and installation cost of extremely large wind turbines and the severe structure dynamic loads on the blades and the tower, as well as the need to develop each component, including the blade, bearing, generator, gearbox, etc. to be suitable for the large-scale single rotor turbines, which includes risk in cost and quality. Also, it includes a risk if a failure occurs; then the whole wind turbine will shut down and no power will be produced until the failed part is fixed or replaced [3].

The multi-rotor systems (MRS) is a technology with a long history that goes back to 1930, but it has fallen out of consideration for its structure complexity, while large scale single rotor wind turbines have become technically feasible [4]. However, with current advanced materials technology, the materials used to construct the rotors have a higher strength to weight ratio. With those advances in the materials technology, MRS is a promising alternative to large-scale wind turbines. The main advantages of MRS are the standardization of the wind turbine components, ease of transportation (since the rotors are of a small scale compared to large-scaled wind turbines), ease of installation, and the cost and the reliability of wind turbines, since multi rotors ensure that if there is a failure in one rotor, then the other rotors will still produce power. The major challenges in an MRS are the complexity of the supporting structure, the yawing system, and the aerodynamic interaction between rotors placed closely to each other.

As MRS has become an interesting field of research, researchers from many countries have made attempts to issue very interesting research points. Some attempts were interested in studying aerodynamics and the aerodynamic interaction between the rotors of an MRS. Experiments made by Goltenbott et al. [5] have shown that two and three di ffuser augmented rotor configurations can increase the power produced per rotor by 5% and 9% respectively, compared to a single rotor. Also, the computational fluid dynamic (CFD) simulations made by Chasapogiannis et al. [6] on a seven rotor system have shown a power increase of 3% per rotor. The coherence e ffect on the produced power and tower loads on a seven rotor MRS has been studied by Yoshida et al. [7]; Wind models with three di fferent coherences were used in the simulation and showed that larger coherence implies higher power production ye<sup>t</sup> increases the collective loads. MRS was also found to improve the wake recovery; the wakes were found to recover faster for MRS compared to a single-rotor configuration and showed a smaller turbulence intensity in the wake [8].

One of the leading research institutes showing grea<sup>t</sup> interest in MRS is the Technical University of Denmark (DTU). The DTU constructed a four rotor wind turbine at the Risø campus. They conducted experiments as well as simulations for the four rotor wind turbine, and both agreed that the interaction between the rotors improved the power performance by 1.8 ± 0.2%, which can increase the annual power production by 1.5 ± 0.2% [9].

Downscaling the design and cost of wind turbine rotors to replace a single large rotor with multiple smaller rotors has been done by Verma et al. [10]. A 5 x 1 MW multi-rotor turbine was compared to the National Renewable Energy Laboratory (NREL) 5MW single rotor turbine. The scaled down multi-rotor configuration has shown a 37% reduction in weight and a 25% reduction in cost compared to a single rotor producing the same amount of power. As an extension of Verma's work, Mate et al. [11] have designed a support structure for the five rotor configuration proposed in their colleague's work, in addition to other configurations that Mate proposed himself. A finite element approach was used for modeling the support structure. However, these simulations did not include a study of the aeroelastic behavior of either the blades or the support structure and the wind conditions.

Aeroelastic analyses for wind turbines are doable for single rotors using either a Computaional Fluid Dynamics (CFD) model like the Fluid Structure Interaction (FSI) made by Bazilevs et al. [12] or Halawa et al. [13] or using deterministic models. NREL's tool, FAST, is one of the most used aeroelastic tools for modeling wind turbines, it is based on models derived from fundamental theory of aerodynamics and structure analysis, which are more time-e fficient compared to CFD models [14].

So far there has been no research or open-source tools proposed that introduce aeroelastic analysis for the multi-rotor concept. In this work, the support structure of an MRS is being aero-elastically analyzed, so that structural problems can be addressed in further research studies. It is the first attempt to develop an in-house aeroelastic tool for an MRS support structure. In this work, the present tool is used to model a twin-rotor wind turbine, with two coplanar rotors placed on a T-shaped tower. This tool can be later extended to model support structures for di fferent configurations of MRS.

The theory used in this work is the Blade Element Momentum (BEM) theory, which is used to calculate the aerodynamic loads, and the virtual work method with a modal approach to calculate the structural deformations of blades and tower. Combining the two theories creates an aeroelastic interface between the blades and tower on one hand, and wind on the other.

The results of the present tool are verified by comparing results of a single rotor wind turbine, to NREL's FAST results of the same turbine model. Then, the results for the twin-rotor configuration are introduced.

## **2. Mathematical Model**

## *2.1. Aerodynamic Model*

Unsteady BEM was used in this work to estimate the aerodynamic loads tangential to and normal on each section along the blade. The blade coordinate system was used for the governing equations, as shown in Figure 1.

**Figure 1.** Blade Coordinate System [15].

Initially, the relative wind velocity components on each blade section, in the rotor plane, tangential to the blade width (*yB*-axis), and normal to it (*xB*-axis) as shown in Figure 2 are as follows:

$$V\_{rel,x\_B} = V\_{0,x\_B} + w\_{xy\_B} \tag{1}$$

$$V\_{\text{rel},y\_B} = V\_{0,y\_B} - \alpha r + \mathfrak{w}\_{\text{yB}} \tag{2}$$

where *V*0 is the inflow velocity, ω is the rotational speed of the rotor, *r* is the blade section position, and *wxB*, *wyB* are the induced velocities in *xB* and *yB* directions.

**Figure 2.** Velocity Triangle on a Blade Section [16].

Prandtl's tip loss factor and Glauert correction are used in the calculation of induced velocities; then, depending on the relative velocity components on each blade section, for every time step, the coefficients of lift and drag forces can be interpolated from the blade airfoils' data tables. Next, lift and drag forces can be calculated for each section, and hence the aerodynamic loads in the plane of rotation (*LyB,aero*) and normal to it (*LxB,aero*) can also be calculated. Dynamic wake, dynamic stall, and yaw misalignment effects were ignored in this study. With the tangential and normal load distributions known, the rotor aerodynamic thrust, torque, and power can be calculated [16].

## *2.2. Structure Model*

The principle of virtual work was used to calculate the structure dynamics parameters. This principle is a method that helps to set up the correct matrices for a discretized mechanical

system, such as Newton's second law shown in Equation (3). This method is well suited for a chained multi-body system like wind turbines:

$$\mathbf{M}\ddot{\mathbf{x}} + \mathbf{C}\dot{\mathbf{x}} + \mathbf{K}\mathbf{x} = \mathbf{F}\_{\mathcal{X}} \tag{3}$$

where; *M* is the mass matrix, *C* is the damping matrix, *K* is the sti ffness matrix, and *Fg* is the generalized forces array. *x*, . *x*, and .. *x*, are the generalized coordinates for the three modes used in the model, their first, and second time derivatives, respectively.

The degrees of freedom used in this model are the uncoupled modal shapes for both the blades and the tower. For the blades, the first and second flap-wise, and the first edge-wise modes are used, denoted by *uf* 1, *uf* 2, and *ue*<sup>1</sup> respectively. For the tower, the first and second fore-aft, and the first side-side modes are used, denoted by *ufa*1, *ufa*2, and *uss*1. These modes were chosen because the orthogonality constraint between the eigen modes turns the mass matrix into a diagonal matrix, and hence into the damping and sti ffness matrices that are dependent on the mass matrix. This results in uncoupled di fferential equations that can be solved by a fourth order Runge-Kutta numerical technique.

Solving this system of equations results in the values of *x*, . *x*, and .. *x*, on which the deformation of either the blade or the tower is assumed to be dependent, as a linear combination of the three modes. For instance, in the case of the blade, with a sti ff tower, the deformation distribution along the blade section position (*r*), denoted by *uxB* in the flap-wise direction and *uyB* in the edge-wise direction, will be as follows:

$$\mu\_{\rm xB}(r) = \mathbf{x}\_1 \boldsymbol{u}\_{\rm xB}^{f1}(r) + \mathbf{x}\_2 \boldsymbol{u}\_{\rm xB}^{e1}(r) + \mathbf{x}\_3 \boldsymbol{u}\_{\rm xB}^{f2}(r) \tag{4}$$

$$
\mu\_{yB}(r) = \mathbf{x}\_1 \boldsymbol{\mu}\_{yB}^{f1}(r) + \mathbf{x}\_2 \boldsymbol{\mu}\_{yB}^{e1}(r) + \mathbf{x}\_3 \boldsymbol{\mu}\_{yB}^{f2}(r). \tag{5}
$$

The velocity and acceleration distribution along the blade are calculated the same way, except for using . *x*, and .. *x* instead of *x* [17].

## *2.3. Aeroelastic Coupling*

The deformation of the blade, together with the velocity and acceleration of its vibration, will result in change in the loads and hence structural deformation in the next time-step. The relative wind velocity components on the blade sections can now be updated to include the blade vibrations, with the blade velocity distribution always opposing the wind direction. Equations (1) and (2) can be updated as follows:

$$V\_{\text{ref}, \text{xB}} = V\_{0, \text{xB}} + \text{w}\_{\text{xB}} - \dot{\mathbf{u}}\_{\text{xB}} \tag{6}$$

$$V\_{rel,yB} = V\_{0,yB} - \alpha r + w\_{yB} - \dot{u}\_{yB} \tag{7}$$

where . *uxB*, and . *uyB* are the blade sections' vibrational velocities in the flap-wise and edge-wise directions.

The loads on the blade sections are also updated to include gravitational loads and inertia loads due to the blade vibrations. The total load distribution along the blade sections in the normal and tangential directions to the plane of rotation will be as follows:

$$L\_{\rm xB}(r) = L\_{\rm xB\,\mu\sigma\nu}(r) - m'(r)\ddot{u}\_{\rm xB}(r) + m'(r)g\sin(\theta\_l + \sin\theta\_\varepsilon)\cos\theta\_A \tag{8}$$

$$L\_{yB}(r) = L\_{yB, \text{arro}}(r) - m'(r)\bar{u}\_{yB}(r) + m' \sin \theta\_A \tag{9}$$

where *m* is the mass density distribution along the blade length, .. *uxB* and .. *uyB* are the blade sections' vibrational accelerations, θ*t* is the tilt angle, θ*c* is the cone angle, and θ*A* is the azimuth angle. The loads in the direction of the blade length are neglected.

The updated relative wind velocity and loads are used for the next time step, guaranteeing that both the aerodynamics and structure models a ffect each other to create an aeroelastic model.

The tower is modelled independently from the blades, and the rotor loads are transmitted to the tower, including rotor thrust, weight, and torque, taking into consideration the aeroelastic behaviour of those loads. Also, the aerodynamic load on the tower itself is calculated, considering the tower vibrations in the relative wind speed to the tower, and the tower inertial loads.

#### **3. Verification of the Single Rotor Model**

The single rotor wind turbine simulated for validation is the NREL 5MW wind turbine. In order to ensure the reliability of the developed aeroelastic tool, the results of aerodynamic analysis, dynamic structural analysis, and aeroelastic analysis are compared to results of the NREL 5MW results in the definition report of the wind turbine [18] and FAST simulation for the aeroelastic case.

## *3.1. Model Outline*

The geometric and material properties of the turbine blades and tower are fully defined in the NREL definition report. General rotor specifications are shown below:


Figure 3 shows a CAD model for the turbine blade.

**Figure 3.** NREL 5MW wind turbine blade [19].

The blade and tower structural specifications are shown in Figure 4.

**Figure 4.** *Cont*.

**Figure 4.** Blade structural properties: (**a**) Blade mass density distribution; (**b**) Blade stiffness distribution; (**c**) Tower mass density distribution; and (**d**) Tower stiffness distribution.

## *3.2. Rotor Aerodynamics*

The NREL 5MW wind turbine is a pitch-controlled turbine. Over the rated wind speed of 11.4 m/s, the pitch control is applied to maintain the rated power output and rotor speed. So far, pitch control has not been applied in the present tool. Accordingly, the full power region cannot be simulated in the current version of the tool, but it will be added later. The BEM in the developed tool is run for wind speeds starting from the cut-in speed of 3 m/s, to the rated speed of 11.4 m/s. The rotor power, thrust, and torque are compared. The results are shown in Figure 5.

**Figure 5.** Steady state responses: comparison between results of the developed tool to NREL reference file.

Figure 5 shows good agreemen<sup>t</sup> between the calculated steady state values of the power, torque, and thrust of the developed tool, compared to the turbine definition document's results. This agreemen<sup>t</sup> proves that the first part of the tool concerned with calculating the aerodynamic responses, without considering the aeroelastic behavior of the blades, is reliable, and the next step would be the dynamic responses of the blades and tower.

## *3.3. Modal Analysis*

In this section, the free vibration modes of the blades and tower are investigated. The structural model is based on the uncoupled modes of vibrations, and it is important that the mode shapes and natural frequencies will be correct, before proceeding to the aeroelastic coupling between the aerodynamics and structural behavior of both the blade and the tower.

The NREL tool Modes [20] was used to calculate the uncoupled mode shapes and natural frequencies of the blades and the tower. Mass and stiffnesses distribution along the blade and tower are available in the NREL 5MW definition document and were used as input for Modes. The mode shapes and natural frequencies are also available in the definition document and are shown for comparison with the results of simulation.
