CC BY
4
УДК 517.958+532.5.032
DOI: 10.18698/1812-3368-2021-3-46-63
ИСПОЛЬЗОВАНИЕ МОДИФИЦИРОВАННОГО МЕТОДА LS-STAG ДЛЯ РАСЧЕТА ПЛОСКОГО ТЕЧЕНИЯ ВЯЗКОУПРУГОЙ ЖИДКОСТИ В КАНАЛЕ С ВНЕЗАПНЫМ СУЖЕНИЕМ 4:1
И.К. Марчевский1, 2 iliamarchevsky@bmstu.ru
В.В. Пузикова2 vvp@dms-at.ru
1 МГТУ им. Н.Э. Баумана, Москва, Российская Федерация
2 ИСП РАН, Москва, Российская Федерация
Аннотация
Разработана модификация метода погруженных границ ЬБ-БТАС для моделирования течений вязкоупругих жидкостей, описываемых линейными и квазилинейными моделями скоростного типа. Данный численный метод реализован в разрабатываемом авторами программном комплексе для моделирования течений вязкой несжимаемой среды методом ЬБ-БТАС. Программный комплекс позволяет проводить расчеты течений вязкоупругих жидкостей, описываемых моделями Максвелла, Джеффри, Джонсона — Сигельмана, Макс-велла-А, Олдройда-Б, Олдройда-А, верхней конвективной моделью Максвелла. Описано построение дискретных аналогов конвективных производных (Олдройда, Коттера — Ривлина, Яумана — Зарембы — Нолла). К трем разнесенным сеткам базового метода ЬБ-БТАС добавлена четвертая сетка, ячейки которой являются контрольными объемами для построения дискретного аналога уравнения для расчета касательных неньютоновских вязкоупругих напряжений. В результате нормальные вязкоупругие напряжения вычисляются в центрах ячеек основной сетки, а касательные — в их углах. Интегрирование по времени получающейся после ЬБ-БТАС-дискретизации дифференциально-алгебраической системы проводится с помощью метода, основанного на схеме предиктор-корректор первого порядка. Рассмотрена модельная задача о плоском течении вязкоупругой жидкости, описываемой моделью Олдройда-Б, в канале с внезапным сужением 4:1. Определено предельное значение числа Вайсенберга. Результаты вычислительных экспериментов хорошо согласуются с приведенными в литературе данными
Ключевые слова
Метод погруженных границ, метод LS-STAG, вязкоупругая жидкость скоростного типа, несжимаемая среда, канал с внезапным сужением, конвективная производная
Поступила 18.08.2020 Принята 15.09.2020 © Автор(ы), 2021
Работа выполнена при финансовой поддержке РНФ (проект РНФ № 17-79-20445)
Введение. Описание многих промышленных процессов, например литьевого формования, требует рассмотрения течений неньютоновских вяз-коупругих [1, 2] полимеров в областях сложной формы. Для оптимизации подобных процессов требуется высокая точность моделирования течения; численные методы решения такого типа задач активно разрабатываются. Одним из основных вопросов является проблема моделирования течений, характеризуемых высоким значением числа Вайсенберга (high Weissenberg number problem, HWNP), заключающаяся в отсутствии сходимости для таких течений. Как показывает практика, ее решение осложняется тем, что предельное значение числа Вайсенберга весьма чувствительно к выбору способа дискретизации уравнений [3].
Метод погруженных границ [4] LS-STAG [5] c функциями уровня [6] позволяет снизить вычислительные затраты на выполнение моделирования в областях сложной формы. Метод LS-STAG и его модификации могут быть применимы к решению широкого класса задач, включая сопряженные задачи гидроупругости [7]. Разработана также модификация метода LS-STAG для моделирования течений вязкоупругих жидкостей в случае неподвижных погруженных границ [8, 9]. Этот метод реализован в разрабатываемом авторами программном комплексе [8].
Цель работы — верификация созданной модификации метода LS-STAG для расчета течений вязкоупругих жидкостей, описываемых моделями скоростного типа. Модельная задача о плоском течении в канале с внезапным сужением 4:1 часто используется для тестирования алгоритмов, поскольку в области резкого изменения формы канала наблюдаются высокие уровни напряжений. Это означает, что такая задача хорошо подходит для оценки критического значения числа Вайсенберга, характерного для тестируемого численного метода. Кроме того, подобная конфигурация области широко распространена в технологии формования полимерных изделий на этапах транспортировки.
Постановка тестовой задачи. Плоское течение вязкоупругой несжимаемой жидкости в канале с внезапным сужением 4:1 рассматривается в прямоугольной расчетной области Q = [0;32D]x[0;8D] (рис. 1). Расстояние от входной границы до сужения 16D. В безразмерных переменных постановка задачи имеет вид
ду 1 ~
V-v = 0, — + (у-V)v + Vp--Av = V-Te, Te + XDxe =2ve(S + VDS);
8t Re (1)
v(x, y, 0) = vо (x, y), Te (x, y, 0) = t0 (x, y), (x, y) e Q;
v |г2=Vinflow = -4D (8D " 4D]ex' v 'Г1=v 'Гз=v K=
c)v
дП
Г4
= ^ dn
= 0;
г и k
|r1 U Г2 U Г3 = Te'bc(x, y, t), Xex't = 2Xv(
(du ^
dy.
(1)
re,bc =
yy —
0,
re ,bc = xy
ve£. f |k =°, £ dy an
= °.
Г4
Здесь t — время; x, у — координаты; p — давление; п — внешняя нормаль к границе области течения K = К иK2; V = у(х, у, t) = = u ■ е Л + v ■ е у — скорость; Re — число Рейнольдса; ve — вязкость вязко-упругой составляющей; Хг — время запаздывания; X — время релаксации; те — вязкоупругая составляющая тензора напряжений [10]; Б = (Уу + [Уу]т) /2 — тензор скоростей деформации; Б — оператор дифференцирования по времени. Модель скоростного типа определяется видом оператора Б и значением Хг. Так, если под дифференцированием по времени подразумевается частная производная по времени, то в зависимости от наличия константы Хг получаем линейные модели Максвелла [11] и Джеффри [12]. Если при дифференцировании по времени используются конвективные производные [2], то получаем квазилинейные модели. Так, в верхней конвективной модели Максвелла [11] и в модели Олдройда-Б [13] задействована верхняя конвективная производная
DX = I = ® dt
(v -V)X - X(Vv) - (Vv )т X;
(2)
в модели Максвелла-А [11] и в модели Олдройда-А [13] — нижняя
А яу
БX = X =-+ (V • У)Х + Х(Уу ) + (УУ )т X; (3)
dt
в модели Джонсона — Сигельмана [14] — вращательная
о 1 (V А
БХ = X = 21 X + X
Далее будет рассматриваться жидкость Олдройда-Б, так как для нее все константы отличны от нуля и в конвективной производной присутствуют все члены.
Vinflow г3 |з£» 1
—А
Г4
J ву 7 Гх
4D
Рис. 1. Расчетная область
Модификация ЬБ-БТАв-сетки. В области О строим «основную» прямоугольную сетку с ячейками - = (хг_ьхг' )х (у-у-), площади которых
равны V',- = Ахг-Ау-. Центрам ячеек соответствуют радиус-векторы xci ■ = = (хС, уС). Границы этих ячеек обозначим Гг-,- (рис. 2). Контрольный объем Ог,- используется для построения дискретных аналогов уравнения неразрывности и уравнений для нормальных неньютоновских вязкоупругих напряжений. Для построения дискретных аналогов уравнения импульса
в проекциях на Ох и Оу используются контрольные объемы ОЦ- = = (х\, хС+1) X (у-, у--+1) и - =(хг_1,хг )х (у С, у^+1) с границами ГЦ - и Ц -
У+1
Рис. 2. Разнесенные сетки
(см. рис. 2). Кроме того, для расчета касательных вязкоупругих напряжений
необходимо рассмотреть контрольные объемы = (х- ,хгс+1)х (ус_1, ус),
которые в результате составят дополнительную сетку, называемую ху-сеткой. Если каждой ячейке Ог-,^ «основной» сетки присвоить вес
0, если
(4)
«г, j = Н
твердая ячейка, г ,j — треугольная ячейка, 1/4 в остальных случаях,
г j
1/3, если Q
то площадь ячейки О.ХУу можно будет вычислять по формуле МХу =
Граница Г,ь твердого тела 0!ь, форма которого в общем случае может быть произвольной, задается функцией расстояния (функцией уровня)
Ф(г) [6]:
ф(г)< 0, г еП/ = П\{Ой иРь},
ф(г) = 0, г еГЙ,
ф (г) > 0, г .
Граница Г-Ь для всех усеченных ячеек О.-,^ ЬБ-БТАв-сетки аппроксимируется отрезком прямой. Для определения координат концов каждого такого отрезка используется линейная интерполяция значений функции уровня ф (х{, у-}). Коэффициенты заполнения ячеек Щ^, е[0,1] соответствуют доле жидкости на восточной и северной границах ячейки и позволяют однозначно идентифицировать тип усеченной ячейки. Отметим, что при решении плоских задач существует всего три типа усеченных ячеек: 1) трапециевидные; 2) треугольные; 3) пятиугольные ячейки. При этом принадлежность ячейки той или иной группе однозначно определяет местоположение точек вычисления неизвестных (рис. 3).
Построение дискретных аналогов конвективных производных. Верхняя (2) и нижняя (3) конвективные производные тензора
X =
Хх X
xy
X X
xy
yy.
содержат слагаемые
X(Vv) =
Xxx du ' Xxy du Xxx dv Xxy dv >
öx dy dx dy
Xxy du ~ X yy du Xyy dv ' Xxy dv
dx dy dy dx j
(
du
du
(Vv )т X =
Xxx + Xxy
dx dy
X 5v ^ X dv_
Xxx + Xxy
ox dy
X Ё1 ± X Ё1
Xxy + Xyy
ox dy X 5v ^ X av
Xyy + Xxy
dy dx,
dxjij, dy\i,j
Рис. 3. Северная трапециевидная ячейка (а), северо-западная треугольная ячейка (б) и северо-западная пятиугольная ячейка (б) ЬБ-БТЛО-сетки с отмеченными точками вычисления неизвестных
С учетом уравнения неразрывности получаем
( ( 2
X(Vv) + (Vv )т X =
X - + X ^
xxx ^ "Г" -Л-xy -
V
Xxx
du dx dv
+ X
dx f Y
x
v Yxy
yy
dy du dy Y Л
xy
Xxx ^ + Xyy _
X
dv dx dv
yy
dy
+ X
du >) dy dv
xy
dx
lyy у
Дискретные аналоги нормальных
du dx
dv
i, j 8у
и касательных напря-
U j
жении
du dy
j
dv dx
необходимые, в частности, для вычисления компо-
U j
нент тензора скоростей деформации, построены в [8, 9] для всех типов ячеек ЬБ-БТАв-сетки. Значения Ххх, Хуу, Ухх и Ууу, как и нормальные напряжения, вычисляются на «основной» сетке, а Хху и Уху, как и каса-
тельные напряжения, — на дополнительной ху-сетке. Поэтому при вычислении Ухх и Ууу следует использовать средние по О,,у значения
Хху — и Хху —, которые обозначим / Хху —\ и \Хху —) . С уче-
ду дх \ ду/.г} \ дх\{,
том (4) среднее по О,,у значение величины ы с ху-сетки равно (ы),. = а,-,]Ыц,] + ы,,н + у + ^¿-1,;_1).
При этом на усеченных ячейках О.,,у часть значений ы,,у_1, у, ы,,у, Ыг_1,берется из граничных условий. При интегрировании Уху необходимо учитывать различие в сетках:
I учйу=(Ххх)*у |х +(Хуу|у
xy j dxitj " 77/1,j ay
оу.
i, j
Эта формула выполняется для всех типов ячеек 0,ху. Здесь (Xхх}ху,
(Хуу)ху — средние по ^у значения Ххх и Хуу. С использованием (4)
среднее по 0.ху значение величины г с «основной» сетки вычисляется следующим образом:
= а-, №, г, ^ +а,, у, )-1 +а,_1, ¡У1 _1, у г,-1, ^ +а,_1, у-1^-1,у-1,
Кроме того, для квазилинейных моделей вязкоупругих жидкостей необходимо построить дискретные аналоги интегралов от слагаемых (у -У)Х, входящих во все конвективные производные. Для уравнений импульса дискретизация конвективных членов приводила к кососиммет-ричным матрицам [5] вследствие использования центральных разностей. Однако рассматриваемые уравнения переноса вязкоупругих напряжений являются гиперболическими, поэтому использование центральных разностей при построении дискретного аналога конвективного члена приведет к вычислительной неустойчивости [15]. Возможным решением указанной проблемы является использование направленных разностей против потока [15].
Для Хху с учетом формулы Остроградского — Гаусса получаем
J (v ■W)XxydV = J Xxyv -ndS.
nxy г xy
t, j t, j
Для заполненной жидкостью ячеИки Qxy имеем
i, j
J (v • n)X xydS ~ C^y[U]w(i, j)Xxy |i-1,j + Cexy[U]e(i, j)Xxy |i+1,j +
rxy i, j
+ сху[и ]р (/, -)Хху |/, - + сху[и ] (/, -)Хху |г,-_1 + С^у[и ] (/, -)Хху |г,-+1. (5)
Для Пгху (рис. 4) конвективный член можно записать как сумму потоков через четыре ее элементарные грани:
I (v • n)XxydS = - | (v • eх)Ххуйу + | (v • eх)Ххуйу-
rxy i, j
1 i, j+1ui i, j
Tsw,e , rnw,e 1 i+1, j+1ui i+1, j
{ (v • e y )Xxydx + j (v • e y )Xxydx.
(6)
Tne,s ,rnw,s
4j uli+1,j
Tse,n rsw,n 1 i, j +1i+1, j +1
а б
Рис. 4. Контрольный объем Q*y (а) и ячейка Qi,j «основной» сетки (б)
Поток массы через жидкую грань Г?- (см. рис. 4) равен щ,-
yj
= J v • exdS = J u(xi, y)dy. Поскольку точка вычисления ui,j располага-
re. i,)
yj-1
ется в центре жидкой части грани (см. рис. 3), имеем
Г 1
ui, j = u
xi, yj_1 н—^uujAyj |. Применяя квадратурную формулу централь-
ных прямоугольников, получаем и/,- -. Разностный аналог по-
тока, относящегося к твердой грани, имеет вид
-ib
Ui» uibj [nxAS]* + vfbj [nyASfbj.
ib
Здесь [«хДЭДЬ = у - Цу) Дуу; [«уАУ^у = у_1 - у) Ах,; скорость у¿Ьу = (и^, У^) на Гу вычисляется на усеченной ячейке в концах отрезка Гпо граничным условиям. Тогда разностный аналог (6) можно вычислить по формуле:
rxy
1 ( ui-1, j + ui, j ^ ui, j+1 + ui-1, j+1 { u : + u
i (v• n)XxydS Zj^ + j+1 ' -1,j+1 j [XxyVj +
[ Xxy ]f, j -
xy Ji, j
2 у 2 2 J
i, j
Ui, j + Ui+1, j Ui+1, j+1 + Ui, j +1
+ -
2 1
V
f V
Vi, j1 + Vi, ^ + Vi+1, j +Vi+1, j-1 1 [ Xxy 1 j +
2 ч 2
1 + -2 ( V- ■ Vi, j + Vi, j+1 |
\ 2
1,;+12 [Хху «, (7)
[ Х ]Ы = 1 + Щу + и,_1,у + и,_1,у+1 + Ц,у+Х) Х |
[Хху],,у = 2 Хху |г-1,у +
1 - и,,у + иг-1,у +иг-1,у+1 + и,,у+1) Х | + 2 Хху |г,у,
[ Х ]е = 1 - и,,у + и,+1,у +й;+1,у+1 + и,,у+1) Х |
[Хху],,у = 2 Хху |г+1,у +
1 + и,у + иг+1,у +иг+1,у +1 +и,,у+1) Х |
+ 2 Хху |;,у,
[ Х ]5 = 1 + чу + Уг,у-1 + у+1,]-1 +1,у ) Х |
[Хху ],,у = 2 Хху у-1 +
1 - У,у + У>,у_1 + У{+1гу_1 + Ущ,у ) Х | + 2 Хху I,,у,
[ Х « = 1 - в§п(У,,у + у+1 + У'+1,у+1 +1,у ) Х |
[Хху ],,у = 2 Хху у+1 "+"
1 + У,у + У,,у+1 + У,+1,у+1 + У{+1гу ) Х
+ 2 Хху ,,} •
Сравнивая (7) и (5), получаем формулы для вычисления элементов Сеу[и ], например:
сху[ик а, ;) =
1
( ui_1, j + ui, j ui, j+1 + ui_1, j+1 ^
1 + sgn(ui_1, j + ui, j + ui, j+1 + ui-1, j+1)
Для усеченных ячеек в формулах следует учитывать значения у/ь (х, у)
—гЬ г п
и потоков и/,- по аналогии со случаем ньютоновской жидкости [5].
Для компонент Ххх и Хху разностные аналоги конвективных членов должны строиться по другим формулам, поскольку в отличие от Хху эти компоненты вычисляются на «основной» сетке. Построим разностный аналог интегралов величины {у-У)г для г с «основной» сетки. С учетом формулы Остроградского — Гаусса | (у • n)zdV = | ху ■ ndS.
О- Г- •
г, - ', -
Для заполненной жидкостью ячейки О/,- имеем
I (у■ n)zdS « Сг[и]^(г, -)х/-1,- + С?[и]в(г,-)х/+1,- +
Г/, -
+ Се [и]р(г, -)хг,- + Се [и]5 (г, -)х/,-_1 + Се [иN (г, ))хи-+1. (8)
Для о/ - (см. рис. 4) дискретизация состоит в записи конвективного члена как суммы потоков через четыре его элементарные грани:
| (у • n)zdS = - | (У • ех )zdy + | (У • ех )zdy -
г,, - г". ге.
'' t, - г,-
- | (У • еу )zdx + | (У • еу )zdx. (9)
уП
Тогда разностный аналог (9) можно вычислить по формуле
{ (У • п^Б « -ц-1,- [х]у- + йи- {х}? - - ъ,--1{х}1 - + Ъ,- {Х}п-, (10)
{х }?, - = {х }1+1, -, {х }п - = {х -+1,
= 1 + 5§П Щ-1,- х + 1 - 5§Лм/'-1,- х {Х}г, - = 2 1, - 2 Х/, -,
{Х}г, - = 2 ,2 Х/, -.
Сравнивая (10) и (8), получаем формулы для Се [и], например Се [ик (г, -) = -Щ.
Интегрирование по времени. После построения описанных дискретных аналогов имеем дифференциально-алгебраическую систему, которую необходимо численно проинтегрировать по времени. Для этого
построим метод на основе схемы предиктор-корректор первого порядка точности. Шаг предиктора заключается в решении разностного аналога уравнения Гельмгольца для прогноза скорости и в момент времени:
м(й~йп)+с [ип ]ип+$ь'с'п - в т(рп - тент) -
М
2 Re
,П
- Dxryn ~ (KU + Sib'v) = 0.
Здесь Рп — вектор, состоящий из значений ре, п — верхний индекс, указывающий на принадлежность к моменту времени tn = nAt; Т£'(£т — вектор
е \п е \п те п - е 1п
с компонентами тХх . и туу I ; ±Ху — вектор, состоящий из тХЛ .;
n,} n,} п,)
С [и] и К — матрицы, полученные в результате построения дискретных аналогов конвективных и вязких потоков; ип — вектор с компонентами ип] и уп^; 8'Ь,с,п и Б1Ь,У — векторы, соответствующие источниковым членам, которые появляются вследствие граничных условий; м — диагональная матрица, формируемая из значений площадей ячеек и О.У^;
—Вт — матрица, возникающая в результате построения дискретного аналога оператора градиента; Вт — матрица, которая соответствует дискретному аналогу оператора дивергенции, построенному на ху-сетке.
Шаг корректора соответствует решению разностного аналога уравнения Пуассона:
АФ = Бй + йЬ'п+1, Л = -ВМ"1Вт.
В результате становится возможным вычислить скорости и давление при ^+1 =(п + l)At:
йп+1 = й + М"ВтФ, Рп+1 = Рп + (Ф / Дt).
Затем выполняется расчет вязкоупругих напряжений:
М (те,п+1 _ те,п )
МеТе,п+1 + X МеТ-=
At
= 2ЧеМеЗп+1 + 2^Ме(5п+1 " ^ ) + СВ(йп+1, Хп ).
At
Здесь Ме — диагональная матрица, формируемая из значений площадей ячеек О.], ^ и О.ХУ; ^ — вектор, составленный из значений компонент тен-
зора скоростей деформации в соответствующих точках; Хп = = 2уДг5п+1 -ХТе,п; СО(ип+1, Хп) соответствует дискретному аналогу
конвективной производной Хп, из которой вычтена частная производная по времени.
Вычислительные эксперименты. При решении задачи (1) полагаем, что характерный размер О = 1, характерная скорость Vх> =1. Тогда имеем
У^О = 1 = 1 = _р_
V У5 + Уе V*
ЛАТ XVcc ,
We =-= а,
D
Re = -
v
Здесь We — число Вайсенберга; — вязкость ньютоновской жидкости-растворителя; Р = 1/9 — доля вязкости растворителя. Таким образом,
1 -р
v* =•
Re
и X r = pWe.
Расчеты проводились на двух неравномерных сетках размерами 96 х 88 и 192 х 176, которые далее обозначим М1 и М2 соответственно. Обе сетки содержат 39,2 % твердых ячеек. В окрестности углов используются блоки равномерной сетки с шагом И (рис. 5).
уЮ
|111111111111111|||111111Ш||1-: \
О 5 10 15 20 25 30 */£>
Рис. 5. Неравномерная сетка М1 (твердые ячейки показаны серым цветом)
Результаты сравнения предельных значений Wecr числа Вайсенберга, определенных при верификации модифицированного метода ЬБ-БТАв и указанных в работах [3, 16], в которых расчеты проводились с помощью различных реализаций метода конечных элементов (МКЭ), приведены в табл. 1. Для всех методов и сеток указаны значения шагов по пространству и времени, поскольку при измельчении сетки значение Wecr может уменьшаться [3]. Во всех работах использовались неравномерные по пространству сетки, причем наиболее подробные блоки сеток применялись в окрестности угла, выступающего в месте сужения. Это вызвано тем, что численная неустойчивость при расчете высоковязких течений связана именно с этой особенностью области. Поэтому в табл. 1 представлены
значения к / Л, используемые в окрестности этой точки. При расчетах течений, характеризуемых значением числа Рейнольдса меньше единицы, наиболее высокое значение Weсг демонстрирует метод, используемый в [16]. Однако при увеличении числа Re значение Weсг для этого метода начинает уменьшаться, в то время как для метода ЬБ-БТАв оно остается без изменений.
Таблица 1
Результаты сравнения предельных значений числа Вайсенберга
Параметр Сетка
M1 М2 UM3 [3] Mesh III [16]
h / D 0,06250 0,03125 0,00500 0,01500
At 0,00500 0,00500 0,00050 0,00100
We cr (Re < 1) 3,00000 3,00000 2,50000 3,60000
WeCr (Re > 1) 3,00000 3,00000 - 2,50000
Результаты измерения длины вихря Ь, образующегося в углу области перед сужением (рис. 6), представлены в табл. 2. Расчеты проводились в диапазоне чисел Вайсенберга 0-3 и при значениях числа Рейнольдса 0,010; 0,100 и 1,000. В случаях когда значение Ь не определено в результате достижения Weсг, вместо числовых значений в соответствующих ячейках
14,0 14,5 15,0 15,5 16,0 16,5 17,0 х/Б
Рис. 6. Пример определения безразмерной длины вихря Ь на сетке М1 при \¥е = 3,000 и 11е = 0,010
табл. 2 указана аббревиатура HWNP. При значениях числа Вайсенберга, не превышающих 2,500, рассмотренные методы демонстрируют очень близкие результаты: разность результатов не превышает 3 %. Однако при значениях We от 2,500 и более разброс результатов сильно увеличивается и в рассмотренных случаях составляет от 7 до 16 %. Однако модификация метода
ЬБ-БТАв позволяет смоделировать монотонное уменьшение значения Ь при увеличении числа We, о котором идет речь в [3], в то время как в [16] данный эффект не был воспроизведен. Это объясняется высокой схемной искусственной вязкостью используемого численного метода, которая сильно отражается на точности расчета при высоких значениях числа Вайсен-берга.
Таблица 2
Результаты сравнения безразмерной длины вихря
We Re LS-STAG МКЭ
M1 М2 [3] [16]
0 0,010 1,720 1,490 - 1,503
0,100 1,660 1,450 1,448 -
1,000 1,350 1,240 - 1,245
0,500 0,010 1,490 1,450 - 1,454
0,100 1,360 1,340 1,339 -
1,000 1,220 1,180 - 1,167
2,000 0,010 1,360 1,350 - 1,357
0,100 1,120 1,100 1,073 -
1,000 0,880 0,860 - 0,852
2,500 0,010 1,190 1,320 - 1,448
0,100 0,800 0,900 0,966 -
1,000 0,790 0,840 - 0,851
3,000 0,010 1,400 1,300 - 1,543
0,100 0,820 0,840 HWNP -
1,000 0,780 0,830 - HWNP
Заключение. Описаны основные особенности модификации метода погруженных границ ЬБ-БТАв. Такая модификация позволяет моделировать течения высоковязких жидкостей. Получены формулы для вычисления разностных аналогов верхней, нижней и вращательной производных. В целях верификации разработанного численного метода рассмотрена модельная задача о моделировании плоского течения жидкости Олдройда-Б в канале с внезапным сужением 4:1. Вычислительные эксперименты показали, что на этой задаче предельное значение числа Вайсенберга для разработанного численного метода равняется 3, что является хорошим показателем для теста. Кроме того, измеренные значения безразмерной длины смоделированного вихря хорошо согласуются с результатами расчетов, полученными другими авторами.
ЛИТЕРАТУРА
[1] Owens R.G., Phillips T.N. Computational rheology. Imperial College Press, 2002.
[2] Galdi G.P., Robertson A.M., Rannacher R., et al. Hemodynamical flows: modeling, analysis and simulation. Oberwolfach Seminars, vol. 37. Birkhauser Basel, Springer, 2008. DOI: https://doi.org/10.1007/978-3-7643-7806-6
[3] Kim J.M., Kim C., Kim J.H., et al. High-resolution finite element simulation of 4:1 planar contraction flow of viscoelastic fluid. J. Nonnewton. Fluid Mech., 2005, no. 129, iss. 1, pp. 23-37. DOI: https://doi.org/10.1016/j.jnnfm.2005.04.007
[4] Mittal R., Iaccarino G. Immersed boundary methods. Annu. Rev. Fluid Mech., 2005, vol. 37, pp. 239-261. DOI: https://doi.org/10.1146/annurev.fluid.37.061903.175743
[5] Cheny 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, vol. 229, iss. 4, pp. 1043-1076. DOI: https://doi.org/10.1016/j.jcp.2009.10.007
[6] Osher S., Fedkiw R.P. Level set methods and dynamic implicit surfaces. Applied Mathematical Sciences, vol. 153. New York, Springer, 2002.
DOI: https://doi.org/10.1007/b98879
[7] Marchevsky I.K., Puzikova V.V. Numerical simulation of wind turbine rotors autorotation by using the modified LS-STAG immersed boundary method. Int. J. Rotating Mach, 2017, vol. 2017, art. 6418108. DOI: https://doi.org/10.1155/2017/6418108
[8] Пузикова В.В. Модификация метода погруженных границ LS-STAG для моделирования течений вязкоупругих жидкостей. Труды ИСП РАН, 2017, т. 29, № 1, с. 71-84. DOI: https://doi.org/10.15514/ISPRAS-2017-29(1)-5
[9] Botella O., Cheny Y., Nikfarjam F., et al. Application of the LS-STAG immersed boundary/cut-cell method to viscoelastic flow computations. Commun. Comput. Phys., 2016, vol. 20, iss. 4, pp. 870-901. DOI: https://doi.org/10.4208/cicp.080615.010216a
[10] Уилкинсон У.Л. Неньютоновские жидкости. М., Мир, 1964.
[11] Maxwell J.C. On the dynamical theory of gases. Philos. Trans. R. Soc. Lond., 1867, vol. 157, pp. 49-88. DOI: https://doi.org/10.1098/rstl.1867.0004
[12] Jeffreys H. The Earth. Its origin, history and physical constitution. Cambridge Univ. Press, 1929.
[13] Oldroyd J.G. On the formulation of rheological equations of state. Proc. Roy. Soc. London, 1950, vol. 200, iss. 2063, pp. 523-541.
DOI: https://doi.org/10.1098/rspa.1950.0035
[14] Johnson M.W. Jr., Segalman D. A model for viscoelastic fluid behavior which allows non-affine deformation. J. Nonnewton. Fluid Mech., 1977, vol. 2, iss. 3, pp. 255270. DOI: https://doi.org/10.1016/0377-0257(77)80003-7
[15] Роуч П. Вычислительная гидродинамика. М., Мир, 1980.
[16] Li X., Han X., Wang X. Numerical modeling of viscoelastic flows using equal low-order finite elements. Comput. Methods Appl. Mech. Eng., 2010, vol. 199, iss. 9-12, pp. 570-581. DOI: https://doi.org/10.1016/jxma.2009.10.010
Марчевский Илья Константинович — канд. физ.-мат. наук, доцент кафедры «Прикладная математика» МГТУ им. Н.Э. Баумана (Российская Федерация, 105005, Москва, 2-я Бауманская ул., д. 5, корп. 1); старший научный сотрудник ИСП РАН (Российская Федерация, 109004, Москва, ул. А. Солженицына, д. 25).
Пузикова Валерия Валентиновна — канд. физ.-мат. наук, старший научный сотрудник ИСП РАН (Российская Федерация, 109004, Москва, ул. А. Солженицына, д. 25).
Просьба ссылаться на эту статью следующим образом:
Марчевский И.К., Пузикова В.В. Использование модифицированного метода LS-STAG для расчета плоского течения вязкоупругой жидкости в канале с внезапным сужением 4:1. Вестник МГТУ им. Н.Э. Баумана. Сер. Естественные науки, 2021, № 3 (96), с. 46-63. DOI: https://doi.org/10.18698/1812-3368-2021-3-46-63
THE MODIFIED LS-STAG METHOD APPLICATION FOR PLANAR VISCOELASTIC FLOW COMPUTATION IN A 4:1 CONTRACTION CHANNEL
I.K. Marchevsky1' 2 iliamarchevsky@bmstu.ru
V.V. Puzikova2 vvp@dms-at.ru
1 Bauman Moscow State Technical University, Moscow, Russian Federation
2 Institute for System Programming, Russian Academy of Sciences, Moscow, Russian Federation
Abstract
In this study we present the modification of the LS-STAG immersed boundary cut-cell method. This modification is designed for viscoelastic fluids. Linear and quasilinear viscoelastic fluid models of a rate type are considered. The obtained numerical method is implemented in the LS-STAG software package developed by the author. This software is created for viscous incompressible flows simulation both by the LS-STAG method and by it developed modifications. Besides of this, the software package is designed to compute extra-stresses for viscoelastic Maxwell, Jeffreys, upper-convected Maxwell, Maxwell-A, Oldroyd-B, Oldroyd-A, Johnson — Segalman fluids on
Keywords
Immersed boundary method, the LS-STAG method, rate type viscoelastic flow, incompressible flow, contraction channel, convective derivative
the LS-STAG mesh. The construction of convective derivatives discrete analogues is described for Oldroyd, Cotter — Rivlin, Jaumann — Zaremba — Noll derivatives. The centers of base LS-STAG mesh cells are the locations for shear non-Newtonian stresses computation. The corners of these cells are the positions for normal non-Newtonian stresses computation. The first order predictor-corrector scheme is the basis for time-stepping numerical algorithm. Benchmark solutions for the planar flow of Oldroyd-B fluid in a 4:1 contraction channel are presented. A critical value of Weissenberg number is Received 18.08.2020 defined. Computational results are in good agreement Accepted 15.09.2020 with the data known in the literature © Author(s), 2021
This work was supported by the Russian Science Foundation (RSFproject no. 17-79-20445)
REFERENCES
[1] Owens R.G., Phillips T.N. Computational rheology. Imperial College Press, 2002.
[2] Galdi G.P., Robertson A.M., Rannacher R., et al. Hemodynamical flows: modeling, analysis and simulation. Oberwolfach Seminars, vol. 37. Birkhäuser Basel, Springer, 2008. DOI: https://doi.org/10.1007/978-3-7643-7806-6
[3] Kim J.M., Kim C., Kim J.H., et al. High-resolution finite element simulation of 4:1 planar contraction flow of viscoelastic fluid. J. Nonnewton. Fluid Mech., 2005, iss. 129, no. 1, pp. 23-37. DOI: https://doi.org/10.1016/j.jnnfm.2005.04.007
[4] Mittal R., Iaccarino G. Immersed boundary methods. Annu. Rev. Fluid Mech., 2005, vol. 37, pp. 239-261. DOI: https://doi.org/10.1146/annurev.fluid.37.061903.175743
[5] Cheny 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, vol. 229, iss. 4, pp. 1043-1076. DOI: https://doi.org/10.1016/j.jcp.2009.10.007
[6] Osher S., Fedkiw R.P. Level set methods and dynamic implicit surfaces. Applied Mathematical Sciences, vol. 153. New York, Springer, 2002.
DOI: https://doi.org/10.1007/b98879
[7] Marchevsky I.K., Puzikova V.V. Numerical simulation of wind turbine rotors autorotation by using the modified LS-STAG immersed boundary method. Int. J. Rotating Mach., 2017, vol. 2017, art. 6418108. DOI: https://doi.org/10.1155/2017/6418108
[8] Puzikova V.V. The LS-STAG immersed boundary method modification for viscoelastic flow computations. Proceedings of ISP RAS, 2017, vol. 29, no. 1, pp. 71-84 (in Russ.). DOI: https://doi.org/10.15514/ISPRAS-2017-29(1)-5
[9] Botella O., Cheny Y., Nikfarjam F., et al. Application of the LS-STAG immersed boundary/cut-cell method to viscoelastic flow computations. Commun. Comput. Phys., 2016, vol. 20, iss. 4, pp. 870-901. DOI: https://doi.org/10.4208/cicp.080615.010216a
[10] Wilkinson W.L. Non-Newtonian fluids. Pergamon Press, 1960.
[11] Maxwell J.C. On the dynamical theory of gases. Philos. Trans. R. Soc. Lond., 1867, vol. 157, pp. 49-88. DOI: https://doi.org/10.1098/rstl.1867.0004
[12] Jeffreys H. The Earth. Its origin, history and physical constitution. Cambridge Univ. Press, 1929.
[13] Oldroyd J.G. On the formulation of rheological equations of state. Proc. Roy. Soc. London, 1950, vol. 200, iss. 2063, pp. 523-541.
DOI: https://doi.org/10.1098/rspa.1950.0035
[14] Johnson M.W. Jr., Segalman D. A model for viscoelastic fluid behavior which allows non-affine deformation. J. Nonnewton. Fluid Mech., 1977, vol. 2, iss. 3, pp. 255270. DOI: https://doi.org/10.1016/0377-0257(77)80003-7
[15] Roache P.J. Computational fluid dynamics. Hermosa, 1998.
[16] Li X., Han X., Wang X. Numerical modeling of viscoelastic flows using equal low-order finite elements. Comput. Methods Appl. Mech. Eng., 2010, vol. 199, iss. 9-12, pp. 570-581. DOI: https://doi.org/10.1016/j.cma.2009.10.010
Marchevsky I.K. — Cand. Sc. (Phys.-Math.), Assoc. Professor, Department of Applied Mathematics, Bauman Moscow State Technical University (2-ya Baumanskaya ul. 5/1, Moscow, 105005 Russian Federation); Senior Researcher, Institute for System Programming, Russian Academy of Sciences (A. Solzhenitsyna ul. 25, Moscow, 109004 Russian Federation).
Puzikova V.V. — Cand. Sc. (Phys.-Math.), Senior Researcher, Institute for System Programming, Russian Academy of Sciences (A. Solzhenitsyna ul. 25, Moscow, 109004 Russian Federation).
Please cite this article in English as:
Marchevsky I.K., Puzikova V.V. The modified LS-STAG method application for planar viscoelastic flow computation in a 4:1 contraction channel. Herald of the Bauman Moscow State Technical University, Series Natural Sciences, 2021, no. 3 (96), pp. 46-63 (in Russ.). DOI: https://doi.org/10.18698/1812-3368-2021-3-46-63
