Nonlinear dynamics of a third-order reduced model of thermomechanically coupled plate under different thermal excitations

Thermoelastic analysis of a shear deformable reduced model of laminated plates with von K´arm´an nonlinearities and cubic temperature along the thickness is presented. Parametric investigation of the response is accomplished by means of bifurcation diagrams, phase portraits and planar cross sections of the four-dimensional basins of attraction, in order to describe the local and global dynamical behavior of the model. Diﬀerent scenarios induced by the boundary conditions are highlighted, along with the eﬀect of a different modeling of temperature distribution (constant/ dome-shape) on the system response. The eﬀects of di-verse (boundary/body) thermal sources are also com-paratively explored. In all examined cases, proper consideration of system global dynamics is shown to be es-sential for reliably unveiling transient to steady eﬀects due to thermomechanical coupling.


Introduction
Investigation of the effects of thermomechanical coupling on the response of materials and structures in a nonlinear dynamic environment is a research topic of major interest, also in connection with the increased importance of multiphysics phenomena occurring in a va- Bending heat source inside the plate riety of technological contexts, which range from aerospace engineering, to civil engineering, up to microelectro-mechanics. Such effects are generally addressed according to a one-way (i.e. partial) coupling approach, in which changes in the mechanical response are produced by a thermal environment due to the presence of temperature-dependent coupling terms in the equations of motion, but no thermal effects are assumed to be entailed by the structural motion. In this respect, considering a constant temperature field obtained through the solution of the steady-state heat transfer equation, and using finite element and/or reduced-order models (ROMs) [3,4,5,6,7,8,9], interesting results have been obtained, e.g., for plates, as regards variations of frequency-response curves, bifurcation scenarios, transition to chaos, and mode involvement due to nonlinear and thermal couplings. Based on the reasonable assumption that thermal dynamics evolves over a much slower time-scale than structural dynamics, the oneway coupling approach is commonly used also in the case of time-varying thermal environment, although the need to account for vibration modes capable to adapt to instantaneous changes of temperature distribution has recently been highlighted in the framework of model reduction [10].
However, consideration of two-way (i.e., full) thermomechanical coupling, in which the thermal energy equation is coupled with the governing mechanical equations via the presence of additional mechanical (strainrate or velocity) terms, may be important for catching meaningful effects at both material and structure levels. Two-way coupling in elastic materials subjected to quasi-static thermal-mechanical loading was discussed already in earlier studies (e.g. [11]), but the inclusion of inelastic deformation and possible dynamic effects caused by impact loading has recently been seen to entail more meaningful coupling in both monolithic [12,13] and composite [14,15] materials, as well as in composite structures addressed via a multiscale methodology combining global structural analysis with cells micromechanics theory [16]. As regards composite beams and plates, several recent investigations have highlighted the role of full coupling in both static and free vibration analyses (e.g., [17,18,19,20]).
In nonlinear structural dynamics within a steady thermal environment, two-way coupled reduced models of plates employing single-term Galerkin expansions for the transverse displacement and two thermally-induced stress resultants have been used to describe thermoelastic damping of the response, bifurcation behavior, chaotic dynamics, and the influence of various material parameters on the decay of nonlinear vibration amplitudes [21,22,23,24,25]. Recently, a two-way coupled ROM with three mechanical and two thermal variables has been used to investigate the role played by mechanical and thermomechanical impedances in forecasting bifurcation events in a beam, with also experimental validation, as well as the influence of slowly changing thermal loads on the beam post-buckled and snapthrough vibrations [26]. A general theory to study nonlinear dynamics of beam-plate structures for MEMS devices, taking into account the coupling of deformation and temperature fields, has also been proposed, however with the coupling being neglected in the particular case of quasi-static formulation considered as application [27]. Overall, getting rid of the complicatedness generally occurring in the analysis and interpretation of nonlinear phenomena when using richer (e.g., finite element) models, low-order models preserving the main features of the underlying continuum formulations may allow easier analyses and deeper understanding of the basic, yet possibly involved, effects of coupling on the finite amplitude vibrations of geometrically nonlinear structures.
Still in the two-way thermomechanical coupling perspective, a minimal reduced model of shear-indeformable composite plate with von Kármán nonlinearities and assumed linear temperature distribution along the thickness, with one mechanical and two thermal generalized variables, labelled CTC [1], has been used for a systematic investigation of the general features of local and global response, under either passive [28] or active [29,30] thermal conditions. In the latter case, depending on the kind of considered thermal excitation (of either membrane or bending nature), local bifurcation analysis of thermomechanically coupled response has highlighted a quite variable and involved transition from monostable to multistable dynamics, with the latter exhibiting a variety of buckled solutions [31] to be possibly avoided or induced via thermal excitation. In turn, global dynamics has shown to be of major importance for reliably describing the meaningful effects on the structure overall response entailed just by the much slower time scale over which thermal phenomena develop with respect to mechanical ones.
The present numerical study is in-line with the investigation performed in [29], however considering a minimal model which is much more refined than the CTC one used therein, in both mechanical and thermal aspects, while being equivalent in computational terms. Indeed, the new reduced model accounts for third-order shear deformability and a consistently assumed cubic variation of the temperature along the thickness, which still allow to end up to one mechanical and two thermal generalized variables and three corresponding ordinary differential equations (ODEs) [2]. The model is named TTC just for being based on a Third-order theory with Thermomechanical Coupling, instead of the Classical theory with Thermomechanical Coupling in the background of CTC. It has been derived from the twodimensional (2D) continuum formulation of the thermomechanical problem of laminated plates, developed in [2] within the framework of Tonti's unified approach to physical theories [32], via a proper dimension reduction accomplished through kinematic condensation of both in-plane displacements and shear angles [2,33], and single-mode modal reductions of the transverse displacement and the two independent membrane and bending temperatures.
Besides the greater mechanical richness of the somehow 'standard' assumption of third-order shear deformability, the meaningfully major thermal richness of the new reduced model stands in the possibility to account for a much wider variety of thermal boundary conditions, to be possibly prescribed on the plate external (upper and lower) surfaces, than with the previous CTC model. Indeed, besides the sole condition of free heat exchange with the environment, mathematically allowed by the latter, boundary conditions on the external surfaces of the TTC model may include temperature prescribed, heat flow prescribed, and thermal insulation, along with all relevant combinations. Preliminary investigations of the nonlinear response under different boundary conditions have also included the case of a time-varying (impulsive) heat flow [2].
Here, the TTC model is used for an extended numerical investigation of the phenomenological effects produced by a variety of steady thermal sources, of boundary and body type, on the nonlinear dynamic response of a thermomechanical plate. On one side, systematic analysis of local dynamics through bifurcation diagrams with a varying thermal parameter and ensuing localized phase portraits is accomplished. Depending on the considered thermal excitation, it will allow to grasp the overall scenarios of mechanical and thermal responses, with a variable occurrence of multistable periodic solutions, thermally-induced buckling, and nonregular responses. Features of the TTC model in accounting for the variable distribution of temperature through the thickness entailed by different thermal excitations will be addressed, too. On the other side, global dynamics of the two-way coupled model will be investigated via properly selected 2D cross-sections of the actual 4D basins of attraction in the mechanical phaseplane, since one more objective of the paper is to detect the influence of thermal excitations on the mechanical response. Systematic comparison of 2D cross-sections with the actual 2D basins of attraction of the uncoupled mechanical model, in which thermal contributions are taken into account through the mean steady values of the temperature actual evolution, will allow to highlight the pervasive and systematic effects produced on the steady mechanical response just by the slow transient evolution of the temperature, when the structure is subjected to an active thermal excitation in addition to a mechanical one.
The paper is organized as follows. The main features of the TTC model are summarized in Sect. 2, where the variable structure of the three equations governing the thermomechanical problem under different thermal excitations is also presented. The following three sections are the core of the paper, and deal with a com-parative analysis of the effects of the different excitations on the nonlinear dynamic response. In particular, Section 3 considers different thermal boundary conditions, distinguishing between free heat exchange between plate and environment, and constant prescribed temperatures on the plate external surfaces. Section 4 presents the meaningfully different outcomes of two alternative assumptions for the prescribed temperature on the surfaces, which also highlight the capability of the reduced model to account for the effects of a different modelling of the same kind of boundary excitation. Further dwelling with the reduced model capability to account for variable thermal sources, Section 5 investigates the dynamic responses provided by two excitations, of boundary and body type, which, though being thermally equivalent, produce non-trivially different local and global dynamics. A conclusion section ends the paper.

