Divergence Structural Analysis of an Aircraft Tapered-Unswept Wing Using Aerodynamic Strip Theory for the Case of Flight at Sea Level

This study investigates the utilize of the iteration of matrix numerically to get the divergence speed of an aircraft tapered-unswept wing via using method of strip theory (without finite span correction and with finite span correction). The 2-D fluid flow across the wing airfoil section is assumed by the aerodynamic strip theory and then integrated over wing span. In the present work, first the tapered unswept wing is divided into four Multhopp’s stations. Then, elastic properties of this wing in relation to torsional influence coefficients related with this four Multhopp’s stations have been evaluated. Second, equations for equilibrium are constructed as integral equations. Then, an appropriate aerodynamic theory in the form of strip theory was addressed, as well as the technique of solution for determining the divergence speed. Finally, using strip theory, the integral equation has been expressed in matrix form. Iteration of matrix has been done using MATLAB environment to reach for the solution to converge. Also, an increase of 15% in torsional stiffness of aircraft wing has been considered to illustrate its effects on torsional divergence speed of aircraft wing. The attained results are going to be helpful in understanding of wing instability for modern aircraft designer.

Introduction

Aero-Elasticity Description

Aero-elasticity is the study of the interactions of inertia, elastic and aerodynamic forces in a flexible structure and the phenomenon that might occur as a result [1]. The effect of inertial forces is introduced in dynamics. Elasticity predicts the form of an elastic body under a given load. classical aerodynamic approaches estimate the forces working on a body of a certain shape [2]. This subject is best summarized by the classic Collar aeroelastic triangle [3], shown in Fig. 1, which illustrates how the main fields of stability and control, structural dynamics, and static aeroelasticity each result from the interaction of two of the three types of force. Nevertheless, for dynamic aeroelastic effects to occur, all three forces must interact. Aeroelastic phenomenon are often divided into two categories: static and dynamic. The nonoscillatory effects of aerodynamic forces acting on the flexible aircraft structure are taken into account by static aeroelasticity. The oscillatory effects of aeroelastic interactions are studied in dynamic aeroelasticity [1]. The static aeroelastic phenomena involve divergence and control reversal, for instance. While flutter and buffeting are dynamic aeroelasticity.

Fig. 1. The aeroelastic triangle of Collar [1].

Divergence Description

Divergence is a static aeroelastic phenomena that includes the interplay of aerodynamic and elastic forces. The most prevalent type of divergence is wing torsional divergence, which happens when the local angle of the wing grows to the point where structural collapse occurs [4]. This occurs when a wing experiences a torsional moment produced by lifting forces, which causes the wing to twist about its elastic axis. The local angle of attack of the wing is increased as a result of this twist. additional lift is created as a result, which in turn causes additional twist. The elastic reactive forces and the aerodynamic forces will eventually either achieve equilibrium or the reactive forces will be overcome. The local angle of attack, as well as the aerodynamic forces, will continue to rise until structural collapse occurs. Thus, at the divergence speed, the increase in aerodynamic torsional moment caused by a twist angle increase is equal to the equivalent increase in elastic restoring torque. If the airplane surpasses this speed, the increase in aerodynamic moment is larger than the rise in elastic restoring torque [4].

Wing Geometry of the Case Study

Because the wing in consideration is tapered, the chord changes over the wing span as shown in Fig. 2.

Fig. 2. Wing Multhopp’s stations of the case study.

The wing root chord length, C_r = 5.588 m.

The wing tip chord length, C_t = 2.794 m.

The span at any span length, y_i = L × cos ϕ.

The wing span, y = 22.86 m.

• All dimensions are in meter.
• Aerodynamic center = 25% chord.
• Elastic axis = 35% chord.
• From relation: ,

$$\begin{array}{r}{\mathrm{\varnothing }}_1={22.5}^{\circ },{\mathrm{\varnothing }}_2={45}^{\circ },{\mathrm{\varnothing }}_3={67.5}^{\circ },{\mathrm{\varnothing }}_4={90}^{\circ }\end{array}$$

• The span at any span length, y_i = L × cos ϕ

Then:

• Point 4 at zero from wing root.
• Point 3 at 4.8601 m from wing root.
• Point 2 at 8.9803 m from wing root.
• Point 1 at 11.7333 m from wing root.

Structure Properties of the Wing

It is assumed that the wing under investigation is perfectly elastic. Means that when the external loads are removed, the wing structure retains its original shape. Experiments on airplane structures demonstrated that, within certain limitations, force and deflection are linearly connected. Even if the material that makes up the wing structure is stressed at a low level, elastic buckling in the skin of an aircraft wing structure can cause a discontinuity in the force deflection diagram [5]. As a result, the wing structure’s elastic properties are presented in the range under the elastic buckling point.

