考虑二维区域 $\Omega$内的简单闭合的沉浸边界 $\Gamma$与黏性流体内的耦合运动, 流体的流动由Navier-Stokes方程描述, 而沉浸边界 $\Gamma$假定是由弹性材料组成的薄膜, 以保证其在静态下被拉伸或放松时, 能在流体上产生回复力.若 $\Gamma$是简单闭的, 则其可将区域 $\Omega$分成两个子区域, 曲线外部设为 $\Omega^+$, 曲线内部设为 $\Omega^-$.若 $\Gamma$是多纽结闭的, 仍可同样处理, 只是复杂度增加而已.假定在 $\Omega^+$与 $\Omega^-$内, 液体的属性(如密度、凝固点等)是一样的, 但由于区域之间有相互作用力, 这会使得在 $\Gamma$上产生速度梯度, 且在流体压力下会产生断裂.

流体运动的Navier-Stokes方程为

$ \begin{equation} \begin{cases} \dfrac{\partial \mathrm{\mathbf{u}}}{\partial t}+\mathrm{\mathbf{u}}\cdot\nabla \mathrm{\mathbf{u}}=-\nabla p+\mu\nabla^2 u+\mathrm{\mathbf{F}}, \\ \nabla\cdot \mathrm{\mathbf{u}}=0, \end{cases} \end{equation} $

(1.1)

其中 $\mathrm{\mathbf{u}}=(u_1, u_2)$表示流体的速度, $p$为压强, $\mu$为流体的黏性系数, 一般设为常数; $\mathrm{\mathbf{F}}=(F_1, F_2)$为面际压力, 其为一种弹性张力, 产生于边界闭曲线 $\Gamma$的伸缩运动, 该运动用 $\Gamma$上的实质坐标 $\alpha$来表示.流体的密度设为常数, 为方便计算, 通常取为1个单位.

在计算物理领域, 为了得到精确的数值计算方法, 须在界面上强加某些边界条件或跳转条件.这一点, 目前已得到该领域研究学者的共识.如Fedkiw等在文[1]中设计了一种“虚拟流体”的方法, 以求捕获非黏性Euler方程的接触间断点的边界条件.同年, Fedkiw等在文[2]中将该方法推广应用于一般的可压缩流的间断点, 如电击、爆炸与燃爆等.受“虚拟流体”方法的启发, Liu X D等^[3]提出一种文本边界条件捕获的方法, 用以求解界面上的系数Poisson方程.该方法具有鲁棒性且易在三维空间实施.后来, Liu X D等^[4]在弱表达式上补充了一个正式的收敛性证明.而在文[5]中, “虚拟流体”的方法被进一步推广应用于二相不可压缩流, 以讨论界面附近的不连续量的数值模糊拖尾效应.但是, 文[6-8]认为, 文[5]中所提出的 $\delta$函数模糊技巧甚至会在基本计算中导致 $O(1)$误差.要克服这一缺陷, 或许我们还得回头去参考Leveque等^[9]所提出的沉浸边界的方法.该方法结构复杂, 不过好在Li等^[10]提供了一种简单版本, 该版本与文[3]的边界条件捕获法很相似.

为计算方便, 方程(1.1)一般假定双周期边界条件.在 $t$时刻, 实质点 $\alpha$的当前位置记为 $X(\alpha, t)$, $\Gamma$上的弧长 $l$与张力 $\mathrm{\mathbf{F}}=(F_1, F_2)$满足

$ \begin{equation} F_i(x, t)=\int_0^L f_i(l, t)\delta\left[x-X(\alpha(l), t)\right]\mathrm{d}l, \quad i=1, 2, \end{equation} $

(1.2)

其中 $f_i$是 $l$对应弧段处所受到的压力, $\delta$是二维的 $\delta$函数.由于力的平衡, 张力 $\mathrm{\mathbf{f}}$满足

$ \begin{equation} \mathrm{\mathbf{f}}(l, t)]=\dfrac{\partial}{\partial l}[T(l, t)\tau(l, t)], \end{equation} $

(1.3)

其中物质被拉伸时, 张量 $T(l, t)=T_0\left(\left|\frac{\partial X}{\partial\alpha}\right|-1\right)$, 且系数 $T_0$依赖于界面的弹性属性, 一般假设为常数; $\tau(l, t)$为 $\Gamma$的单位切向量, 即 $\tau(l, t)=\frac{\partial X}{\partial l}=\frac{\partial X/\partial\alpha}{|\partial X/\partial\alpha|}$.由于松弛状态下 $|\partial X/\partial\alpha|=1$, 故此时张量消失.

