Sage Journals: Discover world-class research

Abstract

Camber morphing aerofoils have the potential to significantly improve the efficiency of fixed and rotary wing aircraft by providing significant lift control authority to a wing, at a lower drag penalty than traditional plain flaps. A rapid, mesh-independent and two-dimensional analytical model of the fish bone active camber concept is presented. Existing structural models of this concept are one-dimensional and isotropic and therefore unable to capture either material anisotropy or spanwise variations in loading/deformation. The proposed model addresses these shortcomings by being able to analyse composite laminates and solve for static two-dimensional displacement fields. Kirchhoff–Love plate theory, along with the Rayleigh–Ritz method, are used to capture the complex and variable stiffness nature of the fish bone active camber concept in a single system of linear equations. Results show errors between 0.5% and 8% for static deflections under representative uniform pressure loadings and applied actuation moments (except when transverse shear exists), compared to finite element method. The robustness, mesh-independence and analytical nature of this model, combined with a modular, parameter-driven geometry definition, facilitate a fast and automated analysis of a wide range of fish bone active camber concept configurations. This analytical model is therefore a powerful tool for use in trade studies, fluid–structure interaction and design optimisation.

Keywords

Morphing wings variable camber plate theory Rayleigh–Ritz penalty method stiffness discontinuities orthogonal polynomials composite plates

Introduction

Due to the nature of aircraft flight operations, the lift coefficient required to sustain flight varies depending on the flight stage. Factors such as aircraft weight (variable during flight), altitude, atmospheric conditions and manoeuvres determine the required lift coefficient at a specific flight stage. There are currently two main approaches to generate these variations in lift coefficient: changing the angle of attack of the whole aircraft or by varying camber distribution. The latter approach can be achieved with traditional control surfaces and high lift devices (i.e. ailerons, elevator, trailing edge flaps), which consist of hinged panels attached to the central wing. Although the deflection of these hinged panels causes an effective variation in lift force, the sharp change in aerofoil section produces a significant drag penalty and noise. While the drag associated with manoeuvres does not typically make a significant impact on fuel burn over a mission, the significant drag penalty of traditional plain flaps makes them an unattractive option for more ambitious attempts to actively and continuously control the camber of the entire wing as a means of directly affecting the efficiency of the aircraft.

Camber morphing allows aerofoil properties to be modified with a much lower drag penalty, as the changes in camber occur in a smooth and continuous way. An ideal solution to this issue is to generate a smooth and continuous change in camber distribution. These not only make camber morphing devices attractive to replace trailing edge hinged plain flaps, but also to be actively used during flight to optimise the spanwise lift distribution.

As a consequence, there has been a significant motivation to develop wings that are capable of continuously adapting their geometries to different flight scenarios. These types of wings are known as morphing wings, and they are used as a method to improve the overall aerodynamic performance of aircraft. These research efforts are summarised by Thill et al. (2008), Barbarino et al. (2011) and Valasek (2012); specifically, all authors highlight that one of the most commonly pursued and researched method to optimise aerodynamic performance, at several flight stages, is active changes in the camber distribution of aerofoils.

One of the first designs that targeted variable camber without surface discontinuities was presented in the NASA’s F-111 Mission Adaptive Wing Aircraft, where a variable camber wing is implemented to optimise aerodynamic performance at different flight scenarios (Larson, 1986). Another variable camber concept is the DARPA Smart Wing (Kudva, 2004), where a centre laminate with honeycomb core and silicone skin give the trailing edge the flexibility needed to create changes in camber (Figure 1).

Figure 1.

DARPA smart wing variable camber concept. Reproduced from Kudva (2004), with author’s permission.

Some research groups have studied variable camber using leading edge devices: Vasista et al. (2017) designed and manufactured a compliant ‘droop-nose’ morphing leading edge device using superelastic materials. Furthermore, combining leading and trailing edge morphing within the same wing has also been investigated: Kota et al. (2003) performed a study on the feasibility of variable camber combining both leading and trailing edge morphing devices using both voice coil motors and piezoelectric actuators, while De Gaspari et al. (2014) designed, manufactured and tested a variable camber wing with morphing leading and trailing edge devices based on compliant-based mechanisms. Also, Monner et al. (2009) proposed a composite monolithic compliant morphing trailing edge device; however, no actuation mechanism was selected. Moreover, De Gaspari and Ricci (2010) performed a two-level optimisation that first obtained the ‘best’ aerofoil configuration for a given condition, followed by obtaining the best internal structural configuration capable of achieving such aerofoil geometry.

Furthermore, other concepts achieved variable camber by focusing the camber changes on the trailing edge, for example, by embedding actuators within the wing skin (Bilgen et al., 2011; Molinari et al., 2011), exploiting bistability in non-symmetric fibre-reinforced composite laminates (Daynes et al., 2010; Diaconu et al., 2008) and also by active actuation of the internal load-bearing structural members (Barbarino et al., 2009; Grohmann et al., 2008; Milojević and Pavlović, 2016; Previtali et al., 2015).

One concept which has shown promise is the fish bone active camber (FishBAC), introduced by Woods and Friswell (2012). This compliance-based morphing aerofoil device consists of a biologically inspired internal skeletal structure covered with a flexible skin. The skeleton consists of a central bending beam spine with a series of stringers branching off to support the skin, which is made from pre-tensioned elastomeric materials (Figure 2). This structure is highly anisotropic by design, with a low chordwise bending stiffness but high spanwise bending stiffness. The resulting bending moments due to actuation loads then deform the structure in a manner, which produce smooth and continuous changes in camber distribution, with the highly anisotropic nature of the structure focusing the deformation on the desired camber shape.

Figure 2.

FishBAC morphing trailing edge concept (white), attached to a ‘rigid’ wing (black) and actuated by a pair of antagonist tendons (yellow).

Preliminary wind tunnel tests of this structural design show promising results in terms of lift control authority ( $Δ C_{L} \approx 0.72$ ), which is similar to traditional trailing edge plain flaps, but with a considerably higher (i.e. 20%–25% improvement) lift-to-drag ratio that is also less sensitive to angle of attack (Figure 3, Woods et al., 2014). The initial design, analysis and prototyping of the FishBAC concept were performed using an analytical model based on Euler–Bernoulli beam theory (Woods and Friswell, 2014), which accounts for structural discontinuities, due to stringers, and the presence of the skin by varying the bending stiffness distribution (i.e. $EI (x)$ ) across the chord length. One advantage of modelling the out-of-plane behaviour with a beam model is that a one-dimensional (1D) net aerodynamic pressure distribution can be directly integrated and the actuation loads can be accounted for by adding local point forces and moments. However, a beam model is not able to evaluate a two-dimensional pressure distributions (i.e. distributed loads with chordwise and spanwise variations), which are required to capture the three-dimensional aerodynamics of finite wings. Finally, a beam model is not capable of capturing any spanwise variation of displacements due to difference in actuation inputs or material anisotropy.

Figure 3.

Lift-to-drag ratio versus angle of attack of the isotropic FishBAC. Experimentally obtained by wind tunnel testing (Woods et al., 2014).

However, it is precisely these variations in spanwise deformation that are of great interest as they can be exploited to also optimise the spanwise load distribution of the wing. The ability to have significant control over the spanwise distribution of lift on a three-dimensional (3D) wing would provide a number of potential induced drag, structural, control and aeroelastic benefits. The presented model addresses these shortcomings by modelling the FishBAC using a two-dimensional structural formulation based on Kirchhoff–Love plate theory. This approach is more suitable for the intended application as variations in out-of-plate displacement due to several factors, such as bend–twist coupling, differential actuation inputs and 3D aerodynamic effects (Anderson, 2010), can be captured and eventually be exploited to optimise lift distributions in both chord and spanwise directions. Also, another advantage of modelling the structure using plate theory is to exploit material anisotropy by including classical laminate theory (CLT) formulations within the plate’s differential equation. Derivations of the equations of motions and strain energy of composite laminated plates are well established (Whitney, 1987), and consequently, they can be implemented to design a composite FishBAC prototype.

The objective of this work is to develop an analytical model of the FishBAC concept that captures thein-plane and out-of-plane displacement due to both aerodynamic pressure (i.e. transverse loading distribution) and actuation loads (i.e. point moments at the trailing edge). This model is able to analyse the discontinuous geometry of the FishBAC by modelling its variable stiffness as multiple individual plates that are joined by penalty springs, capturing its complexity within a single system of linear equations. The scope of this work not only represents a more capable approach to modelling the FishBAC, but also shows an advance — beyond existing techniques in literature — in modelling complex discontinuous plate structures. Currently, there is evidence in the literature of neither the use of penalty springs to join a significant number of individual plate partitions (115 in this case) nor the use of Rayleigh–Ritz method to model static deflections under transverse pressure and line moments of such complex assembly.

The analytical nature of the model not only allows local stiffness properties to be defined for each individual plate ‘partition’, but also to rapidly modify the geometric parameters (e.g. stringer spacing, spine and skin thickness, wing dimensions, among others) and material properties. Finally, unlike finite element method (FEM), the analytical model is ‘mesh independent’, which means that its convergence only depends on the number of the assumed shape functions that are used.

