/**
# Free surface flow of a Bingham fluid with Navier Stokes solver
An example of 2D complex flow over a plate with a free surface is presented here.
The configuration is periodic what is injected to the left comes from the right.
On the bottom solid wall there is a no slip boundary condition, on the top fixed free surface a slip condition is imposed.
We use the centered Navier-Stokes solver with regularization for viscosity.
## Equations for a Bingham fluid
The Bingham fluid is a non Newtonian fluid.
The Bingham viscosity that we will put in the Navier-Stokes solver is such that if
$||\tau|| \le \tau_y$ then there is no motion $D=0$
if the stress is high enough
$||\tau|| > \tau_y$
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 06, which is the Frobenius norm.
then stress tensor is linked to the shear rate tensor by
$$\tau_{ij} = 2 \mu_0 D_{ij} + \tau_y \frac{D_{ij}}{||D||}$$
where $D_{ij}$ is the shear strain rate tensor
(tenseur de taux de dÃ©formation) $D_{ij}=(u_{i,j}+u_{j,i})/2$,
the components in 2D:
$D_{11}=\frac{\partial u}{\partial x}$,
$D_{12} =\frac{1}{2}( \frac{\partial u}{\partial y}+ \frac{\partial v}{\partial x})$,
$D_{21} =D_{12} =\frac{1}{2}( \frac{\partial u}{\partial y}+ \frac{\partial v}{\partial x})$,
$D_{22}=\frac{\partial v}{\partial y}$
in the Euclidian norm we have:
$$||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)
hence
$$D_2^2= D_{ij}D_{ij}= ( \frac{\partial u}{\partial x})^2 + (\frac{\partial v}{\partial y})^2 + \frac{1}{2}( \frac{\partial u}{\partial y}+ \frac{\partial v}{\partial x})^2$$
and we have obviously
$||D|| = D_2/\sqrt{2}$
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.
Factorising with $2 D_{ij}$ to obtain a equivalent viscosity
$$\tau_{ij} = 2( \mu_0 + \frac{\tau_y}{2 ||D|| } ) D_{ij}=2( \mu_0 + \frac{\tau_y}{\sqrt{2} D_2 } ) D_{ij} $$
as defined by Balmforth
$$\tau_{ij} = 2 \mu_{eq} D_{ij} $$
with
$$\mu_{eq}= \mu_0 + \frac{\tau_y}{\sqrt{2} D_2 }$$
## 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)$,
this gives (mind square root of 2):
$D_2=\sqrt{D_{ij}D_{ij}} = \frac{1}{\sqrt{2}}\frac{\partial u}{\partial y}$.
$$ \tau_{12} = 2 \mu D_{12} + 2 \tau_y \frac{D_{12}}{ \sqrt{2}D_2} =
\mu \frac{\partial u}{\partial y} + \tau_y $$
Equilibrium between pressure gradient and viscosity (writting $\tau$ for a shorthand of $\tau_{12}$)
$$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
$$ \tau = (-\frac{\partial p}{\partial x})(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 :
$Y 0 && du < 1e-5)
return 1; /* stop */
}
/**
## Implementation of the Bingham viscosity
*/
event bingham(i++) {
scalar muB[];
/** Compute
$D_{11}=\frac{\partial u}{\partial x}$,
$D_{12} =\frac{1}{2}( \frac{\partial u}{\partial y}+ \frac{\partial v}{\partial x})$,
$D_{21} =D_{12} =\frac{1}{2}( \frac{\partial u}{\partial y}+ \frac{\partial v}{\partial x})$,
$D_{22}=\frac{\partial v}{\partial y}$
And where the second invariant is $D_2=\sqrt{D_{ij}D_{ij}}$
hence
$$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
$$\mu_{eq}= \mu_0 + \frac{\tau_y}{\sqrt{2} D_2 }$$
the final one is the min of of $\mu_{eq}$ and a large $\mu_{max}$, then the fluid flows always, it is not a solid, but a very viscous fluid.
*/
foreach() {
double D2 = 0.;
foreach_dimension() {
double dxx = u.x[1,0] - u.x[-1,0];
double dxy = (u.x[0,1] - u.x[0,-1] + u.y[1,0] - u.y[-1,0])/2.;
D2 += sq(dxx) + sq(dxy);
}
if (D2 > 0.) {
D2 = sqrt(D2)/(2.*Delta);
D2c[] = D2;
muB[] = min(tau_y/(sqrt(2.)*D2) + mu_0 , mumax ); }
}
boundary ({muB});
foreach_face() {
muv.x[] = (muB[] + muB[-1,0])/2.;
}
boundary ((scalar *){muv});
}
/**
Save profiles
*/
event profiles (t += .1)
{
FILE * fp = fopen("xprof", "w");
scalar shear[];
foreach()
shear[] = (u.x[0,1] - u.x[0,-1])/(2.*Delta);
boundary ({shear});
for (double y = 0.; y < 1.0; y += 1./pow(2.,LEVEL))
fprintf (fp, "%g %g %g %g\n", y, interpolate (u.x, L0/2, y), interpolate (shear, L0/2, y),
interpolate (D2c, L0/2, y));
fclose (fp);
}
/**
We adapt according to the error on the velocity field.
*/
event adapt (i++) {
// adapt_wavelet ({u}, (double[]){3e-2,3e-2}, 7, 4);
}
/**
## Results and plots
To run the program
~~~bash
qcc -g -O3 -o bingham_simple bingham_simple.c -lm
./bingham_simple
lldb bingham_simple
~~~
Plot of the velocity and shear
~~~gnuplot Velocity and stress profiles for Bingham flow
set xlabel "y"
set ylabel "u, shear"
Y = 1-.25
U = Y*Y/2/1
p[:][:]'xprof' u 1:2 t'U computed' w lp ,''u 1:3 t'tau comp.',(x