数学物理学报, 2022, 42(4): 1238-1255 doi:

论文

非线性感应加热问题的全离散有限元方法

黄媛1, 支越3, 康彤1,2, 王然,1,2, 张红4

1 中国传媒大学数据科学与智能媒体学院 北京 100024

2 中国科学院大学计算地球动力学重点实验室 北京 100049

3 北京市育才学校通州分校 北京 101101

4 北京汇文中学朝阳垂杨柳分校 北京 100021

Fully Discrete Finite Element Scheme for a Nonlinear Induction Heating Problem

Huang Yuan1, Zhi Yue3, Kang Tong1,2, Wang Ran,1,2, Zhang Hong4

1 School of Data Science and Intelligent Media, Communication University of China, Beijing 100024

2 Key Laboratory of Computational Geodynamics, University of Chinese Academy of Sciences, Beijing 100049

3 Tongzhou Branch of Beijing Yucai School, Beijing 101101

4 Chaoyang Chuiyangliu Branch of Beijing Huiwen Middle School, Beijing 100021

通讯作者: 王然, E-mail: wangr@ucas.ac.cn

收稿日期: 2021-03-30  

基金资助: 国家重点研发计划项目.  2020YFA0713401
国家自然科学基金项目.  U2039207
国家自然科学基金项目.  42074108
国家自然科学基金项目.  41904067
中国传媒大学中央高校基本科研业务费专项资金

Received: 2021-03-30  

Fund supported: the National Key Research and Development Program of China.  2020YFA0713401
the National Science Foundation of China.  U2039207
the National Science Foundation of China.  42074108
the National Science Foundation of China.  41904067
the Fundamental Research Funds for the Central Universities

Abstract

It studies an induction heating model described by Maxwell's equations coupled with a heat equation. In the induction heating model, it assumes a nonlinear relation between the magnetic field and the magnetic induction field, and the electric conductivity is temperature dependent. It presents a fully discrete H-based finite element scheme in time and space and discusses its solvability. Moreover, it proves the fully discrete solution converges to a weak solution of the continuous problem. Finally, theoretical results are supported by some numerical experiments.

Keywords: Induction heating ; Maxwell's equations ; Nonlinear ; Finite elements ; Convergence

PDF (1749KB) 元数据 多维度评价 相关文章 导出 EndNote| Ris| Bibtex  收藏本文

本文引用格式

黄媛, 支越, 康彤, 王然, 张红. 非线性感应加热问题的全离散有限元方法. 数学物理学报[J], 2022, 42(4): 1238-1255 doi:

Huang Yuan, Zhi Yue, Kang Tong, Wang Ran, Zhang Hong. Fully Discrete Finite Element Scheme for a Nonlinear Induction Heating Problem. Acta Mathematica Scientia[J], 2022, 42(4): 1238-1255 doi:

1 引言

电磁感应加热是一种非接触的加热导体材料的方法, 该方法常用于工业应用中, 例如金属硬化和锻造预热. 感应加热系统的基本组成包括感应线圈(电感器), 交流电源和工件(见图 1). 交变电流通过电感器会产生交变磁场, 该磁场使工件产生涡流, 从而耗散涡流的能量并产生焦耳热. 关于感应加热系统的研究, 通常依赖于一系列复杂, 昂贵且漫长的实验, 因此, 数学理论分析和数值模拟可以在设计感应加热过程中发挥重要作用.

图 1

图 1   圆盘工件和感应线圈(电感器)


感应加热模型由Maxwell方程和热传导方程组成. 对于线性电磁材料, 文献[13, 17]提出了各种数值计算方案, 文献[4, 7, 9, 10, 1921]研究了耦合问题的适定性并给出了理论结果. 据我们所知, 很少有同时考虑非线性关系$ \mathit{\boldsymbol{B}} = \mathit{\boldsymbol{B}} ( \mathit{\boldsymbol{H}} )$和电导率依赖于温度函数的研究. 例如, 文献[16]研究了导体中基于非线性$ \mathit{\boldsymbol{H}} $的耦合系统, 并讨论了其可解性. 文献[6, 11]研究了区域中包含导体和非导体部分的非线性模型, 给出了模型的矢量—标量势公式. 这些文献仅涉及对非线性感应加热问题时间半离散格式的研究. 目前为止, 很少有文献研究全离散有限元方法.

本文的目的是对非线性感应加热问题提出全离散$ \mathit{\boldsymbol{H}} $有限元方法, 讨论全离散格式解的收敛性. 文章组织结构如下: 第二章将建立非线性感应加热模型; 第三章给出基于$ \mathit{\boldsymbol{H}} $形式的时间半离散格式和一些收敛结果; 第四章提出全离散有限元解耦格式, 证明全离散格式解的存在唯一性, 给出稳定性估计; 第五章将讨论全离散格式解的收敛性. 最后, 给出数值实验来支持理论结果.

2 非线性感应加热模型

下面考虑一个简化的感应加热模型. $ \Omega\subset{\Bbb R}^3 $是一个带有连续边界$ \partial\Omega $的区域, 被电磁材料所填充. 电磁场方程由涡流方程描述

