One finite difference scheme

Question

There is PDE:

$$\frac{\partial u(r,\varphi,\psi,t)}{\partial t}=\operatorname{div}A(r,\varphi,\psi)\nabla u +f(r,\varphi,\psi,t) $$

We solve numerically IBVP for the ball $B_{1}(0)\subset \mathbb{R}^3$ using the numerical scheme of this PDE. Numerical scheme is the following:

$$\frac{u^{n+1}_{i,j,k}-u^{n}_{i,j,k}}{\Delta t}=\frac{\oint_{\partial V_{i,j,k}}(A_{i,j,k}\nabla u^{n}_{i,j,k},\overline{n})ds}{V_{i,j,k}}+\overline{f}_{i,j,k}^{n}, $$ where $V_{i,j,k}$ is a cell of grid with center $(i,j,k)$ and vertexes coordinates $(i\pm 1/2, j\pm 1/2, k\pm 1/2)$; $\overline{f}_{i,j,k}^{n}$ is a mean value of $f$ for $V_{i,j,k}$. This scheme is equivalent to

$$u^{n+1}_{i,j,k}=u^{n}_{i,j,k}+\frac{\Delta t}{V_{i,j,k}}\oint_{\partial V_{i,j,k}}(A_{i,j,k}\nabla u^{n}_{i,j,k},\overline{n})ds+\frac{\Delta t}{V_{i,j,k}}\int\limits_{V_{i,j,k}}f(r,\varphi,\psi,t)dr d\varphi d\psi $$

Questions:

1) is this scheme conditionally stable?

2) To integrate over $\partial V_{i,j,k} $ we must know values of $u_{i\pm 1/2,j,k},u_{i,j\pm 1/2,k},u_{i,j,k\pm 1/2}$ on every time step. What are the methods to find this values on every time step?

3) How is this scheme called and where can I find more information about it?

Answer 1

The scheme you are showing is a fairly standard finite volume formulation coupled with explicit euler time integration.

Finite Volume Derivation

(Note: I noticed after a bit of writing that I had started to ramble but didn't want to scrap what I had. If you don't want to see the equation in the derived skip to the 2nd bold header)

I will step through the discretization of the unsteady heat equation with source (note the similarity to your equation):

$\frac{\partial u(t,\vec{r})}{\partial t} = \nabla \cdot (k(\vec{r}) \nabla u(t,\vec{r})) + S(\vec{r})$

First you discretize in space. Your unknowns are defined as the volumetric average of each cell (cells indexed with $i$). Integrating over each cell you get:

$\iiint_{V_i}\frac{\partial u_i(t,\vec{r})}{\partial t}dV = \iiint_{V_i} \nabla \cdot (k(\vec{r}) \nabla u_i(t,\vec{r})) + S(\vec{r}) dV$

Break the r.h.s. in two:

$\iiint_{V_i}\frac{\partial u_i(t,\vec{r})}{\partial t}dV = \iiint_{V_i} \nabla \cdot (k(\vec{r}) \nabla u_i(t,\vec{r}))dV + \iiint_{V_i}S(\vec{r}) dV$

The unsteady term, in the simplest form, can be estimated by taking the volume of the cell times the time derivative at the center point ($\vec{r}_i$). Then discretize using Euler:

${V_i}\frac{u_i^{n+1} - u_i^n}{\Delta t} = \iiint_{V_i} \nabla \cdot (k(\vec{r}) \nabla u_i(t,\vec{r}))dV + \iiint_{V_i}S(\vec{r}) dV$

The only step left is to apply the divergence thereom to the first term on the rhs:

$\iiint_V \nabla \cdot \vec{f} dV = \iint_{\partial V} \nabla \vec{f} \cdot \hat{n} d \partial V$

which implies

${V_i}\frac{u_i^{n+1} - u_i^n}{\Delta t} = \iint_{\partial V_i} k(\vec{r}) \nabla u_i(t,\vec{r}) \cdot \hat{n}dA + \iiint_{V_i}S(\vec{r}) dV$

Which is your discretization (minus some trivial algebra).

Spatial Discretization

Now to actually attempt to answer your questions. The method is conditionally stable, with the condition depending on the answer to your second question. The important thing to note is that you don't need to know $u$ at the boundary but $\nabla u \cdot \hat{n}$ which is the gradient normal to the surface. If you choose your boundary between cells to be half way between the two points, then you can find

$\iint_{V_i} \nabla u \cdot \hat{n} dV = \Sigma_{j \in N_i} \frac{u_j - u_i}{||\Delta \vec{r}||} A_{i,j}$

where $N_i$ is the set of neighbors for node i, and $A_{i,j}$ is the area of their interface. So while the flux terms in the equation are written as at the one half nodes, you write that flux as a function of the values at the nodes. Of course this method can be changed, to increase both order of accuracy and stability of the numerical scheme. Some of these schemes require defining the boundary between cells differently as well. Depending on what scheme you choose you will have to look for its stability characteristics.