Thermodynamic Consistency of Plate and Shell Mathematical Models in the Context of Classical and Non-Classical Continuum Mechanics and a Thermodynamically Consistent New Thermoelastic Formulation

Inclusion of dissipation and memory mechanisms, non-classical elasticity and thermal effects in the currently used plate/shell mathematical models require that we establish if these mathematical models can be derived using the conservation and balance laws of continuum mechanics in conjunction with the corresponding kinematic assumptions. This is referred to as thermodynamic consistency of the mathematical models. Thermodynamic consistency ensures thermodynamic equilibrium during the evolution of the deformation. When the mathematical models are thermodynamically consistent, the second law of thermodynamics facilitates consistent derivations of constitutive theories in the presence of dissipation and memory mechanisms. This is the main motivation for the work presented in this paper. In the currently used mathematical models for plates/shells based on the assumed kinematic relations, energy functional is constructed over the volume consisting of kinetic energy, strain energy and the potential energy of the loads. The Euler’s equations derived from the first variation of the energy functional for arbitrary length when set to zero yield the mathematical model(s) for the deforming plates/shells. Alternatively, principle of virtual work can also be used to derive the same mathematical model(s). For linear elastic reversible deformation physics with small deformation and small strain, these two approaches, based on energy functional and the principle of virtual work, yield the same mathematical models. These mathematical models hold for reversible mechanical deformation. In this paper, we examine whether the currently used plate/shell mathematical models with the corresponding kinematic assumptions can be derived using the conservation and balance laws of classical or non-classical continuum mechanics. The mathematical models based on Kirchhoff hypothesis (classical plate theory, CPT) and first order shear deformation theory (FSDT) that are representative of most mathematical models for plates/shells are investigated in this paper for their thermodynamic consistency. This is followed by the details of a general and higher order thermodynamically consistent plate/shell thermoelastic mathematical model that is free of a priori consideration of kinematic assumptions and remains valid for very thin as well as thick plates/shells with comprehensive nonlinear constitutive theories based on integrity. Model problem studies are presented for small deformation behavior of linear elastic plates in the absence of thermal effects and the results are compared with CPT and FSDT mathematical models.


