Altair Taiwan Blog

SPH是Smooth Particle Hydrodynamics光滑粒子流體動力學的簡稱，它是基於插值理論的無網格數值方法。

它是用許多與自己周邊粒子相互作用的粒子（質點）來描述物體（連續體），遵從連續介質力學的力學定律。

將物體的材料特性、速度和品質以一定的權重分配給每個粒子。在計算過程中，連續性的改變是由於物體的變形引起的。

SPH是將連續體動力學的守恆定律，以偏微分方程的形式，通過核函數W(r,h)將其轉化為積分方程。

儘管粒子通常顯示為離散點或球體，但在數學上，它們是連續的。

圖 1 HV中粒子顯示和數學上認為的粒子

核函數W(r,h)是用於描述粒子與其周圍半徑r範圍內的粒子的影響度。

h是光滑長度，它有用戶定義，那麼粒子與其周圍半徑為2h內的粒子產生相互作用的影響。h定義的越大，那麼粒子就與越多的周邊粒子有相互作用，計算精度越高，但是計算資源就越多消耗。

所以需要優化粒子分佈和定義適當的h，在精度和性能之間達到最佳折衷。

圖 2 使用於SPH的核函數

SPH方法實際上也是一種Lagrangian方法，與傳統的Lagrangian方法相比較，使用SPH在計算時不會出現由於單元畸變而引起的計算中沒有沙漏能的問題、沒有負體積問題、出現更少的時間步長問題。

圖 3 球穿擊板的SPH模擬方法和Lagrangian模擬方法

SPH特別適用於模擬具有非常大變形的現象，因為此種情況使用Lagrangian網格，單元的變形是非常大的，進而計算時間步長會變得非常小，計算變得非常沒有效率，此時就需要換用SPH的方法。

鳥撞，液體晃動，波浪工程，彈道學，氣體流動等都可以用SPH方法。

Radioss中的SPH方法與大多數函數相容。例如，可以類比兩個物體相互作用，一個由Lagrangian有限元離散，另一個由SPH粒子離散。

使用Hypermesh可以輕鬆創建SPH粒子。可以在1D -> sph進入，也可以從Mesh-> Create-> SPH進入，基於surface和volume創建SPH粒子。

圖 4 HyperMesh中的SPH網格劃分功能塊

Radioss中的粒子可以有三種分佈方式：

 · 一種是六邊形密網分佈法 HCP，

 · 一種是面心立方堆積Face-Centered Cubic FCC

 · 以及還有一種是立方體頂點分佈法。

由於HCP和FCC的分佈方式的粒子最為緊湊，所以推薦使用這兩種分佈方式，而且這兩種分佈方式得到的結果也是非常接近的。

在HyperMesh中(如上圖面板中)可以直接得到FCC的分佈方式。而HCP的分佈方式需要從Altair Connect上下載一個HyperMesh的腳本，然後可以生成相應的HCP粒子分佈。

圖 5 SPH粒子的HCP/FCC分佈方法

例如下面就是在HyperMesh中設置間距為1mm的FCC分佈的粒子：

圖 6 HyperMesh中設置SPH粒子的面板

建議一個模型中選用一種類型的粒子分佈方式，並且所有粒子的間距（分佈密度）要一樣，不同尺寸的粒子之間的相互作用可能導致異常，有些情況下接觸會發生問題。

粒子的單元屬性使用/PROP/TYPE34(SPH)。

如果在HyperMesh中的粒子設置面板（如圖 6 HyperMesh中設置SPH粒子的面板）設置了物體的密度”material density”，那麼HyperMesh會根據產生的粒子數目自動計算每個粒子的品質mp。

它是這樣計算的：

（這裡n是粒子數目）

並且自動計算初始的光滑長度h0，也就粒子與其最近的粒子的距離。這個h0是這樣計算的：

如果是FCC分佈的粒子那麼h取值稍微大一些的√2h0，最大可以取2h0。

所有下面的卡片在HyperMesh中定義好粒子間距和構件密度後自動產生的：

圖 7 SPH粒子的單元屬性卡片

在計算過程中實際的光滑長度h值會隨著材料變形引起的密度變化而變化（遵從 p · h = 常量）。使用/H3D/SPH/ DIAMETER輸出計算過程中每個時刻的光滑長度h值。

在粒子單元屬性卡片中的alpha_cs這個參數是用於保證SPH粒子構件受拉時的穩定性。

alpha_cs的預設值是基於FCC粒子分佈時自動設置的，建議受拉時設置alpha_cs=0.1 （也就是10%），但是這樣可能會導致能量的吸收，建議在在HyperGraph中檢查HE能量。

另外當流體受到負壓時可以選用參數xi_stab來控制來穩定計算，這個參數使用時一般建議設置為0.3。

在粒子單元屬性卡片中還有qa,qb(Artificial viscosity)人工粘度的參數，它是用於正確類比衝擊和高壓縮區域粒子之間相互作用。一般默認使用：

qa = 2

qb = 1

這是用於使用於類比衝擊波/波傳播情況，比如衝擊速度接近聲速在材料中傳播的速度時。

