光伏變電站用氧化鋅避雷器多物理場耦合分析

張天翊，梁志堅，朱瑞

(1. 廣西大學 電氣工程學院，廣西 南寧 530004；2. 國網能源研究院有限公司，北京 102209)

0 引言

光伏系統及變電站作為新能源發電、系統能量轉換與傳輸的重要組成部分，受到了極大關注，其安全及可靠性顯得尤為重要[1-2]。然而，光伏變電站通常安裝在露天區域，尤其是在雷電災害頻發的地區，雷擊已經成為威脅系統設備的主要安全隱患之一[3]，因此需要對光伏變電站采取防雷措施加以保護。近年來，氧化鋅避雷器以其殘壓低、無續流、動作時延小、通流容量大等優點，在高壓及特高壓輸變電系統中得到廣泛應用[4]。但是，隨著直流輸電工程的快速發展與線路電壓等級的提高，避雷器出現了損耗較大、溫升較高等一系列問題[5]。而溫升的大小和分布對氧化鋅閥片及避雷器的使用壽命有重大影響，關乎整個電力系統安全運行[6]。因此，研究避雷器溫度場分布，特別是在老化等故障情況下，就顯得尤為重要[7]。

目前，很多學者對電力設備多物理場進行分析，主要的方法總結為：試驗法、解析法和有限元法等[8]。其中，試驗方法需要加工一個實際模型，導致成本過高。此外，一些特殊工況的實驗條件也過于苛刻，因此，不適用于模型的初步研究。文獻[9]提出場路耦合法對500 kV避雷器進行溫度場計算。光伏變電站配備的220 kV避雷器采用的是兩節元件結構，在高、低壓段間的電場分布發生了躍變現象。因此采用常規電磁解析方法計算損耗，再結合傳熱學方法對模型溫度場進行分析的方法，并不是溫度場分析的最優選擇。此外，場路耦合方法雖然從物理概念上描述比較清晰，但是無法考慮流熱耦合計算中的集膚效應、渦流效應和流體流動等影響，難以全面而準確地描述各個場的具體分布[10]。目前，有限元法采用對稱正定的稀疏矩陣，可以在綜合考慮對流及輻射換熱效應的基礎上，直接求解流固耦合方程[11]，因此能夠比較全面地對電力設備的多物理場進行耦合分析。文獻[12]采用有限元法計算10 kV氧化鋅避雷器的溫升，但是計算的是二維簡化模型，忽略了熱量在垂直方向的傳導過程。文獻[13]采用商業軟件計算動車組車頂避雷器的溫升，但是忽略了熱輻射和空氣對流換熱對溫度計算的影響。雖然有限元法得到了成功的應用，但是該方法是建立在區域變分原理和泛函計算的基礎上，一旦網格剖分結束后就不易進行修正，因此不具備較強的自適應能力[14]。特別是在前處理過程中的網格劃分十分繁瑣，會占用大量時間。同時，網格扭曲或畸變現象導致計算精度差、收斂困難等問題也亟待解決[15]。因此迫切需要一種能保障計算收斂速度，滿足計算精度，并且不依賴網格劃分的新方法。

近年來，一種只需節點信息，無須劃分網格即可構造高階連續函數，具有較高計算精度的無網格法得到廣泛應用[16-17]。但是，無網格法采用徑向基函數法構造函數的近似表達式，導致計算得到的一般是節點參數而非真實節點值。因此函數不滿足Kronecker特性，不能像有限元法一樣便捷地施加邊界條件。這也使得針對溫度場計算邊界的處理十分麻煩[18]。而有限元法處理邊界問題比較方便，因此將有限元法與無網格法結合起來，可以較好地解決電力設備的多物理場耦合計算問題。但是，在有限元區和無網格交界面處的節點值是真實的，而無網格法的節點值是虛假的，無法直接匹配。因此，文獻[19]引入由有限元與無網格法形函數組成的界面單元，雖然保證了方程求解的連續性，但忽略了交界面函數導數的連續性，而對函數及其空間導數的近似計算恰恰是無網格法求解偏微分方程的關鍵所在。目前，考慮魯棒性的針對電力設備有限元與無網格法耦合計算的文獻還未見公開報道。

