CC BY
17
УДК 519.63+65:533+539.3
ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ВЗАИМОДЕЙСТВИЯ СВЕРХЗВУКОВОГО ПОТОКА ГАЗА С ДЕФОРМИРУЕМЫМ СОПЛОМ
КУЗЬМИН ИМ., ТОНКОВ Л.Е., ЧЕРНОВА А. А.
Институт механики Уральского отделения РАН, 426067, г. Ижевск, ул. Т. Барамзиной, 34
АННОТАЦИЯ. Рассматриваются численные методы, алгоритмы и схемы, применяемые для решения задачи взаимодействия образованной в момент запуска ракетного двигателя ударной волны и упругого деформируемого корпуса сопла. В статье сопоставляются результаты, полученные с использованием различных численных схем. Проводится анализ различных факторов, влияющих на картину течения и деформирования корпуса сопла.
КЛЮЧЕВЫЕ СЛОВА: численное моделирование, сопряженная задача, газовая динамика, упругое деформирование, ударная волна, сопло.
ВВЕДЕНИЕ
Сопло в силу своего конструктивного предназначения является одним из наиболее нагруженных узлов ракетных двигателей, что связано, прежде всего, с высокими температурами, уносом материала, эрозионным воздействием, напряжениями в околокритических областях. Кроме того, на начальном этапе работы большинству конструкций ракетных двигателей первой ступени присуща проблема нарушения осесимметричности течения и возникновения боковых нагрузок [1], что обусловлено, прежде всего, высокой геометрической степенью расширения сопла.
При рассмотрении процессов, протекающих в начальный момент работы сверхзвукового сопла ракетного двигателя можно выделить несколько характерных особенностей, присущих достаточно широкому классу конструкций. Прежде всего, это формирование системы скачков уплотнения в сверхзвуковой части сопла, обусловленное перерасширенным режимом течения. При значительной величине градиента давления вблизи стенки сопла возникает характерная Л.-конфигурация скачков уплотнения и отрыв потока. Неизбежно существующая в любой реальной газодинамической системе неоднородность параметров потока приводит к тому, что на некоторых режимах реализуется особый тип отрыва, имеющий существенно трехмерный и нестационарный характер, что и приводит к возникновению динамических высокочастотных боковых нагрузок при возрастании давления в камере сгорания на этапе запуска двигателя. Данные боковые нагрузки, в свою очередь, приводят не только к колебаниям конструкции, но и к изменению формы сопла, нарушению его осесимметричности.
Общая длительность процесса запуска сопла на порядок больше по отношению к промежутку времени, в течение которого развивается и достигает максимума боковая составляющая тяги. В работе [2] запуск смоделирован полностью с учетом процесса подачи компонентов топлива в камеру и кинетики горения топлива. Здесь авторами для сокращения объема вычислений сделана попытка смоделировать процесс распространения скачка уплотнения и его взаимодействия с корпусом сопла при помощи задания разрывных начальных условий. При этом поверхность разрыва совпадает с критическим сечением сопла, значения параметров приведены в таблице.
Существующие подходы к решению комплексных сопряженных задач и особенности их реализации описаны в [3, 4]. С целью учета специфики подзадач газодинамики и механики деформирования твердых тел используется последовательный метод решения. Рассмотренные в статье алгоритмы и численные схемы верифицированы [5] на примере одной из немногих известных экспериментальных работ в данной области.
ПОСТАНОВКА ЗАДАЧИ
Рассмотрим нестационарную задачу о взаимодействии сверхзвукового потока газа с корпусом сопла [6]. Схема сопла с основными размерами приведена на рис. 1.
а) б)
Рис. 1. Расчетная область: а) схема сопла; б) расположение контрольных точек
Свойства рабочего тела
Таблица
Параметры р, МПа Т, К р , кг/м3
В сопле (X > -0,2 м) 0,101253 288 1,225
В камере сгорания (X £ -0,2 м) 26 299 306,25
В рассматриваемом диапазоне температур и плотностей целесообразно использовать уравнения состояния, полученные для реальных газов, но в данной работе, как и в [6], сделано допущение о свойствах рабочего тела - воздух считается совершенным газом с показателем адиабаты у = 1,41 и газовой постоянной Я = 287 Дж/(кг-К), учет вязкости не
производится. Корпус сопла, изготовленный из стали плотностью р = 8400 кг/м3 с модулем Юнга Е = 20,7 ГПа , полагался изотропным и однородным.
В задачах взаимодействия газа с деформируемым телом на контактных поверхностях необходимо сформулировать дополнительные кинематические и динамические граничные условия. Полная формулировка условий совместимости приведена в [3, 4]. Рассматриваемое течение газа с учетом сделанных допущений может быть описано следующей системой уравнений сохранения
Эр / д( + У-(р я и я) = 0,
ф/ъ + ,и,и,) = -У^-у.а,,
дрк8 / у-Ре и ,) = - Р8 (У - и, )-У - ц-У-(а я-уи ,)
(1)
1
2
где рг - плотность газа, рг - давление, иг - вектор скорости, Ег = —Ц, + еуТ - полная энергия, а = - /л- (уи (Уи ) )+ — /л- (Уи ) -1 - тензор вязких напряжений, Т - темпера-
тура, я - тепловой поток. Система (1) замыкается уравнением состояния идеального газа
Р,= Р
Расчетная область представляет собой внутреннее пространство сопла (рис. 1, а). Граничные условия определяются обычным образом для непроницаемой твердой теплоизолированной стенки: и = 0 на Г , Эц / Эп = 0 на Г. На границе Г скорость газа совпадает со скоростью движения границы и = иГ . Начальные данные соответствуют
условиям запуска двигателя, принятым в [6].
Динамика напряженно-деформированного состояния корпуса сопла описывалась системой уравнений: 6л\аs+ рЬ = р^и^, где ав= ,Я1ге(и^ )1 + 2^е(и^) - тензор напряжений
Коши; е(и^) = 1/2(Уи У ит) - тензор малых деформаций; ря - плотность тела. Более подробные сведения о постановках физических подзадач и методах их решения можно найти в [7].
МАТЕМАТИЧЕСКИЕ МОДЕЛИ И АЛГОРИТМЫ РАСЧЕТА
Рассматриваются разделенные подходы для учета взаимодействия потока газа и деформируемого корпуса сопла, использующие различные методы дискретизации уравнений газодинамики (метод конечных объемов, метод Эйлера со стабилизацией по методу Петрова-Галеркина и методом сквозного счета) и механики деформируемого твердого тела (методом конечных элементов) по пространству.
В работе [6] задачи газовой динамики и деформирования решались последовательно, с помощью исследовательского ПО, основанного на решении уравнений Эйлера методом сквозного счета со стабилизацией по противопоточному методу Петрова-Галеркина (8ИРО), и модуля программы РЕТ8е-БЕМ соответственно, использовались согласованные расчетные сетки (рис. 2, а), а двухсторонний обмен данными на границе сопряжения был реализован в рамках постановки слабого связывания.
В данной работе решение задачи деформирования реализуется методом конечных элементов в программном модуле БЕЗШёю [7], решение газодинамической задачи основано на решении уравнений Навье-Стокса методом конечных объемов с применением ТУБ-схемы второго порядка точности в программе ОрепБОАМ [4], использовались несогласованные на границе раздела сетки (рис. 2, б), двухсторонний обмен данными обеспечивался исследовательскими алгоритмами в рамках постановки слабого связывания.
Рис. 2. Расчетные сетки: а) согласованные [6]; б) несогласованные
Согласованные расчетные сетки [6] содержат 334700 тетраэдров для разрешения области газа и 59600 призм - для стенок сопла. Используемые несогласованные неструктурированные расчетные сетки содержат 443573 тетраэдров для разрешения области газа и 35483 тетраэдров - для стенок сопла.
При решении задачи интерполяции значений нагрузок и перемещений на границе сопряжения применялись методы на основе использования функций радиального базиса [8]. Основной идеей рассматриваемого подхода является поиск необходимого интерполянта сЖх^) в виде линейной комбинации, образованной функциями радиального базиса.
В общем случае, для задачи интерполяции данных и деформирования сетки ОР, интерполянт запишется следующим образом:
С°] (х) = ]] 1ф( X - х5 )+ х), ( (X) = {и, рп}
3=1 1
здесь х
5 <
координаты узлов, в которых заданы величины на сетке для задачи
деформирования, коэффициенты полинома V и 1. определяются из условия ( (х5 ) = ^ и
требования ]j¡S1Яjb\xS. )= 0 для всех полиномов Ь степени меньше или равной степени
полинома V. Введение дополнительного полинома V не является необходимым и применяется только в случае использования радиальных базисных функций определенного вида.
Пусть для определенности ии5 - известные перемещения узлов границы Г5. Тогда система уравнений в матричной форме для неизвестных коэффициентов 1 и л :
Ф в*' и 5
_ вТ 0 _ 0 _
где Ф = Ф\
с 3 -й строкой вида
(3)
) - матрица размерности N х N, В5 - матрица размерности N х 4 \хв у8 . Новое положение узлов в области ОР найдем как
1
=[Ф Е8вЕ ]
л
(4)
где ЖР - приращения координат узлов сетки ОР. Алгоритм может использовать интерполирующие функции радиального базиса (РБФ), то есть функции, значение которых зависит только от нормы аргумента и имеющие различные классы гладкости, что в свою очередь позволяет управлять в широких пределах свойствами интерполянта (2). Полностью аналогичным способом осуществляется интерполяция значений давления (нагрузки).
РЕЗУЛЬТАТЫ ЧИСЛЕННОГО МОДЕЛИРОВАНИЯ
Как показано в [3] при использовании достаточно подробной вычислительной сетки, аппроксимаций пространственных производных второго порядка даже в ламинарном приближении осевая симметрия течения нарушается.
Взаимодействие системы скачков уплотнения и отрывных зон приводит к существенным колебаниям давления, что, в свою очередь, наиболее наглядно проявляется в картине деформаций корпуса сопла (рис. 3). Однако, основной целью данной работы, помимо предварительного тестирования реализованных алгоритмов связывания, является оценка влияния на решение газодинамической задачи особенностей моделируемого процесса обусловленных деформацией области интегрирования.
Рис. 3. Фронт ударной волны и деформация корпуса сопла
По результатам расчетов в контрольных точках 1-5 (рис. 1) построены зависимости изменения давления от времени, как без учета газоупругого взаимодействия (рис. 4), и с учетом деформирования сопла (рис. 5).
Здесь необходимо отметить, что в работе [6] не указано для каких именно точек на стенке сопла построены рассматриваемые зависимости. Точки 1-5 на рис. 1, что очевидно, определяют окружности на стенке сопла. Поэтому авторами при построении графиков давления использовалось то меридиональное сечение, в котором на окружности, соответствующей контрольной точке 1, достигался абсолютный максимум давления.
Рис. 4. Графики зависимости давления от времени без учета газоупругого взаимодействия: точка 1: 1 - [6], 2 - расчет; точка 2: 3 - [6], 4 - расчет; точка 3: 5 - [6], 6 - расчет; точка 4: 7 - [6], 8 - расчет; точка 5: 9 - [6], 10 - расчет
Из рис. 4 видно, что рассчитанное время достижения максимума давления хорошо согласуется с результатами [6], однако, полученные в рамках данной работы зависимости характеризуются как большей интенсивностью нагрузок (локальные максимумы давлений в 1,5 раза превышают аналогичные величины, полученные в [6]), так и более резким ростом. Повышение интенсивности ударной волны и, как следствие, повышение интенсивности нагрузок на стенки сопла, может быть связано как со способом задания начальных условий, так и с размерами и качеством используемых расчетных сеток. Отметим общее снижение интенсивности ударной волны при ее движении к срезу сопла. Снижение значений нагрузок, полученных с использованием второго подхода, имеет, как и в [6], практически линейный характер, но характеризуется большей скоростью, вследствие чего в точке 5 (рис. 1, б) локальные максимумы давлений, полученные с использованием обоих подходов, совпадают (рис. 4, кривые 9,10).
Как видно из рис. 5, учет газоупругого взаимодействия приводит к незначительному уменьшению максимальных значений давления, относительно рассчитанных без учета деформирования сопла. Также отмечается появление второго локального максимума (рис. 5, кривые 2, 4, 6, 8, 10), что, вероятнее всего, непосредственно связано с переходом свободного отрыва в присоединенный в результате взаимодействия скачка уплотнения с деформированной стенкой сопла. Так же как и в предыдущем расчетном случае отмечается превышение интенсивности скачка давления на стенке, относительно данных [6], но величины максимумов давлений отличаются не более чем в 1,2 раза для первой и второй контрольных точек и соответствуют максимальным значениям давлений, полученным в [6], для 3 и 4 точек. В контрольной точке 5 (рис. 5, кривые 9, 10) локальные максимумы кривой давления 10, как и средние значения давления, меньше соответствующих значений [6] на 8 %.
800 000
600 ООО
го
с..400 ООО
CL
200 ООО
О
О 0,002 0,004 0,006 0,008 0,01
Рис. 5. Графики зависимости давления от времени с учетом деформирования сопла: точка 1: 1 - [6], 2 - расчет; точка 2: 3 - [6], 4 - расчет; точка 3: 5 - [6], 6 - расчет; точка 4: 7 - [6], 8- расчет; точка 5: 9 - [6], 10 - расчет
Из рис. 5 видно, что полученные в рамках данной работы кривые давления схожи с кривыми [6]. Однако, учет деформирования сопла, предложенный в данной работе, позволяет оценить изменение характера нагрузок при деформации корпуса сопла.
ЗАКЛЮЧЕНИЕ
В результате проведенного численного эксперимента показано, что учет существенно нестационарных деформаций корпуса, несмотря на малость их абсолютных величин, оказывает заметное и воспроизводимое от расчета к расчету влияние на газодинамическую картину течения, что обнаруживается хотя бы по изменению кривых давления на стенке сопла. Реализованные методы и алгоритмы решения сопряженных задач позволяют получить качественно согласующиеся результаты, однако, говорить о точности вычислений как газодинамических параметров, так и деформаций корпуса сопла на данный момент не представляется возможным. Несмотря на выполненное ранее авторами решение тестовых задач для верификации моделей, методов и алгоритмов как слабого, так и сильного связывания, остается актуальным сравнение получаемых результатов с экспериментальными данными и результатами аналогичных вычислительных экспериментов именно для рассматриваемой задачи о воздействии на сверхзвуковое сопло нестационарного несимметричного потока газа.
Работа выполнена при поддержке грантов РФФИ № 14-01-00055-a, 14-08-00064-а. СПИСОК ЛИТЕРАТУРЫ
1. Tomita T., Sakamoto H., Onodera T., Sasaki M., Takahashi M., Tamura H. Experimental Evaluation of Side Load Characteristics on TP, CTP, and TO Nozzles // AIAA Paper. 2004. № 04-3678. P. 1-9.
2. Wang T.S. Transient three-dimensional startup side load analysis of a regeneratively cooled nozzle // Shock Waves. 2009. V. 19. P. 251-264.
3. Тонков Л.Е., Чернова А.А. Моделирование возникновения боковой составляющей тяги сверхзвукового сопла // Десятая Междунар. конф. «Сеточные методы для краевых задач и приложения». Казань : Изд-во Казанского ун-та, 2014. С. 598-602.
4. Копысов С.П., "Гонков Л.Е., Чернова A.A. Сравнение слабо- и жесткосвязных постановок при моделировании взаимодействия сверхзвукового потока и деформируемой панели // Химическая физика и мезоскопия. 20l2. T. l4, № 4. С. 545-552.
5. Копысов С.П., 1онков Л.Е., Чернова A.A. Двухстороннее связывание при моделировании взаимодействия сверхзвукового потока и деформируемой пластины. Сравнение численных схем и результатов эксперимента // Вычислительная механика сплошных сред. 20l3. T. б, № l. С. 78-8S.
6. Luciano Garellia, Rodrigo R. Paza, Mario A. Storti Fluid-Structure interaction study of the start-up of a rocket engine nozzle // Coumputers and Fluids. 2009. V. 39. P. l208-l2l8.
7. Караваев A.C, Копысов С.П., Кузьмин И.М., Недожогин H.С., Новиков A.K, Рычков В.Н., Сагдеева ЮА., Сармакеева A.C, 1онков Л.Е., Чернова A.A. Конечно-элементные вычисления в программном комплексе FEStudio на кластерах гибридной архитектуры // В сб. «Механика и физико-химия гетерогенных сред, наносистем и новых материалов». Ижевск : Изд-во ИМ УрО PA^ 20lS. С. l0S-l29.
В. Караваев A. С., Копысов С. П., Кузьмин И. М. Метод консервативной интерполяции на нестыкующихся поверностных сетках // Вестник Удмуртского университета. Математика. Механика. Компьютерные науки. 20l4. № 4. С. 64-7S.
NUMERICAL MODELING OF INTERACTION OF A SUPERSONIC GAS FLOW WITH THE NOZZLE CASE
Kuzmin I.M., Tonkov L.E., Chernova A.A.
Institute of Mechanics, Ural Branch of the Russian Academy of Sciences, Izhevsk, Russia
SUMMARY. The numerical methods, algorithms and schemes used for solution of a problem of interaction of the shock wave formed at launch of a rocket engine and the elastic deformable nozzle are considered. The article given a discussion of the revealed delay of a shock wave, in relation to the analytically calculated. The numerical results received with use of various numerical schemes are compared. The analysis of the various factors affecting the pattern of the flow and deformation of the case of a nozzle is carried out.
KEYWORDS: the mathematical modeling, the interfaced task, gas dynamics, elastic deformation, shock wave, nozzle.
Кузьмин Игорь Михайлович, младший научный сотрудник ИМ УрО РАН, e-mail: imkuzmin@gmail. com
Тонкое Леонид Евгеньевич, кандидат физико-математических наук, доцент, старший научный сотрудник ИМ УрО РАН, e-mail: tnk@udman. ru
Чернова Алена Алексеевна, кандидат технических наук, научный сотрудник ИМ УрО РАН, e-mail: alicaaa@,gmail. com
