Метод расчета обтекания тела вращения любой формы при произвольном движении в идеальной жидкости Текст научной статьи по специальности «Физика»

Том I

УЧЕНЫЕ ЗАПИСКИ ЦАГИ 1970 “'

Л® 2

УДК 533.6.011.32:532.582.33

МЕТОД РАСЧЕТА ОБТЕКАНИЯ ТЕЛА ВРАЩЕНИЯ ЛЮБОЙ ФОРМЫ ПРИ ПРОИЗВОЛЬНОМ ДВИЖЕНИИ В ИДЕАЛЬНОЙ ЖИДКОСТИ

Л. А. Маслов •

Предлагается метод расчета распределения скорости, давление и потенциала на поверхности, а также в любой точке пространства вокруг тела вращения, совершающего произвольное движение в идеальной жидкости. По сравнению с известными методами в данном случае на форму тела вращения не налагается никаких ограничений и достаточная точность расчетов достигается при значительно меньших затратах времени ЭВМ. Примеры расчетов сравниваются с известными точными решениями и с экспериментальными значениями давлений на поверхности различных тел.

В настоящее время потенциальное обтекание тел вращения любой формы при произвольном движении может быть рассчитано лишь с помощью метода [1]. Однако метод [1] является весьма трудоемким, и более эффективным оказывается метод расчета течения около тела вращения при произвольном движении, изложенный в работах [2] и [3]. Метод [2], кроме того, легко распространяется на случай расчета трехмерного тела при произвольном движении [4]. Однако при использовании его на поверхность накладываются некоторые ограничения.

В настоящей статье предлагается метод расчета, представляющий собой обобщение метода [2]. При этом способ задания поверхности тела позволяет рассматривать тела вращения любой формы при произвольном движении. '

1. Основные соотношения. Если V1, у2> Уз — проекции вектора у поступательной скорости полюса А твердого тела, а £2з, £2з — проекции вектора угловой скорости £2, на связанные с телом оси координат хуг, то для потенциала возмущенных скоростей жидкости можно записать ,

в

у, г, t)^'£лVi{f)Фl(x>y, z), •

/=1 -

где введены обозначения v^ = lQ1\ ъ5 — 1&2-, г/е =/23 (/ —длина тела). ....

1

Для определения каждого из шести единичных потенциалов ФДл:, у, г) требуется решить внешнюю задачу Неймана

ДФ(. = 0;

дФг

дп

V;-tl

V,

(* = 1,2,...,6)

(1.1)

с нулевыми условиями на бесконечности. Здесь ъ1 — скорость точек поверхности тела в соответствующем простом движении; О)—поверхность, ограничивающая твердое тело; п — нормаль к поверхности м, направленная внутрь жидкости. Если решение искать в виде потенциала простого слоя

ГГ и<0

ФДР) = -.].Ь((3)'7Г’ <!-2)

ш

то граничное условие (1.1) приводится к интегральному уравнению Фредгольма второго рода относительно интенсивности слоя ^

2 *МР)+ (1-3)

О)

где Р — произвольная расчетная точка; <3 — текущая точка поверхности «>;/?= О.Р.

Уравнение (1.3) имеет единственное решение, если поверхность ш принадлежит к классу поверхностей Ляпунова [5].

Продольная ось х декартовой системы координат хуг с ортами 1}к совпадает с осью симметрии тела вращения. Начало координат помещается в носке тела (фиг. 1). Наряду с декартовой системой

Фиг. 1

рассматривается система цилиндрических координат хгЬ {у = rcosG, z = rsin0). Расчетной точке Р присваиваются координаты хгв,

