时域热传导仿真方法及存储介质与流程

1.本发明涉及热学仿真技术领域，特别是涉及一种基于高阶四面体单元的时域热传导仿真方法，以及一种用于执行所述时域热传导仿真方法中步骤的计算机可读存储介质。


背景技术：

2.20世纪40年代早期，a.hrennikoff和r.courant针对建筑结构和飞行器结构的设计的力学分析提出发展了现在称为有限元法(finite element method)的数值方法。
3.1965年，winslow将fem应用于电气工程分析中。1966年，kane yee首先提出了基于时域有限差分法(finite difference time domain,fdtd)的电磁学时域仿真方法。
4.1969年，silvester首先将fem运用到高频领域，应用于微波波导设计。
5.1997年，jin-fa lee等人提出并详细介绍了基于时域有限元法(finite element time domain,fetd)的电磁学时域数值仿真方法。
6.2016年，jian-ming jin等人基于fem，fetd实现了电热耦合仿真，利用feti提升计算效率，应用于射频器件的热分析。
7.2019年，zhang hao-xuan等人采用并行计算技术优化射频元器件多物理仿真的效率。
8.2020年，yan su，jian-ming jin等人发表“an enhanced transient solver with dynamic p-adaptation and multirate time integration for electromagnetic and multiphysics simulations”。基于时域离散伽辽金法(discontinuous galerkin time domain，dgtd)开发了多物理(电磁学与等离子体)的仿真解算器，针对多尺度问题(multi-scale problem)利用动态p优化适应和多时间步的方法提升计算效率。
9.关于有限元的算法研究已经比较完善，但是，对于基于高阶四面体单元的热学时域有限元求解方法还比较少。部分研究在进行电热多物理仿真时，只采用了单项耦合的模型，计算了电磁场的热效应而忽略了温度对于电磁场的影响。目前，对各个物理场都采用时域分析(即全时域)的电热双向耦合仿真方法的研究是很少的。因此，基于高阶四面体单元的热学时域有限元求解方法对于多物理问题有很大的作用。


技术实现要素：

