基于坐標(biāo)插值映射法的變壓器電磁熱場(chǎng)模擬

張建瓴，周壯廣，可欣榮

(華南農(nóng)業(yè)大學(xué) 工程學(xué)院，廣東 廣州 510642)

變壓器工程領(lǐng)域常常會(huì)遇到多物理場(chǎng)的耦合仿真問(wèn)題，例如磁場(chǎng)和結(jié)構(gòu)場(chǎng)的耦合、電場(chǎng)和磁場(chǎng)的耦合以及溫度場(chǎng)和流場(chǎng)的耦合等。在電磁場(chǎng)計(jì)算時(shí)，網(wǎng)格剖分需著重考慮導(dǎo)磁材料的集膚效應(yīng)，對(duì)熱流場(chǎng)則需著重考慮換熱的邊界層剖分。若對(duì)不同的物理場(chǎng)使用統(tǒng)一的網(wǎng)格剖分模型，就需要同時(shí)兼顧各物理場(chǎng)對(duì)網(wǎng)格形狀和精細(xì)度的要求，勢(shì)必使問(wèn)題復(fù)雜化且浪費(fèi)大量的計(jì)算資源；若不同物理場(chǎng)使用不同的網(wǎng)格剖分模型，則不同物理場(chǎng)之間的網(wǎng)格節(jié)點(diǎn)不同，相互耦合的物理場(chǎng)之間就需要進(jìn)行網(wǎng)格數(shù)據(jù)映射，以實(shí)現(xiàn)耦合變量的準(zhǔn)確傳遞。

目前已發(fā)展出多種針對(duì)這種異構(gòu)網(wǎng)格節(jié)點(diǎn)數(shù)據(jù)映射問(wèn)題的解決方法，如快速殼法[1]、反距離移動(dòng)平均法[2]、快速映射法[3]、徑向基函數(shù)法[4]及混合法[5]等。劉剛等[5-8]利用快速映射法、徑向基函數(shù)法及混合法等方法在大型超高壓變壓器的多物理場(chǎng)耦合問(wèn)題中做了許多工作，但其映射算法和物理場(chǎng)計(jì)算都僅限于二維模型。

前述各種方法中，網(wǎng)格數(shù)據(jù)信息傳遞的載體主要是節(jié)點(diǎn)，其數(shù)據(jù)映射計(jì)算較為復(fù)雜。如快速映射法需要遍歷已知物理場(chǎng)所有單元來(lái)搜索找到主單元；徑向基函數(shù)法需要求解大型方程組，而且將物理量從稀疏網(wǎng)格映射到密集網(wǎng)格時(shí)異常區(qū)域節(jié)點(diǎn)的映射精度較差；快速殼法只適用于節(jié)點(diǎn)在凸殼內(nèi)的情況；反距離移動(dòng)平均法誤差較大；混合法計(jì)算精度高，但步驟繁雜、計(jì)算耗時(shí)多以及工程實(shí)用性不強(qiáng)。

本文將網(wǎng)格數(shù)據(jù)傳輸?shù)妮d體定義為坐標(biāo)，提出了一種新的異構(gòu)網(wǎng)格數(shù)據(jù)映射方法——坐標(biāo)插值映射法，并以S9-1600/10/0.4型油浸式電力變壓器為例，介紹坐標(biāo)插值映射方法將三維電磁場(chǎng)中結(jié)構(gòu)件的電磁損耗密度分布映射到熱流場(chǎng)的單向耦合過(guò)程，并與Comsol提供的異構(gòu)網(wǎng)格數(shù)據(jù)映射方法進(jìn)行對(duì)比；最后利用Comsol軟件計(jì)算該變壓器的電磁熱場(chǎng)，對(duì)比分析使用均勻熱源及非均勻熱源時(shí)的變壓器溫度場(chǎng)情況。

1 坐標(biāo)插值數(shù)據(jù)映射法

