Some Stable Solutions

Processing math: 0%

Some Stable Solutions

$\def \half {\frac{1}{2}} \def \slechts {\quad \Longleftrightarrow \quad}$ A Finite Difference scheme is associated with any (normed) tri-diagonal system of equations in the following manner:

$- b.T_{i-1} + T_i - a.T_{i+1} = 0$ The equations, when put in this form, can be solved in a direct fashion. Let's try the most simple of all possible solutions: a constant. Substitution of

$T_i = C$ gives:

$- b.C + C - a.C = 0 \quad \Longrightarrow \quad a + b = 1$ Assuming that

$C \ne 0$ (which would result in an even more trivial solution). It is seen that a constant solution can be obtained only if (minus) the sum of the off-diagonal elements is one.
The next trivial solution of the finite difference equation would be a straight line. Substitute

$T_i = C.i + D$ with

$D \ne 0$ and

$C \ne 0$ , giving:

$- b.\left[ C.(i-1) + D \right] + \left[ C.i + D \right] - a.\left[ C.(i+1) + D \right] =$

$C.i ( - b + 1 - a ) + D ( - b + 1 - a ) + ( b - a ) = 0$ Resulting in the same condition

$a + b = 1$ , augmented with

$b - a = 0$ . It is concluded that a {\em linear solution} can only be obtained for:

$a = b = \half$ Let another possible solution be given as:

$T_i = \left( \frac{b}{a} \right)^i$ And substitute. Then:

$- b \left( \frac{b}{a} \right)^{i-1} + \left( \frac{b}{a} \right)^i - a \left( \frac{b}{a} \right)^{i+1} = 0 \quad \Longrightarrow$

$( b / a )^i ( - b . a / b + 1 - a . b / a ) = ( b / a )^i ( - b + 1 - a ) = 0$ Resulting again in the condition

$a + b = 1$ . A more general solution may be cast into the form:

$T_i = C \left( \frac{b}{a} \right)^i + D$ Suppose that we have boundary conditions, like:

$T_1 = 1 \quad \mbox{and} \quad T_N = 0$ Then the constants

$C$ and

$D$ can be calculated explicitly, resulting in:

$T_i = \frac{ (b/a)^i - (b/a)^N }{ (b/a) - (b/a)^N } = \frac{ (b/a)^{i-1} - (b/a)^{N-1} }{ 1 - (b/a)^{N-1} }$ This result is more general than it seems at first sight, because it may be multiplied with an arbitrary constant and another arbitrary constant may be added to it, thus adapting it to quite arbitrary boundary conditions:

$T_i := (T_1 - T_N) . T_i + T_N$ It is recognized that the condition

$a + b = 1$ plays a rather predominant role in its relationship to the solutions which are found so far. The meaning of the condition is that the sum of the row coefficients in the (normed) matrix equals zero. Things are clarified further by writing the equations in a slightly different manner:

$T_i = b.T_{i-1} + a.T_{i+1}$ Knowing that

$b \ge 0$ ,

$a \ge 0$ and

$a+b = 1$ . It is concluded therefrom that

$T_i$ is a {\em weighted mean} of its neighbouring values, which means that it will always lie in between

$T_{i-1}$ and

$T_{i+1}$ :

$T_{i-1} \le T_{i} \le T_{i+1} \quad \mbox{or} \quad T_{i+1} \le T_{i} \le T_{i-1}$ Hence, irrespective of any boundary conditions, the function behaviour of the numerical solution will be continuously decreasing or continuously increasing. It shows no "wiggles". Such a solution

$T_i$ is commonly called {\em stable}. Thus, in fact, we are in search here for Stable Solutions of the tri-diagonal equations. It is remarked, in addition, that stability conditions, like

$a+b=1$ for the 1-D case, are well known and commonly accepted, as they belong to the "Four Basic Rules" in Patankar's book (1980).
Having established stability, our Stable numerical Solution is recalled:

$T_i = \frac{ (b/a)^{i-1} - (b/a)^{N-1} }{ 1 - (b/a)^{N-1} }$ Because, with help of the "Quotient Function", it can be cast now into the following form:

$T_i = \frac{ e^{P.(i-1).dx} - e^{P.(N-1).dx} }{ 1 - e^{P.(N-1).dx} } \slechts T(x) = \frac{ e^{P.x} - e^{P.L} }{ 1 - e^{P.L} }$ Where

