Isogeometric analysis of secondary thermal buckling in composite laminates using Reddy’s higher-order shear deformation theory

Abstract

This study conducts a systematic analysis of the secondary buckling behavior of composite laminated plates under thermal loads based on Reddy’s Higher-Order Shear Deformation Theory (HSDT) and the Isogeometric Analysis (IGA) method. By introducing the von Kármán large-deformation theory and initial geometric imperfections, a mechanical model considering the higher-order shear effects and geometric nonlinearity is established. The Non-Uniform Rational B-Spline (NURBS) basis functions are employed for spatial discretization, effectively achieving the unification of geometrically exact description and high-order continuity. In the model, the in-plane, bending, and initial imperfection stiffness matrices are derived in detail, and the nonlinear equilibrium equations are solved using the Newton-Raphson iterative method. This research proposes a critical load prediction method based on the minimum eigenvalue of the tangent stiffness matrix. By combining the linear interpolation technique, the secondary buckling loads and corresponding modes are accurately captured. Numerical examples verify the effectiveness and accuracy of this method in analyzing the thermal buckling paths, initial imperfection sensitivity, and higher-order mode evolution of composite laminated plates. The results show that initial imperfections significantly reduce the critical buckling load of the structure. Moreover, the combination of IGA and HSDT can accurately characterize the complex deformation behavior of thick plates, providing a theoretical basis for the thermal stability design of composite structures in engineering.

Introduction

In the rapid development process of contemporary engineering technology, composite laminates have been widely applied in numerous fields such as aerospace, automotive, and marine industries, thanks to their remarkable advantages, including high specific strength, high specific modulus, and strong designability. Take the aerospace field as an example, composite laminates are extensively used in critical components of aircraft, like wings and fuselages, to reduce the structural weight, enhance fuel efficiency, and improve flight performance. However, these structures often encounter complex thermal environments during service. For instance, when an aircraft is flying at high speed, frictional heating with the air and solar radiation can cause a sharp increase in the structural temperature. In a thermal environment, due to the anisotropy of composite laminates, there are differences in thermal expansion coefficients in different directions, which in turn generate thermal stresses. When the thermal stress reaches a certain level, thermal buckling of the structure occurs. Traditional research has mainly focused on the primary thermal buckling. In reality, after the primary thermal buckling and with continuous temperature increase, the structure may experience a more complex and more hazardous secondary thermal buckling phenomenon. Currently, the research on the secondary thermal buckling of composite laminates is still insufficient. Conducting in-depth research on this topic is of both urgent practical necessity and significant theoretical importance for ensuring the safety and reliability of engineering structures in thermal environments and optimizing structural design.

To date, a substantial amount of research has been conducted on the thermal buckling of composite laminates. Chen et al.^1–3 systematically investigated the thermal buckling behavior of composite laminates under uniform/non-uniform temperature fields using the finite element method and the Galerkin method. Their research revealed that the critical temperature under clamped boundary conditions is 20%–40% higher than that under simply-supported boundary conditions, and the temperature-dependence of material properties requires iterative solutions. Tauchert et al.,^4,5 based on the Reissner-Mindlin theory, confirmed that the transverse shear effect reduces the critical temperature of moderately thick plates by 15%–30%. The buckling resistance can be improved through the optimization of ply angles and the number of layers. Chang^6,7 developed a high-order displacement theory and pointed out that the transverse normal strain needs to be considered for anti-symmetric angle-ply laminates. The prediction results of this theory deviate from those of the first-order shear theory by up to 25%. Thangaratnam et al. ⁸ and Abramovich ⁹ analyzed and showed that in cross-ply laminates, the modulus ratio (E₁/E₂) and the thermal expansion coefficient ratio (α₂/α₁) dominate the instability modes. Avci et al. ¹⁰ indicated that a 30% increase in the hole size leads to a drastic 35% drop in the critical temperature. In terms of numerical methods, Chen and Liu ¹¹ proposed a Levy-type analytical method for efficient handling of specific boundary conditions. Moradi and Mansouri ¹² applied the differential quadrature method to improve the computational efficiency for arbitrary boundary conditions by 50%. Cetkovic ¹³ combined the Reddy layerwise theory with isogeometric analysis to achieve unified geometric and mechanical modeling. The current research trends focus on global high-order theories^14,15 and simplified high-order shear deformation theories ¹⁶ to balance accuracy and efficiency.

Isogeometric analysis, leveraging its advantages of precise geometric description and high-order continuity, has emerged as a crucial tool in the research on the thermal buckling of composite materials. Mirzaei and Kiani ¹⁷ within the framework of IGA based on NURBS and combining it with the first-order shear theory, investigated the thermal buckling behavior of graphene-reinforced laminates. Their research revealed the significant impacts of the functionally graded distribution and temperature-dependent material properties. Tran et al. ¹⁸ proposed a quasi-three-dimensional six-variable model. By satisfying the C¹ continuity requirement through IGA, they confirmed the crucial role of the transverse normal strain in the thermal buckling of thick plates. Yu et al. ¹⁹ and Farzam and Hassani ²⁰ developed an approach that couples IGA with the modified couple stress theory for functionally graded plates and carbon nanotube-reinforced plates, respectively. This approach resolves the shear locking problem under thermo-mechanical coupled loads and quantifies the influences of parameters such as gradient indices and boundary conditions. Do et al.^21–25 systematically advanced the high-order theories of IGA. The nth order shear deformation theory optimizes the prediction of the nonlinear buckling of temperature-dependent functionally graded material plates. The quasi-three-dimensional IGA model characterizes the thickness expansion effect using only four variables. For the post-buckling analysis of graphene-reinforced plates, the complete instability path is traced by combining the modified Newton iteration method. Mohammadi ²⁶ extended IGA to folded plate structures. By using the bending strip method to handle multi-patch connections, the control mechanism of the graphene distribution pattern on the buckling temperature was verified. Research indicates that IGA, through the high-order continuity of NURBS basis functions, significantly enhances the accuracy of buckling prediction under thermal gradients, material nonlinearities, and complex geometries, thus providing a powerful tool for the design of advanced composite material structures.

