Восстановление параметров движения летательного аппарата по данным их дискретной регистрации. Ч. 1. Способы, не использующие регуляризацию Текст научной статьи по специальности «Математика»

УДК 629.762

Восстановление параметров движения летательного аппарата по данным их дискретной регистрации. Ч. 1. Способы, не использующие регуляризацию

© А.В. Плюснин МГТУ им. Н.Э. Баумана, Москва, 105005, Россия

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

Ключевые слова: газодинамический выброс, летательный аппарат, восстановление параметров движения, кусочно-полиномиальная интерполяция.

Моделирование движения летательного аппарата (ЛА) в контейнере представляет собой сложную задачу, при решении которой в общей постановке необходимо учитывать сопряженные процессы движения газовой среды, обтекающей ЛА, движения самого ЛА, теплообмена в контейнере, а также термомеханические процессы деформирования корпуса контейнера и ЛА. Решению сопряженных задач взаимодействия корпуса ЛА с окружающей газовой и жидкой средой посвящены работы [1-5].

При экспериментальной отработке систем газодинамического выброса ЛА из контейнера [6, 7] обычно измеряют основные параметры функционирования этой системы: ускорение ЛА, давление в контейнере, тепловые параметры и т. п. Во всех случаях важно определить по результатам испытаний динамические параметры, характеризующие продольное движение ЛА в контейнере. Этими параметрами являются путь, скорость и ускорение ЛА, давление в свободном пространстве контейнера за кормой ЛА, массовый расход и тяга стартовых двигателей. В конкретных ситуациях объем измеряемых параметров может быть неполным, тогда для определения (восстановления) некоторых из приведенных параметров применяют математические методы.

В работе [8] представлены примеры восстановления параметров при динамически неполных измерениях, когда при испытаниях достоверно регистрировалось либо ускорение ЛА, либо давление в контейнере, а другой из этих параметров не измеряли или же при записи произошла существенная потеря информации. В этом случае возможно удовлетворительное восстановление недостающих динамических параметров с использованием уравнения продольного движения ЛА в контейнере.

Временной масштаб физических процессов, сопровождающих движение ЛА внутри контейнера, мал по сравнению с последующими траекторными участками. Если бортовые системы регистрации параметров движения ЛА не «подстраиваются» под участок движения в контейнере, то сохраненные данные оказываются малопригодными для содержательного количественного анализа параметров газо-динамического выброса. В частности, по зафиксированным значениям ускорения ЛА невозможно восстановить интерполяцией пик ускорения. В то же время бортовые системы также фиксируют пройденный путь и скорость движения ЛА, что дает дополнительную информацию. Такой объем измерений параметров движения ЛА назовем дискретной регистрацией кинематических параметров [8]. Метод работы [8] предполагает одновременную регистрацию значений пути, скорости и ускорения ЛА и при условии, что данные записаны не реже, чем, допустим, с интервалом 0,1 с, обеспечивает удовлетворительное восстановление всей совокупности параметров движения ЛА в контейнере. В данной работе указанный подход применен для общего случая возможной неодновременной регистрации кинематических параметров.

Рассматривая вопрос методически, зададим полный набор кинематических данных: зависимости от времени V пути Ь ЛА в контейнере, скорости V и ускорения А ЛА при движении в контейнере [6, 7]. Методы реконструкции будем применять к выборкам из этого набора. Наличие полных исходных данных позволяет контролировать точность восстановления параметров. Указанные зависимости приведены на рис. 1 в безразмерном виде. Время отсчитывается от нуля и отнесено к длительности рассматриваемого в статье интервала регистрации tmsк, путь — к максимальному значению Ь^, скорость и ускорение — к отношениям Ь^/^^ и Ь^х/^х соответственно. «Полные» данные задаются с шагом 5^олн = 0,00625, в промежутках восполняются линейной интерполяцией.

Спектральный анализ и реконструкция графика ускорения ЛА. Проанализируем предварительно график ускорения А (см. рис. 1) с точки зрения его спектральных свойств. Будем считать, что А(V) = 0 вне

отрезка [0; Т ], Т = 1. Метод Фурье позволяет представить А(V) в виде и

ряда и интеграла Фурье. Выбрав интегральное представление, вычислим

2 т

, где 5(А)(ш) = |А(т)е"1 <х ёх [9], при-

0

чем вследствие вещественности А(V) имеем 5(А)(-ю) = 5(А)(ю). Поскольку спектр быстро затухает при (рис. 2), его можно ограничить некоторым значением <с = 2 п /с.

энергетический спектр

5(А)(<)

1,4

7 1,2

0

^ 1,0 <Ц

S

1 0,8

S

о

^ 0,6

§ 0,4

I

0,2

ч

|> о с

-0,2 -0,4

/ _---—

/ m \

----

f % • \ : V Lit) \

■ Л

« Ait)

4 > /

0, i о, 2 0, 3 0 -t. _0 5 0 6 0 7 0 8 0 1 9

Время I

Рис. 1. Полный набор кинематических параметров Тогда обратное преобразование Фурье для усеченного спектра

5</» H,Sffl< '

|œ|>œc

дает приближенное восстановление исходной зависимости

œc

Ac(t) = к-1 Re \ S(A)(œ)e*œt dœ« A(t).

о

На рис. 2 выделены некоторые характерные точки энергетического спектра. Выбирая приведенные значения в качестве линейной частоты fc усечения спектра, получаем ряд приближений к исходной зависимости A(t ) (рис. 3).