Introduction
The plate/shell theories and their mathematical models used currently are natural extensions of beam bending theories. The plate bending theories based on Kirchhoff hypothesis (Classical Plate Theory, CPT), first order shear deformation assumption (First order Shear Deformation Theory, FSDT) and the Higher order Shear Deformation Theory (HSDT) are in fact derived using the concepts used in Euler-Bernoulli beam theory, Timoshenko beam theory and higher order beam theories [1] Poisson (1829) and  in references [6]- [12]. A formal presentation of the plate bending theory now referred to as classical plate theory (CPT) is due to  in references [11] [12]. This plate bending theory is also referred to as plate theory based on Kirchhoff hypothesis (kinematic assumptions) that consists of: 1) Normals perpendicular to the middle surface of the plate before deformation remain normal after deformation to the deformed middle surface.
2) The normals perpendicular to the middle surface are assumed inextensible.
3) Based on (1), the normals to the middle surface before deformation rotate during deformation such that they remain normal to the deformed middle surface.
These assumptions are incorporated in the description of displacements ( ) 2) its failure to accommodate deformation physics of thick plates etc. is well known.
In an attempt to correct the shortcomings of CPT many papers, books and monographs have appeared [13]- [21]. The first significant improvement was due to what is now referred as FSDT in which 3 1 u x ∂ ∂ and 3 2 u x ∂ ∂ in the kinematic assumptions for 1 u and 2 u are replaced by unknown rotations 1 φ and 2 φ (similar to Timoshenko beam theory). This resulted in nonzero transverse shear stresses that are constant across the plate cross-section (this is well known to be non physical). The HSDT and their many variations are now available in the published works that alleviate the problems associated with CPT and FSDT.
Other published works related to the plate bending theories can be found in ref- erences [22] [23] [24] [25]. In a recent paper, Surana et al. [26] investigated thermodynamic consistency of the currently used beam mathematical models that are based on kinematic assumptions and are derived either using energy methods or using the principle of virtual work. The authors concluded that the currently used beam models cannot be derived using the conservation and balance laws of CCM or NCCM. Hence, these mathematical models are thermodynamically inconsistent. A serious consequence of thermodynamic inconsistency is that any further enhancement of these models cannot be done using the conservation and balance laws of CCM or NCCM. Almost all currently used mathematical models for plates and shells are derived based on specific kinematic assumptions defining displacements u 1 , u 2 , u 3 as functions of x 1 , x 2 , x 3 and t in conjunction with the following: 1) By constructing a functional (I) consisting of kinetic energy, strain energy due to bending and potential energy of the loads.
2) The first variation of this functional is set to zero ( 0 I δ = ). This is a necessary condition for an extremum of the functional I in (1).
3) The Euler's equations derived from 0 I δ = (using either fundamental lemma of the calculus of variations or the fourth basic lemma ( [27] [28] [29]) constitute the mathematical model (generally referred to as equations of equilibrium) of the deforming plate or the shell. The mathematical model for most plates and shell theories is derived using this approach. The CPT, FSDT, etc. mathematical models are derived using this approach.
Majority of the published works on the mathematical models for plates and shells [13]- [21] either follow this approach described in (1) - (3) or principle of virtual work. One could show that when the deformation due to mechanical work is reversible, both approaches yield identically the same mathematical models. The fundamental differences in various plate and shell theories and the corresponding mathematical models arise due to: sentative mathematical models for establishing thermodynamic consistency of all plate and shell mathematical models.
4) The conservation and balance laws and the associated constitutive theories of continuum mechanics for isotropic and homogeneous solid continua in 3  (such as plate/shell) must be used to describe physics of deformation to ensure that the deforming continua is always in thermodynamic equilibrium. Thus, to establish thermodynamic consistency (or lack of it) of the mathematical models of CPT and FSDT, we must begin with the conservation and balance laws of continuum mechanics and then incorporate the assumed kinematic relations of CPT and FSDT in these to arrive at the final mathematical models. The mathematical models so derived using laws of thermodynamics are obviously thermodynamically consistent and will honor the kinematic assumptions of CPT and FSDT. These mathematical models are compared with currently used mathematical models for CPT and FSDT. a) If the mathematical models of CPT and FSDT derived in this paper using conservation and balance laws along with their corresponding kinematic assumptions are exactly the same as those used currently, then the currently used mathematical models for CPT and FSDT are thermodynamically consistent.
b) If we find that the currently used mathematical models for CPT and FSDT are not the same as the mathematical models derived here in 4(a), then the currently used mathematical models for CPT and FSDT are thermodynamically inconsistent. In this case any further enhancement (such as incorporating dissipation and memory mechanisms) in the currently used mathematical models for CPT and FSDT are not possible using the conservation and balance laws of continuum mechanics. We keep in mind that the mathematical model derived using the conservation and balance laws of continuum mechanics may be invalid too due to enforcing the corresponding kinematic assumptions in their derivations, but this does not effect the issue of determining thermodynamic consistency of the current CPT and FSDT mathematical models. c) If we find that 4(b) holds, then obviously there is a need for thermodynamically consistent mathematical model describing physics of bending of plates/shells based on the conservation and balance laws of continuum mechanics that can describe thin as well as thick plates/shells without violating thermodynamic consistency. This is accomplished in the new formulation presented in the paper.
5) The outline of the work presented in this paper is given in the following a) The derivations of currently used CPT and FSDT mathematical models using energy functional are described and the complete mathematical models are presented and discussed. b) Equations resulting form the conservation and balance laws and constitutive theories are presented for i) Classical continuum mechanics. ii) Non-classical continuum mechanics using internal rotation [31] [32]. iii) Non-classical continuum mechanics using internal and Cosserat rotations [33] [34]. American Journal of Computational Mathematics c) The conservation and balance laws of classical and non-classical continuum mechanics (CCM, NCCM) are used in conjuction with kinematic assumptions of the CPT and FSDT to derive the resulting mathematical models for CPT and FSDT based on CCM as well as NCCM conservation and balance laws. These mathematical models honor the conservation and balance laws and kinematic assumptions.
d) The mathematical models derived in 5(c) are compared with the currently used mathematical models of CPT and FSDT. e) Derivation of the new kinematic assumption free thermodynamically consistent plate/shell mathematical model is presented that is free of kinematic assumptions and is capable of describing the thin as well as thick plate/shell deformation physics.
f) Model problem studies comparing results from currently used CPT, FDST plate models and from the new formulation presented in this paper are given. g) Summary and conclusions are given in the last section of the paper. h) "Appendix A" contains complete mathematical models based on CCM and NCCM i.e., the equations resulting from the CCM as well as NCCM conservation and balance laws, the basic definitions of rotations (both internal and Cosserat) and the constitutive theory considerations for CCM as well as NCCM based on internal rotations and Cosserat rotations.  Figure 1). In this paper we use the same notations as used in the currently published works on plates and shells. Some of these notations are not consistent with the notations used in continuum mechanics. Based on continuum mechanics notations, ij σ refers to Cauchy stress on the i plane in the j direction. Fortunately, this notation is same in plate and shell mathematical models. Likewise ij m is the component of the Cauchy moment tensor acting on the i plane about the j axis. Unfortunately, this notation is not followed in the derivations of plate and shell models. We will illustrate the differences when discussing the various mathematical models for plates and shells. To avoid confusion, we present currently used plate and shell mathematical models using the same notations as used currently in the literature.

Classical Plate Theory (CPT) Based on Energy Functional (or the Principle of Virtual Work)
The classical plate theory is due to Kirchhoff hypothesis [11] [12] and is based on the following kinematic assumptions for the displacements 1 2 3 , , u u u at an arbitrary point Based on (2), the strain energy is only due to 11 is derived by considering energy functional (I) consisting of kinetic energy, strain energy and the potential energy of loads over length L and the area of cross-section A i.e. integral over the volume (V) of the plate and time [0, t]. Using Lagrangian description, we can write the following for I: We use (1) -(3) in (4) and expand each term in (4). For an extremum of I, we set first variation of I to zero ( 0 I δ = ). Using integration by parts all differentiation is transferred from 1 The shear forces Q 13 and Q 23 are obviously zero. The complete mathematical model for CPT consists of (5), (6) and (7) in 1 2 , u u   and u 3 in which M 11 , M 22 and M 12 are functions of u 3 and are defined by (9). N 11 , N 22 and Q 12 can be easily obtained using (8) and (3). For pure bending without inplane load and deformation in x 1 x 2 plane, the mathematical model consists of (7) and (9). If we only consider stationary process (boundary value problem) and if we substitute (9) in (7), then we can obtain a single fourth order partial differential equation in u 3 describing pure bending of plates.

