Применение кубичного полинома Вольтерра к моделированию динамики теплообмена Текст научной статьи по специальности «Математика»

онные системы управления на транспорте. - Иркутск: ИрГУПС, 2002. - №10. - С, 176-179,

3. Мухопад Ю.Ф., Бадмаева Т.С, Синтез автоматов по декомпозированной схеме алгоритма / Информационные системы контроля и управления в промышленности и на транспорте, - Иркутск: ИрГУПС, 2005, - Вып. 13, - С, 14-25.

4. Мухопад Ю.Ф., Рудковский В.П, Структурный синтез самодиагностируемых функциональных преобразователей информации / Информационные системы контроля и управления в промышленности и на транспорте. - Иркутск: ИрГУПС, 2005. - Вып. 13, - С. 45-49.

5. Мухопад А.Ю, Структурная организация автоматов с контролем / Информационные системы контроля и управления в промышленности и на транспорте, - Иркутск: ИрГУПС, 2005, - Вып. 13. - С. 75-78,

6, Мухопад Ю.Ф.1 Мухопад А.Ю., Бадмаева Т.С. Структурная организация системы самодиагностируемых автоматов / Современные технологии, Системный анализ. Моделирование. - Иркутск: ИрГУПС.. 2005, - №1. - С, 81-85.

7, Сапожников В.В., Кравцов ЮЛ, Сапожников Вл.В. Теория дискретных устройств ж,д. автоматики и связи, - М.: Транспорт, 2001. - 304 с.

8, Щербаков Н.С, Достоверность работы цифровых устройств, - М.: Машиностроение, 1989, - 224 с,

9, Мухопад Ю.Ф. Микроэлектронные информационно-управляющие системы реального времени / Современные технологии. Системный анализ, Моделирование. - Иркутск: ИрГУПС, 2003, - Вып. 1, - С, 18-24,

С.В.Солодуша, В.А.Спиряев, М.С.Щербинин

Применение кубичного полинома Вольтерра к моделированию динамики теплообмена*

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

N ' ' _т

у{0 = - ' е (!)

т=\ о 0 '=1

((1) соответствует случаю, когда система стационарна, т.е. ядра Вольтерра Кт не зависят явно от времени; при этом из скалярности х(7) вытекает симметрия Кт по всем переменным),

Рассмотрим в качестве «эталонной» динамической системы математическую модель переходного процесса в элементе теплообменного аппарата (теплообменнике), предложенную в [1], Согласно [1], зависимость возмущения энтальпии А/(г) (кДж/кг) теплоносителя на выходе теплообменника от возмущения расхода теплоносителя (кг/с) имеет следующий вид:

Г / л

-1{\п{е)с1е -Л210(£)с1с

е " - е "

АИ-2 О

drj,te[0,T] (2)

8 (2) £>0 и <2о ~ стационарные значения расхода воды (кг/с) и теплоподвода (кВт) соответственно, £)(/) = 1)0 + , \ и Я^ - корни характеристического уравнения для некоторой системы двух обыкновенных

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

В работах [2]-[8] изучены алгоритмы построения квадратичных моделей эталонных динамических систем. Цель данной статьи - рассмотреть теоретические и прикладные аспекты построения кубичного полинома Вольтерра применительно к задаче исследования нелинейных процессов теплообмена, введенной в [7], [8].

1. Центральная проблема при аппроксимации непрерывного отображения скалярного входного сигнала х(/) в выходной у(?) с помощью (1) заключается в идентификации многомерных переходных характеристик нелинейной

динамической системы типа «вход-выход».

В последнее десятилетие методам восстановления ядер Вольтерра было посвящено большое количество работ.

* Работа поддержана грантами РФФИ 05- 01-00336 и NATO NR RIG 9812876.

Современное состояние исследований в этой области отражено, например, в [9]. Один из наиболее известных подходов к решению данной задачи основан на задании импульсных входных возмущений. Однако, как отмечено в [10], его