This article is structured as follows: first, an introduction of the FishBAC morphing concept and the initial geometric modelling assumptions is presented, followed by a brief description of the Rayleigh–Ritz method, assumed shape functions, boundary conditions and model implementation. Finally, a convergence study and comparison with FEM results are introduced as a validation to the developed analytical model.

Structural configuration

As described by Woods and Friswell (2012), the FishBAC’s main load-bearing member is a central ‘spine’ that follows the aerofoil camber line. Attached to it, a series of perpendicular stringers of varying heights support the skin, maintain the aerofoil thickness distribution and also increase the spanwise stiffness of the whole structure. The height of the stringers decreases in accordance with the aerofoil geometry. A pre-tensioned elastomeric skin is attached to the stringers and acts as the aerodynamic surface of the FishBAC. Pre-tensioning the skin significantly reduces the out-of-plane deformations under aerodynamic loading and also prevents skin buckling, which would otherwise occur due to compressive strains induced by camber morphing.

The structure is actuated by a pair of antagonistic tendons each driven by independent actuators that are attached to the top and bottom of the trailing edge, respectively. The antagonist pair transfers the torque generated by the actuators, located at the root of the FishBAC, to the trailing edge of the aerofoil section (Figure 4). In the structural analysis, the actuation loads are modelled as applied moments at the tendon-solid trailing edge joints.

Figure 4.

Schematic of the structural configuration of the fish bone active camber morphing trailing edge concept (Woods and Friswell, 2014).

Finally, the FishBAC trailing edge device integrates readily with the existing structural solution for the primary load-bearing member that carries the majority of the aerodynamic loads. While the cross-sectional geometry of the primary load-bearing spars may vary depending on the particular application (e.g. box section, D-spar, channel section), the FishBAC can be readily adapted to suit. Given its use with traditional high stiffness wing geometries, one of the initial assumptions is that the deflection in the primary spar is negligible, which is to say that a ‘rigid’ front portion of the wing is assumed in this initial portion of the work. As a consequence, the analytical model presented in this article only focuses on the deformation of the compliant morphing trailing edge section, which is assumed to be ‘clamped’ to this rigid front portion of the wing. Therefore, the FishBAC is modelled as a cantilever plate with three free edges. A detailed explanation of the modelling assumptions and boundary conditions is given in the next sections of this article.

Modelling assumptions: geometry and materials

The analytical model developed here assumes that the spine’s initial geometry (Figure 5(a)) is a flat plate with no initial curvature (Figure 5(b)). For symmetric aerofoils, this assumption has no effect, but for cambered aerofoils, it flattens out the small amount of curvature that exists in the camber line over the morphing region. Due to this initial assumption, a number of ‘simplifications’ are applied to the geometry, such as the skin is flat and parallel to the spine and each skin section is located at an equivalent height (from the spine) that is calculated by estimating the equivalent contribution in second moment of area of the curved skin in the original design. Finally, the solid trailing edge section is ‘discretised’ in five sections of constant thickness.

Figure 5.

FishBAC’s geometry (a), simplified geometry (b) and modelling assumption (c). In the FishBAC geometry (a), the spine follows the camber line of the aerofoil section, while in the simplified geometry (b), the spine presents no curvature. Furthermore, the stiffness of each partition is ‘condensed’ at the local mid-plane, and each partition is joined using a series of artificial torsional springs (c).

Furthermore, in order to capture the significant stiffness discontinuities caused by the stringers, the structure is divided into several partitions of uniform thickness distribution and composite stacking sequence (chordwise and spanwise) and the plate’s energy balance is solved in each partition, individually. Each one of these individual ‘plates’ are then joined together by artificial torsional springs at each local boundary. This method is known as the Courant’s penalty method (Ilanko et al., 2015) and is explained in detail in the following sections of this article. Finally, the stiffness of each partition is ‘condensed’ at its mid-plane (Figure 5(c)) using CLT. This means that displacement and rotation compatibility is only enforced at the mid-plane and not at the stringers–skin joints. This assumption is reasonable and expected to be valid for the FishBAC due to the compliance of the skin. As the materials used for the stringers and spine are at least three orders of magnitudes stiffer than the elastomeric skin, there is no risk of structural penetration at the skin–stringer contacts if compatibility is not enforced at these locations.

Figure 5 shows a comparison between the FishBAC’s actual geometry (Figure 5(a)), the assumed geometry in the analytical model (Figure 5(b)) and the ‘condensed’ stiffness assumption at the mid-plane (Figure 5(c)). Note that the geometry of finite element model that is used to validate the analytical model corresponds to the actual geometry (Figure 5(a)), which allows to validate simultaneously the underlying geometry assumptions and the implementation of the methods. Further details about the FEM model and the validation process are discussed later.

For a model of the size of a Unmanned Aerial Vehicle (UAV) wing or helicopter rotor blade, three types of materials are used throughout this analysis: high-strength carbon fibre–reinforced polymer (CFRP), acrylonitrile butadiene styrene (ABS), 3D printed plastic and silicone rubber. The stringers and solid trailing edge sections are modelled using isotropic ABS plastic, while the skin is modelled as isotropic silicone. This material selection for stringers, solid trailing edge and skin was performed based on the materials that were used for manufacturing the isotropic FishBAC wind tunnel prototype (Woods and Friswell, 2014) and also because these also may be used to manufacture the first FishBAC composite prototype. Finally, the composite spine is modelled using high-strength CFRP. Table 1 presents the stiffness values of each one of the material definitions.

Table 1.

Material properties of CFRP (Datoo, 1991), ABS (Cambridge University Engineering Department (2003)– Materials Data Book) and silicone (Woods and Friswell, 2013).

Material	CFRP	ABS	Silicone
$E_{11}$	140 GPa	2.9 GPa	3.18 MPa
$E_{22}$	10 GPa
$G_{12}$	5 GPa	N/A	N/A
$ν_{12}$	0.3	0.35	0.425
$t$	0.125 mm	N/A	N/A

ABS: acrylonitrile butadiene styrene; CFRP: carbon fibre–reinforced polymer.

The shear modulus of both ABS and silicone are obtained using the isotropic Young’s modulus and Poisson’s ratio.

Analytical model

The following section introduces the fundamentals of plate theory that are used to model the behaviour of the FishBAC, as well as the specific procedure that is followed to obtain the displacement fields, including an introduction to the Rayleigh–Ritz method for structural analysis (Whitney, 1987). Furthermore, this section introduces the assumed shape functions and the global boundary conditions that are implemented in the analytical formulation.

Rayleigh–Ritz method

The Rayleigh–Ritz method is a variational method that can be used to approximate solutions to partial differential equations based on energy formulations. Its foundation lies in the principle of conservation of total energy in a closed system. From a mechanics point of view, this implies that the sum of the strain energy of the body and the kinetic and potential energies due to external loads is a stationary value (Ilanko et al., 2015; Whitney, 1987). This approach assumes that no frictional losses exist, which is a reasonable simplifying assumption for many types of structures.

For an initially flat plate, these energy formulations can be written in terms of a total energy expression that is a function of the plate displacement’s $u^{o}, v^{o} and w$ (equation (1))

\begin{array}{l} Π (u^{o}, v^{o}, w) & = U + W + V - T = constant \end{array},

(1)

where $U$ refers to the strain energy of the body, $V$ and $W$ are the potential energies due to transverse and in-plane loads, respectively, and $T$ is the kinetic energy.

Since the scope of this work is to analyse the static displacement of composite FishBAC structures, kinetic energy is neglected for the time being, although it can still be added in later if dynamics are of interest.

Energy definitions

For an elastic body, the total strain energy is defined as the integral of the sum of the products of stresses and strains across the volume of the body as seen in equation (2)

\begin{array}{l} U = \frac{1}{2} ∭ (σ_{x} ϵ_{x} + σ_{y} ϵ_{y} + σ_{z} ϵ_{z} + σ_{x z} ϵ_{x z} \\ + σ_{y z} ϵ_{y z} + σ_{x y} ϵ_{x y}) d x d y d z . \end{array}

(2)

As mentioned earlier, the presented analytical model is based on Kirchhoff–Love plate theory. Therefore, through-thickness and transverse shear strains are neglected (i.e. $ϵ_{z} = ϵ_{xz} = ϵ_{yz} = 0$ ), as stated by Whitney (1987). Furthermore, the in-plane strains of the laminate can be obtained in terms of the plate’s displacement, as seen in equation (3)

\begin{array}{l} ϵ_{x} = \frac{\partial u^{o}}{\partial x} - z \frac{\partial^{2} w}{\partial x^{2}} \\ ϵ_{y} = \frac{\partial v^{o}}{\partial y} - z \frac{\partial^{2} w}{\partial y^{2}} \\ ϵ_{x y} = \frac{\partial u^{o}}{\partial y} + \frac{\partial v^{o}}{\partial x} - 2 z \frac{\partial^{2} w}{\partial x \partial y .} \end{array}