Remarks.
1) As is well known, in this mathematical model the transverse shear stresses 23 σ and 31 σ (hence the corresponding shear forces) are zero. This is obviously non physical.
2) Because of (1), this mathematical model does not account for strain energy due to these shear stresses and the corresponding shear strains, hence we expect this model to underestimate the actual transverse displacement u 3 .
3) We note that moment M 11 due to 11 σ is about x 2 axis, hence based on standard continuum mechanics notation it should have been M 12 (on face x 1 about x 2 axis). The same holds for other moments as well. We remark that these moments are not components of a moment tensor. We continue with this nota-tion for the currently used mathematical models presented in this paper for the sake of transparency between the material in this paper and the published works.

First Order Shear Deformation Theory (FSDT) Based on Energy Functional (or the Principle of Virtual Work)
The mathematical model for the first order shear deformation theory employs the following kinematic assumptions for the displacements u 1 , u 2 and u 3 at an arbitrary location 1 2 3 , , x x x in the plate for an arbitrary time t [14] [18].
1) We note that the transverse shear forces Q 13 and Q 23 are not zero in this mathematical model.
2) This mathematical model when compared with CPT is expected to yield higher transverse deflection u 3 due to shear deformation (additional strain energy due to transverse shear strains 13 ε , 23 ε and shear stresses 13 σ and 23 σ ).
3) We note from (17) that transverse shear forces Q 13 and Q 23 are not functions of x 3 i.e., these are constant along the thickness of the plate where as their actual distributions are parabolic functions of x 3 . Thus, the strain energy due gives rise to three partial differential equations in u 3 , 1 φ and 2 φ , containing up to second order derivatives of u 3 , 1 φ and 2 φ .

Kinematic Assumption Free Methodology for Bending of Plates and Shells
The derivations of the mathematical models for bending of plates and shells that are free of kinematic assumptions should always begin with the use of conservation and balance laws as this is necessary for thermodynamic consistency of the resulting mathematical models. We present a brief summary of the conservation and balance laws of classical continuum mechanics (CCM) as well as non-classical continuum mechanics (NCCM). In case of CCM, for small strain, small deformation and isothermal reversible mechanical deformation physics, we only need to consider balance of linear momenta, balance of angular momenta and the linear constitutive theory for the Cauchy stress tensor ("Appendix A", Section A.1). However, in case of NCCM, additional considerations are necessary as shown in "Appendix A", Section A.2 and Section A.3. In the following sections we present explicit forms of the conservation and the balance laws in 1 2 3 , , x x x space and time t as needed in the derivations of the CPT and FSDT plate bending mathematical models.

Classical Continuum Mechanics
Complete details of the mathematical model consisting of balance of linear mo-American Journal of Computational Mathematics menta, balance of angular momenta, linear strain measures and the linear constitutive theory for the Cauchy stress tensor are given in "Appendix A", Section A.1.

Non-Classical Continuum Mechanics (NCCM)
In non-classical continuum mechanics the description of the physics of deformation of solid continua considers additional aspects of deformation physics over and beyond classical continuum mechanics. In the present work, we consider two non-classical continuum theories that are directly applicable in the derivation of the mathematical models of the CPT and FSDT due to the kinematic assumptions used.

Non-Classical Continuum Mechanics with Internal Rotations
In this continuum theory the entire displacement gradient tensor, i.

Non-Classical Continuum Mechanics with Internal and Cosserat Rotations
In this continuum theory, we also consider internal rotations ( ) Linear constitutive theories for s σ , m , a σ and q are given by [ Equations (A2), (26), (28) -(30) are a system of twenty one partial differential equations in twenty one dependent variables:

Thermodynamic Consistency of the Currently Used Plate and Shell Mathematical Models
In this section we attempt derivations of CPT and FSDT mathematical models using the conservation and balance laws and constitutive theory of classical continuum mechanics as well as non-classical continuum mechanics based on internal rotations due to antisymmetric part of the displacement gradient tensor as well as internal rotations and Cosserat rotations in conjunction with the kinematic assumptions used in CPT and FSDT mathematical models.

CPT Mathematical Model Derivation Using CCM
We initiate the derivation by considering the conservation and balance laws and the constitutive theory in Section (A2) of "Appendix A". These are valid for the deformation of continuous matter in 3  and ensure thermodynamic consistency. Thus, use of these is justified for bending of plates and shells in 3  . The currently used mathematical model for CPT requires that we incorporate kinematic assumptions (1) in the conservation and balance laws of CCM. From (1) the expressions for the strain measures 11 ε , 22 ε and 12 ε are obtained as shown in (2) and we also find that the 33 0 ε = , 23 We substitute from (2) and (3) in the balance of linear momenta (A2), and consider integral over the volume The resulting integral is valid over arbitrary area A, hence the integrand must be zero, which gives us the following equations: where, 2  2  2   11  11  3  22  22  3  12  12  3  2  2  2   /2  2  2  33  33  2  3  3  31  32  2  2  3 3 Equations (32) are three partial differential equations in 1 2 , u u   and u 3 . (32) is obviously not the same as the currently used mathematical model (5) - (7). Thus, we can conclude that the mathematical model of Section 2.1 (currently used) is thermodynamically inconsistent based on CCM as it cannot be derived using the conservation and balance laws of CCM.

Remarks 1) This mathematical model
2) If we only consider plane stress deformation in the plane of the plate as considered in almost all currently used theories, then 33 0 N = . For this case, the balance of linear momenta (third equation in (32)) contains no deformation physics related to the internal stress field.
3) Transverse shear stresses 23 σ and 31 σ and hence the corresponding transverse shear forces Q 23 and Q 31 are zero (as also in the case of currently used models) for CPT. This is obviously non physical. (7) except the first and the last terms, the remaining terms cannot be justified based on the derivation using CCM (third equation in (32)).

