Enhanced Constitutive Theories for Classical Thermoviscoelastic Polymeric Fluids

Karan S. Surana^*, Thomas J. Ezell

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

DOI: 10.4236/am.2022.1311058                

Abstract

This paper presents ordered rate nonlinear constitutive theories for thermoviscoelastic fluids based on Classical Continuum Mechanics (CCM). We refer to these fluids as classical thermoviscoelastic polymeric fluids. The conservation and balance laws of CCM constitute the core of the mathematical model. Constitutive theories for the Cauchy stress tensor are derived using the conjugate pair in the entropy inequality, additional desired physics, and the representation theorem. The constitutive theories for the Cauchy stress tensor consider convected time derivatives of Green’s strain tensor or the Almansi strain tensor up to order n and the convected time derivatives of the Cauchy stress tensor up to order m. The resulting constitutive theories of order (m, n) are based on integrity and are valid for dilute as well as dense polymeric, compressible, and incompressible fluids with variable material coefficients. It is shown that Maxwell, Oldroyd-B, and Giesekus constitutive models can be described by a single constitutive theory. It is well established that the currently used Maxwell and Oldroyd-B models predict zero normal stress perpendicular to the flow direction. It is shown that this deficiency is a consequence of not retaining certain generators and invariants from the integrity (complete basis) in the constitutive theory and can be corrected by including additional generators and invariants in the constitutive theory. Similar improvements are also suggested for the Giesekus constitutive model. Model problem studies are presented for BVPs consisting of fully developed flow between parallel plates and lid-driven cavities utilizing the new constitutive theories for Maxwell, Oldroyd-B, and Giesekus fluids. Results are compared with those obtained from using currently used constitutive theories for the three polymeric fluids.

Keywords

Polymeric Fluids, Maxwell, Oldroyd-B, Giesekus, Mathematical Models, Nonlinear Rate Constitutive Theories, Representation Theorem, Integrity, Generators, Invariants

1. Literature Review and Scope of Work

Constitutive theories for polymers date back to Maxwell [1], and Oldroyd-B [2] based primarily on phenomenological approaches. The observed motion (under microscope) of polymeric fluids is cast into appropriate mathematical models using empirical, theoretical, or experimental approaches. An organized treatment of these approaches can be found in Bird et al. [3] [4]. These mathematical models have enjoyed large success in predicting the reasonable behavior of polymeric fluids. In the last six decades, groundbreaking work in continuum mechanics based on a sound thermodynamic foundation has provided an incentive to examine these constitutive models more closely to determine if the constitutive theories for polymeric fluids can be derived strictly using principles of CCM, axioms of constitutive theories, the entropy inequality, and the representation theorem. Such a framework based on sound principles, when possible, has provided an incomparable mechanism for further enhancement of the constitutive theories for more complex physics. Surana et al. and Surana [5] [6] [7] used Conservation and Balance Laws (CBLs) of CCM and the conjugate pairs in the entropy inequality to determine constitutive tensors and their argument tensors augmented for additional physics to derive the constitutive theory for the deviatoric Cauchy stress tensor using the representation theorem. These theories utilized convected time derivatives of Green’s strain tensor (or the Almansi strain tensor) up to order n and the convected time derivatives of the deviatoric Cauchy stress tensor up to order m. Thus, these constitutive theories are referred to as ordered rate constitutive theories of orders m and n. The authors showed that: 1) the Maxwell model is a linear constitutive model corresponding to $m=1$ and $n=1$, 2) the Oldroyd-B model is a quasilinear or nonlinear simplified constitutive model that utilizes $m=1$ and $n=2$, and 3) the Giesekus constitutive model is a nonlinear constitutive model based on $m=1$ and $n=1$. These constitutive models are more simplified linear or nonlinear forms of the general constitutive theory of orders m and n that are obtained by discarding the generators and invariants from the list of combined generators and invariants based on the argument tensors of the deviatoric Cauchy stress tensor.

There are numerous other constitutive theories available in the published works that are mostly phenomenological and are derived based on 1D springs and 1D dashpots in series and/or parallel configurations. Such constitutive theories may be useful for the application in hand, but are difficult to extend for more complex applications due to a lack of sound thermodynamic and mathematical foundation. The incentive in references [5] [6] [7] [8] [9] was to show that currently used mathematical models for polymeric fluids such as Maxwell, Oldroyd-B and Giesekus have a thermodynamic foundation based on CCM.

It has been numerically verified that Maxwell and Oldroyd-B models predict zero normal stress in the direction perpendicular to the flow. For example, in the case of fully developed flow between parallel plates, these two models predict reasonable normal stress in the direction of flow, but the normal stress perpendicular to the direction of flow from these models is zero. We can view these models as simplified forms of constitutive theories based on rates m and n using integrity (complete basis). Then it is rather obvious that this problem of zero normal stress perpendicular to the direction of flow is a consequence of neglecting generators and/or invariants in the rate model of orders m, n which should have been retained. This forms the basis for enhancing the current models and improving their performance. The choice of which additional generators and/or invariants (over and beyond those that are already present in the current models) should be included is not trivial. Details are discussed in the following sections. This paper is organized into the sections:

(I) Conservation and balance laws (without derivation);

(II) Derivation of ordered rate constitutive theories of orders m, n;

(III) Simplification of the constitutive theory of order m, n into a single constitutive theory for dilute and dense polymeric fluids;

(IV) Obtaining currently used mathematical models using the constitutive theory in (III);

(V) Enhancement of the constitutive theory in (III) for improved prediction of normal stress in the direction normal to the flow;

(VI) Model problem studies using fully developed flow between parallel plates and fully developed flow in a lid-driven cavity with comparisons to the currently used models;

(VII) Summary and conclusions.

2. Mathematical Model

Conservation and Balance Laws (CBLs) of Classical Continuum Mechanics (CCM) in Eulerian description are used in fluid mechanics, hence for polymeric fluids as well. The CBLs in Eulerian description can be expressed purely in terms of velocities, hence by choosing velocities as observable quantities, displacements of material points can be ignored. Thus, displacements of material points are neither present in the CBLs nor can they be obtained using velocities as the velocities are observable quantities.

In Eulerian descriptions, stress measures must be considered using the deformed tetrahedron in the current configuration, hence, they can be contravariant or covariant, and the corresponding strain rate measures must be covariant or contravariant, respectively. This situation is more complex than in the case of solids. In this paper, we consider notation used in references [8] [9], i.e. quantities with overbar imply Eulerian description or their values in the current configuration. Thus, $x_i$ and ${\bar{x}}_i$ are coordinates of a material point in the reference and current configurations, respectively. Other than this, we use standard Einstein notation, index notation, or matrix and vector notation (whichever maintains more clarity of presentation). Quantities with subscripts are covariant measures and those with superscripts are contravariant measures. Since in this paper we consider convected time derivatives, we introduce the following notation [8] [9], if $\bar{\sigma }$ is the Cauchy stress tensor in the current configuration, then instead of $\bar{\sigma }$ we write ${\bar{\sigma }}^{\left(0\right)}$ or ${\bar{\sigma }}_{\left(0\right)}$. ${\bar{\sigma }}^{\left(0\right)}$ is the contravariant Cauchy stress tensor, in which the superscript (0) refers to a convected time derivative of order zero, i.e. the tensor itself. Thus, we have

${\bar{\sigma }}^{\left(i\right)}\mathrm{<mo>\; ;\; </mo>}i=\mathrm{0,1,}\cdots \mathrm{,}m$ (1)

${\bar{\sigma }}_{\left(j\right)}\mathrm{<mo>\; ;\; </mo>}i=\mathrm{0,1,}\cdots \mathrm{,}m$ (2)

In Equations (1) and (2), we have convected time derivatives of Cauchy stress tensor up to orders m. Equation (1) is a contravariant measure while Equation (2) is a covariant measure. Parentheses imply convected time derivatives. Square brackets are used to denote material derivatives. Thus, ${\epsilon }_{\left[0\right]}$ is a material derivative of order zero of the Green’s strain tensor. Let ${\gamma }_{\left(i\right)}\mathrm{<mo>\; ;\; </mo>}i=\mathrm{1,2,}\cdots \mathrm{,}n$ and ${\gamma }^{\left(j\right)}\mathrm{<mo>\; ;\; </mo>}j=\mathrm{1,2,}\cdots \mathrm{,}n$ be convected time derivatives of Green’s and Almansi strain tensors ${\epsilon }_{\left[0\right]}$ and ${\bar{\epsilon }}^{\left[0\right]}$ up to orders n, then the following constitute the rate of work conjugate pairs.

$\begin{array}{l}{\bar{\sigma }}^{\left(i\right)};i=0,1,\cdots ,m\\ {\gamma }_{\left(j\right)};j=1,2,\cdots ,n\end{array}$ (3)

$\begin{array}{l}{\bar{\sigma }}_{\left(i\right)};i=0,1,\cdots ,m\\ {\gamma }^{\left(j\right)};j=1,2,\cdots ,n\end{array}$ (4)

We can also consider Jaumann measures [3] [4] [8] [9], but we avoid these here due to their limitations [8] [9]. Thus, the CBLs and constitutive theories for polymeric fluids can be derived using measures in Equations (3) or (4). Such a derivation would be basis dependent as the measures in Equations (3) and (4) are basis dependent. To make the derivation basis independent, instead of using Equations (3) or (4), we consider the following notations for the convected time derivatives of the Cauchy stress tensor and strain tensor in the derivation of the CBLs and the constitutive theories.

$\begin{array}{l}{}^{\left(i\right)}\bar{\sigma }\mathrm{<mo>\; ;\; </mo>}i=\mathrm{0,1,}\cdots \mathrm{,}m\\ {}^{\left(j\right)}\gamma \mathrm{<mo>\; ;\; </mo>}j=\mathrm{1,2,}\cdots \mathrm{,}n\end{array}$ (5)

The conjugate pairs in Equation (5) can be chosen based on either Equation (3) or Equation (4).

2.1. Conservation and Balance Laws

Conservation and Balance Laws: Conservation of Mass (CM), Balance of Linear Momentum (BLM), Balance of Angular Momentum (BAM), First Law of Thermodynamics (FLT), and the Second Law of Thermodynamics (SLT) in Eulerian description can be written as [8] [9] (in fixed x-frame)

$\frac{D\bar{\rho }}{Dt}+\bar{\rho }\left(\bar{\nabla }\cdot \bar{v}\right)=0\left(\text{CM}\right)$ (6)

$\bar{\rho }\frac{D\bar{v}}{Dt}-\bar{\rho }{\bar{F}}^b-\bar{\nabla }\cdot \left({}^{\left(0\right)}\bar{\sigma }\right)=0\left(\text{BLM}\right)$ (7)

${є}_{ijk}{}^{\left(0\right)}\bar{\sigma }{}_{ij}=0\left(\text{BAM}\right)$ (8)

$\bar{\rho }\frac{D\left(\bar{e}\right)}{Dt}+\bar{\nabla }\cdot \bar{q}-{}^{\left(0\right)}\bar{\sigma }:\bar{D}=0\left(\text{FLT}\right)$ (9)

$\bar{\rho }\left(\frac{D\bar{\Phi }}{Dt}+\bar{\eta }\frac{D\bar{\theta }}{Dt}\right)+\frac{\bar{q}\cdot \bar{g}}{\bar{\theta }}-{}^{\left(0\right)}\bar{\sigma }:\bar{D}\le 0\left(\text{SLT}\right)$ (10)

in which $\bar{\rho }$ is density, $\bar{v}$ are velocities, ${\bar{F}}^b$ are body forces per unit mass, $є$ is the permutation tensor, $\bar{e}$ is specific internal energy, $\bar{q}$ is the heat vector, $\bar{D}$ is the symmetric part of velocity gradient tensor $\bar{L}$ (gradients of $\bar{v}$ with respect to $\bar{x}$ ), $\bar{\Phi }$ is Helmholtz free energy density, $\bar{\eta }$ is entropy density, $\bar{\theta }$ is absolute temperature, and $\bar{g}$ are the gradients of temperature $\bar{\theta }$ in x-frame.

This constitutive model contains $\bar{\rho }\left(1\right)$, $\bar{v}\left(3\right)$, ${}^{\left(0\right)}\bar{\sigma }\left(6\right)$, $\bar{q}\left(3\right)$, $\bar{\theta }\left(1\right)$, fourteen dependent variables in CM(1), BLM(3), FLT(1), and five partial differential equations. Thus, we need nine additional equations for closure. These are obtained from the constitutive theories for ${}^{\left(0\right)}\bar{\sigma }\left(6\right)$ and $\bar{q}\left(3\right)$. We note that (shown later) $\bar{e}=\bar{e}\left(\bar{\rho },\bar{\theta }\right)$, $\bar{\Phi }=\bar{\Phi }\left(\bar{\rho },\bar{\theta }\right)$ and $\bar{\eta }=\bar{\eta }\left(\bar{\rho },\bar{\theta }\right)$ ; hence $\bar{e}$, $\bar{\Phi }$, and $\bar{\eta }$ are not dependent variables in mathematical model.

2.2. Constitutive Theories

In the derivation of the constitutive theories, the entropy inequality aids in the initial choice of constitutive tensors and their argument tensors. Total deformation in a compressible polymeric fluid consists of volumetric deformation that results in change of volume without change in shape and distortional deformation that results in change of shape without change in volume. Thus, clearly the volumetric and distortional deformations are mutually exclusive, hence a single constitutive theory for the stress tensor cannot possibly describe both deformation physics. This necessitates additive decomposition of the Cauchy stress tensor ${}^{\left(0\right)}\bar{\sigma }$ into equilibrium ( ${}_e{}^{\left(0\right)}\bar{\sigma }$ ) and deviatoric ( ${}_d{}^{\left(0\right)}\bar{\sigma }$ ) Cauchy stress tensors.

${}^{\left(0\right)}\bar{\sigma }={}_e{}^{\left(0\right)}\bar{\sigma }+{}_d{}^{\left(0\right)}\bar{\sigma }$ (11)

