最小二乘有限元法求解非定常應(yīng)力的Navier?Stokes方程

孫晨揚(yáng)，李啟良，楊志剛（同濟(jì)大學(xué)上海地面交通工具風(fēng)洞中心，上海 201804）

文章編號：1001?246X（2015）01?0013?07

最小二乘有限元法求解非定常應(yīng)力的Navier?Stokes方程

孫晨揚(yáng)，李啟良，楊志剛
（同濟(jì)大學(xué)上海地面交通工具風(fēng)洞中心，上海 201804）

為精確求解非定常層流問題，發(fā)展一種非定常速度－應(yīng)力－壓力的方法.采用牛頓法對非線性對流項(xiàng)進(jìn)行線性化處理和預(yù)處理共軛梯度法，實(shí)現(xiàn)了非定常應(yīng)力形式Navier?Stokes方程的求解.方腔層流流動比較發(fā)現(xiàn)，非定常應(yīng)力形式比渦量形式與試驗(yàn)結(jié)果更加吻合，精度更高.該方法有效地解決亞格子應(yīng)力項(xiàng)的問題，實(shí)現(xiàn)基于最小二乘有限元法的湍流求解.比較方腔湍流流動的試驗(yàn)與仿真結(jié)果，證明本文的方法具有可行性，為湍流大渦模擬計(jì)算打下基礎(chǔ).

最小二乘有限元法；非定常速度－應(yīng)力－壓力；大渦模擬；方腔流

0 引言

流體流動均為非定常流動，包括層流和湍流.為有效模擬流體流動，常用的有限差分法和有限體積法均發(fā)展了相應(yīng)的非定常形式［1－3］.最小二乘有限元法作為與之相對應(yīng)的方法，已發(fā)展了非定常的速度－渦量－壓力形式，并在層流流動中有所應(yīng)用［4］.

近年來基于有限體積法的大渦模擬在工程和研究中得到廣泛應(yīng)用.然而，這些計(jì)算結(jié)果均與試驗(yàn)結(jié)果差距較大［5－6］.主要原因是有限體積法的精度限制在二階.而最小二乘有限元法具有通用性強(qiáng)、計(jì)算效率高、收斂性好和計(jì)算準(zhǔn)確并且能夠使用高階精度等特點(diǎn)［4］，成為獲得精確湍流結(jié)果的新方法.

由于傳統(tǒng)的速度－渦量－壓力形式無法處理大渦模擬出現(xiàn)的亞格子應(yīng)力項(xiàng)，本文首先發(fā)展非定常的速度－應(yīng)力－壓力形式，在此基礎(chǔ)上，引入亞格子應(yīng)力，并由此發(fā)展可以求解湍流的方法.通過方腔流動的數(shù)值計(jì)算與試驗(yàn)對比，評估本文方法的準(zhǔn)確性，為最小二乘有限元方法在湍流計(jì)算中的應(yīng)用提供支撐.

1 數(shù)值方法

1.1 非定常應(yīng)力形式的N?S方程

非定常不可壓無量綱應(yīng)力形式的N?S方程，

式中，u為速度，p為壓力，Sij為應(yīng)力張量，f為體積力，Re代表雷諾數(shù).二維問題有6個未知量，三維問題有10個未知量.

1.2 大渦模擬控制方程及亞格子應(yīng)力模型

按照大渦模擬理論［7－8］，經(jīng)濾波函數(shù)處理后的不可壓縮大渦模擬控制方程為

為求解大渦模擬控制方程，引入亞格子應(yīng)力模型，即給出體現(xiàn)小尺度渦對求解方程影響的亞格子應(yīng)力項(xiàng)ij的表達(dá)式.本文采用Smagorinsky亞格子模型

式中，νt是亞格子尺度的湍動粘度，Δi代表沿i軸方向的網(wǎng)格尺寸，Cs稱為Smagorinsky常數(shù)，暫取0.1.

1.3 時間離散形式

時間離散方式對計(jì)算流體力學(xué)的精度和收斂性有較大影響［7－8］.綜合考慮精度和效率，采用經(jīng)典的θ方法.時間離散后如

