变压器是电力系统中极为关键的设备之一, 起到电压变换、电能传输与分配、电压隔离与安全保护以及电压调节与系统稳定等作用. 没有变压器, 就无法实现电力的远距离高效输送Li2020. 因此保证变压器的安全稳定运行是重中之重, 我们需要能够即时对变压器状态进行监测. 多物理场耦合是指电磁场、温度场、流体场等物理场之间相互作用与反馈的综合效应, 对变压器而言, 其内部结构复杂, 物理场之间相互影响, 进行单一场分析往往无法预测其真实状态, 而多物理场耦合建模与数值计算能够反映变压器内部的能量转换与传递. 文献[2]采用三维有限元方法建立了电磁–热耦合模型, 分析了油道结构对变压器内部温度分布的影响, 为冷却结构设计与局部过热风险提供了理论依据. 文献[3]通过建立并计算多物理场耦合模型, 揭示了电磁损耗产生热量后导热与散热的传递规律, 对冷却结构优化与绝缘寿命评估具有重要意义. 文献[4]系统回顾了油浸式变压器中共轭传热过程的建模方法, 指出电–热–流多物理场耦合是准确预测内部温度分布与热点形成机制的关键.

因此建立与计算多物理场耦合模型十分重要, 文献[4]系统性地对多物理场模型建模方法进行了介绍. 多物理场耦合模型计算十分复杂, 不仅源自物理场之前的相互影响, 而且也源于单场计算的不易, 电磁场方程的计算涉及非线性Maxwell方程, curl-curl算子病态性强. 并且由于变压器内部材料特性不同, 还有磁滞效应、涡流效应等存在, 求解Maxwell方程需要进行间接迭代求解, 每次迭代都需要求解大型稀疏系统. 此外铁心具有非线性的B-H曲线, 使得数值迭代容易不收敛. 文献[5]采用了3D时谐有限元方法, 建立了精确的漏磁场-涡流损耗模型, 用以分析三相电力变压器中产生的涡流损耗. 文献[6]构建了三维变压器电磁场模型, 求解了绕组内部电流与磁场分布, 计算绕组所受径向力与轴向力. 此外, 在文献[7, 8]均考虑了电力变压器漏磁场与损耗计算. 变压器内温度场计算必须要与流场进行耦合, 导致计算必须要采用稳态或者瞬态耦合求解, 又由于热源分布不均匀, 变压器材料参数强依赖温度, 导致计算中必须引入更多非线性, 提高计算的不稳定性, 并且由于冷却方式众多, 每种冷却方式都对应不同油流分布与温升模式, 导致系统更加复杂. 文献[9]基于实验测温与有限元导热模型相结合的方法, 系统研究了三相叠片式软磁变压器的热特性, 验证简化模型在热分析中的有效性. 文献[10]提出了一种基于物理约束神经网络(PINN)的油浸式变压器温度场预测方法, 通过引入物理约束使网络能够重构变压器内部的温度分布, 为数字孪生、状态监测和智能运维提供了一种高效可行的温度场重建手段. 文献[11]提出了一个基于降阶建模(reduced-order modeling, ROM)的快速温度场预测方法, 通过从高保真CFD/FEM热场数据中提取主导模式, 建立可在秒级速度推算热点温度和内部温度分布的高效模型.

相较于单场计算, 多物理场耦合计算除了遭遇单场会遇见的难题, 还需要考虑物理场之间的相互影响, 其耦合关系高度非线性、跨尺度、跨材料、跨区域. 并且由于变压器主要材料的物性参数均随温度显著变化, 导致PDE系统系数随解不断更新, 计算更加困难. 文献[3]建立了电磁场-流-热场双向耦合模型, 通过分析绕组在不同油速的温升变化, 发现电磁-热双向耦合会改变热点位置与绕组铜损密度, 而单向模型则不足以刻画. 文献[12]建立电路-磁场-热场耦合模型, 能够准确计算绕组短路导致的损耗与温升. 文献[13]考虑谐波情况下, 构建3D电磁–热–流体耦合仿真模型, 通过有限元获取绕组和铁芯损耗, 再作为热源驱动油流与温度场模拟, 分析热点位置和温度在不同载流条件下的变化. 文献[14]构建有限元模型包含电磁、热传导以及热应力耦合, 分析三相不平衡下绕组和铁芯损耗. 目前多物理场耦合计算尚不成熟, 但准确估计变压器运行状态对变压器技术更新、日常维护以及故障检测等都有十分重要的作用.

有限元方法(Finite Element Method, FEM)是一种求解偏微分方程(PDE)及变分问题的数值方法, 其核心思想是：将复杂区域划分为由简单几何形状组成的网格(如三角形、四边形、四面体等), 并在每个单元上构造局部的低维近似空间, 从而将原连续问题转化为有限维代数方程组求解. 有限元方法具有清晰的数学理论基础、良好的可扩展性以及对复杂几何、高阶精度与多物理场耦合问题的天然适应性, 因此在弹性力学、电磁场、流体力学、声学以及材料科学等领域得到广泛应用. 在本文中, 我们研究油浸式变压器的电磁-热-流耦合模型的有限元方法. 我们采用间接耦合的方法, 借助电磁场的AVA格式, 将电磁-热-流场耦合模型简化为电磁场与热-流体场的迭代模型.

在变压器内部的复杂几何结构中, 各材料的物理参数往往存在显著差异. 在不同介质界面处, PDE解通常呈现较强的奇异性, 此时数值解的误差会显著依赖于材料参数之间的比例：材料对比度越高, 误差越大. 为了有效克服由界面奇性引起的数值困难, 我们拟采用自适应有限元方法, 以在局部自动加密网格、提高计算精度并降低整体计算代价. 自适应有限元方法在求解椭圆型方程时表现出卓越的性能, 能够有效捕捉解的奇异结构并在局部进行网格加密, 从而实现误差与自由度之间的最优收敛性质. 对于抛物型方程, 自适应策略不仅要求在空间上进行细化或粗化, 还需要对时间步长进行自适应调整, 以同时控制时间离散误差和空间离散误差, 实现整体计算效率与精度的平衡Chen2011, Chen2015.

在本文中, 第二节将介绍具体构筑的多物理场耦合模型, 包括几何模型与对应多物理场方程, 而后介绍耦合模型的解耦方法并给出对应基本算法框架. 第三节介绍方程的弱形式以及针对计算的一些处理. 第四节采取有限元近似方法与向后Euler方法对方程进行离散, 基于离散格式给出方程对应的具体的自适应算法, 而后将算法结合, 给出多物理场耦合模型的自适应算法, 最后进行总结.

本节我们将介绍变压器耦合模型的数学方程以及相应的解耦格式.

令$ \Omega $是有界区域, 引入标准Sobolev空间$ W^{k,p}(\Omega) $ (参考[17]), 记$ H^p(\Omega)=W^{2,p}(\Omega) $和$ L^p(\Omega)=W^{0,p}(\Omega) $. 我们用$ (\cdot,\cdot) $表示$ L^2(\Omega) $空间上的内积, 即$ (u,v)_{\Omega} = \int_{\Omega} uv {\rm d}x $, 下标用以表示积分区域, 用$ \langle\cdot,\cdot\rangle_{\partial\Omega} $表示边界$ \partial\Omega $上的内积. 令$ H(\operatorname{curl};\Omega) = \{v\in L^2(\Omega)^3 : \operatorname{curl} v \in L^2(\Omega)^3\}, H_0(\operatorname{curl};\Omega) = H(\operatorname{curl};\Omega) \cap C_0^\infty(\Omega)^3. $我们引入标准的Bochner空间[18], 比如$ L^2(0,T;X) $:

$ \begin{equation} L^2(0,T;X) = \{ u:(0,T)\to X :\int^T_0 ||u(t)||^2_X {\rm d} t < \infty \} \notag. \end{equation} $

当$ X $是标准Sobolev空间时, 譬如$ X = L^2(\Omega) $,

$ \begin{equation} L^2(0,T;L^2(\Omega)) = \{u:(0,T)\to L^2(\Omega) :\int^T_0 ||u(t)||^2_{L^2(\Omega)} {\rm d} t < \infty \}. \notag \end{equation} $

对应时间上的$ H^1 $型Sobolev空间表示为:

$ \begin{equation} H^1(0,T;H^{-1}(\Omega)) = \{u \in L^2(0,T;H^{-1}(\Omega)): u_t \in L^2(0,T;H^{-1}(\Omega))\}. \notag \end{equation} $

为了方便叙述, 在不影响理解且不关心函数关于时间$ t $的光滑性的前提下, 我们会直接使用Sobolev空间的符号表示Bochner空间, 而省略时间部分的描述, 比如：

$ \begin{align*} &u \in H^1(\Omega) \;\;\mbox{表示}\;\; u \in \left\{v \in L^2(0,T;H^1(\Omega)) : \partial_t v \in L^2(0,T;H^{-1}(\Omega))\right\}, \\ &{{\boldsymbol{u}}} \in H({\rm curl},\Omega) \;\;\mbox{表示}\;\; {{\boldsymbol{u}}} \in \left\{{\boldsymbol{v}} \in L^2(0,T;H({\rm curl},\Omega) : \partial_t {{\boldsymbol{v}}} \in L^2(0,T;H^{-1}({\rm div}, \Omega))\right\}. \end{align*} $

我们基于COMSOL提供的单相模型[19], 考虑单相变压器物理模型, 其几何模型的平面剖分图如图 1, 其中区域$ \Omega \in \mathbb{R}^3 $由铁芯$ \Omega_{\rm core} $、绕组$ \Omega_{\rm coil} $以及油域$ \Omega_{\rm oil} $组成.

油浸式变压器多物理场耦合模型的方程组可写为：

$ \begin{equation} \begin{aligned} &\frac{\partial \mathit{\boldsymbol{B}}}{\partial t} + \nabla \times \mathit{\boldsymbol{E}} = 0, \\ &\nabla \times \mathit{\boldsymbol{H}} = \mathit{\boldsymbol{J}}_{\rm ext} + \mathit{\boldsymbol{D}}, \\ &\nabla \cdot \mathit{\boldsymbol{B}} = 0, \;\;\; {\mathit{\boldsymbol{B}}} = \mu(T) {\mathit{\boldsymbol{H}}}, \\ &\nabla \cdot \mathit{\boldsymbol{D}} = 0, \;\;\; {\mathit{\boldsymbol{D}}} = \sigma(T) {\mathit{\boldsymbol{E}}}, \\ &\nabla \cdot \mathbf{u} = 0, \\ &\rho\left(\frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u}\cdot\nabla)\mathbf{u}\right) = -\nabla p + \mu \nabla^2 \mathbf{u} + \rho \mathbf{g} \beta(T-T_{\rm ref}), \\ &\rho c_p \left(\frac{\partial T}{\partial t} + \mathbf{u}\cdot \nabla T\right) = \nabla\cdot(k\nabla T) + Q({\mathit{\boldsymbol{E}}}). \end{aligned} \end{equation} $