Интегральное преобразование Фурье (ИПФ) позволяет анализировать заданную зависимость на языке частоты œ, непрерывно изменяющейся вдоль вещественной оси. Существуют и другие представления со своими «языками». Например, в численном анализе широко применяется дискретное преобразование Фурье (ДПФ) [10, 11]. Рассмотрим вариант ДПФ, в котором упорядоченный набор чисел

[ai S взаимно однозначно преобразуется в упорядоченный набор

чисел

Ис ) С

/

*( D ) ({a, }?=-01 ))С) С

запишем это как

V

со-

гласно следующим формулам преобразований:

N-1 . N-1

,( D )

= N-1 2 aj exp(-2кijkN-1 ); aj = 2 c[D)exp(2кijkN-1 ).

j=0

k = 0

>я

л

<о »

О

3,2

2,8 2,4

2,0 1,6 1,2 0,8 0,4

0

1 1 1

1 1 г' 5"i

1 1

1 1 I

1 1

1 1 I

/V7 1 1

1 1 |

12

2 4 6 8 10

Линейная частота /

Рис. 2. Безразмерная спектральная характеристика зависимости A(t ):

1 — энергетический сектор A(t); 2-5 — / равно соответственно 1,323, 2,580,

6,864 и 13,37

14

Рис. 3. Восстановление зависимости A(t ) по усеченным спектрам:

1 — экспериментальная зависимость; 2-5 — зависимости, восстановленные в диапазонах частот соответственно 0... 1,323, 0.. .2,580, 0.. .6,864 и 0.. .13,37

Полагая aj = A (( 8т), 8т=T (N-1) 1, 7 = 0,1,..., N -1, вычислим 5(А)(ш) по формуле прямоугольников:

Т+0,5 8т т N-1

5(А)(ш) = | А(т)е-'шт ёт*-£ а! ехР -' ш7Т(-1)1

-0,58т N-1 7 = 0 1

Видно, что А} (ок)* N^(Д)(|А((8т)=0 ) лишь при дискрет-

2 кк Ш -1)

ных значениях частоты шк =----, к = 0,1.. N -1. Следовательно,

к TN

однократное применение ДПФ дает существенно меньший объем информации по сравнению с ИПФ. Для вычисения ИПФ с шагом

2к(-1) /, , „ ч

8ш«—ТN—придется применять ДПФ многократно (I = 1, 2,...),

модифицируя при этом исходный набор значений:

ТЫ (( Г _1 N-1 ^

5(А) (сок +18ш)*-) А(78т)ехр -II8ш 7Т^-1)

N-1

\

„(Л) = 0 ______________ к> fcTN

V

V 7 = 0

Однако, приняв ск ) = 0 при условии к>—- (шк >шс) и вы-

N—1

полнив обратное ДПФ, получим намного более эффективное решение задачи (см. рис. 3). Таким образом, ДПФ, малоэффективное для вычисления теоретико-функциональных аналогов-интегралов и рядов Фурье [10], является в то же время мощным аппаратом решения

задач, сформулированных на языке дискретных частот Шк = 2 ккТ-1 (к = 0,1.. N-1).

Еще одним важным результатом спектрального анализа является теорема отсчетов Уиттекера — Котельникова — Шеннона для функций с ограниченным спектром [9]. По построению функция Ас (') имеет ограниченный спектр:

+от +от

| Ас(т)е-'шт ёт = к-1 | ^е | ^(А)(л)ёе-'шт ёт =

-от -от ^ 0 J

+от Г+от 1

= (2к)-1 11 | А)(л)е'л'ёл[е-'шт ёт = А)(ш).

-от 1^-от J

Однако ее носителем является уже вся вещественная ось, а не отрезок [0;Т], как у функции А('). По теореме отсчетов функция Ас ('), спектр

которой ограничен круговой частотой qc = 2 п fc, однозначно определена своими значениями в точках tj = j 5т ( ( = 0, ±1, ± 2,...), заданных с шагом 5 т = (2 fc ) 1, и вычисляется как сумма ряда Котельникова [9]:

, ч sin qc (t -17) Ac (t) = 2 Ac (tj)-. (1)

^ W/ Qc (t - tj )

j

—от

Функция A(t) задана на конечном отрезке [0; T] и не может иметь ограниченного спектра, но при достаточно большом значении частоты усечения qc на отрезке [0; T ] справедливо равенство A(t)- Ac (t), откуда Aj = A (tj )- Ac (tj ) . Подставляя в ряд (1) значения Aj вместо значений Ac (tj ), получаем приближенное восстановление зависимости A(t):

sin Qc (t -1 i)

A(t)-Aкот(t) - 2 Aj-f^^. (2)

j e[0;T] Qc (t - tj )

В этом состоит сущность обобщенной теоремы отсчетов [9].

Для данных A на рис. 1 пик устойчиво воспроизводится при выборе fc >31...32. Результаты расчета для fc = 32 (интервал между точками отсчета составляет 64-1 - 0,0156) приведены на рис. 4.

Рис. 4. Восстановление зависимости А(^) на основании теоремы отсчетов (кружочками обозначены точки отсчета): 1 — экспериментальная зависимость; 2 — сумма ряда Котельникова

Очевидно, что при количестве точек отсчета, при котором исходная зависимость А(') удовлетворительно описывается рядом Котельникова [9], линейная интерполяция также оказывается удовлетворительной. При наличии у графика ускорения негладких участков она даже более предпочтительна. Но, поскольку линейная интерполяция при дискретной реконструкции лишь одного ускорения ЛА не подходит, то и рассмотренный в статье спектральный подход не годится. Удовлетворительное восстановление кинематических параметров ЛА при движении в контейнере в этой ситуации становится невозможным.

Дискретная реконструкция кинематических параметров ЛА при совпадающих точках регистрации. При дискретной регистрации зависимостей от времени пути, скорости и ускорения ЛА в контейнере известны значения (), Ь} = Ь (Ь))), ((), У] = V ('())),

('(А), А7 = А ('(А))) с некоторым постоянным шагом 8', который является слишком большим для успешного применения обычной интерполяции. В работе [8] рассмотрен метод восстановления по этой информации полных кинематических данных, ограниченный условием