以節(jié)點(diǎn)為數(shù)據(jù)傳輸載體的異構(gòu)網(wǎng)格數(shù)據(jù)映射方法的數(shù)據(jù)傳遞示意如圖1所示。圖1中，為獲得未知物理場(chǎng)中待求節(jié)點(diǎn)的物理量，需要先從和其耦合的已知物理場(chǎng)中找到該節(jié)點(diǎn)對(duì)應(yīng)的主單元，然后在已知物理場(chǎng)中進(jìn)行數(shù)據(jù)插值計(jì)算。

圖1 以節(jié)點(diǎn)為數(shù)據(jù)傳輸載體的數(shù)據(jù)映射示意圖Fig.1 Data mapping diagram of heterogeneous grid with node as data transmission medium

圖2為本文所提坐標(biāo)插值映射法的數(shù)據(jù)傳遞示意圖。首先在已知物理場(chǎng)中，利用已知節(jié)點(diǎn)數(shù)據(jù)對(duì)映射域中的物理量及其坐標(biāo)X、Y和Z進(jìn)行插值處理，將得到的插值結(jié)果和已知節(jié)點(diǎn)數(shù)據(jù)一起構(gòu)成一個(gè)新的數(shù)據(jù)集，并定義在和未知物理場(chǎng)坐標(biāo)系重合的一個(gè)坐標(biāo)空間中；然后，利用求解場(chǎng)域的空間坐標(biāo)不變性，根據(jù)未知物理場(chǎng)中的物理量節(jié)點(diǎn)坐標(biāo)，利用所定義的坐標(biāo)空間數(shù)據(jù)集進(jìn)行插值求得。

坐標(biāo)插值映射法中需要進(jìn)行2次插值運(yùn)算。第一次插值計(jì)算是在已知物理場(chǎng)中進(jìn)行，由已知節(jié)點(diǎn)數(shù)據(jù)求出計(jì)算域中坐標(biāo)點(diǎn)(X,Y,Z)的物理量，形成坐標(biāo)映射數(shù)據(jù)集，并將該數(shù)據(jù)集定義在和未知物理場(chǎng)坐標(biāo)系重合的坐標(biāo)空間中；第二次插值計(jì)算時(shí)依據(jù)待求點(diǎn)坐標(biāo)在未知物理場(chǎng)(即所定義的坐標(biāo)空間)中進(jìn)行。

圖2 坐標(biāo)插值映射法的數(shù)據(jù)映射示意圖Fig.2 Data mapping diagram of heterogeneous grid of the coordinate interpolation mapping method

以三維Delaunay三角剖分得到的四面體為例，第一次插值時(shí)，在其內(nèi)部可利用重心坐標(biāo)進(jìn)行線性插值[9]。假設(shè)m點(diǎn)是以A、B、C和D為頂點(diǎn)的四面體中的任意一點(diǎn)，如圖3所示。

圖3 四面體單元示意圖Fig.3 Diagram of a tetrahedron element

圖3中四面體4個(gè)頂點(diǎn)的直角坐標(biāo)為(Xi,Yi,Zi)(i=1, 2, 3, 4)，m點(diǎn)的直角坐標(biāo)為(X,Y,Z)，其重心坐標(biāo)(α1,α2,α3,α4)是以m為頂點(diǎn)、以四面體ABCD各底面為底所構(gòu)成的新四面體體積與原四面體ABCD的關(guān)系為：

(1)

式中：VABCD為原四面體體積；VmBCD、VmACD、VmABD、VmABC分別為m點(diǎn)與相應(yīng)的A、B、C、D組成的四面體體積。

點(diǎn)m的直角坐標(biāo)(X,Y,Z)與重心坐標(biāo)之間滿足：

(2)

在已知物理場(chǎng)中，m點(diǎn)處的物理量值為：

Sm=α1SA+α2SB+α3SC+α4SD.

(3)

式中Sm、SA、SB、SC和SD分別為m、A、B、C和D點(diǎn)處的物理量值。

