油浸式變壓器餅式繞組電磁-溫度-流體場雙向耦合數(shù)值建模與分析

劉慧娟，杜晉文，張振洋，劉 威

(北京交通大學(xué) 電氣工程學(xué)院，北京 100044)

配電變壓器作為保障電網(wǎng)正常工作的關(guān)鍵設(shè)備，其性能和壽命一直以來受到變壓器制造商與應(yīng)用者的廣泛關(guān)注，而變壓器的絕緣材料老化與熱點溫度直接相關(guān)[1],其老化速度與熱點溫度近似成正比例關(guān)系，因此更加合理地設(shè)計與優(yōu)化變壓器冷卻系統(tǒng)變得尤為重要.

目前國內(nèi)外研究學(xué)者針對在工程中應(yīng)用較為廣泛的餅式繞組溫升特性的研究，主要以電磁場、流體場與溫度場順序單向耦合或者電磁場與溫度場單向耦合而流體場與溫度場雙向耦合兩種方式來展開.Zhang等[2]不考慮電磁場、溫度場對流體場的影響作用，只通過流體場分析餅式繞組內(nèi)部變壓器油的壓強與速度分布情況，考慮回流和沿徑向截面積的不一致等因素分析繞組冷卻油道對變壓器油流動的影響，建立了計算更加快速的流體網(wǎng)絡(luò)模型；廖才波等[3-5]針對油浸式三相變壓器與立體卷鐵芯式變壓器的溫升與油流分布，忽略餅式繞組水平油道冷卻作用與導(dǎo)線匝間絕緣等因素，對整塊低壓繞組進行建模，基于三維電磁-溫度-流體場之間的間接耦合計算分析了變壓器結(jié)構(gòu)內(nèi)部的熱特性；謝裕清等[6-7]未考慮電磁場中損耗的計算，主要針對餅式繞組提出采用最小二乘有限元法求解流體場，將計算所得繞組中油的速度分布單向耦合到溫度場，運用迎風(fēng)有限元法分析求解餅式繞組中相互耦合的溫度場，并基于COMSOL多物理場仿真軟件就影響其溫升分布的重要因素進行了參數(shù)化分析；周利軍等[8-9]忽略絕緣材料對傳熱以及溫度場對電磁場中導(dǎo)線電阻率的影響，在溫度場與流體場中利用多面體網(wǎng)格剖分對簡化的變壓器三維計算模型溫升特性進行了研究分析,并通過與實驗對比，驗證了多面體網(wǎng)格剖分求解的有效性并指出多面體網(wǎng)格的計算優(yōu)勢；Torriano等[10-11]忽略溫度場對電磁場的影響作用，針對變壓器餅式繞組溫升的六種數(shù)值計算模型進行了比較分析并著重分析了油流入口溫度邊界條件對繞組溫升分布的影響，并通過二維與三維仿真的對比分析提出了針對二維數(shù)值計算的修正方法.

綜上所述，針對油浸式變壓器餅式繞組溫升特性的研究，考慮繞組匝間絕緣，對其數(shù)值仿真過程中電磁場、流體場和溫度場之間不同耦合程度的比較研究在文獻中較少報道.

本文作者通過分析電磁場、溫度場與流體場三場的控制方程及三場間耦合變量的關(guān)系，提出通過緊耦合求解溫度場與流體場，而在電磁場與溫度場間以松耦合的方式計算.該方法旨在實現(xiàn)電磁模型與流固耦合傳熱模型的雙向耦合仿真，交互共享三場間耦合的變量參數(shù)，實現(xiàn)電磁場、溫度場與流體場間的完全耦合計算.然后針對一臺實際的變壓器餅式繞組，對建模過程中影響數(shù)值仿真精確性的因素展開分析，確定變壓器油的對流方式、邊界條件和離散計算的算法，且以Boussinesq假設(shè)計及自然對流作用，測試導(dǎo)線匝間絕緣對傳熱計算影響的程度，并對其剖分網(wǎng)格的敏感性進行分析.仿真計算該變壓器繞組的溫度場分布和變壓器油的運動特性，并通過比較三種不同程度的耦合仿真方式，結(jié)合三場間耦合變量的關(guān)系式，對仿真計算結(jié)果進行了合理性分析驗證.

