CC BY
49
УДК 532.5 ББК 22.253
А.В. ЧЕЧНЕВ
ЧИСЛЕННОЕ ИССЛЕДОВАНИЕ НЕСТАЦИОНАРНЫХ ГИДРОДИНАМИЧЕСКИХ ЗАДАЧ С ТВЕРДЫМИ И СВОБОДНЫМИ ГРАНИЦАМИ МЕТОДОМ ГРАНИЧНЫХ ПОТЕНЦИАЛОВ
Ключевые слова: моделирование, гидродинамические процессы, численные методы.
Предложен численный метод решения нестационарных задач со свободными и твердыми границами. Он основан на представлении искомой функции на каждом шаге по времени в виде суммы потенциалов простого и двойного слоев. В результате из граничных условий получается система линейных интегральных уравнений Фредгольма 2-го рода, которая аппроксимируется затем СЛАУ путем замены входящих в нее интегралов по методу Н.М. Крылова - Н.Н. Боголюбова. Нормаль скорости точек свободных границ вычисляется по формуле, вытекающей из интегралов представленного потенциала поля скоростей путем его дифференцирования по направлениям нормалей свободных границ. Численные эксперименты показали, что предложенный метод является устойчивым по времени и позволяет получать достаточно точные численные решения.
A. CHECHNEV
NUMERICAL STUDY OF NON-STATIONARY HYDRODYNAMIC PROBLEMS WITH RIGID AND FREE BOUNDARIES BY METHOD OF BOUNDARY POTENTIALS
Key words: modeling, hydrodynamic processes, numerical methods. The article offers a numerical method of solving non-stationary problems with free and rigid boundaries. It is based on the representation of the unknown function to be found at each time step as the sum of the potentials of simple and double layers. The boundary conditions result in a system of linear integral Fredholm equations of the 2nd kind, which is afterwards approximated by the system of linear algebraic equations (SLAE) by replacing its integral constituents according to N.M. Krylov - N.N. Bogoliubov method. The normal rate of free boundary points is calculated according to the formula which is a consequence of the presented velocity field potential integration by its differentiation according to the directions of free boundaries normals. Numerical experiments showed that the proposed method is stable in time and allows obtaining quite accurate numerical solutions.
Как известно, всякая нестационарная гидродинамическая задача с твердыми и свободной границами в модели идеальной несжимаемой жидкости может быть сведена к последовательности смешанных задач для уравнения Лапласа [1-3, 5, 7]. Последние обычно решаются методом интегральных уравнений, позволяющим свести их к решению систем линейных алгебраических уравнения (СЛАУ). Именно от типа используемой системы интегральных уравнений зависит степень обусловленности аппроксимирующей ее СЛАУ [6], а значит, и качество получаемых численных решений. Поэтому при конструировании новых численных методов нужно стремиться к тому, чтобы решение соответствующей системы интегральных уравнений представляло собой корректно поставленную задачу, что является необходимым условием равномерной ограниченности числа обусловленности аппроксимирующей ее СЛАУ относительно времени t е [0, Т] при всех достаточно малых шагах по времени и пространству.
Ниже описывается численный метод, удовлетворяющий указанному условию. Он основан на нахождении потенциала поля скоростей ф в форме суммы потенциалов простого и двойного слоев, распределенных, соответственно, вдоль твердых и свободной границ:
т
ф(го,0 = -£ Iрк(г,01ф -го|(15г - |#(г,0<1ю(г), (1)
к=Т V) Г(0
где г0 е 0{(), В(?) - область, занимаемая жидкостью, в момент времени t е [0, Т]; Т- заданное время эксперимента; Тк, к = 1,т - твёрдые границы;
Г(0 - свободная граница области О^); <ю(г0,г) = 005ф202^ ; ф^- угол меж-
р - го\
ду векторам пг и г0 — г, п2 - орт внешней нормали к кривой Г(0 в точке гладкости г0 е Г(0, рк(г,£),к = 1,т, и q(z,t)- искомые функции.
В качестве примера рассчитано вертикальное движение кругового цилиндра, находившегося в начальный момент на некоторой глубине.
1. Постановка общей гидродинамической задачи с подвижными твердыми и свободной границами. Пусть в момент времени t = 0 покоящаяся жидкость занимает некоторую область 0(0), ограниченную горизонтальной прямой Г(0) и т замкнутыми кусочно-гладкими жордановыми контурами Тк (0), к = 1, т . При t > 0 последние движутся с заданными скоростями Ук ^), к = 1, т .
Требуется найти для всех t е [0, Т] положения границ Г(^ и Тк, к=1, т , и непрерывно-дифференциируемую в замкнутой области ), ограниченную кри-
выми Г(^ и Тк, к=1,т, и гармоническую ограниченную в области О^) функцию ф(г, 1), удовлетворяющую следующим граничным и начальным условиям:
|ф (г, 0=( Ук (ОД), при г е Тк , к = 1^; (2)
дпг
' ^ 1
■^ф(= ^Тф(2А2 - gУ, при г е Г(t); (3)
ф(г,0) = ф0 (г), при г е Г(0);
ёг
(4)
(5)
■ = 8га<ф(г, t), при г е Г^);
[г(0) = /0(1)+/^0(х), -те (00,60);
' t _
гк (т, t) = гк (т,0) + | Ук (t Щ т е [аьЬ ], к = 1, т;
< 0 _
г к (т,0) = /к (т) + igk (т), те[ак,Ьк ],к = 1, т,_
где Ук ^) - заданные скорости движения контуров, Тк (t), к = 1, т ; ф0(г) - начальное распределение потенциала скоростей вдоль свободной границы Г(0); х = ф0(т), у = g0(т), т е (а0, Ь0) - параметрические уравнения кривой Г(0); х = /к(т), у = gk(т), т е [ак, Ьк], к = 1,т - параметрические уравнения твердых границ, Тк(0),к = 1,т; g - ускорение свободного падения. Жидкость предполагается идеальной и несжимаемой.
2. Численный метод решения задачи (2) - (5). Решение задачи (2) - (5) для уравнения Лапласа в области В(?) будем искать приближённо в дискретные моменты времени гп = пА;, п = 0, 1, 2, ..., на дискретном множестве точек
дВ(1п) = : ] = Щ.
Предположим, что в некоторый момент времени гп > 0 нам известны положения границ области ^(;п) и значения искомой функции ф(г, ;п) во всех точках г е дО(;п).
Нам нужно найти их положения и значения функции ф^, ;п+1) в точках дД;п+1) в момент времени ;п+1 = гп + А; (при заданном малом промежутке времени А;).
Положения твердых границ Тк(;п+1), к = 1,т в соответствии с (5) находятся путем их перемещения на векторы.
Ас(п) = |Ук (;)Л, к = 1, т.
Для определения положения свободной границы Г(;п+1) нам достаточно найти, как это следует из условия (4), скорости точек г е Г(;п) в момент времени ;п.
Проекции векторов и (г, гп) на орты хг касательных к кривой Г(;п) в точке г е Г(;п) могут быть вычислены как производные от функции ф(г, ;п) по направлениям хг:
. . дф . их(г,^) = ^, г е ПЛ). дх г
Проекции же векторов на орты нормалей к кривой в точках должны находиться путем численного решения смешанной задачи для уравнения Лап-
д
ласа в области 0(гп) по известным значениям ф(г, гп) и -(г, гп), соответст-
дп2
венно, в точках кривых Г(;п) и Тк(гп), к = 1,т .
Для решения этой задачи будем искать функцию ф(г, гп) в области 0(;п) в виде (1); тогда неизвестные функции рк(г) = рк(г, ;п) и q(z) = q(г, ;п) в силу граничных условий должны удовлетворять следующей системе линейных интегральных уравнений:
... _ соч \1/ _
%Рп(го) + Е | Р] (г)~,-£01 Жг - |q(z)dц(zo,z) = /к(г0), г0 еТк(^), к = 1,т;
1=1Т] рп ) |г - го| Щп )
т .
-^(го) + ] q(z)dю(zo,z)] Р] (г)1п| г - го | = q(zo), го е Д^),
Г(;п) 1=1Т1 (гп)
где
/к(го) = (V(;п)Ао), го еТк(^),к = 1,т, q(zo) = ф(гоЛ), го е Г(^),
(Яго- го) , ч соч©^г соч \гог =—^-^—, ц(го, г) = ---,
I г - го\ \г - го\
соч © го г = ^ = .
\г - го| \г - го|
(6)
Разобьем кривые Г(У„) и Tk(tn), k = 1, m, соответственно, на Nr и Nk, k = 1, m, ячеек и заменим систему (6) следующей приближенной схемой типа Н.М. Крылова - Н.Н. Боголюбова
N Аю i, Nr
Л Рг + Z Pj^T1- Asj -Z 4j A^j = f i=1'Nk
j=1 Asг j=1
Nr Nk -
Mi + ZЧгAojj -ZPjEj = gi, г = 1,Nr,
j=1 j=1
где N - суммарное число разбиений твердых границ;
(7)
f = (^t(t„),nzo), i = 1,N, k = 1,m, gi = ф(z*,tn), i = 1,Nr,
cos ©i,j+1 cos ©i,j Tz*, Zj - Z* . — . —— АЦгу = 1----,cos©i,j =-d-, i = 1,N, j = 1,Nr,
di, j+1 di,j dij
Аюij = arcsin —^—, dij = |zj - z*|, duj+1 = |zj+1 - z*|,
dijdi, j+1
Aij = (xj -x*)(y+ -y*)-(Xj+1 -x*)(yj - j,*), j = 1,N + Nr,
zj+1
Eij = Jln|z - z*|dSz, ASj = |zj+i - zj |,
zi + zi+1 ^ zi - zi+1
i = 1,Nr , j = 1,N,,T* =
TU Ф + Ф+1
Fi+1 - z,
Яг. = (-7>* ,1%), г = 1, Ы, ф? = " Фг = Ф(^, („), / = 1, Ыг.
После решения СЛАУ (7) вычисляем касательные и нормальные составляющие векторов скоростей в точках Хг е Г(^„), Хг = 1, ЫГ+1, координаты и потенциалы этих точек в момент времени ¿„+1.
3. Тестовая задача. При реализации данного алгоритма на ЭВМ на первом шаге по времени автоматически решается соответствующая «ударная задача», которая в случае кругового цилиндра имеет точное аналитическое решение. Результаты её численного решения при V = (0;1) и И = 2 (И - расстояние от центра круга до горизонтальной свободной поверхности жидкости) даны в табл. 1.
Таблица 1
Распределение потенциала ф вдоль единичной окружности (угол 0, град., отсчитывается от оси 0^).
0, град. Точное решение Численное решение
9 -0,60207 -0,60004
27 -0,53538 -0,53161
45 -0,39733 -0,39332
63 -0,19090 -0,18785
81 +0,067052 +0,068587
99 +0,34938 +0,34931
117 +0,62474 +0,62322
135 +0,86264 +0,85991
153 +1,0370 +1,0334
171 +1,1292 +1,1250
4. Численное решение нестационарной задачи. В качестве примера проведён расчёт вертикального движения кругового цилиндра, находившегося в начальный момент на глубине двух его радиусов. Результаты расчёта представлены в табл. 2. В ней дано распределение потенциала поля в момент времени ^ = 2,5 при Я = 1, V = (0;1). Для сравнения там же приведены соответствующие результаты, полученные ранее другим методом в [4].
Таблица 2
Распределение потенциалов вдоль контура цилиндра в момент времени = 2,5 (угол 6 отсчитывается от горизонтали)
0, град. Данный метод Метод [7]
-87 1,0939 1,0956
-75 1,0630 1,0646
-63 0,99259 0,99402
-51 0,88628 0,88743
-39 0,74990 0,75068
-27 0,59129 0,59163
-15 0,42034 0,42022
-3 0,24909 0,24850
9 0,091405 0,090418
21 -0,038746 -0,040009
33 -0,13189 -0,13335
45 -0,18794 -0,18964
57 -0,21545 -0,21744
69 -0,22579 -0,22807
81 -0,22836 -0,23081
87 -0,22860 -0,23109
Как видно из табл. 2, численные решения практически совпадают (то же самое имело место и во все другие моменты времени t е [0; 5]), что свидетельствует о высокой точности предлагаемого численного метода. Это объясняется тем, что СЛАУ (7) имеет малые числа обусловленности, вследствие чего алгоритм является устойчивым.
Литература
1. Белоцерковский С.М., Дворак А.В., Тесемкин Д.А. К моделированию проникновения в жидкость // Доклады академии наук СССР. 1987. Т. 296, № 6. С. 1320-1323.
2. Лаврентьев М.А., Шабат Б.В. Проблемы гидродинамики и их математические модели. М.: Наука, 1977. 408 с.
3. Терентьев А.Г., Чечнев А.В. Численное исследование входа пластины и дисков в сжимаемую жидкость // Известия академии наук СССР. Механика жидкости и газа. 1985. № 2. С. 104-107.
4. Чечнев А.В. Численное решение задачи о вертикальном движении симметричного плоского тела к свободной поверхности жидкости методом сведения к серии задач Неймана // Математические модели и их приложения: сб. науч. тр. Чебоксары: Изд-во Чуваш. ун-та, 2013. С. 62-67.
5. Шифрин Э.Г., Шубников Г.В. Численные методы решения задачи нестационарного течения жидкости с перемещающимися границами // Вычислительная математика и математическая физика. 1982. Т. 22, № 1. С. 163-170.
6. Brebbia C.A., Domínguez J. A two dimension element code for potential problems using quadratic elements. Software for Engineering Workstation, 1988, vol. 4, pp. 134-144.
7. Hirt C. W., Nichols B.D. Volume of fluid (VOF) method for the dynamics of free boundaries. Journal of Computational Physics, 1981, vol. 39, pp. 201-225.
References
1. Belotserkovskii S.M., Dvorak A.V., Tesemkin D.A. K modelirovaniyu proniknoveniya v zhidkost' [To the simulation of penetration into the liquid]. Doklady akademii nauk SSSR [Reports of the Academy of Sciences of the USSR], 1987, vol. 296, no. 6, pp. 1320-1323.
2. Lavrent'ev M.A., Shabat B.V. Problemy gidrodinamiki i ikh matematicheskie modeli [Problems of hydrodynamics and their mathematical models]. Moscow, Nauka Publ., 1977, 408 p.
3. Terent'ev A.G., Chechnev A.V. Chislennoe issledovanie vkhoda plastiny i diskov v szhi-maemuyu zhidkost' [Numerical study entry plate and disks in a compressible fluid]. Izvestiya akademii nauk SSSR. Mekhanika zhidkosti i gaza [Proceedings of the Academy of Sciences of the USSR], 1985, no. 2, pp. 104-107.
4. Chechnev A.V. Chislennoe reshenie zadachi o vertikal'nom dvizhenii simmetrichnogo ploskogo tela k svobodnoi poverkhnosti zhidkosti metodom svedeniya k serii zadach Neimana [Numerical solution of the problem of the vertical motion of a symmetric flat body to the free surface of the liquid by the reduction method to a series of Neumann problems]. Matematicheskie modeli i ikh prilozheniya: sbornik nauchnykh trudov [Mathematical models and their applications: collection of scientific papers]. Cheboksary, Chuvash State University Publ., 2013, pp. 62-67.
5. Shifrin E.G., Shubnikov G.V. Chislennye metody resheniya zadachi nestatsionarnogo tech-eniya zhidkosti s peremeshchayushchimisya granitsami [Numerical methods for solving the problem of unsteady fluid flow with moving boundaries]. Vychislitel'naya matematika i matematicheskaya fizika [Computational Mathematics and Mathematical Physics], 1982, vol. 22, no. 1, pp. 163-170.
6. Brebbia C.A., Dominguez J. A two dimension element code for potential problems using quadratic elements. Software for Engineering Workstation, 1988, vol. 4, pp. 134-144.
7. Hirt C.W., Nichols B.D. Volume of fluid (VOF) method for the dynamics of free boundaries. Journal of Computational Physics, 1981, vol. 39, pp. 201-225.
ЧЕЧНЕВ АЛЕКСАНДР ВЛАСОВИЧ - кандидат физико-математических наук, доцент кафедры компьютерных технологий, Чувашский государственный университет, Россия, Чебоксары (dikt@chuvsu.ru).
CHECHNEV ALEXANDR - candidate of physical and mathematical sciences, assistant professor of Computer Technology Chair, Chuvash State University, Russia, Cheboksary.