Thermomechanical models
The thermomechanical plate model here used is derived within a unified two dimensional/zero-dimensional (2D/ 0D) modeling framework that integrates the mechanical and thermal aspects, by starting from the three dimensional (3D) physics problems. Such framework is based on the Tonti's unified approach to physical theories [32] and is presented in [1,28,2], which the reader should refer to for the formulation details. Here, a brief summary of the main modelling features is provided, together with a concise description of the formulation steps that allow us to define the sets of the reduced 0D equations used for the analyses.
At the 3D level, the thermoelasticity coupled problem [34] and the nonlinear Green's strain tensor simplified under the assumption of small strains and moderate rotations (von Kármán strains) [35] are considered.
At the 2D level, the basic assumptions consist of considering third order displacements with von Kármán nonlinearities [35], along with a correspondingly consistent third order variation of the temperature along the thickness. By properly combining the configuration, phenomenological and balance multiphysics relations, seven (five mechanical and two thermal) partial differential equations (PDEs) of motion are obtained in terms of seven unknown 2D configuration variables (mid-plane displacements, u, v, w; transverse shear rotations, φ 1 , φ 2 ; membrane and bending temperatures, T 0 , T 1 ), together with the corresponding boundary conditions [2]. In particular, the 3D temperature T (x, y, z, t) is related to the 2D components T 0 (x, y, t), T 1 (x, y, t) as follows: f a (z) = r 1 + r 2 z + r 3 z 2 + r 4 z 3 f b (z) = r 5 + r 6 z + r 7 z 2 + r 8 z 3 f c (z) = r 9 + r 10 z + r 11 z 2 + r 12 z 3 The r i coefficients depend on the kind of imposed thermal boundary conditions, and are reported in Appendix A for the three cases considered in the following.
At the 0D level, in the context of a minimal Galerkin discretization to be pursued under conditions of no internal resonance between the transverse modes of the laminate, single-mode dome-shape approximations are assumed for the transverse displacement w, the transverse shear rotations, φ 1 , φ 2 , and also the temperatures T 0 , T 1 , if considering a condition of cold edges. Looking at the latter, the 0D thermal membrane T R0 and bending T R1 variables are introduced: The in-plane displacement components u, v are expressed in terms of w, T 0 , T 1 through static condensation of the corresponding (in-plane) PDEs of motion (under the usual assumption of a frequency of the in-plane vibration much higher than the frequency of the transverse vibration). After the discretization, a further static condensation of the shear rotations reduces the unknowns to one mechanical (W ) and two thermal (T R0 , T R1 ) variables. The formulation ends up as three coupled nonlinear, nondimensional ordinary differential equations (ODEs). Thanks to the richness and flexibility of the underlying continuum formulation, the model allows to account for a variety of thermomechanical assumptions, excitations and boundary conditions, thus representing a substantial improvement of the thermomechanical model previously investigated by the authors in the nonlinear dynamics regime [28,29,30], characterized by shear indeformability and linear temperature variation along the thickness (CTC model). As regards the mechanical aspect, in the present investigations we consider a simply supported plate, with movable edges subjected to uniform stretching forces of magnitudes p x = p y = p in x and y direction, and transversal resonant harmonic forcing (Fig. 1). As regards the thermal aspect, the plate has isothermal edges and is subjected to the boundary conditions on the up/down outer surfaces and thermal loads as specified below.
The first thermal boundary condition considered consists of a free heat exchange on the up/down outer surfaces (FHE). A temperature difference between plate and surrounding medium or an internal heat source act as thermal loads. Referring to Eqs.(70) of [2], nondimensionalization with respect to time and plate thickness allows to obtain the following set of coupled governing equations (mechanical, thermal membrane and thermal bending, respectively): in terms of the unknown 0D configuration nondimensional reduced variables W (deflection of the center of the plate), T R0 (membrane temperature), T R1 (bending temperature). T ∞ is the time-constant difference between the absolute temperature of the surrounding medium and the reference temperature, while e 0 and e 1 are the membrane and bending heat sources inside the plate. The second thermal boundary condition considered is associated with the condition of prescribed temperature on the external, upper and lower, surfaces of the plate, with the temperature distribution which is assumed to be constant on each of the surfaces (PPC). Referring to Eqs.(72) of [2], nondimensionalization allows to obtain the following set of governing equations: where T up and T down are the prescribed constant temperatures on the upper and lower external surfaces. Comparing Eqs. (4) with Eqs. (3), the different physical process activated by the two types of thermal boundary condition causes different coefficients values, the replacement of the T ∞ term in the membrane equation (4b) by a substantially equivalent thermal term and the addition, into the bending temperature equation (4c), of a new term related to the difference between the upper and lower temperatures on the plate surfaces. As a consequence, when the two faces have different temperatures, both the membrane and the bending thermal variables are activated, differently from the free heat exchange condition which, considering a unique temperature for the external environment, determines the triggering of the sole membrane thermal variable.
About the previous thermal boundary condition, a dome-shape (instead of constant) prescribed temperature on the external surfaces of the plate is also considered (PPD). In this case Eqs. (4)(2) become: where T up and T down are the central values of the domeshape temperature (see Eq. (2)) prescribed on the upper and lower external surfaces. Comparing Eqs. (5) with Eqs. (4), the presence of two new terms related to the T up and T down parameters into the mechanical equation is pointed out, modifying the linear mechanical stiffness and adding a constant external excitation. The coefficients of the PPD model are equal to those of the PPC model, with exception of those related to T up and T down . The dynamical behavior of the thermomechanical model is investigated by considering an epoxy/carbon fiber composite plate of dimensions a = b = 1 m and h = 0.01 m. The material's elastic and thermal properties, which are assumed to be independent of the temperature, are taken from [18], and read: Y 2 = 6.91 · 10 9 N m 2 , G 12 = 3.45 · 10 9 N m 2 , where Y 1 , Y 2 , G 12 are longitudinal modulus of rigidity in x and y direction and shear modulus, respectively; ν 12 is the Poisson's ratio; ρ and δ are mass density and damping coefficient; λ 11 , λ 22 , λ 33 are the thermal conductivities along the x, y and z directions; α 1 , α 2 are the thermal expansions along x and y directions; c v is the specific heat at constant strain, and H is the boundary conductance. The subsequent nondimensional numerical coefficients of Eqs. (3)-(5), used for the numerical analyses, are reported in Appendix A. The mechanical loads applied to the plate consist of an in-plane prestressing force p = 2.51, corresponding to a incipient buckling condition [29], and a transversal weak (f = 1) resonant harmonic forcing, contained in the a 13 and a 17 coefficients, respectively.