1 多物理場耦合計算方法

1.1 電磁場控制方程

變壓器低壓繞組的損耗依賴于變壓器內(nèi)部的磁場分布，根據(jù)電磁場理論[12]，基于磁矢量位的電磁場控制方程式可表示為

(1)

式中：Js為繞組的源電流密度，A/m2；A為矢量磁位，Wb/m；t為時間，s；μ為磁介質(zhì)的磁導(dǎo)率，H/m；為繞組的電導(dǎo)率，S/m.

繞組的損耗P繞組為電阻損耗和渦流損耗之和，表示為總電流密度對繞組體積的積分形式

(2)

式中：J為繞組中的總電流密度，A/m2;V為繞組的體積，m3.

1.2 流體場-溫度場控制方程

變壓器油為不可壓縮的牛頓型流體，在變壓器油箱內(nèi)部的流動主要受質(zhì)量守恒定律與動量守恒定律的制約，其控制方程可表示為

(3)

(4)

式中：ρoil、μoil分別為變壓器油的密度、動力黏度，單位分別為kg/m3、N·S/m；Uoil、p為控制方程求解的未知量，分別表示油流的流速、壓力，單位分別為m/s、N；g為重力加速度，m·s-2；ρref為變壓器油密度變化的參考值，本文中取為22 ℃時對應(yīng)的變壓器油密度值.

繞組內(nèi)部熱量通過熱傳導(dǎo)的形式傳到繞組表面，在繞組表面與變壓器油發(fā)生對流換熱，換熱過程中遵循能量守恒定律，根據(jù)變壓器的傳熱機制，低壓繞組為固體區(qū)域，故將繞組內(nèi)部區(qū)域的傳熱控制方程中流速項U設(shè)為零.同理，由于變壓器油內(nèi)部不產(chǎn)熱，即不含熱源，故將變壓器油的熱源密度SE設(shè)為零，則整個計算域內(nèi)的傳熱控制方程可統(tǒng)一表示為

(5)

式中：ρ、cp、λ分別為傳熱區(qū)域?qū)?yīng)的密度、定壓比熱容、導(dǎo)熱系數(shù)，單位分別為kg/m3、J/(kg·K)、W/(m·K)；T為對應(yīng)區(qū)域的溫度，K；SE為熱源密度，為式(2)繞組損耗P繞組對繞組體積的平均值.

1.3 電磁場-流體場-溫度場耦合計算

由于特征區(qū)域內(nèi)變壓器磁場的漏磁場較弱，且變壓器的導(dǎo)線采用多根并聯(lián)繞制的方式，因此不考慮繞組中的渦流損耗.變壓器繞組采用銅導(dǎo)線，在仿真過程中對導(dǎo)線匝間絕緣進行單獨建模，因此確定其繞組損耗隨溫度的變化關(guān)系[13]為

(6)

式中：T繞組為繞組的溫度;P0為溫度為T0時繞組的損耗值.

變壓器油采用Nynas Lynx環(huán)烷烴礦物油[14]，熱膨脹系數(shù)β為0.000 643，單位為1/℃，密度ρoil、導(dǎo)熱系數(shù)λoil、比熱容cpoil、動力黏度μoil隨溫度的變化關(guān)系為

(7)

式中：Toil為變壓器油的溫度，℃；變壓器油的物性參數(shù)除比熱容之外，密度、導(dǎo)熱系數(shù)、動力黏度系數(shù)都隨著溫度的升高而減小.

