Finite Element Processes Based on GM/WF in Non-Classical Solid Mechanics ()
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
(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
, we have constitutive theories for
and
. Choosing
,
,
,
we can write the following for balance laws and the constitutive theories
.
(1)
(2)
(3)
(4)
(5)
(6)
(7)
(8)
(9)
(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 (
) on all quantities and variables indicating that they all have their usual dimensions in terms of length (
), force (
), and time (
). If we choose
and
as reference length, force and time, then the dimensionless length, force and time (
and
) are defined as
(11)
If we consider
,
,
,
,
,
,
,
and choose
,
, then we obtain the dimensionless form of Equations (1)-(10). In these
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
,
,
,
,
,
,
,
,
,
and
. We substitute
from (3) into (1) and (2).
(12)
(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
and
consists of the decomposition in (3) and the constitutive theories for symmetric and antisymmetric Cauchy moment tensors are given by (5)-(8) with
defined by (9), hence there are also functions of the gradients of u and v. Let
be the discretization of
, domain of definition of mathematical model in which
is a typical finite element.
and
in (12) and (13) are differential operators that act on u and v (
and
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
(14)
in which
and
are approximations of u and v over
. Then, based on the fundamental lemma of the calculus of variations [46] - [60] we construct scalar products of (1) and (2) with test functions
and
over the discretization
and set them to zero.
(15)
or
(16)
In which
and
where
,
are local approximations of u and v over an element e with domain
. We consider
and
in which
,
and substitute for
and
from (12) and (13)
(17)
(18)
Integration by parts once for each term in (17) and (18) yields (noting that
,
being closed boundary of
)
(19)
(20)
Integrating by parts again for the moment terms in (19) and (20)
(21)
(22)
Let
(23)
Using (23) we can write (21) and (22) as follows
(24)
(25)
From (24) and (25) we can conclude that the primary and the secondary variables (PV and SV) are
Substituting for
,
,
from (4) into (24) and (25) and
from (9) into (5)-(8) and then (5)-(8) into (24) and (25) we can write
(26)
(27)
in which
(28)
(29)
(30)
(31)
Functionals
,
,
,
are linear in all of their arguments. We note that
,
are concomitants resulting only because of integration by parts. We can represent (28)-(31) in matrix and vector form
(32)
(33)
in which
(34)
(35)
(36)
(37)
The functions
have the following properties.
(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
and
will contain symmetric element coefficient matrix.
Let
be a nine node p-version hierarchical element with local approximation in higher order scalar product space
[57] [58] [59] [60] . Consider
, a map of
in
space in a two unit square. Then, we can write
(39)
(40)
and
are local approximation functions and
and
are corresponding nodal degrees of freedom for u and v. Using (39) and (40)
(41)
Let the total degrees of freedom for an element e be
(42)
Substituting from (39)-(41) into (34)-(37), we can write
(43)
(44)
(45)
(46)
Using (43)-(46) we can write (33) as follows
(47)
in which
(48)
(49)
(50)
(51)
We note that
(52)
From (52) we can conclude that
in (47) is symmetric. For the entire discretization we can write
(53)
Hence
(54)
or
(55)
in which
(56)
(57)
(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
are admissible in (12) and (13) and
i.e. local approximation of class
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
and
of class
but in doing so we rely on weak convergence of the solutions of class
to class
and eventually to class
needed for the mathematical model.
3) Numerical values of the coefficients of
are obtained using Gauss quadrature.
4) Solution is computed using assembled equations (55) for
after imposing boundary conditions.
5) Linearity of the algebraic system and symmetry
and
are due to the fact that the differential operator in (12) and (13) is linear in displacements u and v and the adjoint
of the differential operator
is same as the operator
(when the mathematical model is expressed in displacements u and v)
6) In the study of the model problem we chose
(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
(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.
(59)
(60)
(61)
(62)
In (59)-(62),
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
,
,
,
,
,
,
,
,
,
and
. 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
,
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
i.e. local approximations of class
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
local approximations with weak convergence to
.
6. Model Problems
In this section we can consider three model problems in
: 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
in all studies [45] which implies that
and
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,
and
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
is minimally conforming approximation space (i.e. solutions of class
in x and y) for the integral forms for which all integrals over the spatial discretization are Riemann. On the other hand for
i.e
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
(i.e.
local approximations in x and y with p-level of 5) and also with
, i.e.
local approximations with p-level of 7. In case of square plate with a hole we consider
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
i.e. solutions of class
with p-level of nine as used in references [40] [45] . For
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
inches with width
inches and thickness
inches. With
inches, the dimensionless plate is
. 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
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,
psi, hence
, and
with
.
corresponds to classical continuum theory. Progressively increasing values of
produces progressively more pronounced polar physics (resistance to deformation). Numerical solutions are calculated for the following:
1) For GM/WF we consider
(solutions of class
in x and y) with
. For this choice of k, integrals over the spatial discretization are Riemann.
2) For GM/WF we also consider
(solutions of class
in x and y) with
. 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
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
in x and y with p-level of nine [40] [45] .
Results
GM/WF:
Figures 3(a)-(c) shows plots of v,
and
versus x at
(center line of the plate) for
. For
we have classical continuum behavior. Progressively increasing values of
results in progressively increasing resistance to deformation, hence progressively reducing displacement v, reducing rotation
but increasing moment
. Similar graphs of v,
and
versus x at
for fixed-fixed plate are shown in Figures 4(a)-(c) for
. We observe similar behaviors of v,
and
Figure 1. Model Problem 1 and 2: Schematics, BCs and loading (dimensionless): GM/WF. (a) Simply suppported plate (
); (b) fixed-fixed plate (
).
Figure 2. Model Problem 1 and 2: Schematics, BCs and loading (dimensionless): LSP. (a) Simply suppported plate (
); (b) fixed-fixed plate (
).
with increasing
values. Due to clamped boundaries the displacement v is significantly reduced and
and
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
at p-level 9 are
Figure 3.
and
versus
: Simply supported plate (GM/WF). (a) v versus x at
; (b)
versus x at
; (c)
versus x at
.
Figure 4.
and
versus
: Fixed-fixed plate (GM/WF). (a) v versus x at
; (b)
versus x at
; (c)
versus x at
.
Figure 5. Simply supported plate: comparison of LSP and GM/WF. (a) v versus x at
; (b)
versus x at
; (c)
versus x at
.
compared with those obtained using GM/WF (
solutions at
or
solutions at
). In the LSP the residual functional for the discretization is of the order
. 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,
and
versus x at
for
, 0.001, and 0.1 (
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,
and
obtained using GM/WF and LSP are in perfect agreement with each other for all three values of
.
Figure 6. Fixed-fixed plate: comparison of LSP and GM/WF. (a) v versus x at
; (b)
versus x at
; (c)
versus x at
.
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
. 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
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
. 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
i.e
in x and y space (for distorted elements in
, 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
in the weak sense.
Results and Discussion
The stresses
and
are normalized using
, i.e.
and
. It is well known that based on classical continuum theory (when
) we have stress concentration of 3.0 at E (Figure 7(a)) i.e. in this case we expect
at point E (Figure 7(a)). With increasing values of
increasing presence of internal polar physics is present that results in progressively increasing resistance to deformation. As a consequence stresses increase as
increases. Figure 8(a) shows plots of
Figure 8. Normalized stress versus y at
(a) Normalized stress versus y at
; (b) Exploded view in the vicinity of hole.
versus y along the the edge of ED of the plate for different values of
. At D,
as expected. As we approach E from D, stress
increases. The exploded view in the vicinity of point E shown in Figure 8(b) confirms that when
i.e. the classical theory,
is indeed 3.0. With increasing
(for the range considered here),
as high as 3.5 is obtained. From Figure 8(a) we note that the for progressively increasing values of
,
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)
; (c)
for the largest value of
(
) used in the studies. Figures 9(a)-(c) show carpet plots of
for
. With progressively increasing values of
, higher values of
in the entire quarter of the plate are observed compared to classical theory (
), most significant increase in
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
) 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
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
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
in this mathematical model when expressed purely in terms of displacements is such that the adjoint
of the differential operator
is same as
. 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
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
or of class
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
.
6) Presence of increasing polar physics with increasing
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
. With progressively increasing
, the stress concentration at E increases from 3 to 3.5 for the largest value of
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.