Unified 2D continuous and reduced order modeling of thermomechanically coupled laminated plate for nonlinear vibrations

A unified formulation of the thermomechanical problem of laminated plates with von Karman nonlinearities, undergoing finite amplitude vibrations, is presented. It integrates mechanical and thermal aspects, by addressing them in parallel via the introduction of generalized 2D variables and equations also for the latter. The formulation virtually embeds a multitude of possible models, resulting from different assumptions about the plate mechanical and thermal configurations. The obtained continuous model is then subjected to a minimum reduction via Galerkin procedure. Some analyses of free and forced nonlinear vibrations under variable mechanical and/or thermal excitations are also carried out, to get some hints on the importance of different thermal aspects associated with membrane and bending dynamics, and on the possibility to catch them via variably simplified models.


Introduction
In contemporary literature on nonlinear structural mechanics, multiphysics analyses tailored to reliably account for the involved coupling effects which ensue from the structural system being embedded in a multifield environment are cutting edge topics. The fundamental interest is still concerned with the evaluation of typical mechanical quantities (displacements, vibration, stresses, etc.). Yet, the way in which they are affected by the variable level of interaction with characteristic quantities of coexisting physical fields, of either mechanical (fluids, other solids) or non-mechanical (thermodynamics, magnetism, piezoelectricity, atomic interaction, etc.) nature, is of major importance to the aim of a reliable analysis and design of many modern engineering systems.
Compared to the great interest in these analyses, the non-mechanical aspects are often poorly integrated within the formulation of the 2D structural mechanics theories (beams, plates, shells, etc.) which are now consolidated in engineering research and practice. For example, in the context of dynamic problems of thermoelasticity [1], classical 2D plate structural theories (purely mechanical) [2][3][4][5] are usually flanked by occasional and non-systematic treatments [6][7][8] of thermoelastic coupling, by merely adapting the 3D relationships of heat transfer and thermoelasticity to the problem at hand.
The main objective of the present work is to address such a lack of formulation and to propose a 2D approach which integrates mechanical and thermal aspects by identifying generalized 2D variables and governing equations also for the thermal aspects of the problem. The basic idea consists of using the classical Tonti diagram [9] for physical theories as a suitable framework to implement a unified, consistent and controllable 2D formulation of the thermomechanical problem for geometrically nonlinear laminated plates undergoing finite amplitude vibrations.
The proposed formulation virtually embeds a multitude of possible plate models resulting from different assumptions about mechanical and thermal configurations. So, it represents an advantageous tool for pursuing multiphysics analyses of the nonlinear dynamics exploiting (i) its systematic approach, (ii) the awareness about theoretical assumptions, (iii) the easiness in comparing different models, and (iv) the effective context for computer implementation. In this last respect, general and versatile formulations already present in the literature within a similar unifying perspective are aimed at dealing with an even larger variety of structural theories mostly in view of effective finite element implementations [10][11][12].
The laminated plate in a thermoelastic environment is a scientifically and technically interesting system to be addressed in the unified 2D perspective. Thermal effects in nonlinear dynamics of laminated plates have been recently addressed in a number of papers by considering different discretization procedures and solution techniques. The analyses include thermomechanically uncoupled or coupled formulations. In many real world applications, the vibrations of a plate are associated with a change of plate temperature: a deformation varying in time entails temperature changes and, conversely, temperature changes induce strains. Uncoupled thermomechanical vibration analyses of structural elements [13][14][15][16][17][18][19][20][21] are often carried out by neglecting the interaction phenomenon: knowing the distribution of temperature resulting from the solution of the heat conduction equation, the motion equations with known terms involving temperature gradients are solved. In contrast, coupled vibration analyses [6][7][8][22][23][24] take into account the actual thermomechanical interaction entailed by the presence of displacement (temperature) unknowns in the thermal (motion) equations.
Within a generally limited availability of studies accounting for actual coupling, even less is the number of those using low-order models allowing to understand its basic effects on system nonlinear vibrations. Chang and Wan [22], Ĉ ukić [23] and Yeh et al. [6] investigated the dynamic behavior of thermomechanically coupled isotropic elastic plates by employing a single-term Galerkin expansion. They indicated that the main effect of coupling is to introduce a thermoelastic damping. Yeh et al. [6] and Yeh [7] also analyzed the aperiodic motion of a rectangular isotropic and orthotropic plate, respectively, and reported that thermoelastic damping tends to stabilize the motion of the plate. Yeh [8] parametrically analyzed the influence of various material parameters on the nonlinear free vibration of the orthotropic plate.
Within the pursued unified perspective, another goal of the present work is to associate to the underlying partial differential continuous modeling a minimal system reduction still accounting for the main features of geometric nonlinearity and thermomechanical coupling. It ends up with three coupled nonlinear ordinary differential equations extending to the case of symmetric laminates a reduction previously accomplished for the orthotropic plate [7,8]. It is useful to also allow us to describe and compare the plate nonlinear vibrations under different modeling assumptions and/or simplifications, as well as physical conditions. In fact, the minimal model is used to identify variable contributions of thermal effects, to drive the numerical investigation, and to evaluate the importance of thermomechanical coupling, via preliminary analyses of free and forced nonlinear vibrations. The outcomes of the complete reduced model and of simplified models are compared, and the features of some main thermomechanical phenomena are investigated.
The paper is organized as follows. The unified formulation is sketched in Sect. 2. The modeling of the thermoelastically coupled laminated plate with von Karman nonlinearities is presented in Sect. 3. The minimal dimension reduction of the continuous model, and the analytical framework to support the numerical investigation are addressed in Sect. 4. Analysis of thermal aspects in nonlinear vibrations is carried out in Sect. 5. The paper ends with some conclusions.

