專利名稱：用于噴墨模擬的2d中心差分水平集投影方法
技術領域：
本發(fā)明涉及一種模擬和分析來自壓電打印頭的噴墨的改進模型和伴隨算法。對模擬模型和算法的改進包括用于兩相流的耦合水平集(level set)投影方法的中心差分方案的開發(fā)，以及還有可用于中心差分方案的多重網(wǎng)格兼容(multi-grid-compatible)離散投影算子的構造。該模擬模型和算法可以以軟件形式實施并在計算機上運行，同時可在附隨的監(jiān)視器上觀看隨時間流逝的模擬。
背景技術：
相關申請資料本申請與序列號為10/390,239的申請相關，該申請于2003年3月14日提交，并且題目為“Coupled Quadrilateral Grid Level SetScheme for Piezoelectric Ink-Jet Simulation(用于壓電噴墨模擬的耦合四邊形網(wǎng)格水平集方案)”。在此并入該相關申請的公開內容來作為參考，并且在下面的說明中將其稱作“相關申請”。
上述的相關申請公開了在四邊形網(wǎng)格上的耦合水平集投影方法。在那個發(fā)明中，速度分量和水平集位于單元中心，同時壓力在網(wǎng)格點處。Navier-Stokes方程在不考慮不可壓縮性的條件下首先在每個時間步長(time step)上被求解。然后，獲得的速度被“投影”到無發(fā)散(divergence-free)的場的空間中。在四邊形網(wǎng)格上的Godunov逆風(upwind)方案被用于評估Navier-Stokes方程中的對流項和水平集對流方程。在時間和空間上進行泰勒展開以獲得邊緣速度和水平集。由于外推法可從垂直單元邊緣的左側和右側被實施，以及可從水平邊緣的上側和下側被實施，因此在每個單元邊緣有兩個外推值(參見相關申請的方程式(50)-(52))。然后，Godunov逆風過程被使用以決定采用哪個外推值(參見相關申請的方程式(53)和(54))。在Godunov逆風過程中，局部Riemann問題實際上被解決。對于Newtonian流體，局部Riemann問題簡化為一維Burger方程式，它具有如由相關申請的方程式(53)和(54)給出的非常簡單的解。因此，基本上，局部速度方向決定了哪個外推值用于Newtonian流體。然而，如果流體不是Newtonian流體，則局部Riemann問題不可能具有簡單的解，并且因此Godunov逆風過程將非常難于編程。用于平流項(advection term)的這種逆風方案的優(yōu)點是較低的網(wǎng)格粘性和較高的數(shù)值分辨率。它的缺點是增加的代碼復雜性，該代碼復雜性是在時間和空間上進行泰勒外推和Godunov逆風過程所需要的。

