/** # 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