The constitutive theory for ${}_e{}^{\left(0\right)}\bar{\sigma }$ describes volumetric deformation while the constitutive theory for ${}_d{}^{\left(0\right)}\bar{\sigma }$ addresses distortional deformation physics. From the entropy inequality (as well as other CBLs), initial choices of $\bar{\Phi }$, $\bar{\eta }$, ${}^{\left(0\right)}\bar{\sigma }$, and $\bar{q}$ as constitutive tensors are justified. The rate of work conjugate pair ${}^{\left(0\right)}\bar{\sigma }\mathrm{<mo>\; :\; </mo>}\bar{D}$ and compressible thermoviscoelastic physics suggest $\bar{\rho }$, $\bar{D}$, and $\bar{\theta }$ as possible argument tensors of ${}^{\left(0\right)}\bar{\sigma }$. Likewise, $\bar{\rho }$, $\bar{g}$ and $\bar{\theta }$ are possible argument tensors of $\bar{q}$. The choices of $\bar{\rho }$ and $\bar{\theta }$ as argument tensors of $\bar{\Phi }$ and $\bar{\eta }$ is obvious, while others can be initially considered based on the principle of equipresence. This gives us

${}^{\left(0\right)}\bar{\sigma }={}^{\left(0\right)}\bar{\sigma }\left(\bar{\rho }\mathrm{,}\bar{D}\mathrm{,}\bar{\theta }\right)$ (12)

$\bar{q}=\bar{q}\left(\bar{\rho }\mathrm{,}\bar{g}\mathrm{,}\bar{\theta }\right)$ (13)

$\bar{\Phi }=\bar{\Phi }\left(\bar{\rho }\mathrm{,}\bar{D}\mathrm{,}\bar{g}\mathrm{,}\bar{\theta }\right)$ (14)

$\bar{\eta }=\bar{\eta }\left(\bar{\rho }\mathrm{,}\bar{D}\mathrm{,}\bar{g}\mathrm{,}\bar{\theta }\right)$ (15)

The principle of equipresence for ${}^{\left(0\right)}\bar{\sigma }$ and $\bar{q}$ is ruled out based on conjugate pairs ${}^{\left(0\right)}\bar{\sigma }:\bar{D}$ and $\frac{\bar{q}\cdot \bar{g}}{\bar{\theta }}$ in the entropy inequality. We note that the phsyics described by ${}_e{}^{\left(0\right)}\bar{\sigma }$ and ${}_d{}^{\left(0\right)}\bar{\sigma }$ suggest the following [9].

${}_e{}^{\left(0\right)}\bar{\sigma }={}_e{}^{\left(0\right)}\bar{\sigma }\left(\bar{\rho }\mathrm{,}\bar{\theta }\right)$ (16)

${}_d{}^{\left(0\right)}\bar{\sigma }={}_d{}^{\left(0\right)}\bar{\sigma }\left(\bar{\rho }\mathrm{,}\bar{D}\mathrm{,}\bar{\theta }\right)$ (17)

Thus, based on Equations (16) and (17), Equation (12) is valid. We substitute Equation (11) into Equation (10) to obtain

$\bar{\rho }\left(\frac{D\bar{\Phi }}{Dt}+\bar{\eta }\frac{D\bar{\theta }}{Dt}\right)+\frac{\bar{q}\cdot \bar{g}}{\bar{\theta }}-{}_e{}^{\left(0\right)}\bar{\sigma }:\bar{D}-{}_d{}^{\left(0\right)}\bar{\sigma }:\bar{D}\le 0$ (18)

This inequality is the starting point for deriving constitutive theories.

2.2.1. Constitutive Theory for Equilibrium Stress ${}_e{}^{(\; 0\; )}\bar{\sigma }$

Using Equation (14), we can write the following using the chain rule of differentiation

$\frac{D\bar{\Phi }}{Dt}=\frac{\partial \bar{\Phi }}{\partial \bar{\rho }}\stackrel{\dot{}}{\bar{\rho }}+\frac{\partial \bar{\Phi }}{\partial \bar{D}}:\stackrel{\dot{}}{\bar{D}}+\frac{\partial \bar{\Phi }}{\partial \bar{g}}\cdot \stackrel{\dot{}}{\bar{g}}+\frac{\partial \bar{\Phi }}{\partial \bar{\theta }}\stackrel{\dot{}}{\bar{\theta }}$ (19)

Substituting from Equation (19) into Inequality (18)

$\bar{\rho }\left(\frac{\partial \bar{\Phi }}{\partial \bar{\rho }}\stackrel{\dot{}}{\bar{\rho }}+\frac{\partial \bar{\Phi }}{\partial \bar{D}}:\stackrel{\dot{}}{\bar{D}}+\frac{\partial \bar{\Phi }}{\partial \bar{g}}\cdot \stackrel{\dot{}}{\bar{g}}+\frac{\partial \bar{\Phi }}{\partial \bar{\theta }}\stackrel{\dot{}}{\bar{\theta }}+\bar{\eta }\stackrel{\dot{}}{\bar{\theta }}\right)+\frac{\bar{q}\cdot \bar{g}}{\bar{\theta }}-{}_e{}^{\left(0\right)}\bar{\sigma }:\bar{D}-{}_d{}^{\left(0\right)}\bar{\sigma }:\bar{D}\le 0$ (20)

From continuity (Equation (6)), we have (compressibility condition)

$\stackrel{\dot{}}{\bar{\rho }}=-\bar{\rho }\left(\bar{\nabla }\cdot \bar{v}\right)=-\bar{\rho }{\bar{D}}_{ki}{\delta }_{ik}=-\bar{\rho }\delta :\bar{D}$ (21)

Substituting Equation (21) into Inequality (20) and regrouping terms

$\begin{array}{l}\left(-{\bar{\rho }}^2\frac{\partial \bar{\Phi }}{\partial \bar{\rho }}\delta -{}_e{}^{\left(0\right)}\bar{\sigma }\right):\bar{D}-{}_d{}^{\left(0\right)}\bar{\sigma }:\bar{D}+\bar{\rho }\frac{\partial \bar{\Phi }}{\partial \bar{D}}:\stackrel{\dot{}}{\bar{D}}+\bar{\rho }\frac{\partial \bar{\Phi }}{\partial \bar{g}}\cdot \stackrel{\dot{}}{\bar{g}}\\ +\bar{\rho }\left(\bar{\eta }+\frac{\partial \bar{\Phi }}{\partial \bar{\theta }}\right)\stackrel{\dot{}}{\bar{\theta }}+\frac{\bar{q}\cdot \bar{g}}{\bar{\theta }}\le 0\end{array}$ (22)

The entropy inequality in Inequality (22) is satisfied for arbitrary but admissible choices of $\stackrel{\dot{}}{\bar{D}}$, $\stackrel{\dot{}}{\bar{g}}$, and $\stackrel{\dot{}}{\bar{\theta }}$ if their coefficients are set to zero, giving the following.

$\bar{\rho }\frac{\partial \bar{\Phi }}{\partial \bar{D}}=0\Rightarrow \bar{\Phi }\ne \bar{\Phi }\left(\bar{D}\right)$ (23)

$\bar{\rho }\frac{\partial \bar{\Phi }}{\partial \bar{g}}=0\Rightarrow \bar{\Phi }\ne \bar{\Phi }\left(\bar{g}\right)$ (24)

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

Equations (23) and (24) imply that $\bar{\Phi }$ is not a function of $\bar{D}$ and $\bar{g}$. Based on Equation (25), $\bar{\eta }$ is not a constitutive tensor as it is deterministic using Helmholtz free energy density, $\bar{\Phi }$. Using Equations (23)-(25), Inequality (22) can be written as

$\left(-{\bar{\rho }}^2\frac{\partial \bar{\Phi }}{\partial \bar{\rho }}\delta -{}_e{}^{\left(0\right)}\bar{\sigma }\right):\bar{D}-{}_d{}^{\left(0\right)}\bar{\sigma }:\bar{D}+\frac{\bar{q}\cdot \bar{g}}{\bar{\theta }}\le 0$ (26)

We note that since the constitutive theory for ${}_e{}^{\left(0\right)}\bar{\sigma }$ can only describe volumetric deformation physics, thus we can write Equation (16), i.e. ${}_e{}^{\left(0\right)}\bar{\sigma }={}_e{}^{\left(0\right)}\bar{\sigma }\left(\bar{\rho }\mathrm{,}\bar{\theta }\right)$. The constitutive tensors and their argument tensors in Equations (14) and (15) can now be modified and we can write

${}_e{}^{\left(0\right)}\bar{\sigma }={}_e{}^{\left(0\right)}\bar{\sigma }\left(\bar{\rho }\mathrm{,}\bar{\theta }\right)$ (27)

${}_d{}^{\left(0\right)}\bar{\sigma }={}_d{}^{\left(0\right)}\bar{\sigma }\left(\bar{\rho }\mathrm{,}\bar{D}\mathrm{,}\bar{\theta }\right)$ (28)

$\bar{q}=\bar{q}\left(\bar{\rho }\mathrm{,}\bar{g}\mathrm{,}\bar{\theta }\right)$ (29)

$\bar{\Phi }=\bar{\Phi }\left(\bar{\rho }\mathrm{,}\bar{\theta }\right)$ (30)

Based on Inequality (26) and Equations (27)-(30), we need to derive constitutive theories for ${}_e{}^{\left(0\right)}\bar{\sigma }$, ${}_d{}^{\left(0\right)}\bar{\sigma }$ and $\bar{q}$.

2.2.2. Constitutive Theory for ${}_e{}^{\left(0\right)}\bar{\sigma }$ for Compressible Polymeric Fluids

Based on Equation (30), we can set the coefficient of $\bar{D}$ in Inequality (26) to zero to obtain ${}_e{}^{\left(0\right)}\bar{\sigma }$ as a function of Helmholtz free energy (this in fact means that Inequality (26) is satisfied for arbitrary but admissible $\bar{D}$ if the coefficient of $\bar{D}$ in the first term of Inequality (26) is set to zero)

${}_e{}^{\left(0\right)}\bar{\sigma }=-{\bar{\rho }}^2\frac{\partial \bar{\Phi }}{\partial \bar{\rho }}\delta =\bar{p}\left(\bar{\rho }\mathrm{,}\bar{\theta }\right)\delta$ (31)

$\bar{p}\left(\bar{\rho }\mathrm{,}\bar{\theta }\right)=-{\bar{\rho }}^2\frac{\partial \bar{\Phi }}{\partial \bar{\rho }}$ (32)

and the entropy inequality reduces to

$-{}_d{}^{\left(0\right)}\bar{\sigma }:\bar{D}+\frac{\bar{q}\cdot \bar{g}}{\bar{\theta }}\le 0$ (33)

Equation (31) is the constitutive theory for the equilibrium Cauchy stress tensor for compressible polymeric fluids. $\bar{p}$ is generally referred to as equation of state. Experimental, empirical or analytical expressions for $\bar{p}\left(\bar{\rho }\mathrm{,}\bar{\theta }\right)$ are admissible as long as $\bar{p}$ is continuous and differentiable in $\bar{\rho }$ and $\bar{\theta }$. Equation (23) is called the reduced form of the entropy inequality.

2.2.3. Constitutive Theory for ${}_e{}^{\left(0\right)}\bar{\sigma }$ for Incompressible Polymeric Fluids

When the polymeric fluid is incompressible, $\bar{\rho }\left(\bar{x}\mathrm{,}t\right)=\rho \left(x\mathrm{,}t\right)={\rho }_0$, i.e. density remains constant. In this case

$\stackrel{\dot{}}{\bar{\rho }}=-\bar{\rho }\left(\bar{\nabla }\cdot \bar{v}\right)=0\left(\text{CM}\right)$ (34)

and

$\frac{\partial \bar{\Phi }\left(\bar{\rho }\mathrm{,}\bar{\theta }\right)}{\partial \bar{\rho }}=\frac{\partial \bar{\Phi }\left({\rho }_0\mathrm{,}\bar{\theta }\right)}{\partial \bar{\rho }}=0$ (35)

Hence, the constitutive theory for ${}_e{}^{\left(0\right)}\bar{\sigma }$ cannot be derived using Equation (31). Using Equation (35), Inequality (26) reduces to

$-{}_e{}^{\left(0\right)}\bar{\sigma }:\bar{D}-{}_d{}^{\left(0\right)}\bar{\sigma }:\bar{D}+\frac{\bar{q}\cdot \bar{g}}{\bar{\theta }}\le 0$ (36)

In order to derive constitutive theory for ${}_e{}^{\left(0\right)}\bar{\sigma }$ for incompressible polymeric fluids, we must introduce the incompressibility condition using conservation of mass in the Eulerian description (Equation (34)) into Inequality (36)

$\bar{\nabla }\cdot \bar{v}={\bar{D}}_{jj}={\bar{D}}_{ji}{\delta }_{ij}=\delta :\bar{D}=0$ (37)

We note that when Equation (37) holds, the following also holds

$\bar{p}\left(\bar{\theta }\right)\delta :\bar{D}=0$ (38)

$\bar{p}\left(\bar{\theta }\right)$ is a Lagrange multiplier (function of temperature $\bar{\theta }$ ). Adding Equation (38) to Inequality (36) and regrouping terms

$\left(\bar{p}\left(\bar{\theta }\right)\delta -{}_e{}^{\left(0\right)}\bar{\sigma }\right):\bar{D}-{}_d{}^{\left(0\right)}\bar{\sigma }:\bar{D}+\frac{\bar{q}\cdot \bar{g}}{\bar{\theta }}\le 0$ (39)

Inequality (39) holds for arbitrary but admissible $\bar{D}$ if the coefficient of $\bar{D}$ in the first term of Inequality (39) is set to zero, giving

${}_e{}^{\left(0\right)}\bar{\sigma }=\bar{p}\left(\bar{\theta }\right)\delta$ (40)

Equation (40) is the constitutive theory for incompressible, non-isothermal polymeric fluids. If the physics is isothermal, then $\bar{p}\left(\bar{\theta }\right)$ simply reduces to $\bar{p}$ in Equation (40). Inequality (39) now reduces to (reduced form)

$-{}_d{}^{\left(0\right)}\bar{\sigma }:\bar{D}+\frac{\bar{q}\cdot \bar{g}}{\bar{\theta }}\le 0$ (41)

2.2.4. Constitutive Theory for Deviatoric Stress ${}_d{}^{(\; 0\; )}\bar{\sigma }$

From Inequality (41), the rate of work conjugate pair ${}_d{}^{\left(0\right)}\bar{\sigma }:\bar{D}$ suggests that (as in Equation (28))

${}_d{}^{\left(0\right)}\bar{\sigma }={}_d{}^{\left(0\right)}\bar{\sigma }\left(\bar{\rho }\mathrm{,}\bar{D}\mathrm{,}\bar{\theta }\right)$ (42)

