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

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

關鍵詞：光滑粒子動力學 非牛頓 數(shù)值模擬 張力不穩(wěn)定

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

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

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

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

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

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

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

2.1 數(shù)值精度低、穩(wěn)定性差

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

2.2 壓力振蕩嚴重

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

2.3 存在張力不穩(wěn)定

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

3 解決方案

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

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

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

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

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

4 結語

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

參考文獻

[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] 劉謀斌，常建忠.光滑粒子動力學方法中粒子分布與數(shù)值穩(wěn)定性分析[J].物理學報，2010（59）：3654-3662.

[6] 蔣濤.聚合物自由面問題的修正對稱SPH方法模擬研究[D].博士研究生學位論文，西安：西北工業(yè)大學，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）方法是一種純無網(wǎng)格粒子方法，目前已廣泛應用于非牛頓流體流動問題。該文在探討了光滑粒子動力學方法應用于非牛頓流體流動問題的研究現(xiàn)狀的基礎上，分析了該方法數(shù)值模擬中存在精度低和張力不穩(wěn)定的問題及產(chǎn)生原因，并提出了相應的解決方案。

關鍵詞：光滑粒子動力學 非牛頓 數(shù)值模擬 張力不穩(wěn)定

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

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

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

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

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

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

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

2.1 數(shù)值精度低、穩(wěn)定性差

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

2.2 壓力振蕩嚴重

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

2.3 存在張力不穩(wěn)定

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

3 解決方案

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

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

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

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

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

4 結語

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

參考文獻

[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] 劉謀斌，常建忠.光滑粒子動力學方法中粒子分布與數(shù)值穩(wěn)定性分析[J].物理學報，2010（59）：3654-3662.

[6] 蔣濤.聚合物自由面問題的修正對稱SPH方法模擬研究[D].博士研究生學位論文，西安：西北工業(yè)大學，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）方法是一種純無網(wǎng)格粒子方法，目前已廣泛應用于非牛頓流體流動問題。該文在探討了光滑粒子動力學方法應用于非牛頓流體流動問題的研究現(xiàn)狀的基礎上，分析了該方法數(shù)值模擬中存在精度低和張力不穩(wěn)定的問題及產(chǎn)生原因，并提出了相應的解決方案。

關鍵詞：光滑粒子動力學 非牛頓 數(shù)值模擬 張力不穩(wěn)定

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

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

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

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

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

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

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

2.1 數(shù)值精度低、穩(wěn)定性差

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

2.2 壓力振蕩嚴重

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

2.3 存在張力不穩(wěn)定

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

3 解決方案

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

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

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

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

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

4 結語

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

參考文獻

[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] 劉謀斌，常建忠.光滑粒子動力學方法中粒子分布與數(shù)值穩(wěn)定性分析[J].物理學報，2010（59）：3654-3662.

[6] 蔣濤.聚合物自由面問題的修正對稱SPH方法模擬研究[D].博士研究生學位論文，西安：西北工業(yè)大學，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