CC BY
44
УДК 519.6
ЧИСЛЕННОЕ ИССЛЕДОВАНИЕ ПРОЦЕССОВ ВХОДА И ПОГРУЖЕНИЯ ПЛОСКОГО ЦИЛИНДРА В НЕСЖИМАЕМУЮ ЖИДКОСТЬ МЕТОДОМ СГЛАЖЕННЫХ ЧАСТИЦ
Р. С. Макарчук
NUMERICAL SIMULATION OF ENTRY AND SUBMERGING OF A 2D CYLINDER INTO INCOMPRESSIBLE LIQUID USING SMOOTHED PARTICLE HYDRODYNAMICS
R. S. Makarchuk
В статье исследуются процессы входа, погружения и последующего всплытия плоского кругового цилиндра в бассейне, заполненном однородной несжимаемой жидкостью, методом сглаженных частиц. Особое внимание уделяется вычислению гидродинамических сил, действующих на цилиндр, и нагрузок, испытываемых твердыми вертикальными стенками бассейна.
The article is devoted to numerical simulation of processes of a 2D circular cylinder entry into an incompressible liquid filling the basin, its subsequent submergence and exit with Smoothed Particle Hydrodynamics. Special attention is paid to computation of hydrodynamic forces on the cylinder and loads on the vertical walls of the basin.
Ключевые слова: метод сглаженных частиц, свободная граница, несжимаемая жидкость, гидродинамические нагрузки, погруженное в жидкость тело, схема расщепления
Keywords: SPH, free surface, incompressible fluid flows, hydrodynamic loads, immersed body, split-step scheme.
Изучение особенностей течений жидкости со свободными границами является актуальной задачей и имеет большое прикладное значение, поскольку такого рода течения встречаются повсеместно. Достаточно привести в качестве примеров такие сложные физические процессы, как накат волны на наклонный берег, взаимодействие волн с береговыми и донными сооружениями, погруженными в жидкость телами: задачи глиссирования, посадки гидросамолетов на поверхность водоемов и т. д. Исследованию процессов взаимодействия твердых и упругих тел с жидкостью посвящены работы многих отечественных и зарубежных ученых: Л. И. Седова, Г. В. Логвиновича, В. В. Пухна-чева, А. А. Коробкина, А. Г. Терентьева, Э. И. Григо-люка, А. Г. Горшкова, М. В. Норкина, В. И. Юдовича, Г. Г. Шахверди, H. Wagner, R. Zhao, O. Faltinsen, M. Greenhow, W. M. Lin, X. Zhu и др. Практический интерес в таких задачах представляют как кинематические характеристики течений - поле скоростей, положение свободной границы, так и динамические - поле давления, величины гидродинамических нагрузок на твердые границы области течения или погруженные тела. Использование классических численных методов (как эйлерового, так и лагранжевого типа) для проведения вычислительных экспериментов связано с определенными трудностями. Основная проблема первых состоит в сложности постановки граничных условий в связи с заранее неизвестным положением свободной границы. Недостаток вторых - аварийный останов программы при пересечениях ребер ячеек сетки, возникающих в случае больших деформаций области расчета. В связи с этим в последнее время все большую популярность стали набирать бессеточные лагранжевые методы, одним из наиболее известных представителей которых является метод сглаженных частиц [1, с. 1013
- 1024; 2, с. 375 - 389].
Существует множество работ, в которых представлены кинематические картины течений, полученные в ходе численного моделирования с применением метода сглаженных частиц, и лишь немногие из них посвящены вычислению поля давления, поскольку именно получение динамических характеристик является одной из основных слабых сторон метода. В настоящей работе особое внимание уделяется расчету поля давления и последующему вычислению гидродинамических сил, действующих на погруженное в жидкость тело, и нагрузок на твердые границы области течения. Инструментом проведения вычислительного эксперимента служит метод сглаженных частиц [3, с. 8417 - 8436], использующий модель несжимаемой жидкости и схему расщепления по физическим факторам [4, с. 745 - 762] для интегрирования по времени закона движения сплошной среды - системы уравнений Навье-Стокса. Использование модели несжимаемой жидкости для определения пиковых нагрузок на затупленное тело (круговой цилиндр в настоящей статье) в момент его входа в жидкость, как отмечается в работе [5, с. 159 - 185], не позволяет получить достоверные результаты, однако дальнейший процесс погружения описывается весьма хорошо.
Рассматриваемая в статье задача формулируется следующим образом. В бассейн длиной 1 м, заполненный жидкостью с плотностью р = 1000 кг/м3 и коэффициентом кинематической вязкости V = 10-6 м2/с до уровня 0,3 м, падает абсолютно твердый плоский круговой цилиндр плотностью рь < р с радиусом основания г = 0,055 м. Высота падения цилиндра определяется положением его центра масс хс=(0,5; 0,855) м в начальный момент времени 1 = 0 с. Задача обладает свойством симметрии относительно прямой, проходящей через центр масс цилиндра ортогонально покоящейся сво-
бодной поверхности жидкости в начальный момент времени. Симметричность постановки задачи означает, что моменты сил, действующих на цилиндр, равны нулю. Тем не менее задача решается в полной двумерной постановке.
Движение жидкости (занимающей область О) в бассейне подчиняется следующему векторному уравнению:
сій 1
— =-------Ур +уАа + f,
С р
(3)
ется как р (ґ) = |х є Rс |р(х,., х,) < к|, р - некоторая
метрика в Яё. Функции формы строятся по данному множеству частиц, которое в разные моменты времени в общем случае различно. Используя множество ближайших соседей в качестве набора узлов интегрирования, приходим к дискретному представлению формулы (5):
(1) /(х) = | f(XX)W(х - х, ЩСВ *£/(х,. ) W(х - х,, к) АО,
где и(х, 1) - вектор скорости, х(1;) - радиус-векторы точек области О, р(х, 1) - давление, р - массовая плотность жидкости, V - коэффициент кинематической вязкости, f - вектор плотности массовых сил. Помимо уравнения (1) для несжимаемых сред поле скоростей должно так же соответствовать условию соленоидаль-ности:
V и = 0. (2)
Для решения системы уравнений (1)-(2) необходимо задать начальные и граничные условия. Пусть Г = Г1 иГ2 иГ3 иГ4 - граница области О, где Г1 -свободная граница (поверхность жидкости в бассейне), Г2и Г3 - вертикальные стенки бассейна, Г4 - дно бассейна. Тогда граничные и начальные условия записываются следующим образом:
|Р„ (х^) 1хеГ1 = *) • П - РатмП
\ и(Х’ -) 1хеГ2 иГ3 иГ4 =
и(х> 0|,=0 = Ио(х)- (4)
Здесь п - вектор внешней нормали к свободной поверхности, Т(х, 1) - тензор вязких напряжений, Рп(х, 1)
- вектор напряжений на свободной поверхности в направлении внешней нормали, ратм - атмосферное давление.
К системе уравнений (1)-(2) с граничными и начальными условиями (3)-(4) для перемещения малых объемов жидкости (частиц) необходимо добавить следующее уравнение и начальные условия:
ёх
~Г = и.
Л
х(-) I=0 = хо •
Систему (1)-(2) будем решать численно методом сглаженных частиц. Основу метода составляет формула усреднения функции по Стеклову:
I(х) = | I(X)№(х - X, Н)ёО, (5)
О
где х, X е О с Rrf , d - размерность задачи. Здесь и далее предполагается зависимость области интегрирования В от переменной X. Функция W называется функцией ядра, кЬ - радиус ее носителя, Ь называется радиусом сглаживания или сглаживающей длиной, а значение коэффициента к зависит от конкретного вида функции.
Проведение вычислительного эксперимента требует перехода к дискретным аналогам уравнений динамики жидкости (1)-(2). С этой целью вся область расчета представляется набором расчетных узлов, между которыми нет жестких топологических связей, характерных, например, для конечно-элементных сеток. Вместо этого используется понятие ближайших соседей. Множество ближайших соседей 1-й частицы с центром х1 определя-
,=1
где АО. - связанный с )-й частицей объем, п - число
3 ^ у
частиц в области В. Ввиду отсутствия связной сетки, в методе сглаженных частиц задача определения объема частицы не может быть решена геометрически, подобно тому как это делается в бессеточных методах, основанных на слабой форме уравнений движения (МРЕМ, КЕМ и др.). Вместо этого объем частиц принято вычислять по формуле АО. = т. / р , где массы mJ задаются в качестве начального условия и сохраняют постоянные значения в течение всего времени расчета.
Функция ядра W является одним из ключевых понятий метода сглаженных частиц, и от ее выбора напрямую зависят результаты численных расчетов. Функции ядра, получившие наибольшее распространение, при расчете задач со свободными границами часто являются причиной образования групп (кластеров) частиц по всей области расчета, что оказывает серьезное влияние как на вычисляемые кинематические, так и динамические характеристики изучаемых процессов. Причина заключается в немонотонном поведении первых производных функции ядра, значения которых стремятся к нулю по мере приближения к центру области носителя, что приводит к ослабеванию сил отталкивания (градиент давления) между парами частиц по мере их сближения, обеспечивая благоприятные условия для их кластеризации. Решением этой проблемы может стать использование для аппроксимации градиента давления такой функции ядра, поведение первой производной которой на всей области ее носителя будет монотонным, что, соответственно, обеспечит монотонное возрастание отталкивающих сил между парами частиц по мере их сближения. Такой подход был предложен в работе [6]. В связи с тем, что в настоящей работе для проведения численных расчетов используется сплайн 4-й степени с радиусом носителя, равным 5/2Ь, для вычисления давления в работе предлагается функция ядра подобного вида, нормализованная для плоского случая:
W (д) =
60
931жп
(5/2-д)4 -(3/2-д)4, 0 < д < 3/2; (5/2-д)4, 3/2< д < 5/2;
0, д > 5/2,
где д = х - х /к
Несмотря на то, что вычисление градиента давления на основе предлагаемой функции ядра позволяет регу-ляризовать расчетную "сетку" внутри области течения, во многих задачах возникают все же проблемы, связанные с близким расположением частиц свободной поверхности к внутренним частицам среды. Такое расположение обусловливается недостаточным количеством ближайших соседей с внешней, по отношению к жидкости, стороны свободной границы, и сказывается на
точности вычисления характеристик течения для этих частиц, что приводит к осцилляциям давлений в приграничных областях. С целью избегания подобных эффектов предлагается алгоритм корректировки свободной границы, заключающийся в добавлении в правую часть уравнений движения дополнительных сил отталкивания, действующих на частицы свободной поверхности. В работе предлагается использование потенциала упругих шаров:
1/2К(г -г0)2, г < г0;
0, г > г0
где К - коэффициент жесткости взаимодействующих упругих шаров (в нашем случае частиц жидкости), г -расстояние между ними, г0 - радиус действия потенциала. Предложенный подход позволяет стабилизировать вычисления поля давления и получить достоверные кривые гидродинамических нагрузок. Результаты численного решения ряда модельных задач методом сглаженных частиц с использованием вышеупомянутых модификаций приведены в работе [7, с. 447 - 464].
Для определения перемещения цилиндра используется второй закон Ньютона. Внешние силы (за исключением силы тяжести), действующие на цилиндр, определяются по формуле:
FJ = -| рпйТ + | Т • пй?Г ,
г г
где первый интеграл в правой части - поверхностные силы давления, второй интеграл - силы вязкого трения, Г - поверхность цилиндра, п - вектор внешней нормали к поверхности цилиндра, р - давление, Т - тензор вязких напряжений. Интегралы вычисляются по формулам прямоугольников.
В вычислительных экспериментах принимают участие цилиндры плотности рь = 250 кг/м3,
рь = 500 кг/м3 и рь = 750 кг/м3. Общее количество частиц, использованных в численных расчетах, составило 14515, из которых 11940 частиц служили для представления заполняющей бассейн жидкости, 439 граничных и 1796 виртуальных частиц применялись для постановки условий на неподвижных твердых границах бассейна, 68 граничных и 272 виртуальные частицы моделировали поверхность цилиндра. Полученные методом сглаженных частиц результаты вычислительных экспериментов сравниваются с результатами численных расчетов методом С1Р, представленными в работе [8], и данными лабораторных экспериментов [9]. На рисунке 1 приведены картины течений при входе и погружении в жидкость цилиндра плотности рь =500 кг/м3 в момент
времени 1 = 0,03 с, на рисунке 2 - в момент времени 1 =
0,12 с. Момент касания цилиндром свободной поверхности жидкости в бассейне берется за начальный момент времени 1 = 0 с.
|||Щ11Ш^ в) I
Рис. 1. Картины течения в момент времени 1=0.03 с: а) результаты автора,
б) результаты метода С1Р, в) лабораторные снимки
Ли
а)
...... . ■ . V ........
........ - .. .
‘ - % 1 '' ж
кжуйШг' - Д
в)
Рис. 2. Картины течения в момент времени 1=0.12 с: а) результаты автора, б) результаты метода С1Р, в) лабораторные снимки
Рис. 3. Картины течения в момент времени 1=0.11 с для цилиндра плотностью рь = 1000 кг/м : а) результаты метода С1Р, б) лабораторные снимки
Из рисунка можно увидеть, что метод сглаженных частиц достаточно адекватно описывает процесс погружения тела в жидкость, однако, в отличие от результатов метода С1Р и снимков экспериментов, струи, возникающие по обе стороны от цилиндра и образующие стенки каверны, имеют склонность к направлению в сторону цилиндра. Тем не менее подобную ситуацию можно также наблюдать и в численных расчетах методом С1Р, но уже в случае цилиндра плотностью рь = 1000 кг/м3 (рис. 3).
На рисунке 4 представлен график зависимости глубины погружения цилиндра в жидкость от времени. Здесь глубина Н есть разность начального невозмущенного положения свободной поверхности жидкости в бассейне (у = 0.3 м) и вертикальной координаты нижней точки поверхности цилиндра. Пунктирная кривая - результат, полученный методом сглаженных частиц авторами статьи, непрерывная кривая - методом С1Р в работе [8], квадратами на графике отмечены результаты эксперимента [9]. Можно отметить определенное расхождение в результатах, полученных различными методами, увеличивающееся со временем, которое можно объяснить тем, что в математической модели не учитывается сопротивление воздуха при падении цилиндра.
цилиндра от времени
На рисунке 5 приведены хронограммы гидродинамических нагрузок (а) и уровень жидкости (б) на правой стенке бассейна: кривая 1 соответствует цилиндру плотностью ръ = 250 кг/м3 , кривая 2 - цилиндру плотностью рь = 500 кг/м3 и кривая 3 - цилиндру плотностью ръ = 750 кг/м3. Максимальные значения нагрузок на стенки бассейна (на рисунках обозначенные max1 и max2) соответствуют моментам максимального запле-ска на них волн, а уровень заплеска тем выше, чем больше масса цилиндра, погружение которого способствовало формированию этих волн.
Рис. 5. Хронограммы гидродинамических нагрузок (а) и уровень жидкости (б) на правой стенке бассейна
Рассмотрим процессы погружения и всплытия цилиндров в бассейне с жидкостью более подробно. За время падения цилиндр плотностью ръ = 250 кг/м3 не
приобретает импульса, достаточного для его полного погружения в жидкость. В процессе погружения цилиндра образуется каверна, ограниченная с обеих сторон струями жидкости. Из этих струй впоследствии формируются волны, расходящиеся от цилиндра в направлении стенок бассейна. На рисунке 5 (б) этим волнам соответствуют первые максимумы (обозначенные шах^. Одновременно с этим цилиндр, достигнув своего предельного нижнего положения, начинает подниматься. После достижения цилиндром наивысшей точки подъема (рис. 6, а) цилиндр повторно погружается в жидкость, в результате чего происходит формирование вторых волн, движущихся в направлении вертикальных стенок бассейна. Этим волнам соответствуют вторые максимумы на рисунке 5 (б), обозначенные шах2.
В случае цилиндра плотностью ръ = 500 кг/м3 приобретаемого им импульса достаточно для его полного погружения в жидкость - каверна, образующаяся в процессе погружения, схлопывается над цилиндром, формируя небольшую кумулятивную струю. Схлопы-вание каверны происходит после достижения цилиндром предельного нижнего положения на первом цикле погружения. Как и в предыдущем случае, из образовавшихся в процессе входа цилиндра в жидкость струй происходит формирование первых волн, которые движутся в сторону стенок, ограничивающих бассейн. В процессе последующего подъема по обе стороны от цилиндра происходит формирование вторых волн. Свой вклад в это вносит так же жидкость, которая поднимается вместе с цилиндром и затем стекает с его поверхности. Повторное погружение цилиндра завершает этот процесс.
Механизм формирования вторых волн в случае входа в жидкость цилиндра плотностью ръ = 750 кг/м3 отличен от двух ранее рассмотренных случаев. Приобре
тенного цилиндром импульса достаточно для его достаточно глубокого погружения в жидкость. За это время успевают произойти процессы распада струй, образующихся при входе цилиндра в жидкость, и схлопы-вания каверны над цилиндром, который в этот момент еще продолжает погружение. В результате распада струй образуются первые волны, движущиеся в направлении стенок бассейна. Схлопывание каверны является причиной образования кумулятивной струи, которая, распадаясь, формирует вторые волны. В ранее рассмотренных случаях вторые волны формировались в процессе повторного погружения цилиндра в жидкость.
На рисунке 6 представлены зависимости от времени положения центров масс цилиндров (а) и гидродинамических сил, действующих на цилиндры (б). Следует обратить внимание на то, что в отличие от сил, действующих на цилиндры с плотностью рь =250 кг/м3 и ръ =500 кг/м3, на интервале примерно от 0,2 с до 0,8 с гидродинамическую силу, действующую на цилиндр плотностью рь = 750 кг/м3, можно считать приблизительно постоянной. Это связано с тем, что цилиндр на данном промежутке времени максимально погружен в жидкость, однако процесс схлопывания каверны над цилиндром не дает ему возможности начать всплытие. Гидродинамические силы достигают своих максимальных значений в моменты максимального погружения цилиндров в жидкость (помимо момента удара цилиндра о поверхность жидкости), минимальных значений -в моменты наибольшего их подъема.
Результаты вычислительных экспериментов демонстрируют эффективность применения метода при решении задач, где основной интерес представляют не пиковые нагрузки на погружающиеся тела, а процессы волнообразования и гидродинамические нагрузки на конструкции, расположенные вблизи мест их входа в жидкость.
Рис. 6. Зависимость от времени положения центров масс цилиндров (а) и гидродинамической силы, действующей на цилиндры (б)
Литература
1. Lucy, L. B. A numerical approach to the testing of the fission hypothesis / L. B. Lucy // Astron. J. - 1977. - 82 (12).
2. Gingold, R. A. Smoothed particle hydrodynamics: theory and application to non-spherical stars / R. A. Gingold, J. J. Monaghan // Mon. Not. R. Astr. Soc. - 1977. - 181.
3. Lee, E.-S. Comparisons of weakly compressible and truly incompressible algorithms for the SPH mesh free particle method / E.-S. Lee, C. Moulinec, R. Xu, D. Violeau, D. Laurence, P. Stansby // J. Comput. Phys. - 2008. - 227.
4. Chorin, A. J. Numerical solution of the Navier-Stokes equations / A. J. Chorin // Math. Comput. - 1968. - 22.
5. Korobkin, A. A. Initial stage of water impact / A. A. Korobkin, V. V. Pukhnachov // Ann. Rev. Fluid Mech. - 1988. -20.
6. Johnson, G. R. SPH for high velocity impact computations / G. R. Johnson, R. A. Stryk, S. R. Beissel // Comput. Method. Appl. M. - 1996. - 139.
7. Afanas'ev, K. E. Calculation of hydrodynamic loads at solid boundaries of the computation domain by the ISPH method in problems with free boundaries / K. E. Afanas'ev, R. S. Makarchuk // Russ. J. Numer. Anal. - M., 2011. - Vol. 26. -№ 5.
8. Zhu, X. Application of the CIP Method to Strongly Nonlinear Wave-Body Interaction Problems / X. Zhu // Doctoral thesis for the degree of doctor ingenior. - 2006. - Режим доступа: http://ntnu.diva-portal.org/smash/get/diva2:122822/ FULLTEXT01 (дата обращения: 10.08.2012).
9. Greenhow, M. Non-linear free surface effects: Experiments and theory / M. Greenhow, W. M. Lin // Rep. No. 83-19, Dept. of Ocean Engineering, MIT, Cambridge, MA. - 1983. - Режим доступа: http://www.dtic.mil/cgi-bin/ GetTRDoc?AD=ADA161079 (дата обращения: 10.08.2012).
Информация об авторе:
Макарчук Роман Сергеевич - старший преподаватель кафедры ЮНЕСКО по НИТ КемГУ, 8 (384-2) 58-32-89, e-mail: mak@kemsu.ru.
Makarchuk Roman Sergeevich - senior lecturer at the Department of UNESCO on NIT Kemerovo State University.
