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

Karan S. Surana, Sri Sai Charan Mathi

Department of Mechanical Engineering, University of Kansas, Lawrence, Kansas.

DOI: 10.4236/ajcm.2020.102010                

Abstract

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.

Keywords

Plate and Shell Mathematical Models, Energy Functional, Thermodynamic Consistency, Classical Continuum Mechanics, Non-Classical Continuum Mechanics, Internal Rotations, Cosserat Rotations, Principle of Virtual Work

1. 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] [2] [3] [4] [5]. 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 [6] - [12]. A formal presentation of the plate bending theory now referred to as classical plate theory (CPT) is due to Kirchhoff (1850) 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 $u_1\left(x_1,x_2,x_3,t\right)$ and $u_2\left(x_1,x_2,x_3,t\right)$ at an arbitrary point using the rotations

$\frac{\partial u_3}{\partial x_1}$ and $\frac{\partial u_3}{\partial x_2}$ while displacement u₃ remains a function of x₁, x₂ 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 [13] - [21]. The first significant improvement was due

to what is now referred as FSDT in which $\frac{\partial u_3}{\partial x_1}$ and $\frac{\partial u_3}{\partial x_2}$ in the kinematic

assumptions for $u_1$ and $u_2$ are replaced by unknown rotations ${\varphi }_1$ and ${\varphi }_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 references [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₁, u₂, u₃ as functions of x₁, x₂, x₃ 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 ( $\delta I=0$ ). This is a necessary condition for an extremum of the functional I in (1).

3) The Euler’s equations derived from $\delta I=0$ (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:

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.

Remarks

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 ${ℝ}^2$ and ${ℝ}^3$ are well known [30].

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 x₁x₂ 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 ${ℝ}^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].

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 $a\times b\times h$ rectangular plate with thickness h with its middle surface in $x_1{x}_2$ plane, hence $x_3$ normal to the middle plane of the plate/shell. Let $u_1\left(x_i\mathrm{,}t\right)$, $u_2\left(x_i\mathrm{,}t\right)$ and $u_3\left(x_i\mathrm{,}t\right)$ be the displacements in $x_1,x_2,x_3$ directions of a point $x_i=\left(x_1,x_2,x_3\right)\in \left[\mathrm{0,}a\right]\times \left[\mathrm{0,}b\right]\times \left[-h/2\mathrm{,}h/2\right]$ (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, ${\sigma }_{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 $m_{ij}$ 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 [11] [12] and is based on the following kinematic assumptions for the displacements $u_1\mathrm{,}{u}_2\mathrm{,}{u}_3$ at an arbitrary point $x_1,x_2,x_3\in \left(\left[\mathrm{0,}a\right]\times \left[\mathrm{0,}b\right]\times \left[-h/2\mathrm{,}h/2\right]\right)$ for an arbitrary value of time t.

$\begin{array}{l}{u}_1\left(x_1,x_2,x_3,t\right)={\tilde{u}}_1\left(x_1,x_2,0,t\right)-x_3\frac{\partial u_3\left(x_1,x_2,0,t\right)}{\partial x_1}\\ u_2\left(x_1,x_2,x_3,t\right)={\tilde{u}}_2\left(x_1,x_2,0,t\right)-x_3\frac{\partial u_3\left(x_1,x_2,0,t\right)}{\partial x_2}\\ u_3=u_3\left(x_1,x_2,0,t\right)\end{array}$ (1)

Using (1) we can derive the linear strain measures ( ${\epsilon }_{ij}$ symmetric part of the displacement gradient tensor):

$\begin{array}{l}{\epsilon }_{11}=\frac{\partial {\tilde{u}}_1}{\partial x_1}-x_3\frac{{\partial }^2{u}_3}{\partial x_1^2};{\epsilon }_{22}=\frac{\partial {\tilde{u}}_2}{\partial x_2}-x_3\frac{{\partial }^2{u}_3}{\partial x_2^2};{\epsilon }_{33}=0\\ {\epsilon }_{23}={\epsilon }_{31}=0;{\epsilon }_{12}=\frac{1}{2}\left(\frac{\partial {\tilde{u}}_1}{\partial x_2}+\frac{\partial {\tilde{u}}_2}{\partial x_1}\right)-x_3\frac{{\partial }^2{u}_3}{\partial x_2\partial x_1}\end{array}$ (2)

Based on (2), the strain energy is only due to ${\sigma }_{11}$, ${\epsilon }_{11}$; ${\sigma }_{22}$, ${\epsilon }_{22}$; ${\sigma }_{12}$, ${\epsilon }_{12}$ as ${\epsilon }_{33}$, ${\epsilon }_{23}$ and ${\epsilon }_{31}$ 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)

$\begin{array}{l}{\sigma }_{11}=D_{11}{\epsilon }_{11}+D_{12}{\epsilon }_{22};{\sigma }_{22}=D_{21}{\epsilon }_{11}+D_{22}{\epsilon }_{23};{\sigma }_{12}=D_{66}\left(2{\epsilon }_{12}\right)\\ D_{11}=\frac{E}{1-{\nu }^2}=D_{22};D_{12}=\nu D_{11}=D_{21};D_{66}=G=\frac{E}{2\left(1+\nu \right)}\end{array}$ (3)

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:

$\begin{array}{l}I={\displaystyle \underset{0}{\overset{t}{\int }}}{\displaystyle \underset{V}{\int }}(\frac{1}{2}{\rho }_0\left({\left(\frac{\partial u_1}{\partial t}\right)}^2+{\left(\frac{\partial u_2}{\partial t}\right)}^2+{\left(\frac{\partial u_3}{\partial t}\right)}^2\right)\\ \begin{array}{c}\\ \end{array}-{\sigma }_{11}{\epsilon }_{11}-{\sigma }_{22}{\epsilon }_{22}-{\sigma }_{12}\left(2{\epsilon }_{12}\right)+{\rho }_0{u}_i{F}_i^b)\text{d}V\text{d}t\end{array}$ (4)

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 ( $\delta I=0$ ). Using integration by parts all differentiation is transferred from $\delta {\tilde{u}}_1,\delta {\tilde{u}}_2$ and $\delta u_3$ to their coefficients. After grouping the coefficients of $\delta {\tilde{u}}_1,\delta {\tilde{u}}_2$ and $\delta u_3$ and performing integration with respect

to x₃ using limits $\left[\frac{-h}{2},\frac{h}{2}\right]$, the coefficients of $\delta {\tilde{u}}_1,\delta {\tilde{u}}_2$ and $\delta u_3$ 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):

$h{\rho }_0\frac{{\partial }^2{\tilde{u}}_1}{\partial t^2}-\frac{\partial N_{11}}{\partial x_1}-\frac{\partial Q_{12}}{\partial x_2}-h{\rho }_0{F}_1^b=0$ (5)

$h{\rho }_0\frac{{\partial }^2{\tilde{u}}_2}{\partial t^2}-\frac{\partial Q_{12}}{\partial x_1}-\frac{\partial N_{22}}{\partial x_2}-h{\rho }_0{F}_2^b=0$ (6)

$\begin{array}{c}h{\rho }_0\frac{{\partial }^2{u}_3}{\partial t^2}-\frac{h^3}{12}\frac{{\partial }^4{u}_3}{\partial t^2\partial x_1^2}-\frac{h^3}{12}\frac{{\partial }^4{u}_3}{\partial t^2\partial x_2^2}-\frac{{\partial }^2{M}_{11}}{\partial x_1^2}-2\frac{\partial M_{12}}{\partial x_1\partial x_2}-\frac{{\partial }^2{M}_{22}}{\partial x_2^2}-h{\rho }_0{F}_3^b=0\end{array}$ (7)

$\forall \left(x_i\mathrm{,}t\right)\in {\Omega }_{xt}={\Omega }_x\times {\Omega }_t$

in which

$N_{11}={\displaystyle \underset{-h/2}{\overset{h/2}{\int }}}{\sigma }_{11}\text{d}{x}_3\mathrm{<mo>\; ;\; </mo>}{N}_{22}={\displaystyle \underset{-h/2}{\overset{h/2}{\int }}}{\sigma }_{22}\text{d}{x}_3\mathrm{<mo>\; ;\; </mo>}{Q}_{12}={\displaystyle \underset{-h/2}{\overset{h/2}{\int }}}{\sigma }_{12}\text{d}{x}_3$ (8)

$\begin{array}{l}{M}_{11}={\displaystyle \underset{-h/2}{\overset{h/2}{\int }}}{x}_3{\sigma }_{11}\text{d}{x}_3=\frac{-h^3}{12}\left(D_{11}\frac{{\partial }^2{u}_3}{\partial x_1^2}+D_{12}\frac{{\partial }^2{u}_3}{\partial x_2^2}\right)\\ M_{22}={\displaystyle \underset{-h/2}{\overset{h/2}{\int }}}{x}_3{\sigma }_{22}\text{d}{x}_3=\frac{-h^3}{12}\left(D_{21}\frac{{\partial }^2{u}_3}{\partial x_1^2}+D_{22}\frac{{\partial }^2{u}_3}{\partial x_2^2}\right)\\ M_{12}={\displaystyle \underset{-h/2}{\overset{h/2}{\int }}}{x}_3{\sigma }_{11}\text{d}{x}_3=\frac{-h^3}{12}\left(D_{66}\frac{{\partial }^2{u}_3}{\partial x_2\partial x_1}\right)\end{array}$ (9)

The shear forces Q₁₃ and Q₂₃ are obviously zero. The complete mathematical model for CPT consists of (5), (6) and (7) in ${\tilde{u}}_1,{\tilde{u}}_2$ and u₃ in which M₁₁, M₂₂ and M₁₂ are functions of u₃ and are defined by (9). N₁₁, N₂₂ and Q₁₂ can be easily obtained using (8) and (3). For pure bending without inplane load and deformation in x₁x₂ 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₃ describing pure bending of plates.

$\frac{{\partial }^4{u}_3}{\partial x_1^4}+\frac{{\partial }^4{u}_3}{\partial x_2^2\partial x_1^2}+\frac{{\partial }^4{u}_3}{\partial x_2^4}=\frac{q}{D}$ (10)

in which $q={\rho }_{0}h{F}_3^b$ and $D=\frac{E{h}^3}{12\left(1-{\nu }^2\right)}$.

Remarks.

1) As is well known, in this mathematical model the transverse shear stresses ${\sigma }_{23}$ and ${\sigma }_{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) We note that moment M₁₁ due to ${\sigma }_{11}$ is about x₂ axis, hence based on standard continuum mechanics notation it should have been M₁₂ (on face x₁ about x₂ 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 u₁, u₂ and u₃ at an arbitrary location $x_1,x_2,x_3$ in the plate for an arbitrary time t [14] [18].

$\begin{array}{l}{u}_1\left(x_1,x_2,x_3,t\right)={\tilde{u}}_1\left(x_1,x_2,0,t\right)-x_3{\varphi }_1\left(x_1,x_2,0,t\right)\\ u_2\left(x_1,x_2,x_3,t\right)={\tilde{u}}_2\left(x_1,x_2,0,t\right)-x_3{\varphi }_2\left(x_1,x_2,0,t\right)\\ u_3=u_3\left(x_1,x_2,0,t\right)\end{array}$ (11)

In which ${\varphi }_1$ and ${\varphi }_2$ are the rotations of x₁x₃ and x₂x₃ sections of the plate with respect to center plane. These rotations ${\varphi }_1$ and ${\varphi }_2$ are in fact rotations about x₂ and x₁ axes once again inconsistency of notations, but we continue with (11) in ${\varphi }_1$ and ${\varphi }_2$ as stated above. Using (11), we can derive the linear strain measures

$\begin{array}{l}{\epsilon }_{11}=\frac{\partial {\tilde{u}}_1}{\partial x_1}-x_3\frac{\partial {\varphi }_1}{\partial x_1};{\epsilon }_{22}=\frac{\partial {\tilde{u}}_2}{\partial x_2}-x_3\frac{\partial {\varphi }_2}{\partial x_2};{\epsilon }_{33}=0\\ {\epsilon }_{23}\left(x_1,x_2,0,t\right)={\epsilon }_{32}\left(x_1,x_2,0,t\right)=\frac{1}{2}\left(-{\varphi }_2+\frac{\partial u_3}{\partial x_2}\right)\\ {\epsilon }_{31}\left(x_1,x_2,0,t\right)={\epsilon }_{13}\left(x_1,x_2,0,t\right)=\frac{1}{2}\left(-{\varphi }_1+\frac{\partial u_3}{\partial x_1}\right)\\ {\epsilon }_{12}\left(x_1,x_2,x_3,t\right)={\epsilon }_{21}\left(x_1,x_2,x_3,t\right)=\frac{1}{2}\left(\frac{\partial {\tilde{u}}_1}{\partial x_2}+\frac{\partial {\tilde{u}}_2}{\partial x_1}\right)-\frac{x_3}{2}\left(\frac{\partial {\varphi }_2}{\partial x_2}+\frac{\partial {\varphi }_1}{\partial x_1}\right)\end{array}$ (12)

Based on (12) only ${\sigma }_{33}$, ${\epsilon }_{33}$ do not contribute to strain energy. If we assume linear constitutive theory with plane stress behavior in the plane(x₁x₂ plane), then

$\begin{array}{l}{\sigma }_{11}=D_{11}{\epsilon }_{11}+D_{12}{\epsilon }_{22};{\sigma }_{22}=D_{21}{\epsilon }_{11}+D_{22}{\epsilon }_{22}\\ {\sigma }_{23}=D_{44}\left(2{\epsilon }_{23}\right);{\sigma }_{31}=D_{55}\left(2{\epsilon }_{31}\right);{\sigma }_{12}=D_{44}\left(2{\epsilon }_{12}\right)\\ D_{11}=D_{22}=\frac{E}{1-{\nu }^2};D_{12}=D_{21}=\nu D_{11}\\ D_{44}=D_{55}=D_{66}=G=\frac{E}{2\left(1+\nu \right)}.\end{array}$ (13)

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.

$\begin{array}{l}I={\displaystyle \underset{0}{\overset{t}{\int }}}{\displaystyle \underset{V}{\int }}(\frac{1}{2}{\rho }_0\left({\left(\frac{\partial u_1}{\partial t}\right)}^2+{\left(\frac{\partial u_2}{\partial t}\right)}^2+{\left(\frac{\partial u_3}{\partial t}\right)}^2\right)-{\sigma }_{11}{\epsilon }_{11}\\ \begin{array}{c}\\ \end{array}-{\sigma }_{22}{\epsilon }_{22}-{\sigma }_{23}\left(2{\epsilon }_{23}\right)-{\sigma }_{31}\mathrm{<mo\; stretchy="false">\; (\; </mo>2}{\epsilon }_{31}\mathrm{<mo\; stretchy="false">\; )\; </mo>}-{\sigma }_{12}\left(2{\epsilon }_{12}\right)+{\rho }_0{u}_i{F}_i^b)\text{d}V\text{d}t\end{array}$ (14)

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 ( $\delta I=0$ ). Using integration by parts all differentiation from $\delta {\tilde{u}}_1$, $\delta {\tilde{u}}_2$, $\delta u_3$, $\delta {\varphi }_1$ and $\delta {\varphi }_2$ is transferred to their coefficients. After grouping the coefficients of $\delta {\tilde{u}}_1$, $\delta {\tilde{u}}_2$, $\delta u_3$, $\delta {\varphi }_1$ and $\delta {\varphi }_2$ and

performing integration with respect to x₃ using limits $\left[\frac{-h}{2},\frac{h}{2}\right]$ the coefficients

of $\delta {\tilde{u}}_1$, $\delta {\tilde{u}}_2$, $\delta u_3$, $\delta {\varphi }_1$ and $\delta {\varphi }_2$ 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):

$\begin{array}{l}h{\rho }_0\frac{{\partial }^2{\tilde{u}}_1}{\partial t^2}-\frac{\partial N_{11}}{\partial x_1}-\frac{\partial Q_{12}}{\partial x_2}-h{\rho }_0{F}_1^b=0\\ h{\rho }_0\frac{{\partial }^2{\tilde{u}}_2}{\partial t^2}-\frac{\partial Q_{12}}{\partial x_1}-\frac{\partial N_{22}}{\partial x_2}-h{\rho }_0{F}_2^b=0\\ h{\rho }_0\frac{{\partial }^2{u}_3}{\partial t^2}-\frac{\partial Q_{13}}{\partial x_1}-\frac{\partial Q_{23}}{\partial x_2}-h{\rho }_0{F}_3^b=0\\ {\rho }_0\frac{h^3}{12}\frac{{\partial }^2{\varphi }_1}{\partial t^2}-\frac{\partial M_{11}}{\partial x_1}-\frac{\partial M_{12}}{\partial x_2}+Q_{13}=0\\ {\rho }_0\frac{h^3}{12}\frac{{\partial }^2{\varphi }_2}{\partial t^2}-\frac{\partial M_{12}}{\partial x_2}-\frac{\partial M_{22}}{\partial x_2}+Q_{23}=0\end{array}$ (15)

in which

$\begin{array}{l}{N}_{11}={\displaystyle \underset{-h/2}{\overset{h/2}{\int }}}{\sigma }_{11}\text{d}{x}_3=h\left(D_{11}\frac{\partial {\tilde{u}}_1}{\partial x_1}+D_{12}\frac{\partial {\tilde{u}}_2}{\partial x_2}\right)\\ N_{22}={\displaystyle \underset{-h/2}{\overset{h/2}{\int }}}{\sigma }_{22}\text{d}{x}_3=h\left(D_{21}\frac{\partial {\tilde{u}}_1}{\partial x_1}+D_{22}\frac{\partial {\tilde{u}}_2}{\partial x_2}\right)\end{array}$ (16)

$\begin{array}{l}{Q}_{23}=Q_{32}={\displaystyle \underset{-h/2}{\overset{h/2}{\int }}}{\sigma }_{23}\text{d}{x}_3={\displaystyle \underset{-h/2}{\overset{h/2}{\int }}}{D}_{44}\left(2{\epsilon }_{23}\right)\text{d}{x}_3=h{D}_{44}\left(-{\varphi }_2+\frac{\partial u_3}{\partial x_2}\right)\\ Q_{13}=Q_{31}={\displaystyle \underset{-h/2}{\overset{h/2}{\int }}}{\sigma }_{13}\text{d}{x}_3={\displaystyle \underset{-h/2}{\overset{h/2}{\int }}}{D}_{55}\left(2{\epsilon }_{13}\right)\text{d}{x}_3=h{D}_{55}\left(-{\varphi }_1+\frac{\partial u_3}{\partial x_1}\right)\\ Q_{12}=Q_{21}={\displaystyle \underset{-h/2}{\overset{h/2}{\int }}}{\sigma }_{12}\text{d}{x}_3={\displaystyle \underset{-h/2}{\overset{h/2}{\int }}}{D}_{66}\left(2{\epsilon }_{12}\right)\text{d}{x}_3=h{D}_{66}\left(\frac{\partial {\tilde{u}}_1}{\partial x_2}+\frac{\partial {\tilde{u}}_2}{\partial x_1}\right)\end{array}\}$ (17)

and

