CC BY
10
Механика
УДК 517.968+550.311 DOI: 10.14529/mmph180104
ВОССТАНОВЛЕНИЕ ПАРАМЕТРОВ ТЕЧЕНИЯ ВЯЗКОЙ ТЕПЛОПРОВОДНОЙ ЖИДКОСТИ ПО НЕКОТОРЫМ ИЗМЕРЕНИЯМ НА ЕЕ ПОВЕРХНОСТИ
А.И. Короткий, И.А. Цепелев
Институт математики и механики им. Н.Н. Красовского УрО РАН, г. Екатеринбург, Российская Федерация
E-mails: korotkii@imm.uran.ru, tsepeiev@imm.uran.ru
Определяются физические характеристики установившегося течения вязкой теплопроводной несжимаемой жидкости по измерениям температуры и потока тепла на ее дневной поверхности. Основными искомыми характеристиками являются температура и скорость жидкости во всей модельной области. Задача формализуется как обратная граничная задача для модели течения естественной тепловой конвекции высоковязкой несжимаемой жидкости. Математическая модель течения такой жидкости описывается стационарными уравнениями Навье-Стокса для ньютоновской реологии среды в приближении Буссинеска в поле силы тяжести, уравнением несжимаемости среды, стационарным уравнением сохранения энергии с соответствующими граничными условиями. Плотность и вязкость жидкости нелинейно зависят от ее температуры. Рассматриваемая обратная задача является некорректной и не обладает свойством устойчивости, малое возмущение исходных данных на участке границы, доступной для измерений, приводит к неконтролируемым ошибкам в определении искомых величин. Для численного решения неустойчивых задач требуется разработка специальных методов. Цель данной работы состоит в построении методов и алгоритмов конструктивного устойчивого численного моделирования решения рассматриваемой обратной задачи. Для реализации этой цели предлагается воспользоваться вариационным методом, который основан на сведении исходной задачи к некоторой экстремальной задаче на минимум подходящего целевого функционала и его устойчивой минимизации каким-либо подходящим способом. При такой стратегии организуется итерационный процесс последовательного численного решения краевых задач граничного управления, которые представляют собой системы дифференциальных уравнений с частными производными с полностью определенными граничными условиями. Для минимизации функционала качества применяется метод сопряженных градиентов в реализации Ролака-Рибьера. Градиент этого функционала и шаг спуска определяются аналитически, что позволяет существенно сократить объем вычислений. Метод конечных объемов применяется для интегрирования систем дифференциальных уравнений с частными производными с различными типами граничных условий. Построенные алгоритмы численного моделирования реализованы в пакете вычислений OpenFOAM. Проведен расчет модельного примера.
Ключевые слова: тепловая конвекция; вязкая жидкость; обратная граничная задача; вариационный метод; численное моделирование.
Введение
Рассматривается задача определения температуры и поля скоростей жидкости, движущейся в некоторой области W с R2, по некоторым заданным условиям на границе области W . Математическая модель стационарной тепловой конвекции вязкой несжимаемой неоднородной жидкости включает в себя стационарное уравнение Навье-Стокса вместе с уравнением несжимаемости для
определения поля скоростей жидкости и стационарное уравнение теплового баланса для определения ее температуры. На границе Г = ЭН пространственной области задаются соответствующие граничные условия для скорости и температуры. При этом некоторые граничные условия задаются только на части границы (условия с недостатком), а на других участках границы задаются дополнительные условия (условия с избытком). Соответствующую краевую задачу для определения поля скоростей и температуры в области Н назовем обратной граничной задачей (или просто обратной задачей).
Рассматриваемая обратная задача является некорректной в том, что и не обладает свойством устойчивости [1, 2], малое возмущение исходных данных на участке границы с переопределением приводит к неконтролируемым ошибкам в реконструкции искомых величин. Для численного решения неустойчивых задач требуется разработка специальных методов и алгоритмов. В данной работе будет представлен конструктивный метод устойчивого численного решения рассматриваемой обратной задачи. Этот метод основан на сведении исходной некорректной задачи к вариационной задаче на минимум целевого функционала среднеквадратичной невязки [1-3]. Для решения экстремальной задачи выбрана одна из реализаций метода сопряженных градиентов. При таком алгоритме решения исходная задача заменяется на последовательное решение аналитически полученных краевых задач граничного управления, для которых соответствующие граничные данные определены полностью.
Обратные задачи восстановления тепловых режимов в геофизических процессах составляют новое направление в моделировании различных технологических и природных процессов (см., например, [4-9]). Необходимо отметить большой интерес к решению подобных задач среди специалистов в области моделирования геодинамических процессов, процессов химического синтеза и многих других отраслей знаний [7-9]. Этому способствует совершенствование методов обработки данных и накопление статистической информации о природном процессе или явлении. Решение обратных задач на сегодняшний день довольно наукоемко и сопряжено с большими затратами ресурсов ЭВМ. Благодаря сочетанию аналитических методов теории управления с алгоритмами численного решения систем дифференциальных уравнений с частными производными, в данной работе удалось получать результаты хорошего качества при сравнительно небольших затратах на компьютерную реализацию. Результаты данного численного моделирования могут дополнить уже известные методы исследования природных процессов распространения потоков вулканической лавы. Некоторые вопросы разрешимости задач граничного управления, необходимые и достаточные условия оптимальности и процедуры их численного решения исследованы
в [7, 8].
Проведены серии вычислительных экспериментов, показавших работоспособность и эффективность предложенного метода. Работа продолжает исследования [9, 11].
Постановка задачи
В области Не М2 с границей Г = ЭН рассматривается установившееся конвективное движение вязкой несжимаемой неоднородной теплопроводной жидкости, находящейся в поле силы тяжести. Математическая модель представлена следующими обезразмерными уравнениями в приближении Буссинеска [4]:
- Яе( и, V) и + У-(^(Уи + Уиг )) = Ур - ЯаТ е2, (1)
V- и = 0, (2)
V- (к уТ)= (и,УТ), (3)
где хеН - точка пространства; и = (их(х);и2(х)) - вектор скорости; Т = Т(х) - температура; р = р(х) - давление; ¡и = М(Т) - вязкость; к(Т) = к(Т)/(ср*) - коэффициент температуропроводности; г* - температурно независимая плотность; к = к(Т) - коэффициент теплопроводности; с - удельная теплоемкость; Яа - число Рэлея [3, с. 7]; Яв - число Рейнольдса [4, с. 87]; е2 = (0,1) - единичный вектор; т - транспонирование тензора; (•,•) - скалярное произведение векторов; V - градиент; V» - дивергенция.
На границе Г = П и Г2 иГэ и Г4 (см. рис. 1) рассматриваются следующие граничные усло-
вия:
Г1 Г 2
Г 3
Г
Т = т,
(V Т, п) = 0, (V Т, п) = 0,
и = 0, ^и, п) = 0,
(V?, п) = 0, (V?, п) = 0, р = 0,
(4)
т = Т4, -к(vT,п) = р, (и,п) = 0, оп -(он,п)п = 0,(V?,п) = 0,
где о°т(уи + VuT) - тензор вязких напряжений. На границе Г1 жидкость втекает с априори
известной температурой Г1? и ее скорость на этой границе считается неизвестной. На границе Г3 она вытекает, напряжения равны нулю, и теплообмен осуществляется только с потоком жидкости. На границе Г4 проводятся измерения температуры Т4 и потока тепла р, для скорости заданы условия непротекания и идеального скольжения. На нижней границе Г2 области заданы условия теплоизолированно-сти для температуры и прилипания для скорости.
Наша конкретная обратная граничная задача состоит в том, чтобы по результатам измерений на Г4 температуры и теплового потока найти температуру и поле скоростей жидкости во всей области й, т. е. требуется найти решение задачи (1)-(4).
Задача (1)-(4) обладает свойством неустойчивости [11]. Если эту задачу решать численно, то малые ошибки в измерениях исходных данных задачи, в частности, ошибки при измерении наблюдаемых величин температуры Т4 и теплового потока р = -к (Т)ЭТ / Эп на Г4 и промежуточные вычислительные погрешности приводят к неконтролируемым вычислительным ошибкам при определении температуры и других искомых величин в исходной задаче. Поэтому для этой задачи разрабатываются специальные устойчивые численные методы.
Обратной задаче (1)-(4) поставим в соответствие вспомогательную задачу граничного управления, которую будем называть прямой задачей. Эта задача в расчетной области описывается системой уравнений (1)-(3) со следующими граничными условиями:
Рис. 1. Расчетная область й
Г
Г 2 Г 3 Г
Т = Т,
(5)
и =Х, (V?, п) = 0,
(V Т, п) = 0, и = 0, (Ур, п) = 0,
(V Т, п) = 0, ^и, п) = 0, р = 0,
4 Т = Т4, (и, п) = 0, оп -(оп, п)п = 0,( V?, п) = 0.
С содержательной стороны данная прямая задача состоит в определении температуры и скорости вязкой жидкости в заданной области определения. Известны все соответствующие граничные данные в модели (5). Также искомыми характеристиками в этой задаче являются давление, вязкость и коэффициент температуропроводности. Вопросы математической корректности данной прямой задачи исследованы в работах [4, 11].
Метод решения обратной задачи
Для построения алгоритма решения задачи перейдём к вариационной постановке обратной задачи [6-9]. Пусть наблюдаемый поток тепла р на Г4 в обратной задаче соответствует некоторому заранее неизвестному распределению скорости (управляющему параметру) и = X на грани*
це Г1 в прямой задаче (1)-(3), (5). Обозначим и - компоненту скорости решения прямой задачи (Т*,и*,Р*) в области й при заданном условии на границе Г1 в (5) таком, что р = -к(Т*)ЭТ* /Эп
на Г4 . Пусть V - некоторое множество допустимых распределений скоростей на Г1, включающий в себя элемент X. Рассмотрим функционал качества
J (X) = J (k (Tx) (vtx, n) -j)2
d Г,
Xe V,
где Tx - компонента решения прямой задачи с условием u = X на Гх в (5). Для искомого решения задачи данный функционал должен принимать нулевое значение при X = X*, J(X*) = 0 . Искомое граничное управление X на Г является минимизирующим элементом в вариационной задаче
J(X) ® min: Xe V . (6)
Таким образом, от решения обратной задачи можно перейти к решению вариационной задачи (6). Для решения данной задачи применим устойчивый итерационный метод сопряженных градиентов Полака-Рибьера [10]. Метод сопряженных градиентов является одним из наиболее эффективных градиентных методов, если рассматривать такие критерии, как точность результатов, затраты компьютерных ресурсов и реализацию его на современных ЭВМ. Последовательность аппроксимаций задачи (6) строится по правилу:
X(n+1) =X(n) + d(n), n =1,2,3,..., (7)
d(») = i -VJ(X^ n =1,
[-VJ(X(n)) + ß(n) d(n-1), n = 2,3,..., b(n) =(vJ(X(n)),VJ(X(n))-vJ(X(n-1)))r ||vJ(X(n-1))|r2, n=2,3,... ,
g(n) = d(n), VJ(X(n)))Г ||VA(X(n))d(n^|Г2 d(n), (8)
VJ(X) = qn - m(TX) (Vw + VwT) | г
VAXC = (VWn) Iг4, где (w,q) - компонента решения (w,z,q) сопряженной задачи
- Re (VuXT w-VwuX) + V-(m(TX) (Vw + VwT )) = V q + z vTX,
V-w = 0,
V(k(TX)Vz) + (ux, Vz) + Ra( w, e2) = ß'(Tx) с граничными условиями:
Vw + VwT, Vux + Vux
+
k(Tx)(vTx, vz); (9)
Г1 U Г2 : Г3:
Г4:
w = 0, (vq,П = 0, (Vw, n) = 0, q = 0, k (Tx) = 0,
z = 2(k(T£)(VTx,n) -j), w = 0, (Vq,n) = 0,
z = 0, z = 0,
(10)
и w - компонента решения (v, w, p) следующей задачи:
- Re ((v, V u x + (ux, V v) + V - (mx <7) + V - (m'(Tx)w(Vux + Vux )) = Vp - Rawe2 V- v = 0,
(v,VTx) + (ux,vW = V -(kxVw) + V (k(Tx)wvTx); с граничными условиями:
Г1: w = 0, v =x, (vp, n) = 0,
Г2 Г3 Г„
(Vw, n) = 0, v = 0, (Vp, n) = 0,
(Vw, n) = 0, (Vv, n) = 0, p = 0,
w = 0, (v,n) = 0, <7n-<n,n)n = 0, (Vp,n) = 0, <7 = Vv + VvT.
(11)
(12)
Г
4
Короткий А.И., Восстановление параметров течения вязкой теплопроводной
Цепелев И.А. жидкости по некоторым измерениям на ее поверхности
Вывод сопряженной задачи основан на линеаризации функций k(Tx+z) = k(TX) + k'ÇT^T + o(|||||) и m(Tx+x) = mT) + m(Tx)T + o(||\), где Tx+zT компоненты решений (Tx+\, ux+\, px+\) и (X u^,p^) прямой задачи (1)-(3), (5) для граничных условий u = X + C и u = X на границе Г1 соответственно (см. также [8]). Для вычисления шага спуска в
(8), нелинейный оператор A^ :£е ¿|(Г1) ® k(Tç)(vT£,п^е L2(r4) для u = X + C на Г1 аппроксимируется выражением A^+Z(X+|) » A^(X) + VA^(X)|.
Таким образом, решение исходной неустойчивой задачи (1)-(4) заменяется на устойчивый итерационный метод сопряженных градиентов решения задачи (7). Прямая задача (1)-(3), (5) и сопряжённая к ней задача (9), (10) позволяют аналитически выписать и достаточно просто вычислить градиент функционала для метода сопряженных градиентов, а задача (11), (12) позволяет вычислить шаг спуска в методе. Основная вычислительная нагрузка в реализации метода лежит на достаточно качественном и быстром решении задач граничного управления. Вопросы разрешимости таких задач (математическая корректность) исследованы в работе [11]. При достаточно большом числе Рэлея (например, большем 106) итерационный метод решения прямой задачи (совместного определения поля скоростей, давления и температуры) в (1)-(3), (5) или аналогично в
(9), (10) и (11), (12) замедляет сходимость и на некоторых исходных данных (например, при скачке вязкости от 1 до 106) демонстрирует расходимость. Этот факт накладывает на метод соответствующие ограничения в его применимости.
Численное моделирование
Для тестирования приведенного алгоритма были разработаны программные коды в пакете OpenFOAM (Open Source Field Operation And Manipulation). Данный пакет - это объектно-ориентированная платформа, в которой на языке программирования С++ реализован широкий набор эффективных процедур аппроксимации дифференциальных операторов с различными типами граничных условий, полученных в моделях механики жидкости и газа, а также алгоритмов решения систем линейных алгебраических уравнений (СЛАУ), которые возникают после дискретизации краевых задач в пространственных областях с различной геометрией. В программные коды легко вносить изменения при изменении исходной математической модели и замене структуры вычислительного кластера. Информация о возможностях пакета представлена на официальном сайте разработчика [12].
Задачи граничного управления дискретизировались методом конечных объемов. Демонстрируемый пример рассчитывался на сетке из 104 полиэдральных ячеек. Для определения поля скоростей и давления при заданном распределении температуры применялся SIMPLE-алгоритм (Semi-Implicit Method for Pressure-Linked Equations) [13]. Тестовые расчеты показывают, что полиэдральная сетка имеет преимущества в качестве перед гексаэдральной (что оказывает существенное влияние на скорость сходимости SIMPLE метода) и преимущества в размерности дискретных задач перед тетраэдральной сеткой (что значительно сокращает время расчетов). Для реализации данного алгоритма решались СЛАУ с положительно определенными и симметричными матрицами с реализацией многосеточного метода GAMG (Geometric Agglomerated Algebraic Multigrid) [14]. После дискретизации уравнения теплового баланса для решения СЛАУ применялся бинаправленный метод сопряженных градиентов [15]. Линейная схема Гаусса с коррекцией потока выбиралась везде для аппроксимации операторов Лапласа для аппроксимации конвективного оператора, который является доминирующим - TVD схема с ограничителем min-mod [16], что обеспечивает монотонность вычисляемых аппроксимаций температуры. Для расчетов использовалось одно ядро CPU Intel Core i5, 2,6 GHz, RAM 16 MB, OS Linux. Тестовый пример рассчитывается примерно 5 мин.
Для проведения вычислительного теста фиксируем следующие характерные значения параметров в модели, соответствующие усредненным значениям в потоке вулканической лавы [8, 9, 17, 18]: a = 10-5 K-1, g = 9,8 м/с2, h = 10 м, pef = 2950 кг/м3, ßeef = 3,5-106 Па с, Tef =
300 K, T* = 1473 K, aT = 100 K, Kef =10-6 м2/с, cref = 1200 J/кг K, uref = 10-3 м/c, вычисляем значения Ra »27000, Rc » 2,75-10-9.
Зависимость вязкости от температуры учитывает кристаллизацию жидкости при понижении температуры
m(T) = min (106, exp(n(T* - T))f(f(T) / f)),
где exp(n(T* - T)) учитывает зависимость вязкости жидкости от температуры, и f(f(T)/ f*) учитывает ее кристаллизацию (аналитическая зависимость рассматривалась в [17]).
Зависимость коэффициента теплопроводности от температуры выражается формулой [18]
k (T) =
11,15 + 5,9 • 10-7( T - ff, T < Г, 11,15 + 9,7 • 10-6( T - T)2, T > T.
Зависимость вязкости и коэффициента теплопроводности от температуры показана на рис. 2. Для проверки качества вычисления градиента функционала в задаче проводился следующий
тест [19]. Для некоторого значения управления X на Г1 (Х^Х*) вычислим вектор
evj(Х)/||VJ(X)||, где ее [10-9,0 ] и рассмотрим функцию
v(e) =
J (X + evJ (Х)/| |VJ (X)||) - J (X)
: 1 + O(e).
(Х)||
На рис. 3 функция п(е) принимает значения близкие к единице в диапазоне параметра ее [10-7,10^ , что демонстрирует качественное вычисление градиента функционала невязки через решение прямой и сопряженной задач.
Рис. 2. Зависимости вязкости (пунктирная линия) и коэффициента теплопроводности (сплошная линия) от температуры
Рис. 3. Тест качества вычисления градиента функционала
Ранее в работе [9] рассматривалась аналогичная задача восстановления теплового поля в модели (1)-(4). Данная методика применялась в задаче восстановления характеристик лавового потока без учета кристаллизации вещества. В качестве управляющего параметра была выбрана температура жидкости на нижней границе области. Было установлено, что в случае доминирования диффузии в модели метод решения из [9] позволяет получить качественное восстановление искомых характеристик. Однако при увеличении магнитуды скорости наведенной конвекции точность вычисления градиента функционала падает (что проявляется в том, что не удается стабилизировать функцию n(e) близко к единице), и в методе сопряженных градиентов не удается гарантировать устойчивую сходимость аппроксимаций и получить качественное решение задачи. Подход, предлагаемый в данной работе, является альтернативным, позволяет проводить расчеты в широком диапазоне величины скорости течений жидкости и дополняет ранее изложенные результаты.
При решении обратной задачи методами сопряженных градиентов было замечено, что невязка на первых итерациях убывает довольно быстро, но затем, когда устойчивые аппроксимации
довольно близко подходят к искомому решению, изменения невязки резко падает. Поведение итерационного процесса, демонстрируемого на рис. 4, является типичным для такого рода задач. Учитывая тот факт, что (8) является решением нелинейной экстремальной задачи
/п) = а^тт { J(Х(п+1)):/п) > 0} , для контроля качества итерационного процесса мы рассматривали неравенство J(Х)) > J(Х(+1)) и останавливали процесс при выполнении равенства J(Х+1))/ J(Х)) — е0 , где е0 выбирается близким к 1, например, 0,99. При этом приемлемое решение задачи достигалось за 20 итераций. В тестовом расчете принималось Х*(х) — 1000е1, Х(0)(х) = 100е1, х еГ1. На рис. 5 изображены погрешности восстановления температуры
Jт —
Т (N) _ т *
¿2( П )
^2 ( Г1)
и скорости JU —
,(N)
¿2( п)
Х
в расчетной области.
Г1)
Заключение
В работе был предложен подход к численному решению некорректной обратной граничной задачи конвективного течения вязкой несжимаемой жидкости, устойчивый к погрешностям вычислений. Алгоритм решения основан на сочетании аналитических методов исследования задач граничного управления и эффективных устойчивых методов решения экстремальных задач. Построен численный метод совместного решения стационарного уравнения Навье-Стокса, уравнения несжимаемости и уравнения теплового баланса с нелинейными коэффициентами, и различными типами граничных условий на неструктурированных сетках. Реализация алгоритма в пакете ОрепБОАМ позволила создать программные коды для решения рассматриваемой задачи, соответствующие современному уровню развития вычислительной техники, в том числе на параллельных вычислительных кластерах. Проведены численные расчеты характерного примера обратной задачи.
Работа выполнена при поддержке Комплексной программы фундаментальных научных исследований УпО РАН (тоект 18-1-1-8).
Номер итерации
Рис. 4. Функционал качества (сплошная линия) и его градиент (пунктирная линия)
Рис. 5. Качество восстановления температуры JT (пунктирная линия) и скорости Ju (сплошная линия)
/
/
Литература
1. Тихонов, А.Н. Методы решения некорректных задач / А.Н. Тихонов, В.Я. Арсенин. -М.: Наука, 1979. - 285 с.
2. Иванов, В.К. Теория линейных некорректных задач и её приложения / В.К. Иванов, В.В. Васин, В.П. Танана. - М.: Наука, 1978. - 206 с.
3. Васильев, Ф.П. Методы оптимизации / Ф.П. Васильев. - М.: Факториал Пресс, 2002. -824 с.
4. Chandrasekhar, S. Hydrodynamic and hydromagnetic stability / S. Chandrasekhar. - Oxford: Clarendon Press, 1961. - 654 p.
5. Ландау, Л.Д. Гидродинамика / Л.Д. Ландау, Е.М. Лифшиц. - М.: Наука, 1986. - 736 с.
6. Короткий, А.И. Реконструкция граничных режимов в обратной задаче тепловой конвекции высоковязкой жидкости / А.И. Короткий, Д.А. Ковтунов // Тр. ИММ УрО РАН. - 2006. - Т. 12, № 2.- С. 88-97.
7. Короткий, А.И. Моделирование прямых и обратных граничных задач для стационарных моделей тепломассопереноса / А.И. Короткий, Ю.В. Стародубцева. - Екатеринбург: Издательство Уральского университета, 2015. - 168 с.
8. Ismail-Zadeh, A. Data-driven numerical modelling in geodynamics: methods and applications / A. Ismail-Zadeh, A. Korotkii, I. Tsepelev. - Berlin: Springer International Publishing, 2016. - 105 p.
9. Quantitative reconstruction of thermal and dynamic characteristics of volcanic lava from surface thermal measurements / A. Korotkii, D. Kovtunov, A. Ismail-Zadeh et al. // Geophysical Journal International. - 2016. - Vol. 205. - Issue 3. - P. 1767-1779.
10. Nocedal, J. Numerical optimization / J. Nocedal, S.J. Wright. - New York: Springer, 1999. -664 p.
11. Короткий, А.И. Оптимальное граничное управление системой, описывающей тепловую конвекцию / А.И. Короткий, Д.А. Ковтунов // Тр. ИММ УрО РАН. - 2010. - Т. 16, № 1. - С. 76101.
12. http://www.openfoam.org/
13. Issa, R.I. Solution of implicitly discretised fluid flow equations by operator-splitting / R.I. Issa // J. Comput. Phys. - 1986. - Vol. 62. - Issue 1. - P. 40-65.
14. Wesseling, P. Geometric multigrid with applications to computational fluid dynamics / P. Wesseling, C.W. Oosterlee // Journal of Computational and Applied Mathematics. - 2001. - Vol. 128. - Issue 1-2. - P. 311-334.
15. Van der Vorst, H.A. BI-CGSTAB: A fast and smoothly converging variant of BI-CG for the solution of nonsymmetric linear systems / H.A. van der Vorst // SIAM Journal on Scientific and Statistical Computing. - 1992. - Vol. 13. - № 2. - P. 631-644.
16. Sweby, P.K. High resolution schemes using flux limiters for hyperbolic conservation laws / P.K. Sweby // SIAM Journal on Numerical Analysis. - 1984. - Vol. 21. - Issue 5. - P. 995-1011.
17. Costa, A. A model for the rheology of particle-bearing suspensions and partially molten tocks/ A. Costa, L. Caricchi, N. Bagdassarov // Geochemistry, Geophysics, Geosystems. - 2009. - Vol. 10, № 3. - Q03010.
18. Hidaka, M. VTFS project: Development of the lava flow simulation code LavaSIM with a model for three-dimensional convection, spreading, and solidification / M. Hidaka, A. Goto, S. Umino, E. Fujita // Geochemistry, Geophysics, Geosystems. - 2005. - Vol. 6, № 7. - Q07008.
19. Navon, I.M. Variational data assimilation with an adiabatic version of the NMC spectral model / I.M. Navon, X. Zou, J. Derber, J. Sela // Monthly Weather Review. - 1992. - Vol. 120, № 7. - P. 14331446.
Поступила в редакцию 3 октября 2017 г.
Bulletin of the South Ural State University Series "Mathematics. Mechanics. Physics" _2018, vol. 10, no. 1, pp. 27-36
DOI: 10.14529/mmph180104
RECOVERY OF FLOW PARAMETERS OF VISCOUS HEAT-CONDUCTING FLUID BY SOME CHANGES AT ITS SURFACE
A.I. Korotkii, ITsepelev
N.N. Krasovskii Institute of Mathematics and Mechanics of the Ural Branch of the Russian Academy of Sciences, Yekaterinburg, Russian Federation E-mails: korotkii@imm.uran.ru, tsepelev@imm.uran.ru
Physical characteristics of steady motion of viscous heat-conducting incompressible fluid by changes of temperature and heat flow on its daylight surface are determined in the article. The main desired characteristics are temperature and fluid velocity in the entire simulation area. The problem is formalized as an inverse boundary problem for the flow model of natural thermal convection of highly viscous incompressible fluid. The mathematical flow model of such fluid is described by a stationary Na-vier-Stokes equations for Newtonian rheology of a medium in Boussinesq approximation in the field of gravity, by the medium incompressibility equation, stationary equation of energy conservation with corresponding boundary conditions. Density and viscosity of the fluid have non-linear dependence on its temperature. The considered inverse problem is incorrect and does not possess the property of stability; small perturbation of initial data on the section of the boundary available for measurement leads to uncontrolled errors in determining the desired values. For numerical solution of unstable problems special methods should be developed. The goal of the article is in developing methods and algorithms of constructive sustainable numerical simulation of the considered inverse problem's solution. In order to implement this goal, the use of variational method, which is based on reduction of the initial problem to some extremum problem on the minimum of the appropriate objective functional and its sustainable minimization by some appropriate technique, is proposed. Using this strategy, an iteration process of sequential numerical solution of boundary problems of boundary control, which are systems of differential equations with partial derivatives with completely determined boundary conditions, is organized. In order to minimize quality functional, the Polac-Ribiere conjugate gradient method is used. This func-tional's gradient and the descent step are determined analytically which allows significantly decreasing the volume of calculations. The method of finite volumes is used for integrating the systems of differential equations with partial derivatives with various types of boundary conditions. The developed algorithms of numerical simulation are implemented in the OpenFOAM calculations package. Calculation of the simulated example is carried out.
Keywords: thermal convection; viscous fluid; inverse boundary problem; variational method; numerical simulation.
References
1. Tikhonov A.N., Arsenin V.Y. On the solution of ill-posed problems. New York: John Wiley and Sons, 1977, 258 p.
2. Ivanov V.K., Vasin V.V., Tanana V.P. Teoriya lineynykh nekorrektnykh zadach i eye priloz-heniya (Theory of linear ill-posed problems and its applications). Moscow, Nauka Publ., 1978, 206 p. (in Russ.).
3. Vasiliev F.P. Metody optimizatsii (Optimization methods). Moscow, Faktorial Publ., 2002, 824 p. (in Russ.).
4. Chandrasekhar S. Hydrodynamic and hydromagnetic stability. Oxford, Clarendon Press, 1961, 654 p.
5. Landau L.D., Lifshitz E.M. Fluid Mechanics. Pergamon, 1987, 554 p.
6. Korotkii A.I., Kovtunov D.A. Rekonstruktsiya granichnykh rezhimov v obratnoy zadache teplo-voy konvektsii vysokovyazkoy zhidkosti (Reconstruction of boundary regimes in an inverse problem of
thermal convection of a high viscous fluid). Trudy Instituta Matematiki i Mehaniki UrO RAN, 2006, Vol. 12, no. 2, pp. 88-97. (in Russ.).
7. Korotkii A.I., Starodubtseva Yu.V. Modelirovanie pryamykh i obratnykh granichnykh zadach dlya statsionarnykh modeley teplomassoperenosa (Modelling of direct and inverse problems for models of stationary heat and mass transfer). Yekaterinburg, Izdatelstvo Uralskogo universiteta Publ., 2015, 168 p. (in Russ.).
8. Ismail-Zadeh A., Korotkii A., Tsepelev I. Data-driven numerical modelling in geodynamics: methods and applications. Berlin: Springer International Publishing, 2016, 105 p. DOI:10.1007/978-3-319-27801-8
9. Korotkii A., Kovtunov D., Ismail-Zadeh A., Tsepelev I., Melnik O. Quantitative reconstruction of thermal and dynamic characteristics of volcanic lava from surface thermal measurements. Geophysical Journal International, 2016, Vol. 205, Issue 3, pp. 1767-1779. D0I:10.1093/gji/ggw117
10. Nocedal J., Wright S.J. Numerical optimization. New York, Springer, 1999, 664 p.
11. Korotkii A.I., Kovtunov D.A. Optimal boundary control of a system describing thermal convection. Proceedings of the Steklov Institute of Mathematics (Supplementary issues), 2011, 272, suppl. 1, pp.74-100.
12. http://www.openfoam.org/
13. Issa R.I. Solution of implicitly discretised fluid flow equations by operator-splitting. Journal of Computational Physics, 1986, Vol. 62, Issue 1, pp. 40-65. DOI: 10.1016/0021-9991(86)90099-9
14. Wesseling P., Oosterlee C.W. Geometric multigrid with applications to computational fluid dynamics. Journal of Computational and Applied Mathematics, 2001, Vol. 128, Issue 1-2, pp. 311-334. DOI: 10.1016/S0377-0427(00)00517-3
15. Van der Vorst H.A. BI-CGSTAB: A fast and smoothly converging variant of BI-CG for the solution of nonsymmetric linear systems. SIAM Journal on Scientific and Statistical Computing, 1992, Vol. 13, no. 2, pp. 631-644. DOI: 10.1137/0913035
16. Sweby P.K. High resolution schemes using flux limiters for hyperbolic conservation laws. SIAM Journal on Numerical Analysis, 1984, Vol. 21, Issue 5, pp. 995-1011. DOI: 10.1137/0721062
17. Costa A., Caricchi L., Bagdassarov N. A model for the rheology of particle-bearing suspensions and partially molten tocks. Geochemistry, Geophysics, Geosystems, 2009, Vol. 10, no. 3, Q03010. DOI: 10.1029/2008GC002138
18. Hidaka M., Goto A., Umino S., Fujita E. VTFS project: Development of the lava flow simulation code LavaSIM with a model for three-dimensional convection, spreading, and solidification. Geochemistry, Geophysics, Geosystems, 2005, Vol. 6, no. 7, Q07008. DOI:10.1029/2004GC000869
19. Navon I.M., Zou X., Derber J., Sela J. Variational data assimilation with an adiabatic version of the NMC spectral model. Monthly Weather Review, 1992, Vol. 120, no. 7, pp. 1433-1446. DOI: 10.1175/1520-0493(1992)120<1433:VDAWAA>2.0.CO;2
Received October 3, 2017
