模擬流體中的粒子運動時，選擇合適的公式以提升計算效率

學時習

2023年9月27日 10:39

瀏覽：3658 收藏：10

當你第一次嘗試使用 COMSOL 軟件的粒子追蹤模塊模擬流體中非常小的粒子（通常直徑為幾十微米或更小的粒子）時，可能會發現瞬態求解器使用的時步比平常要短得多。這通常是由于粒子的運動方程表現出數值剛度而導致的。在這篇文章中，我們將介紹與粒子仿真有關的剛度的概念，并提供一些基于粒子大小選擇正確方程的指南。

示例：小球形粒子的重力沉降

以一個小的球形粒子為例，當它掉落在一個以速度 u （SI 單位：m/s）勻速流動的流體中時，遵循牛頓第二運動定律

(1)

其中，

m_p（SI 單位：kg）是粒子的質量
q（SI 單位：m）是粒子的位置矢量
F_t（SI 單位：N）是作用在粒子上的凈力或總力

對于一個在流體中下沉的粒子，總力是重力 F_g 和曳力 F_D 的總和，

（2）

其中，

ρ_p（SI 單位：kg/m^3）是粒子的密度
ρ（SI 單位：kg/m^3）是周圍流體的密度
g（SI單位：m/s^2）是重力引起的加速度（海平面以下約為 9.8m/s^2）
μ（SI 單位：Pa s）是周圍流體的動力黏度
d_p（SI 單位：m）是粒徑
u（SI 單位：m/s）是周圍流體的速度
v（SI 單位：m/s）為粒子的速度（v≡dq/dt）

重力表達式中的（ρ_p-ρ）/ρ_p 項代表浮力。當粒子（例如，空氣中的固體粒子）比其置換的流體重得多時，浮力的值接近 1。當粒子和周圍的流體具有相同的密度時，其值接近于零，在這種情況下，粒子被稱為懸浮粒子。

此處使用的曳力表達式來自斯托克斯曳力定律（Stokes drag law）。當粒子的相對雷諾數非常小時，此曳力定律較為適用：

這對于較小的粒子更有效。

假設粒子沒有改變大小（d_p 和 m_p 為常數），則球形粒子的質量是

(3)

結合方程1–3，我們得到粒子運動方程的簡化表達式：

(4)

這里，引入了常數 τ_p，

τ_p具有時間單位，并且通常被稱為拉格朗日時間尺度（Lagrangian time scale）或粒子速度響應時間（particle velocity response time），下面我們將解釋其原因。

進一步簡化方程，假定周圍的流體是靜止的（?=0），并且所述粒子初始是靜止的（q=0, v=0, 時間 t= 0）。假設我們對齊坐標系，以使重力矢量指向 –y 方向。然后，根據方程4，粒子位置的 y 分量方程變為

(5)

當給定初始條件為 q_y=0 和 v_y=0時，方程5 的精確解或解析解為

其中，v_t 是自由沉降速度，

轉換為無量綱變量

為了更好的理解粒子在 τ_p 之前最初一段時間是如何加速的，我們可以用相應的無量綱量(t?, q_y?, v_y?)代替時間、位置和速度（t, q_y, v_y），定義為

將這些無量綱變量代入解析解中，得到

在下圖中，無量綱的位置和速度被繪制成為無量綱時間t‘的函數。該繪圖表明粒子速度逐漸接近自由沉降速度，粒子加速主要發生在拉格朗日時間尺度τp最初一段時間。在此初始加速期之后，粒子位置似乎呈線性變化。

模擬流體中的粒子運動時，選擇合適的公式以提升計算效率的圖11 從靜止開始，經歷重力沉降的粒子的無量綱位置和速度的繪圖。

一些典型粒徑的時間尺度

為了更好地了解粒子加速所涉及的時間尺度，假設粒子為密度約為 2200 kg/m^3 的石英玻璃珠。下表列出了不同粒徑的粒子在空氣和水中的一些拉格朗日時間尺度值。

τ_p 和直徑平方呈線性關系意味著大粒子比小粒子具有更長的速度響應時間和更大的自由沉降速度。這會產生兩個主要結果：

大粒子掉落到地面的速度比小粒子快得多。
當大粒子以一定的初始速度射入流體時，會沿著入射軌跡，能夠在曳力使它們減速之前行進相當長的距離。相反，較小的粒子將更快地匹配流體速度。當它們散開時，很可能是由于周圍流體的湍流擴散所致。

數值粒子追蹤仿真

在上一節中，很幸運我們由方程4 得到一個精確的解析解。精確解僅可能在引入許多簡化假設時得到，尤其是各處的流體速度 u 均為零。但在大多數實際情況中，周圍流體的速度不僅不為零，而且在空間上是不均勻的，因此僅靠公式不可能找到精確解。

