不可壓縮流動(dòng)的非結(jié)構(gòu)網(wǎng)格自適應(yīng)方法研究*

姜曉坤 李廷秋

(武漢理工大學(xué)交通學(xué)院 武漢 430063)

0 引 言

自適應(yīng)策略設(shè)計(jì)的基本思想是基于現(xiàn)有數(shù)據(jù)對(duì)數(shù)值解做出誤差估計(jì)，根據(jù)得到的誤差指示器調(diào)整網(wǎng)格單元或基函數(shù)階數(shù).調(diào)整網(wǎng)格單元密度的自適應(yīng)方法被稱為h自適應(yīng)方法，保持了基函數(shù)的階不變，得到代數(shù)的收斂階，具有較快的收斂性. 與彈性連續(xù)介質(zhì)中的有限元方法不同，不可壓縮流動(dòng)的數(shù)值計(jì)算需要處理壓力速度耦合問(wèn)題. 分裂格式可以避免在計(jì)算對(duì)流-擴(kuò)散，流體流動(dòng)等問(wèn)題中求解耦合方程.Selim等[1]討論了基于增量壓力修正分裂格式 (incremental pressure correction scheme, IPCS)的自適應(yīng)有限元法求解不可壓縮流動(dòng)，觀察到時(shí)間步長(zhǎng)與動(dòng)量方程計(jì)算誤差呈線性關(guān)系，與連續(xù)性方程的計(jì)算誤差呈平方關(guān)系. Shen等[2]利用離散微分算子推導(dǎo)了一致分裂格式在流場(chǎng)計(jì)算中壓力-速度的最優(yōu)誤差估計(jì)，但將一致分裂格式應(yīng)用在自適應(yīng)方法還未見研究. 不同于IPCS這類基于壓力/速度修正的分裂格式，一致分裂格式?jīng)]有分裂誤差，在渦量和壓力求解上可完全準(zhǔn)確求解. 一致分裂格式具有無(wú)條件穩(wěn)定的特點(diǎn)，較能滿足海洋工程復(fù)雜流場(chǎng)計(jì)算中高精度、高穩(wěn)定性的要求.

文中基于一致分裂格式，針對(duì)物體邊界網(wǎng)格密度不足的問(wèn)題，開發(fā)了基于誤差指示器的非結(jié)構(gòu)網(wǎng)格自動(dòng)加密模塊，應(yīng)用在有限元流場(chǎng)求解程序上.討論了網(wǎng)格加密比率、最大加密迭代次數(shù)以及初始網(wǎng)格密度在圓柱繞流問(wèn)題上的誤差分析. 依照不同的受力構(gòu)建了誤差指示器. 最后討論了方形體繞流計(jì)算中不同分裂格式的計(jì)算效率與計(jì)算精度情況.

1 求解器數(shù)值驗(yàn)證

計(jì)算誤差η由空間離散誤差ηh，時(shí)間離散誤差ηc及數(shù)值誤差ηk組成.

|η|=|ηh+ηk+ηc|

(1)

文中主要討論通過(guò)加密網(wǎng)格減小空間離散誤差ηh的手段提高計(jì)算精度，首先驗(yàn)證一致分裂格式的求解器在求解繞流問(wèn)題上具有一定的數(shù)值精度，即在保證了數(shù)值精度的條件下，可以將ηk近似為0.

待求解的N-S方程與連續(xù)性方程可簡(jiǎn)要表達(dá)為

ut-νΔu+uu+p=0,▽·u=0

(2)

式中:u為速度矢量;ν為粘性系數(shù);p為壓力.

基于文獻(xiàn)[3]中的DFG 2D-2基準(zhǔn)測(cè)試驗(yàn)證有限元求解器在非結(jié)構(gòu)化網(wǎng)格下的數(shù)值精度. 計(jì)算域見圖1，來(lái)流速度定義為

(3)

式中：y為速度入口處縱向坐標(biāo)；Re=100.

圖1 驗(yàn)證算例的計(jì)算域

圖2 驗(yàn)證算例t=8 s時(shí)刻圓柱附近速度云圖與網(wǎng)格劃分

圖3 驗(yàn)證算例圓柱受力系數(shù)歷程曲線

計(jì)算值參考值Cd max/min3.245 7/2.954 83.212 2/3.133 3Cd mean3.162 23.164 7Cl max/min1.052 2/-1.084 50.9839/-1.018 2Cl mean-0.002 1-0.0172Sr0.297 70.3030

注：Cvd-阻力系數(shù)；Cl-升力系數(shù);max,min,mean-最大值、最小值與平均值.

驗(yàn)證算例表明，基于一致分裂格式的流場(chǎng)求解器有較好的數(shù)值精度與穩(wěn)定性.

2 基于一致分裂格式的自適應(yīng)方法

2.1 一致分裂格式

一致分裂格式由Guermond等[4]提出，其基本思想是基于空間的連續(xù)性推導(dǎo)出壓力變分式，并結(jié)合不可壓縮連續(xù)性方程求解壓力的近似值，然后對(duì)壓力進(jìn)行插值并帶入動(dòng)量方程更新速度. 一致分裂格式的流程如下.

