Consistency and Validity of the Mathematical Models and the Solution Methods for BVPs and IVPs Based on Energy Methods and Principle of Virtual Work for Homogeneous Isotropic and Non-Homogeneous Non-Isotropic Solid Continua

Karan S. Surana, Emilio N. Alverio

Mechanical Engineering, University of Kansas, Lawrence, KS, USA.

DOI: 10.4236/am.2020.117039                

Abstract

Energy methods and the principle of virtual work are commonly used for obtaining solutions of boundary value problems (BVPs) and initial value problems (IVPs) associated with homogeneous, isotropic and non-homogeneous, non-isotropic matter without using (or in the absence of) the mathematical models of the BVPs and the IVPs. These methods are also used for deriving mathematical models for BVPs and IVPs associated with isotropic, homogeneous as well as non-homogeneous, non-isotropic continuous matter. In energy methods when applied to IVPs, one constructs energy functional (I) consisting of kinetic energy, strain energy and the potential energy of loads. The first variation of this energy functional (δI) set to zero is a necessary condition for an extremum of I. In this approach one could use δI = 0 directly in constructing computational processes such as the finite element method or could derive Euler’s equations (differential or partial differential equations) from δI = 0, which is also satisfied by a solution obtained from δI = 0. The Euler’s equations obtained from δI = 0 indeed are the mathematical model associated with the energy functional I. In case of BVPs we follow the same approach except in this case, the energy functional I consists of strain energy and the potential energy of loads. In using the principle of virtual work for BVPs and the IVPs, we can also accomplish the same as described above using energy methods. In this paper we investigate consistency and validity of the mathematical models for isotropic, homogeneous and non-isotropic, non-homogeneous continuous matter for BVPs that are derived using energy functional consisting of strain energy and the potential energy of loads. Similar investigation is also presented for IVPs using energy functional consisting of kinetic energy, strain energy and the potential energy of loads. The computational approaches for BVPs and the IVPs designed using energy functional and principle of virtual work, their consistency and validity are also investigated. Classical continuum mechanics (CCM) principles i.e. conservation and balance laws of CCM with consistent constitutive theories and the elements of calculus of variations are employed in the investigations presented in this paper.

Keywords

Energy Methods, Principle of Virtual Work, Calculus of Variations, Euler’s Equation, Mathematical Model, Classical and Non-Classical, Continuum Mechanics

1. Introduction

1.1. Literature Review

Even though this paper does not consider specific aspects of formulation and solution of composite mechanics problems, but such materials are non-homogeneous and non-isotropic, and the energy methods are almost exclusively employed in the study of BVPs and IVPs associated with such materials. It is for this reason that we present a review of published work related to composites based on either energy methods or principle of virtual work. Widespread use of synthesized composite materials is incentive to examine the consistency and legitimacy of the mathematical models for composites as well as the techniques of obtaining their solutions. Broadly, we can categorize composite material in two groups: 1) those that consist of a matrix with impregnated particles of a different material. Ceramic composites with dispersed metal inclusions fall into this category; 2) those that are made by stacking laminas. Each lamina consisting of a matrix (resin) may have short fibers dispersed in the matrix or may contain continuous fibers oriented in specific directions in the plane of each lamina. Laminas are glued and cured with the same epoxy as the matrix of the laminas. Upon curing, the interfaces between the laminas are assumed to have the same strength and properties as the matrix of the laminas. Regardless of which composite we consider; the composites are obviously non-homogeneous and non-isotropic. This is the major source of difficulty in deriving the mathematical models for composites that describe their deformation physics. In the following we present a brief literature review consisting of the development of mathematical models using energy methods and principle of virtual work and their numerical solutions primarily using finite element method.

The classical plate theory is based on Kirchhoff hypothesis (kinematic assumptions) for the deformation of the entire laminated plate cross section [1] [2]. Some analytical solutions of the laminated plate problems are represented based on series solutions in references [3] [4] [5] [6]. It is well known that due to pronounced shear deformation effects in composites, the theories based on Kirchhoff hypothesis do not work well. The laminated composite models incorporating transverse shear deformation effects either by using an assumed stress field [7] [8] [9] or by considering an assumed displacement field [10] [11] [12] have been reported. The extension of the displacement based mathematical models using higher order functions in thickness coordinate have been reported in references [13] - [23]. A summary of many of these works can be found in reference [24]. Since the advent of p-version of finite element method (earlier works can be found in references [24] - [35] ), its application for solving laminated composite problems grew. More recently the p-version hierarchical finite element method has been applied to plates and shells as well as beams including laminated composites in references [36] [37] [38] [39]. In these works, hierarchical p-approximation for displacements is used in the direction transverse to the plane of the laminate. We remark that the mathematical models currently used for composite mechanics are:

1) Mostly derived based on energy methods or principle of virtual work.

2) Only valid for reversible mechanical deformation.

3) Based on kinematic assumptions for the deformation of the cross section (except references [28] - [35] ).

In a recent paper Surana et al. [40] showed that all currently used beam mathematical models for isotropic and homogeneous matter are thermodynamically inconsistent. That is, these mathematical models can neither be derived using the conservation and balance laws of classical continuum mechanics (CCM) nor non-classical continuum mechanics (NCCM). The authors pointed out that the root cause of thermodynamic inconsistency is the use of a priori kinematic assumptions. The authors also presented a kinematic assumptions free formulation for bending of beams that is accurate for slender as well as deep beams and is free of shear correction factors that plague all currently used mathematical models. In another recent paper Surana et al. [41] also showed that all plate/shell mathematical models for isotropic homogeneous matter used currently that are based on kinematic assumption(s) and derived using energy methods also are thermodynamically inconsistent. The authors also presented a thermodynamically consistent plate/shell formulation with complete 3D deformation physics. This new formulation is free of kinematic assumption(s), is accurate for very thin as well as thick plate/shells and does not require any shear correction factors. Thus, we see that introduction of a priori kinematic assumptions followed by the derivation of the beam, plate/shell mathematical models (generally using energy methods) result in thermodynamically inconsistent mathematical models. We remark that in the work presented here, the objective is not to embark on developing better or more sophisticated methods for obtaining the solutions of the mathematical models for isotropic and non-isotropic matter or development of new approaches for deriving mathematical models based on energy methods. This work focuses on determining consistency and legitimacy of the mathematical models themselves for isotropic as well as non-isotropic matter that are derived based on energy methods and principle of virtual work for BVPs and IVPs as well as the solution methods based on these two approaches.

1.2. Scope of Present Work

Based on the principles of classical continuum mechanics (CCM) i.e. the conservation and balance laws, consistent derivation of constitutive theories based on the entropy inequality and elements of the calculus of variations, we establish some guidelines that are helpful in the investigations presented in this paper.

1) Mathematical models consisting of the conservation and balance laws of CCM and non-classical continuum mechanics (NCCM) with consistent constitutive theories based on the entropy inequality ensure thermodynamic equilibrium during the entire evolution. These mathematical models are meritorious, hence preferred if these can be derived for the deformation physics under consideration. This would be the preferred choice for the mathematical models for composites.

2) Calculus of variations considers extremum of functionals and provides a consistent and organized mathematical framework for obtaining solutions of the mathematical models of the initial value problems (IVPs) and boundary value problems (BVPs) regardless of how these models are derived and regardless of the nature of the differential operators in these models [42] [43].

3) We first consider elements of calculus of variations as applied to BVPs. We construct a functional I corresponding to the BVP or corresponding to the physics of the BVP. If I is differentiable in its arguments, then δI is unique and δI = 0 is a necessary condition for an extremum of I. If ${\delta }^{2}I$ provides a unique extremum principle (sufficient condition, $>0,=0,<0$ ), a solution obtained from δI = 0 is unique, hence yields a unique extremum of I. The solution obtained from δI = 0 also satisfies Euler’s equation, a differential or partial differential equation(s) extracted from δI = 0 using fundamental lemma or the fourth lemma. Thus, a correspondence between the extremums of functional and solution of BVP (Euler’s equation) is established. In structural mechanics for isotropic as well as non-isotropic matter one constructs an energy functional I consisting of strain energy and potential energy of loads. δI = 0 is a necessary condition for an extremum I. At this stage we could use δI = 0: 1) to construct a method of approximation such as finite element method to find an approximate solution that yields extremum of energy functional I and/or 2) we could derive Euler’s equation from δI = 0 using fundamental lemma or fourth lemma. This would be the differential form of the mathematical model associated with the energy functional I. In this paper we investigate whether this approach of using an energy functional consisting of strain energy and potential energy to derive the differential form of mathematical model for BVPs is consistent and valid way to derive differential mathematical models for homogeneous, isotropic as well as for non-homogeneous, non-isotropic matter. We also investigate the approach based on principle of virtual work.

4) Next we consider elements of calculus of variations as applied to IVPs. In this case we construct a space-time functional I corresponding to the IVP or corresponding to the physics of the IVP. If I is differentiable in its arguments, then δI is unique and δI = 0 is a necessary condition for an extremum of I. If ${\delta }^{2}I$ provides a unique extremum principle (sufficient condition, $>0,=0,<0$ ), then a solution obtained from δI = 0 also satisfies Euler’s equation, a partial differential equation in space and time extracted from δI = 0 using fundamental lemma or the fourth lemma. Thus, we have a correspondence between the extremums of space-time functionals and the solutions of the IVPs (Euler’s equation in space and time). In structural mechanics for isotropic, homogeneous as well as for non-isotropic, non-homogeneous matter one constructs energy functional I consisting of the kinetic energy, strain energy and potential energy of loads. δI = 0 is a necessary condition for an extremum of I. At this stage (as in the case of BVPs) we could use δI = 0: 1) to construct a method of approximation such as space-time finite element method to find an approximate solution that yields extremum of the energy functional and/or 2) we could derive Euler’s equation from δI = 0 using fundamental or fourth lemma. This would be the differential form of the mathematical model associated with the energy functional of I. In this paper we investigate whether this approach of deriving mathematical models for IVPs in differential form from the energy functional consisting of kinetic energy, strain energy and potential energy of loads yields a consistent and valid mathematical model in the differential form for homogeneous, isotropic as well as for non-homogeneous, non-isotropic matter. Similar investigation is also presented for principle of virtual work.

5) Investigations parallel to 3, and 4 are also presented using principle of virtual work.

6) Homogenization processes, their validity and their impact on energy methods and the mathematical models derived from energy methods are also equally important aspects to consider.

2. Consistency of the Methodologies of Deriving Mathematical Models and the Consistency of the Resulting Mathematical Models

We investigate the consistency and validity of following three different approaches of deriving mathematical models, as well as the consistency and validity of resulting mathematical models for initial value problems (IVPs) and boundary value problems (BVPs) for homogeneous, isotropic matter and non-isotropic, non-homogeneous matter.

1) Methodologies based on energy methods.

2) Methods based on conservation and balance laws of CCM and NCCM.

3) Methods based on principle of virtual work.

It can be easily shown that methods under 1) and 2) are a subset of calculus of variations, studies in applied mathematics for obtaining extremums of functionals as well as solutions of IVPs and BVPs [42] [43]. Thus, it is perhaps fitting to discuss methods 1) and 2) within the context of calculus of variations, a much broader area of study.

2.1. Calculus of Variations, Residual and Energy Functionals: BVPs

To be able to apply elements of calculus of variations to all BVPs rigorously and uniformly for obtaining their solutions or to be able to derive mathematical description of the associated boundary value problem, it is essential to mathematically classify all differential operators appearing in the BVPs into three possible categories: self-adjoint operators (linear and symmetric), non-self-adjoint operators (linear but not symmetric), and non-linear differential operators (neither linear nor symmetric). This mathematical classification allows us to treat all BVPs rigorously regardless of their origin or the field of application. We can use the elements of calculus of variations: 1) to obtain the solutions of BVPs if the mathematical model is known or 2) to derive the mathematical model if the energy functional associated with the BVP of interest is known. Traditionally the calculus of variations begins with

1) Existence of a functional I (generally by construction)

2) If I is differentiable in its arguments then δI is unique and δI = 0 is a necessary condition for an extremum of I.

3) ${\delta }^{2}I$ ( $>0,=0,<0$ ) must yield unique extremum principle or sufficient condition. This ensures that a solution obtained for δI = 0 yields a unique extremum of I.

A solution obtained from δI = 0 also satisfies Euler’s equation (PDE(s)) extracted from δI = 0 using fundamental lemma or the fourth lemma. We consider the boundary value problem (1) in which the differential operator A is linear.

$A\varphi -f=0\forall x\in {\Omega }_x$ (1)

