МЕТОД ПОГРУЖЕННОЙ ГРАНИЦЫ ДЛЯ РАСЧЕТА СВЕРХЗВУКОВОГО ОБТЕКАНИЯ ЗАТУПЛЕННЫХ ТЕЛ НА ПРЯМОУГОЛЬНЫХ СЕТКАХ Текст научной статьи по специальности «Физика»

УДК 532.516.5

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

*

затупленных тел на прямоугольных сетках

В.В. Винников, Д.Л. Ревизников

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

Введение

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

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

Алгоритм метода погруженной границы реализован в виде программы для решения нестационарной внешней задачи аэродинамики и является составной частью программного комплекса, моделирующего воздействие сверхзвукового гетерогенного потока на обтекаемое тело. Решению этой сложной задачи в различных постановках посвящено множество работ [10-18]. * Работа выполнена при поддержке РФФИ (проект № 05-08-01478-а).

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

Постановка задачи

Рассматривается двумерная стационарная задача сверхзвукового поперечного обтекания кругового цилиндра радиуса ^о потоком невязкого сжимаемого идеального газа. Расчетная область, приведенная на рис.1., представляет собой прямоугольник шириной ^ и высотой h . Центр цилиндра расположен в нижнем правом углу расчетной области и имеет координаты (^,0) .

Область покрывается прямоугольной равномерной сеткой с шагами 5х = (Их — 1),

¿у = И/ (Ну — 1). Центры ячеек сетки имеют координаты (г— 0.5)-¿х, (у — 0.5)-8у, где г = 0,1,..., Ых, у = 0,1,..., Му .

В начальный момент времени область заполнена газом плотности р0 под давлением Р0 . Газ, занимающий подобласть , имеет горизонтальную скорость и0 . Газ в подобласти находится в состоянии покоя.

Решение задачи производится методом установления. Система уравнений Эйлера, описывающая течение газа, в декартовой системе координат имеет вид 8q , 8! (ч) + 8О(д)

дг

8х

8у

= 0

где вектор консервативных переменных и потоки задаются выражениями

ч =

Г рЛ

ри ру

рЕ

! =

( ри Л

2

ри + р риу риН

О =

( ру л

риу

2

ру + р руН

Система уравнений дополняется соотношениями для полной энергии и энтальпии:

Е =

Р

р(7— 1) 2'' ' р

Для замыкания системы уравнений на границах расчетной области ставятся следующие краевые условия. На левой границе (х = °) задается краевое условие первого рода, определяющее

1 (и2 + у2), Н = Е

+

Р

сверхзвуковой невозмущенный поток. Значения физических переменных на правой = w) и верхней (у = Ь) границах находятся из однородного условия Неймана, аппроксимирующего условие сверхзвукового выхода из области. На нижней границе (у = 0) задается условие симметрии. На поверхности цилиндра ставится условие непротекания.

3.0

Метод расчета

Для численного решения уравнений Эйлера имеется богатый набор разнообразных методов, среди которых ни один не обладает подавляющими преимуществами над остальными. Из всего множества сеточных методов для решения гиперболических систем были выбраны TVD-монотонизированные варианты методов Куранта-Изаксона-Риса (CIR) и Хартена-Лакса-ван Лира (HLL), см., например, [19]. Эти методы позволяют производить расчет задачи сверхзвукового обтекания затупленных тел, в то время как более точные методы типа Roe или Хартена-Лакса-ван Лира с выделением контактных разрывов (HLLC) менее универсальны, поскольку проявляют на ряде задач численную неустойчивость [20,21] и требуют привлечения дополнительной процедуры энтропийной коррекции. В ходе проведенного сравнения было получено, что для рассматриваемого класса задач метод Хартена-Лакса-ван Лира существенно превосходит по производительности метод Куранта-Изаксона-Риса, не уступая последнему по точности получаемого решения.

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

К

их

i+1/2

К

ЛЛ-ЛЛ (дк —дь)

4 — Ль К»,

0 <АГ

Л < 0 < Л Л < 0

где ^ = К д), ^ = К (дк).

