CC BY
31
МАТЕМАТИЧЕСКАЯ ФИЗИКА. МАТЕМАТИЧЕСКОЕ
МОДЕЛИРОВАНИЕ MATHEMATICAL PHYSICS. MATHEMATICAL MODELING
УДК 519-7
ОСОБЕННОСТИ РЕШЕНИЯ МАТЕМАТИЧЕСКИХ МОДЕЛЕЙ ВЗАМОДЕЙСТВИЯ ЖИДКОСТИ И ДЕФОРМИРУЕМЫХ ТВЕРДЫХ ТЕЛ
FEATURES OF THE SOLUTION FOR FLUID-STRUCTURE INTERACTION
MATHEMATICAL MODELS
О.А. Гальцева, О.В. Гальцев О.А. Galtseva, O.V. Galtsev
Белгородский национальный исследовательский университет, Россия, 308015, г. Белгород, ул. Победы, 85 Belgorod National Research University, 85 Pobedy St, Belgorod, 308015, Russia E-mail:galtseva@bsu.edu.ru, galtsev_o@bsu.edu.ru
Аннотация
В данной работе рассматриваются известные подходы к решению задач однородных однофазных жидкостей. Проводится их модификация и адаптация для использования в решении задач взаимодействия двухфазных жидкостей и упругих тел. С помощью преобразования рассматриваемых уравнений в конечно-элементный вид и проведения пространственной и временной дискретизации было разработано научно -исследовательское программное обеспечение, позволяющее как оценивать напряжение на общей границе «жидкость-упругое тело», так и отслеживать поведение свободной границы раздела жидкостей с учетом структурных деформаций.
Abstract
In the paper, we consider well-known approaches for solving problems of homogeneous single-phase liquids. Its modification and adaptation for use in solving two-phase fluid-structure interaction problems is carried out. Transforming equations into a finite-element form and carrying out spatial and temporal discretization, was developed a research software, that allows to estimate the tension on the common boundary "liquid-elastic body", and to trace the behavior of the free interface taking into account the structural deformations.
Ключевые слова: гидродинамическое моделирование, флюидопотоки, жидкостно-структурное взаимодействие. Keywords: hydrodynamic modeling, fluid flows, fluid-structure interaction.
Введение
При численном моделировании процесса вытеснения одной жидкости другой в упругой пористой среде широко использовалось понятие кинематики. Под кинематикой понимается геометрическое описание движения. Тело определяется как совокупность материальных точек, где точка есть абстрактный геометрический объект с некоторыми физическими свойствами. Определим тело р в некотором состоянии к как рк, где текущее состояние есть состояние тела до деформации. Этот объект р можно определить с помощью функции деформации х. Положение материальной точки в к обозначим как X, а положение точки в деформированном состоянии (пространственная точка) как х. Функция деформации х переносит точку X в пространственную току х (рис. 1).
Можно определить движение по времени как
X(X,t) = x(X,t) =X + w, (1)
где w - вектор перемещения:
w = x(X,t)-X. (2)
Так как решаемая задача содержит жидкость и твердое тело, то необходимо учесть особенности вектора перемещений. Область жидкостно - структурного взаимодействия (Fluid-
Structure Interaction, FSI) можно определить как единую область, состоящую из твердой и жидкой подобласти fl = fls U Hf, где fls - твердое тело, flf - жидкость. Поверхность между ними - Г. Вектор перемещений w = ws б fls и w = wm б flf, где ws - перемещение твердого тела, wm - перемещение жидкой подобласти. Эти перемещения равны на общей границе «жидкость-твердое тело».
Рис. 1. Состояния тела Fig. 1. Body conditions
Градиент деформации определим следующим образом
F = V' x = V' w+II , (3)
где II - единичный тензор. Оператор дифференцирования по X обозначим как V', а градиент по х как V. Детерминант (определитель) F есть
J = det(F) . ^ (4)
Подавляющее большинство численных исследований жидкостно-структурного взаимодействия базируются на лагранжевом, эйлеровом и лагранжево-эйлеровом (Arbitrary Lagrangian-Eulerian, ALE) [3, 4] описании тел. В лагранжевом описании (обычно используется для описания упругого тела) отслеживается движение материальных точек в пространстве через границы эйлеровой ячейки (рис. 2).
Материальная точка Точка сетки
-Движение тела
----Движение области
f t
Рис. 2. Одномерная демонстрация лагранжевого описания. Точки сетки двигаются вместе
с материальными точками Fig. 2. One-dimensional demonstration of the Lagrangian description. The grid points move along
with the material points
В эйлеровом описании (обычно используется для описания жидкости) определяется, насколько материальные точки сдвинулись с течением времени по отношению к узлам сетки (рис. 3).
Материальная точка Точка сетки
-Движение тела
----Движение области
ut
Рис.3. Одномерная демонстрация эйлерового описания. Точки сетки остаются неподвижными в то время, как материальными точками перемещаются Fig. 3. One-dimensional demonstration of the Euler description. The grid points remain stationary
while the material points are moving
Задачи жидкостно-структурного взаимодействия требуют общего описания для уравнений, описывающих движение жидкости и твердого тела. Нами была выбрана третья формулировка, которая объединяет две предыдущие и позволяет расчетной сетке деформироваться в ответ на структурные изменения, а главным преимуществом является способность сохранять высокую точность вблизи границы раздела упругого тела и жидкости (рис. 4).
Рис. 4. Одномерная демонстрация произвольного лагранжево-эйлерового описания. Точки сетки перемещаются произвольно вместе с материальными точками Fig. 4. One-dimensional demonstration of an arbitrary Lagrangian-Euler description. The grid points move arbitrarily
along with the material points
В этом подходе некоторые точки могут совпадать с материальными точками, некоторые остаются без изменений, а остальные находятся между ними.
Основные уравнения основываются на законе баланса или сохранения с некоторыми дополнениями для конкретных веществ и кинематических отношений. Эти законы применяются для некоторой пространственной области , изменение массы в области равно массе на входе и выходе, изменение линейных импульсов в области равно результирующей силе, действующей на область.
Первый закон - это закон баланса (сохранения) массы:
J + V-(pv) = 0,
где р (х, 1) - плотность вещества, V (х, 1) - его скорость. Второй закон - баланс импульса:
DV Г7 Т* I
pm = V'T+Pg'
где g(x,t) - массовая сила, Dv/Dt - производная по времени
Dv 3v „
— = — + V • Vv.,
Dt at
(5)
(6)
(7)
Т(хД) = 2цЮ)(у) - рП.
В эйлеровой формулировке
Т = Тт.
Твердое вещество рассматривается как упругое тело, а жидкость есть вязкая и несжимаемая.
Рассматриваемые уравнения
Основное уравнение для упругого тела записывается в лагранжевой формулировке [2]:
Я ГО (ш)-р I ) + рg = 0,Xens, (8)
где - поле скоростей, а - симметричный тензор второго порядка вида
ГО (ш) = - 07'ш + (У'ш)т) . (9)
Введение второго поля (перемещения) требует условия на скорость и перемещение, чтобы замкнуть систему уравнений
dw
at
- v = 0,Х £ ns
Подставив (io) в (8), получим уравнение для w:
Р^Г - V • (2AJD)(w) - рП) + pg = 0,Х £ О,
(10)
(11)
Для численного решения задачи совместного движения жидкости и упругого тела используются обе формулировки, как формулировка двух функций (w,v) - уравнений (8) и (10), так и формулировка одной функции (w) - уравнение ( 11 ).
Приведем две формулировки уравнений для жидкой компоненты. Первая - эйлерово (естественное) описание жидкости, вторая - свободная лагранжево-эйлерова формулировка (ALE).
ALE-описание необходимо для моделирования жидкостно-структурнош взаимодействия и позволяет связать граничные условия жидкости и упругого тела.
Для вязкой несжимаемой жидкости система уравнений стостоит из уравнения Стокса в эйлеровой формулировке [1]:
Ц О (V) - р I) + pg = 0,х 6 Г (I) , (12)
где Г (1) - текущая деформируемая область жидкости, V - поле скоростей жидкости, р - давление, ц - вязкость, а О (V) - симметричный тензор второго порядка вида
D(v)=-(Vv+(Vv)T)
(13)
и условия несжимаемости
V • v = 0, х 6 nf (t) . (14)
Как говорилось ранее, произвольная лагранжево-эйлерова формулировка необходима для описания жидкостно-структурного взаимодействия в виде единой системы. Для этого перепишем уравнения (12) и (14), используя следующие отношения:
Vv = V' vk F " i (15)
и
£ = ^"Vm-(V' Vk F " i) , (16)
где k - индекс, обозначающий переменные в ALE формулировке, vm - скорость области жидкости, которая определяется как vm = дх/ дt. Таким образом, уравнение (12) примет вид
р^- V • (2ц D(vk) - pkI) + pgk = о,х е nf(0),
(17)
и
вид
( Р - ^ ' ^) = 0 ,Х6Г ( 0 ) , (18)
где Г ( 0 ) - рассматриваемая область, занятая жидкостью. Тензор напряжения будет иметь
0(vk)=i(V'vkF-1 + F-T(V'vk)T).
(19)
Теперь, когда система уравнений Стокса записана в произвольной лагранжево-эйлеровой формулировке, возможно совместить ее с уравнениями для упругого тела и решить задачу жидкостно-структурного взаимодействия.
Двухфазное течение можно реализовать различными способами, например, методом установления уровня жидкости (level-set метод) или методом объема жидкости (volume of fluid). В обоих методах решается транспортное уравнение:
Эф
at
+ v • Vc|) = 0, х 6 nf,
(20)
где ф - транспортная функция. В методе установления уровня поверхность между двумя жидкостями определяется как заданная изоповерхность. Следовательно на каждом шаге по времени фаза в каждой расчетной ячейке определяется простым скалярным сравнением.
Для уравнения (20), как и для уравнения (17), необходима ALE-формулировка для поиска значений деформированной расчетной области
^+(vk-vm)^ = o, xefif, (21)
где - скорость расчетной сетки и есть шаг по времени.
Следует отметить,что уравнения (17), (18) и (21) представляют собой связанную систему, где существует взаимное влияние level-set функции на поле скоростей жидкости и наоборот. Если эта задача дополняется третьей фазой (упругое тело), то жидкости будут влиять на перемещения упругого тела [5], а тело в свою очередь будет действовать на сетку жидкой подобласти (рис. 5).
Рис. 5. Взаимодействие двухфазной жидкой компоненты и упругого тела. Fluid - жидкость^гасШге - упругое тело, Mesh - сетка, Interface - граница раздела жидкостей
Fig. 5. Interaction of a two-phase liquid component and an elastic body. Fluid - the liquid, Structure - the elastic body, Mesh - the mesh, Interface - the interface of liquids
Граничные и начальные условия
В области П = П^ и и у, где у - граница раздела областей П^ и Г^т , численно решалась задача, состоящая из системы уравнений (12)-(14),(2о) в Г^ (рис. 6).
D.J D.J
а.
У У
Of \ / ilf
г
Рис. 6. Рассматриваемая область,граничные условия и системы уравнений (8),(9) в П s, дополненных граничными условиями Fig.6. The considered region, boundary conditions and systems of equations (8), (9) in П supplemented by the boundary conditions
(2 ц ГО (V) - р I ) п = - р±п, х 6 Б*, (22)
( 2 цГО (V) - р I ) п = (ЯГО (ш) - р I ) п, х 6 Г, (23)
V • п = 0, п • (2 ц ГО - р {I ) п = 0 ,х 6 Б0, (24)
где - граница симметрии.
Система (8), (9), (12),(14), (20) дополняется начальными условиями:
ш | 4=0 = 0, х 6 П, V | 4=0 = 0 , х 6 П, ф | 4=0 = 0, х 6 П. (25)
Перемещение упругого тела определяет перемещение расчетной сетки на поверхности раздела жидкость-упругое тело:
шт = ш^хег, (26)
где шт - перемещение расчетной сетки.
Метод конечных элементов
Метод конечных элементов [6, 7, 8] используется для решения начально-краевых задач, записанных в частных производных, предварительно преобразованных в алгебраические уравнения, которые аппроксимируют начальные уравнения. Метод взвешенных невязок используется для преобразования дифференциальных уравнений в интегральный вид (слабая формулировка). Пространственная дискретизации слабой формулировки позволяет получить систему обыкновенных дифференциальных уравнений во времени (полудискретная форма). Временная дискретизация этих полудискретных уравнений приводит к полностью дискретной системе алгебраических уравнений [9, 10]. Если она нелинейна, то ее необходимо линеаризовать.
Для записи слабой формулировки необходимо умножить рассматриваемые уравнения на произвольную пробную функцию вида и проинтегрировать каждое слагаемое по
области.
Таким образом слабая форма уравнения упругости будет иметь следующий вид:
/п5 р"^ • ш - /п8 ^ (2^ ГО (ш) - р о • ш - /п5 рg • ш = 0. (27)
Используя интегрирование по частям, тензор напряжений примет вид:
/_ ^ (2 ц ГО (ш)-р !)-ш = ; (ГО (ш)-р I )^ш, (28)
где t = ( ГО (w) — р I ) • п. Подставив (28) в (27), получим:
Jns pS • w + Jns (2 n Го (w) — P 1 ): Го (w) — ;ns pg • w — /^ t • w = 0, (29)
Проделаем такую же процедуру для уравнений (8) и (10), взяв пробные функции v и w соответственно:
Jns PÍ • v + /ns (2 ^ ГО (w) — P 1 ) : Го (v) — /fls Р sg • v — /Ms t • v = 0. (30)
ins?r^w—/nsv^w=0, (31)
В слабой формулировке система уравнений Стокса будет иметь вид
/* PÍ • v + 2 ^ (v) : ГО (v) — р div v = pg • v + t • v (32)
(33)
где p - пробная функция давления.
Для формы записи ALE процедура будет такая же, за исключением того, что необходимо учесть переход от эйлерового вида в ALE форму
W« a^Vk (х) ы ^ ^ (34)
Далее . Взяв такую же пробную функцию, как и для уравнений в эйлеровой
формулировке, получим слабую формулировку уравнений (17) и (18)
/* I Р1Т •'v + к J (2 ^ГО (vk) — Pk 1 ) F "Т : Го (v) — J pgk • v — /anf • v = 0 (35)
(J F" 1 :V' vT) • p = 0. (36)
Пространственная дискретизация
Смысл дискретизации состоит в разбиении ограниченной области на отдельные кусочно-непрерывные подобласти. Совокупность таких элементов определяется как сетка или триангуляция, которая является моделью исходной области. Основные переменные аппроксимируются узловыми значениями, умноженными на базисный вектор, который необходим для интерполяции внутри элемента. Наша аппроксимация будет иметь следующий вид:
w * 1 w Wj, v * 2f= 1 v Vj, P * EJL 1 PjPj , (37)
где и - -е базисные вектора, - количество узлов, а и - -е значения в узлах,
связанные с первичными переменными. В методе Галеркина набор базисных векторов берется такой же, как и набор взвешенных векторов. Таким образом, взвешенные функции можно определить следующим образом:
w = 2 n= 1 w w^T = 2 n= 1 vV , A = 2 n= 1 p P. (38)
Узловые значения, связанные со взвешеными функциями, являются произвольными, поэтому ими можно пренебречь.
Временная дискретизация
Для аппроксимации по времени рассматриваемых уравнений используются конечные разности. Общее уравнение после пространственной дискретизации будет иметь вид:
^ = g (t) . (39)
Интегрируя обе части по временному шагу
í+1í?dt = /tt;+1 g (t) dt. (40)
Первый интеграл можно записать как
/п+1^dt = vn+^ vn. (41)
■'tn dt
а второй интеграл можно аппроксимировать только использую среднее значение для некоторого времени :
(42)
где . Тогда явная левая и неявная правая эйлерова аппроксимация будут иметь вид
vn+i_vn vn+i_vn
= g (tn+1), = g (tn) . (43)
Для того, чтобы совместить левый и правый эйлеровы методы, воспользуемся - методом:
„П+1..П
= 9g (tn+1 ) + (1 — в) g(tn ) , (44)
где или . Особый случай, когда , можно аппроксимировать как
= 0 . 5 (g (tn+1) + g(tn ) ) , (45)
Д1
чтобы удовлетворить условие второго порядка точности по времени.
Используя 0-метод, перепишем уравнения 3 0 и 3 1 в следующем виде:
^n+^n Tn+l.v+f tn+l.v)
Jns г At v Jns R Jans )
г л
L р-
+ (1 _ 0) J"n¡. Tr • V + fgn^ tn • v),
i+i_wn
• w = 0 / Vn+1 • w + (1 - 0) / V" • w,
[Wsl "RWS "
Vs vs
Vf RVf
Pf RPf
Wm D wm
Jns P At ' ' "■'Os ' ■ ' ' ( " 0 )Jns ' ■ ' ' (47)
где TR = 2 цШ (v) — p I . Здесь для краткости не дискретизированны по времени, однако производные по времени находятся с использованием 0 -метода.
Уравнения жидкостно-структурного взаимодействия
Систему уравнений жидкостно структурного взаимодействия можно представить в следующем виде:
rjwsws JwsVs JwsVf JwsPf Jwswm Ivsws Ivsvs Ivsvf Ivspf Ivswm
JvfWs K'|VS JvfVf Jv|P| JvfWm
Ipfws IpfVs Ipfvf Ipfpf IpfWm Iwmws Iwmvs JwmVf Iwmpf Iwmwm
где переменные упругой компоненты обозначаются индексом s , жидкой компоненты - f, а сеточные переменные т. Нелинейные члены линеаризуются с помощью метода Ньютона-Рафсона, основанного на определении якобиана:
Jhk - (48)
где есть -й остаток и - некоторая неизвестная.
Эта система возникает путем нахождения остатков дискретизованных уравнений (ошибок дискретизации) и применения разложения в ряд Тейлора этих остатков.
Твердое тело исследуется в формулировках двух переменных уравнений упругости, где остатки находятся как
R1ws = /ns ir • ws — Jns v • wS, (49)
R1vs = J"ns P^ 5 • Vs + JOs (2 H Ш (ws) — p I ) : Ш (V) — Jr ts • vi (50)
Положение границы <<жидкость-упругое тело>> определяется нормальным напряжением
ts = ( Hf(V' VkF - 1 + FJT V' v£) — pk I) • (J F- Tn) . (51)
Жидкость описывается системой уравнений Стокса в ALE формулировке со следующими остатками:
RVf = Jo. ÍJ Pf^Vr • V) + Jn„ (J ( Hf(V' Vk F - 1 + F -TV' vT) — Pk I ) F -T) V' vj, (52)
Rpf = /Of (J F - W vT)^1, (53)
(54)
и
где индекс k указывает на принадлежность к области ALE. Остаток для расчетной сетки:
R1wm = ¡nf ^m(V'wm + V'wJ1)V'wjn-Используя эти определения и применив уравнение (48), можно определить якобианы. Многие якобианы в общей системе равны нулю:
Jwsws
Jvsws
о о о
Jwsvs
Jvsvs
о о о
о
Ivsvf Ivfvf Ipfvf
о
о
JvsPf
Ivfpf
о о
о
Jvswm
Ivfwm Ipfwm
Jwmwm
rwsi RWS "
Vs RVS
Vf RVf
Pf RPf
Wm D wm
Для удобства обозначим () ' = ^ () / При оценке якобиана необходимы следующие тождества:
г = 11Г (Р 'Р -1) , (55)
( Р -т) ' = ( (Р - ^ ' ) т, (56)
( Р - ^ ' = - Р - 1 Р 'Р - 1 , (57)
Р ' = V' ^ (58)
Аппроксимируем скорость движения сетки в жидкой компоненте как Ут = (х —хп) / следовательно
Якобианы твердой части будут иметь вид:
jL
At
I4 = Г Р -
где
I I = (60)
= (б1)
1 = х, V V*, (б2)
V, V =/п Р (6з)
Т^ = V Ш (2 ц5Ш (ш5) + р 5 I ) (I + (2 ц5Ш (Ш) + р 5 I ) . (64)
Якобианы жидкости будут иметь следующий вид:
4г = -/пг "Т V' V, (65)
(66)
^ = 4 1РГ-5 • ^ + /пг МС^ЧР"1) • ч + (I ^Р - -Т) V' V + (I К Р -Т V' V Р-Т) V' V. (67)
Якобиан расчетной сетки
1™ = (?'Шт + V' ШтТ) V' «С (68)
И якобианы взаимодействия жидкости и расчетной сетки будут иметь вид = /Пг Ур^ ■ V} - /Пг (ГркР"Т)^ - /Пг 0Рк(Р—+
;Пр (I ^ ( р -Т) 'V Чт Р -т) V' V+;Пр о ^ Р -Т V ЧТ (Р -т) ' ) V' V, (69)
Ст = 4 (1'Р - ^ V' чТ) • Р; + (I (Р - Х) '= V' чТ) • р. (70)
Ключевым моментом здесь является то, что пробные функции также должны быть непрерывными, так как они определены с использованием таких же базисных векторов, как и для скоростей и перемещений. Это означает, что vf — | г = 0. Перемещение твердого тела определяет смещение сетки на границе <<жидкость - твердое тело>>
Шш—ш5 | г = 0, ^ (71)
и скорость твердого тела определяет скорость жидкости на этой же границе
(72)
Численные алгоритмы решения задачи жидкостно0структурного взаимодействия в упругих пористых средах
Реализовать описанную выше конечно-элементную формулировку можно четырьмя различными способами:
• однонаправленной связью;
• частичной связью;
• монолитной связью;
• частично-монолитной связью упругой и жидкой компонент. Алгоритм частичной связи упругой и жидкой компонент
В данном алгоритме жидкая компонента определяется после задания начальных и граничных условий (рис. 7). Условие 7 2 означает, что на границе раздела Г скорости должны быть равными. Если это первая FSI-итерация по времени, скорость упругого тела берется с последнего временного шага, в противном случае - с предыдущей итерации. Далее находятся остатки и якобианы. Это необходимо для выполнения условия нелинейности несжимаемости жидкости. Для первой итерации Ньютона оценка берется с предыдущего шага по времени. Норма невязки 1°° рассчитывается для определения сходимости и, при необходимости, система решается повторно для обновления решения. Здесь для сходимости необходимо 1° < е — 1 4. Если решение не сходится, система уравнений Стокса решается с помощью прямого линейно-алгебраического решателя. Действия повторяются, пока одно из условий сходимости не будет удовлетворено.
Рис. 7. Блок-схема жидкостно-структурного взаимодействия с частичной связью упругой и жидкой компонент Fig. 7. Block diagram of the liquid-structural interaction with the partial coupling of the elastic and liquid components
Далее происходит переход к решению системы уравнений Ламе, которое также начинается с применения граничных условий. Условия для упругой компоненты на границе раздела «жидкость-твердое тело» соблюдаются в слабой форме и их необходимо добавить в остатки. Остатки и якобианы упругой компоненты получаются с границы Г из уравнения ( 5 1 ). Так как используется нелинейная деформация, то необходимо использовать итерации Ньютона для достижения сходимости, как для жидкой компоненты. Условие для этого базируется на разности перемещений упругого тела на границе Г между текущей и предыдущей FSI-итерацией. Для первой общей FSI-итерации решение предыдущего временного шага используется в нахождении решения предыдущей FSI-итерации. Далее берется норма разности Ь2 по границе. Если Ь2 > е — 8 , то решение считается не сходящимся и необходимы дополнительные FSI-итерации. Если необходимо, то можно воспользоваться методом релаксации для обеспечения сходимости.
Далее, следуя алгоритму, происходит определение перемещения расчетной сетки. У движения сетки есть строгие условия на поверхности , что перемещение сетки и упругого тела здесь равны (71). В отличие от жидкости и упругого тела, движение сетки основано на системе линейных уравнений и это означает, что данное решение сходится на первой итерации каждой общей итерации. Далее процесс повторяется, пока не будет достигнуто определенное условие окончания расчетов.
Алгоритм монолитной связи упругой и жидкой компонент
Отличительной чертой данного подхода является то, что граничные условия здесь определены глобально (рис. 8). Также, как можно видеть на рисунке 8, нахождение неизвестных для жидкости, упругой части и расчетной сетки происходит последовательно и, если необходимо удовлетворить условие сходимости, произвести итерации Ньютона. И уже после достижения сходимости можно перейти на следующий временной слой.
Алгоритм частично-монолитной связи упругой и жидкой компонент
Используемый нами алгоритм частично-монолитной связи упругой и жидкой компонент получается из двух вышеописанных. Граничные условия так же, как и в алгоритме с монолитной связью, определяются глобально. Этот алгоритм решает уравнения Стокса и Ламе только монолитно, но после нахождения решения и проверки на сходимость, как в алгоритме частичной
связью жидкости и упругого тела. Эти два процесса повторяются в общих FSI-итерациях, пока параметры общей сходимости не будут достигнуты (рис. 9).
Рис. 8. Блок-схема жидкостно-структурного взаимодействия с монолитной связью
упругой и жидкой компонент Fig. 8. Block diagram of the liquid-structural interaction with the monolithic relation of the elastic and liquid components
Рис. 9. Блок-схема жидкостно-структурного взаимодействия с частично-монолитной связью упругой и жидкой компонент Fig. 9. Block diagram of the liquid-structural interaction with a partially-monolithic relation of the elastic and liquid components
Апробация разработанных алгоритмов решения математических моделей совместного движения неоднородных жидкостей
Описанные выше алгоритмы использовались при численном моделировании процесса вытеснения нефти водой в упругом поровом пространстве.
На рисунке 10 и 11 представлены графики, отражающие изменение позиции границы раздела жидкостей с течением времени для различных коэффициентов упругости Ламе без учета поверхностного натяжения.
! ! ! : : : j 1 : j 'i ft : f it
............. .............
; ; .............:.............:............ j
: ____________i____________1____________ _____________ : : _____________!...„.....J____________ // / /У //.....
; ; i /У ,............
: : .............1.............1.............j............ I ! : .............1............. Л Л f
: : : ____________i____________;_____________•____________ _____________
I У ......... \ .............
..........i....... .............
j .............\............. ! .............
-л - 15 -1 = 50 к = 2D0
О s 10 15 2D 25 30 35 Л0 45 30
Игле
Рис. 10. График изменения позиции свободной границы с течением времени
для различных значений X Fig. 10. The graph of changing the position of the free boundary with the passage of time for different values of X
Рисунок 11. График изменения позиции свободной границы с течением времени
для различных значений X Figure 11. The graph of changing the position of the free boundary with the passage of time
for different values of X
Разработанный алгоритм был протестирован для геометрии порового пространства в виде несвязных элементов упругого скелета (рис. 12 и рис. 13).
Рис.12. Вытеснение нефти водой в упругом пористом скелете в виде несвязных элементов Fig. 12. The displacement of oil by water in an elastic porous skeleton in the form
of disconnected elements
Рис. 13. Обтекание водой упругого элемента пористой среды Fig. 13. Water flow around an elastic element of a porous medium
Заключение
В настоящей работе разработаны новые алгоритмы моделирования процесса вытеснения одной жидкости другой с учетом взяимодействия с упругим телом. Используя полученные алгоритмы, разработано научно-исследовательское программное обеспечение, позволяющее решать совместно уравнения Стокса для жидкостей и Ламе для упругого скелета в поровом пространстве различной геометрией: капилляры и несвязные элементы упругого скелета. При
этом сами уравнения могут быть как стационарными, так и нестационарными, и описывать как сжимаемую, так и несжимаемую жидкость.
Список литературы References
1. Batchelor G. K. 1967. An Introduction to Fluid Dynamics, Cambridge University Press.
2. Meirmanov A. 2014. Mathematical models for poroelastic flows, Atlantis Press, Paris.
3. Dahong Chen. 1997. "Numerische Simulation von Strömungsvorgängen mit der "Arbitrary Lagrangian-Eulerian method" (ALE-Methode)" Aachen: Mainz, XVII, 139 S. ISBN 3-89653-610-9.
4. Hirt, C.W., Amsden, A.A., and Cook, J.L. 1974. An Arbitrary Lagrangian-Eulerian Computing Method for all Flow Speeds, J. Comp. Phys., 14: 227.
5. Gong, K., Shao, S., Liu, H. et al. 2016. (2 more authors) . Two-phase SPH Simulation of Fluid-Structure Interactions. Journal of Fluids and Structures, 65: 155-179.
6. Галлагер Р. 1984. Метод конечных элементов. Основы: Пер. с англ. — М.: Мир.
Gallager R. 1984. Metod konechnykh elementov. Osnovy: Per. s angl. — M.: Mir.
7. Деклу Ж. 1976. Метод конечных элементов: Пер. с франц. М.: Мир, 1976.
Deklu Zh. 1976. Metod konechnykh elementov: Per. s frants. M.: Mir.
8. Зенкевич О. 1975. Метод конечных элементов в технике. М.: Мир.
Zenkevich O. 1975. Metod konechnykh elementov v tekhnike. M.: Mir.
9. Зенкевич О. 1986. Морган К. Конечные элементы и аппроксимация: Пер. с англ. М.: Мир.
Zenkevich O., Morgan K. 1986. Konechnye elementy i approksimatsiya: Per. s angl. — M.: Mir.
10. Сегерлинд Л. 1979. Применение метода конечных элементов — М.: Мир: 392.
Segerlind L. 1979. Primenenie metoda konechnykh elementov. M.: Mir: 392.
