引言

气液流体相变是自然界和化学化工领域常见的现象。van der Waals流体方程组是描述一类气液流体相变的模型，其分析与计算在流体动力学中占有重要的地位。在数值模拟时，如果按照原始模型计算，由于气液流体在相变区的不稳定性，模型的解会出现剧烈的振荡。为解决这一问题，通常在连续性方程中加入人工黏性项以抑制振荡。近年来，对于带人工黏性的非理想流体，其初边值问题解的性态分析有不少研究。文献[1-3]讨论了类van der Waals流体周期边界下解的存在性、有界性和渐近稳定性；文献[4-8]讨论了van der Waals流体多种情况下解的渐近稳定性。

以上文献中的黏性系数一般取为常数，但是在分析和计算时，人工黏性应主要在不稳定的椭圆区(相变区)起作用，而当人工黏性系数为常数时，在整个求解区域内都用到了同样程度的人工黏性，这就意味着在稳定的非椭圆区中加入了多余的黏性，从而增大了非椭圆区内的误差。针对这一问题，本文将人工黏性系数取为依赖于质量密度(或比容)的函数，使密度在非椭圆区时的人工黏性较小，在不稳定的椭圆区时的人工黏性增大，以减缓震荡，从而在一定程度上减少误差，提高数值模拟的准确性。

本文针对Lagrange坐标下具有变人工黏性系数的一维van der Waals等温模型的周期初边值问题，通过构造泛函进行相关估计，给出了其解只有平凡解和存在非平凡解且容易判别的充分条件，并对理论结果进行了相关数值验算，同时也对初边值问题进行了计算。分析和计算结果表明，即使非相变区域的人工黏性系数很小，也可得到与人工黏性系数较大的常数时类似的抑制相变区解震荡的效果。

1 模型的建立与分析

考虑如下Lagrange坐标系下的一维van der Waals等温流体流动模型的周期初边值问题

