ОГЛАВЛЕНИЕ

СПИСОК СОКРАЩЕНИЙ, УСЛОВНЫХ ОБОЗНАЧЕНИЙ И СИМВОЛОВ

ВВЕДЕНИЕ

1. ОБЩИЕ СВЕДЕНИЯ О РАКЕТЕ 3М-14

2. РАСЧЕТ ДИНАМИКИ ВЫХОДА РАКЕТЫ ИЗ ШАХТНОЙ ПУСКОВОЙ УСТАНОВКИ (аналитический расчет)

2.1 Исходные данные

2.2 Постановка задачи

2.3 Проектирование ШПУ

2.4 Расчет параметров выхода ракеты при заданной схеме

Выводы

Анализ полученных результатов и заключение

3. МЕТОДИКА ЧИСЛЕННОГО ЭКСПЕРИМЕНТА

3.1 Постановка задачи

3.2 Общие сведения о проведении численного эксперимента газовой динамики

3.3 Компьютерные пакеты для численного решения задач газовой динамики

3.4 Общие сведения о методе контрольного объема

Выводы

4. РАСЧЕТ ДИНАМИКИ ВЫХОДА РАЕТЫ ИЗ ШАХТНОЙ ПУСКОВОЙ УСТАНОВКИ (численный эксперимент)

4.1 Постановка задачи

4.2 Стратегия постановки численного эксперимета

4.3 Анализ полученных результатов

Выводы

5. ОПРЕДЕЛЕНИЕ АЭРОДИНАМИЧЕСКИХ НАГРУЗОК

Введение

5.1 Постановка задачи

5.2 Подготовка и проведение численного эксперимента

5.3 Анализ полученных результатов

5.3.1 Вертикальная составляющая потока ()

5.3.2 Поперечная составляющая потока ()

5.3.3 Сложение составляющих потока (;)

5.4 Определение моментов

Выводы

ВЫВОДЫ

ЗАКЛЮЧЕНИЕ

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

ПРИЛОЖЕНИЕ

СПИСОК СОКРАЩЕНИЙ, УСЛОВНЫХ ОБОЗНАЧЕНИЙ И СИМВОЛОВ

ШПУ – шахтная пусковая установка;

АПЛ – атомная подводная лодка;

СУ – система управления;

П.П.П. – программный прикладной пакет;

ЭВМ – электронная вычислительная машина;

ПЭВМ – персональная электронная вычислительная машина;

МКО – метод контрольного объема;

ЛА – летательный аппарат;

АДХ – аэродинамические характеристики;

АДК – аэродинамические коэффициенты;

p – давление в подракетном пространстве;

Re – число Рейнольдса;

T_к – температура в камере сгорания;

T₂ – температура в «подракетном» пространстве;

F_x , F_y , F_z – проекция силы на ось X, Y, Z, соответственно;

X – сила лобового сопротивления;

Y – подъемная сила;

Z – боковая сила;

С_x – коэффициент подъемной силы;

С_y – коэффициент нормальной силы;

х_ц.д. – положение центра давления;

С_p – коэффициент давления.

ВВЕДЕНИЕ

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

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

Перед исполнителем поставлена задача создания схемы расчета динамики выхода ракеты из шахты и определения действующих на нее нагрузок. Решение поставленной задачи ведется на примере ракеты 3М-14.

В найденных литературных источниках тема исследования динамики выхода ракеты из ШПУ недостаточно освещена, ввиду специфичности тематики.

На основе вышеизложенного материала и по мере выполнения задания, в работе сформулированы следующие задачи:

1. Решение задачи динамики выхода ракеты из пусковой установки (газодинамическая задача);

2. Определение А.Д.Х. при поперечном обтекании корпуса ракеты под действием ветровой нагрузки.

1. ОБЩИЕ СВЕДЕНИЯ О РАКЕТЕ 3М-14

Ракета 3М-14 входит в состав системы «Клаб». «Клаб» (Club) – ракетная система с крылатыми ракетами надводного, подводного, наземного и авиационного базирования.

• История создания

Ракетная система «Клаб» (Club) разработана и производится (основные элементы) ОКБ «Новатор» (г. Екатеринбург). Первый испытательный пуск противокорабельной ракеты (ПКР) системы, по данным СМИ, состоялся с атомной подводной лодки (АПЛ) на Северном флоте в марте 2000 г., второй - в июне того же года с дизельной подводной лодки (ДПЛ) проекта 877 Балтийского флота. Оба пуска были признаны успешными.

Первым основным элементом системы является универсальная ракета «Альфа», которая была продемонстрирована в 1993 г. (через 10 лет после начала ее разработки) на выставке вооружений в Абу-Даби и на международном авиакосмическом салоне МАКС-93 в г. Жуковский. В том же году она была принята на вооружение. [1]

