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

УЧЕНЫЕ ЗАПИСКИ Ц А Г И Т о м VII 197 6

№ 3

УДК 533.6.011.55

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

А. Н. Минайлос

Методом сквозного счета рассчитано обтекание пластины прямоугольной формы в плане в сверхзвуковом потоке под углами атаки до 30°. Решение построено без учета отрыва потока от боковых кромок пластины (условие Н. Е. Жуковского не использовалось). Проанализированы полученные результаты расчетов и схемы течения. Приведены аэродинамические характеристики.

Автором совместно с И. В. Иерусалимской и А. П. Косых разработана программа для расчета сверхзвуковых пространственных течений газа. Двумерный вариант примененного метода расчета описан в [1]. Метод и программа будут подробно изложены в отдельной работе. В этой статье рассмотрены некоторые результаты исследования нелинейных конических течений предложенным методом.

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

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

Система координат хуг связана с телом: ось х направлена вдоль его оси. Система газодинамических уравнений должна быть такой, чтобы продольный компонент скорости и (вдоль оси х) был больше скорости звука.

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

7—Ученые записки

97

проводится на каждом шаге hx по оси х. Размеры этого шага диктуются критерием устойчивости применяемой конечно-разностной схемы.

В случае расчета течения около прямоугольной пластины она помещается между двумя слоями ячеек в плоскости у =1 и аппроксимация поверхности не нужна (ее результаты просто дают точную форму поверхности). В направлениях у и z сетка имеет шаги Лу = Лг = 0,05 (все линейные размеры отнесены к половине длины пластины в направлении оси г).

Используемая конечно-разностная схема является стационарным аналогом схемы Р. П. Колгана [2]. Схема [2] имеет первый порядок аппроксимации по временной переменной и .почти второй”—по пространственной. В нашем случае это означает, что хотя в целом схема имеет первый порядок, в поперечных плоскостях х = const поля течения описываются точнее [I], чем при использовании метода С. К. Годунова [3, 4].

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

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

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

разными знаками больших ошибок в вычислении давления и плотности при описании течения типа течения Прандтля—Майера в окрестности точки разворота потока.

При использовании схемы [1], в которой величины разрывов на границах ячейки меньше, чем в схеме [3], влияние первой причины ослабевает.

Вторая причина обусловлена потерей аппроксимации схемы вблизи точки разворота потока. Для уменьшения фиктивного энтропийного слоя, увеличивающего энтропийную функцию в областях расширения на 50—70%, используется модификация метода, автоматически включающая в расчет в областях расширения условия изоэнтропийности вместо проекции уравнения движения на ось х. Цель такой замены состоит в том, чтобы в области, где схема теряет аппроксимацию, ошибки в значениях давления и плотности, сохраняя свою величину (уменьшить их по модулю практически не удается), имели бы одинаковый знак, что существенно уменьшит ошибку в значении энтропийной функции.

Ошибку в норме С энтропийной функции удается уменьшить с 71 до 2,3%, а в числе М с 19 до 2,4%.

Ошибки при вычислении интеграла Бернулли и при решении уравнения расхода в поле течения очень малы.

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

Расчеты и обработка результатов. Расчеты проведены в диапазоне изменения чисел 2 <; <; 8 и углов атаки 0<а^30° на сетке 65 X 50 узлов на

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

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

скорости Vmax, давление —к р^ 1^ах.

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

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

Положение дискретного разрыва можно определить по профилям распределения функций, например, плотности, давления, энтропийной функции, числа М. Было принято считать местом положения дискретного разрыва точку со значением функции, соответствующим середине высоты размытого разрыва, в случае определения положения скачка уплотнения по значениям конического числа М — точку Мкон = 1*. Последний вариант позволяет просто определять положение ударных волн (в программе обработки результатов есть блоки расчета изолиний в плоскости уг) и поэтому использован наиболее широко. Отметим, однако, что условие Мкон = 1 не является ни необходимым, ни достаточным условием существования скачка в данной точке пространства.

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

Дтах = И13Х {Д,Утах> Д^тах} = 5%.

Схема течения. Результаты расчетов подтверждают предположение [5] о существовании над поверхностью прямоугольного крыла „висячих* конических ударных волн.

