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

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

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

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

Движение троса рассматривается по отношению к инерциальной системе координат \(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.\]

где \(m\) – масса узла, \(\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\) одинаковы для всех отрезков.

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

Для троса закрепленного в точке, определяемой радиус-вектором \(r_0\) \[l_1 = | \vec{r}_1 - \vec{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\).

Код

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

В файл-функцию передается время \(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

2021

2020

2019

2018


© 2020. All rights reserved.

Powered by Hydejack v9.1.4