CC BY
32
УДК 517.968+550.311
DOI: 10.14529/mmph160402
ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ТЕЧЕНИЯ ВЯЗКОЙ ЖИДКОСТИ ПО ТЕПЛОВЫМ ИЗМЕРЕНИЯМ НА ЕЕ ПОВЕРХНОСТИ
А.И. Короткий, И.А. Цепелев
Институт математики и механики им. Н.Н. Красовского УрО РАН, г. Екатеринбург, Российская Федерация E-mail: korotkii@imm.uran.ru
Определяются характеристики течения вязкой теплопроводной жидкости по измерениям температуры и потока тепла на ее дневной поверхности. Искомыми характеристиками являются температура и скорость жидкости. Задача рассматривается в стационарной постановке и формализуется как обратная граничная задача для модели высоковязкой несжимаемой жидкости. Задача является некорректной и решается вариационным методом. Проведены расчеты модельных примеров.
Ключевые слова: вязкая жидкость; обратная граничная задача; вариационный метод; численное моделирование.
Введение
Рассматриваются две взаимно связанные между собой задачи. Одна из этих задач условно называется прямой, другая - обратной. Охарактеризуем содержательные стороны этих задач. Прямая задача состоит в определении характеристик потока вязкой жидкости, протекающей в некоторой известной области. Считается, что математическая модель движения жидкости известна, известны также соответствующие исходные (граничные) данные. Искомыми характеристиками в прямой задаче являются поле температур жидкости в известной области, поле скоростей жидкости и переменная вязкость жидкости. Искомыми характеристиками могут быть также распределения вязкости и плотности в известной области или какие-нибудь другие параметры. Обратная задача состоит в определении характеристик потока жидкости по некоторым «избыточным» исходным данным, заданным на части границы известной области. При этом математическая модель движения жидкости считается известной, а на оставшейся части границы области исходные данные заданы с «недостатком» или вовсе не заданы. Один из методов решения обратной задачи состоит в определении «корректных» граничных данных на участке границы с «недостатком» с помощью модели движения и известным граничным данным на участке границы с «избыточными» данными. Затем, с помощью модели движения и всей совокупности граничных данных определяются в результате решения прямой задачи искомые характеристики потока жидкости во всей известной области. Существуют методы, с помощью которых искомые характеристики могут находиться непосредственно по модели движения и граничным данным на участке границы с «избытком». С физической точки зрения «избыточность» граничных данных на некотором участке границы означает доступность этой части границы, возможность делать на ней дополнительные измерения необходимых величин. «Недостаточность» граничных данных на другом участке границы означает недоступность этой части границы, невозможность делать на ней измерения необходимых величин. Для простоты и определенности будут рассматриваться стационарные (установившиеся) варианты задач.
Рассматриваемая обратная задача является, как правило, некорректной и не обладает свойством устойчивости [1, 2], малое возмущение исходных данных на участке границы с переопределением приводит к неконтролируемым ошибкам в определении искомых величин. Обычные классические численные методы не пригодны для решения неустойчивых задач, поэтому для их численного решения требуется применение или разработка специальных методов и алгоритмов. Одна из основных целей данной работы состоит в построении методов и алгоритмов конструктивного устойчивого численного моделирования решения рассматриваемой обратной задачи. Для решения задачи предлагается воспользоваться вариационным методом [1, 2], который основан на сведении исходной задачи к некоторой экстремальной задаче на минимум подходящего целевого функционала (невязки) и минимизации этого функционала каким-либо подходящим методом (например, модифицированными градиентными методами или подходящими квазиньютоновски-
ми методами [3-6]). При минимизации целевого функционала организуется итерационный процесс, который фактически сводит исходную задачу к серии устойчивых корректно поставленных прямых задач.
Привлечение моделей механики вязкой жидкости для моделирования различных технологических и природных процессов получило широкое распространение, как в России, так и за рубежом (см., например, [7-12]). В частности, такие модели используются для моделирования различных процессов в земной коре и недрах, в вулканологии, металлургии, стекольном деле и многих других областях [10-12]. Совершенствование измерительных приборов и методов обработки данных позволяет получать более детальную информацию о процессе или явлении. Решение обратных задач восстановления характеристик среды по доступным прямому измерению данным с помощью методов решения некорректных задач является перспективным направлением исследований. Интерес к таким задачам возрастает в связи с увеличением производительности вычислительных мощностей современных ЭВМ и расширяющимся кругом приложений. Совместный анализ данных и результатов численного моделирования приводит к более глубокому пониманию природы явлений. Дополнительные мотивации к теме исследований и анализ состояния вопроса описаны, например, в [10-12].
В данном исследовании разработаны алгоритмы, реализующие вариационный метод для рассматриваемой обратной задачи. Проведены серии вычислительных экспериментов, показавших работоспособность предложенных метода и алгоритмов.
Постановка задачи
Математическая модель установившегося движения вязкой жидкости области Ос К с границей Г = Г1 иГ2 иГз и Г4 (рис. 1) в приближении Буссинеска представлена обезразмеренными уравнениями [7, 8]
У-(тСО(Уи + Уит )) = Яв{и,V) и + УР-ЯаТе2, хеО, (1)
V- и = 0, х еО, (2)
У-(Я(Т) УТ) = (и, VТ), х еО, (3)
где хе О, - точка пространства с декартовыми координатами (х1;х2) ; и = (м1(х);м2(х)) - вектор скорости движения жидкости; Р = Р(х) - давление; Т = Т(х) - температура; ¡1 = №(Т) - вязкость; 1 = 1(Т) = к(Т) / (ер*) - коэффициент температуропроводности; к = к(Т) - коэффициент теплопроводности; е - удельная теплоемкость; р* - темпе-ратурно независимая плотность; Яв - число Рей-нольдса [8, с. 87]; Яа - число Рэлея [7, с. 7]; е2 = (0;—1) - единичный вектор; V - градиент; т -
Рис. 1. Область W и ее граница Г
операция транспонирования; ^ •,•/ - скалярное произведение векторов; V* - дивергенция.
Прямая задача состоит в нахождении решения (Т, и, Р) системы уравнений (1)-(3), удовлетворяющего граничным условиям
Г1: Т = Т1, и = и1; Г2 : Т = Т2, и = 0 ;
Р = 0 , ап = 0 ; (4)
(и,п) = 0, ап — (ап,п)п = 0;
где S = m(Vu + VuT) - тензор вязких напряжений; Т1, Т2, Т3, Т4 , и1 - заданные функции.
Обратная задача состоит в нахождении решения (Т,и, Р) системы уравнений (1)-(3), удовлетворяющего граничным условиям
Г3: T = Тъ. Г4 : T = T4; '
Короткий А.И., Численное моделирование течения вязкой жидкости
Цепелев И.А. по тепловым измерениям на ее поверхности
G : T = T1, u = Uj; Г2 : u = 0 ;
Г3 : T = T3, P = 0, sn = 0;
Г4 : T = T4, kЭT / Эn = j, (u,n) = 0, sn -(sn,n)n = 0;
функции T1, T3, T4 , j, Uj заданы (T4 , j - результаты измерений температуры и потока тепла).
Метод решения обратной задачи
Дадим вариационную формулировку обратной задаче [9, 10]. Пусть наблюдаемый поток тепла j на Г4 в обратной задаче соответствует некоторому заранее неизвестному распределению
температуры T = T2 = v* на границе Г2 . Пусть T* - компонента решения (T*,u*,P*) прямой задачи при условии T = T2 = v* на границе Г2 в (4), тогда j = k(T*)ЭT* / Э n на Г4 . Пусть V - некоторое множество допустимых распределений температуры на Г2, содержащее в себе элемент v*. Для v е V рассмотрим функционал качества
J(v) = f (k(Tv) ^ -j)2 d Г, ;■! Э n
1 4
где Tv - компонента решения (Tv,uv,Pv) прямой задачи с условием T = T2 = v на Г2 в (4). Функ-
* * ч
ционал качества должен принимать нулевое значение при v = v , J(v ) = 0 . Искомое граничное
температурное распределение v* является минимизирующим элементом в вариационной (экстремальной) задаче
J(v) ® min : ve V . (5)
Таким образом, от решения обратной задачи можно перейти к решению вариационной (экстремальной) задачи (5).
Для минимизации целевого функционала часто необходим градиент этого функционала. Градиент находится по правилу
Э z I
VJ (v) = k (Tv) — Г2, Э n 1 2
где z - компонента решения (z,w, q) сопряженной задачи
V-(u(Tv)(Vw + VwT )) = Vq + zVTv + Re(VuvTw-Vwuv), xe W;
V-w = 0, x eW;
V-(l(Tv )Vz) + <uv, Vz) + Ra( w, e^ = ju'(Tv )2 [Vw + VwT, Vuv + VuvT ] +X\TV)(VTv, Vz), x eW;
Г1 и Г2 : z = 0, w = 0; Г3 : z = 0, q = 0, sn = 0 (s = h(Tv)(Vw + VwT));
Г4 : z = 2 ^ k (Tv) Эр-jj, (w, n) = 0, sn-(sn, n)n = 0.
Функционал качества минимизируется методом сопряженных градиентов в реализации По-лака-Рибьера [4-6].
Численное моделирование
Для демонстрации работы приведенного алгоритма была разработана программа в пакете инженерных вычислений OpenFOAM (Open Source Field Operation And Manipulation). Он представляет собой открытую объектно ориентированную платформу, реализованную на языке программирования С++, обладает большой функциональностью для решения задач механики жидкости и газа, опускает реализацию приложений на мультипроцессорных кластерах. В частности в нем реализованы процедуры решения дифференциальных уравнений с частными производными второго порядка в произвольной пространственной области.
Реализация в пакете OpenFOAM использует широкий набор эффективных процедур аппроксимации дифференциальных операторов, различных типов граничных условий и решения систем линейных алгебраических уравнений, которые возникают после дискретизации краевых задач в пространственных областях с различной геометрией. Программные коды легко адаптировать к изменениям в исходной математической модели и архитектуре вычислителя. Полную информацию о возможностях данного пакета можно найти в [13].
Прямая и сопряженная задачи дискретизировались методом конечных объемов. Тестовый пример рассчитывался на сетке из 1500 гексаэдральных ячеек. Для определения поля скоростей и давления при заданном распределении температуры применялся SIMPLE-алгоритм (Semi-Implicit Method for Pressure-Linked Equations) [14]. Для реализации данного алгоритма решались системы линейных алгебраических уравнений (СЛАУ) с положительно определенными и симметричными матрицами с реализацией метода сопряженных градиентов с подходящим предобуславливателем. Для уравнения теплового баланса соответствующие СЛАУ решалась бинаправленным методом сопряженных градиентов [15]. Для аппроксимации оператора Лапласа всюду выбиралась линейная схема Гаусса с коррекцией потока, для аппроксимации конвективного оператора - TVD схема с ограничителем minmod [16]. Для расчетов использовалось одно ядро CPU Intel Core i5, 2,6 GHz, RAM 16 MB, OS X 10.10. Одна итерация метода сопряженных градиентов рассчитывалась примерно 15 мин. (решение прямой и сопряженной задач и определение шага спуска в методе сопряженных градиентов).
Уточним модельную область W . Граница области состоит из следующих частей: Г1 есть отрезок прямой линии, соединяющий точки (0; 2,5) и (0; 1,5); Г2 - дуга окружности, соединяющая точки (0; 1,5), (1,5; 0,5) и (3, 0;0); Г3 - есть отрезок прямой линии, соединяющий точки (3,0; 0) и (3,0; 0,5); Г4 - дуга окружности, соединяющая точки (3,0; 0,5), (1,5; 1,2) и (0; 2,5).
Фиксируем следующие характерные значения параметров в модели, соответствующие усредненным значениям в потоке лавы [11, 12, 17]: a = 10-5 K-1, g = 9,8 м/с2, h = 10 м, prf = 3000
кг/м3, mef = 3,5-109 Па с, Tef = 300 K, T* = 1473 K, Tc = 1273 K, AT = T* -Tr
ref' ref -8
1r =10-6 м2/с,
ref
= 1200 1/кг К, игв^ = 10 м/с, откуда получаем значения Яа = 100, Яе = 8,5-10 В модельных расчетах принимается зависимость вязкости от температуры
т(Т) = 1 + 2103 (1 + 1апЬ(—100(Тге/Т — Тс))). Принимается зависимость коэффициента теплопроводности от температуры [16]
ТгеГТ — Т*)2, Тгв/Т
k (T) = <
1,15 + 5,9 -10-7(TrefT - T*)2, TeefT < T*
1,15 + 9,7 -10-6 (TrefT - T*)2, TrefT > T*.
(21)
(22)
lref-
ref
На границе Г1 задается температура Т[(х1, х2) = 5,0 — 0,5(х2 — х^), х2 е[х^, х2, ], и скорость иДх2) = и(х2)п1, где пх = (л/2/2;—>/2/2) и и(х2) есть парабола, проходящая через три точки: и (х£) = 10, и (х21) = 0, и (0,5( х2, + х1)) = 7,25. Температура Т3(х1, х2) = 3,5 — 2(х2 — ),
х2 е [х^,х^], на Г3 и Т4(х1,х2) = 4,5 — 2(хх — х^4)/3, хх е [ х4, х1Е ], на Г4 .
Интерес представляет область повышенной вязкости, так как она «тормозит» течения и тепловой поток. Увеличение вязкости жидкости определяется по пороговым значениям вязкости как подмножество О с О, в точках х е G которого
т(Т(х)) >р£ , где т=Апах —е(ттах —mmin), 0 1.
Р ис. 2. Область повышенной вязкости (черная)
Короткий А.И., Численное моделирование течения вязкой жидкости
Цепелев И.А. по тепловым измерениям на ее поверхности
На рис. 2 показано местоположение области повышенной вязкости, которую необходимо восстановить. На рис. 3 иллюстрируется точность восстановления области повышенной вязкости на шагах итерации п = 0, 5, 10, 30 при е = 0,5 .
1000
Рис. 4. Распределение температуры в расчетной области
На рис. 4 показано распределение температуры в области Н; на рис. 5 - точность
восстановления температуры
T(n) (x) — T* (x) на шагах итерации n = 0, 5, 10, 30.
На рис. 6 показано поле распределения скоростей в области W; на рис. 7 - точность восстановления поля скоростей u(n)(x) — u*(x) на шагах итерации n = 0 и n = 30.
Результаты расчетов показывают, что последовательность приближений метода сопряженных градиентов минимизирует функционал невязки устойчивым образом. Для данной аппроксимации расчетной области относительные точности решения систем линейных алгебраических уравнений в SIMPLE методе и методе решения уравнения теплового баланса составляли 10-2, что
согласуется с точностью аппроксимации дифференциальных операторов на заданной неструктурированной сетке. Для тестового примера приемлемый уровень точности решения обратной задачи в определении температуры и поля скоростей достигается примерно за 30 итераций.
Рис. 6. Поле скоростей в расчетной области
n = 0
n = 30
Рис. 7. Точность восстановления поля скоростей в области W
Короткий А.И., Цепелев И.А.
Численное моделирование течения вязкой жидкости по тепловым измерениям на ее поверхности
Заключение
Разработан алгоритм численного решения неустойчивой обратной граничной задачи механики вязкой жидкости, устойчивый к вычислительным погрешностям. Данный алгоритм основан на сочетании аналитических методов исследования математической модели и эффективных устойчивых методов решения экстремальных задач. Реализация алгоритма в пакете OpenFOAM позволила создать программные коды для решения рассматриваемой задачи, соответствующие современному уровню развития вычислительной техники и программного обеспечения для этой техники.
Работа выполнена при поддержке РФФИ (проект № 14-01-00155) и Комплексной программы фундаментальных научных исследований УрО РАН (проект 15-16-1-10).
Литература
1. Тихонов, А.Н. Методы решения некорректных задач / А.Н. Тихонов, В.Я. Арсенин. - М.: Наука, 1979. - 288 с.
2. Иванов, В.К. Теория линейных некорректных задач и её приложения / В.К. Иванов, В.В. Васин, В.П. Танана. - М.: Наука, 1978. - 206 с.
3. Васильев, Ф.П. Методы оптимизации / Ф.П. Васильев. - М.: Факториал Пресс, 2002. -824 с.
4. Nocedal, J. Numerical optimization / J. Nocedal, S.J. Wright. - New York: Springer, 1999. -664 p.
5. Polak, E. Computational methods in optimization: a unified approach / E. Polak. - New York: Academic Press, 1971. - 329 p.
6. Floudas, Ch.A. Encyclopedia of optimization / Ch.A. Floudas, P.M. Pardalos. - New York: Springer, 2009. - 4626 p.
7. Chandrasekhar, S. Hydrodynamic and hydromagnetic stability / S. Chandrasekhar. - Oxford: Clarendon Press, 1961. - 654 p.
8. Ландау, Л. Д. Гидродинамика / Л.Д. Ландау, Е.М. Лифшиц. - М.: Наука, 1986. - 736 с.
9. Короткий, А.И. Реконструкция граничных режимов в обратной задаче тепловой конвекции высоковязкой жидкости / А.И. Короткий, Д.А. Ковтунов // Тр. ИММ УрО РАН. - 2006. - Т. 12, № 2.- С. 88-97.
10. Короткий, А.И. Моделирование прямых и обратных граничных задач для стационарных моделей тепломассопереноса / А.И. Короткий, Ю.В. Стародубцева. - Екатеринбург: Издательство Уральского университета, 2015. - 168 с.
11. 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. DOI 10.1007/978-3-319-27801-8
12. Korotkii, A. Quantitative reconstruction of thermal and dynamic characteristics of volcanic lava from surface thermal measurements / A. Korotkii, D. Kovtunov, A. Ismail-Zadeh, I. Tsepelev, O. Melnik // Geophysical Journal International. - 2016. - Vol. 205. - Issue 3. - P.1767-1779. DOI: 10.1093/gji/ggw117.
13. http://www.openfoam.org/
14. Issa, R.I. Solution of implicitly discretised fluid flow equations by operator-splitting / R.I. Issa // J. Comput. Phys. - 1986. - Vol. 62. - P. 40-65.
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 J. Sci. Stat. Comp. - 1992. -Vol. 13. - № 2. - P. 631-644.
16. Sweby, P.K. High resolution schemes using flux limiters for hyperbolic conservation laws / P.K. Sweby // J. Numer. Anal. - 1984. - Vol. 21. - P. 995-1011.
17. 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 // Geochem. Geophys. Geosyst. - 2005. - № 6. - Q07008, doi:10.1029/2004GC000869.
Поступила в редакцию 27 июня 2016 г.
Bulletin of the South Ural State University Series "Mathematics. Mechanics. Physics" _2016, vol. 8, no. 4, pp. 17-25
DOI: 10.14529/mmph160402
NUMERICAL SIMULATION OF VISCOUS FLUID FLOW BASED ON THERMAL MEASUREMENTS AT ITS SURFACE
A.I. Korotkiy, ITsepelev
Institute of Mathematics and Mechanics of the Ural Branch of the Russian Academy of Sciences (RAS), Yekaterinburg, Russian Federation E-mail: korotkii@imm. uran.ru
The viscous heat-conducting fluid flow characteristics are determined based on temperature measurements and heat flow on its daylight surface. The desired characteristics are temperature and fluid velocity in the whole model area. The problem is considered in a stationary setting and formalized as an inverse boundary problem for the model of high-viscosity incompressible fluid. A mathematical model of this fluid flow is described with the help of the Navier-Stokes equations for a Newtonian fluid in the Boussinesq approximation in a gravity field, incompressible fluid equation, and equation of the energy conservation with the appropriate boundary conditions. The fluid density and viscosity depend on the temperature. The considered problem is incorrect and does not possess the property of stability. Therefore, a small perturbation of the initial data on the accessible part of the border leads to uncontrolled errors in the determination of the unknown quantities in the model area. Conventional classical numerical methods are not suitable for solving the problem, which is why a variation method is used for its numerical solution, which reduces the solution of the original inverse problem to a series of solutions for stable problems. The Polak-Ribiere conjugate gradient method is used to minimize a merit functional in a variation method. This method steadily solves a corresponding extremal problem. The gradient of merit functional is defined analytically as a sequential solution of the direct and conjugate boundary problems. Direct and conjugate problems are numerically solved by the classical method of finite volumes. Constructed algorithms of numerical simulation are implemented in OpenFOAM software. The calculations of model problems are done.
Keywords: viscous fluid; inverse boundary problem; variation method; numerical simulation.
References
1. Tikhonov A.N., Arsenin V.Ya. Metody resheniya nekorrektnykh zadach (Methods of solution for ill-posed problems). Moscow, Nauka Publ., 1979, 288 p. (in Russ.).
2. Ivanov V.K., Vasin V.V., Tanana V.P. Teoriya lineynykh nekorrektnykh zadach i eye priloz-heniya (The theory of linear ill-posed problems and its applications). Moscow, Nauka Publ., 1978, 206 p. (in Russ.).
3. Vasil'ev, F.P. Metody optimizatsii (Optimization methods). Moscow, Faktorial Press Publ., 2002, 824 p. (in Russ.).
4. Nocedal J., Wright S.J. Numerical optimization. New York, Springer, 1999, 664 p. DOI: 10.1007/b98874
5. Polak E. Computational methods in optimization: a unified approach. New York, Academic Press, 1971, 329 p.
6. Floudas Ch.A., Pardalos P.M. Encyclopedia of optimization. New York, Springer, 2009, 4626 p.
7. Chandrasekhar S. Hydrodynamic and hydromagnetic stability. Oxford, Clarendon Press, 1961. 652 p. DOI: 10.1017/s0022112062210592
8. Landau L.D., Lifshits E.M. Gidrodinamika (Hydrodynamics). Moscow, Nauka, 1986, 736 p. (in Russ).
9. Korotkiy A.I., Kovtunov D.A. Rekonstruktsiya granichnykh rezhimov v obratnoy zadache teplo-voy konvektsii vysokovyazkoy zhidkosti (Reconstruction of boundary regimes in the inverse problem of
Короткий А.И., Цепелев И.А.
Численное моделирование течения вязкой жидкости по тепловым измерениям на ее поверхности
thermal convection of a highly viscous fluid). Tr. IMM UrO RAN, 2006, Vol. 12, no. 2, pp. 88-97. (in Russ).
10. Korotkiy A.I., Starodubtseva Yu.V. Modelirovanie pryamykh i obratnykh granichnykh zadach dlya statsionarnykh modeley teplomassoperenosa (Simulation of direct and inverse boundary problems for stationary models of heat and mass transfer). Ekaterinburg: Izdatel'stvo Ural'skogo universiteta Publ., 2015, 168 p. (in Russ.).
11. 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_1
12. 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. DOI: 10.1093/gji/ggw117.
13. http://www.openfoam.org/
14. Issa R.I. Solution of implicitly discretised fluid flow equations by operator-splitting. J. Comput. Phys, 1986, Vol. 62, pp. 40-65. D0I:10.1016/0021-9991(86)90099-9
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 J. Sci. Stat. Comp., 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. J. Numer. Anal., 1984, Vol. 21, pp. 995-1011. DOI: 10.1137/0721062
17. 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. Geochem. Geophys. Geosyst., 2005, no. 6., Q07008. DOI: 10.1029/2004GC000869
Received June 27, 2016
