Методика аппроксимации аэродинамических характеристик баллистических ракет на активном участке траектории в зависимости от числа Маха и угла атаки

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

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

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

 (1)

где V_X, V_Y, V_Z, R_X, R_Y, R_Z, P_X, P_Y, P_Z, g_X, g_Y, g_Z – составляющие скорости, полной аэродинамической силы, силы тяги и ускорения от силы притяжения соответственно в проекциях на оси стартовой системы координат, m – масса БР,

Для оценки влияния аэродинамических коэффициентов на скорость можно воспользоваться зависимостью:

 (2)

В уравнении (2) рассмотрим пассивный (атмосферный) участок траектории. Исходя из этого масса принимается постоянной m = const. Из [1] известно, что зависимость для плотности от высоты можно расписать следующим образом:

 (3)

где  – плотность на уровне земли, H_k – текущая высота,

В свою очередь, зависимость высоты от времени можно выразить:

H_k = H_{k – 1} – V_{в(k – 1)}t, (4)

где H_{k – 1} – начальная высота БЦ в рассматриваемый момент времени t,

V_{в(k – 1)} = V_{k – 1}sin θ_вх – вертикальная скорость в рассматриваемый момент времени t,

θ_вх = 15° – угол входа в атмосферу.

При этом принимается следующее допущение: полагаем, что на протяжении 30-секундного полета на пассивном атмосферном участке скорость V_в за малый промежуток времени Δt = 1 сек и угол θ_вх на протяжении 30 сек изменяется незначительно.

С учетом (3) и (4) уравнение (2) принимает вид:

 (5)

Введем следующие обозначения:

 (6)

 (7)

Проинтегрировав уравнение (5), получаем:

 (8)

 (9)

где Δt – шаг, k ∈ N – дискретное значение времени полета, с.

Таким образом, получена зависимость текущей скорости от аэродинамического коэффициента лобового сопротивления. Зависимость (9) представляет собой рекуррентную формулу, в которой на вход в качестве скорости V_{k – 1} будет подставляться скорость V_k, полученная на предыдущем шаге, и с соответствующим пересчетом высоты аналогично формуле (4).

На рисунке 1 представлена зависимость изменения скорости БЦ от времени с шагом, равным Δt = 1, при различных значениях c_x.

Как видно из графика на рисунке 1, при изменении коэффициента лобового сопротивления на 0,1 разница в скорости на 30-й секунде на высоте приблизительно 30 км может достигать больше 200 м/с. Таким образом, для того чтобы ошибка оценки текущей скорости БЦ не превышала 10 м/с, необходимо, чтобы среднеквадратичная ошибка аппроксимации была менее 0,005.

Из уравнения (1) видно, что аэродинамические силы и моменты оказывают влияние на параметры траектории на большем участке полета.

Известны методы аппроксимации зависимости аэродинамических характеристик (АДХ) от угла атаки и числа Маха, в том числе:

− метод, основанный на представлении указанной зависимости в виде ряда Тейлора [2];

 (10)

− метод, основанный на использовании упрощенных зависимостей АДХ от угла атаки [3]:

C_x = C_x0 + C_xi(α);
C_y = C_y^α α; 
C_z = C_z^β β.                   (11)

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

Кроме того, данный метод применяется только при малых нестационарных возмущениях [2] и при малых углах атаки, когда действителен линейный закон изменения аэродинамических характеристик от кинематических параметров. Поскольку в процессе полета БР пространственный угол атаки может достигать значений 15°, то при использовании данного способа аппроксимации будут возникать большие отклонения.

На рисунке 2 приведены средние квадратические отклонения на интервале чисел Маха (0; 8) при использовании метода, основанного на использовании упрощенных зависимостей АДХ от угла атаки.

Из рисунка видно, что использование упрощенной зависимости приводит к большим ошибкам (от 15 до 30 %) на значительном интервале чисел Маха. Таким образом, ошибки оценки торможения БР могут составлять доли g, что скажется на точности моделирования траектории БР.

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

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

− непрерывность аппроксимации коэффициента C_xa в зависимости от угла атаки и числа Маха;
− погрешность аппроксимации коэффициента C_xa не должна превышать 5 %;
− параметры аппроксимационной зависимости должны определяться на основе данных о значениях коэффициента C_xa БЦ при определенных значениях угла атаки и числа Маха.

При этом дискретные данные о значениях коэффициента C_xa могут определяться с использованием существующего математического аппарата, например на основе программного комплекса SolidWorks Flow Simulation. При расчете дискретных значений коэффициента C_xa параметры атмосферы принимаются в соответствии с [4].

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