將所有經(jīng)式(3)插值運(yùn)算求得的物理量及其坐標(biāo)，和已知節(jié)點(diǎn)物理量及其坐標(biāo)一起構(gòu)成一個(gè)數(shù)據(jù)集，并將其定義在某個(gè)坐標(biāo)系和未知物理場(chǎng)坐標(biāo)系重合的坐標(biāo)空間中。

第二次插值時(shí)，在所定義的坐標(biāo)空間中對(duì)數(shù)據(jù)集使用三維分片線性插值，以圖4所示六面體插值為例，插值函數(shù)可定義為[10]：

(4)

式中：x、y和z為插值點(diǎn)的坐標(biāo)；f(x,y,z)為插值點(diǎn)函數(shù)值(即插值點(diǎn)處的物理量值)；ai(i=0, 1, 2, …, 7)為插值系數(shù)，滿足矩陣方程：

MAT=QT.

(5)

式中：A=[a0a1…a7]；Q=[Q000Q100Q010Q110Q001Q101Q011Q111]，其中Q中各分項(xiàng)為六面體各頂點(diǎn)坐標(biāo)的物理量值；M為插值變換矩陣，且

其中，(xi,yj,zk)(i,j,k=0, 1)為六面體的頂點(diǎn)坐標(biāo)，如圖4所示。求解式(5)可得到系數(shù)ai(i=0, 1, 2, …, 7)的解析表達(dá)式。

圖4 六面體單元及其頂點(diǎn)示意圖Fig.4 Schematic diagram of a hexahedral element and its vertex coordinate

經(jīng)過(guò)2次插值運(yùn)算后，由式(4)即可得到未知物理場(chǎng)中插值點(diǎn)的物理量。第一次插值時(shí)插值點(diǎn)取的越多，可提供已知物理場(chǎng)量的分布信息就越多；在第二次插值時(shí)，插值點(diǎn)(待求點(diǎn))可取為插值域內(nèi)的任意點(diǎn)，包括映射前后網(wǎng)格節(jié)點(diǎn)和非節(jié)點(diǎn)、網(wǎng)格內(nèi)點(diǎn)和網(wǎng)格外點(diǎn)等[11]。網(wǎng)格數(shù)據(jù)映射誤差主要和2次插值時(shí)的算法及第一次插值所提供的數(shù)據(jù)集有關(guān)。

坐標(biāo)插值映射方法具體實(shí)施步驟包括：①由有限元法的后處理程序在整個(gè)映射域內(nèi)進(jìn)行第一次插值運(yùn)算和平滑處理；②提取經(jīng)第一步運(yùn)算后整個(gè)映射域的全部數(shù)據(jù)集(包括坐標(biāo)值及物理量值)，導(dǎo)入到未知物理場(chǎng)坐標(biāo)系所在的坐標(biāo)空間中；③根據(jù)未知物理場(chǎng)中節(jié)點(diǎn)坐標(biāo)，在所定義的坐標(biāo)空間中由有限元法前處理程序進(jìn)行第二次插值運(yùn)算。坐標(biāo)插值映射法避免了快速映射法搜索主單元的過(guò)程，也無(wú)須像徑向基函數(shù)法那樣求解大型方程組；同時(shí)，因第二次插值是在所定義的數(shù)據(jù)集中進(jìn)行，所以對(duì)網(wǎng)格外點(diǎn)無(wú)需特殊處理，相比混合法而言實(shí)現(xiàn)起來(lái)步驟要簡(jiǎn)單，而且對(duì)瞬態(tài)的電磁熱流耦合計(jì)算[12]，坐標(biāo)插值映射法可直接將三維線性插值擴(kuò)充到四維，直接導(dǎo)入瞬態(tài)熱源，使計(jì)算更簡(jiǎn)潔高效。

2 仿真算例

2.1 物理場(chǎng)計(jì)算方法

變壓器的金屬結(jié)構(gòu)件在交變電磁場(chǎng)作用下產(chǎn)生損耗，利用渦流場(chǎng)模型可求出變壓器內(nèi)部電磁場(chǎng)和渦流，得到各結(jié)構(gòu)件損耗分布。因中小型變壓器的磁通密度較小[13]，故忽略鐵心的磁滯損耗；同時(shí)油箱的渦流損耗占總損耗的比重較小[14]，計(jì)算中也忽略油箱的損耗影響。

