The aim of this paper is to present analytical and exact expressions for the frequency and buckling of large amplitude vibration of the symmetrical laminated composite beam (LCB) with simple and clamped end conditions. The equations of motion are derived by using Hamilton's principle. The influences of axial force, Poisson effect, shear deformation, and rotary inertia are taken into account in the formulation. First, the geometric nonlinearity based on the von Karman's assumptions is incorporated in the formulation while retaining the linear behavior for the material. Then, the displacement fields used for the analysis are coupled using the equilibrium equations of the composite beam. Substituting this coupled displacement fields in the potential and kinetic energies and using harmonic balance method, we obtain the ordinary differential equation in time domain. Finally, applying first order of homotopy analysis method (HAM), we get the closed form solutions for the natural frequency and deflection of the LCB. A detailed numerical study is carried out to highlight the influences of amplitude of vibration, shear deformation and rotary inertia, slenderness ratios, and layup in the case of laminates on the natural frequency and buckling load.

1. Introduction

The uses of composite materials in structural components are increasing due to their attractive properties such as high strength/stiffness, light weight, and facility to materialize fiber orientation, material and stacking sequence. Beam-like structures have numerous applications in the industries such as aerospace and robotics. The increased use of the LCB requires a better understanding of large amplitude vibration aspects of them.

Nonlinear dynamic analysis of the LCB on the basis of classical lamination theory (CLT) has received a good amount of attention in the literature [1ā€“6] while relatively little investigations [7, 8] have been performed for dynamic investigation of such beams with low slenderness ratio (length of the beam to the radius of gyration of the cross-section i.e. = šæ/š‘Ÿ) which have to be taken into account the shear deformation and rotary inertia. It should be noted that for the later, the equation of motion is complex, and even obtaining an approximate solution is much more difficult.

In our previous work in [1], we presented analytical expressions for large amplitude-free vibration and postbuckling analysis of unsymmetrically LCB on elastic foundation by using variational iteration method. We used the CLT to obtain differential equations of motion. A differential quadrature approach for nonlinear free vibration analysis of symmetric angle-ply laminated thin beams on nonlinear elastic foundation with elastically restrained against rotation edges has been studied by Malekzadeh and Vosoughi [2]. Patel et al. [3] have extended a new three-nodded shear flexible beam element based on consistency approach to analyze the nonlinear flexural vibration and postbuckling behavior of simply supported laminated orthotropic composite beams on a two-parameter elastic foundation. The axial and transverse deformations of a geometrically nonlinear composite beam were considered to obtain a closed form solution for the postbuckling behavior of composite beams and to exactly solve for the free vibrations that take place around a buckled position by Emam and Nayfeh [4]. Kapania and Raciti [5] presented a two-nodded Timoshenko beam element with 10 degrees of freedom per node to examine the large amplitude free vibrations of LCB. The nonlinear free flexural vibrations of isotropic/laminated orthotropic straight/curved beams have been studied using a cubic B-spline shear flexible curved element [6]. The nonlinear governing equations are solved by employing Newmarkā€™s numerical integration scheme coupled with modified Newton-Raphson iteration technique. Kiran et al. [7] had investigated the large amplitude free vibrations of generally layered laminated composite beams by developing a finite element model. A one-dimensional finite element formulation based on a higher-order displacement model has been developed by Obst and Kapania [8]. The model accounts for geometric nonlinearities, a parabolic shear strain distribution through the thickness, and satisfies the shear stress free boundary conditions at the upper and lower surfaces of the beam. The model has been applied to the nonlinear static and transient analysis, free vibration analysis, and impact analysis of laminated beams.

To the best author's knowledge, there is no analytical expressions in the literature on the large amplitude free vibrations analysis of the symmetrically LCB on the basis of first shear deformation theory (FSDT). Thus, the contribution of this paper is to expand our previous work in [1] by studying the more complicated and important problem to obtain the closed-form expressions for natural frequency and deflection.

2. Equations of Motion and Boundary Condition

The considered LCB with simple end conditions is shown in Figure 1 with the coordinates Ģ‚š‘„ along the axis of the beam and Ģ‚š‘§ along the thickness direction. īĢ‚Ģ‚š‘”šæ,š‘ƒ,š‘¤(Ģ‚š‘„,š‘”), are length, compressive axial load, deflection of the beam, and time, respectively.

