1.3.1　复合棒材稳态热传导的有限元计算

已知一个直径为20mm复材圆棒，如图1-3-1所示。其中长度为100mm的棒材热传导系数为100 W·m ⁻¹ ·℃ ⁻¹ ，长度为200mm的棒材热传导系数为20 W·m ⁻¹ ·℃ ⁻¹ ，长度为400mm的棒材热传导系数为80 W·m ⁻¹ ·℃ ⁻¹ ，圆棒周向绝热，左侧输入热流4000 W·m ⁻² ，右侧保持温度为80℃，求稳态条件下棒材的温度分布。

注意

为保证3个零件模型共节点，在Spaceclaim目录树上单击SYS，在其下的Properties中展开Analysis，将Share Topology设置为Merge。

1．建立单元热传导矩阵

对于一维（仅 x 向）、稳态、截面积一致和无内热源的模型，其第一类边界条件为（给定壁面温度）：

其中， ，即热流出口温度与入口温度之差； ，即单元长度； A 为单元截面积； k 为热传导系数。

对于热量入口，上式变换为：

根据能量守恒定律，热量出口为：

列成单元矩阵为：

2．建立整体单元热传导矩阵

简化计算，将每种材料的棒材仅划分为一个单元，从左至右节点编号依次为1、2、3、4，则每个单元的热传导矩阵为：

整体热传导矩阵由各单元热传导矩阵叠加而成：

对于本例，模型内部没有热源和能量损失，则 ，上面的矩阵简化为：

计算可得 T ₁ =144℃、 T ₂ =140℃、 T ₃ =100℃、 q ₄ =−4000 W·m ⁻² 。

提示

矩阵计算用Excel也可以轻易完成，不需要另装其他软件。以本计算为例，矩阵形式为：

简化为：

如图1-3-2所示采用Excel完成矩阵计算，先计算逆矩阵，框选G91:I93格，在插入函数中输入=Minverse(B91:D93)，按Ctrl+Shift+Enter组合键即可得到逆矩阵；再框选M91:M93格，在插入函数中输入=MMULT(G91:I93,K91:K93)，按Ctrl+Shift+Enter组合键即可得到 T ₁ 、 T ₂ 、 T ₃ 的结果。同理计算MMULT(B91:E94,M91:M94)，即可得到 q ₄ 的结果。

注意

本计算是基于每种材料的棒材仅划分为一个单元，如果划分为两个单元，同理计算可得 T ₁ =144℃、 T ₂ =142℃（位于原1单元中间）、 T ₃ =140℃、 T ₄ =120℃（位于原2单元中间）、 T ₅ =100℃、 T ₆ =90℃（位于原3单元中间），此结果正确；但是 q ₂ =−4.5E−11、 q ₃ =−2.3E−11、 q ₄ =0、 q ₅ =0、 q ₆ =1.13E−11，这与模型内部没有热源和能量损失有关。

3．Steady-State Thermal模块计算

如图1-3-3所示，分别创建a、b、c三种材料，其Isotropic Thermal Conductivity分别定义为100、20、80。

单击Geometry→Geom\SYS\a、b、c，在下方的Details窗口中展开Material选项，将Assignment分别对应设置为a、b、c，即选择定义的a、b、c材料。

Mesh设置保持默认。

边界条件设置如图1-3-4所示。其中选择a模型的左端面定义加载Heat Flux，数值为4000 W·m ⁻² ；选择c模型的右端面定义加载Temperature，数值为80℃。

计算完成后，分别插入Temperature（Geometry设置为All Bodies）、Temperature2（Geometry设置为a、b模型交界面）、Temperature3（Geometry设置为b、c模型交界面）和Total Heat Flux四项，后处理的计算结果如图1-3-5所示。

由图可知，系统最高温度为144℃，位于模型最左侧；a、b模型交界面温度为140℃；b、c模型交界面温度为100℃，总热通量为4000.1 W·m ⁻² 。与整体热传导矩阵计算结果相比，仅在总热通量上有0.0025%的误差。

4．伽辽金法求解稳态热传导分析

对于物理问题一般都可以采用微分方程进行表达，但是并不是所有方程都可以通过泛函的极值方法进行求解。当采用加权余量法时，如果试函数与形函数一致，即为伽辽金法。伽辽金法有广泛的适用性，有限元软件一般都基于该方法进行计算。下面采用二阶单元对本模型进行计算。

由于二阶单元为3个节点，基于等参化建立插值函数（不是单元的形函数）为：

利用插值条件，可得：

求出 a i 、 b i 、 c i ，代入插值函数，则：

令 ，则 （取 ）、 （取 ）、 （取 ），插值函数为：

根据伽辽金法，热传导矩阵中的任意一项为：

i,m =1, 3

又因为 ，对于任一单元，其中单元长度 ， ，则：

，

热传导矩阵中任一项为：

i,m =1,2,3，即

每种材料的棒材仅划分为一个单元，则a材料的热传导矩阵为：

同理，b、c材料的热传导矩阵分别为：

整体热传导矩阵为（每个材料单元的最后一个节点与下一个材料单元的第一个节点共点）：

则：

其中已知： Q ₁ /A =4000， Q _2-6 =0， T ₇ =80。

计算可得： T ₁ =144℃、 T ₂ =142℃、 T ₃ =140℃、 T ₄ =120℃、 T ₅ =100℃、 T ₆ =90℃（均无误）、 q ₂ =0 W·m ⁻² 、 q ₃ =−1.7E−11 W·m ⁻² 、 q ₄ =−1.1E−11 W·m ^-2 、 q ₅ =4.2E−12 W·m ⁻² 、 q ₆ =0 W·m ⁻² 、 q ₇ =−4000 W·m ⁻² 。对比线性单元，每个材料棒材划分两个单元的结果： q ₂ =−4.5E−11、 q ₃ =−2.3E−11、 q ₄ =0、 q ₅ =0、 q ₆ =1.13E−11，不论量级和数值，二次单元精度均高于线性单元。Excel计算如图1-3-6所示。其中逆矩阵为对传导矩阵前6×6矩阵求逆，热通量第五项、第六项计算基于高斯消元法（O79=−G79*H81，O80=−G80*H81），温度为逆矩阵与热通量矩阵之积。