龔小權(quán) 吳曉軍 唐靜 李明 張健

(中國空氣動(dòng)力研究與發(fā)展中心 計(jì)算空氣動(dòng)力研究所, 綿陽 621000)

間斷Galerkin(discontinuous Galerkin,DG)有限元[1-2]方法作為一種適用于非結(jié)構(gòu)網(wǎng)格的高精度方法,具有計(jì)算精度高,計(jì)算模板緊致,適合處理復(fù)雜外形及邊界、初值問題,并行能力強(qiáng),自適應(yīng)計(jì)算方便等優(yōu)點(diǎn)。

與其他高精度格式類似,DG 方法在模擬包含強(qiáng)間斷流場時(shí)存在收斂性和魯棒性差的問題。目前,大部分學(xué)者通過添加人工黏性[3-4]或限制器[5-6]來抑制激波附近的震蕩。 人工黏性方法作為一種先驗(yàn)方法,在控制方程中添加代表黏性的二階導(dǎo)數(shù)項(xiàng),該項(xiàng)數(shù)值與計(jì)算網(wǎng)格、模擬問題相關(guān),具有一定經(jīng)驗(yàn)性,其對數(shù)值解的干擾沒有統(tǒng)一的評估標(biāo)準(zhǔn)。 限制器作為一種后驗(yàn)方法,通過限制或重構(gòu)多項(xiàng)式高階自由度抑制數(shù)值震蕩,其存在殘差收斂滯止、結(jié)果收斂性差的問題。 目前,DG 方法的激波捕捉技術(shù)已成為其發(fā)展的瓶頸之一,國內(nèi)外在限制器方面開展了大量研究,包括Barth-Jespersen 限制器[7]、Venkatakrishnan[8]限制器,以及Shu 等[9]按照WENO 重構(gòu)思想在三角形網(wǎng)格上提出的HWENO 限制器。

激波前后密度、壓力變化大,在數(shù)值模擬前由于無法預(yù)測激波位置,初始網(wǎng)格一般為均勻分布網(wǎng)格,激波附近網(wǎng)格分辨率不足進(jìn)一步加劇DG方法在模擬強(qiáng)間斷問題時(shí)的收斂性和魯棒性問題。 基于以上原因,基于網(wǎng)格自適應(yīng)和限制器技術(shù)開展DG 方法的間斷流場數(shù)值模擬是重要研究方向,能夠提高計(jì)算精度,改善激波模擬分辨率。

耦合流動(dòng)特征、動(dòng)態(tài)調(diào)整網(wǎng)格單元分布的r型、h 型網(wǎng)格自適應(yīng)技術(shù)能夠改善氣動(dòng)特性計(jì)算精度,提高復(fù)雜流動(dòng)模擬精度[10-11],其在誤差修正[12]、網(wǎng)格無關(guān)性[13]、驗(yàn)證與確認(rèn)[14]、非定常流動(dòng)模擬[15]等方面都開展了應(yīng)用研究。 Karen等[16]在數(shù)值模擬飛行器再入過程中采用基于網(wǎng)格點(diǎn)溫度的r 型網(wǎng)格點(diǎn)移動(dòng)方法,引入網(wǎng)格點(diǎn)周圍單元的最小內(nèi)切圓半徑及基于雷諾數(shù)的網(wǎng)格點(diǎn)目標(biāo)長度作為控制參數(shù)。 從結(jié)果看到,該方法實(shí)現(xiàn)了r 型網(wǎng)格自適應(yīng),在強(qiáng)激波附近聚集2 ～3 層網(wǎng)格。 吳澤艷[17]、孫強(qiáng)[18]等將h 型網(wǎng)格自適應(yīng)方法應(yīng)用在二維DG 的激波捕捉方法中,采用HWENO 限制器捕捉激波,有效實(shí)現(xiàn)了激波區(qū)域的網(wǎng)格自適應(yīng)加密。 王利和周偉江[19]基于伴隨方法開展了網(wǎng)格自適應(yīng)技術(shù)的二維DG 方法研究,提高了計(jì)算精度。 Woopen[20]、Fidkowski[21-22]等在DG 方法上建立了基于伴隨方法發(fā)展網(wǎng)格自適應(yīng)技術(shù),并將其用于阻力預(yù)測,取得了較好的效果。 r 型網(wǎng)格自適應(yīng)方法[23-24]作為一種激波捕捉和誤差估計(jì)方法,通過移動(dòng)網(wǎng)格點(diǎn)使網(wǎng)格聚集在流動(dòng)變化劇烈區(qū)域,實(shí)現(xiàn)對流動(dòng)特征區(qū)域網(wǎng)格加密,改善網(wǎng)格分辨率,提高流場模擬精度和氣動(dòng)特性求解精度。 r 型網(wǎng)格自適應(yīng)的最大優(yōu)點(diǎn)是網(wǎng)格單元的連接關(guān)系和網(wǎng)格的空間拓?fù)浣Y(jié)構(gòu)在自適應(yīng)迭代過程中始終保持不變,流場解算器在運(yùn)行過程中不需要調(diào)整網(wǎng)格數(shù)據(jù)結(jié)構(gòu),不必考慮并行動(dòng)態(tài)負(fù)載平衡問題。

