Аналитическое решение уравнений Навье – Стокса для описания неоднородного течения Куэтта с квадратичным профилем в слое с проницаемыми границами

Введение. Точное управление структурой течения в тонких слоях критически важно для современных технологий — от микрофлюидных устройств до систем смазки. Однако в подобных случаях некоторые классические гидродинамические модели не позволяют учесть комплекс технологических параметров: широкий диапазон вязкостей рабочих жидкостей, проницаемость стенок, а также их пространственную неоднородность (шероховатость, распределенные источники) [1]. Пример такой «недостаточной» модели — течение Куэтта между непроницаемыми пластинами [2]. Представленное исследование актуально в плане развития перечисленных ниже прикладных направлений.

1. Биомедицинская микрофлюидика и так называемые лаборатории на чипе. Эффективное смешение реагентов или адресная доставка клеток предполагают, во-первых, инжекцию жидкости через пористые мембраны (V_w≠ 0). Во-вторых, необходимы неоднородные сдвиговые течения, формируемые за счет микрорельефа стенок (что моделируется граничными коэффициентами A и B) [3].
2. Энергоэффективные мембранные технологии (опреснение, разделение газов). В этой сфере производительность определяется взаимодействием продольного потока и поперечной фильтрации через проницаемую стенку, а также вязкостными свойствами разделяемой среды [4].
3. Трибология пористых и текстурированных подшипников. Распределение давления, трение и износ напрямую зависят от течения смазочного материала в микроскопическом зазоре со сложным рельефом (A, B) и возможной фильтрацией (V_w) через пористый вкладыш [5].

Анализ литературы позволил систематизировать подходы к описанию указанных факторов и выявил существенный пробел. И для модельных [6], и для прикладных задач [7] детально исследовались эффекты проницаемости границ, определяющие нелинейные (экспоненциальные) профили скорости. Хорошо описано влияние сложных (в том числе полиномиальных) граничных условий на непроницаемых стенках, что связано с задачами устойчивости, тепломассопереноса и моделирования шероховатости [8]. Глубоко изучено влияние реологических свойств, включая широкий диапазон вязкостей реальных жидкостей [9]. Для решения сопряженных задач активно применяются современные численные [10] и аналитические методы [11]. Однако перечисленные выше факторы традиционно изучались вне системного подхода. Аналитические решения либо описывают проницаемость при однородных граничных условиях [12], либо учитывают нелинейность границ только для непроницаемых стенок [13]. Таким образом, в литературе не представлено комплексное аналитическое решение, которое в явном виде объединяет ключевые безразмерные параметры:

• динамическую вязкость (μ);
• проницаемость (α = ρV_w/μ);
• линейную (A) и квадратичную (B) поправки к граничному профилю скорости.

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

Цель работы — создание и анализ точного решения для обобщенного стационарного течения Куэтта ньютоновской жидкости между проницаемыми пластинами с квадратичным профилем скорости на границе. Научная новизна заключается в получении замкнутого аналитического решения, впервые явно и комплексно учитывающего синергию параметров α, A, B и μ. Основное преимущество такого подхода перед численным моделированием — возможность мгновенного получения решения и проведения прямого анализа физических зависимостей, не опосредованного сеточными аппроксимациями.

Для достижения цели решаются четыре задачи.

1. Вывод и строгая верификация аналитического решения, включая проверку предельных переходов к известным частным случаям.
2. Исследование влияния знака и величины граничных коэффициентов A и B на пространственное распределение скорости и завихренности.
3. Количественный анализ влияния динамической вязкости реальных жидкостей (вода, моторное масло, нефть) на профиль скорости и толщину сдвигового слоя.
4. Обсуждение прикладной значимости результатов, ограничений модели и перспективных направлений ее обобщения.

Материалы и методы. Рассмотрим установившееся (стационарное), ламинарное, изотермическое течение ньютоновской вязкой и несжимаемой жидкости в канале, образованном двумя бесконечными, параллельными, плоскими пластинами (рис. 1).