(3)

This assumption leads to a mathematical expression in terms of the plate’s displacements and its derivatives and the stiffness terms, which in this case are expressed in terms of the ABD matrix, when CLT is implemented (equation (4)). The ABD matrix describes the stiffness of the composite laminate; it combines both material and geometric stiffness in a single expression (Hyer, 2014).

Furthermore, the net aerodynamic pressure distribution acting on the FishBAC (which is found separately using an aerodynamic solver, for example, panel methods or computational fluid dynamics (CFD)) can be treated as a transverse pressure distribution, on the plate, with both variations in x and y. The potential energy due to transverse pressure loading (i.e. force per unit area) is defined as the integral of the pressure times the transverse displacement across the surface area (equation (5), Whitney, 1987). The potential energy due to the actuation moments at the trailing edge depends on the first derivative of the transverse displacement along the bending direction and the applied moment intensity M_x (equation (6), Mansfield, 1989). These three expressions can be found as follows:

\begin{array}{l} U_{i j} = \frac{1}{2} \iint A_{11} {(\frac{\partial u^{o}}{\partial x})}^{2} + 2 A_{12} \frac{\partial u^{o}}{\partial x} \frac{\partial v^{o}}{\partial y} + A_{22} {(\frac{\partial v^{o}}{\partial y})}^{2} \\ + 2 (A_{16} \frac{\partial u^{o}}{\partial x} + A_{26} \frac{\partial v^{o}}{\partial y}) (\frac{\partial u^{o}}{\partial y} + \frac{\partial v^{o}}{\partial x}) \\ + A_{66} {(\frac{\partial u^{o}}{\partial y} + \frac{\partial v^{o}}{\partial x})}^{2} - 2 B_{11} \frac{\partial u^{o}}{\partial x} \frac{\partial^{2} w}{\partial x^{2}} - 2 B_{22} \frac{\partial v^{o}}{\partial y} \frac{\partial^{2} w}{\partial y^{2}} \\ - 2 B_{12} (\frac{\partial v^{o}}{\partial y} \frac{\partial^{2} w}{\partial x^{2}} + \frac{\partial u^{o}}{\partial x} \frac{\partial^{2} w}{\partial y^{2}}) - 4 B_{66} \frac{\partial^{2} w}{\partial x \partial y} (\frac{\partial u^{o}}{\partial y} + \frac{\partial v^{o}}{\partial x}) \\ - 2 B_{16} [\frac{\partial^{2} w}{\partial x^{2}} (\frac{\partial u^{o}}{\partial y} + \frac{\partial v^{o}}{\partial x}) + 2 \frac{\partial u^{o}}{\partial x} \frac{\partial^{2} w}{\partial x \partial y}] \\ - 2 B_{26} [\frac{\partial^{2} w}{\partial y^{2}} (\frac{\partial u^{o}}{\partial y} + \frac{\partial v^{o}}{\partial x}) + 2 \frac{\partial v^{o}}{\partial y} \frac{\partial^{2} w}{\partial x \partial y}] + D_{11} {(\frac{\partial^{2} w}{\partial x^{2}})}^{2} \\ + D_{22} {(\frac{\partial^{2} w}{\partial y^{2}})}^{2} + 4 D_{66} {(\frac{\partial^{2} w}{\partial x \partial y})}^{2} + 2 D_{12} \frac{\partial^{2} w}{\partial x^{2}} \frac{\partial^{2} w}{\partial y^{2}} \\ + 4 (D_{16} \frac{\partial^{2} w}{\partial x^{2}} + D_{26} \frac{\partial^{2} w}{\partial y^{2}}) \frac{\partial^{2} w}{\partial x \partial y} d x d y \end{array}

(4)

V_{ij} = - \int \int q (x, y) w (x, y) dx dy

(5)

W_{ij} = - \int M_{x} \frac{\partial w (a_{i}, y)}{\partial x} dy .

(6)

Displacement fields and shape functions

The energy definitions presented in the previous subsection are all in terms of the plate’s displacement and its derivatives, the material and geometric stiffness represented by the ABD matrix terms (equation (4)) and external loads (equations (5) and (6)). In this case, both material properties and external loads are known and treated as inputs, while the displacements are unknown and, therefore, their shapes need to be determined. Within the context of Rayleigh–Ritz method, all three displacements (i.e. u^o, v^o and w) are normally defined in the form of three sets of double summations of terms in the x- and y-direction that satisfy compatibility conditions (equation (7))

u_{ij}^{o} = \sum_{m = 0}^{M} \sum_{n = 0}^{N} Q_{mn}^{ij} X_{m}^{i} (x) Y_{n}^{j} (y)

(7a)

v_{ij}^{o} = \sum_{m = 0}^{M} \sum_{n = 0}^{N} R_{mn}^{ij} X_{m}^{i} (x) Y_{n}^{j} (y)

(7b)

w_{ij} = \sum_{m = 0}^{M} \sum_{n = 0}^{N} S_{mn}^{ij} X_{m}^{i} (x) Y_{n}^{j} (y),

(7c)

where $Q_{mn}$ , $R_{mn}$ and $S_{mn}$ are the displacement amplitudes (unknown constants) and $X_{m} (x)$ and $Y_{n} (y)$ are the assumed shape functions, in x- and y-directions, respectively. Note that since the analytical model performs the energy balance in each one of the plate’s partitions individually, different sets of shape functions need to be defined for each one of the partitions. The subscripts $i$ and $j$ in equation (7) refer to the individual partitions in the chordwise and spanwise directions, respectively.

Previous studies have considered several types of shape function. A common approach in plate mechanics for solving the plate’s differential equation is to assume that the displacement occurs in a periodic form, which makes the use of cosine and sine Fourier series expansions convenient as it allows for closed-form solutions to the differential equation. Several examples of using periodic functions are presented in Green (1944), Fo-Van (1980), Bhaskar and Kaushik (2004) and Khalili et al. (2005), among others.

Another alternative to periodic functions is using orthogonal polynomials. They present better convergence rates when deflections do not occur in a periodic way as they can capture localised features using less expansion terms. In the context of plate mechanics, successful examples of implementing generic orthogonal polynomials are presented by Bhat (1986) and Rango et al. (2013) performed static analysis of fibre-reinforced composite plates using orthogonal polynomials.

Furthermore, a specific set of polynomials, known as the ‘Jacobi family’, are commonly used in structural mechanics. This ‘family’ includes the Gegenbauer polynomials, which is a special case of the Jacobi polynomials and the Chebyshev and Legendre polynomials, which themselves are a special case of the Gegenbauer set (Abramowitz and Segun, 1968; Boyd and Petschek, 2014).

From the structural mechanics point of view, Legendre polynomials have been successfully implemented in several cases, for example, in predicting buckling of highly anisotropic plates (Wu et al., 2012) and discontinuous panels with variable stiffness (Coburn et al., 2014), analysing displacements of variable stiffness beams and plates under transverse pressure loading (O’Donnell and Weaver, 2017) and also in capturing step changes in thickness (Vescovini and Bisagni, 2012). However, their integrals have an exact value of zero when integrated across their normalised domain. This would imply that there is zero net work when an uniform transverse pressure distribution and external moments are applied (equations (5) and (6)). On the other hand, Chebyshev polynomials do not present this property. Therefore, the shape functions that are implemented in this analytical model are Chebyshev polynomials of the first kind, defined as (equation (8))

T (ζ) = \frac{1}{2} [{(ζ - \sqrt{ζ^{2} - 1})^{n} + (ζ + \sqrt{ζ^{2} - 1})}^{n}],

(8)

where n corresponds to the polynomial order.

One important aspect to consider during the analysis is that since Chebyshev polynomials are normalised and defined from $[- 1, 1]$ , the numerical analysis must be performed within this domain. Consequently, two non-dimensional variables ( $ζ$ and $η$ $\in [- 1, 1]$ ) are defined to relate the physical dimensions of each plate to a normalised frame (equation (9)). These two non-dimensional variables are used to replace $x$ and $y$ in equations (4) to (7), which changes the integration bounds from physical boundaries to $[- 1, 1]$ . Figure 6 shows the relationship between physical and normalised frames

ζ_{i} = \frac{2 x_{i}}{a_{i}} and η_{j} = \frac{2 y_{j}}{b_{j}} .

(9)

Figure 6.

Coordinate transformation, from physical $x_{i}, y_{j}$ to normalised $ζ_{i}, η_{i}$ coordinate frames, performed in each one of the partitions.

Global boundary conditions

The FishBAC morphing trailing edge section is modelled as a cantilever plate (Figure 7) clamped to the rigid forward section of the wing. This constraint implies that the root of the FishBAC’s plate must present zero displacement and rotation. Since the Chebyshev polynomials do not naturally meet this condition, the expansion in the chordwise polynomial functions must be modified to enforce the clamped boundary at the first row of plate partitions at the root. Jaunky et al. (1995) introduced the concept of using a circulation function to enforce boundary conditions at any location $ζ_{c}$ (equation (10))

