Finite difference simulation of three-dimensional micropolar fluid flow in a cubic cavity with a moving upper wall

Full Text

ВВЕДЕНИЕ

Развитие численных методов для моделирования течений микроструктурированных жидкостей представляет значительный интерес для современной гидродинамики [17; 37]. Классические уравнения Навье – Стокса оказываются недостаточными для описания жидкостей со сложной внутренней структурой, таких как полимерные растворы, суспензии, жидкости с микрочастицами и биологические жидкости [31]. Теория микрополярных жидкостей, предложенная А. К. Эрингеном [14; 15], расширяет классическую гидродинамику, вводя независимые вращательные степени свободы микроэлементов среды, что позволяет более адекватно описывать перенос момента количества движения [4; 13].

Численное моделирование трёхмерных течений микрополярных жидкостей сопряжено с рядом методических трудностей [31]. Уравнения движения включают дополнительные уравнения для микрополярной скорости, которые связаны с уравнениями Навье – Стокса через нелинейные члены. Это требует разработки специальных алгоритмов, обеспечивающих устойчивость и точность решения. Особую сложность представляет обеспечение сходимости решения при наличии конвективных членов и сильной связи между уравнениями [18; 32].

Классическая задача Куэтта о течении между параллельными пластинами, одна из которых движется с постоянной скоростью, служит фундаментальным тестовым случаем для верификации численных методов [6; 21]. Для ньютоновской жидкости это течение имеет простое аналитическое решение, что позволяет оценить точность численного алгоритма. Для микрополярной жидкости аналитические решения существуют лишь для упрощённых одномерных постановок [5; 10], что делает задачу верификации трёхмерных решений особенно актуальной. Однако в большинстве практических приложений, связанных с микрофлюидными устройствами, течение происходит не в бесконечных каналах, а в замкнутых полостях конечных размеров, где наличие боковых стенок существенно влияет на структуру течения.

В существующей литературе основное внимание уделяется разработке численных методов для двумерных течений микрополярных жидкостей [3; 20]. Трёхмерные постановки исследованы значительно меньше, а вопросы верификации численных алгоритмов для трёхмерных течений освещены недостаточно полно. Отсутствие стандартизированных тестовых случаев для трёхмерных микрополярных течений затрудняет сравнение различных численных методов и оценку их точности.

Разработанные ранее алгоритмы для микрополярных жидкостей часто используют упрощённые подходы к дискретизации конвективных членов или применяют методы, не гарантирующие сохранения свойств исходных уравнений [30; 28]. Это может приводить к физически некорректным результатам или ухудшению сходимости при определённых значениях параметров микрополярности.

Новизна настоящей работы заключается в следующем. Впервые предложен и обоснован конечно-разностный алгоритм для трёхмерного течения микрополярной жидкости в замкнутой кубической полости с подвижной верхней стенкой, в котором:

– использована устойчивая комбинация метода проекций для поля скорости и явной схемы для уравнений микровращения с сохранением второго порядка точности по пространству;

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

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

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

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

Задачи исследования:

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

Математическая модель

Уравнения движения микрополярной жидкости записываются в безразмерной форме с использованием следующих обозначений: u_i – компоненты вектора скорости, p – давление, ω_i – компоненты вектора микрополярной скорости. Система уравнений включает уравнение неразрывности, уравнения баланса количества движения и уравнения баланса момента количества движения [14; 15; 23]:

$\frac{\partial u_i}{\partial x_i}=0$ (1)

$\frac{\partial u_i}{\partial t}+u_j\frac{\partial u_i}{\partial x_j}=-\frac{\partial p}{\partial x_i}+\frac{1+N}{Re}\frac{{\partial }^2{u}_i}{\partial x_j\partial x_j}+\frac{N}{Re}{ϵ}_{ijk}\frac{\partial {\omega }_k}{\partial x_j}$ (2)

$\frac{\partial {\omega }_i}{\partial t}+u_j\frac{\partial {\omega }_i}{\partial x_j}=\frac{1}{Re\cdot m^2}\frac{{\partial }^2{\omega }_i}{\partial x_j\partial x_j}+\frac{2N}{Re\cdot m^2}\left({ϵ}_{ijk}\frac{\partial u_k}{\partial x_j}2{\omega }_i\right)$ (3)