$ \left\{ {\begin{array}{*{20}{l}} {{v_t} - {u_x} = {{\left( {{\varepsilon _1}(v){v_x}} \right)}_x}}\\ {{u_t} - \sigma {{(v)}_x} = {\varepsilon _2}{u_{xx}}} \end{array}\;\;\;\;(x,t) \in {\bf{R}} \times {{\bf{R}}_ + }} \right. $

(1)

$ {\left. {(v,u)} \right|_{t = 0}} = \left( {{v_0},{u_0}} \right)(x),x \in ( - \infty ,\infty ) $

(2)

$ \begin{array}{l} v,u(x,t) = (v,u)(x + 2L,t),(x,t) \in ( - \infty ,\\ \infty ) \times [0,\infty ) \end{array} $

(3)

式中，u为流速，ρ为质量密度，v=1/ρ为比容, σ(v)为压强函数，ε₁(v)>0为依赖于比容的人工黏性系数，ε₂>0为流体黏性系数。

在van der Waals流体模型中，σ(v)的函数关系为$ -\sigma(v)=\frac{R \theta}{v-b}-\frac{a}{v^{2}}$，其中R, θ, a, b都是正常数，且v>b。在有气液相变时需满足$ \frac{R \theta b}{a}<\frac{8}{27} $，此时σ(v)有两个临界点v₁和v₂(v₁＜v₂)，并有σ′(v₁)=σ′(v₂)=0。当v∈(b, v₁)∪(v₂, ∞)时，σ′(v)>0，当v∈(v₁, v₂)时，σ′(v)＜0。(b, v₁)为液态区域，(v₂, ∞)区域为气态区域，(v₁, v₂)为气液共存区域，即相变区域。在数学上，(b, v₁)∪(v₂, ∞)是方程(2)的稳定的非椭圆区，而(v₁, v₂)是不稳定的椭圆区。

利用周期边界条件，对式(1)在v₀(x)=v₀(x+2L), u₀(x)=u₀(x+2L)上积分，可以得到

$ \int_0^{2L} {v(x,t){\rm{d}}x} = \int_0^{2L} {{v_0}(x){\rm{d}}x} $

$ \int_0^{2L} {u(x,t){\rm{d}}x} = \int_0^{2L} {{u_0}(x){\rm{d}}x} $

令

$ \left\{ {\begin{array}{*{20}{l}} {{m_0} = \frac{1}{{2L}}\int_0^{2L} {{v_0}} (x){\rm{d}}x}\\ {{m_1} = \frac{1}{{2L}}\int_0^{2L} {{u_0}} (x){\rm{d}}x} \end{array}} \right. $

则

$ \left\{ {\begin{array}{*{20}{l}} {\int_0^{2L} {\left[ {v(x,t) - {m_0}} \right]{\rm{d}}x} = 0}\\ {\int_0^{2L} {\left[ {u(x,t) - {m_1}} \right]{\rm{d}}x} = 0} \end{array}} \right. $

相应地，当u, v与t无关时，有以下稳态的周期边值问题

$ \left\{ {\begin{array}{*{20}{l}} { - {u_x} = {{\left( {{\varepsilon _1}(v){v_x}} \right)}_x}}\\ { - \sigma {{(v)}_x} = {\varepsilon _2}{u_{xx}}}\\ {(v,u)(x) = (v,u)(x + 2L)}\\ {\frac{1}{{2L}}\int_0^{2L} v (x){\rm{d}}x = {m_0}}\\ {\frac{1}{{2L}}\int_0^{2L} u (x){\rm{d}}x = {m_1}} \end{array}} \right. $

(4)

对式(4)中第一、二个方程在[0, x]上积分，得到

$ \left\{ {\begin{array}{*{20}{l}} { - u = {\varepsilon _1}(v){v_x} + {c_1}}\\ { - \sigma (v) = {\varepsilon _2}{u_x} + {c_2}} \end{array}} \right. $

(5)

式中，$ c_{1}=-\frac{1}{2 L} \int_{0}^{2 L} u(x) \mathrm{d} x, c_{2}=-\frac{1}{2 L} \int_{0}^{2 L} \sigma(v(x)) \mathrm{d} x $。

再将式(4)第一个方程代入式(5)第二个方程，得到边值问题

$ \left\{ {\begin{array}{*{20}{l}} {{\varepsilon _2}{{\left( {{\varepsilon _1}(v){v_x}} \right)}_x} = \sigma (v) - c}\\ {v(x) = v(x + 2L)}\\ {\int_0^{2L} v (x){\rm{d}}x/2L = {m_0}} \end{array}} \right. $

(6)

式中，$ c = \frac{1}{{2L}}\int_0^{2L} \sigma (v(x)){\rm{d}}s $。

假设V(x)=v(x)－m₀，则定解问题式(6)可写成

$ \left\{ {\begin{array}{*{20}{l}} {{\varepsilon _2}{{\left( {{\varepsilon _1}(v){V_x}} \right)}_x} = \bar \sigma (V) - \bar a}\\ {V(x) = V(x + 2L)}\\ {\int_0^{2L} V (x){\rm{d}}x = 0} \end{array}} \right. $

(7)

式中, σ(V)=σ(V+m₀)－σ(m₀)，$ \bar a = \frac{1}{{2L}}\int_0^{2L} {\bar \sigma } (V(x)){\rm{d}}x $。

定义周期的Sobolev空间为

$ H_{{\rm{per}},0}^1 = \left\{ {v|v \in H_{{\rm{per}}}^1(R),\int_0^{2L} v (y){\rm{d}}y = 0} \right\} $

并定义Sobolev空间上的一个泛函

$ G(v) = \int_0^{2L} {\left[ {\frac{{{\varepsilon _1}\left( {v + {m_0}} \right){\varepsilon _2}}}{2}v_x^2 + H(v)} \right]{\rm{d}}x} $

(8)

式中

$ \begin{array}{l} H(v) = \int_0^v {\bar \sigma } (s){\rm{d}}s = - \frac{{a{v^2}}}{{m_0^2\left( {v + {m_0}} \right)}} + R\theta \left( {\frac{v}{{{m_0} - b}} - } \right.\\ \left. {\ln \left( {1 + \frac{v}{{{m_0} - b}}} \right)} \right) \end{array} $

参照文献[1]，有以下引理。

引理1   V是定解问题式(7)的解，等价于V是G(v)在H_{per, 0}¹上的临界点。

引理2   若V∈H_{per, 0}¹是式(7)的解，则

$ {\left\| V \right\|_{{L^2}}} \le \frac{L}{{\rm{ \mathsf{ π} }}}{\left\| {{V_x}} \right\|_{{L^2}}} $

引理1、2的证明可参照文献[1]，故此处从略。由引理1、2可以得到定理1和定理2。

定理1

若$ B = \frac{{\varepsilon _1^*{\varepsilon _2}{{\rm{ \mathsf{ π} }}^2}}}{{{L^2}}} + \frac{1}{b}\left[ {\frac{{2R\theta }}{{{m_0} - b}} - \left( {\frac{a}{{m_0^2}} + \frac{a}{{{m_0}b}}} \right)} \right] > 0 $，$ \varepsilon _1^* = \mathop {\inf }\limits_{v > b} {\varepsilon _1}(v) $，则式(7)只有平凡解V(x)≡0，即稳态问题式(4)只有解v(x)≡m₀和u(x)≡m₁。

证明：因为

$ \begin{array}{l} 0 = \int_0^{2L} {\left[ { - {\varepsilon _2}{{\left( {{\varepsilon _1}(v){V_x}} \right)}_x} + \bar \sigma (V) - \bar a} \right]V{\rm{d}}x} = \\ \int_0^{2L} {\left( {{\varepsilon _2}{\varepsilon _1}(v)V_x^2 + \bar \sigma (V)V} \right){\rm{d}}x} \end{array} $

又因为

$ \begin{array}{l} \bar \sigma (V)V = \left( {\sigma \left( {V + {m_0}} \right) - \sigma \left( {{m_0}} \right)} \right)V = \\ \left( {\frac{{R\theta }}{{\left( {{m_0} - b} \right)(v - b)}} - \frac{{a\left( {v + {m_0}} \right)}}{{{v^2}m_0^2}}} \right){V^2} \ge m\left( {{m_0}} \right){V^2} \end{array} $

式中，$ m\left(m_{0}\right)=\inf \limits_{v>b}\left(\frac{R \theta}{\left(m_{0}-b\right)(v-b)}-\frac{a\left(v+m_{0}\right)}{v^{2} m^{2}}\right) $。

根据引理1、2和文献[1]中引理2.3的证明可得，当$ \frac{{\varepsilon _1^*{\varepsilon _2}{{\rm{ \mathsf{ π} }}^2}}}{{{L^2}}} + m\left( {{m_0}} \right) > 0 $时，式(7)只有平凡解V(x)≡0。

对m(m₀)进行估计。

令$ g(v)=\frac{R \theta}{\left(m_{0}-b\right)(v-b)}-\frac{a\left(v+m_{0}\right)}{v^{2} m_{0}^{2}} $，设$ \rho = \frac{1}{v} $，则$ 0 < \rho < \frac{1}{b} $，且有

$ g(v) = h(\rho ) = \rho \left( {\frac{{R\theta }}{{{m_0} - b}} \times \frac{1}{{1 - b\rho }} - \frac{a}{{m_0^2}} - \frac{a}{{{m_0}}}\rho } \right) $

因为$ \frac{R \theta}{m_{0}-b}-\frac{a}{m_{0}^{2}}>0, \frac{1}{1-b \rho}>1+b \rho $，所以有

$ \begin{array}{l} h(\rho ) > \rho \left( {\frac{{R\theta }}{{{m_0} - b}} - \frac{a}{{m_0^2}} + \left( {\frac{{R\theta b}}{{{m_0} - b}} - \frac{a}{{{m_0}}}} \right)\rho } \right) \buildrel \Delta \over = \\ q(\rho ) \end{array} $

若$ \frac{{R\theta b}}{{{m_0} - b}} - \frac{a}{{{m_0}}} \ge 0 $，即$ {m_0} \le \frac{{ab}}{{a - R\theta b}} $，则有h(ρ)>0，同时有m(m₀)≥0，此时式(7)只有平凡解V(x)≡0。

若$ \frac{{R\theta b}}{{{m_0} - b}} - \frac{a}{{{m_0}}} < 0 $，则当$ 0 < \rho < \frac{1}{b} $时，有

$ \begin{array}{l} q(\rho ) \ge \min \left\{ {0,\frac{1}{b}\left[ {\frac{{R\theta }}{{{m_0} - b}} - \frac{a}{{m_0^2}} + \frac{1}{b}\left( {\frac{{R\theta b}}{{{m_0} - b}} - } \right.} \right.} \right.\\ \left. {\left. {\left. {\frac{a}{{{m_0}}}} \right)} \right]} \right\} = w\left( {{m_0}} \right) \end{array} $

进而有m(m₀)≥w(m₀)，即当

$ \frac{{\varepsilon _1^*{\varepsilon _2}{{\rm{ \mathsf{ π} }}^2}}}{{{L^2}}} + \frac{1}{b}\left[ {\frac{{2R\theta }}{{{m_0} - b}} - \left( {\frac{a}{{m_0^2}} + \frac{a}{{{m_0}b}}} \right)} \right] > 0 $时，式(7)只有平凡解V(x)≡0。综上，定理1得证。

定理2   若$ C = \frac{{{\varepsilon _1}\left( {{m_0}} \right){\varepsilon _2}{{\rm{ \mathsf{ π} }}^2}}}{{2{L^2}}} - \frac{a}{{m_0^3}} + \frac{{R\theta }}{{2{{\left( {{m_0} - b} \right)}^2}}} < 0 $，则式(7)除了有平凡解V(x)≡0外，还有非平凡解，即稳态问题式(4)还有非常数解。

证明：因为V是式(7)的解，则V是泛函G(v)在H_{per, 0}¹上的临界点，V(x)≡0对应于G(0)=0。若在一定条件下存在v_*(x)∈H_{per, 0}¹使得G(v_*)＜0，则式(7)至少存在一个非平凡解。

取$ {v_*} = \eta \frac{1}{{\sqrt L }}\sin \left( {\frac{{{\rm{ \mathsf{ π} }}x}}{L}} \right) $，其中η>0是一个小正数，则v_*∈H_{per, 0}¹。记$ {v_1}(x) = \frac{1}{{\sqrt L }}\sin \left( {\frac{{{\rm{ \mathsf{ π} }}x}}{L}} \right) $，则有

$ \left\{ {\begin{array}{*{20}{l}} {\int_0^{2L} {v_1^\prime } {{(x)}^2}{\rm{d}}x = \frac{{{{\rm{ \mathsf{ π} }}^2}}}{{{L^2}}}}\\ {\int_0^{2L} {v_1^2} (x){\rm{d}}x = 1}\\ {\int_0^{2L} {v_1^3} (x){\rm{d}}x = 0} \end{array}} \right. $

因为$ \mathop {\lim }\limits_{\eta \to {0^ + }} {\varepsilon _1}\left( {{v_*} + {m_0}} \right) = {\varepsilon _1}\left( {{m_0}} \right) $，所以对于任意给定的大整数N存在η_N>0，使得η＜η_N时有$ {\varepsilon _1}\left( {{v_*} + {m_0}} \right) < {\varepsilon _1}\left( {{m_0}} \right) + \frac{1}{N} $。同时，当η>0充分小且η＜η_N时有$ \left| {\frac{{{v_*}}}{{{m_0} - b}}} \right| < 1 $，此时有

$ \left\{ {\begin{array}{*{20}{l}} {\frac{1}{{{v_*} + {m_0}}} = \frac{1}{{{m_0}}} \times \frac{1}{{1 + \frac{{{v_*}}}{{{m_0}}}}} \ge \frac{1}{{{m_0}}}\left( {1 - \frac{{{v_*}}}{{{m_0}}}} \right)}\\ {\ln \left( {1 + \frac{{{v_*}}}{{{m_0} - b}}} \right) \ge \frac{{{v_*}}}{{{m_0} - b}} - \frac{1}{2}{{\left( {\frac{{{v_*}}}{{{m_0} - b}}} \right)}^2}} \end{array}} \right. $

由此可得

$ H\left( {{v_*}} \right) \le - \frac{a}{{m_0^3}}v_*^2 + \frac{a}{{m_0^4}}v_*^3 + \frac{{R\theta }}{{2\left( {{m_0} - b} \right)}}v_*^2 $

当η>0充分小时，有

$ \begin{array}{l} G\left( {{v_*}} \right) = \int_0^{2L} {\left[ {\frac{{{\varepsilon _1}\left( {{v_*} + {m_0}} \right){\varepsilon _2}}}{2}{{v'}_ * }{{(x)}^2} + } \right.} \\ \left. {H\left( {{v_*}} \right)} \right]{\rm{d}}x \le \left[ {\frac{{\left( {{\varepsilon _1}\left( {{m_0}} \right) + \frac{1}{N}} \right){\varepsilon _2}{{\rm{ \mathsf{ π} }}^2}}}{{2{L^2}}} - \frac{a}{{m_0^3}} + } \right.\\ \left. {\frac{{R\theta }}{{2{{\left( {{m_0} - b} \right)}^2}}}} \right]{\eta ^2} \end{array} $

所以当$ \frac{{{\varepsilon _1}\left( {{m_0}} \right){\varepsilon _2}{{\rm{ \mathsf{ π} }}^2}}}{{2{L^2}}} - \frac{a}{{m_0^3}} + \frac{{R\theta }}{{2{{\left( {{m_0} - b} \right)}^2}}} < 0 $时，只要η>0充分小，就有G(v_*)＜0，即式(7)至少存在一个非平凡解。定理2得证。

2 数值模拟

先利用算例由边值问题式(6)佐证定理1和定理2。

在式(6)中，令$ q(x) = \int_0^x v (s){\rm{d}}s, w(x) = {\varepsilon _1}(v){\mathit{v}_x} $，则式(6)化为

$ \left\{ \begin{array}{l} {v_x} = \frac{w}{{{\varepsilon _1}(v)}}\\ {w_x} = \frac{{\sigma (v) - c}}{{{\varepsilon _2}}}\\ {q_x} = v\\ v(0) = v(2L),w(0) = w(2L),q(0) = 0,\\ \;\;\;\;\;q(2L) = 2L{m_0} \end{array} \right. $

(9)

式(9)即为含未知参数c的两点边值问题，可用MATLAB的bvp4c命令进行如下数值求解。

令

$ {\varepsilon _1}(v) = 0.1{{\rm{e}}^{ - 5{{(v - 1)}^2}}} $

(10)

同时取$ a = 3, b = \frac{1}{3}, R = \frac{8}{3}, \theta = 0.9, {\varepsilon _2} = 0.1, L = 0.5 $，其中a、b、R的值为van der Waals无量纲简化形式的对应值。此时稳定区为(1/3，0.718 6)∪(1.528 5，+∞)，不稳定区为(0.718 6，1.528 5)。σ(v)和ε₁(v)的函数图像如图 1、2所示。

当式(10)成立时，ε₁(v)在不稳定区取值较大，从而会有较好的抑制解震荡的作用；而在稳定区取值较小时，可使解在稳定区误差尽量小。

当m₀=0.489时，B=0.047 565>0时满足定理1的条件。由定理1可知式(9)只有平凡解v≡m₀、w≡0和q≡m₀x。在bvp4c命令中，在合适的范围内取多个不同的初始值时，式(9)的解均为平凡解，如图 3所示。

周期边界条件式(3)可等价转化为[0, 2L]上的周期边界条件

$ \left\{ {\begin{array}{*{20}{l}} {{{\left. {(v,u)} \right|}_{x = 0}} = {{\left. {(v,u)} \right|}_{x = 2L}}}\\ {{{\left. {\left( {{v_x},{u_x}} \right)} \right|}_{x = 0}} = {{\left. {\left( {{v_x},{u_x}} \right)} \right|}_{x = 2L}}} \end{array}} \right. $

(11)

取式(10)的ε₁(v)和参数及m₀=0.489，再取

$ \left\{ {\begin{array}{*{20}{l}} {{v_0}(x) = 0.489 + 0.1\sqrt 2 \sin (2{\rm{ \mathsf{ π} }}x)}\\ {{u_0}(x) = 0.1\cos (2{\rm{ \mathsf{ π} }}x)} \end{array}} \right. $

然后对式(1)、(2)、(11)组成的模型进行离散编程计算，结果如图 4。由图可知，当时间趋于无穷大时，初边值问题式(1)~(3)的解趋于平凡的稳态解。

当m₀=1时，C=-0.102 608＜0，此时参数满足定理2的条件。取如式(12)的初始值时

$ \begin{array}{*{20}{l}} {{v_0}(x) = 1 + 0.2\sqrt 2 \sin (0.7{\rm{ \mathsf{ π} }}x)}\\ {{w_0}(x) = 0.14\sqrt 2 \pi \cos (0.7{\rm{ \mathsf{ π} }}x)}\\ {{q_0}(x) = x + \sin (2{\rm{ \mathsf{ π} }}x)} \end{array} $

(12)

边值问题式(9)有非平凡解，如图 5所示。

经多次类似的数值模拟得出，当式(13)成立时

$ \frac{{\varepsilon _1^*{\varepsilon _2}{{\rm{ \mathsf{ π} }}^2}}}{{{L^2}}} + \frac{1}{b}\left[ {\frac{{2R\theta }}{{{m_0} - b}} - \left( {\frac{a}{{m_0^2}} + \frac{a}{{{m_0}b}}} \right)} \right] > 0 $

(13)

对应的稳态问题只有平凡解。

当式(14)成立时，对应的稳态问题有非平凡解。

$ \frac{{{\varepsilon _1}\left( {{m_0}} \right){\varepsilon _2}{{\rm{ \mathsf{ π} }}^2}}}{{2{L^2}}} - \frac{a}{{m_0^3}} + \frac{{R\theta }}{{2{{\left( {{m_0} - b} \right)}^2}}} < 0 $

(14)

3 实例验证

选用一些计算实例来考察当人工黏性系数ε₁是v的函数并且在不稳定区域取值较大、在稳定区域取值较小时，其对不稳定区域解的震荡的抑制作用，并与ε₁为常数的情形作比较。参数a, b, R, θ, L, ε₂取值同式(10)，取m₀=0.7, ε₁(v)分别为0.49e^－5(v－1)2、0.50e^－5(v－1)2、0.33和0.34，同时设v₀(x)=m₀+0.1sin (2πx), u₀(x)=0.1cos (2πx)，得到初边值问题式(1)~(3)的计算结果如图 6所示。

从图 6中可以看出，当ε₁(v)=0.49e^－5(v－1)2时，初始阶段解有明显的震荡，而当ε₁(v)=0.50e^－5(v－1)2时，不稳定区解的震荡已被抑制而很快趋于稳态的平凡解，这种变化趋势与ε₁(v)=0.34的情形类似，而ε₁(v)=0.50e^－5(v－1)2时在不稳定区取值与0.34相近，在稳定区远小于0.34。其他的算例也有类似的情况。这表明，当人工黏性系数ε₁取为v的函数时，若其在不稳定区域取值较大、而在稳定区域取值较小，则其对不稳定区域解的震荡的抑制作用与在不稳定区域取值相近的常数人工黏性系数类似。