Static analysis of composite beams on variable stiffness elastic foundations by the Homotopy Analysis Method

New analytical solutions for the static deflection of anisotropic composite beams resting on variable stiffness elastic foundations are obtained by the means of the Homotopy Analysis Method (HAM). The method provides a closed-form series solution for the problem described by a non-homogeneous system of coupled ordinary differential equations with constant coefficients and one variable coefficient reflecting variable stiffness elastic foundation. Analytical solutions are obtained based on two different algorithms, namely conventional HAM and iterative HAM (iHAM). To investigate the computational efficiency and convergence of HAM solutions, the preliminary studies are performed for a composite beam without elastic foundation under the action of transverse uniformly distributed loads considering three different types of stacking sequence which provide different levels and types of anisotropy. It is shown that applying the iterative approach results in better convergence of the solution compared with conventional HAM for the same level of accuracy. Then, analytical solutions are developed for composite beams on elastic foundations. New analytical results based on HAM are presented for the static deflection of composite beams resting on variable stiffness elastic foundations. Results are compared to those reported in the literature and those obtained by the Chebyshev Collocation Method in order to verify the validity and accuracy of the method. Numerical experiments reveal the accuracy and efficiency of the Homotopy Analysis Method in static beam problems.


Introduction
Laminated composite structures resting on elastic foundations are increasingly used in aerospace, marine, civil, biomedical, and other engineering applications due to their high strength-to-weight ratio, improved damage tolerance nature, and corrosion resistance. Composite beams are fundamental structural elements, and understanding of their static deflection behaviour is essential for engineering design and modelling purposes. Analytical solutions offer suitable tools to assist in this process and provide an insight into the governing physics of the problem allowing critical characteristics of composite structures on variable elastic foundations to be analysed and predicted. Different analytical approaches have been applied to analyse static and dynamic structural behaviour of beams resting on elastic foundations, for example, solutions such as that following Navier [1][2][3][4][5][6][7], Fourier series [8][9][10][11][12], Power series [13][14][15], Variational Iteration Method [16][17][18][19][20], and Homotopy Perturbation Method [19,[21][22][23]. However, these techniques suffer from significant drawbacks. For example, Navier's solution and Fourier series methods, in which all variables and their derivatives are expressed by trigonometric series expansions, are limited to specific simply supported boundary conditions. When applied to other boundary conditions, Fourier series tends to become slowly convergent, and its derivatives are not easily obtained and could even be discontinuous. These limitations can be overcome by modified Fourier series enriched with auxiliary polynomials noting that obtaining them poses additional difficulties. A power series solution involves discretisation of variables which usually requires large amounts of computer memory and processing time, leading to rounding errors, and thus results in loss of accuracy. The Variational Iteration Method requires finding the general Lagrange multiplier which is a challenging task in complex problems. Results provided by the Homotopy Perturbation Method highly depend on the good choice of a linear operator and initial guess. The Homotopy Analysis Method was proposed by Liao in [24] and further developed in [25] as an efficient analytical technique for solving nonlinear differential equations. The method employs the concept of the homotopy in topology and generates a solution in a form of convergent series. The main tenet of HAM is to replace a nonlinear equation by a system of linear ordinary differential equations that can be solved easily with the help of symbolic software packages such as MATLAB or Maple. In contrast to previously mentioned methodologies, HAM is an easy-to-use analytical tool which is general in terms of loading or boundary conditions. It does not involve discretisation or any physical parameters, small or large. Unlike other analytical approaches, HAM allows the convergence region and rate of convergence to be controlled by using an auxiliary parameter. Flexibility in the choice of this auxiliary parameter as well as linear operator and initial guess is another great advantage of the method. The Homotopy Analysis Method has been successfully applied to solve many problems in biology [26][27][28], chemistry [29][30][31], physics [32][33][34], optimal control theory [35][36][37], fluid mechanics [38][39][40][41], and solid mechanics [42][43][44]. More specifically, HAM successfully solved static and dynamic problems of isotropic beams. Wang et al. [45] used HAM to solve the bending problem of a cantilever beam under tip load resulting in large deformation. An accurate analytical expression for the rotation angle at the free end was obtained and then used to calculate the vertical and horizontal displacements of the beam. Kimiaeifar et al. [46] exploited HAM to solve analytically the large deflection and rotation of a nonlinear cantilever beam under the action of a co-planar follower static load. HAM results were compared with those obtained from a fourth-order Runge-Kutta method to prove the efficiency of the method. Later, Kimiaeifar et al. [47] used HAM to obtain an approximate solution for the governing equation of large deflection and rotation of a nonlinear Euler-Bernoulli beam with variable flexural rigidity loaded by circulatory forces. Maleki et al. [48] investigated the nonlinear large deformation of Euler-Bernoulli beams subject to arbitrary distributed loads by means of HAM. Deflection, slope, and bending moment were obtained for the case of constant and periodic distributed loads. Using HAM, Kimiaeifar et al. [49] derived an analytical solution for the large deflection analysis of an isotropic cantilever beam under the action of end point and uniformly distributed loads. Comparison of results with those obtained from the Finite Element Method (FEM) demonstrated the accuracy and computational efficiency of HAM. Liao [50] employed HAM to obtain convergent series solutions to eigenvalues and eigenfunctions for buckling of Euler-Bernoulli beams with uniform cross section under the action of axial load. Using the Galerkin method, Pirbodaghi et al. [51] transformed a nonlinear partial differential equation in space and time into an ordinary differential equation and then applied HAM to obtain analytical expressions for nonlinear vibration of Euler-Bernoulli beams subject to axial loads. By comparing HAM results against others in the literature, the accuracy of the method was demonstrated for simply supported and clamped beams. Hoseini et al. [52] employed HAM to obtain an accurate analytical solution for the fundamental nonlinear natural frequency and corresponding displacement of tapered beams presented by nonlinear differential equations of second order. Sedighi et al. [53] used HAM to derive analytical solutions for transversal vibration of hinged-hinged flexible beams subject to a constant load at their tips described by a fifth-order nonlinear differential equation. Excellent accuracy of HAM results was confirmed by comparing them against Runge-Kutta methods.
In the context of composite beams, Jafari-Talookolaei et al. [54] applied HAM to derive analytical expressions for the large amplitude free vibration of an unsymmetrically laminated composite beam on an elastic foundation under the action of axial load. An elastic foundation with quadratic and cubic nonlinearities was considered. Recently, Tang et al. [55] used HAM to find a closed-form solution for nonlinear free vibration of Euler-Bernoulli beams made of bi-directional functionally graded materials expressed by a third-order nonlinear partial differential equation which was transformed to ordinary differential equations (ODE) by means of the Galerkin method. The accuracy of the analytical results was confirmed by comparing them with those obtained by the Runge-Kutta method. Lin et al. [56] employed HAM to derive analytical expressions for the large deformation of axially functionally graded cantilever beams subject to tip loads. Results were compared with those from FEM to verify efficiency and accuracy.
Despite its advantages, HAM is a relatively new analytical approach in the context of composite structures, and the full potential of this method is not well known for solving various problems arising in this field. Mathematically, the coupled deflection of composite beams is a boundary value problem presented by a linear system of four coupled governing differential equations, reflecting four degrees of freedom, namely bending in two principal directions, axial elongation and twist [57][58][59][60]. In [61], this model was extended to composite beams resting on constant stiffness elastic foundations. A corresponding exact analytical solution was provided. In the current work, the model is further elaborated to composite beams resting on variable stiffness Winkler elastic foundations presented by a system of four coupled ordinary differential equations one of which has variable coefficients, noting that there is no exact solution for this newly developed problem. Thus, the main purpose of this work is to apply HAM, as a fresh attempt, to obtain a converged series solution for this problem and to investigate its accuracy, efficiency as well as factors affecting its convergence. To implement HAM, two different algorithms are used, one based on the traditional approach proposed in [24] and the other based on an iterative approach [62]. The accuracy and computational efficiency of both algorithms are evaluated by comparing their convergence and computational time required to obtain the results of the prescribed level of accuracy. In addition, to demonstrate the potential of HAM, results obtained from conventional and iterative approaches are compared with those obtained by the Chebyshev Collocation Method (CCM), which has proved to be very effective in solving beam problems [63][64][65][66].
The rest of the paper is organised as follows. In Sect. 2, a mathematical formulation of static deflection of composite beams resting on elastic foundations is given. In Sect. 3, a brief outline of the Homotopy Analysis Method is presented. In Sect. 4, the implementation of the method to the particular problem is explained. In Sect. 5, the validity and efficiency of the method are illustrated with some examples. A brief analysis of the obtained results is given. Finally, some discussions and conclusions are presented in Sect. 6. Figure 1 illustrates a slender composite beam of length measured along the x coordinate axis with the cross section in the y − z plane resting on a Winkler elastic foundation described by a set of mutually independent spring elements with modulus k w . Due to imperfections and material non-homogeneity of the foundation, the elastic coefficient of the Winkler foundation can vary along the beam span, i.e. k w = k w (x). The x-axis is assumed to be along the beam reference line.

Problem statement
Displacement of a generic point on the cross-section can be expressed as where U x , U y , and U z are the components of displacement vector; u, v, and w represent displacements of the beam reference line in x, y, and z directions; and ϕ, θ y , and θ z are the rotations of the beam cross-section about x, y, and z axes, respectively.
Using Eq. (1.1), the strain measures of the beam can be written as According to Euler-Bernoulli beam theory, the cross-section of the beam remains orthogonal to the reference axis after deformation; therefore, The internal work of the beam can be obtained by where σ xy , σ xz , and σ x x , are stresses in x, y, and z directions, respectively. Using definitions of the strain measures x x , γ xy , and γ xz given by Eqs. (2.1), Eq. (4) can be rewritten as (5) Internal forces and moments of a beam can be expressed as Considering Eqs. (3) and (6.1), Eq. (5) can be further modified as follows The principle of virtual work for a beam resting on an elastic foundation is written as where δW int , δW f , and δW ext are the variations in internal, elastic foundation, and external works, respectively. The variation in internal work of the beam can be obtained from Eq. (7) as follows: where vector of strains and curvatures ε, vector of internal forces and moments N, and stiffness matrix S can be defined as where x is a strain, κ x , κ y , and κ z are curvatures in x, y, and z directions, respectively, E A is the extensional stiffness, G J is the twist stiffness, E I y is the out-of-plane bending stiffness, E I z is the in-plane bending stiffness, S ET is the coupling between axial elongation and twist, S E F is the coupling between out-of-plane bending and axial elongation, S E L 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 F L is the coupling between out-of-plane and in-plane bending. It is worth mentioning that in Eq. (9) a linear relation between internal forces and moments N and strains and curvatures ε is assumed, i.e. N = Sε.
For more details of this assumption, see [60]. The variation of external work is presented by where the vector of displacements and rotations U and the vector of external forces and moments Q are expressed by where q x , q y , q z are functions of x representing distributed loads and q ϕ is the function of x representing distributed torque. Further, the variation of work due to the Winkler elastic foundation can be written as [3] 0 Substituting Eqs. (9), (11), and (13) into Eq. (8), applying the integration by parts and collecting the coefficients of δu, δϕ, δw, δv, δw , and δv , the set of governing equations can be obtained: At x = 0 and x = , the boundary conditions can be expressed as follows where () denotes the derivative with respect to x. The distribution of the variable stiffness Winkler foundation along the axial direction is assumed to be linear and varies according to the following rule: where k 0 is the stiffness of the foundation, 1 − ax is a function of the spatial coordinate along the beam length, and a is a variation parameter of the elastic foundation.

An overview of the Homotopy Analysis Method
Consider the differential equation where N is a nonlinear operator, x is an independent variable, and u(x) is an unknown function. According to Liao [25], the zeroth-order deformation equation can be written as where q ∈ [0, 1] is an embedding parameter, L is an auxiliary linear operator, φ(x; q) is the unknown function, u 0 (x) is the initial approximation of u(x), andh = 0 is an auxiliary convergence control parameter. Setting the embedding parameter q = 0 leads to and when q = 1, Eq. (19) is equal to the original governing equation meaning that as q increases from 0 to 1, the solution φ(x; q) varies from the initial guess u 0 (x) to the solution u(x). Expanding φ(x; q) in a Maclaurin series with respect to q and using Eq. (20) leads to where HAM provides great freedom in choosing the auxiliary parameterh. For an appropriately chosen value of h, the series presented by Eq. (22) converges when q = 1, and thus, where the unknown functions u n (x) can be obtained by using the nth-order deformation equation where and The nth-order deformation equation (25) is derived by differentiating the zeroth-order deformation equation (19) n times with respect to q, dividing by n! and then setting q = 0. Finally, applying the inverse operator L −1 , the expression for the unknown function can be obtained as follows: Substituting Eq. (29) into Eq. (24) leads to a series solution of Eq. (18).

Implementation of the Homotopy Analysis Method
In order to solve the static deflection of the composite beam described by the system of coupled Eq. (14) by means of HAM, consider the auxiliary linear operators of the form This selection is based on the method of highest order differential matching which implies the use of the highest order derivative appearing in the differential equation. In turn, inverse operators are integral operators given by where c i , i = 1, 2, . . . , 12, are arbitrary constants of integration to be determined from boundary conditions. Defining nonlinear operators as

4)
nth-order deformation equations for the problem can be written as follows: The next important step in employing the Homotopy Analysis Method is to choose the initial solution. Appropriate selection of the initial approximation (i.e. as close to the exact solution as possible) results in better accuracy and convergence of the approximate solution. A suitable choice for the static deflection problem of composite beams is to use the solution of uncoupled equations of isotropic beams as for different types of boundary equations these solutions are expressed by polynomial functions (as discussed in Sect. 5). The polynomial functions are easy to integrate, which in turn simplifies the process of deriving the nthorder deformation equations and reduces the computational time required to obtain the approximate solutions. However, as the order of deformation equation increases, the computational complexity grows rapidly. To overcome this issue and to increase the convergence and accuracy of HAM, Liao [62] suggested to update the initial guess at each iteration by using the solution of the previous step as the initial guess of the current step. This approach is known as iterative HAM. In Sect. 5, both traditional HAM and iterative HAM are used to obtain the analytical solutions of static deflection of composite beams with and without elastic foundations. The efficacy of iHAM over HAM in solving beam static deflection problems is demonstrated, and then, the convergence and accuracy of the approach are investigated by comparing convergence rates and CPU time required to obtain results of specific levels of precision.

