Finite Element Processes Based on GM/WF in Non-Classical Solid Mechanics

Show more

1. Introduction Literature Review, and Scope of Work

In Lagrangian description of deforming matter, the Jacobian of deformation is a fundamental quantity of the measure of deformation of the solid continua. In general, the Jacobian of deformation varies between material points, i.e. it varies between a material point and its neighbors. Polar decomposition of the Jacobian of deformation at material points into stretch (left of right) and pure rotation shows that if the Jacobian of deformation varies between a material point and its neighbors so do the rotations. We could also consider the decomposition of the displacement gradient tensor into symmetric and skew symmetric tensors. The skew symmetric tensor is a measure of pure rotations while the symmetric tensor is a measure of strains. Strain measures (such as Green’s strain) are purely a function of stretch tensor or alternatively symmetric part of the displacement gradient tensor. In these measures, rotation tensor plays no role.

In non-polar continuum theories, only conjugate stress and strain tensors contribute to the stored energy in the deforming solid continua. Likewise, the dissipation mechanism is purely due to stress tensor and rates of conjugate strain tensor. In such theories, the influence of rotations and the influence of the rates of rotations on the mechanism of energy storage and dissipation is not considered. In the present work, we consider solid continua in which the rotations and the rates of rotations that exist between neighboring material points are resisted by the constitution of the matter, hence result in energy storage and energy dissipation. Thus, the continuum theory used here for solid continua in Lagrangian description incorporates new physics associated with varying internal rotations and their conjugate moments. This physics is completely absent in the currently used continuum theories for isotropic, homogeneous solid continua. As established in the abstract the theory presented here is indeed a polar continuum theory that incorporates internal varying rotations and conjugate moments in the derivation of conservation and balance laws.

The theory used here is a continuum theory in Lagrangian description for polar continuum and should not be confused with micropolar continuum theories [1] - [11] that are designed to accommodate effects at scales smaller than the continuum scale. Micropolar continuum theories require definitions of additional strain measures [6] related to micromechanics. The polar continuum theory used here incorporates standard measures of strains as used currently in non-polar continuum theories. In the polar continuum theory used here, the motivation is to account for the influence of varying rotations at neighboring material points that arise during evolution as these may result in additional energy storage in some solid continua. Polar decomposition of the Jacobian of deformation at neighboring material points clearly substantiates this. An important point to note is that the theory considered here can only account for local rotation effects due to deformation at material points; hence the theory used here is intrinsically a local polar continuum theory, thus cannot account for nonlocal effects.

In the following we present a brief literature review on micropolar theories, nonlocal theories and stress couple theories. A comprehensive treatment of micropolar theories can be found in the works by Eringen [1] - [9] . The concept of couple stresses is presented by Koiter [10] . Balance laws for micromorphic materials are presented in reference [11] . The micropolar theories consider micro deformation due to micro constituents in the continuum. In references [12] [13] [14] by Reddy et al. and reference [15] by Zang et al. nonlocal theories are presented for bending, buckling and vibration of beams, beams with nanocarbon tubes and bending of plates. The nonlocal effects are believed to be incorporated due to the work presented by Eringen [6] in which definition of a nonlocal stress tensor is introduced through integral relationship using the product of macroscopic stress tensor and a distance kernel representing the nonlocal effects. The polar continuum theory for solid continua presented in this paper is strictly local and non-micropolar. The concept of couple stresses was introduced by Voigt in 1881 by assuming a couple or moment per unit area on the oblique plane of the deformed tetrahedron in addition to the stress or force per unit area. Since the introduction of this concept many published works have appeared. We cite some recent works, most of which are related to micropolar stress couple theories. Authors in reference [16] report experimental study of micropolar and couple stress elasticity of compact bones in bending. Conservation integrals in couple stress elasticity are reported in reference [17] . A microstructure- dependent Timoshenko beam model based on modified couple stress theories is reported by Ma et al. [18] . Further account of couple stress theories in conjunction with beams can be found in references [19] [20] [21] . Treatment of rotation gradient dependent strain energy and its specialization to Von Kármán plates and beams can be found in reference [22] . Other accounts of micropolar elasticity and Cosserat modeling of cellular solids can be found in references [23] [24] [25] . We remark that in references [16] - [25] , Lagrangian description is used for solid matter, however the mathematical descriptions are purely derived using strain energy density functional and principle of virtual work. This approach works well for elastic solids in which mechanical deformation is reversible. Extension of these works to thermoviscoelastic solids with and without memory is not possible. In such materials the thermal field and mechanical deformation are coupled due to the fact that the rate of work results in rate of entropy production. In references [26] - [37] various aspects of the kinematics of micropolar theories, stress couple theories, etc. are discussed and presented including some applications to plates and shells.

If the varying rotations and their rates result in energy storage and dissipation, then their energy conjugate moment (shown later in the paper) must exist in the deforming matter. This necessitates the existence of moment (per unit area) on the oblique plane of the deformed tetrahedron. Thus, at the onset, we consider average force per unit area and displacements, and average moment per unit area and the rotations on the oblique plane of the deformed tetrahedron. The work used here [38] - [44] follows a strictly thermodynamic approach using these i.e., for polar solid continua we consider: (i) Conservation of mass and present reasons for not deriving conservation of inertia (ii) Balance of linear momenta (iii) Balance of angular momenta (iv) Balance of moments of moments (or couples) (v) First law of thermodynamics and (vi) Second law of thermodynamics in Lagrangian description in which stress and strain, moment and rotations are energy conjugate pairs. The mathematical description for polar solid continua used here is applicable to polar thermoelastic solids for small deformation and small strain.

In the present work the mathematical model derived by Surana, et al. [38] - [44] in Lagrangian description for thermoelastic polar solids incorporating internal rotations with small strain, small deformation, isothermal case is used to derive balance of linear momenta equations purely in terms of displacements for boundary value problems. 2D plane stress case is used to present details. These equations are then used to construct finite element formulations using GM/WF. Merits and advantages of this approach over least squares finite element formulation based on mathematical model consisting of a first order system of equations are illustrated in terms of formulation details as well as through three plane stress model problems.

2. Mathematical Model

For non-classical elastic solid matter with internal rotation and conjugate moment physics undergoing small deformation and small strain, the mathematical model for BVPs has been presented by Surana et al. [38] - [44] . In present work we assume isothermal deformation process i.e. no entropy production due to mechanical work, hence the mathematical model in Lagrangian description consists of balance of linear momenta, balance of angular momenta, balance of moments of moments (as a balance law or its absence) [38] - [45] and the constitutive theories for: symmetric part of Cauchy stress tensor, symmetric part of Cauchy moment tensor and antisymmetric part of Cauchy moment tensor (if balance of moments of moments is not used as a balance law). We have the following dimensionless form of the mathematical model in ${\mathbb{R}}^{2}$ (neglecting body forces) assuming that balance of moments of moments is not a balance law [45] . Using the decomposition of the Cauchy moment tensor into symmetric and antisymmetric tensors $m={}_{s}m+{}_{a}m$ , we have constitutive theories for ${}_{s}m$ and ${}_{a}m$ . Choosing $x={x}_{1}$ , $y={x}_{2}$ , $u={u}_{1}$ , $v={u}_{2}$ we can write the following for balance laws and the constitutive theories $\forall x\mathrm{,}y\in {\Omega}_{xy}$ .

$\frac{{\partial}_{s}{\sigma}_{xx}}{\partial x}+\frac{{\partial}_{s}{\sigma}_{yx}}{\partial y}+\frac{{\partial}_{a}{\sigma}_{yx}}{\partial y}=0$ (1)

$\frac{{\partial}_{s}{\sigma}_{xy}}{\partial x}+\frac{{\partial}_{s}{\sigma}_{yy}}{\partial y}-\frac{{\partial}_{a}{\sigma}_{yx}}{\partial x}=0$ (2)

$\begin{array}{l}\frac{\partial {m}_{xz}}{\partial x}+\frac{\partial {m}_{yz}}{\partial y}+2\left({}_{a}{\sigma}_{yx}\right)=0\\ {m}_{xz}={}_{s}m{}_{xz}+{}_{a}m{}_{xz}\\ {m}_{yz}={}_{s}m{}_{yz}+{}_{a}m{}_{yz}\end{array}$ (3)

$\begin{array}{l}{}_{s}{\sigma}_{xx}={D}_{11}\frac{\partial u}{\partial x}+{D}_{12}\frac{\partial u}{\partial y}\\ {}_{s}{\sigma}_{yy}={D}_{21}\frac{\partial u}{\partial x}+{D}_{22}\frac{\partial v}{\partial y}\\ {}_{s}{\sigma}_{xy}={D}_{33}\left(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\right)\end{array}$ (4)

${}_{s}{m}_{xz}=\left(\frac{{E}_{0}}{{m}_{0}{L}_{0}}\right)\underset{\u02dc}{\alpha}\frac{\partial \left({}_{i}{\Theta}_{z}\right)}{\partial x}$ (5)

${}_{s}{m}_{yz}=\left(\frac{{E}_{0}}{{m}_{0}{L}_{0}}\right)\underset{\u02dc}{\alpha}\frac{\partial \left({}_{i}{\Theta}_{z}\right)}{\partial y}$ (6)

${}_{a}{m}_{xz}=-\left(\frac{{E}_{0}}{{m}_{0}{L}_{0}}\right)\underset{\u02dc}{\beta}\frac{\partial \left({}_{i}{\Theta}_{z}\right)}{\partial x}$ (7)

${}_{a}{m}_{yz}=-\left(\frac{{E}_{0}}{{m}_{0}{L}_{0}}\right)\underset{\u02dc}{\beta}\frac{\partial \left({}_{i}{\Theta}_{z}\right)}{\partial y}$ (8)

${}_{i}{\Theta}_{z}=\frac{1}{2}\left(\frac{\partial u}{\partial y}-\frac{\partial v}{\partial x}\right)$ (9)

${D}_{11}={D}_{22}=\frac{E}{1-{\nu}^{2}};\text{\hspace{1em}}{D}_{12}={D}_{21}=\frac{\nu E}{1-{\nu}^{2}};\text{\hspace{1em}}{D}_{33}=G=\frac{E}{2\left(1+\nu \right)}$ (10)

The Cauchy stress tensor has also been decomposed into symmetric and antisymmetric tensors. In order to obtain the dimensionless Equations (1)-(10), we first write these with hat ( $\stackrel{^}{}$ ) on all quantities and variables indicating that they all have their usual dimensions in terms of length ( $\stackrel{^}{L}$ ), force ( $\stackrel{^}{F}$ ), and time ( $\stackrel{^}{t}$ ). If we choose ${L}_{0}\mathrm{,}{F}_{0}$ and ${t}_{0}$ as reference length, force and time, then the dimensionless length, force and time ( $L\mathrm{,}F$ and $t$ ) are defined as

$L=\frac{\stackrel{^}{L}}{{L}_{0}},\text{\hspace{1em}}F=\frac{\stackrel{^}{F}}{{F}_{0}},\text{\hspace{1em}}t=\frac{\stackrel{^}{t}}{{t}_{0}}$ (11)

If we consider $\stackrel{^}{E}=E{E}_{0}$ , $\stackrel{^}{x}=x{L}_{0}$ , $\stackrel{^}{y}=y{L}_{0}$ , $\stackrel{^}{m}=m{m}_{0}$ , ${m}_{0}=\frac{{\tau}_{0}}{{L}_{0}}$ , ${F}_{0}=\frac{{\tau}_{0}}{{L}_{0}}$ , $\underset{\u02dc}{\stackrel{^}{\alpha}}=\underset{\u02dc}{\alpha}{E}_{0}$ , $\underset{\u02dc}{\stackrel{^}{\beta}}=\underset{\u02dc}{\beta}{E}_{0}$ and choose ${L}_{0}$ , ${E}_{0}$ , then we obtain the dimensionless form of Equations (1)-(10). In these $\frac{{E}_{0}}{{m}_{0}{L}_{0}}$ is in fact unity but has been left in

the constitutive theories for the moment tensor for sake of clarity. Equations (1)-(9) are a system of eleven partial differential equations in eleven dependent variables $u$ , $v$ , ${}_{s}{\sigma}_{xx}$ , ${}_{s}{\sigma}_{yy}$ , ${}_{s}{\sigma}_{xy}$ , ${}_{a}{\sigma}_{yx}$ , ${}_{s}{m}_{xz}$ , ${}_{s}{m}_{yz}$ , ${}_{a}{m}_{xz}$ , ${}_{a}{m}_{yz}$ and ${}_{i}{\Theta}_{z}$ . We substitute ${}_{a}{\sigma}_{yx}$ from (3) into (1) and (2).

