Computational Fluid Dynamics - Intermediate (b)

2-D Poisson's equation

This section focuses on the numerical approximation of the solution to a boundary value problem concerning a two-dimensional, steady-state Poisson's equation which is an elliptical PDE. To find the solution, the equation is discretized using a second-order centered finite difference scheme in both directions. An algorithm is used to calculate the results over different domains, then the quality and accuracy of our approximations are analyzed by comparing them to the exact solution of the equation and evaluating the resulting errors.

The equation and its boundary conditions are given as:
\begin{equation} \frac{d^2 u}{dx^2} + \frac{d^2 u}{dy^2} + 2\pi\sin(\pi x)\cos(\pi y) = 0 \end{equation} \begin{equation} \begin{split} (a) \hspace{5mm} (x,0)&=\sin (\pi x)\\ (b) \hspace{5mm} (x,1)&=-\sin (\pi x)\\ (c) \hspace{5mm}(0,y)&=0\\ (d) \hspace{5mm}(1,y)&=0\\ \end{split} \end{equation} Eq. (1) is to be discretized using a 2nd-order finite differences scheme and then solved over 2 different domains: Case 1, $ (x,y) \in [0,1]\times [0,1] $, and Case 2, $ (x,y) \in [-2,2]\times [-2,2] $. It should be noted that the boundary conditions given above are for case 1, and for case 2, it should be updated considering the exact solution. Due to the fact that the function \(u\) is dependent on 2 variables, the discretization results in 2 separate indices which need to be converted into a single index using an appropriate method. Then the problem can be represented in matrix-form as \([L]\{u\}=\{f\}\) and subsequently be solved by an algorithm. The results are to be plotted and compared to the exact solution of the PDE, given in the assignment as:
\begin{equation} u=\sin(\pi x) \cos(\pi y) \end{equation}

Discretization

The aim is to find an approximation of the solution to Eq. (1) using a numerical algorithm. To do this, first we have to define discrete points in the domain over which we intend to solve the problem (as an easy example $(x, y) \in [0,1] \times [0,1]$ as shown in Fig. b.1).

Fig. b.1: An example of the problem discretization.
Since there are 2 independent variables $x$ and $y$, each node $u_{i,j}$ needs to be assigned by 2 separate indices. Using a centered difference scheme, we have:
\begin{equation} \frac{u^j_{i-1} -2u^j_i +u^j_{i+1}}{\Delta x^2} + \frac{u^{j-1}_i -2u^j_i +u^{j+1}_i}{\Delta y^2} = - 2 \pi^2 \sin(\pi x) \cos (\pi y) \end{equation} This equation could already be turned into matrix-form $[L]\cdot \{u\}=\{f\}$, but in order to be able to use a flexible algorithm it is more convenient to have the nodes continuously numbered over the whole domain. For this, the following index transformation is employed (also shown in Fig. b.2):
\begin{equation} k=i+(j-1)I \end{equation}

Fig. b.2: Mapping the 2-indices to 1-index notation.
Where $I$ and $J$ are the total number of nodes along the $x$ (index $i$) and $y$ (index $j$) axis respectively. Applied to equation Eq. (4), this turns the left hand side into:
\begin{equation} \begin{split} \nabla^2 \vec u &= \frac{u_{k-1} - 2u_k + u_{k+1}}{\Delta x^2} + \frac{u_{k-I} - 2u_k + u_{k+I}}{\Delta y^2}\\ &=\frac{u_{k-1} + u_{k+1}}{\Delta x^2} + \frac{u_{k-I} + u_{k+I}}{\Delta y^2} - 2 u_k \left(\frac{1}{\Delta x^2} + \frac{1}{\Delta y^2}\right)\\ &= \beta (u_{k-1} + u_{k+1}) + \gamma (u_{k-I} + u_{k+I}) + \alpha u_k \end{split} \end{equation} in which:
\begin{equation} \alpha = - 2 \left(\frac{1}{\Delta x^2} + \frac{1}{\Delta y^2}\right),\hspace{5mm} \beta = \frac{1}{\Delta x^2},\hspace{5mm} \gamma = \frac{1}{\Delta y^2} \end{equation} Written the problem in a matrix form $[L]\cdot \{u\}=\{f\}$, the resulting scheme is as follows (the distance between the $\alpha$ and the corresponding $\gamma$ is equal to $I$):
\begin{equation} \begin{bmatrix} \alpha & \beta & 0 & 0 & 0 & \gamma & 0 & \cdots \\ \beta & \alpha & \beta & 0 & 0 & 0 & \gamma & \cdots \\ 0 & \beta & \alpha & \ddots & & & & \ddots \\ 0 & 0 & \ddots & \ddots \\ 0 & 0 \\ \gamma & 0\\ 0 & \gamma \\ \vdots & \vdots & \ddots \\ \end{bmatrix} \begin{pmatrix} u_1\\ u_2\\ u_3 \\ \vdots \\ \\ \\ \\ \vdots \\ u_n \end{pmatrix} = \begin{pmatrix} f_1 \\ f_2 \\ f_3 \\ \vdots \\ \\ \\ \\ \vdots \\ f_n \end{pmatrix} \end{equation}