with some boundary condition, ${\bar{\Omega }}_x={\Omega }_x\cup \Gamma$ is closure of ${\Omega }_x$, $\Gamma$ being closed boundary of open domain ${\Omega }_x$.

2.2. Residual Functional

The functional I can be the residual functional constructed using (1). If ${\varphi }_n$ is an approximation of $\varphi$ over ${\bar{\Omega }}_x$, then

$E=A{\varphi }_n-f\forall x\in {\bar{\Omega }}_x$ (2)

is the residual function and the residual functional I is constructed using E followed by δI and ${\delta }^{2}I$.

$I={\left(E\mathrm{,}E\right)}_{{\bar{\Omega }}_x}$ ; existence of I (by construction) (3)

$\delta I=2{\left(E,\delta E\right)}_{{\bar{\Omega }}_x}=0$ ; necessary condition for an extremum of I (4)

${\delta }^{2}I=2{\left(\delta E,\delta E\right)}_{{\bar{\Omega }}_x}>0$ ; when A is linear, unique extremum principle (5)

Thus, the integral form (4) is variationally consistent (VC) i.e., unconditionally stable. This process can used to obtain solutions of BVP (1). This of course requires that we know the mathematical model (1). It is obvious that in this case the Euler’s equation obtained δI = 0 (4) is indeed the BVP (1). A solution obtained from δI = 0 also satisfies Euler’s equation, hence is the solution of the BVP. One could show that this method based on residual functional can also be used for nonlinear differential operators [42] with unique extremum principle. Thus, this method based on residual functional for obtaining solutions of BVP is always unconditionally stable (variationally consistent, VC) regardless of the type of differential operator. This method requires that we know the mathematical model of the BVP, thus can always be used for isotropic, homogeneous matter for which the mathematical model can be derived using CCM or NCCM. In case of non-homogeneous and non-isotropic matter, we first need to derive the mathematical model before we can proceed with this residual functional approach described here. It is important to note that the residual functional is not energy functional. This approach of obtaining solutions of BVPs may not allow us to study specific physics of interest that may require concepts of stiffness, mass, damping, etc.

2.3. Energy Functional

The functional I can also be used to represent the energy functional consisting of strain energy and potential energy of loads. We consider self-adjoint, non-self-adjoint and non-linear operators A in (1). First let us consider the case when the differential operator A is known i.e. the mathematical model is known. Then, we can construct the integral form of (1) over ${\bar{\Omega }}_x$ using fundamental lemma of the calculus of variations.

${\left(A{\varphi }_n-f,v\right)}_{{\bar{\Omega }}_x}=0$ (6)

This integral form represents necessary condition from which the solution ${\varphi }_n$ is determined. v in (6) is the test function. Thus, at this stage (6) represents first variation of some functional I set to zero i.e. we assume that, there exists a functional I such that

$\delta I={\left(A{\varphi }_n-f,v\right)}_{{\bar{\Omega }}_x}=0$ (7)

${\varphi }_n$ and $v$ belong to appropriate scalar product space $V\subset H^{k\mathrm{,}p}\left({\bar{\Omega }}_x\right)$. The test function $v=0$ on boundary ${\Gamma }^{\mathrm{<mo>\; *\; </mo>}}$ if ${\varphi }_n={\varphi }_0$ (given or known) on ${\Gamma }^{\mathrm{<mo>\; *\; </mo>}}$. We have assumed existence of functional I and have constructed δI or integral form using fundamental lemma, but extremum principle still needs to be established. It is well known [42] that using (7), various methods of approximation can be considered.

1) Galerkin method (GM):

In this method we choose $v=\delta {\varphi }_n$ in (7), thus

$\delta I={\left(A{\varphi }_n-f,v\right)}_{{\bar{\Omega }}_x}=0;v=\delta {\varphi }_n$ (8)

or

${\left(A{\varphi }_n,v\right)}_{{\bar{\Omega }}_x}={\left(f,v\right)}_{{\bar{\Omega }}_x}$ (9)

or

$B\left({\varphi }_n\mathrm{,}v\right)=l\left(v\right)$ (10)

and

$\begin{array}{c}{\delta }^{2}I=\delta \left(B\left({\varphi }_n,v\right)-l\left(v\right)\right)=B\left(v,v\right)\\ =\left(Av,v\right)>0,=0,<0\forall v\in Visnotensured.\end{array}$ (11)

Hence, the integral form (11) in GM is variationally inconsistent (VIC) [42]. For this case we cannot construct functional I using functionals $B\left(\cdot ,\cdot \right)$, and $l(\cdot )$ in (14) or in any other way.

2) Petrov Galerkin Method (PGM) or weighted residual method (WRM):

In this method we also consider the integral form (10) but $v\ne \delta {\varphi }_n$ ; however we still choose v such that $v=0$ on ${\Gamma }^{\mathrm{<mo>\; *\; </mo>}}$ if $\varphi ={\varphi }_0$ on ${\Gamma }^{\mathrm{<mo>\; *\; </mo>}}$. Thus, in this method $v=\psi \ne \delta {\varphi }_n$, $\psi =0$ on ${\Gamma }^{\mathrm{<mo>\; *\; </mo>}}$ if $\varphi ={\varphi }_0$ on ${\Gamma }^{\mathrm{<mo>\; *\; </mo>}}$ and (10) becomes

$B\left({\varphi }_n\mathrm{,}\psi \right)=l\left(\psi \right)$ (12)

and

$\begin{array}{c}{\delta }^{2}I=\delta \left(B\left({\varphi }_n,\psi \right)-l\left(\psi \right)\right)=B\left(v,\psi \right)\\ =\left(Av,\psi \right)>0,=0,<0\forall \psi \in Visnotensured\text{.}\end{array}$ (13)

Thus, the integral form (13) is VIC. In this case we cannot construct functional I using $B\left(\cdot ,\cdot \right)$, and $l(\cdot )$ in (11) or any other way.

3) Galerkin method with weak form (GM/WF)

In 1) and 2) we have assumed the operator A to be linear. In this section we further assume that adjoint $A^{\mathrm{<mo>\; *\; </mo>}}$ of A is the same as A. The assumption of the operator A being linear and that $A^{\mathrm{<mo>\; *\; </mo>}}=A$ defines differential operators appearing in small deformation, small strain physics in structural mechanics in which the mechanical deformation is reversible. Such operators contain even order derivatives of the dependent variable with respect to spatial coordinates. We consider the following based on fundamental lemma.

${\left(A{\varphi }_n-f\mathrm{,}v\right)}_{{\bar{\Omega }}_x}=0;inwhichv=\delta {\varphi }_n$ (14)

or

${\left(A{\varphi }_n\mathrm{,}v\right)}_{{\bar{\Omega }}_x}={\left(f\mathrm{,}v\right)}_{{\bar{\Omega }}_x}$ (15)

We transfer half of the differentiation from ${\varphi }_n$ to v using integration by parts.

$B\left({\varphi }_n\mathrm{,}v\right)+\langle A{\varphi }_n\mathrm{,}v\rangle ={\left(f\mathrm{,}v\right)}_{{\bar{\Omega }}_x}$ (16)

in (16), $B\left({\varphi }_n\mathrm{,}v\right)$ contains all terms that have both ${\varphi }_n$ and v and $\langle A{\varphi }_n\mathrm{,}v\rangle$ is concomitant [42]. We can write (16) as

$B\left({\varphi }_n\mathrm{,}v\right)=l\left(v\right)$ (17)

$l\left(v\right)={\left(f\mathrm{,}v\right)}_{{\bar{\Omega }}_x}-\langle A{\varphi }_n\mathrm{,}v\rangle$ (18)

In Equation (17) $B\left({\varphi }_n\mathrm{,}v\right)$ is bilinear and symmetric i.e. $B\left({\varphi }_n\mathrm{,}v\right)=B\left(v,{\varphi }_n\right)$ and $l\left(v\right)$ is linear. For this case the functional I can be constructed.

1) $I\left({\varphi }_n\right)=\frac{1}{2}B\left({\varphi }_n\mathrm{,}{\varphi }_n\right)-l(\; v\; )$

2) $\delta I\left({\varphi }_n\right)=0$ gives $B\left({\varphi }_n,v\right)=l(\; v\; )$

3) ${\delta }^{2}I\left({\varphi }_n\right)=\delta \left(I\left({\varphi }_n\right)\right)=B\left(v,v\right)$

Since $B\left({\varphi }_n\mathrm{,}v\right)$ is symmetric $B\left(v,v\right)>0\forall v\in V$, hence, we have a unique extremum principle and the integral form (17) is variationally consistent (VC). $B\left(v,v\right)>0$ implies that $I\left({\varphi }_n\right)$ is minimized by a ${\varphi }_n$ obtained from the weak form (17).

Remarks

1) If the BVP (1) describes a structural mechanics problem with small deformation, small strain and reversible mechanical work, then $\frac{1}{2}B\left({\varphi }_n\mathrm{,}{\varphi }_n\right)$ in

$I\left({\varphi }_n\right)$ represents strain energy and $l\left({\varphi }_n\right)$ represents potential energy of loads.

2) For such problems, $I\left({\varphi }_n\right)$ can be called as energy functional. This in fact forms the basis for energy methods. One constructs $I\left({\varphi }_n\right)$ directly using strain energy and potential energy of loads (without the knowledge of mathematical model). Then, $\delta I\left({\varphi }_n\right)=0$ will yield (17), the weak form and the Euler’s equation(s) extracted from $\delta I\left({\varphi }_n\right)=0$ will be the corresponding mathematical model (1).

3) We need to ensure that construction of energy functional $I\left({\varphi }_n\right)$ using strain energy and potential energy of loads and the extraction of Euler’s equation (the mathematical model) holds for isotropic, homogeneous matter as well as for non-homogeneous and non-isotropic matter (such as composites).

4) We have shown that the energy methods employing energy functional can only be used for linear reversible mechanical deformation in which the mathematical models contain self-adjoint differential operators. This implies small deformation, small strain and reversible mechanical work deformation physics only. Thus, the usefulness of energy methods for designing computational processes (such as finite element method) or for establishing the mathematical model using the energy functional is limited to self-adjoint operators. We present illustrative examples in a subsequent section.

5) Based on (4), it is clear that the energy methods cannot be used to derive mathematical models for finite deformation and finite strain physics even for homogeneous isotropic matter as in this case, the differential operators are non-linear, hence the existence of energy functional cannot be proven or its construction is not possible using calculus of variations.

6) Use of energy functional approach for non-homogeneous and non-isotropic matter in deriving differential mathematical models needs to be investigated.

2.4. Calculus of Variation, Residual and Energy Functionals: IVPs

IVPs contain mathematical descriptions in which the dependent variables exhibit simultaneous dependence on space and time and the domain of definition of IVP is naturally a space-time domain. The space-time differential operator can be mathematically classified as: non-self-adjoint (linear but $A^{\mathrm{<mo>\; *\; </mo>}}\ne A$ ) and non-linear. When applying calculus of variation to obtain approximate solutions of IVP we naturally construct space-time functionals $I\left({\varphi }_n\left(x\mathrm{,}t\right)\right)=I\left({\varphi }_n\mathrm{,}t\right)$ and follow the same steps as in case of BVPs.

1) Existence of $I\left({\varphi }_n\mathrm{,}t\right)$ ; by construction.

2) If $I\left({\varphi }_n\mathrm{,}t\right)$ is differentiable in its arguments, then $\delta I\left({\varphi }_n\mathrm{,}t\right)$ is unique and $I\left({\varphi }_n,t\right)=0$ is a necessary condition for an extremum of $I\left({\varphi }_n\mathrm{,}t\right)$.

3) ${\delta }^{2}I\left({\varphi }_n,t\right)\left(>0,=0,<0\right)$ must yield a unique extremum principle or sufficient condition. This ensures that a ${\varphi }_n$ obtained from $\delta I\left({\varphi }_n,t\right)=0$ yields a unique extremum of $I\left({\varphi }_n\mathrm{,}t\right)$.

A ${\varphi }_n$ obtained from $\delta I\left({\varphi }_n,t\right)=0$ also satisfies Euler’s equation extracted from $\delta I\left({\varphi }_n,t\right)=0$ using fundamental lemma or fourth lemma of the calculus of variations.

2.5. Residual Functional

The functional $I\left({\varphi }_n\mathrm{,}t\right)$ can represent the residual functional constructed using IVP.

$A\varphi \left(x,t\right)-f\left(x,t\right)=0\forall x,t\in {\Omega }_{xt}={\Omega }_x\times {\Omega }_t$ (19)