Secondary buckling, as a high-order instability phenomenon in post-buckling behavior, has witnessed crucial breakthroughs in the research of various types of structures regarding its mechanism and prediction methods. Chien ²⁷ pioneered the use of the continuation method combined with the local perturbation technique. Numerically, they revealed that the secondary buckling of thin plates exhibits asymmetric wrinkling characteristics. Wu^28,29 based on the Lyapunov-Schmidt reduction, proved that an elastic strut with a rectangular cross-section will inevitably trigger secondary bifurcation under the critical section ratio. They also developed the augmented system method to directly solve for the double buckling loads and modes, avoiding the main path tracking and thus enhancing the computational efficiency by more than 40%. Hunt and Everall^30,31 introduced the Arnold tongue theory to quantify the modal competition, revealing that compound bifurcations reduce the post-buckling reserve by 30%–50%, and the stability margin of long plates is twice as high as that of short plates. In the field of composite materials, Tiwari and Hyer^32,33 verified through the first-order shear finite element method and experiments that a 20% increase in the normal displacement constraint on the unloaded edge leads to a 15% decrease in the critical load. Zhong et al.^34,35 established a strict geometric nonlinear theory to capture the modal jumps of hygrothermal laminates, where the displacement amplitude at the bifurcation point is three times that at the primary buckling point. Tao and Dai ³⁶ based on the meshless RPIM method, confirmed that the sensitivity of the secondary instability of porous S-FGM plates to simply-supported boundaries is 40% higher than that to fixed-supported boundaries. Li and Huang ³⁷ revealed that when the ratio of the shear strength of the core material to the tensile strength of the face sheet in sandwich plates is lower than the threshold of 1.8, a chain reaction of secondary buckling is triggered. Regarding the mechanism of special structures, Knoche and Kierfeld ³⁸ derived the critical conditions for polygonal wrinkling of spherical shells through the Pogorelov model, proving that the critical volume reduction rate is linearly related to the bending stiffness. Parry et al., ³⁹ based on the elastic hinged rod-contact model, predicted that the critical aspect ratio for the “ruled-bubble” transition of thin films is 2.5 ± 0.3. Zhang and Gu ⁴⁰ used the asymmetric grid perturbation method to confirm that the fluctuation of the secondary buckling load of composite cylindrical shells is less than 2%, establishing it as a core design criterion. In terms of experimental verification, Alvarez and Bisagni^41,42 developed a low-thermal-expansion fixture and successfully captured the thermally induced modal jumps of laminated plates. Chen et al.^43,44 found that the modal migration frequency of thermally-buckled plates reaches 70% of the fundamental frequency, and the energy of the anti-symmetric mode surges by 30% after the jump, thus destroying the stability. Murakami and Yao ⁴⁵ quantified the snap-through effect of rectangular plates: an initial imperfection in primary buckling increases the snap-through load by 25%, while an initial imperfection in secondary buckling reduces it by 20%.

The thermoplastic matrix endows the composite material with unique plastic deformation ability, significantly enhancing the post-buckling damage tolerance. Shah et al. ⁴⁶ revealed through post-impact compression tests that the reduction in the buckling load of thermoplastic 3D woven composites (44.5%) is far lower than that of thermosetting epoxy-based composites (84.5%). This can be attributed to the local dissipation of impact energy by the plastic deformation of the matrix and better interfacial adhesion. Gaudenzi further confirmed that under pure shear loading, thermoplastic laminates can delay the propagation of buckling through the plastic flow of the matrix. ⁴⁷ The error between the non-linear finite element model and the experimental results is less than 8%. In the study of structures with holes, Yapici ⁴⁸ and Yazici ⁴⁹ found that the sensitivity of the critical buckling load of metal-reinforced thermoplastic laminates to holes increases with the enhancement of boundary constraints. A square hole under fixed-support boundaries leads to a 35% decrease in the load, while under simply-supported boundaries, the decrease is only 15%. Regarding the prediction of buckling in a thermal environment, Barton ⁵⁰ proposed a closed-form approximation method to solve the thermal buckling temperature of symmetrically laminated plates with one-side fixed support. The computational efficiency is 50% higher than that of the Rayleigh-Ritz method. Hieu and Tung systematically developed the theory of thermal post-buckling for carbon nanotube-reinforced structures. For a cylindrical shell with a gradient distribution of CNTs on an elastic foundation, the critical temperature increases by 40%. ⁵¹ Tangential constrained edges can increase the post-buckling bearing capacity of sandwich panels by 30%.^52,53 A non-uniform temperature field (sinusoidal/linear distribution) induces local mode jumps. ⁵⁴ Kumar ⁵⁵ combined the Reddy higher-order theory with the stochastic perturbation method and found that a 30% increase in the diameter of the circular hole in a functionally graded plate causes a sudden 35% drop in the critical temperature. Nam/Phuong et al.^56,57 innovatively designed a spider-web-stiffened graphene-reinforced porous spherical shell. The thermo-mechanical coupling buckling load of this shell is 2.1 times higher than that of the non-stiffened structure. Wan et al. ⁵⁸ conducted a study on the sealing performance of the buckling joints of glass-fiber-reinforced thermoplastic pipes. Finite element analysis shows that a 5-mm buckling amount optimizes the distribution of sealing stress, but a buckling amount greater than 8 mm will cause excessive plastic deformation of the fibers, leading to failure. The critical design parameters are a tooth height of 3 mm and a combination of tooth angles of 55°/25°, which can ensure safe operation under a working pressure of 6.4 MPa. The results of this paper can provide a reference for the design of thermal buckling loads of thermoplastic composites.

From the above research, it can be observed that although scholars have a thorough understanding of the primary thermal buckling of composite laminates, their knowledge of the secondary thermal buckling is relatively limited. This limitation is insufficient for fully exploiting the load-bearing capacity of composite laminates. Therefore, in this paper, based on Reddy’s HSDT, the von Kármán large-deformation theory and initial geometric imperfections are introduced to establish a mechanical model that takes into account both higher-order shear effects and geometric nonlinearity. Spatial discretization is carried out using NURBS basis functions. The effectiveness of the method proposed in this paper is verified by comparing the calculation results with the experimental results and finite element analysis results reported in Refs. 41 and 42. Finally, through parametric analysis, the influences of initial geometric imperfections, aspect ratios, ply angles, elastic moduli, and thermal expansion coefficients on the secondary critical buckling temperature under different boundary conditions are thoroughly investigated.

Theory and formulas

Problem statement

Consider a composite laminate of the N_L layer with a length of a, a width of b, and a thickness of h. The ply-stacking angles are θ _k (k = 1, 2, …, N_L). As shown in Figure 1, this laminate is placed in a uniform thermal environment with a temperature change of ΔT. A Cartesian coordinate system O-xyz is established, where the xy-plane lies on the geometric mid-surface of the laminate.

Figure 1.

Schematic diagram of the composite laminate.

Basic equations

Based on Reddy’s HSDT, the displacement field is assumed as follows:

$\begin{array}{l}u\left(x,y,z\right)=u_0\left(x,y\right)+z{\theta }_x\left(x,y\right)-\frac{4z^3}{3h^2}\left({\theta }_x+\frac{\partial w}{\partial x}\right),\\ v\left(x,y,z\right)=v_0\left(x,y\right)+z{\theta }_y\left(x,y\right)-\frac{4z^3}{3h^2}\left({\theta }_y+\frac{\partial w}{\partial y}\right),\\ w\left(x,y,z\right)=w_0\left(x,y\right)+\mathrm{\Delta }{w}_0\left(x,y\right),\end{array}$

(1)

Here, u₀, v₀, and w₀ represent the displacement components of the mid-surface. θ _x and θ_y denote the transverse shear rotation angles. z is the coordinate along the thickness direction, with its range being [−h/2, h/2]. Δw_0​ represents the initial geometric imperfection, which is characterized by a sine function.

