Simplified and Accurate Stiffness of a Prismatic Anisotropic Thin-Walled Box

RESEARCH ARTICLE Simplified and Accurate Stiffness of a Prismatic Anisotropic ThinWalled Box Giacomo Canale, Felice Rubino, Paul M. Weaver, Roberto Citarella and Angelo Maligno Rolls-Royce plc, Moor Lane, Derby, UK Department of Industrial Engineering, University of Salerno, Giovanni Paolo II 132, 84084 Fisciano (SA), Italy Department of Aerospace Engineering, University of Bristol, Queens Building, University Walk, Bristol BS8 1TR, UK Institute for Innovation in Sustainable Engineering, University of Derby, Quaker Way, Derby DE1 3EE, UK


INTRODUCTION
Beam models have proven themselves to be effective in a variety of engineering structures.In aeroelastic studies, for example, beam models have been used to reveal important trends and generally, rectangular or trapeze cross sections have been used to model wing boxes [1 -4].
One example of prismatic thin-walled box of rectangular cross section is shown in Fig. (1).Its global coordinate system, shown in Fig. (2), is denoted with letters X, Y, Z, and it can be located in any point of the cross section.The origin of the frame is denoted by the letter O.For aeroelastic purposes, the structures analysed are often slender and long, so that three essential stiffness constants are required to evaluate the deformations.As suggested by Weisshaar [5], at any cross section along the reference x-axis of the beam, a linear equation describes the structural behaviour:  Once the internal loads M(X) and T(X) are calculated from the external distribution of forces and moments, Eq. 1 can be inverted and the deformation of the beam can be evaluated.
In order to obtain appropriate deformations, the evaluation of stiffness EI, GK and K must be as accurate as possible.A correct evaluation of the stiffness parameters of a composite beam is not a trivial problem as confirmed by the number of different formulations presented in literature [6 -15].Such models often show a lack of precision when cross sections with different geometries (even simple rectangular cross sections) and unbalanced lay-ups are analysed.This effect is more accentuated when a high level of anisotropy exists.It is the case of reinforcements made by fibres with the same orientation.
Concerning bend-twist coupling stiffness, K, the analytical predictions of Lemanski and Weaver [12] is relatively accurate, simple and based on physical reasoning.Their results have been investigated and reproduced as part of this work.Examples, in fact, are shown in section 3.No further discussion on the stiffness K is therefore provided.
The same cannot be said for bending and torsional stiffness, since formulae presented in the literature are not able to provide sufficient accuracy when cross sections of different geometries and lay-ups are analysed.Therefore, models able to predict accurate results are required.
The approach used by Lemanski and Weaver to calculate K is extended in this paper to evaluate the bending stiffness EI of a symmetric composite thin-walled box.A new analytical formula is also proposed to evaluate GK.

Evaluation of EI
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: 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.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 [10]:   N xy shear force per unit of length of horizontal laminate they are evaluated by using Classical Lamination Theory [16].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 Δ 1 and Δ 2 are added for the following reasons: Lateral deformation ε y of the horizontal laminate is constrained by the vertical walls and depends not only on the 1.
characteristics of the laminate itself but also on the elastic properties of the vertical walls.Shear deformation γ xy of the laminate is present; otherwise no bend/twist coupling effect exists.

2.
The strain field of the horizontal laminate is re-written as where A ij is the membrane constants of the top laminate.The terms Δ 1 and Δ 2 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   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 [12].The displacement, δ, is that of a cantilever beam [14]: where N y is found directly from Classical Lamination Theory, If Eq. 7 is substituted in to Eq. 6 and combined with Eq. 5, then: The second expression can be deduced from the equilibrium of tangential forces, as suggested in reference [10]: where G v is the shear modulus of the vertical wall.
An algebraic system comprising two equations (Eq.8 and 9) in two unknowns Δ 1 and Δ 2 is obtained and is readily solved.Once Δ 1 and Δ 2 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).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, Eq. 10 can be divided into two components, those from the horizontal and vertical walls.The contribution from the top and bottom laminates is: Where: 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: The contribution of the vertical walls to the bending stiffness is, consequently: where t v is the thickness of vertical walls.The final expression of bending stiffness EI is:

Evaluation of Torsional Stiffness GK
In this section, a model to predict the torsional stiffness GK is presented.The starting point is the formula developed by Librescu and Song [7].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 [10], 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: There are two contributions to the torsional stiffness: the contribution given by the vertical walls: the contribution given by the top and bottom laminates: 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 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).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 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 The stiffness GK can be therefore evaluated with Bredt's formula [17]: Where N is an average expression of the shear flow along all the contour found from 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.

RESULTS
Two different cross-sections have been analysed.The first one is with rectangular shape.Even if the geometry is simple, to the authors' knowledge the analytical model proposed in this paper is the only one providing good accuracy for structures with high level of anisotropy.The second one is a trapeze cross section, with realistic laminates.It has been analysed to show that the model is able to provide good results with more realistic structures.

Rectangular Cross Section
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: 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.Finite element analysis (FE) using Patran/Nastran [18] 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).The following steps were done to obtain bending and torsional stiffness.
A unitary bending moment is applied to the tip of the beam.1.
Bending and twisting deformations are measured in a cross section located at the middle span of the beam: 2.
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).
A unitary twisting moment is applied to the tip of the wing with no bending moment.

3.
Twisting deformation is measured and consequently twisting compliance.

4.
Once all the compliances are known, they can be inverted.The inverse of the compliance matrix is the stiffness 5. matrix of Eq.1.Results for EI of unbalanced boxes are shown in Figs.(11)(12)(13).The analytical model presented in this paper has been compared with two other different models (Kollar and Pluzsik theory [10] and Librescu and Song's theory [7] and also with FE.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.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)(16)(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.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 [12].The development of a new theory is therefore not needed.In Figs.(18)(19)(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.

Trapeze Cross Section
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.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.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

(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 H p
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.

DISCUSSION
With regards to the bending stiffness EI, the current analytical formulation gives the same level of accuracy as the model of Kollar and Pluzsik but its formulation is simpler.The current model does not evaluate the contribution of vertical and horizontal walls separately.The contribution of the top and bottom walls to EI is considered as a function of the stiffness of the vertical walls and vice versa.
With regards to the torsional stiffness, GK, the current analytical model is able to give accurate results for several geometries (rectangular and trapeze cross sections) and lay-ups (high anisotropy or realistic laminates) while the other models underestimate or overestimate the stiffness.The model of Librescu and Song has been enhanced in two steps.In the first one, a parabolic distribution of the maximum shear flow has been applied.In the second one, the average shear flow has been calculated in order to apply the Bred's formulation of thin walled closed cross sections.

CONCLUSION
An analytical model to evaluate bending and torsional stiffness of a composite box has been presented.It is relatively accurate, simple and based on physical reasoning.It can be also easily extended to trapeze cross sections.For cross sections with high level of anisotropy, the accuracy of the proposed formulation is within 2% for EI and within 6% for GK whilst other models can give errors of ca.40% and 100% respectively.These percentages of error of the proposed model have been approximately found constant for all the geometries analysed.The percentage of error is even further reduced when more realistic (and less anisotropic) laminates are used.Analyses performed on EI and GK show that the formulation proposed is able to give highly accurate results for different dimensions, length ratios of horizontal and vertical walls and different lay-ups.

K
the bend-twist coupling stiffness [Nm 2 ] M(X) the bending moment at the span wise coordinate x [Nm] T(X) the twisting moment at the span wise coordinate x [Nm] X the span wise coordinate along the main dimension of the beam [m]

1 ] 1 ]
Engineering Journal, 2018, Volume 12 3 ϑ(X) the twisting angle at the span wise coordinate x [m φ(X) the flexural angle and the span wise coordinate x [m - 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] strain νxy Poisson's ratio of composite laminate
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.
of the membrane matrix of the vertical wall Ω the area enclosed by the contour of the cross section

H 1
the height of the left wall H 1 the height of the right wall0

1 ]
of the vertical wall of the box [m] K = The bend-twist coupling stiffness [Nm 2 ] M(X) = The bending moment at the span wise coordinate x [Nm] at the span wise coordinate x [Nm] t h = Thickness of the horizontal wall of the box [m] t v = Thickness of the vertical wall of the box [m] w = Length of the horizontal wall [m] X = The span wise coordinate along the main dimension of the beam [m] at the span wise coordinate x [m and the span wise coordinate x [m -1 ] Ω = Area enclosed by the cross section [m 2 ]