Прикладная математика и механика, 2023, T. 87, № 2, стр. 211-225

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

Е. И. Роменский 1*, И. М. Пешков 2**

1 Институт математики им. С.Л. Соболева СО РАН
Новосибирск, Россия

2 Университет Тренто
Тренто, Италия

* E-mail: evrom@math.nsc.ru
** E-mail: Ilya.Peshkov@unitn.it

Поступила в редакцию 02.02.2023
После доработки 05.03.2023
Принята к публикации 05.03.2023

Полный текст (PDF)

Аннотация

Представлена модель двухфазного течения сжимаемых несмешивающихся жидкостей, вывод которой основан на использовании теории симметрических гиперболических термодинамически согласованных систем. Модель является расширением предложенной ранее термодинамически согласованной модели сжимаемых двухфазных течений за счет включения новых переменных состояния среды, связанных с силами поверхностного натяжения. Определяющие уравнения модели образуют гиперболическую систему дифференциальных уравнений первого порядка и удовлетворяют законам термодинамики (сохранение энергии и возрастание энтропии). Исследованы свойства уравнений модели, показано, что закон капиллярного давления Юнга–Лапласа выполнен в асимптотическом приближении на континуальном уровне.

Ключевые слова: двухфазное течение, поверхностное натяжение, гиперболические уравнения

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

В настоящее время, начиная с работы [1] (см. [2, 3]), широко распространенным является подход, в котором поверхность раздела жидкостей моделируется градиентом функции цвета. Определяющие уравнения в таком подходе содержат параболические члены, что создает определенные трудности при их численной реализации. Недавно в работах [46] была предложена гиперболическая переформулировка модели, которая по существу основана на использовании в уравнениях переменной состояния в виде вектора, соответствующего градиенту функции цвета. Модифицированная модель дает те же самые результаты, что и [1]. Следует, однако, отметить, что неясна возможность применения этого подхода для описания дисперсных смесей жидкостей с поверхностным натяжением на континуальном уровне.

Целью работы является применение формализма симметрических гиперболических термодинамически согласованных (СГТС) систем для создания модели двух несмешивающихся жидкостей с учетом поверхностного натяжения. Класс СГТС-уравнений включает в себя многие известные уравнения механики сплошной среды [7, 8] и, кроме этого, на основе принципов формулировки уравнений этого класса, можно строить гиперболические модели сред с усложненными свойствами, в частности, многофазных сред [9, 10]. Формулировка определяющих уравнений основана на расширении СГТС-модели многофазных сжимаемых течений жидкостей, предложенной в работе [11]. Расширение приводит к нелокальной модели и производится за счет введения новой переменной состояния – градиента объемной доли одной из фаз. При этом соответствующая объемная доля играет роль функции цвета как это принято, например, в [2]. В соответствии с формализмом СГТС-моделей для введенной новой переменной должна быть определена парная переменная состояния (которая характеризует инерциальное изменение объемной доли) и дополнительный закон сохранения для нее. Для описания диссипативных и дисперсионных эффектов вводятся правые части уравнений релаксационного типа, при этом важную роль для моделирования поверхностного натяжения играет антисимметричная их часть, не приводящая к возрастанию энтропии.

Полученная система определяющих дифференциальных уравнений является гиперболической и ее решения удовлетворяют законам термодинамики (сохранение энергии и возрастание энтропии). Проанализированы характеристические скорости в одномерном случае, исследованы релаксационные пределы для малых параметров релаксации и показано, что модель на континуальном уровне описывает закон поверхностного натяжения Юнга–Лапласа.

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

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

2.1. Порождающая система для обратимых сжимаемых двухфазных течений. Начнем с формулировки общих термодинамически согласованных уравнений без учета диссипации. Система порождающих дифференциальных уравнений для сжимаемых двухфазных течений сформулирована в [11, 12] и имеет следующий вид

$\begin{gathered} \frac{{\partial \rho \alpha }}{{\partial t}} + \frac{{\partial \rho \alpha {{u}_{k}}}}{{\partial {{x}_{k}}}} = 0 \\ \frac{{\partial \rho }}{{\partial t}} + \frac{{\partial \rho {{u}_{k}}}}{{\partial {{x}_{k}}}} = 0 \\ \end{gathered} $
(2.1)
$\begin{gathered} \frac{{\partial \rho {{u}_{l}}}}{{\partial t}} + \frac{\partial }{{\partial {{x}_{k}}}}\left( {\rho {{u}_{l}}{{u}_{k}} + {{\rho }^{2}}{{E}_{\rho }}{{\delta }_{{lk}}} + \rho {{w}_{l}}{{E}_{{{{w}_{k}}}}}} \right) = 0 \\ \frac{{\partial \rho c}}{{\partial t}} + \frac{\partial }{{\partial {{x}_{k}}}}\left( {\rho {{u}_{k}}c + \rho {{E}_{{{{w}_{k}}}}}} \right) = 0 \\ \end{gathered} $
$\begin{gathered} \frac{{\partial {{w}_{k}}}}{{\partial t}} + \frac{\partial }{{\partial {{x}_{k}}}}\left( {{{u}_{l}}{{w}_{l}} + {{E}_{c}}} \right) + {{u}_{l}}\left( {\frac{{\partial {{w}_{k}}}}{{\partial {{x}_{l}}}} - \frac{{\partial {{w}_{l}}}}{{\partial {{x}_{k}}}}} \right) = 0 \\ \frac{{\partial \rho S}}{{\partial t}} + \frac{{\partial \rho S{{u}_{k}}}}{{\partial {{x}_{k}}}} = 0 \\ \end{gathered} $

Для описания течения использованы следующие переменные состояния: $\alpha = {{\alpha }_{1}}$ объемная концентрация (объемная доля) фазы с номером 1, удовлетворяющая условию насыщения, что означает, что объемная доля фазы с номером 2 вычисляется как ${{\alpha }_{2}} = 1 - {{\alpha }_{1}}$, $\rho = {{\alpha }_{1}}{{\rho }_{1}} + {{\alpha }_{2}}{{\rho }_{2}}$ – массовая плотность смеси, ${{\rho }_{1}}$, ${{\rho }_{2}}$ – массовые плотности фаз, ${{u}_{i}} = {{c}_{1}}u_{i}^{1} + {{c}_{2}}u_{i}^{2}$ – скорость смеси, где $u_{i}^{1}$, $u_{i}^{2}$ – скорости фаз, $c = {{c}_{1}} = {{\alpha }_{1}}{{\rho }_{1}}{\text{/}}\rho $, ${{c}_{2}} = 1 - c$ = ${{\alpha }_{2}}{{\rho }_{2}}{\text{/}}\rho $ – массовые доли фаз $\left( {{{c}_{1}} + {{c}_{2}} = 1} \right)$. Далее, ${{w}_{i}} = u_{i}^{1} - u_{i}^{2}$ – относительная скорость фаз, $S$ – энтропия смеси (подчеркнем, что при формулировке модели используется приближение одной энтропии, что допустимо при небольших вариациях температур фаз). E – обобщенная внутренняя энергия смеси, а индексы внизу обозначают производные по соответствующим переменным состояния.

