Неявная схема переменных направлений и промежуточные граничные условия

Для решения многомерных задач, приводящих к уравнениям в частных производных параболического типа, широкое применение получил метод переменных направлений (ADI – alternate directions implicit method) [1]. Несмотря на то, что метод имеет давнюю историю и хорошо описан в учебной литературе, попытки его реализации порой оказываются неверными или неточными [2]. Неточность проявляется тогда, когда при учете граничных условий пренебрегается их задание на промежуточных временных шагах. Это пренебрежение может становиться причиной возникновения неустойчивостей даже в том случае, когда сама используемая схема является безусловно устойчивой по спектральному признаку [3]. Тому, как производить корректный учет промежуточных граничных условий для схемы Писмана - Рекфорда (Peaceman-Rachford), посвящены соответствующие разделы в [4, 5].

Ниже мы рассмотрим учет промежуточных граничных условий для схемы переменных направлений Дугласа - Рекфорда (Douglas-Rachford) [3]:
, (1) 
, (2) 
, (3) 

где - оператор Лапласа, - искомая функция, определенная на пространственно-временной сетке:

(4) где
, (5)
, (6)

i, j, s - номера узлов пространственной сетки вдоль направлений Y, X и Z соответственно. - числа узлов пространственной сетки вдоль направлений X, Y и Z соответственно. - величины пространственных шагов вдоль направлений X, Y и Z соответственно. - величина временного под-шага: шаг по времени равен .

Пусть на гранях прямоугольной области заданы некоторые граничные условия (не будем уточнять какие и какого рода):

, (7)

, (8)

, (9)

, (10)

, (11)

, (12)

Перепишем выражения (1), (2) и (3) в такой форме, чтобы неизвестные , и , соответствующие (1), (2) и (3), находились в левой части:

, (13)

, (14)

, (15)

и выясним, какие величины не являются нам известными явным образом. Для этого рассмотрим области, в которых решаются уравнения (13) – (15).

Очевидно, что поскольку в выражение (13) входят операторы и и рассчитать выражения , , и (сохраняя присущий схеме Дугласа - Рекфорда порядок аппроксимации ) оказывается невозможным, областью определения выражения (13) является (см. Рисунок 1):

. 16)

Рис. 1 – Схематичное изображение области моделирования. Внутренний прямоугольник изображает область, заданную соотношением (16), в которой решается уравнение (13). Красными линиями очерчены контуры плоскостей, в точках которых полученные значения не имеют физического смысла (см. обсуждение ниже).

В правых частях выражений (14) и (15) стоят конечно-разностные операторы, которые относятся лишь к тем пространственным направлениям, вдоль которых решается квази-одномерная задача. Поэтому областью определения для (14) и (15) является вся область моделирования:

. (17)

Таким образом, решая уравнение (13), находим (, единица над означает, что величина найдена из первого уравнения Дугласа-Рекфорда) с индексами из области (16).

Однако здесь возникает вопрос – что следует принимать в качестве граничных условий на гранях и в промежуточный момент времени ? Ведь именно величины при и при необходимы для расчета по формуле (14) – см. Рисунок 2.

Рис. 2 – Схематичное изображение области моделирования. Красными линиями очерчены контуры плоскостей, в точках которых необходимо знать истинные (корректные) значения граничных условий для проведения расчета по формуле (14).

При этом, (здесь очень важный момент!) величины при и при , которые рассчитаны из выражения (13) – не являются истинными граничными условиями в момент времени !

Это следует из того факта, что выражение (13) дает решение , в котором учтена временная эволюция только вдоль оси Y. И эта величина, выше мы ее обозначили через , не содержит в себе корректную информацию о граничных условиях в момент времени .

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

Выразим из выражения (14) :

, (18)

а из выражения (15) :

. (19)

Используя (19) можно найти значения при всех

, (20)

поскольку граничные условия являются известными на всех основных (не промежуточных) временных слоях, в соответствии с (4).

Таким образом, зная граничные значения для ,можно по формуле (18) рассчитать “истинные” граничные значения при и при . А эти граничные значения, в свою очередь, определят значения внутри области моделирования. И далее, используя (15), можно определить для всех внутренних точек, переходя, таким образом, на следующий временной слой (см. (4)). Далее все повторяется.

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

Рассмотрим это подробнее.

Пусть требуется решить уравнение теплопроводности, аппроксимируемое следующей схемой Дугласа-Рекфорда на сетке (4):

 
, (21)
, (22)
, (23)
где - теплоемкость, зависящая от температуры.

К этой схеме можно применить описанную выше процедуру нахождения граничных условий на промежуточных временных слоях. Однако сложная функциональная зависимость теплоемкости от температуры приводят к появлению в решаемой системе уравнений (21) – (23) слагаемых следующего вида:

1. в уравнении (22):

(24)

2. в уравнении (23):

, (25)

где I, J, S - номера граничных узлов вдоль осей Y, Х и Z.

Таким образом, необходимо разрешать выражения, содержащие слагаемые вида (24) – (25) относительно и . Данная задача, в принципе, может быть решена численно при условии задания хорошего начального приближения. Однако решение может быть сопряжено со значительными временными затратами и необходимостью производить анализ полученных корней уравнений (таких корней может быть несколько и необходимо будет осуществлять выбор основываясь на каких-либо критериях).

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

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

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

Автор выражает благодарность профессору T.I. Lakoba (Unversity of Vermont) за полезные обсуждения.

Литература:

1. Ch. Gao , Y. Wang “A general formulation of Peaceman and Rachford ADI method for the N-dimensional heat diffusion equation.” Int. Comm. Heat Mass Transfer, Vol. 23, No. 6, pp. 845 – 854 (1996)
2. W.Y. Yang, W. Cao, T.-S. Chung, J. Morris “Applied numerical methods using MATLAB”. Wiley-Interscience (2005).
3. Peaceman D.W. Fundamentals of Numerical Reservoir Simulation. - Amsterdam : Elsevier, 1977.
4. Strikwerda, J. (2004). Finite difference schemes and partial differential equations. Madison: SIAM.
5. Lakoba, T. (13 04 2012 r.). Lecture notes on Numerical Differential Equations. Vermont, USA: University of Vermont. – internet resource.

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

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


*

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