Расстояние между пластинами постоянно и равно h. Нижняя пластина неподвижна и расположена в плоскости z = 0. Верхняя расположена в плоскости z = h и движется с постоянной скоростью W в положительном направлении оси x.

Ключевая особенность данной задачи: проницаемость обеих пластин обеспечивает фильтрацию жидкости в направлении, нормальном к их поверхностям [6]. Предполагается, что эта нормальная скорость (V_w) постоянна по всей поверхности каждой пластины и одинакова для обеих пластин. Если V_w > 0, жидкость втекает в канал через нижнюю пластину и вытекает через верхнюю; если V_w < 0, направление потока через стенки обратное.

Продольная компонента скорости V_x задается в виде квадратичного разложения по поперечной координате y. Такое обобщение классического профиля Куэтта позволяет учесть более сложные граничные условия на верхней пластине [5]:

(1)

Здесь функции U(z), u₁(z) и u₂(z) представляют собой соответственно основную (осевую), линейную и квадратичную моды профиля скорости, зависящие только от координаты z.

Граничные условия формулируются на основе условия прилипания. На неподвижной нижней пластине (z = 0) задается полное отсутствие движения:

(2)

На верхней пластине (z = h) задается обобщенный профиль скорости:

(3)

где W — скорость поступательного движения пластины; A — коэффициент, определяющий линейный градиент скорости вдоль y; B — коэффициент, определяющий кривизну профиля:

Неизвестные функции U(z), u₁(z) и u₂(z) можно найти по уравнению Навье – Стокса для x‑компоненты скорости. В предположении стационарности оно имеет вид [2]:

(4)

Рассматриваемое течение — обобщение чистого течения Куэтта, для которого характерно отсутствие градиента давления вдоль канала (∂p/∂x = 0). Учтем это условие и то, что V_x не зависит от x, V_y = 0 и V_z = V_w = const [3]. В данном случае (4) значительно упрощается:

(5)

Подстановка (1) в (5) и вычисление соответствующих производных дает уравнение, которое должно выполняться для всех значений y. Приравняем коэффициенты при одинаковых степенях y и получим систему трех независимых обыкновенных дифференциальных уравнений (ОДУ) [11]:

(6)

(7)

(8)

Выражения (6) и (7) — линейные однородные ОДУ второго порядка с постоянными коэффициентами. Для их решения вводится безразмерный параметр проницаемости α = ρV_w / μ, который физически представляет собой отношение инерционных сил к вязким силам в нормальном направлении. Разделив (6) и (7) на μ, получаем каноническую форму:

(9)

Общие решения этих уравнений имеют вид:

(10)

Применение граничных условий (2) и (3) позволяет определить константы интегрирования. Для функции u₁(z) условия u₁(0) = 0 и u₁(h) = A дают:

(11)

Аналогично для u₂(z) из условий u₂(0) = 0 и u₂(h) = B находим:

(12)

Окончательные выражения для линейной и квадратичной мод профиля скорости:

(13)

(14)

Уравнение (8) для основной моды U(z) является неоднородным. Подставим в него найденное выражение для u₂(z) (14) и разделим на μ:

(15)

Для уравнения (15) применяется метод неопределенных коэффициентов [13]. Его общее решение — это сумма общего решения однородного уравнения и частного решения неоднородного уравнения. После определения констант интегрирования из граничных условий (2) и (3) окончательное выражение для основной моды принимает вид [14]:

(16)

Таким образом, подстановка (13), (14) и (16) в исходный анзац (1) дает определение полного поля продольной компоненты скорости. Полученное аналитическое решение удовлетворяет как уравнению движения (5), так и всем граничным условиям, т. е. (2) и (3).