It is well known [3] [4] that in order to describe fading or short term memory in polymeric fluids, the constitutive theory for the deviatoric stress tensor must be a differential equation in time, otherwise existence of the memory modulus cannot be established. This necessitates that we must at the very least consider ${}_d{}^{\left(0\right)}\bar{\sigma }$ and ${}_d{}^{\left(1\right)}\bar{\sigma }$ in which ${}_d{}^{\left(1\right)}\bar{\sigma }$ is the constitutive tensor and ${}_d{}^{\left(0\right)}\bar{\sigma }$ is its argument tensor (in addition to others). If we consider convected time derivatives of ${}_d{}^{\left(0\right)}\bar{\sigma }$ of up to orders m, i.e. ${}_d{}^{\left(k\right)}\bar{\sigma }$ ; $k=0,1,\cdots ,m$ (basis independent notation), then we can generalize the choices of constitutive tensors and their argument tensors.

Likewise, since $\bar{D}$ is the first convected time derivative of Green’s strain tensor ${\epsilon }_{\left[0\right]}$ and also the first convected time derivative of the Almansi strain tensor [8] [9], we can also generalize the choice of $\bar{D}$ in Equation (42) by replacing it with ${}^{\left(j\right)}\gamma$ ; $j=1,2,\cdots ,n$, the convected time derivatives of the strain tensor (in basis independent notation). Thus, Equation (42) is replaced with

${}_d{}^{\left(m\right)}\bar{\sigma }={}_d{}^{\left(m\right)}\bar{\sigma }\left(\bar{\rho }\mathrm{,}{}^{\left(j\right)}\gamma \mathrm{,}{}_d{}^{\left(k\right)}\bar{\sigma }\mathrm{,}\bar{\theta }\right)\mathrm{<mo>\; ;\; </mo>}j=\mathrm{1,2,}\cdots \mathrm{,}n;k=\mathrm{0,1,}\cdots \mathrm{,}m-1$ (43)

We note that ${}_d{}^{\left(k\right)}\bar{\sigma }$ ; $k=\mathrm{0,1,}\cdots \mathrm{,}m$ and ${}^{\left(j\right)}\gamma$ ; $j=\mathrm{1,2,}\cdots \mathrm{,}n$ are all symmetric tensors of rank two, $\bar{\rho }$ and $\bar{\theta }$ are tensors of rank zero. The constitutive theory for ${}_d{}^{\left(m\right)}\bar{\sigma }$ can be derived using representation theorem [10] - [21].

Let ${}^{\sigma }\underset{˜}{G}{}^i$ ; $i=\mathrm{1,2,}\cdots \mathrm{,}N$ be the combined generators of the argument tensors of ${}_d{}^{\left(m\right)}\bar{\sigma }$ in Equation (43) that are symmetric tensors of rank two and let ${}^{\sigma }\underset{˜}{I}{}^j$ ; $j=1,2,\cdots ,M$ be the combined invariants of the same argument tensors. Then tensors $I$, ${}^{\sigma }\underset{˜}{G}{}^i$ ; $i=\mathrm{1,2,}\cdots \mathrm{,}N$ constitute the integrity, i.e. complete basis of the space of the constitutive tensor ${}_d{}^{\left(m\right)}\bar{\sigma }$, hence we can express ${}_d{}^{\left(m\right)}\bar{\sigma }$ as a linear combination of this basis.

${}_d{}^{\left(m\right)}\bar{\sigma }={}^{\sigma }\underset{˜}{\alpha }{}^{0}I+{\displaystyle \underset{i=1}{\overset{N}{\sum }}}{}^{\sigma }\underset{˜}{\alpha }{}^i\left({}^{\sigma }\underset{˜}{G}{}^i\right)$ (44)

in which

${}^{\sigma }\underset{˜}{\alpha }{}^i={}^{\sigma }\underset{˜}{\alpha }{}^i\left(\bar{\rho }\mathrm{,}{}^{\sigma }\underset{˜}{I}{}^j\mathrm{,}\bar{\theta }\right)\mathrm{<mo>\; ;\; </mo>}i=\mathrm{0,1,}\cdots \mathrm{,}N\mathrm{<mo>\; ;\; </mo>}j=\mathrm{1,2,}\cdots \mathrm{,}M$ (45)

In Equation (44), ${}^{\sigma }\underset{˜}{\alpha }{}^i$ ; $i=\mathrm{0,1,}\cdots \mathrm{,}N$ are coefficients of the linear combination (and not the material coefficients). Material coefficients are derived by considering Taylor series expansion of ${}^{\sigma }\underset{˜}{\alpha }{}^i$ ; $i=\mathrm{0,1,}\cdots \mathrm{,}N$ (based on the axiom of smooth neighborhood) in ${}^{\sigma }\underset{˜}{I}{}^j$ ; $j=\mathrm{1,2,}\cdots \mathrm{,}M$ about a known configuration $\underset{\_}{\Omega }$ and retaining only up to linear terms in ${}^{\sigma }\underset{˜}{I}{}^j$ ; $j=\mathrm{1,2,}\cdots \mathrm{,}M$ (for simplicity). The Taylor series expansion of ${}^{\sigma }\underset{˜}{\alpha }{}^i$, $i=\mathrm{0,1,}\cdots \mathrm{,}N$ in ${}^{\sigma }\underset{˜}{I}{}^j$, $j=\mathrm{1,2,}\cdots \mathrm{,}M$ about $\underset{\_}{\Omega }$ gives

${}^{\sigma }\underset{˜}{\alpha }{}^i={{}^{\sigma }\underset{˜}{\alpha }{}^i|}_{\underset{\_}{\Omega }}+{\displaystyle \underset{j=1}{\overset{M}{\sum }}}\left({\frac{\partial {}^{\sigma }\underset{˜}{\alpha }{}^i}{\partial {}^{\sigma }\underset{˜}{I}{}^j}|}_{\underset{\_}{\Omega }}\right)\left({}^{\sigma }\underset{˜}{I}{}^j-{{}^{\sigma }\underset{˜}{I}{}^j|}_{\underset{\_}{\Omega }}\right);i=0,1,\cdots ,N$ (46)

Substituting ${}^{\sigma }\underset{˜}{\alpha }{}^i$ ; $i=\mathrm{0,1,}\cdots \mathrm{,}N$ from Equation (46) into Equation (44) and collecting coefficients of $I$, ${}^{\sigma }\underset{˜}{I}{}^j$, ${}^{\sigma }\underset{˜}{G}{}^i$, and ${}^{\sigma }\underset{˜}{I}{}^j\left({}^{\sigma }\underset{˜}{G}{}^i\right)$ ; $i=1,2,\cdots ,N$, $j=\mathrm{1,2,}\cdots \mathrm{,}M$ and regrouping, we can obtain the following.

${}_d{}^{\left(m\right)}\bar{\sigma }={{\underset{˜}{\bar{\sigma }}}^0|}_{\underset{\_}{\Omega }}I+{\displaystyle \underset{i=1}{\overset{M}{\sum }}}{}^{\sigma }\underset{˜}{a}{}_j\left({}^{\sigma }\underset{˜}{I}{}^j\right)I+{\displaystyle \underset{i=1}{\overset{N}{\sum }}}{}^{\sigma }\underset{˜}{b}{}_i\left({}^{\sigma }\underset{˜}{G}{}^i\right)+{\displaystyle \underset{i=1}{\overset{N}{\sum }}}{\displaystyle \underset{i=1}{\overset{M}{\sum }}}{}^{\sigma }\underset{˜}{c}{}_{ij}\left({}^{\sigma }\underset{˜}{I}{}^j\right)\left({}^{\sigma }\underset{˜}{G}{}^i\right)$ (47)

in which

$\begin{array}{c}\begin{array}{l}{{\underset{˜}{\bar{\sigma }}}^0|}_{\underset{\_}{\Omega }}={{}^{\sigma }\underset{˜}{\alpha }{}^0|}_{\underset{\_}{\Omega }}-{\displaystyle \underset{j=1}{\overset{M}{\sum }}}{\frac{\partial {}^{\sigma }\underset{˜}{\alpha }{}^0}{\partial {}^{\sigma }\underset{˜}{I}{}^j}|}_{\underset{\_}{\Omega }}\left({{}^{\sigma }\underset{˜}{I}{}^j|}_{\underset{\_}{\Omega }}\right)\\ {}^{\sigma }\underset{˜}{a}{}_j={\frac{\partial {}^{\sigma }\underset{˜}{\alpha }{}^0}{\partial {}^{\sigma }\underset{˜}{I}{}^j}|}_{\underset{\_}{\Omega }}\\ {}^{\sigma }\underset{˜}{b}{}_i={{}^{\sigma }\underset{˜}{\alpha }{}^i|}_{\underset{\_}{\Omega }}-{\displaystyle \underset{j=1}{\overset{M}{\sum }}}{\frac{\partial {}^{\sigma }\underset{˜}{\alpha }{}^i}{\partial {}^{\sigma }\underset{˜}{I}{}^j}|}_{\underset{\_}{\Omega }}\left({{}^{\sigma }\underset{˜}{I}{}^j|}_{\underset{\_}{\Omega }}\right)\\ {}^{\sigma }\underset{˜}{c}{}_{ij}={\frac{\partial {}^{\sigma }\underset{˜}{\alpha }{}^i}{\partial {}^{\sigma }\underset{˜}{I}{}^j}|}_{\underset{\_}{\Omega }}\end{array}\end{array}\};\begin{array}{c}i=1,2,\cdots ,N\\ j=1,2,\cdots ,M\end{array}$ (48)

Equation (47) is the constitutive theory for the deviatoric Cauchy stress tensor based on the integrity, i.e. complete basis in which ${}^{\sigma }\underset{˜}{a}{}_j$, ${}^{\sigma }\underset{˜}{b}{}_i$, ${}^{\sigma }\underset{˜}{c}{}_{ij}$ ; $i=1,2,\cdots ,N$ ; $j=1,2,\cdots ,M$ are $\left(N+M+NM\right)$ material coefficients defined in a known configuration $\underset{\_}{\Omega }$. The material coefficients can be functions of the invariants ${}^{\sigma }\underset{˜}{I}{}^j$ ; $j=1,2,\cdots ,M$, $\bar{\rho }$, and $\bar{\theta }$ (based on Equation (45)).

2.3. Derivation of Currently Used Constitutive Theories for Polymeric Fluids

Equation (47) contains many material coefficients, some of which may not be significant for specific types of polymers. Nonetheless, we point out that Equation (47) represents the totality of all possible constitutive theories for polymeric fluids as it is based on a complete basis of the space of the constitutive tensor. In this section, we first show derivations of Maxwell, Oldroyd-B, and Giesekus constitutive theories using Equation (47) by selective choices of generators and invariants. This is followed by enhancement of these theories by incorporating additional generators and/or invariants (that are permissible based on Equation (47)) in the existing constitutive theories based on Equation (47) to remedy some obvious deficiencies in them, especially for dilute polymeric fluids. The currently used Maxwell, Oldroyd-B and Giesekus constitutive theories for the stress tensor are listed below [3] [4] [8] [9] for the compressible case. We refer to these as model A in the model problem studies.

${}_d{}^{\left(0\right)}\bar{\sigma }+\lambda \left({}_d{}^{\left(1\right)}\bar{\sigma }\right)=2\eta \left({}^{\left(1\right)}\gamma \right)+\kappa \text{tr}\left({}^{\left(1\right)}\gamma \right)I\left(\text{Maxwell}\right)$ (49)

${}_d{}^{\left(0\right)}\bar{\sigma }+\lambda \left({}_d{}^{\left(1\right)}\bar{\sigma }\right)=2\eta \left({}^{\left(1\right)}\gamma \right)+\kappa \text{tr}\left({}^{\left(1\right)}\gamma \right)I+2\eta {\lambda }_2\left({}^{\left(2\right)}\gamma \right)\left(\text{Oldroyd-B}\right)$ (50)

${}_d{}^{\left(0\right)}\bar{\sigma }+\lambda \left({}_d{}^{\left(1\right)}\bar{\sigma }\right)=2\eta \left({}^{\left(1\right)}\gamma \right)+\kappa \text{tr}\left({}^{\left(1\right)}\gamma \right)I+\frac{\lambda }{\eta }\alpha {\left({}_d{}^{\left(0\right)}\bar{\sigma }\right)}^2\left(\text{Giesekus}\right)$ (51)

Remarks

1) From Equation (49), we note that the Maxwell constitutive model only contains the first convected time derivatives of the deviatoric stress tensor and the strain tensor (Green’s or Almansi). Thus, this constitutive model can be obtained by using Equation (47) with $n=m=1$, and by deleting the terms in it other than those in Equation (49).

2) From the Oldroyd-B constitutive model in Equation (50), we note that it contains only up to first convected time derivatives of the deviatoric stress tensor but contains up to second convected time derivatives of the strain tensor. Hence, the constitutive model is a subset of the constitutive theory in Equation (47) for $n=2$ and $m=1$, thus can be obtained from the general form in Equation (47) for $n=2$ and $m=1$ by deleting terms other than those that appear in Equation (50).

3) The Giesekus constitutive model is also a subset of Equation (47) for $n=1$ and $m=1$. This model is the same as the Maxwell model except for the nonlinear term in ${}_d{}^{\left(0\right)}\bar{\sigma }$.

4) From Equations (1)-(3), we conclude that Maxwell, Oldroyd-B and Giesekus constitutive models are a subset of the general constitutive model in Equation (47) for $n=2$ and $m=1$. We remark that the complete constitutive model (based on the integrity for $n=2$ and $m=1$ ) contains many more generators and invariants of the argument tensors than those appearing in Equations (49)-(51). The generators and invariants in Equation (47) that are not considered in Equations (49)-(51) provide the basis for enhancing constitutive theories in Equations (49)-(51).

5) In the following, we first present a single constitutive theory using Equation (47) with $n=2$ and $m=1$ that describes all three constitutive models (Equations (49)-(51)) used currently.