основной недостаток заключается в том, что многие физические процессы не допускают импульсных входов. В работах [2]-[6] развита оригинальная методика, основанная на задании (т -1) - параметрического семейства кусочно-

постоянных тестовых сигналов с амплитудами а, i = 1,ти. При этом задача идентификации Кт сводится к решению линейных интегральных уравнений Вольтерра I рода, допускающих явные формулы обращения.

Отметим однако, что для целей моделирования реакции системы на внешнее воздействие знание самих ядер Вольтерра, вообще говоря, избыточно. Достаточно уметь подсчитывать многомерные свертки, входящие в (1). Применим для их приближенного вычисления Product Integration Method (см„ например, [11]}, согласно которому при выбранном шаге сетки h в одномерном случае

Jh

if

^K(s)x(ih-s)ds « j [ K(s)ds9 х м = х /-.;' + - \h

О У=1 '~J*2(j-\)h "/+2 Л .

т

\

(3)

(этот метод особенно эффективен, если - сильно осциллирующая функция; в (3) I = 1 пк~Т), Тогда задача идентификации сводится к нахождению не самих ядер, а соответствующих интегралов от них.

Замечание 1. Данный подход реализован в работе [12] на примере эталонной динамической системы вида (1) с

К,п = — , N - 3,4, так что ш\

,п=1 т\о

Подчеркнем, что исследование [12] существенно опиралось на предположение о принадлежности симметричных

многомерных ядер Вольтерра к специальным классам, допускающим применение Р1 - аппроксимаций (аппроксимаций на основе Product Integration Method) одномерных сверток.

Естественно теперь поставить задачу идентификации полиномов Вольтерра вида (1) с использованием pi- аппроксимаций многомерных сверток

Л

, te [0,7].

(4)

ih ih

J... JK„ (J,,.., sm) J1 x(ih - sk Щ,. ,dsn 0 0 k—l

0 0 f

hh

LA

(5)

J

1 - JKis^sjds^ds,

И ''.„=1 ^ ¿ (/,-l)A (/,„-])A

\im = l,n , m - \,N) и задании многопараметрических семейств тестовых входных сигналов, введенных в [2]-[4]. Рассмотрим кубичную модель (в (1) N = 3)

t ii 2 1 t I з

Усиь (0 = \К1 (s)x(t - s)ds + J \к2 (5,, ,s2 )fl x(t - s, Щ + f J \кз (s} , s2, j3 )fj - j, )ds, . (ó)

o o

ООО

Сеточный аналог кубичного полинома Вольтерра (6), построенный на базе квадратуры (5), имеет вид

I / /

у,. = +

№

+ Л Z, 4j*JXi-j*l/2Xi-W2Xi-M/2 > / = 1> я, я/? - 7\

>! М

(7)

где

/А

т:

('- 1)А ih jh kh

ih ¡h í J

(/-l)A(./-l)A

J O] . PÍJ = J jX > ¿2

= / | С> > , /, у, к = 1, П .

(М)А(у-1)А(Аг-1)А

Замечание 2. Согласно [6], для реализации процедуры разделения откликов системы при тестовых возмущениях на составляющие, обусловленные влиянием конкретного интегрального слагаемого в (1), приходится проводить допол-

нительные тесты. В частности, при идентификации модели (6) фактически требуется -/троенное число тестовых сигналов вида

(0 = а, (е(0 - 2 • - ^ ) + - ¿ц - <»2)), 0 < щ ,со2<(<1\ (8)

с амплитудами ап / = 1,3, а{ Ф а2 Ф а3 Ф 0. Это неоправданно усложняет всю методику, если иметь в виду ее

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

Для снятия указанного недостатка будем использовать следующую тактику выбора входных возмущений, Использование откликов эталонной модели, например (2) или (4), на тестовые сигналы с шириной к

(9)

(в (8) сох — к, со2 =0), позволяет восстановить т1 и диагональные элементы д1п, В (9) е(Т) - дискрет-

ный аналог функции Хевисайда на равномерной сетке (/, = гк, г = 1,п, пк = Т), а а Ф /3 Ф 0 - вещественные числа. Прямой подсчет дает, что идентификацию р, ] и недиагональных элементов д. . к обеспечивают отклики эталонной модели на тестовый сигнал

(0 = а №) - 2) + )), ■= Iк, 1 < у + к < I < п .