Mathematical Approaches

Influence Coefficients Apprehensible

The influence coefficients conception is used throughout the study, which takes into account the wing structural deflections caused by various loads. I’ll use the Fig. 3 general situation as an example. According to this method, total angular and linear deformation of each point can be represented as the summation of deformations at that point caused by individual forces and moments. This can be expressed via the superposition principle, which serves as the foundation of linear analysis of systems [5]. On Fig. 3, sign Q is attributed to the arbitrary force or moment known as “generalized forces”. Likewise, sign q corresponds to angular or linear deflection of the point of application of each generalized force and is referred to as “generalized coordinates”. Generalized coordinates are quantities that illustrates the system’s possible independent displacements. As a result, they must not violate the system’s geometric restrictions. Because of this requirement, deformation of the wing at point of connection with the fuselage has to be zero.

Fig. 3. A number of discrete forces and moments are applied to an elastic body [5].

According to superposition concept, the deflection of point of application i^th generalized force caused by n generalized forces provided by:

(1) $$\begin{array}{r}{\text{q}}_{\text{i}}={\sum }_{\text{j}=1}^{\text{n}}{\complement }_{\text{ij}}{\text{Q}}_{\text{j}},(\text{i}=1,2,3,\dots .\text{n})\end{array}$$

where C_ij referred to the coefficients of flexibility influence. The aforementioned equation can be written as follows in matrix notation:

$$\begin{array}{r}\left[\begin{array}{c}{q}_1\\ q_2\\ .\\ .\\ .\\ q_n\end{array}\right]=\left[\begin{array}{c}{\complement }_{11}{\complement }_{12}...{\complement }_{1n}\\ {\complement }_{21}{\complement }_{22}...{\complement }_{2n}\\ ......\\ ......\\ ......\\ {\complement }_{n1}....{\complement }_{nn}\end{array}\right]\left[\begin{array}{c}{Q}_1\\ Q_2\\ .\\ .\\ .\\ Q_n\end{array}\right]\end{array}$$

It can be written as follows in short matrix notation:

$$\begin{array}{r}\left\{q\right\}=[C]\left\{q\right\}\end{array}$$

Symmetry is an important feature of influence coefficients and associated matrices. This characteristic is expressed as:

$$\begin{array}{r}{C}_{ij}=C_j\end{array}$$

Consider matrix [C]’s properties, which are relevant to Fig. 3. The following four matrices, each having a different type of influence coefficient, can be created when n = 4:

$$\begin{array}{r}\left[\complement \right]=\left[\begin{array}{cccc}{\complement }_{11}^{\delta \delta }& {\complement }_{12}^{\delta \delta }& {\complement }_{13}^{\delta \alpha }& {\complement }_{14}^{\delta \alpha }\\ {\complement }_{21}^{\delta \delta }& {\complement }_{22}^{\delta \delta }& {\complement }_{23}^{\delta \delta }& {\complement }_{11}^{\delta \delta }\\ {\complement }_{31}^{\alpha \delta }& {\complement }_{32}^{\delta \delta }& {\complement }_{33}^{\alpha \alpha }& {\complement }_{11}^{\delta \delta }\\ {\complement }_{41}^{\alpha \delta }& {\complement }_{42}^{\delta \delta }& {\complement }_{43}^{\delta \delta }& {\complement }_{44}^{\alpha \alpha }\end{array}\right]\end{array}$$

There are four categories of elements:

∁^δδ: Linear deformation at i as a result of unit force at j.

∁^αδ: Linear deflection at i caused by unit moment at j.

∁^αα: Linear deformation at i caused by unit moment at j.

∁^δα: Linear deflection at i as a result of unit force at j.

The reciprocal theorem of Betti has to be applied for the matrix [C] to be symmetric:

$$\begin{array}{r}{\complement }_{ij}^{\delta \delta }={\complement }_{ji}^{\delta \delta },{\complement }_{ij}^{\alpha \alpha }={\complement }_{ji}^{\alpha \alpha },{\complement }_{ij}^{\delta \alpha }={\complement }_{ji}^{\alpha \delta }\end{array}$$

Strain Energy Expression According to Influence Coefficients

In order to apply the energy approach to an aeroelastic system, the following strain energy formulas in relation to influence coefficients can be constructed with reference to Fig. 3:

(2) $$\begin{array}{r}U=\frac{1}{2}{\sum }_{i=1}^n{Q}_i{q}_i\end{array}$$

Replacing (1) into (2) provides the strain energy expression with respect to flexibility influence coefficients:

(3) $$\begin{array}{r}U=\frac{1}{2}{\sum }_{i=1}^n{\sum }_{j=1}^n{\complement }_{ij}{Q}_i{Q}_j\end{array}$$

