Лазарев Ю.Ф. Mатематическое моделирование физических процессов и технических систем в MATLAB

Подождите немного. Документ загружается.

381

Рис. 9. 6. Содержимое раздела 3DoF

Первый блок позволяет моделировать продольное движение путем численного интегрирования уравне-

ний продольного движения ЛА. Второй вычисляет по заданным проекциям скорости ЛА на оси

и угол

атаки Alpha и величину вектора воздушной скорости.

На рис. 9.7. показано окно настраивания блока Equations of Motion.

Рис. 9. 7. Окно настраивания блока Equations of Motion

С его помощью вводятся следующие параметры, необходимые для численного интегрирования дифференци-

альных уравнений продольного движения.

Initial velocity (m/s)

Начальная скорость

Initial body attitude (rad)

Начальный угол подъема вектора скорости над плоско-

стью горизонта

382

Initial incidence (rad)

Начальный угол атаки

Initial body rotate rate (rad/sec)

Начальная угловая скорость тангажа

Initial position [x z] (m)

Начальное положение центра масс

Mass (kg)

Масса ЛА

Inertia (kg m^2)

Момент инерции ЛА относительно поперечной оси

Acceleration due to gravity (m/s/s)

Ускорение силы тяжести

Входы блока перечислены ниже.

Fx (N)

текущее значение силы (в ньютонах), действующей на ЛА вдоль его продольной оси

Fz (N)

текущее значение силы (в ньютонах), действующей на ЛА вдоль его нормальной оси

M (Nm)

текущее значение момента сил (в ньютонах на метр), действующего на ЛА вокруг его попе-

речной оси

Выходами блока являются нижеуказанные величины.

Altitude

(rad)

текущее значение угла между вектором скорости ЛА и плоскостью горизонта (в радиа-

нах)

q (rad/s)

текущее значение проекции угловой скорости ЛА на его поперечную ось (в радианах в

секунду)

qdot

(rad/s^2)

текущее значение проекции углового ускорения ЛА на его поперечную ось (в радианах в

секунду в квадрате)

Xe, Ze (m)

текущие координаты центра масс ЛА в продольной плоскости в земной системе коорди-

нат (в метрах)

U, w (m/s)

текущие значения проекций скорости ЛА соответственно на продольную и нормальную

оси ЛА (в метрах в секунду)

Ax, Az

(m/s^2)

текущие значения проекций ускорения ЛА соответственно на продольную и нормаль-

ную оси ЛА (в метрах в секунду в квадрате)

Рассмотренные блоки являются основными для моделирования движения ЛА. В них сосредоточены

программы, осуществляющие численное интегрирование дифференциальных уравнений. Но для их работы не-

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

на ЛА в течение его полета. Эти

силы и моменты можно разделить на три группы:

- силы и моменты, действующие на ЛА со стороны маршевого двигателя;

- силы и моменты, действующие на корпус ЛА со стороны окружающей его среды; сюда относятся си-

лы и моменты аэродинамического сопротивления движению ЛА в атмосфере и

сила тяжести ЛА;

- силы и моменты, накладываемые на ЛА со стороны управляющих органов, таких как рули, элероны,

элевоны, закрылки и т. п.

В общем случае эти силы и моменты настолько различаются для разных типов ЛА, что нельзя создать

универсальные блоки для их вычисления. Поэтому в библиотеке Aerospace имеются лишь некоторые

блоки,

осуществляющие наиболее часто встречающиеся зависимости в законах образования некоторых сил и момен-

тов.

Раздел Environment

В этом разделе библиотеки содержатся три группы блоков, позволяющих учитывать влияние парамет-

ров окружающей среды на силы и моменты, действующие на ЛА в его полете (рис. 9.8), - Atmosphere (Атмо-

сфера), Gravity (Гравитация) и Wind (Ветер).

Рис. 9. 8. Содержимое раздела Environment

383

В первой группе (Atmosphere) размещены блоки (рис. 9) ISA Atmosphere Model (ISA-модель атмосфе-

ры) и COESA Atmosphere Model (COESA-модель атмосферы).

Рис. 9. 9. Содержимое группы Atmosphere

Эти блоки рассчитывают параметры атмосферы на текущей высоте полета. Входным параметром обоих

блоков является значение текущей высоты полета ЛА. Выходные параметры приведены ниже.

