The strategy used by Lemanski and Weaver to evaluate bend-twist coupling stiffness K is now extended to the evaluation of the bending stiffness EI.

The key point is the definition of the stiffness K: it is the twisting moment arising in a cross section when a unitary bending curvature is applied.

When a unitary bending curvature is applied, (without twisting, i.e. a torsional moment is obtained from Eq. 1 as:

(2)

The numerical value of the stiffness therefore coincides with the twisting moment arising in the walls of the cross section. Similar reasoning can be used to evaluate EI: as the bending moment arising in the cross section when the same unitary bending deformation is applied.

Consider a hollow box shown in Fig. (3).

The box is symmetrical, such that the top and bottom laminates are identical. Moreover, the vertical walls are made by balanced laminates, as unbalanced laminates are not necessary to induce bend-twist coupling of the box. Let us denote its geometrical characteristics by:

w the length of the horizontal wall [m]

2h the length of the vertical wall [m]

t_y the thickness of the vertical wall [m]

t_h the thickness of the horizontal wall [m]

Now apply a unitary bending deformation. On the top wall of the box, if the material is orthotropic and the stiffness of vertical walls negligible, then strains would be [10L.P. Kollar, and A. Pluzsik, "Analysis of thin-walled composite beams with arbitrary layup", J. Reinf. Plast. Compos., vol. 21, no. 16, pp. 1423-1465.
[http://dx.doi.org/10.1177/0731684402021016928] ]:

(3)

Where:

ε_x axial strain

ε_xy transversal strain

V_xy shear strain

νxy Poisson’s ratio of composite laminate

Fig. (3) Geometric characteristics of a box cross section.

The strains of Eq. 3 are written with respect to the local coordinates x, y of the laminate, shown in Fig. (4).

Fig. (4) Local coordinates in a laminate wall.

The forces per unit of length corresponding to these strains are:

N_x axial force per unit of length of horizontal laminate

N_y lateral forces per unit of length of horizontal laminate

N_xy shear force per unit of length of horizontal laminate

they are evaluated by using Classical Lamination Theory [16Z. Gurdal, R.T. Haftka, and P. Hajela, Design and optimization of Laminated composite materials., Wiley Interscience: New York, .]. As the box is made of laminated composite materials and the presence of vertical walls constrains the deformation of horizontal laminate, the strain field, Eq. 3, requires modification.

Two corrective terms Δ₁ and Δ₂ are added for the following reasons:

1. Lateral deformation ε_y of the horizontal laminate is constrained by the vertical walls and depends not only on the characteristics of the laminate itself but also on the elastic properties of the vertical walls.
2. Shear deformation γ_xy of the laminate is present; otherwise no bend/twist coupling effect exists.

The strain field of the horizontal laminate is re-written as

(4)

where A_ij is the membrane constants of the top laminate. The terms Δ₁ and Δ₂ are unknown. Two algebraic equations are needed to determine them. The first equation can be written evaluating displacement δ of Node 1, shown in Fig. (5).

Node 1 can be thought as a part of the horizontal laminate. Therefore, the displacement δ can be written as

(5)

Fig. (5) Displacement δ of Node 1.

Fig. (6) An idealization of the deflection of the vertical wall.

On the other hand, Node 1 is also part of the vertical wall. Its deflection can be calculated by modelling half of the vertical wall as a cantilevered beam (Fig. 6), as suggested by Lemanski and Weaver [12S.L. Lemanski, and P.M. Weaver, "Flap-torsion coupling in sandwich beams and filled box sections", Thin-walled Struct., vol. 43, no. 6, pp. 923-955.
[http://dx.doi.org/10.1016/j.tws.2004.12.004] ]. The displacement, δ, is that of a cantilever beam [14C. Zhang, and S. Wang, "Structure mechanical modelling of thin walled closed section composite beams, Part 1: Single cell cross section", Compos. Struct., vol. 113, pp. 12-22.
[http://dx.doi.org/10.1016/j.compstruct.2014.03.005] ]:

(6)

where N_y is found directly from Classical Lamination Theory,

(7)

If Eq. 7 is substituted in to Eq. 6 and combined with Eq. 5, then:

(8)

The second expression can be deduced from the equilibrium of tangential forces, as suggested in reference [10L.P. Kollar, and A. Pluzsik, "Analysis of thin-walled composite beams with arbitrary layup", J. Reinf. Plast. Compos., vol. 21, no. 16, pp. 1423-1465.
[http://dx.doi.org/10.1177/0731684402021016928] ]:

(9)

where G_v is the shear modulus of the vertical wall.

An algebraic system comprising two equations (Eq. 8 and 9) in two unknowns Δ₁ and Δ₂ is obtained and is readily solved. Once Δ₁ and Δ₂ are calculated, they are substituted in to Eq. 4 and the strain field of the top laminate is completely defined.

Concerning vertical walls, their local coordinates x_v, y_v are shown in Fig. (7).

Fig. (7) Local coordinates of vertical walls.

Consequently, x_v and x, are coincident and they have the same direction of the global axis X. Forces per unit of length of vertical walls are:

N_xv the axial force per unit of length.

N_vv the lateral force per unit of length

N_xvy is the shear force per unit of length

While corresponding strains are given by:

ε_xv the axial strain

ε_vv the lateral strain

γ_xvy the shear deformation

As EI is the bending moment arising in the cross section when a unitary bending curvature is applied, it can be calculated as the bending moment with respect to the global Y axis,

(10)

Where:

dC the infinitesimal element of the contour

the axial force per unit of length in the global reference. It includes contributions from both Nx and Nxv.

Eq. 10 can be divided into two components, those from the horizontal and vertical walls. The contribution from the top and bottom laminates is:

(11)

Where:

(12)

And the strains of Eq. 12 are found from Eq. 4.

Now consider half of the vertical wall, in order to evaluate the second component. The strain ε_xv, for simple bending, is a linear function of the global coordinate Z:

ε_xv = Z

(13)

The contribution of the vertical walls to the bending stiffness is, consequently:

(14)

where t_v is the thickness of vertical walls. The final expression of bending stiffness EI is:

(15)

In this section, a model to predict the torsional stiffness GK is presented. The starting point is the formula developed by Librescu and Song [7L. Librescu, and O. Song, "On the static aeroelastic tailoring of composite aircraft swept wings modelled as thin-walled beam structures", Compos. Eng., vol. 2, no. 5-7, pp. 497-512.
[http://dx.doi.org/10.1016/0961-9526(92)90039-9] ]. It has been chosen among several formulations because it shows two good characteristics:

• It is quite accurate, especially if compared to the other models investigated.
• Its formulation is relatively straightforward.

Other models, such as that due to Kollar and Pluzsik [10L.P. Kollar, and A. Pluzsik, "Analysis of thin-walled composite beams with arbitrary layup", J. Reinf. Plast. Compos., vol. 21, no. 16, pp. 1423-1465.
[http://dx.doi.org/10.1177/0731684402021016928] ], for example, contain more information and are more involved to implement.

The initial formula proposed by Librescu and Song can be re-arranged as follows:

(16)

Where:

A_ij^v the terms of the membrane matrix of the vertical wall

Ω the area enclosed by the contour of the cross section

(17)

There are two contributions to the torsional stiffness:

• the contribution given by the vertical walls:
• the contribution given by the top and bottom laminates:

(18)

(19)

Using Lemanski and Weaver’s approach, the stiffness GK can be thought as the twisting moment arising in a cross section when a unitary twisting curvature is applied. The forces per unit of length arising on the vertical and horizontal walls are derived from Eq.16, as

(20)

To ensure dimensional accuracy, forces per unit length in Eq. 20 should be divided by . These forces per unit of length are usually not equal and they are constant along the walls. This fact implies a discontinuity of the tangential forces per unit of length at the corners of the cross sections, as shown in the example of Fig. (8).

Fig. (8) Discontinuity of tangential forces in the model of Librescu and Song.

The shear continuity at the corners is imposed by assuming a parabolic distribution in the horizontal laminates, of force per unit of length, with the maximum value between N_xyv and N_xy. For example, when N_xy is greater than N_xyv, the shear flow of the top/bottom laminate will be assumed to vary parabolically: Its value at the corner will be equal to N_xyv, and its maximum value, at the middle point of the wall, will be equal to N_xy (Fig. 9). In other words, the shear flow of the horizontal laminates is now a parabolic function N_xy(y), where y is the local axis of the laminate, shown in Fig. (4). An analogous distribution occurs when N_xyv is greater than N_xy. In this case, the shear flow of the vertical laminates is assumed to vary parabolically. The average shear flow in the horizontal laminate is

(21)

The stiffness GK can be therefore evaluated with Bredt’s formula [17 Pitagora Ed. Erasmo Viola, Theory of Constructions, Theory of Elasticity, Vol. 1 , 1990, p. 628. ISBN: 978-8837104931]:

(22)

Where N is an average expression of the shear flow along all the contour found from

(23)

This formulation provides excellent results for rectangular geometries and it can be easily extended to the trapeze cross sections, as it will be shown in the section 3.2.

Fig. (9) Parabolic correction of Nxy (in red) on the Librescu formulation (in black).

Three rectangular cross sections representing three different boxes, have been analyzed. Top and bottom laminates are made with one single layer whose fibres orientation α can vary from 0 to 90 degrees with respect to the local frame represented in Fig. (4). This assumption has been made to maximise the effects of anisotropy. Vertical walls are orthotropic. They are made with one single layer with fibres oriented at 0 degrees. Appropriate elastic properties for both vertical and horizontal laminates are:

E₁ = 181 GPa

E₂ = 10.3 GPa

G₁₂ = 4.55 GPa

ν₁₂ = 0.28

Elastic properties of α- oriented laminates are calculated by using the constants reported above and the classical lamination theory.

These three cross sections have different wall lengths, but the same shape and the same area enclosed by the contour. Data are reported in Table 1.

Table 1
Geometric properties of three different boxes.

Finite element analysis (FE) using Patran/Nastran [18 MSC software, MSC Nastran 2016. Newport Beach, CA, 2016. Available from: http://documents.mscsoftware.com/ sites/default/files/ro_nastran-2016_ltr_w.pdf] was performed to validate the results. The displacements and rotations are set to zero on one side of the beam whilst MPC of type RBE2 have been used to apply the end tip loads as shown in Fig. (10).

Fig. (10) An example of FE model used to validate the analytical results.

The following steps were done to obtain bending and torsional stiffness.

1. A unitary bending moment is applied to the tip of the beam.
2. Bending and twisting deformations are measured in a cross section located at the middle span of the beam: Sufficiently far from the tip and root, in order to avoid the effects of local deformations. The ensuing bending and twisting deformations, due to a unitary bending moment, represent bending and bend/twist coupling compliances (Eq. 1).
3. A unitary twisting moment is applied to the tip of the wing with no bending moment.
4. Twisting deformation is measured and consequently twisting compliance.
5. Once all the compliances are known, they can be inverted. The inverse of the compliance matrix is the stiffness matrix of Eq.1.

Fig. (11) EI for the “Tall” unbalanced composite wing box. Fibres angle vary from 0 to 90 degrees.

Results for EI of unbalanced boxes are shown in Figs. (11-13). The analytical model presented in this paper has been compared with two other different models (Kollar and Pluzsik theory [10L.P. Kollar, and A. Pluzsik, "Analysis of thin-walled composite beams with arbitrary layup", J. Reinf. Plast. Compos., vol. 21, no. 16, pp. 1423-1465.
[http://dx.doi.org/10.1177/0731684402021016928] ] and Librescu and Song’s theory [7L. Librescu, and O. Song, "On the static aeroelastic tailoring of composite aircraft swept wings modelled as thin-walled beam structures", Compos. Eng., vol. 2, no. 5-7, pp. 497-512.
[http://dx.doi.org/10.1016/0961-9526(92)90039-9] ] and also with FE.

Fig. (12) EI for the “Square” unbalanced composite wing box. Fibres angle vary from 0 to 90 degrees.

Fig. (13) EI for the “Representative” unbalanced composite wing box. Fibres angle vary from 0 to 90 degrees.

The analytical model is shown to be accurate and it gives approximately the same results as Kollar and Pluzsik theory, but its formulation is simpler.

The model has also been tested with structures whose top and bottom walls are made of balanced composite materials (laminates with + α and -α fibres). Results obtained for the square wing box are shown in Fig. (14). In this case, all the models give approximately the same results.

Fig. (14) EI for the “Square” balanced composite wing box. Fibres angle vary from 0 to 90 degrees.

Results for the evaluation of GK have also been produced. The new formulation has been compared with the models of Librescu and Song and Kollar and Pluzsik, respectively Figs. (15-17). It is evident that the model proposed in this paper is the only one (to our knowledge) which works sufficiently well for all three different geometries.

Fig. (15) GK for the “representative” unbalanced composite wing box. Fibres angle vary from 0 to 90 degrees.

Fig. (16) GK for the “square” unbalanced composite wing box. Fibres angle vary from 0 to 90 degrees.

Fig. (17) GK for the “square” unbalanced composite wing box. Fibres angle vary from 0 to 90 degrees.

The main reason lies on the fact that the proposed way to calculate stiffness is based on its physical meaning (deformation under unitary loads) rather than the integration of stiffness along the segments of the cross section (furthermore Librescu’s models offer a slightly simpler formulation compared to Kollar’s model and the results could be slightly less accurate when dealing with high level of anisotropy).

Concerning bend/twist coupling stiffness K, good results have been obtained by using Lemanski-Weaver's model [12S.L. Lemanski, and P.M. Weaver, "Flap-torsion coupling in sandwich beams and filled box sections", Thin-walled Struct., vol. 43, no. 6, pp. 923-955.
[http://dx.doi.org/10.1016/j.tws.2004.12.004] ]. The development of a new theory is therefore not needed. In Figs. (18-20), a study case on boxes described in Table 1 is shown. The models of Lemanski-Weaver and Kollar-Pluzsik have the same level of accuracy, but the implementation of the model of Lemanski and Weaver, also in this case, is simpler.

Fig. (18) Bend/twist coupling stiffness K for a representative wing box.

Fig. (19) Bend twist coupling stiffness K for a square wing box.

Fig. (20) Bend twist coupling stiffness K for the “tall” wing box.

A more complex geometry is studied in this section. Rear and front webs have different lengths. The cross section, in other words, has a trapeze shape (Fig. 21). This geometry can be used as a simplified model of wing boxes. Their stiffness can be used for aeroelastic applications. It is also important to remark that such stiffness parameters are independent from the angle of sweep of the wing.

Fig. (21) Trapeze cross section: geometrical properties.

Top and bottom walls are built with laminates having the following volume fractions: 40% 0 degrees fibres, 20% 90 degrees fibres and 40% α degrees fibres, with α varying between 0 and 90 degrees. Vertical walls are orthotropic. They are made by 20% 90 degrees fibres, 40% 0 degrees fibres and 40% 45 degrees fibres. Appropriate elastic properties are those reported in section 3.1.

Geometrical properties are described in Table 2.

Table 2
Geometrical properties of the trapeze box.

The analytical model for EI described in section 1 can be applied once the structural symmetry of the cross section is found, even if the structure in not geometrically symmetric. In other words, the analytical model can be applied once the neutral axis (only for bending along the Y-axis) is found, as shown in Fig. (22). The neutral axis is such that

Fig. (22) Neutral axis of a trapeze cross section.

(24)

Where:

Z_n the distance of the point of the cross section from the neutral axis.

The position of this neutral axis varies slightly with the angle α. However, when the difference between the lengths of left and right walls is not large, and it is often the case of wing box models, the position of the neutral axis with the respect to the axis shown in Fig. 22 can be calculated by using the simplified formula

(25)

Where:

p the height of the neutral axis (Fig. 21).

H₁ the height of the left wall

H₁ the height of the right wall

Bending stiffness is calculated by applying the analytical formulation of section 1 to the semi-rectangular cross sections placed below the neutral axis.

Models to evaluate GK and K can be applied to trapeze cross sections without difficulties. Numerical results are shown in Figs. (23, 24 and 25). The analytical model has been compared, also in this case, with the models of Kollar and Librescu and with FE results, obtained with the same technique described in section 3.1. Concerning EI, all the models produce almost the same result. Regarding K and GK, the analytical model is the only one with a good level of accuracy. It is important to remark that the model of Librescu underestimates the torsional stiffness, even if no parabolic correction is applied. The model of Canale and Weaver, on the other hand, predicts good results because of the assumption resumed in Eq.23.

Fig. (23) Bending stiffness of a trapeze cross section with realistic laminates.

Fig. (24) Torsional stiffness of a trapeze cross section with realistic laminates.

Fig. (25) Bend-twist coupling stiffness of a trapeze cross section with realistic laminates.