A semi-analytical approach based on the variational iteration method for static analysis of composite beams

Using the Variational Iteration Method (VIM) the 3D static deflection problem of composite beams subject to concentrated tip and uniformly distributed loads is analysed, resulting in a system of coupled non‐ homogeneous ordinary differential equations. Using a general Lagrange multiplier, identified by variational theory, a special type of functional is constructed. By making an initial approximation in the form of a Maclaurin series and by using successive iterations, the solution in the form of convergent series is obtained. The results based on VIM are compared against those of the exact solution and Chebyshev Collocation Method (CCM) for different layups and boundary conditions and good agreement is observed between them. These results show the applicability and effectiveness of VIM for the static analysis of composite beams.


Introduction
Due to their nature, analytical solutions which are expressed symbolically and in closed-form, are indisputably important for understanding physical effects as well as for design purposes. Analytical solutions quantify the importance of model parameters directly, compared to numerical methods, while also being free from numerical instabilities.
The main advantage of VIM is the solution of a problem in closed form obtained without discretising the variables or using small parameters in the differential equation. As a result, the solution is not affected by computational rounding errors. In addition, VIM does not require restrictive assumptions for the nonlinear terms such as being differentiable with respect to the dependent variable which allows complicated analytic calculations to be avoided. For most linear cases the first iteration of the bespoke functional leads to the exact solution using the exact Lagrangian multiplier. One of its main merits is the ability to provide a solution in symbolic form, but this feature can also become a disadvantage if excessive symbolic calculations are involved, making calculations difficult and time consuming. VIM can also have difficulties with evaluation of its Lagrange multiplier and the uncontrollability of non-zero boundary conditions. VIM has been used for the linear and nonlinear static analysis of beams. Ozer [20] used VIM to solve the deflection problem of simply supported Euler-Bernoulli beam with an internal hinge on a uniform elastic foundation under a uniform distributed load and clamped-clamped column with two internal hinges, both subject to jump discontinuities. Results were compared with analytical solutions available in the literature and shown to be in a good agreement. Martin [21] applied VIM to obtain the quasi-static exact solution for simply supported isotropic linear viscoelastic beams under a uniform distributed load. By applying a modified variational iteration method which involves Laplace transforms, the dynamic case of the problem was also analysed. Salehi et al. [22] used VIM to study large deflections of a cantilever beam subject to a tip loads. Results were compared with those obtained from the Differential Transformation Method (DTM) and the Adomian Decomposition Method (ADM). Baghani [23] applied the modified VIM using Laplace transforms to obtain an analytical solution for the size-dependent static nonlinear deflection of cantilever Euler-Bernoulli microbeams under transverse distributed electrical forces. The results were compared with experimental data and with results obtained from the Finite Difference Method and found to be in good agreement. Ghaffarzadeh and Nikkar [24] employed a modified version of VIM which involves the use of Laplace transforms to obtain an explicit analytical expression for the large deformation of a cantilever beam under a tip load. The accuracy and convergence of the method were investigated and compared with those of the fourth-to fifthorder Runge-Kutta-Fehlberg method as well as VIM and Homotopy Analysis Method. Eroglu [25] used VIM to study large in-plane deflections of planar curved beams made of functionally graded materials with material properties considered as an arbitrary function of the position on the cross-section of beam. Results obtained from VIM were compared with those existing in the literature.
In the context of composite beams, VIM has been applied for solving eigenvalue problems involving buckling and free vibration analyses [26][27][28][29][30][31][32][33][34][35]. However, the applicability of this method for the deflection of fully coupled composite beams has not yet been examined. Such types of deflection are described by boundary value problems, governed by a system of coupled fourth-order ordinary differential equations [36,37]. Depending on the type of applied load, this system of equations can be homogeneous or non-homogeneous. Usually obtaining the exact solution of such coupled systems of equations is mathematically demanding and a time consuming process. Besides, exact results can be difficult to interpret. In this situation obtaining fast semi-analytical solutions of reliable accuracy in a convenient form is desired. Consequently, the current main goals are to apply VIM as a fresh attempt and to investigate its accuracy and benefits for the static analysis of fully coupled Euler-Bernoulli composite beams. The accuracy of VIM is investigated by comparing the obtained results with an exact solution as well as the numerical solutions based on the Chebyshev Collocation Method which has been shown to be accurate and efficient in beam problems [38][39][40][41]. In the current paper, the accuracy of VIM is shown to depend on the combination of coupling terms in the stiffness matrix and also on the type of loading conditions. This article is organised as follows: Section 2 presents the governing equations and corresponding boundary conditions of composite beam; Section 3 provides a clear overview of the Variational Iteration Method; Section 4 describes how VIM is applied to the problem of static deflection of composite beam; Sections 5 and 6 demonstrate the effectiveness of the method with various examples before a brief conclusion is given in Section 7. To keep the formulation as general as possible, the stiffness matrix of this beam is expressed in engineering constants as follows ; ð1Þ

