# One-Dimensional Steady-State Convection and Diffusion

(Difference between revisions)
 Revision as of 05:39, 21 July 2010 (view source)← Older edit Current revision as of 02:11, 27 July 2010 (view source) (14 intermediate revisions not shown) Line 1: Line 1: {{Comp Method for Forced Convection Category}} {{Comp Method for Forced Convection Category}} - Exact Solution + ==Exact Solution== - The objective of this subsection is to introduce various discretization schemes of the convection-diffusion terms through discussion of the one-dimensional steady state convection and diffusion problem. For a one-dimensional steady-state convection and diffusion problem, the governing equation is + The objective of this article is to introduce various discretization schemes of the convection-diffusion terms through discussion of the one-dimensional steady state convection and diffusion problem. For a one-dimensional steady-state convection and diffusion problem, the governing equation is {| class="wikitable" border="0" {| class="wikitable" border="0" Line 90: Line 90: which will be used as a criterion to check the accuracy of various discretization schemes. which will be used as a criterion to check the accuracy of various discretization schemes. - [[Central Difference Scheme]]
+ ==Central Difference Scheme== - Integrating [[One-Dimensional_Steady-State_Convection_and_Diffusion#equation_.282.29|the governing equation]] over the control volume P (shaded area in the figure to the right), one obtains + Integrating eq. (1) over the control volume P (shaded area in the figure to the right), one obtains [[Image:Fig4.17.png|thumb|400 px|alt=Control volume for one-dimensional problem | Control volume for one-dimensional problem.]] [[Image:Fig4.17.png|thumb|400 px|alt=Control volume for one-dimensional problem | Control volume for one-dimensional problem.]] Line 101: Line 101: - |{{EquationRef|(1)}} + |{{EquationRef|(11)}} |} |} - The right-hand side of eq. (1) can be obtained by assuming the distribution of $\varphi$ between any two neighboring grid points is piecewise linear, i.e., + The right-hand side of eq. (11) can be obtained by assuming the distribution of $\varphi$ between any two neighboring grid points is piecewise linear, i.e.,
$\left( \Gamma \frac{d\varphi }{dx} \right)_{e}=\Gamma _{e}\frac{\varphi _{E}-\varphi _{P}}{(\delta x)_{e}}$
$\left( \Gamma \frac{d\varphi }{dx} \right)_{e}=\Gamma _{e}\frac{\varphi _{E}-\varphi _{P}}{(\delta x)_{e}}$
Line 110: Line 110:
$\left( \Gamma \frac{d\varphi }{dx} \right)_{w}=\Gamma _{w}\frac{\varphi _{P}-\varphi _{W}}{(\delta x)_{w}}$
$\left( \Gamma \frac{d\varphi }{dx} \right)_{w}=\Gamma _{w}\frac{\varphi _{P}-\varphi _{W}}{(\delta x)_{w}}$
- where Γe and Γw are the diffusivities at the faces of the control volume. To ensure that the flux of $\varphi$ across the faces of the control volume is continuous, the harmonic mean diffusivity at the faces should be used. To evaluate the left hand side of eq. (1), it is necessary to know the values of $\varphi$ at the faces of the control volume. If the piecewise linear profile of $\varphi$ is chosen, it follows that + where Γe and Γw are the diffusivities at the faces of the control volume. To ensure that the flux of $\varphi$ across the faces of the control volume is continuous, the harmonic mean diffusivity at the faces should be used. To evaluate the left hand side of eq. (11), it is necessary to know the values of $\varphi$ at the faces of the control volume. If the piecewise linear profile of $\varphi$ is chosen, it follows that
$\varphi _{e}=\frac{\varphi _{E}+\varphi _{P}}{2}$
$\varphi _{e}=\frac{\varphi _{E}+\varphi _{P}}{2}$
Line 117: Line 117:
$\varphi _{w}=\frac{\varphi _{P}+\varphi _{W}}{2}$
$\varphi _{w}=\frac{\varphi _{P}+\varphi _{W}}{2}$
- Therefore, eq. (1) becomes + Therefore, eq. (11) becomes {| class="wikitable" border="0" {| class="wikitable" border="0" Line 125: Line 125: - |{{EquationRef|(2)}} + |{{EquationRef|(12)}} |} |} Defining the mass flux and diffusive conductance Defining the mass flux and diffusive conductance Line 134: Line 134: $F=\rho u,\text{ }D=\frac{\Gamma }{\delta x}$ $F=\rho u,\text{ }D=\frac{\Gamma }{\delta x}$ - |{{EquationRef|(3)}} + |{{EquationRef|(13)}} |} |} - eq. (2) can be rearranged as + eq. (12) can be rearranged as {| class="wikitable" border="0" {| class="wikitable" border="0" Line 146: Line 146: - |{{EquationRef|(4)}} + |{{EquationRef|(14)}} |} |} where where Line 155: Line 155: $a_{E}=D_{e}-\frac{1}{2}F_{e}$ $a_{E}=D_{e}-\frac{1}{2}F_{e}$ - |{{EquationRef|(5)}} + |{{EquationRef|(15)}} |} |} Line 163: Line 163: $a_{W}=D_{w}+\frac{1}{2}F_{w}$ $a_{W}=D_{w}+\frac{1}{2}F_{w}$ - |{{EquationRef|(6)}} + |{{EquationRef|(16)}} |} |} Line 171: Line 171: $\begin{matrix}{}\\\end{matrix}a_{P}=a_{W}+a_{E}+(F_{e}-F_{w})$ $\begin{matrix}{}\\\end{matrix}a_{P}=a_{W}+a_{E}+(F_{e}-F_{w})$ - |{{EquationRef|(7)}} + |{{EquationRef|(17)}} |} |} - This scheme is termed the central difference scheme because the values of $\varphi$ at the faces of the control volume are taken as the averaged value between two grid points. The continuity equation requires that $F_{e}=F_{w}$ and therefore, eq. (7) reduces to + This scheme is termed the central difference scheme because the values of $\varphi$ at the faces of the control volume are taken as the averaged value between two grid points. The continuity equation requires that $F_{e}=F_{w}$ and therefore, eq. (17) reduces to $a_{P}=a_{W}+a_{E}$ $a_{P}=a_{W}+a_{E}$ - To evaluate the performance of the central difference scheme, let us consider the case of a uniform grid, i.e., $(\delta x)_{e}=(\delta x)_{w}=\delta x$, for which case eq.  (2) can be rearranged as + To evaluate the performance of the central difference scheme, let us consider the case of a uniform grid, i.e., $(\delta x)_{e}=(\delta x)_{w}=\delta x$, for which case eq.  (12) can be rearranged as {| class="wikitable" border="0" {| class="wikitable" border="0" Line 184: Line 184: $\varphi _{P}=\frac{1}{2}\left[ \left( 1-\frac{\text{Pe}_{\Delta }}{2} \right)\varphi _{E}+\left( 1+\frac{\text{Pe}_{\Delta }}{2} \right)\varphi _{W} \right]$ $\varphi _{P}=\frac{1}{2}\left[ \left( 1-\frac{\text{Pe}_{\Delta }}{2} \right)\varphi _{E}+\left( 1+\frac{\text{Pe}_{\Delta }}{2} \right)\varphi _{W} \right]$ - |{{EquationRef|(8)}} + |{{EquationRef|(18)}} |} |} where where Line 193: Line 193: $\text{Pe}_{\Delta }=\frac{\rho u\delta x}{\Gamma }=\frac{F}{D}$ $\text{Pe}_{\Delta }=\frac{\rho u\delta x}{\Gamma }=\frac{F}{D}$ - |{{EquationRef|(9)}} + |{{EquationRef|(19)}} |} |} is the Peclet number using grid size as the characteristic length, which is referred to as the grid Peclet number. The grid Pe is a ratio of the strength of convection over diffusion. To ensure stability of the discretization scheme, the value of $\varphi _{P}$ should always fall between $\varphi _{E}$ and $\varphi _{W}$, which requires that the coefficients, $\varphi _{E}$ and $\varphi _{W}$, are positive, i.e., is the Peclet number using grid size as the characteristic length, which is referred to as the grid Peclet number. The grid Pe is a ratio of the strength of convection over diffusion. To ensure stability of the discretization scheme, the value of $\varphi _{P}$ should always fall between $\varphi _{E}$ and $\varphi _{W}$, which requires that the coefficients, $\varphi _{E}$ and $\varphi _{W}$, are positive, i.e., Line 202: Line 202: $\left| \text{Pe}_{\Delta } \right|\le 2$ $\left| \text{Pe}_{\Delta } \right|\le 2$ - |{{EquationRef|(10)}} + |{{EquationRef|(20)}} |} |} - This is the criterion for stability of the central difference scheme. It can be demonstrated that the central difference becomes unstable if eq. (10) is violated. The fact that the central difference scheme is stable under small grid Peclet number indicates that the central difference scheme is accurate only if the convection is not very significant. + This is the criterion for stability of the central difference scheme. It can be demonstrated that the central difference becomes unstable if eq. (20) is violated. The fact that the central difference scheme is stable under small grid Peclet number indicates that the central difference scheme is accurate only if the convection is not very significant. - [[Upwind Scheme]]
+ ==Upwind Scheme== - The central difference scheme assumes that the effects of the values of $\varphi$ at two neighboring grid points on the value of $\varphi$ at the face of the control volume are equal. This assumption is valid only if the effect of diffusion is dominant. If, on the other hand, the convection is dominant, one can expect that the effect of the grid point upwind is more significant than that of the point downwind. If we can assume that the value of $\varphi$ at the face of the control volume is dominated by the value of $\varphi$ at the grid point at the upwind side and that the effect of the value of $\varphi$ at the downwind side can be neglected, the two terms on the left hand side of eq. (4.211) can be expressed as + The central difference scheme assumes that the effects of the values of $\varphi$ at two neighboring grid points on the value of $\varphi$ at the face of the control volume are equal. This assumption is valid only if the effect of diffusion is dominant. If, on the other hand, the convection is dominant, one can expect that the effect of the grid point upwind is more significant than that of the point downwind. If we can assume that the value of $\varphi$ at the face of the control volume is dominated by the value of $\varphi$ at the grid point at the upwind side and that the effect of the value of $\varphi$ at the downwind side can be neglected, the two terms on the left hand side of eq. (11) can be expressed as - $(\rho u\varphi )_{e}=\left\{ \begin{matrix} + [itex](\rho u\varphi )_{e}=\left\{ \begin{matrix} F_{e}\varphi _{P},\text{ }F_{e}>0 \\ F_{e}\varphi _{P},\text{ }F_{e}>0 \\ F_{e}\varphi _{E},\text{ }F_{e}<0 \\ F_{e}\varphi _{E},\text{ }F_{e}<0 \\ - \end{matrix} \right.$ + \end{matrix} \right.
- $(\rho u\varphi )_{w}=\left\{ \begin{matrix} + [itex](\rho u\varphi )_{w}=\left\{ \begin{matrix} F_{w}\varphi _{W},\text{ }F_{w}>0 \\ F_{w}\varphi _{W},\text{ }F_{w}>0 \\ F_{w}\varphi _{P},\text{ }F_{w}<0 \\ F_{w}\varphi _{P},\text{ }F_{w}<0 \\ - \end{matrix} \right.$ + \end{matrix} \right.
The above two equations can be expressed in the following compact form: The above two equations can be expressed in the following compact form: - $(\rho u\varphi )_{e}=\varphi _{P}\left[\!\left[ F_{e},0 \right]\!\right]-\varphi _{E}\left[\!\left[ -F_{e},0 \right]\!\right]$ +
$(\rho u\varphi )_{e}=\varphi _{P}\left[\!\left[ F_{e},0 \right]\!\right]-\varphi _{E}\left[\!\left[ -F_{e},0 \right]\!\right] - [itex](\rho u\varphi )_{w}=\varphi _{W}\left[\!\left[ F_{w},0 \right]\!\right]-\varphi _{P}\left[\!\left[ -F_{w},0 \right]\!\right]$ +
$(\rho u\varphi )_{w}=\varphi _{W}\left[\!\left[ F_{w},0 \right]\!\right]-\varphi _{P}\left[\!\left[ -F_{w},0 \right]\!\right] - where the operator [itex]\left[\!\left[ A,B \right]\!\right]$ denotes the greater of A and B (Patankar, 1980). Substituting the above expression into the left hand side of eq. (4.211) and using central difference for the right hand side of eq. (4.211), the discretized equation becomes + where the operator $\left[\!\left[ A,B \right]\!\right]$ denotes the greater of ''A'' and ''B'' Patankar, S.V., 1980, Numerical Heat Transfer and Fluid Flow, Hemisphere, Washington, DC. . Substituting the above expression into the left hand side of eq. (11) and using central difference for the right hand side of eq. (11), the discretized equation becomes {| class="wikitable" border="0" {| class="wikitable" border="0" Line 236: Line 236: - |{{EquationRef|(1)}} + |{{EquationRef|(21)}} |} |} where where Line 245: Line 245: $a_{E}=D_{e}+\left[\!\left[ -F_{e},0 \right]\!\right]$ $a_{E}=D_{e}+\left[\!\left[ -F_{e},0 \right]\!\right]$ - |{{EquationRef|(2)}} + |{{EquationRef|(22)}} |} |} Line 253: Line 253: $a_{W}=D_{w}+\left[\!\left[ F_{w},0 \right]\!\right]$ $a_{W}=D_{w}+\left[\!\left[ F_{w},0 \right]\!\right]$ - |{{EquationRef|(3)}} + |{{EquationRef|(23)}} |} |} Line 261: Line 261: $\begin{matrix}{}\\\end{matrix}a_{P}=a_{W}+a_{E}+(F_{e}-F_{w})$ $\begin{matrix}{}\\\end{matrix}a_{P}=a_{W}+a_{E}+(F_{e}-F_{w})$ - |{{EquationRef|(4)}} + |{{EquationRef|(24)}} |} |} - The above scheme is referred to as the upwind scheme because the value of $\varphi$ at the grid point on the upwind side was used as the value of $\varphi$ at the face of the control volume to discretize the convection term. The upwind scheme ensures that the coefficients in eq. (4.221) are always positive so that a physically unrealistic solution can be avoided. + The above scheme is referred to as the upwind scheme because the value of $\varphi$ at the grid point on the upwind side was used as the value of $\varphi$ at the face of the control volume to discretize the convection term. The upwind scheme ensures that the coefficients in eq. (21) are always positive so that a physically unrealistic solution can be avoided. - [[Hybrid Scheme]]
+ ==Hybrid Scheme== - The upwind scheme uses the value of $\varphi$ from the grid point at the upwind side as the value of $\varphi$ at the face of the control volume regardless of the grid Peclet number. While this treatment can yield accurate results for cases with high Peclet number, the result will not be accurate for cases where the grid Peclet number is near zero; for which cases the central difference scheme can produce better results. Spalding (1972) proposed a hybrid scheme that uses the central difference scheme when $\left| \text{Pe}_{\Delta } \right|\le 2$ and the upwind scheme when $\left| \text{Pe}_{\Delta } \right|>2$. + The upwind scheme uses the value of $\varphi$ from the grid point at the upwind side as the value of $\varphi$ at the face of the control volume regardless of the grid Peclet number. While this treatment can yield accurate results for cases with high Peclet number, the result will not be accurate for cases where the grid Peclet number is near zero; for which cases the central difference scheme can produce better results. Spalding Spalding, D.B., 1972, “A Novel Finite-difference Formulation for Differential Expressions Involving both First and Second Derivatives,” Int. J. Num. Methods Eng., Vol. 4, pp. 551-559. proposed a hybrid scheme that uses the central difference scheme when $\left| \text{Pe}_{\Delta } \right|\le 2$ and the upwind scheme when $\left| \text{Pe}_{\Delta } \right|>2$. - To observe the difference between the central difference and upwind schemes, the coefficient for the east neighboring grid point, eqs. (4.215) and (4.222), can be rewritten as + To observe the difference between the central difference and upwind schemes, the coefficient for the east neighboring grid point, eqs. (15) and (22), can be rewritten as {| class="wikitable" border="0" {| class="wikitable" border="0" |- |- Line 274: Line 274: $a_{E}/D_{e}=1-\frac{1}{2}\text{Pe}_{\Delta e},\text{ Central difference scheme}$ $a_{E}/D_{e}=1-\frac{1}{2}\text{Pe}_{\Delta e},\text{ Central difference scheme}$ - |{{EquationRef|(1)}} + |{{EquationRef|(25)}} |} |} Line 282: Line 282: $a_{E}/D_{e}=1+\left[\!\left[ -\text{Pe}_{\Delta e},0 \right]\!\right],\text{ Upwind scheme}$ $a_{E}/D_{e}=1+\left[\!\left[ -\text{Pe}_{\Delta e},0 \right]\!\right],\text{ Upwind scheme}$ - |{{EquationRef|(2)}} + |{{EquationRef|(26)}} |} |} The hybrid scheme can then be expressed as The hybrid scheme can then be expressed as - $a_{E}/D_{e}=\left\{ \begin{matrix} + [itex]a_{E}/D_{e}=\left\{ \begin{matrix} -\text{Pe}_{\Delta e}\text{ Pe}_{\Delta e}<-2 \\ -\text{Pe}_{\Delta e}\text{ Pe}_{\Delta e}<-2 \\ 1-\frac{1}{2}\text{Pe}_{\Delta e}\text{ }-\text{2}\le \text{ Pe}_{\Delta e}\le 2 \\ 1-\frac{1}{2}\text{Pe}_{\Delta e}\text{ }-\text{2}\le \text{ Pe}_{\Delta e}\le 2 \\ 0\text{ Pe}_{\Delta e}>2 \\ 0\text{ Pe}_{\Delta e}>2 \\ - \end{matrix} \right.$ + \end{matrix} \right.
which can be rewritten in the following compact form which can be rewritten in the following compact form Line 299: Line 299: $a_{E}/D_{e}=\left[\!\left[ -\text{Pe}_{\Delta e},1-\frac{1}{2}\text{Pe}_{\Delta e},0 \right]\!\right]$ $a_{E}/D_{e}=\left[\!\left[ -\text{Pe}_{\Delta e},1-\frac{1}{2}\text{Pe}_{\Delta e},0 \right]\!\right]$ - |{{EquationRef|(3)}} + |{{EquationRef|(27)}} |} |} The coefficient for the west neighbor grid point can be obtained using a similar approach. The coefficient for the west neighbor grid point can be obtained using a similar approach. Line 308: Line 308: $a_{W}/D_{w}=\left[\!\left[ \text{Pe}_{\Delta w},1+\frac{1}{2}\text{Pe}_{\Delta w},0 \right]\!\right]$ $a_{W}/D_{w}=\left[\!\left[ \text{Pe}_{\Delta w},1+\frac{1}{2}\text{Pe}_{\Delta w},0 \right]\!\right]$ - |{{EquationRef|(4)}} + |{{EquationRef|(28)}} |} |} - The above hybrid scheme combines the advantages of the central difference and upwind schemes to yield better results for cases where $\left| \text{Pe}_{\Delta } \right|\to \infty$ or $\left| \text{Pe}_{\Delta } \right|\sim 0$. However, there is still room for improvement of the solution when $\left| \text{Pe}_{\Delta } \right|$ is near 2 (see Problem 4.23). + The above hybrid scheme combines the advantages of the central difference and upwind schemes to yield better results for cases where $\left| \text{Pe}_{\Delta } \right|\to \infty$ or $\left| \text{Pe}_{\Delta } \right|\sim 0$. However, there is still room for improvement of the solution when $\left| \text{Pe}_{\Delta } \right|$ is near 2. - [[Exponential and Power Law Schemes]]
+ ==Exponential and Power Law Schemes== - Since the exact solution of eq. (4.201) exists, one can reasonably expect that an  accurate scheme can be derived if the result of the exact solution, eq. (4.210), is utilized. Equation (4.201) can be rewritten as + Since the exact solution of eq. (1) exists, one can reasonably expect that an  accurate scheme can be derived if the result of the exact solution, eq. (10), is utilized. Equation (1) can be rewritten as {| class="wikitable" border="0" {| class="wikitable" border="0" Line 320: Line 320: $\frac{d}{dx}\left( \rho u\varphi -\Gamma \frac{d\varphi }{dx} \right)=0$ $\frac{d}{dx}\left( \rho u\varphi -\Gamma \frac{d\varphi }{dx} \right)=0$ - |{{EquationRef|(1)}} + |{{EquationRef|(29)}} |} |} Defining the total flux of $\varphi$ due to convection and diffusion Defining the total flux of $\varphi$ due to convection and diffusion Line 329: Line 329: $J=\rho u\varphi -\Gamma \frac{d\varphi }{dx}$ $J=\rho u\varphi -\Gamma \frac{d\varphi }{dx}$ - |{{EquationRef|(2)}} + |{{EquationRef|(30)}} |} |} - eq. (4.229) becomes + eq. (29) becomes {| class="wikitable" border="0" {| class="wikitable" border="0" Line 338: Line 338: $\frac{dJ}{dx}=0$ $\frac{dJ}{dx}=0$ - |{{EquationRef|(3)}} + |{{EquationRef|(31)}} |} |} - Integrating eq. (4.231) over the control volume P (shaded area in Fig. 4.17), yields + Integrating eq. (31) over the control volume P (shaded area in first figure), yields {| class="wikitable" border="0" {| class="wikitable" border="0" Line 347: Line 347: $\begin{matrix}{}\\\end{matrix}J_{e}=J_{w}$ $\begin{matrix}{}\\\end{matrix}J_{e}=J_{w}$ - |{{EquationRef|(4)}} + |{{EquationRef|(32)}} |} |} - Instead of assuming piecewise linear distribution of $\varphi$ as with central difference scheme or assuming $\varphi$ at the face of the control volume is equal to the value of $\varphi$ at the grid point on the upwind side in the upwind scheme, the distribution of $\varphi$ between grid points can be taken as that obtained from the exact solution, eq.  (4.210). Applying eq. (4.210) between grid points E and P, we have + Instead of assuming piecewise linear distribution of $\varphi$ as with central difference scheme or assuming $\varphi$ at the face of the control volume is equal to the value of $\varphi$ at the grid point on the upwind side in the upwind scheme, the distribution of $\varphi$ between grid points can be taken as that obtained from the exact solution, eq.  (10). Applying eq. (10) between grid points E and P, we have {| class="wikitable" border="0" {| class="wikitable" border="0" Line 356: Line 356: $\frac{\varphi (x)-\varphi _{P}}{\varphi _{E}-\varphi _{P}}=\frac{\exp [\text{Pe}_{\Delta \text{e}}(x-x_{P})/(\delta x)_{e}]-1}{\exp (\text{Pe}_{\Delta \text{e}})-1}$ $\frac{\varphi (x)-\varphi _{P}}{\varphi _{E}-\varphi _{P}}=\frac{\exp [\text{Pe}_{\Delta \text{e}}(x-x_{P})/(\delta x)_{e}]-1}{\exp (\text{Pe}_{\Delta \text{e}})-1}$ - |{{EquationRef|(5)}} + |{{EquationRef|(33)}} |} |} - Substituting eq. (4.233) into eq. (4.230) and evaluating the result at $x=x_{e}$, the total flux of $\varphi$ at the face of control volume becomes + Substituting eq. (33) into eq. (30) and evaluating the result at $\begin{matrix} + {} & x={{x}_{e}} \\ + \end{matrix}$, the total flux of $\varphi$ at the face of control volume becomes {| class="wikitable" border="0" {| class="wikitable" border="0" Line 365: Line 367: $J_{e}=F_{e}\left[ \varphi _{P}+\frac{\varphi _{P}-\varphi _{E}}{\exp (\text{Pe}_{\Delta e})-1} \right]$ $J_{e}=F_{e}\left[ \varphi _{P}+\frac{\varphi _{P}-\varphi _{E}}{\exp (\text{Pe}_{\Delta e})-1} \right]$ - |{{EquationRef|(6)}} + |{{EquationRef|(34)}} |} |} Similarly, the total flux at the west face of the control volume is Similarly, the total flux at the west face of the control volume is Line 374: Line 376: $J_{w}=F_{w}\left[ \varphi _{W}+\frac{\varphi _{W}-\varphi _{P}}{\exp (\text{Pe}_{\Delta w})-1} \right]$ $J_{w}=F_{w}\left[ \varphi _{W}+\frac{\varphi _{W}-\varphi _{P}}{\exp (\text{Pe}_{\Delta w})-1} \right]$ - |{{EquationRef|(7)}} + |{{EquationRef|(35)}} |} |} - Substituting eqs. (4.234) and (4.235) into eq. (4.232) and rearranging the resulting equation yields + Substituting eqs. (34) and (35) into eq. (32) and rearranging the resulting equation yields {| class="wikitable" border="0" {| class="wikitable" border="0" Line 384: Line 386: - |{{EquationRef|(8)}} + |{{EquationRef|(36)}} |} |} where where Line 393: Line 395: $a_{E}=\frac{F_{e}}{\exp (\text{Pe}_{\Delta e})-1}$ $a_{E}=\frac{F_{e}}{\exp (\text{Pe}_{\Delta e})-1}$ - |{{EquationRef|(9)}} + |{{EquationRef|(37)}} |} |} Line 401: Line 403: $a_{W}=\frac{F_{w}\exp (\text{Pe}_{\Delta w})}{\exp (\text{Pe}_{\Delta w})-1}$ $a_{W}=\frac{F_{w}\exp (\text{Pe}_{\Delta w})}{\exp (\text{Pe}_{\Delta w})-1}$ - |{{EquationRef|(10)}} + |{{EquationRef|(38)}} |} |} Line 409: Line 411: $\begin{matrix}{}\\\end{matrix}a_{P}=a_{W}+a_{E}+(F_{e}-F_{w})$ $\begin{matrix}{}\\\end{matrix}a_{P}=a_{W}+a_{E}+(F_{e}-F_{w})$ - |{{EquationRef|(11)}} + |{{EquationRef|(39)}} |} |} - Equations (4.237) and (4.238) can be rewritten in a format similar to that of eqs. (4.225) – (4.228), i.e., + Equations (37) and (38) can be rewritten in a format similar to that of eqs. (25) – (28), i.e., {| class="wikitable" border="0" {| class="wikitable" border="0" Line 418: Line 420: $a_{E}/D_{e}=\frac{\text{Pe}_{\Delta e}}{\exp (\text{Pe}_{\Delta e})-1}$ $a_{E}/D_{e}=\frac{\text{Pe}_{\Delta e}}{\exp (\text{Pe}_{\Delta e})-1}$ - |{{EquationRef|(12)}} + |{{EquationRef|(40)}} |} |} Line 426: Line 428: $a_{W}/D_{w}=\frac{\text{Pe}_{\Delta w}\exp (\text{Pe}_{\Delta w})}{\exp (\text{Pe}_{\Delta w})-1}$ $a_{W}/D_{w}=\frac{\text{Pe}_{\Delta w}\exp (\text{Pe}_{\Delta w})}{\exp (\text{Pe}_{\Delta w})-1}$ - |{{EquationRef|(13)}} + |{{EquationRef|(41)}} |} |} - The comparison of $a_{E}/D_{e}$ for different schemes is shown in Fig. 1. It can be seen that the hybrid scheme can be viewed as an envelope of the exponential scheme. The hybrid scheme is a good approximation if the absolute value of the grid Peclet number is either very large or near zero. + [[Image:Fig4.18.png|thumb|400 px|alt=Comparison among different schemes |Comparison among different schemes.]] - [[Image:Fig4.18.png|thumb|400 px|alt=Comparison among different schemes |Figure 1: Comparison among different schemes.]] + The comparison of $\begin{matrix}{} & {{a}_{E}}/{{D}_{e}} \\\end{matrix}$ for different schemes is shown in the figure to the right. It can be seen that the hybrid scheme can be viewed as an envelope of the exponential scheme. The hybrid scheme is a good approximation if the absolute value of the grid Peclet number is either very large or near zero. - While the exponential scheme is accurate, the computational time is much longer than for the central difference, upwind or hybrid schemes. Patankar (1981) proposed a power law scheme that has almost the same accuracy as the exponential scheme but a substantially shorter computational time. The coefficient of the neighbor grid point on the east side can be obtained by + While the exponential scheme is accurate, the computational time is much longer than for the central difference, upwind or hybrid schemes. Patankar Patankar, S.V., 1981, “A Calculation Procedure for Two-Dimensional Elliptic Situations,” Numerical Heat Transfer, Vol. 4, pp. 409-425. proposed a power law scheme that has almost the same accuracy as the exponential scheme but a substantially shorter computational time. The coefficient of the neighbor grid point on the east side can be obtained by - + - $a_{E}/D_{e}=\left\{ \begin{matrix} + [itex]a_{E}/D_{e}=\left\{ \begin{matrix} -\text{Pe}_{\Delta e}\text{ Pe}_{\Delta e}<-10 \\ -\text{Pe}_{\Delta e}\text{ Pe}_{\Delta e}<-10 \\ (1+0.1\text{Pe}_{\Delta e})^{5}-\text{Pe}_{\Delta e}\text{ }-1\text{0}\le \text{Pe}_{\Delta e}<0 \\ (1+0.1\text{Pe}_{\Delta e})^{5}-\text{Pe}_{\Delta e}\text{ }-1\text{0}\le \text{Pe}_{\Delta e}<0 \\ (1-0.1\text{Pe}_{\Delta e})^{5}\text{ 0}\le \text{Pe}_{\Delta e}\le 10 \\ (1-0.1\text{Pe}_{\Delta e})^{5}\text{ 0}\le \text{Pe}_{\Delta e}\le 10 \\ 0\text{ Pe}_{\Delta e}>10 \\ 0\text{ Pe}_{\Delta e}>10 \\ - \end{matrix} \right.$ + \end{matrix} \right.
which can be rewritten in the following compact form which can be rewritten in the following compact form Line 449: Line 450: $a_{E}/D_{e}=\left[\!\left[ 0,\left( 1-0.1\left| \text{Pe}_{\Delta e} \right| \right)^{5} \right]\!\right]+\left[\!\left[ 0,-\text{Pe}_{\Delta e} \right]\!\right]$ $a_{E}/D_{e}=\left[\!\left[ 0,\left( 1-0.1\left| \text{Pe}_{\Delta e} \right| \right)^{5} \right]\!\right]+\left[\!\left[ 0,-\text{Pe}_{\Delta e} \right]\!\right]$ - |{{EquationRef|(14)}} + |{{EquationRef|(42)}} |} |} - [[A Generalized Expression of Discretization Schemes]]
+ ==A Generalized Expression of Discretization Schemes== - The above discretization schemes can be expressed in a single generalized form. The total flux J at the interface between two grid points that were defined in eq. (4.230)  can be used to define: + The above discretization schemes can be expressed in a single generalized form. The total flux ''J'' at the interface between two grid points that were defined in eq. (30)  can be used to define: {| class="wikitable" border="0" {| class="wikitable" border="0" Line 460: Line 461: $J^{*}=\frac{J}{\Gamma /\delta x}=\text{Pe}_{\Delta }\varphi -\frac{d\varphi }{d(x/\delta x)}$ $J^{*}=\frac{J}{\Gamma /\delta x}=\text{Pe}_{\Delta }\varphi -\frac{d\varphi }{d(x/\delta x)}$ - |{{EquationRef|(1)}} + |{{EquationRef|(43)}} |} |} - which relates to the values of $\varphi$ at grid points i and i+1 (see Fig. 4.19). The first term on the right side of eq. (4.243) will be related to some weighted average of  $\varphi _{i}$ and $\varphi _{i+1}$, and the second term will be related to the difference between  $\varphi _{i}$ and $\varphi _{i+1}$. Thus, one can express + - $J^{*}$ the total flux as (Patankar, 1980) + [[Image:Fig4.19.png|thumb|400 px|alt=Total flux at the interface between grid points ''i'' and ''i''+1  | Total flux at the interface between grid points ''i'' and ''i''+1 .]] + + which relates to the values of $\varphi$ at grid points ''i'' and ''i''+1 (see figure to the right). The first term on the right side of eq. (43) will be related to some weighted average of  $\varphi _{i}$ and $\varphi _{i+1}$, and the second term will be related to the difference between  $\varphi _{i}$ and $\varphi _{i+1}$. Thus, one can express + $J^{*}$ the total flux as {| class="wikitable" border="0" {| class="wikitable" border="0" Line 470: Line 474: $J^{*}=B(\text{Pe}_{\Delta })\varphi _{i}-A(\text{Pe}_{\Delta })\varphi _{i+1}$ $J^{*}=B(\text{Pe}_{\Delta })\varphi _{i}-A(\text{Pe}_{\Delta })\varphi _{i+1}$ - |{{EquationRef|(2)}} + |{{EquationRef|(44)}} |} |} - where A and B are dimensionless coefficients that are functions of the grid Peclet number. If the field of $\varphi$ is uniform, we will have $d\varphi /dx=0$ and  eq. (4.243) becomes + where ''A'' and ''B'' are dimensionless coefficients that are functions of the grid Peclet number. If the field of $\varphi$ is uniform, we will have $d\varphi /dx=0$ and  eq. (43) becomes {| class="wikitable" border="0" {| class="wikitable" border="0" Line 479: Line 483: $J^{*}=\text{Pe}_{\Delta }\varphi _{i}=\text{Pe}_{\Delta }\varphi _{i+1}$ $J^{*}=\text{Pe}_{\Delta }\varphi _{i}=\text{Pe}_{\Delta }\varphi _{i+1}$ - |{{EquationRef|(3)}} + |{{EquationRef|(45)}} |} |} - Comparing eqs. (4.244) and (4.245) yields + Comparing eqs. (44) and (45) yields {| class="wikitable" border="0" {| class="wikitable" border="0" Line 488: Line 492: $\begin{matrix}{}\\\end{matrix}B(\text{Pe}_{\Delta })-A(\text{Pe}_{\Delta })=\text{Pe}_{\Delta }$ $\begin{matrix}{}\\\end{matrix}B(\text{Pe}_{\Delta })-A(\text{Pe}_{\Delta })=\text{Pe}_{\Delta }$ - |{{EquationRef|(4)}} + |{{EquationRef|(46)}} |} |} - [[Image:Fig4.19.png|thumb|400 px|alt=Total flux at the interface between grid points i and i+1  |Figure 1: Total flux at the interface between grid points i and i+1 .]] + For the grid system shown in the above figure, if we reconsider the problem in a reversed coordinate system $x'$ ($x'=-x$), the grid Peclet number will become $-\text{Pe}_{\Delta }$ and $J^{*}$ becomes - + - For the grid system shown in Fig. 1, if we reconsider the problem in a reversed coordinate system $x'$ ($x'=-x$), the grid Peclet number will become $-\text{Pe}_{\Delta }$ and $J^{*}$ becomes + {| class="wikitable" border="0" {| class="wikitable" border="0" Line 500: Line 502: $J^{*}=B(-\text{Pe}_{\Delta })\varphi _{i+1}-A(-\text{Pe}_{\Delta })\varphi _{i}$ $J^{*}=B(-\text{Pe}_{\Delta })\varphi _{i+1}-A(-\text{Pe}_{\Delta })\varphi _{i}$ - |{{EquationRef|(5)}} + |{{EquationRef|(47)}} |} |} - The symmetric properties of A and B can be obtained by comparing eqs. (4.244) and (4.247), i.e., + The symmetric properties of ''A'' and ''B'' can be obtained by comparing eqs. (44) and (47), i.e., {| class="wikitable" border="0" {| class="wikitable" border="0" Line 509: Line 511: $\begin{matrix}{}\\\end{matrix}A(-\text{Pe}_{\Delta })=B(\text{Pe}_{\Delta })$ $\begin{matrix}{}\\\end{matrix}A(-\text{Pe}_{\Delta })=B(\text{Pe}_{\Delta })$ - |{{EquationRef|(6)}} + |{{EquationRef|(48)}} |} |} Line 517: Line 519: $\begin{matrix}{}\\\end{matrix}B(-\text{Pe}_{\Delta })=A(\text{Pe}_{\Delta })$ $\begin{matrix}{}\\\end{matrix}B(-\text{Pe}_{\Delta })=A(\text{Pe}_{\Delta })$ - |{{EquationRef|(7)}} + |{{EquationRef|(49)}} |} |} - For the exponential schemes discussed above, one can obtain $J^{*}$ from eq(4.234)or (4.235), i.e., + For the exponential schemes discussed above, one can obtain $\begin{matrix}{} & {{J}^{*}} \\\end{matrix}$ from eq(34)or (35), i.e., - \begin{align} + [itex]\begin{align} & J^{*}=\text{Pe}_{\Delta }\left[ \varphi _{i}+\frac{\varphi _{i}-\varphi _{i+1}}{\exp (\text{Pe}_{\Delta })-1} \right] \\ & J^{*}=\text{Pe}_{\Delta }\left[ \varphi _{i}+\frac{\varphi _{i}-\varphi _{i+1}}{\exp (\text{Pe}_{\Delta })-1} \right] \\ & =\frac{\exp (\text{Pe}_{\Delta })\text{Pe}_{\Delta }}{\exp (\text{Pe}_{\Delta })-1}\varphi _{i}-\frac{\text{Pe}_{\Delta }}{\exp (\text{Pe}_{\Delta })-1}\varphi _{i+1} \\ & =\frac{\exp (\text{Pe}_{\Delta })\text{Pe}_{\Delta }}{\exp (\text{Pe}_{\Delta })-1}\varphi _{i}-\frac{\text{Pe}_{\Delta }}{\exp (\text{Pe}_{\Delta })-1}\varphi _{i+1} \\ - \end{align} + \end{align}
- Comparing the above expression with eq. (4.244), one obtains + Comparing the above expression with eq. (44), one obtains - $B=\frac{\exp (\text{Pe}_{\Delta })\text{Pe}_{\Delta }}{\exp (\text{Pe}_{\Delta })-1},\text{ }A=\frac{\text{Pe}_{\Delta }}{\exp (\text{Pe}_{\Delta })-1}$ +
$B=\frac{\exp (\text{Pe}_{\Delta })\text{Pe}_{\Delta }}{\exp (\text{Pe}_{\Delta })-1},\text{ }A=\frac{\text{Pe}_{\Delta }}{\exp (\text{Pe}_{\Delta })-1} - It can be verified that the above A and B satisfy eqs. (4.246), and (4.248) – (4.249). + It can be verified that the above ''A'' and ''B'' satisfy eqs. (46), and (48) – (49). - The implication of the above properties of A and B is that if the function A(PeΔ) for the case that [itex]\text{Pe}_{\Delta }>0$ is known, the expressions of A and B for all $\text{Pe}_{\Delta }$ can be obtained. For example, if $\text{Pe}_{\Delta }<0$, eq. (4.246) can be used to obtain + The implication of the above properties of ''A'' and ''B'' is that if the function ''A''(PeΔ) for the case that $\text{Pe}_{\Delta }>0$ is known, the expressions of ''A'' and ''B'' for all $\text{Pe}_{\Delta }$ can be obtained. For example, if $\text{Pe}_{\Delta }<0$, eq. (46) can be used to obtain - $A(\text{Pe}_{\Delta })=B(\text{Pe}_{\Delta })-\text{Pe}_{\Delta }$ +
$\begin{matrix}{} & {} \\\end{matrix}A(\text{Pe}_{\Delta })=B(\text{Pe}_{\Delta })-\text{Pe}_{\Delta } - Substituting eq. (4.248) into the above equation yields + Substituting eq. (48) into the above equation yields - [itex]A(\text{Pe}_{\Delta })=A(-\text{Pe}_{\Delta })-\text{Pe}_{\Delta }$ +
$\begin{matrix}{} & {} \\\end{matrix}A(\text{Pe}_{\Delta })=A(-\text{Pe}_{\Delta })-\text{Pe}_{\Delta } - Considering [itex]-\text{Pe}_{\Delta }=\left| \text{Pe}_{\Delta } \right|$ for the case that $\text{Pe}_{\Delta }<0$, the above expression can be rewritten as + Considering $-\text{Pe}_{\Delta }=\left| \text{Pe}_{\Delta } \right|$ for the case that $\begin{matrix}{} & {} \\\end{matrix}\text{Pe}_{\Delta }<0$, the above expression can be rewritten as - $A(\text{Pe}_{\Delta })=A(\left| \text{Pe}_{\Delta } \right|)+\left| \text{Pe}_{\Delta } \right|\text{ for Pe}_{\Delta }<0$ +
$A(\text{Pe}_{\Delta })=A(\left| \text{Pe}_{\Delta } \right|)+\left| \text{Pe}_{\Delta } \right|\text{ for Pe}_{\Delta }<0 - Since [itex]A(\text{Pe}_{\Delta })=A(\left| \text{Pe}_{\Delta } \right|)\text{ for Pe}_{\Delta }>0$ the following expression for A under any grid Peclet number can be expressed as + Since $A(\text{Pe}_{\Delta })=A(\left| \text{Pe}_{\Delta } \right|)\text{ for Pe}_{\Delta }>0$ the following expression for ''A'' under any grid Peclet number can be expressed as {| class="wikitable" border="0" {| class="wikitable" border="0" Line 550: Line 552: $A(\text{Pe}_{\Delta })=A(\left| \text{Pe}_{\Delta } \right|)+\left[\!\left[ -\text{Pe}_{\Delta },0 \right]\!\right]$ $A(\text{Pe}_{\Delta })=A(\left| \text{Pe}_{\Delta } \right|)+\left[\!\left[ -\text{Pe}_{\Delta },0 \right]\!\right]$ - |{{EquationRef|(8)}} + |{{EquationRef|(50)}} |} |} - Similarly, the expression of B for any grid Peclet number can be expressed as (see Problem 4.24). + Similarly, the expression of ''B'' for any grid Peclet number can be expressed as. {| class="wikitable" border="0" {| class="wikitable" border="0" Line 559: Line 561: $B(\text{Pe}_{\Delta })=A(\left| \text{Pe}_{\Delta } \right|)+\left[\!\left[ \text{Pe}_{\Delta },0 \right]\!\right]$ $B(\text{Pe}_{\Delta })=A(\left| \text{Pe}_{\Delta } \right|)+\left[\!\left[ \text{Pe}_{\Delta },0 \right]\!\right]$ - |{{EquationRef|(9)}} + |{{EquationRef|(51)}} |} |} - Therefore, different discretization schemes for the convection-diffusion terms can be characterized by different + Therefore, different discretization schemes for the convection-diffusion terms can be characterized by different $A(\left|text{Pe}_{\Delta } \right|)$. - $A(\left| \text{Pe}_{\Delta } \right|)$ + - . + To derive the generalized formula for different discretization schemes, let us begin from eq. (32), i.e., - To derive the generalized formula for different discretization schemes, let us begin from eq. (4.232), i.e., + {| class="wikitable" border="0" {| class="wikitable" border="0" Line 571: Line 572: $J_{e}^{*}\text{D}_{e}=J_{w}^{*}\text{D}_{w}$ $J_{e}^{*}\text{D}_{e}=J_{w}^{*}\text{D}_{w}$ - |{{EquationRef|(10)}} + |{{EquationRef|(52)}} |} |} - The total fluxes at the faces of the control volumes can be obtained from eq. (4.244), i.e., + The total fluxes at the faces of the control volumes can be obtained from eq. (44), i.e., {| class="wikitable" border="0" {| class="wikitable" border="0" Line 580: Line 581: $J_{e}^{*}=B(\text{Pe}_{\Delta e})\varphi _{P}-A(\text{Pe}_{\Delta e})\varphi _{E}$ $J_{e}^{*}=B(\text{Pe}_{\Delta e})\varphi _{P}-A(\text{Pe}_{\Delta e})\varphi _{E}$ - |{{EquationRef|(11)}} + |{{EquationRef|(53)}} |} |} Line 588: Line 589: $J_{w}^{*}=B(\text{Pe}_{\Delta w})\varphi _{W}-A(\text{Pe}_{\Delta w})\varphi _{P}$ $J_{w}^{*}=B(\text{Pe}_{\Delta w})\varphi _{W}-A(\text{Pe}_{\Delta w})\varphi _{P}$ - |{{EquationRef|(12)}} + |{{EquationRef|(54)}} |} |} - Substituting the above expressions into eq. (4.252) and rearranging the resulting equation yields + Substituting the above expressions into eq. (52) and rearranging the resulting equation yields - $\left[ D_{e}B(\text{Pe}_{\Delta e})+D_{w}A(\text{Pe}_{\Delta w}) \right]\varphi _{P}=D_{e}A(\text{Pe}_{\Delta e})\varphi _{E}+D_{w}B(\text{Pe}_{\Delta w})\varphi _{W}$ +
$\left[ D_{e}B(\text{Pe}_{\Delta e})+D_{w}A(\text{Pe}_{\Delta w}) \right]\varphi _{P}=D_{e}A(\text{Pe}_{\Delta e})\varphi _{E}+D_{w}B(\text{Pe}_{\Delta w})\varphi _{W} which can be rearranged as which can be rearranged as Line 601: Line 602: [itex]a_{P}\varphi _{P}=a_{E}\varphi _{E}+a_{W}\varphi _{W}$ $a_{P}\varphi _{P}=a_{E}\varphi _{E}+a_{W}\varphi _{W}$ - |{{EquationRef|(13)}} + |{{EquationRef|(55)}} |} |} where where Line 610: Line 611: $a_{E}=D_{e}A(\text{Pe}_{\Delta e})=D_{e}\left\{ A(\left| \text{Pe}_{\Delta e} \right|)+\left[\!\left[ -\text{Pe}_{\Delta e},0 \right]\!\right] \right\}$ $a_{E}=D_{e}A(\text{Pe}_{\Delta e})=D_{e}\left\{ A(\left| \text{Pe}_{\Delta e} \right|)+\left[\!\left[ -\text{Pe}_{\Delta e},0 \right]\!\right] \right\}$ - |{{EquationRef|(14)}} + |{{EquationRef|(56)}} |} |} Line 618: Line 619: $a_{W}=D_{w}B(\text{Pe}_{\Delta w})=D_{w}\left\{ A(\left| \text{Pe}_{\Delta w} \right|)+\left[\!\left[ \text{Pe}_{\Delta w},0 \right]\!\right] \right\}$ $a_{W}=D_{w}B(\text{Pe}_{\Delta w})=D_{w}\left\{ A(\left| \text{Pe}_{\Delta w} \right|)+\left[\!\left[ \text{Pe}_{\Delta w},0 \right]\!\right] \right\}$ - |{{EquationRef|(15)}} + |{{EquationRef|(57)}} |} |} Line 625: Line 626: |- |- | width="100%" |
| width="100%" |
- $a_{P}=a_{E}+a_{W}+(F_{e}-F_{w})$ + $\begin{matrix}{} & {} \\\end{matrix}a_{P}=a_{E}+a_{W}+(F_{e}-F_{w})$
- |{{EquationRef|(16)}} + |{{EquationRef|(58)}} |} |} - [[Image:Fig4.20.png|thumb|400 px|alt=Comparison of A(|PeΔ|) for different schemes |Figure 2: Comparison of A(|PeΔ|) for different schemes.]] + [[Image:Fig4.20.png|thumb|400 px|alt=Comparison of A(|PeΔ|) for different schemes | Comparison of $\text{Pe}_{\Delta }$ for different schemes.]] - In arriving at eqs. (4.256) and (4.257), A and B were obtained from eqs. (4.250) and (4.251). At this point, it is apparent that different discretization schemes can be characterized by different expressions for A(|PeΔ|). By comparing eqs. (4.256) and (4.257) with different expressions of aE and aW for different schemes, the corresponding A(|PeΔ|) for different schemes can be summarized in Table 1 and plotted in Fig. 2. It should be noted that the difference between the power law and exponential scheme is exaggerated for clear presentation. The generalized formula represented by eqs. (4.255) – (4.258) will be very helpful to develop a generalized computer code for all schemes. A special module or subroutine can be written for different schemes. + In arriving at eqs. (56) and (57), ''A'' and ''B'' were obtained from eqs. (50) and (51). At this point, it is apparent that different discretization schemes can be characterized by different expressions for A(|PeΔ|). By comparing eqs. (56) and (57) with different expressions of ''a''E and ''a''W for different schemes, the corresponding A(|PeΔ|) for different schemes can be summarized in the following table and plotted in the figure to the rightFaghri, A., Zhang, Y., and Howell, J. R., 2010, Advanced Heat and Mass Transfer, Global Digital Press, Columbia, MO.. It should be noted that the difference between the power law and exponential scheme is exaggerated for clear presentation. The generalized formula represented by eqs. (55) – (58) will be very helpful to develop a generalized computer code for all schemes. A special module or subroutine can be written for different schemes.
- '''Table 1''' Summary of A(|PeΔ|) for different schemes + '''Table''' Summary of A(|PeΔ|) for different schemes {| class="wikitable" border="1" {| class="wikitable" border="1" | align="center" style="background:#f0f0f0;" width="30%" | Scheme | align="center" style="background:#f0f0f0;" width="30%" | Scheme Line 657: Line 658: |} |}
+ + ==References== + {{Reflist}}