Γ_{x} (ζ) = (ζ - ζ_{c})^{n},

(10)

where the value of n is set depending on the nature of the boundary condition of the $ζ_{c}$ location (Table 2). Note that as long as the circulation term is multiplied to every term of the Chebyshev polynomial expansion, the set remains orthogonal.

Figure 7.

Global coordinate frame of the FishBAC, including global dimensions and locations at where actuation moments are applied. Note that the location where the actuation moments correspond to the tendon-spine points of contact (tendons not shown).

Table 2.

Boundary conditions as implemented by circulation function in equation (10) (Coburn, 2015).

Boundary condition at $ζ_{c}$	n	Displacement	Rotation
Free edge (F)	0	Free	Free
Simply supported (SS)	1	0	Free
Clamped (C)	2	0	0

Stiffness discontinuities and local boundary conditions

The stiffness of the FishBAC is inherently discontinuous due to the presence of the stringers, which implies that the energy balance presented in equation (1) has to be calculated in each section of uniform stiffness as the ABD matrix terms in equation (4) vary significantly between regions with and without stringers. Note that since independent shape functions are used for each individual section (equation (7)), a coordinate transformation from the physical to the normalised frame has to be performed in each partition, individually. Hence, a local coordinate system is defined at the centre of each element and then individually mapped to local $ζ_{i}$ and $η_{j}$ variables (equation (9)) as observed in Figure 6.

In structures with stiffness discontinuities, shear force and bending moments at each ‘joint’ must be continuous when approached from either side of the boundary. However, due to the ‘step’ change in both geometric and material stiffness, curvatures are not continuous. These types of structures are known as ‘C¹-continuous’, where displacement and rotations at local boundaries must be continuous, but higher order derivatives do not. Since Chebyshev polynomials do not inherently meet this type of structural continuity at local boundaries, these have to be enforced by other means (Coburn, 2015).

There are two common approaches for ensuring displacement and rotation continuities: Lagrange multiplier method or Courant’s penalty method (Ilanko et al., 2015). The former one consists of deriving a set of constraint equations that are scaled by unknown coefficients, called Lagrange multiplier, that represent the exact value that the constraints need to be weighted by to enforce continuity. The latter approach consists of using a penalty energy term, analogous to joining each section with a torsional/displacement springs and accounting for the spring energy that is needed to enforce displacement and rotation compatibility (i.e. continuous displacements and rotations). Due to the number of equations and separate Lagrange multipliers that would be needed, the Courant’s penalty method in form of spring penalty energies is is selected in this opportunity. It is worth noting that while these are not the only methods to enforce the ‘C¹-continuity’ (e.g. each individual polynomial set could be modified so they naturally meet this condition), these two approaches are by far the most common in the literature as modifying each polynomial set would be difficult to setup and computationally expensive.

Courant’s penalty method

As mentioned in the previous subsection, each one of the plates sections is assumed to be joined with an artificial penalty spring with a stiffness equal to $k_{k}$ . Given the relevant degrees of freedom (DOFs) between partitions in this analysis, a set of penalty equations for displacements (equation (11)) and out-of plane rotations (equation (12)) are defined as

U_{pu, kl} = \frac{k_{k}}{2} \int_{- b_{j} / 2}^{b_{j} / 2} {(u_{k} (x_{kl}^{(+)}, y_{j}) - u_{l} (x_{kl}^{(-)}, y_{j}))}^{2} dy

(11a)

U_{pv, kl} = \frac{k_{k}}{2} \int_{- b_{j} / 2}^{b_{j} / 2} {(v_{k} (x_{kl}^{(+)}, y_{j}) - v_{l} (x_{kl}^{(-)}, y_{j}))}^{2} dy

(11b)

U_{pw, kl} = \frac{k_{k}}{2} \int_{- b_{j} / 2}^{b_{j} / 2} {(w_{k} (x_{kl}^{(+)}, y_{j}) - w_{l} (x_{kl}^{(-)}, y_{j}))}^{2} dy

(11c)

U_{p w_{x}, kl} = \frac{k_{k}}{2} \int_{- b_{j} / 2}^{b_{j} / 2} {(\frac{\partial w_{k} (x_{kl}^{(+)}, y_{j})}{\partial x} - \frac{\partial w_{l} (x_{kl}^{(-)}, y_{j})}{\partial x})}^{2} dy

(12a)

U_{p w_{y}, kl} = \frac{k_{k}}{2} \int_{- b_{j} / 2}^{b_{j} / 2} {(\frac{\partial w_{k} (x_{kl}^{(+)}, y_{j})}{\partial y} - \frac{\partial w_{l} (x_{kl}^{(-)}, y_{j})}{\partial y})}^{2} dy,

(12b)

where $k$ and $l$ correspond to two adjacent partitions with different stiffness and $x_{kl}$ to the location where they meet. When the spring stiffness $k_{k}$ is ‘large’, the energy is minimised when the difference in displacements and rotations in equations (11) and (12) are minimal. Furthermore, a similar set of penalty expressions are implemented — to account for spanwise discontinuities — where integration with respect to $x$ at locations $y_{kl}$ is performed, instead of with respect to $y$ at $x_{kl}$ .

As mentioned earlier, the values of $k_{k}$ need to be ‘large’ to avoid discontinuities. However, if they are ‘too large’, numerical errors arise due to ill-conditioning of the coefficient matrix. The condition number of a matrix is defined as the product of the Euclidean norm of the matrix and its inverse (equation (13)). The resulting condition number is a measure of how close to singular the matrix is. A low condition number represents that the matrix is ‘well-conditioned’, and hence, inverting it is stable, while a high condition number indicates ‘ill-conditioning’. An ill-conditioned system is much more sensitive to changes in response due to small changes in stimuli (Cline et al., 1979; Groh, 2015), introducing numerical error and decreasing solution stability

κ (A) = ∥ A ∥ ∥ A^{- 1} ∥ .

(13)

Previous studies have selected the stiffness of penalty springs based on convergence studies of their models. Coburn (2015) performed a convergence study based on percentage difference with respect to FEM and estimated that the model was accurate for a penalty stiffness between $k = 1 \times 10^{5} - 1 \times 10^{12} N / m$ . On the other hand, Vescovini and Bisagni (2012) estimated that a penalty stiffness of $k = 1 \times 10^{8} N / m$ was ‘sufficient’ for convergence.

Convergence studies for selecting the magnitude of both chordwise and spanwise penalty springs were performed using the FishBAC’s geometry. These studies showed that for four different spine composite ply stacking sequences, a value of $k = 1 \times 10^{8} N / m$ provides stable results. Figure 8 shows an example of these convergence studies for a $[45 / 45 / 45]_{S}$ spine stacking sequence.

Figure 8.

Convergence study for selecting the magnitude of the penalty springs for Chebyshev terms of $M = N = 5$ primary axis (solid) presents the RMS percentage error with respect to FEM, while the secondary axis (dashed) presents the condition number of the coefficient matrix of the analytical model.

Furthermore, in order to mitigate numerical errors due to high condition number, the coefficient matrix is normalised by dividing each individual row K_T by its root mean square (RMS), as defined in equation (14) (Groh, 2015). In this particular application, this normalisation reduces the condition number of the coefficient matrix by at least four orders of magnitude

\tilde{K_{T}} = \frac{K_{T_{j}}}{\sqrt{\sum K_{T_{j}}^{2}}} .

(14)

Principle of minimum potential energy

As previously stated, the Rayleigh–Ritz method is based on the assumption of conservation of total energy in a closed system. This approach implies that the sum of energies defined in equation (1) has a stationary value. Therefore, differentiating the total energy formulation with respect to any of the unknown constant shape function amplitudes, $Q_{mn}$ , $R_{mn}$ and $S_{mn}$ , leads to a state of minimum energy (Coburn, 2015; Ilanko et al., 2015). Since the total potential energy is quadratic in terms of the unknown coefficients (Whitney, 1987), a system of $3 \times (M \times N) \times (i \times j)$ linearly independent equations can be constructed in each individual FishBAC’s partition (equation (15)), where $i \times j$ refers to the total number of partitions in chord and span directions. Consequently, a single system of linear equations is solved to estimate the value of all the unknown amplitudes of the entire structure, regardless of the number of individual partitions that are used to approximate the FishBAC geometry. This approach is computationally convenient, as the static behaviour of the highly discontinuous geometry can be captured with a single coefficient matrix