Problem statement
where EA is the extensional stiffness, GJ is the twist stiffness, EI y is the out-of-plane bending stiffness, EI z is the in-plane bending stiffness, S ET is the coupling between axial elongation and twist, S EF is the coupling between out-of-plane bending and axial elongation, S EL is the coupling between in-plane bending and axial elongation, S FT is the coupling between out-of-plane bending and twist, S LT is the coupling between in-plane bending and twist, and S FL is the coupling between out-ofplane and in-plane bending. The set of governing equations for the beam can be written as [36] À EAu 00 À S ET φ 00 þ S EF w 000 À S EL v 000 ¼ q x À S ET u 00 À GJφ 00 þ S FT w 000 À S LT v 000 ¼ q φ At x ¼ 0 and x ¼ ', boundary conditions can be expressed as follows where ðÞ0 denotes the derivative with respect to x; the displacements of the beam reference line in x; y and z directions are denoted as u; v and w, and φ is the twist of beam cross-section about x; q x ; q y ; q z are the distributed loads and q φ is the distributed torque; b f x ; b f y ; b f z are the tip loads in the x; y and z direction respectively, and b m x is the tip torque, and b m y ; b m z are the tip moments about the y and z axes respectively. The exact solution of Eqs.
(2) is provided in Appendix A.

The concept of Variational Iteration Method
To understand the concept of Variational Iteration Method, consider the general nonlinear equation where L is a linear operator, N is a nonlinear operator and g x ð Þ is an inhomogeneous source term.
A functional for Eq. (4) is written as [15] where λ is a general Lagrange multiplier, the subscript n denotes the nth approximation and u ∼ n is restricted variation meaning δu ∼ n ¼ 0. He [1] originally called this a "correction functional" but we prefer to call it the VIM functional to recognise its importance as a special kind of augmented functional.
The identification of the Lagrange multiplier λ is a key step in VIM. Variational theory is usually used for this purpose. To start the iteration process, an initial solution for u 0 x ð Þ is proposed. Usually, it is presented by any arbitrary continuous function with some unknown parameters which satisfy initial and/or boundary conditions. Usually, the type of polynomial selected for this purpose affects the convergence process but not the final result and Maclaurin series is used to expand the initial solutions [13]  where M is the degree of Eq. (4). Once the Lagrange multiplier and initial solution have been identified, iteration formula Eq. (5) can be applied and the solution for Eq. (4) can be obtained as In practice, n is a finite integer chosen depending on the required accuracy of result.

