Solving Nonlinear Parabolic Equations by a Strongly Implicit Finite-Difference Scheme

TL;DR

This paper presents a stable, second-order accurate finite-difference scheme based on Crank--Nicolson for solving nonlinear parabolic equations modeling viscous gravity currents, with applications to thin-film flows and validation against analytical solutions.

Contribution

The authors develop and implement a strongly implicit, second-order finite-difference scheme for nonlinear parabolic PDEs, specifically applied to viscous gravity currents, demonstrating its stability and accuracy.

Findings

Scheme achieves second-order accuracy in space and time.

Numerical results confirm strong stability and mass conservation.

Applicable to various geometries and non-Newtonian fluids.

Abstract

We discuss the numerical solution of nonlinear parabolic partial differential equations, exhibiting finite speed of propagation, via a strongly implicit finite-difference scheme with formal truncation error $\mathcal{O}\left[(\Delta x)^2 + (\Delta t)^2 \right]$. Our application of interest is the spreading of viscous gravity currents in the study of which these type of differential equations arise. Viscous gravity currents are low Reynolds number (viscous forces dominate inertial forces) flow phenomena in which a dense, viscous fluid displaces a lighter (usually immiscible) fluid. The fluids may be confined by the sidewalls of a channel or propagate in an unconfined two-dimensional (or axisymmetric three-dimensional) geometry. Under the lubrication approximation, the mathematical description of the spreading of these fluids reduces to solving the so-called thin-film equation for the current's shape $h(x,t)$. To solve such nonlinear parabolic equations we propose a finite-difference scheme based on the Crank--Nicolson idea. We implement the scheme for problems involving a single spatial coordinate (i.e., two-dimensional, axisymmetric or spherically-symmetric three-dimensional currents) on an equispaced but staggered grid. We benchmark the scheme against analytical solutions and highlight its strong numerical stability by specifically considering the spreading of non-Newtonian power-law fluids in a variable-width confined channel-like geometry (a "Hele-Shaw cell") subject to a given mass conservation/balance constraint. We show that this constraint can be implemented by re-expressing it as nonlinear flux boundary conditions on the domain's endpoints. Then, we show numerically that the scheme achieves its full second-order accuracy in space and time. We also highlight through numerical simulations how the proposed scheme accurately respects the mass conservation/balance constraint.

Tables

Table 1. Table 1: Selected models of the propagation of viscous gravity currents herein simulated by a finite-difference scheme.

Case/Variable	$A$ [m${}_{}{}^{1-p+q}\cdot$s-1]	$p$ [–]	$q$ [–]	$\psi$ [m]
\svhline Newtonian fluid, fixed-width HS cell: $b\left(x\right)=b_1$. (see Huppert and Woods (1995))	$\frac{\Delta \rho g{b}_1^2}{12\mu }$	0	0	$h$
Newtonian fluid, variable-width HS cell: $b\left(x\right)=b_1{x}^n$. (see Zheng et al (2014))	$\frac{\Delta \rho g{b}_1^2}{12\mu }$	$n$	$3n$	$h$
Power-law fluid: $\mu ={\mu }_0{\dot{\gamma }}^{r-1}$, variable-width HS cell: $b\left(x\right)=b_1{x}^n$. (see Di Federico et al (2017); Longo (2017))	$\left(\frac{r}{2r+1}\right){\left(\frac{\Delta \rho g}{{\mu }_0}\right)}^{1/r}{\left(\frac{b_1}{2}\right)}^{\left(r+1\right)/r}$	$n$	$n\left(\frac{2r+1}{r}\right)$	$h{\left|\frac{\partial h}{\partial x}\right|}^{\left(1-r\right)/r}$
Newtonian fluid, 2D porous medium, variable porosity: $\varphi \left(x\right)={\varphi }_1{x}^m$, variable permeability: $k\left(x\right)=k_1{x}^n$. (see Zheng et al (2014))	$\frac{\Delta \rho g{k}_1}{\mu {\varphi }_1}$	$m$	$n$	$h$
Power-law fluid: $\mu ={\mu }_0{\dot{\gamma }}^{r-1}$, 2D porous medium, variable porosity: $\varphi \left(x\right)={\varphi }_1{x}^m$, variable permeability: $k\left(x\right)=k_1{x}^n$. (see Ciriello et al (2016))	$\begin{array}{ccc}{2}^{\left(3r+1\right)/2}{\left(\frac{r}{3r+1}\right)}^{1/r}& & \\ \times {\left(\frac{k_1}{{\varphi }_1}\right)}^{\left(r+1\right)/2r}{\left(\frac{\Delta \rho g}{{\mu }_0}\right)}^{1/r}& & \end{array}$	$m$	$\frac{m\left(r-1\right)+n\left(r+1\right)}{2r}$	$h{\left|\frac{\partial h}{\partial x}\right|}^{\left(1-r\right)/r}$

Equations

$$\frac{\partial h}{\partial t}=\frac{A}{x^p}\frac{\partial }{\partial x}\left(x^q\psi \frac{\partial h}{\partial x}\right).$$

$$\mathcal{V}(t)\equiv {\int }_{\ell }^{L}h(x,t)b_1\varphi (x)\mathrm{d}x={\mathcal{V}}_0+{\mathcal{V}}_{\mathrm{in}}{t}^{\alpha },$$