電磁場、溫度場和流體場的控制方程均為高度非線性的方程，直接聯(lián)立各場域方程耦合求解計算成本較高且更容易發(fā)散，因此為了更加全面地考慮各場域之間的耦合作用，同時保證整個耦合計算更加容易收斂，本文在溫度場與流體場之間的耦合基于有限體積法求解，電磁場運用有限元法計算繞組的損耗，并在一個電周期內(nèi)對單匝導(dǎo)線損耗取平均值傳遞到溫度場，作為流體場與溫度場計算的熱載荷，流體場與溫度場緊耦合計算每匝導(dǎo)線對應(yīng)的溫度及變壓器油流的分布特性，并將每匝導(dǎo)線對應(yīng)的溫度值取平均反饋至電磁場(松耦合)，三場反復(fù)交替迭代計算直至收斂.

耦合迭代計算過程中，溫度場中繞組的溫度由式(6)不斷修正繞組的電阻值，并反饋至電磁場，而變壓器油的密度、比熱容、導(dǎo)熱系數(shù)和動力黏度系數(shù)由式(7)修正，并反饋至流體場.具體的耦合流程如圖1所示.

2 計算模型構(gòu)建

2.1 物理模型

由于低壓繞組散熱條件較高壓繞組差等因素限制，油浸式三相變壓器的熱點溫度總是位于低壓繞組中上部[4-5].本文以66MVA-225/26.4 kV ONAN/ONAF配電變壓器中的低壓繞組[10]作為研究對象，對其多物理場間耦合作用的不同簡化方式展開研究.該變壓器低壓繞組結(jié)構(gòu)如圖2所示，采用 “之字形”冷卻結(jié)構(gòu)，由78段線餅組成，每段線餅由18匝扁銅導(dǎo)線沿輻向繞制而成，低壓繞組內(nèi)側(cè)與對稱軸間的距離為316.2 mm，左右側(cè)豎直油道寬均為6.4 mm，水平油道寬為4.1 mm，由墊圈將整塊餅式繞組分隔成5個區(qū)域，第一區(qū)域由2段線餅構(gòu)成，其他分區(qū)均由19段線餅構(gòu)成，上一分區(qū)的油流出口作為下一分區(qū)的油流入口，繞組內(nèi)外側(cè)均由絕緣筒所圍，油流從低壓繞組底部進入，通過墊圈及絕緣筒引導(dǎo)變壓器油流過水平、豎直油道，與餅式繞組發(fā)生熱交換.文獻[10]中通過比較分析二維與三維的數(shù)值計算結(jié)果對二維數(shù)值模型進行修正，驗證了二維數(shù)值計算的工程適用性，圖2為繞組的特征流體存在的區(qū)域，二維數(shù)值計算占用較少的計算資源，因此根據(jù)低壓繞組分區(qū)的幾何對稱性，通過二維軸對稱數(shù)值計算探究多場之間的耦合作用.

2.2 數(shù)值仿真模型精確性分析

2.2.1 網(wǎng)格剖分及對流方式

數(shù)值分析方法計算的準確與否和網(wǎng)格剖分的質(zhì)量直接相關(guān)，為了保證數(shù)值計算精度，通過特征區(qū)域內(nèi)最高溫度與對應(yīng)剖分網(wǎng)格數(shù)量之間的關(guān)系驗證網(wǎng)格的收斂性[15]，如圖3所示.

由圖3可知，對比網(wǎng)格單元數(shù)為539與394 415時計算結(jié)果，二者計算所得最高溫度相對差值高達38.84%，印證了網(wǎng)格可靠性驗證的必要性.當網(wǎng)格數(shù)量大于100 000時，最高溫度的計算值變化小于0.24%，結(jié)合實際計算機資源，最終選擇網(wǎng)格單元數(shù)為108 686，節(jié)點數(shù)為135 492進行計算.

變壓器油在近壁面處存在黏滯效應(yīng)[16]，為了捕捉變壓器油邊界層的溫度梯度與速度分布梯度特征，在所建立的數(shù)值計算模型中對接近壁面區(qū)域的邊界層加密剖分，第一層網(wǎng)格高度選為0.001 mm，以1.1的比率向油道中心區(qū)域逐漸增大，共剖分10層邊界網(wǎng)格，結(jié)果如圖4所示.

