Numerical Modelling of Static Aeroelastic Deformations of Slender Wing in Aerodynamic Design

The VLM determines wing lifting force and induced drag. The VLM is based on the potential flow theory. By taking into account the Prandtl–Glauert correction, this method has a satisfactory level of accuracy for flight speeds with a Mach number lower than 0.6. The method requires only a coarse definition of the aircraft geometry and the flight state, the simplest general method for 3D aerodynamic analyses of aircraft. Due to few input parameters, analyses can be set up with little effort, and analyses are computationally inexpensive. A detailed formulation of the steady-state and the unsteady variants of VLM can be found in Ref. [13].

The AVL program has been implemented into the aeroelastic platform XAVEL. AVL was created by Mark Drela from MIT Aero & Astro and Harold Youngren [14] and is an open-source program for the aerodynamic and flight-dynamic analysis of rigid aircraft of arbitrary configuration. It employs an extended vortex lattice model for the lifting surfaces, together with a slender-body model for fuselages and nacelles. Extensions of the classic VLM in AVL include accounting for profile drag, calculated in wing sections based on local induced angles of attack and local profile characteristics obtained by other means, experimental or computational.

Another module used in the XAVEL calculations is X-Foil. The X-Foil is an open-source program created by Drela [15] to estimate profile aerodynamic coefficients, including viscous drag. With the assumption of two-dimensional flow and modelling of boundary layer flow, XFOIL allows for the computation of a total airfoil drag value for the range of angles of attack where flow is attached up to the occurrence of moderate regions of flow separation. This translates into much more accurate predictions of wing performance parameters such as lift-to-drag ratio or energy function.

In the proposed approach to the aeroelastic analysis of the wing, the decision was made to use a simplified structural model—a beam model. Beam models used in numerical calculations are usually based on two popular theories—the classical Euler–Bernoulli beam theory and the Timoshenko theory.

An Euler–Bernoulli beam model assumes that sections perpendicular to the beam’s axis remain as such after beam deformation, and it is suitable for medium and high-aspect-ratio beams (aspect ratio higher than 8).

The Euler–Bernoulli beam model is well described in the literature [16-19]. The relationships describing bending and torsional deformations formed the basis for the numerical implementation of the beam model using the C programing language. The finite difference method is used for the beam model’s numerical implementation.

The low-fidelity method uses a two-way interaction between the aerodynamic and structural models. Loading the wing with an aerodynamic force and moment causes the wing to deform, changing the aerodynamic load. Therefore, such a problem must be solved iteratively until the load and deformation match. The implementation of purely one-way coupling often misses the point, as the twisting of the wing directly affects the angle of attack and often changes the wing loading significantly. Figure 1 shows the diagram of two-way fluid-structure interaction implemented in XAVEL.

Figure 1.

There are many ways to monitor the convergence of the aforementioned iterative process. In the presented case, the changes of wing angle of attack in successive iterations are observed. A change value below the assumed threshold ends the iterative process. Maintaining a constant lift force coefficient rather than the angle of attack (which is a typical case) has a huge advantage with this type of aeroelastic analysis. It allows to balance the aircraft’s weight and precisely determine envelope points and corresponding deformations.

The aerodynamic model described in the previous chapter, based on the VLM method, is restricted to moderate angles of attack and low-to-moderate Mach numbers, which is a serious applicability constraint. For this reason, a more advanced model involving solutions of Reynolds-Averaged Navier–Stokes equations and mesh morphing was prepared to be used with a model for computing deformations of a wing structure. The wider applicability range of numerical solutions of Navier–Stokes equations has a price, being the dependence on a computational mesh for flow volume, which needs to follow the deformations of the solid structure. In the presented work, the main effort on the CFD model concerned the method of modification of the computational mesh in the vicinity of the deformed solid structure. For this purpose, the mesh deformation module using RBF functions was implemented in the ANSYS Fluent software using Dynamic Mesh macros available in the UDF’s library. This necessitated the application of the Unsteady Navier–Stokes equations (URANS) solver, although the computational cases used to validate the method consisted of purely steady-flow.

