Конвективное слагаемое в схеме Дугласа – Рекфорда

1. Введение

При моделировании тепловых полей в грунте необходимо учитывать конвективный теплообмен, т.е. перенос тепла посредством переноса вещества. В грунте конвективный теплообмен происходит из-за наличия фильтрационных процессов, причиной которых, например, являются метеорологические осадки. Распределение температуры описывается дифференциальным уравнением теплопроводности, в котором, в случае наличия конвекции, присутствует, так называемое, конвективное слагаемое. Ввиду произвольности формы расчетной области уравнение теплопроводности решают численно, например, конечно разностными методами. Одним из таких методов является схема переменных направлений Дугласа – Рекфорда. Модификации данной схемы при наличии конвекции в расчетной области посвящена настоящая заметка.

2. Монотонная схема для одномерного уравнения
теплопроводности с конвективным слагаемым

Рассмотрим одномерное уравнение теплопроводности с конвективным слагаемым

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

где k=k(x,t)  – коэффициент теплопроводности, v=v(x,t)  – скорость жидкости, f(x,t) – функция тепловых источников  и u=u(x,t) – искомая температура в точке с координатой x в момент времени t. Пусть для уравнения (1) задано граничное условие первого рода

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

а также начальное распределение температуры

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

Для численного решения задачи (1) – (3) воспользуемся методом конечных разностей. Через \omega_{h}=\{x_{i}=ih|i=\overline{0,N},\; h=a/N\} обозначим пространственную сетку с равномерным шагом, а через \omega_{\tau}=\{t_{j}=j\tau|j=\overline{0,M},\; \tau=T/M\} – временную сетку. Заменив в уравнении (1) дифференциальные операторы разностными, в частности, производную в конвективном члене – центральной разностью, получим конечноразностное уравнение

(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},

где

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})

и \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}}

аппроксимирующее исходное уравнение (1) с точностью O(\tau+h^{2}) [1]. Однако приведенная разностная схема не является монотонной [2]. Вместо разностной схемы (4) на практике используют следующую монотонную схему [2], обладающую тем же порядком точности,

(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},

где

\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. Конвективный член в схеме Дугласа – Рекфорда

Модифицируем схему Дугласа – Рекфорда для расчета методом переменных направлений тепловых полей в грунте при наличии фильтрации. Для этого поступим аналогично тому, как это сделано в случае одномерного уравнения (см. пункт 2). Разностная схема Дугласа – Рекфорда примет следующий вид:

(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}^{*}

где

(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});

и

(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\}

В формулах (6) – (12) через C_{eff}(T_{i,j,s}^{n}) обозначена эффективная теплоемкость, через k(T_{i,j,s}^{n}) – эффективная теплопроводность, через \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} – вектор скоростей фильтрации, T_{i,j,s}^{n} – приближенное значение температуры в узле (i,j,s) пространственной сетки \omega_{h}=\{(x_{i},y_{j},z_{s})|i=\overline{0,N_{x}}, j=\overline{0,N_{y}}, s=\overline{0,N_{z}}\} в момент времени t_{n}\in{\omega_{\tau}}=\{t_{n}|n=\overline{0, N_{t}-1}\}. Через C_{w}  в формуле (13) обозначена объёмная теплоемкость воды.

4. Тестовые расчеты

При использовании описанного в пункте 3 подхода учета конвективного теплообмена нет необходимости изменять критерий устойчивости разностной схемы. Приведем пример задачи, в которой с одной стороны задана большая скорость фильтрации, а с другой – время расчета методом переменных направлений сравнимо с временем решения задачи без фильтрации.

Рассмотрим расчетную область (см. рис. 1), которая содержит 3 добывающие скважины и 9 различных грунтов (инженерно-геологических элементов). Геометрические размеры расчетной области следующие: длина – 60 м, ширина – 20 м, высота – 100 м. Пусть в инженерно-геологическом элементе A (см. рис. 1), задана достаточно большая априорная скорость фильтрации \vec{v}(x,y,z;T)=(0;10^{-5}\times(1-w(x,y,z;T));0)^{T} м\с, где через w(x,y,z;T) обозначено значение льдистости в точке с координатами (x,y,z) и температурой T=T(x,y,z).

Расчет на 100 дней (расчетная область содержала более 1 300 000 узлов) однопоточным решателем на основе метода переменных направлений занял 27 мин, в то время как расчет задачи без фильтрации – 24 мин.

3D модель расчетной области

Рисунок 1 – 3D модель расчетной области

Результаты расчетов задач при наличии фильтрации и без нее приведены на рисунках 2 – 5.

Температурные изолинии вокруг скважины через 100 дней Влияние фильтрации на тепловое распределение в грунте
Рисунок 2 – Температурные изолинии около центральной скважины через 100 дней (срез плоскостью y=70) Рисунок 3 – Температурные изолинии при наличии фильтрации около центральной скважины через 100 дней (срез плоскостью y=70)

Отметим, что решателю на основе явной конечно-разностной схемы, в котором для аппроксимации конвективного слагаемого использовалась центральная разность, потребовалось более чем в 10 раз больше времени для решения задачи с заданным в расчетной области вектором \vec{v}(x,y,z;T) фильтрации. С другой стороны, шаг по времени при расчете таким решателем значительно меньше, чем при расчете рассмотренным выше решателем, поэтому точность результатов выше.

Расчет задачи теплопереноса в грунте вокруг скважины Теплоперенос в породе вокруг скважины с учетом фильтрации
Рисунок 4 – Температурное распределение через 100 дней (срез плоскостью x=0) Рисунок 5 – Температурное распределение при наличии фильтрации через 100 дней (срез плоскостью x=0)

5. Литература

1. Самарский, А.А. Введение в теорию разностных схем / А.А. Самарский. – М: Наука, 1971. –552 с.
2. Полевиков, В.К. Численные методы математической физики: курс лекций / В.К. Полевиков. – Минск: БГУ, 2003. – 104 с.

Добавить комментарий

Ваш e-mail не будет опубликован. Обязательные поля помечены *


*

Можно использовать следующие HTML-теги и атрибуты: <a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <strike> <strong>