具有切換邏輯的局部/局部和混合的局部/全局插值的製作方法

2023-07-14 02:38:56

專利名稱：具有切換邏輯的局部/局部和混合的局部/全局插值的製作方法
技術領域：
本發明涉及用於對列印頭的噴墨進行建模、模擬和分析的系統與方法。更具體的，本發明的實施例包括用於基於有限差分的噴墨模擬的四邊形網格，其中，設計了一種算法來求解一組用於四邊形網格上的兩相流的偏微分公式。可以用軟體、硬體或它們的組合來實現上述仿真模型，並且，可以在計算機或其他處理器控制的設備上實施該模型。
背景技術：
噴墨列印頭是一種列印設備，該設備通過將墨滴噴到印刷媒質上來產生圖象。對確保任何由列印頭創建的產品的質量而言，對噴墨過程以及隨之產生的墨滴進行控制都是必不可少的。為實現這樣的控制，重要的是對列印和噴墨過程進行準確和有效的模擬。模擬這個過程包括對至少兩種流體(即墨水和空氣)和這些流體之間的界面進行建模。現有技術的方法已使用計算流體動力學、有限元分析、有限差分分析和水平集方法來對該行為進行建模。
水平集方法是一種用於捕獲界面(如列印頭噴嘴中的空氣和墨水之間的界面)的有效技術。為保持使用水平集方法的模擬的準確性和穩定性，應周期性地停止模擬和對水平集進行重新設距(re-distanced)。現有技術的方法使用雙三次插值法，簡化的雙三次插值法和三角剖分的快速行進法來對水平集進行重新設距。當該界面具有銳角轉角(sharp corner)時，現有技術的方法難以進行上述的重新設距。

發明內容
本發明的一個目標是提供一種解決了上述問題的用於對噴墨進行模擬和分析的方法。從而，便可以更精確的控制墨滴的大小和形狀。
本發明是一種用於模擬經過通道和離開通道的流體流的系統或方法。該流體流包含第一流體和第二流體之間的界面。創建表示該通道的物理空間和通道周圍的一部分物理空間的網格。創建包括一組值的水平集。每個值與網格中的一點相關。該組值中的每一個值與相關點至上述界面的最短距離成比例。求解一組描述了上述第一流體、第二流體和界面的各個方面的方程。使用以下重新設距方法中的一種或更多種方法對水平集中的特定值進行重新設距雙三次插值法、全局插值法和局部線性插值法。對於每個特定值，根據一個或多個切換規則來在上述重新設距方法之間進行切換。通過結合附圖的以下說明和權利要求，可以更清楚地理解本發明的其他目標和業績，並可以完全理解本發明。