Preliminary studies
Preliminary studies of the convergence and computational efficiency of both approaches of HAM, conventional and iterative, are performed in this Subsection. For that purpose, a cantilever laminated fibre-reinforced slender beam without elastic foundation (i.e. k w (x) = 0) with boundary conditions expressed by Eq. (15.1) at x = 0 and Eq. (16.1) at x = is considered. To investigate the effect of material properties on the convergence of the solution, three types of the layup introducing different combinations of coupling terms are assumed, namely symmetric [45 3 ] s with bend-twist coupling, cross-ply [0 3 /90 3 ] with axial-bend coupling, and nonsymmetric [60 3 /30 3 ] with bend-twist, axial-twist, and axial-bend coupling. The corresponding stiffness matrices, obtained by using closed-form expressions developed by Yu and Hodges [67], are presented in Table 1.
The following material properties are assumed: E 11 = 135.64 GPa, E 22 = 10.14 GPa, G 12 = 5.86 GPa, and ν 12 = 0.29. While the proposed mathematical model presented by Eq. (14) allows three types of distributed loads along the beam axis and torque, in the numerical examples only transverse load q z acting in the vertical direction is considered and used to normalise the deflections as follows: The following solutions of the uncoupled equations of isotropic beam are used as initial guesses: Before comparing the two HAM approaches, an appropriate order of deformation equation for iHAM should be established. For that purpose, iHAM results are obtained for a symmetric cantilever beam under the action of uniformly distributed load for different orders of deformation equation, namely N = 2, N = 3, and N = 4, using Eq. (35.1) as initial guesses andh = −1. Using these results and those obtained from the exact solution [60], the relative error was calculated as Figure 2 demonstrates that for both degrees of freedom available for symmetric layup, the log of error is decreasing rapidly when N = 2, while increasing the order of deformation equations up to N = 3 or N = 4 leads to divergent results. Thus, using the second-order deformation equation is a natural choice, and iHAM results throughout the paper are obtained based on this assumption.
With this in mind, the convergence of both algorithms of HAM can be investigated. As shown in Fig. 3, both traditional and iterative approaches provide converged results for all three types of stacking sequence. Generally, it is observed that better accuracy is obtained for symmetric and cross-ply layups compared to the non-symmetric layup which can be attributed to the occurrence of more coupling terms resulting in higher complexity. For example, the 10th-order deformation equation is sufficient to match exact results up to three decimal places in the case of the symmetric layup. In contrast, 25th-order and 35th-order deformation equations are required to attain the same level of accuracy in the cross-ply and non-symmetric layups, respectively. It is also worth noting that the convergence of iHAM results exhibits more stable behaviour compared to HAM results.
In order to demonstrate the computational efficiency of conventional HAM and iHAM, CPU time versus log of relative error is depicted in Fig. 4. It is worth noting that the calculations are performed using a 1.8 GHz Intel Core i5 system. Both conventional and iterative approaches allow results up to a specific degree of precision to be achieved. However, it is observed that for the same level of accuracy iHAM needs less CPU time than the traditional HAM. It is noted that iHAM converges approximately three times as fast as HAM.
The nondimensionalised transverse bending deflections for symmetric, cross-ply, and non-symmetric layups are shown in Fig. 5. Table 2 provides numerical results for the tip deflections of all three types of layups. It is shown that iHAM results are in good agreement with those obtained from the exact solution.