5)
We observe that in this derivation we began with perfectly valid balance laws of CCM, a valid mathematical model for continuous deformation in 3  , but these are altered and damaged due to kinematic assumptions, finally resulting in a mathematical model that is not valid for the physics under consider-American Journal of Computational Mathematics ation.

Derivation of Mathematical Model for CPT Using NCCM Based on Internal Rotations
The motivation for the consideration of NCCM in deriving the CPT mathematical model becomes clear if we decompose displacement gradient tensor d J into symmetric and skew symmetric components. (34) in which components of d a J are in fact related to the internal rotations First, we define the internal rotations Using the kinematic assumptions (1), we obtain following for the internal rotations defined in (35) The components of d a J can now be defined using (36) The components s σ are same as those defined by (2). Substituting the total stress tensor σ from (38) in the balance of linear momenta and integrating with respect to 3 , we obtain the following: 2  1  11  21  31  0  0  1  2  1  2  3   2  2  12  22  32  0  0  2  2  1  2  3   2  3  13  23  33  0  0  3 12 21 31 13 23 , , , , Q Q Q Q Q and 32 Q are defined later.

The internal rotation gradient tensor
is decomposed into symmetric and skew symmetric tensors.   (46) and are same as those in (3). 33 33 0 s σ σ = ≠ if we consider the deformation of plate bending in 3  (discussed earlier). Using (47) and (42), explicit expressions for the constitutive theory for ij m can be obtained in terms of the symmetric part of the rotation gradient tensor.
Using these expressions for ij m in BAM (44) and dividing them by 2 and then integrating with respect to 3 We note that 4) The moment tensor m is symmetric due to balance of moment of moments balance law [37] [38] necessary in NCCM.
5) It is rather obvious that the physics considered in the present derivation is consistent with the kinematic assumptions. This is not the case in derivation based on the energy method or the principle of virtual work. 6) We comment that the moment tensor m in the present derivations is Cauchy moment tensor. This is not the case for

Derivation of Mathematical Model for FSDT Using CCM
In this derivation we need to incorporate kinematic assumptions (11), strains (12) and associated stresses (13) in the conservation and balance laws of CCM ("Appendix A").
Substituting (13) in the conservation and balance laws of CCM ("Appendix A") and integrating in 3 We note that , ,0, Q Q x x t = are functions of 1 2 , x x and t only.
Hence, this mathematical model is not valid.
2) Equations four and five in (15)  3) This mathematical model is obviously not the same as the currently used FSDT mathematical model presented in Section 2.2, and is also invalid due to lack of closure.

FSDT Model Derivation Using NCCM Based on Internal and Cosserat Rotations
The The total rotations Based on the kinematic assumptions (11), strains are defined by (12) and the components of the symmetric stress tensor are defined by (13) if we assume plane stress behavior in the plane of the plate. However, if we consider the plate deformation in 3  , then coefficients ij D are modified and are due to linear constitutive theory in 3  giving rise to 33 0 σ ≠ and American Journal of Computational Mathematics ( ) 1  2  33  31  32  33  1 2 3  1  2 , , , The stress tensor σ is decomposed into symmetric and skew symmetric tensors s σ and a σ : Definitions of the components of the s σ is same as in (13) 12 21 31 13 23 , , , , Q Q Q Q Q and 32 Q appearing in (39) are defined by (67) -(71). (53) and using (57):   (59) in which where µ , λ are Lame's constants and µ     3) This mathematical model is not only different from the currently used FSDT mathematical model, but is also invalid due to lack of closure. 3) From (1) and (2), all currently used plate/shell mathematical models that are based on kinematic assumptions and energy methods or principle of virtual work are thermodynamically inconsistent. That is, when these mathematical models are used to study deformation of plates/shells, thermodynamic equilibrium is not ensured during the evolution of the deformation.

4)
Since laws of thermodynamics cannot be used in deriving the currently used mathematical models, the dissipation, memory mechanism and any further enhancements in these mathematical models can only be done phenomenologically. As well known [30], the phenomenologically incorporated dissipation and memory mechanism are only possible in 1  and their extensions to 2  and 3  (required for plates/shells) is not possible. Due to the absence of energy equation and the entropy inequality thermodynamically consistent treatment of dissipation and memory mechanisms is not possible in currently used plate/shell mathematical models.

