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

УЧЕНЫЕ ЗАПИСКИ ЦАГН

Том XXII І99І №6

УДК 533.6.011.5: 532.582.3

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

С. В. Михайлов,^ Н. С. Яцкевич

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

На практике достаточно часто- реализуются режимы обтекания заостренных носовых частей ЛА с отсоединенной головной ударной волной. Связано это с тем, что любой сверхзвуковой летательный аппарат проходит диапазон чисел М от нуля до значения, соответствующего максимальной скорости полета. Режимы с отошедшей ударной волной характеризуйся резким ростом коэффициента сопротивления. Возникающая в связи с этим задача расчета аэродинамических коэффициентов сравнительно точно решается при помощи методов расчета потенциальных течений (линейных н панельных методов). Тем не менее, для определения полной картины течения в возмущенной области (что имеет практическое значение, так как дает возможность определять условия на входе силовой установки) необходимо использовать численные методы, реализующие интегрирование полной системы уравнений Эйлера. ^метим, что задача расчета течения с отошедшей ударной волной в случае заостренного тела имеет принципиальное отличие от аналогичной задачи в случае затупленного тела. Это отличие связано с невозможностью построения на заостренном теле регулярной расчетной сетки. Сложностью вопроса обусловливается и сравнительно небольшое число работ, посвященных его решению. Так, в работе - [1] представлены примеры расчета комбинации конус — цилиндр на режимах с отсоединенной ударной волной. В этой работе используется разностная схема Годунова [5]. Головная ударная волна в работе [1] представляется как область больших градиентов параметров. В работах [2, 3] рассматривается задача об обтекании комбинации конус — сфера. Расчеты выполнены по модифицированной схеме Воскресенского — Русанова — Бабенко [6] первого порядка аппроксимации. Существенной деталью работ [2, 3] является использование аналитического разложения решения и при.менение специального вида разностной схемы вблизи особенностей течения (ось симметрии, острый носок тела). Головная удар-

ная волна в работах [2, 3] выделяется в виде ; разрыва. Форма волны определяется путем численного решения дифференциального уравнения.

Авторами работы [4] предложен интересный подход к расчету обтекания осесимметричных и плоских воздухозаборников на режимах с выбитой ударной волной. В их' работах используется схема Годунова в совокупности с явным выделением газодинамических разрывов. Расчет проводится на многоблочной адаптированной сетке; Из-за ограниченности ресурсов ЭВМ расчет п^ктранственного течения (осесимметричный воздухозаборник под углом атаки) проводится с привлечением полуэмпирического метода интерполирования решения по угловой переменной. Это позволяет авторам добиться хорошего качества решения при малом количестве ячеек по угловой переменной.

Целью настоящей работы является разработка методики для расчета обтекания заостренных носовых частей ЛА. Задача решается на основе численного интегрирования нестационарной системы уравнений Эйлера при помощи метода Годунова — Колгана — Родионова второго порядка аппроксимации [11, 12]. Расчет ведется с явным выделением головной ударной волны.

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

Одно из требований, предъявляемых авторами к методам расчета обтекания носовых частей — это универсальность и надежность. Исходя из этого, предлагаемая в настоящей работе методика не требует учета априорной информации о рассчитываемом течении. В частности — не постулируется ни местоположение, ни даже факт существования точки торможения. Естественно, это приводит к загрублению метода по отношению к подходам, использующим априорную информацию о виде решения (например — [2, 3]). В данной работе на примере задач об обтекании комбинаций конус — цилиндр и конус — сфера рассмотрен вопрос о точности получаемого решения. Показано, что: 1) предложенная методика позволяет без специальных усилий получать достаточно точное для практики решение даже в окрестности носка конуса; 2) получаемое решение обладает сходимостью по шагу расчетной сетки; 3) рассчитанные по предлагаемой методике поля течения с точнос-

Рис. 1

тью не хуже 1 — 3% соответствуют результатам работ [2, 3, 6]. Подробно рассмотрен вопрос модификации алгоритма выделения скачка уплотнения [13, 14]. Достигнуто практически полное совпадение положений ударной волны, определенных по модифицированному алгоритму [13, 14] и по методике [6].

