非連續(xù)函數(shù)逼近的SPH方法

周學(xué)君

(1. 黃岡師范學(xué)院 數(shù)理學(xué)院，湖北 黃州 438000；2. 河海大學(xué) 力學(xué)與材料學(xué)院，江蘇 南京 210098)

?

非連續(xù)函數(shù)逼近的SPH方法

周學(xué)君1,2

(1. 黃岡師范學(xué)院 數(shù)理學(xué)院，湖北 黃州 438000；2. 河海大學(xué) 力學(xué)與材料學(xué)院，江蘇 南京 210098)

摘要利用光滑粒子流體動(dòng)力學(xué)(SPH)方法研究非連續(xù)函數(shù)的逼近問題。通過改進(jìn)傳統(tǒng)SPH方法，推導(dǎo)出解決非連續(xù)問題的DSPH近似公式，并成功地應(yīng)用于非連續(xù)函數(shù)的數(shù)值逼近實(shí)例。算例結(jié)果表明，DSPH不僅能準(zhǔn)確確定函數(shù)的非連續(xù)位置，而且近似效果令人滿意。

關(guān)鍵詞非連續(xù)函數(shù)；逼近；光滑粒子流體動(dòng)力學(xué)；數(shù)值模擬

函數(shù)逼近問題一直是函數(shù)論及其他數(shù)學(xué)分支研究的熱點(diǎn)，并在諸多相關(guān)學(xué)科的基礎(chǔ)理論和技術(shù)應(yīng)用中占有重要地位。對(duì)于連續(xù)函數(shù)的近似，已經(jīng)有很多成熟的理論和方法，如插值法、Fourier方法、小波方法和機(jī)器學(xué)習(xí)算法等；但對(duì)于非連續(xù)函數(shù)的逼近，目前研究的方法有限，大都從算子逼近的角度去開展，如高福洪等利用Kantorovich算子對(duì)非連續(xù)函數(shù)進(jìn)行近似[1]，張曉紅等運(yùn)用Bernstein算子對(duì)非連續(xù)函數(shù)逼近[2]，郭順生使用Feller算子逼近具有第一類間斷點(diǎn)的非連續(xù)函數(shù)[3]，劉智新利用廣義Bernstein-Bezier多項(xiàng)式對(duì)非連續(xù)函數(shù)逼近[4]。相比于連續(xù)函數(shù)，非連續(xù)函數(shù)存在不連續(xù)點(diǎn)即間斷點(diǎn)，在間斷點(diǎn)處對(duì)函數(shù)進(jìn)行近似存在著困難，于是如何處理好間斷點(diǎn)成為非連續(xù)函數(shù)逼近的關(guān)鍵。

本文考慮利用一種稱為光滑粒子流體動(dòng)力學(xué)(Smoothed Particle Hydrodynamics，簡稱SPH)的數(shù)值方法來研究非連續(xù)函數(shù)的逼近問題，該方法是由Lucy(1977)[5]、Gingold和Monaghan(1977)[6]提出的一種Lagrangian形式的無網(wǎng)格粒子計(jì)算方法。SPH方法最初提出是用于解決三維開放空間的天體物理學(xué)問題，至今已經(jīng)被廣泛研究和擴(kuò)展。該方法早期主要應(yīng)用在流體力學(xué)的相關(guān)領(lǐng)域，如彈性流、磁流體動(dòng)力學(xué)、多相流、多孔介質(zhì)流、熱傳導(dǎo)等，并取得了很大的成功；隨著對(duì)SPH方法研究的不斷深入，應(yīng)用領(lǐng)域也從計(jì)算流體力學(xué)問題擴(kuò)展到計(jì)算固體力學(xué)問題，如脆性固體斷裂、金屬成型、高速碰撞、炸藥爆炸等大變形和沖擊荷載問題。

傳統(tǒng)的SPH方法在解決連續(xù)性問題時(shí)比較順利，但當(dāng)所解決的問題具有非連續(xù)性，傳統(tǒng)SPH方法遇到瓶頸，需要對(duì)傳統(tǒng)SPH作一些修正，這里把解決非連續(xù)性問題的SPH方法稱為DSPH(Discontinuous Smoothed Particle Hydrodynamics)[7]。閆蕊等已經(jīng)對(duì)DSPH方法的有效性進(jìn)行了驗(yàn)證[8]，Liu等提出一維的DSPH方法，并用來模擬一維的激波管取得了成功[9]，趙燕等利用DSPH方法模擬界面不連續(xù)的物理量，并驗(yàn)證方法的數(shù)值有效性[10]。本文只考慮DSPH模擬非連續(xù)函數(shù)的純數(shù)學(xué)問題，沒有涉及實(shí)際問題的工程或?qū)嶒?yàn)背景。