本文針對DG 方法在數(shù)值模擬強(qiáng)間斷方面存在的計(jì)算收斂性差、魯棒性低的問題,建立了一種基于網(wǎng)格點(diǎn)壓力和網(wǎng)格點(diǎn)相對位移變化為控制參數(shù)的r 型網(wǎng)格自適應(yīng)方法及Venkatakrishnan 限制器來提高DG 方法模擬強(qiáng)間斷流場的能力。 重點(diǎn)介紹了r 型網(wǎng)格自適應(yīng)獨(dú)特的網(wǎng)格點(diǎn)移動(dòng)驅(qū)動(dòng)策略;Venkatakrishnan 限制器考慮了DG 方法高斯數(shù)值積分的特殊性。 首先,給出了DG 方法,包括控制方程、空間離散、時(shí)間離散和Venkatakrishnan限制器;其次,介紹了r 型網(wǎng)格自適應(yīng)方法,詳細(xì)給出了網(wǎng)格移動(dòng)策略;再次,采用NACA0012 翼型跨聲速繞流驗(yàn)證DG 方法及限制器;然后,數(shù)值模擬了并列NACA0012 翼型超聲速流動(dòng),驗(yàn)證了二維流動(dòng)下本文r 型網(wǎng)格自適應(yīng)和限制器技術(shù);最后,模擬了三維雙半球-圓柱流動(dòng)干擾流動(dòng),驗(yàn)證了基于r 型網(wǎng)格自適應(yīng)的DG 方法在激波捕捉方面的能力。

1 間斷Galerkin 有限元方法

1.1 控制方程

1.2 DG 方法空間離散

空間離散前,整個(gè)求解區(qū)域Ω被劃分為互不重疊的非結(jié)構(gòu)單元Ωk,對于三維問題,Ωk包括三棱柱、四面體、六面體和金字塔。

將方程(1)兩端乘以測試函數(shù)?j(x,y,z),再在單元內(nèi)積分,并使用分部積分方法得到DG 的弱形式。

式中:n為單元邊界面的外法向;?Ω為單元Ω的邊界;Γ為單元Ω的邊界?Ω的微分符號。

假設(shè)變量在單元內(nèi)具有以下多項(xiàng)式分布:

1.3 DG 方法時(shí)間離散

對方程(4)采用一階歐拉后差:

本文采用LU-SGS 隱式離散[25]求解該線性系統(tǒng),求解過程中引入局部時(shí)間步長、變CFL 數(shù)等加速收斂技術(shù)。

1.4 Venkatakrishnan 限制器

為保證計(jì)算魯棒性,在間斷區(qū)域及網(wǎng)格質(zhì)量較差區(qū)域,采用限制器來抑制震蕩DG 方法的限制器可能帶來計(jì)算精度降低、殘值收斂滯止問題。DG 方法作為一種內(nèi)自由度方法,具有模板緊致特性,這要求限制過程也能保持緊致特性。 有限體積的斜率類限制器通過改進(jìn)后能用于DG 方法,其核心思想是:限制單元變量梯度,使單元任何一點(diǎn)變量值不超過相鄰單元的最值,防止出現(xiàn)“過沖”現(xiàn)象。 Barth 和Jespersen[7]建立了非結(jié)構(gòu)網(wǎng)格中常用的一類Barth-Jespersen 限制器,該限制器在計(jì)算限制因子過程中引入不可微函數(shù)帶來殘值收斂性問題。 Venkatakrishnan[8]通過近似min(1,y) 函數(shù)避免不可微過程,建立了Venkatakrishnan 限 制 器。 Venkatakrishnan 限 制 器 由于其較好的殘值收斂性、魯棒性和計(jì)算精度而得到了廣泛應(yīng)用。 本文參考有限體積法中這類限制器的構(gòu)造原理,得到適合DG 方法的斜率限制器。