По западной классификации ракета получила обозначение SS-N-27 Sizzler (от «sizzle» - шипящий звук, издаваемый кипящим на сковороде маслом). В России и за рубежом (по данным различных СМИ, справочников серии Jane's и пр.) она обозначалась как Klub, «Бирюза» (Biryuza) и «Альфа» (Alpha или Alfa).

• Назначение

Ракетная система «Клаб» предназначена для поражения надводных кораблей и подводных лодок противника всех типов при ведении боевых действий в условиях сильного радиоэлектронного и огневого противодействия.

• Состав

В состав ракетной системы «Клаб» (Club) входят ракетные (ударные ракетные) комплексы (РК) «Клаб-Н» (Club-N) и «Клаб-С» (Club-S), которые устанавливаются на надводных кораблях и подводных лодках соответственно в качестве ударного ракетного оружия.

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

Ракеты системы во многом унифицированы между собой, но, в зависимости от предназначения и базирования, имеют различные наименования и некоторые отличия. [1]

• Ракета 3М-14

Двухступенчатая крылатая ракета для поражения наземных (береговых) целей подводного (ЗМ-14Э) и надводного (3М-14ТЭ) базирования по внешнему виду, аэродинамической схеме, габаритным характеристикам и двигательной установке аналогична ПКР ЗМ-54Э1 и имеет сходство со стратегической КР ракетного комплекса РК-55 «Гранат» (дальность стрельбы до 3000 км). Отличается фугасной (вместо проникающей) БЧ, подрыв которой осуществляется в воздухе для нанесения максимального ущерба объекту и активной радиолокационной головкой самонаведения АРГС-14Э (ОАО «Радар ММС», г. Санкт-Петербург) с высокоэффективной системой наведения ракеты на цель на конечном участке траектории полета. По этим показателям она превосходит зарубежные аналоги, в т.ч. и американскую Tomahawk, которой может быть поставлена помеха в системе спутниковой навигации GPS. При стартовой массе 2000 кг (БЧ 450 кг) и скорости полета до 240 м/с способна поражать цели на дальности до 300 км. Впервые была показана в февраля 2004 г. на 3 Международной выставке сухопутных и военно-морских вооружений «Defexpo India» (г. Дели). При ее разработке в качестве прототипа была использована стратегическая крылатая ракета "Гранат" (код НАТО SS-N-21 Sampson), предназначенной для вооружения атомных подводных лодок проекта 971, 945, 671РТМ, 667АТ и др. [1]

Таблица 1.1. Основные характеристики 3М54Э1/ТЭ1

Характеристика	Показатель
Длина, м	6,200/8,916
Диаметр, м	0, 533/0, 645
Дальность стрельбы максимальная, км	300/275
Высота полета, м на маршевом участке на конечном участке	10-20 менее 10
Максимальная скорость, М на маршевом участке на конечном участке	0,6-0,8 0,6-0,8
Масса, кг: стартовая БЧ	1780/1505 400
Система управления и наведения	инерциальная+активная СН

• Пусковые установки

Пусковыми установками для пуска всех ракет комплекса «Клаб-С» (Club-S), независимо от их назначения, являются штатные 533-мм торпедные аппараты, а ракет комплекса «Клаб-Н» (Club-N) – подпалубные унифицированные ВПУ (типа ЗС-14Э, «КБСМ», г. Санкт Петербург) модульной конструкции с индивидуальными бронированными крышками (для строящихся) или наклонные ПУ (для модернизируемых) кораблей.

• Модификации

«Клаб-М» (Club-M) наземного базирования разработан для поражения наземных и надводных целей ракетами 3М-14КЭ и 3М-54КЭ соответственно. Дальность обнаружения целей: в пассивном (активном) режиме достигает 450 (250) км, дальность стрельбы по надводным (наземным) целям - до 220 (275) км.

В состав комплекса входят самоходные пусковые установки, ракеты для поражения наземных и надводных целей, машины связи и управления, машины технического обслуживания, комплекс наземного технологического оборудования, учебно-тренировочные средства. Пусковая установка была представлена на МАКС-2007. [1]

2. РАСЧЕТ ДИНАМИКИ ВЫХОДА РАЕТЫ ИЗ ШАХТНОЙ ПУСКОВОЙ УСТАНОВКИ (АНАЛИТИЧЕСКИЙ РАСЧЕТ)

2.1 Исходные данные

Для расчета были приняты следующие исходные данные:

Стартовый вес – 1225кг;

Длина со стартовым ускорителем – 6,4м;

Диаметр корпуса – 0,53м;

Тяга разгонного двигателя – 13800Н;

Давление в камере сгорания – 7 Па;

Диаметр критического сечения сопла – 0,039м;

Диаметр выходного сечения сопла – 0,15м;

Давление атмосферное – 99600Па.

Так как данные о топливе ракеты 3М-14 отсутствуют, то было выбрано топливо на основе перхлората аммония (68%), полиуретана (17%) и алюминия (15%).

Характеристики топлива:

Удельный импульс – 2460Н;

Газовая постоянная – 290 Дж/(кг*К);

Показатель адиабаты – 1,16;

Температура горения – 3300К. [2]

2.2 Постановка задачи

На основе схемы, изображенной на рис.2.1, разработать и определить схему расчета параметров выхода ракеты из ШПУ. Для подтверждения качества математической модели и отработки расчетного алгоритма определить один из геометрических параметров ШПУ, а именно ширину отвода газа из «подракетного» пространства и сделать выводы о реальности данного параметра.

Рисунок 2.1. Схема ШПУ с ракетой

2.3 Проектирование ШПУ

Для проектирования шахтной пусковой установки (ШПУ), требуется определить допустимые проектные параметры, которые лимитируются максимальными перегрузками на приборы системы управления и прочностью твердотопливного заряда ракеты с РДТТ.

Примем максимальную перегрузку с равной 15 [2], тогда

.

Таким образом, с учетом этих параметров, определим максимально допустимое значение давление в «подракетном» пространстве :

– Площадь критического сечения сопла

– Массовый расход РДТТ (примем на старте )

кг/с;

где [2].

– Площадь донного среза

.

– Скорость истечения газа из двигателя

;

где – давление в «подракетном» пространстве [3].

– Общая тяга ракеты с учетом давления в «подракетном» пространстве

.

Составим систему уравнений

;

отсюда, допустимое давление .

Составленная система уравнений и последующие уравнения решены в математическом пакете MathCAD.

Теперь, рассмотрим случай, когда объем «подракетного» пространства замкнут, т.е. выход газа из объема отсутствует (рис.2.2).

Рисунок 2.2. «Замкнутая» схема

– Начальный объем «подракетного» пространства

.

– Теплоемкость газа примем

;

– Температура газа в «подракетном» пространстве

,

где =+[2].

– Давление будем определять из условия

[4];

где , т.к. система замкнутая - примем , V – скорость выхода ракеты из ШПУ.

– Из уравнения суммы всех сил действующих на ракету , определим V

;

– Соответственно, высоту вертикального выхода ракеты из ШПУ будем определять из уравнения

;

– Составим систему уравнений

Примем .

Решая данную систему уравнений, построим график изменения давления от времени (рис.2.3) и определим максимальное значение давления.

Рисунок 2.3. Диаграмма

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

Рисунок 2.4. Диаграмма

Рисунок 2.5. Диаграмма

Рисунок 2.6. Диаграмма

Рисунок 2.7. Диаграмма

Рисунок 2.8. Диаграмма

Рисунок 2.9. Диаграмма

Из графика видно, что при t=0.0135секунд будет достигнуто максимальное давление в ШПУ .

Из расчета было определено, что , поэтому необходимо снизить давление за счет установки дополнительного выхода газа из «подракетного» пространства (рис.2.10).

Рисунок 2.10. Схема ШПУ с отводом газа

Проведем расчет при t=0.0135с, добавив в уравнение

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

.

В результате решения получим, что допустимый массовый расход отвода кг/с.

Определив допустимый массовый расход отвода, спроектируем оптимальное сечение отвода.

Расчет будем проводить также, при t=0.0135с.

Исходные данные для расчета, полученные при нахождении :

все остальные параметры указаны в общих исходных данных.

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

– Определим среднюю плотность газа из уравнения Менделеева-Клайперона

.

– Динамический коэффициент вязкости из формулы Милликена

[5];

– Кинематический коэффициент вязкости

.

– Площадь выходного сечения отвода газа

.

– Для определения числа Re, найдем осредненную по всему объему скорость газа

.

– Число Рейнольдса Re

.

– Из закона сопротивления Блазиуса, найдем параметр трения

.

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

, [6]

где – коэффициент сопротивления, отнесенный к гидравлическому диаметру; U – смоченный периметр.

Решая систему уравнений (1) – (5), получим, что d₂= 0.586м, т.е. зазор равен 8 мм.

Выводы:

Необходимо отметить, что расчет был проведен в программе MathCAD в интегральной форме с помощью блока Given-Find. Выявлено, что расчет данной модели достаточно неустойчив, т.к. число сложных нелинейных уравнений в системе стремится к числу, при котором MathCAD может давать сбои (5-7 уравнений) [24]. Поэтому, необходим подбор параметров расчетной схемы, таких как точность (TOL) и, в особенности, начальные приближения искомых величин. При расчете удалось добиться небольшой погрешности из-за вычислительной схемы Δ_max=0.05% (рис.2.11).

Для сравнения был проведен расчет в дифференциальной форме используя метод Рунге-Кутта, который, как и предполагалось, показал еще большую неустойчивость расчета и, как результат, большую погрешность (Δ_max=23%) (рис.2.12). Так как в системе решаются сложные дифференциальные уравнения с быстро изменяющимися и осцилирующими функциями, то здесь проявилась сложность с подбором расчетного шага. Поэтому, данная расчетная схема в экстремумах функции дала максимальные погрешности.

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

Рисунок 2.11. Диаграмма погрешности вычислений в интегральной модели

Рисунок 2.12. Диаграмма погрешности вычислений в дифференциальной модели

Примечание: Погрешность вычислений Δ, определялась из уравнения давлений

– для интегральной схемы

;

– для дифференциальной схемы

.

Из рис.2.12 видно, что при 0,02 секунды погрешность максимальна (23%), что говорит о том, что параметр p(t) в данной точке не отражает действительность (рис.2.13), а значит, говорит о не пригодности дифференциальной схемы в данном случае.

Рисунок 2.13. Дианрамма , рассчитанная в дифференциальной схеме

2.4 Расчет параметров выхода ракеты при заданной схеме

Теперь, определим параметры выхода ракеты для схемы с отводом (рис.2.10) с учетом вычисленных требований по давлению в «подракетном» пространстве.

Составим систему уравнений движения ракеты в ШПУ с учетом отвода газа.

;

На основе полученных данных при решении системы уравнений, построим графики зависимости параметров p, P, V, T₂, U_ист., от времени. Численные результаты приведены в таблице 2.1.

Рисунок 2.14. График p(t)

Рисунок 2.15. График V(t)

Рисунок 2.16. График P(t)

Рисунок 2.17. График T₂(t)

Рисунок 2.18. График (t)

Рисунок 2.19. График H(t)

Рисунок 2.20. График U_ист(t)

Рисунок 2.21. График n_x(t)

Из рис.2.14 – рис.2.21 видно, что процессы, протекающие в подракетном пространстве после достижения максимума давления, меняются менее значительно, чем в случае с замкнутым пространством.

Таблица 2.1. Сводная таблица параметров старта

t, с p, МПа V, м/с P, Па U, м/с T, ºК , кг/с H, м nx
0	9.96E+04	0	3.53E+04	2.50E+03	1.86E+03	0.00E+00	0	2.935
0.01	6.70E+05	1.194	1.58E+05	1.97E+03	2.41E+03	3.999	0.012	13.169
0.02	7.38E+05	2.629	1.73E+05	1.93E+03	2.45E+03	4.515	0.053	14.398
0.03	7.59E+05	4.055	1.78E+05	1.92E+03	2.46E+03	4.677	0.122	14.778
0.04	7.66E+05	5.46	1.79E+05	1.92E+03	2.46E+03	4.734	0.218	14.914
0.05	7.68E+05	6.841	1.80E+05	1.92E+03	2.46E+03	4.749	0.342	14.947
0.06	7.67E+05	8.199	1.79E+05	1.92E+03	2.46E+03	4.741	0.492	14.93
0.07	7.65E+05

9.535

1.79E+05

1.92E+03

2.46E+03

4.722

0.667

14.885

0.08

7.61E+05

10.848

1.78E+05

1.92E+03

2.46E+03

4.696

0.868

14.823

0.09

7.57E+05

12.141

1.77E+05

1.93E+03

2.45E+03

4.665

1.093

14.751

0.1

7.53E+05

13.413

1.76E+05

1.93E+03

2.45E+03

4.632

1.341

14.673

0.11

7.49E+05

14.666

1.75E+05

1.93E+03

2.45E+03

4.597

1.613

14.591

0.12

7.44E+05

15.901

1.74E+05

1.93E+03

2.45E+03

4.562

1.908

14.507

0.13

7.39E+05

17.118

1.73E+05

1.93E+03

2.45E+03

4.526

2.225

14.422

0.14

7.35E+05

18.317

1.72E+05

1.94E+03

2.44E+03

4.49

2.564

14.337

0.15

7.30E+05

19.5

1.71E+05

1.94E+03

2.44E+03

4.454

2.925

14.252

0.16

7.25E+05

20.667

1.70E+05

1.94E+03

2.44E+03

4.418

3.307

14.167

0.17

7.21E+05

21.819

1.69E+05

1.94E+03

2.44E+03

4.382

3.709

14.083

0.18

7.16E+05

22.956

1.68E+05

1.95E+03

2.44E+03

4.347

4.132

14

0.19

7.11E+05

24.079

1.67E+05

1.95E+03

2.43E+03

4.313

4.575

13.918

0.2

7.07E+05

25.188

	1.66E+05	1.95E+03	2.43E+03	4.279	5.038	13.838
0.21	7.03E+05	26.284	1.65E+05	1.95E+03	2.43E+03	4.246	5.52	13.758
0.2273	6.95E+05	28.149	1.64E+05	1.96E+03	2.43E+03	4.189	6.398	13.626

Выводы

Составлена приблизительная схема расчета динамики выхода ракеты из ШПУ с элементами проектирования пусковой установки. В результате расчета определены параметры газа, находящегося в «подракетном» пространстве и параметры движения ракеты от начала движения до ее полного выхода из ШПУ. Кроме того, определены некоторые геометрические характеристики пусковой установки.

Анализ полученных результатов и заключение

В первом приближении расчетное время выхода ракеты из ШПУ составило примерно 0.23 секунды. Для более точного решения требуется уточнить характеристики стартового РДТТ, включить в расчет действие аэродинамических сил при выходе ракеты из шахты. Но следует отметить, что при более тщательном рассмотрении вопроса увеличится число уравнений в системе и, следовательно, сложность вычислений.

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

Таким образом, результаты расчета могут быть применены при определении аэродинамических нагрузок, нагрузок на ШПУ, а также деформаций, углов отклонений конструкции ракеты, которые могут возникать при выходе из пусковой установки.

3. МЕТОДИКА ЧИСЛЕННОГО ЭКСПЕРИМЕНТА

3.1 Постановка задачи

Определить и уточнить методику численного эксперимента для решения поставленных задач, в том числе: динамика выхода ракеты из ШПУ и определения АДХ ракеты.

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

3.2 Общие сведения о проведении численного эксперимента газовой динамики

Первым этапом проведения эксперимента является возникновение или постановка физической задачи. Правильная постановка задачи в значительной мере предопределяет успех всей проводимой работы. Моделирование опыта можно рассматривать как второй этап. Это может быть либо численное моделирование, либо моделирование физических процессов. Каждый из этих методов имеет свои особенности. Численное моделирование дает количественное выражение закономерностей, присущих математической модели. Напротив, с помощью моделирования физических процессов наблюдается аналог действительности или даже сама действительность. Таким образом, полезность расчета ограниченна обоснованностью математической модели [7].

Численное моделирование задач газовой динамики осуществляется по следующему алгоритму:

1. Построение геометрии модели (тело + расчетная область).

Построение геометрии осуществляется в специализированных пакетах таких, как Solid Works, KOMPAS, Unigraphics, ProEngineer и др. Построенная объемная модель, специального формата, импортируется в сеточный генератор.

2. Построение блочной сетки.

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

3. Подготовка модели (задание свойств среды, начальных условий, граничных условий, параметров решателя).

4. Запуск решателя.

Здесь необходимо проверить “сходимость” решения.

5. Просмотр результатов в модуле Post.

Численное решение, производимое в решателе, представлено на рис.3.1.

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

На практике дискретизация производных по времени осуществляется почти исключительно с использованием разностных методов [8]. При дискретизации пространственных производных используется, как правило, метод конечных разностей, конечных элементов, конечных объемов [8, 9, 10, 11].

Рисунок 3.1. Структура численного решения

При применении данных методов алгебраические уравнения связывают между собой значения искомых переменных в группе соседних узловых точек (сеточных узлов) [12, 8]. Также подразумевается, что сетка, состоящая из дискретных точек, распределена по всей вычислительной области во времени и в пространстве.

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

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

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

Если учесть, что для большинства задач газовой динамики определяющие уравнения являются нелинейными, процесс построения численного решения обычно ведется посредством итераций (используются методы Ньютона, многосеточные методы, метод сопряженных градиентов). Иначе говоря, решение для каждого искомого переменного в каждой узловой точке последовательно корректируется посредством обращения к дискретизованным уравнениям [9, 14, 8, 10, 15, 13].

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

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

• Основные уравнения газовой динамики

Упомянутые выше основные уравнения, являются уравнениями газовой динамики. Газовая динамика – это раздел механики сплошных сред, описывающий движение жидкостей и газов в рамках модели сплошной среды. Последнее означает, что рассматриваются масштабы явлений, значительно превосходящие длину свободного пробега молекул. В рамках данного подхода все физические законы, а также свойства являются общими как для макрообъектов, так и бесконечно малых объемов.

В наиболее общем случае для задачи газовой динамики требуется решить систему из четырех независимых уравнений, которая носит название системы уравнений Навье-Стокса:

1. Уравнение неразрывности (сохранения массы)

. 9

2. Уравнение количества движения (сохранения импульса)

, 9

где

- тензор напряжений, записываемый в виде

; 9

- дельта-функция Кронекера

. 9

3. Уравнение энергии (сохранения энергии)

, 9

где

, 9

. 9

4. Уравнение состояния

9

Для записи соотношений - использованы следующие обозначения: - давление; - плотность; - скорость; - температура; - время; - полная энтальпия; - статическая энтальпия;   источниковый член для импульса; - источниковый член для энергии; - коэффициент динамической вязкости; - коэффициент теплопроводности; - оператор Гамильтона (набла); - обозначает векторную величину.

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

При рассмотрении некоторых основных дифференциальных уравнений гидродинамики -, можно сделать вывод, что основные переменные подчиняются обобщенному закону сохранения [4, 8]. Если обозначить зависимую переменную , то обобщенное дифференциальное уравнение можно записать в следующем виде:

9

где

- коэффициент диффузии;

- источниковый член.

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

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

При историческом развитии динамики жидкости в рассмотрение был введен ряд классов течений, описываемых значительно более простыми системами, чем указанная выше -. Эти различные классы возникают при пренебрежении или ограничении некоторых свойств течений. Для течений, представляющих практический интерес, соответствующая классификация приведена в таблице 2.1. Классификация проведена по двум параметрам - вязкости и плотности. Несжимаемые течения, как правило, ассоциируются с течениями, скорость которых мала по сравнению со скоростью звука (). Наоборот, для сжимаемых течений (, либо разница температур в потоке велика) требуется рассмотреть полное уравнение неразрывности и учитывать полное уравнение энергии.

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

Таблице 2.1. Классификация течений

Вязкость	Плотность
	Несжимаемые ()	Сжимаемые ()
Невязкие течения ()	Потенциальные течения ()	Газовая динамика ()
Течение в пограничных слоях (вязкость существенна вблизи поверхности)	Ламинарные течения (очень малые числа ) Турбулентные течения (большие числа )	Существенен перенос тепла
Отрывные течения (вязкость существенна везде)	Ламинарные течения (малые числа ) Турбулентные течения (очень большие числа )	Существенен перенос тепла

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

У плохо обтекаемых тел (например, автомобиля) на подветренной стороне возникают области отрывных течений, в которых существенны эффекты вязкости. Если числа Рейнольдса не слишком малы, течения в таких зонах являются турбулентными и часто нестационарными. Обычно для описания отрывных течений необходимо решать полную систему уравнений Навье-Стокса для сжимаемой и несжимаемой жидкостей.

Приведенные выше соображения позволяют сильно упростить решение задачи при соответствующих обоснованных допущениях. Однако даже в подобных идеализированных случаях точное математическое решение существует только для простых тел (пластина, сфера, цилиндр, клин).

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

3.3 Компьютерные пакеты для численного решения задач газовой динамики

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

Естественным шагом в эволюции численного моделирования динамики жидкости и газа стало создание расчетных пакетов (CFD-пакетов или комплексов), ориентированных на широкую аудиторию пользователей – научных работников, студентов, инженеров и т. д. В таком виде математический аппарат, заключенный в численные методы, стал действительно универсальным, а с учетом стремительного развития вычислительной техники и мощным средством в проведении численных расчетов по газовой динамике. Кроме того, при использовании CFD-пакетов становится необязательным обладать глубокими знаниями по численным методам, способам дискретизации и т.п.

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

1. «Тяжелые» - комплексы высокого класса, подходящие как для научных, так и инженерных расчетов, способные решать самые сложные задачи с учетом большого количества эффектов и использованием широкого набора математических подходов, в том числе специфических. К классу «тяжелых» относятся лидеры среди коммерческих CFD-пакетов – ANSYS CFX (ANSYS, Inc.), Star-CD (CD-adapco), FLUENT (ANSYS, Inc.совместно с Fluent, Inc.). Все они содержат большое число моделей турбулентности, способны решать задачи различной сложности с учетом горения, химических реакций, многофазных потоков, поддерживают различные типы сеток и т. д.

2. Среднего класса. Предназначены, главным образом, для расчетов инженерного уровня сложности. Набор используемых моделей также может быть достаточно широким. К этому разряду можно отнести COSMOSFloWorks (Solid Works Co.), STAR-CCM+ (CD-adapco), ANSYS FLOTRAN (ANSYS, Inc.).

3. «Легкие» - CFD-комплексы, использующие алгоритмы невысокой точности (используются, например, в качестве учебно-методических), либо имеющие узкую направленность расчета (специально созданные под определенную проблематику).

Следует отметить тот факт, что подавляющее большинство CFD-кодов, реализованных в программах, основано на использовании МКО в различных вариациях.

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

3.4 Общие сведения о методе контрольного объема

Выбранный для примера программный комплекс Ansys CFX v.10 основан на конечно-объемном методе (МКО) решения уравнений гидродинамики. Понимание метода позволит безошибочно провести численный эксперимент, несмотря на то, что в CFX МКО дополнен и утонен рядом других уравнений. Поэтому, можно рассмотреть классический метод контрольного объема, описанный Патанкаром.

В МКО расчетную область разбивают на N-е число непересекающихся контрольных объемов (рисунок 3.2) таким образом, что каждая узловая точка содержится в одном контрольном объеме. Дифференциальное уравнение интегрируют по каждому контрольному объему. Для вычисления интегралов используют кусочно-непрерывные функции, которые описывают изменение зависимой переменной (например, одной из составляющих скорости) между сеточными узлами. В результате находят дискретный аналог дифференциального уравнения. Дискретный аналог представляет собой алгебраическое уравнение,

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

1. Выбор начального приближения или оценка значений Ф во всех узловых точках.

2. Расчет предварительных значений коэффициентов в дискретном аналоге на основе начального профиля Ф.

3. Решение линейной системы алгебраических уравнений, дающее новые значения Ф.

4. Возврат ко второму этапу и повторение процесса до тех пор, пока дальнейшие приближения (итерации) перестанут давать сколько-нибудь существенные изменения в значениях Ф.

Рис.3.2. Контрольный объём (заштрихованная область), для двумерного случая и обозначения для алгоритма в декартовой сетке

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

Ansys CFX позволяет проводить расчеты на смешанных сетках, состоящих из различных типов элементов: тетраэдров, призм, клиновидных элементов и гексаэдров (рис.3.3).

Рисунок 3.3. Ячейки а) – гекса, б) – тетра, в) – призма, г) - пирамида

