淺談SPH方法在非牛頓自由表面流中的應用

摘 要：光滑粒子動力學（SPH）方法是一種純無網格粒子方法，目前已廣泛應用于非牛頓流體流動問題。該文在探討了光滑粒子動力學方法應用于非牛頓流體流動問題的研究現狀的基礎上，分析了該方法數值模擬中存在精度低和張力不穩定的問題及產生原因，并提出了相應的解決方案。

關鍵詞：光滑粒子動力學 非牛頓 數值模擬 張力不穩定

中圖分類號：O35 文獻標識碼：A 文章編號：1674-098X（2014）07（a）-0250-02

1 光滑粒子動力學方法應用簡介

光滑粒子動力學（Smoothed Particle Hydrodynamics，SPH）方法是一種基于Lagrange描述的純無網格粒子方法。它不需要使用任何網格，而是將連續體離散成有限個粒子，這些粒子攜帶所有物理量（密度、壓力、速度、溫度等）的信息。其次通過核函數將連續介質的控制方程轉化為積分形式。這樣整個流場變成了一系列粒子的表達，每個粒子的物理量都由周圍粒子的物理量插值得到。

SPH方法最初由Lucy[1]、Gingold和Monaghan[2]提出并用于求解三維開放空間天體物理問題。1994年，Monaghan[3]首次將SPH方法用于自由表面流的數值模擬。隨后，許多學者運用SPH方法對牛頓流體的自由表面流動問題[4-5]進行了模擬研究，大大推動了SPH方法在流體力學領域的發展。目前，SPH方法在聚合物自由面方面的研究也取得了一定進展。2005年，Ellero和Tanner[6]運用SPH方法模擬了低Reynolds數下粘彈性熔體剪切流動問題。2007年、2008年Rafiee[6]運用SPH方法模擬了非牛頓熔體的自由表面流問題。上述研究均是對簡單流動進行研究。2006年，Fang等人[6]首次采用SPH方法研究了粘彈性液滴落到剛性固壁上的復雜流動問題，討論了SPH方法中的人工應力參數的選取。2012年，本文作者等[6]對傳統SPH方法進行改進的基礎上研究了聚合物自由面問題。2010年，Fan等人[6]首次運用SPH方法研究了高粘性冪律流體的充模問題。

總之，隨著計算機技術的快速發展，SPH方法在聚合物流動問題中的應用越來越廣泛。

2 光滑粒子動力學方法中存在的缺點

與其他方法相比，SPH方法也是計算流體力學領域中一種較新的數值方法。隨著應用范圍越來越廣，SPH方法的一些缺點也逐漸顯現出來。SPH方法的缺點主要體現以下幾方面[4，7-8]。

2.1 數值精度低、穩定性差

SPH方法在粒子分布均勻的區域內具有二階精度，但在粒子分布不均勻區域或邊界處至多具有一階精度。經研究發現，SPH方法精度低的主要原因是：在粒子分布不均勻區域，某流體粒子影響域內的粒子數比較少；而對于邊界附近，流體粒子的影響域被邊界截斷，邊界外不存在粒子，從而出現比較嚴重的粒子缺失現象。而流動過程中必然出現粒子分布不均勻的情況。提高SPH方法的精度成為促進SPH方法進一步發展的重要因素。

2.2 壓力振蕩嚴重

SPH方法雖然在模擬自由表面流動問題中具有優勢，尤其是對于大變形問題。然而，由SPH方法得到的壓力卻存在嚴重振蕩。這是因為：除了上述SPH方法數值精度比較低的原因，另一個原因就是弱可壓SPH方法中始終采用壓力狀態方程來求解壓力。雖然有學者針對此問題提出了改進措施，如Colagrossi[6]采用密度初始化方法、耗散項方法等在一定程度上改善了壓力振蕩的情況，但是上述方法始終采用壓力狀態方程來求解壓力，無法從根本上解決壓力振蕩問題。

2.3 存在張力不穩定

張力不穩定性或稱“拉伸不穩定性”，該術語來源于傳統SPH方法求解固體力學中材料拉伸問題時出現的一種非物理粒子簇集現象[7]。許多學者[6]在固體變形問題的數值模擬中對其出現的原因進行了研究和分析，但迄今為止關于張力不穩定性產生的本質原因沒有一個統一的理論分析，也未能夠提出普遍應用的克服方法。目前，應用比較廣泛的克服張力不穩定性的方法是Monaghan和Gray等人[6]提出的人工應力法。