Verification studies
In this Subsection, the validity and potential of the Homotopy Analysis Method in solving the static deflection problem of beams resting on elastic foundations are verified. Since numerical results for bending of fully coupled composite beams on variable elastic foundations are not available in the literature, static deflection analysis of homogeneous isotropic and composite beams on constant stiffness elastic foundations considered in several well-known studies is performed, and the results are compared with analytical and numerical results existing in the literature.
First, values of mid-span deflection of a uniformly loaded isotropic beam with clamped-clamped and simply supported-simply supported boundary conditions obtained from the present HAM and iHAM solutions, exact analytical solution [61], analytical solution based on Green's functions [68], and the DQM solution [69] are presented in Table 3. For convenience, the mid-span transverse deflection and Winkler elastic foundation parameter are normalised as follows: where E I y is the flexural rigidity and the length-to-thickness ratio / h = 120. The analysis has been performed for three different values of k w resulting in an elastic foundation of different stiffnesses, namely k w = 0 meaning there is no elastic foundation; k w = 10 and k w = 100 describing stiffer foundations. The results obtained from different theories for different values of the elastic foundation parameter k w are in close agreement. It is worth mentioning that HAM results were obtained using N = 20, while iHAM results were obtained using second-order deformation equations with 8 iterations. In both approaches,h = −1 was assumed. Next, nondimensional mid-span deflection of a uniformly loaded clamped laminated composite beam resting on a constant stiffness elastic foundation is obtained considering symmetric [45 3 ] s , cross-ply [0 3 /90 3 ], and unsymmetric [60 3 /30 3 ] stacking sequences. Results obtained based on the conventional and iterative algorithms of HAM, CCM, and from the exact solution [61] are given in Table 4 for normalised Winkler elastic foundation parameter k w taking values 0, 10, and 100. It is noted thath i = −1, i = u, ϕ, w, v, was used for both HAM approaches. The order of deformation equation required to obtain conventional HAM results of reasonable accuracy for a beam with symmetric stacking sequence was N = 12. With the increasing complexity of the combination of coupling terms in cross-ply and unsymmetric layups, the order of deformation equation was increased up to N = 20 and N = 26, respectively. Similarly, the number of iterations in iHAM was growing from 6 for the case of symmetric layup up to 10 and 12 for cross-ply and unsymmetric cases, respectively. It is worth mentioning that for iHAM the length of the homotopy approximations increases at each iteration exponentially which may affect the computational efficiency of the method. To address this issue, Liao [62] proposed that an approximate solution could be truncated. However, this strategy was not applied in the current work since CPU times required were acceptable. The results obtained from the exact, analytical, and numerical solutions are in good agreement.

