Дозвуковое обтекание профилей Текст научной статьи по специальности «Физика»

(l + т2с2 ))(i'/+1) - у(| < 2т2 (1 + п)k(и, u)q''-1||(р^1 - (p(j-1)\ . Заметим, что при 1 + 2л < jq условие (5) выполняется при любом Т, т.е., начиная с временного слоя j = 8, в любом случае может быть указано 0 < q < 1, независимо от Т. Для последующих временных слоёв значение q может быть уменьшено, что означает ускорение сходимости. Это полностью подтверждается практикой численного моделирования. На начальных итерациях процесса установления невязка может возрастать, затем наблюдается монотонное убывание, причём от временного слоя к слою убывание ускоряется.

БИБЛИОГРАФИЧЕСКИЙ СПИСОК

1. Васильев В.С., Сухинов А.И. Прецизионные двумерные модели мелких водоемов // Математическое моделирование. - 2003. - Т. 15, № 10. - С. 17-34.

2. Филатов Н.Н. Гидродинамика озер. - СПб.: Наука, 1991. - 200 с.

3. Богданов К.Т. Приливные явления в Тихом океане. - М.: Наука, 1994. - 144 с.

Статью рекомендовал к опубликованию д.т.н., профессор Я.Е. Ромм.

Васильев Владислав Сергеевич - Федеральное государственное автономное образовательное учреждение высшего профессионального образования «Южный федеральный университет»; e-mail: sunspire@mail.ru; 347928, г. Таганрог, пер. Некрасовский, 44; тел.: 88634371606; кафедра высшей математики; к.ф.-м.н.; доцент.

Vasiliev Vladislav Sergeevich - Federal State-Owned Autonomy Educational Establishment of Higher Vocational Education "Southern Federal University"; e-mail: sunspire@mail.ru; 44, Nekrasovsky, Taganrog, 347928, Russia; phone: +78634371606; the department of higher mathematics; cand. of. eng. sc.; associate professor.

УДК 534.22

И.Э. Гамолина, В.В. Дурягина, В.В. Семенистый ДОЗВУКОВОЕ ОБТЕКАНИЕ ПРОФИЛЕЙ

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

Аэродинамика; дозвуковой поток; система уравнений Навье-Стокса; схема расщепления; метод дробных шагов; крыловой профиль.

I.E. Gamolina, V.V. Duryagina, V.V. Semenistyi SUBSONIC PROFILES FLOW

The numerical modeling of an in viscid flow of plane wing profiles in up to sound operational modes is considered. The complexity of aerodynamic processes with possible circulating zones of a separation of a boundary layer requires qualitative calculation with the purpose of deriving the improved designer wing performances. The used finite difference scheme with decomposition of a difference operator on space and physical processes precisely approximates the ini-

tial continuous task. The numerical score is carried on in physical space with a consequent transformation of coordinates in calculated rectangular area before installation on time and exit on a fixed flow regime. The results of numerical experiment for pressure coefficient are considered.

Aerodynamic; subsonic flow; Navier—Stokes equations; difference scheme with decomposition fractional step method; wing profile.

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

Аэродинамика различных частей летательных аппаратов (крыльев, оперений и т.д.) в значительной мере определяется формой их сечений в направлении набегающего потока, т.е. формой их профилей [1]. Это связано с тем, что аэродинамические характеристики профилей крыльев и оперений должны обеспечить летательному аппарату при заданной подъемной силе наименьшее сопротивление, необходимую устойчивость и управляемость при всех режимах полета.

Остановимся более подробно на постановке задачи и ее конечно-разностном алгоритме решения.

Для описания течения возле крылового профиля рассматривалась начально-

краевая задача для полной системы уравнений Навье-Стокса:

+ (1)

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

Искомый вектор / выбирался в переменных р,и1,и2,и3 и Т.

После преобразования координат х± = x1(q1,q2,q3); Х2 ~ Х2 {ЧъЧг.Чз); хз = хз (Л, Ч2,Чз) физическая область течения переходит в стандартную математическую область. В математической области [0,1] хТ (Т = {а^ < qt < biti = 1,2,3} — параллелепипед) с равномерными узлами по пространственным переменным и временному (итерационному) параметру т для аппроксимации системы (1) использовалась следующая конечно-разностная схема (предложенная в работе [2]):

fn+l_fn J

+ + (1-х)/"] = Fn. (2)

В уравнении (2) разностный оператор, аппроксимирующий дифференциальный оператор, представлен в виде суммы расщепленных операторов, содержащих конвективные и вязкие члены (i = 1,3,5), и расщепленных операторов, содержащих члены с давлением и члены вида div V в уравнениях неразрывности и энергии (i = 2,4,6) по каждому координатному направлению. Выбор параметра (веса) х осуществляется из условия устойчивости схемы (2).

Преобразуя схему (2) к каноническому виду и затем приближенно фактори-зуя оператор (/ + тх £f=1 получаем схему дробных шагов с тем же порядком 0(тт + hk) аппроксимации, что и уравнение (2):

+ TxQ![h)^n+Le = I = 1.....6,

(3)

[fn+1=fn + T^n+1, п = 0,1,2,..., N. Здесь вектор fn = — Q.khfn + Fn, = Е?= i и

I П+4 П+-7 П+-7 n+4 n+4

^ 6 = (СР б, б, б, С„3 б, Ст б) - векторы невязок. На каждом дробном шаге решается система из пяти разностных уравнений. Эти уравнения решаются методом прогонки вдоль заданного направления. Остановимся подробно на реализации первого и второго дробных шагов.

На первом дробном шаге расчет проводится вдоль строк каждого сечения параллелепипеда Т плоскостями q3 = const (горизонтальными плоскостями). Строки параллельны оси q± и перпендикулярны оси q2. Учитывая диагональный вид оператора I + решается следующая система разностных уравнений:

1) для нахождения координаты плотности вектора невязки f в:

Т1+ та t^ + laj

ISр Jij 1 , )

№ +1) i,