Temperature (K) Температура ( в градусах Кельвина)

Speed of sound (m/s) Скорость звука (в метрах в секунду)

Air Pressure (N/m^2) Давление воздуха (ньютон на метр в квадрате)

Air Density (kg/m^3) Плотность воздуха (килограмм на метр

кубический)

Первый блок производит расчеты по международной модели стандартной атмосферы до высоты полета

20 км. Второй – по американской расширенной модели стандартной атмосферы. Здесь возможен учет особен-

ностей атмосферы до высоты 84852 м.

Параметров настраивания у этих блоков нет.

Вторая группа Gravity – содержит только один блок WGS84 Gravity Model (рис. 10), который рассчиты-

вает значение ускорения свободного падения

на текущей высоте полета ЛА и на текущей географической ши-

роте.

Рис. 9. 10. Содержимое группы Gravity

Раздел Propulsion

Раздел Propulsion (Двигатель) содержит единственный блок (рис. 11) – Turbofan Engine System (Систе-

ма турбовентиляторного двигателя).

384

Рис. 9. 11. Содержимое раздела Propulsion

Блок имеет следующие входы:

Throttle position

Текущее положение регулирующего дросселя

Mach

Текущее значение числа Маха

Altitude (m)

Текущее значение высоты полета (в метрах)

и такие выходы:

Thrust (N)

Сила тяги (в ньютонах)

Fuel flow (kg/s)

Расход горючего (в килограммах в секунду)

Окно настраивания блока приведено на рис. 12.

Рис. 9. 12. Окно настраивания блока Turbofan Engine System

В число параметров настраивания блока входят:

Initial thrust source

Источник (внутренний или внешний) начального значения тяги

двигателя

Initial thrust

Начальное значение тяги двигателя

Maximum sea level static

thrust

Максимальное значение силы тяги на уровне океана

Fastest engine time constant

at sea level static (sec)

Постоянная времени двигателя на уровне океана ( в секундах)

Sea level static thrust

specific fuel consumption

Удельный расход топлива на уровне океана

Ratio of installed thrust to

uninstalled thrust

Отношение установленной силы тяги к неустановленной

385

Разделы Actuators и GNC

Разделы Actuators (Исполнительные элементы) и GNС (Регуляторы управления движением) содержат

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

В общем случае система автоматического управления движением состоит из следующих частей:

- измерителей параметров движения ЛА – гироскопических приборов, измеряющих углы поворота

корпуса ЛА и его угловые скорости, измерителей скорости ЛА, углов атаки и дрейфа и т. п.;

- регуляторов – звеньев системы управления, которые формируют закон управления (требуемые зави-

симости регулируемых параметров от измеренных величин);

- силового привода, обеспечивающего поворот рулей, элеронов и элевонов

на требуемые углы, величи-

ны которых определены регулятором;

- исполнительных органов (рулей, элеронов и элевонов), обеспечивающих наложение на ЛА требуемых

моментов сил.

Модели измерителей параметров движения и работы исполнительных органов пользователю необхо-

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

Раздел Actuators содержит блоки, моделирующие движение силового привода рулевых

исполнитель-

ных органов (рис. 13).

Рис. 9. 13. Содержимое раздела Actuators

В нем находятся два блока – Second Order Linear Actuator (Линейный силовой привод второго порядка)

и Second Order Nonlinear Actuator (Нелинейный силовой привод второго порядка). Оба блока моделируют про-

цесс перемещения рулевого органа при подаче на вход силового привода заданного значения этого перемеще-

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

частотой собственных

колебаний и коэффициентом демпфирования.

В обоих блоках входом является текущее требуемое значение (Ac_dem) положения регулирующего ор-

гана, а выходом – действительное (Ac_ac) значение его положения.

Окна настраивания блоков изображены на рис. 14 и 15. У обоих блоков два параметра настройки оди-

наковы:

Natural frequency

Частота собственных колебаний

Damping ratio

Относительный коэффициент затухания

Initial position

Начальное положение регулирующего органа

Рис. 9. 14. Окно настраивания блока Second Order Linear Actuator

386

Рис. 9. 15. Окно настраивания блока Second Order Nonlinear Actuator

Блок нелинейного силового привода содержит еще такие параметры настраивания:

Maximum

deflection

Максимальное отклонение регулирующего органа

Minimum

deflection

Минимальное отклонение регулирующего органа

