Монотонизация модифицированной схемы с эрмитовой интерполяцией для численного решения неоднородного уравнения переноса с поглощением

Введение

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

Разнообразие областей знаний, для развития которых требуется решение уравнения переноса, приводит к ещё большему многообразию способов численного решения уравнения переноса. Предлагаемые разностные схемы различны по своим свойствам, так как требования, которым должны удовлетворять схемы, могут сильно отличаться для разных классов задач.

Общей идеей, определяющей развитие разностных схем в настоящий момент, является стремление построить высокоточные схемы, позволяющие получать правильное решение даже на неподробных сетках. Дополнительным требованием при решении задач переноса излучения является возможность расчёта с крупным шагом по времени, поэтому предпочтение отдаётся неявным безусловно устойчивым схемам. Использование схем с компактным шаблоном упрощает построение решения вблизи границ расчётной области и в окрестностях точек разрыва коэффициентов уравнения.

Развитие схем с порядком аппроксимации выше первого повлекло за собой развитие способов монотонизации получаемого решения, так как немонотонные схемы не обладают свойством положительности, а искомое решение с физической точки зрения представляет собой функцию распределения частиц (плотность в фазовом пространстве), которая положительна по определению. Нефизичное численное решение с отрицательными значениями искомой функции на некотором наборе узлов создаёт целый ряд сложностей. Например, встаёт вопрос о корректном вычислении величин, получаемых осреднением по угловой переменной.

Для построения монотонных схем высокого порядка аппроксимации могут использоваться различные подходы. Созданы различные схемы с ограничителями потоков (TVD схемы – Total Variation Diminishing [4]), схемы с переключениями шаблонов (ENO схемы – Essentially Non-Oscillatory, WENO схемы – Weighted Essentially Non-Oscillatory [5-7]), гибридные схемы, в которых решение представляет собой среднее взвешенное решений по схемам высокого и низкого порядков [8]. Гибридные схемы, сохраняющие интегральные средние по ячейке, были разработаны в [9].

В данной работе рассматривается одна из модификаций схемы с эрмитовой интерполяцией – схема CIP (Cubic Interpolation Polynomial), которая впервые была предложена в [10-16]. Это интерполяционно-характеристическая схема третьего порядка аппроксимации по времени и пространству. Интерполянт Эрмита строится по узловым значениям функции и её производных, поэтому при построении схемы указывается способ пересчёта узловых значений, которые переносятся вдоль характеристики, выпущенной назад, а также способ пересчёта производных. Рассматриваемые модификации схемы CIP отличаются способом пересчёта производных. В работах [10-16] для пересчёта производных схема применяется не только для самого уравнения переноса, но и для его дифференциального следствия. Мы будем рассматривать модификацию схемы, в которой производные вычисляются через интегральные средние и формулу Эйлера–Маклорена [17, 18]. Также используется проекционный способ замыкания схемы CIP [19].

Модифицированная схема CIP была исследована авторами в применении к решению неоднородного уравнения переноса с постоянным и кусочно-постоянным коэффициентом поглощения [20]. В одной из предыдущих работ была предложена консервативная монотонизация модифицированной схемы CIP [21]. В настоящей работе рассматривается другой способ монотонизации, конструируется гибридная схема, в которой расчёт решения производится по схеме CIP и по характеристической схеме первого порядка, после чего строится гибрид.

В работе [20] для вычисления интеграла вдоль характеристики применяется формула Симпсона четвёртого порядка аппроксимации для интеграла Римана, что в случае больших оптических толщин приводит к значительной ошибке численного интегрирования. В данной работе формула Симпсона применяется к интегралу вдоль характеристики, которому придаётся вид интеграла Стилтьеса.

