Local and global nonlinear dynamics of thermomechanically coupled composite plates in passive thermal regime

Uniﬁed 2D continuum formulation of the nonlinear dynamic problem for a von K´arm´an shear indeformable symmetric cross-ply composite plate in a thermomechanical environment is presented, along with the ensuing reduction procedure ending up to a three-mode discretized model with unknown transverse displacement and membrane/bending temperatures. Systematic numerical analyses in the case of thermal dynamics passively entrained by the solely active mechanical excitations allow to unveil the main features of the nonlinear response, while highlighting fundamental aspects associated with the thermomechanical coupling. Local and global dynamics of a single-layer orthotropic plate are investigated under varying in-plane/transverse excitations or thermal property of the material. Comparison with the response provided by partially coupled models and the uncoupled mechanical oscillator enables to identify situations in which thermomechanical coupling affects the nonlinear response even in the solely passive thermal setting and to frame the relevant effects within known literature results.

However, in nonlinear dynamics, low-order models preserving the main features of the underlying continuum formulation may allow easier analyses and in-depth understanding of the basic, yet possibly involved, effects of coupling on the finite amplitude vibrations of geometrically nonlinear laminates. Indeed, they get rid, at least in the initial stage of the thermomechanical investigation, of the complicatedness often ensuing in the interpretation of nonlinear phenomena from the use of much richer models (e.g., finite elements), also possibly implemented within an effective unified framework [7,19].
In such a perspective, and in view of overcoming the mere juxtaposition of a non-systematic treatment of thermoelastic coupling to plate mechanical theories [23,25], a minimal reduced order model (ROM) of thermomechanically coupled plates has recently been developed in the framework of a unified 2D continuum formulation of the problem [28]. This has been done by introducing and developing mechanical and thermal variables, along with the relevant equations, in parallel, based on the classical Tonti diagram [29,30] for physical theories. Specifically, Tonti's approach has been effectively used to implement a unified 2D formulation of Classical von Kármán laminated plates with Thermomechanical Coupling (CTC) and assumed linear temperature variation along the thickness [28], though generally allowing to obtain and compare a multitude of models resulting from different, yet controllable and consistent, assumptions about the mechanical and thermal configurations/quantities [31]. In [28], the CTC model has been used to exemplarily describe the nonlinear vibrations of symmetric cross-ply laminates under various modeling simplifications, focusing on the phenomenon of thermoelastic damping of mechanical vibrations entailed by temperature variations Yet, though being low-dimensional, the CTC ROM is quite rich in terms of material thermal properties and variety of mechanical and thermal excitations to be possibly addressed, thus allowing us to systematically analyze the features of nonlinear dynamic response in the absence of internal resonances and to unveil possibly meaningful phenomena of multifield interaction. In this respect, it is important to distinguish between material and excitation conditions under which the coupling terms in the mechanical and thermal equations entail a solely dragged, i.e. passive, presence of thermal phenomena within the system overall response, and those entailing an actually active role of coupled thermal dynamics. All of this certainly helping to identify situations of possible inadequacy of an uncoupled, though thermally affected, mechanical ROM to grasp the system actual response, and shedding light on the need of comprehensive analyses allowing to understand the rich dynamical scenarios of actually coupled systems.
Being the amount of situations of theoretical or technical interest quite large, and the underlying analyses computationally demanding, in the present paper attention is paid to the system nonlinear response in the sole condition of passive thermal dynamics, leaving the analysis of the effects entailed by the presence of active thermal sources to a companion study. In particular, the role played by the thermomechanical coupling in the common condition of purely mechanical loadings, with no direct thermal excitations, is addressed. Both local and global features of the nonlinear response are analyzed, by also comparing them with those provided by partially coupled models and/or by the underlying uncoupled mechanical oscillator, in the light of known literature results.
The manuscript is organized as follows. The unified continuum formulation of the von Kármán shear indeformable symmetric cross-ply laminated plate in a thermomechanical environment is summarized in Sect. 2, along with the dimension reduction procedure which allows to end up to the three-mode discretized model, whose coupling features are pointed out. Then, local and global nonlinear dynamics of a single-layer orthotropic plate in passive thermal regime is dealt with, by analyzing the effects of varying either inplane and transverse mechanical excitations (Sect. 3.1) or a selected thermal property of the material (Sect. 3.3), and by highlighting the non-trivial role played by the thermomechanical coupling even in a purely passive thermal setting (Sect. 3.2). The paper ends with some conclusions.
2 Unified formulation of the thermomechanically coupled problem