Free heat exchange between plate and environment
The first thermal boundary condition considered for the analyses models a free heat exchange (FHE) between plate and environment, which depends on the T ∞ parameter representing a time-constant temperature difference between plate and surrounding medium. Such condition activates pure convection on the external surfaces, and pure internal conduction. The relevant set of equations describing the motion of the plate is reported in the previous section, Eqs. (3). Looking at Eq. (3b), it can be observed that a non-vanishing T ∞ value directly activates the dynamics of the membrane temperature T R0 , due to the naturally symmetric distribution of the temperature along the thickness produced by the thermal condition. The mechanical response is in turn altered thanks to the presence of the coupling linear term a 16 T R0 W which modifies the mechanical stiffness. The dynamics of the thermomechanical plate with free heat exchange condition has been analyzed by the authors in other papers [28,29], in passive as well in active thermal regime, by referring to a classical shear-indeformable von Kármán plate with prescribed linear temperature variation along its thickness. Despite the different mechanical and thermal assumptions considered in the present work, the equations of motion for the two models are formally the same, notwithstanding the different expressions of the a ij coefficients. However, for a thin plate here considered, the two models prove to furnish essentially the same dynamical response, in terms of both local and global dynamics behavior. Thus, this section summarizes the main dynamical features high- lighted in the previous works [28,29,30], to refer to for all details.
The effect of the thermal boundary condition on the mechanical response is reported in Fig. 2(a), in terms of bifurcation diagram for varying T ∞ parameter. Here, the minimum and maximum response curves are reported to point out the evolution of the solution amplitudes. Starting from a monostable pre-buckling configuration, a cooling of the environment (i.e., T ∞ < 0) causes the preservation of monostability, with the sole 1-period existing response P1 which reduces its amplitude for decreasing temperatures. Conversely, a warming of the environment (i.e., T ∞ > 0) produces the arise of two couples of buckled solutions, oscillating around the varied positive/negative configurations of the plate. Their characterization is reported in Fig. 2(b), at the sample value T ∞ = 190. The mechanical phase plane (left panel of Fig. 2(b)) shows the presence of a couple of low-amplitude buckled responses (orange P1 I /red P1 II ) and a couple of high-amplitude buckled solutions (cyan P1 III /blue P1 IV ), coexisting with the gray P1 prebuckling solution, which in the post-buckling regime becomes a cross-well response oscillating around both the varied equilibria. It is worth noting that the system mechanical response is globally symmetric with respect to the trivial equilibrium, while the phase portraits of the solutions of each couple are mirrored with a radial symmetry of the Poincaré points. The membrane thermal response (center panel of Fig. 2(b)), which is the one directly activated by the thermal condition, shows a self-symmetry of the couples of periodic responses, which oscillate, with very low amplitudes, around the T R0 mean steady state value T R0 = 1.7185. Instead, the bending thermal variable (right panel of Fig. 2(b)) oscillates around the rest position, as it is indirectly activated by the mechanical response by means of the coupling terms. For this reason, it is also possible to recognize the same features of the mechanical phase portraits.
As concerns the global dynamics, the evolution of the system basins of attraction has been investigated in [29] by realizing several planar sections of the 4dimensional basins. As the main objective of the analysis is to grasp the influence of the thermal conditions on the mechanical response of the plate, the global dynamics is mostly studied in the mechanical phase plane, by setting trivial thermal initial conditions, i.e., T R0 = T R1 = 0 at t = 0. This is, in fact, the most natural configuration to conceive, from a physical point of view. Figure 3 shows the basin organization in the mechanical plane, at T ∞ = 190, of the model (3) (Fig. 3(a)) and of the relevant uncoupled mechanical model (Fig. 3(b)), obtained by considering only the mechanical equation (3a) in which the thermal contributions are introduced in terms of mean steady state values of the thermal variables, T R0 = 1.718 and T R1 = 0, derived from the resolution of the (uncoupled) thermal equations (3b) and (3c). Comparing the outcomes furnished by the two models, meaningful differences can be pointed out; in fact, the uncoupled mechanical model exhibits a multistable behavior coherent with the results furnished by the bifurcation diagram of Fig. 2(a), while the coupled system displays an (almost) monostable behavior, with the presence of the sole pre-buckling solution (gray) basin, suggesting the inability of the thermal boundary condition to induce mechanical buckling. As exten-sively discussed in [29,30], this discrepancy is due to the long transient needed by the thermal variables to attain their steady values, which provides a slow contribution to the mechanical stiffness responsible for the buckling occurrence. As a consequence, the mechanical response, which has a much quicker dynamics, settles at its first steps of the temporal evolution to the sole stable solution depicted by the system in the pre-buckling regime, i.e. the P1 gray solution. Since this solution represents a robust attractor also in the post-buckling scenario, the trajectories settled on it do not modify their behavior even when the thermal evolution is completed and all the contribution has been furnished to the mechanical stiffness. Indeed, the apparent incongruity between the results of local and global dynamics for the coupled model can be easily explained by remembering that Fig. 3(a) represents only a planar section of the 4D basins, obtained by properly fixing the thermal initial conditions. Changing them allows us to describe the actual basin spatial organization, up to exactly reproducing the response of the uncoupled model when they are set to the relevant mean steady values, as shown in Fig. 3(c). In this case, in fact, the thermal evolution is precluded and the coupled system reduces to the uncoupled one, which demonstrates to describe only a specific condition of the thermomechanical plate. this kind of thermal condition affects also the bending equation, due to the fact that the temperature can have a generic, also nonsymmetric, distribution along the plate thickness, depending on the temperature values imposed on the upper (T up ) and lower (T down ) surfaces. Moreover, the pure conduction mechanism in the PPC model causes an increase in the numerical values of all coefficients relevant to the thermal equations (4b) and (4c), with respect to the FHE model accounting also for convection. As a consequence, the PPC model under pure mechanical excitations displays the same mechanical response as the FHE system, with however higher amplitudes of the thermal responses activated by the mechanical coupling terms a 24 and a 33 . Looking at the effect of the thermal condition on the system dynamics, the bifurcation diagrams as a function of the temperature on the lower surface T down are shown in Fig. 4. Here, the temperature on the upper surface has been fixed to T up = 100, a value that implies the arise of the two high-amplitude buckled responses P1 III and In addition to the loss of symmetry distinguishing the PPC response with respect to the FHE one, another interesting difference in the response of the two systems can be observed in the monostable region obtained for negative T ∞ (T down ) values in the FHE (PPC) case. Figure 4(a) shows the presence of a wide unstable region for −633 < T down < −356, characterized by the chaotic response shown in Fig. 5, for which the first Lyapunov exponent is +0.05.
In order to analyze the global dynamics of the PPC model, the planar section of the basins of attraction for the coupled system with T up = 100, T down = 320 is reported in Fig. 6(a), while the global response of the relevant uncoupled mechanical system with steady state values of the thermal variables T R0 = 1.93431, T R1 = 2.03161 (marked by the blue lines in Figs. 4(b) and (c)) is displayed in Fig. 6(b). Looking at the latter, the symmetry breaking can be recognized also in terms of basin organization inside the two buckled wells, while the comparison between the responses of the coupled system with thermal initial conditions equal to zero ( Fig. 6(a)) or to the relevant steady state values ( Fig. 6(c)) confirms the crucial role played by the long thermal transient in modifying the stationary response scenario in the mechanical cross-section (see also the mechanical phase portrait in Fig. 6, where the thermal responses are reported for completeness, too).
When a symmetric temperature condition is applied to the external surfaces of the PPC model, the sole membrane thermal equation is activated, thus producing a system response which is similar to the one studied with the FHE model when a temperature difference is imposed between plate and environment. Figure 7 displays the comparison between the response of the coupled (a,b) and uncoupled (c,d) FHE (a,c) and PPC (b,d) models when T ∞ = 100 and T up = T down = 100, respectively. The systems show a similar behavior in terms of both the coupled and the uncoupled dynamics, despite some minor differences, as e.g. a slightly more marked, and strongly fractalized, presence of the buckled basins in the coupled PPC model (Fig. 7(b)), and a reduction of their compact part in the relevant uncoupled response (Fig. 7(d)). As already observed, they are due to the different physical processes involved in the two thermal conditions which reflect on the values of the system coefficients. In particular, the higher value of the membrane temperature stiffness coefficient a 22 in the PPC model is responsible for the sensible shortening of the relevant thermal transient, as reported in the central panel of Fig. 8. However, even if quicker than in the FHE case, the thermal evolution of the PPC model is orders of magnitude slower that the mechanical one, so that it is unable to move the mechanical response to a buckled configuration.