判斷當(dāng)前單元所有邊界面的高斯積分點(diǎn)變量值是否滿足以下條件:

積分點(diǎn)限制因子αj具體形式如下:

其中:K為量級為1 的自由可調(diào)參數(shù),其大小影響殘值收斂性,值越小代表限制越大;h為單元的參考長度。

韓國的風(fēng)水說主要應(yīng)用于墓地、部落等的選址以及都城、寺剎、住宅等建造方面。 就墓地而言，中國道教風(fēng)水學(xué)說傳入朝鮮半島后，古代朝鮮人民也認(rèn)為如果把祖先葬在藏風(fēng)納氣的風(fēng)水寶地上，這種“佳吉?dú)狻本蜁?huì)感應(yīng)到子孫后代身上，讓子孫繁衍昌盛； 如果祖先所葬之地沒有“氣”，則會(huì)出現(xiàn)天壤之別的情形。

2 r 型網(wǎng)格自適應(yīng)方法

r 型網(wǎng)格自適應(yīng)方法通過移動(dòng)網(wǎng)格點(diǎn)使網(wǎng)格聚集在流動(dòng)變化劇烈的區(qū)域,實(shí)現(xiàn)流動(dòng)特征區(qū)域網(wǎng)格加密,改善網(wǎng)格分辨率,進(jìn)而提高流場模擬精度。 網(wǎng)格自適應(yīng)技術(shù)主要包括誤差估計(jì)、網(wǎng)格加密/稀疏、表面網(wǎng)格投影、空間體網(wǎng)格匹配及并行實(shí)現(xiàn)。

本文r 型網(wǎng)格自適應(yīng)方法在自適應(yīng)過程中限制表面網(wǎng)格點(diǎn)移動(dòng),因此不考慮表面網(wǎng)格投影。r 型網(wǎng)格自適應(yīng)在自適應(yīng)過程中網(wǎng)格單元的連接關(guān)系和網(wǎng)格空間拓?fù)浣Y(jié)構(gòu)保持不變,也不考慮并行問題。 本文重點(diǎn)介紹r 型網(wǎng)格自適應(yīng)網(wǎng)格點(diǎn)移動(dòng)策略。

Karen 等[16]采用溫度作為網(wǎng)格點(diǎn)移動(dòng)控制變量適合高超聲速流動(dòng),對于一般的跨聲速或超聲速流動(dòng)基于壓力變化可能更具有普適性。 從文獻(xiàn)[16]算例看到,基于單一的溫度作為控制參數(shù),網(wǎng)格自適應(yīng)區(qū)域的網(wǎng)格聚集數(shù)較少,經(jīng)過100 次自適應(yīng)迭代后,只有2 ～3 排網(wǎng)格點(diǎn)移動(dòng),激波附近網(wǎng)格點(diǎn)移動(dòng)信息向其周圍網(wǎng)格點(diǎn)傳遞不明顯。

本文r 型網(wǎng)格自適應(yīng)網(wǎng)格點(diǎn)移動(dòng)的驅(qū)動(dòng)力由網(wǎng)格點(diǎn)的流場變量值確定,以歸一化網(wǎng)格點(diǎn)的變量差及網(wǎng)格點(diǎn)相對位移量這2 類參數(shù)作為網(wǎng)格點(diǎn)移動(dòng)的權(quán)函數(shù)。 針對激波問題,本文以壓力P作為網(wǎng)格點(diǎn)移動(dòng)控制變量。 網(wǎng)格點(diǎn)移動(dòng)的計(jì)算流場如圖1 所示。

圖1 r 型網(wǎng)格自適應(yīng)網(wǎng)格點(diǎn)移動(dòng)位移計(jì)算流程Fig.1 Calculation process of grid point displacement in r-grid adaptive

以圖2 為例具體說明網(wǎng)格點(diǎn)移動(dòng)控制過程。