由于张力 $\mathrm{\mathbf{F}}$的奇性, 方程(1.1)的解相对 $\Gamma$而言是不光滑的.这一效应可用压力 $p$与速度 $u$的跳转条件来表示.在边界上, 对于任意带跳跃间断点的量 $q$, 其跳跃可公式化为

$ \begin{equation} [q(X, t)]=\lim\limits_{\varepsilon\to0^+}q(X+\varepsilon\mathrm{\mathbf{n}}, t)-\lim\limits_{\varepsilon\to0^-}q(X-\varepsilon\mathrm{\mathbf{n}}, t), \end{equation} $

(1.4)

其中 $\mathrm{\mathbf{n}}$为 $\Gamma$的单位法向量.而压力 $p$与速度 $\mathrm{\mathbf{u}}$的跳转条件分别为

$ \begin{equation} \begin{cases} [p]=f\cdot\overrightarrow{n}, \ \left[\dfrac{\partial p}{\partial \overrightarrow{n}}\right]=\dfrac{\partial}{\partial s}(f\cdot\tau), \\[3mm] [u]=0, \ \mu\left[\dfrac{\partial u}{\partial \overrightarrow{n}}\right]=-(f\cdot\tau)\tau. \end{cases} \end{equation} $

(1.5)

显然, 速度 $\mathrm{\mathbf{u}}=(u_1, u_2)$对于 $\Gamma$是连续的, 且 $\Gamma$随着黏性流体一起运动.

Peskin ^[11-13]设计了形如式(1.4) 的沉浸边界方法, 用以解决流体内的耦合运动问题, 并将之推广后应用于生物力学领域.然而, Peskin的方法只具有一阶精度.而文[14, 15]中的沉浸边界方法虽然具有二阶精度, 但其讨论的不是等高的分界面.利用沉浸边界的方法得到较好的精度的应该是文[16, 17], 其中陈述了解的跳转与该跳转关于边界压力的导数, 这一方法被首次应用于Stokes方程(见文献[16]).将方程(1.1) 中的速度 $u$与压强 $p$分解为Stokes分量与正则分量, 分别用指标 $\mathrm{s}$与 $\mathrm{r}$表示, 即

$ \mathrm{\mathbf{u}}=\mathrm{\mathbf{u}}_s+\mathrm{\mathbf{u}}_r, \ p=p_s+p_r, $

则Stokes方程为

$ \begin{equation} \begin{cases} \nabla p_{s}=\mu\nabla^2 \mathrm{\mathbf{u}}_{s}+F, \\ \nabla\cdot \mathrm{\mathbf{u}}_{s}=0. \end{cases} \end{equation} $

(1.6)

可见, Stokes分量在数学上是由Stokes方程决定, 而在物理学中, 其由奇性面际压力决定. Stokes分量可以利用沉浸界面的方法得到, 该方法具有二阶精度.正则分量在数学上由原Navier-Stokes决定, 而在物理学中, 其为Stokes分量所施加的体压力.正则分量可利用时间步长的方法得到.精确描述Stokes流中的面际运动要比在Navier-Stokes流中容易一些, 且Stokes方程的解可用边界积分来表示(见文献[18]), 但其忽略了流体的正则分量.

比较方程(1.1) 与(1.6), 可得正则分量的方程为

$ \begin{equation} \begin{cases} \dfrac{\partial \mathrm{\mathbf{u}}_{r}}{\partial t}+\mathrm{\mathbf{u}}\cdot\nabla \mathrm{\mathbf{u}}_{r}=-\nabla p_{\mathrm{r}}+\mu\nabla^2 \mathrm{\mathbf{u}}_{\mathrm{r}}+\mathrm{\mathbf{F}}_b, \\ \nabla\cdot \mathrm{\mathbf{u}}_{r}=0, \end{cases} \end{equation} $

(1.7)

其中 $\mathrm{\mathbf{F}}_b$为Stokes速度的实体导数所提供的体压力, 即

$ \begin{equation} \mathrm{\mathbf{F}}_b=-\dfrac{\partial u_{s}}{\partial t}-u\cdot\nabla u_{s}. \end{equation} $

(1.8)