6) The Giesekus constitutive model used currently utilizes additive decomposition of ${}_d{}^{\left(0\right)}\bar{\sigma }$ into ${\left({}_d{}^{\left(0\right)}\bar{\sigma }\right)}_s$ and ${\left({}_d{}^{\left(0\right)}\bar{\sigma }\right)}_p$, the solvent and polymer deviatoric stresses. Newton’s law of viscosity is used to describe constitutive theory for ${\left({}_d{}^{\left(0\right)}\bar{\sigma }\right)}_s$ while ${}_d{}^{\left(0\right)}\bar{\sigma }$ in Equation (51) is replaced with ${\left({}_d{}^{\left(0\right)}\bar{\sigma }\right)}_p$. It has been shown [8] [9] that this decomposition and the use of ${\left({}_d{}^{\left(0\right)}\bar{\sigma }\right)}_s$ and ${\left({}_d{}^{\left(0\right)}\bar{\sigma }\right)}_p$ in the constitutive theories as described above is not supported by classical continuum mechanics. Thus, we do not use their additive decomposition, hence ${}_d{}^{\left(0\right)}\bar{\sigma }$ and ${}_d{}^{\left(1\right)}\bar{\sigma }$ are maintained in Equation (51).

For deriving Equations (49)-(51) constitutive theories, we consider $n=2$ and $m=1$, that is

${}_d{}^{\left(1\right)}\bar{\sigma }={}_d{}^{\left(1\right)}\bar{\sigma }\left(\bar{\rho }\mathrm{,}{}^{\left(1\right)}\gamma \mathrm{,}{}^{\left(2\right)}\gamma \mathrm{,}{}_d{}^{\left(0\right)}\bar{\sigma }\mathrm{,}\bar{\theta }\right)$ (52)

The use of representation theorem [10] - [21] to derive the constitutive theory for ${}_d{}^{\left(1\right)}\bar{\sigma }$, a symmetric tensor of rank two, requires that we must consider the combined generators of the argument tensors of ${}_d{}^{\left(1\right)}\bar{\sigma }$ in Equation (52) that are symmetric tensors of rank two as well as their combined invariants. The generators from each argument tensor ${}^{\left(1\right)}\gamma$, ${}^{\left(2\right)}\gamma$, and ${}_d{}^{\left(0\right)}\bar{\sigma }$ that are symmetric tensors of rank two are ${}^{\left(1\right)}\gamma$, ${\left({}^{\left(1\right)}\gamma \right)}^2$ ; ${}^{\left(2\right)}\gamma$, ${\left({}^{\left(2\right)}\gamma \right)}^2$, and ${}_d{}^{\left(0\right)}\bar{\sigma }$, ${\left({}_d{}^{\left(0\right)}\bar{\sigma }\right)}^2$ and their invariants (principal) are $I_{{}^{\left(1\right)}\gamma }$, $I{I}_{{}^{\left(1\right)}\gamma }$, $II{I}_{{}^{\left(1\right)}\gamma }$ ; $I_{{}^{\left(2\right)}\gamma }$, $I{I}_{{}^{\left(2\right)}\gamma }$, $II{I}_{{}^{\left(2\right)}\gamma }$, and $I_{{}_d{}^{\left(0\right)}\bar{\sigma }}$, $I{I}_{{}_d{}^{\left(0\right)}\bar{\sigma }}$, $II{I}_{{}_d{}^{\left(0\right)}\bar{\sigma }}$. In addition to these, there are combined generators and invariants of the argument tensors considering these tensors in sets, two at a time and three at a time, i.e. ${}^{\left(1\right)}\gamma$, ${}^{\left(2\right)}\gamma$ ; ${}^{\left(2\right)}\gamma$, ${}_d{}^{\left(0\right)}\bar{\sigma }$ ; ${}^{\left(1\right)}\gamma$, ${}_d{}^{\left(0\right)}\bar{\sigma }$, and ${}^{\left(1\right)}\gamma$, ${}^{\left(2\right)}\gamma$, ${}_d{}^{\left(0\right)}\bar{\sigma }$. $\bar{\rho }$ and $\bar{\theta }$ are tensors of rank zero, hence they do not contribute to the combined generators and invariants. The constitutive theories in Equations (49)-(51) do not contain combined generators, hence we need not consider these in this derivation. Thus, a general constitutive theory for ${}_d{}^{\left(1\right)}\bar{\sigma }$ based on generators of ${}^{\left(1\right)}\gamma$, ${}^{\left(2\right)}\gamma$, and ${}_d{}^{\left(0\right)}\bar{\sigma }$ would be

$\begin{array}{c}{}_d{}^{\left(1\right)}\bar{\sigma }={}^{\sigma }\underset{˜}{\alpha }{}^{0}I+{}^{\sigma }\underset{˜}{\alpha }{}^1\left({}^{\left(1\right)}\gamma \right)+{}^{\sigma }\underset{˜}{\alpha }{}^2{\left({}^{\left(1\right)}\gamma \right)}^2+{}^{\sigma }\underset{˜}{\alpha }{}^3\left({}^{\left(2\right)}\gamma \right)+{}^{\sigma }\underset{˜}{\alpha }{}^4{\left({}^{\left(2\right)}\gamma \right)}^2\\ +{}^{\sigma }\underset{˜}{\alpha }{}^5\left({}_d{}^{\left(0\right)}\bar{\sigma }\right)+{}^{\sigma }\underset{˜}{\alpha }{}^6{\left({}_d{}^{\left(0\right)}\bar{\sigma }\right)}^2\end{array}$ (53)

in which coefficients ${}^{\sigma }\underset{˜}{\alpha }{}^i$ : $i=\mathrm{0,1,}\cdots \mathrm{,6}$ in the linear combination in Equation (53) are functions of $\bar{\rho }$, $\bar{\theta }$ and the invariants of ${}^{\left(1\right)}\gamma$, ${}^{\left(2\right)}\gamma$, and ${}_d{}^{\left(0\right)}\bar{\sigma }$. The material coefficients are established by considering a Taylor series expansion of ${}^{\sigma }\underset{˜}{\alpha }{}^i$ : $i=\mathrm{0,1,}\cdots \mathrm{,6}$ in the invariants of ${}^{\left(1\right)}\gamma$, ${}^{\left(2\right)}\gamma$, and ${}_d{}^{\left(0\right)}\bar{\sigma }$ about a known configuration $\underset{\_}{\Omega }$, retaining only up to linear terms in the invariants. This constitutive theory will contain generators (as in Equation (53)) and invariants of ${}^{\left(1\right)}\gamma$, ${}^{\left(2\right)}\gamma$, ${}_d{}^{\left(0\right)}\bar{\sigma }$ as well as the products of the generators and invariants. This constitutive theory is the most comprehensive constitutive theory based on Equation (53). We use this constitutive theory as a guide for modifications of Equations (49)-(51). To derive Equations (49)-(51) as a subset of this single constitutive theory, we retain generators ${}^{\left(1\right)}\gamma$, ${}^{\left(2\right)}\gamma$, and ${}_d{}^{\left(0\right)}\bar{\sigma }$ and ${\left({}_d{}^{\left(0\right)}\bar{\sigma }\right)}^2$ in Equation (53), which reduces Equation (53) to (redefining coefficients in the linear combination)

${}_d{}^{\left(1\right)}\bar{\sigma }={}^{\sigma }\underset{˜}{\alpha }{}^{0}I+{}^{\sigma }\underset{˜}{\alpha }{}^1\left({}^{\left(1\right)}\gamma \right)+{}^{\sigma }\underset{˜}{\alpha }{}^2\left({}^{\left(2\right)}\gamma \right)+{}^{\sigma }\underset{˜}{\alpha }{}^3\left({}_d{}^{\left(0\right)}\bar{\sigma }\right)+{}^{\sigma }\underset{˜}{\alpha }{}^4{\left({}_d{}^{\left(0\right)}\bar{\sigma }\right)}^2$ (54)

${}^{\sigma }\underset{˜}{\alpha }{}^i$ : $i=0,1,\cdots ,4$ are functions of $\bar{\rho }$, $\bar{\theta }$ and invariants of ${}^{\left(1\right)}\gamma$, ${}^{\left(2\right)}\gamma$, and ${}_d{}^{\left(0\right)}\bar{\sigma }$. Substituting the Taylor series expansion of ${}^{\sigma }\underset{˜}{\alpha }{}^i$ : $i=0,1,\cdots ,4$ in the invariants of ${}^{\left(1\right)}\gamma$, ${}^{\left(2\right)}\gamma$, and ${}_d{}^{\left(0\right)}\bar{\sigma }$ about a known configuration $\underset{\_}{\Omega }$ (retaining only up to linear terms in the invariants, for simplicity) and retaining only those generators and invariants that appear in Equations (49)-(51), we can obtain

${}_d{}^{\left(1\right)}\bar{\sigma }={\bar{\sigma }}_{0}I+a_1\left({}^{\left(1\right)}\gamma \right)+a_2\left({}^{\left(2\right)}\gamma \right)+a_3\left({}_d{}^{\left(0\right)}\bar{\sigma }\right)+a_4{\left({}_d{}^{\left(0\right)}\bar{\sigma }\right)}^2$ (55)

By dividing throughout by $a_1$, rearranging terms and defining new coefficients, we can obtain the following

${}_d{}^{\left(0\right)}\bar{\sigma }+\lambda \left({}_d{}^{\left(1\right)}\bar{\sigma }\right)=2\eta \left({}^{\left(1\right)}\gamma \right)+2\eta {\lambda }_2\left({}^{\left(2\right)}\gamma \right)+\frac{\lambda }{\eta }\alpha {\left({}_d{}^{\left(0\right)}\bar{\sigma }\right)}^2+\kappa \text{tr}\left({}^{\left(1\right)}\gamma \right)I$ (56)

in which $\eta$ is viscosity, $\lambda$ is relaxation time, ${\lambda }_2$ is retardation time, $\kappa$ is second viscosity and $\alpha$ is mobility factor. If the polymer is incompressible, $\text{tr}\left({}^{\left(1\right)}\gamma \right)=0$, in which case Equation (56) yields Equations (49)-(51) when

1) ${\lambda }_2=0$, $\alpha =0$ ; Maxwell

2) $\alpha =0$ ; Oldroyd-B

3) ${\lambda }_2=0$ ; Giesekus

we note that Equation (56) is in ${}_d{}^{\left(0\right)}\bar{\sigma }$ and ${}_d{}^{\left(1\right)}\bar{\sigma }$, not ${\left({}_d{}^{\left(0\right)}\bar{\sigma }\right)}_p$ and ${\left({}_d{}^{\left(1\right)}\bar{\sigma }\right)}_p$ as used currently.

Remarks

The decomposition

${}_d{}^{\left(0\right)}\bar{\sigma }={\left({}_d{}^{\left(0\right)}\bar{\sigma }\right)}_s+{\left({}_d{}^{\left(0\right)}\bar{\sigma }\right)}_p$ (57)

suggests that we substitute this into the entropy inequality to determine how to derive constitutive theories for ${\left({}_d{}^{\left(0\right)}\bar{\sigma }\right)}_s$ and ${\left({}_d{}^{\left(0\right)}\bar{\sigma }\right)}_p$. Using the reduced form of Inequality (41), we can write

$-\left({\left({}_d{}^{\left(0\right)}\bar{\sigma }\right)}_s+{\left({}_d{}^{\left(0\right)}\bar{\sigma }\right)}_p\right):\bar{D}+\frac{\bar{q}\cdot \bar{g}}{\bar{\theta }}\le 0$ (58)

At this stage, it is perhaps convenient to conclude that

${\left({}_d{}^{\left(1\right)}\bar{\sigma }\right)}_s={\left({}_d{}^{\left(1\right)}\bar{\sigma }\right)}_s\left(\bar{\rho }\mathrm{,}\bar{D}\mathrm{,}\bar{\theta }\right)={\left({}_d{}^{\left(1\right)}\bar{\sigma }\right)}_s\left(\bar{\rho }\mathrm{,}{}^{\left(1\right)}\gamma \mathrm{,}\bar{\theta }\right)$ (59)

and ${\left({}_d{}^{\left(1\right)}\bar{\sigma }\right)}_p={\left({}_d{}^{\left(1\right)}\bar{\sigma }\right)}_p\left(\bar{\rho }\mathrm{,}{}^{\left(1\right)}\gamma \mathrm{,}{}^{\left(2\right)}\gamma \mathrm{,}{}_d{}^{\left(0\right)}\bar{\sigma }\mathrm{,}\bar{\theta }\right)$ (60)

This assumption may not be reflective of the true physics due to the fact that a polymer is an isotropic, homogeneous fluid which has its own properties, and constitutive theories that relates to the constituents (solvent and polymer) as done in Equations (59) and (60) are not possible. In references [8] [9], model problem studies are presented to demonstrate that the use of Equation (59) and Equation (60) instead of Equation (56) leads to drastically different results.

3. Enhancement of Constitutive Theories Used Currently

It is well known that the mathematical models for polymers utilizing constitutive theory in Equation (56) are deficient in simulating normal stresses perpendicular to the direction of flow, for example, in the case of flow between parallel plates in $x_1$ direction, ${}_d{}^{\left(1\right)}\bar{\sigma }{}_{11}$ is in fairly good agreement with experiments but ${}_d{}^{\left(1\right)}\bar{\sigma }{}_{22}$ is zero in Maxwell and Oldroyd-B models. The constitutive theory in Equation (47) of orders m and n based on the integrity (complete basis) contains all possible generators and invariants due to the argument tensors of ${}_d{}^{\left(m\right)}\bar{\sigma }$ in Equation (43). In the current constitutive theories ( $m=1$, $n=2$ ) for deviatoric Cauchy stress tensor, we consider

${}_d{}^{\left(1\right)}\bar{\sigma }={}_d{}^{\left(1\right)}\bar{\sigma }\left(\bar{\rho }\mathrm{,}{}^{\left(1\right)}\gamma \mathrm{,}{}^{\left(2\right)}\gamma \mathrm{,}{}_d{}^{\left(0\right)}\bar{\sigma }\mathrm{,}\bar{\theta }\right)={}_d{}^{\left(1\right)}\bar{\sigma }\left(\bar{\rho }\mathrm{,}\bar{D}\mathrm{,}{}^{\left(2\right)}\gamma \mathrm{,}{}_d{}^{\left(0\right)}\bar{\sigma }\mathrm{,}\bar{\theta }\right)$ (61)

The combined generators of the argument tensors of ${}_d{}^{\left(1\right)}\bar{\sigma }$ in Equation (61) that are symmetric tensors of rank two are:

$\text{due}\text{to}{}^{\left(1\right)}\gamma \mathrm{<mo>\; :\; </mo>}{}^{\sigma }\underset{˜}{G}{}^1={}^{\left(1\right)}\gamma \mathrm{<mo>\; ;\; </mo>}{}^{\sigma }\underset{˜}{G}{}^2={\left({}^{\left(1\right)}\gamma \right)}^2\text{or}\bar{D}\mathrm{,}{\bar{D}}^2$ (62)

