  When modeling engineering systems, it can be difficult to identify the key parameters driving system behavior because they are often deep within the model. Analytical models can help because they describe systems using mathematical equations, showing exactly how different parameters affect system behavior. Analytical models can review the loads and bending moments on the wing of a small passenger aircraft, determining whether the wing design meets the strength requirements. Models are first derived in the notebook interface in Symbolic Math Toolbox, and then the data management and analysis tools in MATLAB are used to simulate the models for different scenarios to verify that anticipated bending moments are within design limits.

Leading aerospace companies use analytical models to describe the systems and processes that they work with. While this example is specific to aircraft wing analysis, analytical models are useful throughout the aerospace industry; for example, they can be used to describe the trajectory of mechanical elements in a landing gear, or to model pumps, compressors, and other mechanical and electrical systems in an aircraft.

To being, one must take a look at the evaluation of three primary loads that act on the aircraft wing: aerodynamic lift, load due to wing structure weight, and load due to the weight of the fuel contained in the wing. These loads act perpendicular to the wing surface, and their magnitude varies along the length of the wing (Figures 1a, 1b, and 1c).

The analytical model of wing loads can be determined in the Symbolic Math Toolbox notebook interface, which offers an environment for managing and documenting symbolic calculations. The notebook interface provides direct support for the MuPAD language, which is optimized for handling and operating on symbolic math expressions.

Equations for each load component are derived separately and then the individual components are added to obtain total load. Assuming an elliptical distribution for lift across the length of the wing, resulting in the following expression for lift profile:  q1(x) = ka √L2 – x2; where L = length of wing; x = position along wing; and ka = lift profile coefficient. The total lift can be determined by integrating across the length of the wing: Lift = 0∫L ka √L2 – x2 dx.

The notebook interface will be able to define ql(x) and calculate its integral (Figure 2). Incorporated in the notebook interface are math equations, descriptive text, and images into the calculations to clearly document the customers work. Completed notebooks may be published in PDF or HTML format. Figure 2: The notebook that is showing the derivation of lift on an aircraft wing.Through integration, it is found that Lift = ∏ L2 ka / 4. Ka is solved by equating the lift expression that we just calculated with the lift expressed in terms of the aircraft’s load factor. In aircraft design, load factor is the ratio of lift to total aircraft weight: n = Lift / Wto. Load factor equals one during straight and level flight, and is greater than one when an aircraft is climbing or during other maneuvers where lift exceeds aircraft weight.

By equating two lift expressions, Wto n / 2 = ∏ L2 ka / 4, the unknown ka term is solved. The analysis 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 two and the lift on the fuselage or other surfaces are not considered.

Plugging ka into the original ql(x) expression, an expression for lift may be obtained: ql(x) = 2 Wto n √L2 – x2 / L2 ∏.
An analytical model like this helps the customer understand how various parameters affect lift. For example, when lift is directly proportional to load factor (n) and that for a load factor of one, the maximum lift 2 Wto / ∏ L, occurs at the wing root (x=0).

Weight of wing structure
Assuming that the load caused by the weight of the wing structure is proportional to chord length (the width of the wing), which is highest at the wing base (Co) and tapers off toward the wing tip (Ct). Therefore, the load profile can be expressed as qw (x) = kw ((Ct – Co / L)x + Co). To define qw(x) and integrate it across the length of the wing the total load from the wing structure must be calculated: x −> kw (Co – x (Co – Ct) / L)) and L kw (CO + Ct) / 2.

Then the structural load is determined by using the equation with the structural load expressed in terms of load factor and weight of the wing structure (Wws) is: Wws n / 2 = kw L (Co + Ct) / 2, and solve for kw.

Plugging kw into the original qw(x) expression, an analytical expression is obtained for load due to weight of the wing structure. qw(x) = – Wws n (Co – x(Co – Ct) / L) / L (Co +Ct).