$\begin{array}{l}{M}_{11}={\displaystyle \underset{-h/2}{\overset{h/2}{\int }}}{x}_3{\sigma }_{11}\text{d}{x}_3=\frac{-h^3}{12}\left(D_{11}\frac{\partial {\varphi }_1}{\partial x_1}+D_{12}\frac{\partial {\varphi }_2}{\partial x_2}\right)\\ M_{22}={\displaystyle \underset{-h/2}{\overset{h/2}{\int }}}{x}_3{\sigma }_{22}\text{d}{x}_3=\frac{-h^3}{12}\left(D_{21}\frac{\partial {\varphi }_1}{\partial x_1}+D_{22}\frac{\partial {\varphi }_2}{\partial x_2}\right)\\ M_{12}={\displaystyle \underset{-h/2}{\overset{h/2}{\int }}}{x}_3{\sigma }_{12}\text{d}{x}_3={\displaystyle \underset{-h/2}{\overset{h/2}{\int }}}{x}_3{D}_{66}\left(2{\epsilon }_{12}\right)\text{d}{x}_3=\frac{-h^3}{12}{D}_{66}\left(\frac{\partial {\varphi }_1}{\partial x_2}+\frac{\partial {\varphi }_2}{\partial x_1}\right)\end{array}\}$ (18)

Equations (15) constitute the complete mathematical model in ${\tilde{u}}_1,{\tilde{u}}_2,{\tilde{u}}_3,{\varphi }_1$ and ${\varphi }_2$. $N_{11}\mathrm{,}{N}_{22}\mathrm{,}{Q}_{12}\mathrm{,}{M}_{11}\mathrm{,}{M}_{22}\mathrm{,}{M}_{12}$ as functions of ${\tilde{u}}_1,{\tilde{u}}_2,{\tilde{u}}_3,{\varphi }_1$ and ${\varphi }_2$ are defined by (16) - (18).

Remarks.

1) We note that the transverse shear forces Q₁₃ and Q₂₃ are not zero in this mathematical model.

2) This mathematical model when compared with CPT is expected to yield higher transverse deflection u₃ due to shear deformation (additional strain energy due to transverse shear strains ${\epsilon }_{13}$, ${\epsilon }_{23}$ and shear stresses ${\sigma }_{13}$ and ${\sigma }_{23}$ ).

3) We note from (17) that transverse shear forces Q₁₃ and Q₂₃ are not functions of x₃ i.e., these are constant along the thickness of the plate where as their actual distributions are parabolic functions of x₃. Thus, the strain energy due ${\sigma }_{23}$, ${\epsilon }_{23}$ and ${\sigma }_{13}$, ${\epsilon }_{13}$ in this model is not correct(as well known). Shear correction factors [14] [35] 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 (x₁x₂ 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 u₃, Q₁₃, Q₂₃, M₁₁, M₁₂, M₂₂, ${\varphi }_1$ and ${\varphi }_2$. Q₁₃, Q₂₃, M₁₁, M₁₂ and M₂₂ can be eliminated by substituting them in equations three, four and five of (16). This gives rise to three partial differential equations in u₃, ${\varphi }_1$ and ${\varphi }_2$, containing up to second order derivatives of u₃, ${\varphi }_1$ and ${\varphi }_2$.

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 $x_1,x_2,x_3$ 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. [31] [32] [36] 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:

$\begin{array}{l}{\sigma }_{ij}\ne {\sigma }_{ji}\text{for}j\ne i\mathrm{<mo>\; ;\; </mo>}{\sigma }_{ij}={}_s\sigma {}_{ij}+{}_a\sigma {}_{ij}\\ {}_s\sigma {}_{ij}={}_s\sigma {}_{ji}\text{and}{}_a\sigma {}_{ij}=-{}_a\sigma {}_{ji}\text{for}i\ne j\text{and}{}_a\sigma {}_{ij}=0\text{for}j=i\end{array}$ (19)

$m_{mk\mathrm{,}m}=-{ϵ}_{ijk}{\sigma }_{ij}=-{ϵ}_{ijk}\left({}_a{\sigma }_{ij}\right)\left(\text{BAM}\right)$ (20)

$m_{ij}=m_{ji}\left(\text{BMM}\right)\mathrm{.}$ (21)

Linear constitutive theories for ${}_s\sigma$, $m$ and $q$ are given by

${}_s\sigma {}_{ij}=2\mu {\epsilon }_{ij}+\lambda {\epsilon }_{kk}{\delta }_{ij}$ (22)

$m_{kj}=2\underset{˜}{\mu }\left({}_s^{{}_i\Theta }{J}_{\text{}kj}\right)$ (23)

and

$q_i=-k{g}_i\mathrm{.}$ (24)

Equations (A2), (20), (22) and (23) constitute eighteen partial differential equations in eighteen dependent variables $u_i$, ${}_s\sigma {}_{ij}$, ${}_a\sigma {}_{ij}$ and $m_{ij}$, hence the mathematical model has closure. $\mu$, $\lambda$ and $\underset{˜}{\mu }$ are material coefficients.

3.2.2. Non-Classical Continuum Mechanics with Internal and Cosserat Rotations

In this continuum theory, we also consider internal rotations $\left({}_i{\Theta }_{x_1}{,}_i{\Theta }_{x_2}{,}_i{\Theta }_{x_3}\right)$ 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 ( ${}_e{\Theta }_{x_i}$, 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 ${}_t\Theta$, sum of ${}_t\Theta$ and ${}_e\Theta$ 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 $\sigma$ is not symmetric.

Following Section A3 in “Appendix A”, we have the following from the remaining balance laws and the linear constitutive theories:

$\begin{array}{l}{\sigma }_{ij}\ne {\sigma }_{ji}\text{for}j\ne i\mathrm{<mo>\; ;\; </mo>}{\sigma }_{ij}={}_s\sigma {}_{ij}+{}_a\sigma {}_{ij}\\ {}_s\sigma {}_{ij}={}_s\sigma {}_{ji}\text{and}{}_a{\sigma }_{ij}={-}_a{\sigma }_{ji}\text{for}i\ne j\text{and}{}_a{\sigma }_{ij}=0\text{for}j=i\end{array}$ (25)

$m_{mk\mathrm{,}m}=-{ϵ}_{ijk}{\sigma }_{ij}=-{ϵ}_{ijk}\left({}_a{\sigma }_{ij}\right)\left(\text{BAM}\right)$ (26)

$m_{ij}=m_{ji}\left(\text{BMM}\right)\mathrm{.}$ (27)

Linear constitutive theories for ${}_s\sigma$, $m$, ${}_a\sigma$ and $q$ are given by [31] [32] [33] [34] [36]:

${}_s\sigma =2\mu \epsilon +\lambda \left(\text{tr}\epsilon \right)I$ (28)

$m=2\underset{˜}{\mu }\left({}_s^{{}_t\Theta }\text{}\text{}J\right)+\underset{˜}{\lambda }\text{tr}\left({}_s^{{}_t\Theta }\text{}\text{}J\right)$ (29)

${}_a\sigma =2\underset{˜}{\beta }\left({}_a^{t}r\right)$ (30)

$q=-kg\mathrm{.}$ (31)

Equations (A2), (26), (28) - (30) are a system of twenty one partial differential equations in twenty one dependent variables: $u_i$, ${}_s\sigma {}_{ij}$, ${}_a\sigma {}_{ij}$, $m_{ij}$ and ${}_t{\Theta }_{x_i}$.

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 ${ℝ}^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 ${\epsilon }_{11}$, ${\epsilon }_{22}$ and ${\epsilon }_{12}$ are obtained as shown in (2) and we also find that the ${\epsilon }_{33}=0$, ${\epsilon }_{23}=0$ and ${\epsilon }_{31}=0$. 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 ${ℝ}^3$, then ${\sigma }_{33}\ne 0$ even though ${\epsilon }_{33}=0$. On the other hand if we assume plane stress, then ${\sigma }_{33}=0$ and the only nonzero components of stress tensor are ${\sigma }_{11}$, ${\sigma }_{22}$ and ${\sigma }_{12}$. In the present work, we consider plate deformation in ${ℝ}^3$, hence we have ${\sigma }_{33}\ne 0$. In addition to (3), we have ${\sigma }_{33}=D_{31}{\epsilon }_{11}+D_{32}{\epsilon }_{22}$ in which $D_{31}=D_{13}=\nu D_{11}$.

We substitute from (2) and (3) in the balance of linear momenta (A2), and consider integral over the volume $V=A\times \left[-h/2\mathrm{,}h/2\right]$ of the plate. In which A is the area of the middle plane ( $x_3=0$ ) of the plate. We integrate with respect to $x_3\in \left[-h/2\mathrm{,}h/2\right]$. The resulting integral is valid over arbitrary area A, hence the integrand must be zero, which gives us the following equations:

$\begin{array}{l}{\rho }_{0}h\frac{{\partial }^2{\tilde{u}}_1}{\partial t^2}-{\rho }_{0}h{F}_1^b-\frac{\partial N_{11}}{\partial x_1}-\frac{\partial Q_{12}}{\partial x_2}=0\\ {\rho }_{0}h\frac{{\partial }^2{\tilde{u}}_2}{\partial t^2}-{\rho }_{0}h{F}_2^b-\frac{\partial Q_{12}}{\partial x_1}-\frac{\partial N_{22}}{\partial x_2}=0\\ {\rho }_{0}h\frac{{\partial }^2{\tilde{u}}_3}{\partial t^2}-{\rho }_{0}h{F}_3^b-\frac{\partial N_{33}}{\partial x_3}=0\end{array}$ (32)

where,

$\begin{array}{l}{N}_{11}={\displaystyle \underset{-h/2}{\overset{h/2}{\int }}}{\sigma }_{11}\text{d}{x}_3;N_{22}={\displaystyle \underset{-h/2}{\overset{h/2}{\int }}}{\sigma }_{22}\text{d}{x}_3;Q_{12}={\displaystyle \underset{-h/2}{\overset{h/2}{\int }}}{\sigma }_{12}\text{d}{x}_3\\ \frac{\partial N_{33}}{\partial x_3}={\displaystyle \underset{-h/2}{\overset{h/2}{\int }}}\left(-\frac{\partial {\sigma }_{33}}{\partial x_3}\right)\text{d}{x}_3=h\left(D_{31}\frac{{\partial }^2{u}_2}{\partial x_1^2}+D_{32}\frac{{\partial }^2{u}_3}{\partial x_2^2}\right).\end{array}$ (33)

Equations (32) are three partial differential equations in ${\tilde{u}}_1,{\tilde{u}}_2$ and u₃.

Remarks

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 $N_{33}=0$. 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 ${\sigma }_{23}$ and ${\sigma }_{31}$ and hence the corresponding transverse shear forces Q₂₃ and Q₃₁ 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 ${ℝ}^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 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 ${}^{d}J$ into symmetric and skew symmetric components.

${}^{d}J={}_s{}^{d}J+{}_a{}^{d}J$ (34)

in which components of ${}_a{}^{d}J$ are in fact related to the internal rotations ${}_i{\Theta }_{x_j}$.

First, we define the internal rotations ${}_i{\Theta }_{x_j}$

$\begin{array}{l}\nabla \times u=e_1\left({}_i{\Theta }_{x_1}\right)+e_2\left({}_i{\Theta }_{x_2}\right)+e_3\left({}_i{\Theta }_{x_3}\right)\\ \text{or}\nabla \times u=e_1\left(\frac{\partial u_3}{\partial x_2}-\frac{\partial u_2}{\partial x_3}\right)+e_2\left(\frac{\partial u_1}{\partial x_3}-\frac{\partial u_3}{\partial x_1}\right)+e_3\left(\frac{\partial u_2}{\partial x_1}-\frac{\partial u_1}{\partial x_2}\right)\mathrm{.}\end{array}$ (35)

Using the kinematic assumptions (1), we obtain following for the internal rotations defined in (35)

$\begin{array}{l}{}_i{\Theta }_{x_1}=2\frac{\partial u_3}{\partial x_2}\\ {}_i{\Theta }_{x_2}=-2\frac{\partial u_3}{\partial x_1}\\ {}_i{\Theta }_{x_3}=\left(\frac{\partial {\tilde{u}}_2}{\partial x_1}-\frac{\partial {\tilde{u}}_1}{\partial x_2}\right)\mathrm{.}\end{array}$ (36)

The components of ${}_a{}^{d}J$ can now be defined using (36)

$\begin{array}{l}{}_a{}^{d}J{}_{12}=-{}_a{}^{d}J{}_{21}=\frac{-{}_i\Theta {}_{x_3}}{2}\\ {}_a{}^{d}J{}_{13}=-{}_a{}^{d}J{}_{31}=\frac{{}_i\Theta {}_{x_2}}{2}\\ {}_a{}^{d}J{}_{23}=-{}_a{}^{d}J{}_{32}=\frac{-{}_i\Theta {}_{x_1}}{2}\end{array}$ (37)

where ${}_i\Theta {}_{x_1}\mathrm{,}{}_i\Theta {}_{x_2}\mathrm{,}{}_i\Theta {}_{x_3}$ are internal rotations [31] [32] [33] [34] [36] 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 $\sigma$ is decomposed into symmetric ${}_s\sigma$ and skew symmetric ${}_a\sigma$ tensors.

$\sigma ={}_s\sigma +{}_a\sigma$ (38)

The components ${}_s\sigma$ are same as those defined by (2). Substituting the total stress tensor $\sigma$ from (38) in the balance of linear momenta and integrating with respect to $x_3$ $\left(\left[-h/2\mathrm{,}h/2\right]\right)$, we obtain the following:

$\begin{array}{l}h{\rho }_0\frac{{\partial }^2{\tilde{u}}_1}{\partial t^2}-\frac{\partial N_{11}}{\partial x_1}-\frac{\partial Q_{21}}{\partial x_2}-\frac{\partial Q_{31}}{\partial x_3}-{\rho }_{0}h{F}_1^b=0\\ h{\rho }_0\frac{{\partial }^2{\tilde{u}}_2}{\partial t^2}-\frac{\partial Q_{12}}{\partial x_1}-\frac{\partial N_{22}}{\partial x_2}-\frac{\partial Q_{32}}{\partial x_3}-{\rho }_{0}h{F}_2^b=0\\ h{\rho }_0\frac{{\partial }^2{u}_3}{\partial t^2}-\frac{\partial Q_{13}}{\partial x_1}-\frac{\partial Q_{23}}{\partial x_2}-\frac{\partial N_{33}}{\partial x_3}-{\rho }_{0}h{F}_3^b=0\end{array}$ (39)

$\begin{array}{l}{N}_{11}={}_{s}N{}_{11}={\displaystyle \underset{-h/2}{\overset{{}_s}{\int }}}{}_s\sigma {}_{11}\text{d}{x}_3;N_{22}={}_{s}N{}_{22}={\displaystyle \underset{-h/2}{\overset{{}_s}{\int }}}{}_s\sigma {}_{22}\text{d}{x}_3\\ \frac{\partial N_{33}}{\partial x_3}=\frac{\partial \left({}_s{N}_{33}\right)}{\partial x_3}=h\left(D_{31}\frac{{\partial }^2{u}_2}{\partial x_1^2}+D_{32}\frac{{\partial }^2{u}_3}{\partial x_2^2}\right)\end{array}$ (40)

$Q_{12}\mathrm{,}{Q}_{21}\mathrm{,}{Q}_{31}\mathrm{,}{Q}_{13}\mathrm{,}{Q}_{23}$ and $Q_{32}$ are defined later.

The internal rotation gradient tensor ${}^{{}_i\Theta }J$

$\left[{}^{{}_i\Theta }J\right]=\left[\frac{\partial \left\{{}_i\Theta \right\}}{\partial \left\{x\right\}}\right]$ (41)

is decomposed into symmetric and skew symmetric tensors.

${}^{{}_i\Theta }J={}_s{}^{{}_i\Theta }J+{}_a{}^{{}_i\Theta }J$ (42)

in which

$\begin{array}{l}{}_s{}^{{}_i\Theta }J{}_{11}=\frac{\partial \left({}_i{\Theta }_{x_1}\right)}{\partial x_1}=2\frac{{\partial }^2{u}_3}{\partial x_1\partial x_2}\\ {}_s{}^{{}_i\Theta }J{}_{22}=\frac{\partial \left({}_i{\Theta }_{x_2}\right)}{\partial x_2}=-2\frac{{\partial }^2{u}_3}{\partial x_1\partial x_2}\\ {}_s{}^{{}_i\Theta }J{}_{33}=\frac{\partial \left({}_i{\Theta }_{x_3}\right)}{\partial x_3}=0.\\ {}_s{}^{{}_i\Theta }J{}_{12}=\frac{1}{2}\left(\frac{\partial \left({}_i{\Theta }_{x_1}\right)}{\partial x_2}+\frac{\partial \left({}_i{\Theta }_{x_2}\right)}{\partial x_1}\right)=\left(\frac{{\partial }^2{u}_3}{\partial x_2^2}-\frac{{\partial }^2{u}_3}{\partial x_1^2}\right)\\ {}_s{}^{{}_i\Theta }J{}_{13}=\frac{1}{2}\left(\frac{\partial \left({}_i{\Theta }_{x_1}\right)}{\partial x_3}+\frac{\partial \left({}_i{\Theta }_{x_3}\right)}{\partial x_1}\right)=\frac{1}{2}\left(\frac{{\partial }^2{\tilde{u}}_2}{\partial x_1^2}-\frac{{\partial }^2{\tilde{u}}_1}{\partial x_1{x}_2}\right)\\ {}_s^{{}_i\Theta }\text{}\text{}{J}_{23}=\frac{1}{2}\left(\frac{\partial \left({}_i{\Theta }_{x_2}\right)}{\partial x_3}+\frac{\partial \left({}_i{\Theta }_{x_3}\right)}{\partial x_2}\right)=\frac{1}{2}\left(\frac{{\partial }^2{\tilde{u}}_2}{\partial x_1{x}_2}-\frac{{\partial }^2{\tilde{u}}_1}{\partial x_2^2}\right)\end{array}$ (43)

From balance of angular momenta we have

$m_{mk\mathrm{,}m}+{ϵ}_{ijk}{\sigma }_{ij}=0.$ (44)

Balance of moment of moments show that [37] [38]

$m_{ij}=m_{ji}\mathrm{.}$ (45)

The constitutive theories for ${}_s\sigma$ and $m$ are [31] [32] [33] [34] [36]

${}_s\sigma {}_{ij}=2\mu {\epsilon }_{ij}+\lambda {\epsilon }_{kk}{\delta }_{ij}$ (46)

$m_{ij}=2\mu \left({}_s{}^{{}_i\Theta }J{}_{ij}\right)\mathrm{<mo>\; ;\; </mo>}\text{as}\left({}_s{}^{{}_i\Theta }J{}_{kk}=0\right)\mathrm{.}$ (47)

The constitutive theory for ${}_s\sigma {}_{11}={\sigma }_{11}$, ${}_s\sigma {}_{22}={\sigma }_{22}$ and ${}_s\sigma {}_{12}={}_s\sigma {}_{21}$ are obtained using (46) and are same as those in (3). ${}_s\sigma {}_{33}={\sigma }_{33}\ne 0$ if we consider the deformation of plate bending in ${ℝ}^3$ (discussed earlier). Using (47) and (42), explicit expressions for the constitutive theory for $m_{ij}$ can be obtained in terms of the symmetric part of the rotation gradient tensor.

Using these expressions for $m_{ij}$ in BAM (44) and dividing them by 2 and then integrating with respect to $x_3$ over $\left(\left[-h/2\mathrm{,}h/2\right]\right)$ we obtain

$\begin{array}{l}h\left({\underset{˜}{m}}_{\mathrm{11,1}}+{\underset{˜}{m}}_{\mathrm{21,2}}\right)=-{}_{a}Q{}_{23}\left(x_1\mathrm{,}{x}_2\right)\\ h\left({\underset{˜}{m}}_{\mathrm{12,1}}+{\underset{˜}{m}}_{\mathrm{22,2}}\right)=-{}_{a}Q{}_{13}\left(x_1\mathrm{,}{x}_2\right)\\ h\left({\underset{˜}{m}}_{\mathrm{13,1}}+{\underset{˜}{m}}_{\mathrm{23,2}}\right)=-{}_{a}Q{}_{12}\left(x_1\mathrm{,}{x}_2\right)\\ \text{in}\text{which}{m}_{mk\mathrm{,}m}=\frac{1}{2}{\displaystyle \underset{-h/2}{\overset{h/2}{\int }}}{m}_{mk\mathrm{,}m}\text{d}{x}_3\\ \text{and}{}_{a}Q{}_{13}=-{}_{a}Q{}_{31}\mathrm{<mo>\; ;\; </mo>}{}_{a}Q{}_{21}=-{}_{a}Q{}_{12}\mathrm{<mo>\; ;\; </mo>}{}_{a}Q{}_{32}=-{}_{a}Q{}_{23}\end{array}$ (48)

$\begin{array}{l}{Q}_{12}={}_{s}Q{}_{12}+{}_{a}Q{}_{12}\\ {}_{s}Q{}_{12}={\displaystyle \underset{\frac{-h}{2}}{\overset{\frac{h}{2}}{\int }}}{}_s\sigma {}_{12}\text{d}{x}_3\\ {}_{a}Q{}_{12}=-{}_{a}Q{}_{21}={\displaystyle \underset{\frac{-h}{2}}{\overset{\frac{h}{2}}{\int }}}{}_a\sigma {}_{12}\text{d}{x}_3\end{array}$ (49)

$\begin{array}{l}\frac{\partial Q_{31}}{\partial x_3}=\frac{{\partial }_s{Q}_{31}}{\partial x_3}+\frac{{\partial }_a{Q}_{31}}{\partial x_3}=0\\ \frac{\partial Q_{23}}{\partial x_3}=\frac{{\partial }_s{Q}_{23}}{\partial x_3}+\frac{{\partial }_a{Q}_{23}}{\partial x_3}=0\\ \frac{\partial Q_{23}}{\partial x_2}=\frac{{\partial }_s{Q}_{23}}{\partial x_2}+\frac{{\partial }_a{Q}_{23}}{\partial x_2}=\frac{{\partial }_a{Q}_{23}}{\partial x_2}\\ \frac{\partial Q_{31}}{\partial x_1}=\frac{{\partial }_s{Q}_{31}}{\partial x_1}+\frac{{\partial }_a{Q}_{31}}{\partial x_1}=\frac{{\partial }_a{Q}_{31}}{\partial x_1}.\end{array}$ (50)

We note that $\frac{\partial Q_{31}}{\partial x_3}$ and $\frac{\partial Q_{32}}{\partial x_3}$ terms in (39) are zeros as $Q_{31}=Q_{31}\left(x_1\mathrm{,}{x}_2\mathrm{,0,}t\right)$ and $Q_{32}=Q_{32}\left(x_1\mathrm{,}{x}_2\mathrm{,0,}t\right)$.

This mathematical model consists of BLM (39) (3 equations), BAM (44) (3 equations), constitutive theories for ${}_s\sigma$ (or ${}_{s}N{}_{11}\mathrm{,}{}_{s}N{}_{22}\mathrm{,}{}_{s}Q{}_{12}\mathrm{,}{}_{s}N{}_{33}$, Equations (40) and (49)) (4 equations) and constitutive theories for $m$ ( $m_{11}\mathrm{,}{m}_{12}\mathrm{,}{m}_{22}\mathrm{,}{m}_{13}\mathrm{,}{m}_{23}$, Equations (47)) (5 equations). These are a total of fifteen partial differential equations in fifteen variables ${\tilde{u}}_1\mathrm{,}{\tilde{u}}_2\mathrm{,}{u}_3\mathrm{<mo>\; ;\; </mo>}{}_{s}N{}_{11}\mathrm{,}{}_{s}N{}_{22}\mathrm{,}{}_{s}N{}_{33}\mathrm{,}{}_{s}Q{}_{12}\mathrm{,}{}_{a}Q{}_{12}\mathrm{,}{}_{a}Q{}_{31}\mathrm{,}{}_{a}Q{}_{23}$; $m_{11},m_{12},m_{22},m_{23},m_{13}$, hence the mathematical model has closure.

Remarks

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 $\sigma$.

3) ${}_s{Q}_{31}$ and ${}_s{Q}_{23}$ are zero but ${}_a{Q}_{31}$ and ${}_a{Q}_{23}$ are nonzero due to nonzero ${}_a{\sigma }_{31}$ and ${}_a{\sigma }_{23}$. Existence of ${}_a{Q}_{31}$ or ${}_a{\sigma }_{31}$ and ${}_a{Q}_{23}$ or ${}_a{\sigma }_{23}$ is due to Cauchy moment tensor (necessitated due to the presence of internal rotations).

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 $M_{11}$, $M_{22}$ and $M_{12}$ obtained in the currently used mathematical models by integrating moments of ${\sigma }_{11}$, ${\sigma }_{22}$ and ${\sigma }_{12}$ in the $x_3$ direction for $x_3\in \left(-h/2\mathrm{,}h/2\right)$.

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 $x_3$ direction between −h/2 to h/2 we obtain

$Q_{12}\left(x_1,x_2\right)=Q_{21};Q_{31}\left(x_1,x_2\right)=Q_{13};Q_{32}\left(x_1,x_2\right)=Q_{23}$ (51)

$\begin{array}{l}h{\rho }_0\frac{{\partial }^2{\tilde{u}}_1}{\partial t^2}-\frac{\partial N_{11}}{\partial x_1}-\frac{\partial Q_{21}}{\partial x_2}-{\rho }_{0}h{F}_1^b=0\\ h{\rho }_0\frac{{\partial }^2{\tilde{u}}_2}{\partial t^2}-\frac{\partial Q_{12}}{\partial x_1}-\frac{\partial N_{22}}{\partial x_2}-{\rho }_{0}h{F}_2^b=0\\ h{\rho }_0\frac{{\partial }^2{u}_3}{\partial t^2}-\frac{\partial Q_{13}}{\partial x_1}-\frac{\partial Q_{23}}{\partial x_2}-\frac{\partial N_{33}}{\partial x_3}-{\rho }_{0}h{F}_3^b=0\end{array}$ (52)

We note that $\frac{\partial Q_{31}}{\partial x_3}=0$ and $\frac{\partial Q_{32}}{\partial x_3}=0$ as $Q_{31}=Q_{31}\left(x_1\mathrm{,}{x}_2\mathrm{,0,}t\right)$ and $Q_{32}=Q_{32}\left(x_1\mathrm{,}{x}_2\mathrm{,0,}t\right)$ are functions of $x_1,x_2$ and t only. $\frac{\partial N_{33}}{\partial x_3}$ is given in Equation (40).

Forces $N_{11}$, $N_{22}$, $Q_{12}=Q_{21}$, $Q_{31}=Q_{13}$ and $Q_{23}=Q_{32}$ are defined by Equation (16).

Remarks

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 ${\tilde{u}}_1$, ${\tilde{u}}_2$, $u_3$; $N_{11}$, $N_{22}$, $Q_{13}=Q_{31}$, $Q_{12}=Q_{21}$, $Q_{23}=Q_{32}$; ${\varphi }_1$ and ${\varphi }_2$. 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 ${\varphi }_1$ and ${\varphi }_2$ 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) ${}^{d}J={}_s{}^{d}J+{}_a{}^{d}J$ or alternatively expand $\nabla \times u$.

