CC BY
40
Том XX
УЧЕНЫЕ ЗАПИСКИ ЦАГИ 1989
№ 1
УДК 533.6.011.5 : 532.582.2
ЧИСЛЕННОЕ РЕШЕНИЕ ЗАДАЧИ ОБТЕКАНИЯ ТРАПЕЦИЕВИДНОГО КЛИНА СВЕРХЗВУКОВЫМ ПОТОКОМ ИДЕАЛЬНОГО ГАЗА
С. М. Босняков, В. В. Коваленко, С. В. Михайлов, Н. X. Ремеев
Предложенным ранее методом [1, 2] решена задача обтекания трапециевидного клина сверхзвуковым потоком идеального газа с использованием конечно-разностной схемы первого порядка аппроксимации [3]. Проведено сравнение результатов расчета, полученных в рамках этого метода, с результатами,, полученными с использованием схемы второго порядка аппроксимации [4] й результатами проведенного авторами эксперимента.
В работе [4] разработан метод решения задачи обтекания трапециевидного клина и проведено сравнение различных подходов к решению поставленной задачи [5 — 7]. Показано преимущество метода второго порядка аппроксимации. При этом основное внимание уделено течению у наветренной поверхности клина. На практике важно знать распределение параметров также у боковой и подветренной поверхностей, что требует выделения в расчете нижней и верхней боковых кромок клина. Нижняя кромка в методе [4] выделена путем фиксации числа узлов у наветренной поверхности, верхняя кромка в процессе счета не выделяется. Это приводит к изменению числа узлов, приходящихся соответственно на подветренную и боковую поверхности с изменением соотношения h/b, где Л —расстояние между нижней и верхней кромками, b — ширина клина. Опыт [5 — 7] показывает, что подвижная сетка в окрестности острых кромок дает забросы в решении. Особенно чувствительны к этому схемы повышенного порядка аппроксимации.
Используемый метод [1, 2] базируется на идее разбиения расчетной области на подобласти, предложенной в работах [8—11]. Впервые авторы реализовали его применительно к задаче обтекания плоского сверхзвукового воздухозаборника [1] и компоновки крыла с фюзеляжем [2]. В данной статье рассмотрен вопрос точности расчета обтекания трапециевидного клина.
1. Рассчитывается пространственное обтекание клина трапециевидной формы в плане сверхзвуковым потоком невязкого газа. Клин располагается в декартовой системе координат таким образом, что его передняя кромка совпадает с плоскостью Х=0, а строительная горизонталь — с плоскостью У = const. Для поверхностей клина введены условные названия: 1) наветренная (нижняя); 2) боковая; 3) подветренная (верхняя). Угол наклона
MD
Рис. 1
наветренной поверхности к плоскости Y= const обозначается бь угол наклона боковой поверхности к плоскости Z= const обозначается 62, (рис. 1).
2. Решается краевая задача для полной системы безразмерных уравнений Эйлера в декартовой системе координат, дополненных уравнением Бернулли:
(р + ры2) + (quv) + -Jj (р uw) = 0;
17 + її + pt^ + 17 = 0;
(р uw) + JL (ру ш) _|_ JL (р _|_ р w2) = 0;
-t№+-k{pv) + ~t{pw) = °’
и -\-V -\-W2 4-
27
7—1 р
7= 1,4.
const;
В Качестве характерных величин выбраны: а* — скорость звука в критическом сечении; р,*, — плотность набегающего потока, L—линейный размер, равный полуширине клина в сечении Х=0. Давление отнесено к величине р^а*. Система (1) предполагается Л-гиперболической и решается в области, ограниченной поверхностью головной ударной волны. На границах расчетной области выполняются условия: 1) непротекания на твердых поверхностях и плоскости симметрии; 2) Рэнкина — Гюгонио на головной волне. Для расчета используется маршевая процедура. Решение на слое X АХ = const получается по известному решению на слое X = const с привлечением законов сохранения массы, импульса, энергии и условий на границах расчетной области. Шаг АХ выбирается исходя из условия-Куранта — Фридрихса — Леви.
3. Метод [1,2] использует конечно-разностную схему С. К. Годунова [3]. В случае, клина конечной ширины расчетная область в некотором сечении X = const представляет собой криволинейный четырехугольник, одна из сто-' рои которого имеет угловые точки, см. рис. 1. Эта область разбивается на пять элементарных подобластей, в трех из которых строится сетка, близкая к ортогональной, а в двух — веерная сетка, стягивающаяся в точку. Таким образом обеспечивается сгущение узлов сетки в областях высоких градиентов параметров потока и фиксируется постоянное число узлов у каждой стороны клина. Алгоритм построения сетки основан на принципе прямого и обратного отображения криволинейного четырехугольника на некоторый прямоугольник с разбиением последнего на ячейки. Треугольная подобласть рассматривается как частный случай четырехугольной. В каждой подобласти представляются граничные условия и решается краевая задача. По завершении расчетного шага границы подобластей стыкуются, на них обеспечивается непрерывность (с точностью метода [3]) рассчитываемых функций.
Описанный метод реализован в программе для ЭВМ БЭСМ-6. При разработке программы проведена формализация задачи. Для описания структуры расчетной области в ОЗУ ЭВМ введены дополнительные данные о взаимном расположении подобластей, граничных условиях в каждой подобласти, о размерности сетки. С целью адаптации программы к особенностям течения предусмотрена возможность перестройки структуры расчетной области. Под перестройкой структуры понимается изменение типов граничных условий на гранях, изменение количества и взаимного расположения подобластей. Задание начальных данных проводится в рамках двумерных и конических приближений. Такой подход позволяет использовать программу для расчета обтекания широкого класса различных тел.
С целью ускорения счета в программе предусмотрена возможность автоматического объединения наиболее мелких ячеек сетки по заранее заданному критерию малости. Объединение осуществляется автоматически с сохранением всех потоков через грани объединенной ячейки и осреднением параметров внутри нее. При расчете клина процедура объединения реализуется при X близких к нулю и в угловых подобластях.
4. Метод [4] использует конечно-разностную схему Мак-Кормака [12, 13] второго порядка аппроксимации. Расчетная сетка строится с использованием ряда последовательных преобразований Кармана — Треффтца [14]. Это позволяет избавиться от угловых точек на контуре клина в сечениях X = const. В преобразованной плоскости осуществляется переход к полярным координатам и строится центрированная сетка с равномерным разбиением. Все особые точки преобразования [14] размещаются внутри контура клина и якобиан преобразования не имеет особенностей в расчётной области. Однозначность преобразования обеспечивается проведением разрезов, соединяющих пары сопряженных особых точек таким образом, чтобы разрезы проходили внутри клина (см. подробно [4]). Начальные данные задаются в некотором сечении X — const либо путем предварительно установления конического течения, либо заданием параметров невозмущенного потока. Результаты счета обрабатываются монотонизатором [15]. У нижней боковой кромки
клина реализуется алгоритм образования тангенциального разрыва, созданный в соответствии с [16]. Полученная обратным преобразованием сетка представлена на рис. 1.
5. Экспериментальное исследование проведено на модели трапециевидного клина в аэродинамической трубе с размерами рабочей части 500 X 500 мм2. Модель состояла из трех частей: 1) основного прямоугольного в плане клина шириной 100 мм, длиной 220 мм, углом между верхней и нижней поверхностями бкл=14°; 2) двух съемных боковых пирамидальных накладок с углами 61=14°, 62=14°. Клин и одна из накладок дренированы (нижняя и боковая поверхности). В аэродинамической трубе модель установлена на хвостовой державке, закрепленной на стойке сайр механизмами. Для учета влияния прогиба державки проведены контрольные нагружения до 120 кгс и соответствующие поправки на угол Да внесены в результаты эксперимента. Давление измерено в девяноста точках вдоль поверхности клина и в поперечных сечениях ^=0,77; 1,45; 2/7. Абсолютная величина погрешности в измерении статического давления составляла ±4Х Ю2 Па.
6. При обтекании клина на его поверхности можно выделить три области: в области / реализуется двумерное, в области II—коническое, а в области III—пространственное течение. Границами областей являются характеристические конуса от боковых кромок клина. Углы раствора конусов у наветренной и подветренной поверхностей различны и определяются значениями чисел Мь рассчитанными в двумерном приближении у каждой поверхности. Схема расположения областей течения представлена на рис. 1, где изображены подветренная поверхность и разворот боковых поверхностей. Область конического течения II заштрихована. Указанные особенности могут быть использованы для анализа точности расчета.
7. На рис. 2—4 сопоставлены распределения статического давления по поверхности клина в сечениях Л" = 1; 2; 4: Расчет методом [1, 2] проведен на сетках с числом узлов МХМ=9Х48 и Л1Х#= 18X85 (М—число узлов в направлении тело—скачок, а N—вдоль поверхности тела). Расчет методом [4] проведен соответственно на сетках 9X34 и 18X68.
Распределение давления по наветренной и подветренной поверхностям клина представлено вдоль координаты £ = (г^го) / (х\.%Ьъ), где 2о=\. Таким образом, точка |=1 соответствует боковой кромке трапециевидного клина, а точка £кон — границе между областями двумерного и конического течения. На рис. 2 граница отмечена вертикальной чертой и горизонтальными стрелками. Сопоставление эпюр статического давления на наветренной поверхности (см. рис. 2) показывают, что оба метода заметно размазывают эту границу, но при помощи метода второго порядка давление в ее окрестности определяется приблизительно на 1,5% точнее. В области / значения статического давления получены практически точно. Этот результат уже отмечен в [4 — 6]. При Екои<£<1 решение не должно зависеть от значения X. В расчете отличие результатов в сечениях X— 1 и Х=2 для обоих методов достигает 3%. Различие объясняется отсутствием аппроксимации конического течения при малых X и размытием границы областей / и II. Кроме того, метод [4] дает заметные осцилляции решения в окрестности боковой кромки (£= 1). Экспериментальные данные находятся в хорошем соответствии с результатами расчета^ см. рис. 2. Наибольшее отличие «5% наблюдается в окрестности боковой кромки клина и на границе областей двумерного, конического и пространственного течений.
На боковой поверхности клина при Л^4 реализуется чисто коническое течение. Обтекание нижней боковой кромки происходит с разрежением и возможно с отрывом потока. В сечении ^=1 метод [1—2] в основном формирует профиль решения (см. рис. 3), но понижение давления в окрестности нижней боковой кромки не проявляется (0<|<0,2). С увеличением значения координаты X этот недостаток постепенно устраняется на мелкой сетке.'
Рис. 2
Бо к о fast поверхность М-3 а = 5 ,
Рис. 3
Рис. 4
В методе [4] расчет нижней боковой кромки организован специальным образом (см. [16]), что позволяет аппроксимировать разрежение потока при О<£<0,2 уже при Х=1 на крупной сетке. Однако немонотонность решения не позволяет получить сходимости с ростом значения X. Забросы в решении методом [4] обусловлены также переменным числом узлов расчетной сетки, приходящихся на боковую поверхность при разных X. Экспериментальные данные находятся в хорошем соответствии с расчетом при £>0,5 (см. рис. 3). При £ = 0 экспериментальная точка лучше согласуется с точкой, рассчитанной методом [4].
На подветренной поверхности клина происходит взаимодействие двух сверхзвуковых потоков: одного, прошедшего через веер волн разрежения Прандтля — Майера у передней кромки клина и другого — перетекающего через верхнюю боковую кромку. В окрестности боковой кромки поток претерпевает интенсивное разрежение. Так по данным расчета методом [1, 2] вблизи кромки (|= 1 на рис. 4) статическое давление понижается до 0,2 (графики построены по значениям параметров в центрах ячеек). В диапазоне —1<|<0 уровень статического давления повышается до значения, превышающего уровеиЬ давления в области /. Это обусловлено возникновением
и=ц ;« = J°; х= ч,о, tf=fz= if
внутренних скачков уплотнения на подветренной поверхности клина. Метод [4] лучше, чем [1,2] описывает положение границы областей / и //, а также четче выделяет начало разрежения (£>0) на расширяющейся части поверхности клина. В окрестности боковой верхней кромки на подветренной поверхности клина сходимость результатов, полученных методом [4] по сеткам, отсутствует. Отличие значений р/р„ рассчитанных обоими методами на участке |<0 оценивается величинами «3 — 4%. На участке £>0 отличие достигает 30—40%. В этой области оба метода дают неверный результат вследствие отсутствия аппроксимации вихря, образующегося при обтекании острой кромки (см. подробно в п. 8 данной статьи).
На рис. 5 приведены линии равных значений статического давления и чисел М в возмущенном поле у поверхности клина. Внешняя линия соответствует положению скачка уплотнения, переходящего в границу возмущенной области. В целом результаты расчета методами [1, 2] и [4] дают удовлетворительное соответствие друг с другом. Так, у наветренной поверхности клина изолинии М = const рассчитанные методом [1,2] с уменьшением шага сетки сближаются в центральной части клина (рис. 5) с изолиниями, рассчитанными методом [4]. В окрестности нижней кромки клина сходимость результатов отсутствует. Это связано с образованием характерной петли изолинией М = 2,2, рассчитанной методом [4] на сетке MX N= 18X68. Петля образуется вследствие сноса энтропии из точек расчетного поля на поверхность клина при постановке граничного условия непротекания. У боковой поверхности клина изолинии, рассчитанные методом [4], имеют нерегулярный «пилообразный» характер, что обусловлено немонотонностью метода [12—13] и особенностями построения сетки. На подветренной поверхности клина характерный излом изолинии р/роо = 0,8 связан с расчетом границы областей II и III, которая является либо характеристикой, либо слабым скачком уплотнения. Метод [1, 2] довольно сильно «размазывает» эту границу. Метод [4] «размазывания» практически не дает, но в решении при этом имеется заброс. С увеличением числа расчетных ячеек в два раза величина заброса уменьшается также приблизительно в два раза, см. рис. 5. Таким образом, анализ результатов на рис. 2—5 показывает, что наилучших результатов в решении задачи обтекания трапециевидного клина можно достичь с использованием монотонной схемы второго порядка аппроксимации на сетке,выделяющей острые боковые кромки как особенности.
8. Основной задачей экспериментальных исследований было установление степени соответствия расчетных и экспериментальных данных в областях безотрывного течения и определение тех отличий, которые вносят в картину течения влияние вязкости газа, приводящее к отрыву пограничного слоя и образованию вихрей при перетекании потока через острую боковую кромку с одной на другую поверхность клина. Сопоставление расчетных и экспериментальных эпюр статического давления дано на рис. 2, 4, 6. В областях безотрывного обтекания на наветренной (рис. 2), боковой (рис. 4) и подветренной (рис. 6) поверхностях соответствие результатов расчета и эксперимента удовлетворительное, отличие не превышает 3%. На боковой поверхности метод [4] дает результаты, наиболее близкие к эксперименту (рис. 3, пунктирные линии).
Отрыв пограничного слоя с острой боковой кромки на нижней поверхности возникает при углах атаки а< — 5°. С уменьшением угла а область отрыва расширяется и при угле а= —14°, когда нижняя поверхность располагается под нулевым углом к набегающему потоку, отрыв на этой поверхности занимает примерно половину расширяющейся части клина (рис. 6). Минимум давления соответствует месту расположения основного вихря [17]. В районе расположения вихря экспериментальные значения давления по поверхности на 40—50% ниже рассчитанных методом [1, 2]. Как видно из рис. 6, глубина провала в профиле статического давления возрастает с ростом числа М от 1,5 до 2,5. Полученные экспериментальные данные подтверждают конический характер течения в сечениях Х=0,77; 1,45; 2,7.
м=г,5; х= 2,7; 8гіФ'',бг=п0
Хотя отрыв пограничного слоя оказывает заметное влияние на распределение статического давления, область отрыва локализована и к тому же расположена на подветренной поверхности с низким уровнем давления. Это позволяет рассчитывать интегральные характеристики клина сх, су и т. д. с точностью 1,5—2%, что вполне приемлемо для практики.
ЛИТЕРАТУРА
1. Босняков С. М., Михайлов С. В. Программа расчета трехмерных сверхзвуковых течений идеального газа. Труды IX конференции молодых ученых ФАЛТ. Долгопрудный, МФТИ, 1984. — /Депонир. рук., ВИНИТИ, 13.09.84, № 6242 — 84.
2. Босняков С. М., Карпов Е. В., М и х а й л о в С. В. Об одном подходе к решению задачи интерференции крыла и фюзеляжа. — Ученые записки ЦАГИ, 1988, т. 19, № 5.
3. Годунов С. К- Разностный метод численного расчета разрывных решений уравнений гидродинамики. — Матем. сб., 1959, т. 3, № 47.
4. Босняков С. М., Коваленко В. В., Минайлос А. Н. Расчет сверхзвукового обтекания элемента плоского воздухозаборника с выделенной головной волной. — Ученые записки ЦАГИ, 1984, т. 15, № 2.
5. Босняков С. М., Минайлос А. Н., Ремеев Н. X. Обтекание клина конечной ширины сверхзвуковым потоком газа. — Ученые записки ЦАГИ, 1977, т. 8, № 6.
6. Босняков С. М., Минайлос А. Н., Ремеев Н. X. Исследование пространственного обтекания двухступенчатых клиньев конечной ширины сверхзвуковым потоком газа. — Ученые записки ЦАГИ, 1980, т. 11, № 1.
7. Д у г а и о в В. В., И в а н о в М. Я. Сверхзвуковое обтекание боковой кромки половины клина.—Ученые записки ЦАГИ, 1977, т. 8, № 6.
8. Годунов С. К., Жданов 'А. И., Семендяев К. А. Численные методы решения одномерных неустановившихся задач газовой динамики.—Доклад на Всесоюзном съезде по теоретической и прикладной механике.— М.; 1960.
9. Алалыкин Г. Б., Годунов С. К., Киреева И. Л.,
Плинер Л. А. Решение одномерных задач газовой динамики в подвижных сетках. — М.: Наука, 1970. „
10. 3 а р у б и н А. Г. Реализация конечно-разностного метода А. Н. Краи-ко для расчета обтекания комбинации крыла и фюзеляжа сверхзвуковым потоком идеального газа. — Труды ЦАГИ, 1978, вып. 1941.
11. М а к а р о в В. Е. К выделению поверхностей разрывов при численном расчете сверхзвуковых конических течений. — Ж. вычисл. мат. и матем. физ., № 5, 1982.
12. М с С о г m a k R. W. The effect of viscocity in hyper velocity impact cratering. — AIAA P., 1969, N 69— 354.
13. W а г m i n g R. F., В e a m R. M. Upwing second-order difference schemes and applications in aerodynamics flows. — AIAA J., 1976, vol. 14, N 9.
14. Moretti G. Conformal mapping for computations of steady threedimensional supersonic flows. — Numerical Laboratory Computer Methods in Fluid Mechanics. — ASME, 1976.
15. Жмакин А. И., Фурсенко А. А. Об одном классе монотонных разностных схем сквозного счета. — Препринт ФТИ АН СССР им. А. Ф. Иоффе, 1979, № 623.
16. Минайлос А. Н. Расчет сверхзвукового обтекания крыльев с учетом сходящих с кромки тангенциальных разрывоз в рамках модели, использующей систему уравнений Эйлера. —Изв. АН СССР, МЖГ, 1978, № 1.
17. Босняков С. М., Ремеев Н. X. Исследование пространственного обтекания клина конечной ширины сверхзвуковым потоком газа при наличии углов атаки и скольжения. — Ученые записки ЦАГИ, 1981, т. 12, № 6.
Рукопись поступила 8/V 1987
