Квазилинейное уравнение теплопроводности в 3D и задача Стефана в вечномерзлых грунтах в рамках конечно-разностной схемы переменных направлений

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

Аннотация — Произведено исследование квазилинейного уравнения теплопроводности, в котором теплопроводность и теплоемкость зависят от температуры в трех пространственных измерениях, применительно к задаче о фазовом переходе в вечномёрзлых грунтах. Явно сформулированы условия, при которых конечно-разностная схема переменных направлений Дугласа-Рекфорда сохраняет устойчивость. Произведено сравнение полученных численных решений с известным аналитическим решением задачи Стефана в одномерном пространстве.

Ключевые слова — квазилинейное уравнение теплопроводности, задача Стефана, конечные разности, схема переменных направлений, численная устойчивость.

I. Введение

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

Хотя и доказано, что ADI методы являются эффективными и экономичными в отношении временных затрат, и, в большинстве случаев, безусловно устойчивыми, они обладают некоторыми недостатками:

1) Их конечно-разностные формулировки позволяют рассматривать только прямоугольные пространственные области (из-за требований коммутативности, накладываемых на факторизованные и расщепленные операторы) [7];

2) Применение ADI схем к задачам с изменяющимися во времени граничными условиями Неймана и Робина сталкивается с серьезными проблемами из-за необходимости определения этих граничных условий на промежуточных этапах схемы [8];

3) При применении к решению нелинейных уравнений теплопроводности, операторы, составляющие ADI схему, не коммутируют, что приводит к потере безусловной устойчивости схемы.

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

Второй и третий недостатки представляют серьёзную проблему для успешного применения ADI схемы. Насколько нам известно, в настоящее время в литературе отсутствует анализ устойчивости схемы ADI применительно к нелинейному уравнению теплопроводности в трёхмерной пространственной области. Этот факт явился одним из оснований для настоящей работы.

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

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

В настоящей работе мы обсуждаем применение ADI схемы Дугласа - Рекфорда к решению задачи Стефана в пористых вечномёрзлых грунтах. Структура статьи следующая: следующий раздел содержит формулировку задачи и некоторые предположения, которые будут использоваться при доказательстве численной устойчивости ADI схемы, а в разделе 3 представлено само доказательство. В 4 разделе представлены численные результаты, за которыми следуют выводы.

II. Постановка задачи и математический формализм

A. Квазилинейное уравнение теплопроводности

Уравнение, описывающее теплообмен в системе с фазовым переходом имеет следующий вид:

где неизвестная функция (температура: в каждый фиксированный момент , причём является действительным линейным пространством положительных функций), - теплоёмкость, - теплопроводность. Посредством мы обозначим ограниченное, открытое подмножество Евклидова пространства с границей , причем замыкание будет обозначаться . Следуя классификации приведенной в [10], будем называть уравнение (1) квазилинейным уравнением теплопроводности. Начальные и граничные условия принимаются равными:

Мы рассматриваем случаи, когда и , таким образом, уравнение (1) является равномерно параболическим и имеет единственное решение для задачи (3) - (2) [11].

B. Задача Стефана для вечномёрзлых грунтов

Опираясь на результаты [12], сформулируем модель для задачи Стефана без явного использования условия для определения положения фазового фронта. Такой подход оправдан тем, что в случае моделей, использующих явное определение положения фронта в задаче Стефана в засоленных вечномёрзлых грунтах, появляются переохлаждённые области (т.е., некорректное описание края промерзания) [13] (анализ края промерзания и связанных с ним процессов морозного пучения и миграции влаги к фронту промерзания [9] выходят за рамки настоящей работы).
Таким образом, фазовый переход учитывается путём введения эффективной теплоёмкости , которая включает скрытую удельную теплоту фазового перехода :

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

где и - значения теплопроводности в талой и мёрзлой фазах соответственно.

C. Линейное пространство

Для формулировки конечно-разностной схемы мы предлагаем следующую процедуру дискретизации. Пусть вектор имеет положительные координаты и будет множеством всех точек , где индексы являются целыми числами. Две точки и , принадлежащие называются соседними, если . Точки , каждый сосед которых принадлежит обозначаются . Следуя [14], обозначим через точки , обладающие следующим свойством: по крайней мере 1 сосед точек из принадлежит . Таким образом, полная пространственная сетка , причём является множеством граничных точек (вне ). В дальнейшем мы предполагаем, что является однородной (т.е., - константы) кубической областью.

