Некоторые аспекты математического моделирования дозвукового обтекания тел с отрывом потока в донной области Текст научной статьи по специальности «Физика»

_МЕЖДУНАРОДНЫЙ НАУЧНЫЙ ЖУРНАЛ «ИННОВАЦИОННАЯ НАУКА» №12-4/2016 ISSN 2410-6070_

УДК 531.6.011.32:532.582.4:517.958

В.Н. Тимофеев

к.т.н.

МГТУ им. Н.Э. Баумана г. Москва, Российская Федерация

НЕКОТОРЫЕ АСПЕКТЫ МАТЕМАТИЧЕСКОГО МОДЕЛИРОВАНИЯ ДОЗВУКОВОГО ОБТЕКАНИЯ ТЕЛ С ОТРЫВОМ ПОТОКА В ДОННОЙ ОБЛАСТИ

Аннотация

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

Ключевые слова

Математическое моделирование, дозвуковое отрывное обтекание, эквивалентное тело.

Основанная на концепции вязко-невязкого взаимодействия [1] методика математического моделирования пространственного обтекания тел дозвуковым потоком газа была предложена автором в работе [8]. Затем данная методика была применена к изучению отрывного обтекания осесимметричных тел ([9]). В соответствии с одним из основных положений концепции вязко-невязкого взаимодействия давление на поверхности тела, обтекаемого потоком вязкой среды, находилось по результатам расчета невязкого обтекания некоторого эквивалентного тела, называемого также телом вытеснения. Так как реализация концепции вязко-невязкого взаимодействия в полном объеме требует больших вычислительных ресурсов, то в методике использовался подход, основанный на сочетании численного и физического экспериментов с априорным определением конфигурации эквивалентного тела. В данной работе предлагается другой вариант построения эквивалентного тела.

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

обтекании тел предполагается, что линия отрыва потока L0 совпадает с контуром донного среза L^ (рис. 1).

Передний участок X i поверхности эквивалентного тела, простирающийся вниз по потоку от носовой

оконечности обтекаемого тела, за линией стыковки Lc дополняется хвостовым участком X 2, моделирующим влияние спутного следа. В совокупности два указанных участка и образуют поверхность эквивалентного тела

Рисунок 1- Конфигурация обтекаемого и эквивалентного тел.

_МЕЖДУНАРОДНЫЙ НАУЧНЫЙ ЖУРНАЛ «ИННОВАЦИОННАЯ НАУКА» №12-4/2016 ISSN 2410-6070_

Для построения переднего участка Z i эквивалентного тела необходимо сместить соответствующие

точки поверхности Zq обтекаемого тела в направлении нормали к Zq. на расстояния, равные толщине S* вытеснения пограничного слоя. При этом точки, принадлежащие линии стыковки Lc , также окажутся смещенными на расстояния, равные S* , от соответствующих точек линии отрыва потока L0 . Как следует из теории пограничного слоя, давление в произвольной точке поверхности Zq обтекаемого тела должно

равняться давлению в соответствующей точке переднего участка Z i поверхности эквивалентного тела.

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

необходимость задавать конфигурацию поверхности Z 2 в первом приближении, а затем проводить

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

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

Возможные способы снижения трудоемкости расчетов состоят в следующем. Во многих прикладных задачах значения чисел Рейнольдса, рассчитанные по скорости набегающего потока, превышают величину 106. Для таких значений чисел Рейнольдса, величины S* толщин вытеснения пограничного слоя малы по сравнению с поперечными размерами хорошо обтекаемых тел. Данное обстоятельство позволяет не осуществлять поправку на малые толщины вытеснения. Поэтому для режимов течения с линией отрыва

потока близкой к контуру донного среза поверхность Zi переднего участка эквивалентного тела можно считать совпадающей с поверхностью Zq обтекаемого тела (рис. 2).

Рисунок 2 - Эквивалентное тел с полубесконечным хвостовом участком.

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

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

