Вязкие напряжения как источник вихревой неустойчивости в слое поверхностно-активного вещества

Интенсивное изучение процессов на границе раздела сред в последние десятилетия обусловлено как широким распространением многофазных систем в природе, так и их высокой востребованностью в технологических приложениях [1 –4] . Вклад поверхностных сил становится определяющим, когда элементы жидкой системы, из которых она состоит, имеют малый характерный размер, например, пузырьки, капли, плёнки и другое. Интерес к управлению такими системами стремительно растёт в связи с развитием технологий проектирования микрожидкостных устройств [5, 6] , платформ «Lab-on-a-Chip» [7 –9] , проточных химических микрореакторов [10 –12] . Характерный размер перечисленных устройств, как правило, не превышает одного миллиметра. На таких масштабах капиллярные эффекты доминируют над объёмными силами, и межфазная граница является не просто пассивным геометрическим барьером, а становится полноценным активным генератором движения жидкости. Описание поверхностной динамики существенно усложняется в присутствии поверхностно-активных веществ (ПАВ). Молекулы ПАВ модифицируют граничные условия для скорости течения, поскольку создают дополнительные напряжения. Перераспределение напряжений в пространстве порождает градиент поверхностного натяжения и, как следствие, дополнительную касательную силу. При балансе напряжений на границе раздела её движение прекращается, а граничное условие скольжения сменяется условием прилипания [13] . На практике этот эффект приводит к стабилизации конвективных неустойчивостей [14, 15] и замедлению движения включений [16 –18] , что находит применение при решении различных технологических задач. Если ПАВ растворимы в жидкости, динамика дополнительно усложняется обменом вещества через границу раздела с основным объёмом жидкости (субфазой) за счёт адсорбции/десорбции. В этом случае напряжение на поверхности определяется балансом между поверхностной адвекцией и обменной кинетикой ПАВ [3] . Также ситуацию усугубляет способность слоя ПАВ претерпевать фазовые переходы. При этом основной характеристикой служит поверхностное давление Π = σ ₀ - σ , определяющее величину снижения поверхностного натяжения от значения σ 0 для чистой поверхности до значения σ при наличии на ней плёнки.

Управляющим параметром фазового перехода является число упругости E, равное отношению поверхностного давления Π к полному касательному напряжению, приложенному к плёнке [19] . Наконец, ПАВ обладают собственной поверхностной реологией, которая также изменяет граничные условия. Деформация слоя ПАВ вызывает сочетание поверхностных вязких (сдвиговых и дилатационных) напряжений с напряжениями Марангони, что также увеличивает число параметров системы и расширяет спектр возможных режимов движения [20 –23] .

Ввиду разнообразия и сложности механизмов, управляющих динамикой слоя ПАВ, теоретические исследования зачастую идут по пути упрощения геометрии. Удобной с математической точки зрения моделью служит двумерное течение в объёмной фазе с одномерной межфазной границей. Недавняя работа [24] , посвящённая изучению влияния частично загрязнённой границы раздела на термокапиллярный поток в двумерном щелевом канале, показала, что взаимодействие течения со слоем нерастворимого ПАВ может быть исчерпывающе описано единственным безразмерным параметром — числом упругости E. При 0 1 на границе раздела формируются две зоны: свободная от ПАВ область с интенсивным термокапиллярным потоком и так называемая застойная зона, где касательные напряжения термо- и концентрационно-капиллярного происхождения взаимно уравновешиваются. В работах [19, 25] получено аналитическое решение плоской двумерной задачи течения жидкости в прямоугольной полости со свободной поверхностью, покрытой плёнкой ПАВ, в условиях фиксированного градиента температуры. Уравнение состояния плёнки, связывающее поверхностное давление П с концентрацией ПАВ, предполагалось линейным. На поверхности также выявлены две характерные зоны — свободная и застойная. При определённых допущениях авторам удалось вывести аналитические выражения, в частности, формулу для положения границы между зонами как функции параметра E. Полученные результаты подтверждены прямым численным моделированием в ряде исследований, а именно [26 –28] . В работе [29] получены новые результаты для случая, когда ПАВ в пределах границы одновременно существует в двух фазовых состояниях. Результаты также засвидетельствовали, что одна часть поверхности свободна от ПАВ и подвержена интенсивному термокапиллярному течению, а другая является застойной, где скорость поверхности примерно на два порядка ниже.

