Computational Fluid Dynamics - Basic (b)



1-D time dependent Parabolic differential equations


Parabolic differential equations could be written as:
\begin{equation} \frac{\partial u(x,t)}{\partial t} + L[u(x,t)] = f(x,t), \ \ \ \ \ \ x \in [\alpha,\beta], \ t>0 \\ \end{equation}
where $ L $ is an elliptic operator applied to the solution, $ f $ is a known function, and $ u $ is the solution. Before implementing the finite difference schemes, it is important to consider some main features of the parabolic equations:
  - The solution admits just one characteristic line;
  - This line defines the dependence domain;
  - Along the parabolic direction $ t $ , just one “initial condition” is required;
  - Each node at the edge needs a boundary condition;

A generic example of such a problem is:
\begin{equation} \frac{\partial u(x,t)}{\partial t} = -c \frac{\partial u(x,t)}{\partial x} + \sigma \frac{\partial ^2 u(x,t)}{\partial x^2} + f(x,t); \end{equation} and the more specific case is the diffusion equation:
\begin{equation} \frac{\partial u(x,t)}{\partial t} = \sigma \frac{\partial ^2 u(x,t)}{\partial x^2} + f(x,t). \end{equation} Because of the importance of this equation in describing physical phenomena, in the following it is considered as a paradigm of parabolic equations.


Diffusion equation: finite difference discretization


Along the time direction, parabolicity of the equation suggests the use of forward scheme in order to consider the current solution as a function of just the previous time-steps.
Along the spatial direction, the homogeneous mathematical problem is symmetrical, so the best initial choice is to respect the symmetry; this way also guarantees the best convergence speed for constant $ \Delta x $ grid.
Referring to the previous chapter, a first order forward scheme is used for the unsteady term, and the second order centred scheme is employed for the diffusion term as follows:
\begin{equation} \frac{u^{n+1}-u^n}{\Delta t} = \sigma \frac{u_{i+1} - 2 u_i + u_{i-1}}{\Delta x^2} + f_i. \end{equation}

Explicit Euler Scheme


This scheme is said to be the most convenient finite difference scheme for discretizing the ODEs. The diffusion equation can be discretized using this scheme as follows:
\begin{equation} \frac{u_i^{n+1}-u_i^n}{\Delta t} = \sigma \frac{u_{i+1}^n - 2 u_i^n + u_{i-1}^n}{\Delta x^2} + f_i^n. \end{equation}
Fig. b.1. Explicit Euler scheme


Before writing the code for the discretization, it would be better to write down the discretized equation in matrix form:
\begin{equation} \left[ I \right] \left(\{u\}^{n+1} - \{u\}^{n} \right) = \left[L\right] \{u\}^n + \Delta t \{f\}^n \end{equation} where:
\begin{equation} [L] = \left( \begin{array}{ccccc} -2\nu & \nu & 0 & \ldots & 0 \\ \nu & -2\nu & \nu & \ldots & 0 \\ 0 & \nu & -2\nu & \ldots & 0 \\ \vdots & & & \ddots & \vdots \\ 0 & 0 & 0 & \ldots & -2\nu \end{array} \right) \ , \ \nu = \frac{\sigma \Delta t}{\Delta x^2} \end{equation} If $ [S] = [I] + [L] $, then we have:
\begin{equation} \{u\}^{n+1} = \left[S\right]\{u\}^{n} + \Delta t \{f\}^n \end{equation}

Two Dirichlet boundary conditions