在低馬赫時，由於內部計算人工粘度時引入過度阻尼，所以係數必須降低，qa，qb應如下設置為不考慮阻尼的情況：

qa = 2.0e-30

qb = 1.0e-30

比如鳥撞中的SPH鳥模型需要設置這個參數。

這裡要注意真實的物理粘度必須通過材料卡片設置，而不是在這裡設置。

在Radioss中設置粒子的單元屬性時還用一個用於定義全域的單元屬性卡片/SPHGLO。

這裡的Maxsph在v14.0.22版本後就不需要設置了，Radioss自動計算。

Nneigh是定義周邊粒子的最大個數，當在定義的2h半徑範圍內搜索到超過Nneigh的粒子個數，那麼就在這些粒子中選取最靠近的Nneigh個粒子當中的Lneigh個例子計算相互作用，所以Nneigh ≥ Lneigh。通常使用預設值就可以了。

Alpha_sort是用於計算搜索相鄰粒子的安全半徑，也就是每個粒子會搜索到大於安全半徑內的相應粒子數，也就是擴大了（1+Alpha_sort）的搜索半徑，默認的Alpha_sort=0.25。這個參數最好不要隨意改動，增加這個參數，由於有更多的相鄰粒子包含進去，會大大降低計算速度。

在Radioss中SPH粒子法類比的物體可以使用絕大多的可以用於3D實體單元的材料卡片，具體可以參見工具書關於材料卡片相容性的章節。

通常我們可以使用LAW4（HYD_JCOOK）來模擬金屬材料（SPH粒子），常用於穿擊金屬板，切削金屬等類比；可以使用LAW6（HYD_VISC）來模擬流體，比如水，油等。

這裡注意粘滯性是在材料卡片這裡設置，比如在LAW6中設置，而不是在單元屬性卡片中的qa,qb中設置。

那麼各項異性的材料也可以用SPH粒子法，比如加勁混凝土使用SPH粒子法，由於各個方向上的配筋不同，加勁混凝土需要定義為正交各向異性（使用LAW24），此時需要在/PROP/SPH粒子單元屬性的卡片中使用skew_ID定義相應的材料方向,skew中的局部x軸就是材料方向1，如果沒有定義skew_ID那麼就認為材料是各項同性。

另外材料卡片同樣可以配合/EOS卡片來設置狀態，使用/FAIL失效卡片來設置材料的失效。比如下面的使用SPH粒子法的金屬鋁材料卡片（單位制g-mm-mus）。

通常使用SPH粒子法在Radioss中失效並不是用單元刪除，而是停止計算粒子與其周邊粒子的相互作用。

SPH粒子的接觸可以有粒子和粒子之間的接觸，這個接觸不需要特別的設置，它是一直有的。還有一種是粒子與Lagrangian構件之間的接觸，那就可以用/INTER/TYPE7來設置，此時需要注意下面幾點：

 · Lagrangian的殼體或者實體需要設置為主面master

 · SPH粒子的間距小於Lagrangian的殼體或者實體的網格，否則不能呈現光滑的接觸面，而是會出現下面的凹凸不平的現象。

 · 設置接觸時，Gap是粒子中心到殼的中性面的距離，或者粒子中心到實體單元的外表面的距離。

 · 只能設置Istf=0或者1000，這樣接觸剛度就使用Lagrangian的主面材料的剛度。

還有一種tied接觸也可以設置，就是使用/INTER/TYPE2將貼近Lagrangian的主面粒子設置為tied接觸中的從點（slave），而且必須設置Spotflag = 4，這樣就不會傳遞旋轉自由度。

另外還需要將向距離鐵道接觸表面2h範圍（或更大）的粒子添加到/INTER/TYPE7接觸中，以免這些粒子會與Lagrangian構件在變形過程中有穿透。

SPH粒子法的邊界設置的確是與普通的Lagrangian中的邊界設置有所不同。

比如使用/BCS時，Lagrangian中只要設置最外面的一層點的邊界約束就可以，而SPH粒子法中這種約束還不夠，最好約束2h範圍內的粒子才可以。

通常SPH粒子法中使用/SPHBCS的邊界約束。

使用/SPHBCS可以類似/BCS一樣選取最外面一層的點。

Radioss計算時會產生出相應的一下虛擬點（ghost particles）， 比如下圖中所示，黑色點選取作為邊界，那麼使用/BCS會使得內部的點（紅色）在力的作用下，由於力沒有平衡而越過邊界。

但是如果使用/SPHBCS那麼計算時，Radioss就會產生（藍色的）虛擬點，在受力過程中有很好的力的平衡，而不會越過邊界。

但是如果力太大還是越過了邊界，那麼參數Ilev可以用來定義這些越過邊界的點是直接從計算中剔除，或者使用彈性衝擊方程將這些點反彈回來。

有多少個虛擬點會產生那麼由/SPHGLO中的Nghost 參數決定。另外還可在/SPHBCS中通過參數type定義粒子是在約束平面上固結（tied）還是滑移（slide）的。

