CN112836399B - 一种基于有限元算法的非线性接触热阻热分析求解方法 - Google Patents

Abstract

本发明属于三维热传导有限元求解技术领域，具体涉及一种基于有限元算法的非线性接触热阻热分析求解方法。本发明针对目前热分析领域无法处理非线性接触热阻的弊端，将接触热阻问题转化为等效的边值问题，并提出一套新型非线性求解器，使用基于能量函数的牛顿‑拉夫逊方法对得到的有限元方程进行求解，通过在每步NR迭代中寻找一个最优的

的近似值

使得牛顿‑拉夫逊方法的总体迭代次数减少；在不改变原有模型的基础上，高效、规范、精确地实现了非线性接触热阻热分析求解。

Description

一种基于有限元算法的非线性接触热阻热分析求解方法

技术领域

本发明属于三维热传导有限元求解技术领域，具体涉及一种基于有限元算法的非线性接触热阻热分析求解方法。

背景技术

对于两个名义上互相接触的固体表面，实际上固体和固体的直接接触只能发生在一些离散点或微小的面积上，由于间隙介质的导热系数与固体导热系数一般相差很大，因而引起接触面附近热流改变，形成的热流附加阻力，即接触热阻。在航天、机械制造、微电子等领域内，各部件之间接触热阻是热力响应的关键参数，在有些情况下，考虑接触热阻和不考虑接触热阻，热分析的结果甚至可以相差50％。如果无法精确计算出接触热阻带来的影响，对物体进行热分析会产生一定的误差，甚至计算出完全错误的结果。

目前所有的热有限元代码处理接触热阻的方式主要有两种。一种方式是接触薄层法：该方法在两个组件的接触处创建一个厚度非常薄(通常微米量级)且和接触面共形的薄层(避免破坏原来几何结构)；另一种方法是ANSYS、CST等商业有限元代码最新版本中采用的接触边界法：和接触薄层法不同，接触边界法无需构建接触薄层，而只需在接触面上设置接触热阻，然后将接触面作为边界条件施加到有限元热分析中。

当下市面上所有热有限元代码中接触边界法的接触热阻值只能设为一固定值，认为部件工作状态下接触热阻的值是恒定的，但通过前人的研究发现，接触热阻主要是由热、力、材料三种因素耦合而成的，除此之外还受到表面粗糙度、界面载荷、接触体的材料特性等因素的影响。对于航空航天、核能、微电子等领域，由于部件的工作环境的复杂性，在热分析中只考虑线性接触热阻的影响，最终得到的结果可能会产生较大的误差。

发明内容

针对上述存在问题或不足，为解决现有有限元代码的接触边界法在实际中效率低下，误差难以控制，未考虑接触热阻的非线性特征等问题，本发明提供了一种基于有限元算法的非线性接触热阻热分析求解方法，该方法将接触热阻问题转化为有限元边界条件进行处理，高效便捷的解决了现有的技术难题。

一种基于有限元算法的非线性接触热阻热分析求解方法，包括以下步骤：

S1.对欲进行热分析的对象建立对应的几何结构模型。

S2.采用四面体网格划分策略，对S1得到的几何结构模型进行网格划分，获得网格数据。

S3.在几何结构模型的物理接触面上形成数值接触面，将物理接触面上设置的接触热阻转化为边界条件，采用伽辽金方法获得热分析的有限元弱形式。

根据接触热阻的定义，在虚拟数值接触面上可得如下边界条件：

采用伽辽金方法得到热分析的有限元的弱形式：

上式中W为测试函数；L表示慢波结构中物理接触面的总数；上标l表示该变量为第l个接触面中的物理量。

S4.使用叠层基函数对S3获得的热分析有限元弱形式进行离散，得到最终待求解有限元矩阵与右端项。

S5.使用基于能量函数的牛顿-拉夫逊(Newton-Raphson)NR方法对得到的有限元方程进行求解，并通过减少每一步迭代的迭代时间和减少总体迭代次数，加快非线性求解器的求解速度，计算出最终结果；

在本发明中，S4所得的最终待求解矩阵是一个复杂的高度非线性矩阵方程，为了进一步提高求解的速度与精度，这里提出一种新的非线性求解器。

常规的牛顿-拉夫逊方法是基于以下迭代方案获得收敛的解：

{x}^k+1＝{x}^k+α^k{Δx}^k where {Δx}^k＝-([J]^k)^-1{r}^k (12)

对于上述方法，为了加快非线性求解器的求解速度，可以通过两种方式进行处理：1)减少每一步迭代的迭代时间，2)减少总体迭代次数。