Слагаемое ${{u}_{l}}\left( {\frac{{\partial {{w}_{k}}}}{{\partial {{x}_{l}}}} - \frac{{\partial {{w}_{l}}}}{{\partial {{x}_{k}}}}} \right)$ в уравнении для относительной скорости, содержащее ротор относительной скорости $\nabla \times {\mathbf{w}}$, играет очень важную роль в структуре определяющих уравнений (2.1) (см. [11]). Оно обеспечивает существование дополнительного закона сохранения энергии для системы (2.1), а также приведение этой системы к симметрической форме. Можно доказать, что если дополнительный стационарный закон сохранения $\nabla \times {\mathbf{w}} = 0$ выполнен в начальных данных при $t = 0$, то он выполнен для всех последующих моментов времени $t > 0$.

При помощи суммирования всех уравнений, а также стационарного закона сохранения $\nabla \times {\mathbf{w}} = 0$, умноженных на соответствующие множители, можно показать, что решения (2.1) удовлетворяют дополнительному закону сохранения энергии (см. [11])

(2.2)
$\frac{\partial }{{\partial t}}\rho \left( {E + \frac{1}{2}{{u}_{l}}{{u}_{l}}} \right) + \frac{\partial }{{\partial {{x}_{k}}}}\left( {\rho {{u}_{k}}\left( {E + \frac{1}{2}{{u}_{l}}{{u}_{l}} + \frac{p}{\rho }} \right) + \rho {{u}_{l}}{{w}_{l}}{{E}_{{{{w}_{k}}}}} + \rho {{E}_{c}}{{E}_{{{{w}_{k}}}}}} \right) = 0$

Здесь, $E\left( {\alpha ,\rho ,c,{{w}_{k}}} \right)$ – плотность внутренней энергии смеси. Наличие дополнительного закона сохранения позволяет переформулировать (2.1) в терминах так называемого порождающего термодинамического потенциала и порождающих переменных и привести уравнения к симметрической форме [7, 8, 13, 14].

Для (2.1) необходимо иметь только одно замыкающее соотношение, а именно плотность внутренней энергии $E$, которая должна быть указана для каждого конкретного течения смеси. В двухфазных моделях, сформулированных в [11], внутренняя энергия определяется как сумма усредненных по массе внутренних энергий фаз и кинематической энергии относительного движения

(2.3)
$E\left( {\alpha ,\rho ,c,{{w}_{k}}} \right) = c{{e}_{1}}\left( {{{\rho }_{1}},S} \right) + \left( {1 - c} \right){{e}_{2}}\left( {{{\rho }_{2}},S} \right) + c\left( {1 - c} \right)\frac{{{{w}_{l}}{{w}_{l}}}}{2}$

Заметим, что внутренние энергии фаз можно записать в виде зависимостей ${{e}_{i}}\left( {\rho {{c}_{i}}{\text{/}}{{\alpha }_{i}},S} \right)$, $i = 1,2$, которые могут оказаться полезными при изучении выпуклости полной энергии.

Для симметрической гиперболичности по Фридрихсу [1416] необходимо потребовать выпуклость энергии (2.3) по переменным $\alpha ,\rho ,c$, ${{w}_{k}}$. Хотя строгого доказательства выпуклости не имеется, проведенные многочисленные расчеты по этой модели не выявили аномалий в поведении решений, что является подтверждением выпуклости, доказательство которой является предметом дальнейших исследований.

Если внутренняя энергия смеси определена, то все термодинамические силы (производные энергии по переменным состояния) можно вычислить по формулам (см. [11]):

(2.4)
$\begin{gathered} p = {{\rho }^{2}}\frac{{\partial E}}{{\partial \rho }} = {{\alpha }_{1}}{{p}_{1}} + {{\alpha }_{2}}{{p}_{2}},\quad \frac{{\partial E}}{{\partial {{w}_{k}}}} = {{c}_{1}}{{c}_{2}}{{w}_{k}} \\ \frac{{\partial E}}{{\partial c}} = {{e}_{1}} + \frac{{{{p}_{1}}}}{{{{\rho }_{1}}}} - {{e}_{2}} - \frac{{{{p}_{2}}}}{{{{\rho }_{2}}}} + \left( {1 - 2c} \right)\frac{{{{w}_{l}}{{w}_{l}}}}{2},\quad {{E}_{\alpha }} = \frac{{{{p}_{2}} - {{p}_{1}}}}{\rho }, \\ \end{gathered} $
где давления компонент смеси вычисляются как ${{p}_{i}} = \rho _{i}^{2}\frac{{\partial {{e}_{i}}}}{{\partial {{\rho }_{i}}}}$, $i = 1,2$.

2.2. Порождающая система для необратимых сжимаемых двухфазных течений. Типичные необратимые процессы в двухфазных сжимаемых течениях – это релаксация фазового давления к общему значению давления фаз и межфазное трение. Эти процессы учитываются в (2.1) как источниковые алгебраические члены в уравнениях для объемной доли и для относительной скорости. В итоге, определяющие уравнения модели с необратимыми процессами имеют следующий вид

$\frac{{\partial \rho \alpha }}{{\partial t}} + \frac{{\partial \rho \alpha {{u}_{k}}}}{{\partial {{x}_{k}}}} = - {{\lambda }^{{ - 1}}}\rho {{E}_{\alpha }}$
$\frac{{\partial \rho }}{{\partial t}} + \frac{{\partial \rho {{u}_{k}}}}{{\partial {{x}_{k}}}} = 0$
(2.5)
$\begin{gathered} \frac{{\partial \rho {{u}_{l}}}}{{\partial t}} + \frac{\partial }{{\partial {{x}_{k}}}}\left( {\rho {{u}_{l}}{{u}_{k}} + {{\rho }^{2}}{{E}_{\rho }}{{\delta }_{{lk}}} + \rho {{w}_{l}}{{E}_{{{{w}_{k}}}}}} \right) = 0 \\ \frac{{\partial \rho c}}{{\partial t}} + \frac{\partial }{{\partial {{x}_{k}}}}\left( {\rho {{u}_{k}}c + \rho {{E}_{{{{w}_{k}}}}}} \right) = 0 \\ \end{gathered} $
$\frac{{\partial {{w}_{k}}}}{{\partial t}} + \frac{\partial }{{\partial {{x}_{k}}}}\left( {{{u}_{l}}{{w}_{l}} + {{E}_{c}}} \right) + {{u}_{l}}\left( {\frac{{\partial {{w}_{k}}}}{{\partial {{x}_{l}}}} - \frac{{\partial {{w}_{l}}}}{{\partial {{x}_{k}}}}} \right) = - {{\chi }^{{ - 1}}}{{E}_{{{{w}_{k}}}}}$
$\frac{{\partial \rho S}}{{\partial t}} + \frac{{\partial \rho S{{u}_{k}}}}{{\partial {{x}_{k}}}} = Q = \frac{\rho }{{{{E}_{S}}}}\left( {{{\lambda }^{{ - 1}}}{{E}_{\alpha }}{{E}_{\alpha }} + {{\chi }^{{ - 1}}}{{E}_{{{{w}_{k}}}}}{{E}_{{{{w}_{k}}}}}} \right) \geqslant 0$

Можно увидеть, что диссипативные члены (источники в правой части первого и пятого уравнения в системе (2.5)) выбираются пропорциональными термодинамическим силам, соответствующим переменным уравнений, в которые они введены. Предполагая, что в состоянии термодинамического равновесия термодинамические силы $\rho {{E}_{\alpha }}$ и ${{E}_{{{{w}_{k}}}}}$ равны нулю, можно видеть, что в термодинамическом равновесии диссипация отсутствует. Учет необратимости приводит к появлению источника $Q$ производства энтропии в уравнении ее баланса. Заметим, что производство энтропии неотрицательно ($Q \geqslant 0$) в силу выбора диссипативных источников пропорциональными термодинамическим силам. Введение диссипации не влияет на закон сохранения энергии (2.2), и его вид остается тем же самым для системы (2.5)