Continuous 2D modeling
Consider the laminated rectangular plate in Fig. 1, with thickness h and edge lengths a and b in the xand y-directions, respectively. The mid-plane of the plate coincides with the xy plane of an orthogonal Cartesian coordinate system. The plate is subjected to uniform compressive forces of magnitudes p x and p y on the plate edges, distributed mechanical transverse excitation, and thermal loadings.
The various fundamental models for the 2D thermomechanical plate depend on the selected mathematical form of the basic generalizing assumptions 3D → 2D: For the Classical von Kármán model with Thermomechanical Coupling (CTC [28]), Eq. (1a) has the following expression [32]: where u i (x, y, z,t), i = 1, 2, 3 are the components of the 3D displacement along the x, y and z directions, while u(x, y,t), v(x, y,t) and w(x, y,t) are the unknown displacement components of the 2D plate model (mid-plane displacement) along the x, y and z directions. Equations (2) result from the assumptions of a cross section plane normal to the mid-plane (shear-indeformability) and invariable along the z direction. As regards Eq. (1b), the temperature is assumed to vary linearly in the transverse direction, consistent with the mechanical assumption (2), thus writing where T (x, y, z,t) is the 3D temperature, while T 0 (x, y,t) and T 1 (x, y,t) are the unknown membrane and bending temperature components of the 2D model. Thus, the displacement {u v w} and the temperature {T 0 T 1 } are the generalized configuration variables of the CTC model. In the case of symmetric cross-ply laminate, the following governing equations in terms of configuration variables can be derived from the CTC general ones in the Appendix of [28]: (first mechanical equation) A 66 (u ,xy + v ,xx + w ,y w ,xx + w ,x w ,xy ) + A 12 (u ,xy + w ,x w ,xy ) with the mechanical and thermal boundary conditions defined by prescribing the following variables: along x = 0, a : u or N 11 conductance; T ∞ is the constant difference between the absolute temperature of the surrounding medium and the reference temperature; λ (1) 33 is the thermal conductivity along z of the first (and last) lamina.
Equations (4) are obtained starting from a unified scheme for the formulation of the thermomechanical problem of laminated plates, which integrates mechanical and thermal aspects as ensuing from the balance, configuration and phenomenological equations of a generic 2D plate model. Figure 2 summarizes such a scheme filled with the equations, variables and phenomenological quantities whose treatment ends up to Eqs. (4); for the detailed meaning of the various symbols, see [28].