變壓器油的黏度較大而流速很小，根據(jù)流體力學(xué)、傳熱學(xué)理論中雷諾數(shù)Re與理查森數(shù)Ri的實際物理意義以及對應(yīng)的計算公式[16]，確定其流動方式為層流；對流換熱方式為以強迫對流換熱主導(dǎo)的混合對流傳熱.變壓器油的熱膨脹系數(shù)β較小，即在該變壓器的運行過程中始終滿足β(T-T0)?1，結(jié)合Boussinesq假設(shè)的適用范圍，故采用Boussinesq假設(shè)來考慮自然對流作用.

2.2.2 溫度場-電磁場間迭代數(shù)據(jù)映射基本單元

圖5(a)為以單匝導(dǎo)線且包括匝間絕緣作為溫度場與電磁場間數(shù)據(jù)映射的基本單元時繞組的溫度分布云圖；以整塊餅式繞組作為溫度場與電磁場間數(shù)據(jù)映射的基本單元，計算所得繞組的溫度分布如圖5(b).仿真過程中二者除數(shù)據(jù)映射基本單元不一樣之外，其他求解條件均相同.

由圖5可知，兩種模型計算所得的繞組熱點位置沿軸向的位置基本相同，而單匝繞組模型比整塊繞組模型計算所得繞組溫升高約19.64%，且以單匝導(dǎo)線作為迭代計算數(shù)據(jù)映射的基本單元可以更加準確地計算出熱點位置的徑向坐標，在實際應(yīng)用中可以更加有效地指導(dǎo)監(jiān)測所用傳感器的安裝，因此采用整塊繞組作為多物理場間迭代計算數(shù)據(jù)映射基本單元是不合理的.由于“之字形”冷卻油道的流向使得餅式繞組沿徑向溫度分布不均勻，采用單匝導(dǎo)線建?？梢愿忧逦乇憩F(xiàn)出繞組溫度沿徑向的溫度梯度，且可以考慮匝間絕緣對傳熱的影響.因此對繞組單匝導(dǎo)線分別進行建模，并考慮繞組導(dǎo)線匝間絕緣對傳熱的影響，計算模型如圖6所示.

2.2.3 計算參數(shù)及邊界條件

銅導(dǎo)線迭代計算的初始電導(dǎo)率值取為22 ℃時對應(yīng)的電導(dǎo)率，其他物性參數(shù)受溫度的影響忽略不計，具體參數(shù)如表1所示[11].

表1 材料物性參數(shù)

關(guān)于計算區(qū)域中邊界條件的設(shè)定為

1)變壓器散熱片等結(jié)構(gòu)對低壓繞組冷卻的作用，由繞組入口邊界條件等效處理，冷卻油垂直進入繞組內(nèi)部，確定油流入口邊界條件為流量入口邊界，流量為0.002 167 kg/s，入口處溫度為46.7 ℃[10]；油流出口邊界條件為自然壓力出口邊界；旋轉(zhuǎn)軸為軸對稱邊界條件.

2)墊圈與絕緣筒的導(dǎo)熱性能較差，均近似為絕熱邊界；變壓器冷卻油與繞組、絕緣筒及墊圈之間接觸面采用流固耦合無滑移邊界.

電磁場采用有限元法求解繞組損耗，而溫度場與流體場應(yīng)用有限體積法求解繞組穩(wěn)態(tài)溫度以及變壓器油的流動速度.流體場中壓力與速度之間的離散采用Coupled算法，同時求解式(3)與式(4)；壓力、動量與能量項均采用二階逆風(fēng)算法進行空間離散，防止引起數(shù)值振蕩，提高計算精度；所有計算中涉及物理量的殘差均為10-8.通過驗證繞組與變壓器油之間的流固耦合交界面熱流量傳遞的連續(xù)性以及整個區(qū)域內(nèi)變壓器油的質(zhì)量守恒驗證迭代計算結(jié)果的收斂性.

