RE.PUBLIC@POLIMI Buckling Behaviour of Thin and Thick Variable-Stiffness Panels

: The possibility of designing composite panels with non-uniform stiffness properties offers a chance for achieving highly-efﬁcient conﬁgurations. This is particularly true for buckling-prone structures, whose response can be shaped through a proper distribution of the membrane and bending stiffnesses. The thermal buckling behaviour of composite panels is among the aspects that could largely beneﬁt from the adoption of a variable-stiffness design, but, in spite of that, it has rarely been addressed. The paper illustrates a semi-analytical approach for evaluating the thermal buckling response of variable-stiffness plates (VSP) by considering different boundary conditions. The formulation relies upon the method of Ritz and a variable-kinematic approach, leading to a computationally efﬁcient implementation, which is particularly useful for exploring the larger design spaces, typical of variable-stiffness conﬁgurations. Due to the possibility of choosing the underlying kinematic approach as an input of the analysis, the formulation is not restricted to thin plates, but is suitable for analyzing the response of thick plates as well. Novel results are derived, which can be useful for benchmarking purposes and for gathering insight into the mechanical behaviour of variable-stiffness plates. Furthermore, the importance of transverse shear ﬂexibility is illustrated with respect to the boundary conditions as well as the degree of steering of the ﬁbers.


Introduction
Thermal buckling phenomena are of crucial importance in the design of aerospace structures [1]. High-speed aircraft and launch vehicles' structural elements are indeed exposed to aerodynamic heating that may promote elastic instability phenomena. In some cases, such as for cryogenic launch vehicles, the source of internal forces leading to instability is due to low temperatures, and cooling-induced buckling should be considered during the design process. With these motivations, thermal buckling has been the subject of many investigations in the aerospace structural community. In this context, an important role is played by rectangular plates that, despite their simplicity, are often a good approximation, at least in a early design phases, of sub-portions of more complex structural topologies. For instance, the portion of skin between stringers and frames can be assumed as a rectangular plate subjected to ideal constraints along the boundaries. Another advantage of analyzing idealized plate-like structures lies in the possibility of gathering insight into the underlying structural mechanics in a clear and straightforward manner, obtaining design guidelines that can be extended even to more complex structures.
The tailoring opportunities offered by composite materials can be further improved by allowing the orientation of the fibers to vary locally, according to the desired load path to be achieved. This chance has been known for a long time, with the first scientific paper mentioning this topic dating back to almost 50 years ago [14]. The lack of available technologies hindered the progress of this concept until the early 1990s, when the pioneering studies of Leissa and Martin [15], Hyer [16,17] and Gürdal [18] increased interest in the concept of variable-stiffness structures. In the past several years, drastic improvements in the manufacturing techniques allowed variable-stiffness constructions to become a hot topic and suitable candidates for the aerostructures of the next generation. Theoretical investigations and numerical models were carried out for analyzing their buckling [19][20][21] and post-buckling [22][23][24] response, in most cases referring to thin plate theory or FSDT [25]. Refined theories have been less commonly employed: the free vibration response of variable-stiffness plates (VSP) is investigated in [26,27] using finite elements, while the authors presented a variable-kinematics, Ritz-based procedure for buckling and free-vibration analysis [28] of VSP based on Carrera's Unified Formulation (CUF) [29][30][31]. Tornabene and co-workers discussed the application of variable-kinematics approaches for analyzing the modal and static response of variable-stiffness laminates [32,33] and sandwich constructions [34,35].
Recently, the thermal buckling behaviour of variable-stiffness plates was investigated in [36], using a thin-plate finite element formulation, and [37], by considering FSDT. Thermally induced multi-stable configurations were assessed by Haldar et al. [38] using finite elements and the Ritz approach proposed by Dano and Hyer [39].
To the best of the author's knowledge, no previous attempts could be found in the literature to analyze the thermal buckling of VSP using refined theories. The present work aims at filling this gap by presenting a unified Ritz-based modeling framework, where the analyst has the possibility of choosing the underlying plate theory-both as equivalent single layer (ED) or layer-wise (LD)-so that both thin and thick plate configurations can be studied.
Furthermore, novel closed-form solutions are presented for the evaluation of pre-buckling internal stress distribution, which can be successfully employed for gathering further understanding into the mechanical response of VSP, as well as improving the efficiency of the buckling solution procedure.