Imposing Dirichlet boundary conditions at each boundary of the diffusion equation means that a value of the function is specified on that boundary of the domain. Considering the matrix form of the diffusion equation this type of boundary condition is described in the following.
\begin{equation} \left\lbrace\begin{array}{c} u_1\\ u_2\\ \vdots\\ u_i\\ \vdots\\ u_{N-1} \\ u_N \\ \end{array}\right\rbrace^{n+1} = \left( \begin{array}{ccccc} 1 & 0 & 0 & \ldots & 0 \\ & & & & \\ & & & & \\ & & [S] & & \\ & & & & \\ & & & & \\ & & & &\\ 0 & 0 & 0 & \ldots & 1 \\ \end{array} \right) \left\lbrace\begin{array}{c} u_{Dl}\\ u_2\\ \vdots\\ u_i\\ \vdots\\ u_{N-1}\\ u_{Dr} \\ \end{array}\right\rbrace^n + \Delta t\left\lbrace\begin{array}{c} 0\\ f_2\\ \vdots\\ f_i\\ \vdots\\ f_{N-1}\\ 0 \\ \end{array}\right\rbrace^n \end{equation} which can be also presented as:
\begin{equation} \begin{cases} u^{n+1}_1 = u^{n+1}_{Dl} \\ u^{n+1}_i - u^{n}_i = \frac{\sigma \Delta t}{\Delta x^2} \left({u^n_{i-1} -2 u^n_i + u^n_{i+1}}\right) + f^n_i\Delta t \\ u^{n+1}_N = u^{n+1}_{Dr} \end{cases} \end{equation} $ u_{Dl} $ and $ u_{Dr} $ are the value of Dirichlet boundary conditions at nodes $ x_1 $ and $ x_N $ respectively. The highlighted parts of the vectors and the matrix show the implementation of the Dirichlet boundary conditions. For each time-step $ n+1 $, the values of $ u $ are calculated using the given values of the previous time-step $ n $.

Neumann boundary condition on the left


In order to implement the Neumann boundary condition at node $ x_1 $, we have:
\begin{equation} \left.\frac{\partial u}{\partial t}\right|_{x=x_1} = \frac{u_0^{n+1} - u_2^{n+1}}{2\Delta x} = u_{Nl}; \end{equation} where $ u_{N1} $ is the value of Neumann boundary condition at node $ x_1 $. It means that a ghost node $ x_0 $ should be considered. However, using Eq. (5), $ u_2^{n+1} $ can be obtained:
\begin{equation} u_2^{n+1}-u_2^n = \nu \left(u_1^n - 2 u_2^n + u_3^n\right) + \Delta t f_2^n. \end{equation} Consequently, we have:
\begin{equation} u_0^{n+1} = \nu u_1^n + (1-2\nu) u_2^n + \nu u_3^n + 2 u_{Nl}^{n+1} \Delta x + f_2^n \Delta t. \end{equation} As a result, the Neumann boundary condition at $ x_1 $ can be shown as following:
\begin{equation} \left\lbrace\begin{array}{c} u_0\\ u_1\\ \vdots\\ u_i\\ \vdots\\ u_N \\ \end{array}\right\rbrace^{n+1} = \left( \begin{array}{cccccc} 0 & \nu & 1-2\nu & \nu & \ldots & 0 \\ & & & & & \\ & & & & & \\ & & [S] & & & \\ & & & & & \\ & & & & & \\ 0 & 0 & 0 & 0 & \ldots & 1 \\ \end{array} \right) \left\lbrace\begin{array}{c} u_0\\ u_1\\ \vdots\\ u_i\\ \vdots\\ u_{Dr} \\ \end{array}\right\rbrace^n + \Delta t\left\lbrace\begin{array}{c} 2\frac{\Delta x}{\Delta t}u_{Nl} + f_2\\ f_1\\ \vdots\\ f_i\\ \vdots\\ 0 \\ \end{array}\right\rbrace^n \end{equation}

Periodic boundary conditions


For the periodic boundary conditions, the last node is omitted from the solution vector. The following corrections on matrix $ S $ describes the periodic boundary conditions:
\begin{equation} [S] = \left( \begin{array}{cccccc} -2\nu+1 & \nu & 0 & \ldots & \nu & \not{0} \\ \nu & -2\nu+1 & \nu & \ldots & 0 & \not{0} \\ 0 & \nu & -2\nu+1 & \ldots & 0 & \not{0} \\ & & & & & \\ & & & & & \\ & & & & & \\ & & & & & \\ \nu & 0 & 0 & \ldots & -2\nu+1 & \not{\nu} \\ \not{0} & \not{0} & \not{0} & & \not{\nu} & \not{-}\not{2}\not{\nu}\not{+}\not{1}\\ \end{array} \right) \end{equation} As it is shown above, the last column and row are eliminated.