Weight of fuel stored in wing
Defining the load from the weight of the fuel stored in the wing as a piece-wise function, where load is zero when x > Lf, it can be assumed that this load is proportional to the width of the fuel tank. This fuel tank is at its maximum at the base of the wing (Cof), and tapers off as it approaches the tip of the fuel storage tank (Ctf). It may be derived qf(x) in the same way that qw(x) is derived, resulting in an equation of the same form: qf(x) = 0 if Lf < x and {Wf n (Cof – x (Cof – Ctf) / Lf) / Lf (Cof + Ctf) if x ≤ Lf.

Calculating the total load is accomplished by adding the three individual load components. This analytical model gives a clear view of how aircraft weight and geometry parameters affect total load. qt(x) = {n (2 Co Wto √L2 – x2 + 2 Ct Wto √L2 – x2 – ∏ Co L Wws + ∏ Co Wws x – ∏ Ct Wws x) / L2 ∏ (Co + Ct) if Lf < x and 2 Wto n √L2 – x2/ L2 ∏ – Wws n (Ct x – Co x + Co L) / L2 (Co + Ct) – Wf n (Ctf x – Cof x + Cof Lf) / Lf2 (Cof + Ctf) if x ≤ Lf.

After completing the analytical model for wing loads, the customer can use the model to evaluate aircrafts with various wing dimensions and weights. The small passenger aircraft have the following parameters:

• Wto = 4800kg
• total aircraft weight
• Wws = 630kg (weight of wing structure)
• Wf = 675kg (weight of fuel stored in wing)
• L = 7m (length of wing)
• Lf = 4.8m (length of fuel tank within wing)
• Co = 1.8m (chord length at wing root)
• Ct = 1.4m (chord length at wing tip)
• Cof = 1.1m (width of fuel tank at wing root)
• Ctf = 0.85m (width of fuel tank at Lf)

To evaluate load during the climb phase, it is assumed that there is a load factor of 1.5, then plot the individual load components and total load.

It is obvious that lift is the largest contributor to total load and that the maximum load of 545 N*m occurs at the end of the fuel tank. Fuel load also contributes significantly to the total load, while the weight of the wing is the smallest contributor.

While it is useful to visualize wing loads, what becomes a major concern is the shear force and bending moments resulting from these loads. The next step is to determine whether worst-case bending moments experienced by the wing are within design limits.

A bending moment model
The expressions that are derived can be used for load on the wing to calculate bending moment. The first step is to start by integrating total load to determine shear force: V(x) = – ∫ qr(x) dx. The bending moment can then be calculated by integrating shear force: M(x) = ∫ V(x) dx.

A custom function is written in the MuPAD language, CalcMoment.mu, that accepts load profile and returns the bending moment along the length of the wing. Symbolic Math Toolbox includes an editor, debugger, and other programming utilities for authoring custom symbolic functions in the MuPAD language.

Using this function with the aircraft parameters that were previously defined to obtain an expression for bending moment as a function of length along wing (x) and load factor (n).

As with wing loads, the user plots bending moment assuming a load factor of 1.5.

As expected, the bending moment is highest at the wing root with a value of 8.5 kN*m. The wing is designed to handle bending moments up to 40 kN*m at the wing root, but since regulations require a safety factor of 1.5, bending moments exceeding 26.7 kN*m are unacceptable. The user will simulate bending moments for various operating conditions, including conditions where the load factor is greater than 1.5, to ensure that they are not in danger of exceeding the 26.7 kN*m limit.

Simulating bending moment
The bending moment equation is then saved in our notebook as moment. Using the getVar command in Symbolic Math Toolbox, the user will import this variable into the MATLAB workspace as a sym object, which can be operated on using MATLAB symbolic functions included in Symbolic Math Toolbox – bendingMoment = getVar(nb,’moment’).

Too numerically evaluate the bending moments for various conditions, the user must convert the sym object to its equivalent numeric function using the matlabFunction command – h_MOMENT = matlabFunction(bendingMoment).

h_MOMENT is a MATLAB function that accepts load factor (n) and length along wing (x) as inputs. By evaluating the bending moments at the wing root (x=0), load factor becomes the only variable that affects bending moments. As mentioned earlier, load factor is equal to Lift / Wto. Using the standard lift equation, and assuming the aircraft is not banking, load factor is then expressed as: n = p A CL V2 / 2 Wto.

Where p = air density – 1.2 kg/m^3; A = planform area (approximately equal to total surface area of wings) – 23 m^2; CL = lift coefficient (varies with aircraft angle of attack, which ranges from 3 degrees to 12 degrees) – 0.75 to 1.5; V = net aircraft velocity (accounts for aircraft speed and external wind conditions) – 40m/s to 88m/s.

These variables can be defined in MATLAB and evaluate the h_MOMENT for the range of lift coefficients and aircraft velocities listed above. We store the results in a dataset array, available in Statistics Toolbox.

Dataset arrays provide a convenient way to manage data in MATLAB, enabling the user to filter the data to view and analyze the subsets that they are most interested in. Since the user want to determine whether bending moments ever exceed the 26.7 kN*m threshold, they only need the bending moment data where the load factor is greater than 1.5. We filter the dataset and plot the data for these conditions.

moment_filt = moment (moment.loadFactor>1.5,:) x = moment_filt.netVel; y = moment_filt.CL; z = moment_filt.maxMoment; [X,Y] = meshgrid(x,y); Z = griddata(x,y,z,X,Y); surf(X,Y,Z).

The plot shows that the peak bending moment, 19.3 kN*m, occurs when net aircraft velocity and lift coefficient are at their maximum values of 88m/s and 1.5, respectively. This result confirms that bending moments will be safely below the 26.7 kN*m limit, even for worst-case conditions.

The Value
The analytical models gave the user a clear view into how different aircraft parameters and operating conditions affect loads and bending moments on the aircraft wing. They were able to verify that the wing can withstand worst-case loading conditions that it could encounter during the climb phase of flight.

The models discussed in this article were used only for high-level proof-of-concept analysis, but analytical models could also be used for more detailed system modeling tasks – for example, to model the airflow near the leading edge or tip of the wing.