где $Re=\frac{\rho U_{0}L}{\mu }$ – число Рейнольдса, $N=\frac{\kappa }{\mu +\kappa }$ – параметр микрополярности (0 ≤ N < 1), $m=L\sqrt{\frac{\kappa +\mu }{\gamma }}$ – параметр, характеризующий соотношение макро- и микромасштабов, μ – динамическая вязкость, κ – коэффициент вихревой вязкости, γ – коэффициент спин-градиентной вязкости, ϵ_ijk – символ Леви-Чивиты.

Диссипативная функция Φ для микрополярной жидкости определяется выражением [17; 35; 37]:

$\text{Φ}=\left(1+N\right)D_{ij}{D}_{ij}+\frac{1}{m^2}\frac{\partial {\omega }_i}{\partial x_j}\frac{\partial {\omega }_i}{\partial x_j}+N{\left({ϵ}_{ijk}\frac{\partial u_k}{\partial x_j}2{\omega }_i\right)}^2$ (4)

где $D_{ij}=\frac{1}{2}\left(\frac{\partial u_i}{\partial x_j}\text{+}\frac{\partial u_j}{\partial x_i}\right)$ – тензор скоростей деформации.

Рассматривается трёхмерная замкнутая кубическая полость со стороной L. На нижней стенке (z = 0) и боковых стенках (x = 0, L; y = 0, L) задаются условия прилипания для скорости (u = v = w = 0). На верхней стенке (z = L) задаётся скорость u = U₀, v = w = 0. Для микрополярной скорости на всех твёрдых границах принимаются условия Неймана: ∂ω_i /∂n = 0 [11]. Данный выбор обусловлен рядом физических и вычислительных аспектов. В микрополярной жидкости граничные условия для микровращения не являются столь же однозначно определёнными, как условия прилипания для скорости [40]. Условия Неймана соответствуют физическому предположению об отсутствии потока момента количества движения через твёрдую стенку, что часто используется при моделировании течений микрополярных жидкостей в замкнутых областях [29; 40]. Альтернативный подход – условия Дирихле (ω_i=0 на стенке) – соответствует полному «прилипанию» микровращения и также используется в литературе [15], однако в данной работе предпочтение отдано условиям Неймана, поскольку они обеспечивают более устойчивое численное решение и лучше согласуются с постановкой задачи для скорости (верхняя стенка движется, нижняя и боковые неподвижны). Как показано в [11; 15], для течений с подвижными границами условия Неймана позволяют корректно описать распределение микровращения в пристеночных областях.

Численный метод

Для численного решения системы уравнений (1)–(3) разработан алгоритм на основе метода конечных разностей [16; 29]. Пространственная дискретизация выполняется на равномерной сетке с шагами Δx, Δy, Δz. Производные по пространственным переменным аппроксимируются центральными разностями второго порядка точности [24]:

$\frac{\partial f}{\partial x}\approx \frac{f_{i+\mathrm{1,}j,k}-f_{i-\mathrm{1,}j,k}}{2\text{Δ}x},\frac{{\partial }^{2}f}{\partial x^2}\approx \frac{f_{i+\mathrm{1,}j,k}-2f_{i,j,k}+f_{i-\mathrm{1,}j,k}}{\text{Δ}{x}^2}$

Для конвективных членов в уравнениях (2) и (3) используется устойчивая схема против потока, которая обеспечивает монотонность решения при больших числах Рейнольдса [18]. Временная дискретизация осуществляется по явной схеме с шагом Δt, удовлетворяющим условию устойчивости Куранта – Фридрихса – Леви [1].

Алгоритм решения реализует метод проекций, который состоит из трёх основных этапов [12; 16]. На первом этапе вычисляется промежуточное поле скорости без учёта градиента давления. На втором этапе решается уравнение Пуассона для давления:

$\frac{{\partial }^{2}p}{\partial x_i\partial x_i}=\frac{1}{\text{Δ}t}\frac{\partial u_i^{*}}{\partial x_i}$ (5)

где u_i* – промежуточная скорость. Уравнение (5) решается итерационным методом с использованием семиточечного шаблона для трёхмерного лапласиана. На третьем этапе выполняется коррекция скорости с учётом найденного градиента давления.