Theoretical Framework
This work is devoted to the thermal buckling analysis of variable-stiffness composite plates, i.e., flat plates characterized by the stacking of plies with non-straight fibers.
Several laws were proposed in the literature for describing the orientation of the fibers within the panel domain, including linear and nonlinear ones [20,25,40]. Linear variation of the orientation angle is considered here due to its effectiveness in allowing the description of a variable-stiffness ply by means of two parameters only. This choice tends to limit the design space with respect to more general nonlinear descriptions, but it appears useful when trying to condense the behaviour of VSP in design charts or compact summary tables, as done in this study. In other words, the understanding of the mechanical response of VSP, their potential for elastic tailoring, as well as the main features associated with their analysis can be grasped even by restricting the focus to the case of fiber linear variation.
The plates are characterized by planar dimensions a × b, where a denotes the longitudinal direction along the x-axis and b is the transverse direction parallel to the y-axis, as illustrated in Figure 1. The total thickness is denoted with h.
According to the notation adopted here, a ply denoted with < θ 0 |θ 1 > x is characterized by a fiber variation along the x-axis, where θ 0 is the orientation, measured with respect to the x-axis, at the plate center, and θ 1 the one at the plate edge (see Figure 1a). Similarly, steering of the fibers along the direction y is denoted with < θ 0 |θ 1 > y (see Figure 1b). Starting from the fiber passing through the plate center, the other ones are obtained by translation of the reference fiber along the xor y-direction, referring to a strategy sometimes denoted as a shifted method [26,41]. By considering fiber variation along the x-direction, the orientation angle at a generic point is given as: An analogous linear law is obtained when variation takes place along the y-direction. Referring to common design rules of the aerospace industry, the stacking sequence is taken as symmetric with respect to the midplane, i.e., a layer oriented at < θ 0 |θ 1 > at the distance z from the plate midplane implies the presence of another identical layer at the distance −z. In addition, the analysis is restricted to the case of balanced laminates, meaning that a ply at < θ 0 |θ 1 > requires the presence of another ply at − < θ 0 |θ 1 > at any position in the stack.

Variational Formulation and Approximate Solution
A variable-kinematic approach is considered due to the possibility of analyzing wide ranges of plate configurations. Specifically, the underlying framework refers to the well-known Carrera's Unified Formulation (CUF) (see, for instance, [29,42,43]), offering the advantage of allowing a relatively easy implementation, which is independent from the kinematic model to be used for the structural analysis. Indeed, the theory is defined as an input of the problem, and can be chosen among the class of Equivalent Single Layer (ED) and Layer-Wise (LD) formulations.
The formulation is developed in the context of a displacement-based approach, where the three displacement components are the unknowns of the problem. The buckling condition is sought by application of the Trefftz criterion, which is written as: where δ 2 Π is the quadratic part of the Taylor expansion of the total potential energy, whose expression is given as: where Ω denotes the domain [−a/2 a/2] × [−b/2 b/2], and k is the ply-related index. Note that the first five contributions entering the right-hand side of Equation (3) should be intended as variations with respect to the equilibrium condition, whilst the last one, σ k p0 , defines the membrane stress state obtained from the pre-buckling analysis.
The components of the deformation and stress tensors reported in Equation (3), in the case of a generic 3D elasticity formulation of the problem, are split into the in-plane and normal contributions, according to the definitions: while the pre-buckling membrane stresses are: The strain components in work-conjugacy relation with those of Equation (6) are the nonlinear part of the Green-Lagrange strain tensor, which is given by: The functional of Equation (3) can be expressed as a function of the displacement vector u k = u k v k w k T after introducing the constitutive law, in this case taken as the generalized Hooke's Law in 3D for an orthotropic material. The approximate solution is sought here referring to the Ritz method. The procedure is briefly outlined in the following, and the interested reader is referred to [28,[44][45][46][47][48] for more details.
According to the CUF formalism, the displacement components are expanded as the product of thickness functions F τ and generalized displacement components u k τ as: where ξ and η are non-dimensional in-plane coordinates, defined as ξ = Following the classical Ritz scheme, the generalized displacement components u k τ are represented as the linear superposition of an arbitrary number of trial functions: where summation is implied with respect to the index i, whilst r defines the generic displacement component. The coefficients c k rτi are the unknown Ritz amplitudes, and the trial function N rτi reads: where p i are Chebyshev polynomials, while f and g are polynomial boundary functions (for a definition, see [49]).
After substituting Equation (9) into the functional given by Equation (3), and applying the criterion given by Equation (2), the set of discrete equations governing the buckling response is obtained as: where K and G are the linear and the geometric stiffness matrices, respectively. The scalar λ is the buckling multiplier, and is obtained as the smallest magnitude eigenvalue of the eigenproblem established by Equation (11). In the context of thermal buckling analysis, it is generally of interest to identify the smallest positive and negative eigenvalues, in order to obtain the critical temperatures associated with heating and cooling. In some cases, all of the eigenvalues have the same sign, thus instability is possible only for cooling or heating.