本文借助于上述文献的优点, 力图克服其缺点, 而引入一种速度与压强分解的方法, 来求解方程(1.1)-(1.2).该方法在时空上都具有二阶精度.速度与压强分解成两部分, 一部分由稳定的Stokes方程与 $\Gamma$上的面际压力决定, 另一部分也命名为正则部分, 其在面际附近的正则格点上可无需特殊手段地计算.这一分解方法使得流体中变量的间断点可用相对简单的方式来精确描述.此外, 在数值化计算压力时可采用更小的时间步长.至于精度的分析, 本文借鉴文[19]的方法.

选定 $\triangle t>0$, 设 $t^n=n\triangle t$, 表示第 $n$次迭代时间, $n=0, 1, 2, \cdots$.将区域 $\Omega$划分为矩形格子, 网格间距 $h_x, h_y$分别沿着 $x, y$轴计算, 通常可取 $h_x=h_y=h$.设 $i, j=0, 1, 2, \cdots, N$, 则每个格点可以用坐标 $(i\cdot h, j\cdot h)$来表示, 而 $t$时刻沉浸边界的位置可以用一个边界标记集合 $\mathbf{X}_k(t)\equiv (X_k(t), Y_k(t))$来表示, 其中 $k=0, 1, 2, \cdots, N_k$, 且由于 $\Gamma$是简单闭曲线, 则 $\mathbf{X}_0(t)=\mathbf{X}_{N_k}(t)$.第 $k$个边界标记近似为 $\mathbf{X}(\alpha_k, t)$, 且 $\alpha_k=kL_0/N_k$, $L_0$为 $\gamma$在松弛状态的长度.可见, 这些标记是在 $\Gamma$处于松弛状态下被等距间隔的.

Stokes压强 $p_s$、速度 $\mathrm{\mathbf{u}}_s$分别与 $p$、 $u$具有同样跳转条件, 即满足式(1.5).一旦 $\mathrm{\mathbf{u}}_s$可知, 则当前时刻产生于物质的体压力 $\mathrm{\mathbf{F}}_b$与 $\mathrm{\mathbf{u}}_s$的全微分就可求出.这需要用到半Lagrange离散化, 正如文[20]指出, 这一方法一直被应用于流体计算.在半Lagrange离散方案中, 水平对流项被包含于导数中, 且式(1.7) 第一项变为

$ \begin{equation} \dfrac{\mathrm{d}\mathrm{\mathbf{u}}_{r}}{\mathrm{d}t}=-\nabla p_{r}+\mu\nabla^2 \mathrm{\mathbf{u}}_{r}+\mathrm{\mathbf{F}}_b, \end{equation} $

(2.1)

其中 $\mathrm{\mathbf{F}}_b=-\frac{\mathrm{d}u_{s}}{\mathrm{d}t}$.

参照文[20, 21], 采用二阶向后差分, 沿着轨迹将式(2.1) 离散化为

$ \begin{equation} \dfrac{3\mathrm{\mathbf{u}}_{r}^{n+1} -4\widetilde{\mathrm{\mathbf{u}}}_{r}^{n} +\widetilde{\mathrm{\mathbf{u}}}_{r}^{n-1}}{2\triangle t}+\nabla p_{r}^n =\mu\nabla^2\mathrm{\mathbf{u}}_{r}^{n+1}+\mathrm{\mathbf{F}}_b^{n+1}, \end{equation} $

(2.2)

其中 $\widetilde{\mathrm{\mathbf{u}}}_{r}^{n}$与 $\widetilde{\mathrm{\mathbf{u}}}_{r}^{n-1}$分别为上游点 $x^{n}$与 $x^{n-1}$的流动速度, 即

$ \begin{equation} \widetilde{\mathrm{\mathbf{u}}}_{r}^{n} =\mathrm{\mathbf{u}}_{r}(x^{n}, t^{n}), \ \widetilde{\mathrm{\mathbf{u}}}_{r}^{n-1}=\mathrm{\mathbf{u}}_{r}(x^{n-1}, t^{n-1}), \end{equation} $

(2.3)

式(2.1) 右边的 $\mathrm{\mathbf{F}}_b$中还含有 $\mathrm{d}u_{s}$的实体导数, 可类似于式(2.2) 将其离散化.

为了估计式(2.3), 首先得估计坐标 $x^{n}, x^{n-1}$.由于初始问题满足

$ \begin{equation} \dfrac{\mathrm{d}x(t)}{\mathrm{d}t}=\mathrm{\mathbf{u}}_{r}(x, t), \ x(t^{n+1})=x_0, \end{equation} $

(2.4)

