雙向交錯迭代耦合模擬計算方法研究

雙向交錯迭代耦合模擬計算方法研究

水資源是現(xiàn)有能源中唯一的可大規(guī)模開發(fā)能源?，F(xiàn)在已建成的一系列大型節(jié)水設(shè)施吸引了很多關(guān)注和開發(fā)的資源，但大型混流式機械的負荷變化頻繁，運行條件差。國內(nèi)調(diào)試的主要原因是振動大，葉片破裂。機組振動影響水輪機的正常穩(wěn)定運行,嚴重的甚至?xí)斐烧麄€廠房及大壩發(fā)生共振,危及電站安全運行,至今仍沒有很好的解決方案。轉(zhuǎn)輪葉片是水輪機的重要水力部件,轉(zhuǎn)輪葉片受流體自激誘發(fā)振動,會造成整個機組不穩(wěn)定和廠房的振動。本文主要研究轉(zhuǎn)輪葉片在流動激勵下的誘發(fā)振動,研究流體動壓力脈動對轉(zhuǎn)輪葉片振動的影響,以解釋水輪機所存在的一直無法很好解決的水力振動問題。由于轉(zhuǎn)輪葉片在復(fù)雜的黏性流場作用下的振動問題是一個非常復(fù)雜的多場耦合的非線性動力學(xué)問題,它不僅結(jié)構(gòu)復(fù)雜,振動原因很多,而且振動機理尤為復(fù)雜。水輪機的振動是由水力、電氣、機械等多種因素引起的,在水輪機運行過程中,流體激勵振動、機械振動和電磁振動共存,耦合在一起構(gòu)成一個復(fù)雜的大型耦合振動系統(tǒng)。而轉(zhuǎn)輪葉片的耦合誘發(fā)水力振動也是重要因素之一。當激勵渦頻和薄壁的固有頻率接近時發(fā)生共振。本文主要模擬在動水壓脈動作用下轉(zhuǎn)輪葉片耦合界面上的壓力和位移脈動,并對整個湍流流道進行計算,進而為目前水輪機頻頻出現(xiàn)的振動問題提供依據(jù)。1有限元單元位移、加速度和剛度矩陣的計算在流固耦合的固體部分應(yīng)用瞬時彈性模擬,固體的運動動力學(xué)方程為Μ¨q(t)+C˙q(t)+Κq(t)=f(t)(1)w(x?y?z?t)=Ν∑i=1qi(t)φi(x?y?z)(2)式中:f(t)為外部荷載力;q(t),˙q(t),¨q(t)分別為有限元單元的位移、速度和加速度矢量;M,C,K分別為質(zhì)量、阻尼、剛度矩陣;w(x,y,z,t)為固體在任一瞬時的位置;qi(t)為固體在任一瞬時的位移矢量;φi(x,y,z)為固體的模態(tài);N為全部模態(tài)階數(shù)。2流體密度積分不可壓縮牛頓流體運動的控制方程(Navier-Stokes)在笛卡爾坐標下可表示為?ui?t+ui?ui?xj=-1ρ?p?xi+ν?2ui?xj?xj+fi(3)?ui?xj=0(4)式中:ρ為流體密度;ν為流體運動黏性系數(shù);fi為質(zhì)量力強度。有限體積法對流體在控制體上的積分時體積的取法就是目前所面臨的問題。目前通用的是由Thomas和Lombard提出的GCL(geometricconservationlaw)方法。3接口連接處理3.1固體單元es圖1為界面插值荷載算法,圖中χj是fj上的自然坐標,而在流體節(jié)點上的壓力形函數(shù)為Φj=∫ΓF(-p?n+σF?n)Djds(5)式中:Dj為界面節(jié)點上的函數(shù);Φj為界面上節(jié)點的流體壓力;p為正應(yīng)力;n為法矢量;σF為流體部分的應(yīng)力。在有2個節(jié)點的1個固體邊界單元上所受的來自4個流體單元的4個節(jié)點的壓力表達為f1=Ν1Φ1+Ν1Φ2+Ν1Φ3+Ν1Φ4+Ν1Φ5=Φ1+Ν1Φ2+Ν1Φ3(6)f2=Ν2Φ1+Ν2Φ2+Ν2Φ3+Ν2Φ4+Ν2Φ5=Ν2Φ2+Ν2Φ3+Φ4(7)那么,作用于固體單元的荷載和為f1+f2=Φ1+(Ν1+Ν2)Φ2+(Ν1+Ν2)Φ3+Φ4(8)再來看已知結(jié)構(gòu)位移,插值獲得界面流體的位移值,為說明清楚,把耦合界面Γfs(fs表示流固耦合)分為流體側(cè)Γf(f表示流體)和結(jié)構(gòu)側(cè)Γs(s表示固體),并人為地把Γf和Γs拉開一定的距離。在Γf上有一流體網(wǎng)格節(jié)點fj,Ωes是與fj對應(yīng)的Γs上的結(jié)構(gòu)網(wǎng)格。由固體和流體動網(wǎng)格在耦合界面的位移協(xié)調(diào)條件,可得流體節(jié)點對應(yīng)的網(wǎng)格位移xfcj為xfcj=xfc(fj)=usc(χj)=∑Νi(χj)usci(i∈Γs?j∈Γf)(9)式中:i為固體單元Ωes的節(jié)點總數(shù);Ni為固體單元Ωes上節(jié)點i的插值形函數(shù);usci為固體單元Ωes上節(jié)點i的位移值。3.2流體固體動網(wǎng)格的計算σf?n=σs?n=-p?n+σμ?n(10)其中σμ=2μ(s-13Ι??v)在無滑移界面Γ上?uf=?us˙us=v在滑移界面Γ上uf?n=us?nn?˙uf=n?us式中:σs為固體部分的應(yīng)力;σμ為黏性剪應(yīng)力;μ為固體節(jié)點上的位移;?uf為流體邊界上的位移;?us為固體部分的位移;˙uf為流體邊界上的速度;˙us為固體邊界上的速度。式(10)表示界面Γ滿足的平衡條件,流體在固體上的等效節(jié)點力可由式(10)計算。流體固體動網(wǎng)格在Γ上滿足usc=xfc(11)˙usc=˙xfc=ualec(12)式中:usc為固體邊界上的網(wǎng)格位移;xfc為邊界上流體的網(wǎng)格位移;˙usc為邊界上固體的速度;˙xfc為邊界上流體的速度;ualec為ALE描述下的網(wǎng)格速度。式(11)和式(12)表示邊界上固體和流體網(wǎng)格位移和速度相協(xié)調(diào)。4水輪機動力學(xué)模型計算模型為水輪機全流道流動,模擬對象為某水電站水輪機從蝸殼進口到尾水管出水口全流道內(nèi)的流動,并在轉(zhuǎn)輪葉片處設(shè)流體-固體耦合界面(FSI),固體部分楊氏模量Y=250GPa,泊松比為0.3,密度為7850kg/m3,采用非線性彈性瞬態(tài)模擬。圖2為水輪機轉(zhuǎn)輪模型及網(wǎng)格。針對計算空間復(fù)雜的特點,再考慮耦合界面的動網(wǎng)格和網(wǎng)格變形等因素,流體區(qū)域采用適應(yīng)性較強的非結(jié)構(gòu)化四面體網(wǎng)格劃分技術(shù),將流體區(qū)域劃分為1216146個網(wǎng)格單元(圖2(b))。轉(zhuǎn)輪及其網(wǎng)格見圖2(c)(d)。圖3為葉片上測試點位置。測試點1、測試點2和測試點3分別為葉片壓力面上不同位置的監(jiān)控點,而測試點4、測試點5和測試點6為負壓面上對應(yīng)的監(jiān)控點。在小流量工況下,進口速度取4.5m/s,出口壓力設(shè)為1個大氣壓。5流動壓力脈動所誘發(fā)的水力振動本文將固體采用非線性彈性瞬態(tài)模擬,流體為不可壓縮線性牛頓流體,進行雙向交錯迭代耦合計算,分析轉(zhuǎn)輪葉片流動壓力脈動所誘發(fā)的水力振動。在葉片的壓力面上,對比圖4可知,最大網(wǎng)格位移出現(xiàn)在葉片測試點2的位置。從整個曲線變化趨勢可以看出,在t=0.05s時,葉片上網(wǎng)格位移最大,也就是說葉片擺幅較大,這樣就會導(dǎo)致整個轉(zhuǎn)輪偏心,進而引起整個水輪機振動。