2.1.1 變壓器渦流場(chǎng)及損耗計(jì)算

使用A-V-A法[15]計(jì)算變壓器的渦流場(chǎng)，在復(fù)頻域中，鐵心、夾件和油箱為渦流區(qū)，繞組及變壓器油為非渦流區(qū)。頻域似穩(wěn)態(tài)下各求解域的控制方程為：

(6)

式中：μ為磁導(dǎo)率，H/m；σ為電導(dǎo)率，S/m；ω為角頻率，rad/s；Je為繞組電流密度，A/m2；A為矢量磁位，Wb/m；V為標(biāo)量電位，V。

繞組采用均勻化線圈，同時(shí)考慮集膚效應(yīng)。電流密度Je由繞組實(shí)際相電流乘匝數(shù)后再除以繞組截面積求得。利用電流密度及電導(dǎo)率可直接求出繞組的歐姆損耗，而其渦流損耗受繞組磁通密度的影響，在只考慮基波時(shí)繞組的單元渦流損耗密度為[7]

Pwei=(σω2l2|Bi|2)/24.

(7)

式中:l為導(dǎo)線軸向尺寸，m；|Bi|為繞組單元磁通密度模，T。

鐵心采用均勻化材料模型，設(shè)y軸垂直于硅鋼片軋制方向，x和z軸沿硅鋼片軋制方向，則等效電導(dǎo)率矩陣及磁導(dǎo)率張量矩陣[16]為：

(8)

(9)

式中：F為鐵心疊積系數(shù)；a為疊片的寬度，mm；d為疊片厚度，mm；μfx、μfy和μfz分別為沿x、y和z方向的磁導(dǎo)率，H/m；μ0為真空磁導(dǎo)率，H/m。

按變壓器的實(shí)際尺寸[13]，取短軸方向的一半作為渦流場(chǎng)模型的計(jì)算域，采用有限元方法計(jì)算空載及額定負(fù)載下各結(jié)構(gòu)件的損耗密度分布。

2.1.2 變壓器熱流場(chǎng)計(jì)算

將變壓器油近似為不可壓縮層流。質(zhì)量方程為：

ρu=0.

(10)

式中：ρ為材料密度，kg/m3；u為流體速度場(chǎng)，m/s。動(dòng)量守恒方程為[17]

ρ(u)u=[-pI+τ(u+(u)T)]+

(ρ-ρref)g.

(11)

式中：p為油的壓力，Pa；I為單位張量；τ為動(dòng)力粘度，Pas；ρref為參考密度，kg/m3；g為重力加速度，m/s2。能量方程忽略粘性耗散等后，退化為熱對(duì)流擴(kuò)散方程

ρCpuT+(-λT)=Q.

(12)

式中：Cp為恒壓熱容，J/(kgK)；T為溫度場(chǎng)變量，K；為導(dǎo)熱系數(shù)，W/(mK)；Q為熱源，W/m3。

仿真計(jì)算時(shí)，油為流體對(duì)流換熱，結(jié)構(gòu)件為傳導(dǎo)傳熱。設(shè)置變壓器油箱外表面為傳熱學(xué)第三類(lèi)邊界條件，油箱入口為第一類(lèi)邊界條件，出口為自由流出，計(jì)及油箱外表面對(duì)環(huán)境輻射，中間剖面為對(duì)稱(chēng)邊界條件。忽略散熱器的影響，油箱入口油溫取303.15 K(30 ℃)，法向流入速度取0.05 m/s，出口為開(kāi)邊界，多物理場(chǎng)耦合采用流動(dòng)耦合及溫度耦合。變壓器主要材料參數(shù)見(jiàn)表1。

表1 變壓器材料參數(shù)Tab.1 Main material parameters of transformer

2.1.3 求解器配置