Reduction to minimal models
A dimensional reduction of the 2D plate model (4) can be performed starting from the basic assumptions: along the lines reported in [28] which are herein summarized for the sake of completeness. Eqs. (6) express the 2D generalized configuration variables displacement and temperature (remind Eqs. (1)) in terms of 0D reduced configuration variables (through shape functions). In the general case, the mathematical form of Eqs. (6) for the CTC model is: where the sets of shape functions (φ u mn , φ v mn , φ w mn , φ T 0 mn , φ T 1 mn ) and the sets of reduced displacement {U mn V mn W mn } and reduced temperature {T R0mn T R1mn } components appear.
In the context of a Galerkin discretization, static condensation of the in-plane u, v displacement allows the minimal reduction of the CTC model to a system of three coupled nonlinear ODEs [28].
We refer to a simply supported rectangular laminate with immovable and isothermal edges, subjected to uniform compressive forces of magnitude p x and p y in x and y direction ( Fig. 1), respectively, whose boundary conditions read N and M being the membrane and bending components of the 2D tension variable.
For the configuration components w, T 0 and T 1 in Eqs. (7), the following single mode approximations are assumed [25,28]: where sin(πx/a) sin(πy/b) is a shape function satisfying both mechanical and thermal boundary conditions (8), that is a good idealization of also the temperature variables (so-called dome-shaped profile [33]) in many technical cases [34]. In turn, for the u and v in-plane displacements in Eqs. (7), the following basis in terms of the reduced variables W (t) and T R0 (t), ensuing from the static condensation, is used [28]: that satisfy Eqs. (4a)(4b) if neglecting therein inertia terms (which corresponds to considering frequencies of in-plane vibration much higher than frequencies of transverse vibration) and eliminating body forces. For the sake of brevity, the expressions of the coefficients u T 11 , u w 20 , (10) are not reported. Thus, Galerkin procedure can be applied to the left hand side of Eqs. (4c)-(4e). Neglecting rotational inertia and adding a linear viscous term δ w v (δ is the damping coefficient) in Eq. (4c), we end up to the following three coupled nonlinear ordinary differential Fig. 2 Unified scheme for the 2D fundamental plate models [28] filled with the equations, variables and phenomenological quantities whose treatment ends up to Eqs. (4).
Nonlinear coupling force in mechanical oscillator, deriving from the coupling terms N 11 w ,x , N 22 w ,y , N 12 w ,x and N 12 w ,y in Eqs. B1, with N i j expressed through the aliquot Nonlinear coupling force in thermal membrane oscillator, deriving from the term Linear coupling force in thermal bending oscillator, deriving from the terms a (1) = i j } in Eqs. P5 Table 1 Description of the multiphysics coupling "forces" with reference to the underlying continuous terms (Fig. 2) equations (see Eq. 61 in [28], with a convenient renumbering of the coefficients in the second equation): The expressions of the coefficients a i j and the fundamental linear frequency of the laminate ω = a 13 /a 11 are reported in the Appendix (Eqs. (14)) in terms of the lamina properties (Eq. (15)), for the case of an orthotropic single layer plate (later on considered in Sect. 3). To get an overall idea of the unified thermomechanical scheme which underlies the governing ODEs (11) of the 0D reduced model -similar to the unified scheme ( Fig. 2) underlying the PDEs (4) -, it appears useful to refer to the component balance, configuration and phenomenological equations/quantities reported in Fig. 3. A mechanical type body diagram of the reduced model is shown in Fig. 4; it schematizes in mechanical terms also the thermal aspect of the problem, and interprets the effects of the multiphysics coupling as problem-dependent configuration distortions applied to each oscillator. Table 1 describes the multiphysics coupling "forces" of the discrete model by referring to the underlying continuous terms.
The three-mode model (11) represents the simplest approximation of the underlying 2D continuum fundamental equations which still preserves the main thermomechani-cal aspects of the latter, thus enabling systematic investigation of their effects on the nonlinear vibrations of laminated plates in the absence of internal resonances among modes.

Local and global nonlinear dynamics
Model (11) allows us to consider diverse combinations of mechanical and thermal excitations, with an expected richness of nonlinear dynamic response. To get information on the role played by the thermomechanical coupling in the simple and common condition of purely mechanical excitations, no thermal loadings are considered (i.e., T ∞ = E (0) = E (1) = 0 in Eqs. (11a)-(11c)), so that the thermal field variables are entrained in the response by the directly excited mechanical displacement within a merely passive regime. In order to non-dimensionalize the model, the following transformations are introduced [25,28]: Substituting Eqs. (12) into Eqs. (11), and considering an harmonic transverse excitation F (0) 3 = F cos(Ωt), we have (the overbar upon the variables has been dropped for conve-nience): where the expressions of the coefficients a i j are reported in the Appendix (Eqs. (16)) in terms of adimensional parameters (Eq. (17)) for the case of an orthotropic single layer plate. Note that in Eq. (13a) the numbering of coefficients has been shifted with respect to Eq. (11a) because the linear stiffness coefficient a 13 now accounts for both the elastic stiffness (a 13 in Eq. (11a)) and the geometric contribution due to the axial pre-load (a 14 in Eq. (11a)), see also Fig. 4.
Reference is made to an orthotropic single layer plate of dimensions a = b = 1 m and h = 0.01 m, in epoxy/carbon fiber composite with high thermal conductivity and low specific heat, specifically designed to activate thermal processes with no computational criticalities due to the low value of the thermal stiffness in the reduced membrane equation, while also achieving a transverse strength of the considered plate such to entail physically admissible displacements. The material properties are: m 2 ·K , whose meaning is described in the Appendix (Eq. (15)). The non-dimensional coefficients of Eqs. (13a)-(13c) (now written with no overbar) become: a 12 = 0.03/η, a 13 = (1− p)/η, a 14 = 0.73/η 2 , a 15 = −0.1/η 2 , a 16 = −0.003/η 2 , a 17 = − f /η 2 , a 23 = 10 −4 /η, a 24 = 0.1, a 32 = 1.6 · 10 −3 /η, a 33 = 0.83, where η is the frequency ratio Ω /ω, f is the non-dimensional forcing amplitude, and p is the non-dimensional axial load (see Eq. (17) of the Appendix). In the following, numerical analyses are developed in primary resonance conditions, i.e. η = 1, which is known to be the most critical situation for an externally forced dynamical system; AUTO [35] and Mathematica [36] are used to describe the system periodic responses, and a C++ ad-hoc routine is implemented to obtain the cross-sections of the 4D basins of attraction.

