Численный метод решения пространственного уравнения поверхности горения конденсированных систем Текст научной статьи по специальности «Математика»

УДК 519.6:531.57

ЧИСЛЕННЫЙ МЕТОД РЕШЕНИЯ ПРОСТРАНСТВЕННОГО УРАВНЕНИЯ ПОВЕРХНОСТИ ГОРЕНИЯ КОНДЕНСИРОВАННЫХ СИСТЕМ

Ю.М. МИЛЁХИН, Г. Г. МЕДВЕДЕВ, И.Г. ВОРОПАЕВА

ФЦДТ "СОЮЗ", Дзержинский, Россия E-mail: [email protected]

АННОТАЦИЯ, Излагается численный метод решения пространственного уравнения поверхности горения конденсированного вещества для неоднородного нелинейно анизотропного поля скорости горения, основанный на принципе минимального времени выгорания заряда.

ВВЕДЕНИЕ

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

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

Для анализа трудностей, возникающих при разработке численных методов, рассмотрим уравнение горящей поверхности для линейно анизотропного поля скорости горения [1]:

Представленное уравнение относится к классу гиперболических уравнений первого порядка в частных производных и является характеристическим уравнением для волнового уравнения.

= 0.

(1)

Наличие релейного переключателя $1§п(ип) в представленном уравнении исключает возможность использования разностных схем.

Кроме того, весьма затруднителен расчет многосвязных областей и учет самопересечения границ заряда при его выгорании - образование отдельных кусков горящего топлива.

Поэтому для решения пространственного уравнения поверхности горения был разработан специальный метод, основанный на использовании вариационного принципа Ферма - принципа минимального времени выгорания заряда.

Разработанный алгоритм включает в себя следующие сравнительно автономные расчетные блоки:

- блок дискретизации тела заряда (конденсированного вещества) на элементарные объемы с использованием равномерно распределенных в пространстве точек ЛПт-последовательности;

блок расчета траекторий движения точек горящей поверхности на основе вариационного принципа Ферма - принципа минимального времени выгорания заряда;

блок расчета изменения контура исходного заряда по мере его выгорания; блок регуляризации расчетной функции изменения площади горения заряда.

БЛОК ДИСКРЕТИЗАЦИИ ТЕЛА ЗАРЯДА

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

С этой целью тело заряда разбивается на элементарные объемы, центрами которых служат точки ЛПТ- последовательности [2].

Точки ЛПТ - последовательности обладают тем замечательным свойством, что произвольная ограниченная группа точек Рь Рг,..-, Ры при любом количестве точек N представляет собой равномерно распределенную последовательность точек в единичном кубе К3 с декартовыми координатами Р( (х^ , ] = ).

При увеличении N «плотность уплотнения» куба увеличивается, а равномерность распределения точек по объему куба сохраняется.

Подобный выбор элементов позволяет наиболее оптимально заполнить внутреннюю область изделия и обеспечить любую требуемую точность расчета текущей поверхности горения за счет увеличения количества точек ЛПТ- последовательности. Необходимые теоретические сведения.

Равномерно распределенные последовательности в трехмерном кубе. Обозначим через К3 единичный куб в трехмерном пространстве. К3 состоит из всех точек Р с декартовьми координатами Р(х1 , Х2, хз), которые удовлетворяют неравенствам

0<Х|<1, у = 1,Ы]. Введем также 3 - мерные параллелепипеды я с ребрами, параллельными координатным осям. Проекция такого параллелепипеда на ось Ох] представляет собой отрезок, который обозначим Ц. Очевидно, л состоит из всех точек Р = (хь хг, хз), координаты которых удовлетворяют условиям е 1-} при 1 < ] < N. Объем такого 3-мерного параллелепипеда равен произведению У3 = I /1 I • | /2 I • I /31 -

Параллелепипед як называется двоичным, если все определяющие его отрезки /кь ¿кЗ двоичные. Заметим, что сумма всех однотипных двоичных параллелепипедов 7гк составляет весь куб К*5 .

Рассмотрим последовательность точек Р1 , Р2 , Р] , ... , принадлежащих К'1 и обозначим через Бы (в) количество точек с номерами 1 < } < Ы, принадлежащих множеству й.

Последовательность точек Р1, ..., Р], ... называется равномерно распределенной в К (сокращенно р. р.) [2], если для любого п

При больших N количество точек р. р. последовательности, принадлежащих любой области в, пропорционально объему в.

В теории р.р. последовательностей доказана [2] следующая лемма: Для того чтобы {PJ } была р. р., необходимо и достаточно выполнение условия (3) для всех двоичных пк.

о

Отклонение. Рассмотрим в К сетку, состоящую из N произвольных точек Рь Р2, ..., Ры- Каждой точке Р из К поставим в соответствие параллелепипед тср с диагональю ОР (рис.1). Объем Ур этого параллелепипеда равен произведению хь хз координат точки Р.

Отклонением сетки Рь Р2, Ры называется число

Нш 8м(1с)/К = Уя.

(2)

