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

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

Т о м VI

197 5

№ 2

УДК 533.6.013

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

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

Имеется весьма ограниченное число работ, посвященных численному решению пространственных задач газовой динамики с помощью метода установления [1]. Для решения подобных четырехмерных задач кроме большого быстродействия и памяти ЭЦВМ требуется эффективный метод сквозного счета повышенного порядка точности. Целью настоящей работы является применение нового конечноразностого метода, основанного на принципе минимальных значений производной [2], к задаче пространственного обтекания тела под углом атаки.

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

Законы сохранения массы, импульса и энергии могут быть представлены в дивергентном виде:

В. П. Колган

Р(х,У, г) = 0.

(О

■Ж + ~Тх И) + 1} (Р®) + 4^ №) = 0;

"зг (ри> + &+р“2)++ж (Ри=°;

Через р, и, V, ге>, р обозначены соответственно плотность, компоненты вектора скорости V и давление, отнесенные к характерным размерным величинам р«,, | ^оо! и р<» У», взятым на бесконечности.

Рассмотрим граничные условия. На теле должно выполняться условие непротекания

где« — внешняя нормаль к поверхности тела.

На бесконечности зададим однородный невозмущенный поток

тела, Моо — число М на бесконечности.

Пусть в момент времени t = 0 имеется некоторое начальное поле течения около тела. Тогда, применяя метод установления, стационарное решение будем искать как предел нестационарного решения при 09.

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

Начало координат поместим на поверхности тела (фиг. 1). Луч х^-0 направим так, чтобы он целиком лежал внутри тела. Будем рассматривать тела, симметричные относительно плоскости г = 0.

(3)

Рос = 1; «со = соэ а; -у«, = эШ а; те»«, = 0; рх = —5-, где а — угол атаки

. хМ1, .

г

а)

5)

Фиг. 1

Тогда достаточно рассмотреть течение в полупространстве 0. В случае тела вращения ось х совпадает с осью тела. Вокруг оси х вводится угловая координата <р так, что полуплоскость 2=0, у<Ц 0 соответствует углу ср = 0, а полуплоскость z = 0, .у>0 соответствует углу cp = ir. Пространство между этими полуплоскостями разбивается на J частей координатными плоскостями <р = сру (у = 0, 1, 2,.. . /) с требуемым шагом, равномерным или неравномерным. Линия пересечения поверхности тела с плоскостью <р = 0 приближенно аппроксимируется ломаной, состоящей из / отрезков, каждый из которых имеет длину lt и образует с осью х угол Ь. (i = 1, 2, . .., /). Далее в плоскости <р = 0 из каждого угла полученной ломаной линии