$\frac{\partial \left({}_{s}{\sigma}_{xx}\right)}{\partial x}+\frac{\partial \left({}_{s}{\sigma}_{yx}\right)}{\partial y}+\frac{1}{2}\frac{\partial}{\partial y}\left(\frac{\partial \left({m}_{xz}\right)}{\partial x}+\frac{\partial \left({m}_{yz}\right)}{\partial y}\right)={A}_{1}\left(u,v\right)=0,\text{\hspace{1em}}\forall \text{\hspace{0.17em}}x,y\in {\Omega}_{xy}$ (12)

$\frac{\partial \left({}_{s}{\sigma}_{xy}\right)}{\partial x}+\frac{\partial \left({}_{s}{\sigma}_{yy}\right)}{\partial y}-\frac{1}{2}\frac{\partial}{\partial x}\left(\frac{\partial \left({m}_{xz}\right)}{\partial x}+\frac{\partial \left({m}_{yz}\right)}{\partial y}\right)={A}_{2}\left(u,v\right)=0,\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\forall \text{\hspace{0.17em}}x,y\in {\Omega}_{xy}$ (13)

Equations ((12) and (13)) form the basis for finite element formulation based on GM/WF.

3. Finite Element Formulation

The symmetric stress tensor components are defined by (4) and are function of u and v. Likewise ${m}_{xz}$ and ${m}_{yz}$ consists of the decomposition in (3) and the constitutive theories for symmetric and antisymmetric Cauchy moment tensors are given by (5)-(8) with ${}_{i}{\Theta}_{z}$ defined by (9), hence there are also functions of the gradients of u and v. Let ${\stackrel{\xaf}{\Omega}}_{xy}^{T}={{\displaystyle {\cup}_{}}}_{e}\text{\hspace{0.05em}}{\stackrel{\xaf}{\Omega}}_{xy}^{e}$ be the discretization of ${\stackrel{\xaf}{\Omega}}_{xy}$ , domain of definition of mathematical model in which ${\stackrel{\xaf}{\Omega}}_{xy}^{e}={\Omega}_{xy}^{e}\cup {\Gamma}^{e}$ is a typical finite element. ${A}_{1}$ and ${A}_{2}$ in (12) and (13) are differential operators that act on u and v ( ${m}_{yz}$ and ${m}_{xz}$ are functions of the gradients of u and v and so are stresses). The balance of linear momenta equations in the form (12) and (13) are helpful in keeping the derivation of integral form based on GM/WF compact.

Let

$\begin{array}{l}{w}_{1}=\delta {u}_{h}\\ {w}_{2}=\delta {v}_{h}\end{array}$ (14)

in which ${u}_{h}$ and ${v}_{h}$ are approximations of u and v over ${\stackrel{\xaf}{\Omega}}_{xy}^{T}$ . Then, based on the fundamental lemma of the calculus of variations [46] - [60] we construct scalar products of (1) and (2) with test functions ${w}_{1}$ and ${w}_{2}$ over the discretization ${\stackrel{\xaf}{\Omega}}_{xy}^{T}$ and set them to zero.

$\begin{array}{l}{\left({A}_{1}\left({u}_{h},{v}_{h}\right),{w}_{1}\right)}_{{\stackrel{\xaf}{\Omega}}_{xy}^{T}}=0\text{\hspace{0.05em}}\text{\hspace{0.05em}};\text{\hspace{1em}}{w}_{1}=\delta {u}_{h}\\ {\left({A}_{2}\left({u}_{h},{v}_{h}\right),{w}_{2}\right)}_{{\stackrel{\xaf}{\Omega}}_{xy}^{T}}=0\text{\hspace{0.05em}}\text{\hspace{0.05em}};\text{\hspace{1em}}{w}_{2}=\delta {v}_{h}\end{array}$ (15)

or

$\begin{array}{l}\left({A}_{1}\left({u}_{h},{v}_{h}\right),{w}_{1}\right)=\underset{e}{{\displaystyle \sum}}{\left({A}_{1}\left({u}_{h}^{e},{v}_{h}^{e}\right),{w}_{1}\right)}_{{\stackrel{\xaf}{\Omega}}_{xy}^{e}}=0\\ \left({A}_{2}\left({u}_{h},{v}_{h}\right),{w}_{2}\right)=\underset{e}{{\displaystyle \sum}}{\left({A}_{2}\left({u}_{h}^{e},{v}_{h}^{e}\right),{w}_{2}\right)}_{{\stackrel{\xaf}{\Omega}}_{xy}^{e}}=0\end{array}$ (16)

In which ${u}_{h}={\displaystyle {\cup}_{e}{u}_{h}^{e}}$ and ${v}_{h}={\displaystyle {\cup}_{e}{v}_{h}^{e}}$ where ${u}_{h}^{e}$ , ${v}_{h}^{e}$ are local approximations of u and v over an element e with domain ${\stackrel{\xaf}{\Omega}}_{xy}^{e}$ . We consider

${\left({A}_{1}\left({u}_{h}^{e}\mathrm{,}{v}_{h}^{e}\right)\mathrm{,}{w}_{1}\right)}_{{\stackrel{\xaf}{\Omega}}_{xy}^{e}}$ and ${\left({A}_{2}\left({u}_{h}^{e}\mathrm{,}{v}_{h}^{e}\right)\mathrm{,}{w}_{2}\right)}_{{\stackrel{\xaf}{\Omega}}_{xy}^{e}}$ in which ${w}_{1}=\delta {u}_{h}^{e}$ , ${w}_{2}=\delta {v}_{h}^{e}$ and substitute for ${A}_{1}$ and ${A}_{2}$ from (12) and (13)

${\left({A}_{1}\left({u}_{h}^{e}\mathrm{,}{v}_{h}^{e}\right)\mathrm{,}{w}_{1}\right)}_{{\stackrel{\xaf}{\Omega}}_{xy}^{e}}={\left(\frac{\partial \left({}_{s}{\sigma}_{xx}\right)}{\partial x}+\frac{\partial \left({}_{s}{\sigma}_{yx}\right)}{\partial y}-\frac{1}{2}\frac{\partial}{\partial y}\left(\frac{\partial \left({m}_{xz}\right)}{\partial x}+\frac{\partial \left({m}_{yz}\right)}{\partial y}\right),{w}_{1}\right)}_{{\stackrel{\xaf}{\Omega}}_{xy}^{e}}$ (17)

${\left({A}_{2}\left({u}_{h}^{e},{v}_{h}^{e}\right),{w}_{2}\right)}_{{\stackrel{\xaf}{\Omega}}_{xy}^{e}}={\left(\frac{\partial \left({}_{s}{\sigma}_{xy}\right)}{\partial x}+\frac{\partial \left({}_{s}{\sigma}_{yy}\right)}{\partial y}+\frac{1}{2}\frac{\partial}{\partial x}\left(\frac{\partial \left({m}_{xz}\right)}{\partial x}+\frac{\partial \left({m}_{yz}\right)}{\partial y}\right),{w}_{2}\right)}_{{\stackrel{\xaf}{\Omega}}_{xy}^{e}}$ (18)

Integration by parts once for each term in (17) and (18) yields (noting that ${\stackrel{\xaf}{\Omega}}_{xy}^{e}={\stackrel{\xaf}{\Omega}}_{xy}^{e}\cup {\Gamma}^{e}$ , ${\Gamma}^{e}$ being closed boundary of ${\stackrel{\xaf}{\Omega}}_{xy}^{e}$ )

$\begin{array}{l}{\left({A}_{1}\left({u}_{h}^{e},{v}_{h}^{e}\right),{w}_{1}\right)}_{{\stackrel{\xaf}{\Omega}}_{xy}^{e}}\\ ={\displaystyle \underset{{\stackrel{\xaf}{\Omega}}_{xy}^{e}}{\int}}\left(-{}_{s}\sigma {}_{xx}\left(\frac{\partial {w}_{1}}{\partial x}\right)-{}_{s}\sigma {}_{yx}\left(\frac{\partial {w}_{1}}{\partial y}\right)+\frac{1}{2}\left(\frac{\partial \left({m}_{xz}\right)}{\partial x}+\frac{\partial \left({m}_{yz}\right)}{\partial y}\right)\frac{\partial {w}_{1}}{\partial y}\right)\text{d}\Omega \\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}+{\displaystyle \underset{{\Gamma}^{e}}{\oint}}\left({}_{s}{\sigma}_{xx}\left({n}_{x}\right)+{}_{s}\sigma {}_{yx}\left({n}_{y}\right)\right){w}_{1}\text{d}\Gamma -\frac{1}{2}{\displaystyle \underset{{\Gamma}^{e}}{\oint}}\left(\frac{\partial \left({m}_{xz}\right)}{\partial x}+\frac{\partial \left({m}_{yz}\right)}{\partial y}\right){n}_{y}{w}_{1}\text{d}\Gamma \end{array}$ (19)

$\begin{array}{l}{\left({A}_{2}\left({u}_{h}^{e},{v}_{h}^{e}\right),{w}_{2}\right)}_{{\stackrel{\xaf}{\Omega}}_{xy}^{e}}\\ ={\displaystyle \underset{{\stackrel{\xaf}{\Omega}}_{xy}^{e}}{\int}}\left(-{}_{s}\sigma {}_{xy}\left(\frac{\partial {w}_{2}}{\partial x}\right)-{}_{s}\sigma {}_{yy}\left(\frac{\partial {w}_{2}}{\partial y}\right)-\frac{1}{2}\left(\frac{\partial \left({m}_{xz}\right)}{\partial x}+\frac{\partial \left({m}_{yz}\right)}{\partial y}\right)\frac{\partial {w}_{2}}{\partial x}\right)\text{d}\Omega \\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}+{\displaystyle \underset{{\Gamma}^{e}}{\oint}}\left({}_{s}{\sigma}_{xy}\left({n}_{x}\right)+{}_{s}\sigma {}_{yy}\left({n}_{y}\right)\right){w}_{2}\text{d}\Gamma +\frac{1}{2}{\displaystyle \underset{{\Gamma}^{e}}{\oint}}\left(\frac{\partial \left({m}_{xz}\right)}{\partial x}+\frac{\partial \left({m}_{yz}\right)}{\partial y}\right){n}_{x}{w}_{2}\text{d}\Gamma \end{array}$ (20)

Integrating by parts again for the moment terms in (19) and (20)

$\begin{array}{l}{\left({A}_{1}\left({u}_{h}^{e},{v}_{h}^{e}\right),{w}_{1}\right)}_{{\stackrel{\xaf}{\Omega}}_{xy}^{e}}\\ ={\displaystyle \underset{{\stackrel{\xaf}{\Omega}}_{xy}^{e}}{\int}}\left(-{}_{s}\sigma {}_{xx}\left(\frac{\partial {w}_{1}}{\partial x}\right)-{}_{s}\sigma {}_{yx}\left(\frac{\partial {w}_{1}}{\partial y}\right)-\frac{1}{2}\left({m}_{xz}\frac{{\partial}^{2}\left({w}_{1}\right)}{\partial x\partial y}+{m}_{yz}\frac{{\partial}^{2}\left({w}_{1}\right)}{\partial {y}^{2}}\right)\right)\text{d}\Omega \\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}+{\displaystyle \underset{{\Gamma}^{e}}{\oint}}{w}_{1}\left({}_{s}{\sigma}_{xx}\left({n}_{x}\right)+{}_{s}\sigma {}_{yx}\left({n}_{y}\right)-\frac{1}{2}\left(\frac{\partial \left({m}_{xz}\right)}{\partial x}+\frac{\partial \left({m}_{yz}\right)}{\partial y}\right){n}_{y}\right)\text{d}\Gamma \\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}+\frac{1}{2}{\displaystyle \underset{{\Gamma}_{e}}{\oint}}\frac{\partial {w}_{1}}{\partial y}\left({m}_{xz}\left({n}_{x}\right)+{m}_{yz}\left({n}_{y}\right)\right)\text{d}\Gamma \end{array}$ (21)