(2.1)

方程组的前四个方程表示准静态的麦克斯韦(Maxwell) 方程组及相应的本构关系, 用来描述电场$ \mathit{\boldsymbol{E}} $和磁场$ \mathit{\boldsymbol{H}} $, 其中$ \mu $和$ \sigma $分别表示磁导率和电导率, 取决于温度$ T $. $ \mathit{\boldsymbol{J}}_{\rm ext} $是外加电流密度源项, 满足

$ \mathit{\boldsymbol{J}}_{\rm ext} = \left\{ \begin{aligned} &\mathit{\boldsymbol{t}}\frac{n_c}{S_c} I(t) &\text{in }\Omega_{\rm coil}, \\ &0 &\text{otherwise}. \end{aligned}\right. $

其中, $ {\mathit{\boldsymbol{t}}} $为与绕组相切的单位向量, $ n_c $为匝数, $ S_c $为绕组的总横截面积, $ I(t) $为每匝电流. 方程组(2.1) 的第五、六个方程是经典的Navier-Stokes方程, 用来描述变压器油的运动状态. 其中, 第五个方程是质量守恒方程, $ \mathbf{u} $是流体速度向量, 该方程表示流体不可压缩. 第六个方程是动量守恒方程, 其中$ \rho $表示流体密度, $ p $表示流体压强, $ \mu \nabla^2\mathbf{u} $是黏性扩散项, $ \rho {\mathit{\boldsymbol{g}}} \beta(T-T_{\rm ref}) $是采用了Boussinesq近似的浮力项, 这里$ {\mathit{\boldsymbol{g}}} $是重力加速度, $ \beta $是体膨胀系数, $ T $表示流体温度, $ T_{\rm ref} $是参考温度. 方程组(2.1) 的最后一个方程是温度场方程, 其中$ c_p $是定压比热容, $ k $是导热系数, $ Q $表示体积热源, 由电场$ {\mathit{\boldsymbol{E}}} $给出.

变压器耦合模型(2.1) 是一个高度复杂的非线性多物理场系统, 其数值计算面临多方面挑战. 首先, 模型中涉及的物理量众多, 且非线性程度较高, 增加了求解的难度. 其次, 计算区域通常具有复杂几何结构, 对网格生成与求解器性能提出较高要求. 此外, 各物理场存在显著的时间与空间尺度差异, 例如电磁场变化迅速, 而温度场变化相对缓慢, 这导致不同方程组之间的耦合具有强烈的多尺度特性. 最后, 模型中不同材料的物理参数往往具有高对比度, 使得数值解在材料界面附近容易出现奇异性, 从而进一步增加了求解的难度. 在本节, 我们先给出各个物理场所对应的具体数学方程以及边界条件, 并给出耦合模型的解偶方法.

Maxwell方程组本身是全量的矢量偏微分方程, 电场$ {\mathit{\boldsymbol{E}}} $和磁场$ {\mathit{\boldsymbol{H}}} $相互耦合, 为了消除其耦合影响, 工程上常采用AVA格式. AVA格式采用势函数形式, 将部分耦合关系通过势的定义自动满足, 使得实际求解的方程组结构更简单, 减少直接处理旋度耦合带来的困难. 引入磁矢势$ \mathit{\boldsymbol{A}} $和电标势$ V $：$ \mathit{\boldsymbol{B}} = \nabla \times \mathit{\boldsymbol{A}}, \mathit{\boldsymbol{E}} = -\frac{\partial \mathit{\boldsymbol{A}}}{\partial t} - \nabla V. $将其代入Maxwell方程组, 并结合库伦规范$ \nabla \cdot \mathit{\boldsymbol{A}}=0 $, 可得AVA方程：

$ \begin{equation} \left\{ \begin{aligned} &\nabla \times (\nu\nabla\times \mathit{\boldsymbol{A}}) - \nabla(\nu \nabla\cdot\mathit{\boldsymbol{A}}) + \sigma \left(\frac{\partial \mathit{\boldsymbol{A}}}{\partial t} + \nabla V\right) = 0 & &\text{in }\Omega_{\rm core}, \\ &\nabla \cdot \left(\sigma\frac{\partial \mathit{\boldsymbol{A}}}{\partial t} + \sigma\nabla V\right) = 0 & &\text{in }\Omega_{\rm core}, \\ &\mathit{\boldsymbol{n}} \cdot \left(\sigma\frac{\partial \mathit{\boldsymbol{A}}}{\partial t} + \sigma\nabla V\right) = 0 & &\text{on } \partial\Omega_{\rm core}, \\ &\mathit{\boldsymbol{n}}_1 \times (\nu_1 \nabla \times \mathit{\boldsymbol{A}}_1) + \mathit{\boldsymbol{n}}_2 \times (\nu_2 \nabla \times \mathit{\boldsymbol{A}}_2) = 0 & &\text{on }\partial\Omega_{\rm core}, \\ &\nu_1 (\nabla\cdot\mathit{\boldsymbol{A}}_1) \mathit{\boldsymbol{n}}_1 + \nu_2 (\nabla\cdot\mathit{\boldsymbol{A}}_2) \mathit{\boldsymbol{n}}_2 = 0 & &\text{on }\partial\Omega_{\rm core}, \\ &\nabla\times(\nu\nabla\times\mathit{\boldsymbol{A}}) - \nabla(\nu \nabla\cdot\mathit{\boldsymbol{A}}) = \mathit{\boldsymbol{J}}_{\rm ext} & &\text{in }\Omega\setminus\overline{\Omega}_{\rm core}, \\ &\mathit{\boldsymbol{A}}\times \mathit{\boldsymbol{n}} = 0 & &\text{on }\partial\Omega \end{aligned}\right. \end{equation} $

(2.2)

其中, $ \nu=\mu^{-1} $表示电阻率. 在AVA格式中, $ {\mathit{\boldsymbol{A}}} $负责处理磁场旋度结构, $ V $负责处理电场的梯度结构. 在数学模型中显式地加上库伦规范约束, 可以使得$ {\mathit{\boldsymbol{A}}} $的零空间非常小(只剩常数场等与物理无关的部分), 从而不会造成伪解.

通过电磁场方程计算电磁损耗后, 这些损耗以体积热源分布的形式传递给热场, 成为温升分析的输入. 温度场-流场耦合模型由以下不可压缩的Navier-Stokes方程[20] 和对流扩散方程[21]联合描述：

$ \begin{equation} \left\{ \begin{aligned} &\nabla \cdot \mathbf{u} = 0& &\text{in }\Omega_{\rm oil}, \\ &\rho\left(\frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u}\cdot\nabla)\mathbf{u}\right) = -\nabla p + \mu \nabla^2 \mathbf{u} + \rho \mathbf{g} \beta(T-T_{\rm ref}), & &\text{in }\Omega,\\ &\mathbf{u} = 0 & &\text{on }\partial\Omega_{\rm oil}, \\ & \int_{\Omega_{\rm oil}} p \ {\rm{d}} x = 0, \\ &\rho c_p \left(\frac{\partial T}{\partial t} + \mathbf{u}\cdot \nabla T\right) = \nabla\cdot(k\nabla T) + Q({\mathit{\boldsymbol{E}}})& & \text{in }\Omega, \\ &T_{\rm oil} = T_i, & &\text{on }\Gamma_{i}, \\ &-k_{\rm oil}\nabla T_{\rm oil} \cdot {\mathit{\boldsymbol{n}}} = -k_i \nabla T_i \cdot {\mathit{\boldsymbol{n}}} & &\text{on }\Gamma_{i},\\ &-k\frac{\partial T}{\partial n} = h(T^4-T_{\rm ref}^4) + \lambda(T-T_{\rm ref}) & &\text{on }\partial \Omega \end{aligned} \right. \end{equation} $

(2.3)

其中$ \Gamma_1 = \partial\Omega_{\rm oil} \cap \partial\Omega_{\rm core},\,\Gamma_2 = \partial\Omega_{\rm oil} \cap \partial\Omega_{\rm coil} $. 流体速度场$ {\mathit{\boldsymbol{u}}} $只定义在油域, 但注意到其满足无滑移边界条件$ {\mathit{\boldsymbol{u}}}|_{\partial\Omega_{\rm oil}}=0 $, 因此可以延拓定义$ {\mathit{\boldsymbol{u}}}|_{\Omega\setminus\Omega_{\rm oil}} = 0 $, 从而有$ {\mathit{\boldsymbol{u}}} \in H_0^1(\Omega)^3 $. 另外, 温度场在界面$ \Gamma_i\,(i=1,2) $处的方程表示连续性条件, 在外边界$ \partial\Omega $的边界条件由对流项$ \lambda(T-T_{\rm ref}) $和辐射项$ h(T^4-T_{\rm ref}^4) $联合给出.

由于物理量太多, 直接求解耦合模型(2.1) 几乎不可行, 在实际计算时不仅需要占用大量计算资源, 而且计算时间长. 为此, 我们选择采用间接耦合的方法进行求解. 我们采用Picard线性化的方法对方程进行处理, 可以注意到将多物理场耦合方程中非线性项线性化后, 电磁场方程中耦合项自然解耦, 我们能够将多物理场方程解耦为电磁场方程和热-流方程的迭代求解, 见算法1. 其核心思想是在每个热时间步内迭代更新电磁参数, 求解电磁场并计算损耗, 再将其作为热源输入热–流体方程, 直至温度场收敛. 求解电场$ {\mathit{\boldsymbol{E}}}(x,t) $的初始值一般选为上一个热时间步得到的电场, 电场和温度场、流场的边界条件根据实际情况而定.

Algorithm 1  多物理场耦合计算
1:  给定初始温度$ T_0 $以及热场时间最大步数$ N_{\rm th} $, 并记$ j=1 $；
2:  while $ j \leq N_{\rm th} $ do
3:      令$ T_{j}=T_{j-1} $
4:      repeat
5:          令$ T=T_{j} $, 更新电磁参数$ \nu(T) $和$ \sigma(T) $;
6:          求解Maxwell方程组, 得到电场$ {\mathit{\boldsymbol{E}}}(x,t) $；
7:          计算电磁损耗得到热源$ Q({\mathit{\boldsymbol{E}}}) $;
8:          求解温度场–流场方程, 得到新的温度场$ T_j $;
9:      until $ |T-T_j|\leq 0.01$;
10:     $ j=j+1 $;
11:  end while
12:  输出温度场$ T=\{T_j\} $.

本节, 我们首先分别给出电磁场方程的AVA格式(2.2) 和热-流耦合方程(2.3) 的弱形式, 并引入惩罚项或者Lagrange乘子法消除常数自由度, 然后联合得到耦合模型方程的变分问题.

利用Green公式, 取任意函数$ \mathit{\boldsymbol{\phi}} \in H_0({\rm curl},\Omega) $以及$ w \in H^1(\Omega_{\rm core}) $乘上方程(2.2), 并在相应区域积分, 可得

$ \begin{equation*} \begin{aligned} & \int_{\Omega\setminus\partial\Omega_{\rm core}} \nu \big(\nabla\times\mathit{\boldsymbol{A}}\cdot\nabla\times\mathit{\boldsymbol{\phi}} + (\nabla\cdot\mathit{\boldsymbol{A}})(\nabla\cdot\mathit{\boldsymbol{\phi}})\big) {\rm\,d}x + \int_{\Omega_{\rm core}} \sigma \left(\frac{\partial\mathit{\boldsymbol{A}}}{\partial t} + \nabla V\right)\big(\mathit{\boldsymbol{\phi}}+\nabla w\big) {\rm\ d}x \\ =& \int_{\Omega_{\rm coil} }{\mathit{\boldsymbol{J}}}_{\rm ext} \cdot \mathit{\boldsymbol{\phi}} {\rm\ d}x \end{aligned} \end{equation*} $

定义双线性形式

$ \begin{equation*} \begin{aligned} a(T;\mathit{\boldsymbol{A}},V;\mathit{\boldsymbol{\phi}},w) :=&\int_{\Omega\setminus\partial\Omega_{\rm core}} \nu(T) \big(\nabla\times\mathit{\boldsymbol{A}}\cdot\nabla\times\mathit{\boldsymbol{\phi}} + (\nabla\cdot\mathit{\boldsymbol{A}})(\nabla\cdot\mathit{\boldsymbol{\phi}})\big) {\rm{d}}x \\ & +\int_{\Omega_{\rm core}} \sigma(T) \left(\frac{\partial\mathit{\boldsymbol{A}}}{\partial t} + \nabla V\right)\big(\mathit{\boldsymbol{\phi}}+\nabla w\big) {\rm d}x. \end{aligned} \end{equation*} $

为了消除$ (\mathit{\boldsymbol{A}},V) $的常数自由度, 避免线性系统奇异, 在实际计算时需要加上惩罚项. 因此, $ \mathit{\boldsymbol{A}} $-$ V $-$ \mathit{\boldsymbol{A}} $方程(2.2)的弱形式可写为：求$ (\mathit{\boldsymbol{A}},V) \in X:=H_0({\rm curl},\Omega) \times H^1(\Omega_{\rm core}) $使得

$ \begin{equation} \begin{aligned} a(T;\mathit{\boldsymbol{A}},V;\mathit{\boldsymbol{\phi}},w) +\eta \bigg(\int_{\Omega} \mathit{\boldsymbol{A}}\cdot\mathit{\boldsymbol{\phi}} \ {\rm d}x + \int_{\Omega_{\rm core}} Vw {\rm \ d}x\bigg) = \int_{\Omega_{\rm coil}} \mathit{\boldsymbol{J}}_{\rm ext}\cdot\mathit{\boldsymbol{\phi}} {\,\rm d}x , \forall\, (\mathit{\boldsymbol{\phi}},w) \in X. \end{aligned} \end{equation} $

(3.1)

其中, $ \eta $是罚参数. 我们引入罚参数时需要选取$ \eta>0 $为一个较大的数, 虽然这样会改变方程结构, 但从能量泛函的角度来考虑, 新的能量泛函相比于此前, 增加了$ \mathcal{E}_{\eta}(A,V) = \mathcal{E}_0(A,V) + \eta(\|A\|^2 +\| V\|^2) $其中$ \mathcal{E}_0(A,V) = a({\mathit{\boldsymbol{A}}},V;{\mathit{\boldsymbol{A}}},V) $. 从能量最小的角度来考虑, 这会使得惩罚后的解在满足物理条件的解中选择对应范数最小的那一个, 从而消除自由度, 但并不会改变满足约束的最优解. 除了惩罚法, Lagrange乘子法也常用于消除常数自由度. 将上述两个条件视作问题本身的约束, Lagrange乘子法引入新的自由度$ \kappa_{{\mathit{\boldsymbol{A}}}} \in \mathbb{R}^3, \ \kappa_V \in \mathbb{R} $, 则方程(2.2) 的变分问题可写为：求$ (\mathit{\boldsymbol{A}},V,\kappa_{{\mathit{\boldsymbol{A}}}},\kappa_V) \in \tilde X := (X,\mathbb{R}^3,\mathbb{R}) $使得:

$ \begin{equation} \begin{aligned} &a(T;{\mathit{\boldsymbol{A}}},V;\mathit{\boldsymbol{\phi}},w) + \int_{\Omega} \kappa_{{\mathit{\boldsymbol{A}}}} \cdot \boldsymbol \phi \ {\rm d}x + \int_{\Omega} \kappa_V w\ {\rm d}x + \int_{\Omega} {{\mathit{\boldsymbol{A}}}} \cdot r_{{\mathit{\boldsymbol{A}}}} \ {\rm d}x + \int_{\Omega} V r_V\ {\rm d}x \\ =& \int_{\Omega_{\rm coil}} \mathit{\boldsymbol{J}}_{\rm ext}\cdot\mathit{\boldsymbol{\phi}} {\,\rm d}x \qquad \forall\, (\mathit{\boldsymbol{\phi}},w,r_{{\mathit{\boldsymbol{A}}}},r_V) \in \tilde X. \end{aligned} \end{equation} $

(3.2)

惩罚法与Lagrange乘子法相比, 惩罚法实现更为简单, 不会改变矩阵的稀疏结构, 但人工参数$ \eta $的选择十分麻烦, 需要试算或者根据经验选取, 不然会造成矩阵近似奇异或系统病态. 而且, 由于$ \eta $选择较大, 方程结构发生改变, 数值尺度不均, 并且矩阵的条件数会依赖于$ \eta $的取值, 网格尺寸或时间步长都需要重新调整；Lagrange乘子法则能严格满足约束, 无需调整参数, 问题属于标准的鞍点结构, 但由于引入了新的未知量, 并且变为鞍点问题, 实现起来更为复杂.

取任意函数$ \mathit{\boldsymbol{v}} \in (H^1_0(\Omega_{\rm oil}))^3 $, $ q \in L^2_0(\Omega_{\rm oil}) $以及$ \theta \in H^1(\Omega) $乘上热-流体耦合方程(2.3), 其中$ L^2_0(\Omega_{\rm oil}):= \{u \in L^2(\Omega_{\rm oil}): \int_{\Omega_{\rm oil}}u \ {\rm d}x = 0\} $, 在相应区域积分并利用Green公式可得：

$ \begin{equation*} \left\{ \begin{aligned} &\rho(\partial_t\mathbf{u},{\mathit{\boldsymbol{v}}})_{\Omega_{\rm oil}} + \rho((\mathbf{u\cdot \nabla})\mathbf{u},{\mathit{\boldsymbol{v}}})_{\Omega_{\rm oil}} + \mu (\nabla \mathbf{u},\nabla {\mathit{\boldsymbol{v}}})_{\Omega_{\rm oil}} - (p,\nabla \cdot {\mathit{\boldsymbol{v}}})_{\Omega_{\rm oil}} = (\rho {\mathit{\boldsymbol{g}}} \beta(T-T_{\rm ref}),{\mathit{\boldsymbol{v}}})_{\Omega_{\rm oil}}, \\ &(\nabla \cdot {\mathit{\boldsymbol{u}}},q)_{\Omega_{\rm oil}} = 0, \\ &\rho c_p (\partial_t T,\theta)_{\Omega}+\rho c_p(\mathbf{u}\cdot \nabla T,\theta)_{\Omega_{\rm oil}} +(k\nabla T,\nabla \theta)_{\Omega} + \left<h(T^4-T^4_{\rm ref})+\lambda(T-T_{\rm ref}),\theta\right>_{\partial \Omega} = (Q,\theta)_{\Omega}. \end{aligned} \right. \end{equation*} $

这是一个非线性耦合系统, 需要线性化并解偶求解. 类似电磁场方程, 约束条件$ \int_{\Omega_{\rm oil}} p {\rm\ d}x=0 $可通过Lagrange乘子法处理. 定义双线性形式

$ \begin{align*} b(\mathbf{w};\mathbf{u},p;{\mathit{\boldsymbol{v}}}, q) :=& \rho(\partial_t\mathbf{u},{\mathit{\boldsymbol{v}}})_{\Omega_{\rm oil}} + \rho((\mathbf{w}\cdot \nabla)\mathbf{u},{\mathit{\boldsymbol{v}}})_{\Omega_{\rm oil}} + \mu (\nabla \mathbf{u},\nabla {\mathit{\boldsymbol{v}}})_{\Omega_{\rm oil}} \\ & - (p,\nabla \cdot {\mathit{\boldsymbol{v}}})_{\Omega_{\rm oil}} - (\nabla \cdot {\mathit{\boldsymbol{u}}},q)_{\Omega_{\rm oil}}, \\ c(\mathbf{u};w;T,\theta) :=& \rho c_p (\partial_t T,\theta)_{\Omega}+\rho c_p(\mathbf{u}\cdot \nabla T,\theta)_{\Omega_{\rm oil}} +(k\nabla T,\nabla \theta)_{\Omega} \\ & +\left< h w^3+\lambda) T,\theta\right>_{\partial \Omega}, \end{align*} $

则热-流耦合方程(2.3) 的弱形式为：求$ (\mathbf u, p,T) \in Y:= H_0^1(\Omega_{\rm oil})^3 \times L_0^2(\Omega_{\rm oil}) \times H^1(\Omega) $, 使得

$ \begin{equation} \begin{aligned} b(\mathbf{u},\mathbf{u},p;{\mathit{\boldsymbol{v}}},q) &= (\rho {\mathit{\boldsymbol{g}}} \beta(T-T_{\rm ref}),{\mathit{\boldsymbol{v}}})_{\Omega_{\rm oil}}, \\ c(\mathbf{u};T;T,\theta) &= (Q,\theta)_{\Omega} + \left<h T_{\rm ref}^4 + \lambda T_{\rm ref},\theta\right>_{\partial\Omega}. \end{aligned} \end{equation} $

(3.3)

联合(3.1) (或(3.2)) 和(3.3), 可得耦合模型(2.1) 的弱形式, 求$ ({\mathit{\boldsymbol{A}}},V,\mathbf{u},p,T) \in X \times Y $使得以下同时成立：

$ \begin{align} a(T;\mathit{\boldsymbol{A}},V;\mathit{\boldsymbol{\phi}},w) + \eta \left(\int_{\Omega} \mathit{\boldsymbol{A}}\cdot\mathit{\boldsymbol{\phi}} \ {\rm d}x + \int_{\Omega_{\rm core}} Vw {\rm \ d}x\right)= \int_{\Omega_{\rm coil}} \mathit{\boldsymbol{J}}_{\rm ext}\cdot\mathit{\boldsymbol{\phi}} {\,\rm d}x, \end{align} $

(3.4a)

$ \begin{align} b(\mathbf{u},\mathbf{u},p;{\mathit{\boldsymbol{v}}},q) = (\rho {\mathit{\boldsymbol{g}}} \beta(T-T_{\rm ref}),{\mathit{\boldsymbol{v}}})_{\Omega_{\rm oil}}, \end{align} $

(3.4b)

$ \begin{align} c(\mathbf{u};T;T,\theta) = (Q({\mathit{\boldsymbol{E}}}),\theta)_{\Omega} + \left<h T_{\rm ref}^4 + \lambda T_{\rm ref},\theta\right>_{\partial\Omega}, \end{align} $

(3.4c)

其中方程(3.4c) 的热源项$ Q $由电场$ {\mathit{\boldsymbol{E}}} = -\partial_t {\mathit{\boldsymbol{A}}} - \nabla V $以及损耗公式给出.

本节聚焦于变分问题(3.4) 的有限元离散方法, 并引入自适应有限元算法降低由模型几何、参数等带来的奇性.

我们用$ \mathcal{T}_h $表示对区域$ \Omega $以网格大小$ h $的拟一致单纯形剖分网格, 用$ \Omega_h $表示对区域$ \Omega $的近似, 用$ \mathcal{T}^{\rm core}_h $表示$ \mathcal{T}_h $在区域$ \Omega_{\rm core} $的限制, 用$ \Omega^{\rm core}_h $表示$ \Omega_h $在子区域$ \Omega_{\rm core} $的限制. 此外, 我们用$ K $来表示网格中的有限元单元, 用$ P_k(K) $来表示定义在单元$ K $上次数小于等于$ k $的多项式空间. 基于网格剖分, 我们可以定义所需的有限元空间$ X_h := {\mathit{\boldsymbol{N}}}_h \times S_h $:

$ \begin{equation} \begin{aligned} {\mathit{\boldsymbol{N}}}_h &:= \{{\mathit{\boldsymbol{v}}}_h \in H_0(\operatorname{curl};\Omega): {\mathit{\boldsymbol{v}}}_h|_{K} \in (P_1(K))^3 + \text{edge shape functions} , \ \ \forall K \in \mathcal{T}_h \} \\ S_h &:= \{ q_h \in H^1(\Omega_{\rm core}): q_h|_K \in P_1(K), \ \ \forall K \in \mathcal{T}_h^{\rm core} \} \end{aligned} \end{equation} $

(4.1)

详细的Nédélec元有限元空间参考文献[22]. 有限元空间$ {\mathit{\boldsymbol{N}}}_h $与$ S_h $将分别用来近似$ {\mathit{\boldsymbol{A}}} $与$ V $.

记$ D_h = \Omega_h \setminus \Omega^{\rm core}_h $. 根据弱形式(3.1), 我们可得如下半离散有限元格式：给定温度场$ T_h $, 求$ (\mathit{\boldsymbol{A}}_h,V_h) \in X_h $使得

$ \begin{equation} \begin{aligned} &\bigl(\,\nu(T)\,\nabla\times\boldsymbol A_h,\ \nabla\times\boldsymbol\phi_h\,\bigr)_{D_h} + \bigl(\,\nu(T)\,\nabla\!\cdot\boldsymbol A_h,\ \nabla\!\cdot\boldsymbol\phi_h\,\bigr)_{D_h} + \bigl(\ \sigma(T)(\partial_t \boldsymbol A_h+\nabla V_h),\ \boldsymbol\phi_h+\nabla w_h\ \bigr)_{\Omega_h^{\mathrm{core}}} \\ & + \eta\big(\,(\boldsymbol A_h,\boldsymbol\phi_h)_{\Omega_h} + (V_h,w_h)_{\Omega_h^{\mathrm{core}}}\big) = (\mathit{\boldsymbol{J}}_{\rm ext},\mathit{\boldsymbol{\phi}}_h) \quad \forall\, (\mathit{\boldsymbol{\phi}}_h,w_h) \in X_h. \end{aligned} \end{equation} $

(4.2)

记时间区间为$ I=[0,T_{\rm end}] $, 我们采用隐式Euler方法来处理时间项, 即对时间区间$ I $以步长$ \tau_i $做剖分, 从而有$ t^n = \Sigma_{i}^n\tau_i $, 以及$ {\mathit{\boldsymbol{A}}}^n = {\mathit{\boldsymbol{A}}}(x,t^n) $, 从而可得如下全离散有限元格式：

$ \begin{equation} \begin{aligned} &\bigl(\,\nu(T_h)\,\nabla\times\boldsymbol A^{n+1}_h,\ \nabla\times\boldsymbol\phi_h\,\bigr)_{D_h} + \bigl(\,\nu(T_h)\,\nabla\!\cdot\boldsymbol A^{n+1}_h,\ \nabla\!\cdot\boldsymbol\phi_h\,\bigr)_{D_h} \\ & + \bigg(\ \sigma(T_h)\Big( \frac{\boldsymbol A^{n+1}_h-\boldsymbol A^{n}_h}{\tau_n}+\nabla V^{n+1}_h\Big),\ \boldsymbol\phi_h+\nabla w_h\ \bigg)_{\Omega_h^{\mathrm{core}}} + \eta\big(\,(\boldsymbol A_h^{n+1},\boldsymbol\phi_h)_{\Omega_h} + (V^{n+1}_h,w_h)_{\Omega_h^{\mathrm{core}}}\big)\\ =& (\mathit{\boldsymbol{J}}_{\rm ext}^{n+1},\mathit{\boldsymbol{\phi}})\quad \forall\, (\mathit{\boldsymbol{\phi}}_h,w_h) \in X_h. \end{aligned} \end{equation} $

(4.3)

移项且合并同类项后可得

$ \begin{equation} \begin{aligned} &\bigl(\,\nu\,\nabla\times\boldsymbol A^{n+1}_h,\ \nabla\times\boldsymbol\phi_h\,\bigr)_{D_h} + \bigl(\,\nu\,\nabla\!\cdot\boldsymbol A^{n+1}_h,\ \nabla\!\cdot\boldsymbol\phi_h\,\bigr)_{D_h} \\ & + \bigl(\ \sigma( \tau_n^{-1}{\boldsymbol A^{n+1}_h}+\nabla V^{n+1}_h),\ \boldsymbol\phi_h+\nabla w_h\ \bigr)_{\Omega_h^{\mathrm{core}}} + \eta\Bigl[\,(\boldsymbol A_h^{n+1},\boldsymbol\phi_h)_{\Omega_h} + (V^{n+1}_h,w_h)_{\Omega_h^{\mathrm{core}}}\Bigr]\\ =& (\mathit{\boldsymbol{J}}_{\rm ext}^{n+1},\mathit{\boldsymbol{\phi}}) + (\sigma \tau_n^{-1}{\boldsymbol A^{n}_h},\boldsymbol \phi_h + \nabla w_h)_{\Omega^{\rm core}_h}. \end{aligned} \end{equation} $

(4.4)

类似地, 我们考虑热-流耦合方程的有限元方法. 基于同样的网格剖分, 我们用$ \mathcal{T}_h^{\rm oil} $来表示油域上的网格剖分, 以$ \Omega_h^{\rm oil} $来表示对$ \Omega_{\rm oil} $的近似. 定义以下有限元空间：

$ \begin{equation} \begin{aligned} {\mathit{\boldsymbol{V}}}_h &:= \{ {\mathit{\boldsymbol{v}}}_h \in H^1_0({\Omega_{\rm oil}}): {\mathit{\boldsymbol{v}}}_h|_K \in (P_2(K))^3 \ \ \forall K \in \mathcal{T}_h^{\rm oil} \}, \\ Q_h &:= \{ q_h\in L^2_0({\Omega_{\rm oil}}): q_h|_K \in P_1(K) \ \ \forall K \in \mathcal{T}_h^{\rm oil} \}, \\ W_h &:= \{ \theta_h\in H^1(\Omega): \theta_h|_K \in P_1(K) \ \ \forall K \in \mathcal{T}_h \}. \end{aligned} \end{equation} $

(4.5)

有限元空间$ {\mathit{\boldsymbol{V}}}_h $, $ Q_h $与$ W_h $将分别用来近似$ \mathbf{u} $, $ p $与$ T $. 根据(3.3), 可得如下半离散有限元问题：对任意测试函数$ ({\mathit{\boldsymbol{v}}}_h,q_h,\theta_h) \in ({\mathit{\boldsymbol{V}}}_h,Q_h,W_h) $, 求$ (\mathbf{u}_h,p_h,T_h)\in ({\mathit{\boldsymbol{V}}}_h,Q_h,W_h) $使得：

$ \begin{equation} \begin{aligned} &\rho(\partial_t\mathbf{u}_h,{\mathit{\boldsymbol{v}}}_h)_{\Omega_h^{\rm oil}} + \rho((\mathbf{u}_h\cdot \nabla)\mathbf{u}_h,{\mathit{\boldsymbol{v}}}_h)_{\Omega_h^{\rm oil}} + \mu (\nabla \mathbf{u}_h,\nabla v_h)_{\Omega_h^{\rm oil}} - (p_h,\nabla \cdot {\mathit{\boldsymbol{v}}}_h)_{\Omega_h^{\rm oil}} - (\nabla \cdot \mathbf{u}_h,q_h)_{\Omega_h^{\rm oil}} \\ =& \rho g \beta(T_h-T_{\rm ref},{\mathit{\boldsymbol{v}}}_h \cdot {\mathit{\boldsymbol{e}}}_g)_{\Omega_h^{\rm oil}}, \\ &\rho c_p(\partial_t T_h,\theta_h)_{\Omega_h}+\rho c_p(\mathbf{u}_h\cdot \nabla T_h,\theta_h)_{\Omega_h^{\rm oil}}+(k\nabla T_h,\nabla \theta_h)_{\Omega_h} +\left<h(T_h^4-T^4_{\rm ref})+\lambda(T_h-T_{\rm ref}),\theta_h\right>_{\partial \Omega_h}\\ =& (Q,\theta_h)_{\Omega_h}, \end{aligned} \end{equation} $

(4.6)

其中$ {\mathit{\boldsymbol{e}}}_{g} $为重力加速度方向的单位向量.

同样地, 在时间区间$ I = (0,T_{\rm end}) $以时间步长$ \tilde\tau_i $进行剖分, 使得$ \mathbf{u}^{n} = \mathbf{u}(x,\Sigma_i^n\tilde\tau_i) $. 采用隐式Euler方法, 并利用Picard线性化技巧处理非线性项(参考[23]), 可得

$ \begin{equation} \begin{aligned} &\rho \bigg(\frac{\mathbf{u}^{n+1}_h-\mathbf{u}_h^n}{\tilde\tau_n},{\mathit{\boldsymbol{v}}}_h\bigg)_{\Omega_h^{\rm oil}} + \rho((\mathbf{u^n_h\cdot \nabla})\mathbf{u}_h^{n+1},{\mathit{\boldsymbol{v}}}_h)_{\Omega_h^{\rm oil}} + \mu (\nabla \mathbf{u}^{n+1}_h,\nabla v_h)_{\Omega_h^{\rm oil}} \\ & - (p^{n+1}_h,\nabla \cdot {\mathit{\boldsymbol{v}}}_h)_{\Omega_h^{\rm oil}} - (\nabla \cdot \mathbf{u}_h^{n+1},q_h)_{\Omega_h^{\rm oil}} = \rho g \beta(T_h^{n+1}-T_{\rm ref},{\mathit{\boldsymbol{v}}}_h \cdot {\mathit{\boldsymbol{e}}}_g)_{\Omega_h^{\rm oil}}, \\ &\rho c_p \bigg(\frac{T_h^{n+1}-T_h^n}{\tilde\tau_n} ,\theta_h\bigg)_{\Omega_h}+\rho c_p(\mathbf{u}_h^{n}\cdot \nabla T_h^{n+1},\theta_h)_{\Omega_h^{\rm oil}}+(k\nabla T^{n+1}_h,\nabla \theta_h)_{\Omega_h} \\ &+\langle h\big((T^{n}_h)^4+4(T^{n}_h)^3(T^{n+1}_h-T^n_h)-T^4_{\rm ref}\big)+\lambda(T^{n+1}_h-T_{\rm ref}),\theta_h\rangle_{\partial \Omega_h}= (Q^{n+1},\theta_h)_{\Omega_h}. \end{aligned} \end{equation} $

(4.7)

经过移项且合并同类项后我们可以得到热-流耦合问题(3.3) 的有限元格式：

$ \begin{equation} \begin{aligned} &\rho\left(\tilde\tau_n^{-1}{\mathbf{u}^{n+1}_h},{\mathit{\boldsymbol{v}}}_h\right)_{\Omega_h^{\rm oil}} + \rho((\mathbf{u}^n_h\cdot \nabla)\mathbf{u}_h^{n+1},{\mathit{\boldsymbol{v}}}_h)_{\Omega_h^{\rm oil}} + \mu (\nabla \mathbf{u}^{n+1}_h,\nabla v_h)_{\Omega_h^{\rm oil}} \\ & - (p^{n+1}_h,\nabla \cdot {\mathit{\boldsymbol{v}}}_h)_{\Omega_h^{\rm oil}} - (\nabla \cdot \mathbf{u}_h^{n+1},q_h)_{\Omega_h^{\rm oil}}\\ =&\rho g \beta(T^{n+1}_h-T_{\rm ref},{\mathit{\boldsymbol{v}}}_h \cdot {\mathit{\boldsymbol{e}}}_g)_{\Omega_h^{\rm oil}} + \rho\left(\tilde\tau_n^{-1}{\mathbf{u}_h^n},{\mathit{\boldsymbol{v}}}_h\right)_{\Omega_h^{\rm oil}} ,\\ &\rho c_p\left(\tilde\tau_n^{-1}{T_h^{n+1}} ,\theta_h\right)_{\Omega_h}+\rho c_p(\mathbf{u}_h^{n}\cdot \nabla T_h^{n+1},\theta_h)_{\Omega_h^{\rm oil}}+(k\nabla T^{n+1}_h,\nabla \theta_h)_{\Omega_h} +\langle(4h(T^n_h)^3+\lambda) T^{n+1}_h,\theta_h\rangle_{\partial\Omega_h} \\ =&\left(Q^{n+1}+\rho c_p \tilde\tau_n^{-1} {T^n_h},\theta_h\right)_{\Omega_h}+\langle3h(T^{n}_h)^4+hT^4_{\rm ref}+\lambda T_{\rm ref },\theta_h \rangle_{\partial \Omega_h}. \end{aligned} \end{equation} $

(4.8)

需要注意的是, 在采取Picard线性化时, 边界项所含$ T^4 $是一个高度非线性项, 不能直接通过替代来进行线性化, 而是需要通过一阶Taylor展开来进行.

变压器结构复杂, 几何模型在网格剖分时通常呈现明显的尺度不均匀性, 导致不同单元上的有限元误差分布存在显著差异. 同时, 由于构成变压器的不同材料在电磁参数上存在较大差别, 数值解在材料界面附近往往表现出较强奇异性. 基于上述因素, 本研究引入自适应有限元算法, 通过后验误差估计指示子动态调整局部网格密度, 使得自由度分布与误差收敛之间达到近似最优的平衡, 从而提高整体计算效率与精度.

基于[15, 16], 我们实际处理的电磁方程的AVA形式, 在铁芯区域即为标准的抛物型方程, 在油区退化为了椭圆型方程, 因此我们只需要将误差指标的具体表达进行对应调整, 就能将其中抛物型自适应算法应用至AVA方程中, 从而实现计算目标. 我们记双线性形式$ a^{n+1}_h({\mathit{\boldsymbol{A}}}^{n+1}_h,V^{n+1}_h;\boldsymbol \phi_h,w_h) $满足：

$ \begin{equation} \begin{aligned} a^{n+1}_h({\mathit{\boldsymbol{A}}}^{n+1}_h,V^{n+1}_h;\boldsymbol \phi_h,w_h):= &\bigl(\,\nu\,\nabla\times\boldsymbol A^{n+1}_h,\ \nabla\times\boldsymbol\phi_h\,\bigr)_{D_h} + \bigl(\,\nu\,\nabla\!\cdot\boldsymbol A^{n+1}_h,\ \nabla\!\cdot\boldsymbol\phi_h\,\bigr)_{D_h} \\ & + \bigl(\ \sigma( \frac{\boldsymbol A^{n+1}_h}{\tau_n}+\nabla V^{n+1}_h),\ \boldsymbol\phi_h+\nabla w_h\ \bigr)_{\Omega_h^{\mathrm{core}}} \quad\ \\ &+ \eta\Bigl[\,(\boldsymbol A_h^{n+1},\boldsymbol\phi_h)_{\Omega_h} + (V^{n+1}_h,w_h)_{\Omega_h^{\mathrm{core}}}\Bigr]. \end{aligned} \end{equation} $

(4.9)

此前我们已经推导出第n次时间迭代步的近似为：给定$ {\mathit{\boldsymbol{A}}}^n_{h} $, 求解$ ({\mathit{\boldsymbol{A}}}^{n+1}_h,V^{n+1}_h)\in X_h $, 使得：

$ \begin{equation} \begin{aligned} a^{n+1}_h({\mathit{\boldsymbol{A}}}^{n+1}_h,V^{n+1}_h;\boldsymbol \phi_h,w_h)= (\overline{\mathit{\boldsymbol{J}}_{\rm ext}^{n+1}},\mathit{\boldsymbol{\phi}}) + (\sigma\frac{\boldsymbol A^{n}_h}{\tau_n},\boldsymbol \phi_h + \nabla w_h)_{\Omega^{\rm core}_h}, \end{aligned} \end{equation} $

(4.10)

其中$ \overline{\mathit{\boldsymbol{J}}_{\rm ext}^{n+1}} = \frac{1}{\tau_{n+1}}\int_{t^n}^{t^{n+1}} {\mathit{\boldsymbol{J}}_{\rm ext}} {\rm d}t $

我们用$ \mathcal{B}^{n+1} $来表示在$ \Omega $上剖分$ \mathcal{M}^{n+1} $中所有单元之间内部公共边的集合. 而后我们定义残差$ R^{n+1}_{\rm time} $为时间残差：

$ \begin{equation} R^{n+1}_{\rm time} = \overline{\mathit{\boldsymbol{J}}_{\rm ext}^{n+1}}-\sigma\chi_{\rm core}\frac{\boldsymbol A_h^{\,n+1}-\boldsymbol A_h^{\,n}}{\tau_{n+1}}. \end{equation} $

(4.11)

$ R^{n+1}_{\rm ell} = \Sigma_K R^{\,n+1}_{\mathrm{ell},K} $为空间残差, $ R^{\,n+1}_{\mathrm{ell},K} $为局部残差：

$ R^{\,n+1}_{\mathrm{ell},K} = \overline{\boldsymbol J_{\mathrm{ext}}^{\,n+1}} - [\nabla\times\bigl(\nu\,\nabla\times \boldsymbol A^{n+1}_h\bigr) - \nabla\bigl(\nu\,\nabla\!\cdot \boldsymbol A^{n+1}_h\bigr)]\big|_K - \sigma\,\chi_{\mathrm{core}}\, \nabla V_h^{\,n+1}\big|_K . $

$ J^{n+1} = \Sigma_e J^{n+1}_e $, 记$ K\in \mathcal{M}^{n+1} $为剖分中的单元, $ h_K $代表$ K $的外接圆直径, $ h_e $代表$ e \in \mathcal{B}^{n+1} $的尺寸. 我们定义$ e $上的跨越残差：

$ \begin{equation} \begin{aligned} &J^{n+1}_e:=[[\nu\nabla \times {\mathit{\boldsymbol{A}}}^{n+1}_h \times n_e]]_e + [[\nu \nabla \cdot {\mathit{\boldsymbol{A}}}^{n+1}_hn_e]]_e,\\ &e = \partial K_1 \cap \partial K_2. \end{aligned} \end{equation} $

(4.12)

其中$ \nu_e $是从$ K_2 $指向$ K_1 $的关于$ e $的单位法向量. 此外我们定义能量范数$ |||\boldsymbol \varphi|||^2_{\Omega} = \bigl(\,\nu\,\nabla\times\boldsymbol \varphi,\ \nabla\times\boldsymbol\varphi\,\bigr)_{D} + $. 基于以上定义, 我们可以定义时间误差指标$ \eta^n_{\rm time} $空间误差指标$ \eta^n_{\rm space} $以及局部误差指标$ \eta^n_K $：

$ \begin{equation} \begin{aligned} \eta^{n+1}_{\rm time} &= \frac{1}{3}|||{\mathit{\boldsymbol{A}}}^{n+1}_h-{\mathit{\boldsymbol{A}}}^{n}_h|||^2_{\Omega} \approx C \ \tau_{n+1}||R^{n+1}_{\rm time}||^2_{L^2(\Omega)}, \ \ \eta^{n+1}_{\rm space} = \Sigma_{K \in \mathcal{M}^{n+1}}\eta^{n+1}_K, \\ \eta^{n+1}_K&= \frac{1}{2} h_K^2 \| R^{\,n+1}_{\mathrm{ell},K}\|^2_{L^2(K)} + \frac{1}{2}\sum\limits_{e\subset\partial K} h_e \| J^{n+1}_e\|^2_{L^2(e)}, \end{aligned} \end{equation} $

(4.13)

其中$ \Omega_e $拥有公共边$ e \in \mathcal{B}^n $的元素的集合, $ C $为一个大于0的常数.

我们记$ TOL_{\rm time} $为时间相关总的容忍度, 希望满足：

$ \begin{equation} \Sigma^N_{n=1} \tau_n \eta ^n_{\rm time} + 2(\Sigma^N_{n=1}\int^{t^n}_{t^{n-1}}||\mathit{\boldsymbol{J}}_{\rm ext}-\overline{\mathit{\boldsymbol{J}}_{\rm ext}^{n+1}}||_{L^2(\Omega)}{\rm d}t)^2 \leq TOL_{\rm time}. \end{equation} $

(4.14)

达成上式的一个自然的方式是调整时间步长$ \tau_{n+1} $, 使得：

$ \begin{equation} \eta^{n+1}_{\rm time} \leq \frac{TOL_{\rm time}}{2T}, \ \frac{1}{\tau_{n+1}}\int^{t^{n+1}}_{t^{n}}||\mathit{\boldsymbol{J}}_{\rm ext}-\overline{\mathit{\boldsymbol{J}}_{\rm ext}^{n+1}}||_{L^2(\Omega)}{\rm d}t \leq \frac{1}{2T}\sqrt{TOL_{\rm time}}. \end{equation} $

(4.15)

而后我们考虑涉及空间区域的近似, 记$ TOL_{\rm space} $为空间相关忍让度, 常见的关于空间网格自适应算法在时间步$ n+1 $的停机准则为：

$ \begin{equation} \eta^{n+1}_{\rm space} \leq \frac{TOL_{\rm space}}{T}. \end{equation} $

(4.16)

此停机准则适用于网格细化. 但在自适应算法中, 我们同样需要考虑对网格进行粗化, 我们记$ TOL_{\rm coarse} $为粗化相关忍让度, 我们还需要定义粗化误差指标:

$ \begin{equation} \eta^{n+1}_{\rm coarse} = \frac{1}{\tau_n} || I^{n+1}_H {\mathit{\boldsymbol{A}}}^{n+1}_h-{\mathit{\boldsymbol{A}}}^{n+1}_h||^2_{L^2(\Omega)} + |||I^{n+1}_H {\mathit{\boldsymbol{A}}}^{n+1}_h-{\mathit{\boldsymbol{A}}}^{n+1}_h|||^2_{\Omega}. \end{equation} $

(4.17)

其中$ I^{n+1}_H: X^{n+1}_h \to X^{n+1}_H $是一个$ H(\operatorname{curl}) $插值算子.

引理3.3.1  $ \label{lem:ada} $给定$ ({\mathit{\boldsymbol{A}}}^{n-1}_h,V^{n-1}_h) \in X_h $, $ \tau_n>0 $. 令$ {\mathit{\boldsymbol{U}}}^n_*:=({\mathit{\boldsymbol{A}}}^{n}_*,V^{n}_*) $为时间离散后连续问题的解. 记$ \mathcal{M}^n_H $为将$ \mathcal{M}^n $粗化后的网格. 令$ {\mathit{\boldsymbol{U}}}^n_h:=({\mathit{\boldsymbol{A}}}^{n}_h,V^{n}_h) $, $ {\mathit{\boldsymbol{U}}}^n_H := ({\mathit{\boldsymbol{A}}}^{n}_H,V^{n}_H) $分别为问题(4.10)在网格$ \mathcal{M}^n_H $与$ \mathcal{M}^n $上的解. 则有如下的误差估计：

$ \begin{equation} ||{\mathit{\boldsymbol{U}}}^n_*-{\mathit{\boldsymbol{U}}}^n_H||^2_{\tau_n,\Omega} \leq ||{\mathit{\boldsymbol{U}}}^n_*-{\mathit{\boldsymbol{U}}}^n_h||^2_{\tau_n,\Omega} + ||{\mathit{\boldsymbol{U}}}^n_h-I^n_H{\mathit{\boldsymbol{U}}}^n_h||^2_{\tau_n,\Omega}. \end{equation} $

(4.18)

其中$ ||\boldsymbol \varphi||_{\tau_n,\Omega} := a^{n}_h(\boldsymbol \varphi,\boldsymbol \varphi) $, 对任意的$ \boldsymbol \varphi\in X $.

该引理在[16][Theorem 7.15]中有类似证明, 我们就不多做赘述.

结合此上我们定义的误差指标, 我们可以得到Maxwell方程的AVA格式关于时间与空间近似的自适应算法：

Algorithm 2  Maxwell方程AVA形式自适应算法
1:  给定容忍度$ TOL_{\rm time} $, $ TOL_{\rm space} $以及$ TOL_{\rm coarse} $与参数$ \delta_1\in(0,1) $, $ \delta_2 >1 $, 以及$ \theta_{\rm time} \in (0,1) $. 给定来自此前时间步$ n-1 $, 网格$ \mathcal{M}^{n-1} $, 时间迭代步长$ \tau_{n-1} $所得$ U^{n-1}_h $：
2:  $ \mathcal{M}^n := \mathcal{M}^{n-1} $, $ \tau_n:=\tau_{n-1}, \ t^n : = t^{n-1}+\tau_n $
基于向后Euler方法, 将数据$ {\mathit{\boldsymbol{U}}}^{n-1}_h $代入至(4.10)求解得到$ {\mathit{\boldsymbol{U}}}^n_h $并计算出$ \mathcal{M}^n $上的误差估计以及对应误差指标.
3:  $ {\rm while} $(4.15)不满足
$ \tau_n:=\delta_1 \tau_n $, $ t^n:=t^{n-1}+\tau_n $
再次将数据$ {\mathit{\boldsymbol{U}}}^{n-1}_h $代入至(4.10)求解得到$ {\mathit{\boldsymbol{U}}}^n_h $并计算出$ \mathcal{M}^n $上的误差估计以及对应误差指标.
end while
4:  $ {\rm while} $ $ \eta^n_{\rm space}>TOL_{\rm space}/T $
细化网格$ \mathcal{M}^n $, 得到修正过的网格, 同样命名为$ \mathcal{M}^n $
将数据$ {\mathit{\boldsymbol{U}}}^{n-1}_h $代入至(4.10), 在新网格求解得到$ {\mathit{\boldsymbol{U}}}^n_h $并计算出误差估计以及对应误差指标.
while(4.15)不满足
重复步骤3的操作.
end while
end while
5:  if $ \eta^n_{\rm coarse} \leq \frac{TOL_{\rm coarse}}{T} $
then
粗化网格$ \mathcal{M}^n $, 得到修正过的网格, 同样命名为$ \mathcal{M}^n $, 将数据$ {\mathit{\boldsymbol{U}}}^{n-1}_h $代入至(4.10), 在新网格求解得到$ {\mathit{\boldsymbol{U}}}^n_h $.
end if
6:  $ {\rm if} $ $ \eta^n_{\rm time}\leq \theta_{\rm time}\frac{TOL_{\rm time}}{2T},\ \frac{1}{\tau_n}\int^{t^n}_{t^{n-1}}||\mathit{\boldsymbol{J}}_{\rm ext}-\overline{\mathit{\boldsymbol{J}}_{\rm ext}^{n+1}}||_{L^2(\Omega)} {\rm d}t \leq \frac{1}{2T}\sqrt{\theta_{\rm time}TOL_{\rm time}} $
then
$ \tau_n := \delta_2 \tau_n $
end if

在实际计算中, 采用向后Euler方法时, 参数通常选择为$ \delta_1 = 0.5, \ \delta_2 = 2 $以及$ \theta_{\rm time} = 0.5 $. 此外, 我们会令粗化容让度$ TOL_{\rm coarse} $远小于空间容让度$ TOL_{\rm space} $, 并且我们会使粗化容忍度与时间容忍度$ TOL_{\rm time} $的比值也保持在较小水平.

在算法中, 第一二步为代入数据, 计算出对于误差估计与误差指标. 第三步判断整体时间容忍度是否满足, 不满足则缩小时间步长再次计算. 第四步判断空间容忍度是否满足, 用于局部网格细化. 在空间容忍度不满足时, 根据局部误差指标将网格分类, 将不满足局部误差容忍的局部网格进行加细, 加细之后再次计算, 加细之后还需考虑网格加细后对时间容忍度的影响. 第五步判断粗化容忍度是否满足, 用于去掉冗余的细网格, 如果满足, 则对网格全体进行粗化, 引理3.3.1确保在满足粗化容忍度的情况下只需要执行一次网格粗化, 并且不再需要检查粗化后问题的解是否满足之前的其余容忍度. 第六步判断时间容忍度是否充分满足, 若满足, 则减少按比例增大时间步长, 减少冗余. 自此完成一个时间迭代步的自适应算法.

采用Picard线性化与隐式Euler方法结合, 我们能够得到时间步为$ n+1 $时热-流场方程的离散格式：

$ \begin{equation} \begin{aligned} &\rho\bigg(\frac{\mathbf{u}^{n+1}_h}{\tilde\tau_n},{\mathit{\boldsymbol{v}}}_h\bigg)_{\Omega_h^{\rm oil}} + \rho((\mathbf{u}^n_h\cdot \nabla)\mathbf{u}_h^{n+1},{\mathit{\boldsymbol{v}}}_h)_{\Omega_h^{\rm oil}} + \mu (\nabla \mathbf{u}^{n+1}_h,\nabla v_h)_{\Omega_h^{\rm oil}} - (p^{n+1}_h,\nabla \cdot {\mathit{\boldsymbol{v}}}_h)_{\Omega_h^{\rm oil}} \\ =&\rho g \beta(T^{n+1}_h-T_{\rm ref},{\mathit{\boldsymbol{v}}}_h \cdot {\mathit{\boldsymbol{e}}}_g)_{\Omega_h^{\rm oil}} + \rho\bigg(\frac{\mathbf{u}_h^n}{\tilde\tau_n},{\mathit{\boldsymbol{v}}}_h\bigg)_{\Omega_h^{\rm oil}}, \\ &(\nabla \cdot \mathbf{u}_h^{n+1},q_h)_{\Omega_h^{\rm oil}} = 0, \\ &\rho c_p\bigg(\frac{T_h^{n+1}}{\tilde\tau_n} ,\theta_h\bigg)_{\Omega_h}+\rho c_p(\mathbf{u}_h^{n}\cdot \nabla T_h^{n+1},\theta_h)_{\Omega_h^{\rm oil}}+(k\nabla T^{n+1}_h,\nabla \theta_h)_{\Omega_h} +\langle(4h(T^n_h)^3+\lambda) T^{n+1}_h,\theta_h\rangle_{\partial\Omega_h}\\ =&\bigg(\bar{Q}^{n+1}+\rho c_p\frac{T^n_h}{\tilde\tau_n},\theta_h\bigg)_{\Omega_h}+\langle3h(T^{n}_h)^4+hT^4_{\rm ref}+\lambda T_{\rm ref },\theta_h \rangle_{\partial \Omega_h}, \end{aligned} \end{equation} $

(4.19)

其中$ \bar{Q}^{n+1} = \frac{1}{\tilde\tau_{n+1}}\int^{t_{n+1}}_{t_n} Q \ {\rm}d t $

热-流体耦合方程可以采用类似自适应方法进行计算, 我们在此不多作赘述, 主要叙述对应误差指标.

对于油区单元$ K \subset \Omega_{\rm oil} $, 我们可以定义动量残差$ R^{n+1}_{\rm mom,K} $:

$ \begin{equation} \begin{aligned} R^{n+1}_{\rm mom,K} = \rho (\frac{\mathbf{u}^{n+1}_h-\mathbf{u}^n_h}{\tilde\tau_{n+1}}+(\mathbf{u}^n_h\cdot \nabla)\mathbf{u}^{n+1}_h) + \nabla p^{n+1}_h - \mu \Delta \mathbf{u}^{n+1}_h -\rho g \beta(T^{n+1}_h-T_{\rm ref})e_g. \end{aligned} \end{equation} $

(4.20)

方程残差$ R^{n+1}_{\rm div,K}:=\nabla \cdot \mathbf{u}^{n+1}_h $, 能量残差$ R^{n+1}_{T,K} $:

$ \begin{equation} R^{n+1}_{T,K} = \rho c_p (\frac{T^{n+1}_h-T^n_h}{\tilde\tau_{n+1}}+\mathbf{u}^{n}_h \cdot \nabla T^{n+1}_h)-\nabla \cdot (k \nabla T^{n+1}_h)-\bar{Q}^{n+1}. \end{equation} $

(4.21)

同样地, 我们定义跳跃残差, 我们选择公共面$ e = \partial K_1 \cap \partial K_2 $, 法向$ {\mathit{\boldsymbol{n}}}_e $从$ K_2 \to K_1 $.

$ \begin{equation} \begin{aligned} J^{n+1}_{{\rm mom},e}&:=[[\mu(\nabla\mathbf{u}^{n+1}_h){\mathit{\boldsymbol{n}}}_e-p^{n+1}_h {\mathit{\boldsymbol{n}}}_e]]_e, \quad J^{n+1}_{T,e} &:= [[k\nabla T^{n+1}_h \cdot {\mathit{\boldsymbol{n}}}_e]]_e. \end{aligned} \end{equation} $

(4.22)

最后我们定义边界残差$ R^{n+1}_{\rm bc} $：

$ \begin{equation} R^{n+1}_{\rm bc} = -k \nabla T^{n+1}_h \cdot n - (h(T^{n+1}_h)^4-T^4_{\rm ref}) - \lambda (T^{n+1}_h-T_{\rm ref}). \end{equation} $

(4.23)

现在我们可以定义误差指标. 首先我们定义局部误差指标$ \eta ^{n+1}_k $:

$ \begin{equation} \begin{aligned} (\eta ^{n+1}_K)^2 =& h^2_K (||R^{n+1}_{{\rm mom},K}||^2_{L^2(K)} + ||R^{n+1}_{{\rm div},K})||^2_{L^2(K)}+||R^{n+1}_{T,K})||^2_{L^2(K)} \\ &+ \Sigma _{e \subset \partial K} h_e (||J^{n+1}_{{\rm mom},e})||^2_{L^2(e)}+||J^{n+1}_{T,e})||^2_{L^2(e)} + h_K ||R^{n+1}_{\rm bc}||^2_{L^2(\partial \Omega \cap \partial K)}. \end{aligned} \end{equation} $