Работа имеет следующую структуру. В первом разделе дана математическая формулировка рассматриваемой дифференциальной задачи. Второй раздел содержит краткое описание модифицированной схемы CIP и вычислительные формулы. В третьем разделе изложена процедура гибридизации схемы CIP; при этом рассмотрены три варианта монотонизации: локальная, послойная и глобальная. Далее проведено тестирование гибридных схем при различных числах Куранта. В четвёртом разделе приведены результаты тестов для уравнения адвекции, а в пятом – для неоднородного уравнения переноса с поглощением. Шестой раздел содержит описание нового варианта интегрирования вдоль характеристики, использующего интеграл Стилтьеса, для уменьшения ошибки численного интегрирования в случае больших оптических толщин. Результаты тестирования этого варианта интегрирования приведены в седьмом разделе.

1. Постановка задачи

Рассмотрим линейное неоднородное уравнение переноса

в котором $u=u(x,t)$ – искомая функция, зависящая от единственной пространственной переменной $x$ и времени $t$. Правую часть $f=f(x,t)$ и коэффициент поглощения $\kappa =\kappa (x,t)$ будем считать известными функциями. Кроме того, ограничимся случаем положительной скорости переноса $a>0$.

Для постановки начально-краевой задачи дополним уравнение (1) начальными условиями Коши

и классическими граничными условиями при $x=0$:

В дальнейшем наряду с граничными условиями (3) мы будем использовать также периодические граничные условия вида

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

Воспользуемся характеристическими свойствами уравнения (1) и выполним замену переменных

Тогда при $\xi =x-at=\text{const}$ уравнение (1) сведётся к обыкновенному дифференциальному уравнению

Будем считать коэффициент поглощения постоянным: $\kappa (x,t)={\kappa }_0=\text{const}$. В этом случае решение уравнения (5) есть функция

Здесь $u^{\ast }=u({\eta }^{\ast })$, где ${\eta }^{\ast }=x^{\ast }+a{t}^{\ast }$; точка $(x^{\ast },t^{\ast })$ принадлежит характеристике.

В частном случае ${\kappa }_0=0$, $f\equiv 0$ уравнение переноса (1) принимает вид уравнения адвекции

для которого, согласно (6),

то есть величина $u$ переносится вдоль характеристики без изменений.

2. Построение разностной схемы

Для численного решения начально-краевой задачи (1-3) или (1, 2, 4) будем использовать модифицированную схему CIP. Принцип её построения подробно изложен в работе [17]. Ограничимся кратким описанием схемы с указанием основной идеи построения и вычислительных формул.

Введём равномерную сетку с шагами $h$ по пространству и $\tau$ по времени (рис. 1). Пусть число Куранта $\sigma =a\tau /h\le 1$. Из точки A выпустим характеристику $x-at=\text{const}$, которая пересечёт нижнюю границу ячейки в точке $E=(x^{\ast },t^n)$.

На отрезке CB введём барицентрические координаты

и рассмотрим базисные функции Эрмита:

Кубический интерполянт запишем в виде

где $y_m^n$ и $d_m^n$ – сеточные значения искомого решения и его пространственной производной соответственно.

Из формулы (9) можно найти $y^{\ast }$ – значение сеточной функции в точке E. Для уравнения адвекции (7), в силу равенства (8),

При решении неоднородного уравнения (1) для нахождения $y_{m+1}^{n+1}$ необходимо кроме $y^{\ast }$ вычислить интеграл из формулы (6) вдоль отрезка характеристики EA. Поскольку функция $f$ произвольна, будем использовать квадратурную формулу Симпсона для нахождения интеграла из (6), что при умеренных оптических толщинах ячейки не снижает точности расчёта [20].

Поскольку в (9) входят сеточные значения производной $d_m^n$, для замыкания алгоритма схемы CIP необходимо указать способ их вычисления при переходе к новому слою по времени. С этой целью определим интегральное среднее по отрезку BA:

Воспользуемся формулой Эйлера–Маклорена

в которой $g_{m+1}^{n+1}$ и $g_{m+1}^n$ – сеточные значения производной по времени, связанные с сеточными значениями пространственной производной равенствами

