基于離散渦方法的粒子強(qiáng)度交換法模擬流場

基于離散渦方法的粒子強(qiáng)度交換法模擬流場

0基于隨機(jī)渦流的彈性渦流模擬方法,主要特點(diǎn)如下該市的方法是模擬非定常流的拉格朗日方法，可以忽略離散誤差。該方法具有良好的強(qiáng)度特性。這兩個(gè)特點(diǎn)使得渦方法成為模擬流場的有效的工具。二維渦方法發(fā)展最成熟的是由Chorin提出的隨機(jī)走步渦團(tuán)法,它是用隨機(jī)走步來模擬渦量的粘性擴(kuò)散,具有統(tǒng)計(jì)平均意義,該方法已經(jīng)成功地運(yùn)用到了橋梁斷面的流場模擬和抗風(fēng)分析中,如AllanLarsen和Zhou。但是它的缺點(diǎn)不能忽視,對(duì)長時(shí)間歷程問題,誤差急劇增長,隨機(jī)走步模擬擴(kuò)散精度較低,且用新生渦團(tuán)處理無滑移邊界條件有較高的數(shù)值噪聲,確定性渦方法克服了隨機(jī)渦方法的缺點(diǎn),而不改變渦方法作為Lagrange方法的優(yōu)點(diǎn)。過去的二十年里,在渦方法領(lǐng)域內(nèi),已經(jīng)取得了較大的發(fā)展:快速多極子展開的應(yīng)用(Greedgard和Rohklin,Barnes和Hut)用粒子強(qiáng)度交換法精確地模擬粘性擴(kuò)散(Degond和Mas-Gallic),以及精確施加粘性邊界條件(Koumoutsakos,etal,Mas-Gallic,Leonardetal,Ploumhans,etal,Benhaddouch)。所有這些使得渦方法能夠更高精度地模擬具有粘性條件的流場。例如,Koumoutsakos和Leonard,Koumoutsakos和Shiels,用確定性渦方法分別模擬了突然起動(dòng)的圓柱和平板的流場。當(dāng)然,他們這種精確的方法也只應(yīng)用到了特殊的截面(圓形和平板)。目前,粒子強(qiáng)度交換法尚未應(yīng)用到工程中,且國內(nèi)對(duì)這種方法的研究比較少,對(duì)其在工程上的應(yīng)用更是空白。本文對(duì)粒子強(qiáng)度交換法進(jìn)行了研究并程序化,將其初步應(yīng)用到橋梁截面的靜風(fēng)分析中。1確定波方法1.1邊界條件條件渦量場的二維Navier-Stokes方程可以匯總為如下數(shù)學(xué)模型?ω?t+u??ω=ν?2ωinD(1a)?2Ψ=-ωinD(1b)u=?×ΨinD(1c)ω(x,0)=ω0(x)inD(1d)u=Ub(t)+Ω(t)ez×(xs-xb)on?D(1e)u→U∞at∞(1f)?ω?t+u??ω=ν?2ωinD(1a)?2Ψ=?ωinD(1b)u=?×ΨinD(1c)ω(x,0)=ω0(x)inD(1d)u=Ub(t)+Ω(t)ez×(xs?xb)on?D(1e)u→U∞at∞(1f)其中ω為渦量,ν、ρ分別為流體的運(yùn)動(dòng)粘性系數(shù)和流體的密度,Ψ為流函數(shù),固體的平動(dòng)速度為Ub(t),繞其形心xb的角速度為Ω(t),U∞是無窮遠(yuǎn)處的速度,式(1b)為泊松方程,(1c)為不可壓縮條件,(1d)為初始條件,(1e)為固體邊界條件,(1f)為遠(yuǎn)場邊界條件。以上方程和二維Navier-Stokes方程是等效的,渦量-流函數(shù)方程使得方程中壓力項(xiàng)消失,但是對(duì)于有邊界的流體域,它引入了附加的運(yùn)動(dòng)約束,需要將速度邊界條件轉(zhuǎn)變?yōu)闇u量邊界條件。1.2didt、dh、ni本文中的粘性擴(kuò)散采用粒子強(qiáng)度交換法,該方法是用積分算子來估計(jì)拉普拉斯算子,?2ω(x)≈2σ2∫ησ(x-y)(ω(y)-ω(x))dy(2)?2ω(x)≈2σ2∫ησ(x?y)(ω(y)?ω(x))dy(2)其中ησ(x)=(1/σ2)η(|x|/σ)?η(s)=-1sddsζ(s)?ζ(s)=12πexp(-s22),因此dΓidt=2υσ2Ν∑j=1ησ(xi-xj)(ΓjSi-ΓiSj)(3)Γni=Γn-1i+2υσ2Ν∑j=1ησ(xi-xj)(Γn-1jSi-Γn-1Sj)(4)在實(shí)際運(yùn)用中,由于η的迅速降低,只有靠近粒子i的粒子才對(duì)dΓi/dt有明顯的影響,為了減少計(jì)算量,只考慮靠近粒子i的粒子。稱對(duì)dΓi/dt有明顯影響的粒子集為Di,如果采用高斯光滑核,Di包括離粒子xi距離小于5σ所有粒子。Γni?Γn-1i+2υσ2∑j∈Diησ(xi-xj)(Γn-1jSi-Γn-1iSj)(5)當(dāng)所有粒子的不同時(shí),方程(5)變?yōu)?Γni?Γn-1i+2υσ2ij∑j∈Diηij(xi-xj)(Γn-1jSi-Γn-1iSj)(6)其中ηij(xi-xj)=(1/σ2ij)η(|xi-xj|/σij),σ2ij=(σi+σj)/2。1.3幾何標(biāo)準(zhǔn)指數(shù)網(wǎng)格的構(gòu)建粒子分布的開始階段,像線動(dòng)量和角動(dòng)量等全局量一般不會(huì)有誤差產(chǎn)生,隨著時(shí)間的增長流場中的粒子會(huì)失去重疊或重疊太多,由于局部粒子的扭曲不能正確顯示渦結(jié)構(gòu),而粒子的聚集導(dǎo)致一些非物理現(xiàn)象產(chǎn)生,從而不能精確地求解光滑的渦域。因此在粒子強(qiáng)度交換法中,粒子強(qiáng)度重分布是必不可少的一步。采用四階內(nèi)插核,滿足渦量守恒,渦量的一階矩,二階矩和三階矩守恒,Λ3為粒子周圍沒有邊界的情況,Λ′2為邊界位于x=-12?Λ′3邊界位于x=-1,網(wǎng)格間距為1,單位強(qiáng)度的粒子位于-12≤x≤12,新粒子的強(qiáng)度分別為Γ1、Γ2、Γ3、Γ4,Λ3中新粒子的位置為-32、-12、12、32?Λ′2中新粒子的位置分別為0、1、2,Λ′3中新粒子的位置分別為-12、12、32、52,其方程分別為式(7-9)。二維示意圖如圖1所示。在外部流場中,渦量主要位于邊界層和尾流中,渦量的大小隨著渦流向下游逐漸變小。因此需要在固體邊界采用更高的精度而在尾流中可以降低精度。為了實(shí)現(xiàn)這個(gè)目標(biāo),將具有規(guī)則序數(shù)不均勻的網(wǎng)格(i,j)覆蓋精度逐漸變化的物理域上。當(dāng)然可以考慮不同的網(wǎng)格形式,本文采用文獻(xiàn)中提出的指數(shù)網(wǎng)格形式。粒子在流體中的面積為:S=h2?((2πmR0)?exp(2πm(i+12)))2(10)h?(2πmR0)?exp(2πm(i+12))(11)其中m為軸向劃分的個(gè)數(shù),R0為最小圓的半徑,h為網(wǎng)格寬度,S為粒子在流體中的面積。2表面環(huán)量與子強(qiáng)度交換的模擬結(jié)果對(duì)雷諾數(shù)為550的突然起動(dòng)圓柱的流場進(jìn)行模擬,采用指數(shù)網(wǎng)格形式,m取592,其粒子的分布如圖2所示(圖中間距是計(jì)算所取間距的四倍即m=148),時(shí)間步取0.01,表面渦片長度為0.0053,對(duì)流和粒子強(qiáng)度交換模擬擴(kuò)散過程均采用歐拉格式,表面渦片向流場內(nèi)擴(kuò)散過程采用高斯四點(diǎn)積分,將計(jì)算的圓柱上表面產(chǎn)生的環(huán)量dΓup/dt=υR∫π0?ω?tdθ和其他文獻(xiàn)的結(jié)果做了比較,如圖3所示,其中KL為Kou-moutsakos和Leonard方法的計(jì)算結(jié)果,G、H、V分別為Ploumhans提出的一般性方法,圓柱表面加“暈環(huán)”和粒子強(qiáng)度重分配時(shí)使用“隨機(jī)振動(dòng)”,詳細(xì)參考文獻(xiàn),從圖中可以看出本文也取得了較好的結(jié)果。同時(shí)對(duì)阻力系數(shù)和Ix=∫ΩyωdΩ分別與文獻(xiàn)和解析解做了比較,阻力系數(shù)如圖4所示,Ix的比較見圖5。圖6給出了突然起動(dòng)圓柱的渦量等勢線,并于Ploumhans的結(jié)果進(jìn)行了比較,從圖中可以看出計(jì)算結(jié)果比較理想。3阻力系數(shù)的時(shí)程曲線為了進(jìn)一步驗(yàn)證程序的正確性,計(jì)算了風(fēng)攻角為15°的方形流場,計(jì)算方法與圓形的方法相同,只是指數(shù)網(wǎng)格的圓心不在圓點(diǎn),與文獻(xiàn)中一致,取在(-2.5,-0.67),m取為850,計(jì)算的阻力系數(shù)如圖7所示。圖7也給出了文獻(xiàn)中阻力系數(shù)的時(shí)程曲線,盡管阻力系數(shù)的振動(dòng)幅值有所不同,但是其平均值和周期還是比較接近的。圖8給出了本文與文獻(xiàn)計(jì)算的在T=30的渦量等勢線的比較。從阻力系數(shù)和渦量等勢線的比較均可看出本文均取得了較為理想的計(jì)算結(jié)果。4計(jì)算誤差分析大海帶東橋引橋主梁截面為閉口箱形截面,其尺寸如圖9所示,將粒子強(qiáng)度交換法首次應(yīng)用到橋梁抗風(fēng)計(jì)算中,本文選擇大海帶橋進(jìn)行數(shù)值模擬,為了使得網(wǎng)格更為貼體,對(duì)本橋梁斷面的計(jì)算選用了與笛卡爾坐標(biāo)平行的方形網(wǎng)格,網(wǎng)格的間距取為0.01,無量綱時(shí)間步為0.01,雷諾數(shù)為Re=UB/υ=104,其中U為風(fēng)速,B為橋?qū)?υ為運(yùn)動(dòng)粘性系數(shù)。計(jì)算的三分力系數(shù)時(shí)程如圖10所示,并與隨機(jī)渦方法計(jì)算結(jié)果的比較,從圖中可以看出,兩個(gè)計(jì)算的結(jié)果吻合得很好,但是隨機(jī)渦方法計(jì)算的結(jié)果振蕩比較明顯,確定性渦方法計(jì)算的曲線比較光滑。取阻力系數(shù)時(shí)程穩(wěn)定階段進(jìn)行平均,得到平均阻力系數(shù),且對(duì)升力系數(shù)時(shí)程進(jìn)行傅立葉變換得到St=fH/U,其中f為旋渦頻率,H為截面高度。表1給出了確定性渦方法計(jì)算的阻力系數(shù)和St數(shù)與其他數(shù)值方法和實(shí)驗(yàn)計(jì)算結(jié)果的比較,計(jì)算結(jié)果表明確定性渦方法的結(jié)果和實(shí)驗(yàn)結(jié)果比較接近。圖11給出了確定性渦方法計(jì)算的渦量等勢線。圖12為確定性渦方法計(jì)算的局部速度和隨機(jī)渦方法的局部