A generalisation of the integral Maxwell model: the gK-BKZ model—frame invariance and analytical solutions

Abstract

In this work, we present a generalisation of the Maxwell integral model by developing a new relaxation modulus based on the Mittag-Leffler function. The new model provides a better fit to linear experimental data compared to a classical 2-mode Maxwell model. The model is improved to become frame invariant and several analytical solutions are derived to better understand the behaviour of the model.

Semi-analytical solution for the flow between parallel plates

The equations governing the flow of an incompressible gK-BKZ fluid are:

$$\begin{aligned} \nabla \cdot \textbf{u}&=0 \end{aligned}$$

(A1)

$$\begin{aligned} \frac{\partial \textbf{u}}{\partial t}+\nabla \cdot (\textbf{v v})&=-\nabla p+\frac{1}{R e} \nabla ^{2} \textbf{u}+\nabla \cdot \varvec{\Phi }+\frac{1}{F r^{2}} \textbf{g}. \end{aligned}$$

(A2)

In (A2), $\Phi$ is given by

$$\begin{aligned} \varvec{\Phi }=\varvec{\sigma }-\frac{1}{R e} \dot{\gamma }, \end{aligned}$$

(A3)

where $\dot{\gamma }=\nabla \textbf{u}+(\nabla \textbf{u})^{\textrm{T}}$, $Re=\frac{\rho _{0} U L}{\eta _{0}}$ is the Reynolds number, $F r=\frac{U}{\sqrt{L g}}$ is the Froude number, $W i=\lambda _{r e f} \frac{U}{L}$ is the Weissenberg number; U is the mean velocity, L is a characteristic length, $\textbf{g}$ and $\rho _{0}$ are the gravity and density of the fluid, respectively. The zero shear-rate viscosity is given by $\eta _{0}= S_0 \lambda _{r e f}$ with $\lambda _{r e f}$ the mean relaxation time.

The system of governing equation becomes complete by adding Eqs. (37) and (49).

The equations become dimensionless by proceeding as follows. The dimensional variables are rewritten with a horizontal bar over them. Let

$$\begin{aligned} \overline{\textbf{x}}= & {} L \textbf{x}, \quad \bar{t}=\frac{L}{U} t,\\ \overline{\textbf{u}}= & {} U \textbf{u},\quad \bar{p}=\rho _{0} U^{2} p,\quad \overline{\varvec{\sigma }}=\rho _{0} U^{2} \varvec{\sigma },\\ \bar{\lambda }_{k}= & {} \lambda _{r e f} \lambda _{k}, \quad \bar{S}_{0}=\rho _{0} U^{2} S_0. \end{aligned}$$

Equations (49)-(A2) are solved in 2D. Note that $p=$ $p(\textbf{x}, t), \textbf{u}=(u(\textbf{x}, t), v(\textbf{x}, t))^{T}$ and

$$\begin{aligned} \begin{aligned} \varvec{\sigma }(\textbf{x}, t)&=\left[ \begin{array}{ll} \sigma _{x x} &{} \sigma _{x y} \\ \sigma _{x y} &{} \sigma _{y y} \end{array}\right] , \quad \varvec{\Phi }(\textbf{x}, t)=\left[ \begin{array}{cc} \Phi _{x x} &{} \Phi _{x y} \\ \Phi _{x y} &{} \Phi _{y y} \end{array}\right] , \\ \nabla \textbf{u}(\textbf{x}, t)&=\left[ \begin{array}{ll} \displaystyle \frac{\partial u}{\partial x} &{} \displaystyle \frac{\partial v}{\partial x} \\ \\ \displaystyle \frac{\partial u}{\partial y} &{} \displaystyle \frac{\partial v}{\partial y} \end{array}\right] , \quad \textbf{B}_{t^{\prime }}(\textbf{x}, t)=\left[ \begin{array}{ccc} { B_{t^{\prime }_{x x}}} &{} {B_{t^{\prime }_{x y}}} &{} 0 \\ {B_{t^{\prime }_{x y}}} &{} {B_{t^{\prime }_{y y}}} &{} 0 \\ 0 &{} 0 &{} 1 \end{array}\right] . \end{aligned} \end{aligned}$$

Note that to compute the total tensor in 2D, we consider the three-dimensional Finger tensor. This definition takes into account the absence of deformation in the third coordinate ($B_{t^{\prime }_{z z}}$ in three-dimensional flows) and also guarantees the nonzero determinant property for the Finger tensor.

Considering two-dimensional Cartesian flows, the governing Eqs. (A1)-(A3) are discretized as follows [44, 45]