текущей точке Q — координаты Spft; R = (х — S) i + (г cos 6 — р sin Ь) /+ -(— (г sin 0 — р sin 9) k.

Если полюс А выбрать на оси х на расстоянии Ха от носка тела, то для скорости точек поверхности тела при его произвольном движении можно записать

V = (l>5 Г Sin 6 —Vq Г COS 0 —v^i -f- [^e (•* — Ха) — r sin 0 — v2]j -f

+ [z/4rcos0 — vb(x — xA) — v3]k, (1.4)

где линейные размеры отнесены к длине тела I.

Уравнение образующей тела вращения

г = г(х) (1.5)

представляется в параметрическом виде:

г = r(s); x = x(s), (1.6)

где s — длина дуги образующей, отсчитываемая от начала координат. Значение параметра s, характеризующее текущую точку Q поверхности тела, обозначается буквой а.

Для представления дифференциальных элементов поверхности вводятся обозначения

= — (17)

Г ds ’ * — ds ' К ■ ’

Элемент площади поверхности равен

d(o = rdsdb.

В каждой точке поверхности тела вводится связанная с этой точкой прямоугольная система координат, единичные векторы которой равны

п = — г'i -f-x' cos 0/ f xr sin bk;

T = x' i + t' cos ey — f sin 0 k\ b = — sin 0y -j- cos 6 k;

(1.8)

здесь rt — нормаль; вектор т направлен по касательной к образующей в сторону увеличения длины дуги s; вектор b лежит в поперечной плоскости и образует с п и т правую систему координат.

В отличие от метода [2], где поверхность задается в форме (1.5) и используется производная dr/dx, в данном случае форма (1.6) позволяет представить любую кривую, ограничивающую односвязную область.

В результате подстановок основное интегральное уравнение (1.3) принимает вид

2*М«, е) = »/«(«. е)~

L 2тс

\х• r-r'{x — \)-x' Р cos (0 — &)] pda db

j l(x —S)3 + r2 + P2 —2rpcos(0 —ft)p/2 >

Я

о 0

где Ь— полная длина дуги образующей от носовой точки (х = 0, г = 0) до хвостовой (х = 1, г = 0).

Аналогичные выражения получаются для скоростей и потенциалов.

2. Расчетные формулы. Как следует из соображений симметрии, для характеристики произвольного движения тела вращения в идеальной жидкости достаточно знать параметры обтекания лишь для трех простых движений: .

— поступательного вдоль оси х со скоростью ^1=1 (/=1);

— поступательного вдоль оси у со скоростью о2=1 (/==2);

— вращательного вокруг поперечной оси, например параллельной оси г и проходящей через полюс А, с угловой скоростью 1>б= 1 а = 6).

1* 3

Суммарные значения относительной скорости и потенциала на поверхности тела вращения можно представить в виде —1 > —►

и = X [и\ х + М2 т (v2 COS 6 + V3 sin 0) + н6 T (ve cos 0 — vb sin 0)] +

+ b [a2 b (v2 sin 0 — -DgCos 0) + ueb(ve sin Ъ -\-vb cos0) — ru4]; (2.1)

Ф == / [cpj vx + Ъ (v2 cos 0 + u3 sin 0) + cp6 (ve cos 0—^5 sin 0)], (2.2)

где IIIi и uib — безразмерные компоненты относительных скоростей по осям т и Ь;

<р,— безразмерные потенциалы в соответствующем простом движении, которые и подлежат дальнейшему определению.

Решение интегрального уравнения (1.9) для указанных простых движений следует искать в виде

М«> 0) = М«); 0)-=Ms)cos0; Ms, 6) = (j.6 (s) cos 0. (2-3>

Вместо интенсивности источников (s) удобно рассматривать другие неизвестные функции gi(s), связанные с }1г соотношениями

gi = 2*r£L, (2.4)

V I

и все вычисления вести в безразмерном виде, принимая за характерный линейный размер длину тела /.

После подстановки величин (2.3) и (2.4) в уравнение (1.9) и длительных преобразований, аналогичных приведенным в работе [2], интегральные уравнения для каждого простого движения (t=l, 2, 6) приобретают вид

L

Si (s) = fi0 (s) - f gt (о) к, о (s, о) da. (2.5)

°

Здесь для известных функций fl 0 (s), равных произведению нормального компонента переносной скорости на радиус тела в данной точке, используются соотношения

fio = rr'\ /20= — rx'\ f60 = r[rr' + х'(х — хА)]. (2.6)

Ядра интегралов обозначены:

KiQ{s,a) = (BG1+x'G2)Ai, K20{s, a) = Kw(s, а) = (ВН1 + х'Н2)А;

4_.___________1 . р 2г [х' (г — р) г'(х — £)].'

кУ(х-Ц)2 + (г + р)2 ’ (*-&)* + (г-р)а ’

G^Eik*); G2 = K(k*)-E(k2);

(2.7)

