Self-similar solution of the first Stokes problem for non-Newtonian fluids with power-law viscosity

Full Text

Среди нестационарных задач динамики вязкой несжимаемой ньютоновской жидкости с постоянным значением вязкости хорошо известна так называемая первая задача Стокса [1]. В этой задаче рассматривается неустановившееся ламинарное течение жидкости в полупространстве, которое ограничено пластиной, внезапно приведенной в своей плоскости в поступательное прямолинейное движение с постоянной скоростью. Такая задача допускает достаточно простое автомодельное решение [2]. Кроме самостоятельного интереса при проведении моделирования соответствующих гидродинамических процессов, это решение может быть использовано для тестирования различных численных методов [3, 4].

Вместе с тем, широкий круг реальных жидкостей демонстрирует более сложное механическое поведение по сравнению с ньютоновскими жидкостями [5–7]. Вязкость таких сплошных сред изменяется в зависимости от скорости деформирования. В первом приближении эта зависимость, зачастую, может быть аппроксимирована в рамках степенного закона. Именно такой вид был принят в одних из первых и хорошо известных в настоящее время нелинейных реологических моделях Оствальда – де Ваэля [8, 9] и Гершеля–Балкли [10]. Такие модели постоянно используются при рассмотрении самых различных задач гидродинамики [11, 12].

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

Постановка задачи

Рассмотрим одномерное, ламинарное течение вязкой несжимаемой неньютоновской жидкости, возникающее в полупространстве, ограниченном пластиной, которая мгновенно приводится в своей плоскости в поступательное прямолинейное движение с постоянной скоростью U₀. Введем декартову систему координат традиционным образом, сориентировав ось Ox вдоль пластины в направлении движения, а ось Oy – по нормали к ней. В безразмерной форме уравнение, описывающее динамику жидкости в рамках таких допущений, записывается следующим образом:

$\frac{\partial u^{\text{'}}}{\partial t^{\text{'}}}=\frac{\partial {{\tau }^{\text{'}}}_{xy}}{\partial y^{\text{'}}};$$u^{\text{'}}=\frac{u}{U_0};$$t^{\text{'}}=\frac{t{U}_0}{L};$${{\tau }^{\text{'}}}_{xy}=\frac{{\tau }_{xy}}{\rho U_0^2};$${{\tau }^{\text{'}}}_{xy}=\frac{{\tau }_{xy}}{\rho U_0^2};$  (1.1)

где t – время; u = u(y, t) – скорость жидкости в направлении продольной оси Ox ; τ_xy – соответствующая компонента тензора напряжений; ρ – плотность жидкости; L – некоторый характерный, принимаемый в качестве масштабного, линейный размер области течения.

Здесь и далее безразмерные величины отмечены верхними штрихами.

Начальное условие, а также одно из граничных условий (на поверхности движущейся пластины) имеют вид

$t^{\text{'}}=0;$$u^{\text{'}}=0;$$y^{\text{'}}=0;$ (1.2)

$y^{\text{'}}=0;$$u^{\text{'}}=1;$$t^{\text{'}}\ge 0.$ (1.3)

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

Случай течения жидкости Оствальда – де Ваэля

Пусть механическое поведение жидкости удовлетворяет реологической модели Оствальда – де Ваэля, в соответствии с которой в рамках рассматриваемой постановки задачи (1.1)–(1.3) имеем

${\tau }_{xy}={\mu }_{eff}\frac{\partial u}{\partial y};$${\mu }_{eff}=k{\left|\frac{\partial u}{\partial y}\right|}^{n-1},$ (2.1)

где µ_eff – эффективная вязкость; k – коэффициент консистенции; n – индекс течения.

Полагая априори, что

$\frac{\partial u}{\partial y}<\mathrm{0,}$

приходим к следующей безразмерной форме записи выражения (2.1) для касательного напряжения

 ${{\tau }^{\text{'}}}_{xy}=-k^{\text{'}}{\left(-\frac{\partial u^{\text{'}}}{\partial y^{\text{'}}}\right)}^n,$$k^{\text{'}}=\frac{k}{\rho L^n{U}_0^{2-n}}.$ (2.2)

Заметим следующее. Рассматриваемая схема течения не предполагает каких-то конкретных характерных размеров. В этой связи, в качестве L предлагается выбирать величину, получаемую некоторой комбинацией прочих физических (но не геометрических) параметров задачи. Например, принимая во внимание основные положения теории размерностей, в качестве характерного размера может быть принято значение

$L={\left(\frac{k}{\rho U_0^{2-n}}\right)}^{\frac{1}{n}}.$ (2.3)

Тогда с учетом (2.3) из (2.2) получаем, что

$k^{\text{'}}=1.$ (2.4)

Подставляя теперь (2.2) с учетом (2.4) в (1.1), приходим к уравнению

$\frac{\partial u^{\text{'}}}{\partial t^{\text{'}}}=-\frac{\partial }{\partial y^{\text{'}}}{\left(-\frac{\partial u^{\text{'}}}{\partial y^{\text{'}}}\right)}^n.$ (2.5)

Здесь в дополнение к (1.2), (1.3) должно быть поставлено еще одно граничное условие:

 $y^{\text{'}}\to \infty ;$$u^{\text{'}}=0.$ (2.6)

Введем в рассмотрение следующую безразмерную автомодельную переменную:

$\eta =\eta (y^{\text{'}},t^{\text{'}})=y^{\text{'}}{\left(2n(n+1)t^{\text{'}}\right)}^{-\frac{1}{n+1}}.$ (2.7)

Тогда с учетом (2.7) уравнение (2.5) сводится к обыкновенному дифференциальному уравнению относительно u´(η) следующего вида:

$\frac{d^2{u}^{\text{'}}}{d{\eta }^2}-2\eta {\left(-\frac{d{u}^{\text{'}}}{d\eta }\right)}^{2-n}=0.$ (2.8)

При этом исходные начальное (1.2) и граничные условия (1.3), (2.6) преобразуются и записываются в форме

$\eta =0;$$u^{\text{'}}=1;$ (2.9)

$\eta \to \infty ;$$u^{\text{'}}=0.$ (2.10)

Отметим, что при n = 1 соотношения (2.8)–(2.10) приводят к постановке классической первой задачи Стокса для случая течения ньютоновской жидкости с постоянным значением вязкости.

Решение (2.8) с учетом условия (2.9) может быть представлено следующим образом:

$u^{\text{'}}\left(\eta \right)=1-\underset{0}{\overset{\eta }{\int }}{\left[(1-n)\left({\xi }^2+C\right)\right]}^{\frac{1}{n-1}}d\xi .$ (2.11)

Здесь С – неопределенная пока константа интегрирования, которая, принимая во внимание условие (2.10), должна определяться из решения уравнения

$\underset{0}{\overset{\infty }{\int }}{\left[\left(1-n\right)\left({\xi }^2+C\right)\right]}^{\frac{1}{n-1}}d\xi =1.$ (2.12)

Анализ (2.12) приводит к следующему заключению. В случае дилатантного поведения рассматриваемой жидкости (n > 1) это уравнение не имеет решения (теряет смысл), поскольку не существует конечного отрицательного значения константы С, при котором на всем полуинтервале [0, ∞) интегрирования основание степени в квадратных скобках оставалось бы положительным. Иначе говоря, полученное автомодельное решение (2.11) может быть отнесено только к случаю жидкости с псевдопластическим поведением, когда 0 < n < 1.

В общем случае определить из (2.12) константу интегрирования аналитически представляется затруднительным за исключением некоторых частных значений параметра n. Так, при n = 0.5 имеем С = π^2/3. В этой связи решение уравнения (2.12) проводили численно.

В качестве иллюстрации на рис. 1 представлена полученная зависимость значений константы интегрирования С от показателя степени n в реологической модели Оствальда – де Ваэля. Из рассмотрения этого графика следует, что функция C(n) на диапазоне 0.1 < n < 1 наиболее часто встречающихся на практике значений параметра n не является монотонной, достигая в некоторой точке n ≈ 0.271 экстремума типа минимум. При этом значению n = 1 соответствует вертикальная асимптота.

Рис. 1. Зависимость константы интегрирования C в автомодельном решении (2.11) от показателя степени n в реологической модели для жидкости Оствальда – де Ваэля, которая демонстрирует псевдопластическое поведение.

Случай течения жидкости Гершеля–Балкли

Пусть механическое поведение среды удовлетворяет реологической модели вязкопластической жидкости Гершеля–Балкли.

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

Априори можно полагать, что граница раздела этих зон представляет собой параллельную пластине плоскость, которая перемещается в пространстве, удаляясь от пластины, в направлении оси Oy. При этом координата этой границы раздела представляет собой неизвестную заранее функцию времени y_p = y_p(t).

В безразмерной форме записи с учетом допущений и обозначений предыдущего раздела в зоне сдвигового течения выражение для касательного напряжения согласно модели Гершеля–Балкли может быть представлено следующим образом:

${{\tau }^{\text{'}}}_{xy}=-{{\tau }^{\text{'}}}_p-k^{\text{'}}{\left(-\frac{\partial u^{\text{'}}}{\partial y^{\text{'}}}\right)}^n;$${{\tau }^{\text{'}}}_p=\frac{{\tau }_p}{\rho U_0^2}.$ (3.1)

Тогда с учетом (2.4), (3.1) вновь приходим из (1.1) к уравнению вида (2.5).

Учитывая наличие границы раздела между застойной зоной и зоной сдвигового течения, в дополнение к (1.3) здесь должны быть поставлены следующие условия:

$y^{\text{'}}={y^{\text{'}}}_p;$$u^{\text{'}}=0;$$\frac{\partial u^{\text{'}}}{\partial y^{\text{'}}}=0.$  (3.2)

Последнее граничное условие в (3.2) означает, что касательное напряжение на границе раздела зон по модулю принимает значение предела текучести.

Вводя теперь в рассмотрение автомодельную переменную (2.7), вновь приходим к обыкновенному дифференциальному уравнению вида (2.8).

При этом исходные начальное (1.2) и граничные условия (1.3), (3.2) преобразуются и записываются в форме

$\eta =0;$$u^{\text{'}}=1;$ (3.3)

$\eta ={\eta }_p;$$u^{\text{'}}=0;$$\frac{d{u}^{\text{'}}}{d\eta }=0.$  (3.4)

Здесь принято обозначение

${\eta }_p=\frac{{y^{\text{'}}}_p\left(t^{\text{'}}\right)}{\sqrt[n+1]{2n\left(n+1\right)t^{\text{'}}}}.$ (3.5)

Отметим, что начальное условие (1.2), которое сводится к виду (2.10), выполняется в точности, поскольку в застойной зоне при η_p < η < ∞ заведомо имеем u′= 0.

Решение уравнения (2.8) для распределения скорости в зоне сдвигового течения описывается соотношением вида

$u^{\text{'}}\left(\eta \right)=1-\underset{0}{\overset{\eta }{\int }}{\left[(n-1)\left({\eta }_p^2-{\xi }^2\right)\right]}^{\frac{1}{n-1}}d\xi .$ (3.6)

Анализируя выражение (3.6), можно видеть, что в случае 0 < n < 1 не представляется возможным удовлетворить краевым условиям (3.3), (3.4). В этой связи построение автомодельного решения для рассматриваемой задачи возможно лишь в случае n > 1, когда жидкость Гершеля–Балкли демонстрирует дилатантное поведение.

В (3.6) значение η_p автомодельной переменной на границе раздела зоны сдвигового течения и застойной зоны, находится с учетом условия (3.3) в зависимости от значения параметра n из решения следующего уравнения:

$\underset{0}{\overset{{\eta }_p}{\int }}{\left[\left(n-1\right)\left({\eta }_p^2-{\xi }^2\right)\right]}^{\frac{1}{n-1}}d\xi =1.$ (3.7)

Определить аналитически значение η_p из решения (3.7) получается лишь для некоторых частных значений индекса течения. Например, при n = 1.5 имеем точное значение η_p = (7.5)^0.2. Для произвольного же значения индекса течения найти η_p аналитически представляется затруднительным.

В этой связи решение уравнения (3.7) проводили численно. Результаты такого решения представлены на рис. 2 в виде графика зависимости η_p от n. Заметим, что при n ≈ 7.199 график функции η_p(n) имеет экстремум типа минимум. При этом значению n = 1 соответствует вертикальная асимптота.

Рис. 2. Зависимость значения η_p автомодельной переменной, которая определяет положение границы раздела между зоной сдвигового течения и застойной зоной для жидкости Гершеля–Балкли, от показателя степени n.

Определившись с решением η_p(n) уравнения (3.7), из (3.5) приходим к выражению

${y^{\text{'}}}_p\left(t^{\text{'}}\right)={\eta }_p\left(n\right)\sqrt[n+1]{2n(n+1)t^{\text{'}}},$

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

Заключение

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

Рассмотрены две хорошо известные реологические модели такого типа. Применительно к модели Оствальда – де Ваэля показано, что полученное автомодельное решение для распределения скорости распространяется только на случай, когда жидкость демонстрирует псевдопластическое поведение (0 < n < 1).

Что же касается модели жидкости Гершеля–Балкли, то здесь, наоборот, автомодельное решение для распределения скорости может быть получено лишь для случая, когда жидкость проявляет свойство дилатантности (n > 1).

About the authors

V. N. Kolodezhnov

Author for correspondence.
Email: kvn117@mail.ru
Russian Federation, Voronezh, 394064

References

Supplementary files

Supplementary Files

Action

1. JATS XML

Download

2. Fig. 1. Dependence of the integration constant C in the self–similar solution (2.11) on the exponent n in the rheological Ostwald-de Wael fluid model, which demonstrates pseudoplastic behavior.

Download (48KB)

Indexing metadata

3. Fig. 2. Dependence of the np value of the self-similar variable, which determines the position of the interface between the shear flow zone and the stagnant zone for the Herschel–Bulkley liquid, on the exponent n.

Download (47KB)

Indexing metadata