对于减少总体迭代次数，由于常规的牛顿-拉夫逊方法(a^k固定为1)只具有局部收敛性，非线性热有限元分析不仅需要比较精确的初始解向量{x}⁰，而且其高度非线性会使得牛顿-拉夫逊方法迭代步数很长，甚至无法收敛。本方法提出一种泛函最小化技术，通过在每步NR迭代中寻找一个最优的

的近似值

以解决这个困难，使得牛顿-拉夫逊方法的迭代次数大大减少。

具体为：在有限元分析中，解向量{x}需要使得能量泛函F最小化，而在每步NR迭代中，解向量{x}^k+1可由a^k表示，因此最优的

同样应该要使泛函F^k+1最小化。因而在每步NR迭代解得{Δx}^k后，计算泛函F^k+1关于a^k的偏导数。利用有限元方法的Ritz方程以及(12)式，可将该偏导数简化，并且不用F^k+1显式表示，即：

上式上标T表示向量转置。最优的a_ko使得该偏导数为零，即有：

由于F^k+1可近似看成a^k的二次函数，因此

可近似为线性函数。利用这一性质，我们只需通过计算两个固定的a^k值(如0.5和1.0)下的

值，从而获得

近似的线性方程；直接令该线性方程为零就可以解得最优

的近似值

进一步的，所述步骤S5还包括减少每一步迭代的迭代时间，即加快{Δx}^k的求解速度，由于[J]^k为大规模不对称矩阵，我们的非线性求解器采用不对称迭代法即GCR方法求解{Δx}^k并使用三阶p型多重网格预处理技术和非对称ILU分解技术对{Δx}^k的求解进行进一步加速。

本发明通过使用接触边界法，在不构建接触薄层的前提下，通过在接触面上设置接触热阻，将接触面作为边界条件施加到有限元热分析中。并进一步考虑到材料非线性问题与接触热阻非线性问题，通过对接触面上网格单独离散处理，精确模拟接触热阻存在所带来的温度跳变，使得热分析结果更加符合真实情况。最后为了进一步提高热分析矩阵的求解的速度与精度，本发明提出了一种新型的非线性求解器，通过在每步NR迭代中寻找一个最优的

的近似值

使得牛顿-拉夫逊方法的总体迭代次数减少，大大加快了复杂的高度非线性矩阵方程的求解速度。针对目前热分析领域无法处理非线性接触热阻的弊端，本发明采用以上种种技术，将接触热阻问题转化为等效的边值问题，在不改变原有模型的基础上，高效、规范、精确地实现了非线性接触热阻热分析求解。

附图说明

图1是本发明的流程示意图。

图2是本发明实施例网格划分示意图。

图3是本发明实施例物理接触面转化数值接触面的示意图。

图4是实施例接触热阻作为边界条件进行处理的示意图。

图5是实施例获得的温度场分布示意图。

具体实施方式

下面结合附图和具体实施例对本发明作进一步说明。

本发明的流程如图1所示，包括如下步骤：

S1.对欲进行热分析的对象建立对应的几何结构模型。

S2.采用四面体网格划分策略，对S1得到的几何结构模型进行网格划分，获得网格数据。

在步骤S1得到的几何结构模型上，使用四面体划分方法，对有限元求解区域内的空间进行对应的网格划分。

S3.几何结构模型物理接触面上形成数值接触面，将物理接触面上接触热阻转化为边界条件，采用伽辽金方法获得热分析的有限元弱形式。

在使用有限元方法对三维模型进行热分析时，模型内温度T分布由以下边值问题控制：

第一个方程表示在计算域Ω内的稳态导热方程，后续三个方程分别表示热流边界Γ_H，对流冷却边界Γ_C和辐射冷却边界Γ_R上的边界条件方程，其中λ表示热导率，

表示拉普拉斯算子，n表示边界面上的传出法向量，h表示对流传热系数，T表示待求解温度，T_a表示环境温度，σ表示斯蒂芬-波尔兹曼(Stefan-Boltzmann)常数，ε表示辐射边界的辐射系数。根据接触热阻的定义，在虚拟数值接触面上可得如下边界条件：

上式中T₁、T₂式表示接触面两端相邻的区域1、区域2的温度，n₂、n₁表示接触面Γ₁、Γ₂的外向法向量，σ(T₁,T₂)表示接触面上的接触热阻值，λ₁(T₁)、λ₂(T₂)表示区域1与区域2随温度变化的热导率。

