CN107577857B - 一种基于热辐射边界条件的三维有限元模拟方法 - Google Patents

Abstract

本发明属于三维热分析数值求解技术领域，涉及一种基于热辐射边界条件的三维有限元模拟方法。本发明首先对要进行热分析的器件进行建模，然后将辐射边界条件引入热传导问题，采用伽辽金残数加权的方法，得到热辐射边界条件的有限元弱形式。接着采用四面体网格剖分模型，选择二阶叠层基函数，离散有限元弱形式方程，配合Newton‑Raphson迭代方法得到有限元单元矩阵和右端向量，集成最终的方程组，最后运用科学的非线性收敛判据，经过不断迭代，快速准确地得到最终的数值计算结果。

Description

一种基于热辐射边界条件的三维有限元模拟方法

技术领域

本发明属于三维热分析数值求解技术领域，涉及一种基于热辐射边界条件的三维有限元模拟方法。

背景技术

作为卫星大功率核心器件的空间行波管，在工作过程中不仅要受到真空高低温的影响，还要受到空间太阳辐射、红外辐射和真空冷黑背景温度的长期影响，这些热环境将影响空间行波管的电参数和热稳定性，对空间行波管的正常工作构成严重威胁。为了保证空间行波管各个部件保持在正常的工作温度，散热设计变得尤为重要。由于空间行波管暴露在外界真空环境，辐射式散热几乎是其唯一散热方式，因此有必要对空间行波管热辐射问题进行研究。

目前，在数值计算领域，有关的热传导差分和有限元法著作[如数值传热学，作者：陶文铨]以及其他相关文献，大都只涉及线性边界条件(狄利克雷边界、诺曼边界、罗宾边界)的处理方法，而对于非线性边界条件如热辐射只是大略提及，没有阐明具体的解决方法。西安交通大学的张博提出了一种间接迭代的方式来处理，但是没有从有限元的本质上处理辐射边界问题，只是将非线性辐射问题转化成了诺曼边界条件来求解，只能应用在特定领域，而且针对复杂模型也缺乏必要的验证。

发明内容

针对上述存在问题或不足，为解决热辐射边界条件在有限元方法中的应用难题，本发明提供了一种基于热辐射边界条件的三维有限元模拟方法，该方法采用二阶叠层基函数以及Newton-Raphson迭代方法来进行辐射边界条件的有限元求解，可以很快地求得高精度的数值模拟结果。

其具体技术方案，包括以下步骤：

A.对目标器件进行建模，建立对应的几何结构模型；

B.采用伽辽金加权残数法得到热辐射边界条件的有限元弱形式；

C.采用四面体网格剖分求解域；

D.选择基函数，离散步骤B中得到的有限元弱形式，得到热辐射边界条件的有限元方程；

E.对步骤D中的有限元非线性方程组不断地进行迭代，直到其温度值满足收敛规则,收敛规则是||F-S^(q+1)A^(q+1)||＜ε或者

ε为设定的收敛精度值。

进一步优选，所述的步骤D中，在基本的三维四面体有限元基函数基础上，构造了二阶叠层基函数，相比于插值高阶基函数，叠层基函数构造方法更加简便，而且对于后面的有限元处理过程也有极大的好处，提高了有限元求解的精度。此外，步骤D中的迭代方法选用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)是前一次迭代的温度值。因为辐射的高度非线性，矩阵的积分求解过程很难完成，本发明采用了高斯数值积分的方式，采用12阶高斯积分，很方便地求解了非线性积分矩阵项，并且保证了矩阵计算的精度。

