自研复合材料有限元软件进展：局部坐标系定义

关键词：瞬态，热传导，有限元求解器，三角形单元热传递有三种方式：热传导、热对流、热辐射。就热传导问题而言，无论是结构力学还是流体力学都会涉及，两边都没拿它当外人。前面的文章提到过，结构力学的有限元发展得非常成熟，大部分的刚度矩阵在文献里面都推导好了。而流体力学的很多单元类型的有限元方程，可能需要自行推导完成。在热传导问题中，我采用加权余量法进行处理，推导出了符合结构力学有限元文献中给出的刚度矩阵，殊途同归。实际上，传统的结构力学有限元三大控制方程：几何方程、物理方程、平衡方程。几何方程描述位移-应变关系，物理方程描述应力-应变关系，平衡方程描述内应力-外载荷关系。传热问题从控制方程角度，更偏向流体力学（能量方程）。但是热对于结构变形太重要了，因此结构有限元必须要把传热问题解决掉。从结构力学跨到流体力学，在有限元方法中，流体力学控制方程左边的矩阵都可以用刚度矩阵去看待它。控制方程的右边的列阵，都可以用载荷的角度去看待，对于第二类边界条件，则可以分成左侧矩阵的修正+右侧列阵的载荷组合。有些文献上，用所谓的“内部单元方程”、“边界单元方程”的描述，会增加我们的困惑，可以不必纠结在此。控制方程二维瞬态热传导控制方程如下：这个方程里面的常数有密度、比热容、导热系数。三种边界条件：（1） 已知边界温度值，属于第一类边界条件，它的处理就和结构有限元里面的位移以一样，可以用置大数法对方程左边的矩阵进行约束处理。（2） 已知边界热流密度，属于第二类边界条件，作为热源。可以类比到结构有限元里面的均布载荷。（2） 已知边界对流换热系数和接触环境温度，也属于第二类边界条件。这个边界条件在处理的时候，需要进行拆分，一部分放到左侧单元矩阵，一部分作为右侧的载荷。有限元思路这部分在结构有限元教材中介绍的比较多，流程：（1） 根据单元类型，确定插值函数。此时单元温度用权函数表达。（2） 采用伽辽金方法，权函数=插值函数，控制方程与权函数相乘，积分取0。（3） 在每个单元域内，方程转换为权函数的积分形式，最终形成单元矩阵。单元方程使用三角形线性单元对应的插值函数：有些教材中，会把面积项提取出来，写成以下这种形式，所以有的教材上刚度矩阵结果用a、b、c表达的时候，会存在差异，但是本质都是一样的。最终单元方程如下，其中M是热容矩阵，K是传导矩阵，F是热载荷。热容矩阵乘的是温度的导数。在瞬态问题的求解中，导数项可以写成前后时间变量差值与时间间隔的比值：代入后得到如下形式：求解思路在求解过程中，把Tn+1当作未知量，Tn作为已知量。这样在每个时间点，求解方法和结构有限元方法一致。初始时候，可以指定一个温度作为全域已知初始温度，然后在迭代过程中，Tn和Tn+1会逐渐接近，达到收敛状态。案例效果设计案例如下，同时包含对流换热边界条件和热流，时间总长10000s，每步时间间隔50s。自研求解器和商用软件结果对比如下，从结果可以看出，自研求解器结果与商用软件结果一致。自研求解器结果：最终温度分布商用软件结果：最终温度分布自研求解器结果：平均温度时间曲线商用软件结果：平均温度时间曲线来源：静界有限元