Kinematic Assumptions Free Methodology Plates and Shells
The derivations of the mathematical models for currently used plate and shell theories are based on kinematic assumptions that contain internal rotations that arise due to displacement gradient tensor and/or assumed Cosserat rotations American Journal of Computational Mathematics (that are additional unknown degrees of freedom) at a material point. We have shown that if the kinematic assumptions in plate/shell mathematical models are a requirement, then the currently used mathematical models for plates/shells cannot be derived using the conservation and balance laws of CCM or NCCM incorporating internal or the internal and the Cosserat rotations. Inadequacy of the energy methods in deriving the plate/shell mathematical models for enhanced deformation physics have already been discussed and is the strong motivator for the present work. Additional benefit of the new methodology presented in this paper is that a single formulation will be able to address the physics of all plate/shell deformations for linear as well as nonlinear elasticity regardless of the plate/shell thickness, loading and the boundary conditions (and the initial conditions). As is well known the conservation and balance laws of CCM and NCCM always yield thermodynamically consistent mathematical models, hence are always considered as a first step in all of the derivations presented in this paper for establishing thermodynamic consistency of the currently used plate/shell mathematical models. We have shown that it is only after incorporating the kinematic assumptions that the mathematical models lose thermodynamic consistency. We consider the following guidelines: 1) It is perhaps meritorious to consider a methodology in which the kinematic assumptions as they are used in the current theories are not considered (i.e., eliminated) as these are the main cause of thermodynamic inconsistency of the resulting mathematical model. Instead, we begin with mathematical model in 3  consisting of conservation and balance laws and incorporate the desired physics of the kinematics of deformation in a much more general manner so that: a) it maintains thermodynamic consistency, and b) the resulting formulation holds and remains accurate for slender as well as thick plate and shell deformation physics.
2) For continuous media CCM and NCCM are two possible approaches for deriving mathematical models of the deforming continuous matter. The CCM has been the foundation for describing the physics of deformation of continuous media. NCCM incorporates additional physics in the mathematical models over and beyond CCM. This may be beneficial in some instances. Nonetheless for isotropic, homogeneous linear elastic continuous matter, classical continuum mechanics must at least provide some reasonable description of the deformation physics regardless of 1  , 2  and 3  . Thus, we should be able to describe the plate/shell deformation physics using the conservation and balance laws of CCM as well as NCCM if we eliminate the kinematic assumptions at the onset of the derivations of the mathematical model(s). In the following, we consider the conservation and balance laws of CCM in describing the plate and shells deformation physics.
3) The finite element method is perhaps the most general and versatile technique of obtaining solution of the mathematical models consisting of partial differential equations. Such solutions are numerical but piecewise analytical and with higher order global differentiability local approximations, these solutions can truly approach theoretical solutions [28] [29]. Thus, if the finite element technique is eventually to be used to obtain solution of plate or shell mathematical models, then keeping this in mind we can begin with conservation and balance laws of CCM in 3  , followed by the design of a geometric description of the plate or the shell finite element that is identical to the manner in which plates and shells are geometrically described currently i.e., the geometric description of the midplane of the plate or the shell and the nodal vector at each point located at the middle surface given by the difference of the two ends which define the bottom and top surface of the plate or the shell.

4)
We map the geometry of the plate or the shell element in (3) in 3  into a natural coordinate space ξηζ in a two unit cube in which ξ and η are in the middle plane of the element and ζ is the transverse direction to the middle plane. The origin of the coordinate system ξηζ is at the center of the two unit cube ( Figure 2).

5)
We can design p-version hierarchical local approximation in ξ , η and ζ of p-levels p ξ , p η and p ζ with variable choice of the orders of the approximation space to ensure desired global differentiability of the approximations. In the local approximation, displacements 1 u , 2 u and 3 u can be approximated by polynomials of degrees p ξ , p η and p ζ as well as higher orders of differentiability. Choices of p-levels p ξ , p η and p ζ define approximations in the plane of the plate or the shell element as well as deformation of its cross section(s). By choosing p ξ , p η and p ζ and desired orders of global differentiability of the approximations i.e., the orders of the approximation space defining the global differentiability of the approximation for 1 u , 2 u and 3 u desired kinematic behavior of the mid plane of the plate or the shell as well as its cross-sections can be achieved. 6) Thus now we have: a) conservation and balance laws in 3  defining the mathematical model for the plate or the shell. b) Geometric description of the plate or the shell element and the local approximation for the displacements describing the deformation physics. The algebraic equation for the plate or the shell element are derived using the integral form based on Galerkin Method with Weak Form (GM/WF) or residual functional i.e., least squares method [28] [29]. For stationary processes i.e., boundary value problems and for linear elastic reversible deformation of solid continua, GM/WF yields unconditionally stable computational processes [28] [29]. In all other cases, least squares method based on residual functional can be shown to obtain unconditionally stable computations [28] [29]. 7) In the methodology presented here, desired kinematic requirements are achieved through local approximation by progressively increasing p-levels. In this methodology conservation and balance laws of CCM or those of NCCM can be considered. In the present work we consider CCM conservation and balance laws. temperature effects and linear constitutive theory) was originally derived by Surana et al. [39] [40] using principal of virtual work. Exactly similar derivation as in ref [39] [40] is also possible using energy method. This approach was also kinematic assumption free, but has the following major shortcomings and limitations. a) Principle of virtual work as used in reference [39] can only be used for reversible mechanical deformation, hence cannot be extended to include physics of dissipation and memory mechanisms.
b) The constitutive theory for Cauchy stress tensor is assumed to be linear generalized Hooke's law. c) Due to absence of energy equation, rate of work conjugate pairs needed for constitute theories cannot be established. The energy methods and principle of virtual work used in ref [39] [40] assumes a constitutive theory as there is no mechanism to derive it.
9) The methodology used here for thermoelastic plates and shells based on the conservation and balance laws of CCM is free of all these restrictions described in (8), yet results in the same formulation as due to principle of virtual work or energy method for reversible small deformation physics and linear constitutive theory in the absence of thermal effects.