with some boundary conditions and initial conditions. We consider the space-time differential operator to be linear (hence it is non-self-adjoint). If ${\varphi }_n$ is an approximation of $\varphi$ over ${\bar{\Omega }}_{xt}={\Omega }_{xt}\cup \Gamma$, $\Gamma$ being closed boundary of open domain ${\Omega }_{xt}$, then

$E=A{\varphi }_n-f\forall x\mathrm{,}t\in {\bar{\Omega }}_{xt}$ (20)

$I={\left(E,E\right)}_{{\bar{\Omega }}_{xt}}$ ; existence of I (by construction) (21)

$\delta I=2{\left(E,\delta E\right)}_{{\bar{\Omega }}_{xt}}=0$ ; necessary condition for an extremum of I (22)

${\delta }^{2}I=2{\left(\delta E,\delta E\right)}_{{\bar{\Omega }}_{xt}}>0$ ; unique extremum principle (23)

Thus, the space-time integral form (22) is space-time variationally consistent (STVF) [43]. This process can be used to obtain solutions of IVP (19). However, this requires that we know the mathematical model (19). It is straight forward to see that in this case the Euler’s equation obtained from δI = 0 (23) is indeed the IVP (19). STVC space-time integral form (23) ensures unconditionally stable computational process due to unique extremum principle (23). This method requires that we know the mathematical model (19) of the IVP, thus, can be used for isotropic, homogeneous matter for which the mathematical model (19) can be obtained from CCM or NCCM. In case of non-homogeneous, non-isotropic matter we first need to derive mathematical model before we can proceed with the residual functional approach presented here.

2.6. Energy Functional

The functional $I\left({\varphi }_n\mathrm{,}t\right)$ can also represent energy functional consisting of kinetic energy, strain energy and the potential energy of loads. This of course can be constructed directly over space time domain ${\bar{\Omega }}_{xt}$ without the knowledge of the mathematical model.

We recall that the space-time differential operators are either non-self-adjoint or non-linear. Thus, methods of approximation such as, space-time Galerkin method (STGM), STPGM, STGM/WF will yield space-time integral forms that are space-time variationally inconsistent i.e. if we consider (based on fundamental lemma)

$\left(A{\varphi }_n\left(x\mathrm{,}t\right)-f\left(x\mathrm{,}t\right)\mathrm{,}v\left(x\mathrm{,}t\right)\right)=0\forall x\mathrm{,}t\in {\Omega }_{xt}$ (24)

and if $B\left({\varphi }_n\mathrm{,}v\right)=l\left(v\right)$ in this final integral from any of the methods of approximation, then $B\left({\varphi }_n\mathrm{,}v\right)$ is bilinear but not symmetric. Hence, in this case the functional $I\left({\varphi }_n\mathrm{,}t\right)$ constructed using

$I\left({\varphi }_n\right)=\frac{1}{2}B\left({\varphi }_n,{\varphi }_n\right)-l\left({\varphi }_n\right)$ (25)

is not valid because

$\delta I\left({\varphi }_n\right)\ne B\left({\varphi }_n\mathrm{,}v\right)-l\left(v\right)$ (26)

Thus, for the space-time differential operators, if the mathematical model (19) is known, the energy functional $I\left({\varphi }_n\mathrm{,}t\right)$ cannot be constructed due to the fact that the space-time operator is either non-self-adjoint or non-linear. This raises an extremely important question and a serious concern whether the mathematical models derived using energy functional consisting of kinetic energy, strain energy and the potential energy of loads are even valid for isotropic, homogeneous matter. The validity of the energy functional itself and the resulting mathematical model for non-isotropic and non-homogeneous matter also need to be investigated.

3. CCM and NCCM: Mathematical Models for Homogeneous Isotropic and Non-Homogeneous, Non-Isotropic Continua (Composites)

The conservation and balance laws of classical continuum mechanics (CCM) in Lagrangian description: Conservation of mass (CM), balance of linear momenta (BLM) balance of angular momenta (BAM), first law of thermodynamics (FLT), second law of thermodynamics (SLT) provide valid and thermodynamically consistent mathematical models for deforming solid continua. The derivation of the conservation and balance laws in differential form assumes that the continuous matter is homogeneous and isotropic, hence each material point in a volume of matter is identical. This assumption implies that the volume of matter considered in the derivation of the conservation and balance laws is of no consequence. How many material points are present in the volume under consideration in the derivation of a balance law is of no consequence as each material point is identical, thus we may as well consider only one material point. This allows us to set the integrand to zero (at a material point). Thus, all balance laws in the differential form are the rate laws that are valid for a single material point due to the assumptions that the solid continua are homogeneous and isotropic. NCCM permits inclusion of additional physics in the derivation of the conservation and balance laws, but the assumption that the matter is homogeneous and isotropic must hold in order to be able to obtain the balance laws in the differential forms. Derivation of constitutive theory follows from entropy inequality, still based on the assumptions that the matter is isotropic and homogeneous.

Remarks

1) If we examine the conservation and balance laws in the absence of constitutive theory, we can conclude these are independent of the constitutive theories. The stress tensor $\sigma$ and heat vector $q$ are assumed to be present in the deforming continua as a consequence of kinematics of deformation and thermal disturbances applied to the volume of matter. This in fact is true.

2) As explained immediately preceding the remarks, that the assumption of homogeneity and isotropy are intrinsic requirements in deriving the differential form of the conservation and balance laws.

3) These assumptions in 2) must hold regardless of the constitutive theories.

4) Based on 1)-3), we have additional restrictions on constitutive relations. Regardless of whether we are to substitute the constitutive equations in balance laws or if we are to consider them as additional equations augmenting the conservation and balance laws, the constitutive theories must be for the isotropic and homogeneous matter as the conservation and balance laws only hold for isotropic and homogeneous continua. In other words, anisotropic or orthotropic constitutive relations can neither be included as additional equations in the conservation and balance laws of CCM or NCCM nor can these be substituted in them.

In view of the remarks, the phenomenological orthotropic constitutive theories such as in Section 3.1, though possibly verifiable experimentally, but cannot be used in the conservation and balance laws of CCM or NCCM. More specifically, in the stress strain relationship used for non-isotropic matter such as composites, the stresses as a function of strain cannot be substituted in the BLM and FLT (or considered as additional equations) as these balance laws are only valid for isotropic, homogeneous matter.

Thus, it is straightforward to conclude that the conservation and balance laws of CCM (or NCCM) cannot be used as a possible mathematical model for non-isotropic and non-homogeneous matter in conjunction with phenomenological constitutive theory.

At this stage, it is worth mentioning that in many works the differential form of BLM (valid for isotropic and homogeneous matter) is used in the process of determining solutions of model problems in composite mechanics using orthotropic constitutive theories. Based on the details presented here, this is not valid and the resulting solutions are not correct solutions of the actual model problems.

3.1. Material Properties: Non-Homogeneous, Non-Isotropic, Orthotropic (Laminated Composites)

The process of constructing laminas with preferred orientation of the fibers at the centerline of each lamina and the stacking of lamina in a desired sequence of material direction orientations by gluing them with epoxy (same resin as the matrix of the laminas) and the curing to construct a laminate is a well-known procedure. In this type of construction, each lamina contains fibers assumed to be at the center of the lamina with different material properties in two perpendicular directions in the plane of the lamina. Matrix (i.e. resin) is assumed to be isotropic and homogeneous with mechanical properties significantly different from the fibers. However, at the common boundary and in its vicinity between the glued laminas, the material is essentially resin (or matrix) which is assumed to be isotropic and homogeneous. Thus, our view is that the material properties make a transition from the center of the lamina (properties close to fiber) to the transverse boundaries of the lamina (properties close to resin or matrix).

In the published works based on rules of mixtures and volume fractions and other homogenization considerations, the mechanical properties of the fibers and the matrix are used to obtain equivalent homogenized mechanical properties that are called orthotropic in the two orthogonal directions in the plane of the lamina and the third one orthogonal to the lamina plane. These mechanical properties are assumed to hold throughout the lamina. When laminas with these homogenized properties are used to construct a laminate, naturally some discontinuity of stress is realized at the lamina interfaces if a single common displacement approximation is used for the entire laminate as the case is in classical plate theory (CPT), first-order shear deformation theory (FSDT), higher-order shear deformation theory (HSDT) etc. The homogenization process is somewhat empirical and phenomenological and the assumption of the same homogenized properties to be valid throughout the lamina thickness seems not in accordance with what the construction of the laminate suggests. From the point of view of the work presented here, it is sufficient for us to note that from the actual non-homogeneous, non-isotropic matter we obtain homogenized material properties that still remain non-isotropic.

3.2. Constitutive Equations for Orthotropic Materials (Laminated Composites)

CCM and NCCM only provide a thermodynamically consistent framework for deriving constitutive theories for homogeneous and isotropic continua, solid or fluent. For small deformation, small strain, isotropic and homogeneous solid continua with reversible mechanical deformation one could show [44] that the constitutive theory for Cauchy stress tensor $\sigma$ can be derived using

${\sigma }_{kl}={\rho }_0\frac{\partial \Phi \left(\epsilon \right)}{\partial {\epsilon }_{kl}}$ (27)

or

${\sigma }_{kl}={\rho }_0\frac{\partial \pi \left(\epsilon \right)}{\partial {\epsilon }_{kl}}$ (28)

in which $\Phi \left(\epsilon \right)$ is Helmholtz free energy density and $\pi \left(\epsilon \right)$ is strain energy density. The frame invariance requirement of the material coefficients in the constitutive theory necessitates that instead of $\Phi =\Phi \left(\epsilon \right)$ and $\pi =\pi \left(\epsilon \right)$ we must consider $\Phi$ and $\pi$ functions of the invariants of $\epsilon$, either principal invariants $I_{\epsilon }\mathrm{,}I{I}_{\epsilon }\mathrm{,}II{I}_{\epsilon }$ or $i_{\epsilon }\mathrm{,}i{i}_{\epsilon }\mathrm{,}ii{i}_{\epsilon }$ [44]. Since the two sets of invariants are related, we can choose either. Thus, we consider

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

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

Details of the derivation can be found in reference [44]. A linear constitutive theory resulting from (29) or (30) will contain only two material coefficients: Lame’s constants $\mu$ and $\lambda$ or modulus of elasticity E and Poisson’s ratio $\nu$. A constitutive theory with more than two material coefficients will be non-linear constitutive theory in $\epsilon$ and its invariants.

Currently used constitutive theories for anisotropic and orthotropic materials are derived by considering Taylor series expansion of $\Phi \left(\epsilon \right)$ or $\pi \left(\epsilon \right)$ in $\epsilon$ about a known configuration and then substituting it in (27) or (28). A linear constitutive theory resulting from this approach for Cauchy stress tensor $\sigma$ can be written as.

${\sigma }_{ij}=c_{ijkl}{\epsilon }_{kl}$ (31)

In (31) $\sigma$ and $\epsilon$, symmetric tensors of rank two are related through $c_{ijkl}$, a tensor of rank four containing 81 material coefficients for an anisotropic material. Symmetry of $\sigma$, $\epsilon$ and symmetry of $c_{ijkl}$ helps us in reducing the 81 material coefficients to 18. Further consideration of a monoclinic material with three orthogonal planes of symmetry reduces theses 18 coefficients to nine engineering coefficients: Young’s moduli $E_1$, $E_2$, $E_3$ ; Poisson’s ratios ${\nu }_{21}$, ${\nu }_{31}$, ${\nu }_{23}$ and shear moduli: $G_{23}$, $G_{31}$, $G_{12}$. In Poisson’s ratios, the first subscript refers to directions of stress and the second subscript is the direction of strain. Using Voigts notation we can write this constitutive theory as

$\left\{\sigma \right\}=\left[D\right]\left\{\epsilon \right\}$ (32)

or

$\left\{\epsilon \right\}=\left[C\right]\left\{\sigma \right\}\mathrm{<mo>\; ;\; </mo>}\left[D\right]={\left[C\right]}^{-1}$ (33)

in which

${\left\{\sigma \right\}}^{\text{T}}=\left[{\sigma }_{11}\mathrm{,}{\sigma }_{22}\mathrm{,}{\sigma }_{33}\mathrm{,}{\sigma }_{23}\mathrm{,}{\sigma }_{31}\mathrm{,}{\sigma }_{12}\right]$ (34)

${\left\{\epsilon \right\}}^{\text{T}}=\left[{\epsilon }_{11}\mathrm{,}{\epsilon }_{22}\mathrm{,}{\epsilon }_{33}\mathrm{,}{\gamma }_{23}\mathrm{,}{\gamma }_{31}\mathrm{,}{\gamma }_{12}\right]$ (35)

and