(10)

п(п +1) п(п + 1)(п + 2)

Таким образом, все п Н---— ч------------- неизвестных в (7) восстановлены.

2 6

Подстановка (9) и (10) в (7) дает следующую СЛАУ

Г

У и* = а

\

Г

-/+1 ы 1 у

+ а'

к у Г J к

Л

/+1

V м /=1 У V м /=1 у

2 _ „ г г _ „ гЗ = г ,2 _ „ j2

(П)

в принятых обозначениях I. = mj, IJ = , I,Ik = pj k, /, = q,}}, /,/, = q,^ , Ijlk = qjJk, IJJ^q.,,, 1 <j + k<i<n.

С помощью системы Maple 10 [13] получено решение системы (11):

т.

Ри

2а/3(/32 -а2)

ма 4- va

2а'

(12) (13)

а а _ « __ « _ а _ а . а , а

_ Л-У-7,0 Л-1.0./-У Л,/-Д0 i,0J-j УМ,/-/-1,0 У/"—1,0,/—• ■ ПЛ)

Pi J ~~ " А ~> ~ " ' ^ ' I '

4а'

2а2^,о + + Р)У"\,й + W - «) Ai

1а2р{р2-а2)

(15)

hjjt

~ Уц-j+lJ-k Уц-jj-k* 1 + Уц-jJ-k + ^yi-lj-jj-k 2У.-Ы

J-k-l ,

6а'

. i = j, J > к , (16)

а л а , а , О,,« _ _ т;а

yjj-j+\j-k ~ ^Уi,i-j,j-k +yjj-jj-k+1 /^KZ-y-ij-A+l yt-U-jj-k ,

= —-------------■ 1 > J* J - K' i17)

6a'

4

я tr or , о; _ a _,}a , й

-y^j+ij-k ~yiti-jj-k +yj,i-jj-k+1 ^ Уi-\J-j,j-k ~ Уi-lj-jJ-k-l y,-H-j-\J -¿+1 + ■//•-i,>-y-l

12a3 i > j, j> к .

(18)

В формулах (12)-(18) ^ ^ и ^ 0, есть отклики заданной эталонной модели на сигналы (9) и (10) соответствен-

но.

Для контроля правильности алгоритмической и программной реализации предлагаемой методики проведены вычислительные эксперименты на «эталонной» математической модели (4) с известными ядрами (предполагалось

N = 4). Отметим, что новый подход в построении (6) дает нулевое отклонение от «эталона» на тестах (9), (10), с помощью которых эта интегральная модель была идентифицирована.

2, Проиллюстрируем описанную выше технику на примере математической модели теплообменного аппарата (2).

На сигналах Л£>",/7(7,) = вида (9), (10) формула (2) дает требуемый для идентификации кубичного полино-

ма набор откликов = > позволяющий найти переходные характеристики р1], к . Подста-

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

кубичной модели на входное возмущение произвольного вида в полосе [-В, В], где В = ±25%£>0; ± 50%/)0; ±75%О0.

Изложенная схема была внедрена в программный комплекс «Динамика», позволяющий проводить идентификацию и тестирование квадратичной и кубичной моделей применительно к математической модели теплообменного аппарата. Работа правильности программы и построения интегральных моделей проверялась, как и в случае с (4), на наборе сигналов, использованных для построения модели (7). Далее полагаем для всех тестов В = 75%В0, /)0 =0,16 (кг/с), Т = 30, к = 1 (сек), Такой выбор шага связан с реальными данными, полученными в ходе эксперимента на

