Thermomechanical Coupling and Transient to Steady Global Dynamics of Orthotropic Plates

Di ﬀ erent reduced order models of thermomechanically coupled von Kármán shear indeformable plate with prescribed linear temperature along the thickness are comparatively investigated in terms of local and global dynamics exhibited in active thermal regime, under harmonic transverse and constant axial mechanical excitations. Two-d.o.f. one-way coupled models are then used for “economical” yet reliable numerical investigation of the nonlinear dynamic response in terms of the mechanical variable and of the dominant (membrane or bending) thermal variable, under the corresponding thermal excitations inducing variably rich scenarios of buckled mechanical response. Attention is focused on the important role played by global dynamics in unveiling the meaningful e ﬀ ects entailed on the structure steady mechanical response by the variably slow thermal transients taken into account by the thermomechanically coupled model.


Introduction
Thermomechanical coupling in laminated plates has been recently a matter of refined modeling at both the continuum and discretized level Rega, 2014, 2017), with a view to highlighting the relevant effects on the structure nonlinear dynamic response. Dealing with a reduced order model (ROM) of von Kármán shearindeformable plate with one mechanical (mid-plane transverse displacement) and two thermal (membrane and bending) variables, the effects of thermoelastic coupling on mechanical vibrations in passive thermal conditions (i.e., with no thermal excitations) have shown to be relatively minor (Saetta and Rega, 2014;Settimi et al, 2017), although the thermal variables dragged into the structure response by the directly excited mechanical variable exhibit distinct qualitative and quantitative features, with the relevant coupling terms affecting the displacement in a different way (Settimi et al, 2017). The role of different variables/terms of coupling in the mechanical and thermal equations was also addressed (Saetta and Rega, 2014), along with the possibility to refer to variably simplified models for catching some main aspects of the thermally affected mechanical response (Saetta and Rega, 2016). Based on preliminary results, no seemingly significant effects of thermoelastic coupling on mechanical vibrations were detected, even in the presence of a non-vanishing heat flow or of explicit thermal sources giving rise to a condition of active thermal dynamics (Saetta and Rega, 2014).
However, more extended analyses conducted in the regime of full thermomechanical coupling have originated a meaningful set of outcomes where the much slower time scale over which thermal phenomena develop with respect to mechanical ones entails non-trivial steady effects on the structure overall response (Settimi et al, 2018). This influence has been highlighted by means of a systematic numerical investigation of the structure global dynamics, pursued through properly selected two-dimensional (2D) cross-sections of the actual 4D basins of attraction.
Global nonlinear dynamics plays a fundamental role in the analysis, control, and design of all engineering systems, as recently highlighted for a variety of them both with reference to specific mechanical/dynamical issues (Settimi and Rega, 2016b,a;Carvalho et al, 2017) and within a more general perspective (Rega and Lenci, 2015). For the considered thermomechanically coupled plate, even in passive thermal regime global analysis has allowed to detect and understand a substantially different scenario of dynamic response occurring in specific mechanical conditions (i.e., vanishing mechanical stiffness) between coupled and uncoupled models, due to the meaningful role played in the former just by the slow transient evolution of the system thermal dynamics (Settimi et al, 2017). When considering also thermal excitations in addition to mechanical ones, these slow transient effects on the steady structural response are seen to be more pervasive and systematic, as highlighted by the comprehensive investigation in Settimi et al (2018).
The present work aims at overviewing some main effects of thermomechanical coupling in active thermal regime and in the presence of also a prescribed axial mechanical excitation, whose possible buckling effects may combine with thermal ones. Yet, for coupled systems, one main issue consists of preliminarily evaluating the relative importance of the various terms and variables of coupling in the governing equations, to the aim of a satisfactory -and more or less comprehensivedescription of the system response. Dealing in particular with the thermomechanical interaction, fully (Yeh, 2005) or (to a different extent) partially (Ribeiro, 2007;Alijani et al, 2011;Shen and Xiang, 2014;Thanh et al, 2017) coupled analyses can indeed be pursued, with possibly different outcomes depending on the specifically considered situation or with an incomplete description.
Accordingly, the first part of this work (Sect. 2) is devoted to comparatively investigating some main local and global dynamical aspects of the structure nonlinear response as obtained with different ROMs ranging from the original, fully coupled, one with three degrees-of-freedom (d.o.f.) down to the simplest uncoupled mechanical model in which the presence of a coexisting active thermal environment is accounted for by merely including its known steady effect in the relevant equation of motion. This allows to properly identify the actual need and/or the possible added value (if any) of all workable models and to select a suitable one to be used for reliable yet "economical" systematic investigations, based on the structural response aspects which one is interested in. Section 3 highlights the role played by global analysis in unveiling the effects of thermal transient on the structure steady dynamics, by considering two different kinds of thermal excitation entailing variably rich scenarios of buckled mechanical response.