10.在发明内容部分中引入了一系列简化形式的概念,该简化形式的概念均为本领域现有技术简化，这将在具体实施方式部分中进一步详细说明。本发明的发明内容部分并不意味着要试图限定出所要求保护的技术方案的关键特征和必要技术特征，更不意味着试图确定所要求保护的技术方案的保护范围。
11.为了克服上述现有技术的不足，本发明提供了一种基于高阶四面体单元的时域热传导仿真方法。
12.为解决上述技术问题，本发明提供的时域热传导仿真方法，包括以下步骤：
13.s1，利用有限元网格剖分工具几何建模和网格剖分，读取几何信息；
14.s2，设置固体传热学分析中各个区域材料相关参数；
15.s3，构建热学时域有限元的基本方程计算出单元矩阵中每一个元素的值，热学采用基于四面体单元的线性点基标量基函数，获得每个单元的单元矩阵；
16.s3，通过每个单元、节点、边的编号，将单元矩阵中的元素与全局矩阵的元素对应，通过遍历叠加利用单元矩阵，得到全局矩阵；
17.s4，通过热学时域有限元的基本方程对每一个时间步进行迭代运算。
18.可选择的，所述的时域热传导仿真方法，还包括：
19.s5，通过验证工具对热学时域有限元的基本方程结果进行验证。
20.可选择的，所述的时域热传导仿真方法，步骤s1中，通过有限元网格剖分工具(例如trelis)进行几何建模和网格剖分，通过第三方仿真软件(例如matlab)平台读取几何信息。
21.可选择的，所述的时域热传导仿真方法，步骤s2中，各个区域材料相关参数包括：热导率，质量密度，比热容、边界条件以及观测数据所需的域点探针。
22.可选择的，所述的时域热传导仿真方法，热学时域有限元的基本方程表达如下；
[0023][0024]
t为温度，c为材料比热容，t为时间，[c
t
]为阻尼矩阵，[k
t
]为刚度矩阵，{t}为离散形式的温度场标量。
[0025]
可选择的，所述的时域热传导仿真方法，步骤s3中，每个单元的单元矩阵包括质量矩阵、阻尼矩阵和刚度矩阵。
[0026]
可选择的，所述的时域热传导仿真方法，步骤s4中，进行迭代运算采用lupq矩阵分解法。
[0027]
可选择的，所述的时域热传导仿真方法，步骤s5中，通过构建算例采用第三方仿真软件(例如comsol multiphysics)作为验证工具。
[0028]
本发明还提供了一种用于执行上述任意一项所述的时域热传导仿真方法中步骤的计算机可读存储介质。
[0029]
传统的电子器件设计环节中，仅对与电子器件工作相关的电磁场进行分析，经常忽略了力学热学等其他物理场对于器件性能的影响。21世纪以来，电子产品向着大功率，高性能，小型化的方向发展，电子元件、设备的发热问题越来越不可忽视。因此温度场的分析于设计在工业中愈发重要。任何让工程师更好地分析不同物理场之间的联系，成为了高性能产品工业设计中一个十分重要的问题。科学研究上，随着跨领域交叉学科的发展，对于热学时域有限元研究的需求也逐渐增大。
[0030]
本发明利用有限元网格剖分工具几何建模和网格剖分，读取几何信息；基于高阶四面体单元构建热学时域有限元的基本方程计算获得全局矩阵；再通过热学时域有限元的基本方程对每一个时间步进行迭代运算；通过验证工具对热学时域有限元的基本方程结果进行验证进而能获得准确的时域热传导仿真结果。可靠高效的时域热传导方程仿真方法对于高性能电子产品综合设计与科学理论问题的研究验证有着十分重要的意义。
附图说明
[0031]
本发明附图旨在示出根据本发明的特定示例性实施例中所使用的方法、结构和/或材料的一般特性，对说明书中的描述进行补充。然而，本发明附图是未按比例绘制的示意图，因而可能未能够准确反映任何所给出的实施例的精确结构或性能特点，本发明附图不应当被解释为限定或限制由根据本发明的示例性实施例所涵盖的数值或属性的范围。下面结合附图与具体实施方式对本发明作进一步详细的说明：
[0032]
图1是对待分析对象进行网格剖分示意图。
[0033]
图2是fetd热学求解与comsol multiphysics结果对比示意图。
[0034]
图3是fetd热学求解与comsol multiphysics对比相对误差示意图。
[0035]
图4搭载散热片的体热源示意图。
具体实施方式
[0036]
以下通过特定的具体实施例说明本发明的实施方式，本领域技术人员可由本说明书所公开的内容充分地了解本发明的其他优点与技术效果。本发明还可以通过不同的具体实施方式加以实施或应用，本说明书中的各项细节也可以基于不同观点加以应用，在没有背离发明总的设计思路下进行各种修饰或改变。需说明的是，在不冲突的情况下，以下实施例及实施例中的特征可以相互组合。本发明下述示例性实施例可以多种不同的形式来实施，并且不应当被解释为只限于这里所阐述的具体实施例。应当理解的是，提供这些实施例是为了使得本发明的公开彻底且完整，并且将这些示例性具体实施例的技术方案充分传达给本领域技术人员。
[0037]
本发明所涉及的物理场为热学场(温度分布)，对应的对其进行描述的方程为热传导方程。
[0038]
热传导方程描述了固体物质中的热传递现象。
[0039][0040]
其中t为温度，ρ为质量密度，c为材料比热容，λ为热传导系数，q为热功率体密度。
[0041]
在热分析中，采用对流边界条件(conviction boundary)：
[0042][0043]
其中t为温度，tf为流体环境温度，λ为热传导系数，为边界面法向量。
[0044]
数值方法是采用计算机求解数学问题近似解的问题。对于不同的问题，已经提出了很多相关的数值算法。常用于求解偏微分方程组的算法有矩量法，有限差分法以及本发明使用的有限元法。
[0045]
本发明提供一种时域热传导仿真方法，包括以下步骤：
[0046]
s1，利用有限元网格剖分工具几何建模和网格剖分，读取几何信息；
[0047]
示例性的，几何建模：通过编写.jou脚本结合gui在cubit trelis中进行几何建模。trelis是一专用的有限元网格剖分工具，其中内置了几何建模工具，并允许用户对构建的几何模型设置标签信息(如材料编号，几何体编号等)。cubit trelis的几何建模支持从一维到三维的几何结构。遵循基于专业工具的方法，cubit trelis提供了各种不同的网格
划分技术，协作基础架构和算法。cubit trelis支持.jou脚本或python脚本进行几何建模，在批量生成算例时，可以使用脚本建模提升速度。
[0048]
网格剖分：通过编写.jou脚本结合gui在cubit trelis中实现网格剖分。以.gfile文件导出。处理三维模型时cubit trelis的网格剖分支持六面体单元，四面体单元等常用单元类型，在批量生成算例时，可以使用脚本建模提升速度。基函数。在有限元法中，基函数是局部的。有限元法通过网格剖分将大型的几何对象剖分为更小，更简单的单元，基于这些单元，基函数得以构建。
[0049]
导入几何信息：将.g file转换为文本文件，再在matlab平台中编写脚本读取文件中包含的几何信息。
[0050]
s2，设置固体传热学分析中各个区域材料相关参数；
[0051]
示例性的，将读取的信息导入matlab workspace后，可以根据需要进行后处理，一般的做法采用单独的.m对参数信息，即固体传热学分析中各个区域材料的热导率，质量密度，比热容以及边界条件等参数进行设置。此外，对于观测数据所需的域点探针，需要在此步骤设置。
[0052]
s3，构建热学时域有限元的基本方程计算出单元矩阵中每一个元素的值，热学采用基于四面体单元的线性点基标量基函数，获得每个单元的单元矩阵；
[0053]
有限元单元可根据基础单元的几何形状分类。针对三维有限元仿真，常见的形状有：长方体，四面体，五面体，曲面六面体等。
[0054]
长方体单元是比较简单的单元。选用长方体单元时，剖分算法较为简单，对于比较简单的对象甚至可以人工进行划分，获取几何信息时较为直接方便。但是受到阶梯效应的影响，长方体单元在模拟几何结构时性能较低,因此在工业中比较少直接使用简单长方体作为有限元的基本单元。
[0055]
四面体单元单元单元相对于长方体单元而言，在模拟几何外形时更加灵活，产生的误差较小，因此广泛的应用于热学，电磁学，力学有限元分析中。
[0056]
五面体单元相较于四面体而言没有明显的精度上的提升，而每个单元的顶点数和棱边数更多，一定程度上提高了计算量。此外，五面体单元的适用范围不如四面体单元广泛，较难找到合适的参考文献，这一定程度上增大了编写程序实现有限元方针的难度。
[0057]
曲面六面体的几何模拟能力比四面体单元更强，但是曲面六面体单元不能通过顶点简单的描述，这给几何信息的处理带来很大的麻烦。对于具有弯曲边界面的单元，在定义其形状时，需要额外的辅助点来保证单元足够准确地拟合了仿真对象的局部几何形状。选用曲面六面体单元实现有限元，在进行网格剖分以及空间离散化时过程比较复杂，实现难度较高。
[0058]
本发明选取性能与可行性皆较高的四面体单元作为基本单元。
[0059]
有限元的基函数，以函数变量进行区分，可分为矢量基函数与标量基函数两大类。热学分析的场量主要为温度，属于标量，所以实现热学有限元时采用标量基函数进行分析。基函数可以根据有限元单元的顶点或棱边进行定义，对于矢量基函数，采用棱边定义是比较方便的；而对于标量基函数，采用顶点定义比较方便。因此，热学有限元采用点基标量基函数。
[0060]
为了求解得到更加精确的结果，可以应用更高次的单元对求解区域进行剖分。图
中显示了各阶四面体单元。一阶四面体有四个结点，分别在各个顶点处。二阶四面体单元在四面体的六条棱边内各增加了一个结点，三阶四面体在二阶四面体的基础上又在六条棱边内新增加了六个结点，并在四面体的四个面内也增加了一个结点。高阶单元将会包含更多的结点信息，因此采用高阶单元剖分求解区域，有限元计算精度会极大改善。出于计算难度和复杂度的考虑，本发明选用二阶四面体结构来对求解区域进行网格剖分。
[0061]
使用高阶四面体单元对求解域进行网格剖分后，接下来要选择一组具有某些性质的试探函数作为单元中未知解的插值函数。插值函数的选取将直接关系到有限元计算结果的准确性。对于含有n个结点的单元e中，无论该单元是多少阶的，每一个结点将对应一个插值函数。虽然高阶单元的插值函数会比低阶单元的插值公式更复杂，但它的计算精度也会大大提高。
[0062]
对于二阶四面体单元，插值函数分为角结点和棱内结点两种：
[0063]
角结点：ni＝(2l
i-1)li(i＝1，2，3，4)
[0064]
棱内结点：n5＝4l1l2[0065]
n6＝4l1l3[0066]
有限元法中几何和物理场的空间离散化是基于许多微小单元的，故被称作有限元法。在早期，有限元的理论发展是基于变分方法的，在最开始使用有限元法计算仿真时，使用的仿真多用变分法推导其基本方程。因此它在物理领域的应用场景多为由laplace方程和poisson方程描述的场景。有限元法也因此常被用于与功能极值问题密切相关的问题。在之后的研究中，发现可以使用加权余量法(weighted residual method，wrm)中的伽辽金法(galerkin’s method)也可以得出相同形式的有限元基本方程。这一理论上的发展很大程度拓宽了有限元法的应用场景。理论上，任何偏微分方程组描述的物理场都可以用有限元法进行仿真计算。因此，近年来越来越多关于有限元法的研究都改使用伽辽金法作为推导有限元法方程的基础，并把有限元法归为加权余量法的一个分支。在本毕业论文中的理论推导部分，也采用伽辽金法推导。
[0067]
有限元法能用于仿真包括热学、光学、电磁学、流体力学、固体力学在内的等多种物理场，因此在工业界和学术界的各个领域得到了广泛的应用。因为适用于大多数物理场景，有限元也是比较适合用于仿真多物理场耦合问题的算法。
[0068]
有限元是通过在空间尺寸上进行特定的空间离散化来实现数值计算的。具体的过程包括几何网格剖分，基于剖分的与选定基函数的空间离散化与最终的计算(对于时域有限元法而言是迭代计算)三个步骤。其中几何网格剖分为空间离散化提供了基础的几何信息，而空间离散化将原本复杂的问题转换为一个线性方程组以便最终的计算。
[0069]
时域有限元则是有限元法的时域版本，通过某种时间步进方法进行迭代计算，从而得到整个观测时间域上场量随时间变化情况的方法，常见的时间步进方法有前/后向差分法和龙格-库塔法。
[0070]
之前已经给出，热传导分析的域内控制方程为：
[0071][0072]
之前已经给出，热传导分析的边界内控制方程为：
[0073]
[0074]
根据之前给出的加权余量法，可以推导出热传导分析的矩阵形式有限元方程，这里直接给出：
[0075][0076]
上式中，[c
t
]为阻尼矩阵，[k
t
]为刚度矩阵，{t}为离散形式的温度场标量。上式中的矩阵以及列向量中的元素可由以下公式给出：
[0077]
[c
t
]
ij
＝∫vρcninjdv
[0078][0079]
{p
t
}i＝∫vqnidv+∫shnitads
[0080]
式中n是有限元中使用的基函数，由于温度是标量，在热学分析中这些基函数采用点基(nodal based)标量基函数。热学的时间离散方法采用后向差分法(backward difference method，bdm)[19]。针对温度-时间的变化关系，后向差分法方程如下：
[0081][0082]
dt
[0083]
无论的值如何选取，采用后向差分法构建的系统是无条件稳定的。将后向差分法带入原式，可以得到热学时域有限元的最终公式：
[0084][0085]
计算出单元矩阵中每一个元素的值，热学采用基于四面体单元的线性点基标量基函数。获取方程中的几个矩阵(质量矩阵，阻尼矩阵，刚度矩阵)至关重要。由于有限的基函数的自由度是分配给每个单元的边或点上的，所以一个单元中的物理场仅有此单元节点和边的物理场值与基函数决定，与其他单元的基函数无关，一般称这种方法中的基函数为局部基函数(local basis function)。所以构建有限元方程中的矩阵时，一般是根据计算矩阵元素的公式先计算每个单元的单元矩阵(element mattix)，最后再将各个单元的矩阵组装起来，得到最终需要的全局矩阵(global matrix)。
[0086]
s3，通过每个单元、节点、边的编号，将单元矩阵中的元素与全局矩阵的元素对应，通过遍历叠加利用单元矩阵，得到全局矩阵；
[0087]
由于有限元法使用的基函数为局部基函数(local basis function)，此基函数只在包含定义其自由度的点或边的单元中场值不为零，所以不相邻单元中的基函数一般是正交的，根据热学时域有限元的基本方程，利用这种有限元基函数得到的矩阵一般是稀疏矩阵。在matlab中实现有限元矩阵构造时，可以采用matlab中的稀疏矩阵(sparse matrix)的数据格式即相关的数学方法，这样可以减少程序的内存占用并加快程序运行速度。
[0088]
s4，通过热学时域有限元的基本方程对每一个时间步进行迭代运算；
[0089]
在进行矩阵线性运算时，可以通过lu分解等矩阵分解方法来提高程序的运算性能(采用矩阵分解法求解大型线性方程组的效率远高于直接对大矩阵求逆)。在具体毕业设计中，采用lupq矩阵分解法。利用matlab中的lu函数，将稀疏矩阵s分解为一个单位下三角矩阵l、一个上三角矩阵u、一个行置换矩阵p以及一个列置换矩阵q，并满足：
[0090]
psq＝lu
[0091]
热学时域有限元的基本方程都可以转换线性方程组求解的一般形式：
[0092]as
x＝b
[0093]
利用lupq矩阵分解方程：
[0094]
pasq＝lu
[0095]
可以得到用lupq表示的线性方程组：
[0096]
p-1
luq-1
x＝b
[0097]
从而得到x的解：
[0098]
x＝pl-1
u-1
qb
[0099]
经过大量的实验，发现在利用matlab进行线性方程组进行求解时，利用lupq方法的求解效率与不做矩阵分解的直接求解的效率相比有明显的优势，且在绝大多数情况下利用lupq求解方法进行线性问题求解的效率也比其他矩阵分解的效率高。
[0100]
第二实施例；
[0101]
本发明第二实施例是基于上述第一实施例形成，相同部分不再赘述，在上述第一实施例基础上增加下述步骤；
[0102]
s5，通过验证工具对热学时域有限元的基本方程结果进行验证。
[0103]
示例性的，验证和最终验证采用的参考软件为comsol multiphysics。comsol multiphysics最初也是于matlab平台作为应用开发的，在工业界和学术界广泛用于对偏微分方程描述的物理问题通过某些数学方法进行数值求解。comsol multiphysics支持使用的算法繁多，且有许多优化性算法以用于不同场景，但是这些算法本质上都属于有限元方法。作为一款跨平台的有限元分析，求解器和多物理场仿真软件。comsol multiphysics允许通过传统的基于物理学的用户界面进行建模并构造偏微分方程的耦合系统，是一款使用自由度较高的软件。comsol为机械，流体，电气，射频，声学，电化学反应、等离子体和光学等不同领域的仿真研究工作提供ide和统一的工作流程。
[0104]
为验证热学时域有限元求解器的准确性，本发明构建了几个算例，并与comsol multiphysics的结果进行比对。
[0105]
第一个算例为铜块冷却：0.1m
·
0.1m
·
0.1m的铜块，其热学参数：导热系数为377w/(m
·
k)，密度为8900kg/m3,恒压热容为386j/(kg
·
k),固体的边界条件为对流边界条件,对流流体的温度为300k，换热系数为500。
[0106]
为观察其温度变化，在铜块中心设置一域点探针，获取改点的温度。在comsol multiphysics和自行编写的fetd热学求解器中计算该程序，对比结果。其结果参考图2和图3所示。
[0107]
根据计算结果可以看出，自行编写的fetd热学求解器得到的结果与comsol multiphysics中的结果基本吻合，误差在可以接受的范围内。
[0108]
第二个算例为搭载散热片的体热源，参考图4所示。有四个2cm
·
2cm
·
2cm的体热源，与9片厚度为0.5cm的薄散热片，通过厚度为1cm的底座连接在一起。整个物体处于对流边界条件中，由于此模型设计的相关参数较多，因此不再赘述其具体参数，将于附录中给出具体的参数。
[0109]
在comsol multiphysics和自行编写的fetd程序中对此算例进行仿真。对比仿真结果得出，comsol multiphysics和自行编写的fetd程序的计算结果大体一致。通过两个算
例的对比验证，可以得出本发明的仿真结构是可靠的。
[0110]
第三实施例；
[0111]
本发明提供一种用于执行第一实施例或第二实施例任意一项所述的时域热传导仿真方法中步骤的计算机可读存储介质。
[0112]
除非另有定义，否则这里所使用的全部术语(包括技术术语和科学术语)都具有与本发明所属领域的普通技术人员通常理解的意思相同的意思。还将理解的是，除非这里明确定义，否则诸如在通用字典中定义的术语这类术语应当被解释为具有与它们在相关领域语境中的意思相一致的意思，而不以理想的或过于正式的含义加以解释。
[0113]
以上通过具体实施方式和实施例对本发明进行了详细的说明，但这些并非构成对本发明的限制。在不脱离本发明原理的情况下，本领域的技术人员还可做出许多变形和改进，这些也应视为本发明的保护范围。