1光滑粒子流體動(dòng)力學(xué)方法

在SPH方法中，將問題域離散成有限個(gè)具有獨(dú)立質(zhì)量和體積的粒子，每個(gè)粒子都攜帶有速度、密度、應(yīng)力等材料特性；在空間上，每個(gè)粒子通過其影響域(支持域)內(nèi)鄰近粒子進(jìn)行插值近似來計(jì)算它們的材料特性。利用SPH方法近似方程一般分兩步進(jìn)行，第一步是將函數(shù)按照積分加權(quán)的形式進(jìn)行積分近似；第二步是將積分離散化(粒子近似)，應(yīng)用離散化的粒子加權(quán)近似函數(shù)。

對(duì)于函數(shù)f(x)，其在問題域Ω內(nèi)的積分近似表示式為

(1)

式中x,x′為坐標(biāo)向量，且都包含在問題域Ω內(nèi)；dx′表示x處無窮小體元；W稱為光滑函數(shù)或者核函數(shù)；h是光滑長度，表征光滑函數(shù)W的影響區(qū)域。

光滑函數(shù)W在近似過程中，實(shí)際上充當(dāng)權(quán)函數(shù)的作用，它很大程度上決定著SPH方法的精度和穩(wěn)定性。光滑函數(shù)一般都要滿足以下三個(gè)基本條件：

(Ⅰ)正則化條件，即光滑函數(shù)在問題域內(nèi)的積分為1：

(2)

(Ⅱ)當(dāng)光滑函數(shù)的光滑長度趨于零時(shí)，具有delta函數(shù)性質(zhì)：

(3)

(Ⅲ)緊支性條件：當(dāng)|x-x′|>κh時(shí)，W(x-x′,h)=0

其中κ是常數(shù)，確定光滑函數(shù)的有效(非零)范圍，這個(gè)范圍也是點(diǎn)x處光滑函數(shù)的支持域。

光滑函數(shù)的種類有很多，如鐘形函數(shù)、高斯型光滑函數(shù)、三次樣條函數(shù)、四次樣條函數(shù)等。本研究采用現(xiàn)有SPH文獻(xiàn)中最為廣泛應(yīng)用的光滑函數(shù)——三次樣條函數(shù)：

(4)

式中ad是正則化因子，在一維、二維和三維空間中分別為1/h，15/(7πh2)和3/2πh3。R是在點(diǎn)x和x′處兩粒子之間的相對(duì)距離，即

(5)

其中r表示兩點(diǎn)之間的距離。對(duì)于三次樣條函數(shù)而言，常數(shù)κ一般取2。

場函數(shù)f(x)的SPH積分表示式(1)可寫成離散化的粒子近似式：

(6)

這里mj，ρj分別表示粒子j的質(zhì)量和密度，j=1,2,…，N,N為在x處粒子的支持域內(nèi)的粒子總數(shù)。

于是，在粒子i處的場函數(shù)的粒子近似式可寫成

(7)

2非連續(xù)SPH方法

傳統(tǒng)SPH方法建立在對(duì)函數(shù)進(jìn)行Taylor展開的基礎(chǔ)上，要求函數(shù)在整個(gè)支持域上連續(xù)，不能直接應(yīng)用在具有非連續(xù)性的問題中。若對(duì)SPH方法進(jìn)行修正，就得到可以應(yīng)用到非連續(xù)性問題的DSPH方法。

2.1　一維的DSPH公式

關(guān)于一維情形的DSPH公式在文獻(xiàn)[7]中有著詳盡的闡述，這里主要引用該文獻(xiàn)的相關(guān)內(nèi)容。在一維空間上，上述的位置矢量x就是標(biāo)量，于是在本部分中點(diǎn)的位置就用標(biāo)量x表示。考慮在數(shù)軸上，尋求在點(diǎn)xi支持域內(nèi)任意函數(shù)f(x)的積分近似式。假設(shè)支持域?yàn)閰^(qū)間[a,b]，長度為2κh，若函數(shù)f(x)在支持域上的點(diǎn)d是非連續(xù)的，且xi

圖1　在d處非連續(xù)的一維函數(shù)的積分近似

于是在整個(gè)支持域上f(x)與光滑函數(shù)W乘積的積分，可以利用積分區(qū)域的可加性寫成

(8)

將式(8)右邊第一項(xiàng)中的函數(shù)f(x)在點(diǎn)xi處進(jìn)行Taylor展開，第二項(xiàng)中的函數(shù)f(x)在點(diǎn)xk處展開(xk∈[d,b])，可得

(9)

這里r表示Taylor展開式余項(xiàng)。

