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

姜曉坤 李廷秋

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

0 引 言

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

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

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

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

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

(1)

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

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

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

(2)

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

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

(3)

式中：y為速度入口處縱向坐標；Re=100.

圖1 驗證算例的計算域

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

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

計算值參考值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-最大值、最小值與平均值.

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

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

2.1 一致分裂格式

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

(4)

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

(5)

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

[,,q]/Δt

(6)

(7)

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

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

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

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

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

(8)

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

(9)

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

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

(10)

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

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

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

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

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

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

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

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

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

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

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

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

壓差阻力(M1):

壓差升力(M2):

摩擦阻力(M3):

(11)

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

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

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

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

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

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

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

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

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

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

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

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

4 結(jié) 束 語

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

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