Структура поверхностного течения качественно меняется при повышении размерности задачи — при переходе к трёхмерному объёмному течению с двумерной поверхностью. В этом случае появляется возможность нарушения симметрии в изначально однородном объёмном течении, ограниченном поверхностью со слоем ПАВ. В результате на границе раздела возникает периодическая последовательность вихрей [30] . Подобные структуры неоднократно наблюдались в однородном плоском течении в декартовой [31 –34] и цилиндрической [35 –37] геометриях, однако удовлетворительного объяснения механизма их образования долгое время не существовало. Например, в ряде теоретических работ [38, 40] авторы полагали, что имеют дело с чистой поверхностью, а явление вызвано неустойчивостью объёмного течения. Однако эксперименты показали, что очистка поверхности приводит к симметричному течению в объёме, устойчивому во всём исследуемом диапазоне скоростей [30] . Добавление же молекул ПАВ неизменно вызывало нарушение симметрии. При достаточно больших мощностях источника объёмного течения слой ПАВ сжимался, и на поверхности формировалась двухзонная структура, аналогичная наблюдавшейся в экспериментах с одномерной границей. Под свободной поверхностью течение оставалось симметричным, тогда как при наличии дополнительного пространственного измерения в застойной зоне образовывались вихри. Как и в двумерном случае, сжатие слоя, образование двух зон и положение застойной линии определяются значением числа упругости [29] . Тем не менее данный параметр не объясняет причины возникновения вихревой структуры.

Вероятно, основы теоретического анализа кризиса равновесия в такой системе были заложены в работе [40] при исследовании плоского течения Куэтта, сдвигающего слой ПАВ. Рассматривалась устойчивость течения под поверхностью, касательные напряжения на которой задавались градиентом поверхностного натяжения. Дополнительно принималось предположение, что поверхность является плоской и несжимаемой. При этом динамическое влияние ПАВ сводилось к действию потенциальной силы и не могло объяснить вихревую составляющую движения, условия возникновения которой фактически описывались уравнением чистой поверхности. Авторы пришли к выводу о беспороговом развитии неустойчивости, что противоречило эксперименту. Усовершенствованная теория предложена в [41] , где впервые удалось получить пороговую неустойчивость с отличным от нуля критическим волновым числом. Локализация движения вблизи границы позволила провести сравнительно простой анализ при рассмотрении поверхности как самостоятельной гидродинамической сущности. Для этого авторы сформулировали феноменологическое замыкание посредством условий проскальзывания, накладываемых на объёмную и поверхностную скорость. Дивергентная составляющая движения, связанная с градиентом поверхностного давления, учитывалась за счёт допущения сжимаемости слоя ПАВ. Модель дала возможность впервые получить порог неустойчивости слоя. Во всём диапазоне волновых чисел мнимая часть инкремента оказывалась отличной от нуля. Это говорило об осциллирующем характере возникающей структуры и свидетельствовало в пользу эксперимента [32] . Тем не менее критическая скорость набегающего потока на три порядка превышала значения, полученные при натурных испытаниях. Это несоответствие явно указывало на недостаточность модели.

Наконец, в недавних экспериментах в осесимметричной геометрии с различными генераторами радиального движения было показано, что при достижении критической мощности источника в слое нерастворимого ПАВ развивается двумерная вихревая конвекция [42]. Поверхностное течение также имеет двухзонную структуру: осесимметричную дивергентную область вблизи источника и азимутально-периодическое вихревое движение на периферии. Как и ранее, разделение поверхности на две характерные области вызывает сжатие слоя ПАВ, когда касательные напряжения со стороны источника превышают поверхностное давление в слое. В результате слой ПАВ сносится к периферии, формируется кольцевая застойная зона, ширина которой определяется параметром упругости E. Для выяснения роли объёмного течения в формировании вихрей эксперимент предусматривал ограничение движения жидкости различными барьерами. Оказалось, что при установке барьера на пути объёмного течения поверхность полностью теряет подвижность. Если же барьер соприкасался с жидкостью со стороны газовой фазы, на поверхности формировались вихревые структуры. Вихри всегда возникали только в застойной зоне межфазной границы и исчезали по мере распространения вглубь объёмной фазы. Таким образом, убедительно продемонстрировано, что природа неустойчивости не связана с напряжениями Марангони. Единственным источником движения в застойной зоне является поток импульса от течения в объёме, обусловленный вязкими напряжениями.

Итак, целью настоящей работы является аналитическое и численное исследование вихревой неустойчивости в слое поверхностно-активного вещества. Обсуждается возможный физический механизм её развития. Предполагается, что застойная зона при равновесном состоянии системы «слой ПАВ–жидкость» при объёмном течении играет роль твёрдой стенки. При достижении касательными напряжениями критического значения поток импульса через границу раздела претерпевает скачок, в результате этого слой ПАВ приходит в движение. Первая часть работы посвящена математической постановке задачи: выводу уравнений движения слоя ПАВ на основе модели Буссинеска–Скривена с их последующим усреднением поперёк слоя объёмной фазы. Далее излагаются результаты линейного анализа устойчивости: находится дисперсионное соотношение, определяются нейтральные кривые, исследуется зависимость критических параметров от чисел Рейнольдса и Буссинеска. Затем приводятся результаты прямого численного моделирования полной нелинейной системы уравнений, описывается эволюция вихревых структур, их динамика во времени и пространстве. В заключение результаты суммируются, обсуждается их соответствие экспериментальным данным.