Из характеристического свойства (8) уравнения адвекции (7) следует, что интегральное среднее (11) совпадает со средним по отрезку EB:

Интеграл (14) может быть вычислен аналитически или приближённо по формуле Симпсона. В более общем случае уравнения (1) интегральные средние по отрезкам BA и EB не совпадают, поэтому ${\bar{y}}_{m+1}$ из (11) вычисляется также по формуле Симпсона:

Величину $y_{m+1}^{n+1/2}$ можно найти по формуле (6), записанной для отрезка характеристики GF.

Отметим, что для неоднородного уравнения (1) связь пространственных и временных производных принимает вид

Итак, в случае уравнения адвекции (7) после расчёта величины $y_{m+1}^{n+1}$ из соотношения (10) с помощью формул (9) и (14) вычисляется интегральное среднее ${\bar{y}}_{m+1}$. Далее из формулы (12) выражается производная по времени $g_{m+1}^{n+1}$, которая пересчитывается в производную по пространству согласно (13). В случае неоднородного уравнения (1) сначала из формулы (9) находится $\hat{y}$ – значение сеточной функции в точке G. Затем по $\hat{y}$ с помощью формулы (6) для отрезка характеристики GF рассчитывается $y_{m+1}^{n+1/2}$. Далее используется (15) и вычисляется ${\bar{y}}_{m+1}$. Сеточное значение производной $g_{m+1}^{n+1}$ находится из (12), что позволяет вычислить $d_{m+1}^{n+1}$ с помощью (16).

Приведём окончательный вид вычислительных формул модифицированной схемы CIP при $\sigma \le 1$. Для уравнения адвекции (7) имеем:

для неоднородного уравнения (1) формулы принимают вид

При $\sigma >1$ модифицированная схема CIP строится аналогично. В этом случае характеристика EA пересекает левую границу ячейки, а не нижнюю (рис. 2), и поэтому эрмитова интерполяция применяется на отрезке CD, причём используются производные по времени, а не по пространству. Формула Эйлера–Маклорена записывается в виде

Вычислительные формулы схемы CIP для случая $\sigma >1$ можно получить из представленных выше формул для случая $\sigma \le 1$, если заменить $\sigma$ на ${\sigma }^{-1}$, $m$ – на $n$, $n$ – на $m$, а сеточную функцию ${\tilde{d}}_m^n$ заменить функцией ${\tilde{g}}_m^n\equiv \tau g_m^n$. Аргументы функции $f$ видоизменяются так, чтобы соответствовать новым точкам пересечения характеристик с границами ячейки (рис. 2).

Ранее модифицированная схема CIP была применена для решения уравнения адвекции (7) и неоднородного уравнения (1) при $\sigma \le 1$. Был продемонстрирован третий порядок сходимости по времени и пространству на гладких тестах [17]. Подробные результаты тестирования схемы содержатся в работе [20].

3. Монотонизация модифицированной схемы CIP

При численном решении уравнения переноса обычно отдают предпочтение монотонным схемам, поскольку они позволяют избежать нефизичных решений, а для немонотонных схем высокого порядка аппроксимации предлагают способы монотонизации. Ранее была предложена консервативная монотонизация схемы CIP [21], позволяющая получить монотонный профиль решения при выполнении ряда условий на сеточные значения производных. Мы рассматриваем другой вариант монотонизации, который основан на использовании гибридных схем и не требует выполнения дополнительных условий.

Построим гибридную схему, следуя работе [8]. Пусть $y^{(\text{CIP})}$ – решение разностной задачи для уравнения адвекции (7) или неоднородного уравнения (1), полученное с помощью схемы CIP, а $y^{(\text{1})}$ – решение той же разностной задачи по характеристической схеме низкого порядка, построенной на шаблоне явного (при $\sigma \le 1$) или неявного (при $\sigma >1$) противопоточного уголка. Тогда гибридное решение вычисляется по формуле

в которой