= -тт [(1 + Е (*2) ~ 2к" к

кг

Н2 = -1 [(1 + 3/г'2) К (№) - (3 + А'2) Е (й2)],

&

где К (к2) и Е (А2) — полные эллиптические интегралы первого и второго рода с модулем

ьг —___________________4гр___________________ „ — . а

* (х —6)« + (г + Р)* И А . 1 к •

Аналогичным образом получаются расчетные соотношения для компонентов безразмерных скоростей и потенциалов, указанных в формулах (2.1) и (2.2):

£

И/х («) = /, 1 («) + | гI (°) К11 (5, °) С?а;

О

I

*1ь(*) =/| 2 (в) + -у | & (°) к, 2 (5, о) Л;

о

^ /с, 2 (в, о) *;

о

/и=/21 = г'; Л1 = -г'{х- хАу,

/22 — 1; f6г = x — ха\

/Си.(в, а)=^СС1 + -^С2^Л; /С,,^)-^;

(«, з) = *«, (5, а) = (снг +у /У2) Л; /С22 (5, о) = /С62 (5, а) = Н8А,

(2.8)

(2.9)

(2.10)

где кроме величин, определяемых формулами (2.7), вводятся обозначения

^ 2х (л: Е) + г (г р) ^ К(№)ш

(д;_$)2 + (г_р)2 -

Яз = ^[(1+*'2)К(*2)-2ОД].

(2.11)

В концевых точках, гдел: = 0, г = 0, либох=1, г — 0, нетрудно показать, что интегральные уравнения (2.5) всегда имеют нулевое решение и не требуют специального рассмотрения, как это делается в методе [2]. При определении значений скоростей и потенциалов в концевых точках достаточно вычислить компоненты безразмерной скорости при поперечном и вращательном движениях (I = 2,6) вдоль оси у

I

^ Ы°)Р^а

о

и1у — сг

2 [(х - I)2 + р2]3'2

(2.12)

и безразмерный потенциал при продольном обтекании

I

g1 (о) йз

?1

[(.к-^ + рТ2 ’

(2.13)

где £*=1уДЛЯ носовой точки следует положить Х = 0 и С6 = ХА, •Ф5ДЛЯ ХВОСТОВОЙ ТОЧКИ .£=1 и С6 — Ха — 1. Остальные компоненты скоростей и потенциалов в этих точках равны нулю.

Скорости в точках пространства, не принадлежащих поверхности тела, удобно вычислять в цилиндрической системе коорди-

нат лтб. Для компонентов безразмерной скорости в каждом простом движении нетрудно получить ’

При этом формула для и,есовпадает с формулой (2.8) для и1Ьг а для ядер интегралов первых двух формул (2.14) введены обозначения

Здесь г\ = {х — Е)2 + (г — р)2, а для остальных величин использованы обозначения (2.7) и (2.10). Безразмерные компоненты переносной скорости /п(х) и /14(х) равны /„=/„= 1; /„ = /^ = /„ = 0;: /б4 ~ Ха х.

Для расчетов коэффициента давлений в случае произвольного движения следует пользоваться интегралом Лагранжа. При подсчете присоединенных масс необходимо иметь в виду, что для тела вращения независимыми и не равными тождественно нулю являются четыре присоединенные массы: Хп, Х22, Х26 и Х6в. По аналогии с работой [3] нетрудно получить

где ро — массовая плотность жидкости.

3. Методы вычислений. Решение основного интегрального уравнения задачи (2.5) определяется в конечном числе дискретных расчетных точек поверхности тела. Чем сложнее форма тела, тем больше расчетных точек необходимо выбрать для достаточно точного описания образующей. Так, в составленной для ЭВМ «Минск-2» программе, по которой выполнялись примеры расчетов, может использоваться до 160* расчетных точек, произвольно расположенных на образующей тела вращения, Например, для плавных тел достаточная точность достигается при, 50 расчетных точках.

и1х(х, г)=/,з(х) + з(х, г, о)йа\

о

I

и1г(х, г) =/(4 М + / ё1 (о) К14 (х, г, з) йо;

(2.14)

о

I

о

I

и

I

*22 = лРо ( гх'ч2й8\

О

ь

(2.15)'

Ч6 = — *Ро / г \гг' + х' (х — ха)] ?2 <1$;

о

I

О

Предлагаемый метод допускает наличие конечного числа точек излома обводов образующей. В самой точке излома для формального выполнения условий Ляпунова следует предположить малый радиус закругления. В процессе расчетов закругление осуществляется автоматически применением квадратичной интерполяции между расчетными точками. В местах резкого изменения обводов необходимо уплотнение расчетных точек.

Основной трудностью является вычисление несобственных интегралов в уравнениях и формулах для скоростей и потенциалов. По аналогии с работой [2] несобственные интегралы вычисляются с помощью замены переменных:

. 0,5 Lh2 1Ч

о— s = sign (h)----—■ . (3.1)

1 — h2

Обратная зависимость определится соотношением

h (S, о) = Sign (а - S) У |a —1° + ; (3-2)

-i- = ±y>-5|(|a-SH-°,5 Lf- (3.3)

h h L

Начало отсчета переменной А совпадает с особенностью интеграла, находящейся в расчетной точке, координата которой х выступает в этой формуле как параметр. Функции (3.2) и (3.3) непрерывны и отличны от нуля всюду, за исключением самой особой точки, где (3.3) стремится к нулю как |« — з!1/2.

По аналогии с работой [4] замена используется лишь на участке

|* —о | <0,05 А, (3.4)

т. е. вблизи особой точки. Интегралы представляются в виде

I 5^0,05/. I 1/УП

| *(«)**+ ] г(»)№+ | ё{0)К^.{Ъ.Ь)

0 0 9 + 0,05£.

В формуле (3.5) при х 0,05/- первое слагаемое, а при х>0,951 второе слагаемое не учитываются, так как в этих случаях они включены в третье. При этом пределы третьего интеграла соответственно изменяются и вычисляются непосредственно по формуле (3.2), где следует положить 0 = 0 в первом случае и а = Ь, во втором.

Первые два интеграла формулы (3.5) не имеют особенности и вычисляются с помощью правила трапеций по узлам интегрирования, расположенным в расчетных точках. Последний интеграл вычисляется по правилу трапеций с постоянным шагом в новых переменных й. В этом случае замена (3.2) обеспечивает симметричное расположение узлов интегрирования относительно особенности при увеличении плотности расположения узлов в ее окрестности, что соответствует условиям существования главного значения по Коши и позволяет вычислять несобственные интегралы по правилу трапеций с постоянным шагом в новых переменных. Необходимые при этом величины £ (о), е(а) и р (а) в узлах интегрирования, оказывающихся между расчетными точками, определяются, квадратичной интерполяцией по Ньютону [6].

Интегральные уравнения решаются методом последовательных приближений по следующим схемам:

для I = 1

1 Г 1

ёю = Ло I ё1, *—1 *10 ёиг == (ё\ А “Н ^1, *—1)'

для / = 2,6

ё/0— /(01 ёиг~Ао ^ ёь *-1 0 ^а>

где к, — номер приближения.

Решение считается найденным, если

I ё 1к («) - ёи к-1 (я) |тах < 0,005 I

/г 01тах>

т. е. граничные условия (1.3) выполняются с точностью 0,5%. Число приближений колеблется от 4 для плавных тел до 10—15 — для тел сложной формы.

Информация о теле задается таблицами значений х и г в расчетных точках. Производные (1.7) вычисляются при помощи формул Ньютона по значениям л: и г в трех расчетных точках [6]. Интегралы в формулах присоединенных масс (2.15) вычисляются по правилу трапеций. Шаг переменной к был выбран равным А /г = 1 /■»/" 1536.

4. Примеры расчетов. С целью проверки изложенного метода на ЭВМ «Минск-2» для эллипсоидов вращения выполнены расчеты присоединенных масс и относительных скоростей, для которых известны точные аналитические выражения [7]. Применимость метода к определению давлений в реальной жидкости показана на примере расчета тела, имеющего вертикальный участок образующей, для которого не могло быть получено решение по методу [2], а также тела в виде комбинации конуса с цилиндром, имеющего местное резкое изменение обводов.

Фиг. 2

На фиг. 2 и 3 точками нанесены результаты вычисления предлагаемым методом безразмерных скоростей на поверхности эллиптическо-

- а Г» 1

го диска, имеющего соотношение полуосей -^-=0,1, и эллипсоида, у

которого = 9. Сплошными линиями построены точные значения скоростей. Вычисление поля скоростей проверено на примере обтекания шара единичного диаметра поступательным потоком вдоль оси х. В табл. 1 приведены значения скорости «і* в точках вертикального диаметра Ос = 0,5) и в точках горизонтального диаметра (у = 0), расположенных на различных расстояниях от поверхности шара.

Отстояние от поверхности Значения скорости и1х Отстояние от поверхности Значения скорости щх

по предлагаемому методу точные значе- ния по предлагаемому методу точные значения

По вертикальному диаметру 0,01 1,470 1,471 По горизонтальному диаметру 0,01 0,0587 0,577

0,1 1,288 1,289 0,1 0,4213 0,4213

2,0 1,004 1,004 2,0 0,9920 0,9920

Результаты подсчета коэффициентов присоединенных масс для трех эллипсоидов предлагаемым методом даны в сравнении с точными значениями в табл. 2.

Таблица 2

а Кп Кп Кее *22 *ее

Ь по предлагаемому методу точные значения

0,1 6,130 0,0751 6,184 0,0748

1.0 0,499 0,500 — 0,500 0,500 —

9,0 0,0244 0,950 0,860 0,0244 0,954 0,864

Сравнение точных величин скоростей и присоединенных масс с результатами расчетов по предлагаемому методу свидетельствуют о достаточной точности последних. Расхождение точных и расчетных величин составляет менее 1 %.

На фиг. 4 и 5 показано сравнение расчетных и экспериментальных величин коэффициента давления р в носовой части поверхности тел, содержащих вертикальный участок контура. На фиг. 4 показаны обводы тел вращения 1—3 и дано распределение р при нулевом угле атаки *. На фиг. 5 построено распределение р по одному меридиану тела вращения 3 при угле атаки а = —10°, когда этот меридиан является наветренным, и угле атаки а ==10°, когда этот меридиан становится подветренным. Совпадение результатов расчета с экспериментом оказывается достаточно хорошим. Такая же хорошая сходимость расчетных и экспериментальных значений давления

* Тела вращения 1—3 соответствуют гондолам двигателей № 25, 85 и 87, исследованным экспериментально в работе [8].

Р

ff.lt

О

О,*

(\ 1 |

> —предлагаемый ме/пов о тело / ъ » г И » 3

Ц &

г

4* УФ*

V'

X N

'\ ? 4 г—Ч 1

ч— 3

Фиг. 4

&

получена для носовой части тела, представляющей собой комбинации конуса с цилиндром (фиг. 6).

Следует заметить, что результаты, приведенные на фиг. 4, для тела с вертикальными участками образующей при нулевом угле атаки могут быть получены при помощи методов [1] и [9], в то время как случай обтекания под углом атаки может

быть рассчитан лишь по методу [1]. В методе [9] в силу использования функции тока может рассматриваться только продольное движение тела вращения. Примеры расчетов тел вращения, имеющих местные резкие изменения формы обводов, приведенные на фиг. 2 и 6, встречаются лишь у авторов работы [1]. Однако в методе [1] объем вычислений оказывается значительно большим, чем в данном случае.

* *

*

ЛИТЕРАТУРА

1. Hess J. L., S m i t h A. M. O. Calculation of potential flow about arbitrary body shapes. International Symposium on Analogue and Digital Techniques Applied to Aeronautics, Liege, 1963, Bruxelles, 1964.

2. Vandrey F. A method for calculating the pressure distribution of a body of revolution moving in a circular path through a perfect incompressible fluid. ARC R. and М., 1960, No. 3139.

3. Маслов Jl. А. Использование одного метода определения потенциала скоростей обтекания тела вращения идеальной жидкостью для расчета присоединенных масс. «Инженерный журнал», 1965, т. 5, № 4.

4. М а с л о в Л. А. Произвольное движение продолговатого тела в идеальной жидкости. Изв. АН СССР, МЖГ, 1966, № 6.

5. М и х л и н С. Г. Лекции по линейным интегральным уравнениям. Физматгиз, 1959.

6. Березин И. С., Жидков Н. П. Методы вычислений. Т. 1. Физматгиз, 1959.

7. Л а м б Г. Гидродинамика. ОГИЗ, 1948.

8. К а н т е р П. М. Распределение давления по капотам моторов воздушного охлаждения. Технические отчеты ЦАГИ, вып. 5, 1941.

9. Сидоров О. П. Решение задачи об обтекании тела вращения,

Труды КАИ, вып. XXIII, 1949.

Рукопись поступила 21 /V 1969 г.