李浩, 宋艷爭, 劉剛*, 呂金潮, 靳立鵬, 武衛(wèi)革

(1.華北電力大學(保定)電氣與電子工程學院, 保定 071003; 2.國網(wǎng)冀北電力有限公司技能培訓中心, 保定 071051;3.保定天威保變電氣股份有限公司河北省輸變電裝備電磁與結(jié)構(gòu)性能重點實驗室, 保定 071056)

油浸式變壓器是電力系統(tǒng)中最重要的設(shè)備之一,其在電力系統(tǒng)中扮演著重要角色。油浸式變壓器的絕緣老化主要受溫升影響,油浸式變壓器溫升與油流速度是一種強耦合關(guān)系[1],油流速度的大小直接影響著溫升。油浸式變壓器在運行中油流不可避免的存在著氣泡,而油流速度又會影響氣泡移動從而影響變壓器的擊穿特性[2]。當油浸式變壓器使用瓦斯保護時,往往由于油流涌動存在著誤動的情況,對于油浸式變壓器出口油流速度與壓力的監(jiān)控顯得尤為重要[3]。由此可看出,準確計算出油浸式變壓器油流速度對其運行、繼電保護及使用壽命具有重要的研究意義。

目前,流場數(shù)值計算方法主要有有限體積法和有限元法。鄧永清等[4]基于自適應(yīng)網(wǎng)格的有限體積法對油流速度進行了計算,減小了計算誤差并提高了計算效率。魏本剛等[5]基于有限體積法對三維變壓器冷卻問題進行了研究。Mufuta等[6]采用有限體積法求解了大型變壓器盤型繞組速度分布。然而采用有限體積法求解復雜邊界結(jié)構(gòu)時,對流通量與擴散通量會出現(xiàn)擴散現(xiàn)象影響精度[7]。因此也有大量學者采用了有限元法對流場問題求解,而有限元法又由于對控制方程處理離散方法不同,分為伽遼金有限元法、迎風有限元法及最小二乘有限元法。而對于大型求解對流占優(yōu)問題時,伽遼金有限元法存在數(shù)值振蕩,所以往往選擇最小二乘有限元法與迎風有限元法。劉剛等[8]基于三角形網(wǎng)格采用迎風有限法計算了變壓器單分區(qū)模型。Liu等[9]運用最小二乘有限元法對變壓器流場流動進行了計算。最小二乘有限元法(least-squares finite element method,LSFEM)是在有限元法的基礎(chǔ)上,通過最小二乘的思想,將離散方程組的系數(shù)矩陣轉(zhuǎn)化為對稱正定的稀疏矩陣[10],收斂性較好,且不會出現(xiàn)數(shù)值振蕩問題。

但最小二乘有限元法在計算流場時,由于高斯積分公式及方程組超定問題,使用四邊形網(wǎng)格計算會更加準確。然而對于含有角環(huán)、圓形和尖角等復雜部件的電力裝備而言,如果僅采用四邊形網(wǎng)格剖分,會出現(xiàn)網(wǎng)格形狀畸變或者剖分節(jié)點畸變,甚至剖分節(jié)點無法與邊界重合等諸多問題,此時就無法采用最小二乘有限元法進行計算。而三角形網(wǎng)格具有很高的幾何靈敏性,對于復雜結(jié)構(gòu)具有良好的貼合性[11]。三角形網(wǎng)格允許在各種形狀的單元之間生成高質(zhì)量網(wǎng)格,例如,在線圈之間的狹窄間隙內(nèi),或在本模型范圍弧狀結(jié)構(gòu)內(nèi)。此外,它有利于粗四邊形和細四邊形網(wǎng)格區(qū)之間的過渡。張?zhí)斓萚12]用混合網(wǎng)格模擬了復雜結(jié)構(gòu)水下潛器的操作。Alauzet等[13]采用各向異性網(wǎng)格自適應(yīng)對于葉輪旋轉(zhuǎn)進行了模擬計算。針對復雜結(jié)構(gòu)無法簡易剖分問題,目前國內(nèi)外學者主要的研究方法是使用混合網(wǎng)格或者自適應(yīng)網(wǎng)格。