$$\mathcal{V}(t)\equiv {\int }_{\ell }^{L}h(x,t)b(x)\mathrm{d}x={\mathcal{V}}_0+{\mathcal{V}}_{\mathrm{in}}{t}^{\alpha }.$$

$$\begin{gathered} \frac{\partial }{\partial t}{\int }_{\ell }^{L}h(x,t)b_1{x}^n\mathrm{d}x={\int }_{\ell }^L\frac{\partial h}{\partial t}{b}_1{x}^n\mathrm{d}x={\int }_{\ell }^L{b}_1{x}^n\frac{A}{x^n}\frac{\partial }{\partial x}\left(x^q\psi \frac{\partial h}{\partial x}\right)\mathrm{d}x \\ =A{b}_1{\left(x^q\psi \frac{\partial h}{\partial x}\right)|}_{x=\ell }^{x=L}\stackrel{\text{by }\text{eqref}eq:mas{s}_{c}onservn}{=}\frac{\mathrm{d}({\mathcal{V}}_{\mathrm{in}}{t}^{\alpha })}{\mathrm{d}t}=\alpha {\mathcal{V}}_{\mathrm{in}}{t}^{\alpha -1}. \end{gathered}$$

$${\left(x^q\psi \frac{\partial h}{\partial x}\right)|}_{x=\ell }$$

$${\left(\psi \frac{\partial h}{\partial x}\right)|}_{x=L}$$

$${\left(\psi \frac{\partial h}{\partial x}\right)|}_{x=\ell }=-\frac{\alpha B}{A{\ell }^q}{t}^{\alpha -1},\alpha >0.$$

$${\left(x^q\psi \frac{\partial h}{\partial x}\right)|}_{x\to 0}$$

$${\left(\psi \frac{\partial h}{\partial x}\right)|}_{x=L}$$

