Журнал вычислительной математики и математической физики, 2022, T. 62, № 2, стр. 305-319

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

А. М. Блохин^1, 2, *, Б. В. Семисалов^1, 2, **

¹ НГУ
630090 Новосибирск, ул. Пирогова, 1, Россия

² ИМ СО РАН
630090 Новосибирск, пр-т Коптюга, 4, Россия

^* E-mail: blokhin@math.nsc.ru
^** E-mail: vibis@ngs.ru

Поступила в редакцию 20.08.2020
После доработки 20.08.2020
Принята к публикации 17.09.2021

EDN: VHTTEL
DOI: 10.31857/S0044466922020065

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

Аннотация

Проведен численный анализ процесса установления стационарных течений несжимаемой вязкоупругой полимерной жидкости в канале с прямоугольным сечением под действием постоянного перепада давления. Для описания течений применяется реологическая мезоскопическая модель Покровского–Виноградова. При использовании интерполяций с узлами Чебышёва по пространственным переменным и неявной схемы по времени разработан алгоритм решения начально-краевых задач для нестационарных уравнений модели. Аналитически показано, что в стационарном случае модель допускает три решения высокой гладкости. Вопрос о том, какое из этих решений реализуется на практике, исследован с помощью расчетов предельного решения нестационарных уравнений. Установлено, что предельное решение с высокой точностью совпадает с одним из трeх решений стационарной задачи, и рассчитаны значения параметров, при которых происходит переключение с одного решения на другое. Библ. 11. Фиг. 6. Табл. 3.

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

1. ВВЕДЕНИЕ

Работа посвящена численному исследованию устойчивости стационарных режимов течения несжимаемой вязкоупругой полимерной жидкости в канале с прямоугольным поперечным сечением. Такие течения по своим качественным особенностям схожи с решениями Пуазейля уравнений Навье–Стокса, для которых в некоторых простых случаях (например, течение в цилиндрической трубе) известны строгие результаты об экспоненциальной сходимости нестационарных решений к стационарным (см. [1]), что гарантирует устойчивость последних.

В отличие от ньютоновской жидкости течения растворов и расплавов полимерных материалов демонстрируют сложную градиентную зависимость элонгационной и сдвиговой вязкостей от скоростей деформации, эффекты запаздывания и зависимость потока от ориентации и размеров макромолекул (см. [2]–[4]). Вследствие этого получить в аналитическом виде решения уравнений, описывающих течения полимерной жидкости типа Пуазейля, и исследовать их устойчивость не представляется возможным. Далее для описания таких течений использована реологическая мезоскопическая модель Покровского–Виноградова (см. [2], [5]). Анализ квазилинейного уравнения, полученного в [6] для расчета стационарных течений в каналах под действием постоянного перепада давления, показал, что его коэффициенты и правая часть определяются неоднозначно. При определeнных условиях это приводит к возникновению трeх качественно различных решений высокой гладкости. В данной работе численно исследуется вопрос о том, какое из этих решений является предельным для исходной нестационарной постановки, т.е. какое из них может быть реализовано на практике.

Для этого в статье разработан алгоритм, основанный на приближениях без насыщения и подходе к аппроксимации двумерных краевых задач для уравнений эллиптического типа, предложенном в [7]. При использовании этого алгоритма в совокупности с неявной схемой по временной переменной и задании начальных данных определенного вида получено предельное решение нестационарного уравнения, которое с высокой точностью совпадает с одним из трех решений стационарной задачи. Причем при определенных значениях параметров реализуется “переключение” предельного течения с одного решения стационарной задачи на другое. Таким образом, удалось прояснить структуру коэффициентов и правой части уравнения, описывающего устойчивые стационарные течения.

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

2. ПОСТАНОВКА ЗАДАЧИ

Следуя монографиям [1]–[3], [8], [9] и статье [5], запишем уравнения реологической мезоскопической модели Покровского–Виноградова, которая описывает течения несжимаемой вязкоупругой полимерной жидкости. Эти уравнения в безразмерной форме в прямоугольной декартовой системе координат $({{x}_{1}},{{x}_{2}},{{x}_{3}})$ выглядят следующим образом:

(1)

$\operatorname{div} {\mathbf{u}} = 0,$

(2)

$\frac{{d{\mathbf{u}}}}{{dt}} + \nabla P = \frac{1}{{{\text{Re}}}}\operatorname{div} (\Pi ),$

(3)

$\frac{{d{{a}_{{ij}}}}}{{dt}} - \sum\limits_{l = 1}^3 \frac{{\partial {{u}_{i}}}}{{\partial {{x}_{l}}}}{{a}_{{lj}}} - \sum\limits_{l = 1}^3 \frac{{\partial {{u}_{j}}}}{{\partial {{x}_{l}}}}{{a}_{{li}}} - \frac{1}{W}\left( {\frac{{\partial {{u}_{i}}}}{{\partial {{x}_{j}}}} + \frac{{\partial {{u}_{j}}}}{{\partial {{x}_{i}}}}} \right) + {{\mathcal{L}}_{{ij}}} = 0,\quad i = 1,2,3,\quad j = i,...,3.$

Здесь $t$ – время; ${{u}_{1}}$, ${{u}_{2}}$, ${{u}_{3}}$ – компоненты вектора скорости ${\mathbf{u}}$;

– $P$ – давление;

– ${{a}_{{ij}}},i,j = 1,2,3$, – компоненты симметрического тензора анизотропии $\Pi $ второго ранга;

– ${{{\mathbf{a}}}_{1}}$, ${{{\mathbf{a}}}_{2}}$, ${{{\mathbf{a}}}_{3}}$ – столбцы симметрической матрицы $\Pi = ({{a}_{{ij}}}) = ({{{\mathbf{a}}}_{1}},{{{\mathbf{a}}}_{2}},{{{\mathbf{a}}}_{3}})$;

– ${{\left\| {{{{\mathbf{a}}}_{{\mathbf{i}}}}} \right\|}^{2}} = ({{{\mathbf{a}}}_{{\mathbf{i}}}},{{{\mathbf{a}}}_{{\mathbf{i}}}})$, $i = 1,2,3$;

– $\operatorname{div} (\Pi ) = {{(\operatorname{div} ({{{\mathbf{a}}}_{1}}),\operatorname{div} ({{{\mathbf{a}}}_{2}}),\operatorname{div} ({{{\mathbf{a}}}_{3}}))}^{{\text{T}}}}$;

– ${{\mathcal{L}}_{{ij}}} = \left( {\tfrac{{{{a}_{{ij}}}}}{W} + \tfrac{{\bar {k}({{a}_{{11}}} + {{a}_{{22}}} + {{a}_{{33}}})}}{3}{{a}_{{ij}}} + \beta ({{{\mathbf{a}}}_{{\mathbf{i}}}},{{{\mathbf{a}}}_{{\mathbf{j}}}})} \right)$, $i,j = 1,2,3$;

– $\bar {k} = k - \beta $, $k$, $\beta $ ($0 < \beta < 1$) – феноменологические параметры модели, характеризующие вклады, связанные с анизотропией (величина $\beta $ учитывает ориентацию макромолекулярного клубка, число $k$ – его размеры, см. [5], [10]);

– ${\text{Re}} = \tfrac{{\varrho {{u}_{H}}l}}{{\eta _{0}^{ * }}}$ – число Рейнольдса; $\varrho ( = {\text{const}})$ – плотность среды;

– $W = \tfrac{{\tau _{0}^{ * }{{u}_{H}}}}{l}$ – число Вейсенберга;

– $\eta _{0}^{ * }$, $\tau _{0}^{ * }$ – начальные значения сдвиговой вязкости и времени релаксации при $T = {{T}_{0}}$, где ${{T}_{0}}$ – температура окружающей среды.

Далее для координат будем использовать обозначения $x,y,z$ ($x = {{x}_{1}}$, $y = {{x}_{2}}$, $z = {{x}_{3}}$).

Замечание 1. В [5], [10] рассмотрены случаи $k = \beta $, $k = 1.2\beta $. В настоящей работе исследуется более общий случай: $k = {{c}_{k}}\beta $, где ${{c}_{k}} > 0$ – новый феноменологический параметр.

Система (1)–(3) записана в безразмерном виде: время $t$; координаты $x$, $y$, $z$; компоненты вектора скорости ${{u}_{1}}$, ${{u}_{2}}$, ${{u}_{3}}$; давление $P$; компоненты тензора анизотропии ${{a}_{{ij}}}$ отнесены к $\tfrac{l}{{{{u}_{H}}}}$; $l$; ${{u}_{H}}$; $\varrho u_{H}^{2}$; $W{\text{/}}3$ соответственно, где $l$ – характерная длина, ${{u}_{H}}$ – характерная скорость.

Пусть полимерная жидкость течет в канале $\Omega $ с прямоугольным поперечным сечением, лежащим в плоскости $(y,z)$ (фиг. 1):

$\Omega = \{ (y,z): - 1 \leqslant y \leqslant 1,\; - {\kern 1pt} r \leqslant z \leqslant r\} .$

В канале вдоль оси $x$ действует перепад давления $\Delta P(t)$. Будем искать частное решение исходной системы (1)–(3) следующего вида:

(4)

$\begin{gathered} {{u}_{2}} = {{u}_{3}} \equiv 0,\quad {{u}_{1}} = u(t,y,z), \hfill \\ P = \mathcal{P}(t,y,z) - A(t)x, \hfill \\ {{a}_{{ij}}} = {{a}_{{ij}}}(t,y,z),\quad i,j = 1,2,3. \hfill \\ \end{gathered} $

Величина $A(t) = \tfrac{{\Delta P(t)}}{{\varrho u_{H}^{2}h}}$ есть безразмерный перепад давления на отрезке по $x$ длины $h$, $\mathcal{P}(t,y,z)$ – компонента давления, зависящая только от времени и от координаты точки сечения.

Фиг. 1.

Канал с прямоугольным сечением $\Omega $.

Пусть

${{\alpha }_{{ij}}} = \tfrac{{{{a}_{{ij}}}}}{{{\text{Re}}}},\quad i,j = 1,2,3,\quad {{\alpha }_{i}} = {{\alpha }_{{ii}}} + {{\unicode{230} }^{2}},\quad {{\unicode{230} }^{2}} = \tfrac{1}{{W{\text{Re}}}};$

${{K}_{I}} = \operatorname{Re} \left( {{{\unicode{230} }^{2}} + \tfrac{{\bar {k}}}{3}I} \right),\quad I = {{\alpha }_{{11}}} + {{\alpha }_{{22}}} + {{\alpha }_{{33}}},\quad {{\tilde {K}}_{I}} = {{K}_{I}} + \beta \operatorname{Re} I.$