Реконструкция векторов консервативных переменных на гранях ячейки выполнялась согласно выражениям:

дь = д, + °.5 • Н^^д,—1, Дд,), дк = д1+1 — 0.5 • Нткш^Дд^, Ад,),

где Дд, = д, — дг-1, а функция Нткег() покомпонентно применяет одну из функций-ограничителей к паре векторных аргументов. В настоящей работе в качестве основного использовался лимитер minnod(•):

• ,ч ísign (а)• ш1п (а, Ы1 аЬ > 0

тштоо(а, Ь ) = ^ ^ 111

4 7 [ 0, аЬ < 0'

Существенных отличий в результатах решения при использовании других лимитеров отмечено не было.

Собственные значения Ль, на вертикальных гранях ячейки вычислялись по вектору

и 1 у т-* I * * * т-1* \

усредненных значений физических переменных и = Р и V Е ):

Ль = и* — с, Ля = и* + с, с = МЕ* —

И4к!}(,—1),

где

р =

Рь + Ря и * = Рьиь + Ряия V = РЛ+Ря^ Е * = РьЕь + РпЕК

V =

2 Рь +Ря Рь + Ря Рь +Ря

Также можно проводить Roe-усреднение физических переменных:

* /- . иыРь + иял Ря . vьл Рь + МРя

Р = л РьРя , и =-/=-Г=— , v =-Г=-Г=

лРь +^Ря лРь +^Ря

■яМУя ,,* = уь'\1Нь 1 уя'\1Ня Е* = Еьч Рь + Ея"<1 Ря

Вычисление потоков вдоль оси у через грани ячейки по методу Хартена-Лакса-ван Лира осуществлялось аналогичным образом:

О

их

У+12

Ов

0 <лв

Хт°в —Хв°т +Хв1т(Ят—Чв) ^ < 0

'в "в^т

Хт Лв

От, Л < 0

где О в = о(Чв ), От = о(Чт ),

Чв = + 05 - Пткег(Дчу.—1з Дчу ), чт = чу+1 — 0.5 - Пт^ег^Дч^, Дчу),

ДЧу = Чу — Чу—!, Хв = У*— с, Л = у* + с .

Аппроксимация краевых условий

Для аппроксимации краевых условий на криволинейной границе используются три узла (х1, у1), (х2, у2 ), (хО, уО ) и одна точка на границе (х0, у0 ), как показано на рис.2.

Рис.2. Используемые при дискретизации краевых условий узлы расчетной сетки и вспомогательные точки.

На криволинейной границе твердого тела для уравнений Эйлера реализуется условие непротекания. Согласно этому условию нормальная составляющая скорости ип равна 0:

ип0 = и0 со$О + у0 8шО = 0 .

где О - угол наклона нормали к поверхности. Тангенциальная составляющая скорости равна:

и о = —и0 зтО + у0 соъО

Значения величин и0 и у0 в точке на границе (х0, у0 ) связаны билинейным интерполяционным соотношением с известными значениями в двух приграничных узлах (х1, у1), (х2, у2 ) и искомым значением в фиктивном узле (хО, уО ):

uo =(1 XO Уо)

f1 X yл

(1 yo) 1 X У 2

v1 XG Уо

