Автоматика и телемеханика, Л- 5, 2007
РАС Б 47.55.pl). 47.11.-]. 47.63.mf, 02.30.Zz
© 2007 г. А.И. КОРОТКИЙ, д-р физ.-мат. наук, И.А. ЦЕПЕЛЕВ, канд. физ.-мат. наук (Институт математики и механики УрО РАН, Екатеринбург)
ПРЯМЫЕ И ОБРАТНЫЕ ЗАДАЧИ ДИНАМИКИ ВЫСОКОВЯЗКОЙ ЖИДКОСТИ1
Описываются методы и алгоритмы числеппого моделирования трехмерных термокопвективпых движений неоднородной высоковязкой несжимаемой жидкости в прямом и обратном направлениях времени. В отличие от задачи в прямом времени задача в обратном направлении времени является некорректно поставленной, для ее решения привлекается одип из вариантов метода квази-обращепия. Основное внимание уделяется организации эффективных вычислительных процессов па ЭВМ параллельного действия. В качестве примера приводятся результаты моделирования обратной задачи о реконструкции истории развития тепловых мантийных плюмов.
1. Введение
В работе описываются методы и алгоритмы численного моделирования термоконвективных движений неоднородной высоковязкой несжимаемой жидкости. Математическая модель динамики (тепловой конвекции) такой жидкости опирается на общие соотношения механики сплошной среды [1 3]: уравнение сохранения количества движения, уравнение сохранения теплового баланса, уравнения состояния, реологические уравнения. В уравнении движения опускаются инерционные члены (поскольку в силу высокой вязкости для характерных движений изучаемой жидкости числа Рейнольдса и Фруда очень малы). Для получившейся системы дифференциальных уравнений ставятся соответствующие краевые задачи для моделирования движения в прямом и обратном направлениях времени. Краевая задача, поставленная в прямом направлении времени, позволяет составить прогноз развития процесса в будущем. Краевая задача, поставленная в обратном направлении времени, позволяет найти предысторию процесса в прошлом. Задача в обратном направлении времени является, как правило, некорректно поставленной и для построения численных методов решения требует привлечения методов регуляризации [4 7].
Аналитическое решение рассматриваемых краевых задач представляет собой довольно сложную проблему и осуществимо в редких случаях. Альтернативный путь решения рассматриваемой проблемы вычислительный эксперимент и численное моделирование. Численное моделирование подобных задач, равно как и большинства задач математической физики, сопряжено с большими объемами вычислений н привлечением больших вычислительных ресурсов. Для получения качественных
1 Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (гранты № 04-07-90120, № 05-01-00098).
результатов необходимо проводить расчеты на достаточно мелких сетках, что и приводит к некоторым дополнительным и повышенным требованиям к вычислительному устройству. Такие вычисления имеет смысл проводить на многопроцессорных вычислительных системах. Они являются мощным и эффективным инструментом в решении подобных задач. Однако для их использования необходимо разработать соответствующие методы и алгоритмы, которые позволят разбить исходную задачу на несколько независимых частей (подзадач) и решать их параллельно [8 11]. О том. как это можно и удобно сделать для рассматриваемых здесь задач, подробно пойдет речь ниже. Предварительно отметим следующее.
При численных расчетах решений краевых задач будут использоваться метод конечных элементов и метод Галеркииа [12. 13] для решения уравнений движения, а также разностные методы и схемы [13 16] для решения уравнения теплового баланса. Алгоритмы численных расчетов будут ориентированы на применение параллельных вычислительных устройств. Решение обратной задачи будет опираться на один из вариантов метода квазиобращения [7. 16]. В конце работы приводятся результаты численных расчетов по моделированию прямой и обратной ретроспективной задач динамики тепловых плюмов в мантии Земли [17 29].
2. Постановки задач
В параллелепипеде О = (0,^) х (0,12) х (0,/з) С М3 рассматривается движение высоковязкой неоднородной несжимаемой теплопроводной жидкости, находящейся под воздействием гравитационного и теплового полей. В декартовых координатах в приближении Буссинеска движение такой жидкости описывается [1 3]: краевой задачей для определения поля скоростей (она включает в себя уравнение Стокса, уравнение несжимаемости, а также соответствующие граничные условия)
Ур = д\у(^(Т) еу ) + КаТ е3, х е О,
. ёгуи = 0, х е О,
^ ^ и = 0, X е гь
(и, п) = 0, дит/дп = 0, х е Г2,
и краевой задачей для определения поля температур (она включает в себя уравнение теплового баланса с соответствующими граничными условиями)
I дТ/дь + (УТ, и) = АТ + /, X е О, ь е [0,0], 1 ' \ о1 Т + о2 дТ/дп = Ть, X е дО, Ь е [0,0].
Здесь и = (и1,и2,и3) - вектор скорости движения жидкости; е3 - орт вертикальной оси Охз; п - единичный вектор внешней нормали в точках границы дО области О; Р ~ давление; ^ - вязкое ть; еу = дщ/дх^ + ди^/дх^ - тензор скоростей деформаций; У - операция взятия градиента; - операция взятия дивергенции; А - оператор Лапласа; Т - температура; / - плотность внутренних источников тепла; Ка - число Рэлея; 01 и 02 — некоторые кусочно-гладкие функции, не об-
Ть
ловой режим па границе области О. Подбирая функции 01, о2 и Ть соответствующим образом, можно задать желаемое распределение температуры или заданный тепловой поток на тех или иных участках границы. Для поля скоростей рассматриваются два типа граничных условий: на Г1 С дО условие прилипания, па Г2 С дО условие идеального скольжения с непроницаемостью, Г1 и Г2 = дО и Г1 П Г2 = 0. Плотность жидкости р = р(Ь, х, Т) связана с температурой уравнением состояния р(Ь, х,Т) = 1 — а (Т(Ь, х) — 1), где а - коэффициент объемного расширения (входит в
число Рэлея). Вязкость ц = ц(Ь,х,Т) связана с температурой реологическим уравнением /л(г,х,Т) = ехр(д/(Т + С) — д/(0,5 + С)), где Q = 255/ 1п(Д) - 0,25/ 1п(Д), С = 15/ 1п(Д) — 0,5 Д - некоторая положительная константа (отношение вязкости на верхнем слое к вязкости на нижнем слое параллелепипеда П). Движение жидкости рассматривается па конечном отрезке времени 0 ^ £ < $.
Для задачи, решаемой в прямом направлении времени, считается известной температура жидкости То = То(х), х € П, в начальный момент времени £ = 0:
(3) Т(0, х) = То(х), х € П.
Для задачи, решаемой в обратном направлении времени, считается известной температура жидкости Т$ = Т$(х), х € П, в финальный момент времени £ = $:
(4) Т($,х) = Т#(х), х € П.
В данной модели не учитывается генерация тепла в жидкости из-за внутреннего вязкого трепня (считается, что такое тепловыделение мало). Реология среды считается ньютоновской. Физические параметры жидкости (температура, скорость, давление, вязкость, плотность) считаются переменными, зависящими от времени и пространственных координат. Вязкость и плотность зависят от температуры. Соотношения (1) и (2) приведены в безразмерном виде. Переход к безразмерному виду был подробно описан в [21].
Зада ча 1. Прямая задача тепловой конвекции: требуется определить температуру Т = Т(£, х), поле скороетей и = и(£, х) и давление р = р(Ь, х), удовлетворяющие при £ € (0,$] и х € П краевым задачам (1)-(2) и начальному условию (3).
За да ча 2. Обратная ретроспективная задача тепловой конвекции: требуется определить температуру Т = Т(Ь,х), поле скоростей и = и(£, х) и давление р = р(Ь,х), удовлетворяющие при £ € [0,$) и х € П краевым задачам (1)-(2) и финальному условию (4).
Аналитические решения подобных задач удается построить только в наиболее простых частных случаях. Поэтому в данной работе все внимание будет уделено построению устойчивых численных методов решения поставленных задач. Задача тепловой конвекции в прямом направлении времени устойчива в том смысле, что малые ошибки при задании физических величин в допустимых исходных данных н погрешности нормального вычислительного процесса не приводят к неограниченному росту погрешности при численном решении задачи и определении состояния среды па практически любом конечном промежутке времени (см.. например. [2. 13 15]). Задача тепловой конвекции в обратном направлении времени является некорректной в том смысле, что малые ошибки при задании физических величин в допустимых исходных данных и погрешности нормального вычислительного процесса могут привести к неограниченному росту погрешности при численном решении задачи и определении состояния среды на любом конечном промежутке времени (см.. например. [7, 10]). Для построения численных алгоритмов решения такой задачи необходимо разработать соответствующие методы регуляризации. Опишем теперь методы и алгоритмы решения поставленных задач.
3. Методы и алгоритмы численных расчетов
Введем предварительно некоторые понятия и обозначения, которые будут ис-
П
моугольную сетку ш с узлами
ш^к = (х^х^х^), 0 ^ г ^ и\, 0 ^ ] ^ п^, 0 ^ к ^ пз, 0 = х0 <х\ < ... < хп-1 < хп = ¡,, хг+1 — хг = Ни 0 < г < щ — 1, г = 1, 2, 3, Ы = и/щ, N =(щ + 1)(п2 + 1)(пз + 1), щ е М, г = 1, 2, 3.
На отрезке времени [0, 0] будем рассматривать разбиени е точками ¿п, п = 0,...,т,
0 = ¿о <tl < ... < гт-1 <гт = 0, т е N.
Символом г^к будем обозначать значение соответствующей сет очной функции г в узле сетки ш^к е ш в момент времени ¿п е [0,0]. Если функция г не зависит от времени или значение момента времени для нее по каким-то причинам не важно, то будем писать г^к- Когда для функции г важно значение момента времени и неважно значение в узле пространственной сетки, то будем писать гп. Сетка ш будет фиксированной, разбиение отрезка [0, 0] будет выполняться адаптивно, т.е. по ходу расчетов в зависимости от сложившихся значений рассчитываемых параметров.
Для внутренних узлов ш^к сетки ш та примере направления х1 определим одномерные сеточные операторы и выражения
Лх1х1 Т = (Тг+10к — 2Т^к + Тг-10к) /Н2 к д2Т/дх\,
Лх1х1х1х1 Т = (Т-2зк — 4Т-1зк + 6Тг0к — + Тг+2]к)/Н\ к д4Т/дх"4, Со(Т, и) = ^-4Н-+-4Н-^
С1(Т, и) = — Р-) /Н1,
Р+1 =0,5и1^к (Т^^к + ТгГ+1/2зк) — ^Кл^ (Т,1+1/2^к — ТгГ+1/2зк) , Р-1 = ®,5и1,1]к (Тг+-1/2зк + Тг-1/2зк) — (Тг+-1/2зк — Тг-1/2зк) ,
Т-+1/2зк = Тцк + 0,5Ф(,), Т+1/ш = Т^к — 0
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.