An Analytical Study on Thermal Post-Buckling Behaviors of Geometrically Imperfect FRC Laminated Beams Using a Modified Zig-Zag Beam Model

1. Introduction

As a kind of light-weight and high-strength material, fiber-reinforced composites (FRC) are increasingly used in military aircrafts where weight and performance requirements are stringent [1,2,3]. However, the extreme service environment of military aircrafts puts forward higher requirements for the structure safety analysis and design techniques, that is, to be able to predict the nonlinear response of structures under given working conditions, analyze the influence of various factors on the response results and make reliable estimates of the accuracy of prediction results [4,5]. For slender composite structures with immovable boundaries serving in aerospace, both the initial geometric imperfections and the axial compressive internal forces caused by the thermal environment will reduce the axial stiffness and cause structural buckling instability，which requires quantitative and precise analysis. Nowadays, researchers not only focus on how to avoid buckling, but also how to make full use of the bearing capacity of structures in post-buckling due to the requirement of lightweight structure design [6,7]. To achieve this, the accurate prediction of buckling and the accurate analysis of stability and defect sensitivity of the post-buckling equilibrium path is the key point.

Researches on post-buckling behaviors of geometrically imperfect beams have received a good amount of attention in the literatures [8,9,10,11,12] that reveal the important influence of geometric imperfections, support conditions and other factors on structural buckling and post-buckling. Based on Reddy's higher order shear deformation beam theory, the post-buckling behaviors were analyzed of simply supported laminated beams with geometric imperfections considering the temperature-dependent material properties and uniform temperature distribution through the thickness [8,13]. Yaghoobi and Torabi conducted the post-buckling and nonlinear free vibration analysis of geometrically imperfect functionally graded beams resting on nonlinear elastic foundation using the Galerkin’s method and variational iteration method [14]. Wu et al. studied the sensitivity of the post-buckling behaviors of functionally graded carbon nanotube-reinforced composite beams to different types of initial generic geometric imperfections by differential quadrature method and an iteration procedure [9]. Numerical results suggested the post-buckling response is highly sensitive to the imperfection amplitude, especially for the L1-mode imperfection. Mohammadi et al. investigated the effects of nonlocal parameter, geometrical imperfection, and elastic foundation on the static instability of the nonlinear nanobeam in post-buckling, snap through and bifurcation instability based on a size-dependent Euler–Bernoulli beam [15]. Fan and Wang studied the effects of initial thermal bending moment, geometric imperfections and matrix cracking on thermal post-buckling behaviors of beam with carbon nanotube reinforced composite layers and piezoelectric fiber reinforced composite layers based on a higher-order beam theory [16].

Because the accuracy of the results depends largely on the precision of the theory, many scholars have been devoting themselves to the development of the non-linear theories of beams, plates and shells [17,18,19,20]. For decades, the beam models like Euler-Bernoulli beam, Timoshenko beam and higher-order shear deformation beam have been well developed and successfully used for the researches on the static and dynamic characteristics of single-layer homogeneous isotropic beams. However, there will exist some limitations in accuracy if applying these theories to the analysis of laminated structures because these theories do not accurately describe the true displacement distribution of laminated structures that should be continuous in its thickness direction due to the difference of stiffness of the adjacent laminates. In the current terminology, the slope change at the interlayer interface is termed as the zig-zag effect. Researches suggested for laminated structures with large thickness or transverse shear deformation, the free vibration frequencies of structures obtained by the equivalent single-layer theory neglecting the zig-zag effects are quite different from those obtained by the three-dimensional elastic theory [21]. To account for such continuous through-thickness displacement field, a series of layer-wise theories have been proposed which separately assumes displacement variables for each layer and may trace the local variations in each layer more accurately [22]. However, the improvement of accuracy is at the expense of the large number of unknowns proportional to the number of layers. For composite laminated structures with more layers, too many displacement variables make the form of solutions too complex, and the corresponding finite element model is also lack of practical value because of the large number of the nodal DOF. To be refined without increasing calculation costs, the equivalent single-layer zig-zag models have been developed to achieve the accuracy of layer-wise models with a lower number of variables. By adding a zig-zag function to the displacement field of the equivalent single-layer models, Murakami pioneered the development of zig-zag theories in which the continuity condition of displacement could be satisfied and the number of displacement variables could be independent of the number of layers [23]. Carrera modified the zig-zag function to meet the condition that the shear stress component on the upper and lower surfaces of the laminated beam is zero [24]. Numerical results revealed that the multilayered plate and shell theories can be greatly improved by use of this modified zig-zag function. Based on this work, Xie et al. developed a general higher-order shear deformation zig-zag theory for analyzing the aero thermoelastic characteristics of composite laminated panels under supersonic airflows [25]. Inspired by Qu and Carrera E’s work, this study applies the modified higher-order shear deformation zig-zag beam model to the post-buckling analysis of geometrically imperfect FRC laminated beams for the first time, and analyze the accuracy and applicability of the theory considering the zig-zag effect and not considering the zig-zag effect.

A number of analytical and numerical methods have been developed to study the static and dynamic responses of laminated beams. They include, but are not limited to, the variational method [26], Navier type method [27], Galerkin method [28], multiple scales methods [29], Green function method [30], harmonic balance method [31], differential quadrature method [9], finite element method [32], meshless method [33], and perturbation method [34]. For structures modeled by such a refined theory considering the zig-zag effects at the interlayer interface, finite element method is the general method to obtain the numerical solutions [35,36]. For the first time, a two-step perturbation method that has better computation efficiency and is more convenient for parameter analysis is used to obtain the asymptotic analytic solution to the post-buckling behaviors of slender structures modeled by a refined zig-zag theory where the displacement field is represented by a piecewise function.

To the best of the author's knowledge, almost all the researches on the post-buckling analyses of slender structures with geometric imperfections are limited to beams with hinged-hinged ends. In this work, the buckling and post-buckling behaviors of perfect/imperfect hinged-hinged and clamped-clamped FRC laminated beams are systematically studied for the first time. The solution results based on the zig-zag theory and the equivalent single-layer theories are verified and compared in detail, and the applicability of the two kinds of theories are analyzed. Through detailed parameter study, various factors affecting the buckling and post-buckling behaviors of the beam are intuitively presented, and the existence of the equilibrium path of bifurcation type is discussed in detail. It is believed that the current zig-zag models and the used solution methods can be extended to the study on the nonlinear static and dynamic characteristics of laminated or sandwich structures with more precise results.

2. Modeling of FRC Laminated Beams

Consider a laminated composite beam with length L, thickness h, and unit width, consisting of N plies of any kind as shown in Figure 1. The two ends are assumed to be hinged-hinged or clamped-clamped and the in-plane boundary conditions supposed immovable. Considering the large transverse shear deformation of beam and the discontinuity of shear strain at interlayer interface, a zig-zag beam model based on Reddy’s higher-order shear deformation theory is adopted, in which the axial and transverse displacement fields are expressed as

$$\tilde{u}\left(x,z\right)=u\left(x\right)+f\left(z\right){\displaystyle \frac{\partial w\left(x\right)}{\partial x}}+g\left(z\right)\vartheta \left(x\right)+\phi \left(z,k\right)\eta \left(x\right)$$

(1a)

$$\tilde{w}\left(x,z\right)=w\left(x\right)+w^{\ast }\left(x\right)$$

(1b)

where $\tilde{u}$ and $\tilde{w}$ are the displacements in x- and z-directions, respectively, at any material point in the (x,z) plane. $u,w$ represents the axial and transverse displacements of the mid-plane in geometry, with the transverse initial geometric imperfection denoted by $w^{\ast }$. z is the depth of the material point measured from the mid-plane along the positive z-axis. ${\displaystyle \frac{\partial w}{\partial x}}$ is the rotation angle of the mid-plane, $\vartheta$ is the angle between the tangent of the warped section at the mid-plane and the positive z-axis and $\eta$ is a higher-order generalized displacement. $f,g$ are the higher-order shape functions, defined as [17]

$$f(z)=-{\displaystyle \frac{4z^3}{3h^2}},g(z)=z-{\displaystyle \frac{4z^3}{3h^2}}$$

(2)

and $\phi$ is a zig-zag function that has been modified on the basis of Murakami’s zig-zag function to meet the stress boundary conditions [25], that is

$$\phi ={\left(-1\right)}^k\left\{{\displaystyle \frac{2}{h_k}}\left[z-{\displaystyle \frac{1}{2}}\left(z_{k+1}+z_k\right)\right]-{\displaystyle \frac{8z^3}{3h_k{h}^2}}\right\}$$

(3)

in which $z_k$ and $z_{k+1}$ represent the coordinates of the bottom and top of the kth laminate in z direction while $h_k$ and h represent the thickness of the kth laminate and the whole beam, respectively.

According to Marguerre’s nonlinear theory of structures with initial deflection [36], the geometric equations considering the initial imperfection can be derived

$${\epsilon }_{xx}={\displaystyle \frac{\partial \tilde{u}}{\partial x}}+{\displaystyle \frac{1}{2}}{\left({\displaystyle \frac{\partial w}{\partial x}}\right)}^2+{\displaystyle \frac{\partial w}{\partial x}}{\displaystyle \frac{\partial w^{\ast }}{\partial x}}={\epsilon }_{xx}^{\left(0\right)}+f{\epsilon }_{xx}^{\left(1\right)}+g{\epsilon }_{xx}^{\left(2\right)}+\phi {\epsilon }_{xx}^{\left(3\right)}$$

(4a)

$${\gamma }_{zx}={\displaystyle \frac{\partial \tilde{u}}{\partial z}}+{\displaystyle \frac{\partial w}{\partial x}}=\overline{f}{\gamma }_{zx}^{(0)}+\overline{g}{\gamma }_{zx}^{(1)}+\overline{\phi }{\gamma }_{zx}^{(2)}$$

(4b)

in which ${\epsilon }_{xx}^{(i)}$(i=0,1,2,3) are generalized in-plane positive strains, ${\gamma }_{zx}^{(i)}$(i=0,1,2) are generalized shear strains, and $\overline{f}$、$\overline{g}$ and $\overline{\phi }$ are defined as

$$\overline{f}={\displaystyle \frac{\partial f}{\partial z}}+1,\overline{g}={\displaystyle \frac{\partial g}{\partial z}}, \overline{\phi }={\displaystyle \frac{\partial \phi }{\partial z}}$$

(5)

The generalized strain components defined in Eq. (4) can be expressed as

$${\epsilon }_{xx}^{(0)}={\displaystyle \frac{\partial u}{\partial x}}+{\displaystyle \frac{1}{2}}{\left({\displaystyle \frac{\partial w}{\partial x}}\right)}^2+{\displaystyle \frac{\partial w}{\partial x}}{\displaystyle \frac{\partial w^{\ast }\left(x\right)}{\partial x}}, {\epsilon }_{xx}^{(1)}={\displaystyle \frac{{\partial }^{2}w}{\partial x^2}}, {\epsilon }_{xx}^{(2)}={\displaystyle \frac{\partial \vartheta }{\partial x}}, {\epsilon }_{xx}^{(3)}={\displaystyle \frac{\partial \eta }{\partial x}}$$

(6)

$${\gamma }_{zx}^{(0)}={\displaystyle \frac{\partial w}{\partial x}}, {\gamma }_{zx}^{(1)}=\vartheta , {\gamma }_{zx}^{(2)}=\eta$$

(7)

Based on the in-plane bending assumption, stresses in the y direction can be ignored and the constitutive model of the kth layer of the beam under thermal environment can be written as [37]

$${\left[\begin{array}{c}{\sigma }_{xx}\\ 0\\ 0\end{array}\right]}^{(k)}=\left[\begin{array}{ccc}{\overline{Q}}_{11}^{(k)}& {\overline{Q}}_{12}^{(k)}& {\overline{Q}}_{16}^{(k)}\\ {\overline{Q}}_{12}^{(k)}& {\overline{Q}}_{22}^{(k)}& {\overline{Q}}_{26}^{(k)}\\ {\overline{Q}}_{16}^{(k)}& {\overline{Q}}_{26}^{(k)}& {\overline{Q}}_{66}^{(k)}\end{array}\right]\left\{{\left[\begin{array}{c}{\epsilon }_{xx}\\ {\epsilon }_{yy}\\ {\epsilon }_{xy}\end{array}\right]}^{(k)}-\Delta T{\left[\begin{array}{c}{\alpha }_{xx}\\ {\alpha }_{yy}\\ {\alpha }_{xy}\end{array}\right]}^{(k)}\right\}$$