одновременной регистрации параметров: '7 = '7) = 7) = '7А). В этом случае восстановление данных на каждом из отрезков ; '7-+1 ] можно выполнять независимо. Действительно, на концах отрезка ; '7+1 ] известны значения Ь , Ь = V и Ь=А (точками обозначены производные по времени). При кусочно-полиномиальной аппрокси-

5

мации Ь(П)(') г ] = Х а7к*к для Ь(') условия Ь(П) ('7 ) = Ь7,

У7 ; '7+1 ] к = 0

¿™ (ь) = Vj, Ь(П) ('7 ) = А в точках '7 и '7+1 позволяют составить невырожденную систему линейных алгебраических уравнений (СЛАУ) относительно коэффициентов а 7к и обеспечивают непрерывность

функций Ь(П), Ь(П) и Ь(П).

Проиллюстрируем данный подход на примере выборки из полных данных, представленных на рис. 1. Результаты вычислений для выборки, составленной относительно начальной точки '0 = 0 с шагом 8' = 0,125 , приведены на рис. 5. Восстановление данных хорошее, в отличие от сплайновой аппроксимации ускорения по заданной выборке.

Начало отсчета '0 и шаг 8' влияют на качество восстановления данных (рис. 6). Оценивая уклонение реконструированной зависимости /(П)(') от исходной /(') величиной йг = 102 /(П) - / ^ ^\/\\Ь^,

Т

выраженной в процентах, где ||/|\ =|/2 ёт , / е{Ь, V, А) , получаем

при /0 = 0 и Ы = 0,0125...0,25 ошибки «0,002...0,25, dV «0,1...2,4, йА «2,6...28, а при 5/ = 0,125 и /0 = 0...0,125 ошибки «0,01...0,03, dV «0,22...0,57, dA «6,1.13,4. В данном примере удовлетворительная реконструкция ускорения будет в диапазоне 0,016<5/<0,15 (оптимум 5/«0,04...0,05 ).

Рис. 5. Восстановление параметров при дискретной регистрации с кратными узлами (= ) = /(А), = 0, 5/=0,125):

1, 5, 8 — полные данные для А/10, Ь и V; 2—4 — для А/10 соответственно дискретные данные, сплайн по дискретным данным и реконструкция; 6, 7 и 9, 10 — дискретные данные и реконструкция соответственно для Ь и V

С ростом 5/ пик ускорения слишком занижается, а с приближением сверху к шагу дискретности полных данных 5/полн увеличивается

«зашумленность» результатов, вызванная несогласованным восполнением полных данных между известными значениями (см. рис. 6). При хорошо подготовленных данных этот эффект пропадает. Скорость и путь восстанавливаются с высокой точностью.

Маршевый (пошаговый) подход в общем случае дискретной реконструкции. В общем случае произвольно заданных узлов (с одним и тем же шагом дискретности 5/ ) будем строить кусочные полиномы

1,4

1,2

1,0

^

£ 0,8 о '

0 &

0,6

1

о

*—I

0,4

о

<0

I

¡8

о

-0,2

-0,4

8

V / / 7 < ч ■ц

4 К . 1 ( / 1 Л ч

1 \ Г ща * Л

1 ;

2 \ |/ I 111[ 5 /, \

Л А ¡К" м ¿Г 1 1 ! ' мМ? 1 « • МгЛ ¿¡X 6

0 1 с 1,2 0, ,3 0 А ' ' о, 1 ИМ } "Т 111:' II ■■ 5 0: 1 6 0, < 7 0, т г 8 0, ,9 1! ,

М

Время г

Рис. 6. Качество восстановления данных при малых и больших шагах выборки (*0 = 0,02925, Ы=0,01; *0 = 0,03925, Ы=0,1875) :

1, 6 — полные данные для АИ0 и V; 2, 4 — выборка для А/10 с шагом соответственно 0,1875 и 0,01; 3, 5, 7 и 8 — восстановление данных с шагами 0,1875 и 0,01 соответственно для А/10 и V

на отрезках

*(А). *(А) 7 ' 7+1

, чтобы через границы отрезков продолжать

графики пути и скорости, а не график ускорения, который наименее точен (как вторая производная от пути). Взаимную нумерацию узлов регистрации подчиним при этом условиям

*(А) . *(Ь) < *(А) (А) . (V) < (А)

Ч-1/ -'/ ' Ч-1^Ч -Ч ■

(3)

При некратных узлах регистрации число условий на отрезке */А); /+1 снижается с шести до четырех. Нельзя, однако, соответствующим образом понизить до третьей степени порядок кусочных полиномов, поскольку у графиков пути и скорости возникнут разрывы

при переходе через узлы */А). Рассмотрим следующее формальное

обобщение предыдущего подхода, в котором пятый порядок аппроксимирующих полиномов сохраняется. На это обобщение будем в дальнейшем ссылаться как на маршевую (пошаговую) реконструкцию.