A unified scheme for 2D modeling of thermoelastic laminated plate
Initially, the thermo-mechanical problem of the plate may be considered a generic mathematical problem of 3D multiphysical field, schematized as follows: Assigned (i) the geometric region of the field, (ii) the nature of the material that fills the region, (iii) the sources and (iv) the boundary conditions of the field ? determine its configuration.
As it is known, the formulation of such a problem leads to the so-called fundamental model of the field, which is constituted by a set of PDEs establishing a direct relationship between source and configuration variables (Fig. 1).
The structure of many physical theories, aimed to build the various fundamental models, is independent of the underlying physical content, as the classical Tonti diagram [9] highlights; it decomposes the 3D fundamental equations of a generic theory into three sets of configuration, balance and phenomenological equations.
In 2D engineering multiphysics modeling, a similar structure can be obtained starting from the basic generalizing assumption 3D ? 2D for the plate: that expresses the 3D configuration variables in terms of 2D generalized ones (through mathematical function shapes), and generalizing all the equations and variables of Tonti's decomposition. Performing this procedure for the 2D modeling of thermoelastic laminated nonlinear plates, the scheme in Fig. 2 Fig. 2) we obtain the fundamental model of the plate, namely the thermoelastic governing equations which directly link fundamental source (body force and source energy) and configuration (displacement and temperature) variables.
Boundary conditions must be associated with the fundamental model, to specify the value of some variables on the plate edges.
Different assumptions about mechanical and thermal configurations (Fig. 2) lead to a multitude of fundamental plate models, which those available in the literature can be referred to.
As regards the mechanical assumptions, the change of configuration can be described by composing two aspects: (i) the reference-plane deformation, which depends on the assumptions about the displacements of the reference plane; (ii) the crosssection deformation, which depends on the assumptions about the displacements of planes parallel to the reference plane.
As regards the thermal assumptions, the mathematical aspect of the equations in Fig. 2 depends on the hypotheses about the features of temperature variations.
In the next section, the unified scheme of Fig Consider the laminated simply supported rectangular plate in Fig. 3. The plate has a thickness h and edge lengths a and b, respectively, in the x-and y-directions. The reference plane of the plate coincides with the xy plane of an orthogonal Cartesian coordinate system. The plate is subjected to both mechanical lateral and thermal loadings. Moreover, the edges of the plate are subjected to uniform stretching forces of magnitudes p x and p y in the xand y-directions, respectively.

Fundamental configuration variables
The various fundamental models for the 2D plate depend on the selected mathematical form of the basic generalizing assumption 3D ? 2D (1), that for the problem of thermoelastic plate can be decomposed into: In Eq. (2a, b) displacement and temperature variables depend only on the x and y coordinates of the reference plane and on time t, while the shapes govern the dependence on the thickness z coordinate.
For the model CTC, Eq. (2a) has the following expression [2,3]: where the functions u 1 (x, y, z, t), u 2 (x, y, z, t) and u 3 (x, y, z, t) are the components of the 3D displacement variable along the x, y and z directions, while u(x, y, t), v(x, y, t) and w(x, y, t) are the components of the displacement variable of the 2D plate model. The latter represent the unknown displacements of the reference plane along the x, y and z directions. Equation (3) results from the following assumptions: the cross section is (i) plane, (ii) normal to the reference plane (shear-indeformability) and (iii) invariable along the z direction. About Eq. (2b), it is assumed that the temperature varies linearly, consistent with the mechanical assumption (3); so we can write: where the function T(x, y, z, t) is the 3D temperature variable, while T 0 (x, y, t) and T 1 (x, y, t) are the unknown membrane and bending components of the temperature variable of the 2D plate model. Thus, the following components for the configuration variables of the CTC model emerge (Fig. 2): temperature: T 0 T 1 f g ð6Þ