$\left[C\right]=\left[\begin{array}{cccccc}1/E_1& -{\nu }_{21}/E_2& -{\nu }_{31}/E_3& & & \\ & 1/E_2& -{\nu }_{32}/E_3& & \left[0\right]& \\ & & 1/E_3& & & \\ & & & 1/G_{23}& 0& 0\\ & \text{Symm}\text{.}& & & 1/G_{31}& 0\\ & & & & & 1/G_{12}\end{array}\right]$ (36)

in which ${\gamma }_{32}=2{\epsilon }_{32}$, ${\gamma }_{31}=2{\epsilon }_{31}$, and ${\gamma }_{12}=2{\epsilon }_{12}$. Generally, directions 1, 2, 3, refer to the lamina coordinate system. $\left[C\right]$ is transformed from the lamina coordinate system to the global x-frame [44].

Remarks

1) The constitutive theory (32) is purely phenomenological. It is calibrated by conducting experiments to determine nine material coefficients.

2) A serious drawback of assuming $\Phi =\Phi \left(\epsilon \right)$ and $\pi =\pi \left(\epsilon \right)$ is that material coefficients are now functions of $\epsilon$ in a known configuration, hence are not frame invariant, a strict requirement for the material coefficients. However, if the material coefficients are not dependent on deformation (i.e. are constant), then this is not a handicap, but certainly limits them to be constant throughout the evolution.

3) When the arguments of $\Phi$ and $\pi$ are invariant of $\epsilon$, the resulting linear constitutive theory for homogeneous, isotropic solid continua requires only two material coefficients. Clearly the derivation based on Taylor series expansion of $\Phi \left(\epsilon \right)$ and $\pi \left(\epsilon \right)$ in $\epsilon$ about a known configuration is not supported by continuum mechanics (CCM or NCCM). Thus, this constitutive theory for orthotropic matter is purely phenomenological or empirical.

4) We also remark that the orthotropic (or anisotropic in general) constitutive relations can neither be substituted in the balance laws nor the balance laws can be augmented using these constitutive relations as additional equations due to the fact that differential form of the balance laws only holds for homogeneous and isotropic matter.

4. Mathematical Models for BVPs Using Energy Functional

We consider methods of deriving mathematical models based on calculus of variations, more specifically energy method for BVPs for homogeneous, isotropic matter as well as for non-homogeneous, non-isotropic matter (such as composites). We recall that the orthotropic constitutive relations used for composites (in particular laminated composites) are also frequently used in other applications as well. As shown earlier, these constitutive relations are neither supported by CCM nor NCCM as these only consider isotropic homogeneous matter.

In this approach of deriving mathematical models for BVPs considered in this section, one constructs an energy functional I consisting of strain energy and potential energy of loads over the spatial domain ${\bar{\Omega }}_x$ of the deforming volume of matter. We consider small deformation and small strain and furthermore we consider the matter to be isotropic. For this case, we can define the energy functional I given in the following.

$I={\displaystyle \underset{{\bar{\Omega }}_x}{\int }}\left(-\text{tr}\left(\left[\sigma \right]\left[\epsilon \right]\right)+{\rho }_0{\left\{F^b\right\}}^{\text{T}}\left\{u\right\}\right)\text{d}\Omega$ (37)

in which $\left[\sigma \right]$ is the Cauchy stress tensor and $\left[\epsilon \right]$ is the linearized Green’s strain tensor. $\left\{F^b\right\}$ are body forces per unit mass and $\left\{u\right\}$ is a vector of displacements. In (37), $\left[\sigma \right]$ and $\left[\epsilon \right]$ are work conjugate pairs, their origin is from the first law of thermodynamics for isotropic, homogeneous matter. In CCM the second law of thermodynamics, entropy inequality can be written as (in Lagrangian description).

$\left({\rho }_0\frac{D\Phi }{Dt}+\eta \frac{D\theta }{Dt}\right)-\text{tr}\left(\left[\sigma \right]\left[\stackrel{\dot{}}{\epsilon }\right]\right)+\frac{q_i{g}_i}{\theta }\le 0$ (38)

in which $\Phi$ is the Helmholtz free energy density, $\eta$ is the entropy density, $\theta$ is temperature, $\left[\stackrel{\dot{}}{\epsilon }\right]$ is strain rate, $\left\{q\right\}$ is the heat flux and $\left\{g\right\}$ is a vector of temperature gradients. The entropy inequality (38) only holds for homogeneous, isotropic matter in which $\left[\sigma \right]$ and $\left[\stackrel{\dot{}}{\epsilon }\right]$ are rate of work conjugate pair. If the matter is non-isotropic and when the constitutive theory for $\sigma$ in terms of $\epsilon$ is like (32)-(36), then in the energy functional I (37) there is no basis for $\left[\sigma \right]$ and $\left[\stackrel{\dot{}}{\epsilon }\right]$ to be rate of work conjugate pair. Thus, use of (37) as energy functional for non-isotropic matter such as composites is not valid either. Hence, we can only proceed further with (37) if we assume homogeneous and isotropic matter. We do so in the following. We rewrite (37) using Einstein notation ( $\left[\sigma \right]$ and the $\left[\epsilon \right]$ are symmetric tensors of rank two).

$I={\displaystyle \underset{{\bar{\Omega }}_x}{\int }}\left(-{\sigma }_{ij}{\epsilon }_{ji}-{\rho }_0{F}_i^b{u}_i\right)\text{d}\Omega$ (39)

in which

${\epsilon }_{ji}=\frac{1}{2}\left(\frac{\partial u_j}{\partial x_i}+\frac{\partial u_i}{\partial x_j}\right)$ (40)

Thus

$\begin{array}{c}{\sigma }_{ij}{\epsilon }_{ji}=\frac{1}{2}{\sigma }_{ij}\left(\frac{\partial u_j}{\partial x_i}+\frac{\partial u_i}{\partial x_j}\right)=\frac{1}{2}\left({\sigma }_{ij}\frac{\partial u_j}{\partial x_i}+{\sigma }_{ij}\frac{\partial u_i}{\partial x_j}\right)\\ =\frac{1}{2}\left({\sigma }_{ji}\frac{\partial u_i}{\partial x_j}+{\sigma }_{ji}\frac{\partial u_j}{\partial x_i}\right)={\sigma }_{ji}\frac{\partial u_i}{\partial x_j}\end{array}$ (41)

using (41) in (39)

$I={\displaystyle \underset{{\bar{\Omega }}_x}{\int }}\left(-{\sigma }_{ji}\frac{\partial u_i}{\partial x_j}+{\rho }_0{F}_i^b{u}_i\right)\text{d}\Omega$ (42)

We transfer one order of differentiation from $u_i$ to ${\sigma }_{ji}$ in the first term of the integrand in (42).

$I={\displaystyle \underset{{\bar{\Omega }}_x}{\int }}\left(\frac{\partial {\sigma }_{ji}}{\partial x_j}{u}_i+{\rho }_0{F}_i^b{u}_i\right)\text{d}\Omega -{\displaystyle \underset{\Gamma }{\oint }{\sigma }_{ji}{n}_j{u}_i\text{d}\Gamma }$ (43)

If I is differentiable in its arguments, then $\delta I$ is unique and $\delta I=0$ is necessary condition for an extremum of I.

$\delta I={\displaystyle \underset{{\bar{\Omega }}_x}{\int }}\left(\frac{\partial {\sigma }_{ji}}{\partial x_j}\delta u_i+{\rho }_0{F}_i^b\delta u_i\right)\text{d}\Omega -{\displaystyle \underset{\Gamma }{\oint }{\sigma }_{ji}{n}_j\delta u_i\text{d}\Gamma }=0$ (44)

For simplicity we assume that $u_i=0$ on $\Gamma$ (BC) then $\delta u_i=0$ on $\Gamma$ and (44) reduces to

$\delta I={\displaystyle \underset{{\bar{\Omega }}_x}{\int }}\left(\frac{\partial {\sigma }_{ji}}{\partial x_j}+{\rho }_0{F}_i^b\right)\delta u_i\text{d}\Omega =0$ (45)

using fundamental lemma, we obtain the following Euler’s equation from (45)

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

We remark that since the volume ${\bar{\Omega }}_x$ contains homogeneous and isotropic matter, each material particle in ${\bar{\Omega }}_x$ is identical. Thus, size of the volume ${\bar{\Omega }}_x$ is irrelevant, hence permitting us to extract (46) from (45) at a material point. When the matter is non-isotropic and non-homogeneous, the first problem is with the strain energy term in (39) which is only valid for isotropic and homogeneous matter. Even if we could find a more justifiable expression for strain energy for non-isotropic matter, Euler’s equation similar to (46) cannot be extracted from δI = 0 due to the fact that volume $\bar{\Omega }$ contains non-isotropic matter, hence each material point could potentially be different in volume ${\bar{\Omega }}_x$. This forces us to maintain the volume integral form (45) instead of differential form (46) so that we know the specific material point in it. We note that (46) is balance of linear momenta for homogeneous, isotropic matter (for BVPs). Thus, when the deforming matter is isotropic and the mechanical deformation is reversible, the energy functional for BVPs only contains strain energy and potential energy of load in which case the Euler’s equations derived from the energy functional are indeed balance of linear momenta as expected. In this case, the differential operator (in displacements) is linear and its adjoint is the same as the operator itself. An integral from constructed using (46) based on fundamental lemma and GM/WF will yield a weak form

$\left(Au-f,v\right)=B\left(u,v\right)-l\left(v\right)=0$ (47)

in which $B\left(\cdot ,\cdot \right)$ is bilinear and symmetric and $l(\cdot )$ is linear and the energy functional I (Equation (39)) can be constructed using

$I=\frac{1}{2}B\left(u\mathrm{,}u\right)-l\left(u\right)$ (48)

Thus, in this case we see that the calculus of variations, energy methods and CCM yield the same mathematical model. This is only possible for a linear system in which the mechanical deformation is reversible. Strains and deformation are small and the matter is isotropic and homogeneous.

5. Mathematical Models for IVPs Using Energy Functional

The space-time differential operators appearing in IVPs are either non-self-adjoint or non-linear [43]. We can use STGM, STPGM, STWRM and STGM/WF for constructing space-time integral forms based on the fundamental lemma of the calculus of variations for space-time functionals (Surana et al. [43] ). The integral form resulting from any of the methods for both types of operators are space-time variationally inconsistent. That is, corresponding to these integral forms, representing first variation (δI) of some functional I, a unique extremum principle (a sufficient condition) does not exist. In this case, the resulting computations based on the integral form (δI = 0 for some functional I) are not ensured to be unconditionally stable. The bilinear functional (for linear space-time differential operators) in the space-time integral form is not symmetric. Hence, based on the calculus of variations [42] [43], it is not possible to construct a space-time functional using the functionals appearing in the space-time integral form. That is, corresponding to a valid IVP, a valid energy functional cannot be constructed. Conversely, there does not exist a valid energy functional such that the Euler’s equation extracted from its first variation set to zero represents a valid IVP. In the present work we show that if we construct an energy functional consisting of kinetic energy, strain energy and potential energy of loads, then the Euler’s equations extracted from its first variation of this functional set to zero cannot yield a valid mathematical model. This is obviously due to the fact that for a valid mathematical model for an IVP, it is not possible to construct a corresponding valid space-time functional due to the fact that for space-time linear differential operators $A^{\mathrm{<mo>\; *\; </mo>}}\ne A$. Thus, even if we know the PDEs describing IVP, then the functional I corresponding to these cannot be constructed. On the other hand, if we construct energy functional I directly without using space-time differential operator, the Euler’s equations resulting from the first variation of this functional cannot be the valid mathematical model (PDEs). First, we consider a simple IVP to illustrate the main source of problem in extracting mathematical models from energy functional for IVPs.

5.1. Model Problem

Consider the following IVP representing 1D balance of linear momenta in which all coefficients in the dimensionless form are unity. We have considered small deformation, small strain reversible mechanical deformation physics with linear constitutive theory for the Cauchy stress tensor and the matter is isotropic and homogeneous.

$\frac{{\partial }^{2}u}{\partial t^2}-\frac{{\partial }^{2}u}{\partial x_1^2}-F_1^b=0\forall x,t\in {\Omega }_{xt}={\Omega }_x\times {\Omega }_t=\left(0,L\right)\times \left(0,\tau \right)$ (49)

