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

1. Введение

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

Общим требованием к численным методам, которые было бы целесообразно использовать для решения актуальных научно-практических задач, можно считать высокую точность, которую можно было бы достигнуть на не слишком подробных сетках. Стремясь удовлетворить это требование, обычно предлагают два ключевых подхода: расширение шаблона разностной схемы или включение в список неизвестных не только узловых значений искомой величины, но и, например, значений производной в узлах сетки или интегральных средних по ячейке. Каждый из вариантов находит применение в практических задачах. Явные схемы с широким пространственным шаблоном, например, ENO и WENO схемы [1-3], могут использоваться для решения задач газовой динамики, так как уже исходная система уравнений является нелинейной. В задачах переноса излучения большой интерес представляют неявные схемы с компактным шаблоном, что связано с тем, что фотоны движутся со скоростью света, таким образом, применение явных схем привело бы к существенному ограничению на шаг по времени для выполнения условия устойчивости. Примерами таких разностных схем являются бикомпактные схемы Б.В. Рогова [4-10], схема CIP (Cubic Interpolation Polynomial) [11-17]. Бикомпактные схемы Б.В. Рогова строятся методом прямых, и поэтому интегрирование по времени можно осуществить с выбранным порядком аппроксимации. В методе CIP, наоборот, порядки аппроксимации по времени и пространству связаны, так как для построения разностной схемы используется сеточно-характеристический метод (также известный как интерполяционно-характеристический метод). Построение разностной схемы можно разбить на два этапа, первым из которых является построение интерполянта Эрмита по узловым значениям функции и ее первой пространственной производной на нижнем слое по времени. На втором этапе производится интегрирование уравнения переноса вдоль обратной характеристики от точки ее пересечения с границами ячейки до заданной точки нового временного слоя. Для получения значений пространственных производных на новом временном слое может использоваться дифференциальное следствие исходного уравнения переноса, полученное взятием производной по пространственной переменной, однако этот вариант имеет ряд недостатков. В работе [18] был предложен другой вариант, позволяющий получить значения пространственной производной на новом слое по времени, используя интегральные средние и формулу Эйлера–Маклорена. Данная работа продолжает разработку этого варианта.

К достоинствам CIP схемы можно отнести высокий порядок аппроксимации (третий по времени и пространству), малую диссипацию и экстрамалую дисперсию [19]. В [20] предложен способ монотонизации рассматриваемой модификации CIP схемы. Компактность используемого шаблона упрощает ситуацию вблизи границ и при наличии разрывов коэффициентов уравнения. Свойства предложенной в [18] модификации CIP схемы в применении к уравнению переноса с поглощением не были подробно изучены ранее, и в рамках этой работы исследуется способ организации расчета в ячейках, на границу которых попадает разрыв коэффициента поглощения. Исследуются порядки сходимости разностной схемы на аналитических тестах разной гладкости.

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

Рассмотрим задачу нахождения численного решения нестационарного уравнения переноса в плоской одномерной геометрии вида

где $f(x,t)$ является известной функцией. При $\kappa (x,t)=0$ и $f(x,t)=0$ уравнение (1) примет вид

Дополним уравнение (1) начальными условиями

Рассмотрим случай положительной скорости переноса a>0. Поставим на левой границе классические граничные условия вида

В данной работе будет численно решаться задача (1,3-4) для уравнения переноса с поглощением, а также задача (2-4) для линейного однородного уравнения переноса.

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

Метод построения используемой разностной схемы подробно изложен в [18], приведем краткое описание схемы. Пусть введена равномерная сетка с шагом h по пространству и шагом $\tau$ по времени, ограничимся случаем a>0 и $\sigma =\frac{a\tau }{h}\le 1$. Выпустим из точки $(x_{m+1},t^{n+1})$ характеристику $x-at=const$ назад до пересечения с нижней гранью ячейки (Рис.1). Обозначим координату пересечения обратной характеристики с границей ячейки $x*$.

Для построения интерполянта Эрмита на отрезке $[x_m,x_{m+1}]$ на нижнем слое по времени используются $y_m^n$, $y_{m+1}^n$ – сеточные значения искомой функции и $d_m^n$, $d_{m+1}^n$ – сеточные значения пространственных производных искомой функции, поэтому для построения разностной схемы необходимо указать способ расчета в точке $(x_{m+1},t^{n+1})$ значений искомой функции и ее производной по пространству.

Используя характеристическую форму уравнения переноса (1), можно получить в любой точке характеристики

