Модель отделения створки головного обтекателя
Интегрирование уравнений движения совместно с уравнениями связей на примере модели движения створки головного обтекателя ракеты-носителя
Движение створки головного обтекателя ракеты-носителя
Головной обтекатель отделяется от ракеты-носителя после того, как ракета-носитель набрала высоту, на которой действие набегающего потока разреженного воздуха уже невелико, что позволяет сбросить “лишнюю” массу головного обтекателя, который до этого момента защищал выводимый на орбиту космический аппарат при движении ракеты-носителя в плотных слоях атмосферы. Существует несколько способов отделения головного обтекателя. Далее рассматривается способ отделения с разделением обтекателя на две створки и последующим их разворотом [2,3].
Схема системы
Для отделения головного обтекателя производится его механическое разделение на две створки, при этом каждая створка остается шарнирно-закрепленной на ракете-носителе. После этого начинается разворот створок под действием толкателей (по два на створку), каждый из которых закреплен одним концом на створке, а другим – упирается в переходный отсек ракеты-носителя. Толкатели создают момент относительно оси вращения створки, разворачивая её относительно ракеты. При достижении некоторого угла (более 45 градусов) шарнирная связь между створкой и ракетой-носителем разрывается и створка продолжает свободное движение.
Системы координат
Построим модель отделения обтекателя с момента его разделения на две створки до свободного движения створок. Будем считать, что створка – это абсолютно-твердое тело с массой значительно меньшей массы ракеты-носителя, что позволять не учитывать влияние движения створки на движение ракеты-носителя. Будет рассматриваться плоское движение створки, полагая, что все силы действуют на створку в одной плоскости, перпендикулярной оси её вращения, проходящей через центр масс створки.
Поскольку необходимо исследовать движение створки относительно ракеты, для записи уравнений движения удобней использовать неинерциальную систему координат \(O x_0 y_0 z_0\), которая связана с ракетной-носителем, движущейся с постоянным ускорением вдоль оси \(O x_0\). Также при записи уравнений движения будет использоваться система координат \(C x_1 y_1 z_1\), связанная со створкой.
Уравнения движения
Для записи уравнений систем многих тел для получения компактных уравнений движения целесообразно использовать матричную форму записи, которую потом легко перевести на языки компьютерной алгебры. Рассматриваемая механическая система простая и имеет небольшое число степеней свободы. Так при плоском движении на первом этапе движения (этапе разворота) створка имеет одну степень свободы относительно ракеты-носителя, на втором этапе – этапе свободного движения – три степени свободы. Для записи уравнений движения можно использовать уравнения Лагранжа, чтобы получить минимальное количество уравнений движения, но здесь для записи уравнений будут использованы уравнения движения центра масс и уравнения движения вокруг центра масс с учетом действия на створку силы реакции в шарнире. Уравнения движения створки будут решаться совместно с уравнениями связей. Такой подход позволит на всех этапах движения системы использовать один набор обобщенных координат: положение центра масс и поворот створки, а также определить в процессе интегрирования уравнений движения силы реакции связей в шарнире створки.
Уравнение движения центра масс створки
В системе координат \(O x_0 y_0 z_0\) движение центра масс створки описывается одним матричным уравнением, которому соответствует три скалярных уравнения \[m \ddot{\boldsymbol r}^{(0)} = \boldsymbol{P}^{(0)} + \boldsymbol \Phi^{(0)} + \boldsymbol R^{(0)}\]
где \(\boldsymbol{P}^{(0)}\) – координатный столбец вектора силы толкателя, \(\boldsymbol \Phi^{(0)}\) – координатный столбец переносной силы инерции, \(\boldsymbol R^{(0)}\) – координатный столбец силы реакции в шарнире. Верхним индексом в круглых скобках обозначается то, что координатные столбцы векторов в уравнении записаны в системе координат \(O x_0 y_0 z_0\). В неинерциальной системе координат к активным силам и силам реакции необходимо добавить переносную силу инерции, направленную вдоль продольно оси ракеты-носителя. Координатный столбец переносной силы инерции определяется следующим образом: \[\boldsymbol \Phi^{(0)} = - \begin{bmatrix} 1 \\ 0 \\ 0 \end{bmatrix} m g n_x\]
где \(m\) – масса створки, \(g\) – ускорение свободного падения, \(n_x\) – перегрузка, с которой движется ракета-носитель, которая определяет “действующую” на створку переносную силу инерции в неинерциальной системе координат \(O x_0 y_0 z_0\)
Сила толкателя, действующего на створку, пропорциональна расстоянию между точками \(P_0\) и \(P_1\) и направлена вдоль вектора \(P_0 P_1\). Обозначим это вектор буквой \(d\). В системе координат \(O x_0 y_0 z_0\) этот вектор определяется следующим образом \[\boldsymbol d^{(0)} = - r_{p0}^{(0)} + r^{(0)} + \boldsymbol{A}_1 \boldsymbol{\rho}_{p1}^{(1)}\]
Рассматривается плоское движение створки, поэтому действие двух толкателей на створку можно заменить действием одного толкателя с силой равной суммарной силе двух толкателей. Обозначим как \(\boldsymbol{e}_P\) единичный вектор, вдоль которого действует сила толкателя P \[\boldsymbol{e}_P^{(0)} = \frac{\boldsymbol d^{(0)}}{|\boldsymbol d^{(0)}|},\]
тогда координатный столбец силы, действующей на створку, в системе координат \(O x_0 y_0 z_0\) будет иметь вид \[\boldsymbol{P}^{(0)} = \boldsymbol{e}_P^{(0)} P(d),\]
где \(P(d)\) зависимость силы толкателя от его длины (перемещения штока). Сила пружинного толкателя линейно зависит перемещения штока: \[P(d) = P_0 + \frac{P_0-P_K}{h} (d-d_0)\]
где \(P_0\) - начальное усилие толкателя (двух толкателей в рассматриваемой плоской задачи), \(P_K\) - конечное усилие толкателя (это также суммарная сила от действия двух толкателей), \(d_0\) - начальная длина толкателя (расстояние от точки \(p_0\) до \(p_1\) ), \(h\) - ход толкателя.
Уравнение движения вокруг центра масс
Движение створки вокруг центра масс описывается динамическим уравнением Эйлера, которые в матричной форме в системе координат, связанной со своркой, имеют вид: \[\boldsymbol{J}^{(1)} \dot{\boldsymbol{\omega}}^{(1)} + \boldsymbol{\omega}^{(1)} \times \boldsymbol{J}^{(1)} \boldsymbol{\omega}^{(1)} = \boldsymbol{M}^{(1)}(P) + \boldsymbol{M}^{(1)}(R)\]
В правой части уравнений движения вокруг центра масс момент от силы толкателя и момент от силы реакции \(\boldsymbol{R}\). Момент от силы толкателя в системе координат, связанной со створкой записывается следующим образом: \[\boldsymbol{M}^{(1)}(P) = \tilde{\boldsymbol{r}}_{p1}^{(1)} \boldsymbol{A}_1^T \boldsymbol{P}^{(0)}\]
Здесь и далее оператор тильда обозначает кососимметричную матрицу, записанную из компонент координатного столбца вектора \[\tilde{\boldsymbol{a}} = \begin{bmatrix} 0 & - a_z & a_y \\ a_z & 0 & - a_x \\ - a_y & a_x & 0 \end{bmatrix}\]
Кососимметричная матрица используется для матричной записи векторного произведения двух координатных столбцов векторов \[\boldsymbol{a} \times \boldsymbol{b} = \tilde{\boldsymbol{a}} \boldsymbol{b} = - \tilde{\boldsymbol{b}} \boldsymbol{a}\]
Момент от силы реакции относительно центра масс створки: \[\boldsymbol{M}^{(1)}(R) = \tilde{\boldsymbol{c}}_{11}^{(1)} \boldsymbol{A}_1^T \boldsymbol{R}^{(0)}\]
С учётом полученных выражений для моментов, уравнение движения вокруг центра масс принимает вид \[\boldsymbol{J}^{(1)} \dot{\boldsymbol{\omega}}^{(1)} = \tilde{\boldsymbol{r}}_{p1}^{(1)} \boldsymbol{A}_1^T \boldsymbol{P}^{(0)} + \tilde{\boldsymbol{c}}_{11}^{(1)} \boldsymbol{A}_1^T \boldsymbol{R}^{(0)} - \boldsymbol{\omega}^{(1)} \times \boldsymbol{J}^{(1)} \boldsymbol{\omega}^{(1)}\]
В уравнения движения центра масс и движения вокруг центра масс входят неизвестные силы реакции, которые можно определить, дополнив полученные уравнения движения уравнениями связей.
Уравнения связей
Положение шарнирной точки, связанной со створкой, в системе координат \(O x_0 y_0 z_0\) определяется суммой вектора положения центра масс створки и положения шарнирной точки створки относительно центра масс \(\mathbf{r}_1^{(0)} + \mathbf{A}_1 \mathbf{c}_{11}^{(1)}\). Матрица поворота или матрица преобразования координат используется в этом выражении для преобразования координат вектора \(\mathbf{c}_{11}^{(1)}\) из системы координат, связанной со створкой в систему координат \(O x_0 y_0 z_0\). Система координат, связанная со створкой \(C_1 x_1 y_1 z_1\), повёрнута относительно системы координат \(O x_0 y_0 z_0\) на угол \(\varphi\) вокруг оси \(z_0\), поэтому матрица поворота имеет простой вид \[\mathbf{A}_1 = \begin{bmatrix} \cos \varphi & -\sin \varphi & 0 \\ \sin \varphi & \cos \varphi & 0 \\ 0 & 0 & 1 \end{bmatrix}\]
Шарнирная точка створки должна совпадать с соответствующей шарнирной точкой на ракете-носителе, положение которой относительно \(O x_0 y_0 z_0\) определяется неизменным координатным столбцом \(\mathbf{c}_{01}^{(0)}\), поэтому \[\mathbf{c}_{01}^{(0)} = \mathbf{r}_1^{(0)} + \mathbf{A}_1 \mathbf{c}_{11}^{(1)}\]
Полученное уравнение представляет собой уравнение связи, которое говорит о том, что створка и ракета-носитель в процессе движения имеют одну общую точку – шарнирную точку. При развороте створки относительно шарнира ограничивается и угловое движение створки как твёрдого тела относительно двух других осей, перпендикулярных оси шарнира, однако если на створку действует плоская система сил, лежащая в плоскости \(O x_0 y_0 z_0\), в которой также находится и центр масс створки, то дополнительные уравнения связей можно не записывать.
Для интегрирования дифференциальных уравнений движения совместно с алгебраическими уравнениями связей, обычно используют следующую технику решения, которая позволят использовать для интегрирования уравнений движения “стандартные” численные методы решения обыкновенных дифференциальных уравнений. Уравнения связей дважды дифференцируют для получения уравнений связей на ускорения. Учитывая постоянство координатных столбцов \(\mathbf{c}_{11}^{(1)}\) и \(\mathbf{c}_{01}^{(0)}\), после первого дифференцирования уравнения связи получим: \[\mathbf{0} = \mathbf{v}_1^{(0)} + \dot{\mathbf{A}}_1 \mathbf{c}_{11}^{(1)}\]
Производная матрицы поворота определяется выражением [1]: \[\dot{\mathbf{A}}_1 = \mathbf{A}_1 \tilde{\mathbf{\omega}}_1^{(1)}\]
где \(\tilde{\mathbf{\omega}}_1^{(1)}\) – кососимметрическая матрица, составленная из компонент вектора угловой скорости створки в проекциях на оси системы координат, связанной со створкой \[\tilde{\mathbf{\omega}}_1^{(1)} = \begin{bmatrix} 0 & - \omega_z^{(1)} & \omega_y^{(1)} \\ \omega_z^{(1)} & 0 & - \omega_x^{(1)} \\ - \omega_y^{(1)} & \omega_x^{(1)} & 0 \end{bmatrix}\]
С учётом последнего выражения, матричное уравнение связи примет вид \[\mathbf{0} = \mathbf{v}_1^{(0)} + \mathbf{A}_1 \tilde{\mathbf{\omega}}_1^{(1)} \mathbf{c}_{11}^{(1)}\]
Это уравнение, в отличие исходного уравнения связи, уже говорит о том, что не положения, а скорости общих точек створки и ракеты-носителя должны совпадать. При условии того, что начальные условия будет удовлетворять исходному уравнению связи и при точном интегрировании уравнений движения исходное уравнение связи тоже будет выполняться, однако вследствие неизбежной погрешности численного интегрирования возможен так называемый “дрейф” связи – нарушение исходного уравнения связи и удаление шарнирных на створке и ракете-носителе друг от друга. В рассматриваемой задаче, учитывая небольшой интервал времени, на котором рассматривается движение системы, эта погрешность будет мала при разумном выборе точности интегрирования системы дифференциальных уравнений.
После второго дифференцирования уравнение связи принимает вид \[\mathbf{0} = \mathbf{a}_1^{(0)} + \mathbf{A}_1 \tilde{\boldsymbol{\omega}}_1^{(1)} \tilde{\boldsymbol{\omega}}_1^{(1)} \mathbf{c}_{11}^{(1)} + \mathbf{A}_1 \tilde{\boldsymbol{\varepsilon}}_1^{(1)} \mathbf{c}_{11}^{(1)}\]
Полученное уравнение связи требует, чтобы ускорения общих точек створки и ракеты-носителя были одинаковы в процессе интегрирования уравнений движения. Начальные условия должны удовлетворять исходному уравнению связи и уравнению связи, записанному для скоростей точек.
Перепишем уравнение связи, выделив матрицу коэффициентов перед ускорениями, используя тождество \(\tilde{\mathbf a} \mathbf b = - \tilde{\mathbf b} \mathbf a\) \[\mathbf{0} = \mathbf{a}_1^{(0)} + \mathbf{A}_1 \tilde{\boldsymbol{\omega}}_1^{(1)} \tilde{\boldsymbol{\omega}}_1^{(1)} \mathbf{c}_{11}^{(1)} - \mathbf{A}_1 \tilde{\mathbf{c}}_{11}^{(1)} \boldsymbol{\varepsilon}_1^{(1)}\]
Собрав из линейных и угловых ускорений створки столбец \(6 \times 1\), запишем это уравнение в блочном матричном виде \[\begin{bmatrix} \mathbf{E}_{3 \times 3} & - \mathbf{A}_1 \tilde{\mathbf{c}}_{11}^{(1)} \end{bmatrix} \begin{bmatrix} \mathbf{a}_1^{(0)} \\ \boldsymbol{\varepsilon}_1^{(1)} \end{bmatrix} = - \mathbf{A}_1 \tilde{\mathbf{\omega}}_1^{(1)} \tilde{\boldsymbol{\omega}}_1^{(1)} \mathbf{c}_{11}^{(1)}\]
Система дифференциально-алгебраических уравнений в матричном виде
Этап 1
Объединив уравнения движения и уравнения связи в одно блочное матричное уравнение, получим для первого этапа движения (вращение створки) \[\begin{bmatrix} \mathbf{E}_{3 \times 3} & \mathbf{0}_{3 \times 3} & - \mathbf{E}_{3 \times 3} \\ \mathbf{0}_{3 \times 3} & \mathbf{J}_1 & - \tilde{\mathbf{c}}_{11}^{(1)} \mathbf{A}_1^T \\ \mathbf{E}_{3 \times 3} & - \mathbf{A}_1 \tilde{\mathbf{c}}_{11}^{(1)} & \mathbf{0}_{3 \times 3} \end{bmatrix} \begin{bmatrix} \mathbf{a}_1^{(0)} \\ \boldsymbol{\varepsilon}_1^{(1)} \\ \mathbf{R}^{(0)} \end{bmatrix} = \begin{bmatrix} \boldsymbol{P}^{(0)} + \boldsymbol{F}^{(0)} \\ \tilde{\boldsymbol{r}}_{p1}^{(1)} \boldsymbol{A}_1^T \boldsymbol{P}^{(0)} \\ - \mathbf{A}_1 \tilde{\mathbf{\omega}}_1^{(1)} \tilde{\mathbf{\omega}}_1^{(1)} \mathbf{c}_{11}^{(1)} \end{bmatrix}\]
Поскольку \[(\mathbf{A}_1 \tilde{\mathbf{c}}_{11}^{(1)})^T = - \tilde{\mathbf{c}}_{11}^{(1)} \mathbf{A}_1^T,\]
с точностью до знака полученная матрица коэффициентов является симметричной. Умножив уравнение связи на минус 1, получим симметричную матрицу коэффициентов \[\begin{bmatrix} \mathbf{E}_{3 \times 3} & \mathbf{0}_{3 \times 3} & - \mathbf{E}_{3 \times 3} \\ \mathbf{0}_{3 \times 3} & \mathbf{J}_1 & - \tilde{\mathbf{c}}_{11}^{(1)} \mathbf{A}_1^T \\ - \mathbf{E}_{3 \times 3} & \mathbf{A}_1 \tilde{\mathbf{c}}_{11}^{(1)} & \mathbf{0}_{3 \times 3} \end{bmatrix} \begin{bmatrix} \mathbf{a}_1^{(0)} \\ \boldsymbol{\varepsilon}_1^{(1)} \\ \mathbf{R}^{(0)} \end{bmatrix} = \begin{bmatrix} \boldsymbol{P}^{(0)} + \boldsymbol{F}^{(0)} \\ \tilde{\boldsymbol{r}}_{p1}^{(1)} \boldsymbol{A}_1^T \boldsymbol{P}^{(0)} \\ \mathbf{A}_1 \tilde{\mathbf{\omega}}_1^{(1)} \tilde{\mathbf{\omega}}_1^{(1)} \mathbf{c}_{11}^{(1)} \end{bmatrix}\]
Структура этой системы дифференциально-алгебраических уравнений следующая \[\begin{bmatrix} \boldsymbol{M} & \boldsymbol{Q}^T \\ \boldsymbol{Q}^T & \boldsymbol{0} \end{bmatrix} \begin{bmatrix} \ddot{\boldsymbol{X}} \\ \boldsymbol{R} \end{bmatrix}= \begin{bmatrix} \boldsymbol{L} \\ \boldsymbol{B} \end{bmatrix}\]
Эту систему можно представить в виде двух матричных уравнений \[\begin{aligned} \boldsymbol{M} \ddot{\boldsymbol{X}} + \boldsymbol{Q}^T \boldsymbol{R} = \boldsymbol{L} \\ \boldsymbol{Q}^T \ddot{\boldsymbol{X}} = \boldsymbol{B} \end{aligned}\]
Выразив из первого уравнения \(\ddot{\boldsymbol{X}}\) \[\ddot{\boldsymbol{X}} = \boldsymbol{M}^{-1} ( \boldsymbol{L} - \boldsymbol{Q}^T \boldsymbol{R} )\]
и подставив во второе уравнение, получим систему линейный уравнений для определения сил реакций \[(\boldsymbol{Q}^T \boldsymbol{M}^{-1} \boldsymbol{Q}^T) \boldsymbol{R} = \boldsymbol{Q}^T \boldsymbol{M}^{-1} \boldsymbol{L} - \boldsymbol{B},\]
что позволяет найти столбец ускорений \[\ddot{\boldsymbol{X}} = \boldsymbol{M}^{-1} (\boldsymbol{L} - \boldsymbol{Q}^T (\boldsymbol{Q}^T \boldsymbol{M}^{-1} \boldsymbol{Q}^T)^{-1} (\boldsymbol{Q}^T \boldsymbol{M}^{-1} \boldsymbol{L} - \boldsymbol{B}) )\]
Этап 2
На втором этапе движения механическая связь между створкой и ракетной-носителем отсутствует. Для учёта этого уравнение связи можно заменить на матричное уравнение (три скалярных уравнения) \[\mathbf{R}^{(0)} = 0\]
тогда, система уравнений примет вид \[\begin{bmatrix} \mathbf{E}_{3 \times 3} & \mathbf{0}_{3 \times 3} & - \mathbf{E}_{3 \times 3} \\ \mathbf{0}_{3 \times 3} & \mathbf{J}_1 & - \tilde{\mathbf{c}}_{11}^{(1)} \mathbf{A}_1^T \\ \mathbf{0}_{3 \times 3} & \mathbf{0}_{3 \times 3} & \mathbf{E}_{3 \times 3} \end{bmatrix} \begin{bmatrix} \mathbf{a}_1^{(0)} \\ \boldsymbol{\varepsilon}_1^{(1)} \\ \mathbf{R}^{(0)} \end{bmatrix} = \begin{bmatrix} \boldsymbol{P}^{(0)} + \boldsymbol{F}^{(0)} \\ \tilde{\boldsymbol{r}}_{p1}^{(1)} \boldsymbol{A}_1^T \boldsymbol{P}^{(0)} \\ - \mathbf{A}_1 \tilde{\mathbf{\omega}}_1^{(1)} \tilde{\mathbf{\omega}}_1^{(1)} \mathbf{c}_{11}^{(1)} \end{bmatrix}\]
Код
Примеры кода на языках MATLAB и Python, реализующие описанный алгоритм.
Код на языке MATLAB
Файл-скрипт main.m
clc; clear all;
% Положение шарнирной точки отчносительно центра масс створки с ССК
p.c11 = [-4.0; 1.0; 0.0];
% Положение точки закрепления толкателя на створке относительно центра масс
% створки в ССК
p.rp1 = [-3.5; -0.8; 0.0];
% Положение точки закрепления толкателя на РН относительно БСК
p.rp0 = [ 0.0; 0.2; 0.0];
% Ускорение свободного падения
p.g = 9.807;
% Начальное усиле толкателей
p.P0 = 30000;
% Конечное усилие толкателей
p.Pk = 10000;
% Ход толкателя
p.h = 0.4;
% Начальная длина толкателя
p.d0 = 0.5;
% Масса створки
p.m = 1000;
% Тензор инерции створки
p.J = diag([10000,10000,10000]);
%
p.phimax = 60;
% Начальные условия
q0 = [4.0; 1.0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0];
% Первый этап
p.Stage = 1;
opt = odeset('Events', @(t,q) PhiEvent(t,q,p));
[t1,q1] = ode45(@(t, q) dqdt(t, q, p), 0:0.02:2, q0, opt);
% Второй этап
p.Stage = 2;
[t2,q2] = ode45(@(t, q) dqdt(t, q, p), [t1(end):0.02:3], q1(end, :)');
% Склейка результатов двух этапов
t = [t1(1:end - 1); t2];
q = [q1(1:end - 1, :); q2];
Файл-функция правых частей дифференциальных уравнений dqdt.m
function [dq, R]= dqdt(t, q, p)
% Извлекаем компоненты из вектора состояния
% Координатный столбцец центра масс
r = q(1:3);
% Углы поворота створки
% в плоской задаче испольуется только последний, а остальные принимаются
% равными нулю
phi = q(4:6);
% Координатный столбец скорости центра масс створки в СК РН
v = q(7:9);
% Координатный столбец угловой скорости
w = q(10:12);
% в плоской задаче испольуется только последнее значение - wz, а остальные принимаются
% равными нулю
wz = w(3);
% Этап движения
Stage = p.Stage;
% Матрица поворота вокруг оси z
A = @(a) [cos(a), -sin(a) 0; sin(a), cos(a), 0;0 0 1];
% Оператор тильда
tilde = @(a) [0 -a(3) a(2); a(3) 0 -a(1); -a(2) a(1) 0];
% Матрица поворота на угол phi вокруг оси z
A1 = A(phi(3));
% Матрица коэффиицентов при ускорениях в уравнениях связей, ограичивающих
% перемещение шарнирной точки
Q = - A1*tilde(p.c11);
% Матрица коэффициентов при реакциях в уравнении движения створки вокруг
% центра масс
QT = - tilde(p.c11)*A1';
% Вектор от точки закрепления толкателя на РН до точки закрепления
% толкателя на створке
d = - p.rp0 + r + A1*p.rp1;
% Модуль этого вектора
nd = norm(d);
% Единичный вектор этого напарвления
ed = d/nd;
% Если длина толкателя не больше его начальной длины и хода
% вычисляем силу
if nd < p.d0 + p.h
Pin0 = ed*(p.P0 - (p.P0-p.Pk)*(nd-p.d0)/p.h );
% иначе, толкатель уже не работает
else
Pin0 = [0;0;0];
end
% Координатный столбец силы в ССК створки
Pin1 = A1'*Pin0;
% Координатный столбец силы инерции в СК РН (вдоль продольной оси)
InertialForce = [-1;0;0]*p.m*p.g;
E3 = eye(3);
Z3 = zeros(3);
if Stage == 1
% Матрица коэффициентов для первого этапа
M = [p.m*E3 Z3 -E3;
Z3 p.J -QT;
E3 Q Z3];
end
if Stage == 2
% Матрица коэффициентов для второго этапа
% Реакции равны нулю
M = [p.m*E3 Z3 E3;
Z3 p.J -QT;
Z3 Z3 E3];
end
% Матрица-столбец правой части ДАУ
B = [Pin0 + InertialForce;
tilde(p.rp1)*Pin1;
-A1*tilde(w)*tilde(w)*p.c11];
% Находим ускорения и реакции
X = M\B;
R = X(7:9);
% Кинематические уравнения для плоского движения
% Производня последнего угла равна угловой скорости, остальные компоненты
% равны нулю
dphi = [0;0;wz];
% Производня вектора состояния
dq = [v; dphi; X(1:6)];
end
Файл-функция событий системы дифференциальных уравнений dqdt.m
Файл-функция событий передается в качестве параметра в функцию интегрирования дифференциальных уравнений и предназначена для остановки процесса интегрировании при заданном условии, которое задается функцией, изменяющей знак. В рассматриваемой задаче такой функцией будет угол поворота створки.
function [Value, isterminal, direction] = PhiEvent(t, q,p)
% Если угол поврота превышает phimax...
Value = q(6)-deg2rad(p.phimax);
% при увеличении угла,
direction = 1;
% то остановить интегрирование
isterminal = 1;
end
Код на языке Python
Вес код программы можно разместить в одном файле. Удобно для решения подобных задач использовать Jupyter notebook
# Массивы, матрицы
import numpy as np
# Для численного интегрирования дифференциальных уранвений
from scipy.integrate import solve_ivp
# Для построения графиков
import matplotlib.pyplot as plt
import matplotlib as mpl
# Настройки графиков по умолчанию
import matplotlib.pylab as pylab
params = {'legend.fontsize': 14, 'figure.figsize': (10, 7), 'axes.labelsize': 14,
'axes.titlesize':14, 'xtick.labelsize':14,'ytick.labelsize':14}
pylab.rcParams.update(params)
# Для передачи параметров в функцию объеденим их в структуру namedtuple
from collections import namedtuple
# Параметры системы
Parameters = namedtuple('Parameters' , 'm J c11 c01 rp0 rp1 d0 h P0 PK nx g phimax')
param1 = Parameters(m=1000.0,
J=np.diag([5000,10000,10000.]),
c11 = np.array([-4.0, 1.0, 0.0]),
c01 = np.array([ 0.0, 2.0, 0.0]),
rp0 = np.array([ 0.0, 0.4, 0.0]),
rp1 = np.array([-3.5,-0.6, 0.0]),
d0 = 0.5,
h = 0.4,
P0 = 50000,
PK = 40000,
nx = 1.0,
g = 9.807,
phimax = np.deg2rad(50))
Функции матрицы поворота и оператор тильда
# Матрица 0 3x3
Z3 = np.zeros((3, 3))
# Единичная матрица 3x3
E3 = np.identity(3)
# Матрица поворота
def Az(a):
return np.array([[np.cos(a), -np.sin(a), 0],[np.sin(a), np.cos(a), 0],[0, 0, 1.0]]);
# Тильда
def tilde(a):
return np.array([[ 0, -a[2], a[1] ],
[ a[2], 0, -a[0] ],
[ -a[1], a[0], 0 ] ])
Функция правых частей системы дифференциальных уравнений.
def dqdt(t,q,p):
# координатный столбец центра масс
r = q[0:3]
# углы поворота створки
a = q[3:6]
# в рассматриваемой плоской задаче
# используется только один - последний из трёх
phi = a[2]
# координатный столбец скорости центра масс
v = q[6:9]
# координатный столбец угловой скорости створки
w = q[9:12]
# Вычисляем матрицу поворота для текущего значения угла
A1 = Az(phi)
# Коэффициенты при ускорениях в уравнениях связей
Q = np.dot(A1, tilde(p.c11))
# Переносная сила инерции
F = np.array([-1,0,0])*p.m*p.g*p.nx
# Вектор p0p1
d = -p.rp0 + r + np.dot(A1,p.rp1)
# Длина вектора
dn = np.sqrt(np.dot(d,d))
# Едичиный вектор направления силы в базисе x0y0z0
eP = d/dn
# Сила толкателя
P = 0*eP
if dn<p.d0+p.h:
P = (p.P0 + (dn-p.d0)*(p.P0-p.PK)/p.h)*eP
# Матрица коэффициентов
A = np.block([[E3*p.m, Z3, -E3],
[Z3 , p.J, np.transpose(Q)],
[-E3 , Q, Z3]])
if STAGE == 2:
# Матрица коэффициентов
A = np.block([[E3*p.m, Z3, -E3],
[Z3 , p.J, np.transpose(Q)],
[Z3 , Z3, E3]])
# Матрица правой части
B = np.block([ F + P,
-np.cross(w,np.dot(p.J,w)) + np.cross(p.rp1, np.dot(np.transpose(A1),P)),
np.dot(np.dot(A1,np.dot(tilde(w),tilde(w))),p.c11)])
X = np.linalg.solve(A,B)
dq = np.block([v,w,X[0:6]])
return dq
Функция-“детектор”, передаваемая в интегратор (параметр events), для определения времени перехода от этапа 1 к этапу 2 и остановки процесса интегрирования функция возвращает разницу между углом потери связи и текущим углом
def event_phi_eq_phimax(t, q):
return param1.phimax-q[5]
# функция-детектор останавливает процесс интегрирование
event_phi_eq_phimax.terminal = True
# при пересечении функции нуля при убывании
event_phi_eq_phimax.direction = -1.0
Запускаем процесс интегрирования
q0 = [4, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0];
STAGE = 1
s1 = solve_ivp(lambda t, y: dqdt(t, y, param1), [0, 1.5], q0,
events = [event_phi_eq_phimax], rtol = 1e-7, max_step=0.1)
STAGE = 2
q1 = s1.y[:,-1]
s2 = solve_ivp(lambda t, y: dqdt(t, y, param1), [s1.t[-1], 2], q1, rtol = 1e-6, max_step=0.1)
График изменения угла и угловой скорости створки
plt.plot(s1.t,s1.y[5]*180/np.pi,'r-')
plt.plot(s2.t,s2.y[5]*180/np.pi,'r--')
plt.plot(s1.t,s1.y[11]*180/np.pi,'b-')
plt.plot(s2.t,s2.y[11]*180/np.pi,'b--')
plt.xlabel('t, c'), plt.ylabel('$\\varphi$, $\\omega$');
plt.grid(ls='--'), plt.legend(['$\\varphi$ Этап 1','$\\varphi$ Этап 2']);
Источники
- Маркеев А. П. Теоретическая механика: Учебник для университетов. 3-е изд. — М.; Ижевск: РХД, 2007. — 592 с.
- Колесников К.С., Кокушкин В.В., Борзых С.В., Панкова Н.В. Расчет и проектирование систем разделения ступеней ракет: учебное пособие. М.: МГТУ им. Н.Э. Баумана. 2006. 376 с.
- Круглов Г. Е. Аналитическое проектирование механических систем [Электронный ресурс] : учеб. пособие. Самарский государственный аэрокосмический университет им. С. П. Королева. 2000.