3 解決方案

上述SPH方法的缺點嚴重限制了SPH方法的推廣與發展。因此，對SPH方法的改進有著重要的理論意義和應用價值。

針對傳統SPH方法數值精度低的問題，眾所周知，核函數在SPH方法中扮演著舉足輕重的角色。因此，本文借助Taylor展開方法[4，8]，在核函數導數的左側乘以一個修正矩陣，且此修正矩陣中不需要求解核函數導數，以此來提高SPH方法的數值精度，即

針對壓力振蕩嚴重的問題，本文考慮結合不可壓縮SPH方法和弱可壓SPH方法的優點，借助于有限差分方法的思想，對壓力進行迭代求解，直至求得的壓力計算得到的密度與初始密度的誤差在允許范圍之內時，迭代終止。

針對張力不穩定性問題，本文根據Swegle等人[7]提出的張力不穩定性結論，給出了一種容易施加的人工應力方法，其具體離散形式為

其中，為一個正常數（），為流體的應力張量分量，核函數的二階導數分量為。上述人工應力方法離散的形式比已有的人工應力形式簡單，且同時考慮了拉伸與壓縮兩種情況引起的不穩定。

4 結語

本文首先對SPH方法應用于聚合物流動問題的研究現狀進行了闡述，然后根據相關文獻及作者的數值模擬經驗分析了運用SPH方法模擬聚合物流動問題時存在精度低和張力不穩定性問題及產生的原因，最后針對存在的問題，提出了相應的解決方案。

參考文獻

[1] L.B.Lucy.A numerical approach to the testing of the ssion hypothesis[J].Astronomical Journal，1977（83）：1013-1024.

[2] R.A.Gingold，J.J.Monaghan. Smoothed particle hydrodynamics theory and application to non-spherical stars[J].Monthly Notices of the Royal Astronomical Society，1977（181）：375-389.

[3] J.J. Monaghan.Simulating free surface ows with SPH[J].Journal of Computational Physics，1994（110）：399-406.

[4] G.R.Liu，M.B.Liu.Smoothed Particle Hydrodynamics：A Mesh-free Particle Method[M].Singapore：World Scientific，2003.

[5] 劉謀斌，常建忠.光滑粒子動力學方法中粒子分布與數值穩定性分析[J].物理學報，2010（59）：3654-3662.

[6] 蔣濤.聚合物自由面問題的修正對稱SPH方法模擬研究[D].博士研究生學位論文，西安：西北工業大學，2012.

[7] J.Swegle，W.D.L.Hicks， S.W.Attaway.Smooth particle hydrodynamics stability analysis[J].Journal of Computational Physics，1995， （116）：123-134.

[8] Jiang T，Ouyang J，Ren J，Yang B，Xu X.A mixed corrected symmetric SPH（MC-SSPH）method for computational dynamic problems[J].Computer Physics Communications，2012（183）：50-62.endprint

摘 要：光滑粒子動力學（SPH）方法是一種純無網格粒子方法，目前已廣泛應用于非牛頓流體流動問題。該文在探討了光滑粒子動力學方法應用于非牛頓流體流動問題的研究現狀的基礎上，分析了該方法數值模擬中存在精度低和張力不穩定的問題及產生原因，并提出了相應的解決方案。

關鍵詞：光滑粒子動力學 非牛頓 數值模擬 張力不穩定

中圖分類號：O35 文獻標識碼：A 文章編號：1674-098X（2014）07（a）-0250-02

1 光滑粒子動力學方法應用簡介

光滑粒子動力學（Smoothed Particle Hydrodynamics，SPH）方法是一種基于Lagrange描述的純無網格粒子方法。它不需要使用任何網格，而是將連續體離散成有限個粒子，這些粒子攜帶所有物理量（密度、壓力、速度、溫度等）的信息。其次通過核函數將連續介質的控制方程轉化為積分形式。這樣整個流場變成了一系列粒子的表達，每個粒子的物理量都由周圍粒子的物理量插值得到。

