Математическая модель активного управления конвекцией в термосифоне с запаздыванием по времени Текст научной статьи по специальности «Математика»

УДК 517.958:57

Брацун Дмитрий Анатольевич

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

Красняков Иван Васильевич

студент физического факультета

ФГБОУ ВПО «Пермский государственный гуманитарно-педагогический

университет», Пермь, Россия

614990, Пермь, Сибирская, 24, (342) 238-63-64, e-mail: [email protected]; [email protected]

МАТЕМАТИЧЕСКАЯ МОДЕЛЬ АКТИВНОГО УПРАВЛЕНИЯ КОНВЕКЦИЕЙ В ТЕРМОСИФОНЕ С ЗАПАЗДЫВАНИЕМ

____л

ПО ВРЕМЕНИ*

Dmitry A. Bratsun

DS, Head of the Theoretical Physics Department

Ivan V. Krasnyakov

Student of the Faculty of Physics

Federal State Budget Educational Institution of Higher Professional Education «Perm State Humanitarian Pedagogical University» 24, Sibirskaja, 614990, Perm, Russia, e-mail: [email protected]; [email protected]

MATHEMATICAL MODEL OF ACTIVE CONTROL OF CONVECTION IN THERMOSYPHON WITH TIME DELAY

© Брацун Д.А., Красняков И.В., 2014

Работа выполнена при финансовой поддержке Министерства образования Пермского края (грант С-26/244) и Программы стратегического развития ПГГПУ (проект 031-Ф).

Аннотация. В работе теоретически рассматривается задача об автоматическом управлении движением неоднородно нагретой жидкости в термосифоне с помощью подсистемы, которая подавляет конвекцию посредством малых изменений ориентации системы в пространстве. Разработана математическая модель явления. Показано, что чрезмерное усиление обратной связи возбуждает в системе колебания, причина которых кроется в запаздывании работы контроллера. Изучены динамические режимы работы управляемого термосифона.

Ключевые слова: тепловая конвекция, активное управление, запаздывание.

Abstract. The problem of active control of the mechanical equilibrium of an inhomogeneously heated fluid in a thermosyphon is studied theoretically. The control is performed by using a feedback subsystem which inhibits convection by changing the orientation of thermosyphon in space. The mathematical model of thermosyphon dynamics is derived. It is shown that excess feedback leads to the excitation of oscillations which are related to a delay in the controller work. The dynamical regimes of a thermosyphon under external control are studied.

Keywords: thermal convection, active control, time-delay.

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

29

управления конвекцией посредством обратной связи внесла исследовательская группа во главе с Бау. В работах [15, 17] управление с обратной связью, организованное с помощью специального подогрева системы, поддерживало стационарное течение в конвективной петле тороидальной формы при значениях параметров, когда в отсутствии управления реализуется хаотический режим. В более свежих работах интересы этой исследовательской группы переместились к использованию нейросети для управления хаотической конвекцией [18], конвекцией Марангони [9], а также течениями в микрожидкостных системах [16]. Наконец, в работе [14] рассмотрена возможность управления конвекцией Релея-Бенара в плоском слое жидкости, подогреваемом снизу.

В вышеупомянутых работах в качестве управляющего воздействия использовалось изменение температуры на границе системы. Однако такая методика оправдывает себя только при небольших отклонениях параметра Рэлея от критического значения. Более того, колебания температуры на границах могут привести к возбуждению вторичной конвекции. Кроме того, к недостаткам такого метода управления стоит отнести и тот факт, что при лабораторной реализации не удается достичь равномерного изменения температуры из-за влияния конечной теплопроводности реальных жидкостей. Изменение температуры границы приводит к образованию вблизи нее теплового скин-слоя, что существенно осложняет задачу [3]. Исходя из этого, для управляющего воздействия лучше использовать модуляцию силового поля.

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

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

При экспериментальном и теоретическом изучении активного управления конвективной устойчивостью достаточно часто используется термосифон -замкнутая трубка, заполненная жидкостью, находящейся в неизотермических условиях. Такой выбор обусловлен, с одной стороны, относительно простыми для описания режимами течений, возникающими в связанных каналах, а с другой - тем, что эти режимы достаточно легко поддаются управлению. Кроме того, задача об управлении в конвективной петле интересна как модельная, описывающая в некоторых случаях управление вихревым конвективным движением в замкнутом объеме.