$\begin{array}{l}{\left({A}_{2}\left({u}_{h}^{e},{v}_{h}^{e}\right),{w}_{2}\right)}_{{\stackrel{\xaf}{\Omega}}_{xy}^{e}}\\ ={\displaystyle \underset{{\stackrel{\xaf}{\Omega}}_{xy}^{e}}{\int}}\left(-{}_{s}\sigma {}_{xy}\left(\frac{\partial {w}_{2}}{\partial x}\right)-{}_{s}\sigma {}_{yy}\left(\frac{\partial {w}_{2}}{\partial y}\right)+\frac{1}{2}\left({m}_{xz}\frac{{\partial}^{2}\left({w}_{2}\right)}{\partial {x}^{2}}+{m}_{yz}\frac{{\partial}^{2}\left({w}_{2}\right)}{\partial y\partial x}\right)\right)\text{d}\Omega \\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}+{\displaystyle \underset{{\Gamma}^{e}}{\oint}}{w}_{2}\left({}_{s}{\sigma}_{xy}\left({n}_{x}\right)+{}_{s}\sigma {}_{yy}\left({n}_{y}\right)+\frac{1}{2}\left(\frac{\partial \left({m}_{xz}\right)}{\partial x}+\frac{\partial \left({m}_{yz}\right)}{\partial y}\right){n}_{x}\right)\text{d}\Gamma \\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}-\frac{1}{2}{\displaystyle \underset{{\Gamma}_{e}}{\oint}}\frac{\partial {w}_{2}}{\partial x}\left({m}_{xz}\left({n}_{x}\right)+{m}_{yz}\left({n}_{y}\right)\right)\text{d}\Gamma \end{array}$ (22)

Let

$\begin{array}{l}{t}_{x}={}_{s}\sigma {}_{xx}\left({n}_{x}\right)+{}_{s}\sigma {}_{yx}\left({n}_{y}\right)-\frac{1}{2}\left(\frac{\partial}{\partial x}\left({m}_{xz}\right)+\frac{\partial}{\partial y}\left({m}_{yz}\right)\right){n}_{y}\\ {t}_{y}={}_{s}\sigma {}_{xy}\left({n}_{x}\right)+{}_{s}\sigma {}_{yy}\left({n}_{y}\right)+\frac{1}{2}\left(\frac{\partial}{\partial x}\left({m}_{xz}\right)+\frac{\partial}{\partial y}\left({m}_{yz}\right)\right){n}_{x}\\ \text{and}\text{\hspace{1em}}{m}_{n}={m}_{xz}\left({n}_{x}\right)+{m}_{yz}\left({n}_{y}\right)\end{array}$ (23)

Using (23) we can write (21) and (22) as follows

$\begin{array}{l}{\left({A}_{1}\left({u}_{h}^{e},{v}_{h}^{e}\right),{w}_{1}\right)}_{{\stackrel{\xaf}{\Omega}}_{xy}^{e}}\\ ={\displaystyle \underset{{\stackrel{\xaf}{\Omega}}_{xy}^{e}}{\int}}\left(-{}_{s}\sigma {}_{xx}\left(\frac{\partial {w}_{1}}{\partial x}\right)-{}_{s}\sigma {}_{yx}\left(\frac{\partial {w}_{1}}{\partial y}\right)-\frac{1}{2}\left({m}_{xz}\frac{{\partial}^{2}\left({w}_{1}\right)}{\partial x\partial y}+{m}_{yz}\frac{{\partial}^{2}\left({w}_{1}\right)}{\partial {y}^{2}}\right)\right)\text{d}\Omega \\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}+{\displaystyle \underset{{\Gamma}^{e}}{\oint}}\left({w}_{1}{t}_{x}\right)\text{d}\Gamma +{\displaystyle \underset{{\Gamma}^{e}}{\oint}}\frac{\partial {w}_{1}}{\partial y}\left(\frac{{m}_{n}}{2}\right)\text{d}\Gamma \end{array}$ (24)

$\begin{array}{l}{\left({A}_{2}\left({u}_{h}^{e},{v}_{h}^{e}\right),{w}_{2}\right)}_{{\stackrel{\xaf}{\Omega}}_{xy}^{e}}\\ ={\displaystyle \underset{{\stackrel{\xaf}{\Omega}}_{xy}^{e}}{\int}}\left(-{}_{s}\sigma {}_{xy}\left(\frac{\partial {w}_{2}}{\partial x}\right)-{}_{s}\sigma {}_{yy}\left(\frac{\partial {w}_{2}}{\partial y}\right)+\frac{1}{2}\left({m}_{xz}\frac{{\partial}^{2}\left({w}_{2}\right)}{\partial {x}^{2}}+{m}_{yz}\frac{{\partial}^{2}\left({w}_{2}\right)}{\partial y\partial x}\right)\right)\text{d}\Omega \\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}+{\displaystyle \underset{{\Gamma}^{e}}{\oint}}\left({w}_{2}{t}_{y}\right)\text{d}\Gamma +{\displaystyle \underset{{\Gamma}^{e}}{\oint}}\frac{\partial {w}_{2}}{\partial x}\left(-\frac{{m}_{n}}{2}\right)\text{d}\Gamma \end{array}$ (25)

From (24) and (25) we can conclude that the primary and the secondary variables (PV and SV) are

Substituting for ${}_{s}{\sigma}_{xx}$ , ${}_{s}{\sigma}_{xy}$ , ${}_{s}{\sigma}_{yy}$ from (4) into (24) and (25) and ${}_{i}{\Theta}_{z}$ from (9) into (5)-(8) and then (5)-(8) into (24) and (25) we can write

${\left({A}_{1}\left({u}_{h}^{e},{v}_{h}^{e}\right),{w}_{1}\right)}_{{\stackrel{\xaf}{\Omega}}_{xy}^{e}}={B}_{1}^{e}\left({u}_{h}^{e},{v}_{h}^{e};{w}_{1}\right)-{l}_{1}^{e}\left({w}_{1}\right)$ (26)

${\left({A}_{2}\left({u}_{h}^{e},{v}_{h}^{e}\right),{w}_{2}\right)}_{{\stackrel{\xaf}{\Omega}}_{xy}^{e}}={B}_{2}^{e}\left({u}_{h}^{e},{v}_{h}^{e};{w}_{2}\right)-{l}_{2}^{e}\left({w}_{2}\right)$ (27)

in which

$\begin{array}{c}{B}_{1}^{e}\left({u}_{h}^{e},{v}_{h}^{e};{w}_{1}\right)={\displaystyle \underset{{\stackrel{\xaf}{\Omega}}_{xy}^{e}}{\int}}\{-\left({D}_{11}\frac{\partial {u}_{h}^{e}}{\partial x}+{D}_{12}\frac{\partial {v}_{h}^{e}}{\partial y}\right)\frac{\partial {w}_{1}}{\partial x}-{D}_{33}\left(\frac{\partial {u}_{h}^{e}}{\partial y}+\frac{\partial {v}_{h}^{e}}{\partial x}\right)\frac{\partial {w}_{1}}{\partial y}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}-\frac{1}{2}\left(\underset{\u02dc}{\alpha}-\underset{\u02dc}{\beta}\right)\left[\left(\frac{{\partial}^{2}{u}_{h}^{e}}{\partial x\partial y}-\frac{{\partial}^{2}{v}_{h}^{e}}{\partial {x}^{2}}\right)\frac{{\partial}^{2}{w}_{1}}{\partial x\partial y}+\left(\frac{{\partial}^{2}{u}_{h}^{e}}{\partial {y}^{2}}-\frac{{\partial}^{2}{v}_{h}^{e}}{\partial y\partial x}\right)\frac{{\partial}^{2}{w}_{1}}{\partial {y}^{2}}\right]\}\text{d}\Omega \end{array}$ (28)

${l}_{1}^{e}\left({w}_{1}\right)=-{\displaystyle \underset{{\Gamma}^{e}}{\oint}}{w}_{1}{t}_{x}\text{d}\Gamma -{\displaystyle \underset{{\Gamma}^{e}}{\oint}}\frac{\partial {w}_{1}}{\partial y}\left(\frac{{m}_{n}}{2}\right)\text{d}\Gamma $ (29)

$\begin{array}{c}{B}_{2}^{e}\left({u}_{h}^{e},{v}_{h}^{e};{w}_{2}\right)\\ ={\displaystyle \underset{{\stackrel{\xaf}{\Omega}}_{xy}^{e}}{\int}}\{-{D}_{33}\left(\frac{\partial {u}_{h}^{e}}{\partial y}+\frac{\partial {v}_{h}^{e}}{\partial x}\right)\frac{\partial {w}_{2}}{\partial x}-\left({D}_{21}\frac{\partial {u}_{h}^{e}}{\partial x}+{D}_{22}\frac{\partial {v}_{h}^{e}}{\partial y}\right)\frac{\partial {w}_{2}}{\partial y}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+\frac{1}{2}\left(\underset{\u02dc}{\alpha}-\underset{\u02dc}{\beta}\right)\left[\left(\frac{{\partial}^{2}{u}_{h}^{e}}{\partial x\partial y}-\frac{{\partial}^{2}{v}_{h}^{e}}{\partial {x}^{2}}\right)\frac{{\partial}^{2}{w}_{2}}{\partial {x}^{2}}+\left(\frac{{\partial}^{2}{u}_{h}^{e}}{\partial {y}^{2}}-\frac{{\partial}^{2}{v}_{h}^{e}}{\partial x\partial y}\right)\frac{{\partial}^{2}{w}_{2}}{\partial y\partial x}\right]\}\text{d}\Omega \end{array}$ (30)

${l}_{2}^{e}\left({w}_{2}\right)=-{\displaystyle \underset{{\Gamma}^{e}}{\oint}}{w}_{2}{t}_{y}\text{d}\Gamma -{\displaystyle \underset{{\Gamma}^{e}}{\oint}}\frac{\partial {w}_{2}}{\partial x}\left(-\frac{{m}_{n}}{2}\right)\text{d}\Gamma $ (31)

Functionals ${B}_{1}^{e}\left(\cdot \mathrm{,}\cdot \right)$ , ${B}_{2}^{e}\left(\cdot \mathrm{,}\cdot \right)$ , ${l}_{1}^{e}(\cdot )$ , ${l}_{2}^{e}(\cdot )$ are linear in all of their arguments. We note that ${l}_{1}^{e}(\cdot )$ , ${l}_{2}^{e}(\cdot )$ are concomitants resulting only because of integration by parts. We can represent (28)-(31) in matrix and vector form

$\left\{\begin{array}{c}{\left({A}_{1}\left({u}_{h}^{e},{v}_{h}^{e}\right),{w}_{1}\right)}_{{\stackrel{\xaf}{\Omega}}_{xy}^{e}}\\ {\left({A}_{2}\left({u}_{h}^{e},{v}_{h}^{e}\right),{w}_{2}\right)}_{{\stackrel{\xaf}{\Omega}}_{xy}^{e}}\end{array}\right\}=\left\{\begin{array}{c}{B}_{1}^{e}\left({u}_{h}^{e},{v}_{h}^{e};{w}_{1}\right)\\ {B}_{2}^{e}\left({u}_{h}^{e},{v}_{h}^{e};{w}_{2}\right)\end{array}\right\}-\left\{\begin{array}{c}{l}_{1}^{e}\left({w}_{1}\right)\\ {l}_{2}^{e}\left({w}_{2}\right)\end{array}\right\}$ (32)

$\mathrm{=}\left\{\begin{array}{c}{b}_{11}^{e}\left({u}_{h}^{e}\mathrm{,}{w}_{1}\right)+{b}_{12}^{e}\left({v}_{h}^{e}\mathrm{,}{w}_{1}\right)\\ {b}_{21}^{e}\left({u}_{h}^{e}\mathrm{,}{w}_{2}\right)+{b}_{22}^{e}\left({v}_{h}^{e}\mathrm{,}{w}_{2}\right)\end{array}\right\}-\left\{\begin{array}{c}{l}_{1}^{e}\left({w}_{1}\right)\\ {l}_{2}^{e}\left({w}_{2}\right)\end{array}\right\}$ (33)

in which

$\begin{array}{c}{b}_{11}^{e}\left({u}_{h}^{e},{w}_{1}\right)={\displaystyle \underset{{\stackrel{\xaf}{\Omega}}_{xy}^{e}}{\int}}\{-\left({D}_{11}\frac{\partial {u}_{h}^{e}}{\partial x}\frac{\partial {w}_{1}}{\partial x}\right)-{D}_{33}\left(\frac{\partial {u}_{h}^{e}}{\partial y}\frac{\partial {w}_{1}}{\partial y}\right)\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}-\frac{1}{2}\left(\underset{\u02dc}{\alpha}-\underset{\u02dc}{\beta}\right)\left[\left(\frac{{\partial}^{2}{u}_{h}^{e}}{\partial x\partial y}\frac{{\partial}^{2}{w}_{1}}{\partial x\partial y}\right)+\left(\frac{{\partial}^{2}{u}_{h}^{e}}{\partial {y}^{2}}\frac{{\partial}^{2}{w}_{1}}{\partial {y}^{2}}\right)\right]\}\text{d}\Omega \end{array}$ (34)