(2.6)
$\frac{\partial }{{\partial t}}\rho \left( {E + \frac{1}{2}{{u}_{l}}{{u}_{l}}} \right) + \frac{\partial }{{\partial {{x}_{k}}}}\left( {\rho {{u}_{k}}\left( {E + \frac{1}{2}{{u}_{l}}{{u}_{l}} + \frac{p}{\rho }} \right) + \rho {{u}_{l}}{{w}_{l}}{{E}_{{{{w}_{k}}}}} + \rho {{E}_{c}}{{E}_{{{{w}_{k}}}}}} \right) = 0$

3. Порождающие определяющие уравнения для двухфазных течений с поверхностным натяжением.

3.1. Формулировка определяющих уравнений. В этом разделе представлена гиперболическая термодинамически согласованная формулировка модели двухфазного течения с поверхностным натяжением, основанная на расширении модели двухфазного течения, описанной в предыдущих разделах. Одним из ключевых положений в представленной выше двухфазной модели является закон баланса объемных долей, включенный в систему (2.5), который может быть записан в форме

(3.1)
$\frac{{\partial \alpha }}{{\partial t}} + {{u}_{k}}\frac{{\partial \alpha }}{{\partial {{x}_{k}}}} = - \frac{1}{\lambda }{{E}_{\alpha }}$
и описывает перенос значения объемной концентрации вдоль траектории движения со скоростью смеси ${{u}_{k}}$, а также ее изменение за счет релаксации давления обеспечиваемой источником в правой части, см. выражение для ${{E}_{\alpha }}$ в (2.4).

Для моделирования поверхностного натяжения в течении несмешивающихся жидкостей необходимо иметь информацию о поверхностях раздела. Общепринятый подход состоит в использовании некоего скалярного параметра $\phi $ (например, функции цвета – color function) для описания поверхностей раздела, который дает возможность определить континуальные характеристики сил поверхностного натяжения (см., напр., [1, 18], и приведенные там ссылки). Параметр $\phi $ принимает значение 0, если точка в среде принадлежит одной из фаз, и значение 1, если точка принадлежит другой фазе. Распределение поверхностей (плотность площади поверхности границ раздела в единице объема и их усредненная ориентация) на континуальном уровне может характеризоваться градиентом $\nabla \phi $ поля $\phi $, а силы поверхностного натяжения характеризуются тензором, образованным этим вектором. При этом, уравнения соответствующей модели содержат параболические члены со вторыми производными от параметра $\phi $.

Цель статьи – сформулировать гиперболические уравнения, содержащие только пространственные производные первого порядка, и принадлежащие классу симметрических гиперболических термодинамически согласованных систем. При этом, в представленной формулировке, в качестве параметра $\phi $ оказывается естественным выбрать объемную долю $\alpha $, которая как раз принимает значения 0 и 1 в зависимости от принадлежности точки среды к той или иной фазе. Для определения распределения поверхностей в элементе среды, также будем основываться на нелокальном описании сплошной среды с использованием некоторого вектора ${\mathbf{b}}$, пропорционального вектору градиента объемной доли $\nabla \alpha $ и стремящегося к нему в некотором релаксационном пределе.

Рассмотрим уравнение баланса объемной доли в виде

(3.2)
$\frac{{\partial \alpha }}{{\partial t}} + {{u}_{k}}\frac{{\partial \alpha }}{{\partial {{x}_{k}}}} + \mathcal{F} = 0,$
которое, аналогично (3.1), описывает перенос объемной доли вдоль траектории движения смеси и содержит источник $\mathcal{F}$, который мы определим ниже.

Введем новую переменную ${{b}_{n}} = \partial \alpha {\text{/}}\partial {{x}_{n}}$, уравнение для которой может быть получено дифференцированием (3.2) по ${{x}_{n}}$:

(3.3)
$\frac{{\partial {{b}_{n}}}}{{\partial t}} + {{u}_{k}}\frac{{\partial {{b}_{n}}}}{{\partial {{x}_{k}}}} + {{b}_{k}}\frac{{\partial {{u}_{k}}}}{{\partial {{x}_{n}}}} + \frac{{\partial{ \mathcal{F}}}}{{\partial {{x}_{n}}}} = 0$

Заметим, что это уравнение может быть записано в эквивалентном виде

(3.4)
$\frac{{\partial {{b}_{n}}}}{{\partial t}} + \frac{{\partial \left( {{{b}_{k}}{{u}_{k}} + \mathcal{F}} \right)}}{{\partial {{x}_{n}}}} + {{u}_{k}}\left( {\frac{{\partial {{b}_{n}}}}{{\partial {{x}_{k}}}} - \frac{{\partial {{b}_{k}}}}{{\partial {{x}_{n}}}}} \right) = 0$
и, кроме того, в силу определения ${{b}_{n}} = \partial \alpha {\text{/}}\partial {{x}_{n}}$ имеет место стационарное уравнение

(3.5)
$\frac{{\partial {{b}_{n}}}}{{\partial {{x}_{k}}}} - \frac{{\partial {{b}_{k}}}}{{\partial {{x}_{n}}}} = 0$

Это стационарное уравнение может быть также получено как следствие (3.4), если предположить, что (3.5) выполнено для начальных данных.

Уравнение (3.4) будет присутствовать в качестве одного из определяющих уравнений для модели поверхностного натяжения. В соответствии с формализмом термодинамически согласованных систем полная система уравнений имеет парную структуру, или гамильтонову структуру [7, 8] – для каждой переменной состояния (обобщенная координата) имеется парная переменная (обобщенный импульс), а их потоки при выводе закона сохранения энергии образуют произведение в потоке сохранения энергии. При этом поток одной переменной выражается через производную обобщенной энергии по другой, парной, переменной состояния. Таким образом, следует принять предположение о том, что поток $\mathcal{F}$ в уравнении (3.4), зависящий от переменных состояния среды, имеет вид $\mathcal{F} = \partial E{\text{/}}\partial d$. Здесь $d$ – скалярная переменная состояния, которая, как видно из уравнения (3.2), характеризует скорость изменения объемной доли (учет инерции) при изменении термодинамических параметров течения, аналогично тому как в уравнении (3.1) правая часть характеризует релаксацию давлений фаз к единому значению. При таком выборе $\mathcal{F}$, уравнение для объемной доли, записанное в дивергентной форме, принимает вид

(3.6)
$\frac{{\partial \rho \alpha }}{{\partial t}} + \frac{{\partial \rho \alpha {{u}_{k}}}}{{\partial {{x}_{k}}}} = - \rho {{E}_{d}}$

Для введенной новой переменной состояния $d$ естественно выбрать следующее парное к (3.4) определяющее уравнение:

(3.7)
$\frac{{\partial \rho d}}{{\partial t}} + \frac{\partial }{{\partial {{x}_{k}}}}\left( {\rho d{{u}_{k}} + \rho {{E}_{{{{b}_{k}}}}}} \right) = \rho {{E}_{\alpha }} - {{\varepsilon }^{{ - 1}}}\rho {{E}_{d}},$
где $\varepsilon $ предполагается константой. Как будет показано ниже, именно такая форма уравнения для $d$ обеспечивает закон сохранения энергии для полной системы. Заметим, что источники в правых частях уравнений (3.6), (3.7) содержат антисимметричные члены, которые моделируют изменение объемной доли при изменении термодинамических переменных состояния смеси.

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

(3.8)
${{B}_{k}} = \delta {{b}_{k}},\quad D = {{\delta }^{{ - 1}}}d,$
где $\delta $ – безразмерный параметр, который может быть интерпретирован как отношение масштабов длины

(3.9)
$\delta \sim L{\text{/}}\ell $

Здесь, длина $L$ ассоциирована с макромасштабом (например, с разрешающим масштабом инструментов измерения или разрешающим масштабом вычислительной сетки), а $\ell $ – некий характерный микромасштаб ассоциированный с границами раздела фаз.

С использованием этих новых переменных состояния уравнения запишутся следующим образом:

$\frac{{\partial \rho \alpha }}{{\partial t}} + \frac{{\partial \rho \alpha {{u}_{k}}}}{{\partial {{x}_{k}}}} = - \frac{1}{\delta }\rho {{E}_{D}}$
(3.10)
$\frac{{\partial \rho D}}{{\partial t}} + \frac{\partial }{{\partial {{x}_{k}}}}\left( {\rho D{{u}_{k}} + \rho {{E}_{{{{B}_{k}}}}}} \right) = \frac{1}{\delta }\rho {{E}_{\alpha }} - \frac{1}{\varepsilon }\rho {{E}_{D}}$
$\frac{{\partial {{B}_{k}}}}{{\partial t}} + \frac{\partial }{{\partial {{x}_{k}}}}\left( {{{B}_{l}}{{u}_{l}} + {{E}_{D}}} \right) + {{u}_{l}}\left( {\frac{{\partial {{B}_{k}}}}{{\partial {{x}_{l}}}} - \frac{{\partial {{B}_{l}}}}{{\partial {{x}_{k}}}}} \right) = 0$

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

Таким образом, полная система термодинамически согласованных определяющих уравнений для описания течения несмешивающихся жидкостей с поверхностным натяжением может быть получена из системы (2.5), дополненной уравнениями (3.10). При этом в уравнении для импульса в потоке появится дополнительное слагаемое, и, кроме того, изменится формула для производства энтропии $Q$. Полная система будет выглядеть следующим образом:

$\frac{{\partial \rho }}{{\partial t}} + \frac{{\partial \rho {{u}_{k}}}}{{\partial {{x}_{k}}}} = 0$
$\frac{{\partial \rho {{u}_{l}}}}{{\partial t}} + \frac{\partial }{{\partial {{x}_{k}}}}\left( {\rho {{u}_{l}}{{u}_{k}} + {{\rho }^{2}}{{E}_{\rho }}{{\delta }_{{lk}}} + \rho {{w}_{l}}{{E}_{{{{w}_{k}}}}} + \rho {{B}_{l}}{{E}_{{{{B}_{k}}}}}} \right) = 0$
$\frac{{\partial \rho c}}{{\partial t}} + \frac{\partial }{{\partial {{x}_{k}}}}\left( {\rho {{u}_{k}}c + \rho {{E}_{{{{w}_{k}}}}}} \right) = 0$
(3.11)
$\frac{{\partial {{w}_{k}}}}{{\partial t}} + \frac{\partial }{{\partial {{x}_{k}}}}\left( {{{u}_{l}}{{w}_{l}} + {{E}_{c}}} \right) + {{u}_{l}}\left( {\frac{{\partial {{w}_{k}}}}{{\partial {{x}_{l}}}} - \frac{{\partial {{w}_{l}}}}{{\partial {{x}_{k}}}}} \right) = - \frac{1}{\chi }{{E}_{{{{w}_{k}}}}}$
$\frac{{\partial \rho \alpha }}{{\partial t}} + \frac{{\partial \rho \alpha {{u}_{k}}}}{{\partial {{x}_{k}}}} = - \frac{1}{\delta }\rho {{E}_{D}}$
$\frac{{\partial \rho D}}{{\partial t}} + \frac{\partial }{{\partial {{x}_{k}}}}\left( {\rho {{u}_{k}}D + \rho {{E}_{{{{B}_{k}}}}}} \right) = \frac{1}{\delta }\rho {{E}_{\alpha }} - \frac{1}{\varepsilon }\rho {{E}_{D}}$
$\frac{{\partial {{B}_{k}}}}{{\partial t}} + \frac{\partial }{{\partial {{x}_{k}}}}\left( {{{B}_{l}}{{u}_{l}} + {{E}_{D}}} \right) + {{u}_{l}}\left( {\frac{{\partial {{B}_{k}}}}{{\partial {{x}_{l}}}} - \frac{{\partial {{B}_{l}}}}{{\partial {{x}_{k}}}}} \right) = 0$
$\frac{{\partial \rho S}}{{\partial t}} + \frac{{\partial \rho S{{u}_{k}}}}{{\partial {{x}_{k}}}} = Q = \frac{\rho }{{{{E}_{S}}}}\left( {{{\chi }^{{ - 1}}}{{E}_{{{{w}_{k}}}}}{{E}_{{{{w}_{k}}}}} + {{\varepsilon }^{{ - 1}}}{{E}_{D}}{{E}_{D}}} \right)$

Для замыкания системы необходимо задать зависимость плотности внутренней энергии от переменных состояния среды $E\left( {\alpha ,\rho ,c,D,{{w}_{k}},{{B}_{k}},S} \right)$. Ниже мы приведем конкретный вид подходящей формулы для плотности внутренней энергии, обобщающий (2.3).

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

(3.12)
$\begin{gathered} {{q}_{1}} = E + \rho {{E}_{\rho }} - \frac{1}{2}{{u}_{l}}{{u}_{l}} - c{{E}_{c}} - \alpha {{E}_{\alpha }} - D{{E}_{D}} - S{{E}_{S}},\quad {{u}_{l}},\quad n = {{E}_{c}} \\ {{z}_{k}} = \rho {{E}_{{{{w}_{k}}}}},\quad {{q}_{2}} = {{E}_{\alpha }},\quad r = {{E}_{D}},\quad {{j}_{k}} = \rho {{E}_{{{{B}_{k}}}}},\quad \theta = {{E}_{S}} \\ \end{gathered} $

В итоге, получим уравнение

(3.13)
$\begin{gathered} \frac{\partial }{{\partial t}}\rho \left( {E + {{u}_{l}}{{u}_{l}}{\text{/}}2} \right) + \\ + \;\frac{\partial }{{\partial {{x}_{k}}}}\left( {\rho {{u}_{k}}\left( {E + \tfrac{1}{2}{{u}_{l}}{{u}_{l}} + \frac{p}{\rho }} \right) + \rho {{u}_{l}}\left( {{{w}_{l}}{{E}_{{{{w}_{k}}}}} + {{B}_{l}}{{E}_{{{{B}_{k}}}}}} \right) + \rho {{E}_{c}}{{E}_{{{{w}_{k}}}}} + \rho {{E}_{D}}{{E}_{{{{B}_{k}}}}}} \right) = 0 \\ \end{gathered} $