Expressive configuration variables and related equations
As regards the mechanical aspect, the expressive configuration variables are: the velocity with components defined as and the deformation with components defined as [2,3] e ð0Þ The latter are related to von Karman 3D strainsassociated with Eq. (3)-by [2,3]: 12 Þ are the bending strains (known as curvatures), that is, the components of the deformation variable.
As regards the thermal aspect, the expressive configuration variable is the thermal gradient, with components defined by [1]: that are related to 3D thermal gradients-associated with Eq. (4)-by: In Eq. (10), ðg 2 Þ are the membrane and bending thermal gradients, respectively.
Thus, the following relationships for the Configuration block of the CTC model emerge (Fig. 2): C1 velocity$displacement À À À À ÀÀ ! Eq. (7) C2 deformation$displacement À À ! Eq. (8) C3 thermal gradient$temperature À ! Eq. (10) that define, in terms of fundamental variables, the following components of the expressive configuration variables: velocity: deformation: e where f 1 , f 2 and f 3 are the components of the 3D body force vector along the x, y and z directions (including in Eq. (17) also the vertical loads applied on the outer surface but assuming zero the moments) and E is the 3D source energy due to chemical reactions, nuclear fission or passage of electric current through the plate [1].

Expressive balance variables and related equations
As regards the mechanical aspect, the expressive balance variables emerge from the known motion equations of the von Karman plate [2,3]: are the membrane and bending components of the tension variable, ensuing from the 3D stresses, while R i,t are the components of the linear momentum variable.
As regards the thermal aspect, the expressive balance variables and the related equations can be obtained by the 3D thermal balance equation [1] in the case of non-stationary conduction and thermoelastic coupling: where q i is the 3D heat flow along the x, y and z directions, b is the 3D internal energy due to the nonstationary conduction, a is the 3D interaction energy due to the thermoelastic coupling and E is the 3D source energy. The 2D balance consists of two equations, the first is obtained from Eq. (21) performing integrations over the thickness h of the plate, and the second is obtained multiplying Eq. (21) by z and then performing integrations over the thickness h; it gives: where the following 2D quantities are defined as In Eq. (23), q Thus, the following relationships for the Balance block of the CTC model emerge (Fig. 2): B1 tension$linear momentum$displacement$body force À À À À À À À À À À ÀÀ ! Eq.
that link the following components of the expressive balance variables linear momentum: heat flow: q to the fundamental source and mechanical configuration variables. These quantities will be explicited in the following subsection.

Phenomenology block
The 3D phenomenological equations allow to express the 3D balance variables in terms of the 3D configuration variables [9]. The 2D generalization of these equations is obtained integrating over the thickness h of the plate; it brings out the global phenomenological coefficients associated with the entire laminate. Operationally, this procedure is characterized by three aspects.
1) Composite laminates have several (N) layers, each with different orientation of its material coordinates with respect to laminate coordinates. Thus, the following relations can be used to transform phenomenological equations from the material coordinates ð" x; " y; " zÞ of each layer to the plate coordinates (x, y, z) used in the problem description [2,3] (Fig. 4): x y z 8 > < > : 2) Although the 3D configuration variables are continuous through the thickness, 3D balance variables are not, due to the change in material coefficients through the thickness (i.e., each lamina). Hence, the integration of 3D balance variables through the laminate thickness requires lamina-wise integration.
3) The global phenomenological coefficients have components of various types (membrane, bending, bending-membrane, higher order) depending on the initial assumptions about the fundamental displacement and temperature variables, and on the lamination scheme.

Equations aimed to explicit the linear momentum
The components of the linear momentum variable are provided by the Newton's second law generalized to the entire laminate; in the case of CTC model [2,3] it gives: where the inertial quantities I i are defined as: with q (k) mass density of the generic lamina.

Equations aimed to explicit the tension
The thermoelastic linear constitutive relations for the kth orthotropic lamina in the principal material coordinates of a lamina are [2]: where " Q ðkÞ ij are the plane stress-reduced elastic stiffnesses, and " b For a lamina arbitrarily oriented with respect to the plate coordinate system, the previous equation becomes:  x; " y; " z ð Þ and plate (x, y, z) coordinate systems [2] that provide where the quantities (A ij , B ij , D ij ) are respectively the membrane, bending-membrane and bending elastic stiffnesses of the laminate, defined as while the quantities ðb A ij ; b B ij ; b D ij Þ are respectively the membrane, bending-membrane and bending thermoelastic stiffnesses of the laminate, defined as

