Журнал вычислительной математики и математической физики, 2021, T. 61, № 1, стр. 95-107

Монотонные схемы для задач конвекции-диффузии с конвективным переносом в различной форме

П. Н. Вабищевич^1, 2, *

¹ ИБРАЭ РАН
115191 Москва, Б. Тульская ул., 52, Россия

² СВФУ им. М.К. Аммосова
677000 Якутск, ул. Белинского, 58, Россия

^* E-mail: vabishchevich@gmail.com

Поступила в редакцию 10.03.2020
После доработки 18.06.2020
Принята к публикации 18.09.2020

DOI: 10.31857/S0044466920120157

Полный текст (PDF)

Аннотация

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

Ключевые слова: задачи конвекции-диффузии, двухслойные разностные схемы, логарифмическая норма, монотонные схемы.

ВВЕДЕНИЕ

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

При рассмотрении параболических уравнений второго порядка конвективные слагаемые чаще всего берутся в недивергентной (характеристической) или дивергентной (консервативной) формах [4]. В этом случае устойчивость наиболее естественно рассматривается в соответствующих банаховых пространствах: в равномерной норме при использовании недивергентной формы и интегральной норме для дивергентной формы. Большие возможности в вычислительной гидродинамике предоставляет [5] запись конвективных слагаемых в так называемой симметричной форме, когда берется полусумма недивергентной и дивергентной форм. В этом случае обеспечивается безусловная кососимметричность оператора конвективного переноса, его энергетическая нейтральность. Сама краевая задача конвекции-диффузии рассматривается в гильбертовом пространстве, что позволяет ориентироваться на стандартные конечно-элементные аппроксимации по пространству [6], [7] и использование общих результатов теории устойчивости (корректности) операторно-разностных схем [8], [9].

При численном решении стационарных и нестационарных задач конвекции-диффузии большое внимание уделяется монотонным аппроксимациям – выполнению принципа максимума на дискретном уровне, наследованию такого свойства для решений эллиптических и параболических уравнений второго порядка [10]. На основе принципа максимума получены [11] априорные оценки для задач конвекции-диффузии в ${{L}_{\infty }}(\Omega )$, когда конвективные слагаемые записываются в недивергентной форме. Аналогично [4] формулируется принцип максимумадля задач с конвективными слагаемыми в дивергентной форме. В этом случае соответствующие априорные оценки имеют место в пространстве ${{L}_{1}}(\Omega )$.

Безусловно монотонные схемы для задач конвекции-диффузии строятся на основе аппроксимации конвективных слагаемых направленными разностями [8]. Хорошо известен основной недостаток таких схем, который связан с первым порядком аппроксимации. В вычислительной гидродинамике [3], [12] для построения монотонных схем второго порядка аппроксимации применяются различные подходы. В основном, явно или неявно используется идея перехода к схемам с направленными разностями в области, где нарушается условие монотонности схем с обычными центрально разностными аппроксимациями конвективных слагаемых. Подобные регуляризованные разностные схемы для задач конвекции-диффузии с недивергентными и дивергентными конвективными слагаемыми рассмотрены в книге [4]. Второй класс безусловно монотонных схем строится на основе трансформации исходного уравнения конвекции-диффузии, когда конвективный и диффузионный перенос записывается в виде диффузионного переноса вспомогательной величины. В этом случае мы приходим к экспоненциальным схемам, которые предложены в работах [13], [14], и в различных вариантах широко используются в вычислительной практике.

Для задач конвекции-диффузии с конвективными слагаемыми в недивергентной и дивергентной формах естественно ориентироваться на использование банаховых пространств ${{L}_{\infty }}(\Omega )$ и ${{L}_{1}}(\Omega )$) соответственно. Исследование устойчивости в банаховых пространствах ${{L}_{\infty }}(\omega )$ и ${{L}_{1}}(\omega )$ (сеточных аналогах ${{L}_{\infty }}(\Omega )$ и ${{L}_{1}}(\Omega )$) обычно проводится на основе принципа максимума. В работе [15] условия устойчивости двухслойных разностных схем для нестационарных задач конвекции-диффузии формулируются с использованием понятия логарифмической нормы [16], [17]. Устойчивость и монотонность этих схем обеспечивается диагональным преобладанием по строкам (в ${{L}_{\infty }}(\omega )$) или по столбцам (в ${{L}_{1}}(\omega )$).

Аналогичные проблемы устойчивости и монотонности двухслойных схем применительно к задачам конвекции-диффузии, в которых конвективный перенос записан в симметричной форме, обсуждаются в настоящей работе. Рассмотрение ведется в гильбертовом пространстве ${{L}_{2}}(\omega )$, а достаточные условия устойчивости и монотонности формулируются в виде специального варианта диагонального преобладания. Для того чтобы выделить ключевые моменты без усложнения нашего исследования непринципиальными техническими деталями, в качестве основного объекта рассматривается одномерное по пространству уравнение конвекции-диффузии. Общие комментарии относительно обобщения результатов на многомерные задачи даны в конце работы.

1. УРАВНЕНИЯ КОНВЕКЦИИ-ДИФФУЗИИ

Будем рассматривать модельные нестационарные задачи конвекции-диффузии с постоянным (не зависящим от времени, но зависящим от точки расчетной области) коэффициентом диффузионного переноса. В прикладных задачах естественно ориентироваться на случай, когда коэффициенты конвективного переноса зависят также и от времени.

Сформулируем простейшие одномерные задачи конвекции-диффузии. Уравнение конвекции-диффузии с конвективными слагаемыми в недивергентном виде [4] имеет вид:

(1.1)

$\frac{{\partial u}}{{\partial t}} + v(x,t)\frac{{\partial u}}{{\partial x}} - \frac{\partial }{{\partial x}}\left( {k(x)\frac{{\partial u}}{{\partial x}}} \right) = f(x,t)$

при $0 < x < l,$ $t > 0$. Это уравнение дополним простейшими однородными граничными условиями Дирихле и начальным условием

(1.2)

$u(0,t) = 0,\quad u(l,t) = 0,\quad t > 0,$

(1.3)

$u(x,0) = {{u}^{0}}(x),\quad 0 < x < l.$

Вторым важнейшим примером является нестационарное уравнение конвекции-диффузии с конвективным переносом в дивергентной форме:

(1.4)

$\frac{{\partial u}}{{\partial t}} + \frac{\partial }{{\partial x}}(v(x,t)u) - \frac{\partial }{{\partial x}}\left( {k(x)\frac{{\partial u}}{{\partial x}}} \right) = f(x,t).$

Основным объектом нашего рассмотрения являются задачи конвекции-диффузии, в которых конвективный перенос берется в симметричной форме:

(1.5)

$\frac{{\partial u}}{{\partial t}} + \frac{1}{2}\left( {\frac{\partial }{{\partial x}}(v(x,t)u) + v(x,t)\frac{{\partial u}}{{\partial x}}} \right) - \frac{\partial }{{\partial x}}\left( {k(x)\frac{{\partial u}}{{\partial x}}} \right) = f(x,t).$

Будем рассматривать множество функций $u(x,t)$, удовлетворяющих граничным условиям (1.2). Нестационарную задачу конвекции-диффузии запишем в виде дифференциально-операторного уравнения

(1.6)

