Численный анализ разрушения одномерного течения полимерной жидкости с фронтом

Texto integral

Введение

Экспериментальные и теоретические исследования течений растворов и расплавов полимерных материалов со специфическими реологическими свойствами представляют важное направление современной гидродинамики. Их актуальность обусловлена бурным развитием технологий производства электронных устройств на полимерной основе методами печати, экструзии, напыления и др. (см. [1–3]). Обеспечение бездефектного производства в подобных технологиях зачастую связано с реализацией устойчивых ламинарных течений в каналах. В связи с этим вопрос об описании и контроле ламинарно-турбулентных переходов в таких течениях является принципиально важным (см. [4]).

Изучение этого вопроса наталкивается на ряд сложностей и новых эффектов. Во-первых, обнаружены течения вязкоупругих жидкостей, в которых вязкость практически не влияет на переход к турбулентности, а важными являются упругие свойства (см. критерий из [5]). Возникает задача описания непрерывного перехода между режимами такой (упругой) и классической (инерционной) турбулентностей (см. [6; 7]). Во-вторых, до сих пор не существует общепринятой математической модели разрушения ламинарных течений полимерной жидкости с прямыми линиями тока (течения Пуазейля и Куэтта), которое наблюдается в экспериментах при существенно более низких значениях числа Рейнольдса, Re, чем в случае течений ньютоновской жидкости (см. [8; 9]). Исследование устойчивости течения, которое, как правило, удается провести только на линейном уровне, дает качественное объяснение перехода к упруго-инерционной турбулентности в разбавленных растворах полимеров (см. [10; 11]), однако является недостаточным для случая чисто упругой турбулентности при малых Re в концентрированных растворах и расплавах полимеров. Здесь, по всей видимости, принципиальным является вклад нелинейных эффектов (см. [4]).

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

Для описания течений высококонцентрированных растворов и расплавов полимерной жидкости будем использовать реологическую мезоскопическую модель Покровского–Виноградова (см. [12; 13]). Ее важное отличие от трендовых моделей, таких как степенные модели, Oldroyd-B, FENE-P и т.п. (см. обзор [4, разд. III]), состоит в учете микроструктуры (конкретно, размера и ориентации полимерных молекул) при макроскопическом моделировании течения, что может оказаться принципиально важным для анализа разрушения устойчивых ламинарных течений. В настоящей работе мы приводим упрощенный одномерный (1D) вариант мезоскопической модели, аналогичный по своей сути уравнению Бюргерса для системы Навье–Стокса. Отметим, что именно на примере упрощенной 1D модели, предложенной Бюргерсом, были даны первые описания нелинейных волновых и турбулентных процессов в течениях классических жидкостей (см. [14]). В последствии эта модель стала эталонным примером для разработки и тестирования численных методов. Конечно-разностные, проекционные и спектральные методы, полулагранжевы и адаптивные схемы тестировались на примере решения уравнения Бюргерса (см. краткий обзор во введении [15]). Подчеркнем, что спектральные подходы в рамках нашего исследования имеют существенные преимущества над другими, во-первых, потому что в случае принадлежности искомого решения классу аналитических функций они реализуют экспоненциальную асимптотику сходимости, а, во-вторых, они являются чувствительными к особенностям решения: скорость их сходимости строго связана с положением и типом особых точек аналитического продолжения решения в комплексную плоскость.

Мезоскопическая 1D-модель существенно сложнее уравнения Бюргерса, она включает систему из двух квазилинейных уравнений первого порядка с параметрами, характеризующими вязкость жидкости (число Рейнольдса, Re), упругость полимерных молекул (число Вайсенберга, W) и анизотропию течения, b. Численный анализ этой модели для течений типа Пуазейля в канале с сечением, ограниченным двумя софокусными эллипсами (см. [16]), показал, что у решения имеются особенности, влияние которых усиливается с ростом Re, W и b, что в итоге приводит к расходимости нескольких вычислительных алгоритмов, построенных в указанной работе. Однако отдельно вопрос о траекториях движения особых точек в той работе не исследовался.

Отметим, что для ряда других моделей этот вопрос активно обсуждается в работах зарубежных исследователей (см. [17–20]). Методы локализации особенностей, предложенные в этих работах, опираются на результаты комплексного анализа (теорему Коши, разложение Лорана, принцип аргумента и др.), теории приближений (теоремы о сходимости рядов Фурье для ортогональных систем функций) и недавний прогресс в области дробно-рациональных приближений. По поводу последнего отметим работы Г. Шталя, Е. А. Рахманова, С. П. Суетина, Л. Н. Трефетена и их соавторов (см. [21–26]), в которых обсуждаются асимптотики сходимости наилучших дробно-рациональных приближений, приближений Чебышёва–Паде и Эрмита–Паде; а также экспоненциальное сгущение и чередование нулей и полюсов таких аппроксимаций в окрестности особых точек приближаемой функции. Заметим однако, что применение этих результатов для вычисления координат особых точек решений дифференциальных уравнений наталкивается на существенную сложность: алгоритмы поиска нулей и полюсов соответствующих дробно-рациональных приближений при высокой степени числителя и знаменателя, как правило, работают неустойчиво (см. [18, заключение; 26, разд. 7] и заключение данной статьи).

Эффективный способ решения этой проблемы при наличии одной особенности, состоящий в применении аппроксимаций Чебышёва–Паде, где в числителе стоит полином высокой степени, а в знаменателе — полином второй степени, предложен в [27] и развит в [28]. В этих работах особые точки решений начально-краевых задач приближались полюсами аппроксимаций Чебышёва–Паде, а для аппроксимации решения использовались барицентрические интерполяционные формулы. Преимущество последних над классическими интерполяциями с узлами в нулях ортогональных полиномов состоит в возможности адаптации узлов к положениям особенностей приближаемой функции с помощью специальных замен переменной (см. по поводу таких замен [29–31], по поводу адаптивных алгоритмов — [32]).