Using $\nabla \times u$ and the kinematic relations (11) we can obtain

$\begin{array}{l}\nabla \times u=e_1\left({}_t{\Theta }_{x_1}\right)+e_2\left({}_t{\Theta }_{x_2}\right)+e_3\left({}_t{\Theta }_{x_3}\right)\\ \text{in}\text{which}{}_t{\Theta }_{x_1}=\left(\frac{\partial u_3}{\partial x_2}+{\varphi }_2\right)\\ {}_t{\Theta }_{x_2}=\left(-{\varphi }_1-\frac{\partial u_3}{\partial x_1}\right)\\ {}_t{\Theta }_{x_3}=\left(\frac{\partial {\tilde{u}}_2}{\partial x_1}-\frac{\partial {\tilde{u}}_1}{\partial x_2}\right)+x_3\left(-\frac{\partial {\varphi }_2}{\partial x_1}+\frac{\partial {\varphi }_1}{\partial x_2}\right)\end{array}$ (53)

The total rotations ${}_t{\Theta }_{x_1}$, ${}_t{\Theta }_{x_2}$ and ${}_t{\Theta }_{x_3}$ 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 ${\varphi }_1$ and ${\varphi }_2$, 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 ${ℝ}^3$, then coefficients $D_{ij}$ are modified and are due to linear constitutive theory in ${ℝ}^3$ giving rise to ${\sigma }_{33}\ne 0$ and

${\sigma }_{33}=D_{31}\frac{\partial u_1}{\partial x_1}+D_{32}\frac{\partial u_2}{\partial x_2}={\sigma }_{33}\left(x_1,x_2,x_3,t\right)$ (54)

The stress tensor $\sigma$ is decomposed into symmetric and skew symmetric tensors ${}_s\sigma$ and ${}_a\sigma$ :

$\sigma ={}_s\sigma +{}_a\sigma$ (55)

Definitions of the components of the ${}_s\sigma$ is same as in (13). Substituting total stress tensor $\sigma$ into the balance of linear momenta and integrating with respect to $x_3\left(\left[-h/2\mathrm{,}h/2\right]\right)$ we obtain (39) in which

$\begin{array}{l}{N}_{11}={}_{s}N{}_{11}={\displaystyle \underset{-h/2}{\overset{h/2}{\int }}}{}_s{\sigma }_{11}\text{d}{x}_3;N_{22}={}_{s}N{}_{22}={\displaystyle \underset{-h/2}{\overset{h/2}{\int }}}{}_s{\sigma }_{22}\text{d}{x}_3\\ \frac{\partial N_{33}}{\partial x_3}=\frac{\partial \left({}_s{N}_{33}\right)}{\partial x_3}=h\left(D_{31}\frac{{\partial }^2{u}_3}{\partial x_1^2}+D_{32}\frac{{\partial }^2{u}_3}{\partial x_2^2}\right).\end{array}$ (56)

$Q_{12}\mathrm{,}{Q}_{21}\mathrm{,}{Q}_{31}\mathrm{,}{Q}_{13}\mathrm{,}{Q}_{23}$ and $Q_{32}$ appearing in (39) are defined by (67) - (71).

The gradient ${}^{{}_t\Theta }J$ of total rotation ${}_t\Theta$ is defined by

(57)

Components of $\left[{}^{{}_t\Theta }J\right]$ can be obtained using definition of ${}_t{\Theta }_{x_i}$ in (53) and using (57):

$\begin{array}{l}{}^{{}_t\Theta }\text{}\text{}{J}_{11}=\frac{\partial \left({}_t{\Theta }_{x_1}\right)}{\partial x_1}=\left(\frac{{\partial }^2{u}_3}{\partial x_1\partial x_2}+\frac{\partial {\varphi }_2}{\partial x_1}\right)\\ {}^{{}_t\Theta }\text{}{J}_{12}=\frac{\partial \left({}_t{\Theta }_{x_1}\right)}{\partial x_2}=\frac{{\partial }^2{u}_3}{\partial x_2^2}+\frac{\partial {\varphi }_2}{\partial x_2}\\ {}^{{}_t\Theta }\text{}{J}_{13}=\frac{\partial \left({}_t{\Theta }_{x_1}\right)}{\partial x_3}=0\\ {}^{{}_t\Theta }\text{}\text{}{J}_{21}=\frac{\partial \left({}_t{\Theta }_{x_2}\right)}{\partial x_1}=\left(-\frac{\partial {\varphi }_1}{\partial x_1}-\frac{{\partial }^2\text{}{u}_3}{\partial x_1^2}\right)\\ {}^{{}_t\Theta }\text{}\text{}{J}_{22}=\frac{\partial \left({}_t{\Theta }_{x_2}\right)}{\partial x_2}=\left(-\frac{\partial {\varphi }_1}{\partial x_2}-\frac{{\partial }^2{u}_3}{\partial x_1\partial x_2}\right)\end{array}$

$\begin{array}{l}{}^{{}_t\Theta }\text{}\text{}{J}_{23}=\frac{\partial \left({}_t{\Theta }_{x_2}\right)}{\partial x_3}=0\\ {}^{{}_t\Theta }\text{}\text{}{J}_{31}=\frac{\partial \left({}_t{\Theta }_{x_3}\right)}{\partial x_1}=\left(\frac{{\partial }^2{\tilde{u}}_2}{\partial x_1^2}-\frac{{\partial }^2{\tilde{u}}_1}{\partial x_1\partial x_2}\right)+x_3\left(-\frac{{\partial }^2{\varphi }_2}{\partial x_1^2}+\frac{{\partial }^2{\varphi }_1}{\partial x_1\partial x_2}\right)\\ {}^{{}_t\Theta }\text{}\text{}{J}_{32}=\frac{\partial \left({}_t{\Theta }_{x_3}\right)}{\partial x_2}=\left(\frac{{\partial }^2{\tilde{u}}_2}{\partial x_1\partial x_2}-\frac{{\partial }^2{\tilde{u}}_1}{\partial x_2^2}\right)+x_3\left(-\frac{{\partial }^2{\varphi }_2}{\partial x_1\partial x_2}+\frac{{\partial }^2{\varphi }_1}{\partial x_2^2}\right)\\ {}^{{}_t\Theta }\text{}\text{}{J}_{33}=\frac{\partial \left({}_t{\Theta }_{x_3}\right)}{\partial x_3}=\left(-\frac{\partial {\varphi }_2}{\partial x_1}+\frac{\partial {\varphi }_1}{\partial x_2}\right).\end{array}$ (58)

Decomposition of ${}^{{}_t\Theta }J$ into symmetric ( ${}_s{}^{{}_t\Theta }J$ ) and skew symmetric ( ${}_a{}^{{}_t\Theta }J$ ) tensors gives

${}^{{}_t\Theta }J={}_s{}^{{}_t\Theta }J+{}_a{}^{{}_t\Theta }J$ (59)

in which

$\begin{array}{l}{}_s{}^{{}_t\Theta }J=\frac{1}{2}\left({}^{{}_t\Theta }J+{}^{{}_t\Theta }J{}^{\text{T}}\right)\\ {}_a{}^{{}_t\Theta }J=\frac{1}{2}\left({}^{{}_t\Theta }J-{}^{{}_t\Theta }J{}^{\text{T}}\right)\mathrm{.}\end{array}$

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 [31] [32] [36] we have

$m_{mk\mathrm{,}m}+{ϵ}_{ijk}{\sigma }_{ij}=0.$ (60)

From the balance of moment of moments balance law [37] [38], we can show that

$m_{ij}=m_{ji}\mathrm{.}$ (61)

The linear constitutive theories for ${}_s\sigma$, $m$ and ${}_a\sigma$ are given by [31] [32] [33] [34] [36]

${}_s\sigma {}_{ij}=2\mu {\epsilon }_{ij}+\lambda {\epsilon }_{kk}{\delta }_{ij}$ (62)

$m_{ij}=2\underset{˜}{\mu }\left({}_s{}^{{}_t\Theta }J{}_{ij}\right)+\underset{˜}{\lambda }{}_s{}^{{}_t\Theta }J{}_{kk}{\delta }_{ij}$ (63)

${}_a\sigma {}_{ij}=2\underset{˜}{\beta }\left({}_a{}^{{}_t\Theta }J{}_{ij}\right)$ (64)

where $\mu$, $\lambda$ are Lame’s constants and $\underset{˜}{\mu }$, $\underset{˜}{\lambda }$ and $\underset{˜}{\beta }$ 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 $x_3$ $\left(\left[-h/2\mathrm{,}h/2\right]\right)$ to obtain the forces per unit length used in the balance of linear momenta (39). First, from (60):

$\begin{array}{l}{}_{a}Q{}_{12}={\underset{˜}{m}}_{m3,m}=-{}_{a}Q{}_{21}\\ {}_{a}Q{}_{23}={\underset{˜}{m}}_{m1,m}=-{}_{a}Q{}_{32}\\ {}_{a}Q{}_{31}={\underset{˜}{m}}_{m2,m}=-{}_{a}Q{}_{13}\end{array}$ (65)

in which

${\underset{˜}{m}}_{mk,m}=\frac{1}{2}{\displaystyle \underset{-h/2}{\overset{h/2}{\int }}}{m}_{mk,m}\text{d}{x}_3$ (66)

${}_a{Q}_{12}={\displaystyle \underset{-h/2}{\overset{h/2}{\int }}}{}_a{\sigma }_{12}\text{d}{x}_3;{}_a{Q}_{23}={\displaystyle \underset{-h/2}{\overset{h/2}{\int }}}{}_a{\sigma }_{23}\text{d}{x}_3;{}_a{Q}_{31}={\displaystyle \underset{-h/2}{\overset{h/2}{\int }}}{}_a{\sigma }_{31}\text{d}{x}_3.$ (67)

We note that $N_{11}$, $N_{22}$, and $N_{33}$ due to ${\sigma }_{11}={}_s\sigma {}_{11}$, ${\sigma }_{22}={}_s\sigma {}_{22}$, ${\sigma }_{33}={}_s\sigma {}_{33}$ are already defined in (40).

Furthermore,

$Q_{12}={}_{s}Q{}_{12}+{}_{a}Q{}_{12};Q_{23}={}_{s}Q{}_{23}+{}_{a}Q{}_{23};Q_{31}={}_{s}Q{}_{31}+{}_{a}Q{}_{31}$ (68)

in which

$\begin{array}{l}{}_{s}Q{}_{12}={}_{s}Q{}_{21}={\displaystyle \underset{\frac{h}{2}}{\overset{\frac{h}{2}}{\int }}}{}_s{\sigma }_{12}\text{d}{x}_3\\ {}_{s}Q{}_{23}={}_{s}Q{}_{32}={\displaystyle \underset{\frac{h}{2}}{\overset{\frac{h}{2}}{\int }}}{}_s{\sigma }_{23}\text{d}{x}_3\\ {}_{s}Q{}_{31}={}_{s}Q{}_{13}={\displaystyle \underset{\frac{h}{2}}{\overset{\frac{h}{2}}{\int }}}{}_s{\sigma }_{31}\text{d}{x}_3\end{array}$ (69)