_МЕЖДУНАРОДНЫЙ НАУЧНЫЙ ЖУРНАЛ «ИННОВАЦИОННАЯ НАУКА» №12-4/2016 ISSN 2410-6070_

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

Согласно рекомендациям работы [2] эквивалентное тело должно представлять собой полубесконечное тело вытеснения, поперечный размер хвостового участка которого по мере удаления вниз по потоку от донного среза обтекаемого тела плавно изменяется, достигая на бесконечности некоторого конечного значения (см. рис. 2). Именно такая форма эквивалентного тела исследуется в данной работе, в отличие от предлагавшейся ранее [10] конфигурации с хвостовым участком конечной длины (рис. 3).

Рисунок 3 - Эквивалентное тел с хвостовом участком конечной длины.

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

Qu Qv Qw ^

--1---1--= О;

Qx Qy Qz

Qu Qu Qu Qu

--b u--ь v--ь w-=

Qt Qx Qy Qz

Qv Qv Qv Qv --b u--b v--b w- =

Qt Qx Qy Qz

Qw Qw Qw Qw + u--b v--b w ■

1 Qp

p Qx 1 Qp

p Qy'

1 Qp

(1)

Qt

Qx

Эу Эе р Эе

где и, V, ^ - координаты вектора V скорости частиц среды, определяемые в декартовой системе координат, связанной с обтекаемым телом; р и р - плотность и статическое давление.

Поскольку использовалась декартова система координат, связанная с обтекаемым телом, то граничное условия непротекания поверхности эквивалентного тела записывалось в следующем виде:

V (М0)-п (М0)| Е—0 , (2)

где Мо(Хо, Уо, Ео) - произвольная точка на поверхности X ; п (М0)| Е - орт вектора нормали к указанной поверхности.

В соответствии с граничным условием затухания возмущений вектор скорости частиц среды

V (М0) — и -1 + V - j + ^ - k должен стремиться на бесконечности (при К — Х02 + у02 + ^ да ) к вектору скорости набегающего потока

К

M0 ^да

или

u ■

да

i + v ■ j + w ■ k

да J да

<

_МЕЖДУНАРОДНЫЙ НАУЧНЫЙ ЖУРНАЛ «ИННОВАЦИОННАЯ НАУКА» №12-4/2016 ISSN 2410-6070_

Предполагалось, что вне поверхности эквивалентного тела течение является потенциальным. Поэтому

вектор возмущенной скорости, равный разности (V(Mq) — V*) , считался равным градиенту потенциала <P(Mq ) возмущенных скоростей, а вектор скорости потока находился по следующей формуле:

V (Mq) = V0+V^(Mq) , (3)

д - д - д г

где V =-1 Н--j Н--k - оператор Гамильтона.

дг0 ду0 dz0

После определения потенциала P(Mq) , вектора скорости V(Mq) и его модуля статическое давление может быть найдено с помощью интеграла Бернулли.

Для несжимаемой среды уравнение неразрывности записывается в виде divV = О . После замены вектора скорости потока выражением, равным правой части равенства (3), уравнение неразрывности

преобразовывалось в соотношение dlVgradp =Q, представляющее собой уравнение Лапласа относительно потенциала возмущенных скоростей:

Ар = 0 ,

. д2 д2 д2 где А =-2 Н--2 Н--2 - оператор Лапласа.

дхо дУо дг0

Подстановка указанного выражения в условие непротекания (2) позволило получить равенство следующего вида:

(V + gradp) • n (Mq))| 2= 0 .

Наличие взаимосвязи между градиентом и производной по направлению нормали

др(М о)

(gradp•n (М Q)| 2 = ■

дп

2

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

др( Mq)

дп

= — V*- п(Мо)

2

2. (4)

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

Таким образом, потенциал возмущенных скоростей должен являться решением внешней задачи Неймана для уравнения Лапласа:

Ар(Мо ) = 0, (5)

др( Мо )

= —V •П (Мп ) ^ да v о'

2, (6)

дп

<РМ0 ) — 0, М0 -да, (7)