Implicit Euler scheme


This scheme is quite similar to the previous scheme, but differs in that it is an implicit method. The general form of this scheme is:
\begin{equation} \frac{u^{n+1}_i - u^{n}_i}{\Delta t} = \sigma \frac{u^{n+1}_{i-1} -2 u^{n+1}_i + u^{n+1}_{i+1}}{\Delta x^2} + f^{n+1}_i \end{equation}
Fig. b.2. Implicit Euler scheme


The discretized equation in matrix form can be written as:
\begin{equation} \left[ I \right] \left(\{u\}^{n+1} - \{u\}^{n} \right) = \left[L\right] \{u\}^{n+1} + \Delta t \{f\}^{n+1} \end{equation} where:
\begin{equation} [L] = \left( \begin{array}{ccccc} -2\nu & \nu & 0 & \ldots & 0 \\ \nu & -2\nu & \nu & \ldots & 0 \\ 0 & \nu & -2\nu & \ldots & 0 \\ \vdots & & & \ddots & \vdots \\ 0 & 0 & 0 & \ldots & -2\nu \end{array} \right) \ , \ \nu = \frac{\sigma \Delta t}{\Delta x^2} \end{equation} If $ [S] = [I] - [L] $, then we have:
\begin{equation} \{u\}^{n+1} = \left[S\right]^{-1}\left(\{u\}^{n} + \Delta t \{f\}^{n+1}\right) \end{equation}

Dirichlet boundary conditions


The Dirichlet boundary conditions on the left and right can be expressed as:
\begin{equation} \left( \begin{array}{ccccc} 1 & 0 & 0 & \ldots & 0 \\ & & & & \\ & & & & \\ & & [S] & & \\ & & & & \\ & & & & \\ & & & &\\ 0 & 0 & 0 & \ldots & 1 \\ \end{array} \right) \left\lbrace\begin{array}{c} u_1\\ u_2\\ \vdots\\ u_i\\ \vdots\\ u_{N-1} \\ u_N \\ \end{array}\right\rbrace^{n+1} = \left\lbrace\begin{array}{c} u_{Dl}\\ u_2\\ \vdots\\ u_i\\ \vdots\\ u_{N-1}\\ u_{Dr} \\ \end{array}\right\rbrace^n + \Delta t\left\lbrace\begin{array}{c} 0\\ f_2\\ \vdots\\ f_i\\ \vdots\\ f_{N-1}\\ 0 \\ \end{array}\right\rbrace^n \end{equation} which can be also presented as:
\begin{equation} \begin{cases} u^{n+1}_1 = u^{n+1}_{Dl} \\ u^{n+1}_i - u^{n}_i = \frac{\sigma \Delta t}{\Delta x^2} \left({u^{n+1}_{i-1} -2 u^{n+1}_i + u^{n+1}_{i+1}}\right) + f^{n+1}_i\Delta t \\ u^{n+1}_N = u^{n+1}_{Dr} \end{cases} \end{equation}

Neumann boundary condition on the left


In order to implement the Neumann boundary condition at node $ x_1 $, we have:
\begin{equation} \frac{u_0^{n+1} - u_2^{n+1}}{2\Delta x} = u^{n+1}_{Nl} \end{equation} and in matrix form:
\begin{equation} \left( \begin{array}{ccc} 0 & \ldots & 0 \\ & & \\ & & \\ [I] \\ & & \\ & &\\ & & \\ \end{array} \right) \left(\left\lbrace\begin{array}{c} u_0\\ u_1\\ \vdots\\ u_i\\ \vdots\\ u_N \\ \end{array}\right\rbrace^{n+1} - \left\lbrace\begin{array}{c} u_0\\ u_1\\ \vdots\\ u_i\\ \vdots\\ u_N \\ \end{array}\right\rbrace^{n} \right) = \left( \begin{array}{c} -1 \ \ 0 \ \ 1 \ \ \ldots \ \ 0 \\ \\ \\ [L] \\ \\ \\ \\ \end{array} \right) \left\lbrace\begin{array}{c} u_{0}\\ u_1\\ \vdots\\ u_i\\ \vdots\\ u_{N} \\ \end{array}\right\rbrace^{n+1} + \Delta t\left\lbrace\begin{array}{c} 2\frac{\Delta x}{\Delta t}u_{Nl}\\ f_1\\ \vdots\\ f_i\\ \vdots\\ f_{N} \\ \end{array}\right\rbrace^{n+1} \end{equation} It should be mentioned that in this case, the right boundary condition is not set.