Важные экспериментальные исследования были приведены в работах [2, 5]. В отличие от множества работ зарубежных и российских авторов (см., например, [4, 12, 13, 17]), в качестве объекта управления был рассмотрен термосифон прямоугольной формы. Оказывается, несмотря на кажущуюся привлекательность тороидального термосифона, так популярного среди исследователей, практическая реализация механического равновесия в нем даже без управления представляет собой нетривиальную задачу. И действительно, чтобы в круглом термосифоне установилось равновесное состояние, необходимо добиться идеально вертикального градиента температуры. А это значит, что подогревать термосифон нужно весьма специальным образом, причем величина нагрева должна меняться от точки к точке конвективной петли. Видимо, именно по этой причине ни в одной из публикаций, посвященных управлению тороидальным термосифоном, экспериментально получить механическое равновесие (не говоря уж об управлении его устойчивостью) так и не удалось. Как показано в [5] получение состояния механического равновесия жидкости в прямоугольной конвективной петле - дело технически несложное. В этой работе впервые было продемонстрирована экспериментальная возможность управления механическим равновесием жидкости в конвективной системе. Были обнаружены новые явления, которые не нашли своего объяснения в теоретических работах. В частности, было выявлено, что чрезмерное усиление обратной связи возбуждает в системе колебания. Теоретическое объяснение было дано в работе [2]. Оказалось, что эти колебания вызываются запаздыванием управляющей подсистемы вносить коррекцию в состояние управляемой системы. Таким образом, как оказалось, мы сталкиваемся здесь с типичной системой, имеющей отрицательную обратную связь с запаздыванием.

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

популяционная динамика [7], нелинейные С" т В

химические реакции и вопросы генной регуляции [1, 10], механические системы [11] и т.д. Запаздывание может быть обусловлено самыми различными причинами, например, т ограниченностью скорости распространения взаимодействия (светового сигнала), наличием инерционности некоторых элементов (при управлении с обратной связью) или

В

существованием цепочек многоэтапных

последовательных реакций с известным Рис1 Схематическое изображение

прямоугольного термосифона

результатом в конце (при транскрипции генов).

В работе [2] была рассмотрена лишь линейная устойчивость состояния

механического равновесия. Кроме того, в силу ограниченности журнального

формата писем, некоторые важные моменты вывода модели не были освещены

тексте статьи. Настоящая работа посвящена подробному изложению методики

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

полученной модели.

Математическая модель. Для получения уравнений рассмотрим систему

двух связанных параллельных каналов, окружённых массивом высокой

теплопроводности. В массиве задаётся линейный градиент температуры £,

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

тяжести может динамически меняться во времени. Это изменение задаётся

единичным вектором § = -, который составляет угол р с направлением

каналов (рис. 1). Обозначим высоту и ширину термосифона как

Н и Ь соответственно, а поперечные размеры канала определим неравенствами

33

- ^ < х < ^ и - ё2 < х < ё2. Исходными уравнениями для вывода модели является система уравнений в приближении Буссинеска:

1

+ V -Уу = -1 Ур + уАу + gPT у , (1)

д г р

дТ

дТ + у .УТ = х АТ, (2)

д г

У-у = 0, (3)

где у - скорость, T - температура, p - давление. Остальные обозначения стандартны. На границах для скорости ставится стандартное условие прилипания:

V = о. (4)

Что касается температуры, то три стенки канала будем считать высокотеплопроводными, а четвёртую, лицевую к наблюдателю, -теплоизолированной (как в эксперименте [5]):

Т = -Аг, — = 0. (5)

11-3 , ду 4 ( )

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

Вектор у имеет компоненты у: (вт< ),0, сов }) и, таким образом, может вращаться в плоскости (х, г) (см. рис. 1). Угол поворота<(/) задаётся обратной связью, посредством которой и осуществляется внешнее управление конвективной системой. Как видно из уравнений (1-5), состояние равновесия в системе реализуется только в случае, когда сила тяжести и градиент температуры направлены в одну сторону < = 0. При малейшем отклонении термосифона от горизонтального положения в сторону 0 механическое равновесие жидкости без внешнего управления в системе становится невозможным.

Так как ширина каналов мала по сравнению с Н и температура жидкости в пределах поперечного сечения меняется мало, то исходные уравнения конвекции (1-5) могут быть усреднены по сечению канала. Представим поля и температуры в виде следующих аппроксимаций, удовлетворяющих граничным условиям (4,5):

