CC BY
124
УЧЕНЫЕ ЗАПИСКИ Ц А Г И
Т о м IV 1 9 7 3 М3
УДК 533.6.013.2.011
РАСЧЕТ ОБТЕКАНИЯ ОСЕСИММЕТРИЧНЫХ ТЕЛ И НЕСУЩИХ КРЫЛОВЫХ ПРОФИЛЕЙ ТРАНСЗВУКОВЫМ ПОТОКОМ ГАЗА
А. С. Фонарев
Рассматривается задача обтекания осесимметричных тел и несущих крыловых профилей трансзвуковым потоком газа, численные решения выполнены с использованием конечноразностного метода Годунова. Приводятся результаты расчета обтекания осесимметричных тел некоторых классов (цилиндрическое тело со сферическим затуплением, плоский торец цилиндра), а также двух профилей (КАСА 64А-410 и пикообразного) под углами атаки в диапазоне чисел МОС = 0,7-1,3.
Применение численных методов установления к проблеме обтекания крыловых профилей и осесимметричных тел трансзву^ ковым потоком газа позволило решить некоторые задачи в этой области. В работах [1] и [2] с использованием модификации конечноразностной схемы Бабенко, Воскресенского, Любимова и Русанова решены задачи обтекания заостренного симметричного профиля, составленного из дуг окружности (под нулевым углом атаки), и осесимметричного обтекания тел типа сферы и эллипсоида. В работах [3—5] методом Годунова рассмотрена задача дифракции ударной волны на симметричном профиле, сфере и цилиндре с последующим установлением стационарного трансзвукового обтекания. Для решения задач трансзвукового обтекания успешно использовались и другие численные методы — метод крупных частиц и различные модификации метода Лакса — Вендрова [6—8]. В работе [9] с использованием метода Лакса —Вендрова впервые была решена задача обтекания несимметричного профиля под углом атаки. Наряду с методами установления в последнее время стали интенсивно развиваться итерационные методы расчета стационарного обтекания трансзвукового профиля в рамках потенциальной теории малых возмущений [10, 11].
Ниже рассматривается задача обтекания осесимметричных тел и несущего профиля крыла (под углом атаки) трансзвуковым потоком газа, численные решения выполнены конечноразностным методом Годунова. В качестве примера приведены результаты расчетов
обтекания осесимметричных тел некоторых классов — цилиндрического тела со сферическим затуплением, плоского торца цилиндра, а также двух профилей — NACA 64А-410 и пикообразного — под углом атаки в диапазоне чисел Мда = 0,7 1,3.
Постановка задачи. Решение стационарной задачи обтекания будем получать как предельное из следующей нестационарной задачи. Пусть в начальный момент времени неподвижное тело начинает внезапно двигаться с постоянной скоростью, соответствующей заданному числу М0 . Начало прямоугольной системы координат ху поместим в носок рассматриваемого тела, ось х направим по оси симметрии тела или по хорде профиля. Рассмотрим плоский или осесимметричный случай обтекания тела невязким нетеилопроводным совершенным газом с показателем адиабаты 7=1,4, относительно формы контура тела никаких ограничений делать не будем.
В нестационарном течении и предельном установившемся обтекании присутствуют ударные волны, поэтому будем исходить из интегральных законов сохранения для нестационарных уравнений Эйлера [12]. В качестве граничных условий зададим на поверхности тела обычные условия непротекания: нормальная к поверхности тела составляющая вектора скорости равна нулю. В начальный момент времени (поскольку движение начинается внезапно) во всем набегающем на профиль потоке газа параметры постоянны:
р=рл, р = рзо, и — Macacos а, V = Мао Sin а.
Здесь р — давление, р— плотность, и, v — составляющие вектора скорости в направлении осей л и у соответственно, а — скорость звука, а — угол атаки профиля (для осесимметричных тел рассматривается случай а = 0); индекс „оо“ соответствует параметрам набегающего потока.
В качестве характерного линейного размера /0 примем в плоском случае хорду профиля, в осесимметричном — радиус поперечного сечения цилиндрической части тела. Приводимые ниже численные результаты даны в безразмерном виде, в качестве характерных величин приняты параметры невозмущенного потока рм, Роо и характерная длина /0-
Метод численного расчета. Обтекание тел трансзвуковым потоком характеризуется наличием местных сверхзвуковых зон, замыкающихся, как правило, скачками уплотнения. При этом положение звуковой линии и скачка уплотнения, разделяющих дозвуковую и сверхзвуковую области потока, заранее неизвестны. Применим в связи с этим конечноразностный метод С. К. Годунова [12], позволяющий единым образом рассматривать дозвуковые, сверхзвуковые и смешанные нестационарные течения без специального выделения сильных разрывов; опыт расчета [3—5] уже доказал его эффективность применительно к решению такого рода задач.
Так как никаких ограничений на форму плоского или осесимметричного тела не накладывается, то для решения задачи должна использоваться расчетная сетка, легко приспособляемая к различным формам тел. Использовались разнообразные типы сеток. Для простых тел, описываемых одним характерным размером, например сферы, куба и т. д., сетка состояла из радиальных линий и концентрических окружностей, квадратов и т. д. (фиг. 1, а). Для удлиненных тел типа профиля (при небольших изменениях угла наклона тела) одно семейство линий сетки состояло из параллель-
ных вертикальных линий, второе семейство по-прежнему повторяло форму контура (фиг. 1, б). И, наконец, для тел более сложной формы сетка могла быть комбинированной (фиг. 1,в). Как правило, величина поперечного шага расчетной сетки при удалении от тела увеличивается (по закону геометрической прогрессии), что дает возможность последнюю линию (внешнюю границу поля) расположить достаточно далеко от тела, в области практически невозмущенного потока, чтобы при решении задачи на ней можно было
^ов ;.0! 0^7 а,в? о,э±
задать либо условия невозмущенного потока, либо равенство нулю градиентов всех параметров газа в направлении от тела вдоль линий первого семейства [3].
Конечноразностная схема расчета строится аналогично приведенной в работе [12]; расчет параметров распада разрыва проводится либо по точным итерационным соотношениям [13], либо по линеаризованным, в зависимости от величины градиентов соответствующих функций, при этом тип расчета выбирается ЭВМ автоматически.
Обтекание осесимметричных тел. Численное решение проводилось ДЛЯ двух ТИПОВ тел конеч- >,08^0^93^85 0^5
ных размеров: цилиндра со сферическим затуплением радиусом 1 (длина цилиндрической части £=3) и торца цилиндра с плоским срезом (длина цилиндрической части £ = 4). Кормовая часть обоих тел представляет собой плоский срез.
Рассмотрены три режима обтекания: Мэо = 0,8; 1,0; 1,3.
На фиг. 2, а, б показаны линии равных значений давления около рассматриваемых тел (иллюстрируется только случай, когда Моо = 1).
Давление отнесено к статическому давлению невозмущенного по-
тока. Хорошо видно, как происходит разгон потока от условий торможения в критической точке до образования сверхзвуковой области (звуковая линия обозначена пунктиром) с возникновением скачка уплотнения. При этом в случае обтекания торца цилиндра резкое торможение потока (типа скачка) происходит сразу же после расширения в угловой точке, дальнейшее торможение потока происходит довольно плавно (это хорошо видно на фиг. 4).
Во всех рассмотренных случаях в кормовой части тел за срезом образуется область возвратного течения. Типичная картина линий тока приведена на фиг. 2, в (Моо=1)- Штрих-пунктирная линия разделяет основной поток и зону возвратного течения.
Сделаем несколько замечаний относительно образования зон возвратного течения при численном решении уравнений Эйлера.
Эти зоны достаточно устойчивы и появляются независимо от способа установления потока (например, в случае разгона с ускорением, при дифракции ударной волны с последующим установлением обтекания). Зоны возвратного течения наблюдаются и в случае применения других численных методов, например метода „крупных частиц" |6]. Опыт расчета показал, что они образуются в сверхзвуковых и смешанных течениях (там, где есть градиенты энтропии); в чисто дозвуковом потенциальном потоке реализуется безотрывное обтекание. Как правило, точность определения параметров в зонах возвратного течения существенно ниже, чем в основном потоке; обычно при точности решения в основном потоке 1—2% в зонах возвратного течения ошибка в значении константы Бернулли достигает 5—10%.
Эти зоны качественно похожи на реальные срывные зоны в донных областях, их величина зависит от параметров потока. Так, для рассматриваемых тел было обнаружено, что они увеличиваются с ростом числа Мсо. В таблице приведены значения расстояния по оси л: от донного среза до точки разделения потоков (в долях радиуса цилиндрической части тела). Из этих данных следует, что форма носовой части рассматриваемых тел практически не влияет на размеры зоны возвратного течения.
Образование зон возвратного течения при численном решении уравнений Эйлера представляет собой недавно обнаруженное, вообще говоря, парадоксальное явление, которое обычно связывают с присутствием искусственной, или схемной, вязкости в конечноразностных схемах.
На фиг. 3 показано распределение коэффициента давления
р = (р — Рх>) I 7М_0р<х по передней и цилиндрической части тела
со сферическим затуплением; по оси х отложена продольная координата. Кружками нанесены экспериментальные данные, полученные в ЦАГИ. Совпадение экспериментальных и теоретических значений достаточно хорошее, особенно в передней части тела, где происходит интенсивный разгон потока и возникают
Форма тела Мсо Ах
Цилиндр со сфериче- 0.8 1,35
ским затуплением 1,0 1,53
1.3 1.80
Цилиндр с плоским 0.8 1.4!
торцом 0,95 1,56
1,3 1,80
большие аэродинамические силы.^В области скачка уплотнения, замыкающего сверхзвуковую зону, наблюдается довольно заметное расхождение численных и экспериментальных данных, отмечаемое авторами всех численных исследований трансзвукового обтекания тел. Оно обусловлено существенным проявлением эффекта вязкости при взаимодействии скачка уплотнения с пограничным слоем. В случае сверхзвукового потока (Моо=1,3) результаты численных расчетов хорошо согласуются с данными эксперимента.
, Мл О,* Л О 1 1
О \ о о Аі
[ V-" Т-п а V / -о
<7 1 І /7 і
1 п '
\_li_
Фиг. 3
М -. . . V 1 /*«,= 0,35 1,3
V \ \ \
\ V. V Ч—"
<7 0 / і 3 х [ I !
^4 і I
а . . . •Л
Фиг. 4
На фиг. 4 показано распределение коэффициента давления по контуру торца цилиндра, по оси х—отложено расстояние вдоль контура тела. Кривые имеют довольно сложный вид. Пики давления соответствуют точкам излома. Сверхзвуковая область начинается сразу за передней угловой точкой после сильного расширения потока (веер расширения хорошо виден на фиг. 2, б). Кружками нанесены экспериментальные данные, полученные в ЦАГИ для торца цилиндра.
Трансзвуковое обтекание профиля под углом атаки. Будем рассматривать профили с закругленной передней и острой задней кромками. Расчетная сетка, сохраняя принципиально основные черты осесимметричного случая, имеет особенности, связанные с необходимостью расчета как верхней, так и нижней половины, поэтому приведем ее описание.
Разобьем все поле расчета вертикальной линией на две части: область/(см. фиг. 1, г), прилежащую к носку и характеризующуюся большими значениями кривизны профиля, и область II, в которой наклоны образующей профиля невелики. В области / контур профиля разобьем на частей и проведем из точек разбиения лучи с равномерным шагом по углу Дб = 1с/А1. Разбив верхнюю и нижнюю стороны оставшейся части профиля и некоторый отрезок оси * позади профиля соответственно на /г4 частей, проведем
через точки деления вертикальные лучи. Перенумеруем их (первое семейство линий), начиная с крайнего правого луча п — О, идущего вниз по оси х, обходя профиль по часовой стрелке, до крайнего правого луча п = N, уходящего вверх от оси л (из каждой точки разбиения оси х выходит два луча — вниз и вверх). Таким образом, .V = А, + £2 -)- к3 + 2&4. На первом луче (п = 0) отложим последовательно М отрезков, длины которых 1т увеличиваются по закону геометрической прогрессии: 1т — 1т-(т = 2, 3, 4, где д — знаменатель прогрессии. Ломаную линию, сое-
диняющую начальные точки всех лучей, расположенные на оси х и контуре тела, назовем нулевой линией т = 0 второго семейства. Остальные линии второго семейства проведем через точки разбиения первого луча параллельно соответствующим отрезкам нулевой линии. Линия т — М и лучи п = 0 и n — N являются внешними границами расчетной области. Удаленность внешней границы от контура профиля определяется выбором величин /, и <7, а также длиной отрезка оси позади профиля. Предполагая, что граница находится достаточно далеко от тела, так что возмущения, пришедшие к ней от профиля, малы и практически не влияют на поле течения около профиля, зададим условия на границе поля в виде равенства нулю производных газодинамических параметров [3] вдоль лучей, идущих от тела. Ячейки, прилегающие к контуру тела, будем рассчитывать с учетом выполнения условий непроте-кания на теле, а ячейки, примыкающие к оси х позади профиля, будем рассматривать как обычные внутренние ячейки.
По описанной выше методике было рассчитано течение около профилей двух типов—ЫАСА 64А-410 с относительной толщиной 10% и пикообразного профиля, координаты которого задаются формулой [14]
V = —с ) х (1 — х)
4/3
(с — относительная толщина профиля), с расположением максимальной толщины на 1/3 хорды.
При проведении расчетов были выбраны следующие параметры расчетной сетки. Область / с переменным направлением лучей занимала 2,5% хорды в носовой части профиля; число секторов было равно &! — 21. В остальной части течения лучи строились вертикально. Верхняя и нижняя части профиля разбивались на й2 = й3 = 23 части. Число отрезков между лучами в области за профилем &4 = 20. Крайний луч располагался на расстоянии, равном
примерно 3 от задней кромки профиля, шаг между лучами в этой области увеличивался по мере удаления от профиля по арифметической прогрессии, начинающейся с шага 0,025. На поверхности профиля шаг сетки выбирался переменным (максимальное значение шага — 0,05, минимальное ~ 0,0018), обеспечивалось по возможности плавное изменение величины шага при переходе из одной области в другую.
Второе семейство расчетной сетки состояло из 38 линий. Начальный шаг разбиения по вертикали (на луче п = 0, см. фиг. 1, г) был равен /, = 0,002, знаменатель прогрессии <7 = 0,18-т-1,2. Выбранные значения этих величин обеспечивали удаление внешней границы расчетной сетки на расстояние примерно 5—6. Можно ожидать (и это подтвердили расчеты), что при таком расположении внешних границ параметры течения вблизи профиля будут мало отличаться от параметров в случае обтекания безграничным потоком.
Сильная неравномерность шага сетки связана с необходимостью расчета областей профиля с большой и малой кривизной (соответственно с резким и плавным изменением параметров потока). Отметим, что область резкого изменения параметров в носке профиля занимает очень малую часть хорды (—2,5%), тем не менее от точности расчета течения в этой области зависит правильность решения задачи на остальной части профиля. Можно считать, что отношение минимального и максимального шагов сетки на профиле определяется, в основном, радиусом кривизны носка и длиной хорды профиля. Шаг по времени определяется минимальным координатным шагом, и из-за неравномерности координатных шагов он в значительной части течения (~95% хорды) существенно меньше допустимого (для рассмотренных случаев отношение максимального локального шага по времени к минимальному равно примерно 25). Это, в свою очередь, приводит к резкому увеличению потребного машинного времени, в несколько раз превышающего время расчета аналогичной задачи с одним характерным размером (например, трансзвукового обтекания цилиндра). Наблюдалось также заметное увеличение времени установления в трансзвуковом режиме по сравнению со сверхзвуковым.
Результаты расчета приведены на фиг. 5—7.
На фиг. 5, а изображены линии равных давлений в поле течения около профиля ЫАСА 64А-410 при М,;г =0,72, а = 4°, На нижней поверхности профиля реализуется сжатие потока, на верхней — сильное разрежение с образованием сверхзвуковой области, граница которой нанесена пунктирной линией. Сверхзвуковая область замыкается скачком уплотнения, представленным на графике в виде сгущения изобар.
Распределение давления по поверхности профиля показано на фиг. 5,6: сплошная линия—данные настоящей работы, пунктирная — численные результаты из работы [9], кружками нанесены данные, полученные экспериментально [15]. Кривые давления, соответствующие нижней стороне поверхности профиля, хорошо согласуются одна с другой. На верхней стороне профиля наблюдается заметное различие результатов в области расположения скачка уплотнения, особенно различаются экспериментальная и теоретические кривые. Как уже отмечалось, расхождение (теории и эксперимента) в положении скачков уплотнения наблюдается также при расчете другими методами на других телах [1, 3].
р-'О^б О'в 0,75 0,3 0,16 09!
При решении невязкой задачи обтекания профиля обычно
накладывают условие Жуковского о расположении задней критической точки на острой кромке профиля. В численных расчетах данной работы это условие специально не задавалось, тем не менее при небольших углах атаки получалось течение с расположением критической точки на острой кромке, что, по-видимому,
объясняется влиянием „схемной вязкости11, присутствующей в используемом конечноразностном методе.
При увеличении угла атаки характер обтекания изменяется: на верхней стороне профиля в центральной его части образуется область возвратного течения, подобного описанному выше течению в кормовой части осесимметричных тел;при этом перетекания потока с нижней стороны на верхнюю по-прежнему не происходит. Значение угла атаки, при котором происходит переход от одного режима обтекания к другому на рассматриваемом профиле, специально не определялось. Отметим лишь, что при а = 8° устанавливается обтекание без возвратной зоны, при а = '24° эта зона занимает значительную часть верхней поверхности профиля.
03
1.0
1.2
1А
4 \ \ \ \ 1 \ 1
ш Г о » \ ' 1 V \ \
г
1
0.2
0,4 0,5
0)
Фиг. 5
03
ка х
Фиг. 6
Приведем результаты расчетов обтекания описанного выше пикообразного профиля толщиной с = 0,1 при числе Моо=1.
Этот профиль подробно исследован в работе [14], где получено приближенное решение для распределения давления с учетом вязкости методом, аналогичным методу локальной линеаризации Спрейтера [16]. Кроме того, для него имеются экспериментальные результаты. На фиг. 6 представлены линии равных значений чисел М при угле атаки профиля а = 2°. Видно, что разгон потока от критической точки вверх по поверхности профиля происходит весьма интенсивно, так что звуковая линия, обозначенная пунктиром, располагается вблизи передней кромки. В силу симметрии профиля и малости угла атаки на нижней стороне профиля поток также разгоняется. Скачок уплотнения, замыкающий сверхзвуковую зону на верхней и нижней поверхностях, располагается на задней кромке, что согласуется с приближенными расчетными и экспериментальными данными.
На фиг. 7 представлена зависимость коэффициента давления от продольной координаты профиля при угле атаки, равном нулю. Сравниваются три кривые: сплошная кривая — численный расчет, светлые кружки — экспериментальные данные [17] для Моо=1,01 и пунктирная линия — приближенное аналитическое решение [14J; следует отметить неплохое согласование как теоретических данных между собой, так и теоретических данных с экспериментальными.
Там же показано распределение коэффициента давления по хорде профиля при а = 2° и Моо=1. Кривая с обозначением а —— 2° соответствует нижней поверхности. Экспериментальные данные (отмечены черными кружками) взяты из работы [17]. Как видно из приведенных материалов, расчетные результаты удовлетворительно согласуются с экспериментальными.
В заключение автор приносит благодарность В. В. Сычеву и А. И. Голубинскому за полезное обсуждение результатов работы и В. П. Колгану за помощь при проведении некоторых численных расчетов.
ЛИТЕРАТУРА
1. Лифшиц Ю. Б. О расчете трансзвукового обтекания симметричного профиля в свободной струе. „Изв. АН СССР — МЖГ“,
1969, № 1.
2. Л и п н и ц к и й Ю. М., Лифшиц Ю. Б. О расчете обтекания тел вращения трансзвуковым потоком. ПММ, т. 34, № 3, 1970.
<1 V 3. Фонарев А. С. Расчет дифракции ударной волны на профиле с последующим установлением стационарного сверхзвукового и трансзвукового обтекания. „Изв. АН СССР — МЖГ“, 1971, № 4. i. 'J 4. Фонарев А. С., К о л г а й В. П. Численный расчет дифракции ударной волны на сфере и цилиндре и установление стационарного обтекания. Труды ЦАГИ, вып. 1324, 1971.
5. К о л г а н В. П., Фонарев А. С. Установление стационарного обтекания при набегании ударной волны на сферу и цилиндр. .Изв.
АН СССР — МЖГ-, 1972, № 5.
6. Белоцерковский О. М., Давыдов Ю. М. Нестационарный метод крупных частиц для газодинамических расчетов. „Журн. вычисл. матем. и матем. физ.“, т. II, № 1, 1971.
7. С г о s s in а п В., М о г е 11 у Q. Time dependent calculation of transonic flows. AIAA Paper. No 70-1322.
8. Kentzer C. P. Computations of time dependent flows on an in finite domain. AIAA Paper, No 70-45.
9. Magnus R„ Yoshihara H. Inviscid transonic flow over airfoils. AIAA Paper, No 70-47.
10. Mur man E. М., Cole J. D. Calculation of plane steady transonic flows. AIAA Paper, No 70-IS8.
11. Steger J. L., Lomax H. Numerical calculation of transonic flow about two-dimensional airfoils by relaxation procedures. AIAA Paper, No 71-669.
12. Годунов С. К,, 3 а б.р о д и я А. В., Прокопов Г. П. Разностная схема для двумерных нестационарных задач газовой динамики и расчет обтекания с отошедшей ударной волной. „Журн. вычисл. матем. и математической физики*, т. 1, № 6. 1961.
13. Годунов С. К. Разностный метод численного расчета разрывных решений уравнений гидродинамики. .Математический сборник”, 47, вып. 3, 1959.
14. Truit R. Shockless transonic airfoils. AIAA Paper, No 70-187.
15. Stivers L. S. Effects of subsonic Mach number on the forces
and pressure distributions on four NACA 64A-Series airfoil sections at angles of attack as high as 28°. NACA TN 3162, 1954.
16. S p r e i t e r J. R., A 1 k s n e, A I b e r t a J. Thin airfoil theory based
on approximate solution of the transonic flow' equation. NASA TN 3970,1957.
17. S m e i a n a F. О., К n e p p e r D. P. Predicted and measured pressure distributions on two airfoils at Mach number 1. AROD-8130, North
Carolina State University.
Рукопись поступила 1 iIX 1972
