不可壓縮圓柱繞流坐標(biāo)變換SOR方法

1、不可壓縮圓柱繞流摘 要 本文對無粘不可壓圓柱繞流問題進行數(shù)值模擬，并將數(shù)值結(jié)果和理論解進行對比。在給定圓柱固壁和遠場速度的邊界條件下，通過求解流函數(shù)的Laplace方程，得到域內(nèi)的流函數(shù)，進而計算出速度場。數(shù)值模擬的主要步驟有網(wǎng)格劃分，物理域和計算域坐標(biāo)變換，對流函數(shù)方程進行差分離散， SOR迭代方法求解和后處理繪出速度場。本文將先建立本問題的數(shù)學(xué)提法和理論解，隨后基于數(shù)學(xué)提法進行數(shù)值求解，最后將理論和數(shù)值進行對比。理論和數(shù)值得到的速度場表現(xiàn)出了一致性，并且均反映出無粘假設(shè)下的流動和實際流動間的巨大差異。關(guān)鍵詞 不可壓縮；圓柱繞流；坐標(biāo)變換；SOR方法1 問題提出考慮圖 1所示無粘不可壓圓柱繞

2、流，來流速度為u=20m/s,v=0。選取合適的計算域范圍和合理的邊界條件，通過求解流函數(shù)方程或者速度勢方程求圓柱周圍的流場并與理論解相比較。圖 1 無粘不可壓圓柱繞流示意圖2 理論求解在圓柱圓心處取笛卡爾坐標(biāo)，對于不可壓縮二維流動，存在流函數(shù)有 無窮遠處來流u=20m/s，為無旋流動，根據(jù)Kelvin-Lagrange定理可知，本問題中，流場始終無旋。由，代入可得主控方程定解的邊界條件有兩個，一個是圓柱的固壁條件，考慮到固壁的不可穿透性有 上式中為圓柱半徑，第二個邊界條件是無窮遠處速度分布在圓柱坐標(biāo)系中求解上述問題，坐標(biāo)定義如圖 2，可得到坐標(biāo)變換關(guān)系和微分關(guān)系圖 2 坐標(biāo)關(guān)系說明根據(jù)進行坐

3、標(biāo)轉(zhuǎn)換，可得到主控方程和邊界條件如下 通過分離變量法1，可得流函數(shù)的解為根據(jù)，算得理論速度場為 3 數(shù)值求解3.1 網(wǎng)格劃分與坐標(biāo)變換設(shè)計求解物理域為與圓柱同軸的圓環(huán)域，范圍為 rc=20m，r=100m。為保證網(wǎng)格正交性，采用O型網(wǎng)格對物理域進行劃分，如圖 3所示。網(wǎng)格沿徑向有161個節(jié)點，間距為dr=0.5m；沿周向有361個節(jié)點，間距為d=1°，第一個節(jié)點和第361節(jié)點在物理域上是重合的。圖 3 物理域網(wǎng)格取計算域為Cartesian網(wǎng)格，坐標(biāo)系記為-坐標(biāo)，每個網(wǎng)格的大小為1×1，即=1。計算域中點(i,j)對應(yīng)物理域中點(rcos,rsin)，其中r=rc+dr(i

4、-1), =(j-1)d。根據(jù)此映射關(guān)系，可以數(shù)值求得從x-y到-的變換微分關(guān)系、J、g11、g12和g22對于二維平面任意曲線坐標(biāo)系中的測度以及變化Jacobian行列式J的表達式為度量張量為3.2 主控方程離散與迭代求解由于坐標(biāo)變換從x-y到-，因此流函數(shù)方程2=0變換為 引入中間變量 則有 按圖 4所示節(jié)點排列進行差分：圖 4 節(jié)點離散示意圖在點進行差分離散： 其中，中間變量差分為： 整理以上5個式子，即可算出i,j*，利用SOR迭代方法得到新的i,jn+13.3 邊界條件根據(jù)邊界條件和邊界點的坐標(biāo)變換，代入u=20，v=0，可推出邊界處的流函數(shù)梯度和；在圓柱表面，取r=rc=0。據(jù)此即

5、給定了邊界條件。計算域中所有點的流函數(shù)初值均為0，圓柱表面流函數(shù)值始終為0。除了內(nèi)外邊界，其他點均可由主控方程算出數(shù)值。計算時，先算=2時各點的值，依次類推，直到算完=160的各點的值。計算外邊界，即=161各點時，采用如下辦法這樣使得計算外邊界流函數(shù)時，能同時利用到兩個方向的梯度和。3.4 收斂條件計算兩輪迭代的流函數(shù)點陣差值的二范數(shù)，當(dāng)該值小于10-5時，即判定為收斂。4 計算結(jié)果及對比4.1 數(shù)值解圖 5為數(shù)值解的流函數(shù)分布，流函數(shù)等值線即流線，因此，根據(jù)該圖可以很明顯看到不可壓無粘繞流特征：即流體在圓柱前一分為二，在圓柱后重新合二為一，流函數(shù)關(guān)于x=0對稱分布，關(guān)于y=0反對稱分布。這

6、個分布很好的反映了該流動無損失，且圓柱受到合力為0的理論估計。圖 5 數(shù)值仿真所得流函數(shù)分布圖 6、圖 7和圖 8分別為u分量、v分量和速率分布，再次驗證了分流特征和關(guān)于流場對稱性的分析。值得注意的是，圓柱前后均有速度為0的滯止區(qū)域，近圓柱面的繞流先加速，在圓柱上下頂點達到最大值，之后減速離開圓柱。圖 6 數(shù)值解速度u分量分布圖 7 數(shù)值解速度v分量分布圖 8 數(shù)值解速率分布根據(jù)Cauchy-Lagrange積分這里是勢函數(shù)，對二維定常問題，且不計體積力，有假定p=100kPa，=1.2kg/m3可以求得全場的壓力分布，如圖 9所示，首先驗證了圓柱所受合力為零，另外，在前后滯止區(qū)內(nèi)有壓力最大值，在上下頂點有壓力最小值。圖 9 數(shù)值解壓力分布4.2 數(shù)值與理論對比記理論所得流函數(shù)為T,則流函數(shù)偏差分布e=-T,如圖所示。圖 10流函