Pre-Buckling Solutions
The geometric stiffness matrix G depends on the membrane pre-buckling stresses, whose evaluation represents the first step of the overall buckling procedure. To this aim, closed-form solutions are derived for establishing the state of pre-buckling internal forces-to the best of the author's knowledge, they were not previously available in the literature-for four sets of boundary conditions. The conditions under investigation are summarized in Figures 2 and 3. According to the nomenclature adopted hereinafter, the term Case-Txi specifies a configuration characterized by fiber variation along the x-axis, and i is equal to: 1 if all the four edges are restrained against the expansion/contraction along the direction normal to the plate edge; 2 if only the longitudinal edges are prevented from the movement along the direction normal to the plate edges. Similarly, Case-Tyi conditions are relative to two possible cases associated with fiber variation along y.
The approach for deriving the pre-buckling, thermally-induced stresses refers to classical membrane theory of laminated plates. Under the assumptions introduced here in terms of stacking sequence and boundary conditions, the solution are exact for a generic 3D laminate. Similar results, but restricted to the purely elastic case, can be found in [50].

Thermoelastic Strains and Constitutive Relation
Following the classical approach adopted in thermoelastic analysis of laminates [51], the vector of membrane strains can be organized as: where the superscript tot is the total deformation, which is represented as the sum of a contribution due to thermal and mechanical strains, th and m . With respect to Equation (4), the subscript p is now omitted.
Assuming linear thermoelastic behaviour, the thermally-induced strains are expressed as: having indicated with ∆T a temperature variation, taken positive for heating, with respect to the stress-free reference condition, andα the vector collecting the laminate coefficients of thermal expansions, which are obtained as: where A is the well-known membrane stiffness matrix, and its inverse a = A −1 is the membrane compliance. The unit thermal stress resultantN th of Equation (14) can obtained as [51]: Recalling now the constitutive equation of a symmetrically layered membrane: and inverting Equation (16) and using Equation (15), is obtained, which provides an expression for the total strains in terms of membrane forces and temperature variation.

Semi-Inverse Approach
The pre-buckling problem is solved by means of a semi-inverse approach, where the initial assumption regards the shear membrane force per unit length N xy , which is taken identically null inside the plate domain and along its boundaries, i.e., N xy = 0. The closed-form solutions for N xx and N yy are then derived by enforcing the fulfillment of the equilibrium equations and the kinematic and equilibrium conditions at the boundaries. The final step of the approach consists of verifying the compatibility of the solution, providing a proof that the resulting solution is indeed the exact solution of the linear elastic pre-buckling problem.
The in-plane pre-buckling equations, which are uncoupled from the out-of-plane one due to the assumption of symmetric laminate, are given as: Starting from the tentative solution characterized by N xy = 0, it is seen, from Equation (18), that the equilibrium equations reduce to: meaning that the direct membrane forces are a function of one single coordinate. Depending on the essential conditions at the boundaries, which are discussed next, and starting from the equilibrium requirement given by Equation (19), it is possible to determine the pre-buckling solution. The resulting internal force distribution is verified ex-post to be respectful of the compatibility condition, which is given by: xx,yy + yy,xx = γ xy,xy .
Note that Equation (20) can be written in terms of force resultants N ik after expressing the strain components xx , yy and γ xy as function of the laminate membrane compliance coefficients.
It is noted that all the solutions derived next do satisfy equilibrium and compatibility conditions, although the verifications of Equation (20) are omitted here for the sake of brevity.

Case-Tx1
The first set of conditions is illustrated in Figure 2a, and is relative to a VSP with fibers varying along the direction x. The in-plane conditions are such that normal displacements at the boundaries are zero. Recalling Equations (17) and (19), the axial shortening along the x-direction is then: where a ik are the coefficients of the compliance matrix a. The expression of Equation (21) shows that N xx = N xx is constant across the panel domain.
The shortening along the transverse direction is: Substituting Equation (22) into Equation (21) and rearranging the expressions leads to: where:

Case-Tx2
The case illustrated in Figure 2b refers to a plate that is constrained along its transverse edges to prevent the normal displacement along the x-direction. The transverse sides are free to expand or contract, but are constrained to remain straight. The first essential conditions regarding the transverse edges are given by Equation (21), and can be written as: which is clearly verified by ∀y. The transverse displacement is unknown and, in general, different from zero. Referring to Equation (22), its expression is derived as: The second boundary condition regards the vanishing of the force resultant directed along the y-axis, i.e., Substituting Equation (26) into Equation (27) allows for expressing the transverse displacement as a function of the membrane force N xx and the temperature variation ∆T as: where: Substituting back Equation (28) into Equations (26) and (27) and solving for N xx and N yy leads to: where:

Case-Ty1
The case relative to the plate with fiber variation along the y-direction, and subjected to the boundary conditions of Figure 3a is analogous to the one discussed in Case-Tx1. The solution is readily available by modifying Equation (23) to account for the π/2 rotation between the two cases (in other words, this is obtained by flipping the indexes x and y, and 1 and 2): where:

Case-Ty2
Pre-buckling boundary conditions relative to Case-Ty2 are depicted in Figure 3b, and refer to fiber steering along y, with normal displacement prevented at the two transverse edges, and stress free conditions along the longitudinal ones.
The first kinematic condition is imposed as: while the second condition regards the membrane forces at the two edges: which implies that N yy = 0 over the entire plate domain. It follows that the pre-buckling solution is: where: It is straightforward to verify, through the evaluation of the transverse displacement ∆v, that the condition of Equation (35), along with the fiber variation along the y-direction, implies that the longitudinal edges remain straight during the pre-buckling deformation process.

Comparison with Literature Results
A first part of the section is devoted to the validation of the present Ritz-based procedure with results available in the literature and finite element calculations.
The first comparison is presented against the 3D exact solutions derived by Noor and Burton [12] for a set of angle-ply, straight fiber configurations. Despite the focus of this preliminary example not being expressly directed towards the case of variable-stiffness plates, it is of particular interest because the reference solutions are derived by means of an exact three-dimensional approach.
The plates are square and different ratios b/h are considered for illustrating the role played by transverse shear deformability. Pre-buckling constraints Case-Tx2 are considered, i.e., the normal motion of the four edges is prevented, whilst simply supported boundary conditions are considered along the four edges for solving the buckling eigenvalue problem.
An orthotropic material is considered, having the following non-dimensional thermoelastic properties: The lay-up is characterized by the stacking of 10 plies oriented at [±θ] 5 . It should be noted that, for the problem at hand, the exact elasticity solution of the pre-buckling problem (see also [12]) is characterized by a membrane state of internal stresses, where the only not null components of the stress tensor are those reported in Equation (6). It follows that no approximations are introduced in the pre-buckling analysis if Equation (30) is used.
The results are summarized in Table 1 for different values of b/h, ply angles and equivalent single layer kinematic theories of order 2,3 and 4. Based on a preliminary convergence analysis, Ritz results are computed using 12 functions along both directions.
As seen from Table 1, close agreement can be achieved between the 3D results reported in [12] and those calculated with the present Ritz formulation. As expected, thicker plates demand increased orders in the kinematic model to guarantee adequate accuracy of the predictions. When b/h is equal to 100, no substantial improvement is observed by refining the theory from the second to fourth order. On the contrary, the critical temperature for thick plates with b/h of 20/3 can be accurately estimated using the third-or fourth-order theories ED3 and ED4. In general, it can be noted that the discrepancy with the reference solution is higher for increasing values of θ. This behaviour is due to the skewness of the in-plane displacement pattern when the plies are rotated in the off-axis direction. In these cases, the convergence of the solution is slower, and a proper description of the buckling pattern requires a larger number of trial functions [45].
A second test case is taken from [37], where the thermal buckling analysis is presented for variable-stiffness plates subjected to restrained pre-buckling thermal expansion and simply supported boundary conditions. The material properties are those of an orthotropic material, whose non-dimensional thermoelastic parameters are given as: Three stacking sequences are considered, with plies oriented at 4 , where the angles are measured with respect to the x-axis, and the variation is linear with x.
The analysis of [37] is conducted based on finite element models, using FSDT as underlying kinematic theory with a shear factor equal to 5/6.
To further validate the results, finite element computations are performed with Abaqus finite element simulations. To this aim, a mesh size of 50 × 50 S4R shell elements was found adequate to obtain converged results. Note that the mesh is constructed such that each element is characterized by a constant orientation, thus the fiber steering is represented in a stepwise manner. Given the linear variation of the fiber angle, each row of elements at x = const share the same orientation, with a discrete increase of the angle when moving along the x-coordinate.
The results are summarized in Table 2, where the present Ritz formulation is exploited by considering three different kinematic theories. In addition, the convergence of the method is assessed by considering an increasing number of trial functions R = S. The convergence of the Ritz solutions is clearly seen from Table 2, where the critical temperatures get smaller and and smaller as the number of trial function increases. It is observed that the number of integration points for evaluating the Ritz in-plane integrals is taken as P + 5. Due to the numerical integration process, convergence of the temperatures from above cannot be guaranteed a priori, although it is generally experienced [28]. For the case at hand, a few degrees of freedom are necessary, in most cases, to reach convergence up to the first two digits. Overall, the convergence is faster for the straight fiber configuration, whose orientation angle at 15 is responsible for a slight amount of buckling mode skewness. The convergence is slower when variable-stiffness panels are of concern, with increasing need to consider more trial functions as the steering of the fibers is stronger.
As seen by the results in Table 2, a substantial agreement is obtained with the predictions of [37], although in some cases the discrepancies can be non-negligible. For instance, the Ritz results relative to the [± < 15|75 > x ] 4 lay-up are higher with respect to those of [37], irrespective of the kinematic theory. For this reason, Abaqus analyses were performed, demonstrating, on the contrary, very close agreement. Note that a local buckling mode, characterized by several short half waves in proximity of the boundaries, is predicted as the first instability by Ritz and Abaqus analyses for the plate with lay-up [± < 15|75 > x ] 4 and b/h = 10. Local instabilities are not of concern in this study, so the values in the table refer to the first global mode. For the configurations with b/h equal to 40, Abaqus and Ritz results are almost identical. Clearly, no substantial advantages are associated with the adoption of higher-order theories, and an ED2 model is generally sufficient. The matching with finite element results is clear even for thicker configurations. In these cases, the adoption of ED4 or LD2 theories can be beneficial to obtain accurate predictions, whilst a second-order theory ED2 is generally responsible for large errors. It is noted that LD2 results generally lead to smaller critical temperatures with respect to Abaqus analyses based on shell elements, as one may expect. In a few cases, this behaviour is not respected due an over-correction associated with the evaluation of shear corrections factor performed by Abaqus for non-symmetric lay-ups. Similar issues were observed in [46]. It is useful to note that the higher-order theories available in the variable-kinematic formulation do not require the definition of a shear factor, and the effect of transverse shear flexibility is inherently accounted for by the kinematic models themselves.
An extensive set of results for the thermal buckling of thin variable-stiffness panels is provided by Duran and co-workers [36]. The authors analyzed several composite materials, focusing on the sensitivity of the optimal configuration with respect to the material properties. The analyses of [36] are based on finite element analyses using four-node, Kirchhoff finite elements, so they are restricted to thin-plate configurations.
A set of results is presented in Table 3, where square plates are analyzed, with non-dimensional ratio b/h=148. The pre-buckling boundary conditions are those denoted here as Case-Tx1, while the flexural ones refer to the simply supported assumptions. Different materials and lay-ups are considered. For the sake of conciseness, the material properties are not reported here but can be found in [36]. The lay-ups are characterized by a linear variation of the fiber angle along the direction x, with stacking sequence given by [± < θ 0 |θ 1 > x ] s . It is noted that the configurations of Table 3 are those maximizing the critical temperatures according to the optimizations performed in [36] and, for this reason, the angles θ 0 and θ 1 are non-integer values. The results are reported by considering the equivalent single layer theories ED3 and ED4; furthermore, the second-order layer-wise theory LD2 is adopted for furnishing lower bound reference results.
As observed, close agreement is noted between the temperatures obtained by Duran et al. [36] and the present formulation. In most cases, the differences are of a few percent points, thus confirming the correct implementation of the Ritz code. The same cannot be claimed for the Kevlar/Epoxy and Carbon/Epoxy configurations, where the discrepancies are very high. To clarify this aspect, Abaqus results were performed for all the configurations: as seen from Table 3, the differences between Ritz and Abaqus results are very small for all the cases. It is then concluded that even the evaluation of the Kevlar/Epoxy and Carbon/Epoxy plates is correct. As a further verification, the comparison against Abaqus analyses is reported in Figure 4 for the Kevlar/Epoxy plate in terms of pre-buckling stress resultants and buckled shape. The pre-buckling forces are evaluated for a unitary temperature increase, and are characterized by a uniform N xx distribution, and a nonlinear one for N yy , according to Equation (23). The quality of the semi-analytical Ritz predictions can be noticed, also with reference to the skew buckled pattern, which is correctly predicted using a few degrees of freedom.   Following the study of [36], the Graphite/Epoxy configuration is analyzed by evaluating the critical temperatures for all the combinations of θ 0 and θ 1 . between -90 and 90. Note that the current analyses do not consider any technological requirement, nor the effects of strong fiber steering on the constitutive law.
The results are calculated using ED3 theory with 12 × 12 functions, and are reported in the plot of Figure 5.
The contour plot is very similar to the one reported in [36], to which the interested reader is referred. The symmetry of the plot can be noted, due to the plate configuration, which is such that T cr (θ 0 , θ 1 ) = T cr (−θ 0 , −θ 1 ). By considering an angle step of 2.5, the configuration guaranteeing the maximum is characterized by the lay-up [± < 57.5|40.0 > x ] s . The orientation of the angle at the plate center is close to the optimal one reported in Table 3, while a slight difference is noted with respect to the fiber orientation at the plate edge.
Overall, this kind of analysis illustrates the potentialities of a Ritz-based procedure. Indeed, the possibility of exploring the design spaces offered by variable-stiffness configurations is subjected to the availability of computationally efficient analysis tools. For instance, the total number of analyses for generating the plot of Figure 5 is, in this case, beyond 5000. The Ritz approach, due to the few degrees of freedom required to achieve convergence, is particularly suited for this, and the time for the analysis is on the order of few seconds.

