Linear and non-linear finite element analysis of shear-corrected composites box beams

The Updated Lagrangian formulation for non-linear finite element analysis is applied to the problem of thin-walled composites box beams undergoing large displacements. A shear correction factor for thin-closed rectangular sections is introduced into some terms of the variational formulation and its influence in the results is analyzed, in both linear and non-linear problems. The Vlasov’s theory describing the coupled flexural-torsional phenomenon and the implementation of the FSDT theory in thin-walled beams is discussed, as well as the strains, stress resultants and constitutive relationships for composites box beams. The application of the Mechanic of Laminated Beam theory (MLB) to the calculation of the shear correction factors considering more constitutive terms for computing the shear flow in flanges and webs than those ones used in the original theory is also debated. The Updated Lagrangian finite element model applied in this work for non-linear analysis of box beams is described. A comparison of the numerical results with those obtained experimentally, analytically and by the numerical models proposed by other authors is done for both linear and non-linear problems. The assumptions made in this work and the formulation developed only applies to thin-walled beams that undergo large displacements, but small to moderate twist rotations.


INTRODUCTION
The thin-walled box beams are one of the most common industrial applications of composites materials.They can be found currently in the aerospace, naval and construction industries, principally like elements supporting flexural and torsional loads.Most of these elements are produced by pultrusion.
In many practical situations, flexural-torsional coupling and nonlinearities arise in those beams.The finite element method has emerged as an efficient computational numerical technique to analyze that kind of problems.The theoretical foundations of most of the FEA methods for this kind of analysis have been extracted from the theories of Vlasov [1] and Gjelsvik [2].Hodges et al [3] applied the Variational Asymptotic Method (VAM) to non-linear three dimensional analysis of beams, obtaining the well-known VABS method (Variational Asymptotic Beam Section analysis), that reduces the three dimensional problem to a one-dimensional cross-sectional one, leading to the Vlasov beam element shown in the Figure 2 and used in this work.
In order to simplify the analysis of thin-walled beams, many authors have modeled the plate elements of the beams by separately using the CLPT theory [4,5], that can be done if transverse shear effects can be neglected.However, as it was shown by Gu and Chattopadhyay [6], the transverse shear strains have relevant effects in the case box beams made of anisotropic composite materials, even for very thin-walled beams.In the present work, both transverse and warping shear strains are considered in the FEM formulation using the model proposed by Vo and Lee [7], but the warping shear strain energy is neglected in the shear correction factors computations based on the premise of small to moderate twist rotations.
Ones of the most valuable researches about the analytical and numerical treatment of the flexural-torsional behavior of thin-walled beams have been conducted by Jaehong Lee.A general theoretical model based on the FSDT theory to describe the behavior of composites box beams under flexural and torsional loads was developed by this author and his co-workers, for both linear and non-linear problems [7,8].The present work makes use of the aforementioned theoretical model, so it is concisely presented in the Part 2 of the paper.One substantial difference between the Lee's model and the model presented here is the introduction of the shear correction factor into the variational formulation for some special cases of laminates, as it is explained in Section 4. The shearcorrected linearized equation of motion in Updated Lagrangian formulation is presented in the Section 3.
As it was mentioned above, shear correction factors will be introduced into the variational formulation of the problem treated in this work for some special cases of plating.One of the most used methods for the shear correction factor determination comes from the theory of Mechanics of Laminated Beams (MLB) constructed by Barbero et al [9].This theory has been employed by many authors to compute shear correction factors in analysis of pultruded composites beams [10,11] and it will be also used in this work to calculate the shear correction factors added in the stiffness matrix of the finite element model.