We substitute constitutive theory for ${}_s{\sigma }_{12}$, ${}_s{\sigma }_{23}$ and ${}_s{\sigma }_{31}$ in (69).

From (63) after integration of $\frac{1}{2}{m}_{ij}=\frac{2\underset{˜}{\mu }}{2}\left({}_s{}^{{}_t\Theta }J{}_{ij}\right)+\frac{\underset{˜}{\lambda }}{2}\left({}_s{}^{{}_t\Theta }J{}_{kk}\right)$ we obtain

${\underset{˜}{m}}_{ij}=\underset{˜}{\mu }{\displaystyle \underset{\frac{h}{2}}{\overset{\frac{h}{2}}{\int }}}\left({}_s{}^{{}_t\Theta }J{}_{ij}\right)\text{d}{x}_3+\frac{\lambda }{2}{\displaystyle \underset{\frac{h}{2}}{\overset{\frac{h}{2}}{\int }}}\left({}_s{}^{{}_t\Theta }\text{}{J}_{kk}\right)\text{d}{x}_3\mathrm{.}$ (70)

Likewise integrating (64) with respect to $x_3$ between the limits $\left[-h/2\mathrm{,}h/2\right]$ we obtain

$\begin{array}{l}{}_{a}Q{}_{12}=2\beta {\displaystyle \underset{-h/2}{\overset{h/2}{\int }}}\left({}_a{}^{{}_t\Theta }J{}_{12}\right)\text{d}{x}_3\\ {}_{a}Q{}_{23}=2\beta {\displaystyle \underset{-h/2}{\overset{h/2}{\int }}}\left({}_a{}^{{}_t\Theta }J{}_{23}\right)\text{d}{x}_3\\ {}_{a}Q{}_{31}=2\beta {\displaystyle \underset{-h/2}{\overset{h/2}{\int }}}\left({}_a{}^{{}_t\Theta }J{}_{31}\right)\text{d}{x}_3\end{array}$ (71)

We finally rewrite the balance of linear momenta using symmetric and antisymmetric decomposition of shear forces.

$\begin{array}{l}h{\rho }_0\frac{{\partial }^2{\tilde{u}}_1}{\partial t^2}-\frac{\partial N_{11}}{\partial x_1}-\frac{\partial \left({}_s{Q}_{21}+{}_{a}Q{}_{21}\right)}{\partial x_2}-\frac{\partial \left({}_s{Q}_{31}+{}_{a}Q{}_{31}\right)}{\partial x_3}-{\rho }_{0}h{F}_1^b=0\\ h{\rho }_0\frac{{\partial }^2{\tilde{u}}_2}{\partial t^2}-\frac{\partial \left({}_s{Q}_{12}+{}_{a}Q{}_{12}\right)}{\partial x_1}-\frac{\partial N_{22}}{\partial x_2}-\frac{\partial \left({}_s{Q}_{32}+{}_{a}Q{}_{32}\right)}{\partial x_3}-{\rho }_{0}h{F}_2^b=0\\ h{\rho }_0\frac{{\partial }^2{u}_3}{\partial t^2}-\frac{\partial \left({}_s{Q}_{13}+{}_{a}Q{}_{13}\right)}{\partial x_1}-\frac{\partial \left({}_s{Q}_{23}+{}_{a}Q{}_{23}\right)}{\partial x_2}-\frac{\partial N_{33}}{\partial x_3}-{\rho }_{0}h{F}_3^b=0\end{array}$ (72)

The complete mathematical model consists of twenty equations. BLM (3 equations), BAM (3 equations), constitutive theory for ${}_s{\sigma }_{ij}$ defining $N_{11}$, $N_{22}$, $N_{33}$, ${}_s{Q}_{12}$, ${}_s{Q}_{23}$, ${}_s{Q}_{31}$ (through (40) and (69)) (6 equations), constitutive theory for ${}_a{\sigma }_{ij}$ defining ${}_a{Q}_{12}$, ${}_a{Q}_{23}$, ${}_a{Q}_{31}$ (3 equations), constitutive theory for $m_{ij}$ (5 equations) in dependent variables: ${\tilde{u}}_1\mathrm{,}{\tilde{u}}_2$, $u_3$ (3); $m_{ij}$ (5); $N_{11}$, $N_{22}$, $N_{33}$, ${}_s{Q}_{12}$, ${}_s{Q}_{23}$, ${}_s{Q}_{31}$, ${}_a{Q}_{12}$, ${}_a{Q}_{23}$, ${}_a{Q}_{31}$ (9), ${\varphi }_1$ and ${\varphi }_2$ (2), a total of 19.

Remarks

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 ${\varphi }_1$, and ${\varphi }_2$ 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 [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.

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 ${ℝ}^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 $\xi \eta \zeta$ in a two unit cube in which $\xi$ and $\eta$ are in the middle plane of the element and $\zeta$ is the transverse direction to the middle plane. The origin of the coordinate system $\xi \eta \zeta$ is at the center of the two unit cube (Figure 2).

5) We can design p-version hierarchical local approximation in $\xi$, $\eta$ and $\zeta$ of p-levels $p_{\xi }$, $p_{\eta }$ and $p_{\zeta }$ with variable choice of the orders of the approximation space to ensure desired global differentiability of the approximations. In the local approximation, displacements $u_1$, $u_2$ and $u_3$ can be approximated by polynomials of degrees $p_{\xi }$, $p_{\eta }$ and $p_{\zeta }$ as well as higher orders of differentiability. Choices of p-levels $p_{\xi }$, $p_{\eta }$ and $p_{\zeta }$ define approximations in the plane of the plate or the shell element as well as deformation of its cross section(s). By choosing $p_{\xi }$, $p_{\eta }$ and $p_{\zeta }$ 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 $u_1$, $u_2$ and $u_3$ 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.

8) It is important to point out that the formulation of p-version hierarchical $C^0$ plate/shell finite element formulation (for reversible deformation without

(a)(b)(c)(d)(e)(f)

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 $V_3$ defining the top and bottom surfaces at node i in $\xi \mathrm{,}\alpha$ space ( ${\bar{\Omega }}^{\xi \alpha }$ ); (d) Map of the element of figure - (b) in the natural coordinate space $\xi$, $\eta$ and $\zeta$; (e) Nodal configuration in $\zeta$ direction; (f) Map of node i and vector $V_3$.

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.

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 [30] [41] - [57]. Based on the conservation and balance laws (more specifically second law of thermodynamics), the choice of $\Phi$, $\eta$, $\sigma$ and $q$ as constitutive variables at the onset is rather straight forward. Using the conjugate pairs, ${\sigma }_{ij}{\stackrel{\dot{}}{\epsilon }}_{ij}$ and $\frac{q\cdot g}{\theta }$ and thermoelastic deformation, choice of $\epsilon$, $\theta$ as argument tensors of $\sigma$ and $g$ and $\theta$ as argument tensors of $q$ is appropriate. Using the principle of equipresence, we can choose $\epsilon$, $g$ and $\theta$ as argument tensors for $\Phi$, $\eta$. Thus, we have

$\sigma =\sigma \left(\epsilon \mathrm{,}\theta \right)$ (73)

$q=q\left(g\mathrm{,}\theta \right)$ (74)

$\Phi =\Phi \left(\epsilon \mathrm{,}g\mathrm{,}\theta \right)$ (75)

$\eta =\eta \left(\epsilon \mathrm{,}g\mathrm{,}\theta \right)$ (76)

6.1.1. Constitutive Theory of σ

Constitutive theory for Cauchy stress tensor $\sigma$ 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 $\epsilon$ but linear in $\theta$ when based on integrity [30].

1) Constitutive theory for $\sigma$ using Helmholtz free energy density $\Phi$ using (75), we can write

$\frac{D\Phi }{Dt}=\stackrel{\dot{}}{\Phi }=\frac{\partial \Phi }{\partial {\epsilon }_{kl}}{\stackrel{\dot{}}{\epsilon }}_{kl}+\frac{\partial \Phi }{\partial g_i}{\stackrel{\dot{}}{g}}_i+\frac{\partial \Phi }{\partial \theta }\stackrel{\dot{}}{\theta }$ (77)

Substituting (77) in SLT (A5) and rearranging terms, we can write

$\left({\rho }_0\frac{\partial \Phi }{\partial {\epsilon }_{kl}}-{\sigma }_{kl}\right){\stackrel{\dot{}}{\epsilon }}_{kl}+{\rho }_0\left(\eta +\frac{\partial \Phi }{\partial \theta }\right)\stackrel{\dot{}}{\theta }+{\rho }_0\frac{\partial \Phi }{\partial g_i}{g}_i\le 0$ (78)

For (78) to hold for arbitrary but admissible choices of $\stackrel{\dot{}}{\epsilon }$, $\stackrel{\dot{}}{\theta }$ and $\stackrel{\dot{}}{g}$, the following must hold

${\rho }_0\frac{\partial \Phi }{\partial {\epsilon }_{kl}}-{\sigma }_{kl}=0\Rightarrow {\sigma }_{kl}={\rho }_0\frac{\partial \Phi }{\partial {\epsilon }_{kl}}$ (79)

${\rho }_0\left(\eta +\frac{\partial \Phi }{\partial \theta }\right)=0\Rightarrow \eta =-\frac{\partial \Phi }{\partial \theta }$ (80)

${\rho }_0\frac{\partial \Phi }{\partial g_i}{g}_i=0\Rightarrow \Phi \ne \Phi \left(g\right)$ (81)

and SLT (78) reduces to

$\frac{q\cdot g}{\theta }\le 0$ (82)

with (79) - (82), SLT (78) is satisfied.

Remarks

a) Equation (79) can be used to derive constitutive theory for Cauchy stress tensor $\sigma$.

b) Equation (80) implies that $\eta$ is not a constitutive variable as it is defined by $\frac{\partial \Phi }{\partial \theta }$.

c) Equation (82) implies that $\Phi$ is not a function of $g$, thus $\Phi =\Phi \left(\epsilon \mathrm{,}\theta \right)$.

In the following we derive constitutive theory for $\sigma$ using

${\sigma }_{kl}={\rho }_0\frac{\partial \varphi \left(\epsilon \mathrm{,}\theta \right)}{\partial {\epsilon }_{kl}}$ (83)

Since, the constitutive theories must be frame independent, we must consider $\Phi =\Phi \left(I_{\epsilon }\mathrm{,}I{I}_{\epsilon }\mathrm{,}II{I}_{\epsilon }\mathrm{,}\theta \right)$; $I_{\epsilon }\mathrm{,}I{I}_{\epsilon }\mathrm{,}II{I}_{\epsilon }$ being principle invariants of $\epsilon$. Thus (83) becomes

${\sigma }_{kl}={\rho }_0\frac{\partial \varphi \left(I_{\epsilon }\mathrm{,}I{I}_{\epsilon }\mathrm{,}II{I}_{\epsilon }\mathrm{,}\theta \right)}{\partial {\epsilon }_{kl}}$ (84)

or

${\sigma }_{kl}={\rho }_0\left(\frac{\partial \varphi }{\partial I_{\epsilon }}\frac{\partial I_{\epsilon }}{\partial {\epsilon }_{kl}}+\frac{\partial \varphi }{\partial I{I}_{\epsilon }}\frac{\partial I{I}_{\epsilon }}{\partial {\epsilon }_{kl}}+\frac{\partial \varphi }{\partial II{I}_{\epsilon }}\frac{\partial II{I}_{\epsilon }}{\partial {\epsilon }_{kl}}\right)$ (85)

Following the details of similar derivations in reference [30], we can obtain the following from (85) after using the Hamilton-Cayley theorem.

$\left[\sigma \right]={}^{\sigma }\underset{˜}{\alpha }{}^0\left[I\right]+{}^{\sigma }\underset{˜}{\alpha }{}^1\left[\epsilon \right]+{}^{\sigma }\underset{˜}{\alpha }{}^2{\left[\epsilon \right]}^2$ (86)

in which

${}^{\sigma }\underset{˜}{\alpha }{}^i={}^{\sigma }\underset{˜}{\alpha }{}^i\left(I_{\epsilon }\mathrm{,}I{I}_{\epsilon }\mathrm{,}II{I}_{\epsilon }\mathrm{,}\theta \right)\mathrm{<mo>\; ;\; </mo>}i=\mathrm{0,1,2}$ (87)

Material coefficients are determined by expanding ${}^{\sigma }\underset{˜}{\alpha }{}^i\left(\mathrm{.}\right)$; $i=\mathrm{0,1,2}$ in Taylor series in $I_{\epsilon }$, $I{I}_{\epsilon }$, $II{I}_{\epsilon }$ and $\theta$ about a known configuration $\underset{\_}{\Omega }$ and only retaining up to linear terms in the invariants and temperature $\theta$ (for simplicity).

$\begin{array}{c}{}^{\sigma }\underset{˜}{\alpha }{}^i={{}^{\sigma }\underset{˜}{\alpha }{}^i|}_{\underset{\_}{\Omega }}+{}^{\sigma }\underset{˜}{\alpha }{}^i{\frac{\partial {}^{\sigma }\underset{˜}{\alpha }{}^i}{\partial I_{\epsilon }}|}_{\underset{\_}{\Omega }}\left(I_{\epsilon }-{I_{\epsilon }|}_{\underset{\_}{\Omega }}\right)+{}^{\sigma }\underset{˜}{\alpha }{\frac{\partial {}^{\sigma }\underset{˜}{\alpha }{}^i}{\partial I{I}_{\epsilon }}|}_{\underset{\_}{\Omega }}\left(I{I}_{\epsilon }-{I{I}_{\epsilon }|}_{\underset{\_}{\Omega }}\right)\\ +{}^{\sigma }\underset{˜}{\alpha }{}^i{\frac{\partial {}^{\sigma }\underset{˜}{\alpha }{}^i}{\partial II{I}_{\epsilon }}|}_{\underset{\_}{\Omega }}\left(II{I}_{\epsilon }-{II{I}_{\epsilon }|}_{\underset{\_}{\Omega }}\right)+{\frac{\partial {}^{\sigma }\underset{˜}{\alpha }{}^i}{\partial \theta }|}_{\underset{\_}{\Omega }}\left(\theta -{\theta }_{\underset{\_}{\Omega }}\right)\end{array}$ (88)

We substitute (88) in (86) and group the terms and coefficients and if we define

$\begin{array}{l}{}^{\sigma }\underset{˜}{I}{}^1=I_{\epsilon }\mathrm{<mo>\; ;\; </mo>}{}^{\sigma }\underset{˜}{I}{}^2=I{I}_{\epsilon }\mathrm{<mo>\; ;\; </mo>}{}^{\sigma }\underset{˜}{I}{}^3=II{I}_{\epsilon }\\ \text{i}\text{.e}\text{.,}{}^{\sigma }\underset{˜}{I}{}^j\mathrm{<mo>\; ;\; </mo>}j=\mathrm{1,2,}\cdots \mathrm{,}M\mathrm{<mo>\; ;\; </mo>}M=\mathrm{3\; <mo>\; ;\; </mo>}\text{theinvariants}\end{array}$ (89)

and

$\begin{array}{l}\left[{}^{\sigma }{\underset{˜}{G}}^1\right]=\left[\epsilon \right]\mathrm{<mo>\; ;\; </mo>}\left[{}^{\sigma }{\underset{˜}{G}}^2\right]={\left[\epsilon \right]}^2\\ \text{i}\text{.e}\text{.,}\left[{}^{\sigma }{\underset{˜}{G}}^i\right]\mathrm{<mo>\; ;\; </mo>}i=\mathrm{1,2,}\cdots \mathrm{,}N\mathrm{<mo>\; ;\; </mo>}N=\mathrm{2\; <mo>\; ;\; </mo>}\text{thegenerators}\mathrm{.}\end{array}$ (90)

Then (86) can be written as [30]:

$\begin{array}{c}\left[\sigma \right]={{\underset{\_}{\sigma }}^0|}_{\underset{\_}{\Omega }}+{\displaystyle \underset{j=1}{\overset{M}{\sum }}}{}^{\sigma }\underset{\_}{a}{}_j{}^{\sigma }\underset{˜}{I}{}^1\left[I\right]+{\displaystyle \underset{i\mathrm{<mo>\; =\; </mo>1}}{\overset{N}{\sum }}}{}^{\sigma }\underset{\_}{b}{}_i\left[{}^{\sigma }\underset{˜}{G}{}^i\right]+{\displaystyle \underset{i=1}{\overset{N}{\sum }}}{\displaystyle \underset{j=1}{\overset{M}{\sum }}}{}^{\sigma }\underset{˜}{C}{}_{ij}{}^{\sigma }\underset{˜}{I}{}^j\left[{}^{\sigma }\underset{˜}{G}{}^i\right]\\ -{\displaystyle \underset{i=1}{\overset{N}{\sum }}}{}^{\sigma }\underset{\_}{d}{}_i\left(\theta -{\theta |}_{\underset{\_}{\Omega }}\right)\left[{}^{\sigma }\underset{˜}{G}{}^i\right]-{{\underset{\_}{\alpha }}_{tm}|}_{\underset{\_}{\Omega }}\left(\theta -{\theta |}_{\underset{\_}{\Omega }}\right)\left[I\right]\end{array}$ (91)

in which; ${}^{\sigma }\underset{\_}{a}{}_j\mathrm{,}{}^{\sigma }\underset{\_}{b}{}_i\mathrm{,}{}^{\sigma }\underset{˜}{C}{}_{ij}\mathrm{,}{}^{\sigma }\underset{\_}{d}{}_i$; $i=\mathrm{1,2,}\cdots ,N$ and $j=\mathrm{1,2,}\cdots ,M$ are material coefficients that are functions of invariants ${}^{\sigma }\underset{˜}{I}{}^j$; $j=\mathrm{1,2,}\cdots ,M$ i.e., $I_{\epsilon }$, $I{I}_{\epsilon }$, $II{I}_{\epsilon }$ and $\theta$ in the known configuration $\underset{\_}{\Omega }$.

Remarks

1) This constitutive theory for $\sigma$ contains ( $N=2,M=3$ ) fourteen material coefficients and contains up to fifth degree terms of the components of strain tensor $\left[\epsilon \right]$, but is linear in temperature $\theta$.

2) ${{\underset{\_}{\sigma }}^0|}_{\underset{\_}{\Omega }}$ is the initial stress field.

3) A linear constitutive theory in which the products of ${}^{\sigma }\underset{˜}{I}{}^j$, ${}^{\sigma }\underset{˜}{G}{}^i$, $\left(\theta -{\theta |}_{\underset{\_}{\Omega }}\right)$ are neglected and only up to linear terms in $\left[\epsilon \right]$ and $\theta$ are retained is given by

$\left[\sigma \right]={{\underset{\_}{\sigma }}^0|}_{\underset{\_}{\Omega }}\left[I\right]+{}^{\sigma }\underset{˜}{a}{}_1{}^{\sigma }\underset{˜}{I}{}^1\left[I\right]+{}^{\sigma }\underset{˜}{b}{}_1\left[{}^{\sigma }{G}^1\right]-{{\underset{\_}{\alpha }}_{tm}|}_{\underset{\_}{\Omega }}\left(\theta -{\theta |}_{\underset{\_}{\Omega }}\right)\left[I\right]$