• 2.    Математическая формулировка задачи

Рассмотрим несжимаемую однородную жидкость с динамической вязкостью η и плотностью ρ. Движение жидкости ограничено твёрдой стенкой снизу и слоем нерастворимого ПАВ сверху. Свяжем с данной системой декартовы координаты x, y, z и будем считать, что слой ПАВ расположен в плоскости z = 0. Расстояние между граничными поверхностями обозначим h . В условиях, когда поверхность находится в состоянии механического равновесия, в объёмной фазе имеет место заданный плоскопараллельный поток U ₀ = U 0 (z) e x , направленный вдоль оси x (Рис. 1) . Тогда уравнения возмущённого движения объёмной жидкой фазы имеют окончательный вид:

^d v +U o (z) | v + v z dU^ + tv ^V ) v = - ¹ V p+v V 2 v ,                   (1)

∂t       ∂x dz              ρ

V^v = 0, где v = v(x,y,z) — вектор поля скорости, p — давление, v = n/p — кинематическая вязкость. На нижней твёрдой границе выполняется условие прилипания:

z = — h : v = 0.

( a )

Рис. 1. Схематическое изображение области течения под слоем ПАВ: слой в состоянии механического равновесия, характерная скорость U течения под поверхностью ниже критической U c ( а ); слой в состоянии двумерной вихревой неустойчивости, U>U c ( б )

Обсудим подробнее граничные условия, которым должны удовлетворять скорости и силы на границе со слоем ПАВ. В работах [41, 42] предложено считать слой на поверхности жидкости отдельной фазой. В соответствии с данным предположением имеет смысл рассматривать поверхность как самостоятельный объект и моделировать её собственным уравнением движения, включающим связь с полем основного течения в объёме. Буссинеск [20] был первым, кто предложил модель вязкого сопротивления движению поверхности для объяснения аномальной скорости осаждения капель. Позднее Скривеном [21] сделано обобщение данной модели, в которую межфазная граница включена как отдельная двумерная ньютоновская жидкость, обладающая собственными реологическими свойствами.

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

∂ u

Psi _d +( u •^ S ) u

-V s П + n s V2_S u —n^ ∂z

z =0

где p s , П , u (x,y) — поверхностные плотность, давление и скорость, V s — оператор набла, действующий вдоль поверхности, V | = (v x ,v y ) — проекция вектора скорости на плоскость слоя, n s — сдвиговая поверхностная вязкость. При разумной характерной скорости объёмного потока в сжатой застойной зоне нерастворимые монослои фактически являются несжимаемыми, а ширина застойной зоны остаётся неизменной [42] . Таким образом, можно игнорировать переходный режим, связанный со сжатием слоя, и считать, что

VS • и = 0.

Уравнение (2) по сути представляет собой динамическое граничное условие, которому удовлетворяет течение в объёме. Кроме того, непрерывность поля скоростей на поверхности требует задания кинематического условия вида:

z = 0:   v|| = u ,      v z = 0.

Как правило, система (1) - (4) , наряду с p = const, замыкается уравнением состояния, связывающим поверхностную плотность ρ S с давлением Π . Однако эффекты, связанные с поверхностным натяжением и его зависимостью от концентрации, не могут повлиять на вихревую моду неустойчивости, а дивергентная мода исключается уравнением (3) . К тому же существуют экспериментальные подтверждения, что на стадии формирования вихрей в ходе кризиса равновесия сжатой застойной зоны концентрация ПАВ практически остаётся постоянной [42] , и запись уравнения переноса вещества вдоль поверхности не требуется (p s = const).

Приведём систему к безразмерному виду. Выберем в качестве единиц: расстояния — толщину слоя h, скорости — характерную скорость течения в объёмной фазе U, времени — h/U, объёмного давления — ρU ² , поверхностного давления — U η _S /h. Тогда уравнения (1) и (2) перепишутся в виде:

^{d V} +U o ^{d V} +v z dU ⁰ + ( v ^V ) v = -V p+^ V ² v , (5) di dx dz Re

∂ u