## Exact Solution

 $\frac{d(\rho u\varphi )}{dx}=\frac{d}{dx}\left( \Gamma \frac{d\varphi }{dx} \right)$ (1)

where the velocity, u, is assumed to be known. Both density, ρ, and diffusivity, Г, are assumed to be constants. The continuity equation for this one-dimensional problem is

 $\frac{d(\rho u)}{dx}=0$ (2)

Equation (1) is subject to the following boundary conditions:

 $\varphi =\varphi _{0},\text{ }x=0$ (3)
 $\varphi =\varphi _{L},\text{ }x=L$ (4)

By introducing the following dimensionless variables

 $\Phi =\frac{\varphi -\varphi _{0}}{\varphi _{L}-\varphi _{0}},\text{ }X=\frac{x}{L}$ (5)

the one-dimensional steady-state convection and diffusion problem can be nondimensionalized as

 $\text{Pe}\frac{d\Phi }{dX}=\frac{d^{2}\Phi }{dX^{2}}$ (6)
 $\begin{matrix}{}\\\end{matrix}\Phi =0,\text{ X}=0$ (7)
 $\begin{matrix}{}\\\end{matrix}\Phi =1,\text{ }X=1$ (8)

where

 $\text{Pe}=\frac{\rho uL}{\Gamma }$ (9)