Maximum rate

Максимальная скорость отклонения регулирующего органа

Именно наличием указанных ограничений на отклонения регулирующего органа и его скорость отли-

чается нелинейный привод от линейного.

Содержимое раздела GNC показано на рис. 16. В него включены блоки, моделирующие процесс фор-

мирования сигналов, управляющих отклонением рулевых органов ЛА.

Рис. 9. 16. Содержимое раздела GNC

Кроме трех последних вспомогательных блоков, осуществляющих интерполяцию матриц, 13 блоков

этого раздела призваны вырабатывать сигнал, пропорциональный требуемому углу поворота регулирующего

органа. Поэтому выход во всех этих блоках один - u - требуемое текущее положение регулирующего органа.

В 11 блоках, представляющих различного вида регуляторы и наблюдатели, основным входом является

вектор «y» величин, характеризующих

текущее движение объекта и измеряемых на борту ЛА имеющимися

измерительными приборами. Более подробное изучение этих блоков представляет интерес, главным образом,

для специалистов в области проектирования систем автоматического управления движением ЛА и не входит в

задачу предварительного ознакомления с возможностями библиотеки Aerospace.

Раздел Transfomations

Сюда помещены блоки двух групп (рис. 17) – Axes (Оси) и Units (Единицы).

387

Рис. 9. 17. Содержимое раздела Transfomations

Группа Axes (рис. 18) включает 6 блоков, осуществляющих различные преобразования форм представ-

ления углового положения твердого тела в пространстве и 1 блок (3×3 Cross Product) векторного произведения

двух трехкомпонентных векторов.

Рис. 9. 18. Содержимое группы Axes

Представим назначение блоков преобразования координат.

Quaternions to Euler Angles

Преобразует кватернион поворота в вектор трех углов Эйлера

Quaternions to Direction

Cosine Matrix

Преобразует кватернион поворота в матрицу направляющих ко-

синусов

Euler Angles to Quaternions

Преобразует вектор трех углов Эйлера в вектор кватерниона по-

ворота

Direction Cosine Matrix to

Quaternions

Преобразует матрицу направляющих косинусов в кватернион по-

ворота

Euler Angles to Direction

Cosine Matrix

Преобразует вектор трех углов Эйлера

в матрицу направляющих косинусов

Direction Cosine Matrix to

Euler Angles

Преобразует матрицу направляющих косинусов в вектор трех уг-

лов Эйлера

Вторая группа блоков – Units – cодержит (рис. 19) 11 блоков преобразования величин из одной системы единиц

в другую.

388

Рис. 9. 19. Содержимое группы Units

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

9.2. Моделирование свободного углового движения космического

аппарата (КА)

Работу с библиотекой Aerospace начнем с создания простой модели углового движения КА в инерци-

альном пространстве.

Вариант SvDvigKA.mdl блок-схемы этой модели SvDvigKA.mdl представлен на рис. 20.

Рис. 9. 20. Модель свободного углового движения КА

В основу модели можно положить блок 6DoF (Euler Angles), обеспечивающий моделирование движе-

ние тела с шестью степенями свободы.

Если задачей моделирования является исследование свободного движения, то на вход этого блока сле-

дует подать векторные сигналы, соответствующие проекциям внешней силы и внешнего момента сил, равные

нулю. Это можно

обеспечить с помощью двух подсистем - «П/C CИЛА» (рис. 21) и «П/C МОМЕНТ» (рис. 22).

В дальнейшем будем использовать следующие обозначения.

Масса КА

Матрица моментов инерции КА

XYZ0

Вектор проекций координат центра масс КА

Вектор проекций скорости центра масс КА

389

UG0

Вектор углов поворота КА

UgSk0

Вектор проекций угловой скорости КА на оси КА

Рис. 9. 21. Подсистема «П/C CИЛА»

Рис. 9. 22. Подсистема «П/C МОМЕНТ»

На рис. 23 показаны настройки блока 6DoF (Euler Angles).

Рис. 9. 23. Настройки блока 6DoF (Euler Angles)

390

Управление работой модели и выведение результатов осуществим при помощи управляющей програм-

мы SvobDvigKA_upr, текст которой приводится ниже.

% SvobDvigKA_upr

% Управляющая программа для модели SvDvigKA

% Лазарев Ю.Ф. 29-03-2004

clear all, clc