The applied mesh-deformation technique, exploiting Radial Basis Functions, belongs to a category of mesh morphing, which involves morphing of the whole or part of the computational domain following the deformation of its boundary. The most important feature of mesh morphing is that it does not change the number of mesh elements nor the internal structure of the mesh, which means that the structure of the dense layers of elements in the wing boundary layer is preserved after the deformation. In the presented method, morphing is limited to a part of the computational domain, limited by a surface made of points where a mesh deformation coefficient assumes values of zero. The most-often used competing techniques of mesh modification in vicinity of deforming boundary belong to the re-meshing category, which involves destroying parts of the existing mesh and creating new elements. These techniques are used most often in the modelling of store separation and modification of mesh when the domain boundary undergoes discontinuous changes, like the opening of aircraft internal bays. Another technique that could be used in the modelling of deformation of the domain boundary is overset, which involves two meshes in the same domain: the background nondeformable mesh and deforming mesh of smaller size, covering only the vicinity of the deforming boundary. Exchange of information between the two meshes during the solution is conducted in overset interfaces created by the solver. The interfaces are modified if the geometry of the deformable zone changes during the solution. This approach has not been applied in the present work, but is considered for future development of the method.

The RBF re-meshing process involves the following steps:

Define source points: Select a set of source points in the domain that will be used to guide the deformation. These points can be located on the surface of the boundary (wing surface) and within the mesh volume, at a distance from the boundary (e.g., to limit the morphing region).

The RBF functions are suitable for interpolation and extrapolation of the deformations. A set of discrete points where deformation is known must be defined. A scalar RBF interpolating/extrapolating function S is defined as: 1S(x)=∑i=1Nηiφ(∥x−xi∥)+h(x)$$S(x) = \sum\nolimits_{i = 1}^N {{\eta _i}} \varphi \left( {\parallel x - {x_i}\parallel } \right) + h(x)$$ where: N is the number of source points with known deformations, φ is a radial function from a selected basis, being a transformation ℝ → ℝ, x = [x, y, z]^T is a point where S is computed, x_i = [x_i, y_i, z_i]^T is a source point where deformation is known, η_i is an element of a computed vector of coefficients η, and h is a correction polynomial.

The correction polynomial of the order m–1, where m, an order of the φ function, is required for the uniqueness of the solution and enables rigid motion of the structure. The interpolating function belonging to basis φ exists if the coefficients η and elements of polynomial h can be found such that the following two conditions are fulfilled: 2S(x)=G(xi),1≤i≤N$$S(x) = G\left( {{x_i}} \right),1 \le i \le N$$ where G(x_i) is a known value at point x_i. In the presented method, G(x_i) is a component of boundary deformation, namely translation of point x_i. along each one of the coordinate axes.

The second condition is an orthogonality condition for all polynomials p with a degree less or equal to the degree of the polynomial h: 3∑i=1Nηip(xi)=0$$\sum\nolimits_{i = 1}^N {{\eta _i}} p\left( {{x_i}} \right) = 0$$

If the order of function φ is m ≤ 2, then in the three-dimensional space, a polynomial h may be applied, of the form: 4h(x)=β1+β2x+β3y+β4z$$h(x) = {\beta _1} + {\beta _2}x + {\beta _3}y + {\beta _4}z$$

The values of coefficients η and β may be determined from the solution of the system of order N: 5(UPPT0)(ηβ)=(G0)$$\left( {\matrix{ U & P \cr {{P^T}} & 0 \cr } } \right)\left( {\matrix{ \eta \hfill \cr \beta \hfill \cr } } \right) = \left( {\matrix{ G \hfill \cr 0 \hfill \cr } } \right)$$ where: G = {G₁, …, G_N}^T are known values of interpolation function in source points, and β = {β₁, …,β_N}^T are coefficients of polynomial h.

U is interpolation matrix N × N, dependent on the distances between points: Uij=φ(∥xi−xj∥),1≤i≤N,1≤j≤N;$${U_{ij}} = \varphi \left( {\parallel {x_i} - {x_j}\parallel } \right),1 \le i \le N,1 \le j \le N;$$