2.1. Coupling Equations

Based on FSDT, the LCB is characterized by its bending rigidity EI, torsional rigidity GJ, bending-torsion coupling rigidity C, and shear rigidity kAG [9]. The derivation of these rigidities is briefly summarized in Appendix A [10].

The kinetic and potential energies of the LCB and the work done by the static lateral load per unit length of the beam š¹z can be written as1š‘‡=2ī€œšæ0ī‚ƒš¼1īš‘¤,Ģ‚š‘”2+š¼2šœƒ,Ģ‚š‘”2+š¼š‘šœ“,Ģ‚š‘”2ī‚„1š‘‘Ģ‚š‘„,š‘ˆ=2ī€œšæ0ī‚øEIšœƒ,Ģ‚š‘„2ī‚µīš‘¤āˆ’š‘ƒ,Ģ‚š‘„2+š¼š‘š¼1šœ“,Ģ‚š‘„2ī‚¶+2š¶šœ“,šœƒĢ‚š‘„,ī€·īš‘¤Ģ‚š‘„+kAG,ī€øĢ‚š‘„āˆ’šœƒ2+GJšœ“,Ģ‚š‘„2ī‚¹ī€œš‘‘Ģ‚š‘„,š‘Š=šæ0š¹š‘§īš‘¤š‘‘Ģ‚š‘„,(1) where comma denotes derivative with respect to š‘„; šœƒ,šœ“, are the angle of rotation and torsional rotation, respectively; š¼1,š¼2 and š¼p are the first and second moment of inertia about š‘¦-axis and the polar mass moment of inertia per unit length about the š‘„-axis, respectively.

For the subsequent results to be general, we use the following nondimensional variables:š‘„=Ģ‚š‘„šæīš‘¤,š‘¤=š‘ŸĢ‚š‘”,š‘”=ī”š‘š·11/ī€·š¼1šæ4ī€øī‚™,š‘Ÿ=š¼š“(2) which š‘Ÿ,š‘,š¼, and š“ are the radius of gyration of the cross-section, beam width, second moment of area, and the area of the cross-section, respectively, and š·11 is defined in Appendix A.

Thus, the kinetic and potential energy and the work done can be rewritten as follows:š¼š‘‡=1šæ52š‘š·11ī€œ10ī€ŗš¼1š‘Ÿ2š‘¤2,š‘”+š¼2šœƒ2,š‘”+š¼š‘šœ“2,š‘”ī€»1š‘‘š‘„,š‘ˆ=2ī€œ10ī‚øEIšæšœƒ2,š‘„āˆ’š‘ƒšæī‚µš‘Ÿ2š‘¤2,š‘„+š¼š‘š¼1šœ“2,š‘„ī‚¶+2š¶šæšœ“,š‘„šœƒ,š‘„+kAGšæī€·š‘Ÿš‘¤,š‘„ī€øāˆ’šæšœƒ2+GJšæšœ“2,š‘„ī‚„ī€œš‘‘š‘„,š‘Š=10š¹š‘§š‘¤š‘‘š‘„.(3) Applying the principle of minimization of total potential energy (š›æ(š‘ˆāˆ’š‘Š)), we obtain the following equilibrium equations as follows:EIšæ2šœƒī…žī…žī‚€š‘Ÿ+kAGšæš‘¤ī…žī‚+Cāˆ’šœƒšæ2šœ“ī…žī…ž=0,(4)ī‚µš‘ŸkAGšæ2āˆ’1šæšœƒī…žī‚¶āˆ’š‘ƒš‘Ÿšæ2š‘¤ī…žī…ž+š¹š‘§=0,(5)GJšæ2šœ“ī…žī…ž+Cšæ2šœƒī…žī…žāˆ’š‘ƒš¼š‘ƒš¼1šæ2šœ“ī…žī…ž=0(6) in which prime denotes derivative with respect to nondimensional coordinate š‘„. These equations are coupled equations and can be solved to get the coupled displacement fields [11].

2.2. Coupled Displacement Fields

