Static Aeroelastic Response of Wing-Structures Accounting for In-Plane Cross-Section Deformation

International Journal of Aeronautical and Space Sciences.
2013.
Dec,
14(4):
310-323

This is an Open Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License (http://creativecommons.org/licenses/bync/3.0/) which permits unrestricted non-commercial use, distribution, and reproduction in any medium, provided the original work is properly cited.

- Received : October 21, 2013
- Accepted : December 12, 2013
- Published : December 30, 2013

Download

PDF

e-PUB

PubReader

PPT

Export by style

Article

Metrics

Cited by

TagCloud

In this paper, the aeroelastic static response of flexible wings with arbitrary cross-section geometry via a coupled CUFXFLR5 approach is presented. Refined structural one-dimensional (1D) models, with a variable order of expansion for the displacement field, are developed on the basis of the Carrera Unified Formulation (CUF), taking into account cross-sectional deformability. A three-dimensional (3D) Panel Method is employed for the aerodynamic analysis, providing more accuracy with respect to the Vortex Lattice Method (VLM). A straight wing with an airfoil cross-section is modeled as a clamped beam, by means of the finite element method (FEM). Numerical results present the variation of wing aerodynamic parameters, and the equilibrium aeroelastic response is evaluated in terms of displacements and in-plane cross-section deformation. Aeroelastic coupled analyses are based on an iterative procedure, as well as a linear coupling approach for different free stream velocities. A convergent trend of displacements and aerodynamic coefficients is achieved as the structural model accuracy increases. Comparisons with 3D finite element solutions prove that an accurate description of the in-plane crosssection deformation is provided by the proposed 1D CUF model, through a significant reduction in computational cost.
N_{EL}
1D finite elements. A cartesian coordinate system is defined with axes
x
and
z
parallel to the cross-section, whereas
y
represents the longitudinal coordinate. According to the displacementbased framework of CUF
[31]
, the displacement field is assumed to be an expansion of a certain class of functions
F
_{τ}
, which depend on the cross-section coordinates
x
and
z
. Introducing the shape functions
N_{i}
and the nodal displacement vector
q
, the displacement vector
u
with components
u_{x}
,
u_{y}
, and
u_{z}
becomes:
where
q _{τi}
contains the degrees of freedom of the τ-th expansion term corresponding to the
i
-th element node. The compact expression in Eq. 1 is based on Einstein’s notation: repeated subscripts τ and
i
indicate summation. Multivariate Taylor’s polynominals of the
x
and
z
variables are employed here as cross-section funtions
F
_{τ}
, and
N
is defined as the expansion order, which is a free parameter of the formulation. Elements with number of nodes
N_{N}
=4 are formulated in the present work, and named B4, using third-order Lagrange polynomials as shape functions
[19]
. The number of degrees of freedom (DOFs) used through the proposed approach is:
The variational statement employed is the Principle of Virtual Displacements:
where
L_{int}
is the internal strain energy, and
L_{ext}
is the work of external loadings variationally consistent with the present method, and here derived for the case of a generic concentrated load
P
= {
P_{ux} p_{uy} p_{uz}
}
^{T}
, acting on the arbitrary load application point (
x_{p}, y_{p}, z_{p}
), which does not necessarily lie along the 1D finite element mesh, unlike standard 1D FE models. δ stands for the virtual variation. By usuing Eq. 1, δ
L_{ext}
becomes:
where
F
_{τ}
is evaluated in (
x_{p}, z_{p}
) and
N_{i}
is calculated in
y_{p}
. In the case of small displacements with respect to the length
L
, the in-plane (subscript
p
) and out-of-plane (subscript
n
) cross-section stress and strain vectors in Eq. 3 are related to the displacements vector via linear differential matrix operators
D _{p}
,
D _{n}
, and material stiffness matrices
C _{pp}
,
C _{pn}
,
C _{np}
,
C _{nn}
, as follows:
Using Eq. 5, Eq. 3 can be rewritten in terms of virtual nodal displacements:
where, the 3 × 3 and 3 × 1
fundamental nuclei
K ^{ijτs}
and
F _{τi}
are introduced. From Eq. 6, the governing equation of motion can be derived through a finite element assembly procedure:
where
K
is the structural stiffness matrix, and
F
is the vector of equivalent nodal forces. It should be noted that no assumptions on the expansion order have so far been made. Therefore, it is possible to obtain higher-order 1D models without changing the formal expression of the nuclei components, as well as classical beam models, such as Euler-Bernoulli’s and Timoshenko’s. Higher-order models provide an accurate description of the shear mechanics, the in-plane and out-of-plane cross-section deformation, Poisson’s effect along the spatial directions, and the torsional mechanics, in more detail than classical models do. Thanks to the CUF, the present hierarchical approach is invariant with respect th the order of the displacemant field expansion. More details are not reported here, but can be found in the work of
[31]
.
V
(
x, y, z
) is equal to zero:
In this case, the velocity vector
V
(
x, y, z
) can be considered as the gradient of a potential function
ϕ
:
Hence, the analysis of a wing or an airfoil under these conjectures can be performed by potential methods. The potential function describing the flow field around an object can be defined as a combination of singularities, such as doublets, vortices, sources, or uniform flux over the external body surface. According to the detailed exposition in
[36]
, the equation to be used to compute the solution of the aerodynamic problem is Laplace’s equation:
Laplace’s equation describes a potential flow field only if the compressibility effects can be neglected, as occurs for the results presented afterwards, where the free stream velocity is rather low. Otherwise, some corrections, e.g. the Prandtl-Glauert transformation, are necessary, as explained in
[37]
. The assumptions here introduced lead to an intergral-differential equation, which expresses the potential funtion in an arbitrary point of the fluid domain as a combination of singularities. For the sake of brevity, this equation is not reported here, but more details can be found in
[36]
. Among the potential methods, the panel methods can be formulated following a low-order, or a high-order approach. The low-order (first-order) panel method employs triangular or quadrilateral panels having constant values of singularities’ strenth, such as the Hess and Smith approach. The higher-order panel methods instead use higher than first-order panels (e.g. paraboloidal panels) and a varying singularity strength over each panels.
V _{i}
on the generic
i
-th aerodynamic panel with normal vector
n _{i}
is equal to zero:
Further details on this method can be found in
[15]
. The 3D Panel Method schematizes the wing surface as a distribution of doublets and sources. The strength of the doublets and sources is calculated to meet the appropriate boundary conditions (BCs), which may be of Dirichlet- or Neumann-type. According to the creator of the program, after a trial and error process, the best results can be obtained by using just the Dirichlet BC type
[38]
. The 3D Panel Method employs a low-order panel method. The LLT approach is not able to evaluate the pressure coefficients on the wing surface, but only the lifting loads along the lifting line. The VLM is able to analyze the pressure coefficients, but only on the reference surface, which is defined as the mean surface between the upper and the lower wing surfaces. The 3D Panel Method is able to calculate the pressure coefficients on both the upper and the lower wing surfaces. Therefore, this method offers the most realistic description of the aerodynamic field.
N
.
Structural displacements are evaluated in specific sections distributed regularly along the wing span. The cross-section distortion
s
is defined as the in-plane displacemant, i.e. a quantity that express the in-plane difference between the deformed shape and the “undeformed” shape of the airfoil cross-section:
where Δ
u_{x}
and Δ
u_{z}
are the cartesian components of the relative displacement vector Δ
u
along the chord direction
x
and the transversal direction
z
, respectively, between the deformed cross-section and the
base section
. Given a structural model, the base section corresponds to the undeformed cross-section, shifted and rotated in such a way that its leading edge and trailing edge points correspond to the leading edge and trailing edge points, respectively, of the deformed cross-section obtained by such a structural model.
The iterative process in
Fig. 1
is stopped once the convergences of the lifting coefficient
C_{L}
, the moment coefficient
C_{M}
, and the cross-section distortion of the wing sections are achieved simultaneously. The description of a similar iterative process can also be found in
[12]
. The convergence controls are thus:
where
toll
is equal to 10
^{-4}
,
are the lifting coefficient, the moment coefficient, and the average cross-section distortion for the generic i-th iteration, respectively. The average distortion
is defined in Eq. 18. A linear approach is adopted as usual in classical aeroelasticity: for each iteration, the aerodynamic loads computed for the deformed wing configuration are applied to the undeformed configuration, without changing the structural stiffness matrix
K
of Eq. 7.
Aeroelastic iterative procedure, with controllers on aerodynamic coefficients and wing deformation.
pressure, and in this work it is modeled as distributed forces. The generic force acting on the j-th aerodynamic panel is evaluated as:
where
V
_{∞}
is the free stream velocity, and
ρ
_{∞}
is the air density.
A_{j}
is the area of the j-th aerodynamic panel, which the pressure coefficient
refers to. According to XFLR5 notation, normal vectors are considered positive when
is negative, and their verse is outer-pointing. Each aerodynamic force is transferred from the aerodynamic model to the structural model, following the approach described in section 2.1 for the generic concentrated load
P
.
For each iteration, the three-dimensional deformed configuration of the wing is built using 11 airfoils along the half-wing span, at a distance of 0.5 m from each other. The first section lies at the wing root. The wing is discretized through a lattice of quadrilateral aerodynamic panels. Let
be the number of panels along the chord line, and let
be the number of panels along the half-wing span. When the VLM is employed, the total number of aerodynamic panels NAP is equal to
For the 3D Panel Method,
N_{AP}
must be calculated as
where the term
is the number of panels along the wing span, on the upper and lower surfaces of the wing. In addition, the term
is the number of panels placed on the tip lateral cross-sections. For the sake of convenience, only the half-wing is analyzed, since the aerodynamic loads are considered to be symmetric with respect to the wing root.
Geometrical configuration of the straight wing.
Effects of the (a) airfoil thickness and (b) camber line on the spanwise local lifting coefficient Cl of the straight wing, along the y axis. Comparison of the VLM and the 3D Panel Method. V _{∞}=50m/s, ρ _{∞}=1.225kg/m^{3}, α=3°.
choice of
remains the same.
For the present assessment analysis, the free stream velocity is assumed to be
V
_{∞}
=50m/s, such that the compressibility effects can be neglected. The air density is assumed to be
ρ
_{∞}
=1.225kg/m
^{3}
. The angle of attack α of the wing is equal to 3 deg. In all the following analyses, the air density
ρ
_{∞}
and the angle of attack α will be invariable parameters. The results focus on the variation of the spanwise local lifting coefficient
C_{l}
along the wing span, defined as:
where,
c(y)
and
L(y)
are the chord and the Lift Force, respectively, generated by the pressure acting on the panels with span-length 2
e(y)
, placed at the y coordinate. More details can be found in
[32]
. As a first result, the trend of
C_{l}
along the y axis (right half-wing) is reported in
Fig. 3
a. This analysis is carried out, considering the variation of the airfoil thickness. As expected, the VLM is not able to take into account the variation of airfoil thickness, since it computes aerodynamic pressures on the wing reference surface, and underestimates
C_{l}
with respect to the 3D Panel Method. In contrast, the 3D Panel Method is able to evaluate the change of the lifting coefficient as the airfoil thickness increases, as can be seen in
Fig. 3
a.
Figure 3
b reports the trend of the spanwise local lifting coefficient
C_{l}
as the camber line changes. It is evident that both aerodynamic methods are able to analyze the influence of the camber line. Comparing
Figs. 3
a and
3
b, it should be noted that the spanwise local lifting coefficient, and thus the aerodynamic pressures, are affected more by the camber line change than the airfoil thickness change. It can be concluded that the 3D Panel Method is able to provide a more realistic evaluation of the pressure distribution on the wing than the VLM. Moreover, the 3D Panel Method affords pressure loads on the actual wing surface, which are fundamental for an accurate study of the actual wing deformation and airfoil distortion, in lieu of loads applied on a fictitious wing reference surface, as for the VLM case. These reasons make the 3D Panel Method the recommended classical aerodynamic tool for the following aeroelastic wing analyses.
Percent error obtained by different 1D CUF models in the computation of the distortion along the airfoil (a) upper and (b) lower surfaces, at the wing tip cross-section (y=5m). Structural assessment: static wing response to a variable pressure distribution. Reference solution: Nastran solid.
root cross-section (at
y
=0), whereas the tip cross-section is free. The cross-section employed is a 2415 NACA airfoil, with constant thickness equal to 2 mm. A spar with a thickness equal to 2 mm is inserted along the spanwise direction, at 25% of the chord. The isotropic material adopted is aluminum: Young’s modulus
E
=69
GPa
, and Poisson’s ratio
v
=0.33.
Due to the small thickness and the well-known aspect ratio restrictions typical of solid elements, this wing is modeled in MSC Nastran by 214,500 solid Hex8 elements and 426,852 nodes, corresponding to 1,280,556 degrees of freedom (DOFs). The same structure is analyzed through CUF models with a variable expansion order up to
N
=14, and discretized through a 1D mesh of 10 B4 finite elements (31 nodes). The number of DOFs depends on
N
, as expressed in Eq. 2; for instance, with 10 B4 elements and
N
=14, the DOFs are 11,160. However, an analysis of the present structure is also carried out through a Nastran shell FE model, but it is not reported herein, for the sake of brevity. Nonetheless, the error obtained between 1D CUF and shell results is comparable with the error obtained between 1D CUF and
Pressure distribution on the wing along the spanwise direction, for the structural assessment. V _{∞}=50m/s, ρ _{∞}=1.225kg/m^{3}, p_{ref} = ½ρ∞V ∞^{2} = 551.25 Pa .
solid results.
A variable pressure distribution step-like along the spanwise direction is applied to the upper and lower wing surfaces, in order to simulate a real pressure distribution, see
Table 1
. The static structural response of the wing is evaluated in terms of the distortion s at the tip cross-section. For the upper and lower surfaces,
Figs. 4
a and
4
b show the percent error
e
obtained by computing the distortion through 1D CUF models and the Nastran solid model, which is taken as reference:
As depicted in
Figs. 4
a and
4
b, the proposed 1D FEs provide a convergent solution, by gradually approaching the Nastran solid results, as the expansion order increases from 8 to 14, according to the conclusions made in previous CUF works
[39]
. For
N
=14, the maximum percent error is about 3% for the upper surface, and about 2.7% for the lower surface. For the wing configuration considered, the choice of
N
=14 seems hence to be accurate enough to detect the cross-section distortion with an acceptable error with respect to the Nastran 3D results, and with a remarkable reduction in terms of DOFs (about a 91% reduction, 11,160 vs. 1,280.556).
V
_{∞}
=30m/s, via the iterative CUF-XFLR5 procedure. This analysis aims at evaluating the influence of the CUF expansion order N on the aeroelastic behavior of the structure, as the accurate description of the cross-section distortion depends on
N
. The same material and straight wing configuration as those considered in the previous
Convergence of lifting and moment coefficients in the iterative aeroelastic analysis, for structural models with different accuracy. Aerodynamic method: 3D Panel. V _{∞}=30m/s.
assessment are employed here, see
Fig. 2
a. In this case, the cross-section is the NACA 2415 airfoil, depicted in
Fig. 2
b. The spar thickness
t
_{3}
is constant and equal to 2 mm; whereas, the skin thickness of upper and lower surfaces varies gradually, from 2 mm (
t
_{1}
in
Fig. 2
b) to 1 mm (
t
_{2}
in
Fig. 2
b), in the zone between 40% and 45% of the chord. This particular choice is coherent with the purpose of studying a highly-deformable nonclassical cross-section.
The 1D structural mesh consists of 10 B4 elements. For the sake of brevity, a convergent study on the number of mesh elements is not reported here. In fact, the choice of 10 B4 elements yields a good evaluation of displacements for all the points of the structure, as detailed in
[32
,
34]
, where a similar structural case in terms of wing configuration and applied aerodynamic loads was studied via the present structural model, and successfully assessed with a commercial FE solid model.
The aeroelastic analysis is now carried out following the iterative coupled procedure CUF-XFLR5 described in
Fig. 1
, and varying
N
. The convergence process on the lifting and moment coefficients is drawn in
Fig. 5
a, by means of a dimensionless parameter
C_{L}/C_{L}^{CONV}
, and in
Fig. 5
b,
Convergent values of lifting coefficient C_{L}^{conv} and moment coefficient C_{M}^{conv} , for different structural models. V _{∞}=30 m/s, C_{L}^{in} =0.4637, C_{M}^{in} =- 0.1629.
respectively.
C_{L}^{conv}
is the final convergent value of the lifting coefficient, which is different for each expansion order employed, as well as the final convergent moment coefficient
C_{M}^{conv}
, as reported in
Table 2
.
Hence, a different choice of
N
influences the structural response of the wing to the aerodynamic loads, and consequently also affects the aerodynamic analysis, due to the aeroelastic coupling. The higher the expansion order employed, the more difference appears between
C_{L}^{conv}
(
C_{M}^{conv}
) and the initial value
C_{L}^{in}
(
C_{M}^{in}
) evaluated for the undeformed wing. For the cases presented in this work, the number of iterations required to achieve the convergence of the lifting coefficient is the same as that required to achieve the convergence of the moment coefficient. It can be seen that the increase of
N
corresponds to the increasing number of iterations
required to achieve the convergence of aerodynamic coefficients. This tendency will be clearly explained afterwards, as a consequence of the introduction of higher-order terms in the model formulation, which enriches the displacement field.
An average cross-section distortion
is now introduced in order to evaluate the aeroelastic deformation of the cross-section
Spanwise distribution of the average distortion s^{-} of the airfoil cross-sections, for different structural models. V _{∞}=30m/s.
“-” : convergence achieved with a tolerance toll = 10^{-4}
shape along the wing span. Given an airfoil cross-section, the average distiortion
is defined as:
where
l
is the curvilinear coordinate along the external airfoil surface, and s is the distortion of the single point of the external airfoil surface defined in Eq. 12. It is noteworthy that s is a positive quantity, and a null value for the average distortion
means no distortion.
Figure 6
plots the trend of the average distortion along the wing span, showing which are the most in-plane deformed airfoil cross-sections in the static aeroelastic equilibrium response. A remarkable variation in the trend of the average distortion appears, depending on the accuracy of the structural model chosen. Models with an expansion order higher than 9 reveal that the section at
y
=4 m appears to be the most distorted section.
For this cross-section,
Table 3
presents the numerical values of average distortion
in the iterative aeroelastic analysis, for different structural theories. As occurred for the convergence of aerodynamic coefficients, the number
Deformation of the airfoil cross-section at y =4m, computed for structural models with different accuracy. V _{∞}=30m/s.
of iterations
required to achieve the convergence of
increases as
N
, and consequently DOFs, increase. In fact, increasing the expansion order
N
, the structural model becomes, in general, more deformable, approaching the real structural behavior. This means that a complete three-dimensional displacement field, as well as local effects, are evaluated with increasing accuracy, especially for structures with highly-deformable cross-sections, see
Figs. 4
a and
4
b. Since the model accuracy increases, the structural deformation is therefore more sensitive to the variations of aerodynamic loads, which are different for each iteration, following the convergent trend in
Figs. 5
a and
5
b, leading to an increasing
The numerical results in
Table 3
highlight that, given an expansion order, a higher number of iterations is necessary to achieve convergence on structural distortion than convergence on aerodynamic coefficients
although the tolerance employed is
Distortion of the airfoil upper surface of the cross-section at y =4m, computed for different structural models. V _{∞}=30m/s.
Convergence of lifting and moment coefficients in the iterative aeroelastic analysis, for different free stream velocities. Structural model: N =14. Aerodynamic method: 3D Panel.
the same.
For
N
>8, the displacement field becomes accurate enough to relevantly take into account a cross-section distortion for the airfoil case considered, as can also be seen in
Fig. 7
. As previously explained, given a structural model, the distortion is computed by comparing the deformed cross-section to the corresponding base section. For the sake of simplicity, only the base section for
N
=1 is plotted in
Fig. 7
.
As expected, low-order models provide a correct evaluation of the bending and torsional structural behavior, but not an exhaustive description of the in-plane deformation. This conclusion is confirmed by
Fig. 8
, where the airfoil distortion
s
computed by variable kinematic models is depicted along the upper surface, at
y
=4m. The maximum distortion value is reached in the part of the cross-section next to the trailing edge, since the stiffening effect due to the spar at 25% of the chord limits the cross-section distortion. Nonetheless, the chordwise position of the maximum distortion points on the airfoil upper and lower surfaces changes, depending on the accuracy of the structural model, see
Table 4
. As a consequence, it is worth pointing out that the increase of N is relevant, not only for an accurate detection of distortion values, but also of the accurate shape-type deformation.
In general, improvements of the structural theory yield more realistic deformations of the wing, until a good convergence is achieved for
N
=14, according to the conclusions made for
Figs. 4
a and
4
b in the structural assessment. In other words, the difference between the results obtained through the generic (
N
-1)-th and N-th expansion orders decreases and becomes minimal for
N
=14. For this reason, it is possible to consider the fourteenth-order model sufficiently accurate to describe the aeroelastic behavior of the structure here considered.
Deformation of the airfoil cross-section at y =4m, computed for different free stream velocities. Structural model: N =14.
Convergent average distortion s^{-conv} [mm] of the cross-section at y=4m, for different structural models. Values and chordwise positions of the maximum distortions s^{US}_{max} [mm] and s^{LS}_{max} [mm], on the airfoil upper and lower surfaces. V _{∞}=30m/s.
Convergent values of lifting coefficient C_{L}^{conv} moment coefficient C_{M}^{conv} and average distortion s^{-conv} [mm] of cross-section at y=4m, for different free stream velocities V _{∞} [m/s]. Structural model: N =14. C_{L}^{in} =0.4637, C_{M}^{in} =-0.1629
Values and chordwise positions of the maximum distortions s^{US}_{max1} , s^{US}_{max2} , s^{US}_{max1} , s^{US}_{max2} , s^{US}_{max1} [mm] on the airfoil upper and lower surfaces of the cross-section at y =4m, for different free stream velocities V _{∞} [m/s]. Structural model: N =14.
the structural model considered is
N
=14. The free stream velocities considered are 25, 30, and 35 m/s. As in the previous analysis, the aerodynamic convergence process is presented through the dimensionless parameter
C_{L}/C_{L}^{conv}
, as illustrated in
Fig. 9
a. The convergence of the moment coefficient is also shown in
Fig. 9
b, through the parameter
C_{M}/C_{M}^{conv}
.
Distortion of the airfoil upper and lower surfaces of the cross-section at y =4m, computed for different free stream velocities. Structural model: N =14.
In this case,
C_{L}^{conv}
and
C_{M}^{conv}
represent the final convergent values of the lifting and moment coefficient for a given
V
_{∞}
, respectively. As occurred for the previous aeroelastic analysis, the trends do not show any numerical problems, such as oscillations. From
Figs. 9
a and
9
b, it is important to note that the number of iterations
required to achieve the aerodynamic convergence increases as
V
_{∞}
increases, and the final convergent values are much different from the initial values, as summarized in
Table 4
. The reason for this behavior is easily explained by the fact that an increasing free stream velocity means increasing aerodynamic loads, and consequently higher structural deformations, and lastly, a more relevant coupling effect on the aeroelastic response of the wing. In fact, an increasing airfoil distortion for the most deformed cross-section at
y
=4m is obtained with
V
_{∞}
according to the numerical results in
Table 5
and airfoil deformed profiles in
Fig. 10
. Also for velocity values different from 30m/s, a higher number of terations is necessary to achieve convergence on structural distortion than convergence on the aerodynamic lifting coefficients
see
Table 5
.
The limitation of distortion close to the airfoil leading edge, due to the spar, is enhanced for
V
_{∞}
=35m/s. The trends of distortion on the airfoil upper and lower surfaces, which are indicated as US and LS, respectively, are depicted in
Fig. 11
at
y
=4, for different velocities. It is important to note that deformations of the upper and lower surfaces also remarkably differ, because of different aerodynamic pressure distributions.
Table 6
shows that not only the maximum distortion values on the airfoil upper
and lower
surface changes as
V
_{∞}
varies, but also their corresponding chordwise positions. This aspecthighlights the importance of higher-order models, in particular for an accurate evaluation of the in-plane cross-section distortion of highly-deformable structures.
As far as the present hierarchical one-dimensional approach is concerned, the results point out that:
These reasons make the future use of the proposed CUF-XFLR5 approach appear promising for a versatile flight optimization tool.

aeroelasticity
;
higher-order 1D finite elements
;
Carrera Unified Formulation
;
in-plane cross-section deformation

1. Introduction

The in-depth understanding of aeroelastic effects on deformable lifting bodies (LBs), due to steady and unsteady aerodynamic loadings, is a typical challenging issue for the current design of aerospace vehicles
[1]
. Furthermore, with the forthcoming employment of composite materials in nextgeneration aircraft configurations, such as High-Altitude Long Endurance aircraft (HALE)
[2]
, and strut-braced wings
[3]
, accurate evaluation of the aeroelastic response becomes even more crucial
[4]
.
Recently, special attention has been directed to the profitable exploitation of the aeroelastic phenomena comprehension , by studying the concept of morphing wings, which are able to adapt and optimize their shape depending on the specific flight conditions and mission profiles
[5
,
6]
. The smart wing is very flexible and could allow a number of advantages, such as drag reduction and aeroelastic vibrations suppression by means of adaptive control
[7]
and different solutions, such as compliant structures
[8]
, bi-stable laminate composites
[9]
, piezoelectric
[10]
and shape memory alloy actuation
[11]
.
In order to develop aeroelastic tools that are able to work in any regime and with any LB geometry, the literature from the last decades has been widely influenced by research devoted to build reliable methods, to couple computational fluid dynamics (CFD) or classical aerodynamic methods with the finite element method (FEM) for structural modeling
[12]
. Valuable examples are in the review articles by Dowell and Hall
[13]
, and Henshaw et al.
[14]
. Reduced approaches, for instance panel methods, are widely used for the classical aerodynamics of wings under some limitations
[15]
, allowing a sizeable reduction in terms of computational cost
[16
,
17]
. The assumption of undeformable airfoil-crosssections
[18]
is typically not proper for recent configurations, since the weight reduction makes wings more flexible and highly-deformable. Hence, two-dimensional (2D) plate/shell and three-dimensional (3D) solid methods are usually employed for the structural modeling, instead of classical one-dimensional (1D) theories, such as the Euler-Bernoulli, Timoshenko, or Vlasov theories
[19]
.
With the advent of smart wings, detailed structural and aeroelastic models are even more essential to fully exploit the non-classical effects in wing design, due to the properties characterizing advanced composite materials, such as anisotropy, heterogeneity and transverse shear flexibility. Beam-like components can be analyzed by means of refined one-dimensional (1D) formulations, which have the main advantage of a lower computational cost required compared with 2D and 3D models. A detailed review of the recent development of refined beam models can be found in
[20]
. El Fatmi
[21]
improved the displacement field over the beam cross-section, by introducing a warping function, to refine the description of normal and shear stress of the beam. Generalized beam theories (GBT) originated with Schardt’s work
[22]
and improved classical theories, by using a piecewise beam description of thin-walled sections
[23]
. An asymptotic type expansion, in conjunction with variational methods, was proposed by Berdichevsky et al.
[24]
, where a commendable review of prior works on beam theory development was given. An alternative approach to formulating refined beam theories, based on asymptotic variational methods (VABS), has led to an extensive contribution in the last few years
[25]
.
A considerable amount of research activity devoted to aeroelastic analysis and optimization was undertaken in the last decades, by using reduced 1D models. A review was carried out by Patil
[26]
, who investigated the variation of aeroelastic critical speeds with the composite ply layup of box beams, via the unsteady Theodorsen’s theory. A thin-walled anisotropic beam model in-corporating nonclassical effects was introduced by Librescu and Song
[27]
to analyze the sub-critical static aeroelastic response, and the divergence instability of swept-forward wing structures. Qin and Librescu
[28]
developed an aeroelastic model to investigate the influence of the directionality property of composite materials, and non-classical effects on the aeroelastic instability of thin-walled aircraft wings. Among the several composite rotor blades applications, the work done by Jeon and Lee
[29]
concerning the steady equilibrium deflections, via a large deflection type beam theory with small strains, is worth mentioning. An example of the use of a refined beam theory for aeroelastic analysis can be found in
[30]
, where the static and dynamic responses of a helicopter rotor blade are evaluated by means of a YF/VABS model.
Higher-order 1D models with generalized displacement variables, based on the Carrera Unified Formulation, have recently been proposed by Carrera and co-authors, for the static and dynamic analysis of isotropic and composite structures
[31]
. The CUF is a hierarchical formulation, which considers the order of the model as a free-parameter of the analysis. In other words, models of any order can be obtained, with no need for ad hoc formulations, by exploiting a systematic procedure. Structural 1D CUF models were used to analyze the structural response of isotropic aircraft wings, under aerodynamic loads computed through the Vortex Lattice Method (VLM), in
[32]
. The aeroelastic CUF-VLM coupling was preliminarily formulated in
[33]
for isotropic flat plates, and then extended to instability divergence detection and the evaluation of composite material lay-up effects on the aeroelastic response of moderate and high-aspect ratio wing configurations, in
[34]
. Flutter analyses of composite lifting surfaces were also presented in
[35]
, by coupling the CUF approach with the Doublet Lattice Method.
The present work couples a refined one-dimensional finite element model based on CUF to an aerodynamic 3D Panel Method, implemented in the software XFLR5. Two potential methods are here compared: the VLM and the 3D Panel Method. The aeroelastic static response of a straight wing is computed through a coupled iterative procedure, and a linear coupling approach. Particular attention is drawn to the in-plane deformation of the wing airfoil cross-sections, as well as the aeroelastic influence of free stream velocity.
2. Numerical models: refined 1D CUF model and panel methods

- 2.1 Variable kinematic 1D CUF FE model

For the sake of completeness, some details about the formulation of CUF finite elements are here retrieved from previous works
[32
,
34]
. A structure with axial length L and cross-section Ω is discretized through a mesh of
Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

- 2.2 A numerical approach for wing aerodynamic analysis

- 2.2.1 Preliminaries

The evaluation of aerodynamic loads can be typically carried out through a CFD code, which solves, for emample, either Navier-Stokes equations or Euler equations numerically. This kind of analysis has a high computational cost, but under some assumptions it is possible to employ simplified approaches. In the wing cases considered in the present work, the flow field is assumed to be steady, and the fluid viscosity is not decisive since the viscous effects can be confined into a small region (boundary layers and wake regions). The fluid can be thus considered as inviscid, and the flow field is irritational, since the curl of the velocity vector
Lager Image

Lager Image

Lager Image

- 2.2.2 XFLR5: an implementation of aerodynamic potential methods

XFLR5 is a software developed by the Andre Deperrois. It performs viscous and inviscid aerodynamic analysis on airfoils and wings, using three potential methods: the Lifting Line Theory (LLT), the VLM, and the 3D Panel Method. The LLTmethod derives from Prandtl’s wing theory and considers the wing as a linear distribution of vortices. The VLM considers a wing as an infinitely thin lifting surface, via a distribution of vortices placed over a wing reference surface. This method requires the non-penetration condition on the reference surface as a boundary condition. Hence, the normal component of the induced velocity
Lager Image

3. Aeroelastic static response ananlysis via the CUF-XFLR5 approach

In this work the aeroelastic static response of the wing is computed through an iterative procedure, based on a coupled CUF-XFLR5 method. Hence, the aerodynamic analysis in XFLR5, as previously mentioned; whereas, variable kinematic 1D CUF models provide the structural wing deformation with a variable expansion order
- 3.1 Iterative procedure

Figure 1
shows in detail the aeroelastic iterative process, which starts with the evaluation of the pressure coefficients for the underformed wing configuration. The further steps to be repeated for each iteration are:
- 1. post-processing calculation of the aerodynamic forces;
- 2. structural analysis of the wing, subject to the aerodynamic forces previously computed;
- 3. new calculation of the aerodynamic pressure coefficients for the new deformed configuration;
- 4. post-processing evaluation of the wing deformation and cross-section distortion.

Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

- 3.2 Aerodynamic loads computation

The aerodynamic load computed by XFLR5 is a distributed
Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

4. Numerical results

- 4.1 Aerodynamic assessment

Firstly, an aerodynamic assessment of the VLM and the 3D Panel Method, which are able to evaluate the pressure coefficients on the wing surface, is performed analyzing the effects of two typical geometrical parameters: the airfoil thickness and the camber line. A straight wing is considered: the wing span is 10 m, and the airfoil chord is 1 m long, as drawn in
Fig. 2
a, where the right half-wing is depicted. This wing configuration is also used in the following structural and aeroelastic analyses. The effect of the camber line on the aerodynamic field is evaluated, using NACA 2415, 3415 and 4415 airfoils. The analysis of the influence of the airfoil thickness is then carried out, using the symmetric NACA 0005, 0010 and 0015 airfoils. The number of aerodynamic panels is chosen as a compromise between the limit number of panels that can be used in XFLR5 (= 5,000)
[38]
, and the number of panels required, in order to achieve convergence in the aerodynamic results. In the following analyses, the
Lager Image

Lager Image

Lager Image

Lager Image

- 4.2 Structural assessment

In order to validate the results given by the proposed higher-order 1D CUF approach, a comparison of the static structural wing response is here performed, with MSC Nastran. Only the right half-wing of the straight configuration introduced in the previous aerodynamic assessment (see
Fig. 2
a) is considered here, due to loads and structural symmetry. A clamped boundary condition is taken into account for the
Lager Image

Pressure distribution on the wing along the spanwise direction, for the structural assessment.V∞=50m/s,ρ∞=1.225kg/m3,pref= ½ρ∞V∞2= 551.25Pa.

Lager Image

Lager Image

- 4.3 Aeroelastic coupling

This section focuses on the results regarding the equilibrium aeroelastic response of a wing exposed to a free stream velocity
Lager Image

Convergent values of lifting coefficientCLconvand moment coefficientCMconv, for different structural models.V∞=30 m/s,CLin=0.4637,CMin=- 0.1629.

Lager Image

Lager Image

Lager Image

Lager Image

Convergence of the average distortions-[mm] in the iterative aeroelastic analysis, for different structural models. Airfoil cross-section aty=4m.V∞=30m/s.

Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

- 4.4 Free stream velocity influence

This analysis aims at establishing the influence of the free stream velocity on the wing distortion. The wing configuration employed for this analysis is the same as that used in the previous study. According to the previous conclusion,
Lager Image

Convergent average distortions-conv[mm] of the cross-section at y=4m, for different structural models. Values and chordwise positions of the maximum distortionssUSmax[mm] andsLSmax[mm], on the airfoil upper and lower surfaces.V∞=30m/s.

Lager Image

Convergent values of lifting coefficientCLconvmoment coefficientCMconvand average distortions-conv[mm] of cross-section at y=4m, for different free stream velocitiesV∞[m/s]. Structural model:N=14.CLin=0.4637,CMin=-0.1629

Lager Image

Values and chordwise positions of the maximum distortionssUSmax1,sUSmax2,sUSmax1,sUSmax2,sUSmax1[mm] on the airfoil upper and lower surfaces of the cross-section aty=4m, for different free stream velocitiesV∞[m/s]. Structural model:N=14.

Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

5. Conclusions

Variable kinematic 1D finite elements were formulated on the basis of the Carrera Unified Formulation (CUF) and coupled to an aerodynamic 3D panel method, implemented in XFLR5. The aeroelastic static response of a straight wing with a highly-deformable airfoil cross-section was computed through a coupled iterative procedure, for increasing structural accuracy and for different free stream velocities. An aerodynamic assessment confirmed that the 3D Panel Method provides a more realistic evaluation of the pressure distribution on the wing, than the Vortex Lattice Method (VLM). As far as the use of 1D higher-order models is concerned, the following main conclusions can be drawn:
- 1. The introduction of higher-order terms in the displacement field is even more important for the aeroelastic analysis, due to the fluid-structure coupling.
- 2. In the case that the wing is rather flexible, the in-plane cross-section deformation has a great impact on the alteration of the aerodynamic loadings.
- 3. The higher the free stream velocity, the more marked the in-plane distortion effect.

- a. The CUF is an ideal tool to easily compare different higher-order theories, since the model accuracy is a free parameter of the analysis.
- b. The in-plane airfoil cross-section deformation is well-described by the proposed 1D structural model, in good agreement with a three-dimensional FE solution, and with a remarkable reduction in terms of DOFs.
- c. A convergent trend of displacements and aerodynamic coefficients is achieved as the structural model accuracy increases. This proves that the proposed 1D higher-order approach does not introduce additional numerical problems in the aeroelastic analysis of wings with arbitrary cross-section geometry.
- d. A higher number of iterations is necessary to achieve convergence on structural distortion than for convergence on aerodynamic coefficients.

Fung Y.C.
2008
An Introduction to the Theory of Aeroelasticity
Dover Publications

Patil M.J.
,
Hodges D.H.
,
Cesnik C.E.S.
2001
“Nonlinear aeroelasticity and flight dynamics of high-altitude long-endurance aircraft”
Journal of Aircraft
38
(1)
88 -
94
** DOI : 10.2514/2.2738**

Sulaeman E.
,
Kapania R.
,
Haftka R.T.
“Parametric studies of flutter speed in a strut-braced wing”
Proceeding of the 43rd AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference
Denver, CO
April 2002
AIAA Paper 2002-1487

Hodges D.H.
,
Pierce G.A.
2002
Introduction to Structural Dynamics and Aeroelasticity
Cambridge University Press

Sofla A.Y.N.
,
Meguid S.A.
,
Tan K.T.
,
Yeo W.K.
2010
“Shape morphing of aircraft wing: status and challenges”
Materials & Design
31
(3)
1284 -
1292
** DOI : 10.1016/j.matdes.2009.09.011**

Gamboa P.
,
Vale J.
,
Lau F.J.P.
,
Suleman A.
2009
“Optimization of a morphing wing based on coupled aerodynamic and structural constraints”
AIAA Journal
47
(9)
2087 -
2104
** DOI : 10.2514/1.39016**

Suleman A.
,
Costa A.P.
2004
“Adaptive control of an aeroelastic flight vehicle using piezoelectric actuators”
Computers & Structures
82
(17-19)
1303 -
1314
** DOI : 10.1016/j.compstruc.2004.03.027**

Kota S.
,
Hetrick J.A.
,
Osborn R.
,
Paul D.
,
Pendleton E.
,
Flick P.
,
Tilmann C.
“Design and application of compliant mechanisms for morphing aircraft structures”
Proceedings of the SPIE Vol. 5054, Smart Structures and Materials 2003: Industrial and Commercial Applications of Smart Structures Technologies
San Diego, CA
12 August 2003

Diaconu C.G.
,
Weaver P.M.
,
Mattioni F.
2008
“Concepts for morphing airfoil sections using bi-stable laminated composite structures”
Thin-Walled Structures
46
(6)
689 -
701
** DOI : 10.1016/j.tws.2007.11.002**

Lim S.M.
,
Lee S.
,
Park H.C.
,
Yoon K.J.
,
Goo N.S.
2005
“Design and demonstration of a biomimetic wing section using a lightweight piezo-composite actuator (LIPCA)”
Smart Materials and Structures
14
(4)
496 -
503
** DOI : 10.1088/0964-1726/14/4/006**

Elzey D.M.
,
Sofla A.Y.N.
,
Wadley H.N.G.
2005
“A shape memory-based multifunctional structural actuator panel”
International Journal of Solids and Structures
42
(7)
1943 -
1955
** DOI : 10.1016/j.ijsolstr.2004.05.034**

Bhardwaj M.K.
1995
“Aeroelastic analysis of modern complex wings using ENSAERO and NASTRAN”
Technical Report NASA 19980235582

Dowell E.H.
,
Hall K.C.
2001
“Modeling of fluidstructure interaction”
Annual Review of Fluid Mechanics
33
445 -
490
** DOI : 10.1146/annurev.fluid.33.1.445**

Henshaw M.J. de C
2007
“Non-linear aeroelastic prediction for aircraft applications”
Progress in Aerospace Sciences
43
(4-6)
65 -
137
** DOI : 10.1016/j.paerosci.2007.05.002**

Katz J.
,
Plotkin A.
2001
Low-Speed Aerodynamics
Cambridge University Press

Wang Y.
,
Wan Z.
,
Yang C.
2012
“Application of highorder panel method in static aeroelastic analysis of aircraft”
Procedia Engineering
31
136 -
144
** DOI : 10.1016/j.proeng.2012.01.1003**

Yang C.
,
Zhang B.C.
,
Wan Z.Q.
,
Wang Y.K.
2011
“A method for static aeroelastic analysis based on the highorder panel method and modal method”
Science China Technological Sciences
54
(3)
741 -
748
** DOI : 10.1007/s11431-010-4253-4**

Newman J.C. III
,
Newman P.A.
,
Taylor A.C. III
,
Hou G.J.W.
1999
“Efficient nonlinear static aeroelastic wing analysis”
Computers & Fluids
28
(4)
615 -
628
** DOI : 10.1016/S0045-7930(98)00047-4**

Bathe K.J.
1996
Finite Element Procedures
Prentice Hall
Upper Saddle River, New Jersey

Kapania K.
,
Raciti S
1989
“Recent advances in analysis of laminated beams and plates, part II: Vibrations and wave propagation”
AIAA Journal
27
(7)
935 -
946
** DOI : 10.2514/3.59909**

El Fatmi R.
2007
“Non-uniform warping including the effects of torsion and shear forces. Part I: A general beam theory”
International Journal of Solids and Structures
44
(18-19)
5912 -
5929
** DOI : 10.1016/j.ijsolstr.2007.02.006**

Schardt R.
1966
“Extension of the engineer’s theory of bending to the analysis of folded plate structures”
Der Stahlbau
35
161 -
171

Silvestre N.
,
Camotim D.
2002
“Second-order generalised beam theory for arbitrary orthotropic materials”
Thin-Walled Structures
40
(9)
791 -
820
** DOI : 10.1016/S0263-8231(02)00026-5**

Berdichevsky V.L.
,
Armanios E.
,
Badir A.
1992
“Theory of anisotropic thin-walled closed-cross-section beams”
Composites Engineering
2
(5-7)
411 -
432
** DOI : 10.1016/0961-9526(92)90035-5**

Yu W.
,
Volovoi V.V.
,
Hodges D.H.
,
Hong X.
2002
“Validation of the variational asymptotic beam sectional analysis (VABS)”
AIAA Journal
40
(10)
2105 -
2113
** DOI : 10.2514/2.1545**

Patil M.J.
“Aeroelastic tailoring of composite box beams”
Proceedings of the 35th Aerospace Sciences Meeting and Exhibit, AIAA Paper 97-0015
Reno, NV
January 1997

Librescu L.
,
Song O.
1992
“On the static aeroelastic tailoring of composite aircraft swept wings modelled as thinwalled beam structures”
Composites Engineering
2
(5-7)
497 -
512
** DOI : 10.1016/0961-9526(92)90039-9**

Qin Z.
,
Librescu L.
2003
“Aeroelastic instability of aircraft wings modelled as anisotropic composite thinwalled beams in incompressible flow”
Journal of Fluids and Structures
18
(1)
43 -
61
** DOI : 10.1016/S0889-9746(03)00082-3**

Jeon S.M.
,
Lee I.
2001
“Aeroelastic response and stability analysis of composite rotor blades in forward flight”
Composites: Part B
32
(1)
249 -
257
** DOI : 10.1016/S1359-8368(00)00061-5**

Friedmann P.P.
,
Glaz B.
,
Palacios R.
2009
“A moderate deflection composite helicopter rotor blade model with an improved cross-sectional analysis”
International Journal of Solids and Structures
46
(10)
2186 -
2200
** DOI : 10.1016/j.ijsolstr.2008.09.017**

Carrera E.
,
Giunta G.
,
Petrolo M.
2011
Beam Structures: Classical and Advanced Theories
John Wiley & Sons

Varello A.
,
Carrera E.
,
Demasi L.
2011
“Vortex lattice method coupled with advanced one-dimensional structural models”
Journal of Aeroelasticity and Structural Dynamics
2
(2)
53 -
78

Varello A.
,
Demasi L.
,
Carrera E.
,
Giunta G.
“An improved beam formulation for aeroelastic applications”
Proceedings of the 51st AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference
Orlando, FL
12-15, April 2010
AIAA Paper 2010-3032

Carrera E.
,
Varello A.
,
Demasi L.
2013
“A refined structural model for static aeroelastic response and divergence of metallic and composite wings”
CEAS Aeronautical Journal
4
(2)
175 -
189
** DOI : 10.1007/s13272-013-0063-2**

Petrolo M.
2013
“Flutter analysis of composite lifting surfaces by the 1D Carrera Unified Formulation and the doublet lattice method”
Composite Structures
95
539 -
546
** DOI : 10.1016/j.compstruct.2012.06.021**

Maskew B.
1987
“Program VSAERO theory Document: a computer program for calculating nonlinear aerodynamic characteristics of arbitrary configurations”
NASA contractor report CR-4023

Oosthuizen P.H.
,
Carscallen W.E.
1997
Compressible fluid flow
McGraw-Hill

Deperrois A.
2011
Guidelines for XFLR5 v6.03

Carrera E.
,
Varello A.
2012
“Dynamic response of thin-walled structures by variable kinematic onedimensional models”
Journal of Sound and Vibration
331
(24)
5268 -
5282
** DOI : 10.1016/j.jsv.2012.07.006**

Citing 'Static Aeroelastic Response of Wing-Structures Accounting for In-Plane Cross-Section Deformation
'

@article{ HGJHC0_2013_v14n4_310}
,title={Static Aeroelastic Response of Wing-Structures Accounting for In-Plane Cross-Section Deformation}
,volume={4}
, url={http://dx.doi.org/10.5139/IJASS.2013.14.4.310}, DOI={10.5139/IJASS.2013.14.4.310}
, number= {4}
, journal={International Journal of Aeronautical and Space Sciences}
, publisher={The Korean Society for Aeronautical & Space Sciences}
, author={Varello, Alberto
and
Lamberti, Alessandro
and
Carrera, Erasmo}
, year={2013}
, month={Dec}