This MATLAB program solves multiple complex loading conditions on aircraft wing by assuming the wing is made out of two Euler-Bernoulli beams with square cross section. This description includes the scenario given, program's objective, proposed solution, results from MATLAB and lastly a table of content explaining each file in this repository.
[Figure 1] show the orthographic drawing of a flight. As a structure engineer, you are requested to estimate the size of supporting beam for the wing. Based on the design, the wing is supported by two square beams as show in Figure 2 and the material for the beam is Carbon – Graphite with the modulus of elasticity 25 GPA.
Given the dry weight for the Turbine is 1950kg, and the uniform distribution weight for wing is 135 kg/meter. The distribution of fuel weight a along the wing is shown in the [Figure 3].
Let assume the lift force acting on the wing as in [Figure 4], and then the distribution lift force can be calculated in the equation below.
Where,
- a) To determine the cross section of the beam if the maximum deflection of the wing must not larger than 0.5 meter.
- b) To calculate the location and at what condition happen the maximum transverse stress.
- c) To calculate the location and at what condition happen the maximum bending stress.
- d) To plot the 2D graph maximum shear force verse speed for the angle of attack 0 deg., 10 deg. and 15 deg. in same graph.
- e) To plot the 2D graph maximum bending moment verse speed for the angle of attack 0 deg., 10 deg. and 15 deg. in same graph.
- f) To plot the 2D graph deflection verse speed for the angle of attack 0 deg., 10 deg. and 15 deg. in same graph.
- g) To plot the 3D graph for shear force with the step size for speed 50km/h and step size for angle attack 1 deg.
- h) To plot the 3D graph for bending moment with the step size for speed 50km/h and step size for angle attack 1 deg.
- i) To plot the 3D graph for deflection with the step size for speed 50km/h and step size for angle attack 1 deg.
Before constructing any codes in MATLAB, a correct Macaulay equation is very important for any further calculations. Therefore, we need to first draw a free body diagram of the scenario as shown in [Figure 5].
Given :
- Weight of fuel = 15.2kg/m at x=0 , 7.78kg/m at x=15.3
- Weight of turbine = 1950kg at x=3
- Weight of wing = 135kg/m throughout
- Lift force = $ -2/17 \rho xv^2 C_L+2.55\rho v^2 C_L $ N/m throughout
Where
-
$\rho$ = density of air 1.225 kg/m3 -
$v$ = speed (0-750 km/h) -
$C_L$ = lift coefficient depends on angle of attack ($C_L=0.09\rho+0.4$ ) -
$\alpha$ = angle of attack (0-15 deg.) -
$x$ = distance from body
By applying
We get
By applying
We get
In both the equations above, note that there are two unknowns, the speed of the aircraft(v) and the lift coefficient(C_L). Therefore, when we code in MATLAB, looping is required to evaluate the different reaction forces(R) and moments(M) at the fixed support at different speed and angle of attack and stores the values in the form of 2D array.
After that, Macaulay equation can be constructed according to the free body diagram as shown in [Figure 6] below :
We get
Integrate the equation above
Above is the slope equation, to fulfil the requirements of this project, we need the deflection equation where the slope equation will be integrated again.
Above is the deflection that we desire. However, before we code the equation into MATLAB, we have to solve for the unknows C and D.
From slope equation, when x=0, slope at fixed support is 0,
From deflection equation, when x=0, deflection at fixed support is 0,
Therefore, our final deflection equation is as follow :
To ease the coding process, the equation above will be split into few small and simple parts by using superposition theory.
- Fuel weight component =
$-(1/24)(15.2)(9.81)x^4+(7.42)(9.81)/(120)(15.3)\ x^5$ - Turbine component =
$-\left(\frac{1}{6}\right)\left(1950\right)\left(9.81\right)(x-3)^3$ - Wing weight component =
$-\left(\frac{1}{24}\right)\left(135\right)\left(9.81\right)x^4$ - Lift force component =
$+\left(\frac{1}{24}\right)\left(2.55\right)\left(1.225\right)v^2C_L x^4-\left(\frac{2}{17}\right)\left(1.225\right)v^2C_L\frac{x^5}{120}$ - Moment at support =
$+\frac{M}{2}x^2$ - Reaction at support =
$+\frac{R}{6}x^3$
The varying variable in the deflection equation is x, a looping is then made to evaluate the deflection at x=0 to x=15.3. this looping should be within the looping for speed and angle of attack as those will affect the value of R and M for every loop.
To determine the cross-section of the beam with a maximum deflection of 0.5m, we need to calculate the value of area moment of inertia, I.
To get the I for a single beam, we divide the I by 2,
To calculate the cross-sectional area, we apply the formula of area moment of inertia
The cross-sectional area,
To determine the location of position and location of maximum transverse stress and maximum bending stress, the equations of both stresses must be studied.
Transverse stress =
Where
-
$V$ = Shear force -
$Q$ = First moment of area -
$I$ = Second moment of area -
$T$ = Thickness of the cross-section -
$M$ = Bending moment -
$c$ = Radius of interest
From the two equations above, we can conclude that the first moment of area, second moment of area, thickness of the cross-section and the radius of interest is fixed along the beam. The only differentiating factor in both cases is the shear force and the bending moment. Thus, to determine the maximum transverse stress and maximum bending stress, we need to consider only the maximum shear force and maximum bending moment.
Since the values of force and bending moment is stored in 3D arrays, we can simply quote the location of the maximum values by MATLAB command of “ind2sub”. It will return the position of the maximum value in 3 single variables for row, column and page.
To plot graphs in MATLAB, commands like “plot(x,y)” and “surfc(x,y,z)” are used. The size of the array in x,y and z axis must be equal to each other. Commands like “meshgrid”, “xlabel”, “ylabel”, “zlabel”, “title” are also used to provide a better visual representation of the graphs.
The following is the result from MATLAB program attached.
-
For question a)
- The cross section of the beam is 0.609 m^2.
-
For question b)
- The maximum transverse stress occurs at point x=0.00,with the conditions of speed 750 km per h and angle of attack 15 degree.
-
For question c)
- The maximum bending stress occurs at point x=0.00,with the conditions of speed 750 km per h and angle of attack 15 degree.
-
For question d), e), f), g), h), i),
| File | Description |
|---|---|
| Macaulays_method.m | MATLAB program solving the scenario explained |
| Figure1.png | Aircraft orthographic drawing |
| Figure2.png | Beam arrangement of wing |
| Figure3.png | Distribution of fuel weight along the wing |
| Figure4.png | Projected area of the wing |
| Figure5.png | Free body diagram of load distribution on wing |
| Figure6.png | Macaulay's method on the free body diagram |