проводятся лучи, образующие с осью х углы {5,(/==0, 1, 2......../).

Осевому лучу при угле р0 = 0 соответствует отрицательная часть оси х. На этом луче откладывается К отрезков, длина которых hk увеличивается по мере удаления от тела по закону геометрической прогрессии со знаменателем ^>1

hk — haqk (£ = 0, 1, 2,..., AT—1).

В плоскости <р = 0 через точки разбиения нулевого луча проводятся (/С-М) параллельных линий под углом Qj к оси х до пересечения с первым лучом, затем под углом 62 до пересечения со вторым лучом и так далее до последнего луча с номером /. В результате получается расчетная сетка для плоскости 9 = 0.

Через узлы этой сетки от плоскости ср = 0 проводится ряд параллельных линий до пересечения с плоскостью = «р,. Эти линии параллельны плоскости х = 0 и составляют угол а! с осью у. Точки пересечения этих линий с плоскостью <? = ?! являются узлами расчетной сетки для плоскости f = tp1. Далее под углом а2 аналогичным ОбраЗОМ ПрОВОДЯТСЯ ЛИНИИ ОТ ПЛОСКОСТИ ср = cpj до плоскости 9 = <р2 и так далее до последней плоскости ср = тс. Осевой луч Ро = 0 является общим для всех плоскостей ср = Значения ctj (j= 1, 2,..., J) определяют форму поперечного сечения тела.

Итак, построение пространственной сетки закончено. Если требуется .построить сетку, состоящую из JXJXK ячеек, то необходимо задать пять массивов чисел:

<ру (/' = 0, 1, 2,. .. У); /„ 0г (¿= 1, 2,.../);

^ (¿ = 0, 1,2..../); а у (у = 1, 2,... J).

Дополнительные два параметра А0 и q определяют размер сетки в направлении от тела.

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

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

схему сквозного счета, в которой местоположение разрывов определяется автоматически. Одним из таких методов является известный конечноразностный метод, предложенный Годуновым С. К. [3]. Однако точность этого метода недостаточно высока, так как для всех функций применяется кусочно-постоянная аппроксимация по ячейке. В работе [2] предложена новая конечноразностная схема, основанная на применении принципа минимальных значений производной, что позволяет повысить порядок аппроксимации по пространственным переменным с первого порядка до второго.

Этот метод и был положен в основу при составлении программы расчета на БЭСМ-6 пространственного течения газа с использованием описанной выше расчетной сетки.

Дадим описание граничных условий, которые использовались при расчетах на границах расчетной сетки. С внешней стороны граничного слоя ячеек /г = АТвыполнялось условие невозмущенного потока. На последнем слое ячеек i — I, ограничивающем расчетную область вниз по потоку, интерполяция параметров внутри ячеек производилась с использованием значений параметров соседних внутренних ячеек. В ячейках, граничащих с поверхностями непро-текания (тело, плоскости <¡> = 0 и <р = я) интерполяция проводилась также с использованием внутренних ячеек, а условие непротекания обеспечивалось при расчете распада соответствующих разрывов*

Ускорение сходимости решения и экономия памяти ЭЦВМ. Обычно в методах установления величина шага по времени Д£ определяется из соображений устойчивости счета и равна минимуму предельно допустимых шагов по всем ячейкам расчетной сетки. Величина At оказывается пропорциональной размеру самой маленькой ячейки сетки. Поэтому существенная неравномерность координатных шагов вызывает большое увеличение времени расчета, которое используется нерационально: для крупных ячеек шаг повремени слишком мал и параметры в них меняются слишком медленно, т. е. участки с крупными ячейками желательно считать с соответствующим крупным шагом, что невозможно из-за появления неустойчивости решения на участках с мелкими ячейками.

Применяемая в настоящей работе расчетная сетка может обладать значительной неравномерностью. Например, при <7 = 1,2 и ЛГ=20 отношение поперечных размеров большой и маленькой ячеек достигает 32. Чтобы избежать нерационального увеличения времени счета при установлении в настоящей работе расчет ведется с преобразованием времени в зависимости от координатного шага сетки: каждая ячейка рассчитывается со своим шагом по времени, близким к предельно допустимому для этой ячейки, что позволило существенно (примерно на порядок) сократить время решения задачи на ЭЦВМ. При указанном „деформировании“ времени теряется соответствие между физическим и фиктивным временем, фигурирующим в расчете. Но это является несущественным, поскольку ищется предельное решение при Ь ОС.

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

При расчетах давление всюду, кроме вычисления распада разрывов, определялось из уравнения Бернулли. Это позволило сократить на 20% потребную память ЭЦВМ и исключить из рассмотрения уравнение энергии в исходной системе уравнений (2). Численный эксперимент показал, что введение интеграла Бернулли не повлияло на устойчивость счета.

Вычислительная машина БЭСМ-6 производит все расчеты с точностью 12 значащих цифр. Однако точность газодинамических расчетов обычно находится в пределах трех-четырех значащих цифр, не более. Поэтому разумно поставить вопрос о хранении промежуточных результатов расчетов с пониженной точностью, что позволит без потери точности расчетов в целом сэкономить оперативную память машины и не пользоваться внешними запоминающими устройствами. В настоящей работе было применено „упаковывание“ числовой информации, т. е. хранение большого количества чисел, но с меньшим количеством знаков. В процессе счета определялся диапазон изменения каждой функции, который разбивался на 4096 интервалов. В памяти машины хранилось не значение функции, а номер соответствующего интервала, в который попадало значение функции до упаковки. При необходимости произведения с упакованными числами некоторых действий машина производила их „распаковку“, т. е. делала приведение чисел (приближенно) к первоначальному десятичному виду. Указанный прием позволил хранить в одной ячейке памяти по четыре числа, соответствующие параметрам газа в одной ячейке расчетной сетки. Расчеты показали, что при этом не происходит видимой потери точности получаемых результатов.

Примеры расчетов. Приведем для примера результаты расчета обтекания осесимметричного тела — цилиндра со сферическим затуплением. Расчетная сетка имела размеры /Х-^Х Л=30 X 17X20. Общее число неизвестных параметров (исключая давление) равнялось 40800. Начальное поле течения задавалось однородным с параметрами, равными параметрам на бесконечности. Для получения установившегося режима требовалось сделать примерно 800 шагов по времени, что занимало 12 часов машинного времени БЭСМ-6.

На фиг. 2 приведено распределение коэффициента давления

ср = 2 (р—роо)1роэ V2 по поверхности тела, обтекаемого сверхзвуковым потоком с числом Моо= 1,5. На сферическую часть тела приходится 18 расчетных точек с шагом Д0 = 5°. Результаты расчета при угле атаки а = 0 практически везде совпадают (с точностью 0,2%) с данными обтекания сферы, взятыми из работы [4], в которой

течение около лобовой части сферы рассчитано по методу Годунова С. К. с координатным шагом Д0=1,44°, однако вблизи оси тела 6 = 0 расхождение результатов несколько больше и достигает 4—5%. Этот эффект, даваемый разностной схемой, объясняется особенностью выбранной системы координат — ячейки, прилежащие к оси 9 = 0 имеют только пять граней, а не шесть, как все остальные. Крестиками нанесены результаты расчета по настоящему методу, но с применением в ячейках кусочно-постоянной аппроксимации (аналогично методу Годунова С. К.). Погрешность расчета в этом случае существенно больше (—3%).

На фиг. 3 приведены результаты расчета обтекания тела потоком газа, имеющим звуковую скорость на бесконечности М=1. Полученные данные хорошо согласуются с результатами эксперимента^], а при а = 0 и с численными расчетами работы [4].

ЛИТЕРАТУРА

1. Богачевский, Мейтс. Прямой метод расчета течения около затупленного осесимметричного тела при угле атаки. РТК, 1966, № 5.

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

3. Годунов С. К., Забродин А. В., Прокопов Г. П. Разностная схема для двумерных нестационарных задач газовой динамики и расчет обтекания с отошедшей ударной волной. „Журн. вычислит, матем. и матем. физ.“, т. 1, № 6, 1961.

4. К о л г а н В. П., Фонарев'А. С. Установление стационарного обтекания при набегании ударной волны на сферу и цилиндр. „Изв. АН СССР, МЖГ“, 1972, № 5.

5. Л е у т и н П. Г. Распределение давления по лобовой поверхности тел вращения при а = 0 — 10° в трансзвуковом потоке. .Ученые записки ЦАГИ“, т. V, № 2, 1974.

Рукопись поступила 29/ІІІ 1974 г.