высокотемпературном контуре (ВТК) ИСЭМ СО РАН.

Для иллюстрации эффективности использования кубичной модели в статье представлены сигналы следующего вида (рис. 1 и рис. 2):

АД/) = 0,060(0 ~ е(Г -13)) - 0,04(е(Г -13) - е(Г - 25)), t е [0,30], (19)

Д£>(/) = 0,04(е(0 - е{1 -13)) - 0,06(<?(Г -13) -- 25)), Г е [0,30]. (20)

0,06 0.04 0,02 1 1 1

•0,02 - « 0. 1 сЗ 1 04 сЗ Г 1

| -0.04 1 1 1

| -0.05 1

I -0.08 I

Рис. 1. Кусочно-постоянный сигнал (19) Рис. 2. Кусочно-постоянный сигнал (20)

Пусть у щ (/) и усиЬ (7) есть отклики квадратичной и кубичной моделей соответственно, невязки моделирования ;!лм| = (1)\ и ¡А^Сй6|| = тах^^ (?) - уе1 (7)|. Так как в приложениях важную роль играет значе-

ние отклика динамической системы в конце переходного процесса 1 = Т, введем ^ и ¡Л^Д т, которые

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

На рисунке 3 показана реакция двух интегральных моделей на заданный произвольный сигнал.

Таблица 1