SPH方法最初由Lucy[1]、Gingold和Monaghan[2]提出并用于求解三維開放空間天體物理問題。1994年，Monaghan[3]首次將SPH方法用于自由表面流的數值模擬。隨后，許多學者運用SPH方法對牛頓流體的自由表面流動問題[4-5]進行了模擬研究，大大推動了SPH方法在流體力學領域的發展。目前，SPH方法在聚合物自由面方面的研究也取得了一定進展。2005年，Ellero和Tanner[6]運用SPH方法模擬了低Reynolds數下粘彈性熔體剪切流動問題。2007年、2008年Rafiee[6]運用SPH方法模擬了非牛頓熔體的自由表面流問題。上述研究均是對簡單流動進行研究。2006年，Fang等人[6]首次采用SPH方法研究了粘彈性液滴落到剛性固壁上的復雜流動問題，討論了SPH方法中的人工應力參數的選取。2012年，本文作者等[6]對傳統SPH方法進行改進的基礎上研究了聚合物自由面問題。2010年，Fan等人[6]首次運用SPH方法研究了高粘性冪律流體的充模問題。

總之，隨著計算機技術的快速發展，SPH方法在聚合物流動問題中的應用越來越廣泛。

2 光滑粒子動力學方法中存在的缺點

與其他方法相比，SPH方法也是計算流體力學領域中一種較新的數值方法。隨著應用范圍越來越廣，SPH方法的一些缺點也逐漸顯現出來。SPH方法的缺點主要體現以下幾方面[4，7-8]。

2.1 數值精度低、穩定性差

SPH方法在粒子分布均勻的區域內具有二階精度，但在粒子分布不均勻區域或邊界處至多具有一階精度。經研究發現，SPH方法精度低的主要原因是：在粒子分布不均勻區域，某流體粒子影響域內的粒子數比較少；而對于邊界附近，流體粒子的影響域被邊界截斷，邊界外不存在粒子，從而出現比較嚴重的粒子缺失現象。而流動過程中必然出現粒子分布不均勻的情況。提高SPH方法的精度成為促進SPH方法進一步發展的重要因素。

2.2 壓力振蕩嚴重

SPH方法雖然在模擬自由表面流動問題中具有優勢，尤其是對于大變形問題。然而，由SPH方法得到的壓力卻存在嚴重振蕩。這是因為：除了上述SPH方法數值精度比較低的原因，另一個原因就是弱可壓SPH方法中始終采用壓力狀態方程來求解壓力。雖然有學者針對此問題提出了改進措施，如Colagrossi[6]采用密度初始化方法、耗散項方法等在一定程度上改善了壓力振蕩的情況，但是上述方法始終采用壓力狀態方程來求解壓力，無法從根本上解決壓力振蕩問題。

2.3 存在張力不穩定

張力不穩定性或稱“拉伸不穩定性”，該術語來源于傳統SPH方法求解固體力學中材料拉伸問題時出現的一種非物理粒子簇集現象[7]。許多學者[6]在固體變形問題的數值模擬中對其出現的原因進行了研究和分析，但迄今為止關于張力不穩定性產生的本質原因沒有一個統一的理論分析，也未能夠提出普遍應用的克服方法。目前，應用比較廣泛的克服張力不穩定性的方法是Monaghan和Gray等人[6]提出的人工應力法。

3 解決方案

上述SPH方法的缺點嚴重限制了SPH方法的推廣與發展。因此，對SPH方法的改進有著重要的理論意義和應用價值。

針對傳統SPH方法數值精度低的問題，眾所周知，核函數在SPH方法中扮演著舉足輕重的角色。因此，本文借助Taylor展開方法[4，8]，在核函數導數的左側乘以一個修正矩陣，且此修正矩陣中不需要求解核函數導數，以此來提高SPH方法的數值精度，即

針對壓力振蕩嚴重的問題，本文考慮結合不可壓縮SPH方法和弱可壓SPH方法的優點，借助于有限差分方法的思想，對壓力進行迭代求解，直至求得的壓力計算得到的密度與初始密度的誤差在允許范圍之內時，迭代終止。

針對張力不穩定性問題，本文根據Swegle等人[7]提出的張力不穩定性結論，給出了一種容易施加的人工應力方法，其具體離散形式為

其中，為一個正常數（），為流體的應力張量分量，核函數的二階導數分量為。上述人工應力方法離散的形式比已有的人工應力形式簡單，且同時考慮了拉伸與壓縮兩種情況引起的不穩定。

4 結語

本文首先對SPH方法應用于聚合物流動問題的研究現狀進行了闡述，然后根據相關文獻及作者的數值模擬經驗分析了運用SPH方法模擬聚合物流動問題時存在精度低和張力不穩定性問題及產生的原因，最后針對存在的問題，提出了相應的解決方案。