Рассмотрим поле изохор, изображенных в виде тонких сплошных линий на •фиг. 1 (х — 4, = 5, се=15°). У плоскости симметрии плоское течение сжатия —

снизу крыла и течение расширения — сверху. Треугольниками отмечено точное положение ударной волны снизу (течение у клина) и первой возмущенной характеристики течения Прандтля—Майера сверху. Поскольку число М около плоскости симметрии под крылом меньше, чем сверху, возмущения от угловой точки под крылом распространяются в направлении отрицательной оси г дальше и при х = 4 уже достигают плоскости симметрии. Это область разрежения, уменьшающая несущие свойства крыла. Сплошные жирные линии изображают -сечения ударных волн, построенных в результате анализа полей скоростей, изобар, изохор, изоэнтроп и изомах конических. Положение этих волн определено по изомахам и изохорам. Поверхности Мкон = 1 показаны белыми кружками.

Основная ударная волна I далеко отходит от кромки крыла в направлении -г и не замкнута над верхней поверхностью. Интенсивность волны максимальна в районе плоскости симметрии течения. Измеренная перепадом энтропийной функции 5 =/?//?*, она равна = 16% (отнесена к 5^=0,0238). Две .висячие' ударные волны 2 и 3 располагаются над крылом в потоках, прошедших соответственно над боковой и передней кромками. Эти потоки разделяются поверхностью контактного разрыва 4. Положение ее получено как положение изоэн-тропы 5 = .9^. Это стало возможным, несмотря на наличие волны 3, из-за систематической ошибки в значении энтропийной функции при расчете течения Прандтля — Майера над крылом. Эта ошибка, характерная для использованного метода, занижает величину 5 на 2,3% по сравнению с а в волне 3 величина 5 увеличивается на 1%. Таким образом, за ударной волной 3 значение 5“ оказалось меньше, чем 5^. Газ же, прошедший над боковой кромкой, имеет значительно большую энтропию. Таким образом, изоэнтропа 5 = 5^ лежит в пределах размытого контактного разрыва. В целом картина течения для системы ударных волн 2 и 3 и контактного разрыва 4 аналогична картине течения, получающейся при столкновении двух плоских сверхзвуковых потоков с образованием двух ударных волн [4]. Такой элемент течения был получен также в расчетах обтекания треугольных, трапециевидных крыльев и крыльев с изломом передней кромки. Назовем его элементом СКС (скачок — контактный разрыв — скачок).

На фиг. 1 штриховыми линиями показаны следы изоэнтропических поверхностей. Огибая боковую кромку, поток разгоняется, проходит звуковую поверхность, достигает, расширяясь, числа МКОн = 2,5 и тормозится в ударной волне 2 {см. фиг. 1). Эта волна содержит участок конически прямой волны и по интенсивности превосходит волну 1 в 2,6 раза. Полученный аналитически прирост энтропии в прямой волне с Мдф = 2,5 составляет Д52 = 42%. Таким образом, энтропийная функция за волной 2 должна достигать значения

^ = $»(!+ Д5, + АЯг) = 0,040/