Натуральный параметр $m$ из формулы для весового коэффициента $\chi$ должен удовлетворять условию $m\ge m_0-1$, где $m_0=3$ – порядок аппроксимации схемы CIP по времени. Таким образом, $m\ge 2$. Положительная константа $C_1$ зависит от решаемой задачи и в каждом конкретном случае подбирается так, чтобы обеспечить наилучший результат монотонизации.

В формулу для расчёта производной (17) входит величина $y_{m+1}^{n+1}$, в качестве которой можно использовать значение ${\left[y^{(\text{CIP})}\right]}_{m+1}^{n+1}$, рассчитанное по схеме CIP до гибридизации, или значение ${\left[y^{(\text{hyb})}\right]}_{m+1}^{n+1}$, полученное в ходе гибридизации. Первый способ приводит к меньшим погрешностям численного решения, что будет обосновано в следующем разделе.

Гибридизацию можно проводить на разных этапах расчёта. В первом варианте гибридное решение вычисляется по формуле (18) сразу после вычисления $y^{(\text{CIP})}$ и $y^{(\text{1})}$ для текущей ячейки. Второй вариант сводится к вычислению $y^{(\text{CIP})}$ и $y^{(\text{1})}$ для целого временно́го слоя с последующей гибридизацией. Наконец, в третьем варианте гибридизация проводится после предварительного расчёта $y^{(\text{CIP})}$ и $y^{(\text{1})}$ для всех слоёв по времени. Перечисленные варианты будем называть локальной, послойной и глобальной монотонизацией соответственно.

При $\sigma \le 1$ локальная и послойная монотонизации эквивалентны, а при $\sigma >1$ все три варианта различны. Таким образом, при $\sigma \le 1$ будем рассматривать локальную и глобальную монотонизацию, а при $\sigma >1$ – локальную, послойную и глобальную.

4. Тестирование гибридных схем для уравнения адвекции

Проведём тестирование гибридных схем для уравнения адвекции (7), дополненного начальными условиями (2) и граничными условиями (3), выбранными так, что

В этом случае точное решение задаётся функцией

Поскольку список неизвестных в схеме CIP включает сеточные значения производных, для них также зададим начальные и граничные условия. При $\sigma \le 1$ имеем:

аналогично при $\sigma >1$

Для исследования порядков сходимости гибридных схем предложим набор из четырёх тестов.

В тестах 2,3,4 вместо классических граничных условий (3) используются периодические граничные условия (4). При этом точное решение также имеет вид (19).

Сравним результаты численного решения уравнения (7) при использовании двух способов расчёта производной, описанных в предыдущем разделе. В одном варианте расчёта используется гибридное решение, а в другом – нет. Начнём с варианта, в котором $y_{m+1}^{n+1}={\left[y^{(\text{CIP})}\right]}_{m+1}^{n+1}$.

Положим $a=1$, $\sigma =0,5$, $X=1$, $T=1$. Тогда $N=2M$. Погрешности численного решения $\epsilon =\Vert {[u]}_{\tau h}-y\Vert$ и порядки сходимости $p$ схемы CIP с локальной монотонизацией на различных сетках приведены в табл. 1. При вычислениях используется одна из трёх норм: $C$, $L_1$ или $L_2$. Для теста 4 с разрывной функцией определение порядка сходимости в норме $C$ не является корректным, и поэтому для данного теста использованы только нормы $L_1$ и $L_2$.

Во всех тестах мы положили $m=4$, что удовлетворяет условию $m\ge 2$ из предыдущего раздела. Константа ${С}_1$ подбиралась так, чтобы профиль численного решения был близок к монотонному при возможно более высоком порядке сходимости.

Было установлено, что наилучшая монотонизация достигается при $C_1=0,70$ для тестов 2,3,4. В случае теста 1 в силу гладкости начального профиля $u_0(x)$ численное решение по схеме CIP не имело немонотонностей, так что при ${С}_1\in \left(0,100\right]$ величины ошибок и порядки сходимости существенно не изменялись; поэтому для теста 1 мы выбрали ${С}_1=1,00$.