Схема рассматриваемой конфигурации и ее расположение в системе координат представлены на рис. 1. Тело описывается полууглом раствора конуса (8) и длиной конического участка (хКОв) -

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

#+£+-£+-£=0 " о)

где:

р ри ру ' ■ рш

ри р + ри2 риу риш

ру ; а= риу ; 6= р + ру2 ; с = руш

рш риш руш р + рш2

е (е + р)и _ (е+р)у (е+р)ш

где р — давление; р — плотность; и, V, ш — компоненты вектора скорости; е = р (е + (и2 + V2 + ш2)/2) — полная энерг"

— 1)р) —внутренняя энергия единицы объема. Значения скорости отнесены к V* (критической скорости), плотности — к р"" (плотности набегающего потока), давления — к р ""V*. На границах расчетной области ставятся условия:

— условие непротекания на твердых границах и плоскостях симметрии;

— условие Ренкина — Гюгонио на скачках уплотнения;

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

С использованием теоремы Остроградского — Гаусса запишем систему уравнений (1) в интегральной форме:

Шо* ^х^у^г + а • + Б • + <:• ^х^у^/ = О. (2)

Задача Коши решается методом установления при помощи разностной схемы Годунова — Колгана — Родионова [5, 7 —12]. Для обеспечения консервативности схема записывается в декартовой системе координат. [7]. В расчетной области в моменты времени / и / + т строится трехмерная расчетная сетка, состоящая из шестигранных ячеек. Предполагается, что узлы сетки движутся с постоянной скоростью на интервале [/, / + т].

В центре каждой ячейки в момент ' времени / задаются значения газодинамических параметров (ГДП). В каждой ячейке имеет место линейное распределение ГДП по пространственным переменным. Предполагается, что значения ГДП в каждой точке пространства на интервале [/, / + т] изменяются линейно. Значения ГДП в момент времени / + т определяются при помощи разностного аналога уравнения (2):

6 о о о о

О У‘ = 0/У + 2 ( — Уу+т) , (3)

N=1

где 1 = (1 — 1/2, j — 1/2, к — 1/2); индексы снизу означают, что ГДП берутся в момент времени /, сверху в момент времени / + т; знак градуса — в момент времени / + т/2; о, а, В, с — векторы консервативных величин из (1) и V — объем ячейки в моменты времени, определяемые индексами; Ун — приращение объема ячейки, т. е. — объем, заметаемый боковой гранью с номером

N за время т; £",,= (5°^, 5;^, 5°) —вектор внешней нормали к положению грани с номером N в момент времени * + т/2. По абсолютной величине он равен площади грани в момент времени * + т/2.