На рисунке 3 приведены дискретные зависимости коэффициента лобового сопротивления при разных углах атаки от числа Маха.

На этом же рисунке для зависимости при нулевом угле атаки обозначены характерные точки: x₀, x₁, x₂, связанные с точками изменения характера функции C_xa. Анализ рисунка 3 позволяет выявить характерные участки изгибов C_xa в зависимости от числа Маха и выбрать для них соответствующий вид аппроксимации.

Выбор вида аппроксимации осуществлялся исходя из знания характерных точек x₀, x₁, x₂, а также дополнительных условий максимума в точке x₁ и постоянства на участках x < x₀ и x > x₂.

В результате функцию f(x) зависимости коэффициента C_xa от числа Маха можно представить в виде

 (12)

где f ^* (x) – переменная часть функции f(x).

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

 (13)

Анализ показал, что в качестве аппроксимирующей функции f ^*_апп (x^*) может быть использована функция вида

 (14)

где A, B, a, b, c – параметры аппроксимации;

 (15)

где x^* – аргумент аппроксимирующей функции f^*_апп (x^*), выбираемый исходя из удобства оценки параметров аппроксимации, x^*= x – x₀, 

Для нахождения коэффициентов a, b, c необходимо решить следующую систему уравнений

 (16)

Учитывая, что полученная система уравнений является линейной относительно a, b, c, она может быть решена методом Крамера.

После проведения анализа зависимости коэффициента C_xa от числа Маха были проанализированы зависимости параметров аппроксимации от угла атаки. Параметрами аппроксимации являются характерные точки, а также соответствующие им значения функции. Значения характерных точек x₀, x₂ являются константами. Была получена следующая система уравнений (17) для определения зависимостей входных параметров аппроксимации от угла атаки. 

 (17)

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

 (18)

где M^* = M – M₀;
M₀ = x₀ – определяется выражением (17);
A [f₀(α)] = f₀(α) – 1;
B [f₀(α), f₂(α)] = f₂(α) – f₀(α) + 1;
φ [M^* , x^*₁(α)] – определяется выражением (15);
a [f₀(α), f₁(α), f₂(α), x^*₁(α), x^*₂(α)] – определяется из выражения (16);
b [f₀(α), f₁(α), f₂(α), x^*₁(α), x^*₂(α)] – определяется из выражения (16);
c [f₀(α), f₁(α), f₂(α), x^*₁(α), x^*₂(α)] – определяется из выражения (16);
x^*₁(α) = x₁(α) – x₀(α);
x^*₂(α) = x₂(α) – x₀(α);
f₀(α), f₁(α), f₂(α), x₀(α), x₁(α), x₂(α) – определяется выражением (17).

График сравнения исходной и аппроксимированной функции зависимости коэффициента лобового сопротивления от числа Маха для нулевого угла атаки, построенной с использованием функционала (18), представлен на рисунке 4.

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

За ошибку аппроксимации примем величину 

График зависимости величины ошибки Δ от числа Маха при разных углах атаки представлен на рисунке 5. 

Анализ ошибок аппроксимационных кривых на рисунке 5 показал следующее. Наибольшее расхождение между значениями коэффициентов лобового сопротивления, полученных в результате расчета в среде SolidWorks Flow Simulation, и значениями, полученными в результате аппроксимации, наблюдается на участке x₀ ≤ x ≤ x₁, что связано с резким ростом функции на коротком промежутке. Данный участок требует более тщательной проработки в дальнейшем. Величина ошибки на данном участке может достигать 29 %, на остальном участке ошибка не превышает 2,5 %.

Таким образом, анализ существующих методов расчета АДХ показал, что существующие методы аппроксимации имеют ряд недостатков, такие как трудоемкость использования для моделирования ракетно-космической обстановки в режиме реального времени, ограничения в применимости или высокая, до 30 %, ошибка. Сформирована пятипараметрическая функция аппроксимации зависимости коэффициента лобового сопротивления C_xa баллистических целей от числа Маха и угла атаки, которая позволяет на сверхзвуковых скоростях при числах Маха свыше 1,5 добиться ошибки аппроксимации не более 2,5 %, что в 10 раз точнее по сравнению с упрощенным методом. В дальнейших исследованиях целесообразно рассмотреть возможность использования предложенного подхода для других классов БР в различных вариантах их сборки, что позволит в будущем с единых позиций осуществлять моделирование АДХ.