P is a constraint matrix of the form: 6P=(1x1y1z11x2y2z2⋮⋮⋮⋮1xNyNzN).$${\bf{P}} = \left( {\matrix{ 1 & {{x_1}} & {{y_1}} & {{z_1}} \cr 1 & {{x_2}} & {{y_2}} & {{z_2}} \cr \vdots & \vdots & \vdots & \vdots \cr 1 & {{x_N}} & {{y_N}} & {{z_N}} \cr } } \right).$$

In the described implementation of RBF functions, the interpolated variables (scalars S in Eq. (1)) are components of transitions of the nodes of the computational mesh: 7dx=Sx(x)=∑i=1Nηxiφ(∥x−xi∥)+βx1+βx2x+βx3y+βx4zdy=Sy(x)=∑i=1Nηyiφ(∥x−xi∥)+βy1+βy2x+βy3y+βy4zdz=Sz(x)=∑i=1Nηziφ(∥x−xi∥)+βz1+βz2x+βz3y+βz4z.$$\matrix{ {{d_x} = {S_x}(x) = \sum\nolimits_{i = 1}^N {{\eta _{xi}}} \varphi \left( {\parallel x - {x_i}\parallel } \right) + {\beta _{x1}} + {\beta _{x2}}x + {\beta _{x3}}y + {\beta _{x4}}z} \hfill \cr {{d_y} = {S_y}(x) = \sum\nolimits_{i = 1}^N {{\eta _{yi}}} \varphi \left( {\parallel x - {x_i}\parallel } \right) + {\beta _{y1}} + {\beta _{y2}}x + {\beta _{y3}}y + {\beta _{y4}}z} \hfill \cr {{d_z} = {S_z}(x) = \sum\nolimits_{i = 1}^N {{\eta _{zi}}} \varphi \left( {\parallel x - {x_i}\parallel } \right) + {\beta _{z1}} + {\beta _{z2}}x + {\beta _{z3}}y + {\beta _{z4}}z.} \hfill \cr } $$

The node translations described above depend on translations defined in source points, which need not be mesh nodes. It makes the method mesh-independent, which feature in Refs [20, 21], and is called ‘meshless’. The other important feature of the method is the possibility of application in parallel computations on computer nodes when each computational thread has access to only a part of the computational mesh. It requires only sending to each computational process tables of source points, and coefficients η and β which can be determined in the host process of the computations.

Two groups of source points are defined in the described implementation of RBF functions. The first group consists of points on the wing surface and the other group consists of points at some distance from the surface, where it is desired to end the deformable region of the mesh and keep the remaining region intact. This is shown in Fig. 2, where the dense distribution of source points on the wing surface in the center is surrounded by two layers of source points at some radial distance, in order to end the deformable zone.

Figure 2.

The source points, which are located on the wing surface, are distributed in chordwise rows, with density increasing at leading and trailing edges, as shown in Fig. 3. In the conducted analyses, 20–50 points in a row were used, and the number of rows was 10–20.

Figure 3.

The mesh deformations produced by wing deformation must vanish in the space between the layers of the source points external to the wing, visible in Fig. 3. This is accomplished by introducing an additional mesh deformation coefficient d. The function defining its distribution is given in the cylindrical coordinate system with a longitudinal axis along span and spherical at the wing tip: 8d=1ifr→∈(0,r→0),d=0.5⋅(1+cos(r→−r→0r→1−r→0⋅π))ifx→∈(r→0,r→1),d=0ifr→>r→1,$$\matrix{ {d = 1{\rm{if}}\vec r \in \left( {0,{{\vec r}_0}} \right),} \hfill \cr {d = 0.5 \cdot \left( {1 + \cos \left( {{{\vec r - {{\vec r}_0}} \over {{{\vec r}_1} - {{\vec r}_0}}} \cdot \pi } \right)} \right){\rm{if}}\vec x \in \left( {{{\vec r}_0},{{\vec r}_1}} \right){\rm{,}}} \hfill \cr {d = 0{\rm{if}}\vec r > {{\vec r}_1},} \hfill \cr } $$ where r₀ and r₁ are the radii of the two layers of external source points.

The final mesh deformation in a point in the flow domain is the product of the deformation computed with Eq. (7) and the mesh deformation coefficient defined by Eq. (8). This assures the mesh deformation is continuous and decreases to zero in the desired region. The volume where the value of d is higher than zero is thus a part of the computational domain subject to morphing.