除了上面的邊界約束，還可以定義邊界出入口，比如模擬注水和放水，只要在水管入口處使用/SPH/INOUT卡片定就可以。

以水管進出水為例，首先在水管的進出口用殼單元劃分網格，這裡要注意進口處殼單元的法向Normal順著流體進去方向，而出口處殼單元的法向是逆著流體出去的方向。這些殼單元可以用void的材料和單元屬性設置（/MAT/VOID+/PROP/VOID）。

然後在/SPH/INOUT卡片中入口設置Itye=1；而出口設置Itype=2，如果Itype=3那麼就是設置沒有反彈的出口Non-reflective frontiers (NRF)。如果這種Itype=4那麼就是控制流量的出口。

在surf_ID中就使用出入口的殼網格形成的面。

這樣就可以使得入口的位置當注入材料的品質等於粒子的品質時，將創建一個粒子。如何計算注入材料的品質？可以通過定義密度vs時間曲線，以及定義材料沿著入口法向的進入速度曲線來計算。所以需要定義I_Rho這個密度曲線和Iv_n這個速度曲線。

而粒子運動到出口位置，通過計算壓力來決定多少粒子可以出去，所以卡片中需要定義I_ps這個壓力曲線。

如果是選用NRF出口（Itype=3）那麼類似Itype=2需要設置壓力曲線，並且還需要設置lc這個長度，它是用於計算下面的截止頻率fc的，一般取的大一點的數值，比如這個實例中就取整個SPH粒子域裡面的水管模型的長度就可以了。

無論是任何一種Itype的設置都需要Disp_Gap來控制粒子間距，一般Disp_Gap可以設置2h就可以了。

如果需要這個水管進出水的SPH粒子模型可以聯繫我們。

SPH粒子法穩定計算的時間步長同樣也有節點時間步長控制和單元時間步長控制。

單元時間步長計算如下：

通常粒子之間間距越大，穩定計算的時間步長就越大。使用單元時間步長控制（/DT/SPHCEL）時Δtsca一般建議使用0.3。

節點時間步長：

（mi是粒子的品質，

而Ki是粒子之間的剛度有核函式定義）

節點時間步長Δtsca控制時一般建議使用0.67。通常推薦使用節點時間步長控制，它比單元時間步長控制計算更加穩定。

在Radioss中還有一個功能是SOL2SPH，也就是在不失效時它就是普通的實體單元Solid，當實體單元失效時，該實體單元就會變成相應的SPH粒子。這樣的方式比全部使用SPH粒子的方法計算要更加高效。

而且同時具有在大變形區域由於使用了SPH粒子所以單元畸變，沒有負體積等等這些實體單元會碰到的問題，使用SOL2SPH和全部使用SPH的計算精度相當。

那麼在定義了SOLSPH後，到底是如何觸發實體單元到SPH粒子的呢？有下面幾種情況：

1.    當定義了/DT/BRICK/DEL的控制，實體單元的變形導致時間步長得到滿足所定義的最小時間步長那麼這個實體單元不是刪除而是觸發成了SPH粒子。

    

2. 當實體單元滿足所定義的材料失效準則，實體單元不是刪除而是觸發成了SPH粒子。

    

3. 當定義接觸（/INTER/TYPE7）時，當任何一個粒子在它的搜索範圍（詳見前面SPH單元屬性）內搜索到了不同於它的其他部件（part）的沒有觸發的SPH粒子，那麼這個每個觸發SPH粒子的實體單元就被觸發了，裡面的粒子就和這個粒子有了接觸，因此當兩個部件都使用SOL2SPH時，他們直接不需要再特別設置TYPE7的接觸了。

當/SPHGLO中的Isol2sph設置為2時，那就不是不同部件（part）而是不同子集（subset）之間計算粒子接觸而引起的SOL2SPH粒子觸發。

那麼SOL2SPH在模型中的設置也是非常簡單。

以下面為例，如果定義part 2為SOL2SPH，那麼在part 2使用的實體單元屬性（/PROP/SOLID或者/PROP/SOL_ORTH）中定義Ndir和sphpart_ID這兩個參數，在sphpart_ID中定義SPH的part資訊（下面實例中是ID=33的SPH的part），這個SPH的part使用和實體的part一樣的材料ID（下面實例中都使用ID=11的材料卡片）。

在SPH的part中需要定義SPH的單元屬性（下面實例中即ID=15的SPH單元屬性）。只有實體構件的實體單元屬性和SPH單元屬性就定義好了。

這裡實體單元屬性中的Ndir是指實體單元變成SPH粒子時，每個實體單元中每個方向的SPH粒子個數。比如下面的六面體單元中當Ndir=3就是六面體每個方向上（r,s,t）有3個SPH粒子產生。如果是四面體單元當Ndir=2那就是四面體每個方向（r,s,t）有2個SPH粒子產生。

當然目前SOL2SPH的使用時注意一些相容性問題，比如只能用於六面體和四面實體單元，不能用於SPH的inlet和outlet，具體可以參見Radioss的工具書。