$$h_0(x)=\{\begin{array}{ll}a\left({\mathfrak{X}}_0^c-x^c\right),& x\le {\mathfrak{X}}_0,\\ 0,& x>{\mathfrak{X}}_0,\end{array}$$

$$h_0(x)=\{\begin{array}{ll}a\left(-1+b{\mathrm{e}}^{cx}\right),& x\le {\mathfrak{X}}_0,\\ 0,& x>{\mathfrak{X}}_0.\end{array}$$

$$\mathcal{L}[h]\equiv \frac{A}{x^p}\frac{\partial }{\partial x}\left(x^q\psi \frac{\partial h}{\partial x}\right).$$

$$\frac{h^{n+1}_{i}-h^{n}_{i}}{\Delta t}=\frac{1}{2}\big{(}\mathcal{L}_{d}\left[h_{i}^{n+1}\right]+\mathcal{L}_{d}\left[h_{i}^{n}\right]\big{)},$$

$${\mathcal{L}}_d\left[h_i^n\right]$$

$${\mathcal{L}}_d\left[h_i^{n+1}\right]$$

$$\frac{h_i^{n,k+1}-h_i^n}{\Delta t}=$$

$$+$$

$$S_{i+1/2}\equiv {\left(\frac{\partial h}{\partial x}\right)}_{x=x_{i+1/2}}\approx \frac{h_{i+1}-h_i}{\Delta x}.$$

$${\psi }_{i+1/2}^{n+1/2,k}$$

$${\psi }_{i-1/2}^{n+1/2,k}$$

$${\psi }_i\equiv {\left[h{\left|\frac{\partial h}{\partial x}\right|}^{(1-r)/r}\right]}_{x=x_i}\approx h_i{\left|\frac{h_{i+1}-h_{i-1}}{2\Delta x}\right|}^{(1-r)/r}.$$

$${\psi }_0$$

$${\psi }_N$$

$$\begin{gathered} \underset{\text{matrix subdiagonal coefficient}}{\underbrace{\left[-\frac{A\Delta t}{2(\Delta x{)}^2}\frac{x_{i-1/2}^q}{x_i^p}{\psi }_{i-1/2}^{n+1/2,k}\right]}}{h}_{i-1}^{n,k+1} \\ +\underset{\text{matrix diagonal coefficient}}{\underbrace{\left[1+\frac{A\Delta t}{2(\Delta x{)}^2}\left(\frac{x_{i+1/2}^q}{x_i^p}{\psi }_{i+1/2}^{n+1/2,k}+\frac{x_{i-1/2}^q}{x_i^p}{\psi }_{i-1/2}^{n+1/2,k}\right)\right]}}{h}_i^{n,k+1} \\ +\underset{\text{matrix superdiagonal coefficient}}{\underbrace{\left[-\frac{A\Delta t}{2(\Delta x{)}^2}\frac{x_{i+1/2}^q}{x_i^p}{\psi }_{i+1/2}^{n+1/2,k}\right]}}{h}_{i+1}^{n,k+1} \\ =h_i^n+\frac{A\Delta t}{2(\Delta x{)}^2}\left[\frac{x_{i+1/2}^q}{x_i^p}{\psi }_{i+1/2}^{n+1/2,k}(h_{i+1}^n-h_i^n)-\frac{x_{i-1/2}^q}{x_i^p}{\psi }_{i-1/2}^{n+1/2,k}(h_i^n-h_{i-1}^n)\right] \end{gathered}$$

$$\underset{0\le i\le N}{\max }\left|h_i^{n,K}-h_i^{n,K-1}\right|<10^{-8}\underset{0\le i\le N}{\max }\left|h_i^{n,K-1}\right|.$$

$$\left[1+\frac{A\Delta t}{(\Delta x{)}^2}\right]h_i^{n+1}=\frac{A\Delta t}{(\Delta x{)}^2}\left(h_{i-1}^{n+1}+h_{i+1}^{n+1}+h_{i-1}^n+h_{i+1}^n\right)+\left[1+\frac{A\Delta t}{(\Delta x{)}^2}\right]h_i^n.$$

$${\psi }_{1/2}^{n+1/2,k}\frac{1}{\Delta x}\left(h_1^{n,k+1}-h_0^{n,k+1}\right)$$

$$\frac{1}{\Delta x}\left(h_N^{n,k+1}-h_{N-1}^{n,k+1}\right)$$

$$h_1^{n,k+1}-h_0^{n,k+1}$$

$$h_N^{n,k+1}-h_{N-1}^{n,k+1}$$

$$h_1^{n,k+1}-h_0^{n,k+1}$$

$$h_N^{n,k+1}-h_{N-1}^{n,k+1}$$

$$\left|{\int }_{\ell }^{L}h(x,t)b_1{x}^n\mathrm{d}x-({\mathcal{V}}_0+{\mathcal{V}}_{\mathrm{in}}{t}^{\alpha })\right|,$$

$$\frac{1}{2}\left[\mathcal{S}(x_i,t^{n+1},h_i^{n+1})+\mathcal{S}(x_i,t^n,h_i^n)\right]$$

$${\left(\frac{\partial \Psi }{\partial x}\right)}_{x=x_i}\approx \frac{1}{2}\left[\left(\frac{{\Psi }_{i+1}^{n+1}-{\Psi }_{i-1}^{n+1}}{2\Delta x}\right)+\left(\frac{{\Psi }_{i+1}^n-{\Psi }_{i-1}^n}{2\Delta x}\right)\right].$$

Full text

1 Introduction

In his lucid 2015 book Questions About Elastic Waves Engelbrecht (2015), Engelbrecht asks “What is a wave?” and answers “As surprising as it may sound, there is no simple answer to this question.” Indeed, the definition of ‘wave’ depends on the physical context at hand Christov (2014). Although most wave phenomena in classical continuum mechanics are described by hyperbolic (wave) equations, one of the surprises of 20th century research into nonlinear partial differential equations (PDEs) is that certain parabolic (diffusion) equations also yield structures with finite speed of propagation. Two examples are (i) a linear diffusion equation with a nonlinear reaction term Kolmogorov et al (1991); Fisher (1937), and (ii) a diffusion equation that is nonlinear due to a concentration-dependent diffusivity Barenblatt (1952).111More specifically, Barenblatt (1952) (see also (Ostriker et al, 1992, p. 13)) credits the observation of finite-speed of propagation in a nonlinear diffusion equation to a difficult-to-find 1950 paper by Zeldovich and Kompaneets. Indeed, it is known that certain aspects of wave phenomena can be reduced to a problem of solving a parabolic PDE, as gracefully illustrated by Engelbrecht (Engelbrecht, 1997, Ch. 6) through a series of selected case studies; further examples include, but are not limited to: electromagnetic waves propagating along the earth’s surface Vlasov and Talanov (1995), seismic waves Jerzak et al (2002), underwater acoustics Jensen et al (2011), and the classical theory of nerve pulses (Engelbrecht, 1997, §6.4.2), which nowadays has been updated by Engelbrecht et al (2018a) to a nonlinear hyperbolic (wave) model in the spirit of the Boussinesq paradigm Engelbrecht et al (2018b); Christov et al (2007).

Of special interest to the present discussion are physical problems that are modeled by nonlinear parabolic PDEs. These nonlinear problems lack general, all-encompassing solution methodologies. Instead, finding a solution often involves methods that are specific to the nature of the governing equation or the physical problem that it describes (Evans, 2010, Ch. 4) (see also the discussion in Christov (2014) in the context of heat conduction). The classical examples of nonlinear parabolic PDEs admitting traveling wave solutions come from heat conduction (Zel’dovich and Raizer, 2002, Ch. X) (see also Straughan (2011)) and thermoelasticity Straughan (2011); Berezovski and Ván (2017). The sense in which these nonlinear parabolic PDEs admit traveling-wave and ‘wavefront’ solutions now rests upon solid mathematical foundations Gilding and Kersner (2004); Vázquez (2007), including the case of gradient-dependent nonlinearity Tedeev and Vespri (2015) (e.g., the last case in table 2.2 to be discussed below).

A classical example of a nonlinear parabolic PDE governing the finite-speed wave-like motion of a substance arises in the study of an ideal gas spreading in a uniform porous medium Barenblatt (1952). A similar nonlinear parabolic equation was derived for the interface between a viscous fluid spreading horizontally underneath another fluid of lower density ($\Delta\rho>0$ between the fluids) Huppert (1982). The motion of the denser fluid is dictated by a balance of buoyancy and viscous forces at a low Reynolds number (viscous forces dominate inertial forces). Such viscous gravity current flows are characterized by ‘slender’ fluid profiles i.e., they have small aspect ratios ($h/L\ll 1$, where $h$ and $L$ are typical vertical and horizontal length scales, respectively). Therefore, these flows can be modeled by lubrication theory (Leal, 2007, Ch. 6). Generically, one obtains a nonlinear parabolic equation for the gravity current’s shape $h$ as a function of the flow-wise coordinate $x$ and time $t$. The case of the spreading of a fixed mass of Newtonian fluid was originally explored contemporaneously by Didden and Maxworthy (1982) and Huppert (1982).

Being governed by a parabolic (irreversible) equation, these currents ‘forget’ their initial conditions after some time has elapsed; this is Barenblatt’s concept of intermediate asymptotics Barenblatt and Zel’dovich (1972); Barenblatt (1996). Moreover, the PDE (1) can be reduced to an ordinary differential equation (ODE) through a self-similarity transformation. If the similarity variable can obtained by a scaling (dimensional) analysis, then the solution is termed a self-similar solution of the first kind (Barenblatt, 1996, Ch. 3). Specifically, the transformation is $h(x,t)=\mathfrak{C}\,t^{\beta}f(\zeta)$ ($h$ in units222Throughout the chapter, we use SI units for all dimensional quantities. of meters), where $\zeta=x/(\eta_{N}t^{\delta})$ is the similarity variable (dimensionless), $f(\zeta)$ is the self-similar profile to be determined by solving an ODE, and $\mathfrak{C}$ and $\eta_{N}$ are dimensional consistency constants. The exponents $\beta$ and $\delta$ are obtained through scaling (dimensional analysis) of the governing PDE. As a representative example, consider the one-dimensional (1D) spreading of a fixed mass of fluid having an arbitrary ‘wavy’ initial shape, as shown in figure 1. Suppose the fluid’s shape $h(x,t)$ is governed by the linear diffusion equation $\partial h/\partial t=A\partial^{2}h/\partial x^{2}$ (taking $A=1$ m2/s in this example without loss of generality) subject to $(\partial h/\partial x)|_{x=\pm L}=0$. The initial condition (IC) is quickly ‘forgotten,’ and the ultimate asymptotic state (here, flat) is achieved after passing through an intermediate asymptotic regime. It is straightforward to determine the self-similarity transformation: $\beta=-\delta=-1/2$ and $f(\zeta)=\mathrm{e}^{-\zeta^{2}}$ (see, e.g., Barenblatt (1996) and the Appendix). Here, $\mathfrak{C}$ depends on the initial condition and $\eta_{N}=2$. The convergence of the rescaled $h(x,t)$ profiles towards $f(\zeta)$ can be clearly observed in figure 1(b). The IC is forgotten, and the profile converges onto the Gaussian intermediate asymptotic shape Barenblatt (1996). The profile $f(\zeta)$ is termed ‘universal’ because it is independent of $h(x,0)$.

Having illustrated the notion of first-kind self-similarity as intermediate asymptotics, let us summarize its use in studying the gravitational spreading of Newtonian viscous fluids in a variety of physical scenarios. For example, gravity currents arise in geophysical applications associated with flows through porous rocks Woods (2015) such as in ground water extraction Bear (1988), during oil recovery Huppert (2000); Felisa et al (2018), and during CO2 sequestration Huppert and Neufeld (2014). In these examples, $h(x,t)$ represents an interface between two immiscible fluids in the limit of large Bond number (gravity dominates surface tension). There is now an extensive literature featuring a wealth of exact and approximate analytical self-similar solutions for gravity currents in porous media, e.g., Barenblatt (1952); Huppert and Woods (1995); Anderson et al (2003); Lyle et al (2005); Vella and Huppert (2006); Hesse et al (2007); Anderson et al (2010); De Loubens and Ramakrishnan (2011); Ciriello et al (2013); Huppert et al (2013); Zheng et al (2014) amongst many others.

In this chapter, we focus on the propagation of non-Newtonian gravity currents, specifically ones for which the denser fluid obeys a power-law rheology. This tractable model of non-Newtonian rheological response is also known as the Oswald–de Weale fluid Bird et al (1987). In unidirectional flow, the power-law model simply dictates that fluid’s viscosity depends upon a power of the velocity gradient. Di Federico et al (2006) generalized Huppert’s problem Huppert (1982) to power-law fluids, although Gratton et al Gratton et al (1999); Perazzo and Gratton (2003) had also considered some related problems. Even earlier, Kondic et al (1996, 1998) derived the governing equations for power-law fluids under confinement (i.e., in Hele-Shaw cells) using the lubrication approximation. These works have contributed to the use of a modified Darcy law to model the flow of non-Newtonian fluids in porous media using the analogy to flow in Hele-Shaw cells. Aronsson and Janfalk (1992) were perhaps the first to combine a Darcy law for a power-law fluid with the continuity equation to obtain a single PDE, of the kind studied herein, governing the gravity current’s shape. Recently, Lauriola et al (2018) highlighted the versatility of this approach by reviewing the existing literature and extending it to two-dimensional axisymmetric spreading in media with uniform porosity but variable permeability. All these flows are of interest because exact analytical self-similar solutions in closed form have been derived previously Gratton et al (1999); Perazzo and Gratton (2005); Di Federico et al (2012, 2017); Ciriello et al (2016). Specifically, the solution of Ciriello et al (2016) will be used in §4.1 below to verify the truncation error of the proposed numerical method.

For a self-similar solution to exist, both the governing PDE and its boundary conditions (BCs) must properly transform into an ODE in $\zeta$ with suitable BCs. A number of studies have specifically shown that the volume of fluid within the domain can be transient, varying as a power law in time, $\mathcal{V}(t)\propto t^{\alpha}$ ($\alpha\geq 0$), and a self-similar solution still exists (see, e.g., Barenblatt (1952); Lyle et al (2005); Hesse et al (2007); Di Federico et al (2012); Zheng et al (2014); Di Federico et al (2017) and the references therein). However, the nonlinear ODE in $\zeta$ often cannot be integrated exactly in terms of known function, except for $\alpha=0$. In §2.2 below, we discuss how a constraint of the form $\mathcal{V}(t)\propto t^{\alpha}$ can be implemented numerically through flux BCs at the computational domain’s ends.

With increasing complexity of the flow physics incorporated in the model, finding a self-similarity transformation may no longer be possible simply by scaling (dimensional) arguments. Gratton and Minotti (1990) classified a number of such situations, including the so-called ‘focusing’ flows involving fluid axisymmetrically flowing towards the origin on a flat planar surface. Further examples involving confined currents in channels with variable width, and/or in porous media whose permeability and porosity are functions of $x$, were proposed by Zheng et al (2014), as illustrated in figure 2. These gravity currents do enter a self-similar regime, even though a self-similar transformation cannot be obtained by scaling arguments alone. The exponents $\beta$ and $\delta$ in the transformation are unknown a priori, hence this situation represents a self-similarity of the second kind (Barenblatt, 1996, Ch. 4). The governing equation can be transformed to an ODE, following which a nonlinear eigenvalue problem must be solved for $\beta$ and $\delta$ through a phase-plane analysis Gratton and Minotti (1990); Gratton (1991). Alternatively, experiments or numerical simulations are necessary to determine $\beta$ and $\delta$. For example, early numerical simulations were performed to this end by Diez et al (1992). However, a ‘pre-wetting film’ ahead of the current’s sharp wavefront ($x=x_{f}(t)$ where h\big{(}x_{f}(t),t\big{)}=0) was required to avoid numerical instabilities. The scheme therein was also first-order accurate in time only. In this chapter, we propose a modern, high-order-accurate implicit numerical method for use in such problems.

Specifically, we develop and benchmark a strongly implicit and conservative numerical scheme for 1D nonlinear parabolic PDEs arising in the study of gravity currents. We show how the proposed scheme can be used to simulate (with high accuracy and at low computational expense) the spreading of 1D non-Newtonian viscous gravity currents in variable geometries (specifically, Hele-Shaw cells with widths varying as a power law in $x$). To this end, we build upon the work of Zheng et al (2014), which introduced this type of finite-difference scheme for simulating the spreading of a finite mass of Newtonian fluid in a variable-width Hele-Shaw cell. Owing to its accuracy and stability, this finite-difference scheme has been recently applied by Alhashim and Koch (2018) to study hydraulic fracturing of low-permeability rock.

This chapter is organized as follows. In §2, we briefly summarize existing models describing certain flows of viscous gravity currents. Then, we introduce a convenient general notation for such nonlinear parabolic PDEs. In §3.1, we introduce the 1D equispaced but staggered grid upon which the proposed finite-difference scheme is to be implemented. The derivation of the BCs for the PDE, from the mass conservation constraint, is discussed in §2.2. Then, we construct the nonlinear Crank–Nicolson scheme in §3.2 and discuss the discretized form of the nonlinear flux BCs in §3.4. Continuing, in §4.1, the scheme’s accuracy is justified by comparing the numerical solution provided by the finite-difference scheme (up to a specified physical time) against an analytical solution obtained through a self-similar transformation of the PDE. Specifically, this approach involves three validation cases: (i) a symmetric (about $x=0$) lump of fixed fluid mass spreading in two directions (convergence is independent of BCs), (ii) a fixed fluid mass spreading away from the origin ($x=0$) (requires only no flux BCs), and (iii) a variable fluid mass injected at the origin spreading away from it (requires careful implementation of the nonlinear BCs). In all three cases, the scheme is shown to be capable of accurately computing the evolution of gravity current’s shape. In §4.2, we analyze the scheme’s conservation properties by verifying numerically that it respects the mass constraint $\mathcal{V}(t)\propto t^{\alpha}$. We consider two validation cases: (i) release of a fixed fluid mass ($\alpha=0$), and (ii) fluid mass injection into the domain ($\alpha>0$). In both cases, we specifically focus on the challenging case of a non-Newtonian (power-law) displacing fluid in a variable-width channel. As a benchmark, we use previously derived first-kind self-similar solutions from the literature, which are discussed in the Appendix.

2 Preliminaries

In this section, we summarize the mathematical model for viscous gravity currents in a selected set of applications involving Newtonian and non-Newtonian fluids. We study their spreading in a fixed- or variable-width channel geometry (also known as a “Hele-Shaw cell”), as well as flows in heterogeneous porous media with independently variable permeability and porosity. Our goal is to highlight the fact that all these models can be concisely summarized by a single nonlinear parabolic PDE supplemented with a set of nonlinear Neumann (flux) BCs.

2.1 Fluid Domain and Flow Characteristics

The flow domain is assumed to be long and thin. For example, it can be a channel existing in the gap between two impermeable plates, i.e., a Hele-Shaw (HS) cell, which may or may not have variable transverse (to the flow) width as shown in figure 2(a); or, it can be slab of uniform-thickness heterogeneously porous material, as shown in figure 2(c). The viscous gravity current consists of one fluid displacing another immiscible fluid. Therefore, a sharp interface $y=h(x,t)$ separates the two fluids at all times. The present study considers the limit of negligible surface (interfacial) tension (compared to gravitational forces). The density difference $\Delta\rho$ between the two fluids is large compared to the density of the lighter fluid, and the denser fluid flows along the bottom of the cell, which is a horizontal impermeable surface. In doing so, the denser fluid displaces the lighter fluid out of its way. Here, the geometry is considered to be vertically unconfined so that the details of the flow of the upper (lighter) fluid can be neglected.

We are interested in the evolution of the interface $h(x,t)$ between the two fluids. Owing to the vertically unconfined, long and thin geometry of the flow passage, the denser fluid has a slender profile (small aspect ratio), and the fluid flow can be described by lubrication theory. The lubrication approximation also requires that viscous forces dominate inertial forces; this is the limit of small Reynolds number. In this regime of small Reynolds number but large Bond number, the flow is governed by a balance of viscous forces and gravity. Furthermore, the lubrication approximation allows for (at the leading order in the aspect ratio) the variation of quantities across the transverse direction, as well as the vertical velocities of the fluids to be neglected.

As shown in figure 2(a), for the flow in a HS cell, we allow the cell’s width to vary as a power-law of the streamwise coordinate $x$, i.e., $b(x)=b_{1}x^{n}$, where $n\geq 0$ is a dimensionless exponent, and $b_{1}>0$ is a dimensional consistency constant having units m1-n. Since the cell has a variable width, it originates from a cell ‘origin,’ which is always taken to be $x=0$ such that $b(0)=0$. As discussed in Zheng et al (2014), in such a flow geometry, the lubrication approximation may fail when $b(x)$ is an increasing function of $x$ i.e., $\mathrm{d}b/\mathrm{d}x=nb_{1}x^{n-1}>1$. In such quickly-widening cells, the transverse variations of properties become significant. We ensure the validity of the lubrication approximation, and models derived on the basis of it, by only considering $n<1$ such that $\mathrm{d}b/\mathrm{d}x$ remains a decreasing function of $x$.

The porosity can also be varied by filling the HS cell with beads of fixed diameter, as illustrated in figure 2(b). We also consider a gravity current spreading horizontally in a porous slab of constant transverse width ($b(x)=b_{1}=const.$) with heterogeneous porosity $\phi(x)=\phi_{1}x^{m}$ and permeability $k(x)=k_{1}x^{n}$, as shown in figure 2(c). Here, $m,n\geq 0$ are dimensionless exponents and $\phi_{1},k_{1}>0$ are dimensional constants needed for consistency with the definitions of porosity and permeability, respectively; specifically $\phi_{1}$ has units of units of m*-m*, and $k_{1}$ has units of m2-n. These variations are illustrated by the streamwise changes of bead radii in figure 2(c). Now, the point at which the porosity and permeability vanish is the origin of the cell. Another interesting case, that of a medium with vertically heterogeneous porosity, has been explored by Ciriello et al (2016). In this chapter, we limit our discussion to flow in a completely porous (i.e., unobstructed, $\phi=1$) HS cell of variable width as in figure 2(a). However, the numerical scheme developed herein can readily treat any of these cases, taking the appropriate parameter definitions from table 2.2 in §2.2.

We allow the denser fluid to be non-Newtonian. Specifically, it obeys the power-law rheology. In unidirectional flow, the one unique non-trivial shear stress component is given by $\tau=\mu(\dot{\gamma})\dot{\gamma}$, where the dynamic viscosity $\mu$ depends on the shear rate $\dot{\gamma}$ as $\mu(\dot{\gamma})=\mu_{0}\dot{\gamma}^{r-1}$. Here, $\mu_{0}$ is the flow consistency index (units of Pa$\cdot$sr), and $r$ ($>0$) is the fluid’s rheological index. Fluids having $r<1$ are termed shear-thinning (e.g., blood), and fluids with $r>1$ are termed shear-thickening (e.g., dense particulate suspensions). In the special case $r=1$, the power-law model reduces to the Newtonian fluid. As stated above, the flow of the displaced fluid is immaterial to the dynamics of the gravity current, as long as the viscosity and density contrasts are large. This condition is satisfied, e.g., by assuming (for the purposes of this chapter) the displaced fluid is air.

Finally, the volume of the fluid in the cell itself may be either fixed (constant mass) or vary with time (injection). Consistent with the literature, we consider the instantaneous volume of fluid in the cell to increase as a power law in $t$: $\mathcal{V}(t)=\mathcal{V}_{0}+\mathcal{V}_{\mathrm{in}}t^{\alpha}$, where $\mathcal{V}_{0}$ is the initial volume of fluid in the HS cell (measured in m3), $\alpha\geq 0$ is a dimensionless exponent, and $\mathcal{V}_{in}$ is an injection pseudo-rate (in units m3s*-α*), becoming precisely the injection rate for $\alpha=1$. Next, we discuss how this assumption leads to BCs for the physical problem and for the numerical scheme.

2.2 Governing Equation, Initial and Boundary Conditions

The propagation of a viscous gravity current is described by a diffusion equation for the interface $h(x,t)$, which is the shape of profile of the denser fluid. The models are derived either from porous medium flow under Darcy’s law and the Dupoit approximation (Bear, 1988, Ch. 8) or using lubrication theory with no-slip along the bottom of the cell and zero shear stress at the fluid–fluid interface (Leal, 2007, Ch. 6-C). The resulting velocity field is combined with a depth-averaged continuity equation to derive the nonlinear parabolic PDE for $h(x,t)$. We propose to summarize all gravity current propagation along horizontal surfaces through a single ‘thin-film’ Oron et al (1997) equation:

[TABLE]

According to Engelbrecht (2015, Ch. 5), eq. (1) can be classified as an ‘evolution equation.’ The term in the parentheses on the right-hand side of eq. (1), roughly, represents a fluid flux balanced by the change in height on the left-hand side. The multiplicative factor $A/x^{p}$ arises due to (i) geometric variations of the flow passage in the flow-wise direction, (ii) porosity variations in the flow-wise direction, or (iii) from the choice of coordinate system in the absence of (i) or (ii). Here, $A$ is dimensional constant depending on the flow geometry, the domain, and the fluid properties. Additionally, $p$ and $q$ are dimensionless exponents that depend on the flow geometry and fluid rheology. The quantity denoted by $\psi$ represents specifically the nonlinearity in these PDEs. Thus, it is necessarily a function of $h$, and possibly $\partial h/\partial x$ for a non-Newtonian fluid (as in the third and fifth rows of table 2.2).333Interestingly, an ‘$r$-Laplacian’ PDE, similar to eq. (1) for a power-law fluid in a HS cell (third row of table 2.2), arises during fluid–structure interaction between a power-law fluid and an enclosing slender elastic tube Boyko et al (2017). This PDE can also be tackled by the proposed finite-difference scheme.

As stated in §1, several versions of eq. (1) will be explored herein, incorporating geometric variations, porosity variations, non-Newtonian behavior. The pertinent physical scenarios that will be tackled herein (using the proposed numerical scheme) are presented in table 2.2, which lists expressions for $A$, $p$, $q$ and $\psi$. From a dimensional analysis of the PDE (1), it follows that the constant $A$ must have units of m${}^{1+p-q}\cdot$s*-1*, as long as the nonlinearity $\psi$ has units of length (as is the case for all the models summarized in Table 2.2). It is worth noting that in the case of 1D, linear diffusion ($p=q=0$ and $\psi=1$), $A$ becomes the ‘diffusivity’ in units of m2/s.

The PDE (1) is solved on the finite space-time interval $(x,t)\in(\ell,L)\times(t_{0},t_{f}]$. Here, $t_{0}$ and $t_{f}$ represent the initial and final times of the numerical simulation’s run, respectively. An initial condition (IC) $h_{0}(x)$ is specified at $t=t_{0}$, so that $h(x,t_{0})=h_{0}(x)$ is known. Meanwhile, $\ell$ is a small positive value (close to 0). Boundary conditions (BCs) are specified at $x=\ell$ and $x=L$. These involve some combination of $h$ and $\partial h/\partial x$. The reason for taking $x=\ell\neq 0$ becomes clear below.

Thus, let us now discuss such a suitable set of BCs. The BCs are based on the imposed mass conservation/growth constraint. Consider the case of a viscous gravity current in a porous slab with variable porosity $\phi(x)=\phi_{1}x^{m}$, and transverse width $b_{1}=const.$ Then, the conservation of mass constraint (see Zheng et al (2014)) takes the form

[TABLE]

where $\alpha\geq 0$. In the parallel case of a HS cell with variable width $b(x)=b_{1}x^{n}$ and porosity $\phi_{1}=const.$, which can either be set to unity or absorbed into $b(x)$ via $b_{1}$, the mass constraint becomes

[TABLE]

Taking a time derivative of eq. (3) and employing eq. (1), we obtain

[TABLE]

Here, $p=n$ in this case of interest, as described in table 2.2, and $Ab_{1}=const.$ Thus, we have obtained conditions relating $x^{q}\psi\partial h/\partial x$ at $x=\ell$ and $x=L$ to $\alpha\mathcal{V}_{\mathrm{in}}t^{\alpha-1}$. These conditions, if satisfied, automatically take into account the imposed volume constraint from eq. (3). The calculation starting with eq. (2) is omitted as it is identical, subject to proper choice of $p$.

For the case of propagation away from the cell’s origin (i.e., any injection of mass must occur near $x=0$, specifically at $x=\ell$), to satisfy eq. (4), we can require that

[TABLE]

where $B=\mathcal{V}_{\mathrm{in}}/b_{1}$. Recall, the case of $\alpha>0$ represents mass injection. Although eq. (3) and eqs. (5) are equivalent, the imposition of the nonlinear BC in eq. (5a) must be approached with care. It should be clear that to impose a flux near the origin (at $x=0$), we need \big{(}x^{q}\psi\partial h/\partial x\big{)}\big{|}_{x\to 0} to be finite. Then, $\psi\partial h/\partial x=\mathcal{O}(1/x^{q})$ as $x\to 0$. On the spatial domain $x\in(0,L)$, such an asymptotic behavior is possible for $p=q=0$. However, in a variable-width cell ($p,q\neq 0$), the local profile and slope as $x\to 0$ blow up if they are to satisfy $\psi\partial h/\partial x=\mathcal{O}(1/x^{q})$ as $x\to 0$. To avoid this uncomputable singularity issue, we defined the computational domain to be $x\in(\ell,L)$, where $\ell$ is ‘small’ but $>0$. The BC from eq. (5a) at $x=\ell$ can then be re-written as

[TABLE]

It may also be of interest to consider the case of a gravity current released a finite distance away from the origin and then spreading towards $x=0$. In this case, an additional length scale arises in the problem: the initial distance of the current’s edge from the origin, say $x_{f}(0)$. The existence of this extra length scale complicates the self-similarity analysis, leading to solutions of the second-kind (Barenblatt, 1996, Ch. 4), as discussed in §1. However, the numerical scheme can handle this case just as well; in fact, it requires no special consideration, unlike spreading away from the origin. Now, we may simply take $\ell=0$ and consider spreading on the domain $(0,L)$ subject to the following BCs:

[TABLE]

which together allow us to satisfy eq. (4) and, thus, eq. (3) for all $t\in(t_{0},t_{f}]$.

The most significant advantage of defining nonlinear flux BCs, such as those in eqs. (5) or (7), is that a nonlinear nonlocal (integral) constraint, such as that in eq. (2) or (3), no longer has to be applied onto the solution $h(x,t)$. Furthermore, if we start with compact initial conditions, i.e., there exists a nose location $x=x_{f}(t_{0})$ such that h\big{(}x_{f}(t_{0}),t_{0}\big{)}=0, then the finite-speed of propagation property of the nonlinear PDE (1) Gilding and Kersner (2004); Vázquez (2007) ensures that this nose $x_{f}(t)$ exists for all $t>t_{0}$ and h\big{(}x_{f}(t),t\big{)}=0 as well. The proposed fully-implicit scheme inherits this property of the PDE. Therefore, we can solve the PDE on the fixed domain $x\in(\ell,L)$, without any difficulty, instead of attempting to rescale to a moving domain on which $x_{f}(t)$ is one of the endpoints with $h=0$ as the BC applied there. The latter approach proposed by Bonnecaze et al (1993) (and used in more recent works Acton et al (2001) as well) leads to a number of additional variable-coefficient terms arising in the PDE (1), due to the non-Galilean transformation onto a shrinking/expanding domain. From a numerical methods point of view, having to discretize these additional terms is not generally desirable.

Having defined a suitable set of BCs, the last remaining piece of information required to close the statement of the mathematical problem at hand is the selection of pertinent initial conditions (ICs). For the case of the release of a finite fluid mass ($\alpha=0$), an arbitrary polynomial IC may be selected, as long as it has zero slope at the origin ($x=0$), leading to satisfaction of the no-flux boundary condition (5a). To this end, let the IC be given by

[TABLE]

where $\mathfrak{X}_{0}$ is a ‘release-gate’ location defining the initial position of the current’s nose, i.e., $\mathfrak{X}_{0}=x_{f}(t_{0})$ and h\big{(}x_{f}(t_{0}),t_{0}\big{)}\equiv h_{0}(\mathfrak{X}_{0})=0. The constant $c>1$ is an arbitrary dimensionless exponent. Finally, $a$ (units of m1-c) is set by normalizing $h_{0}(x)$ such that the initial volume of fluid corresponds to the selected initial fluid volume, $\mathcal{V}_{0}$, via eq. (3).

The case of the release of a finite mass of fluid is particularly forgiving in how we set the IC, and its slope at $x=0$. In fact, we could even take $c=1$ in eq. (8) and the scheme will provide an initial flux of fluid at $t=t_{0}^{+}$, with $(\partial h/\partial x)_{x=0}=0$ thereafter. On the other hand, the case of mass injection ($\alpha>0$) governed by the nonlinear BCs is not as forgiving. By virtue of the ‘point-source’ mass injection at $x=\ell$, the slope at the origin rises sharply from the moment of mass injection. This very sharp rise has a tendency to introduce unphysical oscillations in the current profile when starting from the IC in eq. (8). To avoid this, we must select a ‘better’ IC, which has a shape more similar to the actual solution’s singularity near $x=0$. Having tested a few different options, we found that an exponential function works well:

[TABLE]

Here, $b$ (dimensionless) and $c$ (units of m*-1*) are positive constants, $\mathfrak{X}_{0}=\frac{1}{c}\ln{\frac{1}{b}}$ ensures that the IC has no negative values and a sharp wavefront, and $a$ (units of m) is set by normalizing $h_{0}(x)$ to the selected intial volume $\mathcal{V}_{0}$ via eq. (3), as above.

Finally, it should be noted that the IC from eq. (8) is not used in the convergence studies for finite initial mass (§LABEL:sec:full_conv and §2.2). Rather, the IC is taken to be the exact self-similar solution of Ciriello et al (2016) for a power-law fluid in a uniform-width ($n=0$) HS cell (see also the Appendix). The reasoning behind this particular choice is further expounded upon in §4.1.