[AB] и [CD]: vz = х(оcos — cos ^, к = о, к = о ;

2d1 2d2 '

y '

(6)

[BC] и [DA]: V = х(0cos—cos ^, v = о, v2 = о ; x 2d 2d2

(7)

жХ

f

[AB] и [CD]: T = &(ç, t ) cos—sin

2dx

Л

жу ж V4d2 4J

- Az ;

(8)

[BC]: T = &(ç, t ) cos — sin

2dx

(

ж у ж

\

V4d2 4j

- AH.

(9)

ЖХ

f

[DA]: T = &(ç, t ) cos-sin

2d1

Л

жу ж V4d2 4 J

(10)

Поперечное сечение конвективной петли определяется как - ^ < х < ^, - ^ < х < ^ на отрезках [АВ] и [СБ] и - ^ < х < ^, - ^ < х < ^ на отрезках [ВС]

и [ОА]. Уравнения (1-3) будем усреднять в смысле следующего скалярного произведения:

di d2 di d2 di d2 d1 ^

+ J J J...dxdyd£ +

AB-d -di BC-d -di CD-d -di DA-d, -d,

"1 "2 "1 "2 "1 "2 "1 "2 J J J ...dxdyd£+JJ J...dzdyd£ + J J J...dxdyd£ + JJ J...dzdyd£ , (11)

где £ пробегает по контуру петли (см. рис. 1).

Изложим подробнее процедуру получения динамической модели методом Галеркина. Подставляя представления (6-10) в уравнения (1-3) и интегрируя согласно (11), в результате получим следующую систему из двух уравнений.

Проекция уравнения движения (1) на базис (6-10). 1. Уравнение по части контура AB

= vAuz + gpTM cos <р (12)

ot

приводит к интегральному уравнению

Л V U1 u2 2 U1 u2

dX ff 2 ЖХ 2 Ж Ж f f 2 ЖХ 2 Ж

—H cos -cos -dxdy =-v—-XH cos -cos -

dt J J 2d 2d, 4d2 J J 2d 2d,

-^1 -Ü2 1 2 -^1 -Ü2 1 2

d1 d2 „„, „„, ^ _Л

.e^cos^f f f ©(<£,t)cos2-cos —T— sin +— d£dxdy - (13)

6 J J J y ' ои ои ^ 4^ 4 )b J v 7

£ -d -d2 2d1 2d2

d1 d2

gPAzH cos ^ J J cos cos dxdy.

d, 2d 2di

i 1 2

В результате расчета интегралов в (13) получим уравнение

—Hddt = -v -%XHd,d2 + gPcos^<f©(£t)d£-gPAzHcos^16^. (14) dt 4d 3ж j ж

2. Уравнение по части контура BC

ди

—x = vAux + gpTBC sin ty (15)

ot

приводит к интегральному уравнению

«v d1 d2 2 dl d2

ОЛ ff 2 TZ 2 Try , , T w ff 2 ж 2 Try . .

— L cos -cos -dzdy =-v—:rXL\ cos -cos -dzdy +

dt J J 2d, 2d, 4d2 J J 2d 2d,

-d1 -d2 1 2 -d1 -d2 1 2

dl d2 m f

g^sin^J J J ©(£, t)cos2-cos sin —— + — d^dzdy - (16)

!;-di -d2 2dl 2d2 V4d2 4 У

d1 d2

- gpAHL sin p J J cos-cos —— dzdy.

-d1 -d2 2d1 2d2

В результате расчёта интеграла (16) получим уравнение

—Ld1d2 = -v —^ XLd1d 2 + gtfsinp8^ <f©(£ t)d£-gPAHL sinp16^. (17) 4d 3— J —

3. Уравнение по части контура CD

до

= vA°z + gpTCD cos p (18)

dt

приводит к интегральному уравнению

d1 d2 2 d1 d2 dX rr 2 — 2 — , , — f f 2 — 2 —

—H cos -cos -dxdy =-v—- XH cos -cos -

dt J J 2d 2d, 4d2 J J 2d 2d,

di d2 щ. m, f

+ . f f Ые — _ —y

g^ cos pJ J J ©(£, t)cos2 —— cos —y- sin —^ + — d^dxdy + (19)

d1 d2

+ gPAzH cos p J J cos cos —— dxdy.

-d1 -d2 2d1 2d2 В результате расчёта интеграла (19) получим уравнение

—Hdd2 =-v Ж%XHdxd2 + gtfcosp8^{©(£,t)d£ + g¿AzHcosp16^. (20) 4d 3ж j ж

4. Уравнение по части контура DA

^ = vAux + gpTDA sin p (21)

dt

приводит к интегральному уравнению

й1 й2 2 й1 й2 дХ гг 2 Ти 2 ту , , т лгг гг 2 т 2 ту . .

— ь соб -СОБ -о— = -V—-Хь СОБ -СОБ -йгйу +

дг $ $ 2й 2й7 4й2 $ $ 2й 2й7

-йг -й2 1 2 -й! -й2 1 2

й й-,

+

ер81пр| | | ©(%,г)соб2 собб1п

(

% -йу -й2

2—х 2й2

ту т

\

(22)

V 4й2 4 J

—%—2—У.

В результате расчёта интеграла (22) получим уравнение

ОХ т2

—ь——7 = -V—-дг 4й2

хь——2 + gPsiпv ( ©(%, г)!% 3т

(23)

В результате складываем вместе четыре уравнения (ЛБ+БС+СВ+БЛ) по частям:

дХ дХ дХ дХ дХ

—И±й7 + —МЖ + —И±й7 + —МЛ? = —йЛ, 2(н + ь) (24) дг 12 дг 12 дг 12 дг 12 дг 121 '

-„-2 2 2 2 т т т т

-V—-ХИ—^ -V—-ХЬ——7 -V—-ХИ——7 -V—-ХЬ——7 = 4й2 4й2 4й2 4й2

(25)

т

= -V-

4й

2 Хй^г 2(н + ь)

еР соб р $ ©(%, г й + еР §1п р 8——^ $ ©(%, г й + еР соб р 8——^ $ ©(%, г )й%

8d1d,,

8d1d,,

3т

3т

3т

+

+

О 7 1 О Л Л ( ^

еР б1п р —1—— $©(%, г)й% = 2еР —^ соб ($©(%, г)й% + $1п р$©{%, г)%

........+ !

3т 3т

% V % %

(26)

16d1d,

- ЕРАгИ соб р-- Е^РЛНЬ б1п р-+ ЕРАгИ соб р

16d1d,,

= -еРЛИЬ Бтр

т 16d1d

т

т

1 2

т

2

(27)

Конечное уравнение движения получается следующим:

2 ( Л

— т— Х + 8еР Ч соб©(%,г)й% + 81пр(©(%,г)% - 81пр.(28)

дг 4й2 3т(И + ь)[ $ К ^ $ т (и + ь) ^ 7

V %

Проекция уравнения теплопроводности (2) на базис (6-10).

%

%

%

%

1. Уравнение по части контура АВ:

а© ?

& &-,

^ &1

дг

-&1 -&2

&1 &2 ^

II(...) - X НО+а {{(...

у-&1 -&2

2 &1 &2 „„_ ж г г 2 же . 2 = I I соб -бш

4&2 -& -&,

2&

-&1 - &2 у

V 4&2 4у

&С&у;

(29)

д©

Г 16&& & I ж , , \

& X -^-+ А—— = -г^= &&©(£,г).

дг 12 к 3ж д£ 3ж ' 12 7

&& - X

у

4&2

(30)

2. Уравнение по части контура ВС:

д© &1

&&

&1 &2

дг

&1 -&2

&1 &2

I I« - X I /(...) + АН I |(...

V-&1 - &2

2 &1 &2 ж г г 2 ж .2 = -/— I I соб —-Бт

4& 2 -& -&

2&

-&1 -&2 у

жу ^ж

V 4&2 4у

&г&у;

(31)

д©

л 1&2 XГд© + АН= г).

дг V 3ж д£ 3ж ■

&& - X

у

4&2

(32)

3. Уравнение по части контура СЭ:

д© &1

&&

Г &1 &2

дг

-& -&л

&1 &2 ^

II(...) - X ||<..) + Аг II(...

2 &1 &2

ж Г 2 2 ЖС . 2

= -% I I соб —— БШ

4& 2 -& -&

2&

-& -& у жу ж

V 4&2 4 у

&С&у;

(33)

д©

1дг

&&2 — X

Г 16&1&2 д© + Аг ^ = ^&&©(£, г).

3ж д£

3ж

4&2

(34)

4. Уравнение по части контура ЭА:

д©

дг

&1 &2

&1 &2

-& -&

&1 &2

I К.)-X / /(.)+оI /<...

V-&l -&2

2 &1 &2 ж 1 2 2 жг 2

= I соб —-БШ

-&1 -&2 у

( т, жу ж

4& 2 -& -&

2&

■ + — V 4&2 4 у

&г&у;

(35)

д©

— X

1б&1&2 а© += &&©с,г).

3ж д£ 3ж

4&2

(36)

V-& -&2

Складываем 4 уравнения (AB+BC+CD+DA) по частям:

д® . . д® . . д® . . д® . . д® , ,

-а,а9 +-а,а9 +-а,а9 +-а,а9 = 4-а,а9

дг 12 дг 12 дг 12 дг 12 дг 12

(37)

_ (1бйхй2 )_Х( 2а21

3ж2 д£

3ж

3ж2 д£

3ж

(1ва,а7 д® , 1 (1ваа7 д® 8а

- X -^ — + —— - х —^— + о 1 2

3ж2 д£

3ж

= -4 хага2

(

1в д®

ж

+ л(^

У V 8

3ж2 д£

3ж

3ж

(38)

х-

ж

аа 2 г )-х^ 2 г )-х^ 2 г)-

ж

4а2

4а2

4а2

22 жж

^ 2 ®(£, г ) = -4хж= йхй 2 ®(£, г).

4а2 4а2

(39)

Редуцированное уравнение теплопроводности получается следующим:

д® ж _/-е д „ дг 4а

1в д®

+ л(£)-

V3ж2 д^ 3жу

(40)

где Л(^) = -А для 0<£<Ни л(^) = А для Н +1<£<2Н +1.

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

дХ ж2

— = -у—^Х + —^—^ дг 4а2 3ж(Н +1)

СОБ ®(£, г № + бш ®(£, г -^щ^^ ^,(41)

2 X

дг х 4а2 ' 13ж д£ ^

3ж

(42)

Здесь равновесное состояние присутствует, его описывают последние слагаемые в уравнениях (41, 42). После усреднения система (41, 42) остаётся в частных производных, так как амплитуда температуры в (9, 10) остаётся функцией координат. Упростим эту систему с помощью приближённого решения методом Галёркина. Для этого амплитудная функция г) аппроксимируется двумя первыми модами

2

8

8

©(£, г ) = У (г )вт-

ж\С +

н+ь

+ 2 (г)

ж\с+

сов

ь

Н + ь

(43)

которые удовлетворяют граничным условиям (4, 5).

Подставляем свёртку в уравнение, где есть частные производные и решаем:

н

сов У (г )вт

н+ь

( (

-1

2УН+Ь ас, = - сов <У

( (р

ж

сов-

2

Н + Ь

= - сов <У

н + ь

ж

ж

сов-

н + ь

v_

Н + ь

ж

- сов -

ь 2

v

Л

н + ь

v

(

= - сов <У

= - сов <У

н + ь

ж

н + ь

и ь ь ь ь Л

жн + ж—+ ж— жн + ж—+ж

- 2вт

22

ж

2(н + ь)

ж жн

вт-

22

2(н + ь)

2в1п —

2 2(н + ь

= 2 сов <У

н+ь

ж

(44)

н+ь

жк+

| 2(г)с

сов< I 2(г)сов

н+ь

ь 1

2• 7н+ь

ас, = вт <2

ж\с +

ж

вт-

2

н+ь

н+ь

= вт <2

н + ь

= вт <2

= вт <2

ж

н + ь

ж\ н + ь +

ь

ж\ н +

вт-

н + ь

ж

н + ь

- вт -

ь

ь

н+ь

ь

ь

(45)

жн + жь + ж--жн-ж— жн + жь + ж—+жн + ж

ь Л

- 2в1п-

2

2

2(н + ь)

сов-

2

2(н + ь)

ж

2вт

жь

Л

2(н + ь)

сов ж

н + ь . жь

1

= 2 вт <2-вт , ч

ж 2(н + ь)

= вт <2ь.

н

0

0

н

н

2

2

Подставляем полученные значения в уравнение вместо компонент с частными производными:

дX ж2 V 16Р -Г ^РЬ 7 ЧРАЛЬ I

— = -V—-X + У собр + —^—т 7-Л. (46)

дг 4&2 Зж V Зж(Н + Ь) ж2 (Н + Ь)) к 7

Уравнение для теплопроводности распадается на две части:

— =-хж= У (г)+ /6 ч 7Х + AX, (47) дг Х 4&2 w Зж(Н + Ь) Зж2 v '

— = -х-^ 7 (г)+ /6 ч ух . (48)

дг Х 4&2 и Зж(Н + Ь) ( )

Получили систему из трёх дифференциальных уравнений:

дХ ж2 16Р .Г ^РЬ 7 ^0АНЬ I

— = -у—^Х + у собр + Б1пр —^—т 7-И, (49)

дг 4&2 Зж2 V Зж(Н + Ь) ж2 (Н + Ь)У

дУ ж2 \ 16 „„ 32 /СЛч

— = -Х^= У (г)+—7-ч 7Х + —- АХ, (50)

дг Х 4&2 ^^ Зж(Н + Ь) Зж2 v 7

д7 ж2 /ч 16 /с1\

а =~Х 4&27(г)+ЗГ(Н+Ь)УХ. (51)

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

рг) = -кУ (г -г), (52)

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

Далее, обезразмериваем систему, вводя следующие единицы измерения:

4&2

[г] = - время, (53) ж Х

[х] = Ж - длина, (54)

[X] = 9ж (н + - компонента скорости, (55)

64Ж2

[т] = (Я + 1)х - температура. (56)

4096Ж2 Ж2 ^^

Подставив безразмерные величины (53-56) в систему уравнений (4951), получим новую безразмерную систему из трёх уравнений:

ЖХ

— = Рг(- X + У ^(КУ (Г - г}) + Бт(КУ(Г - - цЯ)), (57)

Ж

— = -У + Х2 + ХЯ, (58)

Ж

— = -1 + ХУ . (59)

Ж

Появились следующие безразмерные параметры:

Рг = —- - число Прандтля, (60)

хй2

Я =--- число Рэлея, (61)

9ж8ух

К = ^ 4096 да2 сТ - параметр управления, (62)

ЗлНЬ 7тН ,

11 = ^Я—¿У ' а = ^я—I) - геометрические параметры термосифона. (63)

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

Линейный анализ устойчивости квазиравновесия. Рассмотрим сначала случай обратной связи без запаздывания. Линейный анализ устойчивости равновесия механического равновесия жидкости X* = 0, У* = 0, 2* = 0 по отношению к бесконечно малым возмущениям приводит к следующему характеристическому уравнению для инкремента 2:

(-1 - Л)(Л2 + (1 + Р)Л + Р - ЯР(1 - а К Я)) = 0 .

(64)

Как видно из (64), система (57-59) при отсутствии запаздывания допускает неустойчивость только монотонного типа, которая определяется простой формулой для нейтральной кривой, не зависящей от числа Прандтля:

К =

Я -1 аЯ 2

(65)

Рис. 2. Нейтральные кривые устойчивости механического равновесия жидкости на плоскости числа Релея Я и коэффициента усиления обратной связи К.

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

44

Из соотношения (65) хорошо видно, что в отсутствии внешнего управления к = о механическое равновесие жидкости становится неустойчивым при я = 1. Таким образом, критерий подобия я, введенный в (62), по своему смыслу явлется относительным числом Релея, которое выражается в единицах надкритичности. Предположим теперь, что запаздывание в управляющей системе имеется т ф о. Линеаризуя систему уравнений (57-59) около положения равновесия х * = о, г * = о, 2* = о, получим:

Мода, амплитуда которой описывается уравнением (68), всегда устойчива и может быть исключена из рассмотрения. Уравнения (66, 67) описывают динамические процессы, которые являются существенно разными по временным масштабам. Так как большинство экспериментов проводилось с трансформаторным маслом (р = 3 1о2), уравнение (66) описывает медленную эволюцию амплитуды х на фоне быстрых изменений г. Поэтому можно положить:

ах(г) = Р (-Х (г) + г (г) -а к яг (г - т)) ,

аг

(66)

(67)

(68)

X (г)« Г (г) -а к яг (г - т) .

(69)

Тогда вместо (67) получим одно замкнутое уравнение:

аг(г) = (Я - 1)г (г) - акя2 г (г - т) .

аг

(70)

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

для нормальных возмущений, приходим к следующему аналитическому выражению для инкремента неустойчивости х = и +(:

и = ^(А(-таКЯ2ет(1-Я)))+Я -1. (71)

т

( = 1Im(А(-т аКЯ2ет(1-Я))) , (72)

т

где А(х) - функция Ламберта, которая по определению удовлетворяет выражению А( x)exp(А(x)) = х. В частном случае Я = 1 уравнение (70) может быть решено точно, - бифуркационное значение параметра управления в этом случае равно к* = п/2та. Период колебаний равен 4т. В общем случае выражения (71) и (72) определяют целое семейство нейтральных кривых для бифуркации Хопфа (и = 0 0), параметрически зависящих от значения запаздывания т (рис. 2). Как видно из рисунка, с ростом времени запаздывания управляющей подсистемы неустойчивость спускается ниже, существенно ограничивая возможности управления механическим равновесием жидкости. Попытки управления конвективной системой при т> 10 приводят к её дестабилизации даже при тех значених параметров, когда система была устойчива без всякого управления. Таким образом, качество управления сильно зависит от способности управляющей подсистемы быстро реагировать на происходящие в термосифоне изменения.

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

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

Нелинейная динамика в области неустойчивости. Рассмотрим динамические свойства системы уравнений (57-59) с запаздывающей обратной связью (56). Как показал анализ линейной устойчивости основного состояния X* = 0, у* = 0, 2* = 0, отвечающего механическому равновесию жидкости в термосифоне, в системе присутствуют критические возмущения двух типов: монотонные и колебательные (см. рис. 2). Необходимо отметить, что с ростом числа Релея я нейтральные кривые для этих возмущений сливаются. Это достаточно необычное поведение для нейтральных кривых, так как в автономных динамических системах, как правило, колебания возникают при размерности системы не ниже двух. Соответственно нейтральная кривая для колебательной неустойчивости двукратно выражена и может перезамыкаться только с двумя монотонными кривыми. Однако здесь мы имеем дело с динамической системой с запаздыванием. В таких системах колебания могут возникать даже в одномерном случае. Таким образом, между областями монотонной и колебательной неустойчивости выше точки слиния нейтральных кривых жесткой границы нет.

Рассмотрим два характерных среза К = 0,5 и К = 3 параметрического пространства (К, Я). Вертикальная линия К = 3 пересекает нейтральную кривую в области колебательных возмущений примерно при Я « 3. На рис. 3 в проекции на плоскость (X,У) приведена эволюция динамических режимов, возникающих внутри области неустойчивости при Я > 3.

Рис.3. Эволюция динамических режимов системы (57-59) с ростом числа Релея я в области колебательной неустойчивости при т = о,3 , к = 3. Фазовые портреты приведены в проекции на плоскость (X, г)

Как видно из рис. 3, выше бифуркации Хопфа в системе возникает один

предельный цикл (см. рис. 3, я = 3,5). В точке бифуркации период колебаний

равняется т«1.66. С увеличением числа Релея амплитуда колебаний растет,

сами колебания приобретают все более сложную форму (см. рис. 3, я = 4,8).

Примерно при г«4,9 происходит вилочная бифуркация цикла, нарушающая

отражательную симметрию уравнений X , г ^-г, 2 ^ 2. В результате

48

рождается сразу два несимметричных цикла (на рис. 3 показан лишь один), а старый цикл становится неустойчивым. После этого в системе реализуется сценарий перехода к хаотическому поведению через удвоение периода. На рис. 3 приводятся циклы удвоенного (Я = 4,96) и учетверённого периода (Я = 5,0). Результат такой эволюции - рождение странного аттрактора, имеющего фрактальную структуру (см. рис. 3, Я = 5,3).

Рис.4. Эволюция фазового портрета динамической системы (57-59) с ростом числа Релея Я в области монотонной неустойчивости при т = 0,3 , К = 0,5 . Фазовые портреты приведены в проекции на плоскость (x, у). Крестиками отмечены начальные условия расчёта траекторий

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

Рассмотрим теперь, что происходит внутри области монотонной неустойчивости. Зафиксируем К = 0,5 и проследим эволюцию системы с ростом числа Релея я . Вертикальная линия К = 0,5 пересекает нейтральную кривую для монотонных возмущений примерно в точке Я «1,1. На рис. 4 приведена серия фазовых портретов системы в проекции на плоскость (X,у). Крестиками обозначены начальные условия для интегрирования соответствующей фазовой траектории. Как видно из рис. 4, ниже нейтральной кривой в системе существует единственное стационарное решение о, которое притягивает к себе все траектории (см. рис. 4, Я = 0,5). Это решение соответствует механическому равновесию жидкости X* = 0, У* = 0, 2* = 0. При пересечении нейтральной кривой в системе в силу указанной выше дискретной симметрии происходит вилочная бифуркация, в результате которой ответвляются еще два решения -о1 и о2 . Каждое из них описывает конвективную циркуляцию жидкости в термосифоне определенного направления (по или против часовой стрелке). Основное состояние о при этом становится неустойчивым (см. рис. 4, Я = 2,5). С ростом числа Релея фазовое пространство вокруг состояний равновесия о1 и о2 становится всё более закрученным. Это говорит об осциллирующем характере затухания возмущений (см. рис. 4, г = 4,0 и г = 5,0). В конце концов это приводит к вторичной бифуркации при Я«5,2, которая является бифуркацией Хопфа. От состояний равновесия о1 и о2 ответвляются устойчивые предельные циклы А1 и А2 соответственно (см. рис. 4, Я = 5,4). Период колебаний Т «1,66 совпадает с периодом колебаний предельного цикла на рис. 3, что говорит о единой природе появления этих осцилляций, связанных

с запаздыванием в системе. При дальнейшем увеличении числа Релея циклы Ä1 и Ä2 сливаются в один предельный цикл Ä (рис. 4, R = 5,8). После этого с циклом Ä повторяется сценарий, описанный для колебательной области неустойчивости (см. рис. 3).

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

Список литературы

1. Брацун Д.А., Захаров А.П. Моделирование пространственно-временной динамики циркадианных ритмов Neurospora crassa // Компьютерные исследования и моделирование. - 2011. - Т. 3, № 2. - С. 191-213.

2. Брацун Д.А., Зюзгин А.В., Половинкин К.В., Путин Г.Ф. Об активном управлении равновесием жидкости в термосифоне // ПЖТФ. - 2008. -Т. 34. - Вып. 15. - С. 36-42.

3. Гетлинг А.В. Концентрация конвективных движений у границы горизонтального слоя жидкости с неоднородным по высоте неустойчивым градиентом температуры // Изв. АН СССР, МЖГ. -1975. - № 5. - C. 45-52.

4. Дроздов С.М. Экспериментальное исследование конвекции жидкости в замкнутом тороидальном канале // Известия РАН, МЖГ. - 1995. - № 4. - C. 20-28.

5. Зюзгин А.В., Путин Г.Ф. Динамическое управление устойчивостью механического равновесия конвективной системы // Гидродинамика. -1998. - C. 123-139.

6. Келлер И.О., Тарунин Е.Л. Управление устойчивостью конвективного равновесия жидкости, подогреваемой снизу // Изв. АН СССР. МЖГ. -1990. - № 4. - C. 6-11.

7. Мюррей Дж. Математическая биология. - М.: ИКИ-РХД, 2009. - Т.1. - 774 с.

8. Чернатынский В.И., Шлиомис М.И. Конвекция вблизи критических чисел Релея при почти вертикальном градиенте температуры // Изв. АН СССР, МЖГ. - 1973. - № 1. - C. 64-70.

9. Bau H. Control of Marangoni-Benard Convection // Int. J. Heat Mass Transfer. - 1999. - Vol. 42. - P. 1327-1341.

10. Bratsun D., Volfson D., Hasty J., Tsimring L. Non-Marcovian processes in Gene Regulation // Noise in Complex Systems and Stochastic Dynamics III / edited by Laszlo B.Kish, Katja Lindenberg, Zoltan Gingl, Proceeding of SPIE. - 2005. - Vol. 5845. - P. 210-219.

11. Bratsun D.A. Effect of unsteady forces on the stability of non-isothermal particulate flow under finite-frequency vibrations // Microgravity Sci. Technol. - 2009. - Vol. 21. - P.153-158.

12. Erhard P., Muller U. Dynamical behavior of natural convection in a singlephase loop // J. Fluid Mech. - 1990. - Vol. 217. - P. 487-518.

13. Ott E., Grebogi C, Yorke J.A. Controlling chaos // Phys. Rev. Lett. - 1990. - Vol. 64. - P. 1196-1199.

14. Remillieux M., Zhao H., Bau H. Suppression of Rayleigh-Benard convection with proportional-derivative controller // Phys. Fluid. - 2007. -Vol. 19. - P. 017102-017111.

15. Singer J., Bau H. Active control of convection // Phys. Fluids A. - 1991. -Vol. 3. - P. 2859-2865.

16. Wang J., Chen Z., Qian S., Bau H. Thermally-Actuated, Phase-Change Flow Control for Microfluidic Systems // Lab on Chip. - 2005. - Vol. 5. -P. 1277-1285.

17. Wang Y., Singer J., Bau H. Controlling chaos in a thermal convection loop // J. Fluid. Mech. - 1992. - Vol. 237. - P. 479-498.

18. Yuen, P. K., Bau H. Controlling Chaotic Convection Using Neural Nets -Theory and Experiments // Neural Networks. - 1998. - Vol. 11. - P. 557569.