Переход к моменту времени * + т осуществляется в два шага (предиктор — корректор) [11, 12]. Шаги отличаются способом определения параметров на боковых гранях ячейки. Для того, чтобы сделать шаг предиктор, величины на гранях (в правой части формулы (3» вычисляются при помощи значений функций и градиентов в каждой отдельно взятой ячейке. Законы сохранения консервативных величин на шаге предиктор не выполняются. При помощи этого шага вычисляются приближенные значения ГДП на момент времени * + т. Они используются для определения ГДП по разные стороны граней в формуле (3) на шаге корректор. Значения ГДП на грани для шага корректор определяются путем решения задачи о распаде произвольного разрыва [5]. Величина шага по времени т вычисляется из условия Куранта — Фридрихса — Леви [5]. При определении градиентов ГДП используется модификация принципа минимальных значений производной В. П. Кол-гана [8, 9]. Обобщение на трехмерную неравномерную сетку проведено в соответствии с рекомендациями [10]. Значения градиентов считаются постоянными на интервале [*, * + т].

2. Расчет проводится на сетке, состоящей из двух блоков. Блоком называется участок сетки, которому в пространстве индексов соответствует параллелепипед. В каждом из блоков сетка строится отдельно. На границах блоков (подобластей) обеспечивается стыковка расчетных сеток. Граница блоков, а также некоторые сечения сетки представлены на рис. 1. В случае осесимметричного течения расчетная область, представленная на рис. 1, может сокращаться до сектора, образованного плоскостью г = О и плоскостью, образующейся в результате поворота плоскости г = О на угол, определяемый необходимой точностью расчета. Для построения сетки используется алгоритм, предложенный в работе [5] для плоского криволинейного четырехугольника, а также процедура, позволяющая использовать этот алгоритм для построения сетки внутри криволинейного шестигранника. Для задания границ подобласти необходимо построить поверхностные сетки, соответствующие различным типам граничных условий (твердое тело, скачок уплотнения, поверхность стыковки подобластей, плоскость симметрии, граница сверхзвукового выхода и т. д.). Остановимся подробно лишь на алгоритме выделения в процессе установления скачка уплотнения.

При разработке предлагаемого в данной работе алгоритма выделения скачка уплотнения для трехмерной нестационарной задачи использованы идеи работ [13, 14]. Представленные ниже рассуждения применимы в случае выделения скачка уплотнения в равномерном набегающем (внешнем) потоке. Пусть фронт ударной волны в момент времени * аппроксимируется сеткой из четырехугольных ячеек. Рассмотрим один из узлов сетки (точка О на рис. 2, а). В этом узле сходятся грани, принадлежащие четырем соседним ячейкам. Известны параметры по обе стороны от скачка уплотнения. Пусть с внешней стороны от скачка расположен равномерный поток. Задача состоит в том, чтобы определить положение узла О в момент времени * + т.

Обозначим четыре сходящиеся в узле грани как 1, 2, 3, 4. Для каждой из них имеется единичный вектор внешней нормали 5і и набор ГДП с внутренней стороны от скачка уплотнения (р, р, и, и, ш). Набор ГДП с внешней стороны от скачка для всех четырех граней одинаков и соответствует набегающему потоку.

Перейдем в подвижную систему координат. Скорость этой системы равна скорости набегающего потока. Тогда первоначальные наборы ГДП заменятся на следующие:

(р, р, (и, и, ш)°),- = (р, р, (и-иоо,и-иоо, ш-ш„))(;

(Р, р, (и, и, ш)°)оо = (р, р, (О, О, 0)00.

, Рис. 2

В дальнейшем изложении вплоть до возврата в неподвижную систему координат знаки градуса у величин опускаются. ■

Определим нормальные (М) и тангенциальные (Г,) компоненты вект°р°в скорости по отношению к соответствующей грани:

01 = (и, V, - да),-; N = 0(- • 8,-; М = №• 8,-;

_ N"",= 0; 1 = 0,-М; 1"", = 0.

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

о(» Р&>1’ РОО,-

истец ется

определим следующим образом:

В выбранной подвижной системе координат , qt всегда положительно- Будем считать, Что узел О смещается по направлению средней нормали 8, которую

8°= 28,.; 8=8°/|8°|.

,=1

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

Рис. 3

диусами, равными скорости крайнего возмущения, распространяющегося в соответствующую сторону. Рассмотрим для простоты аналогичную двумерную нестационарную задачу. На рис. 2, виг изображены возмущения от двух соседних отрезков фронта волны. Пусть 3! — вектор, соединяющий центр отрезка номер I и рассматриваемый узел. Тогда требуемое смещение будет вычисляться по следующим соотношениям:

{ dS,■ = ца, если &) > О;

I ^5, = цт:/ (5(-^5), если 5;) <О.

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

Для дальнейшего решения задачи необходимо предложить алгоритм вычисления смещения узла О вдоль средней нормали. Авторы работ [13, 14] рекомендуют для этого выбирать максимально удаленную точку. В данной работе также была предпринята попытка установить ударную волну, выбирая максимальное смещение. Однако эта попытка не имела успеха. На рис. 3 в верхней части представлены положения ударной волны, характерные для процесса установления при выборе максимального смещения. Внимательное рассмотрение позволяет выяснить причину такого поведения решения. На рис. 2, б представлены два различных способа определения смещения фрагмента ударной волны АВ: смещение ЛтахВтах — по максимуму, ЛтМВ^п— по минимуму смещения вдоль средней нормали. Заштрихованный участок соответствует области, для которой при выборе смещения ЛтахВтах не может быть указан источник возмущений. По существу, в этой области нарушается условие Куранта — Фридрихса — Леви. Сказанное выше заставляет авторов данной работы использовать минимальное смещение вдоль средней нормали:

dS = MIN{dSi} •

Окончательно, смещение узла будет равно сумме смещения вдоль средней нормали с выбранной скоростью в подвижной системе координат и

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

+ ^5. 8 + (и, V, т) 00 • т.

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

Для оценки точности определения положения скачка уплотнения проведено сравнение с результатами работы [6]. Для этого выполнены расчеты обтекания комбинации конус — цилиндр на нескольких режимах, характеризующихся присоединенной ударной волной. На всех режимах отличие в положении ударной волны и в распределении параметров не превышает 0,5% от эталонного решения [6]. Этот результат говорит о том, что изложенный выше алгоритм выделения скачка уплотнения может использоваться при численном интегрировании нестационарной системы уравнений Эйлера.

3. Известно, что получение решения в ячейках с вырожденными границами не может быть гарантировано. Связано это с тем, что теоретически расчетный шаг для таких ячеек равен нулю. Для преодоления этого затруднения применяются различные меры. Так, например: в работах [2, 3] при применении дифференциальной схемы используется специальная асимптотическая запись уравнений вблизи носка тела и оси симметрии; по данным работы [10] при расчете обтекания внешних двугранных углов сверхзвуковым потоком с применением веерных сеток и схемы Годунова удовлетворительно работает процедура интерполяции или даже сноса ГДП в ячейках, имеющих вырожденные границы.

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

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

Распределения ГДП построены для расчетов, выполненных на трех различных сетках:

10 Х 10 — сетка, содержащая 10 ячеек по вертикали и по горизонтали в каждой подобластн;

20 Х 20 — сетка, содержащая 20 ячеек по вертикали и по горизонтали в каждой подобласти;

30 Х 30 — сетка, содержащая 30 ячеек по вертикали и по горизонтали в каждой подобласти.

Полуугол раствора конуса составлял 45°, число М набегающего потока — 1,5, угол атаки — 0°.

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

Рис. 4

тела, приближаются к нулю. Получаемая в расчете эпюра статического давления практически не зависит- от шага сетки.

, Критерием точности получаемых результатов могут' служить потери полного давления, обусловленные погрешностью схемы. Для решения, полученного на сетке 30 Х 30, эти потери на отрезке оси симметрии от ударной ^шны до носка тела сбавляют менее 1%, на сетке 10 Х 10 —не более 2%. Реализующееся в решении значение полного давления на оси симметрии для всех сеток с точностью в 0,5% соответствует потерям в прямом скачке.

4. На практике основное требование к методике расчета обтекания носовых частей ЛА — это возможность расчета аэродинамических коэффициентов носовой части. В случае, если с некоторого сечения х = const ведется маршевый счет, методика. должна обеспечивать достаточную точНОсть параметров в этом сечении. Для оценки точности получаемого решения произведено сравнение с результатами работы [3]. Выполнен расчет обтекания комбинации конус — сфера. Полуугол раствора конуса — 45°. Число Маха набегающего потока — 2. На рис. 5 представлено распределение давления и х-компо-ненты скорости по поверхности тела н в выходном сечении х = const. При построении графиков использовано обезразмеривание и координаты из работы [3]. Видно, что для давления отличие нигде не превышает 1 %. Для плотности сравнение дает аналогичные результаты. Точность расчета х-ком-поненты скорости в выходном сечении — не ниже 3%. Для значений х-компо-ненты скорости вблизи заострения (не более 5—10% длины тела) имеет место значительное расхождение. Однако на остальной части тела точность опред^фния х-компоненты скорости не ниже 3%. Для оценки точности определения коэффициента сопротивления произведено сравнение с результатами работы [6] на примере расчета обтекания конуса с полууглом раствора 45° на режимах с присоединенной ударной волной. Расчеты произведены на сетке 10 Х 10 (см. п. 3). В диапазоне 2,2 <М <4,1 точность определения коэффициента с* не ниже 1—2%. При М>4,1 возможно применение маршевых методов.

5. В качестве примера практического использования предлагаемой методики на рис. 6 представлены диаграммы коэффициента сопротивления и относительного расстояния отхода ударной волны. Данные получены при расчете на нулевом угле атаки комбинации конус — цилиндр в случае отошедшей ударной волны. В расчетах варьировались полуугол раствора конуса (0 = 35°, 40°, 45°) и число Маха набегающего потока (1,01 <М<МКр(0)). Использовалась сетка 1 О Х 10. .

Авторы благодарят С. М. Боснякова за ряд полезных замечаний, высказанных в процессе подготовки-работы.

ЛИТЕРАТУРА

1. И в а н о в М. Я. К решению двумерных и пространственных задач обтекания тел околозвуковым потоком.— Ж. вычисл. матем. н матем. фнз.,

1975, т. 15, М 5.

2. И в а н о в а В. Н., Ра д в о г н н Ю. Б. Численный метод расчета трехмерных обтеканий головной части заостренных тел с отошедшей ударной волной.— Препринт ИПМ, 1980, № 126.

3. И в а н о в а В. Н., Ра д в о г и нЮ.Б. Численное нсследование трехмерного обтекания заострен'!ых тел с отошедшей ударной волной.— Препринт ИПМ, 1981, М 28.

4. М и л е ш н н В. И., Т нлл я е в а Н. И. Сравненне расчетных и экспериментальных данных по обтеканию осесимметричных воздухозаборников на режимах с выбнтой ударной волной.— Ученые записки ЦАГИ, 1981, т. 13,

М 2.

5. Г о д у н о в С. К., 3 а б р о д н н А. В., И в а н о в М. Я., К р а й-к о к Н., Прок о п о в Г. П. Численное решенне многомерных задач газовой динамики.— М.: Наука, 1976.

6. Б а б е н к о .К. И., В о с к р е с е н с к и й Г. И., Л ю б н м о в А. Н.,

Ру с а н о в В. В. Пространственное обтекание гладких тел идеальным газом.— М.: Наука, 1^И.

7. 3 а р у б и н А. Г. О точности метода Годунова в различных сис-темах.координат.— Ученые записки ЦАГИ, 1977, т. 8, № 4.

8. К о л г а н В. П. Применение принципа мннимальных значений произ-

водной к построению конечно- разностной схемы для расчета разрывных решений газовой динамики.— Ученые запискн ЦАГИ, 1972, т. 3, М 6. •

9. К о л г а н В. П. Конечно-разностная схема для расчета разрывных течений нестационарной газовой дннамнки.—— Ученые запискн ЦАГИ, 1975, т. 6, № 1.

10. Т и л л я е в а Н. И. Обобщение модифицнрованной схемы С. К. Годунова на произвольные нерегулярные сетки.— Ученые запискн ЦАГИ, 1986. т. 17, №> 2.

11. Р о д н о н о в А. В. Монотонная схема второго порядка аппроксимации для сквозного расчета неравновесных течений.— Ж. вычнсл. матем. и матем. физ., 1987, т. 27, № 4.

12. Р о д н о н о в А. В. Повышение порядка аппроксимацни схемы

С. К. Годунова.— Ж. вычисл. матем. и матем. физ., 1987-,

13. К райк о А. Н., М а к а р о в В. Е., Т и л л я е в а Н. И. К чнс-ленному построению ударных волн.— Ж. вычисл. матем. и матем. фнз., 1980. т. 20, №> 3.

14. М а к а р о в В. Е. К выделению поверхностей разрыва прн численном расчете сверхзвуковых коннческих течений.— Ж. вычисл. матем. и

матем. физ., 1982, т. 22, М 5.

Рукопись поступила /8^/1 /990 г.