Viscoplastic channel flow

The mosolov example solves the problem this library was written for: steady, isothermal, laminar flow of a Bingham (viscoplastic) fluid down a rectangular channel. It is the numerical scheme of

J. C. Robertson and R. C. Kerr, Isothermal dynamics of channeled viscoplastic lava flows and new methods for estimating lava rheology, J. Geophys. Res. 117, B01202 (2012), doi:10.1029/2011JB008550.

The physical problem

A Bingham fluid behaves as a rigid solid until the shear stress exceeds a yield strength $\tau_y$, then flows with plastic viscosity $\mu$. In a channel inclined so that gravity drives the flow along $x$, the only non-zero velocity component is $u(y, z)$ over the cross-section. Non-dimensionalising lengths on the flow depth $H$ and velocity on $gH^2/\nu$, the momentum balance becomes

\[ \nabla^2 u - B\,\nabla\!\cdot\!\left(\frac{\nabla u}{|\nabla u|}\right) = -1 \quad\text{where the fluid yields}, \qquad \nabla u = 0 \quad\text{in the plugs}, \]

with two dimensionless groups:

Symbol	Meaning	Code
$B = \dfrac{\tau_y}{\rho g H}$	Bingham number -- dimensionless yield strength	`MosolovSettings::binghamNumber`
$b = W/H$	channel aspect ratio -- channel width over depth	`MosolovSettings::channelAspectRatio`

channelAspectRatio is metadata, not geometry: the multigrid solve runs on a square half-channel with multigridSettings.aspectRatio = 1 (square cells, see discretisation). The value of $b$ feeds critical_bingham, the output file-name stem (A{b}B{B}) and the netCDF aspect_ratio attribute.

The flow fills the channel with a central plug (moving at the maximum velocity) and corner plugs welded to the walls, separated by a yielded layer. As $B$ rises the plugs grow; at a critical Bingham number $B^\*(b)$ they meet and the flow stops. For a rectangular channel (Appendix A of the paper, eq. 19):

\[ B^\*(b) = \frac{2 + b - \sqrt{4 + b^2 + (2\pi - 4)\,b}}{4 - \pi}, \]

which is exactly critical_bingham(aspect) in mosolov.hpp. For $b = 2$ this gives $B^\* = 2/(2 + \sqrt{\pi}) \approx 0.530$.

The boundary conditions (eq. 6) are no-slip on the wall ($u = 0$ at $y = b/2$) and floor ($u = 0$ at $z = 1$), and no-shear on the symmetry plane ($\partial u/\partial y = 0$ at $y = 0$) and the free surface ($\partial u/\partial z = 0$ at $z = 0$). These map onto the example boundary pattern: zero-Dirichlet right/bottom, zero-Neumann left/top.

The augmented-Lagrangian scheme

The $|\nabla u|$ in the denominator is singular at the yield surface, so the problem is not solved directly. Following Glowinski & Le Tallec (1989), the strain-rate vector $\dot\gamma$ is introduced as an independent unknown, constrained to equal $\nabla u$ by a Lagrange multiplier $\lambda$, and an Uzawa iteration decouples the nonlinearity (paper eqs. 8--12):

Initial guess. Solve the Poisson problem $(1 + \alpha)\,\nabla^2 u_0 = -1$; set $\lambda_1 = 0$.
For $n \ge 1$, with $s_n = \alpha\,\nabla u_{n-1} + \lambda_n$:

$$ \dot\gamma_n = \begin{cases} \left(1 - \dfrac{B}{|s_n|}\right)\dfrac{s_n}{\alpha}, & |s_n| \ge B, \[2ex] 0, & \text{otherwise.} \end{cases} $$

Velocity update. Solve the linear problem

$$ (1 + \alpha)\,\nabla^2 u_n = \nabla!\cdot!\bigl(\alpha\,\dot\gamma_n - \lambda_n\bigr) - 1. $$

Multiplier update. $\lambda_{n+1} = \lambda_n + \alpha\,(\nabla u_n - \dot\gamma_n)$.
Convergence. Stop when $\lVert \nabla u_n - \dot\gamma_n \rVert / \lVert \nabla u_n \rVert < \delta$.

$\alpha$ is the augmenting parameter (MosolovSettings::augmentingParameter): too small and the iteration crawls, too large and the plug regions go unstable.

Mosolov::solve is this loop. Step 1 and step 3 are each a single LinearMultigrid::multigrid full-multigrid solve of the 5-point Poisson problem -- which is why the whole library is built around a fast linear solver. Iteration counts run from 1 (Newtonian, $B = 0$) to ${\sim}1000$ near $B^\*$, with a median around 100.

What the output contains

Mosolov::write records the fields and the parameters that define the regime:

netCDF entry	Quantity
variable `velocity`	$u(y, z)$ on the finest grid
variable `strain_rate`	$\lvert\nabla u\rvert$ -- zero in the plugs, so it maps the yield surfaces
variable `log_residual`	$\log_{10}$ of the PDE residual
attribute `bingham_number`	$B$
attribute `aspect_ratio`	$b$
attribute `total_flux`	$q$, the down-channel volumetric flux per unit width (`FDArray::calculate_flux`)

Validation

The paper gives closed-form limits that the code is checked against:

Newtonian ($B = 0$). The problem reduces to Poisson channel flow with the Fourier-series solution in eq. 15. tests/test_mosolov.cpp checks that the $B = 0$ Mosolov solve matches a direct LinearMultigrid Poisson solve.
Flux vs. Bingham number. $q(B)$ must decrease monotonically to zero as $B \to B^\*(b)$; tests/test_mosolov.cpp checks the monotone decrease and that $q$ has collapsed well below its Newtonian value by $0.9\,B^\*$.
Wide-channel limit. As $b \to \infty$ the solution approaches the 2-D film flow of eq. 13, $u(z) = \tfrac12(1-B)^2$ for $0 \le z < 1-B$.

Symbol	Meaning	Code
\(B = \dfrac{\tau_y}{\rho g H}\)	Bingham number -- dimensionless yield strength	`MosolovSettings::binghamNumber`
\(b = W/H\)	channel aspect ratio -- channel width over depth	`MosolovSettings::channelAspectRatio`

netCDF entry	Quantity
variable `velocity`	\(u(y, z)\) on the finest grid
variable `strain_rate`	\(\lvert\nabla u\rvert\) -- zero in the plugs, so it maps the yield surfaces
variable `log_residual`	\(\log_{10}\) of the PDE residual
attribute `bingham_number`	\(B\)
attribute `aspect_ratio`	\(b\)
attribute `total_flux`	\(q\), the down-channel volumetric flux per unit width (`FDArray::calculate_flux`)