Reduced order models and their comparative outcomes
The thermomechanically coupled, geometrically nonlinear model here considered refers to a rectangular laminated plate in Fig. 1, of thickness h, and edge lengths a and b in the x and y directions, respectively. The plate is subjected to uniform compressive forces p x and p y on the plate edges, to distributed harmonic transverse mechanical excitation, and to thermal loadings. The model is derived within a unified modeling framework integrating mechanical and thermal aspects which is presented in Saetta and Rega (2014); Settimi et al (2017), to be referred for the formulation details. By considering classical displacements with von Kármán nonlinearities, along with a correspondingly consistent linear variation of the temperature along the thickness, and in the contest of a minimal Galerkin discretization in conditions of no internal resonance between the transverse modes of the laminate, single-mode approximations are assumed for the transverse displacement w and the temperatures T 0 and T 1 , while the in-plane displacement components are statically condensed. The subsequent nondimensionalization allows to obtain three coupled nonlinear ordinary differential equations for a simply supported orthotropic single-layer plate with immovable and isothermal edges: in terms of the nondimensional reduced variables W (deflection of the center of the plate), T R0 (membrane temperature) and T R1 (bending temperature). Coefficients a i j are defined in Settimi et al (2017) and Settimi et al (2018), while the thermal excitations are represented by the constant thermal difference between plate and environment T ∞ , and by the membrane and bending excitations e 0 and e 1 derived from body source thermal energies, whose distribution along the plate thickness is reported in Fig. 1.
In passive thermal regime, under specific conditions, the mechanical activation of thermal variables thanks to the coupling terms in the relevant equations may furnish contributions back to the mechanical equation able to modify the basins of attraction organization, highlighting the need to consider the two-way model to grasp the actual mechanical dynamics (Settimi et al, 2017). In the active thermal regime considered here, there are hints from Saetta and Rega (2014) about the low role played in the mechanical response by both the non-directly excited thermal variables/equations (which govern the number of essential thermal degrees of freedom) and the mechanical coupling terms in thermal equations (which govern the two-way coupling). To deeper investigate this aspect, and in order to detect the simplest model to be used to grasp the main thermomechanical features of the system response, it is of interest to develop a critical analysis of the various models derivable from the general system (1). To this aim, reference is made to an epoxy/carbon fiber composite plate of dimensions a = b = 1 m and h = 0.01 m. The relevant mechanical and thermal properties are reported in Settimi et al (2018), and furnish the following values of the equations coefficients: a 12 = 0.0592,ā 13 = (1 − p),ā 14 = 0.6827,ā 15 = −0.3674,ā 16 = −0.9658, a 17 = − f ,ā 22 = −1.4507,ā 23 = 9.1137 · 10 −5 ,ā 24 = 1.01 · 10 −4 , a 25 = −0.997719,ā 32 = 7.8735 · 10 −4 ,ā 33 = 8.8714 · 10 −4 ,ā 34 = −12 obtained in primary resonance conditions, which represents the most critical situation for an externally forced system, as a function of the precompression p and of the harmonic forcing amplitude f . With respect to the latter, preliminary analyses in presence of the sole mechanical excitations (Settimi et al, 2017) show the possibility to catch the main dynamical aspects of the system response also for low values of the amplitude f , which is therefore fixed to f = 1, and the ability of the precompression in inducing post-buckling behaviors for values higher than p = 2.52. In order to critically discuss the effects of the thermal excitations on the mechanical behavior of the plate, yet dealing with parameter values corresponding to physically admissible quantities, the precompression load is set to p = 2.51, i.e. a mechanical incipient but not yet triggered buckling configuration.
With the same objective, the thermal excitation is selected to reproduce the condition able to most affect the mechanical response of the system, i.e. to more quickly provide its contribution into the mechanical equation thanks to the coupling termsā 15 orā 16 . Observing the linear thermal stiffnessesā 23 andā 32 , it results that the membrane temperature evolution is much slower than the bending one (ā 23 /ā 32 ∼ = 0.12), thus slowing down the influence of the membrane excitations (T ∞ and e 0 ) on the mechanical dynamics (Settimi et al, 2018). For these reasons, the model comparison is developed by considering the presence of a bending excitation e 1 with a linear variation along the thickness which reproduces a cooling of an external surface and an equivalent warming of the other.
The definition of the various, gradually simplified, models to be compared moves from the most general fully coupled, or in different words two-way (from thermal to mechanical, and from mechanical to thermal) coupled, 3 d.o.  (1a) (1 model) without membrane contributions (ā 16 = 0) and with the bending thermal excitation taken into account by setting the T R1 value corresponding to the mean steady state response obtained from linearization of Eq. (1c): Local and global dynamics of the five above described models are compared in terms of bifurcation diagrams, basins of attraction and temporal evolution of selected trajectories, as shown in Figs. 2-7. The outcomes of the bifurcation diagrams in terms of minimum and maximum value of the mechanical displacement W , obtained with the fully coupled 3-2 model and reported in Fig. 2, show that, starting from a pre-buckling scenario for e 1 = 0 with the sole existing P1 (thin black) cross-well solution, the addition of a positive/negative bending excitation e 1 proves to be able to generate buckled responses, though confined around only one positive/negative equilibrium (i.e., in one of the two potential wells, in global dynamics terms), depending on the sign of e 1 . Moreover, consistent with the physically expected effect of changing the sign of e 1 , the overall scenario of the system mechanical response is antisymmetric, with the two stable low-amplitude buckled solutions (P1 I thick/P1 II thin light gray curves) existing in complementary ranges of the excitation, and with the high-amplitude ones (P1 III thick/P1 IV thin gray curves) being confined in closer complementary ranges. As a consequence, for sufficiently high values of e 1 , the system displays period-1 buckled solutions oscillating around one sole varied configuration, together with cross-well period-2 solutions P2 (thick black curves) confined in limited ranges of the positive/negative excitation values, thus highlighting the possibility to exploit the e 1 excitation to force the mechanical buckling around a selected equilibrium, according to its sign. The described behavior is exactly reproduced by all the four simplified models 3-1, 2-2, 2-1, 1, pointing out that the steady dynamics of the periodic responses, or in other words the system attractors, is insensitive to both the possibly different thermomechanical transient and coupling.
Instead, differences among the models responses can be caught if looking at the bifurcation diagrams in terms of the bending and membrane thermal variables T R1 and T R0 , displayed in Fig. 3. In fact, while the response of the directly activated bending variable T R1 (Fig. 3a) cannot be followed by the mechanical 1 model (though its value, as known parameter, can be deduced from Eq. (2)), the behavior of the dragged membrane temperature T R0 (Fig. 3b) can be described by the sole 3-2 and 3-1 models which include both the thermal variables. Moreover, with respect to the latter, the fully coupled 3-2 model demonstrates to be the only one able to correctly grasp its behavior, since the 3 d.o.f one-way model, although considering the membrane variable, furnishes identically null solutions due to the non-activation of Eq. (1b), which remains uncoupled and unforced. It is worth noting that the response in terms of membrane temperature is not affected by the sign of the excitation, i.e. it is symmetric with respect to the e 1 = 0 line, coherently with T R0 describing the thermal behavior in the plate mid-plane. Differently, the response of the bending thermal variable T R1 is organized along a straight line, due to the almost linear nature of the relevant thermal equation.
When considering also possible variations of the mechanical initial conditions, the comparison between coupled (i.e. 3-2,3-1,2-2,2-1) and uncoupled 1 models highlights meaningful discrepancies. Figure 4 displays the basins of attraction (or the relevant cross sections in the mechanical plane, depending on the system d.o.f.'s) realized, with trivial values of the thermal variables, after the application of a bending excitation e 1 = 3 × 10 −4 , corresponding to a power density linearly distributed along the thickness z of 8687.3 × 10 3 × z kW/m 3 . The coupled models, whose outcomes are coinciding and reported in Fig. 4a, provide discordant indications with respect to those furnished by the bifurcation diagrams, with a monostable behavior characterized by the sole (thin gray) high-amplitude buckled P1 IV solution, without any evidence of the expected basins of attraction of the thin light gray P1 II and of the thick black P2 responses. Conversely, the results obtained from the mechanical 1 model coherently identify all the three different basins with dominance of the P1 II response and a markedly fractal organization. The reasons of this disagreement have to be sought in the effect of the thermomechanical coupling together with the contemporary presence of slow and fast dynamics. To better understand the behavior of the coupled systems, in fact, it is worth noticing that the slowness of the bending thermal transient (which is taken into account by all the four models due to the combined presence of Eq. (1c) and termā 15 ) causes the contribution of the e 1 excitation to be supplied gradually into the mechanical equation by means of the coupling term related to T R1 . On the other hand, the mechanical vibration is much faster than the thermal one and its transient, needed for reaching a stable solution, is very short. From a phenomenological viewpoint, it appears possible to neglect the mechanical transient and to look only at the attractor of the system, whose evolution with increasing values of the thermal bending excitation from zero to the selected value e 1 = 3 · 10 −4 can be followed in the bifurcation diagram of Fig. 2. For low values of e 1 , where the thin black P1 attractor is stable, the mechanical response is cross-well at least in its first initial steps, before possibly jumping to the coexisting gray buckled solution P1 IV upon its onset. But  Fig. 3 Bifurcation diagrams of the thermal variables T R1 and T R0 as a function of thermal excitation e 1 . a Bifurcation diagrams of the bending temperature T R1 furnished by the 3-2, 3-1, 2-2, 2-1 models. b Bifurcation diagrams of the membrane temperature T R0 furnished by the 3-2 and 3-1 models. Circle: saddle-node bifurcation; Diamond: period-doubling bifurcation.
when the bending excitation reaches a value (e 1 ∼ = ±5.6 · 10 −5 ) providing a bending thermal variable T R1 in the mechanical equation such to instabilize the P1 response via period doubling bifurcation (see Fig. 2a), the mechanical trajectories are more likely to swiftly jump onto the P1 IV buckled response. Since this solution is stable also after the arise of the further (thin light gray) buckled solution P1 II , the system response steadily remains on it in the whole considered cross-section of mechanical initial conditions. In contrast, the thermal transient is totally ignored in model 1, in which only the final value of the bending temperature evolution is added into the mechanical equation, and the system is able to display the basins of all the three solutions detected by the bifurcation diagram.
The effects of the various thermomechanical coupling terms can be further discussed by observing the temporal evolution of a single trajectory with fixed trivial initial conditions (W =Ẇ = T R0 = T R1 = 0) followed by considering the five proposed models. Looking at the mechanical response in Fig. 5a, the outcomes confirm the coinciding behavior of the four coupled models which reach the steady high-amplitude gray P1 IV solution, in contrast with the low-amplitude light gray P1 II oscillation detected by the uncoupled 1 system, due to the fact that the selected mechanical initial condition belongs to different basins of attraction, as shown comparing the graphics of Fig. 4.
The evolution of the directly activated bending variable displayed in Fig. 5b, possibly followed only by the four coupled models, allows to point out the role played by the two-way coupling. In fact, despite the congruent steady value reached by all models, the periodic motion can be grasped by the sole gray 3-2 and red 2-2 ones, thanks to the oscillatory contribution provided by the mechanical response by means of the coupling term related toā 33 , while the one-way 3-1 and 2-1 models show a stationary solution. Moreover, due to the coinciding outcomes of the 3 d.o.f. Fig. 4 Basins of attraction in the mechanical plane, for p = 2.51 and e 1 = 3 × 10 −4 . a Planar section of the basins of attraction for the 3-2, 3-1, 2-2, 2-1 models. b Basins of attraction for the 1 model. Light gray: basin of the low-amplitude buckled solution P1 II ; Gray: basin of the highamplitude buckled solution P1 IV ; Black: basin of the period-2 solution P2. and of the 2 d.o.f. (3-2≡2-2, 3-1≡2-1) models, it can be observed that the membrane thermal variable does not affect the behavior of the bending temperature.
Finally, as previously stated (see Fig. 3b), the dynamics of the dragged membrane temperature can be followed only if employing the fully coupled 3-2 model (see Fig. 5c), which is able to originate a non-null response by means of theā 24 coupling term, while the other model that also contemplates the presence of the T R0 variable, i.e. the 3-1 model, provides null results due to the non-activation of Eq. (1b).
The characterization of the detected responses in the mechanical and thermal planes is reported in Fig. 6 in terms of phase portraits and Poincaré maps; here, the amplitude difference between the mechanical response identified by the coupled models and that determined by the mechanical 1 system is highlighted (Fig. 6a), and the very small amplitudes of the thermal variables are described (Figs. 6b,c).
In the light of the obtained results, and with the aim to focus the attention on the mechanical behavior of the thermomechanical plate, it can be thus convenient to work with the reduced 2 d.o.f. one-way 2-1 model, which is capable to correctly describe the response in terms of both transversal displacement and stationary bending temperature, while requiring a minor computational effort and a dimensionally Time histories of displacement and thermal variables, for p = 2.51, e 1 = 3 × 10 −4 and trivial mechanical and thermal initial conditions. a Temporal evolution of the mechanical displacement W for 3-2, 3-1, 2-2, 2-1, 1 models. b Temporal evolution of the bending temperature T R1 for 3-2, 3-1, 2-2, 2-1 models. c Temporal evolution of the membrane temperature T R0 for 3-2, 3-1 models. Gray: 3-2 model; Green: 3-1 model; Red: 2-2 model; Blue: 2-1 model; Cyan: 1 model. lowered mathematical system with respect to the fully coupled one (1). Clearly, with this model, information about the oscillatory nature of the bending temperature, as well as any indication about the behavior of the membrane one, are neglected, which have however proved to be negligible.
As a conclusive analysis, it is of interest to comparatively summarize the transient dynamics of the mechanical response displayed by all the considered models, whose outcomes are presented in Fig. 7, which underline the influence of the coupling terms on the length of the non-stationary evolution. As shown, the shortest transient is relevant to the mechanical cyan 1 model, where the thermal contribution furnished by the bending excitation is entirely provided into the mechanical equation at the beginning of its time history. The two one-way models (green 3-1 and blue 2-1) display the same transient duration, due to the fact that the membrane thermal variable T R0 , although considered in the 3-1 model, is not activated by the excitation or by the coupling term, therefore having no influence on the displacement evolution. Instead, when the coupling terms are present also into the thermal equations, i.e., when considering the two-way gray 3-2 and red 2-2 models, the non-stationary dynamics lengthens, with the fully coupled 3-2 model which exhibits the longest transient thanks to the very slow contribution supplied by also the dragged membrane variable (whose slowness can be grasped by comparing the three graphics of Fig. 5). Fig. 6 Phase portraits and Poincaré maps in the mechanical and thermal planes, for p = 2.51, e 1 = 3 × 10 −4 and trivial mechanical and thermal initial conditions. a Phase portraits in the mechanical plane (W,Ẇ ) for 3-2, 3-1, 2-2, 2-1, 1 models. b Phase portraits in the thermal plane (T R1 ,Ṫ R1 ) for 3-2, 3-1, 2-2, 2-1 models. c Phase portraits in the thermal plane (T R0 ,Ṫ R0 ) for 3-2, 3-1 models. Gray: 3-2 model; Green: 3-1 model; Red: 2-2 model; Blue: 2-1 model; Cyan: 1 model.

Bending
As shown in the previous section, the cross sections of the basins of attraction furnish, for the coupled models, contrasting results with respect to the outcomes provided by the bifurcation diagrams, due to the effect of the thermal transient into the evolution of the mechanical displacement.
Referring to the reduced 2 d.o.f. one-way 2-1 model formerly introduced, the matter can be better understood by also looking at different planar cross-sections of the 3D basins of attraction. In fact, the basin considered in Fig. 4a for null values of the thermal variable is certainly the reference natural one if thinking in purely mechanical terms; however, a more comprehensive description of the basins organization in the state space can only be obtained by considering planar cross-sections for also non-trivial values of the thermal variables. This is necessary mostly if being interested in grasping the final outcomes of the dynamics started with given initial conditions of the variable specifically governing the system response, which is here the bending temperature.
With application of the excitation value considered in Fig. 4, and looking at the final outcomes of the dynamics started with given T R1 initial conditions, Figure 8a points out that relevant increasing values (Fig. 8b) succeed in catching the presence of also other basins of attraction, up to reproducing the response of the mechanical 1 model when the initial condition is set to the relevant regime value, as highlighted by the coincidence of results between Fig. 4b and the right panel of Fig. 8a. This is due to the progressive shortening of the transient dynamics as the thermal initial W t Fig. 7 Transient dynamics of the mechanical response, for p = 2.51, e 1 = 3 × 10 −4 and trivial mechanical and thermal initial conditions. Gray: 3-2 model; Green: 3-1 model; Red: 2-2 model; Blue: 2-1 model; Cyan: 1 model. conditions become closer to the steady value to be attained, as shown in Fig. 8b, which corresponds to reduce the gap that the bending thermal variable has to cover.
To summarize this behavior, a cross section of the basins of attraction in the (T R1 ,Ẇ ) plane is presented in Fig. 8c, which is obtained by fixing the mechanical initial condition to W = 1.75, within the buckled P1 II light gray basin at T R1 = 3 (see the mid panel of Fig. 8a). For initial T R1 lower than 3, the response of the coupled model always settles on the high-amplitude buckled P1 IV gray basin, due to the slow thermal contribution into the mechanical equation which is insufficient to move the response towards the other, yet existing, P2 and P1 II attractors. For 3 < T R1 < 5, the arise and enlargement of the low-amplitude P1 II light gray basin is highlighted, while for high T R1 initial conditions the latter becomes the only existing solution for the system, according to the outcomes of the bifurcation diagrams of Fig. 2 showing the P1 II response as the only possible one for high values of the thermal excitation. In fact, setting a high T R1 initial condition corresponds to provide, in the first step of the system dynamical evolution, a high contribution into the mechanical equation by means of the coupling term related toā 15 . As a consequence, the displacemenṫ response is initially moved into the monostable range characterized by the P1 II solution, which therefore attracts all the trajectories regardless of the chosen initial velocity. Due to the fact that such solution is stable also when the thermal variable lowers and stabilizes around its steady value, it represents the only possible response of the system in the considered range.

Membrane
To critically discuss the effects of a membrane thermal excitation on the dynamical behavior of the system, two different kinds of thermal forcing can be considered according to the model, i.e. a time constant thermal difference T ∞ between plate and surrounding medium, which activates pure thermal convection on the external surfaces and pure internal thermal conduction, or a time-independent thermal excitation e 0 , constant along the thickness, which can be physically obtained by exploiting the Joule effect due to the current passage in the plate (thanks, for example, to the insertion of conductive metallic fibers or carbon nanotubes into epoxy matrix, aimed at enhancing its electrical conductivity). From linearization of Eq. (1b), and the subsequent substitution into Eq. (1a), a direct relation can be obtained among the two thermal excitations and also the mechanical precompression p: It can be observed that e 0 and T ∞ produce the same effect on both the thermal and the mechanical vibrations, and, due to the coupling into the mechanical equation by means of a linear T R0 term, their action on the displacement variable corresponds to that of the precompression p. Such correspondence is stressed in the bifurcation diagram of the W displacement reported in Fig. 9a, obtained as a function of the two different (thermal e 0 and mechanical p) excitations, and, as done in the previous sections, by matching the zero of the thermal excitation with a p value corresponding to the incipient buckling (p = 2.51).
The results show the ability of e 0 in inducing a multistable post-buckling scenario, with the passage from the sole pre-buckling thin black P1 solution to the contemporary presence of two couples of buckled responses oscillating around varied positive/negative configurations of the plate, corresponding to low-amplitude buckled solutions (thick P1 I /thin P1 II light gray curves) and to high-amplitude solutions (thick P1 III /thin P1 IV gray curves) coexisting in a relatively wide range of the control parameters. The bifurcation diagram highlights also a globally symmetric behavior of the solutions of each couple with respect to the trivial equilibrium, coherently with the constant distribution of the membrane excitation along the plate thickness, and in contrast with what obtained from the application of the bending excitation previously described. Moreover, with respect to the results of Fig. 2, it is worth underlining that in this case the pre-buckling solution P1 is stable also in the post-buckling regime, where it becomes a high-amplitude cross-well solution oscillating around both the varied plate equilibria.
In particular, when a positive thermal excitation of e 0 = 10 −4 is applied, the system steady response is characterized by a multistable behavior including two pairs of buckled solutions, P1 III /P1 IV and P1 I /P1 II , coexisting with the pre-buckling P1 response. However, the cross-section of the 3D basins of attraction with trivial values of the thermal variables displays a completely different scenario, as shown in Fig. 9b, in which the system behavior results to be monostable. The discrepancies between expected and obtained behaviors of the model can be understood by considering again the long thermal transient needed by the membrane temperature to attain its final steady value.
In fact, it entails, in the mechanical equation, a slow thermal contribution to the system overall stiffness up to reaching the value necessary to achieve the buckled configuration. As a consequence, the mechanical response in the very first initial steps of its temporal evolution falls onto the pre-buckling solution (black basin), which represents the only stable response for the system when no thermal excitation is applied, irrespective of the chosen mechanical initial conditions. Since it represents a robust attractor in the whole range of parameters here considered, the trajectories already settled on it do not modify their behavior also when other buckled responses arise in the 3D state space.
Also in this case, as seen before, changing the thermal initial conditions causes the arise and enlargement of the buckled P1 III and P1 IV basins of attraction (gray basins), and the subsequent birth of the low-amplitude buckled P1 I and P1 II basins (light gray basins), as shown in Fig. 10a for increasing values of T R0 . This behavior is well represented in Fig. 10c, where the section in the (T R0 -Ẇ ) plane is reported. The outcomes confirm the dependency of the mechanical response on the thermal initial values, and highlight a response similar to that obtained when analyzing the effects of the bending excitation (see Fig. 8c). However, here three regions can be detected, corresponding to qualitatively different responses of the system. For low values of T R0 , the mechanical displacement settles on the pre-buckling (black) P1 solution, while the basins of the high-amplitude buckled solutions P1 III and P1 IV appear and enlarge their compact part for 0.15 < T R0 < 1, anyway coexisting with the pre-buckling one and showing a markedly fractal organization. For T R0 higher than 1, the low-amplitude buckled responses P1 I and P1 II (light gray basins) also arise and widen, with evident dominance of the one oscillating around the positive configuration, according to the chosen positive mechanical fixed initial condition W = 2.0. The region of high T R0 values (i.e.,T R0 > 1.8) is finally governed by the sole P1 II basin, confirming the outcomes of the relevant bifurcation diagram of Fig. 9a which displays the presence of the sole buckled P1 I /P1 II couple for high values of the thermal dynamics.
As a final comment, it is worth noting that the fractal arrangement of the basins herein detected is not present in the previously analyzed bending forced case, due tȯ the different kind of mechanical buckling (symmetric vs antisymmetric) achievable with the two thermal excitations. When the membrane variable is activated, in fact, the coexistence of several buckled solutions possibly reachable by the trajectories causes the fractalization of the relevant basins, while the antisymmetric behavior of the bending variable strongly reduces the multistability region and the basins are organized in a more compact way.

Conclusions
Thermomechanical coupling of a reduced order model of von Kármán shear indeformable plate in active thermal regime and in the presence of an harmonic transverse mechanical excitation and a constant axial excitation has been addressed, focusing on both modeling aspects and phenomenological features of the transient and steady nonlinear dynamic response.
The first part of the work has been devoted to comparatively investigating some main local and global dynamical aspects of the response, as obtained with different reduced order models ranging from the original, fully coupled, model with one mechanical (mid-plane transverse displacement) and two thermal (membrane and bending) variables, down to the simplest uncoupled mechanical model in which the presence of an active thermal environment is accounted for by merely including its known steady effect in the relevant equation of motion. Both the two-way (from thermal to mechanical, and from mechanical to thermal) and the one-way (from thermal to mechanical) coupling have been considered, along with intermediate (two-d.o.f.) models accounting for only the directly excited thermal equation/variable. This has allowed to evaluate the relative importance of the various terms and variables of coupling, in view of a satisfactory, and more or less comprehensive, description of the system dynamics. Based on the considered (membrane or bending) thermal excitation, the alternative two-d.o.f one-way coupled models obtained by skipping either one (bending or membrane, respectively) of the two thermal equations/variables and by neglecting the mechanical contribution into the remaining (membrane or bending) equation have been identified as the most "economical" ones to be used for a systematic numerical investigation aimed at reliably describing the response in terms of mechanical variable and of the dominant (membrane or bending) thermal one, although caught in its solely steady features.
In the second part of the work, the two alternative models have been used to investigate the plate dynamics under a bending or membrane thermal excitation, each one of them inducing a variably rich scenario of buckled mechanical responses. Attention has been focused on the important role played by global dynamics investigations in unveiling the meaningful effects entailed on the plate steady mechanical response by the variably slow thermal transients taken into account by the (though simplified) thermomechanically coupled model. For each of the considered excitations, these effects have been highlighted via a combined use of properly selected cross-sections of the 4D basin of attraction and of the transient time history of the corresponding thermal variable, by looking at the relevant outcomes against the local bifurcation scenario of the mechanical displacement.