CN110059327A - 一种基于辐射换热的三维有限元模拟方法 - Google Patents
一种基于辐射换热的三维有限元模拟方法 Download PDFInfo
- Publication number
- CN110059327A CN110059327A CN201811432872.XA CN201811432872A CN110059327A CN 110059327 A CN110059327 A CN 110059327A CN 201811432872 A CN201811432872 A CN 201811432872A CN 110059327 A CN110059327 A CN 110059327A
- Authority
- CN
- China
- Prior art keywords
- heat transfer
- finite element
- radiation heat
- iteration
- equation
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Pending
Links
- 230000005855 radiation Effects 0.000 title claims abstract description 41
- 238000000034 method Methods 0.000 title claims abstract description 35
- 238000004088 simulation Methods 0.000 title claims abstract description 9
- 230000006870 function Effects 0.000 claims abstract description 32
- 239000011159 matrix material Substances 0.000 claims abstract description 27
- 239000013598 vector Substances 0.000 claims abstract description 6
- 238000002940 Newton-Raphson method Methods 0.000 claims description 10
- 238000003475 lamination Methods 0.000 claims description 6
- 238000005516 engineering process Methods 0.000 abstract description 3
- 238000002076 thermal analysis method Methods 0.000 abstract description 2
- 238000004458 analytical method Methods 0.000 abstract 1
- 230000008569 process Effects 0.000 description 7
- 238000004364 calculation method Methods 0.000 description 4
- 230000014509 gene expression Effects 0.000 description 4
- 230000008878 coupling Effects 0.000 description 3
- 238000010168 coupling process Methods 0.000 description 3
- 238000005859 coupling reaction Methods 0.000 description 3
- 238000010586 diagram Methods 0.000 description 2
- 239000000463 material Substances 0.000 description 2
- 230000008901 benefit Effects 0.000 description 1
- 230000005540 biological transmission Effects 0.000 description 1
- 238000004422 calculation algorithm Methods 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 238000004891 communication Methods 0.000 description 1
- 238000010276 construction Methods 0.000 description 1
- 230000007812 deficiency Effects 0.000 description 1
- 238000004519 manufacturing process Methods 0.000 description 1
- 230000009466 transformation Effects 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
- G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
- G06F2119/08—Thermal analysis or thermal optimisation
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Computer Hardware Design (AREA)
- Evolutionary Computation (AREA)
- Geometry (AREA)
- General Engineering & Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Complex Calculations (AREA)
Abstract
本发明属于三维热分析数值求解技术领域，涉及一种基于辐射换热的三维有限元模拟方法。本发明首先对要进行热分析的器件进行建模，然后将辐射热平衡方程与热传导控制方程耦合求解，采用伽辽金残数加权的方法，得到热辐射边界条件的有限元弱形式。接着采用四面体网格剖分模型，选择二阶叠层基函数，离散有限元弱形式方程，配合Newton‑Raphson迭代方法得到有限元单元矩阵和右端向量，集成最终的方程组，最后运用科学的非线性收敛判据，经过不断迭代，快速准确地得到最终的数值计算结果。
Description
技术领域
本发明属于三维热分析数值求解技术领域，具体涉及一种基于辐射换热的三维有限元模拟方法。
背景技术
行波管是微波器件中高功率、宽频带放大器的主要管型，大量用作卫星通信、电子对抗和雷达系统的放大级或发射功率源。随着新材料、新工艺和新技术的不断出现，行波管的性能和技术指标有了很大提高，但随之也带来了高的热负载。电子枪作为行波管的核心部件，对于行波管寿命和可靠性都有重要的影响。电子枪在工作状态下，其阴极要发射电子注，必须要达到足够高的温度，一般高达1000摄氏度。在如此高的温度下，加上电子枪内部的真空环境，辐射换热是一个极其重要的热传递过程。
目前，在数值计算领域，有关的辐射换热理论主要分为两种：一种是封闭环境内的辐射热平衡方程与热传导控制方程迭代求解的显式计算方法，该种方式计算方便快捷，但是对于某些特定的情况(比如边界条件缺失)会导致问题不可求解；一种是辐射热平衡方程与热传导控制方程的耦合求解，这是一种隐式求解方法，有的文献耦合求解会采用Newton-Raphson等非线性迭代方法，但是其耦合矩阵是非对称的，而且采用的是低阶基函数，在网格加密的情况下造成极大的内存和计算负担。
发明内容
针对上述存在问题或不足，为解决现有的辐射热平衡方程与热传导控制方程的耦合求解效率低的问题，本发明提供了一种基于辐射换热的三维有限元模拟方法，该方法从辐射能量的本质出发，经过巧妙的变换，使得最终的耦合矩阵变为对称矩阵，并且采用二阶叠层基函数以及牛顿-拉夫森(Newton-Raphson)迭代方法来进行非线性的有限元求解，可以很快地求得高精度的数值模拟结果。
其具体技术方案，包括以下步骤：
A.对目标器件进行建模，建立对应的几何结构模型；
B.从辐射换热的本质出发，得到辐射换热问题的有限元弱形式；
C.采用四面体网格剖分求解域；
D.选择叠层基函数，离散B中得到的有限元弱形式，得到辐射换热问题的有限元方程组；
E.对步骤D中的有限元非线性方程组不断地进行迭代，直到其温度值满足设定的收敛规则。
进一步优选，所述步骤D中选择二阶叠层基函数，相比于插值高阶基函数，叠层基函数构造方法更加简便，而且对于后面的有限元处理过程也有极大的好处，提高了有限元求解的精度。此外，步骤E中的迭代方法选用Newton-Raphson非线性迭代方法，通过这种方法的使用，使得辐射换热的高度非线性难题迎刃而解，并且可以很快地达到收敛，极大地提高了求解速度。
Newton-Raphson方法最终的非线性方程组迭代形式如下：
J(q)ΔA(q)＝F-S(q)A(q) (1)
其中J(q)是雅可比矩阵，是本发明需要求解的最重要矩阵，ΔA(q)是前后两次迭代的温度差值，F是右端向量，q是迭代次数，S(q)是采用Newton-Raphson方法之前的有限元初始矩阵项，A(q)是前一次迭代的温度值。
所述收敛规则是||F-S(q+1)A(q+1)||＜ε或者(ε是人为设定的收敛精度值，S(q+1)是采用Newton-Raphson方法之后的有限元初始矩阵项，A(q+1)是后一次迭代的温度值)。
本发明针对辐射换热在有限元方法中的应用问题，将辐射热平衡方程与热传导控制方程联立，得到有限元弱形式，并采用高阶叠层基函数离散成非线性方程组，与Newton-Raphson迭代方法结合，通过设定的收敛判据，不断地迭代，直到求得最终的温度值。
与现有技术相比，本发明可以准确、快速地求解高度非线性辐射换热问题，并且解决了现有技术关于这方面描述不足以及一些特定解决方法效率低下的问题。
附图说明
图1是本发明的流程图；
图2是辐射换热示意图；
图3是四节点四面体单元示意图；
图4是Newton-Raphson非线性迭代方法的流程图；
具体实施方式
下面结合附图和具体实施例来详细描述本发明的技术方案。
参照图1，一种基于辐射换热的三维有限元模拟方法，包括以下步骤：
A.对目标器件进行建模，建立对应的几何结构模型。
B.从辐射换热的本质出发，得到辐射换热问题的有限元弱形式。
要解决辐射换热问题，首先要构建封闭包壳。如图2所示，考虑封闭包壳中的某个面a，定义其有效辐射
其中Ra表示有效辐射，σ是Stefan-Boltzmann常数，εa是热辐射发射率，表示面a的平均温度，Ha表示入射辐射。同样的，可以定义净辐射
其中qa表示净辐射热流密度。由(2)、(3)得：
此外，入射辐射还满足如下的能量方程：
其中Sa、Sb表示辐射面的面积，Fab、Fba表示辐射面之间的角系数，f表示封闭包壳内的辐射面总数。由(2)和(5)可得：
本发明的关键是对(6)式进行巧妙的处理，在方程两侧分别乘以以保证后续的有限元矩阵是对称矩阵，考虑到角系数的性质，变换后的矩阵形式如下：
其中：
将式(4)和式(7)结合，可以得到：
其中：
并且Lab是Gab的逆矩阵。
求得了净辐射热流密度qa，就可以与热传导控制微分方程结合，进行有限元求解。先构建边值问题，包括热传导的控制微分方程和热辐射净热流密度两部分，具体如下方程所示：
其中u是求解域内的温度值，k是热传导系数，Q是内部产热量，ρ是密度，c是比热容，t是时间，n表示法向方向。采用伽辽金残数加权法，可以得到(11)式、(12)式的加权残数表达式如下：
其中表示残差项，Ω表示求解域，Γ表示辐射面。从而进一步可以得到
其中v1,v2是权函数，c1是任意实数。
这里定义面积分和体积分如下
(u,v)Ω＝∫Ω(u,v)dV (16)
<u,v>Γ＝∫Γ(u,v)dS (17)
其中u、v表示任意两个函数，V表示体积，S表示面积。
对于(15)式中的通过格林定理，可以展开为
由于v1,v2的任意性，令c1＝-1,v1＝v2，可得
为了利用弱形式得到问题的近似解，首先应该选择试函数Ni(x,y,z)来代替真实解，如下式，且必须满足必要的边界条件。
其中C0、Ci是任意实数，Ni为简单函数，例如低阶多项式。在伽辽金方法中，直接采用试函数自身作为权函数，即
v1(x,y,z)＝Ni(x,y,z) (21)
后面的有限元过程考虑稳态热传导，不考虑控制方程中的时间项，所以三维热传导方程弱形式写成
C.采用四面体网格剖分求解域；
采用四面体网格剖分求解域，剖分后的求解域被分割为三维四面体网格，从而将连续的几何结构空间转化为离散的网格空间。
D.选择叠层基函数，离散B中得到的有限元弱形式，得到辐射换热问题的有限元方程组；
如图3所示四面体单元中i,j,k,l代表四个顶点的编号，我们首先得到四个最基本的基函数：
式中
将(27)式、(28)式、(29)式和(30)式中的i,j,k,l轮换，得到aj,ak,al，bj,bk,bl，cj,ck,cl，dj,dk,dl。V为四面体的体积。
对于标量叠层高阶基函数的选择有如下规则叠层规则
式(32)(33)中，Wp表示所有基函数的集合，p表示选取基函数的阶数，Dim(Wp)表示基函数的个数，i,j,k,l表示上标。所以二阶叠层标量基函数，选择N1,N2,N3,4N1N2,4N1N3,4N2N3这六个基函数，系数4考虑计算方便。对于有限元过程来说，把域Ω离散为M个单元之后，如同(22)式所示的弱形式定积分，可以通过将每个单元的积分贡献简单相加，即
对于每一个单元来说，跟有限元方程组右端项有关系的的求解在很多有限元基础材料中都有介绍，在此不再赘述。本发明的关注点是用Newton-Raphson迭代方法来进行非线性辐射边界的处理，因此重点考虑的是有限元方程组左端矩阵的求解。
假设非线性方程组的形式为：
SA＝F (35)
其中S为左端矩阵，A为待求解向量，F为右端向量。
对于Newton-Raphson法是一种梯度算法，将(35)式的有限元方程组写成如下形式：
f(A)＝SA-F＝0 (36)
其中f(A)为非线性函数。应用Newton-Raphson方法则有
A(q+1)＝A(q)-[f′(A(q))]-1f(A(q)) (37)
其中上标q表示迭代次数，令J(q)＝f′(A(q)),并化简有
J(q)ΔA(q)＝F-S(q)A(q) (38)
其中J(q)为非线性函数f(A)对A的导数矩阵[称为雅可比(Jacobi)矩阵]的第q次迭代值，ΔA(q)前后两次迭代的温度差值，满足
ΔA(q)＝A(q+1)-A(q) (39)
其中A(q+1)为后一次求解的值，A(q)为前一次求解的值。通过求解线性方程组的方法求解方程组(38)式，可以得到ΔA(q)，然后根据(39)式可以获得A(q+1)，其流程图如图4所示。
Newton-Raphson迭代方法需要计算雅可比矩阵J,对于有限元计算来说，雅可比矩阵J可以由各单元的雅可比矩阵Je叠加组成。由(34)式，定义单元矩阵Sij：
令Sij＝Kij+Mij，其中
难点在于如何处理进行M矩阵的求解，下面具体阐述。首先考虑平均温度的定义：
其中Ta表示节点插值系数(包括三个顶点和三个中点)。所以二阶叠层基函数的平均温度表达式为：
考虑式(42)，为了实现降维，令所以：
所以：
同样地，雅克比矩阵为：
E.对步骤D中的有限元非线性方程组不断地进行迭代，直到其温度值满足设定的收敛规则。
在步骤D中已经完成所需矩阵J和S的求解，接下来只需要按照Newton-Raphson迭代方法所述，完成非线性方程组的构建。对于非线性方程组的迭代求解过程，本发明采用的收敛规则是||F-S(q+1)A(q+1)||＜ε或者(ε是人为设定的收敛精度值，S(q+1)是采用Newton-Raphson方法之后的有限元初始矩阵项，A(q+1)是后一次迭代的温度值)，前者的收敛判据更加准确，但是处理过程可能相对后者更加复杂，后者判据则相对简单，但精度可能会有所降低。
综上所述，本发明针对目前热辐射换热问题在有限元方法中的应用难题，提出了一套通用的有限元解决方法，采用高阶叠层基函数、Newton-Raphson非线性迭代方法以及高斯数值积分的使用，可以准确、快速地求解高度非线性辐射换热问题，并且解决了现有技术关于这方面描述不足以及一些特定解决方法的适用性不足的难题。
Claims (3)
1.一种基于辐射换热的三维有限元模拟方法，包括以下步骤：
A.对目标器件进行建模，建立对应的几何结构模型；
B.从辐射换热的本质出发，得到辐射换热问题的有限元弱形式；
C.采用四面体网格剖分求解域；
D.选择叠层基函数，离散B中得到的有限元弱形式，得到辐射换热问题的有限元方程组；
对于标量叠层高阶基函数的选择有如下规则叠层规则
式(1)(2)中，Wp表示所有基函数的集合，p表示选取基函数的阶数，Dim(Wp)表示基函数的个数，i,j,k,l表示上标；
E.对步骤D中的有限元非线性方程组不断地进行迭代，直到其温度值满足设定的收敛规则。
2.如权利要求1所述基于辐射换热的三维有限元模拟方法，其特征在于：所述步骤D中选择二阶叠层基函数。
3.如权利要求1所述基于辐射换热的三维有限元模拟方法，其特征在于：
所述步骤E中的迭代方法选用Newton-Raphson非线性迭代方法，Newton-Raphson方法最终的非线性方程组迭代形式如下：
J(q)ΔA(q)＝F-S(q)A(q) (3)
其中J(q)是雅可比矩阵，ΔA(q)是前后两次迭代的温度差值，F是右端向量，q是迭代次数，S(q)是采用Newton-Raphson方法之前的有限元初始矩阵项，A(q)是前一次迭代的温度值。
所述步骤E中收敛规则是||F-S(q+1)A(q+1)||＜ε或者ε是人为设定的收敛精度值，S(q+1)是采用Newton-Raphson方法之后的有限元初始矩阵项，A(q+1)是后一次迭代的温度值。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201811432872.XA CN110059327A (zh) | 2018-11-28 | 2018-11-28 | 一种基于辐射换热的三维有限元模拟方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201811432872.XA CN110059327A (zh) | 2018-11-28 | 2018-11-28 | 一种基于辐射换热的三维有限元模拟方法 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN110059327A true CN110059327A (zh) | 2019-07-26 |
Family
ID=67315751
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201811432872.XA Pending CN110059327A (zh) | 2018-11-28 | 2018-11-28 | 一种基于辐射换热的三维有限元模拟方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN110059327A (zh) |
Cited By (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110826266A (zh) * | 2019-10-01 | 2020-02-21 | 复旦大学 | 基于变换热辐射和热传导理论设计的热旋转器装置 |
| CN112836399A (zh) * | 2020-11-20 | 2021-05-25 | 电子科技大学 | 一种基于有限元算法的非线性接触热阻热分析求解方法 |
| CN113688495A (zh) * | 2021-07-01 | 2021-11-23 | 复旦大学 | 基于温度依赖变换热电场理论的热电转换器及其设计方法 |
| CN113761762A (zh) * | 2021-08-03 | 2021-12-07 | 西北核技术研究所 | 用于有限元数值模拟后验误差估计的平衡通量构造方法 |
Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN106354954A (zh) * | 2016-08-31 | 2017-01-25 | 电子科技大学 | 一种基于叠层基函数的三维力学模态仿真模拟方法 |
| CN107577857A (zh) * | 2017-08-28 | 2018-01-12 | 电子科技大学 | 一种基于热辐射边界条件的三维有限元模拟方法 |
-
2018
- 2018-11-28 CN CN201811432872.XA patent/CN110059327A/zh active Pending
Patent Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN106354954A (zh) * | 2016-08-31 | 2017-01-25 | 电子科技大学 | 一种基于叠层基函数的三维力学模态仿真模拟方法 |
| CN107577857A (zh) * | 2017-08-28 | 2018-01-12 | 电子科技大学 | 一种基于热辐射边界条件的三维有限元模拟方法 |
Non-Patent Citations (2)
| Title |
|---|
| PENG XIE 等: "《Thermal analysis of electron gun using finite element method》", 《2017 EIGHTEENTH INTERNATIONAL VACUUM ELECTRONICS CONFERENCE (IVEC)》 * |
| 谈和平等: "计算热辐射学的进展", 《科学通报》 * |
Cited By (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110826266A (zh) * | 2019-10-01 | 2020-02-21 | 复旦大学 | 基于变换热辐射和热传导理论设计的热旋转器装置 |
| CN110826266B (zh) * | 2019-10-01 | 2023-05-30 | 复旦大学 | 基于变换热辐射和热传导理论设计的热旋转器装置 |
| CN112836399A (zh) * | 2020-11-20 | 2021-05-25 | 电子科技大学 | 一种基于有限元算法的非线性接触热阻热分析求解方法 |
| CN112836399B (zh) * | 2020-11-20 | 2022-11-08 | 电子科技大学 | 一种基于有限元算法的非线性接触热阻热分析求解方法 |
| CN113688495A (zh) * | 2021-07-01 | 2021-11-23 | 复旦大学 | 基于温度依赖变换热电场理论的热电转换器及其设计方法 |
| CN113688495B (zh) * | 2021-07-01 | 2024-04-26 | 复旦大学 | 基于温度依赖变换热电场理论的热电转换器及其设计方法 |
| CN113761762A (zh) * | 2021-08-03 | 2021-12-07 | 西北核技术研究所 | 用于有限元数值模拟后验误差估计的平衡通量构造方法 |
| CN113761762B (zh) * | 2021-08-03 | 2023-10-20 | 西北核技术研究所 | 用于电场/温度有限元数值解的后验误差估计方法 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN110059327A (zh) | 一种基于辐射换热的三维有限元模拟方法 | |
| Ma et al. | Partially‐coupled least squares based iterative parameter estimation for multi‐variable output‐error‐like autoregressive moving average systems | |
| Toro et al. | ADER schemes for scalar non-linear hyperbolic conservation laws with source terms in three-space dimensions | |
| CN107577857A (zh) | 一种基于热辐射边界条件的三维有限元模拟方法 | |
| CN105468909B (zh) | 基于sod‑ps‑r算法的时滞电力系统机电振荡模式计算方法 | |
| CN109684740B (zh) | 一种基于混合网格及时间步长的电磁学多尺度计算方法 | |
| CN103412988B (zh) | 基于相移降阶模型周期结构的三维电磁场仿真模拟方法 | |
| CN107153721A (zh) | 一种运动目标下的辛时域有限差分电磁仿真方法 | |
| CN116484586A (zh) | 一种计算磁约束等离子体湍流特性的方法、系统、设备和存储介质 | |
| Wang et al. | A new family of exponential-based high-order DGTD methods for modeling 3-D transient multiscale electromagnetic problems | |
| Zingale et al. | Improved coupling of hydrodynamics and nuclear reactions via spectral deferred corrections | |
| Li et al. | Hermite spectral method for Fokker-Planck-Landau equation modeling collisional plasma | |
| Cai et al. | Optimal guidance for hypersonic reentry using inversion and receding horizon control | |
| CN117439731B (zh) | 基于同态加密的隐私保护大数据主成分分析方法及系统 | |
| CN108808706B (zh) | 基于sod-ps-ii-r算法的时滞电力系统机电振荡模式计算方法 | |
| CN112836399B (zh) | 一种基于有限元算法的非线性接触热阻热分析求解方法 | |
| Xue et al. | Efficient transient thermal analysis based on spectral element time domain method with curvilinear hexahedrons | |
| Liu et al. | A finite element domain decomposition combined with algebraic multigrid method for large-scale electromagnetic field computation | |
| Gamzina et al. | Thermo-mechanical stress in high-frequency vacuum electron devices | |
| Lin et al. | A discontinuous Galerkin method for two-temperature plasmas | |
| Sulaiman et al. | NUMERICAL SOLUTIONS OF NONLINEAR SECOND-ORDER TWO-POINT BOUNDARY VALUE PROBLEMS USING HALF-SWEEP SOR WITH NEWTON METHOD. | |
| CN107526869B (zh) | 一种基于函数逼近自适应三维微波管输入输出窗模型降阶的数值方法 | |
| CN113065201B (zh) | 一种考虑滑移修正的辐射平衡温度计算方法 | |
| Malek | Applications of nonstandard finite difference methods to nonlinear heat transfer problems | |
| Bihlo | Invariant meshless discretization schemes |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| RJ01 | Rejection of invention patent application after publication | ||
| RJ01 | Rejection of invention patent application after publication |
Application publication date: 20190726 |