楊賀奎，丁　寧，徐江榮

(杭州電子科技大學(xué)能源研究所，浙江 杭州 310018)

?

基于函數(shù)梯度無(wú)網(wǎng)格自適應(yīng)方法研究

楊賀奎，丁寧，徐江榮

(杭州電子科技大學(xué)能源研究所，浙江 杭州 310018)

摘要：以計(jì)算一維泊松方程為例，通過(guò)研究節(jié)點(diǎn)分布對(duì)無(wú)網(wǎng)格法數(shù)值模擬的影響，提出了一種基于函數(shù)梯度大小調(diào)整節(jié)點(diǎn)密度的自適應(yīng)方法.首先利用無(wú)網(wǎng)格法求出函數(shù)的近似梯度函數(shù)，然后分析該函數(shù)的梯度值，在梯度較大的位置加密計(jì)算節(jié)點(diǎn)，以提高計(jì)算精度，減少數(shù)值模擬的誤差.利用該方法計(jì)算二維圓柱繞流N-S方程，模擬計(jì)算的結(jié)果基本符合圓柱繞流現(xiàn)象，表明了基于函數(shù)梯度進(jìn)行節(jié)點(diǎn)加密的方法達(dá)到了很好的效果.

關(guān)鍵詞：梯度函數(shù)；無(wú)網(wǎng)格方法；節(jié)點(diǎn)分布；RPIM形函數(shù)

0引言

在解決工程數(shù)值分析領(lǐng)域中，傳統(tǒng)網(wǎng)格法的運(yùn)用出現(xiàn)了許多難于求解的問(wèn)題.特別是對(duì)于氣固兩相流運(yùn)動(dòng)界面不連續(xù)的問(wèn)題，模擬時(shí)往往產(chǎn)生不能求解網(wǎng)格極端扭曲或難于跟蹤自由表面等難題.網(wǎng)格法在求解計(jì)算時(shí)，需要不斷地重新劃分網(wǎng)格，使得計(jì)算量增大，計(jì)算精度降低，穩(wěn)定性變差.為了解決這一難題，文獻(xiàn)[1]提出利用無(wú)網(wǎng)格法來(lái)解決相關(guān)問(wèn)題.無(wú)網(wǎng)格法在數(shù)值計(jì)算中不需要將節(jié)點(diǎn)連接成網(wǎng)格單元進(jìn)行劃分和重構(gòu)，而是按照一些任意分布的坐標(biāo)點(diǎn)構(gòu)造插值函數(shù)離散控制方程，能方便地模擬各種復(fù)雜形狀的流場(chǎng)[2].因此無(wú)網(wǎng)格法為計(jì)算流體力學(xué)的自適應(yīng)分析、界面追蹤、多相流動(dòng)等方面開(kāi)辟了一個(gè)極具誘惑的新途徑.然而運(yùn)用無(wú)網(wǎng)格法計(jì)算有一個(gè)很大的問(wèn)題，其計(jì)算量幾乎和計(jì)算節(jié)點(diǎn)數(shù)成正比，隨著節(jié)點(diǎn)數(shù)的增加，計(jì)算量迅速增大[3].如何保證在計(jì)算總節(jié)點(diǎn)數(shù)一定的情況下，提高計(jì)算精度成為一個(gè)難題.本文以計(jì)算一維泊松方程為例，通過(guò)研究節(jié)點(diǎn)分布對(duì)無(wú)網(wǎng)格法流場(chǎng)數(shù)值模擬的影響，提出了一種基于函數(shù)梯度大小調(diào)整節(jié)點(diǎn)密度的自適應(yīng)方法.

1基于函數(shù)梯度大小進(jìn)行節(jié)點(diǎn)加密分析

無(wú)網(wǎng)格法可以任意地在計(jì)算區(qū)域內(nèi)布點(diǎn)，不同的布點(diǎn)方式給模擬計(jì)算精度帶來(lái)不同的影響.一般情況下，函數(shù)梯度反映了函數(shù)變化的趨勢(shì)，在高梯度的位置函數(shù)的變化比較大，因此在這些位置應(yīng)該多分布一些節(jié)點(diǎn)，如果分布的節(jié)點(diǎn)比較少必然會(huì)增大數(shù)值模擬的誤差[4].

1.1一維泊松方程

下面以計(jì)算一維泊松方程為例，使用配點(diǎn)法模擬，研究如何進(jìn)行節(jié)點(diǎn)加密，使得在計(jì)算總節(jié)點(diǎn)數(shù)一定的情況下，得到較高的模擬計(jì)算精度.其解析函數(shù)梯度絕對(duì)值分布如圖1所示.