Первый узел /0А всегда можно отнести к начальному состоянию

Ь

(/0А) (/0А) )=а(А) )=0.

(4)

Добавляя к этим условиям три условия в узлах Ь), ^), /1°А), имеем СЛАУ для вычисления коэффициентов полинома Ь(ПЩ г (А) (А)п =

у 0 ;' 1 \

5

= 2 аок/к. Решив ее, находим затем значения Ь((А) = Ь(П) (А) )г( А) (А) -,

к = 0 |/0 \

(А) = Ь(П) (А) )| (А) (А)-, что позволяет продолжить вычисления на

А) = Ь

1[/0А); 4А) \

следующем отрезке и т. д.

К сожалению, результаты реконструкции данным методом оказываются совершенно неудовлетворительными. Если расчет для случая кратных узлов регистрации обладал свойством локальности (вычисления на каждом из отрезков Г/-; //+1 \ выполнялись независимо),

\(А). .(А)"

1- ; Ч+1

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

выполня-

ются с использованием «начальных данных» Ь-А), V-А), включивших

в себя погрешности, накопленные при предыдущих вычислениях.

Покажем, что решающее значение имеет накопление ошибок аппроксимации, а не погрешностей машинной арифметики. Пусть,

например, мы имеем выборку из точных данных Ь(/) = /0,5, V (/) = 0,5 / -0,5, А(/) = -0,25 / -1,5, соответствующую расстановке узлов

5/ = 0,1, = 0,15, ^)= 0,105, /0А) = 0,085, к которой добавим начальные условия, аналогичные условию (4), но вычисленные по

заданным точным зависимостям. Полином Ь(П)(/) г (А) (А), предста-

{_'( ;\

вим таким образом, чтобы обеспечить выполнение условий на левой границе отрезка:

Ь(П)(/) ............2 3,2

[' (А);/ (А1 \

=Ь(А)-

-V]А) т + 0,5 А- т2 + т3

2 а -к тк

к = 0

где т = / -/у4"*. Легко выписать СЛАУ относительно оставшихся трех

неизвестных, методом исключения понизить ее порядок до второго и выписать точные формулы для коэффициентов по правилу Крамера. Тогда по результатам вычислений на отрезке /0А\ А) заданные

условия будут удовлетворены полиномом Ь(П)(/) г (А) (А)- точно или

с уровнем погрешности порядка машинного эпсилона [12]. В то же время начальные данные Ь^, ^^^ для вычислений на следующем

отрезке приобретают погрешности относительно точных значений — примерно 0,04 и 0,7 % соответственно, причем эти значения устойчивы относительно возмущений коэффициентов а 0 к.

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

Глобальные кусочно-полиномиальные интерполирующие функции. В пошаговом методе непрерывность реконструированных

графиков пути и скорости обеспечивалась переносом значений Ь(А), V-А), вычисленных на отрезке 'jA) , на следующий отрезок

t<jA\ в качестве левых граничных значений. Условия непре-

рывности можно автоматически выполнить, если значения Ь(А), V- А)

включить в состав узловых значений кусочно-полиномиальной интерполяции. В случае когда в узлах задается не только функция, но и ее производные, интерполяция называется эрмитовой [13]. При побудем для общ-

г (А); 1 (А) V ; 1-+1

строении такой интерполяции на отрезке

ности и для удобства записи делать различие между зарегистрированными значениями Ь-+1, V-+1, А-, А-+1 и соответствующими узловыми значениями интерполяции Ь^ , Г-+1), А-А), А((+1 . Кроме того, узловыми значениями интерполяции будут также Ь(А), V-А), Ь</A+(,

V (А) и ^) у-+1 и Ь-+1.

Каждому узлу интерполяции поставим в соответствие его «степени свободы» V, согласно заданию в этом узле значения функции

(у = 0) и(или) ее производных соответствующего порядка. Например,

XV)

если узел г- +1 является внутренним и с ним связывается лишь значение скорости Ку+С, то этот узел имеет одну степень свободы v = 1. Каждый из граничных узлов /-А) и /1 будет иметь по три степени свободы (V = 0,1,2). Для степеней свободы узлов и соответствующих

узловых значений введем краткие обозначения ^ ) и ' соответственно. Подробная расшифровка обозначений приведена в таблице:

Таблица

Индекс к 0 1 2 3 4 5 6 7 8

Узел ¿(к) 1 (А) '7 1 (А) '7 1 (А) 17 < (А) ¿7 + 1 ^ (А) ^ (А) ¿7 + 1 ^ (т) ¿7 + 1 t (V) ¿7 + 1 (V) + 1

Степень свободы V 0 1 2 0 1 2 0 1 0

Параметр д/к) Т( А) ¿7 у( А) А) Т( А) ¿7+1 А А1 т( ¿) ¿7+1 VI? ¿(V) ¿7+1

В данных обозначениях будем рассматривать интерполяции с узловыми значениями к = 0,..., п , где п = 5 (интерполяция пятой степени по крайним узлам), п = 7 (интерполяция седьмой степени по четырем

узлам, не учитывающая узловое значение ¿7+1 ) и п = 8 (полный случай). При п = 7 узел + имеет степень свободы у = 1, но не имеет степени свободы у = 0, что означает неполную эрмитову интерполяцию [14]. Если какие-либо из узлов регистрации совпадут, появятся повторы в узловых значениях, которые должны быть устранены.