針對上述復雜結(jié)構(gòu)二維流場計算的需求,現(xiàn)提出一種三角形網(wǎng)格再處理為四邊形的方法,使其適用于最小二乘有限元法。為驗證方法有效性,針對方腔頂蓋驅(qū)動模型,用所提方法的計算結(jié)果與Fluent仿真結(jié)果進行對比。最后將所提方法應(yīng)用于具有角環(huán)、圓角等復雜結(jié)構(gòu)的換流變壓器單分區(qū)模型的油流分布計算。

1 無量綱最小二乘有限元法計算油流速度的控制方程

變壓器油常常被近似視為不可壓縮流體,滿足質(zhì)量守恒與動量守恒定律,穩(wěn)態(tài)情況下的納維斯托克斯方程可表示為如下形式。

?U=0

(1)

(ρU·?)U+?p=μ?2U+f

(2)

式中:U為速度,m/s;ρ為流體密度,kg/m3;μ為動力黏度,Pa·s;f為外應(yīng)力,N;p為壓強,Pa;?為拉普拉斯算子。其中式(2)中等式左邊第一項為對流項,第二項為壓力梯度項。等式右邊第一項為黏性動量,第二項為外應(yīng)力。

由于式(2)壓力-速度耦合方程是速度的二階偏微分方程,無法直接使用最小二乘有限元法。為此引入渦量ω=?×U,將上述方程轉(zhuǎn)換為一階的壓力-速度-渦量耦合方程可得

(ρU·?)U+?p=μ?×ω+f

(3)

ω-?×U=0

(4)

?·U=0

(5)

對每個變量引入?yún)⒖剂?將各變量與常量無量綱化[14],在直角坐標系中式(3)～式(5)可展開成如下形式。

(6)

將式(6)一階偏微分方程描述成矩陣方程形式為

(7)

將式(7)中各變量表示成矩陣形式為

(8)

(9)

對矩陣形式的式(7)通過殘差公式,以及最小二乘法的基本思想可化為

[A(u+tv]-f)dΩ

(10)

對式(9)進行處理,代入插值函數(shù),便可以得到離散方程組。

其中系數(shù)矩陣A為對稱正定的稀疏矩陣,最小二乘有限元法的詳細推導過程請參考文獻[14]。

2 基于三角形網(wǎng)格的LSFEM實現(xiàn)過程

對于圖1所示的半圓形仿真模型,在采用四邊形網(wǎng)格剖分時,如果采用自由剖分往往會出現(xiàn)網(wǎng)格單元畸變,即出現(xiàn)部分三角形網(wǎng)格的現(xiàn)象,如圖1的局部放大圖所示。

圖1 規(guī)則四邊形網(wǎng)格畸變?yōu)榉墙Y(jié)構(gòu)三角形

然而對于大型電力裝備而言,由于部件多、形狀比圖1所示模型復雜得多,在網(wǎng)格剖分時很難實現(xiàn)僅用四邊形網(wǎng)格全場域的剖分。此時,如果采用三角形剖分,則具有很好的靈活性,且剖分效果較好,如圖2所示。

圖2 非結(jié)構(gòu)三角形剖分弧線結(jié)構(gòu)

因此,為了在復雜結(jié)構(gòu)的二維流體場分析中使用無量綱最小二乘有限元法,將通過對三角形網(wǎng)格再處理得到適合最小二乘有限元法使用的四邊形網(wǎng)格,過程如下。

步驟1將仿真模型用三角形進行剖分;再對三角形網(wǎng)格單元節(jié)點進行再處理得到四邊形網(wǎng)格,處理流程如圖3所示。

圖3 八節(jié)點四邊形網(wǎng)格形成過程

步驟2在三角形剖分得到圖3(a)的網(wǎng)格后,將單元關(guān)聯(lián)矩陣和節(jié)點坐標矩陣分別儲存。初始三角形含有6個節(jié)點,其節(jié)點關(guān)聯(lián)矩陣的節(jié)點序號如表1所示。

表1 三角形6節(jié)點單元節(jié)點關(guān)聯(lián)矩陣

步驟3在圖3(b)中,先形成一個新節(jié)點,即取三角形3個頂點橫、縱坐標的平均值作為新節(jié)點的坐標,從而形成新的節(jié)點關(guān)聯(lián)矩陣和坐標矩陣。