$\mathrm{\Delta }w=\phi h\hspace{0.17em}\sin \left(\frac{\pi x}{a}\right)\sin \left(\frac{\pi y}{b}\right)$

(2)

Here, φ represents the amplitude of the initial imperfection.

Taking into account the von Kármán large-deformation theory and initial imperfections, the total strain is divided into linear and nonlinear parts.

$\begin{array}{l}{\epsilon }_x=\frac{\partial u}{\partial x}+\frac{1}{2}{\left(\frac{\partial w}{\partial x}\right)}^2-\frac{1}{2}{\left(\frac{\partial \mathrm{\Delta }{w}_0}{\partial x}\right)}^2-z\frac{{\partial }^2\left(w-\mathrm{\Delta }{w}_0\right)}{\partial x^2},\\ {\epsilon }_y=\frac{\partial v}{\partial y}+\frac{1}{2}{\left(\frac{\partial w}{\partial y}\right)}^2-\frac{1}{2}{\left(\frac{\partial \mathrm{\Delta }{w}_0}{\partial y}\right)}^2-z\frac{{\partial }^2\left(w-\mathrm{\Delta }{w}_0\right)}{\partial y^2},\\ {\gamma }_{xy}=\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}+\frac{\partial w}{\partial x}\frac{\partial w}{\partial y}-\frac{\partial \mathrm{\Delta }{w}_0}{\partial x}\frac{\partial \mathrm{\Delta }{w}_0}{\partial y}-2z\frac{{\partial }^2\left(w-\mathrm{\Delta }{w}_0\right)}{\partial x\partial y},\\ {\gamma }_{xz}=\frac{\partial u}{\partial z}+\frac{\partial w}{\partial x}={\theta }_x-\frac{4z^2}{h^2}\left({\theta }_x+\frac{\partial w}{\partial x}\right)+\frac{\partial \mathrm{\Delta }{w}_0}{\partial x},\\ {\gamma }_{yz}=\frac{\partial v}{\partial z}+\frac{\partial w}{\partial y}={\theta }_y-\frac{4z^2}{h^2}\left({\theta }_y+\frac{\partial w}{\partial y}\right)+\frac{\partial \mathrm{\Delta }{w}_0}{\partial y},\end{array}$

(3)

Where ε _x and ε _y are the in-plane normal strains, γ _xy is the in-plane shear strain, and γ _xy and γ _xy are the transverse shear strains. $\frac{\partial \mathrm{\Delta }{w}_0}{\partial x}$ , $\frac{\partial \mathrm{\Delta }{w}_0}{\partial y}$ represents the gradient term of the initial imperfection, which is used to modify the geometric nonlinear strain.

For a composite single-layer plate, the constitutive relationship in the material principal axis (1-2 coordinate system) is as follows:

$\left[\begin{array}{l}{\sigma }_1\\ {\sigma }_2\\ {\tau }_{12}\\ {\tau }_{13}\\ {\tau }_{23}\end{array}\right]=\left[\begin{array}{lllll}{Q}_{11}& Q_{12}& 0& 0& 0\\ Q_{12}& Q_{22}& 0& 0& 0\\ 0& 0& Q_{66}& 0& 0\\ 0& 0& 0& Q_{44}& 0\\ 0& 0& 0& 0& Q_{55}\end{array}\right]\left(\left[\begin{array}{l}{\epsilon }_1\\ {\epsilon }_2\\ {\gamma }_{12}\\ {\gamma }_{13}\\ {\gamma }_{23}\end{array}\right]-\left[\begin{array}{l}{\alpha }_1\\ {\alpha }_2\\ 0\\ 0\\ 0\end{array}\right]\mathrm{\Delta }T\right)$

(4)

Among them, ${\sigma }_1,{\sigma }_2,{\tau }_{12},{\tau }_{13},{\tau }_{23}$ are the stress components in the material principal axis; ${\epsilon }_1,{\epsilon }_2,{\gamma }_{12},{\gamma }_{13},{\gamma }_{23}$ are the strain components in the material principal axis. Q _ij represents the elastic constants in the material principal axis:

$\begin{array}{l}{Q}_{11}=\frac{E_1}{1-{\nu }_{12}{\nu }_{21}},Q_{22}=\frac{E_2}{1-{\nu }_{12}{\nu }_{21}},Q_{12}=\frac{{\nu }_{12}{E}_2}{1-{\nu }_{12}{\nu }_{21}},\\ Q_{66}=G_{12},Q_{44}=G_{13},Q_{55}=G_{23}.\end{array}$

(5)

And α₁ and α₂ are the thermal expansion coefficients in the directions of the material principal axes. $\mathrm{\Delta }T=T-T_0$ represents the temperature change, where T₀ is the reference temperature.

Rotate the material principal axes (1-2) by an angle θ to the global coordinate system (x-y), and define the transformation matrix T:

$\mathbf{T}\left(\theta \right)=\left[\begin{array}{lllll}{\cos }^2\hspace{0.17em}\theta & {\sin }^2\hspace{0.17em}\theta & \sin \hspace{0.17em}\theta \hspace{0.17em}\cos \hspace{0.17em}\theta & 0& 0\\ {\sin }^2\hspace{0.17em}\theta & {\cos }^2\hspace{0.17em}\theta & -\sin \hspace{0.17em}\theta \hspace{0.17em}\cos \hspace{0.17em}\theta & 0& 0\\ -2\hspace{0.17em}\sin \hspace{0.17em}\theta \hspace{0.17em}\cos \hspace{0.17em}\theta & 2\hspace{0.17em}\sin \hspace{0.17em}\theta \hspace{0.17em}\cos \hspace{0.17em}\theta & {\cos }^2\hspace{0.17em}\theta -{\sin }^2\hspace{0.17em}\theta & 0& 0\\ 0& 0& 0& \cos \hspace{0.17em}\theta & \sin \hspace{0.17em}\theta \\ 0& 0& 0& -\sin \hspace{0.17em}\theta & \cos \hspace{0.17em}\theta \end{array}\right]$

(6)

The relationship between the elastic matrix in the global coordinate system and the elastic matrix in the material principal axes is as follows:

${\overline{\mathbf{Q}}}_{ij}={\mathbf{T}}^{-1}\left(\theta \right){\mathbf{Q}}_{ij}\mathbf{T}\left(\theta \right)$

(7)

The constitutive equation in the global coordinate system is as follows:

$\left[\begin{array}{l}{\sigma }_x\\ {\sigma }_y\\ {\tau }_{xy}\\ {\tau }_{xz}\\ {\tau }_{yz}\end{array}\right]=\left[\begin{array}{lllll}{\overline{Q}}_{11}& {\overline{Q}}_{12}& {\overline{Q}}_{16}& 0& 0\\ {\overline{Q}}_{12}& {\overline{Q}}_{22}& {\overline{Q}}_{26}& 0& 0\\ {\overline{Q}}_{16}& {\overline{Q}}_{26}& {\overline{Q}}_{66}& 0& 0\\ 0& 0& 0& {\overline{Q}}_{44}& {\overline{Q}}_{45}\\ 0& 0& 0& {\overline{Q}}_{45}& {\overline{Q}}_{55}\end{array}\right]\left(\left[\begin{array}{l}{\epsilon }_x\\ {\epsilon }_y\\ {\gamma }_{xy}\\ {\gamma }_{xz}\\ {\gamma }_{yz}\end{array}\right]-\left[\begin{array}{l}{\alpha }_x\\ {\alpha }_y\\ {\alpha }_{xy}\\ 0\\ 0\end{array}\right]\mathrm{\Delta }T\right)$

(8)

A brief of the NURBS basis functions

In the concept of IGA, one-dimensional basis functions are typically constructed using an open knot vector $\mathrm{\Xi }=\left\{{\xi }_1,{\xi }_2,{\xi }_3,\dots {\xi }_{i,\dots ,}{\xi }_{n+p+1}\right\}$ , where ${\xi }_i\le {\xi }_{i+1}$ . The multiplicities of the first knot ξ₁ and the last knot ξ_n+p+1 are p+1, where p and n denote the order and the number of basis functions, respectively. The parameter-space interval [ξ₁, ξ_n+p+1] is referred to as a patch, which serves to map the physical domain to the parameter space. Meanwhile, the interval [ξ _i , ξ_i+1) represents a knot span, which corresponds to the parameter space of an element. The B-spline basis functions associated with a given knot vector are constructed in a recursive manner.