Under varying external parameters
To analyze the effect of both the considered external loads, i.e. the transverse harmonic excitation and the axial load, a response chart summarizing the main solutions in the p − f plane is presented in Fig. 5. It shows in colors the sole regions of dominant and possibly competing 1-period solutions, with the relevant overlapping zones, and in white the regions of chaotic responses, which become dominant for  high values of axial loads and transverse forcing, while it excludes zones of possible coexistence with solutions of higher periodicity for the sake of clarity. The chart shows a quite rich response scenario for increasing values of the axial load, with multistability regions characterized by several couples of 1-period responses. A few of them, bifurcated from the fundamental cross-well solution P1 for increasing amplitude of the transverse excitation (first panel in Fig. 6(a)), are exemplarily illustrated in Fig. 6(b). Bifurcation diagrams of the displacement and of the two thermal variables are reported in Fig 6(a) in the absence of axial load, with the former pointing out the need to look at the admissible W values to be considered for avoiding the crisis of the material. Relevant strength analyses based on the maximum stress failure criterion [37] furnish the maximum non-dimensional displacement value W = 3.5 (corresponding to 3.5 cm), which entails reducing the admissible range of transverse forcing amplitude to low values. Accordingly, the following analyses are carried out for f = 1, also in order to investigate the rich multistable scenario depicted by the system for increasing values of the axial load.
The bifurcation diagram for increasing p is reported in Fig. 7, with the characterization of the periodic responses in the state plane being shown in Fig. 8(a,b), for p = 4. The results highlight that for values higher than p ∼ = 2.51 the system displays dual buckled solutions of different kinds, starting with the P1 III and P1 IV high-amplitude oscillations around the plate deformed (positive and negative) configu-(a) (b) Fig. 6 Bifurcation diagrams of displacement W , membrane temperature T R0 , and bending temperature T R1 , with increasing transverse excitation f , for p = 0. Circle: saddle-node bifurcation; Square: transcritical bifurcation (a); phase portraits of the P1 V /P1 VI (green) and P1 VII /P1 VIII (yellow) solutions, at f = 3.5 (b). rations (cyan/blue orbits in Fig. 8(a)). As the axial load increases, after the value p ∼ = 3.48, other two, low-amplitude, buckled solutions P1 I and P1 II raise up (orange/red orbits in Fig. 8(a)) and remain stable for high values of p. All these responses are generated by saddle-node bifurcations, as other 2-period and 4-period solutions that exist in narrow ranges of the control parameter. The fundamental 1-period solution P1 (gray line) representing the sole pre-buckling response of the system, exists for the whole considered range of p, becoming, in the post-buckling region, a cross-well solution oscillating around the two buckled configurations. From Fig. 7 and 8(a,b) it can be observed that the postbuckled response of the system in terms of the bending thermal component T R1 is qualitatively similar to the mechanical one, with distinct branches of periodic solutions of each pair in the bifurcation diagrams and mirrored phase-plane orbits with radially-symmetric Poincaré points, while the membrane thermal component T R0 exhibits coinciding branches and orbits, as per the self-symmetric response in the plate mid-plane. Moreover, it is one order of magnitude smaller than the bending one, likely owed to the numerical values of the coefficients a 24 and a 33 of the relevant coupling forces (see Table 1) in the thermal Eqs. (13b)-(13c): in fact, T R1 is linear with respect to the mechanical velocity, with a ratio of almost 1, while T R0 is related to the WẆ product, reduced by one order of magnitude.
To investigate the robustness of the various solutions as regards possible changes in the initial conditions of the system variables, basins of attraction are analyzed via an ad-hoc routine. Since the basins are 4-dimensional in state space, their evaluation (representation) is a computationally (topologically) demanding task, so that they are investigated through cross-sections in the plane of the mechanical variables (W,Ẇ ) by considering naturally vanishing initial conditions of the thermal variables, i.e. T R0 (0) = T R1 (0) = 0.0, as per their solely passive role in the response. Fig. 8(c) shows the system basins of attraction at f = 1, p = 4. Here it is evident that the most robust solutions are the P1 I and P1 II responses, which display compact basins of attraction (orange/red in Fig. 8(c)) around the buckled configurations. The other solu- Fig. 9 Comparison of the bifurcation diagrams as a function of the axial load between the thermomechanical model (red curves) and the purely mechanical one (black curves). Circle: saddle-node bifurcation; Square: transcritical bifurcation; Diamond: period-doubling bifurcation.
tions, conversely, exhibit basins with markedly fractal characteristics, corresponding to an almost zero robustness.
It is also of interest to compare the results obtained by the coupled thermomechanical model (13a)-(13c) with those furnished by the underlying mechanical model with no thermal contributions (represented by Eq. (13a) with a 15 = a 16 = 0), which, after buckling, corresponds to an uncoupled doublewell Duffing oscillator. Fig. 8(d) shows the basins of attraction of the latter, to be compared with the cross-section in Fig. 8(c) relevant to the thermomechanical model, while Figure 9 compares the bifurcation diagrams of the two systems (red for the thermomechanical model and black for the purely mechanical one) for increasing axial load. An almost exact coincidence of the local and global responses of the two models is noticed, in terms of both stability of the main periodic solutions and robustness of the basins of attraction, meaning that in the passive thermal regime the coupling with temperatures does not alter the steady mechanical response but only introduces non-null steady values of the thermal variables in the overall system response.