% Установка параметров КА

J=[3400 300 -200;300 2200 100;-200 100 1400]; % Матрица моментов инерции КА

m=2000; % Масса КА

% Установка начальных условий

XYZ0=[0 0 0]; % Начальное положение КА

V0=[0 0 0]; % Начальные скорости КА

UG0=[0 0 0]; % Начальные углы КА

UgSk0=[1 0.1 0]; % Начальные угловые скорости КА

% Установка параметров интегрирования

TK=300; % Конечное время интегрирования

hi=0.1; % Шаг интегрирования

% Запуск модели

sim('SvDvigKA');

% Запись результатов интегрирования

FI=yout(:,1); TE=yout(:,2); PSI=yout(:,3);

omx=yout(:,4); omy=yout(:,5); omz=yout(:,6);

t=tout;

% Графическое представление результатов

subplot(2,2,1)

plot(TE,PSI), grid

axis('equal'); set(gca,'FontSize',12)

title('Движение оси Хb в пространстве');

xlabel('Угол \theta (градусы)'); ylabel('Угол \psi (градусы)');

subplot(2,2,3)

plot(t,TE,t,PSI,'.'), grid

set(gca,'FontSize',12); title('Углы во времени');

xlabel('Время (с)'); ylabel('Углы (градусы)'); legend(' \theta ',' \psi ',0);

subplot(2,2,4)

plot(t,omy,t,omz,'.'), grid

set(gca,'FontSize',12); title('Угловые скорости во времени');

xlabel('Время (с)'); ylabel('Угловые скорости (рад/с)'); legend(' omy ',' omz ',0);

subplot(2,2,2)

axis('off');

h=text(-0.3,1.1,'Свободное угловое движение космического аппарата','FontSize',14);

h=text(0.1,0.9,'| ','FontSize',12); h=text(0.2,0.9,num2str(J(1,1)),'FontSize',12);

h=text(0.4,0.9,num2str(J(1,2)),'FontSize',12); h=text(0.6,0.9,num2str(J(1,3)),'FontSize',12);

h=text(0.8,0.9,'| ','FontSize',12); h=text(-0.1,0.8,'J = ','FontSize',12);

h=text(0.1,0.8,'| ','FontSize',12); h=text(0.2,0.8,num2str(J(2,1)),'FontSize',12);

h=text(0.4,0.8,num2str(J(2,2)),'FontSize',12); h=text(0.6,0.8,num2str(J(2,3)),'FontSize',12);

h=text(0.8,0.8,'| ','FontSize',12); h=text(0.1,0.7,'| ','FontSize',12);

h=text(0.2,0.7,num2str(J(3,1)),'FontSize',12); h=text(0.4,0.7,num2str(J(3,2)),'FontSize',12);

h=text(0.6,0.7,num2str(J(3,3)),'FontSize',12); h=text(0.8,0.7,'| ','FontSize',12);

h=text(-0.1,0.5,'Начальные углы (градусы)','FontSize',12);

h=text(0.1,0.4,['\psi0 = ',num2str(UG0(3)*180/pi)],'FontSize',12);

h=text(0.4,0.4,['\theta0 = ',num2str(UG0(2)*180/pi)],'FontSize',12);

h=text(0.7,0.4,['\phi0 = ',num2str(UG0(1)*180/pi)],'FontSize',12);

h=text(-0.1,0.2,'Начальные угловые скорости (рад/с)','FontSize',12);

h=text(0.1,0.1,['omx0 = ',num2str(UgSk0(1))],'FontSize',12);

h=text(0.4,0.1,['omy0 = ',num2str(UgSk0(2))],'FontSize',12);

h=text(0.7,0.1,['omz0 = ',num2str(UgSk0(3))],'FontSize',12);

h=text(-0.1,-0.05,'---------------------------------------------------------------------------------------------');

h=text(-0.1,-0.1,'Программа SvobDvigKA-upr Лазарев Ю. Ф. 29-03-2004');

h=text(-0.1,-0.15,'---------------------------------------------------------------------------------------------');

Далее (рис. 24…26) приведены результаты моделирования для трех случаев различного распределения

масс КА: 1) матрица моментов инерции является диагональной (динамически симметричный КА); 2) матрица

моментов инерции является диагональной с разными моментами инерции; 3) матрица моментов инерции явля-

ется КА произвольна.