(4)

步驟2將壓力近似值帶入動(dòng)量方程的變分形式求解當(dāng)前速度.

(5)

式中：v為速度場(chǎng)至希爾伯特空間的投影向量；σ為科氏應(yīng)力張量；ε為速度的對(duì)稱梯度，定義為

[,,q]/Δt

(6)

(7)

與速度/壓力修正算法相比，一致分裂格式在壓力上沒有分裂誤差. 一致分裂格式在壓力收斂上速度更快，可以達(dá)到二階的收斂速度，而增量壓力修正方法壓力收斂速度僅為一階。

2.2 誤差指示器與自適應(yīng)算法

誤差先驗(yàn)估計(jì)是一個(gè)對(duì)偶問(wèn)題，相關(guān)的數(shù)學(xué)推導(dǎo)可參考文獻(xiàn)[5]，本文僅對(duì)誤差表達(dá)與指示器做一般性闡述.

首先，根據(jù)伴隨方程的解φ得到誤差目標(biāo)函數(shù)M為

M(e)=r(e,φ)=r(u,φ)-r(U,φ)=-r(U,φ)

(8)

式中:誤差e可表示為e=u-U，u為對(duì)偶問(wèn)題中的測(cè)試函數(shù)(test function)解，U為變分方程近似解. 將目標(biāo)函數(shù)表達(dá)為空間離散形式：

(9)

式中：k為每次迭代步中需要加密的網(wǎng)格比例. 需要加密的網(wǎng)格單元誤差指示器為εk=|r(U,φ)k|. 為了防止自適應(yīng)全局加密時(shí)指示器r(U,φ)=0而局部的誤差指示器不為0的情況，將誤差指示器修改為

εk=|[r(U,φ),Z-πhZ]k|

(10)

式中：πh為矢量速度場(chǎng)與壓力在半離散空間上的插值.

自適應(yīng)算法流程為：

步驟1選擇初始的網(wǎng)格作為起始猜測(cè)值.

步驟2基于一致分裂格式求解t=j時(shí)刻的速度與壓力.

步驟3利用求解弱對(duì)偶問(wèn)題的殘差絕對(duì)值與速度/壓力之間進(jìn)行插值，求得每一個(gè)單元的誤差指示器εk.

步驟4使用Rivara加密方法[6]對(duì)網(wǎng)格比例k的網(wǎng)格進(jìn)行加密，當(dāng)小于最大迭代次數(shù)或不滿足收斂準(zhǔn)則時(shí)，返回步驟3.

在構(gòu)建誤差指示器時(shí)，得到插值πh時(shí)需要分析耦合方程的誤差. 使用一致分裂格式無(wú)需討論分裂誤差項(xiàng)，誤差指示器εk計(jì)算更加簡(jiǎn)潔.

3 自適應(yīng)方法數(shù)值分析

3.1 圓柱體繞流問(wèn)題的自適應(yīng)

根據(jù)以上算法，采用與驗(yàn)證算例相同的計(jì)算域，在流體完全流入計(jì)算域的時(shí)刻，開始網(wǎng)格迭代自動(dòng)加密. 圓柱繞流在迭代的非結(jié)構(gòu)網(wǎng)格自適應(yīng)示意圖見圖4.使用全局加密的方法，可以觀察到自適應(yīng)模塊在圓柱附近自動(dòng)加密了網(wǎng)格. 與手工加密網(wǎng)格的結(jié)果(見圖2)相比，自適應(yīng)網(wǎng)格具有一定的非對(duì)稱性，且加密區(qū)域集中在圓柱分離點(diǎn)處. 該示意圖說(shuō)明非結(jié)構(gòu)網(wǎng)格自適應(yīng)模塊很好的處理了初始網(wǎng)格過(guò)于稀疏的問(wèn)題，并且在圓柱尾流處加密幅度較小，較能符合區(qū)域自適應(yīng)的要求.

圖4 三個(gè)加密迭代步下的圓柱繞流網(wǎng)格自適應(yīng)

在自適應(yīng)誤差估計(jì)時(shí)，需要誤差目標(biāo)函數(shù)，分別定義為

壓差阻力(M1):

壓差升力(M2):

摩擦阻力(M3):

(11)

文獻(xiàn)[7]對(duì)圓柱后方的尾渦進(jìn)行了自動(dòng)加密，但在海洋工程實(shí)際應(yīng)用中，更關(guān)心的是圓柱體的受力尤其是升力系數(shù)的計(jì)算精度，基于圓柱壓力構(gòu)造的目標(biāo)函數(shù)更符合工程實(shí)際需求.