Equations (4) and (6) are used to couple the rotation šœƒ, torsion šœ“, and deflection of the beam š‘¤ [11]. These equations can be rewritten as šœ“ī…žī…ž=āˆ’š“1šœƒī…žī…ž,šœƒī…žī…ž+š“2š“šœƒ=2šœ‚š‘¤ī…ž,(7) whereš“1=š¶ī€·GJāˆ’š‘ƒš¼š‘ƒ/š¼1ī€ø,š“2=kAGšæ2Cš“1āˆ’EI,(8) and šœ‚ is the slenderness ratio, that is, šœ‚=šæ/š‘Ÿ.

2.2.1. Simple End Conditions

The boundary conditions for the simple end conditions beam are as follows:(š‘¤=0,š‘€=0,š‘‡=0)atš‘„=0,1āŸ¹ī€·š‘¤=0,šœƒī…ž=0,šœ“ī…žī€ø=0atš‘„=0,1.(9) The transverse displacement š‘¤ is assumed asš‘¤(š‘„,š‘”)=š‘¢(š‘”)sin(šœ‹š‘„),(10) where š‘¢(t) is the amplitude parameter and only is a function of š‘”. Substituting the function š‘¤ in (7), we obtain the coupled rotation šœƒ and torsion šœ“ asšœƒ(š‘„,š‘”)=š·1š‘¢(š‘”)cos(šœ‹š‘„),šœ“=āˆ’š“1š·1š‘¢(š‘”)cos(šœ‹š‘„),(11) whereš·1=šœ‹š“2šœ‚ī€·š“2āˆ’šœ‹2ī€ø.(12)

2.2.2. Clamped End Conditions

In the case of the clamped-clamped beam, the boundary conditions are as follows:(š‘¤=0,šœƒ=0,šœ“=0)atš‘„=0,1.(13) The admissible function for the š‘¤ is taken as follows:š‘¤=š‘¢(š‘”)(1āˆ’cos(2šœ‹š‘„)).(14) From (7) and (14), the coupled displacement fields for šœƒ and šœ“ are obtained as follows:šœƒ=š·2š‘¢(š‘”)sin(2šœ‹š‘„),šœ“=āˆ’š“1š·2š‘¢(š‘”)sin(2šœ‹š‘„),(15) where:š·2=2šœ‹š“2šœ‚ī€·š“2āˆ’4šœ‹2ī€ø.(16)

2.3. Large Amplitude Vibrations

Consider that there is no external load on the beam, that is, š¹š‘§=0. The total energy for free oscillation of a system without damping is constant (energy conservation principle):š‘‡+š‘ˆ+š‘Š1=Constant,(17) where š‘Š1 is the work done by the tension developed in the LCB due to the large amplitudes and is given by [11, 12]š‘Š1=EIš‘Ÿ28šæ3ī‚µī€œ10š‘¤2,š‘„ī‚¶š‘‘š‘„2.(18) Substituting the expressions for š‘¤,šœƒ, and šœ“, the expressions for š‘ˆ,š‘‡, and š‘Š1 can be obtained asš‘‡=š¶1Ģ‡š‘¢2ī€·š¶,š‘ˆ=2āˆ’š¶3š‘ƒī€øš‘¢2,š‘Š1=š¶4š‘¢4(19) in which coefficients š¶š‘–(š‘–=1āˆ’4) for simple and clamped end conditions are expressed as