$ \begin{equation} \left\{ \begin{array}{ll} & \partial_t\mathit{\boldsymbol{B}} + \nabla \times\mathit{\boldsymbol{E}} =0, \\ & \nabla \times\mathit{\boldsymbol{H}} = \mathit{\boldsymbol{J}}+\mathit{\boldsymbol{J}}_s, \end{array} \right. \end{equation} $

其中$ \mathit{\boldsymbol{E}} $是电场, $ \mathit{\boldsymbol{H}} $是磁场, $ \mathit{\boldsymbol{B}} $是磁通密度, $ \mathit{\boldsymbol{J}}_s $$ \Omega $外部电感器中的源电流密度. 假设$ \mathit{\boldsymbol{J}}_s $是无散度的, 那么一个交变磁场$ \mathit{\boldsymbol{H}}_s $使得

使用Biot-Savart定律可以直接计算出$ \mathit{\boldsymbol{H}} _s $如下

这里的$ \Omega $是一个足够大的关于$ \mathit{\boldsymbol{y}} $的积分区域. 假设集合$ \{ \mathit{\boldsymbol{y}}: \mathit{\boldsymbol{J}}_s \not\equiv {\bf{0}} \} $到边界$ \partial \Omega $的距离足够远, 使得$ \mathit{\boldsymbol{H}}_s $$ \partial \Omega $上近似满足$ \mathit{\boldsymbol{H}}_s \times \mathit{\boldsymbol{n}} = {\bf{0}} $. 交变磁场在$ \Omega $中产生涡流, 即$ \mathit{\boldsymbol{J}}=\sigma\mathit{\boldsymbol{E}} $, 其中电导率$ \sigma = \sigma(u) $依赖于温度函数$ u $. $ \sigma $是严格正的且有界, 即存在正的常数$ \sigma_* $$ \sigma^* $, 使得

$ \begin{equation} 0<\sigma_*\leq \sigma(s)\leq \sigma^*<\infty, \; \; s>0. \end{equation} $

引入$ \sigma $的反函数$ \gamma $, 同样$ \gamma $也是有界的, 即$ 0 <\gamma_*\leq \gamma\leq \gamma^*<\infty $.

研究$ \Omega $$ \mathit{\boldsymbol{H}} $$ \mathit{\boldsymbol{B}} $之间呈非线性关系的这类模型是至关重要的. 此外, 铁磁体的磁化曲线具有单调性. $ \mathit{\boldsymbol{B}} $$ \mathit{\boldsymbol{H}} $之间的非线性关系为$ \mathit{\boldsymbol{B}} = \mathit{\boldsymbol{B}}(\mathit{\boldsymbol{H}}) $. 假设矢量场$ \mathit{\boldsymbol{B}} $为有势场, Lipschitz连续的和严格单调的. 仿照文献[18], 定义$ \mathit{\boldsymbol{B}} $的势为$ \Phi_{\mathit{\boldsymbol{B}}} $, 即$ \mbox{grad}\; \Phi_{{\mathit{\boldsymbol{B}}}}=\mathit{\boldsymbol{B}} $, 满足下列条件

$ \begin{equation} \begin{array}{ll} (\mathit{\boldsymbol{B}}(\mathit{\boldsymbol{x}})-\mathit{\boldsymbol{B}}(\mathit{\boldsymbol{y}}))\cdot (\mathit{\boldsymbol{x}}-\mathit{\boldsymbol{y}})\geq b_*|\mathit{\boldsymbol{x}}-\mathit{\boldsymbol{y}}|^2, \; \; b_*>0\; \; \; \; \; &\forall\mathit{\boldsymbol{x}}, \mathit{\boldsymbol{y}}\in {\Bbb R}^3, \\ \mathit{\boldsymbol{B}}{({\bf{0}})}={\bf{0}}, \\ \mathit{\boldsymbol{B}}^{-1}(\mathit{\boldsymbol{x}})\cdot(\mathit{\boldsymbol{x}}-\mathit{\boldsymbol{y}})\geq \Phi_{\mathit{\boldsymbol{B}}^{-1}}(\mathit{\boldsymbol{x}})-\Phi_{\mathit{\boldsymbol{B}}^{-1}}(\mathit{\boldsymbol{y}})&\forall\mathit{\boldsymbol{x}}, \mathit{\boldsymbol{y}}\in{\Bbb R}^3, \\ \Phi_{{\mathit{\boldsymbol{B}}}^{-1}}(\mathit{\boldsymbol{B}}(\mathit{\boldsymbol{x}}))\geq C_0|\mathit{\boldsymbol{x}}|^2 &\forall\mathit{\boldsymbol{x}}\in{\Bbb R}^3. \end{array} \end{equation} $

$ \mathit{\boldsymbol{B}} $的严格单调性可知, $ \mathit{\boldsymbol{B}} $为可逆的, $ \mathit{\boldsymbol{B}}^{-1} $是Lipschitz连续的. 根据文献[18, 定理5.1], $ \mathit{\boldsymbol{B}} $的势$ \Phi_{\mathit{\boldsymbol{B}}} $是严格凸的. 由文献[18, 定理8.4]得到

$ \begin{equation} \begin{array}{ll} \mathit{\boldsymbol{B}}^{-1}(\mathit{\boldsymbol{x}})\cdot(\mathit{\boldsymbol{x}}-\mathit{\boldsymbol{y}})\geq \Phi_{\mathit{\boldsymbol{B}}^{-1}}(\mathit{\boldsymbol{x}})-\Phi_{\mathit{\boldsymbol{B}}^{-1}}(\mathit{\boldsymbol{y}})\hskip 1cm \forall\mathit{\boldsymbol{x}}, \mathit{\boldsymbol{y}}\in{\Bbb R}^3, \end{array} \end{equation} $

$ \begin{equation} \frac{\rm d}{{\rm d}t}\Phi_{{\mathit{\boldsymbol{B}}}^{-1}}\big(\mathit{\boldsymbol{B}}(\mathit{\boldsymbol{x}})\big)={\mathit{\boldsymbol{B}}}^{-1}\big(\mathit{\boldsymbol{B}}(\mathit{\boldsymbol{x}})\big)\cdot \frac{{\rm d}\mathit{\boldsymbol{B}}(\mathit{\boldsymbol{x}})}{{\rm d}t}=\mathit{\boldsymbol{x}}\cdot \frac{{\rm d}\mathit{\boldsymbol{B}}(\mathit{\boldsymbol{x}})}{{\rm d}t}. \end{equation} $

考虑下述磁场的初始边值问题

$ \begin{equation} \left\{ \begin{array}{ll} \partial_t\mathit{\boldsymbol{B}}(\mathit{\boldsymbol{H}}'+\mathit{\boldsymbol{H}}_s) + \nabla \times\big(\gamma(u)\nabla \times\mathit{\boldsymbol{H}}'\big) = {\bf{0}} & \; (\mathit{\boldsymbol{x}}, t)\in \Omega\times (0, T), \\ \mathit{\boldsymbol{H}}'\times\mathit{\boldsymbol{n}}={\bf{0}} & \; (\mathit{\boldsymbol{x}}, t)\in \partial\Omega\times (0, T), \\ \mathit{\boldsymbol{H}}'(0) = \mathit{\boldsymbol{H}}_0 & \; \mathit{\boldsymbol{x}}\in \Omega, \; t=0, \end{array} \right. \end{equation} $

其中$ {\mathit{\boldsymbol{n}}} $$ \partial\Omega $上的单位外法线向量. 初值$ \mathit{\boldsymbol{H}}_0 $$ \partial \Omega $上满足$ \mathit{\boldsymbol{H}}_0 \times \mathit{\boldsymbol{n}} = {\bf{0}} $.

本文主要研究如何求解$ \mathit{\boldsymbol{H}}' $. 为了简便起见, 下述内容依然使用$ \mathit{\boldsymbol{H}} $代替$ \mathit{\boldsymbol{H}}' $.

由于$ \mathit{\boldsymbol{B}} $的散度为零, 设

$ \Omega $中电流产生的局部焦耳热为

为了处理热传导方程中的焦耳热项, 对其引入截断函数[16]

其中$ r $是正的常数. 温度函数$ u $的变化由下面带有初始条件和边界条件的热传导方程表示

$ \begin{equation} \left\{ \begin{array}{ll} \partial_t u-\nabla \cdot(\lambda\nabla u)={\cal R}_r\big(\gamma(u)|\nabla \times\mathit{\boldsymbol{H}}|^2\big) & \; (\mathit{\boldsymbol{x}}, t)\in \Omega\times (0, T), \\ { } -\lambda\frac{\partial u}{\partial\mathit{\boldsymbol{n}}}=0 & \; (\mathit{\boldsymbol{x}}, t)\in \partial\Omega\times (0, T), \\ u(0)=u_0& \; \; \mathit{\boldsymbol{x}}\in \Omega, \; t=0, \end{array} \right. \end{equation} $

其中导热系数$ \lambda(\mathit{\boldsymbol{x}}, t) $是严格正的, 有界的$ (0<\lambda_*\leq \lambda\leq \lambda^*<\infty) $, 满足

$ \begin{equation} |\lambda(\mathit{\boldsymbol{x}}, t_2)-\lambda(\mathit{\boldsymbol{x}}, t_1)|\leq C_\lambda|t_2-t_1|\hskip 0.5cm \forall\mathit{\boldsymbol{x}}\in\Omega, \; \forall t_1, t_2\in[0, T]. \end{equation} $

通过(2.6)式中的$ \gamma(u) $项和(2.7)式中的焦耳热项, 磁场方程和热传导方程被耦合在一起.

3 时间离散化变分公式

首先, 给出一些贯穿全文的概念和符号. 定义$ L^2(\Omega) $空间的內积及范数

定义$ H^m(\Omega) := \{v\in L^2(\Omega): D^\xi v \in L^2(\Omega), |\xi|\leq m\} $空间的范数

其中$ \xi $表示非负三重指数. 使用黑体表示向量值空间, 如$ \mathit{\boldsymbol{L}}^2(\Omega) := \big(L^2(\Omega)\big)^3 $. 定义Hilbert空间

相应的內积和范数为

$ \mathit{\boldsymbol{n}} $表示$ \partial\Omega $上的单位外法线向量.

应用上述定义, (2.6)和(2.7)式的变分形式为

$ \begin{equation} \begin{array}{ll} \big(\partial_t\mathit{\boldsymbol{B}}(\mathit{\boldsymbol{H}}+\mathit{\boldsymbol{H}}_s), \mathit{\boldsymbol{Q}}\big)+\big(\gamma(u)\nabla \times\mathit{\boldsymbol{H}}, \nabla \times\mathit{\boldsymbol{Q}}\big)=0, \quad \forall\mathit{\boldsymbol{Q}}\in\mathit{\boldsymbol{H}}_0(\mathit{\boldsymbol{curl}}, \Omega) \end{array} \end{equation} $

$ \begin{equation} \big(\partial_t u, v\big)+\big(\lambda\nabla u, \nabla v\big)=\Big({\cal R}_r\big(\gamma(u) |\nabla \times\mathit{\boldsymbol{H}}|^2\big), v\Big), \quad \forall v\in H^1(\Omega). \end{equation} $

下面, 基于向后欧拉法, 给出(3.1)–(3.2)式的时间离散格式. 设$ n $为正整数, 将时间$ [0, T] $均匀剖分成$ n $份, 时间步长$ \tau = T/n $, 则$ \{t_i = i\tau:\; i = 0, \cdots, n\} $. 对于任意函数$ z $, 有

时间离散格式为: 对于$ i=1, \cdots, n $,

$ \begin{equation} \begin{array}{ll} \big(\delta\mathit{\boldsymbol{B}}(\mathit{\boldsymbol{h}}_i+\mathit{\boldsymbol{H}}_{s, i}), \mathit{\boldsymbol{Q}}\big)+\big(\gamma(u_{i-1})\nabla \times\mathit{\boldsymbol{h}}_i, \nabla \times\mathit{\boldsymbol{Q}}\big)=0, \quad \forall\mathit{\boldsymbol{Q}}\in\mathit{\boldsymbol{H}}_0(\mathit{\boldsymbol{curl}}, \Omega), \end{array} \end{equation} $

$ \begin{equation} \big(\delta u_i, v\big)+\big(\lambda_i\nabla u_i, \nabla v\big)=\Big({\cal R}_r\big(\gamma(u_{i-1})|\nabla \times\mathit{\boldsymbol{h}}_i|^2\big), v\Big), \quad \forall v\in H^1(\Omega). \end{equation} $

使用单调算子引理[13, 18], 证明在每个时间步上, (3.3)式解的存在唯一性.

定理3.1  令(2.2), (2.3), (2.8)式均满足, 设$ \mathit{\boldsymbol{H}}_0\in \mathit{\boldsymbol{L}}^2(\Omega) $, $ u_0\in L^2(\Omega) $, $ \mathit{\boldsymbol{H}}_{s}\in L^2((0, T); $$ \mathit{\boldsymbol{L}}^2(\Omega)) $, 则对于任意的$ i=1, \cdots, n $, (3.3)和(3.4)式存在唯一的解$ \mathit{\boldsymbol{h}}_i\in\mathit{\boldsymbol{H}}_0(\mathit{\boldsymbol{curl}}, \Omega) $, $ u_i\in H^1(\Omega) $.

  定义算子$ {\cal L}_{\gamma, i}:\mathit{\boldsymbol{H}}_0(\mathit{\boldsymbol{curl}}, \Omega)\rightarrow \mathit{\boldsymbol{H}}^{-1}(\mathit{\boldsymbol{curl}}, \Omega) $满足

为了在时间步$ t_i $上得到唯一的解$ \mathit{\boldsymbol{h}}_i $, 需要满足

由于上式右端已知, 算子$ {\cal L}_{\gamma(u_{i-1}), i} $是半连续的, 严格单调的和强制的. 根据文献[18, 定理18.2], 能够证明(3.3)式解的存在性.

接下来给出唯一性. 设$ \mathit{\boldsymbol{h}}_1 $$ \mathit{\boldsymbol{h}}_2 $是两个不同的解, 有

因此, 可以得到$ \mathit{\boldsymbol{h}}_1=\mathit{\boldsymbol{h}}_2 $.

显然, 由Lax-Milgram定理可知, 线性问题(3.4)具有唯一的解.

此外, 构造时间上的分片常数插值和分片线性函数插值, 对于任意的$ i=1, \cdots, n $, 函数$ f_{i} $表示为

在文献[16]中证明了下述定理, 该定理给出了时间离散格式(3.3)–(3.4)解$ (\mathit{\boldsymbol{h}}_i, u_i) $的收敛结果, 证明了(3.1)–(3.2)式弱解的存在性.

定理3.2  令(2.2)–(2.5), (2.8)式均满足, 设$ u_0\in H^1(\Omega) $, $ \mathit{\boldsymbol{H}}_{s}\in H^1((0, T);\mathit{\boldsymbol{H}}(\mathit{\boldsymbol{curl}}, \Omega)) $, $ \gamma(s) $是全局Lipschitz连续函数, 使得

4 全离散有限元解耦格式

本节将给出(3.1)–(3.2)式的全离散有限元格式. 设$ {\cal T}^h $$ \Omega $上的正则的四面体三角剖分, 网格尺寸为$ h $. 引入Nédélec $ \mathit{\boldsymbol{H}}(\mathit{\boldsymbol{curl}}, \Omega) $—协调有限元空间[14], 定义如下

其中$ \mathit{\boldsymbol{a}}_K $$ \mathit{\boldsymbol{b}}_K $为常数向量. 已知任意函数$ \mathit{\boldsymbol{v}}^h\in \mathit{\boldsymbol{V}}_h $都是由每个元素$ K\in {\cal T}^h $的集合$ M_E(\mathit{\boldsymbol{v}}) $中的自由度唯一确定的

其中$ \tau $是边缘上的单位切向量. 定义有限元空间

其中$ P_1(K) $是线性多项式空间. 令$ \mathit{\boldsymbol{\pi}}^h:\mathit{\boldsymbol{H}}(\mathit{\boldsymbol{curl}}, \Omega)\rightarrow \mathit{\boldsymbol{V}}_h $是投影算子, 定义为

$ \begin{equation} (\nabla \times\mathit{\boldsymbol{\pi}}^h \mathit{\boldsymbol{U}}, \nabla\times\mathit{\boldsymbol{Q}}^h)+(\mathit{\boldsymbol{\pi}}^h \mathit{\boldsymbol{U}}, \mathit{\boldsymbol{Q}}^h)=(\nabla \times \mathit{\boldsymbol{U}}, \nabla\times\mathit{\boldsymbol{Q}}^h)+(\mathit{\boldsymbol{U}}, \mathit{\boldsymbol{Q}}^h)\; \; \forall \mathit{\boldsymbol{Q}}^h\in \mathit{\boldsymbol{V}}_h, \end{equation} $

$ \Pi^h:H^1(\Omega)\rightarrow W_h $是Galerkin投影算子, 有

$ \begin{equation} (\nabla \Pi^h u, \nabla v^h)+( \Pi^h u, v^h) =(\nabla u, \nabla v^h)+(u, v^h)\; \; \forall v^h\in W_h. \end{equation} $

全离散有限元格式如下: 给定$ \mathit{\boldsymbol{h}}_0^h=\mathit{\boldsymbol{\pi}}^h \mathit{\boldsymbol{H}}_0 $$ u_0^h=\Pi^h u_0 $, 对于任意的$ i=1, \cdots, n $, 求解$ \mathit{\boldsymbol{h}}_i^h\in \mathit{\boldsymbol{V}}_h^0 $$ u_i^h\in W_h $, 使得

$ \begin{equation} \big(\delta\mathit{\boldsymbol{B}}(\mathit{\boldsymbol{h}}_i^h+\mathit{\boldsymbol{\pi}}^h\mathit{\boldsymbol{H}}_{s, i}), \mathit{\boldsymbol{Q}}^h\big)+\big(\gamma(u_{i-1}^h)\nabla \times\mathit{\boldsymbol{h}}_i^h, \nabla \times\mathit{\boldsymbol{Q}}^h\big)=0, \quad \forall \mathit{\boldsymbol{Q}}^h\in \mathit{\boldsymbol{V}}_h^0, \end{equation} $

$ \begin{equation} \big(\delta u_i^h, v^h\big)+\big(\lambda_i\nabla u_i^h, \nabla v^h\big)=\Big({\cal R}_r\big(\gamma(u_{i-1}^h)|\nabla \times\mathit{\boldsymbol{h}}_i^h|^2\big), v^h\Big), \quad \forall v^h\in W_h. \end{equation} $

在下面的定理中, 将证明离散问题(4.3)的可解性. 使用Lax-Milgram定理, 容易得到(4.4)式解是存在唯一的.

定理4.1  令定理3.1的假设条件均满足. 对$ \forall i=1, \cdots, n $, 全离散格式(4.3)存在唯一的解$ \mathit{\boldsymbol{h}}_i^h\in \mathit{\boldsymbol{V}}_h^0 $.

  在每个时间步上需要求解$ \mathit{\boldsymbol{h}}^h\in \mathit{\boldsymbol{V}}_h^0 $

$ \begin{equation} \begin{array}{ll} \frac{1}{\tau}\big(\mathit{\boldsymbol{B}}(\mathit{\boldsymbol{h}}^h+\mathit{\boldsymbol{\pi}}^h\mathit{\boldsymbol{H}}_{s})-\mathit{\boldsymbol{B}}(\mathit{\boldsymbol{\pi}}^h\mathit{\boldsymbol{H}}_{s}), \mathit{\boldsymbol{Q}}^h\big)+\big(\gamma\nabla\times \mathit{\boldsymbol{h}}^h, \nabla\times\mathit{\boldsymbol{Q}}^h\big)=\big(\mathit{\boldsymbol{f}}, \mathit{\boldsymbol{Q}}^h\big), \end{array} \end{equation} $

$ \forall\mathit{\boldsymbol{Q}}^h\in \mathit{\boldsymbol{V}}_h^0 $, 这里$ \mathit{\boldsymbol{f}} $是已知函数. 基函数表示为$ \mathit{\boldsymbol{\varphi}}_i $, 则$ \mathit{\boldsymbol{V}}_h^0=span\{\mathit{\boldsymbol{\varphi}}_1, \mathit{\boldsymbol{\varphi}}_2, \dots, \mathit{\boldsymbol{\varphi}}_k\} $, $ k $是网格边数. 记$ \mathit{\boldsymbol{h}}^h = \sum\limits_{i=1}^{k} \alpha_i\mathit{\boldsymbol{\varphi}}_i $, 其中$ \mathit{\boldsymbol{\alpha}} = \big(\alpha_1, \alpha_2, \dots, \alpha_k\big) $. 定义非线性算子$ {\cal L} $$ {\Bbb R}^k $$ {\Bbb R}^k $的映射: $ {\cal L}(\mathit{\boldsymbol{\alpha}})=\big({\cal L}_1(\mathit{\boldsymbol{\alpha}}), \dots, {\cal L}_{k}(\mathit{\boldsymbol{\alpha}})\big) $,

其中$ j=1, \dots, k $. 这样问题就简化成对非线性代数方程$ \mathit{\boldsymbol{L}}(\mathit{\boldsymbol{\alpha}})={\bf{0}} $的求解. 对于算子$ \mathit{\boldsymbol{L}} $, 有

根据文献[18, 引理18.2], 我们得出结论, 在集合$ \{\mathit{\boldsymbol{\alpha}}\in {\Bbb R}^k:\; |\mathit{\boldsymbol{\alpha}}| \le r\} $中, 如果$ r > \|\mathit{\boldsymbol{f}}\| $, 则方程$ \mathit{\boldsymbol{L}}(\mathit{\boldsymbol{\alpha}})={\bf{0}} $至少有一个解.

下面, 讨论(4.5)式解的唯一性. 设$ \mathit{\boldsymbol{h}}^h_1 $$ \mathit{\boldsymbol{h}}^h_2 $是两个不同的解, 有

$ \mathit{\boldsymbol{Q}}^h=\mathit{\boldsymbol{h}}^h_1-\mathit{\boldsymbol{h}}^h_2 $, 有

因此, 可以得到$ \mathit{\boldsymbol{h}}^h_1=\mathit{\boldsymbol{h}}^h_2 $.

下面的两个引理将给出$ (\mathit{\boldsymbol{h}}_i^h, u_i^h) $的稳定性估计.

引理4.1  假设$ u_0\in H^1(\Omega) $, 则存在一个正的常数$ C_r $, 取决于截断函数$ R_r $中的参数$ r $, 使得

  (ⅰ) 在(4.4)中取$ v^h=\delta u_i^h\tau $, 对$ i=1, \cdots, j $$ (1\leq j\leq n ) $求和,

$ \begin{equation} \sum\limits_{i=1}^{j}\big\|\delta u_i^h\big\|^2\tau+\sum\limits_{i=1}^{j}\big(\lambda_i\nabla u_i^h, \nabla u_i^h-\nabla u_{i-1}^h\big) =\sum\limits_{i=1}^{j}\Big({\cal R}_r\big(\gamma(u_{i-1}^h)|\nabla \times\mathit{\boldsymbol{h}}_i^h|^2\big), \delta u_i^h\Big)\tau. \end{equation} $

对(4.6)式左端第二项处理(参考文献[6, 引理5, (ⅰ)]), 有

由于$ {\cal R}_r $是有界函数, 因此使用Hlöder不等式和Young不等式, 有

整理上述估计, 应用Gronwall不等式, 能够得到(ⅰ).

(ⅱ) 由上述(ⅰ)和$ u_j^h=u_0^h+ \sum\limits_{i=1}^{j}\delta u_i^h\tau $, 易得

(ⅲ) 对于任意的$ v^h\in W_h $, 有

证毕.

类似于引理4.1中(ⅲ)的证明, 在时间离散格式(3.4)中, 能够推导出

得到

假设函数$ \sigma $满足

因此

于是

这里需要注意的是, 对偶空间中的估计对于离散场是无效的, 因为变分形式的试探函数仅在$ H^1(\Omega) $的子空间$ W_h $中, 现在假设

$ \begin{equation} \begin{array}{ll} |\gamma(v^h_1)-\gamma(v^h_2)|\leq C\frac{\big|\big(v^h_1-v^h_2, v^h\big)\big|}{\big\|v^h\big\|_{H^1(\Omega)}}, \, \forall v^h_1, v^h_2\in W_h. \end{array} \end{equation} $

使用引理4.1 (ⅲ) 得到

引理4.2  令(2.2)–(2.4)及(4.7)式的假设均满足, 设$ \mathit{\boldsymbol{H}}_{s}\in H^1((0, T);\mathit{\boldsymbol{H}}(\mathit{\boldsymbol{curl}}, \Omega)) $, 则存在一个正的常数$ C $, 使得

  (ⅰ) 在(4.3)式中取$ \mathit{\boldsymbol{Q}}^h=\tau\delta\mathit{\boldsymbol{h}}_i^h $, 对$ i=1, \dots, j $求和,

结合(2.3)式, 对上式左端第一项估计,

对左端第二项估计, 由引理4.1中的(ⅲ)和(4.7)式, 易得

利用矢量场$ \mathit{\boldsymbol{B}} $的Lipschitz连续性及$ \mathit{\boldsymbol{H}}_{s}\in H^1((0, T);\mathit{\boldsymbol{H}}(\mathit{\boldsymbol{curl}}, \Omega)) $, 对右端处理, 得到

整理上述估计, 应用Gronwall不等式, 得到

(ⅱ) 在(4.3)式中取$ \mathit{\boldsymbol{Q}}^h=\tau\mathit{\boldsymbol{h}}_i^h $, 对$ i=1, \dots, j $$ (1\leq j\leq n ) $求和,

对于上式左端第一项, 由(2.3)和(2.4)式可以得到

对于左端第二项, 有

使用Cauchy不等式, 结合假设条件$ \mathit{\boldsymbol{H}}_{s}\in H^1((0, T);\mathit{\boldsymbol{H}}(\mathit{\boldsymbol{curl}}, \Omega)) $及矢量场$ \mathit{\boldsymbol{B}} $的Lipschitz连续性, 得到

因此, 有

证毕.

5 离散格式解的收敛性

本节的目的是证明当$ \tau, h\rightarrow 0 $时全离散格式解的收敛性. 首先确定时间步长$ \tau $, 令网格尺寸$ h $趋向于零. 这样得到时间离散格式(3.3)–(3.4)式的解$ (\mathit{\boldsymbol{h}}_i, u_i) $. 具体过程见定理5.1和5.2. 第二步, 令$ \tau\rightarrow 0 $, 证明问题(3.3)–(3.4)的解$ (\mathit{\boldsymbol{h}}_i, u_i) $收敛于(3.1)–(3.2)式的弱解$ (\mathit{\boldsymbol{H}}, u) $. 收敛结果已在定理3.2完成.

不失一般性, 设$ {\cal T}_1\prec {\cal T}_2\prec \cdots {\cal T}_k\prec \cdots $$ \Omega $中均匀且形状规则的网格序列, 使得$ \lim\limits_{k\rightarrow \infty} h_k =0 $, $ {\cal T}_{k+1} $表示网格$ {\cal T}_k $的加密. 为了明确离散解对$ {\cal T}_k $的依赖程度, 对离散解标注一个上标, 记为$ \mathit{\boldsymbol{h}}_i^{(k)} :=\mathit{\boldsymbol{h}}_i^h, $$ u_i^{(k)} :=u_i^h. $ 定义空间$ \mathit{\boldsymbol{V}}_h $, $ \mathit{\boldsymbol{V}}_h^0 $$ W_h $, 这取决于网格划分$ {\cal T}_k $, 分别为$ \mathit{\boldsymbol{V}}_h^k $, $ \mathit{\boldsymbol{V}}_h^{0, k} $$ W_h^k $. (4.1)和(4.2)分别定义了投影算子$ \mathit{\boldsymbol{\pi}}^k:\mathit{\boldsymbol{H}}(\mathit{\boldsymbol{curl}}, \Omega)\rightarrow \mathit{\boldsymbol{V}}_h^k $$ \Pi^k:H^1(\Omega)\rightarrow W_h^k $. 使用这些符号, 全离散问题(4.3)–(4.4)的格式如下

$ \begin{equation} \big(\delta\mathit{\boldsymbol{B}}(\mathit{\boldsymbol{h}}_i^{(k)}+\mathit{\boldsymbol{\pi}}^k\mathit{\boldsymbol{H}}_{s, i}), \mathit{\boldsymbol{Q}}^{(k)}\big)+ \big(\gamma(u_{i-1}^{(k)})\nabla \times\mathit{\boldsymbol{h}}_i^{(k)}, \nabla \times\mathit{\boldsymbol{Q}}^{(k)}\big)=0, \end{equation} $

$ \begin{equation} \big(\delta u_i^{(k)}, v^{(k)}\big)+\big(\lambda_i\nabla u_i^{(k)}, \nabla v^{(k)}\big)=\Big({\cal R}_r\big(\gamma(u_{i-1}^{(k)})|\nabla \times\mathit{\boldsymbol{h}}_i^{(k)}|^2\big), v^{(k)}\Big), \end{equation} $

其中, $ \forall\mathit{\boldsymbol{Q}}^{(k)}\in \mathit{\boldsymbol{V}}_h^{0, k} $, $ \forall v^{(k)}\in W_h^k $.

利用非线性矢量场$ \mathit{\boldsymbol{B}} $的单调性和Minty-Browder方法[8, 18]克服非线性项收敛的困难, 可以证明定理5.1. 下面部分中的序列仍然是子序列, 使用相同的符号表示.

定理5.1  令引理4.1–4.2的假设均满足. 设$ \gamma(s) $是全局Lipschitz连续函数, 对于$ 1\leq i\leq n $, 当$ k\rightarrow \infty $时, 有

  (ⅰ) 由引理4.1可知

因此, 根据$ H^1 (\Omega) $的自反性和Sobolev空间的紧嵌入, 存在一个子列和$ u_i \in H^1(\Omega) $, 使得

利用函数$ \gamma $的Lipschitz连续性, 得到

所以$ \gamma\big(u_i^{(k)}\big)\rightarrow \gamma\big(u_i\big) $ a.e. 在$ \Omega $中.

(ⅱ) 根据投影算子(4.1)的定义, 有

(ⅲ) 利用引理4.2得到

因此, 可以得到一个子列$ \{\mathit{\boldsymbol{h}}^{(k)}_i\}_{k=1}^\infty $, 其弱收敛于某个$ \mathit{\boldsymbol{h}}_i\in\mathit{\boldsymbol{H}}_0(\mathit{\boldsymbol{curl}}, \Omega) $. 利用$ \mathit{\boldsymbol{H}}_0(\mathit{\boldsymbol{curl}}, \Omega) $$ \mathit{\boldsymbol{L}}^2(\Omega) $的紧单射, 以及$ \mathit{\boldsymbol{L}}^2(\Omega) $的紧性, 能够找到一个子列(参考文献[12]), 使得

(ⅳ) 利用矢量场的Lipschitz连续性和引理4.2中的稳定性估计, 得到

$ \begin{equation} \|\mathit{\boldsymbol{B}}(\mathit{\boldsymbol{h}}_i^{(k)}+\mathit{\boldsymbol{\pi}}^k\mathit{\boldsymbol{H}}_{s, i})\|\leq C\|\mathit{\boldsymbol{h}}_i^{(k)}+\mathit{\boldsymbol{\pi}}^k\mathit{\boldsymbol{H}}_{s, i}\|\leq C. \end{equation} $

故存在一个子列和$ \mathit{\boldsymbol{U}} \in \mathit{\boldsymbol{L}}^2(\Omega) $, 使得$ k\rightarrow \infty $时, 有

根据矢量场$ \mathit{\boldsymbol{B}} $的单调性, 下述不等式对于任意的$ \mathit{\boldsymbol{w}}\in L^2((0, T);\mathit{\boldsymbol{L}}^2(\Omega)) $及任意非负$ \zeta\in C^\infty_0(\Omega) $成立,

$ \begin{equation} \big(\mathit{\boldsymbol{B}}(\mathit{\boldsymbol{h}}_i^{(k)}+\mathit{\boldsymbol{\pi}}^k\mathit{\boldsymbol{H}}_{s, i})-\mathit{\boldsymbol{B}}(\mathit{\boldsymbol{w}}+\mathit{\boldsymbol{H}}_{s, i}), \zeta((\mathit{\boldsymbol{h}}_i^{(k)}+\mathit{\boldsymbol{\pi}}^k\mathit{\boldsymbol{H}}_{s, i})-(\mathit{\boldsymbol{w}}+\mathit{\boldsymbol{H}}_{s, i}))\big) \geq 0. \end{equation} $

(5.4)式可以分解成四项

$ P_1 $可以写成如下形式

因此

对其余各项取极限$ k\rightarrow \infty $, 有

返回到(5.4)式, 有

由于$ \mathit{\boldsymbol{w}} $是任意选取的, 令$ \mathit{\boldsymbol{w}}=\mathit{\boldsymbol{h}}_i+\varepsilon\mathit{\boldsymbol{h}}' $, 其中$ \varepsilon>0 $, $ \mathit{\boldsymbol{h}}'\in \mathit{\boldsymbol{L}}^2(\Omega) $, 得到

$ \varepsilon $趋向于$ 0 $, 则

由于$ \mathit{\boldsymbol{h}}' $的选取是任意的, 令$ \mathit{\boldsymbol{h}}'=-\mathit{\boldsymbol{h}}' $, 上面的不等式反向成立, 即

对于任意的$ \mathit{\boldsymbol{h}}'\in \mathit{\boldsymbol{L}}^2(\Omega) $和非负$ \zeta\in C^\infty_0(\overline{\Omega}) $成立. 因此, 得到在$ \Omega $中, 几乎处处$ \mathit{\boldsymbol{U}}=\mathit{\boldsymbol{B}}(\mathit{\boldsymbol{h}}_i+\mathit{\boldsymbol{H}}_{s, i}) $, 即在$ \mathit{\boldsymbol{L}}^2(\Omega) $$ \mathit{\boldsymbol{B}}_i^k\rightharpoonup\mathit{\boldsymbol{B}}(\mathit{\boldsymbol{h}}_i+\mathit{\boldsymbol{H}}_{s, i}) $.

(ⅴ) 对于任意的$ \mathit{\boldsymbol{Q}}\in \mathit{\boldsymbol{H}}_0(\mathit{\boldsymbol{curl}}, \Omega) $, 有

显然在$ \mathit{\boldsymbol{L}}^2(\Omega) $$ \nabla \times\mathit{\boldsymbol{h}}_i^{(k)}\rightharpoonup \nabla \times\mathit{\boldsymbol{h}}_i $.

接下来, 需要讨论如下不等式

$ \begin{eqnarray} 0&\leq&\gamma_*\int_{\Omega}\big|\nabla \times\mathit{\boldsymbol{h}}_i^{(k)}-\nabla \times\mathit{\boldsymbol{h}}_i\big|^2{\rm d}\mathit{\boldsymbol{x}}\leq\int_{\Omega}\gamma(u_{i-1}^{(k)})\big|\nabla \times\mathit{\boldsymbol{h}}_i^{(k)}-\nabla \times\mathit{\boldsymbol{h}}_i\big|^2{\rm d}\mathit{\boldsymbol{x}}{}\\ &=&\big(\gamma(u_{i-1}^{(k)})\nabla \times\mathit{\boldsymbol{h}}_i^{(k)}, \nabla \times\mathit{\boldsymbol{h}}_i^{(k)}\big)+\big(\gamma(u_{i-1}^{(k)})\nabla \times\mathit{\boldsymbol{h}}_i, \nabla \times\mathit{\boldsymbol{h}}_i\big){}\\ &&-2\big(\gamma(u_{i-1}^{(k)})\nabla \times\mathit{\boldsymbol{h}}_i^{(k)}, \nabla \times\mathit{\boldsymbol{h}}_i\big). \end{eqnarray} $

根据本定理中的(ⅰ), (ⅱ), 应用Lebesgue控制定理, 有(参考文献[16, 定理1 (ⅹⅲ)])

$ \begin{equation} \lim\limits_{k\rightarrow \infty}\big(\gamma(u_{i-1}^{(k)})\nabla \times\mathit{\boldsymbol{h}}_i, \nabla \times\mathit{\boldsymbol{h}}_i\big)=\big(\gamma(u_{i-1})\nabla \times\mathit{\boldsymbol{h}}_i, \nabla \times\mathit{\boldsymbol{h}}_i\big) \end{equation} $

$ \begin{equation} \lim\limits_{k\rightarrow \infty}\big(\gamma(u_{i-1}^{(k)})\nabla \times\mathit{\boldsymbol{h}}_i^{(k)}, \nabla \times\mathit{\boldsymbol{h}}_i\big)=\big(\gamma(u_{i-1})\nabla \times\mathit{\boldsymbol{h}}_i, \nabla \times\mathit{\boldsymbol{h}}_i\big). \end{equation} $

对于固定的时间步长$ \tau $, 利用(5.1), (5.3), (2.3)和本定理中的(ii)–(iii), 得到

$ \begin{eqnarray} &&\lim\limits_{k\rightarrow \infty}\big(\gamma(u_{i-1}^{(k)})\nabla \times\mathit{\boldsymbol{h}}_i^{(k)}, \nabla \times\mathit{\boldsymbol{h}}_i^{(k)}\big) {}\\ &=&-\lim\limits_{k\rightarrow \infty}\big(\delta\mathit{\boldsymbol{B}}(\mathit{\boldsymbol{h}}_i^{(k)}+\mathit{\boldsymbol{\pi}}^k\mathit{\boldsymbol{H}}_{s, i}), \mathit{\boldsymbol{h}}_i^{(k)}\big) {}\\ &=&-\lim\limits_{k\rightarrow \infty}\big(\delta\mathit{\boldsymbol{B}}(\mathit{\boldsymbol{h}}_i^{(k)}+\mathit{\boldsymbol{\pi}}^k\mathit{\boldsymbol{H}}_{s, i}), \mathit{\boldsymbol{h}}_i^{(k)}-\mathit{\boldsymbol{h}}_i\big) -\big(\delta\mathit{\boldsymbol{B}}(\mathit{\boldsymbol{h}}_i+\mathit{\boldsymbol{H}}_{s, i}), \mathit{\boldsymbol{h}}_i\big) {}\\ &&-\lim\limits_{k\rightarrow \infty}\big(\delta\mathit{\boldsymbol{B}}(\mathit{\boldsymbol{h}}_i^{(k)}+\mathit{\boldsymbol{\pi}}^k\mathit{\boldsymbol{H}}_{s, i})-\delta\mathit{\boldsymbol{B}}(\mathit{\boldsymbol{h}}_i+\mathit{\boldsymbol{H}}_{s, i}), \mathit{\boldsymbol{h}}_i\big) {}\\ &=&-\big(\delta\mathit{\boldsymbol{B}}(\mathit{\boldsymbol{h}}_i+\mathit{\boldsymbol{H}}_{s, i}), \mathit{\boldsymbol{h}}_i\big)=\big(\gamma(u_{i-1})\nabla \times\mathit{\boldsymbol{h}}_i, \nabla \times\mathit{\boldsymbol{h}}_i\big), \end{eqnarray} $

整理(5.5)–(5.8)式, 有

$ \begin{equation} \lim\limits_{k\rightarrow \infty}\int_{\Omega}\big|\nabla \times\mathit{\boldsymbol{h}}_i^{(k)}-\nabla \times\mathit{\boldsymbol{h}}_i\big|^2{\rm d}\mathit{\boldsymbol{x}}=0, \end{equation} $

故在$ \mathit{\boldsymbol{L}}^2(\Omega) $$ \nabla \times\mathit{\boldsymbol{h}}_i^{(k)}\rightarrow \nabla \times\mathit{\boldsymbol{h}}_i $.

下面的定理表明, 当固定时间步长$ \tau $并使网格尺寸$ h $趋于零时, 得到Maxwell方程与热传导方程耦合的时间离散模型.

定理5.2  令定理5.1的假设均满足. 对于任意的$ i=1, \cdots, n $, $ \mathit{\boldsymbol{h}}_i\in \mathit{\boldsymbol{H}}_0(\mathit{\boldsymbol{curl}}, \Omega) $$ u_i\in H^1(\Omega) $$ (3.3) $$ (3.4) $式的唯一解.

  首先证明$ (\mathit{\boldsymbol{h}}_i, u_i) $求解$ (3.3) $式. 利用(5.1)式及$ \mathit{\boldsymbol{V}}_h^{0, l}\subset\mathit{\boldsymbol{V}}_h^{0, k} $, 对于$ l\leq k $, 有

对上式取极限$ k\rightarrow \infty $, 有

$ l\rightarrow \infty $时, 由$ \mathit{\boldsymbol{V}}_h^{0, l} $$ \mathit{\boldsymbol{H}}_0(\mathit{\boldsymbol{curl}}, \Omega) $的稠密性, 得到

$ \begin{equation} \begin{array}{ll} \big(\delta\mathit{\boldsymbol{B}}(\mathit{\boldsymbol{h}}_i+\mathit{\boldsymbol{H}}_{s, i}), \mathit{\boldsymbol{Q}}\big)+ \big(\gamma(u_{i-1})\nabla \times\mathit{\boldsymbol{h}}_i, \nabla \times\mathit{\boldsymbol{Q}}\big)=0\; \; \forall \mathit{\boldsymbol{Q}}\in \mathit{\boldsymbol{H}}_0(\mathit{\boldsymbol{curl}}, \Omega). \end{array} \end{equation} $

接下来, 为了证明$ \mathit{\boldsymbol{h}}_i $$ u_i $求解(3.4)式, 对于$ l\leq k $, 有

$ \begin{equation} \big(\delta u_i^{(k)}, v^{(l)}\big)+\big(\lambda_i\nabla u_i^{(k)}, \nabla v^{(l)}\big)=\Big({\cal R}_r\big(\gamma(u_{i-1}^{(k)})|\nabla \times\mathit{\boldsymbol{h}}_i^{(k)}|^2\big), v^{(l)}\Big)\; \; \forall v^{(l)}\in W_h^l . \end{equation} $

利用定理5.1和Lebesgue控制收敛定理, 注意到函数$ {\cal R}_r $是连续且有界的, 得到

对(5.11)式取极限, 有

$ l\rightarrow \infty $时, 由$ W_h^l $$ H^1(\Omega) $中的稠密性, 得到

证毕.

6 数值实验

本节将给出非线性感应加热问题的数值实验.

实验6.1  本次实验的目的是在固定时间步长下, 检验当网格尺寸$ h $减小时误差的变化. 计算区域$ \Omega $是一个单位立方体, 包含有限数量的标准四面体.

选取满足条件的非线性函数

其中磁场$ \hat{ \mathit{\boldsymbol{H}}}= \mathit{\boldsymbol{H}} + \mathit{\boldsymbol{H}}_s $. 初始值和参数为

选取源场为

将时间$ [0, 1] $均分成$ 100 $份, 使用COMSOL软件中使用2阶Nédélac单元求解$ \mathit{\boldsymbol{H}} $, 使用2阶Lagrange单元求解$ u $. 在每个时间步上使用MUMPS求解器来求解整个系统. 将立方体的每条边均匀剖分成64份, 在这种情况下, 得到的解作为参考解, 分别表示为$ \mathit{\boldsymbol{H}}_{{\rm ref}} $$ u_{{\rm ref}} $.

为了表明文中格式解的收敛性, 计算$ h = 1/2, 1/4, 1/8, 1/16, 1/32 $的数值解$ \mathit{\boldsymbol{H}}_{{\rm n}} $$ u_{{\rm n}} $, 将数值解与$ \mathit{\boldsymbol{H}}_{{\rm ref}} $$ u_{{\rm ref}} $进行比较. 选取$ \Omega $中的特定测量点以及特定时间$ t_i=0.01i $, 其中$ i=1, \cdots, 10 $. $ \mathit{\boldsymbol{H}}_n $$ \mathit{\boldsymbol{H}}_{{\rm ref}} $, $ u_n $$ u_{{\rm ref}} $之间的相对误差为

其中$ P $是特定测量点的集合. 随着网格尺寸$ h $的减小, 相对误差的变化如图 2所示. 图 2(a)中的回归线表示为$ \log_2 ( {\rm rel} \, \mathit{\boldsymbol{H}}_n ) = 1.784 \log_2 ( \tau ) + 2.391 $, 图 2(b)中的回归线表示为$ \log_2 ( {\rm rel} \, u_n ) = 2.202 \log_2 ( \tau ) -2.485 $, 故数值解收敛于参考解.

图 2

图 2   网格尺寸$ h $与相对误差的对数标度图


实验6.2  实验模型见图 1. 本次实验将给出一些数值模拟, 演示当改变电流频率$ \omega $时, 涡流和温度的分布情况.

交流电源为

通过线圈的交流电源产生一个通过工件的交变磁场$ \mathit{\boldsymbol{H}}_s $, 该磁场在工件中产生涡流和焦耳热. 计算区域$ \Omega $是工件. 时间为$ [0, 1] $, 时间步长$ \tau = 0.2 \pi/ \omega $.

$ t=1.0 $时, 在不同系数$ \omega $下的涡流密度$ \nabla \times \mathit{\boldsymbol{H}} $的变化见图 3, 图 4图 5, 温度$ u $的变化如图 6, 图 7图 8所示.

图 3

图 3   涡流密度$ \omega = \pi, \; t = 1.0 $


图 4

图 4   涡流密度$ \omega = 10\pi, \; t = 1.0 $


图 5

图 5   涡流密度$ \omega = 100\pi, \; t = 1.0 $


图 6

图 6   温度分布$ \omega = \pi, \; t = 1.0 $


图 7

图 7   温度分布$ \omega = 10\pi, \; t = 1.0 $


图 8

图 8   温度分布$ \omega = 100\pi, \; t = 1.0 $


数值模拟结果表明, 涡流的大小随着与工件表面距离的增大而减小. 当电流频率$ \omega $增加时, 涡流的分布集中在工件的外侧. 此外, 温度的分布也受$ \omega $的影响. 工件的整体温度随着$ \omega $的增加而升高, 并且工件内侧与外侧之间的温度差逐渐增大. 这表明温度的变化与涡流的分布是一致的.

7 结论

本文研究了非线性Maxwell方程与热传导方程耦合的全离散有限元解耦格式. 证明了全离散格式解的存在唯一性. 同时也证明了全离散格式的解收敛于连续耦合问题的解. 通过两组数值实验来支持文中的理论结果, 实验结果表明, 理论结果与感应加热模型的实际研究结果一致.

参考文献

Akrivis G , Larsson S .

Linearly implicit finite element methods for the time dependent Joule heating problem

BIT Numerical Mathematics, 2005, 45: 429- 442

DOI:10.1007/s10543-005-0008-1      [本文引用: 1]

Barglik J , Doležel I , Karban P , et al.

Modelling of continual induction hardening in quasi-coupled formulation

Compel-the International Journal for Computation and Mathematics in Electrical and Electronic Engineering, 2005, 24 (1): 251- 260

DOI:10.1108/03321640510571273     

Bermúdez A , Gómez D , Muniz M C , et al.

Transient numerical simulation of a thermoelectrical problem in cylindrical induction heating furnaces

Advances in Computational Mathematics, 2007, 26: 39- 62

DOI:10.1007/s10444-005-7470-9      [本文引用: 1]

Bień M .

Global solutions of the non-linear problem describing Joule's heating in three space dimensions

Mathematical Methods in the Applied Sciences, 2005, 28 (9): 1007- 1030

DOI:10.1002/mma.559      [本文引用: 1]

Bossavit A , Emson C , Mayergoyz I D . Méthodes Numériques en Électromagnétisme. Paris: Eyrolles, 1991

Chovan J , Geuzaine C , Slodička M .

$\bm{A}$-$\phi$ formulation of a mathematical model for the induction hardening process with a nonlinear law for the magnetic field

Computer Methods in Applied Mechanics and Engineering, 2017, 321: 294- 315

DOI:10.1016/j.cma.2017.03.045      [本文引用: 2]

Elliott C M , Larsson S .

A finite element model for the time-dependent Joule heating problem

Mathematics of Computation, 1995, 62 (212): 1433- 1453

[本文引用: 1]

Evans L C. Graduate Studies in Mathematics: Partial Differential Equations. Second ed. New York: American Mathematical Society, 2010

[本文引用: 1]

Hömberg D .

A mathematical model for induction hardening including mechanical effects

Nonlinear Analysis: Real World Applications, 2004, 5 (1): 55- 90

DOI:10.1016/S1468-1218(03)00017-8      [本文引用: 1]

Hömberg D , Petzold T , Rocca E .

Analysis and simulations of multifrequency induction hardening

Nonlinear Analysis: Real World Applications, 2005, 22: 84- 97

[本文引用: 1]

Kang T , Wang R , Zhang H .

Potential field formulation based on decomposition of the electric field for a nonlinear induction hardening model

Communications in Applied Mathematics and Computational Science, 2019, 14 (2): 175- 205

DOI:10.2140/camcos.2019.14.175      [本文引用: 1]

Kang T , Wang R , Zhang H .

Fully discrete $\bm{T}$-$\psi$ finite element method to solve a nonlinear induction hardening problem

Numerical Methods for Partial Differential Equations, 2021, 37: 546- 582

DOI:10.1002/num.22540      [本文引用: 1]

Nečas J . Introduction to the Theory of Nonlinear Elliptic Equations. Chichester: John Wiley & Sons Ltd, 1986

[本文引用: 1]

Nédélec J C .

Mixed finite elements in $R^3$

Numerische Mathematick, 1980, 35: 315- 341

DOI:10.1007/BF01396415      [本文引用: 1]

Roubíček T . Nonlinear Partial Differential Equations with Applications. Berlin: Birkhöuser, 2005

Slodička M , Chovan J .

Solvability for induction hardening including nonlinear magnetic field and controlled Joule heating

Applicable Analysis, 2017, 96: 2780- 2799

DOI:10.1080/00036811.2016.1243661      [本文引用: 4]

Sun D , Manoranjan V S , Yin H M .

Numerical solutions for a coupled parabolic equations arising induction heating processes

Conference Publications, 2007, 2007 (Special): 956- 964

[本文引用: 1]

Vajnverg M . Variational Method and Method of Monotone Operators in the Theory of Nonlinear Equations. New York: John Wiley, 1973

[本文引用: 7]

Yin H M .

On a nonlinear Maxwell's system in quasi-stationary electromagnetic fields

Mathematical Models and Methods in Applied Sciences, 2004, 14 (10): 1521- 1539

DOI:10.1142/S0218202504003787      [本文引用: 1]

Yin H M .

Regularity of weak solution to Maxwell's equations and applications to microwave heating

Journal of Differential Equations, 2004, 200 (1): 137- 161

DOI:10.1016/j.jde.2004.01.010     

Yin H M , Wei W .

Regularity of weak solution for a coupled system arising from a microwave heating model

European Journal of Applied Mathematics, 2014, 25: 117- 131

DOI:10.1017/S0956792513000326      [本文引用: 1]

/