где $\xi =x-at$, $\eta =x+at$. Численное значение $y*$ получается по формулам интерполяции Эрмита, интеграл вдоль характеристики $\xi =const$, входящий в (5), рассчитывается по формуле Симпсона.

Для получения пространственной производной $d_{m+1}^{n+1}$ вычисляется интеграл вдоль правого ребра ячейки на отрезке от $(x_{m+1},t^n)$ до $(x_{m+1},t^{n+1})$ и используется формула Эйлера–Маклорена. Это позволяет найти производную искомой функции по времени $g_{m+1}^{n+1}$, которая уже может быть пересчитана в пространственную производную $d_{m+1}^{n+1}=\left(f_{m+1}^{n+1}-{(\kappa y)}_{m+1}^{n+1}-g_{m+1}^{n+1}\right)/a$ из уравнения (1).

4. Результаты численного тестирования схемы для однородного уравнения переноса

Приведем результаты тестирования схемы CIP, примененной для решения уравнения (2) с начальными условиями (3) и граничными условиями (4). Если задано начальное распределение $u_0(x)$, то решением задачи (2-4) будет функция $u(x,t)=u_0(x-at)$, если краевое условие (4) согласовано с начальным условием (3), то есть $\varphi (t)=u_0(-at)$. Пусть a=1, X=1, $\sigma =0.5$. Для исследования сходимости схемы использовался следующий набор тестов:

В тестах 2,3,4 используются периодические граничные условия.

Третий порядок сходимости схемы по времени и пространству в нормах С, L₁, L₂ демонстрируют результаты теста 1 с бесконечно дифференцируемой функцией в качестве начального распределения (см. Табл. I и Рис. 1А). При понижении степени гладкости функции порядок сходимости уменьшается. В тесте 2 с функцией с разрывом второй производной порядок сходимости равен 1.73, 2.44 и 2.1 в нормах С, L₁ и L₂ соответственно (см. Табл. I и Рис. 1Б). В тесте 3 с функцией с разрывом первой производной порядок сходимости равен 0.74, 1.56 и 1.21 в нормах С, L₁ и L₂ соответственно (см. Табл. I и Рис. 1В). Для теста 4 с разрывной функцией определение порядка сходимости в норме С не является корректным. В нормах L₁ и L₂ получены порядки сходимости 0.78 и 0.42 соответственно (см. Табл. I и Рис. 1Г).

5. Результаты численного тестирования схемы для неоднородного уравнения переноса при постоянном коэффициенте поглощения $\kappa$ = const

Следующим этапом тестирования схемы CIP является применение схемы для решения неоднородного уравнения переноса с поглощением (1). Пусть коэффициент поглощения не зависит от времени и координаты, то есть $\kappa =const$. Для анализа порядка сходимости использованы функции из тестов 1,2,3,4, функция правой части задается в виде $f(x,t)=\kappa \cdot u_0(x-at)$, $\kappa =1$.

Порядок сходимости схемы CIP по времени и пространству в применении к неоднородному уравнению переноса (1) не отличается существенно от порядка сходимости схемы в применении к однородному уравнению (2) для всех проведенных тестов. Тест 1 демонстрирует третий прядок сходимости по времени и пространству во всех рассматриваемых нормах (см. Табл. II и Рис. 2А), тест 2 демонстрирует порядки сходимости 1.71, 2.43 и 2.08 в нормах С, L₁ и L₂ соответственно, в тесте 3 получены порядки сходимости 0.75, 1.59 и 1.21 в нормах С, L₁ и L₂ соответственно (см. Табл. II и Рис. 2Б,2В). В тесте 4 порядки сходимости в нормах L₁ и L₂ равны 0.79 и 0.41 соответственно (см. Табл. II и Рис. 2Г).

6. Выбор значения коэффициента поглощения $\kappa$ в точке разрыва

При построении схемы подразумевалось, что значение $\kappa$ постоянно в пределах одной ячейки сетки, это допускает использование схемы в случае, когда коэффициент $\kappa$ имеет разрыв на границе ячейки. Пусть коэффициент поглощения $\kappa (x,t)\text{}$ на отрезке $0\le x\le 1$ является кусочно-постоянной функцией вида

В точке $x=0.5$ функция (6) терпит разрыв, в ячейке слева от точки разрыва значение коэффициента поглощения равно ${\kappa }_1$, а в ячейке справа – ${\kappa }_2$. Введем обозначение $\kappa (0.5,t)={\kappa }_c$. Одним из этапов расчета по модифицированной схеме CIP является расчет интегрального среднего на отрезке от $(x_{m+1},t^n)$ до $(x_{m+1},t^{n+1})$ и пересчет производной по времени в производную по пространству с использованием уравнения (1) при $x=0.5$. Это вызывает вопрос, какое значение коэффициента поглощения ${\kappa }_c$ использовать при $x=0.5$.