本方法直接采用其他研究者提出的热分析模型常用材料热导率精确计算公式进行计算：

λ(T)＝λ₀(1+β₁T+β₂T²+...+β_pT^p) (4)

由上式可知模型材料热导率随温度变化明显，导致接触热阻值也随接触面两边组件的温度变化明显，因此σ可记为σ(T₁,T₂)，使用的接触热阻计算公式如下：

上式中sr₁、sr₂表示接触面两端材料的表面粗糙度，H表示接触面两端材料微硬度的小值，P_F为接触面上的接触压强。部件采用不同装配工艺其P_F值不同，需要事先通过试验测得。一旦P_F值测得，当类似结构采用相同的装配工艺时，即使采用不同的材料也可认为接触处的压强P_F保持不变。和热导率λ相比，sr、H和P_F随温度的变化很小，可以忽略，因此在(5)中它们对于温度是常数。

将模型中所有物理接触面都按照上述过程进行操作，则在每对数值接触面上都可得到类似的边界条件。将这些数值接触面上的接触边界条件添加到热分析边值问题中就构成能够处理接触热阻的热分析边值问题。

通过使用伽辽金残数法，能很容易地将(1)-(3)式结合，得到对应的非线性热分析的有限元的弱形式：

上式中W表示伽辽金残数法的测试函数；L表示热分析模型中物理接触面的总数；上标l表示该变量为第l个接触面中的物理量。

S4.使用叠层基函数对S3获得的热分析有限元弱形式进行离散，得到最终待求解有限元矩阵与右端项。

考虑到四面体单元具有良好的边界拟合性，本方法选择基于四面体单元的高阶叠层标量基函数进行有限元弱形式的空间离散。在离散过程中，在每个物理接触面上无需做任何处理，就能生成一套三角形网格。然后将每个物理接触面上网格都虚拟成完全一样的两套网格分别作为数值接触面

和

上的网格。离散后，(6)式中的T,

用一系列基函数按照如下形式展开：

上式中N_i表示四面体单元的第i个基函数，x_i表示N_i对应的插值系数，D表示四面体单元内基函数的总数,

表示接触面上的温度基函数形式。

由于在数值接触面

和

上有各自的网格，因此

和

上的基函数也不同，其中第i个基函数分别记为

和

和

为其对应的展开系数，D_l为

或

的总数。需用特别指出的是，正是这一重要的性质使得

和

可以用不同基函数和系数展开(如(7)式所示)以保证

和

取不同的值，从而能精确模拟接触热阻存在所带来的温度跳变。在本方法中，N_i、

和

的构造方式都采用三阶标量叠层方式构造。最后，将每一个N_i、

和

分别作为测试函数W、

和

代入到(6)式中，就可以将热分析转换成求解如下的矩阵方程：

其中上式中各个矩阵和右端项的具体项由下给出：

由于数值接触面

和

上的每一个三角形网格都属于一个四面体单元，因此

和

是解向量{x}的一部分。

S5.使用基于能量函数的牛顿-拉夫逊(Newton-Raphson)方法对得到的有限元方程进行求解，计算出最终结果。

对于S4得到待求解有限元矩阵，由于热导率λ(T)和接触热导σ(T₁,T₂)都是温度函数，因此都使得矩阵方程(8)为非线性方程。当两者同时存在时，(8)式中不仅体积分矩阵[S]为非线性，而且面积分矩阵如

等也为非线性，因此(8)式为一个复杂的高度非线性矩阵方程。为了精确快速求解(8)所得到的非线性方程，本方法采用一种基于能量函数的牛顿-拉夫逊方法。

常规的牛顿-拉夫逊方法是基于以下迭代方案获得收敛的解：

{x}^k+1＝{x}^k+α^k{Δx}^k where {Δx}^k＝-([J]^k)^-1{r}^k (12)

其中[J]^k和{r}^k分别代表第k步迭代过程中的雅可比矩阵与残差，对于本方法而言，根据(8)-(11)，具体的雅克比矩阵和残差可由下式得到：

在(12)式中{Δx}^k的计算是通过求解线性矩阵方程得到:

[J]^k{Δx}^k＝-{r}^k (17)

由于[J]^k为大规模不对称矩阵，我们的非线性求解器采用不对称迭代法即GCR方法求解线性矩阵方程(17)，并使用三阶p型多重网格预处理技术和非对称ILU分解技术对上式的求解进行进一步加速，以上技术属于本领域公知，这里不再赘述。