Equations aimed to explicit the heat flow
The Fourier law for the kth orthotropic lamina in the principal material coordinates of a lamina is [1]: For a lamina arbitrarily oriented with respect to the plate coordinate system, the previous equation becomes: where the transformation relations of the thermal conductivities k where the quantities ðk A ij ; k B ij ; k D ij Þ are respectively the membrane, bending-membrane and bending thermal conductivities of the laminate, defined as Remark Recalling Eq. (11), Eq. (40) provides 33 T 1 ; therefore the heat flow component q 3 varies through the thickness due to the change of conductivity k ðkÞ 33 through the thickness (i.e., each lamina), but not due to the temperature, which has constant gradient along z. All that implies q 3,z = 0 in 3D thermal balance (21): it is a situation not realistic but acceptable in many cases, which shows the physical inconsistence of the linear thermal configuration assumption (4). With regards to the 2D thermal balance (22), the energy rates due to the heat flow component q 3 is however accounted by the global quantities Q (0) and Q (1) . They are obtained from the convective boundary condition: heat flow due to conduction must be in balance with the heat flow due to convection (Sect. 3.4). Therefore it is a thermally balanced but not congruent theory. This thermal approximation appears quite analogous to the Kirchhoff mechanical approximation (linear displacement assumption (3)); in fact, with regards to the mechanical aspect, in Eq. (19c) the transverse shearing forces are present only for the purposes of balance (although they are not directly visible for being expressed in terms of moments [2,3]), but they do not appear in the generalized variables of the problem (as the 2D heat flow along z direction). Instead, it is noted that a cubic thermal configuration assumption would describe a more realistic situation, because it involves a certain consistence of thermal configuration, and appears quite analogous to the Reddy mechanical approximation [2].

Equations aimed to explicit the internal energy
The internal energy for the kth lamina is expressed in terms of temperature by [1,7]: where C (k) = q (k) c v (k) is the thermal capacity of a lamina, function of mass density q (k) and of specific heat at constant strain c ðkÞ v . Then the components of the internal energy variable, recalling the (43), (23) and (4), are: where the quantities (C A , C B , C D ) are respectively the membrane, bending-membrane and bending thermal capacities of the laminate, defined as

Equations aimed to explicit the interaction energy
The interaction energy for the kth orthotropic lamina is expressed in terms of strains by the 3D thermoelasticity theory [1,7], but omitting the term b zz e zz due to the usual mechanical hypothesis for the plate, e zz = e 33 = 0, implicitly assumed with Eq. (9); so it is: where the thermoelastic stiffnesses b ij have been already specified at the Eq. (34), Then the components of the interaction energy variable, recalling Eqs. (47), (23) and (9), are: 12 þ z e where the quantities ðb A ij ; b B ij ; b D ij Þ are respectively the membrane, bending-membrane and bending thermoelastic stiffness of the laminate, already defined by Eq. (37).
Thus, the following relationships for the Phenomenology block of the model CTC ensue (Fig. 2): P1 linear momentum$velocity À À À À À À ! Eq. (31) P2 tension$deformation$temperature ÀÀ ! Eq. (35) P3 heat flow$thermal gradient À À À À À À À ! Eq. (41) P4 internal energy$temperature À À À À À À ! Eq. (44) P5 interaction energy$deformation À À ÀÀ ! Eq. (48) quantitatively governed by the following global phenomenological coefficients: inertia 3.4 Energy rates of the out-of-plane heat flow The energy rates Q (0) and Q (1) in Eq. (24), due to the heat flow q 3 in z direction, can be expressed through the convective boundary conditions on the upper and lower surfaces of the plate [1,7]: where the first member is the heat flow in z direction q 3 at the outer surfaces of the plate; the second member is its expression in terms of thermal gradient and conductivities; the third member is the heat flow exchanged by convection between the outer surfaces of the plate and the surrounding medium; the fourth member is its expression by the Newton's law of cooling [1], where H is the boundary conductance, T ? is the constant difference between the reference temperature and the absolute temperature of the surrounding medium, (T) h/2 and (T) -h/2 are the temperature values at the outer surfaces of the plate. Thus, performing the integration lamina by lamina along the thickness of the laminate and using Eq. (49), Eq. (24) becomes (omitting mathematical steps):  (19) and (22)], we obtain the fundamental CTC model of the thermoelastic laminated plate. These equations are expressed in terms of fundamental variables displacement and temperature and are given in Eq. 82 in Appendix. In the particular case of orthotropic plate, they reduce to the analogous equations in [7].