Профили численных решений для разных тестов показаны на рис. 3. Для тестов 2 и 3 получены монотонные профили решений (рис. 3а-3б, 3в-3г). В тесте 4 гибридизация не позволила устранить все немонотонности (рис. 3ж): немонотонность профиля в полуокрестностях точек разрыва аналитического решения, где решение было отлично от нуля, сохранялась при любых значениях параметров $m$ и ${С}_1$.

По данным табл. 1 построены графики зависимостей $\epsilon =\epsilon (M)$ в логарифмических координатах (рис. 5). Пунктиром на этих графиках изображены прямые с коэффициентами, полученными методом наименьших квадратов (МНК). На рис. 5а для наглядности проведена прямая с угловым коэффициентом $k=-3$. Порядки сходимости численного решения по гибридной схеме, найденные по МНК, приведены в табл. 2. На гладком тесте продемонстрирован третий порядок сходимости во всех нормах (табл. 1, рис. 5а). При снижении гладкости функции в тестах 2,3,4 порядок сходимости оказывается меньше порядка аппроксимации (табл. 1, рис. 5б, 5в, 5г), как и ожидалось [20]. Приведённые в работе [20] порядки сходимости схемы CIP без монотонизации не отличаются от результатов из табл. 2.

Тест 1 с гладкой функцией был проведён также для схемы CIP с глобальной монотонизацией. Полученные погрешности и порядки сходимости представлены в табл. 3, которая демонстрирует, что порядок сходимости гибридной схемы равен 2 в норме $C$ и 2,47 в норме $L_2$ уже на гладком решении. Этот вывод оказывается справедливым при любых $m$ и $C_1$. В то же время схема с локальной монотонизацией на гладком тесте демонстрирует третий порядок сходимости во всех нормах (табл. 2).

На рис. 6в, 6г представлен профиль решения по схеме с глобальной монотонизацией для теста 2. Из рисунка видно, что пиковое значение решения существенно уменьшается уже за один период. Также в гибридной схеме проявляется значительная диссипация слева и справа от пика. Из сказанного следует, что локальная монотонизация более эффективна по сравнению с глобальной. Аналогичный вывод справедлив для тестов 3,4.

Рассмотрим теперь вариант расчёта производных с использованием гибридного решения, положим $y_{m+1}^{n+1}={\left[y^{(\text{hyb})}\right]}_{m+1}^{n+1}$ в формуле (17). В тесте 1 порядок сходимости равен 3, погрешности численного решения не отличаются существенно от погрешностей, приведённых в табл. 1. На рис. 4 представлен профиль численного решения для теста 2. Наблюдается значительное размытие профиля: пиковое значение решения уменьшается примерно на 70% за один период, размытие профиля оказывается даже более сильным, чем при глобальной монотонизации с $y_{m+1}^{n+1}={\left[y^{(\text{CIP})}\right]}_{m+1}^{n+1}$ (рис. 6в). Кроме того, в широкой области изменения параметров сетки и гибридизации сходимость отсутствует. Это обосновывает выбор в пользу расчёта производных по значению ${\left[y^{(\text{CIP})}\right]}_{m+1}^{n+1}$. Все дальнейшие тесты будут выполняться при $y_{m+1}^{n+1}={\left[y^{(\text{CIP})}\right]}_{m+1}^{n+1}$ в формуле (17).

Перейдём теперь к тестированию гибридных схем при $\sigma >1$. Воспользуемся теми же тестами 1,2,3,4, отказавшись от периодических граничных условий; таким образом, во всех четырёх случаях будем использовать классические граничные условия (3), согласованные с начальными условиями (2). Пусть $a=1$, $\sigma =2$, $T=1$. Отрезок по координате выберем вдвое бо́льшим по сравнению с тем же отрезком при $\sigma \le 1$: $X=2$, – чтобы распространение нетривиальной части профиля численного решения можно было наблюдать на единичном отрезке: от $x=0,5$ до $x=1,5$. При выбранных значениях параметров $M=4N$. Результаты тестирования схемы CIP с локальной монотонизацией приведены в табл. 4, 5 и на рис. 8. Порядки сходимости гибридной схемы существенно не отличаются от представленных в табл. 1 для случая $\sigma \le 1$.

