A large-scale layered stationary convection of an incompressible viscous fluid under the action of shear stresses at the upper boundary. Velocity field investigation

Full Text

Введение. Известно, что в потоке жидкости могут существовать точки, где нулевое значение принимает одна или несколько компонент вектора скорости. Такие точки называются застойными, или критическими [1]. Рядом с такими точками, как правило, образуются застойные зоны, то есть области с обратным течением. Данное утверждение справедливо для решений уравнений Навье-Стокса и их модификаций, являющихся непрерывными функциями координат. Другие точные решения, обладающие застойными точками, приведены в [2-4]. Противотечения возникают, например, в следе за цилиндром при небольших числах Рейнольдса и Фруда в жидкости с небольшим периодом плавучести [5]. Экспериментальные исследования показали, что противотечения наблюдаются в поверхностных, подповерхностных и глубинных слоях практически во всех районах Мирового океана [6-9]. Проблема учета и исследования застойных зон оказывается актуальной не только в задачах, связанных с обтеканием движущихся в жидкой или газообразной среде тел применительно, в частности, к судам и летательным аппаратам. Она возникает, например, при использовании скважин в нефтедобыче [10], когда необходимо определить характер течения в окрестности точки изменения направления потока. Данная проблема также имеет место при обтекании ветром зданий и является определяющей для расчета выбросов вредных веществ техногенных систем. Теоретическое изучение этих и многих других процессов, связанных с исследованием противотечений, в настоящее время далеко от завершения. В связи с этим отыскание новых точных решений уравнений гидродинамики, описывающих появление застойных зон в жидкости, является актуальной задачей. Первым изотермическим точным решением, посвященным исследованию застойных точек, является течение Хименца [11]. Исследование противотечений в океане было начато Экманом [12]. В [12] было представлено точное решение, описывающее течение вязкой несжимаемой жидкости при балансе вязких и кориолисовых сил. Давление в этом случае полагалось однородным по горизонтальным координатам. На основе течений Экмана [12], Чарни [13] и Стоммела [14] были предложены различные математические модели описания экваториальных противотечений. Дальнейшие исследования в этой области опирались на вычислительный эксперимент. Отметим, что в [15-20] построены широкие классы точных решений для описания жидкостей при вращении. Несмотря на большое число работ, в которых учитывается планетарная завихренность, особый интерес представляют течения жидкости в отсутствие поля сил Кориолиса. В этом случае для описания течений можно использовать «бесконечно замкнутый» горизонтальный слой. В [21-23] приведены точные решения крупномасштабных изотермических течений при постоянном давлении и конвективных движений различной природы. В статьях [24-35] приведены точные решения для конвективных термокапиллярных течений, Б у р м а ш е в а Н. В., П р о с в и р я к о в Е. Ю. которые являются обобщениями течения Остроумова-Бириха [36, 37]. В данной работе мы рассматриваем движение вязкой несжимаемой жидкости под действием заданных на верхней границе тангенциальных сил, не являющихся результатом капиллярного эффекта. 1. Постановка задачи. Рассмотрим задачу описания конвективных движений в слоистых крупномасштабных течениях вязкой несжимаемой жидкости в стационарном случае. Наиболее распространенным подходом в таких задачах является использование уравнений тепловой конвекции в приближении Буссинеска, предполагающего аппроксимацию плотности как линейную функцию температуры. С математической точки зрения данная модель [38] состоит из - уравнения Навье-Стокса (V · ∇)V = -∇P + ν V + gβT k, (1) - уравнения теплопроводности V · ∇T = χ T, (2) ∇ · V = 0. (3) - уравнения несжимаемости Здесь V (x, y, z) = {Vx , Vy , Vz } - вектор скорости, причем Vz ≡ 0, так как рассматривается слоистое течение; P - отклонение давления от гидростатического, деленное на постоянную среднюю плотность ρ жидкости; T - отклонение от средней температуры; ν, χ - коэффициенты кинематической вязкости и температуропроводности жидкости соответственно; k - орт оси z, направленной вертикально вверх; ∇ - оператор Гамильтона; - двумерный оператор Лапласа. Заметим, что данная система суть система из пяти скалярных уравнений относительно четырех неизвестных (компонент скорости Vx , Vy и физических полей T и P ). То есть система переопределена. Следовательно, чтобы решать такую систему, надо либо убедиться в непротиворечивости системы условий, либо искать такое решение, на котором одно из имеющихся условий удовлетворится тождественно. Кроме того, отметим, что члены конвективной производной из уравнения Навье-Стокса и уравнения несжимаемости тождественно равны нулю, но они сохраняются в уравнении теплопроводности. Поэтому, строго говоря, решается нелинейная задача. 2. Условие разрешимости переопределенной системы. Для того чтобы разрешить переопределенную исходную систему (1)-(3) предлагается пойти вторым из предложенных выше путей и искать решение в классе, предложенном в [20, 24]: Vx = u(z), Vy = v(z). В этом случае уравнение неразрывности (3) обращается в тождество: ∂Vx ∂Vy ∂Vz ∂u(z) ∂v(z) + + = + ≡ 0. ∂x ∂y ∂z ∂x ∂y 182 (4) Крупномасштабная слоистая стационарная конвекция вязкой несжимаемой жидкости. . . Далее будем искать давление и температуру в виде специальных линейных форм по горизонтальным координатам: T = T0 (z) + T1 (z)x + T2 (z)y, P = P0 (z) + P1 (z)x + P2 (z)y. (5) Подставим выбранный класс (4), (5) в уравнение Навье-Стокса (1): Vx u(z) ∂Vy ∂Vx ∂Vz ∂P + Vy + Vz =- + ν Vx , ∂x ∂y ∂z ∂x ∂v(z) ∂Vz ∂P ∂u(z) + v(z) +0 =- + ν Vx , ∂x ∂y ∂z ∂x ∂2 ∂2 ∂2 ∂ P0 (z) + P1 (z)x + P2 (z)y = ν + + u(z), ∂x ∂x2 ∂y 2 ∂z 2 откуда согласно методу неопределенных коэффициентов получаем P1 = νu . Здесь и всюду далее штрихом обозначено дифференцирование по переменной z. Аналогично из второго уравнения Навье-Стокса (относительно компоненты Vy ) имеем P2 = νv . Из третьего уравнения сохранения импульсов в жидкости имеем Vx ∂Vy ∂Vx ∂Vz ∂P + Vy + Vz =- + ν Vz + gβT, ∂x ∂y ∂z ∂z u(z) ∂u(z) ∂v(z) ∂Vz ∂P + v(z) +0 =- + gβT, ∂x ∂y ∂z ∂z ∂P = gβT ∂z или ∂ P0 (z) + P1 (z)x + P2 (z)y = gβ T0 (z) + T1 (z)x + T2 (z)y , ∂z откуда согласно методу неопределенных коэффициентов получаем три скалярных соотношения: ∂Pi = gβTi , i = 0, 1, 2. ∂z И, наконец, подставим выбранный класс решений (4), (5) в уравнение теплопроводности (2): Vx u ∂T ∂T ∂T ∂2 ∂2 ∂2 + Vy + Vz =χ + + T, ∂x ∂y ∂z ∂x2 ∂y 2 ∂z 2 ∂T0 (z) + T1 (z)x + T2 (z)y ∂T0 (z) + T1 (z)x + T2 (z)y +v = ∂x ∂y 183 Б у р м а ш е в а Н. В., П р о с в и р я к о в Е. Ю. =χ ∂2 ∂2 ∂2 + + ∂x2 ∂y 2 ∂z 2 T0 (z) + T1 (z)x + T2 (z)y , uT1 + vT2 = χT0 (z) + χT1 (z)x + χT2 (z)y. Отсюда согласно методу неопределенных коэффициентов имеем uT1 + vT2 = χT0 (z), T1 (z) = T2 (z) = 0. Таким образом, окончательно получаем нелинейную систему уравнений: T1 = 0, ∂P1 = gβT1 , ∂z νu = P1 , T2 = 0. ∂P2 = gβT2 ∂z νv = P2 , χT0 = uT1 + vT2 , ∂P0 = gβT0 , ∂z (6) (7) (8) (9) (10) 3. Краевая задача. Система обыкновенных дифференциальных уравнений (6)-(10) является системой тринадцатого порядка. Следовательно, чтобы корректно определить краевую задачу, необходимо иметь тринадцать условий, заданных на границе слоя, которые будут использованы для определения констант, возникающих в ходе интегрирования рассматриваемой системы. Будем далее считать, что нижняя граница слоя жидкости является абсолютно твердой и неподвижной. Верхнюю границу будем полагать свободной. В качестве краевых условий рассмотрим следующие: - пусть на нижней границе (z = 0) слоя жидкости выполняется условие прилипания, а температура задается функцией T = Ax + By; - на верхней границе (z = h) действует постоянное давление S, а температура задается функцией T = ϑ + Cx + Dy; - кроме того, на свободной поверхности z = h заданы напряжения: η du = ξ1 , dz η dv = ξ2 . dz Таким образом, приходим к следующей системе граничных условий: T0 (0) = 0, T0 (h) = ϑ, P0 (h) = S, 184 u(0) = v(0) = 0, T1 (0) = A, T1 (h) = C, P1 (h) = 0, du η (h) = ξ1 , dz T2 (0) = B, T2 (h) = D, P2 (h) = 0, Крупномасштабная слоистая стационарная конвекция вязкой несжимаемой жидкости. . . η dv (h) = ξ2 . dz Без ограничения общности можно считать S = 0, тем самым ведя отсчет приведенного давления от уровня, задаваемого на верхней границе. 4. Решение системы уравнений. Рассмотрим частный случай поставленной выше краевой задачи - задание температурного возмущения только на верхней границе (A = B = 0). Начнем последовательно интегрировать основную систему уравнений с уравнения (6). С учетом соответствующих граничных условий получаем z T1 = C , h z T2 = D . h Далее полученные выражения подставляем в уравнения (7) и находим градиенты P1 , P2 приведенного давления: P1 = gβC z 2 - h2 , 2h P2 = gβD z 2 - h2 . 2h Затем из уравнений (8) находим компоненты вектора скорости: u(z) = z gβD 8ξ1 (2h2 z - z 3 ) + , 8 hν η v(z) = z gβC 8ξ2 (2h2 z - z 3 ) + . 8 hν η Далее найденные выражения компонент скорости и градиентов температуры подставляем в уравнение (9), откуда T0 (z) = z(-h3 (Cξ1 + Dξ2 ) + 12χηϑ + (Cξ1 + Dξ2 )z 3 ) + 12χηh CDgβz(-16h6 + 21h2 z 4 - 5z 6 ) + , 840χh2 ν и, наконец, из уравнения (10) определим составляющую P0 давления: P0 (z) = βg[3h5 (Cξ1 + Dξ2 ) - 60χηh2 ϑ - 5(h3 (Cξ1 + Dξ2 ) - 12χηϑ)z 2 ] + 120χηh 2βg(Cξ1 + Dξ2 )z 5 β 2 CDg 2 (41h8 - 64h6 z 2 + 28h2 z 6 - 5z 8 ) + + . 120χηh 6720χh2 ν 5. Переход к термокапиллярной конвекции. Рассматриваемая краевая задача может быть сведена к задаче исследования термокапиллярной конвекции (конвекции Бенара-Марангони), описанной в [24]. Для этого достаточно уточнить условия на нижней (твердой) и верхней (свободной) поверхностях, положив T0 (h) = 0, η du (h) = ξ1 = -σT1 (h), dz η dv (h) = ξ2 = -σT2 (h). dz 185 Б у р м а ш е в а Н. В., П р о с в и р я к о в Е. Ю. 6. Исследование компонент скорости. Функции компонент скорости u, v имеют ноль в начале координат (множитель z выделяется явно). Изучим вопрос о том, сколько еще нулей может быть у этих функций. Этот вопрос является довольно принципиальным, поскольку наличие нулей у этих функций говорит о существовании противотечений (расслоении) в конвективном потоке вязкой несжимаемой жидкости. Пронормируем найденные в п. 4 выражения компонент u(z), v(z) вектора скорости V на величину C, приведя их к безразмерному виду. Для этого введем новые параметры системы: ∗ = D/C, δ = h/l, где h - характерный вертикальный размер слоя, а l - характерный горизонтальный размер слоя. Кроме этого, перейдем к безразмерной координате z → z/h. Для этого разделим u(z), v(z) на gβl3 /ν. В результате получим u(z) = zδ 3 8 ∗ (2z - z 3 ) + 8W1 , v(z) = zδ 3 2z - z 3 + 8W2 . 8 Здесь Wi = νξi gβCηh2 - числа Вебера, посчитанные для значений ξi , i = 1, 2. Очевидно, что W1 /W2 = ξ1 /ξ2 . Отметим также, что после обезразмеривания координата z теперь меняется в диапазоне [0, 1]. Рассмотрим сначала функцию u(z). Случай ∗ = 0 опустим из-за тривиальности. Поэтому далее представим функцию u(z) в мультипликативном виде: 8W1 z ∗δ3 2z - z 3 + ∗ u(z) = 8 и введем в рассмотрение функцию u1 (z) = 2z - z 3 + 8W1 ∗ . Эта функция является полиномом третьей степени относительно z. Чтобы определить, сколько нулей у функции u1 (z), рассмотрим вспомогательную функцию f (z) = 2z - z 3 , график которой представлен на рис. 1. Отметим, у функции f (z) на [0, 1] есть всего один ноль - точка z0 = 0. Кроме этого, глобальный максимум функции f на интересующем нас отрезке [0, 1] есть точка z∗ = 2/3, значение функции в которой равно √ 4 2 fmax = f (z∗ ) = √ > 1. 3 3 А на правой границе z = 1 она принимает значение f1 = f (1) = 1. 186 Крупномасштабная слоистая стационарная конвекция вязкой несжимаемой жидкости. . . Рис. 1. График функции f (z) [Figure 1. The graph of f (z)] Производя параллельный перенос графика этой функции вдоль вертикальной оси, то есть рассматривая семейство функций вида fc (z) = f (z) + c, можно добиться того, что функция fc (z) пересечет ось z на отрезке [0, 1] один или два раза, или же вообще не будет пересекать ось z (не имея при этом общих точек или только единственную общую точку - точку касания). Различные положения функции fc представлены на рис. 2. Таким образом, возвращаясь к функции u1 , получаем следующие оценки: 1) - 8W1 ∗ 2) 1 = - < 1 - единственный ноль, одна смена знака; 8W1 - два различных нуля, одна смена знака; √ 8W1 4 2 3) 1 < - ∗ < √ - два различных нуля, две смены знака; 3 3 √ 4 2 8W1 4) √ = - ∗ - один ноль, функция u1 (z) не меняет знак; 3 3 ∗ Рис. 2. Семейство функций fc (z) = f (z) + c для c = -f (1) + 0.05 (линия 1), для c = -(f (z∗ ) + f (1))/2 (линия 2), для c = -f (z∗ ) (линия 3), для c = -f (z∗ ) - 0.05 (линия 4) [Figure 2. The family of functions fc (z) = f (z) + c for c = -f (1) + 0.05 (line 1), for c = -(f (z∗ ) + f (1))/2 (line 2), for c = -f (z∗ ) (line 3), for c = -f (z∗ ) - 0.05 (line 4)] 187 Б у р м а ш е в а Н. В., П р о с в и р я к о в Е. Ю. √ 4 2 8W1 5) √ < - ∗ - нулей нет, функция u1 (z) не меняет знак. 3 3 В соответствии с каждым из перечисленных случаев для функции u(z) можем получить следующие картины (без ограничения общности взято δ = 1, ∗ = 1) (см. рис. 3). ∗ 3 Рис. 3. Семейство функций u(z) = z 8 δ 2z - z 3 + 8W∗1 для W1 = -(f (1) - 0.05) (линия 1), для W1 = -(f (z∗ )+f (1))/2 (линия 2), для W1 = -f (z∗ ) (линия 3), для W1 = -(f (z∗ )+0.05) (линия 4), когда δ = 1, ∗ = 1 ∗ 3 [Figure 3. The family of functions u(z) = z 8 δ 2z - z 3 + 8W∗1 for W1 = -(f (1) - 0.05) (line 1), for W1 = -(f (z∗ ) + f (1))/2 (line 2), for W1 = -f (z∗ ) (line 3), for W1 = -(f (z∗ ) + 0.05) (line 4), when δ = 1, ∗ = 1] Проводя аналогичный анализ для функции v1 (z), можно получить следующую систему условий, определяющую расположение этой кривой относительно оси z: 1) -8W2 < 1 - единственный ноль, одна смена знака; 2) 1 = -8W2 - два √ различных нуля, одна смена знака; 4 2 3) 1 < -8W2 < √ - два различных нуля, две смены знака; 3 3 √ 4 2 4) √ = -8W2 - один ноль, функция v1 (z) не меняет знак; 3√3 4 2 5) √ < -8W2 - нулей нет, функция v1 (z) не меняет знак. 3 3 Отметим, что если вернуться к изначальным параметрам задачи, то получим W1 ∗ = νξ1 , gβDηh2 W2 = νξ2 . gβCηh2 Таким образом, наличие противотечений в слое жидкости, их число и поνξ1 ложение застойных точек зависят от значения комбинаций параметров gβDηh 2 и νξ2 . gβCηh2 7. Годограф скорости. Особый интерес в этой ситуации представляет изучение траектории, описываемой концом вектора-годографа скорости. Чисто теоретически может быть не более 25 различных форм годографа (5 вариантов для кривой u(z) и 5 вариантов для кривой v(z)), но есть ситуации, которые порождают один и тот же годограф (например, при ∗ = 1 188 Крупномасштабная слоистая стационарная конвекция вязкой несжимаемой жидкости. . . 1 годографом будет биссекв случаях W1 = W2 = -1 и W1 = W2 = - 3√ 6 триса первой четверти). Кроме того, не следует забывать про связь между параметрами задачи (W1 ξ2 = W2 ξ1 ), которая приводит к дополнительным ограничениям и, следовательно, к снижению числа принципиально разных траекторий, описываемых концом годографа. Итак, если функция u1 (z) такова, что имеет единственный ноль и меняет знак один раз, то годограф имеет вид, представленный на рис. 4, a. Если функция u1 (z) имеет два нуля и меняет знак один раз, то годограф имеет вид, представленный на рис. 4, b. Если функция u1 (z) имеет два нуля, но меняет знак два раза, то годограф имеет вид, представленный на рис. 4, c. Если функция u1 (z) имеет один ноль, но не меняет знак, то годограф имеет вид, представленный на рис. 4, d. Если у функции u1 (z) нулей нет (функция не меняет знак), то годограф имеет вид, представленный на рис. 4, e. 8. Исследование завихренности и спиральности. Выпишем основную количественную меру завихренности - псевдовектор вихря или завихренность: i j k ∂ ∂ ∂ Ω = rot V = ∂x , ∂y ∂z u v 0 или покоординатно: Ωx = - δ3 ∂v = - (2W2 + z - z 3 ), ∂z 2 Ωy = ∂u δ3 = 2W1 + ∂z 2 ∗ (z - z 3 ) , Ωz = 0. Особый интерес здесь вызывают комбинации параметров, при которых завихренность обращается в ноль хотя бы в нескольких точках. Принципиально таких случая два. Случай первый. При ∗ = W1 = W2 = 0 завихренность на интересующем нас интервале обращается в ноль в точках z = 0 и z = 1. В этом случае u(z) ≡ 0, v(z) = δ3 2 z (2 - z 2 ). 8 То есть противотечений в слое не возникает (поскольку нули функции v(z) лежат вне интервала (0, 1)). Жидкость в слое движется в одном строго зафиксированном направлении и имеет место «волновой» эффект изменения модуля скорости движения жидкости по высоте z рассматриваемого слоя. Случай второй. При ∗ = 0, ∗ W2 = W1 . Так как ξ1 W2 = ξ2 W1 , последнее условие равносильно условию ξ2 = ξ1 . В этом случае псевдовектор вихря становится нулевым в точках, являющихся корнями уравнения 2W2 + z - z 3 = 0. Число действительных корней этого уравнения зависит от W2 . Экстре2 мум функции 2W2 + z - z 3 равен 3√ и достигается он в точке z∗∗ = √13 . 3 189 Рис. 4. Вид годографа в плоскости u(z) - v(z) [Figure 4. Hodograph plots in the u(z) - v(z) plane] Б у р м а ш е в а Н. В., П р о с в и р я к о в Е. Ю. 190 Крупномасштабная слоистая стационарная конвекция вязкой несжимаемой жидкости. . . 1 , то у приведенного выше Таким образом, если W2 > 0 или W2 < - 3√ 3 1 уравнения корней на отрезке [0, 1] нет. Если W2 = - 3√ , то существует 3 единственный действительный корень на отрезке [0, 1] - точка z∗∗ . Если 1 - 3√ W2 0, то таких корня два - один левее z∗∗ , другой правее. 3 Кроме того, u(z) = zδ 3 8 ∗ 2z - z 3 + 8W1 ∗ = ∗ zδ 8 3 (2z - z 3 + 8W2 ) = ∗ v(z), то есть имеет место линейная зависимость между компонентами u и v вектора скорости V . В этом случае годограф скорости описывает своим концом прямую, проходящую через начало координат, тангенс угла наклона которой равен 1/ ∗ . А значит, если координата z такова, что удовлетворяет уравнению 2W2 + z - z 3 = 0, то направление движения жидкости на этой высоте слоя оказывается заранее известным. Таким образом, в зависимости от сочетаний конкретных значений параметров задачи можно попасть в ситуации, когда вихрь меняет свое направление. Исследуем теперь спиральность [39]. По определению, спиральность есть число H = V · Ω = V · rot V = -u ∂u δ6 ∂v +v = (W1 - ∂z ∂z 8 ∗ W2 )z 2 (-2 + 3z 2 ). Очевидно, что спиральность оказывается равной нулю при W1 - ∗ W2 = 0, а это есть одно из условий обращения завихренности в ноль. Это выглядит логичным, учитывая, что завихренность Ω входит множителем в определение спиральности H. Однако есть два случая, когда спиральность обращается в ноль, а завихренность может и не быть нулевой - это случай z = 0 и z = 2/3 (нас интересуют только точки слоя 0 z 1). В точке z = 0 вектор скорости V оказывается равен нулю (условие прилипания), а компоненты завихренности определяются выражениями Ωx = -W2 δ 3 , Ωy = W1 δ 3 . Если W1 = W2 = 0, то завихренность наряду со спиральностью становится нулевой при z = 0. Аналогичной подстановкой точки z = 2/3 легко можно ∗ показать, что если W1 = - 3√6 , то обращаются в ноль u( 2/3) и Ωy , а если 1 W2 = - 3√ , то также в ноль обращаются v( 2/3) и Ωx . 6 Таким образом, даже если спиральность стала нулевой, то в зависимости от значений констант в граничных условиях завихренность может в ноль и не обращаться. Конкурирующие интересы. Мы заявляем, что у нас нет конфликта интересов в отношении авторства и публикации этой статьи. Авторская ответственность. Мы несем полную ответственность за предоставление окончательной рукописи в печать. Каждый из нас одобрил окончательную версию рукописи. 191

About the authors

Natalya V Burmasheva

Ural Federal University named after the First President of Russia B. N. Yeltsin; Institute of Engineering Science, Urals Branch, Russian Academy of Sciences

Email: nat_burm@mail.ru
Cand. Techn. Sci.; Associate Professor; Institute of Mathematics and Computer Science, Dept. of Mechanics and Mathematical Modeling1 ; Researcher; Laboratory of Mechanics of Deformations 2 19, Mira st., Ekaterinburg, 620002, Russian Federation; 34, Komsomolskaya st., Ekaterinburg, 620049, Russian Federation

Evgeny Yu Prosviryakov

Institute of Engineering Science, Urals Branch, Russian Academy of Sciences

Email: evgen_pros@mail.ru
Cand. Phys. & Math. Sci.; Leading Researcher; Laboratory of Applied Mechanics 34, Komsomolskaya st., Ekaterinburg, 620049, Russian Federation

References

Supplementary files