В настоящей работе предложено обобщение алгоритма из [27] на случай начально-краевой задачи для системы 1D квазилинейных уравнений мезоскопической модели, проверена его сходимость при различных значениях Re и найдены соотношения между параметрами Re и W, обеспечивающие переход к режиму с разрушением гладкого решения задачи. Под разрушением здесь и далее мы понимаем расходимость используемых приближений, связанную с тем, что точка аналитического продолжения решения в $ℂ.$ достигает области задачи. В работе также исследованы зависимости времени существования решения и координаты точки, в которой происходит разрушение, от значений Re и W.

1. Одномерная модель течения полимерной жидкости и постановка задачи

В настоящей работе предложена 1D-модель течения несжимаемой вязкоупругой полимерной жидкости, учитывающая эффекты релаксации и размеры полимерных молекул. Она основана на мезоскопическом подходе к описанию динамики полимерных сред (см. [12; 13]). Для анализа разрушения течений в рамках мезоскопического подхода запишем безразмерные уравнения из [33] в 1D-случае:

$\begin{array}{l}\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}=\frac{1}{\mathrm{Re}}\frac{\partial a}{\partial x},\\ \frac{\partial a}{\partial t}+u\frac{\partial a}{\partial x}=2\left(a+W^{-1}\right)\frac{\partial u}{\partial x}-W^{-1}a-\beta a^2,\end{array}$(1)

где $u\left(t,x\right)$ — скорость, $a\left(t,x\right)$ характеризует напряжение в жидкости; t — время, x — пространственная координата; $\mathrm{Re}=\frac{\rho u_{H}l}{{\eta }_0^{*}}$ — число Рейнольдса, $\rho \left(=const\right)$ — плотность среды, u_H — характерная скорость, l — характерная длина; $W=\frac{{\tau }_0^{*}{u}_H}{l}$ — число Вайсенберга; ${\eta }_0^{*}$, ${\tau }_0^{*}$ — начальные значения сдвиговой вязкости и времени релаксации (см. [12; 13]); b характеризует размер и ориентацию молекул полимера (см. [13]).

Система (1) записана в безразмерном виде: время t, координата x, скорость u и характеристика напряжения a отнесены к $T_C=l/u_H$, l, u_H и $\mathrm{Re}W/3$ соответственно, где T_C — характерный временной масштаб модели.

Отметим, что в случае $\beta =0$ систему (1) можно вывести из модели Oldroyd-B при нулевой вязкости растворителя, а в случае $W\to 0$ из второго уравнения (1) имеем $a\to \partial u/\partial x,$ и вся система сводится к уравнению Бюргерса

$\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}=\frac{1}{\mathrm{Re}}\frac{{\partial }^{2}u}{\partial x^2}.$ (2)

Заметим, что в модели из [33] функция $a\left(t,x\right)$ является компонентой тензора анизотропии среды, поэтому система (1) учитывает лишь 1D-проекцию сложной пространственной картины движения полимерных молекул, что, конечно, является значительным упрощением. Однако мы будем полагать, что такое 1D-упрощение является качественно верным. Подтверждением этому служат предельные случаи, указанные выше.

Для (1) будем рассматривать начально-краевую задачу с условиями

$u\left(\mathrm{0,}x\right)=\exp \left(-8x^2\right),a\left(\mathrm{0,}x\right)\equiv \mathrm{0,}x\in [-\mathrm{1,1}];u\left(t,\pm 1\right)=\mathrm{0,}t>0.$ (3)

Отметим, что аналогичная постановка для уравнения (2) рассматривалась в [28]. Для описания процесса формирования больших градиентов и разрывов решения (1), (3) используем подход из [27; 28], где показано, что возникновение больших градиентов и разрушение гладких решений уравнения Бюргерса при предельно малых 1 / Re можно связать с наличием у решения особых точек в комплексной плоскости $ℂ.$ Эти точки (полюса, точки ветвления и т.п.) могут возникать в определенные моменты времени и двигаться по сложным траекториям в окрестности отрезка действительной оси, где ищется решение. Выход точки на отрезок ведет к разрыву решения или его производных.

Наш основной интерес состоит в изучении связей между реологическими параметрами жидкости Re, W, b и траекториями движения особых точек решения задачи (1), (3) в $ℂ.$ Подчеркнем, что именно характер движения особых точек, если они существуют, определяет формирование фронтов, сильных и слабых разрывов и возникновение решений типа blow-up.

2. Численный алгоритм

2.1. Аппроксимация решения по пространственной переменной

Для приближения неизвестных функций по переменной x используем барицентрическую форму записи интерполяционного многочлена с узлами Чебышёва–Лобатто $x_k,k=\mathrm{0,}\dots ,N$ (см. [34]). Например, для функции u(t,x) это приближение имеет вид

$u\left(t,x\right)\approx P_N\left[u\right]\left(t,x\right)=\frac{{\displaystyle \sum _{k=0}^N\frac{{\lambda }_k}{x-x_k}u\left(t,x_k\right)}}{{\displaystyle \sum _{k=0}^N\frac{{\lambda }_k}{x-x_k}}},x_k=\cos \frac{k\pi }{N},$ (4)