(8a)

$${\left[\begin{array}{c}0\\ {\sigma }_{zx}\end{array}\right]}^{(k)}=\left[\begin{array}{cc}{\overline{Q}}_{44}^{(k)}& {\overline{Q}}_{45}^{(k)}\\ {\overline{Q}}_{45}^{(k)}& {\overline{Q}}_{55}^{(k)}\end{array}\right]{\left[\begin{array}{c}{\epsilon }_{yz}\\ {\epsilon }_{zx}\end{array}\right]}^{(k)}$$

(8b)

where the transformed thermal coefficients ${\alpha }_{xx}$, ${\alpha }_{yy}$ and ${\alpha }_{xy}$ in Eq. (8a) are defined as

$$\left[\begin{array}{c}{\alpha }_{xx}\\ {\alpha }_{yy}\\ {\alpha }_{xy}\end{array}\right]=\left[\begin{array}{cc}{m}^2& n^2\\ n^2& m^2\\ 2mn& -2mn\end{array}\right]\left[\begin{array}{c}{\alpha }_1\\ {\alpha }_2\end{array}\right]$$

(9)

and ${\overline{Q}}_{ij}^{(k)}$ are the transformed stiffness coefficients, defined by [38]

$${\left\{\begin{array}{c}\begin{array}{c}\begin{array}{c}{\overline{Q}}_{11}\\ {\overline{Q}}_{22}\end{array}\\ {\overline{Q}}_{12}\end{array}\\ {\overline{Q}}_{66}\\ {\overline{Q}}_{16}\\ {\overline{Q}}_{26}\end{array}\right\}}^{\left(k\right)}={\left[\begin{array}{cccc}\begin{array}{c}\begin{array}{c}{m}^4\\ n^4\end{array}\\ m^2{n}^2\end{array}& \begin{array}{c}{n}^4\\ m^4\\ m^2{n}^2\end{array}& \begin{array}{c}2m^2{n}^2\\ 2m^2{n}^2\\ m^4+n^4\end{array}& \begin{array}{c}4m^2{n}^2\\ 4m^2{n}^2\\ -4m^2{n}^2\end{array}\\ m^2{n}^2& m^2{n}^2& -2m^2{n}^2& {\left(m^2-n^2\right)}^2\\ m^{3}n& -m{n}^3& m{n}^3-m^{3}n& 2\left(m{n}^3-m^{3}n\right)\\ m{n}^3& -m^{3}n& m^{3}n-m{n}^3& 2\left(m^{3}n-m{n}^3\right)\end{array}\right]}^{\left(k\right)}{\left\{\begin{array}{c}{Q}_{11}\\ Q_{22}\\ Q_{12}\\ Q_{66}\end{array}\right\}}^{\left(k\right)}$$

(10a)

$${\left\{\begin{array}{c}{\overline{Q}}_{44}\\ {\overline{Q}}_{45}\\ {\overline{Q}}_{55}\end{array}\right\}}^{\left(k\right)}={\left[\begin{array}{cc}\begin{array}{c}{m}^2\\ -mn\end{array}& \begin{array}{c}{n}^2\\ mn\end{array}\\ n^2& n^2\end{array}\right]}^{\left(k\right)}{\left\{\begin{array}{c}{Q}_{44}\\ Q_{55}\end{array}\right\}}^{\left(k\right)}$$

(10b)

in which

$$Q_{11}^{\left(k\right)}={\displaystyle \frac{E_1^{\left(k\right)}}{1-{\mu }_{12}^{\left(k\right)}{\mu }_{21}^{\left(k\right)}}},Q_{12}^{\left(k\right)}={\displaystyle \frac{{\mu }_{12}^{\left(k\right)}{E}_2^{\left(k\right)}}{1-{\mu }_{12}^{\left(k\right)}{\mu }_{21}^{\left(k\right)}}},Q_{22}^{\left(k\right)}={\displaystyle \frac{E_2^{\left(k\right)}}{1-{\mu }_{12}^{\left(k\right)}{\mu }_{21}^{\left(k\right)}}}$$

(11a)

$$Q_{66}^{\left(k\right)}=G_{12},Q_{44}^{\left(k\right)}=G_{23},Q_{55}^{\left(k\right)}=G_{13}$$

(11b)

and $m=\cos \theta , n=\sin \theta$, where $\theta$ is the lamination angle with respect to the beam axial direction.

On the basis of Eqs. (8), the axial normal stress and transverse shear stress at any point in the k-layer of the beam can be further written, according to the model in [39]

$${{\sigma }_{xx}}^{(k)}={\tilde{Q}}_{11}^{\left(k\right)}{{\epsilon }_{xx}}^{(k)}-\left[\begin{array}{ccc}{\overline{Q}}_{11}^{(k)}& {\overline{Q}}_{12}^{(k)}& {\overline{Q}}_{16}^{(k)}\end{array}\right]\left[\begin{array}{c}{\alpha }_{xx}\\ {\alpha }_{yy}\\ {\alpha }_{xy}\end{array}\right]\Delta T$$

(12a)

$${{\tau }_{zx}}^{(k)}={\tilde{Q}}_{55}^{\left(k\right)}{{\gamma }_{zx}}^{(k)}$$

(12b)

where

$${ \mathrm{Q} ˜}_{11}^{\left(\mathrm{k}\right)}={ \mathrm{Q} ¯}_{11}^{\left(\mathrm{k}\right)}-{\left[\begin{array}{c}{ \mathrm{Q} ¯}_{12}^{\left(\mathrm{k}\right)}\\ { \mathrm{Q} ¯}_{16}^{\left(\mathrm{k}\right)}\end{array}\right]}^{\mathrm{T}}{\left[\begin{array}{cc}{ \mathrm{Q} ¯}_{22}^{\left(\mathrm{k}\right)}& { \mathrm{Q} ¯}_{26}^{\left(\mathrm{k}\right)}\\ \mathrm{sym}& { \mathrm{Q} ¯}_{66}^{\left(\mathrm{k}\right)}\end{array}\right]}^{-1}\left[\begin{array}{c}{ \mathrm{Q} ¯}_{12}^{\left(\mathrm{k}\right)}\\ { \mathrm{Q} ¯}_{16}^{\left(\mathrm{k}\right)}\end{array}\right]$$

(13a)

$${ \mathrm{Q} ˜}_{55}^{\left(\mathrm{k}\right)}=\left[{ \mathrm{Q} ¯}_{55}^{\left(\mathrm{k}\right)}-{\displaystyle \frac{1}{{ \mathrm{Q} ¯}_{44}^{\left(\mathrm{k}\right)}}}{\left({ \mathrm{Q} ¯}_{45}^{\left(\mathrm{k}\right)}\right)}^2\right]$$

(13b)

By integrating stresses in thickness direction, the generalized internal forces and moments are expressed as

$$\left[\begin{array}{c}{N}_x\\ M_x\\ P_x\\ T_x\end{array}\right]=\left[\begin{array}{cccc}{A}_{11}& B_{11}& E_{11}& F_{11}\\ & D_{11}& G_{11}& H_{11}\\ & & I_{11}& J_{11}\\ sym& & & R_{11}\end{array}\right]\left[\begin{array}{c}{\epsilon }_{xx}^{\left(0\right)}\\ {\epsilon }_{xx}^{\left(1\right)}\\ {\epsilon }_{xx}^{\left(2\right)}\\ {\epsilon }_{xx}^{\left(3\right)}\end{array}\right]-\left[\begin{array}{c}{A}_{11}^T\\ B_{11}^T\\ E_{11}^T\\ F_{11}^T\end{array}\right]\Delta T$$

(14a)

$$\left[\begin{array}{c}{Q}_{zx}\\ P_{zx}\\ T_{zx}\end{array}\right]=\left[\begin{array}{ccc}{D}_{55}& G_{55}& H_{55}\\ & I_{55}& J_{55}\\ sym& & R_{55}\end{array}\right]\left[\begin{array}{c}{\gamma }_{zx}^{\left(0\right)}\\ {\gamma }_{zx}^{\left(1\right)}\\ {\gamma }_{zx}^{\left(2\right)}\end{array}\right]$$

(14b)

where $N_x$ and $Q_{zx}$ are axial positive and transverse shear forces, $M_x$ is bending moment per unit length, $P_x$ and $T_x$ are higher-order bending moments, and $P_{zx}$ and $T_{zx}$ are higher-order shear forces, respectively. The beam reduced stiffness coefficients are expressed by

$$\begin{array}{l}\left(A_{11},B_{11},D_{11},E_{11},F_{11},G_{11},H_{11},I_{11},J_{11},R_{11}\right)\\ ={\displaystyle \sum }_{k=1}^{N_l}\underset{z_k}{\overset{z_{k+1}}{{\displaystyle \int }}}\left(1,f,f^2,g,\phi ,fg,f\phi ,g^2,g\phi ,{\phi }^2\right){\tilde{Q}}_{11}^{\left(k\right)}dz\end{array}$$

(15a)

$$(D_{55},G_{55},H_{55},I_{55},J_{55},R_{55})={\displaystyle \sum }_{k=1}^{N_l}\underset{z_k}{\overset{z_{k+1}}{{\displaystyle \int }}}\left({\overline{f}}^2,\overline{f}\overline{g},\overline{f}\overline{\phi },{\overline{g}}^2,\overline{g}\overline{\phi },{\overline{\phi }}^2\right){\tilde{Q}}_{55}^{\left(k\right)}dz$$

(15b)

and

$$\left(A_{11}^T,B_{11}^T,E_{11}^T,F_{11}^T\right)={\displaystyle \sum }_{k=1}^{N_l}\underset{z_k}{\overset{z_{k+1}}{{\displaystyle \int }}}\left(1,f,g,\phi \right){\tilde{\alpha }}_x^{\left(k\right)}dz$$

(16)

in which

$${\tilde{\alpha }}_x^{\left(k\right)}=\left[\begin{array}{ccc}{\overline{Q}}_{11}^{\left(k\right)}& {\overline{Q}}_{12}^{\left(k\right)}& {\overline{Q}}_{16}^{\left(k\right)}\end{array}\right]\cdot \left[\begin{array}{c}{\alpha }_x\\ {\alpha }_y\\ {\alpha }_{xy}\end{array}\right]$$

(17)

The virtual work done by the external force is ignored and the Hamilton energy-variational principle for such a beam is

$${\displaystyle {\int }_{t_1}^{t_2}\left(\delta T-\delta U\right)dt=}0$$

(18)

where $\delta U$ is the virtual strain energy, $\delta T$ is the virtual kinetic energy, and

$$\begin{array}{l}U={\displaystyle \frac{1}{2}}{\displaystyle \sum }_{k=1}^{N_l}\underset{z_k}{\overset{z_{k+1}}{{\displaystyle \int }}}\underset{0}{\overset{L}{{\displaystyle \int }}}\left({{\sigma }_{xx}}^{\left(k\right)}{{\epsilon }_{xx}}^{\left(k\right)}+{{\tau }_{xz}}^{\left(k\right)}{{\gamma }_{xz}}^{\left(k\right)}\right)dxdz\\ ={\displaystyle \frac{1}{2}}\underset{0}{\overset{L}{{\displaystyle \int }}}\left(N_x{\epsilon }_{xx}^{\left(0\right)}+M_x{\epsilon }_{xx}^{\left(1\right)}+P_x{\epsilon }_{xx}^{\left(2\right)}+T_x{\epsilon }_{xx}^{\left(3\right)}+Q_{xz}{\gamma }_{xz}^{\left(0\right)}+P_{xz}{\gamma }_{xz}^{\left(1\right)}+T_{xz}{\gamma }_{xz}^{\left(2\right)}\right)dx\end{array}$$

(19a)

$$T={\displaystyle \frac{1}{2}}{\displaystyle \sum }_{k=1}^{N_l}\underset{z_k}{\overset{z_{k+1}}{{\displaystyle \int }}}\underset{0}{\overset{L}{{\displaystyle \int }}}{\rho }^{\left(k\right)}\left[{\left(\dot{u}+f{\displaystyle \frac{\partial \dot{w}}{\partial x}}+g\dot{\vartheta }+\phi \dot{\eta }\right)}^2+{\dot{w}}^2\right]dxdz$$