Например, если = то Т^} = ¿¿1, и число узловых значений

должно уменьшиться на единицу вместе со степенью интерполяции.

Далее рассматриваются некратные узлы регистрации, но применяемые подходы без затруднений можно использовать в случаях

полностью или частично совпадающих узлов. Случай ^) = ) = ¿(А)

проще рассматривать изложенным выше методом.

Имеется ряд способов явной записи интерполяционных полиномов через узловые значения [13, 14]. Удобен подход, применяемый в методе конечных элементов (МКЭ) [15, 16], которым воспользуемся в данной работе.

Подобно локальным треугольным (тетраэдральным) координатам для плоских (пространственных) конечных элементов (КЭ) введем на

КЭ ¿(А); ¿у^! линейные координаты, причем в рассматриваемом одномерном случае достаточно только одной из них:

х=(-^л;1. (5)

Для имеющихся условий шаг регистрации 8t является постоянным для всех кинематических параметров, следовательно ( ( = 0,1... N) :

Л = г(А) -,(А) =&. = г(т) - г(¿) = г(V) - (V) "/-'7+1 +1 7 " 7 + 1 7 '

Для того чтобы оставаться в пределах данных рис. 1, номер N должен быть выбран с учетом условия < Т = 1. Определим значения

^ -X, =( -)-1; =( -*0Л))-1; Xа =1, (6)

характеризующие взаимное расположение моментов регистрации различных параметров, причем из условий (3) следуют неравенства

0<\ь <1; 0<\у <1.

Интерполяционный полином на КЭ

(Л). ,(Л)'

Ч ; Ч+1

запишем в ло-

кальных координатах:

Ь( П )(, )

= У а (]) Xк Г,(Л).,(Л) 1" 2-1 к .

к=о

(7)

Используя выражение (7), вводим обозначения для полинома

!(П)(0

Г ] ау

и его производных порядка V: п к!

/)(.)Г (Л) (Л) -ЬГУ-^аУ)Хк^=У /к(>)Хк , (8) ^ ] k=V(k-V)! к ¿0 к

(9)

^^^а^, (к = 0,1,..

где п^,= п - V .

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

,(Л). ,(Л) '0 > ^

записи уравнений МКЭ в формулировке по методу Бубнова — Га-леркина эти функции определены в глобальном смысле, так как они являются одновременно и пробными функциями. Тем не менее явные выражения этих функций обычно не требуются и для работы на КЭ их удобно нумеровать локально. В соответствии с приведенными выше обозначениями функции формы для степеней свободы узлов

КЭ Г4Л); 4Л+11 обозначим как ^) (,).

Функции формы имеют своим носителем объединение КЭ, которым принадлежит соответствующий узел, во всех остальных точках

они тождественно равны нулю. Ограничения функций формы на КЭ

(А) ¥(А) ~

Iу , +1 определяются как эрмитовы интерполяционные полиномы степени не выше п, удовлетворяющие условиям нормировки [13-16]

К(к) ()=5Н, (к,/=0..п), (10)

где 5 к/ — символ Кронекера. Эрмитовы функции формы в отличие от лагранжевых могут иметь степень ниже п. Например, кубический полином /(х) = -27/32(х + 1) (х-1) удовлетворяет пяти условиям нормировки /(-1) = /'(-1) = 0, /(1/3) = 1, /'(1/3) = 0, /(1) = 0. Подобные примеры можно получить из общего выражения для эрмитова полинома [13].

При пересечении границы КЭ первая производная от функций формы претерпевает разрыв (поэтому при решении уравнений второго порядка прибегают к слабым формулировкам). При условии своего

существования функции формы

образуют на КЭ I /уА); tI

базис п-мерного линейного пространства полиномов п-й степени. В случае полной эрмитовой интерполяции, в том числе разнопорядковой, функции формы существуют и определяются единственным образом [13]. Матрица, которую при этом приходится обращать, является обобщенной (конфлюэнтной [17]) матрицей Вандермонда. Неполная эрмитова интерполяция не обязательно осуществима, как, например, при поиске линейной функции с произвольно заданными значениями производной на концах отрезка. При п = 7 определитель

обобщенной матрицы Вандермонда равен -ХУй-7 (1 -Ль)3 х х(1-ЛУ)2 (7Л2 -4ЛУ -6ХЬХУ + )и обращается в нуль на линиях

Лу =1/7(2 + 3Л£ ±у]-9Л£ + 4 ) открытого квадрата 0<ХЬ, Лу < 1.

При возникновении проблемы можно повысить порядок интерполяционного полинома, добавляя недостающие степени свободы. Добавляя узловое значение , получаем п = 8, и определитель матрицы Ван-

дермонда 4Л1Л6 к]7 (1 -Х )3 (1 -Лу )6 (Л у -Л) не обращается в нуль

(рассмотрение данного определителя имеет смысл только при несовпадающих узлах регистрации). Следует отметить, что качество интерполяции по заданным узловым значениям может при этом ухудшиться из-за появления ложных колебаний, как показывают приведенные на

рис. 7 результаты расчетов для выборки 5t = 0,125, А) = 0,06875, Л ь = 0,65, Лу = 0,65 (откуда ^ц = 0,025, } =-0,03125).

1,4 1Д

ь 10

о >

»в

0,8 0,6 0,4

и N

8 0,2 р. 8

* о

-0,2 -0,4

-л • ■ —? 4 ' '

/ 9 ч

ч к * А \ 10 ч

\ 1 у 5 \ у < 1