Conservation and Balance Laws of CCM and Constitutive Theories
For thermoelastic deformation physics, the conservation and balance laws are given in the "Appendix A", Section A.1. We consider thermoelastic solid continua reversible mechanical deformation and general constitutive theories based on integrity [30] [41]- [57]. Based on the conservation and balance laws (more specifically second law of thermodynamics), the choice of Φ , η , σ and q as constitutive variables at the onset is rather straight forward. Using the conjugate pairs, ij ij σ ε and θ ⋅ q g and thermoelastic deformation, choice of ε , θ as ar-American Journal of Computational Mathematics gument tensors of σ and g and θ as argument tensors of q is appropriate.
Using the principle of equipresence, we can choose ε , g and θ as argument tensors for Φ , η . Thus, we have

Constitutive Theory of σ
Constitutive theory for Cauchy stress tensor σ can be derived either based on Helmholtz free energy density or using representation theorem. We consider both approaches here. The resulting constitutive theories from both approaches are same and are nonlinear in the components of ε but linear in θ when based on integrity [30].
1) Constitutive theory for σ using Helmholtz free energy density Φ using (75), we can write Substituting (77) in SLT (A5) and rearranging terms, we can write For (78) to hold for arbitrary but admissible choices of  ε , θ  and g  , the following must hold and SLT (78) reduces to Following the details of similar derivations in reference [30], we can obtain the following from (85) We substitute (88) in (86) and group the terms and coefficients and if we de- i.e., ; 1,2, , ; 2; the generators.
Then (86) can be written as [30]:  (86), the rest of the derivation remains same as in case of the derivation using Φ given above.

Constitutive Theory of q
We begin with ( ) ,θ = q q g based on entropy inequality and use representation theorem [30] [41]- [57] which gives us the following constitutive theory for q ( ) ( ) This constitutive theory is based on integrity, hence uses complete basis. From (94), we can also derive a linear constitutive theory where 1 2 , , κ κ κ are material coefficients defined in a known configuration Ω .
These can be functions of ( ) Ω ⋅ g g (invariant of g ) and θ Ω .

Final Mathematical Model (Thermoelastic)
The entropy inequality is identically satisfied with the choices we have made in deriving the constitutive theories.
Remarks 1) Use of this mathematical model (96) -(100) in 3  will result in a comprehensive nonlinear plate and shell formulation that describes evolution in a) The constitutive theory for σ is based on integrity.
b) The thermal effects are incorporated using the mechanism consistent with the conservation and balance laws.
2) In the formulation of reference [39] [40] based on principle of virtual work, the resulting Euler's equations representing the associated mathematical model In this approach, even though principle of virtual work permits non-linear reversible mechanical deformation, there is no mechanism of the constitutive theory like (99). Secondly, consideration of thermal effects requires consideration of energy equation (98) which is not possible in energy methods or principle of virtual work. Nonetheless, if we limit the physics to conform to what has been considered in reference [39] [40], then the mathematical model (96) -(100) obviously reduces to (101) and (102).
3) We clearly see that enhancement of mathematical model used in reference [39] [40] to (96) -(100) using energy methods or principle of virtual work is not always possible (as shown here) due to limitations of energy methods and principle of virtual work. 4) In this paper, CPT and FSDT mathematical models for BVPs are used as representative plate/shell mathematical models that are based on reversible mechanical small deformation without thermal effects. Thus, to compare the work presented here with CPT and FSDT, we need to consider the mathematical model consisting of (101) and (102) in 3  (simplified form of (96) -(100)).

Complete Mathematical Model
Expanded form of the mathematical model in 3  ((101), (102)) are given in the following:  ; ; The nonzero elements of (6 × 6) [D] material coefficient matrix are given by: The boundary conditions (BCs) are defined by:

Integral Form and Element Formulation
If we substitute stresses from (104) in balance of linear momenta (103), then we obtain a system of three partial differential equations in 1 u , 2 u and 3 u . One could easily show that the differential operator A and its adjoint * A in this system of differential equations are same, hence for this system of Equations x Ω is a typical element "e". Let 1 2 3 , , w w w be test functions such that  T  T   1  1  2  2  3  3 . , 0, . , 0, . , 0 We substitute from (103) 21  31  1  1  0 1  1  1  2  3   12  22  32  2  2  0 2  2  1  2  3   13  23  33  3  3  0 3  3  1  2  3 . , , Transfer one order of differentiation from the stresses to the test functions 1 w , 2 w and 3 w in (113) using integration by parts (IBP) to obtain  2  3   3  3  3  3  1  13  23  33  0 3  3  1  2  3   3  13  23 From (114), 1 2 3 , , u u u are primary variables (PVs) and the coefficients of  T  T  T  T   1  2  3 , , in which is the nodal load vector due to body forces. Assembly of the element equations, imposition of boundary condition and solution of resulting linear simultaneous algebraic equations follows American Journal of Computational Mathematics usual procedure [28] [29]. We note that T e e K K     =     , a consequence of the fact that adjoint of the differential operator in (103) is same as the operator itself. In order to be able to calculate