Для подтверждения физической состоятельности решения рассмотрим предельный переход V_w → 0, что эквивалентно α → 0. В этом пределе, используя разложение экспоненты в ряд Тейлора, получаем для квадратичной моды:

(17)

Аналогично для линейной моды:

(18)

Для основной моды:

(19)

Численное моделирование. На данном этапе выполняются:

• верификация и анализ полученного решения;
• исследование влияния физических параметров на структуру течения.

Для этого авторы представленной статьи выполнили комплексное численное моделирование в среде Matlab. Основные цели программной реализации:

• визуальное представление полей скорости;
• количественное сравнение характеристик течения для различных типов жидкостей и граничных условий.

Центральный элемент кода — функция solve_couette_full¹, которая реализует формулы (13), (14), (16).

На входе будут:

• физические свойства жидкости (плотность ρ и динамическая вязкость μ);
• геометрические параметры задачи (расстояние между пластинами h);
• кинематические характеристики (W, A, B, V_w);
• координатная сетка по z.

На выходе — базовые функции U(z), u₁(z), u₂(z) и их производные по z для последующего расчета таких важнейших гидродинамических величин, как напряжение сдвига и завихренность.

Моделирование проходило в два этапа. На первом выполнили базовый расчет для стандартного набора параметров. Это позволило проверить выполнение граничных условий. Максимальное отклонение скорости от заданных значений на границах z = 0 и z = h не превысило 10^–12, что подтвердило высокую точность аналитического подхода и корректность программной реализации.

На втором этапе код расширили для анализа двух ключевых аспектов. Во-первых, детально исследовали влияние знаков и величин коэффициентов A и B в граничном профиле скорости на верхней пластине. Рассматривались четыре комбинации: A > 0, B > 0; A > 0, B < 0; A < 0, B > 0; A < 0, B < 0. Это позволило смоделировать несколько физических сценариев (от ускорения потока к центру канала до его торможения у стенок), а также переход от выпуклых к вогнутым профилям скорости. Для каждой комбинации строились двумерные карты распределения скорости V_x(y, z) и профили скорости вдоль нормали к пластинам при фиксированных значениях поперечной координаты y. Это дало наглядное представление о том, как изменение знака квадратичного члена B кардинально меняет форму изолиний скорости, а знак линейного коэффициента A определяет асимметрию профиля относительно плоскости y = 0.

Во-вторых, рассчитали течение для различных жидкостей. В модель заложили физические свойства воды, моторного масла SAE 30² и сырой нефти при температурах 20 °C, 40 °C или 60 °C [14]. Использование справочных данных для динамической вязкости и плотности позволило перейти от абстрактной математической модели к практически значимым инженерным расчетам. Для каждой жидкости автоматически определялось число Рейнольдса Re = Wh / ν, что дало возможность оценить режим течения. Как и ожидалось, для воды при комнатной температуре число Рейнольдса достигало Re ≈ 1000. При этом течение оставалось ламинарным, что согласуется с известной устойчивостью классического течения Куэтта. Для масла и нефти течение оставалось глубоко ламинарным (Re < 20).

Особое внимание уделили расчету производных величин. Для каждого значения координаты y аналитически вычислили напряжение сдвига τ_xz = μ∂V_x/∂z и завихренность ω_y = −∂V_x/∂z. Это позволило избежать ошибок численного дифференцирования [15] и обеспечило высокую точность при построении распределений указанных величин как вдоль зазора между пластинами, так и в зависимости от поперечной координаты y.

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

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

Функции u₁(z), u₂(z), U(z) параметрически исследовали методом численного моделирования в Matlab на основе аналитического решения, выведенного из (13), (14) и (16). Расчеты выполнены при фиксированных параметрах (рис. 2–5): h = 1,0 м, W = 0,3 м/с, V_w = 0,001 м/с, μ = 0,01 Па·с для четырех комбинаций знаков: A = ±0,03 с^–1 и B = ±0,005 м^–1·с^–1.

