Aerostructural optimization of a morphing wing for airborne wind energy applications

The morphing wing concept investigated in this study, shown in figure 2, is developed from the compliant structure introduced by Molinari et al [7–10], and attains variations of the local lift coefficient (c_l) by means of chordwise morphing. Differently from the concept introduced by Molinari et al, in the proposed design the mechanical energy required to morph the structure is provided by electromechanical linear actuators instead of active piezoelectric elements. These actuators are chosen due to their higher static load-carrying capability (capacity of maintaining a fixed position without using energy), their energy density, and their favorable force-displacement characteristics [11, 12]. The deformation introduced by the actuators is guided by a numerically optimized distributed compliance rib. The actuators and ribs are organized in rib-actuator sets, consisting of a rib pair and a single actuator. To achieve the desired shape changes, the possibility to elongate part of the compliant skin is crucial, therefore, a corrugated skin with highly anisotropic stiffness is introduced on the lower side of the wing. This corrugated skin allows the desired in-plane elongation to occur, while contributing to withstand the large bending loads experienced by the structure [13–15].

Zoom In Zoom Out Reset image size

Overall, the wing consists of nine evenly distributed rib-actuator sets per half span. A stiff D-nose wing box carries the majority of the bending, shear and torsional loads. The material used for the wing and for the internal compliant rib is carbon fiber reinforced plastic (CFRP) owing to its anisotropy and favorable lightweight characteristics. The layup is concurrently optimized with the structure and aerodynamic shape, and is presented in section 3.1.

Owing to the high wing loading and to the compliance of the entire wing—and particularly of its compliant section—it is necessary to consider aeroelastic interactions to accurately predict the aerodynamic and structural behavior and to assess the resulting performance. The response of the wing on actuation inputs and aerodynamic loads is thus assessed by using a two-way weakly coupled 3D static aeroelastic analysis tool [7]. The aerodynamic behavior of the wing is simulated with a 3D panel method [16, 17] and a nonlinear extended lifting-line technique [17, 18], the latter using XFOIL [19, 20] to calculate the sectionwise 2D aerodynamic characteristics. The structural behavior is assessed by means of a 3D finite elements model. The wing skin and internal structure are modeled by plate elements, the stringers and actuators are modeled, respectively, as beam and as thermally-expanding rod elements, and the corrugated skin is modeled using flat elements and a substitute plate model [21, 22]. The aeroelastic analysis tool, validated through high fidelity numerical simulations and wind tunnel tests [7, 10], iteratively assesses the aerodynamic characteristics (lift and drag) and the structural displacements caused by actuation inputs and aerodynamic loads.

In this study, analytical models are used to assess the energy extraction capabilities of the AWE wing. The power P produced by an AWE system in the traction phase at optimal reel-out speed (the reel-in phase is neglected in this investigation, as it is strongly dependent upon the followed trajectory) can be modeled as a function of the lift and effective drag coefficient c_L and c_D, the wing area $A$, wind speed V_W, and angle between tether and ground $\gamma$ [1, ch 1]:

$$\begin{eqnarray}&&P=\displaystyle \frac{2}{27}\,\rho \,A\,\displaystyle \frac{{c}_{{\rm{L}}}^{3}}{{c}_{{\rm{D}}}^{2}}\,{V}_{{\rm{W}}}^{3}\,{\cos }^{3}(\gamma ).\end{eqnarray} \tag{ 1 }$$

The effective drag coefficient consists of the drag of the wing and also of the drag of the tether. Houska et al [23] suggest the following approximation:

$$\begin{eqnarray}&&{c}_{{\rm{D}}}={c}_{{\rm{D}},{\rm{wing}}}+\displaystyle \frac{1}{4}\displaystyle \frac{{l}_{{\rm{tether}}}\cdot {d}_{{\rm{tether}}}}{{A}_{{\rm{wing}}}}{c}_{{\rm{D}},{\rm{tether}}}.\end{eqnarray} \tag{ 2 }$$

At a given wind speed, the flight speed $V$ at optimal reel-out speed (${V}_{{\rm{ro}},\mathrm{op}{\rm{t}}}$ = $1/3$V_W) can be calculated as [1, ch 1]:

$$\begin{eqnarray}&&V=\displaystyle \frac{2}{3}\,\displaystyle \frac{{c}_{{\rm{L}}}}{{c}_{{\rm{D}}}}\,{V}_{{\rm{W}}}\cos (\gamma ).\end{eqnarray} \tag{ 3 }$$

The influence of the mass of the aircraft ${m}_{{\rm{A}}}={m}_{\mathrm{Wing}}\,+{m}_{\mathrm{Actuator}}+{m}_{\mathrm{Fuselage}}$ on the power output of the AWE system can be considered in the analysis by increasing the tether elevation angle by ${\rm{\Delta }}{\gamma }_{{\rm{mass}}}$ (cosine losses), and by rolling the aircraft into the circular trajectory (of radius r) to counteract centrifugal forces [1, ch 3].

$$\begin{eqnarray}&&P=\displaystyle \frac{2}{27}\,\rho \,A\,\displaystyle \frac{{{c}}_{{\rm{L}}}^{3}}{{c}_{{\rm{D}}}^{2}}\,{V}_{{\rm{W}}}^{3}\,{\cos }^{3}(\gamma +{\rm{\Delta }}{\gamma }_{\mathrm{mass}})\cdot {\left(1-{\left(\displaystyle \frac{2{m}_{{\rm{A}}}}{A\rho {c}_{{\rm{L}}}r}\right)}^{2}\right)}^{\tfrac{3}{2}}.\end{eqnarray} \tag{ 4 }$$

The mass of the wing ${m}_{\mathrm{Wing}}$ is obtained from the FE model, and the mass of the fuselage is set to ${m}_{\mathrm{Fuselage}}=10\,\mathrm{kg}$, to account for the mass of the wing-to-tail connection, the tail empennage, and the electronics. The mass of the actuators is determined as a function of the required power using the relation ${m}_{{\rm{Actuator}}}=0.1317\cdot {({P}_{{\rm{Actuator}}})}^{0.8251}$ [11]. The required actuators power output can be calculated by the force acting on the actuator, assessed within the aeroelastic simulation, the maximum stroke of the actuator, and the maximum time required to fully activate the actuator, which is set to ${t}_{\mathrm{Actuator}}=0.5\,{\rm{s}}$.

From formulas (1)–(4), the maximum power output of an AWE wing with known aerodynamic properties (lift and drag coefficient) can be calculated as a function of the wind speed. To predict the average annual power production in the following optimization, the performance of each individual wing is assessed at four flight speeds and corresponding angles of attack, where V₃ is the design flight speed, V₂ and V₁ are, respectively, the intermediate and the slow flight speed, and V₄ is the design speed plus gust speed assuming a horizontal gust¹ .

$$\begin{eqnarray}{V}_{1} & = & 0.5{V}_{3},\\ {V}_{2} & = & 0.75{V}_{3},\\ {V}_{4} & = & {V}_{3}+{\rm{\Delta }}{V}_{\mathrm{gust}}.\end{eqnarray} \tag{ 5 }$$

The aerodynamic coefficients (c_L and c_D) at the other flight speeds are determined by linear interpolation and constant extrapolation from V₁ to ${V}_{\mathrm{stall}}$. The power can then be assessed for flight velocities ranging from ${V}_{\mathrm{stall}}$ to ${V}_{3}$. At ${V}_{3}$, the maximum tether force and wing loading ${F}_{\mathrm{aero},\max }$ is reached. After V₃ is reached, the reel-out speed is increased to limit the force acting on the wing, while still extracting power. Apart from limiting tether and wing loading, the power of an AWE system is also constrained by the nominal generator power. However, as the focus of this study lies on the wing and not on the overall system, no constraint on the maximum power is enforced. As a result, the power at velocities exceeding ${V}_{3}$ is calculated by:

$$\begin{eqnarray}&&P(V\gt {V}_{3})={F}_{\mathrm{aero},\max }\cdot {V}_{{\rm{ro}}}={P}_{V3}\displaystyle \frac{{V}_{{\rm{W}}}(V)}{{V}_{{\rm{W}}}({V}_{3})}.\end{eqnarray} \tag{ 6 }$$

With this set of formulas, the average power output of the AWE system can be calculated by considering the wind distribution probability f_W as:

$$\begin{eqnarray}&&\bar{P}=\sum {f}_{{\rm{W}}}\cdot P({V}_{{\rm{W}}}).\end{eqnarray} \tag{ 7 }$$

At V₄, the morphing capability at high loads is investigated. The compliant section of the wing is fully deflected upwards over the whole span and the ability to alleviate gust loads is studied. The time to fully deflect the actuator is set to ${t}_{\mathrm{Actuator}}=0.5\,{\rm{s}}$, shorter than the rise time of an extreme gust [24]. The gust velocity is calculated after the IEC 61400-1 standard [25]. The detailed integration into the optimization procedure is explained in section 3.2.

The wind distribution probability used in this study is the average wind distribution at the locations of all conventional offshore wind parks with a minimum capacity of $200\,\mathrm{MW}$ at altitudes of 50, 100, and $200\,{\rm{m}}$ [26, 27]. The rationale for the choice of these locations is that AWE systems placed offshore show greater potential due to the stronger and more uniform winds compared to systems located onshore. Additionally, visual and acoustic impacts of offshore systems are smaller compared to their inland counterparts [28].

The calculated average offshore wind distribution is shown in figure 3.

Zoom In Zoom Out Reset image size

3.1.1. Aerodynamic shape

3.1.2. Structure

The compliant rib structure, shown schematically in figure 4, is described by a Voronoi-based parametrization [7, 31]. The main advantages of this parametrization are the small number of variables to describe the structure, its stability properties (small variable changes lead to small shape changes), and the fact that the obtained structure has no unconnected nodes and segments. In addition to the coordinates of the Voronoi-sites, the thickness of the segments is considered in the optimization, as well as the width of the rib. Having twelve Voronoi-sites, a total of 37 design variables per rib are used. The actuator is defined by four design variables per rib: three describe its position (vertical position at spar y_AF, vertical and horizontal position at trailing edge y_AR and x_AR), and one defines the maximum stroke of the actuator. At the design flight speed, the actuators are unactuated. At both speeds ${V}_{1}$ and ${V}_{2}$, the actuators are activated at an optimized level (p₁ and p₂) ranging from −1 to 1, and at gust speed, the actuators are fully actuated at a uniform level ${p}_{4}=-1$. Additionally, the position of the spar and the position of the corrugation is defined as a design variable (horizontal position spar x_S, horizontal position corrugation x_C, length corrugation ${l}_{{\rm{C}}}=1/4({x}_{{\rm{C}}}-{x}_{{\rm{S}}})$), and the position of the outer tether attachment point (TAP) can be varied along the spanwise direction.

Zoom In Zoom Out Reset image size

3.1.3. Wing skin and spar layup

The wing skin and spar composite layup ply thickness is concurrently optimized with the wing structure and the aerodynamic shape. The wing skin is split in four sections along its span, and in five sections in chordwise direction, leading to a total of 20 skin-sections, shown in figure 5. The five sections in chordwise direction are the wing box, the top intermediate and top rear section, the bottom section, and the corrugation. The wing box section has a [0°, 90°, ± 45°, 90°, 0°] layup, of which the thickness of the individual CFRP Unidirectional (UD) plies are considered as design variables² . The thickness of the composite layers of the top rear, top intermediate, and bottom segments are parametrized with two design variables, and feature a [0°, 90°, 90°, 0°] CFRP UD layup. The skin at the location of each compliant rib is reinforced with an additional 90° CFRP UD ply, resulting in one additional design variable per rib-actuator set. Likewise, the internal compliant rib structure consists of 90° CFRP UD. The corrugation layup is constant, whereas the upper and lower stringer thickness and cross-section and the spar thickness can be varied, resulting in the last five design variables.

Zoom In Zoom Out Reset image size

This parametrization leads to a total of 604 design variables, summarized in table 1.

Table 1.  Design variables.

Group	#DV	Design variables
Aerodynamic	161	15 × 9 airfoil, 1 span
shape		9 chord, 8 sweep, 8 twist
Structure	356	36 × 9 Voronoi, 9 corrugation
		9 rib width, 9 spar, 4 stringer
		1 tether attachment point
