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     . The earlier presentations of the concepts and formulations related to the bending of plates can be traced back to Chladni et al. (1802), Germain (1826), Lagrange et al. (1828), Cauchy (1828), Poisson (1829) and Kirchhoff (1850) in references  - . A formal presentation of the plate bending theory now referred to as classical plate theory (CPT) is due to Kirchhoff (1850) in references  . 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 and at an arbitrary point using the rotations
and while displacement u3 remains a function of x1, x2 and t only. The problems associated with this plate theory are:
1) zero transverse shear stress (non physical)
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  - . The first significant improvement was due
to what is now referred as FSDT in which and in the kinematic
assumptions for and are replaced by unknown rotations and (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 references    . In a recent paper, Surana et al.  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 u1, u2, u3 as functions of x1, x2, x3 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 ( ). This is a necessary condition for an extremum of the functional I in (1).
3) The Euler’s equations derived from (using either fundamental lemma of the calculus of variations or the fourth basic lemma (   ) 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  -  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:
1) Kinematic assumptions describing the deformation of the plate cross-section and its middle surface. This leads to specific strain, hence stress fields and specific expressions for the moments obtained by integrating the moments of the stresses over the thickness of the plate or the shell.
2) The second significant difference arises due to specific chosen physics related to strain energy in the functional I.
1) Different choices in (1) - (2) lead to different mathematical models.
2) The approaches described above (energy functional or principle of virtual work) for deriving mathematical models for plates and shells can only be used for isothermal processes with reversible mechanical deformation and the energy equation is not part of the mathematical models. Thus, the mathematical models for the plate and shell theories requiring consideration of dissipation and memory mechanisms, nonlinear considerations and thermal effects cannot be derived using this approach.
3) Dissipation and memory mechanisms in the current plate and shell theories can only be addressed using phenomenological approaches employing 1D springs and 1D dashpots. Shortcomings of this approach and the difficulties associated with its extensions to and are well known .
4) The consideration of kinetic energy and potential energy of loads in the energy functional I is rather straight forward. However, the strain energy consideration requires rate of work conjugate pairs which are only known from the derivation of the energy equation (first law of thermodynamics). Since stress is a consequence of strain through material response, moment results due to curvature through material response, thus it is perhaps fitting to say the stress and the strain rate, moment and curvature rate are rate of work conjugate pairs. This in fact is the basis for strain energy considerations in the functional I. It is rather obvious that this approach relies heavily on clear and precise understanding of physics and may work well in simple deformation physics, but in more complex situations this approach may not be as straight forward and may even lead to erroneous choices.
Various Considerations, Scope and Outline of Work
In the following we outline the basic approach used and the scope of work presented in this paper.
1) In establishing thermodynamic consistency of the plate and shell mathematical models it suffices to use bending of plates only as shell formulations are their extensions. We consider x1x2 plane to be the middle surface of the plate of thickness h.
2) In the work presented here we consider small strain, small deformation and reversible mechanical deformation. The conservation and balance laws and the constitutive theory(ies) for such deformation are summarized in the “Appendix A” for classical as well as non-classical continuum mechanics.
3) The CPT and FSDT mathematical models used currently contain all features of the majority of the plate and shell mathematical models used currently, hence in this paper we only consider these two mathematical models as representative 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 (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  .
iii) Non-classical continuum mechanics using internal and Cosserat rotations  .
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.
2. Currently Used Mathematical Models for CPT and FSDT
For simplicity consider bending of a rectangular plate with thickness h with its middle surface in plane, hence normal to the middle plane of the plate/shell. Let , and be the displacements in directions of a point (see 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, refers to Cauchy stress on the i plane in the j direction. Fortunately, this notation is same in plate and shell mathematical models. Likewise 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
Figure 1. Schematic of the plate.
currently used plate and shell mathematical models using the same notations as used currently in the literature.
2.1. Classical Plate Theory (CPT) Based on Energy Functional (or the Principle of Virtual Work)
The classical plate theory is due to Kirchhoff hypothesis   and is based on the following kinematic assumptions for the displacements at an arbitrary point for an arbitrary value of time t.
Using (1) we can derive the linear strain measures ( symmetric part of the displacement gradient tensor):
Based on (2), the strain energy is only due to , ; , ; , as , and are zero, hence make no contribution to strain energy. The mathematical model based on (1) and (2) and the linear constitutive theory for Cauchy stress tensor (isotropic, homogeneous matter)
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 ( ). Using integration by parts all differentiation is transferred from and to their coefficients. After grouping the coefficients of and and performing integration with respect
to x3 using limits , the coefficients of and for arbitrary
time, length and the remaining spatial limits are set to zero. These are the Euler’s equations constituting the mathematical model for CPT based on (1) - (4):
The shear forces Q13 and Q23 are obviously zero. The complete mathematical model for CPT consists of (5), (6) and (7) in and u3 in which M11, M22 and M12 are functions of u3 and are defined by (9). N11, N22 and Q12 can be easily obtained using (8) and (3). For pure bending without inplane load and deformation in x1x2 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 u3 describing pure bending of plates.
in which and .
1) As is well known, in this mathematical model the transverse shear stresses and (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 u3.
3) We note that moment M11 due to is about x2 axis, hence based on standard continuum mechanics notation it should have been M12 (on face x1 about x2 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 notation 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.
2.2. 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 u1, u2 and u3 at an arbitrary location in the plate for an arbitrary time t  .
In which and are the rotations of x1x3 and x2x3 sections of the plate with respect to center plane. These rotations and are in fact rotations about x2 and x1 axes once again inconsistency of notations, but we continue with (11) in and as stated above. Using (11), we can derive the linear strain measures
Based on (12) only , do not contribute to strain energy. If we assume linear constitutive theory with plane stress behavior in the plane(x1x2 plane), then
The mathematical model based on (11) - (13) is also (as in the case of CPT) derived by considering energy functional I consisting of kinetic energy, strain energy and the potential energy of loads over the volume V of the plate and time [0, t]. Using Lagrangian description we can write the following for I.
We use (11) - (13) in (14) and expand each term in (14). For an extremum of I, we set first variation of I to zero ( ). Using integration by parts all differentiation from , , , and is transferred to their coefficients. After grouping the coefficients of , , , and and
performing integration with respect to x3 using limits the coefficients
of , , , and for arbitrary time and spatial limits are set to zero. These are the Euler’s equations constituting the mathematical model for FSDT based on (11) - (14):
Equations (15) constitute the complete mathematical model in and . as functions of and are defined by (16) - (18).
1) We note that the transverse shear forces Q13 and Q23 are not zero in this mathematical model.
2) This mathematical model when compared with CPT is expected to yield higher transverse deflection u3 due to shear deformation (additional strain energy due to transverse shear strains , and shear stresses and ).
3) We note from (17) that transverse shear forces Q13 and Q23 are not functions of x3 i.e., these are constant along the thickness of the plate where as their actual distributions are parabolic functions of x3. Thus, the strain energy due , and , in this model is not correct(as well known). Shear correction factors   are used to correct this situation. In the present work, we do not consider shear correction factors, instead we use this final mathematical model (15) - (18) that is consistent with the kinematic assumptions (11) when obtaining solutions for comparison with the new formulation.
4) In case of pure bending (in the absence of in plane (x1x2 plane) loads) we only need to consider equations three through five in (15), equations one and two in (17) and equations in (18). This mathematical model for pure bending consists of eight first order partial differential equations in eight dependent variables u3, Q13, Q23, M11, M12, M22, and . Q13, Q23, M11, M12 and M22 can be eliminated by substituting them in equations three, four and five of (16). This gives rise to three partial differential equations in u3, and , containing up to second order derivatives of u3, and .
3. 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 space and time t as needed in the derivations of the CPT and FSDT plate bending mathematical models.
3.1. Classical Continuum Mechanics
Complete details of the mathematical model consisting of balance of linear momenta, 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.
3.2. 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.
3.2.1. Non-Classical Continuum Mechanics with Internal Rotations
In this continuum theory the entire displacement gradient tensor, i.e. symmetric part (strain measures) as well as the antisymmetric part (measures of internal rotations) are considered in the derivation of the conservation and balance laws (see Surana et al.    for details), and we have the following. Balance of linear momenta yields the same equations as in CCM (Equation (A2), “Appendix A”). Following Section A.2 in “Appendix A”, we have the following from the remaining balance laws including the linear constitutive theories:
Linear constitutive theories for , and are given by
Equations (A2), (20), (22) and (23) constitute eighteen partial differential equations in eighteen dependent variables , , and , hence the mathematical model has closure. , and are material coefficients.
3.2.2. Non-Classical Continuum Mechanics with Internal and Cosserat Rotations
In this continuum theory, we also consider internal rotations due to antisymmetric part of the displacement gradient tensor in addition to strain measures (hence we consider the displacement gradient tensor in its entirety) as well as three Cosserat rotations ( , additional degrees of freedom) at a material point both acting about the axes of the same triad with its axes parallel to fixed x-frame. Thus, now we have total rotations , sum of and at each material point. In the following we present conservation and balance laws and the constitutive theories that consider this physics. The balance of linear momenta remains the same as in CCM (Equation (A2) in “Appendix A”) except that Cauchy stress tensor is not symmetric.
Following Section A3 in “Appendix A”, we have the following from the remaining balance laws and the linear constitutive theories:
Linear constitutive theories for , , and are given by     :
Equations (A2), (26), (28) - (30) are a system of twenty one partial differential equations in twenty one dependent variables: , , , and .
4. 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.
4.1. 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 and ensure thermodynamic consistency. Thus, use of these is justified for bending of plates and shells in . 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 , and are obtained as shown in (2) and we also find that the , and . The expressions for the components of the Cauchy stress tensor are established using linear constitutive theory as shown in (3). If we assume deformation in , then even though . On the other hand if we assume plane stress, then and the only nonzero components of stress tensor are , and . In the present work, we consider plate deformation in , hence we have . In addition to (3), we have in which .
We substitute from (2) and (3) in the balance of linear momenta (A2), and consider integral over the volume of the plate. In which A is the area of the middle plane ( ) of the plate. We integrate with respect to . The resulting integral is valid over arbitrary area A, hence the integrand must be zero, which gives us the following equations:
Equations (32) are three partial differential equations in and u3.
1) This mathematical model (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.
2) If we only consider plane stress deformation in the plane of the plate as considered in almost all currently used theories, then . 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 and and hence the corresponding transverse shear forces Q23 and Q31 are zero (as also in the case of currently used models) for CPT. This is obviously non physical.
4) We note that in Equation (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 , but these are altered and damaged due to kinematic assumptions, finally resulting in a mathematical model that is not valid for the physics under consideration.
4.2. 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 into symmetric and skew symmetric components.
in which components of 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 can now be defined using (36)
where are internal rotations      at a material point about the axis of a triad with axes parallel to x-frame. The positive rotations in (35) are the counterclockwise sense. From (37) we note that kinematic assumptions in (1) employs internal rotations which are not supported by CCM, hence the motivation for considering NCCM with internal rotations for the derivations considered in this section. Details of strain tensor components remain the same as defined by (2). The stress tensor is decomposed into symmetric and skew symmetric tensors.
The components 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 , we obtain the following:
and are defined later.
The internal rotation gradient tensor
is decomposed into symmetric and skew symmetric tensors.
From balance of angular momenta we have
Balance of moment of moments show that  
The constitutive theories for and are     
The constitutive theory for , and are obtained using (46) and are same as those in (3). if we consider the deformation of plate bending in (discussed earlier). Using (47) and (42), explicit expressions for the constitutive theory for can be obtained in terms of the symmetric part of the rotation gradient tensor.
Using these expressions for in BAM (44) and dividing them by 2 and then integrating with respect to over we obtain
We note that and terms in (39) are zeros as and .
This mathematical model consists of BLM (39) (3 equations), BAM (44) (3 equations), constitutive theories for (or , Equations (40) and (49)) (4 equations) and constitutive theories for ( , Equations (47)) (5 equations). These are a total of fifteen partial differential equations in fifteen variables ; , hence the mathematical model has closure.
1) Obviously, the mathematical model used currently consisting of (5) - (9) is not the same as derived here using NCCM consisting of Equations (39), (40), (48) - (50).
2) The balance of linear momenta (39) contain usual terms related to the gradients of the stress tensor .
3) and are zero but and are nonzero due to nonzero and . Existence of or and or is due to Cauchy moment tensor (necessitated due to the presence of internal rotations).
4) The moment tensor is symmetric due to balance of moment of moments balance law   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 in the present derivations is Cauchy moment tensor. This is not the case for , and obtained in the currently used mathematical models by integrating moments of , and in the direction for .
7) The mathematical model derived here is consistent with the conservation and balance laws of NCCM and incorporates the kinematic assumptions of CPT, but is obviously different than the mathematical model used currently for CPT. Thus, we conclude the CPT mathematical model used currently is thermodynamically inconsistent based on NCCM.
8) Since the kinematic assumptions (1) in CPT does not contain unknown rotations (Cosserat rotations), there is no motivation or need for undertaking the derivation of CPT mathematical model based on NCCM that considers both internal and Cosserat rotations.
4.3. 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 direction between −h/2 to h/2 we obtain
We note that and as and are functions of and t only. is given in Equation (40).
Forces , , , and are defined by Equation (16).
1) This mathematical model does not have closure. We have Equations (52) (3 equations) and Equations (16) (5 equations), a total of eight equations in ten dependent variables , , ; , , , , ; and . Hence, this mathematical model is not valid.
2) Equations four and five in (15) are not possible in the derivation considered here as the time derivatives of and do not exist in this derivation. These only appear in the derivation based on energy methods or in the methods based on principle of virtual work.
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.
4.4. FSDT Model Derivation Using NCCM Based on Internal and Cosserat Rotations
The motivation for considering NCCM based on internal and Cosserat rotations in deriving the mathematical model becomes clear if we decompose displacement gradient into symmetric and antisymmetric tensors (as in Section 4.2) or alternatively expand .
Using and the kinematic relations (11) we can obtain
The total rotations , and are the rotations at a material point about the axes of a triad with its axes parallel to x-frame. These rotations consist of internal rotations as well as the unknown Cosserat rotations and , both of which are not supported in CCM. This suggests that we should perhaps undertake the derivations of FSDT plate theory using conservation and balance laws and the constitutive theories of NCCM that considers internal and Cosserat rotations (Section A.3 in “Appendix A”).
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 , then coefficients are modified and are due to linear constitutive theory in giving rise to and
The stress tensor is decomposed into symmetric and skew symmetric tensors and :
Definitions of the components of the is same as in (13). Substituting total stress tensor into the balance of linear momenta and integrating with respect to we obtain (39) in which
and appearing in (39) are defined by (67) - (71).
The gradient of total rotation is defined by
Components of can be obtained using definition of in (53) and using (57):
Decomposition of into symmetric ( ) and skew symmetric ( ) tensors gives
The nonzero components of the symmetric and the skew symmetric tensors in (59) can be easily obtained using (58).
From the balance of angular momenta    we have
From the balance of moment of moments balance law  , we can show that
The linear constitutive theories for , and are given by     
where , are Lame’s constants and , and are material coefficients related to constitutive theory for non-classical physics due to internal and Cosserat rotations. We integrate (60), (62), (63) and (64) with respect to to obtain the forces per unit length used in the balance of linear momenta (39). First, from (60):
We note that , , and due to , , are already defined in (40).
We substitute constitutive theory for , and in (69).
From (63) after integration of we obtain
Likewise integrating (64) with respect to between the limits we obtain
We finally rewrite the balance of linear momenta using symmetric and antisymmetric decomposition of shear forces.
The complete mathematical model consists of twenty equations. BLM (3 equations), BAM (3 equations), constitutive theory for defining , , , , , (through (40) and (69)) (6 equations), constitutive theory for defining , , (3 equations), constitutive theory for (5 equations) in dependent variables: , (3); (5); , , , , , , , , (9), and (2), a total of 19.
1) This mathematical model does not have closure as we have twenty equations but only nineteen dependent variable. Thus, this mathematical model derived based on the kinematic assumptions of FSDT with conservation and the balance laws of NCCM incorporating internal and Cosserat rotations is invalid.
2) As in case of the derivation in Section 4.3, here also we observe that the time derivative of , and present in the currently used FSDT model derived based on the energy method (or the principle of virtual work) does not exist in the derivation presented here using NCCM.
3) This mathematical model is not only different from the currently used FSDT mathematical model, but is also invalid due to lack of closure.
5. General Remarks
1) In Section 4 we have shown that currently used mathematical models for CPT and FSDT cannot be derived using the conservation and balance laws of classical as well as the non-classical continuum mechanics based on internal and/or Cosserat rotations.
2) Currently used mathematical model for CPT requires use of internal rotation(s) and the currently used FSDT model requires internal as well as Cosserat rotations. Furthermore, all currently used mathematical models for plate bending (and their extensions to shells) either require use of internal rotation(s) or internal and Cosserat rotation(s). Thus, use of CPT and FSDT as typical representative mathematical model for investigating thermodynamic consistency of all currently used plate/shell mathematical models is justified. Hence, the remarks and the inferences presented for CPT and FSDT mathematical models hold in general for all currently used plate/shell mathematical models.
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 , the phenomenologically incorporated dissipation and memory mechanism are only possible in and their extensions to and (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.
6. 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 (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 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 , and . 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  . 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 , 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 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 , and with variable choice of the orders of the approximation space to ensure desired global differentiability of the approximations. In the local approximation, displacements , and can be approximated by polynomials of degrees , and as well as higher orders of differentiability. Choices of p-levels , and define approximations in the plane of the plate or the shell element as well as deformation of its cross section(s). By choosing , and 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 , and 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 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  . For stationary processes i.e., boundary value problems and for linear elastic reversible deformation of solid continua, GM/WF yields unconditionally stable computational processes  . In all other cases, least squares method based on residual functional can be shown to obtain unconditionally stable computations  .
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.
8) It is important to point out that the formulation of p-version hierarchical plate/shell finite element formulation (for reversible deformation without
Figure 2. Curved shell element geometry, mappings and nodal configurations. (a) A 3D solid element with nine node configuration on its opposite faces (top and bottom); (b) A 3D curved shell (or plate) element with nodes on the middle surface and nodal vector ends defining top and bottom faces; (c) Node i on the middle surface of the shell and vector defining the top and bottom surfaces at node i in space ( ); (d) Map of the element of figure - (b) in the natural coordinate space , and ; (e) Nodal configuration in direction; (f) Map of node i and vector .
temperature effects and linear constitutive theory) was originally derived by Surana et al.   using principal of virtual work. Exactly similar derivation as in ref   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  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   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.
6.1. 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   - . Based on the conservation and balance laws (more specifically second law of thermodynamics), the choice of , , and as constitutive variables at the onset is rather straight forward. Using the conjugate pairs, and and thermoelastic deformation, choice of , as argument tensors of and and as argument tensors of is appropriate. Using the principle of equipresence, we can choose , and as argument tensors for , . Thus, we have
6.1.1. 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 .
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 , the following must hold
and SLT (78) reduces to
with (79) - (82), SLT (78) is satisfied.
a) Equation (79) can be used to derive constitutive theory for Cauchy stress tensor .
b) Equation (80) implies that is not a constitutive variable as it is defined by .
c) Equation (82) implies that is not a function of , thus .
In the following we derive constitutive theory for using
Since, the constitutive theories must be frame independent, we must consider ; being principle invariants of . Thus (83) becomes
Following the details of similar derivations in reference , we can obtain the following from (85) after using the Hamilton-Cayley theorem.
Material coefficients are determined by expanding ; in Taylor series in , , and about a known configuration and only retaining up to linear terms in the invariants and temperature (for simplicity).
We substitute (88) in (86) and group the terms and coefficients and if we define
Then (86) can be written as :
in which; ; and are material coefficients that are functions of invariants ; i.e., , , and in the known configuration .
1) This constitutive theory for contains ( ) fourteen material coefficients and contains up to fifth degree terms of the components of strain tensor , but is linear in temperature .
2) is the initial stress field.
3) A linear constitutive theory in which the products of , , are neglected and only up to linear terms in and are retained is given by
using the notation (Lame’s constants) and realizing that and we can write (91) as follows
this of course is generalized Hooke’s law.
b) Constitutive theory for can also be derived using and the representation theorem   - , this leads to
in which , and are combined generators of the argument tensors and that are symmetric tensors of rank two. Since (93) is same as (86), the rest of the derivation remains same as in case of the derivation using given above.
6.1.2. Constitutive Theory of q
We begin with based on entropy inequality and use representation theorem   -  which gives us the following constitutive theory for
This constitutive theory is based on integrity, hence uses complete basis. From (94), we can also derive a linear constitutive theory
where are material coefficients defined in a known configuration . These can be functions of (invariant of ) and .
7. Final Mathematical Model (Thermoelastic)
The entropy inequality is identically satisfied with the choices we have made in deriving the constitutive theories.
1) Use of this mathematical model (96) - (100) in will result in a comprehensive nonlinear plate and shell formulation that describes evolution in which [(a)]
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   based on principle of virtual work, the resulting Euler’s equations representing the associated mathematical model are
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  , then the mathematical model (96) - (100) obviously reduces to (101) and (102).
3) We clearly see that enhancement of mathematical model used in reference   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 (simplified form of (96) - (100)).
8. Finite Element Formulation
8.1. Complete Mathematical Model
Expanded form of the mathematical model in ((101), (102)) are given in the following:
in which is the closed spatial domain in . is the closed boundary of . The linear constitutive theory for Cauchy stress tensor and in terms of linear strain is given in the following (using Voigt’s notation)
The nonzero elements of (6 × 6) [D] material coefficient matrix are given by:
The boundary conditions (BCs) are defined by:
where, , a closed boundary of , , and are direction cosines of a unit exterior normal to the boundary . E and are modulus of elasticity and Poisson’s ratio. Equations (103) and (104), nine equations in nine unknowns and constitute the complete mathematical model for the plate or the shell in . are known displacement boundary conditions on and , and are known tractions on the boundary .
8.2. 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 , and . One could easily show that the differential operator and its adjoint in this system of differential equations are same, hence for this system of Equations (103) and (104) in whatever form, GM/WF would yield an integral form that is variationally consistent (VC)  , thus the resulting computational processes are ensured to be unconditionally stable for all choices of computational parameters (h, p, k) and the physical parameters in the mathematical model. Construction of the finite element formulation using GM/WF is facilitated using the mathematical model in the form of (103) and (104). Let be discretization of , then in which is a typical element “e”. Let be test functions such that
We construct integral form over using (103) and (110) based on fundamental lemma of the calculus of variations   .
We substitute from (103) in (111) for , and and transfer one order of differentiation from the stress derivative terms to the test functions.
Let , and be approximations of and over and , and be the approximations of and over an element “e” with domain such that where are the approximations of over .
Consider ; over
Transfer one order of differentiation from the stresses to the test functions , and in (113) using integration by parts (IBP) to obtain
From (114), are primary variables (PVs) and the coefficients of are secondary variables (SVs). Boundary conditions (108) are essential boundary conditions (EBCs) and those defined by (109) are natural boundary conditions (NBCs). Let
be equal order, equal degree local approximations for and in which are approximation functions and , and are nodal degrees of freedom for , and . From (115), we obtain
First we define secondary variables , and as
We substitute strains (106) in (104), then (104) in (114) to obtain
: a vector of secondary variables at the nodes.
The matrix is defined by
is the element stiffness matrix, 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 usual procedure  . We note that , 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 and we need to determine the local approximation functions in in (115).
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(a). We connect all nine pairing nodes between the bottom and top surfaces of the element and define vectors connecting nodes and . The surface containing the middle points of vectors are the shell element nodes and are assumed to represent the middle surface of the plate/shell element. The surfaces containing the bottom and top ends of these vectors define the bottom and top surfaces of the shell element (Figure 2(b)). Details at a typical node i of the plate/shell element are shown in (Figure 2(c)). The coordinates of nodes (Figure 2(b)) and the vectors completely define the geometry of shell element in . 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 is mapped into , being the map of the element in the two unit cube. In the map of the element, the element nodes are located in plane at . The natural coordinates and define bottom and top faces of the plate/shell element. Thus, each is mapped into a two unit length in the direction (perpendicular to plane at ). Let , , ; be the nodal coordinates of the shell element nodes located in plane. Then, the coordinates of an arbitrary point say for arbitrary but are defined by
in which ; for . The coordinates of a point at an arbitrary location with respect to the middle plane ( ) are given by
Thus, the coordinates at any arbitrary point in the plate/shell element are given by ; . Using (123) and (124) we obtain the following for the coordinates at an arbitrary point at a location in the natural coordinate map .
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 for , and . Such approximations are in a scalar product space ; , with p-levels , and . Extensions of the approximations to higher classes defining local approximations of higher order global differentiability can be considered based on ref    - . To illustrate the details, we consider to be the dependent variable for which we wish to establish local approximation functions such that
The p-version local approximation for in the direction for the three node configuration (with nodes 1, 2 and 3 located at , and can be written as  .
Similarly in the direction for the three node configuration (with nodes 1, 2 and 3 located at , and ) we can write the following for p-version hierarchical local approximation.
In the direction, the shell element nodes are located at (Figure 2(b), Figure 2(d)), hence the Lagrange (or others) interpolations (local approximations) in the direction for , etc. corresponding to 2, 3, 4, …, etc. nodal configurations in direction respectively need to be reduced down to a single node located at .
We consider a map of at each node in coordinate system pointing in the same direction as in which is mapped into at each node (say i) (Figure 2(e)). This mapping is given by
in which , length of at node i of the shell element. Derivations of the one dimensional p-version hierarchical functions in or direction is given in the following. Consider a typical plate/shell element node i ( ) shown in Figure 2(c), Figure 2(d), Figure 2(f). For p-level of one ( ) in the direction the Lagrange interpolation for the two node configuration of Figure 2(e) can be written as
Substituting from (130) in (131)
using (134) and (135) in (131)
Likewise for , for the three node configuration in Figure 1, we can write (Lagrange interpolation)
Substituting for from (130)
differentiating with respect to twice and evaluating the derivatives at and substituting these back in (138) and then substituting for from (130) we obtain the following for the shell node i
Following the same procedure, we can derive the following for a node i located at for p-level of .
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.
In the one dimensional approximations (127), (128) and (140) we separate the the approximation functions and the nodal variable operators (from the dofs)  .
In direction: for nodes 1, 2 and 3
Nodal variable operators:
In direction: for nodes 1, 2 and 3
Nodal variable operators:
In direction: for node i at or
Nodal variable operators:
By taking tensor products of the 1D nodal variable operators in and and letting them act on we obtain the degrees of freedom for or over (hence over ) and the corresponding approximation functions are obtained by taking the tensor product of the 1D approximation functions in and and we can write
in which are the degrees of freedom. Using (141) for and we can write the following for local approximation for and of and (assuming and to be same for and ).
In which are the degrees of freedom for and for an element e.
From the geometry defined by (125), we define the Jacobian of mapping between and master element (using Murnaghan’s notation)  .
The volumes in and frames and the derivatives of the approximation functions are related through
Using (145), matrix in (122) is completely defined. Coefficients of are calculated using Gauss quadrature with element map . Details can be found in references  . Assembly of the element equations and their solutions follows standard procedure  .
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 and . While p-levels and primarily control the deformation in the plane, controls transverse deformation.
3) With local approximation described here, the integral over are in Lebesgue sense. Smoothness of the analytical solution considered in the model problems ensures convergence of solution to class 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   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   (as ) results in more comprehensive mathematical model as presented here.
9. Solutions of Model Problems
We consider a square plate ( ) simply supported on all four boundaries. Midplane of the plate lies in plane and the origin of fixed x-frame is at the bottom left corner. Plate thickness is “h”. is normal to - plane pointing upwards. Top surface of the plate ( ) is subjected to uniformly distributed load acting in the negative direction. A schematic of the plate is shown in Figure 3. We nondimensionalize density , modulus and distance (or displacement) by using reference density , reference modulus and reference length . Reference velocity and reference time are and . Schematic of the plate, coordinate system, material properties, reference quantities and the dimensionless quantities are shown in Figure 3.
The values of the distributed load q are chosen such that it produces a deflection at point A ( ) of 0.1 for all thicknesses using CPT. This gives us the following values of q for thicknesses of and 10.0 (Table 1).
Figure 3. Model problem; schematic, discretization, material properties and boundary conditions.
Table 1. Distributed load “q” values based on CPT model.
Material Properties Reference Quantities Dimensionless Quantities
In the numerical studies we choose hence . We consider a 4 × 4 uniform discretization for the entire plate. We consider four different thicknesses: and 10.0. Due to the choice of , the dimensionless distances, coordinates, thicknesses, deflections etc., are same as those with dimensions. In all studies, we keep the discretization fixed (4 × 4 uniform) and increase p-levels in the and direction ( and ) to obtain converged displacements, strains and stresses. The local approximations are considered in spaces with , hence solutions of class . For this choice of local approximations, the integrals over (discretization) are in Lebesgue sense. Due to the fact that solutions of the model problems are smooth, with p-refinement we expect solutions of class to converge to in the weak sense. Solutions of the mathematical models for CPT and FSDT given in the following are also computed using finite element method based on residual functional (least squares method) with p-version hierarchical local approximations. For all four thicknesses of the plate a (2 × 2) discretization with p-levels of 4, 4 or 5, 5 was sufficient to obtain converged displacements for CPT and FSDT. These solutions of displacement for CPT and FSDT are compared with the new formulation for all four thicknesses.
(a) (b) (c) (d)
Figure 4. Displacement of the centerline ( ) versus axial distance for . (a) Displacement versus ( ); (b) Displacement versus ( ); (c) Displacement versus ( ); (d) Displacement versus ( ).
Figure 5. Mating element edges at .
(a) (b) (c) (d)
Figure 6. Shear stress versus distance at for . (a) Shear stress versus at ( ); (b) Shear stress versus at ( ); (c) Shear stress versus at ( ); (d) Shear stress versus at ( ).
(a) (b) (a) (b)
Figure 7. Normal stress versus distance at for . (a) Normal stress versus at ( ); (b) Normal stress versus at ( ); (c) Normal stress versus at ( ); (d) Normal stress versus at ( ).
Discussion of Results
Figures 4(a)-(d) show plots of transverse displacement of the middle plane of the plate as a function of at for and 10.0. In the case of new formulation, and yield virtually the same displacement (shown in figures). In case of , the plate is thin, hence CPT and FSDT are expected to produce almost identical results (Figure 4(a)). The new formulation is 3D elasticity formulation, hence contains more comprehensive and complete deformation physics compared to CPT and FSDT. Plot of versus for the new formulation are in quite good agreement with those from CPT and FSDT. For , CPT under estimates (Figure 4(b)) due to lack of shear deformation. For this thickness ( ), FSDT results are quite reasonable as the shear deformation contribution is not very significant. The converged versus from the new formulation and those from FSDT are in good agreement (Figure 4(b)). For and , the shear deformation is progressively more dominant. Figure 4(c) and Figure 4(d) show versus for CPT, FSDT and the new formulation. In both figures, we observe higher values of versus for FSDT compared to CPT. The new formulation yields higher values of compared to FSDT as well as CPT and the difference between CPT and new formulation as well as the difference between FSDT and the new formulation increases for compared to . These results from the mathematical models are in accordance with the physics of deformation.
Next we consider the behavior of at as a function of coordinate (obtained using the new formulation). At , four vertical edges are coincident. Since the local approximations are of class , we expect discontinuity of (function of displacement gradients) at lower p-level. However, with progressively increasing p-levels we expect convergence of (weak convergence of solutions to the class ).
Figures 6(a)-(d) show plots of versus at for
and 10.0. From Figure 6(a) for , we note that at and , there is a discontinuity of . We have two graphs of versus even though at , there are four element edges coincident. Figure 5 shows locations with four element edges marked that are coincident. Let be shear stress along
the edges . Since ,
due to discontinuity of even though is continuous. Similarly , however and . Thus, at graphs of and versus are coincident. Likewise graphs of and versus are coincident as well but have different values than those of and . Thus, in Figure 6(a) we only observe two sets of graphs. The discontinuity of is due to discontinuity of along the edge, a consequence of local approximation for displacements at and . At , , the discontinuity of diminishes and at , and 11-14 versus is almost converged (in weak sense) and has unique values from all the edges. Graphs of versus at for and 10.0 shown in Figures 6(b)-(d) exhibit exactly similar behavior at lower p-levels. We observe discontinuity of but with progressively increasing p-levels, the discontinuity of diminishes and for p-levels in the range of 5 to 9, we obtain converged solutions for . Higher thicknesses of plate naturally require higher p-levels. versus at behaves in exactly similar fashion; graphs and the details are omitted for the sake of brevity.
Figures 7(a)-(d) show plots of versus at for and 10.0. On the top surface of the plate ( ) a uniform pressure load of is applied and on the bottom surface of the plate ( ), . This holds regardless of the thickness of the plate. For all thickness values, p-levels and are needed to obtain converged values of , but no difficulties are encountered in doing so. In case of CPT, FSDT and HSDT plate mathematical models, . This is obviously non physical especially for thick plates. In case of the new formulation presented here is calculated accurately regardless of the plate thickness.
10. 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 due to the absence of the energy equation and the entropy inequality. This approach cannot be extended to and . 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 : 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 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 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 GM/WF in which the integral form is variationally consistent  , hence the resulting computations are unconditionally stable. The local approximation for displacements in constructed in space with i.e., solutions of class but with higher p-levels (p) in all three directions . In this approach, local approximations with different p-levels choices result in the corresponding kinematic descriptions, hence these vary depending on the choice of the p-levels. The formulation presented here is based on CCM in . It addresses physics of deformation of extremely thin as well as thick plates/shells as demonstrated in the model problem studies. Extension of this formulation for plates/shells with dissipation and memory mechanism is consistent and systematic through entropy inequality and free of any phenomenological considerations.
It is worth remarking that the new formulation for plates/shells is a truly a formulation in . in all plate/shell mathematical models is generally zero, but this is not physical. The new formulation presented here simulates for thin as well as thick plates without any difficulty as shown in the model problems. We remark that even though we have considered a flat square plate as a model problem, the formulation is naturally valid for curved shells with curved geometry. This is obvious from the geometric description of the shell and the use of conservation and balance laws in constituting the mathematical model. Lastly, we point out that this single formulation presented in this paper addresses all plate/shells deformation physics regardless of their thickness or geometry. Studies for curved shells will be presented in a subsequent paper.
First author is grateful for his endowed professorship and the department of mechanical engineering of the University of Kansas for providing financial support to the second author. The computational facilities provided by the Computational Mechanics Laboratory of the mechanical engineering departments are also acknowledged.
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 .
In which is the mass density in the reference configuration, , and are body force per unit mass in and directions, is Cauchy stress tensor, is linear strain tensor, e is specific internal energy, is heat flux vector, is temperature gradient vector, is Helmholtz free energy density, is entropy density and is temperature.
From the conjugate pairs and , 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 .
A.2. Non-Classical Continuum Theory Incorporating Internal Rotations
Displacement gradient tensor and its decomposition into symmetric and antisymmetric tensors can be written as
in which and are symmetric and antisymmetric tensor, thus and are strain and rotation tensors. contains rotations , and about , and axes. Alternatively
with this definitions of , and in (9), we have
all others are zero. The rotations , and are about the axes of the triad located at a material point. If we define the rotations as a vector , , the gradient of i.e. and its decomposition into symmetric and skew-symmetric tensor and can be written as
We have the following conservation and balance laws   
Additionally we have used decomposition of Cauchy stress tensor into symmetric tensor and antisymmetric tensor
The Cauchy moment tensor is symmetric due to balance of moment of moments balance law, an additional balance law needed in non-classical continuum theories    due to new physics associated with rotations. From the conjugate pairs in the entropy inequality, at the very least, the following must hold (thermoelastic solid continua)
Choice of the argument tensors for and can be based on the principle of equipresence . The derivation of the constitutive theory for is same as for in case of CCM (section). A linear constitutive theory for can be written as
where is the material coefficient related to the constitutive theory for the Cauchy moment tensor.
A.3. Non-Classical Continuum Theory Incorporating Internal Rotations and Cosserat Rotations
Let be external or Cosserat rotations (unknown) about the axes of the same triad at a material point about which internal rotations act, then the total rotations are given by
Gradient of , and its decomposition into symmetric and antisymmetric tensors gives
The CM, BLM, BAM and BMM balance laws in this case are same as in Section A.2. The FLT and the SLT   are given by
The decomposition of (Section A.2) is used here as well. The Cauchy moment tensor is symmetric in this case also (balance of moment of moments balance law  ). From the conjugate pairs in the entropy inequality (28), at the very least the following must hold (thermoelastic solid continua)
Argument tensors for and at this stage can be established using principle of equipresence. Constitutive theories for and in the absence of reduce to  
 Han, S.M., Benaroya, H. and Wei, T. (1999) Dynamics of Transversely Vibrating Beams Using Four Engineering Theories. Final Version, Academic Press, Cambridge.
 Kirchhoff, G. (1850) Uber das Gleichgewicht und die Bewegung einer elastischen Scheibe. Journal für die reine und angewandte Mathematik, ger, 40, 51-88.
 Surana, K.S., Mysore, D. and Reddy, J.N. (2019) Thermodynamic Consistency of Beam Theories in the Context of Classical and Non-Classical Continuum Mechanics and a Thermodynamically Consistent New Formulation. Continuum Mechanics and Thermodynamics, 31, 1283-1312.
 Surana, K.S., Joy, A.D. and Reddy, J.N. (2017) Non-Classical Continuum Theory for Solids Incorporating Internal Rotations and Rotations of Cosserat Theories. Continuum Mechanics and Thermodynamics, 29, 662-698.
 Surana, K.S., Mysore, D. and Reddy, J.N. (2019) Non-Classical Continuum Theories for Solid and Fluent Continua and Some Applications. International Journal of Smart and Nano Materials, 10, 28-89.
 Reddy, J.N. (2015) An Introduction to Nonlinear Finite Element Analysis: With Applications to Heat Transfer, Fluid Mechanics and Solid Mechanics. Oxford University Press, Oxford.
 Surana, K.S., Mohammadi, F. and Reddy, J.N. (2016) Ordered Rate Constitutive Theories for Non-Classical Internal Polar Thermoviscoelastic Solids without Memory. International Journal of Mathematics, Science, and Engineering Applications, 10, 99-131.
 Yang, F., Chong, A.C.M. and Lam, D.C.C. (2002) Couple Stress Based Strain Gradient Theory for Elasticity. International Journal of Solids and Structures, 39, 2731-2743.
 Surana, K.S., Shanbhag, R. and Reddy, J.N. (2018) Necessity of Law of Balance of Moment of Moments in Non-Classical Continuum Theories for Solid Continua. Meccanica, 53, 2939-2972.
 Surana, K.S. and Sorem, R.M. (1990) Higher Order Hierarchical Plates and Curved Shell Elements Based on p-Approximation in the Thickness Direction. Mathematical and Computer Modelling, 14, 849-854.
 Surana, K.S. and Sorem, R.M. (1991) Completely Hierarchical p-Version Curved Shell Element for Laminated Composite Plates and Shells. Computational Mechanics, 7, 237-251.
 Todd, J.A. (1948) Ternary Quadratic Types. Philosophical Transactions of the Royal Society of London. Series A: Mathematical and Physical Sciences, 241, 399-456.
 Smith, G.F. (1970) On a Fundamental Error in Two Paper of C. C. Wang, “On Representation for Isotropic Functions, Part I and Part II”. Archive for Rational Mechanics and Analysis, 36, 161-165.
 Smith, G.F. (1971) On Isotropic Functions of Symmetric Tensors, Skew-Symmetric Tensors and Vectors. International Journal of Engineering Science, 9, 899-916.
 Spencer, A.J.M. and Rivlin, R.S. (1959) The Theory of Matrix Polynomials and Its Application to the Mechanics of Isotropic Continua. Archive for Rational Mechanics and Analysis, 2, 309-336.
 Spencer, A.J.M. (1971) Chapter 3 Treatise on Continuum Physics, I. In: Erigen, A.C., Ed., Theory of Invariants, Academic Press, New York, 239-353.
 Zheng, Q.S. (1993) On the Representations for Isotropic Vector-Valued, Symmetric Tensor-Valued and Skew-Symmetric Tensor-Valued Functions. International Journal of Engineering Science, 31, 1013-1024.
 Zheng, Q.S. (1993) On Transversely Isotropic, Orthotropic and Relatively Isotropic Functions of Symmetric Tensors, Skew-Symmetric Tensors, and Vectors. International Journal of Engineering Science, 31, 1399-1453.
 Surana, K.S., Petti, S.R., Ahmadi, A.R. and Reddy, J.N. (2004) On p-Version Hierarchical Interpolation Functions for Higher-Order Continuity Finite Element Models. International Journal of Computational Engineering Science, 2, 653-673.
 Maduri, R.K., Surana, K.S., Romkes, A. and Reddy, J.N. (2009) Higher Order Global Differentiability Local Approximations for 2-D Distorted Triangular Elements. International Journal for Computational Methods in Engineering Science and Mechanics, 10, 20-26.
 Ahmadi, A., Surana, K.S., Maduri, R.K., Romkes, A. and Reddy, J.N. (2009) Higher Order Global Differentiability Local Approximations for 2-D Distorted Quadrilateral Elements. International Journal for Computational Methods in Engineering Science and Mechanics, 10, 1-19.
 Surana, K.S., Ahmadi, A.R. and Reddy, J.N. (2002) The k-Version of Finite Element Method for Self Adjoint Operators in BVP. International Journal of Computational Engineering Science, 2, 155-218.
 Surana, K.S., Ahmadi, A.R. and Reddy, J.N. (2003) The k-Version of Finite Element Method for Non-Self-Adjoint Operators in BVP. International Journal of Computational Engineering Science, 4, 737-812.
 Surana, K.S., Ahmadi, A.R. and Reddy, J.N. (2004) The k-Version of Finite Element Method for Nonlinear Operators in BVP. International Journal of Computational Engineering Science, 5, 133-207.