參考文獻

[1] L.B.Lucy.A numerical approach to the testing of the ssion hypothesis[J].Astronomical Journal，1977（83）：1013-1024.

[2] R.A.Gingold，J.J.Monaghan. Smoothed particle hydrodynamics theory and application to non-spherical stars[J].Monthly Notices of the Royal Astronomical Society，1977（181）：375-389.

[3] J.J. Monaghan.Simulating free surface ows with SPH[J].Journal of Computational Physics，1994（110）：399-406.

[4] G.R.Liu，M.B.Liu.Smoothed Particle Hydrodynamics：A Mesh-free Particle Method[M].Singapore：World Scientific，2003.

[5] 劉謀斌，常建忠.光滑粒子動力學方法中粒子分布與數值穩定性分析[J].物理學報，2010（59）：3654-3662.

[6] 蔣濤.聚合物自由面問題的修正對稱SPH方法模擬研究[D].博士研究生學位論文，西安：西北工業大學，2012.

[7] J.Swegle，W.D.L.Hicks， S.W.Attaway.Smooth particle hydrodynamics stability analysis[J].Journal of Computational Physics，1995， （116）：123-134.

[8] Jiang T，Ouyang J，Ren J，Yang B，Xu X.A mixed corrected symmetric SPH（MC-SSPH）method for computational dynamic problems[J].Computer Physics Communications，2012（183）：50-62.endprint

摘 要：光滑粒子動力學（SPH）方法是一種純無網格粒子方法，目前已廣泛應用于非牛頓流體流動問題。該文在探討了光滑粒子動力學方法應用于非牛頓流體流動問題的研究現狀的基礎上，分析了該方法數值模擬中存在精度低和張力不穩定的問題及產生原因，并提出了相應的解決方案。

關鍵詞：光滑粒子動力學 非牛頓 數值模擬 張力不穩定

中圖分類號：O35 文獻標識碼：A 文章編號：1674-098X（2014）07（a）-0250-02

1 光滑粒子動力學方法應用簡介

光滑粒子動力學（Smoothed Particle Hydrodynamics，SPH）方法是一種基于Lagrange描述的純無網格粒子方法。它不需要使用任何網格，而是將連續體離散成有限個粒子，這些粒子攜帶所有物理量（密度、壓力、速度、溫度等）的信息。其次通過核函數將連續介質的控制方程轉化為積分形式。這樣整個流場變成了一系列粒子的表達，每個粒子的物理量都由周圍粒子的物理量插值得到。

SPH方法最初由Lucy[1]、Gingold和Monaghan[2]提出并用于求解三維開放空間天體物理問題。1994年，Monaghan[3]首次將SPH方法用于自由表面流的數值模擬。隨后，許多學者運用SPH方法對牛頓流體的自由表面流動問題[4-5]進行了模擬研究，大大推動了SPH方法在流體力學領域的發展。目前，SPH方法在聚合物自由面方面的研究也取得了一定進展。2005年，Ellero和Tanner[6]運用SPH方法模擬了低Reynolds數下粘彈性熔體剪切流動問題。2007年、2008年Rafiee[6]運用SPH方法模擬了非牛頓熔體的自由表面流問題。上述研究均是對簡單流動進行研究。2006年，Fang等人[6]首次采用SPH方法研究了粘彈性液滴落到剛性固壁上的復雜流動問題，討論了SPH方法中的人工應力參數的選取。2012年，本文作者等[6]對傳統SPH方法進行改進的基礎上研究了聚合物自由面問題。2010年，Fan等人[6]首次運用SPH方法研究了高粘性冪律流體的充模問題。

總之，隨著計算機技術的快速發展，SPH方法在聚合物流動問題中的應用越來越廣泛。

2 光滑粒子動力學方法中存在的缺點

與其他方法相比，SPH方法也是計算流體力學領域中一種較新的數值方法。隨著應用范圍越來越廣，SPH方法的一些缺點也逐漸顯現出來。SPH方法的缺點主要體現以下幾方面[4，7-8]。

2.1 數值精度低、穩定性差