Composite beam on variable stiffness elastic foundation
In the previous Subsections, the efficiency of HAM in solving static problems of beams resting on elastic foundations was demonstrated. In addition, the superiority of an iterative approach over a conventional one was shown. Based on these results, in this Subsection iHAM analysis is extended to a composite beam resting on variable stiffness Winkler elastic foundation. The stiffness of this foundation varies linearly according to the rule described by Eq. (17) and is normalised as  equations of isotropic beams are obtained and used as initial guesses, To study the effect of Winkler elastic coefficient k 0 on the static behaviour of a composite beam, the deflections of a beam with various stacking sequences were analysed assuming the linear variation parameter a = 0 while k 0 was allowed to vary from 0 to 100, noting k 0 = 0 describes the beam without elastic foundation, whereas k 0 = 0 corresponds to a Winkler elastic foundation of constant stiffness. Table 5 gives a comparison between the analytical iHAM results and numerical approximations obtained by CCM. The nondimensionalised values of deflections obtained from both methods are in good agreement. It is worth noting that to obtain iHAM results for all values of k 0 , auxiliary parametersh u ,h ϕ , andh v equal to −1 were used, whileh w was assumed to be equal to −1 for k 0 = 0 and k 0 = 10, and to −1/2 for k 0 = 50 and k 0 = 100, respectively. It is worth mentioning that there are different approaches for finding appropriate values ofh, for example, plotting the so-calledh-curve which allows the valid region of the convergence control parameter to be identified, or making use of the traditional square residual error employing integration or corresponding discrete form (for more details see [62]). However, for complicated coupled problems plottingh-curves may be computationally inefficient, while the exact integration approach involved in the square residual error may not be possible or computational cost of its discrete form may not be acceptable. Thus, in the current paper the trial-and-error approach is used to find the appropriate convergence control parameterh. Also it should be noted that more iterations of iHAM are required as the value of k 0 increases. The deformed configurations of the beam for all three types of the stacking sequence for different values of k 0 are depicted in Fig. 6. It is observed that the level of deformation depends on the type of the layup, consistently increasing from symmetric to non-symmetric and then to cross-ply stacking sequences. Additionally, the inverse relationship between the value of k 0 and the amplitude of the beam deflection can be observed.
To investigate the effect of the linear variation parameter a on the deflection of the composite beam resting on a variable stiffness elastic foundation, Winkler elastic coefficient k 0 was assumed to be 10 while a was allowed to vary from 200 to 600. Using iHAM and CCM, the nondimensional deflections were calculated for each  Table 7. The results obtained from both methods are generally in good agreement. It is worth noting that the level of accuracy of iHAM results presented in Table 7 is significantly influenced by the choice of auxiliary parameters, and more iterations are required to obtain the accurate iHAM results as the value of the linear variation parameter a increases. For the current example values ofh i , i = u, ϕ, w, v, corresponding to different degrees of freedom, are presented in Table 6. The deformed configurations of simply supported composite beams resting on linearly variable elastic foundations for different types of the layup are shown in Fig. 7. It is shown that the amplitude of the deflection depends inversely on the value of the coefficient of the variable term, namely linear variation parameter a. In addition, the value of this parameter affects the symmetry of the beam deformation, meaning the higher the value of a the more asymmetric the shape.

