Burgers' equation: numerical solution - Dirichlet+Neumann boundary conditions

The numerical solution considering Dirichlet boundary conditions on all edges are explained in the previous sub-section. Here the numerical results for the same problem, but with a Neumann boundary condition on one edge (other edges remain the same) will be obtained. First the formulae derivation for Neumann boundary conditions will be explained. Then generating the solution using Newton's method will be discussed.

Formulae derivation for Neumann boundary conditions

As shown in Fig. a.3.1, the Neumann boundary conditions is implemented on the right edge of the 2-D domain. The same procedure can be employed to derive the Neumann boundary conditions on the other edges. Considering the right edge in Fig. a.3.1, the Neumann boundary conditions for values of velocity in $ x $ and $ y $ directions can be expressed as some functions of $ x $:
\begin{equation} \frac{\partial u}{\partial y} = g^u(x) \ , \ \frac{\partial v}{\partial y} = g^v(x) \end{equation} Before, applying Dirichlet boundary conditions, the values on the edge were given; but for Neumann boundary condition, the values on the edge are unknown and only their derivative in the normal direction to the edge are given.

In order to derive the formulae for Neumann boundary conditions, so called ghost nodes are added to the right of the edge. Then the 2nd order approximation of derivatives on the edge can be written as:
\begin{equation} \begin{cases} g^u(i,N) = \frac{u(i,N+1)-u(i,N-1)}{2 \Delta y} \\ \\ g^v(i,N) = \frac{v(i,N+1)-v(i,N-1)}{2 \Delta y} \end{cases} \end{equation}

Fig. a.3.1: Implementation of Neumann boundary condition on the right edge.
On the other hand, from the previous sub-section we have:
\begin{equation} \begin{cases} \!\begin{aligned} F_k = \frac{c_1}{4} & \left(u_{k+1}^2 - u_{k-1}^2\right) + \frac{c_2}{2} v_k \left(u_{k+I} - u_{k-I}\right) \\ & - \nu \left[ c_3 \left(u_{k+1} - u_{k-1}\right) + c_4 \left(u_{k+I} - u_{k-I}\right) - 2c_5 u_k \right] \end{aligned} \\ \\ \!\begin{aligned} G_k = \frac{c_1}{2} & u_k \left(v_{k+1} - v_{k-1}\right) + \frac{c_2}{4} \left(v_{k+I}^2 - v_{k-I}^2\right) \\ & - \nu \left[ c_3 \left(v_{k+1} - v_{k-1}\right) + c_4 \left(v_{k+I} - v_{k-I}\right) - 2c_5 v_k \right] \end{aligned} \end{cases} \end{equation} To be consistent with the previous notations, we assume $ I = N $. However, on the right edge, the values of $ u_{k+I} $ and $ v_{k+I} $ do not exist and have to be replaced. According to Eq. (2):
\begin{equation} u_{k+1} = u_{k-1} + 2 \Delta y \cdot g_k^u \ , \ v_{k+1} = v_{k-1} + 2 \Delta y \cdot g_k^v \end{equation} Substituting Eq. (4) into Eq. (3) results in:
\begin{equation} \begin{cases} \!\begin{aligned} F_{k,\text{right}} = \frac{c_1}{4} & \left(u_{k+1}^2 - u_{k-1}^2\right) + v_k g_k^u \\ & - \nu \left[ c_3 \left(u_{k+1} - u_{k-1}\right) + 2 c_4 \left(u_{k-I} + \Delta y g_k^u\right) - 2c_5 u_k \right] \end{aligned} \\ \\ \!\begin{aligned} G_{k,\text{right}} = \frac{c_1}{2} & u_k \left(v_{k+1} - v_{k-1}\right) + g_K^v \left(v_{k-I} + \Delta y g_k^v\right) \\ & - \nu \left[ c_3 \left(v_{k+1} - v_{k-1}\right) + 2 c_4 \left(v_{k-I} + \Delta y g_k^v\right) - 2c_5 v_k \right] \end{aligned} \end{cases} \end{equation}

Generating the solution using Newton's method

Expressions for $ F_{k,\text{right}} $ and $ G_{k,\text{right}} $ on the right edge are obtained in Eq. (5). The derivatives of the set of equations can be added to the Jacobian matrix of the Newton's method (see details here), which contains a set of equations' derivatives for the internal nodes.

The same procedure can be employed to apply Neumann boundary conditions on the other edges.

Using the provided Matlab code, the numerical results can be obtained for different grid sizes. For different grid sizes, it is recommended to modify the under-relaxation factor, $ \Omega $, if the code does not converge (line 115 in Matlab code). Fig. a.3.2 shows the exact together with the numerical results for case 1 and the grid size of $ 21 \times 21 $ and in Fig. a.3.3, the same but for case 2 and the grid size of $ 51 \times 51 $.