${\lambda }_k=\frac{1}{{{\Theta }^{\text{'}}}_N^{}\left(x_k\right)}=\frac{2^{N-1}}{N}{\left(-1\right)}^k{\delta }_k,{\Theta }_N\left(x\right)={\prod }_{j=0}^N\left(x-x_j\right),{\delta }_k=\left\{\begin{array}{c}\frac{1}{2},k=\mathrm{0,}N,\\ 1èíà÷å.\end{array}\right.$

Числа l_k называются весами интерполяционной формулы. Приближения (4) обладают свойством прямой численной устойчивости, которое заключается в том, что разница интерполяционных полиномов, вычисленных в точной арифметике и при наличии ошибок округления, ограничена величиной, зависящей от приближаемой функции f и медленно растущей с ростом N (см. оценки констант Лебега из [35]). Кроме того, расчет значения $P_N\left[u\right]\left(t,x\right)$ в любой точке $x\ne x_k$ требует порядка O(N) арифметических операций, что делает формулу (4) экономичной.

Положим, что функция $u\left(t,x\right)$ аналитична и ее продолжение в $ℂ$ (z — координата комплексной плоскости) имеет особую точку $z^{*}=\delta +i\epsilon$, расположенную на малом расстоянии от отрезка $[-\mathrm{1,1}]$. Для адаптации (4) к этой особенности будем использовать следующую модификацию барицентрического приближения (см. [29]):

$\begin{array}{c}u\left(t,x\right)\approx R_N\left[u\right]\left(t,x\right)=\frac{{\displaystyle \sum _{k=0}^N\frac{{\lambda }_k}{x-{\xi }_k}u\left(t,{\xi }_k\right)}}{{\displaystyle \sum _{k=0}^N\frac{{\lambda }_k}{x-{\xi }_k}}},\\ {\xi }_k=g\left(x_k\right)=g\left(x_k,\delta ,\epsilon \right).\end{array}$ (5)

Отметим, что в отличие от $P_N\left[u\right]\left(t,x\right)$ приближение $R_N\left[u\right]\left(t,x\right)$ является дробно-рациональным. Здесь $g:D_1\to D_2$ — конформное отображение, $D_1,D_2$ — области в , содержащие отрезок $[-\mathrm{1,1}]$, $g\left(x\right)$ переводит $[-\mathrm{1,1}]$ в себя, а обратное отображение $g^{-1}\left(\xi \right)$ переводит $z^{*}$ в точку, лежащую значительно дальше от отрезка $[-\mathrm{1,1}]$. Это приводит к увеличению размера эллипса Бернштейна, соответствующего исходному полиномиальному приближению, и, как следствие, к ускорению сходимости. Оценка функционала погрешности приближения (5) получена в [29]. Напомним, что в теории приближений эллипс Бернштейна аналитической функции f, определенной на отрезке $[-\mathrm{1,1}]$, — эллипс с фокусами в точках $\pm 1$, в который f может быть аналитически продолжена. Этот эллипс содержит на своем контуре ближайшую к отрезку $[-\mathrm{1,1}]$ особую точку f. При этом скорость сходимости полиномиальных приближений функции f оценивается как $O\left({\tilde{\rho }}^{-N}\right)$, где $\tilde{\rho }$ — сумма полуосей эллипса Бернштейна, N — степень приближающего полинома.

Далее будем использовать отображение

$$\begin{gathered} g\left(x,\delta ,\epsilon \right)=\delta +\epsilon \sinh \left\{\left(a^{-}+a^{+}\right)\frac{x-1}{2}+a^{-}\right\}, \\ {a}^{\pm }={\sinh }^{-1}\frac{1\pm \delta }{\epsilon }. \end{gathered}$$(6)

Детали, касающиеся построения функции такого вида, описаны в [27].

Покажем сначала преимущества приближения (5) на сетке x_k, k = 0, ... ,N, построенной с использованием (6), на примере функции, предложенной К. Рунге,

$f\left(x\right)=\frac{1}{1+25x^2},x\in \left[-\mathrm{1,1}\right].$

Помимо дробно-рационального приближения (5) будем использовать интерполяционные полиномы на сетке с узлами в нулях многочленов Чебышёва (4) (при равномерном расположении узлов, как известно, наблюдается расходимость) и кубические сплайны с заданными на концах точными значениями вторых производных приближаемой функции.

На фиг. 1а величина N для (4) и (5) обозначает число узлов, а для кубических сплайнов — количество подынтервалов разбиения. Видно, что для получения абсолютной погрешности $\parallel E_a{\parallel }_{\infty }$ порядка 10^–15 для дробно-рациональных приближений требуется около 60 узлов, а для алгебраических полиномов — 170. Зависимость логарифма погрешности интерполяции кубическими сплайнами от N является выпуклой вниз функцией, что говорит о сходимости с конечным порядком. Действительно, порядок сходимости здесь равен четырем независимо от N. В этом случае для получения точности порядка 10^–15 потребуется уже около 16 500 подынтервалов.

 

Фиг. 1. (а) Десятичный логарифм погрешности для дробно-рационального приближения на адаптивной сетке (1), для интерполяционного полинома с чебышёвскими узлами (2), для интерполяции кубическими сплайнами на равномерной сетке (3); (б) эллипсы Бернштейна для приближаемой функции f (штрихпунктир) и $f\circ g$ (пунктир).

 

Особенность функции Рунге находится в точках $A_{\mathrm{1,2}}=\pm i/5$. Обратное к (6) отображение переводит их приблизительно в $B_{\mathrm{1,2}}=\pm 0.68i$. На фиг. 1б приведены эллипсы Бернштейна функций f и $f\circ g$. Для f сумма полуосей эллипса Бернштейна ${\tilde{\rho }}_f\approx 1.22$, для $f\circ g$ — ${\tilde{\rho }}_{f\circ g}\approx 1.89$, при этом тангенс угла наклона линии десятичного логарифма погрешности на фиг. 1а для интерполяционного полинома с чебышёвскими узлами (линия 2) с хорошей точностью равен $\lg {\tilde{\rho }}_f$, а для дробно-рационального приближения (линия 1) — $\lg {\tilde{\rho }}_{f\circ g}$, что согласуется с оценками погрешности из [29].

По сравнению с особыми точками решения начально-краевой задачи (1), (3) особенность функции Рунге лежит далеко от действительной оси. Тем не менее даже в этом примере дробно-рациональные интерполяционные формулы (5) со специальной адаптацией к особенности функции (6) существенно выигрывают у других методов приближения. Отметим, что аналогичные преимущества такой аппроксимации наблюдаются при решении начально-краевых задач для уравнения, описывающего тепловой взрыв, и уравнения Бюргерса в [27], а также задачи Дирихле для нелинейного эллиптического уравнения с особенностью в комплексной плоскости в двумерном случае (см. [32]). Далее будем применять такие приближения для решения задачи.

Эффективность приближений (5) при решении начально-краевых задач связана также с возможностью применения матричных аппроксимаций операторов дифференцирования. Пусть x — вектор-столбец длины $N+1$, состоящий из узлов интерполяции, $u\left(x\right)$ и $u^{\text{'}}\left(x\right)$ — вектор-столбцы значений функции $u\left(t_0,x\right)$ и ее производной по x в узлах интерполяции в некоторый момент времени t₀, $d{R}_N\left[u\right]/dx$ — вектор значений производной $R_N\left[u\right]\left(t_0,x\right)$ по x. Для матричной аппроксимации операции дифференцирования в работе [36] на основе (4) выведена следующая формула:

$u^{\text{'}}\left(x\right)\approx \frac{d{R}_N\left[u\right]}{dx}=D^{1}u\left(x\right),D_{ij}^1=\left(\begin{array}{c}\frac{{\lambda }_j}{{\lambda }_i\left(x_i-x_j\right)},åñëèi\ne j,\\ -\sum _{k=\mathrm{0,}k\ne j}^N{D}_{ik}^1,åñëèi=j,\end{array}\right.$

которая справедлива также для случая приближения (5) при замене x_j на x_j. Здесь $D^1=\left(D_{ij}^1\right)$ — матрица, аппроксимирующая оператор производной первого порядка.

2.2. Дискретизация задачи по времени

Дискретизация задачи (1), (3) по времени описана на примере первого уравнения системы (1):

$\frac{\partial u}{\partial t}=f\left(t,x,u,\frac{\partial u}{\partial x},\frac{\partial a}{\partial x}\right)=\frac{1}{\mathrm{Re}}\frac{\partial a}{\partial x}-u\frac{\partial u}{\partial x}.$

Будем искать значения решения $u\left(t,x\right)$ в моменты времени $t=t^n$ в узлах Чебышёва–Лобатто, $x=x_l$, $l=\mathrm{0,}\dots ,N$. Пусть $u^{n+1},u^n,$ $a^{n+1},a^n$ — векторы значений $u\left(t,x\right)$ и $a\left(t,x\right)$ в точках $x_l$ на $\left(n+1\right)$-м и $n$-м слоях по времени соответственно, $x$ — вектор, состоящий из точек $x_l$, $f\left(t^n,x,u^n,D^1{u}^n,D^1{a}^n\right)$ — вектор значений $f\left(t^n,x_l,u_l^n,{\sum }_{m=0}^N{D}_{lm}^1{u}_m^n,{\sum }_{m=0}^N{D}_{lm}^1{a}_m^n\right)$, $u_l^n=u\left(t^n,x_l\right)$, $a_l^n=a\left(t^n,x_l\right)$. Для перехода со слоя $n$ на слой $n+1$ применим формулы метода Рунге–Кутты:

$\begin{array}{c}{u}^{n+1}=u^n+{\tau }_n\sum _{j=1}^s{b}_j{k}_j^n,\\ k_1^n=f\left(t^n,x,u^n,D^1{u}^n,D^1{a}^n\right),\\ k_j^n=f\left(t^n+c_j{\tau }_n,x,u^n+\sum _{l=1}^{j-1}{\alpha }_{jl}{k}_l^{n-1},D^1\left(u^n+\sum _{l=1}^{j-1}{\alpha }_{jl}{k}_l^{n-1}\right),D^1\left(a^n+\sum _{l=1}^{j-1}{\alpha }_{jl}{k}_l^{n-1}\right)\right).\end{array}$ (7)

В последнем выражении $j=\mathrm{2,}\dots ,s$; ${\tau }_n$ — шаг по времени, заданный на n-м слое по времени, s — количество стадий метода, которое в расчетах задано равным шести. Значения ${\alpha }_{jl},c_j,b_j$ для схем 4-го и 5-го порядков можно найти в [37].

2.3. Локализация особой точки решения

Приближение (5) использует отображение (6), которое включает зависимость от координат особой точки приближаемой функции. Следуя идее из [27], будем полагать, что решение задачи (1), (3) также имеет особую точку с неизвестными априори координатами, которые необходимо вычислить. Для аппроксимации координат особенности используем на каждом временном шаге приближение Чебышёва–Паде неизвестной функции и вычислим нули его знаменателя. Это приближение имеет вид (см. [38, ч. 2, п. 1.6])

${\tilde{u}}^n\approx r_{M,S}\left(x\right)=\frac{p\left(x\right)}{q\left(x\right)}=\frac{{\displaystyle \sum _{k=0}^M{b}_k{T}_k\left(x\right)}}{1+{\displaystyle \sum _{k=1}^S}{c}_k{T}_k\left(x\right)},$ (8)

где $T_k\left(x\right)$ — полином Чебышёва I рода степени k, $b_k,c_k$ — коэффициенты, подлежащие определению, ${\tilde{u}}^n$ — приближенное решение, полученное на n-м шаге по времени на сетке x_k, $k=\mathrm{0,}\dots ,N$ (см. (5)). Следуя [27], здесь в качестве M мы выбирали $N/2$, округленное до целого числа, и задавали $S=2$. Далее использовали разложение Фурье–Чебышёва функции ${\tilde{u}}^n\left(x\right)$:

${\tilde{u}}^n\left(x\right)=\sum _{k=0}^{\infty }{a}_k{T}_k\left(x\right).$ (9)

Здесь коэффициенты a_k можно вычислить, зная значения ${\tilde{u}}^n\left(x\right)$ в узлах Чебышёва–Лобатто $x={\tilde{x}}^j=\cos \left(j\pi /\nu \right)$, $j=\mathrm{0,}\dots ,\nu$, $\nu =M+S$, с помощью быстрого преобразования Фурье. Умножая (8) на функцию $q\left(x\right)$ и собирая подобные члены при $T_k\left(x\right)$, получаем

${\tilde{u}}^n\left(x\right)q\left(x\right)-p\left(x\right)=O\left(T_{\nu +1}\right).$ (10)

Подставляя сюда разложения Фурье–Чебышёва функций ${\tilde{u}}^n\left(x\right)$, $q\left(x\right)$, $p\left(x\right)$, приходим к системе линейных алгебраических уравнений (СЛАУ) для определения коэффициентов $b_k,c_k$:

$\left(\begin{array}{cccccc}1& & & -d_{\mathrm{0,1}}& \dots & -d_{\mathrm{0,}S}\\ & \ddots & & ⋮& ⋮& ⋮\\ & & 1& -d_{M\mathrm{,1}}& \dots & -d_{M,S}\\ & & & -d_{M+\mathrm{1,1}}& \dots & -d_{M+\mathrm{1,}S}\\ & & & ⋮& ⋮& ⋮\\ & & & -d_{M+S\mathrm{,1}}& \dots & -d_{M+n,S}\\ & & & & & \end{array}\right)\left(\begin{array}{c}{b}_0\\ ⋮\\ b_M\\ c_1\\ ⋮\\ c_S\end{array}\right)=\left(\begin{array}{c}{a}_0\\ ⋮\\ a_M\\ a_{M+1}\\ ⋮\\ a_{M+S}\end{array}\right).$

Здесь элементы матрицы $d_{j,k}$ находятся из условия ${\sum }_{j=0}^{\infty }{d}_{j,k}{T}_j\left(x\right)={\sum }_{j=0}^{\infty }{a}_j{T}_j\left(x\right)T_k\left(x\right)$. Для приближенного вычисления координат особой точки необходимо найти только коэффициенты $c_k,k=\mathrm{1,}\dots ,S$, так как ее положение определяется нулями знаменателя (8). В случае $S=2$ элементы матрицы СЛАУ для вычисления c_k можно выразить явно:

$d_{M+\mathrm{1,1}}=a_{M+1}/2+a_{M+3}/\mathrm{2,}{d}_{M+\mathrm{1,2}}=a_M/2+a_{M+4}/\mathrm{2,}$

$d_{M+\mathrm{2,1}}=a_{M+2}/2+a_{M+4}/\mathrm{2,}{d}_{M+\mathrm{2,2}}=a_{M+1}/2+a_{M+5}/2.$

После решения системы необходимо отыскать корни знаменателя (8). Они являются собственными числами сопровождающей матрицы (см. [39, гл. 18]). Однако при $S=2$ несложно записать знаменатель в виде полинома второй степени и найти точные значения его корней.

Отметим, что для более точной локализации особой точки в алгоритме используется информация о ее положении на предыдущем шаге по времени (см. [27; 28]). Пусть $\hat{\delta },\hat{\epsilon }$ — координаты особой точки на предыдущем шаге, тогда разложение (9) на текущем шаге делается не на отрезке $[-\mathrm{1,1}]$, а на суженном отрезке $[\max \left(\widehat{\delta }-10\widehat{\epsilon },-1\right),\min \left(\widehat{\delta }+10\widehat{\epsilon }\mathrm{,1}\right)]$.

2.4. Схема численного решения задачи (1), (3)

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

Шаг 1. Задаем сетку из узлов Чебышёва–Лобатто x_k, $k=\mathrm{0,}\dots ,N$, задаем начальный шаг по времени t и находим значения $u\left(\mathrm{0,}{x}_k\right)$, $a\left(\mathrm{0,}{x}_k\right)$.

Шаг 2. Делаем переход на следующий шаг по времени по формулам Рунге–Кутты вида (7).

Шаг 3. Находим приближенно координаты (d,e) особой точки функции $u\left(t^n,x\right)$ на текущем шаге с помощью аппроксимации Чебышёва–Паде.

Шаг 4. Адаптируем сетку к положению особенности, используя отображение $g\left(x_k,\delta ,\epsilon \right)$. Переинтерполируем текущее решение на новую сетку, используя (5).

Шаг 5. Повторяем шаги 2–4, пока не дойдем до заданного момента времени $t=t_{\max }.$

Отметим, что в ряде случаев (об этом сказано ниже) на шаге 3 вместо особых точек $u\left(t^n,x\right)$ отыскивались особые точки функции $a\left(t^n,x\right)$. В качестве $t_{\max }$ задавали время разрушения гладкого решения, которое находилось экспериментально. Для случаев, когда решение не разрушается, значения $t_{\max }$ определены для каждого эксперимента отдельно.

3. Результаты расчетов

3.1. Верификация предложенного алгоритма

Для тестирования алгоритма задача (1), (3) решалась при $W={10}^{-4}$, $\beta =0$, $\mathrm{Re}=100.$ Расчет велся до момента времени $t_{\max }=0.6$. Для анализа сходимости метода рассчитаны относительные отклонения $\parallel E_r^u{\parallel }_{\infty },$ $\parallel E_r^a{\parallel }_{\infty }$ решений задачи (1), (3), полученные при увеличении числа узлов пространственной сетки:

${||E_r^u||}_{\infty }=\frac{\underset{x\in Q}{{\displaystyle \max }}\left|{\tilde{u}}_{N_1}\left(x\right)-{\tilde{u}}_{N_2}\left(x\right)\right|}{\underset{x\in Q}{{\displaystyle \max }}\left|{\tilde{u}}_{N_2}\left(x\right)\right|},{||E_r^a||}_{\infty }=\frac{\underset{x\in Q}{{\displaystyle \max }}\left|{\tilde{a}}_{N_1}\left(x\right)-{\tilde{a}}_{N_2}\left(x\right)\right|}{\underset{x\in Q}{{\displaystyle \max }}\left|{\tilde{a}}_{N_2}\left(x\right)\right|},$ (11)

где ${\tilde{u}}_{N_1},{\tilde{u}}_{N_2}$ — приближенные решения для скорости при $t=t_{\max },$ полученные на сетках из $N_1+\mathrm{1,}{N}_2+1$ узлов соответственно (аналогично для ${\tilde{a}}_{N_1},{\tilde{a}}_{N_2}$), $N_2>N_1$, $Q$ — множество из 1000 точек, равномерно распределенных на отрезке $[-\mathrm{1,1}]$.

В тестовом расчете показано, что относительные отклонения решений при $t_{\max }=0.6$ и $N_2=N_1+10$ с ростом N₁ стремительно уменьшаются и при $N_1=50$ составляют примерно 10^–8. Этот расчет также показал, что алгоритм вычисления координаты особой точки устойчив, только если точка лежит достаточно близко к отрезку $[-\mathrm{1,1}]$.

Установлено, что применение алгоритма с адаптацией сетки к особенности функции $a\left(t,x\right)$ приводит к тем же результатам, что и применение алгоритма с адаптацией к особенности функции $u\left(t,x\right)$. В данном случае при $t\in \left[\mathrm{0.1,0.6}\right]$ приближенные значения координат особых точек функций u и a вычислялись устойчиво и отличались друг от друга на величину порядка ${10}^{-2}$. Так, при $t=0.6$ координата особой точки функции $u\left(t,x\right)$ — $z^{*}\approx 0.607+0.15i$, функции $a\left(t,x\right)$ — $z^{*}\approx 0.605+0.136i$. Отметим, что при использовании аппроксимации (8) при $S=2$ в рамках разработанного алгоритма мы всегда получаем пару комплексно-сопряженных полюсов. Здесь и далее без ограничения общности полагаем, что особая точка лежит в верхней комплексной полуплоскости.

В качестве верификации предложенных модели и алгоритма численные решения задачи (1), (3) при W = 10^–4, Re = 1000 и b = 0 сопоставлялись с решениями уравнения Бюргерса (2) с начальными и краевыми данными (3). Напомним, что при $W \to 0$ система (1) преобразуется в уравнение (2). Решения задач (1), (3) и (2), (3) имеют по одному фронту. Относительные отклонения решений при t = 0.6 на фронтах в этом расчете составили 0.167. Вне фронтов величины отклонений имеют порядок 10^–3. Отклонения траекторий особых точек по порядку не превышают 10^–2. Это свидетельствует о том, что качественные свойства модели воспроизводятся в расчетах верно.

Далее была исследована сходимость численных решений с увеличением числа узлов N при различных Re. При этом решения задачи (1), (3) найдены в момент времени $t=t_{\max }=0.6$, и построены графики логарифмов их отклонений (11) в зависимости от $N=N_1$ (фиг. 2). В каждом из проведенных расчетов мы полагали $N_2=N_{\max }+10$, где $N_{\max }$ соответствует максимальному значению $N=N_1$, представленному на графиках. Из анализа графиков на фиг. 2 следует, что характер сходимости является экспоненциальным, а вклад погрешностей, связанных с округлением действительных чисел в памяти компьютера (минимальные значения отклонений для каждого из четырех расчетов), с ростом Re возрастает незначительно.

Отметим, что при увеличении N координаты особых точек изменялись незначительно.

 

Фиг. 2. Десятичный логарифм относительных отклонений решений u(t,x) (а) и a(t,x) (б) при t = 10^–5 / 2, W = 10^–4 и различных Re, равных 10 (1), 100 (2), 1000 (3), 10⁴ (4).

 

Например, в случае Re = 1000, $t=t_{\max }$ отклонения не превышали $2\times {10}^{-3}$. Эти результаты свидетельствуют о том, что особые точки у искомых решений существуют и их координаты c использованием разработанного алгоритма определяются устойчиво.

3.2. Численный анализ течений полимерной жидкости

В ходе численного решения задачи (1), (3) при вариации параметров Re, W и b установлено, что наибольшее влияние на картину течения оказывает параметр W, а также его соотношение с Re. Исследования проведены для достаточно большого диапазона значений W от 10^–5 до 10. Эти значения включают, как случаи малого влияния упругих эффектов ($W\ll 1$), в которых реализуется классическая инерционная и упруго-инерционная турбулентности, так и случаи, когда вклад упругости является существенным ($W\gtrsim 1$). Для последних речь может идти о переходе к чисто упругой турбулентности.

В указанном диапазоне для W при умеренных значениях Re обнаружены два качественно различных режима течения с пороговым значением $W\approx 0.5$. В обоих режимах у скорости формируется большой градиент, что связано с приближением $z^{*}$ к $[-\mathrm{1,1}]$. Однако в первом режиме при W > 0.5 и любых значениях остальных параметров траектория $z^{*}\left(t\right)$ становится практически вертикальной в окрестности $t=0.4$ и особая точка вплотную приближается к $[-\mathrm{1,1}]$ (фиг. 3, линии 3–5). В случае W = 1 при вариации N и t нам удавалось получать корректные результаты вплоть до $t=t^{*}=0.42$. В случаях W = 5 и W = 10 — вплоть до $t=t^{*}=0.415$. Далее при увеличении t у решений наблюдались осцилляции в окрестности $x={\delta }^{*}=0.5$, характерные для явления Гиббса, что свидетельствует о наличии разрыва, причем рвется функция u(t,x) (фиг. 4а, градиенты линий 3–5), а у функции a(t,x) возникает разрыв производной (фиг. 4б, изломы линий 3–5).

 

Фиг. 3. Траектории движения особых точек скорости u(t,x) (а) и функции a(t,x) (б) при Re = 100, b = 0.42 и различных W, равных 0.01 (-.-.-., 1), 0.1 (- - -, 2), 1 (темно-серый цвет, 3), 5 (черный цвет, 4), 10 (светло-серый цвет, 5).

 

Во втором режиме при W < 0.5 картина течения существенно зависит от Re. При значениях Re порядка сотни и меньше величина Im $z^{*}$ достигает минимума и далее практически не изменяется (фиг. 3, линии 1, 2). Моменты времени, при которых достигается минимум, отмечены на фиг. 3. Аналогичные траектории особых точек для u(t,x) наблюдались в [28] при решении задачи (2), (3) для уравнения Бюргерса. Функция a(t,x) имеет пик в окрестности $x=\mathrm{0.5.}$ При этом вершина пика на фиг. 4б (см. линии 1, 2) сглажена. Другая ситуация реализуется, когда W < 0.5 и Re принимает значение больше нескольких сотен. В таком случае особая точка достигает отрезка [–1,1], как и в первом режиме, только несколько позже. При этом момент времени, при котором точка попадает на отрезок [–1,1], $t^{*}$ существенно зависит от Re: чем меньше Re, тем больше $t^{*}$.

 

Фиг. 4. Профили u (0.4, x) (а) и a (0.4, x) (б) при Re = 100, b = 42 и различных W, равных 0.01 (1), 0.1 (2), 1 (3), 5 (4), 10 (5).

 

Отдельно следует сказать о случае $W\ll 1$ и $\mathrm{Re}\gg 1000$. Например, в табл. 1 при Re = 10⁴, $\tilde{\beta }=0.42$ и разных W приведены приближенные значения моментов времени $t^{*}$, когда особая точка достигает отрезка [–1,1], и координат точек отрезка ${\delta }^{*}\in [-\mathrm{1,1}]$, в которые попадает особая точка при $t=t^{*}$. Прочерки в табл. 1 означают, что особая точка не достигает отрезка и в определенный момент времени начинает двигаться параллельно действительной оси в $ℂ.$ Чтобы убедиться в этом, мы продлили расчет до момента времени $t=t_{\max }=1.$

 

Таблица 1. Приближенные значения момента времени появления разрыва решения и координаты разрыва при Re = 10⁴ и различных W

W	tmax	d0
0.1	0.413	0.5005
0.05	0.4158	0.5021
0.01	0.4303	0.5109
0.005	0.443	0.5184
0.001	0.504	0.5539
0.0006	–	–

 

На фиг. 5 приведены траектории движения особых точек при Re = 10⁴, b = 0.42 и при различных W. Пороговое значение W, при котором характер движения особой точки меняется, приблизительно равно 10^–3 (линия 5 на фиг. 5). В этом случае у траектории появляется точка перегиба, характерная для режимов, когда особенность не достигает действительной оси. Однако приближенное решение все равно разрушается со временем, при этом положение особой точки сильно колеблется.

 

Фиг. 5. Траектории особых точек функций u (а) и a (б) при Re = 10⁴, b = 0.42 и разных W, равных 0.1 (черный цвет, 1), 0.05 (темно-серый цвет, 2), 0.01 (серый цвет, 3), 0.005 (светло-серый цвет, 4), 0.001 (-.-.-., 5), 0.0005 (- - -, 6).

 

При дальнейшем увеличении Re наблюдается аналогичная ситуация, но пороговое значение W уменьшается по закону близкому к обратной пропорциональности со значением Re. Например, при Re = 10⁵ и W = 10^–4 особая точка достигает действительной оси, а при W = 10^–5 она начинает двигаться параллельно оси, но вскоре численное решение разрушается.

Число Рейнольдса также оказывает существенное влияние на характер течения. При увеличении Re фронт скорости u становится более крутым и градиент a меняется более резко (фиг. 6). При этом в случае W = 0.5 возникают два режима. При умеренных значениях Re < 100 особая точка остается на небольшом расстоянии от отрезка [–1,1] и разрывов не наблюдается. При Re > 100 особая точка резко приближается к [–1,1], что ведет к скачку значения функции u и производной a. Важно отметить, что время формирования разрыва остается прежним: $t\approx 0.4$.

 

Фиг. 6. Профили u(0.4,x) (а) и a(0.4,x) (б) при W = 0.5, b = 0.42 и различных Re, равных 10 (1), 100 (2), 10⁴ (3).

 

На фиг. 7 приведены профили решения при W = 0.5, Re = 100, $t_{\max }=0.4$ и различных b, которые для скорости u визуально совпадают. У профилей решения для a при разных b отличаются максимальные значения. При этом траектории движения особых точек во всех случаях достигали действительной оси в один и тот же момент времени. Таким образом, вариация параметра b оказывает наименьшее влияние на эволюцию решений.

 

Фиг. 7. Профили u(0.4,x) (а) и a(0.4,x) (б) при W = 0.5, Re = 100 и различных b, равных 0 (1), 0.14 (2), 0.42 (3), 0.7 (4), 1.26 (5).

 

4. Заключение

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

Вопросы о существовании, единственности и типе особой точки, а также о непрерывной зависимости ее координат от числа степеней свободы используемых аппроксимаций выше не обсуждались. Отметим, что необходимых и достаточных условий существования-единственности особой точки приближаемой функции в рамках теории дробно-рациональных аппроксимаций до сих пор установлено не было. Вместе с тем из приведенных графиков видно, что решение имеет один фронт, и барицентрические интерполяции, построенные с учетом одной особой точки, демонстрируют экспоненциальный характер сходимости. Это свидетельствует о существовании и единственности особой точки аналитического продолжения искомого решения в малой окрестности области задачи в $ℂ.$ Следовательно, для ее локализации достаточно использовать приближения Чебышёва–Паде с многочленом второй степени в знаменателе. Алгоритм построения таких приближений является экономичным (см. [27]), что принципиально для решения нестационарных задач.

Важным вопросом остается определение траекторий особенностей в случае, когда их несколько. Для решения этой проблемы можно использовать различные подходы. Повышение степени S знаменателя аппроксимаций Паде вида (8) может привести к плохой обусловленности, численной неустойчивости и ложным полюсам, так называемым дублетам Фруассара (Froissart doublets) (см. [18; 27; 38]). Тем не менее в рассмотренной задаче такие аппроксимации демонстрируют сгущение нулей и полюсов в окрестности единственной найденной нами особой точки, что является еще одним подтверждением ее существования и единственности. Заметим, что для строгого анализа распределения нулей и полюсов помимо решения открытых проблем теории приближения дробно-рациональными функциями требуются неулучшаемые оценки погрешности разработанного метода и теоремы о влиянии этой погрешности на положения нулей и полюсов, что является темой дальнейших исследований. В связи с этим графики распределений нулей и полюсов приближений решения здесь мы не обсуждаем. Отметим только еще один подход к локализации нескольких особенностей в $ℂ,$ состоящий в построении дробно-рациональных приближений с полиномами второй степени в знаменателе (8) на малых подынтервалах исходного отрезка (см. [28]). Предложенное в указанной статье отображение позволяет строить адаптивную сетку для произвольного числа особенностей. Подобная идея обсуждается также в [30], где она реализована с подачи классика теории барицентрических приближений, проф. Жан-Поля Берру (Jean-Paul Berrut).

С точки зрения механики полученные результаты позволяют констатировать качественное различие течений ньютоновской и полимерной жидкостей при достаточно больших W. Прежде всего это касается разрушения решений задач с гладкими начальными данными при умеренных Re. Пороговое значение W, при котором наблюдается этот эффект, есть W ≈ 0.5, а характерное безразмерное время существования решения $t\approx 0.42$. Важно обратить внимание на близость этих величин и на то, что по определению W и t есть соответственно время релаксации и время эволюции течения, нормированные на временной масштаб задачи T_C (см. комментарии к формуле (1)), т.е. если время релаксации превышает 0.5T_C, то время существования гладкого решения меньше 0.5T_C.

Заметим, что режим $W\ll 1$ и $\mathrm{Re}\gg 1000$ представляет течения с большим вкладом вязких и малым вкладом упругих эффектов, в которых возможен переход от классической инерционной к инерционно-упругой турбулентности. В работах [6; 9] показано, что турбулентное течение в таком случае возникает при определенном соотношении между критическими значениями параметров вида Re ~ E(1 –b))^–1/2, где b — отношение вязкостей растворителя и полимерного раствора, E = W / Re — характеристика эластичности (см. [6, рис. 4 и комментарии к нему]). Нетрудно видеть, что эта связь соответствует соотношению $\mathrm{Re}~W^{-1},$ установленному для этого режима и в нашей работе.

Отметим, что, как и в классическом случае для уравнения Бюргерса, увеличение Re приводит к росту градиентов решения, но при этом для достаточно больших W возникают пороговые значения Re, при превышении которых гладкие решения разрушаются, что уже нехарактерно для уравнения Бюргерса с ненулевой вязкостью. Для большей размерности, конечно, вопрос о существовании гладких решений уравнений Навье–Стокса остается открытым и составляет так называемую задачу тысячелетия. Есть все основания полагать, что соответствующая проблема для течений вязкоупругих жидкостей является еще более сложной. В этой связи подчеркнем, что предложенный метод расчета траекторий особых точек, являясь аппроксимационным, не дает необходимых и достаточных условий их существования. Однако он представляет важный шаг на пути к созданию аппарата для численного анализа указанных проблем и может послужить импульсом к развитию новых теоретических подходов. В частности, на основе асимптотического и комплексного анализа с применением теорем о сходимости полиномиальных и дробно-рациональных приближений авторами исследуются возможности вывода систем дифференциальных уравнений, описывающих эволюцию особой точки. Анализ таких систем даст ответы на многие вопросы о существовании гладких решений задач гидродинамики.

Sobre autores

Л. Брындин

Новосибирский государственный университет; Институт теоретической и прикладной механики им. С. А. Христиановича СО РАН

Autor responsável pela correspondência
Email: l.bryndin@g.nsu.ru
Rússia, 630090 Новосибирск, ул. Пирогова, 2; 630090 Новосибирск, ул. Институтская, 4/1

Б. Семисалов

Новосибирский государственный университет; Институт математики им. С. Л. Соболева СО РАН

Email: vibis@ngs.ru
Rússia, 630090 Новосибирск, ул. Пирогова, 2; 630090 Новосибирск, пр-т Акад. Коптюга, 4

В. Беляев

Новосибирский государственный университет; Институт теоретической и прикладной механики им. С. А. Христиановича СО РАН

Email: belyaevasily@mail.ru
Rússia, 630090 Новосибирск, ул. Пирогова, 2; 630090 Новосибирск, ул. Институтская, 4/1

В. Шапеев

Новосибирский государственный университет; Институт теоретической и прикладной механики им. С. А. Христиановича СО РАН

Email: shapeev.vasily@mail.ru
Rússia, 630090 Новосибирск, ул. Пирогова, 2; 630090 Новосибирск, ул. Институтская, 4/1

Bibliografia

Arquivos suplementares