Повторим те же тесты с послойной монотонизацией вместо локальной. Результаты приведены в табл. 6, 7 и на рис. 9. Сравнение с табл. 4, 5 показывает, что изменение способа гибридизации слабо влияет на величины ошибок и порядки сходимости схемы. В тесте 2 профиль решения по схеме с послойной монотонизацией является монотонным, но оказывается асимметричным. Наблюдается более сильное размытие профиля справа от пика (см. рис. 7в, 7г). Отсюда следует, что предпочтение следует отдавать локальной монотонизации, при которой профиль численного решения монотонен и симметричен. Аналогичный результат наблюдается в тесте 3.

Как и в случае $\sigma \le 1$, переход к гибридным схемам не устраняет всех немонотонностей на тесте 4. Кроме того, в тесте 1 при $\sigma >1$ решения по схеме CIP и схеме низкого порядка на каждом шаге по времени отличаются более значительно, чем при $\sigma \le 1$, поэтому наблюдается более сильная зависимость порядка сходимости гибридной схемы от параметров гибридизации; третий порядок сходимости достигается только при небольших значениях C₁.

5. Тестирование гибридных схем для неоднородного уравнения переноса

Перейдём к тестированию гибридных схем для неоднородного уравнения переноса (1). Будем использовать тот же набор тестов, что и в предыдущем разделе. Функцию правой части выберем в виде

Тогда точное решение, как и в случае уравнения адвекции (7), задаётся формулой (19).

Начальные и граничные условия на производную, использующиеся при численном решении неоднородного уравнения по схеме CIP, вообще говоря, отличаются от условий на производную (20) и (21) для уравнения адвекции. В случае $\sigma \le 1$ вместо (20) имеем условия

а в случае $\sigma >1$ вместо (21) – условия

Изменения связаны с тем, что при выводе начальных и граничных условий на функции ${\tilde{d}}_m^n$ и ${\tilde{g}}_m^n$ используется связь пространственных и временны́х производных, которая задаётся решаемым уравнением. Однако при выборе правой части в виде (22) указанное различие пропадает, так что в наших тестах начальные и граничные условия к уравнению адвекции (7) и неоднородному уравнению переноса (1) оказываются одинаковыми.

Как показано в предыдущем разделе, при решении уравнения адвекции наилучшие результаты даёт схема с локальной монотонизацией, поэтому, переходя к решению неоднородного уравнения, мы ограничиваемся этим вариантом гибридной схемы. В качестве схемы низкого порядка при вычислении гибридного решения будем использовать сеточно-характеристическую схему с линейной интерполяцией, в которой для расчёта интеграла из (6) применяется формула прямоугольников. Таким образом, схема низкого порядка (при $\sigma \le 1$) имеет вид:

Аналогичная формула справедлива и в случае $\sigma >1$:

Результаты тестов для случая $\sigma \le 1$ при периодических граничных условиях приведены в табл. 8, 9 и на рис. 10. В этих тестах $a=1$, ${\kappa }_0=1$, $\sigma =0,5$, $X=1$, $T=1$, а $N=2M$. Погрешности численного решения и порядки сходимости не отличаются существенно от представленных в табл. 1, 2.

В случае $\sigma >1$ полагаем $a=1$, ${\kappa }_0=1$, $\sigma =2$, $X=2$, $T=1$, при этом $M=4N$ и вместо периодических граничных условий (4) используются классические условия (3). Результаты тестов приведены в табл. 10, 11 и на рис. 11. Существенных отличий от результатов из табл. 4, 5 также не наблюдается.

6. Схема CIP для случая больших оптических толщин