求解變壓器渦流場(chǎng)時(shí)，首先采用MUMPS求解器分析線圈幾何尺寸，再利用PARDISO直接求解器求解得到場(chǎng)量A、V及各繞組的電流。

對(duì)熱流場(chǎng)進(jìn)行分離式求解，非線性方法使用阻尼牛頓法，溫度T的阻尼因子為0.8，速度u和壓力p的阻尼因子為0.5。線性求解部分，溫度T使用PARDISO求解器，速度和壓力則使用代數(shù)多重網(wǎng)格。

2.2 網(wǎng)格剖分及數(shù)據(jù)映射

圖5為變壓器渦流場(chǎng)及熱流場(chǎng)仿真時(shí)繞組和鐵心的網(wǎng)格剖分情況，網(wǎng)格單元總數(shù)和網(wǎng)格映射特征見(jiàn)表2。

圖5 渦流場(chǎng)及熱流場(chǎng)計(jì)算時(shí)網(wǎng)格剖分情況Fig.5 Mesh generation of winding and iron core in eddy current field and thermal field simulation

表2 結(jié)構(gòu)件的網(wǎng)格單元總數(shù)統(tǒng)計(jì)Tab.2 Total number of grid units for the structural parts

在渦流場(chǎng)模型中，取拉格朗日點(diǎn)為計(jì)算點(diǎn)，在所有位置進(jìn)行平滑處理(閾值為1)。為了提高平滑效果，提取到與物理實(shí)際更相符的數(shù)據(jù)，損耗密度的分辨率取為常規(guī)，隨后導(dǎo)出繞組和鐵心的體積損耗密度值及其坐標(biāo)形成的映射數(shù)據(jù)集。在熱流場(chǎng)模型中，導(dǎo)入上述映射數(shù)據(jù)集，定義熱源類(lèi)型為廣義熱源，各節(jié)點(diǎn)的發(fā)熱功率可根據(jù)其坐標(biāo)在映射數(shù)據(jù)集中進(jìn)行三維線性插值計(jì)算得到。在Comsol軟件中通過(guò)定義三維插值函數(shù)可以將數(shù)據(jù)集導(dǎo)入熱流場(chǎng)模型中，隨后在熱源的設(shè)置中調(diào)用該插值函數(shù)，實(shí)現(xiàn)未知和已知物理場(chǎng)坐標(biāo)系的重合。對(duì)不同的結(jié)構(gòu)件均采用相同的方法。

為了驗(yàn)證本文方法的有效性，將本文的方法和Comsol軟件提供的異構(gòu)網(wǎng)格映射方法進(jìn)行了對(duì)比。該方法通過(guò)求解以下公式進(jìn)行平滑和映射：

P-P0=δh2(P).

(13)

式中：P為待求映射變量；P0為已知映射變量；δ為可調(diào)參數(shù)，取0.01；h為網(wǎng)格大小。在軟件中使用弱形式偏微分方程對(duì)式(13)進(jìn)行求解。

2.3 結(jié)果分析

圖6為額定負(fù)載時(shí)鐵心的損耗密度分布在映射前后的對(duì)比。由圖6可知：進(jìn)行數(shù)據(jù)映射后鐵心的損耗密度值及其分布與渦流場(chǎng)計(jì)算的結(jié)果(映射前)基本相同；相比Comsol提供的映射方法，采用坐標(biāo)插值映射法的計(jì)算結(jié)果與映射前更接近，在鐵心末級(jí)鐵拐角高損耗區(qū)映射精度更高。

圖6 額定負(fù)載時(shí)損耗密度分布映射前后對(duì)比Fig.6 Loss density distribution of low voltage winding and iron core before and after mapping at rated load

利用坐標(biāo)插值映射法，額定負(fù)載時(shí)高低壓繞組和鐵心損耗密度映射前后的數(shù)值范圍對(duì)比見(jiàn)表3。由表3可以看出：映射后物理量的數(shù)值范圍均有所縮小；與鐵心的數(shù)值相比，映射前后高、低壓繞組的數(shù)值分布范圍相對(duì)較小。