THE VLASOV WARPING THEORY AND THE MODEL OF LEE
The axial displacement of the cross section due to the torsion is called warping and it causes the deplanation of the transverse planes.When a torsion-loaded thin-walled beam does not warp because of the form of the cross-sections (circular, X-shaped, T-shaped, L-shaped, among others) or when the warp displacement is not constrained in this element, the classical Saint Venant torsion theory describes well the torsion phenomenon.In that case, it is considered that all torsional torque originates only shear stresses [12].However, when the cross section is restrained against axial displacement at some locations along the member or when a torsion moment loading is imposed in any position of the cross section, it appears a warping moment and the normal stresses in the longitudinal direction arise [13].In that case, Saint Venant theory does not apply anymore and Vlasov warping theory emerges to describe this complex stress state [12,13,14].
Latin American Journal of Solids and Structures 10(2013) 647 -673 In the Vlasov warping theory, for a free warping section, the axial displacement is compute by: u !z = −ω s .ϕ !(z), where "z" represents the longitudinal axis of the beam and "s" is called the contour coordinate of the section.For a restrained warping condition, normal stresses appear and they can be expressed by: σ !!,! s, z = −E z .ω s .ϕ !! z [14].In the last expressions, ϕ is the twist angle and ω(s) is called the warping function.
To a better understanding of the Vlasov warping theory, two coordinate systems are defined (Figure 1): The Cartesian coordinate system (x, y, z), and the local plate coordinate system (n, s, z).There is a point, P, with respect to which all global displacements and rotations are taken; this point is called the pole, shear center or center of twists and the axis passing through it and parallel to the "z" axis is called the pole axis; the angle of rotation of the cross section about that pole axis is the twist angle (ϕ).Both coordinate systems are related between them by an angle θ.
Taking into account the Figure 1, the general expression for the warping function, ω(s), can be written as [15]: Where F(s) is the Saint-Venant circuit shear flow and t(s) is the thickness of each point of the cross section.In the particular case of a box beams, the explicit expressions for the Saint-Venant circuit shear flow, F(s), and the warping functions, ω(s), can be found in [16].
According to Vlasov theory, when the warping is restrained, a warping moment or bi-moment is originated in the cross-section.This moment is considered in the present work and its equation and relationship with primary Saint Venant torsional moment can be found in [14,17].
The present formulation is based on the one-dimensional Vlasov beam element [8,[18][19][20][21], in which seven degrees of freedom are considered per master node (Figure 2), namely one axial dis- placement, two transverse displacements, one torsional rotation, two flexural rotations and one warping rotation.The last rotation is a measure of the warping intensity and its rate, of the magnitude of the warping moment [22].In the Vlasov beam element, transverse shear strains γ !" and γ !" and warping shear strain γ ! are included and it is assumed that they are uniform in each cross section of the beam element.According to this consideration and to the assumptions of the FSDT theory, Vo and Lee [7,8,23] obtained the following relationships (See Figure 1 for reference): • Displacement field of the contour midline in terms of displacements of the pole point P: • Displacement field of any generic point C in the section in terms of displacements of the contour midline (FSDT theory): Where"n !" is the distance along the normal direction between the midline contour and the generic point C and ψ !* (s, z) and ψ !* (s, z) denote the rotations of transverse plane about the "z" and "s" axis of the local plate coordinate system, respectively.The expressions for those rotations can be found in [7,8].

STRAINS, STRESS RESULTANTS AND CONSTITUVE RELATIONSHIPS
In the UL formulation, the work-conjugated pair used for the variational statement is the Almansi strain-Cauchy Stress tensor.In agreement with many authors that have analysed large deformation in Timoshenko`s beams [7,8,24,25], only nonlinear terms product of the variation in the longitudinal direction of transverse displacements will be considered in this work (Von Karman strains), leading to the following simplification of the Almansi strain tensor for any generic C point on the contour: A more complex analysis of composites beams by FSDT considering the whole terms of the strain tensor can be found in [26], where the shear locking effects are reduced significantly, but the computational cost is high.In this case, for sake of numerical simplicity, shear locking is dealt with the application of a reduced integration scheme in the components of the linear stiffness matrix containing A !! , A !" and A !! , which are shear extensional stiffnesses implicit in some coefficients of the constitutive matrix of the composites box beam.
Following the same procedure employed by Machado and Cortínez [27], the following expressions are gotten for the Von Karman strains of a generic C point in terms of the strains and curvatures of the midsurface of the wall: Where: are the biaxial curvatures in the midsurface.
is the high order biaxial curvature in midsurface of the wall.
Both the explicit expressions for the strains and curvatures in the midsurface in terms of the seven degrees of freedom of a Vlasov beam element (Figure 2) and the expressions for the stress resultant for flexural-torsional analysis of thin-walled beams can be found in [8] and those expressions are fitted to the UL formulation in the present work.

Latin American Journal of Solids and Structures 10(2013) 647 -673
Based on the constitutive in-plane and out-plane relationships for a lamina [28], direct constitutive equations relating the stress resultants and the global strains and curvatures of the cross section can be obtained, as it is detailed in [29].In general, these expressions can be represented in UL formulation as: Where: The explicit values for E !" can be gotten from [8,29].

SHEAR CORRECTION FACTOR
Like in most of the works about thin-walled composite beams considering shear deformability [7,22,30], here it is assumed that shear transverse strains (γ !" ,γ !" ) and warping shear strains (γ ! ) are uniform over the entire cross sections of the beam.This assumption leads to uniform shear strain energy in the global coordinate system (x, y, z), which is not necessarily true in most of cases.
To account for the error involved in this assumption, they are introduced shear correction factors into the constitutive relationships of shear stress and shear strains [28].
There are several methods to determine the shear correction factors in cross-section composites beams [31][32][33][34][35][36].In this work, an energy-based method coming from the MLB theory of Barbero et al [9] is employed and it is briefly exposed in the next lines.
Considering small to moderate twist angles, the contribution of the warping shear strain energy to the average shear strain energy can be neglected according to [37], leading to the next expression for the average shear strain energy per unit length in a cross section of the beam (superscripts and subscripts indicating the configuration are suppressed for sake of simplicity): Following the same procedure used by Barbero et al [9] and dividing the energy contribution into webs energy and flanges energy, the shear strain energy due to the actual shear strain distribution can be represented by the next equations: Latin American Journal of Solids and Structures 10(2013) 647 -673 The expressions for local laminate stress resultants appearing in equation ( 8) are: Local laminate load in the contour direction (9a) Local laminate load in the normal direction (9b) Where "i" represents both the flanges and the webs.Now, considering thin-walled hollow beams, they are neglected the shear deformations in the flanges and webs caused by the local laminate load in the normal direction, V !",! * , [9] what is equivalent to say that the shear stresses σ !",! are not considered in this case for the computing of the shear correction factor.According to these assumptions, the equilibrium equations for any point in the cross section are: Applying the constitutive relationships of the CLPT theory to each wall, neglecting the warping effects for the shear correction factors calculation, based on the supposition of small to moderate twist rotation, and converting the midsurface shell strains and curvatures into global beam strains and curvatures [38], the following expression for the local laminate axial load is gotten: Latin American Journal of Solids and Structures 10(2013) 647 -673 The equations ( 6) and ( 11) are used to obtain an explicit expression for N !* in terms of the resultants N !, M !, M !, M !, V ! and V !(bimoment, secondary torsional moment and high order normal force are not considered).The shear flow variation in the "i" wall is then obtained by using equation (10a) and from the equilibrium equations for a beam subjected to flexion and Saint Venant torsion.The final expressions are: • For the flanges: • For the webs: Where "w !(z)" and "w !(z)" are the distributed loads and "e !" and "e !" are the eccentricities of those loads.
The explicit expressions for the constants of the equations (12a) and (12b) can be found in Appendix A, where "α" is the shear-corrected constitutive matrix (See equation 16).
The integration of equations (12a) and (12b) to obtain expressions for the shear flow in the flanges and webs is done as it is described in [39] for closed sections.
According to [28], the shear-extension and bending-twisting coupling can be suppressed for laminates with off-axis plies balanced symmetrically (antisymmetric angle-ply laminate, for example).If the laminate is symmetric, not necessarily balanced, the coupling matrix, B, is equal to zero.In this work, only the former case will be considered for the shear correction factor calculations; the second case will be considered in further works.So, neglecting high order biaxial curvatures and considering all assumptions mentioned so far, the constitutive relationships of equations ( 9) are reduced to: As the layers of the laminate are thin enough, the influence of the coupling term B !" is small enough too, and in this case, the biaxial curvature κ !" (!) = −2 ∂ !u ∂s ∂z is neglected for the calculation of the shear correction factors under the hypothesis of small to moderate twist rotation.So, the equations ( 8) are significantly reduced to: Latin American Journal of Solids and Structures 10(2013) 647 -673 The shear correction factors are introduced into the constitutive relationships for V ! and V !, as it is indicated in the next equations (sum convention): The last equations lead to a modified constitutive matrix, α, defined as: The final expressions to find the shear correction factors are gotten by equating the average shear strain energy to the actual shear strain energy: Equation for Ksx: Equation for Ksy: Latin American Journal of Solids and Structures 10(2013) 647 -673 In linear analysis, the resultants appearing in Equations ( 17) are known in advance from the equilibrium equations in each section.In non-linear analysis, they are taking the resultants from the last configuration.
In the numerical examples presented here, they are reported the averaged values of the shear correction factors computed for the cross sections of the beam considered in the finite element analysis.

FINITE ELEMENT FORMULATION
This analysis is limited to box beams with equal-thickness flanges and equal-thickness webs, in such a way that the pole point coordinate ( !,  ! ) is located in the centroid of the beam.The basic steps of the continuum mechanics incremental decomposition for UL formulation as stated in [4] are applied to this analysis and it is obtained the following linearized equation of motion, where the shear correction factor as calculated on the Section 4 is introduced into the constitutive relationships of the shear stress resultants,  ! and  ! .
Latin American Journal of Solids and Structures 10(2013) 647 -673 From last equation, explicit expressions for the linear stiffness matrix can be stated.The non-linear stiffness matrix components are established in terms of the resultant loads in the current configuration.
Considering only the transverse loads and the torsional moments, the vector of external forces for an element has the next form: For the distributed loads, q (distributed transverse load) and t (distributed torsional moment), the load assigned to each node, "i", is defined as: .
For punctual loads, Q (punctual transverse load) and T (punctual torsional load), the load assigned to each node, "i", is defined as: T , where ξ ! is the point of application of the punctual loads, Qand T, in the isoparametric coordinate system.
The final non-linear finite element formulation is represented by: R( ∆ Where:
6 RESULTS AND DISCUSSIONS

Result comparison with analytical models
Both Euler-Bernoulli and Timoshenko beam models are used to estimate the vertical displacements in several locations of a cantilever composite box beams supporting a central distributed load.The whole sides of the box beam have the same laminate, which is done of chopped strand mats of E fiberglass with vinylester resin, having the following mechanical properties [42]:E != E ! = 14236 MPa, ν !" = 0.3287 and G !" = G !" = G !" = 5357 MPa.Some comparisons between these analytical models and the present numerical model are shown in Figure 4.As it can be appreciated, the present numerical solution is closer to the analytical solution of Timoshenko when the beam is shorter, and this matching is improved by the introduction of the shear correction factor, K !" , into the variational formulation.Despite the introduction of a reduced integration scheme into some terms of the stiffness matrix, the shear locking is yet appreciated, principally when beam is large, where stiffening effects in the results are more notorious.In this particular problem, the shear correction factor is K !",!"#$#%& = 0.4712 for L=0.25 m; the difference with the one computed by Bank [33] is 12.27 % (K !",!"#$ = 0.4197) and 8.50% is the difference with the factor calculated by Cowper [32] (K !",!"#$%& = 0.4355).

Results comparison with other models for nonlinear analysis of composites box beam
A non-linear problem proposed in [7] and [43] will be considered for comparison purposes.The problem consists in a cantilever composite box beam with the next data (see Table 1) The updated Lagrangian finite element model exposed in Section 5 is used, but the shear correction factor calculated as proposed here cannot be used in this example because the off-axis plies in the stacking sequence are not balanced.The number of iterations for each load step is determined in function of a minimum allowable value of E(δΔ) = 10 !! , as it is showed in Figure 5, where the error is plotted as function of the iterations to obtain the equilibrium configuration, considering θ = 45° in the flanges.It can appreciate that this particular problem is very stable, because the increase in the number of iteration leads to a reduction in the error, which is not ever the case.For θ = 45°and θ = 90°, the results obtained from various authors for the vertical displacement, "v", the horizontal displacement, "u", and the twist angle, "∅", the tip of the cantilever beam, are showed in Table 2.   [44] 48.50 cm 53.40 cm Vo and Lee (FSDT) [7] 55.00 cm 56.5 cm Vo and Lee (CLPT) [43] 51.3 cm 52.9 cm Present 58.8 cm 59.12 cm Horizontal displacement, u.

Reference θ=45° θ=90°
Stemple and Lee [44] 0.09 rad's 0 rad's Vo and Lee (FSDT) [7] 0.072 rad's 0 rad's Vo and Lee (CLPT) [43] 0.065 rad's 0 rad's Present 0.045 rad's 0 rad's The Table 2 indicates that the present formulation predicts larger vertical displacements of the box beam than other well-known models developed up to now.As in the models developed by Vo and Lee [7,43], this UL model predicts very closed values of the vertical displacement of the tip of the beam for stacking sequences in the flange with θ = 45° and θ = 90°, from which it can be inferred a weak influence of the angle θ in the vertical displacement of the beam in this range, what can be confirmed in the mentioned works [7,43].Contrary to the vertical displacement, the present UL model is stiffer than the other ones considered here for the horizontal displacement and the twist angle.
In Figure 6, it is showed the change of the three degrees of freedom considered in this example (v, u, ϕ) in the longitudinal direction of the beam for θ = 45°, solving the problem linearly and nonlinearly using the UL formulation.The differences among the results of these approaches are significant, what means that this problem is geometrically nonlinear and cannot be treated as a linear one by analytical or numerical methods.
Latin American Journal of Solids and Structures 10(2013) 647 -673 Figure 6 Differences between linear and UL non-linear solutions for v,u and ∅.

Result comparison with experimental data
For verification purposes, a cantilever short composite box beam was subjected to several loads in the tip.The whole experimental setup, which allows both applying instantaneous load and measuring directly the respective deflection, is shown in Figure 7.The setup consists of a 2-ton bridge crane, a sample holder, an electrical hoist, a 1-ton capacity load sling and open drum with lifting accessories.The upper part of setup also includes a DSLR camera at the same level of the composite box beam behind which is installed a measurement pattern to record (by a picture) the deflection at the moment when each load is applied.
The sample was prepared with a piece of steel angle at the tip to avoid load-sling sliding.In the other side a solid aluminum block was tightly inserted into the interior of the box beam.At the moment when sample is installed and screwed firmly, this internal insert ensures that fixed end condition and original squared geometry remains.The sample before load is shown in Figure 8.The camera focuses the measurement tapes, which were located in two directions in the test: the horizontal one to check the length where the load was applied and the longitudinal location of measurement stations and the vertical ones installed each 2 inches, to measure the vertical displacement directly when the picture was taken.
Several punctual loads were applied at the tip of the box beam with the objective to evaluate the vertical displacements at determined locations.The first step was calculating the drum´s tare and takes it into account in the mass.Then, at the floor level without the load sling connected to the drum, the recipient was filled until a predetermined level, taking care to measure exactly the quantity of liquid added; this mass of water plus the drum's tare and load sling´s weight are summed to calculate the first load that was applied to the tip of the beam.The mass was elevated using the electric hoist, then the load sling was connected to the drum and the load was applied obtaining a specific vertical displacement.
After the vertical displacement is measured digitally in the picture, the box beam is completely unloaded using the electric hoist.More water is added and the load process in the box beam is repeated several times, like it was explained above.The properties of the material of the box beam were obtained from manufacturer's information and they are exposed in Table 3.In the Figure 10, they are plotted the vertical displacements of the short beam obtained from the present model, considering the not shear and shear corrected (K !" = 0.415) variational formulation.It can be appreciated that the shear-corrected model is less stiff than the other one.In the Table 4 and Table 5, these displacements are compared with those obtained from the test described above.The relative error between numerical and experimental data can be attributed to the sum of three kinds of errors; numerical errors due to shear locking effects, experimental errors coming from the approximation of measurements and the lack of parallelism in the camera (despite special careful was taken in the visual alignment, no parallelism verification devices at milimetric scale employed) and errors in the characterization of the properties of the material coming from the manufacturer.The average relative errors point to conclude that the addition of the shear correction factor into the variational formulation originates numerical solutions better approximate to experimental data for the case of linear analysis, which is the one considered in this section.

Case study
The geometry, boundary conditions and lamina material of the composite box beam employed in Section 5.2 will be used in this case study and they are taken into account the following stacking sequences for both flanges and webs: [(θ/-θ)2]s (balanced and symmetrical), [(θ/-θ)4] (balanced and non-symmetrical) and [0/θ2/0]s (unbalanced and symmetrical).The torsion and warping rotations are analyzed as function of the orientation of the fibers and it is also examine the influence of the shear correction factor in the angle of twist.
In Figure 11, it is showed the variation of the twist angle for several values of the angle θ, considering a distributed load of q = −1000 N/m.As it can be seen, for most of fiber orientations, θ, a torsional coupling appears for the balanced and symmetrical laminate and the angles of twist predicted can be initially attributed to geometric nonlinearities and to the twisting-flexural coupling, D !" .In the case of the angle-ply antisymmetrical laminate, the results obtained are only slightly larger than the ones of the balanced and symmetrical plating.This means that in this particular case study, the material coupling terms D16 (present in the balanced and symmetrical laminate) and B16 (present in the angle-ply antisymmetrical laminate), have not a strong influence at inducing twist angles; therefore, in the first two stacking sequence analyzed, [(θ/-θ)2]sand [(θ/-θ)4], these twist angles can be essentially attributed to a bending-twisting coupling coming from the high order coupling terms, F !" , originated from geometric nonlinearities [8,29].The most important coupling torsion effects can be seen when plies are oriented in the ±15° directions, which is consistent with results of [7].For the symmetrical non-balanced stacking sequence considered in this case study, a stronger twisting-flexural and shear-extensional coupling can be observed, if the results are compared with the former two cases mentioned, being the maximum values of twist angles reached again in the ±15° directions, corroborating again the results of [7].In this case, all terms of coupling matrix, B, are zero, but the shear-extensional (A16) and twisting-flexural (D !" ) coupling terms still remains, and when compare Figure 11c with Figures 11a and 11b, it can be appreciated that the non-balancing of off-axis plies in the laminate brings to a significant influence of those material coupling terms in the twist angle.
In Figure 11, it can be observed that the change of the angle of twist is not constant for some fiber θ orientations.It is an indication of the presence of variable warping effects, and this can be confirmed in the Figure 12, where the variation of the angle of warping in the longitudinal direction of the beam is plotted, for the unbalanced symmetrical plating [0/θ2/0]s.For θ = 0°and θ = 90°, no significant warping angles are predicted by the model, what is consistent with the small values of the twist angles predicted in those fiber orientations.However, for other fiber orientations, warping angles are important.The longitudinal location of the point supporting the maximum warping effects is variable with the angle of orientation of fibers and the most relevant warping effects are presented when fibers are oriented at θ = 15°.From the point of view of the induced warping effects, the behavior of the plot indicates that the most critical section for the cantilever beam is among the fixed support, where warping moment is the maximum, and approximately 0.55 meters from that support, where warping deformation is the highest.
In order to examine the influence of the shear correction factor, K !" , into the results of the present UL formulation, the twist angles computed for the beam not shear and shear-corrected are plotted in Figure 13 for the balanced non-symmetrical stacking sequences considered in this case study, [(θ/-θ)4].In this figure, it is clearly appreciated that the introduction of K !" into the variational formulation leads to a very unstable solution for the twist angles; this solution also overpredicts these twist angles.It can be concluded that the introduction of the shear correction factor as proposed here does not give reliable solutions for the induced-torsion phenomenon in this UL non-linear analysis.This can be explained studying the assumptions implicit in the shear correction factors calculations.The material decoupling assumptions are related to the off-axis balanced symmetrically stacking sequence and to the small thickness of each lamina of the laminate, what is complied in this particular case.Another important assumption is the non-consideration of the shear stresses in the normal direction of the local plate coordinate system; this assumption is valid for walls with small thickness, what is complied in this case too.The other important assumption is to neglect the warping and high order shear strain energies; the small contributions of these energies to the total shear strain energy cannot be guaranteed in this case and the error associated to that assumption is probably the main cause of the instability of the shear-corrected solution for the twist angles, what suggests that those shear strain energies shall be considered in the shear correction factor computation for this kind of non-linear problems.

CONCLUDING REMARKS AND FUTURE WORKS
Linear and UL non-linear analysis of thin-walled composites box beams were carried out in the present work.A simple method for computing the shear correction factors of this kind of beams based on the MLB theory was proposed too.In the case of linear problems, the introduction of a shear correction factor into the variational formulation leads to more accurate results when compared them with analytical and experimental ones, essentially for short beams.In the case of non-linear analysis, the UL formulation developed shows to be in agreement with other well-known numerical models; however, it is important to mention that the present formulation is less stiff than those models for the vertical displacements, on the contrary to what happen with the horizontal displacements and the angle of twist.The introduction of the shear correction factor into the linear terms of the UL non-linear variational formulation brings to unstable solutions for the angle of twist, what can be attributed to the errors caused by the non-consideration of warping and high order shear strain energies into the calculation of the shear correction factors; these strain energies can have an important contribution in some problems of composites box beams undergoing large displacements according to the results obtained in the present work.
Future works will be focused on the development of shear-corrected UL finite element models to predict large twist rotations, where the model for the calculation of the shear correction factors incorporated into the variational formulation, accounts for the warping and high order shear strain energies.

Appendix
Explicit expressions for the constants for shear flow calculation in webs and flanges.

Figure 1
Figure 1 Coordinate systems of the Vlasov Warping Theory.

Figure 4
Figure 4 Comparison between analytical models and the actual numerical model for different lengths of beam.

Figure 5
Figure 5 Change of the error with the number of iterations for a determined load step.
the third iteration Latin American Journal of Solids and Structures 10(2013) 647 -673

Figure 7 Figure 8
Figure 7 Schematic setup configuration

Figure 10 Figure 9
Figure 10 Vertical displacements of the box beam using the not shear and shear corrected present model.

Figure 11 Figure 12 Figure 13
Figure 11 Variation of angle of twist in the longitudinal direction for several stacking sequences.

Table 1
Data of cantilever composite box beam

Table 2
Vertical displacement for the cantilever composite box beam problem being considered

Table 3
Specifications of the composites box beam.

Table 3
Comparison between numerical results an experimental data for the not shear model.

Table 4
Comparison between numerical results an experimental data for the shear corrected model.