Тогда с учетом (4) запишем систему (1)–(3) в виде двух систем:

(5)

$\begin{gathered} {{u}_{t}} - {{({{\alpha }_{{12}}})}_{y}} - {{({{\alpha }_{{13}}})}_{z}} = A(t), \hfill \\ {{({{\alpha }_{{12}}})}_{t}} - {{\alpha }_{2}}{{u}_{y}} - {{\alpha }_{{23}}}{{u}_{z}} + {{{\tilde {K}}}_{I}}{{\alpha }_{{12}}} + \beta \operatorname{Re} ({{\alpha }_{{13}}}{{\alpha }_{{23}}} - {{\alpha }_{{12}}}{{\alpha }_{{33}}}) = 0, \hfill \\ {{({{\alpha }_{{13}}})}_{t}} - {{\alpha }_{{23}}}{{u}_{y}} - {{\alpha }_{3}}{{u}_{z}} + {{{\tilde {K}}}_{I}}{{\alpha }_{{13}}} + \beta \operatorname{Re} ({{\alpha }_{{12}}}{{\alpha }_{{23}}} - {{\alpha }_{{13}}}{{\alpha }_{{22}}}) = 0, \hfill \\ \end{gathered} $

(6)

$\begin{gathered} {{({{\alpha }_{{11}}})}_{t}} - 2({{\alpha }_{{12}}}{{u}_{y}} + {{\alpha }_{{13}}}{{u}_{z}}) + {{K}_{I}}{{\alpha }_{{11}}} + \beta \operatorname{Re} (\alpha _{{11}}^{2} + \alpha _{{12}}^{2} + \alpha _{{13}}^{2}) = 0, \hfill \\ {{({{\alpha }_{{22}}})}_{t}} + {{{\tilde {K}}}_{I}}{{\alpha }_{{22}}} + \beta \operatorname{Re} (\alpha _{{12}}^{2} - {{\alpha }_{{11}}}{{\alpha }_{{22}}} + \alpha _{{23}}^{2} - {{\alpha }_{{22}}}{{\alpha }_{{33}}}) = 0, \hfill \\ {{({{\alpha }_{{33}}})}_{t}} + {{{\tilde {K}}}_{I}}{{\alpha }_{{33}}} + \beta \operatorname{Re} (\alpha _{{13}}^{2} - {{\alpha }_{{11}}}{{\alpha }_{{33}}} + \alpha _{{23}}^{2} - {{\alpha }_{{22}}}{{\alpha }_{{33}}}) = 0, \hfill \\ {{({{\alpha }_{{23}}})}_{t}} + {{{\tilde {K}}}_{I}}{{\alpha }_{{23}}} + \beta \operatorname{Re} ({{\alpha }_{{12}}}{{\alpha }_{{13}}} - {{\alpha }_{{11}}}{{\alpha }_{{23}}}) = 0. \hfill \\ \end{gathered} $

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

${{\mathcal{P}}_{y}} = {{({{\alpha }_{{22}}})}_{y}} + {{({{\alpha }_{{23}}})}_{z}},\quad {{\mathcal{P}}_{z}} = {{({{\alpha }_{{23}}})}_{y}} + {{({{\alpha }_{{33}}})}_{z}}.$

Дифференцируя первое уравнение (5) по $t$, второе уравнение по $y$ и третье уравнение по $z$, а затем складывая их, приходим к уравнению второго порядка для функции $u(t,y,z)$:

(7)

${{u}_{{tt}}} - \tilde {\Delta }u - \tilde {A}(t,y,z){{u}_{y}} - \tilde {B}(t,y,z){{u}_{z}} + {{\tilde {K}}_{I}}[{{u}_{t}} - A(t)] + F(t,y,z) = A{\kern 1pt} '(t),$

где

$\tilde {\Delta }u = {{\alpha }_{2}}{{u}_{{yy}}} + 2{{\alpha }_{{23}}}{{u}_{{yz}}} + {{\alpha }_{3}}{{u}_{{zz}}},\quad \tilde {A}(t,y,z) = {{({{\alpha }_{2}})}_{y}} + {{({{\alpha }_{{23}}})}_{z}},\quad \tilde {B}(t,y,z) = {{({{\alpha }_{{23}}})}_{y}} + {{({{\alpha }_{3}})}_{z}},$

$\begin{gathered} F(t,y,z) = {{({{{\tilde {K}}}_{I}})}_{y}}{{\alpha }_{{12}}} + {{({{{\tilde {K}}}_{I}})}_{z}}{{\alpha }_{{13}}} + \beta \operatorname{Re} \text{[}{{({{l}_{{22}}})}_{z}} + {{({{l}_{{33}}})}_{y}}], \\ {{l}_{{22}}} = {{\alpha }_{{12}}}{{\alpha }_{{23}}} - {{\alpha }_{{13}}}{{\alpha }_{{22}}},\quad {{l}_{{33}}} = {{\alpha }_{{13}}}{{\alpha }_{{23}}} - {{\alpha }_{{12}}}{{\alpha }_{{33}}}. \\ \end{gathered} $

Дополним уравнения (5)–(7) граничными условиями (условиями прилипания жидкости для скорости)

(8)

$u(t,y,z) = 0\quad {\text{при}}\quad t \geqslant 0,\quad (y,z) \in \partial \Omega ,$

и начальными данными в предположении, что в начальный момент времени $t = 0$ жидкость в канале покоится:

(9)

${{\alpha }_{{ij}}}(0,y,z) = 0,\quad u(0,y,z) = 0,\quad {{u}_{t}}(0,y,z) = 0\quad \forall (y,z) \in \Omega {\kern 1pt} .$

3. СТАЦИОНАРНЫЙ СЛУЧАЙ

Следуя [5] (см. также [6]), заметим, что при реализации условий (4) в стационарном случае, когда переменные задачи не зависят от времени: $u = u(y,z)$, ${{a}_{{ij}}} = {{a}_{{ij}}}(y,z)$, $i,j = \overline {1,3} $, $\mathcal{P} = \mathcal{P}(y,z)$, $A$ = = const, система уравнений (1)–(3) в канале с прямоугольным сечением сводится к одному разрешающему уравнению для определения продольной скорости течения:

(10)

$\hat {a}{{u}_{{yy}}} - 2\hat {b}{{u}_{{yz}}} + \hat {c}{{u}_{{zz}}} = - \operatorname{Re} A\tilde {\mathcal{K}},$

где

$\hat {a} = 1 - u_{y}^{2}\mathcal{L},\quad \hat {c} = 1 - u_{z}^{2}\mathcal{L},\quad \hat {b} = {{u}_{y}}{{u}_{z}}\mathcal{L},$

$\mathcal{L} = \frac{{1 - \Delta }}{{{{\lambda }^{2}}}},\quad \Delta = \frac{\lambda }{{\mu ({{{\tilde {\mathcal{K}}}}_{\mu }}\mu + \tilde {\mathcal{K}})}},\quad \lambda = \sqrt {u_{y}^{2} + u_{z}^{2}} ,$

$\tilde {\mathcal{K}} = \frac{{1 + W\sigma + 2W(\bar {k}{\text{/}}3 + \beta )(\sigma (1 + W\sigma ) + W{{\mu }^{2}})}}{{{{{(1 + W\sigma )}}^{2}}}},\quad \sigma = {{a}_{{22}}} + {{a}_{{33}}},\quad \mu = \sqrt {a_{{12}}^{2} + a_{{13}}^{2}} ,$

${{\tilde {\mathcal{K}}}_{\mu }} = \frac{{W{{\sigma }_{\mu }}}}{{{{{(1 + W\sigma )}}^{2}}}}[2(\bar {k}{\text{/}}3 + \beta ) - 1] + \frac{{4{{W}^{2}}\mu }}{{{{{(1 + W\sigma )}}^{3}}}}(\bar {k}{\text{/}}3 + \beta )[W\sigma - W\mu {{\sigma }_{\mu }} + 1]{\kern 1pt} .$

Используя формулы для $\tilde {\mathcal{K}}$ и выражение $\sigma (\mu )$, величину $\mu $ можно определить как решение нелинейного уравнения

$\tilde {\mathcal{K}}(\mu )\mu = \lambda .$

Функция $\sigma = \sigma (\mu )$ является решением кубического уравнения

(11)

$\begin{gathered} a{{\sigma }^{3}} + b{{\sigma }^{2}} + c\sigma + d = 0,\quad \sigma \ne - 1{\text{/}}W, \\ a = \widetilde \beta W,\quad b = \widetilde \beta + 1,\quad c = {{W}^{{ - 1}}} + W{{\mu }^{2}}\widetilde \beta ,\quad d = \beta {{\mu }^{2}},\quad \widetilde \beta = \frac{{2\bar {k}}}{3} + \beta . \\ \end{gathered} $

Анализ этого уравнения, проведенный в [6], привел к заключению, что в прямоугольной области при достаточно малых градиентах давления $A$ это уравнение имеет три действительных решения – функции ${{\sigma }_{1}}(\mu )$, ${{\sigma }_{2}}(\mu )$, ${{\sigma }_{3}}(\mu )$. Приведем далее выражения для этих решений. Следуя формулам Кардано для решения кубических уравнений, введем обозначения

$p = \frac{{3ac - {{b}^{2}}}}{{3{{a}^{2}}}},\quad q = \frac{{2{{b}^{3}} - 9abc + 27{{a}^{2}}d}}{{27{{a}^{3}}}},\quad \mathbb{Q} = {{(p{\text{/}}3)}^{3}} + {{(q{\text{/}}2)}^{2}}.$

После преобразований получаем $p = {{\mu }^{2}} - \mathbb{A}$, $q = \mathbb{B}{{\mu }^{2}} + \mathbb{C}$, где

$\mathbb{A} = \frac{{1 - \widetilde \beta + \mathop {\widetilde \beta }\nolimits^2 }}{{3\mathop {\widetilde \beta }\nolimits^2 {{W}^{2}}}},\quad \mathbb{B} = \frac{{3\beta - \widetilde \beta - 1}}{{3W\widetilde \beta }},\quad \mathbb{C} = \frac{{(1 + \widetilde \beta )(\widetilde \beta - 2)(2\widetilde \beta - 1)}}{{27{{W}^{3}}\mathop {\widetilde \beta }\nolimits^3 }}.$

Указанные три решения имеют вид

${{\sigma }_{1}} = 2cos\alpha \left| {\sqrt[6]{{\frac{{{{q}^{2}}}}{4} - \mathbb{Q}}}} \right| - \frac{b}{{3a}},\quad {{\sigma }_{{2,3}}} = - \left| {\sqrt[6]{{\frac{{{{q}^{2}}}}{4} - \mathbb{Q}}}} \right|( - cos\alpha \pm \sqrt 3 sin\alpha ) - \frac{b}{{3a}},\quad \alpha = \frac{{arctan(2\sqrt { - \mathbb{Q}} {\text{/}}q)}}{3}.$