\frac{\partial Π}{\partial Q_{m n}^{i j}}, \frac{\partial Π}{\partial R_{m n}^{i j}}, \frac{\partial Π}{\partial S_{m n}^{i j}} = 0 (\begin{matrix} m = 1, 2, \dots, M \\ n = 1, 2, \dots, N \end{matrix} .

(15)

Model implementation

The structural model described in the previous sections was implemented using MATLAB^® R2016a, on an Intel^® Core™ i7-4790 3.60 GHz CPU processor, using a 64-bit OS with 16 GB of physical memory. The geometric dimensions of the model were selected based on a 600 × 900 mm NACA 2510 wing tunnel model that was used for previous experiments. Out of those 600 mm of chord length, the last 140 mm corresponds to the FishBAC, which implies that the morphing device starts at 76.6% of chord length. In other words, the first 460 mm of chord length from the leading edge corresponds to the rigid section of the wing, which is not being modelled in this work.

The derivatives of the shape functions are computed analytically, while the integrals in equations (4) and (5) are computed using ‘integral2’ a two-dimensional adaptive quadrature MATLAB built-in function. Furthermore, the integrals at the partition boundaries are needed for calculating the work due to actuation moments (equation (6)) and penalty energy terms (equations (11) and (12)). However, the analytical derivative of the Chebyshev polynomial recurrence formula (equation (8)) cannot be directly integrated at these locations as it is indeterminate due to a zero denominator. This limitation is solved by calculating the limit of the Chebyshev polynomial as its two boundaries are approached (i.e. at $ζ, η = - 1 and 1$ ), followed by integrating the resulting expressions along the reference length using ‘integral’, a one-dimensional adaptive quadrature MATLAB built-in function.

It is important to mention that all the integrals were performed in the non-dimensional reference frame, defined from $[- 1, 1]$ in both in-plane directions. Since each plate partition is rectangular, and it is assumed that the relationship between the normalised and physical frames (equation (9)) do not vary with increasing displacements, the normalised-to-physical scaling of the integrals can be performed ‘outside’ the integrals as the Jacobian matrix of each partition is constant (and can therefore be factored outside the integral). There are two reasons why this is important: first, the integrals do not need to be computed for every individual partition as they can be later scaled for each individual element, and second, this approach allows all integrals of interest to be computed and stored before performing the structural analysis, which has a significant positive impact in both computation times and RAM memory use.

Furthermore, the FishBAC was divided into a total of 115 partitions: 23 in the chordwise direction and 5 in the spanwise direction. This many chordwise partitions allows all stiffness discontinuities due to stringers to be captured and also to ‘discretise’ the solid trailing edge thick section to avoid steep changes in thickness. In terms of physical dimensions, the spine has a uniform stiffness of 0.75 mm, while each stringer height varies in accordance with the aerofoil thickness distribution, having values that range from 14 to 7 mm. Finally, a uniform skin thickness of 0.5 mm maintains the NACA 2510 aerofoil section. It is important to mention that even though this structural model was developed around the FishBAC concept, it can be easily adapted to model any plate-based structure, regardless of its level of discontinuity, stiffness properties, boundary conditions or dimensions.

Finally, the authors consider that a targeted accuracy between 5% and 10% is acceptable for initial design and elastic tailoring. However, the authors are also aware that this target may need to change depending on the sensitivity of the aerodynamics to structural deflections. Future wind tunnel tests and fluid–structure interaction models will give a better insight of the accuracy needed for modelling this morphing device.

A top view diagram of the partitions, with their respective dimensions, can be found in Figure 9. It is important to mention that the number of partitions and dimensions, stringer spacing and so on can be easily modified by manually changing the input parameters of the MATLAB script file. Also, since the integrals are not computed for each partition individually, increasing the number of partitions does not have a significant impact on the computational cost of the simulation.

Figure 9.

Top view of the analytical model geometry (with local dimensions). A total of 23 chordwise and 5 spanwise partitions, respectively, are used to model the complex geometry of the FishBAC. All dimensions are in millimetre.

Validation: finite element model

A FEM model of the FishBAC was developed, using ABAQUS/CAE^® version 6.14-1, to validate the analytical model. As this FEM model is based on the true FishBAC geometry, comparing the analytical model to it will simultaneously test both the simplified geometry (Figure 5) and mathematical assumptions that are implemented in this analytical model. The FEM model consists of a combination of shell, continuum shell and solid elements that are joined together as a single part. The composite plate, stringers and the skin of the compliant section are modelled using four-node shell (S4R) elements. Furthermore, the spine’s material is defined as a composite laminate (i.e. shell elements with material defined on a ply-by-ply basis), while both the skin and the stringers are modelled as isotropic (Table 1).

The ‘thick’ section of the trailing edge is modelled using a combination of solid 8-node (C3D8R) elements, for the isotropic parts, and continuum shell (SC8R) elements, for the laminated composite parts (Figure 10). This is due to the fact that this region contains a section of the composite spine that is located in between solid isotropic material. A geometrically non-linear analysis was performed for several different spine stacking sequences under a range of uniform pressure distributions, ranging from 20 to 500 Pa. The nodal displacements along the free edges were tracked and extracted to allow for comparison to the analytical model. A fully clamped boundary condition was applied to the root of the FishBAC. Furthermore, the FishBAC’s skin is pre-stretched by 10% to minimise out-of-plane deflections under aerodynamic loading and to avoid buckling in compression. In order to simulate this, the skin is prestressed in the FEM model by applying a prescribed uniform, in-plane predefined stress field equal to the Young’s modulus of the skin times 10% strain in the chordwise direction. Spanwise stress/strain due to Poisson’s ratio effect is not considered, as during the manufacture of actual FishBAC skins, the skin is free to contract in the spanwise direction before it is bonded to the structure. A convergence study was performed to set the element size used for comparison by varying the global element size of the mesh and calculating the percentage difference for each increment. In this study, convergence is achieved when the average percentage difference in tip displacements vary less than 0.5% for two consecutive increments in mesh density.

Figure 10.

Schematic of the type of elements that are used to model the behaviour of the FishBAC’s static displacement. Solid and continuum shell elements are displayed in ‘green’ colour, while shell elements are displayed in ‘pink’ colour.

Loading condition

One of the assumptions in the structural analysis of the FishBAC is that the aerodynamic pressure distributions on the upper and lower skins are translated directly to the spine. This equivalent pressure is found by subtracting the top surface pressure distribution from the bottom one (equation (16)). This approach has one disadvantage: it does not allow the out-of-plane deformation of the skin, under aerodynamic loads, to be modelled. However, for the level of fidelity and purpose of this tool (robust and efficient fluid–structure interaction for design and optimisation), these skin deformations would not be captured with the aerodynamic tools that will be used. A full CFD analysis is needed to estimate the skin deformation between stringers, which would add significant computational expense. A potential alternative is to separately model the skin deflections using CFD and superimpose the results with the current structural model. Finally, as previously mentioned, pre-tensioning the skin mitigates the effects of skin deformation on the aerodynamic loads, as it increases its stiffness, reduces out-of-plane skin deflections and prevents buckling when in compression

Furthermore, the actuating loads are applied as point moments at the location where the tendons meet with the spine. Figure 11 displays the type of loading applied in both analytical and FEM models, while Figure 7 shows the locations where the external actuation moments are applied, including their moment arms with respect to the clamped edge. For the purpose of this initial investigation, $P_{down} (x, y)$ and $P_{up} (x, y)$ are assumed to be uniform, which results in a net positive uniform pressure distribution acting on the spine. Using uniform pressure distributions facilitates the direct comparison between the analytical and FEM models, without having to interpolate to match the location of the FEM nodes and the analytical partitions. Under ‘real’ aerodynamic loads, the rate of deflection would be shallower along the chord than under uniform pressure. The highest pressure would occur at the root of the morphing device and gradually reducing, up to a value of zero at the trailing edge of the device.

Figure 11.

External loads applied in both analytical and finite element models in the global in the global x–y coordinate frame. Dimensions and specific locations of the moment loads are presented in Figures 7 and 9.

The magnitude of the transverse pressure loadings was selected based on preliminary FEM simulations, targeting similar maximum deflections of 4–5 mm for all four spine material configurations, allowing for direct error comparison between all four cases. Note that the highest pressure value applied to each case is similar or higher to the values that the composite wind tunnel model of similar dimensions would experience at Mach 0.15

q_{net} (x, y) = (P_{down} (x, y) - P_{up} (x, y)) .

(16)

Results and discussion

The analytical model presented in this article is compared against a non-linear finite element analysis of the FishBAC under uniform pressure distribution and actuating moments. In addition to the resulting displacement fields, a convergence study in terms of Chebyshev polynomial terms is presented. In these sets of results, two types of errors are reported: the maximum absolute value percentage difference and the RMS percentage difference between analytical and FEM results.

Polynomial term convergence

A convergence study of Chebyshev polynomial terms against the converged FEM results is presented for four different spine material configurations: isotropic ABS, $[0 / 90 / 0]_{S}$ CFRP, $[45 / 45 / 45]_{S}$ CFRP and $[90 / 90 / 90]_{S}$ CFRP (where 0° is aligned with the global x-axis). The spine’s thickness is kept constant at $t = 0.75 mm$ for all four cases, and the dimension, position and material properties of the stringers and skin remain constant. The structure is loaded under uniform transverse pressure and maximum absolute and RMS percentage errors were computed. Furthermore, the study also involved comparing these errors with both computation time and condition number, as observed in Figure 12. Table 3 shows a summary of the convergence study.

Figure 12.

Convergence study: comparison between analytical and geometrically non-linear FEM results for a spine’s stacking sequence of [45/45/45]_S.

Table 3.

Comparison of analytical and geometrically non-linear FEM results in function of Chebyshev polynomials terms (polynomial order).

Laminate	Material	Polynomial terms (M = N)	Max. error (abs. value) (%)	RMS error (%)	DOF
Isotropic	ABS plastic	2	7.978	7.847	3015
		3	7.141	7.846	5520
		4	5.054	4.912	8625
		5	2.332	2.185	12,420
		6	0.791	0.732	16,905
		7	3.634	3.570	22,080
[0/90/0]_S	CFRP	2	7.814	7.711	3015
		3	6.973	6.874	5520
		4	4.150	4.039	8625
		5	0.691	0.580	12,420
		6	3.308	3.129	16,905
		7	6.933	6.745	22,080
[45/45/45]_S	CFRP	2	14.27	9.927	3015
		3	13.28	9.136	5520
		4	10.80	6.903	8625
		5	7.757	4.864	12,420
		6	7.753	4.702	16,905
		7	11.48	6.675	22,080
[90/90/90]_S	CFRP	2	7.825	7.691	3015
		3	7.077	6.974	5520
		4	4.675	4.530	8625
		5	1.730	1.578	12,420
		6	1.642	1.575	16,905
		7	4.715	4.644	22,080

RMS: root mean square; DOF: degree of freedom; ABS: acrylonitrile butadiene styrene; CFRP: carbon fibre–reinforced polymer.

Results show that for all four spine configurations, a RMS error of less than 4.9% can be achieved when five Chebyshev polynomial terms in both chordwise and spanwise direction (i.e. M = N = 5) are used. In general, it is observed that adding a sixth term has a positive effect in further reducing the RMS percentage error, at the expense of doubling the computation time. On the other hand, there is a significant increase in error when adding a seventh polynomial term. This also corresponds to an increase in condition number of about seven orders of magnitude in all four cases, as observed in Figure 12. It can be concluded that the proposed analytical model has a limit of six Chebyshev polynomial terms, as further increasing the number of polynomial terms yields severe ill-conditioning. Finally, it can be observed that results within ≈7% RMS error can be obtained when four polynomial terms are used, which represents an attractive case for future optimisation studies as this case represents a 66% reduction in computational time than when five terms are used. It is also attractive in terms of computational effort with respect to the FEM, as it requires less than 8% of the total number of DOFs as the converged FEM. As a commercial software package, ABAQUS/CAE (FEM) benefits from a highly optimised solution process for its equations of motion, but further optimisation of the solution methods for the process shown here is likely to yield significant improvements in solution time. This approach is in addition to the benefits of simple, parameter-driven geometry/material definition (ideal for optimisation and aeroelasticity) provided by the proposed method.

Regarding the type of error, it is observed that the RMS and the maximum absolute value percentage error are very similar in the isotropic, $[0 / 90 / 0]_{S}$ and $[90 / 90 / 90]_{S}$ . However, these two significantly differ from the most anisotropic case (i.e. $[45 / 45 / 45]_{S}$ ). This implies that the maximum absolute value percentage error in this last case is driven by differences in localised deformation, rather than from the global ones. Since the application of this structural model is to predict deformations due to aerodynamic loads, which are mainly dependent on global deformations, the RMS error is more useful for comparing the model with FEM results and is used subsequently as a measure of comparison to FEM results for the rest of this article.

Finally, the difference in percentage error (Table 3) among the four material configuration is primarily due to the degree of material anisotropy and the presence of transverse shear in the $[45 / 45 / 45]_{S}$ material configuration. A comparison between linear and non-linear FEM was performed, showing no significant differences between the two of them.

Uniform pressure loading

Uniform transverse pressure was applied to the same four laminates presented in the previous subsection. The RMS percentage difference along the spanwise $x = a / 2$ edge is used to measure the difference between the analytical and FEM results.

Table 4 shows that the RMS error remains stable as the load increases. Furthermore, it can be observed that the most anisotropic layup (i.e. $[45 / 45 / 45]_{S}$ ) consistently presents a higher error than the other three cases. The corresponding displacement field is displayed in Figure 13, while Figures 14 and 15 show a comparison between analytical and FEM results along the spanwise $x = a / 2$ and chordwise $y = - b / 2$ edge, respectively.

Table 4.

Comparison, in terms of RMS percentage difference between analytical ( $M = N = 5$ terms) and FEM, of the four different configurations under uniform transverse pressure loading.

Laminate	Material	Pressure (Pa)	RMS error (%)
Isotropic	ABS plastic	20	1.815
		40	1.929
		60	2.047
		80	2.165
[0/90/0]_S	CFRP	100	0.518
		200	0.536
		300	0.553
		400	0.567
		500	0.580
[45/45/45]_S	CFRP	50	4.798
		100	4.823
		150	4.846
		200	4.864
[90/90/90]_S	CFRP	20	1.224
		40	1.308
		60	1.398
		80	1.488
		100	1.577

RMS: root mean square; ABS: acrylonitrile butadiene styrene; CFRP: carbon fibre–reinforced polymer.

Figure 13.

Comparison between analytical (solid) and FEM ( $\circ$ ) displacement field for a $[45 / 45 / 45]_{S}$ spine’s stacking sequence under uniform transverse pressure.

Figure 14.

Analytical (solid) versus FEM ( $\circ$ ) displacement along spanwise edge $x = a / 2$ for a $[45 / 45 / 45]_{S}$ spine’s stacking sequence under uniform transverse pressure.

Figure 15.

Analytical (solid) versus FEM ( $\circ$ ) displacement along chordwise edge $y = - b / 2$ for a $[45 / 45 / 45]_{S}$ spine’s stacking sequence under uniform transverse pressure.

Even though localised discrepancies can be observed along the spanwise edge, the global deformation of the FishBAC is properly captured (Figure 13). As mentioned earlier, properly capturing the global deformation is a priority for this application, as this would drive any change in aerodynamic loads.

On the other hand, the $[0 / 90 / 0]_{S}$ case consistently presents the lowest percentage error out of all four cases. Figure 16 displays the two-dimensional displacement field, while Figures 17 and 18 show a comparison between analytical and FEM results along the chordwise $y = - b / 2$ and spanwise $x = a / 2$ edge, respectively. In this case, both local and global features are well captured by the analytical model. Unlike in the $[45 / 45 / 45]_{S}$ case (Figure 14), it can be observed that there is no significant spanwise variation in out-of-plane displacement in the $[0 / 90 / 0]_{S}$ (Figure 17). This symmetry in spanwise deflection is expected as cross-ply laminates present no bend-twist coupling – $D_{16}$ and $D_{26}$ terms in the ABD matrix are zero.

Figure 16.

Comparison between analytical (solid) and FEM (∘) displacement field, for a $[0 / 90 / 0]_{S}$ spine’s stacking sequence, under uniform transverse pressure.

Figure 17.

Analytical (solid) versus FEM (∘) displacement along the spanwise edge $x = a / 2$ for a $[0 / 90 / 0]_{S}$ spine’s stacking sequence under uniform transverse pressure.

Figure 18.

Analytical (solid) versus FEM (∘) displacement along the chordwise edge $y = - b / 2$ for a $[0 / 90 / 0]_{S}$ spine’s stacking sequence under uniform transverse pressure.

Actuation moments

A comparison of analytical versus FEM displacements, under input moments, was performed for the two most compliant laminates in the chordwise direction (i.e. isotropic and $[90 / 90 / 90]_{S}$ ). The magnitude of the maximum applied moments was selected based on preliminary non-linear FEM results, such that the deflections in both stacking sequence cases had similar values. Three different scenarios were analysed: the first one consisted of applying two positive moments of the same magnitude, the second one, applying moments of equal magnitude but opposite directions, and the third one, applying a negative moment to just one of the actuators. The RMS percentage difference between analytical and FEM results is computed at the spanwise $x = a / 2$ edge, and results are summarised in Tables 5 and 6 for the isotropic and $[90 / 90 / 90]_{s}$ spine laminates, respectively, where a total of six Chebyshev polynomial terms were used in both chordwise and spanwise directions.

Table 5.

Comparison, in terms of RMS percentage difference, between analytical ( $M = N = 6$ terms) and FEM for the isotropic spine configuration under all three moment loading scenarios.

Case	$M_{x_{1}}$ (N m)	$M_{x_{2}}$ (N m)	RMS error (%)
Case 1 (symmetricinput)	0.14	0.14	3.927
	0.28	0.28	5.098
	0.42	0.42	6.204
	0.56	0.56	7.236
	0.7	0.7	8.187
Case 2 (differential input)	0.1	−0.1	43.74
	0.2	−0.2	43.74
	0.3	−0.3	43.72
	0.4	−0.4	43.70
	0.5	−0.5	43.38
Case 3 (singleinput)	0	−0.14	5.550
	0	−0.28	6.371
	0	−0.42	7.342
	0	−0.56	8.406
	0	−0.70	9.526

RMS: root mean square.

Table 6.

Comparison in terms of RMS percentage difference between analytical ( $M = N = 6$ terms) and FEM of the CFRP $[90 / 90 / 90]_{s}$ spine configuration under all three moment loading scenarios.

Case	$M_{x_{1}}$ (N m)	$M_{x_{2}}$ (N m)	RMS error (%)
Case 1 (symmetricinput)	0.14	0.14	3.252
	0.28	0.28	4.110
	0.42	0.42	4.946
	0.56	0.56	5.755
	0.7	0.7	6.534
Case 2 (differential input)	0.1	−0.1	44.09
	0.2	−0.2	44.08
	0.3	−0.3	44.07
	0.4	−0.4	44.05
	0.5	−0.5	44.02
Case 3 (single input)	0	−0.14	6.028
	0	−0.28	6.356
	0	−0.42	6.812
	0	−0.56	7.370
	0	−0.70	8.008

RMS: root mean square.

When comparing displacement fields, it is observed that Case 1 and Case 3 show similar behaviours to FEM, in terms of chordwise and spanwise edge displacement and slopes, with an overall RMS error of less than 8.6%. Figure 19 shows the displacement field of the $[90 / 90 / 90]_{s}$ spine under uniform actuation, while Figure 20 shows the same laminate under single actuation.

Figure 19.

Analytical (solid) versus FEM (∘) for a [90/90/90]_S spine’s laminate under uniform positive actuation inputs.

Figure 20.

Analytical (solid) versus FEM (∘) for a [90/90/90]_S spine’s laminate under a single negative $M_{x_{2}}$ input. Legend corresponds to the moment magnitude of the right actuator.

On the other hand, Case 2 (i.e. differential moment inputs) shows significant discrepancies in terms of RMS error along the spanwise free edge and the deformed shape. Tables 5 and 6 show RMS errors as high as 44%, for the isotropic and $[90 / 90 / 90]_{s}$ cases, respectively. In terms of displacement magnitudes, it can be observed that the analytical model is stiffer than the FEM model (Figure 21) by a factor of approximately 2. After further investigation, it was concluded that these discrepancies are caused by the presence of significant transverse shear deformations along the y–z plane. Since the stringers divide the plate in several partitions, the effective chordwise in-plane global dimension of the FishBAC is reduced; hence, local thickness-to-width ratios of approximately $t / a \approx 1 / 13$ exist. As a consequence, through-thickness normal lines to mid-plane can no longer be assumed to remain straight and normal when significant transverse shear exists.

Figure 21.

Analytical (solid) versus FEM (∘) deflection along chordwise edge y =–b/2, for a [90/90/90]_S spine’s laminate, under differential moment inputs.

Unlike the other load cases that have been analysed, the differential moment input case does present transverse shear as these moment inputs generate a net torque along the spine’s x-axis. As a consequence, shear flow is induced, resulting in transverse shear stresses and strains (i.e. $τ_{yz}$ and $γ_{yz}$ , respectively) that cannot be captured using Kirchhoff–Love plate theory. Further FEM estimations were performed to analyse this case and they confirm this limitation. A comparison between the original FishBAC model and a case in which transverse shear in the y–z direction is artificially ‘suppressed’ (by increasing the CFRP’s $G_{yz}$ shear modulus by four orders of magnitude) is performed. Figure 22 shows the variation along the span of the transverse shear stress $τ_{yz}$ between stringers four and five. It is observed that for the case of realistic $G_{xz}$ , transverse shear stress does exist and varies along the span. Finally, to verify that transverse shear is the source of discrepancy in deformed shapes, the displacement along the chordwise edge $y = - b / 2$ of the analytical model is compared with the FEM ‘suppressed’ transverse shear case. It is observed that similar deformations are achieved, with error of less than 3% (Figure 23). A potential solution to this limitation is to expand the analytical model to account for transverse shear deformation using first-order shear deformation theory (FSDT), also known as Mindlin–Reissner plate theory, in the context of plate mechanics.

Figure 22.

Transverse shear stress $τ_{yz}$ between stringers four and five and along the span. The solid line corresponds to the FishBAC with ‘real’ material properties, while the line with cross markers (+) is obtained when transverse shear deformation in the y–z plane is artificially suppressed.

Figure 23.

Analytical (solid) versus FEM (o) displacement along chordwise edge $y = - b / 2$ , when transverse shear deformations in FEM are ‘artificially’ suppressed.

Computational comparison: model size and convergence study

A final comparison, in terms of DOFs between analytical and FEM model, is presented in this section (Table 7). The FEM’s DOFs correspond to a converged mesh and include all boundary conditions and constraints applied during the analysis. It is observed that among all cases, the analytical model reduces the DOF that need to be solved by at least 84%, compared to a FEM converged mesh. Note that for the $[45 / 45 / 45]_{S}$ case, the FEM model presents about seven times more DOFs than the other three cases. This means that a finer mesh was needed for FEM convergence, which is consistent with the fact that this is the most anisotropic case and experiences more localised displacement features and variations along the span than the isotropic or orthotropic cases.

Table 7.

Analytical versus FEM degrees of freedom.

Laminate	Analytical (DOF)		FEM (DOF)
	(Chebyshev terms)
	M, N = 5	M, N = 6
Isotropic	12,420	16,905	109,346
[0/90/0]_S	12,420	16,905	109,346
[45/45/45]_S	12,420	16,905	821,906
[90/90/90]_S	12,420	16,905	109,346

DOF: degree of freedom; FEM: finite element method.

Finite element information corresponds to converged meshes and includes all applied boundary conditions and constraints.

A reduction in DOFs is not only a measure of the efficiency of the analytical model, as it can obtain converged solutions (except when transverse shear due to differential moment inputs exist) with a smaller linear system, but also significantly decreases the amount of RAM memory required to obtain static deformations.

Conclusion

An efficient and mesh-independent two-dimensional analytical model for predicting the static behaviour of the FishBAC trailing edge device has been developed. The model is capable of analysing fully anisotropic FishBAC geometry and material configurations and predicting the in-plane and out-of-plane displacement fields under external two-dimensional transverse pressure loading and applied actuation moments. It achieves the static modelling by condensing all of the geometric and material features of the FishBAC into a single system of linear equations, obtained using Rayleigh–Ritz method.

Its novelty lies in capturing the complex structure of the FishBAC by discretizing its geometry into a series of individual plates, each one with an equivalent stiffness at the mid-plane, that are joined together using a series of artificial penalty springs. Furthermore, it allows for rapid modification of the structural and material parameters of the FishBAC (e.g. spine stacking sequence, dimensions, stringer spacing), and due to its mesh-independent and analytical nature, it can compute converged displacement fields by solving a fixed number of linear equations that only vary with the polynomial order of the assumed shape functions.

Results show that under uniform pressure loading, the analytical model converges in five Chebyshev polynomial terms, with a percentage error under 4.8% with respect to FEM, using 84% less DOFs. Furthermore, errors in predicting large deflections due to actuation loads range from 3% to 8%, except when a differential moment input is applied. This load case results in a net torque on the FishBAC’s spine, causing transverse shear deformations along the y–z plane. Significant discrepancies exist in this specific load case, and they are due to Kirchhoff–Love plate theory’s inability to capture any through-thickness strains.

Although the analytical model has been developed around the FishBAC morphing trailing edge concept, this approach can be used to model any plate structure, regardless of its level of discontinuity or material properties. Furthermore, this model opens the design space for future design iterations of the FishBAC, not only allowing for the use of composite laminates, but also of core materials (e.g. sandwich configuration). Also, since it is built around Rayleigh–Ritz method, dynamic analysis can be introduced by accounting for kinetic energy, which will be required in the near future to model deflections under unsteady aerodynamic load cases. Also, Rayleigh–Ritz Method can be used to perform both static and dynamic aeroelastic studies.

In summary, the analytical model represents a powerful, robust and fast tool for future design and optimisation of the FishBAC as the structural and material parameters can be easily modified. Future work will address the limitations to accurately predict deflections when transverse shear exist; methods, such as FSDT, will be considered. Finally, this model will act as the foundation of a fluid–structure interaction model that will be able to couple deflections with aerodynamic loads.

Footnotes

Acknowledgements

The first author would like to thank Dr Matt O’Donnell for his support on the numerical implementation of Rayleigh–Ritz method in MATLAB and Mr Sergio Minera for his assistance on transformations between coordinate systems. The authors would also like to thank the reviewers for their constructive suggestions to improve the quality of this article.

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was supported by the Engineering and Physical Sciences Research Council through the EPSRC Centre for Doctoral Training in Advanced Composites for Innovation and Science (grant number EP/L016028/1). The second author would like to acknowledge the Royal Society for the Royal Society Wolfson Merit award and also the Science Foundation Ireland for the award of a Research Professor grant. The third author would like to acknowledge the Royal Academy of Engineering for the Research Professorship award.

Data Access Statement

All underlying data are provided in full within this paper.

ORCID iD

Andres E Rivero

References

Abramowitz

Segun

(eds) (1968) Handbook of Mathematical Functions: With Formulas, Graphs, and Mathematical Tables. 5th ed. New York: Dover Publications.

Anderson

(2010) Fundamentals of Aerodynamics. New York: McGraw-Hill.

Barbarino

Bilgen

Ajaj

et al . (2011) A review of morphing aircraft. Journal of Intelligent Material Systems and Structures 22(9): 823–877.

Barbarino

Pecora

Lecce

et al . (2009) A novel SMA-based concept for airfoil structural morphing. Journal of Materials Engineering and Performance 18(5–6): 696–705.

Bhaskar

Kaushik

(2004) Simple and exact series solutions for flexure of orthotropic rectangular plates with any combination of clamped and simply supported edges. Composite Structures 63(1): 63–68.

Bhat

(1986) Plate deflections using orthogonal polynomials. Journal of Engineering Mechanics 111: 1301–1309.

Bilgen

Friswell

Kochersberger

et al . (2011) Surface actuated variable-camber and variable-twist morphing wings using piezocomposites. In: Proceedings of the 52nd AIAA/ASME/ASCE/AHS/ASC structures, structural dynamics, and materials conference, Denver, CO, 4–7 April.

Boyd

Petschek

(2014) The relationships between Chebyshev, Legendre and Jacobi polynomials: the generic superiority of Chebyshev polynomials and three important exceptions. Journal of Scientific Computing 59(1): 1–27.

Cambridge University Engineering Department (2003) Materials data book. Technical Report, Cambridge University, Cambridge. Available at: http://www-mdp.eng.cam.ac.uk/web/library/enginfo/cueddatabooks/materials.pdf

10.

Cline

Moler

Stewart

et al . (1979) An estimate for the condition number of a matrix. SIAM Journal on Numerical Analysis 16(2): 368–375.

11.

Coburn

(2015) Buckling of stiffened variable stiffness panels. PhD Thesis, University of Bristol, Bristol.

12.

Coburn

Weaver

(2014) Buckling analysis of stiffened variable angle tow panels. Composite Structures 111(1): 259–270.

13.

Datoo

(1991) Mechanics of Fibrous Composites. London: Elsevier.

14.

Daynes

Nall

Weaver

et al . (2010) Bistable composite flap for an airfoil. Journal of Aircraft 47(1): 334–338.

15.

De Gaspari

Ricci

(2010) A two levels approach for the optimal design of morphing airfoils. In: Proceedings of the 13th AIAA/ISSMO multidisciplinary analysis optimization conference, Fort Worth, TX, 13–15 September.

16.

De Gaspari

Ricci

Riccobene

(2014) Design, manufacturing and wind tunnel validation of an active camber morphing wing based on compliant structures. In: Proceedings of the 25th international conference on adaptive structures and technologies, The Hague, 6–8 October.

17.

Diaconu

Weaver

Mattioni

(2008) Concepts for morphing airfoil sections using bi-stable laminated composite structures. Thin-Walled Structures 46(6): 689–701.

18.

Fo-Van

(1980) Bending of uniformly cantilever rectangular plates. Applied Mathematics and Mechanics 1(3): 371–383.

19.

Green

(1944) Double fourier series and boundary value problems. Mathematical Proceedings of the Cambridge Philosophical Society 40(3): 222–228.

20.

Groh

(2015) Non-classical effects in straight-fibre and tow-steered composite beams and plates. Available at: http://aerospaceengineeringblog.com/wp-content/uploads/2013/03/PhD_Thesis_DigitalVersion-1.pdf

21.

Grohmann

Maucher

Prunhuber

et al . (2008) Multidisciplinary design and optimization of active trailing edge for smart helicopter rotor blade. Mechanics of Advanced Materials and Structures 15(3–4): 307–324.

22.

Hyer

(2014) Stress Analysis of Fiber-Reinforced Composite Materials. New Delhi, India: McGraw-Hill.

23.

Ilanko

Monterrubio

Mochida

(2015) The Rayleigh-Ritz Method for Structural Analysis (Iste Series). London and New York: John Wiley & Sons.

24.

Jaunky

Knight

Ambur

(1995) Buckling of arbitrary quadrilateral anisotropic plates. AIAA Journal 33(5): 938–944.

25.

Khalili

Malekzadeh

Mittal

(2005) A new approach to static and dynamic analysis of composite plates with different boundary conditions. Composite Structures 69(2): 149–155.

26.

Kota

Hetrick

Osborn

et al . (2003) Design and application of compliant mechanisms for morphing aircraft structures. Available at: https://www.spiedigitallibrary.org/conference-proceedings-of-spie/5054/0000/Design-and-application-of-compliant-mechanisms-for-morphing-aircraft-structures/10.1117/12.483869.short

27.

Kudva

(2004) Overview of the DARPA smart wing project. Journal of Intelligent Material Systems and Structures 15(4): 261–267.

28.

Larson

(1986) Flight control system development and flight test experience with the F-111 mission adaptive wing aircraft. Technical Report, NASA Ames Research Center, Edwards, CA.

29.

Mansfield

(1989) The Bending and Stretching of Plates. 2nd ed. Cambridge: Cambridge University Press.

30.

Milojević

Pavlović

(2016) Development of a new adaptive shape morphing compliant structure with embedded actuators. Journal of Intelligent Material Systems and Structures 27(10): 1306–1328.

31.

Molinari

Quack

Dmitriev

et al . (2011) Aero-structural optimization of morphing airfoils for adaptive wings. Journal of Intelligent Material Systems and Structures 22(10): 1075–1089.

32.

Monner

Kintscher

Lorkowski

et al . (2009) Design of a smart droop nose as leading edge high lift system for transportation aircrafts. In: Proceedings of the 50th AIAA/ASME/ASCE/AHS/ASC structures, structural dynamics, and materials conference, Palm Springs, CA, 4–7 May.

33.

O’Donnell

Weaver

(2017) RAPID analysis of variable stiffness beams and plates: Legendre polynomial triple-product formulation. International Journal for Numerical Methods in Engineering 112: 86–100.

34.

Previtali

Molinari

Arrieta

et al . (2015) Design and experimental characterisation of a morphing wing with enhanced corrugated skin. Journal of Intelligent Material Systems and Structures 27(2): 1–15.

35.

Rango

Bellomo

Nallim

(2013) A general Ritz algorithm for static analysis of arbitrary laminated composite plates using first order shear deformation theory. Journal of Engineering Research 10(2): 1–12.

36.

Thill

Etches

Bond

et al . (2008) Morphing skins. Aeronautical Journal 112(1129): 117–139.

37.

Valasek

(ed.) (2012) Morphing Aerospace Vehicles and Structures (Aerospace Series). London: John Wiley & Sons.

38.

Vasista

Riemenschneider

Van

Kamp

et al . (2017) Evaluation of a compliant droop-nose morphing wing tip via experimental tests. Journal of Aircraft 54(2): 519–534.

39.

Vescovini

Bisagni

(2012) Buckling analysis and optimization of stiffened composite flat and curved panels. AIAA Journal 50(4): 904–915.

40.

Whitney

(1987) Structural Analysis of Laminated Anisotropic Plates. Lancaster, PA: Technomic Publishing.

41.

Woods

BKS

Friswell

(2012) Preliminary investigation of a fishbone active camber concept. In: Proceedings of the ASME 2012 conference on smart materials, adaptive structures and intelligent systems. Available at: http://michael.friswell.com/PDF_Files/C327.pdf

42.

Woods

BKS

Friswell

(2013) Fluid-structure interaction analysis of the fish bone active camber mechanism. In: Proceedings of the 54th AIAA/ASME/ASCE/AHS/ASC structures, structural dynamics, and materials conference. Boston, MA, 8–11 April, pp. 1–15. Reston, VA: American Institute of Aeronautics and Astronautics.

43.

Woods

BKS

Friswell

(2014) Structural characterization of the fish bone active camber morphing airfoil. In: Proceedings of the 22nd AIAA/ASME/AHS adaptive structures conference, National Harbor, MD, 13–17 January.

44.

Woods

BKS

Bilgen

Friswell

(2014) Wind tunnel testing of the fish bone active camber morphing concept. Journal of Intelligent Material Systems and Structures 25(7): 772–785.

45.

Raju

Weaver

(2012) Comparison of variational, differential quadrature, and approximate closed-form solution methods for buckling of highly flexurally anisotropic laminates. Journal of Engineering Mechanics 139(8): 1073–1083.

Parametric structural modelling of fish bone active camber morphing aerofoils

Abstract

Keywords

Introduction

Structural configuration

Modelling assumptions: geometry and materials

Analytical model

Rayleigh–Ritz method

Energy definitions

Displacement fields and shape functions

Global boundary conditions

Stiffness discontinuities and local boundary conditions

Courant’s penalty method

Principle of minimum potential energy

Model implementation

Validation: finite element model

Loading condition

Results and discussion

Polynomial term convergence

Uniform pressure loading

Actuation moments

Computational comparison: model size and convergence study

Conclusion

Footnotes

Acknowledgements

Declaration of conflicting interests

Funding

Data Access Statement

ORCID iD

References