Уф(М0) — 0, М0 — да, (8)

Внешняя задача Неймана имеет единственное решение, которое записывалось в виде потенциала

двойного слоя

*Mo) = f Jr-^g<M)dГ , (9)

г

где Х0, У0, - точка, в которой вычисляется потенциал возмущенных скоростей и скорость потока; п(М) и g(М) - орт вектора нормали к поверхности X и поверхностная плотность потенциала

двойного слоя, определяемые в точке М(х, у, г) , расположенной на элементе поверхности d X; Г - вектор, соединяющий точку Мс точкойМо; Г - модуль указанного вектора, вычисляемый по следующей формуле:

= V(x - xo)2 +(У - Уо)2 +(z - zo)2 .

Потенциал двойного слоя является решением уравнения Лапласа (5) и удовлетворяет граничным условиям (7)-(8) затухания возмущений на бесконечности. Для выполнения граничного условия непротекания (6) поверхностная плотность g(М) потенциала двойного слоя должна являться решением

следующего интегрального уравнения [4] :

( Э г(7-п(М)л ^

J_

4п

Q № g (м )d z

= -Vда■ n (Mo) .

Здесь М0 - текущая точка, принадлежащая поверхности X , а нормальную производную — следует

Эп

понимать в смысле, определенном, в работе [5].

Поверхностная плотность потенциала двойного слоя находилась с помощью численного моделирования. Для этого поверхность X эквивалентного тела аппроксимировалась конечным числом

панелей X^ ; к — 1,...,N , причем панели имели форму многоугольников. Суммарное число панелей,

аппроксимирующих поверхность X , обозначалось через N .В данной работе рассматривалось эквивалентное тело с полубесконечным хвостовым участком. Предполагалось, что на некотором удалении

вниз по потоку от донного среза поперечное сечение поверхности X 2 хвостового участка будет сохранять

постоянную форму и размеры. Поэтому у некоторых из многоугольников X^ одна из сторон располагалась

на бесконечном удалении вниз по потоку от обтекаемого тела. На каждой из панелей X^ поверхностная

плотность потенциала двойного слоя g(М) считалась постоянной и равной величине , к — 1,...,N .

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

N 1 сГ - п(М)

*м.)=Г * 4- J ^ z

*—1 { г

и формулу для вычисления скорости потока

N 1 г г -п(М)

V (мо) = Уда+Г g, (f VJ -^M>dZ). (10)

k=i 4ж Z г

Градиент потенциала двойного слоя, имеющего постоянную поверхностную плотность на панели

X ^ , равен скорости, индуцированной замкнутой вихревой нитью ^к , расположенной на границе ЭХк

указанной панели , при условии, что циркуляция Г ^ этой вихревой нити, равняется величине —gк . Скорость, индуцированная замкнутой вихревой нитью, определяется с помощью формулы Био-Савара. Поэтому

градиент потенциала двойного слоя панели 2 ^ с постоянной поверхностной плотностью gk был преобразован следующим образом:

g^i7 - п (М £l \ds Х Г 4ж ^ r3 4ж} r3 '

После введения в рассмотрение вектора функции скорости, равного вектору скорости,

индуцированной в точке Mq k-ой вихревой нитью L с единичной циркуляцией

—,-к, ч 1 f ds Х 7

Wk (Мо) = — J —— ,

Г r

формула (10) преобразовывалась к виду, позволяющему применять метод дискретных вихрей [6]:

_ _ N _

V (Mq) = V*+Z rkWk (MQ). (11)

k=1

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

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

Для определения неизвестных циркуляций Г^ , k = 1,..,N требовалось, чтобы граничные условия непротекания поверхности эквивалентного тела Z выполнялись в контрольных точках Cv, v = 1,.., N,

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

_ _ N _

= (V (С) - Уда) • п(С) = X (С) • п (С),

потенциала двойного слоя непрерывны, поэтому с учетом соотношений

д¥с С)

дп

граничное условие непротекания представлялось в следующем виде:

К(С)• п(С)) = п(С) , у = l,..,N , (12)