Для решения системы уравнений микровращения (3) применяется явная схема с использованием тех же аппроксимаций для конвективных и диффузионных членов. Связь между уравнениями скорости и микровращения учитывается через члены, содержащие символ Леви-Чивиты.

Критерием сходимости стационарного решения является выполнение условий:

$max\mid u^{n+1}-u^n\mid <\epsilon , max\mid {\omega }^{n+1}-{\omega }^n\mid <\epsilon$ (6)

где ε = 10^-5 – заданная точность, n – номер итерации. Расчёты прекращаются при выполнении условий (6) или при достижении максимального числа итераций N_max=500.

Алгоритм реализован в среде MATLAB [19]. Для интегрирования диссипативной функции (4) по объёму используется метод прямоугольников с весом, равным объёму элементарной ячейки сетки.

Верификация разработанного численного метода выполнена по нескольким направлениям.

Сравнение с аналитическим решением для ньютоновской жидкости (N=0). В центральном сечении полости (вдали от боковых границ) профиль скорости должен быть близок к линейному u=U₀ z/L, характерному для течения Куэтта. Расчёты для сетки 31×31×31 при Re = 5, U₀ = 0.2 дали относительную погрешность не выше 2.4×10^-4, что подтверждает корректность дискретизации [6; 21].

Исследование сходимости при сгущении сетки. Проведены расчёты на сетках 16³, 31³, 61³. Установлен второй порядок сходимости по пространственным переменным, что согласуется с теорией метода конечных разностей [24; 30]. Норма погрешности для скорости убывает пропорционально h².

Выполнение закона сохранения массы. Расходимость поля скорости в стационарном решении не превышает 10^-8, что свидетельствует о корректной реализации уравнения неразрывности [32].

Сходимость итерационного процесса. На рисунке 1 представлена зависимость невязок для скорости и микровращения от номера итерации. Невязки монотонно уменьшаются и достигают ε=10^-5 за 300–400 итераций. Скорость сходимости для микровращения несколько ниже, что обусловлено наличием дополнительных нелинейных членов в уравнении (3) [7; 26].

 

Рисунок 1. Сходимость численного решения: зависимости невязок для скорости и микровращения от номера итерации

 

Физические законы сохранения. Проверка баланса кинетической энергии и момента количества движения показала, что диссипация энергии соответствует работе внешних сил с точностью до 0.1 % [2].

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

Базовый набор параметров для верификационных расчётов: Re = 5, N = 0.2, m = 0.3, U₀ = 0.2, размер области L = 1.

На рисунке 2 представлены профили скорости в центральном сечении полости (x = 0.5L, y = 0.5L) для различных значений параметра микрополярности N. При N = 0 профиль практически линейный, что соответствует течению Куэтта. С увеличением N профиль становится нелинейным: вблизи стенок градиент скорости возрастает, в центральной части уменьшается. Это обусловлено влиянием микровращения на перенос количества движения [2; 31]. При N=0.9 отклонение от линейного профиля достигает 18 % на расстоянии z/L=0.5.

 

Рисунок 2. Профили скорости в центральном сечении полости (x = 0.5L, y = 0.5L) для различных значений параметра микрополярности N

 

На рисунке 3 показано распределение модуля микровращения $\mid \omega \mid =\sqrt{{\omega }_x^2+{\omega }_y^2+{\omega }_z^2}$ в сечении z = 0.5L. Максимальные значения микровращения наблюдаются вблизи верхней и нижней стенок, а также в угловых зонах, где градиенты скорости наибольшие. В центральной области канала микровращение минимально. Такое распределение соответствует физическому смыслу микровращения как меры вращения микроэлементов жидкости под действием градиентов скорости [27; 33]. Изолинии имеют регулярную структуру, что свидетельствует об отсутствии численных осцилляций в решении.

 

Рисунок 3. Распределение модуля микровращения |ω| в сечении z = 0.5L

 

