В этом примере показано, как загрузить данные полета и оценить силы G во время полета.
Записанные данные содержат следующие параметры полета:
угол атаки (альфа) в радианах,
угол бокового выступа (бета) в радианах,
приборная скорость (МСУ) в узлах,
угловые скорости тела (омега) в радианах/с,
нисходящие и перекрестные положения в футах, и
высота (alt) в футах.
load('astflight.mat');
Переменные MATLAB ® создаются для угла атаки (альфа), угла боковины (бета), угловых скоростей тела (омега) и высоты (альт) из записанных данных. convangvel функция используется для преобразования угловых скоростей тела из радиан в секунду (рад/с) в градусы в секунду (град/с).
alpha = fltdata(:,2); beta = fltdata(:,3); omega = convangvel( fltdata(:,5:7), 'rad/s', 'deg/s' ); alt = fltdata(:,10);
В этом наборе полетных данных была зафиксирована приборная скорость (IAS). Индикация приборной скорости (IAS) индицируется в приборе кабины. Для выполнения расчетов обычно используется истинная воздушная скорость (TAS) - скорость без ошибок измерения.
Погрешности измерений вводятся по индикаторам ПВД, по которым определяют скорость. Этими ошибками измерения являются погрешность плотности, погрешность сжимаемости и погрешность калибровки. Применение этих ошибок к истинной скорости воздуха приводит к указанной скорости.
Ошибка плотности возникает из-за меньшей плотности воздуха на высоте. Эффект заключается в том, что индикатор скорости на больших высотах имеет значение ниже истинной скорости. Когда разность или погрешность в плотности воздуха на высоте от плотности воздуха в стандартные сутки на уровне моря применяется к истинной скорости воздуха, это приводит к эквивалентной скорости воздуха (EAS). Эквивалентная скорость воздуха - это истинная скорость, измененная с изменением плотности атмосферы, которая влияет на показатель скорости.
Ошибка сжимаемости возникает из-за ограниченной способности воздуха противостоять сжатию. Эта способность снижается увеличением высоты, увеличением скорости или ограниченным объёмом. В индикаторе скорости воздуха имеется некоторое количество захваченного воздуха. При полете на больших высотах и выше воздушных скоростей калиброванная воздушная скорость (КАС) всегда выше эквивалентной. Калиброванная воздушная скорость представляет собой эквивалентную воздушную скорость, модифицированную с эффектом сжимаемости воздуха, влияющим на индикатор скорости.
Погрешность калибровки характерна для данной конструкции самолета. Погрешность калибровки является результатом положения и размещения статического (ых) вентиляционного (ых) отверстия (ей) для поддержания давления, равного атмосферному давлению внутри индикатора скорости воздуха. Положение и размещение статического вентиляционного отверстия вместе с углом атаки и скоростью самолета будут определять давление внутри указателя скорости воздушного движения и, таким образом, величину погрешности калибровки указателя скорости воздушного движения. Таблица калибровки обычно приводится в руководстве по эксплуатации летчика (POH) или в других спецификациях самолета. По данной градуировочной таблице по калиброванной скорости воздуха определяют указанную скорость (МСУ), изменяя ее с погрешностью калибровки индикатора скорости.
Ниже приведена таблица калибровки скорости для индикатора скорости самолета с нулевым отклонением закрылка. Таблица калибровки скорости воздуха преобразует указанную скорость (IAS) в калиброванную скорость (CAS) путем удаления ошибки калибровки.
flaps0IAS = 40:10:140; flaps0CAS = [43 51 59 68 77 87 98 108 118 129 140];
По приборной скорости (IAS) из таблицы калибровки полета и скорости воздуха определяют калиброванную скорость (CAS) полета.
CAS = interp1( flaps0IAS, flaps0CAS, fltdata(:,4) );
Атмосферные свойства, температуру (Т), скорость звука (а), давление (Р) и плотность (ро) определяют на высоте для стандартных суток с помощью atmoscoesa функция.
[T, a, P, rho]= atmoscoesa( alt );
Как только калиброванная воздушная скорость (CAS) и атмосферные свойства определены, истинная воздушная скорость (Vt) может быть вычислена с использованием correctairspeed функция.
Vt = correctairspeed( CAS, a, P, 'CAS', 'TAS' );
Используйте datcomimport для переноса цифровых данных DATCOM в MATLAB. Единицами измерения для этой аэродинамической информации являются футы и градусы.
data = datcomimport( 'astflight.out', true, 0 );
Его можно увидеть в выходном файле Digital DATCOM и проверить импортированные данные, которые
имеют данные только в первом альфа-значении. По умолчанию для отсутствующих точек данных установлено значение 99999. Отсутствующие точки данных заполняются значениями для первого альфа-канала, поскольку эти точки данных предназначены для использования для всех значений альфа-канала.
aerotab = {'cyb' 'cnb' 'clq' 'cmq'};
for k = 1:length(aerotab)
for m = 1:data{1}.nmach
for h = 1:data{1}.nalt
data{1}.(aerotab{k})(:,m,h) = data{1}.(aerotab{k})(1,m,h);
end
end
end
Стабильность и динамические производные в цифровой структуре DATCOM представляют собой 3-D таблицы, которые являются функциями числа Маха, угла атаки в градусах и высоты в футах. Для выполнения 3-D линейной интерполяции (interp3), индексы для производных таблиц должны быть монотонными и плоскими. Индексы этой формы генерируются meshgrid функция.
[mnum, alp, h] = meshgrid( data{1}.mach, data{1}.alpha, data{1}.alt );
Так как угловые единицы производных в градусах, единицы угла атаки (альфа) преобразуются из радианов из градусов функцией convang.
alphadeg = convang( alpha, 'rad', 'deg' );
Числа Маха для полета вычисляются функцией machnumber используя скорость звука (а) и воздушную скорость (Vt).
Mach = machnumber( convvel( [Vt zeros(size(Vt,1),2)], 'kts', 'm/s' ), a );
Закольцовывание по условиям полета, позволяет interp3 используется для линейной интерполяции таблиц производных для нахождения статических и динамических производных в этих условиях полета.
for k = length(alt):-1:1 cd(k,:) = interp3( mnum, alp, h, data{1}.cd, Mach(k), alphadeg(k), alt(k), 'linear'); cyb(k,:) = interp3( mnum, alp, h, data{1}.cyb, Mach(k), alphadeg(k), alt(k), 'linear'); cl(k,:) = interp3( mnum, alp, h, data{1}.cl, Mach(k), alphadeg(k), alt(k), 'linear'); cyp(k,:) = interp3( mnum, alp, h, data{1}.cyp, Mach(k), alphadeg(k), alt(k), 'linear'); clad(k,:) = interp3( mnum, alp, h, data{1}.clad, Mach(k), alphadeg(k), alt(k), 'linear'); end
После нахождения производных для условий полета можно вычислить аэродинамические коэффициенты.
Эталонные длины и площади, используемые при вычислении аэродинамических коэффициентов, извлекаются из цифровой структуры DATCOM.
cbar = data{1}.cbar;
Sref = data{1}.sref;
bref = data{1}.blref;
Угловые единицы производных в градусах, поэтому единицы угла бокового скольжения (бета) преобразуются из радианов из градусов функцией convang.
betadeg = convang( beta, 'rad', 'deg' );
Для вычисления аэродинамических коэффициентов необходимо задать угловые скорости тела (омега) в осях устойчивости, подобно производным. Функция dcmbody2wind формирует матрицу направления косинуса для осей тела к осям устойчивости (Tsb), когда угол бокового выступа (бета) равен нулю.
Tsb = dcmbody2wind( alpha, zeros(size(alpha)) );
Скорость изменения угла атаки (alpha_dot) также необходима для нахождения угловых скоростей в оси устойчивости (omega_stab). Функция diff используется на альфа в градусах, деленных на время выборки данных (0,50 секунды) для аппроксимации скорости изменения угла атаки (alpha_dot).
alpha_dot = diff( alphadeg/0.50 );
Последнее значение alpha_dot сохраняется для поддержания длины alpha_dot в соответствии с другими массивами в этом вычислении. Это необходимо, потому что diff функция возвращает массив, значение которого короче входного значения
alpha_dot = [alpha_dot; alpha_dot(end)];
Для полетных данных вычисляют угловые скорости по оси устойчивости (omega_stab). Угловые скорости преобразуются в матрицу 3-D, которую умножают на матрицу 3-D для матрицы направления косинуса для осей тела к осям устойчивости (Tsb).
omega_temp = reshape((omega - [zeros(size(alpha)) alpha_dot zeros(size(alpha))])',3,1,length(omega)); for k = length(omega):-1:1 omega_stab(k,:) = (Tsb(:,:,k)*omega_temp(:,:,k))'; end
Вычислите коэффициент лобового сопротивления (CD), коэффициент боковой силы (CY) и коэффициент подъема (CL). convvel функция используется для получения единиц скорости (Vt), соответствующих единицам скорости производных.
CD = cd; CY = cyb.*betadeg + cyp.*omega_stab(:,1)*bref/2./convvel(Vt,'kts','ft/s'); CL = cl + clad.*alpha_dot*cbar/2./convvel(Vt,'kts','ft/s');
Для вычисления аэродинамических сил используются аэродинамические коэффициенты сопротивления, боковой силы и подъемной силы.
Динамическое давление необходимо для расчета аэродинамических сил. Функция dpressure вычислить динамическое давление от скорости (Vt) и плотности (rho). convvel функция используется для приведения единиц скорости (Vt) в соответствие с единицами плотности (rho).
qbar = dpressure( convvel( [Vt zeros(size(Vt,1),2)], 'kts', 'm/s' ), rho );
Для нахождения сил в осях тела необходима матрица направления косинуса для осей устойчивости к осям тела (Tbs). Матрица направления косинуса для осей устойчивости к осям тела (Tbs) представляет собой транспонирование матрицы направления косинуса для осей тела к осям устойчивости (Tsb). Чтобы выполнить транспонирование массива 3-D, permute используется функция.
Tbs = permute( Tsb, [2 1 3] );
При закольцовывании через точки полетных данных вычисляются аэродинамические силы и преобразуются из осей остойчивости в оси кузова. convpres функция используется для получения единиц динамического давления (qbar), соответствующих единицам эталонной области (Sref).
for k = length(qbar):-1:1 forces_lbs(k,:) = Tbs(:,:,k)*(convpres(qbar(k),'Pa','psf')*Sref*[-CD(k); CY(k); -CL(k)]); end
В осях корпуса оценивается постоянная тяга.
thrust = ones(length(forces_lbs),1)*[200 0 0];
Оценка постоянной тяги добавляется к аэродинамическим силам, и единицы измерения преобразуются в метрические.
forces = convforce((forces_lbs + thrust),'lbf','N');
Используя рассчитанные силы, оцените силы G во время полета.
Ускорения оцениваются с использованием расчетных сил и массы, преобразованных в килограммы с помощью convmass. Ускорения преобразуются в силы G с помощью convacc.
N = convacc(( forces/convmass(84.2,'slug','kg') ),'m/s^2','G''s'); N = [N(:,1:2) -N(:,3)];
Силы G строятся над полетом.
h1 = figure; plot(fltdata(:,1), N); xlabel('Time (sec)') ylabel('G Force') title('G Forces over the Flight') legend('Nx','Ny','Nz','Location','Best')
close(h1);