ВЕСТНИК ЮЖНОГО НАУЧНОГО ЦЕНТРА РАН Том 5, № 4, 2009, стр. 23-26
= МАТЕМАТИКА И МЕХАНИКА
УДК 532.546.2.013.4:536.24
ЧИСЛЕННЫЙ МЕТОД ИССЛЕДОВАНИЯ КОНВЕКЦИИ МНОГОКОМПОНЕНТНОЙ ЖИДКОСТИ В ПОРИСТОЙ СРЕДЕ
© 2009 г. А.Д. Немцев1, В.Г. Цибулин2
Для исследования режимов гравитационной конвекции многокомпонентной жидкости развита численная схема смещенных сеток. Разработана среда моделирования, позволяющая проводить расчет и визуализацию трехмерных режимов, продолжать их по параметрам задачи, исследовать переходы от стационарных к нестационарным движениям, вычислять спектр устойчивости стационарных режимов. Представлены результаты расчета конвекции двухкомпонентной жидкости в пористом параллелепипеде с двумя противоположными теплоизолированными боковыми гранями и линейным распределением температуры по высоте для других граней.
Ключевые слова: конвекция, пористая среда, многокомпонентная жидкость, схема смещенных сеток, семейство стационарных режимов.
1. ПОСТАНОВКА ЗАДАЧИ
Математическое моделирование фильтрационной конвекции представляет большой интерес в связи с задачами геофизики и реализацией ряда технологических приложений. Из-за нелинейности уравнений, описывающих конвекцию в пористой среде, необходимо применение численных методов и развитие средств вычислительного эксперимента. Особые требования предъявляются к численным схемам, в частности, они должны наследовать существенные свойства рассматриваемых задач. К таким свойствам относится сильная неединственность решения, обнаруженная в плоской задаче фильтрационной конвекции на основе гипотезы Дарси [1]. В этой задаче от состояния механического равновесия ответвляется однопараметрическое семейство стационарных режимов, вызванное ко-симметрией системы. При определенных условиях в трехмерной задаче фильтрационной конвекции также возможно существование семейства плоских стационарных движений [2, 3]. Целью настоящей работы является разработка метода расчета и численное исследование переходов в трехмерной задаче с учетом многокомпонентности жидкости, насыщающей пористую среду.
Уравнения фильтрационной конвекции многокомпонентной жидкости с учетом массовых сил, источников тепла и примесей на основе модели
1 Южный научный центр Российской академии наук, Ростов-на-Дону, 344006, Ростов-на-Дону, ул. Чехова, 41, апет1Беу@ mail.ru.
2 Южный федеральный университет, Ростов-на-Дону, 344006, Ростов-на-Дону, ул. Большая Садовая, 105, 18уЪиИп@ math.rsu.ru.
Дарси выведены в [4]. При отсутствии массовых сил и внутренних источников тепла уравнения в безразмерных переменных записываются в виде
е dv = -dp - v+/ m r ei, (i)
dt 7=1
V ■ v = 0, (2)
яег
ьг —+v■ ver = кгдеr+v■ k, r =1,2,...,5. (3)
dt
Здесь v(x, y, z, t) - векторное поле скорости, p(x, y, z, t) - давление, i1(x, y, z, t) - температура, отсчитываемая от среднего значения, S - 1 - число примесей, ir(x, y, z, t), r = 2, 3, ..., S - отклонения массовых концентраций от средних значений, (x, y, z) - точка области D, занятой жидкостью, t -время, k (0,0,1) - орт, противоположный направлению силы тяжести, е - коэффициент пористости, Ar, кг - фильтрационные параметры Рэлея и коэффициенты диффузии для температуры и каждой компоненты примеси.
Задача рассматривается в параллелепипеде D = [0, Lx] х [0, Ly] х [0, Lz]. На границе для скорости принимается условие непротекания, для температуры на двух противолежащих боковых гранях d1D = {y = 0} U {y = Ly} ставится условие отсутствия теплового и концентрационных потоков. На остальных гранях d2D = dD\d1D поддерживается равновесное распределение температуры и концентраций примесей:
n ■ v = 0, (x, y, z) e SD, (4)
дyer =0, (x, y, z) e Ö1D, er =0, (x, y, z) e 52D, r = 1, 2, ..., S. (5)
24
А.Д. НЕМЦЕВ, В.Г. ЦИБУЛИН
Начальные условия задаются для скорости V, температуры и компонент примесей 6г:
Цх, у, г, 0) = Vо(х, у, 2), 6Г(х, у, г,0) = {0(х, у, 2), г = 1,2,...,
2. МЕТОД СМЕЩЕННЫХ СЕТОК
Дискретизация краевой задачи (1)-(5) осуществляется методом смещенных сеток [5], для чего по пространственным координатам вводятся основные и вспомогательные сетки:
Ьх
XI = гЬх, г = 0,1,..., Ых + 1, Их
Nx + 1
i + 1/2
xi + xi + 1 2
ih%
L
yj = ~ — +jhy, j = 0,1,..., Ny + 1, hy = n
y j + 1/2 :
yj + yj+1
zk = khz, k = 0,1,..., Nz + 1, hz =
Lz
Nz + 1
6k+1/2
zk+zk+1 2
d 1 6,' + 1/2, j, k d1 6i+1/2, j, k =
d2 6i, j+1/2, k '
6i + 1, j, k - 6i, j, k
hx
6i + 1, j, k + 6i, j, k
2
6i, j + 1, k - 6i, j, k
hy
d2 6i, j + 1/2, k '
d3 6i, j, k+1/2 :
d3 6i, j, k+1/2
6i, j+1, k + 6i, j, k
2
6i, j, k+1 - 6i, j, k
К
6i, j, k+1 + 6i, j, k
(6)
2
Для аппроксимации конвективных членов используется формула [5]
з / з
-(6r, ^ ) i, j, k:
J / dm dm (6 ' П dn V
m=1 y nФm
2
3
3
+ J / dm П dn (00 6rdm Vm)
(7)
т = 1 п ф т
В расчетах применяется метод искусственной сжимаемости: в системе (1)-(3) вместо уравнения неразрывности рассматривается уравнение
д р + с-1 V ■ V = 0. (8)
Аппроксимация уравнений (1), (3), (8) записывается с помощью разностных операторов (6), (7):
[V1 - е-1 (-й 1 р - V!)] г,у+1/2, к+1/2 = 0, [V2 - е-1(-й2р - V2)]+1/2,у, к+1/2 = 0,
Здесь Ых, Ыу и Ы2 - число внутренних узлов основной сетки по соответствующей координате. В узлах (хг, у, гк) определяется значение температуры 61 и концентрации примесей 6г (г, у, к), а давление рассчитывается в узлах (хг-1/2, у—1/2> гк_1/2). На вспомогательных сетках, смещенных относительно точек (хг-1/2, у-_т, гк1/2) на полшага по соответствующим координатам, вычисляются компоненты вектора
1 2 3
скорости V1, V и V3.
Такая дискретизация позволяет автоматически удовлетворить краевым условиям для температуры, примесей и скорости на границе д2В. С помощью законтурных узлов для температуры и скорости V2 реализуются краевые условия второго рода на границе дхВ.
Разностные аналоги дифференциальных операторов первого порядка и операторы усреднения по координатам введены на двухточечных шаблонах:
V3-
f-1 (з p - V 3+ / К d1 d2 6 Г
0, (9)
г+1/2,у+1/2, к
[р) + С-1(й 1 V 1 + й2 V 2 + й3 У 3)] г+1/2,у + 1/2, к+1/2 = 0
[Ьг6к - \гАИ6г - 8182 V3 + 7(6г, V)] г,к = 0.
В результате получается система порядка (4 + 5)ЫхЫуЫ2 обыкновенных дифференциальных уравнений, ее интегрирование проводилось при помощи метода Рунге-Кутта.
3. ЧИСЛЕННЫЕ РЕЗУЛЬТАТЫ
На основе описанной схемы была разработана программа Diff_m в среде Delphi для расчета фильтрационной конвекции многокомпонентной жидкости в параллелепипеде. Для визуализации выводятся трехмерные изображения поля скорости, распределения температуры и примеси в различных сечениях области. При выводе изображений используется открытая графическая библиотека OpenGL. Программа позволяет вычислять спектр стационарных режимов при помощи вызова функций OLE-сервера MATLAB.
Для получения сосуществующих стационарных режимов предусмотрены различные способы зада-
2
ИССЛЕДОВАНИЕ КОНВЕКЦИИ МНОГОКОМПОНЕНТНОЙ ЖИДКОСТИ
25
е1 Lz
Рис. 1. Изолинии температуры и концентрации примеси для трех плоских стационарных режимов Ьу = 0,5
Рис. 2. Изолинии температуры для нескольких трехмерных стационарных режимов Ly = 1
ния начальных условий. В частности, отклонения от средних значений температуры и концентраций можно вводить по формуле
йг,. . м АГ . irx ■ jry ■ krz i (i, j, k)= A sin-sin-sin-,
Lx Ly Lz
r = 1, 2, ..., S, i, j, k = 1, 2, ...,
где Ar - амплитуда отклонения. Для проведения эксперимента по трехмерной устойчивости (неустойчивости) плоских движений реализован пользовательский интерфейс, позволяющий использовать данные программы расчета плоской задачи Дарси [6].
На рисунках 1 и 2 представлены результаты расчета конвекции двухкомпонентной жидкости (S = 2, к1 = 1, к2 = 0,2) в параллелепипеде с Lx = 2, Lz = 1 при фиксированных числах Рэлея А1 = 120, m2 = 10. Вычисления для параллелепипедов умеренной длины показали, что в зависимости от глубины области реализуются два сценария ответвления
конвективных движении от теряющего устойчивость механического равновесия v = 0, 6r = 0. При достаточно малой глубине ответвляется семейство плоских стационарных режимов, а при глубине, сравнимой с высотой, реализуются трехмерные изолированные стационарные режимы.
На рисунке 1 приведены некоторые режимы из семейства плоских (v2 = 0) стационарных режимов. Аналогично случаю фильтрационной конвекции однокомпонентной жидкости [3, 5], в спектре всех режимов присутствует практически нулевое собственное значение (10-7 ~ 10-8). Для всех сечений области D, параллельных плоскости xOz, распределения температуры и концентраций примеси одинаковы.
На рисунке 2 представлены изолинии температуры для трех изолированных режимов. Эти конвективные движения имеют сложную трехмерную структуру, в их спектре нет нулевых собственных значений. Режимы I и II являются устойчивыми, а режим III - метастабильным, в результате длительного установления он переходит в режим I. В силу имеющихся в задаче симметрий существуют конвективные движения, отличающиеся от приведенных направлением вращения жидкости и преобразованным температурным полем.
Работа выполнена при финансовой поддержке АВЦП «Развитие научного потенциала высшей школы» (№ 2.1.1/6095).
СПИСОК ЛИТЕРАТУРЫ
1. Юдович В.И. Косимметрия, вырождение решений операторных уравнений, возникновение фильтрационной конвекции // Мат. заметки. Т. 49. 1991. № 5. С. 142-148.
2. Брацун Д.А., Любимов Д.В., Теплое В.С. Трехмерные конвективные движения в пористом цилиндре конечной длины // Гидродинамика. Пермь, 1998. Вып. 11. С. 58-77.
3. Немцев А.Д., Цибулин В.Г. Численное исследование первого перехода в трехмерной задаче фильтрационной конвекции // Известия РАН. Механика жидкости и газа. 2007. № 4. С. 144-150.
4. Юдович В.И. Косимметрия и конвекция многокомпонентной жидкости в пористой среде // Изв. вузов. Северо-Кавказ. регион. Естеств. науки. Спецвыпуск. Математ. моделир. 2001. С. 174-178.
5. KarasozenB., NemtsevA.D., Tsybulin V.G. Staggered grids discretization in three-dimensional Darcy convection // Comput. Phys. Commun. 2008. Vol. 178. P. 885-893.
6. Кантур О.Ю., Цибулин В.Г. Численное исследование
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.