本論文研究如何將無網格方法應用到電力設備的流-熱耦合計算中，重點研究有限元求解的流體場與無網格求解的溫度場在交接界面的數據傳遞。首先建立1臺經過簡化處理的某光伏變電站220 kV氧化鋅避雷器三維流體場模型。采用商業有限元軟件Ansys分別對其額定工作狀態和老化故障下的流速分布進行計算，得到各壁面的對流換熱系數。其次，通過在有限元法與無網格法交界面之間構造過渡區域，采用徑向基函數法構造無網格法的形函數，推導出基于無網格法的傳熱學控制方程。在此基礎上，建立了2種工況下的無網格溫度場模型，得到了相應的溫度分布。然后，將無網格-有限元法間接耦合計算結果與純粹有限元流體-溫度直接耦合仿真結果進行對比。最后通過實驗，在證明了本文耦合計算方法準確性的同時，還能發現無網格計算的速度遠高于有限元法。本文所述方法也可以為其他電力設備多物理場耦合分析提供一定的參考。

1 避雷器結構參數與耦合計算流程

1.1 模型主要參數

圖1 避雷器三維仿真模型Fig. 1 Three dimensional simulation model of surge arrester

表1 避雷器主要參數Table 1 Major parameters of surge arrester

表2 材料主要參數Table 2 Major parameters of materials

為了便于分析，對避雷器做以下假設。

（1）在絕緣套管與空氣之間為大空間輻射換熱，因此將氣體和固體交界面設置為無滑移邊界[20]；

（2）流體被視為不可壓縮氣體，忽略浮力影響，同時忽略溫度對空氣粘度的影響，計算得到的避雷器壁面的對流換熱系數均為平均值；

（3）假設初始溫度為室溫，忽略接觸熱阻，所有介質均可視為各向同性介質。

1.2 流體-溫升耦合計算流程

無網格法與有限元法結合的氧化鋅避雷器多物理場耦合分析流程如圖2所示。首先，建立避雷器三維流體場有限元模型，并有如下簡化處理：忽略避雷器中的一些細微結構，如不影響溫升計算的退刀槽、螺紋孔、螺栓等；另外將曲率較小弧形的導角均按照直角處理；其次，對模型進行網格劃分，并進行求解；然后，將流體場求解出的對流換熱系數作為邊界，采用有限元法對避雷器進行共軛傳熱計算；最后，通過在有限元與無網格計算區域中間建立過渡函數，實現流體場與溫度場的數據傳遞。在同樣的激勵與邊界條件的基礎上，利用無網格法再次進行溫度場計算，得出不同工況下避雷器溫升的分布。

圖2 避雷器多物理場耦合計算流程Fig. 2 Process of fluid-thermal coupling analysis

2 避雷器多物理場耦合數學模型

2.1 有限元法流體場控制方程

在避雷器流體場有限元計算中，主要的散熱方式為自然冷卻，因此氣體流通模式主要為層流，滿足的流體場控制方程[21]如下。

2.2 有限元法溫度場控制方程

在溫度場計算中，同時考慮了熱傳導、對流傳熱和熱輻射的影響，其控制方程與邊界條件[23]為

借助準滑移面方法研究RT位移模式非極限狀態被動土壓力，結果表明：繞墻頂轉動擋土墻非極限狀態被動土壓力為凹曲線分布，越靠近墻底，被動土壓力增大越快；轉動越大，被動土壓力越大，其合力作用點越低，合力作用點在墻高的下三分點以下。

氧化鋅避雷器運行時會通過套管向周圍不斷輻射熱量，計算公式[24]表示為

式中：εt為表面輻射率，對于絕緣套管表面取為0.95，而對于氧化鋅閥片則為0.71；σ為斯忒藩-玻耳茲曼(Stefan-Boltzmann)常數，文中取為5.7×10-8W/(m2?K4)；Sto為壁面輻射面積。

最終溫度場有限元離散形式為

式中：KFEA為剛度矩陣；N為由四面體單元組成的形函數。

2.3 無網格法溫度場控制方程

2.4 有限元法與無網格法耦合

3 多物理場耦合計算結果

3.1 流體-溫度場有限元計算

圖3為避雷器的網格剖分，可看出網格質量良好，符合計算精度要求?？偟墓濣c數為1 815 403，網格數為8 841 828。因此，網格劃分需要占用大量的計算時間和資源，這也是本文所述無網格法優勢所在。

圖3 流體計算網格Fig. 3 Meshes of fluid calculation

額定運行與上節閥片老化故障時避雷器流體場粒子分布如圖4所示?？梢钥闯?，在2種工況下，流體的分布區別不大，最大值分別為0.18 m/s和0.181 m/s，都出現在套管外部正上方。

圖4 避雷器流速分布Fig. 4 Flow velocity distribution of surge arrester

求出的對流換熱系數可以用以下方程表示。

對于避雷器側壁面，可表示為

對于避雷器頂面和底面，則分別為

式中：Cf為阻力系數；Ra為瑞利數。