$$\begin{aligned} \begin{aligned}&\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y} =0, \\&\frac{\partial u}{\partial t}+\frac{\partial \left( u^{2}\right) }{\partial x}+\frac{\partial (u v)}{\partial y}=-\frac{\partial p}{\partial x}+\frac{1}{{\text {Re}}}\left[ \frac{\partial ^{2} u}{\partial x^{2}}+\frac{\partial ^{2} u}{\partial y^{2}}\right] \\&\quad +\frac{\partial \Phi _{x x}}{\partial x}+\frac{\partial \Phi _{x y}}{\partial y}+\frac{1}{F r^{2}} g_{x}, \\&\frac{\partial v}{\partial t}+\frac{\partial (u v)}{\partial x}+\frac{\partial \left( v^{2}\right) }{\partial y} =-\frac{\partial p}{\partial y}+\frac{1}{R e}\left[ \frac{\partial ^{2} v}{\partial x^{2}}+\frac{\partial ^{2} v}{\partial y^{2}}\right] \\&\quad +\frac{\partial \Phi _{x y}}{\partial x}+\frac{\partial \Phi _{y y}}{\partial y}+\frac{1}{F r^{2}} g_{y}, \end{aligned} \end{aligned}$$

where

$$\begin{aligned} \begin{aligned} \Phi _{x x}&=\sigma _{x x}-\frac{2}{R e}\left( \frac{\partial u}{\partial x}\right) , \\ \Phi _{x y}&=\sigma _{x y}-\frac{1}{R e}\left( \frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\right) , \\ \Phi _{y y}&=\sigma _{y y}-\frac{2}{R e}\left( \frac{\partial v}{\partial y}\right) . \end{aligned} \end{aligned}$$

$\varvec{\sigma }$ is given by:

$$\begin{aligned} \begin{aligned} \sigma _{x x}&=\int _{-\infty }^{t} m(t-t')\\&\quad \frac{{a}}{{a}-3+{b} I_{1}+(1-{b}) I_{2}} {B_{t^{\prime }_{x x}}} d t^{\prime } \\ \sigma _{x y}&=\int _{-\infty }^{t} m(t-t')\\&\quad \frac{{a}}{{a}-3+{b} I_{1}+(1-{b}) I_{2}} {B_{t^{\prime }_{x y}}} d t^{\prime } \\ \sigma _{y y}&=\int _{-\infty }^{t} m(t-t')\\&\quad \frac{{a}}{{a}-3+{b} I_{1}+(1-{b}) I_{2}} {B_{t^{\prime }_{y y}}} d t^{\prime } \end{aligned} \end{aligned}$$

with

$$\begin{aligned} m(t-t')= & {} \frac{S_0\Gamma (\beta )}{W i\;\alpha } \left( E_{\alpha ,\alpha +\beta -1}\left( \frac{{t'}-t}{W i }\right) \right. \nonumber \\{} & {} \left. +(1-\beta ) E_{\alpha ,\alpha +\beta }\left( \frac{{t'}-t}{W i }\right) \right) . \end{aligned}$$

(A4)

The Finger tensor components satisfies (see 37):

$$\begin{aligned} \begin{aligned} \frac{\partial }{\partial t} {B_{t^{\prime }_{x x}}}&=-\frac{\partial }{\partial x}\left( u {B_{t^{\prime }_{x x}}}\right) -\frac{\partial }{\partial y}\left( v {B_{t^{\prime }_{x x}}}\right) \\&\quad +2\left[ \frac{\partial u}{\partial x} {B_{t^{\prime }_{x x}}}+\frac{\partial u}{\partial y} {B_{t^{\prime }_{x y}}}\right] \\ \frac{\partial }{\partial t} {B_{t^{\prime }_{x y}}}&=-\frac{\partial }{\partial x}\left( u {B_{t^{\prime }_{x y}}}\right) -\frac{\partial }{\partial y}\left( v {B_{t^{\prime }_{x y}}}\right) \\&\quad +\frac{\partial u}{\partial y} {B_{t^{\prime }_{y y}}}+\frac{\partial v}{\partial x} {B_{t^{\prime }_{x x}}} \\ \frac{\partial }{\partial t} {B_{t^{\prime }_{y y}}}&=-\frac{\partial }{\partial x}\left( u {B_{t^{\prime }_{y y}}}\right) -\frac{\partial }{\partial y}\left( v {B_{t^{\prime }_{y y}}}\right) \\&\quad +2\left[ \frac{\partial v}{\partial x} {B_{t^{\prime }_{x y}}}+\frac{\partial v}{\partial y} {B_{t^{\prime }_{y y}}}\right] \end{aligned} \end{aligned}$$

We assume the flow is fully developed and look for velocity, pressure and stress tensors of the form $u=u(y)$, $v(y)=0$, $p=p(x, y)$, $\sigma _{x y}=\sigma _{x y}(y)$, $\sigma _{x x}=\sigma _{x x}(y)$, $\sigma _{y y}=\sigma _{y y}(y)$, for $y \in [0,1]$.

The components of the Finger tensor are defined for simple shear flows by:

$$\begin{aligned} \textbf{B}_{t'}(\textbf{x},t)) = \begin{pmatrix} 1+ \dot{\gamma }^2(t-t')^2 &{} \dot{\gamma }(t-t') &{} 0\\ \dot{\gamma }(t-t') &{} 1 &{} 0\\ 0 &{} 0 &{} 1 \end{pmatrix} \end{aligned}$$

leading to the invariants $I_{1}$ and $I_{2}$, needed to compute $H\left( I_{1}, I_{2}\right)$ in the Papanastasiou equation.

$$\begin{aligned} I_{1}=I_{2}=3+\dot{\gamma }^{2}\left( t-t^{\prime }\right) ^{2} \end{aligned}$$

Here, $\dot{\gamma }=\displaystyle \frac{\textrm{d} u}{\,\textrm{d} y}$.

The momentum equation becomes,

$$\begin{aligned} -\frac{\partial p}{\partial x} + \frac{\partial \sigma _{x y}}{\partial y}= & {} 0 \end{aligned}$$

(A5)

$$\begin{aligned} -\frac{\partial p}{\partial y} + \frac{\partial \sigma _{y y}}{\partial y}= & {} 0 \end{aligned}$$

(A6)

and the stress tensor components are given by

$$\begin{aligned} \sigma _{x y}= & {} \int _{-\infty }^{t} m(t-t')\frac{\dot{\gamma }\left( t-t^{\prime }\right) }{{a}+\dot{\gamma }^{2}\left( t-t^{\prime }\right) ^{2}} d t^{\prime } \end{aligned}$$

(A7)

$$\begin{aligned} \sigma _{x x}= & {} \int _{-\infty }^{t} m(t-t')\frac{\left[ 1+\dot{\gamma }^{2}\left( t-t^{\prime }\right) ^{2}\right] }{{a}+\dot{\gamma }^{2}\left( t-t^{\prime }\right) ^{2}} d t^{\prime }, \end{aligned}$$

(A8)

$$\begin{aligned} \sigma _{y y}= & {} \int _{-\infty }^{t} m(t-t')\frac{1}{{a}+\dot{\gamma }^{2}\left( t-t^{\prime }\right) ^{2}} d t^{\prime } \end{aligned}$$

(A9)

with $m(t-t')$ defined by (A4).

From Eq. (A6) we have that

$$\begin{aligned} p(x, y)=\sigma _{y y}(y)+F(x), \end{aligned}$$

(A10)

therefore, Eq. (A5) can be rewritten as

$$\begin{aligned} \frac{\partial \sigma _{x y}}{\partial y}=F^{\prime }(x). \end{aligned}$$

(A11)

We can argue that the left side of Eq. (A11) is only a function of y, leading to $F^{\prime }$ being a constant,

$$\begin{aligned} C=\frac{\partial p}{\partial x}. \end{aligned}$$

Therefore,

$$\begin{aligned} \sigma _{x y}(y)=C y+h(x). \end{aligned}$$

Since $\sigma _{x y}(1 / 2)=0$, by symmetry, we have that $h(x)=-\frac{1}{2} C$, leading to

$$\begin{aligned} \sigma _{x y}(y)=(y-1 / 2) C. \end{aligned}$$

(A12)

Case $\alpha \ne 1$ or $\beta \ne 1$:

The memory function is given by (A4), then using (A12) Eq. (A7) can be written as,

$$\begin{aligned} \begin{aligned} (y-1 / 2) C&=\frac{S_0 \Gamma (\beta ){a}}{\alpha W i} \int _{-\infty }^{t} \left[ E_{\alpha ,\alpha +\beta -1}\left( \frac{{t'}-t}{W i }\right) \right. +\\&\quad \left. + (1-\beta ) E_{\alpha ,\alpha +\beta }\left( \frac{{t'}-t}{W i }\right) \right] \\&\quad \frac{\dot{\gamma }\left( t-t^{\prime }\right) }{{a}+\dot{\gamma }^{2}\left( t-t^{\prime }\right) ^{2}} d t^{\prime }. \end{aligned} \end{aligned}$$

The change of variable $s=t-t^{\prime }$ leads to

$$\begin{aligned} (y-1 / 2) C= & {} \frac{S_0 \Gamma (\beta ){a}}{\alpha W i} \int _{0}^{\infty } \left( E_{\alpha ,\alpha +\beta -1}\left( \frac{-s}{W i }\right) \right. \nonumber \\{} & {} \left. +(1-\beta ) E_{\alpha ,\alpha +\beta }\left( \frac{-s}{W i }\right) \right) \nonumber \\{} & {} \frac{\dot{\gamma } s}{ {a}+\dot{\gamma }^{2} s^{2}} d s \end{aligned}$$

(A13)

Let us define $s=\left( \varvec{a}^{1 / 2} / \dot{\gamma }\right) \mathcal {T}$. Then, we have that

$$\begin{aligned} (y-1 / 2) C= & {} K \int _{0}^{\infty } \left( E_{\alpha ,\alpha +\beta -1}\left( -\phi \mathcal {T}\right) \right. \nonumber \\{} & {} \left. +(1-\beta ) E_{\alpha ,\alpha +\beta }\left( -\phi \mathcal {T}\right) \right) \frac{\mathcal {T}}{1+\mathcal {T}^{2}} d \mathcal {T} \end{aligned}$$

(A14)

where

$$\begin{aligned} \phi =\frac{{a}^{1 / 2}}{\dot{\gamma } W i}, \quad K=\frac{S_0 \Gamma (\beta ){a}}{\dot{\gamma } \alpha W i}. \end{aligned}$$

(A15)

The integral in Eq. (A14) can be solved numerically by a higher order quadrature method.

Taking the boundary conditions at the channel entrance, where the velocity $u^{i n}(y)$ is defined by

$$\begin{aligned} u^{i n}(y)=1-4\left( y-\frac{1}{2}\right) ^{2}, \quad 0 \le y \le 1 \end{aligned}$$

we have

$$\begin{aligned} \int _{0}^{1} u^{i n}(y) d y=\frac{2}{3} \end{aligned}$$

and by mass conservation, the solution u(y) should verify the following equation:

$$\begin{aligned} \int _{0}^{1} u(y) d y=\frac{2}{3}. \end{aligned}$$

(A16)

Integrating Eq. (A16) by parts leads to

$$\begin{aligned} \int _{0}^{1} y \dot{\gamma }(y) d y+\frac{2}{3}=0. \end{aligned}$$

(A17)

To obtain the solutions $u(y),\; p(x, y),\; \sigma _{x y}(y),\; \sigma _{x x}(y)$ and $\sigma _{y y}(y)$, first, calculate $\dot{\gamma }$ and C by the following procedure:

• P1- Estimate a value for C; which can initially be taken to be the corresponding Newtonian value;

• P2- Discretize the interval $\left( \frac{1}{2}, 1\right]$, corresponding to half the channel, considering N points, $y_j=\displaystyle \frac{1}{2}+j \frac{1}{2 N},\, j=1,2, \cdots , N$.

• P3- Calculate the profile $\dot{\gamma }_{j}$, finding the zeros of the function (A13), $j=1,2, \cdots , N;$

  • − Note that we need to approximate the integral on the right side of (A13). For that we need to truncate the domain of integration and use an appropriate quadrature rule.

• P4- Check if $\dot{\gamma }_{j}$ satisfies the Eq. (A17);

• P5- If the Eq. (A17) is not satisfied, provide another value for C and redo steps P3 and P4.

After obtaining the values for C and $\dot{\gamma }\left( y_{j}\right)$ and using the steps P1-P5, the value of $\sigma _{x y}\left( y_{j}\right)$ is obtained from (A12). The values of $\sigma _{x x}\left( y_{j}\right)$ and $\sigma _{y y}\left( y_{j}\right)$ can be obtained from (A8) and (A9), respectively.

Assuming the change of variables $s=t-t^{\prime }$ and taking $s={a}^{1 / 2} \mathcal {T} / \dot{\gamma }\left( y_{j}\right)$, these integrals simplify as

$$\begin{aligned} \sigma _{x x}\left( y_{j}\right)= & {} K \int _{0}^{\infty } \left( E_{\alpha ,\alpha +\beta -1}\left( \frac{-s}{W i }\right) \right. \nonumber \\{} & {} \left. +(1-\beta ) E_{\alpha ,\alpha +\beta }\left( \frac{-s}{W i }\right) \right) \frac{1+{a} \mathcal {T}^{2}}{1+\mathcal {T}^{2}} d \mathcal {T} \end{aligned}$$

(A18)

$$\begin{aligned} \sigma _{y y}\left( y_{j}\right)= & {} K \int _{0}^{\infty } \left( E_{\alpha ,\alpha +\beta -1}\left( \frac{-s}{W i }\right) \right. \nonumber \\{} & {} \left. +(1-\beta ) E_{\alpha ,\alpha +\beta }\left( \frac{-s}{W i }\right) \right) \frac{1}{1+\mathcal {T}^{2}} d \mathcal {T} \end{aligned}$$

(A19)

where K and $\phi$ are given by (A15). The integrals in (A18) and (A19) can be solved numerically.