圖2 網(wǎng)格點(diǎn)分布示意圖Fig.2 Distribution diagram of grid points

1) 采用距離權(quán)計(jì)算網(wǎng)格點(diǎn)0 及周圍網(wǎng)格點(diǎn)1、2、3、4、5 的壓力。

7) 限制特殊點(diǎn)的位移,如對稱面上網(wǎng)格點(diǎn)y方向位移為零(對稱面在y=0 平面)。

8) 移動(dòng)網(wǎng)格后檢查網(wǎng)格質(zhì)量,如果出現(xiàn)負(fù)體積,則修改移動(dòng)量Δr0,Δ?r0=C·Δr0,C為限制系數(shù),本文取0.5,并再次檢查網(wǎng)格質(zhì)量,如果出現(xiàn)負(fù)體積,再按Δ?r0=C·Δr0減小移動(dòng)量,直到體積正常。

3 數(shù)值算例

本文所有算例都在中國空氣動(dòng)力研究與發(fā)展中心的計(jì)算機(jī)群上完成。 首先,采用NACA0012翼型驗(yàn)證本文DG 方法和Venkatakrishnan 限制器;然后,采用典型超聲速二維及三維流動(dòng)算例展示r 型網(wǎng)格自適應(yīng)方法在DG 方法的激波捕捉技術(shù)中的應(yīng)用能力。 算例表明,本文限制器方法具有較好的魯棒性,能夠處理各向異性網(wǎng)格及網(wǎng)格尺度急劇變化的網(wǎng)格。

3.1 DG 方法及其限制器驗(yàn)證

圖3 NACA0012 翼型計(jì)算網(wǎng)格Fig.3 Computational grid of NACA0012 airfoil

圖4 不同方法壓力分布比較Fig.4 Comparison of pressure distribution between different methods

3.2 并列兩段NACA0012 翼型超聲速繞流

本節(jié)數(shù)值模擬了2 個(gè)并列放置的NACA0012翼型的超聲速繞流。 計(jì)算的初始網(wǎng)格如圖5 所示,遠(yuǎn)場為30 倍弦長,共25 982 個(gè)三棱柱單元,展向拉伸一排網(wǎng)格。 其中一個(gè)翼型相對另外一個(gè)翼型上移和后移0. 5 m,計(jì)算來流馬赫數(shù)Ma=2.0,來流迎角0°,來流靜壓P=26 500 Pa,來流靜溫T=223.25 K,自適應(yīng)迭代11 次,在自適應(yīng)迭代中保持翼型表面網(wǎng)格點(diǎn)不移動(dòng)。

翼型錯(cuò)位放置,翼型間產(chǎn)生復(fù)雜的激波干擾,上部翼型的脫體激波在下部翼型的上表面產(chǎn)生反射激波后反射到上部翼型的下表面。 圖5 ～圖9分別給出了初始網(wǎng)格、第1 次自適應(yīng)、第3 次自適應(yīng)、第5 次自適應(yīng)、第11 次自適應(yīng)后的網(wǎng)格分布及壓力系數(shù)空間分布。 可以看到,初始網(wǎng)格和自適應(yīng)后的網(wǎng)格都捕捉到脫體激波及翼型間的反射激波;r 型網(wǎng)格自適應(yīng)方法準(zhǔn)確判別到流場劇烈變化區(qū)域,并將網(wǎng)格沿著激波方向各向異性自適應(yīng)聚合,自適應(yīng)后的激波捕捉更清晰和銳利;隨著自適應(yīng)迭代次數(shù)的增加,激波附近的網(wǎng)格逐步被自適應(yīng)加密,網(wǎng)格更緊密聚合于激波處,約有4 ～5排網(wǎng)格,捕捉激波的網(wǎng)格寬度更小,拉伸比更大。自適應(yīng)方法同樣也準(zhǔn)確識別翼型之間的反射激波,并沿激波方向進(jìn)行了自適應(yīng)加密。

圖7 NACA0012 翼型第3 次自適應(yīng)后計(jì)算網(wǎng)格及流場Fig.7 Computational grid and flowfield of NACA0012 airfoil after third r-adaptation

圖8 NACA0012 翼型第5 次自適應(yīng)后計(jì)算網(wǎng)格及流場Fig.8 Computational grid and flowfield of NACA0012 airfoil after fifth r-adaptation

