(1)Three-dimensional aeroelasticity
program for finite thin-wings in incompressible and steady air flow
by
Christophe Paget
and
Jérôme Matrat
Introduction
A wing in airflow faces an aerodynamic pressure applied on a structural wing. Those two different fields constitute the aeroelasticity. In the world-wide literature, aeroelasticity was modelled in several different way, such as aerodynamic forces applied in basic 1D or 2D mechanical structures. In those cases, it is difficult to view the deformation field of such wing. Therefore the idea was to create a program which combines both aerodynamic part and 3D structural part. For this latter, a finite element method was used. As first approach, the program will be studying steady flow. For educational purposes, the program fully uses the graphic user interfaces. The program has been written in Matlab code combining a Fortran executable for the finite element code.
1. Aerodynamic theory
Aeroelasticity, as the name explains, is the combination of aerodynamic and elasticity theories. Let us start with the aerodynamics. In the aerodynamic point of view, the flow over the finite wing was considered as incompressible and steady. For incompressible and steady flow, several theories were developed in the early sixties, such as Kernel function method (KFM) [1], Panel Method [2] or even Theodorsen method [3].
The method used in this program is the panel method. Several methods were available for 3-D models, such as Lifting-line solution by Horseshoe elements, or Lifting-surface solution by Vortex ring elements. For cambered wings, the most appropriate method was the Lifting-surface solution by Vortex ring elements.
As the name explains, the theory involves vortex ring elements, which have the same shape as the top surfaces of the elements of the top layer. The main advantage of this element is in the simple programming effort that it requires (although its computational efficiency can be further improved). Additionally, the exact boundary conditions will be satisfied on the actual wing surface, which can have camber and various planform shapes.
The singularity element is based on the vortex line solution of the incompressible continuity equation. The boundary condition that must be satisfied by the solution is the zero normal flow across the thin wing’s solid surface.
(1)
In order to solve the lifting-surface problem numerically, the wing is divided into elements containing vortex ring singularities. The solution procedure is as follows:
1.1. Induced velocity formulation
The velocity induced by such a vortex line of circulation G , according to the Biot-Savart law is:
(2)
where dl is the infinitesimal portion of the filament.
For a straight vortex line segment (between the points P1 and P2) at a point P, the velocity induced becomes:
(3)
where
(4)
(5)
(6)
(7)
(8)
(9)
In the program, the function was call VORTXL, such as:
[u,v,w]=VORTXL(xp, yp, zp, x1, y1, z1, x2, y2, z2, G ) (10)
To obtain the velocity induced by the four segments (Figure 1) of the parallelepiped vortex ring with circulation G , the velocity induced by each segment was defined as follow:
|
[u1,v1,w1]=vortxl(xp, yp, zp, x1, y1, z1, x2, y2, z2, G ) |
|
|
[u2,v2,w2]=vortxl(xp, yp, zp, x2, y2, z2, x3, y3, z3, G ) |
|
|
[u3,v3,w3]=vortxl(xp, yp, zp, x3, y3, z3, x4, y4, z4, G ) |
|
|
[u4,v4,w4]=vortxl(xp, yp, zp, x4, y4, z4, x1, y1, z1, G ) |
|
|
|
|
|
Then, |
|
|
(u,v,w)=(u1,v1,w1)+(u2,v2,w2)+(u3,v3,w3)+(u4,v4,w4) |
|
|
|
|
|
In the program, this routine was called VORING. |
|
|
Figure 1: Vortex ring element |
1.2. Choice of singularity element
The method by which the thin-wing planform is divided into panels is shown in Figure 2. The leading segment of the vortex ring is placed on the panel’s quarter chord line and the collocation point is at the centre of the three-quarter chord line, as shown in Figure 2. The normal vector n is defined at this point, too. A positive G is defined here according to the right-hand rotation rule (for the leading segment).
Figure 2: Nomenclature for the vortex ring elements
From the numerical point of view these vortex ring elements are stored in rectangular patches (arrays) with i, j indexing as shown by Figure 2. The velocity induced at an arbitrary point P(x ,y , z), by a typical vortex ring at the location i, j was computed under the function:
[u, v, w]=VORING(x, y, z, i, j, G )
The use of this subroutine can considerably shorten the programming effort; however, for the vortex segment between two such rings the induced velocity is computed twice. For the sake of simplicity this routine will be used for this problem, but more advanced programming can easily correct this loss of computation effort.
At this stage, the velocity induced by the trailing vortex segments only (the vortex lines parallel to the free-stream). This information is needed for the induced drag computations and if done at this phase will only slightly increase the computational effort. The influence of the trailing segments is obtained by simply omitting the (u1,v1,w1)+(u3,v3,w3) part from the final expression of the induced velocity.
Therefore the velocity induced by the trailing vortex segments is:
(u,v,w)*=(u2,v2,w2)+(u4,v4,w4) (11)
1.3. Discretisation and grid generation
Because only the wing semispan is modelled, the mirror-image method was used to account for the other semispan. For the pressure distribution calculations the local circulation is needed, which for the leading edge panel is equal to G i but for all the elements behind it is equal to the difference G i-G i-1. In the case of increasing surface curvature the above-described vortex rings will not be placed exactly on the lifting surface, and a finer grid needs to be used, or the wing surface can be redefined accordingly. By placing the leading segment of the vortex ring at the quarter chord line of the panel the two-dimensional Kutta condition is satisfied along the chord (confer the lumped-vortex element [2]). Also, along the wing trailing edges, the trailing vortex of the last panel row (which actually simulates the starting vortex) must be cancelled to satisfy the three-dimensional trailing edge condition:
g
T.E.=0 (12)For steady-state flow this is done by attempting to align the wake vortex panels parallel to the local streamlines, and their strength is equal to the strength of the shedding panel at the trailing edge (G T.E. = G W for each row).
For the example of Figure 2, the chord is divided equally into M=3 panels and the semi-span is divided equally into N=3 panels. Therefore, the chordwise counter i will have values from 1 è M and the spanwise counter j will have values between 1 è N. Also, geometrical information such as the vortex ring corner points, panel area SK, normal vector nk, and the coordinates of collocation points are calculated at this phase (note the panel sequential counter i will have values between 1 and M x N). A simple and fairly general method for evaluating the normal vector is shown in Figure 3. The panel opposite corner points define two vectors Ak and Bk and their vector product will point in the direction of nk:
|
|
|
|
The results of the grid generating |
|
|
phase are shown schematically in |
|
|
Figure 3. For more information |
|
|
about generating panel corner points, |
|
|
collocation points, area and normal |
|
|
vector, refer to the subroutine PANEL. |
|
|
Figure 3: Definition of the wing outward normal |
1.4. Influence coefficient
Let us assume that the collocation point scanning has started at i=1 and j=1. The velocity induced by the first vortex ring is then
[ui, vi, wi]11=VORING(x, y, z, i=1, j=1, G =1.0)
and from its image on the other semispan
[uii, vii, wii]11=VORING(-x, y, z, i=1, j=1, G =1.0)
and the velocity induced by the unit strength G 1 and its image at collocation point 1 is:
(u,v,w)11=(ui-uii,vi+vii,wi+wii) (14)
Note that the subscript ( )11 represents the influence of the first vortex at the first collocation point, and both counters can have values from 1 to M x N. Also, a unit strength vortex is used in the process of evaluating the influence coefficient a11, which is:
(15)
To scan all the vortex rings influencing this point, an inner scanning loop is needed with the counter L=1 è N x M. Thus, at this point, the K counter is at point 1, and the L counter will scan all the vortex rings on the wing surface, and all the influence coefficients a1L are computed (also, in eq. (14) the ( )11 index means K=1, L=1):
(16)
When a particular vortex ring is at the trailing edge, a "free wake" vortex ring with the same strength is added to cancel the spanwise starting vortex line. Therefore, when the influence of such a trailing-edge panel vortex is calculated, the contribution of this segment is added. For example, in Figure 2 the first wake panel is encountered when i=3 (or L counter is equal to 7). If the wake grid is added into the M+1 corner point array, then the velocity due to the i=3, j=1 (or L=7) panel is
[u, v, w]19=VORING(x1, y1, z1, i=3, j=1, G =1.0)
and due to the attached wake
[u, v, w]19W=VORING(x1, y1, z1, i=3+1, j=1, G =1.0)
When the wing is symmetric as in this case and only the right half wing is panelled, then the 
velocity components of the trailing edge and wake panels include the influence of the left-hand side image (as in eq. (14)). The corresponding influence coefficient is
(17)
As mentioned before, parallel to the computation of the aKL coefficients, the normal velocity component induced by the streamwise segment can also be computed by using the 
portion as in eq. (11). For the first element then
(18)
This procedure continues until all the collocation points have been scanned.
1.5. Establish the right-hand side vector (RHS)
The RHS vector is computed as before by scanning each of the collocation points on the wing:
(19)
where Q¥ has its components in the axis attacked to the wing:
(20)
1.6. Solve linear set of equations
Once the computations of the influence coefficients and the right-hand side vector are completed, the zero normal flow boundary condition on each of the collocation points will result in the following set of algebraic equations:
(21)
Therefore
(22)
Here K is the vertical collocation point counter and L is the horizontal vortex ring counter and the order of this matrix is m = M x N.
The induced downwash is defined by the linear equation:
(23)
where again m = N x M.
1.7. Secondary computations
The solution of the above set of equations results in the vector (G 1, …, G K, …,G m). At this stage, the counter K is replaced back to the original i, j counters. The lift of each bound vortex segment is obtained by using the Kutta-Joukowski theorem:
i > 1 (24)
and when the panel is at the leading edge then
i = 1 (25a)
Therefore the lift coefficient is:
(25b)
where 
is the panel area.
The pressure difference across this panel is
(26)
where is the panel area.
The induced drag is defined in the same manner as the lift:
i > 1 (27)
i = 1 (28)
where the induced downwash at each collocation point i, j is computed by summing up the velocity induced by all the trailing vortex segments.
Therefore the drag coefficient is:
where is the panel area.
Further, the load condition of the finite element model is provided of the local pressure difference corresponding for each element.
2. Aerodynamic results
In order to validate the aerodynamic loads obtained with the lifting-surface solution by vortex ring element, some results were compared with test provided in the world-wide literature.
To start with, the effect of the aspect ratio on the lift coefficient slope of the untapered planar wing with zero sweep angle is compared with some experimental data [4] in Figure 4. In this example, the wing is a planar wing planform, where the leading, trailing and side edges are made of straight lines and the wing has no camber (NACA 0020).
The Figure 4 shows good agreements between the vortex ring element method and the experimental data.
Figure 4: Effect of the aspect ratio on the lift coefficient slope of untapered planar wing, for a zero sweep angle; Solid line corresponds to the present and the various dots to experimental data [ref. 4]
Figure 5: Effect of the drag coefficient on the lift coefficient of an untapered wing (NACA 23004), for AR=1; Solid line corresponds to AeroWing calculations and the dots to experimental data [ref. 5]
Figure 6: Effect of the lift coefficient on the angle of attack of an untapered wing (NACA 4415), for AR=12.07, Reynolds number of 6.0 x 106; Solid line corresponds to AeroWing calculations and the triangles to experimental data [ref. 6]
The Figures 4 to 6 shows that the program AeroWing gives sufficient accuracy on the aerodynamics of finite wings.
3. Finite element code
In the scope of solving static equilibrium problems in small strain solid elasticity, a three-dimensional Finite Element program, so-called 20nodes, was done. Including twenty noded isoparametric elements described in the three dimensional space, this Fortran code calculates the deform shape, strains and stresses from the pressure distribution previously evaluated on a finite wing subjected to an incompressible and steady air flow.
The following will describe the code after a small introduction of the FE theory. Then, graphs will be presented to validate the program comparatively to an in-house FE code developed at FFA for several types of linear elastic problems.
FE formulation of three-dimensional elasticity
Considering an arbitrary part of a body, the forces acting on the part are given by the traction vector t along the boundary surface S and the body force b in the region V. Equilibrium requires that
Inserting the stress tensor in the expression of the forces on the external surface S gives
Using Gauss’ divergence theorem, one obtains
, for arbitrary volumes V
Thus, 
, which can be written in the compact matrix form :
(1)
where s contains all the stress components and Ñ is a matrix differential operator
Derive (1), multiply it by a virtual displacement v and integrate it over the volume V to perform integration by parts using the Green-Gauss theorem.
Finally the expression reduces to
which stands for the virtual work equation (2)
This weak form of the equilibrium equations is straightforward to formulate the FE equations for three-dimensional elasticity.
In the FE technique, the continuous variable u, i.e. the displacements, is approximated in terms of its nodal values, through functions of the space variable called shape functions, N. Consequently, the virtual displacement v can be expressed
and
where 
which, by inserting those expressions in (2), yields
Thus,
(3)
For elastic materials, 
, where D is the constitutive matrix.
Under the assumption of small displacement gradients, the strains can be derived as shown in the following equation
Then, (3) becomes
(4)
The boundary conditions can be also defined as
on Sf
and 
on Sg
where the load vector f is known along the boundary Sf and the displacement vector g is known along the boundary Sg.
Assuming in addition that we only consider external forces on the body, (4) takes the form
(5)
where KM is called the stiffness matrix.
The goal of the finite element method is to numerically solve this system of equations.
Isoparametric 20 noded Brick element
In this paragraph, the finite element formulation of a 20 noded brick element is examined.
Element definition
Two coordinate schemes are used in the element formulation:
The global coordinate system (x,y,z) defining the whole finite element model and the local coordinate system for each element (x,h,c).
Shape functions or interpolation functions are associated with each node so that the assumed displacement field and the global coordinates are defined at any point (x,h,c) within the element by the expressions:
Strain-displacement matrix B
The strains are defined in terms of the nodal displacements and shape functions derivatives as demonstrated previously in the paragraph "FE formulation of three-dimensional elasticity", by the expression:
Stiffness matrix evaluation
We now have all of the information necessary to calculate the element stiffness matrix KM where, as previously shown,
but dxdydz = dxdhdc.detJ
with 
and
J is the Jacobian matrix.
Then, the following transformation of integral is obtained
Numerical integration
In practice, Gauss-Legendre integration is used almost exclusively within isoparametric FE formulations. With this technique, weights and also positions of the integration points are determined so that a given polynomial is integrated exactly.
In this problem, as 20noded bricks element were used, 3 integration points were chosen at positions xI with corresponding weight wi.
|
xi |
wi |
|
+(Ö15)/5 |
5/9 |
|
0 |
8/9 |
|
-(Ö15)/5 |
5/9 |
Then the three-dimensional Gaussian quadrature formulae is written as:
for an integration over the volume of one brick element, in its local coordinate system (x,h,c)
Consistent nodal forces
To represent a distributed pressure load on one face of the element in terms of nodal forces, the consistent nodal force formulae shown in (5) was used.
For a distributed pressure P, the force Fi on node I can be expressed:
The integral for Fi is calculated using the Gaussian integration, as seen in the previous paragraph.
Sparse matrices
The sparse matrices are a special class of matrices that contain a significant number of zero-valued elements. This property allows reducing the amount of memory required for data storage by storing only the nonzero elements of the matrix together with their indices.
Moreover by eliminating operations on zero elements, computation time is significantly improved.
Choleski back-substitution
Factorisation of the system matrix:
Since K is symmetric and positive definite, it can be transformed using Cholesky factorisation into K = RTR, where R is upper triangular. This further reduces the storage space since only R has to be stored.
RTy=p
Ru=y
Finite element program structure
Based on the implementation of the FE method as described in Smith and Griffiths "Programming the Finite Element Method", this program adopted a modular approach in the sense that separate subroutines perform the various main finite element operations.
Organisation of the program:
|
Program description |
Subroutine used |
|
INPUT AND INITIALISATION |
|
|
Store the nodes into array COORD Store the 20noded brick elements into array ELE Store the material data in array MSTIFF Store the boundary conditions in array NF Find the maximum bandwidth for each freedom to determine the size of vector KV |
|
|
Store the weights and sampling points for the Gauss-Legendre quadrature (Numerical integration) into array SAMP |
GAUSS |
|
ELEMENT STIFFNESS INTEGRATION AND ASSEMBLY |
|
|
For each element IE equals 1 to NE Store the nodes coordinates of element IE into array NCOORD |
|
|
Full vector G with the nodes degree of freedom |
FORMG |
|
Null the stiffness matrix KM |
NULL |
|
For each Gauss integration points IG,JG,KG equals 1 to 3 |
|
|
Form the 3D shape functions FUN and their derivatives DER in local coordinates |
FMQUA3 |
|
Calculate the Jacobian matrix JAC = DER * NCOORD |
MATMUL |
|
Invert the Jacobian matrix and calculate its determinant DET, JAC1 = JAC -1 |
TREEX3 |
|
Calculate the shape functions derivatives DERIV in global coordinates DERIV = JAC1 * DER |
MATMUL |
|
Form the BEE strain – displacement matrix |
FORMB3 |
|
Store the material data associated with the element IE in MDATA |
|
|
Form the matrix of elastic constants DEE |
FORMD3 |
|
DBEE = DEE * BEE |
MATMUL |
|
Transpose the BEE matrix in BT |
MATRAN |
|
BTDB = BT * DBEE |
MATMUL |
|
Calculate the coefficients of the summation QUOT = DET * (the Gaussian weights SAMP) |
|
|
BTDB = QUOT * BTDB |
MSMULT |
|
Form the new stiffness matrix KM = BTDB + KM for the element IE |
MATADD |
|
Store and assemble the global stiffness matrix as a sparse matrix in vector KV |
FSPARV |
|
NODAL LOADS CALCULATION |
|
|
Read the pressure distribution PRESS on element ILOW to element IHIGH Find the face NFACE of the element where the pressure is applied For each element IE equals ILOW to IHIGH Store the nodes coordinates of element IE face NFACE into array NCOORDSURF |
|
|
Null the vector NODALOAD |
NULVEC |
|
For each Gauss integration points IG,JG equals 1 to 3 |
|
|
Form the 2D shape functions FUNSURF and their derivatives DERSURF |
FMQUAD |
|
Calculate the Jacobian matrix JACSURF = DERSURF * NCOORDSURF |
MATMUL |
|
Calculate the Jacobian matrix determinant DET |
DET2 |
|
Calculate the coefficients of the summation QUOT = DET * PRESS * (the Gaussian weight SAMP) |
|
|
FUNSURF = QUOT * FUNSURF |
VSMULT |
|
Form the load vector NODALOAD = NODALOAD + FUNSURF, for element IE |
VECADD |
|
Store the global load vector LOADS for each nodes with a non zero degree of freedom LOADS = LOADS + NODALOAD |
|
|
EQUATION SOLUTION |
|
|
Solve the matrix system KV * u = LOADS for the displacement vector u, using the Choleski factorisation with sparse matrices. |
SPARIN SPABAC |
|
The resulting displacement are stored in vector LOADS and printed to the output file |
PRINTDISP |
|
RECOVER ELEMENT STRAINS AND STRESSES AT BRICK CENTROIDS |
|
|
For each element IE equals 1 to NE Store the nodes coordinates of element IE into array NCOORD |
|
|
Full vector G with the nodes degree of freedom |
FORMG |
|
Store the displacements from LOADS into ELD and incorporate the nodes locked, Which were not considered for the resolution of the matrix system. For one Gauss integration point at the brick centroid |
|
|
Form the 3D shape functions FUN and their derivatives DER in local coordinates |
FMQUA3 |
|
Calculate the Jacobian matrix JAC = DER * NCOORD |
MATMUL |
|
Invert the Jacobian matrix and calculate its determinant DET, JAC1 = JAC -1 |
TREEX3 |
|
Calculate the shape functions derivatives DERIV in global coordinates DERIV = JAC1 * DER |
MATMUL |
|
Form the BEE strain – displacement matrix |
FORMB3 |
|
Form the strain vector VEPS = BEE * ELD |
MVMULT |
|
Store the material data associated with the element IE in MDATA |
|
|
Form the matrix of elastic constants DEE |
FORMD3 |
|
Form the stress vector VSIGMA = DEE * VEPS |
MVMULT |
|
Store the strains and stresses at each element centroid into EPS and SIGMA |
|
|
Print the results to the output file |
PRINTRES |
Description of the input-file
The input-file for the program 20nodes is divided in several sequences or cards that are terminated by an X typed in column 1.
The two first lines of the input are used to label it and may contain arbitrary information.
TITLE
TEXT
The nodes and their coordinates are then typed in columns.
NODE
NN Xn Yn Zn
NN Node number
Xn X-coordinate for node N
Yn Y-coordinate for node N
Zn Z-coordinate for node N
X
ENDNODES
The following card defines the 20-node isoparametric volume elements, numbered as shown
VOL20
NE N1 N2 N3 N4 N5 N6 N7 N8 -1 MT
NE N9 N10 N11 N12 N13 N14 N15 N16 N17 N18 N19 N20
NE Element number
Ni I=1,8 Corner nodes numbered as shown previously
Nj j=9,20 Mid-edge nodes
MT 3-D material index used in the element
X
ENDELEM
The 3D elastic material data are defined in the next card.
VOLMA
MT E1 E2 E3 G12 G23 G13
MT P21 P31 P32
MT 3-D material index
E1, E2, E3 Young's moduli in the global coordinate system X, Y, Z
G12, G23, G13 Shear moduli in planes XY, YZ, XZ
P21, P31, P32 Poisson's ratios
X
Locked degrees of freedom are applied to nodes as boundary conditions in card LOCK.
LOCK
ID NLOW NHIGH INCRE I1-I6
ID Identification integer
Index of node to which restraint applies from NLOW to NHIGH with index increment INCR
I1-I6 Flag to specify restrictions as a function of the direction (I1 = X-dir, I2 = Y-dir, I3 = Z-dir)
0 Free to move
1 Locked
X
Applied pressure is distributed onto the nodal points which belong to the face EFACE of the element considered. The faces are numbered as shown:
LOAD
LTEXT
P 6 ELOW DIR EFACE EHIGH INCRE
PRESSURE
LTEXT Arbitrary label
Pressure applied on the element face EFACE from element number ELOW to EHIGH with increment INCRE. The direction of the pressure load is in the global coordinate system DIR = 1 to 3 with 1 = X, 2 = Y and 3 = Z.
X
X
4. Validation of the code
Several calculations were performed both with this code and Stripe, an in-house finite element program developed at FFA. Good agreements were found between the calculations of the two programs. Figure 7 displays the percentage of change between the displacement field calculated by Stripe and the one calculated by the FE code, on each node of a wing. As shown in the legend, three different types of wings were considered (Airbus, Boeing 777, Rafale fighter as described in AeroWing program manual) under a pressure distribution due to a low velocity steady flow. In this example, a maximum relative difference in the results of approximately 3% was calculated. The accuracy of the results is highly dependent on the aspect ratio of the elements which composed the wing. With an appropriate meshing of the wing, with almost square shaped element (especially in the plane of the wing), displacements, stresses and strains can be calculated with a good accuracy.
Figure 7:
Differences in the displacement field (U, V, W) calculated by this FE code and Stripe, at each node.
5. Aeroelasticity results
After verifying that the module, i.e. the aerodynamic module and the finite element code module, gives coherent results, the whole program, combining aerodynamic and structure is verified in this section.
To achieve this goal a reference needed to be found to compare our results with. Being aware of the lack of data in the literature about the deflection of a finite wing induced by a steady airflow, the easiest way to proceed was to use a commercial finite element program including an aeroelastic module. MSC Nastran version 70.7 was used for this purpose.
Two types of wing design were investigated: a swept wing (15°) and a straight wing.
In both cases, a Mach number of 0.45 was chosen. For the swept wing, the angle of attack was set to 6° and 12° for the straight wing.
Figure 8: Deflections along several lines parallels to the leading edge for a straight wing
Figure 9: Deflections along several lines parallels to the leading edge for a swept wing
Figures 8 and 9 show the deflection along the axis of the wing parallel to the leading edge for 5 different depths equally distribute along the wing chord. Both Nastran and AeroWing (after one iteration) calculation deflection are plotted against the distance to the wing attach.
The straight wing appeared to give the reasonable results in terms of accuracy when compared with Nastran results. Indeed, a maximum relative error of 9% was found at the largest deflection. This error was increased by a factor two in the case of the swept wing. Improvements of the deflection calculation accuracy can be done by selecting an iterative resolution in AeroWing.
6. Further works
As first approximation, the wing section does not have any thickness distribution. Indeed, the thickness was constant with amplitude equal to the maximum amplitude of the airfoil profile. Note the camber line was respected in the model.
To model an airfoil, such as NACA profile, more elements are needed, which means that more CPU time and memory are necessary. In this case special algorithms to reduce the bandwidth of the stiffness matrix need to be implemented in the finite element code. Besides, wedge shape are expected to take place to improve the wing profile, therefore the finite element code needs to handle elements of 15 nodes.
The program can be rewritten for other codes that Matlab, such as SciLab.
In the case of composite wing, the filler made of advanced aluminium needs to be tailed in order to stick to the reality with reinforcements instead of massive material.
The composite section does handle only the unidirectional fibre, and other types of reinforcement forms can be needed for various purposes, and then need to be implemented.
Finally, the program handles only steady flow at low velocity. The reason came from the fact that the present finite element code did not handle dynamic load. Therefore the finite element program and the aerodynamic program need to be modified to take into account the dynamic case for the finite element code and the unsteady case for the aerodynamics.
Reference
[1] D. Tourde, " Design of a smart wing," FFA TN 1999-05, FFA, Sweden, 1999
[2] J. Katz and A. Plotkin, Low-Speed Aerodynamics: From Wing Theory to Panel Methods, McGraw-Hill, 1991
[3] T. Theodorsen, "General theory of aerodynamic instability and the mechanism of flutter," Report 496, National Advisory Committee for Aeronautics, 1935
[4] R. T. Jones and D. Cohen, "High Speed Wing Theory," Princeton Aeronautical Paperback, No. 6, Princeton University Press, Princeton, N. J., 1960
[5] L. Prandtl, "Applications of Modern Hydrodynamics to Aeronautics," NACA Report No. 116, 1921
[6] I. H. Abbott, A. E. Von Doenhoff and L. S. Stivers, "Summary of Airfoil Data," NACA Report 824, p. 142, 1945
[7] E. Hinton and D. R. J. Owen, Finite Element Programming, Academic Press, Department of civil engineering, University College of Swansea, UK, 1977
[8] I.M. Smith and D. V. Griffiths, Programming the Finite Element method, John Wiley & Sons, 2nd edition, University of Manchester UK, 1988
[9] Niels Saabye Ottosen and Hans Petersson, Introduction to the Finite Element Method, Prentice Hall, University of Lund, Sweden, 1992
(13)