The spatial distribution of coefficient d computed for an exemplary rectangular wing of aspect-ratio 5 is shown in Fig. 4. For demonstration purposes, the wing has been subject to a deformation resulting in wing-tip translation in the normal direction of 5% span. The spatial distribution of node translation in the direction normal to the wing planform of the same wing is shown in Fig. 5 as isolines of constant translation. The translations are presented in two perpendicular planes: one plane along the wing span and another along the wing chord, in approximately 0.7 half-span. In Fig. 6 is presented the crosssection of the domain in a plane located in the wing half-chord, and presenting details of the mesh. The resulting distribution of node translations shows that the goal of ensuring continuous, smooth distribution of node deformation, vanishing at a prescribed distance from the deforming structure, has been reached.

Figure 4.

Figure 5.

Figure 6.

The mesh deformation is initiated by the import to the solver of the set of source points located on the deformed surface. Then the system of Eq. (5) is solved in the host process of ANSYS Fluent in order to determine the vectors of coefficients η and β. These vectors are then sent to the parallel processes operating on fragments of the domain. The deformation, described by Eqs. (7) and (8), is conducted by the solver executing a user-defined function using the macro DEFINE_GRID_MOTION of ANSYS Fluent solver during the prescribed number of time steps.

A finite-element beam model, described in Ref. [23], was used to discretize the structure of the wing and to determine the deformations under aerodynamic loads. The beam model includes the following parameters at the nodes of the elements: –

longitudinal displacements u₁, u₂,

The stiffness matrix of an element has size 12 by 12, and its elements are dependent on tensile stiffness, EA, bending stiffness in two planes: EJ_y, EJ_z, and torsional stiffness, GK, which may be determined analytically, e.g., by methods described in [23]. The stiffness matrix of the beam element has the form [21]: 9[ EAl012EJzl30012EJyl3000GKsl00004EJyl000004EJzl−EAl00000EAlSymmetry0−12EJzl3000−6EJzl2012EJzl300−12EJyl30000012EJyl3000−GKsl00000GKsl00−6EJyl202EJyl000004EJyl06EJzl20002EJzl0−6EJzl20004EJzl ]$$\left[ {\matrix{ {{{EA} \over l}} \hfill & {} \hfill & {} \hfill & {} \hfill & {} \hfill & {} \hfill & {} \hfill & {} \hfill & {} \hfill & {} \hfill & {} \hfill & {} \hfill \cr 0 \hfill & {{{12E{J_z}} \over {{l^3}}}} \hfill & {} \hfill & {} \hfill & {} \hfill & {} \hfill & {} \hfill & {} \hfill & {} \hfill & {} \hfill & {} \hfill & {} \hfill \cr 0 \hfill & 0 \hfill & {{{12E{J_y}} \over {{l^3}}}} \hfill & {} \hfill & {} \hfill & {} \hfill & {} \hfill & {} \hfill & {} \hfill & {} \hfill & {} \hfill & {} \hfill \cr 0 \hfill & 0 \hfill & 0 \hfill & {{{G{K_s}} \over l}} \hfill & {} \hfill & {} \hfill & {} \hfill & {} \hfill & {} \hfill & {} \hfill & {} \hfill & {} \hfill \cr 0 \hfill & 0 \hfill & 0 \hfill & 0 \hfill & {{{4E{J_y}} \over l}} \hfill & {} \hfill & {} \hfill & {} \hfill & {} \hfill & {} \hfill & {} \hfill & {} \hfill \cr 0 \hfill & 0 \hfill & 0 \hfill & 0 \hfill & 0 \hfill & {{{4E{J_z}} \over l}} \hfill & {} \hfill & {} \hfill & {} \hfill & {} \hfill & {} \hfill & {} \hfill \cr {{{ - EA} \over l}} \hfill & 0 \hfill & 0 \hfill & 0 \hfill & 0 \hfill & 0 \hfill & {{{EA} \over l}} \hfill & {} \hfill & {} \hfill & {{\rm{Symmetry}}} \hfill & {} \hfill & {} \hfill \cr 0 \hfill & {{{ - 12E{J_z}} \over {{l^3}}}} \hfill & 0 \hfill & 0 \hfill & 0 \hfill & {{{ - 6E{J_z}} \over {{l^2}}}} \hfill & 0 \hfill & {{{12E{J_z}} \over {{l^3}}}} \hfill & {} \hfill & {} \hfill & {} \hfill & {} \hfill \cr 0 \hfill & 0 \hfill & {{{ - 12E{J_y}} \over {{l^3}}}} \hfill & 0 \hfill & 0 \hfill & 0 \hfill & 0 \hfill & 0 \hfill & {{{12E{J_y}} \over {{l^3}}}} \hfill & {} \hfill & {} \hfill & {} \hfill \cr 0 \hfill & 0 \hfill & 0 \hfill & {{{ - G{K_s}} \over l}} \hfill & 0 \hfill & 0 \hfill & 0 \hfill & 0 \hfill & 0 \hfill & {{{G{K_s}} \over l}} \hfill & {} \hfill & {} \hfill \cr 0 \hfill & 0 \hfill & {{{ - 6E{J_y}} \over {{l^2}}}} \hfill & 0 \hfill & {{{2E{J_y}} \over l}} \hfill & 0 \hfill & 0 \hfill & 0 \hfill & 0 \hfill & 0 \hfill & {{{4E{J_y}} \over l}} \hfill & {} \hfill \cr 0 \hfill & {{{6E{J_z}} \over {{l^2}}}} \hfill & 0 \hfill & 0 \hfill & 0 \hfill & {{{2E{J_z}} \over l}} \hfill & 0 \hfill & {{{ - 6E{J_z}} \over {{l^2}}}} \hfill & 0 \hfill & 0 \hfill & 0 \hfill & {{{4E{J_z}} \over l}} \hfill \cr } } \right]$$