Отметим, что производство энтропии в вышеприведенной системе (3.11) неотрицательно в силу выбора источников в правых частях уравнений. Как оказалось, антисимметричные части описывающих источники членов в правых частях уравнений для $\alpha $ и $D$ при выводе уравнения энергии взаимно сократились. Тем не менее, они должны оказывать существенное влияние на волновые поля, описываемые моделью и, по-видимому, описывают дисперсионные эффекты.

Энергия $E$ может быть взята в форме

(3.14)
$E = {{c}_{1}}{{e}_{1}}\left( {{{\rho }_{1}},S} \right) + {{c}_{2}}{{e}_{2}}\left( {{{\rho }_{2}},S} \right) + {{c}_{1}}{{c}_{2}}\frac{{{{w}_{l}}{{w}_{l}}}}{2} + \frac{1}{2}\frac{\Sigma }{\rho }{{B}_{l}}{{B}_{l}} + \frac{1}{2}\frac{{{{G}^{2}}}}{\rho }{{D}^{2}},$
где $G$ – модуль, характеризующий инерционные эффекты, связанные с физикой поверхностного натяжения, $\Sigma $ называют коэффициентом капиллярности , и его связь с традиционным коэффициентом поверхностного натяжения $\sigma $ обсуждается в параграфе 3.1. Отметим только, что физические размерности $\Sigma $ и $\sigma $ соотносятся как $\left[ \Sigma \right] = \left[ \sigma \right] \cdot L$, где $L$ – некий масштаб длины.

Для энергии в форме (3.14) термодинамические силы вычисляются по формулам

$p = {{\rho }^{2}}\frac{{\partial E}}{{\partial \rho }} = {{\alpha }_{1}}{{p}_{1}} + {{\alpha }_{2}}{{p}_{2}} - \frac{1}{2}\frac{\Sigma }{{{{\rho }^{2}}}}{{B}_{l}}{{B}_{l}} - \frac{1}{2}\frac{{{{G}^{2}}}}{{{{\rho }^{2}}}}{{D}^{2}}\quad \left( {{{p}_{i}} = \rho _{i}^{2}\frac{{\partial {{e}_{i}}}}{{\partial {{\rho }_{i}}}},\;i = 1,2} \right)$
(3.15)
$\frac{{\partial E}}{{\partial c}} = {{e}_{1}} + \frac{{{{p}_{1}}}}{{{{\rho }_{1}}}} - {{e}_{2}} - \frac{{{{p}_{2}}}}{{{{\rho }_{2}}}} + \left( {1 - 2c} \right)\frac{{{{w}_{l}}{{w}_{l}}}}{2},\quad {{E}_{\alpha }} = \frac{{{{p}_{2}} - {{p}_{1}}}}{\rho }$
$\frac{{\partial E}}{{\partial {{w}_{k}}}} = {{c}_{1}}{{c}_{2}}{{w}_{k}},\quad \frac{{\partial E}}{{\partial D}} = \frac{{{{G}^{2}}}}{\rho }D,\quad \frac{{\partial E}}{{\partial {{B}_{k}}}} = \frac{\Sigma }{\rho }{{B}_{k}}$

3.2. Асимптотический анализ и закон Юнга–Лапласа. В этом разделе мы покажем, что при формальном асимптотическом переходе при $\delta \to 0$, на решениях предложенной модели, в первом приближении, выполняется закон Юнга–Лапласа, связывающий кривизну поверхности раздела фаз и скачок давления при переходе через эту границу.

Поскольку физические единицы параметров $\varepsilon $ и ${{\lambda }^{{ - 1}}}$ совпадают, то в дальнейшем удобно предполагать, что

(3.16)
$\varepsilon = {{\lambda }^{{ - 1}}}{{\delta }^{2}}$

Тогда пятое–седьмое уравнения системы (3.11) могут быть записаны как

$\rho \frac{{d\alpha }}{{dt}} = - \frac{1}{\delta }\rho {{E}_{D}}$
(3.17)
$\frac{{d{{B}_{k}}}}{{dt}} + {{B}_{j}}\frac{{\partial {{u}_{j}}}}{{\partial {{x}_{k}}}} + \frac{{\partial {{E}_{D}}}}{{\partial {{x}_{k}}}} = 0$
$\rho \frac{{dD}}{{dt}} + \frac{{\partial \rho {{E}_{{{{B}_{k}}}}}}}{{\partial {{x}_{k}}}} = - \frac{\lambda }{\delta }\left( {\frac{1}{\delta }\rho {{E}_{D}} - \frac{1}{\lambda }\rho {{E}_{\alpha }}} \right),$
где для краткости мы используем материальную производную $d{\text{/}}dt$ = $\partial {\text{/}}\partial t$ + ${{u}_{l}}\partial {\text{/}}\partial {{x}_{l}}$.

Определим поведение решения при $\delta \to 0$ при постоянном $\lambda $. Разложим решение в формальный ряд по степеням малого параметра $\delta $:

$\alpha = {{\alpha }_{0}} + \delta {{\alpha }_{1}} + {{\delta }^{2}}{{\alpha }_{2}} + \ldots $
(3.18)
$\begin{gathered} \rho = {{\rho }_{0}} + \delta {{\rho }_{1}} + {{\delta }^{2}}{{\rho }_{2}} + \ldots \\ D = {{D}_{0}} + \delta {{D}_{1}} + {{\delta }^{2}}{{D}_{2}} + \ldots \\ \end{gathered} $
${{B}_{k}} = {{B}_{{k,0}}} + \delta {{B}_{{k,1}}} + {{\delta }^{2}}{{B}_{{k,2}}} + \ldots $

То же самое проделаем с термодинамическими силами:

${{E}_{\alpha }} = {{E}_{{\alpha ,0}}} + \delta {{E}_{{\alpha ,1}}} + {{\delta }^{2}}{{E}_{{\alpha ,1}}} + \ldots $
(3.19)
${{E}_{D}} = {{E}_{{D,0}}} + \delta {{E}_{{D,1}}} + {{\delta }^{2}}{{E}_{{D,1}}} + \ldots $
${{E}_{{{{B}_{k}}}}} = {{E}_{{{{B}_{k}},0}}} + \delta {{E}_{{{{B}_{k}},1}}} + {{\delta }^{2}}{{E}_{{{{B}_{k}},2}}} + \ldots ,$
для которых мы должны предположить
(3.20)
${{E}_{{\alpha ,0}}} = 0,$
что означает равенство давлений в элементе среды в равновесии, см. (3.15).

Уравнение для $D$:

Подставляя (3.18) и (3.19) в уравнение для $D$, и, собирая коэффициенты при одинаковых степенях $\delta $, получим, что

(3.21)
${{E}_{{D,0}}} = 0,\quad {{E}_{{D,1}}} = 0$

Далее, сравнивая коэффициенты при ${{\delta }^{2}}$, получим

(3.22)
${{E}_{{D,2}}} = \frac{1}{{\lambda {{\rho }_{0}}}}\left[ {{{\rho }_{0}}{{E}_{{\alpha ,1}}} - \frac{\partial }{{\partial {{x}_{k}}}}\left( {{{\rho }_{0}}{{E}_{{{{B}_{k}},0}}}} \right)} \right]$

Уравнение для ${{B}_{k}}$:

Рассмотрим далее разложение в уравнении для ${{B}_{k}}$. При нулевой и первой степени по $\delta $ получим