$x = (i-1).dx$ is the 1-D coordinate and

$L$ is the total length of the mesh. Herewith, the general Stable (:

$a+b=1$ ) numerical Solution becomes quite alike an Analytical one:

$T(x) = (T_0 - T_L) \frac{ e^{P.x} - e^{P.L} }{ 1 - e^{P.L} } + T_L$ This solution is applicable to 1-D problems of Convection and Diffusion. It is concluded therefrom that any tri-diagonal system of linear equations, where the off-diagonal elements are negative and the rows sum up to zero, is actually a blueprint for Convection and Diffusion.
The number

$P.L$ in the exponents is recognized as the {\em P\'echlet number}. In real life, the dimensionless P\'echlet number has the following meaning:

$Pe = \rho.c.v.L/\lambda$ , where

$\rho =$ density

$(kg/m^3)$ ,

$c =$ heat capacity

$(J/kg/K)$ ,

$v =$ velocity

$(m/s)$ and

$\lambda =$ thermal conductivity

$(J/s/m/K)$ .
For the problem of 1-D Convection and Diffusion, the off-diagonal coefficients

$a$ and

$b$ can be written in a form which will be useful in the subsequent work. Define

$\epsilon$ by:

$\epsilon = b - a \qquad \mbox{while knowing} \qquad a + b = 1$ By addition and subtraction we find:

$b = \half (1 + \epsilon) \quad \mbox{and} \quad a = \half (1 - \epsilon) \quad \mbox{where} \quad -1 \le \epsilon \le +1$ It is remarked that the range of

$\epsilon$ is quite restricted in the first place. Furthermore it can be argued that, as the grid spacing becomes smaller and smaller with successive refinements of the grid, then

$\epsilon$ also must become smaller and smaller:

$\frac{1-\epsilon}{1+\epsilon} = \frac{a}{b}(dx \rightarrow 0) \rightarrow 1$ With help of a clever trick, we can express

$\epsilon$ solely in quotients of the off-diagonal coefficients. Divide namely nominator and denominator by

$\sqrt{a.b}$ in:

$\epsilon = \frac{b - a}{b + a} = \frac{\sqrt{b/a} - \sqrt{a/b}}{\sqrt{b/a} + \sqrt{a/b}} = \frac{e^{+P/2.dx}-e^{-P/2.dx}}{e^{+P/2.dx}+e^{-P/2.dx}}$ The latter expression is recognized as a Hyperbolic Tangent:

$tanh(P/2.dx)$ , as will be seen in one of the subsequent paragraphs. A sketch of the

$tanh$ function reveals that

$\epsilon$ will have the following properties:

$\lim_{P \rightarrow -\infty} \epsilon(P) = -1 \qquad \epsilon(0) = 0 \qquad \lim_{P \rightarrow +\infty} \epsilon(P) = +1$ When cast in Finite Element form, the analytical Convection-Diffusion scheme reads as follows:

$\left[ \begin{array}{cc} + a & - a \\ - b & + b \end{array} \right] = \left[ \begin{array}{cc} + \half - \half \epsilon & - \half + \half \epsilon \\ - \half - \half \epsilon & + \half + \half \epsilon \end{array} \right]$ The trace of any such element matrix is

$1$ . Double check for pure diffusion and for pure convection, respectively:

$P = 0 \quad \Longrightarrow \quad \left[ \begin{array}{cc} + a & - a \\ - b & + b \end{array} \right] = \left[ \begin{array}{cc} + 1/2 & - 1/2 \\ - 1/2 & + 1/2 \end{array} \right]$

$P = \infty \quad \Longrightarrow \quad \left[ \begin{array}{cc} + a & - a \\ - b & + b \end{array} \right] = \left[ \begin{array}{cc} 0 & 0 \\ - 1 & + 1 \end{array} \right]$

$P = - \infty \quad \Longrightarrow \quad \left[ \begin{array}{cc} + a & - a \\ - b & + b \end{array} \right] = \left[ \begin{array}{cc} + 1 & - 1 \\ 0 & 0 \end{array} \right]$ The latter expressions are demanding for a boundary condition at the inlets

$x=0$ and

$x=L$ respectively. Otherwise, the temperature

$T(x)$ at these places would be undefined.