# Planar Couette flow of Generalized Newtonian Fluid

This code extends the method used in /sandbox/M1EMN/Exemples/bingham_simple.c and generalizes it for any Power Law fluid (using regularization method). Another difference between the two is that this code calculates the second invariant of deformation tensor at the face-centers of the cells instead of the cell centers.

## Mathematical Formulations

Unlike the Newtonian fluids, non-Newtonian fluids do not have a linear stress-strain rate relationship. One way to represent the relationship is using the Generalized Newtonian fluid method: \displaystyle \tau = \tau_y + 2\mu_0D_{ij}^n The fluid is such that if \|\tau\| \le \tau_y then there is no motion D_{ij}=0\:\forall\:(i,j) if the stress is high enough \|\tau\| > \tau_y then there is motion

Note: that \|\tau\| is the modulus defined as the Euclidian norm \sqrt{\frac{1}{2}{\tau_{ij} \tau_{ij}}}. It is not \sqrt{\tau_{11}^2 + \tau_{12}^2} as in Balmorth et al. (2006), which is the Frobenius norm.

D_{ij} is the shear strain rate tensor (or the deformation tensor)

D_{ij}=(u_{i,j}+u_{j,i})/2: the components in 2D: \displaystyle D_{11}=\frac{\partial u}{\partial x} \displaystyle D_{12} =\frac{1}{2}\left(\frac{\partial u}{\partial y}+ \frac{\partial v}{\partial x}\right) \displaystyle D_{21} =D_{12} =\frac{1}{2}\left( \frac{\partial u}{\partial y}+ \frac{\partial v}{\partial x}\right) \displaystyle D_{22}=\frac{\partial v}{\partial y}

In the Euclidian norm we have: \displaystyle \|D\|=\sqrt{\frac{D_{ij}D_{ij}}{2}} The second invariant defined by D_2=\sqrt{D_{ij}D_{ij}} (this is the Frobenius norm) is given by: \displaystyle D_2^2= D_{ij}D_{ij}= \left( \frac{\partial u}{\partial x}\right)^2 + \left(\frac{\partial v}{\partial y}\right)^2 + \frac{1}{2}\left( \frac{\partial u}{\partial y}+ \frac{\partial v}{\partial x}\right)^2 and we have obviously \|D_{ij}\| = D_2/\sqrt{2}

## Numerical regularization

\displaystyle \tau_{ij} = \tau_y\left(\frac{D_{ij}}{\|D_{ij}\|}\right) + 2\mu_0\|D_{ij}\|^{n-1}D_{ij}^n Factorising with 2D_{ij} to obtain a equivalent viscosity \displaystyle \tau_{ij} = 2\left(\mu_0 \|D_{ij}\|^{n-1} + \frac{\tau_y}{2 \|D_{ij}\|}\right)D_{ij} \displaystyle \tau_{ij} = 2 \mu_{eq}D_{ij} \displaystyle \mu_{eq} = \mu_0\|D_{ij}\|^{n-1} + \frac{\tau_y}{2\|D_{ij}\|} \mu is the min of \mu_{eq} and a large \mu_{max} so that the viscosity does not blow up. \displaystyle \mu = \text{min}\left(\mu_{eq}, \mu_{max}\right) Note: We present here the formulation in Balmforth, he uses \dot{\gamma} which is by his definition \sqrt{\frac{1}{2}\dot{\gamma_{ij}}\dot{\gamma_{ij}}} and as \dot{\gamma_{ij}}=2 D_{ij} then \dot{\gamma} is \sqrt{2}D_2, that is why we have a \sqrt{2} in the equations.

## Exact solution in the proposed case

We look at an unidirectional flow, a pure shear flow u(y), v=0, so D_{11}=D_{22}=0 and D_{12}=D_{21}=(1/2)(\partial u/\partial y). \displaystyle \tau_{12} = 2\mu D_{12}^n + \tau_y = 2^{1-n}\mu\left(\frac{\partial u}{\partial y}\right)^n + \tau_y

