Дискретная модель троса в MATLAB

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

Трос рассматривается как система материальных точек (узлов). Узлы соединяются между собой невесомыми пружинами, которые имитируют упругие свойства троса. Предположим, что масса троса \(M\) распределяется по длине троса равномерно. В этом случае при разбиении троса на узлы масса одного узла \(m\) будет определяться выражением: \[m = \frac{M}{n},\]

где n – количество точек (узлов), на которые разбивается трос.

Движение троса рассматривается относительно неподвижной системы координат \(Ox_0 y_0\), находящейся в однородном поле силы тяжести. Ускорение свободного падения \(g\) направлено в противоположном направлении оси \(Oy_0\). Предположим, что один конец троса закреплен. Уравнение движения узла троса имеет вид: \[m \ddot{\vec{r}}_i = - \vec{F}_i + \vec{F}_{i+1} + \vec{F}_{ie}, \quad i = 1,\ldots,n-1.\]

где \(\vec{F}_i\), \(\vec{F}_{i+1}\) – силы, действующие со стороны смежных узлов, вызванные растяжением троса, \(\vec{F}_{ie}\) – внешняя сила, действующая на \(i\)-ый узел.

На последний узел действует только одна сила упругости со стороны предыдущего смежного узла, поэтому уравнение движения последнего узла имеет вид: \[m \ddot{\vec{r}}_n = - \vec{F}_n + \vec{F}_{ne}.\]

сила \(\vec{F}_k\), действующая между узлом \(k-1\) и узлом \(k\) определяется выражением \[\vec{F}_ k = c (l_k - l_0) \vec{e}_k, \quad k =1, \ldots, n\]

где \(c\) – жесткость участка троса между двумя узлами, зависящая от свойств материала троса, его толщины и свободной длины участка троса, соединяющего два узла: \[c = \frac{EF}{l_0}\]

где \(EF\) – произведение модуля упругости материала троса на площадь его поперечного сечения, \(l_0\) – свободная длина троса, соединяющего два узла. Если трос длиной \(L_0\), разбивается узлами на отрезки одинаковой длины, то значения \(l_0\) одинаковы для всех отрезков: \[l_0 = \frac{L_0}{n}\]

Расстояние между узлом \(k-1\) и узлом \(k\) – \(l_k\) определяется следующим образом: \[l_k = | \vec{r}_k - \vec{r}_{k-1} |, \quad k=2,\ldots,n\]

Длина первого участка троса, соединяющего неподвижную точку троса с первым узлом: \[l_1 = | \vec{r}_1 - \vec{r}_{0} |.\]

где \(r_0\) – радиус вектор точки закрепления.

Если длина \(l_k\) будет больше \(l_0\), то между узлами \(k-1\) и \(k\) будет действовать упругая сила \(\vec{F}_k\). Единичный вектор \(\vec{e}_k\), определяет положительное направление силы \(\vec{F}_k\), действующей со стороны узла \(k\) на узел \(k-1\). \[\vec{e}_k = \frac{\vec{r}_k - \vec{r}_{k-1}}{l_k}.\]

На узел \(k\), действует та же по модулю сила противоположного направления.

При движении троса в поле силы тяжести, как показано на рисунке, на каждый узел будет действовать сила тяжести: \[\vec{F}_{ke} = - m g \vec{e}_y\]

где \(\vec{e}_y\) – единичный вектор оси \(O y_0\).

Код

Программа моделирования движения троса будет состоять из двух файлов:

• файл-функция правых частей дифференциальных уравнений dq_tether.m;
• файл-скрипт для запуска процесса интегрирования main.m.

Файл-функция правых частей.

В файл-функцию передается время \(t\), столбец состояния системы в момент \(t\), представляющей собой последовательность пар координат узлов и скоростей узлов: \[q = [x_1, y_1, x_2, y_2, \ldots, x_n, y_n, \; \dot{x}_1, \dot{y}_1, \dot{x}_2, \dot{y}_2, \ldots, \dot{x}_n, \dot{y}_n]^T,\]

и структуру с параметрами системы p, которая будет объявлена в главном файл-скрипте.

function dq = dq_tether(t, q, p)    
    % Количество точек (узлов)
    n  = size(q,1)/4;        
    % Формируем из q матрицу столбцов координат узлов 2xn [r1, r2, ..., rn]
    r  = reshape(q(1:2*n),2,[]);    
    % Добавим точку закрепления троса в начало матрицы r
    r  = [[0;0] r];
    % Матрица радиус-векторов, соединяющих смежные точки (от узла k к узлу k+1)
    dr = diff(r, 1, 2);
    % Расстояния между узлами
    l  = sqrt(sum(dr.^2,1));    
    % Удлинение участка троса
    delta = l - p.L0;    
    % Матрица единичных векторов направления от узла k к узлу k+1
    e = dr./repmat(l,2,1);    
    % Сила между узлами. Учитывается, что только растяжение 
    % вызывает возникновение упругой силы
    F = delta*p.c.*(delta>0);    
    % Вектор силы
    Fvec = repmat(F,2,1).*e;
    % На каждый узел кроме последнего, действуют две силы    
    % Одна со знаком минус, другая со знаком плюс 
    Fp = -Fvec + [Fvec(:, 2:n) [0;0]];    
    % Ускорение    
    ap = Fp./repmat(p.point_masses,2,1);    
    % В направлении минус Y действует сила тяжести
    ag = repmat([0; -p.g],1,n);    
    % Суммарное ускорение точки (узла)
    ap = ap + ag;
    % Результат работы функции - производная вектора состояния
    dq = [q(2*n+1:end);ap(:)];    
end

Главный файл-скрипт

clc;
% Трос длиной 3 км 
L = 3000; 
% диаметром 3 мм, 
d = 0.003;
% Площадь поперечного сечения
F  = pi*d^2/4;
% плотность материала троса г/см3 -> кг/м3 (кевлар)
rho = 1.44*(1e3); 
m   = rho*pi*d^2/4*L;
fprintf('Масса троса %3.1f кг', m);
% Количество точек
n = 50;
% Массы точек (узлов) троса - массив
p.point_masses = repmat(m/n,1,n);
% На конце троса масса 10 кг
p.point_masses(n) = 10;
% Модуль упругости 83 ГПа
E   = 83e9;
% Свободная длина троса
p.L0 = L/n;
% Жесткость отрезка троса, соединяющего две точечные массы
p.c = E*F/p.L0;
% Ускорение свободного падения
p.g = 9.807;
% Формируем начальные условия
x0 = (1:n)*p.L0;  
y0 = zeros(1,n);
r0 = [x0; y0];
v0 = zeros(2,n);
% Интегрируем методом Рунге-Кутты
[t,q] = ode45(@(t,q) dq_tether(t,q,p), 0:0.5:100, [r0(:); v0(:)]);
% Анимация
for i=1:size(t,1)
    cla;
    % Положения точек
    r = [ [0;0], reshape(q(i,1:2*n),2,[])] ;
    % Точки, соединенные отрезками 
    plot(r(1,:), r(2,:),'.-');    
    % Включить сетку
    grid on;
    % Границы графика должны быть неизменны
    xlim([-L,L]); ylim([-L,L*0.1]);    
    % Масштаб по всем осям одинаков для исключения искажений    
    daspect([1 1 1]);
    % Сохранить кадр
    getframe;
end