Boundary conditions

To introduce the boundary conditions, the boundary nodes have to be recognized in the $ u $ vector and the corresponding rows in the $L$ matrix have to be replaced so that the value of each boundary node turns into $u_i = bc_i$. Taking this into consideration, the example with $I = 5$ becomes:
\begin{equation} \begin{matrix} &\left.\begin{matrix} \\ \\ I+1\text{ rows}\\ \text{replaced}\\ \\ \\ \end{matrix} \right\{ \\ &\left.\begin{matrix} I-2 \text{ rows}\\ \text{unchanged} \\ \\ \end{matrix} \right\{ \\ &\left.\begin{matrix} 2 \text{ rows}\\ \text{replaced} \end{matrix} \right\{ \\ &\left.\begin{matrix} I-2 \text{ rows}\\ \text{unchanged} \end{matrix} \right\{ \\ &\vdots \end{matrix} \begin{bmatrix} 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \\ 0& 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \cdots \\ 0 & 0 & 1 & 0 & 0 & 0 &0 & 0 & 0 & 0 & 0 & 0 & \\ 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0& 0 & 0\\ 0 & 0 & 0 & 0& 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & \gamma & 0 & 0 & 0 & \beta & \alpha & \beta & 0 & 0 & 0 & \gamma\\ 0 & 0 & \gamma & 0 & 0 & 0 & \beta & \alpha & \beta & 0 & 0 & 0 \\ 0 & 0 & 0 & \gamma & 0 & 0 & 0 & \beta & \alpha & \beta & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & \gamma & 0 & 0 & 0 & \beta & \alpha \\ & \vdots & & & & & & \ddots & & & & & \ddots \\ \\ \end{bmatrix} \begin{pmatrix} u_1 \\ u_2 \\ u_3 \\ u_4 \\ u_5 \\ u_6 \\ u_7 \\ u_8 \\ u_9 \\ u_{10} \\ u_{11} \\ u_{12} \\ \vdots \end{pmatrix} = \begin{pmatrix} bc_1 \\ bc_2 \\ bc_3 \\ bc_4 \\ bc_5 \\ bc_6 \\ f_7 \\ f_8 \\ f_9 \\ bc_{10} \\ bc_{11} \\ f_{12} \\ \vdots \\ \end{pmatrix} \end{equation}

Numerical and exact solutions

By plotting Eq. (3) over the 2-dimensional domains, the exact solutions for cases 1 and 2 can be obtained, which are shown in Fig. b.3 and b.4.

Fig. b.3: Numerical solution of case 1.

Fig. b.4: Numerical solution of case 2.

Validation of the code

To validate the code the root mean square error is plotted vs. the number of grid nodes. As shown in Fig. b.5, the scheme has a second order accuracy.

Fig. b.5: Second order convergence of the root mean square error.

The relevant Matlab codes can be downloaded using the following link:
Matlab Codes Bank

Special thanks to Julian Apfelthaler and Daniel Lange for their contributions to prepare this post.