При расчете стационарных вариантов процесс итерации по времени завершается при достижении требуемого уровня сходимости, определенного пользователем. Для расчета переходного режима итерационная процедура обновляет нелинейные коэффициенты на каждом временном шаге (цикл для коэффициентов), в то время как внешний цикл приближается к решению по времени. [12,16] Для наглядности, схема описанного алгоритма изображена на рис.3.4.

Рисунок 3.4. Схема алгоритма

Выводы

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

Для проведения численного эксперимента был выбран наиболее оптимальный программный прикладной пакет – Ansys CFX, основанный на методе контрольного объема.

4. РАСЧЕТ ДИНАМИКИ ВЫХОДА РАЕТЫ ИЗ ШАХТНОЙ ПУСКОВОЙ УСТАНОВКИ (ЧИСЛЕННЫЙ ЭКСПЕРИМЕНТ)

4.1 Постановка задачи

На основе исходных данных и расчетах, проведенных в разделе 2, разработать стратегию численного эксперимента газодинамической задачи выхода ракеты из шахты. На основе проведенной работы сделать выводы.

4.2 Стратегия постановки численного эксперимета

Газодинамическая задача является временной, поэтому решение ее будет проводиться в зависимости от времени. Для этого необходимо построить несколько моделей, в которых будет меняться объем «подракетного» пространства, в зависимости от движения ракеты. При достижении давления, которое обеспечивает движение ракеты на величину H (величина Н выбирается из условия точности расчета), решается новая модель, начальные условия в которой соответствуют параметрам в последний момент предыдущего расчета. При таком подходе определяются только параметры газа, а движение ракеты определяется из уравнения

