## Symbolic Math Toolbox |

This example shows how to derive an analytical model for aircraft wing loads using symbolic computation. The model will include the major force components on an aircraft wing, including the aerodynamic lift, load from the weight of the wing structure, and load from the weight of the fuel stored in the wing.

We start by modeling the aerodynamic lift, which is depicted below.

The lift profile can be estimated using an elliptic distribution, . We calculate the total lift by integrating across the length of the wing:

q_l:= x --> ka*sqrt(L^2 - x^2)

LiftTotal := int(q_l(x), x=0..L)

We would like to express ka in terms of meaningful parameters. We can do this by equating the above lift expression with the lift expressed in terms of the aircraft's load factor, which is the ratio of lift to total aircraft weight, . We equate the 2 lift expressions and solve for ka. Note that our calculation assumes that lift forces are concentrated on the two wings of the aircraft, which is why the left-hand side of the equation is divided by 2. We do not consider lift on the fuselage or other surfaces.

eqn := n*Wto/2 = LiftTotal

ka := (solve(eqn,ka) assuming L>0)[1]

We plug the resulting ka term into our expression for .

q_l(x)

In addition to lift, there are also downward loads resulting from the weight of the wing. One downward load is due to the weight of the wing structure itself, while the other is due to the weight of the fuel stored in the wing. We derive analytical models for these loads separately.

__Wing structure__

The load due to the weight of the wing structure is depicted below.

We assume that this load is proportional to chord length, which is maximum at the base of the wing (Co), and tapers off linearly as you approach the tip of the wing (Ct). The load profile can therefore be expressed by the following equation: . We define our expression for and integrate it across the length of the wing to calculate the total load due to the weight of the wing structure.

q_w := x --> kw*(x*(Ct - Co)/L + Co)

structLoadTotal := int(q_w(x), x=0..L)

As we did when we calculated aerodynamic lift, we equate the resulting load expression with the structural load expressed in terms of load factor and weight of the wing structure (Wws), and solve for kw.

eqn := n*Wws/2 = structLoadTotal

kw := -(solve(eqn,kw) assuming L>0 and Co + Ct > 0)[1]

We plug the resulting kw term into our expression for.

q_w(x)

__Fuel stored in wing__

Placing fuel in aircraft wings is a common design practice because the weight of the fuel counteracts the lift force, reducing the net load on the wing. The fuel load is depicted below.

The magnitude of this load is calculated in the same way as the structural load, resulting in a solution with the same form. However, because the fuel storage does not extend to the wing tip, the load profile ends midway through the wing. Therefore, we define the fuel load as a piecewise function.

q_f:= x --> piecewise([x <= Lf, -(Wf*n)/((Ctf + Cof)*Lf)*(x*(Ctf-Cof)/Lf + Cof)],[x > Lf, 0])

We calculate total load by adding the aerodynamic lift, structural load, and fuel load.

q_t:= x --> q_l(x) + q_w(x) + q_f(x)

simplify(q_t(x))

This analytical model gives a clear view of how aircraft weight and geometry parameters affect total load.

We would like to evaluate total load for our specific aircraft, which has the following parameters.

Wto := 4800:

Wws := 630:

Wf := 675:

L := 7:

Lf := 2.4:

Co := 1.8:

Ct := 1.4:

Cof := 1.1:

Ctf := 0.85:

Wws := 630:

Wf := 675:

L := 7:

Lf := 2.4:

Co := 1.8:

Ct := 1.4:

Cof := 1.1:

Ctf := 0.85:

We plot the individual load components and total load on the aircraft wing, assuming a load factor of 1.5.

P1 := plot::Function2d(q_l(x) | n=1.5, x=0..L,

LineColor=RGB::Red,Legend="Lift",LineStyle=Dashed):

P2 := plot::Function2d(q_w(x) | n=1.5, x=0..L,

LineColor=RGB::Blue,Legend="Wing Structure",LineStyle=Dashed):

P3 := plot::Function2d(q_f(x) | n=1.5, x=0..L,

LineColor=RGB::Green,Legend="Fuel",LineStyle=Dashed):

P4 := plot::Function2d(q_t(x) | n=1.5, x=0..L,

LineColor=RGB::Black,Legend="Total load",

LineStyle=Solid,LineWidth=1):

plot(P1,P2,P3,P4,

AxesTitles = ["Length", "Load (N/m)"],Width=160,Height=120)

LineColor=RGB::Red,Legend="Lift",LineStyle=Dashed):

P2 := plot::Function2d(q_w(x) | n=1.5, x=0..L,

LineColor=RGB::Blue,Legend="Wing Structure",LineStyle=Dashed):

P3 := plot::Function2d(q_f(x) | n=1.5, x=0..L,

LineColor=RGB::Green,Legend="Fuel",LineStyle=Dashed):

P4 := plot::Function2d(q_t(x) | n=1.5, x=0..L,

LineColor=RGB::Black,Legend="Total load",

LineStyle=Solid,LineWidth=1):

plot(P1,P2,P3,P4,

AxesTitles = ["Length", "Load (N/m)"],Width=160,Height=120)

While it is useful to visualize wing loads, what we are most interested in are the shear forces and bending moments resulting from these loads. We calculate shear force by integrating total load along the length of the wing, and then calculate bending moment by integrating shear force. We wrote custom MuPAD procedures called CalcShear and CalcMoment that calculate shear force and bending moment for a given load profile.

__Shear Force__

Shear force is calculated by integrating total load across the length of the wing.

read(NOTEBOOKPATH."CalcShear.mu");

shear_force := CalcShear(q_t(x))

shear_force := CalcShear(q_t(x))

plotfunc2d(shear_force | n=1.5, x=0..L, AxesTitles = ["Length", "Shear (N)"])

__Bending Moment__

Bending moment is calculated by integrating shear force across the length of the wing.

read(NOTEBOOKPATH."CalcMoment.mu");

bending_moment := CalcMoment(q_t(x))

bending_moment := CalcMoment(q_t(x))

plotfunc2d(bending_moment | n=1.5, x=0..L, AxesTitles = ["Length", "Moment (N*m)"])

The maximum bending moment is approximately 10 kN*m, and occurs at the base of the wing. We must ensure that our wing design can handle this bending moment.