$\frac{{du}}{{dt}} + \mathcal{A}u = f(t),\quad \mathcal{A} = \mathcal{A}(t) = \mathcal{C}(t) + \mathcal{D},$

где $\mathcal{C}(t)$ – оператор конвективного переноса, а $\mathcal{D}$ – оператор диффузионного переноса. Рассматривается задача Коши для эволюционного уравнения (1.6), т.е. уравнение дополняется условием

(1.7)

$u(0) = {{u}^{0}}.$

Приведем простейшие априорные оценки для рассматриваемых задач конвекции-диффузии (1.1)–(1.3) и (1.2)–(1.4) и (1.2), (1.3), (1.5). Рассматриваются пространства ${{L}_{\infty }}(0,l)$, ${{L}_{1}}(0,l)$ и ${{L}_{2}}(0,l)$, нормы в которых есть

${{\left\| v \right\|}_{\infty }} = \mathop {max}\limits_{0 < x < l} \left| {v(x)} \right|,\quad {{\left\| v \right\|}_{1}} = \int\limits_0^l {{\text{|}}v(x){\kern 1pt} {\text{|}}dx} ,\quad {{\left\| v \right\|}_{2}} = {{\left( {\int\limits_0^l {{{v}^{2}}(x)dx} } \right)}^{{1/2}}}.$

Для решения нестационарной задачи конвекции-диффузии (1.1)–(1.3) (конвективный перенос в недивергентной форме) справедлива [18] априорная оценка в ${{L}_{\infty }}(0,l)$, которая следует из принципа максимума:

(1.8)

${{\left\| {u(x,t)} \right\|}_{\infty }} \leqslant {{\left\| {{{u}^{0}}(x)} \right\|}_{\infty }} + \int\limits_0^t {{{{\left\| {f(x,\theta )} \right\|}}_{\infty }}d\theta } .$

Для решения задачи (1.2)–(1.4) имеет (см. [4]) место аналогичная априорная оценка, но в ${{L}_{1}}(0,l)$:

(1.9)

${{\left\| {u(x,t)} \right\|}_{1}} \leqslant {{\left\| {{{u}^{0}}(x)} \right\|}_{1}} + \int\limits_0^t {{{{\left\| {f(x,\theta )} \right\|}}_{1}}d\theta } .$

Наиболее просто устанавливается оценка для задачи конвекции-диффузии с конвективными слагаемыми в симметричной форме (1.2), (1.3), (1.5) в гильбертовом пространстве ${{L}_{2}}(0,l)$:

(1.10)

${{\left\| {u(x,t)} \right\|}_{2}} \leqslant {{\left\| {{{u}^{0}}(x)} \right\|}_{2}} + \int\limits_0^t {{{{\left\| {f(x,\theta )} \right\|}}_{2}}d\theta } .$

Для доказательства (1.10) привлекаются свойства положительности самосопряженного оператора диффузионного переноса ($\mathcal{D} = \mathcal{D}{\kern 1pt} * > 0$) и кососимметричности оператора конвективного переноса ($\mathcal{C} = - \mathcal{C}{\kern 1pt} *$ ). Мы хотим иметь аналогии априорных оценок (1.8)–(1.10) при рассмотрении дискретных задач.

2. УСТОЙЧИВОСТЬ И МОНОТОННОСТЬ ДВУХСЛОЙНЫХ СХЕМ

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

Ищется решение следующей системы линейных обыкновенных уравнений первого порядка:

(2.1)

$\frac{{d{{w}_{i}}}}{{dt}} + \sum\limits_{j = 1}^m \,{{a}_{{ij}}}(t){{w}_{j}} = {{\varphi }_{i}}(t),\quad i = 1,2, \ldots ,m,\quad t > 0.$

Полагая $w = w(t) = \{ {{w}_{1}},{{w}_{2}}, \ldots ,{{w}_{m}}\} ,$ $A = \{ {{a}_{{ij}}}\} $, запишем (2.1) в матричном (операторном) виде:

(2.2)

$\frac{{dw}}{{dt}} + A(t)w = \varphi (t).$

Будем строить разностные схемы для приближенного решения задачи Коши, когда (2.2) дополняется начальными условиями

(2.3)

$w(0) = {{u}^{0}}.$

Нас интересует устойчивость разностного решения задачи (2.2), (2.3) в пространствах ${{L}_{\infty }}$, ${{L}_{1}}$ и ${{L}_{2}}$. Для нормы вектора и подчиненной ей нормы матрицы [19] в ${{L}_{\infty }}$ имеем

(2.4)

${{\left\| w \right\|}_{\infty }} = \mathop {max}\limits_{1 \leqslant i \leqslant m} {\text{|}}{{w}_{i}}{\kern 1pt} {\text{|}},\quad {{\left\| A \right\|}_{\infty }} = \mathop {max}\limits_{1 \leqslant i \leqslant m} \sum\limits_{j = 1}^m \,{\text{|}}{{a}_{{ij}}}{\kern 1pt} {\text{|}}.$

Аналогично в ${{L}_{1}}$

(2.5)

${{\left\| w \right\|}_{1}} = \sum\limits_{i = 1}^m \,{\text{|}}{{w}_{i}}{\kern 1pt} {\text{|}},\quad {{\left\| A \right\|}_{1}} = \mathop {max}\limits_{1 \leqslant j \leqslant m} \sum\limits_{i = 1}^m \,{\text{|}}{{a}_{{ij}}}{\kern 1pt} {\text{|}}.$

Для ${{L}_{2}}$ имеем

(2.6)

${{\left\| w \right\|}_{2}} = {{\left( {\sum\limits_{i = 1}^m \,{\text{|}}{{w}_{i}}{\kern 1pt} {{{\text{|}}}^{2}}} \right)}^{{1/2}}},\quad {{\left\| A \right\|}_{2}} = \mathop {max}\limits_{1 \leqslant i \leqslant m} {\kern 1pt} {\kern 1pt} {{({{\lambda }_{i}}(AA{\kern 1pt} *))}^{{1/2}}},$

где ${{\lambda }_{i}}(G)$ – собственные значения симметричной матрицы $G$.

Устойчивость решения задачи Коши (2.2), (2.3) удобно формулировать с привлечением понятия логарифмической нормы. Логарифмическая норма матрицы $A$ есть (см. [20], [21]) число

$\mu (A) = \mathop {lim}\limits_{\delta \to 0 + } \frac{{\left\| {I + \delta A} \right\| - 1}}{\delta },$

где $I$ – единичная матрица. Для решения задачи (2.2), (2.3) имеет место априорная оценка

(2.7)

$\left\| {w(t)} \right\| \leqslant exp(\mu ( - A)t)\left( {\left\| {{{u}^{0}}} \right\| + \int\limits_0^t \,exp( - \mu ( - A)(t - \theta ))\left\| {\varphi (\theta )} \right\|d\theta } \right),$

которая обеспечивает устойчивость по начальным данным и правой части.

Для логарифмической нормы матрицы в ${{L}_{\infty }}$ (согласованной с (2.4)), в ${{L}_{1}}$ (согласованной с (2.5)) и в ${{L}_{2}}$ (согласованной с (2.6)) имеют место выражения

${{\mu }_{\infty }}(A) = \mathop {max}\limits_{1 \leqslant i \leqslant m} \left( {{{a}_{{ii}}} + \sum\limits_{j \ne i,j = 1}^m \,{\text{|}}{{a}_{{ij}}}{\kern 1pt} {\text{|}}} \right),$