и .

• Построение твердотельной модели

Как уже было сказано, модель можно построить в любом графическом редакторе, в том числе KOMPAS, Solid Works, Unigraphics, ProEngineer. В данном случае выбран пакет Solid Works, в котором была создана модель и расчетная область. Так как модель симметрична, то можно ограничиться построением сектора (рис.4.1).

Конечным результатом построения должны стать сохраненные файлы формата Parasolid (текстовые).

Рисунок 4.1. Схемы твердотельных моделей (нижняя часть)

• Построение сетки

Одним из самых сложных этапов моделирования, является построение сетки расчетной области. Качество сетки влияет на сходимость решаемой задачи и самое главное, на правдивость полученных результатов. Необходимо помнить, что сетка должна иметь больших различий по соотношениям объемов ячеек. При уменьшении объема ячеек – увеличивается точность, но при этом и увеличивается время расчета. Но, зная, что наиболее важной характеристикой точности является величина ячеек на границе с телом, то объем ячеек сетки будем уменьшать от краев расчетной области к границам рассматриваемого тела.

Сетка для схемы изображенной на рис.4.1 создавалась в специализированном сеточном генераторе ICEM CFD компании ANSYS. Для быстроты построения сетки были использованы тетраэдрические ячейки (рис.4.2). Но, несмотря на быстроту построения этот метод обладает плохой сходимостью при расчете и снижает скорость самого расчета.