$\begin{array}{c}{b}_{12}^{e}\left({v}_{h}^{e},{w}_{1}\right)={\displaystyle \underset{{\stackrel{\xaf}{\Omega}}_{xy}^{e}}{\int}}\{-\left({D}_{12}\frac{\partial {v}_{h}^{e}}{\partial y}\frac{\partial {w}_{1}}{\partial x}\right)-{D}_{33}\left(\frac{\partial {v}_{h}^{e}}{\partial x}\frac{\partial {w}_{1}}{\partial y}\right)\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}-\frac{1}{2}\left(\underset{\u02dc}{\alpha}-\underset{\u02dc}{\beta}\right)\left[-\left(\frac{{\partial}^{2}{v}_{h}^{e}}{\partial {x}^{2}}\frac{{\partial}^{2}{w}_{1}}{\partial x\partial y}\right)-\left(\frac{{\partial}^{2}{v}_{h}^{e}}{\partial y\partial x}\frac{{\partial}^{2}{w}_{1}}{\partial {y}^{2}}\right)\right]\}\text{d}\Omega \end{array}$ (35)

$\begin{array}{c}{b}_{21}^{e}\left({u}_{h}^{e},{w}_{2}\right)={\displaystyle \underset{{\stackrel{\xaf}{\Omega}}_{xy}^{e}}{\int}}\{-\left({D}_{33}\frac{\partial {u}_{h}^{e}}{\partial y}\frac{\partial {w}_{2}}{\partial x}\right)-{D}_{21}\left(\frac{\partial {u}_{h}^{e}}{\partial x}\frac{\partial {w}_{2}}{\partial y}\right)\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+\frac{1}{2}\left(\underset{\u02dc}{\alpha}-\underset{\u02dc}{\beta}\right)\left[\left(\frac{{\partial}^{2}{u}_{h}^{e}}{\partial x\partial y}\frac{{\partial}^{2}{w}_{2}}{\partial {x}^{2}}\right)+\left(\frac{{\partial}^{2}{u}_{h}^{e}}{\partial {y}^{2}}\frac{{\partial}^{2}{w}_{2}}{\partial y\partial x}\right)\right]\}\text{d}\Omega \end{array}$ (36)

$\begin{array}{c}{b}_{22}^{e}\left({v}_{h}^{e},{w}_{2}\right)={\displaystyle \underset{{\stackrel{\xaf}{\Omega}}_{xy}^{e}}{\int}}\{-\left({D}_{33}\frac{\partial {v}_{h}^{e}}{\partial x}\frac{\partial {w}_{2}}{\partial x}\right)-{D}_{22}\left(\frac{\partial {v}_{h}^{e}}{\partial y}\frac{\partial {w}_{2}}{\partial y}\right)\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+\frac{1}{2}\left(\underset{\u02dc}{\alpha}-\underset{\u02dc}{\beta}\right)\left[-\left(\frac{{\partial}^{2}{v}_{h}^{e}}{\partial {x}^{2}}\frac{{\partial}^{2}{w}_{2}}{\partial {x}^{2}}\right)-\left(\frac{{\partial}^{2}{v}_{h}^{e}}{\partial x\partial y}\frac{{\partial}^{2}{w}_{2}}{\partial y\partial x}\right)\right]\}\text{d}\Omega \end{array}$ (37)

The functions ${b}_{ij}^{e};i,j=1,2$ have the following properties.

$\begin{array}{l}{b}_{11}^{e}\left({u}_{h}^{e},{w}_{1}\right)={b}_{11}^{e}\left({w}_{1},{u}_{h}^{e}\right)\Rightarrow {b}_{11}^{e}\left(\cdot ,\cdot \right)\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{is}\text{\hspace{0.17em}}\text{symmetric}\\ {b}_{22}^{e}\left({v}_{h}^{e},{w}_{2}\right)={b}_{22}^{e}\left({w}_{2},{v}_{h}^{e}\right)\Rightarrow {b}_{22}^{e}\left(\cdot ,\cdot \right)\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{is}\text{\hspace{0.17em}}\text{symmetric}\\ {b}_{12}^{e}\left({w}_{2},{u}_{h}^{e}\right)={b}_{21}^{e}\left({u}_{h}^{e},{w}_{2}\right)\\ {b}_{21}^{e}\left({w}_{1},{v}_{h}^{e}\right)={b}_{12}^{e}\left({v}_{h}^{e},{w}_{1}\right)\end{array}\}$ (38)

Equations (32) are the weak form of the mathematical model (12) and (13). Equations (38) imply that the element equations constructed from (38) using local approximation ${u}_{h}^{e}$ and ${v}_{h}^{e}$ will contain symmetric element coefficient matrix.

Let ${\stackrel{\xaf}{\Omega}}_{xy}^{e}$ be a nine node p-version hierarchical element with local approximation in higher order scalar product space ${H}^{k\mathrm{,}p}\left({\stackrel{\xaf}{\Omega}}_{xy}^{e}\right)$ [57] [58] [59] [60] . Consider ${\stackrel{\xaf}{\Omega}}_{xy}^{e}\to {\stackrel{\xaf}{\Omega}}_{\xi \eta}^{e}=\left[-\mathrm{1,1}\right]\times \left[-\mathrm{1,1}\right]$ , a map of ${\stackrel{\xaf}{\Omega}}_{xy}^{e}$ in $\xi \mathrm{,}\eta $ space in a two unit square. Then, we can write

${u}_{h}^{e}\left(\xi ,\eta \right)={\displaystyle \underset{i=1}{\overset{{n}^{u}}{\sum}}}\text{\hspace{0.05em}}\text{\hspace{0.05em}}{N}_{i}^{u}\left(\xi ,\eta \right)\left({}^{u}{\delta}_{i}^{e}\right)$ (39)

${v}_{h}^{e}\left(\xi ,\eta \right)={\displaystyle \underset{i=1}{\overset{{n}^{v}}{\sum}}}\text{\hspace{0.05em}}\text{\hspace{0.05em}}{N}_{i}^{v}\left(\xi ,\eta \right)\left({}^{v}{\delta}_{i}^{e}\right)$ (40)

${N}_{i}^{u}\left(\xi \mathrm{,}\eta \right)$ and ${N}_{i}^{v}\left(\xi \mathrm{,}\eta \right)$ are local approximation functions and $\left({}^{u}{\delta}_{i}^{e}\right)$ and $\left({}^{v}{\delta}_{i}^{e}\right)$ are corresponding nodal degrees of freedom for u and v. Using (39) and (40)

$\begin{array}{l}{w}_{1}=\delta {u}_{h}^{e}={N}_{j}^{u}\left(\xi ,\eta \right);\text{\hspace{1em}}j=1,2,\cdots ,{n}^{u}\\ {w}_{2}=\delta {v}_{h}^{e}={N}_{k}^{v}\left(\xi ,\eta \right);\text{\hspace{1em}}k=1,2,\cdots ,{n}^{v}\end{array}$ (41)

Let the total degrees of freedom for an element e be $\left\{{\delta}^{e}\right\}$

$\left\{{\delta}^{e}\right\}=\left\{{}^{u}{\delta}^{e}\right\}\cup \left\{{}^{v}{\delta}^{e}\right\}=\left\{\begin{array}{c}\left\{{}^{u}{\delta}^{e}\right\}\\ \left\{{}^{v}{\delta}^{e}\right\}\end{array}\right\}$ (42)

Substituting from (39)-(41) into (34)-(37), we can write

$\begin{array}{l}{b}_{11}^{e}={\displaystyle \underset{{\stackrel{\xaf}{\Omega}}_{xy}^{e}}{\int}}\{-{D}_{11}\left({\displaystyle \underset{i=1}{\overset{{n}^{u}}{\sum}}}\frac{\partial {N}_{i}^{u}}{\partial x}\left({}^{u}{\delta}_{i}^{e}\right)\right)\frac{\partial {N}_{j}^{u}}{\partial x}-{D}_{33}\left({\displaystyle \underset{i=1}{\overset{{n}^{u}}{\sum}}}\frac{\partial {N}_{i}^{u}}{\partial y}\left({}^{u}{\delta}_{i}^{e}\right)\right)\frac{\partial {N}_{j}^{u}}{\partial y}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.05em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.05em}}\text{\hspace{0.17em}}-\frac{1}{2}\left(\underset{\u02dc}{\alpha}-\underset{\u02dc}{\beta}\right)[\left({\displaystyle \underset{i=1}{\overset{{n}^{u}}{\sum}}}\frac{{\partial}^{2}{N}_{i}^{u}}{\partial x\partial y}\left({}^{u}{\delta}_{i}^{e}\right)\right)\frac{{\partial}^{2}{N}_{j}^{u}}{\partial x\partial y}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.05em}}\text{\hspace{0.05em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}+\left({\displaystyle \underset{i=1}{\overset{{n}^{u}}{\sum}}}\frac{{\partial}^{2}{N}_{i}^{u}}{\partial {y}^{2}}\left({}^{u}{\delta}_{i}^{e}\right)\right)\frac{{\partial}^{2}{N}_{j}^{u}}{\partial {y}^{2}}]\}\text{d}\Omega \text{\hspace{0.05em}};\text{\hspace{0.17em}}\text{\hspace{0.17em}}j=1,2,\cdots ,{n}^{u}\end{array}$ (43)

$\begin{array}{l}{b}_{12}^{e}={\displaystyle \underset{{\stackrel{\xaf}{\Omega}}_{xy}^{e}}{\int}}\{-{D}_{12}\left({\displaystyle \underset{i=1}{\overset{{n}^{v}}{\sum}}}\frac{\partial {N}_{i}^{v}}{\partial y}\left({}^{v}{\delta}_{i}^{e}\right)\right)\frac{\partial {N}_{j}^{u}}{\partial x}-{D}_{33}\left({\displaystyle \underset{i=1}{\overset{{n}^{v}}{\sum}}}\frac{\partial {N}_{i}^{v}}{\partial y}\left({}^{v}{\delta}_{i}^{e}\right)\right)\frac{\partial {N}_{j}^{u}}{\partial y}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}+\frac{1}{2}\left(\underset{\u02dc}{\alpha}-\underset{\u02dc}{\beta}\right)[\left({\displaystyle \underset{i=1}{\overset{{n}^{v}}{\sum}}}\frac{{\partial}^{2}{N}_{i}^{v}}{\partial {x}^{2}}\left({}^{v}{\delta}_{i}^{e}\right)\right)\frac{{\partial}^{2}{N}_{j}^{u}}{\partial x\partial y}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}+\left({\displaystyle \underset{i=1}{\overset{{n}^{v}}{\sum}}}\frac{{\partial}^{2}{N}_{i}^{v}}{\partial y\partial x}\left({}^{v}{\delta}_{i}^{e}\right)\right)\frac{{\partial}^{2}{N}_{j}^{u}}{\partial {y}^{2}}]\}\text{d}\Omega \text{\hspace{0.17em}};\text{\hspace{1em}}j=1,2,\cdots ,{n}^{u}\end{array}$ (44)

$\begin{array}{l}{b}_{21}^{e}={\displaystyle \underset{{\stackrel{\xaf}{\Omega}}_{xy}^{e}}{\int}}\{-{D}_{33}\left({\displaystyle \underset{i=1}{\overset{{n}^{u}}{\sum}}}\frac{\partial {N}_{i}^{u}}{\partial y}\left({}^{u}{\delta}_{i}^{e}\right)\right)\frac{\partial {N}_{j}^{v}}{\partial x}-{D}_{21}\left({\displaystyle \underset{i=1}{\overset{{n}^{u}}{\sum}}}\frac{\partial {N}_{i}^{u}}{\partial x}\left({}^{u}{\delta}_{i}^{e}\right)\right)\frac{\partial {N}_{j}^{v}}{\partial y}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}+\frac{1}{2}\left(\underset{\u02dc}{\alpha}-\underset{\u02dc}{\beta}\right)[\left({\displaystyle \underset{i=1}{\overset{{n}^{u}}{\sum}}}\frac{{\partial}^{2}{N}_{i}^{u}}{\partial x\partial y}\left({}^{u}{\delta}_{i}^{e}\right)\right)\frac{{\partial}^{2}{N}_{j}^{v}}{\partial {x}^{2}}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}+\left({\displaystyle \underset{i=1}{\overset{{n}^{u}}{\sum}}}\frac{{\partial}^{2}{N}_{i}^{u}}{\partial {y}^{2}}\left({}^{u}{\delta}_{i}^{e}\right)\right)\frac{{\partial}^{2}{N}_{j}^{v}}{\partial y\partial x}]\}\text{d}\Omega \text{\hspace{0.17em}};\text{\hspace{0.17em}}\text{\hspace{1em}}j=1,2,\cdots ,{n}^{v}\end{array}$ (45)