using the notation ${}^{\sigma }\underset{˜}{b}{}_1={2\mu |}_{\underset{\_}{\Omega }},{}^{\sigma }\underset{˜}{a}{}_1={\lambda |}_{\underset{\_}{\Omega }}$ (Lame’s constants) and realizing that $\left[{}^{\sigma }{\underset{˜}{G}}^1\right]=\left[\epsilon \right]$ and ${}^{\sigma }\underset{˜}{I}{}^1=\text{tr}\left[\epsilon \right]$ we can write (91) as follows

$\left[\sigma \right]={{\underset{\_}{\sigma }}^0|}_{\underset{\_}{\Omega }}\left[I\right]+{2\mu |}_{\underset{\_}{\Omega }}\left[\epsilon \right]+{\lambda }_{\underset{\_}{\Omega }}\left(\text{tr}\left[\epsilon \right]\right)\left[I\right]-{{\underset{\_}{\alpha }}_{tm}|}_{\underset{\_}{\Omega }}\left(\theta -{\theta |}_{\underset{\_}{\Omega }}\right)\left[I\right]$ (92)

this of course is generalized Hooke’s law.

b) Constitutive theory for $\sigma$ can also be derived using $\sigma =\sigma \left(\epsilon ,\theta \right)$ and the representation theorem [30] [41] - [57], this leads to

$\left[\sigma \right]={}^{\sigma }\underset{˜}{\alpha }{}^0\left[I\right]+{}^{\sigma }\underset{˜}{\alpha }{}^1\left[\epsilon \right]+{}^{\sigma }\underset{˜}{\alpha }{}^2{\left[\epsilon \right]}^2$ (93)

in which $\left[I\right]$, $\left[\epsilon \right]$ and ${\left[\epsilon \right]}^2$ are combined generators of the argument tensors $\left[\epsilon \right]$ and $\theta$ 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 $\Phi$ given above.

6.1.2. Constitutive Theory of q

We begin with $q=q\left(g\mathrm{,}\theta \right)$ based on entropy inequality and use representation theorem [30] [41] - [57] which gives us the following constitutive theory for $q$

$q=-{\kappa |}_{\underset{\_}{\Omega }}g-{{\kappa }_1|}_{\underset{\_}{\Omega }}\left(g\cdot g\right)g-{{\kappa }_2|}_{\underset{\_}{\Omega }}\left(\theta -{\theta |}_{\underset{\_}{\Omega }}\right)g$ (94)

This constitutive theory is based on integrity, hence uses complete basis. From (94), we can also derive a linear constitutive theory

$q=-{\kappa |}_{\underset{\_}{\Omega }}g$ (95)

where $\kappa \mathrm{,}{\kappa }_1\mathrm{,}{\kappa }_2$ are material coefficients defined in a known configuration $\underset{\_}{\Omega }$. These can be functions of ${\left(g\cdot g\right)}_{\underset{\_}{\Omega }}$ (invariant of $g$ ) and ${\theta |}_{\underset{\_}{\Omega }}$.

7. Final Mathematical Model (Thermoelastic)

${\rho }_0\frac{D^2{u}_i}{D{t}^2}-{\rho }_0{F}_i^b-\frac{\partial {\sigma }_{ji}}{\partial x_j}=0\left(\text{BLM}\right)$ (96)

${\sigma }_{ij}={\sigma }_{ji}\left(\text{BAM}\right)$ (97)

${\rho }_0\frac{De}{Dt}+\frac{\partial q_i}{\partial x_i}-{\sigma }_{ij}{\stackrel{\dot{}}{\epsilon }}_{ij}=0\left(\text{FLT}\right)$ (98)

$\begin{array}{c}\left[\sigma \right]={{\underset{\_}{\sigma }}^0|}_{\underset{\_}{\Omega }}+{\displaystyle \underset{j=1}{\overset{M}{\sum }}}{}^{\sigma }\underset{\_}{a}{}_j{}^{\sigma }\underset{˜}{I}{}^1\left[I\right]+{\displaystyle \underset{i\mathrm{<mo>\; =\; </mo>1}}{\overset{N}{\sum }}}{}^{\sigma }\underset{\_}{b}{}_i\left[{}^{\sigma }\underset{˜}{G}{}^i\right]+{\displaystyle \underset{i=1}{\overset{N}{\sum }}}{\displaystyle \underset{j=1}{\overset{M}{\sum }}}{}^{\sigma }\underset{˜}{C}{}_{ij}{}^{\sigma }\underset{˜}{I}{}^j\left[{}^{\sigma }\underset{˜}{G}{}^i\right]\\ -{\displaystyle \underset{i=1}{\overset{N}{\sum }}}{}^{\sigma }\underset{\_}{d}{}_i\left(\theta -{\theta |}_{\underset{\_}{\Omega }}\right)\left[{}^{\sigma }\underset{˜}{G}{}^i\right]-{{\underset{\_}{\alpha }}_{tm}|}_{\underset{\_}{\Omega }}\left(\theta -{\theta |}_{\underset{\_}{\Omega }}\right)\left[I\right]\end{array}$ (99)

$q=-{\kappa |}_{\underset{\_}{\Omega }}g-{{\kappa }_1|}_{\underset{\_}{\Omega }}\left(g\cdot g\right)g-{{\kappa }_2|}_{\underset{\_}{\Omega }}\left(\theta -{\theta |}_{\underset{\_}{\Omega }}\right)g$ (100)

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 which [(a)]

a) The constitutive theory for $\sigma$ 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 are

$-{\rho }_0{F}_i^b-\frac{\partial {\sigma }_{ji}}{\partial x_j}=0$ (101)

${\sigma }_{ij}=2\mu {\epsilon }_{ij}+\lambda {\epsilon }_{kk}{\delta }_{ij}$ (102)

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)).

8. Finite Element Formulation

8.1. Complete Mathematical Model

Expanded form of the mathematical model in ${ℝ}^3$ ((101), (102)) are given in the following:

$\begin{array}{l}{A}_1\left(\mathrm{.}\right)=-\frac{\partial {\sigma }_{11}}{\partial x_1}-\frac{\partial {\sigma }_{21}}{\partial x_2}-\frac{\partial {\sigma }_{31}}{\partial x_3}-{\rho }_0{F}_1^b=0\hfill \\ A_2\left(\mathrm{.}\right)=-\frac{\partial {\sigma }_{12}}{\partial x_1}-\frac{\partial {\sigma }_{22}}{\partial x_2}-\frac{\partial {\sigma }_{32}}{\partial x_3}-{\rho }_0{F}_2^b=0\forall \left(x_1\mathrm{,}{x}_2\mathrm{,}{x}_3\right)\in {\Omega }_x\hfill \\ A_3\left(\mathrm{.}\right)=-\frac{\partial {\sigma }_{13}}{\partial x_1}-\frac{\partial {\sigma }_{23}}{\partial x_2}-\frac{\partial {\sigma }_{33}}{\partial x_3}-{\rho }_0{F}_3^b=0\hfill \end{array}$ (103)

in which ${\bar{\Omega }}_x={\Omega }_x\cup \Gamma$ is the closed spatial domain in ${ℝ}^3$. $\Gamma ={\Gamma }_1\cup {\Gamma }_2$ is the closed boundary of ${\Omega }_x$. The linear constitutive theory for Cauchy stress tensor $\sigma$ and in terms of linear strain is given in the following (using Voigt’s notation)

$\left\{\sigma \right\}=\left[D\right]\left\{ϵ\right\}$ (104)

$\begin{array}{l}{\left\{\sigma \right\}}^{\text{T}}=\left[{\sigma }_{11}\mathrm{,}{\sigma }_{22}\mathrm{,}{\sigma }_{33}\mathrm{,}{\sigma }_{23}\mathrm{,}{\sigma }_{31}\mathrm{,}{\sigma }_{12}\right]\\ {\left\{\epsilon \right\}}^{\text{T}}=\left[{\epsilon }_{11}\mathrm{,}{\epsilon }_{22}\mathrm{,}{\epsilon }_{33}\mathrm{,}{\epsilon }_{23}\mathrm{,}{\epsilon }_{31}\mathrm{,}{\epsilon }_{12}\right]\end{array}$ (105)

$\begin{array}{l}{\epsilon }_{11}=\frac{\partial u_1}{\partial x_1};{\epsilon }_{22}=\frac{\partial u_2}{\partial x_2};{\epsilon }_{33}=\frac{\partial u_3}{\partial x_3}\\ {\epsilon }_{23}=\frac{1}{2}\left(\frac{\partial u_2}{\partial x_3}+\frac{\partial u_3}{\partial x_2}\right);{\epsilon }_{31}=\frac{1}{2}\left(\frac{\partial u_1}{\partial x_3}+\frac{\partial u_3}{\partial x_1}\right);{\epsilon }_{12}=\frac{1}{2}\left(\frac{\partial u_1}{\partial x_2}+\frac{\partial u_2}{\partial x_1}\right)\end{array}$ (106)

The nonzero elements of (6 × 6) [D] material coefficient matrix are given by:

$\begin{array}{l}{D}_{11}=D_{22}=D_{33}=\frac{E\left(1-\nu \right)}{\left(1+\nu \right)\left(1-2\nu \right)}\\ D_{12}=D_{21}=D_{13}=D_{31}=D_{23}=D_{32}=\frac{\nu E}{\left(1+\nu \right)\left(1-2\nu \right)}\\ D_{44}=D_{55}=D_{66}=2G\\ G=\frac{E}{2\left(1+\nu \right)}\end{array}$ (107)

The boundary conditions (BCs) are defined by:

$u_1={\underset{˜}{u}}_1\mathrm{,}{u}_2={\underset{˜}{u}}_2\text{and}{u}_3={\underset{˜}{u}}_3\text{on}{\Gamma }_1$ (108)

$\begin{array}{l}{\sigma }_{11}{n}_{x_1}+{\sigma }_{21}{n}_{x_2}+{\sigma }_{31}{n}_{x_3}=t_1\hfill \\ {\sigma }_{12}{n}_{x_1}+{\sigma }_{22}{n}_{x_2}+{\sigma }_{32}{n}_{x_3}=t_2\hfill \\ {\sigma }_{13}{n}_{x_1}+{\sigma }_{23}{n}_{x_2}+{\sigma }_{33}{n}_{x_3}=t_3\hfill \end{array}\}\text{on}{\Gamma }_2$ (109)

where, $\Gamma ={\Gamma }_1\cup {\Gamma }_2$, a closed boundary of ${\Omega }_x$, $n_{x_1}$, $n_{x_2}$ and $n_{x_3}$ are direction cosines of a unit exterior normal to the boundary ${\Gamma }_2$. E and $\nu$ are modulus of elasticity and Poisson’s ratio. Equations (103) and (104), nine equations in nine unknowns $u_1,u_2,u_3$ and ${\sigma }_{ij}$ constitute the complete mathematical model for the plate or the shell in ${ℝ}^3$. ${\underset{˜}{u}}_1,{\underset{˜}{u}}_2,{\underset{˜}{u}}_3$ are known displacement boundary conditions on ${\Gamma }_1$ and ${\underset{˜}{t}}_1$, ${\underset{˜}{t}}_2$ and ${\underset{˜}{t}}_3$ are known tractions on the boundary ${\Gamma }_2$.

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 $u_1$, $u_2$ and $u_3$. One could easily show that the differential operator $A$ and its adjoint $A^{\mathrm{<mo>\; *\; </mo>}}$ 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) [28] [29], 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 ${\bar{\Omega }}_x^{\text{T}}$ be discretization of ${\bar{\Omega }}_x$, then ${\bar{\Omega }}_x^{\text{T}}={\displaystyle \underset{e}{\cup }{\bar{\Omega }}_x^e}$ in which ${\bar{\Omega }}_x^e$ is a typical element “e”. Let $w_1,w_2,w_3$ be test functions such that

$w_1=\delta u_1\mathrm{<mo>\; ;\; </mo>}{w}_2=\delta u_2\mathrm{<mo>\; ;\; </mo>}{w}_3=\delta u_3\mathrm{.}$ (110)

We construct integral form over ${\bar{\Omega }}_x^{\text{T}}$ using (103) and (110) based on fundamental lemma of the calculus of variations [27] [28] [29].

${\left(A_1(.),w_1\right)}_{{\bar{\Omega }}_x^{\text{T}}}=0,{\left(A_2(.),w_2\right)}_{{\bar{\Omega }}_x^{\text{T}}}=0,{\left(A_3(.),w_3\right)}_{{\bar{\Omega }}_x^{\text{T}}}=0$ (111)

We substitute from (103) in (111) for $A_1(.)$, $A_2(.)$ and $A_3(.)$ and transfer one order of differentiation from the stress derivative terms to the test functions.

Let ${\left(u_1\right)}_h$, ${\left(u_2\right)}_h$ and ${\left(u_3\right)}_h$ be approximations of $u_1\mathrm{,}{u}_2$ and $u_3$ over ${\bar{\Omega }}_x^{\text{T}}$ and ${\left(u_1\right)}_h^e$, ${\left(u_2\right)}_h^e$ and ${\left(u_3\right)}_h^e$ be the approximations of $u_1\mathrm{,}{u}_2$ and $u_3$ over an element “e” with domain ${\bar{\Omega }}_x^e$ such that ${\left(u_i\right)}_h={\displaystyle \underset{e}{\cup }{\left(u_i\right)}_h^e}$ where $i=1,2,3$ are the approximations of $u_1\mathrm{,}{u}_2\mathrm{,}{u}_3$ over ${\bar{\Omega }}_x^{\text{T}}$.

$\begin{array}{l}{\left(A_1(.),w_1\right)}_{{\bar{\Omega }}_x^{\text{T}}}={\displaystyle \underset{e}{\sum }}{\left(A_1^e(.),w_1\right)}_{{\bar{\Omega }}_x^e}=0\\ {\left(A_2(.),w_2\right)}_{{\bar{\Omega }}_x^{\text{T}}}={\displaystyle \underset{e}{\sum }}{\left(A_2^e(.),w_2\right)}_{{\bar{\Omega }}_x^e}=0\\ {\left(A_3(.),w_3\right)}_{{\bar{\Omega }}_x^{\text{T}}}={\displaystyle \underset{e}{\sum }}{\left(A_3^e(.),w_3\right)}_{{\bar{\Omega }}_x^e}=0\end{array}$ (112)

Consider ${\left(A_i^e\left(\mathrm{.}\right)\mathrm{,}{w}_i\right)}_{{\bar{\Omega }}_x^e}$; $i=1,2,3$ over ${\bar{\Omega }}_x^e$

$\begin{array}{l}{\left(A_1^e(.),w_1\right)}_{{\bar{\Omega }}_x^e}={\left(-\frac{\partial {\left({\sigma }_{11}\right)}_h^e}{\partial x_1}-\frac{\partial {\left({\sigma }_{21}\right)}_h^e}{\partial x_2}-\frac{\partial {\left({\sigma }_{31}\right)}_h^e}{\partial x_3}-{\rho }_0{F}_1^b,w_1\right)}_{{\bar{\Omega }}_x^e}\\ {\left(A_2^e(.),w_2\right)}_{{\bar{\Omega }}_x^e}={\left(-\frac{\partial {\left({\sigma }_{12}\right)}_h^e}{\partial x_1}-\frac{\partial {\left({\sigma }_{22}\right)}_h^e}{\partial x_2}-\frac{\partial {\left({\sigma }_{32}\right)}_h^e}{\partial x_3}-{\rho }_0{F}_2^b,w_2\right)}_{{\bar{\Omega }}_x^e}\\ {\left(A_3^e(.),w_3\right)}_{{\bar{\Omega }}_x^e}={\left(-\frac{\partial {\left({\sigma }_{13}\right)}_h^e}{\partial x_1}-\frac{\partial {\left({\sigma }_{23}\right)}_h^e}{\partial x_2}-\frac{\partial {\left({\sigma }_{33}\right)}_h^e}{\partial x_3}-{\rho }_0{F}_3^b,w_3\right)}_{{\bar{\Omega }}_x^e}\end{array}$ (113)

Transfer one order of differentiation from the stresses to the test functions $w_1$, $w_2$ and $w_3$ in (113) using integration by parts (IBP) to obtain

$\begin{array}{l}{\left(A_1^e(.),w_1\right)}_{{\bar{\Omega }}_x^e}={\displaystyle {\int }_{{\bar{\Omega }}_x^e}}\left(\frac{\partial w_1}{\partial x_1}{\left({\sigma }_{11}\right)}_h^e+\frac{\partial w_1}{\partial x_2}{\left({\sigma }_{21}\right)}_h^e+\frac{\partial w_1}{\partial x_3}{\left({\sigma }_{31}\right)}_h^e-\left({\rho }_0{F}_1^b\right)w_1\right)\text{d}\Omega \\ -{\displaystyle {\oint }_{\Gamma }}{w}_1\left({\left({\sigma }_{11}\right)}_h^e{n}_{x_1}+{\left({\sigma }_{21}\right)}_h^e{n}_{x_2}+{\left({\sigma }_{31}\right)}_h^e{n}_{x_3}\right)\text{d}\Gamma \\ {\left(A_2^e(.),w_1\right)}_{{\bar{\Omega }}_x^e}={\displaystyle {\int }_{{\bar{\Omega }}_x^e}}\left(\frac{\partial w_2}{\partial x_1}{\left({\sigma }_{12}\right)}_h^e+\frac{\partial w_2}{\partial x_2}{\left({\sigma }_{22}\right)}_h^e+\frac{\partial w_2}{\partial x_3}{\left({\sigma }_{32}\right)}_h^e-\left({\rho }_0{F}_2^b\right)w_2\right)\text{d}\Omega \\ -{\displaystyle {\oint }_{\Gamma }}{w}_2\left({\left({\sigma }_{12}\right)}_h^e{n}_{x_1}+{\left({\sigma }_{22}\right)}_h^e{n}_{x_2}+{\left({\sigma }_{32}\right)}_h^e{n}_{x_3}\right)\text{d}\Gamma \\ {\left(A_3^e(.),w_1\right)}_{{\bar{\Omega }}_x^e}={\displaystyle {\int }_{{\bar{\Omega }}_x^e}}\left(\frac{\partial w_3}{\partial x_1}{\left({\sigma }_{13}\right)}_h^e+\frac{\partial w_3}{\partial x_2}{\left({\sigma }_{23}\right)}_h^e+\frac{\partial w_3}{\partial x_3}{\left({\sigma }_{33}\right)}_h^e-\left({\rho }_0{F}_3^b\right)w_3\right)\text{d}\Omega \\ -{\displaystyle {\oint }_{\Gamma }}{w}_3\left({\left({\sigma }_{13}\right)}_h^e{n}_{x_1}+{\left({\sigma }_{23}\right)}_h^e{n}_{x_2}+{\left({\sigma }_{33}\right)}_h^e{n}_{x_3}\right)\text{d}\Gamma \end{array}$ (114)

From (114), $u_1\mathrm{,}{u}_2\mathrm{,}{u}_3$ are primary variables (PVs) and the coefficients of $w_1\mathrm{,}{w}_2\mathrm{,}{w}_3$ are secondary variables (SVs). Boundary conditions (108) are essential boundary conditions (EBCs) and those defined by (109) are natural boundary conditions (NBCs). Let

