How should boundary conditions be applied when using finite-volume method?

Question

Following from my previous question I am trying to apply boundary conditions to this non-uniform finite volume mesh,

I would like to apply a Robin type boundary condition to the l.h.s. of the domain ($x=x_L)$, such that,

$$ \sigma_L = \left( d u_x + a u \right) \bigg|_{x=x_L} $$

where $\sigma_L$ is the boundary value; $a, d$ are coefficients defined on the boundary, advection and diffusion respectively; $u_x = \frac{\partial u}{\partial x}$, is the derivative of $u$ evaluated at the boundary and $u$ is the variable for which we are solving.

Possible approaches

I can think of two ways to implement this boundary condition on the above finite volume mesh:

1. A ghost cell approach.

   Write $u_x$ as a finite difference including a ghost cell.$$ \sigma_L = d \frac{u_1 - u_0}{h_{-}} + a u(x_L)$$

   A. Then use linear interpolation with points $x_0$ and $x_1$ to find the intermediate value, $u(x_L)$.

   B. Alternatively find $u(x_L)$ by averaging over the cells, $u(x_L) = \frac{1}{2}(u_0 + u_1)$

   In either case, the dependence on ghost cell can be eliminated in the usual way (via substitution into the finite volume equation).

2. An extrapolation approach.

   Fit a linear (or quadratic) function to $u(x)$ by using the values at points $x_1, x_2$ ($x_3$). This will provide the value at $u(x_L)$. The linear (or quadratic) function can then be differentiated to find an expression for the value of the derivative, $u_x(x_L)$, at the boundary. This approach does not use a ghost cell.

Questions

• Which approach of the three, (1A, 1B or 2) is "standard" or you would recommend?
• Which approach introduces the smallest error or is the most stable?
• I think I can implement the ghost cell approach myself, however, how can the extrapolation approach be implemented, does this approach have a name?
• Are there any stability difference between fitting a linear function or a quadratic equation?

Specific equation

I wish to apply this boundary to the advection-diffusion equation (in conservation form) with non-linear source term,

$$ u_t = -au_x + du_{xx} + s(x,u,t) $$

Discretising this equation on the above mesh using the $\theta$-method gives,

$$ w_{j}^{n+1} - \theta r_a w_{j-1}^{n+1} - \theta r_b w_{j}^{n+1} - \theta r_c w_{j+1}^{n+1} = w_j^n + (1-\theta) r_a w_{j-1}^n + (1-\theta) r_b w_j^n + (1-\theta) r_c w_{j+1}^n + s(x_j,t_n) $$

However for the boundary point ($j=1$) I prefer to use a fully implicit scheme ($\theta=1$) to reduce the complexity,

$$ w_{1}^{n+1} - r_a w_{0}^{n+1} - r_b w_{1}^{n+1} - r_c w_{2}^{n+1} = w_1^n + s_1^n $$

Notice the ghost point $w_0^{n+1}$, this will be removed by applying the boundary condition.

The coefficients have the definitions,

$$ r_a = \frac{\Delta t}{h_j}\left( \frac{ah_j}{2h_{-}} + \frac{d}{h_{-}} \right) $$

$$ r_b = - \frac{\Delta t}{h_j}\left( \frac{a}{2}\left[ \frac{h_{j-1}}{h_{-}} - \frac{h_{j+1}}{h_{+}} \right] + d\left[-\frac{1}{h_{-}} - \frac{1}{h_{+}} \right]\right) $$

$$ r_c = \frac{\Delta t}{h_j}\left(- \frac{ah_j}{2h_{+}} + \frac{d}{h_{+}} \right) $$

All the "$h$" variables are defined as in the above diagram. Finally, $\Delta t$ which is the time step (N.B. this is a simplified case with constant $a$ and $d$ coefficients, in practice the "$r$" coefficients are slightly more complicated for this reason).

Answer 1

This is rather a general remark on FVM than an answer to the concrete questions. And the message is that there shouldn't be the need for such an adhoc discretization of the boundary conditions.

Unlike in FE- or FD-methods, where the starting point is a discrete ansatz for the solution, the FVM approach leaves the solution untouched (at first) but averages on a segmentation of the domain. The discretization of the solution comes into play only when the obtained system of balance equations is turned into an algebraic equation system by approximating the fluxes across the interfaces.

In this sense, in view of the boundary conditions, I advise to stick to the continuous form of the solution as long as possible and to introduce the discrete approximations only at the very end.

Say, the equation $$ u_t = -au_x + du_{xx} + s(x,u,t) $$ holds on the entire domain. Then it holds on the subdomain $[0,h_1)$, and an integration in space gives \begin{align} \int_0^{h_1}u_t \text{d}x &=& \int_0^{h_1} \partial_x(-au + du_{x})\text{d}x &+& \int_0^{h_1} s(x,u,t)\text{d}x \\ &=& (-au + du_{x})|_{x=h_1}-(-au + du_{x})|_{x=0}&+&\int_0^{h_1} s(x,u,t)\text{d}x, \end{align} which is the contribution of the first cell to the equation system. Note that, apart from taking only averages, there has been no discretization of $u$.

But now, to turn this into an algebraic equation, one typically assumes that on cell $C_i$ the function $u$ is constant in space, i.e. $u(t,x)|_{C_i} = u_i(t)$. Thus, having associated $u(x_i)\approx u_i$, one can express $u_x|_{h_i}$ at the cell boarders via the difference quotient in $u_{i}$ and $u_{i+1}$. To express $u$ at the cell boarders one can use interpolation (i.e. central differences or upwind schemes).

What to do at the boundary? In the example, it is all about approximating $(-au + du_{x})|_{x=0}$, no matter what has been done to $u$ so far.

• Given $u|_{x=0}=g_D$ one can introduce a ghost cell and the condition that an interpolant between $u_0$ and $u_1$ is equal to $g_D$ at the boarder.

• Given ${u_x|}_{x=0} = g_N$ one can introduce a ghost cell and the condition that an approximation to the derivative between $u_0$ and $u_1$ matches $g_N$ at the boarder

• If the flux itself is prescribed: $(-a u + du_x )|_{x=0}= g_R$, there is no need for a discretization.

However, I am not sure, what to do in the case that there are Robin type bc's that do not match the flux directly. This, will need some regularization because of the discontinuity of the advection and diffusion parameters.

---

===> Some personal thoughts on FVM <===

• FVM is not a scaled FDM, as examples of 1D Poisson's equations on a regular grid often suggest
• There shouldn't be a grid in FVM, there should be cells with interfaces and, if necessary, centers
• That's why I think that a stencil formulation of the discretization is not suitable
• Assembling of the equation system should be done according to the discretization approach, i.e. by iterating over the cells rather than defining an equation for every unknown. I mean to think of the $i$-th row of the coefficient matrix as the part of the problem posed on cell $\Omega_i$, rather than of the equation that is associated with $u_i$.
• This is particularly important for 2D or 3D problems but may also help to have a clear notation in 1D: Make a difference between the volume (in 1D: length) of the cell, here $h_i$, and the distance between the centers, maybe in 1D: $d_i:=d_{i,i+1}=|x_i-x_{i+1}|$.

---