$\begin{array}{l}{b}_{22}^{e}={\displaystyle \underset{{\stackrel{\xaf}{\Omega}}_{xy}^{e}}{\int}}\{-{D}_{33}\left({\displaystyle \underset{i=1}{\overset{{n}^{v}}{\sum}}}\frac{\partial {N}_{i}^{v}}{\partial x}\left({}^{v}{\delta}_{i}^{e}\right)\right)\frac{\partial {N}_{j}^{v}}{\partial x}-{D}_{22}\left({\displaystyle \underset{i=1}{\overset{{n}^{v}}{\sum}}}\frac{\partial {N}_{i}^{v}}{\partial y}\left({}^{v}{\delta}_{i}^{e}\right)\right)\frac{\partial {N}_{j}^{v}}{\partial y}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}-\frac{1}{2}\left(\underset{\u02dc}{\alpha}-\underset{\u02dc}{\beta}\right)[\left({\displaystyle \underset{i=1}{\overset{{n}^{v}}{\sum}}}\frac{{\partial}^{2}{N}_{i}^{v}}{\partial {x}^{2}}\left({}^{v}{\delta}_{i}^{e}\right)\right)\frac{{\partial}^{2}{N}_{j}^{v}}{\partial {x}^{2}}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}+\left({\displaystyle \underset{i=1}{\overset{{n}^{v}}{\sum}}}\frac{{\partial}^{2}{N}_{i}^{v}}{\partial x\partial y}\left({}^{v}{\delta}_{i}^{e}\right)\right)\frac{{\partial}^{2}{N}_{j}^{v}}{\partial x\partial y}]\}\text{d}\Omega \text{\hspace{0.17em}};\text{\hspace{0.17em}}\text{\hspace{1em}}j=1,2,\cdots ,{n}^{v}\end{array}$ (46)

Using (43)-(46) we can write (33) as follows

$\begin{array}{c}\left\{\begin{array}{c}{\left({A}_{1}\left({u}_{h}^{e},{v}_{h}^{e}\right),{w}_{1}\right)}_{{\stackrel{\xaf}{\Omega}}_{xy}^{e}}\\ {\left({A}_{2}\left({u}_{h}^{e},{v}_{h}^{e}\right),{w}_{2}\right)}_{{\stackrel{\xaf}{\Omega}}_{xy}^{e}}\end{array}\right\}=-\left[{K}^{e}\right]\left\{{\delta}^{e}\right\}+\left\{{P}^{e}\right\}\\ =-\left[\begin{array}{cc}\left[{}^{11}{K}^{e}\right]& \left[{}^{12}{K}^{e}\right]\\ \left[{}^{21}{K}^{e}\right]& \left[{}^{22}{K}^{e}\right]\end{array}\right]\left\{\begin{array}{c}\left\{{}^{u}{\delta}^{e}\right\}\\ \left\{{}^{v}{\delta}^{e}\right\}\end{array}\right\}+\left\{\begin{array}{c}\left\{{P}_{1}^{e}\right\}\\ \left\{{P}_{2}^{e}\right\}\end{array}\right\}\end{array}$ (47)

in which

$\begin{array}{l}{}^{11}{K}_{ij}^{e}={\displaystyle \underset{{\stackrel{\xaf}{\Omega}}_{xy}^{e}}{\int}}\{{D}_{11}\frac{\partial {N}_{i}^{u}}{\partial x}\frac{\partial {N}_{j}^{u}}{\partial x}+{D}_{33}\frac{\partial {N}_{i}^{u}}{\partial y}\frac{\partial {N}_{j}^{u}}{\partial y}+\frac{1}{2}\left(\underset{\u02dc}{\alpha}-\underset{\u02dc}{\beta}\right)[\frac{{\partial}^{2}{N}_{i}^{u}}{\partial x\partial y}\frac{{\partial}^{2}{N}_{j}^{u}}{\partial x\partial y}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}+\frac{{\partial}^{2}{N}_{i}^{u}}{\partial {y}^{2}}\frac{{\partial}^{2}{N}_{j}^{u}}{\partial {y}^{2}}]\}\text{d}\Omega \text{\hspace{0.17em}};\text{\hspace{1em}}i,j=1,2,\cdots ,{n}^{u}\end{array}$ (48)

$\begin{array}{l}{}^{12}{K}_{ij}^{e}={\displaystyle \underset{{\stackrel{\xaf}{\Omega}}_{xy}^{e}}{\int}}\{{D}_{12}\frac{\partial {N}_{i}^{u}}{\partial x}\frac{\partial {N}_{j}^{v}}{\partial y}+{D}_{33}\frac{\partial {N}_{i}^{u}}{\partial y}\frac{\partial {N}_{j}^{v}}{\partial x}-\frac{1}{2}\left(\underset{\u02dc}{\alpha}-\underset{\u02dc}{\beta}\right)[\frac{{\partial}^{2}{N}_{i}^{u}}{\partial x\partial y}\frac{{\partial}^{2}{N}_{j}^{v}}{\partial {x}^{2}}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}+\frac{{\partial}^{2}{N}_{i}^{u}}{\partial {y}^{2}}\frac{{\partial}^{2}{N}_{j}^{v}}{\partial y\partial x}]\}\text{d}\Omega \text{\hspace{0.17em}};\text{\hspace{1em}}i=1,2,\cdots ,{n}^{u};\text{\hspace{0.17em}}j=1,2,\cdots ,{n}^{v}\end{array}$ (49)

$\begin{array}{l}{}^{21}{K}_{ij}^{e}={\displaystyle \underset{{\stackrel{\xaf}{\Omega}}_{xy}^{e}}{\int}}\{{D}_{33}\frac{\partial {N}_{i}^{v}}{\partial x}\frac{\partial {N}_{j}^{u}}{\partial y}+{D}_{21}\frac{\partial {N}_{i}^{v}}{\partial y}\frac{\partial {N}_{j}^{u}}{\partial x}-\frac{1}{2}\left(\underset{\u02dc}{\alpha}-\underset{\u02dc}{\beta}\right)[\frac{{\partial}^{2}{N}_{i}^{v}}{\partial {x}^{2}}\frac{{\partial}^{2}{N}_{j}^{u}}{\partial x\partial y}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}+\frac{{\partial}^{2}{N}_{i}^{v}}{\partial y\partial x}\frac{{\partial}^{2}{N}_{j}^{u}}{\partial {y}^{2}}]\}\text{d}\Omega \text{\hspace{0.17em}};\text{\hspace{1em}}i=1,2,\cdots ,{n}^{u};\text{\hspace{0.17em}}j=1,2,\cdots ,{n}^{v}\end{array}$ (50)

$\begin{array}{l}{}^{22}{K}_{ij}^{e}={\displaystyle \underset{{\stackrel{\xaf}{\Omega}}_{xy}^{e}}{\int}}\{{D}_{33}\frac{\partial {N}_{i}^{v}}{\partial x}\frac{\partial {N}_{j}^{v}}{\partial x}+{D}_{22}\frac{\partial {N}_{i}^{v}}{\partial y}\frac{\partial {N}_{j}^{v}}{\partial y}+\frac{1}{2}\left(\underset{\u02dc}{\alpha}-\underset{\u02dc}{\beta}\right)[\frac{{\partial}^{2}{N}_{i}^{v}}{\partial {x}^{2}}\frac{{\partial}^{2}{N}_{j}^{v}}{\partial {x}^{2}}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}+\frac{{\partial}^{2}{N}_{i}^{v}}{\partial x\partial y}\frac{{\partial}^{2}{N}_{j}^{v}}{\partial x\partial y}]\}\text{d}\Omega \text{\hspace{0.17em}};\text{\hspace{0.17em}}\text{\hspace{1em}}i,j=1,2,\cdots ,{n}^{v}\end{array}$ (51)

We note that

$\begin{array}{l}{}^{11}{K}_{ij}^{e}={}^{11}K{}_{ji}^{e}\text{\hspace{0.17em}};\text{\hspace{1em}}i,j=1,2,\cdots ,{n}^{u},\text{\hspace{0.17em}}\text{hence}\text{\hspace{0.17em}}\left[{}^{11}{K}^{e}\right]\text{\hspace{0.17em}}\text{is}\text{\hspace{0.17em}}\text{symmetric}\\ {}^{12}{K}_{ij}^{e}={}^{21}K{}_{ji}^{e}\text{\hspace{0.17em}};\text{\hspace{0.17em}}\left[{}^{12}{K}^{e}\right]={\left[{}^{21}{K}^{e}\right]}^{\text{T}}\\ {}^{22}{K}_{ij}^{e}={}^{22}K{}_{ji}^{e}\text{\hspace{0.17em}};\text{\hspace{1em}}i,j=1,2,\cdots ,{n}^{u},\text{\hspace{0.17em}}\text{hence}\text{\hspace{0.17em}}\left[{}^{22}{K}^{e}\right]\text{\hspace{0.17em}}\text{is}\text{\hspace{0.17em}}\text{symmetric}\end{array}$ (52)

From (52) we can conclude that $\left[{K}^{e}\right]$ in (47) is symmetric. For the entire discretization we can write

$\begin{array}{c}\left\{\begin{array}{c}{\left({A}_{1}\left({u}_{h}^{e},{v}_{h}^{e}\right),{w}_{1}\right)}_{{\stackrel{\xaf}{\Omega}}_{xy}^{e}}\\ {\left({A}_{2}\left({u}_{h}^{e},{v}_{h}^{e}\right),{w}_{2}\right)}_{{\stackrel{\xaf}{\Omega}}_{xy}^{e}}\end{array}\right\}={\displaystyle \underset{e}{\sum}}\left\{\begin{array}{c}{\left({A}_{1}\left({u}_{h}^{e},{v}_{h}^{e}\right),{w}_{1}\right)}_{{\stackrel{\xaf}{\Omega}}_{xy}^{e}}\\ {\left({A}_{2}\left({u}_{h}^{e},{v}_{h}^{e}\right),{w}_{2}\right)}_{{\stackrel{\xaf}{\Omega}}_{xy}^{e}}\end{array}\right\}\\ =-{\displaystyle \underset{e}{\sum}}\left(\left[{K}^{e}\right]\left\{{\delta}^{e}\right\}+\left\{{P}^{e}\right\}\right)=0\end{array}$ (53)

Hence

$\underset{e}{\sum}}\left(\left[{K}^{e}\right]\left\{{\delta}^{e}\right\}\right)={\displaystyle \underset{e}{\sum}}\left\{{P}^{e}\right\$ (54)

or

$\left[K\right]\left\{\delta \right\}=\left\{P\right\}$ (55)

in which

$\left[K\right]={\displaystyle \sum}\left[{K}^{e}\right];\text{\hspace{1em}}\text{assemblyofelementequations}$ (56)

$\left\{\delta \right\}={\displaystyle {\cup}_{e}\left\{{\delta}^{e}\right\}}$ (57)

$\left\{P\right\}={\displaystyle \underset{e}{\sum}}\left\{{P}^{e}\right\}$ (58)

4. Approximation Spaces and Some Remarks

1) Since the mathematical model ((12) and (13) contains up to fourth order derivatives of the displacements, the approximation functions in spaces ${V}_{h}\subset {H}^{k\mathrm{,}p}\left({\stackrel{\xaf}{\Omega}}_{xy}^{e}\right),k\ge 5$ are admissible in (12) and (13) and $k=5$ i.e. local approximation of class ${C}^{4}\left({\stackrel{\xaf}{\Omega}}_{xy}^{e}\right)$ corresponds to minimally conforming space.

2) Weak form (32) resulting from GM/WF only contains derivatives of up to order two of u and v, hence it is tempting to use ${u}_{h}^{e}$ and ${v}_{h}^{e}$ of class ${C}^{1}\left({\stackrel{\xaf}{\Omega}}_{xy}^{e}\right)$ but in doing so we rely on weak convergence of the solutions of class ${C}^{1}$ to class ${C}^{2}$ and eventually to class ${C}^{4}$ needed for the mathematical model.