3 仿真結(jié)果分析

將溫度場與流體場間緊耦合計算，而電磁場與溫度場之間松耦合計算，實現(xiàn)電磁-溫度-流體場三場完全耦合的計算的方法標記為EFT，由該方法迭代計算繞組溫升及變壓器油速度分布，計算過程中能量與速度變化的殘差曲線如圖7所示，與電磁場之間共迭代6次.

三場耦合計算達到收斂共迭代2 368步，其中電磁場與溫度場間以松耦合方式計算，迭代過程中計算的收斂速度隨著迭代次數(shù)的增加而加快.圖8為迭代計算過程中繞組總損耗值的變化曲線.

繞組初始的電導(dǎo)率選擇為繞組22 ℃時的值，根據(jù)式(6)可得，在忽略繞組渦流損耗的情況下，繞組損耗與繞組溫度成線性關(guān)系，因此在第一次迭代計算過程雙向耦合對繞組損耗值的修正幅度最大.設(shè)定兩次迭代計算的誤差小于5%，迭代6次計算即收斂，最終計算所得繞組的總損耗值收斂于1 437.1 W.

3.1 特征區(qū)域溫升分布

將同時忽略溫度場對電磁場與流體場中材料參數(shù)的影響，三場單向耦合的方式簡記為 FT1；將忽略溫度場對電磁場中繞組電阻的影響，而溫度場與流體場以緊耦合計算的方式簡記為FT2.圖9為通過三種耦合方式計算得到的繞組溫度分布云圖，將繞組從下至上，依次編號為1,2，…，19.

由圖9可知，在3種耦合計算方式下，繞組的溫度分布趨勢基本相同，當變壓器油從繞組內(nèi)側(cè)流入時，入口處與出口處溫度均較低，分區(qū)內(nèi)繞組的最高溫度均出現(xiàn)在中間區(qū)域第10餅沿徑向靠外側(cè)，即在三種耦合方式下計算所得最高溫度的位置是相同的.靠近出口處由于溫度邊界層的作用均表現(xiàn)出“熱斑”現(xiàn)象，油道中心區(qū)域溫度相較兩側(cè)邊界溫度高.

比較圖9中(a)、(b)，由FT1方式計算所得特征區(qū)域內(nèi)繞組最高溫度為60.804 ℃，F(xiàn)T2計算的最高溫度為64.240 ℃，相差5.66%.同理，由圖9(b)、(c)計算結(jié)果可得，F(xiàn)T2與EFT兩種方式計算的繞組最高溫度相差4.78%.為了更進一步說明三場迭代耦合對繞組溫升計算的影響，選擇特征區(qū)域中溫度最高的第10餅繞組，比較繞組水平中心線沿徑向的溫度分布，如圖10所示.

由圖10(a)可知，溫度變化平緩處為繞組單匝導(dǎo)線的所在位置，第10餅繞組沿徑向的最高溫度和最低溫度差值達到17.52 ℃，而由圖10(b)可知，單匝導(dǎo)線沿徑向溫度最大差值僅僅為0.5 ℃，因此選擇單餅繞組作為迭代映射的基本單元，將其溫度取平均值反饋至電磁場計算誤差較大，將單匝導(dǎo)線的溫度取平均值反饋至電磁場進行迭代耦合計算更為合理.

三種不同耦合程度下，計算所得繞組溫度分布沿繞組徑向方向變化趨勢基本一致，繞組的最高溫度均出現(xiàn)在繞組中心靠油流出口一側(cè)，且繞組溫度沿徑向呈先升高后降低趨勢.FT1、FT2與EFT三種方式計算所得的溫度依次增大，且對于溫度越高的區(qū)域，計算差值越大.由于溫度與損耗成線性關(guān)系，通過電磁場與溫度場間的松耦合計算，可以迭代修正兩場中的損耗值與溫度值，對于溫度越高的區(qū)域，修正的幅度越大.由式(7)可得，溫度影響變壓器油的材料參數(shù)，比如該變壓器油的導(dǎo)熱系數(shù)隨著溫度升高而減小，F(xiàn)T1未能考慮溫度升高引起的導(dǎo)熱系數(shù)等材料參數(shù)的變化，因此計算所得的溫度整體較低，而在FT2中溫度場與流體場間緊耦合計算，變壓器油的導(dǎo)熱系數(shù)等材料參數(shù)隨著變壓器油的溫度值實時更正，實現(xiàn)二者間的完全耦合計算.

