CC BY
8
МЕХАНИКА
Челябинский физико-математический журнал. 2021. Т. 6, вып. 4- С. 440-448.
УДК 621.455+629.76.085.5 Б01: 10.47475/2500-0101-2021-16404
ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ
КАВИТАЦИОННОГО ОБТЕКАНИЯ ТЕЛА ВРАЩЕНИЯ
В. И. Пегов
Южно-Уральский научный центр УрО РАН, Миасс, Россия
Государственный ракетный центр имени академика В. П. Макеева, Миасс, Россия src@makeyev.ru
Статья посвящена численному моделированию осесимметричного кавитационного обтекания тела вращения. Метод решения задачи основан на замене кавитационного обтекания тела обтеканием гидродинамических источников, непрерывно распределённых по поверхностям тела и каверны. Решение заключается в нахождении формы каверны, распределения давления и присоединённых масс по поверхности тела.
Ключевые слова: гидродинамика, каверна, кавитатор, замыкатель, число кавитации, метод потенциалов, свободная граница.
Разработка численного решения кавитационного обтекания тела вращения является актуальной в связи с выбором удобообтекаемых осесимметричных форм подводных аппаратов [1; 2].
Расчёт кавитационного течения будем проводить при следующих предположениях: жидкость идеальная невесомая и несжимаемая, вызванное течение жидкости — потенциальное [2-6]. При этих предположениях потенциал скорости Ф должен удовлетворять уравнению Лапласа
ДФ = 0 (1)
и следующим граничным условиям:
- непротекания на смоченной поверхности тела шт
дФ
к к = V ■ п; (2)
- постоянства давления на границе каверны шк
Р к = Рк; (3)
- потенциал Ф при удалении в бесконечность стремится к нулю:
Ф ^ 0 при л/х2 + у2 ^ то. (4)
В случае безвихревого движения интеграл Коши — Лагранжа служит для выражения давления в каверне Рк через кинематические элементы Ф, ит, V:
дФ + и! + Рк = Р^ + V!
дЬ + 2+ р р + 2 ,
где Ф = Ух—Фабс. — потенциал течения при обращении движения; ит — касательная скорость на границе каверны; У — скорость потока; — давление в жидкости. Введём безразмерные величины
- Уг _ х I = —, х = —
Лн Лн
безразмерное число кавитации
ит = Ф
ит = —, Ф =-
т У , У Л н
а = —
Лн дУ
У2 ~дь'
а =
РР
1 то 1 к
Здесь Лн есть радиус насадки. Тогда с помощью интеграла Коши получим выражение производной от потенциала по времени
Лагранжа
д Ф
Ж
(1 + а + 2аФ — и?) .
При установившемся кавитационном обтекании тела производная от потенциала по времени равна нулю = ^, и из последнего уравнения можно найти точное значение касательной скорости ит на границе каверны
ит =
7ТГ
а.
(5)
Далее будем опускать черту над буквенными обозначениями, считая, что все величины безразмерные.
Задача о кавитационном течении сводится к смешанной краевой задаче, так как на поверхности тела задана производная от потенциала (2), а на поверхности каверны может быть задано значение потенциала. Действительно, поскольку известно
. л- дФ
значение касательной скорости на свободной границе ит = у 1 + а и —— = ит, то с
дв
помощью интегрирования ит по контуру в найдём значение потенциала на границе каверны
Ф(в) = Ф(в^) + / ит (а)(а; вн < в < вр.
^ вм
(6)
Здесь в^ и вР — длины контура в в точках схода струй с кавитатора и с замыкателя.
При расчёте кавитационных течений в невесомой жидкости применим схему Ря-бушинского, обладающую продольной симметрией относительно плоскости, проходящей через максимальный диаметр каверны (рис. 1).
Рис. 1. Симметричная схема Рябушинского
в
При этом число кавитации, форму кавитатора будем считать заданными, а форму границ каверны будем определять из решения задачи.
Распределим по поверхности тела и свободной поверхности каверны источники с поверхностной интенсивностью ^(Р), которая является функцией точки Р поверхности тела и свободной поверхности шт + шк (рис. 1). Выражение для абсолютного потенциала движения жидкости в цилиндрической системе координат х, г, в в этом случае имеет вид
Фабс (P) = -
I R
' Шт+Шк R
где R = (x — £)г + (r cos в — р cos id)j + (r sin в — p sin <d)k — расстояние от расчётной точки P(x, r, в) до текущей точки Q(£,p,$).
С помощью преобразований граничное условие непротекания (2) приводится к интегральному уравнению Фредгольма II рода для интенсивности источников g(s) на поверхности тела [2]:
Г sl
g(s) = rr' — g(s)Kw(s,a)da; 0 ^ s ^ sn; sp ^ si ^ si. (7)
J 0
Для нахождения интенсивности источников на свободной поверхности используем полученное из условия постоянства давления (3) выражение для потенциала на свободной поверхности (6). С учётом этого получим уравнение [1]
Г sl
Ф^) = x — g(a)Ki2(s,a)da; sN ^ s ^ sP. (8)
0
Уравнение (8) будем рассматривать как интегральное уравнение для определения g(s) на границе каверны, которое является интегральным уравнением Фредгольма I рода.
В результате математическое описание задачи сводится к двум интегральным уравнениям Фредгольма I и II рода (8), (7), которые эквивалентны уравнению Лапласа (1) c граничными условиями (2)-(4). Для решения задачи эту систему уравнений необходимо дополнить уравнениями для определения формы границ каверны. Сначала задаётся приближённо форма каверны и затем с помощью метода установления проводится её уточнение. Схема расчёта реализуется в виде итерационного процесса, где на каждом шаге итераций проводится уточнение координат свободной поверхности.
Согласно кинематическому условию свободную поверхность можно рассматривать как непроницаемую поверхность [4], на которой выполняются условия непротекания, и интегральное уравнение (7) можно распространить при этом на всю поверхность (смоченную поверхность тела и границу каверны):
g(s) = rr' — Г g(s)Kio(s, a)da; 0 ^ s ^ si. (9)
0
По полученной при решении этого уравнения интенсивности д можно рассчитать потенциал Ф и скорость UT. При приближённом же задании формы каверны д не удовлетворяет интегральному уравнению (8) и производная от потенциала по времени
дФ 1 / 2л
= 2 V1 + ° — ит) ; sn ^ s ^ sp,
как и в нестационарном случае, будет отличаться от нуля. Следуя методу уста-
дф дъ
дФ г,
новления, с помощью интегрирования по времени с шагом п уточним значения
потенциала на границе каверны
Ф*(в) = Ф(в) + дФП; вм ^ в ^ вр. (10)
Затем значение интенсивности $(в), соответствующее динамическому граничному условию (3), найдём из совместного решения интегральных уравнений (7), (8):
гвь
д(в) = гг' I д(в)Кю(в, а)(а; 0 ^ в ^ вм; вр ^ в ^ вг;
'о
св1
Ф*(в) = х(в) — / ^(а)К12(в, а)(а; в^ ^ в ^ вР. о
(11)
По полученному из решения системы интегральных уравнений (11) значению интенсивности $*(в) рассчитываются составляющие скорости и£, иГ и ит на свободной границе, которые удовлетворяют динамическому граничному условию. Потребуем, чтобы при полученных значениях иХ, иГ и ит выполнялось на свободной поверхности также и кинематическое граничное условие, т. е.
дх их дг иг (12)
дв ит' дв ит
Составляющие скорости uтj, uxj и urj = N N + 1,... , Р) в расчётных точках каверны вычисляются по формулам
Г в1 г' Г в1
и^ = xj + / ,а)(а; Uxj = 1 — # — + / ^(а)^^,а)(а;
Xj
./0 ' ./0
Х' Г'
Urj = --+ ,а)(а.
г0
Выражения для ядер К10(в,а), К11(в,а), К12(в,а) приведены в [2], ядра К13(в,а) и К14(в, а) имеют вид
К1з(в, а) = 2(Х — $) ^А; К14(в, а) = Г2?^ + —) А.
г12 г12 г
Интегрируя систему дифференциальных уравнений (12), обеспечиваем выполнение динамического граничного условия на этом временном шаге. Многократное повторение вычислений по формулам (9)-(12) позволяет определить форму каверны.
Интегральные уравнения Фредгольма I и II рода (9) и (11) решаются с помощью замены интегральных уравнений конечной системой алгебраических уравнений методом Фредгольма [2]. В результате получим систему М уравнений с М неизвестными. Решив эту систему алгебраических уравнений каким-либо численным методом, например методом Гаусса, найдём значения интенсивностей в расчётных точках <7ъ д, ..., дм.
Разработанный численный метод расчёта осесимметричных кавитационных течений был реализован в виде программы для ЭВМ. Основным блоком в составленной программе является блок формирования матриц, соответствующих системе линейных алгебраических уравнений. Составляющие скорости их, иг, ит и потенциал Ф вычисляются в расчётных точках с использованием блока формирования матриц.
По схеме Рябушинского проведены расчёты кавитационного обтекания диска (в = 180°) и острых конусов с углом раствора (в = 90° и 120°) при изменении числа кавитации от 0.02 до 0.10. При вычислении каждого отдельного профиля каверн требуется от трёх до шести итераций. Рассчитанные профили каверн приведены в виде графиков на рис. 2. Из приведённых графиков следует, что длина 1к и максимальный диаметр RK каверны зависят как от числа кавитации, так и от угла раствора кавитатора: при уменьшении числа кавитации и увеличении угла раствора 1к и RK возрастают. Эта закономерность подтверждается экспериментальными данными [4].
С7=0,02
О 10 20 30 40 50 60 70 80 90 (X-Xh)/RH
а=0,06
о=0,1
-/3=180°;--/3=120°; . . . 0=90°
Рис. 2. Расчётные профили каверн за диском и острыми конусами в невесомой жидкости а = 0.02; 0,06; 0,10
Проведены также более подробные исследования по распределению давления и потенциала течения по поверхности кавитатора, которые недостаточно представлены в имеющихся данных [3; 5; 6].
Значения коэффициента давления Р = 2(Р — Р^)/рУ2 на смоченной поверхности кавитатора в зависимости от числа кавитации а и координаты т/Ян приведены в табл. 1 для диска (в = 180°) и острых конусов (в = 90° и в = 120°).
Таблица 1
Диск Острые конусы
в 180° 120° 90°
а 0.02 0.06 0.10 0.02 0.06 0.10 0.02 0.06 0.10
Т/Кн Р
0 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
0.2 0.99 0.99 0.99 0.91 0.91 0.91 0.78 0.78 0.78
0.4 0.97 0.97 0.97 0.83 0.83 0.82 0.67 0.67 0.66
0.5 0.95 0.95 0.95 0.78 0.77 0.77 0.63 0.63 0.63
0.6 0.92 0.92 0.92 0.72 0.71 0.71 0.56 0.56 0.56
0.7 0.88 0.88 0.87 0.66 0.65 0.65 0.48 0.48 0.48
0.8 0.81 0.81 0.80 0.58 0.56 0.56 0.39 0.39 0.39
0.9 0.68 0.68 0.66 0.49 0.44 0.42 0.27 0.27 0.27
0.95 0.55 0.54 0.52 0.42 0.34 0.32 0.18 0.18 0.16
0.98 0.43 0.42 0.39 0.29 0.28 0.26 0.15 0.11 0.04
0.99 0.39 0.38 0.34 0.05 0.01 —0.04 0.06 0.02 —0.02
1.00 —0.02 —0.06 —0.10 —0.02 —0.06 —0.10 —0.02 —0.06 —0.10
Значения абсолютного потенциала течения Фабс., равного Фабс. = Ф — х на смоченной поверхности кавитатора в зависимости от числа кавитации а и координаты т/Ян приведены в табл. 2 для диска (в = 180°) и острых конусов (в = 90° и в = 120°).
Таблица 2
Диск Острые конусы
в 180° 120° 90°
а 0.02 0.06 0.10 0.02 0.06 0.10 0.02 0.06 0.10
Т/Кн Фабс.
0 1.44 1.42 1.41 1.01 0.98 0.91 0.74 0.70 0.66
0.2 1.43 1.42 1.40 1.08 1.04 0.98 0.82 0.78 0.74
0.4 1.41 1.39 1.38 1.11 1.06 1.01 0.88 0.84 0.80
0.5 1.39 1.37 1.36 1.12 1.07 1.02 0.89 0.85 0.81
0.6 1.36 1.35 1.33 1.12 1.07 1.02 0.90 0.86 0.82
0.7 1.33 1.32 1.30 1.12 1.07 1.02 0.90 0.86 0.82
0.8 1.29 1.28 1.26 1.10 1.06 1.00 0.90 0.85 0.82
0.9 1.25 1.23 1.21 1.08 1.03 0.98 0.89 0.85 0.81
0.95 1.21 1.20 1.18 1.07 1.02 0.97 0.89 0.84 0.80
0.98 1.20 1.18 1.16 1.06 1.01 0.96 0.89 0.84 0.80
0.99 1.19 1.17 1.15 1.05 1.00 0.95 0.89 0.84 0.79
1.00 1.18 1.16 1.14 1.04 1.99 0.94 0.88 0.83 0.79
В табл. 3 приведены значения длины 1к и максимального радиуса Як каверны, рассчитанные по разработанному методу для диска (в = 180°) и острых конусов (в = 90° и в = 120°). Для диска получаются наибольшие значения 1к и Як. Здесь же представлены результаты расчётов 1к и Як для диска по формулам
1к = Ян (1-92 — 3а); Як = Ян,/084Щ (13)
а V с
полученным с использованием опытных данных [4]. Результаты расчётов хорошо согласуются с аппроксимационными зависимостями (13). Таблица 3
в 180° 120° 90°
а 0.02 0.06 0.10 0.02 0.06 0.10 0.02 0.06 0.10
lK/Rn расчёт 98.4 28.3 15.8 86 24.4 13.8 75.6 21.5 12.1
по фор. (13) 93 29 16.2 - - - - - -
Rk/Rh расчёт 6.80 3.96 3.20 5.85 3.34 2.78 5.07 2.96 2.44
по фор. (13) 6.55 3.85 3.04 - - - - - -
Коэффициент сопротивления кавитатора Cx вычисляется по формуле
2 fSN
Cx = о- + ^ гг'(! - )ds. Jo
При внезапном ускорении тела, обтекаемого в режиме кавитации, в предположении, что кавитационная зона не подвергается внезапному изменению, коэффициент продольной присоединённой массы кавитатора Ax вычисляется по формуле
. м 2п fSN
А' = PR = -r>L ГГ(Ф - x)ds-
Для рассмотренных форм кавитатора результаты расчёта коэффициентов сопротивления Cx и Ax приведены в табл. 4. Здесь же приводится сравнение рассчитанных значений Cx с их значениями, определёнными по формуле [4]
Cx = Cxo (1 + о) (14)
при значениях Cx0 = 0.82; 0.607; 0.465 соответственно для диска и конусов (в = 180°; 120°; 90°). Видим, что уменьшение угла раствора конического кавитатора позволяет значительно уменьшить его сопротивление.
Таблица 4
в 180° 120° 90°
а 0.02 0.06 0.10 0.02 0.06 0.10 0.02 0.06 0.10
СХ расчёт 0.839 0.876 0.909 0.630 0.652 0.677 0.470 0.492 0.509
Сх по фор. (14) 0.836 0.869 0.902 0.619 0.643 0.668 0.474 0.493 0.512
Ax = ß/pRl 4.39 4.23 4.04 2.90 2.86 2.83 2.50 2.46 2.42
Из проведённого сравнения результатов расчётов по симметричной схеме Рябу-шинского с зависимостями (13), (14), выведенными на основании обобщения экспериментальных данных [4], можно заключить, что разработанный метод подтверждается опытными данными.
Список литературы
1. Пантов Е.М., Махин Н.Н., Шереметов Б. Б. Основы теории движения подводных аппаратов. Л. : Судостроение, 1973.
2. ДегтярьВ.Г., Пегов В. И. Гидродинамика подводного старта ракет. М. : Машиностроение/Машиностроение — Полет, 2009.
3. Ш^епеленко В. Н. К расчёту кавитационных течений // Журн. приклад. механики и тех. физики. 1968. № 5. C. 100-105.
4. Логвинович Г. В. Гидродинамика течений со свободными границами. Киев : Науко-ва думка, 1969.
5. Амромин Э. Л., Иванов А. Н. Осесимметричное обтекание тел в режиме развитой кавитации // Изв. АН СССР. Механика жидкости и газа. 1975. № 3. C. 37-42.
6. Пегов В. И., Мошкин И. Ю. Расчёт гидродинамики кавитационного способа старта ракет // Челяб. физ.-мат. журн. 2018. Т. 3, вып. 4. С. 476-485.
Поступила в 'редакцию 18.10.2021. После переработки 06.11.2021.
Сведения об авторе
Пегов Валентин Иванович, доктор технических наук, профессор, ведущий научный сотрудник, Южно-Уральский научный центр УрО РАН; главный научный сотрудник, Государственный ракетный центр имени академика В. П. Макеева, Миасс, Россия; e-mail: src@makeyev.ru.
Chelyabinsk Physical and Mathematical Journal. 2021. Vol. 6, iss. 4. P. 440-448.
DOI: 10.47475/2500-0101-2021-16404
NUMERICAL SIMULATION OF CAVITATIONAL FLOW AROUND A ROTARY BODY
V.I. Pegov
South Ural Research Center of the Ural Branch of RAS, Miass, Russia Academician V.P. Makeyev State Rocket Centre, Miass, Russia src@makeyev.ru
A numerical simulation of an axisymmetric cavitational flow around a rotary body is performed. The problem solving method is based on substitution of the cavitational flow around a body for the flow around hydrodynamic sources continuously distributed over the body and the cavity surfaces. The solution consists in defining a cavity shape, a distribution of the pressure and apparent masses over the body surface.
Keywords: hydrodynamics, cavity, cavitator, closer, cavitation number, method of potentials, free boundary.
References
1. PantovЕ.N., MakhinN.N., Sheremetov B.B. Osnovy teorii dvizheniya podvodnykh apparatov [Foundations of the submersible motion theory]. Leningrad, Sudostroyeniye Publ., 1973. (In Russ.).
2. DegtiarV.G., Pegov V.I. Gidrodinamika podvodnogo starta raket [Hydrodynamics of rocket launches from underwaters]. Moscow, Mashinostroeniye / Mashinostroeniye — Polyot Publ., 2009. (In Russ.).
3. Shepelenko V.N. Calculations for cavitating flows. Journal of Applied Mechanics and Technical Physics, 1968, vol. 9, iss. 5, pp. 598-601.
4. Logvinovich G.V. Gidrodinamika techeniy so svobodnymi granitsami [Hydrodynamics of flows with free boundaries]. Kiev, Naukova dumka Publ., 1969. (In Russ.).
5. AmrominE.L., IvanovА.N. Axisymmetric flow around bodies in the developed cavitation mode. Fluid Dynamics, 1975, vol. 10, iss. 3, pp. 396-400.
6. Pegov V. I., MoshkinI.Yu. Computation of the rockets launch hydrodynamics. Chelyabinsk Physical and Mathematical Journal, 2018, vol. 3, iss. 4, pp. 476-485.
Accepted article received 18.10.2021. Corrections received 06.11.2021.