3) Numerical values of the coefficients of $\left[{K}^{e}\right]$ are obtained using Gauss quadrature.

4) Solution is computed using assembled equations (55) for ${\stackrel{\xaf}{\Omega}}_{xy}^{T}$ after imposing boundary conditions.

5) Linearity of the algebraic system and symmetry $\left[{K}^{e}\right]$ and $\left[K\right]$ are due to the fact that the differential operator in (12) and (13) is linear in displacements u and v and the adjoint ${A}^{\mathrm{*}}$ of the differential operator $A$ is same as the operator $A$ (when the mathematical model is expressed in displacements u and v)

6) In the study of the model problem we chose $\beta =0$ (based on the material presented in [45] ) i.e. we consider balance of moments of moments as a balance law, hence the Cauchy moment tensor is symmetric.

5. A Least Squares Formulation in ${\mathbb{R}}^{2}$ (Plane Stress) Based on Residual Functional

We consider the following mathematical model (obtained using (1)-(10)) in the dimensionless form (in the absence of balance of moments of moments as a balance law [45] ) consisting of first order partial differential equations.

$\begin{array}{l}\frac{{\partial}_{s}{\sigma}_{xx}}{\partial x}+\frac{{\partial}_{s}{\sigma}_{yx}}{\partial y}+\frac{{\partial}_{a}{\sigma}_{yx}}{\partial y}=0;\text{\hspace{1em}}\frac{{\partial}_{s}{\sigma}_{xy}}{\partial x}+\frac{{\partial}_{s}{\sigma}_{yy}}{\partial y}-\frac{{\partial}_{a}{\sigma}_{yx}}{\partial x}=0\\ \frac{\partial {m}_{xz}}{\partial x}+\frac{\partial {m}_{yz}}{\partial y}+2\left({}_{a}{\sigma}_{yx}\right)=0\end{array}$ (59)

$\begin{array}{l}{}_{s}{\sigma}_{xx}={D}_{11}\frac{\partial u}{\partial x}+{D}_{12}\frac{\partial u}{\partial y};{\text{\hspace{1em}}}_{s}{\sigma}_{yy}={D}_{21}\frac{\partial u}{\partial x}+{D}_{22}\frac{\partial v}{\partial y};\\ {}_{s}{\sigma}_{xy}={D}_{33}\left(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\right)\end{array}$ (60)

$\begin{array}{l}{}_{s}{m}_{xz}=\left(\frac{{E}_{0}}{{m}_{0}{L}_{0}}\right)\underset{\u02dc}{\alpha}\frac{\partial \left({}_{i}{\Theta}_{z}\right)}{\partial x};{\text{\hspace{1em}}}_{s}{m}_{yz}=\left(\frac{{E}_{0}}{{m}_{0}{L}_{0}}\right)\underset{\u02dc}{\alpha}\frac{\partial \left({}_{i}{\Theta}_{z}\right)}{\partial y}\\ {}_{a}{m}_{xz}=-\left(\frac{{E}_{0}}{{m}_{0}{L}_{0}}\right)\underset{\u02dc}{\beta}\frac{\partial \left({}_{i}{\Theta}_{z}\right)}{\partial x};{\text{\hspace{1em}}}_{a}{m}_{yz}=-\left(\frac{{E}_{0}}{{m}_{0}{L}_{0}}\right)\underset{\u02dc}{\beta}\frac{\partial \left({}_{i}{\Theta}_{z}\right)}{\partial y}\end{array}$ (61)

$\begin{array}{l}{D}_{11}={D}_{22}=\frac{E}{1-{\nu}^{2}};\text{\hspace{1em}}{D}_{12}={D}_{21}=\frac{\nu E}{1-{\nu}^{2}};\text{\hspace{1em}}{D}_{33}=G=\frac{E}{2\left(1+\nu \right)}\\ {}_{i}{\Theta}_{z}=\frac{1}{2}\left(\frac{\partial u}{\partial y}-\frac{\partial v}{\partial x}\right)\end{array}$ (62)

In (59)-(62), $\frac{{E}_{0}}{{m}_{0}{L}_{0}}$ is in fact one, but it has been left in constitutive theory

for the moment tensors for sake of clarity. Equations (59)-(62) are a system of eleven first order linear coupled differential equations in eleven dependent variables $u$ , $v$ , ${}_{s}{\sigma}_{xx}$ , ${}_{s}{\sigma}_{yy}$ , ${}_{s}{\sigma}_{xy}$ , ${}_{a}{\sigma}_{yx}$ , ${}_{s}{m}_{xz}$ , ${}_{s}{m}_{yz}$ , ${}_{a}{m}_{xz}$ , ${}_{a}{m}_{yz}$ and ${}_{i}{\Theta}_{z}$ . A least squares formulation (LSF) of (59)-(62) is constructed using residual functionals [57] [61] - [66] resulting from each of the eleven equations when the local approximations for the dependent variables are substituted in them. The local approximations considered in higher order scalar product space ${H}^{k\mathrm{,}p}\left({\stackrel{\xaf}{\Omega}}_{xy}^{e}\right)$ , ${\stackrel{\xaf}{\Omega}}_{xy}^{e}$ being an element of the discretization which are p-version hierarchical with higher order global differentiability. Since (59)-(62) are a system of first order equations $k=2$ i.e. local approximations of class ${C}^{1}\left({\Omega}_{xy}^{e}\right)$ for each variable constitute minimally conforming space of approximations [57] . However, for the model problems considered here the solutions are sufficiently smooth, thus permitting the use of ${C}^{0}\left({\stackrel{\xaf}{\Omega}}_{xy}^{e}\right)$ local approximations with weak convergence to ${C}^{1}\left({\stackrel{\xaf}{\Omega}}_{xy}^{e}\right)$ .

6. Model Problems

In this section we can consider three model problems in ${\mathbb{R}}^{2}$ : 1) Simply supported thin plate with transverse in plane loading. 2) fixed-fixed thin plate with transverse in plane loading. 3) a square plate with a circular hole at the center subjected to uniaxial uniform loading.

Remarks.

・ In all numerical studies we consider both formulations, GM/WF as well as LSP.

・ We choose $\underset{\u02dc}{\beta}=0$ in all studies [45] which implies that ${}_{a}m=0$ and $m={}_{s}m$ implying that balance of moments of moments is a balance law. This is necessary for incorporating correct polar physics due to internal rotations in the mathematical model [45] . Thus, ${}_{a}{m}_{xz}=0$ and ${}_{a}{m}_{yz}=0$ in (59)-(62) and the mathematical model reduces to nine partial differential equations in nine dependent variables.

・ We note that the integral form in GM/WF contains upto second order derivatives of u and v, hence $k=3$ is minimally conforming approximation space (i.e. solutions of class ${C}^{2}$ in x and y) for the integral forms for which all integrals over the spatial discretization are Riemann. On the other hand for $k=2$ i.e ${C}^{1}$ approximations in x and y, the integrals over the spatial discretization are Lebesgue. For simply supported and fixed-fixed plate we consider numerical studies with $k=3,p=5$ (i.e. ${C}^{2}$ local approximations in x and y with p-level of 5) and also with $k=2,p=7$ , i.e. ${C}^{1}$ local approximations with p-level of 7. In case of square plate with a hole we consider $k=2$ with p-level of 7.

・ Computations for least squares formulation are only performed and compared with those from GM/WF for the simply supported and fixed-fixed plate to ensure that the solutions obtained using GM/WF in fact have the desired accuracy for the choices of k and p. In these studies we choose $k=1,p=9$ i.e. solutions of class ${C}^{0}$ with p-level of nine as used in references [40] [45] . For $k=1$ integrals over the discretization is in Lebesgue sense.

6.1. Simply Supported and Fixed-Fixed Plate: Model Problems 1 and 2

We consider a thin plate of length $\stackrel{^}{l}=20$ inches with width $\stackrel{^}{h}=0.5$ inches and thickness $\stackrel{^}{t}=0.1$ inches. With ${L}_{0}=10$ inches, the dimensionless plate is $L\times h\times t=2\times 0.05\times 0.01$ . Figure 1(a) and Figure 1(b) show schematics of the plate, boundary conditions and loading for the formulation based on GM/WF for both simply supported and fixed-fixed boundary conditions. The load is applied over a length of $b=0.4$ as forces at the nodes that corresponds to uniform stress in the y-direction (see Figure 1). Figure 2(a) and Figure 2(b) show same schematics with BCs and loading used in least squares formulation. In all numerical studies the plates are discretized using a 20 element uniform discretization (10 elements along the length and two elements width b) using a nine node p-version hierarchical higher order global differentiability finite elements. In all computations we choose Poisson’s ratio of 0.3, $\stackrel{^}{E}={E}_{0}=30\times {10}^{6}$ psi, hence $E=1$ , and $0\le \underset{\u02dc}{\alpha}\le 0.25$ with $\underset{\u02dc}{\beta}=0.0$ . $\underset{\u02dc}{\alpha}=0$ corresponds to classical continuum theory. Progressively increasing values of $\underset{\u02dc}{\alpha}$ produces progressively more pronounced polar physics (resistance to deformation). Numerical solutions are calculated for the following:

1) For GM/WF we consider $k=3$ (solutions of class ${C}^{2}$ in x and y) with $p=5$ . For this choice of k, integrals over the spatial discretization are Riemann.

2) For GM/WF we also consider $k=2$ (solutions of class ${C}^{1}$ in x and y) with $p=7$ . For this choice of k integrals over the spatial discretization are in Lebesgue sense.

3) Since the solution for LSP yields residual functional values of the order of $O\left({10}^{-15}\right)$ or lower, comparing the computed solutions from (1) and (2) we confirm that when both solutions are almost indistinguishable from each other, the solution from GM/WF has good accuracy.

4) For Least squares formulation we consider solutions of class ${C}^{0}$ in x and y with p-level of nine [40] [45] .

Results

GM/WF:

Figures 3(a)-(c) shows plots of v, ${}_{i}\Theta {}_{z}$ and ${}_{s}{m}_{xz}$ versus x at $y=0.025$ (center line of the plate) for $0\le \underset{\u02dc}{\alpha}\le 0.25$ . For $\underset{\u02dc}{\alpha}=0$ we have classical continuum behavior. Progressively increasing values of $\underset{\u02dc}{\alpha}$ results in progressively increasing resistance to deformation, hence progressively reducing displacement v, reducing rotation ${}_{i}{\Theta}_{z}$ but increasing moment ${}_{s}{m}_{xz}$ . Similar graphs of v, ${}_{i}{\Theta}_{z}$ and ${}_{s}{m}_{xz}$ versus x at $y=0.025$ for fixed-fixed plate are shown in Figures 4(a)-(c) for $0\le \underset{\u02dc}{\alpha}\le 0.25$ . We observe similar behaviors of v, ${}_{i}{\Theta}_{z}$ and ${}_{s}{m}_{xz}$

Figure 1. Model Problem 1 and 2: Schematics, BCs and loading (dimensionless): GM/WF. (a) Simply suppported plate ( $b=0.4$ ); (b) fixed-fixed plate ( $b=0.4$ ).

Figure 2. Model Problem 1 and 2: Schematics, BCs and loading (dimensionless): LSP. (a) Simply suppported plate ( $b=0.4$ ); (b) fixed-fixed plate ( $b=0.4$ ).

with increasing $\alpha $ values. Due to clamped boundaries the displacement v is significantly reduced and ${}_{i}{\Theta}_{z}$ and ${}_{s}{m}_{xz}$ follow accordingly.

Comparison of results: GM/WF and LSP:

The numerical solutions obtained from LS formulation for exactly same BCs and loading (Figure 2) using local approximations of class ${C}^{0}$ at p-level 9 are

Figure 3. $v\mathrm{,}{}_{i}\Theta {}_{z}$ and ${}_{s}{m}_{xz}$ versus $y=0.025$ : Simply supported plate (GM/WF). (a) v versus x at $y=0.025$ ; (b) ${}_{i}{\Theta}_{z}$ versus x at $y=0.025$ ; (c) ${}_{s}{m}_{xz}$ versus x at $y=0.025$ .

Figure 4. $v\mathrm{,}{}_{i}\Theta {}_{z}$ and ${}_{s}{m}_{xz}$ versus $y=0.025$ : Fixed-fixed plate (GM/WF). (a) v versus x at $y=0.025$ ; (b) ${}_{i}{\Theta}_{z}$ versus x at $y=0.025$ ; (c) ${}_{s}{m}_{xz}$ versus x at $y=0.025$ .