где п (Су ) - орт вектора нормали к панели 2к в контрольной точке Су.

В соответствии с равенствами (12) циркуляции Г^ определялись из системы линейных алгебраических уравнений

N

X аукГк = К, у = 1,.., N , (13)

к=1

в которой коэффициенты и правые части вычислялись следующим образом:

ак = ^к(С) •п(С X у = 1,.., N ; к = 1,.., ж ; К=-К^п(С) , у = 1,.., N

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

составленная из коэффициентов аук , получается вырожденной, поэтому система (13) решалась с использованием регуляризирующей переменной [7].

После нахождения неизвестных циркуляций Гд скорость потока в точках, расположенных вне

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

(C , (C )■

AV(C)=^P4(CV)) , (14)

от от.

2

Og(C) 0g(Cv)

причем частные производные - и - должны вычисляться по направлениям двух

дтх дт2

неколлинеарных ортов Т- (C ) и Т2 (Cv), лежащих в касательной плоскости, проходящей через

контрольную точку C . Так как на каждой панели 2 v поверхностная плотностью g(C ) была постоянной и равнялась величине (—г ), то в формуле (14) частные производные определялись численным дифференцированием по значениям циркуляций V-ого и соседних с ним вихревых многоугольников. После процедуры вычисления значений указанных частных производных и нахождения векторов AV(C v)

скорость потока в контрольных точках определялась как предельное значение градиента потенциала

двойного слоя

— — " - 1 -

V (C ) = F +Х rkwk (C ) + - AV (C ).

k=1 2

Статическое давление находилось из интеграла Бернулли

p = Р» + 2 А, V — V2) , (15)

где V - модуль вектора скорости в рассматриваемой точке течения; Vx - модуль вектора скорости

набегающего потока. В заключение определялись безразмерные коэффициенты давления

2^

С, = 2(p -px)l(pJ2) ,

которые в соответствии с интегралом Бернулли рассчитывались по следующей формуле:

с, = 1 - (V / V )2. (16)

Изучалось обтекание осесимметричных тел при нулевом угле. При таких условиях эквивалентное тело оказывалось осесимметричным (рис. 4).

dm -d

i X

- I

h и

Рисунок 4 - Геометрические размеры эквивалентного тела.

Исследовалось обтекание цилиндрических тел, снабженных головными частями оживальной формы. Характерным линейным размером считался диаметр миделевого сечения цилиндрической части dm . Длина

головной оживальной части обтекаемого тела обозначалась через Ig , а основной цилиндрической части -

через 1с . Суммарная длина обтекаемого тела 1Ь равнялась сумме 1ё + 1с .

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

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

от донного среза обтекаемого тела, диаметр хвостового участка <(х) уменьшался пропорционально

косинусу разности Х — ¡ъ от значения, равного диаметру миделевого сечения dm , до величины ds , которая

оставалась постоянной при дальнейшем увеличении координаты Х.

Такой подход позволил характеризовать форму полубесконечного хвостового участка эквивалентного

тела всего двумя параметрами: длиной I а и минимальным диаметром ds. Для получения замкнутой расчетной методики эти параметры определялись следующим образом. При известном значение минимального диаметра ds использование формулы Хорнера [12] и приемов, описанных в работе [11], позволяло находить параметр как функцию суммарной длины обтекаемого тела 1Ь и числа Рейнольдса

Re, определяемого по скорости набегающего поток и по длине 1Ъ .В свою очередь, диаметр ds

рассматривался как размер эквивалентного тела на бесконечности, что позволяло (см. [2]) выразить его через аэродинамический коэффициент суммарного сопротивления, составляющими которого являлись коэффициенты сопротивления давления и трения, а также коэффициент донного сопротивления.

Проведенные расчеты показали, что в исследованных диапазонах значений Iъ = 6...12 и Яе—10й...2*10/ можно приближенно принять величину минимального диаметра ds равной 0,3.