SPH方法在粒子分布均勻的區域內具有二階精度，但在粒子分布不均勻區域或邊界處至多具有一階精度。經研究發現，SPH方法精度低的主要原因是：在粒子分布不均勻區域，某流體粒子影響域內的粒子數比較少；而對于邊界附近，流體粒子的影響域被邊界截斷，邊界外不存在粒子，從而出現比較嚴重的粒子缺失現象。而流動過程中必然出現粒子分布不均勻的情況。提高SPH方法的精度成為促進SPH方法進一步發展的重要因素。

2.2 壓力振蕩嚴重

SPH方法雖然在模擬自由表面流動問題中具有優勢，尤其是對于大變形問題。然而，由SPH方法得到的壓力卻存在嚴重振蕩。這是因為：除了上述SPH方法數值精度比較低的原因，另一個原因就是弱可壓SPH方法中始終采用壓力狀態方程來求解壓力。雖然有學者針對此問題提出了改進措施，如Colagrossi[6]采用密度初始化方法、耗散項方法等在一定程度上改善了壓力振蕩的情況，但是上述方法始終采用壓力狀態方程來求解壓力，無法從根本上解決壓力振蕩問題。

2.3 存在張力不穩定

張力不穩定性或稱“拉伸不穩定性”，該術語來源于傳統SPH方法求解固體力學中材料拉伸問題時出現的一種非物理粒子簇集現象[7]。許多學者[6]在固體變形問題的數值模擬中對其出現的原因進行了研究和分析，但迄今為止關于張力不穩定性產生的本質原因沒有一個統一的理論分析，也未能夠提出普遍應用的克服方法。目前，應用比較廣泛的克服張力不穩定性的方法是Monaghan和Gray等人[6]提出的人工應力法。

3 解決方案

上述SPH方法的缺點嚴重限制了SPH方法的推廣與發展。因此，對SPH方法的改進有著重要的理論意義和應用價值。

針對傳統SPH方法數值精度低的問題，眾所周知，核函數在SPH方法中扮演著舉足輕重的角色。因此，本文借助Taylor展開方法[4，8]，在核函數導數的左側乘以一個修正矩陣，且此修正矩陣中不需要求解核函數導數，以此來提高SPH方法的數值精度，即

針對壓力振蕩嚴重的問題，本文考慮結合不可壓縮SPH方法和弱可壓SPH方法的優點，借助于有限差分方法的思想，對壓力進行迭代求解，直至求得的壓力計算得到的密度與初始密度的誤差在允許范圍之內時，迭代終止。

針對張力不穩定性問題，本文根據Swegle等人[7]提出的張力不穩定性結論，給出了一種容易施加的人工應力方法，其具體離散形式為

其中，為一個正常數（），為流體的應力張量分量，核函數的二階導數分量為。上述人工應力方法離散的形式比已有的人工應力形式簡單，且同時考慮了拉伸與壓縮兩種情況引起的不穩定。

4 結語

本文首先對SPH方法應用于聚合物流動問題的研究現狀進行了闡述，然后根據相關文獻及作者的數值模擬經驗分析了運用SPH方法模擬聚合物流動問題時存在精度低和張力不穩定性問題及產生的原因，最后針對存在的問題，提出了相應的解決方案。

參考文獻

[1] L.B.Lucy.A numerical approach to the testing of the ssion hypothesis[J].Astronomical Journal，1977（83）：1013-1024.

[2] R.A.Gingold，J.J.Monaghan. Smoothed particle hydrodynamics theory and application to non-spherical stars[J].Monthly Notices of the Royal Astronomical Society，1977（181）：375-389.

[3] J.J. Monaghan.Simulating free surface ows with SPH[J].Journal of Computational Physics，1994（110）：399-406.

[4] G.R.Liu，M.B.Liu.Smoothed Particle Hydrodynamics：A Mesh-free Particle Method[M].Singapore：World Scientific，2003.

[5] 劉謀斌，常建忠.光滑粒子動力學方法中粒子分布與數值穩定性分析[J].物理學報，2010（59）：3654-3662.

[6] 蔣濤.聚合物自由面問題的修正對稱SPH方法模擬研究[D].博士研究生學位論文，西安：西北工業大學，2012.

[7] J.Swegle，W.D.L.Hicks， S.W.Attaway.Smooth particle hydrodynamics stability analysis[J].Journal of Computational Physics，1995， （116）：123-134.

[8] Jiang T，Ouyang J，Ren J，Yang B，Xu X.A mixed corrected symmetric SPH（MC-SSPH）method for computational dynamic problems[J].Computer Physics Communications，2012（183）：50-62.endprint