Burgers' equation: numerical solution - Dirichlet boundary conditions
In the previous section, the exact solution for the steady Burgers' equation is shown. In this chapter, the numerical solution is presented using Newton's method.
As it explained before, the 2-D viscous, steady Burgers' equation in conservative form can be written as:
\begin{equation} \begin{cases} 0.5 \partial_{x}u^2+v\partial_{y}u = \nu \left(\partial_{xx}u^2+\partial_{yy}u^2\right) \\ u\partial_{x}v+0.5 \partial_{y}v^2 = \nu \left(\partial_{xx}v^2+\partial_{yy}v^2\right) \end{cases} \end{equation} which are to be solved with Dirichlet boundary conditions on $ u $ and $ v $. For discretisation, the three-point centered difference scheme is employed using a uniform grid:
\begin{equation} \begin{cases} \frac{1}{2} \frac{u_{i+1,j}^2 - u_{i-1,j}^2}{2\Delta x} + v_{i,j} \frac{u_{i,j+1} - u_{i,j-1}}{2\Delta y} = \nu \left(\frac{u_{i+1,j} - 2u_{i,j} + u_{i-1,j}}{\Delta x^2} + \frac{u_{i,j+1} - 2u_{i,j} + u_{i,j-1}}{\Delta y^2}\right) \\ \\ u_{i,j} \frac{v_{i+1,j} - v_{i-1,j}}{2\Delta x} + \frac{1}{2} \frac{v_{i,j+1}^2 - v_{i,j-1}^2}{2\Delta y} = \nu \left(\frac{v_{i+1,j} - 2v_{i,j} + v_{i-1,j}}{\Delta x^2} + \frac{v_{i,j+1} - 2v_{i,j} + v_{i,j-1}}{\Delta y^2}\right) \end{cases} \end{equation} In order to simplify the code, it is useful to map the 2-indices notation to a 1-index notation: $ u_{i,j} \longleftrightarrow u_k $ in which:
\begin{equation} k = i + (j-1)I \end{equation} where $ I $ is the number of nodes in $ x $ direction. Then the equations can be written as:
\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] = 0 \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] = 0 \end{aligned} \end{cases} \end{equation} where:
\begin{equation} c_1 = \frac{1}{\Delta x} \ , \ c_2 = \frac{1}{\Delta y} \ , \ c_3 = \frac{1}{\Delta x^2} \ , \ c_4 = \frac{1}{\Delta y^2} \ , \ c_5 = c_3 + c_4. \end{equation} The boundary values of $ u $ and $ v $ are taken from the exact solution described before and the velocity values of all internal nodes, $ 2(I-2)(J-2) $, have to be solved. In order to solve the system of equations using Newton's method, the following equation should be iteratively solved:
\begin{equation} \Delta \boldsymbol{u}^{i+1} = \boldsymbol{u}^{i+1} - \boldsymbol{u}^i = - \boldsymbol{J}^{-1}(u^i) \boldsymbol{F}(u^{i}), \end{equation} where $ \boldsymbol{u}^i $ is the vector of solutions at $ i $-th iteration, $ \boldsymbol{F} $ contains $ F_k $ and $ G_k $, and $ \boldsymbol{J} $ is the Jacobian matrix which can be shown as follows:
Each block contains terms like (as an example for the first block):
\begin{equation} \frac{\partial F_k}{\partial u_{k-1}} \ , \ \frac{\partial F_k}{\partial u_k} \ , \ \frac{\partial F_k}{\partial u_{k+1}} \ , \ \frac{\partial F_k}{\partial u_{k-I}} \ , \ \frac{\partial F_k}{\partial u_{k+I}}, \end{equation} for all values of $ k $ and similar terms for other blocks.
All the formulation is implemented in the Matlab code which can be downloaded here.
The numerical results for case 1 together with the respective exact solution is plotted in Fig. a.2.1 and a.2.2. The grid size is $ 100 \times 100 $.
Similarly, the results for case 1 are illustrated in Fig. a.2.3 and a.2.4.
Convergence: Newton's method vs. numerical solution
It is a common mistake to mix up the convergence rate of Newton's method and the convergence of the numerical solution to the exact solution. When the Newton's method iteration converges (e.g. with the convergence rate of $ 10^{-8} $), it does not mean that the numerical solution error converges with the same rate. The numerical solution convergence also depends on the grid size. This is illustrated for case 1 in Fig. a.2.5. It is clear that the Newton's method (red curve) converges steadily to the defined convergence criteria ($ 10^{-8} $), but the numerical solution convergence (blue curve) saturates to a much larger error value ($ 2.3 \times 10^{-3} $) for grid size of $ 100 \times 100 $.
Implementing Neumann boundary conditions will be discussed in the next sub-section.
The relevant Matlab codes can be downloaded using the following link:
Matlab Codes Bank
No comments:
Post a Comment