i+lj + № + 1) f

i+2j.

a1+la11

hi

№ +1)

+ № + i)[f

¡-и

i-2j

= №

(4)

2) для нахождения координат скорости вектора невязки f

та; J a^ + laj

Л

n+i

_ _| I +fc n+- ,

[Eli

(5)

Выписано уравнение для нахождения первой координаты скорости вектора

1

п+-

невязки £6. Уравнения для нахождения других координат скорости аналогичные. Здесь ^ и Я = ц' — - коэффициенты вязкости.

При построении разностных уравнений использовались оператор центральной разности Лх и направленные разностные операторы

Л; =---L-L (к = 1,2),

' 2hj v ' ^

]

где Т+

j операторы сдвига вправо и влево по j-й координате.

fPi

н

Параметр а,- = + — и^ + —иу (a, ft,у = 1,2,3) влияет на выбор знака

разностной производной Л,+к (из соображений устойчивости);

3) для нахождения координаты температуры вектора невязки f

+ТХ

\ a-L+laJ

-к

71+-

Сfr)y 6 + +fc

П+ п _|„ I тк П+- 1

! "J 2 1 Ч

„ „ HkTA +

(6)

+ ^ Aj ^ ЯктЛ + ^ Л* ^ ЯктЛ

Здесь kT - коэффициент теплопроводности; Ht - коэффициенты Ламе (Я = Я1Я2Я3) и = ^ - компоненты Якобиана преобразования.

Система уравнений (4)-(6) решается с помощью пятидиагональной прогонки (при изменении знака выражения аъ благодаря наличию циркуляционных зон течения) или трехдиагональной прогонкой при постоянстве знака а±. Расчет первого дробного шага завершается.

Численное решение на втором дробном шаге системы (3) проводится снова вдоль строк сечения q3 = const области Т. Решается второе из разностных уравнений системы (3): (/ + = fn+6.

И+-) , wt

В покоординатно операторно-разностном виде система уравнений для второго дробного шага имеет вид

foH + ТХ^Л!-^)".^ + ТХ^Л!-^)".^ +

1

+ТХ

н н Z31P

н

Н п+-