Вид сигнала 1К! \Кь\\ Кь[=т

Рис. 1 16,2434 5,9679 0,9286 0,0726

Рис. 2 54,8668 8,2455 3,3307 0,3354

Рис. 4 37,7924 6,8341 2,2246 0,3795

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

Анализ результатов, полученных для различных входных сигналов х;/(Г) , у е [-В,В], показал явное преимущество кубичной модели в сравнении с квадратичной при В > 0, превышающем 75% от О0.

Рис, 3. Отклик на возмущение, изображенное на рисунке 4

Замечание 3. Дальнейшее развитие работы связано с обобщением данной техники на случай векторных входных возмущений разной физической природы и оптимальном в некотором смысле выборе а, ß в тестовых сигналах, обеспечивающих идентификацию полиномов Вольтерра в случаях N = 2,3, наиболее важных в приложениях.

Заключение. В работе на примере эталонной модели теплообмена (2) была впервые построена кубичная модель на базе метода Product Integration, развитого на случай многомерных сверток. Были аналитически получены локальные формулы для нахождения переходных характеристик модели (6). Проведено сравнение точности моделирования квадратичным и кубичным полиномами Вольтерра, построенными с помощью метода Product integration. Эффективность работы кубичной модели проиллюстрирована на типовых примерах.

Библиографический список

1. Таиров Э.А, Нелинейное моделирование динамики теплообмена в канале с однофазным теплоносителем // Изв. АН СССР. Энергетика и транспорт. - 1989. - № 1, - С. 150-156.

2. Apartsyn A.S, Mathematical modelling of the dynamic systems and objects with the help of the Volterra integral series // EPRI-SEI Joint seminar, - Beijing, China, 1991. - P, 117-132,

3. Apartsyn A.S, On some identification method for nonlinear dynamic systems II 1SEMA-92. - Shenzhen, China, 1992. - P. 288-292.

4. Апарцин А,С, О новых классах линейных многомерных уравнений I рода типа Вольтерра II Изв, вузов, Математика. - 1995. -№ 11,-С, 28-41.

5. Апарцин А.С, Теоремы существования и единственности решений уравнений Вольтерра I рода, связанных с идентификацией нелинейных динамических систем (скалярный случай). - Иркутск: СЭИ СО РАН, 1995. - № 9. - 30 с. - Препринт.

6. Апарцин А,С, Новый алгоритм моделирования нелинейных динамических систем на базе полиномов Вольтерра II Оптимизация, управление, интеллект. - 2000. - Ч. 1. - №5. - С. 26-31.

7. Апарцин A.C., Таиров Э.А., Солодуша С,В., Худяков Д.В. Применение интегро-степенных рядов Вольтерра к моделированию динамики теплообменников II Изв. РАН, Энергетика. - 2000. - № 3. - С. 138-145.

8. Apartsyn A.S. and Solodusha S.V. Mathematical Simulation of Nonlinear Dynamic Systems by Volterra Series // Engineering Simulation, - 2000. - Vol, 17. - P. 143-153.

9. Doyle III F„ Pearson R„ Ogunnaike B, Identification and Control Using Volterra Models. Springer - Verlag, - 2002,

10. Льюнг А Идентификация систем, Теория для пользователя. М: Наука, 1991, - 431 с.

11. Linz P. Product integration method for Volterra integral equations of the first kind II BIT, 1971. - Vol. 11. - P. 413-421,

12, Апарцин А.С., Спиряев В,А. Об одном подходе к идентификации полиномов Вольтерра // Оптимизация, управление, интеллект. - 2005, - № 2(10). - С. 109-117.

13. MAPLE 10 II Licensed to: Energy Systems Institute of the SB RAS, Serial number: 5GG8Z57RMT8DA5ZE.

B.B.Суржик

Моделирование динамики экранопданов

Проектирование экранопланов какой либо компоновочной схемы существенно отличается от создания самолетов такого же взлетного веса, Эти отличия связаны не столько с конструктивными особенностями, а с особыми свойствами околоэкранной эродинамики, При движении несущей поверхности (крыла) над экранирующей поверхностью (над землей, водой, льдом, снегом и т.д.) растет его подъемная сила и уменьшается сопротивление, Кроме этого, при изменении высоты крыла над экраном, существенно меняются и его динамические характеристики. В частности, кроме аэродинамического фокуса на крыле появляется экранный фокус, положение которого зависит от высоты положения крыла над экранирующей поверхностью или от отстояния крыла. Расположение этих фокусов аэродинамического и экранного в конкретной компоновочной схеме и определяет динамические свойства и устойчивость экраноплана. При изменении высоты крыла над экраном, эти виртуальные точки меняют свое место приложения, а при колебании экраноплана вокруг положения равновесия при воздействии внешнего возмущения происходит так называемая "пляска" фокусов. Эта пляска фокусов есть изменение взаимного местоположения аэродинамического и экранного фокусов в конкретной компоновочной схеме, Она и определяет устойчивость экраноплана,

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

Для крыла, движущегося над экраном, и находящегося на некотором расстоянии от центра тяжести (ЦТ) схемы, функциональная зависимость его аэродинамического коэффициента подъемной силы представим в виде:

Су = +с;да, +с;дй, +1)Д5±с;(х„ ±х,„) —^ . (1)

/и а {

где С" а, - стационарное значение коэффициента подъемной силы ¡-ого крыла;

СуАа{ - изменение коэффициента подъемной силы ¡-ого крыла при изменении угла атаки на А а,;

СкуАк, - изменение коэффициента подъемной силы 1-ого крыла при изменении его отстояния на А/г,; Су. +1)А$ - изменение С при изменении угла тангажа на величину А

Су (хт,±хРа)--- - изменение Су при вращении экраноплана вокруг поперечной оси 02с связанной

1 <№_ ¡л (И

системы координат с угловой скоростью —;

си

2т -,и =- - относительная плотность компоновочной схемы экраноплана, о а - средняя аэродинамическая

Р$2Ъах хорда крыла;

х - относительная координата носика профиля до ЦТ схемы; х Ра - относительная координата аэродинамического фокуса крыла,

Функциональная зависимость коэффициента сопротивления крыла Сх на основании аналогичных соображений имеет вид:

С, = С" а + С" Да, + С* Дй* + С* (х +1)Д(2)

/и а !