步驟4在圖3(c)中,需要將7節(jié)點三角形單元形成3個新的四節(jié)點四邊形單元,即提取原先的一個頂點與相鄰的兩個邊上的節(jié)點與四邊形中點形成四邊形單元,如表2所示。

表2 四邊形四節(jié)點單元節(jié)點關(guān)聯(lián)矩陣

步驟5在將單元編號如表1向表2轉(zhuǎn)化的過程中,需要注意形成單元的節(jié)點關(guān)聯(lián)矩陣時,節(jié)點序號的編號需按照逆時針順序排序,否則代入插值函數(shù)時計算會出現(xiàn)錯誤。在圖3(d)中,再取每個四邊形的每條邊的中點,在四邊形四節(jié)點基礎(chǔ)上形成新的節(jié)點編號與相應(yīng)節(jié)點坐標,從而得到四邊形八節(jié)點單元,修改單元信息與節(jié)點信息如表3所示,即可完成三角形網(wǎng)格到八節(jié)點四邊形網(wǎng)格的轉(zhuǎn)變。

表3 八節(jié)點四邊形單元編號

步驟6在使用所提方法得到四邊形八節(jié)點網(wǎng)格后,即可以使用無量綱化最小二乘有限元法計算流場速度,其流場的計算流程如圖4所示。

圖4 基于三角形網(wǎng)格的油流速度計算流程圖

3 方法有效性驗證

驗證模型為如圖5所示的頂蓋驅(qū)動方腔模型,其邊長為0.1 m,頂蓋水平方向速度為1 m/s,豎直方向速度為0 m/s,剩余三邊水平速度與豎直速度均設(shè)置為0 m/s,在底邊中點設(shè)置壓力參考點為0 Pa。

u為水平方向速度;v為豎直方向速度;P為壓力

對圖5模型分別采用四邊形和三角形剖分,分別得到1 225個四邊形網(wǎng)格和1 210個三角形,再處理成四邊形網(wǎng)格,如圖6所示。對于三角形網(wǎng)格采用第二部分介紹的方法進行二次處理,然后分別采用有量綱和無量綱最小二乘有限元法對流體場進行計算。有量綱最小二乘法計算時,流體密度為1 000 kg/m3,動力黏度為0.25 Pa·s。無量綱最小二乘有限元計算時,可計算得到雷諾數(shù)為400。計算獲得的流體場云圖如圖7所示。

圖6 方腔剖分網(wǎng)格

圖7 頂蓋驅(qū)動方腔模型速度云圖

從圖7速度云圖可以看出,基于三角形網(wǎng)格的方法與四邊形結(jié)果基本一致,僅存在微小差別。為了進一步驗證方法的有效性,又取中心線上各點的速度,x=0.05 m,繪制在圖8中,同時給出了采用Fluent軟件的仿真結(jié)果。

圖8 中心線速度對比

從圖8中可以看出,所提方法有無量綱中心線上的流場速度計算結(jié)果與四邊形網(wǎng)格的結(jié)果基本重合。故只需將本文方法與Fluent中心線水平、豎直方向速度對比,如圖9所示。

圖9 中心線分向速度對比

經(jīng)過比較,驗證了所提基于三角形網(wǎng)格再處理成四邊形網(wǎng)格的方法的適用性。

4 變壓器單分區(qū)模型油流速度計算

4.1 變壓器單分區(qū)模型,計算步驟及邊界條件

由于網(wǎng)格剖分原因,文獻[15]在流-熱耦合分析時將靜電環(huán)和角環(huán)等復雜結(jié)構(gòu)假設(shè)成“橫平豎直”的結(jié)構(gòu),對計算結(jié)果可能帶來一定的影響。為此,采用三角形網(wǎng)格對文獻[15]中換流變壓器模型的單分區(qū)模型進行剖分,然后分析其油流分布。