在附圖中，相似的附圖標記表示相似的部分圖1示出了典型的噴墨列印頭噴嘴；圖2示出了可用於噴墨模擬的邊界經過擬合的四邊形網格；圖3的一系列模擬結果示出了墨滴的噴射情況；圖4示出了由模擬形成的贗像(artifact)；圖5A和5B示出了均勻網格上的一些相對的點，在這些點處，可計算各個變量；圖6示出了一種數值算法的流程圖，在其中實施了本發明的實施例；圖7是關於一部分均勻網格的圖，它示出了根據雙三次插值法而使用的單元；圖8A-8E是關於偏移四邊形網格的一些部分的圖，該網格包括某界面的一些部分，在這些圖中，示出了一些局部插值的實例；圖9是一部分偏移四邊形網格的圖；圖10是一種執行局部線性插值的方法的流程圖；圖11是一種使用輪廓繪製算法的結果來對水平集的值進行重新設距的方法的流程圖；圖12是一部分偏移四邊形網格的圖，該網格包括某界面的一部分，在該界面中，可以使用給出輪廓繪製算法的結果的重新設距的實例；圖13是一部分四邊形網格的圖，該圖中示出了逼近界面的一部分的分段1-7；這些分段與它們之前的分段的夾角以及叉積結果，所述叉積結果表示分段相對於其之前的分段的朝向；以及圖14的框圖示出了可用於實施本發明的各個方面的典型系統。
具體實施例方式
I.引言圖1示出了典型的噴墨列印頭噴嘴100，該噴嘴包括墨水102以及墨水102與空氣之間的界面104。對墨水102施加壓力脈衝引起在界面104處形成墨滴。可通過向經壓力板連接到墨水102的壓電(PZT)致動器施加動態電壓來產生上述壓力脈衝。
在設計列印頭噴嘴100時，用計算流體動力學(CFD)代碼模擬墨滴的產生是有用的。CFD代碼通過求解一組偏微分控制方程(即用於兩相流的不可壓縮的納維-斯託克斯公式)得出流體速度、壓力和界面位置來模擬墨滴的產生。納維-斯託克斯公式是示例性的，也可以使用描述墨滴行為的其他公式而不改變要求保護的本發明的範圍。圖2示出了經過體擬合的四邊形網格200的實例，在本發明的一個實施例中，CFD代碼可以用該網格來表示其中產生墨滴的物理空間。
圖3的一系列快照表示CFD代碼的結果。時間點t0處的快照302表示列印頭噴嘴100的初始狀態。時間點t1處的快照304表示列印頭噴嘴100開始產生墨滴。快照306表示正在生成墨滴時的稍後時間點t2處的列印頭噴嘴100，快照306包含贗像308。贗像308並不表示墨水102的行為，而是由CFD代碼產生的贗像。圖4是贗像308的放大圖。本發明的一個目的是減小這類贗像的尺寸或消除這類贗像的出現。
界面104由水平集函數φ的零水平集Γ表示。將水平集函數初始化為到界面104的帶符號距離。即水平集的值是墨水102側的到界面的最短距離，且是空氣側的最短距離的負值。在網格200的每個節點處計算了水平集函數φ。
II.控制公式用於兩相流的控制公式包含連續性公式、納維-斯託克斯公式和水平集對流公式(1)，在附錄中和以下討論中引用的其他編號的公式中說明了這些公式。在這些公式中，u是速度矢量、t是時間、ρ是相對密度、p是壓力、μ是相對動力學粘度、Re是雷諾數、We是韋伯數、Fr是弗魯德數、κ是曲率、δ是狄拉克delta函數，D是形變率張量。相對密度ρ、相對動力學粘度μ以及曲率κ都是用水平集函數φ定義的。
典型的列印頭噴嘴100是旋轉對稱的。因此，在柱坐標系中對列印頭噴嘴100進行建模是有利的。假設結果與方位角無關是合理的。因此，可以將模型從三維簡化至軸對稱物理空間X＝(r，z)。在該物理空間X中創建了經過體擬合的四邊形網格200。可通過變換Φ來對該物理空間進行變換，其中，該變換將物理空間X中的節點映射到計算空間Ξ＝(ξ，η)。在該計算空間Ξ中，可以用有限差分分析來求解網格200上的控制公式。可以在附錄中找到該變換的雅可比矩陣J(2)和變換矩陣T(3)。對於軸對稱坐標系，g＝2πr。求算計算空間Ξ中的連續性公式、納維-斯託克斯公式以及水平集對流公式(1)包括藉助雅可比矩陣J(2)和變換矩陣T(3)將公式(1)變換成公式(4)。在計算空間Ξ中進行求算時，由公式(5)給出納維-斯託克斯公式(4)右側的粘度項。
III.數值算法數值算法是在四邊形網格中形成的。在下文中，上標n(或n+1)表示時間步長。對於給出的量un，pn和φn，算法的目的是獲得滿足控制公式的un+1，pn+1和φn+1。本發明的實施例可包含在時間上一階準確和在空間上二階準確的算法。本發明不限於包含這種準確度的算法，而是可以包含其他形式。
A.界面的模糊在典型的系統中，兩種流體的密度和粘度是非常不同的。穿越界面104的密度和粘度的突然變化會造成模擬上的困難。因此，這種突然的變化將由平滑的函數取代。這種平滑發生在以界面104為中心的、厚度為2ε的空間區域內。從而，ε表示界面模糊的範圍。這種平滑函數的實例是公式(6)中表示的經過平滑的階躍函數和經過平滑的狄拉克delta函數。從而，將穿越界面104的突然變化的相對密度平滑為(7)，其中，ρ2/ρ1是第二流體(空氣)與第一流體(墨水)的密度之比。ε的典型值可以是網格200中的四邊形單元的平均大小的1.7至2.5倍。
B.水平集、速度場和壓力場更新可以用公式(8)計算水平集函數φn+1。可以用僅需時間點tn處得到的數據的顯式預測校正法來求算公式(8)中的以時間為中心(time-centerd)的平流項。
可以用由公式(4)得到的公式(9)計算中間值u*。在此之後，可以使用公式(10)計算壓力場pn+1。可以使用公式(11)計算新的速度矢量場un+1。
C.水平集的重新初始化為正確地捕獲界面104和準確地計算界面張力，水平集需要保持為至界面104的帶符號的距離函數。然而，如果通過公式(8)來更新該水平集，則該水平集不能保持為這種狀態。因此，周期性地停止模擬和重新創建新的水平集函數φ。下面將給出關於這種重新初始化的更為詳細的描述。
D.空間離散化如圖5A和5B所示，變換X＝Φ(Ξ)使得計算空間中的計算網格500由單位正方形組成，Δξ＝Δη＝1。圖2的經過體擬合的四邊形網格被映射到具有單位正方形的均勻計算網格500。
在以下的所有討論中，變量是在一個或多個網格的節點上計算的。節點的相對位置由圓括號中的或表示為下標的變量對(Z(x，y)Zx，y)表示。
圖5A和5B示出了均勻計算網格500的一部分。如圖5A所示，在整數時間步長處，在均勻計算網格的各個單元的中心處計算了速度分量un(i，j)和水平集的值φn(i，j)。同樣如圖5A所示，在整數時間步長處，在均勻計算網格的節點處計算了壓力pn(i，j)。圖5B示出了在半時間步長處和各邊緣的中點處計算的以時間為中心的邊緣速度un+1/2(i+1/2，j)和水平集預測值φn+1/2(i+1/2，j)。
注意，在該算法中僅需要變換(或映射)的局部定義。全局變換X＝Φ(Ξ)的存在或準確形式不重要。
E.平流項用於平流項的算法基於John B.BELL等人在「用於可變密度流的二階投影方法(A second-order Projection Method for Variable-DensityFlows)」(Journal of Computational Physics，101，334-348，1992，用於兩相(兩個恆定密度)流)中描述的非拆分的二階Godunov方法。它是一種以單元為中心的預測校正方法。
為求算上述平流項，使用了公式(12)和(13)。在這些公式中，通過空間和時間上的泰勒展開得到邊緣速度和邊緣水平集。速度的時間導數由納維-斯託克斯公式代替。水平集的時間導數由水平集對流公式代替。外插是從上述邊緣的兩側進行的，然後，用Godunov類型的上風法選擇使用哪個外插結果。
使用以下方法計算了項un+1/2(i+1/2，j)。從左進行外插，得到公式(14)，其中，在公式(15)中給出了Fn(i，j)。從右進行外插，得到公式(16)。下一步是Godunov上風法。平流速度由公式(17)和(18)決定。其他以時間為中心的邊緣值可以用類似方式進行計算。
F.流程圖可以將本發明的實施例實施成通過圖6中的流程圖600舉例說明的數值算法的一部分。在步驟601中，讀取噴嘴的幾何形狀。然後，在步驟602中，讀取控制參數。典型的控制參數是模擬的結束時間；界面模糊的範圍和水平集重新初始化的頻率。給定噴嘴的幾何形狀，在步驟603中創建了如圖2中的經過體擬合的四邊形網格200。在步驟604中，計算了變換矩陣T和雅可比矩陣J。在步驟605中，將時間、當前的時間步長數目和初始的流體速度設為零。給定模糊參數，在步驟606中設定了界面的厚度。可將最初的墨水至空氣界面104假設為平坦的，相應地，在步驟607中初始化水平集。
在步驟608中，通過檢查時間t是否少於循環的設計的運行時間量tend，開始了該循環。如果是這樣，在步驟609中，計算確保代碼的穩定性的時間步長。否則停止模擬。在步驟610中，更新時間。在步驟611中，用時間步長和墨水流速確定墨水的流入速度。一旦確定了該流入速度，則用CFD代碼求解控制墨水行為的偏微分方程。在步驟612中計算了水平集。在步驟613中確定是否在步驟614中對水平集進行重新設距。可以在特定數目的時間步長之後對水平集進行重新設距，或使用其他標準。在步驟615中用水平集的值計算了粘度和密度。在步驟616中，可以計算速度預測值。在步驟617中，將速度的預測值投影到無散度空間來獲得新的壓力和不可壓縮的速度場。在步驟618中，計算了新的墨水流速。在步驟619中，增加時間步長，且循環回到步驟608。
IV.水平集的重新設距如之前指出的，當使用現有技術描述的方法來解決這類問題時，會產生如圖4所示的那類贗像。為解決這個問題，發明人開發了以下對水平集進行重新設距的方法。這些方法可以在步驟614中執行。這些方法也可以在業內熟知的其他算法的操作中執行。
本發明的實施例可包括局部插值法；全局插值法；以及指導在兩種方法之間進行選擇的切換邏輯。可以將全局插值方法實施為具有準確的距離計算的二維輪廓繪製算法。上述局部插值法可以是一種雙三次插值法。在與界面的距離處於一個單元內的各節點處，可使用雙三次插值法或輪廓繪製方法計算兩個新的水平集值。當在這兩種方法之間進行選擇時，切換邏輯可使用輪廓繪製方法的結果。
本發明的另一個實施例可包括雙三次插值法、局部線性插值法以及指導在這兩種方法之間進行選擇的切換邏輯。又一個實施例包括使用三角剖分的快速行進法。可以將三角剖分的快速行進法用於不靠近界面104的單元或與界面104的距離大於一個單元的單元。不靠近界面104的其他單元未被重新設距。
A.雙三次插值法圖7是偏移計算網格700。網格700的節點與圖5A所示的φn(i，j)的位置重合。參考圖7，對於中心單元i，j，建立了局部插值函數.f(r，z)，且該函數被用於對水平集進行重新設距。這種局部插值函數的一個實例是在申請號為10/957349的、由本申請的同一發明人於2004年10月1日提交的美國專利申請中描述的雙三次插值函數，此處通過引用將其包含於本文之中。在中心單元i，j的節點處，使用來自單元i，j和它的8個最近的相鄰單元的信息計算了水平集的值。對某些節點使用雙三次插值法是無效的。例如，如果缺少某些最近的鄰居(如鄰近模擬邊界的單元i，j)，則使用雙三次插值法是無效的。在本申請中使用的雙三次插值應指雙三次插值方法或化簡的雙三次插值方法。
B.局部線性插值法當雙三次插值法不適用時，也可以使用局部線性插值法。圖8A示出了偏移四邊形網格800。如網格700一樣，四邊形網格800偏離網格200，因為網格800的節點集中於計算水平集φ的各點。圖8A也示出了在其中可使用局部插值法的那些點。圖8A示出了這樣的一個實例節點i，j的四個相鄰節點中，僅有一個節點處於界面104的另一側。在這裡，可以將節點處i，j處的新的水平集值視為從節點x(i，j)，y(i，j)至界面104與網格800的網格線相交的點a1，b1的帶符號距離s，如公式(19)至(21)所示，可以用φ(i，j)與φ(i，j-1)計算距離s。或者，可以通過計算從上述節點到界面104的最短距離來找到新的水平集值。上述計算是以四邊形網格為背景描述的。如業內熟知的，也可以在均勻網格中進行這些計算。
圖8B示出了節點i，j，該節點的四個相鄰節點中的兩個處於界面104的另一側。在界面與網格線相交的點中，處於右側的點的坐標為a2，b2，處於下方的點的坐標為a1，b1。通過以與公式(19)和(20)類似的方式對舊的水平集值φ(i，j)，φ(i+1，j)和φ(i，j-1)進行線性插值得到了這些坐標。可以以如公式(22)中所示的參數格式表示從連接a1，b1到a2，b2的分段至節點i，j的最近點的坐標。可使用公式(23)-(25)計算φ(i，j)的新值。
圖8C示出了節點i，j，該節點的相鄰節點中的三個處於界面104的另一側。如上所述，通過對舊的水平集值φ(i，j)，φ(i+1，j)，φ(i，j+1)和φ(i，j-1)進行線性插值可以容易地獲得界面與網格線相交的處於右、下、上側的點的坐標a1，b1，a2，b2，a3，b3。最近的點可能處於連接a1，b1和a2，b2的線上，或處於連接a3，b3與a2，b2的線上。因此，我們重複圖8B中提出的方法，並選擇可能的新水平集值中較小者。在這種情況下，用於求解新的水平集值的公式由公式(26)-(30)給出。
圖8D示出了節點i，j，該節點的兩個相對的相鄰節點處於界面104的另一側。新的φ(i，j)可以是從x(i，j)，y(i，j)到a1，b1的距離或從x(i，j)，y(i，j)到a2，b2的距離中的較短者。或者，新的φ(i，j)可以是從x(i，j)，y(i，j)到界面104的上部或下部的距離中的較短者。界面104可包括多個未連接的界面。
圖8E示出了節點i，j，該節點的所有四個相鄰節點均處於界面104的另一側。再一次地，通過對舊的水平集值φ(i，j)，φ(i+1，j)，φ(i，j+1)，φ(i，j-1)進行線性插值，能容易地獲得界面與網格線相交的處於左、右、上、下側的點的坐標a2，b2，a4，b4，a1，b1，a3，b3。最近的點可能處於連接a1，b1和a2，b2的線上或處於連接a3，b3和a4，b4的線上。因此，我們重複上述過程，並選擇可能的新的水平集值中的較小者。
圖9是關於網格800的另一個圖，在其中，指明了水平集值φ(i，j)和各個水平集值φ(i，j)的笛卡爾坐標x(i，j)，y(i，j)。圖9也示出了四個最近的鄰居和它們的笛卡爾坐標。
圖10的流程圖1000示出了可在本發明的實施例內執行的局部線性插值。如果邊界104在節點i，j與節點i，j的四個最近鄰居中的一個或多個之間通過，則可以對特定的水平集值進行局部線性插值。在步驟1010中，開始FOR循環，這使得一組偏移指數m，n＝{i-1，j；i+1，j；i，j-1；i，j+1}從步驟1030至1060進行循環。在步驟1030中，將φ(i，j)的符號與φ(i+m，j+n)的符號進行比較。對給定的偏移指數m，n而言，如果上述符號相同，則跳過步驟1040和1050，並使用下一個偏移指數。如果上述符號相反，則在步驟1040中，用線性插值來確定邊界與φ(i，j)和φ(i+m，j+n)之間的線的交點。在步驟1050中，將這一點存儲為φtemp(m，n)。在步驟1060中，確定是否檢查了每個偏移指數。如果沒有，則重複步驟1030至1050。在所有的偏移指數經過循環之後，終止FOR循環。在步驟1070中，臨時將由φtemp(m，n)形成的組中絕對值最小者指派給φtemp(i，j)。在步驟1080中，將φnew(i，j)的新值確定為φ(i，j)的符號乘以φtemp(i，j)。可以對每個水平集值φ(i，j)或假設靠近界面104的水平集值的子集進行如流程圖1000中所述的局部線性插值。
步驟1070和1080的一種替代方法包括使用一個或多個近似的界面分段來逼近界面104，這些分段與來自φtemp(m，n)的那組點中的一對或更多對點相交。可找到一條或多條與節點(i，j)相交並垂直於上述近似的界面分段之一的線。可以將水平集φnew(i，j)的新值設為最短垂直線的長度。
上述的計算是在物理空間X中進行的。或者，這些計算也可以在計算空間Ξ中進行。為清楚起見，在流程圖1000中未描述異常處理。為正確執行本發明，有必要實施異常處理。需要得到處理的異常的例子是模擬邊界，且此時界面104與某節點相交而使得水平集值φ為零。此外，可以容易地將上述方法推廣到三維或更多維。
C.全局插值法處於本發明的背景中的全局插值法是指用於對水平集的特定節點進行重新設距的方法。當執行全局插值時，首先尋求零水平集Γ的某種近似。這種近似的一個實例是由沿零水平集Γ的點表示的一系列線段。可通過找到從該特定節點到這一系列線段的最短距離來計算出特定節點的新的水平集值。
1.用輪廓繪製算法找到零水平集Γ可以用輪廓繪製算法找到水平集函數φ的零水平集Γ。對可視化目的而言，零水平集Γ可能值得關注，或者，當對水平集函數φ進行重新設距時，可以使用零水平集Γ。可以將對水平集函數的近似存儲為水平集值φ(i，j)的二維數組。輪廓繪製算法的目標是找到連接單元節點的網格線上的「零點」，即零水平集Γ與網格的網格線的交點。
輪廓繪製算法使用兩個標誌數組lf(i，j)與mf(i，j)來標記已就與界面104的相交進行檢查的網格線。數組lf和mf中的每個元素可以為單個位，而這些數組的維數可以與水平集φ(i，j)的維數相同。標誌數組lf的元素(i，j)表示連接節點(i，j)和節點(i，j-1)的第一網格線。類似的，數組mf的元素(i，j)表示連接節點(i，j)和節點(i-1，j)的第二網格線。可將這些標誌數組初始化成全為零，這表示網格線均未得到檢查。一旦檢查了一條網格線，則將相關標誌數組中的相應元素設為1。可使用變量來存儲已找到的零點的當前數目。
可通過從上至下搜索每條網格線來完成對零點的搜索。或者，可以在域邊界上進行對零點的搜索，然後在內部的網格線上進行這種搜索。一旦找到界面104的第一個零點，便可以跟隨界面104，直到在域邊界或第一個零點處結束為止。重複上述步驟，直到每條網格線均得到檢查為止。將零點存儲到2xn維數組z(xm，ym)，其中m是從1至n+1的指數。
上述的跟隨界面104包括檢查特定零點的北方、南方、東方和西方的相鄰網格線，以確定界面104是否與這些網格線相交。跟隨界面104也包括檢查特定的網格線，以確定界面104與該網格線的交點是否多於一個。跟隨界面104也包括確保界面104不與其自身相交。
上述方法僅僅是使用輪廓繪製算法找到界面104的其中一種方法。另一種用於找到界面104的方法是計算水平集函數的一階和/或二階梯度，並使用該信息來找到和跟隨界面104。業內存在一些熟知的用於從水平集φ找到界面104的方法。
2.從輪廓圖計算新的水平集圖11的流程圖1100示出了一種在給定表示界面104的一組零點z(xm，ym)的情況下對水平集φ的特定節點(x0，y0)進行重新設距的方法。給出這種方法僅僅是為了說明的目的，其他當界面104的輪廓可得到時計算特定節點處的水平集φ的方法在業內是熟知的。理想地，可將這種方法用於與界面104的距離處於一個單元之內的水平集值。圖12示出了包括界面104的一部分的四邊形網格。圖12也示出了流程圖1100中引用的某些變量。
如上所述，數組z是逼近界面104的一系列點。一般來說，界面104由一條或多條連續線組成。這些連續線在模擬空間的邊界處開始和結束。或者，這些連續線中的一條或多條線可形成封閉的表面(如包圍墨滴的表面)。可以對數組z中的上述一系列點進行排序，以使得形成界面104的各分段的每對點在數組中彼此相鄰。如果有必要，可以添加複製的點。可以將其他值添加到數組中，以指明界面104的終點或其他異常情況。
步驟1110初始化FOR循環，這使得當沿界面104檢查所有分段時使用指數m。根據公式(31)將開始於點xm，ym並繼續達到點xm+1，ym+1的分段進行參數化。在步驟1120中，根據公式(31)對給定的節點和給定的分段m求解參數t。在步驟1130中，如果t大於1，將t設置為1。在步驟1140中，如果t小於0，將t設置為0。執行步驟1130和1140，以將t限於分段m上的點。圖12示出了分段m和m+1的點t。在步驟1150中，根據公式(31)計算沿分段m的參數t的笛卡爾坐標(xt，yt)。在步驟1160中，根據公式(31)計算從節點(x0，y0)至(xt，yt)的距離sm。圖12示出了對應於線段m和m+1的距離sm和sm+1。在步驟1170中，存儲值sm，以供後面使用。在步驟1180中，終止FOR循環，或增加指數m並重複步驟1120至1170。FOR循環的輸出可以是由距離s形成的數組。如果終止FOR循環，則在步驟1190中，從步驟1170中存儲的由距離s形成的數組中的最小值計算臨時的水平集值。在步驟1195中，通過將節點(x，y)處的舊水平集值的符號乘以該臨時的水平集值計算了節點(x，y)處的新的水平集值。步驟1195的一種替代方案是與舊的水平集值無關地計算新的水平集值的符號。步驟1190和1170的一種替代方案包括僅存儲小於計算的直至那一點的最小距離sm。
從而，基於對界面104的全局求算計算了特定節點的新的水平集值。本發明的另一個實施例包括對界面104的全局求算，然而，未將整個界面104用於計算新的水平集值。相反，僅使用了靠近節點(x0，y0)的單元中的有限的一組界面分段來計算新的水平集值。
為簡單起見，未在流程圖1100中示出異常處理。本發明的實施例可包括異常處理，如邊界條件和其他異常情況。已經以物理空間X中的二維四邊形網格為背景描述了全局插值。經過適當調整，上述方法也可以在計算空間Ξ中的矩形網格上執行。
D.切換邏輯對水平集φ進行重新設距包括使用多種插值方法，這些方法包括雙三次插值法、全局插值法、局部線性插值法。在對水平集進行重新設距時，本發明的實施例可使用這些方法中的兩種或三種方法。在本發明中，將在對水平集的特定節點進行重新設距時選擇使用何種插值方法稱為切換邏輯。使用切換邏輯的動機包括確保代碼穩定性、減少贗像和提高速度。
雙三次插值法可以是預設的插值方法。在某些情況下也可以使用三角剖分的快速行進方法。如果某特定節點是鄰近域邊界的單元的一部分，則可以對該特定節點使用全局插值法或局部線性插值法。如果某特定節點處或其鄰近區域中的流體加速度高於閾值，則也可以對該節點使用全局插值法或局部插值法。另外，如果某特定節點是包括具有銳角轉角的界面104的單元的一部分，則可以對該特定節點使用全局插值法或局部插值法。而且，如果某特定節點是包括具有大於閾值的曲率的界面104的單元的一部分，則可以對該特定節點使用全局插值法或局部插值法。一般來說，對大多數噴墨模擬情形而言，僅需為有限的一組節點而從雙三次插值起簡化插值。因此，在界面附近，對水平集進行重新設距的準確度仍為三階。
1.鄰近域邊界如果界面104通過鄰近域邊界的單元，則不應使用雙三次插值法。為建立雙三次插值基於的雙三次多項式，雙三次插值需要16個條件。對鄰近域邊界的單元而言，不存在16個條件。從而不應使用雙三次插值法。另外，使用雙三次插值法對鄰近邊界的單元處的水平集進行重新設距具有將界面104推離域邊界的趨勢。這將導致在模擬中引入贗像。例如，在緊接對稱軸的位置使用雙三次插值法常常導致模擬產生一個非常長的未夾斷(pinch off)的墨滴。這樣的結果是模擬中的贗像，並違背了毛細不穩定性的物理原理。
2.高的局部流體加速度本申請中描述的水平集投影方法在時間上是明確的。時間步長受Courant-Friedrichs-Lewy(CFL)條件、表面張力、粘度和局部流體加速度的約束。如公式(32)所表示h是最小單元的尺寸(對於軸對稱坐標，h＝minimum(Δr，Δz))；F在公式15中定義；Re是雷諾數；We是韋伯數；ρ是相對密度；ρ1和ρ2是第一種流體和第二種流體的密度；μ是相對動力學粘度；u是徑向速度；v是垂直速度。
可以在模擬開始時確定來自表面張力和粘度的約束，因為它們僅取決於網格大小和流體性質。應該在每個時間步長處對每個單元檢查CFL和流體加速度。為了實施的目的，如公式33所示，可根據表面張力和粘度約束計算初始的時間步長Δt0。可以在每個時間步長處計算公式34中的Δt』。用這兩種約束中較小者作為實際的時間步長大小，即Δt＝minimum(Δt0，Δt』)。對於給定的單元，如果Δt』＜Δt0，那麼應當用全局插值法或局部線性插值法來對該單元的水平集值進行重新設距。
3.銳角轉角當界面104在單元中或接近單元處具有銳角轉角時，不應使用雙三次插值法，因為它具有平滑銳角轉角的傾向。對於在一個方向上連續彎曲的界面104的部分，雙三次插值法能很好的工作。但是，當界面104開始在新的方向上彎曲時，雙三次插值法不能很好的工作。至少因為這個原因，如果兩個分段之間的角小於特定的閾值，不足以宣稱不應使用雙三次插值法。必須使用更為精細的方法。
用於執行全局插值的界面104的近似可用來分辨界面是否具有銳角轉角。圖13示出了網格800上的界面104的一部分。界面104在分段1和分段2之間具有銳角轉角。從而，不應將雙三次插值法用於包含分段1和2的單元，而可以將其用於包含分段3，4，5，6和7的單元。
假設我們考慮如圖13所示的具有端點xm，ym和xm+1，ym+1的界面104的第m分段，其中，分段1是第m分段。各分段可以由公式35中所示的矢量(Δxm，Δym)表示。應計算由第m分段與前一分段以及下一分段形成的夾角αm和αm+1。這可以使用如公式36所示的餘弦定理來完成。圖13示出了若干個可以通過這種方式計算的夾角(αm-1…αm+4)。如圖13所示，夾角αm是由第m分段和之前分段以及連接這兩個分段的線形成的三角形的內角。其他用於計算或逼近由兩條線形成的夾角的方法在業內是熟知的。
用於確定單元是否包含銳角轉角的另一信息是分段彎曲的相對方向。可通過計算上述之前分段與第m分段的叉積得出該信息，即dirm＝(Δxm-1，Δym-1)×(Δxm，Δym)。使用的來自該計算的唯一信息是dirm的符號。該叉積的結果是圖13中所示的分段m-1至分段m+4。指示離開頁面的方向的矢量用點表示，而指示進入頁面的方向的矢量用叉表示。
如果包含分段m的單元滿足以下兩個條件之一，則不應使用雙三次插值法，而應使用全局插值法或局部線性插值法dirm·dirm+1≥0且0＜αm≤40°且0＜αm+1≤40°；或dirm·dirm+1＜0且(0＜α≤ 75°或0＜αm+1≤75°對於第一個條件，第一上限可以是40°，儘管也可以使用諸如30°，50°或60°的其他角度，而不至於超出本發明的範圍。對於第二個條件，第二上限可以是75°，儘管也可以使用諸如60°或80°的其他角度，而不至於超出本發明的範圍。在本發明的實施例中，第一上限小於第二上限。
4.界面的曲率當界面104在特定的單元內或其附近具有銳角轉角時，求算的該特定單元的節點處的水平集的曲率κ將傾向於比較高。在四邊形網格上，可通過公式37計算曲率。通過使用中心差分公式，可以在計算空間中以數值計算的方式計算上述公式中的所有導數。計算導數的其他方法在業內是熟知的。
如果單元中的任意節點的曲率高於閾值，則認為該單元的曲率是高的。該閾值與上面提到的表示界面104的模糊程度的參數ε成反比。如果曲率大於該閾值，則不應使用雙三次插值法而應使用全局或局部線性插值法。
本說明書已描述了本發明的細節，現在，將結合圖14描述用於實施本發明的一個或多個方面的示範性系統1400。如圖14所示，該系統包括提供計算資源和控制計算機的中央處理單元(CPU)1401。CPU1401可以用微處理器等實現，也可以包含圖形處理器和/或用於數學計算的浮點協處理器。系統1400也可包括隨機存取存儲器(RAM)和只讀存儲器(ROM)形式的系統存儲器1402，。
如圖14所示，也可以提供若干個控制器和外圍設備。輸入控制器1403表示至鍵盤、滑鼠或指示筆等各種輸入設備1404的接口。也可設置與掃描儀1406通信的掃描儀控制器1405。系統1400也可包括與一個或多個存儲設備1408連接的存儲控制器1407，這些存儲設備1408中的每個設備均包括磁帶、磁碟或光學媒體之類的存儲介質，這些存儲介質可用於記錄作業系統的指令程序、包括實施本發明的各個方面的程序實例在內的公用程序與應用程式。根據本發明，也可以用存儲設備1408存儲處理後的數據或待處理的數據。系統1400也可包括用於向顯示設備1411提供接口的顯示控制器1409，所述顯示設備1411可以是陰極射線管(CRT)或薄膜電晶體(TFT)顯示器。系統1400也可包括與印表機1413通信的印表機控制器1412。通信控制器1414與一個或多個通信設備1415進行連接，這些通信設備使得系統1400能通過包括網際網路、區域網(LAN)、廣域網(WAN)在內的各種網絡或包括紅外信號在內的任何合適的電磁載波信號與遠程設備相連。
在所示的系統中，所有主要的系統部件均連接到可表示多於一條物理總線的總線1416。然而，各種系統部件彼此之間可以相鄰或不相鄰。例如，可以將輸入數據和/或輸出數據從一個物理位置遠程地發送到另一個物理位置。同樣，可通過網絡從遠程位置(如伺服器)訪問實施本發明的各個方面的程序。這些數據和/或程序可以通過各種計算機可讀媒體進行傳送，所述媒體包括磁帶、磁碟、光碟、網絡信號或包括紅外信號在內的任何合適的電磁載波信號。
可以方便的用軟體實現本發明。然而，也可以存在其他實施方式，包括硬體實施方式或軟體/硬體實施方式。可使用ASIC、數位訊號處理電路等等來實現任何用硬體實現的功能。相應的，權利要求中的「裝置」既包括軟體實施方式，又包括硬體實施方式。類似的，本文中使用的術語「計算機可讀媒體」包含軟體、其上包含硬接線形式的指令程序的硬體以及所述軟體和硬體的組合。考慮到存在這些替代的實施方式，應當懂得，附圖和說明書提供了本領域技術人員編寫程序代碼(即軟體)或製造電路(即硬體)來執行要求的處理所需的功能信息。
根據本發明的其他方面，可以將任何上述方法或步驟體現為存儲在、或傳送到計算機或其他處理器控制的設備以供執行的指令程序(即軟體)。或者，可以用功能上等同的硬體(如專用集成電路(ASIC)、數位訊號處理電路等)或軟體與硬體的組合來實施任何所述的方法或步驟。
儘管結合若干個具體的實施例描述了本發明，但是，對本領域的技術人員而言，根據上述說明，其他替代形式、修改、變更和應用是顯而易見的。其中，這些替代形式的例子包括更高維數；其他坐標系；其他流體；其他邊界條件；其他通道；其他控制公式以及具有多於兩種流體的系統。因此，本文描述的發明意在包括所有落在附錄的權利要求的精神和範圍內的替代形式、修改、變更和應用。
附錄·u＝0，DuDt=-1p+1Re(2D)-1We-1Frez,...(1)]]>t+u=0.]]>J=gdet=gdetrrzz....(2)]]>T=g-1J[]-1=z-T-zr....(3)]]>Ξ·u＝0，u＝qTu， -gJ2We(gTTT|TT|)(TT),...(4)]]>t+J-1u=0.]]> =gJRe[TT][gJ-1TTu+(gJ-1TTu)T]...(5)]]>+JRe{g2J-1TTTu}+Re-ur20.]]> =21+(1-21)H....(7)]]>n+1=n-tJ[u]n+1/2....(8)]]>
(gTu*)=(g2t(n+1/2)JTTTpn+1)....(10)]]>un+1=u*-gt(n+1/2)JTTpn+1....(11)]]>[(u)u]i,jn+1/2=ui+1/2,jn+1/2+ui-1/2,jn+1/22(ui+1/2,jn+1/2-ui-1/2,jn+1/2)]]>+i,j+1/2n+1/2+i,j-1/2n+1/22(ui,j+1/2n+1/2-ui,j-1/2n+1/2)....(12)]]>[(u)]i,jn+1/2=ui+1/2,jn+1/2+ui-1/2,jn+1/22(i+1/2,jn+1/2-i-1/2,jn+1/2)]]>+i,j+1/2n+1/2+i,j-1/2n+1/22(i,j+1/2n+1/2-i,j-1/2n+1/2)....(13)]]>ui+1/2,jn+1/2,L=ui,jn+12u,i,jn+t2ut,i,jn]]>=ui,jn+(12-t2Ji,jui,jn)u,i,jn-t2Ji,j(u)i,jn+t2Fi,jn....(14)]]> ui+1/2,jn+1/2,R=ui+1,jn-12u,i+1,jn+t2ut,i+1,jn]]>=ui+1,jn-(12+t2Ji+1,jui+1,jn)u,i+1,jn-t2Ji+1,j(u)i+1,jn+t2Fi+1,jn....(16)]]>
i,j-1+sy(i,j-i,j-1)=0,y=(xi,j-xi,j-1)2+(yi,j-yi,j-1)2...(19)]]>s=i,j-1i,j-1-i,jy....(20)]]>(φi，j)新＝sign(φi，j)舊×s. (21)(a1+tΔa，b1+tΔb)，(Δa，Δb)＝(a2-a1，b2-b1).(22)((x，y)i，j-(a1+tΔa，b1+tΔb))·(Δa，Δb)＝0. (23)t=a(xi,j-a1)+b(yi,j-b1)a2+b2....(24)]]> t1=(a2-a1)(xi,j-a1)+(b2-b1)(yi,j-b1)(a2-a1)2+(b2-b1)2....(26)]]>|i,j1|=[xi,j-a1-t1(a2-a1)]2+[yi,j-b1-t(b2-b1)]2....(27)]]>t2=(a2-a3)(xi,j-a3)+(b2-b3)(yi,j-b3)(a2-a3)2+(b2-b3)2....(28)]]>|i,j2|=[xi,j-a3-t2(a2-a3)]2+[yi,j-b3-t(b2-b3)]2....(29)]]> xt＝xm+t(xm+1-xm)yt＝ym+t(ym+1-ym)，t=(xm+1-xm)(x0-xm)+(ym+1-ym)(y0-ym)(xm+1-xm)2+(ym+1-ym)2,...(31)]]>sm=(x0-xt)2+(y0-yt)2.]]>
tmini,j[ru,z,We1+28h3/2,Re2nn(1r2+1z2)-1,2h|F|]....(32)]]>t0mini,j[We1+28h3/2,Re2nn(1r2+1z2)-1]....(33)]]>tmini,j[ru,z,2h|F|]....(34)]]>(Δxm，Δym)＝(xm+1，ym+1)-(xm，ym).(35)m=arccos(xm-1xm+ym-1ym)xm-12+ym-12xm2+ym2,]]>m+1=arccos(xm+1xm+ym+1ym)xm+12+ym+12xm2+ym2....(36)]]>=J-1(gTTT|TT|)....(37)]]>
權利要求
1.一種用於模擬經過通道和離開通道的流體流的方法，該流體流包含第一流體和第二流體之間的界面，所述方法包含以下步驟創建表示所述通道的物理空間和所述通道周圍的一部分物理空間的網格；創建包括一組值的水平集，每個值與所述網格中的一個點相關，這組值中的每一個值與從所述相關點至所述界面的最短距離成比例；求解描述所述第一流體、第二流體以及所述界面的各個方面的一組公式；使用以下重新設距方法中的兩種或更多種方法對所述水平集中的特定值進行重新設距雙三次插值法、全局插值法或局部線性插值法；以及對每個特定值而言，根據一個或多個切換規則在雙三次插值法與全局插值法或局域線性插值法之間進行切換。
2.如權利要求1所述的方法，其中，所述第一流體是墨水，所述第二流體是空氣，且所述通道包括作為噴墨頭的一部分的噴墨噴嘴。
3.如權利要求1所述的方法，其中，所述網格是四邊形網格，且所述方法還包括以下步驟將所述網格變換成計算空間中的均勻的正方形網格；將所述的一組公式從所述物理空間變換到所述計算空間；在所述計算空間中求解該組公式。
4.如權利要求1所述的方法，其中，求解這組公式包括使用有限差分方法。
5.如權利要求1所述的方法，其中，所述切換規則包括對一個特定的點而言，如果該特定的點是包括所述界面的單元的一部分，且所述單元的所述界面的一部分具有銳角轉角，則切換到所述全局插值法或所述局部插值法。
6.如權利要求5所述的方法，其中，如果在所述界面的某一特定轉角處連接的兩個連續分段之間的角小於四十五度，則該特定轉角是銳角轉角。
7.如權利要求5所述的方法，其中，如果三個連續分段朝彼此彎曲，且這三個分段之間的兩個角均小於或等於四十度，則所述界面的某一特定轉角是銳角轉角。
8.如權利要求5所述的方法，其中，如果三個連續分段朝彼此遠離的方向彎曲，且這三個分段之間的兩個角之一小於或等於七十五度，則所述界面的某一特定轉角是銳角轉角。
9.如權利要求5所述的方法，其中，如果由包括所述界面的某一特定轉角的或靠近該轉角的單元的四個節點組成的組的水平集的曲率大於臨界曲率值，則所述界面的所述轉角是銳角轉角。
10.如權利要求9所述的方法，其中，所述臨界曲率值與界面模糊的範圍成反比。
11.如權利要求1所述的方法，其中，所述切換規則包括對某一特定的點而言，如果該特定點是與域的邊界相鄰的單元的一部分，則切換到一般的插值法。
12.如權利要求1所述的方法，其中，所述切換規則包括對某一特定的點而言，如果該特定點處的局部流體加速度大於閾值，則切換到一般的插值法。
13.如權利要求1所述的方法，其中，點與所述網格上的節點進行搭配。
14.如權利要求1所述的方法，還包括求解所述公式組來描述一部分所述第一流體從所述通道的噴出的步驟。
15.如權利要求1所述的方法，其中，如果與所述水平集中的某一特定值相關的點處的流體是所述第一流體，則該值具有第一符號，如果該相關點處的流體是所述第二流體，則該值具有相反的符號。
16.如權利要求1所述的方法，其中，在3維坐標系中對所述流體流進行模擬。
17.如權利要求16所述的方法，其中，所述坐標系是軸對稱坐標系，且未模擬沿方位角的流體流。
18.如權利要求1所述的方法，還包括使用三角剖分的快速行進法將所述水平集重新設距於距所述界面大於一個單元的各個點處。
19.一種包括用於執行如權利要求1所述的方法的模塊的裝置。
20.一種包括一組用於指示裝置來執行如權利要求1所述的方法的指令的計算機可讀媒體。
全文摘要
用於模擬經過通道和離開通道的流體流的系統和方法。該流體流包括第一流體和第二流體之間的界面。創建表示上述通道的物理空間和該通道周圍的一部分物理空間的網格。創建表示上述界面的水平集。求解了一組描述第一流體、第二流體和上述界面的各個方面的方程。使用以下重新設距方法中的兩種或更多種方法對水平集中的特定值進行了重新設距雙三次插值法、全局插值法或局部線性插值法。對每個特定值而言，根據一個或多個切換規則在上述重新設距方法之間進行切換。
文檔編號G06F17/50GK101049769SQ20071009682
公開日2007年10月10日 申請日期2007年4月4日 優先權日2006年4月6日
發明者於峻德 申請人:精工愛普生株式會社