/** # Planar Couette flow of Generalized Newtonian Fluid This code extends the method used in [/sandbox/M1EMN/Exemples/bingham_simple.c](/../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: $$\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: $$D_{11}=\frac{\partial u}{\partial x}$$ $$D_{12} =\frac{1}{2}\left(\frac{\partial u}{\partial y}+ \frac{\partial v}{\partial x}\right)$$ $$D_{21} =D_{12} =\frac{1}{2}\left( \frac{\partial u}{\partial y}+ \frac{\partial v}{\partial x}\right)$$ $$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) is given by: $$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 $$\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 $$\tau_{ij} = 2\left(\mu_0 \|D_{ij}\|^{n-1} + \frac{\tau_y}{2 \|D_{ij}\|}\right)D_{ij}$$ $$\tau_{ij} = 2 \mu_{eq}D_{ij}$$ $$\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. $$\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)$. $$\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}$) $$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 = \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 : $$\left\{\tau<\tau_y, \frac{\partial u}{\partial y} = 0,\:\&\:u=U\:\forall\:Y 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$$-\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$$D_{11} = \frac{\partial u}{\partial x}D_{12} = \frac{1}{2}\left( \frac{\partial u}{\partial y}+ \frac{\partial v}{\partial x}\right)D_{21} = \frac{1}{2}\left( \frac{\partial u}{\partial y}+ \frac{\partial v}{\partial x}\right)D_{22} = \frac{\partial v}{\partial y}$$The second invariant is D_2=\sqrt{D_{ij}D_{ij}} (this is the Frobenius norm)$$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\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.$$ \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 ~~~bash #!/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](https://www.dropbox.com/sh/8at1yk9vigovdg7/AAAr-Td7p106Kt_3cIK4mg_ia?dl=0) # Bibliography * [Same example in Basilisk using the calculation of D2 at cell centers](http://basilisk.fr/sandbox/M1EMN/Exemples/bingham_simple.c) and its application to [1D Collapse](http://basilisk.fr/sandbox/M1EMN/Exemples/bingham_collapse_noSV.c) * [Related example in Gerris](http://gerris.dalembert.upmc.fr/gerris/tests/tests/couette.html) * [Related example with augmented Lagrangian](http://basilisk.fr/sandbox/popinet/poiseuille-periodic.c) * K. F. Liu and C. C. Mei Liu, K.F. and Mei, C.C., 1990. Approximate equations for the slow spreading of a thin sheet of Bingham plastic fluid. Physics of Fluids A: Fluid Dynamics, 2(1), pp.30-36.; [doi: 10.1063/1.857821](https://doi.org/10.1063/1.857821) * The Theoretical Formulations: Bird, R.B., 1987. Armstrong and RC Hassager, O.,“Dynamics of Polymeric Liquids”. v.1. * Balmforth, N.J., Craster, R.V., Rust, A.C. and Sassi, R., 2006. Viscoplastic flow over an inclined surface. Journal of Non-Newtonian Fluid Mechanics, 139(1-2), pp.103-127. [doi: 10.1016/j.jnnfm.2006.07.010](https://doi.org/10.1016/j.jnnfm.2006.07.010) */