Convective Term in the Douglas-Rachford ADI Scheme

1. Introduction

For the mathematical simulation of thermal field distribution in the ground, it is necessary to account for convective heat transfer (heat transfer by means of mass transfer). Convective heat transfer is caused by the filtration of water into the ground resulting, for example, from precipitation. The temperature distribution is described by the partial differential equation of heat conduction where, in the case of convection, there is a so-called convective term. Since the computational domain is arbitrary, the heat equation is solved numerically by using finite-difference methods (FDM). The Douglas – Rachford ADI scheme is one of these methods. In this paper, we focus on modification of the scheme to account for convection in the computational domain.

2. Unconditionally stable scheme for one-dimensional heat equation with convective term

Let’s consider a one-dimensional heat equation with convective term (1):

(1) $\begin{equation*} u'_{t}=(ku'_{x})'_{x}+vu'_{x}+f(x,t),\;\;\; x\in{[0,a]},\; t\in{[0,T]} \end{equation*}$

where $k=k(x,t)$ is the conduction coefficient, $v=v(x,t)$ is the fluid velocity, $f(x,t)$ is the heat source or sink density function and $u=u(x,t)$ is the sought temperature value at the point $x$ at time $t$ . A boundary condition of the first type is specified for equation (1)

(2) $\begin{equation*} u(0,t)=\mu_{0}(t),\;\;\;u(a,t)=\mu_{1}(t), \end{equation*}$

as well as the initial temperature distribution:

(3) $\begin{equation*} u(x,0)=u_{0}(x). \end{equation*}$

For the numerical solution of problem (1) – (3), we employ the FDM. We denote the spatial mesh with uniform step by $\omega_{h}=\{x_{i}=ih|i=\overline{0,N},\; h=a/N\}$ , and the time mesh by $\omega_{\tau}=\{t_{j}=j\tau|j=\overline{0,M},\; \tau=T/M\}$ . If we substitute the convective term derivative in equation (1) with central difference, then we obtain the finite-difference equation:

(4) $\begin{equation*} \frac{y_{i}^{j+1}-y_{i}^{j}}{\tau}=\frac{k_{i,j+1}^{+0.5}\Lambda^{+}y_{i}^{j+1}-k_{i,j+1}^{-0.5}\Lambda^{-}y_{i}^{j+1}}{h}+v_{i,j+1}\Lambda^{0}y_{i}^{j+1}+f_{i,j}, \end{equation*}$

$y_{0}^{j}=\mu_{0}(t_{j}), \; y_{N}^{j}=\mu_{1}(t_{j}), \; y_{i}^{0}=u_{0}(x_{i}), \; i=\overline{0,N}, j=\overline{0,M},$

where

$k_{i,j+1}^{\pm0.5}=(k(x_{i},t_{j+1})+k(x_{i\pm1},t_{j+1}))/2, \; v_{i,j+1}=v(x_{i},t_{j+1}), \; f_{i,j}=f(x_{i},t_{j})$

and $\Lambda^{-}y_{i}^{j+1}=\frac{y_{i}^{j+1}-y_{i-1}^{j+1}}{h}, \; \Lambda^{+}y_{i}^{j+1}=\frac{y_{i+1}^{j+1}-y_{i}^{j+1}}{h}}, \; \Lambda^{0}y_{i}^{j+1}=\frac{y_{i+1}^{j+1}-y_{i-1}^{j+1}}{2h}}$

approximating the differential equation (1) with accuracy $O(\tau+h^{2})$ [1]. Unfortunately, the above difference scheme is conditionally stable [2]. Therefore, instead of the difference scheme (4), the following unconditionally stable scheme [2] with the same order of accuracy is used:

(5) $\begin{equation*} \frac{y_{i}^{j+1}-y_{i}^{j}}{\tau}=\mu_{i,j+1}\frac{k_{i,j+1}^{+0.5}\Lambda^{+}y_{i}^{j+1}-k_{i,j+1}^{-0.5}\Lambda^{-}y_{i}^{j+1}}{h}+ \end{equation*}$

$+\tilde{v}_{i,j+1}^{+}k_{i,j +1}^{+0.5}\Lambda^{+}y_{i}^{j+1}+\tilde{v}_{i,j+1}^{-}k_{i,j +1}^{-0.5}\Lambda^{-}y_{i}^{j+1}+f_{i,j},$

where

