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

УДК 531.6.011.32:532.582.4:517.958

Численное моделирование дозвукового отрывного обтекания осесимметричных тел

© В.Н. Тимофеев, А.Ю. Бушуев МГТУ им. Н.Э. Баумана, Москва, 105005, Россия

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

Ключевые слова: дозвуковое отрывное обтекание, математическое моделирование, метод дискретных вихрей.

Введение. Из множества возможных подходов к изучению отрывных течений [1-4] авторами данной работы предпочтение было отдано концепции вязко-невязкого взаимодействия [5], на базе которой сформулированы основные положения методики математического моделировании обтекания тел с отрывом дозвукового потока в донной области [6]. Они состоят в следующем. Согласно одному из основных положений теории пограничного слоя давление в произвольной точке поверхности Х0 обтекаемого тела заданной формы вниз по потоку от носовой оконечности до линии отрыва можно считать равным давлению в соответствующей точке поверхности Х1,

расположенной на нормали к поверхности X 0 и удаленной от Х0 на

*

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

Концепция вязко-невязкого взаимодействия предлагает считать поверхность Х1 передней частью поверхности так называемого эквивалентного тела или тела вытеснения. За линией отрыва потока зад-

няя часть Е 2 поверхности эквивалентного тела, моделирующая влияние спутного следа, должна представлять собой продолжение поверхности ее передней части. Таким образом, поверхность эквивалентного тела Е формируется как объединение поверхностей Е1 и

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

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

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

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

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

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

дп ^ду ^ дм о дх ду дг

дп дп дп дп 1 др

--ъ п--ъ V--ъ м— =---;

дХ дх ду дг р дх

ду ду ду ду 1 др

--ъ п--ъ у--ъ м— =---;

дХ дх ду дг р ду

дм дм дм дм 1 др

--ъ п--ъ у--ъ м— =---,

дХ дх ду дг р дг

(1)

здесь п, у и м - координаты вектора скорости V = п ■ I + у ■ ] + м • к в декартовой системе координат; р и р - плотность и давление.

Граничными условиями являются условие непротекания поверхности эквивалентного тела Е и условие затухания возмущений на бесконечности. Первое состоит в том, что в произвольной точке М0(х0, у0, г0) на поверхности Е эквивалентного тела должна равняться нулю нормальная составляющая вектора скорости потока:

V,, (М0) = V (М0) • п(М0)\2 = 0, (2)

где п (М0)| 2 - вектор нормали к поверхности Е в точке М0.

Согласно условию затухания возмущений на бесконечности возмущения, вызванные телом, должны затухать, т. е. при М0 ^ 0 скорость V должна стремиться к скорости набегающего потока V0.

Предполагалось, что всюду вне поверхности эквивалентного тела течение потенциальное, и существует потенциал ф(М0) возмущенных скоростей, а скорость потока вычисляется по формуле

V (М 0) = + Уф( М 0), (3)

где V00 - вектор скорости набегающего потока; V - оператор Гамильтона; Уф(М0) - градиент потенциала возмущенных скоростей. Первое уравнение системы (1) может быть записано в виде

= 0.

Подставляя в это уравнение выражение (3) получим

ё^гаёф = 0,

т. е. уравнение Лапласа:

Лф = 0,

где Л - оператор Лапласа.

Подстановка выражения (3) в граничное условие непротекания (2) приводит к следующему соотношению на поверхности эквивалентного тела:

дф(М0)

дп

-V,- п (М 0)| 2. (4)

Учитывая граничные условия непротекания поверхности эквивалентного тела X, записанные в форме (4), и граничные условия затухания возмущений на бесконечности получаем, что потенциал возмущенных скоростей ф(М0) является решением внешней задачи Неймана для уравнения Лапласа:

Дф(М0) = 0, (5)

дф(М 0)

= -К0О- п (М0)| 2, (6)

I

дп

ф(М0) ^ 0, М0 ^ о, (7)

Уф(М0) ^ 0, М0 ^о. (8)

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