simple end conditions: š¶1=š¼1šæ5ī€·š¼1š‘Ÿ2+š¼š‘ƒš“12š·12+š¼2š·12ī€ø4š‘š·11,š¶2=šœ‹2ī€·EIš·12āˆ’2š¶š·12š“1ī€ø+ī€·4šækAGšæ2š·12āˆ’2kAGš‘Ÿšœ‹šæš·1+kAGš‘Ÿ2šœ‹2+GJš“12š·12šœ‹2ī€ø,š¶4šæ3=šœ‹2ī€·š¼š‘ƒš“12š·12+š‘Ÿ2š¼1ī€ø4šæš¼1,š¶4=EIš‘Ÿ2šœ‹432šæ3,(20) clamped end conditions:š¶1=š¼1šæ5ī€·3š¼1š‘Ÿ2+š¼š‘ƒš“12š·22+š¼2š·22ī€ø4š‘š·11,š¶2=šœ‹2ī€·EIš·22āˆ’2Cš·22š“1ī€øšæ+ī€·kAGšæ2š·22āˆ’4kAGš‘Ÿšœ‹šæš·2+4kAGš‘Ÿ2šœ‹2+4GJš“12š·22šœ‹2ī€ø,š¶4šæ3=šœ‹2ī€·š¼š‘ƒš“12š·22+š‘Ÿ2š¼1ī€øšæš¼1,š¶4=EIš‘Ÿ2šœ‹42šæ3.(21) Substituting the expressions for š‘‡,š‘ˆ, and š‘Š1 in (17), we get the following equation:š¶1Ģ‡š‘¢2+ī€·š¶2āˆ’š¶3š‘ƒī€øš‘¢2+š¶4š‘¢4=Constant.(22) Differentiating (22), we obtain the nonlinear ordinary differential equations as follows:š¶Ģˆš‘¢+2āˆ’š¶3š‘ƒš¶1š‘¢+2š¶4š¶1š‘¢3=0ā‰”Ģˆš‘¢+š›¼š‘¢+š›½š‘¢3=0.(23) In order to obtain the postbuckling load-deflection relation, one can set all time-derivative terms in (23) equal to zero which yieldsš‘ƒNL=2š¶4š‘¢2+š¶2š¶3.(24) Neglecting the contribution of š‘¢ in (24), the critical buckling load can be determined asš‘ƒcr=š¶2š¶3.(25)

3. Method of Solution: Implementation of the HAM in Beam Vibrations

For convenience of the reader, we present a brief description of the HAM in Appendix B.

Under the transformation šœ=šœ”š‘” and š‘Š(šœ)=š‘¢(š‘”), (23) becomes as follows:šœ”2Ģˆš‘Š+š›¼š‘Š+š›½š‘Š3=0.(26) The zero-order deformation equation can be written as below ī€ŗšœ‘(1āˆ’š‘ž)šæ(šœ;š‘ž)āˆ’š‘Š0ī€»[šœ‘](šœ)=š‘žā„Žā„(šœ)š‘(šœ;š‘ž)(27) in whichš‘[]šœ‘(šœ;š‘ž)=šœ”2šœ•2šœ‘(šœ;š‘ž)šœ•šœ2+š›¼šœ‘(šœ;š‘ž)+š›½šœ‘(šœ;š‘ž)3=0.(28) We chose the following auxiliary linear operator as:šæ[]šœ‘(šœ;š‘ž)=šœ”02ī‚øšœ•2šœ‘(šœ;š‘ž)šœ•šœ2ī‚¹+šœ‘(šœ;š‘ž).(29) We employ Taylor expansion series for šœ‘(š‘”;š‘ž) andšœ”(š‘ž) asšœ‘(šœ;š‘ž)=šœ‘(šœ;0)+āˆžī“š‘š=11šœ•š‘š!š‘ššœ‘(šœ;š‘ž)šœ•š‘žš‘š||||š‘ž=0Ɨš‘žš‘š=š‘Š0(šœ)+āˆžī“š‘š=1š‘Šš‘š(šœ)š‘žš‘š,(30)šœ”(š‘ž)=šœ”0+āˆžī“š‘š=11šœ•š‘š!š‘ššœ”(š‘ž)šœ•š‘žš‘š||||š‘ž=0š‘žš‘š=šœ”0+āˆžī“š‘š=1šœ”š‘šš‘žš‘š.(31) In order to satisfy the initial conditions, the initial guess of š‘Š(šœ) is chosen as follows:š‘Š0(šœ)=š‘Šmaxcos(šœ).(32) In our case, to obtain the first-order approximation, the function of š‘Š1(t) can be expressed as (see Appendix B)šæī€ŗš‘Š1ī€»[šœ™]||(š‘”)=ā„ā„Ž(š‘”)š‘(š‘”;š‘)š‘=0,(33)š‘Š1(0)=0,š‘‘š‘Š1(0)š‘‘š‘”=0.(34) Assuming ā„=āˆ’1,ā„Ž(t)=1 and after substituting (32) in (33), one would getšœ”02ī€·Ģˆš‘Š1+š‘Š1ī€ø=š‘Šmaxī‚€šœ”cos(šœ)023āˆ’š›¼āˆ’4š›½š‘Š2maxī‚āˆ’š›½š‘Š3max4cos(3šœ),(35)š‘Š1Ģ‡š‘Š(0)=0,1(0)=0.(36) Eliminating the secular term, we havešœ”0=ī‚™3š›¼+4š›½š‘Š2max.(37) Solving (35) and (36), the š‘Š1(šœ) is obtained as follows:š‘Š1šæ(šœ)=18šœ”02(cos(šœ)āˆ’cos(3šœ)),(38) wherešæ11=āˆ’4š›½š‘Š3max.(39) Thus the first-order approximation of the š‘Š(šœ) yields to,š‘Š(šœ)=š‘Š0(šœ)+š‘Š1(šœ)(40) in whichšœ=šœ”š‘”,šœ”=šœ”0.(41)