is the Peclet number that reflects the relative level of convection and diffusion. Pe becomes zero for the case of pure diffusion and becomes infinite for the case of pure advection. The exact solution of eqs. (6) - (8) can be obtained as

 $\Phi =\frac{\varphi -\varphi _{0}}{\varphi _{L}-\varphi _{0}}=\frac{\exp (\text{Pe}X)-1}{\exp (\text{Pe})-1}$ (10)

which will be used as a criterion to check the accuracy of various discretization schemes.

## Central Difference Scheme

Integrating eq. (1) over the control volume P (shaded area in the figure to the right), one obtains

Control volume for one-dimensional problem.
 $(\rho u\varphi )_{e}-(\rho u\varphi )_{w}=\left( \Gamma \frac{d\varphi }{dx} \right)_{e}-\left( \Gamma \frac{d\varphi }{dx} \right)_{w}$ (11)

The right-hand side of eq. (11) can be obtained by assuming the distribution of $\varphi$ between any two neighboring grid points is piecewise linear, i.e.,

$\left( \Gamma \frac{d\varphi }{dx} \right)_{e}=\Gamma _{e}\frac{\varphi _{E}-\varphi _{P}}{(\delta x)_{e}}$

$\left( \Gamma \frac{d\varphi }{dx} \right)_{w}=\Gamma _{w}\frac{\varphi _{P}-\varphi _{W}}{(\delta x)_{w}}$