${{\mu }_{1}}(A) = \mathop {max}\limits_{1 \leqslant i \leqslant m} \left( {{{a}_{{ii}}} + \sum\limits_{j \ne i,j = 1}^m \,{\text{|}}{{a}_{{ji}}}{\kern 1pt} {\text{|}}} \right),$

${{\mu }_{2}}(A) = \mathop {max}\limits_{1 \leqslant i \leqslant m} \left( {{{\lambda }_{i}}\left( {\frac{{A + A{\kern 1pt} *}}{2}} \right)} \right).$

При рассмотрении задач конвекции-диффузии мы ориентируемся на априорные оценки (1.8)–(1.10). Принимая во внимание (2.7), это соответствует тому, что

(2.8)

${{\mu }_{\alpha }}( - A) \leqslant 0,\quad \alpha = \infty ,1,2.$

Сформулируем достаточные условия выполнения (2.8), которые имеют место для $\alpha = \infty ,1$.

Задачу (2.2), (2.3) будем рассматривать при следующих ограничениях. Диагональные элементы матрицы $A$ предполагаются неотрицательными и имеется диагональное преобладание в том или ином варианте. Будем считать, например, что справедливы неравенства

(2.9)

${{a}_{{ii}}} \geqslant \sum\limits_{j \ne i,j = 1}^m \,{\text{|}}{{a}_{{ij}}}{\kern 1pt} {\text{|}},\quad i = 1,2, \ldots ,m,$

т.е. имеет место нестрогое диагональное преобладание по строкам. Второй вариант связан с нестрогим диагональным преобладанием по столбцам:

(2.10)

${{a}_{{ii}}} \geqslant \sum\limits_{j \ne i,j = 1}^m \,{\text{|}}{{a}_{{ji}}}{\kern 1pt} {\text{|}},\quad i = 1,2, \ldots ,m.$

Непосредственно убеждаемся, что при выполнении (2.9) или (2.10) имеет место неравенство (2.8) при $\alpha = \infty $ или $\alpha = 1$ соответственно.

При $\alpha = 2$ логарифмическая норма определяется по симметричной части $A$, причем

${{\mu }_{2}}( - A) = \mathop {max}\limits_{1 \leqslant i \leqslant m} \left( {{{\lambda }_{i}}\left( { - \frac{{A + A{\kern 1pt} *}}{2}} \right)} \right) = - \mathop {min}\limits_{1 \leqslant i \leqslant m} \left( {{{\lambda }_{i}}\left( {\frac{{A + A{\kern 1pt} *}}{2}} \right)} \right).$

Тем самым ${{\mu }_{2}}( - A) \leqslant 0$, если $A$ – неотрицательно определенная матрица:

(2.11)

$A \geqslant 0.$

Матрица с нестрогим диагональным преобладанием является неотрицательно-определенной. Это достаточное условие для выполнения (2.11) запишем для симметричной части матрицы $A$:

(2.12)

${{a}_{{ii}}} \geqslant \frac{1}{2}\sum\limits_{j \ne i,j = 1}^m \,{\text{|}}{{a}_{{ij}}} + {{a}_{{ji}}}{\kern 1pt} {\text{|}},\quad i = 1,2, \ldots ,m.$

В силу этого при выполнении (2.12) для решения задачи (2.2), (2.3) справедлива оценка (1.10).

Для приближенного решения задачи (2.2), (2.3) будем использовать двухслойные схемы с весами. Обозначим приближенное решение на момент времени ${{t}^{n}} = n\tau $, где $\tau $ – шаг по времени, через ${{y}^{n}}$. Оно определяется из уравнения

(2.13)

$\frac{{{{y}^{{n + 1}}} - {{y}^{n}}}}{\tau } + A(\sigma {{y}^{{n + 1}}} + (1 - \sigma ){{y}^{n}}) = {{\varphi }^{n}},$

где, например, $A = {{A}^{{n + \sigma }}} = A(\sigma {{t}^{{n + 1}}} + (1 - \sigma ){{t}^{n}})$, при начальном условии

(2.14)

${{y}^{0}} = {{u}^{0}}.$

Достаточные условия устойчивости разностной схемы (2.13), (2.14) в ${{L}_{\infty }}$ и в ${{L}_{1}}$ формулируются в виде следующего утверждения.

Теорема 1. Для задачи Коши (2.2), (2.3) с матрицей $A$, удовлетворяющей условиям (2.9) (или (2.10)), разностная схема с весами (2.13), (2.14) безусловно устойчива в ${{L}_{\infty }}$ (или в ${{L}_{1}}$) при $\sigma \geqslant 1$. Если матрица $A$ удовлетворяет условиям (2.12), то схема (2.13), (2.14) безусловно устойчива в ${{L}_{2}}$ при $\sigma \geqslant 0.5$. При этом для разностного решения верна априорная оценка

(2.15)

$\left\| {{{y}^{{n + 1}}}} \right\| \leqslant \left\| {{{u}^{0}}} \right\| + \sum\limits_{k = 0}^n \,\tau \left\| {{{\varphi }^{k}}} \right\|.$

Доказательство. При рассмотрении устойчивости в ${{L}_{\infty }}$ и ${{L}_{1}}$ будем следовать [15]. Схему с весами (2.13), (2.14) удобно рассматривать как чисто неявную схему

(2.16)

$\frac{{{{y}^{{n + \sigma }}} - {{y}^{n}}}}{{\sigma \tau }} + A{{y}^{{n + \sigma }}} = {{\varphi }^{n}}$

для функции ${{y}^{{n + \sigma }}} = \sigma {{y}^{{n + 1}}} + (1 - \sigma ){{y}^{n}}$. Из (2.16) непосредственно следует

$\left\| {(I + \sigma \tau A){{y}^{{n + \sigma }}}} \right\| \leqslant \left\| {{{y}^{n}}} \right\| + \tau \sigma \left\| {{{\varphi }^{n}}} \right\|.$

Для логарифмической нормы имеют место оценки

$\left\| {Au} \right\| \geqslant - \mu (A)\left\| u \right\|,\quad \left\| {Au} \right\| \geqslant - \mu ( - A)\left\| u \right\|.$

В силу этого с учетом того, что при диагональном преобладании справедливо (2.8), имеем

$\left\| {(I + \sigma \tau A){{y}^{{n + \sigma }}}} \right\| \geqslant - \mu ( - I - \sigma \tau A)\left\| {{{y}^{{n + \sigma }}}} \right\| = (1 + \sigma \tau \mu ( - A))\left\| {{{y}^{{n + \sigma }}}} \right\| \geqslant \left\| {{{y}^{{n + \sigma }}}} \right\|.$

Тем самым приходим к оценке

(2.17)

$\left\| {{{y}^{{n + \sigma }}}} \right\| \leqslant \left\| {{{y}^{n}}} \right\| + \tau \sigma \left\| {{{\varphi }^{n}}} \right\|.$

При ограничениях $\sigma \geqslant 1$ имеем

$\left\| {{{y}^{{n + \sigma }}}} \right\| \geqslant \sigma \left\| {{{y}^{{n + 1}}}} \right\| - (\sigma - 1)\left\| {{{y}^{n}}} \right\|.$

В этих условиях из (2.17) получим неравенство

(2.18)