(4.24)

类似地, 空间误差指标即为$ \eta^{n+1}_{\rm space}:=\Sigma_{K \in \mathcal{M}^{n+1}}\eta^{n+1}_k $, 我们记时间误差指标为$ \eta ^{n+1}_{\rm time} $, 满足：$ \eta ^{n+1}_{\rm time} = ||\mathbf{u}^{n+1}_h - \mathbf{u}^n_h||^2_{L^2(\Omega)} + ||T^{n+1}_h-T^n_h||^2_{L^2(\Omega)}. $最后我们定义粗化误差指标: $ \eta^{n+1}_{\rm coarse}:= \frac{1}{\tilde\tau_{n+1}} || I^{n+1}_H {\mathit{\boldsymbol{U}}}^{n+1}_h-{\mathit{\boldsymbol{U}}}^{n+1}_h||^2_{L^2(\Omega)} + || I^{n+1}_H {\mathit{\boldsymbol{U}}}^{n+1}_h-{\mathit{\boldsymbol{U}}}^{n+1}_h||^2_{H^1(\Omega)}. $

定义误差指标后, 热-流体场上单步自适应算法与此前的自适应算法几乎相同, 只是对应计算方程与误差指标有所更新, 因此我们暂不详细叙述.

结合此前我们提出的Maxwell方程AVA格式的自适应算法以及热-流体场耦合的自适应算法, 首先我们进行定义, 电磁场时间步长为$ \tau_n $, 对应时间为$ t^n = \Sigma^n_i \tau_i $, 热-流体场时间步长为$ \tilde\tau_n $, 对应时间为$ t_n = \Sigma _i^n \tilde\tau_n $. 此处我们给出在热场时间步$ n $处的多物理场耦合的自适应算法：

