摘要 油浸式电力变压器绕组温升研究对于保障电力系统的安全运行具有重要意义。针对经典Galerkin有限元法在求解对流扩散方程时存在的非物理振荡，以及基于结构化交错网格的有限体积法在处理流固耦合分界面时计算复杂度高的问题，该文提出了一种将二阶迎风有限体积法（SOUFVM）与流线迎风有限元法（SUFEM）结合的流热耦合计算方法。为了避免交错网格中的额外插值操作，采用基于结构化同位网格的二阶迎风有限体积法求解流场控制方程，以获得油流速度分布。同时，为了有效避免传统有限体积法中因确保热流密度在流固耦合界面的连续性而需单独处理分界面的问题，采用流线迎风有限元法求解油流与绕组耦合传热的统一温度场方程，以获得温升分布。以一台额定容量为321.1 MV·A，额定电压为530/
kV的油浸式电力变压器为研究对象，建立二维单分区不分匝绕组仿真模型，并采用所提方法对其油流与温升分布进行计算。将计算结果与一阶迎风有限体积法和伽辽金有限元法相结合（FOUFVM-GFEM）的方法及Fluent软件计算结果进行对比验证。结果表明，该方法在数值精度上优于FOUFVM-GFEM方法；与Fluent结果相比，油流速度峰值的绝对误差为0.000 8 m/s，相对误差为1.38%；热点温度的绝对误差为2.06 K，相对误差为0.6%。总体而言，该方法的相对误差均控制在2%以内，充分验证了其在油流与温升计算中的有效性，具有较高的工程应用价值。
关键词：油浸式电力变压器 绕组温升 有限体积法 有限元法 流场 温度场
变压器绕组温升是影响其运行性能的重要指 标[1-2]。研究表明，绕组热点温升过高是导致变压器故障的主要原因之一[3-4]。因此，准确掌握绕组温升分布及热点温度对于保障设备安全运行至关重要。然而，变压器内部复杂的多物理场耦合使得传统的单一物理场分析方法难以全面反映其流热特性。为此，本文提出一种高效的多物理场耦合求解方法，为国产高压电力装备的多物理场仿真软件提供关键技术支撑。
目前，针对油浸式电力变压器绕组温升的研究方法主要分为直接测量法和间接计算法两类[5]。直接测量法通过在变压器内部埋设传感器，直接获取热点温度。该方法具有较高的测量精度，能够实现实时监控，并且由于传感器埋设在变压器内部，外界环境因素对测量的影响较小，因此测量稳定性较高[6-7]。然而，直接测量法难以提供变压器内部各部分状态的全面信息[8]。与此相比，间接计算方法通过运用计算公式或仿真模型来估算变压器的热点温度，主要包括经验公式法、热路模型法和数值模拟法等。间接计算方法能够保持变压器结构的完整性，但其计算精度依赖输入数据的精度和计算模型的准确性。如果输入数据或计算模型存在误差，计算结果也会受到影响[9]。在这些方法中，IEEE STD C57.91—1995推荐的变压器绕组热点温度经验模型被广泛应用[10]。该模型通过计算油浸式电力变压器的环境温度、顶部油温与环境温度之间的温差，以及顶部油温与绕组热点温度之间的温差，来推算绕组的热点温度。尽管该方法计算简便，但其结果往往存在较大的误差。此外，文献[11]提出了基于热电类比的变压器顶层油温模型，该热路模型无法准确反映变压器内部的温度分布，因此在实际应用中存在较大的局限性。由于经验公式法和热路模型法均未考虑变压器内部的温度场分布及绕组结构与散热介质的相互作用，因此，近年来学者们开始转向基于多物理场耦合的数值模拟方法，这些方法能够更全面地模拟变压器内部的流动与传热特性。
随着计算机技术的快速发展，数值模拟已成为分析变压器绕组温升分布的重要手段。常用的数值方法包括有限元法[12-13]、有限体积法[14-15]以及有限差分法[16]等。有限体积法因其严格满足局部守恒性，在流动与传热问题中得到了广泛应用[17]，其常用的中心差分格式在对流主导问题中稳定性较差，而一阶迎风格式则存在数值耗散，限制了模拟精度。为了提高计算稳定性与精度，二阶迎风有限体积法（Second-Order Upwind Finite Volume Method, SOUFVM）被提出，通过采用高阶迎风格式，在保证守恒性的同时有效抑制数值振荡，克服了传统方法在精度上的局限。此外，为了避免因交错网格导致的额外插值操作，采用结构化同位网格布置方式，在同一节点上同时定义速度、压力与温度变量，从而简化了数值实现流程，降低了计算复杂度[18-20]。另外，有限元法因其良好的几何适应性和高阶逼近能力，在处理复杂边界和多层介质问题中具有较高的精度[21]。但经典Galerkin有限元法在求解不可压缩流体问题时，可能引发非物理振荡，从而影响计算结果的准确性[22]。为此，流线迎风有限元法（Streamline Upwind Finite Element Method, SUFEM）通过在流线方向引入人工黏性，增强了数值稳定性，尤其在处理对流-扩散型主导的温度场问题中表现出良好性能。与传统有限体积法需在流固界面对热流密度进行插值处理不同，SUFEM能够在离散过程中自然满足温度与热流密度在界面上的连续性，从而显著简化边界处理过程[23]。
综合利用SOUFVM与SUFEM的各自优势，本文提出了一种流-热耦合求解策略：采用基于结构化同位网格的SOUFVM求解流场控制方程，以获得绕组内部油流的速度分布；同时，采用SUFEM求解油流与绕组耦合传热的统一温度场控制方程，以获得温升分布。最后，通过顺序迭代策略交替求解流场与温度场方程，从而获得整体场域的最终稳态温度分布。
以一台额定容量为321.1 MV·A、额定电压为530/ kV的油浸式电力变压器为研究对象，建立二维单分区不分匝绕组仿真模型，采用所提方法对其油流与温升分布进行计算，将计算结果与一阶迎风有限体积法和伽辽金有限元法相结合（First- Order Upwind Finite Volume Method-Galerkin Finite Element Method, FOUFVM-GFEM）的方法以及基于有限体积法的Fluent软件进行对比，以验证所提方法在油流与温升分析中的有效性。
变压器油通常被近似为不可压缩牛顿流体。根据流体动力学，油浸式电力变压器内变压器油的流动遵循质量守恒定律和动量守恒定律，具体的方程表达式[24]为
（1）
（2）
式中，U为速度矢量（m/s）；
为流体密度（kg/m3）；
为动力黏度（N·s/m2）；p为压强（Pa）；f为外力密度矢量（N/m3）。在式（2）中，左端第一项为对流项，第二项为扩散项，第三项为压力梯度项。
采用有限体积法对式（2）进行离散，得到控制方程为
（3）
式中，V为控制体积。
流体流动的边界条件通常包含壁面边界、入口边界和出口边界，其处理方式如下[25]。
1）壁面边界条件
当壁面为非渗透性时，壁面边界条件为速度分量u=0和v=0，即无滑移边界条件，可表示为
（4）
式中，n为法线方向单位矢量。
2）入口边界条件
当入口边界条件为速度边界条件时，需要指定入口处的法向和切向速度，可表示为
（5）
式中，u0和v0分别为已知的法向和切向速度。
3）出口边界条件
在计算过程中，通常根据具体的计算模型假设出口流体的性质，并据此设定出口边界条件。本文采用一种常见的压力出口边界条件，其具体形式为
（6）
式中，p0为已知的出口压力值。
基于上述离散方法和边界条件，式（3）离散后形成线性方程组为
（7）
式中，K为刚度矩阵（包含对流与扩散项）；z为待求解的速度向量（包含所有计算单元在水平和竖直方向的速度分量组成的列向量）；F为源项（包含压力梯度项、边界条件项及修正项）。
对流项的一个重要特征是具有迁移性。因此，数值计算时需要采用考虑这些物理特性的离散格式，以确保计算的稳定性和精度。迎风格式便是针对该特性提出的一种方法，其核心思想是从上游方向获取信息以离散对流项[26]。
根据迎风格式的思想，其一阶、二阶离散格式[27]为
（8）
式中，
和
分别为采用一阶和二阶迎风格式计算得到的界面速度；Uf、UC分别为面f中心的速度和单元C质心的速度；dCf为单元C质心至面f中心的距离。
二阶迎风格式相比一阶迎风格式在数值精度方面具有优势，其理论截断误差可由O(Dx)提高至O((Dx)2)，同时具备更低的数值耗散特性。然而，已有研究表明，二阶格式在处理强对流或复杂几何条件下的流动问题时，容易引发非物理振荡，导致数值不稳定[28]。为了解决上述问题，本文在一阶迎风格式的基础上，提出一种融合延迟修正技术的混合格式方法。该方法以一阶格式为基准结构，引入二阶迎风项作为修正补偿，在保持算法稳定性的同时有效提升数值精度，具体构建方式为
（9）
其中
mf=rUfSf
式中，mf为对流质量流系数；Sf为控制体面f的法向面积。为了兼顾数值稳定性与计算精度，采用隐式处理方式构建一阶迎风格式的对流通量，并在此基础上，通过显式引入二阶迎风修正项，形成稳定且高精度的混合离散格式。
然而，在实际计算中，采用二阶迎风格式虽然能提高精度，但也会增加计算成本。特别是在方程组装过程中，传统方法通常需要逐一计算每个单元的通量，这会导致相邻单元界面处的冗余计算，从而降低效率并增加计算开销。为了优化这一过程，本文在内部面计算对流通量时，以一阶迎风格式为基础，仅对每个界面执行一次质量流计算，并结合界面法向方向自动区分上下游单元，从而避免重复计算。同时，为了兼顾精度和稳定性，引入延迟修正思想，将二阶迎风格式作为显式修正项叠加于一阶格式之上，形成结构清晰、计算高效的混合离散格式。一阶迎风格式示意图如图1所示。在结构化网格中，由于拓扑关系规则，相邻单元共享统一的界面信息，因此单元C和单元D在分界面的通量分别为rSfUfnCUC和rSfUfnDUC。其中，UC为单元C质心的速度，nC、nD分别为单元C和单元D分界面上的法线方向单位矢量。
图1 一阶迎风格式示意图
Fig.1 Schematic of the first-order upwind scheme
在对流项的计算中，边界通量的处理至关重要。本文采用已知对流质量流系数和边界值计算边界面上的对流通量，并将计算结果添加到源项。具体处理方法为
（10）
式中，Fconv为对流通量；Ubdy为边界法向速度值。
此外，在数值计算中，除了对流项的处理，扩散项的离散格式同样影响计算精度与稳定性。因此，合理的扩散项离散方法也是实现高效计算的关键。对于扩散项，若各个面采用单点积分，则在结构化网格上[27]可表示为
（11）
式中，
为单元C质心与单元D质心间的距离；UD为单元D质心的速度。
在传统有限体积法中，通量计算与数据结构通常基于单元构建。然而，在处理内部面通量及复杂几何边界时，此类结构常导致重复计算和刚度矩阵组装效率低下。针对这一局限，本文提出了一种基于面-单元的数据结构优化方法。该方法以面为循环控制单元，避免了对相同面的重复访问，有效降低了通量计算的冗余。计算结果可直接分配至相邻单元，简化刚度矩阵的组装过程。
对于第一个单元，其计算表示为
（12）
而对于相邻单元，其计算则可表示为
（13）
在结构化网格的边界面上，由于法向方向与梯度计算方向一致，可直接沿梯度方向计算扩散通量，并将计算结果添加到源项。具体处理方法为
（14）
式中，Fdiff为扩散通量；mf为面上的动力黏度；ddf为相邻单元间的距离。
除了对流项和扩散项外，压力梯度项的离散化同样是数值计算中的关键步骤。在其离散化过程中，单元质心与单元面梯度的计算至关重要。对于结构化网格，二阶迎风格式要求精确计算这些梯度。然而，格林-高斯方法在灵活性和精度方面存在一定的局限性[27]。为此，本文采用最小二乘法计算单元质心的梯度，并通过插值方法将其传递至单元面，以获得更加准确的面梯度。这一改进有效提高了压力梯度项离散化的精度和稳定性。
此外，为了进一步优化压力梯度项的离散化，利用前一迭代步骤的计算结果来求解各单元质心的梯度，具体计算方法[27]为
（15）
（16）
式中，pC和pD分别为单元C和单元D质心的值；
为单元C质心处的梯度；NB(C)为单元C的相邻单元数量；wi为权重因子；pi为第i个相邻单元质心处的值；Dpi为第i个相邻单元与当前单元之间的压力差；ri和rC、rD分别为第i个相邻单元和单元C、D质心的空间坐标位置。
为了精确计算质心梯度，本文采用一种优化程序，通过最小化目标函数GC来实现梯度的求解，具体实现步骤[27]为
（17）
式（17）可以表示为矩阵形式，有
（18）
式中，nb为相邻单元的数量；∂p/∂x、∂p/∂y分别为质心处沿x和y方向的压力梯度分量；Dxi和Dyi分别为第i个相邻单元和当前单元之间在x和y方向的坐标差。该方法不仅提高了梯度计算的精度，还增强了其在结构化网格中的适应性。
在此基础上，采用跨界面相邻单元质心梯度的加权平均值来计算单元面的梯度，从而有效捕捉界面变化。通过合理估计相邻单元间的梯度，该方法确保在结构化网格中实现更高精度的计算。具体实现方法[27]为
（19）
式中，gC为界面f相对于质心C位置相关的几何加权因子；
为单元D质心处的梯度；
为质心D与界面f之间的距离。
在前述对流项、扩散项及压力梯度项离散化方法的基础上，进一步展开式（7），并将其转化为有限体积法的标准离散形式。具体为
（20）
式中，KC为单元C的刚度矩阵系数，表示单元C自身对速度的影响，其主要由对流项和扩散项的离散化贡献；KCD为相邻单元D对单元C速度影响的耦合系数矩阵，其来源于对流项与扩散项的离散化；L为压力梯度项，表示压力对动量方程的贡献；B为边界条件项，考虑边界速度或外部力的影响（本文不考虑外力）；R为修正项；D～NB(C)表示单元C的邻接单元集合。
式（20）可进一步拆分为x和y方向上的分量，其具体表达式为
（21）
式中，N为计算单元的总数；uC、uD和vC、vD分别为单元C和单元D在x和y方向上的速度分量；
、
和
、
分别为在x和y方向上，采用二阶与一阶迎风格式计算得到的界面速度分量；gN为扩散通量系数；mfC为单元C的对流质量流系数；mfD为单元D的对流质量流系数；Bxn、Byn分别为x、y方向的边界条件项。
式（21）中各项代表不同的物理作用，其中等式左边第一项为主对角项，表示单元C的刚度矩阵系数KC，用于控制单元C处的速度UC。gN表示扩散项对单元C速度的影响，其计算方法为
（22）
此外，mfC用于描述对流项在单元C处的影响，具体计算方法为
（23）
式中，
为单元面上的流速投影。
在式（21）中，等式左边第二项为邻接单元项，表示单元C受到邻接单元D速度的影响（KCDUD）；gN为扩散项对单元D的贡献，计算方式与前述相同。mfD具体计算方法为
（24）
此外，式（21）中的等式右边第一项为压力梯度项，表示压力梯度对速度的驱动，具体表述为
（25）
最后，式（21）中等式右边第二项为边界条件项，其具体形式取决于边界类型。等式右边第三项和第四项为二阶迎风格式的延迟修正项，分别在上游单元C和下游单元D进行修正。
综上所述，将式（21）所有单元的离散方程写成矩阵形式，有
（26）
式中，Kx、Ky分别为x、y方向的刚度矩阵（包含对流项和扩散项）；zx、zy分别为x、y方向的速度分量组成的列向量；Lx、Ly分别为x、y方向的压力梯度项；Bx、By分别为x、y方向的边界条件项；Rx、Ry分别为x、y方向的二阶迎风格式延迟修正项。
接下来，将式（26）展开，具体表达式为
（27）
式中，
为所有相邻单元的扩散贡献求和；
为从相邻单元流入单元C的贡献。
在此基础上，最终zx、zy通过并行稀疏直接求解器（Parallel sparse Direct Solver, PARDISO）求解线性方程组，从而获得速度场的解为
（28）
基于上述理论，设计了典型二维方腔算例用于对比实验，评估同位网格与交错网格有限体积法在计算时间与总迭代步数等关键性能指标上的差异。为了确保分析的全面性，统一设置了网格数量、初始条件及收敛准则等参数。结果表明，本文提出的同位网格方法因避免交错网格中界面变量的插值操作，形成了更为简洁的数值实现结构，在保证解精度的同时显著降低了计算成本，详见图2和表1。
图2 模型中心线速度分量分布
Fig.2 Distribution of velocity components along the centerline of the model
表1 计算结果对比
Tab.1 Comparison of computational results
方法网格数量迭代次数总耗时/s 交错网格37×3750829.62 同位网格37×3732415.43
由图2可见，基于同位网格方法计算所得的竖直中心线上的水平速度分布与U. Ghia采用的涡量-流函数法结果更为吻合，进一步验证了本文方法的数值准确性[29]。
由表1可见，同位网格相比交错网格的迭代次数减少了184次，计算时间缩短了14.19 s。可见，该方法不仅优化了数值实现结构，在计算效率与收敛性方面也表现出更优的性能，验证了其在标准算例中的有效性。因此，本文进一步将该方法拓展应用于油浸式电力变压器的油流与温升分析，以评估其在复杂工程问题中的有效性与可行性。
油浸式变压器内流体和固体的传热过程遵循能量守恒定律，稳态温度场的控制方程[30]可表示为
（29）
式中，cp为比定压热容[J/(kg·K)]；k为导热系数[W/(m·K)]；HT为单位体积热源的产热率（W/m3）；T为温度（K）。
当流体速度为零时，式（29）中的能量守恒方程将退化为固体的热传导方程。因此，式（29）能够统一描述固体与流体的传热过程。在求解变压器内部的传热问题时，首先通过式（1）和式（2）计算油道中的油流速度，同时在固体区域设置速度为零，利用式（29）求解流体和固体区域的温度分布。
针对流固耦合传热问题，传统方法通常采用分区求解策略，即分别构建流体与固体区域的传热模型，并通过迭代调整耦合边界条件以实现区域间的信息交换。相比之下，本文采用整场求解方法，将流体与固体区域统一纳入一个传热控制方程中，从而避免分区方法中对耦合边界的重复处理。场域的边界条件及流固耦合分界面条件[31]为
（30）
式中，G1、G2、G3分别为第一、二、三类边界条件；∂T/∂n为温度T沿界面法线方向n的梯度；Ta为流体远场温度；h为表面传热系数；qs为表面热流密度；S1、S2分别为流固分界面上固体侧与流体侧的边界；S12、S21为流固分界面；k1与k2分别为固体材料与流体材料的导热系数。式（30）中第三、四式为分界面处的耦合边界条件。
针对温度场控制方程，采用加权余量法构造有限元方程，并对求解区域进行离散，具体形式为
（31）
式中，W为求解率；Wi为权函数；n为求解区域内有限元网格节点数。
在经典Galerkin有限元法中，对流项通常被处理为对称形式，尽管该形式简化了节点间的耦合关系，但容易引发数值振荡甚至导致求解失败。针对这一问题，本文在此基础上引入流线迎风有限元法，以改进式（31）所描述的流体传热问题的对流项处理。该方法在标准Galerkin权函数的基础上引入一个与流线方向相关的修正项，从而增强了数值稳定性，其具体形式[32]为
（32）
其中
（33）
式中，Pe为佩克莱特数；Ni为网格节点温度的插值基函数；
为修正系数；l为单元特征长度。
本文采用四节点四边形单元的形状函数计算T的近似解，在单元内的表示形式[33]为
（34）
式中，
为单元e中第j个节点的温度。
根据格林公式，将式（31）中的二阶偏微分项转化为一阶偏微分项，其具体推导为
（35）
将式（35）代入式（31），可得
（36）
式中，G为整体场域的边界。根据式（30）第四式可知，式（36）中第三项积分为零，这表明有限元法能够自动确保两种介质分界面上热流密度的连续性。在离散整体场域时，分界面处仅有一个节点，且该节点上的温度值唯一。因此，耦合分界面处的温度连续性条件式（30）第三式得以自动满足。
将式（30）第二式中的第二、三类边界条件代入式（36）中，可得
（37）
将式（34）代入式（37），得到温度场控制方程的流线迎风有限元离散形式为
（38）
式中，ne为整个场域内的单元总数；be为满足第二类和第三类边界条件的单元总数。
由于式（38）已经包含第二、三类边界条件，因此无需额外处理。对于第一类边界条件，温度为常数，式（36）中第一类边界条件对应的边界积分设为零。在构造有限元刚度矩阵后，通过强加边界条件的方法处理第一类边界条件，从而求取场域的温度分布。
稳态流场与温度场采用间接耦合计算方法，即依次顺序迭代求解流场和温度场的控制方程，直至收敛得到稳定的数值解。流场与温度场耦合计算流程如图3所示。
初始化速度、压力、界面流速和对流质量流系数，采用压力耦合方程的半隐式方法（Semi-Implicit Method for Pressure Linked Equations, SIMPLE）求解动量方程，从而计算出预测速度场u*、v*（u*、v*为未修正的值），有
（39）
采用Rhie-Chow插值法计算单元表面的速度，具体表达式[27]为
图3 流场与温度场耦合计算流程
Fig.3 The flow and temperature field coupling calculation flowchart
（40）
式中，
为初步计算的界面速度；
为用于调整界面速度计算的系数；
为相邻单元界面f上的压力梯度；
为相邻单元界面f上压力梯度的平均值。
将式（40）转化为矢量形式，以便适用于多维压力修正方程，具体形式为
（41）
其中
（42）
式中，pC、pD分别为界面两侧单元C和D的压力值；eCD为从单元C指向单元D的单位向量。
相比传统的线性插值方法，Rhie-Chow插值法在网格类型适应性方面具有显著优势，特别是在捕捉复杂流场中界面流速变化方面。该方法不仅能够有效避免棋盘效应，还能确保计算结果的物理合理性，并在复杂边界和流动条件下提高界面流速的准确性与稳定性。
动量方程求解得到的速度分量u*、v*尚未满足连续性方程。为了保证质量守恒，本文对速度场和压力场进行修正。设
和
分别为速度场和压力场的修正量，其关系为
（43）
式中，n表示第n次迭代。
根据式（41），将式（43）代入连续性方程，可得
（44）
式中，
为界面f上的速度；S为单元面的法向面积；
为通过动量方程求得的未修正的质量流系数。
为了提高计算效率，SIMPLE算法忽略式（44）右侧括号中的项，从而将速度修正与压力修正直接关联。由于该方程为修正方程，修正项在迭代过程中将逐渐趋于零，因此省略这些项不会影响最终解的精度。其具体形式为
（45）
修正后的方程简化为仅包含扩散项和源项的压力泊松方程，从而在保持计算精度的同时提高计算效率。
在SIMPLE算法中，忽略某些方程项可能导致较大的压力修正，进而减缓收敛速度或导致发散。为了提高稳健性和收敛性，本文引入了显式欠松弛策略，该策略仅对压力修正进行欠松弛处理，而速度更新则无需处理，因为通过调整压力修正，能够有效保证质量守恒。具体的修正方程为
（46）
式中，VC为控制体积的大小；aC为单元C的离散化方程系数；
为欠松弛因子[27]。
在扩散问题中，扩散系数通常依赖未知因变量
，并可能产生较大的交叉扩散项[27]。为了解决这一问题，本文引入延迟修正方法。然而，在迭代过程中，的显著变化可能引起源项和系数的较大波动，导致迭代发散。这种发散源于系数及交叉扩散项的非线性，进而使源项受到尚未收敛的解场影 响[27]。为了促进收敛并稳定迭代过程，采用欠松弛技术减缓的变化。
标准的离散化方程通常采用形式[27]为
（47）
式中，aC、aD分别为单元C与邻接单元D相关的系数；
、
分别为单元C与邻接单元D处的未知量；bC为与单元C相关的源项。
通过使用前一次迭代的值进行修正并利用松弛因子
进行调整，式（47）可表示为
（48）
式中，
为前一次迭代的值。
通过变换，式（48）可表示为
（49）
在流场-温度场耦合计算中，可以结合有限体积法和有限元法，充分发挥两者在流动和传热模拟中的优势。在流场计算中，有限体积法通过在单元质心定义物理量，确保物理量的局部守恒性。为了提高节点物理量的计算精度，采用反距离加权法，利用相邻单元质心物理量的加权平均值来估算节点变量[34]。为了减少误差累积并提高数值稳定性，对反距离加权法进行了优化，通过归一化处理权重，使其总和为1。这一优化不仅提高了计算精度，还显著增强了在复杂网格中的稳定性。在温度场计算中，有限元法通常在单元节点定义物理量，并通过插值函数近似单元内的物理量分布。在实际应用中，通过加权平均四个节点的温度值来估算单元质心的温度，以提高计算精度。结合有限体积法的守恒特性与有限元法的高精度插值方法，能够更准确地模拟实际工程中的流体流动和热传递过程。
文献[29]研究表明，绕组的不分匝模型与分匝模型在流场和温度场的分布上基本一致。考虑到绕组内部的挡板结构，这些挡板将绕组划分为若干个结构相似的分区，其中油流自入口进入每个分区，沿S形路径与绕组充分接触后通过出口流出。本文以一台额定容量为321.1 MV·A、额定电压为530/
kV的油浸式电力变压器为研究对象，采用二维单分区不分匝绕组传热模型，并将计算结果与FOUFVM- GFEM方法及Fluent进行对比分析，以验证SOUFVM- SUFEM方法在绕组稳态温升计算中的有效性与准确性。
程序开发在VS Code环境中进行，计算使用的设备配置包括12th Gen Intel(R) Core(TM) i9-12900 2.40 GHz处理器和128 GB内存。油浸式电力变压器绕组结构与温升计算模型如图4所示。
图4 绕组结构与温升计算模型
Fig.4 Schematic diagram of the winding structure and the temperature rise calculation model
图4所示二维绕组温升计算模型的尺寸为100 mm×125 mm，包含8个尺寸为88 mm×10 mm的线饼。线饼之间的水平油道宽度为5 mm，油流的竖直进出口油道宽度为6 mm。该模型采用结构化网格划分，共有12 500个单元和12 726个节点。网格划分与油流流线如图5所示。
图5 网格划分与油流流线
Fig.5 Grid generation and oil flow streamline
在油浸式电力变压器绕组温升模拟中，需要根据绕组金属导体和流体流动特性设置相应参数。模型中的绕组和流动区域分别采用铜和变压器油的物理性质进行定义，各部分材料物性参数详见表2。
在求解绕组温升分布时，首先将计算模型离散为不重叠的结构化网格。并设定相应的边界条件：
（1）油流入口：u=0 m/s，v=0.05 m/s。
（2）初始油温：330 K。
（3）油流出口：u=0 m/s，p=0 Pa。
（4）由于绝缘筒具有较小的导热系数，因此视为绝热表面。绕组与油流分界面及绝缘筒的表面均设置为无滑移壁面边界条件，即u=0 m/s，v=0 m/s。
表2 材料物性参数[35]
Tab.2 Material properties parameters[35]
材 料参数类型数 值 变压器油密度/(kg/m3)1 098.72-0.712T 动力黏度/(N·s/m2)0.084 6-4×10-4T+5×10-7T2 比定压热容/[J/(kg·K)]807.2+3.58T 导热系数/[W/(m·K)]0.150 9-7.1×10-5T 铜密度/(kg/m3)8 900 比定压热容/[J/(kg/K)]381 导热系数/[W/(m·K)]387.6
（5）绕组温升的主要热源为电流通过绕组时产生的电阻损耗。考虑到温度对绕组电阻的影响，电阻损耗应视为温度的函数[36-37]，有
（50）
式中，T0为初始温度，T0=273 K；H0为在该温度下的绕组损耗，在本算例中参考了文献[38]取值为227 kW/m3；
为导体温度系数，对于铜导体，= 0.003 93 K-1。由于绕组体积对温度的依赖较小，因此式（50）可视为绕组损耗密度与温度的函数。
图6展示了采用本文所提SOUFVM-SUFEM方法与Fluent软件计算的油流速度分布云图。从图中可以看出，本文方法与Fluent在稳态条件下的计算结果基本一致，油流速度峰值的绝对误差为0.000 8 m/s，相对误差为1.38%，这些结果验证了所提计算方法在油流速度计算中的有效性。
图6 模型油流速度分布云图
Fig.6 The flow velocity distribution of the oil model
为了进一步验证计算方法的有效性，提取了竖直中心线上的水平速度值进行分析。由图7可见，FOUFVM-GFEM方法与Fluent的计算结果相比，油流速度的最大绝对误差为0.001 083 m/s，相对误差为7.9%。而本文方法得到的竖直中心线水平速度与Fluent的计算结果相比，油流速度的最大绝对误差为0.000 88 m/s，相对误差为6.3%。这一结果进一步验证了所提计算方法的有效性，同时也表明，本文方法相较于FOUFVM-GFEM方法具有更高的精度。
图7 模型竖直中心线油流速度分量分布
Fig.7 The velocity component distribution map of oil flow along the vertical centerline of the model
图8展示了采用本文所提SOUFVM-SUFEM方法与Fluent计算得到的模型温度分布云图。从图中可以看出，本文方法计算得到的稳态热点温度为341.59 K，热点位置坐标为 (82, 71)，位于第4号线饼处；Fluent计算得到的热点温度为343.65 K，热点位置坐标为 (73, 72)，同样位于第4号线饼处。两者之间的绝对误差为2.06 K，相对误差为0.6%，这些结果验证了所提方法在温度分布计算中的有效性。
图8 模型温度分布云图
Fig.8 The temperature distribution contour map of the model
为了进一步验证计算方法的有效性，本文提取了竖直中心线上的线饼温度进行分析。模型竖直中心线线饼温度分布如图9所示。由图9可见，本文方法得到的竖直中心线线饼温度与Fluent计算结果存在一定误差，最大绝对误差为3.42 K，相对误差为1.01%。这一差异主要源于两种方法在离散原理上的本质区别：SUFEM通过流线迎风加权的有限元形式能够自然保证流固耦合界面热流密度的连续性，而Fluent基于有限体积法，其对流固耦合界面的处理依赖显式插值，这可能导致界面附近温度场出现系统性偏差。尽管数值结果存在差异，但相对误差仍控制在允许范围内，并且从线饼温度的总体变化趋势来看，二者基本一致。与FOUFVM-GFEM相比，本文方法在计算精度上有所提高，进一步验证了所提计算方法的有效性。
图9 模型竖直中心线线饼温度分布
Fig.9 The temperature distribution map of the disc along the vertical centerline of the model
图10展示了模型右侧竖直油道中心线的温度分布。从图中可以看出，FOUFVM-GFEM与Fluent的计算结果之间存在较大的数值误差，而本文方法与Fluent计算结果的误差较小，这一结果验证了所提方法在有效性与准确性方面的优势。
图10 模型右侧竖直油道中心线温度分布
Fig.10 The temperature distribution along the vertical oil passage centerline on the right side of the model
图11展示了模型顶部水平油道中心线的温度分布。可以看出，FOUFVM-GFEM方法在模型左侧竖直油道与水平油道交界处出现轻微的数值振荡，并且与Fluent的计算结果相比存在较大误差。而本文方法与Fluent计算结果的误差较小，有效避免了非物理振荡，并提高了计算精度和稳定性。
图11 模型顶部水平油道中心线温度分布
Fig.11 Temperature distribution along the horizontal oil passage centerline at the top of the model
本文提出了一种结合二阶迎风有限体积法（SOUFVM）与流线迎风有限元法（SUFEM）的流-热耦合计算方法，以提高油浸式电力变压器绕组温升计算的数值精度与稳定性，并降低计算复杂度。
1）Fluent作为基于有限体积法（FVM）开发的通用流体力学与传热模拟软件，采用交错网格与SIMPLE算法实现压力-速度耦合，虽能有效抑制非物理振荡并保证守恒性，但在流固耦合中，由于变量位置不一致，需要进行显式插值操作，以保证界面热流密度与温度的连续性，从而增加了数值实现的复杂度。相比之下，本文提出的方法采用同位网格布置，在同一网格节点上离散速度、压力与温度等变量，避免了交错网格方法中因变量位置不一致而需进行插值操作，从而简化了数据结构并提高了计算效率。通过SUFEM对温度场的离散，该方法无需显式插值即可自然实现温度与热流密度的连续性，特别适用于复杂边界条件下的计算。此外，结合SOUFVM对动量方程的高阶离散，所构建的流-热耦合方法在数值精度、稳定性和实现效率上均具显著优势，适用于油浸式电力变压器温升等多物理场耦合计算。
2）以一台额定容量为321.1 MV·A、额定电压为
kV的油浸式电力变压器为研究对象，建立二维单分区不分匝绕组仿真模型，并采用所提方法对其油流速度与温升分布进行计算。将计算结果与FOUFVM-GFEM方法及Fluent软件的结果进行对比，结果表明：所提方法在计算精度上优于FOUFVM-GFEM方法，且在油流速度和热点温度的计算中，与Fluent计算结果的相对误差分别为1.38%和0.6%，均处于合理范围内。这些结果验证了所提方法在油流与温升计算中的有效性与准确性。
3）本文提出的方法有效克服了传统方法在处理复杂几何结构与流固耦合问题时的局限性，为油浸式电力变压器绕组温升分析提供了一种高效可靠的计算工具。此外，该方法为流热耦合计算方法的优化与应用提供了新的思路，同时为实际工程中的绕组温升评估与故障预测提供了关键技术支持，具有显著的工程应用潜力。
参考文献
[1] Li Yan, Meng Tiannan, Li Ziwei, et al. Research on the influence of uneven temperature distribution of transformer windings on short circuit strength[J]. IEEE Transactions on Applied Superconductivity, 2024, 34(8): 5501604.
[2] Liu Xingmou, Yang Yongming, Yang Fan, et al. Numerical research on the losses characteristic and hot-spot temperature of laminated core joints in transformer[J]. Applied Thermal Engineering, 2017, 110: 49-61.
[3] Mikhak-Beyranvand M, Faiz J, Rezaei-Zare A, et al. Electromagnetic and thermal behavior of a single- phase transformer during Ferroresonance considering hysteresis model of core[J]. International Journal of Electrical Power & Energy Systems, 2020, 121: 106078.
[4] 闻新, 战泓廷, 王戬, 等. 基于改进PSO-LSTM模型的变压器绕组热点温度预测研究[J]. 电气工程学报, 2025, 20(5): 352-361.
Wen Xin, Zhan Hongting, Wang Jian, et al. Research on transformer winding hot spot temperature pre- diction based on improved PSO-LSTM model[J]. Journal of Electrical Engineering, 2025, 20(5): 352- 361.
[5] Feng Xinyan, Zhang Da, Zhao Tingzhi, et al. Magnetic-thermal-flow field coupling simulation of oil-immersed transformer[J]. IEEE Access, 2024, 12: 65462-65470.
[6] Abdali A, Abedi A, Masoumkhani H, et al. Magnetic- thermal analysis of distribution transformer: vali- dation via optical fiber sensors and thermography[J]. International Journal of Electrical Power & Energy Systems, 2023, 153: 109346.
[7] Nicola M, Nicola C I, Sacerdoţianu D, et al. Monitoring system for power transformer windings hot spot temperature using fiber optic sensors, Kalman filter and integration in SCADA system[J]. American Journal of Signal Processing, 2018, 8(2): 33-44.
[8] Tang Pengfei, Zhang Zhonghao, Tong Jie, et al. Predicting transformer temperature field based on physics-informed neural networks[J]. High Voltage, 2024, 9(4): 839-852.
[9] Chandran L R, Ajith Babu G S, Nair M G, et al. A review on status monitoring techniques of transformer and a case study on loss of life calculation of distribution transformers[J]. Materials Today: Pro- ceedings, 2021, 46: 4659-4666.
[10] IEEE Guide for Loading Mineral-Oil-Immersed Transformers: C57.91—1995[S]. Institute of Electrical and Electronics Engineers, 1995.
[11] Chen Weigen, Pan Chong, Yun Yuxin. Power transformer top-oil temperature model based on thermal-electric analogy theory[J]. European Transa- ctions on Electrical Power, 2009, 19(3): 341-354.
[12] 刘刚, 胡万君, 郝世缘, 等. 油浸式变压器绕组瞬态温升降阶快速计算方法[J]. 电工技术学报, 2024, 39(3): 643-657.
Liu Gang, Hu Wanjun, Hao Shiyuan, et al. Reduced order calculation method of steady temperature rise of oil immersed power transformer[J]. Transactions of China Electrotechnical Society, 2024, 39(3): 643- 657.
[13] 谭又博, 余小玲, 臧英, 等. 谐波电流对换流变压器绕组损耗及温度分布特性的影响[J]. 电工技术学报, 2023, 38(2): 542-553.
Tan Youbo, Yu Xiaoling, Zang Ying, et al. The influence of harmonic current on the loss and temperature distribution characteristics of a converter transformer winding[J]. Transactions of China Elec- trotechnical Society, 2023, 38(2): 542-553.
[14] 邓永清, 阮江军, 董旭柱, 等. 基于流线分析的10kV油浸式变压器绕组热点温度反演模型建立及验证研究[J]. 中国电机工程学报, 2023, 43(8): 3191- 3204.
Deng Yongqing, Ruan Jiangjun, Dong Xuzhu, et al. Establishment and verification of 10kV oil immersed transformer winding hot spot temperature inversion model based on streamline analysis[J]. Proceedings of the CSEE, 2023, 43(8): 3191-3204.
[15] Xiao Jiayi, Zhang Zhanlong, Hao Yuefeng, et al. Method for measuring thermal flow field distribution in oil-immersed transformer using dynamic heat transfer coefficient[J]. IEEE Transactions on Instru- mentation and Measurement, 2024, 73: 9002011.
[16] 王永强, 马伦, 律方成, 等. 基于有限差分和有限体积法相结合的油浸式变压器三维温度场计算[J]. 高电压技术, 2014, 40(10): 3179-3185.
Wang Yongqiang, Ma Lun, Lü Fangcheng, et al. Calculation of 3D temperature field of oil immersed transformer by the combination of the finite element and finite volume method[J]. High Voltage Engin- eering, 2014, 40(10): 3179-3185.
[17] Star S K, Sanderse B, Stabile G, et al. Reduced order models for the incompressible Navier-Stokes equations on collocated grids using a ‘discretize-then-project’ approach[J]. International Journal for Numerical Methods in Fluids, 2021, 93(8): 2694-2722.
[18] Nie J H, Li Zengyao, Wang Qiuwang, et al. A method for viscous incompressible flows with a simplified collocated grid system[C]//Proceeding of Proceedings of Symposium on Energy Engineering in the 21st Century, Hong Kong, China, 2023: 177-183.
[19] Ferziger J H, Perić M. Computational methods for fluid dynamics[M]. New York: Springer, 2002.
[20] Bharadwaj A S. Addressing the checkerboard problem in an Eulerian meshless method for incompressible flows[J]. arXiv Preprint, 2023: 2307.09778.
[21] 刘刚, 靳艳娇, 马永强, 等. 基于非平均热源的油浸式变压器2维温度场分析[J]. 高电压技术, 2017, 43(10): 3361-3370.
Liu Gang, Jin Yanjiao, Ma Yongqiang, et al. Two- dimensional temperature field analysis of oil- immersed transformer based on non-uniformly heat source[J]. High Voltage Engineering, 2017, 43(10): 3361-3370.
[22] Hendriana D, Bathe K J. On upwind methods for parabolic finite elements in incompressible flows[J]. International Journal for Numerical Methods in Engineering, 2000, 47(1/2/3): 317-340.
[23] 彭丽丹. 电力变压器温度场数值计算研究[D]. 北京: 华北电力大学, 2016.
Peng Lidan. The numerical calculation research on the temperature field of power transformer[D]. Beijing: North China Electric Power University, 2016.
[24] 刘刚, 郝世缘, 胡万君, 等. 基于子循环自适应串行交错时间匹配算法的油浸式变压器绕组瞬态温升计算[J]. 电工技术学报, 2024, 39(4): 1185-1197.
Liu Gang, Hao Shiyuan, Hu Wanjun, et al. Transient temperature rise calculation of oil immersed trans- former winding based on sub cyclic adaptive staggered time matching algorithm[J]. Transactions of China Electrotechnical Society, 2024, 39(4): 1185- 1197.
[25] 谢裕清. 油浸式电力变压器流场及温度场耦合有限元方法研究[D]. 北京: 华北电力大学, 2017.
Xie Yuqing. Study on flow field and temperature field coupling finite element methods in oil-immersed power transformer[D]. Beijing: North China Electric Power University, 2017.
[26] 明平剑, 张文平. 计算多物理场: 有限体积方法应用[M]. 北京: 北京航空航天大学出版社, 2015.
[27] Moukalled F, Mangani L, Darwish M. The Finite Volume Method in Computational Fluid Dynamics[M]. Cham: Springer International Publishing, 2016.
[28] Khosla P K, Rubin S G. A diagonally dominant second-order accurate implicit scheme[J]. Computers & Fluids, 1974, 2(2): 207-209.
[29] Ghia U, Ghia K N, Shin C T. High-Re solutions for incompressible flow using the Navier-Stokes equations and a multigrid method[J]. Journal of Computational Physics, 1982, 48(3): 387-411.
[30] 刘刚, 高成龙, 胡万君, 等. 基于鲸鱼优化算法超参数优化的径向基函数响应面模型的油浸式变压器绕组挡板结构优化[J]. 电工技术学报, 2024, 39(17): 5331-5343.
Liu Gang, Gao Chenglong, Hu Wanjun, et al. Optimization of winding block washer structure for oil immersed transformers based on radial basis function response surface model with whale optimi- zation algorithm hyper-parameters optimization[J]. Transactions of China Electrotechnical Society, 2024, 39(17): 5331-5343.
[31] 陶文铨. 数值传热学[M]. 2版. 西安: 西安交通大学出版社, 2001.
[32] Zienkiewicz O C, Taylor R L, Nithiarasu P. The Finite Element Method for Fluid Dynamics[M]. 7th ed. Amsterdam: Elsevier, 2014
[33] 王泽忠. 简明电磁场数值计算[M]. 北京: 机械工业出版社, 2011.
[34] 王立, 喻高明, 傅宣豪, 等. 基于反距离加权插值法的产量劈分新方法[J]. 断块油气田, 2018, 25(5): 617-621.
Wang Li, Yu Gaoming, Fu Xuanhao, et al. New method for production cleavage by inverse distance weighted interpolation[J]. Fault-Block Oil & Gas Field, 2018, 25(5): 617-621.
[35] 刘刚, 胡万君, 刘云鹏, 等. 降阶技术与监测点数据融合驱动的油浸式变压器绕组瞬态温升快速计算方法[J]. 电工技术学报, 2024, 39(19): 6162- 6174.
Liu Gang, Hu Wanjun, Liu Yunpeng, et al. A fast calculation method for transient temperature rise of oil immersed transformer windings driven by fusion of order reduction technology and monitoring point data[J]. Transactions of China Electrotechnical Society, 2024, 39(19): 6162-6174.
[36] 宋浩永, 黄青丹, 陈于晴, 等. 110 kV环保型天然酯绝缘油变压器绕组的温度场分析[J]. 高压电器, 2024, 60(5): 117-123.
Song Haoyong, Huang Qingdan, Chen Yuqing, et al. Temperature field analysis of windings of natural ester insulating oil-immersed transformer 110 kV environment friendly[J]. High Voltage Apparatus, 2024, 60(5): 117-123.
[37] 刘刚, 郝世缘, 朱章宸, 等. 基于动态模态分解-自适应变步长油浸式电力变压器绕组瞬态温升快速计算方法[J]. 电工技术学报, 2024, 39(12): 3895- 3906.
Liu Gang, Hao Shiyuan, Zhu Zhangchen, et al. Research on rapid calculation method of transient temperature rise of winding of dynamic mode decomposition-adaptive time stepping oil-immersed power transformer[J]. Transactions of China Electro- technical Society, 2024, 39(12): 3895-3906.
[38] 谢裕清, 李琳, 宋雅吾, 等. 油浸式电力变压器绕组温升的多物理场耦合计算方法[J]. 中国电机工程学报, 2016, 36(21): 5957-5965, 6040.
Xie Yuqing, Li Lin, Song Yawu, et al. Multi-physical field coupled method for temperature rise of winding in oil-immersed power transformer[J]. Proceedings of the CSEE, 2016, 36(21): 5957-5965, 6040.
Abstract With the continuous development of power systems, oil-immersed power transformers, as key power equipment, are experiencing rising winding temperatures, which has become a core factor affecting their safe and stable operation. Studies show that excessively high hotspot temperatures in windings can directly cause transformer failures. Therefore, accurately calculating the windings' temperature distribution and hotspot temperatures is crucial. Traditional temperature-rise calculation methods often struggle to balance numerical accuracy, stability, and computational complexity when dealing with the coupled, complex flow and temperature fields inside transformers.
Combined with the second-order upwind finite volume method (SOUFVM) and the streamline upwind finite element method (SUFEM), this paper proposes a coupled flow-heat calculation method to eliminate the non-physical oscillations arising from the classical Galerkin finite element method when solving convection- diffusion equations. To avoid additional interpolation operations in staggered grids, the second-order upwind finite-volume method on structured collocated grids is used to solve the flow-field control equations, thereby obtaining the oil flow velocity distribution. The streamline upwind finite element method is employed to solve the unified temperature field equation for oil flow and winding heat transfer, thereby obtaining the temperature rise distribution. The proposed coupling calculation method is applied to a two-dimensional single-region winding simulation model of an oil-immersed power transformer, with a rated capacity of 321.1 MV·A and a rated voltage of 530/
kV. The distributions of oil flow and temperature rise are calculated. The results indicate that the proposed method achieves higher numerical accuracy than the Galerkin finite element method. Compared with fluent software, the absolute errors (relative errors) for the peak oil flow velocity and the hotspot temperature are 0.000 8 m/s (1.38%) and 2.06 K (0.6%), respectively. Overall, the relative error of this method remains within 2%.
In summary, the coupled flow-heat calculation method proposed in this paper integrates the advantages of finite volume and finite element methods, overcoming the shortcomings of existing numerical methods in handling complex geometries and flow-structure coupling problems. This method effectively solves the problem of calculating temperature rise in transformer windings and also provides strong technical support for other complex multiphysics coupling problems. With the continuous advancement of computer technology, this method has broad application prospects in the multiphysics simulation software for power equipment. It is expected to play an essential role in promoting the intelligent and accurate development of high-voltage power equipment simulation software.
keywords：Oil-immersed power transformer, winding temperature rise, finite volume method, finite element method, flow field, temperature field
DOI: 10.19595/j.cnki.1000-6753.tces.250520
中图分类号：TM411
国家重点研发计划资助项目（2021YFB2401700）。
收稿日期 2025-04-01
改稿日期 2025-04-24
辛纪威 男，1995年生，博士研究生，研究方向为电气设备多物理场耦合计算。E-mail: xin291573287@vip.qq.com（通信作者）
李 琳 男，1962年生，教授，博士生导师，研究方向为电磁场理论及应用与先进输电技术。E-mail: lilin@ncepu.edu.cn
（编辑 崔文静）