$\left\| {{{y}^{{n + 1}}}} \right\| \leqslant \left\| {{{y}^{n}}} \right\| + \tau \left\| {{{\varphi }^{n}}} \right\|,$

из которого следует оценка (2.15).

Для гильбертова пространства ${{L}_{2}}$ безусловная устойчивость схемы с весами (2.13), (2.14) при выполнении (2.11) имеет место [8], [9] при более слабых ограничениях на вес. Домножим скалярно уравнение (2.13) на ${{y}^{{n + \sigma }}}$, что с учетом (2.11) дает

(2.19)

$({{y}^{{n + 1}}} - {{y}^{n}},{{y}^{{n + \sigma }}}) \leqslant \tau ({{y}^{{n + \sigma }}},{{\varphi }^{n}}).$

Для оценки левой части (2.19) воспользуемся следующим вспомогательным результатом [22].

Лемма 1. Пусть $w = \sigma u + (1 - \sigma ){v}$ и постоянная $\sigma \geqslant 0.5$ для $u,\;v$ из некоторого гильбертового пространства ${{H}_{G}}$ ($G = G{\kern 1pt} * > 0$). Тогда имеем

(2.20)

$(G(u - v),w) \geqslant ({{\left\| u \right\|}_{G}} - {{\left\| v \right\|}_{G}}){{\left\| w \right\|}_{G}}.$

В нашем случае $u = {{y}^{{n + 1}}},$ $v = {{y}^{n}},$ $w = {{y}^{{n + \sigma }}},$ $G = I$ и поэтому для левой части (2.19) при $\sigma \geqslant 0.5$ неравенство (2.20) дает

$({{y}^{{n + 1}}} - {{y}^{n}},{{y}^{{n + \sigma }}}) \geqslant \left( {\left\| {{{y}^{{n + 1}}}} \right\| - \left\| {{{y}^{n}}} \right\|} \right)\left\| {{{y}^{{n + \sigma }}}} \right\|.$

С учетом неравенства

$({{y}^{{n + \sigma }}},{{\varphi }^{n}}) \leqslant \left\| {{{y}^{{n + \sigma }}}} \right\|\left\| {{{\varphi }^{n}}} \right\|$

получаем (2.19).

Замечание 1. При малых $\sigma $ мы можем рассчитывать на условную устойчивость. Например, для матриц $A$ с диагональным преобладанием (2.9) (или (2.10)) разностная схема с весами (2.13), (2.14) условно устойчива [15] при $0 \leqslant \sigma < 1$ в ${{L}_{\infty }}$ (в ${{L}_{1}}$), если

(2.21)

$\tau \leqslant \frac{1}{{1 - \sigma }}{{\left( {\mathop {max}\limits_{1 \leqslant i \leqslant m} {{a}_{{ii}}}} \right)}^{{ - 1}}}.$

C приведенными выше условиями устойчивости в виде диагонального преобладания непосредственно связывается свойство монотонности разностного решения задачи (2.13), (2.14) при дополнительном предположении о неположительности внедиагональных элементов матрицы $A$. Докажем следующее утверждение о безусловной монотонности схемы с весами (2.13), (2.14).

Теорема 2. Пусть в схеме (2.13), (2.14) выполнены условия диагонального преобладания (2.9) (или (2.10), или (2.12)) при

(2.22)

${{a}_{{ij}}} \leqslant 0,\quad i \ne j,\quad i,j = 1,2, \ldots ,m,$

и пусть

${{u}^{0}} \geqslant 0,\quad {{\varphi }^{n}} \geqslant 0,\quad n = 0,1, \ldots ,$

тогда

${{y}^{{n + 1}}} \geqslant 0,\quad n = 1,2, \ldots ,$

при любых $\tau > 0$ и $\sigma \geqslant 1$.

Доказательство. Для перехода с одного временного слоя на другой из (2.16) имеем

(2.23)

$B{{y}^{{n + \sigma }}} = {{g}^{n}},\quad n = 0,1, \ldots ,$

где

$B = I + \sigma \tau A,\quad {{g}^{n}} = {{y}^{n}} + \tau {{\varphi }^{n}}.$

Доказательство неотрицательности решения проводится по индукции. Предположим, что ${{y}^{n}} \geqslant 0$ (при $n = 0$ это имеет место в силу предположений теоремы). Покажем, что при этом неотрицательным будет и ${{y}^{{n + 1}}} \geqslant 0$.

Сначала установим, что в условиях теоремы матрица системы линейных алгебраических уравнений (2.23) является M-матрицей. Для элементов матрицы $B$ имеем

${{b}_{{ii}}} = 1 + \sigma \tau {{a}_{{ii}}},\quad {{b}_{{ij}}} = \sigma \tau {{a}_{{ij}}},\quad j = 1,2, \ldots ,m,\quad i = 1,2, \ldots ,m.$

Тем самым обеспечена положительность диагональных и неположительность внедиагональных элементов матрицы.

При (2.9) имеем строгое диагональное преобладание по строкам:

${{b}_{{ii}}} > \sum\limits_{j \ne i,j = 1}^m \,{\text{|}}{{b}_{{ij}}}{\kern 1pt} {\text{|}},\quad i = 1,2, \ldots ,m.$

Аналогично, при ограничениях (2.10) имеем строгое диагональное преобладание по столбцам:

${{b}_{{ii}}} > \sum\limits_{j \ne i,j = 1}^m \,{\text{|}}{{b}_{{ji}}}{\kern 1pt} {\text{|}},\quad i = 1,2, \ldots ,m.$

В силу этого (см. [19]) матрица $B$ является М-матрицей.

Отдельно рассматривается случай ограничений (2.12). Сформулируем следующее вспомогательное утверждение [23] (см. также [24], [25]) для матриц с $\alpha $-диагональным преобладанием.

Лемма 2. Пусть элементы матрицы $A$

${{a}_{{ii}}} > 0,\quad {{a}_{{ij}}} \leqslant 0,\quad i \ne j,\quad i,j = 1,2, \ldots ,m,$

и имеет место строгое $\alpha $-диагональное преобладание

(2.24)

${{a}_{{ii}}} > \alpha {{R}_{i}}(A) + (1 - \alpha ){{S}_{i}}(A),\quad i = 1,2, \ldots ,m,$

где

${{R}_{i}}(A) = \sum\limits_{j \ne i,j = 1}^m \,{\text{|}}{{a}_{{ij}}}{\kern 1pt} {\text{|}},\quad {{S}_{i}}(A) = \sum\limits_{j \ne i,j = 1}^m \,{\text{|}}{{a}_{{ji}}}{\kern 1pt} {\text{|}},\quad i = 1,2, \ldots ,m.$

Тогда матрица $A$ является М-матрицей при всех $\alpha \in [0,1]$.

В условиях (2.12) и (2.22) для матрицы имеем

${{b}_{{ii}}} > \frac{1}{2}({{R}_{i}}(B) + {{S}_{i}}(B)),\quad i = 1,2, \ldots ,m,$

т.е. она является матрицей с $\alpha $-диагональным преобладанием при $\alpha = 0.5$. В силу леммы 2 матрица $B$ является М-матрицей.