4. Numerical Results and Discussion

To validate the presented relations, the results of the natural frequency for the single isotropic lamina and the LCB with simple and clamped end conditions are compared with the available results in the literatures.

The variables used in the tables and plots are the nonlinear to linear frequency ratio, šœ”NL/šœ”L, the nondimensional amplitude of the vibration, š‘Šmax, and the postbuckling load to critical load, š‘ƒNL/š‘ƒcr, (henceforth referred to as the frequency ratio, the amplitude ratio, and the buckling load ratio, resp.).

The materials steel and AS4/3501 Graphite-Epoxy, having the following mechanical properties, are used for the study.


AS4/3501 graphite/epoxy:šø11=144.8Gpa,šø22šŗ=9.65Gpa,12=4.14Gpa,šŗ13šŗ=4.14Gpa,(43)23=4.14Gpa,šœ12=0.3,šœŒ=1389.23kg/m3.(44)

The shear correction factor š‘˜ is taken as 5/6, as commonly used in the literature. All the layers are of equal thickness. Fiber orientation is measured from š‘„-axis.

The validity of the results is established by Tables 1 and 3 for simple end conditions and Table 2 for clamped end conditions. Table 1 shows the frequency ratio for the single isotropic beam along with the results from literatures. It can be seen that the frequency obtained from two different beam models (i.e., Bernoulli-Euler [1, 13, 14] and Timoshenko [15] beams) yields the same value provided that the beamā€™s slenderness ratio be at least 100. This finding is also reported by [16] in which only the linear analysis was studied. The same results from [17] are included in Table 2 for the sake of comparison, and the agreement of the values of frequency ratio is good for various values of the amplitude and slenderness ratios.

Also, comparison of the frequency ratio for four-layer symmetric cross-ply, angle-ply and general lay-up composite beams with the finite element solutions from [7] for different amplitude ratio is reported in Table 3. As can be seen from these tables, the agreement between the results is quite good.

Figure 2 illustrates the effects of shear deformation and rotary inertia on the frequency ratio for three different amplitude ratios. It is clear that the CLT predicts the frequency ratio lower than the one obtained by FSDT. For thin beams and based on the two considered beam models, namely, CLT and FSDT, almost no difference is seen for different slenderness ratio.

The same layers, that is, layers oriented at 0, 90, 90, and 0, in two different configurations are considered in Figure 3 to examine the stacking sequence on the frequency ratio. As can be observed, the frequency ratio for [0/90/90/0] configuration is higher than the [90/0/0/90] configuration. Furthermore, as the amplitude ratio increases, the frequency ratio increases for two considered lamination schemes.

Figure 4 displays the frequency ratio versus the angle of orientation for the symmetric angle-ply (šœƒ/āˆ’šœƒ/āˆ’šœƒ/šœƒ) thin and moderately thick LCB and for three different amplitude ratios and two considered boundary conditions. It is evident that all figures for thin and thick beams reach their maximum values at the orientation around (15Ā°ā€“20Ā°). This result is also reported by [18] in which only the linear study was conducted. In the intermediate region (60āˆ˜ā‰¤šœƒā‰¤90āˆ˜), the frequency ratio remains almost unchanged.

Variation of the buckling load ratio versus amplitude ratio is shown in Figure 5. As can be observed for thick beams, there are significant differences between buckling load ratio predicted by CLT and FSDT. This difference is more considerable at larger amplitudes.

5. Conclusions

Analytical expressions for the frequency and deflection equations of the LCB at large amplitude on the basis of FSDT with simple end conditions have been derived. The applicability of the theory is demonstrated by numerical results, which show good agreement with published results.

