相关申请资料
本申请与序列号为10/390,239的申请相关，该申请于2003年3 月14日提交，并且题目为“Coupled Quadrilateral Grid Level Set Scheme for Piezoelectric Ink-Jet Simulation(用于压电喷墨模 拟的耦合四边形网格水平集方案)”。在此并入该相关申请的公开内 容来作为参考，并且在下面的说明中将其称作“相关申请”。
上述的相关申请公开了在四边形网格上的耦合水平集投影方法。 在那个发明中，速度分量和水平集位于单元中心，同时压力在网格点 处。Navier-Stokes方程在不考虑不可压缩性的条件下首先在每个时 间步长(time step)上被求解。然后，获得的速度被“投影”到无发 散(divergence-free)的场的空间中。在四边形网格上的Godunov 逆风(upwind)方案被用于评估Navier-Stokes方程中的对流项和水 平集对流方程。在时间和空间上进行泰勒展开以获得边缘速度和水平 集。由于外推法可从垂直单元边缘的左侧和右侧被实施，以及可从水 平边缘的上侧和下侧被实施，因此在每个单元边缘有两个外推值(参 见相关申请的方程式(50)-(52))。然后，Godunov逆风过程被使 用以决定采用哪个外推值(参见相关申请的方程式(53)和(54))。 在Godunov逆风过程中，局部Riemann问题实际上被解决。对于 Newtonian流体，局部Riemann问题简化为一维Burger方程式，它具 有如由相关申请的方程式(53)和(54)给出的非常简单的解。因此， 基本上，局部速度方向决定了哪个外推值用于Newtonian流体。然而， 如果流体不是Newtonian流体，则局部Riemann问题不可能具有简单 的解，并且因此Godunov逆风过程将非常难于编程。用于平流项 (advection term)的这种逆风方案的优点是较低的网格粘性和较高 的数值分辨率。它的缺点是增加的代码复杂性，该代码复杂性是在时 间和空间上进行泰勒外推和Godunov逆风过程所需要的。

考虑到前述内容，本发明的目的是提供一种改进的模型和伴随的 算法，以模拟和分析来自压电打印头的喷墨。
另一目的是包括在改进的模型和算法中用于耦合水平集投影方法 的中心差分方案，以及可选择地进一步包括可用于中心差分方案的多 重网格兼容离散投影算子。
在本发明的中心差分方案中，不考虑不可压缩的条件，在每个时 间步长中再次首先求解Navier-Stokes方程。然后，获得的速度被“投 影”到无发散的场的空间中。但是，与相关申请的逆风方案对比，这 里的中心差分方案使用交错的单元平均，并且局部速度的方向是不重 要的。在时间上的泰勒外推仍被用来计算以时间为中心的(time- centered)速度和水平集预测算子。但由于它仅仅是在时间上的外推， 因此代码执行容易得多。该方案的一个特征是，如果在时间步长t＝tn 时速度和水平集位于单元中心处，则它们将在下一个时间步长t＝tn+1 时位于网格点处。在时间步长t＝tn+2时，速度和水平集将再次位于单元 中心处。在整个时间步长系列上，速度和水平集交替布置，这就是为 什么该方法被描述为使用“交错的”网格的原因。
为了改善该中心差分方案的性能，两个离散投影算子被构造以将 速度投影到无发散的空间中。从投影算子得到的线性系统可通过多重 网格方法被快速求解。投影算子包括有限差分投影和有限元投影，它 们被用在不同的条件下。当速度位于网格点处而压力位于单元中心处 时，使用有限差分投影算子。当速度位于单元中心处而压力位于网格 点处时，使用有限元投影算子。当采用有限差分投影算子时，单元密 度在时间上滞后半个时间步长。这个单元密度时间延迟被称为“延迟” 密度。
如本发明所打算的，该中心差分方案的优点包括(i)容易执行； (ii)较低的计算机存储要求；以及(iii)更容易推广应用至non- Newtonian流体流动。
依据本发明的一个方面，提供了一种用于模拟和分析来自通道 (channel)的流体喷射的方法。优选地，该通道包括作为压电喷墨头 的一部分的喷墨喷嘴。在流过该通道的第一流体(例如墨水)和第二 流体(例如空气)之间存在一个边界。该方法包括下列步骤：(a)用 公式表示(formulate)在均匀矩形网格上基于中心差分的离散化；(b) 使用基于中心差分的离散化来在均匀矩形网格上执行有限差分分析， 以求解控制通过该通道的至少第一流体的流动的方程；和(c)基于执 行的有限差分分析，模拟第一流体通过该通道的流动和从该通道的喷 射。
优选地，在执行步骤(b)中，水平集方法被用于收集通道中的边 界的特征。
优选地，均匀矩形网格包括多个单元，并且用公式表示基于中心 差分的离散化，以使对于每一单元，存在用于确定第一流体速度的速 度向量和用于收集通道中边界的特征的水平集值。对于每一单元，对 应的速度向量和水平集值在一个时间点位于该单元的近似中心，且在 下一个时间点位于该单元的网格点。
优选地，步骤(a)包括构造多重网格兼容的投影算子，该算子包 括有限差分投影算子和有限元投影算子。
优选地，在步骤(b)中，使用基于中心差分的离散化在均匀矩形 网格上执行有限差分分析以求解方程，该方程控制在每个时间点在各 个单元中通过通道的至少第一流体的流动，在步骤(b)后施加的该有 限差分投影算子在速度向量和水平集值位于网格点的情况下被执行， 并且在步骤(b)后施加的该有限元投影算子在速度向量和水平集值位 于近似单元中心的情况下被执行。
在另一方面，本发明包括一种用于模拟和分析来自通道的流体喷 射的装置，在流过该通道的第一流体和第二流体之间有一个边界。该 装置包括一个或多个被配置成执行上述处理的部件或模块。这些处理 的一些或全部可被以软件、硬件或其组合的形式实现的指令程序来方 便地规定。所述部件之一优选地包括用于能够可视化观察模拟的显示 器。
依据本发明的另一方面，本发明的多个方面可以指令程序的形式 被实现，该指令程序可以是软件的形式，该软件可以被存储或传送至 计算机或其他受处理器控制的设备以用于执行。可替换地，指令程序 可直接以硬件(例如专用集成电路(ASIC)、数字信号处理电路等) 来实现，或者这种指令可被实现为软件和硬件的组合。
通过参考结合附图的下面的说明和权利要求书，更全面地理解本 发明的其它目的和成果将变得清楚和明白。

在附图中，相同的附图标记指的是相同的部件：
图1是显示在t＝tn处速度和压力的位置的网格的示意图；
图2是在本发明的实施例中使用的交错网格的示意图；
图3是说明依据本发明的实施例的算法的流程图；
图4是说明在t＝tn+1处非无发散(non-divergence-free)速度校 正算子和压力的位置的示意图；
图5是说明在t＝tn+2处非无发散速度校正算子和压力的位置的示意 图；
图6是对流体泡落在较轻的周围流体中的模拟；
图7是在喷墨模拟中被施加在喷嘴流入上的动态压力；
图8是2D小滴喷射的模拟；以及
图9是说明可用于执行本发明的多个方面的示例性系统的框图。