( \ \ 2 / * > Л \ 7

( J I. ..(. „ 4

Г 0 I < 3 /

0 1 0 2 0 3 0 4 0 ,5 0 ,6 с ,7 С ,8 С ,9

V

Время I

Рис. 7. Эрмитова интерполяция по узловым значениям (выборка Ы = 0,125, ^А) = 0,06875, X Ь = 0,65, XV = 0,65)

1, 5, 8 — полные данные для А/10, Ь и V; 2, 6, 9 — дискретные данные для А/10, Ь и V; 3, 4 — эрмитова интерполяция для А/10 при п = 7 и п = 8; 7, 10 — эрмитова интерполяция (п = 8) для Ь и V

Условия (10) и линейная независимость функций формы позволяют записать выражение (7) в следующем виде:

Ь(п ^ )

7 ; 13+1

^ akP Xк = £ Щ") 4

(") 3

(11)

к = 0

к = 0

Коэффициенты а") находят из соотношений

а (?> V = 0,1,2;

= V = 0,1,2;

■=«Г', V = 0,1;

=аГ>, V = 0,

(12)

образующих в развернутом виде СЛАУ с обобщенной матрицей Ван-дермонда С.

Примем В = С 1 и представим решение соотношений (12) в виде

Л ;) -

=Х В«ч ч).

(13)

I = 0

В результате подстановки выражения (13) в формулу (11) получаем

хе^* (к)=хк<к)а) ).

к = 01 = 0

к = 0

Вследствие независимости параметров с]) находим отсюда функции формы в локальных координатах:

К(к V) = х Х1В1к.

I = 0

Таким образом, коэффициентами функций формы при полиномиальном разложении в локальных координатах являются столбцы матрицы В. При п = 5 матрицы С и В имеют соответственно вид:

С =

В =

1 0 0 0 0 0

0 V1 0 0 0 0

0 0 щ2 0 0 0

1 1 1 1 1 1 9

0 ]1 щ1 щ1 4] %

0 0 2 щ2 ш;2 2°л; -2

1 0 0 0 0 0

0 Л] 0 0 0 0

0 0 0,5й2 0 0 0

-10 -6к] Ч^2 10 -4Л; 0,5Л 2

15 Щ 1,5Л2 -15 -7Л; -"2

-6 -Щ -0,5л2 6 -3Л; 0,5Л2

В данном случае С является нижней блочно-треугольной матрицей [11], несложно обращаемой.

Используя выражение (13), запишем формулу (9) в виде ( V)-*,-у(к + у)!

1Г

к!

■ЕВк+у,1Ч-{1), (к-0,1,..Пу).

(14)

I — о