The system with nearly zero mechanical stiffness
Yet, a non-trivial role of the thermomechanical coupling can be highlighted even in a purely passive thermal setting, i.e. in the absence of whatever thermal excitation. To this aim, it is worth analyzing the global features of the response occurring around p = 1.0, which entails zeroing of the a 13 coefficient in Eq. (13a), i.e. vanishing of the mechanical linear stiffness.
Comparing the cross-sections at T R0 = T R1 = 0.0 of the 4D basins of attraction of the coupled model ( Fig. 10(a)) and the 2D basins of the uncoupled Duffing oscillator (Fig. 10(b)), meaningful differences can be observed. They pertain mostly to the two P2 solutions, which -though playing an important role already in the global response scenario of the symmetric Duffing oscillator, where they strongly compete with the P1 solution at perfect vanishing of the linear stiffnessbecome definitely more robust for the coupled model: this is clearly witnessed by their basins filling nearly completely the right and left portions of the considered window of initial conditions in Fig. 10(a), while also affecting the window core.
To better grasp the rationale behind these differences, the effect of the two thermal contributions present in Eq. (13a) on the global response scenario is investigated separately (Fig. 11). Setting a 16 = 0 (i.e., artificially imposing vanishing of the T R0 membrane coupling, right panels) entails no differences with respect to the scenario obtained with the fully coupled model (left panels), thus making apparent that the periodic modulation of the linear stiffness of the mechanical oscillator induced by the parametric coupling term a 16 T R0 W is actually of minor significance, even where (i.e., for p = 1.0) it represents the sole thermally-induced stiffening effect. The situation is completely different if a 15 is set to zero, namely, if vanishing of T R1 bending contribution is imposed, which entails a global scenario of response (mid-right panels) substantially identical to the one exhibited by the purely mechanical oscillator (mid-left panels), with the dual P2 solutions already existing for p = 0.8 (as also suggested in [38]) playing a major role for p = 1, which corresponds to the Ueda oscillator [39]. In this case, it is important to also look at the time histories of the mechanical and thermal components of the response for the fully coupled model (red curves in Fig. 12) and the coupled one with no bending contribution in the mechanical oscillator (gray curves in Fig. 12), obtained starting from fixed mechanical initial conditions (−3.4, 0), which belong to the green crosssection basin of the P2 solution ( Fig. 11(b), left panel) and to the gray basin of the P1 solution ( Fig. 11(b), mid-right panel), respectively. In the latter case (a 15 = 0), the steady coupled response oscillates with period one as in the uncoupled Duffing oscillator (Fig. 11(b), mid-left panel, and gray response in Fig. 10(c)), notwithstanding the modulation of the linear stiffness herein entailed by the T R0 coupling term, while in the former case (a 15 = 0), the steady dynamics is 2-periodic (green response in Fig. 10(c)).
The investigation of the transient part of the time histories reported in Fig. 12(b) reveals that the different behavior of the two systems has to be ascribed to the very initial steps of the responses evolution, when the T R1 component of the response, which is passively activated by the coupling with the mechanical velocity, oscillates around a non-vanishing mean value, thus providing -in the fully coupled modela feedback time-varying contribution a 15 T R1 (t) in the mechanical oscillator equation which includes a constant indirect excitation, overall entailing a more robust P2 response of the latter.
The crucial role of the T R1 mean value is confirmed when analyzing the purely mechanical Duffing oscillator with varying constant excitation, i.e. when a 16 = 0 and T R1 (t) = T R1 = c in Eq. (13a). The relevant bifurcation diagram in Fig. 13 reveals that the coexistence of the P1 solution with the two P2 responses occurs for −1.58 < c = T R1 < 1.58, while the sole P2 solutions last outside this range. Such outcomes are in agreement with what already highlighted in the literature about the symmetric Duffing oscillator with no linear stiffness [39], for which the addition of a constant excitation to the harmonic external one makes the mechanical oscillator markedly asymmetric [40,41], thus enhancing the role of the P2 solutions. This is equivalent to what happens, more generally, for a Helmholtz-Duffing oscillator [40,42] close to the stability limit of 1 2 -subharmonic resonance, where the P2 solution becomes the governing one.
However, for the coupled thermomechanical model under analysis, these results have to be framed in the transient evolution of the system response, since at the steady regime the constant excitation provided by the bending thermal variable to the mechanical oscillator is certainly negligible (i.e. the T R1 component settles on oscillations around a nearly zero equilibrium, as shown in the right panel of Fig. 12(a)). In this sense, the T R1 contribution to the mechanical variable during the temporal evolution of the system can be interpreted as a decreasing constant excitation starting with a value close to the mechanical initial condition (which directly affects the bending thermal initial process) and ending up to zero. The comparison between the bifurcation diagram of the Duffing system and the crosssection basin of the coupled model reported in Fig. 13 helps to clarify this conjecture. If the starting mechanical initial condition (and consequently the T R1 initial equilibrium position, i.e. the value of the constant excitation) belongs to the range of coexistence of the P1 and P2 solutions, the coupling with the bending thermal component is manifested in the slight modification of the basins organization, without anyway qualitatively altering their overall topology. In contrast, when the mechanical initial condition implies the shift of the transient bending thermal oscillations around highly varied configurations, the associated (time decreasing) constant excitation to the mechanical oscillator forces the trajectories to settle on one of the P2 attractors (depending on the sign of the initial condition). Once the trajectories reach the P2 attractors, they remain on them as the T R1 contribution decreases, since they represent the sole periodic solutions that the system exhibits for both high and zero constant T R1 values. Of course, the Duffing system depicts also other periodic solutions for high constant excitations, e.g. the stable P4 solutions arising from the period doubling of the P2 responses, represented with magenta curves in Fig. 13. For even higher constant excitations, ensuing chaotic responses occur for both the Duffing oscillator [43] and the fully coupled thermomechanical model, however in the latter case affecting the solely transient regime. Indeed, such solutions do not last when the c constant excitation approaches the zero value in time, so that, during the temporal evolution of the T R1 variable in the thermomechanical model, the relevant attractors disappear and the trajectories move to the enduring P2 responses which become the steady ones for the coupled system.