При этом выражения для производных $\tfrac{{d{{\sigma }_{{1,2,3}}}}}{{d\mu }}$, фигурирующих в формуле для ${{\tilde {\mathcal{K}}}_{\mu }}$, имеют вид

$\frac{{\partial {{\sigma }_{1}}}}{{\partial \mu }} = \sigma _{1}^{'} = \frac{{{{2}^{{2/3}}}}}{{3\sqrt { - \mathbb{Q}} \left| {{{{({{q}^{2}} - 4\mathbb{Q})}}^{{5/6}}}} \right|}}((qsin\alpha - 2\sqrt { - Q} cos\alpha )\mathbb{Q}{\kern 1pt} {\text{'}} + (q\sqrt { - \mathbb{Q}} cos\alpha - 2\mathbb{Q}sin\alpha )q{\kern 1pt} '),$

$\begin{gathered} \frac{{\partial {{\sigma }_{{2,3}}}}}{{\partial \mu }} = \sigma _{{2,3}}^{'} = \frac{{ - 1}}{{3\sqrt[3]{2}\sqrt { - \mathbb{Q}} \left| {{{{({{q}^{2}} - 4\mathbb{Q})}}^{{5/6}}}} \right|}}( \pm (\sqrt 3 cos\alpha \pm sin\alpha )(2\mathbb{Q}q{\kern 1pt} {\text{'}} - q\mathbb{Q}{\kern 1pt} ') + \\ + \;(cos\alpha \mp \sqrt 3 sin\alpha )(2\mathbb{Q}{\kern 1pt} {\text{'}} - qq{\kern 1pt} ')\sqrt { - \mathbb{Q}} ), \\ \end{gathered} $

где для ${{\sigma }_{2}}$ выбираем верхний из знаков “$ \pm $", "$ \mp $", для ${{\sigma }_{3}}$ – нижний;

$q{\kern 1pt} ' = \frac{{\partial q}}{{\partial \mu }} = - 2\mu \frac{{1 - 3\beta + \widetilde \beta }}{{3\widetilde \beta W}},$

$\begin{gathered} \mathbb{Q}{\kern 1pt} ' = \frac{{\partial \mathbb{Q}}}{{\partial \mu }} = \frac{\mu }{{27{{{\tilde {\beta }}}^{4}}{{W}^{4}}}}\left\{ {27{{\beta }^{2}}{{{\tilde {\beta }}}^{2}}{{\mu }^{2}}{{W}^{2}} + \widetilde \beta (1 + \widetilde \beta {{\mu }^{2}}{{W}^{2}})[\widetilde \beta (\widetilde \beta (6{{\mu }^{2}}{{W}^{2}} - 1) + 4) - 1] - } \right. \\ - \;\left. {\beta (1 + \widetilde \beta )[\widetilde \beta (2\widetilde \beta (9{{\mu }^{2}}{{W}^{2}} - 1) + 5) - 2]} \right\}. \\ \end{gathered} $

Стационарное уравнение (10) дополняется граничными условиями

(12)

$u(y,z) = 0\quad {\text{при}}\quad (y,z) \in \partial \Omega .$

В [5] описан алгоритм решения краевой задачи (10), (12), основанный на методе установления, приближениях без насыщения и методе коллокаций (см. также [7]). Этот алгоритм использован в настоящей работе для поиска решений стационарной задачи.

Возникает естественный вопрос: какое из трех возможных решений задачи (10), (12) реализуется на практике? Нестационарная формулировка модели, полученная в разд. 2, имеет целью исследование этого вопроса. Для этого в модели будет задан постоянный градиент давления и реализованы вычисления вплоть до выхода на стационарный режим течения при некотором $t = {{t}_{s}}$. Сопоставление скорости $u({{t}_{s}},y,z)$ нестационарного течения в момент времени t = t_s и решения задачи (10), (12), полученного при использовании разных зависимостей $\sigma (\mu )$, покажет, какой режим реализуется на практике.

4. ОПИСАНИЕ ВЫЧИСЛИТЕЛЬНОГО АЛГОРИТМА

Для реализации расчетов введем сетку по времени с шагом $\tau $ и узлами ${{t}_{n}} = n\tau $, $n = 0,1,...$ . Обозначим ${{u}^{n}} = {{u}^{n}}(y,z) = u({{t}_{n}},y,z)$, $\alpha _{{ij}}^{n} = \alpha _{{ij}}^{n}(y,z) = {{\alpha }_{{ij}}}({{t}_{n}},y,z)$ и аппроксимируем в уравнениях (5)–(7) производные по времени разностными отношениями вида

${{u}_{t}} \approx \frac{{{{u}^{{n + 1}}} - {{u}^{n}}}}{\tau },\quad {{u}_{{tt}}} \approx \frac{{{{u}^{{n + 1}}} - 2{{u}^{n}} + {{u}^{{n - 1}}}}}{{{{\tau }^{2}}}},\quad n = 1,2,...\;.$

4.1. Линеаризация уравнений модели

Для реализации вычислений необходимо линеаризовать по Ньютону уравнения (5)–(7) относительно неизвестных функций ${{\alpha }_{{ij}}}$ и $u$ и на каждом шаге по времени осуществлять итерации по нелинейности. В расчетах невязка итераций по нелинейности была установлена равной ${{10}^{{ - 10}}}$. Запишем линеаризованные уравнения для функций ${{\alpha }_{{ij}}}$ в виде

$\frac{{\alpha _{{11}}^{{n + 1}} - \alpha _{{11}}^{n}}}{\tau } = - {\text{Re}}\left[ {\left( {\frac{{\bar {k}}}{3} + 2\beta } \right)\alpha _{{11}}^{n} + \frac{{\bar {k}}}{3}{{I}^{n}} + {{\unicode{230} }^{2}}} \right]\alpha _{{11}}^{{n + 1}} + $

$ + \;2(\alpha _{{12}}^{n}u_{y}^{n} + \alpha _{{13}}^{n}u_{z}^{n}) + \frac{{\operatorname{Re} \bar {k}}}{3}{{(\alpha _{{11}}^{n})}^{2}} - \beta \operatorname{Re} [{{(\alpha _{{12}}^{n})}^{2}} + {{(\alpha _{{13}}^{n})}^{2}} - {{(\alpha _{{11}}^{n})}^{2}}],$

$\frac{{\alpha _{{22}}^{{n + 1}} - \alpha _{{22}}^{n}}}{\tau } = - \left[ {K_{I}^{n} + {\text{Re}}\left( {\frac{{\bar {k}}}{3} + 2\beta } \right)\alpha _{{22}}^{n}} \right]\alpha _{{22}}^{{n + 1}} + {\text{Re}}\left( {\frac{{\bar {k}}}{3}{{{(\alpha _{{22}}^{n})}}^{2}} + \beta [{{{(\alpha _{{22}}^{n})}}^{2}} - {{{(\alpha _{{12}}^{n})}}^{2}} - {{{(\alpha _{{23}}^{n})}}^{2}}]} \right),$

(13)

$\frac{{\alpha _{{33}}^{{n + 1}} - \alpha _{{33}}^{n}}}{\tau } = - \left[ {K_{I}^{n} + {\text{Re}}\left( {\frac{{\bar {k}}}{3} + 2\beta } \right)\alpha _{{33}}^{n}} \right]\alpha _{{33}}^{{n + 1}} + {\text{Re}}\left( {\frac{{\bar {k}}}{3}{{{(\alpha _{{33}}^{n})}}^{2}} + \beta [{{{(\alpha _{{33}}^{n})}}^{2}} - {{{(\alpha _{{13}}^{n})}}^{2}} - {{{(\alpha _{{23}}^{n})}}^{2}}]} \right),$

$\frac{{\alpha _{{12}}^{{n + 1}} - \alpha _{{12}}^{n}}}{\tau } = [\beta \operatorname{Re} \alpha _{{33}}^{n} - \tilde {K}_{I}^{n}]{\kern 1pt} \alpha _{{12}}^{{n + 1}} + \alpha _{2}^{n}u_{y}^{n} + \alpha _{{23}}^{n}u_{z}^{n} - \beta \operatorname{Re} \alpha _{{13}}^{n}\alpha _{{23}}^{n},$

$\frac{{\alpha _{{13}}^{{n + 1}} - \alpha _{{13}}^{n}}}{\tau } = [\beta \operatorname{Re} \alpha _{{22}}^{n} - \tilde {K}_{I}^{n}]\alpha _{{13}}^{{n + 1}} + \alpha _{{23}}^{n}u_{y}^{n} + \alpha _{3}^{n}u_{z}^{n} - \beta \operatorname{Re} \alpha _{{12}}^{n}\alpha _{{23}}^{n},$

$\frac{{\alpha _{{23}}^{{n + 1}} - \alpha _{{23}}^{n}}}{\tau } = [\beta \operatorname{Re} \alpha _{{11}}^{n} - \tilde {K}_{I}^{n}]{\kern 1pt} \alpha _{{23}}^{{n + 1}} - \beta \operatorname{Re} \alpha _{{12}}^{n}\alpha _{{13}}^{n},$

где ${{I}^{n}} = a_{{11}}^{n} + a_{{22}}^{n} + a_{{33}}^{n}$, $K_{I}^{n} = {\text{Re}}\left( {{{\unicode{230} }^{2}} + \tfrac{{\bar {k}}}{3}{{I}^{n}}} \right)$, $\tilde {K}_{I}^{n} = K_{I}^{n} + \beta \operatorname{Re} {{I}^{n}}$, $\alpha _{i}^{n} = \alpha _{{ii}}^{n} + {{\unicode{230} }^{2}}$, $i = 2,3$, $u_{y}^{n}$, $u_{z}^{n}$ – производные функции $u(t,y,z)$ при $t = {{t}_{n}}$.

Для поиска функции ${{u}^{{n + 1}}}(y,z)$ из (7) несложно вывести следующее линеаризованное уравнение:

(14)