對于更一般的問題，我們可以通過數值仿真來獲得近似解。其主要思想是，在初始時間 t=0 時，給定初始粒子位置 q_0 和速度 v_0，我們可以使用數值時間長算法來計算一組離散的時步 t_1，t_2，t_3，……的解。為此，設計了各種各樣不同的時間步長算法，其中有許多是在 COMSOL Multiphysics? 軟件中可用的。

用數值方法求解一組微分方程會引入一定量的誤差,即實際粒子運動與計算得到的數值解之間的差異。雖然通常不能指望從數值仿真中獲得一個完美的解，但更現實的目標是，當時間間隔（t_1，t_2–t_1，t_3–t_2等）減小時，模擬的粒子運動應變得更加精確。

需要權衡的是，如果時間步較小，則需要花更多的時間步才能達到相同的輸出時間。最終，這可能會導致實際運行時間顯著增加，這是仿真完成的時間。進行數值仿真的工程師必須始終在解精度和執行時間之間尋求合理的平衡。

COMSOL Multiphysics?中的粒子追蹤模塊提供了一個流體流動顆粒跟蹤接口，該接口通過數值求解牛頓第二定律來模擬周圍流體中單個粒子的運動。基本上，此接口可求解方程1，同時允許我們向方程右側添加各種不同的力。它還包括用于設置初始粒子位置和速度以及檢測和處理粒子與周圍幾何中的表面的碰撞的各種選項。

處理小粒子和長時間尺度

在許多實際應用中，粒子追蹤模型的需要求解時間的范圍遠大于拉格朗日時間尺度 τ_p。例如，假設我們要在 1s 的總仿真時間內追蹤水中直徑約 20μm 的石英玻璃顆粒的運動。從上述表格我們可知，水中這樣的小粒子的拉格朗日響應時間約為 5×10^-5 s，所以總仿真時間約為 2000τ_p。如果我們想在幾分鐘或幾小時的跨度內追蹤更小的粒子，那么我們的總仿真時間可能比 τp 大幾百萬倍。

下面的截圖顯示了瞬態求解器在跟蹤這些 20μm 粒子時所采取的時間步日志。在步驟1 中輸出時間的范圍：瞬態節點已被設置為 range(0,0.1,1)，這意味著它將僅以 0.1s 的倍數存儲輸出。但是，這并不妨礙求解器在必要時采取更短的時間步來獲得精確的解。如這里所顯示的，求解器先從采取 1ms 或更小的時間步開始，然后在粒子接近其最終速度時逐漸采取更大的時間步。

如下面的步驟24 所示，在 COMSOL Multiphysics 中，粒子追蹤物理場接口通常使用嚴格的時間步算法，該算法至少要求求解器所采取的某些步長與輸出時間一致。但這并不是所有物理場的普遍要求；對于某些物理場接口，可以通過在求解器所采用的最近步長之間進行插值來獲得輸出時間。

在研究接近尾聲時因為粒子幾乎不再加速，所以時間步可能會很大。最終，求解器需要 24 個時間步才能在 0.1s 達到第一個輸出時間，但是只需要再增加 12 個時間步就能在 1s 到達最終時間。

經歷重力沉降的粒子運動方程是剛性常微分方程（或剛性 ODE）的一個示例。大多數粒子追蹤模型中使用的默認時間步進方法被稱為廣義 α，這是一種二階隱式時間步方案，非常適合用于處理剛性問題。如果需要額外穩定性，則可以在瞬態求解器設置中調整一個被稱為放大高頻的數值阻尼項。因此，隨著粒子速度接近自由沉降速度，時間步變得更大。（相比之下，顯式 Runge–Kutta 方法 RK34 采取 7425 個步長來求解相同的問題！）

但是，如果粒子在幾個不同的釋放時間進入仿真域，或者背景流體速度在空間上不均勻（這樣粒子在以后的研究中仍會加速），則求解器可能直到最終時間都會一直采用如此小的時間步。因為求解器可能需要成千上萬甚至數百萬的時間步，所以如果我們試圖在很長的仿真時間內追蹤非常小的粒子，則最終這些研究將需要大量的執行時間才能完成。

有一個與此密切相關的現象，可能會使 COMSOL? 軟件新用戶感到困惑，當使用入口邊界條件將粒子釋放到仿真域中，并且假設這些粒子被分配了指向仿真域的初始速度。請注意，從之前的截圖中可以看出，初始時步大小（總仿真時間為1秒）為 1 毫秒。如果初始時步仍然遠遠大于 τp，則曳力可能會過度補償，導致粒子速度短暫改變方向并指向入口邊界。如果發生這種情況，粒子可能會錯誤地檢測到與入口邊界的碰撞，從而使它們卡在此處。

粒子追蹤模型中的數值剛度處理

有兩種主要方法可以求解流體中的粒子運動的數值剛度模型，即輸出時間間隔比 τ_p 大幾個數量級的模型。

第一種是我們所說的“強力”方法：只需告訴求解器采取更小的時步即可。如果不想產生大量的輸出（可能會創建大量文件大小），那么可以不考慮輸出時間，而是在求解器序列的瞬態求解器設置中指定一個較小的步長或最大步長。