$\text{due}\text{to}{}^{\left(2\right)}\gamma \mathrm{<mo>\; :\; </mo>}{}^{\sigma }\underset{˜}{G}{}^3={}^{\left(2\right)}\gamma \mathrm{<mo>\; ;\; </mo>}{}^{\sigma }\underset{˜}{G}{}^4={\left({}^{\left(2\right)}\gamma \right)}^2$ (63)

$\text{due}\text{to}{}_d{}^{\left(0\right)}\bar{\sigma }\mathrm{<mo>\; :\; </mo>}{}^{\sigma }\underset{˜}{G}{}^5={}_d{}^{\left(0\right)}\bar{\sigma }\mathrm{<mo>\; ;\; </mo>}{}^{\sigma }\underset{˜}{G}{}^6={\left({}_d{}^{\left(0\right)}\bar{\sigma }\right)}^2$ (64)

$\begin{array}{l}\text{due}\text{to}{}^{\left(1\right)}\gamma \text{and}{}^{\left(2\right)}\gamma \mathrm{<mo>\; :\; </mo>}{}^{\sigma }\underset{˜}{G}{}^7={}^{\left(1\right)}\gamma \cdot {}^{\left(2\right)}\gamma +{}^{\left(2\right)}\gamma \cdot {}^{\left(1\right)}\gamma =\bar{D}\cdot {}^{\left(2\right)}\gamma +{}^{\left(2\right)}\gamma \cdot \bar{D}\\ {}^{\sigma }\underset{˜}{G}{}^8={\left({}^{\left(1\right)}\gamma \right)}^2\cdot {}^{\left(2\right)}\gamma +{}^{\left(2\right)}\gamma \cdot {\left({}^{\left(1\right)}\gamma \right)}^2={\bar{D}}^2\cdot {}^{\left(2\right)}\gamma +{}^{\left(2\right)}\gamma \cdot {\bar{D}}^2\\ {}^{\sigma }\underset{˜}{G}{}^9={}^{\left(1\right)}\gamma \cdot {\left({}^{\left(2\right)}\gamma \right)}^2+{\left({}^{\left(2\right)}\gamma \right)}^2\cdot {}^{\left(1\right)}\gamma =\bar{D}\cdot {\left({}^{\left(2\right)}\gamma \right)}^2+{\left({}^{\left(2\right)}\gamma \right)}^2\cdot \bar{D}\\ \text{due}\text{to}{}^{\left(1\right)}\gamma \text{and}{}_d{}^{\left(0\right)}\bar{\sigma }\mathrm{<mo>\; :\; </mo>}{}^{\sigma }\underset{˜}{G}{}^{10}\mathrm{,}{}^{\sigma }\underset{˜}{G}{}^{11}\mathrm{,}{}^{\sigma }\underset{˜}{G}{}^{12}\left(\text{similar}\text{to}{}^{\sigma }\underset{˜}{G}{}^7\mathrm{,}{}^{\sigma }\underset{˜}{G}{}^8\mathrm{,}{}^{\sigma }\underset{˜}{G}{}^9\right)\\ \text{due}\text{to}{}^{\left(2\right)}\gamma \text{and}{}_d{}^{\left(0\right)}\bar{\sigma }\mathrm{<mo>\; :\; </mo>}{}^{\sigma }\underset{˜}{G}{}^{13}\mathrm{,}{}^{\sigma }\underset{˜}{G}{}^{14}\mathrm{,}{}^{\sigma }\underset{˜}{G}{}^{15}\left(\text{similar}\text{to}{}^{\sigma }\underset{˜}{G}{}^7\mathrm{,}{}^{\sigma }\underset{˜}{G}{}^8\mathrm{,}{}^{\sigma }\underset{˜}{G}{}^9\right)\end{array}$ (65)

$\text{due}\text{to}{}^{\left(1\right)}\gamma \mathrm{,}{}^{\left(2\right)}\gamma \text{and}{}_d{}^{\left(0\right)}\bar{\sigma }$ : can be obtained from [8] [9]

and the combined invariants of the same argument tensors are

$\text{due}\text{to}{}^{\left(1\right)}\gamma \mathrm{<mo>\; :\; </mo>}{}^{\sigma }\underset{˜}{I}{}^1=\text{tr}\left(\bar{D}\right)\mathrm{<mo>\; ;\; </mo>}{}^{\sigma }\underset{˜}{I}{}^2=\text{tr}\left({\bar{D}}^2\right)\mathrm{<mo>\; ;\; </mo>}{}^{\sigma }\underset{˜}{I}{}^3=\text{tr}\left({\bar{D}}^3\right)$ (66)

$\text{due}\text{to}{}^{\left(2\right)}\gamma \mathrm{<mo>\; :\; </mo>}{}^{\sigma }\underset{˜}{I}{}^4=\text{tr}\left({}^{\left(2\right)}\gamma \right)\mathrm{<mo>\; ;\; </mo>}{}^{\sigma }\underset{˜}{I}{}^5=\text{tr}{\left({}^{\left(2\right)}\gamma \right)}^2\mathrm{<mo>\; ;\; </mo>}{}^{\sigma }\underset{˜}{I}{}^6=\text{tr}{\left({}^{\left(2\right)}\gamma \right)}^3$ (67)

$\begin{array}{l}\text{due}\text{to}{}^{\left(1\right)}\gamma \mathrm{,}{}^{\left(2\right)}\gamma \left(\text{or}\bar{D}\mathrm{,}{}^{\left(2\right)}\gamma \right)\mathrm{<mo>\; :\; </mo>}{}^{\sigma }\underset{˜}{I}{}^7=\text{tr}\left(\bar{D}\cdot {}^{\left(2\right)}\gamma \right)\mathrm{,}{}^{\sigma }\underset{˜}{I}{}^8=\text{tr}\left({\bar{D}}^2\cdot {}^{\left(2\right)}\gamma \right)\\ {}^{\sigma }\underset{˜}{I}{}^9=\text{tr}\left(\bar{D}\cdot {\left({}^{\left(2\right)}\gamma \right)}^2\right)\mathrm{,}{}^{\sigma }\underset{˜}{I}{}^{10}=\text{tr}\left({\bar{D}}^2\cdot {\left({}^{\left(2\right)}\gamma \right)}^2\right)\\ {}^{\sigma }\underset{˜}{I}{}^{11}=\text{tr}\left(\bar{D}\cdot {}^{\left(2\right)}\gamma +{}^{\left(2\right)}\gamma \cdot \bar{D}\right)\\ {}^{\sigma }\underset{˜}{I}{}^{12}=\text{tr}\left(\bar{D}\cdot {}^{\left(2\right)}\gamma -{}^{\left(2\right)}\gamma \cdot \bar{D}\right)\end{array}$ (68)

$\begin{array}{l}\text{due}\text{to}{}^{\left(1\right)}\gamma \mathrm{,}{}_d{}^{\left(0\right)}\bar{\sigma }\mathrm{<mo>\; :\; </mo>}{}^{\sigma }\underset{˜}{I}{}^j\mathrm{<mo>\; :\; </mo>}j\mathrm{<mo>\; =\; </mo>13,14,}\cdots \mathrm{,18}\left(\text{similar}\text{to}{}^{\sigma }\underset{˜}{I}{}^j\mathrm{<mo>\; :\; </mo>}j=\mathrm{7,8,}\cdots \mathrm{,12}\right)\\ \text{due}\text{to}{}^{\left(2\right)}\gamma \mathrm{,}{}_d{}^{\left(0\right)}\bar{\sigma }\mathrm{<mo>\; :\; </mo>}{}^{\sigma }\underset{˜}{I}{}^j\mathrm{<mo>\; :\; </mo>}j\mathrm{<mo>\; =\; </mo>19,20,}\cdots \mathrm{,24}\left(\text{similar}\text{to}{}^{\sigma }\underset{˜}{I}{}^j\mathrm{<mo>\; :\; </mo>}j=\mathrm{7,8,}\cdots \mathrm{,12}\right)\end{array}$

$\text{due}\text{to}{}^{\left(1\right)}\gamma \mathrm{,}{}^{\left(2\right)}\gamma \text{and}{}_d{}^{\left(0\right)}\bar{\sigma }$ : can be obtained from [8] [9]

In the currently used constitutive theories, we use generators: ${}^{\sigma }\underset{˜}{G}{}^1=\bar{D}$, ${}^{\sigma }\underset{˜}{G}{}^3={}^{\left(2\right)}\gamma$ and ${}^{\sigma }\underset{˜}{G}{}^6={\left({}_d{}^{\left(0\right)}\bar{\sigma }\right)}^2$. First invariants of tensors $\bar{D}$, ${}^{\left(2\right)}\gamma$, and ${}_d{}^{\left(0\right)}\bar{\sigma }$, i.e. $\text{tr}\left(\bar{D}\right)$, $\text{tr}\left({}^{\left(2\right)}\gamma \right)$, and $\text{tr}\left({}_d{}^{\left(0\right)}\bar{\sigma }\right)$ could have been used but are neglected in Equation (56). $\text{tr}\left(\bar{D}\right)$ in Equation (56) does not appear due to the incompressibility assumption, but must be included in Equation (56) if the polymer is compressible.

The constitutive theory in Equation (56) can be enhanced by using generators and invariants in Equations (62)-(68). Since $\bar{D}$ is a fundamental measure of deformation rate, our first choice must include the generators and invariants related to $\bar{D}$ that are not present in Equation (56). This suggests that we must consider the addition of generators ${\bar{D}}^2$ and invariant $\text{tr}\left({\bar{D}}^2\right)$ in Equation (56). Thus, we must explore the new constitutive theory for Maxwell, Oldroyd-B, and Giesekus polymeric fluids,

${}_d{}^{\left(0\right)}\bar{\sigma }+\lambda \left({}_d{}^{\left(1\right)}\bar{\sigma }\right)=2\eta \bar{D}+2\eta {\lambda }_2\left({}^{\left(2\right)}\gamma \right)+\frac{\lambda }{\eta }\alpha {\left({}_d{}^{\left(0\right)}\bar{\sigma }\right)}^2+{\eta }_1{\left(\bar{D}\right)}^2+{\eta }_3\text{tr}\left({\bar{D}}^2\right)I$ (69)

Determination of ${\eta }_1$ and ${\eta }_3$, i.e. calibration of Equation (69) requires experiments. In the following we present two model problem studies: 1) fully developed flow between parallel plates and 2) a lid-driven square cavity.

4. Complete Mathematical Model and Its Dimensionless Form