Рисунок 4.2. Тетраэдрическая сетка

• Подготовка решателя расчетного комплекса и проведение расчета

Для экспортированной в Ansys CFX сеточной модели были определены граничные условия. Определение границ изображено на рис. 4.3.

Рисунок 4.3. Определение границ расчетной области

В качестве параметров расчета были использованы следующие:

1. Модель газа: модель газа, приближенного к параметрам реального воздуха при температуре 25 градусов - Air at 25C.

2. Модель турбулентности: SST с решением полного уравнения энергии.

3. Модель стенки: адиабатическая стенка (Adiabatic Wall) с учетом поверхностных напряжений трения (без проскальзывания) – No Slip.

4. Схема решения уравнений: Временной.

5. Параметры сходимости решения: максимальное число итераций – 30.

6. Временной шаг – 0,005 секунд

7. Граничные условия задавались следующим образом:

8. Граница входа Inlet: статическая температура (Static Temperature) – 2500 К, массовый расход через вход – 5,36 кг/с

9. Свободный выход Open: статическое давление на выходе (Static Pressure for Entrainment) - 99600 Па, статическая температура (Static Temperature) – 298 К.

10. Граница симметрии – Symmetry.

4.3 Анализ полученных результатов

В результате расчета были получены значения параметров газа в «подракетном» пространстве. Картины параметров газа для t=0.02 c представлены на рис.4.4. и рис.4.5.

Рисунок 4.4. Распределение давлений в «подракетном» пространстве

Рисунок 4.5. Линии тока газа в «подракетном» пространстве

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

Таблица 4.1. Результаты расчетов

t, c	Давление в «подракетном» пространстве (p), Па		Δ, %
	Аналитический расчет	Численный эксперимент
0	99600	99600	0
0.005	4,00Е+05	3,80Е+05	5
0.01	6,70E+05	5,13Е+05	23
0.015	7,10Е+05	4,70Е+05	33
0.02	7,38E+05	3,23Е+05	45

Такие неудовлетворительные результаты могли возникнуть из-за ряда причин, в том числе:

1. использовалась модель воздуха;

2. скорость истечения (340м/с) отличается от расчетной (2500м/с);

3. количество уточнений (30) на каждом промежутке времени недостаточно для используемой сеточной модели.

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

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

Выводы

Разработана схема численного эксперимента газодинамической задачи выхода ракеты из пусковой шахты. В дальнейшем, после доработки расчетной модели, можно будет получить результаты которые можно сравнивать с ручным расчетом и с данными физических экспериментов.

5. ОПРЕДЕЛЕНИЕ АЭРОДИНАМИЧЕСКИХ НАГРУЗОК

Введение

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

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

Таблица 5.1. Характер и значения скорости ветра

Сила ветра по шкале Бофорта			Визуальная оценка действия ветра
Баллы	Скорость, м/с	Характерис-тика ветра
0	0 – 0,2	Штиль	Дым из труб поднимается отвесно
1	0,3 – 1,5	Тихий	Дым слегка отклоняется
2	1,6 – 3,3	Легкий	Движение воздуха ощущается лицом – шелестят листья; начинают шевелиться флаги, флюгер
3	3,4 – 5,4	Слабый	Колеблются тонкие ветки; развеваются флаги; начинается легкий перенос снега по поверхности покрова
4	5,5 – 7,9	Умеренный