$\begin{array}{ll}\text{BCs}:\hfill & u\left(0,t\right)=u_0\hfill \\ \hfill & u\left(L,t\right)=u_L\hfill \end{array}(\begin{array}{c}\forall t\in \left(0,\tau \right)\end{array}$ (50)

$\begin{array}{ll}\text{ICs}:\hfill & u\left(x,0\right)=0\hfill \\ \hfill & \frac{\partial u}{\partial t}\left(x,0\right)=0\hfill \end{array}(\begin{array}{c}\forall x\in \left(0,L\right)\end{array}$ (51)

Figure 1 shows the space time domain ${\bar{\Omega }}_{xt}={\Omega }_{xt}\cup \Gamma$, $\Gamma ={\displaystyle \underset{i=1}{\overset{4}{\cup }}{\Gamma }_i}$ in which $\Gamma$ is the closed boundary of ${\bar{\Omega }}_{xt}$.

Let $n_x$ and $n_t$ be unit exterior normals to each ${\Gamma }_i$. Their values on ${\Gamma }_i$ ; $i=1,2,\cdots ,4$ are shown in Figure 1. We can write (49) as

$Au\left(x\mathrm{,}t\right)=f\left(x\mathrm{,}t\right)\forall x\mathrm{,}t\in {\Omega }_{xt}$ (52)

in which the space-time differential operator A is given by

$A=\frac{{\partial }^2}{\partial t^2}-\frac{{\partial }^2}{\partial x_1^2}\text{and}f\left(x\mathrm{,}t\right)=F_1^b$ (53)

Let $V\subset H^{k\mathrm{,}p}\left({\bar{\Omega }}_{xt}\right)$ be the scalar product space of functions admissible in (49).

Linearity of A:

Clearly operator A is linear as $A\left(\alpha u_1+\beta u_2\right)=\alpha Au+\beta A{u}_2$ holds for $\forall u_1\mathrm{,}{u}_2\in V$ and $\forall \alpha \mathrm{,}\beta \in R$.

Adjoint $A^{\mathrm{<mo>\; *\; </mo>}}$ of A:

Consider scalar product of $Au$ with v over ${\bar{\Omega }}_{xt}$ in which $v=\delta u$.

${\left(Au,v\right)}_{{\bar{\Omega }}_{xt}}={\left(\frac{{\partial }^{2}u}{\partial t^2}-\frac{{\partial }^{2}u}{\partial x_1^2},v\right)}_{{\bar{\Omega }}_{xt}}={\displaystyle \underset{{\bar{\Omega }}_{xt}}{\int }}\left(\frac{{\partial }^{2}u}{\partial t^2}v-\frac{{\partial }^{2}u}{\partial x_1^2}v\right)\text{d}x\text{d}t$ (54)

Figure 1. Space-time domain ${\bar{\Omega }}_{xt}$, boundaries ${\Gamma }_i$, BCs, ICs, and unit exterior normals to ${\Gamma }_i$.

we transfer all differentiation from u to v in both terms in the integral using integration by parts

$\begin{array}{c}{\left(Au,v\right)}_{{\bar{\Omega }}_{xt}}={\displaystyle \underset{{\bar{\Omega }}_{xt}}{\int }}\left(\frac{{\partial }^{2}u}{\partial t^2}v-\frac{{\partial }^{2}u}{\partial x_i^2}v\right)\text{d}\Omega +{\displaystyle \underset{\Gamma }{\oint }}\frac{\partial u}{\partial t}v{n}_t\text{d}\Gamma -{\displaystyle \underset{\Gamma }{\oint }}\frac{\partial v}{\partial t}u{n}_t\text{d}\Gamma \\ -{\displaystyle \underset{\Gamma }{\oint }}\frac{\partial u}{\partial x_1}v{n}_{x_1}\text{d}\Gamma +{\displaystyle \underset{\Gamma }{\oint }}\frac{\partial v}{\partial x_1}u{n}_{x_1}\text{d}\Gamma \end{array}$ (55)

We use BCs (50) and ICs (51) and the following

$\begin{array}{l}v\left(0,t\right)=0,v\left(L,t\right)=0\forall t\in \left[0,\tau \right]\\ v\left(x,0\right)=0,v\left(x,0\right)\forall x\in \left[0,L\right]\end{array}$ (56)

to simplify the integrals over closed boundary $\Gamma$ in (55) and obtain the following.

${\displaystyle \underset{\Gamma }{\oint }}\frac{\partial u}{\partial t}v{n}_t\text{d}\Gamma ={\displaystyle \underset{{\Gamma }_4}{\int }}{\frac{\partial u}{\partial t}|}_{\tau }v\left(n,\tau \right)n_t\text{d}\Gamma ;n_t=1;\text{d}\Gamma =\text{d}{x}_1:{\Gamma }_4:\left[0,L\right]$ (57)

${\displaystyle \underset{\Gamma }{\oint }}\frac{\partial v}{\partial t}u{n}_t\text{d}\Gamma ={\displaystyle \underset{{\Gamma }_4}{\int }}{\frac{\partial v}{\partial t}|}_{\tau }u\left(n,\tau \right)n_t\text{d}\Gamma ;n_t=1;\text{d}\Gamma =\text{d}{x}_1:{\Gamma }_4:\left[0,L\right]$ (58)

$\begin{array}{l}{\displaystyle \underset{\Gamma }{\oint }}\frac{\partial u}{\partial x_1}v{n}_{x_1}\text{d}\Gamma =0\\ {\displaystyle \underset{\Gamma }{\oint }}\frac{\partial v}{\partial x_1}u{n}_{x_1}\text{d}\Gamma ={\displaystyle \underset{{\Gamma }_1}{\int }}{\frac{\partial v}{\partial x_1}|}_{x_1=0}u\left(0,t\right)\left(-1\right)\text{d}\Gamma +{\displaystyle \underset{{\Gamma }_3}{\int }}{\frac{\partial v}{\partial x_1}|}_{x_1=L}u\left(L,t\right)\left(1\right)\text{d}\Gamma \end{array}$ (59)

$\text{d}\Gamma =\text{d}t\mathrm{<mo>\; ;\; </mo>}{\Gamma }_1\mathrm{<mo>\; :\; </mo>}\left[\mathrm{0,}\tau \right]\mathrm{<mo>\; ;\; </mo>}{\Gamma }_3\mathrm{<mo>\; :\; </mo>}\left[\mathrm{0,}\tau \right]$ (60)

We substitute (57)-(60) in (55) and rearrange terms

$\begin{array}{c}{\left(Au,v\right)}_{{\bar{\Omega }}_{xt}}={\displaystyle \underset{{\bar{\Omega }}_{xt}}{\int }}\left(\frac{{\partial }^{2}v}{\partial t^2}u-\frac{{\partial }^{2}v}{\partial x_1^2}u\right)\text{d}\Omega +{\displaystyle \underset{{\Gamma }_4}{\int }}{\frac{\partial u}{\partial t}|}_{\tau }v\left(x,\tau \right)\text{d}\Gamma \\ +{\displaystyle \underset{{\Gamma }_4}{\int }}{\frac{\partial v}{\partial t}|}_{\tau }v\left(x,\tau \right)\text{d}\Gamma +\langle Au,v\rangle \end{array}$ (61)

in which the concomitant is given by $\langle Au,v\rangle$ is given by [43]

$\langle Au,v\rangle ={\displaystyle \underset{{\Gamma }_1}{\int }}{\frac{\partial v}{\partial x_1}|}_{x_1=0}u\left(0,t\right)\left(-1\right)\text{d}\Gamma +{\displaystyle \underset{{\Gamma }_3}{\int }}{\frac{\partial v}{\partial x_1}|}_{x_1=L}u\left(L,t\right)\left(1\right)\text{d}\Gamma$ (62)

In (62) $\langle Au,v\rangle$ is the concomitant. We can write (61) as

$\begin{array}{c}{\left(Au,v\right)}_{{\bar{\Omega }}_{xt}}={\left(u,\frac{{\partial }^{2}v}{\partial t^2}-\frac{{\partial }^{2}v}{\partial x_1^2}\right)}_{{\bar{\Omega }}_{xt}}+{\left(v\left(x,\tau \right),{\frac{\partial u}{\partial t}|}_{\tau }\right)}_{{\Gamma }_4}\\ +{\left({\frac{\partial v}{\partial t}|}_{\tau },u\left(x,\tau \right)\right)}_{{\Gamma }_4}+\langle Au,v\rangle \end{array}$ (63)

Clearly the differential operator $A^{\mathrm{<mo>\; *\; </mo>}}$ (adjoint of A) from the right hand side of (63) is not the same as the differential operator A.

Remarks

1) Thus, in the IVP (for isotropic matter) the space-time differential operator is linear but its adjoint $A^{\mathrm{<mo>\; *\; </mo>}}$ is not the same as the operator A itself. This is primarily due to the additional terms appearing as a consequence of open boundary ${\Gamma }_4$ when performing integration by parts in time. On the open boundary

neither u nor $\frac{\partial u}{\partial t}$ are known, hence the terms resulting from integration by parts (in time) on ${\Gamma }_4$ must be considered in defining adjoint $A^{\mathrm{<mo>\; *\; </mo>}}.$

2) Since in this case A is linear but $A^{\mathrm{<mo>\; *\; </mo>}}\ne A$, the space-time integral form based on fundamental lemma using STGM/WF yields a bilinear functional (containing both u and the test function v) that is not symmetric. The consequence of this is that functional $B\left(\cdot \mathrm{,}\cdot \right)$ and $l(\cdot )$ in the integral cannot be used

to construct a functional $I\left(u\right)=\frac{1}{2}B\left(u\mathrm{,}u\right)-l\left(u\right)$ such that $\delta I\left(u\right)=0$ gives the integral from $B\left(u,v\right)=l\left(v\right)$.

3) Due to non-symmetry of functional $B\left(\cdot \mathrm{,}\cdot \right)$ in the integral form, an extremum principle ( ${\delta }^{2}I$ i.e. first variation of the integral form, $>0,=0,<0$ ) is not possible either.

4) We illustrate 1)-3) in the following.

5.2. Integral Form for Model Problem Using Fundamental Lemma and STGM/WF

Let $v=\delta u$, be the test function, then based on fundamental lemma we can write the following using (49) over the space-time domain ${\bar{\Omega }}_{xt}$

${\left(Au+F_1^b\mathrm{,}v\right)}_{{\bar{\Omega }}_{xt}}=0$ (64)

or

${\displaystyle \underset{{\bar{\Omega }}_{xt}}{\int }}\left(\frac{{\partial }^{2}u}{\partial t^2}v-\frac{{\partial }^{2}u}{\partial x_1^2}v+F_1^{b}v\right)\text{d}\Omega =0$ (65)

We transfer one order of differentiation from u to v for each term in the integrand of (65)

$\begin{array}{l}{\left(Au+F_1^b,v\right)}_{{\bar{\Omega }}_{xt}}\\ ={\displaystyle \underset{{\bar{\Omega }}_{xt}}{\int }}\left(-\frac{\partial v}{\partial t}\frac{\partial u}{\partial t}+\frac{\partial v}{\partial x_1}\frac{\partial u}{\partial x_1}\right)\text{d}\Omega +{\displaystyle \underset{\Gamma }{\oint }}v\frac{\partial u}{\partial t}\text{d}\Gamma -{\displaystyle \oint }v\frac{\partial u}{\partial x_1}{n}_{x_1}\text{d}\Gamma +{\displaystyle \underset{{\bar{\Omega }}_{xt}}{\int }}{F}_1^{b}v\text{d}\Gamma \\ =0\end{array}$ (66)

Simplification of the boundary integrals (see Equation (57) and (59)) using BCs, ICs and their variations gives

$\begin{array}{c}{\displaystyle \underset{\Gamma }{\oint }}\frac{\partial u}{\partial t}v{n}_t\text{d}\Gamma ={\displaystyle \underset{{\Gamma }_4}{\int }}{\frac{\partial u}{\partial t}|}_{\tau }v\left(n,\tau \right)n_t\text{d}\Gamma \\ ={\left(v\left(x,\tau \right),{\frac{\partial u}{\partial t}|}_{\tau }\right)}_{{\Gamma }_4};n_t=1,\text{d}\Gamma =\text{d}x;{\Gamma }_4:\left[0,L\right]\end{array}$ (67)

and

${\displaystyle \underset{\Gamma }{\oint }}\frac{\partial u}{\partial x_1}v{n}_{x_1}\text{d}\Gamma =0$ (68)

using (67) and (68) in (66)

