Computational Fluid Dynamics - Intermediate (c)

2-D Helmholz equation

This section focuses on the numerical approximation of the solution to a boundary value problem concerning a two-dimensional, steady-state Helmholz' equation which is an elliptical PDE. Similar to 2D Poisson's equation, to find the solution, a second-order centered finite difference scheme in both directions is implemented.

The equation and its boundary conditions are given as:
\begin{equation} \frac{d^2 u}{dx^2} + \frac{d^2 u}{dy^2} + \sigma u = f(x,y) \end{equation} \begin{equation} \begin{split} (a) \hspace{5mm} u(x_{in},y)&=D_{left} \hspace{5mm} or \hspace{5mm} -\frac{d u}{dx}(x_{in},y)&=N_{left}\\ (b) \hspace{5mm} u(x,y_{in})&=D_{down} \hspace{5mm} or \hspace{5mm} -\frac{d u}{dy}(x,y_{in})&=N_{down}\\ (c) \hspace{5mm} u(x_{fin},y)&=D_{right} \hspace{5mm} or \hspace{5mm} \frac{d u}{dx}(x_{fin},y)&=N_{right}\\ (d) \hspace{5mm} u(x,y_{fin})&=D_{top} \hspace{5mm} or \hspace{5mm} \frac{d u}{dy}(x,y_{fin})&=N_{top}\\ \end{split} \end{equation} Eq. (1) will be numerically solved on a square domain $ (x,y) \in [x_{in},x_{fin}]\times [y_{in},y_{fin}] $ with a generic combination of Dirichlet's and Neumann's boundary conditions.

Discretization

Referring to the case discretization technique used for solving the 2D Poisson's equation, the centred difference scheme results:
\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} + \sigma u^{j}_i = f^{j}_i \end{equation} The indices mapping employed for the equation is:
\begin{equation} k=i+(j-1)I \end{equation} Where $ I $ is the total number of nodes along the $ x $ (index $ i $) axis. It is then possible to write the problem in a matrix form ( $ [L]\cdot \{u\}=\{f\} $ ) as follows:
\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} in which the distance between the $\alpha$ and the corresponding $\gamma$ is equal to $I$ and $\alpha$ , $\beta$ and $\gamma$ are:
\begin{equation} \alpha = - 2 \left(\frac{1}{\Delta x^2} + \frac{1}{\Delta y^2}\right) + \sigma,\hspace{5mm} \beta = \frac{1}{\Delta x^2},\hspace{5mm} \gamma = \frac{1}{\Delta y^2} \end{equation} The problem is finally completed including the boundary conditions. In the case of Neumann's boundary conditions, a first order scheme is used to discretize the first order derivative of $ u $. Such a choice allows not to introduce ghost nodes outside the domain, needful whether we wanted to continue using a centred scheme, but reduces the convergence order of the solver close to the Neumann's boundaries. Similar to what done for 2nd order ODEs, a forward decentred scheme is implemented on left and bottom boundaries and a backward decentred scheme is chosen for the right and top ones.

Validation

To validate the code the root mean square error is plotted vs. the number of nodes per edge. As shown in Fig. c.1, the scheme has a second order accuracy if we apply Dirichlet's boundary conditions for all the boundaries.

Fig. c.1: Second order convergence of the root mean square error for Dirichlet's boundary conditions.
On the other hand, when we apply Neumann's boundary conditions, it results evident the loose of convergence due to the first order scheme implemented just close to the boundaries. Fig. c.2 shows such a phenomenon considering the case in which left and bottom edges have Neumann's boundary conditions.

Fig. c.2: Loss of convergence for mixed Dirichlet's and Neumann's boundary conditions.

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