(1 X1 y1 ^ -1 ( u11 ( Г X1 У1' -1 ( v1 N

1 X2 у 2 u2 , vo =(1 Xo Уо J1 X2 У2 v2

V1 xg yG j v ug j V1 xg yg j v vg j

( «Л

Ы-,

v «G J

(v Л

cos в +

V VG J

sine

= 0.

Пусть

X y^"

(b1 b2 bg)=(1 XO Уо)

1 X2 У 2

V1 xg yg

тогда выражение можно переписать в виде

uG cos в + vG sin в = -—(u1 cos в + v1 sin в) -—(u2 cos в + v2 sin в)

bG bG

Для однозначного определения искомых значений скоростей uG , vG в фиктивном узле (xG, yG ), необходимы дополнительные уравнения.

Для замыкания системы используем краевое условие = 0 на криволинейной границе в точке

дп

O. Аппроксимация этого краевого условия также осуществляется с помощью процедуры билинейной интерполяции:

b2

и.

sin в - vG cosв = (1 Пусть

(X d 2 do )=(1 xg yg

yG )

(1 X1 y1

v

1 X2 y 2

0 cos в sin в

^ 1( u1sinв- v1cosв> u2 sin в - v2cosв

(дит /дп )

o j

1 x1 1 x2

y1

y2

v0 cos в si^

тогда

в - vG cos в = d1 (u1 sin в - v cos в) + d2 (u2 sin в - v2 cos в) + dO (дит/дп)о .

uG sin в - v G ■

Введем вспомогательные обозначения:

b

RHS1 = —L (u1 cos в + v1 sin в)—- (u2 cosв + v2 sin в)

bG bG

RHS 2 = d1 (u1 sin в - v1 cosв) + d 2 (u2 sin в - v2 cos в).

Получим два уравнения:

1

1

v

2

1

uG+1cos0+vG+1sin9 = Щ,

G k+1

G k+1

sinfl-v^1 cos0 =rhs2.

Откуда можно получить:

4t1=M,cos0+RHS2sinS, vlG " = RHSlsinS- RHS,cosS.

Величина pa определяется из краевого условия на границе дп = о, которое

аппроксимируется линейным соотношением:

Pg =(1 xg yG J 1

X У1

X2 y2 0 cos 9 sin 9

Pi P2 {dp/dn )(

о у

Значение Еа определяется из краевого условия для давления на границе др/дп = {p{uт^fУR :

Eg =

Pg

+

Pg (Ï- 1) 2

1 (uG+vG ),

f 1 x1 У ^ -1 f P1 ^

Pg =(1 xg yg ) 1 X2 У 2 P2

v0 cos 9 sin9y v(dp/dn)o j

Величину uTQ в точке O на границе можно получить из уравнения:

uo =(1 xo Уо )

f1 X У1 1 X2 У 2 0 cos 9 sin 9

^ 1f u1sin9- v1cos9>| u2 sin 9 - v2cos9 (duг/ dn)c

^WO^/ Oil! KS J

а величина pO может быть определена из уравнения Сл „ л-1/ „ л

Ро =(1

xr

Уо )

У1

У 2

v0 cos 9 sin9y v

Р1 Р 2 (dp/dn )

о

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

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

(la ^ P ]

Soa = Ua — Uo , где Ua = ua = uo

va vo

E V^a у E V^ o у

Вводится дополнительная точка-образ I с координатами {x:, yj), которые определяется

выражением:

Г Xj ]

V У: у

= 2

Г x ^

V Уо у

Г XG ^

v yg у

Далее на отрезке 10 с помощью процедуры интерполяции вычисляется вектор наклонов физических переменных по их значениям в узлах ,

Ух ), (х2> У 2 ) (х5> У 5 ) :

.Г 1 1 1 1 1 ТУ 0 >

(

SIO =

Pi р2 Рз Ра

u, u1

v1 v2

V Ei E2

u

u„

4 E4

P5

U5

v5

E,

5 у

X,

X4

x5

У, У 2 Уз Уа У 5

X

X2 X3

X4 X5

22222 V У, У2 Уз Уа У5 у

XO Xj

Уо — У1 22 XO — Xj 22 V Уо — У: у

Полученные наклоны Soa и SI0 используются для восстановления окончательного значения вектора физических переменных иа в фиктивном узле (ха, уа ): Н^М^О > SI0 ).

Ua = Uo

v

Результаты расчетов

Решалась модельная задача поперечного обтекания кругового цилиндра единичного радиуса. Число Маха M о набегающего свободного потока задавалось равным 6. Ширина и высота расчетной области составляли w = 1-5 •R и h = 3 0 • R0 соответственно. В ходе проведения вычислительного эксперимента определялись поля плотности, скорости и давления, обезразмеренные относительно величин р0, *Jp0 /р0 , p0 соответственно. Расчет производился на равномерной прямоугольной сетке 128 х 256 узлов. Изолинии полей плотности и чисел Маха представлены на рис.3-а,Ь. По распределению плотности также численно смоделирована теневая фотография, приведенная на рис.3-с. Нормированная яркость каждого зерна теневой фотографии

задавалась выражением, приведенным в работе [22]: BSchlieren = {l — |Vp|/max|Vp))15.

a) b) c)

Рис.3. a) - Изолинии полей плотности; b) - изолинии чисел Маха; с) - Симуляция теневой фотографии (Schlieren image) и приближенно-аналитическое положение головной ударной волны.

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

x(y) = -R0 - А + RC cot(

(

1 +

У

tan2 (()

,0.5

-1

где

sin(() =

y J Mo '

RC = R0 -1.386exp

1.8

(Mo -1)0 '

А = R0 • 0.386exp

f 4.67 ^

V M0 у

Результаты сравнительного анализа численных расчетов с эталонными распределениями, приведенными в работе [23], показаны на рис.4-6. Поскольку расчетные поля были получены на прямоугольных сетках, была выполнена процедура переинтерполяции в цилиндрическую систему координат. Это внесло некоторую погрешность в значения численных полей вблизи ударной волны. Кроме того, было отмечено несовпадение численных распределений с эталонными вблизи поверхности цилиндра. Это объясняется грубостью расчетной сетки. Расчеты, проведенные на

более подробной сетке 512 х 1024, хорошо согласуются с эталонными распределениями на поверхности цилиндра.

Рис.5. Радиальные сечения расчетного и эталонного полей давления.

Расчет: ■

0o; 20o;

45o

60o

90o

M

Эталон: □ 0o; о 20o; a 45o; v 60o; о 90o

\

—i— -3.0

—1— -2.5

-2.0

-1.5

-1.0

Рис.6. Радиальные сечения расчетного и эталонного чисел Маха.

Заключение

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

5

4

3

2

Список литературы

1. Koh E.P.C., Tsai H.M., Liu F. Euler Solution Using Cartesian Grid with a Gridless Least-Squares Boundary Treatment // AIAA journal, 2005, Vol. 43, № 2, pp. 246-255.

2. Kirshman D.J., Liu F. A gridless boundary condition method for the solution of the Euler equations on embedded Cartesian meshes with multigrid // Journal of Computational Physics, 2004, Vol. 201, № 1, pp. 119-147.

3. Kirshman D.J., Liu F. Flutter prediction by an Euler method on non-moving Cartesian grids with gridless boundary conditions // Comput. fluids, 2006, Vol. 35, № 6, pp. 571-586.

4. Luo H., Baum J.D., Lohner R. A hybrid Cartesian grid and gridless method for compressible flows // Journal of Computational Physics, 2006, Vol. 214, № 2, pp. 618-632.

5. Piquemal A.-S., Chiavassa G., Donat R. A Brinkman-penalization method for compressible viscous flows // Submitted to Computers and Fluids, 2004.

6. Colella P., Graves D.T., Keen B.J., Modiano D. A Cartesian grid embedded boundary method for hyperbolic conservation laws // Lawrence Berkeley National Laboratory. 2004, Paper LBNL-56420.

7. Винников В.В., Ревизников Д.Л. Применение декартовых сеток для решения уравнений Навье-Стокса в областях с криволинейными границами // Математическое Моделирование,

2005, т.17, №8, с. 15-30.

8. Винников В.В., Поликша И.В., Ревизников Д.Л. Применение метода погруженной границы к решению задач теплообмена с подвижным фронтом фазового перехода // Труды Четвертой Российской национальной конференции по теплообмену: в 8 томах. Т.7. Радиационный и сложный теплообмен. Теплопроводность, теплоизоляция. - М.: Издательский дом МЭИ,

2006, с. 183-186.

9. Винников В.В., Ревизников Д.Л. Компьютерный практикум по численному решению задач математической физики в областях с криволинейными границами // Сборник статей "Информационные технологии и программирование", Выпуск 1 (13). - М.:, МГИУ, 2005, с. 5-20.

10. Михатулин Д.С., Полежаев Ю.В., Ревизников Д.Л. Исследование разрушения стеклопластика при полете в запыленной атмосфере // Теплофизика высоких температур, 2001, Том 39, № 4, с. 640-648.

11. Михатулин Д.С., Полежаев Ю.В., Ревизников Д.Л. Исследование разрушения углеродного теплозащитного материала при полете в запыленной атмосфере // Теплофизика высоких температур, 2003, Том 41, № 1, с. 98-105.

12. Михатулин Д.С., Ревизников Д.Л., Способин А.В., Шехтер Ю.Л. Особенности теплоэрозионного разрушения материалов в сверхзвуковом полидисперсном потоке // Труды Четвертой Российской национальной конференции по теплообмену: в 8 томах. Т.6. Дисперсные потоки и пористые среды. Интенсификация теплообмена. - М.: Издательский дом МЭИ, 2006, с. 87-90.

13. Волков А.Н., Циркунов Ю.М. Влияние дисперсной примеси на течение и теплообмен при поперечном обтекании цилиндра сверхзвуковым потоком запыленного газа // Изв. РАН. МЖГ. 2005. № 4. с. 68-85.

14. Стасенко А.Л. Физическая механика многофазных потоков. - М.: Изд-во МФТИ, 2004. -136 с.

15. Осипцов А.Н., Егорова Л.А., Сахаров В.И. Аэродинамическая фокусировка полидисперсных частиц при обтекании тел запыленным газом // Доклады Российской Академии наук, 2004, т. 395, № 6, с. 767-772.

16. Головачев Ю.П., Шмидт А.А. Обтекание затупленного тела сверхзвуковым потоком запыленного газа // Изв. АН СССР МЖГ, 1982, № 3, с. 73-77.

17. Chang S. S.-H. Nonequilibrium phenomena in dusty supersonic flow past blunt bodies of revolution // Physics of Fluids, Vol. 18, № 4, 1975, pp. 446-452.

18. Elangovan R., Cao H.V. Dusty supersonic viscous flow over a two-dimensional blunt body // Journal of Thermophysics and Heat Transfer, Vol. 4, № 4, 1990, pp. 529-533.

19. Куликовский А.Г., Погорелов Н.В., Семенов А.Ю. Математические вопросы численного решения гиперболических систем уравнений. - М.: ФИЗМАТЛИТ, 2001.- 608 с.

20. Dumbser M, Moschetta J.M, Gressier J. A matrix stability analysis of the carbuncle phenomenon // Journal of Computational Physics 2004; Vol. 197, № 2, pp. 647-670.

21. Chauvat Y., Moschetta J.-M., Gressier J. Shock wave numerical structure and the carbuncle phenomenon // Int. J. Numer. Meth. Fluids 2004, Vol. 47, № 8-9, pp 903-909.

22. Hagen T.R., Henriksen M.O., Hjelmervik J.M., Lie K.-A. How to Solve Systems of Conservation Laws Numerically Using the Graphics Processor as a High-Performance Computational Engine // submitted for publication in "Geometric Modelling, Numerical Simulation and Optimization: Industrial Mathematics at SINTEF", Springer Verlag.

23. Любимов А.Н., Русанов В.В. Течение газа около тупых тел: в 2 ч. Ч. 1. - М.: Наука, 1970. -287 с.

Винников Владимир Владимирович, научный сотрудник кафедры вычислительной математики и программирования Московского авиационного института (государственного технического университета), к.ф.-м.н.; e-mail: vvinnikov@list.ru; контактный телефон: 158-48-94;

Ревизников Дмитрий Леонидович, профессор кафедры вычислительной математики и программирования Московского авиационного института (государственного технического университета), д.ф.-м.н; e-mail: reviz@k806.mainet.msk.su; контактный телефон: 158-48-94.