$\begin{array}{l}{\left(Au+F_1^b,v\right)}_{{\bar{\Omega }}_{xt}}\\ ={\displaystyle \underset{{\bar{\Omega }}_{xt}}{\int }}\left(-\frac{\partial v}{\partial t}\frac{\partial u}{\partial t}+\frac{\partial v}{\partial x_1}\frac{\partial u}{\partial x_1}\right)+{\left(v\left(x,\tau \right),{\frac{\partial u}{\partial t}|}_{\tau }\right)}_{{\Gamma }_4}+{\displaystyle \underset{{\bar{\Omega }}_{xt}}{\int }}{F}_1^{b}v\text{d}\Gamma \\ =0\end{array}$ (69)

${\left(Au+F_1^b,v\right)}_{{\bar{\Omega }}_{xt}}=B\left(u,v\right)-l\left(v\right)=0$ (70)

in which

$B\left(u,v\right)={\displaystyle \underset{{\bar{\Omega }}_{xt}}{\int }}\left(-\frac{\partial v}{\partial t}\frac{\partial u}{\partial t}+\frac{\partial v}{\partial x_1}\frac{\partial u}{\partial x_1}\right)\text{d}\Omega +{\left(v\left(x,\tau \right),{\frac{\partial u}{\partial t}|}_{\tau }\right)}_{{\Gamma }_4}$ (71)

$l\left(v\right)=-{\displaystyle \int }{F}_1^{b}v\text{d}\Omega$ (72)

Remarks

1) $B\left(u\mathrm{,}v\right)$ is bilinear

2) $B\left(u\mathrm{,}v\right)\ne B\left(v\mathrm{,}u\right)$ i.e. $B\left(\cdot \mathrm{,}\cdot \right)$ is not symmetric

3) In this case if we construct a functional

$I=\frac{1}{2}B\left(u\mathrm{,}u\right)-l\left(u\right)$ (73)

$=\frac{1}{2}{\displaystyle \underset{{\bar{\Omega }}_{xt}}{\int }}\left(-{\left(\frac{\partial u}{\partial t}\right)}^2+{\left(\frac{\partial u}{\partial t}\right)}^2\right)\text{d}\Omega +\frac{1}{2}{\left(u\left(x,\tau \right),{\frac{\partial u}{\partial t}|}_{\tau }\right)}_{{\Gamma }_4}+{\left(F_1^b,u\right)}_{{\bar{\Omega }}_{xt}}$ (74)

Then

$\begin{array}{c}\delta I={\displaystyle \underset{{\bar{\Omega }}_{xt}}{\int }}\left(-\frac{\partial v}{\partial t}\frac{\partial u}{\partial t}+\frac{\partial v}{\partial x_1}\frac{\partial u}{\partial x_1}\right)\text{d}\Omega +\frac{1}{2}{\left(v\left(x,\tau \right),{\frac{\partial u}{\partial t}|}_{\tau }\right)}_{{\Gamma }_4}\\ +\frac{1}{2}{\left(u\left(x,\tau \right),{\frac{\partial v}{\partial t}|}_{\tau }\right)}_{{\Gamma }_4}+{\left(F_1^b,v\right)}_{{\bar{\Omega }}_{xt}}\end{array}$ (75)

which is not the same as (69) or (70)-(72).

4) Thus, in this case a functional I cannot be constructed using $B\left(\cdot \mathrm{,}\cdot \right)$ and $l(\cdot )$ from the integral form (weak form).

5.3. Derivation of Mathematical Model for IVP from the Energy Functional

In this section, we construct an energy functional I consisting of kinetic energy, strain energy and potential energy of loads that is believed to correspond to the initial value problem (49). We consider isotropic and homogeneous matter. Since (49) is in dimensionless form, we can construct an energy functional as follows

$I={\displaystyle \underset{{\bar{\Omega }}_{xt}}{\int }}\left(\frac{1}{2}{\left(\frac{\partial u}{\partial t}\right)}^2-{\sigma }_{11}{\epsilon }_{11}+F_1^{b}u\right)\text{d}\Omega$ (76)

${\epsilon }_{11}=\frac{\partial u}{\partial x_1}$ (77)

substituting for (77) in (76)

$I={\displaystyle \underset{{\bar{\Omega }}_{xt}}{\int }}\left(\frac{1}{2}{\left(\frac{\partial u}{\partial t}\right)}^2-{\sigma }_{11}\frac{\partial u}{\partial x_1}+F_1^{b}u\right)\text{d}\Omega$ (78)

We take first variation of I and set it to zero, necessary condition for an extremum of I.

$\delta I={\displaystyle \underset{{\bar{\Omega }}_{xt}}{\int }}\left(-\frac{\partial u}{\partial t}\frac{\partial \delta u}{\partial t}-{\sigma }_{11}\frac{\partial \delta u}{\partial x_1}+F_1^b\delta u\right)\text{d}\Omega =0$ (79)

We transfer differentiation from $\delta u$ to $\frac{\partial u}{\partial t}$ and ${\sigma }_{11}$ using integration by parts

$\delta I={\displaystyle \underset{{\bar{\Omega }}_{xt}}{\int }}\left(-\frac{{\partial }^{2}u}{\partial t^2}\delta u+\frac{\partial {\sigma }_{11}}{\partial x_1}\delta u+F_1^b\delta u\right)\text{d}\Omega +{\displaystyle \underset{\Gamma }{\oint }}\frac{\partial u}{\partial t}\delta u{n}_t\text{d}\Gamma -{\displaystyle \oint }{\sigma }_{11}{n}_x\delta u\text{d}\Gamma =0$ (80)

or

$\delta I={\displaystyle \underset{{\bar{\Omega }}_{xt}}{\int }}\left(-\frac{{\partial }^{2}u}{\partial t^2}+\frac{\partial {\sigma }_{11}}{\partial x_1}+F_1^b\right)\delta u+{\displaystyle \underset{\Gamma }{\oint }}\frac{\partial u}{\partial t}\delta u{n}_t\text{d}\Gamma -{\displaystyle \underset{\Gamma }{\oint }}{\sigma }_{11}{n}_x\delta u\text{d}\Gamma =0$ (81)

it has been shown

${\displaystyle \underset{\Gamma }{\oint }}\frac{\partial u}{\partial t}\delta u{n}_t\text{d}\Gamma ={\left({\delta u|}_{\tau }\mathrm{,}{\frac{\partial u}{\partial t}|}_{\tau }\right)}_{{\Gamma }_4}and{\displaystyle \underset{\Gamma }{\oint }}{\sigma }_{11}{n}_{x_1}\delta u\text{d}\Gamma =0$ (82)

using (82) in (81)

$\delta I={\displaystyle \underset{{\bar{\Omega }}_{xt}}{\int }}\left(-\frac{{\partial }^{2}u}{\partial t^2}+\frac{\partial {\sigma }_{11}}{\partial x_1}+F_1^b\right)\delta u+{\left({\delta u|}_{\tau },{\frac{\partial u}{\partial t}|}_{\tau }\right)}_{{\Gamma }_4}=0$ (83)

In the absence of the scalar product term over ${\Gamma }_4$, the Euler’s equation would have been (when the matter is homogeneous and isotropic)

$\frac{{\partial }^{2}u}{\partial t^2}-\frac{\partial {\sigma }_{11}}{\partial x_1}-F_1^b=0$ (84)

or

$\frac{{\partial }^{2}u}{\partial t^2}-\frac{{\partial }^{2}u}{\partial x_1^2}-F_1^b=0$ (85)

which is the same as the IVP. However, because of the open boundary ${\Gamma }_4$ the term ${\left(\mathrm{,}\right)}_{{\Gamma }_4}$ in (76) is not zero, hence (84) or (85) is not the Euler’s equation (mathematical model) corresponding to (83).

Remarks

1) We clearly see that there is no correspondence between the energy functional I (used in structural mechanics) and the initial value problem (85) even for isotropic and homogeneous matter. More specifically the mathematical models used currently that are believed to be based on energy methods for isotropic and homogeneous matter obviously cannot be derived using the energy functional consisting of kinetic energy, strain energy and the potential energy of load without neglecting the influence of the open boundary.

2) When deforming solid continua is non-isotropic and non-homogeneous, then

a) Strain energy terms have no basis (in ${ℝ}^3$ ) in the energy functional I.

b) Derivation of Euler’s equation is not possible due to the fact that ${\bar{\Omega }}_{xt}$ does not contain identical material particles, hence the volume integral must be maintained (83).

c) Even for isotropic, homogeneous matter, in the derivation of mathematical model (85) from (83) we must neglect the boundary integral over open boundary ${\Gamma }_4.$

5.4. Energy Functional in ${ℝ}^3$ and the Mathematical Model for IVP

The energy functional I consisting of kinetic energy, strain energy and potential energy of loads in $R^3$ can be written as (not nondimensionalized) follows for homogeneous, isotropic matter.

$I={\displaystyle \underset{{\bar{\Omega }}_{xt}}{\int }}\left(\frac{1}{2}{\rho }_0\left({\left(\frac{{\partial }^2{u}_1}{\partial t^2}\right)}^2+{\left(\frac{{\partial }^2{u}_2}{\partial t^2}\right)}^2+{\left(\frac{{\partial }^2{u}_3}{\partial t^2}\right)}^2\right)-{\sigma }_{ji}{\epsilon }_{ij}+{\rho }_0{F}_i^b{u}_i\right)\text{d}\Omega$ (86)

we note that (84) is only valid for the isotropic case, as for non-isotropic case ${\sigma }_{ji}{\epsilon }_{ij}$ as strain energy has no theoretical basis. Following Section 4, we can write (86) as

$I={\displaystyle \underset{{\bar{\Omega }}_{xt}}{\int }}\left(\frac{1}{2}{\rho }_0\left({\left(\frac{{\partial }^2{u}_1}{\partial t^2}\right)}^2+{\left(\frac{{\partial }^2{u}_2}{\partial t^2}\right)}^2+{\left(\frac{{\partial }^2{u}_3}{\partial t^2}\right)}^2\right)-{\sigma }_{ji}\frac{\partial u_i}{\partial x_i}+{\rho }_0{F}_i^b{u}_i\right)\text{d}\Omega$ (87)

Assuming that I is differentiable in its arguments, hence δI is unique, then $\delta I=0$ is a necessary condition for an extremum of I.

$\begin{array}{c}\delta I={\displaystyle \underset{{\bar{\Omega }}_{xt}}{\int }}\left({\rho }_0\left(\frac{\partial u_1}{\partial t}\frac{\partial \delta u_1}{\partial t}-\frac{\partial u_2}{\partial t}\frac{\partial \delta u_2}{\partial t}+\frac{\partial u_3}{\partial t}\frac{\partial \delta u_3}{\partial t}\right)-{\sigma }_{ji}\frac{\partial \delta u_i}{\partial x_j}+{\rho }_0{F}_i^b\delta u_i\right)\text{d}\Omega \\ =0\end{array}$ (88)

we transfer differentiation from $\delta u_1\mathrm{,}\delta u_2\mathrm{,}\delta u_3$ and $\delta u_i$ to their coefficients.

$\begin{array}{c}\delta I={\displaystyle \underset{{\bar{\Omega }}_{xt}}{\int }}\left({\rho }_0\left(-\frac{{\partial }^2{u}_1}{\partial t^2}\delta u_1-\frac{{\partial }^2{u}_2}{\partial t^2}\delta u_2-\frac{{\partial }^2{u}_3}{\partial t^2}\delta u_3\right)+\frac{\partial {\sigma }_{ji}}{\partial x_j}\delta u_i+{\rho }_0{F}_i^b\right)\delta u_i\text{d}\Omega \\ +{\displaystyle \underset{\Gamma }{\oint }}{\rho }_0\frac{\partial u_1}{\partial t}{n}_t\delta u_1\text{d}\Gamma +{\displaystyle \underset{\Gamma }{\oint }}{\rho }_0\frac{\partial u_2}{\partial t}{n}_t\delta u_2\text{d}\Gamma \\ +{\displaystyle \underset{\Gamma }{\oint }}{\rho }_0\frac{\partial u_3}{\partial t}{n}_t\delta u_3\text{d}\Gamma -{\displaystyle \underset{\Gamma }{\oint }}{\sigma }_{ji}{n}_{x_j}\delta u_i\text{d}\Gamma \end{array}$ (89)

in (89) even if the last boundary integral is zero (we assume so in the following), the other three boundary terms would yield integrals on open boundary at $t=\tau$ (say ${\Gamma }_4$ ). Thus, we can write (89) as