$\begin{array}{l}{\left(u_1\right)}_h^e=\left[N\right]\left\{{\delta }^{u_1}\right\}={\displaystyle \underset{i=1}{\overset{n}{\sum }}}{N}_i{\delta }_i^{u_1}\\ {\left(u_2\right)}_h^e=\left[N\right]\left\{{\delta }^{u_2}\right\}={\displaystyle \underset{i=1}{\overset{n}{\sum }}}{N}_i{\delta }_i^{u_2}\\ {\left(u_3\right)}_h^e=\left[N\right]\left\{{\delta }^{u_3}\right\}={\displaystyle \underset{i=1}{\overset{n}{\sum }}}{N}_i{\delta }_i^{u_3}\end{array}$ (115)

be equal order, equal degree local approximations for $u_1\mathrm{,}{u}_2$ and $u_3$ in which $\left[N\right]$ are approximation functions and $\left\{{\delta }^{u_1}\right\}$, $\left\{{\delta }^{u_2}\right\}$ and $\left\{{\delta }^{u_3}\right\}$ are nodal degrees of freedom for ${\left(u_1\right)}_h^e$, ${\left(u_2\right)}_h^e$ and ${\left(u_3\right)}_h^e$. From (115), we obtain

$\begin{array}{l}{w}_1=\delta \left({\left(u_1\right)}_h^e\right)=N_j\hfill \\ w_2=\delta \left({\left(u_2\right)}_h^e\right)=N_j;j=1,2,3\hfill \\ w_3=\delta \left({\left(u_3\right)}_h^e\right)=N_j\hfill \end{array}$ (116)

First we define secondary variables $P_1^e$, $P_2^e$ and $P_3^e$ as

$\begin{array}{l}{P}_1^e={\left({\sigma }_{11}\right)}_h^e{n}_{x_1}+{\left({\sigma }_{21}\right)}_h^e{n}_{x_2}+{\left({\sigma }_{31}\right)}_h^e{n}_{x_3}\\ P_2^e={\left({\sigma }_{12}\right)}_h^e{n}_{x_1}+{\left({\sigma }_{22}\right)}_h^e{n}_{x_2}+{\left({\sigma }_{32}\right)}_h^e{n}_{x_3}\\ P_3^e={\left({\sigma }_{13}\right)}_h^e{n}_{x_1}+{\left({\sigma }_{23}\right)}_h^e{n}_{x_2}+{\left({\sigma }_{33}\right)}_h^e{n}_{x_3}\end{array}$ (117)

We substitute strains (106) in (104), then (104) in (114) to obtain

$\left\{\begin{array}{c}{A}_1^e\left(\mathrm{.}\right)\\ A_2^e\left(\mathrm{.}\right)\\ A_3^e\left(\mathrm{.}\right)\end{array}\right\}=\left[K^e\right]\left\{{\delta }^e\right\}-\left\{f^e\right\}-\left\{P^e\right\}$ (118)

in which

$\left[K^e\right]={\displaystyle {\int }_{{\bar{\Omega }}^e}}{\left[B\right]}^{\text{T}}\left[D\right]\left[B\right]\text{d}\Omega$ (119)

${\left\{{\delta }^e\right\}}^{\text{T}}=\left[{\left\{{\delta }_1^u\right\}}^{\text{T}}\mathrm{,}{\left\{{\delta }_2^u\right\}}^{\text{T}}\mathrm{,}{\left\{{\delta }_3^u\right\}}^{\text{T}}\right]$ (120)

$\left\{f^e\right\}=\left\{\begin{array}{c}\left\{f^{u_1}\right\}\\ \left\{f^{u_2}\right\}\\ \left\{f^{u_3}\right\}\end{array}\right\}$ (121)

in which

$f_i^{u_1}={\displaystyle \underset{{\bar{\Omega }}_x^e}{\int }}{N}_i{\rho }_0{F}_1^b\text{d}\Omega ;f_i^{u_2}={\displaystyle \underset{{\bar{\Omega }}_x^e}{\int }}{N}_i{\rho }_0{F}_2^b\text{d}\Omega ;f_i^{u_3}={\displaystyle \underset{{\bar{\Omega }}_x^e}{\int }}{N}_i{\rho }_0{F}_3^b\text{d}\Omega ;i=1,2,\cdots ,n$

$\left\{P^e\right\}$ : a vector of secondary variables at the nodes.

The matrix $\left[B\right]$ is defined by

$\left\{\epsilon \right\}=\left[B\right]\left\{{\delta }^e\right\}$ (122)

$\left[K^e\right]$ is the element stiffness matrix, $\left\{f^e\right\}$ 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 [28] [29]. We note that $\left[K^e\right]={\left[K^e\right]}^{\text{T}}$, 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 $\left[K^e\right]$ and $\left\{f^e\right\}$ we need to determine the local approximation functions in $\left[N\right]$ 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 $V_3^i\mathrm{,}i=\mathrm{1,2},\cdots \mathrm{,9}$ connecting nodes $i_{bottom}$ and $i_{top}$. The surface containing the middle points of vectors $V_3^i\mathrm{,}i=\mathrm{1,2},\cdots \mathrm{,9}$ 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 $V_3^i\mathrm{,}i=\mathrm{1,2},\cdots \mathrm{,9}$ 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 $i=\mathrm{1,2},\cdots \mathrm{,9}$ (Figure 2(b)) and the vectors $V_3^i\mathrm{,}i=\mathrm{1,2},\cdots \mathrm{,9}$ completely define the geometry of shell element in ${ℝ}^3$. The shell element of Figure 2(b) is mapped into a natural coordinate space $\xi ,\eta ,\zeta$ in a two unit cube with the origin of the coordinate system at the center of the cube (Figure 2(d)).

Symbolically ${\bar{\Omega }}^e$ is mapped into ${\bar{\Omega }}^m$, ${\bar{\Omega }}^m$ being the map of the element ${\bar{\Omega }}^e$ in the two unit cube. In the map ${\bar{\Omega }}^m$ of the element, the element nodes $i=\mathrm{1,2},\cdots \mathrm{,9}$ are located in $\xi \mathrm{,}\eta$ plane at $\zeta =0$. The natural coordinates $\zeta =-1$ and $\zeta =1$ define bottom and top faces of the plate/shell element. Thus, each $V_3^i$ is mapped into a two unit length in the $\zeta$ direction (perpendicular to $\xi \eta$ plane at $\zeta =0$ ). Let $x_1^i$, $x_2^i$, $x_3^i$; $i=\mathrm{1,2},\cdots \mathrm{,9}$ be the nodal coordinates of the shell element nodes located in $\zeta =0$ plane. Then, the coordinates of an arbitrary point say $P_m\left(x_1\mathrm{,}{x}_2\mathrm{,}{x}_3\right)$ for arbitrary $\xi \mathrm{,}\eta$ but $\zeta =0$ are defined by

$\left\{\begin{array}{c}{x}_1\\ x_2\\ x_3\end{array}\right\}={\displaystyle \underset{i=1}{\overset{9}{\sum }}}{\tilde{N}}_i\left(\xi ,\eta \right){\left\{\begin{array}{c}{x}_1^i\\ x_2^i\\ x_3^i\end{array}\right\}}_{mid}$ (123)

in which $x_j^i={\left(x_j^i\right)}_{mid}$; $j=\mathrm{1,2,3}$ for $i=\mathrm{1,2},\cdots \mathrm{,9}$. The coordinates of a point at an arbitrary location $\xi \mathrm{,}\eta \mathrm{,}\zeta$ with respect to the middle plane ( $\zeta =0$ ) are given by

$\left\{\begin{array}{c}{}^0{x}_1\\ {}^0{x}_2\\ {}^0{x}_3\end{array}\right\}={\displaystyle \underset{i=1}{\overset{9}{\sum }}}{\tilde{N}}_i\left(\xi \mathrm{,}\eta \right)\left(\frac{\zeta }{2}{\left\{\begin{array}{c}{x}_1^i\\ x_2^i\\ x_3^i\end{array}\right\}}_{top}-\frac{\zeta }{2}{\left\{\begin{array}{c}{x}_1^i\\ x_2^i\\ x_3^i\end{array}\right\}}_{bot}\right)={\displaystyle \underset{i=1}{\overset{9}{\sum }}}{\tilde{N}}_i\left(\xi \mathrm{,}\eta \right)\frac{\zeta }{2}{V}_3^i$ (124)

Thus, the coordinates at any arbitrary point $P\left(\xi \mathrm{,}\eta \mathrm{,}\zeta \right)$ in the plate/shell element are given by $x_i=x_i^m+{}^{0}x{}_i$; $i=1,2,3$. Using (123) and (124) we obtain the following for the coordinates $x_i$ at an arbitrary point at a location $\xi \mathrm{,}\eta \mathrm{,}\zeta$ in the natural coordinate map ${\bar{\Omega }}^m$.

$\left\{\begin{array}{c}{x}_1\\ x_2\\ x_3\end{array}\right\}={\displaystyle \underset{i=1}{\overset{9}{\sum }}}{\tilde{N}}_i\left(\xi \mathrm{,}\eta \right)\left(\left\{\begin{array}{c}{x}_1^i\\ x_2^i\\ x_3^i\end{array}\right\}+\frac{\zeta }{2}{V}_3^i\right)$ (125)

The local approximation functions for the shell element are derived using tensor product of 1D functions in $\xi \mathrm{,}\eta$ and $\zeta$ 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 $\xi \mathrm{,}\eta$ and $\zeta$.

For the sake of simplicity we illustrate the process of deriving details of element local approximation using approximations of class $C^0\left({\bar{\Omega }}^e\right)$ for $u_1$, $u_2$ and $u_3$. Such approximations are in a scalar product space $V\subset H^{k\mathrm{,}p}\left({\bar{\Omega }}^e\right)$; $k=1$, with p-levels $p_{\xi }$, $p_{\eta }$ and $p_{\zeta }$. Extensions of the approximations to higher classes defining local approximations of higher order global differentiability can be considered based on ref [28] [29] [58] - [63]. To illustrate the details, we consider $\psi$ to be the dependent variable for which we wish to establish local approximation functions $\left[N\left(\xi \mathrm{,}\eta \mathrm{,}\zeta \right)\right]$ such that

${\psi }_h^e\left(\xi \mathrm{,}\eta \mathrm{,}\zeta \right)=\left[N\left(\xi \mathrm{,}\eta \mathrm{,}\zeta \right)\right]\left\{{\delta }^{\psi }\right\}$ (126)

The $C^0$ p-version local approximation for $\psi$ in the $\xi$ direction for the three node configuration (with nodes 1, 2 and 3 located at $\xi =-1$, $\xi =1$ and $\xi =1$ can be written as [28] [29].

$\begin{array}{l}\psi \left(\xi \right)=\left(\frac{1-\xi }{2}\right){\psi }_1+\left(\frac{1+\xi }{2}\right){\psi }_3+{\displaystyle \underset{i=1}{\overset{p_{\xi }}{\sum }}}\left(\frac{{\xi }^i-a}{i!}\right){\frac{{\partial }^i\psi }{\partial {\xi }^i}|}_{\xi =0}\\ a=\{\begin{array}{l}1;\text{if}i\text{iseven}\hfill \\ \xi ;\text{if}i\text{isodd}\hfill \end{array}\end{array}$ (127)

Similarly in the $\eta$ direction for the three node configuration (with nodes 1, 2 and 3 located at $\eta =-1$, $\eta =0$ and $\eta =1$ ) we can write the following for $C^0$ p-version hierarchical local approximation.

$\begin{array}{l}\psi \left(\eta \right)=\left(\frac{1-\eta }{2}\right){\psi }_1+\left(\frac{1+\eta }{2}\right){\psi }_3+{\displaystyle \underset{i=1}{\overset{p_{\eta }}{\sum }}}\left(\frac{{\eta }^i-b}{i!}\right){\frac{{\partial }^i\psi }{\partial {\eta }^i}|}_{\eta =0}\\ b=\{\begin{array}{l}\mathrm{1\; <mo>\; ;\; </mo>}\text{if}i\text{iseven}\hfill \\ \eta \mathrm{<mo>\; ;\; </mo>}\text{if}i\text{isodd}\hfill \end{array}\end{array}$ (128)

In the $\zeta$ direction, the shell element nodes are located at $\zeta =0$ (Figure 2(b), Figure 2(d)), hence the Lagrange (or others) interpolations (local approximations) in the $\zeta$ direction for $p_{\zeta }=\mathrm{1,2},\cdots$, etc. corresponding to 2, 3, 4, …, etc. nodal configurations in $\zeta$ direction respectively need to be reduced down to a single node located at $\zeta =0$.

We consider a map of $-1\le \zeta \le 1$ at each node in $\alpha$ coordinate system pointing in the same direction as $\zeta$ in which $\zeta \in \left[-\mathrm{1,1}\right]$ is mapped into $\alpha \in \left[-h_i/2\mathrm{,}{h}_i/2\right]$ at each node (say i) (Figure 2(e)). This mapping is given by

$\alpha =\frac{1-\zeta }{2}\left(\frac{-h_i}{2}\right)+\frac{1+\zeta }{2}\left(\frac{h_i}{2}\right)$ (129)

or

$\zeta =\frac{\alpha }{h_i/2}\mathrm{<mo>\; ;\; </mo>}\left(\text{ateachnode}i\text{oftheelement}\right)$ (130)

in which $h_i=\Vert V_3^i\Vert$, length of $V_3^i$ at node i of the shell element. Derivations of the one dimensional p-version hierarchical functions in $\zeta$ or $\alpha$ direction is given in the following. Consider a typical plate/shell element node i ( $i=\mathrm{1,2},\cdots \mathrm{,9}$ ) shown in Figure 2(c), Figure 2(d), Figure 2(f). For p-level of one ( $p_{\zeta }=1$ ) in the $\zeta$ direction the Lagrange interpolation for the two node configuration of Figure 2(e) can be written as

$\psi \left(\zeta \right)=\left(\frac{{\psi }_2+{\psi }_1}{2}\right)+\zeta \left(\frac{{\psi }_2-{\psi }_1}{2}\right)$ (131)

Substituting from (130) in (131)

$\psi \left(\alpha \right)=\left(\frac{{\psi }_2+{\psi }_1}{2}\right)+\frac{\alpha }{h_i/2}\left(\frac{{\psi }_2-{\psi }_1}{2}\right)$ (132)

$\frac{\partial \psi }{\partial \alpha }=\frac{1}{h_i/2}\left(\frac{{\psi }_2-{\psi }_1}{2}\right)={\frac{\partial \psi }{\partial \alpha }|}_{\alpha =0}$ (133)

Hence,

$\left(\frac{{\psi }_2-{\psi }_1}{2}\right)=\frac{h_i}{2}{\frac{\partial \psi }{\partial \alpha }|}_{\alpha =0}$ (134)

and

$\left(\frac{{\psi }_2+{\psi }_1}{2}\right)=\psi \text{forthenode}i\text{at}\zeta =0\text{or}\alpha =0$ (135)

using (134) and (135) in (131)

$\psi \left(\zeta \right)=\psi +\left(\frac{\zeta h_i}{2}\right){\frac{\partial \psi }{\partial \alpha }|}_{\alpha =0}$ (136)

Likewise for $p_{\zeta }=2$, for the three node configuration in Figure 1, we can write (Lagrange interpolation)

$\psi \left(\zeta \right)={\psi }_2+\frac{\zeta }{2}\left(\frac{{\psi }_3-{\psi }_1}{2}\right)+\frac{{\zeta }^2}{2}\left(\frac{{\psi }_1-2{\psi }_2+{\psi }_3}{2}\right)$ (137)

Substituting for $\zeta$ from (130)

$\psi \left(\zeta \right)={\psi }_2+\frac{\alpha }{h_i}\left(\frac{{\psi }_3-{\psi }_1}{2}\right)+\frac{{\alpha }^2}{{\left(\frac{h}{2}\right)}^2}\frac{1}{2}\left(\frac{{\psi }_1-2{\psi }_2+{\psi }_3}{2}\right)$ (138)

differentiating $\psi \left(\alpha \right)$ with respect to $\alpha$ twice and evaluating the derivatives at $\alpha =0$ and substituting these back in (138) and then substituting for $\alpha$ from (130) we obtain the following for the shell node i

$\psi \left(\zeta \right)=\psi +\left(\frac{\zeta h_i}{2}\right){\frac{\partial \psi }{\partial \alpha }|}_{\alpha =0}+{\left(\frac{\zeta h_i}{2}\right)}^2\frac{1}{2!}{\frac{{\partial }^2\psi }{\partial {\alpha }^2}|}_{\alpha =0}$ (139)

Following the same procedure, we can derive the following for a node i located at $\zeta =0$ for p-level of $p_{\zeta }$.

$\psi \left(\zeta \right)=\psi +\left(\frac{\zeta h_i}{2}\right){\frac{\partial \psi }{\partial \alpha }|}_{\alpha =0}+{\left(\frac{\zeta h_i}{2}\right)}^2\frac{1}{2!}{\frac{{\partial }^2\psi }{\partial {\alpha }^2}|}_{\alpha =0}+\cdots +{\left(\frac{\zeta h_i}{2}\right)}^{p_{\zeta }}\frac{1}{p_{\zeta }!}{\frac{\partial {}^{p_{\zeta }}\psi }{\partial {\alpha }^{p_{\zeta }}}|}_{\alpha =0}$ (140)

we note that $\psi \left(\zeta \right)$ 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) [28] [29].

In $\xi$ direction: for nodes 1, 2 and 3

Nodal variable operators: $\mathrm{1\; <mo>\; ;\; </mo>}\frac{{\partial }^2}{\partial {\xi }^2}\mathrm{,}\cdots \mathrm{,}\frac{{\partial }^{p_{\xi }}}{\partial {\xi }^{p_{\xi }}}\mathrm{<mo>\; ;\; </mo>1}$

Approximation functions: $\frac{1-\xi }{2};\frac{{\xi }^i-a}{i!};i=2,3,\cdots ,p_{\xi };\frac{1+\xi }{2}$

In $\eta$ direction: for nodes 1, 2 and 3

Nodal variable operators: $\mathrm{1\; <mo>\; ;\; </mo>}\frac{{\partial }^2}{\partial {\eta }^2}\mathrm{,}\cdots \mathrm{,}\frac{{\partial }^{p_{\eta }}}{\partial {\eta }^{p_{\eta }}}\mathrm{<mo>\; ;\; </mo>1}$

Approximation functions: $\frac{1-\eta }{2}\mathrm{<mo>\; ;\; </mo>}\frac{{\eta }^i-b}{i\mathrm{<mo>\; !\; </mo>}}\mathrm{<mo>\; ;\; </mo>}i=\mathrm{2,3,}\cdots \mathrm{,}{p}_{\eta }\mathrm{<mo>\; ;\; </mo>}\frac{1+\eta }{2}$

In $\zeta$ direction: for node i at $\zeta =0$ or $\alpha =0$

Nodal variable operators: $1,{\frac{\partial }{\partial \alpha }|}_{\alpha =0},{\frac{{\partial }^2}{\partial {\alpha }^2}|}_{\alpha =0},\cdots ,{\frac{{\partial }^{p_{\zeta }}}{\partial {\alpha }^{p_{\zeta }}}|}_{\alpha =0}$

Approximation functions: $\mathrm{1,}\frac{\zeta h_i}{2}\mathrm{,}{\left(\frac{\zeta h_i}{2}\right)}^2\frac{1}{\mathrm{2\; <mo>\; !\; </mo>}}\mathrm{,}\cdots \mathrm{,}{\left(\frac{\zeta h_i}{2}\right)}^{p_{\zeta }}\frac{1}{p_{\zeta }\mathrm{<mo>\; !\; </mo>}}$

By taking tensor products of the 1D nodal variable operators in $\xi \mathrm{,}\eta$ and $\zeta$ and letting them act on $\psi$ we obtain the degrees of freedom for ${\psi }_h^e\left(\xi \mathrm{,}\eta \mathrm{,}\zeta \right)$ or $\psi \left(\xi \mathrm{,}\eta \mathrm{,}\zeta \right)$ over ${\bar{\Omega }}^m$ (hence over ${\bar{\Omega }}^e$ ) and the corresponding approximation functions are obtained by taking the tensor product of the 1D approximation functions in $\xi \mathrm{,}\eta$ and $\zeta$ and we can write

${\psi }_h^e\left(\xi \mathrm{,}\eta \mathrm{,}\zeta \right)=\left[N\left(\xi \mathrm{,}\eta \mathrm{,}\zeta \right)\right]\left\{{\delta }^{\psi }\right\}$ (141)

in which $\left\{{\delta }^{\psi }\right\}$ are the degrees of freedom. Using (141) for $u_1,u_2$ and $u_3$ we can write the following for local approximation for ${\left(u_1\right)}_h^e,{\left(u_2\right)}_h^e$ and ${\left(u_2\right)}_h^e$ of $u_1,u_2$ and $u_3$ (assuming $p_{\xi },p_{\eta }$ and $p_{\zeta }$ to be same for $u_1,u_2$ and $u_3$ ).

$\left\{{\varphi }_h^e\right\}=\left\{\begin{array}{c}{\left(u_1\right)}_h^e\\ {\left(u_2\right)}_h^e\\ {\left(u_2\right)}_h^e\end{array}\right\}=\left[\begin{array}{ccc}\left[N\left(\xi \mathrm{,}\eta \mathrm{,}\zeta \right)\right]& \left[0\right]& \left[0\right]\\ \left[0\right]& \left[N\left(\xi \mathrm{,}\eta \mathrm{,}\zeta \right)\right]& \left[0\right]\\ \left[0\right]& \left[0\right]& \left[N\left(\xi \mathrm{,}\eta \mathrm{,}\zeta \right)\right]\end{array}\right]\left\{\begin{array}{c}\left\{{\delta }^{u_1}\right\}\\ \left\{{\delta }^{u_2}\right\}\\ \left\{{\delta }^{u_3}\right\}\end{array}\right\}$ (142)

In which $\left\{{\delta }^{u_1}\right\}$ $\left\{{\delta }^{u_2}\right\}$ $\left\{{\delta }^{u_3}\right\}$ are the degrees of freedom for $u_1,u_2$ and $u_3$ for an element e.

From the geometry defined by (125), we define the Jacobian of mapping $\left[J_c\right]$ between ${\bar{\Omega }}_x^e$ and master element ${\bar{\Omega }}^m$ (using Murnaghan’s notation) [28] [29].

$\left[J_c\right]=\left[\frac{x_1\mathrm{,}{x}_2\mathrm{,}{x}_3}{\xi \mathrm{,}\eta \mathrm{,}\zeta }\right]$ (143)

The volumes in $x_1\mathrm{,}{x}_2\mathrm{,}{x}_3$ and $\xi \mathrm{,}\eta \mathrm{,}\zeta$ frames and the derivatives of the approximation functions are related through

$\text{d}x\text{d}y\text{d}z=\det \left[J_c\right]\text{d}\xi \text{d}\eta \text{d}\zeta$ (144)

$\left\{\begin{array}{c}\frac{\partial N_i}{\partial x_1}\\ \frac{\partial N_i}{\partial x_2}\\ \frac{\partial N_i}{\partial x_3}\end{array}\right\}={\left[J_c^{\text{T}}\right]}^{-1}\left\{\begin{array}{c}\frac{\partial N_i}{\partial \xi }\\ \frac{\partial N_i}{\partial \eta }\\ \frac{\partial N_i}{\partial \zeta }\end{array}\right\};i=1,2,\cdots ,n$ (145)

Using (145), $\left[B\right]$ matrix in (122) is completely defined. Coefficients of $\left[K^e\right]$ are calculated using Gauss quadrature with element map ${\bar{\Omega }}^m$. Details can be found in references [28] [29]. Assembly of the element equations and their solutions follows standard procedure [28] [29].

Remarks

1) In this formulation p-levels in $\xi \mathrm{,}\eta$ and $\zeta$ 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 $\xi \mathrm{,}\eta$ and $\zeta$ result in higher degree approximations of $u_1,u_2$ and $u_3$. While p-levels $p_{\xi }$ and $p_{\eta }$ primarily control the deformation in the plane, $p_{\zeta }$ controls transverse deformation.

3) With $C^0\left({\bar{\Omega }}_x^e\right)$ local approximation described here, the integral over ${\bar{\Omega }}_x^{\text{T}}$ are in Lebesgue sense. Smoothness of the analytical solution considered in the model problems ensures convergence of $C^0\left({\bar{\Omega }}_x^e\right)$ solution to class $C^1\left({\bar{\Omega }}_x^e\right)$ 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^{\mathrm{<mo>\; *\; </mo>}}=A$ ) results in more comprehensive mathematical model as presented here.