三種誤差目標(biāo)函數(shù)下網(wǎng)格加密比率k值以及最大迭代次數(shù)對(duì)計(jì)算誤差的影響見圖5，其中，refined表示初始網(wǎng)格加密后的結(jié)果. 可以觀察到：①網(wǎng)格加密比率k對(duì)計(jì)算誤差沒有顯著影響，當(dāng)最大迭代次數(shù)達(dá)到一定次數(shù)后，計(jì)算精度有逐漸收斂的過(guò)程；②應(yīng)用不同受力構(gòu)造的誤差目標(biāo)函數(shù)，其計(jì)算誤差有較大差異，且摩擦阻力構(gòu)造的誤差函數(shù)的網(wǎng)格加密對(duì)圓柱繞流數(shù)值精度有顯著提高.

圖5 圓柱繞流不同誤差目標(biāo)函數(shù)隨網(wǎng)格加密比率及最大迭代次數(shù)計(jì)算精度曲線

3.2 方型體繞流問(wèn)題的自適應(yīng)

方形體繞流算例主要討論尖銳界面下的自適應(yīng)策略及自適應(yīng)方法在非對(duì)稱繞流問(wèn)題下的計(jì)算情況，計(jì)算域設(shè)置與文獻(xiàn)[1]一致.

基于均勻的非結(jié)構(gòu)網(wǎng)格，方型體繞流在五個(gè)迭代步下的自適應(yīng)見圖6，可以觀察到尖銳點(diǎn)處的網(wǎng)格得到了很好加密，對(duì)梯度變化較慢的區(qū)域如角落處的網(wǎng)格加密較小，同時(shí)自適應(yīng)模塊在方型體頂部后方流速較快的區(qū)域也加密了網(wǎng)格，對(duì)應(yīng)的速度云圖也得到了一定的光順.

圖6 方型體繞流迭代自適應(yīng)的速度云圖與網(wǎng)格變化

不同誤差目標(biāo)函數(shù)隨網(wǎng)格變化比率k的變化曲線見圖7，和圖5a)不同，壓差升力的計(jì)算精度明顯高于摩擦阻力和壓差阻力，且隨著網(wǎng)格加密比率的增大，計(jì)算精度反而有減小的情況. 貼壁的方型體繞流與圓柱繞流問(wèn)題不同在于：流體在迎流面會(huì)向上流動(dòng)，作用于物體一個(gè)向上的升力，而將阻力作為誤差目標(biāo)函數(shù)并不能反映出這一受力. 圓柱繞流的升力一直在零值處作周期性振蕩，網(wǎng)格自適應(yīng)如果對(duì)頻率較快的波動(dòng)流場(chǎng)進(jìn)行網(wǎng)格調(diào)整，計(jì)算耗時(shí)較長(zhǎng).

圖7 圓柱繞流不同誤差目標(biāo)函數(shù)隨網(wǎng)格加密比率計(jì)算精度曲線

Chorin投影法、增量壓力修正，以及一致分裂格式隨網(wǎng)格數(shù)的增加其計(jì)算時(shí)間與計(jì)算精度曲線見圖8.

圖8 方體繞流不同分裂格式的計(jì)算時(shí)間及計(jì)算誤差隨網(wǎng)格數(shù)變化曲線

由圖8可知，一致分裂格式在網(wǎng)格數(shù)量增加后相比與經(jīng)典的分裂格式在計(jì)算效率和計(jì)算精度上均有一定的優(yōu)勢(shì)，但差異不大. 一致分裂格式在網(wǎng)格數(shù)量較小時(shí)計(jì)算精度的波動(dòng)性較小，隨著網(wǎng)格數(shù)量的增大，其計(jì)算效率有顯著的提升.

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

本文主要對(duì)繞流問(wèn)題下非結(jié)構(gòu)網(wǎng)格自適應(yīng)策略進(jìn)行了討論，得出結(jié)論有：最大加密迭代次數(shù)、初始網(wǎng)格密度以及網(wǎng)格加密比率對(duì)計(jì)算誤差并無(wú)顯著影響；使用不同的受力特征構(gòu)造的誤差目標(biāo)函數(shù)計(jì)算誤差影響較大，一致分裂格式與增量壓力修正格式與投影法相比，在計(jì)算速度上有一定優(yōu)勢(shì).應(yīng)用到自適應(yīng)方法時(shí)，一致分裂格式無(wú)誤差項(xiàng)，穩(wěn)定性較好.

關(guān)于自適應(yīng)算法在計(jì)算流體力學(xué)的應(yīng)用值得進(jìn)一步研究的問(wèn)題有:①瞬態(tài)問(wèn)題的先驗(yàn)誤差估計(jì),即時(shí)間步長(zhǎng)的自適應(yīng)調(diào)整[9]；②一致分裂格式對(duì)周期性邊界條件的適應(yīng)性較差, 而其改進(jìn)形式或其他現(xiàn)代分裂格式在自適應(yīng)方法中的應(yīng)用需要進(jìn)一步探討[10]；③運(yùn)動(dòng)物體以及多體問(wèn)題的網(wǎng)格自適應(yīng)問(wèn)題，其核心在于網(wǎng)格自動(dòng)生成的方式，基于笛卡爾網(wǎng)格以及切割網(wǎng)格法的h塊自適應(yīng)方法[11]在這一領(lǐng)域具有良好的研究?jī)r(jià)值.