Boundary conditions
The mechanical and thermal boundary conditions associated with the fundamental CTC model are concerned with: properly expressed in terms of displacements.
that for the problem of thermoelastic plate can be decomposed into: Equations (52) and (53) express the 2D generalized configuration variables (remind Eqs. (1) and (2a, b)) in terms of 0D reduced configuration variables (through shape functions); in Eq. (53a, b) the reduced displacement and temperature variables depend only on time t, while the shapes govern the dependence on x and y in the reference plane.
In the general case, the mathematical form of Eq. (53a, b) for the CTC model is: where the generic sets of shape functions u mn and the following sets of components for the reduced configuration variables appear: reduced displacement: U mn V mn W mn f g ð55aÞ reduced temperature: In the context of a Galerkin discretization, the minimal reduction for model CTC will lead to a system of three coupled nonlinear ODEs of variable order, depending on the multiphysics character of the problem.
To simplify the reduction procedure, a static condensation [25,26] of the in-plane u, v displacement is accomplished. This corresponds to what done also in the case of thermoelastic orthotropic plate [7] where, however, the treatment makes use of the Airy function. We refer to a rectangular symmetric cross-ply laminate with simply supported, immovable and isothermal edges, subjected to uniform stretching forces of magnitudes p x = p y = p in x and y direction. It should satisfy the following boundary conditions: As in [7], for the configuration components w, T 0 and T 1 in Eq. (54), the following single mode approximations are assumed: T 0 ¼ T R0 ðtÞ sin px a sin py b ð58bÞ where sin(px/a)sin(py/b) is a shape function satisfying both mechanical and thermal boundary conditions. In turn, for the configuration components u and v (inplane displacements) in Eq. (54), the following basis, ensuing from static condensation and satisfying identically the kinematic boundary conditions (56) and (57), is used: Reducing Eq. (82a, b) in Appendix for the case of symmetric cross-ply laminate, neglecting inertia terms and eliminating body force, and substituting Eqs.
Reducing Eq. (82c-e) for the case of symmetric cross-ply laminate, neglecting rotational inertia and adding a linear viscous term dw v (where d is the damping coefficient) in Eq. (82c), and substituting Eqs. (58) and (59), Galerkin integration provides the following equations: Again, for reasons of brevity, the expressions of the coefficients a ij are not reported.
In order to non-dimensionalize the system equations, the following transformations are introduced [7]: where a medium value for the thermal expansion is considered: Substituting Eq. (62) into Eq. (61), and considering an harmonic transverse excitation F ð0Þ 3 ¼ f cos xt and two harmonic thermal loads E ð0Þ ¼ e ð0Þ cos xt and E ð1Þ ¼ e ð1Þ cos xt, we have (the overbar upon the variables is dropped for convenience): with the adimensional parameters in Eq. (65) being defined in Eq. (83) in Appendix. In the particular case of orthotropic plate, Eqs. (64) and (65) coincide with the analogous equations in [7] (the correspondence between the dimensionless parameters in Eq. (65) and the dimensionless parameters in [7] is shown in Eq. (84) in Appendix). In [7] it has been shown that the one-mode Galerkin expansion may represent a sufficiently accurate approximation for the governing equation of the simply supported thermomechanically coupled rectangular orthotropic plate undergoing large deflections; comparisons are made with [22,[27][28][29].

Thermal contributions to the mechanical equation
The following numerical investigations will be aimed at preliminarily highlighting the influence of thermal dynamics on plate nonlinear vibrations. To this purpose, an analytical framework useful to drive the investigation and analyze the results can be obtained by algebraically solving (64b) and (64c) for T R0 and T R1 , respectively a 32 cosðsÞ ð 66bÞ and by inserting them into (64a), thus ending up with an ''expanded'' mechanical oscillator : where From Eqs. (67) and (68), it is observed that the thermal equations in (64) are reflected in the mechanical equation through the following contributions: In general, the time-varying unknown coefficients in Eq. (67) are obtained by coupling Eq. (67) with Eq. (66a, b) (i.e., by considering the whole system (64)). Moreover, the following particular cases can be considered.
3) Removing the whole thermal dynamics, while still keeping a difference between environmental and reference temperatures (T 1 6 ¼ 0), (66a) and (66b) reduce to: and Eq. (67) reads: governing by itself the nonlinear plate vibrations. 4) Removing all thermal aspects, Eq. (67) reduces to: that is the classical Duffing equation for nonlinear plate vibrations, directly achievable through a single-mode discretization of the sole mechanical variable W under the assumption of heat-insulated system. 5) Several intermediate cases could be considered, such as the one in which thermal dynamics is activated and sustained by thermal loads, while still removing mechanical coupling terms in the thermal equations (64b, c); in this case Eq. (67) loses the coefficients " b 31 and " b 23 ðsÞ and governs itself the mechanical problem, once obtained the residual time-varying unknown coefficients by separately solving the thermal problem.
It is observed that the thermo-mechanical coupling is induced by the coefficients " b 23 ðsÞ and " b 31 . Of course, if the adimensional parameters " b A 2 and " b D 2 (Eq. 83 in Appendix) are zero, this implies " a 25 ¼ " a 33 ¼ 0 in Eq. (66) and then " b 23 ðsÞ ¼ " b 31 ¼ 0 in Eq. (67) (cf. Eq. 68), and the coupling is physically canceled.