Equilibrium between pressure gradient and viscosity (writting \tau for a shorthand of \tau_{12}) \displaystyle 0=-\frac{\partial p}{\partial x} + \frac{\partial \tau}{\partial y} as there is no stress at the free surface y=h, the stress is \displaystyle \tau = \left(-\frac{\partial p}{\partial x}\right)(h-y) the stress \tau increases from the free surface, as long as \tau<\tau_y, we are under the threshold, so shear is zero: \frac{\partial u}{\partial y} =0, hence velocity is constant, say it is U. Let us define Y=h-\tau_y/(-\frac{\partial p}{\partial x}), where \tau=\tau_y.

So : \displaystyle \left\{\tau<\tau_y, \frac{\partial u}{\partial y} = 0,\:\&\:u=U\:\forall\:Y<y<h\right\} Then going down: 0<y<Y we have \tau = 2^{1-n}\mu\left(\frac{\partial u}{\partial y}\right)^n + \tau_y. This gives: \displaystyle \tau_y + 2^{1-n}\mu\left(\frac{\partial u}{\partial y}\right)^n = \left(-\frac{\partial p}{\partial x}\right)(Y-y) After some straight forward manipulations: \displaystyle \left(\frac{\partial u}{\partial y}\right) = \left(\frac{1}{2^{1-n}\mu}\right)^{1/n}\left(-\frac{\partial p}{\partial x}\right)^{1/n}(Y-y)^{1/n} = A_{p\mu}(Y-y)^{1/n} and this allows to solve for the velocity profile \displaystyle u = \frac{nA_{p\mu}}{n+1}\left(Y^{\frac{n+1}{n}} - \left(Y-y\right)^{\frac{n+1}{n}}\right) which is indeed zero in y=0, and for y=Y, we have the plug flow u=U of value: \displaystyle U= \frac{n}{n+1}A_{p\mu}Y^{\frac{n+1}{n}} \displaystyle A_{p\mu} = \left(\frac{1}{2^{1-n}\mu}\right)^{1/n}\left(-\frac{\partial p}{\partial x}\right)^{1/n} For the present case -\frac{\partial p}{\partial x} = 1, \mu = 1, h = 1, which gives: A_{p\mu} = \frac{1}{2^{(1-n)/n}}, Y = 1 - \tau_y

# Code

#include "navier-stokes/centered.h"
char filename[80];
double tauy,mu_0,mumax;
double n;
int imax = 1e4;
#define dtmax (1e-3)
int main(int a, char const *arguments[])
{
sprintf (filename, "%s", arguments[1]);

init_grid (1<<6);
L0 = 1.0;
origin (0.0, 0.0);
DT = dtmax;
stokes = true;
TOLERANCE = 1e-5;

Values of yeild stress, viscosity, and coefficient.
Newtonian: \mu_0 = 1.0; \tau_y = 0. and n = 1
Power law \mu_0 = 1.0; \tau_y = 0. and n = 0.5
Herschel-Bulkley \mu_0 = 1.0; \tau_y = 0.25 and n = 0.5
Bingham \mu_0 = 1.0; \tau_y = 0.25 and n = 1

  mu_0 = 1.0;
tauy= 0.0;
n = 1.0;
if (a >= 3){
mu_0 = atof(arguments[2]);
}
if (a >= 4){
tauy = atof(arguments[3]);
}
if (a >= 5){
n = atof(arguments[4]);
}

the regularisation value of viscosity

  mumax=1000;

Right - left boundaries are periodic

  periodic (right);

slip at the top

  u.t[top] = neumann(0);
u.n[top] = neumann(0);
uf.n[top] = neumann(0);

no slip at the bottom

  u.n[bottom] = dirichlet(0);
uf.n[bottom] = dirichlet(0);
u.t[bottom] = dirichlet(0);

presure conditions are neumann 0.0

  p[top] = neumann(0);
pf[top] = neumann(0);
p[bottom] = neumann(0);
pf[bottom] = neumann(0);

run();
}