式中，上標(biāo)n＋1代表當(dāng)前時間步的物理量，n代表上一時間步的物理量，Δt代表時間步長，θ是時間離散方式的參數(shù).

參數(shù)θ的取值通常有兩種方法：θ＝1對應(yīng)于全隱式時間方案，具有一階的時間精度，既可用于時間推進(jìn)的計(jì)算，也可用于獲得穩(wěn)態(tài)流動結(jié)果；θ＝0.5對應(yīng)于Crank?Nicolson方案，該方案具有二階的時間精度，常用于計(jì)算非定常流動情況，本文采用Crank?Nicolson方案進(jìn)行時間離散.

1.4 最小二乘有限元解法

文獻(xiàn)［9］給出定常應(yīng)力形式層流的最小二乘有限元解法的基本步驟，在此不一一列舉.本文的最大差別主要體現(xiàn)在時間項(xiàng)和亞格子應(yīng)力項(xiàng)的處理上，具體表現(xiàn)為剛度矩陣各系數(shù)的不同，如式（5）所示.

結(jié)合文獻(xiàn)［9］中理論，使用C＋＋語言編寫計(jì)算程序，便可以開展基于最小二乘有限元法的流體計(jì)算.計(jì)算程序的關(guān)鍵在于亞格子湍動粘度的計(jì)算和迭代求解算法的編寫.計(jì)算效率由系數(shù)矩陣決定，最小二乘有限元法中系數(shù)矩陣為對稱正定矩陣，可以采用預(yù)處理共軛梯度迭代求解，迭代收斂性很好.

2 方腔流算例驗(yàn)證

方腔流動作為流體力學(xué)中經(jīng)典流動，同時是對數(shù)值算法進(jìn)行驗(yàn)證的經(jīng)典問題.本文采用Prasad［10］等人的試驗(yàn)參數(shù)和結(jié)果.參照該試驗(yàn)參數(shù)建立幾何模型.對于邊界條件，除了頂蓋給定驅(qū)動速度1m·s－1外，其它面均設(shè)為壁面.幾何模型，如圖1所示.由于靠近壁面處流動變化比較劇烈，且考慮到大渦模擬對近壁處的網(wǎng)格要求，采取兩邊密，中間疏的網(wǎng)格形式，網(wǎng)格數(shù)取30×30×30，近壁面第一個網(wǎng)格處y＋約為5，基本滿足大渦模擬計(jì)算要求.劃分網(wǎng)格如圖2所示.

圖1 幾何模型Fig.1 Geometrymodel

圖2 網(wǎng)格模型Fig.2 Meshmodel

根據(jù)Prasad等人的試驗(yàn)結(jié)果，對圖1中正方形對稱面的兩條中線上的速度按時間記錄，并最終比較它們的平均值與均方根值跟試驗(yàn)之間的誤差，下文比較圖中，點(diǎn)表示試驗(yàn)結(jié)果，線表示計(jì)算結(jié)果.

2.1 Re＝3 200

雷諾數(shù)Re＝3 200時，流動屬于層流［10］.設(shè)定時間步長為0.1 s，每個時間步長內(nèi)迭代10次，使之達(dá)到收斂.整個仿真計(jì)算了4 000個時間步長，取后2 000個時間步長用于數(shù)據(jù)統(tǒng)計(jì).渦量和應(yīng)力形式的平均值與試驗(yàn)結(jié)果的比較，如圖3所示.

從平均速度比較結(jié)果可以看出應(yīng)力形式的結(jié)果比渦量形式的結(jié)果與試驗(yàn)值更加接近，具有更高的精確性.具體表現(xiàn)為，應(yīng)力形式的速度平均值與試驗(yàn)值的變化趨勢相同，并且在數(shù)值上也十分接近.除個別點(diǎn)（x/B＝－0.1處，v試驗(yàn)＝0.017，v仿真＝0.029和y/D＝0.1處，u試驗(yàn)＝－0.018，u仿真＝－0.008）外，在－0.45≤x/B≤0.45（B為方腔的寬度）和－0.45≤y/D（D為方腔的深度）≤0.45內(nèi)，試驗(yàn)與仿真的誤差都在15%以下；在靠近壁面附近，速度梯度較大的區(qū)域，應(yīng)力形式的仿真結(jié)果也能與試驗(yàn)結(jié)果保持一致．而渦量形式的仿真結(jié)果在0.3≤｜x/B｜≤0.5（B為方腔的寬度）和0.3≤｜y/D｜（D為方腔的深度）≤0.5均有不同程度的偏差，個別點(diǎn)的偏差超過100%.