(19b)

Through variational operations, the following equations of motion can be obtained

$$\begin{array}{l}\delta u:{\displaystyle \frac{\partial N_x}{\partial x}}={\rho }_0\ddot{u}+{\rho }_1{\displaystyle \frac{\partial \ddot{w}}{\partial x}}+{\rho }_2\ddot{\vartheta }+{\rho }_3\ddot{\eta }\hfill \\ \delta w:{\displaystyle \frac{\partial Q_{xz}}{\partial x}}-{\displaystyle \frac{{\partial }^2{M}_x}{\partial x^2}}+N_x{\displaystyle \frac{{\partial }^2\left(w+w^{\ast }\right)}{\partial x^2}}={\rho }_0\ddot{w}-{\displaystyle \frac{\partial }{\partial x}}\left({\rho }_1\ddot{u}+{\rho }_4{\displaystyle \frac{\partial \ddot{w}}{\partial x}}+{\rho }_5\ddot{\vartheta }+{\rho }_6\ddot{\eta }\right)\hfill \\ \delta \vartheta :{\displaystyle \frac{\partial P_x}{\partial x}}-P_{xz}={\rho }_2\ddot{u}+{\rho }_5{\displaystyle \frac{\partial \ddot{w}}{\partial x}}+{\rho }_7\ddot{\vartheta }+{\rho }_8\ddot{\eta }\hfill \\ \delta \eta :{\displaystyle \frac{\partial T_x}{\partial x}}-T_{xz}={\rho }_3\ddot{u}+{\rho }_6{\displaystyle \frac{\partial \ddot{w}}{\partial x}}+{\rho }_8\ddot{\vartheta }+{\rho }_9\ddot{\eta } \text{(20a-d)}\hfill \end{array}$$

where

$$\left({\rho }_0,{\rho }_4,{\rho }_5,{\rho }_6,{\rho }_7,{\rho }_8,{\rho }_9\right)={\displaystyle \sum }_{k=1}^{N_l}\underset{z_k}{\overset{z_{k+1}}{{\displaystyle \int }}}{\rho }^{\left(k\right)}\left(1,f^2,fg,f\phi ,g^2,g\phi ,{\phi }^2\right)dz$$

(21)

3. Approximate Analytical Solutions of Post-Buckling Response

Because the thermal expansion of structure is regarded as a quasi-static process, the equilibrium equations governing the static response of the beam can be obtained from Eqs. (21) by eliminating the time-dependent terms, those are

$$\begin{array}{l}\delta u:{\displaystyle \frac{\partial N_x}{\partial x}}=0\\ \delta w:{\displaystyle \frac{\partial Q_{xz}}{\partial x}}-{\displaystyle \frac{{\partial }^2{M}_x}{\partial x^2}}+N_x{\displaystyle \frac{{\partial }^2\left(w+w^{\ast }\right)}{\partial x^2}}=0\\ \delta \vartheta :{\displaystyle \frac{\partial P_x}{\partial x}}-P_{xz}=0\\ \delta \eta :{\displaystyle \frac{\partial T_x}{\partial x}}-T_{xz}=0 \text{(22a-d)}\end{array}$$

Substituting Eq. (6), Eq. (7) and Eq. (14) into Eq. (22) leads to the governing equations expressed by generalized displacements

$$\begin{array}{l}{\displaystyle \frac{\partial }{\partial x}}\left\{A_{11}\left[{\displaystyle \frac{\partial u}{\partial x}}+{\displaystyle \frac{1}{2}}{\left({\displaystyle \frac{\partial w}{\partial x}}\right)}^2+{\displaystyle \frac{\partial w}{\partial x}}{\displaystyle \frac{\partial w^{\ast }}{\partial x}}\right]+B_{11}{\displaystyle \frac{{\partial }^{2}w}{\partial x^2}}+E_{11}{\displaystyle \frac{\partial \vartheta }{\partial x}}+F_{11}{\displaystyle \frac{\partial \eta }{\partial x}}-A_{11}^T\Delta T\right\}=0\\ {\displaystyle \frac{\partial }{\partial x}}\left(D_{55}{\displaystyle \frac{\partial w}{\partial x}}+G_{55}\vartheta +H_{55}\eta \right)+N_x{\displaystyle \frac{{\partial }^2\left(w+w^{\ast }\right)}{\partial x^2}}-{\displaystyle \frac{{\partial }^2}{\partial x^2}}\left(D_{11}{\displaystyle \frac{{\partial }^{2}w}{\partial x^2}}+G_{11}{\displaystyle \frac{\partial \vartheta }{\partial x}}+H_{11}{\displaystyle \frac{\partial \eta }{\partial x}}\right)\\ -B_{11}{\displaystyle \frac{{\partial }^2}{\partial x^2}}\left[{\displaystyle \frac{\partial u}{\partial x}}+{\displaystyle \frac{1}{2}}{\left({\displaystyle \frac{\partial w}{\partial x}}\right)}^2+{\displaystyle \frac{\partial w}{\partial x}}{\displaystyle \frac{\partial w^{\ast }}{\partial x}}\right]=0\\ {\displaystyle \frac{\partial }{\partial x}}\left\{G_{11}{\displaystyle \frac{{\partial }^{2}w}{\partial x^2}}+I_{11}{\displaystyle \frac{\partial \vartheta }{\partial x}}+J_{11}{\displaystyle \frac{\partial \eta }{\partial x}}+E_{11}\left[{\displaystyle \frac{\partial u}{\partial x}}+{\displaystyle \frac{1}{2}}{\left({\displaystyle \frac{\partial w}{\partial x}}\right)}^2+{\displaystyle \frac{\partial w}{\partial x}}{\displaystyle \frac{\partial w^{\ast }}{\partial x}}\right]\right\}\\ -\left(G_{55}{\displaystyle \frac{\partial w}{\partial x}}+I_{55}\vartheta +J_{55}\eta \right)=0\\ {\displaystyle \frac{\partial }{\partial x}}\left\{H_{11}{\displaystyle \frac{{\partial }^{2}w}{\partial x^2}}+J_{11}{\displaystyle \frac{\partial \vartheta }{\partial x}}+R_{11}{\displaystyle \frac{\partial \eta }{\partial x}}+F_{11}\left[{\displaystyle \frac{\partial u}{\partial x}}+{\displaystyle \frac{1}{2}}{\left({\displaystyle \frac{\partial w}{\partial x}}\right)}^2+{\displaystyle \frac{\partial w}{\partial x}}{\displaystyle \frac{\partial w^{\ast }}{\partial x}}\right]\right\}\\ -\left(H_{55}{\displaystyle \frac{\partial w}{\partial x}}+J_{55}\vartheta +R_{55}\eta \right)=0 \text{(23a-d)}\end{array}$$

One notes that Eq. (23a) may be solved for the axial displacement u, and hence it can be eliminated from the other three equations. By integrating Eq. (23) with respect to the spatial coordinate x yields, we can obtain

$$A_{11}\left[{\displaystyle \frac{\partial u}{\partial x}}+{\displaystyle \frac{1}{2}}{\left({\displaystyle \frac{\partial w}{\partial x}}\right)}^2+{\displaystyle \frac{\partial w}{\partial x}}{\displaystyle \frac{\partial w^{\ast }}{\partial x}}\right]+B_{11}{\displaystyle \frac{{\partial }^{2}w}{\partial x^2}}+E_{11}{\displaystyle \frac{\partial \vartheta }{\partial x}}+F_{11}{\displaystyle \frac{\partial \eta }{\partial x}}=C_1=N_x+A_{11}^T\Delta T$$

(24)

In Eq. (24), $C_1$ is a constant of integration. Integrate Eq. (24) once more, and one has

$$u\left|\begin{array}{c}L\\ 0\end{array}\right.=-{\displaystyle {\int }_0^L\left[\frac{1}{2}{\left(\frac{\partial w}{\partial x}\right)}^2+\frac{\partial w}{\partial x}\frac{\partial w^{\ast }}{\partial x}\right]}dx-{\displaystyle \frac{B_{11}}{A_{11}}}{\displaystyle \frac{\partial w}{\partial x}}\left|\begin{array}{c}L\\ 0\end{array}\right.-{\displaystyle \frac{E_{11}}{A_{11}}}\vartheta \left|\begin{array}{c}L\\ 0\end{array}\right.-{\displaystyle \frac{F_{11}}{A_{11}}}\eta \left|\begin{array}{c}L\\ 0\end{array}\right.+{\displaystyle \frac{C_1}{A_{11}}}x\left|\begin{array}{c}L\\ 0\end{array}\right.$$

(25a)

and

$$\left(N_x+A_{11}^T\Delta T\right)L={\displaystyle {\int }_0^L\left\{A_{11}\left[\frac{\partial u}{\partial x}+\frac{1}{2}{\left(\frac{\partial w}{\partial x}\right)}^2+\frac{\partial w}{\partial x}\frac{\partial w^{\ast }}{\partial x}\right]+B_{11}\frac{{\partial }^{2}w}{\partial x^2}+E_{11}\frac{\partial \vartheta }{\partial x}+F_{11}\frac{\partial \eta }{\partial x}\right\}}dx$$

(25b)

For a FRC beam with immovable ends, the boundary conditions for central axial displacement can be expressed as $u\left(0\right)=u\left(L\right)=0$, and thus Eqs. (25) can be further written

$$C_1={\displaystyle \frac{1}{L}}\left\{{\displaystyle {\int }_0^L{A}_{11}\left[\frac{1}{2}{\left(\frac{\partial w}{\partial x}\right)}^2+\frac{\partial w}{\partial x}\frac{\partial w^{\ast }}{\partial x}\right]}dx+B_{11}{\displaystyle \frac{\partial w}{\partial x}}\left|\begin{array}{c}L\\ 0\end{array}\right.+E_{11}\vartheta \left|\begin{array}{c}L\\ 0\end{array}\right.+F_{11}\eta \left|\begin{array}{c}L\\ 0\end{array}\right.\right\}$$

(26a)

$$N_x={\displaystyle \frac{1}{L}}{\displaystyle {\int }_0^L\left\{\begin{array}{l}{A}_{11}\left[\frac{1}{2}{\left(\frac{\partial w}{\partial x}\right)}^2+\frac{\partial w}{\partial x}\frac{\partial w^{\ast }}{\partial x}\right]\\ +B_{11}\frac{{\partial }^{2}w}{\partial x^2}+E_{11}\frac{\partial \vartheta }{\partial x}+F_{11}\frac{\partial \eta }{\partial x}\end{array}\right\}}dx-A_{11}^T\Delta T$$

(26b)

By substituting Eq. (26a) back into Eq. (24), we can obtain

$$\begin{array}{l}{A}_{11}\left[{\displaystyle \frac{\partial u}{\partial x}}+{\displaystyle \frac{1}{2}}{\left({\displaystyle \frac{\partial w}{\partial x}}\right)}^2+{\displaystyle \frac{\partial w}{\partial x}}{\displaystyle \frac{\partial w^{\ast }}{\partial x}}\right]\\ ={\displaystyle \frac{1}{L}}\left\{\begin{array}{l}{\displaystyle {\int }_0^L{A}_{11}\left[\frac{1}{2}{\left(\frac{\partial w}{\partial x}\right)}^2+\frac{\partial w}{\partial x}\frac{\partial w^{\ast }}{\partial x}\right]}dx\\ +B_{11}{\displaystyle \frac{\partial w}{\partial x}}\left|\begin{array}{c}L\\ 0\end{array}\right.+E_{11}\vartheta \left|\begin{array}{c}L\\ 0\end{array}\right.+F_{11}\eta \left|\begin{array}{c}L\\ 0\end{array}\right.\end{array}\right\}-\left(B_{11}{\displaystyle \frac{{\partial }^{2}w}{\partial x^2}}+E_{11}{\displaystyle \frac{\partial \vartheta }{\partial x}}+F_{11}{\displaystyle \frac{\partial \eta }{\partial x}}\right)\end{array}$$

(27)

