用于喷墨模拟的2d中心差分水平集投影方法

专利名称：用于喷墨模拟的2d中心差分水平集投影方法
技术领域：
本发明涉及一种模拟和分析来自压电打印头的喷墨的改进模型和伴随算法。对模拟模型和算法的改进包括用于两相流的耦合水平集(level set)投影方法的中心差分方案的开发，以及还有可用于中心差分方案的多重网格兼容(multi-grid-compatible)离散投影算子的构造。该模拟模型和算法可以以软件形式实施并在计算机上运行，同时可在附随的监视器上观看随时间流逝的模拟。
背景技术：
相关申请资料本申请与序列号为10/390,239的申请相关，该申请于2003年3月14日提交，并且题目为“Coupled Quadrilateral Grid Level SetScheme 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)ρ(φ)DuDt=-▿p+▿·(2μ(φ)D)-σκ(φ)δ(φ)▿φ.---(2)]]>在这些方程中，形变张量的速率和流体速度分别由方程(3)确定。
D=12[▿u+(▿u)T],---(3)]]>u＝ue1+υe2，DDt=∂∂t+(u·▿)]]>是拉格朗日时间导数，ρ是相对密度，p是压力，μ是相对动态粘性，σ是表面张力系数，к是曲率，δ是迪拉克δ函数，以及φ是水平集。
水平集公式用于定义在两个流体譬如墨水和空气之间的界面，因此相对密度、相对动态粘性和曲率全部被定义在如方程(4)所述的水平集φ的项中。
这里，水平集函数φ被初始化为到界面的带符号的距离，即水平集值在流体侧是到界面的最短距离，而在空气侧是最短距离的负值。与该界面正交、从流体2被吸入流体1的单位以及界面的曲率可以如方程(5)中φ的项来表示。
n=▿φ|▿φ||φ=0,]]>κ=▿·(▿φ|▿φ|)|φ=0.---(5)]]>为了使控制方程无量纲，选择如在方程(6)中所述的下面定义。
x＝Lx′，y＝Ly′，u＝Uu′，t=LUt′,---(6)]]>p＝ρ1U2p′，ρ＝ρ1ρ′，μ＝μ1μ′，最初量是无量纲的，并且L、U、ρ1、μ1分别表示特征长度、特征速度、流体1的密度和流体1的动态粘性。特征速度和特征长度可被任意选择，因为它们不影响模拟的结果。
将方程(6)的关系代入方程(1)和(2)，并降低该最初值，得到方程(7)和(8)，其中相对密度比率ρ(φ)、相对粘性比率μ(φ)、雷诺数Re和韦伯数均由方程(9)定义。
·u＝0， (7)∂u∂t+(u·▿)u=-1ρ(φ)▿p+1ρ(φ)Re▿·(2μ(φ)D)-1ρ(φ)Weκ(φ)δ(φ)▿φ,---(8)]]>ρ(φ)=1ifφ≥0ρ1/ρ2ifφ<0,]]>μ(φ)=1ifφ≥0μ2/μ1ifφ<0,---(9)]]>Re=ρ1ULμ1,]]>We=ρ1U2Lσ.]]>由于界面随流体移动，所以水平集的演变由方程(10)控制。
∂φ∂t+u·▿φ=0.---(10)]]>由于方程(5)、(7)和(8)以向量记法的项来表示，所以它们在笛卡儿坐标系和轴对称坐标系中采用相同的形式。
为了说明本发明，我喜欢使用对流项的守恒形式。因此，方程(8)和(10)被分别重写为方程(11)和(12)。
∂u∂t+▿·(uu)=-1ρ(φ)▿p+1ρ(φ)Re▿·(2μ(φ)D)-1ρ(φ)Weκ(φ)δ(φ)▿φ,---(11)]]>∂φ∂t+▿·(uφ)=0.---(12)]]>方程(8)、(10)和(11)、(12)的等价性可由扩展对流项来确认。例如，在方程(12)中的扩展在方程(13)中被说明。只要流体是不可压缩的，即只要·u＝0，则保持最终的等价性。
▿·(uφ)=∂(uφ)∂x+∂(υφ)∂y---(13)]]>=∂(uφ)∂x+∂(υφ)∂y]]>=φ(u,x+υ,y)+uφ,x+υφ,y]]>=φ▿·u+(u·▿)φ]]>=(u·▿)u.]]>B.边界条件在固体壁上，假定速度的法线和切线分量变为零，尽管在三相点处这个假定必须被改变。在流入和流出处，本发明的公式表示允许我们规定速度(14)或压力(15)的边界条件。在方程(15)中，符号n表示正交于流入或流出边界的单位。
u＝uBC(14)p=pBC,∂u∂n=0,---(15)]]>对于喷墨模拟，与时间相关的流入条件由模拟电荷驱动机构的等效电路模型来提供，该电荷驱动机构迫使墨水从槽中进入喷嘴。
C.接触模型在三相点处，其中空气和墨水在固体壁处相遇，如在2002年3月22日提交的且题目为“A Slipping Contact Line Model andMass-Conservative Level Set Implementation for Ink-JetSimulation(执行用于喷墨模拟的滑动接触线模型和质量守恒水平集)”的序列号为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(φ)=0ifφ<-∈12[1+φ∈+1πsin(πφ/∈)]if|φ|≤∈1ifφ>∈---(17)]]>δ(φ)=dH(φ)dφ.---(18)]]>通常选择参数ε与如方程(19)中所述的单元的平均大小成比例，其中Δx是矩形单元的平均大小。我通常选择α在1.7和2.5之间。
ε＝αΔx，(19)因此，密度和粘性变得如方程(20)所述。
ρ(φ)=ρ2ρ1+H(φ)(1-ρ2ρ1),---(20)]]>μ(φ)=μ2μ1+H(φ)(1-μ2μ1).]]>B.中心方案网格由大小为Δx×Δy的均匀矩形单元组成。这些单元以(xi＝(i-1/2)Δx，yj＝(j-1/2)Δy)为中心。在时间步长tn，假定离散速度un＝(ui，j，vi，j)和水平集φi，j均位于单元中心，而压力pi，jn在网格点处，如图1所示。目的是获得在t＝tn+1处的速度、压力和水平集。
B.1插值第一级是从位于单元中心处的离散速度和水平集值线性地重构该场，如方程(21)所述。
un(x,y)=ui,jn+ui,j′Δx(x-xi)+ui,j′Δy(y-yj),]]>φn(x,y)=φi,jn+φi,j′Δx(x-xi)+φi,j′Δy(y-yj).---(21)]]>
在这些方程中， 和 和 分别是在x和y方向上的离散速度(水平集)斜率。有几种计算这些斜率的方法。例如，一种方法可使用如在上述相关申请中的方程(60)-(62)中讨论的二阶单调受限斜率。可替换地，如在此优选的，使用在方程(22)中所述的简单的中心差分来计算这些斜率。
ui,j′Δx=ui+1,jn-ui-1,jn2Δx,ui,j′Δy=ui,j+1n-ui,j-1n2Δy,]]>φi,j′Δx=φi+1,jn-φi-1,jn2Δx,φi,j′Δy=φi,j+1n-φi,j-1n2Δy.---(22)]]>B.2预测算子第二级是获得在t＝tn+1/2处速度和水平集的预测算子。这可如按照方程(23)的由在时间上的泰勒展开式容易地实现。
ui,jn+1/2=ui,jn+Δt2∂u∂t|t=tn,]]>φi,jn+1/2=φi,jn+Δt2∂φ∂t|t=tn.---(23)]]>通过代入方程(11)和(12)得到方程(24)来估计在(23)右侧的时间导数。在单元中心处获得水平集预测算子φi，jn+1/2之后，使用方程(20)来计算以时间为中心的密度(ρ(φn+1/2))和粘性(μ(φn+1/2))。注意，这些密度和粘性在t＝tn+1/2处仍然是以单元为中心的。
ui,jn+1/2=ui,jn+Δt2{-∂pi,jn∂x-(2ui,j′Δxui,jn+ui,j′Δyυi,jn+υi,j′Δyui,jn)]]>+1ρ(φn)Re[(2μ(φ)u,x),x+(μ(φ)(u,y+υ,x)),y]t=tn]]>-1ρ(φn)Weκ(φn)δ(φn)φ,xn},]]>υi,jn+1/2=υi,jn+Δt2{-∂pi,jn∂y-(ui,j′Δxυi,jn+υi,j′Δxui,jn+2υi,j′Δyυi,jn)---(24)]]>+1ρ(φn)Re[(μ(φ)(u,y+υ,x)),x+(2μ(φ)υ,y),y]t=tn]]>-1ρ(φn)Weκ(φn)δ(φn)φ,yn},]]>φi,jn+1/2=φi,jn-Δt2{1Δx(ui,j′φi,jn+ui,jnφi,j′)+1Δy(υi,j′φi,jn+υi,jnφi,j′)}.]]>B.3校正算子第三级是将在tn处的分段线性近似(21)引中至tn+1。新的速度场和水平集首先通过如方程(25)所述的它们的交错单元平均来实现。该方案的中心差分特性被牢固地链接至在(25)中的交错的平均，其应当被集成在如图2所示的控制箱 中。以水平集对流方程为例，将方程(26)代入方程(25)的第一方程中以得到方程(27)。对于在方程(27)右侧的第一项，插入线性重构的水平集(即在(21)中的第二方程)。对于第二和第三项，使用在空间和时间上的中间点求积分。以第二项为例，该代入被显示在方程(28)中，其中以时间为中心的量是那些在前述级获得的速度和水平集预测算子，以及导数算子Dx+和平均算子μy+被定义为如方程(29)中所示。在下面，我同样使用定义的导数和平均算子比如Dy+、Dx-、μx+等，这些算子类似于(27)被定义。
φ‾i+1/2,j+1/2n+1=1ΔxΔy∫Ci+1/2,j+1/2φ(x,y,tn+1)dxdy,]]>u‾i+1/2,j+1/2n+1=1ΔxΔy∫Ci+1/2,j+1/2u(x,y,tn+1)dxdy.---(25)]]>φ=φn-∫▿·(uφ)dt---(26)]]>φ‾i-1/2,j+1/2n+1=1ΔxΔy∫Ci+1/2,j+1/2φ(x,y,tn)dxdy---(27)]]>-1ΔyDx+∫τ=tntn+1∫y∈Jj+1/2(uφ)dydτ-1ΔxDy+∫τ=tntn+1∫x∈Ii+1/2(υφ)dxdτ.]]>1ΔyDx+∫τ=tntn+1∫y∈Jj+1/2(uφ)dydτ---(28)]]>=Δt[Dx+μy+(ui,jn+1/2φi,jn+1/2)],]]>Dx+ui,j=ui+1,j-ui-1,j2Δx,μy+ui,jui,j+1+ui,j-12.---(29)]]>在插入线性重构场和预测算子之后，在方程(27)中的交错水平集平均变得如在方程(30)中所示。同样地，速度场如在方程(31)和(32)中所述。
φ‾i+1/2,j+1/2n+1=μx+μy+φi,jn-Δx8Dx+μy+φi,j′-Δy8Dy+μx+φi,j′---(30)]]>-Δt[Dx+μy+(ui,jn+1/2φi,jn+1/2)+Dy+μx+(υi,jn+1/2φi,jn+1/2)].]]>
u‾i+1/2,j+1/2n+1=1ΔxΔy∫Ci+1/2,j+1/2u(x,y,tn)dxdy]]>-1ΔyDx+∫τ=tntn+1∫y∈Jj+1/2u2dydτ-1ΔxDy+∫τ=tntn+1∫x∈Ii+1/2υudxdτ]]>+1ΔxΔy∫τ=tntn+1∫Ci+1/2,j+1/21ρ(φ)Re{(2μ(φ)u,x),x+[μ(φ)(u,y+υ,x)],y}dxdydτ]]>-1ΔxΔy∫τ=tntn+1∫Ci+1/2,j+1/21ρ(φ)Weκ(φ)δ(φ)φ,xdxdydτ]]>=μx+μy+ui,jn-Δx8Dx+μy+ui,j′-Δy8Dy+μx+ui,j′]]>-Δt[Dx+μy+(ui,jn+1/2ui,jn+1/2)+Dy+μx+(υi,jn+1/2ui,jn+1/2)]]]>+Δtμx+μy+ρ(φn+1/2)Re{(2μ(φ)u,x),x+[μ(φ)(u,y+υ,x)],y}i,jn+1/2]]>-Δtμx+μy+{1ρ(φ)Weκ(φ)δ(φ)φ,x}i,jn+1/2,---(31)]]>和υ‾i+1/2,j+1/2n+1=1ΔxΔy∫Ci+1/2,j+1/2υ(x,y,tn)dxdy]]>-1ΔyDx+∫τ=tntn+1∫v∈Jj+1/2uυdydτ-1ΔxDy+∫τ=tntn+1∫x∈Ii+1/2υ2dxdτ]]>+1ΔxΔy∫τ=tntn+1∫Ci+1/2,j+1/21ρ(φ)Re{[μ(φ)(u,y+υ,x)],x+(2μ(φ)υ,y),y}dxdydτ]]>-1ΔxΔy∫τ=tntn+1∫Ci+1/2,j+1/21ρ(φ)Weκ(φ)δ(φ)φ,ydxdydτ]]>=μx+μy+υi,jn-Δx8Dx+μy+υi,j′-Δy8Dy+μx+υi,j′]]>-Δt[Dx+μy+(ui,jn+1/2υi,jn+1/2)+Dy+μx+(υi,jn+1/2υi,jn+1/2)]]]>+Δtμx+μy+ρ(φn+1/2)Re{[μ(φ)(u,y+υ,x)],x+(2μ(φ)υ,y),y}i,jn+1/2]]>-Δtμx+μy+{1ρ(φ)Weκ(φ)δ(φ)φ,y}i,jn+1/2....(32)]]>注意，这些在t＝tn+1时的水平集和速度校正算子位于网格点处。由于不可压缩的条件没有被加强，所以由方程(31)和(32)给出的速度场不是无发散的。
B.4用于un+1的投影方程(33)中所述的在t＝tn+1时速度场un+1是不可压缩的。如果我们将发散算子应用至(33)的两边并且注意·un+1＝0，那么将获得在方程(34)中的结果。
un+1=u‾n+1-1ρn+1/2▿pn+1---(33)]]>▿·(Δtρ(φn+1/2)▿pn+1)=▿·u‾n+1.---(34)]]>投影方程(34)是椭圆形的。如果密度比率ρ(φn+1/2)是常数，则它简化为泊松方程。投影方程的有限元公式表示被显示在方程(35)中，∫Ωu*·▿ψdx=∫ΩΔtρ(φn+1/2)▿pn+1·▿ψdx+∫Γ1ψuBC·ndS,---(35)]]>其中，Γ1表示其中给出流入或流出速度uBC的全部边界。发散理论可以证明，在Γ1处的隐含的边界条件如方程(36)中所示。
Δtρ(φn+1/2)∂pn+1∂n=(u*-uBC)·n.---(36)]]>值得注意的是，如果在流入和流出处边界压力都被给定，则在(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)所述的穿越时间(crossingtime)问题。
φt′+F|φ|＝0， (37)其中F是给定的法向速度。值得注意的是，t′已被用于方程(37)中以强调它是伪时间变量，并且该方程完全为了重新初始化的目的而被求解。在F＝1的情况下，在时间上向前和向后流过界面，并计算时间t′，在时间t′处水平集函数对于各个单元改变符号。穿越时间(正的和负的)等于带符号的距离，这两个简单的方法中任何一个均适合供本发明使用。两者均已被尝试，其中在模拟结果中没有显著差别。
用于本发明的更好的重新初始化方案在序列号为10/729,637的申请中进行了描述，该申请于2003年12月5日提交并且题目为“Selectively Reduced Bi-cubic Interpolation for Ink-JetSimulations on Quadrilateral Grids(在四边形网格上用于喷墨模拟的有选择地减少双三次插值)”。该申请的公开内容在此被结合以供参考。这个重新初始化方案展现了好得多的质量守恒特性。
D.关于时间步长的约束条件由于时间积分方案是用于对流项的二阶显式以及用于粘性的一阶显式，所以关于时间步长Δt的约束条件由CFL条件、表面张力、粘性和总加速度确定，如方程(38)所示，
Δt<mini,j[Δxu,Δyυ,Weρ1+ρ28πh3/2,Re2ρnμn(1Δx2+1Δy2)-12h|Fn|],---(38)]]>其中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 BackPressure 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，jn+1/2和pi，j+1/2n+1/2是由在单元中心平均密度预测算子来近似的。例子在方程(40)中给出。注意，通过使用由泰勒展开式计算的水平集预测算子φi，jn+1/2来获得密度预测算子。由于在方程(39)中的压力和速度均用于t＝tn+1，所以用于投影中的密度滞后半个时间步长。因此，该密度被称为“延迟密度”。
ΔtΔx2[1ρi+1/2,jn+1/2(pi+1,jn+1-pi,jn+1)-1ρi-1/2,jn+1/2(pi,jn+1-pi-1,jn+1)]]]>+ΔtΔy2[1ρi,j+1/2n+1/2(pi,j+1n+1-pi,jn+1)-1ρi,j-1/2n+1/2(pi,jn+1-pi,j-1n+1)]---(39)]]>=12Δx(u‾i+1/2,j+1/2n+1+u‾i+1/2,j-1/2n+1-u‾i-1/2,j+1/2n+1-u‾i-1/2,j-1/2n+1)]]>+12Δy(u‾i+1/2,j+1/2n+1+u‾i-1/2,j+1/2n+1-u‾i+1/2,j-1/2n+1-u‾i-1/2,j-1/2n+1).]]>ρi+1/2,jn+1/2=ρi,jn+1/2+ρi+1,jn+1/22.---(40)]]>在离散方程(39)中，边界条件是容易实现的。例如，如果固体壁边界通过单元(i，j)和(i+1，j)之间的边缘，则需要对压力实施Neumann条件dp/dn＝0。这可仅仅使用(41)来进行。注意，在计算速度预测算子和校正算子时，已经考虑了速度的非滑动条件。因此，在投影中没有考虑非滑动条件。
pi+1,jn+1=pi,jn+1.---(41)]]>在流入处，给出了流入压力，需要满足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是给定的的所有边界。
∫ΩΔtρn+2▿pn+2·▿ψdx=∫Ωu‾n+2·▿ψdx-∫Γ1ψuBC·ndS.---(43)]]>这里，需要对可使用多重网格方法工作的方程(43)进行离散化。由于在图5中压力位于网格点而速度和密度在单元中心，因此压力可被近似为标准的FEM分段双线性函数，以及速度可被近似为分段常数函数。测试函数还被选择为分段双线性，并且离散值被选择以与压力相搭配。这些选择等价于也就是说，在网格点(i+1/2，j+1/2)处的连续FEM算子是(44)。
∫ΩΔtρn+2▿pn+2·▿ψi+1/2,j+1/2dx∫Ωψi+1/2,j+1/2dx=∫Ωu‾n+2·▿ψi+1/2,j+1/2dx∫Ωψi+1/2,j+1/2dx-∫Γ1ψi+1/2,j+1/2uBC·ndS∫Ωψi+1/2,j+1/2dx.---(44)]]>注意，(44)是在本发明的多重网格代码中用来离散化的实际方程。在离散化中，方程(43)的直接使用将使多重网格求解程序发散。
在均匀矩形网格中，(44)的左侧可被离散化为(45)，同时在右侧的第一项可被写作(46)。
Δt6ρi,jn+2[(2Δx2+2Δy2)pi+1/2,j+1/2n+2-(2Δx2-1Δy2)pi-1/2,j+1/2n+2]]>+(1Δx2-2Δy2)pi+1/2,j-1/2n+2-(1Δx2+1Δy2)pi-1/2,j-1/2n+2]]]>+Δt6ρi+1,jn+2[(2Δx2+2Δy2)pi+1/2,j+1/2n+2+(1Δx2-2Δy2)pi+1/2,j-1/2n+2]]>-(2Δx2-1Δy2)pi+3/2,j+1/2n+2-(1Δx2+1Δy2)pi+3/2,j-1/2n+2]]]>+Δt6ρi,j+1n+2[(2Δx2+2Δy2)pi+1/2,j+1/2n+2+(1Δx2-2Δy2)pi+1/2,j+3/2n+2---(45)]]>+(2Δx2-1Δy2)pi-1/2,j+1/2n+2-(1Δx2+1Δy2)pi-1/2,j+3/2n+2]]]>+Δt6ρi+1,j+1n+2[(2Δx2+2Δy2)pi+1/2,j+1/2n+2+(1Δx2-2Δy2)pi+1/2,j+3/2n+2]]>-(2Δx2-1Δy2)pi+3/2,j+1/2n+2-(1Δx2+1Δy2)pi+3/2,j+3/2n+2]]]>12Δx(u‾i,j+1n+2+u‾i,jn+2-u‾i+1,j+1n+2-u‾i+1,jn+2)+12Δy(υ‾i+1,jn+2+υ‾i,jn+2-υ‾i+1,j+1n+2-υ‾i,j+1n+2)---(46)]]>如果没有规定的流入/流出速度，则右侧的第二项变为零。对于给出流入速度的情况，第二项是容易在数值上计算的。例如，如果流入速度是不变的，以及流入边界是水平的且位于求解域的底部(即在单元(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.8GHzCPU和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、数字信号处理电路等等来实现。因此，在权利要求书中的术语“部件或模块”用来包括软件和硬件实施。类似地，如这里使用的“机器可读介质”包括软件、在其上具有指令硬连线的程序的硬件、或者它们的组合。考虑这些实施方式替换的情况，可以理解的是，附图和相应的说明提供了功能性信息，本领域技术人员将需要该功能性信息来写出程序代码(即软件)或制作电路(即硬件)以执行所要求的处理。
尽管结合几个具体实施例已经说明了本发明，根据前面的说明，另外的替换、改进、变化和应用对于本领域技术人员将是显而易见的。因此，这里说明的本发明用来包括可以处在所附权利要求书的精神和范围内的所有这种替换、改进、变化和应用。
权利要求
1.一种用于模拟和分析来自通道的流体喷射的方法，在流过该通道的第一流体和第二流体之间有一个边界，该方法包括下列步骤(a)在均匀矩形网格上用公式表示基于中心差分的离散化；(b)使用基于中心差分的离散化对均匀矩形网格执行有限差分分析，以求解控制通过该通道的至少第一流体的流动的方程；和(c)基于所执行的有限差分分析，模拟第一流体通过该通道的流动和从该通道的喷射。
2.如权利要求1所述的方法，其中第一流体是墨水，第二流体是空气，以及该通道包括喷墨喷嘴，该喷墨喷嘴是压电喷墨头的一部分。
3.如权利要求1所述的方法，其中在执行步骤(b)中，水平集方法被用于收集通道中的边界的特征。
4.如权利要求1所述的方法，其中均匀矩形网格包括多个单元，用公式表示基于中心差分的离散化，从而对于每一单元都有用于确定第一流体速度的速度向量和用于收集通道中的边界的特征的水平集值。
5.如权利要求4所述的方法，其中，对于每一单元，对应的速度向量和水平集值在一个时间点位于该单元的近似中心，并且在下一个时间点位于该单元的网格点。
6.如权利要求5所述的方法，其中步骤(a)包括构造多重网格兼容的投影算子，所述算子包括有限差分投影算子和有限元投影算子。
7.如权利要求6所述的方法，其中在步骤(b)中，使用基于中心差分的离散化对均匀矩形网格执行有限差分分析以求解方程，该方程控制在每个时间点在各个单元中通过通道的至少第一流体的流动，在步骤(b)后施加的该有限差分投影算子在速度向量和水平集值位于网格点的情况下被执行，并且在步骤(b)后施加的该有限元投影算子在速度向量和水平集值位于近似单元中心的情况下被执行。
8.一种用于模拟和分析来自通道的流体喷射的装置，在流过该通道的第一流体和第二流体之间有一个边界，该装置包括一个或多个部件或模块，所述部件或模块被配置用来在均匀矩形网格上用公式表示基于中心差分的离散化；使用基于中心差分的离散化对均匀矩形网格执行有限差分分析，以求解控制通过该通道的至少第一流体的流动的方程；和基于所执行的有限差分分析，模拟第一流体通过该通道的流动和从该通道的喷射。
9.如权利要求8所述的装置，其中一个或多个部件或模块的处理由以软件、硬件或其组合实现的指令的程序来规定。
10.如权利要求8所述的装置，其中一个或多个部件或模块包括用于可视化观察该模拟的显示器。
11.如权利要求8所述的装置，其中第一流体是墨水，第二流体是空气，以及该通道包括喷墨喷嘴，该喷墨喷嘴是压电喷墨头的一部分。
12.一种具有用于指导机器执行一种方法的指令程序的机器可读介质，该方法用于模拟和分析来自通道的流体喷射，在流过该通道的第一流体和第二流体之间有一个边界，该指令程序包括(a)用于在均匀矩形网格上用公式表示基于中心差分的离散化的指令；(b)用于使用基于中心差分的离散化对均匀矩形网格执行有限差分分析，以求解控制通过该通道的至少第一流体的流动的方程的指令；和(c)用于基于所执行的有限差分分析来模拟第一流体通过该通道的流动和从该通道的喷射的指令。
13.如权利要求12所述的机器可读介质，其中在执行指令(b)中，水平集方法被用于收集通道中的边界的特征。
14.如权利要求12所述的机器可读介质，其中均匀矩形网格包括多个单元，用公式表示基于中心差分的离散化，从而对于每一单元都有用于确定第一流体速度的速度向量和用于收集通道中的边界的特征的水平集值。
15.如权利要求14所述的机器可读介质，其中，对于每一单元，对应的速度向量和水平集值在一个时间点位于该单元的近似中心，并且在下一个时间点位于该单元的网格点。
16.如权利要求15所述的机器可读介质，其中指令(a)包括构造多重网格兼容的投影算子的指令，所述算子包括有限差分投影算子和有限元投影算子。
17.如权利要求16所述的机器可读介质，其中在执行指令(b)中，使用基于中心差分的离散化对均匀矩形网格执行有限差分分析以求解方程，该方程控制在每个时间点在各个单元中通过通道的至少第一流体的流动，在指令(b)后施加的该有限差分投影算子在速度向量和水平集值位于网格点的情况下被执行，并且在指令(b)后施加的该有限元投影算子在速度向量和水平集值位于近似单元中心的情况下被执行。
全文摘要
一种耦合水平集投影方法被结合入基于有限差分的喷墨模拟方法，以使该模拟更容易和更快速地执行并且减少对存储器的要求。该耦合水平集投影方法是基于在均匀矩形网格上的中心差分离散化。以有限差分投影算子或有限元投影算子的形式构造的数值投影算子可被用在中心差分方案中，在该情况中得到的线性系统可通过多重网格方法来求解。当流体速度位于网格点而压力位于单元中心时，使用有限差分投影算子，而当流体速度位于单元中心而压力位于网格点时，使用有限元投影算子。
文档编号G06F17/50GK1755702SQ200510106490
公开日2006年4月5日 申请日期2005年9月30日 优先权日2004年10月1日
发明者于峻德 申请人:精工爱普生株式会社