Для реальных жидкостей расчеты выполнялись при h = 0,02 м, W = 0,05 м/с, A = 0,1 с^–1,
B = 0,02 м^–1·с^–1, V_w = 0,0005 м/с. Использовались: вода при 20 °C (μ = 1,002·10^–3 Па·с) и 40 °C (μ = 0,653·10^–3 Па·с); нефть при 20 °C (μ = 9,15·10^–3 Па·с) и 40 °C (μ = 4,72·10^–3 Па·с); масло SAE 30 при 20 °C (μ = 0,290 Па·с) и 60 °C (μ = 0,050 Па·с).

Результаты исследования. На рис. 2 представлены поля завихренности ω_y. При A > 0 максимальные положительные значения ω_y находятся в области положительных координат y. При A < 0 максимумы ω_y смещены в область y < 0. Амплитуда значений ω_y варьируется с изменением коэффициента B.

На рис. 3 представлены поле модуля скорости |V| и векторное поле для двух комбинаций знаков параметров A и B. При A > 0 векторы длиннее в области y > 0, что указывает на смещение максимума скорости в положительную область поперечной координаты. При A < 0 наблюдается обратная картина: векторы длиннее в области y < 0, что свидетельствует о смещении максимума в отрицательную область. Это проявление асимметрии течения, обусловленной линейным членом A в граничном условии. Форма изолиний |V| различается для случаев B > 0 и B < 0, что связано с влиянием квадратичного члена B.

На рис. 4 представлены изолинии продольной скорости V_x и линии тока. Линии тока отклоняются от направления оси канала. Направление отклонения определяется знаком A. Изолинии V_x при B > 0 и B < 0 имеют различную кривизну.

На рис. 5 показаны профили V_x(y) на разных уровнях, т. е. при z, равном 0, h, h/4, h/2 и 3h/4. При z = 0 скорость V_x равна нулю. С ростом z вариация V_x по y увеличивается и достигает максимума при z = h, где профиль определяется параметрами A и B. При z = h/2 профиль близок к линейному.

Расчеты выполнены для реальных жидкостей при h = 0,02 м, W = 0,05 м/с, A = 0,1 с^–1,
B = 0,02 м^–1·с^–1, V_w = 0,0005 м/с (рис. 6–8). Использовались: вода при 20 °C (μ = 1,002·10^–3 Па·с) и 40 °C (μ = 0,653·10^–3 Па·с); нефть при 20 °C (μ = 9,15·10^–3 Па·с) и 40 °C (μ = 4,72·10^–3 Па·с); масло SAE 30 при 20 °C (μ = 0,290 Па·с) и 60 °C (μ = 0,050 Па·с).

На рис. 6 приведены профили V_x(z) при y = 0. Все профили построены в диапазоне z ∈ [0, 20] мм, соответствующем толщине слоя.

На рис. 7 отображены изолинии |V| и линии тока. Для воды линии тока минимально отклоняются от горизонтали, а для масла при 20 °C — значительно искривляются вблизи верхней стенки.

На рис. 8 представлены двумерные поля завихренности ω_y для различных жидкостей.

Обсуждение. Согласно данным рис. 6, для воды при температуре 20 °C характерно слабое отклонение профиля V_x(z) от линейного закона, что связано с ее малой вязкостью. При повышении температуры до 40 °C это отклонение становится еще менее выраженным, а общий наклон профиля несколько сокращается. Для масла SAE 30 (20 °C) наблюдается умеренная нелинейность профиля, которая при нагреве до 60 °C значительно ослабевает. В данном случает профиль приближается к линейному, что соответствует снижению вязкости при повышении температуры. Для нефти (20 °C) профиль характеризуется выраженной нелинейностью. При увеличении температуры до 40 °C кривая становится ближе к линейной, что также согласуется с уменьшением вязкости.

