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

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

Том XXII 1991 № 6

УДК 533.6.011.5

МЕТОД ДВУХ ПОВЕРХНОСТЕЙ ПОСТРОЕНИЯ

расчетной сетки в задачах

СВЕРХЗВУКОВОГООБТЕКАНИЯ

Н. Н. Бабаева, И. В. Орлов ,

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

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

Одним из способов. которые позволяют осуществить не только простую адаптацию расчетной сетки к геометрии течения, но и выполнить ряд высоких требований к сетке в смысле ее'ортогональности вблизи границ области.- управления расстановкой узлов, а также уменьшения затрат машинного времени, является метод трансфинитной интерполяции. Он описан в работах [1—3].

Частным случаем метода трансфинитной интерполяции является метод двух- поверхностей. В этом случае метод трансфинитной интерполяции выглядит следующим образом.

Заданы две непересекающиеся поверхности Р,(5, и) и ^2(5, и) как векторные функции:

Р, = [Х| (5, и), у, (5, и), 2, (5, и)], Р2 = [Х2(5, и), 02(5, и), 22(5, и)]

двух параметров 5, и, изменяющихся в интервалах 5иач ^ 5 ^ 5К0Н, О ^ и ^ 1. Параметр 5 в методе двух поверхностей играет формальную роль и может быть определен пользоватМем произвольно. Требуется построить функцию

Х(5. и. /) х(Р,(5, и),Р2(5. н),/)

II 3 со С У(5. и. /) = у(Р,(5, и),Р2(5, н),/)

2(5, и. /) 2(Р,(5, и),Р2(5, Н),/)

осуществляющую отображение вычислительного пространства (s, и, t) на некую физическую область . с 'двумя заданными границами Pi и Р2 Остальные границы области получаются в результате самого отображения. Предполагается, что переменная t также нормирована и изменяется в единичном интервале: О ^ t ^ 1, причем

F(s, и, О) = Pi (s, и), F(s, и, 1) = P2(s, и).

Программный модуль строит сетку внутри четырехугольной области, две стороны которой прямолинейны, две другие заданы дискретным набором точек или аналитически. Точки сетки между двумя криволинейными границами соответствуют заданным значениям параметра / ’из интервала (О; 1). В качестве значений параметра s взяты порядковые номера точек вдоль криволинейных границ.

В аналитическом виде границы задаются следующим образом:

Zi = Zi(s), Yi = Yi(s), Z2=Z2(s). Y2 = Y2(s), l<S <S.

Преобразование имеет следующий вид: -

Z(s, /) = ./1|(t)Zi(s) + /12(t)Z2(s) + /13(t)TiD|z(s) + /14(t)'t2D2z(s), y(s, t) = /11 (t)Yi (s) + /12(t)Y2(s) + H3(t)'tiD|y(s) + /14(t)'t2D2y(s).

В этих формулах /1<(t) — функции смешения, D| и П2 — векторы, по направлению которых координатные линии S = const покидаюг поверхности 1 и 2 соответственно.

В ■ общем случае D;=0M + (1 — В) Т, где Nj и T — векторы нормали и касательной к границе, 0 — свободный параметр, задаваемый пользователем. При' 0 = О линии S = const подходят по касательной к поверхностям 1 и 2, а при 0=1 — по нормали. 'ti, 'Т2 — коэффициенты, контролирующие длину векторов нормали на границах области.

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

Частными случаями метода двух поверхностей являются линейная, квадратичная и кубическая интерполяция.

Для линейной интерполяции используется набор функций смешения:

/11=1 -- t, /12=/, /13= /14 = о.

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

Для квадратичной интерполяции выбирается один из двух наборов функций смешения:

/11=1— Л /12= Л /13 =/(1—0, /14=0

или

/1i=(1—t)2, /12 = 2/ — /2, /13 = 0, /14= — /(l—t).

Для кубической интерполяции берется следующий набор:

/11 = 2t3 — 3t2 + 1, /12= — 2t3 + 3t2,

/13=t (l —t)2, /14 ==t2(t—1).

При ненулевых ТА: и 0= 1 сетка, получающаяся в результате такого преобразования, ортогональна к поверхностям Pk(s, и) . Увеличение ТА: увеличивает ортогональность в ячейках, расположенных в глубине расчетной области.

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

В этом случае границы заданы двумя таблицами:

P.(S) = [xi(s),yi(s)]J,

P2(s) — [X2(S), Y2(S)]J, (i=-l,O, . . . ,/,/ + 1).

Необходимые для построения сетки производные можно вычислить с помощью центральных разностей.

(4т) U/W. «'—0 1................1

Аналогично вычисляются > ’ (“IT") '

Преобразования записываются теперь в виде:

x(s, t) — J.I.,(t)x,(s) + J.l.2(t)*2(s) + J.l.3(t)T, (— -M) + J.l.4(t)T2 (—"Ter) ’ y(s, t) — (t) у i (s) + J.l.2(t)Y2(s) + 1 ('11')+

Одним из преимуществ дискретного задания границ является свобода в размещении точек на границах при выбранных значениях параметра s.

В методе заложена также возможность контроля за расположением узлов сетки на линиях S — const. Этот контроль заключается в определении таких значений параметра t, при которых расположение вычисленных узлов сетки удовлетворяло бы поставленным условиям. В качестве таких условий можно взять требование равномерного расположения точек или расположения с заданным заранее отношением длин дуг между точками.

В качестве примеров построено несколько расчетных сеток. На рис. 1, а построена сетка в поперечном сечении области между треугольным крылом с затупленными передними кромками и ударной волной.

Для этого примера взяты значения т\ = Т2 = 3. Сетка ' имеет сгущение точек у поверхности тела.

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

Опыт показывает, что подбор коэффициентов Т|, Т2 не составляет особой трудности. В то же время необходимо обратить особое внимание на это, поскольку неудачный подбор коэффициентов Т|, Т2 может привести к самопересечениям сетки.

-2.0 -1,0 о

1.0 2,0

Рис. 2

".0 5;0

2. Достоинства изложенной методики получения ' расчетных сеток демонстрируются на примере задачи о течении идеального совершенного газа С. показателем адиабаты Пуассона х = 1,4 в плоском канале переменного сечения. Профилирование канала ' заключается в том, что часть прямолинейного участка нижней стенки (рис. 2) заменяется препятствием, Вид которого задается одним периодом гармонической функции. ' Отношение высоты препятствия к ширине канала составляет 1/4.

Если предположить, что сетка будет иметь постоянную размерность 30: 100 и неизменное положение точек на границах, то естественно ожидать изменения сеточных линий только в ограниченной области большой кривизны профиля. Поэтому на рис. 2 приведен характерный фрагмент двух расчетных сеток: А — примитивная интерполяция для линейных функций смешения; В — усложненная интерполяция для кубических функций. Можно видеть, что сетка В отвечает повышенным требованиям к ортогональности сторон расчетных ячеек. При этом ожидается, что данные изменения сеткй приводят к улучшению всего численного решения основной задачи о течении газа.

Решения на полученных сетках находились с использованием конечноразностной схемы С. К. Годунова (СГ) [4]. В качестве же критерия оценки полученных результатов использовались аналогичные данные для модификации СГ для кусочно-линейной аппроксимации распределения по пространству газодинамических величин [5] (МСГ) , имеющей повышенный порядок точности. Применение методов распадного типа к полным нестационарным уравнениям Эйлера равносильно решению системы интегральных законов сохранения для каждой ячейки с конечной площадью 5 и контуром Г, который обходится противчасовой стрелки:

\pdxdy +ф рУп= 0;

(1)

^ригіхгіу + ф (ри^„ + р!!.у) = О; (2)

5 I’

^ругіхгіу + ф (руУ„ — Р!!.х) = О; (3)

5 ('

і р( е + + ф р( * + = 0, (4)

я Г

где V,, = и!!.у — и!!.х.

Здесь и ниже р, р, е, *, V — давление. плотность, удельные энергия и энтальпия газа, модуль вектора скорости, а и, V — его проекции на оси х, у.

Систему (1) — (4) необходимо дополнить уравнением состояния газа, согласно которому выполняются равенства:

р н хР

е (ж — 1)р ’ (*— 1)р ‘

Из предварительных экспериментов было выяснено, что, начиная с определенного числа М> 4, не реализовалось запирание канала в критическом сечении, а отход скачка уплотнения не превышал размеры расчетной области. Вследствие этого скорость невозмущенного потока соответствовала

М"" = 5. Было проведено три численных эксперимента: 1, 2, применяя СГ на двух сетках, и 3 — для МСГ на сетке В.

На рис. 3, а, 6 приведены линии тока и изобары для эксперимента 3 и на рис. 3, в, г изобары соответственно для 2 и 1 эксперимента. Поведение кривых хорошо иллюстрирует основные физические особенности течения. Так, препятствие формирует ■ в набегающем газе отошедший скачок уплотнения, отражающийся от верхней стенки вторичной волной сжатия. Последующее расширение за критическим сечением приводит к отрыву газа от поверхности и формированию в донной области возвратного течения, замыкание которого осуществляется при его взаимодействии с отраженной волной сжатия.

Известно, что на течения отрывного характера сильное влияние оказывает вязкость, в рассматриваемой задаче аппроксимационная. Сравнительный анализ рис. 3,6, в, г показывает, что пониженный уровень схемной вязкости, присущий схемам повышенного порядка, ускоряет момент отрыва, а также качественно ослабляет «размывание» газодинамических разрывов по пространству и удерживает их фактически на двух расчетных ячейках, не вызывая паразитных осцилляций. При этом изменяются точки отрыва и присоединения нулевой линии тока возвратного течения, уточняются местоположения ударной волны и ее первичного отражения. Можно также видеть на рис. 3, а еще одно отражение волны сжатия, которая после , дифракции на отрыве практически не фиксируется на рис. 3, ви г.

Вопрос о правдоподобии моделирования чисто вязких явлений в рамках уравнений Эйлера не рассматривается в данной статье. Хотя этот вопрос долго дискутируется в литературе, начиная с ранних работ [7], [8] и заканчивая [9—12], он уже не носит научного характера.

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

Рис. 3

Так на рис. 4, а—в приведены зависимости давления газа, его плотности и числа М на' теле от координаты х, для удобства индексация кривых и номер эксперимента совпадают. Можно отметить, используя данные на кривой 3 как эталонные, что точки кривой 2 (для улучшенной сетки) в основной массе расположены заметно ближе, хотя максимальные отклонения во всех экспериментах невелики и для величины давления составили приблизительно 0,5%, а для величины плотности и числа М — 2,0%. Используя тот факт, что из уравнений

*+ д^(рУ) = О,

др(

а;

+ д^(рЯУ) = О.

следует сохранение в стационарном пределе вдоль линии тока полной энтальпии (например\ вдоль твердой границы — нулевой линии тока), то по ве-

уу _ ^

личине —-—22- вдоль тела можно оценить погрешность установившегося ре-

J8

25

28

(5

10

5

J____L

i l I

0.)

2.0

<6

(,'1

1,2

11.8

8.6

0.4

-Ofi -0.5 -0,4 -0.3 -0,2 -0.1 Й f) -Ofi -0,5 -0,4 -0,3 -0,2 -0,1 a

- , . 4.040 - /Л

- A 0,035 - / \

_ в/ 0,030 ms : ^

- /у 0,020 -

0,015

~\I I I I 0,010 -BOOS - _i_ i i i i i

3) -0,0 -0.5 -Oft -0,3 -0.2 -0,1 0 'gj -Ofi -0,5 -0.4 -0,3 -0,2 -0,1 0

а) — зависимость давления на теле от х; б) — зависимость плотности на теле от х; в) — зависимость числа Маха на теле от х; г) — зависимость погрешности расчета энтальпии торможения на теле от х

Рис. 4

шения. На рис. 4, г приведены соответствующие зависимости такой характеристики.

Для схемы первого порядка применение сеток с улучшенными параметрами ортогональности, как видно по зависимостям 1 и 2, приводит К снижению погрешности расчета, все же заметно уступая схеме второго порядка, ошибка которой фактически ограничивается 1,0%.

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

ЛИТЕРАТУРА

1. G о r d оп W. J. and Т h i е I L. С. Transfinite mappings and their applications to grid generation. Оп «Numerical Grid Generations ed. Ьу J. F. Thompson.

2. Е r i k ss о n L. Е. Generation of boundary-conforming grids around wing-body configurations using transfinite interpolation.— AIAA J. vol. 20, N 10, October 1982.

3. S m i t h R. Е. Algebraic grid generation. On «Numerical Grid Generation:. ed. Ьу J. F. Thompson.

4. Г о д у н о в С. К, 3 а б р о д и н А. В, И в а н о в М. Я., К р а й-

к о А. Н, Про к о п о в Г. П. Численное решение многомерных задач газо-

вой дннамикн.— М.: Наука, 1976.

5. К о л г а и В. П. Коиечно-разностная схема для расчета двумерных разрывных решений нестационарной газовой динамики.— Ученые записки

АГИ, 1975, т. 6, № 1.

6. Т и л л я е в а Н. И. Обобщение модифицированной схемы Годунова

на произвольной нерегулярной сетке.— Ученые записки ЦАГИ, 1986, т. 17, № 2.

7. Фо н а рев А. С., Кол г а н В. П. Численный расчет дифракции ударной волиы на сфере и цилиндре и установление стационарного обте-каиия.— Труды ЦАГИ, 1971, вып. 1324.

8. Фо н а р е в А. С. Расчет обтекания осесимметричных тел и несущих крыловых профилей трансзвуковым потоком газа.— Ученые записки ЦАГИ, 1973, Т. 4, Nt 3.

9. Ш м ы г л е в с к и й Ю. Д., Щ у р ш а л о в Л. В. О непригодности подмены решений уравнений Навье — Стокса численным решением уравнений Эйлера.— Ж. вычисл. матем. и матем. физ, 1989, т. 29 № 10.

10. К u m а r А, S а 1 а s М. D. Euler and Navier — Stokes solutions for supersonic shear flow past а circular cylinder.— AIAA J., vol. 23, N 4, 1985.

11. N е w s о те S. R. А comparison of Euler and Navier — Stokes solution for suрегsопlc flow over а соп1са1 delta wing.— AIAA j., vol. 24, N 4, 1986.

12. К э д и л О. А, Ч ж у а н А. X. Влияние численной диссипации на результаты pac'lera вихревых течений на основе уравнений Эйлера.— Аэрокосмическая техника, 1988, № 11.

Рукопись поступила 3//V /990 г.