Были рассмотрены различные варианты и проведен численный эксперимент с использованием теста 1 с бесконечно дифференцируемой функцией в качестве начального профиля решения. Функция правой части задавалась в виде $f(x,t)=\kappa (x,t)\cdot u_0(x-at)$, так как коэффициент поглощения выбран в виде разрывной функции, то и функция правой части тоже является разрывной. Расчет проводился при a=1, $\sigma =0.5$, $N_x=100$, $N_t=2N_x$, T=1. Значения коэффициентов ${\kappa }_1$ и ${\kappa }_2$ выбраны величинами одного порядка, чтобы вклад ошибки численного интегрирования по формуле Симпсона в ошибку метода был одинаковым по порядку во всей расчетной области.

Из Табл. III можно видеть, что к наименьшей ошибке в норме С приводит использование при $x=0.5$ коэффициента поглощения из ячейки слева, то есть ${\kappa }_c={\kappa }_1$, причем на результат не влияет то, какое из значений справа или слева от точки разрыва больше по величине. Исследовалась ошибка вычисления производной в точке разрыва $x=0.5$, так как ее величина меньше нормы погрешности вычисления производной $\left|\right|{[{u^{\prime }}_x]}_{\tau h}-d|{|}_C$ во всей расчетной области. Отличие вариантов из Табл. III также невозможно обнаружить по погрешности вычисления самой функции $\left|\right|{[u]}_{\tau h}-y|{|}_C$. Так как в данном тесте было выбрано a>0, то в общем случае оптимальным является выбор ${\kappa }_c$ равным коэффициенту поглощения в ячейке, через которую проходит характеристика, выпущенная из узла с неизвестным значением.

7. Результаты численного тестирования схемы для неоднородного уравнения переноса при кусочно-постоянном коэффициенте поглощения $\kappa$

Применим схему CIP для решения неоднородного уравнения переноса с поглощением (1), коэффициент поглощения зададим в виде разрывной функции (6). Функцию правой части зададим в виде $f(x,t)=\kappa (x,t)\cdot u_0(x-at)$.

Результаты применения схемы в тестах 1,2,3,4 при ${\kappa }_1=1$, ${\kappa }_2=2$ представлены в Табл. IV и на Рис. 3. Порядок сходимости схемы CIP по времени и пространству в применении к неоднородному уравнению переноса (1) с разрывом коэффициента поглощения также не отличается существенно от порядка сходимости схемы в применении к однородному уравнению (2) для всех проведенных тестов. Тест 1 демонстрирует третий прядок сходимости по времени и пространству во всех рассматриваемых нормах (см. Табл. IV и Рис. 3А), тест 2 демонстрирует порядки сходимости 1.72, 2.43 и 2.09 в нормах С, L₁ и L₂ соответственно, в тесте 3 получены порядки сходимости 0.76, 1.62 и 1.22 в нормах С, L₁ и L₂ соответственно (см. Табл. IV и Рис. 3Б,3В). В тесте 4 порядок сходимости в норме С равен нулю, порядки сходимости в нормах L₁ и L₂ равны 0.79 и 0.41 соответственно (см. Табл. IV и Рис. 3Г).

При больших значениях коэффициента поглощения вычисление интеграла, входящего в (5), по формуле Симпсона становится существенно менее точным, погрешность схемы определяется ошибками численного интегрирования. В тесте 1 при ${\kappa }_1=1$, ${\kappa }_2=100$ величина ошибок увеличилась по сравнению с тестом 1 при ${\kappa }_1=1$, ${\kappa }_2=2$ (см. Табл. IV и Рис. 3А), порядок сходимости равен 4, что соответствует порядку аппроксимации формулы Симпсона численного интегрирования (см. Табл. V, Рис. 4).

Предложим еще один тест для исследования предложенной схемы.

где $f_1=3$, $f_2=2$. Наличие разрывов коэффициента поглощения и функции правой части в тесте 5 соответствует распространенной ситуации при решении задач переноса нейтронов.

Результаты применения схемы в тесте 5 демонстрирует Табл. VI. Порядок сходимости схемы CIP по времени и пространству равен 0.79 в норме C и примерно равен 1 в нормах L₁ и L₂.

8. Заключение

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

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