圖3 無量綱速度平均值Fig.3 Non?dimensionalmean velocity，（a）vorticity formulation，（b）stress formulation

圖4所示為無量綱速度均方根的試驗(yàn)值與應(yīng)力形式代碼計(jì)算結(jié)果的比較.同樣能夠發(fā)現(xiàn)其與試驗(yàn)結(jié)果吻合的較好，仿真結(jié)果與試驗(yàn)結(jié)果的變化趨勢相同，在數(shù)值上也比較接近，誤差在大部分點(diǎn)都在30%以下，對于無量綱速度均方根來說，這一結(jié)果是可以接受的.渦量形式得到的無量綱速度均方根在所有測點(diǎn)均小于0.01，與試驗(yàn)結(jié)果差距較大.頂蓋驅(qū)動的方腔流動是由頂蓋以剪切力的形式帶動整個方腔流體流動.它是典型的剪切力驅(qū)動的流動.正因如此，與渦量形式相比，應(yīng)力形式由于能夠在邊界和整個計(jì)算區(qū)域直接給出或計(jì)算出剪切應(yīng)力，因而具有更高精度.渦量形式則需要獲得速度梯度后才能計(jì)算出剪切應(yīng)力，受網(wǎng)格的影響，必然會有相應(yīng)的數(shù)值誤差.這些誤差在速度梯度較大的地方更為明顯，從而導(dǎo)致在上述速度變化劇烈的地方，與試驗(yàn)差距較大.

圖4 無量綱速度均方根（Urms＝10 U－′2/U0，Vrms＝10 V－′2/U0）Fig.4 Non?dimensional root?mean square velocity

2.2 Re＝10 000

當(dāng)雷諾數(shù)Re增大到10 000時，此時方腔流動進(jìn)入到湍流狀態(tài)［10］.湍流耗散起著重要的作用.為了有效模擬方腔的湍流流動，本文發(fā)展了基于最小二乘有限元的大渦模擬方法.應(yīng)該指出的是，整個計(jì)算設(shè)置均與上述設(shè)置相同.圖5和圖6分別給出了無量綱速度平均值和均方根值.

從圖中可以看出，雷諾數(shù)增大改變了無量綱速度平均值的變化規(guī)律.首先它使靠近壁面的區(qū)域速度變化更加陡峭，更加符合湍流邊界層內(nèi)速度分布規(guī)律.另外，對于x/B處，它的最大值有所減少，最小值有所增加；對于y/D，它的最小值有所增加.與層流相比，湍流粘性的作用使速度脈動更加強(qiáng)烈，表現(xiàn)出速度均方根數(shù)值更大，變化更明顯.觀察圖6的數(shù)值曲線可以發(fā)現(xiàn)，基于最小二乘有限元的大渦模擬計(jì)算結(jié)果具有可信性；各測點(diǎn)的數(shù)值與試驗(yàn)值吻合較好.

圖5 無量綱速度平均值Fig.5 Non?dimensionalmean velocity

圖6 無量綱速度均方根（Urms＝10 U－′2/U0，Vrms＝10 V－′2/U0）Fig.6 Non?dimensional root?mean square velocity

3 結(jié)論

發(fā)展了基于最小二乘有限元法的速度－應(yīng)力－壓力的非定常形式，并在此基礎(chǔ)上引入大渦模擬方法成功實(shí)現(xiàn)湍流的求解.通過方腔流動的計(jì)算發(fā)現(xiàn)：

1）在層流流動中，傳統(tǒng)速度－渦量－應(yīng)力形式不如新發(fā)展的速度－應(yīng)力－壓力形式.對于頂蓋驅(qū)動方腔流動，應(yīng)用速度－應(yīng)力－壓力形式在整個監(jiān)控測點(diǎn)區(qū)域均能獲得與試驗(yàn)相一致的計(jì)算結(jié)果；