Consider the wing section depicted in Fig. 4 to formulate the strain energy as an integral form. The wing span-wise coordinate y and the tangential coordinate s are used to identify a point on the cross section. Tangential coordinates are determined positive in anticlockwise direction for the peripheral skin and positive with the positive direction of Z-axis for internal webs [6].

Fig. 4. Loaded cross section of shell beam [6].

As shown in Fig. 4, a torque T(y) acting in the positive direction acts up on the cross section at y distance from the origin. At section’s shear center, T(y)’s application point is located. Shear flow, which is represented as q, is positive with S-direction.

If the beam is free to deform during the application of twisting moments, the strain energy is totally attributable to shear stress and provided as:

(4) $$\begin{array}{r}U=\frac{1}{2G}{\int }_0^l{\int }_{C.S}{q}^2\frac{ds}{t}d\lambda \end{array}$$

Twist angle of the beam caused by specified distributed applied torque T(y) is determined by using Castigliano’s theorem to (4) as follows: where (λ, s) is distribution of shear flow caused by applied torques.

(5) $$\begin{array}{r}\theta \left(y\right)=\frac{\mathrm{\partial }U}{\mathrm{\partial }T}=\frac{1}{G}{\int }_0^l{\int }_{C.S}q(\lambda ,s)\frac{\mathrm{\partial }q}{\mathrm{\partial }T}\frac{ds}{t}d\lambda \end{array}$$

$\mathrm{\partial }q/\mathrm{\partial }T$ is distribution of shear flow caused by unit torque.

T = 1, performed on spanwise section.

Assuming for the shear flow distribution caused by T = 1 is represented as (λ, s), then: replacing (6) into (5), produce: where $J={\displaystyle \frac{1}{{\int }_{c.s}{\vartheta }^{2}ds/t}}$ is the Cross-sectional rigidity of the beam in torsion. Additionally, GJ is known as the beam’s torsional stiffness.

$$\begin{array}{r}q(\lambda ,s)=T\left(\lambda \right)\vartheta (\lambda ,s)\end{array}$$