Limited numerical studies are conducted to examine the effect of slenderness ratio, fiber orientation, and amplitude of vibration on the vibration and buckling characteristics of the LCB. The present study can serve as a quick and accurate reference to predict the vibration and postbuckling response of the composite beam at any amplitude.


A. Derivation of Bending (EI), Torsional (GJ), Bending-Torsion Coupling (C), and Shear (kAG) Rigidities

Based on CLT, the constitutive equations of the laminate relating the stress resultant and the curvatures can be written asāŽ§āŽŖāŽØāŽŖāŽ©š‘€š‘„š‘€š‘¦š‘€š‘„š‘¦āŽ«āŽŖāŽ¬āŽŖāŽ­=āŽ”āŽ¢āŽ¢āŽ¢āŽ£š·11š·12š·16š·12š·22š·26š·16š·26š·66āŽ¤āŽ„āŽ„āŽ„āŽ¦āŽ§āŽŖāŽØāŽŖāŽ©š‘˜š‘„š‘˜š‘¦š‘˜š‘„š‘¦āŽ«āŽŖāŽ¬āŽŖāŽ­,(A.1) where š·š‘–š‘— are anisotropic stiffness coefficients which are defined byš·š‘–š‘—=ī€œš‘”/2āˆ’š‘”/2š‘„š‘–š‘—š‘§2š‘‘š‘§,(š‘–,š‘—=1,2,6)(A.2) where š‘„š‘–š‘— and š‘” are the transformed material constants [19] and thickness of the beam. The moment š‘€š‘¦ can be considered to be zero. Then š‘˜š‘¦ can be eliminated from (A.1) to giveīƒÆš‘€š‘„š‘€š‘„š‘¦īƒ°=āŽ”āŽ¢āŽ¢āŽ¢āŽ¢āŽ¢āŽ£š·11āˆ’š·212š·22š·16āˆ’š·12š·26š·22š·16āˆ’š·12š·26š·22š·66āˆ’š·226š·22āŽ¤āŽ„āŽ„āŽ„āŽ„āŽ„āŽ¦īƒÆš‘˜š‘„š‘˜š‘„š‘¦īƒ°.(A.3) The bending š‘€š‘„ and twisting š‘€š‘„š‘¦ moments can be related to the resultant bending (š‘€) and torsional (š‘‡) moments as follows [20]:š‘€=āˆ’š‘š‘€š‘„,š‘‡=2š‘š‘€š‘„š‘¦.(A.4) Also, the curvatures š‘˜š‘„ and š‘˜š‘„š‘¦ can be related to the bending slope šœƒī…ž and twist rate šœ“ī…ž so that [20]š‘˜š‘„=āˆ’šœƒī…ž,š‘˜š‘„š‘¦=2šœ“ī…ž.(A.5) Thus, the expressions for š‘€ and š‘‡ are given byīƒÆš‘€š‘‡īƒ°āŽ”āŽ¢āŽ¢āŽ¢āŽ¢āŽ¢āŽ£š·=š‘11āˆ’š·212š·222ī‚µš·16āˆ’š·12š·26š·22ī‚¶2ī‚µš·16āˆ’š·12š·26š·22ī‚¶4ī‚µš·66āˆ’š·226š·22ī‚¶āŽ¤āŽ„āŽ„āŽ„āŽ„āŽ„āŽ¦īƒÆšœƒī…žšœ“ī…žīƒ°.(A.6) The common relationship between the bending š‘€ and torsional š‘‡ moments with the bending slope šœƒī…ž and twist rate šœ“ī…ž are as follows:īƒÆš‘€š‘‡īƒ°=īƒ¬šœƒEIš¶š¶GJīƒ­īƒÆī…žšœ“ī…žīƒ°.(A.7) An element-by-element comparison of (A.6) and (A.7), we obtainī‚µš·EI=š‘11āˆ’š·212š·22ī‚¶ī‚µš·,š¶=2š‘16āˆ’š·12š·26š·22ī‚¶,ī‚µš·GJ=4š‘66āˆ’š·226š·22ī‚¶.(A.8) The transverse shear force-strain relation for the LCB can be expressed asš‘„š‘„š‘§ī‚µī€œ=kbš‘”/2āˆ’š‘”/2š‘„55ī‚¶š›¾š‘‘Ģ‚š‘§š‘„š‘§=kAGš›¾š‘„š‘§,(A.9) where:ī€œkAG=kbš‘”/2āˆ’š‘”/2š‘„55š‘‘Ģ‚š‘§.(A.10)