$\begin{gathered} \frac{{{{u}^{{n + 1}}} - 2{{u}^{n}} + {{u}^{{n - 1}}}}}{{{{\tau }^{2}}}} + \operatorname{Re} {{\unicode{230} }^{2}}\left( {\frac{{{{u}^{{n + 1}}} - {{u}^{n}}}}{\tau }} \right) - {{\unicode{230} }^{2}}(u_{{yy}}^{{n + 1}} + u_{{zz}}^{{n + 1}}) = \alpha _{{22}}^{n}u_{{yy}}^{n} + \alpha _{{33}}^{n}u_{{zz}}^{n} + \\ + \;2\alpha _{{23}}^{n}u_{{yz}}^{n} - {\text{Re}}\left( {\frac{{\bar {k}}}{3} + \beta } \right){{I}^{n}}\left( {\frac{{{{u}^{n}} - {{u}^{{n - 1}}}}}{\tau }} \right) + \tilde {K}_{I}^{n}A({{t}_{n}}) + {{F}^{n}}, \\ \end{gathered} $

где

${{F}^{n}} = [{{(\alpha _{2}^{n})}_{y}} + {{(\alpha _{{23}}^{n})}_{z}}]{\kern 1pt} u_{y}^{n} + [{{(\alpha _{{23}}^{n})}_{y}} + {{(\alpha _{3}^{n})}_{z}}]u_{z}^{n} - {\text{Re}}\left( {\frac{{\bar {k}}}{3} + \beta } \right)[I_{y}^{n}\alpha _{{12}}^{n} + I_{z}^{n}\alpha _{{13}}^{n}] - \beta {\text{Re}}[{{(l_{{22}}^{n})}_{z}} + {{(l_{{33}}^{n})}_{y}}] + A{\kern 1pt} '({{t}_{n}}),$

$I_{{y,z}}^{n} = {{(\alpha _{{11}}^{n})}_{{y,z}}} + {{(\alpha _{{22}}^{n})}_{{y,z}}} + {{(\alpha _{{33}}^{n})}_{{y,z}}}$ – производные от функции ${{I}^{n}}$, ${{(l_{{22}}^{n})}_{z}} = {{(\alpha _{{12}}^{n}\alpha _{{23}}^{n} - \alpha _{{13}}^{n}\alpha _{{22}}^{n})}_{z}}$, ${{(l_{{33}}^{n})}_{y}} = {{(\alpha _{{13}}^{n}\alpha _{{23}}^{n} - \alpha _{{12}}^{n}\alpha _{{33}}^{n})}_{y}}$, $A{\kern 1pt} '({{t}_{n}}) = {{\left. {\tfrac{{dA(t)}}{{dt}}} \right|}_{{t = {{t}_{n}}}}}$. Отметим также, что при записи (14) использована аппроксимация

${{\tilde {K}}_{I}}{{u}_{t}} = \operatorname{Re} {{\unicode{230} }^{2}}{{u}_{t}} + {\text{Re}}\left( {\frac{{\bar {k}}}{3} + \beta } \right){{I}^{n}}{{u}_{t}} \approx \operatorname{Re} {{\unicode{230} }^{2}}\left( {\frac{{{{u}^{{n + 1}}} - {{u}^{n}}}}{\tau }} \right) + {\text{Re}}\left( {\frac{{\bar {k}}}{3} + \beta } \right){{I}^{n}}\left( {\frac{{{{u}^{n}} - {{u}^{{n - 1}}}}}{\tau }} \right).$

4.2. Аппроксимация неизвестных функций и их производных

На каждом шаге итерационной схемы по времени будем искать решения уравнений (13), (14) с граничными условиями (8) в классе достаточно гладких в области $\Omega $ функций. Далее для записи вычислительной схемы необходимо аппроксимировать в уравнениях (13), (14) функции ${{u}^{n}}(y,z)$ и $\alpha _{{ij}}^{n}(y,z)$ и их производные по переменным $y$ и $z$ элементами конечномерных функциональных пространств. В качестве таких пространств используем пространства алгебраических полиномов, а в качестве способа аппроксимации – интерполяционные полиномы в форме Лагранжа с узлами в нулях многочленов Чебышёва. Конкретно, для приближения ${{u}^{n}}(y,z)$ и $\alpha _{{ij}}^{n}(y,z)$ будем использовать прямые (тензорные) произведения таких полиномов в области $\Omega $. Свойства таких приближений, а именно отсутствие их насыщения, логарифмический рост констант Лебега, а также оценки погрешности подробно обсуждаются в [11, гл. 3].

Введем в области $\Omega $ сетку с узлами $({{y}_{k}},{{z}_{m}})$, где

${{y}_{k}} = cos\tfrac{{(2k - 1)\pi }}{{2K}}$, ${{z}_{m}} = rcos\tfrac{{(2m - 1)\pi }}{{2M}}$, $k = \overline {1,K} $, $m = \overline {1,M} $.

Сделаем для удобства преобразование координаты $z$: $\tilde {z} = z{\text{/}}r$, $\mathop {\tilde {z}}\nolimits_m = {{z}_{m}}{\text{/}}r$, где $\tilde {z},{{\tilde {z}}_{m}} \in [ - 1,1]$. Пусть ${{(\alpha _{{ij}}^{n})}_{{km}}} = \alpha _{{ij}}^{n}({{y}_{k}},{{z}_{m}})$, $u_{{km}}^{n} = {{u}^{n}}({{y}_{k}},{{z}_{m}})$, тогда

(15)

$\alpha _{{ij}}^{n}(y,z) \approx P(\alpha _{{ij}}^{n},y,z) = \sum\limits_{k = 1}^K \,\sum\limits_{m = 1}^M \,{{l}_{{kK}}}(y){{l}_{{mM}}}(\tilde {z}){{(\alpha _{{ij}}^{n})}_{{km}}},$

(16)

${{u}^{n}}(y,z) \approx {{P}_{\omega }}({{u}^{n}},y,z) = \sum\limits_{k = 1}^K \,\sum\limits_{m = 1}^M \,\frac{{\omega (y,\tilde {z})}}{{\omega ({{y}_{k}},\mathop {\tilde {z}}\nolimits_m )}}{{l}_{{kK}}}(y){{l}_{{mM}}}(\tilde {z})u_{{km}}^{n},$

где

${{l}_{{kK}}}(y) = \frac{{{{T}_{K}}(y)}}{{(y - {{y}_{k}})T_{K}^{'}({{y}_{k}})}},\quad {{l}_{{mM}}}(\tilde {z}) = \frac{{{{T}_{M}}(\tilde {z})}}{{(\tilde {z} - {{{\tilde {z}}}_{m}})T_{M}^{'}({{{\tilde {z}}}_{m}})}}$

суть фундаментальные многочлены интерполяции, ${{T}_{K}}(y)$, ${{T}_{M}}(\tilde {z})$ – многочлены Чебышёва степеней $K$ и $M$ соответственно, множитель $\omega (y,\tilde {z}) = (1 - {{y}^{2}})(1 - {{\tilde {z}}^{2}})$ в формуле (16) обеспечивает автоматическое выполнение граничных условий (8) (см. также [11, гл. 9, § 5]).

Для аппроксимации уравнений (13), (14) далее будет использован метод коллокаций с узлами $({{y}_{k}},{{z}_{m}})$, поэтому нам потребуются значения производных полиномов (15), (16) в этих узлах. Приведем формулы дифференцирования по переменной $y$ (по переменной $z$ формулы строятся аналогично):

(17)