For isothermal, incompressible flow, conservation of mass, balance of linear momentum, and the constitutive theories are given by (using contravariant Cauchy stress tensor, ${\bar{\sigma }}^{(\; 0\; )}$

$\frac{\partial \widehat{\stackrel{¯}{\rho }}}{\partial t}+\widehat{\stackrel{¯}{\rho }}\left(\widehat{\stackrel{¯}{\nabla }}\cdot \widehat{\stackrel{¯}{v}}\right)=0$ (70)

$\widehat{\stackrel{¯}{\rho }}\frac{D\widehat{\stackrel{¯}{v}}}{Dt}+\widehat{\stackrel{¯}{\rho }}{\bar{F}}^b+\widehat{\stackrel{¯}{\nabla }}\cdot \widehat{\stackrel{¯}{p}}-{\left({}_d{\widehat{\stackrel{¯}{\sigma }}}^{\left(0\right)}\right)}^{\text{T}}\cdot \widehat{\stackrel{¯}{\nabla }}=0$ (71)

${}_d\widehat{\stackrel{¯}{\sigma }}{}^{\left(0\right)}+\hat{\lambda }\left({}_d\widehat{\stackrel{¯}{\sigma }}{}^{\left(1\right)}\right)=2\hat{\eta }\widehat{\stackrel{¯}{D}}+2\hat{\eta }{\hat{\lambda }}_2\left({\hat{\gamma }}_{\left(2\right)}\right)+\frac{\hat{\lambda }}{\hat{\eta }}\alpha {\left({}_d\widehat{\stackrel{¯}{\sigma }}{}^{\left(0\right)}\right)}^2+{\hat{\eta }}_1{\left(\widehat{\stackrel{¯}{D}}\right)}^2+{\hat{\eta }}_3\text{tr}\left({\widehat{\stackrel{¯}{D}}}^2\right)I$ (72)

Hat (^{^}) over all quantities indicate that they have their usual dimensions (or units). We choose the following reference quantities (with zero subscript) and dimensionless variables (without hat).

$\begin{array}{l}\bar{x}=\frac{\widehat{\stackrel{¯}{x}}}{L_0}\mathrm{<mo>\; ;\; </mo>}\bar{v}=\frac{\widehat{\stackrel{¯}{v}}}{v_0}\mathrm{<mo>\; ;\; </mo>}\bar{\rho }=\frac{\widehat{\stackrel{¯}{\rho }}}{{\rho }_0}\mathrm{<mo>\; ;\; </mo>}{}_d\bar{\sigma }{}^{\left(0\right)}=\frac{{}_d\widehat{\stackrel{¯}{\sigma }}{}^{\left(0\right)}}{{\tau }_0};\bar{p}=\frac{\widehat{\stackrel{¯}{p}}}{p_0}\\ \eta =\frac{\hat{\eta }}{{\eta }_0}\mathrm{<mo>\; ;\; </mo>}{\eta }_1=\frac{{\hat{\eta }}_1}{{\eta }_0}\mathrm{<mo>\; ;\; </mo>}{\eta }_3=\frac{{\hat{\eta }}_3}{{\eta }_0};t_0=\frac{L_0}{v_0}\mathrm{<mo>\; ;\; </mo>}\bar{F}=\frac{{\bar{F}}^b}{F_0}\end{array}$ (73)

Using Equation (73) in Equations (70)-(72), we can obtain the following dimensionless form for the conservation of mass, balance of linear momentum, and the constitutive theories.

$\frac{\partial \bar{\rho }}{\partial t}+\bar{\rho }\left(\bar{\nabla }\cdot \bar{v}\right)=0$ (74)

$\bar{\rho }\frac{D\bar{v}}{Dt}+\left(\frac{L_0{F}_0}{v_0^2}\right)\bar{\rho }{\bar{F}}^b+\left(\frac{p_0}{{\rho }_0{v}_0^2}\right)\bar{\nabla }\cdot \bar{p}-\left(\frac{{\tau }_0}{{\rho }_0{v}_0^2}\right){\left({}_d\bar{\sigma }{}^{\left(0\right)}\right)}^{\text{T}}\cdot \bar{\nabla }=0$ (75)

$\begin{array}{c}{}_d\bar{\sigma }{}^{\left(0\right)}+De\left({}_d\bar{\sigma }{}^{\left(1\right)}\right)=\left(\frac{{\eta }_0{v}_0}{{\tau }_0{L}_0}\right)2\eta \bar{D}+\left(\frac{{\eta }_0}{L_0{v}_0{\rho }_0}\right)2\eta D{e}_2\left({\gamma }_{\left(2\right)}\right)\\ +\left(\frac{v_0}{{\eta }_0{L}_0}\right)\frac{De}{\eta }\alpha {\left({}_d\bar{\sigma }{}^{\left(0\right)}\right)}^2+\left(\frac{{\eta }_0{v}_0^2}{{\tau }_0{L}_0^2}\right){\eta }_1{\left(\bar{D}\right)}^2\\ +\left(\frac{{\eta }_0{v}_0^2}{{\tau }_0{L}_0^2}\right){\eta }_3\text{tr}\left({\bar{D}}^2\right)I\end{array}$ (76)

If we define

$\begin{array}{l}Re=\frac{{\rho }_0{v}_0{L}_0}{{\eta }_0}\mathrm{<mo>\; ;\; </mo>}De=\frac{{\hat{\lambda }}_1{v}_0}{L_0}\mathrm{<mo>\; ;\; </mo>}D{e}_2=\frac{{\hat{\lambda }}_2{v}_0}{L_0}\mathrm{<mo>\; ;\; </mo>}{\eta }_{10}={\eta }_1\left(\frac{v_0}{L_0}\right)\mathrm{<mo>\; ;\; </mo>}{\eta }_{30}={\eta }_3\left(\frac{v_0}{L_0}\right)\\ p_0={\tau }_0={\rho }_0{v}_0^2\left(\text{characteristickineticenergy}\left(\text{CKE}\right)\right)\end{array}$ (77)

Then using Equation (77), we can write Equations (74)-(76) as

$\frac{\partial \bar{\rho }}{\partial t}+\bar{\rho }\left(\bar{\nabla }\cdot \bar{v}\right)=0$ (78)

${\rho }_0\frac{D\bar{v}}{Dt}+\left(\frac{L_0{F}_0}{v_0^2}\right)\bar{\rho }{\bar{F}}^b+\bar{\nabla }\cdot \bar{p}-{\left({}_d\bar{\sigma }{}^{\left(0\right)}\right)}^{\text{T}}\cdot \bar{\nabla }=0$ (79)

$\begin{array}{c}{}_d\bar{\sigma }{}^{\left(0\right)}+De\left({}_d\bar{\sigma }{}^{\left(1\right)}\right)=\frac{2\eta }{Re}\bar{D}+\frac{2\eta D{e}_2}{Re}\left({\gamma }_{\left(2\right)}\right)+\left(\frac{v_0}{{\eta }_0{L}_0}\right)\frac{De}{\eta }\alpha {\left({}_d\bar{\sigma }{}^{\left(0\right)}\right)}^2\\ +\frac{{\eta }_{10}}{Re}{\bar{D}}^2+\frac{{\eta }_{30}}{Re}\text{tr}\left({\bar{D}}^2\right)I\end{array}$ (80)

Equations (78)-(80) constitute the complete enhanced mathematical model for incompressible, isothermal flow of polymeric fluids in ${ℝ}^3$ and are used in the model problem studies. We refer to this model as model B in the model problem studies.

5. Model Problems

In this section we consider two boundary value problems: fully developed flow between parallel plates and a square lid-driven cavity. Solutions are presented for Maxwell, Oldroyd-B, and Giesekus fluids using the new, enhanced constitutive theories presented in this paper. Results obtained using the new constitutive theory are compared with the constitutive models used currently to demonstrate the benefits of using enhanced constitutive theories for the deviatoric Cauchy stress tensor derived in this paper. In both model problems, we use the following material coefficients and reference values.

Maxwell and Oldroyd-B Fluids:

$L_0=0.015\text{m},v_0=0.015325\text{m}/\text{s},{\rho }_0=\hat{\rho }=998.2\text{kg}/{\text{m}}^{\text{3}},$

$\begin{array}{l}{\eta }_0=\hat{\eta }=1.002\times {10}^{-3}\text{Pa}\cdot \text{s},{\hat{\eta }}_s=9.018\times {10}^{-4}\text{Pa}\cdot \text{s}\mathrm{,}\\ {\hat{\eta }}_p=1.002\times {10}^{-4}\text{Pa}\cdot \text{s}\mathrm{,}{\hat{\lambda }}_1=0.1\text{s}\mathrm{,}{\hat{\lambda }}_2=0.05\text{s}\left(0\text{for}\text{Maxwell}\right)\mathrm{,}\\ \alpha =0,t_0=\frac{L_0}{v_0}=0.97879\text{s}\mathrm{,}Re=\mathrm{229,}De=\mathrm{0.10217,}\\ D{e}_2\mathrm{<mo>\; =\; </mo>0.051085}\left(0\text{for}\text{Maxwell}\right)\mathrm{,}{p}_0={\tau }_0={\rho }_0{v}_0^2(\; CKE\; )\end{array}$

(Parallel plates)

${\hat{\eta }}_1={\hat{\eta }}_3$ are 0%, 5%, 10%, 15%, and 20% of $\hat{\eta }={\hat{\eta }}_s+{\hat{\eta }}_p=1.002\times {10}^{-3}\text{Pa}\cdot \text{s}$ (0.0, 0.501, 1.002, 1.503, 2.004) × 10⁻⁴ Pa·s

(Lid-driven cavity)

${\hat{\eta }}_1={\hat{\eta }}_3$ are 0%, 0.2%, 1%, 5%, and 20% of $\hat{\eta }={\hat{\eta }}_s+{\hat{\eta }}_p=1.002\times {10}^{-3}\text{Pa}\cdot \text{s}$ (0.0, 2.004, 10.02, 50.1, 200.4) × 10⁻⁶ Pa·s

Giesekus Fluid:

$\begin{array}{l}{\rho }_0=\hat{\rho }=800\text{kg}/{\text{m}}^{\text{3}},{\eta }_0=\hat{\eta }=1.426\text{Pa}\cdot \text{s},{\hat{\eta }}_s=0.002\text{Pa}\cdot \text{s}\mathrm{,}\\ {\hat{\eta }}_p=1.424\text{Pa}\cdot \text{s}\mathrm{,}{\hat{\lambda }}_1=0.06\text{s}\mathrm{,}{\hat{\lambda }}_2=\mathrm{0,}\alpha =0.15\end{array}$

(Parallel plates) $L_0=0.015\text{m}$, $v_0=0.015325\text{m}/\text{s}$

${\hat{\eta }}_1={\hat{\eta }}_3$ are 0%, 5%, 10%, 15%, and 20% of $\hat{\eta }={\hat{\eta }}_s+{\hat{\eta }}_p=1.426\text{Pa}\cdot \text{s}$ (0.0, 0.0713, 0.1426, 0.2139, 0.2852) Pa·s

$t_0=\frac{L_0}{v_0}=0.97879\text{s}\mathrm{,}Re=\mathrm{1.2896,}De=\mathrm{0.613,}{p}_0={\tau }_0={\rho }_0{v}_0^2(\; CKE\; )$

(Lid-driven cavity) $L_0=0.1\text{m}$, $v_0=0.025\text{m}/\text{s}$

${\hat{\eta }}_1={\hat{\eta }}_3$ are 0%, 2.5%, 5%, 7.5%, and 10% of $\hat{\eta }={\hat{\eta }}_s+{\hat{\eta }}_p=1.426\text{Pa}\cdot \text{s}$ (0.0, 0.03565, 0.0713, 0.10695, 0.1426) Pa·s

$t_0=\frac{L_0}{v_0}=4\text{s},Re=1.403,De=0.025,p_0={\tau }_0={\rho }_0{v}_0^2(\; CKE\; )$

5.1. Mathematical Model in ${ℝ}^2$

The expanded forms of the conservation of mass, balance of linear momentum, and the constitutive theories in ${ℝ}^2$ for boundary value problems are given by (in the absence of body forces)

$\bar{\rho }\left(\frac{\partial \bar{u}}{\partial \bar{x}}+\frac{\partial \bar{v}}{\partial \bar{y}}\right)=0$ (81)

$\bar{u}\frac{\partial \bar{u}}{\partial \bar{x}}+\bar{v}\frac{\partial \bar{u}}{\partial \bar{y}}+\frac{\partial \bar{p}}{\partial \bar{x}}-\frac{\partial \left({}_d{\bar{\sigma }}_{11}^{\left(0\right)}\right)}{\partial \bar{y}}-\frac{\partial \left({}_d{\bar{\sigma }}_{12}^{\left(0\right)}\right)}{\partial \bar{y}}=0$ (82)

$\bar{u}\frac{\partial \bar{v}}{\partial \bar{x}}+\bar{v}\frac{\partial \bar{v}}{\partial \bar{y}}+\frac{\partial \bar{p}}{\partial \bar{y}}-\frac{\partial \left({}_d{\bar{\sigma }}_{12}^{\left(0\right)}\right)}{\partial \bar{y}}-\frac{\partial \left({}_d{\bar{\sigma }}_{22}^{\left(0\right)}\right)}{\partial \bar{y}}=0$ (83)

Explicit expressions for three constitutive equations for ${}_d{\bar{\sigma }}_{11}^{\left(0\right)}$, ${}_d{\bar{\sigma }}_{22}^{\left(0\right)}$, and ${}_d\bar{\sigma }{}_{12}^{\left(0\right)}={}_d\bar{\sigma }{}_{21}^{\left(0\right)}$ can be obtained using Equation (80) and the following.

$\left[{}_d{\bar{\sigma }}^{\left(0\right)}\right]=\left[\begin{array}{cc}{}_d{\bar{\sigma }}_{11}^{\left(0\right)}& {}_d{\bar{\sigma }}_{12}^{\left(0\right)}\\ {}_d{\bar{\sigma }}_{12}^{\left(0\right)}& {}_d{\bar{\sigma }}_{22}^{\left(0\right)}\end{array}\right];\left[\bar{L}\right]=\left[\begin{array}{cc}\frac{\partial \bar{u}}{\partial \bar{x}}& \frac{\partial \bar{u}}{\partial \bar{y}}\\ \frac{\partial \bar{v}}{\partial \bar{x}}& \frac{\partial \bar{v}}{\partial \bar{y}}\end{array}\right]$ (84)

$\left[{}_d{\bar{\sigma }}^{\left(1\right)}\right]=\frac{D\left[{}_d{\bar{\sigma }}^{\left(0\right)}\right]}{Dt}+\left[\bar{L}\right]\left[{}_d{\bar{\sigma }}^{\left(0\right)}\right]+\left[{}_d{\bar{\sigma }}^{\left(0\right)}\right]{\left[\bar{L}\right]}^{\text{T}}$ (85)

$\frac{D\left[{}_d{\bar{\sigma }}^{\left(0\right)}\right]}{Dt}=\bar{u}\frac{\partial \left[{}_d{\bar{\sigma }}^{\left(0\right)}\right]}{\partial \bar{x}}+\bar{v}\frac{\partial \left[{}_d{\bar{\sigma }}^{\left(0\right)}\right]}{\partial \bar{y}}$ (86)

$\left[\bar{D}\right]=\left[{\gamma }_{\left(1\right)}\right]=\frac{1}{2}\left(\left[\bar{L}\right]+{\left[\bar{L}\right]}^{\text{T}}\right)=\left[\begin{array}{cc}\frac{\partial \bar{u}}{\partial \bar{x}}& \frac{1}{2}\left(\frac{\partial \bar{u}}{\partial \bar{y}}+\frac{\partial \bar{v}}{\partial \bar{x}}\right)\\ \frac{1}{2}\left(\frac{\partial \bar{u}}{\partial \bar{y}}+\frac{\partial \bar{v}}{\partial \bar{x}}\right)& \frac{\partial \bar{v}}{\partial \bar{y}}\end{array}\right]$ (87)

$\left[{\gamma }_{\left(2\right)}\right]=\left(\bar{u}\frac{\partial }{\partial \bar{x}}+\bar{v}\frac{\partial }{\partial \bar{y}}\right)\left[{\gamma }_{\left(1\right)}\right]+{\left[\bar{L}\right]}^{\text{T}}\left[{\gamma }_{\left(1\right)}\right]+\left[{\gamma }_{\left(1\right)}\right]\left[\bar{L}\right]$ (88)

This mathematical model is used to present numerical studies for a lid-driven cavity.

Solutions of the model problems are obtained using finite element methods based on the residual functional (least squares finite element method) [22] in which the local approximations are in $H^{k\mathrm{,}p}\left({\bar{\Omega }}^e\right)$ higher order scalar product spaces permitting a higher degree of local approximation as well as desired higher order global differentiability [23].

5.2. Mathematical Model in ${ℝ}^1$

We consider fully developed flow between parallel plates. If $\bar{x}$ is the direction of flow, then for fully developed flow, the flow is independent of the $\bar{x}$ coordinate, hence the gradients of the dependent variables in the $\bar{x}$ direction are zero. The complete mathematical model in ${ℝ}^1$ can be obtained using Equations (78)-(80). We can write the following (neglecting body forces) for the balance of linear momentum and the constitutive theories.

$\begin{array}{l}\frac{\partial \bar{p}}{\partial \bar{x}}-\frac{\partial \left({}_d{\bar{\sigma }}_{12}^{\left(0\right)}\right)}{\partial y}=0\\ \frac{\partial \bar{p}}{\partial \bar{y}}-\frac{\partial \left({}_d{\bar{\sigma }}_{22}^{\left(0\right)}\right)}{\partial y}=0\end{array}$ (89)

and the constitutive equations are given by

$\begin{array}{l}{}_d{\bar{\sigma }}_{11}^{\left(0\right)}+De\frac{\partial \bar{u}}{\partial \bar{y}}\left(2\left({}_d{\bar{\sigma }}_{12}^{\left(0\right)}\right)\right)=\left(\frac{v_0}{{\eta }_0{L}_0}\right)\frac{De}{\eta }\alpha \left({\left({}_d{\bar{\sigma }}_{11}^{\left(0\right)}\right)}^2+{\left({}_d{\bar{\sigma }}_{12}^{\left(0\right)}\right)}^2\right)\\ +\frac{1}{4}\frac{{\eta }_{10}}{Re}{\left(\frac{\partial \bar{u}}{\partial \bar{y}}\right)}^2+\frac{1}{2}\frac{{\eta }_{30}}{Re}{\left(\frac{\partial \bar{u}}{\partial \bar{y}}\right)}^2\end{array}$

$\begin{array}{l}{}_d{\bar{\sigma }}_{12}^{\left(0\right)}+De\frac{\partial \bar{u}}{\partial \bar{y}}\left({}_d{\bar{\sigma }}_{22}^{\left(0\right)}\right)=\frac{\eta }{Re}\left(\frac{\partial \bar{u}}{\partial \bar{y}}\right)\\ +\left(\frac{v_0}{{\eta }_0{L}_0}\right)\frac{De}{\eta }\alpha ({}_d\bar{\sigma }{}_{11}^{\left(0\right)}{}_d\bar{\sigma }{}_{12}^{\left(0\right)}+{}_d\bar{\sigma }{}_{12}^{\left(0\right)}{}_d\bar{\sigma }{}_{22}^{(\; 0\; )})\end{array}$

$\begin{array}{l}{}_d{\bar{\sigma }}_{22}^{\left(0\right)}=2\frac{D{e}_2\eta }{Re}{\left(\frac{\partial \bar{u}}{\partial \bar{y}}\right)}^2+\left(\frac{v_0}{{\eta }_0{L}_0}\right)\frac{De}{\eta }\alpha \left({\left({}_d{\bar{\sigma }}_{12}^{\left(0\right)}\right)}^2+{\left({}_d{\bar{\sigma }}_{22}^{\left(0\right)}\right)}^2\right)\\ +\frac{1}{4}\frac{{\eta }_{10}}{Re}{\left(\frac{\partial \bar{u}}{\partial \bar{y}}\right)}^2+\frac{1}{2}\frac{{\eta }_{30}}{Re}{\left(\frac{\partial \bar{u}}{\partial \bar{y}}\right)}^2\end{array}$ (90)

6. Model Problem Studies

In this section, we present converged finite element solutions of fully developed flow between parallel plates and lid-driven cavities for currently used Maxwell, Oldroyd-B and Giesekus models as well as using the enhanced constitutive model presented in this paper. In the finite element method based on the residual functional (least squares method) used here, when the approximation spaces are minimally conforming (or of orders higher than minimally conforming), the proximity of the L₂-norm of the residual functional is an absolute measure of the accuracy and convergence of the computed solutions. In all numerical studies presented here, the L₂-norm of the residuals of $O\left({10}^{-4}\right)$ or lower is achieved, ensuring convergence of the computed solutions to the true solution of the BVP. Since the mathematical model consists of nonlinear partial differential equations, the solution of the nonlinear algebraic system of equations resulting from the residual functional formulation is obtained using Newton’s linear method with line search described in the following.

6.1. Solution Procedure for Nonlinear Boundary Value Problems

An unconditionally stable (Variationally Consistent (VC) [22]) finite element formulation of nonlinear BVPs can be constructed using the residual functional. For simplicity, we illustrate the details for a single nonlinear differential equation describing the BVP. Let

$A\varphi -f=0\forall \bar{x}\in \bar{\Omega }\subset ℝ$

be a BVP. Let ${\bar{\Omega }}^T={\displaystyle \underset{e}{\cup }}{\bar{\Omega }}^e$ be the discretization of $\bar{\Omega }$ in which ${\bar{\Omega }}^e={\Omega }^e\cup {\Gamma }^e$ is a finite element e with a closed boundary ${\Gamma }^e$. Let ${\varphi }_h^e$ be the approximation of $\varphi$ over ${\bar{\Omega }}^e$ (local approximation) and ${\varphi }_h$

${\varphi }_h={\displaystyle \underset{e}{\cup }}{\varphi }_h^e$ (91)

be the approximation of $\varphi$ over ${\bar{\Omega }}^T$. Then the residual function E is defined by

$E=A{\varphi }_h-f\forall x\in {\bar{\Omega }}^T$ (92)

The residual functional $I\left({\varphi }_h\right)$ can be written as

$\begin{array}{l}I\left({\varphi }_h\right)={\left(E\left({\varphi }_h\right)\mathrm{,}E\left({\varphi }_h\right)\right)}_{{\bar{\Omega }}^T}\mathrm{<mo>\; =\; </mo>}{\displaystyle \underset{e}{\sum }}{\left(E^e\left({\varphi }_h^e\right)\mathrm{,}{E}^e\left({\varphi }_h^e\right)\right)}_{{\bar{\Omega }}^e}={\displaystyle \underset{e}{\sum }}{I}^e\left({\varphi }_h^e\right)\\ \text{In}\text{which}{E}^e=A{\varphi }_h^e-f\forall \bar{x}\in {\bar{\Omega }}^e\end{array}$ (93)

If $I\left({\varphi }_h\right)$ is differentiable in its arguments (i.e. ${\varphi }_h$ ), then $\delta I\left({\varphi }_h\right)=0$ is a necessary condition for an extremum of $I\left({\varphi }_h\right)$.

$\delta I\left({\varphi }_h\right)=2{\left(E,\delta E\right)}_{{\bar{\Omega }}^T}={\displaystyle \sum }\delta I^e\left({\varphi }_h^e\right)=2{\displaystyle \underset{e}{\sum }}\left(E^e,\delta E^e\right)=2{\displaystyle \underset{e}{\sum }}{g}^e=2g=0$ (94)

From Equation (94), we can confirm that Euler’s equation from $\delta I\left({\varphi }_h\right)=0$ is in fact the BVP. Thus, a function of ${\varphi }_h$ that yields the extremum of $I\left({\varphi }_h\right)$ is also a solution to the BVP.

When the differential operator is nonlinear, then g in Equation (94) is a nonlinear function of ${\varphi }_h$. We must find a solution ${\varphi }_h$ iteratively that satisfies Equation (94). This is accomplished using Newton’s linear method with line search [22]. The final result is that if ${\varphi }_h^0$ is the assumed starting solution then the improved solution ${\varphi }_h$ is given by

${\varphi }_h={\varphi }_h^0+\alpha \Delta {\varphi }_h$ (95)

$\Delta {\varphi }_h=-\frac{1}{2}{\left({\delta }^{2}I\left({\varphi }_h\right)\right)}_{{\varphi }_h^0}^{-1}{g}_{{\varphi }_h^0}$ (96)

in which

${\delta }^{2}I\left({\varphi }_h\right)\simeq 2\left(\delta E\mathrm{,}\delta E\right)=2{\displaystyle \underset{e}{\sum }}\left(\delta E^e,\delta E^e\right)>0$ (97)

$0<\alpha <2\text{is}\text{such}\text{that}I\left({\varphi }_h\right)\le I\left({\varphi }_h^0\right)$ (98)

and ${\varphi }_h$ is considered to be converged when

$\underset{i}{\max }\left|g_i\left({\varphi }_h\right)\right|\le \Delta$ (99)

where Δ is a preset tolerance for computed zero (generally O(10⁻⁶) or lower). If the tolerance in Equation (99) is not satisfied, ${\varphi }_h^0$ is set to ${\varphi }_h$ and Equations (95)-(99) are repeated until Equation (99) is satisfied. Since Newton’s linear method has quadratic convergence, an accuracy of Δ = O(10⁻⁶) is generally achieved in less than five iterations.

This approach for one differential equation can be easily extended to m differential equations. Let $E_i;i=1,2,\cdots ,m$ be the residual functions resulting from each partial differential equation. Then the residual functional I for ${\bar{\Omega }}^T$ can be written as

$I={\displaystyle \underset{i=1}{\overset{m}{\sum }}}{\left(E_i,E_i\right)}_{{\bar{\Omega }}^T}={\displaystyle \underset{i=1}{\overset{m}{\sum }}}{\displaystyle \underset{e}{\sum }}{\left(E_i^e,E_i^e\right)}_{{\bar{\Omega }}^e}={\displaystyle \underset{i=1}{\overset{m}{\sum }}}{\displaystyle \underset{e}{\sum }}{I}_i^e$ (100)

Remaining details follow the details given above for one equation.

6.2. Model Problem I: Fully Developed Flow between Parallel Plates

We consider fully developed flow between parallel plates separated by a dimensionless distance of $H=2$ ( $\hat{H}=3$ cm). Figure 1 shows a schematic, discretization, and boundary conditions.

Figure 1. Fully developed flow between parallel plates.

The origin is located halfway between the centers of the plates. The upper half of the domain ( $0\le \bar{x}\le 1$ ) is discretized using ten three node p-version hierarchical elements. Since the mathematical models for all three fluids consist of a system of first order partial differential equations, the local approximation of class $C^1\left({\bar{\Omega }}^e\right)$ ensure integrals are in the Riemann sense for the discretizations ${\bar{\Omega }}^T$. Initial studies show that a p-level of 5 (p= 5) is sufficient to yield a residual functional I for ${\bar{\Omega }}^T$ of O(10⁻⁶) or lower. Newton’s linear method with line search converges in approximately 2 to 4 iterations with a tolerance of $\Delta ={10}^{-6}$ on $\underset{i}{\max }\left|g_i\right|$. We choose ${\hat{\eta }}_1={\hat{\eta }}_3=0\%$, 5%, 10%, 15%, and 20% of ${\eta }_0$ giving ${\hat{\eta }}_1={\hat{\eta }}_3=0.0$, 0.0000501, 0.0001002, 0.0001503, and 0.0002004 for Maxwell and Oldroyd-B fluids and ${\hat{\eta }}_1={\hat{\eta }}_3=0.0$, 0.0713, 0.1426, 0.2139, and 0.2852 for Giesekus fluid. The flow is pressure driven with $\partial \bar{p}/\partial \bar{x}=-0.0133$ for Maxwell and Oldroyd-B fluids and $\partial \bar{p}/\partial \bar{x}=-0.2$ for Giesekus fluid. Figures 2(a)-(d) show plots of velocity $\bar{u}$, stresses ${}_d{\bar{\sigma }}_{12}^{\left(0\right)}$, ${}_d{\bar{\sigma }}_{11}^{\left(0\right)}$, and ${}_d{\bar{\sigma }}_{22}^{\left(0\right)}$ versus $\bar{y}$ for different values of ${\hat{\eta }}_1={\hat{\eta }}_3$ for Maxwell fluid. Velocity $\bar{u}$ reduces with increasing values of ${\hat{\eta }}_1={\hat{\eta }}_3$. Shear stress ${}_d{\bar{\sigma }}_{12}^{\left(0\right)}$ remains unaffected as it only depends on $\partial \bar{p}/\partial \bar{x}$. ${}_d{\bar{\sigma }}_{22}^{\left(0\right)}$ is zero when ${\hat{\eta }}_1={\hat{\eta }}_3=0$ (standard Maxwell model), but progressively increasing values of ${\hat{\eta }}_1={\hat{\eta }}_3$ yield progressively increasing values of ${}_d{\bar{\sigma }}_{22}^{\left(0\right)}$ for $0\le \bar{y}\le 1$. ${}_d{\bar{\sigma }}_{11}^{\left(0\right)}$ has nonzero values for ${\hat{\eta }}_1={\hat{\eta }}_3=0$ (as expected). The ${}_d{\bar{\sigma }}_{11}^{\left(0\right)}$ values also increase along $0\le \bar{y}\le 1$ for progressively increasing values of ${\hat{\eta }}_1={\hat{\eta }}_3$.

Figures 3(a)-(d) show plots of $\bar{u}$, ${}_d{\bar{\sigma }}_{12}^{\left(0\right)}$, ${}_d{\bar{\sigma }}_{11}^{\left(0\right)}$, and ${}_d{\bar{\sigma }}_{22}^{\left(0\right)}$ versus $\bar{y}$ for different values of ${\hat{\eta }}_1={\hat{\eta }}_3$ (same as those used for Maxwell model) for Oldroyd-B model. Velocity $\bar{u}$ versus $\bar{y}$ and ${}_d{\bar{\sigma }}_{12}^{\left(0\right)}$ versus $\bar{y}$ plots are almost identical to those of the Maxwell model. ${}_d{\bar{\sigma }}_{22}^{\left(0\right)}$ versus $\bar{y}$ in Figure 3(d) is exactly identical to that of the Maxwell model. This is expected as ${}_d{\bar{\sigma }}_{22}^{\left(0\right)}$ is zero for Maxwell as well as Oldroyd-B models when ${\hat{\eta }}_1={\hat{\eta }}_3=0$, thus ${}_d{\bar{\sigma }}_{22}^{\left(0\right)}$ in both models is only due to ${\hat{\eta }}_1={\hat{\eta }}_3\ne 0$ and has the same mechanism. ${}_d{\bar{\sigma }}_{11}^{\left(0\right)}$ in the Oldroyd-B model is lower than that of the Maxwell model for all values of ${\hat{\eta }}_1={\hat{\eta }}_3$ due to additional dissipation.

Figure 2. Fully developed flow between parallel plates: Maxwell.

Figure 3. Fully developed flow between parallel plates: Oldroyd-B.

Figure 4. Fully developed flow between parallel plates: Giesekus.

The currently used Giesekus model naturally produces ${}_d{\bar{\sigma }}_{22}^{\left(0\right)}\ne 0$ when ${\hat{\eta }}_1={\hat{\eta }}_3=0$ (Figure 4(d)), but ${}_d{\bar{\sigma }}_{22}^{\left(0\right)}$ values are enhanced (increased) for progressively increasing ${\hat{\eta }}_1={\hat{\eta }}_3$. ${}_d{\bar{\sigma }}_{11}^{\left(0\right)}$ (Figure 4(c)) in this model is an order of magnitude higher than Maxwell or Oldroyd-B models, thus we do not observe an appreciable change in ${}_d{\bar{\sigma }}_{22}^{\left(0\right)}\ne 0$ versus $\bar{y}$ for increasing values of ${\hat{\eta }}_1={\hat{\eta }}_3$. We note that for an order of magnitude higher pressure gradient in the Giesekus model as compared to Maxwell and Oldroyd-B models, the velocity $\bar{u}$ is an order of magnitude lower. An exploded plot of $\bar{u}$ versus $\bar{y}$ in Figure 4(a) shows progressively decreasing velocity $\bar{u}$ versus $\bar{y}$ for progressively increasing ${\hat{\eta }}_1={\hat{\eta }}_3$.

6.3. Model Problem II: Lid-Driven Square Cavity

Figure 5(a) shows a schematic of the square 1 × 1 (dimensionless) lid-driven cavity with the lid moving at velocity $\bar{u}=1.0$. Boundary conditions are also shown in Figure 5(a). Figure 5(b) shows a graded discretization of the cavity using 36 p-version hierarchical elements with a higher order global differentiability local approximation. ${\hat{\eta }}_1={\hat{\eta }}_3=0\%$, 0.2%, 1%, 5%, and 20% of ${\eta }_0$ giving ${\hat{\eta }}_1={\hat{\eta }}_3=0.0$, 0.000002004, 0.00001002, 0.0000501, and 0.0002004 for Maxwell and Oldroyd-B fluids and ${\hat{\eta }}_1={\hat{\eta }}_3=0\%$, 2.5%, 5%, 7.5%, and 10% of ${\eta }_0$ giving ${\hat{\eta }}_1={\hat{\eta }}_3=0.0$, 0.03565, 0.0713, 0.10695, and 0.1426 for Giesekus fluid are used in the calculations. For solutions of class $C^{11}\left({\bar{\Omega }}^e\right)$ with $p_{\zeta }=p_{\eta }=5$, Newton’s linear method with line search yielded a residual functional I for ${\bar{\Omega }}^T$ of at most O(10⁻⁴) within 12 iterations for these values of ${\hat{\eta }}_1={\hat{\eta }}_3$.

We only present results at $\bar{x}=0.5$ as a function of $\bar{y}$ (vertical center line) for the sake of brevity. Figures 6(a)-(f) show plots of velocities $\bar{u}$, $\bar{v}$, stresses ${}_d{\bar{\sigma }}_{11}^{\left(0\right)}$, ${}_d{\bar{\sigma }}_{22}^{\left(0\right)}$, ${}_d{\bar{\sigma }}_{12}^{\left(0\right)}$, and pressure $\bar{p}$ as a function of $\bar{y}$ at $\bar{x}=0.5$. Velocity $\bar{v}$, stress ${}_d{\bar{\sigma }}_{22}^{\left(0\right)}$, and pressure $\bar{p}$ show the most dependence on progressively increasing values of ${\hat{\eta }}_1={\hat{\eta }}_3$. Similar plots for Oldroyd-B and Giesekus fluids are shown in Figures 7(a)-(f) and Figures 8(a)-(f) respectively. Oldroyd-B model results parallel to those of the Maxwell model. In the case of the Giesekus model, we see significant dependence of velocity $\bar{v}$, stresses ${}_d{\bar{\sigma }}_{11}^{\left(0\right)}$, ${}_d{\bar{\sigma }}_{22}^{\left(0\right)}$, and pressure $\bar{p}$ on progressively increasing values of ${\hat{\eta }}_1={\hat{\eta }}_3$.

Remarks

Based on the numerical studies presented for the two model problems, we can make the following remarks.

1) In the model problem study involving fully developed flow between parallel plates, the new constitutive theory for ${}_d\bar{\sigma }{}^{\left(0\right)}$ produces nonzero, progressively increasing ${}_d{\bar{\sigma }}_{22}^{\left(0\right)}$ with progressively increasing ${\hat{\eta }}_1={\hat{\eta }}_3$ values. This is more dramatic for Maxwell and Oldroyd-B models in which currently used constitutive theories yield ${}_d{\bar{\sigma }}_{22}^{\left(0\right)}=0$.