$\begin{array}{c}\delta I={\displaystyle \underset{{\bar{\Omega }}_{xt}}{\int }}\left({\rho }_0\left(-\frac{{\partial }^2{u}_1}{\partial t^2}\delta u_1-\frac{{\partial }^2{u}_2}{\partial t^2}\delta u_2-\frac{{\partial }^2{u}_3}{\partial t^2}\delta u_3\right)+\left(\frac{\partial {\sigma }_{ji}}{\partial x_j}\delta u_i+{\rho }_0{F}_i^b\delta u_i\right)\right)\text{d}\Omega \\ +{\left({\delta u_1|}_{\tau },{\frac{\partial u_1}{\partial t}|}_{n_t}\right)}_{{\Gamma }_4}+{\left({\delta u_2|}_{\tau },{\frac{\partial u_2}{\partial t}|}_{n_t}\right)}_{{\Gamma }_4}+{\left({\delta u_3|}_{\tau },{\frac{\partial u_3}{\partial t}|}_{n_t}\right)}_{{\Gamma }_4}\\ =0\end{array}$ (90)

The integrals over ${\Gamma }_4$ are not zero as they are on the open boundary ${\Gamma }_4$ (at $t=\tau$ ) on which $u_i$ or $\frac{\partial u_i}{\partial t}$ are not known. In the absense of integrals over ${\Gamma }_4$, the Euler’s equation resulting for $\delta I=0$ in (90) are

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

which is balance of linear momenta.

However, since the integrals over ${\Gamma }_4$ are not zero, Euler’s equations (91) cannot be derived from (90).

Remarks

1) In structural mechanics for homogeneous, isotropic matter as well as non-homogeneous, non-isotropic matter such as composites, the mathematical model (91) is used and is derived using the energy functional (86). The derivation presented here clearly shows this involves approximation of neglecting the boundary integrals over ${\Gamma }_4$ (open boundary at $t=\tau$ ) in (90).

2) Another important point to note is that, in the derivation of the Euler’s equations that they are valid at a material point and for every material point in the volume, there is the inherent assumption that the continua are homogeneous and isotropic. Only then the choice of volume is arbitrary as every material point in the volume is identically the same, hence we can set the integrand to zero which gives us the differential model at a material point.

3) We finally conclude that the derivation presented is only possible for isotropic matter and that (91) cannot be derived from the functional I in (86) as used commonly in the published works for both isotropic and non-isotropic solid continua.

6. Principle of Virtual Work

Principle of virtual work has been used in CCM as well as NCCM. The basic concepts remain the same in both cases. Thus, it suffices to consider CCM only in the following (in Langrangian description).

Definition: When a deformed body (volume of matter) is in stable equilibrium due to actions and reactions, then the work done by all forces due to virtual displacements is zero.

The virtual displacements at the material points are hypothetical or imaginary displacements such that due to the application of these, the boundary conditions and the loads as well as their points of application remain unaltered, hence ensuring that the equilibrium of the deformed body is not affected due to virtual displacements.

We consider a solid continuum with deformed volume $\bar{V}$ bounded by $\partial \bar{V}$. We consider small deformation, small strain, hence the undefined and defined coordinates as approximately same ( $x\approx \bar{x}$ ), thus no change in volume, implying that ( $V\approx \bar{V}$ ) and ( $\partial V\approx \partial \bar{V}$ ). Let the volume V be subjected to body forces,

inertial forces and tractions P on $\partial V$. Then, ${\rho }_0\frac{{\partial }^2{u}_i}{\partial t^2}\mathrm{,}{\rho }_0{F}_i^b$ and $\frac{\partial {\sigma }_{ji}}{\partial x_j};i=1,2,3$

are the inertial forces, body forces and the forces due to traction P on $\partial V$ (after using divergence theorem and Cauchy principle) per unit volume at a material point x or $\bar{x}$. Let $\delta u_i;i=1,2,3$ be virtual displacements at the material point x or $\bar{x}$ then the total virtual work for the volume V is (no sum over i, but sum over j) is given by the following and is zero.

${\displaystyle \underset{{\bar{\Omega }}_{xt}}{\int }}\left({\rho }_0\frac{{\partial }^2{u}_i}{\partial t^2}\delta u_i-{\rho }_0{F}_i^b\delta u_i-\frac{\partial {\sigma }_{ji}}{\partial x_j}\delta u_i\right)\text{d}\Omega =0;i=1,2,3$ (92)

where ${\bar{\Omega }}_{xt}={\bar{\Omega }}_x\times {\bar{\Omega }}_t={\bar{\Omega }}_x\times \left[0,\tau \right]$ ; ${\bar{\Omega }}_x$ ( ${\bar{\Omega }}_{xt}={\Omega }_{xt}\cup \Gamma$, $\Gamma$ being the closed boundary of ${\Omega }_{xt}$ ) being volume V and $\tau$ being final value of time.

or

${\displaystyle \underset{{\bar{\Omega }}_{xt}}{\int }}\left({\rho }_0\frac{{\partial }^2{u}_i}{\partial t^2}-{\rho }_0{F}_i^b-\frac{\partial {\sigma }_{ji}}{\partial x_j}\right)\delta u_i\text{d}\Omega =0$ (93)

Let $v_i=\delta u_i;i=1,2,3$ (94)

where $v_i$ are test functions, then we can write (93) as

${\displaystyle \underset{{\bar{\Omega }}_{xt}}{\int }}\left({\rho }_0\frac{{\partial }^2{u}_i}{\partial t^2}-{\rho }_0{F}_i^b-\frac{\partial {\sigma }_{ji}}{\partial x_j}\right)v_i\text{d}\Omega =0$ (95)

we note that $v_i=\delta u_i=0$ on ${\Gamma }^{\mathrm{<mo>\; *\; </mo>}}$ if $u_i=u_i^0$ (specified) on ${\Gamma }^{\mathrm{<mo>\; *\; </mo>}}$. We can also write (95) as

${\left({\rho }_0\frac{{\partial }^2{u}_i}{\partial t^2}-{\rho }_0{F}_i^b-\frac{\partial {\sigma }_{ji}}{\partial x_j}\mathrm{,}{v}_i\right)}_{{\bar{\Omega }}_{xt}}=0$ (96)

Let

$A_i\left(u\right)={\rho }_0\frac{{\partial }^2{u}_i}{\partial t^2}-\frac{\partial {\sigma }_{ji}}{\partial x_j}and{f}_i={\rho }_0{F}_i^b$ (97)

Then, (96) can be written as

${\left(A_i\left(u\right)-f_i\mathrm{,}{v}_i\right)}_{{\bar{\Omega }}_{xt}}=0$ (98)

using fundamental lemma of the calculus of variations we conclude from (98) that the Euler’s equations are

$A_i\left(u\right)-f_i=0\forall x,t\in {\Omega }_{xt};i=1,2,3$ (99)

which in fact are the balance of linear momenta.

Remarks

1) We see that starting with the statement of virtual work (92), we are able to obtain the corresponding differential form of the mathematical model using the fundamental lemma of the calculus of variations.

2) It is obvious that choice of the work terms in (92) is crucial. In this example we considered time dependent deformation of a simple solid for which the inertial forces, body forces and those due to tractions P on $\partial V$ are well understood and known. This may not be the case in more complex deformation physics.

3) Equivalence of the virtual work statement (92) and the integral form (98) resulting from it using the fundamental lemma is important to note. This confirms that (92) is in fact the integral statement based on fundamental lemma i.e. the integral form based on principle of virtual work is the same as the integral form resulting from the fundamental lemma when we know the mathematical model. It is necessary to know the mathematical model because it is only in this case that the statement of virtual work like (92) will contain all of the terms that are present in the mathematical model.

4) The conservation and balance laws of CCM and NCCM and the consistent constitutive theories are well recognized and accepted as valid mathematical descriptions of solid and fluent continua. In view of this and remark (3), the usefulness of principle of virtual work in deriving mathematical models is very minimal to none.

5) Once we have the mathematical model from CCM or NCCM, calculus of variations provides means for obtaining their solution (approximate), thus there is no need to resort to principle of virtual work as used currently in mechanics.

6) Principle of virtual work is the same as the integral statement resulting from fundamental lemma. Thus, the principle of virtual work is obviously not limited to solid mechanics. Any valid mathematical model based on CCM or NCCM regardless of the field of application or level of complexity can be used in conjunction with the fundamental lemma to construct integral form(s), that are indeed statements of virtual work (or energy). The converse is difficult to achieve. That is, a correct statement of virtual work without the knowledge of a mathematical model is quite difficult and in some cases, may not be possible due to lack of knowledge of precise nature of work or energy conjugate pairs.

7) From the published works, we note that principle of virtual work is primarily used for solid mechanics due to the fact that virtual work implies forces and virtual displacements. We have shown that principle of virtual work is not limited to solid mechanics. Once we know the mathematical model, fundamental lemma is all that is needed to construct integral form that is precisely that the corresponding statement of virtual work.

8) In view of the remarks (1)-(7), it is straight forward to conclude that the concept of principle of virtual work used presently is a more restricted view of the calculus of variations.

9) Mathematical classification of differential operators for BVPs and IVPs and the integral forms based on fundamental lemma for different choices of test function is a much more comprehensive and rigorous approach [42] [43] [45] as opposed to energy methods or principle of virtual work.

10) We only have transparent correspondence between energy methods, principle of virtual work and calculus of variations for self-adjoint differential operators in BVPs for which all three approaches yield the same energy functional, integral statement and mathematical model.

7. Summary and Conclusion

In the following we summarize the work presented in this paper and draw some conclusions.

1) The objective of this work is to investigate the consistency and legitimacy of the mathematical models for BVPs and IVPs for homogeneous and isotropic matter as well as for non-homogeneous and non-isotropic matter that are derived using energy functional approach. For BVPs, the energy functional consists of strain energy and the potential energy of loads. For IVPs, the energy functional consists of kinetic energy, strain energy and the potential energy of loads.

2) We consider BVPs in the following:

a) For linear, reversible mechanical deformation with small strain and small deformation, calculus of variations, energy methods and CCM, all yield the same mathematical model for isotropic, homogeneous matter. This of course is due to the fact that the differential operators in this case are linear and their adjoint is the same as the operator itself.

b) For non-homogeneous and non-isotropic matter with orthotropic or anisotropic constitutive relations, there is no basis for using ${\sigma }_{ji}{\epsilon }_{ij}$ as strain energy. We note that work conjugate pair ${\sigma }_{ji}{\epsilon }_{ij}$ is due to entropy inequality containing ${\sigma }_{ji}{\stackrel{\dot{}}{\epsilon }}_{ij}$ as rate of work conjugate pair. When the matter is not isotropic, there is no basis to assume ${\sigma }_{ji}{\epsilon }_{ij}$ as work conjugate pair. Even if we assume this to be true, the volume ${\bar{\Omega }}_x$ in this case contains non-isotropic matter, hence the differential model (Euler’s equation) cannot be derived as we must maintain the volume integral containing the different material points of non-isotropic matter. Thus, for non-isotropic matter regardless of whether the problem is linear or non-linear (but reversible mechanical deformation), a valid energy functional needs to be established using valid strain energy expression. Even if this is possible to accomplish, the derivation of the differential form of the mathematical model is not possible due to lack of isotropy.

3) In case of IVPs:

a) The space time operators resulting from CCM (BLM) are either non-self-adjoint or non-linear. For linear operators (as in structural mechanics) adjoint $A^{\mathrm{<mo>\; *\; </mo>}}$ of the space-time operator is not the same as the operator A. Thus, for IVPs use of the fundamental lemma and STGM/WF yields functionals $B\left(\cdot \mathrm{,}\cdot \right)$ and $l(\cdot )$ in which $B\left(\cdot \mathrm{,}\cdot \right)$ is bilinear but not symmetric and $l(\cdot )$ is linear. These functionals cannot be used to construct a energy functional

$I\left(u\right)=\frac{1}{2}B\left(u\mathrm{,}u\right)-l\left(v\right)$ such that $\delta I\left(u\right)=B\left(u\mathrm{,}v\right)-l\left(v\right)=0$ holds.

b) We have shown that using the energy functional (for isotropic matter) I consisting of kinetic energy, strain energy and potential energy of loads, it is not possible to derive a mathematical model (the space-time Euler’s equation) that represents BLM (in CCM). This is primarily due to the open boundary (at $t=\tau$ ) on which the boundary integral does not vanish. On the open boundary at $t=\tau$ the dependent variable and its time derivation are not known. Thus, even for isotropic, homogeneous reversible deformation physics, the differential form of the mathematical model cannot be derived from the energy functional consisting of kinetic energy, strain energy and the potential energy of loads.

c) When the deforming matter is non-isotropic, the first issue is the validity of ${\sigma }_{ji}{\epsilon }_{ij}$ as strain energy. The second issue is the same as in the case of isotropic matter i.e. presence of boundary integrals corresponding to the open boundary and third issue is that, in this case we cannot set the integrand to zero in the $\delta I=0$ due to the fact that the volume ${\bar{\Omega }}_x$ now contains different material points due to lack of isotropy.

