# One finite difference scheme

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?

• This actually looks like a fairly basic finite volume scheme (due to the divergence theorem usage). You have your new value equal to the old plus what flows over the boundaries plus the volumetric source. – Godric Seer Jan 24 '14 at 18:18
• Thank you, @GodricSeer. I had problems with translation of this term, so I could not find any information in global internet. Could you tell anything about stabitity of this scheme or about one-half-nodes? – cool Jan 24 '14 at 18:33

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.

• As I understood, $j$ is an index of a set of vertexes $V_{i,j,k}$ (there are 8 vertexes), $i$ are neighbours of vertex $j$ (there are 3 neighbours), $A_{i,j}$ is an area of that face that contains vertexes $i,j$ and $||r||$ is a distance between vertexes $i,j$, am I right? – cool Jan 25 '14 at 13:43
• I apologize, in your question, the equations are formulated on a 3D structured grid. So $i,j,k$ index into a this grid. In my question, i formulated it as unstructured, so $j$ is an index to a neighbor of $i$. $V_i$ is then the volume of cell $i$. – Godric Seer Jan 25 '14 at 15:44
• But there was no $V_i$ in the last formula for surface integral – cool Jan 25 '14 at 17:04
• I see i had an error, it was supposed to be dA not dV. That is where the $A_{i,j}$ comes from. – Godric Seer Jan 25 '14 at 20:06
• I didn't understand how we choose $A_{i,j}$ for vertexes $i,j$. Could you tell it once more? – cool Jan 28 '14 at 15:15