基于運(yùn)動(dòng)嵌套網(wǎng)格的槳葉cfd網(wǎng)格生成

基于運(yùn)動(dòng)嵌套網(wǎng)格的槳葉cfd網(wǎng)格生成

翼價(jià)計(jì)算水流是直升機(jī)空氣動(dòng)力學(xué)的一個(gè)重要研究方向。網(wǎng)格生成技術(shù)在旋翼CFD研究中占有重要的地位。導(dǎo)致目前旋翼CFD發(fā)展比固定翼滯后的原因之一就是在旋翼周圍生成貼體且正交的網(wǎng)格比在固定翼上復(fù)雜得多。這是因?yàn)闃?gòu)成旋翼系統(tǒng)的每一片槳葉與機(jī)翼相比較,不是作單純的直線運(yùn)動(dòng),而是作旋轉(zhuǎn)運(yùn)動(dòng)。無(wú)論是鉸接式旋翼還是無(wú)鉸接式旋翼,槳葉除旋轉(zhuǎn)運(yùn)動(dòng)外,同時(shí)還存在變距、揮舞和擺振運(yùn)動(dòng)。槳葉的揮舞、變距運(yùn)動(dòng)是旋翼區(qū)別于固定翼的重要特征之一,它們對(duì)旋翼的空氣動(dòng)力學(xué)和動(dòng)力學(xué)特性均有著十分重要的影響。若考慮槳葉的彈性,在上述運(yùn)動(dòng)的基礎(chǔ)上還必須疊加槳葉的彈性變形。這些給如何實(shí)時(shí)地在槳葉表面生成貼體網(wǎng)格提出了嚴(yán)峻的挑戰(zhàn)。因此,早期的旋翼CFD研究都未曾計(jì)入槳葉的上述運(yùn)動(dòng)。直到近十幾年,Webster和Boniface等人才在旋翼網(wǎng)格生成中考慮槳葉的變距和揮舞運(yùn)動(dòng),它們都采用了動(dòng)態(tài)網(wǎng)格的思想。其中,前者將坐標(biāo)系選在槳尖平面,這樣做的好處是在此平面中只要考慮槳葉的周期性變距。先在槳葉周圍生成初始網(wǎng)格,并依據(jù)到槳葉表面的距離遠(yuǎn)近將槳葉周圍的流場(chǎng)分成三個(gè)區(qū)域,靠近槳葉表面的區(qū)域隨槳葉作一致的運(yùn)動(dòng),外部網(wǎng)格為不動(dòng)網(wǎng)格,中間網(wǎng)格則按一定的權(quán)函數(shù)作變形運(yùn)動(dòng)以銜接內(nèi)外網(wǎng)格;后者的基本思想也類似。然而這種方法的主要缺點(diǎn)就是中間網(wǎng)格必須變形,當(dāng)變形嚴(yán)重時(shí)會(huì)影響計(jì)算效率和計(jì)算精度。此外,這種方法的初始網(wǎng)格只能選取同一種結(jié)構(gòu)網(wǎng)格。為彌補(bǔ)上述方法的不足,本文引入嵌套網(wǎng)格,并吸收動(dòng)態(tài)網(wǎng)格的思想,給出了一種運(yùn)動(dòng)嵌套網(wǎng)格,可較好地計(jì)入槳葉實(shí)際的揮舞和變距運(yùn)動(dòng)。同時(shí),對(duì)嵌套網(wǎng)格中存在的瓶頸問題——貢獻(xiàn)單元的搜索進(jìn)行了改進(jìn),給出了一種新的偽貢獻(xiàn)單元搜索法,克服了在非結(jié)構(gòu)網(wǎng)格中如何進(jìn)行信息傳遞的困難,將只能對(duì)結(jié)構(gòu)網(wǎng)格有效的“InverseMap”法推廣到也可適合于非結(jié)構(gòu)網(wǎng)格。1c型網(wǎng)格的生成在生成旋翼網(wǎng)格時(shí),最初從事旋翼CFD的研究者們只考慮槳葉的旋轉(zhuǎn)運(yùn)動(dòng),沒有考慮槳葉的揮舞和變距運(yùn)動(dòng)。旋翼的三維網(wǎng)格可簡(jiǎn)單地通過下面方法得到。首先將槳葉沿展向分成若干個(gè)剖面,在槳葉的每個(gè)剖面通過求解Poisson方程生成貼體的C型網(wǎng)格,槳尖外的網(wǎng)格通過槳尖處的網(wǎng)格用代數(shù)方法在法向保持相同的網(wǎng)格間距結(jié)構(gòu)而生成,旋翼整體網(wǎng)格為C-H網(wǎng)格。本文的C型網(wǎng)格生成分兩步進(jìn)行:一是用文的幾何生成方法生成槳葉剖面的C型網(wǎng)格;二是將上一步的網(wǎng)格作為初始網(wǎng)格來(lái)求解Poisson方程,最終得到貼體且近似正交的C型網(wǎng)格。文的方法是在物理平面內(nèi)用兩組拋物線相交的辦法構(gòu)筑接近于正交的繞翼型的C型網(wǎng)格,其中一組拋物線過翼型表面點(diǎn),另一組由選取的一組大于零的數(shù)Fj(j=1,…,J,且F1<F2<…<FJ)來(lái)確定。網(wǎng)格點(diǎn)坐標(biāo)(x,y)由下面兩式得到xi,j=Ai2?(Si+Fj)2yi,j=2Ai(Si+Fj)(1)xi,j=Ai2-(Si+Fj)2yi,j=2Ai(Si+Fj)(1)式中,Ai和Si是與翼型表面點(diǎn)有關(guān)的參數(shù),可由這些點(diǎn)的坐標(biāo)解析表出。為了使生成的網(wǎng)格較貼近翼型表面,本文將Fj更改為Fi,j,則F不僅在j方向(近似垂直于翼型表面的方向)上發(fā)生變化,在i方向(近似平行于翼型表面的方向)上也要發(fā)生變化,更有助于調(diào)整網(wǎng)格間距。有了上述生成的C型網(wǎng)格,就可用來(lái)求解Poisson方程。設(shè)計(jì)算平面與物理平面的變換關(guān)系為a=a(x,y),Z=Z(x,y),其中a,Z滿足橢圓型方程。為了便于進(jìn)行離散化數(shù)值求解,將物理平面內(nèi)的橢圓形方程變換到計(jì)算平面內(nèi),即該方程的邊界條件較容易確定。其右端源項(xiàng)有多種取法,為了使求解更簡(jiǎn)單,本文將P,Q取值為式中:a,c為給定的正數(shù),Zmin為邊界面的Z值,這樣可以將網(wǎng)格線向翼型表面靠攏。用中心差分法將(2)式離散,得到三對(duì)角線性方程組。以第一步得到的網(wǎng)格為初始網(wǎng)格,用追趕法求解此線性方程組,可求得繞翼面的C型網(wǎng)格。為了加速收斂,對(duì)翼型網(wǎng)格的外邊界采用內(nèi)點(diǎn)外插的方法來(lái)進(jìn)行更新。實(shí)踐證明,這一套方法很有效。圖1(a,b)為圍繞NACA0012翼型,分別采用第一步方法和第二步方法得到的網(wǎng)格,但此時(shí)生成的網(wǎng)格與翼面貼近程度還不夠;于是采用文中的指數(shù)函數(shù)加密法得到任意貼體的C型網(wǎng)格,如圖2(a)所示。從圖中可以看出,在翼面附近,網(wǎng)格的貼體性很好,但其正交性沒有得到很好的保證。針對(duì)這種情況,本文將物理平面中內(nèi)邊界點(diǎn)取在翼型表面,采用上述第一步、第二步方法生成繞翼型的C型網(wǎng)格,再進(jìn)行加密處理,最終得到在翼面附近貼體且近似正交的C型網(wǎng)格,見圖2(b)。2結(jié)構(gòu)網(wǎng)格的網(wǎng)格點(diǎn)若考慮槳葉的周期揮舞和變距運(yùn)動(dòng),須先給出槳葉在每一個(gè)方位角和每一個(gè)徑向位置的揮舞、變距角的大小,由下式給出式中:j(t)=Kt,表示槳葉所處的方位角;t為時(shí)間;U0,U1c和U1s分別為旋翼?yè)]舞運(yùn)動(dòng)系數(shù);θ0為總距;θ1c和θ1s分別為縱、橫向周期變距;Δθ為槳葉負(fù)扭度。為了能夠確定動(dòng)態(tài)網(wǎng)格中計(jì)及槳葉變距和揮舞運(yùn)動(dòng)的區(qū)域,在槳葉開始旋轉(zhuǎn)前,先生成一套初始網(wǎng)格Ω0,并在此網(wǎng)格中定義三個(gè)子域式中:為槳葉表面及其附近區(qū)域,定義t>0時(shí),網(wǎng)格隨槳葉作一致的揮舞和變距運(yùn)動(dòng),不發(fā)生變形;為過渡網(wǎng)格區(qū)域,也稱中間網(wǎng)格,定義t>0時(shí),網(wǎng)格開始出現(xiàn)變形;為槳葉的整體旋轉(zhuǎn)運(yùn)動(dòng)網(wǎng)格,定義t>0時(shí),網(wǎng)格隨槳葉同步旋轉(zhuǎn),若在槳葉固定的參考系內(nèi),該區(qū)域中的網(wǎng)格不作運(yùn)動(dòng)。上述方法定義的三個(gè)子域在動(dòng)態(tài)網(wǎng)格中所起的作用不同,域可保持網(wǎng)格在翼型表面的貼體性,該區(qū)域中網(wǎng)格點(diǎn)坐標(biāo)變化的表達(dá)式為式中Δθ和ΔU分別表示槳葉變距和揮舞角的改變量。實(shí)現(xiàn)子域中的網(wǎng)格變形是使用動(dòng)態(tài)網(wǎng)格的關(guān)鍵,文給出了一種依據(jù)到外邊界面的距離而定義權(quán)重函數(shù)的方法來(lái)實(shí)現(xiàn)子域中網(wǎng)格的變形,其網(wǎng)格點(diǎn)坐標(biāo)更新公式是x′new(i,j,k)=x(i,j,k)+w(i,j,k)?(xnew(i,j,k)?x(i,j,k))y′new(i,j,k)=y(i,j,k)+w(i,j,k)?(ynew(i,j,k)?y(i,j,k))z′new(i,j,k)=z(i,j,k)+w(i,j,k)?(znew(i,j,k)?z(i,j,k))(7)x′new(i,j,k)=x(i,j,k)+w(i,j,k)?(xnew(i,j,k)-x(i,j,k))y′new(i,j,k)=y(i,j,k)+w(i,j,k)?(ynew(i,j,k)-y(i,j,k))z′new(i,j,k)=z(i,j,k)+w(i,j,k)?(znew(i,j,k)-z(i,j,k))(7)式中:xnew,ynew和znew由式(6)確定,w為權(quán)重函數(shù)。然而,當(dāng)槳葉揮舞和變距角比較大時(shí),子域中的網(wǎng)格會(huì)發(fā)生很大的扭曲變形。3偽貢獻(xiàn)單元法和psde動(dòng)態(tài)網(wǎng)格中會(huì)出現(xiàn)網(wǎng)格扭曲,如果扭曲過大,將會(huì)導(dǎo)致嚴(yán)重的后果;同時(shí)動(dòng)態(tài)網(wǎng)格的信息傳遞只能在各子域的交界面處進(jìn)行,這樣傳遞的信息量很少,對(duì)整個(gè)流場(chǎng)的計(jì)算精度和速度有一定的影響。如果采用嵌套網(wǎng)格,就可以解決上述兩個(gè)問題。具體方法為,在計(jì)算域中劃分子域,針對(duì)每個(gè)子域生成網(wǎng)格,并允許子域網(wǎng)格相互有重疊。在應(yīng)用嵌套網(wǎng)格時(shí)主要存在兩方面問題:(1)洞邊界和洞點(diǎn)的確定;(2)貢獻(xiàn)單元的搜尋。洞邊界和洞點(diǎn)的確定相比于貢獻(xiàn)單元的搜尋要簡(jiǎn)單許多,可直接采用常用的HoleMap法。目前,對(duì)于貢獻(xiàn)單元的搜尋,基本都采用Meakin提出的InverseMap法。然而,此法只對(duì)結(jié)構(gòu)網(wǎng)格有效,對(duì)非結(jié)構(gòu)網(wǎng)格則不可行。同時(shí)即使對(duì)于結(jié)構(gòu)網(wǎng)格,對(duì)于生成InverseMap的疏密程度有一定要求,且還要單獨(dú)處理可能出現(xiàn)的奇異面(線)問題。為了增加InverseMap方法的通用性,本文找到了一種不依賴于網(wǎng)格節(jié)點(diǎn)序號(hào)的方法,被作者稱為“偽貢獻(xiàn)單元的搜尋法(Pseudo-searchingschemeofdonorelements,PSDE)”。圖3(a)為圍繞槳葉剖面的二維非結(jié)構(gòu)網(wǎng)格,P為需要在F中確定其貢獻(xiàn)單元的點(diǎn),FM為覆蓋于F上的InverseMap(見圖3(b)),是均勻的正交直線網(wǎng)格。下面以二維三角形非結(jié)構(gòu)網(wǎng)格為例說明該方法的基本步驟:(1)角形單元面積具體做法為,首先尋找包含F(xiàn)M網(wǎng)格點(diǎn)PM的F網(wǎng)格單元,記為I;然后再確定PM在該網(wǎng)格單元中的相對(duì)位置,計(jì)為ai=Ai/A(i=1,2)。其中A1和A2分別表示PM點(diǎn)與I網(wǎng)格單元中2,3兩點(diǎn)和3,1兩點(diǎn)所構(gòu)成三角形的面積(如圖4所示),A為第I個(gè)三角形單元面積。此處用面積坐標(biāo)作為插值函數(shù),該公式為式中:ui(i=1,2,3)為三角形單元3個(gè)頂點(diǎn)上的函數(shù)值;u為PM點(diǎn)函數(shù)值。(2)偽貢獻(xiàn)單元abcd設(shè)P點(diǎn)在槳葉剖面當(dāng)?shù)刈鴺?biāo)系中的坐標(biāo)為(xP,yP)。包含P的FM網(wǎng)格單元可通過兩個(gè)一維搜索容易得到,計(jì)該單元為ABCD(見圖5)。此單元即為需要搜索的貢獻(xiàn)單元,因?yàn)椴皇侵苯釉贔中找到的包含P的單元,故稱為偽貢獻(xiàn)單元。P在該單元中的相對(duì)位置為其中點(diǎn)(xC,yC)為單元ABCD的中心點(diǎn)。(3)abcd單元先求出偽貢獻(xiàn)單元4個(gè)頂點(diǎn)的函數(shù)值,可直接由第一步計(jì)算得到,接著根據(jù)P點(diǎn)在此貢獻(xiàn)單元中的相對(duì)位置,通過雙線性插值獲得P點(diǎn)的函數(shù)值,即式中,ui為ABCD單元頂點(diǎn)上的函數(shù)值。通過上述三步方法便可以完成嵌套網(wǎng)格的信息傳遞。與文方法相比,本文方法更通用,對(duì)InverseMap沒有嚴(yán)格的疏密要求,視精度而定。同時(shí),不必單獨(dú)處理奇異線(面)問題,因?yàn)閭呜暙I(xiàn)單元搜尋法對(duì)網(wǎng)格節(jié)點(diǎn)的序號(hào)沒有要求。4網(wǎng)格內(nèi)邊界點(diǎn)的變化下面基于嵌套網(wǎng)格,并結(jié)合動(dòng)態(tài)網(wǎng)格的思想,給出了一種運(yùn)動(dòng)嵌套網(wǎng)格。具體方法是,先將旋翼周圍的物理空間分成幾個(gè)區(qū)域,區(qū)域間允許有重疊;接著讓槳葉表面的貼體網(wǎng)格隨槳葉一起作周期揮舞和變距運(yùn)動(dòng),網(wǎng)格點(diǎn)的變化按式(6)進(jìn)行。在每一個(gè)時(shí)間站位重新確定邊界點(diǎn)(如果相互重疊的網(wǎng)格都發(fā)生了運(yùn)動(dòng)),即確定相互重疊的網(wǎng)格間的人工內(nèi)邊界點(diǎn)和人工外邊界點(diǎn);然后,再通過偽貢獻(xiàn)單元搜尋法來(lái)實(shí)現(xiàn)不同區(qū)域網(wǎng)格間信息的傳遞,從而可實(shí)現(xiàn)整個(gè)流場(chǎng)的計(jì)算。圖6給出了在槳葉變距和揮舞平面內(nèi)的運(yùn)動(dòng)嵌套網(wǎng)格圖,圖7則給出了運(yùn)動(dòng)嵌套網(wǎng)格的總體示意圖。5偽貢獻(xiàn)單元搜尋法(1)結(jié)合嵌套網(wǎng)格和運(yùn)動(dòng)網(wǎng)格的特點(diǎn),給出了一種