Пусть и . Далее, временная сетка определяется следующим образом:

Приближённые решения уравнения (1), , определяются на и принимают значения в вещественном конечномерном линейном пространстве , размерность которого равна числу точек в . Функция называется допустимой, если и на . Как было отмечено в [14], значения не известны в . Таким образом, должно быть сделано следующее предположение: существует нуль-последовательность (которая стремится к нулю) пространственных шагов, такая что и всегда принадлежит последовательности . Это предоположение гарантирует, что допустимая функция однозначно задаётся в , исходя из начальных и граничных значений .

Для каждой определим разностные операторы первого порядка:

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

Таким образом, являются значениями на фиктивных промежуточных узлах сетки при .

Из уравнений (8)–(9) получается:

где

и соответствует значению в и определяется как . Введем также в -внутреннее произведение и соответственно индуцированную норму по :

Максимум нормы определяется как:

Мы определили в модифицированное -внутреннее произведение и индуцированную норму посредством

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

D. Конечно-разностная схема переменных направлений

Теперь поставим в соответствие уравнению (1) следующую конечно-разностную схему:

где даются формулой (15), , при . Принимая во внимание уравнения (8)–(11) и (15), уравнение (16) может быть переписано в следующем виде:

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

Схема (17) была описана в [1] и называется схемой Дугласа-Рекфорда. Дуглас в [1] отметил, что схема (17) с , определённым уравнением (11), является локально согласованной схемой второго порядка по пространственным координатам до тех пор, пока граничные условия не понизят точность. Из уравнения (16) можно видеть, что уравнение (1) аппроксимируется с точностью в . Граничные и начальные условия (3)–(2) сформулированы на сетке следующим образом:

где .
Хорошо известно, что схема (16) безусловно устойчива, когда и в уравнении (15) являются константами (т.е., и ) и являются конечно-разностными формами оператора Лапласа. Это напрямую видно из выражений (16) и (14), если заметить, что являются самосопряжёнными положительно-определёнными операторами, . В этом случае можно использовать следующую теорему, доказанную в [15]:

Теорема 1:

Если операторы и в уравнении (16) являются самосопряжёнными, не изменяются во времени и - положительно определённая, тогда условие

необходимо и достаточно для проведения следующей оценки:

где - норма, индуцированная -внутренним произведением (14) в случае когда и - константы, - начальное условие, и - численное решение, полученное на -ой итерации.

Действительно, если и являются константами, операторы определённые выражением (15) коммутируют попарно, и соотношение оператора (19) справедливо для любых малых . Соотношение (20) означает, что схема (16) является численно устойчивой, когда и - константы.

Введение и , зависящих от , делает операторы , определённые (15), нелинейными, т.е. не является верным то, что . Этот факт делает невозможным доказательство численной устойчивости в единственной норме так, как это сделано в уравнении (20). Ниже мы немного изменим определение оператора и сделаем следующие предположения:

Предположение 1:

В ADI схеме (16) операторы действуют на векторы следующим образом:

Это предположение приводит к рассмотрению линейной конечно-разностной схемы типа (16) для уравнения теплопроводности с переменными коэффициентами и , формы которых зависят от температурного поля , вычисленного на предыдущем временном слое. Такое предположение широко применяется при рассмотрении нелинейных задач, зависящих от времени (например, [16]).

Перепишем (21) в следующем эквивалентном виде:

где

Предположение 2:

Функции , и (определённые согласно (4)–(6), (9) и (24)) отображают элементы из в пространство последовательностей , элементы которого удовлетворяют следующим соотношениям:

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

III. Численная устойчивость ADI схемы, применяемой к нелинейным уравнениям теплопроводности

В этом разделе мы используем предположения 1 и 2 и получаем условия, при которых схема (16) с выражением (21) могут быть успешно применены к решению нелинейного уравнения теплопроводности (1). Основная проблема, которую необходимо рассмотреть, следующая: получить некоторое соотношение между энергетическими нормами операторов , которое может гарантировать численную устойчивость (16). Для того чтобы продолжить, примем следующее определение численной устойчивости, приведённое в [15]:

Определение 1:

Пусть будет положительно определённым самосопряжённым оператором, действующим в вещественном Гильбертовом пространстве H. Пусть будет константой, так что для любой из (7), - константа, которая не зависит от . Схема (16) называется устойчивой в Гильбертовом пространстве (со скалярным произведением и нормой, определёнными по (14)), если для любой и для любых совершенно малых и , для решения уравнения имеет место следующее соотношение:

где - это тождественный оператор.

Для того чтобы доказать, что схема (16) устойчива в указанном выше смысле, нам понадобятся некоторые вспомогательные результаты, доказанные в [15]:

Лемма 1:

Пусть - Гильбертово пространство (вещественное или комплексное). Пусть - линейный оператор заданный на ,так что . Далее соотношение

) является достаточным

для его выполнения.

Лемма 2:

Пусть - зависящий от времени (для из(7)), положительно определенный и самосопряжённый оператор. Тогда, соотношение

эквивалентно соотношению

где - константа, которая не зависит от .
Теперь, на основании результатов [15], мы можем сформулировать основную теорему:

Теорема 2:

Пусть операторы (определённые согласно (21)) удовлетворяют следующему условию для любого

где - константа. Тогда соотношение

из (7), ) достаточно для проведения следующей оценки:

где - решения уравнения (16).

Доказательство:

Операторы (определённые по (21)) являются положительно определёнными и самосопряжёнными, и, таким образом, можно определить и . Перепишем уравнение (16) в следующей форме:

где

Теперь, используя Лемму 1 и (34), мы получаем:

Таким образом, для и имеем:

Учитывая Лемму 2 и первое уравнение из (38), получаем:

Отсюда,

И, наконец, возвращаясь к скалярному произведению, получим (35).

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

где

max() и min() берутся для всех
Из (44) - (48) следует, что для выполнения неравенства (34), достаточным условием является:

Или более строго

Из (50) достаточный критерий устойчивости (который связывает и и гарантирует, что (34) имеет место на каждом временном уровне ) может быть непосредственно выражен следующим образом:

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

IV. Численные эксперименты

В этом разделе мы демонстрируем результаты применения ADI схемы (16) к решению задачи Стефана в прямоугольном параллелепипеде. Мы сравниваем полученные численные результаты с известными аналитическими выражениями решений для бесконечных областей [17]. Рассмотрим следующую задачу: одномерная полубесконечная область содержит материал в жидком состоянии при и в твёрдом состоянии при , причём является точкой, которая разделяет жидкую и твёрдую фазы. Температура в точке в момент времени обозначается . Температура фазового перехода считается равной градусов Цельсия. Мы также предполагаем нулевые тепловые потоки на концах и кусочно-постоянное начальное условие:

Поставленная задача, как известно, имеет следующее решение [17]:

где константа получается путём решения уравнения:

где erfc() является дополнительной функцией ошибок, - скрытая теплота фазового перехода, и являются температуропроводностями в твёрдой и жидкой фазах соответственно и являются постоянными во времени и пространстве.
В нашем численном эксперименте мы рассмотрели пространственную область с размерами , метров. Пространственная сетка имеет следующее число узлов: 4 x 4 x 172 узлов. На всех гранях параллелепипеда были введены граничные условия Неймана с нулевыми потоками. Начальное условие было взято равным (52) с метров. Значения теплопроводностей и теплоёмкостей талой и мёрзлой фаз, которые содержатся в (4)–(6) были взяты: ,
, . Численная устойчивость для задачи с этими параметрами была достигнута тогда, когда шаг по времени и пространственный шаг связывались следующим соотношением: . Это дало нам число итераций: 6486.

Рис.1 Зависимость температуры от Z-координаты: сравнение результатов расчетов (синие звезды и зеленые треугольники), полученных со схемой (16), и аналитического решения (53) (красные кружки). Момент времени: 28 дней. Вставка показывает те же кривые при температурах ниже Tph = 0 градусов Цельсия. Красная вертикальная линия внутри параллелепипеда обозначает на множество точек, для которых представлены результаты.
На Рис. 1 показаны численные результаты вдоль линии метров (т.е. вдоль центральной линии параллелепипеда) вместе с соответствующим аналитическим результатом, полученным согласно (53). Можно увидеть как качественное, так и количественное совпадение численных и аналитических результатов, причем точность оказывается в пределах 0,5 градуса по Цельсию. Как и ожидалось, лучшее согласие достигается при больших значениях параметра сглаживания в (5), который отвечает за учёт резкости фазового перехода.
Хорошее согласование между численными и аналитическими результатами также наблюдается для временной эволюции температуры в заданной точке. Таким образом, на рисунке 2 показана зависимость температуры от времени в точке с координатами метров.

