CFD Online Discussion Forums - Problem with Poisson equation with Neumann BC

Hi everybody !

In my CFD code, I use the pressure correction method and the package Fishpack to solve the pressure poisson equation. I tried to make my own solver for 2D Poisson equation with non homogeneous Neumann BC (FD scheme with

)

$\frac{\partial^2 p}{\partial x^2}+\frac{\partial^2 p}{\partial y^2}=f(x,y) \qquad \forall (x,y) \in [a,b]\times[c,d]$
$\frac{\partial p}{\partial x}=h \qquad \forall y \quad x=a,b$
$\frac{\partial p}{\partial y}=g \qquad \forall x \quad y=c,d$

I realized that it was not as easy as Dirichlet BC : I don't really know how to construct my matrix A. For example for a grid of 3x3 points, here is my matrix A : ( in $f_{i,j}$ i represtents x and j represents y)

$\begin{pmatrix} 1 & -\frac{1}{2} & 0 & -\frac{1}{2} & 0 & 0 & 0 & 0 & 0\\ -\frac{1}{2} & 2 & -\frac{1}{2} & 0 & -1 & 0 & 0 & 0 & 0\\ 0 & -\frac{1}{2} & 1 & 0 & 0 & -\frac{1}{2} & 0 & 0 & 0\\ -\frac{1}{2} & 0 & 0 & 2 & -1 & 0 & -\frac{1}{2} & 0 & 0\\ 0 & -1 & 0 & -1 & 4 & -1 & 0 & -1 & 0\\ 0 & 0 &-\frac{1}{2} &0 & -1 & 2 & 0 & 0 & -\frac{1}{2}\\ 0 & 0 & 0 & -\frac{1}{2} & 0 & 0 & 1 & -\frac{1}{2} & 0\\ 0 & 0 & 0 & 0 & -1 & 0 & -\frac{1}{2} & 2 & -\frac{1}{2}\\ 0 & 0 & 0 & 0 & 0 & -\frac{1}{2} & 0 & -\frac{1}{2} & 1\\ \end{pmatrix}$

My vector of unknowns p :
$\begin{pmatrix} p_{1,1}\\ p_{2,1}\\ p_{3,1}\\ p_{1,2}\\ p_{2,2}\\ p_{3,2}\\ p_{3,1}\\ p_{3,2}\\ p_{3,3}\\ \end{pmatrix}$

and the vector B :
$\begin{pmatrix} \frac{1}{4}(-dx^2f_{1,1}+g+h)\\ \frac{1}{2}(-dx^2f_{2,1}+g)\\ \frac{1}{4}(-dx^2f_{3,1}+g+h)\\ \frac{1}{2}(-dx^2f_{1,2}+h)\\ -dx^2f_{2,2}\\ \frac{1}{2}(-dx^2f_{3,2}+h)\\ \frac{1}{4}(-dx^2f_{3,1}+g+h)\\ \frac{1}{2}(-dx^2f_{3,2}+g)\\ \frac{1}{4}(-dx^2f_{3,3}+g+h)\\ \end{pmatrix}$

What do you think ? Is that correct ? What method should I use to solve this system ? SOR method ?