程 麗 張 超, 肖道林, 張趙威
(1、沈陽大學(xué)機(jī)械工程學(xué)院,遼寧 沈陽110044 2、中國科學(xué)院沈陽自動(dòng)化研究所,遼寧 沈陽110016)
由于科技的發(fā)展和人類對(duì)未知領(lǐng)域探索的需要,力學(xué)仿真廣泛存在于人類科學(xué)研究中。但是傳統(tǒng)力學(xué)仿真所用網(wǎng)格法具有一定局限性。[3-4]
近年,新興的光滑無網(wǎng)格的方法經(jīng)過已經(jīng)成熟,光滑粒子流體力學(xué)(SPH)是由核近似和粒子近似兩步組成的無網(wǎng)格拉格朗日粒子方法[5]。由于SPH 是在不同時(shí)刻計(jì)算該區(qū)域任意分布粒子上承載的物理量,所以具有極強(qiáng)的自適應(yīng)性,并且各個(gè)物理量均為矢量?;赟PH 方法此種特性,對(duì)其一維進(jìn)行研究可以為高維的工程問題提供參考。
在SPH 中,通過光滑核函數(shù)引進(jìn)積分表示式[1]。W(xi-xj,h)稱為光滑核函數(shù),該函數(shù)具有狄拉克函數(shù)的性質(zhì)。其中h 表示光滑長度,? 為積分域,在SPH 中也叫支持域。A(x)j表示物理量A 在坐標(biāo)xj處的值。x 是空間坐標(biāo)向量,i 與j 表示粒子序號(hào)。
實(shí)際工程中流體粒子運(yùn)動(dòng)會(huì)遇到無法穿透的邊界,模擬時(shí)為防止粒子穿透邊界,采用一種Lennard-Jones 邊界力模型
式中,i 和j 分別指邊界粒子和內(nèi)部流體粒子,n1n2為常數(shù),通常取12 和4。D 為最大速度的平方,r0為截止距離,通常取初始距離一半,粒子運(yùn)動(dòng)到邊界粒子r0范圍內(nèi)將受到邊界力的作用。
為了求得粒子在每一個(gè)時(shí)間步運(yùn)動(dòng)狀態(tài),引入蛙跳法。根據(jù)蛙跳積分的方案,在時(shí)間步n 有
圖1 0.005s 粒子運(yùn)動(dòng)分布狀態(tài)
圖2 0.01975s 粒子運(yùn)動(dòng)分布狀態(tài)
圖3 0.025s 粒子分布狀態(tài)
圖4 速度隨時(shí)間變化曲線
仿真采用顆粒密度1.5kg/m,將每一個(gè)粒子看為半徑為0.01m。每一個(gè)時(shí)間步取得到廣泛應(yīng)用的5*10-5s,共計(jì)算500 步。仿真結(jié)果如下,粒子呈現(xiàn)三角形為向右移動(dòng),呈現(xiàn)圓形為向左移動(dòng),正方形為邊界。
經(jīng)過1、圖2、圖3 可以觀測到,0.005s 時(shí)刻所有粒子還是保持向右移動(dòng),0.01975s 開始出現(xiàn)由于邊界力作用而向左移動(dòng)的粒子,最終500 步時(shí)即為0.025s 接近一半粒子向左移動(dòng),另一半粒子還未運(yùn)動(dòng)到邊界繼續(xù)向右移動(dòng),圖4 對(duì)以上加以驗(yàn)證。
最終得出SPH 可以模擬一維顆粒狀樣本變化,面對(duì)高維的工程問題只需提升維度如加上環(huán)境所需要的外力項(xiàng)與粘性項(xiàng)即可。