// un is used to search for a stationary solution
scalar un[];
// muv will be used as the face vector for viscosity
face vector muv[];

## Initialization event

event init (t = 0) {
// preparing viscosity to be used as Non-Newtonian fluid
mu = muv;

presure gradient mdpdx \displaystyle -\frac{dp}{dx} = 1

  const face vector mdpdx[] = {1.0,0.0};
a = mdpdx;

Initialy at rest

  foreach() {
u.x[] = 0;
u.y[] = 0;
}
foreach(){
un[] = u.x[];
}
dump (file = "start");
}

We look for a stationary solution.

event logfile (i += 500; i <= imax) {
double du = change (u.x, un);
fprintf(ferr, "i = %d: err = %g\n", i, du);
if (i > 0 && du < 1e-6){
dump (file = filename);
return 1; /* stop */
}
if (i==imax){
dump (file = filename);
}
}

event properties(i++) { // Overloading the properties event

## Implementation of generalized Newtonian viscosity

\displaystyle D_{11} = \frac{\partial u}{\partial x} \displaystyle D_{12} = \frac{1}{2}\left( \frac{\partial u}{\partial y}+ \frac{\partial v}{\partial x}\right) \displaystyle D_{21} = \frac{1}{2}\left( \frac{\partial u}{\partial y}+ \frac{\partial v}{\partial x}\right) \displaystyle D_{22} = \frac{\partial v}{\partial y} The second invariant is D_2=\sqrt{D_{ij}D_{ij}} (this is the Frobenius norm) \displaystyle D_2^2= D_{ij}D_{ij}= D_{11}D_{11} + D_{12}D_{21} + D_{21}D_{12} + D_{22}D_{22} the equivalent viscosity is \displaystyle \mu_{eq}= \mu_0\left(\frac{D_2}{\sqrt{2}}\right)^{N-1} + \frac{\tau_y}{\sqrt{2} D_2 } Note: \|D\| = D_2/\sqrt{2}

Finally, mu is the min of of \mu_{eq} and a large \mu_{max}.

The fluid flows always, it is not a solid, but a very viscous fluid. \displaystyle \mu = \text{min}\left(\mu_{eq}, \mu_{max}\right)

  double muTemp = mu_0;
foreach_face() {
double D11 = (u.x[] - u.x[-1,0]);
double D22 = ((u.y[0,1]-u.y[0,-1])+(u.y[-1,1]-u.y[-1,-1]))/4.0;
double D12 = 0.5*(((u.x[0,1]-u.x[0,-1])+(u.x[-1,1]-u.x[-1,-1]))/4.0 + (u.y[] - u.y[-1,0]));
double D2 = sqrt(sq(D11)+sq(D22)+2.0*sq(D12))/(Delta);
if (D2 > 0.0) {
double temp = tauy/(sqrt(2.0)*D2) + mu_0*exp((n-1.0)*log(D2/sqrt(2.0)));
muTemp = min(temp, mumax);
} else {
if (tauy > 0.0 || n < 1.0){
muTemp = mumax;
} else {
muTemp = (n == 1.0 ? mu_0 : 0.0);
}
}
muv.x[] = fm.x[]*(muTemp);
}
boundary ((scalar *){muv});
}

## Running the code

Use the following run.sh script

#!/bin/bash
qcc -O2 -Wall Couette_NonNewtonian.c -o Couette_NonNewtonian -lm
./Couette_NonNewtonian lastNewt 1.0 0.0 1.0
./Couette_NonNewtonian lastShThn 1.0 0.0 0.5
./Couette_NonNewtonian lastHB 1.0 0.25 0.5
./Couette_NonNewtonian lastBing 1.0 0.25 1.0

# Output and Results

The post-processing codes and simulation data are available at: PostProcess