Для правой части (2.23) имеем ${{g}^{n}} \geqslant 0$. Поэтому для решения уравнения (2.23) с М-матрицей $B$ имеем ${{y}^{{n + \sigma }}} \geqslant 0$. Принимая во внимание, что $\sigma {{y}^{{n + 1}}} = {{y}^{{n + \sigma }}} + (\sigma - 1){{y}^{n}}$, получаем ${{y}^{{n + 1}}} \geqslant 0$ при $\sigma \geqslant 1$, что и завершает доказательство теоремы.

Замечание 2. При $0 \leqslant \sigma < 1$ монотонность обеспечивается при малых шагах по времени. Например, для матрицы $A$ с внедиагональными неположительными элементами и диагональным преобладанием по строкам или столбцам монотонность имеет место [15] при ограничениях (2.21).

Теперь мы можем использовать установленные результаты для исследования безусловной устойчивости и монотонности разностных схем для нестационарных задач конвекции-диффузии при записи конвективного переноса в различных формах.

3. РАЗНОСТНЫЕ СХЕМЫ ДЛЯ УРАВНЕНИЙ КОНВЕКЦИИ-ДИФФУЗИИ

Будем использовать на отрезке $[0,l]$ равномерные сетки, когда

$\bar {\omega } \equiv \omega \cup \partial \omega = \{ x\,{\text{|}}\,x = {{x}_{i}} = ih,\;i = 0,1, \ldots ,N,\;Nh = l\} ,$

и $\omega $ есть множество внутренних узлов:

$\omega = \{ x\,{\text{|}}\,x = {{x}_{i}} = ih,\;i = 1,2, \ldots ,N - 1,\;Nh = l\} .$

При граничных условиях Дирихле (1.2) аппроксимация уравнения конвекции-диффузии (1.1) (или (1.4), (1.6)) и начальных условий (1.3) дает задаче (2.2), (2.3), при этом $m = N - 1$ и приближенное решение ${{w}_{i}}(t) = w(x,t),$ $x \in \omega $.

Будем использовать стандартные аппроксимации на трехточечном шаблоне [8] при рассмотрении задач конвекции-диффузии [4]. Разностный оператор диффузии определим, например, в виде

(3.1)

$Dw = - \frac{1}{{{{h}^{2}}}}k(x + 0.5h)(w(x + h,t) - w(x,t)) + \frac{1}{{{{h}^{2}}}}k(x - 0.5h)(w(x,t) - w(x - h,t)),\quad x \in \omega ,$

причем

(3.2)

$w(x,t) = 0,\quad x \in \partial \omega .$

Наиболее интересная возможность аппроксимации со вторым порядком $h$ членов конвективного переноса связана с заданием скорости ${v}(x,t)$ в полуцелых узлах сетки $\bar {\omega }$. Оператор конвективного переноса в недивергентной форме (уравнение (1.1)) на множестве сеточных функций (3.2) зададим в виде

(3.3)

${{C}_{\infty }}w = \frac{1}{{2h}}v(x + 0.5h,t)(w(x + h,t) - w(x,t)) + \frac{1}{{2h}}v(x - 0.5h,t)(w(x,t) - w(x - h,t)),\quad x \in \omega .$

Аналогично сеточный оператор конвективного переноса в дивергентной форме (уравнение (1.4)) возьмем в виде

(3.4)

${{C}_{1}}w = \frac{1}{{2h}}v(x + 0.5h,t)(w(x + h,t) + w(x,t)) - \frac{1}{{2h}}v(x - 0.5h,t)(w(x,t) + w(x - h,t)),\quad x \in \omega .$

Для оператора конвективного переноса из уравнения (1.5) положим

${{C}_{2}} = \frac{1}{2}({{C}_{1}} + {{C}_{\infty }}),$

так что

(3.5)

${{C}_{2}}w = \frac{1}{{2h}}(v(x + 0.5h,t)w(x + h,t) - v(x - 0.5h,t)w(x - h,t))),\quad x \in \omega .$

После аппроксимации по пространству придем к задаче (2.2), (2.3), когда

(3.6)

$A = C + D,$

при задании $D$ согласно (3.1) на множестве функций (3.2), а $C = {{C}_{\infty }},{{C}_{1}},{{C}_{2}}$ – согласно (3.3), (3.4), (3.5) соответственно.

Сформулируем условия устойчивости и монотонности схем с весами (2.13), (2.14) для задачи (2.2), (2.3), (3.6). Свое рассмотрение начнем с задачи конвекции-диффузии (1.1)–(1.3) (аппроксимации (3.1)–(3.3), (3.6)).

Для того чтобы воспользоваться теоремами 1, 2, выпишем явно элементы матрицы $A$. Для (3.1)–(3.3), (3.6) ($C = {{C}_{\infty }}$) для диагональных элементов получим

${{a}_{{ii}}} = \frac{1}{{{{h}^{2}}}}({{k}_{{i + 1/2}}} + {{k}_{{i - 1/2}}}) - \frac{1}{{2h}}{{v}_{{i + 1/2}}} + \frac{1}{{2h}}{{v}_{{i - 1/2}}},$

а для внедиагональных имеем

(3.7)

${{a}_{{i,i - 1}}} = - \frac{1}{{{{h}^{2}}}}{{k}_{{i - 1/2}}} - \frac{1}{{2h}}{{v}_{{i - 1/2}}},\quad {{a}_{{i,i + 1}}} = - \frac{1}{{{{h}^{2}}}}{{k}_{{i + 1/2}}} + \frac{1}{{2h}}{{v}_{{i + 1/2}}}$

при использовании обозначений ${{k}_{{i \pm 1/2}}} = k(x \pm 0.5h),x \in \omega $.

Ключевое для монотонности условие неположительности внедиагональных элементов (2.22) дает

(3.8)

$\frac{1}{{{{h}^{2}}}}{{k}_{{i - 1/2}}} + \frac{1}{{2h}}{{v}_{{i - 1/2}}} \geqslant 0,\quad \frac{1}{{{{h}^{2}}}}{{k}_{{i + 1/2}}} - \frac{1}{{2h}}{{v}_{{i + 1/2}}} \geqslant 0.$

При неположительности внедиагональных элементов матрицы $A$ нестрогое диагональное преобладание по строкам (см. (2.9)) выполнено всегда. Неравенства (3.8) запишем в виде ограничений на шаг сетки по пространству:

(3.9)

$\frac{{h\left| {v(x \pm 0.5h,t)} \right|}}{{k(x \pm 0.5h)}} \leqslant 2,\quad x \in \omega .$

В случае конвективного переноса в дивергентной форме (3.1), (3.2), (3.4), (3.6) ($C = {{C}_{1}}$) для диагональных элементов матрицы $A$ имеем

${{a}_{{ii}}} = \frac{1}{{{{h}^{2}}}}({{k}_{{i + 1/2}}} + {{k}_{{i - 1/2}}}) + \frac{1}{{2h}}{{v}_{{i + 1/2}}} - \frac{1}{{2h}}{{v}_{{i - 1/2}}}.$

Внедиагональные элементы снова определяются согласно (3.7), а условия их неотрицательности (3.8) дают ограничения (3.9). Но теперь имеет место слабое диагональное преобладание по столбцам (см. (2.10)).

В наиболее интересном случае конвективного переноса в симметричной форме (3.1), (3.2), (3.5), (3.6) ($C = {{C}_{2}}$) для диагональных элементов матрицы $A$ получим выражение