Periodic boundary conditions


For the periodic boundary conditions, the last node is omitted from the solution vector. The following corrections on matrix $ S $ describes the periodic boundary conditions:
\begin{equation} [S] = \left( \begin{array}{cccccc} 2\nu+1 & -\nu & 0 & \ldots & -\nu & \not{0} \\ -\nu & 2\nu+1 & -\nu & \ldots & 0 & \not{0} \\ 0 & -\nu & 2\nu+1 & \ldots & 0 & \not{0} \\ & & & & & \\ & & & & & \\ & & & & & \\ & & & & & \\ -\nu & 0 & 0 & \ldots & 2\nu+1 & \not{-}\not{\nu} \\ \not{0} & \not{0} & \not{0} & & \not{-}\not{\nu} & \not{2}\not{\nu}\not{+}\not{1}\\ \end{array} \right) \end{equation}


Crank-Nicolson scheme


The Crank-Nicolson scheme can be presented as follows:
\begin{equation} \frac{u^{n+1}_i - u^{n}_i}{\Delta t} = \frac{1}{2}\sigma\left( \frac{u^n_{i-1} -2 u^n_i + u^n_{i+1}}{\Delta x^2} + f^n_i\right) + \frac{1}{2}\sigma \left(\frac{u^{n+1}_{i-1} -2 u^{n+1}_i + u^{n+1}_{i+1}}{\Delta x^2} + f^{n+1}_i\right) \end{equation}
Fig. b.3. Crank-Nicolson scheme


In order to write the method in matrix form, it is better to divide the explicit ($ EE $) and implicit ($ IE $) terms:
\begin{equation} u^{n+1}_i - u^{n}_i = \frac{1}{2} rhs_{EE} + \frac{1}{2} rhs_{IE} \end{equation} So we have:
\begin{equation} \left[ I \right] \left(\{u\}^{n+1} - \{u\}^{n} \right) = \frac{1}{2}\left[L_{EE}\right] \{u\}^{n} + \frac{1}{2}\left[L_{IE}\right] \{u\}^{n+1} + \Delta t \left(\frac{1}{2}\{f\}^{n}+\frac{1}{2}\{f\}^{n+1}\right) \end{equation} If $ [S_{EE}] = [I] + \frac{1}{2}[L_{EE}] $ and $ [S_{IE}] = [I] - \frac{1}{2}[L_{IE}] $ we can rewrite the Eq. 27 as:
\begin{equation} \{u\}^{n+1} = \left[S_{IE}\right]^{-1} \left(\left[S_{EE}\right]\{u\}^{n} + \frac{\Delta t}{2} \left(\{f\}^{n}+\{f\}^{n+1}\right)\right) \end{equation}

Dirichlet boundary conditions