表3 額定負(fù)載下各結(jié)構(gòu)件損耗密度映射數(shù)值范圍Tab.3 Loss density range of the structural parts before and after mapping at rated load W/m3

以渦流場(chǎng)計(jì)算得到的各結(jié)構(gòu)件的總損耗值為基準(zhǔn)，計(jì)算其映射后的相對(duì)誤差

(14)

式中：Prl為映射后熱流場(chǎng)中熱源在各構(gòu)件域中的體積分，W；Pwl為映射前渦流場(chǎng)中體積損耗密度在各構(gòu)域中的體積分，W。計(jì)算結(jié)果見(jiàn)表4。對(duì)于坐標(biāo)插值映射法，繞組的映射是由稀疏網(wǎng)格向密集網(wǎng)格進(jìn)行，由于繞組損耗的數(shù)值范圍較小，因此映射誤差相對(duì)較??；鐵心損耗數(shù)值范圍大，但其映射是由密集網(wǎng)格向稀疏網(wǎng)格進(jìn)行，故映射誤差不大，且在額定負(fù)載時(shí)鐵心損耗的數(shù)值范圍大于空載時(shí)的數(shù)值范圍，故額定負(fù)載時(shí)鐵心的映射誤差比空載時(shí)大。對(duì)鐵心的總損耗相對(duì)誤差，利用坐標(biāo)插值映射法所得的結(jié)果比利用Comsol方法所得的結(jié)果略大。

表4 各結(jié)構(gòu)件損耗映射相對(duì)誤差Tab.4 Relative errors of loss mapping for the structural components %

本文所提方法最主要的耗時(shí)在于映射數(shù)據(jù)的提取和第二次插值，額定負(fù)載時(shí)這2步總耗時(shí)12.77 s，空載時(shí)耗時(shí)為10.32 s；而Comsol的方法最主要的耗時(shí)在于對(duì)式(13)的求解，在額定負(fù)載和空載時(shí)，這一步耗時(shí)均為71 s。因此本文應(yīng)用所提方法耗時(shí)相對(duì)較少。

圖7為額定負(fù)載下，按均勻熱源和非均勻熱源計(jì)算得到的變壓器油箱(中間剖面)和鐵心的穩(wěn)態(tài)溫度分布。由圖7可以看出：熱源非均勻時(shí)主要對(duì)鐵心的溫度場(chǎng)分布產(chǎn)生影響；使用非均勻熱源時(shí)鐵心上端外頂角處和下鐵軛的溫度較使用均勻熱源時(shí)低，B相心柱中間偏上部位溫度有所升高。

3 結(jié)論

本文提出了一種用于異構(gòu)網(wǎng)格數(shù)據(jù)映射的坐標(biāo)插值映射新方法，并得出以下結(jié)論：

a)坐標(biāo)插值映射法有較高的映射精度，計(jì)算步驟簡(jiǎn)單且無(wú)需考慮網(wǎng)格外點(diǎn)問(wèn)題。

b)在本文所提的計(jì)算條件下，額定負(fù)載時(shí)高低壓繞組和鐵心的損耗映射相對(duì)誤差分別為0.029 4 %、0.093 8 %、1.050 1 %；空載時(shí)，鐵心損耗的映射相對(duì)誤差為0.714 0 %。

c)與Comsol的映射方法相比，坐標(biāo)插值映射法在物理量分布上更精確、耗時(shí)短；但對(duì)于映射方向由稀疏網(wǎng)格向密集網(wǎng)格、數(shù)值變化范圍較大、網(wǎng)格單元差異較大的情況時(shí)相對(duì)誤差較大。

圖7 額定負(fù)載下油箱和鐵心在均勻熱源與非均勻熱源時(shí)的溫度分布對(duì)比Fig.7 Temperature distribution of tank and iron core with uniform and non-uniform heat sources at rated load

d)對(duì)于中小型變壓器，額定負(fù)載時(shí)非均勻熱源對(duì)鐵心的溫度場(chǎng)分布影響主要體現(xiàn)在鐵心柱B相和鐵軛的溫度分布上。