Figure 5. Simply supported plate: comparison of LSP and GM/WF. (a) v versus x at $y=0.025$ ; (b) ${}_{i}{\Theta}_{z}$ versus x at $y=0.025$ ; (c) ${}_{s}{m}_{xz}$ versus x at $y=0.025$ .

compared with those obtained using GM/WF ( ${C}^{2}$ solutions at $p=5$ or ${C}^{1}$ solutions at $p=7$ ). In the LSP the residual functional for the discretization is of the order $O\left({10}^{-15}\right)$ . This ensures that the computed solutions satisfy the governing differential equations in the pointwise sense, hence the computed solutions are virtually same as the theoretical solution. Comparison of these solutions with GM/WF provides a check on the accuracy of the solutions obtained using GM/WF as in GM/WF there is no direct measure of accuracy in the method itself.

Figures 5(a)-(c) show the plots of v, ${}_{i}{\Theta}_{z}$ and ${}_{s}{m}_{xz}$ versus x at $y=0.025$ for $\underset{\u02dc}{\alpha}=0$ , 0.001, and 0.1 ( $\underset{\u02dc}{\alpha}=0$ being the classical theory) obtained using GM/WF and a comparison with least squares method for simple supported plate. Similar results for GM/WF and a comparison with LSP for fixed-fixed plate are shown in Figures 6(a)-(c). In both Figure 5 and Figure 6, v, ${}_{i}{\Theta}_{z}$ and ${}_{s}{m}_{xz}$ obtained using GM/WF and LSP are in perfect agreement with each other for all three values of $\underset{\u02dc}{\alpha}$ .

Figure 6. Fixed-fixed plate: comparison of LSP and GM/WF. (a) v versus x at $y=0.025$ ; (b) ${}_{i}{\Theta}_{z}$ versus x at $y=0.025$ ; (c) ${}_{s}{m}_{xz}$ versus x at $y=0.025$ .

6.2. A Square Plate with a Circular Hole: Model Problem 3

We consider a 6" × 6" square plate of thickness 0.1" with a 0.48" diameter circular hole at the center. We use ${L}_{0}=1.5"$ . The material properties, reference quantities etc. used here are same as for model problems 1 and 2. Poisson’s ratio of 0.3 is used. This gives rise to a $4\times 4\times 0.06$ dimensionless plate with a hole diameter of 0.16 (Figure 7(a)). The plate is subjected to uniform displacement of 0.01 (dimensionless) on its vertical faces that creates a uniform dimensionless stress field of ${\left({\sigma}_{xx}\right)}_{0}=0.0048$ . The details of the BCs and loading for quarter plate are shown in Figure 7(a). Figure 7(c) shows a graded discretization of the quarter plate. The plate is divided in four bicubic patches (Figure 7(b)). In each patch a 3 × 3 uniform discretization of nine-node p-version hierarchical elements with higher order global differentiability local approximation [59] [67] is used giving a total of 36 elements for the quarter of the plate. Computations are performed only using the formulation based on GM/WF with local approximation of class ${C}^{1}$ i.e $k=2$ in x and y space (for distorted elements in ${\mathbb{R}}^{2}$ , xy-space

Figure 7. Schematic of a quarter of the plate with a circular hole, BCs, loading and finite element discretization.

[57] [67] ) with p-level of 7. For this choice of k, the integrals over the discretization are Lebesgue, but due to smoothness of the solution we can expect these solutions to converge to class ${C}^{2}$ in the weak sense.

Results and Discussion

The stresses ${}_{s}{\sigma}_{xx}$ and ${}_{s}{\sigma}_{yy}$ are normalized using ${\left({\sigma}_{xx}\right)}_{0}$ , i.e. ${\left({}_{s}{\sigma}_{xx}\right)}_{n}=$ ${}_{s}{\sigma}_{xx}/{\left({\sigma}_{xx}\right)}_{0}$ and ${\left({}_{s}{\sigma}_{yy}\right)}_{n}={}_{s}{\sigma}_{yy}/{\left({\sigma}_{xx}\right)}_{0}$ . It is well known that based on classical continuum theory (when $\underset{\u02dc}{\alpha}=0$ ) we have stress concentration of 3.0 at E (Figure 7(a)) i.e. in this case we expect ${\left({}_{s}{\sigma}_{xx}\right)}_{n}=3.0$ at point E (Figure 7(a)). With increasing values of $\underset{\u02dc}{\alpha}$ increasing presence of internal polar physics is present that results in progressively increasing resistance to deformation. As a consequence stresses increase as $\underset{\u02dc}{\alpha}$ increases. Figure 8(a) shows plots of

Figure 8. Normalized stress versus y at $x=\mathrm{0.0.}$ (a) Normalized stress versus y at $x=0.0$ ; (b) Exploded view in the vicinity of hole.

${\left({}_{s}{\sigma}_{xx}\right)}_{n}$ versus y along the the edge of ED of the plate for different values of $\underset{\u02dc}{\alpha}$ . At D, ${\left({}_{s}{\sigma}_{xx}\right)}_{n}=1$ as expected. As we approach E from D, stress ${\left({}_{s}{\sigma}_{xx}\right)}_{n}$ increases. The exploded view in the vicinity of point E shown in Figure 8(b) confirms that when $\underset{\u02dc}{\alpha}=0$ i.e. the classical theory, ${\left({}_{s}{\sigma}_{xx}\right)}_{n}$ is indeed 3.0. With increasing $\underset{\u02dc}{\alpha}$ (for the range considered here), ${\left({}_{s}{\sigma}_{xx}\right)}_{n}$ as high as 3.5 is obtained. From Figure 8(a) we note that the for progressively increasing values of $\underset{\u02dc}{\alpha}$ , ${\left({}_{s}{\sigma}_{xx}\right)}_{n}$ progressively increases from a value of 1.0 at D to upto 3.5 at E

Figure 9. Contour plots showing quarter of a plate with a circular hole under uniaxial loading. (a) classical; (b) $\underset{\u02dc}{\alpha}=0.0001$ ; (c) $\underset{\u02dc}{\alpha}=\mathrm{0.000185.}$

for the largest value of $\underset{\u02dc}{\alpha}$ ( $\underset{\u02dc}{\alpha}=0.000185$ ) used in the studies. Figures 9(a)-(c) show carpet plots of ${\left({}_{s}{\sigma}_{xx}\right)}_{n}$ for $\underset{\u02dc}{\alpha}=0.0,0.0001,0.000185$ . With progressively increasing values of $\underset{\u02dc}{\alpha}$ , higher values of ${\left({}_{s}{\sigma}_{xx}\right)}_{n}$ in the entire quarter of the plate are observed compared to classical theory ( $\underset{\u02dc}{\alpha}=0.0$ ), most significant increase in ${\left({}_{s}{\sigma}_{xx}\right)}_{n}$ being at point E as expected.

7. General Remarks

In this section we make some remarks related to the two finite element formulations (GM/WF, LSP) in context with the numerical studies presented here.

1) It is obvious that for the model problems (in ${\mathbb{R}}^{2}$ ) the GM/WF has only two dependent variables u and v whereas LSP based on first order system of PDEs has nine dependent variables resulting in enormous computational inefficiency but permitting flexibility that permits ${C}^{0}$ local approximations.

2) In LSP there is no concept of secondary variables as in GM/WF, hence there are no self equilibrating quantities in LS finite element formulation. As a result, all dependent variables pertaining to the known physics, even those that are zero, must be specified on the boundaries of the domain. For example in GM/WF stress free boundaries are automatically satisfied due to sum of secondary variables being zero at a node. Same is true for moment free boundaries. However, in LSP all boundary information must be defined in the problem data. Figure 1(a), Figure 1(b) for GM/WF and Figure 2(a), Figure 2(b) for LSP containing schematics and BCs for model problem 1 and 2 clearly illustrate this. Due to the necessity of defining all dependent variable specifications on the boundaries of the domain in LSP as BCs, often the specifications become cumbersome and result in redundancies in their specifications. GM/WF is completely free of such problems. In case of square plate with a circular hole, the situation is much more difficult in LSP as in this case the stress and moment components normal to the hole are zero while the tangential components need to be computed. Definition of such BCs require either constrained equations or a rotated local coordinate system on the hole boundary that is normal and tangent to the hole boundary. In GM/WF normal tractions (both stress and moment) are secondary variables, hence their sum on the hole boundary in naturally zero at each node, thus this boundary condition is automatically satisfied.

3) In GM/WF as well as in LSP we have considered integrals over the spatial discretization in Lebesgue sense, but there is no issue of their convergence. LS residual functional $O\left({10}^{-15}\right)$ and perfect match of GM/WF results with LSP for model problems 1 and 2 confirm that both GM/WF and LSP results are sufficiently converged to be as good as theoretical solutions.

8. Summary and Conclusions

In this paper the mathematical model consisting of conservation and balance laws in Lagrangian description for non-classical continuum theory for elastic solids (small strain small deformation physics without dissipation and memory) incorporating internal rotation physics due to displacement gradient tensor is considered (derived in reference [45] ). In such solids the deformation physics due to mechanical work is reversible; hence the differential operator $A$ in this mathematical model when expressed purely in terms of displacements is such that the adjoint ${A}^{\mathrm{*}}$ of the differential operator $A$ is same as $A$ . Thus, in such mathematical models GM/WF is ideal for the finite element formulation of the corresponding BVPs. We make the following specific remarks and observations and draw some conclusions from the work presented in this paper.

1) GM/WF is ideal for reversible processes as in the present case. In such mathematical models ${A}^{\mathrm{*}}=A$ holds.

2) LSP with first order system of PDEs is computationally non competitive with GM/WF. In the work presented here GM/WF has only two dependent variables whereas LSP has nine.

3) An important question is “could we have used LSP” for the mathematical model in displacements u and v derived for GM/WF. Of course we could but: (1) this would require solutions of class ${C}^{4}$ or of class ${C}^{3}$ for sure (2) stress and moment boundary conditions (zero or non zero) are extremely difficult to define as the stresses and moments are not dependent variables any more in the mathematical model. (3) Due to lack of secondary variable, zero stress and moment boundary conditions also need to be specified. Due to these difficulties it is perhaps more convenient to use mathematical models consisting of first order PDEs in LSP which are computationally not competitive with GM/WF.

4) Numerical studies for the three model problems clearly demonstrate superiority of GM/WF over LSP in almost all aspects.

5) Numerical solutions computed using GM/WF and those using LSP satisfy PDEs, almost in the pointwise sense as the residual functional for the discretization is $O\left({10}^{-15}\right)$ .

6) Presence of increasing polar physics with increasing $\underset{\u02dc}{\alpha}$ is clearly demonstrated in model problem 1 and 2 (also shown in references [40] [45] using finite element formulation based on LSP) using both finite element formulations.

7) The third model problem is rather difficult to study using finite element formulation based on LSP due to the difficulty of specifying zero boundary conditions of stress and moment normal to the hole boundary. In GM/WF the secondary variables and their sum being zero on free boundaries automatically satisfy these BCs.

8) In the plate problem with a circular hole the stress concentration at point E of 3.0 is predicted correctly when $\underset{\u02dc}{\alpha}=0$ . With progressively increasing $\underset{\u02dc}{\alpha}$ , the stress concentration at E increases from 3 to 3.5 for the largest value of $\underset{\u02dc}{\alpha}=0.000185$ used here.

9) The finite element formulation based on GM/WF for non-classical continuum models is a valuable approach in which the deformation due to mechanical work is reversible. The finite element formulation based on GM/WF for non-classical continuum models presented here is superior and meritorious in all aspects when compared to the finite element formulations based on least square processes. The only disadvantage one could possibly point out is the presence of up to fourth order derivatives of displacements in the mathematical models. In view of the research work on k-version [58] [59] [60] this is hardly of any consequence.

Acknowledgements

The support provided by the first and the third author’s university distinguised professorships is gratefully acknowledged. The financial support provided to the second author by the department of mechanical engineering and school of engineering is also acknowledged. The computational infrastructure of CML of the mechanical engineering department have been instrumental in performing the numerical studies presented in the paper.

References

[1] Eringen, A.C. (1964) Mechanics of Micromorphic Materials. In: Gortler, H., Ed., Applied Mechanics, 131-138.

[2] Eringen, A.C. (1968) Mechanics of Micromorphic Continua. In: Kroner, E., Ed., Mechanics of Generalized Continua, 18-35.

