CC BY
175. Иванова Е.А., Кривцов A.M., Морозов Н.Ф., Фирсова А.Д. Теоретическая механика. Определение эквивалентных упругих характеристик дискретных систем. СПб.: Изд-во Политехи, ун-та, 2004.
6. Подольская Е.А., Кривцов A.M. Описание геометрии кристаллов с гексогональной плотноупакованной структурой на основе парных потенциалов взаимодействия // Физ. твердого тела. 2012. 54, вып. 7. 13271334.
Поступила в редакцию 26.10.2016
УДК 532.594
ИНТЕГРАЛЬНЫЕ СООТНОШЕНИЯ ПОГРАНИЧНОГО СЛОЯ В ТЕОРИИ ВОЛНОВЫХ ТЕЧЕНИЙ КАПИЛЛЯРНЫХ ПЛЕНОК
В. Я. Шкадов1, А. Н. Белоглазкин2
Изложен обобщенный метод вывода модельных уравнений для волновых режимов течения в стекающих пленках. Применяются уравнения вязкой жидкости на основе интегральных соотношений пограничного слоя с весовыми функциями. Построена последовательность систем дифференциальных уравнений эволюционного типа с целочисленным параметром п, которым задается номер весовой функции. При п = 0 получается классическая модель IBL, при п ^ 1 — ее модификации, названные моделями WIBL. Проведены численные расчеты решений в линейном и нелинейном приближениях при п = 0,1, 2. Сопоставление с экспериментальными данными и численными решениями исходных дифференциальных уравнений гидродинамики показало, что, как правило, для течений пленок на вертикальной и наклонной твердых поверхностях наиболее точные результаты получаются при использовании модели п = 0, а с увеличением п точность решений падает, что говорит о преимуществе классической модели IBL перед моделями WIBL.
Ключевые слова: нелинейные волны, неустойчивость, пленка, капиллярность.
A generalized method of deriving the model equations is considered for wave flow regimes in falling liquid films. The viscous liquid equations are used on the basis of integral boundary-layer relations with weight functions. A family of systems of evolution differential equations is proposed. The integer parameter n of these systems specifies the number of a weight function. The case n = 0 corresponds to the classical IBL model. The case n > 1 corresponds to its modifications called the WIBL models. The numerical results obtained in the linear and nonlinear approximations for n = 0,1, 2 are discussed. The numerical solutions to the original hydrodynamic differential equations are compared with experimental data. This comparison leads us to the following conclusions: as a rule, the most exact solutions are obtained for n = 0 in the case of film flows on vertical and inclined solid surfaces and the accuracy of solutions decreases with increasing n. Hence, the classical IBL model has an advantage over the WIBL models.
Key words: nonlinear waves, instability, film, capillarity.
1. Введение. Началом экспериментальных исследований гравитационно-капиллярных волн в жидких пленках послужила статья П.Л. Капицы и С.П. Капицы [1], в которой обнаружена возможность существования не только периодического режима, но и режима, состоящего из одиночных волн, представляющих "новое явление, установленное опытом". Первый опыт теоретического исследования волновых движений в пленке изложен в работе [2], где показано, что применение приближения пограничного слоя может быть успешным в случае периодического режима. Что касается случая одиночных волн, то там же сделано заключение: "Количественное описание свойств одиночных волн в тонком слое текущей жидкости при наличии вязкости и поверхностного натяжения будет
1 Шкадов Виктор Яковлевич — доктор физ.-мат. наук, проф. каф. аэромеханики и газовой динамики мех.-мат. ф-та МГУ, e-mail: shkadovQmech.math.msu.su.
2 Белоглазкин Александр Николаевич — канд. физ.-мат. наук, доцент каф. аэромеханики и газовой динамики мех.-мат. ф-та МГУ, e-mail: belQmech.math.msu.su.
сопряжено со столь серьезными математическими трудностями, что преодоление их едва ли можно предвидеть. Поэтому изучение таких одиночных волн придется вести экспериментально". Решить эту проблему оказалось возможным методами математики и гидромеханики. В работе [3] в развитие идеи о целесообразности использования приближения пограничного слоя построена математическая модель нестационарного волнового течения в пленке, решения которой позволили истолковать все имевшиеся экспериментальные результаты и наметить пути к новым содержательным экспериментам.
В течение многих лет математическая модель статьи [3] служила основным средством для исследования новых, часто неожиданных свойств волновых движений в тонких слоях вязких жидкостей. Важным полезным свойством этой модели является ее универсальность, позволяющая включать в рассмотрение новые разнообразные механические и физико-химические факторы. К настоящему времени на основе этого подхода рассмотрены волновые течения в пленке ньютоновских, неньютоновских и сильновязких жидкостей, течения с растворенным поверхностно-активным веществом, волны в пленке на тонком цилиндре, вращающемся диске и на поверхности с периодическим и ступенчатым рельефом, а также неустойчивость течений и нелинейные волны в однослойных и двухслойных пленках. Изучены воздействия на пленочное течение граничного воздушного потока и внешнего электрического поля, влияние вибраций и нагрева удерживающей поверхности. В рамках подхода статьи [3] рассмотрены сопряженные задачи волновой динамики и тепло- и массопереноса в пленках и капиллярных струях.
Полученные к настоящему времени результаты по гидродинамике и физико-химическим свойствам в волновых пленках жидкости, а также исследуемые новые задачи представлены в широком и активно пополняемом потоке научных публикаций. Построенная впервые в [3] модельная система эволюционных дифференциальных уравнений, получившая название "модель Шкадова", с ее многочисленными приложениями приведена в монографиях [4-6]. Подробный анализ и развитие метода сведения краевых задач для жидких пленок к модельным эволюционным дифференциальным уравнениям, получившего название "метод Капицы-Шкадова", представлены в монографиях [7, 8]. Этапы развития новых идей и достигнутые результаты, касающиеся различных внешних воздействий и процессов тепло- и массопереноса на волновые течения, отражены в обзорах [9-13].
При разумных значениях управляющих внешних параметров, при которых корректность вывода приближенной модели статьи [3] гарантирована, расчеты нелинейных волн по этой модели мало отличаются от имеющихся единичных расчетов, выполненных на основании полной системы Навье-Стокса [14], и от соответствующих экспериментальных данных [5]. Преимущество этой модели связано, в частности, с тем, что с ее помощью можно выполнить расчеты большого количества вариантов и учесть в них влияние всех внешних существенных параметров задачи. Это позволяет не только истолковать отдельные варианты волновых потоков, но и вывести из их анализа заключения о свойствах решений в целом.
В настоящей работе вывод модельных дифференциальных уравнений для волновых пленок основывается на применении обобщенного интегрального соотношения пограничного слоя, предложенного В.В. Голубевым [15]. В этом соотношении выполняется осреднение уравнения движения по поперечной координате с использованием весовых функций. Выводится последовательность модельных уравнений, содержащая произвольный целочисленный параметр п. Ставится цель исследовать точность аппроксимации исходной краевой задачи модельными уравнениями при п = 0,1, 2, 3,....
2. Приближение пограничного слоя и метод интегральных соотношений. Течения тонкого слоя несжимаемой жидкости по наклонной стенке (угол наклона в) под действием сил тяжести, вязкости и поверхностного натяжения описываются системой уравнений пограничного слоя с самоиндуцированным давлением
55 (щ + иих + ьиу) = Iгххх - екх + 1 + ^ иуу,
6 (1)
их + уу = 0
и граничными условиями на твердой удерживающей поверхности и на свободной поверхности жидкости
и = 0, V = 0, у = 0;
иу = 0, кг + икх = у, у = к{х,1). Здесь и, и — компоненты скорости, к — толщина слоя жидкости.
Вывод краевой задачи (1), (2) в приближении пограничного слоя из исходной полной системы уравнений Навье-Стокса и соответствующих граничных условий и его обоснование представлены в работе [3]. Выражения для безразмерных параметров 5, е и оценки допустимых значений этих величин даются, в частности, в [16, 17] и ряде других статей. Отметим здесь, что при приведении исходной задачи к безразмерной форме (1), (2) принимается предположение, что на рассматриваемых классах течений силы тяжести, вязкости и поверхностного натяжения имеют одинаковый порядок и что отношение характерных масштабов длин 1У/1Х = к есть малая величина. Эти соображения выражаются следующими соотношениями:
к2 3 1 1
We кЫе кЕг 55
(3)
Здесь \¥е = , Ег = -^Ц Ые = — числа Вебера, Фруца, Рейнольдса; /го, Щ — характерные
а дпо V
. 15с^
значения толщины слоя и осредненнои скорости жидкости, д = д$е =-о, р — плотность
Ке
жидкости, а — коэффициент поверхностного натяжения, и — коэффициент кинематической вязкости.
Соотношениями (3) вводятся в рассмотрение два параметра 5 и к:
А КП/9 к2/9 Р Г 4-4-1/3
о =-"тт, к = —7тт> = ЗКе, 7 = — [и д) — число Капицы.
4571/з 71/з р
Если 7^1, что имеет место для жидкостей малой вязкости типа воды, то при конечных, но умеренных расходах жидкости в слое получаем к2 <С 1. В этом случае при выводе системы (1) из полной задачи Навье-Стокса можно обоснованно отбросить члены порядка о(к2) и получить систему пограничного слоя (1) с одним безразмерным параметром 5 для слоя на вертикальной стенке и с дополнительным параметром е для слоя на наклонной стенке. Для сильновязких жидкостей условие 7^1 может заметно ослабляться и возникает необходимость учитывать члены порядка о(к2) при выводе (1). Соответствующая система, обобщающая (1), выведена в [16], там же исследованы свойства линейных и нелинейных волновых решений в зависимости от двух основных параметров 5 и к для слоя на вертикальной стенке (е = 0).
Применим к исследованию решений системы (1) метод интегральных соотношений, который по существу соответствует первому шагу метода Галеркина-Петрова. Краевую задачу для поля скоростей и(у, х, /) запишем в виде
Ци) = 0,
(4)
и = 0, у = 0; иу = 0, у = к{х, /). Зависимость решения и(у,х,/) от переменной у представляется одной базисной функцией (р(г)):
У
й(у, х, /) = а(х, ¿Мг?), Г1 = -Ц- (5)
Здесь предполагается, что <р{г]) удовлетворяет граничным условиям (2). Множитель а(х,£) находится из условия ортогональности невязки к весовой функции 1р(у,х,£). После подстановки (5) в (4), умножения на ф и интегрирования получающегося выражения в пределах [0, Щ приходим к уравнению для определения а(х,Ь):
н
У Цу)4>{у-1х,1)(1у = 0. (6)
о
[н
Введем в рассмотрение локальный расход жидкости д = / ийу. Интегрируя уравнение нераз-
рывности в пределах [0,/г], получаем уравнение
/о
Ы + Чх = о.
(7)
Искомую функцию а(х, можно выразить через ¡1, д:
') \
<р<1у\ ■ а{х,г) = д. (8)
Теперь из (6)-(8) получается система двух дифференциальных уравнений для определения функций Л,(ж,£), <?(ж,£).
Эта схема решения с одной базисной функцией впервые была реализована в [3], причем принималось 1р(у, х, ¿) = 1. В настоящей работе эта схема взвешенных невязок с одной базисной функцией используется для различных весовых функций ф(у, х, ¿). При полной реализации метода Галеркина-Петрова предполагается представление решения в виде ряда по системе базисных функций ср^ (в частности, при гр^ = ср^ получается метод Галеркина). К сожалению, реализация разложений с несколькими функциями ср^ приводит вследствие нелинейности задачи (1), (2) к очень громоздким уравнениям для а^, Л-ь Чк-, как эт0 продемонстрировано в [18, 19].
Применительно к задачам пограничного слоя подход с использованием одной базисной функции аналогичен методу интегральных соотношений. Эффективность метода, включая точность решений, наглядность полученных результатов в значительной степени зависят от выбора как базисной функции (р, так и интегрального соотношения, определяемого функцией ф.
Умножая первое уравнение системы (1) на ип (п = 0,1, 2,...) и интегрируя в пределах 0 ^ у ^ /г(ж,£), получим множество интегральных соотношений
А
п + 1
н н
1ип+Чу^ +(^!ип+2<1у о 1 о
н н
= Я J ип dy + ^ J ипиуу (1у, (9)
о о
в которых приняты обозначения <5 = 1 + Ьххх — еНх, А = 55.
В теории пограничного слоя соответствующие аналогичные уравнения рассматривались Т. Карманом и К. Польгаузеном (при п = 0), Л.С. Лейбензоном (при п = 1), В.В. Голубевым (при про-
1 (1р
извольном целом п = 2,3,...) с заменой ЫхЛ) на толщину пограничного слоя 5о(х) и О на ——
р ах
(см. [15, с. 565-566]).
Примем для продольной скорости и(у, ж, ¿) простейшее выражение, согласующееся с граничными условиями (2) и решением краевой задачи (1), (2) в стационарном случае для слоя постоянной толщины:
I у
и{у,х,г) = а(х,г)<р(г]), 1р(г]) = Г] - -Г]2, = (10)
Пользуясь (10), находим
н н 1
<7 = J ийу = ^ аН, а = иуу = — ^, J ип+1 г1у = ап+1Н J <рп+1 йг].
О 0 0
Введем обозначение / (рп(1г] = Фп. Теперь уравнение (9) можно записать в виде
./о
А
[Фп+1 {ап+1Ь)г + Фга+2 {ап+2К)Х] = Фга (апЪ) {Я ~ •
Проводя дифференцирование в скобках по £ и по ж и сокращая обе части уравнения на {апК),
получим
/ „ „2 \ ,
(П)
/¿- / \ П" /
_ п + 3(п + 2)ап Ап — —
п + 1
Ф0 = 1, =
Здесь
= вп\ я)н> п = 0,1, 2,3,...
= Зап, ап Фп+2 Вп = Фп
~ Фп+1 ' ЗФ„+1
2 1 15' Ф3 = 2 = -, Ф4 35' 8 ~ 315
При выводе уравнения (11) из системы уравнений пограничного слоя (1) по существу применяется метод Галеркина-Петрова взвешенных невязок по переменной у с одной базисной функцией аср(г]) и весовыми функциями (сир)п. Таким образом, каждое уравнение, которое следует из соотношения (11) при п = 0,1,2,..., получено методом взвешенных невязок при использовании одной базисной функции <р и весовых функций 1, 1р, 1р2,....
Полагая в (11) п = 0,1,2,..., получаем последовательность уравнений для q, к. Первые три
уравнения метода взвешенных невязок имеют следующий вид:
а) + = " =
б) Д^ + Н!*-!^) =!(«-£)», „ = 1; (12)
в) А + Щ- удх - ^ %]гх ) = ^ - к, п = 2.
9 h4X 3 h? V 9 V^ h?
При n = 0 система (7), (12,а) представляет базовую модель Шкадова [3]. Иногда ее называют моделью IBL (Integral Boundary Layer). Решение систем (7), (12,а) составляет фундамент современной теории волновых течений капиллярных пленок. При п = 1 получается модифицированная модель IBL, называемая также моделью WIBL. Уравнение (12,5) эквивалентно уравнению (12,а) и лишь незначительно отличается от него величинами числовых коэффициентов. Эта модель, ранее выведенная в [20], используется наряду с моделью IBL при решении конкретных задач. При п = 2 модификация модели IBL — уравнения (7), (12,в) — выведена впервые. Другие новые системы можно получить из (7), (11), полагая п = 3,4,.... Для обоснования преимуществ и недостатков трех моделей (7), (12) выполнены их численные решения при значениях управляющих параметров 5, е, а, соответствующих реальным нелинейным волнам. Далее проводится сравнение результатов.
3. Линейная неустойчивость. Сравним результаты исследования линейной неустойчивости, полученные на основе трех моделей (7), (12). Основное течение описывается формулами ho = 1, qo = 1. Полагая
h = 1 + heia{x~ct), q = 1 + qeia(x~ct)
и проводя линеаризацию в системах (7), (12), находим дисперсионные соотношения. Приравнивая к нулю определители возникающих линейных систем для h, q, получим
/ 2 12 6 \ iaA с--—с + -— 7 =3 — с, п = 0;
V 5 5 J
■ л ( 2 17 9 5 Л 5 ,0
гаА с--сН----7 = - (3 - с), п = 1;
V 7 7 6 7 6 v
д / о 22 4 7 \ 7
Здесь
л 9Д
7=д [а+е), е = nctgO, K =
Условием применимости систем (7), (12) является к2 <С 1. В экспериментах [1] с жидкостями малой вязкости значения 5 заключены в интервале [0,06-0,5], при этом к2 ~ 0,1. Для жидкостей повышенной вязкости в [16] выведена модельная система, которая обобщает систему (7), (12,а) и в которой учитываются члены о(п2). Там же даны оценки соответствующих интервалов допустимых значений 5 и к, Re и 7. Аналогичное обобщение с учетом членов о(п2) для системы (7), (12,5) проведено
в [18].
На рис. 1, а представлены нейтральные кривые, построенные для систем (7), (12) при п = 0,1, 2 для пленки, стекающей по вертикальной поверхности (в = 7г/2). Там же приведена нейтральная кривая 4, полученная прямым численным решением исходной линеаризованной задачи. Отметим, что кривая 1 соответствует модели п = 0, а кривая 2 — модели п = 1. О точности, обеспечиваемой этими моделями в линейном приближении, можно судить по степени отклонения кривых 1-3 от кривой 4- Легко видеть, что в интервале рассматриваемых конечных значений 5 модель п = 0
значительно точнее воспроизводит границу области неустойчивости. Имеет место общая закономерность: степень различия с точным решением монотонно растет с увеличением п при каждом 8. Отмеченная закономерность не выполняется вблизи критической точки нейтральной кривой 5 = 0, а = 0. В этой области различия между кривыми 1 3 и точным решением 4 малы и пригодна любая из моделей п = 0,1,2, причем наименьшее отличие от точного решения имеет место для модели п = 2. Однако эти свойства решений в области малых значений 8, а представляют интерес только с точки зрения теории. Дело в том, что волны в этой области малых значений 8, а имеют малую амплитуду и большую длину, проявляют свойства сильной неустойчивости и в экспериментах не фиксируются как четкие нелинейные регулярные волны. В работах [1, 2| вводилось критическое число Ые/г = 0,617~3/11, такое, что нелинейные волны наблюдаются и фиксируются при Ые > Кед.. Для воды при температуре 15°С это соответствует <5д. = 0,047. Данный вопрос исследовался в работе [21] на основе численного решения системы (7), (12,а). Установлено, что при малых 8 имеет место самопроизвольный распад первоначального длинноволнового возмущения и переход к последовательности нестационарных волн с меньшей длиной характерной волны. В [5] в результате анализа этих расчетов и значительного количества измерений выводится заключение, что существует такая граница атт, что нелинейные рмулярные волны наблюдаются только при а>о;т;п. При этом сктт слабо зависит от 8. В [22] на основе анализа системы (7), (12а) установлено, что при малых значениях волнового числа расстояние между главными горбами в периодической последовательности одиночных волн самопроизвольно сокращается вследствие гидродинамической неустойчивости. Найдено критическое значение а^, такое, что при а < ад. существуют бифуркации тина а —>■ 2а, а —>■ За. Отсюда следует, что физически значимые реальные волны, доступные физическим измерениям, существуют лишь на части области линейной неустойчивости, удовлетворяющей условиям 8 > 8ь, а > ехк- Из рассмотрения исключается, в частности, окрестность угловой точки нейтральной кривой 8 = 8С, а = 0.
Рис. 1. Нейтральные кривые для пленки, стекающей по вертикальной (в = 7г/2) поверхности (о): 1 п = 0: 2 1; 3 2: 4 решение уравнения Орра Зоммерфельда. и по наклонной к горизонту (в = 8°) поверхности (б): 1 п = 0: 2 1:3 2: 4 граница области неустойчивости /с: 5 граница области вторичной неустойчивости /2; б граница образования многогорбых волн /я: 7 частота наиболее
растущих по линейной теории волн /„,
На рис. 1,6 представлены нейтральные кривые 1 3, построенные для систем (7), (12) при п = 0,1,2 для пленки, стекающей но наклонной к горизонту поверхности (9 = 8°). Кривая /с представляет теоретическую границу, рассчитанную но линейной теории, и экспериментальную границу области неустойчивости [23]. Две другие кривые, проведенные но экспериментальным точкам, граница области вторичной неустойчивости /2 и граница образования многогорбых волн характеризуют нелинейные свойства наблюдаемых и измеряемых волн. Обе эти кривые начинаются и продолжаются вне окрестности угловой точки нейтральной кривой 8 = 8С, а = 0. Вне этой малой окрестности четко фиксируется основное свойство последовательности моделей (7), (12): наименьшее отличие от точной границы устойчивости /с обеспечивает модельная система п = 0, точность описания гра-
иицы неустойчивости другими модельными системами монотонно уменьшается с ростом п. Детали взаимного расположения нейтральных кривых 1 4 в окрестности 5 = ôc, а = 0 не представляют практического интереса, так как волны в этой зоне характеризуются малыми значениями амплитуды и коэффициента усиления и существенно неустойчивы. Экспериментальные точки /с, /2, fs, fm (частота наиболее растущих по линейной теории волн) расположены вне этой малой окрестности.
Обратимся теперь к сопоставлению нелинейных решений модельных систем (7), (11), соответствующих регулярным пространственно-периодическим волнам в стекающей пленке при конечных значениях 5.
Для численного исследования систем (7), (11) будут использоваться те же подходы, что и при изучении нелинейных волновых решений системы (7)* (12,а). Это касается прежде всего неединственности решений и связанных с этим бифуркаций регулярных волн и аттракторов нестационарных решений [24]. Подробный анализ бифуркаций и аттракторов исходной системы (1) при 5 —> 0 дан в работе [25]. Алгоритм отбора решений для сопоставления их с экспериментальными данными основан на использовании свойств оптимальности и неустойчивости волновых режимов, введенных в рассмотрение в [3, 26].
4. Установившиеся периодические и уединенные нелинейные волны. Рассмотрим краевую задачу на собственные значения, следующую из (7), (11), в предположении, что волны установившиеся:
h = h(0, Q = q(0, £ = x-ct.
Поскольку в этом случае из (7) следует q = qo + c(h — ho), из (11) при ho = 1 получаем уравнение для
h3 (h'" - eh' + 1) - q0 + c(l - h) =
= "7Г - с)2 + (2Bn - An)(qo - c)ch + (Bn - An + 1 )c2h2] h', ^
Un
где qo — постоянная интегрирования. Для периодической волны qo — средний по периоду расход в слое, для уединенной волны-солитона qo = 1 — толщина невозмущенного слоя.
Уравнение (13) имеет тривиальное решение h(x) = 1, qo = 1. В зависимости от краевых условий существуют также другие решения уравнения (13), в частности: в случае периодических волн с длиной волны \ = 2тт/а
h{ 0) = h{\), h'{ 0) = h'{\), h"{ 0) = h"{ A);
в случае уединенных волн-солитонов
£ ->■ ±00 : h^l, h' ->■ 0, h" ->■ 0.
Исследование установившихся периодических волн сводится к нелинейной задаче на собственные значения. Решение представляется конечным Фурье-разложением:
N
h= hk = -hk, /io = l.
k=—N
В этой задаче S, а — заданные внешний и внутренний свободные параметры; q, с — собственные значения. Методом Галеркина выводится система нелинейных алгебраических уравнений для hk, с, qo-
Для уединенных волн в силу того, что h = h — 1—>0 при >■ ± оо, уравнение (13) можно линеаризовать в окрестности стационарной точки (1,0,0):
h'" + nnh' + (3 - с)h = 0, Qn = [Вп{ 1 - С)2 + (2Вп - Ап){ 1 - с)с + (Вп -Ап + 1)с2] - е.
Характеристическое уравнение (14) имеет три корня: ai, 02, 03. Для быстрых положительных солитонов при с > 3
ai > 0, <Т2,з = —a + ib, а> 0.
Асимптотическое поведение при >■ ± оо заднего и переднего фронтов движущейся уединенной волны-солитона дается формулами
h = l + MeaiC,
, (15)
h = 1 + Ne-aCsm(b(), £ ->■ оо,
где = £ + {о, М, N, ^о — константы. При численном решении (13) строится профиль волны, согласующийся с асимптотическими представлениями (15) подбором констант М, N, с.
Для расчета нелинейных периодических регулярных волн и построения полной картины бифуркаций волновых семейств применялся метод инвариантного погружения [24]. При каждом фиксированном значении 5 и переменном параметре а существует несколько решений задачи (13) в виде регулярных волн. С увеличением 5 число таких решений быстро возрастает. Все семейства волн возникают с малой, но конечной амплитудой в точках бифуркации, расположенных вблизи значений ат = ас/т, т = 2,3,... , М, как почти гармонические волны с периодом 2тг/т и непрерывно продолжаются до s = 0, где s = а/ас [27].
Первое основное семейство 71 [3] начинается в точке Si = 1; волны этого семейства имеют наименьшую фазовую скорость среди всех других решений, существующих при каждом s в интервале 0 < s ^ 1. Второе основное семейство 72 [24] начинается в точке бифуркации sm', эти волны имеют наибольшую фазовую скорость при каждом s в интервале 0 < s < sm, М ^ 2. Существуют также промежуточные семейства медленных волн, точки бифуркаций которых s* расположены в интервале sm < s* < 1- При малых значениях 5 существуют только первое и второе основные семейства (М = 2) [24]. Промежуточные семейства появляются с ростом 5, причем М > 2 растет вместе с 5.
Волновые режимы, имеющие при фиксированных значениях 5, s наибольшие средний расход до, амплитуду hmax и фазовую скорость с, называются доминирующими.
Если зафиксировать 5, то получаем множество доминирующих волн в интервале 0 < s ^ 1. Это множество кусочно-непрерывное с разрывами в точках s^, в которых имеет место переход с одного семейства на другое. Как правило, sk ф s*, но — s*| — малая величина. В этих малых интервалах волновых чисел волны неустойчивы, наблюдаются двухпериодические и даже стохастические режимы. Как следует из определения, множество доминирующих волн при любом значении 5 включает часть волн первого 71 и второго 72 семейств, а также части промежуточных семейств, наблюдающихся при выбранном 5.
В плоскости (5, а) существуют линия I, выше которой и до а = ас имеется единственное решение — нелинейная одногорбая волна первого семейства, и линия II, ниже которой и до а = а* имеется нелинейная многогорбая волна второго семейства. В области между линиями I и II наблюдаются нелинейные волны промежуточных бифуркаций, количество которых растет с увеличением 5.
В качестве теста для сравнения моделей (7), (12) выберем зависимости с = c(hmax), где с — фазовая скорость установившейся волны, hmax — наибольшая высота главного горба на длине волны. На рис. 2, а представлены точки (с, hmах) для последовательности ж-периодических волн, порождаемых на вертикальной поверхности при 5 = 0,1; 9 = 7Г/2 изменениями волнового числа а в пределах ас ^ а > ск*. Маркерами 4 показаны результаты измерений с, hmax в пленках на вертикальной поверхности. Кривые 1-3 построены по результатам численных решений модельных систем (7), (12). Кривые А соответствуют режимам одногорбых волн первого семейства, возникающих при olc ^ а ^ ai, где I — линия первой бифуркации. Кривые В соответствуют режимам многогорбых волн, которые возникают при ai > а > ск*. Ближе всего к экспериментальным точкам расположены результаты расчетов по модели п = 0. Видно, что расчеты по моделям и = 1 и в = 2 не улучшают расчеты по базовой модели IBL: соответствующие кривые 2, 3 расположены дальше от экспериментальных точек.
На рис. 2, б представлены в плоскости (с, hmax) точки, соответствующие быстрым многогорбым волнам в пленке на наклонной поверхности при 5 = 0,226; 9 = 8°. Эти точки получены из экспериментов [23] и из численного решения модельных систем (7), (12) (кривые 13). Как и в случае волновых движений в пленке на вертикальной поверхности, на рис. 2, б прослеживаются следующие закономерности. К прямой, проведенной по экспериментальным точкам, ближе всего расположена кривая 1 (п = 0). Это означает, что модельная система (7), (12,а) дает наиболее точное решение исходной задачи (1), (2). Точность решений модельных систем (7), (12) уменьшается с увеличением номера п. В частности, заметно меньшую точность дает решение (7), (12,в) при п = 2.
Применим теперь модельные системы (7), (12) к исследованию уединенных волн-солитонов. Волновые режимы с волнами этого типа принципиально отличаются от режимов с периодической
последовательностью уединенных волн, рассмотренных выше. Основные теоретические результаты о свойствах волн-солитонов впервые получены в рамках подхода п = 0, 9 = ж/2 в [17, 24, 28]. Метод склеивания асимптотических разложений представлен в [17]. Там же указаны расчеты для медленных волн-солитонов (с < 3) и введен в использование параметр 8. Быстрые волны-солитоны (с > 3) при конечных значениях 8 построены в [29]. Основополагающие результаты этих работ получили к настоящему времени существенное развитие: ставится вопрос о том, в какой степени полезными для изучения волн-солитонов могут быть другие модели-модификации подхода 1ВЬ, представленные уравнением (13), при п = 1,2,....
Рис. 2. Зависимость фазовой скорости с от наибольшой высоты главного горба волны /?тах в пленке на вертикальной (6 = 0,1, в = -к/2) поверхности (а): 1 п = 0: 2 1: 3 2: 4 данные работы [5], и на наклонной к горизонту (6 = 0,226, в = 8°) поверхности (б): 1 п = 0: 2 1: 3 2: 4 данные работы [23]
Рис. 3. Зависимости фазовой скорости с одногорбой волны-солитона от параметра 6: 1 п = 0: 2 1: 3 2: 4 данные
На рис. 3 приведены зависимости фазовой скорости одногорбой волны-солитона от параметра 8, рассчитанные но указанному выше алгоритму решения краевой задачи (13), (15) при п = 0,1,2. Для сравнения маркерами показаны результаты, полученные прямым решением уравнений пограничного слоя [30]. Можно отметить хорошее соответствие расчетов уравнений пограничного слоя расчетам модельной задачи (13), (15) по методу 1ВЬ (п = 0). Имеется заметное отклонение от результатов теории пограничного слоя при использовании модифицировавших) метода 1ВЬ (п = 1). Это сравнение методов ранее было проведено в [31], однако схютвететвующего истолкования в [31] не обнаружено. Отклонение от точных результатов растет с увеличением номера модели (п = 2).
5. Нестационарное развитие волн и аттракторы. Обратимся теперь к вопросу о математическом моделировании множества нространственно-нериоди че-ских нестационарных волн, развивающихся из произ-
работы [30]
вольных малых начальных возмущений вследствие гидродинами ческой неустойчивости.
Рассмотрим задачу Коши о возбуждении нелинейных волн с помощью малых возмущений в начальный момент времени. Решение системы (7), (11) для нестационарной нелинейной ж-иериодической волны разыскивается в виде рядов Фурье, коэффициенты которых зависят от времени I:
/0М) =
N
Е
к=-М
лмх
Iк = ¡-к, / = (Л', д), Л'О = 1, С = а{х - с{1),
где а волновое число, С/ скорость соответствующей линейной волны.
После подстановки (16) в (11) методом Галеркина формируется динамическая система для /г,д., дк- Для вычисления нелинейных правых частей уравнений применяются псевдоспектральный метод и процедуры быстрых Фурье-преобразований, как в работах [27, 32]. В начальный момент задаются условия
¿ = 0: 1гк = 1гк( 0), Як = Як{ 0).
Динамическая система интегрируется численно до установления асимптотического предельного решения. Как правило, в начальный момент задается малая амплитуда первой гармоники /г. Каждая пара значений 6, е порождает серию вариантов с параметром а € [а*,ас].
Для вычисления фазовых скоростей гармоник нелинейной волны использовалось соотношение
1с1, 1т Д.
^ — ^ _ _ _ с1ГС1-|" _—
ка (И Ке/д.'
которое получается из (16) дифференцированием но I и по х. Каждое допустимое значение а соответствует определенной временной частоте возмущений, которую можно вычислить по формуле
1 ír¡2\l/'i / = —кН' 6"7г V V )
Основные известные к настоящему времени результаты по нестационарному развитию пространственно-периодических волн и структуре аттракторов получены в численных решениях системы (7), (12,а) метода 1ВЬ [27, 32].
На основе систематических численных экспериментов для фиксированных значений внешних параметров 8, е установлено существование при а* < а < ас глобального аттрактора ж-периодиче-ских решений системы (7), (11), составленного из частей множества доминирующих волн. При любых малых гармонических или стохастических возмущениях основного состояния /г = 1, д = 1 предельными при I —> оо состояниями оказываются, как правило, доминирующие волны. Дополнительно к ним обнаружены осциллирующие режимы на малых интервалах значений а в окрестностях точек бифуркаций. На этих интервалах значений а множество предельных доминирующих воли терпит разрывы. Таким образом, проекция глобальншх) аттрактора на плоскость (с, а) представляет собой кусочно-непрерывную кривую, причем куски ее принадлежат различным семействам рехуляр-ных волн, появляющимся при изменении а. Эти семейства возникают в точках бифуркации по мере уменьшения а: первое семейство 71; промежуточные семейства, количество которых увеличивается с ростом 6; второе семейство 72. Периодические одиночные волны, наблюдаемые в экспериментах, соответствуют решениям вторшх) семейства 72 •
Рис. 4. Глобальный аттрактор (о расход </о- б фазовая скорость с) при 6 = 0,1: в = и/2:
1 п = 0: 2 1: 3 2
На рис. 4 в качестве примера показан глобальный аттрактор кусками сплошной линии при 8 = 0,1; 9 = ж/2. Продвижение по ветви второго семейства с уменьшением а соответствует переходу к периодическим последовательностям одиночных волн с возрастающим расстоянием между их
главными горбами. Фазовая скорость волн монотонно растет с уменьшением а. На рис. 2 представлены точки (с, 6), относящиеся именно к этой части глобального аттрактора, при 6 = 0,1 и 6 = 0,226 соответственно. Средний расход до вначале возрастает с уменьшением а, однако имеется точка а*, в которой до достигает максимума и затем начинает уменьшаться. В расчетах нестационарных решений при п = 0 обнаружен новый эффект, состоящий в том, что при малых значениях а < а* притягивающие свойства второго семейства ослабляются и глобальный аттрактор в этой области значений а разрушается. Наблюдается самопроизвольный переход от волны с длиной 2тт/а к более коротким волнам с длинами 2ж/2а, 2ж/3а [27].
Подробное численное исследование нестационарных волн и аттракторов в пленке на твердой поверхности при углах наклона 9 = 90; 8; 6,4; 4,6° и различных значениях 6 выполнено в [33] на основе решения системы (7), (12,а). Там же даны примеры успешного сравнения проведенных расчетов и экспериментов [23, 34].
Рассмотрим теперь некоторые результаты о нестационарных волнах и аттракторах, полученные на основе численных решений систем (7), (11) при п = 0,1,2. На рис. 5 показаны в плоскостях (до, а), (с, а) следы глобальных аттракторов для случая 6 = 0,226, 9 = 8°. На первом интервале непрерывности, соответствующем первому основному семейству 71, все три кривые качественно правильно воспроизводят зависимость до (ск). Однако точность решений монотонно падает с ростом п относительно варианта п = 0. Самая большая потеря точности имеет место при переходе от варианта п = 0 к варианту п = 1. На интервале промежуточных бифуркаций и двухпериодичееких волновых движений кусочно-непрерывные части глобального аттрактора заметно завышены по сравнению с результатами модели п = 0. Расчеты при п = 1, п = 2 требуют больших затрат времени для разделения нестационарных двухпериодичееких и стохастических движений. Наконец, на интервале второго основного семейства 72 модели п = 1, п = 2 значительно завышают результат модели п = 0, точность заметно падает с ростом номера. Сложную для расчета область, включающую локальный максимум до при малом а, они не воспроизводят в противоположность модели п = 0. Отметим, что зависимости с(а) качественно одинаковы для всех моделей п = 0,1,2 и слабо различаются коли чеетвенно.
Рис. 5. Глобальный аттрактор (о расход </о, б фазовая скорость с) при (5=0,229: в = 8°:
1 п = 0: 2 1: 3 2
Таким образом, модели п = 1, п = 2 не открывают никаких новых возможностей при вычислении нестационарных волн и глобального аттрактора по сравнению с испытанной моделью п = 0, более того, повторяют эту модель в качественном отношении при заметных количественных отклонениях. Степень отклонений от модели п = 0 увеличивается с ростом п.
6. Заключение. Выведена последовательность модельных систем эволюционных уравнений для локальной толщины пленки ]г(х,1) и расхода д(х,1). К уравнениям пограничного слоя применен метод Галеркина Петрова но переменной у = 1гг] с использованием базисной функции и=г]—г]2/2 и системы функций ип(у,х,1) для ортогонализации невязки. Выведенные модели идентифицируются но значению номера п, независимы между собой, и каждая из них может применяться к исследованию волновых режимов в пленке.
Получены численные решения этих систем при п = 0,1,2, и по ним вычислены такие характеристики нелинейных пространственно-периодических и уединенных волновых структур, которые
можно сопоставить между собой и сравнить с экспериментальными данными либо с точными численными решениями соответствующих краевых задач для уравнений Навье-Стокса.
Показано, что во всех рассмотренных случаях нелинейных волн наилучшее соответствие с экспериментальными и точными результатами дает система с параметром п = 0 (метод IBL). Все последующие модельные системы начиная с п = 1 (метод WIBL) дают лишь качественное соответствие с результатами по методу IBL при значительных количественных ошибках. Как правило, величина ошибки во всех рассмотренных задачах растет монотонно с увеличением п при п = 1,2,....
Работа выполнена при финансовой поддержке РФФИ, грант № 15-014)5186.
СПИСОК ЛИТЕРАТУРЫ
1. Капица П.Л., Капица С.П. Волновые течения тонких слоев вязкой жидкости // Журн. эксперим. и теор. физ. 1949. 19, № 2. 105-120.
2. Капица П.Л. Волновые течения тонких слоев вязкой жидкости // Журн. эксперим. и теор. физ. 1948. 18, № 1. 3-28.
3. Шкадов В.Я. Волновые режимы течения тонкого слоя вязкой жидкости под действием силы тяжести // Изв. АН СССР. Механ. жидкости и газа. 1967. № 1. 43-51.
4. Холпанов Л.П., Шкадов В.Я. Гидродинамика и теплообмен с поверхностью раздела. М.: Наука, 1990.
5. Алексеенко С.В., Накоряков В.Е., Покусаев Б.Г. Волновое течение пленок жидкости. Новосибирск: Наука, Сиб. отд-е, 1992.
6. Chang Н.-С., Demekhin Е.А. Complex Wave Dynamics on Thin Films. Amsterdam: Elsevier, 2002.
7. Nepomnyashchy A. A., Velarde M.G., Colinet P. Interfacial phenomena and convection. Boca Raton: Chapman & Hall/CRC, 2002.
8. Kalliadasis S., Ruyer-Quil C., Scheid В., Velarde M.G. Falling Fiquid Films. Fondon: Springer, 2011.
9. Зейтунян P.X. Проблема термокапиллярной неустойчивости Бенара-Марангонп // Успехи физ. наук. 1998. 168, № 3. 259-286.
10. Shkadov V Y. Hydrodynamics of slopped falling films // Interfacial phenomena and the Marangoni effect / Ed. by M. Velarde, R.K. Zeytounian. New York: Springer-Verlag, 2002. 191-224.
11. Шкадов В.Я., Демехин Е.А. Волновые движения пленок жидкости на вертикальной поверхности (теория для истолкования экспериментов) // Успехи механики. 2006. 4, № 2. 3-65
12. Shkadov V. Y., Sisoev G.M. Wavy falling films: theory and computations instead of physical experiment // IUTAM Symposium on Nonlinear Waves in Multi-Phase Flow. Vol. 57 / Ed. by H.-C. Chang. Notre Dame. USA; Dordrecht: Kluwer, 2000. 1-10.
13. Shkadov V. Ya., Velarde G.M., Shkadova V.P. Falling films and the Marangoni effect // Phys. Rev. E. 2004. 69, 056310.1. 1-15.
14. Trifonov Y. Y. Stability and bifurcations of the wavy film flow down a vertical plate: the results of integral approaches and full-scale computations // Fluid Dyn. Res. 2012. 44. 031418 (19pp).
15. Кочин H.E., Кибель И.А., Розе H.B. Теоретическая гидромеханика, II. М.: ГИФМЛ, 1963.
16. Шкадов В.Я. Двухпараметрическая модель волновых режимов течения пленок вязкой жидкости // Вестн. Моск. ун-та. Матем. Механ. 2013. № 4. 56-61.
17. Шкадов В.Я. Уединенные волны в слое вязкой жидкости // Изв. АН СССР. Механ. жидкости и газа. 1977. № 1. 63-66.
18. Ruyer-Quil С., Manneville P. Improved modeling of flows down inclined planes // Eur. Phys. J. B. 2000. 15. 357-369.
19. Актершев С.П., Алексеенко С.В. Модель волнового течения стекающей пленки вязкой жидкости // Прикл. механ. и техн. физ. 2013. 54, № 2. 21-31.
20. Ruyer-Quil С., Manneville P. Modeling film flows down inclined planes // Eur. Phys. J. B. 1998. 6, N 2. 277-292.
21. Демехин E.A., Токарев Г.Ю., Шкадов В.Я. О существовании критического числа Рейнольдса для стекающей под действием веса пленки жидкости // Теор. основы хим. технол. 1987. 21, № 4. 555-559.
22. Сисоев P.M., Шкадов В.Я. Развитие доминирующих волн из малых возмущений в стекающих пленках вязкой жидкости // Изв. РАН. Механ. жидкости и газа. 1997. № 6. 30-41.
23. Liu J., Gollub J.P. Solitary wave dynamics of film flows // Phys. Fluids. 1994. 6. 1702-1712.
24. Вунов А.В., Демехин E.A., Шкадов В.Я. О неединственности нелинейных волновых решений в вязком слое // Прикл. матем. и механ. 1984. 48, № 4. 691-696.
25. Demekhin Е.А., Tokarev G.Y., Shkadov V. Y. Hierarchy of bifurcations of spaceperiodic structures in a nonlinear model of active dissipative media / / Physica D. 1991. 5, N 2. 338-361.
26. Шкадов В.Я. К теории волновых течений тонкого слоя вязкой жидкости // Изв. АН СССР. Механ. жидкости и газа. 1968. № 2. 20-25.
27. Сисоев Г.М., Шкадов В.Я. Доминирующие волны в стекающих пленках вязкой жидкости // Докл. РАН. 1997. 357, № 4. 483-486.
28. Бунов А.В., Демехин Е.А., Шкадов В.Я. Бифуркации уединенных волн в стекающем слое жидкости // Вести. Моск. ун-та. Матем. Механ. 1986. № 2. 73-78.
29. Демехин Е.А., Шкадов В.Я. К теории солитонов в системах с диссипацией // Изв. АН СССР. Механ. жидкости и газа. 1986. № 3. 91-97.
30. Chang Н.-С., Demekhin Е.А., Kalaidin Е. Simulation of noise-driven wave dynamics on a falling film // AIChE J. 1996. 42, N 6. 1553-1568.
31. Ruyer-Quil C., Manneville P. On the speed of solitary waves running down a vertical wall //J. Fluid Mech. 2005. 531. 181-190.
32. Демехин E.A., Токарев Г.Ю, Шкадов В.Я. Двумерные нестационарные волны на вертикальной пленке жидкости // Теор. основы хим. технол. 1987. 21, № 2. 177-183.
33. Тушканов Д.А., Шкадов В.Я. Нелинейные волны в пленке жидкости на почти горизонтальной поверхности // Изв. РАН. Механ. жидкости и газа. 2006. № 3. 11-24.
34. Liu J., Paul J.D., Gollub J.P. Measurements of the primary instabilities of film flows //J. Fluid Mech. 1993. 250. 69-101.
Поступила в редакцию 25.11.2016
УДК 531.396
О ПРИМЕНИМОСТИ МЕТОДА ШАТРОВ В ЗАДАЧЕ СТАБИЛИЗАЦИИ ПЕРИОДИЧЕСКОГО ПРОГРАММНОГО ДВИЖЕНИЯ
К. А. Тузов1
Рассматривается задача поиска оптимальных параметров стабилизации линейных периодических систем. Системы приводятся к стационарному виду, после чего применяется метод шатров В.Г. Болтянского. Работоспособность процедуры демонстрируется на примере.
Ключевые слова: метод шатров В.Г. Болтянского, мультипликаторы Флоке, матричное уравнение Ляпунова.
The search problem for optimal parameters is considered to stabilize linear periodic control systems. Such a system is reduced to a stationary form; then, Boltyanski's tent method is applied. The efficiency of this procedure is illustrated by an example.
Key words: Boltyanski's tent method, Floquet characteristic multipliers, Lyapunov matrix equation.
1. Введение. Метод шатров В.Г. Болтянского [1] является одним из основных методов решения задач экстремального рода. В настоящей работе рассматривается возможность его применения для оптимальной стабилизации периодических линейных систем.
2. Постановка задачи. Задача, исследуемая в работе, появилась в результате обнаружения ошибок при реализации алгоритма имитации невесомости на центрифуге с управляемым кардано-вым подвесом, более подробно описанного в книге [2]. В течение всего этапа имитации орбитального полета консоль центрифуги вращается с угловой скоростью
0J = Wo + sin(2"7rz/t), (1)
где Wo, ш, v — некоторые константы, а внешнее полукольцо и кабина поворачиваются по заданным периодическим законам. Благодаря этому алгоритму происходит изменение модуля и направления
1 Тузов Кирилл Александрович — студ. каф. прикладной механики и управления мех.-мат. ф-та МГУ, e-mail: ktuzovQgmail .com.