ф( М0) = -М ЩМЪ (М V X. (9)

4л X г

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

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

1 д г( Т - п(М)

^^1 g (М * ъ-у^-п (М0).

4-дпМ0 IV г

ъ

Для реализации численного метода поверхность эквивалентного тела аппроксимируется конечным числом панелей - многоугольников Ек, к = 1, ..., Ы, с постоянной плотностью потенциала двойного

слоя gk. Из аддитивности и линейности двойного интеграла и линейности градиента следует, что

7Т N Г 1 ^ Г Г• п(МЛ Л

V(М0) = gk —Ум0 ]-

к= I 4Л Ек Г

(10)

Градиент потенциала двойного слоя, размещенного на панели Ек с постоянной плотностью gk, равен скорости, индуцированной замкнутой вихревой нитью Ьк, расположенной на границе 5Ек и имеющей циркуляцию Гк, равную -gk, следовательно, справедливо следующее соотношение:

4л Мо Е г3 4л} г3

Ек Чк

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

Указанный вектор назовем функцией скорости и вычислим по закону Био - Савара:

1 ^ г Г • п (М) _ 1 г & х г

*мо)=-о г^е=-ч

3

4л V г3 4л; г

Ек Чк

после чего формулу (8) представим в виде

_ _ N _

V (М о) = ГkWk (М о).

к=1

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

Для определения неизвестных циркуляций Гк граничные условия непротекания поверхности эквивалентного тела Е удовлетворяются в контрольных точках Су, У = 1,..., N, расположенных в гео-

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

= (V(Cv) - к,) • n(Cv) = £ rkWk(Cv) • n(Cv) dn k=l

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

N ___

t Гк (Wk (Cv )• n (Cv)) = -V^ n (Cv),

k=1

где n (Cv) - орт нормали к панели Zk в контрольной точке Cv. Поэтому циркуляции определяются из системы линейных алгебраических уравнений

N

t «vk Г = bv, v = 1,..., N, (11)

k=1

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

avk = Wk (Cv )• n (Cv X bv=-Vo • n (Cv ).

Для замкнутых поверхностей t в методе дискретных вихрей матрица ||avk|| получается вырожденной, поэтому система (11) решается с использованием регуляризирующей переменной [9].

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

вне поверхности тела определяются из соотношения (11). Так как контрольные точки расположены на поверхности тела, то в них скорость как градиент потенциала двойного слоя терпит разрыв

(Cv) ^гт , dg(Cv )■

AV (Cv) = ) + ^-^T2(Cv),

причем частные производные вычисляются по направлениям двух ортогональных ортов т1(Су) и т2(Су), лежащих в касательной плоскости. Учитывая, что gV равняется - Г у, эти частные производные определяются численно по значениям циркуляций V-го и соседних с

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

_ _ N _ _

V (С) = ГЛ (Cv) + - АV (Cv).

к=1 2

Давление р определялось из интеграла Бернулли:

р = р«+ 2- р«И - V2). (12)

В соответствии с (12) безразмерные коэффициенты давления

Ср = 2(р - р« )/(рм^)

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

Ср = 1 - (V / V«, )2.

Результаты численных расчетов. По предложенной методике произведено моделирование обтекания осесимметричных тел различной формы. На данном этапе исследований было признано целесообразным осуществлять построение эквивалентного тела наиболее простой конфигурации. Поэтому поправка на толщину вытеснения пограничного слоя не производилась и передний участок поверхности эквивалентного тела совмещался с поверхностью обтекаемого тела. За линией отрыва потока, которая считалась совпадающей с линией донного среза, хвостовой участок поверхности эквивалентного тела, моделирующий влияние спутного следа, представлял собой продолжение поверхности обтекаемого тела. Длина указанного участка поверхности эквивалентного тела обозначалась как 1сл.

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

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

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

Для выявления особенностей расчетной схемы и отработки деталей алгоритма было проведено численное исследование обтекания цилиндрического тела, снабженного головной частью оживальной формы и хвостовой сужающейся частью конической формы. В качестве характерного линейного размера этого составного тела принимался диаметр миделевого сечения ё, т. е. диаметр цилиндрической части. Длина головной оживальной части составляла ё, средней цилиндрической - 2ё, а хвостовой части - 0,5ё. Полуугол сужения хвостовой конической части равнялся 10 град. В соответствии со сделанными предположениями для тела оживально-цилиндрическо-конической формы хвостовой участок эквивалентного тела строился как продолжение конической поверхности хвостовой части обтекаемого тела и имел полуугол сужения, равный 10 град. При такой конфигурации длина хвостового участка эквивалентного тела составила 1сл = 2,3 3ё.

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

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

Рис. 1. Зависимость наибольшего по абсолютному значению разрежения ср . на обтекаемом теле от суммарного числа панелей N

Для данного оживально-цилиндрическо-конического тела на рис. 2 представлена зависимость коэффициента давления ср на образующей обтекаемого тела от безразмерной продольной координаты х = х / ё, где х - абсцисса той точки поверхности, в которой определяется величина ср. Суммарное число панелей N равнялось 1 856.

Расчетное распределение коэффициента давления ср на поверхности

исследуемого оживально-цилиндрическо-конического тела (сплошная линия на рис. 2) согласуется с экспериментальными данными [10] (прямоугольники соответственно на рис. 2).

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

Большой практический интерес представляет исследование обтекания осесимметричных тел, не имеющих сужающихся или расширяющихся хвостовых частей. Поэтому на основе представленной методики было проведено численное исследование обтекания цилиндрических тел, снабженных головной частью оживальной формы, но не имеющих хвостовой части. В ходе вычислительного эксперимента варьировались длины носового оживального и основного цилиндрического участков поверхности обтекаемого тела. В качестве примера ниже представлены результаты математического моделирования осесимметричного обтекания тела, у которого длина головной оживальной части составляла ё, длина основной цилиндрической - 2ё, а общая длина тела 3ё. Сформулированным выше дополнительным условиям соответствует, например, выбор конфигурации заднего участка эквивалентного тела в виде сужающейся оживальной поверхности. На рис. 3 и 4 представлены соответственно зависимости, безразмерной скорости V = V / Уот и коэффициента давления с р на поверхности эквивалентного тела от безразмерной продольной координаты х и длины 1сл хвостового участка эквивалентного тела. Значениям 1сл, равным 2ё, 2,5ё, 3ё и 4ё на рис. 3