因此，電磁-溫度-流體場耦合建模方法EFT結(jié)合有限元法與有限體積法各自的計算優(yōu)勢，綜合運用松耦合與緊耦合的計算方式，可以同時考慮溫度的變化對繞組電阻及變壓器油的材料參數(shù)的影響，且以單匝導(dǎo)線作為松耦合計算中數(shù)據(jù)映射的基本單元可以更進一步提高松耦合的計算精度，考慮繞組匝間絕緣材料對傳熱的影響.

3.2 變壓器油流速分布

三種不同耦合仿真計算方法所得自繞組底部沿其豎直中心線變壓器油的速度分布如圖11所示.

由圖11可知，3種耦合方式計算所得油流分布基本一致，靠近出口一側(cè)水平油道的流速略高于靠近入口側(cè)水平油道的流速，而FT1方式由于流體場與溫度場之間單向耦合不考慮溫度場對流體場的影響，流速略低于另兩種耦合方式.三種耦合方式下，變壓器油流速均呈現(xiàn)靠近入口、出口處流速高而繞組中間區(qū)域流速低，結(jié)合繞組的溫度分布趨勢，由于中間區(qū)域冷卻油流速較低而使得該區(qū)域繞組的溫度較高.三場耦合仿真靠近出口處油流的分布云圖如圖12所示，箭頭代表油流的速度方向.

圖12很清晰地表現(xiàn)出變壓器油與繞組接觸的邊界層流動細節(jié)，從繞組表面向油道中心表現(xiàn)出很明顯的速度變化梯度，由于無滑移邊界條件，繞組與變壓器油接觸面油流動速度最低，驗證了邊界層剖分的有效性.變壓器油均由內(nèi)側(cè)豎直油道通過水平油道流向油道，與繞組外側(cè)豎直油道油發(fā)生匯合，在流動過程中與繞組發(fā)生對流換熱，為繞組提供冷卻.

4 結(jié)論

1)數(shù)值仿真模型中網(wǎng)格的數(shù)量直接影響多物理場耦合仿真計算結(jié)果的精度，當網(wǎng)格數(shù)量為539和394 415時，二者計算所得最高溫度差達38.84%.針對本文模型的計算結(jié)果的分析表明，當網(wǎng)格數(shù)量大于100 000時，計算的溫度差小于0.24%，因此，多物理場耦合仿真必須首先進行網(wǎng)格獨立性驗證，才能確保仿真計算的精度.

2)對電磁場、溫度場和流體場間不同耦合程度的仿真結(jié)果進行對比分析，并仿真對比迭代耦合過程中兩種數(shù)據(jù)映射的基本單元，采用單匝繞組模型的繞組溫度比單餅繞組模型的高19.64%，以單餅繞組作為迭代計算數(shù)據(jù)映射基本單元不能預(yù)測繞組溫度沿徑向的分布且單餅繞組溫度沿徑向分布相差較大，將單餅繞組溫度取平均值反饋至電磁場進行迭代耦合計算存在較大誤差，因此在繞組的溫升仿真計算過程中，必須考慮繞組匝間絕緣材料對傳熱的影響，且不能將單餅繞組作為電磁場與溫度場間迭代耦合計算的基本單元.

3)結(jié)合松耦合與緊耦合建立繞組溫升計算多物理場仿真模型，可以考慮繞組匝間絕緣對傳熱的影響和變壓器油及繞組的材料屬性隨溫度的變化，且在5%的殘差指標下，整體計算達到收斂只需迭代6次，有效地保證了計算效率和精度.