where Γe and Γw are the diffusivities at the faces of the control volume. To ensure that the flux of $\varphi$ across the faces of the control volume is continuous, the harmonic mean diffusivity at the faces should be used. To evaluate the left hand side of eq. (11), it is necessary to know the values of $\varphi$ at the faces of the control volume. If the piecewise linear profile of $\varphi$ is chosen, it follows that

$\varphi _{e}=\frac{\varphi _{E}+\varphi _{P}}{2}$

$\varphi _{w}=\frac{\varphi _{P}+\varphi _{W}}{2}$

Therefore, eq. (11) becomes

 $(\rho u)_{e}\frac{\varphi _{E}+\varphi _{P}}{2}-(\rho u)_{w}\frac{\varphi _{P}+\varphi _{W}}{2}=\Gamma _{e}\frac{\varphi _{E}-\varphi _{P}}{(\delta x)_{e}}-\Gamma _{w}\frac{\varphi _{P}-\varphi _{W}}{(\delta x)_{w}}$ (12)

Defining the mass flux and diffusive conductance

 $F=\rho u,\text{ }D=\frac{\Gamma }{\delta x}$ (13)

eq. (12) can be rearranged as

 $a_{P}\varphi _{P}=a_{E}\varphi _{E}+a_{W}\varphi _{W}$ (14)

