![m inimize m(R ) | R ↦→ {r(x),R(x)} | x ∈ [0,L]
R
subj∀ecxtto rmin ≤ r(x) ≤ rmin + Tmin
R - T ≤ R (x) ≤ R
max min max
r(x)- R (x) ≤∘- Tmin-----
E(σ(x,ξ))+ 6 V ar(σ (x,ξ))- σmax ≤ 0](TemplateICATT0x.png)
ABSTRACT
This paper examines the design of a wing for a new long-endurance aircraft that will be used as a telecommunications platform. In particular, the problem of optimizing the shape of the wing’s spar, acting as the aircraft’s main structural support, in order to minimize its mass, is analyzed and solved. In doing so, a robust object-orientated programming approach is taken to transcribe the problem and flexibly find optimal solutions from user defined design guesses. Using a complex step finite difference method and stochastic collocation, the geometry of the wing’s spar is optimized with respect to its mass, under uncertain loading conditions. It is shown that a design mass of %70 percent lower than the nominal is feasible.
Index Terms— Wing Structure, Shape Optimization
1.
1.1.
Preceding the construction of the optimization problem, it is necessary to define the constant metrics that characterize the nature of the problem. For the case of an annular spar, the problem parametrization requires several design metrics to be taken into account, namely: wing semi-span L, material density ρ, Young’s modulus E, ultimate tensile/compressive strength σmax, aircraft operational mass M, minimum spar thickness Tmin, minimum inner wall radius rmin, and maximum outer wall radius Rmax.
Additionally, the operational conditions that the aircraft will be subjected to need to be taken into account, namely: gravitational acceleration g, and maximum operational G-force G. In this particular investigation, attention will be directed towards the case of a carbon fiber spar, whose specifications are enumerated in Table 1. One can instantiate a UAV of this nature through an object-oriented approach, as shown in Listing 1. The full source code of the UAV class is shown in Listing 5.
| Parameter | Symbol | Value |
| Semi-Span | L | 7.5 m |
| Density | ρ | 1600 kg∕m3 |
| Young’s Modulus | E | 70 GPa |
| Ultimate Strength | σmax | 600 MPa |
| Aircraft Mass | M | 500 kg |
| Minimum Thickness | Tmin | 2.5 mm |
| Minimum Inner Radius | rmin | 1 cm |
| Maximum Outer Radius | Rmax | 5 cm |
| Gravity | g | 9.81 m∕s2 |
| G-Force | G | 2.5 |
The objective of this optimization problem is to minimize the total mass of the spar through structural alteration, while not allowing its structural stress σ at any spanwise location to exceed its material’s ultimate strength σmax during its maximal operational loading conditions. In this particular optimization process, the uncertainty of the wing’s load distributions is taken into account. That is, the stress constraint of the spar becomes the wing’s average stress plus six standard deviations. This statistical constraint characterization is discussed further in Section 1.4.1. Additionally, there are also manufacturing constraints that dictate the spar’s inner radius r and outer radius R to satisfy r ≥ rmin and R ≤ Rmax, while satisfying a thickness constraints R - r ≥ Tmin. In a more formal manner, this problem’s formulation is shown as:
1.3.
Because finite element methods will be used to solve this
optimization problem, it is required to discretize the spar’s
geometry into a series of spanwise finite elements. Herein, the
wing’s spar will be represented by a series of nodes, indexed by i,
where i = 1 and i = n signify the locations of the wing’s root and
tip respectively. As such, it can be deduced that there are n nodes
and n- 1 finite elements along the span of the spar. Hence, the set
of node indexes 
and finite element indexes 
are defined as
follows
![N = [i,i+ 1,...,n]
E = [j,j + 1,...,n- 1]](TemplateICATT1x.png)
1.3.1.
It is appropriate to begin the optimization process with the formulation of the design vector R, composed of inner radii r and outer radii R. Abiding to the form of a column vector, it is necessary to apply a change of variables and define the design vector as follows,
![R = [ξk,ξk+1,...,ξ2n]T
{ξ2k-1,ξ2k} ← {ri,Ri}∀i ∈ I](TemplateICATT2x.png)
1.3.2.
Having deduced the number of nodes n from the user supplied design vector, the spanwise location along the spar can be discretized in a similar manner as,
![x = [xi,xi+1,...,xn] | x1 = 0,xn = L](TemplateICATT3x.png)
1.3.3.
At each node i along the wing’s span, the spar’s cross section is represented as a circular annulus of inner and outer radii, ri and Ri, respectively, whose area is defined by Ai = π(Ri2 - ri2). Assuming this cross sectional area characterization is consistent across all nodes, the inner and outer surfaces of each finite element j are linearly interpolated between their corresponding nodes i = j and i = j + 1. Hence the nodes’ cross sectional areas are integrated along the interpolations, yielding the volume of each finite element as follows
1.4.
As formulated in Section 1.2, the structural stress at any spanwise
location along the spar must not exceed its material’s ultimate
strength, σi - σmax ≤ 0 ∀i ∈
. Modeling the wing’s spar as a
cantilever beam, Euler-Bernoulli beam theory can be exploited to
compute the spar’s nodal stress.
1.4.1.
In this study, a linear load distribution is no longer assumed. Rather, the force distribution experienced by the aircraft’s spar is formulated as the nominal linear force distribution plus a probabilistic perturbation given by
. Hence the force along the wingspan
can be formulated as
In order to calculate the displacement of the spar at any location along its span, the second moment of inertia at each node must first be calculated. For a circular annulus describing the cross section of any node along the spar’s span, the second moment of area is formulated as such
1.4.3.
Assuming a transversely applied load distribution and cross sectional symmetry, a cubic Hermite finite-element basis is used to solve the Euler-Bernoulli beam equations. The vertical and angular displacement of the spar at each node along the beam are calculated, using the aforementioned formulations of the load distribution and second moment of area. The beam’s vertical displacement is formulated as follows δ = [δi,δi+1,…,δn]T. Using the spar’s displacements, its nodal tensile stresses are computed and formulated as such σ = [σi,σi+1,…,σn]T. Hence the constraint function must satisfy the following
The unadjusted abscissas and weights corresponding to the m
collocation points are computed through the Gauss-Hermite
subroutine given in Listing 2.
For further background on the derrivation of this method, the reader is advised to consult [?]. 1.5.
As required by the optimization problem’s specifications, the inner and outer radii of each annular cross section cannot be less than Tmin apart, the inner radius cannot be smaller than rmin, and the outer radius cannot be larger than Rmax.
1.5.1.
In formulating the range to which the inner radii ri and outer radii Ri are bounded, it is necessary to formally define lower and upper bound vectors L and U, respectively, that satisfy L ≤ R ≤ U, corresponding to the form of the transformed design vector R described in Section 1.3.1. The lower and upper bounds for the radii, corresponding to the constraints described in Section 1.2, are formulated as such
![L = [λk,...,λ2n]T | {λ2k-1,λ2k} ← {rmin,rmin +Tmin}
U = [υ ,...,υ ]T | {υ ,υ } ← {R - T ,R }
k 2n 2k- 1 2k max min max](TemplateICATT14x.png)
1.5.2.
In order to enforce the constraint that the annular spar’s wall
thickness may not be smaller than Tmin and that the outer radii
must always be greater than the inner radii Ri > ri ∀i ∈, linear
inequality constraints must be instantiated. In doing so, it is
required to formulate a n × 2n matrix A and n × 1 column
vector b that satisfy AR ≤ b. The matrices are defined as
follows
![⌊ ⌋
αi,k ⋅⋅⋅ αi,2n
A = |⌈ ... ... ... |⌉ | {αi,2i-1,αi,2i} ← {1,- 1}
α ⋅⋅⋅ α
n,k n,2n T
b = [βi,...,βn] | βi ← - Tmin](TemplateICATT15x.png)
2.
2.1.
Sequential quadratic programming (SQP) has proven itself as one of the most successful methods for solving nonlinearly constrained optimization problems. Essentially, SQP aims to model the nonlinear programming problem (NLP) at a given approximate solution Rh by a quadratic programming subproblem, then use the subproblem to construct a better approximation Rh+1. This process is iterated until convergence to a solution R* is achieved. SQP is not a feasible point method; neither the point nor any of the subsequent iterates are required to be feasible. However, the SQP method does satisfy the bounds L and U at every iteration, and is ideal for small to medium scale problem such as in the case of this paper’s optimization problem, where optimality returns are generally diminished for n > 40. This method of nonlinear constrained optimization is handled natively by Matlab’s fmincon function [?], and described in more detail in [?].
2.2.
In order to achieve an optimal design vector R* that satisfies
m(R*) ≤ m(R)∀R
{ri,Ri}∀i ∈, both the necessary and
sufficient conditions for optimality need to be realized. It is
necessary for an optima at a stationary point to have a
first derivative of zero ∇ = 0. It is sufficient to call the
stationary point a local minimum R* if ∇2R* > 0. For further
information on conditions of optimality, the reader should consult
[?].
2.2.1.
Among the simplest methods of finite difference derivative approximations is the forward step method
(hi) represents the truncation error, and hi represents the
step size. Herein, hi represents the positive distance from
|ξi| to the next larger floating-point number of the same
precision as ξi. It should be noted that f may represent either the
problem’s objective function m(R) or constraint function
c(R).
2.2.2.
Similarly to the forward difference method, the approximate derivative of either the objective or constraint functions may be approximated with a perturbation to the design vector R. Under the central difference paradigm, the finite difference derivative approximation becomes
(hi2). Therefor the central difference method is more
accurate than the forward difference method due to its smaller
truncation error.
2.2.3.
Among the best of the finite difference methods is complex-step. It has shown to be very accurate, robust, and easy to implement, while maintaining a reasonable computational cost. The relationship between the real and imaginary parts of the function at hand is exploited to yield a derivative approximation of the following form
![√---
ℑ[f(R-+-hiei---1)] 2
∇ ξif(R) = hi + O(hi)](TemplateICATT19x.png)
The complex step method converges quadratically with decreasing step size. This method is particularly advantageous because it is insensitive to small step sizes and eventually achieves the accuracy of its function evaluations. The reader can learn more about its derivation in [?]. Herein, focus will be directed towards the use of the complex step method to approximate the gradients of both the objective function m(R) and constraint function C(R). The complex step method can be formulated to approximate the functions’ Jacobian matrices as shown in Listing 3. As an example, to approximate the gradient of the objective function, one should execute obj.Complex_Jacobian(@obj.Stress_Constraints, R).
Having instantiated an aircraft with the specifications
enumerated in Table 1, one can supply an initial guess to
the optimizer. The initial guess of the 2n × 1 design vector
R0 will characterize the nominal design of the spar’s radii,
such that {ri,Ri}←{0.0415,0.05}∀i ∈, yielding an
initial total spar mass of m0 = 29.32 kg. With the use of the
complex step finite difference method and a varying amount of
collocation points, the optimizer settled upon a final mass of
m = 8.59221 kg, which is a %70.7026 decrease from the
nominal mass, as shown through the convergence of the
optimizer in Figure 4. The resulting spar geometries are seen in
Figure 6. As one can see from Figure 8, the average stress
and standard deviations are markedly different from the
corresponding nominal plot in Figure 2. It is also seen throughout
the figures that varying the number of collocation points did not
affect the final optimal mass of the spar. It is seen that the
new constraints enforced by the uncertainty model have
made their effect known, as it is seen in Figure 8 that the six
standard deviations above the mean stress have settled upon the
constraint.
Complex Step Convergence with 25 Nodes and 3 Collocation Points
Complex Step Convergence with 25 Nodes and 5 Collocation Points
Optimized Design with 25 Nodes and 3 Collocation Points
Optimized Design with 25 Nodes and 5 Collocation Points
Optimized Stress Statistics with 25 Nodes and 3 Collocation Points.
Optimized Stress Statistics with 25 Nodes and 5 Collocation Points.