Plate and Shell Element Geometry and Local Approximation
We consider three dimensional hexahedral geometry shown in Figure 2(a). The plate/shell element geometry and the local approximation considered here are derived using element in Figure 2(a). Consider nine node configuration on the bottom and top surfaces (planes) of the hexahedron shown in Figure 2  . The shell element of Figure 2(b) is mapped into a natural coordinate space , , ξ η ζ in a two unit cube with the origin of the coordinate system at the center of the cube (Figure 2(d)).
Symbolically e Ω is mapped into m Ω , m Ω being the map of the element Thus, the coordinates at any arbitrary point ( ) ( ) The local approximation functions for the shell element are derived using tensor product of 1D functions in , ξ η and ζ directions keeping in mind that all dofs at the element nodes may not be the same, the tensor product process necessitates that we take tensor product of approximation functions and the nodal variable operators separately. After taking the tensor product, the resulting 3D nodal variable operators when act on the dependent variable(s) produce the corresponding nodal degrees of freedom that correspond to the approximation functions generated by taking the tensor product 1D approximation functions in , ξ η and ζ .
For the sake of simplicity we illustrate the process of deriving details of element local approximation using approximations of class The 0 C p-version local approximation for ψ in the ξ direction for the three node configuration (with nodes 1, 2 and 3 located at Similarly in the η direction for the three node configuration (with nodes 1, 2 and 3 located at

( )
In the ζ direction, the shell element nodes are located at We consider a map of 1 1 ζ − ≤ ≤ at each node in α coordinate system pointing in the same direction as ζ in which at each node (say i) (Figure 2(e)). This mapping is given by ψ ψ ψ ψ ψ ζ ζ Substituting from (130) in (131) ( ) Hence, and 2 1 for the node at 0 or 0 2 i ψ ψ ψ ζ α using (134) and (135) in (131) ( ) Likewise for 2 p ζ = , for the three node configuration in Figure 1, we can write (Lagrange interpolation) ( ) Substituting for ζ from (130) ( ) Following the same procedure, we can derive the following for a node i located at 0 ζ = for p-level of p ζ .
( ) we note that ( ) ψ ζ in (140) is hierarchical i.e. lower p-level approximation functions are a complete subset of the higher p-level approximation functions and the same holds for nodal degrees of freedom.
[ ] 1 2 3 , , , , The volumes in 1 2 3 , , 1, 2, , Using (145) [29]. Assembly of the element equations and their solutions follows standard procedure [28] [29]. Remarks 1) In this formulation p-levels in , ξ η and ζ define the plate/shell deformation behavior in the plane of the plate/shell as well as in the transverse direction to the shell middle surface.
2) Increasing p-levels in , ξ η and ζ result in higher degree approximations of 1 2 , u u and 3 u . While p-levels p ξ and p η primarily control the deformation in the plane, p ζ controls transverse deformation. C Ω in the weak sense. 4) In the derivation details presented here we have intentionally used curve geometry to emphasize that the formulation presented here is valid for flat plates/shells as well as shells of arbitrary geometry. 5) We remark that when the deformation is reversible and when the matter is isotropic and homogeneous, the mathematical model derived using principle of virtual work or energy methods in Lagrangian description (in the absence of kinematic assumptions) are same as those obtained using balance of linear momenta. Thus, for such deformation physics one could derive the finite element formulation directly using principle of virtual work or energy method. This in fact has been done in reference [39] [40] for three dimensional laminated composite curved shell elements. A serious drawback of this approach is that, in this approach only the equations corresponding to balance of linear momenta can be derived. Precise nature of the rate of work conjugate pairs is not possible either due to lack of second law of thermodynamics, hence derivations of a comprehensive constitutive theory is not possible. In the approach presented here using balance of linear momenta, derivation of constitutive theory is based on integrity and the finite element formulation based on VC integral form using GM/WF [28] [29] (as * A A = ) results in more comprehensive mathematical model as presented here.

Solutions of Model Problems
We consider a square plate ( a a × ) simply supported on all four boundaries. Midplane of the plate lies in 1 2 x x plane and the origin of fixed x-frame is at the bottom left corner. Plate thickness is "h". 3 x is normal to 1 2 x x -plane pointing upwards. Top surface of the plate ( 3 2 x h = ) is subjected to uniformly distributed load acting in the negative 3 x direction. A schematic of the plate is shown in Figure 3. We nondimensionalize density ρ , modulus Ê and distance (or displacement) by using reference density 0ρ ρ = , reference modulus     ). , the plate is thin, hence CPT and FSDT are expected to produce almost identical results (Figure 4(a)).