Implementation of Variational Iteration Method
The VIM solution procedure starts with constructing the VIM functionals for the system of governing equations, described by Eqs.
where λ u ; λ φ ; λ w and λ v are general Lagrange multipliers, and u Taking variations of Eqs. (8) with respect to u n ; φ n ; w n and v n respectively leads to and further to Integrating Eq. (10a) by parts, then Thus, collecting like-terms for δu n , the stationary conditions for this case can be written as Integrating Eq. (12a) twice, the following expression for λ u ξ ð Þ can be obtained: where constants of integration c 1 and c 2 can be found by imposing conditions given by Eqs. (12b) and (12c) which results in Evidently, λ φ ¼ λ u . In a similar manner, integrating Eq. (10c) by parts and keeping in mind that δw n 0 ð Þ ¼ 0 results in Collecting like-terms for δw n , corresponding equations for stationary condition can be obtained δw n x ð Þ : Integrating Eq. (16a) results in where c i ; i ¼ 1 . . . 4, are constants of integration which can be obtained by applying conditions given by Eqs. (16b)-(16e). Thus, the Lagrange multiplier can be identified in the following form Substituting the identified Lagrange multipliers, presented by Eqs. (14) and (18), into Eqs. (8) results in the following iteration formulae To apply these iteration formulae, initial solutions are required. Using Eq. (6), u 0 and φ 0 can be obtained by using Maclaurin series with the first two terms while w 0 and v 0 can be obtained by using Maclaurin series with the first four terms being where a i ; b i ; i ¼ 1; . . . ; 4; c j and d j ; j ¼ 3; 4, are unknown constants to be determined according to the applied boundary conditions. Once iteration formulae and initial solutions are established, VIM can be applied to analyse the static deflection of fully coupled Euler-Bernoulli composite beams.
In the following sections results for composite beams under the action of concentrated tip loads and uniformly distributed loads subject to various boundary conditions are presented.

