Optimizing Aircraft Wing Material Distribution using Ant Colony Optimization (ACO)
Overview
The solution is written fully in Python and it aimes to optimize the material distribution of a simplified aircraft wing cross-section to minimize its weight while satisfying structural constraints. The wing is modeled as a thin-walled rectangular section subjected to an aerodynamic lift load.
It is essential to select the thickness and fiber orientation of composite material elements (skin and spars) to ensure the structure can withstand stress, buckling, and displacement constraints, and maintain a minimum natural frequency.
Given Data:
- Wing Geometry:
- Chord length: \( c = 2.0 \, \text{m} \)
- Thickness (height): \( h = 0.1 \, \text{m} \) (5% of chord)
- Span: \( b = 10.0 \, \text{m} \) (but we analyze a 2D cross-section)
- Two spars located at \( x = 0.5 \, \text{m} \) and \( x = 1.5 \, \text{m} \)
- Material Properties (default composite material):
- Longitudinal Young’s modulus: \( E_1 = 135 \, \text{GPa} \)
- Transverse Young’s modulus: \( E_2 = 10 \, \text{GPa} \)
- Shear modulus: \( G_{12} = 5 \, \text{GPa} \)
- Poisson’s ratio: \( \nu_{12} = 0.3 \)
- Density: \( \rho = 1600 \, \text{kg/m}^3 \)
- Yield stress: \( \sigma_y = 1200 \, \text{MPa} \)
- Loads:
- Total lift force: \( F_{\text{lift}} = 3000 \, \text{N} \), distributed parabolically over the chord.
- Constraints:
- Maximum stress: \( \sigma_{\text{max}} \leq \sigma_y = 1200 \, \text{MPa} \)
- Maximum displacement: \( \delta_{\text{max}} \leq 0.15 \, \text{m} \)
- Minimum natural frequency: \( f_{\text{min}} \geq 1.0 \, \text{Hz} \)
- Buckling constraint for spars and skin panels.
- Optimization Variables:
- Thickness options for skin: \( t_{\text{skin}} = \{30, 35, 40, 45\} \, \text{mm} \)
- Thickness options for spars: \( t_{\text{spar}} = \{100, 120, 140, 160\} \, \text{mm} \)
- Fiber angle options: \( \theta = \{0^\circ, 45^\circ, 90^\circ\} \)
- Objective:
- Minimize the total weight of the wing cross-section.
Wing has been simplified into a 2D cross-section divided into four elements:
- Element 1: Skin (top surface, \( 0 \leq x \leq 2 \, \text{m} \), \( y = 0.1 \, \text{m} \))
- Element 2: Skin (bottom surface, \( 0 \leq x \leq 2 \, \text{m} \), \( y = 0 \, \text{m} \))
- Element 3: Spar at \( x = 0.5 \, \text{m} \), \( 0 \leq y \leq 0.1 \, \text{m} \)
- Element 4: Spar at \( x = 1.5 \, \text{m} \), \( 0 \leq y \leq 0.1 \, \text{m} \)
Each element has a thickness \( t_i \) and fiber angle \( \theta_i \). It is assumed that the material properties are uniform across elements for simplicity, except for thickness and fiber angle.
Define the Objective Function
The objective is to minimize the total weight of the wing cross-section. The weight is proportional to the mass, given by:
\[ W = \rho \sum_{i=1}^4 A_i t_i \]
where:
- \( \rho = 1600 \, \text{kg/m}^3 \) is the material density.
- \( A_i \) is the cross-sectional area of element \( i \) (per unit span).
- \( t_i \) is the thickness of element \( i \).
Calculate Areas:
- Skin Elements (1 and 2):
- Length along chord: \( L_{\text{skin}} = 2.0 \, \text{m} \).
- Area per unit span: \( A_{\text{skin}} = L_{\text{skin}} = 2.0 \, \text{m} \).
- Spar Elements (3 and 4):
- Height: \( h = 0.1 \, \text{m} \).
- Area per unit span: \( A_{\text{spar}} = h = 0.1 \, \text{m} \).
Thus:
- \( A_1 = A_2 = 2.0 \, \text{m} \)
- \( A_3 = A_4 = 0.1 \, \text{m} \)
The weight becomes:
\[ W = 1600 \cdot (2.0 \cdot t_1 + 2.0 \cdot t_2 + 0.1 \cdot t_3 + 0.1 \cdot t_4) \cdot 1 \, \text{kg} \]
Assuming a unit span (1 m), mass is computed in kg, which is numerically equivalent to weight in kg for simplicity (ignoring gravitational constant for this context).
Model the Load Distribution
The total lift force \( F_{\text{lift}} = 3000 \, \text{N} \) is distributed parabolically over the top skin (Element 1). The load distribution is:
\[ q(x) = q_0 \sqrt{1 - \left(\frac{x}{c}\right)^2} \]
where \( c = 2.0 \, \text{m} \), and \( q_0 \) is a constant to be determined. The total force is:
\[ F_{\text{lift}} = \int_0^c q(x) \, dx = \int_0^{2.0} q_0 \sqrt{1 - \left(\frac{x}{2}\right)^2} \, dx = 3000 \]
To solve for \( q_0 \):
\[ \int_0^2 \sqrt{1 - \left(\frac{x}{2}\right)^2} \, dx \]
Substitute \( u = x/2 \), \( x = 2u \), \( dx = 2 \, du \):
\[ \int_0^1 q_0 \sqrt{1 - u^2} \cdot 2 \, du = 2 q_0 \int_0^1 \sqrt{1 - u^2} \, du \]
The integral is a quarter-circle:
\[ \int_0^1 \sqrt{1 - u^2} \, du = \frac{\pi}{4} \]
Thus:
\[ 2 q_0 \cdot \frac{\pi}{4} = 3000 \]
\[ q_0 = \frac{3000 \cdot 4}{2 \pi} = \frac{6000}{\pi} \approx 1909.86 \, \text{N/m} \]
So, the distributed load is:
\[ q(x) = \frac{6000}{\pi} \sqrt{1 - \left(\frac{x}{2}\right)^2} \, \text{N/m} \]
This load is applied downward on the top skin (Element 1).
Finite Element Analysis
To analyze the wing, it is modelled as a structure with four elements and apply a simplified finite element approach. Each element has 2 degrees of freedom (DOF) per node (displacement in \( x \) and \( y \)). A coarse mesh with nodes at key points is assumed.
Nodes:
- Node 1: \( (0, 0) \)
- Node 2: \( (0.5, 0) \)
- Node 3: \( (1.5, 0) \)
- Node 4: \( (2.0, 0) \)
- Node 5: \( (0, 0.1) \)
- Node 6: \( (0.5, 0.1) \)
- Node 7: \( (1.5, 0.1) \)
- Node 8: \( (2.0, 0.1) \)
Elements:
- Element 1: Top skin (Nodes 5-6-7-8)
- Element 2: Bottom skin (Nodes 1-2-3-4)
- Element 3: Spar 1 (Nodes 2-6)
- Element 4: Spar 2 (Nodes 3-7)
Boundary Conditions:
- Fixed at \( x = 0 \): Nodes 1 and 5 have zero displacement (\( u_1 = v_1 = u_5 = v_5 = 0 \)).
Stiffness Matrix:
For each element, the stiffness matrix based on its material properties and fiber angle is computed. The stiffness matrix for a composite material in the global coordinate system depends on the fiber angle \( \theta \).
The compliance matrix in the material coordinate system is:
\[ S = \begin{bmatrix}
\frac{1}{E_1} & -\frac{\nu_{12}}{E_1} & 0 \\
-\frac{\nu_{21}}{E_2} & \frac{1}{E_2} & 0 \\
0 & 0 & \frac{1}{G_{12}}
\end{bmatrix} \]
Where:
- \( E_1 = 135 \, \text{GPa} \)
- \( E_2 = 10 \, \text{GPa} \)
- \( G_{12} = 5 \, \text{GPa} \)
- \( \nu_{12} = 0.3 \)
- \( \nu_{21} = \nu_{12} \frac{E_2}{E_1} = 0.3 \cdot \frac{10}{135} \approx 0.0222 \)
Numerically:
\[ S = \begin{bmatrix}
7.407 \times 10^{-12} & -2.222 \times 10^{-12} & 0 \\
-2.222 \times 10^{-12} & 1.0 \times 10^{-10} & 0 \\
0 & 0 & 2.0 \times 10^{-10}
\end{bmatrix} \, \text{Pa}^{-1} \]
To transform to the global coordinate system for fiber angle \( \theta \), we use the transformation matrix:
\[ T = \begin{bmatrix}
\cos^2 \theta & \sin^2 \theta & 2 \cos \theta \sin \theta \\
\sin^2 \theta & \cos^2 \theta & -2 \cos \theta \sin \theta \\
-\cos \theta \sin \theta & \cos \theta \sin \theta & \cos^2 \theta - \sin^2 \theta
\end{bmatrix} \]
The global compliance matrix is:
\[ S_{\text{global}} = T^T S T \]
The stiffness matrix is:
\[ Q = S_{\text{global}}^{-1} \]
\[ T(\theta = 0) = \begin{bmatrix}
1 & 0 & 0 \\
0 & 1 & 0 \\
0 & 0 & 1
\end{bmatrix} \]
\[ S_{\text{global}} = S \]
\[ Q = S^{-1} \]
Calculate the inverse:
\[ S = \begin{bmatrix}
7.407 \times 10^{-12} & -2.222 \times 10^{-12} & 0 \\
-2.222 \times 10^{-12} & 1.0 \times 10^{-10} & 0 \\
0 & 0 & 2.0 \times 10^{-10}
\end{bmatrix} \]
The 2x2 submatrix for \( \sigma_x, \sigma_y \):
\[ S_{2x2} = \begin{bmatrix}
7.407 \times 10^{-12} & -2.222 \times 10^{-12} \\
-2.222 \times 10^{-12} & 1.0 \times 10^{-10}
\end{bmatrix} \]
Determinant:
\[ \text{det} = (7.407 \times 10^{-12}) \cdot (1.0 \times 10^{-10}) - (-2.222 \times 10^{-12})^2 \]
\[ = 7.407 \times 10^{-22} - 4.937 \times 10^{-24} \approx 7.358 \times 10^{-22} \]
Inverse:
\[ Q_{2x2} = \frac{1}{\text{det}} \begin{bmatrix}
1.0 \times 10^{-10} & 2.222 \times 10^{-12} \\
2.222 \times 10^{-12} & 7.407 \times 10^{-12}
\end{bmatrix} \]
\[ Q_{2x2} \approx \begin{bmatrix}
1.359 \times 10^{11} & 3.021 \times 10^9 \\
3.021 \times 10^9 & 1.006 \times 10^{10}
\end{bmatrix} \]
\[ Q_{33} = \frac{1}{2.0 \times 10^{-10}} = 5.0 \times 10^9 \]
\[ Q = \begin{bmatrix}
1.359 \times 10^{11} & 3.021 \times 10^9 & 0 \\
3.021 \times 10^9 & 1.006 \times 10^{10} & 0 \\
0 & 0 & 5.0 \times 10^9
\end{bmatrix} \, \text{Pa} \]
For other angles,\( T \) would be calculated and transformation will be performed, but for simplicity, it will proceed with \( \theta = 0^\circ \) and test other angles later.
Element Stiffness:
For a 2D element, the stiffness matrix is computed using:
\[ K_e = t \int_A B^T Q B \, dA \]
Where \( B \) is the strain-displacement matrix. For simplicity, assumptions is that each element behaves as a 1D bar or plate under tension/compression for stress calculations, and deflections using beam theory are calculated.
Apply Loads and Boundary Conditions
The distributed load \( q(x) \) is applied to the top skin. The load on nodes 5, 6, 7, and 8 by integrating over segments is approximated:
\[ F_i = \int_{x_i}^{x_{i+1}} q(x) \, dx \]
Divide the chord into segments:
- Node 5 to 6: \( x = 0 \) to \( 0.5 \)
- Node 6 to 7: \( x = 0.5 \) to \( 1.5 \)
- Node 7 to 8: \( x = 1.5 \) to \( 2.0 \)
Calculate the force per segment:
\[ F = \int_a^b \frac{6000}{\pi} \sqrt{1 - \left(\frac{x}{2}\right)^2} \, dx \]
For segment \( x = 0 \) to \( 0.5 \):
\[ F_{5-6} = \frac{6000}{\pi} \int_0^{0.25} \sqrt{1 - u^2} \cdot 2 \, du \]
This integral is complex, so the average load will be approximated:
\[ q_{\text{avg}} = \frac{1}{0.5} \int_0^{0.5} \frac{6000}{\pi} \sqrt{1 - \left(\frac{x}{2}\right)^2} \, dx \]
For simplicity, assume uniform distribution for calculation:
\[ F_{\text{total}} = 3000 \, \text{N} \]
Distribute equally among nodes 6, 7, 8 (Node 5 is fixed):
\[ F_6 = F_7 = F_8 = \frac{3000}{3} = 1000 \, \text{N} \text{ (downward)} \]
Force vector (y-direction, negative):
\[ F = [0, 0, 0, 0, 0, 0, 0, -1000, 0, -1000, 0, -1000]^T \]
Simplified Structural Analysis
To make the problem tractable manually, wing is modelled as a cantilever beam with two spars and skins, and compute stresses, displacements, buckling, and frequency.
Displacement:
Treat the wing as a beam with fixed end at \( x = 0 \). The moment of inertia \( I \) depends on the thicknesses:
\[ I = \sum_i \frac{t_i h_i^3}{12} + t_i h_i d_i^2 \]
For simplicity, assume the spars dominate bending resistance:
\[ I \approx t_3 \cdot \frac{h^3}{12} + t_4 \cdot \frac{h^3}{12} \]
\[ I = (t_3 + t_4) \cdot \frac{0.1^3}{12} = (t_3 + t_4) \cdot \frac{0.001}{12} \]
Maximum deflection for a distributed load:
\[ \delta_{\text{max}} = \frac{q L^4}{8 E I} \]
Approximate \( q = 1500 \, \text{N/m} \) (total 3000 N over 2 m):
\( L = 2.0 \, \text{m} \), \( E = E_1 = 135 \times 10^9 \, \text{Pa} \)
\[ \delta_{\text{max}} = \frac{1500 \cdot 2^4}{8 \cdot 135 \times 10^9 \cdot I} = \frac{1500 \cdot 16}{1.08 \times 10^{12} \cdot I} = \frac{24000}{1.08 \times 10^{12} \cdot I} \]
Constraint: \( \delta_{\text{max}} \leq 0.15 \, \text{m} \).
Stress:
Maximum bending stress in the spars:
\[ \sigma = \frac{M c}{I} \]
Maximum moment at the fixed end:
\[ M_{\text{max}} = \frac{q L^2}{2} = \frac{1500 \cdot 2^2}{2} = 3000 \, \text{N·m} \]
\[ c = \frac{h}{2} = 0.05 \, \text{m} \]
\[ \sigma = \frac{3000 \cdot 0.05}{I} = \frac{150}{I} \, \text{Pa} \]
Constraint: \( \sigma \leq 1200 \times 10^6 \, \text{Pa} \).
Buckling (Spars):
Critical buckling load for a spar (fixed-fixed column):
\[ P_{\text{cr}} = \frac{\pi^2 E I_{\text{spar}}}{L^2} \]
\[ I_{\text{spar}} = t_i \cdot \frac{h^3}{12} = t_i \cdot \frac{0.1^3}{12} = t_i \cdot \frac{0.001}{12} \]
\[ L = h = 0.1 \, \text{m} \]
\[ P_{\text{cr}} = \frac{\pi^2 \cdot 135 \times 10^9 \cdot t_i \cdot \frac{0.001}{12}}{0.1^2} \]
\[ = \frac{9.8696 \cdot 135 \times 10^9 \cdot t_i \cdot 0.00008333}{0.01} \]
\[ = 1.111 \times 10^{10} t_i \, \text{N} \]
Load on each spar (approximate):
\[ P = \frac{3000}{2} = 1500 \, \text{N} \]
\[ P_{\text{cr}} \geq 1500 \]
\[ 1.111 \times 10^{10} t_i \geq 1500 \]
\[ t_i \geq 1.35 \times 10^{-7} \, \text{m} \]
This is very small, indicating buckling is not a limiting factor for reasonable thicknesses, but panel buckling for skin is checked.
Skin Buckling:
For a plate under compression:
\[ \sigma_{\text{cr}} = \frac{k \pi^2 E}{12 (1 - \nu^2)} \left(\frac{t}{b}\right)^2 \]
Assume \( k = 4 \), \( \nu = 0.3 \), \( b = 0.5 \, \text{m} \) (distance between spars):
\[ \sigma_{\text{cr}} = \frac{4 \pi^2 \cdot 135 \times 10^9}{12 (1 - 0.3^2)} \cdot \frac{t^2}{0.5^2} \]
\[ = \frac{4 \cdot 9.8696 \cdot 135 \times 10^9}{12 \cdot 0.91} \cdot \frac{t^2}{0.25} \]
\[ = 1.953 \times 10^{11} t^2 \, \text{Pa} \]
\[ \sigma \leq 1.953 \times 10^{11} t^2 \]
Frequency:
Fundamental frequency of a cantilever beam:
\[ f = \frac{1}{2\pi} \sqrt{\frac{k}{m}} \]
Approximate stiffness \( k \):
\[ k = \frac{3 E I}{L^3} \]
Mass:
\[ m = \rho \sum A_i t_i \cdot L \]
Manual Optimization
Combinations of thicknesses and fiber angles to minimize weight while satisfying constraints will be tested.
Trial 1:
- \( t_1 = t_2 = 0.030 \, \text{m} \) (30 mm)
- \( t_3 = t_4 = 0.100 \, \text{m} \) (100 mm)
- \( \theta_1 = \theta_2 = \theta_3 = \theta_4 = 0^\circ \)
Weight:
\[ W = 1600 \cdot (2.0 \cdot 0.030 + 2.0 \cdot 0.030 + 0.1 \cdot 0.100 + 0.1 \cdot 0.100) \]
\[ = 1600 \cdot (0.06 + 0.06 + 0.01 + 0.01) = 1600 \cdot 0.14 = 224 \, \text{kg} \]
Moment of Inertia:
\[ I = (0.100 + 0.100) \cdot \frac{0.001}{12} = 0.2 \cdot \frac{0.001}{12} \approx 1.667 \times 10^{-5} \, \text{m}^4 \]
Displacement:
\[ \delta_{\text{max}} = \frac{24000}{1.08 \times 10^{12} \cdot 1.667 \times 10^{-5}} \]
\[ = \frac{24000}{1.8 \times 10^7} \approx 0.00133 \, \text{m} \]
\[ 0.00133 \leq 0.15 \text{ (Satisfied)} \]
Stress:
\[ \sigma = \frac{150}{1.667 \times 10^{-5}} \approx 8.998 \times 10^6 \, \text{Pa} = 8.998 \, \text{MPa} \]
\[ 8.998 \leq 1200 \text{ (Satisfied)} \]
Skin Buckling:
\[ \sigma_{\text{cr}} = 1.953 \times 10^{11} \cdot (0.030)^2 = 1.953 \times 10^{11} \cdot 0.0009 \]
\[ = 1.758 \times 10^8 \, \text{Pa} = 175.8 \, \text{MPa} \]
\[ 8.998 \leq 175.8 \text{ (Satisfied)} \]
Frequency (approximate):
\[ m = 1600 \cdot 0.14 \cdot 2 = 448 \, \text{kg} \]
\[ k = \frac{3 \cdot 135 \times 10^9 \cdot 1.667 \times 10^{-5}}{2^3} = \frac{6.75 \times 10^6}{8} = 8.4375 \times 10^5 \, \text{N/m} \]
\[ f = \frac{1}{2\pi} \sqrt{\frac{8.4375 \times 10^5}{448}} \approx \frac{1}{6.283} \sqrt{1883.93} \]
\[ \sqrt{1883.93} \approx 43.41 \]
\[ f \approx \frac{43.41}{6.283} \approx 6.91 \, \text{Hz} \]
\[ 6.91 \geq 1.0 \text{ (Satisfied)} \]
All constraints are satisfied. Weight = 224 kg.
Trial 2:
- \( t_1 = t_2 = 0.035 \, \text{m} \)
- \( t_3 = t_4 = 0.120 \, \text{m} \)
- \( \theta = 0^\circ \)
Weight:
\[ W = 1600 \cdot (2.0 \cdot 0.035 + 2.0 \cdot 0.035 + 0.1 \cdot 0.120 + 0.1 \cdot 0.120) \]
\[ = 1600 \cdot (0.07 + 0.07 + 0.012 + 0.012) = 1600 \cdot 0.164 = 262.4 \, \text{kg} \]
Trial 3:
- \( t_1 = t_2 = 0.030 \, \text{m} \)
- \( t_3 = t_4 = 0.100 \, \text{m} \)
- \( \theta_1 = \theta_2 = 45^\circ \), \( \theta_3 = \theta_4 = 0^\circ \)
Recalculate stiffness for skin elements with \( \theta = 45^\circ \):
\[ \theta = 45^\circ, \cos \theta = \sin \theta = \frac{\sqrt{2}}{2} \approx 0.7071 \]
\[ T = \begin{bmatrix}
0.5 & 0.5 & 1 \\
0.5 & 0.5 & -1 \\
-0.5 & 0.5 & 0
\end{bmatrix} \]
Compute \( S_{\text{global}} = T^T S T \), then \( Q \). This is complex, so we note that \( \theta = 45^\circ \) reduces effective stiffness, potentially increasing stress and displacement. Since Trial 1 satisfies constraints, it is selected.
Select Optimal Solution
After testing combinations, Trial 1 yields:
- Thicknesses: \( t_1 = t_2 = 30 \, \text{mm} \), \( t_3 = t_4 = 100 \, \text{mm} \)
- Fiber Angles: \( \theta_1 = \theta_2 = \theta_3 = \theta_4 = 0^\circ \)
- Weight: 224 kg
- Constraints:
- Stress: 8.998 MPa < 1200 MPa
- Displacement: 0.00133 m < 0.15 m
- Buckling: Satisfied
- Frequency: 6.91 Hz > 1.0 Hz
This solution minimizes weight while satisfying all constraints.
Solution
The optimal material distribution for the wing cross-section is:
- Top and Bottom Skin: Thickness = 30 mm, Fiber Angle = 0°
- Spars: Thickness = 100 mm, Fiber Angle = 0°
- Total Weight: 224 kg
All structural constraints are satisfied, ensuring the wing can withstand the applied lift load without excessive stress, displacement, buckling, or insufficient stiffness.
Repository