(3.23)
$\frac{{d{{B}_{{k,0}}}}}{{dt}} + {{B}_{{j,0}}}\frac{{\partial {{u}_{j}}}}{{\partial {{x}_{k}}}} = 0,\quad \frac{{d{{B}_{{k,1}}}}}{{dt}} + {{B}_{{j,1}}}\frac{{\partial {{u}_{j}}}}{{\partial {{x}_{k}}}} = 0,$
а при второй степени:

(3.24)
$\frac{{d{{B}_{{k,2}}}}}{{dt}} + {{B}_{{j,2}}}\frac{{\partial {{u}_{j}}}}{{\partial {{x}_{k}}}} + \frac{{\partial {{E}_{{D,2}}}}}{{\partial {{x}_{k}}}} = 0$

С учетом (3.22) последнее уравнение запишется как

(3.25)
$\frac{{d{{B}_{{k,2}}}}}{{dt}} + {{B}_{{j,2}}}\frac{{\partial {{u}_{j}}}}{{\partial {{x}_{k}}}} + \frac{1}{\lambda }\frac{{\partial {{E}_{{\alpha ,1}}}}}{{\partial {{x}_{k}}}} - \frac{1}{\lambda }\frac{\partial }{{\partial {{x}_{k}}}}\left( {\rho _{0}^{{ - 1}}\frac{{\partial {{\rho }_{0}}{{E}_{{{{B}_{j}},0}}}}}{{\partial {{x}_{j}}}}} \right) = 0$

Уравнение для $\alpha $:

Наконец, рассмотрим уравнение для $\alpha $. При первой степени по $\delta $, получаем

(3.26)
$\frac{{d{{\alpha }_{0}}}}{{dt}} = 0,$
a при второй степени:

(3.27)
${{\rho }_{0}}\frac{{d{{\alpha }_{1}}}}{{dt}} - \frac{1}{\lambda }\frac{{\partial {{\rho }_{0}}{{E}_{{{{B}_{k}},0}}}}}{{\partial {{x}_{k}}}} = - \frac{1}{\lambda }{{\rho }_{0}}{{E}_{{\alpha ,1}}}$

Таким образом, при старших степенях $\delta $, определяющие уравнения для $\alpha $ и ${{B}_{k}}$ запишутся как

(3.28)
$\begin{gathered} {{\rho }_{0}}\frac{{d{{\alpha }_{1}}}}{{dt}} - \frac{1}{\lambda }\frac{{\partial {{\rho }_{0}}{{E}_{{{{b}_{{k,0}}}}}}}}{{\partial {{x}_{k}}}} = - \frac{1}{\lambda }{{\rho }_{0}}{{E}_{{\alpha ,1}}} \\ \frac{{d{{B}_{{k,2}}}}}{{dt}} + {{B}_{{j,2}}}\frac{{\partial {{u}_{j}}}}{{\partial {{x}_{k}}}} + \frac{1}{\lambda }\frac{{\partial {{E}_{{\alpha ,1}}}}}{{\partial {{x}_{k}}}} - \frac{1}{\lambda }\frac{\partial }{{\partial {{x}_{k}}}}\left( {\rho _{0}^{{ - 1}}\frac{{\partial {{\rho }_{0}}{{E}_{{{{B}_{j}},0}}}}}{{\partial {{x}_{j}}}}} \right) = 0 \\ \end{gathered} $

В частности, видно, что на самом деле второе уравнение является следствием первого, если к нему применить оператор градиента. Другими словами, в пределе при $\delta \to 0$, вектор ${{B}_{k}}$ стремится к градиенту $\nabla \alpha $ в старших членах:

(3.29)
${{B}_{{k,2}}} = \frac{{\partial {{\alpha }_{1}}}}{{\partial {{x}_{k}}}},$
если так было в начальный момент времени.

Закон Юнга–Лапласа как асимптотический предел. Рассмотрим границу раздела между двумя не смешивающимися фазами. Предположим, что в начальных условиях выполняется равенство ${{B}_{k}} = \delta \frac{{\partial \alpha }}{{\partial {{x}_{k}}}}$. Далее предположим, что такая среда находится в равновесии ($d{\text{/}}dt = 0$ и ${{u}_{k}} = 0$), тогда первое из определяющих уравнений (3.28) превращается в

(3.30)
${{\rho }_{0}}{{E}_{{\alpha ,1}}} = \frac{{\partial {{\rho }_{0}}{{E}_{{{{B}_{k}},0}}}}}{{\partial {{x}_{k}}}}$

Так, если энергия взята в форме (3.14), последнее равенство перепишется в виде

(3.31)
${{\rho }_{0}}{{E}_{{\alpha ,1}}} = \Sigma \frac{{\partial {{B}_{{k,0}}}}}{{\partial {{x}_{k}}}},$
и далее, учитывая, что для начальных данных ${{B}_{{k,0}}} = \delta \frac{{\partial {{\alpha }_{0}}}}{{\partial {{x}_{k}}}}$, получаем
(3.32)
${{p}_{{2,1}}} - {{p}_{{1,1}}} = \delta \Sigma \frac{\partial }{{\partial {{x}_{k}}}}\frac{{\partial {{\alpha }_{0}}}}{{\partial {{x}_{k}}}},$
где ${{p}_{{1,1}}}$ и ${{p}_{{2,1}}}$ – первые члены в разложении давлений фаз ${{p}_{1}}$ и${{p}_{2}}$ по степеням $\delta $:

(3.33)
${{p}_{1}} = {{p}_{{1,0}}} + \delta {{p}_{{1,1}}} + \ldots ,\quad {{p}_{2}} = {{p}_{{2,0}}} + \delta {{p}_{{2,1}}} + \ldots $

Если учесть интерпретацию (3.9) безразмерного параметра $\delta \sim L{\text{/}}\ell $ и сравнить уравнение (3.32) с континуальным аналогом закона Юнга–Лапласа, см., напр. [3],

(3.34)
$\begin{array}{*{20}{r}} {\left[ p \right] = \sigma \nabla \cdot {\mathbf{n}},\quad {\mathbf{n}} = \frac{{\nabla \alpha }}{{\left\| {\nabla \alpha } \right\|}},} \end{array}$
где $\sigma $ – коэффициент поверхностного натяжения, $\left[ p \right]$ – скачок давления через поверхность раздела, ${\mathbf{n}}$ – единичная нормаль к этой поверхности, то можно заключить, что для согласованности (3.32) и (3.34), необходимо выбирать капиллярный коэффициент $\Sigma $ как
(3.35)
$\sigma = \frac{\Sigma }{\ell },$
а масштаб $L$ надо ассоциировать, например, с масштабом $\Delta x$ вычислительной сетки. При этом подразумевается, что $L < \ell $, поскольку параметр $\delta $ предполагается малым.

Мгновенная релаксация давлений. Если рассмотреть далее предел бесконечно быстрой релаксации давлений фаз, то есть $\lambda \to 0$ в (2.4) и (3.1), то из первого уравнения (3.28), записанного в виде

(3.36)
${{\rho }_{0}}\frac{{d{{\alpha }_{1}}}}{{dt}} = - \frac{1}{\lambda }\left( {{{\rho }_{0}}{{E}_{{\alpha ,1}}} - \frac{\partial }{{\partial {{x}_{k}}}}\left( {{{\rho }_{0}}{{E}_{{{{b}_{{k,0}}}}}}} \right)} \right),$
следует закон Юнга–Лапласа (3.30)–(3.32):

