第25卷第2期刊Vol 25 No. 22004年6月CHINESE QUARTERLY OF MECHANICSJune 2004求解定常勢(shì)流正命題的全時(shí)空變分有限元法陳池,陶毅,劉高聯(lián)(上海大學(xué)力學(xué)所,上海200072)摘要:非定常流動(dòng)變分原理的建立使得用有限元法來求解多工況點(diǎn)的設(shè)計(jì)問題成為可能。本文在劉高聯(lián)的非定常變分理論的基礎(chǔ)上,對(duì)定常變分問題進(jìn)行時(shí)間相關(guān)有限元求解。但由于可壓縮非定常位勢(shì)流動(dòng)的控制方程是雙曲型的簡單地把時(shí)間當(dāng)作同空間一樣的物理維來求解是不可行的。而現(xiàn)有的時(shí)空有限元法極其復(fù)雜增加了計(jì)算復(fù)雜度,使其很難用于工程設(shè)計(jì)中。為此,文[2、3]提出了求解一維非定常問題的新型時(shí)空有限元法。本文把該方法推廣到二維流動(dòng)用它求解二維彎管內(nèi)的流動(dòng)和翼型繞流問題。計(jì)算結(jié)果與用定常方法求得的結(jié)果幾乎重合,說明該方法可以用于多維時(shí)間相關(guān)求解關(guān)鍵詞:非定常;變分;時(shí)空有限元;計(jì)算流體力學(xué)中圖分類號(hào):0355文獻(xiàn)標(biāo)識(shí)碼:A文章編號(hào):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非定常問題是當(dāng)今研究的熱門問題之一,但到目前為止,多數(shù)非定常問題都是用有限差分或有限體積法求解的。有限元法以其能適應(yīng)復(fù)雜的幾何區(qū)域而逐漸在計(jì)算流體力學(xué)領(lǐng)域中得到應(yīng)用。但是,用有限元方法求解非定常問題的文獻(xiàn)也主要是基于 Galerkin變分理論的基礎(chǔ)上進(jìn)行的),對(duì)時(shí)間偏導(dǎo)數(shù)的處理大多數(shù)是采用的差分方法。這一方面是因?yàn)榉嵌ǔ５淖兎衷淼慕⑹菢O其困難的,因?yàn)榉嵌ǔ栴}的控制方程是雙曲線型的,在時(shí)間方向只能給初始條件,而不能給終止條件。劉高聯(lián)根據(jù) Hamilton原理成功地建立了非定常問題的變分原理,為非定常問題的變分求解建立了完整周密的理論基礎(chǔ)。另一方面是因?yàn)椴罘址ㄏ鄬?duì)于時(shí)空有限元法要容易。但我們注意到,有很多與時(shí)間相關(guān)的問題是很難用差分法解決的,例如降落傘的降落問題但如果用時(shí)空有限元法就容易處理°。變分法以其有強(qiáng)大的變域變分工具,可以為葉輪機(jī)進(jìn)行多工況點(diǎn)反設(shè)計(jì)而逐漸得到重視。遺憾的是,迄今為止,還沒有發(fā)現(xiàn)非定常變分計(jì)收稿日期:2003409-23基金項(xiàng)目:國家自然科學(xué)基金項(xiàng)目(50136030,10372055)作者簡介:陳池(1973-),男,四川廣安人,博土力學(xué)季刊第25卷算的有效方法。文獻(xiàn)[7對(duì)二維振蕩翼型進(jìn)行了求解,但是文中對(duì)時(shí)間維的處理方式完全和空間維的處理方式等同。我們知道,非定常絕熱位勢(shì)流動(dòng)的控制方程是雙曲型方程,在時(shí)間方向絕對(duì)不能采取同空間方向一樣的處理。為此,我們必須采用時(shí)空有限元法,但是現(xiàn)有的時(shí)空有限元法使用很不方便,公式也很復(fù)雜,而且也不適合變分計(jì)算。為此,我們需要提出新的求解方法。受雙曲型方程特征的啟發(fā),劉高聯(lián)提出在處理單元內(nèi)的插值函數(shù)時(shí)對(duì)時(shí)間方向上采用展開處理方式的設(shè)想。經(jīng)過一維的計(jì)算試驗(yàn),證明是有效的,本文將用它來對(duì)二維的絕熱勢(shì)流問題進(jìn)行時(shí)間相關(guān)求解。1二維非定常勢(shì)流的控制方程及其變分理論維非定常勢(shì)流其相應(yīng)的無量綱(以滯止參數(shù)無量綱化)氣動(dòng)方程為:8t+·()=0(1)中=A(y-1)(atP(3)上述公式中,、A、、P、y分別為無量綱密度速度、位勢(shì)壓力和絕熱指數(shù)。有量綱量的定義見文獻(xiàn)[7]計(jì)算區(qū)域如圖1所示初邊值條件為(1)在初始邊給定勢(shì)函數(shù)∮和密度值P作為初始條件,穩(wěn)定解與初始值關(guān)系不大。即φ=f2P=f3。(2)遠(yuǎn)場邊界:應(yīng)用 Riesman不變量來計(jì)算。(3)物型面上:元=0(4)激波面上:以g表示其法向分速,則有 Rankine-Hugniot激波關(guān)系。+P(A -g)[lA=0|p(An-g)=0圖1計(jì)算塊[|A,]=0Fig 1 A slab[|H門]=(y-1)g[|An](5)尾流線上下:an=0(6)尾緣點(diǎn): Kutta無載荷條件。用變分推導(dǎo)的系統(tǒng)性途徑可以得到以下的變分原理:上述非定常問題的解將使a=0,其中(1-(2∵()+2(1ydadydt+‖(pAn),dtp陽dady(φ-f2)-2∫3∮dxc第2期陳池,等:求解定常勢(shì)流正命題的全時(shí)空變分有限元法2控制泛函的離散對(duì)勢(shì)函數(shù)在有限單元內(nèi)作如下的離散:1-[中k(,n,5=22+a1(1+51+§1-m+a2(1+5]+251274+a2(1+1-1[+a4(1+5](6)引入二維空間單元的形函數(shù)替代上式中含、?的各項(xiàng),即有(,,3)=N[+a:(1+5](遵守求和約定)其中N=-41=2,M=9+41=2,N2=9+9+2,M,=99+2定義從空間-時(shí)間坐標(biāo)(x,y,t)到有限單元體積坐標(biāo)(、7之間的 Jacobian轉(zhuǎn)換矩陣如下atdx ay atasas as從而有單元內(nèi)勢(shì)函數(shù)對(duì)空間-時(shí)間的偏導(dǎo)數(shù)8沖沖其中J1為 Jacobian矩陣的逆矩陣,J1、J2、J32分別為由 Jacobian矩陣的逆矩陣的第一、二、三列元素組成的行向量。v綽為勢(shì)函數(shù)對(duì)體積坐標(biāo)的散度,分別為a285[+a,(1+5]a7a71[+a,(1+5]a=Na,(遵守求和約定),1=12,4(下面對(duì)六k同)。由此我們可得到單元內(nèi)泛函對(duì)任意a;的偏導(dǎo)數(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力學(xué)季刊第25卷+‖(aA,)N(1+5)]·E3dd+‖pN(1+5·E1dln(9)上列諸式中J|為 Jacobian矩陣的行列式,E12、E13為面積分的轉(zhuǎn)換系數(shù)為了把(9)式中含a的項(xiàng)線性化處理我們把和中的a提出來,從而可以得到單元內(nèi)寫成矩陣形式的代數(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式中對(duì)k有求和約定合并各單元的上述方程可以得到整個(gè)網(wǎng)格上的關(guān)于a的方程。初始密度假設(shè)為來流的密度,初始勢(shì)函數(shù)為來流速度與當(dāng)?shù)貦M坐標(biāo)的乘積,即:=A2。由于我們的計(jì)算方法是隱式的,可以取較大的時(shí)間步長。另外,為了消除上述線性方程組的奇異性,我們必須給定某些點(diǎn)上的a值。由此我們以特征線(或面)公式得到周圍邊界r上的速度和密度,從而可以算出當(dāng)?shù)氐闹怠?yīng)用上述變分泛函求解,不需要其它邊界條件,因?yàn)樵诜汉形覀円呀?jīng)考慮了所有的邊界條件3數(shù)值結(jié)果3.1一維問題求解為了驗(yàn)證這種方法我們首先對(duì)一維非定常的流動(dòng)進(jìn)行了求解。我們的例子選自文[5],即要計(jì)算如圖2中的大容器當(dāng)閥門打開后管內(nèi)的流動(dòng)狀況。用本文的方法所得的結(jié)果和文獻(xiàn)[5]用特征線所求得的結(jié)果見圖3和4二者幾乎完全重合,說明我們的方法用于計(jì)算一維非定常流動(dòng)是很成功的3.2二維非定常求解用上述方法我們對(duì)二維的可壓縮繞流進(jìn)行了非定常計(jì)算。計(jì)算的處理方式同一維問題求解的過程。第2期陳池,等:求解定常勢(shì)流正命題的全時(shí)空變分有限元法161計(jì)算管沒圖2計(jì)算管段及相關(guān)參數(shù)Fig 2 Calculated pipe and related parameters10000100000圖3壓力比較(用符號(hào)衰示的為文[5]的結(jié)果圖4速度比較(用符號(hào)示的為文[5]的結(jié)果Fig 3 Pressure comparison (results denotedFig 4 Velocity comparison(results denotedby character taken from paper [5])by character taken from paper [5])在這里我們做了兩個(gè)算例。其中一個(gè)算例選自NACA的研究報(bào)告,文中是用流函數(shù)和勢(shì)函數(shù)計(jì)算的定常流動(dòng)。計(jì)算區(qū)域及網(wǎng)格如圖5所示。進(jìn)口無量綱速度是0.3993。我們的計(jì)算經(jīng)過12個(gè)周期達(dá)到穩(wěn)定。得到的結(jié)果( Time-Marching)和研究報(bào)告的結(jié)果( Stream-Potential)見圖6和圖7所示?？梢钥闯鏊鼈兪俏呛系孟喈?dāng)好的。勢(shì)流函數(shù)法“魯一喜一著一著一一看圖5二維管道網(wǎng)格圖6下管壁面速度Fig 6 Velocity on the lower wall of pipe另一個(gè)算例是計(jì)算的NACA0012翼型的繞流計(jì)算區(qū)域和網(wǎng)格見圖8。計(jì)算經(jīng)過40多個(gè)時(shí)間步(單步步長為10個(gè)無量綱時(shí)間)計(jì)算基本達(dá)到穩(wěn)定(勢(shì)函數(shù)絕對(duì)誤差總和小于105),所得到的結(jié)果如圖9所示。為了比較我們用定常的公式計(jì)算了同樣的問題把它們的結(jié)果畫在了一起?？梢悦黠@看出,本方法所得的結(jié)果和定常計(jì)算的結(jié)果吻合得相當(dāng)好。實(shí)際的計(jì)算顯示它們之間勢(shì)函數(shù)差值的總和在10以內(nèi)。力學(xué)季刊第25卷特流函數(shù)法圖7上管壁面速度圖8翼型網(wǎng)格Fig. 7 Velocity on the upper wall of pipeFig8 Airfoil grid勢(shì)流函數(shù)法圖9壓力系數(shù)Fig 9 Pressure coefficient4結(jié)論針對(duì)非定常絕熱流動(dòng)的控制方程的雙曲型特點(diǎn),本文采用時(shí)空變分有限元來計(jì)算定常的變分問題。從所得到的結(jié)果可以看出,本方法用于計(jì)算二維時(shí)間相關(guān)的定常流動(dòng)是切實(shí)可行的。通過上面的推導(dǎo)過程還可以看出,本方法可以推廣應(yīng)用于求解多維問題。參考文獻(xiàn):[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)陶毅,李孝偉陳池新型全時(shí)空有限元法[C]水動(dòng)力學(xué)學(xué)術(shù)會(huì)議論文集深圳,2003,[3]陶毅劉高聯(lián)李孝偉陳池計(jì)算非定常流動(dòng)的全時(shí)空有限元法[C]水動(dòng)力學(xué)學(xué)術(shù)會(huì)議論文集深圳,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)二維振蕩機(jī)翼跨聲速非定常繞流的變域變分有限元解空氣動(dòng)力學(xué)學(xué)報(bào)[1996,14(4);217-222[8]劉高聯(lián)輛合熱彈性動(dòng)力學(xué)的統(tǒng)一的變分原理族[].力學(xué)學(xué)報(bào),1999,31(2):165172