соответствуют линии 1, 2, 3 и 4.

Анализ представленных зависимостей показал, что на большей части поверхности обтекаемого тела, вплоть до значений безразмерной продольной координаты х < 2,5, характер распределения безразмерной скорости и коэффициента давления слабо зависит от длины 1сл во всем исследованном диапазоне этого параметра. Если же длина хвостового участка эквивалентного тела превышает величину 4ё, то область влияния параметра 1сл на распределение безразмерной скорости и коэффициента давления вдоль образующей обтекаемого тела уменьшается еще сильнее и это влияние проявляется лишь непосредственно в окрестности донного среза.

Г ч___

/ 1 '/// 2 3 4

/

О 1 2 3 4 х

Рис. 3. Зависимость безразмерной скорости V на поверхности эквивалентного тела от безразмерной продольной координаты х и длины 4л хвостового участка эквивалентного тела:

1 - 1сл = га ; 2 - 2,5а; 3 - за; 4 - ла

I

\

\

\;

1 2 3 4

О 1 2 3 4 X

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

1 - 1сл = га; 2 - 2,5а; 3 - за; 4 - ла

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

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

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

Для фюзеляжей с плоским донным срезом при Яе = 5 • 106 в работе [6] приведены экспериментальные значения аэродинамического коэффициента донного сопротивления, лежащие в диапазоне 0,09...0,13. Если отбросить крайние из представленных значений, можно ориентироваться на величину коэффициента донного сопротивления приблизительно равную 0,1, которой в свою очередь соответствует значение коэффициента донного давления сРдон = _0,1

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

принимает значение, равное коэффициенту донного давления ср . Поэтому зависимости ср = ср (х, 1сл), аналогичные представленной

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

тела реализовалось бы равенство ср = срдон.

Для иллюстрации полученных результатов можно указать, например, что при сохранении длины головной части, равной ё, варьирование длины цилиндрического участка в пределах 3ё ...5ё приводит к изменению длины 1сл заднего участка эквивалентного тела в диапазоне 2,3ё ...2,6ё.

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

ЛИТЕРАТУРА

[1] Андерсон Д., Таннехил Дж., Плетчер Р. Вычислительная гидромеханика и теплообмен. В 2 т. Москва, Мир, 1990. Т. 1, 384 с., Т. 2, 392 с.

[2] Флетчер К. Вычислительные методы в динамике жидкости. В 2 т. Москва, Мир, 1991. Т. 1, 504 с., Т. 2, 552 с.

[3] Чжен П. Отрывные течения. В 3 т. Москва, Мир, 1972. Т. 1, 299 с., Т. 2, 280 с., Т. 3, 333 с.

[4] Четверушкин Б.Н., Шильников Е.В. Вычислительный и программный инструментарий для моделирования трехмерных течений вязкого газа на многопроцессорных системах. Журнал вычислительной математики и математической физики, 2008, т. 48, № 2, с. 309-320.

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

[6] Тимофеев В.Н., Бушуев А.Ю. Математическое моделирование обтекания тел с отрывом дозвукового потока в донной области. Вестник МГТУ им Н.Э. Баумана. Сер. Естественные науки. Спец. Вып. № 4. Математическое моделирование. 2012, с. 94-99.

[7] Девнин С.И. Аэрогидромеханика плохообтекаемых конструкций. Справочник. Ленинград, Судостроение, 1983, 320 с.

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

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

[10] Woodward F.A. An Improved Method for the Aerodynamic Analysis of Wing-Body-Tail Configurations in Subsonic and Supersonic Flow. Washington,1973, NASA CR-2228, рр. 1-130.

Статья поступила в редакцию 27.06.2013

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

Тимофеев В.Н., Бушуев А.Ю. Численное моделирование дозвукового отрывного обтекания осесимметричных тел. Инженерный журнал: наука и инновации, 2013, вып. 9. URL: http://engjournal.ru/catalog/mathmodel/aero/967.html

Тимофеев Валерий Николаевич окончил МВТУ им. Н.Э. Баумана в 1974 г. Канд. техн. наук, доцент кафедры «Вычислительная математика и математическая физика». Автор более 40 научных и методических работ. Область научных интересов: математическое моделирование, численные методы, механика жидкости и газа, аэродинамика, методы оптимизации. e-mail: v_n_1951@mail.ru

Бушуев Александр Юрьевич окончил МВТУ им. Н. Э. Баумана в 1974 г. Канд. техн. наук, доцент кафедры «Вычислительная математика и математическая физика». Автор 20 научных работ. Область научных интересов: математическое моделирование в технике, методы оптимизации и принятия решений, численные методы. A.Ju.Bushuev@ya.ru