* Здесь И ниже Мкон =-[(Ув+ У*ЖХР/Я)11*2> ^8 и У — компоненты полного вектора скорости в конической системе координат с центром в угловой точке крыла.

Однако в расчетах энтропия оказалась значительно большей. В этой коническк дозвуковой части потока при численном решении существует вихрь, центр которого отмечен на фиг. 1 белым ромбом. Для зоны вихря характерны пониженные значения давления, плотности, конического числа М и повышенные значения энтропийной функции (5 достигает величины 0,0506). Обозначим прирост энтропии по сравнению с значением через Д5д. По-видимому, он вызван ошибками аппроксимации конечно-разностной схемы и может характеризовать интенсивность вихря. В рассматриваемом варианте Д5д = 44,5%. Область, в которой 5>52, на фиг. 1 заштрихована. Приращение энтропии возрастает при увеличении угла атаки а и числа (например, при = 5 и а= 10° Д5д=8,4%).

Сделаем два замечания по поводу возрастания энтропийной функции.

Во-первых,’ объясняя увеличение энтропии влиянием ошибок аппроксимации, мы получаем неожиданно большие значения Д5д (ведь в направлениях у и г диссипативные члены присутствуют в разностных уравнениях, начиная с члена с четвертой производной) Отметим, что разрешающей способности расчетной сетки уже недостаточно, чтобы описать слой течения между волной 2 и вихрем. Для детализации и уточнения картины течения в области вихря необходима более густая сетка.

Во-вторых, необходимо помнить, что расчеты проведены для невязкого газа, и в краевых условиях на кромке крыла не выполнено условие Н. Е. Жуковского, т. е. эти условия для системы уравнений Эйлера не удовлетворяют предельному состоянию течения, полученному из уравнений Навье — Стокса при числе Рейнольдса, стремящемся к бесконечности; и в области над крылом в реальном потоке вязкость может существенно изменить картину течения.

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

Таким образом, полученная схема течения состоит из одной конически эллиптической зоны и двух конически гиперболических и содержит два элемента над поверхностью крыла, один из которых (СКС) образован столкновением потоков и содержит разрывы 2, 3 и 4, а другой является вихрем.

Результаты расчетов. Зависимость картины течения (ударные волны и положение вихря) от угла атаки при числе = 5 иллюстрируется фиг. 2. При угле атаки а = 20° картина течения получена при х = Ъ (при х = 4 коничность течения уже нарушена), а затем перестроена с учетом коничности потока. Ось вихря при увеличении угла атаки движется по кривой. При этом интенсивности вихря и скачков возрастают. При а^>15° область эллиптичности делится на две отдельные области.

Фиг. 3

Положение ударной волны в плоскости симметрии у нижней поверхности хорошо соответствует точному значению угла наклона ударной волны около клина (крестики на фиг. 2) и отражает немонотонность изменения этого угла с изменением угла атаки.

На фиг. 3 приведены кривые распределения давления по размаху крыла на нижней и верхней поверхностях (М^ = 3, х = 1) при различных углах атаки. Стрелками показано точное положение первых характеристик течения разрежения, рассчитанное по числу М плоского потока снизу крыла. Штриховая линия указывает значение давления в невозмущенном потоке. Увеличение давления на верхней поверхности при г = 0,7 0,9 соответствует конически дозву-

ковой зоне между волнами 2 и 3.

Силовые характеристики приведены на фиг. 4. Слева на фигуре показана зависимость коэффициента сп от угла атаки при различных значениях числа Мда. Для сравнения с данными эксперимента [6] пластина была выбрана длиной х = 1,87. Экспериментальные значения получены при Ма = 6,9 и показаны на фигуре светлыми точками. В эксперименте исследовалась не пластина, а клин с углом верхней поверхности 6,63°. Поэтому при а = 0 экспериментальное значение сп отрицательно. На графике звездочкой показан угол, начиная с которого верхняя поверхность находится в .ньютоновской* аэродинамической тени. Отличие экспериментальных данных от расчетных при этом значении а составляет около 30%. Только при углах а>10°-4-12° влияние верхней поверхности, ослабевает и отличие уменьшается.

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

ЛИТЕРАТУРА

1. Косых А. П., Минайлос А. Н. Исследование методов сквозного счета для задач сверхзвуковой аэродинамики. .Ученые записки ЦАГИ“, т. 7, № 1, 1976.

2. К о л г ан В. П. Применение принципа минимальных значений производной к построению конечно-разностных схем для расчета разрывных решений газовой динамики. „Ученые записки ЦАГИ“, т. 3, № 6, 1972.

3. Годунов С. К. Разностный метод численного расчета разрывных решений уравнений гидродинамики. „Матем. сборник', т. 3, № 47, 1959.

4. Иванов М. Я., Край ко А. Н., Михайлов Н. В. Метод сквозного счета для двумерных и пространственных сверхзвуковых течений. „Журн. вычислит, матем. и матем. физ.“, т. 12, № 2, 1972.

5. Булах Б. М. Нелинейные конические течения газа. М., .Наука", 1970.

6. Penland I. A. Maximum lift-drag-ratlo characteristics of rectangular and delta wings at Mach 6,9. NASA TN D-2925, 1965.

Рукопись поступила 8jX 1974 г.