雙向交錯(cuò)迭代耦合模擬計(jì)算方法研究

雙向交錯(cuò)迭代耦合模擬計(jì)算方法研究

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