Non-stationary behavior of a solid propellant charge for nozzleless solid rocket motors under gas-dynamic load.pdf Одним из основных требований, предъявляемых к современным интегральным твердотопливным ракетно-прямоточным двигателям (РПД) является отсутствие отбрасываемых элементов конструкции при переходе с режима разгона на маршевый режим полета ракеты. В связи с этим, в качестве стартово-разгонных ступеней нашли применение бессопловые ракетные двигатели (БРД), в конструкции которых нет традиционного соплового блока и соответственно нет необходимости его сбрасывания. Скрепленный с корпусом заряд твердого топлива бессопловой стартовой ступени размещается в камере сгорания прямоточного контура маршевой ступени интегрального РПД. При этом оптимизируются массово-габаритные характеристики и упрощается конструкция ДУ. Однако существуют проблемы, связанные с большими перепадами давления по каналу заряда и уменьшением удельного импульса по сравнению с сопловыми двигателями. Для выполнения тягово-энергетических требований, предъявляемых к двигательной установке (ДУ), возникает необходимость профилирования канала заряда таким образом, чтобы выходная его часть представляла собой имитатор сопла, разгоняющий поток продуктов сгорания твердого топлива (ТТ) до необходимых параметров. В подобных конструкциях критические условия течения газовой смеси (скорость звука) достигаются в сечении канала, предшествующем расширяющейся «сопловой» части. При этом на малом участке длины канала реализуется сверхкритический перепад давления (порядка 0.5Pmax), что ведет к большим напряжениям внутри заряда и деформациям канала, способным разрушить заряд в начале квазистационарного участка работы ДУ. Широкий температурный диапазон эксплуатации РПД (-50, +70) °С предъявляет противоречивые требования к физико-механическим свойствам твердого топлива. Уменьшение модуля упругости ТТ при высоких температурах способствует снижению максимальных напряжений в заряде, но ведет росту опасных деформаций, и наоборот. Одной из задач, возникающих при оптимизации формы заряда и прогноза внутрибаллистических характеристик БРД, является определение максимально возможных перемещений поверхности заряда ТТ после приложения газодинамической нагрузки. С одной стороны, максимальные осевые перемещения определяют размеры торцевых зазоров между корпусом БРД и поверхностью заряда ТТ, с другой стороны, радиальные перемещения поверхности заряда изменяют как величину площади критического сечения сопла, так и площадь поверхности горения, в результате чего изменяется давление в камере сгорания (КС), которое в свою очередь изменяет нагрузку на заряд ТТ. В этом смысле задача определения величины перемещения поверхности ТТ является сопряженной задачей. Решение сопряженной задачи для классического РДТТ с традиционным сопловым блоком и прочноскрепленным зарядом типа «зонтик» в квазистационарной постановке рассматривалось в работе [1]. В работе [2] для модельного двигателя с одноканальным прочноскрепленным зарядом определялись частоты колебаний корпуса, сопла и газовой полости в продольном и радиальном направлениях на основе упруго-массовой модели РДТТ. Для рассмотрения резонансных явлений в РДТТ предполагалось, что возникающие деформации в заряде ТТ не влияют на положение поверхности горения. Упругие свойства современных твердых топлив сильно зависят от температуры термостатирования. Увеличение температуры, как правило, ведет к снижению модуля упругости ТТ, в результате чего оно становится более податливым внешним нагрузкам. Поэтому для разных температур термостатирования заряда ТТ точки его поверхности будут вести себя по-разному под действием одной и той же приложенной нагрузки. Современные ТТ (композиция на основе HTPB и HMX, HTPB и RDX) представляют собой гиперупругий материал, который характеризуется относительно низким модулем упругости и высоким объемным модулем [3]. У таких материалов функция плотности упругой энергии деформации является скалярной функцией тензоров деформации, чьи производные по компонентам деформации определяют соответствующие компоненты напряжения. Поведение таких материалов для малых и умеренных деформаций можно описывать моделью Neo-Hookean (нео-Гука), моделью Огдена или трехпараметрической моделью Муни - Ривлина, а для больших деформаций - использовать модель Arruda - Boyce (Арруда и Бойса) [4]. Целью настоящей работы является исследование нестационарного поведения осесимметричного заряда ТТ бессоплового РДТТ под действием сил давления со стороны продуктов сгорания и оценка максимальных перемещений точек поверхности ТТ на основе численного решения уравнений механики деформируемого твердого тела и газовой динамики продуктов сгорания. Физико-математическая постановка задачи Рассматривается модельный осесимметричный бессопловой РДТТ, состоящий из стального корпуса 1, термоизоляционного вкладыша 2, резиновой прослойки 3, заряда ТТ 4. Область 5 соответствует КС, рис. 1. С поверхности твердого топлива продукты сгорания поступают в канал и вытекают через выходное отверстие наружу. Под действием давления продуктов сгорания заряд твердого топлива деформируется, что приводит как к изменению И.Г. Воропаева, А.А. Козулин, Л.Л. Миньков, Э.Р. Шрагер 50 площади поверхности горения и тем самым к изменению газоприхода, так и к изменению площади критического сечения канала, что в конечном итоге изменяет давление внутри канала. Рис. 1. Разрез бессоплового РДТТ: 1 - корпус, 2 - вкладыш, 3 - прослойка, 4 - заряд, 5 - камера сгорания Fig. 1. Section of a nozzleless solid rocket motor (SRM): 1 - body, 2 - insert, 3 - interlayer, 4 - charge, and 5 - combustion chamber Для исследования протекающих процессов в БРД выделим область, представляющую собой трехмерный цилиндрический сектор, который состоит из подобластей для заряда, вкладыша, прослойки, корпуса и канала, рис.1 Внутрикамерные процессы в ДУ описываются нестационарными уравнениями Навье - Стокса для вязких сжимаемых продуктов сгорания ТТ, уравнениями напряженно-деформированного состояния для упругих и гиперупругих тел и условиями, сопрягающими решения задач на границах разнородных материалов. Материалы, из которых изготовлены прослойка и заряд, предполагаются гиперупругими, а материалы, из которых выполнены корпус и вкладыш, - упругими. Система уравнений для НДС Уравнение сохранения массы d pd dt + Pd Vk U = 0. Уравнение движения частицы среды d d2uf dt2 = Vk, (2) где Uj - перемещение частицы среды. Тензор конечных деформаций 4 = 2 (+VA+VU V A). Связь между напряжениями и деформациями: - для упругого тела: a d = 2ц 4 + v 1 - 2v d 8y Skk - для гиперупругого тела: ad = 1 FXFfl- (3) (4) (5) Нестационарное поведение заряда ТТ бессоплового РДТТ 51 где Fj =dXjldXj - компоненты тензора градиентов деформации, SiJ- - компоненты напряжения, выражающиеся через энергию деформаций W в виде S! = 2 W (6) С Здесь Cij = FkiFkj - правый тензор деформаций Коши - Грина; J = det [ F ] • Согласно модели Neo-Hookean, энергия деформации имеет вид где Wd = Ц +Л ( j -1)2. 2 - ' dd I1 = ijJ2/3 , I1 =Xj2 + X2 + X2, J = XjX2X3; () i1 - первый инвариант матрицы Cy; Xj2, X2, X2 - собственные значения матрицы С; цd = Ed 2 (1 + ѵ d) начальный модуль сдвига материала; dd = 2/Kd параметр несжимаемости материала; K d = - начальный объемный модуль Ed 3 (1 - 2ѵ d) упругости; E - модуль Юнга; vd - коэффициент Пуассона. Индекс d соответствует номеру материала: 1 - корпус, 2 - вкладыш, 3 - прослойка, 4 - топливо. Граничные условия Граница раздела поверхность твердого тела - газ: ^nn = - Р + Tnn - равенство нормальных напряжений в топливе и газе; CTd = Т • CTd = т ПТі ПХ^ nT2 nT2 - равенство касательных напряжения в топливе и газе. Граница контакта между разнородными материалами: _d = _d+1. _d = _d+1. _d = _d+1 °nn~ ^nn ; an*1 ~ ^яі1 ; °ni2 ~ ^n*2 равенство нормальных и касательных напряжений; „,d _ „,d+1 . „.d _„,d+1. „.d _„,d+1 Unn = Unn ; un%1 = UnX1 ; Un%2 = Un%2 - равенство нормальных и касательных перемещений; d = 1, 2, 3. Плоскость симметрии: = о равенство нулю нормальных перемещений. Условия закрепления: uf = о равенство нулю перемещений во всех направлениях. 52 И.Г. Воропаева, А.А. Козулин, Л.Л. Миньков, Э.Р. Шрагер Уравнения течения вязкого сжимаемого газа Уравнение неразрывности ЯГ + = РЛі РѴKer5(xk ^ Rk ) • dt dXfo Уравнение сохранения количества движения dPu: + dPu:uk =-dp I dTik dt dxk dx. dxk () (9) Уравнение сохранения энергии dpe dpukh д (, dT -k-+ k =-1 X-+ Tiku. dt dxk dxk ^ dxk + CpTpPTU01 PVKerS(xk ; Rk ), (10) где Rk , 1 p V2 функции, описывающие поверхность горения; e =-- +---полная k -1 р 2 внутренняя энергия; h = e + - - полная энтальпия; V2 =У uj - квадрат вектора Р k=i скорости; Ker -коэффициент, учитывающий эффект эрозии, задается в виде [5] K =| 11_при 1 < hr еГ W1 + ker (I - Ikr ) 1 при 1 * hr • 1 где I = - - параметр Вилюнова; Ikr - пороговое значение параметра Ptm01 p Вилюнова; Xr - коэффициент сопротивления; U - скорость продуктов сгорания, осредненная по поперечному сечению канала; % - тензор вязких напряжений; (du. duk ---S ; t - время; x - пространственная координата; p - давлеk dx. з dx, k I J J ние; u - скорость; k - отношение удельной теплоемкости при постоянном давлении к удельной теплоемкости при постоянном объеме C-ICv ; R - газовая по ' + ' dxk dx. ik стоянная. Коэффициенты Ikr и ker определялись полуэмпирическим способом из работы [6] и принимались равными 5.6 и 0.08 соответственно. Уравнение состояния: p = pRT (11) Граничные условия Недеформируемая твердая стенка: uk = 0. Деформируемая твердая стенка, топливо: dt Нестационарное поведение заряда ТТ бессоплового РДТТ 53 Выход из канала: Р = Рн - при дозвуковом истечении давление равно давлению окружающей среды. Плоскость симметрии: ип = 0, ^ = 0, Ф = (р, от1,от2,р} . дп В данной постановке горящая поверхность заменена твердой стенкой, а приход массы и энергии газов моделируется через источниковые члены в соответствующих уравнениях (8) и (10) [7]. Метод решения Твердотельные подобласти разбиваются на конечные элементы, а подобласть канала - на конечные объемы. В окружном направлении задавался один конечный элемент (конечный объем) в соответствии с приближением осесимметричности. Для численного решения уравнений (1) - (5) используется метод конечных элементов построенный на линейных функциях формы для конечных элементов, представляющих собой шестигранники. Для численного решения системы уравнений (8) - (11) используется метод Па-танкара [8] и противопоточная разностная схема второго порядка аппроксимации для конвективных слагаемых. Дельта-функция Дирака 5( xk; Rk) в уравнениях (8) и (10) аппроксимируется отношением площади соприкасающейся с поверхностью горения грани ячейки к объему этой ячейки. Результаты решения задачи (1) - (11) были получены с помощью пакета Ansys (Fluent и Transient Structural). При решении сопряженной задачи давление со стороны газа (нагрузка) прикладывалась не сразу, а разбивалась на 2 этапа: на первой итерации к поверхности заряда прикладывалась половина избыточной нагрузки, а на втором этапе - полная избыточная нагрузка. Количество итераций на каждом этапе сопряжения давление - перемещение составляло не более 5 раз. Обсуждение результатов В экспериментальном исследовании [9] показано, что в результате компенсационного эффекта факторов эрозии и зависимости от давления скорости горения, разница между скоростями горения в разных сечениях канала не превышает 10%. Распределения относительных значений скорости горения U = ub/umax и давления P = P/Pmax вдоль оси симметрии заряда ТТ (x = x/xmax ), полученные в результате предварительных расчетов газодинамических параметров, без учета влияния деформации заряда ТТ, для разных моделей скорости горения, представлены на рис. 2 и 3. Среднеинтегральные значения скорости горения топлива на дозвуковом участке канала заряда определяют газоприход и соответственно средний уровень давления в камере (рис. 3). Как следует из рис. 2 (кривая 1), расчетные данные хорошо согласуются с [9] на участке поверхности заряда, ограничивающем дозвуковую часть газодинамического тракта канала. Участок резкого спада давления (рис. 3) соответствует переходу от дозвуковой части газодинамического тракта к сверхзвуковой. 54 И.Г. Воропаева, А.А. Козулин, Л.Л. Миньков, Э.Р. Шрагер 0 0.2 0.4 0.6 0.8 x Рис. 2. Распределение скорости горения по поверхности ТТ: 1 - Ub = U01P VKer ; 2 - ub = U01P V Fig. 2. Burning rate distribution over the propelant surface: 1 Ub - U01P K-er and 2 Ub - U01P 0 0.2 0.4 0.6 0.8 x Рис. 3. Распределение давления по каналу для различных моделей скорости горения: 1 - Ub = U01pVKer ; 2 - Ub = U01pV ; 3 - Ub = U01P0V Fig. 3. Channel pressure distribution for different burning rate models: 1 Ub - U01-P Ker; 2 Ub - U01P , and 3 Ub - U01-P Поскольку точное определение параметров эрозионного горения топлива на этапе проектирования ДУ проблематично, то, как показывает рис. 3, для предварительных оценок ВБХ на основном участке канала целесообразно использовать закон скорости горения Ub - U01P0V (кривая 3), где Р0 - давление у переднего дна КС. При этом прогноз давления будет примерно на 5 % превышать значения, рассчитанные по модели с учетом эрозионной составляющей скорости горения. Для решения сопряженной задачи были приняты следующие значения параметров. Механические свойства материалов приведены в табл. 1. Характеристики твердого топлива приведены в табл. 2, а теплофизические параметры продуктов сгорания ТТ - в табл. 3. Термодинамическая оценка влияния начальной температуры в диапазоне от -50 °С до + 70 °С на теплофизические параметры продуктов сгорания показывает, что ее учет приводит к отклонениям в параметрах, не превышающим 1 %. Нестационарное поведение заряда ТТ бессоплового РДТТ 55 Т аблица 1 Механические свойства материалов БРД Топливо Резина Вкладыш Сталь Плотность, кг/м3 1710 1244 1479 7850 Модуль упругости, Па 5-106 при t = 70°C 107, при t = -50°C 9.807-106 1.442-1010 2-1011 Коэффициент Пуассона 0.4999 0.499 0.2 0.3 Т аблица 2 Характеристики твердого топлива Вариант 1 Вариант 2 Начальная температура ТТ, °С -50 +70 Скорость горения при 1 атм, иП1, м/с 4.48-10-3 6.4-10-3 Показатель степени в законе горения, ѵ 0.428 0.41 Температура горения, K 2960 2960 Т аблица 3 Теплофизические параметры продуктов сгорания ТТ Теплоемкость, Дж/(кг-К) 1800 Молярная масса, кг/моль 25.4 Динамическая вязкость, Па-с 1.79-10-5 Теплопроводность, Вт/(м-К) 0.0242 При численном решении сопряженной задачи давление в камере сгорания задавалось равным атмосферному. Задача о зажигании поверхности ТТ не рассматривалась, и предполагалось, что вся поверхность ТТ начинает гореть одновременно. На рис. 4 показаны зависимости изменения максимального давления в камере сгорания от времени для скоростей горения ТТ, соответствующих разным начальным температурам заряда ТТ (см. табл. 2). • •••** • 1 1 ■ ■ ■ . ■ввнВввввв • ' ■ ■ ■ ■ ■ • и ■ .1 *2 • 1- ■ 0 0.005 0.01 0.015 t, c P, МПа 100 II 80 • 60 , ' 1 *2 40 " 20 Рис. 4. Зависимость максимального давления в КС от времени: 1 - вариант 1; 2 - вариант 2 Fig. 4. Dependence of the maximum pressure in a combustion chamber on time: 1 - variant 1 and 2 - variant 2 И.Г. Воропаева, А.А. Козулин, Л.Л. Миньков, Э.Р. Шрагер 56 Из рисунка видно, что квазиустановление происходит примерно через 6 мс, после чего начинают происходить колебания максимального давления, причем амплитуда колебания тем выше, чем выше скорость горения ТТ (кривая 2). Период первого колебания составляет около 5 мс. Затем частота колебания увеличивается и период колебания уменьшается до 2-3 мс, что соответствует времени прохождения звуковой волны по продуктам сгорания вдоль канала туда и обратно. Характерная скорость распространения поперечных волн по заряду ТТ, оцененная по формуле т/Ej[2(1 + v)р] , составляет 31 и 44 м/с для варианта 1 и 2 соответственно. Толщина заряда в окрестности критического сечения (наиболее узкой части канала) равна 4 см. Тогда время распространение поперечных волн в радиальном направлении (туда и обратно) находится в диапазоне от 2 до 3 мс. Для оценки амплитуды колебания давления в камере сгорания, вызванного колебаниями критического сечения сопла, воспользуемся формулой Бори. Согласно этой формуле, давление в камере сгорания РДТТ po на стационарном участке обратно пропорционально радиусу критического сечения сопла Ro в степени 2/(1 -v). Если радиус критического сечения сопла незначительно изменяется со временем, например по периодическому закону, то и давление также будет изменяться со временем: p(t) ~ R(t)-2/(1-v). Отсюда имеет место равенство p(t )/ Ро =(Rq/ R(t) )2/(1-v). Пусть изменение радиуса критического сечения сопла является осциллирующим и имеет вид R(t) = R0 + Ar • exp(iot), где a - частота колебания, а Ar - амплитуда колебания радиуса. Тогда Ml = f1 + Ar . exp (а*О"2'™ -1, Ро I Ro J где Ap(t) = p(t) - po. A 2 A Если -R

Скачать электронную версию публикации