9. Solutions of Model Problems

We consider a square plate ( $a\times a$ ) simply supported on all four boundaries. Midplane of the plate lies in $x_1{x}_2$ plane and the origin of fixed x-frame is at the bottom left corner. Plate thickness is “h”. $x_3$ is normal to $x_1{x}_2$ - plane pointing upwards. Top surface of the plate ( $x_3=h/2$ ) is subjected to uniformly distributed load acting in the negative $x_3$ direction. A schematic of the plate is shown in Figure 3. We nondimensionalize density $\hat{\rho }$, modulus $\hat{E}$ and distance (or displacement) by using reference density ${\rho }_0=\hat{\rho }$, reference modulus $E_0=\hat{E}$ and reference length $L_0=1\text{'}\text{'}$. Reference velocity $v_0$ and reference time $t_0$ are $v_0=\sqrt{E_0/{\rho }_0}$ and $t_0=L_0/v_0$. 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 ( $x_1=x_2=10,x_3=0$ ) of 0.1 for all thicknesses using CPT. This gives us the following values of q for thicknesses of $h=\mathrm{0.1,1.0,5.0}$ 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

$\hat{E}=30\times {10}^6\text{psi}$ $E_0=\hat{E}$ $E=1$

$\nu =0.3$ ${\rho }_0=\hat{\rho }$ $\rho =1$

$\hat{\rho }=0.289018\text{lbm}/{\text{in}}^{\text{3}}$ ${Ł}_0=1$ $L=\hat{L}$

$\begin{array}{l}\{\begin{array}{l}a\times b\times h\hfill \\ 20\times 20\times 0.1\hfill \\ 20\times 20\times 1.0\hfill \\ 20\times 20\times 5.0\hfill \\ 20\times 20\times 10.0\hfill \end{array}\hfill \end{array}$

In the numerical studies we choose $\hat{a}=20\text{'}\text{'}$ hence $a=20$. We consider a 4 × 4 uniform discretization for the entire plate. We consider four different thicknesses: $h=\mathrm{0.1,1.0,5.0}$ and 10.0. Due to the choice of $L_0=1\text{'}\text{'}$, 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 $\xi \mathrm{,}\eta$ and $\zeta$ direction ( $p_{\xi }\mathrm{,}{p}_{\eta }$ and $p_{\zeta }$ ) to obtain converged displacements, strains and stresses. The local approximations are considered in $H^{k\mathrm{,}p}\left({\bar{\Omega }}_x^e\right)$ spaces with $k=1$, hence solutions of class $C^0$. For this choice of local approximations, the integrals over ${\bar{\Omega }}_x^{\text{T}}$ (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 $C^0$ to converge to $C^1$ 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 $u_3$ of the centerline ( $y=10$ ) versus axial distance $x_1$ for $h=\mathrm{0.1,1,5,10}$. (a) Displacement $u_3$ versus $x_1$ ( $L\times b\times h=20\times 20\times 0.1$ ); (b) Displacement $u_3$ versus $x_1$ ( $L\times b\times h=20\times 20\times 1$ ); (c) Displacement $u_3$ versus $x_1$ ( $L\times b\times h=20\times 20\times 5$ ); (d) Displacement $u_3$ versus $x_1$ ( $L\times b\times h=20\times 20\times 10$ ).

Figure 5. Mating element edges at $x_1=x_2=5$.

(a) (b) (c) (d)

Figure 6. Shear stress ${}_s{\sigma }_{13}$ versus distance $x_3$ at $x_1=x_2=5.0"$ for $h=\mathrm{0.1,1,5,10}$. (a) Shear stress ${}_s{\sigma }_{13}$ versus $x_3$ at $x_1=x_2=5.0$ ( $L\times b\times h=20\times 20\times 0.1$ ); (b) Shear stress ${}_s{\sigma }_{13}$ versus $x_3$ at $x_1=x_2=5.0$ ( $L\times b\times h=20\times 20\times 1$ ); (c) Shear stress ${}_s{\sigma }_{13}$ versus $x_3$ at $x_1=x_2=5.0$ ( $L\times b\times h=20\times 20\times 5.0$ ); (d) Shear stress ${}_s{\sigma }_{13}$ versus $x_3$ at $x_1=x_2=5.0$ ( $L\times b\times h=20\times 20\times 10.0$ ).

(a) (b) (a) (b)

Figure 7. Normal stress ${}_s{\sigma }_{33}$ versus distance $x_3$ at $x_1=x_2=10.0$ for $h=\mathrm{0.1,1,5,10}$. (a) Normal stress ${}_s{\sigma }_{13}$ versus $x_3$ at $x_1=x_2=10.0$ ( $L\times b\times h=20\times 20\times 0.1$ ); (b) Normal stress ${}_s{\sigma }_{13}$ versus $x_3$ at $x_1=x_2=10.0$ ( $L\times b\times h=20\times 20\times 1$ ); (c) Normal stress ${}_s{\sigma }_{13}$ versus $x_3$ at $x_1=x_2=10.0$ ( $L\times b\times h=20\times 0.20\times 5.0$ ); (d) Normal stress ${}_s{\sigma }_{13}$ versus $x_3$ at $x_1=x_2=10.0$ ( $L\times b\times h=20\times 20\times 10.0$ ).

Discussion of Results

Figures 4(a)-(d) show plots of transverse displacement $u_3$ of the middle plane of the plate as a function of $x_1$ at $x_2=10$ for $h=\mathrm{0.1,1.0,5.0}$ and 10.0. In the case of new formulation, $p_{\xi }=p_{\eta }=3$ and $p_{\zeta }=3\text{-}14$ yield virtually the same displacement $u_3$ (shown in figures). In case of $h=0.1$, 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 $u_3$ versus $x_1$ for the new formulation are in quite good agreement with those from CPT and FSDT. For $h=1.0$, CPT under estimates $u_3$ (Figure 4(b)) due to lack of shear deformation. For this thickness ( $h=1.0$ ), FSDT results are quite reasonable as the shear deformation contribution is not very significant. The converged $u_3$ versus $x_1$ from the new formulation and those from FSDT are in good agreement (Figure 4(b)). For $h=5$ and $h=10$, the shear deformation is progressively more dominant. Figure 4(c) and Figure 4(d) show $u_3$ versus $x_1$ for CPT, FSDT and the new formulation. In both figures, we observe higher values of $u_3$ versus $x_1$ for FSDT compared to CPT. The new formulation yields higher values of $u_3$ 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 $h=10$ compared to $h=5$. These results from the mathematical models are in accordance with the physics of deformation.

Next we consider the behavior of ${\sigma }_{13}$ at $x_1=x_2=5.0$ as a function of coordinate $x_3$ (obtained using the new formulation). At $x_1=x_2=5$, four vertical edges are coincident. Since the local approximations are of class $C^0$, we expect discontinuity of ${\sigma }_{13}$ (function of displacement gradients) at lower p-level. However, with progressively increasing p-levels we expect convergence of ${\sigma }_{13}$ (weak convergence of $C^0$ solutions to the class $C^1$ ).

Figures 6(a)-(d) show plots of ${\sigma }_{13}$ versus $x_3$ at $x_1=x_2=5.0$ for

$h=\mathrm{0.1,1.0,5.0}$ and 10.0. From Figure 6(a) for $h=0.1$, we note that at $p_{\xi }=p_{\eta }=2$ and $p_{\zeta }=3$, there is a discontinuity of ${\sigma }_{13}$. We have two graphs of ${\sigma }_{13}$ versus $x_3$ even though at $x_1=x_2=5.0$, there are four element edges coincident. Figure 5 shows locations $x_1=x_2=5.0$ with four element edges marked $A_i{B}_i;i=1,2,3,4$ that are coincident. Let ${\left({\sigma }_{13}\right)}_{A_i{B}_i}$ be shear stress along

the edges $A_i{B}_i;i=1,2,3,4$. Since ${\sigma }_{13}=G\left(\frac{\partial u_1}{\partial x_3}+\frac{\partial u_3}{\partial x_1}\right)$, ${\left({\sigma }_{13}\right)}_{A_1{B}_1}\ne {\left({\sigma }_{13}\right)}_{A_2{B}_2}$

due to discontinuity of $\partial u_3/\partial x_1$ even though $\partial u_1/\partial x_3$ is continuous. Similarly ${\left({\sigma }_{13}\right)}_{A_3{B}_3}\ne {\left({\sigma }_{13}\right)}_{A_4{B}_4}$, however ${\left({\sigma }_{13}\right)}_{A_1{B}_1}={\left({\sigma }_{13}\right)}_{A_3{B}_3}$ and ${\left({\sigma }_{13}\right)}_{A_2{B}_2}={\left({\sigma }_{13}\right)}_{A_4{B}_4}$. Thus, at $x_1=x_2=5.0$ graphs of ${\left({\sigma }_{13}\right)}_{A_1{B}_1}$ and ${\left({\sigma }_{13}\right)}_{A_3{B}_3}$ versus $x_3$ are coincident. Likewise graphs of ${\left({\sigma }_{13}\right)}_{A_2{B}_2}$ and ${\left({\sigma }_{13}\right)}_{A_4{B}_4}$ versus $x_3$ are coincident as well but have different values than those of ${\left({\sigma }_{13}\right)}_{A_1{B}_1}$ and ${\left({\sigma }_{13}\right)}_{A_3{B}_3}$. Thus, in Figure 6(a) we only observe two sets of graphs. The discontinuity of $\left({\sigma }_{13}\right)$ is due to discontinuity of $\partial u_3/\partial x_1$ along the edge, a consequence of $C^0$ local approximation for displacements at $p_{\xi }=p_{\eta }=3$ and $p_{\zeta }=3$. At $p_{\xi }=p_{\eta }=4$, $p_{\zeta }=3$, the discontinuity of ${\sigma }_{13}$ diminishes and at $p_{\xi }=p_{\eta }=5$, $p_{\zeta }=3\text{-}7$ and 11-14 ${\sigma }_{13}$ versus $x_3$ is almost converged (in weak sense) and has unique values from all the edges. Graphs of ${\sigma }_{13}$ versus $x_3$ at $x_1=x_2=5.0$ for $h=\mathrm{1.0,5.0}$ and 10.0 shown in Figures 6(b)-(d) exhibit exactly similar behavior at lower p-levels. We observe discontinuity of ${\sigma }_{13}$ but with progressively increasing p-levels, the discontinuity of ${\sigma }_{13}$ diminishes and for p-levels in the range of 5 to 9, we obtain converged solutions for ${\sigma }_{13}$. Higher thicknesses of plate naturally require higher p-levels. ${\sigma }_{23}$ versus $x_3$ at $x_1=x_2=10.0$ behaves in exactly similar fashion; graphs and the details are omitted for the sake of brevity.

Figures 7(a)-(d) show plots of ${\sigma }_{33}$ versus $x_3$ at $x_1=x_2=5.0$ for $h=\mathrm{0.1,1.0,5.0}$ and 10.0. On the top surface of the plate ( $x_3=h/2$ ) a uniform pressure load of ${\sigma }_{33}=-q$ is applied and on the bottom surface of the plate ( $x_3=h/2$ ), ${\sigma }_{33}=0$. This holds regardless of the thickness of the plate. For all thickness values, p-levels $p_{\xi }=p_{\eta }=5\text{-}7$ and $p_{\zeta }=4\text{-}9$ are needed to obtain converged values of ${\sigma }_{33}$, but no difficulties are encountered in doing so. In case of CPT, FSDT and HSDT plate mathematical models, ${\sigma }_{33}=0$. This is obviously non physical especially for thick plates. In case of the new formulation presented here ${\sigma }_{33}$ 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 ${ℝ}^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 $\theta$. 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 $\sigma$ and $q$ if so desired. When the deformation process is isothermal, $\theta$ 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 [28] [29], hence the resulting computations are unconditionally stable. The local approximation for displacements in constructed in $H^{k\mathrm{,}p}\left({\bar{\Omega }}^e\right)$ space with $k=1$ i.e., solutions of class $C^0$ but with higher p-levels (p) in all three directions $\left(\xi \mathrm{,}\eta \mathrm{,}\zeta \right)$. 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 ${ℝ}^3$. 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 ${ℝ}^3$ 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 ${ℝ}^3$. ${\sigma }_{33}$ in all plate/shell mathematical models is generally zero, but this is not physical. The new formulation presented here simulates ${\sigma }_{33}$ 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 ${ℝ}^3$ 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.

Acknowledgements

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.

Appendix

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].

${\rho }_0\left(x\right)=\left|J\right|\rho \left(x\mathrm{,}t\right)\left(\text{CM}\right)$ (A1)

${\rho }_0\frac{D{v}_i}{Dt}-{\rho }_0{F}_i^b-\frac{\partial {\sigma }_{ji}}{\partial x_j}=0\left(\text{BLM}\right)$ (A2)

$\sigma ={\sigma }^{\text{T}}\left(\text{BAM}\right)$ (A3)

${\rho }_0\frac{De}{Dt}+\nabla \cdot q-{\sigma }_{ji}{\stackrel{\dot{}}{\epsilon }}_{ij}=0\left(\text{FLT}\right)$ (A4)

${\rho }_0\left(\frac{D\varphi }{Dt}+\eta \frac{D\theta }{Dt}\right)+\frac{q\cdot g}{\theta }-{\sigma }_{ji}{\stackrel{\dot{}}{\epsilon }}_{ij}\le 0\left(\text{SLT}\right)$ (A5)

In which ${\rho }_0$ is the mass density in the reference configuration, $F_1^b$, $F_2^b$ and $F_3^b$ are body force per unit mass in $x_1,x_2$ and $x_3$ directions, ${\sigma }_{ij}$ is Cauchy stress tensor, ${\epsilon }_{ij}$ is linear strain tensor, e is specific internal energy, $q$ is heat flux vector, $g$ is temperature gradient vector, $\varphi$ is Helmholtz free energy density, $\eta$ is entropy density and $\theta$ is temperature.

From the conjugate pairs ${\sigma }_{ij}\mathrm{,}{\stackrel{\dot{}}{\epsilon }}_{ij}$ and $\frac{q\cdot g}{\theta }$, we can conclude that at the very least the following must hold (thermoelastic solid continua)

$\sigma =\sigma \left(\epsilon \mathrm{,}\theta \right)$ (A6)

$q=q\left(g\mathrm{,}\theta \right)$ (A7)

Choices of the argument tensors for $\varphi$ and $\eta$ are discussed in section 4 in conjunction with the derivation of the constitutive theories for $\sigma$ and $q$.

A.2. Non-Classical Continuum Theory Incorporating Internal Rotations

Displacement gradient tensor $\left[{}^{d}J\right]$ and its decomposition into symmetric and antisymmetric tensors can be written as

$\left[\frac{\partial \left\{u\right\}}{\partial \left\{x\right\}}\right]=\left[{}^{d}J\right]=\left[{}_s{}^{d}J\right]+\left[{}_a{}^{d}J\right]=\left[\epsilon \right]+\left[{}_a{}^{i}r\right]$ (A8)

in which $\left[{}_s{}^{d}J\right]$ and $\left[{}_a{}^{d}J\right]$ are symmetric and antisymmetric tensor, thus $\left[\epsilon \right]$ and $\left[{}_a{}^{i}r\right]$ are strain and rotation tensors. $\left[{}_a{}^{i}r\right]$ contains rotations ${}_i{\Theta }_{x_1}$, ${}_i{\Theta }_{x_2}$ and ${}_i{\Theta }_{x_3}$ about $x_1$, $x_2$ and $x_3$ axes. Alternatively

$\nabla \times u=e_1\left({}_i{\Theta }_{x_1}\right)+e_2\left({}_i{\Theta }_{x_2}\right)+e_3\left({}_i{\Theta }_{x_3}\right)$ (A9)

with this definitions of ${}_i{\Theta }_{x_1}$, ${}_i{\Theta }_{x_2}$ and ${}_i{\Theta }_{x_3}$ in (9), we have

$\begin{array}{l}{}_a{}^{i}r{}_{12}={}_i\Theta {}_{x_3};{}_a{}^{i}r{}_{13}=-{}_i\Theta {}_{x_2};{}_a{}^{i}r{}_{23}={}_i\Theta {}_{x_1}\\ {}_a{}^{i}r{}_{21}=-{}_a{}^{i}r{}_{12};{}_a{}^{i}r{}_{31}={-}_a^i\text{}{r}_{13};{}_a{}^{i}r{}_{32}=-{}_a{}^{i}r{}_{23}\end{array}$ (A10)

all others are zero. The rotations ${}_i{\Theta }_{x_1}$, ${}_i{\Theta }_{x_2}$ and ${}_i{\Theta }_{x_3}$ are about the axes of the triad located at a material point. If we define the rotations as a vector $\left\{{}_i\Theta \right\}$, ${\left\{{}_i\Theta \right\}}^{\text{T}}=\left[{}_i\Theta {}_{x_1},{}_i\Theta {}_{x_2},{}_i\Theta {}_{x_3}\right]$, the gradient of $\left\{{}_i\Theta \right\}$ i.e. $\left[{}^{{}_i\Theta }J\right]$ and its decomposition into symmetric and skew-symmetric tensor $\left[{}_s{}^{{}_i\Theta }J\right]$ and $\left[{}_a{}^{{}_i\Theta }J\right]$ can be written as

$\left[\frac{\partial \left\{{}_i\Theta \right\}}{\partial \left\{x\right\}}\right]=\left[{}^{{}_i\Theta }J\right]=\left[{}_s{}^{{}_i\Theta }J\right]+\left[{}_a{}^{{}_i\Theta }J\right]$ (A11)

We have the following conservation and balance laws [31] [32] [36]

${\rho }_0\left(x\right)=\left|J\right|\rho \left(x\mathrm{,}t\right)\left(\text{CM}\right)$ (A12)

${\rho }_0\frac{D{v}_i}{Dt}-{\rho }_0{F}_i^b-\frac{\partial {\sigma }_{ji}}{\partial x_j}=0\left(\text{BLM}\right)$ (A13)

$m_{mk\mathrm{,}m}+{ϵ}_{ijk}{\sigma }_{ij}=0\left(\text{BAM}\right)$ (A14)

${\rho }_0\frac{De}{Dt}+\nabla \cdot q-{}_s\sigma {}_{ij}{\stackrel{\dot{}}{\epsilon }}_{ij}-m_{ij}\left({}_s{}^{{}_i\Theta }J{}_{ij}\right)=0\left(\text{FLT}\right)$ (A15)

${\rho }_0\left(\frac{D\varphi }{Dt}+\eta \frac{D\theta }{Dt}\right)+\frac{q\cdot g}{\theta }{-}_s{\sigma }_{ji}{\stackrel{\dot{}}{\epsilon }}_{ij}-m_{ij}\left({}_s{}^{{}_i\Theta }J{}_{ij}\right)\le 0\left(\text{SLT}\right)$ (A16)

Additionally we have used decomposition of Cauchy stress tensor $\sigma$ into symmetric tensor ${}_s\sigma$ and antisymmetric tensor ${}_a\sigma$

$\sigma ={}_s\sigma +{}_a\sigma$ (A17)

The Cauchy moment tensor $m$ is symmetric due to balance of moment of moments balance law, an additional balance law needed in non-classical continuum theories [31] [32] [36] 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)