Discussion of Results
The new formulation is 3D elasticity formulation, hence contains more comprehensive and complete deformation physics compared to CPT and FSDT. Plot   Thus, at 1 2 5.0 x x = = graphs of ( ) σ .
Thus, in Figure 6(a) we only observe two sets of graphs. The discontinuity of ( ) 13 σ is due to discontinuity of actly similar behavior at lower p-levels. We observe discontinuity of 13 σ but with progressively increasing p-levels, the discontinuity of 13 σ diminishes and for p-levels in the range of 5 to 9, we obtain converged solutions for 13 σ . Higher thicknesses of plate naturally require higher p-levels. 23 σ versus 3 x at  33 σ , but no difficulties are encountered in doing so. In case of CPT, FSDT and HSDT plate mathematical models, 33 0 σ = . This is obviously non physical especially for thick plates. In case of the new formulation presented here 33 σ is calculated accurately regardless of the plate thickness.

Summary and Conclusions
Using currently used CPT and FSDT mathematical models as representative mathematical models for all plate/shell mathematical models, that are derived using energy methods or principle of virtual work and are based on kinematic assumptions, the thermodynamic consistency of all plate/shell mathematical models has been evaluated using the conservation and balance laws of CCM as well as NCCM based on internal, internal and Cosserat rotations. CPT utilizes internal rotations whereas FSDT uses Cosserat rotations in the kinematic assumptions but requires both internal and Cosserat rotations in the derivation of the mathematical model. Use of internal and/or Cosserat rotations in the kinematic assumptions is typical of all plate/shell mathematical models derived using the energy methods or the principle of virtual work. Thus, the findings in this work using CPT and FSDT mathematical models hold for all plate/shell mathematical models used currently.
It has been shown that CPT and FSDT mathematical models for plate bending with their respective kinematic assumptions cannot be derived using the conservation and balance laws of CCM or NCCM solely based on internal rotations and both internal and Cosserat rotations. Thus, we can conclude that all currently used plate/shell mathematical models are thermodynamically inconsistent i.e., the deformation resulting from the solutions of these mathematical models is in violation of the principle of thermodynamics of CCM as well as NCCM.
The currently used plate/shell mathematical models are derived using the energy methods and the principle of virtual work, so only reversible mechanical deformation can be included in these mathematical models without resorting to phenomenological approaches. Inclusion of dissipation and memory mechanisms in these mathematical models is only possible using phenomenological approach in 1  due to the absence of the energy equation and the entropy inequality. This approach cannot be extended to 2  and 3  . From the derivations of CPT and FSDT presented in the paper using CCM and NCCM, it is obvious that the main source of problem is a priori assumption of kinematic relations. These kinematic assumptions result in conflict of some type or other when used in conjunction with the conservation and balance laws of CCM or NCCM.
We have presented a kinematic assumption free thermoelastic plate/shell formulation in which the mathematical models consist of the conservation and balance laws of CCM in 3  : balance of linear momenta, balance of angular momenta (implies Cauchy stress tensor is symmetric) and the energy equation.
The constitutive theory for the Cauchy stress tensor is based on integrity (complete basis) and includes thermal effects. This constitutive theory is a non linear constitutive theory in terms of strain tensor components (contains terms up to fifth degree terms) but linear in temperature θ . The constitutive theory for the heat vector q is also based on integrity and is a non linear constitutive theory containing up to third degree terms of the temperature gradients. These balance laws are derived for evolutions (IVPs) but naturally reduce to mathematical models for BVPs including linear constitutive theories for σ and q if so desired. When the deformation process is isothermal, θ is no longer a dependent variable, energy equation is not needed and we only have BLM.
A finite element formulation is constructed using this mathematical model based on CCM for BVPs in which the constitutive theory is linear and the deformation process is isothermal so that the results from this formulation can be compared with CPT and FSDT. The finite element formulation is based on Strains, rotations and their gradients, conjugate pairs and constitutive theories In this appendix we present the conservation and balance laws of CCM as well as NCCM, some details of the strain tensor, rotation tensor, gradients of the rotation tensor and its decomposition, the conjugate pairs and linear constitutive theories for classical as well as non-classical continuum mechanics based on internal rotations and the non-classical continuum mechanics incorporating internal and the Cosserat rotations. The material presented here considers small deformation, small strain, linear as well as non linear reversible mechanical deformation.

A.1. Classical Continuum Mechanics (CCM)
A material point has only three translational degrees of freedom. The symmetric part of the displacement gradient tensor constitutes strain measures and the antisymmetric part of the displacement gradient tensor (a measure of internal rotations) is neglected in the derivation of the conservation and balance laws. We have conservation of mass, balance of linear momenta (BLM), first law of thermodynamics (FLT), the second law of thermodynamics (SLT) and considerations for constitutive theories [30]. x directions, ij σ is Cauchy stress tensor, ij ε is linear strain tensor, e is specific internal energy, q is heat flux vector, g is temperature gradient vector, φ is Helmholtz free energy density, η is entropy density and θ is temperature.
From the conjugate pairs , ij ij σ ε  and θ ⋅ q g , we can conclude that at the very least the following must hold (thermoelastic solid continua) Choices of the argument tensors for φ and η are discussed in section 4 in conjunction with the derivation of the constitutive theories for σ and q .