Pre-Buckling and Buckling Response of Carbon/Epoxy VSP
Having demonstrated the capabilities of the Ritz formulation, further results are discussed in this section. The present study focuses on carbon/epoxy composites, which are commonly employed in advanced constructions such as in the aerospace field. The material properties are those relative to the Graphite/Epoxy material of Table 3, and taken from [36]. For convenience, they are reported here below: It can be noted that the orthotropy ratio is approximately 20; the thermal coefficient along the fiber direction is slightly negative, meaning that a shortening is experienced by the material when subjected to heating, while the value is positive along the transverse direction, i.e., the matrix expands when heated. The distribution of thermal and elastic properties influences the pre-buckling stress distribution, and, in general, can be tailored by allowing the fibers to change their orientation in the panel domain. The pre-buckling distribution of internal membrane forces, in turn, affects the buckling response of the plate. Furthermore, the overall response depends upon the boundary conditions applied to the plate, with in-plane ones having the twofold effect of influencing the pre-buckling response and buckling one. Note that, for symmetric thin laminates, in-plane boundary conditions are relevant only in terms of pre-buckling response, as the out-of-plane buckling equation is uncoupled from the in-plane ones. For thick plates, whenever the assumption of in-plane stress constitutive equation is not adequate to model the plate elastic response, an inherent coupling between in-plane and out-of-plane behaviour is recovered, and the buckling behaviour is affected by the in-plane boundary conditions as well.
A set of panels is studied below here by considering the orientation of the fibers to vary along the y-axis, and the in-plane boundary conditions Case-Ty1 and Case-Ty2. The panels are square, with dimensions equal to 150 mm, and are obtained by the stacking of four plies oriented at [± < θ 0 |θ 1 > y ] s , where the angles are measured with respect to the x axis. To illustrate how the panel response is affected by the fiber path, the angles θ 1 are varied between 0 and 90, while two distinct angles θ 0 are considered, and taken equal to 15 and 75. The first case corresponds to a plate with a relatively high stiffness along the x direction, at least in the middle strips of the plate close to y=0; in the second case the central region is much weaker along the x-direction due to the almost transverse orientation of the fibers.
Firstly, the results are assessed in terms of pre-buckling behaviour; then, the buckling response is illustrated, reporting also the comparison with finite element calculations.

