![m inimize m (R ) | R ↦→ {r(x),R(x)} | x ∈ [0,L]
R
subje∀cxtto rmin ≤ r(x ) ≤ rmin + Tmin
Rmax - Tmin ≤ R(x) ≤ Rmax
r(x)- R (x ) ≤ - Tmin
σ (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 sequential quadratic programming, several numerical differentiation paradigms are investigated in relation to the problem’s objective and constraints. It is shown that an optimal spar mass, over 60% lower than the nominal, is consistently achieved while satisfying all constraints.
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 8.
| 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 while not allowing its structural stress σ at any spanwise location to exceed its material’s ultimate strength σmax during its maximal operational loading conditions. 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 follow:
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, it is assumed that the force distribution in the spanwise direction has an approximately linear distribution, with maximum load at the root and zero load at the tip. Hence the force at each node can be formulated as follows
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
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}
T
U = [υk,...,υ2n] | {υ2k-1,υ2k} ← {Rmax - Tmin,Rmax}](TemplateICATT9x.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,2.n|
A = ⌈ .. .. .. ⌉ | {αi,2i-1,αi,2i} ← {1,- 1}
αn,k ⋅⋅⋅ αn,2n
b = [β ,...,β ]T | β ← - T
i n i min](TemplateICATT10x.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)](TemplateICATT14x.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 2. 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. For the
sake of testing robustness, the initial guess of the 2n × 1 design
vector R0 will characterize the boundary constraints of the spar’s
radii, such that {ri,Ri}←{rmin,Rmax}∀i ∈, yielding an
initial total spar mass of m0 = 90.4779 kg, far beyond the
desired optimal value m* < 5 kg.
One can easily instantiate an initial guess R0 of the design
vector, with n = 20 nodes for instance, as shown in Listing 3.
Once the design vector has been characterized, the user can easily begin the optimization process by executing [Ropt,fval,exitflag,output,lambda,grad,hessian] = Raptor.Optimize(R, step_type);, where step_type can be 'forward', 'central', or 'complex'. With n = 20 nodes, the forward, central, and complex step optimization processes yield the outputs in Listings 4, 5, and 6, respectively. It can quite easily be seen that all three of the finite difference methods yield approximately equal values of the optimal spar mass m* = 4.9069, thereby verifying beyond a doubt that m* is certainly a local minimum of the optimization problem.
Visually, the spar geometry, vertical displacement, and stress can be see for all the finite difference methods in Figures 2, 4, and 6. Additionally, the iterative convergence of the objective function can be see for each method in Figure 8. It can be seen that each finite difference method yields virtually identical results. Each method yielded a final spar mass that is 63% better than the nominal mass of 13.26 kg. However, it was quite apparent that the optimization process using the complex step method was far more expeditious than the other methods.
4.
A virtually identical value of the spar’s optimal mass m* = 4.9069, over 60% better than the nominal mass, has been achieved consistently across three separate finite difference methods. Validity can thereby be concluded for the optimality of the resultant masses and design vectors, as well as the method by which the complex Jacobian matrices were computed in the complex step implementation. Through implementation, as expected, it was seen that the complex step method yielded solution convergence in a fraction of the time it took the forward and central step methods to do so. Lastly, from the geometrical representations of the resulting design vectors, an intuition can be developed, from which a wing’s spar will tend to have greater radii towards its root. It can be seen that the optimization process has indeed solidified this intuition, by yielding geometries which meet the upper constraints towards the root and lower constraints towards the tip.
5.
5.1.