$\frac{{\partial P}}{{\partial y}} = \sum\limits_{k = 1}^K \,\sum\limits_{m = 1}^M \,{{S}_{{kK}}}(y)s_{k}^{2}\left\{ {\frac{{{{T}_{K}}(y)}}{{y - {{y}_{k}}}} - T_{K}^{'}(y)} \right\}{{l}_{{mM}}}(\tilde {z})({{\alpha }_{{ij}}})_{{km}}^{n},$

(18)

$\frac{{\partial {{P}_{\omega }}}}{{\partial y}} = \sum\limits_{k = 1}^K \,\sum\limits_{m = 1}^M \,{{( - 1)}^{{K - 1}}}{{S}_{{kK}}}(y){{\lambda }_{{mM}}}(\tilde {z})\left\{ {\frac{{2y{{y}_{k}} - {{y}^{2}} - 1}}{{y - {{y}_{k}}}}{{T}_{K}}(y) + (1 - {{y}^{2}})T_{K}^{'}(y)} \right\}u_{{km}}^{n}.$

Поскольку уравнение (14) имеет второй порядок, потребуются также аппроксимации второй производной от функции $u(y,z)$:

(19)

$\frac{{{{\partial }^{2}}{{P}_{\omega }}}}{{\partial {{y}^{2}}}} = \sum\limits_{k = 1}^K \,\sum\limits_{m = 1}^M \,{{( - 1)}^{{K - 1}}}{{S}_{{kK}}}(y){{\lambda }_{{mM}}}(\tilde {z})\left\{ {\frac{{2s_{k}^{2} - {{K}^{2}}{{{(y - {{y}_{k}})}}^{2}}}}{{{{{(y - {{y}_{k}})}}^{2}}}}{{T}_{K}}(y) - \frac{{{{y}^{2}} - 3y{{y}_{k}} + 2}}{{y - {{y}_{k}}}}T_{K}^{'}(y)} \right\}u_{{km}}^{n},$

где

${{s}_{k}} = \sqrt {1 - y_{k}^{2}} ,\quad {{S}_{{kK}}}(y) = \frac{{{{{( - 1)}}^{k}}}}{{K(y - {{y}_{k}}){{s}_{k}}}},\quad {{\lambda }_{{mM}}}(\tilde {z}) = \frac{{(1 - {{{\tilde {z}}}^{2}})}}{{1 - \tilde {z}_{m}^{2}}}{{l}_{{mM}}}(\tilde {z}).$

Переходя в (17)–(19) к пределу при $y \to {{y}_{l}}$, $z \to {{z}_{q}}$, $l = \overline {1,K} $, $q = \overline {1,M} $, используя свойство фундаментальных многочленов ${{l}_{{mM}}}(\mathop {\tilde {z}}\nolimits_q ) = {{\delta }_{{mq}}}$ – символ Кронекера и правило Лопиталя, находим

(20)

$\frac{{\partial \alpha _{{ij}}^{n}}}{{\partial y}}({{y}_{l}},{{z}_{q}}) \approx \mathop {lim}\limits_{(y,z) \to ({{y}_{l}},{{z}_{q}})} \frac{{\partial P(y,z)}}{{\partial y}} = \sum\limits_{k = 1,k \ne l}^K \frac{{{{{( - 1)}}^{{l + k}}}{{s}_{k}}}}{{{{s}_{l}}({{y}_{l}} - {{y}_{k}})}}({{\alpha }_{{ij}}})_{{kq}}^{n} - \frac{{{{y}_{l}}}}{{2s_{l}^{2}}}({{\alpha }_{{ij}}})_{{lq}}^{n},$

(21)

$\frac{{\partial {{u}^{n}}}}{{\partial y}}({{y}_{l}},{{z}_{q}}) \approx \mathop {lim}\limits_{(y,z) \to ({{y}_{l}},{{z}_{q}})} \frac{{\partial {{P}_{\omega }}(y,z)}}{{\partial y}} = \sum\limits_{k = 1,k \ne l}^K {{( - 1)}^{{k + l}}}\frac{{{{s}_{l}}}}{{{{s}_{k}}({{y}_{l}} - {{y}_{k}})}}u_{{kq}}^{n} - \frac{{3{{y}_{l}}}}{{2s_{l}^{2}}}u_{{lq}}^{n},$

(22)

$\frac{{{{\partial }^{2}}{{u}^{n}}}}{{\partial {{y}^{2}}}}({{y}_{l}},{{z}_{q}}) \approx \mathop {lim}\limits_{(y,z) \to ({{y}_{l}},\mathop {\tilde {z}}\nolimits_q )} \frac{{{{\partial }^{2}}{{P}_{\omega }}(y,z)}}{{\partial {{y}^{2}}}} = \sum\limits_{k = 1,k \ne l}^K {{( - 1)}^{{k + l - 1}}}\frac{{2s_{l}^{2} + 3{{y}_{l}}({{y}_{l}} - {{y}_{k}})}}{{{{s}_{k}}{{s}_{l}}{{{({{y}_{l}} - {{y}_{k}})}}^{2}}}}u_{{kq}}^{n} - \frac{{({{K}^{2}} + 5)s_{l}^{2} + 3y_{l}^{2}}}{{3s_{l}^{4}}}u_{{lq}}^{n}.$

Пусть

${{\xi }_{{lk}}} = \frac{{{{{( - 1)}}^{{l + k}}}{{s}_{k}}}}{{{{s}_{l}}({{y}_{l}} - {{y}_{k}})}},\quad {{\eta }_{{lk}}} = \frac{{{{{( - 1)}}^{{l + k}}}{{s}_{l}}}}{{{{s}_{k}}({{y}_{l}} - {{y}_{k}})}},\quad {{\zeta }_{{lk}}} = {{( - 1)}^{{l + k - 1}}}\frac{{2s_{l}^{2} + 3{{y}_{l}}({{y}_{l}} - {{y}_{k}})}}{{{{s}_{k}}{{s}_{l}}{{{({{y}_{l}} - {{y}_{k}})}}^{2}}}},$

где $l,k = \overline {1,K} ,l \ne k$,

${{\nu }_{l}} = - \frac{{{{y}_{l}}}}{{2s_{l}^{2}}},\quad {{\hat {\nu }}_{l}} = - \frac{{3{{y}_{l}}}}{{2s_{l}^{2}}},\quad {{\tilde {\nu }}_{l}} = - \frac{{({{K}^{2}} + 5)s_{l}^{2} + 3y_{l}^{2}}}{{3s_{l}^{4}}},\quad l = \overline {1,K} ,$

$\Upsilon _{{ij}}^{n}$, ${{(\Upsilon _{{ij}}^{n})}_{y}}$ суть $K \times M$-матрицы с элементами $\alpha _{{ij}}^{n}({{y}_{k}},{{z}_{m}})$ и $\tfrac{{\partial \alpha _{{ij}}^{n}}}{{\partial y}}({{y}_{k}},{{z}_{m}})$,

${{U}^{n}}$, $U_{y}^{n}$, $U_{{yy}}^{n}$ суть $K \times M$-матрицы с элементами ${{u}^{n}}({{y}_{k}},{{z}_{m}})$, $\tfrac{{\partial {{u}^{n}}}}{{\partial y}}({{y}_{k}},{{z}_{m}})$, $\tfrac{{{{\partial }^{2}}{{u}^{n}}}}{{\partial {{y}^{2}}}}({{y}_{k}},{{z}_{m}})$.

Сформируем $K \times K$-матрицы:

${{A}_{1}} = \left( {\begin{array}{*{20}{c}} {{{\nu }_{1}}}&{{{\xi }_{{12}}}}&{...}&{{{\xi }_{{1K}}}} \\ {{{\xi }_{{21}}}}&{{{\nu }_{2}}}&{...}&{{{\xi }_{{2K}}}} \\ \vdots & \vdots &{...}& \vdots \\ {{{\xi }_{{K1}}}}&{{{\xi }_{{K2}}}}&{...}&{{{\nu }_{K}}} \end{array}} \right),\quad {{\mathcal{A}}_{1}} = \left( {\begin{array}{*{20}{c}} {{{{\hat {\nu }}}_{1}}}&{{{\eta }_{{12}}}}&{...}&{{{\eta }_{{1K}}}} \\ {{{\eta }_{{21}}}}&{{{{\hat {\nu }}}_{2}}}&{...}&{{{\eta }_{{2K}}}} \\ \vdots & \vdots &{...}& \vdots \\ {{{\eta }_{{K1}}}}&{{{\eta }_{{K2}}}}&{...}&{\mathop {\hat {\nu }}\nolimits_K } \end{array}} \right),\quad {{\mathcal{A}}_{2}} = \left( {\begin{array}{*{20}{c}} {{{{\tilde {\nu }}}_{1}}}&{{{\zeta }_{{12}}}}&{...}&{{{\zeta }_{{1K}}}} \\ {{{\zeta }_{{21}}}}&{{{{\tilde {\nu }}}_{2}}}&{...}&{{{\zeta }_{{2K}}}} \\ \vdots & \vdots &{...}& \vdots \\ {{{\zeta }_{{K1}}}}&{{{\zeta }_{{K2}}}}&{...}&{{{{\tilde {\nu }}}_{K}}} \end{array}} \right).$

Переписывая линейные комбинации в правых частях (20)–(22) в виде произведений матриц, для аппроксимации производных по $y$ в уравнениях (13), (14) получим формулы

${{(\Upsilon _{{ij}}^{n})}_{y}} \approx {{A}_{1}}\Upsilon _{{ij}}^{n},\quad U_{y}^{n} \approx {{\mathcal{A}}_{1}}{{U}^{n}},\quad U_{{yy}}^{n} \approx {{\mathcal{A}}_{2}}{{U}^{n}}.$

Дифференцируя приближения (15), (16) по переменной $z$ и проводя выкладки, аналогичные описанным выше, сформируем $M \times M$-матрицы ${{B}_{1}}$, ${{\mathcal{B}}_{1}}$, ${{\mathcal{B}}_{2}}$ для аппроксимации производных по переменной $z$:

${{(\Upsilon _{{ij}}^{n})}_{z}} \approx \Upsilon _{{ij}}^{n}B_{1}^{{\text{T}}},\quad U_{z}^{n} \approx {{U}^{n}}\mathcal{B}_{1}^{{\text{T}}},\quad U_{{zz}}^{n} \approx {{U}^{n}}\mathcal{B}_{2}^{{\text{T}}},\quad {{(\Upsilon _{{ij}}^{n})}_{{yz}}} \approx {{A}_{1}}\Upsilon _{{ij}}^{n}B_{1}^{{\text{T}}},\quad U_{{yz}}^{n} \approx {{\mathcal{A}}_{1}}{{U}^{n}}\mathcal{B}_{1}^{{\text{T}}}.$

Здесь $K \times M$-матрицы ${{(\Upsilon _{{ij}}^{n})}_{z}}$, $U_{z}^{n}$, $U_{{zz}}^{n}$, ${{(\Upsilon _{{ij}}^{n})}_{{yz}}}$, $U_{{yz}}^{n}$ определяются аналогично ${{(\Upsilon _{{ij}}^{n})}_{y}}$, $U_{y}^{n}$.

Далее для построения алгоритмов воспользуемся спектральным разложением матриц ${{\mathcal{A}}_{2}}$, ${{\mathcal{B}}_{2}}$, аппроксимирующих вторые производные:

(23)

${{\mathcal{A}}_{2}} = {{R}_{\mathcal{A}}}{{D}_{\mathcal{A}}}R_{\mathcal{A}}^{{ - 1}},\quad {{\mathcal{B}}_{2}} = {{R}_{\mathcal{B}}}{{D}_{\mathcal{B}}}R_{\mathcal{B}}^{{ - 1}},$

где ${{R}_{\mathcal{A}}}$, ${{R}_{\mathcal{B}}}$ – матрицы собственных вектров ${{\mathcal{A}}_{2}}$ и ${{\mathcal{B}}_{2}}$; ${{D}_{\mathcal{A}}}$, ${{D}_{\mathcal{B}}}$ – диагональные матрицы собственных значений ${{\mathcal{A}}_{2}}$ и ${{\mathcal{B}}_{2}}$ – $d_{\mathcal{A}}^{k}$, $d_{\mathcal{B}}^{m}$.

Замечание 2. Строго говоря, при наличии у матриц ${{\mathcal{A}}_{2}}$ и ${{\mathcal{B}}_{2}}$ пар комплексно-сопряженных собственных чисел, матрицы ${{D}_{\mathcal{A}}}$, ${{D}_{\mathcal{B}}}$ являются блочно-диагональными, содержащими на диагонали блоки размера $2 \times 2$. Однако, как показали численные эксперименты, выбранный базис исключает такую возможность. Важным обстоятельством является также медленный рост обусловленностей матриц ${{R}_{\mathcal{A}}}$, ${{R}_{\mathcal{B}}}$ с ростом $K$, $M$, обеспечивающий устойчивость алгоритма к погрешностям округления (см. [7]).

4.3. Решение задачи линейной алгебры на каждом временном шаге

После подстановки описанных аппроксимаций в уравнения (13), (14) и отбрасывания погрешности приходим к линейным матричным уравнениям для выражения матриц $\Upsilon _{{ij}}^{{n + 1}}$, ${{U}^{{n + 1}}}$ через матрицы $\Upsilon _{{ij}}^{n}$, ${{U}^{n}}$. Запишем для примера одно из уравнений для выражения $\Upsilon _{{11}}^{n}$ (см. первое уравнение в системе (13)):

(24)

$\begin{gathered} {{\mathfrak{M}}_{{11}}} \cdot \Upsilon _{{11}}^{{n + 1}} = \left\{ {\frac{1}{\tau } + \operatorname{Re} \left[ {\left( {\frac{{\bar {k}}}{3} + 2\beta } \right)\Upsilon _{{11}}^{n} + \frac{{\bar {k}}}{3}(\Upsilon _{{11}}^{n} + \Upsilon _{{22}}^{n} + \Upsilon _{{33}}^{n}) + {{\unicode{230} }^{2}}{{E}_{I}}} \right]} \right\} \cdot \Upsilon _{{11}}^{{n + 1}} = \\ = \frac{1}{\tau }\Upsilon _{{11}}^{n} + 2(\Upsilon _{{12}}^{n} \cdot ({{\mathcal{A}}_{1}}{{U}^{n}}) + \Upsilon _{{13}}^{n} \cdot ({{U}^{n}}\mathcal{B}_{1}^{{\text{T}}})) + \frac{{\operatorname{Re} \bar {k}}}{3}(\Upsilon _{{11}}^{n} \cdot \Upsilon _{{11}}^{n}) - \beta \operatorname{Re} (\Upsilon _{{12}}^{n} \cdot \Upsilon _{{12}}^{n} + \Upsilon _{{13}}^{n} \cdot \Upsilon _{{13}}^{n} - \Upsilon _{{11}}^{n} \cdot \Upsilon _{{11}}^{n}). \\ \end{gathered} $

Здесь точка означает поэлементное произведение матриц, ${{E}_{I}}$ – матрица размера $K \times M$, все элементы которой равны единице, начальные данные задаются в соответствии с (9): $\Upsilon _{{11}}^{0} = {{(0)}_{{K \times M}}}$ – нулевая матрица. Аппроксимации остальных пяти уравнений из системы (13) записываются аналогично, при этом в левой части составляются матрицы ${{\mathfrak{M}}_{{ij}}}$.

Аппроксимация уравнения (14) выглядит следующим образом:

(25)

$\begin{gathered} \left( {\frac{1}{{{{\tau }^{2}}}} + \frac{{\operatorname{Re} {{\unicode{230} }^{2}}}}{\tau }} \right){{U}^{{n + 1}}} - {{\unicode{230} }^{2}}{{\mathcal{A}}_{2}}{{U}^{{n + 1}}} - {{\unicode{230} }^{2}}{{U}^{{n + 1}}}\mathcal{B}_{2}^{{\text{T}}} = \mathcal{F}_{U}^{n}: = \frac{1}{{{{\tau }^{2}}}}(2{{U}^{n}} - {{U}^{{n - 1}}}) + \frac{{\operatorname{Re} {{\unicode{230} }^{2}}}}{\tau }{{U}^{n}} + \Upsilon _{{22}}^{n} \cdot ({{\mathcal{A}}_{2}}{{U}^{n}}) + \hfill \\ + \;\Upsilon _{{33}}^{n} \cdot ({{U}^{n}}\mathcal{B}_{2}^{{\text{T}}}) + 2\Upsilon _{{23}}^{n} \cdot ({{A}_{1}}{{U}^{n}}B_{1}^{{\text{T}}}) - \frac{{\operatorname{Re} }}{\tau }\left( {\frac{{\bar {k}}}{3} + \beta } \right)(\Upsilon _{{11}}^{n} + \Upsilon _{{22}}^{n} + \Upsilon _{{33}}^{n}) \cdot ({{U}^{n}} - {{U}^{{n - 1}}}) + \tilde {\mathcal{K}}_{I}^{n}A({{t}_{n}}) + {{\mathcal{F}}^{n}}, \hfill \\ \end{gathered} $

где

$\tilde {\mathcal{K}}_{I}^{n} = \operatorname{Re} {{\unicode{230} }^{2}}{{E}_{I}} + {\text{Re}}\left( {\frac{{\bar {k}}}{3} + \beta } \right)(\Upsilon _{{11}}^{n} + \Upsilon _{{22}}^{n} + \Upsilon _{{33}}^{n}),$

$\begin{gathered} {{\mathcal{F}}^{n}} = [{{A}_{1}}\Upsilon _{{22}}^{n} + \Upsilon _{{23}}^{n}B_{1}^{{\text{T}}}] \cdot ({{\mathcal{A}}_{1}}{{U}^{n}}) + [{{A}_{1}}\Upsilon _{{23}}^{n} + \Upsilon _{{33}}^{n}B_{1}^{{\text{T}}}] \cdot ({{U}^{n}}\mathcal{B}_{1}^{{\text{T}}}) - \\ - \;{\text{Re}}\left( {\frac{{\bar {k}}}{3} + \beta } \right)[({{A}_{1}}(\Upsilon _{{11}}^{n} + \Upsilon _{{22}}^{n} + \Upsilon _{{33}}^{n})) \cdot \Upsilon _{{12}}^{n} + ((\Upsilon _{{11}}^{n} + \Upsilon _{{22}}^{n} + \Upsilon _{{33}}^{n})B_{1}^{{\text{T}}}) \cdot \Upsilon _{{13}}^{n}] - \\ - \;\beta {\text{Re}}{\kern 1pt} [(\Upsilon _{{12}}^{n} \cdot \Upsilon _{{23}}^{n} - \Upsilon _{{13}}^{n} \cdot \Upsilon _{{22}}^{n})B_{1}^{{\text{T}}} + {{A}_{1}}(\Upsilon _{{13}}^{n} \cdot \Upsilon _{{23}}^{n} - \Upsilon _{{12}}^{n} \cdot \Upsilon _{{33}}^{n})] + A{\kern 1pt} '({{t}_{n}}) \\ \end{gathered} $

суть $K \times N$-матрицы.

Для решения матричного уравнения (24) на каждом шаге по времени необходимо выполнять деления элементов матрицы в правой части на элементы матрицы ${{\mathfrak{M}}_{{11}}}$. Составляя по аналогии с (24) матричные уравнения для остальных пяти уравнений системы (13) и осуществляя поэлементные деления их правых и левых частей, получаем приближенные значения решений $\alpha _{{ij}}^{{n + 1}}$ в узлах коллокации.

Для решения (25) используем спектральные разложения (23) матриц ${{\mathcal{A}}_{2}}$, $\mathcal{B}_{2}^{{\text{T}}}$, стоящих в левых частях. Умножим (25) на матрицу $R_{\mathcal{A}}^{{ - 1}}$ слева и на матрицу $R_{\mathcal{B}}^{{ - 1}}$ справа, обозначим ${{V}^{{n + 1}}} = R_{\mathcal{A}}^{{ - 1}}{{U}^{{n + 1}}}R_{\mathcal{B}}^{{ - 1}}$, ${{\mathcal{G}}^{n}} = R_{\mathcal{A}}^{{ - 1}}\mathcal{F}_{U}^{n}R_{\mathcal{B}}^{{ - 1}}$, ${{p}_{\tau }} = \tfrac{1}{{{{\tau }^{2}}}} + \tfrac{{\operatorname{Re} {{\unicode{230} }^{2}}}}{\tau }$ и получим матричное уравнение:

${{p}_{\tau }}{{V}^{{n + 1}}} - {{\unicode{230} }^{2}}{{D}_{\mathcal{A}}}{{V}^{{n + 1}}} - {{\unicode{230} }^{2}}{{V}^{{n + 1}}}{{D}_{\mathcal{B}}} = {{\mathcal{G}}^{n}},$

решение которого относительно элементов ${v}_{{km}}^{{n + 1}}$ матрицы ${{V}^{{n + 1}}}$ дается элементарными формулами

${v}_{{km}}^{{n + 1}} = \frac{{g_{{km}}^{n}}}{{{{p}_{\tau }} - {{\unicode{230} }^{2}}(d_{\mathcal{A}}^{k} + d_{\mathcal{B}}^{m})}},\quad k = \overline {1,K} ,\quad m = \overline {1,M} ,$

где $g_{{km}}^{n}$ – элементы матрицы ${{\mathcal{G}}^{n}}$. В данном случае шаг сетки по времени $\tau $ нужно выбирать, исходя из условия ${{p}_{\tau }} \ne {{\unicode{230} }^{2}}(d_{\mathcal{A}}^{k} + d_{\mathcal{B}}^{m})$ $\forall k,m$ или учитывая, что ${{\unicode{230} }^{2}} = {{(W{\text{Re}})}^{{ - 1}}}$,

$\tau \ne \frac{{{\text{Re}} \pm \sqrt {{\text{R}}{{{\text{e}}}^{2}} + 4W\operatorname{Re} (d_{\mathcal{A}}^{k} + d_{\mathcal{B}}^{m})} }}{{2(d_{\mathcal{A}}^{k} + d_{\mathcal{B}}^{m})}}\quad \forall k,m.$

Зная элементы матрицы ${{V}^{{n + 1}}}$, несложно восстановить значения решения в узлах коллокации: ${{U}^{{n + 1}}} = {{R}_{\mathcal{A}}}{{V}^{{n + 1}}}{{R}_{\mathcal{B}}}$.

Отметим, что в силу нелинейности исходных уравнений для перехода от предыдущего шага по времени к следующему решение линеаризованных матричных уравнений (24), (25) ищется с помощью итераций по нелинейности. Иначе говоря, при фиксированном $t = {{t}_{n}}$ реализуется цикл, на каждом шаге $s$ которого решения $\Upsilon _{{ij}}^{{n,s}}$, ${{U}^{{n,s}}}$ выражаются через решения $\Upsilon _{{ij}}^{{n,s - 1}}$, ${{U}^{{n,s - 1}}}$ в соответствии с формулами, описанными выше, и начальными условиями $\Upsilon _{{ij}}^{{n,0}} = \Upsilon _{{ij}}^{n}$, ${{U}^{{n,0}}} = {{U}^{n}}$. При этом итерации останавливаются при реализации условия

$\frac{{\left\| {{{U}^{{n,s}}} - {{U}^{{n,s - 1}}}} \right\|}}{{\left\| {{{U}^{{n,s}}}} \right\|}} \leqslant {{\varepsilon }_{{NI}}},\quad \frac{{\left\| {\Upsilon _{{ij}}^{{n,s}} - \Upsilon _{{ij}}^{{n,s - 1}}} \right\|}}{{\left\| {\Upsilon _{{ij}}^{{n,s}}} \right\|}} \leqslant {{\varepsilon }_{{NI}}},\quad i = 1,2,3,\quad j = i,...,3.$

Номер $s$ в таком случае обозначим $\sigma $ и положим $\Upsilon _{{ij}}^{{n + 1}} = \Upsilon _{{ij}}^{{n,\sigma }}$, ${{U}^{{n + 1}}} = {{U}^{{n,\sigma }}}$. Здесь и далее норма обозначает максимальный элемент матрицы. В наших расчетах задано ${{\varepsilon }_{{NI}}} = {{10}^{{ - 10}}}$. При переходе на следующий шаг по времени меняются значения функций $A({{t}_{n}})$, $A{\kern 1pt} '({{t}_{n}})$, стоящих в правой части (25) и описывающих зависимость градиента давления в канале от времени.

Отметим, что предложенная аппроксимация уравнений модели в сравнении с классической схемой метода коллокаций приводит к существенному сокращению вычислительных затрат. Действительно, для вычисления $\Upsilon _{{ij}}^{{n + 1}}$ в соответствии с (24) необходимо совершить порядка $KM$ операций (деления элементов матрицы в правой части на элементы матрицы в левой части). Для расчета ${{U}^{{n + 1}}}$ по (25) наиболее затратными являются операции произведения матриц ${{\mathcal{G}}^{n}} = R_{\mathcal{A}}^{{ - 1}}\mathcal{F}_{U}^{n}R_{\mathcal{B}}^{{ - 1}}$, ${{U}^{{n + 1}}} = {{R}_{\mathcal{A}}}{{V}^{{n + 1}}}{{R}_{\mathcal{B}}}$. Порядок числа операций, необходимых для этого, есть $O(MK(M + K))$. Здесь предполагается, что спектральные разложения (23) сделаны заранее, и матрицы этих разложений для различных $K$, $M$ занесены в базу данных. Классическая же схема метода коллокаций в данной задаче ведет к системам с большими разреженными матрицами размера $KM \times KM$, решение которых требует порядка $O({{(KM)}^{3}})$ операций.

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

5. РЕЗУЛЬТАТЫ ЧИСЛЕННЫХ ЭКСПЕРИМЕНТОВ

Разработанный метод использован для решения нестационарной задачи (5)–(9). При этом градиент давления был задан постоянным $A(t) \equiv {{A}_{0}}$. Наряду с нестационарной методом установления решалась также стационарная задача (10), (12) при $A = {{A}_{0}}$. При этом использован алгоритм, предложенный в [7]. Отметим, что сходимость метода установления достигается только при использовании определенных ветвей решения $\sigma (\mu )$ уравнения (11). Какие именно ветви обеспечивают сходимость, и какая ветвь соответствует предельному решению нестационарной задачи (5)–(9), зависит от значений параметров модели.

В табл. 1 приведены описание и диапазоны значений основных параметров модели и алгоритма, использованные в расчетах. Обозначим время выхода течения на стационарный режим величиной $t{\kern 1pt} *$, т.е. $t{\kern 1pt} * = n{\kern 1pt} *{\kern 1pt} \tau $, где $n{\kern 1pt} *$ – номер итерации, на которой впервые реализуются условия стабилизации течения:

$\frac{{\left\| {{{U}^{{n*}}} - {{U}^{{n* - 1}}}} \right\|}}{\tau } \leqslant {{\varepsilon }_{S}},\quad \frac{{\left\| {\Upsilon _{{ij}}^{{n*}} - \Upsilon _{{ij}}^{{n* - 1}}} \right\|}}{\tau } \leqslant {{\varepsilon }_{S}},\quad i = 1,2,3,\quad j = 1,...,i,$

${{\varepsilon }_{S}}$ – малое число. В проведенных экспериментах установлено, что $t{\kern 1pt} *$ практически не зависит от параметров $K$, $M$ и $\tau $. Относительное отклонение нестационарной скорости течения в момент времени $t{\kern 1pt} *$ от решения стационарной задачи (10), (12) обозначим через ${{D}_{u}}(y,z)$. В экспериментах эта величина рассчитывалась как относительная разность двух массивов, содержащих значения соответствующих скоростей в узлах $({{y}_{k}},{{z}_{m}})$ (см. п. 4.2).

Таблица 1.

Параметры модели и численного метода

Параметр	Описание	Значение
$\beta $	Феноменологический параметр	0.01–0.9
${{c}_{k}}$	Феноменологический параметр, $k = {{c}_{k}}\beta $	0.5–20
${{A}_{0}}$	Безразмерный градиент давления	1–330
$W$	Число Вайсенберга	0.01–3
$2r$	Безразмерная ширина канала	0.02–2
$K$, $M$	Число узлов сетки вдоль осей $y$ и $z$	21–31
$\tau $	Шаг сетки по времени	0.0015–0.1
${{\varepsilon }_{S}}$	Погрешность стабилизации течения	10^–8

В численных экспериментах установлено, что в случае ${{c}_{k}} = 1$ и $\beta < 0.5$ (именно такие значения использованы в [2], [10]) сходимость метода установления в стационарной задаче имеет место только при $\sigma = {{\sigma }_{{1,3}}}$, причем предельное решение нестационарной задачи соответствует ветви $\sigma = {{\sigma }_{3}}$. В табл. 2 для этого случая приведены значения $t{\kern 1pt} *$ и $\left\| {{{D}_{u}}(y,z)} \right\|$ при вариации параметров $W$, ${{A}_{0}}$, $r$ в некоторой окрестности предельных значений этих параметров. Под предельными здесь понимаются значения, при переходе через которые стационарный режим течения перестает существовать. Эти значения найдены приближенно в расчетах и приведены в последних строках табл. 2, 3 с индексом lim. В табл. 3 указаны значения $t{\kern 1pt} *$ и $\left\| {{{D}_{u}}(y,z)} \right\|$ при вариации ${{c}_{k}}$ и $\beta $.

Таблица 2.

Значения $t{\kern 1pt} *$ и $\left\| {{{D}_{u}}(y,z)} \right\|$ при $\beta = 0.1$, ${{c}_{k}} = 1$ и вариации параметров $W$, ${{A}_{0}}$, $r$

Вариация $W$ при ${{A}_{0}} = 1$, $r = 1$			Вариация градиента ${{A}_{0}}$ при $W = 0.01$, $r = 1$			Вариация $r$ при ${{A}_{0}} = 200$, $W = 0.05$
$W$	$t{\kern 1pt} *$	$\left\\| {{{D}_{u}}(y,z)} \right\\|$	${{A}_{0}}$	$t{\kern 1pt} *$	$\left\\| {{{D}_{u}}(y,z)} \right\\|$	$r$	$t{\kern 1pt} *$	$\left\\| {{{D}_{u}}(y,z)} \right\\|$
1.5	41.75	$6.23 \times {{10}^{{ - 9}}}$	200	9.1	$1.32 \times {{10}^{{ - 8}}}$	0.025	1.446	$4.49 \times {{10}^{{ - 6}}}$
2.5	61.95	$2.62 \times {{10}^{{ - 8}}}$	300	14.88	$4.62 \times {{10}^{{ - 8}}}$	0.05	1.434	$1.1 \times {{10}^{{ - 5}}}$
3	108.65	$4.64 \times {{10}^{{ - 8}}}$	330	19.09	$5.31 \times {{10}^{{ - 8}}}$	0.1	1.584	$1.32 \times {{10}^{{ - 6}}}$
${{W}^{{{\text{lim}}}}} \approx 3.4$			$A_{0}^{{{\text{lim}}}} \approx 335$			${{r}^{{{\text{lim}}}}} \approx 0.115$

Таблица 3.

Значения $t{\kern 1pt} *$ и $\left\| {{{D}_{u}}(y,z)} \right\|$ при $W = 0.01$, $r = 1$, ${{A}_{0}} = 20$, 200 и вариации параметров $\beta $, ${{c}_{k}}$

Вариация $\beta $ при ${{A}_{0}} = 200$, ${{c}_{k}} = 1$			Вариация ${{c}_{k}}$ при ${{A}_{0}} = 20$, $\beta = 0.37$
$\beta $	$t{\kern 1pt} *$	$\left\\| {{{D}_{u}}(y,z)} \right\\|$	${{c}_{k}}$	$t{\kern 1pt} *$	$\left\\| {{{D}_{u}}(y,z)} \right\\|$
0.15	10.4	$2.58 \times {{10}^{{ - 8}}}$	1	5.1	$3.32 \times {{10}^{{ - 10}}}$
0.2	13.32	$4.39 \times {{10}^{{ - 8}}}$	1.52	5.1	$9.91 \times {{10}^{{ - 5}}}$
0.26	18.12	$6.69 \times {{10}^{{ - 8}}}$	2	5.15	$3.35 \times {{10}^{{ - 10}}}$
${{\beta }^{{{\text{lim}}}}} \approx 0.267$			неограничено

Видно, что при значениях параметров $W$ и ${{A}_{0}}$, близких к предельным, зависимость $t{\kern 1pt} *$ от этих значений носит нелинейный характер. Время выхода на стационарный режим наиболее чувствительно к значению $W$. При $W$, близких к ${{W}^{{{\text{lim}}}}}$, и малых $t$ имеют место колебания значений скорости, и для установления течения требуется достаточно много времени (см. также фиг. 5). При уменьшении ширины канала $r$ течение устанавливается намного быстрее и для бóльших значений $W$ и ${{A}_{0}}$, однако относительное отклонение $\left\| {{{D}_{u}}(y,z)} \right\|$ при этом возрастает (см. также фиг. 2б). При вариации r в окрестности малых значений ($r < 0.05$) время $t{\kern 1pt} *$ практически не зависит от $r$.

Фиг. 2.

Относительные отклонения ${{D}_{u}}(y,z)$ при (a) ${{A}_{0}} = 1$, $\beta = 0.1$, ${{c}_{k}} = 1$, $r = 1$, $W = 1.5$ (стандартная ситуация); (б) ${{A}_{0}} = 200$, $\beta = 0.1$, ${{c}_{k}} = 1$, $W = 0.05$, $r = 0.05$ (узкий канал); (в) ${{A}_{0}} = 20$, $W = 0.01$, $r = 1$, $\beta = 0.37$, ${{c}_{k}} = 1.52$ (переключение $3 \to 1$ установившегося решения нестационарной задачи).

Отметим, что зависимость $t{\kern 1pt} *$ от $\beta $ близка к линейной, и отклонение $\left\| {{{D}_{u}}(y,z)} \right\|$ при приближении $\beta $ к ${{\beta }^{{{\text{lim}}}}}$ мало меняется. Интересная картина наблюдается при вариации ${{c}_{k}}$ (табл. 3): отклонения $\left\| {{{D}_{u}}(y,z)} \right\|$ при изменении ${{c}_{k}}$ в широких диапазонах остаются малыми за исключением окрестности некоторой точки (в данном случае точки ${{c}_{k}} \approx 1.52$), где они резко возрастают (фиг. 2в). Такая особенность связана с переключением установившегося решения нестационарной задачи с решения, соответствующего ветви ${{\sigma }_{3}}$, на решение, соответствующее ветви ${{\sigma }_{1}}$, стационарной задачи. Этот важный эффект, который далее будем именовать “переключение $3 \to 1$”, обсудим чуть ниже.

На фиг. 2 для некоторых режимов течения показаны значения относительных отклонений ${{D}_{u}}(y,z)$ предельных решений нестационарной задачи от решений стационарной задачи. В большинстве случаев отклонения имеют высокочастотные колебания, амплитуда которых возрастает в окрестности границ, достигая значений порядка ${{10}^{{ - 10}}}$–${{10}^{{ - 8}}}$ (фиг. 2а). Проведенные эксперименты позволяют заключить, что такой характер поведения ${{D}_{u}}(y,z)$ связан с вычислительной погрешностью и значениями $K$, $M$ и ${{\varepsilon }_{S}}$. Однако в некоторых случаях наблюдаются конкретные тенденции в поведении ${{D}_{u}}(y,z)$: например, рост этих значений при уменьшении ширины канала в окрестности стенок, имеющих бóльшую длину (фиг. 2б); или рост отклонений в центре канала при переключении $3 \to 1$ (фиг. 2в).

На фиг. 3, 4 показаны зависимости скорости течения в центре канала от времени ${{u}_{c}}(t) = u(t,0,0)$ (штрихи), полученные при решении нестационарной задачи (5)–(9) с вариацией значений $W$ и ${{c}_{k}}$. На этих же графиках серыми линиями отмечены значения решений стационарных уравнений (10), (12) в точке $y = 0$, $z = 0$, полученные при использовании тех ветвей функции $\sigma (\mu )$, которые обеспечивают сходимость метода установления. Значения параметров рассмотренных режимов соответствуют указанным в табл. 2, 3 .

Фиг. 3.

Значения скорости ${{u}_{c}}(t)$, полученные при решении нестационарной задачи (штрихи), и соответствующие значения решений стационарной задачи в точке $y = 0$, $z = 0$ (серые линии) при ${{A}_{0}} = 1$, $\beta = 0.1$, ${{c}_{k}} = 1$, $r = 1$ и различных $W$: (а) $W = 1.5$, (б) $W = 2.5$, (в) $W = 3$.

Фиг. 4.

Значения скорости ${{u}_{c}}(t)$, полученные при решении нестационарной задачи (штрихи), и соответствующие значения решений стационарной задачи в точке $y = 0$, $z = 0$ (серые линии) при $W = 0.01$, $\beta = 0.37$, ${{A}_{0}} = 20$, $r = 1$ и различных ${{c}_{k}}$: (а) ${{c}_{k}} = 1$, (б) ${{c}_{k}} = 1.52$, (в) ${{c}_{k}} = 2$.

Фиг. 5.

Значения скорости $u(t,y,z)$, рассчитанные при ${{A}_{0}} = 1$, $\beta = 0.1$, ${{c}_{k}} = 1$, $r = 1$ и $W = 3$ в различные моменты времени: (а) $t = 2$, (б) $t = 5$, (в) $t = 100$.

Отметим, что колебания значений на фиг. 3 характеризуют сложный процесс установления течения при достаточно больших $W$. Более подробно этот процесс прослеживается на фиг. 5, где приведены распределения скорости $u(t,y,z)$ в различные моменты времени $t$. Видно, что на начальных этапах (при малых $t$) имеют место сильные колебания решения, а также течения в направлении, противоположном градиенту давления (области отрицательных значений скорости видны на фиг. 5a). При $t \approx 8$ сильные колебания затухают и течение плавно выходит на стационарный режим.

На фиг. 4 показано, что количество решений стационарной задачи (10), (12), которые удается получить с помощью метода установления, меняется в зависимости от параметров задачи. Существенное влияние при этом оказывают параметры $\beta $ и ${{c}_{k}}$. Причем при определенном соотношении этих параметров происходит переключение $3 \to 1$.

Более конкретно, при значениях параметров ${{A}_{0}} = 20$, $W = 0.01$, $r = 1$, $\beta = 0.37$ и ${{c}_{k}} < 1.515$ методом установления удается получить два стационарных решения, соответствующих ветвям ${{\sigma }_{1}}$, ${{\sigma }_{3}}$, при этом устойчивое решение соответствует ветви ${{\sigma }_{3}}$ (фиг. 4а). При ${{c}_{k}} \approx 1.515$ метод установления перестает сходиться для $\sigma = {{\sigma }_{3}}$, и происходит переключение $3 \to 1$. В диапазоне $1.515 < {{c}_{k}} < 1.528$ удается получить единственное решение стационарной задачи для $\sigma = {{\sigma }_{1}}$ (фиг. 4б). В диапазоне $1.528 < {{c}_{k}} < 1.532$ – два решения для $\sigma = {{\sigma }_{{1,2}}}$. При ${{c}_{k}} > 1.532$ – три решения для всех $\sigma = {{\sigma }_{{1,2,3}}}$ (фиг. 4в). При этом для ${{c}_{k}} > 1.515$ устойчивым остается решение, соответствующее ветви ${{\sigma }_{1}}$. Подчеркнем, что неравенства, указанные в этом примере, были получены в расчетах и являются приближенными. Заметим, что описанный характер исчезновения и возникновения численных решений стационарной задачи при увеличении ${{c}_{k}}$, а также переключения $3 \to 1$ сохраняется при вариации других параметров задачи.

Обозначим $c_{k}^{{3 \to 1}}$ значение ${{c}_{k}}$, при котором происходит переключение $3 \to 1$. На фиг. 6 показан график зависимости $c_{k}^{{3 \to 1}}$ от $\beta $ при различных значениях параметра $W$. Из фиг. 6 следует, что при малых значениях $W$ кривые $c_{k}^{{3 \to 1}}(\beta )$ практически совпадают (серая и штриховая линии), однако при увеличении $W$ кривая меняется. Аналогичная картина наблюдается при увеличении ${{A}_{0}}$.

Фиг. 6.

Кривые $c_{k}^{{3 \to 1}}(\beta )$ при ${{A}_{0}} = 1$, $r = 1$ и различных $W$: $W = 0.01$ (штрихи), $W = 0.1$ (серая линия), $W = 1$ (черная линия).

Обозначим символом $c_{k}^{{1,2,3}}$ минимальное значение параметра ${{c}_{k}}$, при котором стационарная задача имеет три решения, соответствующие ветвям ${{\sigma }_{{1,2,3}}}$, и рассмотрим зависимость $c_{k}^{{1,2,3}}(\beta )$. Из примера, приведенного выше, можно заключить, что $c_{k}^{{1,2,3}}(0.37) \approx 1.532$. При проведении вычислений установлено, что линия $c_{k}^{{1,2,3}}(\beta )$ практически не меняется при вариации других параметров задачи и с высокой точностью аппроксимируется гиперболой $c_{k}^{{1,2,3}}(\beta ) \approx 3{\text{/}}(4\beta ) - 0.5$. График этой гиперболы визуально неотличим от штриховой линии на фиг. 6. При уменьшении параметров $W$ и ${{A}_{0}}$ кривые перехода $c_{k}^{{3 \to 1}}(\beta )$ стремятся к этой гиперболе снизу. При приближении значений $W$ и ${{A}_{0}}$ к предельным, кривая $c_{k}^{{3 \to 1}}(\beta )$ существенно отдаляется от $c_{k}^{{1,2,3}}(\beta )$ и меняет форму. При малых значениях $W$ и ${{A}_{0}}$ и вариации $r$ в широких диапазонах кривые $c_{k}^{{3 \to 1}}(\beta )$ остаются в окрестности $c_{k}^{{1,2,3}}(\beta )$.

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

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

Лойцянский Л.Г. Механика жидкости и газа. М.: Наука, 1978 [Loitsyanskii L.G. Mechanics of liquids and gases. Oxford: Stewartson Pergamon Press; Moscow: Nauka, 1978].
Алтухов Ю.А., Гусев А.С., Пышнограй Г.В. Введение в мезоскопическую теорию текучести полимерных систем. Барнаул: АлтГПА, 2012 [Altukhov Yu.A., Gusev A.S., Pyshnograi G.V. Introduction to the mesoscopic theory of flowing polymer systems. Barnaul: AltGPA, 2012].
Pokrovskii V.N. The mesoscopic theory of polymer dynamics. 2nd ed. Berlin: Springer, 2010.
Vinogradov G.V., Pokrovskii V.N., Yanovsky Yu.G. Theory of viscoelastic behavior of concentrated polymer solutions and melts in one-molecular approximation and its experimental verification // Rheol. Acta. 1972. V. 7. P. 258–274.
Блохин А.М., Семисалов Б.В. Стационарное течение несжимаемой вязкоупругой полимерной жидкости в канале с эллиптическим сечением // Сиб. журн. индустриал. матем. 2014. Т. XVII. № 4 (60). С. 38–47 [Blokhin A.M., Semisalov B.V. A stationary flow of an incompressible viscoelastic fluid in a channel with elliptic cross section // J. Appl. Industr. Math. 2015. V. 9. 1. P. 18–26].
Блохин А.М., Семисалов Б.В. Расчет стационарных неизотермических МГД течений полимерной жидкости в каналах с внутренними нагревательными элементами // Сиб. журн. индустриал. матем. 2020. Т. 23. № 2. С. 17–40 [Blokhin A.M., Semisalov B.V. Simulation of the stationary nonisothermal MHD flows of polymeric fluids in channels with interior heating elements // J. Appl. Industr. Math. 2020. V. 14. № 2. P. 222–241].
Семисалов Б.В. Быстрый нелокальный алгоритм решения краевых задач Неймана–Дирихле с контролем погрешности // Вычисл. методы и программирование. 2016. Т. 17. № 4. С. 500–522 [Semisalov B.V. A fast nonlocal algorithm for solving Neumann-Dirichlet boundary value problems with error control // Vych. Metody Programm. 2016. V. 17. Is. 4. P. 500–522].
Седов Л.И. Механика сплошной среды. Т. 1. М.: Наука, 1970 [Sedov L.I. Mechanics of continuous media. V. 1. Singapore: World Sci., 1997].
Pai Shih-I. Introduction to the theory of compressible flow. Literary Licensing, LLC. 2013.
Зинович С.А., Головичёва И.Э., Пышнограй Г.В. Влияние молекулярной массы на сдвиговую и продольную вязкость линейных полимеров // Прикладн. механ. и техн. физ. 2000. Т. 41. № 2. С. 154–160 [Golovicheva S.A., Zinovich G.V., Pyshnograi G.V. Effect of the molecular mass on the shear and longitudinal viscosity of linear polymers // J. Appl. Mech. Tech. Phys. 2000. V. 41. P. 347–352].
Бабенко К.И. Основы численного анализа. М.–Ижевск: НИЦ “Регулярная и хаотическая динамика”, 2002 [Babenko K.I. Fundamentals of numerical analysis. Moscow–Izhevsk: Sci. Publ. Center “Regular and chaotic dynamics”, 2002].

Дополнительные материалы отсутствуют.

Инструменты

следующая статья выпуска предыдущая статья выпуска содержание выпуска

Журнал вычислительной математики и математической физики

Архивы выпусков Информация о журнале Отправить рукопись в журнал