Algorithm 3  多物理场耦合自适应算法
1:  给定容忍度$ TOL_{\rm time} $, $ TOL_{\rm space} $以及$ TOL_{\rm coarse} $与参数$ \delta_1\in(0,1) $, $ \delta_2 >1 $, 要求$ \delta_1 \delta_2 = 1 $. 给定初始网格$ \mathcal{M}^0 $以及初始$ {\mathit{\boldsymbol{U}}}^0_{1,h} = ({\mathit{\boldsymbol{A}}}^0,V^0) $, $ {\mathit{\boldsymbol{U}}}^0_{2,h} = (\mathbf{u}^0,p^0,T_{\rm ref}) $, 令$ n=1,\ t_0 = 0 $
2:  基于给定$ T = T_{\rm ref} $, 更新电磁场方程(4.10)中系数$ \sigma(T),\ \nu(T) $. 给定热-流体时间步长为$ \tilde\tau_1 $, 选取电磁场方程初始时间步长$ \tau_1 $, 需要满足$ \tilde\tau_1 = c_1\delta_2^{c_2} \tau_1 $, 其中$ c_1,c_2 $均为正整数. 在时间区间$ [0,\tilde\tau_1] $应用Maxwell方程AVA格式自适应算法2, 我们记求解所得为$ {\mathit{\boldsymbol{U}}}^1_{1,h} = ({\mathit{\boldsymbol{A}}}^1_h, V^1_h) $, 自适应算法后的网格为$ \mathcal{M}^1 $, 以及最后的时间步长序列$ \{ \tau^1_n\} $.
将$ ({\mathit{\boldsymbol{A}}}^1_h, V^1_h) $代入, 求解出对应$ {\mathit{\boldsymbol{B}}}^1_h,\ E^1_h $, 将其代入[24]所得损耗计算公式, 得到瞬时发热功率$ Q^1 $, 对应计算$ \bar{Q}^1 $, 而后将其代入至热-流场耦合方程(4.19), 以时间步长$ \tilde\tau_1 $, 网格$ \mathcal{M}^1 $为初始网格, 进行自适应算法. 倘若需要对时间进行细化, 则记细化后时间为$ \tilde\tau_1' $, 以时间区间$ (0,\tilde\tau_1'), (\tilde\tau'_1, \tilde\tau_1) $进行自适应算法求解. 在执行自适应算法时, 若需要粗化, 则执行算法2粗化判断, 若条件满足, 则粗化. 最后记求得解$ {\mathit{\boldsymbol{U}}}^1_{2,h} = (\mathbf{u}^1_h,p^1_h,T^1_{h}) $, 所得网格记为$ \mathcal{M}^{1} $,
3:  第$ n $次时间迭代：$ \mathcal{M}^n := \mathcal{M}^{n-1} $, $ \tau_n:=\tau_{n-1}, \tilde\tau_n : = \tilde\tau_{n-1} $, 在热时间区间$ (t_{n-1},t_{n-1}+\tilde\tau_n) $, 基于此前所求$ {\mathit{\boldsymbol{U}}}^{n-1}_{1,h} $, $ {\mathit{\boldsymbol{U}}}^{n-1}_{2,h} $, 将$ T^{n-1}_h $代入电磁场方程(4.10), 更新电磁场方程中系数$ \sigma , \nu $. 选取电磁场方程初始时间步长$ \tau_n $, 需要满足$ \tilde\tau_n = c_1\delta_2^{c_2} \tau_n $, 以初始时间为$ t_0 = t_{n-1} $, 终止时间为$ t_{n-1} + \tilde\tau_n $进行Maxwell方程AVA格式自适应算法2, 我们记求解所得为$ {\mathit{\boldsymbol{U}}}^n_{1,h} = ({\mathit{\boldsymbol{A}}}^n_h, V^n_h) $, 自适应算法后的网格为$ \mathcal{M}^n $, 以及最后的时间步长序列$ \{ \tau^n_n\} $.
将$ ({\mathit{\boldsymbol{A}}}^n_h, V^n_h) $代入, 求解出对应$ {\mathit{\boldsymbol{B}}}^n_h,\ E^n_h $, 将其代入[24]所得损耗计算公式, 得到瞬时发热功率$ Q^n $, 对应计算$ \bar{Q}^n $, 而后将其代入至热-流场耦合方程(4.19), 以时间步长$ \tilde\tau_n $, 网格$ \mathcal{M}^n $为初始网格, 进行自适应算法. 倘若需要对时间进行细化, 则记细化后时间为$ \tilde\tau_n' $, 以时间区间$ (t_{n-1},(t_{n-1}+\tilde\tau_n'), (t_{n-1}+\tilde\tau_n', t_{n-1}+\tilde\tau_n) $进行两次自适应算法求解. 在执行自适应算法时, 若需要粗化, 则执行算法2粗化判断, 若条件满足, 则粗化. 最后记求得解$ {\mathit{\boldsymbol{U}}}^n_{2,h} = (\mathbf{u}^n_h,p^n_h,T^n_{h}) $, 所得网格为$ \mathcal{M}^{n} $.
4:  令$ n:=n+1 $, REPEAT

在我们的多物理场耦合自适应算法中, 核心是利用Picard线性化对方程进行解耦与减少计算量, 我们可以首先对整个方程组做Picard线性化, 我们可以发现在电磁场方程中的非线性项只为$ \sigma, \nu $, 因此我们选择上个时间步所计算的温度$ T $, 即能得到对应的$ \sigma, \nu $, 能够将其代入至电磁场方程中, 而后电磁场在完成线性化的同时, 也完成了对方程组的解耦, 使得电磁场方程得以单独摘出, 将方程组解耦成电磁场方程与热-流场耦合方程. 对于热-流场耦合方程而言, 我们在采取Picard线性化时, 除了本身明显的非线性项外, $ Q $是一个与$ {\mathit{\boldsymbol{B}}},\ {\mathit{\boldsymbol{E}}} $相关的非线性项, 因此我们可以选择在同一时间电磁场的计算结果将其代入至$ Q $中, 将其化简. 通过这样的方式我们不仅能够实现对方程的解耦, 也能节约计算资源提高计算效率.

自适应有限元算法的核心是用尽可能少的网格单元, 在给定误差容忍下得到尽量高的精度, 因此我们需要根据局部误差指示, 对误差大的局部网格进行加细, 增加计算精度. 在容忍满足的前提下, 我们需要考虑整体网格是否冗余, 整体网格是否需要退回, 由此节约计算资源, 从而提高整体的计算效率. 在我们的算法中我们不仅考虑了对空间网格的自适应调整, 同样也考虑了对时间步长的自适应调整, 缩小时间步长以达到精度要求, 增加时间步长以此减少计算次数, 减少计算冗余, 同样提高计算效率.

在处理油浸式变压器多物理场耦合模型时, 通过对单相变压器几何模型分析, 我们将变压器划分为铁芯、绕组以及油区三部分区域, 由于计算区域不是一般的凸区域, 解在区域的拐角处会有比较大的奇异性, 导致数值解在区域拐角处的误差比较大, 我们需要对区域拐角局部进行加密以达到计算精度要求. 其次是由于变压器三部分区域的材料不同, 磁导率、电导率等物理参数在不同介质中的数量级差别较大, 这同样会导致数值解在区域界面处的误差比较大, 因此我们还需要在界面处局部加密来保证解的收敛性. 在用有限元方法进行求解, 需要考虑局部网格加密时, 我们自然会选择自适应方法来对问题进行处理.

在本文中, 我们基于对油浸式变压器进行建模, 得到对应的多物理场耦合模型. 针对复杂的耦合模型以及几何条件, 我们利用Picard线性化将电磁-热-流场三场耦合方程简化为电磁场与热-流场耦合方程的迭代求解. 针对电磁场我们采用了AVA形式来保证求解的稳定性与精确度, 引入了Lagrange乘子或是惩罚项来消除常数自由度, 针对热-流场耦合方程, 我们采用向后Euler方法来处理时间导数. 此外, 由于变压器内部区域物质性不连续, 并且几何条件复杂, 导致PDE求解不稳定, 我们采用了自适应有限元算法对问题进行求解, 能够实现对网格的局部细化、对时间步长的具体协调、使得保证精确度的情况下, 尽可能节约计算资源.