0(Р„ Р2, Рм) = 8ир|8ы(тгр)-ЫУ]

(3)

р

Х2 ▲

Р

О

О

Рис.1. К определению объема параллелепипеда п

р

где верхняя грань берется по всем Р е К .

Теорема. Для того, чтобы {Р1 } была р. р., необходимо и достаточно, чтобы при N -> оо

О (Рь Р2, Ры) / И —> 0 . (4)

Неравномерность. Рассмотрим в К^ произвольную сетку, состоящую из N точек Р1, Р2, ... , Ры- Выберем произвольный двоичный параллелепипед пк и перенесем начало координат в центр як. Координатные плоскости разобьют пк на 23 равновеликих гипероктантов (рис.2). Обозначим через У^ совокупность всех положительных, а через

У|7 - всех отрицательных октантов.

Величина (Ук") - 8Ы (ук+) характеризует расположение точек сетки по отношению К 7СК.

Верхняя грань этой величины по всевозможным тск

5ир^(Ук)-8ы(ЧГ) (5)

к

есть целое число, характеризующее расположение точек сетки в К3. Однако принять эту величину за 3 - мерный аналог равномерности нельзя, поскольку хорошие значения

о

(5) еще не гарантируют равномерности распределения точек сетки в К .

Гарантированную оценку можно получить, если рассмотреть проекции точек РЬР2,... ,Ры на различные грани куба К3 и для каждой из них оценить свою величину вида (5).

Обозначим через Кр какую-нибудь б - мерную грань куба К3, где 1 < э < 3. Проекции точек Р1, Р2, ... , Ры на Кр образуют в - мерную сетку в Кр. Для нее можно вычислить свою величину (5).

А к

- +

+ - р

Рис.2. К обоснованию критерия неравномерности

Неравномерностью сетки Pi, Р2, ...Pn называется наибольшая верхняя грань Ф- (Р. > > -Л )= max sup | SN (v")- SN (v;) . (6)

P k

Величина фсо - целое число и отвечает неравенству

1 < Фоэ < N. (7)

Оценка фсо через D выразится

9OO(P„P2,...PN) < 4 "D (P,,P2,...PN) . (8)

Для фоо справедлива следующая теорема:

Для того, чтобы {Pi} была р . р. необходимо и достаточно, чтобы при N со

фос (Рь Р2>—> PN) / N —» 0 . (9)

Как следует из представленных определений точки ЛПт-последовательности отвечают всем вышеприведенным критериям равномерности в отличие от равномерно распределенных величин, полученных другими методами, например, использованием датчиков случайных чисел.

Для иллюстрации использования точек ЛПх-последовательности в блоке дискретизации тела заряда рассмотрим тестовый пример. Выберем заряд осесимметричной формы и заполним его объем точками ЛПт-последовательности. На рис. 3 представлено расчетное сечение заряда.

Как видно из рисунка, точки ЛПт-последовательности равномерно заполняют объем заряда. Каждая точка ЛПт-последовательности характеризует элементарный объем заряда и используется в блоке расчета траекторий движения точек горящей поверх-

Рис.З. Внутренние точки расчетного сечения

ности. Увеличением числа точек ЛПХ - последовательности дает возможность обеспечить любую требуемую точность расчета текущего объема заряда и его поверхности.

БЛОК РАСЧЕТА ТРАЕКТОРИЙ ДВИЖЕНИЯ ТОЧЕК ГОРЯЩЕЙ ПОВЕРХНОСТИ

Для расчета последовательности выгорания элементарных объемов заряда используется прямой метод решения вариационной задачи Ферма - метод Ритца, обладающий свойством консервативности, т.е. при расчете текущей геометрии заряда не происходит накопления погрешностей. Текущее значение наружного слоя заряда определяется по заданному времени его выгорания в форме неравенства <11 < А», где (11* - принятый шаг расчета по времени.

Функционал минимального времени прохождения фронтом горящей поверхности участка дуги [80, 81] для неоднородного анизотропного поля скорости горения имеет вид

и = и (х, у, т) - поле скорости горения топлива, неоднородное по объему заряда;

[х(8), у(8)? 2(8)] - параметрическое задание траектории движения произвольной точки М (х, у, т) фронта горящей поверхности из положения М(80) в положение М^). Анизотропное поле скорости горения задается выражением

где п является вектором внешней нормали к горящей поверхности, а оператор А представляет собой в общем случае нелинейную тензорную функцию.

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

(10)

Здесь: X - время прохождения участка дуги [80, Б]];

9 9 9 1/9

(ёх + с!у + вг ) - элемент длины дуги;

п = А(п),

(П)

БЛОК РАСЧЕТА ИЗМЕНЕНИЯ КОНТУРА ИСХОДНОГО ЗАРЯДА

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

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

БЛОК РЕГУЛЯРИЗАЦИИ РАСЧЕТНОЙ ФУНКЦИИ ИЗМЕНЕНИЯ ПЛОЩАДИ ГОРЕНИЯ ЗАРЯДА

Разработанный численный метод определяет текущее значение сгоревшего объема заряда. Для определения площади горящей поверхности необходимо использовать процедуру численного дифференцирования. Задача численного дифференцирования при неточно заданных исходных данных относится к классу некорректных задач. Для решения таких задач необходимо использовать методы регуляризации, которые сводятся к введению сглаживающего функционала. Известные методы регуляризации приводят к решениям в классе функций с непрерывной первой производной. Такой подход при расчете горящей поверхности может сгладить характерные пики в зависимости изменения площади горящей поверхности от времени, обусловленные в частности выгоранием отдельных кусков заряда. Для сохранения характерных особенностей функции 8(1) в разработанном алгоритме используются следующие подходы:

1) увеличение количества точек ЛПт-последовательности;

