Кристаллография, 2020, T. 65, № 1, стр. 111-118

К вопросу о гидродинамической неустойчивости нематиков в магнитном поле. II. Линейный анализ устойчивости некоторых моделей магнитогидродинамических доменов

А. В. Голованов 1*

1 Институт элементоорганических соединений им. А.Н. Несмеянова РАН
Москва, Россия

* E-mail: gav@ineos.ac.ru

Поступила в редакцию 24.04.2018
После доработки 25.12.2018
Принята к публикации 28.01.2019

Полный текст (PDF)

Аннотация

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

ВВЕДЕНИЕ

К настоящему времени экспериментально и теоретически изучен процесс зарождения и развития магнитогидродинамических (МГД) доменов двух видов, возникающих в термотропных и лиотропных нематиках: полосатых (первого и второго рода) и паркетных. Полосатые домены первого рода представляют собой параллельные цилиндрические вихри, лежащие в плоскости слоя нематика [1], второго рода – вихревые линии с незамкнутыми линиями тока, расположенные перпендикулярно плоскости слоя нематика [24]. Паркетные домены – это цилиндрические вихри, расположенные под углом друг к другу в виде паркета [5, 6]. Появление доменной структуры происходит за порогом гидродинамической неустойчивости [7]. Нахождение порога неустойчивости (критичекого магнитного поля доменообразования) требует совместного решения уравнений движения директора и Навье–Стокса (НС), записанных в линейном приближении, но тип возникающей неустойчивости в этом случае не определяется.

Существует метод, известный как линейный анализ устойчивости (ЛАУ) динамических систем [8, 9]. Этот метод позволяет определять тип устойчивости/неустойчивости, реализующийся в системе при изменении управляющих параметров, и строить соответствующие фазовые траектории (интегральные кривые). Основой метода является оценка устойчивости стационарных решений линейных обыкновенных дифференциальных уравнений, содержащих производную первого порядка по времени, называемых эволюционными уравнениями (ЭУ). Для получения ЭУ в моделях МГД-доменов необходимо учитывать член ${\rho }\frac{{\partial {v}}}{{\partial t}}$ в линеаризованном уравнении НС [10]. В моделях [13, 5, 6] этим членом пренебрегают, за исключением модели для переориентации в случае неограниченной среды [1].

Цель настоящей работы – получение ЭУ и определение типа гидродинамической неустойчивости для моделей, предложенных в [1].

РЕЗУЛЬТАТЫ И ИХ ОБСУЖДЕНИЕ

Модель доменной структуры при переориентации нематика в случае неограниченной среды

В [1] учитывался член ${\rho }\frac{{\partial {v}}}{{\partial t}}$, но система уравнений в частных производных не была преобразована к ЭУ. В рамках модели получены зависимости обратного времени включения доменной структуры ${{s}_{ \pm }}$ от квадрата безразмерного волнового вектора $q_{x}^{2}$ и квадрата напряженности магнитного поля ${{H}^{2}}$, а также зависимость $q_{{xc}}^{2}({{H}^{2}})$.

В рассматриваемой модели неограниченный слой жидкого кристалла ориентируется магнитным полем вдоль оси x, затем поле быстро поворачивается в направлении оси z. Далее исследуется возможность появления вдоль оси x доменной структуры, сопровождаемой переходным потоком вещества. Скорость течения нематика задается одной компонентой ${{{v}}_{z}}$; компонентой ${{{v}}_{x}}$, обеспечивающей непрерывность поля скоростей на границах, пренебрегают. Кроме того, угол θ отклонения директора от положения равновесия и ${{{v}}_{z}}$ считаются не зависящими от z. Уравнения НС и движения директора в линейном приближении записываются в следующем виде:

(1)
$\left\{ \begin{gathered} {\rho }\frac{{\partial {{{v}}_{z}}}}{{\partial t}} = {{{\eta }}_{2}}\frac{{\partial _{{}}^{2}{{{v}}_{z}}}}{{\partial x_{{}}^{2}}} + {{\alpha }_{2}}\frac{\partial }{{\partial x}}\left( {\frac{{\partial {\theta }}}{{\partial t}}} \right), \hfill \\ {{{\gamma }}_{1}}\frac{{\partial {\theta }}}{{\partial t}} = {{{\chi }}_{a}}H_{{}}^{2}{\theta } + {{K}_{3}}\frac{{\partial _{{}}^{2}{\theta }}}{{\partial x_{{}}^{2}}} - {{\alpha }_{2}}\frac{{\partial {{{v}}_{z}}}}{{\partial x}}, \hfill \\ \end{gathered} \right.$
где ${{{\eta }}_{2}}$ – коэффициент вязкости, ${{\alpha }_{2}}$ – коэффициент вязкости Лесли, ${{{\gamma }}_{1}}$ – вращательная вязкость, ${{{\chi }}_{a}}$ – анизотропия диамагнитной восприимчивости, ${{K}_{3}}$ – константа упругости продольного изгиба, ρ – плотность нематика.

В отличие от [1], где решения системы (1) искали в виде

(2)
$\left\{ \begin{gathered} {{{v}}_{z}} = {{{v}}_{0}}\sin ({{q}_{x}}x)\exp (st), \hfill \\ {\theta } = {{{\theta }}_{0}}\cos ({{q}_{x}}x)\exp (st), \hfill \\ \end{gathered} \right.$
где ${{q}_{x}}$ – волновой вектор доменной структуры вдоль оси x, будем искать их как функции
(3)
${{{v}}_{z}} = {{{v}}_{0}}(t)\sin ({{q}_{x}}x),$
(4)
${\theta } = {{{\theta }}_{0}}(t)\cos ({{q}_{x}}x).$
Это позволяет после подстановки (3) и (4) в (1) получить ЭУ задачи
(5)
$\left\{ \begin{gathered} \frac{{d{{{\theta }}_{0}}}}{{dt}} = \frac{1}{{{{{\gamma }}_{1}}}}({{{\chi }}_{a}}H_{{}}^{2} - {{K}_{3}}q_{x}^{2}){{{\theta }}_{0}} - \frac{{{{\alpha }_{2}}}}{{{{{\gamma }}_{1}}}}{{q}_{x}}{{{v}}_{0}}, \hfill \\ \frac{{d{{{v}}_{0}}}}{{dt}} = \frac{{{{\alpha }_{2}}}}{{{\rho }{{{\gamma }}_{1}}}}({{K}_{3}}q_{x}^{2} - {{{\chi }}_{a}}H_{{}}^{2}){{q}_{x}}{{{\theta }}_{0}} + \hfill \\ + \;\frac{1}{{{\rho }{{{\gamma }}_{1}}}}(\alpha _{2}^{2} - {{{\eta }}_{2}}{{{\gamma }}_{1}})q_{x}^{2}{{{v}}_{0}}. \hfill \\ \end{gathered} \right.$
К (5) применяем ЛАУ. Существуют два стационарных решения (5). Первое:
(6)
${{{\theta }}_{{0{\text{ст}}1}}} \ne 0,\quad {{{v}}_{{0{\text{ст}}1}}} = 0$
при
(7)
$q_{x}^{2} = \frac{{{{{\chi }}_{a}}}}{{{{K}_{3}}}}H_{{}}^{2}.$
Выражение (7) является условием, при котором определитель системы линейных однородных уравнений, получаемых из (5), равен нулю. Второе:
(8)
${{{\theta }}_{{0{\text{ст}}2}}} = {{{v}}_{{0{\text{ст}}2}}} = 0$
при

(9)
$q_{x}^{2} \ne \frac{{{{{\chi }}_{a}}}}{{{{K}_{3}}}}H_{{}}^{2}.$

Для проверки на устойчивость решений (6) и (8) задаем их возмущения

(10)
$\left\{ \begin{gathered} {{{\theta }}_{{0{\text{в}}}}} = {{{\theta }}_{{0{\text{ст}}}}} + {\Theta }(t) \hfill \\ {{{v}}_{{0{\text{в}}}}} = {{{v}}_{{0{\text{ст}}}}} + V(t) \hfill \\ \end{gathered} \right.,$
где ${{{\theta }}_{{0{\text{в}}}}}$, ${{{v}}_{{0{\text{в}}}}}$ – возмущенные значения величин ${{{\theta }}_{{0{\text{ст}}}}}$, ${{{v}}_{{0{\text{ст}}}}}$, ${\Theta }(t)$ – возмущение угла ориентации директора, $V(t)$ – возмущение скорости течения нематика, и составляем уравнения для возмущений
(11)
$\left\{ \begin{gathered} \frac{{d{\Theta }}}{{dt}} = \frac{1}{{{{{\gamma }}_{1}}}}({{{\chi }}_{a}}H_{{}}^{2} - {{K}_{3}}q_{x}^{2}){{{\theta }}_{{0{\text{в}}}}} - \frac{{{{\alpha }_{2}}}}{{{{{\gamma }}_{1}}}}{{q}_{x}}{{{v}}_{{0{\text{в}}}}} = {{f}_{1}}, \hfill \\ \frac{{dV}}{{dt}} = \frac{{{{\alpha }_{2}}}}{{{\rho }{{{\gamma }}_{1}}}}({{K}_{3}}q_{x}^{2} - {{{\chi }}_{a}}H_{{}}^{2}){{q}_{x}}{{{\theta }}_{{0{\text{в}}}}} + \hfill \\ + \;\frac{1}{{{\rho }{{{\gamma }}_{1}}}}(\alpha _{2}^{2} - {{{\gamma }}_{1}}{{{\eta }}_{2}})q_{x}^{2}{{{v}}_{{0{\text{в}}}}} = {{f}_{2}}. \hfill \\ \end{gathered} \right.$
Далее, согласно процедуре ЛАУ, переходим к линеаризованным уравнениям для возмущений
(12)
$\left\{ \begin{gathered} \frac{{d{\Theta }}}{{dt}} = {{k}_{{11}}}{\Theta } + {{k}_{{12}}}V, \hfill \\ \frac{{dV}}{{dt}} = {{k}_{{21}}}{\Theta } + {{k}_{{22}}}V, \hfill \\ \end{gathered} \right.$
где
(13)
$\begin{gathered} {{k}_{{11}}} = {{\left( {\frac{{\partial {{f}_{1}}}}{{\partial {{{\theta }}_{{0{\text{в}}}}}}}} \right)}_{{{\text{ст}}}}} = \frac{1}{{{{{\gamma }}_{1}}}}({{{\chi }}_{a}}H_{{}}^{2} - {{K}_{3}}q_{x}^{2}), \\ {{k}_{{12}}} = {{\left( {\frac{{\partial {{f}_{1}}}}{{\partial {{{v}}_{{0{\text{в}}}}}}}} \right)}_{{{\text{ст}}}}} = - \frac{{{{\alpha }_{2}}}}{{{{{\gamma }}_{1}}}}{{q}_{x}}, \\ {{k}_{{21}}} = {{\left( {\frac{{\partial {{f}_{2}}}}{{\partial {{{\theta }}_{{0{\text{в}}}}}}}} \right)}_{{{\text{ст}}}}} = \frac{{{{\alpha }_{2}}}}{{{\rho }{{{\gamma }}_{1}}}}({{K}_{3}}q_{x}^{2} - {{{\chi }}_{a}}H_{{}}^{2}){{q}_{x}}, \\ {\text{ }}{{k}_{{22}}} = {{\left( {\frac{{\partial {{f}_{2}}}}{{\partial {{{v}}_{{0{\text{в}}}}}}}} \right)}_{{{\text{ст}}}}} = \frac{1}{{{\rho }{{{\gamma }}_{1}}}}(\alpha _{2}^{2} - {{{\gamma }}_{1}}{{{\eta }}_{2}})q_{x}^{2}. \\ \end{gathered} $
Находим характеристическое уравнение задачи
(14)
${{s}^{2}} - ({{k}_{{11}}} + {{k}_{{22}}})s + {{k}_{{11}}}{{k}_{{22}}} - {{k}_{{12}}}{{k}_{{21}}} = 0.$
Здесь s – обратное время включения доменной структуры.

Рассмотрим стационарное решение (6). Для нахождения корней (14) исключаем ${{q}_{x}}$ из коэффициентов (13) с помощью (7). Получаем

(15)
$\begin{gathered} {{k}_{{11}}} = 0,\quad {{k}_{{12}}} = - \frac{{{{\alpha }_{2}}}}{{{{{\gamma }}_{1}}}}\sqrt {\frac{{{{{\chi }}_{a}}}}{{{{K}_{3}}}}} H, \\ {{k}_{{21}}} = 0,\quad {{k}_{{22}}} = \left( {\frac{{\alpha _{2}^{2}}}{{{{{\gamma }}_{1}}}} - {{{\eta }}_{2}}} \right)\frac{{{{{\chi }}_{a}}H_{{}}^{2}}}{{{\rho }{{K}_{3}}}}. \\ \end{gathered} $
Отсюда следует, что уравнение (14) имеет два корня: ${{s}_{1}} = 0$ и ${{s}_{2}} = {{k}_{{22}}}$. Определим знак коэффициента ${{k}_{{22}}}$. Здесь и далее в качестве модельной системы используем термотропный нематик n-метоксибензилиден-n-бутиланилин (МВВА), поскольку для него определены все материальные константы, входящие в (15), а именно: ${{{\chi }}_{a}}$ = 0.97 × 10–7 СГС-ед, ρ = = 1.088 г/см3, ${{K}_{3}}$ = 7.5 × 10–7 дин, ${{\alpha }_{2}}$ = –0.775 П, ${{{\gamma }}_{1}}$ = = 0.77 П и ${{{\eta }}_{2}}$ = 1.03 П [11]. Знак коэффициента ${{k}_{{22}}}$ определяет разность, стоящая в скобках. Так как $\frac{{\alpha _{2}^{2}}}{{{{{\gamma }}_{1}}}} < {{{\eta }}_{2}}$, то ${{k}_{{22}}} < 0$ при любых значениях H, отличных от нуля. Следовательно, ${{s}_{2}} < 0$ и закон изменения возмущений в зависимости от времени имеет вид
(16)
$\begin{gathered} {\Theta }(t) = {{c}_{{11}}} + {{c}_{{12}}}\exp ({{s}_{2}}t), \hfill \\ V(t) = {{c}_{{21}}} + {{c}_{{22}}}\exp ({{s}_{2}}t). \hfill \\ \end{gathered} $
Из (16) видно, что с течением времени возмущения угла ориентации директора ${\Theta }(t)$ и скорости V(t) убывают, а при $t \to \infty $ стремятся к начальным возмущениям ${{c}_{{11}}}$ и ${{c}_{{21}}}$. Поэтому стационарное решение (6) является устойчивым. В этом случае тип устойчивости определяется как нейтральная устойчивость [8, 9].

Как известно, исследование устойчивости нелинейных дифференциальных уравнений проводится также методом определения фазовых траекторий (интегральных кривых) в фазовой плоскости [8]. Исключая время в (5) с последующим интегрированием, получаем интегральную кривую вида ${{{v}}_{0}} = \frac{{{{k}_{{22}}}}}{{{{k}_{{12}}}}}{{{\theta }}_{0}} + {\text{const}}$. Она показана на рис. 1а. Стрелками обозначена совокупность интегральных кривых, отличающихся начальными условиями. Их направление указывает на движение системы к состоянию устойчивого равновесия.

Рис. 1.

Интегральные кривые системы: а – для устойчивого стационарного состояния ${{{\theta }}_{{0{\text{ст}}1}}} \ne 0$, ${{{v}}_{{0{\text{ст}}1}}} = 0$: H = 100 Э, t = 10 с; б – для неустойчивого стационарного состояния ${{{\theta }}_{{0{\text{ст}}2}}} = 0$, ${{{v}}_{{0{\text{ст}}2}}} = 0$: H = 100 кЭ, t = 10–6 с. Стрелками обозначены совокупности интегральных кривых, различающихся начальными условиями.

Рассмотрим стационарное решение (8). Для анализа его устойчивости также необходимо исключить ${{q}_{x}}$ из (13). Это можно сделать, исследовав на экстремум функцию, являющуюся решением (14):

(17)
$\begin{gathered} {{s}_{{1,2}}} = \frac{1}{2}\{ A{{H}^{2}} + Bq_{x}^{2} \pm [{{A}^{2}}{{H}^{4}} + \\ + \;2C{{H}^{2}}q_{x}^{2} + Dq_{x}^{4}{{]}^{{1/2}}} \} , \\ \end{gathered} $
где
$A = \frac{{{{{\chi }}_{a}}}}{{{{{\gamma }}_{1}}}},\quad B = \frac{1}{{\rho }}\left( {\frac{{\alpha _{2}^{2}}}{{{{{\gamma }}_{1}}}} - {{{\eta }}_{2}}} \right) - \frac{{{{K}_{3}}}}{{{{{\gamma }}_{1}}}},$
$C = \frac{{{{{\chi }}_{a}}}}{{{{{\gamma }}_{1}}}}\left( {\frac{{\alpha _{2}^{2}}}{{{\rho }{{{\gamma }}_{1}}}} + \frac{{{{{\eta }}_{2}}}}{{\rho }} - \frac{{{{K}_{3}}}}{{{{{\gamma }}_{1}}}}} \right),$
$D = {{\left( {\frac{{\alpha _{2}^{2}}}{{{\rho }{{{\gamma }}_{1}}}} - \frac{{{{{\eta }}_{2}}}}{{\rho }} - \frac{{{{K}_{3}}}}{{{{{\gamma }}_{1}}}}} \right)}^{2}} - \frac{{4{{{\eta }}_{2}}{{K}_{3}}}}{{{\rho }{{{\gamma }}_{1}}}}.$
В результате имеем
(18)
$q_{x}^{2} = a{{H}^{2}},\quad {\text{где}}\quad a = \frac{{C + {{{\left[ {\frac{{{{B}^{2}}({{C}^{2}} - {{A}^{2}}D)}}{{{{B}^{2}} - D}}} \right]}}^{{1/2}}}}}{D}.$
Формулы (7), (17) и (18) совпадают с выражениями для $q_{{xF}}^{2}$, ${{s}_{ \pm }}$ и $q_{{xc}}^{2}$, полученными в [1]. Значения волновых векторов $q_{{xc}}^{2}$, как показано в [1], являются точками максимумов функции ${{s}_{ + }} = f(q_{x}^{2})$ для любых H ≥ 0. Волновые векторы ${{q}_{{xF}}}$ таковы, что ${{s}_{ + }}({{q}_{{xF}}}) = 0$, и определяются через “условие Фредерикса” $H_{{_{C}}}^{2} = q_{{xF}}^{2}\frac{{{{K}_{3}}}}{{{{{\chi }}_{a}}}}$, которое совпадает с (7). Отметим, что из (17) при условии ${{s}_{{1,2}}}$ = 0 следует (7).

На рис. 2 в системе координат ${{q}_{x}}$, $H$, $s$ показаны функции (17), (18) и (7). Поверхность $s = f({{q}_{x}},\;H)$ составлена из кривых $s = f({{q}_{x}})$, параметризованных переменной H. Эти кривые имеют максимумы, через которые проходит прямая ${{q}_{x}} = \sqrt a H$. На рисунке через эти точки проходит плоскость $\sqrt a H - {{q}_{x}} = 0$. В плоскости ${{q}_{x}}$H расположена прямая ${{q}_{x}} = \sqrt {\frac{{{{{\chi }}_{a}}}}{{{{K}_{3}}}}} H$.

Рис. 2.

Функции (7), (17) и (18), представленные в системе координат ${{q}_{x}}$, $H$, $s$.

Подставляя (18) в (17), имеем

${{s}_{{1,2}}} = \frac{1}{2}\{ (A + aB \pm {{[{{A}^{2}} + {{a}^{2}}D + 2aC]}^{{1/2}}}){{H}^{2}} \} ,$
где $A\, = \,1.26 \times 10_{{}}^{{ - 7}}\;\frac{{{\text{см}}\;{\text{с}}}}{{\text{г}}}$, $B\, = \, - {\kern 1pt} 0.25\;\frac{{{\text{с}}{{{\text{м}}}^{{\text{2}}}}}}{{\text{с}}}$, $C\, = \,2.095 \times $ $ \times \;{{10}^{{ - 7}}}\;\frac{{{\text{с}}{{{\text{м}}}^{{\text{3}}}}}}{{\text{г}}}$, $D = 0.053\;\frac{{{\text{с}}{{{\text{м}}}^{{\text{4}}}}}}{{{{{\text{с}}}^{{\text{2}}}}}}$, $a = 4.74 \times {{10}^{{ - 4}}}\;\frac{1}{{{\text{с}}{{{\text{м}}}^{{\text{2}}}}\;{{{\text{Э}}}^{{\text{2}}}}}}$. Вычисления показывают, что ${{s}_{1}} > 0$ и ${{s}_{2}} < 0$ при любом $H > 0$. Решение (12) имеет вид
(19)
$\begin{gathered} {\Theta }(t) = {{c}_{{11}}}\exp ({{s}_{1}}t) + {{c}_{{12}}}\exp ({{s}_{2}}t), \hfill \\ V(t) = {{c}_{{21}}}\exp ({{s}_{1}}t) + {{c}_{{22}}}\exp ({{s}_{2}}t). \hfill \\ \end{gathered} $
Очевидно, что при $t \to \infty $ возмущения увеличиваются с течением времени. Данное стационарное решение является неустойчивым. Это седловая неустойчивость [8, 9]. На рис. 1б показана соответствующая интегральная кривая. Так же, как и на рис. 1а, стрелками обозначена совокупность интегральных кривых, различающихся начальными условиями. Как видно из рисунка, движение вдоль интегральной кривой уводит систему от стационарного состояния.

Анализ показывает, что при любом H > 0 существуют два набора волновых векторов (формулы (7) и (18)), которые соответствуют разным стационарным решениям (5). Из экспериментов [1, 2, 46] следует, что МГД-домены в отличие от электрогидродинамических доменов [10] являются переходными – возникают, развиваются и исчезают, когда поворот директора от положения равновесия (при заданном поле) достигает максимума. Поэтому устойчивое состояние (6) в случае появления доменов в неограниченной среде может возникнуть только после их исчезновения – директор повернулся на угол θ, и течение нематика прекратилось. Следовательно, волновые векторы, определяемые равенством (7), не реализуются, в том числе в силу того, что оно следует из условия s = 0.

Таким образом, реализуются только волновые векторы (18), связанные с неустойчивым состоянием ${{{\theta }}_{{0{\text{ст}}2}}} = {{{v}}_{{0{\text{ст}}2}}} = 0$, поскольку их значения в точности соответствуют точкам максимумов функции (17) (рис. 2).

Кроме того, отметим, что домены возникают при сколь угодно малых значениях магнитного поля, как следует из (18). Это говорит о том, что критического поля доменообразования в этой модели не существует.

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

В данной модели рассматривается планарный слой нематика с положительной анизотропией диамагнитной восприимчивости и бесконечно сильным сцеплением директора на границах. Вектор ${{n}_{0}}$ невозмущенной ориентации нематика направлен вдоль оси x декартовой системы координат, ось z перпендикулярна плоскости слоя, начало координат выбрано в центре слоя толщиной d (рис. 3). Магнитное поле H прикладывается перпендикулярно невозмущенной ориентации директора ${{n}_{0}}$ (вдоль оси z).

Рис. 3.

На верхнем рисунке представлена геометрия задачи, на нижнем – схематичное изображение вихрей в плоскости xz. Оси вихрей параллельны оси y, d – высота вихрей по оси z. Вся структура периодична по оси x. Верхняя и нижняя поверхности нематика определяются уравнением $z = \pm \frac{d}{2}$.

При быстром включении магнитного поля величиной больше ${{H}_{D}}$ (критическое поле доменообразования) через некоторое время в слое жидкого кристалла появляется периодическая структура, состоящая из темных и светлых полос – доменов. Для описания доменов рассматриваются малые возмущения директора и вектора скорости течения нематика

$n \equiv ({\theta }(x,z),0,1),\quad {\mathbf{v}} \equiv ({{{v}}_{x}}(x,z),0,{{{v}}_{z}}(x,z)).$
Для заданной геометрии уравнения движения директора и НС в линейном приближении имеют вид:
(20)
$\left\{ \begin{gathered} {{{\gamma }}_{1}}\frac{{\partial {\theta }}}{{\partial t}} = {{{\chi }}_{a}}{{H}^{2}}{\theta } + {{K}_{1}}\frac{{{{\partial }^{2}}{\theta }}}{{\partial {{z}^{2}}}} + {{K}_{3}}\frac{{{{\partial }^{2}}{\theta }}}{{\partial {{x}^{2}}}} - \hfill \\ - \;{{\alpha }_{2}}\frac{{\partial {{{v}}_{z}}}}{{\partial x}} - {{\alpha }_{3}}\frac{{\partial {{{v}}_{x}}}}{{\partial z}}, \hfill \\ {\rho }\frac{{\partial {{{v}}_{x}}}}{{\partial t}} = - \frac{{\partial p}}{{\partial x}} + {{{\eta }}_{1}}\frac{{{{\partial }^{2}}{{{v}}_{x}}}}{{\partial {{z}^{2}}}} + {{{\eta }}_{4}}\frac{{{{\partial }^{2}}{{{v}}_{z}}}}{{\partial x\partial z}} + \hfill \\ + \;{{{\eta }}_{5}}\frac{{{{\partial }^{2}}{{{v}}_{x}}}}{{\partial {{x}^{2}}}} + {{\alpha }_{3}}\frac{\partial }{{\partial z}}\left( {\frac{{\partial {\theta }}}{{\partial t}}} \right), \hfill \\ {\rho }\frac{{\partial {{{v}}_{z}}}}{{\partial t}} = - \frac{{\partial p}}{{\partial z}} + {{{\eta }}_{2}}\frac{{{{\partial }^{2}}{{{v}}_{z}}}}{{\partial {{x}^{2}}}} + {{{\eta }}_{6}}\frac{{{{\partial }^{2}}{{{v}}_{x}}}}{{\partial x\partial z}} + \hfill \\ + \;{{\alpha }_{4}}\frac{{{{\partial }^{2}}{{{v}}_{z}}}}{{\partial {{z}^{2}}}} + {{\alpha }_{2}}\frac{\partial }{{\partial x}}\left( {\frac{{\partial {\theta }}}{{\partial t}}} \right), \hfill \\ \end{gathered} \right.$
где ${{{\eta }}_{j}}$ – коэффициенты вязкости, выраженные через различные комбинации коэффициентов вязкости Лесли, ${{\alpha }_{i}}$ – коэффициенты вязкости Лесли, ${{{\chi }}_{a}}$ – анизотропия диамагнитной восприимчивости, ${{{\gamma }}_{1}}$ – вращательная вязкость, ρ – плотность нематика, p – давление, ${{K}_{1}}$ – константа упругости поперечного изгиба, ${{K}_{3}}$ – константа упругости продольного изгиба.

Для уменьшения числа переменных в (20) введем в рассмотрение функцию тока Лагранжа ${\psi }(x,\;z,\;t)$:

${{{v}}_{x}} = - \frac{{\partial {\psi }}}{{\partial z}},\quad {{{v}}_{z}} = \frac{{\partial {\psi }}}{{\partial x}}.$
Введенная таким образом функция тока удовлетворяет уравнению непрерывности
$\frac{{\partial {{{v}}_{x}}}}{{\partial x}} + \frac{{\partial {{{v}}_{z}}}}{{\partial z}} = 0.$
Запишем (20) в безразмерном виде
(21)
$\left\{ \begin{gathered} \frac{{\partial {\theta }}}{{\partial{ \tilde {t}}}} = {{h}^{2}}{\theta } + \frac{{{{\partial }^{2}}{\theta }}}{{\partial {{{\tilde {z}}}^{2}}}} + K\frac{{{{\partial }^{2}}{\theta }}}{{\partial {{{\tilde {x}}}^{2}}}} + \frac{{{{\alpha }_{3}}}}{{{{\alpha }_{2}}}}\frac{{{{\partial }^{2}}{\tilde {\psi }}}}{{\partial {{{\tilde {z}}}^{2}}}} - \frac{{{{\partial }^{2}}{\tilde {\psi }}}}{{\partial {{{\tilde {x}}}^{2}}}}, \hfill \\ \frac{\partial }{{\partial{ \tilde {t}}}}\left( {\frac{{{{\partial }^{2}}{\tilde {\psi }}}}{{\partial {{{\tilde {x}}}^{2}}}}} \right) + \frac{\partial }{{\partial{ \tilde {t}}}}\left( {\frac{{{{\partial }^{2}}{\tilde {\psi }}}}{{\partial {{{\tilde {z}}}^{2}}}}} \right) = \frac{{{{{\gamma }}_{{\text{1}}}}{{{\eta }}_{1}}}}{{{\rho }{{K}_{1}}}}\frac{{{{\partial }^{4}}{\tilde {\psi }}}}{{\partial {{{\tilde {z}}}^{4}}}} + \hfill \\ + \;\frac{{{{{\gamma }}_{{\text{1}}}}{{{\eta }}_{7}}}}{{{\rho }{{K}_{1}}}}\frac{{{{\partial }^{4}}{\tilde {\psi }}}}{{\partial {{{\tilde {x}}}^{2}}\partial {{{\tilde {z}}}^{2}}}} - \frac{{{{{\gamma }}_{{\text{1}}}}{{{\eta }}_{{\text{2}}}}}}{{{\rho }{{K}_{1}}}}\frac{{{{\partial }^{4}}{\tilde {\psi }}}}{{\partial {{{\tilde {x}}}^{4}}}} + \hfill \\ + \;\frac{{\alpha _{2}^{2}}}{{{\rho }{{K}_{1}}}}\frac{{{{\partial }^{2}}}}{{\partial {{{\tilde {x}}}^{2}}}}\left( {\frac{{\partial {\theta }}}{{\partial{ \tilde {t}}}}} \right) - \frac{{{{\alpha }_{{\text{2}}}}{{\alpha }_{3}}}}{{{\rho }{{K}_{1}}}}\frac{{{{\partial }^{2}}}}{{\partial {{{\tilde {z}}}^{2}}}}\left( {\frac{{\partial {\theta }}}{{\partial{ \tilde {t}}}}} \right), \hfill \\ \end{gathered} \right.$
где $\tilde {x} = \frac{{\pi }}{d}x$ и $\tilde {z} = \frac{{\pi }}{d}z$ – безразмерные координаты, ${\tilde {\psi }} = \frac{{{{{\alpha }}_{2}}}}{{{{K}_{1}}}}{\psi }$ – безразмерная функция тока, $\tilde {t} = \frac{t}{{{{{\tau }}_{0}}}}$ – безразмерное время, $h = \frac{H}{{{{H}_{F}}}}$, $K = \frac{{{{K}_{3}}}}{{{{K}_{1}}}}$ – безразмерные параметры задачи, ${{{\tau }}_{0}} = \frac{{{{{\gamma }}_{1}}{{d}^{2}}}}{{{{{\pi }}^{2}}{{K}_{1}}}}$ – характерное время задачи, ${{H}_{F}} = \frac{{\pi }}{d}\sqrt {\frac{{{{K}_{1}}}}{{{{{\chi }}_{a}}}}} $ – критическое поле Фредерикса для деформации поперечного изгиба, ${{{\eta }}_{7}} = {{\alpha }_{1}} - {{\alpha }_{2}} + {{\alpha }_{4}} + {{\alpha }_{6}}$.

Граничные условия задачи имеют вид

(22)
${{\left. {\frac{{\partial {\tilde {\psi }}}}{{\partial{ \tilde {x}}}}} \right|}_{{\tilde {z} = \pm \tfrac{{\pi }}{2}}}} = 0;\quad {{\left. {\frac{{{{\partial }^{2}}{\tilde {\psi }}}}{{\partial {{{\tilde {z}}}^{2}}}}} \right|}_{{\tilde {z} = \pm \tfrac{{\pi }}{{\text{2}}}}}} = 0;\quad {{\left. {\theta } \right|}_{{\tilde {z} = \pm \tfrac{{\pi }}{{\text{2}}}}}} = 0.$

Решения уравнений (21) будем искать в виде функций, периодических вдоль оси x:

(23)
$\left\{ \begin{gathered} {\tilde {\psi }} = {{{{\tilde {\psi }}}}_{0}}(\tilde {t})\cos ({{{\tilde {q}}}_{x}}\tilde {x})\cos {\text{(}}{{{\tilde {q}}}_{z}}\tilde {z}), \hfill \\ {\theta } = {{{\theta }}_{0}}(\tilde {t})\cos ({{{\tilde {q}}}_{x}}\tilde {x})\cos {\text{(}}{{{\tilde {q}}}_{z}}\tilde {z}), \hfill \\ \end{gathered} \right.$
где ${{\tilde {q}}_{x}}$ и ${{\tilde {q}}_{z}}$ – безразмерные волновые векторы вдоль осей x и z. В [5] показано, что ${{\tilde {q}}_{z}} = 1$ для случая бесконечно сильного сцепления директора на границах. Поэтому
$\left\{ \begin{gathered} {\tilde {\psi }} = {{{{\tilde {\psi }}}}_{0}}(\tilde {t})\cos ({{{\tilde {q}}}_{x}}\tilde {x})\cos \tilde {z}, \hfill \\ {\theta } = {{{\theta }}_{0}}(\tilde {t})\cos ({{{\tilde {q}}}_{x}}\tilde {x})\cos \tilde {z}, \hfill \\ \end{gathered} \right.$
что удовлетворяет граничным условиям (22).

Далее повторяем процедуру, описанную в предыдущей части, и получаем ЭУ

(24)
$\left\{ \begin{gathered} \frac{{d{{{\theta }}_{0}}}}{{d\tilde {t}}} = (h_{{}}^{2} - K\tilde {q}_{x}^{2} - 1){{{\theta }}_{0}} + (\tilde {q}_{x}^{2} - a){{{{\tilde {\psi }}}}_{0}}, \hfill \\ \frac{{d{{{{\tilde {\psi }}}}_{0}}}}{{d\tilde {t}}} = \frac{{{\nu }(h_{{}}^{2} - K\tilde {q}_{x}^{2} - 1)(\tilde {q}_{x}^{2} - a)}}{{\tilde {q}_{x}^{2} + 1}}{{{\theta }}_{0}} + \hfill \\ + \;\frac{{[{\nu }{{{(\tilde {q}_{x}^{2} - a)}}^{2}} - {{{\nu }}_{1}}(b\tilde {q}_{x}^{4} + c\tilde {q}_{x}^{2} + 1)]}}{{\tilde {q}_{x}^{2} + 1}}{{{{\tilde {\psi }}}}_{0}}, \hfill \\ \end{gathered} \right.$
где ${\nu } = \frac{{{\alpha }_{{\text{2}}}^{{\text{2}}}}}{{{\rho }{{K}_{1}}}}$ и ${{{\nu }}_{1}} = \frac{{{{{\gamma }}_{{\text{1}}}}{{{\eta }}_{1}}}}{{{\rho }{{K}_{1}}}}$ – отношения чисел Эриксена E и Рейнольдса R [12], $a = \frac{{{{\alpha }_{3}}}}{{{{\alpha }_{2}}}}$, $b = \frac{{{{{\eta }}_{2}}}}{{{{{\eta }}_{1}}}}$, $c = \frac{{{{{\eta }}_{7}}}}{{{{{\eta }}_{1}}}}$, и стационарные решения:
(25)
${{{\theta }}_{{{\text{0 ст}}1}}} \ne 0,\quad {{{\tilde {\psi }}}_{{{\text{0 ст}}1}}} = 0$
при
(26)
$\tilde {q}_{x}^{2} = \frac{{{{h}^{2}} - 1}}{K}$
и
(27)
${{{\theta }}_{{{\text{0 ст}}2}}} = {{{\tilde {\psi }}}_{{{\text{0 ст}}2}}} = 0$
при
(28)
$\tilde {q}_{x}^{2} \ne \frac{{{{h}^{2}} - 1}}{K}.$
Линеаризованные уравнения для возмущений
(29)
$\left\{ \begin{gathered} \frac{{d{\Theta }}}{{d\tilde {t}}} = {{k}_{{11}}}{\Theta } + {{k}_{{12}}}{\tilde {\Psi },} \hfill \\ \frac{{d{\tilde {\Psi }}}}{{d\tilde {t}}} = {{k}_{{21}}}{\Theta } + {{k}_{{22}}}{\tilde {\Psi },} \hfill \\ \end{gathered} \right.$
где
(30)
$\begin{gathered} {{k}_{{11}}} = (h_{{}}^{2} - K\tilde {q}_{x}^{2} - 1),\quad {{k}_{{12}}} = (\tilde {q}_{x}^{2} - a), \\ {{k}_{{21}}} = \frac{{{\nu }(h_{{}}^{2} - K\tilde {q}_{x}^{2} - 1)(\tilde {q}_{x}^{2} - a)}}{{\tilde {q}_{x}^{2} + 1}}, \\ {{k}_{{22}}} = \frac{{{\nu }{{{(\tilde {q}_{x}^{2} - a)}}^{2}} - {{{\nu }}_{1}}(b\tilde {q}_{x}^{4} + c\tilde {q}_{x}^{2} + 1)}}{{\tilde {q}_{x}^{2} + 1}}. \\ \end{gathered} $
Рассмотрим стационарное решение (25). Исключая $\tilde {q}_{x}^{2}$ из (30) с помощью (26), получим
(31)
${{k}_{{11}}} = 0,\quad {{k}_{{12}}} = \frac{{(h_{{}}^{2} - 1)}}{K} - a,\quad {{k}_{{21}}} = 0,\quad {{k}_{{22}}} = \frac{{({\nu } - b{{{\nu }}_{1}}){{{(h_{{}}^{2} - 1)}}^{2}}\frac{1}{K} - (2a{\nu } + c{{{\nu }}_{1}})(h_{{}}^{2} - 1) + K({\nu }a_{{}}^{2} - {{{\nu }}_{1}})}}{{h_{{}}^{2} - 1 + K}}.$
Устойчивость/неустойчивость стационарного решения (25), также как и (6), зависит от знака коэффициента ${{k}_{{22}}}$. Поскольку ${{K}_{1}}$ = 6 × 10–7 дин, α3 = = –0.01 П, η1 = 0.24 П, η7 = 1.35 П, то K = 1.25, a = = 0.013, b = 4.29, c = 5.63, ν = 9.2 × 105, ν1 = 1.21 × × 106. Тогда $b{{{\nu }}_{1}} > {\nu }$ и ${{{\nu }}_{1}} > {\nu }a_{{}}^{2}$, а ${{k}_{{22}}} < 0$ при любых значениях ${{h}^{2}} \geqslant 1$. Следовательно, все выводы, сформулированные выше для стационарного решения (6), остаются справедливыми и для (25). Стационарное состояние системы является устойчивым, и тип устойчивости – нейтральная устойчивость. На рис. 4а показана интегральная кривая ${{{\tilde {\psi }}}_{0}} = - \frac{{{{k}_{{22}}}}}{{{{k}_{{12}}}}}{{{\theta }}_{0}} + {\text{const}}$. Как и в первой модели, система, двигаясь вдоль интегральной кривой, стремится к устойчивому состоянию.

Рис. 4.

Интегральные кривые системы: а – для устойчивого стационарного состояния ${{{\theta }}_{{{\text{0 ст}}1}}} \ne 0$, ${{{\tilde {\psi }}}_{{{\text{0 ст}}1}}} = 0$: h2 = 16, t = 10–4 с; б – для неустойчивого стационарного состояния ${{{\theta }}_{{{\text{0 ст}}2}}} = {{{\tilde {\psi }}}_{{{\text{0 ст}}2}}} = 0$: h2 = 14.44, t = 10–5 с. Стрелками обозначены совокупности интегральных кривых, различающихся начальными условиями.

Рассмотрим стационарное решение (27). Решением уравнения (14) для $\tilde {s}$ – обратного безразмерного времени включения доменной структуры – является функция $\tilde {s}({{h}^{2}},\tilde {q}_{x}^{2})$, которая имеет следующий вид:

(32)
${{\tilde {s}}_{{1,2}}} = \frac{1}{2}\left\{ {\frac{{h_{{}}^{2}(\tilde {q}_{x}^{2} + 1) + A\tilde {q}_{x}^{4} - B\tilde {q}_{x}^{2} + {\nu }a_{{}}^{2} - {{{\nu }}_{1}} - 1}}{{\tilde {q}_{x}^{2} + 1}} \pm {{{\left[ \begin{gathered} {{\left( {\frac{{h_{{}}^{2}(\tilde {q}_{x}^{2} + 1) + A\tilde {q}_{x}^{4} - B\tilde {q}_{x}^{2} + {\nu }a_{{}}^{2} - {{{\nu }}_{1}} - 1}}{{\tilde {q}_{x}^{2} + 1}}} \right)}^{2}} + \hfill \\ + \;\frac{{4(h_{{}}^{2} - K\tilde {q}_{x}^{2} - 1)(C\tilde {q}_{x}^{4} + D\tilde {q}_{x}^{2} + {{{\nu }}_{1}})}}{{\tilde {q}_{x}^{2} + 1}} \hfill \\ \end{gathered} \right]}}^{{1/2}}} } \right\},$
где $A = {\nu } - b{{{\nu }}_{1}} - K \approx {\nu } - b{{{\nu }}_{1}}$, $B = 2a{\nu } + c{{{\nu }}_{1}} + K + $ $ + \;1 \approx 2a{\nu } + c{{{\nu }}_{1}}$, $C = b{{{\nu }}_{1}}$, $D = c{{{\nu }}_{1}}$. Из (32) видно, что условие (26) реализуется при $\tilde {s} = 0$.

Зависимость $\;\tilde {q}_{x}^{2}({{h}^{2}})$ для MBBA найдена экспериментально в [1]. Показано, что экспериментальные точки $(h_{i}^{2},\tilde {q}_{{xi}}^{2})$ группируются около кривой, проводимой через точки максимумов семейства кривых $\tilde {s}(\tilde {q}_{x}^{2})$, параметризованных переменной $\;{{h}^{2}}$. В ходе эксперимента изменение значений $\;{{h}^{2}}$ в пределах от 17 до 400 приводило к изменению значений $\;\tilde {q}_{x}^{2}$ в пределах от 4 до 32. Эти данные позволяют определить знаки ${{\tilde {s}}_{{1.2}}}$, не исследуя (32) на экстремум для получения зависимости $\;\tilde {q}_{x}^{2}({{h}^{2}})$. Так как соотношение $\frac{{{{h}^{2}}}}{{\tilde {q}_{x}^{2}}}\;$ лежит в диапазоне значений ∼4–10, то выражение $h_{{}}^{2} - K\tilde {q}_{x}^{2} - 1$ для измеренных значений ${{h}^{2}}$ и$\;\tilde {q}_{x}^{2}$ не отрицательно, поэтому второе слагаемое в (32) больше первого. Следовательно, ${{\tilde {s}}_{1}} > 0$ и ${{\tilde {s}}_{2}} < 0$.

Поскольку решением (29) являются функции

$\begin{gathered} {\Theta }(\tilde {t}) = {{c}_{{11}}}\exp ({{s}_{1}}\tilde {t}) + {{c}_{{12}}}\exp ({{s}_{2}}\tilde {t}), \\ {\tilde {\Psi }}(\tilde {t}) = {{c}_{{21}}}\exp ({{s}_{1}}\tilde {t}) + {{c}_{{22}}}\exp ({{s}_{2}}\tilde {t}), \\ \end{gathered} $
то, как и в предыдущей модели, стационарное решение (27) является неустойчивым – это седловая неустойчивость. Интегральная кривая на рис. 4б иллюстрирует вывод об уходе системы от положения равновесия при движении вдоль интегральной кривой.

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

Стационарное состояние (25) для ограниченной среды может возникнуть в двух случаях. Во-первых, если на границах задана жесткая наклонная ориентация директора, которая при отсутствии внешнего поля или течения задает наклонную ориентацию всего слоя нематика [10]. Во-вторых, оно возникает после исчезновения доменов – течение прекратилось, а директор повернулся на угол θ. Поскольку в модели не задана наклонная ориентация директора на границах, то реализуется второй случай.

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

ЗАКЛЮЧЕНИЕ

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

Учет члена ${\rho }\frac{{\partial {v}}}{{\partial t}}$ в уравнении НС позволяет помимо возможности применять ЛАУ ввести в модель доменной структуры числа Рейнольдса R = = $\frac{{{\rho }{v}l}}{{\eta }}$ и Эриксона E = $\frac{{{\eta }{v}l}}{K}$ – безразмерные соотношения параметров среды (вязкостей η, констант упругости K, плотности ρ нематика), а также характерных скорости ${v}$ и расстояния l модели. Соотношение чисел Рейнольдса и Эриксена ${\mu } = \frac{{{\rho }K}}{{{{{\eta }}^{2}}}}$ (или обратное соотношение ${\nu } = \frac{{{{{\eta }}^{2}}}}{{{\rho }K}}$) исключает из модели характерные величины и оставляет лишь параметры среды (24). Получение значений этих чисел и их соотношений в экспериментах по наблюдению за МГД-доменами позволяет оценить вклад сил различной природы в уравнениях НС и движения директора с целью возможного упрощения моделей.

При построении рис. 1 и 4 использовалась подпрограмма <phaseportrait> математического пакета Maple 17.