https://doi.org/10.1007/978-3-662-30257-6_2

[3] Eringen, A.C. (1968) Theory of Micropolar Elasticity. In: Liebowitz, H., Ed., Fracture, 621-729.

[4] Eringen, A.C. (1970) Balance Laws of Micromorphic Mechanics. International Journal of Engineering Science, 8, 819-828.

[5] Eringen, A.C. (1990) Theory of Thermo-Microstretch Fluids and Bubbly Liquids. International Journal of Engineering Science, 28, 133-143.

[6] Eringen, A.C. (1999) Theory of Micropolar Elasticity. In: Microcontinuum Field Theories, Springer-Verlag, New York, 101-248.

[7] Eringen, A.C. (1966) A Unified Theory of Thermomechanical Materials. International Journal of Engineering Science, 4, 179-202.

[8] Eringen, A.C. (1967) Linear Theory of Micropolar Viscoelasticity. International Journal of Engineering Science, 5, 191-204.

[9] Eringen, A.C. (1972) Theory of Micromorphic Materials with Memory. International Journal of Engineering Science, 10, 623-641.

[10] Koiter, W. (1964) Couple Stresses in the Theory of Elasticity, I & II. Philosophical Transactions of the Royal Society of London B, 67, 17-44.

[11] Oevel, W. and Schröter, J. (1981) Balance Equations for Micromorphic Materials. Journal of Statistical Physics, 25, 645-662.

https://doi.org/10.1007/BF01022359

[12] Reddy, J.N. (2007) Nonlocal Theories for Bending, Buckling and Vibration of Beams. International Journal of Engineering Science, 45, 288-307.

[13] Reddy, J.N. and Pang, S.D. (2008) Nonlocal Continuum Theories of Beams for the Analysis of Carbon Nanotubes. Journal of Applied Physics, 103, Article ID: 023611.

https://doi.org/10.1063/1.2833431

[14] Reddy, J.N. (2010) Nonlocal Nonlinear Formulations for Bending of Classical and Shear Deformation Theories of Beams and Plates. International Journal of Engineering Science, 48, 1507-1518.

[15] Lu, P., Zhang, P.Q., Leo, H.P., Wang, C.M. and Reddy, J.N. (2007) Nonlocal Elastic Plate Theories. Proceedings of Royal Society A, 463, 3225-3240.

https://doi.org/10.1098/rspa.2007.1903

[16] Yang, J.F.C. and Lakes, R.S. (1982) Experimental Study of Micropolar and Couple Stress Elasticity in Compact Bone in Bending. Journal of Biomechanics, 15, 91-98.

[17] Lubarda, V.A. and Markenscoff, X. (2000) Conservation Integrals in Couple Stress Elasticity. Journal of Mechanics and Physics of Solids, 48, 553-564.

[18] Ma, H.M., Gao, X.L. and Reddy, J.N. (2008) A Microstructure-Dependent Timoshenko Beam Model Based on a Modified Couple Stress Theory. Journal of Mechanics and Physics of Solids, 56, 3379-3391.

[19] Ma, H.M., Gao, X.L. and Reddy, J.N. (2010) A Nonclassical Reddy-Levinson Beam Model Based on Modified Couple Stress Theory. Journal of Multiscale Computational Engineering, 8, 167-180.

https://doi.org/10.1615/IntJMultCompEng.v8.i2.30

[20] Reddy, J.N. (2011) Microstructure Dependent Couple Stress Theories of Functionally Graded Beams. Journal of Mechanics and Physics of Solids, 59, 2382-2399.

[21] Reddy, J.N. and Arbind, A. (2012) Bending Relationship between the Modified Couple Stress-Based Functionally Graded Timoshenko Beams and Homogeneous Bernoulli-Euler Beams. Annals of Solid Structural Mechanics, 3, 15-26.

https://doi.org/10.1007/s12356-012-0026-z

[22] Srinivasa, A.R. and Reddy, J.N. (2013) A Model for a Constrained, Finitely Deforming Elastic Solid with Rotation Gradient Dependent Strain Energy and Its Specialization to Von Kármán Plates and Beams. Journal of Mechanics and Physics of Solids, 61, 873-885.

[23] Mora, R.J. and Waas, A.M. (2007) Evaluation of the Micropolar Elasticity Constants for Honeycombs. Acta Mechanica, 192, 1-16.

https://doi.org/10.1007/s00707-007-0446-8

[24] Onck P.R. (2002) Cosserat Modeling of Cellular Solids. C. R. Mecanique, 330, 717-722.

[25] Segerstad, P.H., Toll, S. and Larsson, R. (2009) A Micropolar Theory for the Finite Elasticity of Open-Cell Cellular Solids. Proceedings of Royal Society A, 465, 843-865.

https://doi.org/10.1098/rspa.2008.0267

[26] Grioli, G. (2003) Microstructures as a Refinement of Cauchy Theory. Problems of Physical Concreteness, Continuum Mechanics and Thermodynamics, 15, 441-450.

https://doi.org/10.1007/s00161-003-0122-8

[27] Grekova, E.F. and Maugin, G.A. (2005) Modelling of Complex Elastic Crystals by Means of Multi-Spin Micromorphic Media. International Journal of Engineering Science, 43, 494-519.

[28] Altenbach, J., Altenbach, H. and Eremeyev, V.A. (2010) On Generalized Cosserat-Type Theories of Plates and Shells: A Short Review and Bibliography, Archive of Applied Mechanics, 80, 217-227.

https://doi.org/10.1007/s00419-009-0314-1

[29] Lazar, M. and Maugin, G.A. (2005) Nonsingular Stress and Strain Fields of Dislocations and Disclinations in First Strain Gradient Elasticity. International Journal of Engineering Science, 43, 1157-1184.

[30] Lazar, M. and Maugin, G.A. (2004) Defects in Gradient Micropolar Elasticity? I: Screw Dislocation. Journal of the Mechanics and Physics of Solids, 52, 2263-2284.

[31] Maugin, G.A. (1978) A Phenomenological Theory of Ferroliquids. International Journal of Engineering Science, 16, 1029-1044.

[32] Maugin, G.A. (1981) Wave Motion in Magnetizable Deformable Solids. International Journal of Engineering Science, 19, 321-388.

[33] Maugin, G.A. (1988) Wave Motion in Magnetizable Deformable Solids. Philosophical Transactions of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences, 356, 1367-1395.

https://doi.org/10.1098/rsta.1998.0226

[34] Cosserat, E. and Cosserat, F. (1909) Théorie des corps déformables. Paris.

[35] Green, A.E. (1965) Micro-Materials and Multipolar Continuum Mechanics. International Journal of Engineering Science, 3, 533-537.

[36] Green, A.E. and Rivlin, R.S. (1967) The Relation between Director and Multipolar Theories in Continuum Mechanics. Zeitschrift für angewandte Mathematik und Physik ZAMP, 18, 208-218.

[37] Green, A.E. and Rivlin, R.S. (1964) Multipolar Continuum Mechanics. Archive for Rational Mechanics and Analysis, 17, 113-147.

https://doi.org/10.1007/BF00253051

[38] Surana, K.S., Powell, M.J. and Reddy, J.N. (2015) A More Complete Thermodynamic Framework for Solid Continua. Journal of Thermal Engineering, 1, 1-13.

[39] Surana, K.S., Powell, M.J., Nunez, D. and Reddy, J.N. (2015) A Polar Continuum Theory for Solid Continua. International Journal of Engineering Research and Industrial Applications, 8, 77-106.

[40] Surana, K.S., Powell, M.J. and Reddy, J.N. (2015) Constitutive Theories for Internal Polar Thermoelastic Solid Continua. Journal of Pure and Applied Mathematics: Advances and Applications, 14, 89-150.

[41] Surana, K.S., Joy, A.D. and Reddy, J.N. (2016) A Non-Classical Internal Polar Continuum Theory for Finite Deformation of Solids Using First Piola-Kirchhoff Stress Tensor. Journal of Pure and Applied Mathematics: Advances and Applications, 16, 1-41.

[42] Surana, K.S., Joy, A.D. and Reddy, J.N. (2016) Non-Classical Internal Polar Continuum Theory for Finite Deformation and Finite Strains in Solids. International Journal of Pure and Applied Mathematics, 4, 59-97.

[43] 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, 665-698.

https://doi.org/10.1007/s00161-017-0554-1

[44] Surana, K.S., Joy, A.D. and Reddy, J.N. (2017) Non-Classical Continuum Theory for Fluids Incorporating Internal and Cosserat Rotation Rates. Continuum Mechanics and Thermodynamics, 1-41.

https://doi.org/10.1007/s00161-017-0579-5

[45] Surana, K.S., Shanbhag, R.R. and Reddy, J.N. (2017) Necessity of Balance of Moments of Moments Balance Law in Non-Classical Continuum Theories for Solid Continua. Meccanica. (In Press)

[46] Gelfand, M. and Fomin, S. (2000) Calculus of Variations. Dover Publications.

[47] Hellwig, G. (1967) Differential Operators of Mathematical Physics. Addison-Wesley Publishing Co.

[48] Mikhlin, G. (1964) Variational Methods in Mathematical Physics. Pergamon Press.

[49] Reddy, J.N. (1986) Applied Functional Analysis and Variational Methods in Engineering. McGraw-Hill.

[50] Becker, M. (1964) The Principles and Applications of Variational Methods. MIT Press.

[51] Forray, M. (1968) Variational Calculus in Science and Engineering. McGraw-Hill.

[52] Rektorys, K. (1977) Variational Methods in Mathematics, Science and Engineering. Reidel.

https://doi.org/10.1007/978-94-011-6450-4

[53] Schechter, R.S. (1967) The Variational Methods in Engineering. McGraw-Hill.

[54] Washizu, K. (1982) Variational Methods in Elasticity and Plasticity. Pergamon Press.

[55] Weinstock, R. (1952) Calculus of Variations with Applications to Physics and Engineering. McGraw-Hill.

[56] Reddy, J.N. (1993) An Introduction to the Finite Element Method. McGraw-Hill.

[57] Surana, K.S. and Reddy, J.N. (2016) The Finite Element Method for Boundary Value Problems: Mathematics and Computations. CRC/Taylor & Francis.

[58] 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 Engineering Science, 3, 155-218.

https://doi.org/10.1142/S1465876302000605

[59] 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 Engineering Science, 4, 737-812.

https://doi.org/10.1142/S1465876303002179

[60] Surana, K.S., Ahmadi, A.R. and Reddy, J.N. (2004) The k-Version of Finite Element Method for Non-Linear Operators in BVP. International Journal of Computational Engineering Science, 5, 133-207.

https://doi.org/10.1142/S1465876304002307

[61] Winterscheidt, D. and Surana, K.S. (1993) p-Version Least-Squares Finite Element Formulation for Convection-Diffusion Problems. International Journal for Numerical Methods in Engineering, 36, 111-133.

https://doi.org/10.1002/nme.1620360107

[62] Winterscheidt, D. and Surana, K.S. (1993) p-Version Least-Squares Finite Element Formulation of Burgers’ Equation. International Journal for Numerical Methods in Engineering, 36, 3629-3646.

https://doi.org/10.1002/nme.1620362105

[63] Winterscheidt, D. and Surana, K.S. (1994) p-Version Least-Squares Finite Element Formulation for Two-Dimensional, Incompressible Fluid Flow. International Journal for Numerical Methods in Fluids, 18, 43-69.

https://doi.org/10.1002/fld.1650180104

[64] Bell, B.C. and Surana, K.S. (1994) p-Version Least Squares Finite Element Formulation for Two-Dimensional, Incompressible, Non-Newtonian Isothermal and Non-Isothermal Fluid Flow. International Journal for Numerical Methods in Fluids, 18, 127-162.

https://doi.org/10.1002/fld.1650180202

[65] Bell, B.C. and Surana, K.S. (1994) A Space-Time Coupled p-Version Least-Squares Finite Element Formulation for Unsteady Fluid Dynamics Problems. International Journal for Numerical Methods in Engineering, 37, 3545-3569.

https://doi.org/10.1002/nme.1620372008

[66] Bell, B.C. and Surana, K.S. (1996) A Space-Time Coupled p-Version Least Squares Finite Element Formulation for Unsteady Two-Dimensional Navier-Stokes Equations. International Journal for Numerical Methods in Engineering, 39, 2593-2618.

https://doi.org/10.1002/(SICI)1097-0207(19960815)39:15<2593::AID-NME968>3.0.CO;2-2

[67] Ahmadi, A.R., 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