Pre-Buckling Analysis
Results relative to the boundary conditions Case-Ty1 are presented in Figure 6. It is recalled that the approach developed in the previous section is exact. The comparison against finite element calculations is not presented for the sake of conciseness. However, the correctness of the derivation and its implementation has been verified through an extensive set of comparisons against finite element results. The membrane forces per unit length are calculated for a unitary temperature increase according to Equation (32), and are reported by illustrating the resultant N xx in the left portion of the plot, and the transverse resultant N yy in the right part. Positive values denote traction, while negative ones are associated with compression. Note that the curves are reported for half of the plate's width. The remaining part can be easily recovered by exploiting the symmetry of the internal distribution with respect to the vertical axis.
Referring to Figure 6a, one can note that low values of θ 1 are associated with fiber paths running closer to the longitudinal direction. Given the high value of the ratio |α 22 /α 11 |, the thermal expansion mostly relates to the transverse direction of the laminate. Due to the prescribed null displacement along the longitudinal edges, a larger force N y is introduced with respect to those cases with higher angles θ 1 . At the same time, the compressive force N x tends to increase its value for increasing orientations of the fiber angle θ 1 at the panel's edge. Looking, for instance, at the case with θ 1 equal to 90, it is possible to notice the higher compressive values at the outer edge, where the thermoelastic response is dictated by the response along the matrix direction. The compressive force becomes smaller when moving towards the center of the panel, i.e., when the fibers are progressively rotated along the axial direction x.
The second case of Figure 6b is relative to a panel with fibers oriented at ±75 at the panel's center. In this case, the steering to small values of θ 1 (see, e.g., the case for θ 1 = 0) determines a matrix-dominated response along the y-direction. The thermal dilatation is thus reacted with higher values of N y . As it concerns the longitudinal direction, one can note that compressive forces are milder on the outer regions, as far as the fibers are oriented with an angle θ 1 , which is smaller than θ 0 : in these cases, the steering is such that the matrix contribution to thermal dilatation increases when moving towards the outer regions. An opposite response is observed for θ 1 equal to 90, where the thermal dilation coefficient at the outer part is equal to that of the matrix, and larger with respect to the center, and similarly the force N x .
A second set of results is available in Figure 7, where the in-plane boundary conditions Case-Ty2 are considered. The displacements along the normal direction are prevented at the transverse edges, whilst the longitudinal ones are free to expand or contract. This set of constraints is representative of a design situation where two sides of the plate are restrained with stringers characterized by high in-plane bending stiffness, such that the expansion is prevented. On the contrary, the in-plane bending stiffness of the two other parallel sides is much smaller, and the motion along the normal direction is not restrained. From the results of Equation (36), it is seen that no membrane forces develop along the transverse direction as a response to a uniform heating or cooling. The distribution of internal forces is then purely uniaxial. The behaviour of plates characterized by θ 0 equal to 15 is summarized in Figure 7a. One can note that steering the angle θ 1 from 0 to 90 has the effect of increasing the laminate thermal coefficient, thus increasing the amount of compression introduced at the outer edges. For the case at hand, the smallest values of θ 1 (0, 22.5) are characterized by a state of internal traction: thermal buckling is thus possible only if the plate is cooled. The same conclusion holds for θ 1 equal to 45 and 67.5, as the compression region is relatively small. When θ 1 is equal to 90, an internal state of compression is promoted, thus making the onset of instability phenomena possible due to heating.
The results of Figure 7b are relative to configurations where the fiber angle θ 0 is equal to 75. This set of configurations is particularly interesting, as a mechanical buckling-driven design would suggest configurations where the mid angle is approximately perpendicular to the load direction in order to allow load re-distribution towards the edges. In this case, fiber angles θ 1 equal to 90, 67.5 and 45 are associated with compressive pre-buckling force per unit length N x . In other words, the matrix contribution to thermal expansion is dominant along the direction y, resulting in a compressive force in response to the prevented displacement along y at the two transverse edges. For relatively high values of fiber steering, i.e., when the external fibers are oriented at 0 or 22.5, the small thermal expansion associated with the fiber direction is such that the external parts of the plate experience a traction force resultant, whilst the inner part still is subjected to the matrix-dominated compression forces. Due to the contemporary presence of compressive and tensile internal forces, this latter configuration can undergo instability phenomena due to thermal heating.