\begin{equation} \left( \begin{array}{ccccc} 1 & 0 & 0 & \ldots & 0 \\ & & & & \\ & & & & \\ & & [S_{IE}] & & \\ & & & & \\ & & & & \\ & & & &\\ 0 & 0 & 0 & \ldots & 1 \\ \end{array} \right) \left\lbrace\begin{array}{c} u_1\\ u_2\\ \vdots\\ u_i\\ \vdots\\ u_{N-1} \\ u_N \\ \end{array}\right\rbrace^{n+1} = \left( \begin{array}{ccccc} 0 & 0 & 0 & \ldots & 0 \\ & & & & \\ & & & & \\ & & [S_{EE}] & & \\ & & & & \\ & & & & \\ & & & &\\ 0 & 0 & 0 & \ldots & 0 \\ \end{array} \right)\left\lbrace\begin{array}{c} u_1\\ u_2\\ \vdots\\ u_i\\ \vdots\\ u_{N-1}\\ u_N \\ \end{array}\right\rbrace^n + \frac{\Delta t}{2}\left(\left\lbrace\begin{array}{c} \frac{2 u_{Dl}}{\Delta t}\\ f_2\\ \vdots\\ f_i\\ \vdots\\ f_{N-1}\\ \frac{2 u_{Dr}}{\Delta t} \\ \end{array} \right\rbrace^n + \left\lbrace\begin{array}{c} 0\\ f_2\\ \vdots\\ f_i\\ \vdots\\ f_{N-1}\\ 0 \\ \end{array}\right\rbrace^{n+1}\right) \end{equation}

Neumann boundary condition on the left


\begin{equation} \left( \begin{array}{ccccc} 1 & 0 & -1 & \ldots & 0 \\ & & & & \\ & & & & \\ & & [S_{IE}] & & \\ & & & & \\ & & & & \\ & & & &\\ & & & & \\ \end{array} \right) \left\lbrace\begin{array}{c} u_0\\ u_1\\ \vdots\\ u_i\\ \vdots\\ u_{N-1} \\ u_N \\ \end{array}\right\rbrace^{n+1} = \left( \begin{array}{ccccc} 0 & 0 & 0 & \ldots & 0 \\ & & & & \\ & & & & \\ & & [S_{EE}] & & \\ & & & & \\ & & & & \\ & & & &\\ & & & & \\ \end{array} \right)\left\lbrace\begin{array}{c} u_0\\ u_1\\ \vdots\\ u_i\\ \vdots\\ u_{N-1}\\ u_N \\ \end{array}\right\rbrace^n + \frac{\Delta t}{2}\left(\left\lbrace\begin{array}{c} \frac{4\Delta x u_{Nl}}{\Delta t}\\ f_1\\ \vdots\\ f_i\\ \vdots\\ f_{N-1}\\ f_{N} \\ \end{array} \right\rbrace^n + \left\lbrace\begin{array}{c} 0\\ f_1\\ \vdots\\ f_i\\ \vdots\\ f_{N-1}\\ f_{N} \\ \end{array}\right\rbrace^{n+1}\right) \end{equation}

Periodic boundary conditions


\begin{equation} [S_{IE}] = \left( \begin{array}{cccccc} 2\nu+1 & -\nu & 0 & \ldots & -\nu & \not{0} \\ -\nu & 2\nu+1 & -\nu & \ldots & 0 & \not{0} \\ 0 & -\nu & 2\nu+1 & \ldots & 0 & \not{0} \\ & & & & & \\ & & & & & \\ & & & & & \\ & & & & & \\ -\nu & 0 & 0 & \ldots & 2\nu+1 & \not{-}\not{\nu} \\ \not{0} & \not{0} & \not{0} & & \not{-}\not{\nu} & \not{2}\not{\nu}\not{+}\not{1}\\ \end{array} \right) \end{equation}
\begin{equation} [S_{EE}] = \left( \begin{array}{cccccc} -2\nu+1 & \nu & 0 & \ldots & \nu & \not{0} \\ \nu & -2\nu+1 & \nu & \ldots & 0 & \not{0} \\ 0 & \nu & -2\nu+1 & \ldots & 0 & \not{0} \\ & & & & & \\ & & & & & \\ & & & & & \\ & & & & & \\ \nu & 0 & 0 & \ldots & -2\nu+1 & \not{\nu} \\ \not{0} & \not{0} & \not{0} & & \not{\nu} & \not{-}\not{2}\not{\nu}\not{+}\not{1}\\ \end{array} \right) \end{equation}



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

No comments:

Post a Comment