2) In model problem 1 for the currently used Giesekus model, ${}_d{\bar{\sigma }}_{22}^{\left(0\right)}\ne 0$. The new constitutive theory results in additional ${}_d{\bar{\sigma }}_{11}^{\left(0\right)}$ and ${}_d{\bar{\sigma }}_{22}^{\left(0\right)}$ stresses greater than their values for ${\hat{\eta }}_1={\hat{\eta }}_3=0$ from the currently used constitutive model.

Figure 5. Lid drive square cavity.

Figure 6. Square lid-driven cavity: Maxwell.

Figure 7. Square lid-driven cavity: Oldroyd-B.

Figure 8. Square lid-driven cavity: Giesekus.

3) ${}_d{\bar{\sigma }}_{12}^{\left(0\right)}$ remains unaffected in model problem 1 as it can be theoretically determined using $\partial \bar{p}/\partial \bar{x}$ (as is evident from Figures 2(b)-4(b)).

4) The influence of the new constitutive theory is also evident in lid-driven cavity. However, due to complex flow physics, clear and concise observations (similar to model problem 1) are difficult in this model problem.

5) Fully developed flow between parallel plates serves as a good model problem for calibration, i.e. determination of ${\hat{\eta }}_1={\hat{\eta }}_3$. Measurements of ${}_d{\bar{\sigma }}_{22}^{\left(0\right)}$ for a given $\partial \bar{p}/\partial \bar{x}$ can be simulated numerically using varied values of ${\hat{\eta }}_1={\hat{\eta }}_3$. Values of ${\hat{\eta }}_1={\hat{\eta }}_3$ yielding the same ${}_d{\bar{\sigma }}_{22}^{\left(0\right)}$ as in the experiment are the correctly calibrated values of ${\hat{\eta }}_1={\hat{\eta }}_3$.