B. An Overview of Homotopy Analysis Method (HAM)

For convenience of the reader, we will first present a brief description of the HAM [21]. Consider the following nonlinear homogeneous differential equations:š‘[]š‘¢(š‘”)=0,(B.1) where š‘ is nonlinear operators, š‘” denotes the independent variable, and š‘¢(š‘”)=0 are unknown functions.

Liao constructed the so-called zero-order deformation equation as [21]ī€ŗšœ™(1āˆ’š‘)šæ(š‘”;š‘)āˆ’š‘¢0ī€»[šœ™](š‘”)=š‘ā„ā„Ž(š‘”)š‘(š‘”;š‘),(B.2) where š‘āˆˆ[0,1] is an embedding parameter, ā„ are nonzero auxiliary functions, ā„Ž(t) are nonzero auxiliary function, šæ is an auxiliary linear operator, š‘¢0(š‘”) are initial guesses of š‘¢(t), and šœ™(t;p) are unknown functions. As š‘ increases from 0 to 1, the šœ™(š‘”;š‘) varies from the initial approximation to the exact solution. In other words,šœ™(š‘”;0)=š‘¢0(š‘”) andšœ™(š‘”;1)=š‘¢(š‘”).

The deformation equation (B.2) can provide us with a family of approximation series whose convergence region depends upon the auxiliary parameter ā„ and the auxiliary function ā„Ž(t). More importantly, this provides us with a simple way to adjust and control the convergence regions and rates of approximation series.

Differentiating once more from (B.2) with respect to the embedding parameter š‘ and then setting š‘=0, we obtain the first-order deformation equation asšæī€ŗš‘¢1ī€»[šœ™]=ā„ā„Ž(š‘”)š‘(š‘”;š‘)š‘=0(B.3) which gives the first-order approximation for the š‘¢(š‘”) [21]. The higher order approximation of the solution can be obtained by calculating the š‘š-order (š‘š>1) deformation equation which can be calculated by differentiating (B.3) š‘š times with respect to the š‘ and then setting š‘=0 and finally dividing them by š‘š! [21]. Therefore, the š‘šth-order deformation equation can be obtained as follows:šæī€ŗš‘¢š‘š(š‘”)āˆ’š‘¢š‘šāˆ’1ī€»(š‘”)=ā„ā„Ž(š‘”)š‘…š‘šī€·īƒ‰š‘¢š‘šāˆ’1,ā€¦,īƒ‰š‘¢š‘šāˆ’1ī€ø,(B.4) where (īƒ‰š‘¢š‘šāˆ’1,ā€¦,īƒ‰š‘¢š‘šāˆ’1) and š‘…š‘š(īƒ‰š‘¢š‘šāˆ’1,ā€¦,īƒ‰š‘¢š‘šāˆ’1) are defined asš‘…š‘šī€·īƒ‰š‘¢š‘šāˆ’1,ā€¦,īƒ‰š‘¢š‘šāˆ’1ī€ø=1š‘‘(š‘šāˆ’1)!š‘šāˆ’1š‘[]šœ™(š‘”;š‘)š‘‘š‘š‘šāˆ’1||||š‘ž=0īƒ‰š‘¢š‘šāˆ’1=ī€½š‘¢0,š‘¢1,š‘¢2,ā€¦,š‘¢š‘šāˆ’1ī€¾,(B.5) subjected to the following initial conditions:š‘¢š‘š(0)=0,Ģ‡š‘¢š‘š(0)=0.(B.6) To solve (B.2), we employ Taylor expansion series for šœ™(š‘”;š‘) asšœ™(š‘”;š‘)=šœ™(š‘”;0)+āˆžī“š‘š=11ī‚µšœ•š‘š!š‘ššœ™(š‘”;š‘)šœ•š‘š‘šī‚¶||||š‘=0š‘š‘šāŸ¹š‘¢(š‘”)=š‘¢0(š‘”)+āˆžī“š‘š=1š‘¢š‘š(š‘”)š‘š‘š,(B.7) where š‘¢š‘š(š‘”) is called the š‘šth-order derivative of unknown function šœ™(š‘”;š‘).


The authors are grateful for the encouragement and useful comments of the reviewers.