Рис. 2. Зависимость температуры от времени t в точке метров: сравнение численных результатов, полученных со схемой (16), (синие звезды и зеленые треугольники) и аналитического решения (53) (красные круги) дней.

V. Выводы

Мы показали, что при выполнении Предположений 1 и 2, равномерная численная устойчивость (определение 1) может быть установлена для схемы Дугласа - Рекфорда (17) (которая аппроксимирует уравнение (1)), причем общий вид достаточного критерия устойчивости дается выражением (51). Используемый подход заключается в сведении квазилинейного уравнения к уравнению с переменными коэффициентами (с учетом ограничений (25) - (27)), формы которых изменяются после каждой итерации. Сравнение полученных результатов численного и аналитического решения (рис. (1) - (2)) предполагает, что ADI метод может быть успешно применен к решению задач Стефана в вечномёрзлых грунтах.

REFERENCES

[1] J. Douglas, “Alternating direction methods for three space variables,” Numerische Mathematik, vol. 4, no. 1, pp. 41–63, Dec. 1962.
[2] N. N. Yanenko, The method of fractional steps;: The solution of problems of mathematical physics in several variables, 1st ed. Springer-Verlag, 1971.
[3] G. I. Marchuk, Splitting and Alternating Direction Methods, ser. Finite Difference Methods. North Holland: Elsevier Science Publishers B.V., 1990, vol. 1, pp. 199–462.
[4] K. J. in ’t Hout and B. D. Welfert, “Unconditional stability of secondorder ADI schemes applied to multi-dimensional diffusion equations with mixed derivative terms,” Applied Numerical Mathematics, vol. 59, no. 3-4, pp. 677–692, Mar. 2009.
[5] P. Angot, J. Keating, and P. D. Minev, “A direction splitting algorithm for incompressible flow in complex geometries,” Computer Methods in Applied Mechanics and Engineering, vol. 217-220, pp. 111–120, 2012.
[6] E. Hansen and A. Ostermann, “Dimension splitting for quasilinear parabolic equations,” IMA Journal of Numerical Analysis, vol. 30, no. 3, pp. 857–869, 2010.
[7] G. Birkhoff and R. S. Varga, “Implicit alternating direction methods,” Trans. Amer. Math. Soc., vol. 92, pp. 13–24, 1959.
[8] B. Bialecki and R. Fernandes, “An alternating-direction implicit orthogonal spline collocation scheme for nonlinear parabolic problems on rectangular polygons,” SIAM Journal on Scientific Computing, vol. 28, no. 3, pp. 1054–1077, 2006.
[9] H. R. Thomas, P. Cleall, Y. C. Li, C. Harris, and M. Kern-Luetschg, “Modelling of cryogenic processes in permafrost and seasonally frozen soils,” G´eotechnique, vol. 59, no. 3, pp. 173–184, 2009.
[10] O. A. Ladyˇzenskaja, V. A. Solonnikov, and N. N. Ural’ceva, Linear and quasi-linear equations of parabolic type, ser. Translations of Mathematical Monographs. Providence, RI: American Mathematical Society, 1968, vol. 23.
[11] S. Kamenomostskaya, “On the Stefan problem.” Mat. Sb., N. Ser., vol. 53(95), no. 4, pp. 489 –514, 1961.
[12] A. A. Samarskii and P. N. Vabishchevich, Mathematical Modelling, Volume 1, Computational Heat Transfer, volume 1 ed. Wiley, 1996.
[13] L. Bronfenbrener and E. Korin, “Two-phase zone formation conditions under freezing of porous media,” Journal of Crystal Growth, vol. 198-199, pp. 89–95, 1999.
[14] M. Lees, “Alternating direction and semi-explicit difference methods for parabolic partial differential equations,” Numerische Mathematik, vol. 3, no. 1, pp. 398–412, 1961.
[15] A. A. Samarskii and A. V. Gulin, Stability of Difference Schemes. Moskow: Nauka, 1973.
[16] J. E. Dendy, “Alternating direction methods for nonlinear Time- Dependent problems,” SIAM J. Numer. Anal., vol. 14, no. 2, pp. 313–326, 1977.
[17] E. Javierre, C. Vuik, F. J. Vermolen, and S. van der Zwaag, “A comparison of numerical models for one-dimensional stefan problems,” Journal of Computational and Applied Mathematics, vol. 192, no. 2, pp. 445–459, 2006.

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

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


*

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