对其基于时间段 $[t^{n+1}, t^{n}]$与 $[t^{n+1}, t^{n-1}]$分别向后积分, 可得估计式

$ \begin{align} x^n=x_0-\triangle t\mathrm{\mathbf{u}}_{r}(x^*, t^{n+1/2}),\ x^*=x_0-\dfrac{\triangle t}{2}\mathrm{\mathbf{u}}\left(x_0-\dfrac{\triangle t}{2}\mathrm{\mathbf{u}}^{n+1/2}, t^{n+1/2}\right); \end{align} $

(2.5)

$ \begin{align} x^{n-1}=x_0-2\triangle t\mathrm{\mathbf{u}}_{r}(x^*, t^{n}), \ x^*=x_0-\triangle t\mathrm{\mathbf{u}}(x_0-\triangle t\mathrm{\mathbf{u}}^t, t^n), \end{align} $

(2.6)

其中 $x^*$是根据取中法选取的.

由于流体速度 $\mathrm{\mathbf{u}}$与 $\mathrm{\mathbf{u}}_{r}$在时段 $t^n$与 $t^{n-1}$时只在格点上可知, 而在上游位置 $x^{n}, x^{n-1}$却不可知, 因此, 在空间变量中采用立方体Lagrange插值法, 对 $\mathrm{\mathbf{u}}$与 $\mathrm{\mathbf{u}}_{r}$在其他位置的值进行估计.这一方法具有四阶精度.

在式（2.6）中, $\mathrm{\mathbf{u}}^n$在 $x_0$点估计, 同时在其他位置, $\mathrm{\mathbf{u}}(\cdot, t_n)$可通过方体Lagrange插值方法来求解; 在式(2.5) 中, $\mathrm{\mathbf{u}}^{n+1/2}$在 $x_0$点估计值可利用时间推算近似为 $\frac32\mathrm{\mathbf{u}}^{n}-\frac12\mathrm{\mathbf{u}}^{n-1}$, 同时在其他位置, $\mathrm{\mathbf{u}}(\cdot, t_{n+1/2})$可通过同样的时间推算法与空间插值来求解.一旦 $x^{n}, x^{n-1}$可知, 则可用于求解式(2.2) 中的 $\widetilde{\mathrm{\mathbf{u}}}_{r}^{n}$与 $\widetilde{\mathrm{\mathbf{u}}}_{r}^{n-1}$.同理, 可求得力 $\mathrm{\mathbf{F}}_{b}^{n+1}$.

先采用二阶方案来求解式(2.2) 与式(1.6) 中的第二个散度方程.在式(2.2) 中, 将中间速度 $\mathrm{\mathbf{u}}_{r}^{*}$取代 $\mathrm{\mathbf{u}}_{r}^{n+1}$, 并等价变换为

$ \begin{equation} \left(\dfrac{3}{2\triangle t}-\mu\nabla^2\right)\mathrm{\mathbf{u}}_{r}^{*} =\mathrm{\mathbf{F}}_b^{n+1}-\dfrac{1}{2\triangle t}(-4\widetilde{\mathrm{\mathbf{u}}}_{r}^{n} +\widetilde{\mathrm{\mathbf{u}}}_{r}^{n-1})-\nabla p_{r}^n. \end{equation} $

(2.7)

对式（2.7）两边同时进行快速Fourier变换, 并利用标准的二阶有限差分算子, 则可求出 $\mathrm{\mathbf{u}}_{r}^{*}$.类似于文[14], 将 $\mathrm{\mathbf{u}}_{r}^{*}$近似地映射到自由发散的向量场

$ \begin{equation} \mathrm{\mathbf{u}}_{r}^{n+1}=\mathbb{P}\mathrm{\mathbf{u}}_{r}^{*}, \ \mathbb{P}f=f-\nabla(\nabla^2)^{-1}\nabla\cdot f \end{equation} $

(2.8)

的子空间, 则可定义某个 $\phi$, 使得

$ \mathrm{\mathbf{u}}_{r}^{n+1}=\mathrm{\mathbf{u}}_{r}^{*}-(\triangle t)\nabla\phi. $

由此可得Poisson方程

$ \begin{equation} (\triangle t)\nabla^2\phi=\nabla\cdot\mathrm{\mathbf{u}}_{r}^{*}, \end{equation} $

(2.9)

对式(2.9) 进行Fourier变换并利用差分算子, 则可得出 $\phi$, 并进一步求解 $\mathrm{\mathbf{u}}_{r}^{n+1}$.