${}_s\sigma ={}_s\sigma \left(\epsilon \mathrm{,}\theta \right)$ (A18)

$m=m\left({}_s{}^{{}_i\Theta }J\mathrm{,}\theta \right)$ (A19)

$q=q\left(g\mathrm{,}\theta \right)$ (A20)

Choice of the argument tensors for $\varphi$ and $\eta$ can be based on the principle of equipresence [30]. The derivation of the constitutive theory for ${}_s\sigma$ is same as for $\sigma$ in case of CCM (section). A linear constitutive theory for $m$ can be written as

$m_{ij}=2\underset{˜}{\mu }\left({}_s{}^{{}_i\Theta }J{}_{ij}\right)$ (A21)

where $\underset{˜}{\mu }$ 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 ${}_e\Theta$ be external or Cosserat rotations (unknown) about the axes of the same triad at a material point about which internal rotations ${}_i\Theta$ act, then the total rotations ${}_t\Theta$ are given by

${}_t\Theta ={}_i\Theta +{}_e\Theta$ (A22)

and

$\left[{}_a{}^{t}r\right]=\left[{}_a{}^{i}r\right]+\left[{}_a{}^{e}r\right]$ (A23)

Gradient of ${}_t\Theta$, ${}^{{}_t\Theta }J$ and its decomposition into symmetric and antisymmetric tensors gives

$\left[{}^{{}_t\Theta }J\right]=\left[\frac{\partial \left\{{}_t\Theta \right\}}{\partial \left\{x\right\}}\right]=\left[{}_s{}^{{}_t\Theta }J\right]+\left[{}_a{}^{{}_t\Theta }J\right]=\left[{}_s{}^{{}_t\Theta }J\right]+\left[{}_a{}^{i}r\right]$ (A24)

$\left[{}_s{}^{{}_t\Theta }J\right]=\frac{1}{2}\left(\left[{}^{{}_t\Theta }J\right]+{\left[{}^{{}_t\Theta }J\right]}^{\text{T}}\right)$ (A25)

$\left[{}_a{}^{{}_t\Theta }J\right]=\frac{1}{2}\left(\left[{}^{{}_t\Theta }J\right]-{\left[{}^{{}_t\Theta }J\right]}^{\text{T}}\right)$ (A26)

The CM, BLM, BAM and BMM balance laws in this case are same as in Section A.2. The FLT and the SLT [33] [34] are given by

${\rho }_0\frac{De}{Dt}+\nabla \cdot q-{}_s\sigma {}_{ij}{\stackrel{\dot{}}{\epsilon }}_{ij}-{}_a\sigma {}_{ij}\left({}_a{}^t\stackrel{\dot{}}{r}{}_{ij}\right)-m_{ij}\left({}_s{}^{{}_i\Theta }J{}_{ij}\right)=0\left(\text{FLT}\right)$ (A27)

${\rho }_0\left(\frac{D\varphi }{Dt}+\eta \frac{D\theta }{Dt}\right)+\frac{q\cdot g}{\theta }-{}_s\sigma {}_{ji}{\stackrel{\dot{}}{\epsilon }}_{ij}-{}_a\sigma {}_{ij}\left({}_a{}^t\stackrel{\dot{}}{r}{}_{ij}\right)-m_{ij}\left({}_s{}^{{}_i\Theta }J{}_{ij}\right)\le 0\left(\text{SLT}\right)$ (A28)

The decomposition of $\sigma ={}_s\sigma +{}_a\sigma$ (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 [37] [38]). From the conjugate pairs in the entropy inequality (28), at the very least the following must hold (thermoelastic solid continua)

${}_s\sigma ={}_s\sigma \left(\epsilon \mathrm{,}\theta \right)$ (A29)

$m=m\left({}_s{}^{{}_t\Theta }J\mathrm{,}\theta \right)$ (A30)

${}_a\sigma ={}_a\sigma \left({}_a{}^{t}r\mathrm{,}\theta \right)$ (A31)

Argument tensors for $\Phi$ and $\eta$ at this stage can be established using principle of equipresence. Constitutive theories for $m$ and ${}_a\sigma$ in the absence of $\theta$ reduce to [33] [34]

$m=2\underset{˜}{\mu }\left({}_s{}^{{}_t\Theta }J\right)+\underset{˜}{\lambda }\left(\text{tr}\left({}_s{}^{{}_t\Theta }J\right)I\right)$ (A32)

${}_a\sigma =2\underset{˜}{\beta }\left({}_a{}^{t}r\right)$ (A33)

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Ballarini, R. (2003) The Da Vinci-Euler-Bernoulli Beam Theory? http://www.memagazine.org/contents/current/webonly/webex418.html
[2]	Witmer, E.A. (1991-1992) Elementary Bernoulli-Euler Beam Theory. MIT Unified Engineering Course Notes, 5-114 to 5-164.
[3]	Caresta, M. Vibrations of a Free-Free Beam. https://www.colorado.edu/physics/phys1240/phys1240_fa15/homelabs/Beam_vibration_Mauro %20Caresta_UNS
[4]	Han, S.M., Benaroya, H. and Wei, T. (1999) Dynamics of Transversely Vibrating Beams Using Four Engineering Theories. Final Version, Academic Press, Cambridge. https://doi.org/10.1006/jsvi.1999.2257 https://web.archive.org/web/20060623063248/http://www.memagazine.org/contents/current /webonly/webex418.html
[5]	Timoshenko, S.P. (1953) History of Strength of Materials. McGraw-Hill, New York.
[6]	Chaldni, E.F.F. (1802) Die Akustik. Breitkopf & Härtel, Leipzig.
[7]	Germain, S. (1826) Remarques sur la nature, les bornes et l’étendue de la question des surfaces élastiques et équation générale de ces surfaces. Huzard-Courcier.
[8]	Lagrange, J.L. (1828) Annali di Chimica, 39, 149-207.
[9]	Cauchy, A.L. (1828) Sur l’équilibre et le mouvement d’une plaque solide. Exercises de Mathématiques, 3, 328-355.
[10]	Poisson, S.D. (1829) Mémoire sur l’équilibre et le Mouvement des Corps élastiques. Mémoires de l’Academie des Sciences de l’Institut de France, 8, 257-370.
[11]	Kirchhoff, G. (1850) Uber das Gleichgewicht und die Bewegung einer elastischen Scheibe. Journal für die reine und angewandte Mathematik, ger, 40, 51-88. http://eudml.org/doc/147439 https://doi.org/10.1515/crll.1850.40.51
[12]	Kirchhoff, G. (1850) Uber die uchwingungen einer kriesformigen elastischen scheibe. Annalen der Physik, 81, 58-264. https://doi.org/10.1002/andp.18501571005
[13]	Reddy, J.N. (2004) Mechanics of Laminated Plates and Shells: Theory and Analysis. 2nd Edition, CRC Press, Boca Raton. https://doi.org/10.1201/b12409
[14]	Reddy, J.N. (2007) Theory and Analysis of Elastic Plates and Shells. 2nd Edition, CRC Press, Boca Raton. https://doi.org/10.1201/9780849384165
[15]	Lévy, S. (1942) Bending of Rectangular Plates with Large Deflections. Report No. 737, NACA.
[16]	Wang, C.T. (1948) Bending of Rectangular Plates with Large Deflections. Report No. 162, NACA.
[17]	Lévy, S. (1942) Square plate with Clamped Edges under Pressure Producing Large Deflections. Report No. 847, NACA.
[18]	Timoshenko, S. (1959) Theory of Plates and Shells. McGraw-Hill, New York.
[19]	Leissa, A.W. (1973) Vibration of Shells. NASA SP 288, Washington DC.
[20]	Dym, C.L. (1974) Introduction to the Theory of Shells. Pergamon, New York.
[21]	Ugural, A. (2009) Stresses in Beams, Plates, and Shells. 3rd Edition, CRC Press, Boca Raton.
[22]	Reissner, E. (1944) On the Theory of Bending of Elastic Plates. Journal of Mathematics and Physics, 23, 184-191. https://doi.org/10.1002/sapm1944231184
[23]	Reissner, E. (1945) The Effect of Transverse Shear Deformation on the Bending of Elastic Plates. ASME Journal of Applied Mechanics, 12, A68-A77.
[24]	Reissner, E. (1947) On Bending of Elastic Plates. Quarterly of Applied Mathematics, 5, 55-68. https://doi.org/10.1090/qam/20440
[25]	Reissner, E. (1976) On the Theory of Transverse Bending of Elastic Plates. International Journal of Solids and Structures, 12, 545-555. https://doi.org/10.1016/0020-7683(76)90001-9
[26]	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. https://doi.org/10.1007/s00161-019-00744-8
[27]	Gelfand, M. and Fomin, S. (2017) Calculus of Variations. Dover Publications, Mineola.
[28]	Surana, K.S. and Reddy, J.N. (2017) The Finite Element Method for Boundary Value Problems: Mathematics and Computations. CRC/Taylor and Francis. https://doi.org/10.1201/b22512
[29]	Surana, K. S. and Reddy, J. N. (2017) The Finite Element Method for Initial Value Problems. CRC/Taylor and Francis, Boca Raton. https://doi.org/10.1201/b22512
[30]	Surana, K.S. (2015) Advanced Mechanics of Continua. CRC/Taylor and Francis, Boca Raton. https://doi.org/10.1201/b17959
[31]	Surana, K.S., Powell, M.J. and Reddy, J.N. (2015) A More Complete Thermodynamic Framework for Solid Continua. Journal of Thermal Engineering, 1, 446-459. https://doi.org/10.18186/jte.17430
[32]	Surana, K.S., Reddy, J.N., Nunez, D. and Powell, M.J. (2015) A Polar Continuum Theory for Solid Continua. International Journal of Engineering Research and Industrial Applications, 8, 77-106.
[33]	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. https://doi.org/10.1007/s00161-017-0554-1
[34]	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. https://doi.org/10.1080/19475411.2018.1530700
[35]	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. https://doi.org/10.1093/acprof:oso/9780199641758.001.0001
[36]	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.
[37]	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. https://doi.org/10.1016/S0020-7683(02)00152-X
[38]	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. https://doi.org/10.1007/s11012-018-0851-1
[39]	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. https://doi.org/10.1016/0895-7177(90)90302-4
[40]	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. https://doi.org/10.1007/BF00370038
[41]	Prager, W. (1945) Strain Hardening under Combined Stresses. Journal of Applied Physics, 16, 837-840. https://doi.org/10.1063/1.1707548
[42]	Reiner, M. (1945) A Mathematical Theory of Dilatancy. American Journal of Mathematics, 67, 350-362. https://doi.org/10.2307/2371950
[43]	Rivlin, R.S. and Ericksen, J.L. (1955) Stress-Deformation Relations for Isotropic Materials. Journal of Rational Mechanics and Analysis, 4, 323-425. https://doi.org/10.1512/iumj.1955.4.54011
[44]	Todd, J.A. (1948) Ternary Quadratic Types. Philosophical Transactions of the Royal Society of London. Series A: Mathematical and Physical Sciences, 241, 399-456. https://doi.org/10.1098/rsta.1948.0025
[45]	Rivlin, R.S. (1955) Further Remarks on the Stress-Deformation Relations for Isotropic Materials. Journal of Rational Mechanics and Analysis, 4, 681-702. https://doi.org/10.1512/iumj.1955.4.54025
[46]	Wang, C.C. (1969) On Representations for Isotropic Functions, Part I. Archive for Rational Mechanics and Analysis, 33, 249-267. https://doi.org/10.1007/BF00281278
[47]	Wang, C.C. (1969) On Representations for Isotropic Functions, Part II. Archive for Rational Mechanics and Analysis, 33, 268-287. https://doi.org/10.1007/BF00281279
[48]	Wang, C.C. (1970) On Representations for Isotropic Functions, Part I and Part II. Archive for Rational Mechanics and Analysis, 36, 166-223. https://doi.org/10.1007/BF00272241
[49]	Wang, C.C. (1971) Corrigendum to “Representations for Isotropic Functions”. Archive for Rational Mechanics and Analysis, 43, 392-395. https://doi.org/10.1007/BF00252004
[50]	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. https://doi.org/10.1007/BF00272240
[51]	Smith, G.F. (1971) On Isotropic Functions of Symmetric Tensors, Skew-Symmetric Tensors and Vectors. International Journal of Engineering Science, 9, 899-916. https://doi.org/10.1016/0020-7225(71)90023-1
[52]	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. https://doi.org/10.1007/BF00277933
[53]	Spencer, A.J.M. and Rivlin, R.S. (1960) Further Results in the Theory of Matrix Polynomials. Archive for Rational Mechanics and Analysis, 4, 214-230. https://doi.org/10.1007/BF00281388
[54]	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. https://doi.org/10.1016/B978-0-12-240801-4.50008-X
[55]	Boehler, J.P. (1977) On Irreducible Representations for Isotropic Scalar Functions. Journal of Applied Mathematics and Mechanics, 57, 323-327. https://doi.org/10.1002/zamm.19770570608
[56]	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. https://doi.org/10.1016/0020-7225(93)90109-8
[57]	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. https://doi.org/10.1016/0020-7225(93)90005-F
[58]	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. https://doi.org/10.1142/S1465876301000507
[59]	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. https://doi.org/10.1080/15502280802572270
[60]	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. https://doi.org/10.1080/15502280802572262
[61]	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. https://doi.org/10.1142/S1465876302000605
[62]	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. https://doi.org/10.1142/S1465876303002179
[63]	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. https://doi.org/10.1142/S1465876304002307