For a FRC beam with immovable fixed–fixed ends, the boundary conditions can be expressed as

$$u\left(0\right)=u\left(L\right)={\displaystyle \frac{\partial w}{\partial x}}\left(0\right)={\displaystyle \frac{\partial w}{\partial x}}\left(L\right)=\eta \left(0\right)=\eta \left(L\right)=\vartheta \left(0\right)=\vartheta \left(L\right)=0$$

(28)

and thus Eqs. (26) can be further simplified as

$$C_1={\displaystyle \frac{A_{11}}{L}}{\displaystyle {\int }_0^L\left[\frac{1}{2}{\left(\frac{\partial w}{\partial x}\right)}^2+\frac{\partial w}{\partial x}\frac{\partial w^{\ast }}{\partial x}\right]}dx$$

(29a)

$$N_x={\displaystyle \frac{A_{11}}{L}}{\displaystyle {\int }_0^L\left[\frac{1}{2}{\left(\frac{\partial w}{\partial x}\right)}^2+\frac{\partial w}{\partial x}\frac{\partial w^{\ast }}{\partial x}\right]}dx-A_{11}^T\Delta T$$

(29b)

Substitute Eq. (29a) back into Eq. (24), and we can obtain

$$A_{11}\left[\begin{array}{l}{\displaystyle \frac{\partial u}{\partial x}}+{\displaystyle \frac{1}{2}}{\left({\displaystyle \frac{\partial w}{\partial x}}\right)}^2\\ +{\displaystyle \frac{\partial w}{\partial x}}{\displaystyle \frac{\partial w^{\ast }}{\partial x}}\end{array}\right]={\displaystyle \frac{A_{11}}{L}}{\displaystyle {\int }_0^L\left[\begin{array}{l}\frac{1}{2}{\left(\frac{\partial w}{\partial x}\right)}^2\\ +\frac{\partial w}{\partial x}\frac{\partial w^{\ast }}{\partial x}\end{array}\right]}dx-\left(\begin{array}{l}{B}_{11}{\displaystyle \frac{{\partial }^{2}w}{\partial x^2}}+E_{11}{\displaystyle \frac{\partial \vartheta }{\partial x}}\\ +F_{11}{\displaystyle \frac{\partial \eta }{\partial x}}\end{array}\right)$$

(30)

From the above analysis, it can be seen that the axial force $N_x$ of a beam with two immovable ends can be expressed as