Теперь с учетом соотношения (3), (4) и (6) полные условия в точках регистрации примут вид ( ( — 1... N)

I(П; - 1;0)(ц)—Ц; I(П; - 1;1)(Ху) —У,; I(П; >-1;2)(1)—Л,.

(15)

В отличие от соотношений (12) здесь искомыми величинами являются узловые значения кусочно-полиномиальной интерполяции, а в правой части стоят значения зарегистрированных кинематических параметров. Кроме того, данная СЛАУ может оказаться недоопреде-ленной.

Подставляя в (15) выражения (8) и (14), получаем СЛАУ

Мх — Ъ, (16)

составленную из 3N уравнений относительно 9N узловых значений: ¿Н^ — Ъз,+у, 0<,<N-1, у — 0,1,2, (17)

I—о

Н — Е (к+у)!.уВ

Н у1 — Е-Г1-Ау Вк+У, I,

к — 0

к!

Ц+1, у —0

Ъ3,+у —' У,+1&, у—1;

Л,+1Ы 2, у — 2

Отождествим на время регистрируемые параметры с соответствующими узловыми значениями. Тогда при п — 7 и п — 8 СЛАУ

превращается в систему тождеств, а Ц(Л), у(Л) (и ЦУпри п — 8)

остаются неопределенными. При п — 5 после исключения тожде-

ственных соотношений 'j 1'2((1) = Aj ( будем иметь 2N уравнений относительно 2N искомых значений и VjA). Однако полученная при этом СЛАУ распадается на последовательно решаемые подсистемы и дает тот же результат, что и пошаговый метод. Таким образом, для получения содержательных результатов требуются дополнительные математические средства.

Выводы. В ч. 1 работы поставлена задача реконструкции параметров движения ЛА в контейнере в процессе газодинамического выброса при их регистрации с большой дискретностью, построены кусочно-полиномиальные эрмитовы аппроксимации реконструируемого графика пути, рассмотрены методы непосредственного решения результирующих СЛАУ. В случае общего расположения узлов регистрации эти методы оказались непригодными для решения поставленной задачи. Далее в ч. 2 работы на основе кусочно-полиномиальных аппроксимаций, разработанных в ч. 1, будут рассмотрены подходы, использующие вариационную формулировку задачи совместно с методом регуляризации, что позволит получить практически значимые результаты.

ЛИТЕРАТУРА

[1] Плюснин А.В. Расчет коэффициентов внутреннего и внешнего нестационарного взаимодействия корпуса ЛА с жидкостью методом граничных элементов. Математическое моделирование и численные методы, 2014, вып. 2, с. 77-100.

[2] Плюснин А.В. Учет эффекта вторичного догорания при расчетах систем газодинамического выброса летательного аппарата. Математическое моделирование и численные методы, 2014, вып. 3, с. 55-73.

[3] Димитриенко Ю.И., Захаров А.А., Коряков М.Н., Сыздыков Е.К., Минин В.В. Численное решение сопряженной задачи гиперзвуковой аэродинамики и термомеханики термодеструктирующих конструкций. Инженерный журнал: наука и инновации, 2013, вып. 9. URL: http://engjournal.ru/catalog/mathmodel/ aero/1114.html

[4] Димитриенко Ю.И., Димитриенко И.Д., Коряков М.Н. Комплексное моделирование движения летательных аппаратов и процессов горения в двигателях на основе анизотропных топлив. Инженерный журнал: наука и инновации, 2013, вып. 9. URL: http://engjournal.ru/catalog/mathmodel/technic/1120.html

[5] Димитриенко Ю.И., Захаров А.А., Коряков М.Н., Сыздыков Е.К. Моделирование сопряженных процессов аэрогазодинамики и теплообмена на поверхности теплозащиты перспективных гиперзвуковых летательных аппаратов. Известия высших учебных заведений. Машиностроение. 2014. № 3. с. 23-34.

[6] Ефремов Г.А., Страхов А.Н., Минасбеков Д.А., Горлашкин А.А., Плюснин А.В., Соколов П.М., Бондаренко Л.А., Говоров В.В. Отработка газодинамики подводного старта на наземном газодинамическом стенде предприятия. Ракетные комплексы и ракетно-космические системы — проектирование, экспериментальная отработка, летные испытания, эксплуатация: Труды секции 22 имени академика В.Н. Челомея XXXVIII Академических чтений по космонавтике. Реутов, 2014, с. 65-74.

[7] Плюснин А.В., Бондаренко Л.А., Сабиров Ю.Р. Анализ газогидродинамических процессов и методов их расчета на основе опыта предприятия в отработке подводного минометного старта. Ракетные комплексы и ракетно-космические системы — проектирование, экспериментальная отработка, летные испытания, эксплуатация: Труды секции 22 имени академика В.Н. Челомея XXXIXАкадемических чтений по космонавтике. Реутов, 2015, с. 74-83.

[8] Плюснин А.В. Реконструкция результатов газодинамических испытаний при неполном объеме измерений. Ракетные комплексы и ракетно-космические системы — проектирование, экспериментальная отработка, летные испытания, эксплуатация: Труды секции 22 имени академика В.Н. Челомея XXXIX Академических чтений по космонавтике. Реутов, 2015, с. 84-93.

[9] Харкевич А.А. Спектры и анализ. Москва, Книжный дом «ЛИБРОКОМ», 2009, 240 с.

[10] Jablon C., Simon J.-C. Applications des modèles numériques en physique. Basel: Birkhâuser Verlag, 1978, 236 p.

[11]Demmel J.W. Applied Numerical Linear Algebra. Philadelphia, SIAM, 1997, 430 p.

[12] Dennis J.E., Schnabel R.B. Numerical Methods for Unconstrained Optimization and Nonlinear Equations. New Jersey, Prentice-Hall, Inc., 1983, 440 p.

[13] Калиткин Н.Н. Численные методы. Москва, Наука, 1978, 512 с.

[14]Thompson J.F., Warsi Z.U.A., Mastin C.W. Boundary-Fitted Coordinate Systems for Numerical Solution of Partial Differential Equations — A Review. Journal of Computational Physics, 1982, vol. 47, no. 1, pp. 1-108.

[15]Norrie D.H., de Vries G. An Introduction to Finite Element Analysis. New York, Academic Press, 1978, 304 p.

[16] Ершов Н.Ф., Шахверди Г.Г. Метод конечных элементов в задачах гидродинамики и гидроупругости. Ленинград, Судостроение, 1984, 240 с.

[17] Bjorck A., Pereyra V. Solution of Vandermonde Systems of Equations. Mathematics of Computation, 1970, vol. 24, no. 112, pp. 893-903.

Статья поступила в редакцию 27.01.2016

Ссылку на эту статью просим оформлять следующим образом:

Плюснин А.В. Восстановление параметров движения летательного аппарата по данным их дискретной регистрации. Ч. 1. Способы, не использующие регуляризацию. Математическое моделирование и численные методы, 2016, № 1 (9), с. 68-88.

Плюснин Андрей Владимирович окончил Днепропетровский государственный университет по специальности «Гидроаэродинамика» и аспирантуру кафедры «Волновая и газовая динамика» МГУ им. М.В. Ломоносова. Канд. физ.-мат. наук, доцент кафедры «Вычислительная математика и математическая физика» МГТУ им. Н.Э. Баумана, заместитель начальника отдела АО «ВПК «НПО машиностроения». Автор ряда работ по нестационарным задачам гидрогазодинамики и теории упругости. e-mail: andrey.plusnin@gmail.com

Aircraft motion parameters recovery from the data of their discrete registration. Part 1. Methods without use of regularization

© A.V. Plyusnin Bauman Moscow State Technical University, Moscow, 105005, Russia

The article examines methods of aircraft motion parameters recovery from the data of their low resolution recordings in the gas-dynamic ejection experimental test. Desired conditions were satisfied by the use of Hermitian piecewise polynomial interpolation. Implementation of Tikhonov regularization provides the most flexible approach to the problem under consideration.

Keywords: gas-dynamic ejection, aircraft, motion parameters recovery, piecewise polynomial interpolation.

REFERENCES

[1] Plyusnin A.V. Matematicheskoe modelirovanie i chislennye metody — Mathematical Modeling and Computational Methods, 2014, no. 2, pp. 77-100.

[2] Plyusnin A.V. Matematicheskoe modelirovanie i chislennye metody — Mathematical Modeling and Computational Methods, 2014, no. 3, pp. 55-73.

[3] Dimitrienko Yu.I., Zaharov A.A., Koryakov M.N., Syzdykov E.K., Minin V.V. Inzhenernyy zhurnal: nauka i innovatsii — Engineering Journal: Science and Innovation, 2013, no. 9. Available at: http://engjournal.ru/catalog/mathmodel/ aero/ 1114.html.

[4] Dimitrienko Yu.I., Dimitrienko I.D., Koryakov M.N. Inzhenernyy zhurnal: nauka i innovatsii — Engineering Journal: Science and Innovation, 2013, no. 9. Available at: http://engjournal.ru/catalog/mathmodel/technic/1120.html.

[5] Dimitrienko Yu.I., Zakharov A.A., Koryakov M.N., Syzdykov E.K. Izvestiya vuzov. Ser. Mashinostroeniye — Proceedings of Higher Educational Institutions. Machine Building, 2014, no. 3, pp. 23-34.

[6] Efremov G.A., Strakhov A.N., Minasbekov D.A., Gorlashkin A.A., Plyusnin A.V., Sokolov P.M., Bondarenko L.A., Govorov V.V. Otrabotka gazodinamiki podvodnogo starta na nazemnom gazodinamicheskom stende predpriyatiya [Underwater launch gas dynamic tests on the gas dynamics ground test bench of the enterprise]. Raketnye kompleksy i raketno-kosmicheskiye sistemy — proyektirovaniye, eksperimentalnaya otrabotka, letnye ispytaniya, ekspluatatsiya: Trudy sektsii 22 XXXVIII Akademicheskikh chteniy po kosmonavtike [Proc. of the 22 section of XXXVIII academic space technology readings: Rocket and space-rocket systems — designing, experimental tests, flight tests, exploitation]. Reutov, 2014, pp. 65-74.

[7] Plyusnin A.V., Bondarenko L.A., Sabirov Yu.R. Analiz gazogidrodinamich-eskikh protsessov i metodov ikh rascheta na osnove opyta predpriyatiya v otra-botkye podvodnogo minom'etnogo starta [Analysis of gas and hydro dynamic processes and their calculation methods on the basis of the enterprise experience in underwater mortar launch tests]. Raketnye kompleksy i raketno-kosmicheskiye sistemy — proyektirovaniye, eksperimentalnaya otrabotka, letnye ispytaniya, ekspluatatsiya: Trudy sektsii 22 XXXVIII Akademicheskikh chteniy po kosmonavtike [Proc. of the 22 section of XXXVIII academic space

technology readings: Rocket and space-rocket systems — designing, experimental tests, flight tests, exploitation]. Reutov, 2015, pp. 74-83.

[8] Plyusnin A.V. Reconstruktsiya rezul'tatov gazodinamicheskikh ispytaniy pri nepolnom obyome izmereniy [Gas-dynamic tests results recovery in case of incomplete measurements]. Raketnye kompleksy i raketno-kosmicheskiye sistemy — proyektirovaniye, eksperimentalnaya otrabotka, letnye ispytaniya, ekspluatatsiya: Trudy sektsii 22 XXXVIII Akademicheskikh chteniy po kosmonavtike [Proc. of the 22 section of XXXVIII academic space technology readings: Rocket and spacerocket systems — designing, experimental tests, flight tests, exploitation]. Reu-tov, 2015, pp. 84-93.

[9] Kharkevich A.A. Spektry i analiz [Specters and analysis]. Мoscow, LIBRO-KOM Publ., 2009, 240 p.

[10] Jablon C., Simon J.-C. Applications des modèles numériques en physique. Basel, Birkhäuser Verlag, 1978, 236 p.

[11] Demmel J.W. Applied Numerical Linear Algebra. Philadelphia, SIAM, 1997, 430 p.

[12] Dennis J.E., Schnabel R.B. Numerical Methods for Unconstrained Optimization and Nonlinear Equations. New Jersey, Prentice-Hall, Inc., 1983, 440 p.

[13] Kalitkin N.N. Chislennye metody [Numerical Methods]. Мoscow, Nauka Publ., 1978, 512 p.

[14] Thompson J.F., Warsi Z.U.A., Mastin C.W. Journal of Computational Physics, 1982, vol. 47, no. 1, pp. 1-108.

[15] Norrie D.H., de Vries G. An Introduction to Finite Element Analysis. New York, Academic Press, 1978, 304 p.

[16] Yershov N.F., Shakhverdi G.G. Metod konechnykh elementov v zadachakh gidrodinamiki i gidrouprugosti [Finite Element Method in Hydro Dynamics and Hydro Elasticity Problems]. Leningrad, Sudostroyenie Publ., 1984, 240 p.

[17] Björck A., Pereyra V. Mathematics of Computation, 1970, vol. 24, no. 112, pp. 893-903.

Plyusnin A.V. graduated from the Dnipropetrovsk State University, majoring in Aero-hydrodynamics and obtained a degree in Department of Wave and Gas Dynamics, Lo-monosov Moscow State University. Cand.Sci. (Phys.-Math.), Assoc. Professor of Department of Computational Mathematics and Mathematical Physics, Bauman Moscow State Technical University, Deputy Chief of the Department in JSC MIC Mashinostroy-enia. Author of publications in the field of nonstationary problems of hydraulic gas dynamics and elasticity theory. e-mail: andrey.plusnin@gmail.com