(6) $$\begin{array}{r}\frac{\mathrm{\partial }q}{\mathrm{\partial }T}=\{\begin{array}{l}\vartheta (\lambda ,s)(y>\lambda )\\ 0(y<\lambda )\end{array}\end{array}$$

(7) $$\begin{array}{r}\theta \left(y\right)={\int }_0^y\frac{T\left(\lambda \right)d\lambda }{GJ}\end{array}$$

When we differentiate (7) with regard to y, we obtain the following relationship among the applied torsional moment and rate of twist:

(8) $$\begin{array}{r}d\theta \left(y\right)/dy={\theta }^{{}^{\prime }}=\frac{T}{GJ}\end{array}$$

Coefficients of Torsional Influence

Consider the cantilever wing in Fig. 5, which is subjected to a unit torque force. At a distance η from the origin, a unit torque about the elastic axis is exerted, and the angular deflection resulting at y is denoted as C^θθ (y, η).

Fig. 5. Cantilever wing subjected to unit torque [6].

As previously stated:

$$\begin{array}{r}q(\lambda ,s)=\vartheta (\lambda ,s)(0<\lambda <\eta )\end{array}$$

$$\begin{array}{r}\frac{\mathrm{\partial }q}{\mathrm{\partial }T}=\vartheta (\lambda ,s)(0<\lambda <y)\end{array}$$

Thus, for η > $y$: and for η < $y$:

(9) $$\begin{array}{r}{\complement }^{\theta \theta }(y,\eta )={\int }_0^y\frac{d\lambda }{GJ}\end{array}$$

(10) $$\begin{array}{r}{\complement }^{\theta \theta }(y,\eta )={\int }_0^{\lambda }\frac{d\lambda }{GJ}\end{array}$$

If the distribution of shear flow $\vartheta$ (s, λ) caused by a unit torque is given, the constant of torsional J, stated in (9), (10), can be estimated at any section of the beam. This necessitates knowledge of the wing skin thickness, flange, web and stringer thickness, and so on at each section of the wing, as well as the shear modulus values at each section. According to reference [6], torsional rigidity curve GJ has been calculated and is presented in Fig. 6.

Fig. 6. (a) Curve of shear stiffness, torsional and bending [6]. (b) Curve of torsional stiffness (modified to SI unit) [5].

Matrix of Torsional Influence Coefficients_Numeri- cal Estimation

According to Fig. 2, the wing is segmented into four stations over the semi-span. These specific locations, known as Multhopp’s stations, were chosen to make computing the aerodynamic matrices easier. These stations have been determined to be at y₄ = 0, y₃ = 4.8601, y₂ = 8.9803, and y₁ = 11.7333 meters from the root. Because four wing stations are assumed in this work, the obtained matrix of torsional coefficients is (4 × 4) square matrix, illustrated as follows:

$$\begin{array}{r}\left[C^{\theta \theta }\right]=\left[\begin{array}{c}{C}^{11}{C}^{12}{C}^{13}{C}^{14}\\ C^{21}{C}^{22}{C}^{23}{C}^{24}\\ C^{31}{C}^{32}{C}^{33}{C}^{34}\\ C^{41}{C}^{42}{C}^{43}{C}^{44}\end{array}\right]\end{array}$$

These coefficients’ values have been calculated from (9) and (10) to be:

[$C^{\theta \theta }$]

(11) $$\begin{array}{r}=\left[\begin{array}{c}37.55416.5166.9430\\ 16.51616.5166.9430\\ 6.9436.9436.9430\\ 0000\end{array}\right]\times {10}^{-8}rad.N.m^{-1}\end{array}$$

Divergence Structural Analysis

Divergence, a static instability phenomenon, is caused by the combination of elastic and aerodynamic forces. As a result, calculating the matrix of aerodynamic coefficients and the matrix of structural stiffness is required for the divergence speed analysis. The current study is concerned with the torsional divergence of tapered un-swept wing, which is the most prevalent problem in aero-elasticity.

Take the simple un-swept wing at incidence () with twist center behind the center of aerodynamic as an example. The pressure distribution with major loads situated towards the nose causes the wing to twist in the nose up direction. Due to the structure’s imperfect rigidity, it twists and changes in form in relation to the wing root section. As seen in Fig. 7, it twists about wing flexural center also known torsional axis. This twist increases the effective incidence to (α + θ), which causes a lift increment that acts forward of the torsional axis. As a result, the effect is statically unstable, because, the more it twists, the larger the torsional moment tending to cause it to twist, of course resisted via elastic forces because of twist even more. The twisting propensity is, of course, prevented by elastic forces owing to the structure’s stiffness. This twist resistance grows fast with twist, or strain, till it balances aerodynamic moment of twist and equilibrium is obtained. Nevertheless, as speed rises, the aerodynamic forces (L) rise quickly, which causes the twisting moment (Mo) to rise as well. Speed has no influence on the elastic stiffness, thus as the speed rises, so does the amount of twist. At some point, the elastic restoring torque is only barely enough to counterbalance the twisting moment, and equilibrium is only attained with the wing breaking point. This is a critical speed known as wing torsional divergence speed [5]. Once the wing’s torsional divergence speed is attained, the increase in aerodynamic moment of torsional caused by an arbitrary increase in twist angle equals the increase in elastic restoring torque. The wing is statically unstable if the speed arrives at the torsional divergence speed, in which, the increase in the aerodynamic torsional moment exceeds the increase in the elastic restoring torque. Any rise in speed over this value will lead to structural collapse, and the wing will break off [5].

Fig. 7. Divergence speed estimation of the wing [5].

Equations of Equilibrium

The following assumes are made for simplicity’s sake:

1. Un-swept wings are distinguished via an Elastic-axis that is vertical to the symmetry plane of aircraft.
2. The wing’s chordwise sections maintains rigid; camber bending is to be negligible.

By connecting the rate of twist to the applied torque as previously mentioned, the differential equation of torsional aero-elastic equilibrium of un-swept wing about its elastic axis is represented using (8) as follow:

$$\begin{array}{r}\frac{d\theta \left(y\right)}{dy}={\theta }^{{}^{\prime }}=\frac{T}{GJ}\end{array}$$

You might rewrite this as:

(12) $$\begin{array}{r}\frac{d}{dy}\left(GJ\frac{d\theta }{dy}\right)=\frac{dT}{dy}=-t\left(y\right)\end{array}$$

where (y) is distribution of elastic twist.

Take a slender straight wing that is being affected by aerodynamic and inertial forces, as in Fig. 8.

Fig. 8. Slender straight wing [5].

The utilized torque per unit span t(y) is provided in Fig. 8 by:

(13) $$\begin{array}{r}t(y)=q{c}_{1}ce+q{c}^2{C}_{m_{ac}}-Nmgd\end{array}$$

where $c_l$ is the local lift coefficient.

$C_{m_{ac}}$ is local coefficients of moment about center of aerodynamic.

mg is weight of wing per unit span.

N is the load factor perpendicular to wing surface. At level flight, N equals 1.

When we combine (12) and (13), we get the following differential equation of equilibrium: or

(14) $$\begin{array}{r}\frac{d}{dy}\left(GJ\frac{d\theta }{dy}\right)=Nmgd-q{c}_{l}ce-q{c}^2{C}_{m_{ac}}\end{array}$$

(15) $$\begin{array}{r}\frac{d}{C_{m_{ac}}dy}\left(GJ\frac{d\theta }{dy}\right)+q{c}_{l}ce=Nmgd-q{c}^2\end{array}$$

$\theta \left(0\right)=0;{\theta }^{{}^{\prime }}\left(l\right)=0.$ Are the boundary conditions.

By using Castiglione’s theorem to the energy equation, the wing’s torsional deflection is calculated at any spanwise position y caused by torque t applied at span wise position η.

(16) $$\begin{array}{r}\theta \left(y\right)={\int }_0^l{C}^{\theta \theta }(y,\eta )t\left(y\right)d\eta \end{array}$$

When (14) is added to (16), we get:

(17) $$\begin{array}{r}\theta \left(y\right)={\int }_0^l{C}^{\theta \theta }(y,\eta )\{[c_{l}ce+c^2{C}_{m_{ac}}]q-Nmgd\}d\eta \end{array}$$

The angle of attack may be thought of as a superposition of an elastic twist and a rigid angle.

(18) $$\begin{array}{r}\alpha \left(y\right)={\alpha }^r\left(y\right)+\theta \left(y\right)\end{array}$$

Additionally, local coefficient of lift may be expressed as follows: where $c_l$($y$) is local angle of attack determined from zero lift excepting elastic twist.

(19) $$\begin{array}{r}{c}_l\left(y\right)=c_l^r\left(y\right)+c_l^e\left(y\right)\end{array}$$

$c_l^r\left(y\right)$is distribution of local coefficient of lift caused by rigid twist, αr (y).

$c_l^e\left(y\right)$is distribution of local coefficient of lift due to elastic twist.

The following integral equation is obtained by substituting (19) in (15):

(20) $$\begin{array}{r}\frac{d}{dy}\left(GJ\frac{d\theta }{dy}\right)+qce{c}_l^e=-qce{c}_l^r-q{c}^2{C}_{m_{ac}}+Nmgd\end{array}$$

Similarly, we derive the following integral equation by replacing (19) in (17): where

(21) $$\begin{array}{r}\theta \left(y\right)=q{\int }_0^l{C}^{\theta \theta }(y,\eta )ec{c}_l^{e}d\eta +f\left(y\right)\end{array}$$

$$\begin{array}{r}f\left(y\right)={\int }_0^l{C}^{\theta \theta }(y,\eta )(ec{c}_l^{r}q+q{c}^2{C}_{m_{ac}}-Nmgd)d\eta \end{array}$$

This equation represents the required governing integral equation.

Definition of Aerodynamic Strip Theory

There are several approaches to modeling a wing’s spanwise lift distribution. This section will provide the strip theory for a discretized wing, which is the most basic. According to this approach, the wing is made up of a number of chordwise elements famed “strips,” and it is supposed that the lift coefficient on each of these pieces is proportional to the local angle of incidence α(y), and that the lift on one segment has no effect on other. In its most basic version, root and tip effects, as well as compressibility effects, are neglected. In reality, these presumptions predict a low air speed (M < 0.3) and a high aspect ratio (AR ≥ 6) for the wing. Take note, drag calculations cannot be performed using strip theory [1]. The motion of a fluid can be framed as a 2-D Problems for slender forms in aerodynamic analysis. Setting a given aerodynamic problem into 2-D crosswise direction in comparison with that of longitudinal one. The strip theory is based on the assumption that the 2-dimensional aerodynamic coefficients for added weight for each strip may be calculated and then integrated throughout the length of the body to yield 3-Dimention coefficients [5].

In (21), both θ(y) and $c_l^e\left(y\right)$ are treated as unknown functions. Hence, the situation becomes undefined mathematically. When a relation among the two unknowns is given, the problem becomes mathematically definite. This second relationship is provided by a suitable aerodynamic theory. The aerodynamic theory assumes a linear relationship between the incidence angle and lift distribution, which may be expressed by:

(22) $$\begin{array}{r}(y)=\mathrm{\Theta }(ccl)\end{array}$$

Θ is linear operator used to the lift distribution $c{c}_l$ to produce the required incidence distribution α(y).

According to the definition of strip theory, Θ is just where a₀ is the local 2-D slope of curve of lift coefficient. c is chord of the wing.

(23) $$\begin{array}{r}\mathrm{\Theta }=\frac{1}{a_{0}c}\end{array}$$

Matrix-Solution Method Based on Strip Theory

Equations (21) and (22) constitute the foundation for predicting the elastic twist and distribution of lift of a straight wing include straight elastic axis.

The dynamic pressure q_d with the least Eigen value, which is derived from equilibrium equation in its homogeneous-integral form, are used to calculate the torsional divergence speed of a 3-D wing [5]. Equation (20)’s homogeneous form is:

(24) $$\begin{array}{r}\frac{d}{dy}\left(GJ\frac{d\theta }{dy}\right)+q_{d}ce{c}_l^e=0\end{array}$$

Eq. (21)’s homogeneous form is:

(25) $$\begin{array}{r}\theta \left(y\right)=q_d{\int }_0^l{C}^{\theta \theta }(y,\eta )ec{c}_l^{e}d\eta \end{array}$$

To determine the divergence speed, (24) or (25) can be used with (22). In most practical airplanes problems numerical solutions are required. The integral (25) using matrix techniques is helpful in this case. Equation (25) may be expressed using strip theory as: where ${\displaystyle \frac{d{c}_l}{d\alpha }}$ is the slope of the effective lift curve modified for aspect ratio. Equation (26) may be illustrated on matrix form as follows: where $\left[A\right]={\displaystyle \frac{1}{{\displaystyle \frac{d{c}_l}{d\alpha }}}}\left[1/C\right]$

(26) $$\begin{array}{r}\frac{{\text{c}}_{\text{l}}^{\text{e}}}{\frac{{\text{d}}_{\text{cl}}}{{\text{d}}_{\mathrm{\alpha }}}}={\text{q}}_{\text{d}}{\int }_0^{\text{l}}{\text{C}}^{\mathrm{\theta }\mathrm{\theta }}(\text{y},\mathrm{\eta })\text{ec}{\text{c}}_{\text{l}}^{\text{e}}\text{d}\mathrm{\eta }\end{array}$$

(27) $$\begin{array}{r}\left[A\right]\left\{c{c}_l^e\right\}=q_d\left[E\right]\left\{c{c}_l^e\right\}\end{array}$$

(28) $$\begin{array}{r}E=\left[C^{\theta \theta }\right]\left[e\right]\left[\overline{W}\right]\end{array}$$

where $\left[\overline{W}\right]$ is weighting matrix.

The governing (27) is shown on matrix form to calculating divergence speed numerically.

Numerical Evaluation

Equation (27), which governs the calculation of unswept-tapered wing’s divergence speed, has been formulated and is provided in matrix form. This equation is utilized to calculate divergence speed of wing when matrix $[C^{\theta \theta }]$ has been calculated numerically and is given by (11). Multhopp’s quadratic formula is used to calculate the weighting matrix [$\overline{W}$]. This formula’s use is described in ref [6]. When it is demonstrated that for a specific four-station configuration used in this study, the shape of a diagonal matrix is provided via (31). Equation (27) is solved using a matrix-iteration approach via strip theory with and without corrections of finite span.

Calculation of the Matrix for the Wing Chord

Wing chord at any span $y_i$ may be calculated using the subsequent relationships based on triangle similarity as seen on Fig. 9.

$$\begin{array}{r}\frac{C_i}{C_r}=\frac{y-y_i}{y}\end{array}$$

(29) $$\begin{array}{r}{C}_i=C_r\left(\frac{y-y_i}{y}\right)=C_r(1-\frac{y_i}{y})\end{array}$$

Fig. 9. Geometry of the wing half span platform.

If the previous numerical values are substituted, the following will result:

(29a) $$\begin{array}{r}{C}_i=5.588(1-\frac{y_i}{22.86})\end{array}$$

The following diagonal matrix is obtained by calculating the chord $C_i$ at stations 1, 2, 3, and 4, then making the results into a matrix:

$$\begin{array}{r}\left[C\right]=\left[\begin{array}{c}2.7198000\\ 03.392800\\ 004.39990\\ 0005.588\end{array}\right]\end{array}$$

Eccentricity Calculation

The aerodynamic center (A.C.) of the wing is considered to be one-quarter of the chord taken from the leading edge (0.25 $\times$ chord, C), and (E.A.) located at 0.35 C computed from the leading edge.

The eccentricity, e = 0.35 C − 0.25 C = 0.1 C, which is the distance among the elastic axis and aerodynamic center as in Fig. 10. Because the chord c varies or changes for different parts of the wing, e varies as well.

Fig. 10. Section of an airfoil illustrating the aerodynamic and shear centers.

Therefore,

(30) $$\begin{array}{r}{e}_i=0.1\times C_i\end{array}$$

That results as diagonal matrix as follows:

$$\begin{array}{r}[e]=\left[\begin{array}{c}0.2719000\\ 00.339200\\ 000.43990\\ 0000.5588\end{array}\right]\end{array}$$

Multhopp’s Quadratic Approach

Multhopp’s approximate quadrature is useful when dealing with functions derived from lifting line theory. utilizing the techniques and formula for a wing semi span l = b/2 to a symmetrical lift distribution problem yields the diagonal matrix of weighting values [6]. As follow:

(31) $$\begin{array}{r}\left[\overline{W}\right]=\frac{\pi I}{8}\left[\begin{array}{c}sin{\varnothing }_{1}000\\ 0sin{\varnothing }_{2}00\\ 00sin{\varnothing }_{3}0\\ 000\frac{1}{2}sin{\varnothing }_4\end{array}\right]\end{array}$$

Calculation Divergence Speed of the Wing with No Finite Span Correction Using Aerodynamic Strip Theory

The following results are obtained by multiplying the respective matrix values based on (28), which is:

$$\begin{array}{r}[E]=\left[C^{\theta \theta }\right]\left[e\right]\left[\overline{W}\right]\end{array}$$

$$\begin{array}{r}[E]=\left[\begin{array}{c}0.19490.19760.14070\\ 0.08570.19760.14070\\ 0.03600.08310.14070\\ 0000\end{array}\right]\times {10}^{-6}rad\frac{m}{N}\end{array}$$

We derive the following simplified expression from (27):

$$\begin{array}{r}\frac{1}{\frac{d{C}_l}{d\alpha }}\left[\frac{1}{C}\right]=q_d\left[E\right]\end{array}$$

When we rearrange the previous expression, we obtain:

(32) $$\begin{array}{r}\frac{1}{\frac{d{C}_l}{d\alpha }}=q_d\left[C\right]\left[E\right]\end{array}$$

The matrix product [C] [E] now becomes:

$$\begin{array}{r}\left[\begin{array}{c}0.53000.53730.38270\\ 0.29080.67030.47750\\ 0.15850.36540.61920\\ 0000\end{array}\right]\times {10}^{-6}\frac{\text{rad}\times m^2}{\text{N}}\end{array}$$

By using MATLAB environment to iterate this matrix numerically, we obtained matrix [C][E] = 1.3277$\times {10}^{-6}\text{rad}.{\text{m}}^2/\text{N}$ after nine iterations. We attained the divergence speed of $v_{Div}$ = 472.8420 $\text{m}/\text{s}$ by replacing the dynamic pressure $q_d$ b $1/2\rho v_{Div}^2$.

${\displaystyle \frac{d{C}_l}{d\alpha }}$= 5.5 and ρ = 1.225 kg/${\text{m}}^3$ (at sea level).

Calculation of Wing Divergence Speed with Finite Span Correction Using Aerodynamic Strip Theory

For the strip theory modified for finite span to be useful and accurately forecast the divergence speed of a practical wing, namely, a finite wing, the following calculations must be conducted: where: $\left(AR/AR+2\right)$ is the correction factor of finite span.

$$\begin{array}{r}\frac{d{C}_l}{d\alpha }=a_0\left(\frac{AR}{AR+2}\right)\end{array}$$

wing flat form aspect ratio: AR $={\displaystyle \frac{b}{\overline{c}}}$, b is the span of the wing and $\overline{C}$ is the standard mean chord or geometric.

As seen in the previous equation, the aspect ratio varies with wing span. As a result, the effective lift coefficient curve slope ${\displaystyle \frac{d{C}_l}{d\alpha }}$ will change when the finite wing is taken into account [5].

Let’s now calculate the aspect ratio AR = ${\displaystyle \frac{b}{{\overline{C}}^{.}}}$.

Because the chord changes over the span, we must determine $\overline{c}$. But $\overline{c}$= $S/b$, where S is the area of the wing.

Planform of a wing for the case study is as Fig. 11.

Fig. 11. Planform of a wing.

We know that $C_i=5.588(1-{\displaystyle \frac{y_i}{22.86}})$ from Eq. (29a).

Using $\overline{C}={\displaystyle \frac{S}{b}}={\displaystyle \frac{{\int }_{-b/2}^{b/2}{C}_i\left(y\right)dy}{{\int }_{-b/2}^{b/2}dy}}$

Because of symmetry $\overline{C}=S/b={\displaystyle \frac{2{\int }_0^{b/2}{C}_i\left(y\right)dy}{b}},$

where $b={\int }_{\frac{-b}{2}}^{\frac{b}{2}}dy=22.86\text{m}$

This means:

S $=2{\int }_0^{b/2}{C}_i\left(y\right)$dy

Therefore, $\overline{C}=\frac{S}{b}=\frac{102.5088}{22.86}=4.4842\text{m}.$

$$\begin{array}{r}AR=\frac{b}{\overline{C}}=\frac{22.86}{4.4842}=5.097\end{array}$$

As a result, the strip theory corrected for finite span: where $a_0=5.5$ denotes the lift slope.

$$\begin{array}{r}\frac{d{c}_l}{d\alpha }=a_0\left(\frac{AR}{AR+2}\right)=5.5\times \frac{5.097}{5.097+2}=3.95\end{array}$$

For the finite span correction, the divergence speed is:

The Influence of Increasing the Wing’s Stiffness Rigidity on the Divergence Speed

Assume that the GJ is increased by 15% uniformly along the span of the wing. The matrix $\left[C^{\theta \theta }\right]$ representing the torsional flexibility influence coefficients is as follows: produced matrix of [E] as follow:

$$\begin{array}{r}\left[C^{\theta \theta }\right]=\left[\begin{array}{c}32.6614.366.0370\\ 14.3614.366.0370\\ 6.0376.0376.0370\\ 0000\end{array}\right]\times {10}^{-8}\frac{rad}{N{m}^2}\end{array}$$

$$\begin{array}{r}[E]=\left[C^{\theta \theta }\right]\left[e\right]\left[\overline{W}\right]\end{array}$$

$$\begin{array}{r}=\left[\begin{array}{c}\mathrm{0.16950.17180.12240}\\ \mathrm{0.07450.17180.12240}\\ \mathrm{0.03130.07220.12240}\\ 0000\end{array}\right]\times {10}^{-6}\frac{rad.m^2}{N}\end{array}$$

The matrix product [C][E] now becomes:

$$[C][\text{E}]=\left[\begin{array}{c}0.46100.46720.33280\\ 0.25280.58280.41520\\ 0.13780.31770.53840\\ 0000\end{array}\right]\times {10}^{-6}{\displaystyle \frac{rad.m^2}{N}}$$

By using MATLAB code, after nine iterations we obtained the result of the matrix [C][E] to be 1.1544 $\times {10}^{-6}{\displaystyle \frac{rad.m^2}{N}}$. Replacing the dynamic pressure $q_{d}by1/2\rho v_{Div}^2,\rho =$1.225 kg/${\text{m}}^3$ (at sea level) and ${\displaystyle \frac{d_{cl}}{d\alpha }}$ = 5.5, we obtained the divergence speed of infinite wing to be: $v_{Div}=507.0839\text{m/s}$. Also, for finite span correction with ${\displaystyle \frac{d_{cl}}{d\alpha }}=3.95,$ the divergence speed obtains: $v_{Div}=598.3601\text{m/s}$.

Discussion of the Results

Aerodynamic strip theory has been used throughout this research to determine the divergence speed of unswept-tapered wing (both with and without finite span correction). Divergence speed has been found to increase by approximately 18% when the finite span correction is used as represented on Table. 1. As a result, the divergence speed produced by the strip theory analysis approach without finite span correction is lower than that of a three-dimensional wing, indicating that it is conservative.

Table I. Wing Divergence Speed of Flight at Sea Level

Method of analysis	Divergence speed, VD (m/sec)
Strip theory (No finite span correction)	472.8420
Strip theory (With finite span correction)	557.9546
Strip theory (With increment of 15% of GJ and without finite span correction)	507.0839
Strip theory (With increment of 15% of GJ and with finite span correction)	598.3601

To investigate the influence of increasing up torsional stiffness on divergence speed, a 15% increase in torsional stiffness is tentatively examined, and the divergence speed is predicted via strip theory as well (without and with finite span correction). When the torsional stiffness is raised by 15%, the divergence speed increases by around 18% as illustrated on Table. 1. This is due to the fact that a wing’s divergence speed is exactly related to its torsional stiffness.

Conclusion

To measure torsional deformation of wing under pure torsion, the torsional influence coefficient is utilized. Nevertheless, this matrix was created using the technique of matrix iteration to calculate the wing’s divergence speed. The unknowns in the governing integral equation are the elastic twist distribution and the lift coefficient, making the problem statically-indeterminate. Strip theory used to address this problem. The matrix technique is used to numerically solve the governing integral equation. The matrix has been iterated using MATLAB environment in order to converge on a solution. Strip theory has been used to calculate the wing’s divergence speed with and without finite span corrections. It has also been investigated how the wing’s torsional stiffness affects the speed of torsional divergence. The divergence speed predicted based on strip theory is found to be around 18 percent higher for (3-D wing) in comparison with (2-D wing). Based on strip theory, two-dimensional torsional divergence analysis offers conservative torsional divergence speed. A 15% increase in torsional stiffness caused around 18% improvement on divergence speed of a (3-D wing).

To sum up, this result demonstrates that the divergence speed of a wing is related to the square root of the torsional stiffness [6]. This supports the outcome obtained from the case study presented in reference [5]. The study’s findings are required for the construction of high-performance aircrafts and might also be used by designers for static aeroelastic analysis of aircraft wings.