圖9 NACA0012 翼型第11 次自適應(yīng)后計(jì)算網(wǎng)格及流場Fig.9 Computational grid and flowfield of NACA0012 airfoil after eleventh r-adaptation

圖10 給出了第11 次網(wǎng)格自適應(yīng)后翼型脫體激波附近網(wǎng)格放大圖。 自適應(yīng)后,激波附近網(wǎng)格拉伸比變大,網(wǎng)格尺寸劇烈變化,網(wǎng)格質(zhì)量變差,本文的Venkatakrishnan 限制器能夠很好處理該類網(wǎng)格,保證DG 方法計(jì)算的魯棒性。 圖11 給出了空間截面位置,圖12 和圖13 分別給出了空間截面z=0. 25 m 和z=0. 068 m 處的壓力分布?？梢钥吹?自適應(yīng)前后DG 方法都能夠捕捉激波,隨著r 型網(wǎng)格自適應(yīng)迭代,網(wǎng)格逐步向激波位置聚集,DG 方法捕捉激波更加清晰,主要表現(xiàn)在跨過激波前后壓力峰值更大,激波寬度更窄,激波更銳利。 從圖12 看到,在自適應(yīng)迭代11 次后,數(shù)值模擬的激波變化已經(jīng)很小。

圖10 自適應(yīng)網(wǎng)格局部放大圖Fig.10 Local enlarged view of adaptation grid

圖11 空間截面位置Fig.11 Space section position

圖12 z =0.25 m 處壓力分布及激波處局部放大圖Fig.12 Pressure distribution at section z =0.25 m and local enlarged view of shock wave

圖13 z =0.068 m 處壓力分布及激波處局部放大圖Fig.13 Pressure distribution at section z =0.068 m and local enlarged view of shock wave

圖14 展示了自適應(yīng)計(jì)算過程中翼型阻力系數(shù)CD隨計(jì)算步數(shù)變化曲線。 阻力系數(shù)在每次網(wǎng)格自適應(yīng)后只需較小計(jì)算步就收斂,隨著自適應(yīng)迭代增加,阻力系數(shù)增加量逐漸減小。

圖14 自適應(yīng)過程阻力系數(shù)隨計(jì)算步數(shù)變化曲線Fig.14 Curve of drag coefficient varying with iteration steps in adaptive process

圖15 給出了阻力系數(shù)和Venkatakrishnan 限制器系數(shù)平均值隨自適應(yīng)迭代步變化曲線。 阻力系數(shù)隨網(wǎng)格自適應(yīng)迭代逐步增大,斜率逐步減小,限制系數(shù)平均值及其斜率隨網(wǎng)格自適應(yīng)迭代步數(shù)都逐步減小。 網(wǎng)格自適應(yīng)步數(shù)越大,網(wǎng)格拉伸比越大,網(wǎng)格分布呈現(xiàn)更加不均勻性,因此限制系數(shù)會(huì)逐漸減小。 但由于網(wǎng)格主要在激波附近發(fā)生移動(dòng),而數(shù)值模擬激波時(shí)計(jì)算自身就要降階,此處的網(wǎng)格移動(dòng)對計(jì)算精度影響很小,算例驗(yàn)證了本文r 網(wǎng)格自適應(yīng)方法、DG 方法及Venkatakrishnan 限制器的魯棒性。

圖15 阻力系數(shù)及限制器系數(shù)平均值隨自適應(yīng)步數(shù)變化Fig.15 Resistance and average value of limitation coefficients with adaptive steps

3.3 三維雙半球-圓柱流動(dòng)干擾模擬

為進(jìn)一步驗(yàn)證本文DG 方法的限制器及網(wǎng)格自適應(yīng)技術(shù),本節(jié)計(jì)算了三維雙半球-圓柱流動(dòng)干擾模擬算例。 計(jì)算網(wǎng)格如圖16 所示,共527 612 個(gè)四面體單元。 計(jì)算來流馬赫數(shù)Ma=2.0,來流迎角0°,來流靜壓P= 26 500 Pa,來流靜溫T=223.25 K,自適應(yīng)迭代4 次。 圖16 ～圖18 給出了雙半球-圓柱對稱面及表面初始網(wǎng)格和流場、第2次自適應(yīng)后的網(wǎng)格及流場、第4 次自適應(yīng)后的網(wǎng)格及流場。 隨著自適應(yīng)迭代步的增加,網(wǎng)格逐步向激波附近聚集,用于捕捉激波的網(wǎng)格的寬度減小,拉伸比增大,激波分辨率增高。