由圖1可以看出，前半部分比后半部分的梯度小，這說(shuō)明該函數(shù)值在前半部分變化要比后半部分變化小.因此為了減少計(jì)算誤差，在函數(shù)梯度較大的位置進(jìn)行節(jié)點(diǎn)加密.為了驗(yàn)證這一想法，本文先對(duì)前半部分節(jié)點(diǎn)加密，再對(duì)后半部分節(jié)點(diǎn)加密，比較這兩種情況下數(shù)值模擬的平均誤差.加密的方案是在原來(lái)均勻分布的平均間距上等間距地增加一些節(jié)點(diǎn).

圖1　解析函數(shù)梯度絕對(duì)值分布圖

1)前半部分的節(jié)點(diǎn)加密

前半部分節(jié)點(diǎn)加密后，平均誤差隨節(jié)點(diǎn)個(gè)數(shù)變化如圖2所示.

圖2　前半部分節(jié)點(diǎn)加密引起的平均誤差隨節(jié)點(diǎn)個(gè)數(shù)變化圖

從圖2(a)中可以看到，在前半部分平均間距中，計(jì)算節(jié)點(diǎn)數(shù)增加到8時(shí)，平均誤差迅速猛增.

為了進(jìn)一步體現(xiàn)加密效果，圖2(b)去掉了圖2(a)中Nthick_num=8時(shí)Daverage_error的數(shù)值.可見(jiàn)，在平均間距中增加1～2個(gè)節(jié)點(diǎn)，其平均誤差并沒(méi)有下降，增加1個(gè)節(jié)點(diǎn)還導(dǎo)致平均誤差變大.但繼續(xù)增加節(jié)點(diǎn)，平均誤差就減小了.當(dāng)Nthick_num>6時(shí)，Daverage_error比Nthick_num=0時(shí)增大了很多，說(shuō)明加密過(guò)度使得平均誤差變大，并不能使計(jì)算精度有所提升.因此，要進(jìn)行局部節(jié)點(diǎn)的加密，就需要經(jīng)過(guò)一個(gè)度量的過(guò)程，過(guò)度加密是沒(méi)有意義的.

2)后半部分的節(jié)點(diǎn)加密

后半部分節(jié)點(diǎn)加密后，高梯度處平均誤差隨節(jié)點(diǎn)個(gè)數(shù)變化如圖3所示.

從圖3(a)中可以看出，在后半部分平均間距中，計(jì)算節(jié)點(diǎn)數(shù)增加到5時(shí)，平均誤差迅速猛增.

為了進(jìn)一步體現(xiàn)加密效果，圖3(b)去掉了圖3(a)中Nthick_num=5時(shí)Daverage_error的數(shù)值.當(dāng)Nthick_num>6時(shí)，Daverage_error的數(shù)值比Nthick_num=0時(shí)大，這也說(shuō)明了加密過(guò)度使得平均誤差變大.就本算例而言，在平均間距中增加1～2個(gè)節(jié)點(diǎn)就使得精度有所提升.由圖3(b)可以看出，在平均間距中增加1～2個(gè)節(jié)點(diǎn)效果是最好的，總體效果也比前半部分加密好.這說(shuō)明在梯度比較大的位置加密比在梯度小的位置加密效果更好，精度更高.

圖3　后半部分節(jié)點(diǎn)加密引起的高梯度處平均誤差隨節(jié)點(diǎn)個(gè)數(shù)變化圖

1.2根據(jù)近似函數(shù)計(jì)算函數(shù)梯度方法

圖4　近似函數(shù)梯度絕對(duì)值分布圖

一般情況下，進(jìn)行數(shù)值模擬時(shí)，不知道所求函數(shù)的解析解，也就不知道該函數(shù)的梯度函數(shù)，故而無(wú)法確定函數(shù)值變化的大小.本文提出使用近似函數(shù)的方法，即先在均勻分布的節(jié)點(diǎn)上算出數(shù)值解，再利用形函數(shù)的導(dǎo)數(shù)，近似求出函數(shù)梯度值.

下面驗(yàn)證根據(jù)近似函數(shù)計(jì)算函數(shù)梯度的方法是可行的.近似函數(shù)梯度絕對(duì)值分布如圖5所示.

由圖4可以看出，利用近似函數(shù)求解的梯度函數(shù)的數(shù)值解和解析解基本上重合的，因此可以認(rèn)為這種方法是可行的.

1.3不同位置加密方法的函數(shù)誤差