Поднимается пыль и бумажки; колеблются небольшие сучья; снегопад переходит в метель
5	8 – 10,7	Свежий	Качаются тонкие стволы деревьев; на воде появляются белые «барашки»
6	10,8 – 13,8	Сильный	Качаются толстые сучья; гудят провода; шум ветра слышен в домах
7	13,9 – 17,1	Крепкий	Качаются стволы деревьев; затрудняется движение
8	17,2 – 20,7	Очень крепкий	Ломаются сучья деревьев; колеблются средние деревья; очень трудно идти против ветра
9	20,8 – 24,4	Шторм	Ломаются толстые сучья и небольшие деревья; разрушаются дымовые трубы; сбрасывается черепица; на море высокие волны
10	24,5 – 28,4	Сильный шторм	Разрушаются строения; деревья вырывает с корнем; ломаются телеграфные столбы
11	28,5 – 32,6	Жесткий шторм	Большие разрушения
12	32,7 – 36,9	Ураган	Опустошительные разрушения
17	58,6	–	–

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

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

Обычно транспортные космические системы имеют цилиндрические формы составляющих ее блоков. Поэтому для оценки сил, действующих на ее элементы вблизи земли, следует иметь ввиду явления, присущие поперечному (или близкому к нему) обтеканию цилиндра при дозвуковых скоростях. Дело в том, что благодаря проявлению вязкости воздуха при дозвуковых скоростях могут реализовываться различные режимы обтекания с большим или меньшим сопротивлением в зависимости от скорости ветра и диаметров отдельных частей ракеты. [17], [18]

5.1 Постановка задачи

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

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

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

По результатам расчета, требуется построить эпюры моментов и перерезывающих сил.

5.2 Подготовка и проведение численного эксперимента

• Построение твердотельной модели

Расчетная модель была построена в пакете Solid Works. В данном случае, для исключения влияния стенок на результат, необходимо учитывать площадь заполнения телом расчетной области, требуется, чтоб это параметр был примерно равен 1/8 площади расчетной области.

• Построение блочной сетки

Блочная сетка, для схемы изображенной на рис.5.1, создавалась в специализированном сеточном генераторе ICEM CFD компании ANSYS.

Рисунок 5.1. Схема моделирования

В результате построения сетки с применением блока «O-grid», т.е. грани ячеек внутреннего блока расположены перпендикулярно к касательным генерируемого тела, была создана расчетная модель (рис.5.2)

Рисунок 5.2. Сгенерированная блочная сетка

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

• Подготовка решателя расчетного комплекса и проведение расчета

Для примера примем, что на ракету в перпендикулярном направлении дует «свежий» ветер, эта характеристика соответствует 10 м/с. Скорость и высоту выхода ракеты возьмем из предыдущего расчета. Для данной схемы (рис.5.1) проведем расчеты в П.П.П. Ansys CFX.

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

;

.

Для предварительной оценки проведем шесть расчетов с параметрами, указанными в таблице 5.2. Параметры расчетной среды приведены в таблице 5.3.

Таблица5.2. Расчетные параметры

№	H, м	, м/с	, град
1	0,492	12.93	50.65
2	1,341	16.73	36.7
3	2,225	19.82	30.29
4	3,307	22.95	25.82
5	4,575	26.07	22.55
6	6,398	29.87	19.55

Таблица 5.3. Параметры расчетной среды

Параметр	Значение
давление, Па	99600
температура, K	298
плотность воздуха,	1,2

Данный расчет проведен с использованием SST-модели турбулентности, так как данная модель более точна и надежна для широкого класса потоков (т.е., потоков подверженных градиентам давления, обтекание профилей), чем стандартная k-w модель турбулентности. [16]

Напомним, что жидкость, находящаяся в расчетной области, при M << 1 будет считаться несжимаемой.

• Получение результатов

В результате расчета были получены теневые картины обтекания ракеты. Приведем такие картины для расчетного примера №6 (рис.5.3, рис.5.4).

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

Рисунок 5.3. Теневая картина распределения скоростей потоков (расчетный пример №6)

Рис.5.4. Теневая картина распределения давления (расчетный пример №6)

Полученные результаты следует сравнивать с теоретическими и представленными в литературе подобными экспериментами. Необходимо оценивать схожесть с физическими законами.

5.3 Анализ полученных результатов

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

.

Образуемые от этих потоков силы, складываются векторно, но можно предположить, что аэродинамические эффекты, возникающие от продольных и поперечных сил, будут накладываться и присутствовать на ракете выходящей из шахты. Поэтому, для проверки правильности расчета, проведем дополнительные численные эксперименты с последовательным исключением одной из составляющих потока (,).

5.3.1 Вертикальная составляющая потока ()

На ракету действует под углом атаки поток = 24,079 м/с, равный скорости ракеты для расчетного случая № 4 (H=4,575 м) (I случай). Расчет проводится без учета донного давления.

Данный расчет обладает хорошей сходимостью решения (рис.5.5). Решение было получено через 19 итераций, что составило 20 мин.

Рисунок 5.5. График сходимости для расчета при = 24,079 м/с, =0

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

Рисунок 5.6. Линии тока от продольной составляющей () (I сл.)

Определим аэродинамические коэффициенты сил, действующих в данной расчетной схеме, и построим графики распределения давлений по поверхности ракеты (рис.5.7).

Рисунок 5.7. Диаграмм распределения коэффициента давления по контуру (ось Y) ракеты в продольном направлении

Ввиду отсутствия симметрии в расчетной области наблюдается незначительное неравномерное распределение давлений по контуру ракеты в поперечном направлении (Δ С_p=). Антисимметрия вызвана наличием в граничных условиях вертикальной стенки входа потока (=0) (рис.5.1).

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

Коэффициенты продольной и нормальной сил

где q – скоростной напор, Па;

S – площадь миделя ракеты, м²;

- суммарная продольная сила, Н;

- суммарная нормальная сила, Н.

Замечание: При выходе ракеты из шахты на нее действуют силы в проекции на ось Х от давления в «подракетном» пространстве, от атмосферного давления и от скоростного напора. Так как влияние давлений атмосферы и «подракетного» пространства были учтены при расчете основных параметров, то следует вычитать действие силы от атмосферного давления на ракету в проекции на ось Х из общей силы F_x.

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

Рисунок 5.8. Линии тока от продольной составляющей () (с учетом донной части - IIсл.)

Результаты вычислений и сравнение с моделью без донной части представлены в таблице 5.4

Таблица 5.4. Результаты вычислений

	I сл.	II сл.	Δ, %
Корпус	872,8	862,3 (Н)	1,2%
Донная часть	–	-851 (Н)	–
	872,8	11,3	–

что составляет около 15% различия от физических экспериментов ().

5.3.2 Поперечная составляющая потока ()

Проведем сравнительный численный эксперимент обтекания ракеты, находящейся вертикально к ветровому потоку (=20м/с) (I случай).

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

Характер обтекания и силы, действующие на круговой цилиндр при углах атаки, близких к 90º, является функциями чисел Re и M. (При поперечном обтекании аэродинамические силы относятся к максимальной площади продольного сечения S = l × d.) При небольших скоростях обтекания проявляются силы вязкости воздуха, и основным критерием подобия является число Re, с которым связано так называемое явление «кризиса обтекания» – резкое уменьшение сопротивления при достижении критической величины числа Re_кр ≈ 2,5 ∙, подсчитанного по диаметру цилиндра. В зависимости от скорости ветра и диаметра соответствующего участка ракеты числа Re могут превосходить значения Re = 3 ∙ (при скорости ветра >20м/с). [1]