where

 $a_{E}=D_{e}-\frac{1}{2}F_{e}$ (15)
 $a_{W}=D_{w}+\frac{1}{2}F_{w}$ (16)
 $\begin{matrix}{}\\\end{matrix}a_{P}=a_{W}+a_{E}+(F_{e}-F_{w})$ (17)

This scheme is termed the central difference scheme because the values of $\varphi$ at the faces of the control volume are taken as the averaged value between two grid points. The continuity equation requires that Fe = Fw and therefore, eq. (17) reduces to

aP = aW + aE

To evaluate the performance of the central difference scheme, let us consider the case of a uniform grid, i.e., x)e = (δx)w = δx, for which case eq. (12) can be rearranged as

 $\varphi _{P}=\frac{1}{2}\left[ \left( 1-\frac{\text{Pe}_{\Delta }}{2} \right)\varphi _{E}+\left( 1+\frac{\text{Pe}_{\Delta }}{2} \right)\varphi _{W} \right]$ (18)

where

 $\text{Pe}_{\Delta }=\frac{\rho u\delta x}{\Gamma }=\frac{F}{D}$ (19)

is the Peclet number using grid size as the characteristic length, which is referred to as the grid Peclet number. The grid Pe is a ratio of the strength of convection over diffusion. To ensure stability of the discretization scheme, the value of $\varphi _{P}$ should always fall between $\varphi _{E}$ and $\varphi _{W}$, which requires that the coefficients, $\varphi _{E}$ and $\varphi _{W}$, are positive, i.e.,

 $\left| \text{Pe}_{\Delta } \right|\le 2$ (20)