d) Thus for IVPs, derivation of a valid differential form of mathematical model from the energy functional consisting of kinetic energy, strain energy and potential energy of loads is not possible for isotropic as well as non-isotropic continua (including composites). This is primarily due to the fact that the space-time differential operators are not self-adjoint (may be linear but $A^{\mathrm{<mo>\; *\; </mo>}}\ne A$ ).

4) Based on work presented in this paper, we can conclude that derivation of mathematical models from the energy methods only remains valid for isotropic, homogeneous matter and that a valid mathematical model using energy methods can only be derived for BVPs with linear, isotropic reversible mechanical work with small strain, small deformation physics.

5) We also remark that use of the energy functionals directly in processes such as finite element method employing integral forms over element domain ${\bar{\Omega }}_x^e$ also requires caution:

a) In using ${\sigma }_{ji}{\epsilon }_{ij}$ as strain energy when the matter is non-isotropic as there is no basis for it based on entropy inequality for isotropic matter.

b) In not discarding the terms on the open boundary, if integration by parts is employed in time.

c) To our knowledge, all currently used finite element processes are in violation of one or both of these aspects. When the matter is isotropic and homogeneous, violation of b) occurs and for non-homogeneous and non-isotropic matter, both a) and b) are violated.

6) It is shown that the integral form resulting from the principle of virtual work is exactly the same as the integral form from the fundamental lemma when the mathematical model is known. If the Euler’s equation(s) resulting from the principle of virtual work are not the same as the mathematical models based on CCM or NCCM, then the integral form in the principle of virtual work is flawed. Thus, in view of the thermodynamically consistent mathematical models based on CCM or NCCM and calculus of variations, there is virtually no benefit to pursue principle of virtual work anymore.

Acknowledgements

The first author is grateful for his endowed professorship and to the Department of Mechanical Engineering of the University of Kansas for financial support to the second author. Facilities provided by the Computational Mechanics Laboratory of Mechanical Engineering are also acknowledged. The first author gratefully acknowledges the help provided by his graduate student Sai S.C. Mathi in proof reading, in discussions and in preparation of final manuscript of the paper.

Conflicts of Interest

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

References

[1]	Reissner, E. and Stavsky, Y. (1961) Bending and Stretching of Certain Types of Heterogeneous Aeolotropic Elastic Plates. Journal of Applied Mechanics, 28, 402-408. https://doi.org/10.1115/1.3641719
[2]	Whitney, J.M. and Leissa, A.W. (1969) Analysis of Heterogeneous Anisotropic Plates. Journal of Applied Mechanics, 36, 261-266. https://doi.org/10.1115/1.3564618
[3]	Pagano, N.J. (1970) Exact Solutions for Rectangular Bidirectional Composites and Sandwich Plates. Journal of Composite Materials, 4, 20-34. https://doi.org/10.1177/002199837000400102
[4]	Pagano, N.J. (1970) Influence of Shear Coupling in Cylindrical. Bending of Anisotropic Laminates. Journal of Composite Materials, 4, 330-343. https://doi.org/10.1177/002199837000400305
[5]	Pagano, N.J. and Hatfield, S.J. (1972) Elastic Behavior of Multilayered Bidirectional Composites. AIAA Journal, 10, 931-933. https://doi.org/10.2514/3.50249
[6]	Srinivas, S. and Rao, A.K. (1970) Bending, Vibration and Buckling of Simply Supported Thick Orthotropic Rectangular Plates and Laminates. International Journal of Solids and Structures, 6, 1463-1481. https://doi.org/10.1016/0020-7683(70)90076-4
[7]	Reissner, E. (1944) On the Theory of Bending of Elastic Plates. Journal of Mathematics and Physics, 23, 184-191. https://doi.org/10.1002/sapm1944231184
[8]	Reissner, E. (1945) The Effect of Transverse Shear Deformation on the Bending of Elastic Plates. Journal of Applied Mechanics, 12, A68-A77.
[9]	Reissner, E. (1947) On Bending of Elastic Plates. Quarterly of Applied Mathematics, 5, 55-68. https://doi.org/10.1090/qam/20440
[10]	Hildebrand, F.B., Reissner, E. and Thomas, G.B. (1949) Notes on the Foundations of the Theory of Small Displacements of Orthotropic Shells. National Advisory Committee for Aeronautics, Washington DC.
[11]	Mindlin, R.D. (1951) Influence of Rotatory Inertia and Shear Flexural Motions of Isotropic Elastic Plates. Journal of Applied Mechanics, 18, 31-38.
[12]	Basset, A.B. (1890) On the Extension and Exure of Cylindrical and Spherical Thin Elastic Shells. The Royal Society, London.
[13]	Whitney, J.M. (1969) The Effect of Transverse Shear Deformation on the Bending of Laminated Plates. Journal of Composite Materials, 3, 534-547. https://doi.org/10.1177/002199836900300316
[14]	Whitney, J.M. and Pagano, N.J. (1970) Shear Deformation in Heterogeneous Anisotropic Plates. Journal of Applied Mechanics, 37, 1031-1036. https://doi.org/10.1115/1.3408654
[15]	Yang, P.C., Norris, C.H. and Stavsky, Y. (1966) Elastic Wave Propagation in Heterogeneous Plates. International Journal of Solids and Structures, 2, 665-684. https://doi.org/10.1016/0020-7683(66)90045-X
[16]	Mau, S.T. (1973) A Refined Laminated Plate Theory. Journal of Applied Mechanics, 40, 605-607. https://doi.org/10.1115/1.3423032
[17]	Nelson, R.B. and Lorch, D.R. (1974) A Refined Theory for Laminated Orthotropic Plates. Journal of Applied Mechanics, 41, 177-183. https://doi.org/10.1115/1.3423219
[18]	Lo, K.H., Christensen, R.M. and Wu, E.M. (1977) A High-Order Theory of Plate Deformation Part 1: Homogeneous Plates. Journal of Applied Mechanics, 44, 663-668.
[19]	Lo, K.H., Christensen, R.M. and Wu, E.M. (1977) A High-Order Theory of Plate Deformation Part 2: Laminated Plates. Journal of Applied Mechanics, 44, 560-575.
[20]	Levinson, M. (1980) An Accurate, Simple Theory of the Statics and Dynamics of Elastic Plates. Mechanics Research Communications, 7, 343-350. https://doi.org/10.1016/0093-6413(80)90049-X
[21]	Murthy, M.V.V. (1981) An Improved Transverse Shear Deformation Theory for Laminated Anisotropic Plates. Langley Research Center, National Aeronautics and Space Administration, Scientific and Technical Information Branch.
[22]	Reddy, J.N. (1984) A Simple Higher-Order Theory for Laminated Composite Plates. Journal of Applied Mechanics, 51, 745-752. https://doi.org/10.1115/1.3167719
[23]	Reddy, J.N. (1984) A Refined Nonlinear Theory of Plates with Transverse Shear Deformation. International Journal of Solids and Structures, 20, 881-896. https://doi.org/10.1016/0020-7683(84)90056-8
[24]	Szab, B.A. and Sahrmann, G.J. (1988) Hierarchic Plate and Shell Models Based on Pextension. International Journal for Numerical Methods in Engineering, 26, 1855-1881. https://doi.org/10.1002/nme.1620260812
[25]	Brombolich, L.J. and Gould, P.L. (1969) Finite Element Analysis of Shells of Revolution by Minimization of the Potential Energy Functional. Proceedings for the Conference on Applications of the Finite Element Method in Civil Engineering, Nashville, 279-307.
[26]	Szabo, B.A., Basu, P.K. and Rossow, M.P. (1977) Theoretical Manual and Users’ Guide for COMET-X. U.S. Department of Transportation, Federal Railroad Administration, Washington DC.
[27]	Lamprecht, R.M. and Basu, P.K. (1977) Some Trends in Computerized Stress Analysis. Proceedings of the Seventh Conference on Electronic Computation, St. Louis, August 1979.
[28]	Kalim, P. and Surana, K.S. (1985) Finite Element Formulation for Axisymmetric Shell Heat Conduction with Temperature and Gradients. The Winter Annual Meeting of ASME, Miami Beach, 17-21.
[29]	Surana, K.S. and Kalim, P. (1986) Transition Finite Elements with Temperature and Temperature Gradients as Primary Variables for Axisymmetric Heat Conduction. Computers & Structures, 24, 197-212. https://doi.org/10.1016/0045-7949(86)90279-8
[30]	Surana, K.S. and Phillips, R.K. (1987) Three Dimensional Curved Shell Finite Elements for Heat Conduction. Computers & Structures, 25, 775-785. https://doi.org/10.1016/0045-7949(87)90169-6
[31]	Surana, K.S. and Phillips, R.K. (1990) Three Dimensional Solid-Shell Transition Finite Elements for Heat Conduction. Computers & Structures, 26, 941-950. https://doi.org/10.1016/0045-7949(87)90111-8
[32]	Surana, K.S. and Sorem, R.M. (1989) Higher Order Hierarchical Plates and Curved Shell Elements Based on p-Approximation in the Thickness Direction. 7th International Conference of Mathematical Computational Modeling, Chicago, 849-854. https://doi.org/10.1016/0895-7177(90)90302-4
[33]	Surana, K.S. and Sorem, R.M. (1990) Curved Shell Elements for Elastostatics with p-Version in the Thickness Direction. Computers & Structures, 36, 701-719. https://doi.org/10.1016/0045-7949(90)90085-G
[34]	Surana, K.S. and Sorem, R.M. (1985) p-Approximation Curved Shell Elements for Elastostatics Analysis of Laminated Composite Plates and Shells. ASME Pressure Vessel and Piping Conference, Honolulu.
[35]	Surana, K.S. and Sorem, R.M. (1990) Curved Shell Elements Based on Hierarchical P-Approximation in the Thickness Direction for Linear Static Analysis of Laminated Composites. Internation Journal for Numerical Methods in Engineering, 29, 1393-1420. https://doi.org/10.1002/nme.1620290703
[36]	Surana, K.S. and Sorem, R.M. (1991) Completely Hierarchical p-Version Curved Shell Element for Laminated Composite Plates and Shells. Computational Mechanics, 7, 237-251. https://doi.org/10.1007/BF00370038
[37]	Surana, K.S. and Sorem, R.M. (1990) Higher Order Hierarchical Plates and Curved Shell Elements Based on p-Approximation in the Thickness Direction. Mathematical and Computer Modelling, 14, 849-854. https://doi.org/10.1016/0895-7177(90)90302-4
[38]	Surana, K.S. and Nguyen, S.H. (1990) Higher-Order Shear-Deformable Two Dimensional Hierarchical Beam Elements for Laminated Composites. Mathematical and Computer Modelling, 14, 893-898. https://doi.org/10.1016/0895-7177(90)90310-J
[39]	Surana, K.S. and Nguyen, S.H. (1990) Two-Dimensional Curved Beam Element with Higher-Order Hierarchical Transverse Approximation for Laminated Composites. Computers & Structures, 36, 499-511. https://doi.org/10.1016/0045-7949(90)90284-9
[40]	Surana, K.S., Mysore, D. and Reddy, J.N. (2019) Thermodynamic Consistency of Beam Theories in the Context of Classical and Non-Classical Continuum Mechanics and a Thermodynamically Consistent New Formulation. Continuum Mechanics and Thermodynamics, 31, 1283-1312. https://doi.org/10.1007/s00161-019-00744-8
[41]	Surana, K.S. and Mathi, S.S.C. (2020) Thermodynamic Consistency of Plate and Shell Mathematical Models in the Context of Classical and Non-Classical Continuum Mechanics and a Thermodynamically Consistent New Thermoelastic Formulation. American Journal of Computational Mathematics, 10, 167-220. https://doi.org/10.4236/ajcm.2020.102010
[42]	Surana, K.S. and Reddy, J.N. (2016) The Finite Element Method for Boundary Value Problems: Mathematics and Computations. CRC Press, Boca Raton. https://doi.org/10.1201/9781315365718
[43]	Surana, K.S. and Reddy, J.N. (2017) The Finite Element Method for Initial Value Problems. CRC/Taylor and Francis, Boca Raton. https://doi.org/10.1201/b22512
[44]	Surana, K.S. (2016) Advanced Mechanics of Continua. Applied and Computational Mechanics Series. CRC Press, Boca Raton. https://doi.org/10.1201/b17959
[45]	Gelfand, I.M., Fomin, S.V. and Silverman, R.A. (2000) Calculus of Variations. Dover Books on Mathematics. Dover Publications, Mineola.