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

Развитие микро- и нанотехнологий актуализирует задачи моделирования течений жидкостей со сложной внутренней структурой [1] . Основы континуальной теории, учитывающей вращательные степени свободы материальных точек, заложены братьями Коссера [2] в начале XX века. Эринген [3, 4] развил и систематизировал эту теорию применительно к жидкостям, предложив модель микрополярной жидкости, которая расширяет классическую гидродинамику, вводя представление о континууме, материальные точки которого обладают собственными вращательными степенями свободы. Модель находит применение при описании суспензий, полимерных растворов, жидкостей с микрочастицами и некоторых биологических сред [5, 6] . Математическая формулировка модели микрополярной жидкости включает дополнительные параметры и уравнение для поля микровращений, что приводит к появлению качественно новых эффектов, таких как зависимость локального вращения частиц от градиентов скорости и возникновение связанных с микровращением напряжений [7, 8] .

Классическая задача течения в полости с движущейся крышкой служит базовым тестом для численных методов в гидродинамике и позволяет исследовать отрывные течения, формирование вихрей и процессы диссипации в ограниченной области [9, 10] . В ряде исследований [11 –14] подобное течение рассматривалось для микрополярных жидкостей, однако вопросы, связанные с детальным анализом полной функции диссипации и количественной оценкой вклада различных микрополярных механизмов в рассеяние энергии, остаются недостаточно изученными.

Диссипация механической энергии в вязкой среде является фундаментальным процессом, определяющим динамику течения, тепловыделение и эффективность множества технических устройств [15, 16] . В случае микрополярной жидкости выражение для функции диссипации существенно усложняется, поскольку включает слагаемые, обусловленные градиентами микровращений и их взаимодействием с полем макроскопической завихренности [17, 18] . В работах [19, 20] показано, что вращение даже в классических (ньютоновских) средах оказывает решающее влияние на устойчивость и вихревую структуру адвективных течений в ограниченных областях с твёрдыми границами. Полный учёт этих механизмов необходим для корректного моделирования тепломассопереноса в микрофлюидных системах, при проектировании микромиксеров и теплообменников [21] , а также для понимания поведения неньютоновских сред [22, 23] .

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

Статья опубликована в открытом доступе по лицензии CC BY 4.0

разрабатывается устойчивый численный алгоритм, вычисляются поля всех искомых функций, включая функцию диссипации, и проводится сравнительный анализ с классическим случаем. Актуальность исследования определяется потребностями современной микрофлюидики и технологий, связанных с управлением течениями на микроуровне [24] .

• 2.    Математическая постановка задачи

Рассматривается двумерное стационарное течение несжимаемой микрополярной жидкости внутри ограниченной прямоугольной области, имеющей форму единичного квадрата. Геометрия расчётной области задаётся как ^ = [0,1] х [0,1] , где оси координат x и у лежат в плоскости течения. Такая постановка позволяет сосредоточить внимание на основных физических особенностях микрополярного течения, не усложняя анализ дополнительными эффектами, связанными с трёхмерной геометрией или изменением свойств жидкости в направлении, перпендикулярном плоскости течения.

Микрополярная жидкость отличается от классической ньютоновской тем, что каждая материальная точка континуума обладает не только поступательной скоростью, но и собственным микровращением, что требует введения дополнительного векторного поля микрополярной скорости (микровращения). В двумерной постановке поле скорости макроскопического движения жидкости описывается двумерным вектором v (x,y) = (u(x,y),v(x,y)) , где u и v — компоненты скорости вдоль осей x и y соответственно. Поле микровращений, в свою очередь, обладает единственной ненулевой компонентой, направленной перпендикулярно плоскости течения, и представляется скалярной функцией, интерпретируемой как проекция вектора и на ось z : и (х,у) = (0,0,w(x,y)) .

В рамках линейной теории микрополярной жидкости определяющие соотношения для тензора силовых напряжений t и тензора моментных напряжений m имеют вид [3, 4, 6, 17] :