$$N_x=\left\{\begin{array}{c}{\displaystyle \frac{A_{11}}{L}}{\displaystyle {\int }_0^L\left[\frac{1}{2}{\left(\frac{\partial w}{\partial x}\right)}^2+\frac{\partial w}{\partial x}\frac{\partial w^{\ast }}{\partial x}\right]}dx-A_{11}^T\Delta T\left(\mathrm{for} \mathrm{clamped}-\mathrm{clamped} \mathrm{ends}\right)\\ {\displaystyle \frac{1}{L}}{\displaystyle {\int }_0^L\left\{\begin{array}{l}{A}_{11}\left[\frac{1}{2}{\left(\frac{\partial w}{\partial x}\right)}^2+\frac{\partial w}{\partial x}\frac{\partial w^{\ast }}{\partial x}\right]\\ +B_{11}\frac{{\partial }^{2}w}{\partial x^2}+E_{11}\frac{\partial \vartheta }{\partial x}+F_{11}\frac{\partial \eta }{\partial x}\end{array}\right\}}dx-A_{11}^T\Delta T\left(\begin{array}{c}\mathrm{for} \mathrm{hinged}-\mathrm{hinged} \mathrm{ends}\\ \mathrm{or} \mathrm{hinged}-\mathrm{clamped} \mathrm{ends}\end{array}\right)\end{array}\right.$$

(31)

By the way, Eqs. (29) and Eqs. (30) are still true when the laminated beam is symmetrically layered even if the boundary conditions at both ends are not clamped-clamped.

3.1. Clamped-Clamped Supported

By substituting Eqs. (29b) and Eq. (30) into Eqs (23b-d), the simplified governing equations for clamped-clamped ends can be obtained

$$\begin{array}{l}\left(D_{11}-{\displaystyle \frac{{B_{11}}^2}{A_{11}}}\right){\displaystyle \frac{{\partial }^{4}w}{\partial x^4}}+\left(G_{11}-{\displaystyle \frac{B_{11}{E}_{11}}{A_{11}}}\right){\displaystyle \frac{{\partial }^3\vartheta }{\partial x^3}}+\left(H_{11}-{\displaystyle \frac{B_{11}{F}_{11}}{A_{11}}}\right){\displaystyle \frac{{\partial }^3\eta }{\partial x^3}}\\ -\left\{{\displaystyle \frac{A_{11}}{2L}}\underset{0}{\overset{L}{{\displaystyle \int }}}\left[{\left({\displaystyle \frac{\partial w}{\partial x}}\right)}^2+2{\displaystyle \frac{\partial w}{\partial x}}{\displaystyle \frac{\partial w{}_{*}}{\partial x}}\right]dx-A_{11}^T\Delta T\right\}{\displaystyle \frac{{\partial }^2\left(w+w^{*}\right)}{\partial x^2}}-D_{55}{\displaystyle \frac{{\partial }^{2}w}{\partial x^2}}-G_{55}{\displaystyle \frac{\partial \vartheta }{\partial x}}-H_{55}{\displaystyle \frac{\partial \eta }{\partial x}}=0\\ \left(G_{11}-{\displaystyle \frac{B_{11}{E}_{11}}{A_{11}}}\right){\displaystyle \frac{{\partial }^{3}w}{\partial x^3}}+\left(I_{11}-{\displaystyle \frac{{E_{11}}^2}{A_{11}}}\right){\displaystyle \frac{{\partial }^2\vartheta }{\partial x^2}}+\left(J_{11}-{\displaystyle \frac{E_{11}{F}_{11}}{A_{11}}}\right){\displaystyle \frac{{\partial }^2\eta }{\partial x^2}}-\left(G_{55}{\displaystyle \frac{\partial w}{\partial x}}+I_{55}\vartheta +J_{55}\eta \right)=0\\ \left(H_{11}-{\displaystyle \frac{B_{11}{F}_{11}}{A_{11}}}\right){\displaystyle \frac{{\partial }^{3}w}{\partial x^3}}+\left(J_{11}-{\displaystyle \frac{E_{11}{F}_{11}}{A_{11}}}\right){\displaystyle \frac{{\partial }^2\vartheta }{\partial x^2}}+\left(R_{11}-{\displaystyle \frac{{F_{11}}^2}{A_{11}}}\right){\displaystyle \frac{{\partial }^2\eta }{\partial x^2}}-\left(H_{55}{\displaystyle \frac{\partial w}{\partial x}}+J_{55}\vartheta +R_{55}\eta \right)=0 \text{(32a-c)}\end{array}$$

To facilitate the mathematical treatment of the governing equations, the following dimensionless parameters are defined

$$\begin{array}{l}X=\pi {\displaystyle \frac{x}{L}},W={\displaystyle \frac{w}{L}},\Theta ={\displaystyle \frac{\vartheta }{\pi }},H={\displaystyle \frac{\eta }{\pi }},{\lambda }_T=\Delta T,{\gamma }_T={\displaystyle \frac{A_{11}^T{L}^2}{{\pi }^2{D}_0}},{\gamma }_0={\displaystyle \frac{A_{11}{L}^2}{{\pi }^2{D}_0}},W^{*}={\displaystyle \frac{w^{*}}{L}}\\ \left({\gamma }_1,{\gamma }_2,{\gamma }_3,{\gamma }_4,{\gamma }_5,{\gamma }_6,{\gamma }_7,{\gamma }_8,{\gamma }_9\right)={\displaystyle \frac{1}{D_0}}\left(s_1,s_2,s_3,s_4,s_5,s_6,s_7,s_8,s_9\right)\\ \left({\gamma }_{10},{\gamma }_{11},{\gamma }_{12},{\gamma }_{13},{\gamma }_{14},{\gamma }_{15}\right)={\displaystyle \frac{L^2}{{\pi }^2{D}_0}}\left(D_{55},G_{55},H_{55},I_{55},J_{55},R_{55}\right)\end{array}$$

(33)

in which

$$\begin{array}{l}{s}_1=D_{11}-{\displaystyle \frac{{B_{11}}^2}{A_{11}}},s_2=G_{11}-{\displaystyle \frac{B_{11}{E}_{11}}{A_{11}}},s_3=H_{11}-{\displaystyle \frac{B_{11}{F}_{11}}{A_{11}}}{s}_4=G_{11}-{\displaystyle \frac{B_{11}{E}_{11}}{A_{11}}},s_5=I_{11}-{\displaystyle \frac{{E_{11}}^2}{A_{11}}},\\ s_6=J_{11}-{\displaystyle \frac{E_{11}{F}_{11}}{A_{11}}},s_7=H_{11}-{\displaystyle \frac{B_{11}{F}_{11}}{A_{11}}},s_8=J_{11}-{\displaystyle \frac{E_{11}{F}_{11}}{A_{11}}},s_9=R_{11}-{\displaystyle \frac{{F_{11}}^2}{A_{11}}},D_0={\displaystyle \sum }_{k=1}^{N_l}\underset{z_k}{\overset{z_{k+1}}{{\displaystyle \int }}} {\tilde{Q}}_{11}^{\left(k\right)}{z}^{2}dz\end{array}$$

(34)

Substitute Eqs. (33) into Eqs. (32) and neglect the higher-order small terms, the following dimensionless equations are obtained

$$\begin{array}{l}{\gamma }_1{\displaystyle \frac{{\partial }^{4}W}{\partial X^4}}+{\gamma }_2{\displaystyle \frac{{\partial }^3\Theta }{\partial X^3}}+{\gamma }_3{\displaystyle \frac{{\partial }^{3}H}{\partial X^3}}-\left[{\displaystyle \frac{\pi {\gamma }_0}{2}}\underset{0}{\overset{\pi }{{\displaystyle \int }}}{\left({\displaystyle \frac{\partial W}{\partial X}}\right)}^{2}dX-{\gamma }_T{\lambda }_T\right]{\displaystyle \frac{{\partial }^2\left(W+W^{*}\right)}{\partial X^2}}\\ -{\displaystyle \frac{{\partial }^{2}W}{\partial X^2}}\left(\pi {\gamma }_0\underset{0}{\overset{\pi }{{\displaystyle \int }}}{\displaystyle \frac{\partial W}{\partial X}}{\displaystyle \frac{\partial W^{*}}{\partial X}}dX\right)-{\gamma }_{10}{\displaystyle \frac{{\partial }^{2}W}{\partial X^2}}-{\gamma }_{11}{\displaystyle \frac{\partial \Theta }{\partial X}}-{\gamma }_{12}{\displaystyle \frac{\partial H}{\partial X}}=0\\ {\gamma }_4{\displaystyle \frac{{\partial }^{3}W}{\partial X^3}}+{\gamma }_5{\displaystyle \frac{{\partial }^2\Theta }{\partial X^2}}+{\gamma }_6{\displaystyle \frac{{\partial }^{2}H}{\partial X^2}}-\left({\gamma }_{11}{\displaystyle \frac{\partial W}{\partial X}}+{\gamma }_{13}\Theta +{\gamma }_{14}H\right)=0\\ {\gamma }_7{\displaystyle \frac{{\partial }^{3}W}{\partial X^3}}+{\gamma }_8{\displaystyle \frac{{\partial }^2\Theta }{\partial X^2}}+{\gamma }_9{\displaystyle \frac{{\partial }^{2}H}{\partial X^2}}-\left({\gamma }_{12}{\displaystyle \frac{\partial W}{\partial X}}+{\gamma }_{14}\Theta +{\gamma }_{15}H\right)=0 \text{(35a-c)}\end{array}$$

The dimensionless boundary conditions can be expressed as

$$W\left(0\right)=W\left(\pi \right)={\displaystyle \frac{\partial W}{\partial x}}\left(0\right)={\displaystyle \frac{\partial W}{\partial x}}\left(\pi \right)=H\left(0\right)=H\left(\pi \right)=\Theta \left(\pi \right)=\Theta \left(\pi \right)=0$$

(36)

Although the initial defects of the beam are random, studies have shown that the influence of the waveform defect is the most prominent, so it is assumed that the initial defect of the beam has the same form as the classical solution of the small deflection [10], that is

$$W^{*}=\epsilon A_{10}^{*}\left[1-cos\left(2mX\right)\right]=\epsilon \mu A_{10}^{\left(1\right)}\left[1-cos\left(2mX\right)\right]$$

(37)

where $\mu ={\displaystyle \frac{A_{10}^{*}}{A_{10}^{\left(1\right)}}}$ is the imperfection parameter. It is worth noting that in the post-buckling stage, the defect parameter is a variable rather than a constant.

A two-step perturbation method is employed to solve Eq. (35) and Eq. (36) and determine the thermal post-buckling equilibrium paths [40]. In the present case, thermal load and the generalized displacements are represented in the following series form,

$$\begin{array}{l}{\lambda }_T\left(X,\epsilon \right)={\displaystyle \sum }_{k=0}{\epsilon }^k{\lambda }_T^{\left(k\right)}\left(X\right)，W\left(X,\epsilon \right)={\displaystyle \sum }_{k=1}{\epsilon }^k{w}_k\left(X\right)\\ \Theta \left(X,\epsilon \right)={\displaystyle \sum }_{k=1}{\epsilon }^k{\vartheta }_k\left(X\right)，H\left(X,\epsilon \right)={\displaystyle \sum }_{k=1}{\epsilon }^k{\eta }_k\left(X\right)\end{array}$$

(38)

where $\epsilon$ is a small perturbation parameter with no physical meaning, and k is the number of terms in series.

By substituting Eq. (38) into Eq. (35), a set of the perturbation differential equations are derived. The first-order perturbation equations are expressed as

$O\left({\epsilon }^1\right)$:

$$\begin{array}{l}{\gamma }_1{\displaystyle \frac{{\partial }^4{w}_1}{\partial X^4}}+{\gamma }_2{\displaystyle \frac{{\partial }^3{\vartheta }_1}{\partial X^3}}+{\gamma }_3{\displaystyle \frac{{\partial }^3{\eta }_1}{\partial X^3}}+{\gamma }_T{\lambda }_T^{\left(0\right)}{\displaystyle \frac{{\partial }^2}{\partial X^2}}\left(w_1+{\displaystyle \frac{W^{*}}{\epsilon }}\right)-{\gamma }_{11}{\displaystyle \frac{\partial {\vartheta }_1}{\partial X}}-{\gamma }_{12}{\displaystyle \frac{\partial {\eta }_1}{\partial X}}-{\gamma }_{10}{\displaystyle \frac{{\partial }^2{w}_1}{\partial X^2}}=0\\ {\gamma }_4{\displaystyle \frac{{\partial }^3{w}_1}{\partial X^3}}+{\gamma }_5{\displaystyle \frac{{\partial }^2{\vartheta }_1}{\partial X^2}}+{\gamma }_6{\displaystyle \frac{{\partial }^2{\eta }_1}{\partial X^2}}-\left({\gamma }_{11}{\displaystyle \frac{\partial w_1}{\partial X}}+{\gamma }_{13}{\vartheta }_1+{\gamma }_{14}{\eta }_1\right)=0\\ {\gamma }_7{\displaystyle \frac{{\partial }^3{w}_1}{\partial X^3}}+{\gamma }_8{\displaystyle \frac{{\partial }^2{\vartheta }_1}{\partial X^2}}+{\gamma }_9{\displaystyle \frac{{\partial }^2{\eta }_1}{\partial X^2}}-\left({\gamma }_{12}{\displaystyle \frac{\partial w_1}{\partial X}}+{\gamma }_{14}{\vartheta }_1+{\gamma }_{15}{\eta }_1\right)=0 \text{(39a-c)}\end{array}$$

To solve the above partial differential equations, a set of displacement trial functions that satisfy the boundary conditions is introduced [41]

$$\begin{array}{l}{w}_1\left(X\right)=A_{10}^{\left(1\right)}\left[1-cos\left(2mX\right)\right]\\ {\vartheta }_1\left(X\right)=B_{10}^{\left(1\right)}sin\left(2mX\right)\\ {\eta }_1\left(X\right)=C_{10}^{\left(1\right)}sin\left(2mX\right) \text{(40a-c)}\end{array}$$

Substitute Eqs. (37) and Eq. (40) into Eqs. (39), and a set of linear equations can be obtained.

$$\begin{array}{l}\begin{array}{c}\left[{\gamma }_{10}-{\gamma }_T{\lambda }_T^{\left(0\right)}\left(1+\mu \right)+{\left(2m\right)}^2{\gamma }_1\right]{\left(2m\right)}^2{A}_{10}^{\left(1\right)}\\ +\left[{\gamma }_{11}\left(2m\right)+{\gamma }_2{\left(2m\right)}^3\right]B_{10}^{\left(1\right)}+\left[{\gamma }_{12}\left(2m\right)+{\gamma }_3{\left(2m\right)}^3\right]C_{10}^{\left(1\right)}=0\end{array}\\ \left[{\left(2m\right)}^3{\gamma }_4+2m{\gamma }_{11}\right]A_{10}^{\left(1\right)}+\left[{\left(2m\right)}^2{\gamma }_5+{\gamma }_{13}\right]B_{10}^{\left(1\right)}+\left[{\left(2m\right)}^2{\gamma }_6+{\gamma }_{14}\right]C_{10}^{\left(1\right)}=0\\ \left[{\left(2m\right)}^3{\gamma }_7+2m{\gamma }_{12}\right]A_{10}^{\left(1\right)}+\left[{\left(2m\right)}^2{\gamma }_8+{\gamma }_{14}\right]B_{10}^{\left(1\right)}+\left[{\left(2m\right)}^2{\gamma }_9+{\gamma }_{15}\right]C_{10}^{\left(1\right)}=0 \text{(41a-c)}\end{array}$$

By solving Eqs. (41), the zero-order buckling temperature solutions can be obtained

$${\lambda }_T^{\left(0\right)}={\displaystyle \frac{\left[{\gamma }_{11}+{\gamma }_2{\left(2m\right)}^2\right]\frac{B_{10}^{\left(1\right)}}{A_{10}^{\left(1\right)}}+\left[{\gamma }_{12}+{\gamma }_3{\left(2m\right)}^2\right]\frac{C_{10}^{\left(1\right)}}{A_{10}^{\left(1\right)}}}{\left(2m\right)\left(1+\mu \right){\gamma }_T}}+{\displaystyle \frac{{\left(2m\right)}^2{\gamma }_1+{\gamma }_{10}}{\left(1+\mu \right){\gamma }_T}}$$

(42)

in which

$${\displaystyle \frac{B_{10}^{\left(1\right)}}{A_{10}^{\left(1\right)}}}={\displaystyle \frac{\left|\begin{array}{cc}-\left[{\left(2m\right)}^3{\gamma }_4+2m{\gamma }_{11}\right]& \left[{\left(2m\right)}^2{\gamma }_6+{\gamma }_{14}\right]\\ -\left[{\left(2m\right)}^3{\gamma }_7+2m{\gamma }_{12}\right]& \left[{\left(2m\right)}^2{\gamma }_9+{\gamma }_{15}\right]\end{array}\right|}{\left|\begin{array}{cc}\left[{\left(2m\right)}^2{\gamma }_5+{\gamma }_{13}\right]& \left[{\left(2m\right)}^2{\gamma }_6+{\gamma }_{14}\right]\\ \left[{\left(2m\right)}^2{\gamma }_8+{\gamma }_{14}\right]& \left[{\left(2m\right)}^2{\gamma }_9+{\gamma }_{15}\right]\end{array}\right|}}$$

(43a)

$${\displaystyle \frac{C_{10}^{\left(1\right)}}{A_{10}^{\left(1\right)}}}={\displaystyle \frac{\left|\begin{array}{cc}\left[{\left(2m\right)}^2{\gamma }_5+{\gamma }_{13}\right]& -\left[{\left(2m\right)}^3{\gamma }_4+2m{\gamma }_{11}\right]\\ \left[{\left(2m\right)}^2{\gamma }_8+{\gamma }_{14}\right]& -\left[{\left(2m\right)}^3{\gamma }_7+2m{\gamma }_{12}\right]\end{array}\right|}{\left|\begin{array}{cc}\left[{\left(2m\right)}^2{\gamma }_5+{\gamma }_{13}\right]& \left[{\left(2m\right)}^2{\gamma }_6+{\gamma }_{14}\right]\\ \left[{\left(2m\right)}^2{\gamma }_8+{\gamma }_{14}\right]& \left[{\left(2m\right)}^2{\gamma }_9+{\gamma }_{15}\right]\end{array}\right|}}$$

(43b)

The second order equations can be written as

$O\left({\epsilon }^2\right)$:

$$\begin{array}{l}{\gamma }_1{\displaystyle \frac{{\partial }^4{w}_2}{\partial X^4}}+{\gamma }_2{\displaystyle \frac{{\partial }^3{\vartheta }_2}{\partial X^3}}+{\gamma }_3{\displaystyle \frac{{\partial }^3{\eta }_2}{\partial X^3}}-\left[{\gamma }_{10}-{\gamma }_T{\lambda }_T^{\left(0\right)}\right]{\displaystyle \frac{{\partial }^2{w}_2}{\partial X^2}}-{\gamma }_{11}{\displaystyle \frac{\partial {\vartheta }_2}{\partial X}}-{\gamma }_{12}{\displaystyle \frac{\partial {\eta }_2}{\partial X}}\\ =-{\gamma }_T{\lambda }_T^{\left(1\right)}{\displaystyle \frac{{\partial }^2\left(w_1+W^{*}/\epsilon \right)}{\partial X^2}}\\ {\gamma }_4{\displaystyle \frac{{\partial }^3{w}_2}{\partial X^3}}+{\gamma }_5{\displaystyle \frac{{\partial }^2{\vartheta }_2}{\partial X^2}}+{\gamma }_6{\displaystyle \frac{{\partial }^2{\eta }_2}{\partial X^2}}-\left({\gamma }_{11}{\displaystyle \frac{\partial w_2}{\partial X}}+{\gamma }_{13}{\vartheta }_2+{\gamma }_{14}{\eta }_2\right)=0\\ {\gamma }_7{\displaystyle \frac{{\partial }^3{w}_2}{\partial X^3}}+{\gamma }_8{\displaystyle \frac{{\partial }^2{\vartheta }_2}{\partial X^2}}+{\gamma }_9{\displaystyle \frac{{\partial }^2{\eta }_2}{\partial X^2}}-\left({\gamma }_{12}{\displaystyle \frac{\partial w_2}{\partial X}}+{\gamma }_{14}{\vartheta }_2+{\gamma }_{15}{\eta }_2\right)=0 \text{(44a-c)}\end{array}$$

According to the form of equations, a set of trial functions of second-order generalized displacements satisfying boundary conditions is introduced again.

$$\begin{array}{l}{w}_2\left(X\right)=A_{10}^{\left(2\right)}\left[1-cos\left(2mX\right)\right]\\ {\vartheta }_2\left(X\right)=B_{10}^{\left(2\right)}sin\left(2mX\right)\\ {\eta }_2\left(X\right)=C_{10}^{\left(2\right)}sin\left(2mX\right) \text{(45a-c)}\end{array}$$

Substitute Eqs. (45) and Eq. (37) into Eqs. (44), and a set of linear second-order equations can be obtained.

$$\begin{array}{l}\left[{\gamma }_{10}-{\gamma }_T{\lambda }_T^{\left(0\right)}+{\left(2m\right)}^2{\gamma }_1\right]{\left(2m\right)}^2{A}_{10}^{\left(2\right)}+\left[{\gamma }_{11}\left(2m\right)+{\gamma }_2{\left(2m\right)}^3\right]B_{10}^{\left(2\right)}\\ +\left[{\gamma }_{12}\left(2m\right)+{\gamma }_3{\left(2m\right)}^3\right]C_{10}^{\left(2\right)}=-{\gamma }_T{\lambda }_T^{\left(1\right)}{\left(2m\right)}^2{A}_{10}^{\left(1\right)}\\ \left[{\left(2m\right)}^3{\gamma }_4+2m{\gamma }_{11}\right]A_{10}^{\left(2\right)}+\left[{\left(2m\right)}^2{\gamma }_5+{\gamma }_{13}\right]B_{10}^{\left(2\right)}+\left[{\left(2m\right)}^2{\gamma }_6+{\gamma }_{14}\right]C_{10}^{\left(2\right)}=0\\ \left[{\left(2m\right)}^3{\gamma }_7+2m{\gamma }_{12}\right]A_{10}^{\left(2\right)}+\left[{\left(2m\right)}^2{\gamma }_8+{\gamma }_{14}\right]B_{10}^{\left(2\right)}+\left[{\left(2m\right)}^2{\gamma }_9+{\gamma }_{15}\right]C_{10}^{\left(2\right)}=0 \text{(46a-c)}\end{array}$$

It can be proved that the unique solutions of Eqs. (46) are as follows

$${\lambda }_T^{\left(1\right)}=0$$

(47a)

$$w_2={\vartheta }_2={\eta }_2=0$$

(47b)

The third-order equation can be written as

$O\left({\epsilon }^3\right)$:

$$\begin{array}{l}{\gamma }_1{\displaystyle \frac{{\partial }^4{w}_3}{\partial X^4}}+{\gamma }_2{\displaystyle \frac{{\partial }^3{\vartheta }_3}{\partial X^3}}+{\gamma }_3{\displaystyle \frac{{\partial }^3{\eta }_3}{\partial X^3}}-\left[{\gamma }_{10}-{\gamma }_T{\lambda }_T^{\left(0\right)}\right]{\displaystyle \frac{{\partial }^2{w}_3}{\partial X^2}}-{\gamma }_{11}{\displaystyle \frac{\partial {\vartheta }_3}{\partial X}}-{\gamma }_{12}{\displaystyle \frac{\partial {\eta }_3}{\partial X}}=\\ \left[{\displaystyle \frac{\pi {\gamma }_0}{2}}\underset{0}{\overset{\pi }{{\displaystyle \int }}}{\left({\displaystyle \frac{\partial w_1}{\partial X}}\right)}^{2}dX\right]{\displaystyle \frac{{\partial }^2\left(w_1+\frac{W^{*}}{\epsilon }\right)}{\partial X^2}}+{\displaystyle \frac{{\partial }^2{w}_1}{\partial X^2}}\left(\pi {\gamma }_0\underset{0}{\overset{\pi }{{\displaystyle \int }}}{\displaystyle \frac{\partial w_1}{\partial X}}{\displaystyle \frac{\partial \frac{W^{*}}{\epsilon }}{\partial X}}dX\right)-{\gamma }_T{\lambda }_T^{\left(2\right)}{\displaystyle \frac{{\partial }^2\left(w_1+\frac{W^{*}}{\epsilon }\right)}{\partial X^2}}\\ {\gamma }_4{\displaystyle \frac{{\partial }^3{w}_3}{\partial X^3}}+{\gamma }_5{\displaystyle \frac{{\partial }^2{\vartheta }_3}{\partial X^2}}+{\gamma }_6{\displaystyle \frac{{\partial }^2{\eta }_3}{\partial X^2}}-\left({\gamma }_{11}{\displaystyle \frac{\partial w_3}{\partial X}}+{\gamma }_{13}{\vartheta }_3+{\gamma }_{14}{\eta }_3\right)=0\\ {\gamma }_7{\displaystyle \frac{{\partial }^3{w}_3}{\partial X^3}}+{\gamma }_8{\displaystyle \frac{{\partial }^2{\vartheta }_3}{\partial X^2}}+{\gamma }_9{\displaystyle \frac{{\partial }^2{\eta }_3}{\partial X^2}}-\left({\gamma }_{12}{\displaystyle \frac{\partial w_3}{\partial X}}+{\gamma }_{14}{\vartheta }_3+{\gamma }_{15}{\eta }_3\right)=0 \text{(48a-c)}\end{array}$$

According to the form of the equations, a set of trial functions of third-order generalized displacements satisfying boundary conditions is introduced again.

$$\begin{array}{l}{w}_3\left(X\right)=A_{10}^{\left(3\right)}\left[1-cos\left(2mX\right)\right]\\ {\vartheta }_3\left(X\right)=B_{10}^{\left(3\right)}sin\left(2mX\right)\\ {\eta }_3\left(X\right)=C_{10}^{\left(3\right)}sin\left(2mX\right) \text{(49a-c)}\end{array}$$

By substituting Eqs. (37) and Eq. (49) into Eqs. (48) and solving the obtained set of linear third-order equations in a similar way, the second order buckling temperature solutions can be obtained

$${\lambda }_T^{\left(2\right)}={\displaystyle \frac{{\gamma }_0\left(1+3\mu \right){\left(2m\pi A_{10}^{\left(1\right)}\right)}^2}{4{\gamma }_T\left(1+\mu \right)}}$$

(50)

And

$$w_3\left(X\right)={\vartheta }_3\left(X\right)={\eta }_3\left(X\right)=0$$

(51)

Sufficient accuracy can be obtained by solving up to the third-order perturbation equations, according to [41]. The asymptotic solutions of the dimensionless deflection and thermal load can be written as the form of series expansion.

$$W\left(X,\epsilon \right)={\displaystyle \sum }_{k=1}{\epsilon }^k{w}_k\left(X\right)=\epsilon w_1+O\left({\epsilon }^4\right)=\epsilon A_{10}^{\left(1\right)}\left[1-cos\left(2mX\right)\right]$$

(52a)

$$\begin{array}{l}{\lambda }_T={\lambda }_T^{\left(0\right)}+{\epsilon }^2{\lambda }_T^{\left(2\right)}+O\left({\epsilon }^4\right)={\displaystyle \frac{\left[{\gamma }_{11}+{\gamma }_2{\left(2m\right)}^2\right]\frac{B_{10}^{\left(1\right)}}{A_{10}^{\left(1\right)}}+\left[{\gamma }_{12}+{\gamma }_3{\left(2m\right)}^2\right]\frac{C_{10}^{\left(1\right)}}{A_{10}^{\left(1\right)}}}{\left(2m\right)\left(1+\mu \right){\gamma }_T}}\\ +{\displaystyle \frac{{\left(2m\right)}^2{\gamma }_1+{\gamma }_{10}}{\left(1+\mu \right){\gamma }_T}}+{\displaystyle \frac{{\gamma }_0\left(1+3\mu \right){\left(2m\pi A_{10}^{\left(1\right)}\right)}^2}{4{\gamma }_T\left(1+\mu \right)}}\end{array}$$

(52b)

According to Eq. (52a), we can obtain

$$W_m=W{|}_{X={\displaystyle \frac{\pi }{2m}}}=2\epsilon A_{10}^{\left(1\right)}$$

(53)

Thus, the thermal post-buckling load-deflection relationships are obtained

$$\begin{array}{l}{\lambda }_T={\lambda }_T^{\left(0\right)}+{\epsilon }^2{\lambda }_T^{\left(2\right)}+O\left({\epsilon }^4\right)\\ ={\displaystyle \frac{\left[{\gamma }_{11}+{\gamma }_2{\left(2m\right)}^2\right]\frac{B_{10}^{\left(1\right)}}{A_{10}^{\left(1\right)}}+\left[{\gamma }_{12}+{\gamma }_3{\left(2m\right)}^2\right]\frac{C_{10}^{\left(1\right)}}{A_{10}^{\left(1\right)}}}{\left(2m\right)\left(1+\mu \right){\gamma }_T}}+{\displaystyle \frac{{\left(2m\right)}^2{\gamma }_1+{\gamma }_{10}}{\left(1+\mu \right){\gamma }_T}}+{\displaystyle \frac{{\gamma }_0\left(1+3\mu \right){\left(m\pi W_m\right)}^2}{4\left(1+\mu \right){\gamma }_T}}+\dots \end{array}$$

(54)

in which the expression of ${\displaystyle \frac{B_{10}^{\left(1\right)}}{A_{10}^{\left(1\right)}}}$ and ${\displaystyle \frac{C_{10}^{\left(1\right)}}{A_{10}^{\left(1\right)}}}$ can be referred to Eqs. (43).

It's remarkable that if there are no geometric defects in FRC laminated beams with clamped-clamped ends, the post-buckling thermal load-deflection curve of bifurcation type can be obtained, that is

$$\begin{array}{l}{\lambda }_T={\lambda }_T^{\left(0\right)}+{\epsilon }^2{\lambda }_T^{\left(2\right)}+O\left({\epsilon }^4\right)\\ ={\displaystyle \frac{\left[{\gamma }_{11}+{\gamma }_2{\left(2m\right)}^2\right]\frac{B_{10}^{\left(1\right)}}{A_{10}^{\left(1\right)}}+\left[{\gamma }_{12}+{\gamma }_3{\left(2m\right)}^2\right]\frac{C_{10}^{\left(1\right)}}{A_{10}^{\left(1\right)}}}{\left(2m\right){\gamma }_T}}+{\displaystyle \frac{{\left(2m\right)}^2{\gamma }_1+{\gamma }_{10}}{{\gamma }_T}}+{\displaystyle \frac{{\gamma }_0{\left(m\pi W_m\right)}^2}{4{\gamma }_T}}+\dots \end{array}$$

(55)

Obviously, the equilibrium path of post-buckling is stable, and by letting $W_m=0$ , the critical buckling bifurcation points can be obtained

$${\lambda }_{Tcr}={\displaystyle \frac{\left[{\gamma }_{11}+{\gamma }_2{\left(2m\right)}^2\right]\frac{B_{10}^{\left(1\right)}}{A_{10}^{\left(1\right)}}+\left[{\gamma }_{12}+{\gamma }_3{\left(2m\right)}^2\right]\frac{C_{10}^{\left(1\right)}}{A_{10}^{\left(1\right)}}}{\left(2m\right){\gamma }_T}}+{\displaystyle \frac{{\left(2m\right)}^2{\gamma }_1+{\gamma }_{10}}{{\gamma }_T}}$$

(56)

Note that in analytical formulas, m represents the order of buckling modes. By letting m=1, the real post-buckling equilibrium path of clamped-clamped FRC laminated beams can be obtained.

3.2. Hinged-Hinged Supported

By substituting Eqs. (26b) and Eq. (27) into Eqs (23b-d), and neglecting the higher-order small terms, the simplified governing equations for FRC laminated beams with hinged-hinged ends can be obtained

$$\begin{array}{l}\left(D_{11}-{\displaystyle \frac{{B_{11}}^2}{A_{11}}}\right){\displaystyle \frac{{\partial }^{4}w}{\partial x^4}}+\left(G_{11}-{\displaystyle \frac{B_{11}{E}_{11}}{A_{11}}}\right){\displaystyle \frac{{\partial }^3\vartheta }{\partial x^3}}+\left(H_{11}-{\displaystyle \frac{B_{11}{F}_{11}}{A_{11}}}\right){\displaystyle \frac{{\partial }^3\eta }{\partial x^3}}-D_{55}{\displaystyle \frac{{\partial }^{2}w}{\partial x^2}}-G_{55}{\displaystyle \frac{\partial \vartheta }{\partial x}}-H_{55}{\displaystyle \frac{\partial \eta }{\partial x}}\\ -\left\{{\displaystyle \frac{1}{L}}\underset{0}{\overset{L}{{\displaystyle \int }}}\left\{A_{11}\left[{\displaystyle \frac{1}{2}}{\left({\displaystyle \frac{\partial w}{\partial x}}\right)}^2+{\displaystyle \frac{\partial w}{\partial x}}{\displaystyle \frac{\partial w{}_{*}}{\partial x}}\right]+B_{11}{\displaystyle \frac{{\partial }^{2}w}{\partial x^2}}+E_{11}{\displaystyle \frac{\partial \vartheta }{\partial x}}+F_{11}{\displaystyle \frac{\partial \eta }{\partial x}}\right\}dx-A_{11}^T\Delta T\right\}{\displaystyle \frac{{\partial }^2\left(w+w^{*}\right)}{\partial x^2}}=0\\ \left(G_{11}-{\displaystyle \frac{B_{11}{E}_{11}}{A_{11}}}\right){\displaystyle \frac{{\partial }^{3}w}{\partial x^3}}+\left(I_{11}-{\displaystyle \frac{{E_{11}}^2}{A_{11}}}\right){\displaystyle \frac{{\partial }^2\vartheta }{\partial x^2}}+\left(J_{11}-{\displaystyle \frac{E_{11}{F}_{11}}{A_{11}}}\right){\displaystyle \frac{{\partial }^2\eta }{\partial x^2}}-\left(G_{55}{\displaystyle \frac{\partial w}{\partial x}}+I_{55}\vartheta +J_{55}\eta \right)=0\\ \left(H_{11}-{\displaystyle \frac{B_{11}{F}_{11}}{A_{11}}}\right){\displaystyle \frac{{\partial }^{3}w}{\partial x^3}}+\left(J_{11}-{\displaystyle \frac{E_{11}{F}_{11}}{A_{11}}}\right){\displaystyle \frac{{\partial }^2\vartheta }{\partial x^2}}+\left(R_{11}-{\displaystyle \frac{{F_{11}}^2}{A_{11}}}\right){\displaystyle \frac{{\partial }^2\eta }{\partial x^2}}-\left(H_{55}{\displaystyle \frac{\partial w}{\partial x}}+J_{55}\vartheta +R_{55}\eta \right)=0 \text{(57a-c)}\end{array}$$

On the basis of Eqs (33), the following dimensionless parameters are further introduced.

$${\varsigma }_1={\displaystyle \frac{B_{11}L}{{\pi }^2{D}_0}},{\varsigma }_2={\displaystyle \frac{E_{11}L}{{\pi }^2{D}_0}},{\varsigma }_3={\displaystyle \frac{F_{11}L}{{\pi }^2{D}_0}}$$

(58)

where the expression of $D_0$ can be referred to Eqs. (34).

Consequently, the dimensionless governing equations can be expressed as follows

$$\begin{array}{l}{\gamma }_1{\displaystyle \frac{{\partial }^{4}W}{\partial X^4}}+{\gamma }_2{\displaystyle \frac{{\partial }^3\Theta }{\partial X^3}}+{\gamma }_3{\displaystyle \frac{{\partial }^{3}H}{\partial X^3}}-\left\{\pi \underset{0}{\overset{\pi }{{\displaystyle \int }}}\left[\begin{array}{l}{\displaystyle \frac{{\gamma }_0}{2}}{\left({\displaystyle \frac{\partial W}{\partial X}}\right)}^2+{\varsigma }_1{\displaystyle \frac{{\partial }^{2}W}{\partial X^2}}\\ +{\varsigma }_2{\displaystyle \frac{\partial \Theta }{\partial X}}+{\varsigma }_3{\displaystyle \frac{\partial H}{\partial X}}\end{array}\right]dX-{\gamma }_T{\lambda }_T\right\}{\displaystyle \frac{{\partial }^2\left(W+W^{*}\right)}{\partial X^2}}\\ -{\displaystyle \frac{{\partial }^{2}W}{\partial X^2}}\left(\pi {\gamma }_0\underset{0}{\overset{\pi }{{\displaystyle \int }}}{\displaystyle \frac{\partial W}{\partial X}}{\displaystyle \frac{\partial W^{*}}{\partial X}}dX\right)-{\gamma }_{10}{\displaystyle \frac{{\partial }^{2}W}{\partial X^2}}-{\gamma }_{11}{\displaystyle \frac{\partial \Theta }{\partial X}}-{\gamma }_{12}{\displaystyle \frac{\partial H}{\partial X}}=0\\ {\gamma }_4{\displaystyle \frac{{\partial }^{3}W}{\partial X^3}}+{\gamma }_5{\displaystyle \frac{{\partial }^2\Theta }{\partial X^2}}+{\gamma }_6{\displaystyle \frac{{\partial }^{2}H}{\partial X^2}}-\left({\gamma }_{11}{\displaystyle \frac{\partial W}{\partial X}}+{\gamma }_{13}\Theta +{\gamma }_{14}H\right)=0\\ {\gamma }_7{\displaystyle \frac{{\partial }^{3}W}{\partial X^3}}+{\gamma }_8{\displaystyle \frac{{\partial }^2\Theta }{\partial X^2}}+{\gamma }_9{\displaystyle \frac{{\partial }^{2}H}{\partial X^2}}-\left({\gamma }_{12}{\displaystyle \frac{\partial W}{\partial X}}+{\gamma }_{14}\Theta +{\gamma }_{15}H\right)=0 \text{(59a-c)}\end{array}$$

According to [39] , there are two factors that lead to the initial deformation $W^{*}$for simply supported asymmetric FRC laminated beam, one is the initial deformation $W_T^{*}$ caused by stretching-bending coupling effect, the other $W_G^{*}$ is still caused by geometric defects. The initial deformation of the beam is still assumed to satisfy the hinged-hinged boundary and possess the same shape as the buckling mode, that is

$$W^{*}=W_G^{*}+W_T^{*}=\epsilon A_G^{*}\sin \left(mX\right)+\epsilon A_T^{*}\sin \left(X\right)$$

(60)

Define $\mu ={\displaystyle \frac{A_T^{*}+A_G^{*}}{A_{10}^{\left(1\right)}}}$ as the imperfection parameter. The subscript T and G respectively denote thermal moments and geometric imperfection. Also, notice that the defect parameter is a variable rather than a constant in the post-buckling stage.

The solutions of Eqs. (59a-c) can also be determined by a two-step perturbation technique. The generalized displacements and thermal load can be assumed the same form as Eqs. (38), those are

$$\begin{array}{l}{\lambda }_T\left(X,\epsilon \right)={\displaystyle \sum }_{k=0}{\epsilon }^k{\lambda }_T^{\left(k\right)}\left(X\right), W\left(X,\epsilon \right)={\displaystyle \sum }_{k=1}{\epsilon }^k{w}_k\left(X\right)\\ \Theta \left(X,\epsilon \right)={\displaystyle \sum }_{k=1}{\epsilon }^k{\vartheta }_k\left(X\right), H\left(X,\epsilon \right)={\displaystyle \sum }_{k=1}{\epsilon }^k{\eta }_k\left(X\right)\end{array}$$

(61)

Similar to the methods in section 3.1, by substituting Eqs. (61) into Eqs. (59a-c) and collecting the terms of the same order of $\epsilon$, a set of perturbation equations are obtained and up to the third order perturbation equations are required to solve for post-buckling problem.

Assuming the displacement trial functions step by step and substituting them with Eq. (60) into the perturbation equations of each order for solution, the following forms of asymptotic solutions can be obtained.

$$\begin{array}{l}W=\epsilon w_1+O\left({\epsilon }^4\right)=\epsilon A_{10}^{\left(1\right)}sin\left(mX\right)+O\left({\epsilon }^4\right)\\ \Theta =\epsilon {\vartheta }_1+O\left({\epsilon }^4\right)=\epsilon B_{10}^{\left(1\right)}\cos \left(mX\right)+O\left({\epsilon }^4\right)\\ H=\epsilon {\eta }_1+O\left({\epsilon }^4\right)=\epsilon C_{10}^{\left(1\right)}\cos \left(mX\right)+O\left({\epsilon }^4\right) \text{(62a-c)}\end{array}$$

And

$$\begin{array}{l}{\lambda }_T={\lambda }_T^{\left(0\right)}+\epsilon {\lambda }_T^{\left(1\right)}+{\epsilon }^2{\lambda }_T^{\left(2\right)}+O\left({\epsilon }^4\right)={\displaystyle \frac{{\gamma }_{10}+{\gamma }_1}{{\gamma }_T\left(1+\mu \right)}}+{\displaystyle \frac{\left({\gamma }_{11}+{\gamma }_2\right)B_{10}^{\left(1\right)}}{{\gamma }_T\left(1+\mu \right)A_{10}^{\left(1\right)}}}+{\displaystyle \frac{\left({\gamma }_{12}+{\gamma }_3\right)C_{10}^{\left(1\right)}}{{\gamma }_T\left(1+\mu \right)A_{10}^{\left(1\right)}}}\\ +{\displaystyle \frac{2\pi \left(A_{10}^{\left(1\right)}\epsilon \right)\left({\varsigma }_1+{\varsigma }_2\frac{B_{10}^{\left(1\right)}}{A_{10}^{\left(1\right)}}+{\varsigma }_3\frac{C_{10}^{\left(1\right)}}{A_{10}^{\left(1\right)}}\right)}{{\gamma }_T}}+{\displaystyle \frac{\left(1+3\mu \right){\gamma }_0{\pi }^2{\left(A_{10}^{\left(1\right)}\epsilon \right)}^2}{4\left(1+\mu \right){\gamma }_T}}+\dots \left(m=1\right)\end{array}$$

(63)

In Eq. (63), similar to results of the beams with clamped-clamped ends, m=1 and the real post-buckling equilibrium path of hinged-hinged FRC laminated beams can be obtained. $\left(A_{10}^{\left(1\right)}\epsilon \right)$ is taken as the perturbation parameter relating to the dimensionless maximum deflection $W_m$ and can be obtained by substituting $x={\displaystyle \frac{\pi }{2m}}$ into Eq. (62a), that is

$$W_m=\epsilon A_{10}^{\left(1\right)}$$

(64)

Substitute Eq. (64) into Eq. (63) and the post-buckling equilibrium path (load-deflection curve) is obtained as follows

$$\begin{array}{l}{\lambda }_T={\lambda }_T^{\left(0\right)}+\epsilon {\lambda }_T^{\left(1\right)}+{\epsilon }^2{\lambda }_T^{\left(2\right)}+O\left({\epsilon }^4\right)={\displaystyle \frac{{\gamma }_{10}+{\gamma }_1}{{\gamma }_T\left(1+\mu \right)}}+{\displaystyle \frac{\left({\gamma }_{11}+{\gamma }_2\right)B_{10}^{\left(1\right)}}{{\gamma }_T\left(1+\mu \right)A_{10}^{\left(1\right)}}}+{\displaystyle \frac{\left({\gamma }_{12}+{\gamma }_3\right)C_{10}^{\left(1\right)}}{{\gamma }_T\left(1+\mu \right)A_{10}^{\left(1\right)}}}\\ +{\displaystyle \frac{2\pi W_m\left({\varsigma }_1+{\varsigma }_2\frac{B_{10}^{\left(1\right)}}{A_{10}^{\left(1\right)}}+{\varsigma }_3\frac{C_{10}^{\left(1\right)}}{A_{10}^{\left(1\right)}}\right)}{{\gamma }_T}}+{\displaystyle \frac{{\gamma }_0\left(1+3\mu \right){\left(\pi W_m\right)}^2}{4\left(1+\mu \right){\gamma }_T}}+\dots \end{array}$$

(65)

where ${\displaystyle \frac{B_{10}^{\left(1\right)}}{A_{10}^{\left(1\right)}}}, {\displaystyle \frac{C_{10}^{\left(1\right)}}{A_{10}^{\left(1\right)}}}$ are expressed as

$${\displaystyle \frac{B_{10}^{\left(1\right)}}{A_{10}^{\left(1\right)}}}={\displaystyle \frac{\left|\begin{array}{cc}-\left[{\gamma }_4+{\gamma }_{11}\right]& \left[{\gamma }_6+{\gamma }_{14}\right]\\ -\left[{\gamma }_7+{\gamma }_{12}\right]& \left[{\gamma }_9+{\gamma }_{15}\right]\end{array}\right|}{\left|\begin{array}{cc}\left[{\gamma }_5+{\gamma }_{13}\right]& \left[{\gamma }_6+{\gamma }_{14}\right]\\ \left[{\gamma }_8+{\gamma }_{14}\right]& \left[{\gamma }_9+{\gamma }_{15}\right]\end{array}\right|}}$$

(66a)

$${\displaystyle \frac{C_{10}^{\left(1\right)}}{A_{10}^{\left(1\right)}}}={\displaystyle \frac{\left|\begin{array}{cc}\left[{\gamma }_5+{\gamma }_{13}\right]& -\left[{\gamma }_4+{\gamma }_{11}\right]\\ \left[{\gamma }_8+{\gamma }_{14}\right]& -\left[{\gamma }_7+{\gamma }_{12}\right]\end{array}\right|}{\left|\begin{array}{cc}\left[{\gamma }_5+{\gamma }_{13}\right]& \left[{\gamma }_6+{\gamma }_{14}\right]\\ \left[{\gamma }_8+{\gamma }_{14}\right]& \left[{\gamma }_9+{\gamma }_{15}\right]\end{array}\right|}}$$

(66b)

In particular, if the FRC laminated beams are geometrically perfect and symmetrically layered with respect to the geometric mid-plane, there will be no initial deformation caused by initial thermal moments and geometric defects; hence we will obtain

$$\mu =0$$

(67a)

$${\varsigma }_1={\varsigma }_2={\varsigma }_3=0$$

(67b)

and the post-buckling equilibrium path can be expressed as

$${\lambda }_T={\lambda }_T^{\left(0\right)}+\epsilon {\lambda }_T^{\left(1\right)}+{\epsilon }^2{\lambda }_T^{\left(2\right)}+O\left({\epsilon }^4\right)={\displaystyle \frac{{\gamma }_{10}+{\gamma }_1}{{\gamma }_T}}+{\displaystyle \frac{\left({\gamma }_{11}+{\gamma }_2\right)B_{10}^{\left(1\right)}}{{\gamma }_T{A}_{10}^{\left(1\right)}}}+{\displaystyle \frac{\left({\gamma }_{12}+{\gamma }_3\right)C_{10}^{\left(1\right)}}{{\gamma }_T{A}_{10}^{\left(1\right)}}}+{\displaystyle \frac{{\gamma }_0{\left(\pi W_m\right)}^2}{4{\gamma }_T}}+\dots$$

(68)

It's obvious that the post-buckling thermal load-deflection curve in this case is of bifurcation type, and the post-buckling equilibrium path is stable. By letting $W_m=0$ , the critical buckling bifurcation points can be obtained

$${\lambda }_{Tcr}={\displaystyle \frac{{\gamma }_{10}+{\gamma }_1}{{\gamma }_T}}+{\displaystyle \frac{\left({\gamma }_{11}+{\gamma }_2\right)B_{10}^{\left(1\right)}}{{\gamma }_T{A}_{10}^{\left(1\right)}}}+{\displaystyle \frac{\left({\gamma }_{12}+{\gamma }_3\right)C_{10}^{\left(1\right)}}{{\gamma }_T{A}_{10}^{\left(1\right)}}}$$

(69)

4. Numerical Results and Discussion

In this section, numerical results of buckling and post-buckling responses of perfect/imperfect FRC laminated beams with different material properties, slenderness ratio, laying methods and two different boundaries (H-H, C-C) are presented for validation and parameter analysis respectively. The material characteristics are presented as follows [42]

$$\begin{array}{l}{E}_1/E_2=open,E_2=E_3,G_{12}=G_{23}=0.6E_2,G_{23}=0.5E_2\\ {\mu }_{12}={\mu }_{13}=0.25,L/h=open,{\alpha }_2/{\alpha }_1=open\end{array}$$

(70)

4.1. Validation and Comparison

To verify the correctness of the model and method used in this paper, this section ignores geometric defect factors, and calculates dimensionless critical buckling temperatures with different slenderness ratios, elastic modulus ratios and thermal expansion coefficient ratios. The dimensionless buckling temperature is defined by

$${\overline{\Delta T}}_{cr}=\Delta T_{cr}{\alpha }_1{\left({\displaystyle \frac{L}{h}}\right)}^2$$

(71)

Furthermore, the results obtained in this paper are compared with those obtained by the equivalent single-layer theories and other numerical methods in the literature, and the causes of the results difference are analyzed.

Table 1 and Table 2 calculate the critical buckling loads of laminated beams with different slenderness ratios under two different ply modes (0°/90°) and (0°/90°/0°), respectively. Since there is no buckling value of bifurcation type for simply supported asymmetric (0°/90°) laminated beams, the critical buckling values in this case are not presented in table 1, which will be further explained in section 4.2. It is obvious that the results of this paper and those obtained by using the equivalent single-layer theory in the literature are mutually confirmed. For asymmetric (0°/90°) laminated beams, when the slenderness ratio is large, the error between the results in this paper and those in the literature is less than 0.5%, because the transverse shear deformation of slender beams is small and the results of different beam theories are close. However, when the slenderness ratio is small, the discontinuity of interlaminar shear strain cannot be neglected and the results obtained by using the equivalent single layer theory in the literature are overestimated. For symmetric (0°/90°/0°) laminated beams, the results of this paper are very close to those obtained by using the equivalent single-layer theory in literature. Therefore, the improvement of the accuracy of the zig-zag beam model in this paper is relatively limited for symmetrically laminated beams.

In Table 3 and Table 4, the results of this paper are compared with those modeled by the equivalent single-layer high-order shear deformation theories. For asymmetric (0°/90°) laminated beams, when the ratio of elastic modulus $E_1/E_2$ and thermal expansion coefficient ${\alpha }_2/{\alpha }_1$ is small, the errors between the results obtained in this paper and those obtained by the equivalent single layer theory in literature are very small. Conversely, when and ${\alpha }_2/{\alpha }_1$ are relatively large, the error of the results is relatively large, approaching 10%. Comparatively speaking, the ratio of elastic modulus has a much greater effect on the accuracy of the results than the ratio of thermal expansion coefficient. When $E_1/E_2$ is large, the mechanical properties of materials in two orthogonal directions are quite different and the characteristic of the discontinuity of interlaminar strain cannot be neglected. The results obtained by the equivalent single-layer theory will be relatively high. For symmetric (0°/90°/0°) laminated beams, whether simply or fixedly supported, similar results can be obtained that the results of this paper are very close to those obtained by using the equivalent single-layer theory in literature. Therefore, the improvement of the accuracy of the zig-zag beam model in this paper is relatively limited for symmetrically laminated beams.

4.2. Parameter Analysis

In this section, the material properties of Eqs. (70) are still selected, and numerical examples are given to analyze the various factors affecting the post-buckling characteristics of the FRC laminated beams.

Table 5 reveals that under the assumption of in-plane bending, no matter whether the laminated beams are simply supported at both ends or clamped at both ends, the antisymmetric angle-ply laminated beams have the same bifurcated equilibrium path as the symmetric angle-ply laminated beams. This phenomenon can also be verified by Eqs. (13), that is, the equivalent axial tension-compression stiffness and transverse shear stiffness of the two layers with same absolute values ​​of the ply angle are equal. Therefore, we can draw a conclusion that when the beam is symmetrically or anti-symmetrically laminated, there exists the equilibrium path of bifurcation for simply supported laminated beams, while other ply modes will make the geometric and physical neutral surfaces of the laminated beams do not coincide, resulting in tension-bending coupling effect, so that there is no bifurcated equilibrium path.

Figure 2 reflects the effects of elastic modulus ratio and slenderness ratio on the critical buckling temperature of bifurcation type for perfect FRC laminated beams with immovable clamped-clamped and hinged-hinged ends and different ply-angles. From the above, it is found that for simply supported laminated beams, if the absolute values of the ply angles of the two symmetrical layers in geometry are different, there is no critical buckling value of bifurcation type, so the buckling problem in this case is not considered. As can be seen from the figures, the critical buckling value generally increases first and then decreases with the increase of elastic modulus ratio $E_1/E_2$， and the increase of slenderness ratio L/h will delay the decrease of critical buckling value with the increase of elastic modulus ratio $E_1/E_2$. For the given two boundary conditions, fixed support at both ends has higher buckling strength than simple support at both ends. One can also find that for a given number of common ply modes, when the ratio of elasticity modulus is small, [0°/45°/0°] beams have the best ability to resist instability, otherwise [0°/90°/0°] beams are the best, while [45°/-45°/45°] and [45°/-45°] are the worst ply modes to resist buckling.

Figure 3 shows the effects of slenderness ratio on the critical buckling temperature of bifurcation type for perfect symmetrically layered FRC laminated beams with hinged-hinged and clamped-clamped ends. One can find the boundary conditions have greater influence on the buckling strength than the ply angles in the case of given material parameters and that the critical buckling value of beams may decrease sharply with the increase of slenderness ratio.

Figure 4 and Figure 5 study the effect of thermal expansion coefficient ratio on the critical buckling temperature load and post-buckling behaviors for perfect and imperfect FRC laminated beams under different influencing factors. Numerical results suggest that with the increase of thermal expansion coefficient ratio ${\alpha }_2/{\alpha }_1$, the critical buckling temperature decreases and post-buckling load-bearing capacity decreases. Further, it can be found that the smaller the elastic modulus ratio $E_1/E_2$, the faster the buckling value decreases. It can also be found that the influence of elastic modulus ratio on critical buckling load is not only affected by the paving angles and slenderness ratio, but also by thermal expansion coefficient ratio.

Figure 6 compares the post-buckling equilibrium path curves of (45°/-45°/45°) symmetrically laminated beams with different initial imperfection rates. As can be seen from the figure, the laminated beam structures are not sensitive to initial geometric imperfections in general, and the post-buckling equilibrium path is stable. It can be concluded that for perfect symmetrically paved laminated beams, whether simply supported or fixed at both ends, the whole equilibrium path is of bifurcation type. With the increase of defect rate, the bearing capacity of beams decreases at the initial post-buckling stage, but in the deep post-buckling stage, contrary to what one might suppose, the initial defect rate can enhance the bending stiffness of the beam to a certain extent. Note that that equilibrium path of bifurcation type does not exist for laminated beams with initial imperfections.

Figure 7 investigates the influence of elastic modulus ratio $E_1/E_2$ on the post-buckling equilibrium paths of symmetrically laminated beams with two different kinds of boundaries conditions and ply angles. It is obvious that laminated beams with fixed ends have better post-buckling bearing capacity than simply supported beams at both ends, which can also be verified by Figure 5 and Figure 6. However, the influence of elastic modulus ratio on the post-buckling equilibrium path of beams is also affected by the ply angles and slenderness ratio implied by Figure 2, and it is difficult to simply explain how the elastic modulus ratio influences the post-buckling equilibrium path. It can be also found that (0°/90°/0°) beams have higher post-buckling strength than (30°/60°/30°) beams.

Figure 8 studies and compares the post-buckling equilibrium paths of hinged-hinged and clamped-clamped laminated beams with different ply modes (symmetrically layered and asymmetrically layered). It can be obtained from Figure 8(a) that for perfect antisymmetric laminated beams, whether simply supported or clamped at both ends, the post-buckling equilibrium paths of bifurcation type always exists because there is no tension-bending coupling effect in this case. However, for other asymmetric laminated beams, if the both ends are simply supported, post-buckling equilibrium path of bifurcation type does not exist due to the existence of the thermal bending moment caused by the eccentricity of neutral surface. From the two figures, it can be seen that (45°/-45°/45°) and (45°/-45°) beams always have the worst post-buckling bearing capacity among the given ply angles.

5. Conclusions

In this paper, a modified zig-zag beam model is used to refine the modeling of laminated beams. The thermal buckling and post-buckling behaviors of both clamped-clamped and hinged-hinged FRC laminated beams with geometric imperfections are systematically studied by using a two-step perturbation technique. Through detailed verification and parameter analysis, the main conclusions can be obtained as follows:

1. For symmetrical or anti-symmetrical laminated beams, the results of the zig-zag beam model and the equivalent single-layer theory in this paper are similar. For general asymmetrical laminated beams, when the slenderness ratio is small or the elastic modulus is large, the results obtained by using the equivalent single-layer theory will be higher.

2. There exists the equilibrium path of bifurcation type for perfect fixed supported laminated beams with arbitrary ply modes and ply angles and for perfect simply supported laminated beams when the beam is symmetrically or anti-symmetrically laminated.

3. The post-buckling of FRC laminated beams is insensitive to initial imperfections and the equilibrium path is stable. Also, initial geometric imperfections will reduce the load-carrying capacity of beams in the initial post-buckling stage, but in the deep post-buckling stage, the initial defect can enhance the bending stiffness of the beam to a certain extent.

4. The influence of elastic modulus ratio on the buckling and post-buckling of beam structures is also affected by slenderness ratio, thermal expansion coefficient ratio and pavement angle, but in general, the critical buckling value increases first and then decreases with the increase of elastic modulus ratio $E_1/E_2$.

5. The two-step perturbation technique is also applicable to the higher-order deformation zig-zag beam models where the displacement field is described by piecewise continuous functions.