^Ps ^Re di ^{+( u} '^V s ^{) u} =

-V s П + V ² u --1

Bq dz

.

Сюда входят безразмерные параметры:

P s =   = -,

η S    ν S

η S

Bq= т, ηh

Uh

Re= — ν

.

Число Рейнольдса Re является основным управляющим параметром задачи и определяет скорость движения жидкости в объёме. Величина Р s представляет собой поверхностный аналог числа Прандтля для неизотермических течений и отражает отношение характерных времен диффузии импульса в объёме и на поверхности. Похожую роль выполняет и число Буссинеска Bq. Как P s , так и Bq имеют постоянные значения и используются для описания геометрии течения, поверхностной реологии ПАВ и реологических свойств субфазы.

Очевидно, что (6) является не просто пассивным граничным условием для уравнения (5) , а выступает полноценным динамическим уравнением, представляющим эволюцию слоя ПАВ. Совместное рассмотрение двумерной динамики в плоскости z = 0 и трёхмерной динамики субфазы сталкивается с существенными вычислительными трудностями. Поэтому в случае трёхмерных течений исследования ограничиваются конкретным классом движений с известной симметрией [43] . Поле u (x,y) рассматривается как полностью самостоятельное, связанное со скоростью субфазы v динамическим соотношением. Поскольку ключевые для анализа неустойчивости процессы локализованы вблизи межфазной границы и определяются вектором V | (x,y ,z), основное внимание уделяется далее двумерному течению. Зависимость последнего от поперечной координаты z задаётся явным образом в виде:

v (x,y,z) = U o (z) e x +€(z) v s (x,y)+ 0(z) u (x,y),     V z = 0,                          (7)

где U0(z) — профиль основного плоскопараллельного течения в объёме, в качестве которого далее будет пониматься следующая парабола:

U o (z) = —4z(z +1).

Объёмной модой ξ(z) будем считать первую собственную функцию оператора вязкой диффузии, удовлетворяющую условиям прилипания на обеих границах. Такой выбор при исследовании неустойчивости в объёме не был бы оправдан, но, согласно выводам из работы [42], объёмное течение остаётся устойчивым во всём диапазоне значений числа Re, в котором наблюдается поверхностная неустойчивость. На функцию 0(z) следует наложить дополнительное условие быстрого затухания вглубь объёма субфазы. Так как неустойчивость целиком определяется вязким переносом импульса вблизи поверхности, характерный масштаб затухания свяжем с масштабом возмущений пограничного слоя: 6 ^ Re-1/3. Данная оценка следует из анализа уравнения Эйри, к которому сводится линеаризованное уравнение для возмущений в сдвиговом потоке вблизи границы [44]. Такой выбор является модельным допущением, правильно передающим характерный масштаб затухания возмущений и позволяющим сохранить аналитическую структуру модели. Задача модели состоит не столько в воспроизведении трёхмерных эффектов, сколько в получении минимальной замкнутой системы, качественно согласующейся с экспериментом.

Итак, имеем:

£(z)=sin(n(z +1)),     0(z) = e z/ cos(z/6).                               (8)

Подставим (7) , (8) в уравнения движения и найдём проекции уравнения (5) на ξ(z), а уравнения (6) на θ(z):

( ... ) ^(Xf (z)dz,

- 1

где f = { £,# } . Данная процедура даёт возможность эффективным образом учесть взаимное влияние полей u и v друг на друга, избежав необходимости прямого удовлетворения динамическим граничным условиям. Нетрудно убедиться, что интегралы вида ( £0 ) , ( U 0 £0 ) , ( £'20' ) и ( £0 2 ) оказываются более чем на порядок меньше остальных. Слагаемыми, содержащими их как коэффициенты, можно пренебречь. Кроме того, будем считать, что импульс из объёма к поверхности передаётся исключительно через x-компоненту возмущения скорости. В состоянии равновесия поверхность неподвижна, а касательное напряжение компенсируется градиентом поверхностного давления. Это положение удерживается до тех пор, пока профиль скорости сохраняет свою форму. При наложении возмущения форма профиля вблизи поверхности изменяется. Однако поверхность приходит в движение не при любом возмущении, а лишь тогда, когда кривизна профиля скорости v x (z) у поверхности достигает критического значения и равновесие касательных сил нарушается. Данное условие затрагивает исключительно x-компоненту скорости, поскольку именно она определяет продольное касательное напряжение. Возмущение v y , хотя и порождает поперечное напряжение, не влияет на кривизну профиля v x (z) и, следовательно, не приводит к кризису равновесия. Таким образом, лишь v _x способна инициировать движение поверхности. С учётом всего вышесказанного получаем:

d v S _ _,            _ d v_s    _     1

д^ +N(vs •V'yvs + Li dx^ = —Vp + R (^ vs + L2vs + L3u)’ где

∂u di+(u -Vs)u

-Vs n + VS u -

G i V x + G 2 U X Bq

e x

( £ 3 )     r - ( U o £ )     т         / £

^N ( £ 2 ) ’      ¹     ( £ 2 ) ’      ² ( £ ² )\ dz ²

L 3 =

G 1 =

G 2 =

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

∂Ψ        ∂Ψ       ∂ψ∂ψ

= dy’   v y =^- dx,   u x = dy’   ^{U y} = ^-a?

Здесь Ψ и ψ — функции тока объёмной и поверхностной фаз. Далее, применяя к (9) , (10) операцию взятия ротора, получаем:

tx

PsRef дф + j(ф,^\ = v2ф+ф d2(G1*+G2,   ф = -V4

\ dt        J        Bq dy ²

где Ф = ( Vx v ) _z , ф = ( Vx u ) — вихри скоростей в объёме и на поверхности, J (f,—) =    g — —    gr— — якобиан.

∂x ∂y ∂y ∂x

Уравнения (11) , (12) описывают поля возмущённых объёмного и поверхностного течений.

• 3.    Устойчивость равновесия

Рассмотрим вопрос об устойчивости равновесного состояния поверхности при u = 0, v = U ₀ (z) e _x . Очевидно, что в этом случае на поверхности должно выполняться условие баланса касательных напряжений:

dn ₀ _   1 dU ₀

dx Bq dz выражающее зависимость равновесного поверхностного давления от координаты в направлении течения. Можно видеть, что П0 (x) — линейная возрастающая функция продольной координаты. В основном состоянии слой ПАВ находится в покое, поскольку вязкие напряжения со стороны субфазы уравновешиваются градиентом поверхностного давления.

Считая возмущения малыми, линеаризуем задачу (11) , (12) с учётом условия (13) . Получаем:

W + Li ^ = Г ^2 ^+L 2 ^+L 3 ^ ) , dt dx Re ^v                 '

P S Re V4 „ - ± ^^G^ dt S Bq     dy²

.

Уравнения (14) , (15) являются линейными с постоянными коэффициентами. Решение будем искать в виде:

{ Φ,Ψ,φ,ψ }∼ e i ( αx + βy - ωt ) ,

где α , β — волновые числа вдоль координат x и y , а комплексная величина ω определяет частоту и скорость роста возмущений. Подстановка (16) в (14) , (15) приводит к системе алгебраических уравнений с условием разрешимости:

( i(aL 1 — ^ )Re+(k ² — L ₂ ) β ² G 1 k ² Bq

- L 3

k ^{2 +} ^G - i^P s Re k ² Bq

Здесь k ² = a ² + в ² . Из уравнения видно, что симметрия между продольным (а) и поперечным (в) волновыми числами отсутствует, поскольку они содержатся в дисперсионном соотношении (17) в явном виде. Это является следствием анизотропной связи между объёмом основной жидкости и поверхностью. В классических задачах о конвективной неустойчивости, таких как, например, конвекция Рэлея–Бенара [45] , система изотропна в горизонтальной плоскости, и оба направления оказываются равноправными. Это приводит к вырождению нейтральной кривой по волновому числу: существует бесконечно много мод с одинаковым порогом устойчивости, отличающихся ориентацией волнового вектора. В обсуждаемой здесь задаче изотропия нарушается из-за выделения направления основного течения вдоль оси x .

Далее примем а = 1 и исследуем зависимость порога неустойчивости от волнового числа в в поперечном к основному течению направлении. Значение поверхностного числа Прандтля примем равным P S = 10 ^{- 4} , что соответствует олеиновой и стеариновой кислотам, использованным в качестве ПАВ в экспериментах из работы [42] . Решение системы (17) даёт набор критических значений (Re,Bq,в,^), при которых достигается порог неустойчивости. Прежде всего отметим согласованность применяемой модели с моделью Линде в отношении колебательного характера неустойчивости: при ^ = 0 уравнения (17) не имеют решения, что указывает на принципиально нестационарный, осциллирующий характер потери устойчивости. Ключевым отличием от результатов из работы [32] является то, что вихревая конвекция в слое ПАВ происходит при определённом наборе критических параметров (Рис. 2а ). Таким образом, используемая здесь модель предсказывает неустойчивость поверхности при конечном критическом числе Рейнольдса Re C , соответствующем характерной скорости U C объёмного течения (Рис. 1) . Поверхностная и объёмная вязкости подавляют возмущения и могут полностью остановить их развитие: с ростом числа Буссинеска равновесие стабилизируется (Рис. 2б ).

Область неустойчивости в плоскости (Re,Bq) построена для разных значений поперечного волнового числа в. Можно видеть, что порог Re C монотонно возрастает с увеличением Bq для всего диапазона исследованных

Рис. 2. Нейтральные кривые Re( β ) для набора чисел Буссинеска при α = 1 , P S = 10 ^-4 и в пределе Bq → 0 , отвечающего результатам работы Линде и Шварца [40] , предсказавшим беспороговую неустойчивость ( а ); область неустойчивости в плоскости параметров ( Re , Bq ) ( б )

значений β . Этот факт подтверждает стабилизирующую роль поверхностной вязкости: чем выше сопротивление слоя ПАВ сдвиговым деформациям, тем интенсивнее должно быть объёмное течение для того, чтобы привести поверхность в движение. Наиболее опасной является мода с β ≈ 0.5, для которой критическое число Рейнольдса минимально при достаточно малых значениях Bq. Длинноволновые возмущения (β = 0.1) требуют более высоких чисел Re C , что объясняется слабой связью этих мод с объёмным течением из-за малости градиента напряжения в поперечном направлении. Коротковолновые возмущения (β=5) становятся неустойчивыми при значительно бо´ льших Re, поскольку при меньших Re решающее значение приобретает поверхностная вязкая диссипация, вклад которой пропорционален величине k ² и которая эффективно подавляет эти возмущения. При любом фиксированном Bq зависимость Re C (β) имеет выраженный минимум, положение которого смещается с изменением числа Буссинеска. Таким образом, поверхностная вязкость выделяет наиболее опасную моду.

На рисунке 3 представлены зависимости критического числа Буссинеска Bq _C от поперечного волнового числа β для нескольких значений числа Рейнольдса. Квадратные символы выделяют максимальное значение Bq, при котором для заданной пары (Re, β) существует нейтральная мода с Im(ω) = 0. Другими словами, это наибольшая поверхностная вязкость, допускающая развитие неустойчивости при данных значениях параметров. Каждая кривая Bq _C (β) имеет выраженный максимум, положение которого отвечает наиболее опасной моде — волновому числу β C , для которого неустойчивость сохраняется при самой высокой поверхностной вязкости. Видно, что с ростом числа Рейнольдса максимум смещается: β C увеличивается от значений ∼ 0.5 при малых Re до величин порядка 1 при больших Re. Это означает, что с увеличением объёмного расхода характерный поперечный размер вихревых структур уменьшается: наиболее опасная мода становится более коротковолновой. Высота максимума Bq _C монотонно возрастает с увеличением числа Re: чем интенсивнее объёмное течение, тем более вязкий поверхностный слой оно способно привести в движение.

Рис. 3. Критическое число Буссинеска в зависимости от поперечного волнового числа β при разных значениях числа Рейнольдса; квадратными символами обозначены максимальные значения функции Bq( e ) , определяющие критическую моду β C , для которой развивается неустойчивость

Обнаруженный характер неустойчивости демонстрирует прямую аналогию с классической конвекцией Рэлея–Бенара, которую отмечают авторы в [42] . В задаче Рэлея–Бенара подогрев снизу создаёт вертикальный градиент температуры. При малых числах Рэлея тепло переносится молекулярной теплопроводностью. При достижении потоком тепловой энергии, проходящим через слой, критического значения наступает кризис теплопроводного режима и возникает конвективное движение. В обсуждаемой задаче сдвиговое течение в объёме создаётся потоком x-компоненты импульса к поверхности. При малых числах Рейнольдса импульс переносится вязкой диффузией, и слой ПАВ остаётся неподвижным. При достижении критического значения числа Re C перенос импульса от основного течения к поверхности резко изменяется (перестраивается), что приводит к развитию вихревой конвекции в слое ПАВ. В обоих случаях неустойчивость носит пороговый характер, а критическое значение управляющего параметра зависит от волнового числа возмущения и имеет выраженный минимум. Различие заключается в размерности развивающегося течения. В задаче Рэлея–Бенара конвекция трёхмерна и распространяется по всему объёму. В настоящей задаче — это отсутствие вырождения по волновому числу: нейтральная кривая имеет единственный минимум при β = β C . Данное обстоятельство является, по-видимому, характерной чертой сдвиговых неустойчивостей и отличает рассматриваемый механизм от изотропных конвективных неустойчивостей, индуцированных объёмными силами [44] .

• 4.    Нелинейные режимы движения

  • 4.1.    Метод численного решения

Система уравнений (11), (12) решалась численно методом конечных разностей. В эксперименте [42] использовалась цилиндрическая кювета с источником радиального движения, расположенным в центре. По периферии кюветы течение ограничивалось, и несжимаемая застойная зона представляла собой кольцевую область с фиксированными границами. Для упрощения вычислений расчётная область ограничивалась твёрдыми стенками по оси x, вдоль которой направлено основное течение в объёмной фазе. По оси y задавались периодические граничные условия. Рассматривалась прямоугольная область 0 ⩽ x ⩽ 2π, 0 ⩽ y ⩽ 8π. Пространственные дифференциальные операторы аппроксимировались центральными разностями на равномерной сетке.

Представленные ниже численные результаты получены на сетке с числом узлов 61 х 241. Данная сетка даёт

результат относительной погрешности значения, к которому сходится функция тока при дальнейшем сгущении, в пределах 5%. Нелинейные уравнения решались с использованием явной схемы. Для производной по времени применялась схема Эйлера первого порядка точности. Дискретный аналог функции f (x,y,t) обозначался как

m i,k

и давал значение функции в узле сетки i, k в момент времени t m

уравнений (11) , (12) с учётом принятых обозначений имели вид:

∆t n . Конечно-разностные аналоги

Φ _i ^m _k ⁺¹ =Φ i^mk +∆t m

Re ( V 2 ф т _к + L 2 $ m _k + L з ф m _k ) - NJ (Ф^ к - L 1

∂Φ i ^m k ∂x

m +1 m

φ _i,k = φ _i,k + ∆t _m

P S Re

1 ∂ ² ( G 1 Ψ i ^m k + G 2 ψ i ^m k )

Bq        ∂y ²

- j ^к ,^к )

Опишем алгоритм работы вычислительной схемы. Допустим, что в момент времени t _m известны все данные для полей Ψ _i ^m _k , ψ _i ^m _k , Φ _i ^m _k , φ _i ^m _k . На следующем временном шаге m +1 решаются уравнения Пуассона для функции тока и определяются новые значения Ψ _i ^m _k ⁺¹ , ψ _i ^m _k ⁺¹ . Затем, с использованием явного рекуррентного соотношения обновляются значения полей завихренности Φ _i ^m _k ⁺¹ , φ _i ^m _k ⁺¹ и процедура выходит на следующую итерацию. Уравнения Пуассона для функции тока решаются методом последовательной верхней релаксации. Точность вычислений задаётся значением 10 ^{- 5} . Для обеспечения устойчивости явной схемы на каждой итерации рассчитывается новый шаг по времени в соответствии с критерием Куранта:

∆t =

∆h ²

2(2+max( | ^ | , | Ф | )) ’

Начальным условием служит зашумлённое распределение функции тока с амплитудой 10 ^{- 4} . Завихренность на твёрдых границах вычисляется по формулам Кусковой–Чудова второго порядка точности [46] .

• 4.2.    Решение полной нелинейной задачи

В рамках численного исследования основное внимание уделено области параметров вблизи кризиса равновесия. С этой целью выбраны точки вблизи минимумов нейтральных кривых Re _C ( β ), представленных на рисунке 2а для различных значений числа Буссинеска, и прослежена эволюция структуры течения при изменении числа Рейнольдса в окрестности критического значения, отвечающего заданному Bq.

Прежде всего следует отметить, что в нелинейном расчёте порог неустойчивости оказывается несколько выше предсказываемого линейной теорией. По-видимому, это расхождение обусловлено учётом твёрдых стенок, ограничивающих расчётную область в продольном направлении. Это приводит к дискретизации спектра разрешённых волновых чисел и исключает моды, которые в неограниченной системе развивались бы при меньших значениях Re. В частности, для Bq= 0 . 01 вихревое движение начинается лишь при Re= 10, тогда как линейный анализ даёт Re c и 8. При небольших числах Re поток импульса через границу раздела остаётся недостаточно интенсивным для преодоления стабилизирующего действия поверхностной вязкости; слой ПАВ сохраняет неподвижность, несмотря на наличие объёмного течения. Увеличение интенсивности течения либо уменьшение поверхностной вязкости приводит к возникновению вихревого движения в застойной зоне. Возмущение скорости у поверхности вызывает локальное возрастание вязкой силы, нарушающее баланс между касательным напряжением со стороны объёма и градиентом поверхностного давления; в результате элемент двумерного слоя приходит в движение в направлении основного потока.

На рисунке 4 приведены поля функции тока на межфазной границе для различных чисел Рейнольдса и Буссинеска в фиксированный момент времени. При малых надкритичностях (Re = 20, Рис. 4а ) течение демонстрирует высокую степень упорядоченности: формируется периодическая система вихрей с практически постоянным поперечным волновым числом, близким к значению β C , найденному при линейном анализе. Данный вывод справедлив для всех трёх рассмотренных значений Bq с поправкой на стабилизирующее влияние поверхностной вязкости (см. числовые значения функции тока на изолиниях).

С ростом числа Рейнольдса картина качественно меняется. При Re = 30 (Рис. 4б ) наблюдается постепенное уменьшение характерного поперечного размера вихрей, то есть увеличение доминирующего волнового числа (числа, соответствующего максимуму спектральной плотности функции тока). Кроме того, становится заметным эффект адвективного сноса вихревой структуры вдоль оси x к правой границе расчётной области. Эффект обусловлен переносом возмущений основным течением и существует при любых надкритичностях, однако заметнее проявляется лишь при достаточном удалении от порога устойчивости. Влияние числа Буссинеска на структуру течения остаётся слабым и сказывается главным образом в виде небольшого изменения амплитуды вихрей. При дальнейшем увеличении числа Рейнольдса (Re = 40, Рис. 4в ) тенденция измельчения вихревой структуры сохраняется. Доминирующее волновое число приближается к своему предельному значению в ~ 1, полученному при линейном анализе. Адвективный снос к стенке становится существенным и меняет форму

( a )                                                  ( б )                                                  ( в )

о 50   50   50   50   50   50   50   50   5

X                            X                            X

Рис. 4. Поля функции тока на межфазной поверхности в фиксированный момент времени t = 2 при различных значениях числа Рейнольдса Re : 20 ( a ), 30 ( б ), 40 ( в ); обозначениям в рамках «1», «2», «3» соответствуют значения числа Буссинеска Bq : 0.01, 0.05, 0.1

вихревых образований. Стабилизирующее влияние числа Буссинеска проявляется в том, что при Bq= 0.1 (см. кадр 3) характерная длина волны оказывается несколько больше, чем при меньших значениях Bq. Этот факт согласуется с предсказаниями линейной теории о подавлении коротковолновых мод поверхностной вязкостью.

Рисунок 5 иллюстрирует эволюцию максимального по пространству значения функции тока для разных значений чисел Рейнольдса и Буссинеска, отвечающих кадрам на рисунке 4. Во всём исследованном диапазоне параметров система демонстрирует существенно нестационарное поведение: функция тока испытывает выраженные осцилляции, что качественно согласуется с предсказаниями ранних теоретических работ [40, 41] о колебательном характере потери устойчивости. Можно выделить две характерные составляющие колебаний. Первая — низкоамплитудные осцилляции высокой частоты, обусловленные, по-видимому, адвективным переносом вихревых структур основным течением вдоль оси x, их частота определяется скоростью объёмного потока и слабо зависит от надкритичности. Вторая составляющая — это высокоамплитудные низкочастотные осцилляции, обусловленные процессами перестройки вихревой структуры при слиянии и разрушении вихревых пар. При малых числах Рейнольдса (Рис. 5а ) данный тип колебаний практически отсутствует, то есть течение остаётся упорядоченным, а вихревая структура стабильна (Рис. 4а ). Частоты колебаний различаются приблизительно на порядок. С ростом Re низкочастотные осцилляции усиливаются, что свидетельствует о возрастающей роли нелинейных взаимодействий и постепенной хаотизации течения.

( а )                                                  ( б )                                                  ( в )

Рис. 5. Изменение во времени максимального значения поверхностной функции тока для разных чисел Буссинеска в зависимости от числа Рейнольдса Re : 20 ( а ), 30 ( б ), 40 ( в )

Число Буссинеска оказывает стабилизирующее влияние как на амплитуду установившегося движения, так и на амплитуду колебаний: с увеличением Bq обе величины заметно уменьшаются. Сопоставление кривых на рисунке 5а–в позволяет проследить за динамикой начального роста возмущений. С увеличением числа Рейнольдса функция тока на начальном этапе растёт существенно быстрее: выход на конечное значение происходит за время t ∼ 1 при Re = 20, заt∼0.5 приRe=30изаt∼0.25 приRe=40. Данная тенденция отражает увеличение линейного инкремента с удалением от порога устойчивости и хорошо согласуется с теоретическими предсказаниями. Частота высокоамплитудных колебаний также становится больше по мере роста Re (Рис. 5б, в), что указывает на ускорение процессов перестройки вихревой структуры и, в конечном счёте, на переход системы к хаотическому режиму при достаточно больших надкритичностях.

• 5.    Заключение

В настоящей работе представлена теоретическая модель сдвиговой неустойчивости слоя нерастворимого ПАВ, расположенного на границе раздела «жидкость–газ» в условиях стационарного объёмного течения. Управляющим параметром задачи выступает число Рейнольдса Re. Число Буссинеска Bq выражает отношение поверхностной вязкости к объёмной, служит характеристикой среды и является фактором стабилизации.

Показано, что потеря устойчивости основного состояния происходит пороговым образом и сопровождается возникновением в слое ПАВ периодической последовательности вихрей, ориентированных поперёк основного потока. Критическое значение Re C монотонно увеличивается с ростом Bq, а зависимость Re C (β) при конечном значении поперечного волнового числа имеет выраженный минимум, положение которого зависит от поверхностной вязкости. Полученные нейтральные кривые качественно согласуются с экспериментальными наблюдениями [32, 34, 42] . Система имеет порог устойчивости, существенно меньший предсказанного в теоретической работе [41] и отвечающий интервалу возможных значений Re C , который немного меньше диапазона от 20 до 1000, представленного в эксперименте в работе [42] для Bq ≈ 0.01.

Численное решение полной нелинейной системы уравнений подтверждает основные прогнозы линейной теории. Вблизи порога формируется упорядоченная вихревая структура с волновым числом, близким к β C . С ростом Re наблюдается уменьшение размера вихрей и увеличение значения доминирующего волнового числа. Эта тенденция сохраняется во всём исследованном диапазоне параметров. Система в нелинейном режиме характеризуется ярко выраженным нестационарным поведением: вихревая структура осциллирует, и чем больше её надкритичность, тем выше амплитуда и частота низкочастотных колебаний, обусловленных рождением и аннигиляцией вихрей. При достаточно больших Re система демонстрирует признаки перехода к хаотическому режиму. Обнаружен эффект адвективного сноса вихрей основным течением, проявляющийся заметнее при удалении от порога.

Работа выполнена в рамках крупного научного проекта при финансовой поддержке Минобрнауки России (соглашение № 075-15-2024-535 от 23.04.2024).