發(fā)明內容
考慮到前述內容，本發(fā)明的目的是提供一種改進的模型和伴隨的算法，以模擬和分析來自壓電打印頭的噴墨。
另一目的是包括在改進的模型和算法中用于耦合水平集投影方法的中心差分方案，以及可選擇地進一步包括可用于中心差分方案的多重網(wǎng)格兼容離散投影算子。
在本發(fā)明的中心差分方案中，不考慮不可壓縮的條件，在每個時間步長中再次首先求解Navier-Stokes方程。然后，獲得的速度被“投影”到無發(fā)散的場的空間中。但是，與相關申請的逆風方案對比，這里的中心差分方案使用交錯的單元平均，并且局部速度的方向是不重要的。在時間上的泰勒外推仍被用來計算以時間為中心的(time-centered)速度和水平集預測算子。但由于它僅僅是在時間上的外推，因此代碼執(zhí)行容易得多。該方案的一個特征是，如果在時間步長t＝tn時速度和水平集位于單元中心處，則它們將在下一個時間步長t＝tn+1時位于網(wǎng)格點處。在時間步長t＝tn+2時，速度和水平集將再次位于單元中心處。在整個時間步長系列上，速度和水平集交替布置，這就是為什么該方法被描述為使用“交錯的”網(wǎng)格的原因。
為了改善該中心差分方案的性能，兩個離散投影算子被構造以將速度投影到無發(fā)散的空間中。從投影算子得到的線性系統(tǒng)可通過多重網(wǎng)格方法被快速求解。投影算子包括有限差分投影和有限元投影，它們被用在不同的條件下。當速度位于網(wǎng)格點處而壓力位于單元中心處時，使用有限差分投影算子。當速度位于單元中心處而壓力位于網(wǎng)格點處時，使用有限元投影算子。當采用有限差分投影算子時，單元密度在時間上滯后半個時間步長。這個單元密度時間延遲被稱為“延遲”密度。
如本發(fā)明所打算的，該中心差分方案的優(yōu)點包括(i)容易執(zhí)行；(ii)較低的計算機存儲要求；以及(iii)更容易推廣應用至non-Newtonian流體流動。
依據(jù)本發(fā)明的一個方面，提供了一種用于模擬和分析來自通道(channel)的流體噴射的方法。優(yōu)選地，該通道包括作為壓電噴墨頭的一部分的噴墨噴嘴。在流過該通道的第一流體(例如墨水)和第二流體(例如空氣)之間存在一個邊界。該方法包括下列步驟(a)用公式表示(formulate)在均勻矩形網(wǎng)格上基于中心差分的離散化；(b)使用基于中心差分的離散化來在均勻矩形網(wǎng)格上執(zhí)行有限差分分析，以求解控制通過該通道的至少第一流體的流動的方程；和(c)基于執(zhí)行的有限差分分析，模擬第一流體通過該通道的流動和從該通道的噴射。
優(yōu)選地，在執(zhí)行步驟(b)中，水平集方法被用于收集通道中的邊界的特征。
優(yōu)選地，均勻矩形網(wǎng)格包括多個單元，并且用公式表示基于中心差分的離散化，以使對于每一單元，存在用于確定第一流體速度的速度向量和用于收集通道中邊界的特征的水平集值。對于每一單元，對應的速度向量和水平集值在一個時間點位于該單元的近似中心，且在下一個時間點位于該單元的網(wǎng)格點。
優(yōu)選地，步驟(a)包括構造多重網(wǎng)格兼容的投影算子，該算子包括有限差分投影算子和有限元投影算子。
優(yōu)選地，在步驟(b)中，使用基于中心差分的離散化在均勻矩形網(wǎng)格上執(zhí)行有限差分分析以求解方程，該方程控制在每個時間點在各個單元中通過通道的至少第一流體的流動，在步驟(b)后施加的該有限差分投影算子在速度向量和水平集值位于網(wǎng)格點的情況下被執(zhí)行，并且在步驟(b)后施加的該有限元投影算子在速度向量和水平集值位于近似單元中心的情況下被執(zhí)行。
在另一方面，本發(fā)明包括一種用于模擬和分析來自通道的流體噴射的裝置，在流過該通道的第一流體和第二流體之間有一個邊界。該裝置包括一個或多個被配置成執(zhí)行上述處理的部件或模塊。這些處理的一些或全部可被以軟件、硬件或其組合的形式實現(xiàn)的指令程序來方便地規(guī)定。所述部件之一優(yōu)選地包括用于能夠可視化觀察模擬的顯示器。
依據(jù)本發(fā)明的另一方面，本發(fā)明的多個方面可以指令程序的形式被實現(xiàn)，該指令程序可以是軟件的形式，該軟件可以被存儲或傳送至計算機或其他受處理器控制的設備以用于執(zhí)行?？商鎿Q地，指令程序可直接以硬件(例如專用集成電路(ASIC)、數(shù)字信號處理電路等)來實現(xiàn)，或者這種指令可被實現(xiàn)為軟件和硬件的組合。
通過參考結合附圖的下面的說明和權利要求書，更全面地理解本發(fā)明的其它目的和成果將變得清楚和明白。