На рисунке 4 представлено распределение диссипативной функции Φ в том же сечении. Диссипация максимальна вблизи движущейся верхней стенки и в угловых зонах, где присутствуют значительные градиенты скорости и микровращения [34; 35]. Распределение диссипации несимметрично относительно центра канала, что связано с асимметрией граничных условий (верхняя стенка движется, нижняя неподвижна). Отсутствие отрицательных значений диссипации подтверждает физическую корректность численного решения.

 

Рисунок 4. Распределение диссипативной функции Φ в сечении z = 0.5L

 

Для количественной оценки влияния параметра N вычислены интегральные величины. При увеличении N от 0 до 0.9 максимальная скорость микровращения возрастает от 0 до 3.8×10^-3, а интегральная диссипация увеличивается в 2.8 раза, что указывает на существенное влияние микрополярных эффектов на энергетические потери [25; 38].

Исследована зависимость решения от шага по времени. При уменьшении Δt от 10^-4 до 10^-5 изменение максимальной скорости не превышает 0.01 %. Выбранный шаг Δt = 2×10^-5 обеспечивает устойчивость и достаточную точность [1].

Проведён анализ чувствительности к виду аппроксимации конвективных членов. Центральные разности приводят к численным осцилляциям при Re > 10, тогда как схема против потока даёт монотонное решение во всём диапазоне параметров [9; 16].

Граничные условия выполняются с высокой точностью: на верхней стенке скорость u равна U₀ с точностью до 10^-6, на нижней и боковых стенках условия прилипания – с точностью до 10^-7; для микрополярной скорости нормальные производные на границах близки к нулю с точностью 10^-5.

Метод устойчив при варьировании параметров в широком диапазоне: Re ∈ [1, 50], N ∈ [0,0.9], m ∈ [0.1,0.5]. Время расчёта на сетке 31³ составляет 5–10 минут на процессоре Intel Core i7, что приемлемо для параметрических исследований [27; 39].

Разработанный численный метод имеет второй порядок точности по пространству и первый – по времени. Основной вклад в погрешность вносят аппроксимация конвективных членов схемой против потока и явная схема по времени [16; 22]. Для стационарных течений это не является критическим, так как решение сходится к стационарному состоянию.

Ограничением метода является требование условия Куранта – Фридрихса – Леви, что при больших числах Рейнольдса требует малого шага по времени [1]. Для Re > 100 рекомендуется применение неявных или полунеявных схем [8; 12; 24]. Также метод использует равномерную сетку, что может быть неоптимально для задач с большими градиентами в пристеночных областях. В дальнейшем планируется адаптация алгоритма для неравномерных сеток со сгущением к стенкам [9].

Сравнение с результатами других авторов для двумерных течений микрополярных жидкостей показывает хорошее соответствие [3; 20]. Расхождения не превышают 3 % для скоростей и 5 % для микровращения, что находится в пределах точности различных численных методов.

ЗАКЛЮЧЕНИЕ И ВЫВОДЫ

Разработан и верифицирован конечно-разностный метод для моделирования трёхмерного течения микрополярной жидкости в кубической полости с подвижной верхней стенкой. Метод основан на комбинации схемы против потока для конвективных членов и метода проекций для решения связанной системы уравнений Навье – Стокса и микровращения.

Верификация показала второй порядок сходимости по пространству, высокую точность (относительная погрешность не превышает 2.4×10^-4 при сравнении с аналитическим решением для ньютоновской жидкости) и выполнение законов сохранения с высокой точностью.

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

Метод устойчив в широких диапазонах чисел Рейнольдса и параметров микрополярности, что позволяет использовать его для параметрических исследований микроструктурированных жидкостей в каналах конечных размеров.

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

About the authors

Evgenii Y. Prosviryakov

Ural Federal University; Institute of Engineering Science, Ural Branch of the RAS

Email: evgen_pros@mail.ru

Doctor of Physics and Mathematics, Associate Professor, Professor at the Department of Information Technology and Automation, Head of Sector, Sect. of Nonlinear Vortex Hydrodynamics

Russian Federation, Ekaterinburg; Ekaterinburg

Kristina V. Gubareva

Samara State Technical University

Author for correspondence.
Email: r.kristina2017@mail.ru

Candidate of Engineering Science, Associate Professor at the Department of Industrial Thermal Power Engineering

Russian Federation, Samara

References

Supplementary files