CC BY
31
УДК 532.529:536.24
МЕТОД МОДЕЛИРОВАНИЯ КРУПНЫХ ВИХРЕЙ В ПРИЛОЖЕНИИ К ПРОБЛЕМАМ ВНУТРЕННЕЙ ГАЗОДИНАМИКИ РДТТ
К.Н. ВОЛКОВ, В.Н. ЕМЕЛЬЯНОВ
Балтийский государственный технический университет "ВОЕНМЕХ" им. Д.Ф. Устинова, Санкт-Петербург, Россия
АННОТАЦИЯ. Проводится моделирование крупных вихрей турбулентных течений в каналах с распределенным вдувом. Течения описываются на основе фильтрованных по пространству уравнений Навье-Стокса. Вихревая подсеточная вязкость рассчитывается при помощи И^Ю-модели с поправками на кривизну линий тока. Обсуждаются особенности постановки начальных и граничных условий на проницаемой поверхности канала. Расчеты выполняются для различных скоростей вдува с нижней и верхней стенок канала и сравниваются с данными прямого численного моделирования, решением, полученным на основе к-е модели турбулентности, и данными физического эксперимента.
ВВЕДЕНИЕ
Течение в канале с распределенным вдувом служит моделью течения продуктов разложения твердого топлива в камере сгорания ракетного двигателя на твердом топливе (РДТТ), что отражает наиболее существенную сторону процесса - подвод массы со стороны горящей поверхности заряда. Исследованию течений в каналах с проницаемыми стенками уделяется достаточно большое внимание в литературе [1-14].
Работы [1-4] посвящены экспериментальному изучению режимов течений в каналах со вдувом при различных оформлениях поперечного сечения канала (квадратное, круглое, звездообразное). Исследование [4] ограничивается достаточно малыми числами Рейнольдса, соответствующими ламинарному режиму течения.
Для численных расчетов применяются физико-математические модели различной степени сложности [1, 2, 5-8], реализуемые при помощи конечно-разностных или ко-нечно-объемных методов. Учитывая физические особенности течения, удается добиться существенного упрощения решения задачи [5-7]. Построение упрощенных математических моделей, вычислительная эффективность которых достигается за счет пре-
небрежения влиянием некоторых факторов, обосновывается соответствующими оценками и сравнением различных подходов [7].
Уравнения, описывающие течение вязкой несжимаемой жидкости между двумя параллельными пластинами, с одной из которых производится вдув, а другая является непроницаемой, допускают точное решение [5, 6]. Предполагая, что продольная составляющая скорости зависит от координаты х по линейному закону, а поперечная составляющая скорости не зависит от продольной координаты, расчет распределения скорости сводится к интегрированию уравнения обыкновенного дифференциального уравнения. В случае сильного вдува (Яе—>со) решение получается в конечном виде.
Приближенные модели [1, 2, 7] основаны на построении профиля продольной скорости из модели слоистого течения (решается уравнение Пуассона с постоянным членом в правой части) и последующим определением перепада давления из расчета потери импульса на разгон того расхода, который подведен вдувом.
Решение, описываемое соотношением для сильного вдува, хорошо согласуется с данными измерений [1-4] и допускает обобщение на случай течения вязкой несжимаемой жидкости в круглом канале с проницаемыми стенками, поперечное сечение которого изменяется во времени [8].
Для внутренних течений в РДТТ характерно турбулентное движение рабочего тела, а процессы сложной физической природы протекают на фоне общей газодинамической обстановки в рабочем пространстве РДТТ. Турбулентность сказывается как на газодинамических характеристиках потока, так и выступает в качестве интенсифицирующего фактора теплообмена в предсопловом объеме и сопловом блоке РДТТ. Процессы раздувания или эрозионного горения топлива связаны с особенностями турбулентного переноса вблизи горящей поверхности заряда, а работа внутренней теплозащиты определяется характером приповерхностных течений.
Данные физического и численного эксперимента показывают, что решение для сильного вдува достаточно хорошо описывает распределение скорости в турбулентном режиме при Яе>80...100 [5-7, 9, 10]. Расчет характеристик турбулентности проводится на основе уравнений к-г модели при известном распределении скорости [7]. Вместе с тем, приближение идеальной жидкости приводит к существенным погрешностям при моделировании турбулентных течений в длинных и узких каналах [10]. Невязкое решение неприменимо также для моделирования течений в быстрогорящих каналах [8].
Для замыкания уравнений Рейнольдса привлекаются различные модели турбулентности.
Стандартная к-е модель турбулентности позволяет получить достаточно точные результаты для каналов с двухсторонним вдувом при различных оформлениях поперечного сечения [7].
Для каналов с односторонним вдувом к-е модель не дает удовлетворительного предсказания точки перехода, а также уровня турбулентных пульсаций скорости [11, 12] (их величина в окрестности проницаемой стенки возрастает при увеличении скорости вдува). Вблизи непроницаемой стенки рассогласование расчетных и экспериментальных данных по интенсивности турбулентности достигает 15...20% [13].
Модели турбулентности 3-го и 4-го порядка, например, V2-/модель, позволяют получить результаты, согласующиеся с данными прямого численного моделирования [10, 14] (расчеты проводились для течения в плоском канале при наличии однородного вдува через одну из стенок и отсоса через противоположную). Средние характеристики потока слабо зависят от флуктуаций скорости на проницаемой поверхности [10].
Решение ряда практических задач, в частности, исследование эрозионного горения топлива, устойчивости течений, сформированных вдувом, учет влияния пульсаций скорости на скорость горения, осцилляций параметров потока в канале заряда, связанных с вихревыми структурами, моделирование переноса частиц конденсированной фазы и зашлаковывание участков газодинамического тракта, требует привлечения методов моделирования турбулентных течений, позволяющих рассчитывать не только средние, но и пульсационные характеристики потока.
В данной работе проводится моделирование крупных вихрей турбулентных течений в каналах с распределенным вдувом и обсуждается влияние скорости вдува на средние и пульсационные характеристики потока. Результаты расчетов, полученные для различных отношений скоростей вдува с нижней и верхней стенок канала, сравниваются с данными прямого численного моделирования, решением, полученным на основе к-г модели турбулентности, и данными физического эксперимента.
ОСНОВНЫЕ РАСЧЕТНЫЕ СООТНОШЕНИЯ
Совместим ось х прямоугольной декартовой системы координат с нижней стенкой канала, а оси у и г свяжем с его поперечным сечением (рис. 1). С нижней и верхней стенок канала осуществляется вдув со скоростями VII,) и соответственно.
Нестационарное течение вязкого сжимаемого газа течение описывается фильтрованными по пространству уравнениями Навьсе-Стока, записанными в консервативных переменных, которые формально совпадают с уравнениями Рейнольдса. Для замыкания основных уравнений используется ЯЫй-модель подсеточной вязкости [15]. Для учета эффектов, связанных с кривизной линий тока, подсеточная вязкость умножается на демпфирующую функцию, зависящую от числа Ричардсона [16]. Учет зависимости молекулярной вязкости от температуры производится на основе закона Сазерленда.
С
С
Р
к
А ¿1
£
а*
Рис. 1. Система координат
Ширина фильтра А связывается с размером шага разностной сетки д = у1/3 = (АхАуДя)1'3 , где V - объем ячейки, Ах, Ау, Аг - шаги сетки по координатным
направлениям х, у9 г соответственно.
В начальный момент времени задаются распределения продольной и поперечной составляющих скорости, описываемые решением для сильного вдува. Составляющая скорости в направлении оси 7 полагается равной нулю.
На нижней и верхней стенках канала при выставляются граничные условия нормального вдува (и=м>=О при у= 0 и у=к, у=ут при у=О и v=-vм^2 при у=И). Температура нижней и верхней стенки считается фиксированной (Т=Тт при у= 0 и Т=ТШ при у=И). Левая торцевая стенка канала полагается непроницаемой и теплоизолированной (и=у=м>=0, дТ/дп=0 при х=0).
На границе, через которую газ покидает расчетную область (при х=Ь), ставятся неотражающие граничные условия (условия конвективного переноса).
В направлении оси г задаются периодические граничные условия (условия повторения течения).
Метод моделирования крупных вихрей требует постановки нестационарных граничных условий, обусловленных пульсационным характером турбулентного потока.
Скорость вдува изменяется во времени по гауссовскому закону, но остается постоянной в пространстве [11].
При использовании к-е модели характеристикам турбулентности на проницаемой поверхности канала присваиваются произвольные малые значения. Характеристики турбулентности на непроницаемой стенке находятся при помощи метода пристеночных функций [7].
Дискретизация уравнений проводится при помощи метода контрольного объема на неравномерной сетке [17] и конечно-разностных схем повышенной разрешающей способности по времени и по пространству [18]. Уравнения интегрируются при помощи пятишагового метода Рунге-Кутты [17]. Для дискретизации невязких потоков применяется метод кусочно-параболической реконструкции и схема Чакраварти-Ошера, а
для дискретизации вязких потоков - центрированные конечно-разностные формулы 2-го порядка. Схема Чакраварти-Ошера основана на приближенном решении задачи Ри-мана о распаде произвольного разрыва.
Система разностных уравнений решается многосеточным методом на основе схемы полной аппроксимации (используется V-цикл и 4 уровня сетки). Последовательность вложенных сеток строится при помощи метода схлопывающихся граней (схло-пывание ячейки производится в направлении наиболее короткой грани).
Вычислительная процедура реализована в виде компьютерного кода на языке программирования C/C++. Для распараллеливания вычислительной процедуры применяется интерфейс межпроцессорного взаимодействия MPI.
Расчеты проводились на персональном компьютере типа Pentium IV CPU 2.8 GHz (для решения задачи требуется около 2.5 недель процессорного времени), на кластере типа xSeries Cluster Xeon 2x2.4 GHz (16 узлов), выпускаемом компанией Streamline Computing Ltd, и суперкомпьютере IBM SP/1600 с узлами eServer pSeries 690 на базе процессора Power 4+ 1.7GHz 6.8 GFlops (суперкомпьютерный центр расположен в Central Library of Research Council, Daresbury Laboratory, Warrington, United Kingdom).
РЕЗУЛЬТАТЫ РАСЧЕТОВ
Ширина и длина канала выбираются равными /2=0.01 м и ¿=0.60 м. В качестве рабочей среды принимается воздух.
Расчеты проводились на сетке 400x200x50 со сгущением узлов к переднему днищу (по координате х) и непроницаемой стенке (по координате у) канала по закону геометрической прогрессии. Минимальные и максимальные шаги сетки по координатным направлениям х и у составляли Axmin=1.5-10° м, ДхП1ах=1.2-10"4 м, ДушЬ=1.8-10° м, Дутах=2.0*10"1 м. Шаг сетки по координате z полагался постоянным Az=1.0-10'4 м. Шаг по времени составлял Д/=2.5-10~6 с (делалось 50000 шагов по времени).
Для удобства представления результатов линейные размеры относятся к ширине канала /?, продольная скорость - к максимальной скорости в сечении канала ит, а поперечная скорость - к скорости вдува с верхней стенки vw. Характерным параметром задачи является число Рейнольдса Rc=pvwh/[i.
Картина течения в канале, обработанная в виде в линий равных значений вихря скорости в различных координатных плоскостях, показана на рис. 2 и рис. 3 в момент времени t=0.52 с. Максимумы завихренности соответствуют центрам вихрей.
Характеристики турбулентности зависят от отношения скорости вдува к локальному значению продольной составляющей скорости vw/um. При малых параметрах вдува (vw/um~0.04) в течении доминируют мелкомасштабные вихревые структуры. При более высоких значениях (vw/um~0A) около стенки присутствуют крупномасштабные
х
Рис. 2. Вихревая картина течения в канале со вдувом при ^>-^,=0, vW2=\Л м/с (фрагмент а) и v-И;,=0? уш=5.0 м/с (фрагмент б) в координатной, плоскости
ху
о о и
У У
Рис. 3. Вихревая картина течения в канале со вдувом при vVl^1=0, уш=\.0 м/с (фрагмент а) и vvt,1=0, ^=5.0 м/с (фрагмент б) в координатной плоскости
У*
вихри. Интенсивность вихревого течения в канале усиливается при увеличении скорости вдува. Вихри, образующиеся в расчетной области, имеют вытянутую форму, что свидетельствует об анизотропности турбулентности. В свою очередь, это ставит вопрос о приемлемости моделей турбулентности, основанных на гипотезе вихревой вязкости, для исследования течений в каналах с неравномерным массоподводом.
Распределения продольной составляющей скорости по поперечному сечению канала приведены на рис. 4. Пунктирная линия показывает косинусоидальный профиль продольный скорости вихревого течения невязкой несжимаемой жидкости [5, 6], а штрихпунктирная линия - профиль скорости ламинарного течения вязкой жидкости при 11е=10\ Распределения скорости, описываемые кривыми 1-5, соответствуют различным скоростям вдува с верхней стенки ут=0Л (1); 0.55 (2); 1.0 (3); 2.75 (4); 5.0 м/с (5). Нижняя стенка канала полагается непроницаемой (ут=0). Значки • соответствуют
Рис. 4. Распределения продольной скорости при х//г=40 (фрагмент а) и х!к=50 (фрагмент б)
данным физического эксперимента [11], а значки о - результатам расчета по к-е модели турбулентности при vW2=l м/с (условия для кривой 3).
По мере увеличения числа Рейнольдса максимум продольной составляющей скорости приближается к непроницаемой поверхности, вблизи которой развивается зона вязкого течения, подобная пограничному слою на плоской пластине. При высоких числах Рейнольдса четко выражена зона сингулярного влияния вязкости (через граничное условие прилипания) и зона вихревого течения невязкой жидкости. С увеличением расстояния от левой стенки канала профиль скорости принимает более наполненную форму, чем профиль скорости, полученный на основе к-е модели турбулентности.
Переход ламинарного режима течения в турбулентный определяется по картине поведения коэффициента трения вдоль нижней непроницаемой поверхности канала (рис. 5). Коэффициент поверхностного трения находится из соотношения С^2хм,/(ри2), где Тц; - напряжение поверхностного трения, и - средняя скорость в поперечном сечении. Сплошные линии 1-5 соответствуют различным скоростям вдува с верхней стенки, а значки □, о, • - расчету на основе к-г модели турбулентности для условий, описываемых кривыми 3, 4 и 5. Сравнение полученных результатов показывает, что к-г модель дает удовлетворительное предсказание коэффициента трения лишь при достаточно малых скоростях вдува (кривая 3). При увеличении скорости вдува расхождение результатов возрастает по мере удаления от левой границы.
Резкое увеличение коэффициента трения, начиная с х!Ь~0.2, свидетельствует о ламинарно-турбулентном переходе. Положение точки перехода зависит от параметров задачи, а также оказывается достаточно чувствительным к флуктуациям скорости на проницаемой стенке и при увеличении скорости вдува смещается вниз по потоку. Координата точки перехода удовлетворительно согласуется с данными [13]. На основе полученных результаты можно выделить следующие режимы течения:
Рис. 5. Распределения коэффициента поверхностного трения вдоль непроницаемой стенки канала
- область существенного влияния сил вязкости (х/И<1/Ке);
- ламинарное течение с косинусоидальным профилем скорости (1<х//?<5);
- область перехода (5<х//к10...15);
- турбулентное течение (10... 15<х/к).
Особенностью течения в канале с распределенным вдувом является наличие отрицательного градиента давления, который оказывает существенное влияние на механизм и интенсивность турбулентного переноса.
С увеличением координаты х и скорости вдува происходит смещение максимума кинетической энергии турбулентности от стенки в поток (рис. 6). Профиль кинетической энергии турбулентности имеет два максимума, расположенные вблизи стенок канала. Увеличение уровня турбулентных пульсаций скорости наблюдается в области сильного сдвига на некотором расстоянии от проницаемой стенки, где жидкие частицы, движущиеся по нормали к поверхности, вынуждены развернуться в узкой приповерхностной зоне. Величина максимума, расположенного в окрестности непроницаемой стенки, существенно меньше (почти в 2 раза). Результаты расчетов согласуются с данными физического эксперимента [11] (значки •), исключая пристеночные зоны потока, где расчет переоценивает интенсивность турбулентности [13].
Расчеты по к-г модели (значки о) предсказывают положение максимума кинетической энергии турбулентности на большем расстоянии от проницаемой поверхности канала [8], а также переоценивают уровень турбулентности вблизи его стенок по сравнению с данными измерений и результатами прямого численного моделирования [11, 13].
Рис. 6. Распределения кинетической энергии турбулентности при х//2=40 (фрагмент а) и х/й=50 (фрагмент б)
Интенсивность турбулентности в канале с проницаемыми стенками имеет такой же порядок, что и уровень турбулентности полностью развитого турбулентного течения в канале, когда вязкие и рейнольдсовы напряжения находятся в равновесии с напряжениями, обусловленными давлением [6].
Вклад пульсаций продольной и поперечной составляющих скорости в баланс кинетической энергии турбулентности примерно одинаков (средний клад пульсаций скорости в направлении оси 2 равняется нулю). Распределения нормальных компонент тензора напряжений Рейнольдса в качественном отношении практически не отличается от распределения кинетической энергии турбулентности, в то время как профили касательных компонент меняют знак при переходе через серединную линию канала.
ЗАКЛЮЧЕНИЕ
Проведено моделирование крупных вихрей турбулентных течений в каналах со вдувом для различных комбинаций скоростей вдува с нижней и верхней стенок. Сравнение результатов численного моделирования с расчетами, проведенными на основе к-г модели, и данными физического эксперимента показывает, что метод моделирования крупных вихрей предсказывает более наполненный профиль продольной составляющей скорости по сравнению с решением, полученным на основе модели вихревого течения невязкой жидкости или к-8 модели. При этом максимум кинетической энергии турбулентности располагается ближе к проницаемой поверхности канала, чем это предсказывает к-г модель.
Разработанные средства численной диагностики и полученные результаты являются основой для моделирования движения частиц конденсированной фазы в турбулентном потоке, их инерционной сепарации и межфазного взаимодействия.
СПИСОК ЛИТЕРАТУРЫ
1. Емельянов В.Н. Физическое и вычислительное моделирование трехмерных течений в двигательных установках // Внутрикамериые процессы, горение и газовая динамика дисперсных систем. Санкт-Петербург: Изд-во БГТУ, 1996. С. 124-137.
2. Емельянов В.Н. Внутренние течения сложной структуры // Внутрикамериые процессы, горение и газовая динамика дисперсных систем. Санкт-Петербург: Изд-во БГТУ, 1998. С. 80-91.
3. Бендерский Б.Я., Тененев В.А. Экспериментально-численное исследование течений в осесимметричных каналах сложной формы со вдувом // Известия РАН. МЖГ. 2001. №2. С. 184-188.
4. Cheng Y.C., Hwang G.J. Experimental studies of laminar flow and heat transfer in a one-porous-wall square duct with wall // International Journal of Heat and Mass Transfer. 1995. Vol. 38. No. 18. P. 3475-3484.
5. Райзберг Б.А., Ерохин Б.Т., Самсонов К.П. Основы теории рабочих процессов в ракетных системах на твердом топливе. М.: Машиностроение, 1972. 384 с.
6. Численный эксперимент в теории РДТТ / Под редакцией A.M. Липанова. Екатеринбург: Наука, 1994. 304 с.
7. Волков К.Н., Емельянов В.Н. Математические модели трехмерных турбулентных течений в каналах со вдувом // Математическое моделирование. 2004. Т. 16. № 10. С.41-63.
8. Majdalani J., Vyas А.В., Flandro G.A. Higher mean-flow approximation for a solid rocket motor with radially regressing walls // AIAA Paper. No. 2001-3870. 11 p.
9. Ciucci A., Iaccarino G., Moser R., Najjar F., Durbin P. Simulation of rocket motor internal flows with turbulent mass injection // Center for Turbulence Research. 1998. P. 245-266.
10. Chaouat B. Numerical simulation of channel flows with fluid injection using Reynolds stress model // AIAA Paper. No. 2000-0992. 18 p.
11. Chung Y.M., Sung H.J., Krogstad P. A. Modulation of near-wall turbulence structure with wall blowing and suction // AIAA Journal. 2002. Vol. 40. No. 8. P. 1529-1535.
12. Kourta A. Instability of channel flow with fluid injection and parietal vortex shedding // Computers and Fluids. 2004. Vol. 33. P. 155-178.
13. Никитин H.B., Павельев А.А. Турбулентные течения в канале с проницаемыми стенками. Результаты прямого численного моделирования и трехпараметрической модели // Известия РАН. МЖГ. 1998. № 6. С. 18-26.
14. Vuillot F., Lupoglazoff N. Combustion and turbulent flow effects in 2D unsteady Navier-Stokes simulations of oscillatory solid rocket motors // AIAA Paper. No. 96-0884. 15 p.
15. Horiuti K. Backward scatter of subgrid-scale energy in wall-bounded and free shear turbulence // Journal of Physical Society of Japan. 1997. Vol. 66. No. 1. P. 91-107.
16. Shen S., Ding F., Han J., Lin Y.-L., Arya S.P., Proctor F.H. Numerical modeling studies of wake vortices: real case simulation // AIAA Paper. No. 99-0755. 16 p.
17. Волков К.Н. Применение метода контрольного объема для решения задач механики жидкости и газа на неструктурированных сетках // Вычислительные методы и программирование. 2005. Т. 6. № 1. С. 43-60.
18. Волков К.Н. Дискретизация конвективных потоков в уравнениях Навье-Стокса на основе разностных схем высокой разрешающей способности // Вычислительные методы и программирование. 2004. Т. 5. № 1. С. 129-145.
SUMMARY. Large-eddy simulation of turbulent flows in the ducts with wall injection is developed. The flows are described on the base of filtered Navier-Stokes equations. Sub-grid scale eddy viscosity is calculated with the help of RNG model with a correction on stream line curvature. The features of initial and boundary conditions specification on permeable wall of a duct are discussed. The calculations for different velocities of fluid injection from duct walls are established and compared with the results of direct numerical simulation, solution based on the k-8 model and experimental data.