2）對于湍流流動，通過引入大渦模擬實(shí)現(xiàn)了基于最小二乘有限元法的湍流求解.對于Re＝10 000的方腔流動，本文方法成功模擬出方腔的湍流流動，取得與試驗(yàn)吻合較好的結(jié)果.

本文的方法能夠?qū)崿F(xiàn)非定常層流和湍流模擬.應(yīng)該指出的是，本文方法有待更多算例支撐和評估，從而更加全面比較出最小二乘有限元法的優(yōu)缺點(diǎn).

［1］ HuangW Y，Wei JY，Ge Y B.Finite differencemethod for 3D unsteady incompressible vorticity?velocity formulation Navier?Stokes equations［J］.Journal on Numerical Methods and Computer Applications，2012，4：301－311.

［2］ An J，Sun P，Luo ZD，Huang X M.Stablized discrete finite volumemethod for unsteady Stokes equations［J］.Computational Mathematics，2011，2：213－224.

［3］ Wang F J.Analysis of computational fluid dynamics［M］.Beijing：Tsinghua University Press，2010.

［4］ Tao S.Application of the least?squares finite element method for computational fluid dynamics［D］.Shanghai：Tongji University，2012.

［5］ Xie Z G，Xu C X，Cui G X，Wang Z S.Large eddy simulation of flows around a square cylinder［J］.Chinese Journal of Computational Physics，2007，24（2）：171－180.

［6］ Jia X H，Liu H.Large eddy simulation of flow around two circular cylinders［J］.Chinese Journal of Hydrodynamics Edition：A，2008，6：625－632.

［7］ Gao H，Zhou X J.Implicit iteration time?advancement scheme for compressible Navier?Stokes equations［J］.Chinese Journal of Computational Physics，2008，25（1）：51－57.

［8］ Cheng J，Shu CW.High order schemes for CFD：A Review［J］.Chinese Journal of Computational Physics，2009，26（5）：633－655.

［9］ Smagorinsky J.General circulation experiments with the primitive equations：Part I.The basic experiment［J］.Monthly Weather Review，1963，91：99－164.

［10］ Ghosal S，Moin P.The basic equations for the large eddy simulation of turbulent flows in complex geometry［J］.Journal of Computational Physics，1995，118：24－37.

［11］ Tao S，Yang ZG，Gu W J，et al.Application comparison of least square finite elementmethod and finite volumemethod in CFD［J］.Computer Aided Engineering，2012，21（2）：1－6.

［12］ Prasad A K，Koseff JR.Reynolds number and end?wall effects on a lid?driven cavity flow［J］.Physics of Fluids，1989，A1：208－218.

Least?squares Finite Element M ethod for Unsteady Stress Formulation of Navier?Stokes Equations

SUN Chenyang，LIQiliang，YANG Zhigang （Shanghai AutomobileWind Tunnel Center，Tongji University，201804 Shanghai，China）

To solve unsteady laminar flow problems，a method of velocity?stress?pressure formulation instead of velocity?vorticity?pressure formulation is developed.With Newton's linearizedmethod to linearize convective terms and preconditioned conjugate gradient method to solve equations，unsteady stress formulation of Navier?Stokes equations is solved.Comparison between numerical and experimental results of cavity laminar flow shows that result of stress formulation fits experiment better and has higher accuracy than vorticity formulation.The stress formulation can deal with subgrid stress with least squares finite elementmethod.Comparison with experimental results of cavity turbulent flow reveals feasibility of the method.It lays a firm foundation for large eddy simulation computation.

least squares finite elementmethod；unsteady velocity?stress?pressure；large eddy simulation；cavity flow

TP301.6

A

2013－12－25；

2014－04－01

國家自然科學(xué)基金青年科學(xué)基金（11302153）；同濟(jì)大學(xué)中央高?；究蒲袠I(yè)務(wù)費(fèi)專項(xiàng)資金（20123315）資助項(xiàng)目

孫晨揚(yáng)（1991－），男，江蘇興化人，碩士生，主要研究車身與空氣動力學(xué)，E?mail：sun?chenyang＠163.com

Received date： 2013－12－25；Revised date： 2014－04－01