在泊松方程下不同加密方法的函數(shù)誤差分布如圖5所示.

圖5　泊松方程不同加密方法下函數(shù)誤差分布圖

圖5中，Derror=|u(xi)real-u(xi)h|，即表示每1個(gè)節(jié)點(diǎn)的相對(duì)誤差.由圖5可以看出，后1/4加密函數(shù)誤差相對(duì)較小，前1/4加密函數(shù)誤差相對(duì)較大，所以在函數(shù)梯度大的位置加密效果比較好，在函數(shù)梯度小的位置加密效果不太好.

2基于不同節(jié)點(diǎn)加密影響域PDE的模擬計(jì)算

為了進(jìn)一步驗(yàn)證基于函數(shù)梯度進(jìn)行節(jié)點(diǎn)加密的方法的可行性，采用無(wú)網(wǎng)格徑向基插值法(RPIM)對(duì)二維圓柱繞流N-S方程進(jìn)行模擬驗(yàn)證.在該工況中，圓柱周圍流場(chǎng)的速度變化是最大的，函數(shù)梯度也是相對(duì)較大的，因此要在圓柱周圍進(jìn)行節(jié)點(diǎn)加密.

2.1圓柱繞流N-S方程

2.1.1數(shù)學(xué)模型

圓柱繞流的主控制方程采用原始變量形式的二維非定常不可壓縮流N-S方程，在慣性參考系下其無(wú)量綱形式為.

2.1.2計(jì)算區(qū)域和邊界條件

設(shè)定計(jì)算區(qū)域，水平方向?yàn)閤方向，豎直方向?yàn)閥方向.水平長(zhǎng)度為1 m，豎直長(zhǎng)度為0.5 m.圓柱圓心坐標(biāo)為(0.2,0.2),進(jìn)口速度U=0.01 m/s.流動(dòng)方向與垂直流動(dòng)方向上的均勻速度分別定義為u和v.

進(jìn)口處：給定進(jìn)口速度，u=U,v=0;

出口處：在計(jì)算區(qū)域出口處，設(shè)置自由出口條件.假定在任意時(shí)刻下計(jì)算區(qū)域出口處的流動(dòng)都已充分發(fā)展，那么各速度分量和壓力均取第二類邊界條件，即

上下邊界：為消除固壁帶來(lái)的尺寸效應(yīng)，給定u=U,v=0;

圓柱表面：給定無(wú)滑移，無(wú)穿透邊界條件，u=0,v=0;

二維非定常不可壓縮流N-S方程是非線性方程，含有拋物型和橢圓型兩種性質(zhì).一般的解法是不行的，本文采用數(shù)值方法—速度修正法即分裂步數(shù)法進(jìn)行迭代求解，時(shí)間離散采用差分方法，空間離散采用無(wú)網(wǎng)格徑向基插值法(RPIM)[5].在整個(gè)流場(chǎng)中，只有在圓柱周圍的梯度變化是最大的，因此在圓柱周圍進(jìn)行節(jié)點(diǎn)加密，形成一個(gè)云狀結(jié)構(gòu)，節(jié)點(diǎn)分布如圖6如示.

圖6　圓柱繞流節(jié)點(diǎn)分布圖

無(wú)網(wǎng)格徑向基插值法(RPIM)屬于全局弱式的積分方法.使用該方法需要1個(gè)背景網(wǎng)格，并在每個(gè)網(wǎng)格單元內(nèi)均勻地分布高斯點(diǎn).影響該方法計(jì)算精度的因素主要是影響域半徑、時(shí)間步長(zhǎng)、影響域內(nèi)保留的節(jié)點(diǎn)個(gè)數(shù).這3項(xiàng)數(shù)值的大小不僅影響到數(shù)值計(jì)算的精度，而且很大程度上影響了計(jì)算的時(shí)間.因此要求影響域半徑盡量小，影響域內(nèi)保留的節(jié)點(diǎn)個(gè)數(shù)保持在20以內(nèi).為提高數(shù)值計(jì)算的精度，時(shí)間步長(zhǎng)采用變步長(zhǎng)的形式進(jìn)行求解.時(shí)間步長(zhǎng)序列為

下面對(duì)控制方程進(jìn)行離散：

通過(guò)以上分析，求解圓柱繞流N-S方程，根據(jù)自適應(yīng)分析理論，在圓柱周圍進(jìn)行節(jié)點(diǎn)加密，形成一個(gè)節(jié)點(diǎn)密度較高的云狀結(jié)構(gòu)，而在其他位置節(jié)點(diǎn)分布比較稀疏.

