Дискретная модель троса в 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\) – радиус вектор точки закрепления.

Единичный вектор \(\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;        
    % Столбцы [r1, r2, ..., rn]
    r  = reshape(q(1:2*n),2,[]);    
    % Добавим точку закрепления троса в начало матрицы r
    r  = [[0;0] r];
    % Матрица радиус-векторов, соединяющих смежные точки
    dr = diff(r, 1, 2);
    % Расстояния между точками
    l  = sqrt(sum(dr.^2,1));    
    % Удлинение участка троса
    delta = l - p.L0;    
    % Матрица единичных векторов направления от i к i+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

© 2022. All rights reserved.

Powered by Hydejack v9.1.6