由核函數(shù)W的緊支性，可知|x-xi|≤κh，|x-xk|≤κh，且假設(shè)f′(x)在[a,d)∪(d,b]上有界，這樣式(9)第三個(gè)等號(hào)右邊的第三、四項(xiàng)可以視作Taylor展開式余項(xiàng)r(h)，于是上式即為

(10)

可得

(11)

在非連續(xù)條件下，式(11)中大括號(hào)中的項(xiàng)其實(shí)就是積分近似式的修正項(xiàng)，若函數(shù)f(x)是連續(xù)的，則修正項(xiàng)可以省略。

下一步將進(jìn)行粒子近似。若將xi的支持域利用粒子進(jìn)行離散，因?yàn)閮蓚€(gè)粒子不能位于同一個(gè)位置，所以非連續(xù)點(diǎn)d總是位于兩個(gè)粒子之間，因?yàn)辄c(diǎn)xk的選取是任意的，為方便起見，不妨將xk設(shè)定在距離非連續(xù)點(diǎn)d最近的右端的粒子處，即xi

圖2　在d處非連續(xù)的一維函數(shù)的粒子近似

將式(11)中的積分離散化后，可得f(xi)的粒子近似式為

(12)

值得注意的是，因?yàn)槭?11)大括號(hào)中分式的分子的積分區(qū)域?yàn)閇d,b]，于是式(12)中相應(yīng)部分的粒子近似式的求和也只在該區(qū)域的粒子上進(jìn)行，如圖2所示。

2.2　多維的DSPH公式

假設(shè)Ω為支持域，且Ω=Ω1+Ω2，Γ是Ω1和Ω2的分界面。函數(shù)f(x)在Γ處不連續(xù)，但在Ω的其他地方連續(xù)，xi為Ω1內(nèi)靠近Γ的點(diǎn)，xk為Ω2內(nèi)任意一點(diǎn)，如圖3所示。

圖3　在Γ處非連續(xù)的二維函數(shù)

同一維情形的推導(dǎo)過程類似，可得函數(shù)f(x)與光滑函數(shù)W(x)的乘積在支持域Ω上的積分為

(13)

對(duì)式(13)采用Taylor展開，進(jìn)行適當(dāng)?shù)淖冃魏罂傻煤瘮?shù)f(x)積分近似表達(dá)式

(14)

將上式積分離散化后，可得f(xi)的粒子近似式為

(15)

3算例實(shí)現(xiàn)

對(duì)于具體的非連續(xù)函數(shù)而言，其在所討論區(qū)域的不連續(xù)點(diǎn)(一維情形中的d，多維情形中的Γ)可以很容易確定；而xk的最佳選取位置是距離非連續(xù)點(diǎn)最近的右端的粒子處，可以在分布好問題域的粒子后再確定。DSPH算法的編程是基于Fortran語言編寫的，工作環(huán)境是Visual Studio 2012。另外，在數(shù)值模擬過程中，每個(gè)粒子的支持域內(nèi)只包含有限個(gè)粒子，這些粒子將會(huì)在粒子近似計(jì)算過程中被使用，這里采用全配對(duì)搜索法來確定給定粒子的支持域內(nèi)粒子，關(guān)于全配對(duì)搜索法的更多細(xì)節(jié)請參見文獻(xiàn)[7]。

例1考慮一維函數(shù)

(16)

容易知道f(x)在x=0.5處不連續(xù)，在其余地方都是連續(xù)的。

利用DSPH方法逼近該函數(shù)之前，需要先布置粒子。將區(qū)間[-1,1]等分成40個(gè)子區(qū)間，每個(gè)子區(qū)間的中心分布一個(gè)粒子，一共40個(gè)粒子，粒子間距為0.05，每個(gè)粒子的體積均為0.05，即

(17)

光滑函數(shù)采用三次樣條函數(shù)，κ=2，光滑長度h=0.5。圖4給出了函數(shù)f(x)的DSPH近似效果，圖4中實(shí)線為函數(shù)f(x)的真實(shí)圖像，○表示DSPH方法的模擬結(jié)果，從結(jié)果來看，DSPH方法不僅能準(zhǔn)確地判斷函數(shù)的不連續(xù)點(diǎn)的位置，而且模擬的精度很好。

圖4　例1的DSPH方法近似結(jié)果

例2考慮二維函數(shù)

(18)

直線x-y+0.02=0是二維函數(shù)f(x)的不連續(xù)界線，在二維空間[0,1]×[0,1]均勻布置11×11個(gè)粒子，粒子在橫向和縱向的間距都是0.1，光滑函數(shù)選用三次樣條函數(shù)，κ=2，光滑長度h=0.2，每個(gè)粒子的體積mi/ρj=1/121≈0.0083。