I.引言
本发明涉及在均匀矩形网格上用于尤其是喷墨的二维两相流的模 拟的中心差分水平集投影方法的开发。Navier-Stokes方程在不考虑 不可压缩性的条件下首先在每个时间步长上被求解。然后，获得的速 度被“投影”到无发散的场的空间中。与相关申请的逆风方案对比， 在此的中心差分方案使用交错单元平均，并且局部速度的方向是不重 要的。为了求解Navier-Stokes方程，首先计算速度和水平集斜率。 在时间上的泰勒外推被用于计算以时间为中心的速度和水平集预测算 子。这些预测算子然后被用于计算速度校正算子，该速度校正算子不 是无发散的。
该方案的一个特征是，如果在时间步长t＝tn处速度和水平集位于 单元中心，则在下一个时间步长t＝tn+1处它们将位于网格点上。在时间 步长t＝tn+2处，速度和水平集将再次位于单元中心。速度和水平集在整 个系列时间步长上交替布置是该方法被描述为“交错的”网格的原因。
为了改善该中心差分方案的性能，两个离散投影算子被构造成将 速度投影到无发散空间中。从投影算子得到的线性系统可通过多重网 格方法快速求解。投影算子包括有限差分投影和有限元投影，它们被 用于不同的情况中。当速度位于网格点处且压力在单元中心处时，使 用有限差分投影算子。当速度位于单元中心处且压力在网格点处时， 使用有限元投影算子。当采用有限差分投影算子时，单元密度在时间 上滞后半个时间步长。这个单元密度时间延迟被称为“延迟”密度。
本发明的中心差分算法比使用Godunov逆风方案的模拟算法更 快，并具有较低的计算机开销要求，因为：(1)只使用在时间上的泰 勒外推，而不是在时间和空间上的泰勒外推；(2)没有要求解的局部 Riemann问题；(3)在每个时间步长上只求解一个投影方程(在相关 申请中在每个时间步长上求解两个方程，即MAC投影和FEM投影)。
II.运动方程
下面的说明解释了如何由Kupferman和Tadmor将用于单相流的中 心差分方案扩展至两相流模拟。求解喷墨模拟的控制方程在这部分中 被确定并被表述在附件中(与在这里参考的其它用数字标记的方程一 起)。由于我更喜欢使用用于平流项的守恒形式，因此Navier-Stokes 方程和水平集对流方程与在相关申请中给出的那些方程略有不同。用 于两相流的中心差分方案和说明算法实施例的流程图在部分III中进行 了解释。与快速多重网格方法相兼容的两个离散投影算子在部分IV中 进行了解释。在部分V中给出了模拟的例子。
A.控制方程
用于两相流的控制方程包括连续性方程(1)和Navier-Stokes 方程(2)。
              ·u＝0，                  (1)
$\rho \left(\phi \right)\frac{\mathrm{Du}}{\mathrm{Dt}}=-▿p+▿·\left(2\mu \left(\phi \right)D\right)-\mathrm{\sigma \kappa }\left(\phi \right)\delta \left(\phi \right)▿\phi .---\left(2\right)$
在这些方程中，形变张量的速率和流体速度分别由方程(3)确定。
$D=\frac{1}{2}[▿u+{(▿u)}^T],---\left(3\right)$
u＝ue1+υe2，
$\frac{D}{\mathrm{Dt}}=\frac{\partial }{\partial t}+(u·▿)$ 是拉格朗日时间导数，ρ是相对密度，p是压力，μ是 相对动态粘性，σ是表面张力系数，к是曲率，δ是迪拉克δ函数，以及 φ是水平集。
水平集公式用于定义在两个流体譬如墨水和空气之间的界面，因 此相对密度、相对动态粘性和曲率全部被定义在如方程(4)所述的水 平集φ的项中。

这里，水平集函数φ被初始化为到界面的带符号的距离，即水平集 值在流体侧是到界面的最短距离，而在空气侧是最短距离的负值。与 该界面正交、从流体2被吸入流体1的单位以及界面的曲率可以如方 程(5)中φ的项来表示。
$n=\frac{▿\phi }{|▿\phi |}{|}_{\phi =0},$
$\kappa =▿·\left(\frac{▿\phi }{|▿\phi |}\right){|}_{\phi =0}.---\left(5\right)$
为了使控制方程无量纲，选择如在方程(6)中所述的下面定义。
  x＝Lx′，y＝Ly′，u＝Uu′， $t=\frac{L}{U}{t}^{\prime },---\left(6\right)$
  p＝ρ1U2p′，ρ＝ρ1ρ′，μ＝μ1μ′，
最初量是无量纲的，并且L、U、ρ1、μ1分别表示特征长度、特征 速度、流体1的密度和流体1的动态粘性。特征速度和特征长度可被 任意选择，因为它们不影响模拟的结果。
将方程(6)的关系代入方程(1)和(2)，并降低该最初值，得 到方程(7)和(8)，其中相对密度比率ρ(φ)、相对粘性比率μ(φ)、 雷诺数Re和韦伯数均由方程(9)定义。
         ·u＝0，                      (7)