Influence of intrinsic thermal parameters
Until now, attention has been devoted to analyze the dynamical behavior of the system under variations of external parameters, like the transverse and axial loads, which represent the most easily changeable and controllable properties of the system from a design and operational point of view. However, analyzing the dynamical effects of other, intrinsic, parameters can be also particularly interesting to the aim of properly designing the plate according to the behavior that has to be ensured. Several analyses as a function of the various physical properties have revealed that, among all others, setting the thermal expansion α 2 even beyond the range of values typical of presently known materials produces a substantial effect on the system dynamics, as shown in Fig. 14 and Fig. 15.
In Fig. 14(a) the bifurcation diagrams of the thermomechanical model for increasing α 2 are presented, for f = 1, p = 4. For the α 2 value of the formerly described epoxy/carbon fiber composite, i.e. α 2 = −1.91 · 10 −6 1 K , the response scenario is the one represented in the phase portraits and basins of attraction of Fig. 8, which is characterized by five 1period and two 4-period solutions. Starting from this situation, the increase (or decrease) of the thermal expansion causes a reduction of the multistability which disappears at α 2 = 6.95·10 −3 , where the only stable solutions are the lowamplitude buckled P1 I and P1 II responses. Such local behavior is confirmed by the evolution of the basins of attraction, represented again by cross-sections in the (W,Ẇ ) plane with T R0 = T R1 = 0.0, and reported in Fig. 14(b) for increasing α 2 . The marked fractality which affects all basins of attraction is strongly reduced as the thermal expansion increases, and the phase portrait becomes more and more regularly organized, with a clear separation of the initial conditions leading to the various attractors. Finally, for α 2 > 6.95 · 10 −3 , the sign of the initial displacement determines the buckled response to be achieved.