Conclusions
An explicit analysis of the static deflection of composite beams resting on variable stiffness elastic foundations subject to uniformly distributed loads described by a non-homogeneous system of coupled differential equations with the combination of constant and variable coefficients has been performed by means of the Homotopy Analysis Method. Two HAM algorithms, i.e. conventional and iterative, have been implemented, and the applicability and accuracy of the method were investigated. The main difference between the two algorithms is in the method of achieving more accurate results: in conventional HAM the order of deformation equation increases for this purpose, while in iHAM the initial guess is updated at each iteration. Numerical simulations of the static deflection of composite beams with various coupling terms without elastic foundations were performed to demonstrate the superiority of the iterative approach over the conventional type by the comparison of their convergence rates. It was shown that iHAM converges more rapidly (up to three times faster than HAM) and exhibits more stable behaviour. Comparison of the results obtained by HAM and iHAM with those available in the literature shows the viability of using both approaches in the static analysis of beams resting on elastic foundations. Finally, iHAM was applied to the static analysis of composite beams resting on variable stiffness Winkler elastic foundations with different values of elastic foundation coefficient. The ensuing analytical results were compared against those obtained from the Chebyshev Collocation Method, and excellent agreement between them was observed. It was also noted that the convergence of HAM results correlates with the complexity of combination of coupling terms provided by different stacking sequences Values of the elastic foundation coefficient are another factor significantly influencing the computational performance of the method. It was observed that higher values of these coefficients require a larger number of iterations in obtaining iHAM results of reasonable accuracy. Thus, it can be concluded that the iterative Homotopy Analysis Method is a convenient and efficient method for the static analysis of composite beams resting on variable stiffness elastic foundations.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.