圖16 三維雙半球-圓柱初始計(jì)算網(wǎng)格及流場Fig.16 Initial computational grid and flowfield of 3D double hemispherical cylinder

圖17 雙半球-圓柱第2 次自適應(yīng)計(jì)算網(wǎng)格及流場Fig.17 Computational grid and flowfield of double hemispherical cylinder after second r-adaptation

圖18 雙半球-圓柱第4 次自適應(yīng)計(jì)算網(wǎng)格及流場Fig.18 Computational grid and flowfield of double hemispherical cylinder after fourth r-adaptation

圖19 ～圖21 給出了x=0.22 m 處的初始網(wǎng)格及流場、第2 次自適應(yīng)后的網(wǎng)格及流場、第4 次自適應(yīng)后的網(wǎng)格及流場。 從圖21 看到,r 型網(wǎng)格自適應(yīng)后空間網(wǎng)格向激波附近聚集,來自上部圓柱的脫體激波打到下部圓柱,其激波分辨率提高。

圖19 雙半球-圓柱初始x =0.22 m 處空間網(wǎng)格及流場Fig.19 Spatial grid and flowfield at x =0.22 m section of double hemispherial cylinder

圖22 給出了自適應(yīng)前后對稱面上z=0.04 m截面壓力分布。 自適應(yīng)后激波寬度減小,激波前后壓力峰值更大。 圖23 和圖24 分別給出了自適應(yīng)前后下部圓柱和上部圓柱y=0 截面表面壓力分布。 自適應(yīng)后影響上下圓柱壓力分布,由于本文的r 型網(wǎng)格自適應(yīng)未移動(dòng)表面網(wǎng)格,影響量較小。 圖25 給出了x=0.22 m 截面處下部圓柱的表面壓力分布。 上部圓柱的脫體激波在此處打到下部圓柱,看到自適應(yīng)前后該處的壓力有明顯變化,自適應(yīng)后的表面壓力更低。

圖22 自適應(yīng)前后對稱面z =0.04 m 截面壓力分布Fig.22 Pressure distribution at z =0.04 m section of symmetry plane before and after adaptation

圖23 自適應(yīng)前后下部圓柱y =0 截面表面壓力分布Fig.23 Surface pressure distribution of lower cylinder at y =0 section before and after adaptation

圖24 自適應(yīng)前后上部圓柱y =0 截面表面壓力分布Fig.24 Surface pressure distribution of upper cylinder at y =0 section before and after adaptation

圖25 自適應(yīng)前后下部圓柱x =0.22 m 截面表面壓力分布Fig.25 Surface pressure distribution of lower cylinder at x =0.22 m section before and after adaptotion

4 結(jié) 論

本文針對DG 方法在模擬強(qiáng)間斷流場面臨的問題,發(fā)展了Venkatakrishnan 限制器和r 型網(wǎng)格自適應(yīng)方法,通過二維及三維算例分析不同自適應(yīng)迭代步后的表面和空間壓力分布,展示了本文方法在激波精細(xì)捕捉方面的能力,并能得到如下結(jié)論:

1) r 型網(wǎng)格自適應(yīng)方法具備根據(jù)流場特征開展網(wǎng)格自適應(yīng)能力,能夠沿著激波方向?qū)崿F(xiàn)網(wǎng)格的各向異性自適應(yīng)移動(dòng)加密,具有較好的魯棒性,保證網(wǎng)格不出現(xiàn)負(fù)體積。

2) DG 方法及Venkatakrishnan 限制器針對超聲速問題具有較好的激波捕捉能力,能夠魯棒處理各向異性非結(jié)構(gòu)網(wǎng)格、針對自適應(yīng)后的相鄰網(wǎng)格尺寸差異很大的網(wǎng)格具有較高的計(jì)算魯棒性。

3) 通過在DG 方法上發(fā)展r 型網(wǎng)格自適應(yīng)方法能夠明顯改善數(shù)值模擬強(qiáng)激波能力,提高激波分辨率,進(jìn)一步提高計(jì)算精度。