The effect of thermal aspects on nonlinear vibrations
In the following, to get some hints on the possible importance of the variable contribution of thermal effects, the responses obtained by solving the coupled Eqs. Runge-Kutta method of fourth order is employed to perform numerical investigation in free and forced vibrations. The time step of the numerical integration is chosen to be p/300, the error tolerance is set to be less than 0.0001.

Free vibrations
The main effect of thermomechanical coupling is the thermoelastic damping phenomenon which causes the amplitude of vibration decay [22]. To assess the sole effect of thermomechanical coupling in the free oscillations of the plate, we eliminate from system (64) the mechanical damping (" a 12 ¼ 0), the in-plane mechanical pre-load (" a 14 ¼ 0), the mechanical and thermal forcings (" a 18 ¼ " a 27 ¼ " a 34 ¼ 0) and the difference between environmental and reference temperatures (T 1 ¼ 0). So Eq. (67) becomes: and works with Eq. (66a, b) (also suitably shortened) to evaluate the unknown time-varying coefficients " b 22 ðsÞ, " b 23 ðsÞ and " b 32 ðsÞ. The corresponding model without thermomechanical coupling is obtained from Eq. (75) suitably shortened: Different lamination schemes are considered (Fig. 5): the plate R is orthotropic and characterized by the values of the dimensionless parameters in [8] 1 ; the plate S, also orthotropic, differs from R for the values of some parameters 2 ; the other plates present multilayer lamination schemes consisting of laminae type R and S. Table 2 in Appendix shows the values of the dimensionless parameters in Eq. (83) for each of these plates.
By way of example, for the S-R-S laminate the time responses obtained by the complete coupled model and by the uncoupled one (77) are shown in Fig. 6, starting from the initial deformed configuration Wð0Þ ¼ 1. Figure 6a and b show the plate deflection in the case of coupled and uncoupled model, respectively; Fig. 6c and d show the thermal dynamics of the coupled model. In Fig. 6a, the typical thermoelastic damping causing vibration amplitude decay [22] is observed: it is due to the non-trivial value of T R1 (Fig. 6d), which affects the terms " b 31 _ W and " b 32 ðsÞ in Eq. (76). In contrast, the very small value of T R0 (Fig. 6c) entails a negligible contribution of thermal membrane dynamics, as already noticed in [8] for the orthotropic plate. Indeed, the time histories of the associated terms " b 22 ðsÞ and " b 23 ðsÞ in Eq. (76) are nearly identical but 180°out of phase (Fig. 7, left side), thus producing a very small global contribution (Fig. 7, right side) to the frequency of the oscillator.
Such a situation occurs for all plates of Fig. 5, for which Eq. (76) can be written that works with Eq. (66b) (also suitably shortened) to evaluate the unknown time-varying coefficient " b 32 ðsÞ. The system is then referable to the simplified case 2) of Sect. 4.2 (in the free vibration framework), and the thermomechanical coupling is now induced only by the thermoelastic damping coefficient " b 31 . Table 1 shows the values of the thermoelastic damping coefficient " b 31 for the plates in Fig. 5, while Fig. 8 shows the time evolution of the displacement amplitude peak at the center of each laminated plate.  Fig. 5 Different multilayer lamination schemes consisting of orthotropic laminae type R and S 1 Note that the correspondence between the dimensionless parameters in Eq. (65) and the dimensionless parameters in [7,8] is shown in Eq. (84) in Appendix. 2 The values of " a 1 ¼ " d 1 e " a 3 ¼ " d 3 are three times greater than the corresponding values of plate R, while the values of " b A 2 ¼ " b D 2 are ten times greater.
It can be seen that the highest damping occurs for plate S-S-S-S-S (Fig. 8), whose dimensionless parameters values maximize the thermoelastic damping coefficient " b 31 in Eq. (78) ( Table 1), while the lowest damping is exhibited by the plate R (which has the lowest value of " b 31 ). In this respect, reminding that " b 31 ¼ À " a 33 " a 32 " a 16 (Eq. (68)), it is worth noting that the dimensionless parameters of the system may have a considerably variable importance. For example, it can be numerically verified that " k A 1 has a negligible effect  h=2k A 11 (Appendix, Eq. 83), and also the global conductivity in the x direction does not depend on the lamination scheme, because of its symmetry. Overall, the amount of thermomechanical damping " b 31 may meaningfully differ based on both material properties and lamination scheme.
After the transient, very minor differences are seen to occur among the two vibrational responses, highlighting how the effect of coupled thermal dynamics is definitely poor, as somehow expected. The low membrane thermal contribution results from the low values of T R0 (Fig. 9c); as a consequence, the terms " b 22 ðsÞ and " b 23 ðsÞ barely affect the frequency of the ''expanded'' mechanical oscillator (67) through   Table 1 Fig. 7-like diagrams. As to the bending thermal dynamics, its effect could appear not negligible if looking at the significant amplitude of T R1 , compared to T R0 (Fig. 9d); however, it is actually low as revealed by Eq. (64a), where T R1 is multiplied by the very small coefficient " a 16 ¼ 0:17. Such a low effect can also be realized thinking of the ''expanded'' mechanical oscillator (67). In fact, therein, the thermoelastic damping coefficient ( " b 31 ¼ À2:9) significantly lowers the total damping value (being " a 12 ¼ À0:05), but simultaneously the time-varying coefficient " b 32 ðsÞ adds to the mechanical forcing (Fig. 10, l.s.), and results in a sort of augmented overall excitation (Fig. 10, r.s.); the response outcome of the coupled oscillator is very similar to that of the uncoupled Eq. (75) The case of also a non-vanishing difference between environmental and reference temperatures (T 1 [ 0, i.e. plate temperature higher than environment one), which implies b 21 6 ¼ 0 in Eq. (67) and induces thermal conduction, is considered in Fig. 11.
The influence of the whole thermal dynamics on the mechanical one results in a decrease of the deflection, until the plate reaches the environment temperature, after which a forcing-sustained periodic oscillation is observed (Fig. 12).
The diagrams in Fig. 11a-c are well approximated by the simplified model accounting for the sole membrane thermal dynamics, case 1) of Sect. 4.2 (Eqs. 70-69a) (not plotted because indistinguishable from those of the complete coupled model); thus, the bending thermal dynamics plays a negligible role in the mechanical system response. Moreover, in the membrane thermal dynamics the mechanical coupling term is negligible, as shown below.
The global effect of all thermal membrane coefficients on the frequency in Eq. (67) can be evaluated looking at the time history (Fig. 13) of the following quantities: Fig. 9Þ ð79aÞ Fig: 11Þ ð79bÞ The equal sign in Eq. (79) ensues from comparison of equations (66a)-(68), and highlights that the global   Figure 13 shows these quantities in both short, s 2 0; 500 ½ , and long, s 2 0; 20000 ½ , time ranges. The time histories exhibit small oscillations (Fig. 13, l.s.) near the axis of thermal equilibrium T R0 ¼ 0 (case T 1 ¼ 0), or near a backbone curve (case T 1 [ 0) that lowers with time until the plate reaches the environment temperature (Fig. 13, r.s.). These oscillations are due to the mechanical coupling term " a 25 _ W Á W in thermal equation (64b), while the backbone curve represents the pure thermal conduction due to the term " a 26 T 1 (in the same equation).
Because of the small coupling oscillations observed in Fig. 13, the influence of the cyclic modulation of the frequency (b 23 term in Eq. (67)) appears to be almost negligible compared to the monotonic displacement of the thermal equilibrium point due to pure conduction. This can be seen by comparing the response in Fig. 13 with that of a mechanical single-dof oscillator: where T R0 (s) is obtained (in the case T ? [ 0) by separately integrating the thermal equation (64b) after having neglected the mechanical coupling term " a 25 _ W Á W (this system is referable to the simplified case 5) of Sect. 4.2). The responses of the two models are compared in Fig. 14 into ranges s 2 500; 520 ½ and s 2 10000; 10020 ½ , and show a good agreement. Thus, the main thermal effect on the nonlinear vibrations turns out to be the gradual membrane cooling of the plate produced by pure conduction. This causes a progressive increase of tensile stresses due to thermal contraction, which changes the frequency. The overall outcome is a reduction of the amplitude of mechanical oscillations (see Fig. 11), while the effect of coupled thermal dynamics is still poor.
It is worth noting that the periodic response of the coupled model in Fig. 12 can be well approximated by  the mechanical single-dof oscillator (74) (simplified case 3), with " a 14 ¼ 0, of Sect. 4.2), whose term " b 21 represents an assigned thermal distortion which is the final effect of the conduction phenomenon observed in Fig. 11. Figure 15 shows the effect of an in-plane mechanical pre-load " p ¼ 8 (added to the other excitation sources in Fig. 11), which provides the contribution " a 14 to the frequency in Eq. (67) (added to the thermal coefficients of Eq. (79b)). We see the steady-state response, when the plate has reached the environment temperature. The tensile stresses induced by T 1 [ 0 are not large enough to effectively balance the mechanical pre-load, so the plate vibrates only around the buckled equilibria. The figure also shows the response given by the uncoupled model (74), which shows slight differences with respect to the coupled one.
Finally, Fig. 16 shows the effect of adding the sole thermal forcing to the mechanical transverse forcing. The bending thermal forcing " e ð1Þ ¼ 5 (" e ð0Þ ¼ 0) has practically no effects, the response of the coupled model coinciding with the case " e ð1Þ ¼ 0 (Fig. 9a, b). Instead, a significant difference occurs in the presence of membrane thermal forcing " e ð0Þ ¼ 5 (" e ð1Þ ¼ 0). In this case, the global effect of all membrane coefficients on the frequency in Eq. (67) can be evaluated looking at the time history of the Eq. (79)-like quantity:  which entails a significant cyclic modulation of the frequency (Fig. 16, r.s.). The main contribution to this modulation is provided by the forcing term in the thermal equation (64b), while the mechanical coupling term (in the same equation) is negligible. Again, this may be shown using the simplified case 5) of Sect. 4.2, i.e., using Eq. (80) where T R0 (s) is obtained by separately integrating the thermal equation (64b), after having neglected the mechanical coupling term " a 25 _ W Á W; the response obtained in this way practically coincides with the blue graphs in Fig. 16.
Based on the forced periodic vibrations observed in the variety of the (though limited) number of cases examined, a few summary points can be made as regards the response phenomenology and the modeling in the background.
(i) Also in the presence of a non-vanishing heat flow (T 1 6 ¼ 0) and/or of explicit thermal sources (e ðiÞ 6 ¼ 0), there is no really significant effect of thermoelastic coupling on mechanical vibrations. (ii) The membrane thermal dynamics plays the dominant role (with respect to the bending thermal dynamics) in the mechanical system response, both when associated with a thermal conduction phenomenon and when directly activated by an external thermal source. (iii) In a number of cases, some main aspects of the thermally induced mechanical response can be suitably caught by referring to variably simplified models.
More systematic analyses are of course needed to further substantiate the previous points.

Conclusions
In the context of 2D engineering multiphysics modeling, a unified formulation of the thermomechanical problem of laminated plates undergoing finite amplitude vibrations has been presented. It integrates mechanical and thermal aspects, which are addressed in parallel by defining generalized 2D variables and equations also for the latter. The formulation virtually embeds a multitude of possible models, resulting from different assumptions about the plate mechanical and thermal configurations. It is a suitable framework to be referred to for systematic multiphysics analyses of plate nonlinear dynamics, owing to its favorable features of systematicity, awareness of theoretical assumptions, easiness of comparison of different models, and effective pathway for computer implementation. The scheme has been specialized to the laminated thermoelastic plate with von Karman nonlinearities.
The obtained continuous model has been then subjected to a minimum reduction via Galerkin procedure, ending up with a reduced model of three coupled ODEs. To get first hints on the possible importance of different thermal aspects in free and forced nonlinear vibrations, the numerical responses obtained by the complete reduced order model have been compared with those provided by simplified models variably accounting for the presence of mechanical and/or thermal excitations. The main effect of thermomechanical coupling is the thermoelastic damping causing some amplitude of vibration decay. In this respect, the considered numerical cases highlight that, in free nonlinear vibrations, the thermal membrane dynamics of the coupled model has a negligible influence with respect to thermal bending dynamics. In the case of periodic forced vibrations, for the considered numerical cases the effect of thermomechanical coupling is poor, while typical effects of the decoupled membrane thermal dynamics are recognized. Systematic forced nonlinear dynamics analyses are presently going on. The generality of the presented unified framework can be fruitfully exploited to conduct further studies concerned with many different aspects of system behavior, which range from the comparison of a variety of higher order, geometrically nonlinear, thermoelastic, plate models, up to the inclusion of other multiphysics aspects (e.g. piezoelectric effects, magnetoelastic coupling, and so on) in the overall theoretical treatment.