Так как высота выхода ракеты мала при Re_d < 2,5 ∙, то будем рассматривать случай, когда этот показатель выше.

Максимальное давление поток оказывает в передней критической точке (рис.5.10), после чего давление убывает, так как поток, обтекающий лобовую часть цилиндра, непрерывно разгоняется. При давление становится равным статическому давлению невозмущенного потока. При достижении критических чисел Re пограничный слой потока, обтекающего цилиндр, перед областью отрыва становится турбулентным. Поскольку турбулентный слой содержит больше энергии, чем ламинарный, он получает возможность обтечь контур цилиндра по большей его части. Отрыв потока в ряде случаев доходит до , в результате чего лобовое сопротивление составляет .

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

Данные численного расчета сравнивались с физическим экспериментом (рис.5.8, рис.5.9, рис.5.10). Расчет проводился с = 20 м/с, что соответствует Re = 6,7 ∙.

При расчете были применены несколько моделей турбулентности, в том числе: модель переноса касательных напряжений (Shear-Stress Transport) «SST» («k-ω») и «RNG» модель, основанная на стандартной «k-ε» модели. Из рисунков 5.8 и 5.9 видно, что модель «RNG» наилучшим способом подходит к моделированию данной проблемы, так как данная схема имеет возможность учета эффектов закрученности и вращения, а, в свою очередь, вращение и закрученность осредненного потока оказывает существенное влияние на структуру турбулентности, в частности на турбулентную вязкость[16]. Но надо заметить, что применение модели турбулентности RNG требует расчетного времени в два раза больше чем в случае со схемой SST.

Рисунок 5.8. Распределение давления по поперечному сечению цилиндра при Re_d=

Рис.5.9. Распределение коэффициента давления по поперечному сечению цилиндра при Re_d=

А)

Б)

А) – теоретическая схема обтекания Б) – физический эксперимент

В) численный эксперимент

Рисунок 5.10. Обтекание поперечного сечения цилиндра при Re_d=

Как видно из рис.5.8 – рис.5.10 результаты расчета в CFX сходятся с представленными физическими экспериментами.

Но нужно заметить, что распределение давление по контуру ракеты (рис.5.11) в поперечном направлении, как показано на рис.5.9, можно наблюдать лишь в центральной части ракеты (1м – 3м от Земли), так как в передней части ракеты (3м – 4,575м) большое влияние оказывает закругленная носовая часть, а в нижней (0м – 1м) на распределение давления оказывает стенка расчетной области, имеющая возможность пропускать поток.

Рисунок 5.11. Распределение коэффициента давления по контуру ракеты в продольном направлении (Ось Y)

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

Сравним коэффициенты лобового сопротивления для элемента ракеты (цилиндра) и всей ракеты в целом с представленными в литературе.

Из диаграммы (рис.5.8., рис.5.9), определим коэффициент лобового сопротивления для бесконечно длинного цилиндра:

CFX «RNG»:

CFX «SST»:

физический эксперимент:

Таким образом, расхождение составило 11% («RNG»), 13% («SST»).

Определим лобовое сопротивление всей ракеты

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

5.3.3 Сложение составляющих потока (;)

В итоге, после проведенных в разделах 5.3.1 и 5.3.2 расчетах, были получены удовлетворительные результаты. Попробуем с помощью простейших статистических методов оценки параметров распределения получить картину распределения коэффициента давления по контуру ракеты в продольном направлении (для γ=0 и γ=180) при сложении продольного и поперечного потока (рис.5.12).

Рисунок 5.12. Распределение давления по контуру ракеты в продольном направлении (I сл.)

При сложении применялось уравнение , где – коэффициент давления рассчитанный при продольном обтекании в i-той точке; – коэффициент давления рассчитанный при поперечном обтекании в i-той точке.

Полученный таким образом график С_P не может быть точным, но может показывать характер обтекания. Для уточнения графика необходимо использовать более сложные статистические методы, особенно в точках, которые сильно отклоняются от среднего значения функции, а именно две характерные точки в носовой части – точки максимума и минимума (А и Б).

Теперь, проведем проверочный расчет при =24м/с и = 20м/с () с помощью численного эксперимента.

Данная задача так же, как и случай с поперечным обтеканием ракеты, решается только экспериментальным путем, хотя существуют полуэмпирические методы, дающие большие погрешности при расчете. Опыт показывает, что зависимости подъемной силы от угла сохраняют линейный характер только при малых значениях этих углов. По мере их роста действительные зависимости С_y ()все сильнее отклоняются от линейных.

Степень нелинейности определяется числом маха М, Re и геометрическими формами ЛА. Все это приводит к тому, что при углах приходится подсчитывать с учетом нелинейных составляющих.

В данном случае расчетная схема выглядит следующим образом (рис.5.12):

Рисунок 5.12. Схема обтекания тела вращения под большим углом

Рисунок 5.13. Линии тока обтекания ракеты под углом =40

Рисунок 5.14. Распределение давления по контуру ракеты в продольном направлении полученные наложением и численным экспериментом в CFX

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

Рисунок 5.16. Отрывные течения на подветренной стороне ракеты

Как видно из рис.5.13 – 5.16, на подветренной стороне ракеты образуются сложные вихревые течения. Под действием ветровой нагрузки происходит срыв потока, ее действие максимально в пределах 3 метров от нижней точки ракеты. В носовой части на картину течений оказывает продольная сила, благодаря чему оторвавшийся поток возвращается к контру ракеты. Кроме того, в результате увеличения турбулентности, изменилась точка отрыва потока с поверхности до .

Определим основные аэродинамические коэффициенты:

Запишем формулы для нахождения сил X и Y в поточной системе координат (рис.5.17):

где

- суммарная продольная сила, Н;

- суммарная нормальная сила, Н.

Определим коэффициенты в поточной системе координат:

где

- угол атаки, в градусах;

q – скоростной напор, Па;

S – площадь миделя ракеты, м.

Запишем формулы для нахождения С_x и C_у:

Таким образом, на этапе отработки, нами были получены удовлетворительные результаты, повторяющие, в достаточной степени, физические эксперименты. Поэтому, расчеты подобные описанным выше будем считать также удовлетворительными. Далее, на той же блочной сетки, модели турбулентности и схожих физических параметрах проведем расчет для определения моментов и перерезывающих сил ракеты выходящей из ШПУ.

5.4 Определение моментов

По результатам проведенных расчетов получены силы, действующие на ракету выходящую из ШПУ и центры давления для случаев, указанных в таблице 5.2.

В расчетах этой части была применена модель турбулентности «SST» для повышения скорости расчета, кроме того, как было выяснено в п.5.3.2, эта модель мало влияет на интегральные характеристики, которыми будем оперировать в данном разделе.

Также, для получения полной информации о силах, действующих на ракету, построим графики распределения давление по контуру в продольном направлении (рис.5.18), на которых хорошо видна эволюция влияний внешних факторов на выходящую ракету.

H=0,492м

H=1,341м

H=2,225м

H=3,307м

