CC BY
84
УДК 532.5
С. В. Каштанова, Н. Н. Окулова
МОДЕЛИРОВАНИЕ ТЕЧЕНИЯ ВЯЗКОЙ ЖИДКОСТИ В КАВЕРНЕ МЕТОДОМ КОНТРОЛЬНЫХ ОБЪEМОВ ПРИ ИСПОЛЬЗОВАНИИ СТАБИЛИЗИРОВАННОГО МЕТОДА БИСОПРЯЖEННЫХ ГРАДИЕНТОВ
На примере задачи о моделировании течения вязкой несжимаемой жидкости в области квадратной формы (каверне) исследовано применение метода контрольных объeмов (алгоритм SIMPLE) к задачам расчeта течений в областях с простой геометрией. В численной реализации алгоритма SIMPLE применены два алгоритма решения линейной системы — метод Гаусса и стабилизированный метод бисопряжeнных градиентов c aiLU-предобусловливанием. Проведены вычислительные эксперименты по исследованию скорости счeта. В первом вычислительном эксперименте исследовалась зависимость скорости счeта от параметра а, во втором — от способа хранения матрицы (йельский формат по сравнению с хранением полной матрицы).
E-mail: svetlyachok@gmail.com; nokulova@gmail.com
Ключевые слова: вязкая жидкость, каверна, метод контрольных объемов, стабилизированный метод бисопряженных градиентов.
Постановка задачи. Рассматривается плоская задача о моделировании течения вязкой несжимаемой жидкости в области квадратной формы (каверне) с верхней границей, движущейся в горизонтальном направлении с заданной скоростью U0.
В области течения поля скоростей и давления удовлетворяют уравнениям Навье-Стокса и уравнению неразрывности:
dv 1
_ + (v -V)v = -Vp + Re Av (1)
V-v = 0. (2)
З V ^ R UoL Здесь v = (u, v) — скорость жидкости, p — давление, Re =--
число Рейнольдса, v — коэффициент кинематической вязкости, L — характерный размер, в качестве которого выбирается длина стороны каверны. Система уравнений (1), (2) выписана в безразмерных переменных; величины, имеющие размерность длины, отнесены к L, скорости — к U0, давление — к pU02, где р — плотность жидкости.
Поместим начало координат в левый нижний угол области, тогда граничные условия имеют вид:
(u,v) = (1, 0) при y = 1; (u, v) = (0, 0) на остальной части границы;
dp
— = 0 на всей границе.
on
Будем считать, что при t = 0 жидкость в каверне покоится, а давление во всех точках постоянное (p = p0).
Отметим, что в такой постановке давление в области течения определено с точностью до произвольной постоянной, поэтому будем полагать p = 0 в правом верхнем углу каверны.
Метод решения. Первоначально для решения системы (1) и (2) был опробован метод конечных элементов, реализованный согласно [1, 2]. Однако получаемые поля скоростей оказывались нефизичны-ми, что связано, по всей видимости, с методикой аппроксимации конвективного члена в (1), не обеспечивающей выполнения дискретных аналогов законов сохранения.
В дальнейшем использовался интегро-интерполяционный метод, или, по-другому, метод контрольных объемов (МКО). В данном методе расчетная область разбивается на некоторое число непересекающихся подобластей (контрольных объемов) таким образом, что каждая точка для расчета неизвестной переменной (узловая точка) содержится в одном контрольном объеме. Специфика МКО заключается в том, что дискретные аналоги уравнений (1), (2) получаются как следствия выполнения законов сохранения для каждого контрольного объема.
Во избежание трудностей, указанных в [3, 4] (потеря точности решения или его нефизичность), для каждой переменной следует использовать свою сетку. Расчетная область П = {(0,1) х (0,1)} делится на квадратные ячейки П^ = {(xi-1,xxi) х (yj-1,yj)} — контрольные объемы для давления (основная сетка). Узловые точки для расчета давления pi ,j находятся в центрах соответствующих контрольных объемов. Составляющие скорости рассчитываются в точках, лежащих на гранях контрольных объемов для расчета давления (для u — на гранях, перпендикулярных оси Ox, для v — на гранях, перпендикулярных оси Oy ). Соответствующие контрольные объемы для щ,j и vi,j показаны на рис. 1. Построенная сетка носит название шахматной, или MAC-разнесенной.
Алгоритм решения следующий. На каждом шаге по времени At в узлах основной сетки задаются значения давления p0 j и скоростей Щ,j, v0,j, известные с предыдущего временного слоя или из начальных и граничных условий (u0j и v0j можно вычислить, например, как среднеарифметическое соответствующих компонентов скоростей, найденных в узлах шахматной сетки). Далее выполняется итерационная процедура вычисления давления pi,j и скоростей ui,j, vi,j в момент времени t = t0 + At на шахматной сетке — алгоритм SIMPLE (Semi-Implicit Method for Pressure-Linked Equations) [4].
Рис. 1. Контрольный объем для расчета давления (а), скорости по осям Ох (б)
и Оу (в)
1. Для поля давления р*,3- (найденного на предыдущей итерации) находятся компоненты скоростей и*3, ь*3 на шахматной сетке:
аи,3 и*,3 = ,3 иг-1,3 + ,3 и*+1,3 + аг,3-1иг ,3-1 + аг,3+1иг,3+1 +
+а0,3и0,з + (р*,з — Р*+1з) Ду;
а,'-,3 Ь*,3 = ^«-1,3 Ьг—1,3 + аг+1,3 Ьг+1,3 + аг,3-1Ьг,3-1 + аг,3+1Ьг,3+1 +
„ + (р*,3 — р*,3+0 Дх'
Данные уравнения являются дискретными аналогами уравнений (1), (2), обеспечивающими сохранение массы, импульса и кинетической энергии в каждой ячейке. Значения коэффициентов аи и аУ в правых частях уравнений связаны с влиянием совместных конвективных и диффузионных процессов на гранях контрольных объемов для и и ь соответственно [4]. При этом а0,3 = ДхДу/Д£, Дх = х»-1 — х^, ДУ = У 3-1 — У 3.
2. Предполагается, что искомое давление находится из соотношения р^,3 = ар 'р'г 3 + р*,3, где ар > 0 — коэффициент релаксации. Этот коэффициент подбирается в зависимости от конкретной задачи и влияет на сходимость метода. Поправка давления р^,3 находится из дискретного аналога уравнения неразрывности (2):
р / _ р / i р ' _|_р ' -¡- р ' i
аг,3рг,3 = аг-1,3рг-1,3 + аг+1,3рг+1,3 + а3-1,гр-1,г + а3+1,гр+1,г + Ь,
где коэффициенты ар определяются в результате подстановки поправочных формул для скорости в проинтегрированное по основному контрольному объему уравнение неразрывности, свободный член Ь = (и*,3 — и*-13)Ду — (ь*3 — ь*3-1)Дх. Отметим, что учет граничных условий происходит через приведенное уравнение для поправок давления [4].
3. Рассчитывается суммарное поле давления p,hj = ap p'j + p'j.
4. Рассчитывается суммарное поле скоростей U' j и v' j:
u',j = u*,j+aP(p',j -p'+i,j )Ay/aU j, v',j = v*j+ap(p'j -p'j+1)Ax/avij.
5. В качестве критерия сходимости итерационного процесса берется малость изменения искомых величин на соседних итерациях: max(|u'j — u'j|, |v',j — v'j |) < S, |p'j | < е. Если критерий не выполняется, то алгоритм повторяется, при этом в качестве p'j берется найденное давление p,hj.
Поскольку интерес представляет установившееся течение, то по окончании работы алгоритма проверяем критерий останова по времени; в качестве такого критерия можно взять малость скорости изменения величин на временном шаге. Если критерий останова не выполнен, то переходим на следующий временной слой, если выполнен — найдено установившееся поле скоростей и давления.
Численная реализация. Задача о моделировании течения вязкой жидкости в каверне является известным тестовым примером. В данном случае на этой задаче было протестировано применение метода контрольных объемов (алгоритма SIMPLE) к задачам расчета течения вязкой несжимаемой жидкости.
Первоначально алгоритм SIMPLE был реализован с помощью пакета Wolfram Ma^mat^. Все матрицы хранились полностью, системы линейных алгебраических уравнений решались с помощью встроенного алгоритма метода Гаусса. Расчеты проводились для различных чисел Рейнольдса (Яе = 100, 1000, 2000), во всех случаях полученные результаты хорошо согласовывались с известными в литературе данными [5]. Однако расчет на сетке 10x10 узлов занял 965 с. Поэтому была поставлена задача увеличения скорости счета.
Вторая численная реализация алгоритма SIMPLE была осуществлена в среде C++ ВшШег. В этом варианте оптимизирован алгоритм решения системы линейных уравнений и формат хранения матриц.
Большинство матриц, получающихся в ходе реализации алгоритма SIMPLE, являются разреженными; одним из эффективных способов хранения таких матриц является (использованный в данной работе) йельский формат (Уа1е format). В нем матрица А хранится при использовании следующих массивов:
1) aelem, содержащий ненулевые элементы матрицы А, перечисленные в строчном порядке;
2) jptr, содержащий столько же элементов, сколько ае1ет, и для каждого из них указывает, в каком столбце находится данный элемент;
3) iptr, который имеет длину, равную размерности матрицы, и каждый элемент которого равен числу ненулевых элементов в строке матрицы.
Системы линейных уравнений решаются с помощью стабилизированного метода бисопряженных градиентов c aiLU-предобусловли-ванием. Матрицы приближенно представляются в виде произведения нижнетреугольной и верхнетреугольной матриц (iLU-факторизация). Такое представление позволяет решать систему линейных алгебраических уравнений Ax = b путем выполнения обратных ходов для системы Ly = b и системы Ux = y. При этом портрет матрицы A (множество позиций матрицы с ненулевыми элементами) включает в себя портреты матриц L и U. Алгоритм iLU-факторизации выглядит так [6]:
for k = 1 . ..n
for j = 1 ...k — 1 1 j—1
if (k, j) Pa then Lkj = 0 else Lkj = u-(Afci — L^U,);
j := j + 1; lkk := 1;
for j = k ...n
k—1
if (k, j) </ Pa then Ukj = 0 else Ukj = Akj — ^ LkiUij;
i=1
j := j +1; k := k + 1.
Поскольку iLU-разложение неединственно, можно допустить модификацию формул для вычисления элементов матриц L и U: элементы матрицы L вычисляются по приведенной формуле, а элементы матрицы U — по формуле 1
Ukj = akj — (1 — a) ^ lkjщ.
i=1
Такое разложение называется aiLU-разложением.
Ниже приведен алгоритм стабилизированного метода бисопряженных градиентов с aiLU-предобусловливанием (BiCGStab+aiLU) [6, 7]:
M = LU, Vx0; r0 = b — Ax0, p0 = r0 Vrc : (r 0,r0) = 0; while ||rn|| > в do
yn = M-1pn; an = ; sn = rn - anAyn; = M-1s
(Ayn,r0):
(Azn sn)
Wn = Д n , ) ; xn+1 = xn + anyn + wnzn; rn+1 = sn — wnAzn; (Azn,Azn)
a (rn+1 r0)
вп = n (, n ); pn+1 = rn+1 + вп (pn — WnAyn). Wn (rn,r°)
Алгоритм не требует обращения и умножения матриц. Необходимо один раз вычислить матрицы L и U, а затем запускать для соответствующих систем алгоритм обратного хода метода Гаусса.
Рис. 2. Диаграмма распределения направления скорости
Результаты расчетов. С помощью разработанных алгоритмов были проведены расчеты полей скоростей и давления при различных значениях числа Рейнольдса, построены диаграммы распределения скоростей по направлению и абсолютной величине и диаграммы распределения давления по каверне. Например, результаты расчетов для значений параметров задачи Яе = 2000, ар = 0.01, а = 0,31, 8 = 10-3, е = 10-4, в = 10-7 (параметр критерия останова по времени), АЬ = 0,01, приведены на рис. 2,3. Видно, что вследствие влияния вязкости поле скорости не является симметричным. Результаты, полученные при различных значениях числа Рейнольдса (Яе = 100, 1000, 2000), качественно не отличаются друг от друга. В нижних углах каверны наблюдаются вторичные вихри. Причем чем число Рейнольдса больше, тем сильнее развитие вторичных вихрей. Наибольшие значения скорости достигаются вблизи подвижной крышки и в окрестности правого верхнего угла. Максимальное значение давления наблюдается в правом верхнем углу каверны.
Далее были проведены два вычислительных эксперимента по исследованию скорости счета. В первом исследовалась зависимость скорости счета от параметра а (при условии, что фиксированы остальные параметры). Во втором — от способа хранения матрицы (йельский формат по сравнению с храненем полной матрицы). В обоих случаях использовался алгоритм, реализованный на языке программирования С++.
Исследования по выбору оптимального значения параметра пред-обусловливания а проводились для задачи с фиксированными пара-
метрами Яе = 2000, 8 = 10-3, е = 10-4, 9 = 10-7. Результаты вычислений на сетке 20x20 приведены в табл. 1. Видно, что изменение а в диапазоне от 0,1 до 0,4 мало влияет на скорость счета, а наилучшим шагом по времени является Д£ = 0,15. Аналогичные эксперименты были проведены на других сетках, они показали, что чем подробнее сетка, тем оптимальные шаг по времени и параметр ар меньше. Увеличение этих параметров ведет к численной неустойчивости.
Таблица 1
Зависимость скорости счета от параметра а (для сетки 20 х 20)
а At ар Время счета, с
0,1 0,05 0,05 113,5
0,1 0,05 0,10 Процесс расходится
0,1 0,10 0,05 84,4
0,1 0,10 0,10 Процесс расходится
0,1 0,15 0,05 79,5
0,1 0,15 0,10 Процесс расходится
0,2 0,05 0,05 110,5
0,2 0,10 0,05 87,4
0,2 0,15 0,05 85,2
0,3 0,05 0,05 110,0
0,3 0,10 0,05 85,7
0,3 0,15 0,05 83,0
0,3 0,20 0,05 Процесс расходится
0,4 0,05 0,05 110,6
0,4 0,10 0,05 92,7
0,4 0,15 0,05 82,1
Во втором вычислительном эксперименте исследовалось, какое ускорение счета при решении всей задачи получается, если матрица хранится в йельском формате по сравнению с хранением полной матрицы. Отметим, что на каждой итерации метода BiCGStab+iLU проводится два умножения матрицы на вектор. Расчеты проводились для задачи с фиксированными параметрами Re = 2000, 8 = 10-3, е = 10-4, в = 10-7. Результаты данного вычислительного эксперимента приведены в табл. 2, из которой следует, что использование йель-ского формата хранения матриц для умножения матрицы на вектор значительно (на 30% и более) ускоряет процесс счета. Отметим также, что в алгоритме, реализованном в пакете Wolfram Mathematica (матрицы хранились полностью, решение линейной системы производилось методом Гаусса), расчет на сетке 10x10 узлов занимает, как было указано ранее, 965 с.
Таблица 2
Вычислительный эксперимент № 2
Время счета, с
Сетка At ap а Без использования йельского формата хранения С использованием йельского формата хранения Полученное ускорение, %
10x10 0,20 0,10 0,3 4,0 2,8 31,5
15x15 0,10 0,10 0,3 50,3 32,2 35,9
20x20 0,15 0,05 0,3 132,4 83,0 37,3
25x25 0,05 0,05 0,3 513,7 363,0 29,3
Перспективы дальнейших исследований. Укажем некоторые ограничения рассмотренного метода. Алгоритм SIMPLE предполагает, что основная сетка полностью покрывает область течения. Это означает, что область со сложной геометрией необходимо заменить областью более простой формы, такой, чтобы ее граница совпадала с гранями контрольных объемов основной сетки. Для того чтобы получаемое на модифицированной области решение было хорошим приближением решения в исходной области, необходимо построение достаточно мелкой сетки.
Еще одной особенностью метода SIMPLE является его фактическая применимость только к стационарным задачам, а при решении не стационарных задач переход с одного временного слоя на следующий требует полного выполнения всего итерационного алгоритма SIMPLE. Поэтому время, необходимое для решения нестационарной задачи на достаточно мелкой сетке, будет чрезвычайно велико.
Одним из наиболее эффективных путей решения возникающих сложностей представляется использование метода, предложенного в [8]. Алгоритм метода LS-STAG основан на MAC-методе для разнесе-ных сеток (как и алгоритм SIMPLE), но расчетная сетка не связывается с границей области. Данный метод относится к классу методов погруженных границ — прямоугольные ячейки, которые пересекаются границей области течения, усекаются в соответствии с границей аппроксимированной кусочно-линейной функцией. Таким образом, область решения представляет собой совокупность прямоугольных ячеек, не имеющих пересечений с границей области течения, и усеченных ячеек.
Дискретные аналоги уравнений Навье-Стокса и неразрывности (в интегральной формулировке) в прямоугольных и усеченных ячейках записываются на одном и том же пятиточечном шаблоне, при этом выполняются дискретные аналоги законов сохранения массы, импульса и кинетической энергии.
Представляется целесообразным распространение методики LS-STAG на решение связанных задач тепломассопереноса. В таком случае система уравнений, которую требуется решить, включает в себя уравнение переноса теплоты
dT
_ + (V -V)T = a2AT, (3)
где T = T— температура, а2 — коэффициент температуропроводности.
В интегральном виде система уравнений (1)-(3) имеет вид:
- / u-V + / (,ñ)„dS + / p, ■ ndS - ±1 Vu ■ ndS = „;
п г г
dt / v-V + / (V. ,,v-S + / „ . n-S - ¿jr Vv . n-S = „;
п г г
J V ■ ndS = „; ^JtóV + J(V ■ n)TdS - a2 J VT ■ n dS = „, г п г г
и для ее решения могут быть применены аналогичные идеи и подходы к аппроксимации уравнений, что и при решении задач гидродинамики методом LS-STAG.
СПИСОК ЛИТЕРАТУРЫ
1. Lewis R. W., Nithiarasu P., Seetharamu K. N. Fundamentals of the finite element method for heat and fluid flow. - Chichester: John Wiley and Sons Ltd., 2004. - 335 p.
2. ZenkiewiczO. C., T a y l o r R. L. The finite element method. Vol. 3: Fluid dynamics. - 2000. - 347 p.
3. ПатанкарС. В. Численное решение задач теплопроводности и конвективного теплообмена при течении в каналах. - М.: Изд-во МЭИ, 2003. - 312 с.
4. П а т а н к а р С. В. Численные методы решения задач теплообмена и динамики жидкостей. - М.: Энергоатомиздат, 1984. - 152 с.
5. Г у р о в Д. Б., Е л и з а р о в а Т. Г., Шеретов Ю. В. Численное моделирование течений жидкости в каверне на основе квазигидродинамической системы уравнений // Математическое моделирование. - 1996. - Т. 8, № 7. - C. 33-44.
6. Баладин М. Ю., Шурина Э. П. Методы решения СЛАУ большой размерности. - Новосибирск: Изд-во НГТУ, 2000. - 70 с.
7. П у з и к о в а В. В. Решение систем линейных алгебраических уравнений методом BiCGStab с предобусловливанием // Вестник МГТУ им. Н.Э. Баумана. Сер. Естественные науки. - 2011. Настоящий выпуск. - С. 124-133.
8. C h e n y Y., Botella O. The LS-STAG method: A new immersed boundary/level-set method for the computation of incompressible viscous flows in complex moving geometries with good conservation properties // J. Comput. Phys. -2010. - No. 229. - P. 1043-1076.
Статья поступила в редакцию 03.10.2011