在格点空间, 压力的更新迭代方程满足

$ \begin{equation} \nabla p_{r}^{n+1}=\nabla p_{r}^{n}+\frac32\nabla\phi-\mu\triangle t\nabla^3\phi. \end{equation} $

(2.10)

代入式（2.2）, 则可求出方程

$ \begin{equation} \dfrac{3\mathrm{\mathbf{u}}_{r}^{n+1} -4\widetilde{\mathrm{\mathbf{u}}}_{r}^{n} +\widetilde{\mathrm{\mathbf{u}}}_{r}^{n-1}}{2\triangle t}+\nabla p_{r}^{n+1} =\mu\nabla^2\mathrm{\mathbf{u}}_{r}^{n+1}+\mathrm{\mathbf{F}}_b^{n+1} \end{equation} $

(2.11)

的解.值得注意的是, 由于 $\nabla^2\not=\nabla\cdot\nabla$, 故映射 $P$不是带通常的二阶差分算子满射.

当边界力为刚性力时, 对于计算边界, 除了简化跳转条件的处理方法外, 速度分解法也是个有效途径.数值计算中, 如果利用显式方案, 则时间步长必须足够小才能保持迭代的稳定性, 这将产生非常多的非物理耗散, 花费时间较长.因此, 文[16, 22, 23]中采用隐式或半隐式方案来克服这一困难.隐式方案的优点在于, 我们可放松被刚性强加的时间步长限制.但最大的缺点是, 隐式方案会导致边界上的力产生非局部、非线性问题, 其所涉及的边界点相互作用问题, 将使得方案更为抽象.

本节给出一个计算实例, 以对椭圆的放松或振荡运动进行仿真.这个例子在文[11, 16, 22]等处都曾被用于测试相关文献的数值方法.

取初始边界为一个椭圆, 焦点在横坐标上, 且长半轴与短半轴分别为 $a = 0.7, b = 0.5$, 未拉伸的边界取成半径为 $r_0=0.5$的圆, 张力系数设为0.1.仿真流体分为两类, 一类是粘性较强的, 其扩散系数 $\mu=0.1$; 另一类粘性较弱, 其扩散系数为 $\mu=0.01$; 它们的计算区域都是 $[-1.2, 1.2]^2$, 且流体初始时刻为松弛的, 即 $t = 0$时, $\mathbf{u} = 0$且 $p = 0$.

经过足够长的仿真时间后, 接口界面收敛到一个半径为 $r_e=\sqrt{ab}\approx 0.6124$的圆, 其面积与原椭圆面积相等(这是因为封闭流体的不可压缩性造成的), 而比未拉伸的边界圆要大.在稳定状态下, 流体的速度处处都等于0, 而 $p$在边界内外都为常数值, 且跳跃 $[p] < 0$, 这是因为边界被初始化为一个拉伸状态且流体的速度趋于零.

利用速度分解的方法, 可将公式(1.1)-(1.3) 整合到非量纲时间 $t=10$.对于粘性较强与粘性较弱这两种情况, 在接近终止时, 边界都趋于其稳定状态.所不同的是, 前者的收敛过程一直比较平稳, 而后者早期存在振荡.在该例中, 由于边界力量不是特别刚性, 故分级时间步长的方法不是很实用.有鉴于此, 设 $N_m=1, \triangle t_m=\triangle t$.

图 1显示了Stokes速度 $u_s, v_s$、正则速度 $u_r, v_r$以及整个流体的速度 $u, v$在 $t=1.2, x=0.3$的值; 上图针对于粘性较强的流体类型, 即 $\mu=0.1$; 下图针对于粘性较弱的流体类型, 即 $\mu=0.01$.在早期与后期可观察到, 在这两种类型中, $u_s$与 $u$的法向导数的跳跃间断点可被本文的方法精确捕捉, 而 $u_r$的法向导数在穿过 $\Gamma$时是连续的.当 $\mu=0.1$时, Stokes速度 $u_s$在数量上远远超过正则速度 $u_r$; 而当 $\mu=0.01$时, $u_s$与 $u_r$在数量上都远远超过 $u$, 但方向相反.后期的这一效应与早期的相比, 稍微弱一点.

由于初始速度为0, 且在零时刻 $u_s=-u_r$, 故当流体粘度变小时, 可预计Stokes速度 $\mathbf{u}_s$与Navier-Stokes速度 $\mathbf{u}$的相关性会减弱, 即Reynolds数会增大.