本发明针对热辐射边界在有限元方法中的应用问题，提出了一种通用的处理方法，将热辐射边界条件与热传导控制方程联立，得到有限元弱形式，并采用高阶叠层基函数离散成非线性方程组，与Newton-Raphson迭代方法结合，通过收敛判据(||F-S^(q+1)A^(q+1)||＜ε或

不断地迭代，直到求得最终的温度值。

与现有技术相比，本发明可以准确、快速地求解高度非线性热辐射边界问题，并且解决了现有技术关于这方面描述不足以及一些特定解决方法的适用性不足的难题。

附图说明

图1是本发明的流程图；

图2是四节点四面体单元示意图；

图3是Newton-Raphson非线性迭代方法的流程图；

图4是下标对应关系图。

具体实施方式

下面结合附图和具体实施例来详细描述本发明的技术方案。

参照图1，一种基于热辐射边界条件的三维有限元模拟方法，包括以下步骤：

A.对目标器件进行建模，建立对应的几何结构模型。

B.采用伽辽金加权残数法得到热辐射边界条件的有限元弱形式。

要对热辐射边界条件进行有限元求解，就必须先构建边值问题，包括热传导的控制微分方程和热辐射边界条件两部分，具体如下方程所示：

其中

为拉普拉斯算子，u是求解域内的温度值，k是热传导系数，Q是内部产热量，ρ是密度，c是比热容，t是时间，ε是热辐射发射率，n是法向方向，σ是Stefan-Boltzmann常数，u_f是外界环境温度值。

采用伽辽金残数加权法，可以得到(2)式、(3)式的加权残数表达式如下

其中

表示残差项，Ω表示求解域，Γ表示边界条件。从而进一步可以得到

其中v₁,v₂是权函数，c₁是任意实数。

这里定义面积分和体积分如下

(u,v)_Ω＝∫_Ω(u,v)dV (7)

<u,v>_Γ＝∫_Γ(u,v)dS (8)

其中u、v表示任意两个函数，V表示体积，S表示面积。

对于(6)式中的

通过格林定理，可以展开为

由于v₁,v₂的任意性，令c₁＝-1,v₁＝v₂，可得

为了利用弱形式得到问题的近似解，首先应该选择试函数N_i(x,y,z)来代替真实解，如下式，且必须满足必要的边界条件。

其中C₀、C_i是任意实数，N_i为简单函数，例如低阶多项式。在伽辽金方法中，直接采用试函数自身作为权函数，即

v₁(x,y,z)＝N_i(x,y,z) (12)

后面的有限元过程考虑稳态热传导，不考虑控制方程中的时间项，所以三维热传导方程弱形式写成

C.采用四面体网格剖分求解域；

采用四面体网格剖分求解域，剖分后的求解域被分割为三维四面体网格，从而将连续的几何结构空间转化为离散的网格空间。

D.选择叠层基函数，离散B中得到的有限元弱形式，得到热辐射边界条件的有限元方程；

如图2所示四面体单元中i,j,k,l代表四个顶点的编号，我们首先得到四个最基本的基函数：

式中

将(18)式、(19)式、(20)式和(21)式中的i,j,k,l轮换，得到a_j,a_k,a_l，b_j,b_k,b_l，c_j,c_k,c_l，d_j,d_k,d_l。V为四面体的体积。

标量二阶叠层基函数，选择N₁,N₂,N₃,N₁N₂,N₁N₃,N₂N₃这六个基函数。对于有限元过程来说，把域Ω离散为M个单元之后，如同(13)式所示的弱形式定积分，可以通过将每个单元的积分贡献简单相加，即

对于每一个单元来说，跟有限元方程组右端项有关系的

和

的求解在很多有限元基础材料中都有介绍，在此不再赘述。本发明的核心是用Newton-Raphson迭代方法来进行非线性辐射边界的处理，因此重点考虑的是有限元方程组左端矩阵的求解。

假设非线性方程组的形式为：

SA＝F (24)

其中S为左端矩阵，A为待求解向量，F为右端向量。

对于Newton-Raphson法是一种梯度算法，将(24)式的有限元方程组写成如下形式：

f(A)＝SA-F＝0 (25)

其中f(A)为非线性函数。应用Newton-Raphson方法则有

A^(q+1)＝A^(q)-[f′(A^(q))]^-1f(A^(q)) (26)

其中上标q表示迭代次数，令J^(q)＝f′(A^(q)),并化简有

J^(q)ΔA^(q)＝F-S^(q)A^(q) (27)

其中J^(q)为非线性函数f(A)对A的导数矩阵[称为雅可比(Jacobi)矩阵]的第q次迭代值，ΔA^(q)前后两次迭代的温度差值，满足

ΔA^(q)＝A^(q+1)-A^(q) (28)

其中A^(q+1)为后一次求解的值，A^(q)为前一次求解的值。通过求解线性方程组的方法求解方程组(27)式，可以得到ΔA^(q)，然后根据(28)式可以获得A^(q+1)，其流程图如图3所示。

Newton-Raphson迭代方法需要计算雅可比矩阵J,对于有限元计算来说，雅可比矩阵J可以由各单元的雅可比矩阵J^e叠加组成。由(23)式，定义单元矩阵S_ij：

令S_ij＝K_ij+M_ij，其中

应用Newton-Raphson法可以求得

K_ij矩阵的求解在各类基础有限元书籍中都有介绍，这里不作具体论述,M_ij矩阵含有待求解未知项u，计算比较困难，所以在此考虑高斯数值积分

其中F为函数表达式，S_Δ表示三角形面积，W_k表示权重，m表示采样点的个数，(N_1k,N_2k,N_3k)表示采样点。

接下来介绍M矩阵的求解，令

u＝u₁N_1k+u₂N_2k+u₃N_3k+u₄N_1kN_2k+u₅N_1kN_3k+u₆N_2kN_3k (34)

其中u₁,u₂,u₃,u₄,u₅,u₆表示求得的解向量对应的小三角形辐射面上的顶点和中点的值。当i＜4&&j＜4时有：

其中i,j,k为下标。

当i＜4&&j≥4时有：

其中m,n为下标，并且和j满足图4表格所示的对应关系。

当i≥4&&j≥4时有：

其中a,b为下标，并且和i满足图4表格所示的对应关系，m,n为下标，和j同样满足图4表格所示的对应关系。

接下来介绍J矩阵的求解，用数值积分的方法对(32)式进行展开。

当i＜4&&j＜4时有：

其中u,v,p,q为下标,并且u,v和p满足图4表格的对应关系。

当i＜4&&j≥4时有：

其中u,v,q为下标,并且u,v和p满足图4表格的对应关系；m,n,j为下标,并且m,n和j也满足图4表格的对应关系；

当i≥4&&j≥4时有：

其中u,v,q为下标,并且u,v和p满足图4表格的对应关系；m,n,j为下标,并且m,n和j也满足图4表格的对应关系；a,b,i为下标,并且a,b和i同样满足图4表格的对应关系；E.对步骤D中的有限元非线性方程组不断地进行迭代，直到其温度值满足一定的收敛规则。

在步骤D中已经完成所需矩阵J和S的求解，接下来只需要按照Newton-Raphson迭代方法所述，完成非线性方程组的构建。对于非线性方程组的迭代求解过程，本发明采用的收敛规则是||F-S^(q+1)A^(q+1)||＜ε或者

(ε是人为设定的收敛精度值)，前者的收敛判据更加准确，但是处理过程可能相对后者更加复杂，后者判据则相对简单，但精度可能会有所降低。

综上所述，本发明针对目前热辐射边界条件在有限元方法中的应用难题，提出了一套通用的有限元解决方法，采用高阶叠层基函数、Newton-Raphson非线性迭代方法以及高斯数值积分的使用，可以准确、快速地求解高度非线性热辐射边界问题，并且解决了现有技术关于这方面描述不足以及一些特定解决方法的适用性不足的难题。

Claims (2)

1.一种基于热辐射边界条件的三维有限元模拟方法，应用于空间行波管，包括以下步骤：

A.对目标器件进行建模，建立对应的几何结构模型；

B.采用伽辽金加权残数法得到热辐射边界条件的有限元弱形式；

C.采用四面体网格剖分求解域；

D.选择基函数，离散步骤B中得到的有限元弱形式，得到热辐射边界条件的有限元方程；

E.对步骤D中的有限元非线性方程组不断地进行迭代，直到其温度值满足收敛规则,收敛规则是||F-S^(q+1)A^(q+1)||<ε或者

其中ε为设定的收敛精度值，ΔA^(q)是前后两次迭代的温度差值，F是右端向量，q是迭代次数，S^(q)是采用Newton-Raphson方法之前的有限元初始矩阵项，A^(q)是前一次迭代的温度值；

所述步骤D中的基函数为高阶叠层基函数，迭代方法选用Newton-Raphson非线性迭代方法。

2.如权利要求1所述基于热辐射边界条件的三维有限元模拟方法，其特征在于：所述步骤D中基函数选用二阶叠层基函数。

Images