Abstract

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.

Steel:𝐸=206.8Gpa,𝐺=79Gpa,𝜐=0.3,𝜌=7800kg/m3.(42)

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.

Appendices

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 πœ™(𝑑;𝑝).

Acknowledgment

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