$\mu_{i,j+1}=\frac{1}{1+0.5h|\tilde{v}_{i,j+1}|}, \; \tilde{v}_{i,j+1}=\frac{v_{i,j+1}}{k_{i,j+1}}, \; \tilde{v}_{i,j+1}^{-}=(\tilde{v}_{i,j+1}-|\tilde{v}_{i,j+1}|)/2,$

$\tilde{v}_{i,j+1}^{+}=(\tilde{v}_{i,j+1}+|\tilde{v}_{i,j+1}|)/2.$

3. Convective term in the Douglas-Rachford Scheme

Let us modify the Douglas-Rachford scheme to calculate the thermal field distribution in the ground affected by filtration. Similar to the one-dimensional equation (see section 2), the Douglas-Rachford finite difference scheme is as follows:

(6) $\begin{equation*} C_{eff}(T_{i,j,s}^{n})\frac{T_{i,j,s}^{*}-T_{i,j,s}^{n}}{t_{n+1}-t_{n}}=\frac{2(\mu_{+x}^{n}\Lambda_{x}^{+}-\mu_{-x}^{n}\Lambda_{x}^{-})}{x_{i+1}-x_{i-1}}T_{i,j,s}^{*}+ \end{equation*}$

$+\frac{2(\Lambda_{y}^{+}-\Lambda_{y}^{-})T_{i,j,s}^{n}}{y_{j+1}-y_{j-1}}+\frac{2(\Lambda_{z}^{+}-\Lambda_{z}^{-})T_{i,j,s}^{n}}{z_{s+1}-z_{s-1}}+\Omega_{x}^{n}T_{i,j,s}^{*}$

(7) $\begin{equation*} C_{eff}(T_{i,j,s}^{n})\frac{T_{i,j,s}^{**}-T_{i,j,s}^{*}}{t_{n+1}-t_{n}}=\frac{2(\mu_{+y}^{n}\Lambda_{y}^{+}-\mu_{-y}^{n}\Lambda_{y}^{-})}{y_{j+1}-y_{j-1}}T_{i,j,s}^{**}- \end{equation*}$

$-\frac{2(\Lambda_{y}^{+}-\Lambda_{y}^{-})T_{i,j,s}^{n}}{y_{j+1}-y_{j-1}}+\Omega_{x}^{n}T_{i,j,s}^{*}$

(8) $\begin{equation*} C_{eff}(T_{i,j,s}^{n})\frac{T_{i,j,s}^{n+1}-T_{i,j,s}^{**}}{t_{n+1}-t_{n}}=\frac{2(\mu_{+z}^{n}\Lambda_{z}^{+}-\mu_{-y}^{n}\Lambda_{z}^{-})}{z_{s+1}-z_{s-1}}T_{i,j,s}^{n+1}- \end{equation*}$

$-\frac{2(\Lambda_{z}^{+}-\Lambda_{z}^{-})T_{i,j,s}^{n}}{z_{s+1}-z_{s-1}}+\Omega_{z}^{n}T_{i,j,s}^{*}$

where

(9) $\begin{equation*} \Lambda_{x}^{+}T_{i,j,s}=\frac{k(T_{i,j,s}^{n})+k(T_{i+1,j,s}^{n})}{2(x_{i+1}-x_{i})}(T_{i+1,j,s}-T_{i,j,s}), \end{equation*}$

$\Lambda_{x}^{-}T_{i,j,s}=\frac{k(T_{i-1,j,s}^{n})+k(T_{i+1,j,s}^{n})}{2(x_{i}-x_{i-1})}(T_{i,j,s}-T_{i-1,j,s});$

(10) $\begin{equation*} \Lambda_{y}^{+}T_{i,j,s}=\frac{k(T_{i,j,s}^{n})+k(T_{i,j+1,s}^{n})}{2(y_{j+1}-y_{j})}(T_{i,j+1,s}-T_{i,j,s}), \end{equation*}$

$\Lambda_{y}^{-}T_{i,j,s}=\frac{k(T_{i,j-1,s}^{n})+k(T_{i,j+1,s}^{n})}{2(y_{j}-y_{j-1})}(T_{i,j,s}-T_{i,j-1,s});$

(11) $\begin{equation*} \Lambda_{z}^{+}T_{i,j,s}=\frac{k(T_{i,j,s}^{n})+k(T_{i,j,s+1}^{n})}{2(z_{s+1}-z_{s})}(T_{i,j,s+1}-T_{i,j,s}), \end{equation*}$

