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

孫晨揚，李啟良，楊志剛（同濟大學上海地面交通工具風洞中心，上海 201804）

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

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

孫晨揚，李啟良，楊志剛
（同濟大學上海地面交通工具風洞中心，上海 201804）

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

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

0 引言

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

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

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

1 數值方法

1.1 非定常應力形式的N?S方程

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

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

1.2 大渦模擬控制方程及亞格子應力模型

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

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

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

1.3 時間離散形式

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

式中，上標n＋1代表當前時間步的物理量，n代表上一時間步的物理量，Δt代表時間步長，θ是時間離散方式的參數.

參數θ的取值通常有兩種方法：θ＝1對應于全隱式時間方案，具有一階的時間精度，既可用于時間推進的計算，也可用于獲得穩態流動結果；θ＝0.5對應于Crank?Nicolson方案，該方案具有二階的時間精度，常用于計算非定常流動情況，本文采用Crank?Nicolson方案進行時間離散.

1.4 最小二乘有限元解法

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

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

2 方腔流算例驗證

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

圖1 幾何模型Fig.1 Geometrymodel

圖2 網格模型Fig.2 Meshmodel

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

2.1 Re＝3 200

雷諾數Re＝3 200時，流動屬于層流［10］.設定時間步長為0.1 s，每個時間步長內迭代10次，使之達到收斂.整個仿真計算了4 000個時間步長，取后2 000個時間步長用于數據統計.渦量和應力形式的平均值與試驗結果的比較，如圖3所示.

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

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

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

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

2.2 Re＝10 000

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

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

圖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 結論

發展了基于最小二乘有限元法的速度－應力－壓力的非定常形式，并在此基礎上引入大渦模擬方法成功實現湍流的求解.通過方腔流動的計算發現：

1）在層流流動中，傳統速度－渦量－應力形式不如新發展的速度－應力－壓力形式.對于頂蓋驅動方腔流動，應用速度－應力－壓力形式在整個監控測點區域均能獲得與試驗相一致的計算結果；

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

本文的方法能夠實現非定常層流和湍流模擬.應該指出的是，本文方法有待更多算例支撐和評估，從而更加全面比較出最小二乘有限元法的優缺點.

［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

國家自然科學基金青年科學基金（11302153）；同濟大學中央高校基本科研業務費專項資金（20123315）資助項目

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

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