Summary and future developments
The main steps of a unified 2D continuum formulation of the nonlinear dynamic problem for a von Kármán shear indeformable symmetric cross-ply laminated plate in a thermomechanical environment have been presented, along with the ensuing reduction procedure ending up to a three-mode discretized model. This consists of one mechanical and two thermal ordinary differential equations, with unknown transverse displacement, membrane temperature, and bending temperature, and allows us to account for a variety of excitation conditions and material properties. The reduced model preserves the main nonlinear dynamic features of the underlying continuum formulation, and is suitable for the analysis and understanding of the basic effects of thermomechanical coupling on the finite amplitude vibrations of geometrically nonlinear laminates, in the absence of internal resonances.
In particular, it enables to address both the case of thermal dynamics passively entrained by the coupling in the presence of sole mechanical excitations, and the case of thermal dynamics directly activated by some coexisting thermal source. The present study has focused on the former, and indeed simpler, case, with the aim to analyze the overall features of the system nonlinear response, while also highlighting fundamental aspects associated with the presence of thermomechanical coupling.
The three-mode model has been thus used for systematic investigation of the local and global nonlinear dynamics of a single-layer orthotropic plate, subjected to sole axial loads and transverse harmonically varying excitation. Though the two thermal variables are solely dragged into the system overall response by the directly excited mechanical variable, they exhibit different qualitative and quantitative features, with the relevant coupling terms affecting the mechanical displacement in a different way.
Response scenarios under varying in-plane and/or transverse excitation have been systematically obtained through numerical simulations, along with the effect of a selected thermal property of the material, by also comparing them with the response provided by partially coupled models and/or by the underlying uncoupled mechanical oscillator. Even in the solely passive thermal setting herein considered, for the nearly zero stiffness configuration of the system -also analyzed in the light of literature results for archetypal oscillators -the transient evolution of the thermomechanically coupled response has been shown to play a meaningful role in steadily modifying the system global dynamics with respect to that of the uncoupled oscillator.