7. Summary and Conclusions

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

1) A new, enhanced constitutive theory has been presented for Maxwell, Oldroyd-B, and Giesekus constitutive models for incompressible polymeric fluids. This constitutive theory is designed to correct the major deficiency in Maxwell and Oldroyd-B models of zero normal stress perpendicular to the direction of the flow. These models produce ${}_d{\bar{\sigma }}_{22}^{\left(0\right)}=0$ when $\bar{x}$ (direction 1) is the direction of the flow as shown in the case of fully developed flow between parallel plates.

2) Constitutive theory derivations are initiated using entropy inequality and a general constitutive theory of orders $\left(m\mathrm{,}n\right)$ (in stress and strain rates) is derived based on the integrity (complete basis) and representation theorem.

3) It is shown that the currently used constitutive models for polymeric fluids are a small subset of the general constitutive theory of orders $\left(m\mathrm{,}n\right)$ based on the integrity presented in this paper.

4) The enhancement of the currently used constitutive models is accomplished by retaining additional generators and invariant(s) in the constitutive theory from the integrity. The rationale is presented for retaining the additional generator ${\left[\bar{D}\right]}^2$ and additional invariant $\text{tr}\left({\left[\bar{D}\right]}^2\right)$ in the enhanced constitutive theory.

5) The new, enhanced constitutive theory yields additional nonzero ${}_d\bar{\sigma }{}_{11}^{\left(0\right)}={}_d\bar{\sigma }{}_{22}^{\left(0\right)}$, but shear stress ${}_d{\bar{\sigma }}_{12}^{\left(0\right)}$ remains unaffected. The magnitude of ${}_d{\bar{\sigma }}_{11}^{\left(0\right)}$ and ${}_d{\bar{\sigma }}_{22}^{\left(0\right)}$ depend upon the material coefficients ${\hat{\eta }}_1$ and ${\hat{\eta }}_3$. In the present work, a simple case of ${\hat{\eta }}_1={\hat{\eta }}_3$ is considered. However, ${\hat{\eta }}_1$ and ${\hat{\eta }}_3$ can be two additional material coefficients that can be determined experimentally.

6) We remark that the generator ${\left[\bar{D}\right]}^2$ and invariant $\text{tr}\left({\left[\bar{D}\right]}^2\right)$ are part of the integrity, hence will always be present in a constitutive theory for ${}_d\bar{\sigma }{}^{\left(0\right)}$ if it would have been based on the complete basis or integrity. Extremely simplified forms of currently used Maxwell, Oldroyd-B, and Giesekus models that do not include enough terms from the integrity create this deficiency of normal stress perpendicular to the flow direction being zero.

7) In this paper, we have shown that the inclusions of ${\left[\bar{D}\right]}^2$ and $\text{tr}\left({\left[\bar{D}\right]}^2\right)$ from the integrity help us in restoring nonzero normal stress perpendicular to the flow direction.

8) In the first model problem (parallel plates), we clearly show that nonzero, progressively increasing values of ${}_d{\bar{\sigma }}_{22}^{\left(0\right)}$ are obtained for progressively increasing values of ${\hat{\eta }}_1={\hat{\eta }}_3$. Velocity $\bar{u}$ and stress ${}_d{\bar{\sigma }}_{11}^{\left(0\right)}$ change accordingly while ${}_d{\bar{\sigma }}_{12}^{\left(0\right)}$ remains unaffected as it only depends upon $\partial \bar{p}/\partial \bar{x}$.

9) The influence of the new constitutive theory for lid-driven cavities has also been illustrated for progressively increasing values of ${\hat{\eta }}_1={\hat{\eta }}_3$.

10) Fully developed flow between parallel plates can be used to calibrate the model, i.e. determination of ${\hat{\eta }}_1={\hat{\eta }}_3$ as two material coefficients.

Acknowledgements

The first author is grateful for his endowed professorships and the department of mechanical engineering of the University of Kansas for providing financial support to the second author. The computational facilities provided by the Computational Mechanics Laboratory of the mechanical engineering department are also acknowledged.

Conflicts of Interest

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

References

[1]	Maxwell, J.C. (1867) On the Dynamical Theory of Gases. Philosophical Transactions of the Royal Society of London A, 157, 49-88. https://doi.org/10.1098/rstl.1867.0004
[2]	Oldroyd, J.G. (1950) On the Formulation of Rheological Equations of State. Proceedings of the Royal Society of London A, 200, 523-541. https://doi.org/10.1098/rspa.1950.0035
[3]	Bird, R.B., Armstrong, R.C. and Hassager, O. (1987) Dynamics of Polymeric Liquids, Volume 1, Fluid Mechanics. 2nd Edition, John Wiley and Sons, New York.
[4]	Bird, R.B., Armstrong, R.C. and Hassager, O. (1987) Dynamics of Polymeric Liquids, Volume 2, Kinetic Theory. 2nd Edition, John Wiley and Sons, New York.
[5]	Surana, K.S., Nunez, D., Reddy, J.N. and Romkes, A. (2013) Rate Constitutive Theory for Ordered Thermofluids. Continuum Mechanics and Thermodynamics, 25, 626-662. https://doi.org/10.1007/s00161-012-0257-6
[6]	Surana, K.S., Nunez, D., Reddy, J.N. and Romkes, A. (2014) Rate Constitutive Theory for Ordered Thermoviscoelastic Fluids: Polymers. Continuum Mechanics and Thermodynamics, 26, 143-181. https://doi.org/10.1007/s00161-013-0295-8
[7]	Surana, K.S., Nunez, D. and Reddy, J.N. (2013) Giesekus Constitutive Model for Thermoviscoelastic Fluids Based on Ordered Rate Constitutive Theories. Journal of Research Updates in Polymer Science, 2, 232-260. https://doi.org/10.6000/1929-5995.2013.02.04.5
[8]	Surana, K.S. (2015) Advanced Mechanics of Continua. CRC/Taylor and Francis, Boca Raton. https://doi.org/10.1201/b17959
[9]	Surana, K.S. (2021) Classical Continuum Mechanics. 2nd Edition, CRC/Taylor and Francis, Boca Raton. https://doi.org/10.1201/9781003105336
[10]	Wang, C.C. (1969) On Representations for Isotropic Functions, Part I. Archive for Rational Mechanics and Analysis, 33, 249-267. https://doi.org/10.1007/BF00281278
[11]	Wang, C.C. (1969) On Representations for Isotropic Functions, Part II. Archive for Rational Mechanics and Analysis, 33, 268-287. https://doi.org/10.1007/BF00281279
[12]	Wang, C.C. (1970) A New Representation Theorem for Isotropic Functions, Part I and Part II. Archive for Rational Mechanics and Analysis, 36, 166-223. https://doi.org/10.1007/BF00272241
[13]	Wang, C.C. (1971) Corrigendum to “Representations for Isotropic Functions”. Archive for Rational Mechanics and Analysis, 43, 392-395. https://doi.org/10.1007/BF00252004
[14]	Smith, G.F. (1970) On a Fundamental Error in Two Papers of C.C. Wang, “On Representations for Isotropic Functions, Part I and Part II”. Archive for Rational Mechanics and Analysis, 36, 161-165. https://doi.org/10.1007/BF00272240
[15]	Smith, G.F. (1971) On Isotropic Functions of Symmetric Tensors, Skew-Symmetric Tensors and Vectors. International Journal of Engineering Science, 9, 899-916. https://doi.org/10.1016/0020-7225(71)90023-1
[16]	Spencer, A.J.M. and Rivlin, R.S. (1959) The Theory of Matrix Polynomials and Its Application to the Mechanics of Isotropic Continua. Archive for Rational Mechanics and Analysis, 2, 309-336. https://doi.org/10.1007/BF00277933
[17]	Spencer, A.J.M. and Rivlin, R.S. (1960) Further Results in the Theory of Matrix Polynomials. Archive for Rational Mechanics and Analysis, 4, 214-230. https://doi.org/10.1007/BF00281388
[18]	Spencer, A.J.M. (1971) Theory of Invariants. Chapter 3, Academic Press, Cambridge. https://doi.org/10.1016/B978-0-12-240801-4.50008-X
[19]	Boehler, J.P. (1977) On Irreducible Representations for Isotropic Scalar Functions. Journal of Applied Mathematics and Mechanics, 57, 323-327. https://doi.org/10.1002/zamm.19770570608
[20]	Zheng, Q.S. (1993) On the Representations for Isotropic Vector-Valued, Symmetric Tensor-Valued and Skew-Symmetric Tensor-Valued Functions. International Journal of Engineering Science, 31, 1013-1024. https://doi.org/10.1016/0020-7225(93)90109-8
[21]	Zheng, Q.S. (1993) On Transversely Isotropic, Orthotropic and Relatively Isotropic Functions of Symmetric Tensors, Skew-Symmetric Tensors, and Vectors. International Journal of Engineering Science, 31, 1399-1453. https://doi.org/10.1016/0020-7225(93)90005-F
[22]	Surana, K.S. and Reddy, J.N. (2016) The Finite Element Method for Boundary Value Problems: Mathematics and Computations. CRC/Taylor & Francis, Boca Raton. https://doi.org/10.1201/9781315365718
[23]	Surana, K.S., Carranza, C.H. and Mathi, S.S.C. (2021) k-Version of Finite Element Method for BVPs and IVPs. Mathematics, 9, Article No. 1333. https://doi.org/10.3390/math9121333