Buckling Analysis
Non-dimensional buckling temperatures are reported in Table 4 for the plates discussed in the previous section. In-plane boundary conditions Case-Ty1 and Case-Ty2 are considered, while simply supported assumptions are introduced for analyzing the buckling problem. The results are presented by considering a third-order equivalent single layer theory, ED3, while 14 trial functions are considered along both the orthogonal directions. Table 4. Non-dimensional thermal buckling loads Tα 22 × 10 3 for [± < θ 0 |θ 1 > y ] s laminates subjected to SSSS boundary conditions. Results obtained using ED3 and 14 × 14 trial functions.

Case-Ty1
Case-Ty2 The results illustrate the dependence of the critical temperature with respect to the fiber path, boundary conditions and geometry of the plate. As expected, smaller values of b/h, i.e., thicker plates, lead to higher critical temperatures. In particular, it is noted that the pre-buckling internal forces do not depend on the panel relative thickness b/h, but they are a function of the thermoelastic properties only, i.e., the free thermal strains are not affected by the geometry of the plate. As a consequence, plates with different values of b/h, and heated with a unitary temperature increase, are subjected to the same internal resultants. On the contrary, the buckling condition depends upon the square of the ratio h/b, meaning that thicker plates require higher levels of internal membrane forces to buckle, and thus higher temperature increase. Based on these remarks, one can observe that the ratio between the critical temperatures associated with values of b/h equal to 100 and 50 are slightly smaller than four due to effects of transverse shear flexibility. Similarly, the ratio of the temperatures of plates characterized by b/h equal to 50 and 20 is slightly smaller than 6.25.
Looking at the role played by in-plane boundary conditions, it can be noted that, for a given non-dimensional parameter b/h, all of the critical temperatures fall in a relatively small range of value when the plate is constrained with Case-Ty1 conditions. In particular, the maximum ratio between the highest and the lowest temperatures is equal to 2, approximately. This is due to the inherent biaxiality of the internal pre-buckling forces. As noted, from the plots of Figure 6, the fiber steering has, in general, the effect of mitigating the internal compression along one direction, while increasing that on the mutually orthogonal direction, and vice versa.
A different response is observed for the in-plane conditions identified by Case-Ty2. In this case, the uniaxiality of the internal forces leads to a more pronounced dependence on the fiber path. Some configurations do not buckle under thermal heating, while they are prone to buckling as far as the fiber is steered along the transverse direction at panel edges.
For clarity, a design chart is presented in Figure 8, where the critical temperatures are reported for Case-Ty2 conditions and different stacking sequences. The results for straight fiber configurations display, with close agreement, the behaviour illustrated by Nemeth [4]: heating-induced instability is possible in the range between θ equal to 44 and 90 degrees, while cooling is the only mechanism promoting buckling in the range between 0 and 43 degrees.
The curves of Figure 8 highlight the effect of steering the fibers. One can note that the region of instability due to heating tends to shift on the upper-right part of the plot as θ 0 is reduced. Indeed, resistance against thermal heating is provided by longitudinally directed fibers. On the contrary, the effect of progressively rotating the center plate fibers along the longitudinal direction is that of shifting the cooling stability region to the bottom left of the plot: the matrix-dominated behaviour, at least in the central part of the plate, has the effect of improving the resistance to cooling. It is interesting to highlight that, for some range of values, buckling is possible both for heating and cooling, for temperatures that may reasonably be retained in the field of material linear response. For instance, when θ 0 is equal to 75 and θ 1 is 0, positive and negative critical non-dimensional temperatures are possible with absolute values of 24.74 and 50.68, respectively. The co-existence of positive and negative buckling multipliers is not possible for straight-fiber configurations, where the uniform pre-buckling condition is strictly compressive or tensile. The outlined response is thus peculiar of variable-stiffness plates, where the variability of the internal pre-buckling forces allows for the presence of non-uniform membrane force resultants with compressive and tensile regions. With similar motivations, one can observe that thermal buckling due to heating is always possible for the configurations with θ 0 equal to 75, irrespective of the values of the fiber orientation at the edge θ 1 .
Additional results are presented in Figure 9, where the comparison is presented in terms of pre-buckling stress resultants and buckled shape for three variable-stiffness configurations with lay-up [± < θ 0 |90 > y ] s , and one constant stiffness panel with θ 0 = 90. The internal stress distribution is evaluated referring to Equation (36) by considering a unitary temperature increase. As seen from Figure 9a-d, the internal stress non-uniformity is progressively reduced as the angle θ 0 is increased up to 90. Note that, for clarity, the same color scale is adopted in all of the plots. With the exception of the case θ 0 = 15, which is the buckling configuration associated with the plate cooling, all the other buckled surfaces are characterized by two half-waves, with increased skewness for smaller values of θ 0 due to increased bending-twisting coupling. Despite the similarity of the buckling modes, it is noted that the critical temperatures, available from Figure 8, are drastically affected by the steering of the fibers. The case of clamped conditions is summarized in Figure 10 to quantify the effect of flexural boundary conditions. As expected, the effect of adding constraints is that of raising the values of critical temperatures. It can be noted that the same trends observed for the simply supported case are essentially preserved, the main difference regarding a translation of the stability regions towards higher absolute values of the temperature.
The effects of normal and transverse shear deformability are investigated in Table 5 by presenting the critical temperatures for different values of b/h. All the results are obtained referring to a layer-wise, second-order LD2 theory, which guarantees the possibility of achieving quasi-3D predictions. Aiming at highlighting the role played by higher-order deformability effects, the results are now presented according to the following non-dimensional form [4]: In other words, the non-dimensional temperature is expressed in terms of the buckling coefficient T cr , which, for thin plate theory, eliminates the dependence on the geometric parameter b/h. In the general case of higher-order theories and for fixed plate configuration, any variation ofT cr is only ascribable to normal and transverse shear flexibility effects.  The results provide a clear insight into the importance of properly accounting for shear deformation effects, when thick and moderately thick plates of concern. For the configurations at hand, the buckling coefficients associated with b/h equal to 10 are between 0.50 and 0.78 times the values registered for the thin configurations with b/h equal to 1000.
It is interesting to note that the influence of the fiber orientation on transverse shear effects is relatively weak when the plate is constrained with Case-Ty1 conditions. This behaviour is ascribable to the inherent biaxiality of the internal membrane forces, which tends to reduce transverse shear deformation effects.
When the plate is constrained with Case-Ty2 conditions, the internal forces are uniaxial, and the effects of transverse shear deformability are more pronounced. Referring to Table 5, it can be noted that the relevance of transverse shear effects is higher for increasing intensities of fiber steering. The effect is at a maximum for the lay-up [± < 75|10 > y ] s , where the buckling coefficient associated with b/h=10 is almost half of the one predicted for the thin plate counterpart. In those cases, the need for adopting proper kinematic theories is evident. At the same time, even the analysis of moderately thick plates in the range of b/h between 50 and 20 is susceptible to noticeable approximation if thin plate theory is adopted, especially for variable-stiffness configurations with aggressive steering.

Conclusions
The work illustrated the analysis of the thermal buckling response of VSP using a variable-kinematics Ritz approach. The pre-buckling analysis is performed analytically and, to this aim, a set of novel closed-form solutions was derived for four different sets of in-plane boundary conditions. The availability of those solutions is particularly useful for obtaining a clear insight into the pre-buckling stress redistribution. Furthermore, it facilitates the overall buckling procedure as no linear systems need to be solved to determine the pre-buckling internal stress distributions. Due to the effectiveness of the approach, case studies can be easily conducted and design charts plotted. The results illustrate that highly accurate predictions can be achieved for thin and thick plates with uniform and variable-stiffness. The adoption of high-order theories is particularly relevant when the normal edge displacement is prevented at two parallel edges, and the two others are free to expand. Furthermore, high-steering of fibers has the effect of increasing the role played by transverse shear deformability, thus increasing even more the need for high-order theories.