The rows and columns of the stiffness matrix correspond to the node parameters described above. They may be interpreted as forces and moments created by unit displacements or rotations in the respective nodes.

The stiffness matrix of the beam is composed of matrices for single elements by adding contributions of each element’s stiffness in the rows of the global matrix, corresponding to degrees of freedom in the nodes of the discretized wing. In order to take into account the geometric orientation of elements (a wing may be segmented, with varying dihedral, sweep, twist), each element matrix is first transformed to a beam coordinate system by multiplication by transformation matrices composed of directional cosines of element axis [24].

The finite element method also requires that continuous loads on wing surfaces are reduced to discrete loads in the beam nodes, corresponding to the node degrees of freedom. This is done based on the principle of equality of work of external and internal forces on virtual displacements [24].

Application of boundary conditions in the present case assumes zero nodal linear and angular deflection in particular nodes and directions, e.g., zero deflection and rotation at the wing root. It means that the corresponding rows and columns of the global stiffness matrix are removed from the system.

The wing deformation, consisting of translations and rotations in the element nodes, are determined by the solution of the system: 10{X}=[K]−1⋅{F}$$\{ X\} = {[K]^{ - 1}} \cdot \{ F\} $$ where {X}, deformation vector; [K]^–1, inverse stiffness matrix; and {F}, nodal force vector.

Depending on the problem, the analysis may be conducted in a one-step fashion, by the one-time solution of the system (1), or, in order to account for the nonlinear behavior of the structure, by performing the solution in more steps. The first option is suitable for relatively stiff wings and moderate deformations. In many cases, e.g., high aspect-ratio wings undergoing large deflections, it is necessary to account for nonlinear effects in the deformations. This is conducted in the presented method by application of the predictor–corrector procedure. The stiffness matrix is computed for the nondeformed beam in the predictor step. The system (10) is solved, the deformations are halved, and applied in the nodes, in order to compute the stiffness matrix of the deformed structure for the corrector step. Then the system (10) is solved again, and the deformations are applied to the original, nondeflected beam. As a final step in the structural analysis, the nodal parameters computed for the nodes on the wing (beam) axis are used to compute the displacements in the RBF source points located on the wing surface. This process includes translations of the RBF points along beam axes and rotations of the wing cross-sections as a result of wing bending and twist. This is shown in Fig. 7.

Figure 7.