Actuator	36	4 × 9 actuation
Layup	44	36 skin, 4 rib, 4 spar
General	7	1 velocity, 1 tether length
		2 actuation gains (at V1 and V2)
		3 angles of attack (at V1, V2, V3)
Total	604

In this study, parameters related to different disciplines (structure, aerodynamics) are optimized concurrently. This increases the complexity of the optimization compared to a sequential approach, but will result in the identification of a better-performing individual, due to the possibility of assessing and exploiting interdisciplinary interactions [32]. The objective of the optimization is to maximize the average annual power production per wing area of a morphing AWE system. The goal of maximizing the power per wing area instead of maximizing the overall power prevents the optimizer from simply increasing the area of the wing to improve the objective function. Increases in wing area are still beneficial, as the relative tether drag decreases with increasing wing area. The analytical models discussed in section 2.3 are used to evaluate the power output of the system.

The optimization includes two constraints. The first aims to ensure a wing design with sufficient morphing capacity, able to alleviate extreme gust loads [25]. This constraint is expressed by imposing that the lift generated when flying at gust speed ${V}_{4}={V}_{3}+{\rm{\Delta }}{V}_{\mathrm{gust}}$ with the wing fully morphed to reduce the lift has to be smaller or equal than the lift at the design speed with the actuators inactive, resulting in the constraint:

$$\begin{eqnarray}&&{c}_{{\rm{L}},3}\cdot {V}_{3}^{2}\geqslant {c}_{{\rm{L}},{\rm{gust}}}\cdot {({V}_{3}+{\rm{\Delta }}{V}_{{\rm{gust}}})}^{2}.\end{eqnarray} \tag{ 8 }$$

The second constraint requires the final design to withstand the aeroelastic loads without incurring buckling.

$$\begin{eqnarray}&&{{MS}}_{\mathrm{buckling}}=\displaystyle \frac{{F}_{\mathrm{critical}}}{{F}_{\mathrm{applied}}}-1\geqslant 0.\end{eqnarray} \tag{ 9 }$$

Designs satisfying this constraint were found to simultaneously fulfill strength criteria [9], hence no additional structural constraints are imposed. Dynamic aeroelastic instabilities of the wing are not investigated in the optimization, as preliminary studies showed the flutter speed of the investigated wing to be outside of the flight envelope. To confirm this, the flutter speed of the final wing has been assessed in a subsequent step, confirming the validity of the assumptions and the findings of the preliminary study³ .

The optimization problem can therefore be formulated as follows:

$$\begin{eqnarray}\mathrm{maximize} & & {\rm{f}}({\bf{x}})=\displaystyle \frac{\bar{P}({\bf{x}})}{A}\qquad \qquad \quad {\bf{x}}\in {{\mathbb{R}}}^{n}\\ \mathrm{subject}\ \mathrm{to} & & {\bf{x}}{}_{\mathrm{low}}\leqslant {\bf{x}}\leqslant {{\bf{x}}}_{{\rm{up}}}\qquad \qquad {\bf{x}}{}_{\mathrm{low}},{{\bf{x}}}_{{\rm{up}}}\in {{\mathbb{R}}}^{n}\\ & & {\bf{g}}({\bf{x}})\leqslant 0\\ & & {\bf{g}}({\bf{x}})=\left\{\begin{array}{c}\tfrac{{c}_{{\rm{L}},{\rm{gust}}}\cdot {({V}_{3}+{\rm{\Delta }}{V}_{{\rm{gust}}})}^{2}}{{c}_{{\rm{L}},3}\,\cdot \,{V}_{3}^{2}}-1\\ -{{MS}}_{\mathrm{buckling},{V}3}\\ -{{MS}}_{{\rm{buckling}},{\rm{gust}}}\end{array}\right\}.\end{eqnarray} \tag{ 10 }$$

The objective function f(x) is defined as the average annual power production $\bar{P}$ per wing area. The inequality constraints are related to the morphing capability at gust loads and the buckling load coefficients at V₃ and ${V}_{3}+{\rm{\Delta }}{V}_{\mathrm{gust}}$.

The optimization is performed using the algorithm covariance matrix adaptation-evolution strategy (CMA-ES) [34], previously used by Molinari et al [7], as it satisfies the requirements of being able to treat nonsmooth objective functions and of dealing with a large number of design variables.

