張建瓴，周壯廣，可欣榮

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

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

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

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

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

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

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

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

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

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

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

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

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

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

(1)

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

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

(2)

在已知物理場中，m點處的物理量值為：

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

(3)

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

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

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

(4)

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

MAT=QT.

(5)

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

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

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

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

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

2 仿真算例

2.1 物理場計算方法

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

2.1.1 變壓器渦流場及損耗計算

使用A-V-A法[15]計算變壓器的渦流場，在復(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。

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

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

2.1.2 變壓器熱流場計算

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

ρu=0.

(10)

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

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

(ρ-ρref)g.

(11)

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

ρCpuT+(-λT)=Q.

(12)

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

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

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

2.1.3 求解器配置

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

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

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

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

圖5 渦流場及熱流場計算時網(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)計Tab.2 Total number of grid units for the structural parts

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

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

P-P0=δh2(P).

(13)

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

2.3 結(jié)果分析

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

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

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

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

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

(14)

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

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

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

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

3 結(jié)論

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

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

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

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

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

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