然而，由于常规的牛顿-拉夫逊方法(a^k固定为1)只具有局部收敛性，非线性热有限元分析不仅需要比较精确的初始解向量{x}⁰，而且其高度非线性会使得牛顿-拉夫逊方法迭代步数很长，甚至无法收敛。本方法通过在每步NR迭代中寻找一个最优的

以解决这个困难。

为了找到最优的

本方法提出一种泛函最小化技术，该技术简单高效，其核心思想是：在有限元分析中，解向量{x}是使得能量泛函F最小化，而从(12)式可知，在每步牛顿-拉夫逊方法(NR)迭代中，解向量{x}^k+1可由a^k表示，因此最优的

同样应该要使泛函F^k+1最小化。因而在每步NR迭代解得{Δx}^k后，计算泛函F^k+1关于a^k的偏导数。利用有限元方法的Ritz方程以及(12)式，可将该偏导数简化，并且不用F^k+1显式表示，即：

上式上标T表示向量转置。最优的

应该使得该偏导数为零，即有：

由于F^k+1可近似看成a^k的二次函数，因此

可近似为线性函数。利用这一性质，我们只需通过(18)式计算两个固定的a^k值(如0.5和1.0)下的

值，从而获得

近似的线性方程。直接令该线性方程为零就可以解得最优

的近似值

将

作为第k步NR迭代的a^k值代入到(12)式中，不仅无需精确的初始解向量{x}⁰(本方法中{x}⁰简单设置为零向量)，而且能显著减少(8)式的NR迭代次数，从而大幅提高非线性热分析效率。

综上所述，本方法通过使用基于能量函数的牛顿-拉夫逊方法对得到的有限元方程进行求解，通过在每步NR迭代中寻找一个最优的

的近似值

使得牛顿-拉夫逊方法的总体迭代次数减少，能精确快速地计算出最终结果，得到热分析部件稳态下的温度场分布。本发明针对目前热分析领域无法处理非线性接触热阻的弊端，将接触热阻问题转化为等效的边值问题，并提出一套新型非线性求解器，在不改变原有模型的基础上，高效、规范、精确地实现了非线性接触热阻热分析求解。

Claims (2)

1.一种基于有限元算法的非线性接触热阻热分析求解方法，其特征在于，包括以下步骤：

S1.对欲进行热分析的对象建立对应的几何结构模型；

S2.采用四面体网格划分策略，对S1得到的几何结构模型进行网格划分，获得网格数据；

S3.在几何结构模型的物理接触面上形成数值接触面，将物理接触面上设置的接触热阻转化为边界条件，采用伽辽金方法获得热分析的有限元弱形式；

S4.使用叠层基函数对S3获得的热分析有限元弱形式进行离散，得到最终待求解有限元矩阵与右端项；

S5.使用基于能量函数的牛顿-拉夫逊NR方法对得到的有限元方程进行求解，并通过减少每一步迭代的迭代时间和减少总体迭代次数，加快非线性求解器的求解速度，计算出最终结果；

牛顿-拉夫逊方法基于以下迭代方案获得收敛的解：

{x}^k+1＝{x}^k+α^k{Δx}^k where{Δx}^k＝-([J]^k)^-1{r}^k (12)其中[J]^k和{r}^k分别代表第k步迭代过程中的雅可比矩阵与残差；

所述减少总体迭代次数的方法具体为：

通过在每步NR迭代中寻找一个最优的

的近似值

使得牛顿-拉夫逊方法的总体迭代次数减少；

在有限元分析中，解向量{x}需要使得能量泛函F最小化，而在每步NR迭代中，解向量{x}^k+1可由a^k表示，因此最优的

同样应该要使泛函F^k+1最小化；因而在每步NR迭代解得{Δx}^k后，计算泛函F^k+1关于a^k的偏导数，利用有限元方法的Ritz方程以及(12)式，将该偏导数简化，并且不用F^k+1显式表示，即：

上式上标T表示向量转置，最优的

应该使得该偏导数为零，即有：

由于F^k+1可近似看成a^k的二次函数，因此

可近似为线性函数，利用这一性质，只需通过(18)式计算两个固定的a^k值下的

值，从而获得

近似的线性方程；直接令该线性方程为零以解得最优

的近似值

2.如权利要求1所述基于有限元算法的非线性接触热阻热分析求解方法，其特征在于：所述步骤S5中非线性求解器采用不对称迭代法即GCR方法求解{Δx}^k并使用三阶p型多重网格预处理技术和非对称ILU分解技术进行加速，以减少每一步迭代的迭代时间。

Images