CC BY
86УДК 532.5
12 2 2 3
А.С. Козелков ' , А.А. Куркин , Е.Н. Пелиновский '
МОДЕЛИРОВАНИЕ ПАДЕНИЯ ТЕЛА В ВОДУ В РАЗЛИЧНЫХ УСЛОВИЯХ
НА ОСНОВЕ ЧИСЛЕННОГО РЕШЕНИЯ УРАВНЕНИЙ НАВЬЕ-СТОКСА ПОЛНОСТЬЮ НЕЯВНЫМ МЕТОДОМ
ФГУП «РФЯЦ-ВНИИЭФ», г. Саров,
Нижегородский государственный технический университет им. Р.А. Алексеева,
Институт прикладной физики РАН
Цель: Представлены результаты численного исследования падения тела в воду в различных условиях на основе уравнений Навье-Стокса.
Метод: Для численного моделирования многофазных течений со свободной поверхностью используется модель односкоростного приближения на базе полностью неявной связи скорости и давления. Результаты: Приводятся основные формулы дискретизации уравнений и вид коэффициентов, а также результаты верификации модели на задаче образовании дорожки Кармана при движении тела и падении тела в воду. На базе предложенной модели проведено моделирование как вертикального вхождения, так и под различными углами. Исследованы размеры области возмущения и выявлены закономерности изменения параметров. Область применения: Моделирование падения тел в воду при различных условиях и исследование всех параметров падения.
Ключевые слова: уравнения Навье-Стокса, численное моделирование, неявный метод, жидкость, твердое тело, взаимодействие.
Введение
Вопросу исследования взаимодействия тел с жидкостью посвящены работы многих отечественных и зарубежных ученых [1-5, 10]. В них приводятся аналитические и численные решения, результаты экспериментальных и лабораторных исследований.
Взаимодействие тела с поверхностью жидкости и его дальнейшее погружение определяется многими факторами: скоростью и углом вхождения тела; массой и формой тела; физическим состоянием воздушной и жидкой сред в момент взаимодействия; свойствами тела и многими другими факторами, учет влияния которых существенно усложняет аналитическое и численное исследование [6-9].Сложность рассматриваемых явлений не позволяет в полной математической постановке получить ни аналитические, ни численные решения широкого круга прикладных и научных задач. Даже численное исследование этого процесса представляет большие трудности, возникающие уже на этапе построения численных схем. Эти трудности связанны с существенно трехмерным нестационарным процессом обтекания тела, его турбулентным характером, нелинейностью, а также разрывным характером движения.
В прикладных и научных исследованиях падения тел в воду основными изучаемыми свойствами являются силы, действующие на тело и возникающие волны. Теоретический анализ этих свойств основывается на теории идеальной жидкости без учета влияния воздуха [2]. Однако, даже при сделанных упрощениях, точный расчет действующих сил и конфигурация свободной поверхности практически не возможен. В связи с этим актуальным становится вопрос выявления закономерностей и описания деталей этого явления на простых численных экспериментах с целью выработки теоретических закономерностей. Современное состояние исследований в этом направлении отражено в [2-4, 10].
Процесс падения тела в воду можно условно разбить на три стадии: стадия соударения с образованием каверны на поверхности воды, стадия начала погружения с сопутствующим схлопыванием каверны и стадия непосредственно погружения с сопутствующим схло-пыванием подводного пузыря [1, 3, 8, 10]. Для изучения возникающих в результате падения
© Козелков А.С., Куркин А.А., Пелиновский Е.Н., 2015.
тела поверхностных волн необходим учет влияния всех стадий, что возможно лишь при непосредственном их моделировании. Для задач распространения поверхностных волн зачастую используют параметрическое приближение[7-9,11-14], не учитывающее две последние стадии и представляющее собой каверну, образованную вертикальным падением тела.Самым главным недостатком параметрического источника является то, что он не учитывает угол вхождения тела в воду.
В настоящей работе проведены численные эксперименты вхождения тела в воду под разными углами с целью анализа параметров образующейся каверны. Вхождение моделируется с разными скоростями, отражающими дозвуковой и сверхзвуковой вход под различными углами. Анализируются высоты возникающих при соударении волн, а также три стадии движения тела. Представлена используемая математическая модель и применяемые граничные условия. Приводятся основные формулы дискретизации уравнений и вид коэффициентов матрицы системы. Представлены результаты верификации предложенной модели на задачах, позволяющих оценить правильность расчета сил, действующих на тело, и правильность образования формы источника. Приведены и обобщены результаты исследований вхождения тела под различными углами.
Основные уравнения
Рассмотрим несжимаемые среды, состоящие из произвольного количества гетерогенных фаз, разделенных границей раздела. В рамках односкоростной модели смеси движение данной системы описывается консервативными уравнениями, включающими уравнения неразрывности, сохранения импульса и уравнение для объемных долей фаз [15, 16]:
V-u = 0
dt k
d dt
Здесь u - вектор общей скорости совокупности фаз k ; p(k) - плотность фазы k ; a0~k) - объемная доля фазы k и Xak = 1; p - давление; ju(k) - молекулярная вязкость фазы k; g -
k
ускорение свободного падения.
Основная сложность при численном решении системы (1) заключается в определении связи поля давления с полем скорости. Процедура согласования поля давления с полем скорости должна приводить к одновременному удовлетворению уравнений неразрывности и сохранения импульса. Наиболее распространенными являются методы типа SIMPLE, основанные на процедуре коррекции давления или принципе расщепления неизвестных [17].
Для получения SIMPLE-процедуры численного решения системы (1) опустим уравнение для объемных долей и массовые силы, рассмотрим ячейку P с гранями f = nb(P) и запишем систему уравнений (1) в дискретном виде X unfSf=0
f=nb( P)
Xa(k)p(k) ^^
2 a(* Уk ) u = -V • 2 К )У ) uu) + V • 2 (a(k ) У) Vu) -Vp + pg
k k k
-a( k Уk) +V{a{ k Уk) u ) = 0
(1)
n-1
-V = -2 2a(k Vfk )unflunfSf + 2 2a(k V/ )VunfSf
f=nb(P) k f=nb(P) k
(2)
Здесь п - временной слой; т - временной шаг; - площадь грани /, разделяющей кон-
- 2 pfSf + 2a(k]P(k)gV
f=nb(P) k
k
трольные объемы расчетной сетки (рис. 1); u f — величина скорости на грани (здесь и далее
индекс f означает принадлежность величины к грани); nb(P) — количество граней ячейки (в данном случае ячейки Р).
В случае расчета многофазных течений с помощью алгоритма SIMPLE уравнение для
переноса объемных долей решается для (n — 1) объемной доли a^n= 1 — £ аР^фаз k от-
k ф n
дельно от уравнений неразрывности и сохранения импульса.
Наиболее удобно использовать совмещенный алгоритм решения системы (2), основанный на полностью неявной формулировке. Для его реализации запишем систему (2) в виде
D
f
£ juf
f=nb(P)
(n n — 1 \ U— U )
\vp} -1 ) —(vf )
■Sf = о
V = — £ £a( k)pf) unf~lu nfSf
T f =nb(P) k
(k),,(k V, n v.— V «"V.^V rv(k) r>(k)f
(3)
+ £ £a(k)Mfk )Vu "fSf — £ pfSf + £a( k) p(k) g V
f=nb(P) k
f=nb( P)
Здесь «верхняя черта» означает, что величина на грани получена интерполяцией из соседних контрольных объемов. В уравнении неразрывности в системе (3) была использована поправка Рхи-Чоу [18], нивелирующая разницу приближений градиента давления в уравнениях неразрывности и сохранения импульса, а также связывающая поля скорости и давления при одновременном решении уравнений неразрывности и движения.
В алгебраической форме система имеет вид
£
je{ p,u,v, w}
aPp jp +
£
N=NB (P )
£
je{ p,u,v, w}
aNp N
= bp,
aNjN
(4)
= bp, i = {u,v,w}.
£ Фр + £ £
р,иу,м>} N = ЫЕ(Р) р,иу,
Суммирование по индексу ] для первого уравнения системы (4) - уравнения неразрывности дает коэффициенты общей матрицы системы для вычисления давления в контрольных объемах дискретной модели. Данные коэффициенты имеют вид
(^8/ )■
aPP =
aN =
S f • d pn
a
pp -
p
= — £
a
f=nb( p)
pp
f
aNu = (1 — )Sf, aN = (1 — )Sfy, aNW = (1 — Xf )Sf,
(5)
a
pu p
=£
XfSf
a
pv
4Sf
, apw = £
AfSf
*р = £
/=пЬ( р г /=пЬ( Рг /=пЬ(Р)
При записи данных коэффициентов использовался алгоритм неортогональной коррекции [19] для возможности правильного расчета на произвольных неструктурированных
сетках, а также формула для вычисления давления на грани р / с помощью линейной интерполяции по значениям в центрах ячеек:
Р/ =Л/РР +(! ~Л/ ) Рм • (6)
Для первого уравнения системы (4) правая часть имеет вид
bp = £ D f Vpf ■ S f — D f Vpf
f=nb( p)
Sf -
Sf • Sf S f 'd pn
1pn
(7)
k
Суммирование по индексу / для второго уравнения системы (4) - уравнения сохранения импульса - дает коэффициенты общей матрицы системы для вычисления компонент скорости
Лк) „ _ ^/-Б/ ^ V )Лк)
a^uu — avv — a
4N
4N
к
ff + min(0, Zaf )р(к) Sf) S f ' d PN к
(8)
Первое слагаемое выражений (8) относится к диффузионному слагаемому, и коэффициенты общей матрицы системы будут иметь вид
N = (1 Щ, < = (1 Щ, а^р = (1 Щ,
a
a
up _
P
— Z ff, app — Z ^fSf, app — Z
f/
(9)
f—nb(P) f—nb(P) f—nb(P)
Как и для коэффициентов (5), для записи использовался алгоритм неортогональной коррекции. Второе слагаемое выражений (8) - конвективная составляющая, которая аппроксимируются с помощью любой известной дифференциальной схемы, применимой на произвольных неструктурированных сетках [20-22]. Обычно используется противопотоковая схема (UpwindDifferences, UD). Также используются: противопоточная схема с линейной интерполяцией (LinearUpwindDifferences, LUD), схема QUICK, центрально-разностная схема (Cen-tralDifferences, CD), схемы семейства NVD (NormalizedVariableDiagram), гибридные схемы. Все перечисленные схемы применяются совместно с противопоточной схемой для увеличения монотонности.
Дискретизация нестационарного слагаемого может быть осуществлена по одной из известных неявных схем [23]: схеме Эйлера; схеме Адамса-Бэшфорта. Вклад диффузионного и конвективного слагаемых уравнений сохранения импульса относится и к диагональным коэффициентам общей матрицы системы, которые с учетом нестационарного слагаемого дискретизованного с помощью схемы Эйлера будут иметь вид
Лк)„(k) V
a
p
— - Z
N—NB( P)
a
P
— - Z
N—NB( P)
— - Z
aN +Z< )PPk к T
N + ZaP pP ) к T
Лк)„(к) V
a
(10)
а'р- = - Ъ аТ + К )РР
N=Ж(Р) к Т
Для второго уравнения системы (4) правая часть имеет вид
bu bp
bVp
— Z
f—nb( P)
—Z
f—nb( P)
—Z
f—nb( P)
Vu -Vf - (S f -Vv - Vf - (S f -Vw - Vf - (S f -
S f -S f
S f - d PN
S f -S f
S f - d PN
S f -S f
S f - d PN
d PN )
d PN )
d PN )
+ PP )u
V ^(к'Р(к)gxV,
+ Z
4 p )vpV ^«(к )Р(к) gyV,
(11)
+ ZaP PP) w
V ^
P — + Za т к
(к )р(к) g
Таким образом, совмещенная система линейных алгебраических уравнений полностью неявного алгоритма для моделирования многофазного потока имеет вид
aPpp aPpu aPpv apw ' pp
aPup aPuu aPuv aPw uP - Z N—NBP)
avPp aPuv aPv a? VP
aPwp aPwu aPwv (A p _ wP -
pp pu pv pw
UN UN UN UN
up uu uV uw
UN UN UN UN
aVp aUV aVV aVW
an an an an
wp wu an an
wv ww aN aN
pn u
N
V
N
w
N
bp
bP
Данная система записана для вычисления общей скорости и общего давления многофазного потока, однако может быть обобщена и на случай расчета, где каждая фаза обладает своей скоростью и физическими свойствами, такими как сжимаемость и турбулентность. Эти обобщения будут проведены в следующих работах.
Для моделирования границ раздела фаз, после решения системы (12), решается уравнение переноса объемных долей (третье уравнение системы (1)), которое решается для (п -1) объемной доли фаз. Его дискретизация методом конечных объемов осуществляется
по схеме полностью аналогичной той, что используется для уравнения сохранения импульса. Для аппроксимации конвективного слагаемого уравнения переноса объемных долей используется схема М-С1СБАМ [24], относящаяся к классу сжимающих схем высокого разрешения и обеспечивающая сохранение минимально возможной толщины границы раздела сред, а также сохранение формы распределения объемных долей при параллельном переносе и вращении. Дискретизация нестационарного слагаемого также осуществляется по схемам Эйлера или Адамса-Бэшфорта.
В алгебраической форме данная система уравнений для к -й фазы имеет вид
а^аР + £ 4к)4к) = Ьр) (13)
N=ЛЕ(Р)
Коэффициенты матрицы неявного решения уравнения (13) имеют вид:
а(к > + Г
a™ = min( 0, uf—1) Sf ) a(k) =— £ aN
' N=NB( p)
b (k) = u(n-1) S (a(k) (k)) , V (k),n-1 bp =—uf Sf '(aMC— aUD) + ~ap , (14)
Т
V
Т
где а^кс, аЦкр - значение объемной доли на грани, найденной по схеме МСТСБАМ и по про-
тивопоточной схеме соответственно; ap),n 1 - значение объемной доли на прошлый шаг по
времени. Представленные слагаемые получены с помощью дискретизации методом отложенной коррекции.
Для численного решения итоговая система уравнений должна быть дополнена начальными и граничными условиями. На твердых стенках градиент давления и объемных
долей равен нулю т.е. — = 0, —— = 0, значение скорости равно нулю - u = 0, v = 0, w = 0.
dn dn
На входных границах задается фиксированное значение скорости и объемных долей, гради-
dp
ент давления равен нулю, т.е. — = 0. На выходных границах фиксируется статическое дав-
dn
ление, градиенты скорости и объемных долей равны нулю:
ды = 0, av = 0, ^w = 0, a = 0. dn dn dn dn
Движение тела моделируется с использованием метода погруженных границ (Immersed Boundary Method, IBM) [25]. Метод IBM предполагает выделение в расчетной области ячеек, полностью или частично занятых твердым телом, и введением для данных ячеек, дополнительной силы сопротивления, предложенной в [26]. Такой подход к моделированию подвижного твердого тела достаточно прост, не требует динамического перестроения расчетной сетки и дает хорошие результаты для практических задач, не требующих детального описания течений в пограничном слое вблизи поверхности твердого тела.
Правильное вычисление сил сопротивления, действующих на тело, является одной из самых сложных задач моделирования падения тела. Воспользуемся следующим алгоритмом. Определение силы сопротивления, действующей на тело со стороны жидкости, а также моменты этих сил находятся из следующих интегралов:
F = - | pndS ; Fr =j rtdS
s S
MP = - J p(r - r ) x ndS ; MT = - J (r - r ) x (rt )dS
(15)
б ^
где Гр сила сопротивления давления; ^ - сила сопротивления трения; т - напряжение трения;
п^ - нормальный и тангенциальный к поверхности тела векторы; р - давление в жидкости. Интегралы (15) вычисляются по всей поверхности твердого тела При численном расчете интеграл заменяется суммированием следующим образом:
Fp = - j pndS = -YpnA ; F = \rtds = Yrts ;
S i=1 S i=1
N
MP = - j p(r - r*) x ndS = (r - r*) x niSi ;
(16)
i=1 N
MT=- j (r-r*) x (rt)dS = £ (rt- r*) x (rt)i St ;
S 2=1 i
Основная проблема - нахождение границы твердого тела для определения на ней давления, площади и нормали. Для решения этой проблемы воспользуемся алгоритмом построения изоповерхности. Для нахождения изоповерхности зададим граничный уровень концентрации в ячейке (pcrjt = 0.5. Вначале интерполяцией найдем значение концентрации в узлах расчетной модели. Далее проводим поиск ячеек, содержащих границу твердого тела. Условие содержания границы состоит в том, что ячейка содержит узлы со значением объемной доли твердого тела как менее (pcrjt, так и более. Для ребер ячейки, содержащей границу твердого тела, находим точки, в которых ф = фсг^. Данная операция делается в предположении линейного распределения ф вдоль ребра
.Pcnt-Pi
p = pl+(p2-ply
(17)
где^,Д - узлы рассматриваемого ребра; (рх,(р2 - значение концентраций; Р- точка пересечения.
По имеющимся N точкам строим N треугольников с общей точкой Рс, которая является геометрическим центром найденных точек:
1 N
p =1 у p . с Nу г
Для каждого треугольника находим вектор-площадь по выражениям:
(18)
Sx = 1 x 2
P -P
1 2, y 1 1, y
p -P
1 3, y 11, y
s
= 2
p -p
1 2, x p,x
p -p
1 3,x p,x
p -p
1 2, z 1 1,z p -p
1 3,z 1 1,z p -p
1 2,z p, z p -p
1 3,z p, z
• s = 1
; Sz = 2
p -p
1 2,x p, x
p -p
1 3, x 1 1, x
p -p
1 2, y 1 1, y
p -p
1 3, y 11, y
(19)
(20)
Находим вклад силы и момента сил от каждого треугольника:
/V: (Я /'):
=[(Рс; г.) (.V /')
где Рс,у - центр масс треугольника, ., ^ ., Б . - компоненты вектор-площади; р - давление в ячейке; г* - центр масс всего твердого тела.
В данном виде алгоритм определения результирующей силы и моментов сил пригоден
S
для использования на произвольной неструктурированной сетке, состоящей из многогранников произвольной формы. Представленная модель реализована в пакете программ ЛОГОС -российском программном продукте инженерного анализа, предназначенном для решения сопряженных трехмерных задач конвективного тепломассопереноса, аэродинамики и гидродинамики на параллельных ЭВМ [20, 21, 27]. Пакет программ ЛОГОС успешно прошел верификацию и показал достаточно хорошие результаты на серии различных гидродинамических задач, включая расчеты турбулентных и нестационарных течений [7, 21, 28].
Верификация модели
Представленные здесь верификационные задачи направлены на оценку правильности моделирования сил, действующих на тело в потоке жидкости, и формирования свободной поверхности около тела и вдали от него.
Обтекание цилиндра
Данная задача является одним из самых известных и классических тестов для верификации программ, моделирующих движение жидкости или газа. В ней рассматривается обтекание жидкостью кругового цилиндра. При числах Рейнольдса (Re) выше 30 за цилиндром образуется последовательность вихрей, вращающихся попеременно вправо и влево, известная как «вихревая дорожка Кармана». При достижении критических чисел Рейнольдса (при переходе течения из ламинарного режима в турбулентный) дорожка разваливается. Экспериментальные исследования данной задачи детально представлены в [29]. В эксперименте определяются коэффициент лобового сопротивления цилиндра (Cx) и число Струхаля (St), отражающее временной период колебания (или срыва) вихрей.
Для оценки действующих на цилиндр сил будем рассматривать обратную задачу. Пусть цилиндр движется в канале, заполненном неподвижной вязкой жидкостью. Движение цилиндра соответствует ламинарному режиму с числом Рейнольдса 350. За характерный размер берется диаметр цилиндра d=0,2 м. Движение цилиндра моделируется с помощью метода IBM, а действующие на него силы вычисляются с помощью описанного выше алгоритма. Для рассматриваемого диаметра цилиндра безразмерный коэффициент лобового сопротивления Сх, соответствующий взятому числу Re, равен 1,3.
Моделирование осуществляется на серии сгущающихся блочно-структурированных гексагональных сетках (табл. 1) в геометрической области, соответствующей параллелепипеду с размерами: X Е [0;10]; У € [0;1]; ZE [0;0Д]. Параметры расчетных сеток представлены в табл. 1.
Таблица 1
Расчетные сетки и характерные размеры ячеек
Номер сетки Число ячеек, млн Характерный размер ячейки, м
1 0,1 0,01
2 0,4 0,005
3 0,5 0,0044
4 1,6 0,0025
В расчетах использовалась схема второго порядка точности по времени. Для аппроксимации конвективных слагаемых в уравнении сохранения импульса использовалась проти-вопоточная схема. Шаг по времени выбирается переменным с условием соответствия числу Куранта, равному 0,9. Длина расчетной области позволяет с достаточной точностью определить установившиеся динамические характеристики.
На рис. 1 изображены мгновенные поля скорости в разные моменты времени. Пред-
ставленная картина соответствует классическому описанию отрыва пограничного слоя при обтекании тела с тупой кормовой частью [29]. В табл. 2 приведены результаты сравнения экспериментальных значений со значениями, полученными в результате расчетов.
Таблица 2
Результаты сравнения численных и экспериментальных значений
№ сетки С (расчет) С (эксперимент) Погрешность, % Число Струхаля (расчет) Число Струхаля (эксперимент) Погрешность, %
1 1,59 1,3 22,8 0,22 0,21 5,82011
2 1,40 1,3 8,45 0,22 0,21 5,82011
3 1,39 1,3 7,61 0,21 0,21 2,04082
4 1,37 1,3 5,61 0,21 0,21 2,04082
5 1,35 1,3 4,12 0,21 0,21 0,25063
О 5 10 15 20 х, м
Рис. 1. Поле скорости вокруг цилиндра в разные моменты времени
Результаты численных экспериментов показывают хорошее совпадение между собой и сходимость численного решения при измельчении сетки. Сходимость значений лобового сопротивления подтверждает правильность расчета сил сопротивления, действующих на цилиндр.
Падение параллелепипеда в воду
Верификацию представленного алгоритма на правильное описание образования каверны при падении тела в воду можно осуществить с помощью численного моделирования эксперимента, описанного в [2]. С помощью этой задачи можно также оценить правильность описания распространения волн на свободной поверхности на расстояние.
В этом эксперименте с высоты Н вдоль торцевой стенки бассейна в воду свободно падает прямоугольный параллелепипед с высотой Н1 и длиной I. Прямоугольный бассейн с
горизонтальным дном заполнен водой на высоту И < Н1. В натурном эксперименте длина бассейна составляла 4,3 м, а ширина 0,2 м. В невозмущенном состоянии вода находится в состоянии покоя. Расчетная модель представлена на рис. 2.
1 н1
н ? ? /1 ; ;
дг0=1б х„=31,5 д;0,м Рис. 2. Геометрия задачи
В [2] рассматривается несколько вариантов проведения эксперимента, отличающиеся друг от друга размерами тела, высотой, с которой оно падает, и глубиной воды. В безразмерной форме эти параметры выглядят так:
Н0 = Н/И, Н° = Щ/И, ¡0 = ¡¡И, х0 = хИ, р0 = (р-р)/р
где рх, р - плотность тела и плотность жидкости соответственно. Величины И, Н0, Н0,10, р0, х° меняются в зависимости от постановки.
В первом случае значения этих величин следующие: И = 8 см, . На расстоянии х° = 16 установлен мареограф, фиксирующий колебания водной поверхности. Во втором случае И = 4 см, Н0 = 4.75, Н° = 4.5, 10 = 1.15, р0 = 0.215, а мареограф установлен на расстоянии х° = 31.5.
Для моделирования падения параллелепипеда в воду была построена блочно-структурированная расчетная сетка, состоящая примерно из 1 млн ячеек в области с размерами 4,3х0,4х0,5 м. В месте падения параллелепипеда и установки датчика сетка более мелкая. На левую, правую и нижнюю стенки дискретной модели накладывалось граничное условие прилипания потока, сверху фиксировалось статическое давление.
В расчетах использовалась схема второго порядка точности по времени, для конвективных слагаемых в уравнении сохранения импульса использовалась противопоточная схема, для конвективных слагаемых в уравнении переноса объемных долей - схема М-СГСБАМ. Шаг по времени составлял 0,001 с.
На рис. 5 показано образование волн от падения параллелепипеда для первой постановки задачи. В момент времени I = 0,25 с наблюдается ранняя стадия развития каверны и брызговой струи. В момент I = 0,3 с струя продолжает развиваться. Через 0,1 с дно канала обнажилось, скорость формирования струи уменьшилась. В момент I = 7,5 с показана заключительная стадия схлопывания каверны и генерации брызговой струи. Следует отметить, что в работе [2] не указано, в какие моменты времени были сделаны фотографии. Авторами этой статьи были самостоятельно выбраны моменты времени для построения результатов численного моделирования.
а)
б)
Рис. 3. Образование волн от падения параллелепипеда для первой постановки задачи:
а - результаты численного моделирования на моменты времени 0,25; 0,3; 0,4 и 7,5 с;
б - экспериментальные данные
На рис. 4 приведены зависимости амплитуды колебаний уровня воды в бассейне от времени для второй и третьей постановки задачи соответственно. Цветом выделена кривая, построенная на основе данных, записанных неподвижным волномером.
Рис. 4. Амплитуда колебаний уровня воды в бассейне для второй и третьей постановок задачи соответственно:
А - амплитуда колебаний уровня воды; ^ - время
Сравнение мареографных записей показывает хорошее совпадение для первой приходящей волны. Вычисленная амплитуда вторая волны («отрицательной») также хорошо соответствует экспериментальным данным, однако наблюдается различие для дальнейших колебаний поверхности, которые имеют миллиметровую амплитуду. Это может быть связано как с качеством расчетной сетки, так и со способом записи численных данных в мареографной точке.
Моделирование возмущений при падении тел в воду в различных условиях
Рассмотрим задачу падения тел в воду при различных скоростях падения и с различными углами атаки. Для исследования рассмотрим модельную дискретную область размером
3 „
160x160x50 см , представляющую собой неструктурированную трехмерную (3Б) сетку, состоящую из усеченных многогранников произвольной формы (рис. 5). Расчетная сетка построена средствами препроцессора пакета программ ЛОГОС. В модельной области взята постоянная глубина относительно поверхности 10 м. Высота над нулевым уровнем поверхности воды составляет 40 м.
Верхняя часть рис. 5 справа отражает структуру сетки по вертикальной координате. Расстояние от поверхности воды (нулевой уровень) до максимального верхнего значения границы воздуха составляет 40 м, а до дна 10 м. Относительно нулевого уровня вверх для более детального разрешения основного выброса воды в воздух, было построено несколько сеточных слоев со сгущением.
Рис. 5. Модельная дискретная область (слева) и фрагмент расчетной сетки в месте падения (справа). Зоны сгущения сетки вблизи свободной поверхности и области движения тела (черный цвет)
Характерные размеры ячеек в воздухе: от 0 до 0,2 м - Дг = 10 см; от 0,2 м до 15,5 м -Дг = 25 см; от 15,5 м до 17 м - Дъ = 50 см; от 17 м до 21 м - Дг = 0,75 м; от 21 м до 27 м -Дг = 1,5 м; и от 27 м до 35 м - Дъ = 3 м. Размеры ячеек в воде по вертикали: от 0 до 0.2 м -Дг = 10 см, и от 0,2 м до 10 м - Дг = 25 см.
Для исследований будем рассматривать тело сферической формы диаметром 1 м движущееся со скоростью 10, 30, 60, 100, 150, 200 м/с под углами 0, 20, 30, 45, 60 и 80 град относительно поверхности воды для каждой из скоростей. Для корректного моделирования движения тела в воде проведено сгущение расчетной сетки в цилиндрической области до характерного размера ячеек 10 см так, чтобы на диаметр метеорита приходилось не менее 10 ячеек. Данная область расположена с наклоном, соответствующим углу падения тела (на рис. 5 20 град) (нижняя часть справа).
На рис. 6 - рис. 8 представлены результаты моделирования падения тела под углами 20, 45 и 80 град (тело летит слева направо).
При падении под углом 20 град образовавшаяся каверна имеет несимметричную структуру с выделяющимся всплеском по ходу падения тела (передняя волна) (рис. 6). За телом при малых скоростях (~10 м/с) поднятия поверхности воды (задняя волна) практически не наблюдается, тогда как по ходу падения высота волны достигает 2-4 м. Образовавшаяся каверна имеет клиновидную структуру и ее ширина около момента схлопывания составляет около 8 м. По мере возрастания скорости падения высота основного выброса волны увеличивается. Увеличивается и высота волны, образованной за телом, увеличивая при этом и размер каверны. По мере увеличения скорости амплитуда передней и задней волн отличаются в 4-5
раз. Передняя волна имеет «серповидную» форму, поэтому будет обрушаться как на внутреннюю сторону каверны, так и на внешнюю по ходу падения тела.
-15 -10 0 10 20 3» х, m -15 -10 0 10 20 30 х, ш -15 -10 0 10 20 30 х, m
е)
Рис. 6. Падение тела под углом 200со скоростью:
a - 10 м/с; б - 30 м/с; в - 60 м/с; г - 100 м/с; д - 150 м/с; е - 200 м/с
е)
Рис. 7. Падение тела под углом 450со скоростью:
а - 10 м/с; б - 30 м/с; в - 60 м/с; г - 100 м/с; д - 150 м/с; е - 200 м/с
Рис. 8. Падение тела под углом 800со скоростью
a - 10 м/с; б - 30 м/с; в - 60 м/с; г - 100 м/с; д - 150 м/с; е - 200 м/с
Изменение угла падения до 45 град уменьшает амплитуду переднее волны и увеличивает амплитуду задней, при этом размер полной каверны становится больше (рис. 7). При этом незначительное увеличение угла падения резко меняет форму передней волны, которая от «серповидной» формы переходит к вертикальному всплеску вверх и эта закономерность усиливается при увеличении скорости. Обрушение данной волны будет происходить по схеме падения вертикально вниз с уклоном на внутреннюю часть каверны, практически на захватывая внешней её стороны по ходу падения тела. Во внутреннюю часть каверны будет обрушаться и задняя волна. При таких режимах падения образование внутренней части каверны будет продолжаться даже при обрушении передней и задней волн и в итоге приведет к образованию большой воздушной полости в толще воды. По мере увеличения угла падения такая картина будет сохраняться.
Увеличение угла падения до близкого к вертикальному выравнивает амплитуды выбросов воды по ходу и за ходом падения тела (рис. 8). При этом при увеличении угла каверна начинает схлопывается раньше. При вертикальном падении каверна будет иметь симметричную структуру. В этом случае схлопывание каверны начинается гораздо раньше, чем при меньшем угле падения, а выбросы воды имеют гораздо меньшую амплитуду (может отличаться в разы), чем основной выброс воды при падении под углом.
Оценить закономерность изменения геометрии каверны можно по следующим параметрам (рис. 9): Япп — внутренний радиус каверны; Яш — внешний радиус каверны; Ны — внутреннее углубление каверны; Ноыг, — внешний подъем каверны.
а) б)
Рис. 9. Основные параметры каверны:
а - симметричная каверна; б - косая каверна, образованная падением под углом
В табл. 3 приведены значения этих параметров для различных скоростей и углов падения. Значения для скоростей падения 10-100 м/с представлены на момент времени 0,4 с, а для скоростей 150 и 200 м/с на момент 0,2 с. В эти моменты времени каверна из активной фазы роста начинает переходить в фазу обрушения. Для величины НоШ падения под углом
приводится максимальное значение одного из образовавшихся выбросов воды. Измерения между точками проводились в автоматическом режиме с использованием специальной функции вычисления расстояния постпроцессора пакета программ ЛОГОС (крайние точки выбирались в соответствии с рис. 9 справа). Отметим, что слишком точно определить величины этих параметров весьма проблематично из-за размытости границ, поэтому приводятся приблизительные (оценочные) данные.
Как видно из табл. 3, изменение параметров относительно вертикального падения наблюдается до и после 20 град. При этом внутренний радиус каверны и внутреннее
углубление каверны Нп до 20 градусов увеличиваются, а после 20 градусов начинают уменьшаться постепенно до параметров вертикального падения. С внешним радиусом каверны ЯоШ и внешним подъемом каверны Нои( изменения происходят с точностью до наоборот.
Изменение параметров при увеличении скорости падения происходит более резко, а самое существенное изменение наблюдается при небольших углах падения. Изменение как в сторону увеличения, так и уменьшения параметров происходит по линейному закону (рис. 10).
Таблица 3
Изменение параметров каверны при падении
Скорость падения, м/с Угол падения К, м R ,, м out7 H, м Hout,м
10 0 1.0 2.5 3.3 1.2
20 1.6 2.0 1.3 1.8
30 1.4 1.8 1.9 2.0
45 1.3 2.2 2.4 1.8
60 1.2 2.3 3.0 1.5
80 1.1 2.4 3.1 1.3
30 0 2.5 6.0 9.5 4.6
20 3.6 4.0 3.5 6.2
30 3.3 4.2 5.0 6.0
45 3.0 4.6 7.0 5.6
60 2.7 5.1 8.3 5.2
80 2.6 5.6 8.9 4.7
60 0 3.2 8.0 9.8 7.0
20 5.0 6.0 4.5 10.5
30 4.4 6.5 6.0 10.0
45 4.0 6.9 8.7 8.8
60 3.7 7.4 9.4 8.0
80 3.3 7.8 9.7 7.5
100 0 4.2 10.0 10.0 10.0
20 6.7 7.0 5.8 16.0
30 5.7 7.2 7.9 14.5
45 5.0 8.0 10.0 13.0
60 4.6 9.0 10.0 11.7
80 4.4 9.7 10.0 11.0
150 0 4.7 15.0 10.0 10.2
20 7.5 10.5 6.1 18.0
30 6.3 11.1 7.9 14.0
45 5.6 12.0 10.0 13.7
60 5.1 13.2 10.0 12.2
80 4.8 14.5 10.0 11.1
200 0 5.0 18.0 10.0 11.0
20 8.0 11.0 6.4 19.0
30 7.0 12.0 8.5 14.5
45 6.1 13.3 10.0 14.0
60 5.6 15.8 10.0 13.0
80 5.1 17.0 10.0 12.0
Рис. 10. Графики изменения параметров каверны от угла падения на различных скоростях:
1 - 200 м/с; 2 - 150 м/с; 3 - 100 м/с; 4 - 60 м/с; 5 - 30 м/с; 6 - 10 м/с
Изменение параметров при малых углах происходит по линейному закону. При увеличении угла тенденция сохраняется и приобретает такой же характер при существенном увеличении скорости. Хотя, как видно из графиков, закономерность изменения параметров примерно одинаковая на различных скоростях.
Более резкое изменение параметров на больших скоростях также указывает на то, что для таких режимов падения необходимо учитывать изменение термодинамических свойств жидкости и газа (в частности сжимаемость).
Выводы
В работе представлены результаты численного исследования падения тела в воду при различных углах входа и разных скоростях. Для численных экспериментов используется модель, основанная на уравнениях Навье-Стокса для многофазных течений со свободной поверхностью. Модель основана на односкоростном приближении и базируется на полностью неявной связи скорости и давления. Правильное описание моделью падения тел и их движения в воде показано при численном решении задачи образования дорожки Кармана и падении тела в воду.
Предложенная модель была использована для исследования размеров области возмущения и выявления закономерностей изменения параметров каверны. Показано, что изменение параметров каверны наиболее интенсивно происходит при малых углах падения и под-
чиняется линейному закону. Интенсивность изменения растет по мере увеличения скорости, а тенденция линейного изменения сохраняется.
Библиографический список
1. Лаврентьев, М.А. Проблемы гидродинамики и их математические модели / М.А. Лаврентьев, Б.В. Шабат. - М.: Главная редакция физ.-мат. лит-ры изд-ва «Наука», 1973.
2. Букреев, В.И. Гравитационные волны при падении тела на мелкую воду / В.И. Букреев, А.В.Гусев // Прикладная механика и техническая физика. 1996. Т. 37. №2. С. 90-98.
3. Gekle, S. Supersonic Air Flow due to Solid-Liquid Impact / S. Gekle [et al.] // Physical Review Letters, 0031-9007/10/104(2) /024501(4).
4. Aristoff, J.M. On falling spheres: the dynamics on water entry and descent along a flexible beam // Partial fulfillment of the requirements for the degree of doctor of philosophy at Massachusetts Institute of Technology. - USA: MIT. 2009.
5. Pierazzo, E. Validation of numerical codes for impact and explosion cratering: Impacts on strength-less and metal targets / E. Pierazzo [et al.] // Meteoritics& Planetary Science 43. 2008. № 12. Р. 1917-1938.
6. Комаров, А.А. Падение тела в резервуар с жидкостью и расчет возникающих при этом динамических нагрузок / А.А. Комаров, В.В. Казеннов // Вестник МГСУ. Сер. Гидравлика. Инженерная гидрология // Гидротехническое строительство. 2014. № 5. С. 135-143.
7. Козелков, А.С. Моделирование цунами космогенного происхождения в рамках уравнений Навье-Стокса с источниками различных типов / А.С. Козелков [и др.] // Изв. РАН. Сер. Механика жидкости и газа. 2015. №2. С. 142-150.
8. Козелков, А.С. Эффекты, сопровождающие вхождение метеорита в водную среду // Труды НГТУ им. Р.Е. Алексеева. 2014. № 3(105). С. 48-77.
9. Козелков, А.С. Цунами космогенного происхождения / А.С. Козелков, А.А. Куркин, Е.Н. Пе-линовский // Труды Нижегородского государственного технического университета им. Р.Е. Алексеева. 2014. № 2(104). С. 26-35.
10.Логвинович, Г.В. Гидродинамика течений со свободными границами / Г.В. Логвинович. -Киев: Наукова Думка, 1969.
11. Левин, Б.В. Физика цунами и родственных явлений в океане / Б.В. Левин, М.А. Носов. - М.: Янус-К, 2005.
12.Ward, S.N. Asteroid impact tsunami: A probabilistic hazard assessment / S.N. Ward, E.Asphaug // Icarus. 2000. V. 145. №1. Р. 64-78.
13.Kharif, C. Asteroid impact tsunamis / C. Kharif, E. Pelinovsky // C. R. Physique. 2005. № 6. Р.361-366.
14. Badescu, V. Dynamics and Coastal Effects of Tsunamis Generated by Asteroids Impacting the Black Seа / V. Badescu, D. Isvoranu // Pure Appl. Geophys. 2011. № 168. Р. 1813-1834.
15.Волков, К.Н. Течения газа с частицами / К.Н. Волков, В.Н. Емельянов. - М.: Физматлит, 2008.
16.Darwish, M A coupled finite volume solver for the solution of incompressible flows on unstructured grids / M. Darwish, I. Sraj, F. Moukalled // Journal of Computational Physics. 2009. V. 228. Р.180-201.
17.Ferziger, J. H. Computational methods fluid dynamics / J. H. Ferziger, M. Peric // Springer, 2001.
18.Rhie, C.M. A numerical study of the turbulent flow past an isolated airfoil with trailing edge separation / C.M. Rhie, W.L. Chow // AIAA. 1983. V. 21. Р. 1525-1532.
19. Jasak, H. Error Analysis and Estimation for the finite volume method with applications to fluid flows / H. Jasak. Thesis submitted for the degree of doctor. Department of Mechanical Engineering, Imperial College of Science, 1996.
20. Математические модели и алгоритмы для численного моделирования задач гидродинамики и аэродинамики: учеб. пособие / А.С. Козелков [и др.]; НГТУ им. Р.Е. Алексеева. - Нижний Новгород, 2014. - 166 c.
21.Козелков, А.С. Моделирование течений вязкой несжимаемой жидкости на неструктурированных сетках методом отсоединенных вихрей / А.С. Козелков [и др.] // Математическое моделирование. 2013. Т. 26. №8. С. 81-96.
22.Роуч, П. Вычислительная гидродинамика / П. Роуч. - М.: Мир, 1980.
23.Waclawczyk, T. Remarks on prediction of wave drag using VOF method with interface capturing approach / T. Waclawczyk, T. Koronowicz // Archives of civil and mechanical engineering, 2008. V. 8. Р. 5-14.
24.Mittal, R. Immersed boundary methods / R. Mittal, G. Iaccarino // Ann Rev Fluid Mech. 2005. V. 37. Р. 239-61.
25.Posa, A. Large-eddy simulations in mixed-flow pumps using an immersed-boundary method / A. Posa [et al.] // ComputFluids, 2011. V. 47(1). Р. 33-43.
26.Козелков, А.С. Реализация метода расчета вязкой несжимаемой жидкости с использованием многосеточного метода на основе алгоритма SIMPLE в пакете программ ЛОГОС / А.С. Козелков [и др.] // Журнал ВАНТ. Сер. Математическое моделирование физических процессов. 2013. Вып.. 4. С. 44-56.
27.Козелков, А.С. Минимальный базис задач для валидации методов численного моделирования турбулентных течений вязкой несжимаемой жидкости / А.С. Козелков [и др.] // Труды НГТУ им. Р.Е. Алексеева. 2014. № 4(104). С. 21-69.
28.Шлихтинг, Г. Теория пограничного слоя / Г. Шлихтинг. - М.: Наука, 1974.
Дата поступления в редакцию 15.07.2015
A.S. Kozelkov1'2, A.A. Kurkin2, E.N. Pelinovsky2'3,
MODELING OF A BODY FALLING INTO THE WATER UNDER DIFFERENT CONDITIONS BASED ON THE NUMERICAL SOLUTION OF THE NAVIER-STOKES
EQUATIONS FULLY IMPLICIT METHOD
FSUE «RFNC - VNIIEF»1, Nizhny Novgorod state technical university n.a. R.E. Alexeev , Institute of Applied Physics RAS3
Purpose: In this paper results of numerical investigation of a body falling into the water under different conditions based on the Navier-Stokes equations are discussed.
Method: For the numerical simulation of multiphase flows with a free surface is used a model of a single-speed approach based on fully implicit link of speed and pressure.
Results: The basic formulas and equations sample form of the coefficients, and the results of model verification on the task of education Karman motion of the body and the body fall into the water are presented. On the basis of the proposed model has been simulated as a vertical entry and entry from different angles. The dimensions of the perturbations are studied and regularities of changes of parameters are revealed.
Application domain: Modeling of bodies falling in water under various conditions and study all parameters of fall. Key words: Navier-Stokes equations, numerical simulation, implicit method, fluid, solid, interaction.