$\Lambda_{z}^{-}T_{i,j,s}=\frac{k(T_{i,j,s-1}^{n})+k(T_{i,j,s+1}^{n})}{2(z_{s}-z_{s-1})}(T_{i,j,s}-T_{i,j,s-1});$

and

(12) $\begin{equation*} \mu_{+\alpha}^{n}=\frac{1}{1+0.5(\alpha_{i+1}-\alpha_{i})|\tilde{v}_{i,j,s}^{\alpha}(n)|}, \end{equation*}$

$\mu_{-\alpha}^{n}=\frac{1}{1+0.5(\alpha_{i}-\alpha_{i-1})|\tilde{v}_{i,j,s}^{\alpha}(n)|}, \; \tilde{v}_{i,j,s}^{\alpha}(n)=\frac{v_{i,j,s}^{\alpha}(n)}{k(T_{i,j,s}^{n})}$

(13) $\begin{equation*} \Omega_{\alpha}^{n}T_{i,j,s}=C_{w}(\tilde{v}_{i,j,s}^{+\alpha}(n)\Lambda_{\alpha}^{+}+\tilde{v}_{i,j,s}^{-\alpha}(n)\Lambda_{\alpha}^{-})T_{i,j,s} \end{equation*}$

$\tilde{v}_{i,j,s}^{\pm\alpha}(n)=\frac{(\tilde{v}_{i,j,s}^{\alpha}(n)\pm|\tilde{v}_{i,j,s}^{\alpha}(n)|)}{2}, \; \alpha\in\{x,y,z\}$

In formulas (6) – (12) the $C_{eff}(T_{i,j,s}^{n})$ is the effective heat capacity, $k(T_{i,j,s}^{n})$ is the effective conductivity of ground, $\vec{v}_{i,j,s}(n)=(v_{i,j,s}^{x}(n),v_{i,j,s}^{y}(n),v_{i,j,s}^{z}(n))^{T}$ is the filtration velocity vector, $T_{i,j,s}^{n}$ is the approximate temperature value at the node $(i,j,s)$ of the spatial mesh $\omega_{h}=\{(x_{i},y_{j},z_{s})|i=\overline{0,N_{x}}, j=\overline{0,N_{y}}, s=\overline{0,N_{z}}\}$ at the time $t_{n}\in{\omega_{\tau}}=\{t_{n}|n=\overline{0, N_{t}-1}\}$ . $C_{w}$ in formula (13) defines the volumetric heat capacity of the water.

4. Benckmarks

When using the approach described in section 3 for convective heat transfer, there is no need to change the stability criterion of a difference scheme.
Here is an example of a problem where we compare the computational time using the ADI method with a high filtration velocity setting against the computational time without filtration.

Consider the following computational domain (see Fig. 1): 3 production boreholes and 9 different grounds (layers). The geometric dimensions of the computational domain are: length – 60 m, width – 20 m, height – 100 m. A sufficiently high filtration velocity is manually specified for the ground layer A (see Fig. 1) $\vec{v}(x,y,z;T)=(0;10^{-5}\times(1-w(x,y,z;T));0)^{T}$ m\s, where $w(x,y,z;T)$ is the ice content at point $(x,y,z)$ and temperature $T=T(x,y,z)$ .

The computation for 100 days (computational domain contained more than 1.3 million nodes) was performed on a single-threaded solver employing the ADI method and took 27 minutes, while the computation without filtration took 24 minutes.

Fig. 1: 3D model of computational domain

The results with and without filtration are shown in Figures 2 – 5.


Fig. 2: Temperature isolines around center borehole after 100 days (plane section at $y=70$ )		Fig. 3: Temperature isolines around center borehole with filtration after 100 days (plane section at $y=70$ )

Note that the solver based on the explicit finite-difference scheme (central difference for the approximation of the convective term was used), is more than 10 times slower in solving the problem with filtration vector $\vec{v}(x,y,z;T)$ specified in the computational domain. On the other hand, the time step when running such a solver is much smaller than when using the above solver, so the accuracy of the results is higher.


Fig. 4: Temperature distribution after 100 days (plane section at $x=0$ )		Fig. 5: Temperature distribution with filtration after 100 days (plane section at $x=0$ )

5. References

1. Samarskii A.A., The theory of difference schemes / A.A. Samarskii.– Moscow: Nauka, 1971. –552 p.
2. Polevikov V.K., Numerical methods of mathematical physics: the course of lectures / V.K. Polevikov. – Minsk: BSU, 2003. – 104 p.

Full Name*		E-mail*
Company*		Occupation*

Leave a Reply Cancel reply