(3.37)
${{\rho }_{0}}{{E}_{{\alpha ,1}}} = \frac{\partial }{{\partial {{x}_{k}}}}\left( {{{\rho }_{0}}{{E}_{{{{b}_{{k,0}}}}}}} \right)$

3.3. Характеристические скорости в среде с поверхностным натяжением. Построенная система принадлежит к классу симметрических гиперболических систем, если потенциал энергии $E + {{u}_{l}}{{u}_{l}}{\text{/}}2$ является выпуклой функций параметров состояния. Тем не менее, для построения численных методов и теоретического анализа полезно понимать структуру волн в такой гиперболической системе. Здесь мы приведем только простейшие сведения о характеристических скоростях данной модели. В общем случае, посчитать явно характеристические скорости системы (3.11) вряд ли возможно, и нам удалось найти скорости характеристик только в состоянии равновесия (${{u}_{l}} = 0$, ${{w}_{l}} = 0$, ${{B}_{l}} = 0$), которые оказываются вещественными и вычисляются по формулам

(3.38)
$ \pm {\kern 1pt} \frac{{G\sqrt \Sigma }}{\rho },\quad \pm {\kern 1pt} \sqrt {\frac{{\partial {{p}_{1}}}}{{\partial {{\rho }_{1}}}}} ,\quad \pm {\kern 1pt} \sqrt {\frac{{\partial {{p}_{2}}}}{{\partial {{\rho }_{2}}}}} $

Таким образом, в дополнение к обычным характеристическим скоростям, связанным с компонентами смеси, система обладает характеристиками, связанных с транспортом эффектов поверхностного натяжения, и скорости которых вычисляются по двум новым модулям как $ \pm G\sqrt \Sigma {\text{/}}\rho $.

Отметим также, что в окрестности состояния равновесия, система обладает полным набором собственных векторов.

В общем случае, отличающемся от равновесия, характеристические скорости являются корнями бикубического уравнения, явное решение которого найти не удается.

3.4. Симметрическая форма определяющих уравнений. Определяющие уравнения, сформулированные в предыдущем параграфе, как было показано, являются термодинамически согласованными. Их можно также привести к симметрической форме переписав их в терминах порождающих переменных (3.12). Следуя формализму симметрических гиперболических систем [7, 8], введем порождающий потенциал $L$, производные которого по порождающим переменным являются консервативными переменными уравнений системы (3.12):

$\begin{gathered} {{L}_{{{{q}_{1}}}}} = \rho ,\quad {{L}_{{{{u}_{l}}}}} = \rho {{u}_{l}},\quad {{L}_{n}} = \rho c,\quad {{L}_{{{{z}_{k}}}}} = {{w}_{k}},\quad {{L}_{{{{q}_{2}}}}} = \rho \alpha , \\ {{L}_{r}} = \rho D,\quad {{L}_{{{{j}_{k}}}}} = {{B}_{k}},\quad {{L}_{\theta }} = \rho S \\ \end{gathered} $

Сам потенциал вычисляется с помощью закона сохранения энергии по формуле ${{q}_{1}}{{L}_{{{{q}_{1}}}}} + {{u}_{l}}{{L}_{{{{u}_{l}}}}} + n{{L}_{n}}$ + ${{z}_{k}}{{L}_{{{{z}_{k}}}}} + {{q}_{2}}{{L}_{{{{q}_{2}}}}} + r{{L}_{r}}$ + ${{j}_{k}}{{L}_{{{{j}_{k}}}}} + \theta {{L}_{\theta }} - L$ = ρ$\left( {E + {{u}_{l}}{{u}_{l}}{\text{/}}2} \right)$, откуда получаем $L = {{\rho }^{2}}{{E}_{\rho }}$ + ${{z}_{k}}{{L}_{{{{z}_{k}}}}}$ + ${{j}_{k}}{{L}_{{{{j}_{k}}}}}$. С использованием определенного выше порождающего термодинамического потенциала и порождающих переменных (3.12) потоки системы (3.11) для каждого уравнения запишутся соответственно в виде $\rho {{u}_{k}} = {{\left( {{{u}_{k}}L} \right)}_{{{{q}_{1}}}}}$, $\rho {{u}_{l}}{{u}_{k}}$ + + ${{\rho }^{2}}{{E}_{\rho }}{{\delta }_{{lk}}}$ + $\rho {{w}_{l}}{{E}_{{{{w}_{k}}}}}$ + $\rho {{b}_{l}}{{E}_{{{{b}_{k}}}}}$ = ${{\left( {{{u}_{k}}L} \right)}_{{{{u}_{l}}}}}$ + ${{z}_{k}}{{L}_{{{{z}_{k}}}}}$ + ${{j}_{k}}{{L}_{{{{j}_{k}}}}}$${{\delta }_{{lk}}}{{z}_{\alpha }}{{L}_{{{{z}_{\alpha }}}}}$${{\delta }_{{lk}}}{{j}_{\alpha }}{{L}_{{{{j}_{\alpha }}}}}$, $\rho {{u}_{k}}c + \rho {{E}_{{{{w}_{k}}}}}$ = = ${{\left( {{{u}_{k}}L} \right)}_{n}} + {{z}_{k}}$, ${{u}_{l}}{{w}_{l}} + {{E}_{{{{c}_{1}}}}}$ = ${{u}_{l}}{{L}_{{{{z}_{l}}}}} + n$, $\rho {{\alpha }_{1}}{{u}_{k}}$ = ${{\left( {{{u}_{k}}L} \right)}_{{{{q}_{2}}}}}$, $\rho {{u}_{k}}D + \rho {{E}_{{{{B}_{k}}}}}$ = ${{\left( {{{u}_{k}}L} \right)}_{r}} + {{j}_{k}}$, ${{B}_{l}}{{u}_{l}}$ + + ED = ${{u}_{l}}{{L}_{{{{j}_{l}}}}} + r$, $\rho S{{u}_{k}}$ = ${{\left( {{{u}_{k}}L} \right)}_{\theta }}$.

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