H=4,575м

H=6,398м

Рисунок 5.18 Распределение давлений по контуру ракеты в продольном направлении

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

Момент определялся относительно верхней кромки ШПУ. Ввиду малости участков, центр давления принимается в середине участка по длине. Из суммарного момента определялся общий центр давления ракеты

Полученные в CFX силы экспортировались в Excel и записывались следующим образом (таблица 5.4):

Таблица. 5.4. Табличная запись экспортированных сил

Номер уачстка

	Расстояние верхней кромки ШПУ, мм		, Н	, Н	, Н∙м
. . .	. . .	. . .	. . .	. . .	. . .

Общие результаты вычислений приведены в таблице 5.5.

Таблица 5.5. Результаты вычислений в CFX

№	H, м	, м/с	, град	q, Па	Fx, H	Fу, H	Xц.д., м	M, Н∙м
1	0,492	12.93	-50.65	99.086	-850	-15.007	0.318	-4.77
2	1,341	16.73	-36.7	165.856	-837.8	-37.49	1.017	-38.14
3	2,225	19.82	-30.29	232.89	-847.8	-49.28	1.858	-91.54
4	3,307	22.95	-25.82	312.337	-842.6	-69.202	2.68	-185.52
5	4,575	26.07	-22.55	402.782	-864.6	-91.756	3.83	-352.04
6	6,398	29.87	-19.55	528.767	-833	-112.33	4.2427	-476.59

Для каждого этапа выхода ракеты из ШПУ построим эпюры перерезывающих сил и моментов

Рисунок 5.20. Аэродинамические силы, действующие на ракету (Н=0,492м)

Рисунок 5.21. Аэродинамические силы, действующие на ракету (Н=1,341м )

Замечание: На рис.5.20 - 5.24 эпюры перерезывающих сил носят информативный характер.

Рисунок 5.22. Аэродинамические силы, действующие на ракету (Н=2,225м)

Рисунок 5.23. Аэродинамические силы, действующие на ракету (Н=4,575м)

Рисунок 5.24. Аэродинамические силы, действующие на ракету (Н=6,4м)

Как видно из приведенных выше рисунков и табл., максимальный момент увеличивается с выходом ракеты из ШПУ. В последние секунды выхода ракеты он достигает 476.59 Н∙м.

На рис.5.25 показано распределение моментов по оси ракеты в зависимости от скорости ветра. Видно, что моменты, с увеличением ветра в два раза, увеличиваются в такойже пропорции.

Рисунок 5.25. Распределение моментов по оси ракеты для H=4.575 и скорости ветра 5, 10, 20 м/с

Кроме моментов были определены аэродинамические коэффициенты для ракеты выходящей из пусковой шахты (рис.5.26). При расчете АДК значение аэродинамических сил относилось к площади поперечного сечения корпуса, а коэффициенты моментов – к той же площади и полной длине центрального корпуса.

Рисунок 5.26. АДК в зависимости от угла атаки набегающего потока

Из рис.5.26 – характер изменения АДК соответствует характеру изменения АДК, представленных в литературе [17].

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

Таким образом, данная схема может быть использована для приблизительной оценки нагрузок. Общая погрешность в расчетах относительно физических экспериментов составила около 10-15%.

Выводы

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

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

ВЫВОДЫ

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

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

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

ЗАКЛЮЧЕНИЕ

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

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

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

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

1. URL: WWW.ARMS-EXPO.RU

2. Павлюк Ю.C. Баллистическое проектирование ракет: Учебное пособие для вузов. — Челябинск: Изд. ЧГТУ, 1996. — 92 с.

3. Cидельников Р.В. Теория полета: Краткий конспект лекций. – Челябинск: Изд-во ЮУрГУ, 2003. – 73 с.

4. Есин В.И. Пневмогидравлические системы ПГС и автоматика ракет: Учебное пособие. – Челябинск: ЧПИ, 1988. – 129 с.

5. Сидельников Р.В. Аэрогидрогазодинамика: Краткий конспект лекций. – Челябинск: Изд. ЮУрГУ, 2003. – 103 с.

6. Теория пограничного слоя. Шлихтинг Г., перев. с немецкого, Главная редакция физико-математической литературы издательства «Наука», Москва, 1974, 712 с.

7. Преображенский В.П. Теплотехнические измерения и приборы: Учебник для вузов. – 3-е изд., перераб. – М.: «Энергия», 1978. 704 с.

8. Флетчер К. Вычислительные методы в динамике жидкостей: В 2-х томах: Т. 1: Пер. с англ. – М.: Мир, 1991. – 504 с, ил.

9. Роуч П. Вычислительная гидродинамика. – М.: Мир, 1980. – 616с., ил.

10. Chung T. J. Computational Fluid Dynamic. – Cambridge: Cambridge University Press, 2002. – 1021 p.

11. Ferziger J. H., Perić M. Computatational Methods for fluid Dynamics. – 3., rev. ed. – Berlin: Springer, 2002. – 423 p.

12. Патанкар С. Численные методы решения задач теплообмена и динамики жидкости: Пер. с англ. – М.: Энергоатомиздат, 1984. – 152 с., ил.

13. Versteeg H. K., Malalasekera W. An Introduction to Computational Fluid Dynamics. The finite volume method. – London: Longman Scientific & Technical, 1995. – 252 p.

14. Станкова Е. Н., Затевахин М. А. Многосеточные методы. Введение в стандартные методы. – СПб.: Институт высокопроизводительных вычислений и информационных систем, 2003. – 47 с.

15. Ferziger J. H., Perić M. Computatational Methods for fluid Dynamics. – 3., rev. ed. – Berlin: Springer, 2002. – 423 p.

16. CD Adapco Group. Tutorial (STAR-CD VERSION 3.22.).- 2004

17. Петров К.П. Аэродинамика транспортных космических систем. – М.: Эдиториал УРСС, 2000. – 368 с.

18. Маликов В.Г. Наземное оборудование ракет. М.: Воениздат, 1971. – 304 с.

19. ANSYS ICEM CFD v.11. Release. Help manual.

20. ANSYS CFX v.10. Release. Theory Reference.

21. Гигиенические требования к видеодисплейным терминалам и персональным электровычислительным машинам и организации работ: Санитарные правила и нормы. – М.: Информационно–издательский центр Госкомсанэпиднадзора России, 1996.

22. Стандарт предприятия. Дипломная научно-исследовательская работа студента. Структура и правила оформления. СТП ЮУрГУ 19-2003 / Составители: Т.И. Парубочая, Н.В. Сырейщикова, С.Д. Ваулин, В.Р. Гофман. – Челябинск: Изд. ЮУрГУ, 2003. –19 с.

23. Сидоров А.И., Хашковский А.В., Мирзаева Н.М. Безопасность эксплуатации ЭВМ и микропроцессорной техники в составе автоматизированного производства: Учебное пособие. – Челябинск: ЧГТУ, 1990.

24. Гурский Д.А. Вычисления в MathCAD 12. – СПб.: Питер, 2006. – 544 с.

ПРИЛОЖЕНИЕ

Компьютерный диск CD-диск, содержащий:

1. Электронный вариант пояснительной записки.

2. Электронные варианты использованных источников.

3. Электронные таблицы с экспериментальными данными.