Elliptic PDE finite volume method with Dirichlet boundary condition

Question

I want to discretize the following equation using a Finite Volume Method $$\nabla \cdot (a(x)\nabla u)=f(x)\\x\in \Omega \subset \mathbb{R}^2 \\u_{|\partial\Omega}=g$$ I'm using Voronoi cells here: let $V_i \subset\Omega$, $i=1,\dots,N$, such that $\bigcup\limits_{i=1}^N V_i=\Omega$, $V_i\cap V_j=\partial V_i \cap \partial V_j$ for $i\ne j$ be the mesh cells.

After applying the Divergence theorem and integrating the equation over each $V_i$, obtain

$$\int_{V_i} \nabla \cdot (a(x)\nabla u) dx = \int_{\partial V_i} a(x)\nabla u \cdot n ds=\int_{V_i} fdx$$

Let $x_i$ be the "center" of $V_i$ (this is a cell-cenetered method). Define $n_i\subset \{1,\dots,N\}$ as the set of indices of the cells immediately neighbouring $V_i$. Let the boundary between any two mesh cells be defined as

$$\Gamma_{i,j}:=\partial V_i \cap \partial V_j$$

Define the length of the boundary segment between neighbouring cells $V_i$ and $V_j$ as $$l_{i,j}:=\left|\Gamma_{i,j}\right|, j\in n_i$$

and the distance between two cell "centers" as $$h_{i,j}:=\|x_i-x_j\|_2, j\in n_i$$

Approximate $a(x)$ for each two neighbouring cells $V_i,V_j$ as

$$a_{i,j}:=a\left(\frac{x_i+x_j}{2}\right), j\in n_i$$

So we obtain the following discretization:

$\int_{\partial V_i} a(x)\nabla u \cdot n ds=\sum_{j\in n_i}\int_{\Gamma_{i,j}}a(x)\nabla u(x)\cdot n ds \approx \sum_{j\in n_i} a_{i,j} \frac{u_j-u_i}{h_{i,j}}l_{i,j}=|V_i|f(x_i)\approx \int_{V_i} f(x)dx$

To discretize the Dirichlet BC using centered cells, as far as I understand, one needs to use outer ghost cells. So let $V_{G_j}$ be such a cell, adjacent to $V_j$ for some $j\in\{1,\dots, N\}$, and $x_{G_{j}}\in V_{G_{j}}$ be its "center". But this is where I'm stuck. Where do I go from here?

P.S. I know there is a similar question here, but I couldn't quite follow the notation there.

Answer 1

In a cell centered finite volume setting your are discretising the fluxes over your cell boundary as you indicated. Lets say we look at one face between the cubes A and B, the calculated flux would be:

$$ f_{AB} = \frac{1}{\frac{1}{a_{A}} + \frac{1}{a_{B}} } ~\frac{u_B-u_A}{h} ~ F $$ with the fraction is the harmonic average (!) between the diffusivities in the two cells, the discretized gradient and $F$ being the area.

I don't know what software/code you are using, but whenever you are at a boundary face, simply swap the value u_B with the value you want to enforce at the boundary! Then the flux is just as if the border is another neighbour cell, with the distinction that it does not change. If you want to impose no flux conditions (i.e. du/dx = 0 at the border) then you set the whole flux to zero, as there should be no flux if the derivative is zero!