${{a}_{{ii}}} = \frac{1}{{{{h}^{2}}}}({{k}_{{i + 1/2}}} + {{k}_{{i - 1/2}}}),$

с внедиагональными элементами (3.7). Условие диагонального преобладания выполнено в форме (2.12).

Итогом нашего рассмотрения является следующее утверждение, которое является следствием теорем 1, 2.

Теорема 3. Пусть в схеме (2.13), (2.14), (3.6) с сеточным оператором конвективного переноса $C = {{C}_{\infty }},{{C}_{1}},{{C}_{2}}$, определяемым согласно (3.3)–(3.5), выполнены условия (3.9). Тогда схема монотонна при всех $\tau > 0$, если $\sigma \geqslant 1$, а для разностного решения имеет место априорная оценка (2.15) в пространстве ${{L}_{\infty }}(\omega ),$ ${{L}_{1}}(\omega ),$ ${{L}_{2}}(\omega )$ соответственно.

Ограничение на шаг сетки по пространству можно снять, используя направленные разности для аппроксимации конвективных слагаемых. При этом мы жертвуем точностью: ранее рассмотренные центрально-разностные аппроксимации конвективных слагаемых имели второй порядок аппроксимации по $h$, а аппроксимации конвективного переноса направленными разностями – первый. Введем обозначения

$v(x,t) = {{v}^{ + }}(x,t) + {{v}^{ - }}(x,t),$

${{v}^{ + }}(x,t) = \frac{1}{2}(v(x,t) + \left| {v(x,t)} \right|{\kern 1pt} ) \geqslant 0,\quad {{v}^{ - }}(x,t) = \frac{1}{2}(v(x,t) - \left| {v(x,t)} \right|{\kern 1pt} ) \leqslant 0.$

Ориентируясь на задание скорости в полуцелых узлах, вместо (3.3) положим

(3.10)

${{C}_{\infty }}w = \frac{1}{h}{{v}^{ - }}(x + 0.5h,t)(w(x + h,t) - w(x,t)) + \frac{1}{h}{{v}^{ + }}(x - 0.5h,t)(w(x,t) - w(x - h,t)),\quad x \in \omega .$

Удобно записать (3.10), выделив явно диагональную часть оператора конвективного переноса:

(3.11)