圖5為避雷器額定運行與上節閥片老化故障時的絕緣套管表面溫度分布情況。從圖5 a)可以看出，正常情況下套管表面最大溫升為21.65 °C，主要位于中間法蘭與上節底部外套連接處。由于受到重力加速度對氣體流動的影響，溫度整體呈現出在豎直方向的梯度變化。圖5 b)為上節老化故障情況下的溫度分布，與正常工況對比可以發現，下節套管表面溫升要低一些，溫升最大為22.25 °C。

圖5 基于有限元法的避雷器絕緣套管溫升分布Fig. 5 Thermal distribution of arrester insulating sleeve based on FEA

圖6為氧化鋅閥片額定運行及老化2種狀態下的溫度分布情況?？梢钥闯?，在正常工況下，閥片各節均呈現出中間溫度高而兩端溫度低的分布趨勢，最高溫度為31.42 °C。主要是因為避雷器內部閥片柱兩端的金屬連接件散熱效果較好，所以溫度值較低。當閥片老化故障時，最大溫升為38.45 °C?？梢园l現老化的閥片溫升要明顯降低，故需要特別關注下節閥片的溫升變化。

圖6 基于有限元法的氧化鋅閥片溫升分布Fig. 6 Thermal distribution of ZnO valve plate based on FEA

出現這種現象的主要原因是當氧化鋅閥片在老化時，電阻片吸收了大量的能量。在較大的熱載荷作用下，內部晶界層中的離子活動非?；钴S，增大了離子遷移率，導致閥片內部肖特基勢壘壁高度降低，因此對過電壓產生了抑制效果，有效降低了損耗[27]。此外，在長時間的大電流作用下，閥片會受到較大電磁力和熱應力的影響，導致壓敏電阻片內部結構發生變化，產生明顯的極性效應。因此，損耗分布發生較大變化引起溫升分布不均勻。

3.2 無網格溫度場計算

圖7為避雷器額定運行與上節閥片老化故障時絕緣套管表面溫度分布情況。對比圖5，可以看出二者溫升分布基本相同，最大溫升分別為22.22 °C 和 23.25 °C，相應的誤差分別為 2.6% 和4.5%。

圖7 基于無網格法的絕緣套管溫升分布Fig. 7 Thermal distribution of insulating sleeve based on meshless method

圖8為氧化鋅閥片正常及老化2種狀態下的溫度分布情況。對比圖6，可以看出計算出的溫升分布基本相同。其中，無網格法計算出的閥片正常狀態工況的最大溫升為31.83 °C，誤差為1.3%。而老化后的最大溫升為38.69 °C，誤差為0.6%，可以證明所用無網格計算方法有足夠的精度。

圖8 基于無網格法的氧化鋅閥片溫升分布Fig. 8 Thermal distribution of ZnO valve plate based on meshless method

表3為流體-溫度耦合計算過程中各個環節的用時統計，可以發現，有限元法由于產生大量網格，計算時間遠超無網格法。對比可以發現采用有限元法計算溫度場消耗的時間大約是采用無網格法的210倍。

表3 計算時間統計Table 3 3 Calculation time

4 溫升實驗

為了驗證避雷器多物理場耦合計算的準確性，按照符合國家標準GB18802.1—2011《低壓電涌保護器(SPD)第1部分：低壓配電系統的電涌保護器性能要求和試驗方法》的規定來確定裝備的溫升。溫度由mv2000型溫度傳感器來監測，檢測點的安裝位置如圖1所示。避雷器持續工作20 h，采樣時間為10 min，測量的結果如圖9所示。由圖9可知，各仿真結果與試驗結果基本吻合，滿足精度要求。structure surge arrester[J]. High Voltage Engineering, 2012, 38(8):2129-2136.

圖9 溫升實驗結果Fig. 9 Experimental results of temperature rise

5 結論

本文提出了一種針對光伏變電站用氧化鋅避雷器溫升計算的新方法，構造了斜坡函數，將有限元與無網格法結合，綜合了2種方法的優點。得出如下主要結論。

（1）描述了無網格-有限元耦合法的基本原理，根據溫度場的控制方程，推導出了綜合無網格法與有限元的溫度場矩陣方程；

（2）用有限元法分別計算避雷器額定和老化工況下的流速和溫升分布情況，結果顯示，額定運行時，避雷器溫度場分布較為均勻，而在老化故障下，氧化鋅閥片溫升要明顯降低；

（3）采用無網格-有限元耦合法分析溫度場問題時，不僅能夠保障計算精度，還能夠極大地提高計算效率。此種方法對配置在光伏變電站內的電氣設備數值計算具有一定的參考價值及意義。