∂vi ∂vj tij = -p^ij +М д--+7— 1+ к д--Sijk ^k

∂xj ∂xi∂x

∂ωk      ∂ωi∂ω mij = a^—°j + в д--+ " , .

∂xk      ∂xj∂x

Здесь p — гидродинамическое давление; µ — классическая сдвиговая вязкость; κ — вихревая (вращательная) вязкость; α,β,γ — коэффициенты угловой вязкости; ε ijk — символ Леви–Чивиты; δ ij — символ Кронекера. Для несжимаемой жидкости уравнение неразрывности выражает закон сохранения массы:

∂u ∂v

^V v = dx ⁺ dy =^0.

Подстановка определяющих соотношений в законы сохранения импульса и момента импульса (уравнения Эйлера для сплошной среды) и их последующее обезразмеривание с масштабами длины L (ширина полости) и скорости U (характерная скорость, в качестве которой выбрана скорость движения крышки Ulid), приводят к системе уравнений движения в безразмерной форме. Число Рейнольдса вычисляется по формуле: Re = pUL/^. Уравнение сохранения импульса, обобщающее классическое уравнение Навье–Стокса на случай микрополярных сред, в двумерной постановке записывается как dv .        _ 1+N^2   N Л

di ^{+ (v}•^V) v = ^V R ^V v ⁺-(^Vx . ),

Здесь N = к/р — безразмерный параметр микрополярности, определяющий относительный вклад в движение жидкости вихревой вязкости к по сравнению с классической динамической вязкостью р. Оператор V х и в двумерной постановке даёт двумерный вектор:

Vx w =

что отражает связь между пространственным изменением микровращения и дополнительными объёмными моментами, действующими на жидкость.

Уравнение для поля микровращений (моментная составляющая) в скалярной форме, эквивалентной исходной векторной записи, имеет вид:

^{+ (} v •Яш = —Ц ^V2w + -2^^(z2Ш. dt          Re m 2     Re m 2

В (3) безразмерный параметр m связан с характерной для микрополярной жидкости внутренней длиной l (например, с размером микрочастиц) соотношением: m = l/L , где L — размер полости. Характерная длина l обусловливает величину коэффициентов угловых вязкостей [25] . При l/L = 0.3 внутренний масштаб составляет 30% от размера полости (это типично для микрофлюидных систем, в которых размер канала, как правило, составляет ~ 100 мкм,

• а размер частицы ^ 30 мкм). Параметр Z означает скалярную завихренность макроскопического поля скорости:

  z = (Vx v )• e _г =

  
  

  ∂v ∂u

  ∂x ∂y ^,

  
  

  где e z — единичный вектор, направленный вдоль оси z . Разность (Z — 2ш) интерпретируется как мера несоответствия между локальной завихрённостью макроскопического потока и удвоенным значением микровращения; именно эта величина ответственна за возникновение внутренних моментов, стремящихся уравновесить поступательное и вращательное движениячастиц жидкости. При N = 0 ( к = 0 ) система (2) , (3) тождественно переходит в классические уравнения Навье–Стокса (уравнение (3) становится несвязанным и при нулевых начальных и граничных условиях даёт ш = 0 ). Таким образом, предельный переход к ньютоновской жидкости осуществляется корректно.

  Решение стационарной задачи ищется методом установления — интегрированием нестационарных уравнений (1) - (3) по времени до достижения установившегося состояния. В начальный момент времени t = 0 предполагается, что в области Ω вся жидкость находится в состоянии покоя, а давление и микровращение отсутствуют:

  
  

  ^v(x,y^,0) = ^(0,0), Р^(х,У^,0) = ^{0, м (х,}У^,0) = ⁰.

  
  
  
  

  Особое внимание уделяется рассмотрению диссипативных процессов, связанных с вязким трением как на макро-, так и на микроуровне. Полная функция диссипации энергии для микрополярной жидкости в двумерной постановке записывается в виде [3, 4] :

  
  

  Ф = (1+N)(V v : V v ) + 4 ( V m : V m ) + N (Z— 2ш) 2 . m ²

  
  
  
  

  Первое слагаемое описывает классическую вязкую диссипацию, обусловленную градиентами макроскопической скорости, второе — диссипацию, связанную с пространственными градиентами микровращений, третье — диссипацию, обусловленную несоответствием между макроскопической завихренностью и микровращением. Все три механизма в совокупности определяют полную скорость рассеяния механической энергии в микрополярной среде. Анализ функции диссипации Φ проводится путём оценки локальной скорости рассеяния механической энергии, что в гидродинамике является общепринятым подходом [15, 16] и не требует решения уравнения энергии.

  На верхней стенке ( у = 1 ), движущейся с постоянной скоростью U 0 вдоль оси х , задаются условия прилипания для макроскопической скорости и условие ш = 0 (hyper stick, также известное как no spin), соответствующее полному подавлению собственного вращения микрочастиц на подвижной границе:

  
  

  v =(U o ,0), ш = 0.

  
  
  
  

  Безразмерные переменные вводятся с использованием размера полости L в качестве масштаба длины и вспомогательной скорости U _ref = 2U lid как масштаба скорости и времени, которое при таком обезразмеривании есть L/U _ref. Число Рейнольдса определяется реальной скоростью крышки: Re = pU ljd L/^ . Физика решения полностью соблюдается за счёт использования параметров Re , N , m и не зависит от конкретного выбора масштаба скорости.

  Стенки — нижняя ( у = 0 ) и боковые ( х = 0 и х =1) — полагаются неподвижными:

  
  

  ^v = ^(0,0).

  
  
  
  

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

  
  

  ly=^{0 (}y=⁰⁾,     |X=^{0 (x}=^0,1).

  
  
  
  

  Выбор граничных условий для микровращения обусловлен фундаментальными физическими процессами взаимодействия жидкости со стенками. Так, на движущейся стенке условие ш = 0 соответствует полному подавлению собственного вращения микрочастиц относительно поверхности крышки. На неподвижных стенках условия Неймана дш/дп = 0 физически означают отсутствие внешних микромоментов и свободную адаптацию микровращения к локальной структуре течения. Такая комбинация граничных условий согласуется с рекомендациями в работе [25] .

  Таким образом, полная математическая постановка задачи включает систему уравнений (1) – (3) , начальные условия (4) , граничные условия (6) – (8) , определение функции диссипации (5) и фиксированные значения безразмерных параметров, что обеспечивает однозначное и физически обоснованное описание стационарного микрополярного течения в полости с движущейся крышкой.

  
  

  3. Численное моделирование

  
  

  Для решения задачи стационарного течения микрополярной жидкости в квадратной полости с движущейся крышкой разработан проекционный алгоритм Чорина–Темама [26, 27] с явной временно´ й дискретизацией. Расчётная область Ω покрывается равномерной сеткой. Размер сетки выбран на основе анализа сеточной

  
  

сходимости: выполнены расчёты на сетках с числом узлов 31 х31 , 51 х 51 , 71 х 71 и 101 х 101 для базового набора параметров Re = 10 , N = 0.3 , m = 0.3 . Интегральная диссипация на сетках 51 х 51 и 71 х 71 различается менее чем на 15%, что признаётся приемлемым для качественного анализа. Дальнейшее увеличение сетки до 101 х 101 узлов требует уменьшения шага по времени для получения устойчивого решения и значительно увеличивает время счёта. Таким образом, выбрана сетка 51 х 51 (Дх = Ду = 1 /50 ). Все переменные — компоненты скорости u , v , давление p , микрополярная скорость ω — определяются в центрах ячеек (совмещённая сетка). Временной шаг зафиксирован: ∆t= 10 ^-4 , что обеспечивает устойчивость явной схемы при умеренных значениях числа Рейнольдса. Выбранный временной шаг удовлетворяет необходимому условию устойчивости явной схемы с диффузионными членами.

Уравнение импульса (2) дискретизируется следующим образом: конвективные слагаемые аппроксимируются согласно противопоточной схеме, диффузионные — центральными разностями, производная по времени — явной схемой Эйлера; порядок аппроксимации: первый по времени ( O(∆t) , явная схема Эйлера); второй по пространству для диффузионных членов ( O(∆x ² ) , центральные разности); первый для конвективных членов ( O(∆x) , схема против потока). Граничные условия Дирихле реализуются точно; условия Неймана для функции ω , описыающей поле микровращения, и давления p получаются путём зеркального отражения значений из приграничных узлов, что обеспечивает второй порядок аппроксимации на равномерной сетке.

На каждом временно´ м шаге сначала вычисляется новое поле микровращений ω ^{n +1} ; решается уравнение (3) , записанное в дискретном виде:

_w n +i _n_n ∆t

+( v ⁿ

V h )w ⁿ

—Ц v2 wn+-2N-

Re m ^{2 h} Re m ²

(Z n - 2w n ),

где V _h означает дискретный градиент, а Z n = (dv/dx - du/dy) h есть дискретная скалярная завихрённость. Затем определяются промежуточные компоненты скорости (u ^∗ ,v ^∗ ) без учёта градиента давления. Для восстановления несжимаемости итерационно, методом Якоби с 100 внутренними итерациями, решается уравнение Пуассона для давления:

v h p n ⁺¹ = Д V h • v * .

После этого поле скорости корректируется:

v n ⁺¹ = v * - At V h p ⁿ⁺¹ .

На всех границах задаются условия Неймана ∂p/∂n = 0 . Данное условие возникает естественным образом в проекционном методе при выполнении условия прилипания на твёрдых стенках [26, 27] и является стандартной практикой в расчётах несжимаемых вязких течений.

Граничные условия применяются на каждом подэтапе. На движущейся верхней стенке полагается: u=U 0 , v = 0 , ω = 0 ; на неподвижных стенках u= v = 0 за счёт зеркального отражения значений из приграничных узлов (например, ω 1 ,j = ω 2 ,j) для ω реализуются условия нулевого нормального градиента.

Для повышения устойчивости вводится подрелаксация поля микровращений шп+1 ^ awn+np |(l a)wn с коэффициентом α = 0.3. Здесь ωcno+m1p — значение, полученное непосредственно из разностного уравнения. Итерационный процесс останавливается при достижении максимального изменения всех полей менее ε = 10-6 или после 1000 шагов по времени.

Дополнительно сделана оценка устойчивости решения по параметру m . Расчёты показали, что при m = 0.15 и m = 0.16 численная схема теряет устойчивость (возникает неограниченный рост невязок, а вычисления становятся некорректными — NaN), тогда как при m= 0.18 и выше решение устойчиво. Таким образом, нижняя граница применимости метода составляет m ⩾ 0.18 . В параметрическом исследовании использовались значения m = {0.2,0.4,0.6,0.8} .

Проведена проверка корректности реализации алгоритма: подтверждено точное выполнение граничных условий; сохранение несжимаемости ( |V- v | ^ < 10 ^-6 ); монотонное убывание невязок до уровня 10 ^-5 (Рис. 1) ; качественное совпадение структуры течения с классическим решением при N ^ 0 .

• 4. Результаты и обсуждение

  Численное моделирование показало устойчивую сходимость решения к физически корректному стационарному состоянию. Невязки всех полей монотонно убывали на протяжении 1000 итераций (Рис. 1) , демонстрируя стабильность использованной схемы даже при наличии сильной нелинейности и дополнительной степени свободы — микровращения. Несмотря на то, что расчёт не достиг полной сходимости по критерию ε = 10 ^-6 , достоверность анализа стационарной структуры течения подтверждают малые остаточные невязки ( 7.6 х 10 ^-5 для u ), физическая согласованность распределений всех вычисленных полей и положительная определённость функции диссипации.

Для количественного анализа асимметрии течения с использованием разработанного численного алгоритма выполнены расчёты при трёх значениях числа Рейнольдса Re = {10,500,2000} с фиксированными N = 0.3 , m = 0.3

Рис. 1. Зависимость максимальной невязки компонент скорости u , v и микрополярной скорости ω от номера итерации (базовый случай Re= 10 , N = 0.3 , m = 0.3 )

(см. Рис. 2) . Центр основного вихря, определяемый по минимуму функции тока, при Re= 10 находится в точке (0.50; 0.82). С ростом Re до 500 центр смещается до (0.52; 0.88), а при Re = 2000 достигает точки (0.52; 0.92). Таким образом, увеличение числа Рейнольдса приводит к смещению вихря вверх и вправо, что свидетельствует об усилении асимметрии течения, обусловленной движением верхней стенки. Изолинии функции тока выявляют основной вихрь, смещённый вверх; вторичные вихри в верхних углах не обнаружены. Аналогичные эффекты, связанные с влиянием вращения на структуру течения, обсуждаются в работах [28, 29] .

( а )

Рис. 2. Линии тока при N = 0.3 , m=0.3 и разном числе Рейнольдса: 10 ( а ); 500 ( б ); 2000 ( в )

( в )

Профили горизонтальной скорости u(y) при x = 0.5 для шести значений Re = {10,100,500,1000,1500,2000} приведены на рисунке 3. Как видно, с увеличением числа Рейнольдса градиент скорости у движущейся верхней стенки возрастает, а профиль в ядре течения становится более пологим. Это соответствует классической картине формирования тонкого пограничного слоя при высоких Re и подтверждает физическую корректность выполненных расчётов.

Профили микровращения ω(y) при x = 0.9 (область, где ω достигает максимума) показаны на рисунке 4. Микровращение локализовано вблизи твёрдых стенок и быстро затухает в ядре течения. С увеличением Re максимальное значение ω уменьшается, что объясняется подавлением микровращений возрастающими инерционными силами.

Профили диссипации Φ(y) при x = 0.5 и x = 0.9 представлены на рисунке 5. Диссипация сосредоточена в узких пристеночных слоях, особенно вблизи верхней стенки, где градиенты скорости максимальны. С ростом Re максимальное значение диссипации увеличивается, а область её локализации становится более узкой. Интегральная диссипация для базового случая составляет 2.133 Вт/м. Для количественной оценки влияния микрополярных эффектов проведён расчёт для классической ньютоновской жидкости ( N = 0 ); в этом случае интегральная диссипация равняется 1.812 Вт/м. Таким образом, учёт микровращений при выбранных значениях параметров ( N = 0.3 , m = 0.3 ) выявляет заметное увеличение диссипации.

С целью анализа влияния параметров на диссипацию энергии для течения микрополярой жидкости проведена серия расчётов при следующих значениях параметров:

• -    число Рейнольдса Re = {10,100,500,1000,1500,2000} ;

• —    параметр микрополярности N = {0,0.1,0.3,0.5,0.7} ;

  

  Рис. 3. Профили горизонтальной скорости u(y) при x = 0.5

  
  
  

  Рис. 4. Профили микровращения w(y ) при x = 0.9

  
  
  

  Рис. 5. Профили диссипации Ф(у) при разных значениях координаты х : 0.5 ( а ); 0.9 ( б )

  
  
  
  
  

• —    параметр микромасштаба m = {0.2,0.4,0.6,0.8} .

На рисунке 6 показана зависимость интегральной диссипации от числа Рейнольдса для различных комбинаций (N,m) . Диссипация возрастает с увеличением Re при всех исследованных значениях N и m . Наибольший рост наблюдается в диапазоне Re = 10, ..., 500 , а при дальнейшем увеличении Re до значения 2000 диссипация выходит на плато или слегка снижается, что объясняется насыщением диссипативных механизмов.

Рисунок 7 демонстрирует связь диссипации с параметром микрополярности N . Зависимости, в отличие от свойственных для классических ньютоновских жидкостей, не являются строго монотонными. На большинстве кривых видны локальные максимумы и минимумы (зигзагообразный характер). Например, при Re = 2000 и m = 0.2 диссипация возрастает при увеличении N от 0 до 0.1, достигает максимума, а затем снижается. Такое поведение объясняется нелинейной обратной связью между макроскопической завихренностью ζ и полем микровращения ш в выражении (5) . При увеличении N усиливаются первый (1 + N )(V v : V v ) и третий N (Z — 2ш) ² члены, но одновременно перестраивается поле ш , что приводит к локальному уменьшению разности (Z — 2ш) и, как следствие, к снижению вклада третьего слагаемого, поскольку его квадратичная часть при этом падает и частично компенсирует рост множителя N . Несмотря на локальные немонотонности, общий диапазон изменения диссипации при увеличении N от 0 до 0.7 составляет от 1.12 до 4.94 Вт/м (меняется более чем в 4 раза), что подтверждает ключевую роль параметра N .

На рисунке 8 представлена зависимость диссипации от параметра m . Как и на рисунке 7, здесь наблюдаются немонотонные колебания (зигзаги). Причина та же: изменение m влияет на коэффициент N/m ² во втором слагаемом функции диссипации, что через перестройку поля ω вызывает нелинейную реакцию всей системы. Однако абсолютное изменение диссипации во всём диапазоне m = 0.2,..., 0.8 не превышает 10-15%, тогда как при варьировании Re и N оно достигает сотен процентов. Таким образом, влияние m на интегральные энергетические характеристики является вторичным по сравнению с влиянием Re и N .

—e—^=0.0,/и=0.2 - G-^=0.3,/7/=0.6 —H—^=0.0, ш=0.4 - a- ^=0.3, 777=0.8 —А—^=0.0,т=0.6 —A—^=0.5,//7=0.2 —0—^=0.0, /77=0.8 —0—^=0.5, m=0.4 —V—^0.1, 777=0.2 “V—^=0.5, /77=0.6 —^^=0.1, 777=0.4 -^-^=0.5, 777=0.8 - * -^=0. l,/77=0.6 —«—^=0.7, 777=0.2 - +■ -^=0.1, 777=0.8 -+-^=0.7, 777=0.4 - * -^=o.3,/77=0.2 ■■■■■*■■■■ W=0.7, 777=0.6 - * -^=0.3, 777=0.4 —*-^=0.7, 777=0.8

Re

Рис. 6. Зависимость интегральной диссипации от числа Рейнольдса для различных комбинаций (N,m)

5.5

0      0.1     0.2     0.3     0.4     0.5     0.6     0.7

2V

Рис. 7. Зависимость интегральной диссипации от параметра N для различных комбинаций (Re,m)

Рис. 8. Зависимость интегральной диссипации от параметра m для различных комбинаций (Re,N)

• 5.    Выводы

Проведено численное исследование установившегося двумерного течения микрополярной жидкости в квадратной полости с движущейся крышкой. Разработанный алгоритм на основе проекционного метода показал свою эффективность. Выполнено параметрическое исследование (осуществлено 120 расчётов при разных комбинациях параметров). Обнаружено, что зависимости диссипации от N и m имеют немонотонный (зигзагообразный) характер, что связано с нелинейной обратной связью между макроскопической завихренностью и микровращением. При этом влияние параметра m на диссипацию существенно слабее, чем влияние Re и N . Максимумы диссипации локализованы в пристеночных слоях; распределение микрополярной скорости коррелирует с областями высоких градиентов макроскопической скорости. Полученные авторами данные подтверждают, что микрополярные эффекты существенно влияют на диссипацию энергии; их вклад должен учитываться при анализе вязких течений в замкнутых областях.

Исследование выполнено при поддержке Российского научного фонда (проект № 25-29-00339), 29- 00339/.