Материалы семинара "Методы построения математических моделей космических систем с подвижными элементами"

Авторы: 
Юдинцев В. В.
Тип материала: 
электронное издание
Год опубликования: 
2012
Издательство, журнал, сборник.: 
Материалы семинара кафедры теоретической механики

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

  • системы раскрытия солнечных батарей, антенн, радиаторов;
  • мистемы отделения ступеней отработавших блоков ракет;
  • роботы-манипуляторы;
  • наземные экспериментальные установки;

Анализ таких систем обычно невозможен без использования компьютера: уравнения движения часто настолько сложны многомерны, что найти аналитическое решение невозможно. Этим объясняется и то, что методы формирования уравнений движения этих систем стали разрабатываться с 60х годов 20 века, когда с развитием техники с выходом в космос возникла необходимость построения моделей сложных многотельных механизмов , а также появлилась возможность доведения результата моделирования "до числа". Уже тогда выяснилось, что традиционные методы вывода уравнений приводят к громоздким и неэфективным программам, требуют проведения большой предварительной работы, большого обеъма аналитических операций, труднореализуемых и неэффективных при машинной реализации. Были разработаны методы построения уравнений ориентированные дальнейшее "машинное" численное интегрирование.

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


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

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

Использование абсолютных координат


Разделим маятник на отдельные тела, заменив шарниры реакциями связей и запишем уравненя движения каждого тела в отдельности:
\begin{equation}
\left\{
\begin{aligned}
m_1 \ddot x_1 & =F_{1x} - R_{10x}+{R_{12x}},\\
m_1 \ddot y_1 & =F_{1y}-{R_{10y}}+{R_{12y}},\\
J_{1z} \ddot \varphi_1 & =M_{1z}+ l_1 ({R_{10x}} \cos \varphi_1 + \\
& + {R_{10y}} \sin \varphi_1 + \\
& + {R_{12x}} \cos \varphi_1 +{R_{12y}} \sin \varphi_1) ,\\
m_2 \ddot x_2 & =F_{2x}-{R_{21x}},\\
m_2 \ddot y_2 & =F_{2y}-{R_{21y}},\\
J_{2z} \ddot \varphi_2 & =M_{2z}+{R_{12x}} l_2 \cos \varphi_2 + \\
& +{R_{12y}} l_2 \sin \varphi_2\\
\end{aligned}\right.
\end{equation}

Эта система 6 дифференциальных уравнений содержит 10 неизвестных функций времени: 6 ускорений и 4 реакции связи. Для того, чтобы можно было проинтегрировать эту систему, её необходимо дополнить уравнениями связи, которые "говорят" о том, что кинематические параметры тел не являются независимыми. Система уравнений связи имеет следующий вид:

\begin{equation}
\left\{
\begin{aligned}
\boxed{x_{A_1} = x_{A_0}} = x_1 - l_1 \sin \varphi_1 & = 0 \\
\boxed{y_{A_1} = y_{A_0}} = y_1 + l_1 \cos \varphi_1 & = 0 \\
\boxed{x_{B_1} = x_{B_2}} = x_1 + l_1 \sin \varphi_1 & = \\
= x_2 - l_2 \sin \varphi_2 & \\
\boxed{y_{B_1} = y_{B_2}} = y_1 - l_1 \cos \varphi_1 & = \\
= y_2 - l_2 \cos \varphi_2 & \\
\end{aligned}\right.
\end{equation}

Чтобы решить полученную систему дифференциальных уравнений

\begin{equation}
\left\{
\begin{aligned}
& m_1 \ddot x_1 =F_{1x}-{R_{10x}}+{R_{12x}},\\
& m_1 \ddot y_1 =F_{1y}-{R_{10y}}+{R_{12y}},\\
& J_{1z} \ddot \varphi_1 = M_{1z}+ l_1 ({R_{10x}} \cos \varphi_1 +{R_{10y}} \sin \varphi_1 +
+ {R_{12x}} \cos \varphi_1 +{R_{12y}} \sin \varphi_1) ,\\
& m_2 \ddot x_2 = F_{2x}-{R_{21x}},\\
& m_2 \ddot y_2 = F_{2y}-{R_{21y}},\\
& J_{2z} \ddot \varphi_2 = M_{2z}+{R_{12x}} l_2 \cos \varphi_2 +{R_{12y}} l_2 \sin \varphi_2\\
& x_1 - l_1 \sin \varphi_1 = 0 \\
& y_1 + l_1 \cos \varphi_1 = 0 \\
& x_1 + l_1 \sin \varphi_1 = x_2 - l_2 \sin \varphi_2 \\
& y_1 - l_1 \cos \varphi_1 = y_2 - l_2 \cos \varphi_2 \\
\end{aligned} \right.
\end{equation}

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

$$
\left\{
\begin{aligned}
& m_1 \ddot x_1 =F_{1x}-{R_{10x}}+{R_{12x}},\\
& m_1 \ddot y_1 =F_{1y}-{R_{10y}}+{R_{12y}},\\
& J_{1z} \ddot \varphi_1 = M_{1z}+ l_1 ({R_{10x}} \cos \varphi_1 + {R_{10y}} \sin \varphi_1 +
+{R_{12x}} \cos \varphi_1 +{R_{12y}} \sin \varphi_1) ,\\
& m_2 \ddot x_2 = F_{2x}-{R_{21x}},\\
& m_2 \ddot y_2 = F_{2y}-{R_{21y}},\\
& J_{2z} \ddot \varphi_2 = M_{2z}+{R_{12x}} l_2 \cos \varphi_2 +{R_{12y}} l_2 \sin \varphi_2\\
& \ddot x_1 - \ddot \varphi_1 l_1 \cos \varphi_1 = \dot \varphi_1^2 l_1 \sin \varphi_1\\
& \ddot y_1 - \ddot \varphi_1 l_1 \sin \varphi_1 = \dot \varphi_1^2 l_1 \cos \varphi_1 \\
& \ddot x_1 + \ddot \varphi_1 \cos \varphi_1 - \ddot x_2 + \ddot \varphi_2 l_2 \cos \varphi_2 = \dot \varphi_1^2 l_1 \sin \varphi_1 + \dot \varphi_2^2 l_2 \sin \varphi_2\\
& \ddot y_1 + \ddot \varphi_1 l_1 \sin \varphi_1 - \ddot y_2 - \ddot \varphi_2 l_2 \sin \varphi_2 = \dot \varphi_2^2 l_2 \cos \varphi_2 - \dot \varphi_1^2 l_1 \cos \varphi_1 \\
\end{aligned} \right.
$$

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

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

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


Первое уравнение связи "точка-плоскость"
$$
\boxed{(\pmb \rho_i - \pmb \rho_P) \cdot \mathbf n_i = 0}
$$
ограничивает относительное движение смежных тел так, что неизменная точка одного тела $j$ скользит по плоскости, жестко связанной с другим телом $i$.

Уравнение, ограничивающее относительное вращение требует, чтобы вектор относительной угловой скорости двух тел не имел проекции на заданное направление:
$$
\boxed{ (\pmb \omega_i - \pmb \omega_j) \cdot \mathbf n_i = 0 }
$$

Для уравнения связи точка плоскость:
\begin{equation*}\label{eq:point-plate1}
\boxed{\mathbf n_{i}^{(i)T} \ddot{\pmb{\rho}}_{i}^{(i)} =0}
\end{equation*}
Координаты вектора точки контакта $\pmb \rho_{i}$ в базисах $i$ и $0$
$$
\pmb \rho_{i}^{(i)} = \mathbf{A}^{iT} \pmb \rho_{i}^{(0)}
$$
$$
\pmb \rho _{i}^{(0)}= \mathbf{A}^{j} \pmb \rho _{j}^{(j)} + \mathbf r_j^{(0)} - \mathbf r_i^{(0)}
$$
Скорость точки контакта в базисе $i$
$$
\dot{\pmb{\rho}}_{i}^{(i)} = \dot{\mathbf{A}}^{iT} \pmb \rho_{i}^{(0)} + \mathbf{A}^{iT} \dot{\pmb \rho_{i}}^{(0)}
$$
Производные матрицы поворота
$$\dot{\mathbf{A}}^{iT}=-\tilde{\pmb \omega}_i^{(i)} \mathbf{A}^{iT}, \ \dot{\mathbf{A}}^{i}= \mathbf{A}^{i} \tilde{\pmb \omega}_i^{(i)}
$$
Скорость точки контакта в проекциях на оси инерциального базиса
$$
\dot{\pmb \rho}_{i}^{(0)} = \mathbf{A}^j \tilde{\pmb \omega}_j^{(j)} \pmb \rho_{j}^{(j)}+\dot{\mathbf{r}}_j^{(0)} - \dot{\mathbf{r}}_i^{(0)}
$$
Ускорение точки контакта
$$
\begin{aligned}
\ddot{\pmb \rho}_{i}^{(i)} & = -\dot{\tilde{\pmb{\omega}}}_i^{(i)} \mathbf A^{iT} \pmb{\rho}_{i}^{(0)} + \tilde{\pmb \omega}_i^{(i)} \tilde{\pmb \omega}_i^{(i)} \mathbf{A}^{rT}_i \pmb{\rho}_{i}^{(0)} - \tilde{\pmb \omega}_i^{(i)} \mathbf A^{iT} \dot{\pmb \rho}_{i}^{(0)} - \tilde{\pmb \omega}_i^{(i)} \mathbf{A}^{rT}_i \dot{\mathbf \rho}_{i}^{(0)} + \\ & +\mathbf A^{iT} \ddot{\mathbf \rho}_{i}^{(0)} \\
\ddot{\pmb \rho}_{i}^{(0)} & = \dot{\mathbf{A}}^j \tilde{\pmb \omega}_j^{(j)} {\pmb \rho}_{j}^{(j)}+\mathbf{A}^j \dot{\tilde{\pmb \omega}}_j^{(j)} \rho_{j}^{(j)}+\ddot{\mathbf{r}}_j^{(0)} - \ddot{\mathbf{r}}_i^{(0)}= \\
& = \mathbf{A}^j \tilde{\pmb \omega}_j \tilde{\pmb \omega}_j^{(j)} {\pmb \rho}_{j}^{(j)}+\mathbf{A}^j \dot{\tilde{\pmb \omega}}_j^{(j)} \pmb \rho_{j}^{(j)}+\ddot{\mathbf{r}}_j^{(0)} - \ddot{\mathbf{r}}_i^{(0)} \\
\ddot{\pmb \rho}_{i}^{(i)} & = {\tilde{\pmb \rho}_{i}^{(i)}} \dot{\pmb {\omega}}_i^{(i)} - {\mathbf A^{iT} \mathbf{A}^j \tilde{\pmb \rho}_{j}^{(j)}} \dot{\pmb \omega}_j^{(j)} + {\mathbf A^{iT}} \ddot{\mathbf{r}}_j^{(0)}-{\mathbf A^{iT}} \ddot{\mathbf{r}}_i^{(0)} + \tilde{\pmb \omega}_i^{(i)} \tilde{\pmb \omega}_i^{(i)} \mathbf{A}^{iT} \pmb{\rho}_{i}^{(0)}-\\
& - 2 \tilde{\pmb \omega}_i^{(i)} \mathbf A^{iT} \dot{\pmb \rho}_{i}^{(0)} + \mathbf A^{iT} \mathbf{A}^j \tilde{\pmb \omega}_j^{(j)} \tilde{\pmb \omega}_j^{(j)} \pmb \rho_{j}^{(j)}
\end{aligned}
$$
Матричная форма уравнения связи "точка-плоскость":
\begin{equation*}
\boxed{\mathbf n_{i}^{(i)T} \ddot{\pmb{\rho}}_{i}^{(i)} =0}
\end{equation*}
$$
\mathbf Q_i {\mathbf{X}}_i+ \mathbf Q_j {\mathbf{X}}_j=\mathbf{b}_{ij},
$$
$$
{\mathbf{X}}_i=
\begin{bmatrix}
\ddot{\mathbf r}_i^{(0)} \\
\dot{\pmb{\omega}}_i^{(i)}
\end{bmatrix}, \
{\mathbf{X}}_j=
\begin{bmatrix}
\ddot{\mathbf r}_j^{(0)} \\
\dot{\pmb{\omega}}_j^{(j)}
\end{bmatrix};
$$
$$
\mathbf Q_i = \begin{bmatrix} - \mathbf n_{i}^{(i)T} \mathbf{A}^{iT} & \mathbf n_{i}^{(i)T} \tilde{\pmb{\rho}}_{i}^{(i)} \end{bmatrix},
\
\mathbf Q_j = \begin{bmatrix} \mathbf n_{i}^{(i)T} \mathbf{A}^{iT} & - \mathbf n_{i}^{(i)T}\mathbf A^{iT} \mathbf{A}^j \tilde{\pmb \rho}_{j}^{(j)} \end{bmatrix},
$$
$$
\mathbf{b}_{i}=\mathbf n_{i}^{(i)T}\left(2 \tilde{\pmb \omega}_i^{(i)} \mathbf A^{iT} \dot{\pmb \rho}_{i}^{(0)}-\tilde{\pmb \omega}_i^{(i)} \tilde{\pmb \omega}_i^{(i)} \mathbf{A}^{iT} \pmb \rho_{i} - \mathbf A^{iT} \mathbf{A}^j \tilde{\pmb \omega}_j^{(j)} \tilde{\pmb \omega}_j^{(j)} \pmb \rho_{j}^{(j)}\right).
$$

Для уравнения связи, ограничивающего относительное вращение
$$
\boxed{\mathbf{n}_{i}^{(i)T} \pmb \varepsilon_{ij}^{(i)} = 0}
$$
Относительная угловая скорость:\\
$${\pmb \omega}_{ij}^{(i)} =\mathbf{A}^{iT} \mathbf{A}^{j} {\pmb \omega}_j^{(j)}-\pmb{\omega}_i^{(i)}$$

Относительное угловое ускорение:
$$\pmb \varepsilon_{ij}^{(i)} = -\tilde{\pmb \omega}_i^{(i)} \mathbf{A}^{iT} \mathbf{A}^{j} {\pmb \omega}_j^{(j)}+ \mathbf{A}^{iT} \mathbf{A}^{j} \tilde{\pmb \omega}_j^{(j)} {\pmb \omega}_j^{(j)} +
\mathbf{A}^{iT} \mathbf{A}^{j} {\pmb \varepsilon}_j^{(j)}-
{\pmb \varepsilon}_i^{(i)}$$
$$
\boxed{
(\mathbf{A}^{j} \mathbf{A}^{iT} \mathbf{n}_{i}^{(i)})^T {\pmb \varepsilon}_j^{(j)}-
\mathbf{n}_{i}^{(i)T} {\pmb \varepsilon}_i^{(i)}
=
\mathbf n_{i}^{(i)T}(\tilde{\pmb \omega}_i^{(i)} \mathbf{A}^{iT} \mathbf{A}^{j} {\pmb \omega}_j^{(j)})
}
$$
$$
\mathbf Q_i {\mathbf{X}}_i+ \mathbf Q_j {\mathbf{X}}_j=b_{i},
$$
$$
\mathbf Q_j = \begin{bmatrix} \mathbf 0 & \left(\mathbf{A}^j \mathbf{A}^{iT} \mathbf{n}_{ij}^{(i)}\right)^T \end{bmatrix},
\
\mathbf Q_i = \begin{bmatrix} \mathbf 0 & - \left(\mathbf{n}_{ij}^{(i)}\right)^T\end{bmatrix}, \,
b_{i}=\mathbf n_{i}^{(i)T}(\tilde{\pmb \omega}_i^{(i)} \mathbf{A}^{iT} \mathbf{A}^{j} {\pmb \omega}_j^{(j)}).
$$

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

Нарушение уравнений связи

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

Метод Й. Виттенбурга

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

Метод отдельных тел

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

Рассмотрим механическую систему, состоящую из $n$ последовательно соединенных тел: тело $1$ присоединено к телу $0$, закон движения которого известен. Запишем уравнения движения последнего тела цепи:
\begin{equation}
\mathbf{M}_n \mathbf{w}_n=\mathbf{Q}_n+\mathbf{R}_n, \ \ \ \ (1)
\end{equation}
где $\mathbf{w}_n=\{\mathbf{a}_n,\mathbf{\epsilon}_n\}^T$ - столбец линейных и угловых ускорений тела $n$; $\mathbf{M}_n$ - матрица масс; $\mathbf{Q}_n$ - обобщенная сила, действующая на тело; $\mathbf{R}_n$ - сила и момент реакции в шарнире, приведенные к центру масс. Уравнение движения предыдущего тела цепочки имеет вид
\begin{equation}
\mathbf{M}_{n-1} \mathbf{w}_{n-1}=\mathbf{Q}_{n-1}+\mathbf{R}_{n-1}-\mathbf{C}_{n}^T\mathbf{R}_n \ \ \ \ (2)
\end{equation}
здесь на тело $n-1$ кроме обобщенной силы $\mathbf{Q}_{n-1}$ и реакции $\mathbf{R}_{n-1}$ действует реакция $\mathbf{R}_{n}^*$ со стороны шарнира $n$, определяемая следующим образом:
$$
\mathbf{R}_{n}^*=
\begin{bmatrix}
\mathbf{I} & \mathbf 0 \\ \tilde{\mathbf{\rho}}_n & \mathbf I
\end{bmatrix} (-\mathbf{R}_n)=-\mathbf{C}_n^T \mathbf R_n. \ \ \ \ (3)
$$

Между ускорениями тел $n$ и $n-1$ существует связь, которая в общем виде выражается так
\begin{equation}
\mathbf{w}_n=\mathbf{C}_n \mathbf{w}_{n-1}+\mathbf{S}_n \ddot{\mathbf{q}}_n+\mathbf{w}_n', \ \ \ \ (4)
\end{equation}
где $\mathbf{q}$ - столбец обобщенных координат, задающих положение тела $n$ относительно тела $n-1$, размерность столбца $\mathbf{q}$ определяется числом степеней свободы шарнира между смежными телами. Для идеальных связей выполняется соотношение, следующее из равенства нулю работ на возможном перемещении:
\begin{equation}
\mathbf{S}_n^T \mathbf{R}_n=0. \ \ \ \ (5)
\end{equation}

Учитывая последнее соотношение, умножим (1) слева на $\mathbf{S}_n^T$ для исключения реакции:
\[
\mathbf{S}_n^T \mathbf{M}_n \mathbf{w}_n = \mathbf{S}_n^T \mathbf{Q}_n^T.
\]
В последнее выражение подставим $\mathbf{w}_n$ из (4):
\[
\mathbf{S}_n^T \mathbf{M}_n (\mathbf{C}_n \mathbf{w}_{n-1}+\mathbf{S}_n \ddot{\mathbf{q}}_n+\mathbf{w}_n') = \mathbf{S}_n^T \mathbf{Q}_n^T.
\]
Получим связь между ускорением тела и второй производной от обобщенных координат в шарнире:
\begin{equation}
\ddot{\mathbf{q}}_n = \mathbf{U}_n^{-1} \mathbf{S}_n^T \left( \mathbf{Q}_n - \mathbf{M}_n (\mathbf{C}_n \mathbf{w}_{n-1}+\mathbf{w}_n') \right), \ \ \ \ (6)
\end{equation}
где
\[
\mathbf{U}_n = \mathbf{S}_n^T \mathbf{M}_n \mathbf{S}_n. \ \ \ \ (7)
\]
Далее подставим полученное $\ddot{\mathbf{q}}_n$ в кинематические соотношения (4), подставим $\mathbf{w}_n$ в (1), выразив, таким образом, реакции $\mathbf{R}_n$. Теперь уравнение (2) можно преобразовать к виду
\[
\mathbf{M}_{n-1}^* \mathbf{w}_{n-1}=\mathbf{Q}_{n-1}^*+\mathbf{R}_{n-1},
\]
которое справедливо для любого $n=k$:
\begin{equation}
\mathbf{M}_{k-1}^* \mathbf{w}_{k-1}=\mathbf{Q}_{k-1}^*+\mathbf{R}_{k-1}, \ \ \ \ (8)
\end{equation}
где
\begin{gather}
\mathbf{M}_{k-1}^*=\mathbf{M}_{k-1}+\mathbf{C}_{k}^T \mathbf{M}_{k}^* \mathbf{C}_{k} - \mathbf{C}_{k}^T \mathbf{M}_{k}^* \mathbf{S}_{k} \mathbf{U}_{k}^{-1} \mathbf{S}_{k}^T \mathbf{M}_{k}^* \mathbf{C}_{k}, \ \ \ \ (9)\\
\mathbf{Q}_{k-1}^*=\mathbf{Q}_{k-1}+\mathbf{C}_{k}^T \left( \mathbf{M}_{k}^* \left( \mathbf{S}_{k} \mathbf{U}_{k}^{-1} \mathbf{S}_{k}^T (\mathbf{Q}_{k}^*-\mathbf{M}_{k}^* \mathbf{w}_{k}')+\mathbf{w}_{k}' \right) - \mathbf{Q}_{k}^* \right). \ \ \ \ (10)
\end{gather}

Уравнение (5) можно записать для $k=n-1,n-2,\ldots,1$, при этом в выражения для $\mathbf{M}_{k-1}^*$ и $\mathbf{Q}_{k-1}^*$ будут входить значения этих матриц на предыдущем шаге построенной рекурсии $\mathbf{M}_{k}^*$ и $\mathbf{Q}_{k}^*$. Выполнением этой последовательности операций от тела $n-1$ до тела $0$, ускорение которого известно, реализуется обратный ход алгоритма. Далее, зная ускорение тела $0$, возможно определить ускорение тела $1$ при помощи выражений (4) и (6), так реализуется прямой ход алгоритма. В выражения (4) и (6) будут входить вычисленные при обратном ходе алгоритма матрицы $\mathbf{M}_{k}^*$ и $\mathbf{Q}_{k}^*$.

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

Пример программы на языке MATLAB

Движение цепочки тел под действием силы тяжести.

Источник: Дмитроченко, О. Н. Эффективные методы численного моделирования динамики елинейных систем абсолютно твердых и деформируемых тел, Дис... канд. физ. мат. наук: 01.02.01", 2003г.