Л1 —tfиз) ij

=(ц;-

n+i

п+-

НР 1 ^-'P^l]

Ь^Л^)^

п+-

= )17

(7)

(е = 1,2,3);

п+-

ад. 3 + ТХ

( Н Н 1H2y>uzJi] Н 1H3y>u3Ji] )

1

п+-= 6.

„ 2 18Р 1 ЭР , ЭГ Р deN

Здесь ат = - —, oz = -— и SA = — (— о —).

рдр р дТ де vp ^ Эру

В отличие от первого дробного шага (где система уравнений (4)-(6) может решаться независимо друг от друга) при решении системы (7) сначала выписываются уравнения для нахождения компонент скорости вектора невязки f з, а затем рассчитываются компоненты плотности и температуры вектора невязки.

На третьем и четвертом дробных шагах расчеты проводятся по столбцам сечения q3 = const области Т. Эти расчеты совершенно аналогичны тем, что проводились соответственно на первом и втором шагах. Изменения в формулах связаны только с заменой направления интегрирования.

На пятом и шестом дробных шагах расчеты проводятся вдоль всех столбцов области Т, параллельных q3. Число таких столбцов совпадает с размерностью внутренних узлов сетки в сечении q3 = const. После этого расчет на итерации заканчивается. Восстанавливается вектор решения на новом временном слое и проверяется критерий установления. Если критерий не выполняется, начинается вычисление на новой итерации. Полученное решение составляет правую часть уравнения (3).

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

Л

ql

4 x2 6

Рис. 1. Расчетная и разностная сетки

в

В задаче использовалась сетка, предложенная в работе [3]. Для ее построения минимизируется нелинейная целевая функция вида:

р = ВД1к[аОЯ7иМ) + (1- а)5М(Х;,Щ

Здесь 5М(¿,у, к) - мера гладкости, ОЯТ(1,], к) - мера ортогональности сетки. Метод позволяет строить гладкие, локально-ортогональные сетки. Ниже приведены графики расчета изменения коэффициента давления на профиле крыла в зависимости от числа Маха М, угла атаки а, количества итераций п, коэффициента сопротивления Сх и коэффициента подъемной силы Су.

Рис. 2. Результаты численного эксперимента

БИБЛИОГРАФИЧЕСКИЙ СПИСОК

1. Пинчуков В.И. О численном исследовании трансзвуковых турбулентных течений возле крыла неявными схемами высоких порядков // Вычислительные технологии. - 2006. - № 2. - С. 110-121.

2. Ковеня В.М., Яненко Н.Н. Метод расщепления в задачах газовой динамики. - Новосибирск: Наука, 1981. - 304 с.

3. Пинчуков В.И. Об одном методе построения адаптивных разностных сеток в задачах аэродинамики // Моделирование в механике. - 1988. - № 2. - С. 133-141.

Статью рекомендовал к опубликованию д.ф.-м.н., профессор А.И. Жорник.

Гамолина Ирина Эдуардовна - Федеральное государственное автономное образовательное учреждение высшего профессионального образования «Южный федеральный университет»; e-mail: iegam@rambler.ru; 347928, г. Таганрог, пер. Некрасовский, 44; тел.: +79185190837; к.т.н.; доцент.

Дурягина Вероника Владимировна - e-mail: vepanuka@mail.ru; тел.: +79604426599; старший преподаватель.

Семенистый Владимир Васильевич - e-mail: vlad60sem@gmail.com; тел.: +79282135206; к.ф.-м.н.; доцент.

Gamolina Irina Eduardovna - Federal State-Owned Autonomy Educational Establishment of Higher Vocational Education "Southern Federal University"; e-mail: iegam@rambler.ru; 44, Nekrasovskiy, Taganrog, 347928, Russia; phone: +79185190837; cand. of eng. sc.; associate professor.

Duryagina Veronika Vladimirovna - e-mail: vepanuka@mail.ru; phone: +79604426599; senior lecturer.

Semenistyi Vladimir Vasilievich - e-mail: vlad60sem@gmail.com; phone: +79282135206; cand. of eng. sc.; associate professor.