$\frac{{\partial {{L}_{{{{q}_{i}}}}}}}{{\partial t}} + \frac{\partial }{{\partial {{x}_{k}}}}{{\left( {{{u}_{k}}L} \right)}_{{{{q}_{i}}}}} = 0;\quad i = 1,2$
$\frac{{\partial {{L}_{{{{u}_{l}}}}}}}{{\partial t}} + \frac{\partial }{{\partial {{x}_{k}}}}{{\left( {{{u}_{k}}L} \right)}_{{{{u}_{l}}}}} + {{L}_{{{{z}_{l}}}}}\frac{{\partial {{z}_{k}}}}{{\partial {{x}_{k}}}} - {{L}_{{{{z}_{\alpha }}}}}\frac{{\partial {{z}_{\alpha }}}}{{\partial {{x}_{l}}}} + {{L}_{{{{j}_{l}}}}}\frac{{\partial {{j}_{k}}}}{{\partial {{x}_{k}}}} - {{L}_{{{{j}_{\alpha }}}}}\frac{{\partial {{j}_{\alpha }}}}{{\partial {{x}_{l}}}} = 0$
$\frac{{\partial {{L}_{n}}}}{{\partial t}} + \frac{\partial }{{\partial {{x}_{k}}}}{{\left( {{{u}_{k}}L} \right)}_{n}} + \frac{{\partial {{z}_{k}}}}{{\partial {{x}_{k}}}} = 0$
(3.39)
$\frac{{\partial {{L}_{{{{z}_{l}}}}}}}{{\partial t}} + \frac{\partial }{{\partial {{x}_{k}}}}{{\left( {{{u}_{k}}L} \right)}_{{{{z}_{l}}}}} + {{L}_{{{{z}_{\alpha }}}}}\frac{{\partial {{u}_{\alpha }}}}{{\partial {{x}_{l}}}} - {{L}_{{{{z}_{l}}}}}\frac{{\partial {{u}_{k}}}}{{\partial {{x}_{k}}}} + \frac{{\partial n}}{{\partial {{x}_{l}}}} = 0$
$\frac{{\partial {{L}_{r}}}}{{\partial t}} + \frac{\partial }{{\partial {{x}_{k}}}}{{\left( {{{u}_{k}}L} \right)}_{r}} + \frac{{\partial {{j}_{k}}}}{{\partial {{x}_{k}}}} = 0$
$\frac{{\partial {{L}_{{{{j}_{l}}}}}}}{{\partial t}} + \frac{\partial }{{\partial {{x}_{k}}}}{{\left( {{{u}_{k}}L} \right)}_{{{{j}_{l}}}}} + {{L}_{{{{j}_{\alpha }}}}}\frac{{\partial {{u}_{\alpha }}}}{{\partial {{x}_{l}}}} - {{L}_{{{{j}_{l}}}}}\frac{{\partial {{u}_{k}}}}{{\partial {{x}_{k}}}} + \frac{{\partial r}}{{\partial {{x}_{l}}}} = 0$
$\frac{{\partial {{L}_{\theta }}}}{{\partial t}} + \frac{\partial }{{\partial {{x}_{k}}}}{{\left( {{{u}_{k}}L} \right)}_{\theta }} = 0$

Можно показать, что система (3.39) может быть представлена как система квазилинейных уравнений с симметричными матрицами при временной и пространственных производных. Действительно, матрица при временной производной образована вторыми производными потенциала $L$ по порождающим переменным, матрица при пространственной производной по координате ${{x}_{k}}$ образована вторыми производными от ${{u}_{k}}L$ по порождающим переменным, а остальные слагаемые уравнений также записываются симметричными матрицами.

При условии выпуклости потенциала $L$ система (3.39) является симметрической гиперболической по Фридрихсу [7, 8, 15, 16]. Отметим, что хотя выпуклость, по-видимому, имеет место, строгое доказательство этого для используемого в работе выбора внутренней энергии отсутствует, и является предметом дальнейшего исследования.

Заключение. Таким образом, в работе предложена новая модель течения несмешивающихся сжимаемых жидкостей с поверхностным натяжением на границе раздела фаз, определяющие уравнения которой удовлетворяют законам термодинамики (сохранение энергии и возрастание энтропии) и могут быть представлены в виде симметрической гиперболической системы. Модель является расширением разработанной ранее СГТС-модели сжимаемых двухфазных течений [11] за счет введения новых переменных состояния – векторного поля, близкого по смыслу градиенту объемной доли, и дополнительной парной скалярной переменной, описывающей инерционные свойства поверхностей раздела. Показано, что в приближении малости параметров, входящих в кинетические члены определяющих уравнений, модель описывает на континуальном уровне закон Юнга–Лапласа для капиллярного давления. Дальнейшие исследования будут включать в себя продолжение анализа уравнений модели и численное решение тестовых задач для течений несмешивающихся жидкостей.

Исследования Е.И. Роменского в области разработки определяющих уравнений (разд. 2, 3.1, 3.4) выполнены в рамках государственного задания Института математики им. С.Л. Соболева СО РАН (проект № FWNF-2022-0008), а в области асимптотических разложений (разд. 3.2) выполнены при поддержке Российского научного фонда (грант № 22-11-00104).

Список литературы

  1. Brackbill J.U., Kothe D.B., Zemach C. A continuum method for modeling surface tension // J. Comput. Phys. 1992. V. 100. № 2. P. 335–354.

  2. Perigaud G., Saurel R. A compressible flow model with capillary effects // J. Comput. Phys. 2005. V. 209. № 1. P. 139–178.

  3. Popinet S. Numerical models of surface tension // Annu. Rev. Fluid Mech. 2018. V. 50. № 1. P. 49–75.

  4. Schmidmayer K., Petitpas F., Daniel E., Favrie N., Gavrilyuk S. A model and numerical method for compressible flows with capillary effects // J. Comput. Phys. 2017. V. 334. P. 468–496.

  5. Chiocchetti S., Peshkov I., Gavrilyuk S., Dumbser M. High order ADER schemes and GLM curl cleaning for a first order hyperbolic formulation of compressible flow with surface tension // J. Comput. Phys. 2021. V. 426. P. 109898.

  6. Chiocchetti S., Dumbser M. An exactly curl-free staggered semi-implicit finite volume scheme for a first order hyperbolic model of viscous two-phase flows with surface tension // J. Sci. Comput. 2022. V. 94. P. 24.

  7. Godunov S.K., Romenskii E.I. Elements of Continuum Mechanics and Conservation Laws. Springer US, 2003.

  8. Peshkov I., Pavelka M., Romenski E., Grmela M. Continuum mechanics and thermodynamics in the Hamilton and the Godunov-type formulations // Contin. Mech.&Thermodyn. 2018. V. 30. № 6. P. 1343–1378.

  9. Romenski E., Belozerov A.A., Peshkov I.M. Conservative formulation for compressible multiphase flows // Quart. Appl. Math. 2016. V. 74. P. 113–136.

  10. Romenski E., Reshetova G., Peshkov I. Two-phase hyperbolic model for porous media saturated with a viscous fluid and its application to wavefields simulation // Appl. Math. Model. 2022. V. 106. P. 567–600.

  11. Romenski E., Drikakis D., Toro E. Conservative models and numerical methods for compressible two-phase flow // J. Sci. Comput. 2010. V. 42. № 1. P. 68–95.

  12. Romenski E., Resnyansky A.D., Toro E.F. Conservative hyperbolic formulation for compressible two-phase flow with different phase pressures and temperatures // Quart. Appl. Math. 2007. V. 65. № 2. P. 259–279.

  13. Годунов С.К., Михайлова Т.Ю., Роменский Е.И. Системы термодинамически согласованных законов сохранения, инвариантных относительно вращений // Сиб. матем. ж. 1996. Т. 37. № 4. С. 790–806.

  14. Годунов С.К., Роменский Е.И. Элементы механики сплошных сред и законы сохранения. Новосибирск: Науч. книга. 1998.

  15. Friedrichs K.O. Symmetric positive linear differential equations // Commun. on Pure&Appl. Math. 1958. V. 11. № 3. P. 333–418.

  16. Dafermos K.M. Hyperbolic Conservation Laws in Continuum Physics. Berlin: Springer, 2016.

  17. Dhaouadi F., Dumbser M. A first order hyperbolic reformulation of the Navier–Stokes–Korteweg system based on the GPR model and an augmented Lagrangian approach // J. Comput. Phys. 2022. V. 470. P. 111544.

  18. Dhaouadi F., Gavrilyuk S., Vila J.-P. Hyperbolic relaxation models for thin films down an inclined plane // Appl. Math.&Comput. 2022. V. 433 P. 127378.

Дополнительные материалы отсутствуют.