This is the criterion for stability of the central difference scheme. It can be demonstrated that the central difference becomes unstable if eq. (20) is violated. The fact that the central difference scheme is stable under small grid Peclet number indicates that the central difference scheme is accurate only if the convection is not very significant.

## Upwind Scheme

The central difference scheme assumes that the effects of the values of $\varphi$ at two neighboring grid points on the value of $\varphi$ at the face of the control volume are equal. This assumption is valid only if the effect of diffusion is dominant. If, on the other hand, the convection is dominant, one can expect that the effect of the grid point upwind is more significant than that of the point downwind. If we can assume that the value of $\varphi$ at the face of the control volume is dominated by the value of $\varphi$ at the grid point at the upwind side and that the effect of the value of $\varphi$ at the downwind side can be neglected, the two terms on the left hand side of eq. (11) can be expressed as

$(\rho u\varphi )_{e}=\left\{ \begin{matrix} F_{e}\varphi _{P},\text{ }F_{e}>0 \\ F_{e}\varphi _{E},\text{ }F_{e}<0 \\ \end{matrix} \right.$

$(\rho u\varphi )_{w}=\left\{ \begin{matrix} F_{w}\varphi _{W},\text{ }F_{w}>0 \\ F_{w}\varphi _{P},\text{ }F_{w}<0 \\ \end{matrix} \right.$

The above two equations can be expressed in the following compact form:

$(\rho u\varphi )_{e}=\varphi _{P}\left[\!\left[ F_{e},0 \right]\!\right]-\varphi _{E}\left[\!\left[ -F_{e},0 \right]\!\right]$

$(\rho u\varphi )_{w}=\varphi _{W}\left[\!\left[ F_{w},0 \right]\!\right]-\varphi _{P}\left[\!\left[ -F_{w},0 \right]\!\right]$

where the operator $\left[\!\left[ A,B \right]\!\right]$ denotes the greater of A and B [1]. Substituting the above expression into the left hand side of eq. (11) and using central difference for the right hand side of eq. (11), the discretized equation becomes

 $a_{P}\varphi _{P}=a_{E}\varphi _{E}+a_{W}\varphi _{W}$ (21)

where

 $a_{E}=D_{e}+\left[\!\left[ -F_{e},0 \right]\!\right]$ (22)
 $a_{W}=D_{w}+\left[\!\left[ F_{w},0 \right]\!\right]$ (23)
 $\begin{matrix}{}\\\end{matrix}a_{P}=a_{W}+a_{E}+(F_{e}-F_{w})$ (24)

The above scheme is referred to as the upwind scheme because the value of $\varphi$ at the grid point on the upwind side was used as the value of $\varphi$ at the face of the control volume to discretize the convection term. The upwind scheme ensures that the coefficients in eq. (21) are always positive so that a physically unrealistic solution can be avoided.

## Hybrid Scheme

The upwind scheme uses the value of $\varphi$ from the grid point at the upwind side as the value of $\varphi$ at the face of the control volume regardless of the grid Peclet number. While this treatment can yield accurate results for cases with high Peclet number, the result will not be accurate for cases where the grid Peclet number is near zero; for which cases the central difference scheme can produce better results. Spalding [2] proposed a hybrid scheme that uses the central difference scheme when $\left| \text{Pe}_{\Delta } \right|\le 2$ and the upwind scheme when $\left| \text{Pe}_{\Delta } \right|>2$.

To observe the difference between the central difference and upwind schemes, the coefficient for the east neighboring grid point, eqs. (15) and (22), can be rewritten as

 $a_{E}/D_{e}=1-\frac{1}{2}\text{Pe}_{\Delta e},\text{ Central difference scheme}$ (25)
 $a_{E}/D_{e}=1+\left[\!\left[ -\text{Pe}_{\Delta e},0 \right]\!\right],\text{ Upwind scheme}$ (26)

The hybrid scheme can then be expressed as

$a_{E}/D_{e}=\left\{ \begin{matrix} -\text{Pe}_{\Delta e}\text{ Pe}_{\Delta e}<-2 \\ 1-\frac{1}{2}\text{Pe}_{\Delta e}\text{ }-\text{2}\le \text{ Pe}_{\Delta e}\le 2 \\ 0\text{ Pe}_{\Delta e}>2 \\ \end{matrix} \right.$

which can be rewritten in the following compact form

 $a_{E}/D_{e}=\left[\!\left[ -\text{Pe}_{\Delta e},1-\frac{1}{2}\text{Pe}_{\Delta e},0 \right]\!\right]$ (27)

The coefficient for the west neighbor grid point can be obtained using a similar approach.

 $a_{W}/D_{w}=\left[\!\left[ \text{Pe}_{\Delta w},1+\frac{1}{2}\text{Pe}_{\Delta w},0 \right]\!\right]$ (28)

The above hybrid scheme combines the advantages of the central difference and upwind schemes to yield better results for cases where $\left| \text{Pe}_{\Delta } \right|\to \infty$ or $\left| \text{Pe}_{\Delta } \right|\sim 0$. However, there is still room for improvement of the solution when $\left| \text{Pe}_{\Delta } \right|$ is near 2.

## Exponential and Power Law Schemes

Since the exact solution of eq. (1) exists, one can reasonably expect that an accurate scheme can be derived if the result of the exact solution, eq. (10), is utilized. Equation (1) can be rewritten as

 $\frac{d}{dx}\left( \rho u\varphi -\Gamma \frac{d\varphi }{dx} \right)=0$ (29)

Defining the total flux of $\varphi$ due to convection and diffusion

 $J=\rho u\varphi -\Gamma \frac{d\varphi }{dx}$ (30)

eq. (29) becomes

 $\frac{dJ}{dx}=0$ (31)

Integrating eq. (31) over the control volume P (shaded area in first figure), yields

 $\begin{matrix}{}\\\end{matrix}J_{e}=J_{w}$ (32)

Instead of assuming piecewise linear distribution of $\varphi$ as with central difference scheme or assuming $\varphi$ at the face of the control volume is equal to the value of $\varphi$ at the grid point on the upwind side in the upwind scheme, the distribution of $\varphi$ between grid points can be taken as that obtained from the exact solution, eq. (10). Applying eq. (10) between grid points E and P, we have

 $\frac{\varphi (x)-\varphi _{P}}{\varphi _{E}-\varphi _{P}}=\frac{\exp [\text{Pe}_{\Delta \text{e}}(x-x_{P})/(\delta x)_{e}]-1}{\exp (\text{Pe}_{\Delta \text{e}})-1}$ (33)

Substituting eq. (33) into eq. (30) and evaluating the result at $\begin{matrix} {} & x={{x}_{e}} \\ \end{matrix}$, the total flux of $\varphi$ at the face of control volume becomes

 $J_{e}=F_{e}\left[ \varphi _{P}+\frac{\varphi _{P}-\varphi _{E}}{\exp (\text{Pe}_{\Delta e})-1} \right]$ (34)

Similarly, the total flux at the west face of the control volume is

 $J_{w}=F_{w}\left[ \varphi _{W}+\frac{\varphi _{W}-\varphi _{P}}{\exp (\text{Pe}_{\Delta w})-1} \right]$ (35)

Substituting eqs. (34) and (35) into eq. (32) and rearranging the resulting equation yields

 $a_{P}\varphi _{P}=a_{E}\varphi _{E}+a_{W}\varphi _{W}$ (36)

where

 $a_{E}=\frac{F_{e}}{\exp (\text{Pe}_{\Delta e})-1}$ (37)
 $a_{W}=\frac{F_{w}\exp (\text{Pe}_{\Delta w})}{\exp (\text{Pe}_{\Delta w})-1}$ (38)
 $\begin{matrix}{}\\\end{matrix}a_{P}=a_{W}+a_{E}+(F_{e}-F_{w})$ (39)

Equations (37) and (38) can be rewritten in a format similar to that of eqs. (25) – (28), i.e.,

 $a_{E}/D_{e}=\frac{\text{Pe}_{\Delta e}}{\exp (\text{Pe}_{\Delta e})-1}$ (40)
 $a_{W}/D_{w}=\frac{\text{Pe}_{\Delta w}\exp (\text{Pe}_{\Delta w})}{\exp (\text{Pe}_{\Delta w})-1}$ (41)
Comparison among different schemes.

The comparison of $\begin{matrix}{} & {{a}_{E}}/{{D}_{e}} \\\end{matrix}$ for different schemes is shown in the figure to the right. It can be seen that the hybrid scheme can be viewed as an envelope of the exponential scheme. The hybrid scheme is a good approximation if the absolute value of the grid Peclet number is either very large or near zero.

While the exponential scheme is accurate, the computational time is much longer than for the central difference, upwind or hybrid schemes. Patankar [3] proposed a power law scheme that has almost the same accuracy as the exponential scheme but a substantially shorter computational time. The coefficient of the neighbor grid point on the east side can be obtained by

$a_{E}/D_{e}=\left\{ \begin{matrix} -\text{Pe}_{\Delta e}\text{ Pe}_{\Delta e}<-10 \\ (1+0.1\text{Pe}_{\Delta e})^{5}-\text{Pe}_{\Delta e}\text{ }-1\text{0}\le \text{Pe}_{\Delta e}<0 \\ (1-0.1\text{Pe}_{\Delta e})^{5}\text{ 0}\le \text{Pe}_{\Delta e}\le 10 \\ 0\text{ Pe}_{\Delta e}>10 \\ \end{matrix} \right.$

which can be rewritten in the following compact form

 $a_{E}/D_{e}=\left[\!\left[ 0,\left( 1-0.1\left| \text{Pe}_{\Delta e} \right| \right)^{5} \right]\!\right]+\left[\!\left[ 0,-\text{Pe}_{\Delta e} \right]\!\right]$ (42)

## A Generalized Expression of Discretization Schemes

The above discretization schemes can be expressed in a single generalized form. The total flux J at the interface between two grid points that were defined in eq. (30) can be used to define:

 $J^{*}=\frac{J}{\Gamma /\delta x}=\text{Pe}_{\Delta }\varphi -\frac{d\varphi }{d(x/\delta x)}$ (43)
Total flux at the interface between grid points i and i+1 .

which relates to the values of $\varphi$ at grid points i and i+1 (see figure to the right). The first term on the right side of eq. (43) will be related to some weighted average of $\varphi _{i}$ and $\varphi _{i+1}$, and the second term will be related to the difference between $\varphi _{i}$ and $\varphi _{i+1}$. Thus, one can express J * the total flux as [1]

 $J^{*}=B(\text{Pe}_{\Delta })\varphi _{i}-A(\text{Pe}_{\Delta })\varphi _{i+1}$ (44)

where A and B are dimensionless coefficients that are functions of the grid Peclet number. If the field of $\varphi$ is uniform, we will have $d\varphi /dx=0$ and eq. (43) becomes

 $J^{*}=\text{Pe}_{\Delta }\varphi _{i}=\text{Pe}_{\Delta }\varphi _{i+1}$ (45)

Comparing eqs. (44) and (45) yields

 $\begin{matrix}{}\\\end{matrix}B(\text{Pe}_{\Delta })-A(\text{Pe}_{\Delta })=\text{Pe}_{\Delta }$ (46)

For the grid system shown in the above figure, if we reconsider the problem in a reversed coordinate system x' (x' = − x), the grid Peclet number will become − PeΔ and J * becomes

 $J^{*}=B(-\text{Pe}_{\Delta })\varphi _{i+1}-A(-\text{Pe}_{\Delta })\varphi _{i}$ (47)

The symmetric properties of A and B can be obtained by comparing eqs. (44) and (47), i.e.,

 $\begin{matrix}{}\\\end{matrix}A(-\text{Pe}_{\Delta })=B(\text{Pe}_{\Delta })$ (48)
 $\begin{matrix}{}\\\end{matrix}B(-\text{Pe}_{\Delta })=A(\text{Pe}_{\Delta })$ (49)

For the exponential schemes discussed above, one can obtain $\begin{matrix}{} & {{J}^{*}} \\\end{matrix}$ from eq(34)or (35), i.e.,

\begin{align} & J^{*}=\text{Pe}_{\Delta }\left[ \varphi _{i}+\frac{\varphi _{i}-\varphi _{i+1}}{\exp (\text{Pe}_{\Delta })-1} \right] \\ & =\frac{\exp (\text{Pe}_{\Delta })\text{Pe}_{\Delta }}{\exp (\text{Pe}_{\Delta })-1}\varphi _{i}-\frac{\text{Pe}_{\Delta }}{\exp (\text{Pe}_{\Delta })-1}\varphi _{i+1} \\ \end{align}

Comparing the above expression with eq. (44), one obtains

$B=\frac{\exp (\text{Pe}_{\Delta })\text{Pe}_{\Delta }}{\exp (\text{Pe}_{\Delta })-1},\text{ }A=\frac{\text{Pe}_{\Delta }}{\exp (\text{Pe}_{\Delta })-1}$

It can be verified that the above A and B satisfy eqs. (46), and (48) – (49). The implication of the above properties of A and B is that if the function A(PeΔ) for the case that PeΔ > 0 is known, the expressions of A and B for all PeΔ can be obtained. For example, if PeΔ < 0, eq. (46) can be used to obtain

$\begin{matrix}{} & {} \\\end{matrix}A(\text{Pe}_{\Delta })=B(\text{Pe}_{\Delta })-\text{Pe}_{\Delta }$

Substituting eq. (48) into the above equation yields

$\begin{matrix}{} & {} \\\end{matrix}A(\text{Pe}_{\Delta })=A(-\text{Pe}_{\Delta })-\text{Pe}_{\Delta }$

Considering $-\text{Pe}_{\Delta }=\left| \text{Pe}_{\Delta } \right|$ for the case that $\begin{matrix}{} & {} \\\end{matrix}\text{Pe}_{\Delta }<0$, the above expression can be rewritten as

$A(\text{Pe}_{\Delta })=A(\left| \text{Pe}_{\Delta } \right|)+\left| \text{Pe}_{\Delta } \right|\text{ for Pe}_{\Delta }<0$

Since $A(\text{Pe}_{\Delta })=A(\left| \text{Pe}_{\Delta } \right|)\text{ for Pe}_{\Delta }>0$ the following expression for A under any grid Peclet number can be expressed as

 $A(\text{Pe}_{\Delta })=A(\left| \text{Pe}_{\Delta } \right|)+\left[\!\left[ -\text{Pe}_{\Delta },0 \right]\!\right]$ (50)

Similarly, the expression of B for any grid Peclet number can be expressed as.

 $B(\text{Pe}_{\Delta })=A(\left| \text{Pe}_{\Delta } \right|)+\left[\!\left[ \text{Pe}_{\Delta },0 \right]\!\right]$ (51)

Therefore, different discretization schemes for the convection-diffusion terms can be characterized by different $A(\left|text{Pe}_{\Delta } \right|)$.

To derive the generalized formula for different discretization schemes, let us begin from eq. (32), i.e.,

 $J_{e}^{*}\text{D}_{e}=J_{w}^{*}\text{D}_{w}$ (52)

The total fluxes at the faces of the control volumes can be obtained from eq. (44), i.e.,

 $J_{e}^{*}=B(\text{Pe}_{\Delta e})\varphi _{P}-A(\text{Pe}_{\Delta e})\varphi _{E}$ (53)
 $J_{w}^{*}=B(\text{Pe}_{\Delta w})\varphi _{W}-A(\text{Pe}_{\Delta w})\varphi _{P}$ (54)

Substituting the above expressions into eq. (52) and rearranging the resulting equation yields

$\left[ D_{e}B(\text{Pe}_{\Delta e})+D_{w}A(\text{Pe}_{\Delta w}) \right]\varphi _{P}=D_{e}A(\text{Pe}_{\Delta e})\varphi _{E}+D_{w}B(\text{Pe}_{\Delta w})\varphi _{W}$

which can be rearranged as

 $a_{P}\varphi _{P}=a_{E}\varphi _{E}+a_{W}\varphi _{W}$ (55)

where

 $a_{E}=D_{e}A(\text{Pe}_{\Delta e})=D_{e}\left\{ A(\left| \text{Pe}_{\Delta e} \right|)+\left[\!\left[ -\text{Pe}_{\Delta e},0 \right]\!\right] \right\}$ (56)
 $a_{W}=D_{w}B(\text{Pe}_{\Delta w})=D_{w}\left\{ A(\left| \text{Pe}_{\Delta w} \right|)+\left[\!\left[ \text{Pe}_{\Delta w},0 \right]\!\right] \right\}$ (57)

 $\begin{matrix}{} & {} \\\end{matrix}a_{P}=a_{E}+a_{W}+(F_{e}-F_{w})$ (58)
Comparison of PeΔ for different schemes.

In arriving at eqs. (56) and (57), A and B were obtained from eqs. (50) and (51). At this point, it is apparent that different discretization schemes can be characterized by different expressions for A(|PeΔ|). By comparing eqs. (56) and (57) with different expressions of aE and aW for different schemes, the corresponding A(|PeΔ|) for different schemes can be summarized in the following table and plotted in the figure to the right[4]. It should be noted that the difference between the power law and exponential scheme is exaggerated for clear presentation. The generalized formula represented by eqs. (55) – (58) will be very helpful to develop a generalized computer code for all schemes. A special module or subroutine can be written for different schemes.

Table Summary of A(|PeΔ|) for different schemes

 Scheme A(|PeΔ|) Central difference $1-0.5\left| \text{Pe}_{\Delta } \right|$ Upwind 1 Hybrid $\left[\!\left[ 0,1-0.5\left| \text{Pe}_{\Delta } \right| \right]\!\right]$ Exponential $\left| \text{Pe}_{\Delta } \right|/[\exp (\left| \text{Pe}_{\Delta } \right|)-1]$ Power Law $\left[\!\left[ 0,(1-0.1\left| \text{Pe}_{\Delta } \right|)^{5} \right]\!\right]$

## References

1. 1.0 1.1 Patankar, S.V., 1980, Numerical Heat Transfer and Fluid Flow, Hemisphere, Washington, DC.
2. Spalding, D.B., 1972, “A Novel Finite-difference Formulation for Differential Expressions Involving both First and Second Derivatives,” Int. J. Num. Methods Eng., Vol. 4, pp. 551-559.
3. Patankar, S.V., 1981, “A Calculation Procedure for Two-Dimensional Elliptic Situations,” Numerical Heat Transfer, Vol. 4, pp. 409-425.
4. Faghri, A., Zhang, Y., and Howell, J. R., 2010, Advanced Heat and Mass Transfer, Global Digital Press, Columbia, MO.