2) увеличение шага расчета йи;

3) осреднение зависимости 8(1) при различных положениях заряда в кубе с точками

ЛПх-последовательности.

На рис.5 проиллюстрирован способ регуляризации расчетной зависимости 8(1) за счет вращения тела заряда в кубе с точками ЛПх-последовательности.

Рис. 4. Схема расчета текущего положения горящей поверхности

Рис. 5. Различные положения заряда в кубе с точками ЛПт-последовательности

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

[3,4].

Изложенный алгоритм проиллюстрируем следующими модельными примерами. РЕШЕНИЕ ТЕСТОВЫХ ЗАДАЧ

Пример 1. Заряд торцевого горения. Отличается простотой формы, отсутствием выгорающих кусков и многосвязных областей. Методы расчета геометрии выгорания таких зарядов хорошо разработаны.

На рис.6 представлен исходный контур заряда и показаны характерные положения горящей поверхности по мере его выгорания.

На рис.7 приведен график зависимости площади горения от величины сгоревшего свода.

Рис, 6. Заряд торцевого горения

28000

240001

50 ёо " 70 80 9Г

Рис. 7. Зависимость S(e)

Пример 2. Осесимметричная конструкция заряда, усложненная дополнительными вырезами.

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

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

График зависимости площади горения заряда от величины сгоревшего свода приводится на рис.9.

Пример 3. Заряд с зонтичной проточкой.

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

Осевое сечение заряда и положения поверхности горения в процессе выгорания (отдельные моменты времени) показаны на рис.10.

График зависимости S(e) для данной конструкции изображен на рис. 11.

Рис. 8. Осесимметричный заряд с развитой начальной поверхностью

Рис. 9. Зависимость S(e)

8Ü00Ü

70000

60000

50000

40000

30000

20000

10000

Рис. 10. Динамика изменения поверхности горения заряда с зонтичной проточкой

Рис. 11. Зависимость S(e)

При проведении вычислительного эксперимента прямоугольник, охватывающий сечение тестового заряда, заполнялся различным количеством точек ЛПТ-последовательности. Расчетные варианты представлены в таблице.

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

№ вар-та 1 2 3 4 5

ИлПт 30 000 300 000 1 000 000 2 000 000 3 000 000

5W, % +0.0196 -0.0099 -0.0025 -0.0004 -0.0001

öS, % +0.0568 -0.0085 -0.0052 -0.0030 -0.0006

В таблице приняты следующие обозначения:

Nnrix - количество точек Л ^-последовательности; 5W - погрешность вычисления объема;

5S - погрешность вычисления площади горящей поверхности. Как следует из таблицы, увеличение числа точек ЛПт-последовательности Nnm позволяет существенно повысить точность вычисления как текущего объема заряда, так и площади горящей поверхности. Так при числе точек Нлпт=30 ООО погрешность вычисления объема составила 5W « 0.02%, а площади горящей поверхности » 0.06%. При числе точек Нлпт=3 ООО ООО погрешности вычисления объема и горящей поверхности не превысили, соответственно, 5W » 0.0001% и 5S « 0.0006%.

Таким образом, представленный метод дает возможность рассчитывать с любой наперед заданной точностью текущие значения объема заряда W(t) и поверхности горения S(t). Более того, данный метод позволяет рассчитывать текущую геометрию заряда при наличии анизотропного неоднородного поля скорости горения.

СПИСОК ЛИТЕРАТУРЫ

1. Милёхин Ю.М., Медведев Г.Г., Воропаева И.Г. Тензорное представление скорости горения конденсированных систем. Вывод и анализ пространственного уравнения поверхности горения // Химическая физика и мезоскопия, 2006. Т.8. №1. - С. 5-20.

2. Соболь И.М., Статников Р.Б. Выбор оптимальных параметров в задачах со многими критериями. - М.: Наука, 1981. - 110 с.

3. Соркин P.E. Теория внутрикамерных процессов в ракетных системах на твердом топливе: внутренняя баллистика. - М.: Наука, 1983. - 288 с.

4. Липанов A.M., Алиев A.B. Проектирование РДТТ. - М.: Машиностроение, 1995. -450 с.

SUMMARY. Numerical method of solution of spatial equation of condensed systems burning surface for anisotropic non uniform burning rate field is stated. Method bases on principle of minimum time of condensed system burnout.