圖5和圖6分別給出函數(shù)f(x,y)的真實(shí)圖像和DSPH方法的模擬圖像，可以發(fā)現(xiàn)隨著維數(shù)的增加，DSPH的模擬能力依舊很好，從兩者的形態(tài)來看，DSPH方法對(duì)不連續(xù)處的把握很精準(zhǔn)，同時(shí)圖像的吻合度很高。

圖5　例2的函數(shù)的真實(shí)圖像

圖6　例2的函數(shù)的DSPH方法模擬圖像

傳統(tǒng)SPH方法解決非連續(xù)問題遇到困難，通過對(duì)它進(jìn)行改進(jìn)，推導(dǎo)出適合非連續(xù)問題的DSPH公式，并分別給出一維和多維情形的粒子近似式。結(jié)合兩個(gè)具體的算例，利用DSPH方法對(duì)一維和二維函數(shù)進(jìn)行逼近模擬，結(jié)果表明，DSPH方法不僅可以準(zhǔn)確地找到非連續(xù)點(diǎn)的位置，而且逼近效果很好。

參考文獻(xiàn):

[1]高福洪，王子玉.Kantorovich算子對(duì)不連續(xù)函數(shù)的逼近[J].高等學(xué)校計(jì)算數(shù)學(xué)學(xué)報(bào)，1988，3：274-280.

[2]張曉紅，施禮明.變形Bernstein算子對(duì)不連續(xù)函數(shù)的逼近[J].華東冶金學(xué)院學(xué)報(bào)，1995，12(1)：98-102.

[3]郭順生.用Feller算子逼近第一類間斷點(diǎn)的函數(shù)[J].應(yīng)用數(shù)學(xué)學(xué)報(bào)，1991，14(1)：57-65.

[4]劉智新.廣義Bernstein-Bezier多項(xiàng)式對(duì)非連續(xù)函數(shù)的逼近[J].數(shù)學(xué)季刊，1987，2(3)：100-109.

[5]Lucy L B.Numerical approach to testing of fission hypothesis[J].Astronomical Journal，1977， 82：1013-1024.

[6]Gingold R A, Moraghan J J.Smoothed particle hydrodynamics: Theory and application to non-spherical stars[J].Manual Notebook of Royal Astronomical Society, 1977，181：375-389.

[7]韓旭，楊剛，強(qiáng)洪夫(譯).光滑粒子流體動(dòng)力學(xué)——種無網(wǎng)格粒子法[M].長沙：湖南大學(xué)出版社，2005.

[8]閆蕊，徐緋，張?jiān)狼?DSPH方法的有效性驗(yàn)證及應(yīng)用[J].爆炸與沖擊，2013，33(2)：133-139.

[9]Liu M B，Liu G R，Lam K Y.A one-dimensional meshfree particle figuretion for simulating shock waves[J].Shock Waves，2003，13(3)：201-211.

[10]趙燕，徐緋，李玉龍．一種可考慮界面不連續(xù)的改進(jìn)SPH方法[J]．計(jì)算力學(xué)學(xué)報(bào)，2009，26(6)：928-934.

編輯王菊平

Approximation of discontinuous function based on SPH method

ZHOU Xue-jun1,2

(1.College of Mathematics and Physics, Huanggang Normal University, Huangzhou 438000, Hubei, China;

2.College of Mechanics and Materials, Hohai University, Nanjing 210098, Jiangsu, China)

AbstractThe approximation of discontinuous function based on Smoothed Particle Hydrodynamics (SPH) is discussed. The approximate figures of Discontinuous Smoothed Particle Hydrodynamics (DSPH) for discontinuous problem are derived by improving the original SPH, and DSPH is applied successfully to the numerical approximation examples of discontinuous function. The results show that DSPH can not only find accurately the position of discontinuous, but also achieve satisfactory approximation.

Key wordsdiscontinuous function; approximation; Smoothed Particle Hydrodynamics (SPH); numerical simulation

基金項(xiàng)目湖北省教育廳科學(xué)技術(shù)研究項(xiàng)目(B2015218)，黃岡師范學(xué)院科學(xué)研究項(xiàng)目(2013020103)。

作者簡介周學(xué)君，男，湖北蘄春人，講師，博士生，主要研究方向?yàn)闄C(jī)器學(xué)習(xí)理論及其應(yīng)用、力學(xué)中數(shù)值方法等。

收稿日期2015-08-08

doi10.3969/j.issn.1003-8078.2015.06.04

中圖分類號(hào)O302

文獻(xiàn)標(biāo)志碼A

文章編號(hào)1003-8078(2015)06-0014-05