在附圖中，相同的附圖標記指的是相同的部件圖1是顯示在t＝tn處速度和壓力的位置的網(wǎng)格的示意圖；圖2是在本發(fā)明的實施例中使用的交錯網(wǎng)格的示意圖；圖3是說明依據(jù)本發(fā)明的實施例的算法的流程圖；圖4是說明在t＝tn+1處非無發(fā)散(non-divergence-free)速度校正算子和壓力的位置的示意圖；圖5是說明在t＝tn+2處非無發(fā)散速度校正算子和壓力的位置的示意圖；圖6是對流體泡落在較輕的周圍流體中的模擬；圖7是在噴墨模擬中被施加在噴嘴流入上的動態(tài)壓力；圖8是2D小滴噴射的模擬；以及圖9是說明可用于執(zhí)行本發(fā)明的多個方面的示例性系統(tǒng)的框圖。
具體實施例方式
I.引言本發(fā)明涉及在均勻矩形網(wǎng)格上用于尤其是噴墨的二維兩相流的模擬的中心差分水平集投影方法的開發(fā)。Navier-Stokes方程在不考慮不可壓縮性的條件下首先在每個時間步長上被求解。然后，獲得的速度被“投影”到無發(fā)散的場的空間中。與相關申請的逆風方案對比，在此的中心差分方案使用交錯單元平均，并且局部速度的方向是不重要的。為了求解Navier-Stokes方程，首先計算速度和水平集斜率。在時間上的泰勒外推被用于計算以時間為中心的速度和水平集預測算子。這些預測算子然后被用于計算速度校正算子，該速度校正算子不是無發(fā)散的。
該方案的一個特征是，如果在時間步長t＝tn處速度和水平集位于單元中心，則在下一個時間步長t＝tn+1處它們將位于網(wǎng)格點上。在時間步長t＝tn+2處，速度和水平集將再次位于單元中心。速度和水平集在整個系列時間步長上交替布置是該方法被描述為“交錯的”網(wǎng)格的原因。
為了改善該中心差分方案的性能，兩個離散投影算子被構造成將速度投影到無發(fā)散空間中。從投影算子得到的線性系統(tǒng)可通過多重網(wǎng)格方法快速求解。投影算子包括有限差分投影和有限元投影，它們被用于不同的情況中。當速度位于網(wǎng)格點處且壓力在單元中心處時，使用有限差分投影算子。當速度位于單元中心處且壓力在網(wǎng)格點處時，使用有限元投影算子。當采用有限差分投影算子時，單元密度在時間上滯后半個時間步長。這個單元密度時間延遲被稱為“延遲”密度。
本發(fā)明的中心差分算法比使用Godunov逆風方案的模擬算法更快，并具有較低的計算機開銷要求，因為(1)只使用在時間上的泰勒外推，而不是在時間和空間上的泰勒外推；(2)沒有要求解的局部Riemann問題；(3)在每個時間步長上只求解一個投影方程(在相關申請中在每個時間步長上求解兩個方程，即MAC投影和FEM投影)。
II.運動方程下面的說明解釋了如何由Kupferman和Tadmor將用于單相流的中心差分方案擴展至兩相流模擬。求解噴墨模擬的控制方程在這部分中被確定并被表述在附件中(與在這里參考的其它用數(shù)字標記的方程一起)。由于我更喜歡使用用于平流項的守恒形式，因此Navier-Stokes方程和水平集對流方程與在相關申請中給出的那些方程略有不同。用于兩相流的中心差分方案和說明算法實施例的流程圖在部分III中進行了解釋。與快速多重網(wǎng)格方法相兼容的兩個離散投影算子在部分IV中進行了解釋。在部分V中給出了模擬的例子。
A.控制方程用于兩相流的控制方程包括連續(xù)性方程(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·▿)]]>是拉格朗日時間導數(shù)，ρ是相對密度，p是壓力，μ是相對動態(tài)粘性，σ是表面張力系數(shù)，к是曲率，δ是迪拉克δ函數(shù)，以及φ是水平集。
水平集公式用于定義在兩個流體譬如墨水和空氣之間的界面，因此相對密度、相對動態(tài)粘性和曲率全部被定義在如方程(4)所述的水平集φ的項中。
這里，水平集函數(shù)φ被初始化為到界面的帶符號的距離，即水平集值在流體側是到界面的最短距離，而在空氣側是最短距離的負值。與該界面正交、從流體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的動態(tài)粘性。特征速度和特征長度可被任意選擇，因為它們不影響模擬的結果。
將方程(6)的關系代入方程(1)和(2)，并降低該最初值，得到方程(7)和(8)，其中相對密度比率ρ(φ)、相對粘性比率μ(φ)、雷諾數(shù)Re和韋伯數(shù)均由方程(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)以向量記法的項來表示，所以它們在笛卡兒坐標系和軸對稱坐標系中采用相同的形式。
為了說明本發(fā)明，我喜歡使用對流項的守恒形式。因此，方程(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.邊界條件在固體壁上，假定速度的法線和切線分量變?yōu)榱?，盡管在三相點處這個假定必須被改變。在流入和流出處，本發(fā)明的公式表示允許我們規(guī)定速度(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(執(zhí)行用于噴墨模擬的滑動接觸線模型和質量守恒水平集)”的序列號為10/105,138的申請中給出的滑動接觸線模型被使用。該申請的公開內容在此被結合以供參考。
III.中心差分方案現(xiàn)在用公式表示在均勻矩形網(wǎng)格上的數(shù)值算法。在下面，上標n(或n+1)表示時間步長，即方程(16)等。
un＝u(t＝tn) (16)假定量un、pn、φn，該算法的目的是獲得滿足不可壓縮條件的un+1、pn+1、φn+1。這里說明的顯性模式算法在時間上是一階精確的而在空間上是二階精確的。
A.浸潤(smearing)界面因為由迪拉克δ函數(shù)和由穿越界面的密度和粘性的急劇變化引起的數(shù)值困難，所以Heaviside和迪拉克δ函數(shù)由平滑的函數(shù)取代，即用來浸潤一點界面。Heaviside函數(shù)由方程(17)重新定義。因此，平滑的δ函數(shù)如在方程(18)中所述。
H(φ)=0ifφ<-∈12[1+φ∈+1πsin(πφ/∈)]if|φ|≤∈1ifφ>∈---(17)]]>δ(φ)=dH(φ)dφ.---(18)]]>通常選擇參數(shù)ε與如方程(19)中所述的單元的平均大小成比例，其中Δx是矩形單元的平均大小。我通常選擇α在1.7和2.5之間。
ε＝αΔx，(19)因此，密度和粘性變得如方程(20)所述。
ρ(φ)=ρ2ρ1+H(φ)(1-ρ2ρ1),---(20)]]>μ(φ)=μ2μ1+H(φ)(1-μ2μ1).]]>B.中心方案網(wǎng)格由大小為Δx×Δy的均勻矩形單元組成。這些單元以(xi＝(i-1/2)Δx，yj＝(j-1/2)Δy)為中心。在時間步長tn，假定離散速度un＝(ui，j，vi，j)和水平集φi，j均位于單元中心，而壓力pi，jn在網(wǎng)格點處，如圖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)中討論的二階單調受限斜率?？商鎿Q地，如在此優(yōu)選的，使用在方程(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)的由在時間上的泰勒展開式容易地實現(xiàn)。
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)右側的時間導數(shù)。在單元中心處獲得水平集預測算子φ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)所述的它們的交錯單元平均來實現(xiàn)。該方案的中心差分特性被牢固地鏈接至在(25)中的交錯的平均，其應當被集成在如圖2所示的控制箱 中。以水平集對流方程為例，將方程(26)代入方程(25)的第一方程中以得到方程(27)。對于在方程(27)右側的第一項，插入線性重構的水平集(即在(21)中的第二方程)。對于第二和第三項，使用在空間和時間上的中間點求積分。以第二項為例，該代入被顯示在方程(28)中，其中以時間為中心的量是那些在前述級獲得的速度和水平集預測算子，以及導數(shù)算子Dx+和平均算子μy+被定義為如方程(29)中所示。在下面，我同樣使用定義的導數(shù)和平均算子比如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時的水平集和速度校正算子位于網(wǎng)格點處。由于不可壓縮的條件沒有被加強，所以由方程(31)和(32)給出的速度場不是無發(fā)散的。
B.4用于un+1的投影方程(33)中所述的在t＝tn+1時速度場un+1是不可壓縮的。如果我們將發(fā)散算子應用至(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)是常數(shù)，則它簡化為泊松方程。投影方程的有限元公式表示被顯示在方程(35)中，∫Ωu*·▿ψdx=∫ΩΔtρ(φn+1/2)▿pn+1·▿ψdx+∫Γ1ψuBC·ndS,---(35)]]>其中，Γ1表示其中給出流入或流出速度uBC的全部邊界。發(fā)散理論可以證明，在Γ1處的隱含的邊界條件如方程(36)中所示。
Δtρ(φn+1/2)∂pn+1∂n=(u*-uBC)·n.---(36)]]>值得注意的是，如果在流入和流出處邊界壓力都被給定，則在(35)的右側的第二項將變?yōu)榱?。投影步驟可由離散化方程(34)或(35)來進行。
從方程(34)或(35)解出壓力pn+1之后，速度場un+1可由方程(33)獲得。
B.5交錯特性依據(jù)前述分部分，中心差分方案首先在沒有連續(xù)性的條件下積分Navier-Stokes方程。在t＝tn，速度和水平集位于單元中心處，而壓力在網(wǎng)格點處。使用在時間上的泰勒展開式來首先計算以時間為中心的速度和水平集預測算子。使用以時間為中心的水平集來計算以時間為中心的密度預測算子。然后，這些位于單元中心的預測算子被用于計算在網(wǎng)格點處的速度和水平集校正算子。由于我們想要使壓力和速度之間的相對位置保持相同，所以在投影步驟中獲得的速度un+1位于網(wǎng)格點處，而壓力pn+1在單元中心處。
從t＝tn+1到t＝tn+2的情形是存在的另一種方式。在t＝tn+1，速度位于網(wǎng)格點處，而壓力在單元中心。中心差分方案首先計算在網(wǎng)格點處的以時間為中心的速度和水平集預測算子。然后，這些預測算子被用于在單元中心處獲得校正算子。在進行投影后，不可壓縮的速度un+2位于單元中心處，而壓力pn+2在網(wǎng)格點。由于在我們的數(shù)值方案中離散速度和壓力的位置隨每一個時間步長而改變，所以它具有所謂的“交錯特性”。
C.水平集的重新初始化為了正確地捕獲界面并準確地計算表面張力，水平集需要保持為界面的帶符號的距離函數(shù)。然而，如果水平集由方程(7)更新，就不需要這樣保持。因此，代之以該模擬被周期性地停止，并且作為帶符號的距離函數(shù)的新水平集函數(shù)φ被重新產生，即|φ|＝1，而沒有改變原始水平集函數(shù)的零水平集。
前面已經認識到需要在水平集計算中這樣做，如重新初始化。用于重新初始化的直接而簡單的方法是使用輪廓繪圖儀首先找到界面(零水平集)，然后重新計算從各個單元至界面的帶符號距離。另一簡單的重新初始化選擇是解決如方程(37)所述的穿越時間(crossingtime)問題。
φt′+F|φ|＝0， (37)其中F是給定的法向速度。值得注意的是，t′已被用于方程(37)中以強調它是偽時間變量，并且該方程完全為了重新初始化的目的而被求解。在F＝1的情況下，在時間上向前和向后流過界面，并計算時間t′，在時間t′處水平集函數(shù)對于各個單元改變符號。穿越時間(正的和負的)等于帶符號的距離，這兩個簡單的方法中任何一個均適合供本發(fā)明使用。兩者均已被嘗試，其中在模擬結果中沒有顯著差別。
用于本發(fā)明的更好的重新初始化方案在序列號為10/729,637的申請中進行了描述，該申請于2003年12月5日提交并且題目為“Selectively Reduced Bi-cubic Interpolation for Ink-JetSimulations on Quadrilateral Grids(在四邊形網(wǎng)格上用于噴墨模擬的有選擇地減少雙三次插值)”。該申請的公開內容在此被結合以供參考。這個重新初始化方案展現(xiàn)了好得多的質量守恒特性。
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中的流程圖所示，數(shù)值算法基本上是順序的。代碼首先讀取噴嘴幾何尺寸(步驟301)，還讀取比如趨向(tend)(模擬的結束時間)、α(界面浸潤的范圍)、ifq-reini(水平集應當多久被重新初始化一次)等控制參數(shù)(步驟302)。利用所給出的噴嘴幾何尺寸，代碼產生均勻的矩形網(wǎng)格(步驟303)。當前時間步長的時間和數(shù)量被設置為零，并且初始流體速度在處處被設置為零(步驟304)。利用所給出的浸潤參數(shù)(α)，使用方程(19)來設置界面厚度(步驟305)。然后，通過假定初始墨水-空氣界面是平坦的來初始化水平集(步驟306)。
現(xiàn)在，通過檢查是否t＜趨向來開始時間循環(huán)(步驟307)。如果是的話，根據(jù)序列號為10/652,386的申請的過程來確定一致的反壓力，該申請于2003年8月29日提交，并且題目為“Consistent BackPressure for Piezoelectric Ink-Jet Simulation(用于壓電噴墨模擬的一致的反壓力)”，該申請的公開內容在此結合作為參考(步驟308)。然后通過方程(38)確定時間步長以確保代碼的穩(wěn)定性(步驟309)。更新時間(步驟310)。然后，時間步長和墨水流速(初始流速為零)被傳送至計算當前時間步長的流入壓力的等效電路或類似分析工具(步驟311)。在從等效電路接收流入速度后，CFD代碼進行求解偏微分方程。首先根據(jù)方程(22)計算速度和水平集的斜率(步驟312)。然后，使用方程(24)計算預測算子(步驟313)。一旦獲得水平集預測算子，則還計算以時間為中心的粘性和密度。使用方程(30)至(32)計算速度和水平集校正算子(步驟314)。對于每一個ifq-reini時間步長，水平集還被重新隔開(re-distanced)(步驟315和316)。使用新的水平集值計算新的流體粘性和密度(步驟317)。速度場被投影到無發(fā)散空間中以獲得新的壓力和不可壓縮的速度場(步驟318)。在循環(huán)中所做的最后的事情是計算墨水流動速率(步驟319)，以及更新時間步長的數(shù)量(步驟320)。
IV.多重網(wǎng)格兼容投影在這個部分中，說明了與多重網(wǎng)格方法兼容的離散投影算子的構造。在前面部分中解釋的中心差分方案具有“交錯的”特征；即需要兩個離散投影算子，一個具有位于單元中心的壓力，以及另一個具有位于網(wǎng)格點處的壓力。這里我提出一個有限差分投影和一個有限元投影。它們在模擬中被交替使用。當速度位于網(wǎng)格點而壓力在單元中心處時，使用有限差分投影算子。當速度位于單元中心而壓力在網(wǎng)格點處時，使用有限元投影算子。當采用有限差分投影時，單元密度在時間上滯后半個時間步長，這被稱為“延遲”密度。
A.t＝tn+1的投影依據(jù)前面的部分，如果速度在t＝tn時位于單元中心，則它在t＝tn+1的末尾時位于網(wǎng)格點處。由于我們想要使壓力和速度之間的相對位置保持相同，所以在投影步驟中要解出的壓力即pn+1應當位于單元中心。相對位置在圖4中示出，其中在 中的“否定號”表示這些速度是作為不是不可壓縮的速度校正算子。注意，沒有顯示的密度校正算子也位于網(wǎng)格點。
有限元方法(FEM)投影可用來計算壓力并獲得無發(fā)散的速度。通過讓壓力成為分段線性且讓速度梯度成為分段恒定，可以以直接的方式使FEM投影離散化。IMSL直接求解程序被認為求解線性系統(tǒng)。然而，因為壓力位于單元中心而密度在網(wǎng)格點，所以就算不是不可能，應用該快速多重網(wǎng)格方法來求解線性系統(tǒng)以改善代碼性能也是非常困難的。與多重網(wǎng)格方法相比，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)中，邊界條件是容易實現(xiàn)的。例如，如果固體壁邊界通過單元(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)可容易地實現(xiàn)流入壓力。
pi，0＝2pBC-pi，1. (42)再者，在進行投影時沒有必要擔心流入的速度條件du/dn＝0，因為在預測算子-校正算子的計算中已經考慮到了它。
離散投影方程(39)和邊界條件可通過多重網(wǎng)格方法來求解，這將在接下來的分部分中進行描述。在單元中心壓力被求解后，在網(wǎng)格點處不可壓縮的速度可通過使用方程(33)來計算。
B.t＝tn+2的投影在t＝tn+1，速度un+1和水平集φn+1位于網(wǎng)格點，而壓力pn+1位于單元中心。因此，對于t＝tn+2，速度un+2和水平集φn+2將在單元中心處被構造，而壓力pn+2在網(wǎng)格點處被構造。
FEM投影被采用以獲得不可壓縮的速度和壓力。在校正算子(不是不可壓縮的)和要求解的壓力之間的相對位置在圖5中被示出。如在(43)中所述的FEM投影方程應該用來促進該實施。在(43)中，Γ1表示其中流入或流出速度uBC是給定的的所有邊界。
∫ΩΔtρn+2▿pn+2·▿ψdx=∫Ωu‾n+2·▿ψdx-∫Γ1ψuBC·ndS.---(43)]]>這里，需要對可使用多重網(wǎng)格方法工作的方程(43)進行離散化。由于在圖5中壓力位于網(wǎng)格點而速度和密度在單元中心，因此壓力可被近似為標準的FEM分段雙線性函數(shù)，以及速度可被近似為分段常數(shù)函數(shù)。測試函數(shù)還被選擇為分段雙線性，并且離散值被選擇以與壓力相搭配。這些選擇等價于也就是說，在網(wǎng)格點(i+1/2，j+1/2)處的連續(xù)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)是在本發(fā)明的多重網(wǎng)格代碼中用來離散化的實際方程。在離散化中，方程(43)的直接使用將使多重網(wǎng)格求解程序發(fā)散。
在均勻矩形網(wǎng)格中，(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)]]>如果沒有規(guī)定的流入/流出速度，則右側的第二項變?yōu)榱恪τ诮o出流入速度的情況，第二項是容易在數(shù)值上計算的。例如，如果流入速度是不變的，以及流入邊界是水平的且位于求解域的底部(即在單元(i，1)和(i，0)之間)，則第二項簡化為vBC/2Δy，其中vBC是規(guī)定的流入速度的y分量。
C.多重網(wǎng)格方法的實施在多重網(wǎng)格方法被應用至求解具有快速改變系數(shù)的拉普拉斯方程時，多重網(wǎng)格方法趨于緩慢收斂或者甚至發(fā)散。在噴墨模擬中，多達1000∶1的密度比率必須被適應。因此，系數(shù)在墨水-空氣界面的兩邊是非常快速變化的。在使用多重網(wǎng)格方法來求解從投影得來的線性系統(tǒng)時必須小心進行。由于本發(fā)明不直接適合多重網(wǎng)格方法本身，所以我在編碼多重網(wǎng)格求解程序中采用的選擇在下面簡單地列出●用于最精確級的離散投影算子如在前面兩個分部分所定義。
●采用典型的多重網(wǎng)格V循環(huán)。
●多色Jacobi迭代被用于每一級(除了最近似級之外)來“松馳”線性系統(tǒng)。余數(shù)被限制用在下一個較近似級。所述解被插值以及添加至下一個較精確級。
●對于具有非常高密度比率(例如1000∶1)的問題，在V循環(huán)中，對于每一級進行四個多色Jacobi迭代，以及重復8-10個V循環(huán)來減少余數(shù)7至8個數(shù)量級。
●限制算子和插值算子是雙線性的以及相鄰的。
●多色Jacobi迭代或共軛梯度方法可被應用來在底部(最近似)級求解線性系統(tǒng)。但是在底部級的余數(shù)不必精確至零。而它應當是比所要求的解的精度至少小兩個數(shù)量級。
●在較近似級的數(shù)值投影算子與頂部(最精確)級一樣被精確地定義，除了網(wǎng)格大小(Δx和Δy)增加到兩倍。
●在較近似級中的單元密度可通過平均在較精確級中的這些值來獲得。
V.數(shù)值示例作為第一個例子，考慮由第二較輕流體包圍的2D流體氣泡的下落，正如由圖6中模擬所說明的。假定初始氣泡是直徑為1的圓。它的密度是周圍的第二流體的兩倍高。雷諾數(shù)被選擇為400，而韋伯數(shù)被設定為無窮大(即沒有表面張力)。求解域是5×5，以及256×256的均勻正方網(wǎng)格被用于模擬。時間步長是0.005。如上面所指出的，IMSL直接求解程序或多重網(wǎng)格求解程序可被使用以從投影中求解線性系統(tǒng)。唯一的差別是，如果使用多重網(wǎng)格求解程序，則投影的離散化應當如前面部分中的來進行；否則，多重網(wǎng)格求解程序不可能收斂。注意，僅僅一組結果被繪制，因為使用IMSL和多重網(wǎng)格求解程序所獲得的結果是如此接近，以致于沒有可見的差異。在具有雙Xeon 2.8GHzCPU和266MHz ECC SDRAM的Windows XP工作站上，使用一個CPU和多重網(wǎng)格求解程序對于700個時間步長的模擬花費大約930秒。使用IMSL直接求解程序的同樣模擬花費大約1700秒。比較耗費在求解從投影的離散化得到的線性系統(tǒng)的CPU時間，用于IMSL直接求解程序的數(shù)字是每系統(tǒng)1.38秒，以及用于多重網(wǎng)格求解程序的數(shù)字是每系統(tǒng)0.281秒。
對于第二個簡單的2D例子，考慮通過圖7中的動態(tài)壓力的墨水滴噴射，其中最大壓力是8，而最小壓力是-7。噴嘴直徑在開口處是1，而在底部(流入)處是2。噴嘴開口部分的長度，即其直徑為1的部分的長度是1。傾斜部分的長度是2.2。初始墨水空氣界面被假定是平坦的且從噴嘴開口向下0.1。在這個模擬中使用的其它參數(shù)是雷諾數(shù)等于50；韋伯數(shù)等于31.3。使用32×336矩形網(wǎng)格的模擬結果被繪制在圖8中。該模擬(從t＝0至t＝9)花費900個時間步長完成，并且在每一個時間步長中，代碼求解10750×10750的矩陣系統(tǒng)。在具有266MHz ECC SDRAM的2.8GHz Xeon Windows工作站上，CPU時間是509秒。
VI.應用和效果正如前述所演示的，本發(fā)明提供了一種用于兩相噴墨模擬的耦合水平集投影方法的中心差分方案。本發(fā)明還提供數(shù)值投影算子，以使得到的線性系統(tǒng)可使用多重網(wǎng)格方法來求解。
本發(fā)明的細節(jié)已被說明，現(xiàn)在參考圖9說明可用于實施本發(fā)明一個或多個方面的示例性系統(tǒng)90。如圖9中所示，可以是XP Windows工作站的該系統(tǒng)包括中央處理單元(CPU)91，該中央處理單元(CPU)91提供計算資源并控制計算機。CPU 91可利用微處理器等來實施，以及可代表多于一個的CPU(例如雙Xeon 2.8GHz CPU)，并且還可包括一個或多個諸如圖形處理器之類的輔助芯片。系統(tǒng)90進一步包括可以是隨機存取存儲器(RAM)和只讀存儲器(ROM)的形式的系統(tǒng)存儲器92。
多個控制器和外圍設備還被提供，如圖9所示。輸入控制器93表示諸如鍵盤、鼠標或記錄筆(stylus)之類的各種輸入設備94的接口。存儲控制器95與一個或多個存儲設備96連接，該存儲設備96中的每一個包括存儲介質，比如磁帶或磁盤、或者光學介質，該存儲介質可用來記錄用于操作系統(tǒng)、實用程序(utility)和應用程序的指令的程序，其可包括執(zhí)行本發(fā)明的各個方面的程序的實施例。存儲設備96還可用于存儲依據(jù)本發(fā)明的處理過的或待處理的數(shù)據(jù)。顯示控制器97提供用于觀看模擬的可為任何已知類型的顯示設備98的接口。打印機控制器99還被提供以用于與打印機101通信。通信控制器102與一個或多個通信設備103連接，所述通信設備103使系統(tǒng)90能夠通過包括因特網(wǎng)、局域網(wǎng)(LAN)、廣域網(wǎng)(WAN)的多種網(wǎng)絡中的任一個或通過包括紅外信號的任何合適的電磁載波信號連接至遠程設備。
在所說明的系統(tǒng)中，所有主要系統(tǒng)部件連接至可表示多于一個物理總線的總線104。然而，各種系統(tǒng)部件可以在物理上彼此相鄰或不相鄰。例如，輸入數(shù)據(jù)和/或輸出數(shù)據(jù)可以從一個物理位置遠程傳輸至另一位置。同樣，執(zhí)行本發(fā)明的各個方面的程序可從遠程位置(例如服務器)經網(wǎng)絡來訪問。這種數(shù)據(jù)和/或程序可通過包括磁帶或磁盤或光盤的多種機器可讀介質、網(wǎng)絡信號、或包括紅外信號的任何其它合適的電磁載波信號中任何一個進行傳輸。
本發(fā)明可方便地使用軟件來實施。然而，可替換的實現(xiàn)方式當然也是可行的，包括硬件和/或軟件/硬件實施。硬件實施的功能可使用ASIC、數(shù)字信號處理電路等等來實現(xiàn)。因此，在權利要求書中的術語“部件或模塊”用來包括軟件和硬件實施。類似地，如這里使用的“機器可讀介質”包括軟件、在其上具有指令硬連線的程序的硬件、或者它們的組合?？紤]這些實施方式替換的情況，可以理解的是，附圖和相應的說明提供了功能性信息，本領域技術人員將需要該功能性信息來寫出程序代碼(即軟件)或制作電路(即硬件)以執(zhí)行所要求的處理。
盡管結合幾個具體實施例已經說明了本發(fā)明，根據(jù)前面的說明，另外的替換、改進、變化和應用對于本領域技術人員將是顯而易見的。因此，這里說明的本發(fā)明用來包括可以處在所附權利要求書的精神和范圍內的所有這種替換、改進、變化和應用。
權利要求
1.一種用于模擬和分析來自通道的流體噴射的方法，在流過該通道的第一流體和第二流體之間有一個邊界，該方法包括下列步驟(a)在均勻矩形網(wǎng)格上用公式表示基于中心差分的離散化；(b)使用基于中心差分的離散化對均勻矩形網(wǎng)格執(zhí)行有限差分分析，以求解控制通過該通道的至少第一流體的流動的方程；和(c)基于所執(zhí)行的有限差分分析，模擬第一流體通過該通道的流動和從該通道的噴射。
2.如權利要求1所述的方法，其中第一流體是墨水，第二流體是空氣，以及該通道包括噴墨噴嘴，該噴墨噴嘴是壓電噴墨頭的一部分。
3.如權利要求1所述的方法，其中在執(zhí)行步驟(b)中，水平集方法被用于收集通道中的邊界的特征。
4.如權利要求1所述的方法，其中均勻矩形網(wǎng)格包括多個單元，用公式表示基于中心差分的離散化，從而對于每一單元都有用于確定第一流體速度的速度向量和用于收集通道中的邊界的特征的水平集值。
5.如權利要求4所述的方法，其中，對于每一單元，對應的速度向量和水平集值在一個時間點位于該單元的近似中心，并且在下一個時間點位于該單元的網(wǎng)格點。
6.如權利要求5所述的方法，其中步驟(a)包括構造多重網(wǎng)格兼容的投影算子，所述算子包括有限差分投影算子和有限元投影算子。
7.如權利要求6所述的方法，其中在步驟(b)中，使用基于中心差分的離散化對均勻矩形網(wǎng)格執(zhí)行有限差分分析以求解方程，該方程控制在每個時間點在各個單元中通過通道的至少第一流體的流動，在步驟(b)后施加的該有限差分投影算子在速度向量和水平集值位于網(wǎng)格點的情況下被執(zhí)行，并且在步驟(b)后施加的該有限元投影算子在速度向量和水平集值位于近似單元中心的情況下被執(zhí)行。
8.一種用于模擬和分析來自通道的流體噴射的裝置，在流過該通道的第一流體和第二流體之間有一個邊界，該裝置包括一個或多個部件或模塊，所述部件或模塊被配置用來在均勻矩形網(wǎng)格上用公式表示基于中心差分的離散化；使用基于中心差分的離散化對均勻矩形網(wǎng)格執(zhí)行有限差分分析，以求解控制通過該通道的至少第一流體的流動的方程；和基于所執(zhí)行的有限差分分析，模擬第一流體通過該通道的流動和從該通道的噴射。
9.如權利要求8所述的裝置，其中一個或多個部件或模塊的處理由以軟件、硬件或其組合實現(xiàn)的指令的程序來規(guī)定。
10.如權利要求8所述的裝置，其中一個或多個部件或模塊包括用于可視化觀察該模擬的顯示器。
11.如權利要求8所述的裝置，其中第一流體是墨水，第二流體是空氣，以及該通道包括噴墨噴嘴，該噴墨噴嘴是壓電噴墨頭的一部分。
12.一種具有用于指導機器執(zhí)行一種方法的指令程序的機器可讀介質，該方法用于模擬和分析來自通道的流體噴射，在流過該通道的第一流體和第二流體之間有一個邊界，該指令程序包括(a)用于在均勻矩形網(wǎng)格上用公式表示基于中心差分的離散化的指令；(b)用于使用基于中心差分的離散化對均勻矩形網(wǎng)格執(zhí)行有限差分分析，以求解控制通過該通道的至少第一流體的流動的方程的指令；和(c)用于基于所執(zhí)行的有限差分分析來模擬第一流體通過該通道的流動和從該通道的噴射的指令。
13.如權利要求12所述的機器可讀介質，其中在執(zhí)行指令(b)中，水平集方法被用于收集通道中的邊界的特征。
14.如權利要求12所述的機器可讀介質，其中均勻矩形網(wǎng)格包括多個單元，用公式表示基于中心差分的離散化，從而對于每一單元都有用于確定第一流體速度的速度向量和用于收集通道中的邊界的特征的水平集值。
15.如權利要求14所述的機器可讀介質，其中，對于每一單元，對應的速度向量和水平集值在一個時間點位于該單元的近似中心，并且在下一個時間點位于該單元的網(wǎng)格點。
16.如權利要求15所述的機器可讀介質，其中指令(a)包括構造多重網(wǎng)格兼容的投影算子的指令，所述算子包括有限差分投影算子和有限元投影算子。
17.如權利要求16所述的機器可讀介質，其中在執(zhí)行指令(b)中，使用基于中心差分的離散化對均勻矩形網(wǎng)格執(zhí)行有限差分分析以求解方程，該方程控制在每個時間點在各個單元中通過通道的至少第一流體的流動，在指令(b)后施加的該有限差分投影算子在速度向量和水平集值位于網(wǎng)格點的情況下被執(zhí)行，并且在指令(b)后施加的該有限元投影算子在速度向量和水平集值位于近似單元中心的情況下被執(zhí)行。
全文摘要
一種耦合水平集投影方法被結合入基于有限差分的噴墨模擬方法，以使該模擬更容易和更快速地執(zhí)行并且減少對存儲器的要求。該耦合水平集投影方法是基于在均勻矩形網(wǎng)格上的中心差分離散化。以有限差分投影算子或有限元投影算子的形式構造的數(shù)值投影算子可被用在中心差分方案中，在該情況中得到的線性系統(tǒng)可通過多重網(wǎng)格方法來求解。當流體速度位于網(wǎng)格點而壓力位于單元中心時，使用有限差分投影算子，而當流體速度位于單元中心而壓力位于網(wǎng)格點時，使用有限元投影算子。
文檔編號G06F17/50GK1755702SQ200510106490
公開日2006年4月5日 申請日期2005年9月30日 優(yōu)先權日2004年10月1日
發(fā)明者于峻德 申請人:精工愛普生株式會社