Verification studies
To demonstrate the ability of VIM to predict the static response of composite beams accurately, the comparison of numerical results with some theoretical and experimental developments available in the literature are performed.
First, to verify the numerical results obtained from VIM, a number of test cases of static deflection of laminated composite beams subject  to uniformly distributed load q z are presented. For convenience, the deflection is normalised as follows where b is the width and h is the thickness of rectangular beam crosssection. Deflection is obtained for a slender composite beam ('=h ¼ 50) with the symmetric 0 =90 =0 ½ and anti-symmetric 0 =90 ½ types of layup for various boundary conditions: hinged-hinged (H-H), clamped-hinged (C-H), clamped-clamped (C-C), clamped-free (C-F). VIM results are compared against the Chebyshev Collocation Method (for more details see [36,37]), and some other results available in the literature. The following dimensionless orthotropic material properties are considered: Tables 1,2 demonstrate good agreement between VIM results and the results from the exact solution of Khdeir and Reddy [42] based on Euler-Bernoulli theory, trigonometric series analytical solution of Nguyen et al. [43] and Finite Element Method (FEM) results of Murthy et al. [44] for symmetric and anti-symmetric beams respectively. Now consider a thin-walled composite box-beam with the so-called Circumferentially Uniform Stiffness (CUS) exhibiting axial elongationtwist coupling, subject to clamped-free boundary conditions and a unit torsional tip load (1 in.-lb). Each wall of the box-beam consists of six layers of 15 graphite-epoxy unidirectional prepregs (AS4/3501-6). The stiffness properties of this box-beam are calculated using FEM cross-sectional tools and reported in Table 9 of [45]. A static test was conducted on this beam and repeated twice as described in [46]. The experimental data on twist is compared with the analytical results obtained using the exact solution and one iteration of VIM. As shown in Fig. 2, analytical estimations of twist are within the scatter of the experimental data which demonstrates the efficacy of the proposed approach. Most of the studies of static deflection of beams available in the literature examine isotropic and symmetrically laminated composite beams. However, due to the presence of axial-bending coupling terms in unsymmetrically laminated beams their analysis is significantly more complicated. To demonstrate this fact, symmetric 45 3 ½ s , crossply 0 3 =90 3 ½ and unsymmetric 60 3 =30 3 ½ stacking sequences are considered in this paper. Corresponding stiffness matrices are presented in Table 3. Stiffness terms used in Eq. (2) are obtained using the closed-form expressions based on classical lamination theory presented by Yu and Hodges [47].
The analysis is performed for normalised deflections. In the case when distributed load q z is applied to a composite beam, deflections are defined as In the case of concentrated end load b f z , normalised deflections are Table 4 and Fig. 3 demonstrate how VIM efficiently produces the values for the tip deflection of cantilever composite beams with symmetric, cross-ply and unsymmetric stacking sequences under the action of concentrated load b f z ¼ 0:01 N applied at the free end. For the case of symmetric and cross-ply layups just one iteration of VIM is required to obtain results which are in an excellent agreement with those obtained from the exact solution and CCM with 6 shape functions. Due to the presence of extra coupling terms, five iterations are needed for the case of unsymmetric stacking sequence to obtain reasonably accurate results. In order to verify the efficiency of the variational iteration method in comparison with exact solution and the Chebyshev Collocation Method, relative errors for CCM and VIM are reported as The rapid convergence of VIM solution for all types of layup is depicted in Figs. 4-6. Table 5 and Fig. 7 show the maximum deflection which occurs at the free end of cantilever composite beams with symmetric, crossply and unsymmetric stacking sequences under the action of distributed load q z ¼ 0:01 N=m. In these cases, loads are included in the governing equations and so iteration formulae are more complicated than in the previous example. As a result, 10 iterations should be conducted with VIM in the case of symmetric layup, 20 iterations in the case of cross-ply layup and 35 iterations in the case of unsymmetric stacking sequence to obtain results which are in a good agreement with results from exact and CCM solutions. The convergence of VIM solution for symmetric and unsymmetric cases is graphically presented in Figs. 8-10 respectively. Now consider composite beams with the same geometrical and material properties simply supported at both ends subject to the uniformly distributed load q z ¼ 0:1 N=m. Non-zero deformations of this beam are depicted in Fig. 11. In this case maximum deflection for w occurs at the mid-span of beam, while maximum values for other degrees of freedom can be observed at x≈ 1 4 ' and x≈ 3 4 '. The difference between the 10-iteration VIM (for symmetric layup), 16-iteration VIM (for cross-ply layup) and 35-iteration VIM (for unsymmetric layup), 6shape function of CCM and exact solution are presented in Table 6 and shown to be in a good agreement.

Conclusions
The Variational Iteration Method has been applied to analyse the static deflection of composite Euler-Bernoulli beams described by a system of fully coupled non-homogeneous ordinary differential equations. To perform analysis, a special functional is constructed using general Lagrange multipliers. To identify these multipliers, variational theory is applied. Once the initial approximation with unknown constants is selected by using Maclaurin series, several successive approximations of the VIM functional provide an analytical solution with reasonable accuracy. Analyses have been performed for fully coupled composite beams subject to clampedfree and simply supportedsimply supported boundary conditions under the action of concentrated end load and transverse uniform distributed load for different types of stacking sequence. Results obtained are compared to the exact solution and those provided by the Chebyshev Collocation Method. The convergence of VIM solutions was investigated and it is shown that the accuracy of VIM depends on the loading conditions and the type of coupling terms. For example, just one iteration with VIM is needed to obtain the exact results in the case of symmetric composite beam subject to tip load, while due to the presence of extra coupling terms five iterations are required to obtain the precise results for the case of the beam with unsymmetric layup. The fact that distributed loads are included in the governing equations increases the complexity of iteration formula. As a result, VIM solutions as a sequence of iterations requires a large number of iterations to give reasonable accuracy (10, 20 and 35 iterations for the symmetric, cross-ply and unsymmetric layups respectively). This could lead to a high degree of complexity so the resulting integrations in the relations may be difficult to perform. It is worth noting, that the convergence of method is independent of the choice of boundary conditions. VIM can be considered as a method which provides reasonable analytical approximations to the solution of coupled ordinary differential equations.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.