$\begin{array}{l}{N}_{i,0}\left(\xi \right)=\{\begin{array}{ll}1& {\xi }_i\le \xi <{\xi }_{i+1}\\ 0& \text{otherwise}\end{array}\text{for}\hspace{0.17em}p=0,\\ N_{i,p}=\frac{\xi -{\xi }_i}{{\xi }_{i+p}-{\xi }_i}{N}_{i,p-1}\left(\xi \right)+\frac{{\xi }_{i+p+1}-\xi }{{\xi }_{i+p+1}-{\xi }_{i+1}}{N}_{i,p-1}\left(\xi \right)\hspace{0.17em}\text{for}\hspace{0.17em}p\ge 0.\end{array}$

(9)

Here, i = 1, 2, …, n + p + 1, and the value of 0/0 is assumed to be 0. NURBS basis functions are introduced as a set of B-spline basis functions with associated weight factors. The mathematical definition of one-dimensional NURBS basis functions is as follows:

$R_i^p\left(\xi \right)=\frac{N_{i,p}\left(\xi \right)w_i}{{\sum }_{\widehat{i}=1}^n{N}_{\widehat{i},p}\left(\xi \right)w_{\widehat{i}}}$

(10)

Where w _i > 0 are the weights. If all the weights are set to be the same, the NURBS basis functions degenerate into B-spline basis functions. Two-dimensional NURBS basis functions are obtained by taking the tensor product of two one-dimensional NURBS functions, as shown below:

$R_{i,j}^{p,q}\left(\xi ,\eta \right)=\frac{N_{i,p}\left(\xi \right)M_{j,q}\left(\eta \right)w_{i,j}}{{\sum }_{\widehat{i}=1}^n{\sum }_{\widehat{j}=1}^m{N}_{\widehat{i},p}\left(\xi \right)M_{\widehat{j},q}\left(\eta \right)w_{\widehat{i},\widehat{j}}}$

(11)

Figure 2 illustrates the 2D cubic B-spline basis functions with the open uniform knot vectors $\mathrm{\Xi }=\left\{0,0,0,0,0.2,0.4,0.6,0.8,1,1,1,1\right\}$ and $\mathrm{{\rm H}}=\left\{0,0,0,0,0.2,0.4,0.6,0.8,1,1,1,1\right\}$ . For two-dimensional problems, two-dimensional NURBS surfaces are constructed by performing the following linear combination of NURBS functions and their associated control points:

${\mathrm{\Omega }}_{2D}={\sum }_{i=1}^n{\sum }_{j=1}^m{R}_{i,j}^{p,q}\left(\xi ,\eta \right){\mathbf{B}}_{i,j}$

(12)

Figure 2.

2D cubic B-spline basis functions.

The above equation is the result of the tensor product given the control nets B and knot vectors $\mathrm{\Xi }=\left\{{\xi }_1,{\xi }_2,\dots ,{\xi }_{n+p+1}\right\}$ and $\mathbf{H}=\left\{{\eta }_1,{\eta }_2,\dots ,{\eta }_{m+q+1}\right\}$ .

NURBS discretization

In the physical space, the coordinates (x, y) are defined by the weighted sum of the coordinates of the control points P_i,j=(X_i,j, Y_i,j) and the NURBS basis functions $R_{i,j}^{p,q}\left(\xi ,\eta \right)$ , as follows:

$x\left(\xi ,\eta \right)={\displaystyle \sum }_{i=1}^n{\displaystyle \sum }_{j=1}^m{R}_{i,j}^{p,q}\left(\xi ,\eta \right)X_{i,j},y\left(\xi ,\eta \right)={\displaystyle \sum }_{i=1}^n{\displaystyle \sum }_{j=1}^m{R}_{i,j}^{p,q}\left(\xi ,\eta \right)Y_{i,j}$

(13)

Based on Reddy’s HSDT, the displacement field is discretized using NURBS basis functions.

$\begin{array}{l}\begin{array}{l}{u}_0\left(\xi ,\eta \right)={\sum }_{i,j}{R}_{i,j}{U}_{i,j},\\ v_0\left(\xi ,\eta \right)={\sum }_{i,j}{R}_{i,j}{V}_{i,j},\\ w_0\left(\xi ,\eta \right)={\sum }_{i,j}{R}_{i,j}{W}_{i,j},\\ {\theta }_x\left(\xi ,\eta \right)={\sum }_{i,j}{R}_{i,j}{\mathrm{\Theta }}_{x,i,j},\\ {\theta }_y\left(\xi ,\eta \right)={\sum }_{i,j}{R}_{i,j}{\mathrm{\Theta }}_{y,i,j},\end{array}\\ \mathrm{\Delta }w\left(x,y\right)={\sum }_{i,j}{R}_{ij}\left(u,v\right)\mathrm{\Delta }{\mathbf{P}}_{ij},\end{array}$

(14)

Each control point ${\mathbf{P}}_{i,j}$ corresponds to the degrees of freedom ${\mathbf{u}}_{i,j}={\left[U_{i,j},V_{i,j},W_{i,j},{\mathrm{\Theta }}_{x,i,j},{\mathrm{\Theta }}_{y,i,j}\right]}^{\mathrm{T}}$ .

Then, the control points considering the initial imperfection is:

${\mathbf{P}}_{ij}^{\text{imperfection}}={\mathbf{P}}_{ij}+\mathrm{\Delta }{\mathbf{P}}_{ij}$

(15)

The total potential energy functional Π consists of the strain energy U and the potential energy of external forces W_ext:

$\mathrm{\Pi }=U-W_{\text{ext}}$

(16)

The strain energy U includes linear and nonlinear strain energies, corresponding to the contributions of each stiffness matrix. The potential energy of external forces W_ext is the work done by thermal loads on the displacement field.

The expression for the strain energy is as follows:

$U=\frac{1}{2}{\int }_V\mathbf{\sigma }\hspace{0.17em}\hspace{0.17em}\mathbf{\epsilon }\mathrm{d}V$

(17)

After substituting the stress and strain into equation (16) and rearranging the terms, we obtain:

$\mathrm{\Pi }=\frac{1}{2}{\mathbf{u}}^{\mathrm{T}}\left({\mathbf{K}}_{\mathrm{M}}+{\mathbf{K}}_{\mathrm{B}}+{\mathbf{K}}_{\mathrm{I}}\right)\mathbf{u}+\frac{1}{2}{\mathbf{u}}^{\mathrm{T}}{\mathbf{N}}_1\mathbf{u}+\frac{1}{6}{\mathbf{u}}^{\mathrm{T}}{\mathbf{N}}_2{\mathbf{u}}^3-{\mathbf{u}}^{\mathrm{T}}{\mathbf{F}}_{\text{ext}}$

(18)

where the complete form of ${\mathbf{u}}^{\mathrm{T}}{\mathbf{N}}_2{\mathbf{u}}^3$ is ${\mathbf{u}}^{\mathrm{T}}\left({\mathbf{N}}_2\otimes \mathbf{u}\otimes \mathbf{u}\right)\mathbf{u}$ .

The expressions for each matrix in equation (17) are as follows:

In-plane stiffness matrix K_M.

${\mathbf{K}}_{\mathrm{M}}={\displaystyle \sum }_{k=1}^{N_{\mathrm{L}}}{\int }_{\mathrm{\Omega }}{\mathbf{B}}_{\mathrm{m}}^{\mathrm{T}}{\overline{\mathbf{Q}}}_{\mathrm{m}}^{\left(k\right)}{\mathbf{B}}_{\mathrm{m}}\left|\mathbf{J}\right|\mathrm{d}\xi \mathrm{d}\eta ·{\int }_{z_k}^{z_{k+1}}\left(1,z,z^3\right)\mathrm{d}z$

(19)

Here, B_m is the in-plane strain-displacement matrix, which contains the first-order derivatives of the NURBS basis functions.

${\mathbf{B}}_{\mathrm{m}}=\left[\begin{array}{lllll}\frac{\partial R_i}{\partial x}& 0& -z\frac{{\partial }^2{R}_i}{\partial x^2}& z^3\frac{\partial R_i}{\partial x}& \cdots \\ 0& \frac{\partial R_i}{\partial y}& -z\frac{{\partial }^2{R}_i}{\partial y^2}& z^3\frac{\partial R_i}{\partial y}& \cdots \\ \frac{\partial R_i}{\partial y}& \frac{\partial R_i}{\partial x}& -2z\frac{{\partial }^2{R}_i}{\partial x\partial y}& z^3\left(\frac{\partial R_i}{\partial y}+\frac{\partial R_i}{\partial x}\right)& \cdots \end{array}\right]$

(20)

And ${\overline{\mathbf{Q}}}_{\mathrm{m}}^{\left(k\right)}$ is the in-plane elastic matrix of the kth layer.

The integral term ${\int }_{-h/2}^{h/2}\left(1,z,z^3\right)\mathrm{d}z$ corresponds to the layer-wise contributions of in-plane tension, bending, and high-order shear effects.

Each element of the Jacobian matrix J is the partial derivative of the physical coordinates with respect to the parametric coordinates.

$\mathbf{J}=\left[\begin{array}{ll}\frac{\partial x}{\partial \xi }& \frac{\partial x}{\partial \eta }\\ \frac{\partial y}{\partial \xi }& \frac{\partial y}{\partial \eta }\end{array}\right]$

(21)

where

$\begin{array}{l}\frac{\partial x}{\partial \xi }={\displaystyle \sum }_{i,j}\frac{\partial R_{i,j}}{\partial \xi }{X}_{i,j},\frac{\partial x}{\partial \eta }={\displaystyle \sum }_{i,j}\frac{\partial R_{i,j}}{\partial \eta }{X}_{i,j},\\ \frac{\partial y}{\partial \xi }={\displaystyle \sum }_{i,j}\frac{\partial R_{i,j}}{\partial \xi }{Y}_{i,j},\frac{\partial y}{\partial \eta }={\displaystyle \sum }_{i,j}\frac{\partial R_{i,j}}{\partial \eta }{Y}_{i,j}.\end{array}$

(22)

Bending stiffness matrix K_B.

${\mathbf{K}}_{\mathrm{B}}={\displaystyle \sum }_{k=1}^{N_{\mathrm{L}}}{\int }_{\mathrm{\Omega }}{\mathbf{B}}_{\mathrm{b}}^{\mathrm{T}}{\overline{\mathbf{Q}}}_{\mathrm{b}}^{\left(k\right)}{\mathbf{B}}_{\mathrm{b}}\left|\mathbf{J}\right|\mathrm{d}\xi \mathrm{d}\eta ·{\int }_{z_k}^{z_{k+1}}\left(z^2,z^4,z^6\right)\mathrm{d}z$

(23)

Among them, B_b is the bending strain-displacement matrix, which contains the second-order derivatives of NURBS basis functions.

${\mathbf{B}}_{\mathrm{b}}=\left[\begin{array}{lllll}0& 0& -\frac{{\partial }^2{R}_i}{\partial x^2}& 3z^2\frac{\partial R_i}{\partial x}& \cdots \\ 0& 0& -\frac{{\partial }^2{R}_i}{\partial y^2}& 3z^2\frac{\partial R_i}{\partial y}& \cdots \\ 0& 0& -2\frac{{\partial }^2{R}_i}{\partial x\partial y}& 3z^2\left(\frac{\partial R_i}{\partial y}+\frac{\partial R_i}{\partial x}\right)& \cdots \end{array}\right]$

(24)

And ${\overline{\mathbf{Q}}}_{\mathrm{b}}^{\left(k\right)}$ is the bending elastic matrix, which is isomorphic to the in-plane matrix but has different integration weights.

${\overline{\mathbf{Q}}}_{\mathrm{b}}^{\left(k\right)}={\overline{\mathbf{Q}}}_{\mathrm{m}}^{\left(k\right)}·{\int }_{-h/2}^{h/2}{z}^n\mathrm{d}z,\left(n=2,4,6\right)$

(25)

The integral term ${\int }_{-h/2}^{h/2}\left(z^2,z^4,z^6\right)\mathrm{d}z$ reflects the high-order shear correction of the bending stiffness.

The initial imperfection stiffness matrix K_I.

${\mathbf{K}}_{\mathrm{I}}={\displaystyle \sum }_{k=1}^{N_{\mathrm{L}}}{\int }_{\mathrm{\Omega }}{\mathbf{B}}_{\mathrm{g}}^{\mathrm{T}}{\mathbf{N}}_0{\mathbf{B}}_{\mathrm{g}}\left|\mathbf{J}\right|\mathrm{d}\xi \mathrm{d}\eta \mathrm{d}z$

(26)

Here, B_g is the geometric strain matrix, which includes the initial imperfection gradient.

${\mathbf{B}}_{\mathrm{g}}=\left[\begin{array}{llll}\frac{\partial R_i}{\partial x}& 0& \frac{\partial \mathrm{\Delta }{w}_0}{\partial x}\frac{\partial R_i}{\partial x}& \dots \\ 0& \frac{\partial R_i}{\partial y}& \frac{\partial \mathrm{\Delta }{w}_0}{\partial y}\frac{\partial R_i}{\partial y}& \dots \\ \frac{\partial R_i}{\partial y}& \frac{\partial R_i}{\partial x}& \left(\frac{\partial \mathrm{\Delta }{w}_0}{\partial x}\frac{\partial R_i}{\partial y}+\frac{\partial \mathrm{\Delta }{w}_0}{\partial y}\frac{\partial R_i}{\partial x}\right)& \dots \end{array}\right]$

(27)

And N₀ is the initial stress matrix.

${\mathbf{N}}_0=\left[\begin{array}{ll}{N}_{0x}& N_{0xy}\\ N_{0xy}& N_{0y}\end{array}\right]={\displaystyle \sum }_{k=1}^N{\int }_{z_k}^{z_{k+1}}\left[\begin{array}{ll}{\sigma }_{0x}& {\sigma }_{0xy}\\ {\sigma }_{0xy}& {\sigma }_{0y}\end{array}\right]\mathrm{d}z$

(28)

Among them, σ₀ is the prestress caused by initial imperfections.

Geometric nonlinear terms N₁ and N₂.

$\begin{array}{l}{\mathbf{N}}_1={\displaystyle \sum }_{k=1}^{N_{\mathrm{L}}}{\int }_{\mathrm{\Omega }}{\mathbf{G}}^{\mathrm{T}}{\mathbf{S}}_{\mathrm{L}}\mathbf{G}\left|\mathbf{J}\right|\mathrm{d}\xi \mathrm{d}\eta \mathrm{d}z,\\ {\mathbf{N}}_2={\displaystyle \sum }_{k=1}^{N_{\mathrm{L}}}{\int }_{\mathrm{\Omega }}{\mathbf{G}}^{\mathrm{T}}{\mathbf{S}}_{\mathrm{Q}}\mathbf{G}\left|\mathbf{J}\right|\mathrm{d}\xi \mathrm{d}\eta \mathrm{d}z.\end{array}$

(29)

Here, G is the displacement gradient matrix, which contains the first-order derivatives of NURBS basis functions.

$\mathbf{G}=\left[\begin{array}{llll}\frac{\partial R_i}{\partial x}& 0& \frac{\partial w}{\partial x}\frac{\partial R_i}{\partial x}& \cdots \\ 0& \frac{\partial R_i}{\partial y}& \frac{\partial w}{\partial y}\frac{\partial R_i}{\partial y}& \cdots \\ \frac{\partial R_i}{\partial y}& \frac{\partial R_i}{\partial x}& \left(\frac{\partial w}{\partial x}\frac{\partial R_i}{\partial y}+\frac{\partial w}{\partial y}\frac{\partial R_i}{\partial x}\right)& \cdots \end{array}\right]$

(30)

And S_L is the linear stress matrix, which is composed of the current in-plane stresses.

${\mathbf{S}}_{\mathrm{L}}=\left[\begin{array}{lll}{\sigma }_x& {\tau }_{xy}& 0\\ {\tau }_{xy}& {\sigma }_y& 0\\ 0& 0& 0\end{array}\right]$

(31)

And S_Q is the quadratic stress matrix, which reflects the nonlinear coupling effects.

${\mathbf{S}}_{\mathrm{Q}}=\left[\begin{array}{lll}\frac{1}{2}{\sigma }_x& \frac{1}{2}{\tau }_{xy}& 0\\ \frac{1}{2}{\tau }_{xy}& \frac{1}{2}{\sigma }_y& 0\\ 0& 0& 0\end{array}\right]$

(32)

The expression for the equivalent nodal forces of thermal loads is:

${\mathbf{F}}_{\text{ext}}={\mathbf{F}}_{\text{thermal}}={\displaystyle \sum }_{k=1}^{N_{\mathrm{L}}}{\int }_{\mathrm{\Omega }}{\mathbf{B}}^{\mathrm{T}}{\mathbf{Q}}^{\left(k\right)}{\mathbf{\alpha }}^{\left(k\right)}\mathrm{\Delta }T\left|\mathbf{J}\right|\mathrm{d}\xi \mathrm{d}\eta .{\int }_{z_k}^{{z_k}_{+1}}\mathrm{d}z$

(33)

The strain-displacement matrix B contains the first-order derivatives of NURBS basis functions, which converts nodal displacements into strains.

$\mathbf{B}=\left[\begin{array}{llll}\frac{\partial R_i}{\partial x}& 0& -z\frac{{\partial }^2{R}_i}{\partial x^2}& \cdots \\ 0& \frac{\partial R_i}{\partial y}& -z\frac{{\partial }^2{R}_i}{\partial y^2}& \cdots \\ \frac{\partial R_i}{\partial y}& \frac{\partial R_i}{\partial x}& -2z\frac{{\partial }^2{R}_i}{\partial x\partial y}& \cdots \end{array}\right]$

(34)

The nonlinear equilibrium equations are derived through the variational principle δΠ = 0.

$\left({\mathbf{K}}_{\mathrm{M}}+{\mathbf{K}}_{\mathrm{B}}+{\mathbf{K}}_{\mathrm{I}}+{\mathbf{N}}_1+{\mathbf{N}}_2{\mathbf{u}}^2\right)\mathbf{u}={\mathbf{F}}_{\text{ext}}$

(35)

Within the framework of this research paper, the Newton-Raphson iterative algorithm is utilized to tackle the nonlinear equilibrium equations, the solution process is illustrated in Figure 3. The subsequent content details the specific sequential steps of this solution-seeking process:

(1) Initialization: Given the initial displacement u₀ = 0, set the convergence tolerance ϵ.

(2) Iterative loop (step k):

(a) Calculate the internal forces: ${\mathbf{F}}_{\mathrm{int}}^{\left(k\right)}=\left({\mathbf{K}}_M+{\mathbf{K}}_B+{\mathbf{K}}_I+{\mathbf{N}}_1+{\mathbf{N}}_2{{\mathbf{u}}^2}^{\left(k\right)}\right){\mathbf{u}}^{\left(k\right)},$

(b) Calculate the residual: ${\mathbf{R}}^{\left(k\right)}={\mathbf{F}}_{\text{ext}}-{\mathbf{F}}_{\mathrm{int}}^{\left(k\right)}$ ,

(c) Construct the tangent stiffness matrix: ${\mathbf{K}}_T^{\left(k\right)}={\mathbf{K}}_M+{\mathbf{K}}_B+{\mathbf{K}}_I+{\mathbf{N}}_1+{\mathbf{N}}_2{\mathbf{u}}^{\left(k\right)}\odot {\mathbf{u}}^{\left(k\right)},$

(d) Solve for the displacement increment: ${\mathbf{K}}_T^{\left(k\right)}\mathrm{\Delta }\mathbf{u}={\mathbf{R}}^{\left(k\right)}$ ,

(e) Update the displacement: ${\mathbf{u}}^{\left(k+1\right)}={\mathbf{u}}^{\left(k\right)}+\mathrm{\Delta }\mathbf{u}$ ,

(f) Check for convergence: If $\Vert {\mathbf{R}}^{\left(k\right)}\Vert <ϵ$ , stop the iteration; otherwise, continue.

Figure 3.

Calculation process of critical buckling temperature and mode.

To capture the buckling path, the thermal load λΔT needs to be applied step-by-step,Divide the total load into N incremental steps ${\lambda }_1,{\lambda }_2,\dots ,{\lambda }_N$ . For each load step λi, call the Newton-Raphson iteration to solve for the equilibrium state. At each equilibrium point, calculate the minimum eigenvalue μ _min of the tangent stiffness matrix K_T.

${\mathbf{K}}_{\mathrm{T}}\mathbf{v}=\mu \mathbf{v}$

(36)

If the signs of μ _min for two adjacent load steps λ _i and λ_i+1 are opposite, it indicates that the critical point lies between them. Then, estimate the secondary buckling load λ_cr2 through linear interpolation.

${\lambda }_{\text{cr}2}={\lambda }_i+\frac{\left|{\mu }_{*\hspace{0.17em}\mathit{min},i}\right|}{\left|{\mu }_{*\hspace{0.17em}\mathit{min},i}\right|+\left|{\mu }_{*\hspace{0.17em}\mathit{min},i+1}\right|}\left({\lambda }_{i+1}-{\lambda }_i\right)$

(37)

The mode at the first critical point is the primary buckling mode. At the critical point λ_cr2, the mode v corresponding to the minimum eigenvalue of K_T is the secondary buckling mode.

The two types of boundary conditions used in this paper are as follows.

Four-side simply-supported boundary condition (SSSS):

$\text{At}\hspace{0.17em}x=0,a;\hspace{0.17em}y=0,b:u_0=v_0=w_0=0$

(38)

Four-side clamped boundary condition (CCCC):

$\text{At} x=0,a; y=0,b:u_0=v_0=w_0={\theta }_x={\theta }_y=0$

(39)

Convergence and validation

Convergence

First, a convergence study was conducted on the primary buckling temperature and the secondary buckling temperature. The material used was AS4/8552, ⁴² and its material properties are shown in Table 1. The dimensions of the laminate were 300 mm × 200 mm, and the ply-stacking sequence was [35/-35/10/-10]_s. The boundary conditions were SSSS and CCCC. A cubic NURBS was employed to construct the mesh, as shown in Figure 4. The curves of the primary critical buckling temperature and the secondary critical buckling temperature varying with the mesh are presented in Table 2. The data indicate that as the number of mesh elements increases, the results gradually converge. When the number of mesh elements is greater than 9 × 9, the results converge.

Table 1.

Material properties.

E1 (GPa)	E2 (GPa)	G12 (MPa)	ν 12	α1 (μm/m°C)	α2 (μm/m°C)	tply (mm)
135	9.68	5.6	0.3	0.28	28	0.181

Figure 4.

Cubic IGA meshes (9 × 9) and control net of the composite laminate.

Table 2.

Convergence study of λ_cr1 and λ_cr2 under SSSS and CCCC boundary conditions.

Number of elements	SSSS		CCCC
	λ cr1	λ cr2	λ cr1	λ cr2
3 × 3	24.07	74.25	34.44	100.66
5 × 5	24	69.82	32.35	87.49
7 × 7	23.98	69.17	32.34	85.19
9 × 9	23.98	68.88	31.54	84.36
11 × 11	23.96	68.84	31.5	83.7

Validation

In Ref. 13, Cetkovic employed the finite element method based on the layerwise displacement model to study the thermal buckling of composite laminated plates. He presented the variation of the first initial thermal buckling load with the aspect ratio under the SSSS and CCCC boundary conditions. The comparison between the results of this paper and those in Ref. 13 is shown in Figure 5. The results clearly indicate that there is excellent consistency between the two, which demonstrates that the method proposed in this paper has high accuracy in calculating the initial thermal buckling of composite laminated plates.

Figure 5.

Comparison of the variation of critical buckling temperature with aspect ratio under different boundary conditions.

In Ref. 59, Meyers et al. utilized the variational method and the Rayleigh-Ritz formula to study the thermal post-buckling of composite laminates. The comparison of equilibrium paths between this study and theirs is presented in Figure 6. It is evident that prior to the onset of secondary buckling, there is excellent agreement between the results of this paper and those of Meyers et al. However, since Meyers et al. failed to take into account the possibility of secondary buckling, they overestimated the post-buckling load-bearing capacity of the laminates. This overestimation represents a potential hazard in structural design. Consequently, it is highly imperative to consider secondary buckling when dealing with thin plates.

Figure 6.

Comparison of equilibrium paths with Reference 59.

In Ref. 42, Javier et al. employed experimental and finite element methods to investigate the secondary thermal buckling of composite laminates. The research indicated that after the initial buckling, when the laminate was continuously heated and the temperature reached a certain value, secondary buckling occurred in the laminate. Its load-bearing characteristics changed significantly, and the mode also changed notably. The comparison between the calculated equilibrium path in this paper and the experimental and simulation results in Ref. 42 is shown in Figure 7, the displacement is measured at the center of the plate, which corresponds to point P₀ in Ref. 42. It turns out that the calculation results of this paper are in good agreement with those in the reference. Moreover, compared with the results calculated based on the first-order shear deformation plate theory in Ref. 42, the calculation results of this paper are closer to the experimental results, which demonstrates that the method in this paper has higher accuracy. The comparison between the initial buckling mode and the secondary buckling mode is shown in Figure 8. It can be observed that there is an obvious mode transition in the secondary buckling mode compared with the initial buckling mode, and the number of longitudinal half-waves changes. Additionally, compared with the buckling modes calculated in Ref. 42, the buckling modes calculated in this paper are more consistent with the experimental results.

Figure 7.

Comparison of equilibrium paths.

Figure 8.

Comparison of buckling modes (a) Primary buckling mode-Experiment ⁴² (b) Primary buckling mode-FEM ⁴² (c) Primary buckling mode-Present study (d) Secondary buckling mode-Experiment ⁴² (e) Secondary buckling mode-FEM ⁴² (f) Secondary buckling mode-Present study.

Results and discussion

The critical buckling temperature ratio is defined as the quotient of the secondary critical buckling temperature divided by the primary critical buckling temperature. This ratio serves as a metric for quantifying the load-bearing capacity subsequent to the primary buckling event. It is denoted as:

$k_{\text{cr}}=\frac{{\lambda }_{\text{cr}2}}{{\lambda }_{\text{cr}1}}$

(40)

The equilibrium paths under different imperfection amplitudes for the SSSS and CCCC boundary conditions are presented in Figure 9. The results show that as the imperfection amplitude increases, the inflection points of the equilibrium paths shift downward. This implies that both the primary critical buckling temperature and the secondary critical buckling temperature decrease. However, the degrees of reduction for the two are different. The reduction of the primary critical buckling temperature is relatively significant, while the secondary critical buckling temperature shows lower sensitivity to the initial imperfections.

Figure 9.

Equilibrium paths under different imperfection amplitudes (a) SSSS boundary conditions.

The critical buckling temperature ratios under different imperfection amplitudes are shown in Figure 10. The results indicate that as the imperfection amplitude increases, the critical buckling temperature ratio also increases, yet the rate of growth slows down. When comparing the two boundary conditions, the critical buckling temperature ratio under the SSSS boundary condition is larger than that under the CCCC boundary condition.

Figure 10.

Critical buckling temperature ratios under different imperfection amplitudes.

The variation of the ratio of the secondary critical buckling temperature to the critical buckling temperature with the aspect ratio of the laminated plate under SSSS and CCCC boundary conditions is shown in Figure 11. It can be observed that as the aspect ratio increases, the secondary critical buckling temperature decreases. Moreover, the secondary critical buckling temperature under the CCCC boundary condition is higher than that under the SSSS boundary condition. Regarding the ratio of the critical buckling temperature, when a/b < 1, the ratio remains essentially unchanged. When a/b > 1, the ratio of the critical buckling temperature increases as the aspect ratio increases.

Figure 11.

Influence of aspect ratio on the secondary critical buckling temperature and critical buckling temperature ratio under different boundary conditions.

Assume that the stack sequence of the composite laminated plate is [θ₁/-θ₁/θ₂/-θ₂]_s. The effects of ply angles on the secondary critical buckling temperature and the ratio of the critical buckling temperature of the composite laminated plate under different boundary conditions are shown in Figures 12 and 13. The results suggest that under various boundary conditions, the influence laws of ply angles on the secondary critical buckling temperature and the ratio of the critical buckling temperature of the laminated plate are generally the same. When θ₁ = 0° and θ₂ = 0°, as well as when θ₁ = 90° and θ₂ = 90°, the secondary critical buckling temperature of the laminated plate reaches its minimum. When θ₁ = 0° and θ₂ = 90°, the ratio of the critical buckling temperature is the smallest. The secondary critical buckling temperature of the laminated plate under the CCCC condition is higher than that under the SSSS boundary condition, but the ratio of the critical buckling temperature is smaller than that under the simply-supported boundary condition.

Figure 12.

Effects of ply angles on the secondary buckling of composite laminated plates under SSSS boundary conditions. (a) Secondary critical buckling temperature; (b) Critical buckling temperature ratio.

Figure 13.

Effects of ply angles on the secondary buckling of composite laminated plates under CCCC boundary conditions. (a) Secondary critical buckling temperature; (b) Critical buckling temperature ratio.

The influence of the ratio of the elastic moduli in the two principal material directions of the laminated plate on the secondary critical buckling temperature and the critical buckling ratio under SSSS and CCCC boundary conditions is shown in Figure 14. It can be observed that as the modulus ratio increases, the secondary critical buckling temperature rises, while the critical buckling ratio decreases. This indicates that the elastic modulus ratio has a greater impact on the primary critical buckling temperature than on the secondary critical buckling temperature. Within the scope of this study, the influence of the elastic modulus ratio on the secondary critical buckling temperature is relatively limited, with a maximum increase of no more than 10°C. However, its influence on the critical buckling ratio is quite significant, especially under the SSSS boundary condition.

Figure 14.

Effects of modulus ratio on the secondary critical buckling temperature and the critical buckling temperature ratio under different boundary conditions.

The influence of the ratio of thermal expansion coefficients in the two principal material directions of the laminated plate on the secondary critical buckling temperature and the critical buckling ratio under SSSS and CCCC boundary conditions is shown in Figure 15. The results indicate that as the thermal expansion coefficients ratio increases, both the secondary critical buckling temperature and the critical buckling temperature ratio decrease. However, the thermal expansion coefficient ratio has a more significant impact on the secondary critical buckling temperature, while its influence on the critical buckling ratio is relatively limited. In particular, under the CCCC boundary condition, it has almost no effect.

Figure 15.

Effects of thermal expansion coefficients on the secondary critical buckling temperature and the critical buckling temperature ratio under different boundary conditions.

Conclusion

This paper conducts a systematic study on the secondary buckling behavior of composite laminated plates under thermal loads based on Reddy’s HSDT and the IGA method. By introducing the von Kármán large-deformation theory and initial geometric imperfections, a mechanical model considering high-order shear effects and geometric nonlinearity is established. The NURBS basis functions are employed for spatial discretization, achieving the unity of geometrically accurate description and high-order continuity. The in-plane, bending, and initial imperfection stiffness matrices in the model are derived in detail, and the Newton-Raphson iterative method is used to solve the nonlinear equilibrium equations. Additionally, a critical load prediction method based on the minimum eigenvalue of the tangent stiffness matrix is proposed. Combined with the linear interpolation technique, the secondary buckling load and its corresponding mode are accurately captured.

Numerical examples verify the effectiveness and accuracy of this method in analyzing the thermal buckling path, initial imperfection sensitivity, and high-order mode evolution of composite laminated plates. The research results show that initial imperfections significantly reduce the critical buckling load of the structure, and the combination of IGA and HSDT can accurately characterize the complex deformation behavior of thick plates, providing a theoretical basis for the thermal stability design of composite structures in engineering.

Moreover, this paper also explores the influence of different parameters on the secondary buckling behavior of composite laminated plates. The results indicate that as the aspect ratio of the laminated plate increases, the secondary critical buckling temperature decreases; under different boundary conditions, the influence of the ply angle on the secondary critical buckling temperature and the ratio of critical buckling temperatures follows basically the same pattern; an increase in the elastic modulus ratio leads to an increase in the secondary critical buckling temperature but a decrease in the critical buckling temperature ratio; an increase in the thermal expansion coefficient ratio causes both the secondary critical buckling temperature and the critical buckling temperature ratio to decrease. These findings provide important references for the structural design and optimization of composite laminated plates in thermal environments.

In conclusion, the research methods and conclusions presented in this paper not only enrich the research content of the thermal buckling behavior of composite laminated plates but also offer valuable theoretical support for relevant engineering applications. Future research can further explore the buckling behavior of composite structures in more complex thermo-mechanical coupling environments and develop more efficient numerical analysis methods to handle such problems.