2.2模擬計(jì)算結(jié)果分析

利用Tecplot軟件引入數(shù)值求解圓柱繞流N-S方程，在圓柱周圍進(jìn)行節(jié)點(diǎn)加密，當(dāng)圓柱繞流速度場(chǎng)變化到某個(gè)時(shí)刻時(shí)，圓柱繞流速度矢量發(fā)生的變化情況如圖7所示.

圖7　圓柱繞流速度矢量變化情況

圖7中，當(dāng)圓柱繞流速度場(chǎng)變化40 s時(shí)圓柱后方開(kāi)始形成漩渦；80 s時(shí)已經(jīng)形成了漩渦，而且在圓柱的后方還可以看到脫落的漩渦.120 s時(shí)流場(chǎng)已經(jīng)基本充分發(fā)展，在圓柱的后方已有明顯的小漩渦，脫落的漩渦也很明顯.這個(gè)結(jié)果基本符合圓柱繞流的現(xiàn)象，表明了基于函數(shù)梯度進(jìn)行節(jié)點(diǎn)加密的方法是可行的.

3結(jié)束語(yǔ)

本文研究了基于函數(shù)梯度大小調(diào)整節(jié)點(diǎn)密度的自適應(yīng)方法，利用無(wú)網(wǎng)格方法先求出近似梯度函數(shù)，再根據(jù)梯度函數(shù)來(lái)確定節(jié)點(diǎn)加密的位置.得出在高梯度的位置進(jìn)行節(jié)點(diǎn)加密，可以減少數(shù)值模擬誤差，提高計(jì)算精度，但注意不要過(guò)度加密，否則會(huì)適得其反.并通過(guò)計(jì)算二維圓柱繞流N-S方程進(jìn)行模擬驗(yàn)證，模擬結(jié)果基本符合圓柱繞流的現(xiàn)象，驗(yàn)證了本文方法的有效性.

參考文獻(xiàn)

[1]李九紅,程玉民.無(wú)網(wǎng)格方法的研究進(jìn)展與展望[J].力學(xué)季刊,2006,27(1)：143-152.

[2]LIU Y, ZHANG X, LU M W. Meshless Least-Squares Method for Solving the Steady-State Heat Conduction Equation[J]. Tsinghua Science and Technology, 2005，10(1):61-66.

[3]仇軼,由長(zhǎng)福,祁海鷹,等.無(wú)網(wǎng)格方法中的結(jié)點(diǎn)分布算法[J].數(shù)值計(jì)算與計(jì)算機(jī)應(yīng)用,2006,27(3):176-182.

[4]張雄,劉巖,馬上.無(wú)網(wǎng)格法的理論及應(yīng)用[J].力學(xué)進(jìn)展，2009，39(1)：1-36.

[5]夏茂輝,賈延,劉才.基于徑向基函數(shù)的點(diǎn)插值(RPIM)無(wú)網(wǎng)格法[J].燕山大學(xué)學(xué)報(bào)，2006，32(2):112-117.

The Research of the Adaptive Method about Adjusting the Node Density Based on the Meshless Method

YANG Hekui, DING Ning, XU Jiangrong

(InstituteofEnergy,HangzhouDianziUniversity,HangzhouZhejiang310018,China)

Abstract：This paper is based on calculating the Poisson equation as an example, through the study of node distribution effect on numerical simulation in the flow field by using the meshless method, presents an adaptive method about adjusting the node density based on the function gradient size. It is using meshless method to calculate the approximate gradient functions, and analyze the gradient value of the function. Using the multi-distribution of some nodes in the location of function gradient higher, to improve the calculation accuracy and reduce the error of the numerical simulation. By using this method to calculate the N-S equation of the flow around a circular cylinder, the simulation results are basically in accord with the flow around the cylinder. It is indicated that the method of node encryption based on functional gradient has achieved good results.

Key words：the gradient function; the meshless method; the node distribution; the RPIM shape function

DOI：10.13954/j.cnki.hdu.2016.04.018

收稿日期：2015-10-19

基金項(xiàng)目：國(guó)家自然科學(xué)基金資助項(xiàng)目(51176044)；浙江省自然科學(xué)基金資助項(xiàng)目(Y1111085)

作者簡(jiǎn)介：楊賀奎(1992-)，男，安徽亳州人，碩士研究生，能源機(jī)械.通訊作者：丁寧副教授，E-mail:tinin@hdu.edu.cn.

中圖分類號(hào)：O241

文獻(xiàn)標(biāo)識(shí)碼：A

文章編號(hào)：1001-9146(2016)04-0084-06