Рисунок 5 - Зависимость безразмерной скорости V на поверхности эквивалентного тела от безразмерной продольной координаты Х и безразмерной длины обтекаемого тела 1Ъ : кривая 1 - 1Ь — 6 ; 2 - Iь — 8 ; 3 -

Ть —10 ; 4 - Тъ —12 .

Рисунок 6 - Зависимость коэффициента давления Ср на поверхности эквивалентного тела от безразмерной

продольной координаты Х и безразмерной длины обтекаемого тела Iь : кривая 1 - ¡ь = 6 ; 2 - ¡ь = 8 ; 3 -

Тъ =10; 4 - Ть =12 .

В качестве иллюстрации на рис. 5 и 6 представлены полученные по результатам численного моделирования распределения безразмерной скорости V = V / V, и коэффициента давления Ср по поверхности обтекаемого тела. По оси абсцисс откладывались значения безразмерной продольной координаты X = X / йт . Число Рейнольдса равнялось 5*106 . Безразмерным значениям суммарных длин

обтекаемых тел 1Ъ = 1Ъ / dm , равных 6; 8; 10 и 12, на рис. 5 и 6 соответствуют линии 1; 2; 3 и 4.

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

Список использованной литературы:

1. Гогиш А.В., Степанов Г.Ю. Турбулентные отрывные течения / Глав. ред. физ.-мат.лит. - М.: Наука, 1979. - 368с.

2. Гогиш А.В., Степанов Г.Ю. Отрывные и кавитационные течения - М.: Наука, 1990. - 384с.

3. Лойцянский Л.Г. Механика жидкости и газа / Глав. ред. физ.-мат.лит. - М.: Наука, 1987.-676с.

4. Белоцерковский С.М., Лифанов И.К. Численные методы в сингулярных интегральных уравнениях и их применение в аэродинамике, теории упругости, электродинамике / Глав. ред. физ.-мат. лит. - М.: Наука 1985.-256с.

5. Соболев С.Л. Уравнения математической физики - М.: Наука, 1992. - 432с.

_МЕЖДУНАРОДНЫЙ НАУЧНЫЙ ЖУРНАЛ «ИННОВАЦИОННАЯ НАУКА» №12-4/2016 ISSN 2410-6070_

6. Белоцерковский С.М., Ништ Н.И. «Отрывное и безотрывное обтекание тонких крыльев идеальной жидкостью - М.Наука, 1978, -352с.

7. Лифанов И. К. Метод сингулярных интегральных уравнений и численный эксперимент - М.: ТОО "Янус", 1995. - 520с.

8. Тимофеев В.Н. Математическое моделирование дозвукового пространственного обтекания тел. В Кн: Актуальные направления развития прикладной математики в энергетике, энергоэффективности и информационно-коммуникационных технологиях: сб. трудов международной научной конференции, посвященной 180-летию МГТУ им. Н.Э. Баумана. - М.: 2010. - С. 198-202.

9. Тимофеев В.Н., Бушуев А.Ю. Математическое моделирование дозвукового обтекания тел с отрывом потока в донной области. Вестник Саратовского государственного технического университета" 2012, №1 (64), выпуск 2, -С.11-10. Тимофеев В.Н., Бушуев А.Ю. Численное моделирование дозвукового отрывного обтекания осесимметричных тел. Инженерный журнал: наука и инновации, 2013, вып. 9. URL: http://engj ournal .ru/catalog/mathmodel/aero/967.html

11. Тимофеев В.Н. Математическое моделирование отрывного дозвукового обтекания осесимметричных тел с учетом донного давления. Инженерный журнал: наука и инновации, 2014, вып. 10. URL: http://engj ournal .ru/catalog/mathmodel/aero/1246.html.

12. Hoerner S. F., "Base Drag and Thick Trailing Edges, Journal of the Aeronautical Sciences, Vol. 17, No. 10, 1950, pp.622-628.

© Тимофеев В.Н., 2016