При численном решении неоднородного уравнения переноса по схеме CIP, как и при решении по любой другой интерполяционно-характеристической схеме, необходимо вычислять интеграл из (6). Для этого в схеме CIP используется формула Симпсона четвёртого порядка аппроксимации. При этом важным параметром является безразмерная комбинация $l=\kappa s/a$ ( $s$ – характерный размер задачи), называемая оптической толщиной. До сих пор мы рассматривали случай умеренных оптических толщин, когда ошибки численного интегрирования по формуле Симпсона не превосходят ошибок интерполяции. Однако при больших оптических толщинах основной вклад в погрешность численного решения вносят ошибки численного интегрирования [20, 22].

Предложим использовать в схеме CIP новый вариант интегрирования вдоль характеристики, который позволит уменьшить ошибки численного решения в случае больших оптических толщин. Основная идея этого варианта заключается в том, что выполняется переход от интеграла Римана к интегралу Стилтьеса, в результате которого формула (6) принимает вид:

Здесь ${\gamma }^{\prime }({\eta }^{\prime })=\exp \left[-{\kappa }_0\left(\eta -{\eta }^{\prime }\right)/2a\right]$ – интегрирующая функция, ${\gamma }^{\ast }=\gamma ({\eta }^{\ast })=e^{-{\kappa }_0\tau }$. В данном случае $s=\eta -{\eta }^{\ast }$, что соответствует длине отрезка интегрирования вдоль характеристики EA (рис. 1).

После перехода к интегрированию по Стилтьесу изменяются только формулы для $y_{m+1}^{n+1}$ и $y_{m+1}^{n+1/2}$, в которые входит интеграл вдоль характеристики. При $\sigma \le 1$ эти формулы имеют вид

где ${\gamma }^{\ast }=e^{-{\kappa }_0\tau }$, $\hat{\gamma }=e^{-{\kappa }_0\tau /2}$, ${\gamma }_{\text{I}}=\left(1+{\gamma }^{\ast }\right)/2$, ${\gamma }_{\text{II}}=\left(1+\hat{\gamma }\right)/2$.

Исследуем приведённые формулы схемы CIP в пределе $l\gg 1$, когда ${\gamma }^{\ast }\approx 0$, $\hat{\gamma }\approx 0$. Применение формулы Симпсона к интегралу из (23) даёт

В отличие от формулы Симпсона для интеграла Римана второе слагаемое вычисляется не в середине K отрезка характеристики EA (рис. 1), а в точке $K^{\ast }$, смещённой от K на $\Delta x$ по координате и $\Delta t$ по времени. С увеличением оптической толщины ячейки $l$ точка $K^{\ast }$ сдвигается вдоль характеристики к точке A (поскольку $\underset{l\to \infty }{\lim }L=0$). В пределе $l\to \infty$ обе точки ( $K^{\ast }$и A) совпадают.

Аналогичным образом можно перейти к интегралу Стилтьеса и в случае $\sigma >1$.

Далее для краткости будем использовать следующие обозначения: CIP_R – схема CIP, в которой формула Симпсона применяется к интегралу Римана, CIP_S – схема CIP с интегралом Стилтьеса.

Схемой низкого порядка, как и раньше, является сеточно-характеристическая схема первого порядка на шаблоне уголка, в которой формула прямоугольников теперь применяется к интегралу Стилтьеса. Таким образом, при $\sigma \le 1$ схема низкого порядка имеет вид:

Аналогичную формулу можно получить и для случая $\sigma >1$.

7. Тестирование схемы CIP в случае больших оптических толщин

Начнём с тестирования схемы CIP_S для неоднородного уравнения переноса (1) в случае $\sigma \le 1$. Пусть $a=1$, $\sigma =0,5$, $X=1$, $T=1$.

Табл. 12 демонстрирует третий порядок сходимости схемы CIP_S на гладком тесте при ${\kappa }_0=1$, причём величины ошибок не отличаются существенно от представленных в табл. 8 для схемы CIP_R.