圖10為換流變壓器繞組的單分區(qū)模型。模型大小為108 mm×222 mm,線餅為88 mm×10 mm,油流入口和出口都為10 mm,絕緣筒寬為2 mm,左側(cè)絕緣紙筒長為96 mm,右側(cè)絕緣紙筒長為98 mm,兩個絕緣紙筒與油墊圈的連接內(nèi)側(cè)圓弧半徑均為10 mm,外側(cè)圓弧半徑均為12 mm。線餅與最下層油墊圈的距離為2 mm,其余每個線餅間的距離為5 mm,線餅與角環(huán)的距離也為5 mm,角環(huán)厚度為2 mm,角環(huán)下側(cè)矩形長為74 mm,角環(huán)內(nèi)側(cè)小圓弧半徑為2 mm,角環(huán)外側(cè)小圓弧半徑為4 mm,角環(huán)內(nèi)側(cè)大圓弧半徑為15 mm,角環(huán)外側(cè)大圓弧半徑為17 mm,角環(huán)上側(cè)矩形長48 mm,角環(huán)距油墊圈的距離為2 mm,角環(huán)上側(cè)的墊圈長為68 mm,最上層的墊圈長為88 mm。由于絕緣筒與油墊圈僅僅只是形狀對于油流速度有影響,且暫不考慮其散熱影響。因此可以把模型簡化成圖10右側(cè)的計算模型。

圖10 換流變壓器局部單分區(qū)模型

對圖10計算模型采用三角形剖分,網(wǎng)格如圖11所示,然后采用所提網(wǎng)格再處理方法進行處理。

圖11 換流變壓器單分區(qū)油流區(qū)域網(wǎng)格剖分

采用無量綱最小二乘有限元法計算時,油流入口邊界條件設(shè)置為u=0.03 m/s,油流出口為壓力耦合邊界,繞組、角環(huán)、油墊圈等固體壁面無滑移邊界條件為u=v=0 m/s。變壓器油的密度為863.76 kg/m3,動力黏度為0.028 83 Pa·s。

4.2 換流變壓器油流速度分布結(jié)果

本文方法得到的流體場計算結(jié)果如圖12(a)所示,為了驗證方法的有效性,圖12(b)也給出了仿真軟件Fluent的云圖。

圖12 換流變壓器單分區(qū)內(nèi)部油流分布云圖

由圖12可知,本文方法的結(jié)果與Fluent結(jié)果基本一致,僅在第二油道出口與繞組到角環(huán)處存在著微小差距。根據(jù)流體力學的知識,可以知道速度在這種情況下不可能突變,所以在第二油道處的速度本文方法計算結(jié)果更加準確。油道入口與油道出口附近油流速度相同,滿足能量守恒定律。在圖12(a)速度云圖中,運用本文方法也未出現(xiàn)“死流”現(xiàn)象,驗證了本文方法的有效性。

各油道在x=0.054 m處中心點速度如圖13所示,可以看出,本文方法數(shù)值解基本與Fluent結(jié)果變化趨勢一致。與Fluent結(jié)果相比,本文方法計算的中心點速度在前10個和第12個油道略微偏大,在第11個油道和最后一個油道本文方法油道中心點速度偏小。在最后一個油道中心點速度偏差最大,兩者偏差4.6%。

圖13 各油道中心點速度分布

進一步提取差距最大的第15油道中心線速度并與Fluent計算結(jié)果對比如圖14所示。不難發(fā)現(xiàn),兩者僅僅只是在速度中心點附近存在微小差距,其他位置速度基本一致。造成這些差距的原因主要有兩點,第一是所用方法的不同,Fluent是基于有限體積法,而本文方法是基于最小二乘有限元法,這會造成一定誤差。第二是由于二者采用的網(wǎng)格在節(jié)點和數(shù)量上存在差異,也會造成誤差。綜上所述,驗證了本文方法的有效性。

圖14 第15油道中心線速度分布

5 結(jié)論

提出了一種基于三角形網(wǎng)格再處理為四邊形網(wǎng)格計算二維流體場的方法,并用該方法計算了方腔模型有無量綱速度分布及變壓器單分區(qū)復雜結(jié)構(gòu)的油流速度分布,與商用軟件Fluent結(jié)果對比。分析結(jié)果表明。

(1)通過方腔模型的對比發(fā)現(xiàn),本文方法在有量綱與無量綱兩種情況下計算結(jié)果一致,兩種方法下均適用。

(2)本文方法能正確計算含角環(huán)等復雜部件的變壓器模型的油流分布。

本文方法可以對曲面結(jié)構(gòu)具有良好的剖分特性,在計算含復雜部件的電力裝備的流場時本文所提方法有良好的適用性。雖然是在三角形網(wǎng)格上進行驗證,但是對于含有四邊形和三角形的混合網(wǎng)格,本文方法也適用。