Как показывает рис. 8, для воды при 20 °C и 40 °C завихренность отрицательная и изменяется в диапазоне от –35 до –5 с^–1, что соответствует слабой завихренности по всему каналу. Для масла SAE 30 при 20 °C наблюдается локализованная область повышенной по модулю завихренности
(|ω_y| > 2,3 с^–1), сосредоточенная вблизи верхней стенки. При нагреве до 60 °C область повышенной завихренности расширяется, а ее значения по модулю снижаются, что согласуется с уменьшением вязкости. Для нефти при 20 °C завихренность достигает значений до –4,5 с^–1, а при 40 °C — до –1,5 с^–1, что также соответствует снижению вязкости.

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

Определяющее влияние знака коэффициента A на асимметрию течения (рис. 2–5) объясняется его физическим смыслом как градиента скорости ∂V_x/∂y на границе. Этот градиент, заданный на верхней стенке, за счет вязкой диффузии формирует поперечную составляющую потока во всем объеме, что приводит к смещению линий тока и максимумов завихренности. Коэффициент B, управляющий кривизной профиля (∂²V_x/∂y²), определяет распределение этого поперечного потока по ширине канала, что проявляется в изменении формы изолиний скорости (рис. 3, 4). Ослабление влияния граничных коэффициентов с удалением от верхней стенки (рис. 5) количественно описывается экспоненциальными зависимостями (13) и (14), что согласуется с представлениями о затухании возмущений от границы в вязкой жидкости [16].

Безразмерный параметр проницаемости α зависит от динамической вязкости μ (α = ρV_w / μ). Этим объясняется резкое различие между течением воды и масла (рис. 6, 8).

Для высоковязких жидкостей (малые показатели α) решение вдали от границ стремится к линейному профилю. Однако неоднородность в уравнении (8) для U(z), обусловленная членом с B, приводит к концентрации сдвига в тонком пристенном слое, толщина которого обратно пропорциональна |α|. Таким образом, вязкость — не просто множитель, а параметр, управляющий пространственной локализацией сдвиговой деформации. Это имеет важное практическое значение, например, для расчета трения в подшипниках с пористой смазкой [17].

Полученное решение верифицировано, т. к. корректно воспроизводит известные предельные случаи. При α→0 оно переходит в решение для течения с квадратичным условием на непроницаемой стенке [18]. При A, B → 0 сводится к классическому экспоненциальному профилю для проницаемых стенок [19]. Одновременный переход α, A, B → 0 дает линейный профиль классического течения Куэтта [2].

Явный вид решения позволяет напрямую оценивать влияние каждого параметра на поле течения, что ценно для инженерного проектирования. Например, в микрофлюидном микшере можно варьировать микрорельеф стенки (через A и B) для генерации заданной вихревой структуры, улучшающей перемешивание [20]. В задаче смазки пористого подшипника модель позволяет аналитически связать вязкость масла, скорость просачивания и шероховатость поверхности (через A, B) с распределением напряжения сдвига и диссипацией энергии [21].

Основное ограничение модели — допущения о стационарности, ламинарности и неньютоновском поведении жидкости. Расширение модели для учета турбулентности выходит за рамки данной работы. Это отдельная сложная задача, требующая перехода к осредненным уравнениям Рейнольдса (RANS)³ или моделям крупных вихрей (LES)⁴ [20]. Еще одно перспективное направление — обобщение для неньютоновских жидкостей [21] и нестационарных режимов [22]. Это позволит охватить более широкий класс таких прикладных задач, как пульсирующие течения в биомедицинских устройствах.

Заключение. Впервые получено точное аналитическое решение стационарных уравнений Навье – Стокса, описывающее обобщенное течение Куэтта ньютоновской жидкости между проницаемыми пластинами с квадратичным профилем скорости на границе.

Параметрический анализ выявил, что линейный коэффициент A граничного условия определяет направление асимметрии полей скорости и завихренности, а квадратичный коэффициент B — степень их пространственной нелинейности.

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

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

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