第25卷第2期刊Vol 25 No. 22004年6月CHINESE QUARTERLY OF MECHANICSJune 2004求解定常勢流正命題的全時空變分有限元法陳池,陶毅,劉高聯(lián)(上海大學力學所,上海200072)摘要:非定常流動變分原理的建立使得用有限元法來求解多工況點的設計問題成為可能。本文在劉高聯(lián)的非定常變分理論的基礎上,對定常變分問題進行時間相關有限元求解。但由于可壓縮非定常位勢流動的控制方程是雙曲型的簡單地把時間當作同空間一樣的物理維來求解是不可行的。而現(xiàn)有的時空有限元法極其復雜增加了計算復雜度,使其很難用于工程設計中。為此,文[2、3]提出了求解一維非定常問題的新型時空有限元法。本文把該方法推廣到二維流動用它求解二維彎管內的流動和翼型繞流問題。計算結果與用定常方法求得的結果幾乎重合,說明該方法可以用于多維時間相關求解關鍵詞:非定常;變分;時空有限元;計算流體力學中圖分類號:0355文獻標識碼:A文章編號:02540053(2004)021576A New Time/Space FEM for Steady Potential FlowCHEN Chi, TAO Yi, LIU Gao-LianShanghai University, Shanghai Inst. App. Math. Mech, Shanghai 200072, China)Abstract: On the basis of Liu Gaolian's unsteady variational theories, the time-marching method is used tosolve steady variational problems. Because the governing equations for compressible unsteady potentialflow is hyperbolic, looking time dimension as space dimension in the same way is never appropriate. Unfortunately, the existing time/ space FEM is quite complicated, and difficult to be put into engineeringuse. A new time/space FEM was presented to calculate the one-dimensional unsteady problems in papers[2,3]. The method is now extended to calculate two-dimensional steady flow in a time-marching way. Asexamples, one pipe flow and one airfoil rounding flow are considered the results agree well with those ofthe steady flows, which demonstrates the usefulness of this method in solving multidimensional problemsKey words: unsteady variational calculus; time/ space FEM; CFD非定常問題是當今研究的熱門問題之一,但到目前為止,多數(shù)非定常問題都是用有限差分或有限體積法求解的。有限元法以其能適應復雜的幾何區(qū)域而逐漸在計算流體力學領域中得到應用。但是,用有限元方法求解非定常問題的文獻也主要是基于 Galerkin變分理論的基礎上進行的),對時間偏導數(shù)的處理大多數(shù)是采用的差分方法。這一方面是因為非定常的變分原理的建立是極其困難的,因為非定常問題的控制方程是雙曲線型的,在時間方向只能給初始條件,而不能給終止條件。劉高聯(lián)根據(jù) Hamilton原理成功地建立了非定常問題的變分原理,為非定常問題的變分求解建立了完整周密的理論基礎。另一方面是因為差分法相對于時空有限元法要容易。但我們注意到,有很多與時間相關的問題是很難用差分法解決的,例如降落傘的降落問題但如果用時空有限元法就容易處理°。變分法以其有強大的變域變分工具,可以為葉輪機進行多工況點反設計而逐漸得到重視。遺憾的是,迄今為止,還沒有發(fā)現(xiàn)非定常變分計收稿日期:2003409-23基金項目:國家自然科學基金項目(50136030,10372055)作者簡介:陳池(1973-),男,四川廣安人,博土力學季刊第25卷算的有效方法。文獻[7對二維振蕩翼型進行了求解,但是文中對時間維的處理方式完全和空間維的處理方式等同。我們知道,非定常絕熱位勢流動的控制方程是雙曲型方程,在時間方向絕對不能采取同空間方向一樣的處理。為此,我們必須采用時空有限元法,但是現(xiàn)有的時空有限元法使用很不方便,公式也很復雜,而且也不適合變分計算。為此,我們需要提出新的求解方法。受雙曲型方程特征的啟發(fā),劉高聯(lián)提出在處理單元內的插值函數(shù)時對時間方向上采用展開處理方式的設想。經過一維的計算試驗,證明是有效的,本文將用它來對二維的絕熱勢流問題進行時間相關求解。1二維非定常勢流的控制方程及其變分理論維非定常勢流其相應的無量綱(以滯止參數(shù)無量綱化)氣動方程為:8t+·()=0(1)中=A(y-1)(atP(3)上述公式中,、A、、P、y分別為無量綱密度速度、位勢壓力和絕熱指數(shù)。有量綱量的定義見文獻[7]計算區(qū)域如圖1所示初邊值條件為(1)在初始邊給定勢函數(shù)∮和密度值P作為初始條件,穩(wěn)定解與初始值關系不大。即φ=f2P=f3。(2)遠場邊界:應用 Riesman不變量來計算。(3)物型面上:元=0(4)激波面上:以g表示其法向分速,則有 Rankine-Hugniot激波關系。+P(A -g)[lA=0|p(An-g)=0圖1計算塊[|A,]=0Fig 1 A slab[|H門]=(y-1)g[|An](5)尾流線上下:an=0(6)尾緣點: Kutta無載荷條件。用變分推導的系統(tǒng)性途徑可以得到以下的變分原理:上述非定常問題的解將使a=0,其中(1-(2∵()+2(1ydadydt+‖(pAn),dtp陽dady(φ-f2)-2∫3∮dxc第2期陳池,等:求解定常勢流正命題的全時空變分有限元法2控制泛函的離散對勢函數(shù)在有限單元內作如下的離散:1-[中k(,n,5=22+a1(1+51+§1-m+a2(1+5]+251274+a2(1+1-1[+a4(1+5](6)引入二維空間單元的形函數(shù)替代上式中含、?的各項,即有(,,3)=N[+a:(1+5](遵守求和約定)其中N=-41=2,M=9+41=2,N2=9+9+2,M,=99+2定義從空間-時間坐標(x,y,t)到有限單元體積坐標(、7之間的 Jacobian轉換矩陣如下atdx ay atasas as從而有單元內勢函數(shù)對空間-時間的偏導數(shù)8沖沖其中J1為 Jacobian矩陣的逆矩陣,J1、J2、J32分別為由 Jacobian矩陣的逆矩陣的第一、二、三列元素組成的行向量。v綽為勢函數(shù)對體積坐標的散度,分別為a285[+a,(1+5]a7a71[+a,(1+5]a=Na,(遵守求和約定),1=12,4(下面對六k同)。由此我們可得到單元內泛函對任意a;的偏導數(shù)為e(1+5e[+a(1+5](1+5)JAAn2[+a1(1+5)/、、、1+5+J;18N,(1+)[+a/(1+5(1+J 2anD[+a1(1+]·J2aNo(1+5)NN力學季刊第25卷+‖(aA,)N(1+5)]·E3dd+‖pN(1+5·E1dln(9)上列諸式中J|為 Jacobian矩陣的行列式,E12、E13為面積分的轉換系數(shù)為了把(9)式中含a的項線性化處理我們把和中的a提出來,從而可以得到單元內寫成矩陣形式的代數(shù)方程組:b2b2 ba b2|a2p2(10)b3 b32 b3 b3abb ba b其中的參數(shù)為(1+)aN(1+y)(1+ya(1+5aN(1+5+J,iaN(1+)aN(1+)PJ3 aN(1+5)+j,aNJ,1(1+5NN(1+)(1+JIdedrds0N+‖(a,)N(1+5]·E2d5d3+‖pN(1+t·E1dηp式中對k有求和約定合并各單元的上述方程可以得到整個網格上的關于a的方程。初始密度假設為來流的密度,初始勢函數(shù)為來流速度與當?shù)貦M坐標的乘積,即:=A2。由于我們的計算方法是隱式的,可以取較大的時間步長。另外,為了消除上述線性方程組的奇異性,我們必須給定某些點上的a值。由此我們以特征線(或面)公式得到周圍邊界r上的速度和密度,從而可以算出當?shù)氐闹?。應用上述變分泛函求?不需要其它邊界條件,因為在泛函中我們已經考慮了所有的邊界條件3數(shù)值結果3.1一維問題求解為了驗證這種方法我們首先對一維非定常的流動進行了求解。我們的例子選自文[5],即要計算如圖2中的大容器當閥門打開后管內的流動狀況。用本文的方法所得的結果和文獻[5]用特征線所求得的結果見圖3和4二者幾乎完全重合,說明我們的方法用于計算一維非定常流動是很成功的3.2二維非定常求解用上述方法我們對二維的可壓縮繞流進行了非定常計算。計算的處理方式同一維問題求解的過程。第2期陳池,等:求解定常勢流正命題的全時空變分有限元法161計算管沒圖2計算管段及相關參數(shù)Fig 2 Calculated pipe and related parameters10000100000圖3壓力比較(用符號衰示的為文[5]的結果圖4速度比較(用符號示的為文[5]的結果Fig 3 Pressure comparison (results denotedFig 4 Velocity comparison(results denotedby character taken from paper [5])by character taken from paper [5])在這里我們做了兩個算例。其中一個算例選自NACA的研究報告,文中是用流函數(shù)和勢函數(shù)計算的定常流動。計算區(qū)域及網格如圖5所示。進口無量綱速度是0.3993。我們的計算經過12個周期達到穩(wěn)定。得到的結果( Time-Marching)和研究報告的結果( Stream-Potential)見圖6和圖7所示。可以看出它們是吻合得相當好的。勢流函數(shù)法“魯一喜一著一著一一看圖5二維管道網格圖6下管壁面速度Fig 6 Velocity on the lower wall of pipe另一個算例是計算的NACA0012翼型的繞流計算區(qū)域和網格見圖8。計算經過40多個時間步(單步步長為10個無量綱時間)計算基本達到穩(wěn)定(勢函數(shù)絕對誤差總和小于105),所得到的結果如圖9所示。為了比較我們用定常的公式計算了同樣的問題把它們的結果畫在了一起?？梢悦黠@看出,本方法所得的結果和定常計算的結果吻合得相當好。實際的計算顯示它們之間勢函數(shù)差值的總和在10以內。力學季刊第25卷特流函數(shù)法圖7上管壁面速度圖8翼型網格Fig. 7 Velocity on the upper wall of pipeFig8 Airfoil grid勢流函數(shù)法圖9壓力系數(shù)Fig 9 Pressure coefficient4結論針對非定常絕熱流動的控制方程的雙曲型特點,本文采用時空變分有限元來計算定常的變分問題。從所得到的結果可以看出,本方法用于計算二維時間相關的定常流動是切實可行的。通過上面的推導過程還可以看出,本方法可以推廣應用于求解多維問題。參考文獻:[11 Liu G L. A vital innovation in Hamilton principle and its extension to initial-value problems[C]. Proc 4th Intl Conf Nonlinear MechShanghai Univ Press, 2002, 90-97.[2]劉高聯(lián)陶毅,李孝偉陳池新型全時空有限元法[C]水動力學學術會議論文集深圳,2003,[3]陶毅劉高聯(lián)李孝偉陳池計算非定常流動的全時空有限元法[C]水動力學學術會議論文集深圳,2003[4 Chung T J. Finite element analysis in fluid dynamics[M]. McGraw-Hill Int Book Company, 1978[5 Hirsch C. Numerical computation of internal external flows, vol 2[M). John-wiley, New York, 1990[6] Hughes T J R, Hulbert G M. Space-time finite element methods for elastodynamics, Formulations and error estimates[J]. Comp Meth[7]郭加宏劉高聯(lián)二維振蕩機翼跨聲速非定常繞流的變域變分有限元解空氣動力學學報[1996,14(4);217-222[8]劉高聯(lián)輛合熱彈性動力學的統(tǒng)一的變分原理族[].力學學報,1999,31(2):165172