$\begin{gathered} {{C}_{\infty }}w = \frac{1}{h}({{v}^{ + }}(x - 0.5h,t) - {{v}^{ - }}(x + 0.5h,t))w(x,t) + \\ + \;\frac{1}{h}({{v}^{ - }}(x + 0.5h,t)(w(x + h,t) - {{v}^{ + }}(x - 0.5h,t)w(x - h,t)),\quad x \in \omega . \\ \end{gathered} $

Для сеточного оператора конвективного переноса в дивергентной форме имеем

(3.12)

$\begin{gathered} {{C}_{1}}w = \frac{1}{h}({{v}^{ + }}(x + 0.5h,t) - {{v}^{ - }}(x - 0.5h,t))w(x,t) + \\ + \;\frac{1}{h}({{v}^{ - }}(x + 0.5h,t)(w(x + h,t) - {{v}^{ + }}(x - 0.5h,t)w(x - h,t)),\quad x \in \omega . \\ \end{gathered} $

Сеточный оператор конвективного переноса в симметричной форме аппроксимируется направленными разностями следующим образом:

(3.13)

$\begin{gathered} {{C}_{2}}w = \frac{1}{{2h}}({\kern 1pt} \left| {v(x + 0.5h,t)} \right| + \left| {v(x - 0.5h,t)} \right|{\kern 1pt} )w(x,t) + \\ + \;\frac{1}{h}({{v}^{ - }}(x + 0.5h,t)(w(x + h,t) - {{v}^{ + }}(x - 0.5h,t)w(x - h,t)),\quad x \in \omega . \\ \end{gathered} $

При таких аппроксимациях конвективных слагаемых внедиагональные элементы всегда неположительны, а диагональные – неотрицательны. Условия слабого диагонального преобладания в формах (2.9), (2.10) и (2.12) при аппроксимациях направленными разностями (3.11)–(3.13) проверяются непосредственно.

Теорема 4. Пусть в схеме (2.13), (2.14), (3.6) сеточный оператор конвективного переноса $C = {{C}_{\infty }},{{C}_{1}},{{C}_{2}}$ определяется согласно (3.11)–(3.13). Тогда схема монотонна при всех $\tau > 0$ и $\sigma \geqslant 1$, а для разностного решения имеет место априорная оценка (2.15) в пространстве ${{L}_{\infty }}(\omega ),$ ${{L}_{1}}(\omega ),$ ${{L}_{2}}(\omega )$ соответственно.

4. ЭКСПОНЕНЦИАЛЬНЫЕ СХЕМЫ

Безусловно монотонные схемы наиболее просто конструируются на основе трансформации исходных уравнений конвекции-диффузии с исключением членов конвективного переноса [4]. От (1.1) можно прийти к уравнению

(4.1)

$\frac{{\partial u}}{{\partial t}} - \frac{1}{{\chi (x,t)}}\frac{\partial }{{\partial x}}\left( {k(x)\chi (x,t)\frac{{\partial u}}{{\partial x}}} \right) = f(x,t),$

в котором

(4.2)

$\chi (x,t) = exp\left( { - \int\limits_0^x \,\frac{{{v}(s,t)}}{{k(s)}}ds} \right).$

С использованием функции $\chi (x,t)$ уравнение (1.4) можно записать в виде

(4.3)

$\frac{{\partial u}}{{\partial t}} - \frac{\partial }{{\partial x}}\left( {\frac{{k(x)}}{{\chi (x,t)}}\frac{{\partial (\chi (x,t)u)}}{{\partial x}}} \right) = f(x,t).$

Для уравнения (1.5) имеем

(4.4)

$\frac{{\partial u}}{{\partial t}} - \frac{1}{{2\chi (x,t)}}\frac{\partial }{{\partial x}}\left( {k(x)\chi (x,t)\frac{{\partial u}}{{\partial x}}} \right) - \frac{1}{2}\frac{\partial }{{\partial x}}\left( {\frac{{k(x)}}{{\chi (x,t)}}\frac{{\partial (\chi (x,t)u)}}{{\partial x}}} \right) = f(x,t).$

Отменим некоторые возможности по разностной аппроксимации таких диффузионно-конвективных слагаемых.

Будем отталкиваться от уравнения (4.1). Для сеточных функций, удовлетворяющих (3.2), положим

(4.5)

$\begin{gathered} Aw = - \frac{1}{{{{h}^{2}}\chi (x,t)}}k(x + 0.5h)\chi (x + 0.5h,t)(w(x + h,t) - w(x,t)) + \\ \, + \frac{1}{{{{h}^{2}}\chi (x,t)}}k(x - 0.5h)\chi (x - 0.5h,t)(w(x,t) - w(x - h,t)). \\ \end{gathered} $

Для вычисления $\chi (x \pm 0.5h,t)$ учитываем то, что

$\chi (x \pm 0.5h,t) = \chi (x)exp\left( { - \int\limits_x^{x \pm 0.5h} \,\frac{{{v}(s,t)}}{{k(s)}}ds} \right),$

и сеточные функции $k(x),\;{v}(x,t)$ заданы в полуцелых узлах. С точностью до $\mathcal{O}({{h}^{2}})$ положим

$\chi (x \pm 0.5h,t) = \chi (x)exp( \mp \theta (x \pm 0.5h,t)h)$

при использовании обозначений

$\theta (x,t) = \frac{{{v}(x,t)}}{{2k(x)}}.$

Это позволяет вместо (4.5) использовать следующую аппроксимацию для $A = {{A}_{\infty }}$:

(4.6)

$\begin{gathered} {{A}_{\infty }}w = - \frac{1}{{{{h}^{2}}}}k(x + 0.5h)exp( - \theta (x + 0.5h,t)h)(w(x + h,t) - w(x,t)) + \\ \, + \frac{1}{{{{h}^{2}}}}k(x - 0.5h)exp(\theta (x - 0.5h,t)h)(w(x,t) - w(x - h,t)). \\ \end{gathered} $

По аналогии с (4.5) для уравнения (4.3) используем

$\begin{gathered} Aw = - \frac{{k(x + 0.5h)}}{{{{h}^{2}}\chi (x + 0.5h,t)}}(\chi (x + h,t)w(x + h,t) - \chi (x,t)w(x,t)) + \\ \, + \frac{{k(x - 0.5h)}}{{{{h}^{2}}\chi (x - 0.5h,t)}}(\chi (x,t)w(x,t) - \chi (x - h,t)w(x - h,t)). \\ \end{gathered} $

Упрощая это выражение, получаем для $A = {{A}_{1}}$:

(4.7)

$\begin{gathered} {{A}_{1}}w = - \frac{1}{{{{h}^{2}}}}k(x + 0.5h)exp( - \theta (x + 0.5h,t)h)w(x + h,t) + \frac{1}{{{{h}^{2}}}}k(x + 0.5h)exp(\theta (x + 0.5h,t)h)w(x,t) + \\ \, + \frac{1}{{{{h}^{2}}}}k(x - 0.5h)exp( - \theta (x - 0.5h,t)h)w(x,t) - \frac{1}{{{{h}^{2}}}}k(x - 0.5h)exp(\theta (x - 0.5h,t)h)w(x - h,t). \\ \end{gathered} $

Для уравнения (4.4) используем аппроксимацию для ${{A}_{2}} = 0.5({{A}_{1}} + {{A}_{\infty }})$. С учетом (4.6) и (4.7) получим

(4.8)

$\begin{gathered} {{A}_{2}}w = - \frac{1}{{{{h}^{2}}}}k(x + 0.5h)exp( - \theta (x + 0.5h,t)h)w(x + h,t) + \\ \, + \frac{1}{{2{{h}^{2}}}}k(x + 0.5h)(exp(\theta (x + 0.5h,t)h) + exp( - \theta (x + 0.5h,t)h))w(x,t) + \\ \, + \frac{1}{{2{{h}^{2}}}}k(x - 0.5h)(exp(\theta (x - 0.5h,t)h) + exp( - \theta (x - 0.5h,t)h))w(x,t) - \\ \, - \frac{1}{{{{h}^{2}}}}k(x - 0.5h)exp(\theta (x - 0.5h,t)h)w(x - h,t). \\ \end{gathered} $

Для внедиагональных элементов $A = {{A}_{1}},{{A}_{2}},{{A}_{\infty }}$ имеем

${{a}_{{i,i - 1}}} = - \frac{1}{{{{h}^{2}}}}{{k}_{{i - 1/2}}}exp({{\theta }_{{i - 1/2}}}),\quad {{a}_{{i,i + 1}}} = - \frac{1}{{{{h}^{2}}}}{{k}_{{i + 1/2}}}exp( - {{\theta }_{{i + 1/2}}}),$

т.е. они все неотрицательны. Для диагональных элементов матрицы $A = {{A}_{\infty }}$ из (4.6) следует слабое диагональное преобладание по строкам (2.9). Аналогично, из (4.7) для $A = {{A}_{1}}$ имеем слабое диагональное преобладание по столбцам (2.10). Для $A = {{A}_{2}}$ из (4.8) получим слабое диагональное преобладание в форме (2.12). Тем самым можно сформулировать утверждение.

Теорема 5. Пусть в схеме (2.13), (2.14) сеточный оператор $A = {{A}_{\infty }},{{A}_{1}},{{A}_{2}}$ определяется согласно (4.6)–(4.8). Тогда схема монотонна при всех $\tau > 0$ и $\sigma \geqslant 1$, а для разностного решения имеет место априорная оценка (2.15) в пространстве ${{L}_{\infty }}(\omega ),$ ${{L}_{1}}(\omega ),$ ${{L}_{2}}(\omega )$ соответственно.

Как и для схем с аппроксимацией конвективного переноса направленными разностями (теорема 4) мы имеем безусловно устойчивые и монотонные в различных пространствах. Преимущество экспоненциальных схем состоит в том, что как и для схем с центрально-разностными аппроксимациями (теорема 3), погрешность аппроксимации по пространству имеет второй порядок. Платой за повышение точности является более сложное вычисление коэффициентов разностного оператора.

5. МНОГОМЕРНЫЕ ЗАДАЧИ

Отметим возможности построения монотонных схем для нестационарных уравнений конвекции-диффузии на примере модельных двумерных задач. В прямоугольнике

$\Omega = \{ {\mathbf{x}}\,{\text{|}}\,{\mathbf{x}} = ({{x}_{1}},{{x}_{2}}),\;0 < {{x}_{\alpha }} < {{l}_{\alpha }},\;\alpha = 1,2\} $

рассматривается нестационарное уравнение конвекции-диффузии с конвективными слагаемыми в недивергентном виде

(5.1)

$\frac{{\partial u}}{{\partial t}} + \sum\limits_{\alpha = 1}^2 \,{{{v}}_{\alpha }}({\mathbf{x}},t)\frac{{\partial u}}{{\partial {{x}_{\alpha }}}} - \sum\limits_{\alpha = 1}^2 \,\frac{\partial }{{\partial {{x}_{\alpha }}}}\left( {k({\mathbf{x}})\frac{{\partial u}}{{\partial {{x}_{\alpha }}}}} \right) = f({\mathbf{x}},t).$

Это уравнение дополняется однородными граничными условиями Дирихле

(5.2)

$u({\mathbf{x}},t) = 0,\quad {\mathbf{x}} \in \partial \Omega ,\quad t > 0.$

Кроме того, задается начальное условие

(5.3)

$u({\mathbf{x}},0) = {{u}^{0}}({\mathbf{x}}),\quad {\mathbf{x}} \in \Omega .$

Как и при рассмотрении одномерных задач, вторым примером является нестационарное уравнение конвекции-диффузии с конвективным переносом в дивергентной форме:

$\frac{{\partial u}}{{\partial t}} + \sum\limits_{\alpha = 1}^2 \,\frac{\partial }{{\partial {{x}_{\alpha }}}}({{v}_{\alpha }}({\mathbf{x}},t)u) - \sum\limits_{\alpha = 1}^2 \,\frac{\partial }{{\partial {{x}_{\alpha }}}}\left( {k({\mathbf{x}})\frac{{\partial u}}{{\partial {{x}_{\alpha }}}}} \right) = f({\mathbf{x}},t).$

Основным является уравнение конвекции-диффузии с конвективным переносом в симметричной форме, когда

$\frac{{\partial u}}{{\partial t}} + \frac{1}{2}\sum\limits_{\alpha = 1}^2 \,\left( {\frac{\partial }{{\partial {{x}_{\alpha }}}}({{{v}}_{\alpha }}({\mathbf{x}},t)u) + {{{v}}_{\alpha }}({\mathbf{x}},t)\frac{{\partial u}}{{\partial {{x}_{\alpha }}}}} \right) - \sum\limits_{\alpha = 1}^2 \,\frac{\partial }{{\partial {{x}_{\alpha }}}}\left( {k({\mathbf{x}})\frac{{\partial u}}{{\partial {{x}_{\alpha }}}}} \right) = f({\mathbf{x}},t).$

Операторы конвекции-диффузии в многомерных задачах представляются в виде суммы одномерных операторов конвекции-диффузии. В силу этого при построении монотонных схем для многомерных задач мы можем опираться на рассмотренные выше аппроксимации одномерных операторов конвекции-диффузии. В частности, можно строить экспоненциальные схемы на основе записи уравнений аналогично (4.1), (4.3), (4.4), например, для задачи (5.1)–(5.3) будем использовать уравнение

$\frac{{\partial u}}{{\partial t}} - \sum\limits_{\alpha = 1}^2 \,\frac{1}{{{{\chi }_{\alpha }}({\mathbf{x}},t)}}\frac{\partial }{{\partial {{x}_{\alpha }}}}\left( {k({\mathbf{x}}){{\chi }_{\alpha }}({\mathbf{x}},t)\frac{{\partial u}}{{\partial {{x}_{\alpha }}}}} \right) = f({\mathbf{x}},t),$

где теперь

${{\chi }_{1}}({\mathbf{x}},t) = exp\left( { - \int\limits_0^{{{x}_{1}}} \,\frac{{{{{v}}_{1}}(s,{{x}_{2}},t)}}{{k(s,{{x}_{2}})}}ds} \right),\quad {{\chi }_{2}}({\mathbf{x}},t) = exp\left( { - \int\limits_0^{{{x}_{2}}} \,\frac{{{{{v}}_{2}}({{x}_{1}},s,t)}}{{k({{x}_{1}},s)}}ds} \right).$

При использовании равномерных по каждому направлению сеток исследование устойчивости и монотонности разностных схем для многомерных уравнений конвекции-диффузии в соответствующих пространствах сеточных функций проводится полностью аналогично одномерному случаю.

При расщеплении по направлениям монотонность обеспечивается использованием локально-одномерных схем покомпонентного расщепления [26]. Такое рассмотрение проведено в [15] для многомерных задач конвекции-диффузии, когда конвективный перенос берется в недивергентной или дивергентной формах.

Не принципиальным является также обобщение результатов на случай неравномерных прямоугольных сеток. Более интересные проблемы порождаются использованием общих неструктурированных сеток. В этой связи отметим, что монотонные схемы метода конечных объемов строятся [27] при использовании сеток Делоне и разбиения Вороного, как контрольного объема. В методе конечных элементов мы можем использовать линейные конечные элементы на таких сетках и специальные аппроксимации коэффициента при производной по времени – lumping процедуры (см., например, [7]).

Список литературы

Hundsdorfer W.H. Numerical Solution of Time-Dependent Advection-Diffusion-Reaction Equations. Berlin: Springer, 2003.
Morton K.W. Numerical Solution of Convection-Diffusion Problems. London: Chapman & Hall, 1996.
Wesseling P. Principles of Computational Fluid Dynamics. Berlin: Springer, 2001.
Самарский А.А., Вабищевич П.Н. Численные методы решения задач конвекции–диффузии. М.: URSS, 1999.
Vabishchevich P.N. On the form of the hydrodynamics equations // High Speed Flow Field conference, 19–22, November 2007, Moscow, Russia. Moscow, 2007. P. 1–9.
Brenner S., Scott R.The Mathematical Theory of Finite Element Methods. New York: Springer, 2007.
Thomée V. Galerkin Finite Element Methods for Parabolic Problems. Berlin: Springer, 2006.
Самарский А.А. Теория разностных схем. М.: Наука, 1989.
Самарский А.А., Гулин А.В. Устойчивость разностных схем. М.: Наука, 1973.
Protter M.H., Weinberger H.F. Maximum Principles in Differential Equations. New York: Springer, 1967.
Ладыженская О.А., Солонников В.А., Уральцева Н.Н. Линейные и квазилинейные уравнения параболического типа М.: Наука, 1967.
Tannehill J.C., Anderson D.A., Pletcher R. H. Computational Fluid Mechanics and Heat Transfer. Philadelphia: Taylor & Francis, 1997.
de G.Allen D.N., Southwell R.V. Relaxation methods applied to determine the motion, in two dimensions, of a viscous fluid past a fixed cylinder // Quart. J. Mech. Appl. Math. 1955. V. 8. P. 129–145.
Scharfetter D.L., Gummel H.K. Large-signal analysis of a silicon read diode oscillator // IEEE Trans. Electron Devices. 1969. V. ED-16. P. 4–77.
Afanas’eva N.M., Churbanov A.G., Vabishchevich P.N. Unconditionally monotone schemes for unsteady convection-diffusion problems // Comput. Meth. in Appl. Math. 2013. V. 13. № 2. P. 185–205.
Dekker K., Verwer J.G. Stability of Runge–Kutta Methods for Stiff Nonlinear Differential Equations. Amsterdam: North-Holland, 1984.
Desoer C.A., Vidyasagar M. Feedback Systems: Input-Output Properties. Philadelphia: SIAM, 2008.
Friedman A. Partial Differential Equations of Parabolic Type. Englewood: Prentice-Hall, 1964.
Horn R.A., Johnson C.R. Matrix Analysis. Cambridge: Cambridge University Press, 1990.
Dekker K., Verwer J.G. Stability of Runge–Kutta Methods for Stiff Nonlinear Differential Equations. Amsterdam: North-Holland, 1984.
Hairer E., Norsett S.P., Wanner G. Solving Ordinary Differential Equations. I. Nonstiff Problems. Berlin: Springer, 1987.
Вабищевич П.Н. Потоковые схемы расщепления для параболических задач со смешанными производными // Ж. вычисл. матем. и матем. физ. 2013. Т. 53. № 8. С. 1314–1328.
Sun Y.X. Sufficient conditions for generalized diagonally dominant matrices // Numerical Math. (A Journal of Chinese Universities). 1997. V. 19. № 3. P. 216–223. in Chinese.
Wang G., Hong Z., Gao Z. Sufficient conditions of nonsingular H-matrices // J. of Shanghai University (English Edition). 2004. V. 8. № 1. P. 35–37.
Guo Z.J., Guo Z.J., Yang J.G. A new criteria for a matrix is not generalized strictly diagonally dominant matrix // Applied Mathematical Sciences. 2011. V. 5. № 6. P. 273–278.
Вабищевич П.Н. Аддитивные операторно-разностные схемы (схемы расщепления). М.: URSS, 2013.
Vabishchevich P.N. Finite-difference approximation of mathematical physics problems on irregular grids // Computational Methods in Applied Mathematics. 2005. V. 5. № 3. P. 294–330.

Дополнительные материалы отсутствуют.

Инструменты

следующая статья выпуска предыдущая статья выпуска содержание выпуска

Журнал вычислительной математики и математической физики

Архивы выпусков Информация о журнале Отправить рукопись в журнал