$\frac{\partial u}{\partial t}+(u·▿)u=-\frac{1}{\rho \left(\phi \right)}▿p+\frac{1}{\rho \left(\phi \right)\mathrm{Re}}▿·\left(2\mu \left(\phi \right)D\right)-\frac{1}{\rho \left(\phi \right)\mathrm{We}}\kappa \left(\phi \right)\delta \left(\phi \right)▿\phi ,---\left(8\right)$
$\rho \left(\phi \right)=\left\{\begin{array}{ccc}1& \mathrm{if}& \phi \ge 0\\ {\rho }_1/{\rho }_2& \mathrm{if}& \phi <0\end{array}\right.,$
$\mu \left(\phi \right)=\begin{array}{}\end{array}\left\{\begin{array}{ccc}1& \mathrm{if}& \phi \ge 0\\ {\mu }_2/{\mu }_1& \mathrm{if}& \phi <0\end{array}\right.,---\left(9\right)$
$\mathrm{Re}=\frac{{\rho }_1\mathrm{UL}}{{\mu }_1},$
$\mathrm{We}=\frac{{\rho }_1{U}^{2}L}{\sigma }.$
由于界面随流体移动，所以水平集的演变由方程(10)控制。
$\frac{\partial \phi }{\partial t}+u·▿\phi =0.---\left(10\right)$
由于方程(5)、(7)和(8)以向量记法的项来表示，所以它们 在笛卡儿坐标系和轴对称坐标系中采用相同的形式。
为了说明本发明，我喜欢使用对流项的守恒形式。因此，方程(8) 和(10)被分别重写为方程(11)和(12)。
$\frac{\partial u}{\partial t}+▿·\left(\mathrm{uu}\right)=-\frac{1}{\rho \left(\phi \right)}▿p+\frac{1}{\rho \left(\phi \right)\mathrm{Re}}▿·\left(2\mu \left(\phi \right)D\right)-\frac{1}{\rho \left(\phi \right)\mathrm{We}}\kappa \left(\phi \right)\delta \left(\phi \right)▿\phi ,---\left(11\right)$
$\frac{\partial \phi }{\partial t}+▿·\left(\mathrm{u\phi }\right)=0.---\left(12\right)$
方程(8)、(10)和(11)、(12)的等价性可由扩展对流项来 确认。例如，在方程(12)中的扩展在方程(13)中被说明。只要流 体是不可压缩的，即只要·u＝0，则保持最终的等价性。
$▿·\left(\mathrm{u\phi }\right)=\frac{\partial \left(\mathrm{u\phi }\right)}{\partial x}+\frac{\partial \left(\mathrm{\upsilon \phi }\right)}{\partial y}---\left(13\right)$
$=\frac{\partial \left(\mathrm{u\phi }\right)}{\partial x}+\frac{\partial \left(\mathrm{\upsilon \phi }\right)}{\partial y}$
$=\phi (u_{,x}+{\upsilon }_{,y})+u{\phi }_{,x}+{\mathrm{\upsilon \phi }}_{,y}$
$=\phi ▿·u+(u·▿)\phi$
$=(u·▿)u.$
B.边界条件
在固体壁上，假定速度的法线和切线分量变为零，尽管在三相点 处这个假定必须被改变。在流入和流出处，本发明的公式表示允许我 们规定速度(14)或压力(15)的边界条件。在方程(15)中，符号n 表示正交于流入或流出边界的单位。
            u＝uBC                    (14)
$p=$$p^{\mathrm{BC}},\frac{\partial u}{\partial n}=0,---\left(15\right)$
对于喷墨模拟，与时间相关的流入条件由模拟电荷驱动机构的等 效电路模型来提供，该电荷驱动机构迫使墨水从槽中进入喷嘴。
C.接触模型
在三相点处，其中空气和墨水在固体壁处相遇，如在2002年3月 22日提交的且题目为“A Slipping Contact Line Model and Mass-Conservative Level Set Implementation for Ink-Jet Simulation(执行用于喷墨模拟的滑动接触线模型和质量守恒水平 集)”的序列号为10/105,138的申请中给出的滑动接触线模型被使 用。该申请的公开内容在此被结合以供参考。
III.中心差分方案
现在用公式表示在均匀矩形网格上的数值算法。在下面，上标n (或n+1)表示时间步长，即方程(16)等。
      un＝u(t＝tn)                      (16)
假定量un、pn、φn，该算法的目的是获得满足不可压缩条件的un+1、 pn+1、φn+1。这里说明的显性模式算法在时间上是一阶精确的而在空间 上是二阶精确的。
A.浸润(smearing)界面
因为由迪拉克δ函数和由穿越界面的密度和粘性的急剧变化引起 的数值困难，所以Heaviside和迪拉克δ函数由平滑的函数取代，即用 来浸润一点界面。Heaviside函数由方程(17)重新定义。因此，平 滑的δ函数如在方程(18)中所述。
$H\left(\phi \right)=\left\{\begin{array}{ccc}0& \mathrm{if}& \phi <-\in \\ \frac{1}{2}[1+\frac{\phi }{\in }+\frac{1}{\pi }\sin (\mathrm{\pi \phi }/\in )]& \mathrm{if}& \left|\phi \right|\le \in \\ 1& \mathrm{if}& \phi >\in \end{array}\right.---\left(17\right)$
$\delta \left(\phi \right)=\frac{\mathrm{dH}\left(\phi \right)}{\mathrm{d\phi }}.---\left(18\right)$
通常选择参数ε与如方程(19)中所述的单元的平均大小成比例， 其中Δx是矩形单元的平均大小。我通常选择α在1.7和2.5之间。
          ε＝αΔx，                    (19)
因此，密度和粘性变得如方程(20)所述。
$\rho \left(\phi \right)=\frac{{\rho }_2}{{\rho }_1}+H\left(\phi \right)(1-\frac{{\rho }_2}{{\rho }_1}),---\left(20\right)$
$\mu \left(\phi \right)=\frac{{\mu }_2}{{\mu }_1}+H\left(\phi \right)(1-\frac{{\mu }_2}{{\mu }_1}).$
B.中心方案
网格由大小为Δx×Δy的均匀矩形单元组成。这些单元以 (xi＝(i-1/2)Δx，yj＝(j-1/2)Δy)为中心。在时间步长tn，假定离散速度 un＝(ui，j，vi，j)和水平集φi，j均位于单元中心，而压力pi，j n在网格点处，如图 1所示。目的是获得在t＝tn+1处的速度、压力和水平集。
B.1插值
第一级是从位于单元中心处的离散速度和水平集值线性地重构该 场，如方程(21)所述。
$u^n(x,y)=u_{i,j}^n+\frac{u_{i,j}^{\prime }}{\mathrm{\Delta x}}(x-x_i)+\frac{u_{i,j}^{\prime }}{\mathrm{\Delta y}}(y-y_j),$
${\phi }^n(x,y)={\phi }_{i,j}^n+\frac{{\phi }_{i,j}^{\prime }}{\mathrm{\Delta x}}(x-x_i)+\frac{{\phi }_{i,j}^{\prime }}{\mathrm{\Delta y}}(y-y_j).---\left(21\right)$
在这些方程中，
和 和 分别是在x和y方向上的离散 速度(水平集)斜率。有几种计算这些斜率的方法。例如，一种方法 可使用如在上述相关申请中的方程(60)-(62)中讨论的二阶单调受 限斜率。可替换地，如在此优选的，使用在方程(22)中所述的简单 的中心差分来计算这些斜率。
$\frac{u_{i,j}^{\prime }}{\mathrm{\Delta x}}=\frac{u_{i+1,j}^n-u_{i-1,j}^n}{2\mathrm{\Delta x}},\frac{u_{i,j}^{\prime }}{\mathrm{\Delta y}}=\frac{u_{i,j+1}^n-u_{i,j-1}^n}{2\mathrm{\Delta y}},$
$\frac{{\phi }_{i,j}^{\prime }}{\mathrm{\Delta x}}=\frac{{\phi }_{i+1,j}^n-{\phi }_{i-1,j}^n}{2\mathrm{\Delta x}},\frac{{\phi }_{i,j}^{\prime }}{\mathrm{\Delta y}}=\frac{{\phi }_{i,j+1}^n-{\phi }_{i,j-1}^n}{2\mathrm{\Delta y}}.---\left(22\right)$
B.2预测算子
第二级是获得在t＝tn+1/2处速度和水平集的预测算子。这可如按照 方程(23)的由在时间上的泰勒展开式容易地实现。
$u_{i,j}^{n+1/2}=u_{i,j}^n+\frac{\Delta t}{2}\frac{\partial u}{\partial t}{|}_{t=t^n},$
${\phi }_{i,j}^{n+1/2}={\phi }_{i,j}^n+\frac{\mathrm{\Delta t}}{2}\frac{\partial \phi }{\partial t}{|}_{t=t^n}.---\left(23\right)$
通过代入方程(11)和(12)得到方程(24)来估计在(23)右 侧的时间导数。在单元中心处获得水平集预测算子φi，j n+1/2之后，使用方 程(20)来计算以时间为中心的密度(ρ(φn+1/2))和粘性(μ(φn+1/2))。注 意，这些密度和粘性在t＝tn+1/2处仍然是以单元为中心的。
$u_{i,j}^{n+1/2}=u_{i,j}^n+\frac{\mathrm{\Delta t}}{2}\{-\frac{{\partial p}_{i,j}^n}{\partial x}-(2\frac{u_{i,j}^{\prime }}{\mathrm{\Delta x}}{u}_{i,j}^n+\frac{u_{i,j}^{\prime }}{\mathrm{\Delta y}}{\upsilon }_{i,j}^n+\frac{{\upsilon }_{i,j}^{\prime }}{\mathrm{\Delta y}}{u}_{i,j}^n)$
$+\frac{1}{\rho \left({\phi }^n\right)\mathrm{Re}}[{\left(2\mu \left(\phi \right)u_{,x}\right)}_{,x}+{\left(\mu \left(\phi \right)(u_{,y}+{\upsilon }_{,x})\right)}_{,y}{]}_{t=t^n}$
$-\frac{1}{\rho \left({\phi }^n\right)\mathrm{We}}\kappa \left({\phi }^n\right)\delta \left({\phi }^n\right){\phi }_{,x}^n\},$
${\upsilon }_{i,j}^{n+1/2}={\upsilon }_{i,j}^n+\frac{\mathrm{\Delta t}}{2}\{-\frac{{\partial p}_{i,j}^n}{\partial y}-(\frac{u_{i,j}^{\prime }}{\mathrm{\Delta x}}{\upsilon }_{i,j}^n+\frac{{\upsilon }_{i,j}^{\prime }}{\mathrm{\Delta x}}{u}_{i,j}^n+2\frac{{\upsilon }_{i,j}^{\prime }}{\mathrm{\Delta y}}{\upsilon }_{i,j}^n)---\left(24\right)$
$+\frac{1}{\rho \left({\phi }^n\right)\mathrm{Re}}[{\left(\mu \left(\phi \right)(u_{,y}+{\upsilon }_{,x})\right)}_{,x}+{\left(2\mu \left(\phi \right){\upsilon }_{,y}\right)}_{,y}{]}_{t=t^n}$
$-\frac{1}{\rho \left({\phi }^n\right)\mathrm{We}}\kappa \left({\phi }^n\right)\delta \left({\phi }^n\right){\phi }_{,y}^n\},$
${\phi }_{i,j}^{n+1/2}={\phi }_{i,j}^n-\frac{\mathrm{\Delta t}}{2}\{\frac{1}{\mathrm{\Delta x}}(u_{i,j}^{\prime }{\phi }_{i,j}^n+u_{i,j}^n{\phi }_{i,j}^{\prime })+\frac{1}{\mathrm{\Delta y}}({\upsilon }_{i,j}^{\prime }{\phi }_{i,j}^n+{\upsilon }_{i,j}^n{\phi }_{i,j}^{\prime })\}.$
B.3校正算子
第三级是将在tn处的分段线性近似(21)引中至tn+1。新的速度场 和水平集首先通过如方程(25)所述的它们的交错单元平均来实现。 该方案的中心差分特性被牢固地链接至在(25)中的交错的平均，其 应当被集成在如图2所示的控制箱 中。以水平集对流方 程为例，将方程(26)代入方程(25)的第一方程中以得到方程(27)。 对于在方程(27)右侧的第一项，插入线性重构的水平集(即在(21) 中的第二方程)。对于第二和第三项，使用在空间和时间上的中间点 求积分。以第二项为例，该代入被显示在方程(28)中，其中以时间 为中心的量是那些在前述级获得的速度和水平集预测算子，以及导数 算子Dx +和平均算子μy +被定义为如方程(29)中所示。在下面，我同样 使用定义的导数和平均算子比如Dy +、Dx -、μx +等，这些算子类似于(27) 被定义。
${\bar{\phi }}_{i+1/2,j+1/2}^{n+1}=\frac{1}{\mathrm{\Delta x\Delta y}}{\int }_{C_{i+1/2,j+1/2}}\phi (x,y,t^{n+1})\mathrm{dxdy},$
${\bar{u}}_{i+1/2,j+1/2}^{n+1}=\frac{1}{\mathrm{\Delta x\Delta y}}{\int }_{C_{i+1/2,j+1/2}}u(x,y,t^{n+1})\mathrm{dxdy}.---\left(25\right)$
$\phi ={\phi }^n-\int ▿·\left(\mathrm{u\phi }\right)\mathrm{dt}---\left(26\right)$
${\bar{\phi }}_{i-1/2,j+1/2}^{n+1}=\frac{1}{\mathrm{\Delta x\Delta y}}{\int }_{C_{i+1/2,j+1/2}}\phi (x,y,t^n)\mathrm{dxdy}---\left(27\right)$
$-\frac{1}{\mathrm{\Delta y}}{D}_x^{+}{\int }_{\tau =t^n}^{t^{n+1}}{\int }_{{y\in J}_{j+1/2}}\left(\mathrm{u\phi }\right)\mathrm{dyd\tau }-\frac{1}{\mathrm{\Delta x}}{D}_y^{+}{\int }_{\tau =t^n}^{t^{n+1}}{\int }_{{x\in I}_{i+1/2}}\left(\mathrm{\upsilon \phi }\right)\mathrm{dxd\tau }.$
$\frac{1}{\mathrm{\Delta y}}{D}_x^{+}{\int }_{\tau =t^n}^{t^{n+1}}{\int }_{{y\in J}_{j+1/2}}\left(\mathrm{u\phi }\right)\mathrm{dyd\tau }---\left(28\right)$
$=\mathrm{\Delta t}\left[D_x^{+}{\mu }_y^{+}\left(u_{i,j}^{n+1/2}{\phi }_{i,j}^{n+1/2}\right)\right],$
$D_x^{+}{u}_{i,j}=\frac{u_{i+1,j}-u_{i-1,j}}{2\mathrm{\Delta x}},{\mu }_y^{+}{u}_{i,j}\frac{u_{i,j+1}+u_{i,j-1}}{2}.---\left(29\right)$
在插入线性重构场和预测算子之后，在方程(27)中的交错水平 集平均变得如在方程(30)中所示。同样地，速度场如在方程(31) 和(32)中所述。
${\bar{\phi }}_{i+1/2,j+1/2}^{n+1}={\mu }_x^{+}{\mu }_y^{+}{\phi }_{i,j}^n-\frac{\mathrm{\Delta x}}{8}{D}_x^{+}{\mu }_y^{+}{\phi }_{i,j}^{\prime }-\frac{\mathrm{\Delta y}}{8}{D}_y^{+}{\mu }_x^{+}{\phi }_{i,j}^{\prime }---\left(30\right)$
$-\mathrm{\Delta t}[D_x^{+}{\mu }_y^{+}\left(u_{i,j}^{n+1/2}{\phi }_{i,j}^{n+1/2}\right)+D_y^{+}{\mu }_x^{+}\left({\upsilon }_{i,j}^{n+1/2}{\phi }_{i,j}^{n+1/2}\right)].$
${\bar{u}}_{i+1/2,j+1/2}^{n+1}=\frac{1}{\mathrm{\Delta x\Delta y}}{\int }_{C_{i+1/2,j+1/2}}u(x,y,t^n)\mathrm{dxdy}$
$-\frac{1}{\mathrm{\Delta y}}{D}_x^{+}{\int }_{\tau =t^n}^{t^{n+1}}{\int }_{{y\in J}_{j+1/2}}{u}^2\mathrm{dyd\tau }-\frac{1}{\mathrm{\Delta x}}{D}_y^{+}{\int }_{\tau =t^n}^{t^{n+1}}{\int }_{{x\in I}_{i+1/2}}\mathrm{\upsilon udxd\tau }$
$+\frac{1}{\mathrm{\Delta x\Delta y}}{\int }_{\tau =t^n}^{t^{n+1}}{\int }_{C_{i+1/2,j+1/2}}\frac{1}{\rho \left(\phi \right)\mathrm{Re}}\{{\left(2\mu \left(\phi \right)u_{,x}\right)}_{,x}+[\mu \left(\phi \right)(u_{,y}+{\upsilon }_{,x}){]}_{,y}\}\mathrm{dxdyd\tau }$
$-\frac{1}{\mathrm{\Delta x\Delta y}}{\int }_{\tau =t^n}^{t^{n+1}}{\int }_{C_{i+1/2,j+1/2}}\frac{1}{\rho \left(\phi \right)\mathrm{We}}\kappa \left(\phi \right)\delta \left(\phi \right){\phi }_{,x}\mathrm{dxdyd\tau }$
$={\mu }_x^{+}{\mu }_y^{+}{u}_{i,j}^n-\frac{\mathrm{\Delta x}}{8}{D}_x^{+}{\mu }_y^{+}{u}_{i,j}^{\prime }-\frac{\mathrm{\Delta y}}{8}{D}_y^{+}{\mu }_x^{+}{u}_{i,j}^{\prime }$
$-\mathrm{\Delta t}[D_x^{+}{\mu }_y^{+}\left(u_{i,j}^{n+1/2}{u}_{i,j}^{n+1/2}\right)+D_y^{+}{\mu }_x^{+}\left({\upsilon }_{i,j}^{n+1/2}{u}_{i,j}^{n+1/2}\right)]$
$+\mathrm{\Delta t}\frac{{\mu }_x^{+}{\mu }_y^{+}}{\rho \left({\phi }^{n+1/2}\right)\mathrm{Re}}\{{\left(2\mu \left(\phi \right)u_{,x}\right)}_{,x}+[\mu \left(\phi \right)(u_{,y}+{\upsilon }_{,x}){]}_{,y}{\}}_{i,j}^{n+1/2}$
$-\mathrm{\Delta t}{\mu }_x^{+}{\mu }_y^{+}\{\frac{1}{\rho \left(\phi \right)\mathrm{We}}\kappa \left(\phi \right)\delta \left(\phi \right){\phi }_{,x}{\}}_{i,j}^{n+1/2},---\left(31\right)$
和
${\bar{\upsilon }}_{i+1/2,j+1/2}^{n+1}=\frac{1}{\mathrm{\Delta x\Delta y}}{\int }_{C_{i+1/2,j+1/2}}\upsilon (x,y,t^n)\mathrm{dxdy}$
$-\frac{1}{\mathrm{\Delta y}}{D}_x^{+}{\int }_{\tau =t^n}^{t^{n+1}}{\int }_{{v\in J}_{j+1/2}}\mathrm{u\upsilon dyd\tau }-\frac{1}{\mathrm{\Delta x}}{D}_y^{+}{\int }_{\tau {=t}^n}^{t^{n+1}}{\int }_{{x\in I}_{i+1/2}}{\upsilon }^2\mathrm{dxd\tau }$
$+\frac{1}{\mathrm{\Delta x\Delta y}}{\int }_{\tau =t^n}^{t^{n+1}}{\int }_{C_{i+1/2,j+1/2}}\frac{1}{\rho \left(\phi \right)\mathrm{Re}}\left\{\right[\mu \left(\phi \right)(u_{,y}+{\upsilon }_{,x}){]}_{,x}+{\left(2\mu \left(\phi \right){\upsilon }_{,y}\right)}_{,y}\}\mathrm{dxdyd\tau }$
$-\frac{1}{\mathrm{\Delta x\Delta y}}{\int }_{\tau =t^n}^{t^{n+1}}{\int }_{C_{i+1/2,j+1/2}}\frac{1}{\rho \left(\phi \right)\mathrm{We}}\kappa \left(\phi \right)\delta \left(\phi \right){\phi }_{,y}\mathrm{dxdyd\tau }$
$=$${\mu }_x^{+}{\mu }_y^{+}{\upsilon }_{i,j}^n-\frac{\mathrm{\Delta x}}{8}{D}_x^{+}{\mu }_y^{+}{\upsilon }_{i,j}^{\prime }-\frac{\mathrm{\Delta y}}{8}{D}_y^{+}{\mu }_x^{+}{\upsilon }_{i,j}^{\prime }$
$-\mathrm{\Delta t}[D_x^{+}{\mu }_y^{+}\left(u_{i,j}^{n+1/2}{\upsilon }_{i,j}^{n+1/2}\right)+D_y^{+}{\mu }_x^{+}\left({\upsilon }_{i,j}^{n+1/2}{\upsilon }_{i,j}^{n+1/2}\right)]$
$+\mathrm{\Delta t}\frac{{\mu }_x^{+}{\mu }_y^{+}}{\rho \left({\phi }^{n+1/2}\right)\mathrm{Re}}\{{\left[\mu \left(\phi \right)(u_{,y}+{\upsilon }_{,x})\right]}_{,x}+{\left(2\mu \left(\phi \right){\upsilon }_{,y}\right)}_{,y}{\}}_{i,j}^{n+1/2}$
$-\mathrm{\Delta t}{\mu }_x^{+}{\mu }_y^{+}\{\frac{1}{\rho \left(\phi \right)\mathrm{We}}\kappa \left(\phi \right)\delta \left(\phi \right){\phi }_{,y}{\}}_{i,j}^{n+1/2}....\left(32\right)$
注意，这些在t＝tn+1时的水平集和速度校正算子位于网格点处。由 于不可压缩的条件没有被加强，所以由方程(31)和(32)给出的速 度场不是无发散的。
B.4用于un+1的投影
方程(33)中所述的在t＝tn+1时速度场un+1是不可压缩的。如果我 们将发散算子应用至(33)的两边并且注意·un+1＝0，那么将获得在方 程(34)中的结果。
$u^{n+1}={\bar{u}}^{n+1}-\frac{1}{{\rho }^{n+1/2}}▿p^{n+1}---\left(33\right)$
$▿·(\frac{\mathrm{\Delta t}}{\rho \left({\phi }^{n+1/2}\right)}▿p^{n+1})=▿·{\bar{u}}^{n+1}.---\left(34\right)$
投影方程(34)是椭圆形的。如果密度比率ρ(φn+1/2)是常数，则它 简化为泊松方程。投影方程的有限元公式表示被显示在方程(35)中，
${\int }_{\Omega }{u}^{*}·▿\mathrm{\psi dx}={\int }_{\Omega }\frac{\mathrm{\Delta t}}{\rho \left({\phi }^{n+1/2}\right)}▿p^{n+1}·▿\mathrm{\psi dx}+{\int }_{{\Gamma }_1}{\mathrm{\psi u}}^{\mathrm{BC}}·\mathrm{ndS},---\left(35\right)$
其中，Γ1表示其中给出流入或流出速度uBC的全部边界。发散理论可以 证明，在Γ1处的隐含的边界条件如方程(36)中所示。
$\frac{\mathrm{\Delta t}}{\rho \left({\phi }^{n+1/2}\right)}\frac{\partial p^{n+1}}{\partial n}=(u^{*}-u^{\mathrm{BC}})·n.---\left(36\right)$
值得注意的是，如果在流入和流出处边界压力都被给定，则在 (35)的右侧的第二项将变为零。投影步骤可由离散化方程(34)或 (35)来进行。
从方程(34)或(35)解出压力pn+1之后，速度场un+1可由方程(33) 获得。
B.5交错特性
依据前述分部分，中心差分方案首先在没有连续性的条件下积分 Navier-Stokes方程。在t＝tn，速度和水平集位于单元中心处，而压 力在网格点处。使用在时间上的泰勒展开式来首先计算以时间为中心 的速度和水平集预测算子。使用以时间为中心的水平集来计算以时间 为中心的密度预测算子。然后，这些位于单元中心的预测算子被用于 计算在网格点处的速度和水平集校正算子。由于我们想要使压力和速 度之间的相对位置保持相同，所以在投影步骤中获得的速度un+1位于网 格点处，而压力pn+1在单元中心处。
从t＝tn+1到t＝tn+2的情形是存在的另一种方式。在t＝tn+1，速度位 于网格点处，而压力在单元中心。中心差分方案首先计算在网格点处 的以时间为中心的速度和水平集预测算子。然后，这些预测算子被用 于在单元中心处获得校正算子。在进行投影后，不可压缩的速度un+2 位于单元中心处，而压力pn+2在网格点。由于在我们的数值方案中离散 速度和压力的位置随每一个时间步长而改变，所以它具有所谓的“交 错特性”。
C.水平集的重新初始化
为了正确地捕获界面并准确地计算表面张力，水平集需要保持为 界面的带符号的距离函数。然而，如果水平集由方程(7)更新，就不 需要这样保持。因此，代之以该模拟被周期性地停止，并且作为带符 号的距离函数的新水平集函数φ被重新产生，即|φ|＝1，而没有改变 原始水平集函数的零水平集。
前面已经认识到需要在水平集计算中这样做，如重新初始化。用 于重新初始化的直接而简单的方法是使用轮廓绘图仪首先找到界面( 零水平集)，然后重新计算从各个单元至界面的带符号距离。另一简 单的重新初始化选择是解决如方程(37)所述的穿越时间(crossing time)问题。
φt′+F|φ|＝0，                      (37)
其中F是给定的法向速度。值得注意的是，t′已被用于方程(37)中以 强调它是伪时间变量，并且该方程完全为了重新初始化的目的而被求 解。在F＝1的情况下，在时间上向前和向后流过界面，并计算时间t′， 在时间t′处水平集函数对于各个单元改变符号。穿越时间(正的和负 的)等于带符号的距离，这两个简单的方法中任何一个均适合供本发 明使用。两者均已被尝试，其中在模拟结果中没有显著差别。
用于本发明的更好的重新初始化方案在序列号为10/729,637的 申请中进行了描述，该申请于2003年12月5日提交并且题目为 “Selectively Reduced Bi-cubic Interpolation for Ink-Jet Simulations on Quadrilateral Grids(在四边形网格上用于喷墨模 拟的有选择地减少双三次插值)”。该申请的公开内容在此被结合以 供参考。这个重新初始化方案展现了好得多的质量守恒特性。
D.关于时间步长的约束条件
由于时间积分方案是用于对流项的二阶显式以及用于粘性的一阶 显式，所以关于时间步长Δt的约束条件由CFL条件、表面张力、粘性 和总加速度确定，如方程(38)所示，
$\mathrm{\Delta t}<\underset{i,j}{\min }[\frac{\mathrm{\Delta x}}{u},\frac{\mathrm{\Delta y}}{\upsilon },\sqrt{\mathrm{We}\frac{{\rho }_1+{\rho }_2}{8\pi }}{h}^{3/2},\frac{\mathrm{Re}}{2}\frac{{\rho }^n}{{\mu }^n}{(\frac{1}{\Delta x^2}+\frac{1}{\Delta y^2})}^{-1}\sqrt{\frac{2h}{\left|F^n\right|}}],---\left(38\right)$
其中h＝min(Δx，Δy)及Fn在方程(19)中定义。
E.流程图
如在图3中的流程图所示，数值算法基本上是顺序的。代码首先 读取喷嘴几何尺寸(步骤301)，还读取比如趋向(tend)(模拟的 结束时间)、α(界面浸润的范围)、ifq-reini(水平集应当多久被 重新初始化一次)等控制参数(步骤302)。利用所给出的喷嘴几何尺 寸，代码产生均匀的矩形网格(步骤303)。当前时间步长的时间和数 量被设置为零，并且初始流体速度在处处被设置为零(步骤304)。利 用所给出的浸润参数(α)，使用方程(19)来设置界面厚度(步骤 305)。然后，通过假定初始墨水-空气界面是平坦的来初始化水平集 (步骤306)。
现在，通过检查是否t＜趋向来开始时间循环(步骤307)。如果 是的话，根据序列号为10/652,386的申请的过程来确定一致的反压 力，该申请于2003年8月29日提交，并且题目为“Consistent Back Pressure for Piezoelectric Ink-Jet Simulation(用于压电喷墨 模拟的一致的反压力)”，该申请的公开内容在此结合作为参考(步 骤308)。然后通过方程(38)确定时间步长以确保代码的稳定性(步 骤309)。更新时间(步骤310)。然后，时间步长和墨水流速(初始 流速为零)被传送至计算当前时间步长的流入压力的等效电路或类似 分析工具(步骤311)。在从等效电路接收流入速度后，CFD代码进行 求解偏微分方程。首先根据方程(22)计算速度和水平集的斜率(步 骤312)。然后，使用方程(24)计算预测算子(步骤313)。一旦获 得水平集预测算子，则还计算以时间为中心的粘性和密度。使用方程 (30)至(32)计算速度和水平集校正算子(步骤314)。对于每一 个ifq-reini时间步长，水平集还被重新隔开(re-distanced)(步 骤315和316)。使用新的水平集值计算新的流体粘性和密度(步骤 317)。速度场被投影到无发散空间中以获得新的压力和不可压缩的速 度场(步骤318)。在循环中所做的最后的事情是计算墨水流动速率(步 骤319)，以及更新时间步长的数量(步骤320)。
IV.多重网格兼容投影
在这个部分中，说明了与多重网格方法兼容的离散投影算子的构 造。在前面部分中解释的中心差分方案具有“交错的”特征；即需要 两个离散投影算子，一个具有位于单元中心的压力，以及另一个具有 位于网格点处的压力。这里我提出一个有限差分投影和一个有限元投 影。它们在模拟中被交替使用。当速度位于网格点而压力在单元中心 处时，使用有限差分投影算子。当速度位于单元中心而压力在网格点 处时，使用有限元投影算子。当采用有限差分投影时，单元密度在时 间上滞后半个时间步长，这被称为“延迟”密度。
A.t＝tn+1的投影
依据前面的部分，如果速度在t＝tn时位于单元中心，则它在t＝tn+1 的末尾时位于网格点处。由于我们想要使压力和速度之间的相对位置 保持相同，所以在投影步骤中要解出的压力即pn+1应当位于单元中心。 相对位置在图4中示出，其中在 中的“否定号”表示这些速度是作 为不是不可压缩的速度校正算子。注意，没有显示的密度校正算子也 位于网格点。
有限元方法(FEM)投影可用来计算压力并获得无发散的速度。通 过让压力成为分段线性且让速度梯度成为分段恒定，可以以直接的方 式使FEM投影离散化。IMSL直接求解程序被认为求解线性系统。然而， 因为压力位于单元中心而密度在网格点，所以就算不是不可能，应用 该快速多重网格方法来求解线性系统以改善代码性能也是非常困难 的。与多重网格方法相比，IMSL直接求解程序是非常慢的。
因此，这里提出使用具有“略微延迟的”中心密度预测算子的有 限差分投影(34)。参考图4并考虑中心差分，方程(34)可如方程 (39)所述被离散化。单元边缘密度，象pi+1/2，j n+1/2和pi，j+1/2 n+1/2是由在单元中心 平均密度预测算子来近似的。例子在方程(40)中给出。注意，通过 使用由泰勒展开式计算的水平集预测算子φi，j n+1/2来获得密度预测算子。 由于在方程(39)中的压力和速度均用于t＝tn+1，所以用于投影中的密 度滞后半个时间步长。因此，该密度被称为“延迟密度”。
$\frac{\mathrm{\Delta t}}{\Delta x^2}[\frac{1}{{\rho }_{i+1/2,j}^{n+1/2}}(p_{i+1,j}^{n+1}-p_{i,j}^{n+1})-\frac{1}{{\rho }_{i-1/2,j}^{n+1/2}}(p_{i,j}^{n+1}-p_{i-1,j}^{n+1})]$
$+\frac{\mathrm{\Delta t}}{\Delta y^2}[\frac{1}{{\rho }_{i,j+1/2}^{n+1/2}}(p_{i,j+1}^{n+1}-p_{i,j}^{n+1})-\frac{1}{{\rho }_{i,j-1/2}^{n+1/2}}(p_{i,j}^{n+1}-p_{i,j-1}^{n+1})]---\left(39\right)$
$=\frac{1}{2\mathrm{\Delta x}}({\bar{u}}_{i+1/2,j+1/2}^{n+1}+{\bar{u}}_{i+1/2,j-1/2}^{n+1}-{\bar{u}}_{i-1/2,j+1/2}^{n+1}-{\bar{u}}_{i-1/2,j-1/2}^{n+1})$
$+\frac{1}{2\mathrm{\Delta y}}({\bar{u}}_{i+1/2,j+1/2}^{n+1}+{\bar{u}}_{i-1/2,j+1/2}^{n+1}-{\bar{u}}_{i+1/2,j-1/2}^{n+1}-{\bar{u}}_{i-1/2,j-1/2}^{n+1}).$
${\rho }_{i+1/2,j}^{n+1/2}=\frac{{\rho }_{i,j}^{n+1/2}+{\rho }_{i+1,j}^{n+1/2}}{2}.---\left(40\right)$
在离散方程(39)中，边界条件是容易实现的。例如，如果固体 壁边界通过单元(i，j)和(i+1，j)之间的边缘，则需要对压力实施 Neumann条件dp/dn＝0。这可仅仅使用(41)来进行。注意，在计算 速度预测算子和校正算子时，已经考虑了速度的非滑动条件。因此， 在投影中没有考虑非滑动条件。
$p_{i+1,j}^{n+1}=p_{i,j}^{n+1}.---\left(41\right)$
在流入处，给出了流入压力，需要满足p＝pBC和dp/dn＝0。例如， 如果流入边界位于单元(i，1)和(i，0)之间的边缘，则通过设置(42)可容 易地实现流入压力。
pi，0＝2pBC-pi，1.                         (42)
再者，在进行投影时没有必要担心流入的速度条件du/dn＝0，因为 在预测算子-校正算子的计算中已经考虑到了它。
离散投影方程(39)和边界条件可通过多重网格方法来求解，这 将在接下来的分部分中进行描述。在单元中心压力被求解后，在网格 点处不可压缩的速度可通过使用方程(33)来计算。
B.t＝tn+2的投影
在t＝tn+1，速度un+1和水平集φn+1位于网格点，而压力pn+1位于单 元中心。因此，对于t＝tn+2，速度un+2和水平集φn+2将在单元中心处被 构造，而压力pn+2在网格点处被构造。
FEM投影被采用以获得不可压缩的速度和压力。在校正算子(不是 不可压缩的)和要求解的压力之间的相对位置在图5中被示出。如在 (43)中所述的FEM投影方程应该用来促进该实施。在(43)中，Γ1 表示其中流入或流出速度uBC是给定的的所有边界。
${\int }_{\Omega }\frac{\mathrm{\Delta t}}{{\rho }^{n+2}}▿p^{n+2}·▿\mathrm{\psi dx}={\int }_{\Omega }{\bar{u}}^{n+2}·▿\mathrm{\psi dx}-{\int }_{{\Gamma }_1}\psi u^{\mathrm{BC}}·\mathrm{ndS}.---\left(43\right)$
这里，需要对可使用多重网格方法工作的方程(43)进行离散化。 由于在图5中压力位于网格点而速度和密度在单元中心，因此压力可 被近似为标准的FEM分段双线性函数，以及速度可被近似为分段常数 函数。测试函数还被选择为分段双线性，并且离散值被选择以与压力 相搭配。这些选择等价于也就是说，在网格点(i+1/2，j+1/2)处的连续FEM 算子是(44)。
$\frac{{\int }_{\Omega }\frac{\mathrm{\Delta t}}{{\rho }^{n+2}}▿p^{n+2}·▿{\psi }_{i+1/2,j+1/2}\mathrm{dx}}{{\int }_{\Omega }{\psi }_{i+1/2,j+1/2}dx}=\frac{{\int }_{\Omega }{\bar{u}}^{n+2}·▿{\psi }_{i+1/2,j+1/2}\mathrm{dx}}{{\int }_{\Omega }{\psi }_{i+1/2,j+1/2}\mathrm{dx}}-\frac{{\int }_{{\Gamma }_1}{\psi }_{i+1/2,j+1/2}{u}^{\mathrm{BC}}·\mathrm{ndS}}{{\int }_{\Omega }{\psi }_{i+1/2,j+1/2}\mathrm{dx}}.---\left(44\right)$
注意，(44)是在本发明的多重网格代码中用来离散化的实际方 程。在离散化中，方程(43)的直接使用将使多重网格求解程序发散。
在均匀矩形网格中，(44)的左侧可被离散化为(45)，同时在 右侧的第一项可被写作(46)。
$\frac{\mathrm{\Delta t}}{{6\rho }_{i,j}^{n+2}}[(\frac{2}{\Delta x^2}+\frac{2}{\Delta y^2})p_{i+1/2,j+1/2}^{n+2}-(\frac{2}{\Delta x^2}-\frac{1}{\Delta y^2})p_{i-1/2,j+1/2}^{n+2}$
$+(\frac{1}{\Delta x^2}-\frac{2}{\Delta y^2})p_{i+1/2,j-1/2}^{n+2}-(\frac{1}{\Delta x^2}+\frac{1}{\Delta y^2})p_{i-1/2,j-1/2}^{n+2}]$
$+\frac{\mathrm{\Delta t}}{{6\rho }_{i+1,j}^{n+2}}[(\frac{2}{\Delta x^2}+\frac{2}{\Delta y^2})p_{i+1/2,j+1/2}^{n+2}+(\frac{1}{\Delta x^2}-\frac{2}{\Delta y^2})p_{i+1/2,j-1/2}^{n+2}$
$-(\frac{2}{\Delta x^2}-\frac{1}{\Delta y^2})p_{i+3/2,j+1/2}^{n+2}-(\frac{1}{\Delta x^2}+\frac{1}{\Delta y^2})p_{i+3/2,j-1/2}^{n+2}]$
$+\frac{\mathrm{\Delta t}}{{6\rho }_{i,j+1}^{n+2}}[(\frac{2}{\Delta x^2}+\frac{2}{\Delta y^2})p_{i+1/2,j+1/2}^{n+2}+(\frac{1}{\Delta x^2}-\frac{2}{\Delta y^2})p_{i+1/2,j+3/2}^{n+2}---\left(45\right)$
$+(\frac{2}{\Delta x^2}-\frac{1}{\Delta y^2})p_{i-1/2,j+1/2}^{n+2}-(\frac{1}{\Delta x^2}+\frac{1}{\Delta y^2})p_{i-1/2,j+3/2}^{n+2}]$
$+\frac{\mathrm{\Delta t}}{{6\rho }_{i+1,j+1}^{n+2}}[(\frac{2}{\Delta x^2}+\frac{2}{\Delta y^2})p_{i+1/2,j+1/2}^{n+2}+(\frac{1}{\Delta x^2}-\frac{2}{\Delta y^2})p_{i+1/2,j+3/2}^{n+2}$
$-(\frac{2}{\Delta x^2}-\frac{1}{\Delta y^2})p_{i+3/2,j+1/2}^{n+2}-(\frac{1}{\Delta x^2}+\frac{1}{\Delta y^2})p_{i+3/2,j+3/2}^{n+2}]$
$\frac{1}{2\mathrm{\Delta x}}({\bar{u}}_{i,j+1}^{n+2}+{\bar{u}}_{i,j}^{n+2}-{\bar{u}}_{i+1,j+1}^{n+2}-{\bar{u}}_{i+1,j}^{n+2})+\frac{1}{2\mathrm{\Delta y}}({\bar{\upsilon }}_{i+1,j}^{n+2}+{\bar{\upsilon }}_{i,j}^{n+2}-{\bar{\upsilon }}_{i+1,j+1}^{n+2}-{\bar{\upsilon }}_{i,j+1}^{n+2})---\left(46\right)$
如果没有规定的流入/流出速度，则右侧的第二项变为零。对于给 出流入速度的情况，第二项是容易在数值上计算的。例如，如果流入 速度是不变的，以及流入边界是水平的且位于求解域的底部(即在单 元(i，1)和(i，0)之间)，则第二项简化为vBC/2Δy，其中vBC是规定的流入 速度的y分量。
C.多重网格方法的实施
在多重网格方法被应用至求解具有快速改变系数的拉普拉斯方程 时，多重网格方法趋于缓慢收敛或者甚至发散。在喷墨模拟中，多达 1000∶1的密度比率必须被适应。因此，系数在墨水-空气界面的两边 是非常快速变化的。在使用多重网格方法来求解从投影得来的线性系 统时必须小心进行。由于本发明不直接适合多重网格方法本身，所以 我在编码多重网格求解程序中采用的选择在下面简单地列出：
●用于最精确级的离散投影算子如在前面两个分部分所定义。
●采用典型的多重网格V循环。
●多色Jacobi迭代被用于每一级(除了最近似级之外)来“松 驰”线性系统。余数被限制用在下一个较近似级。所述解被插值以及 添加至下一个较精确级。
●对于具有非常高密度比率(例如1000∶1)的问题，在V循环 中，对于每一级进行四个多色Jacobi迭代，以及重复8-10个V循环 来减少余数7至8个数量级。
●限制算子和插值算子是双线性的以及相邻的。
●多色Jacobi迭代或共轭梯度方法可被应用来在底部(最近似) 级求解线性系统。但是在底部级的余数不必精确至零。而它应当是比 所要求的解的精度至少小两个数量级。
●在较近似级的数值投影算子与顶部(最精确)级一样被精确地 定义，除了网格大小(Δx和Δy)增加到两倍。
●在较近似级中的单元密度可通过平均在较精确级中的这些值来 获得。
V.数值示例
作为第一个例子，考虑由第二较轻流体包围的2D流体气泡的下 落，正如由图6中模拟所说明的。假定初始气泡是直径为1的圆。它 的密度是周围的第二流体的两倍高。雷诺数被选择为400，而韦伯数被 设定为无穷大(即没有表面张力)。求解域是5×5，以及256×256 的均匀正方网格被用于模拟。时间步长是0.005。如上面所指出的，IMSL 直接求解程序或多重网格求解程序可被使用以从投影中求解线性系 统。唯一的差别是，如果使用多重网格求解程序，则投影的离散化应 当如前面部分中的来进行；否则，多重网格求解程序不可能收敛。注 意，仅仅一组结果被绘制，因为使用IMSL和多重网格求解程序所获得 的结果是如此接近，以致于没有可见的差异。在具有双Xeon 2.8GHz CPU和266MHz ECC SDRAM的Windows XP工作站上，使用一个CPU和 多重网格求解程序对于700个时间步长的模拟花费大约930秒。使用 IMSL直接求解程序的同样模拟花费大约1700秒。比较耗费在求解从 投影的离散化得到的线性系统的CPU时间，用于IMSL直接求解程序的 数字是每系统1.38秒，以及用于多重网格求解程序的数字是每系统 0.281秒。
对于第二个简单的2D例子，考虑通过图7中的动态压力的墨水滴 喷射，其中最大压力是8，而最小压力是-7。喷嘴直径在开口处是1， 而在底部(流入)处是2。喷嘴开口部分的长度，即其直径为1的部分 的长度是1。倾斜部分的长度是2.2。初始墨水空气界面被假定是平坦 的且从喷嘴开口向下0.1。在这个模拟中使用的其它参数是：雷诺数等 于50；韦伯数等于31.3。使用32×336矩形网格的模拟结果被绘制在 图8中。该模拟(从t＝0至t＝9)花费900个时间步长完成，并且在 每一个时间步长中，代码求解10750×10750的矩阵系统。在具有 266MHz ECC SDRAM的2.8GHz Xeon Windows工作站上，CPU时间是 509秒。
VI.应用和效果
正如前述所演示的，本发明提供了一种用于两相喷墨模拟的耦合 水平集投影方法的中心差分方案。本发明还提供数值投影算子，以使 得到的线性系统可使用多重网格方法来求解。
本发明的细节已被说明，现在参考图9说明可用于实施本发明一 个或多个方面的示例性系统90。如图9中所示，可以是XP Windows 工作站的该系统包括中央处理单元(CPU)91，该中央处理单元(CPU) 91提供计算资源并控制计算机。CPU 91可利用微处理器等来实施，以 及可代表多于一个的CPU(例如双Xeon 2.8GHz CPU)，并且还可包 括一个或多个诸如图形处理器之类的辅助芯片。系统90进一步包括可 以是随机存取存储器(RAM)和只读存储器(ROM)的形式的系统存储 器92。
多个控制器和外围设备还被提供，如图9所示。输入控制器93表 示诸如键盘、鼠标或记录笔(stylus)之类的各种输入设备94的接口。 存储控制器95与一个或多个存储设备96连接，该存储设备96中的每 一个包括存储介质，比如磁带或磁盘、或者光学介质，该存储介质可 用来记录用于操作系统、实用程序(utility)和应用程序的指令的程 序，其可包括执行本发明的各个方面的程序的实施例。存储设备96还 可用于存储依据本发明的处理过的或待处理的数据。显示控制器97提 供用于观看模拟的可为任何已知类型的显示设备98的接口。打印机控 制器99还被提供以用于与打印机101通信。通信控制器102与一个或 多个通信设备103连接，所述通信设备103使系统90能够通过包括因 特网、局域网(LAN)、广域网(WAN)的多种网络中的任一个或通过 包括红外信号的任何合适的电磁载波信号连接至远程设备。
在所说明的系统中，所有主要系统部件连接至可表示多于一个物 理总线的总线104。然而，各种系统部件可以在物理上彼此相邻或不相 邻。例如，输入数据和/或输出数据可以从一个物理位置远程传输至另 一位置。同样，执行本发明的各个方面的程序可从远程位置(例如服 务器)经网络来访问。这种数据和/或程序可通过包括磁带或磁盘或光 盘的多种机器可读介质、网络信号、或包括红外信号的任何其它合适 的电磁载波信号中任何一个进行传输。
本发明可方便地使用软件来实施。然而，可替换的实现方式当然 也是可行的，包括硬件和/或软件/硬件实施。硬件实施的功能可使用 ASIC、数字信号处理电路等等来实现。因此，在权利要求书中的术语 “部件或模块”用来包括软件和硬件实施。类似地，如这里使用的“机 器可读介质”包括软件、在其上具有指令硬连线的程序的硬件、或者 它们的组合。考虑这些实施方式替换的情况，可以理解的是，附图和 相应的说明提供了功能性信息，本领域技术人员将需要该功能性信息 来写出程序代码(即软件)或制作电路(即硬件)以执行所要求的处 理。
尽管结合几个具体实施例已经说明了本发明，根据前面的说明， 另外的替换、改进、变化和应用对于本领域技术人员将是显而易见的。 因此，这里说明的本发明用来包括可以处在所附权利要求书的精神和 范围内的所有这种替换、改进、变化和应用。