Effect of a different temperature function on the outer surfaces
In this section, attention is devoted to discuss the effects of choosing different temperature functions when imposing the thermal boundary conditions. With reference to the plate with prescribed temperature on the external surfaces, the response of the PPC model, in which the temperature on the external surfaces is assumed to be constant, is compared with the behavior of the model with dome-shaped temperature (PPD -Prescribed Prescribed Dome-shape). The latter assumption, although less intuitive from a physical point of view, results to be coherent with the modal shapes assumed for displacement and temperature variables in the formulation of the reduced order model. Bifurcation diagrams for the PPC (black) and PPD (red) models are reported as a function of the temperature on the upper ( Fig. 9(a), with T down = 0) and lower ( Fig. 9(b), with T up = 100) surfaces. The outcomes point out differences in amplitude and stability regions of the buckled dynamical response, which are mostly relevant to one of the two buckled wells, corresponding to the mechanical configuration bent towards the colder surface. In fact, positive values of the T up parameter in Fig. 9(a) imply a warmer upper surface (i.e. the surface at W < 0, see Fig. 1), while for T down > 100 in Fig. 9(b) the upper surface becomes the colder one.
Looking at the pre-buckling scenario, bifurcation diagrams in Fig. 9(b) highlight another important difference between the two models. In fact, for negative T down values, the unstable region characterized by chaotic response depicted by the PPC model and described in Fig. 5 is replaced, for the PPD model, by the negative high-amplitude buckled response P1 III , reported in Fig. 10, which remains stable in the whole T downnegative range. This entails a monostable periodic scenario lasting for a wide range of T down values.
As regards the behavior of the thermal variables, the same linear trend can be highlighted for both models, with the PPD system reaching lower temperatures with   Fig. 9 Comparison between PPC (black) and PPD (red) models: bifurcation diagrams as a function of the mechanical displacement W for varying temperature on the upper surface Tup, with T down = 0 (a); on the lower surface T down , with Tup = 100 (b); bifurcation diagrams as a function of the membrane temperature T R0 (c) and of the bending temperature T R1 (d) for varying temperature on the lower surface T down , with Tup = 100. respect to the PPC ones. This is due to the numerical values attained by the coefficients a 23 and a 35 relevant to the parameters T up , T down , which are lower for the PPD model with dome-shaped prescribed temperature (see Eqs. (9),(10)), thus reducing the effect of the thermal boundary condition. It is worth noting that all the other coefficients are equal for both models, due to the same mechanical and thermal assumptions. The global dynamics of the coupled (a) and uncoupled (b) PPD model with symmetric thermal boundary condition T up = T down = 100 is presented in Fig. 11, to be compared with the relevant results obtained for the PPC model and presented in the previous Fig. 7(b),(d). Due to the different steady state mean value of the T R0 variable reached by the two models, observable in Fig. 9(c), the response of the uncoupled mechanical models is determined by setting T R0 = 0.920199 and T R0 = 0.567271 for the PPC and for the PPD model, respectively. However, the relevant outcomes (Fig. 7(d) vs Fig. 11(b)) are very similar in terms of number of basins identified and their organization and compactness.
Conversely, the behavior of the basins of attraction of the PPD coupled model (Fig. 11(a)) is characterized by a stronger presence of the fractalized buckled basins with respect to the PPC case ( Fig. 7(b)). Such result might appear counterintuitive, given that in the PPC model the effect of the thermal boundary condition is strengthened by higher values of the relevant coefficients, suggesting a stronger contribution of the thermal dynamics to the mechanical buckling. The explanation has to be sought in the occurrence of the term related to a 18 , present in the PPD mechanical equation (5a) and absent in the PPC one (4a), through which the thermal condition directly modifies the mechanical stiffness responsible for the buckling occurrence.
A different scenario reported in Fig. 12 is obtained if applying an antisymmetric thermal condition T up = 100, T down = −100, which activates the sole bending temperature T R1 . Also in this case, the different temperature function implies a different steady state value of the bending thermal variable for the two models, corresponding to T R1 = −1.84692 for the PPC system, and to T R1 = −1.13818 for the PPD one (see Fig. 9(d) and Fig. 13(a)). Moreover, the different temperature func-tion prescription on the boundary entails qualitatively different dynamical scenarios, corresponding to monostability (PPD) and bistability (PPC), which reflect on evident differences in both the coupled and the uncoupled responses. The monostable P1 III buckled response of the PPD model obviously implies a coincidence of the coupled and uncoupled outcomes (Fig. 12(b) and (d), respectively). As concerns the PPC model, the uncoupled mechanical system (Fig. 12(c)) displays a bistable dynamics with a dominance of the P1 III buckled basin, as well, while in the coupled model ( Fig. 12(a)) the presence of the thermal transient slows down the thermal contribution to the mechanical buckling and moves the system response to a mainly pre-buckling (gray) configuration. As a result, the PPC and PPD coupled models furnish different responses for most of the mechanical initial conditions, as exemplarily shown in Fig. 13.
To clearly understand the reason behind the lower values attained by the thermal variables of the PPD model with respect to the PPC one, it is useful to investigate the temperature evolution along the plate thickness, also separately analyzing the various contributions relevant to the different powers of z described in Eq. (1). In Fig. 14 reference is made to the previously discussed cases of symmetric ((a), (b)) and antisymmetric ((c), (d)) thermal boundary conditions, by looking at the central point (i.e. x = a/2, y = b/2) of the plate. As expected, the symmetric condition activates only the symmetric constant and quadratic contributions, while the antisymmetric one involves the sole linear and cubic antisymmetric contributions. Moreover, the prescribed temperature on the upper and lower surfaces is guaranteed by the fixed values (100/100 or 100/-100) reached at z = −h/2 and z = h/2, respectively. However, looking at the temperature behavior of the PPC model in the middle of the plate (at z = 0), the curve shift with respect to zero (Fig. 14(a,b)) and its angular coefficient (Fig. 14(c,d)), corresponding to T R0 and T R1 in the reduced model respectively, reach evidently higher values with respect to the PPD case. This is due to the discrepancy between the mathematical functions chosen for modeling the thermal variables (dome-shape functions) and for the prescribed boundary condition (constant function) in the PPC model, which enforces temperature to achieve high values in the center of the plate in order to guarantee the prescribed temperature on the external surfaces and concurrently assume an internal dome-shape distribution. Conversely, the PPD model assumes the same dome-shape function for both variables and boundary conditions so that the temperature along the plate thickness can actually develop in an almost purely constant or linear way (i.e., without activating nonlinear contributions), coherently with the symmetric or antisymmetric boundary conditions imposed on the external surfaces. In fact, for the latter model, the results show a negligible contribution furnished by the quadratic and cubic terms, which instead become important to describe the overall temperature behavior of the PPC model. different behavior of the coupled (a) and the uncoupled systems (b), the latter being only resolved when considering in the coupled model the proper steady thermal i.c. (Fig. 15((c)). Moreover, the results underline that the choice of the temperature function is able to qualitatively and quantitatively modify the global dynamics, as well. Looking at the uncoupled systems, major differences in the basin organization can be pointed out inside the well towards the upper colder surface (W < 0), as already observed, in which the basin of the orange P1 I solution disappears in the PPD model and is replaced by that of the P1 III cyan response. Moving to the coupled systems, the PPC model displays the presence of the dominating P1 gray basin coexisting with the fractalized P1 IV blue one, while for the PPD model the P1 IV basin becomes the prevalent response.
Diagrams reported in Fig. 16 confirm the major role played by the constant and linear contributions in determining the steady temperature along the plate thickness for the PPD model ( Fig. 16(a)), with respect to the PPC case in which the nonlinear quadratic and cubic terms actively participate to define the temperature pattern ( Fig. 16(b)).

Effect of an equivalent internal thermal source
The effect of the previously studied thermal boundary excitation is here compared with that produced by a bending thermal body source, which in the 3D model represents a power density possibly due, for example, to the passage of an electric current through the plate, chemical reactions or nuclear fission. In particular, reference is made to the thermomechanical model with prescribed dome-shaped thermal boundary condition described by Eq. (5), in which purely antisymmetric temperatures on the lower and upper surfaces are considered, i.e. T up = −T down . As a consequence, the thermal dynamics develops in terms of the sole bending thermal variable T R1 . In order to determine the bending thermal excitation e 1 (t) able to produce the same effect on the thermal dynamics, the steady uncoupled version of Eq. (5c) is solved to obtain the following boundary condition-body source relation: e 1 = a35 a34 α 1 (T up − T down ).
The evolution of the stationary solutions when alternatively applying T up /T down and e 1 is presented in ics to develop for lower values of the applied thermal source. Comparing blue and red curves in Fig. 17, it can be observed that the response (in terms of existence and stability) furnished by the e 1 excitation (blue) is somehow stretched out with respect to that of the T up /T down boundary condition (red), and the evolution of the various periodic solutions calls for higher thermal energy.
This behavior is confirmed by the analysis of the global dynamics in terms of basins of attraction in the mechanical plane. The response of the PPD model under antisymmetric boundary condition T up = 100 = −T down already presented in Fig. 12(b),(d) can be compared with that produced by an equivalent body source e 1 = −9.83 · 10 −4 , reported in Fig. 18, for the coupled and the uncoupled systems. The outcomes show that the application of the boundary condition causes the instabilization of the P1 solution and the occurrence of a monostable scenario characterized by the buckled P1 III response ( Fig. 12(b),(d)). Differently, when applying the equivalent e 1 excitation the P1 response remains stable and, indeed, appears as the dominant solution for both the coupled and the uncoupled models (Fig. 18), even though with quantitative differences in basin size due to the presence of thermal transient in the coupled case.
To complete the comparison, it is also of interest to analyze the temperature evolution along the plate thickness for the two cases, as reported in Fig. 19. The prescribed temperature on the external surfaces T up = 100 = −T down produces a linear distribution of temperature, already presented in Fig. 14(d), due to the coher-ence of mathematical functions chosen for the boundary conditions and for the temperature variable, which characterizes the PPD model. A significantly different behavior is instead obtained when applying an equivalent body source e 1 , as shown in Fig. 19(a). Note that the equivalence is guaranteed by the same value of the reduced bending thermal variable T R1 attained by the two models. However, the application of the sole e 1 excitation implicitly requires the temperatures on the upper and lower surfaces T up /T down to be fixed to zero value, so that they do not produce contributions into the equations. It is worth remembering that in the free heat exchange (FHE) case the T ∞ parameter representing the thermal boundary condition corresponds to a temperature difference between plate and environment, and setting it to zero implies the thermal equilibrium with the surrounding medium. Conversely, T up and T down parameters describe the actual temperature (in Kelvin) of the plate external surfaces, and fixing them to zero causes the temperature along the thickness to assume zero value at z = ±h/2, thus settling to a cubic evolution.

Conclusions
The behavior of a reduced model of third order shear deformable single-layer orthotropic plate with thermomechanical coupling is investigated in terms of local and global dynamics. Differently from the shear indeformable model with linear temperature distribution along the thickness analyzed in [29,28], the temperature evolution described with the cubic function allows the present model to account for a variety of different thermal boundary conditions, as well as to critically discuss the effect of the assumed temperature functions on the overall response of the model.
As concerns the boundary condition modeling, the dynamics of a plate under free heat exchange with the environment is compared with the condition of prescribed temperature on the plate outer surfaces, in order to stress the effect of thermal constraints on the mechanical response of the system. The different physical phenomena activated by the two conditions have shown to reflect not only on the values of the equation coefficients, but also, and mainly, on the structure of equations themselves, leading to a qualitatively different dynamical behavior. While the free heat exchange condition directly activates only the membrane temperature equation, inducing a symmetric dynamics around the two buckled mechanical configurations, the assignment of a prescribed temperature on the external surfaces of the plate involves also the (antisymmetric) bending thermal variable which is responsible for an overall nonsymmetric mechanical behavior of the model.
Moreover, the choice of the mathematical function describing the temperature distribution on the boundaries is seen to be decisive in determining the mechanical response of the plate. In fact, considering a domeshape distribution, coherent with the mode-shape assumption for the mechanical and thermal variables, implies the addition of a boundary term also into the mechanical equation, and the comparison with the dynamical behavior of a plate with constant prescribed temperature points out differences in terms of existence, stability and amplitude of the possible stationary solutions depicted by the system. In general, the outcomes show that when assuming the same function for inter-nal variables and external conditions, the dynamics of the reduced model settles to lower temperature values, and the temperature distribution along the thickness does not activate the higher quadratic and cubic contributions (although provided by the model), arranging itself in a linear way.
Finally, the effect of different thermal sources is discussed, by comparing the thermal boundary conditions with a body source producing an equivalent effect of the thermal dynamics. Also in this case, the different mechanical scenario furnished by the two cases points out that the choice of a proper model for the thermal excitations is able to strongly influence the response of the system, not only in terms of thermal evolution, as expected, but also as regards the mechanical behavior.
As a concluding comment, the presented results confirm the crucial role played by the thermal transient in determining the steady state behavior of the thermomechanical plate, whose mechanical response demonstrates to be strongly different from that obtained by uncoupled models which intrinsically neglect the thermal evolution, thus stressing the importance of accounting for the thermomechanical coupling to grasp the actual behavior of multiphysics models characterized by slow and fast dynamics.

Acknowledgement
The financial support of Italian research program PRIN 2015 (no. 2015JW9NJT) is gratefully acknowledged.

Compliance with ethical standards
The authors declare that they have no conflict of interest.