Due to the complexity of the problem and the limitations imposed by the constraints, identifying a suitable initial guess is not trivial. To this purpose, the initial wing design is identified through two sequential pre-optimizations. First, a purely aerodynamic optimization, relying on XFOIL, is conducted to identify the initial airfoil shape of the wing. The airfoil shape, common to all wing sections, is parametrized through the Bezier-PARSEC technique previously discussed. Aeroelastic effects are neglected and the airfoil shape is optimized in terms of average annual power production per wing area using the analytical models and wind distributions defined in section 2.3. 3D aerodynamic effects are modeled through analytical approaches to predict the induced drag assuming a constant Oswald factor equal to e = 0.7. With this pre-optimization, first insights about optimal airfoil shapes for AWE wings are gathered.

A preliminary sizing for the various structural components, in order to guarantee the load-carrying capacity and perceivable morphability for the initial guess, is obtained by concurrently optimizing the initial wing structure and composite layup. This second pre-optimization uses the same analysis and optimization tools as the main optimization, however, the objective is to minimize weight whilst achieving the desired morphing capability.

The performance of the initial wing defined by these two pre-optimizations is referred to as 'optimization 0' and is included in the comparison presented in the following section.

In optimization 1, only the 'General' design variables are optimized. The main difference with respect to optimization 0 is that the wing is optimized in terms of power per area and not in terms of weight and morphing capability. Additionally, the angle of attack is optimized at each considered speed, and the wing is actuated at V₁ and V₂, leading to an increase in weighted average annual power production per wing area of 3%.

In optimization 2, the 'General' and 'Aerodynamic shape' variables are optimized simultaneously. In figure 9, the initial and optimized wing planforms are shown. It can be seen that the planform of the wing is drastically changed. The wing chord length is increased at the root and decreased at the tip of the wing, leading to a decreased induced drag and reduced bending moment, ultimately decreasing the structural weight. The total wing area is increased from 5.1 to $5.5\,{{\rm{m}}}^{2}$, decreasing the relative tether drag as previously discussed. In figure 10, the optimized airfoil geometry is shown. Compared to the initial airfoil (optimization 0 and optimization 1), identified through the pre-optimization, the airfoils at the different sections show distinctive properties, being even more strongly cambered and thick at the root (high lift, supported by inner tether), slightly less cambered and thicker around the outer TAP (high bending moment due to outer TAP), and even less cambered towards the tip (to reduce induced drag and bending moment).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As a result of the optimized planform and airfoil geometry, the power harvesting factor is increased, as shown in figure 7(b), resulting in a 35% higher weighted average annual power production per wing area.

Figure 8(a) shows that the structural weight of the wing increases between optimization 1 and optimization 2, due to the larger wing area of the latter individual. The increased weight normally leads to a performance decrease. However, due to the thicker cross-sections of the wing compared to the initial design, the bending stiffness is increased and the flight envelope is enlarged. In detail, it leads to an increase in the design flight speed V₃ (shown in figure 8(b)) from ${V}_{3,\mathrm{opt}1}=91\,{\rm{m}}\,{{\rm{s}}}^{-1}$ at the first stage of the optimization to ${V}_{3,\mathrm{opt}2}=99\,{\rm{m}}\,{{\rm{s}}}^{-1}$ at the second stage, leading to the higher power production.

In the last optimization, the design space is further enlarged from 168 to 604 design variables, accounting for the wing structure and composite layup, and leading to an additional improvement of the objective function by 8%. The final internal structure changes slightly with respect to the initial pre-optimized structure. The airfoil shapes are very similar to the previous optimized shapes, being slightly thicker at the root to additionally increase the bending stiffness. The layup thickness at the different sections is shown in figure 11. It can be seen that the maximum thickness is reached at the wing box of section 1. The thickness at section 2 is slightly decreased compared to section 1. At section 3, the skin needs to be reinforced due to the outer TAP, and thus the thickness increases again.

Zoom In Zoom Out Reset image size

The ply thickness is optimized as a continuous design variable, possibly leading to infeasible designs. Therefore, depending on the manufacturing method and available material, the final ply thickness needs to be adjusted in a post-processing step.