热传导过程有限元分析

第六讲热传导过程有限元分析元计算技术部传热学是研究温差引起的热能传递规律的科学。热力学第二定律指出：凡是有温差存在的地方，就有热能自发地从高温物体向低温物体传递。本讲针对热传导问题从其基本方程、有限元分析、ELAB工程建模等几个方面来介绍其仿真过程。基本方程ELAB模型向导实现有限元脚本文件分析传热过程的基本变量就是温度，它是物体中的几何位置以及时间的函数。根据Fourier传热定律和能量守恒定律，可以建立热传导问题的控制方程，即二维物体的瞬态温度场u(x,y,t)应满足以下方程：()()xyuuuckkQtxxyy其中u表示温度，kx，ky，kz是物体的沿x,y方向的导热系数，ρ是物体的密度（kg/m3），c是物体的比热(J/(kg·K))，Q是物体内热源密度(W/(m·K))传热问题的边界条件有三类：第一类边界条件：0uu第二类边界条件：0xxyyuunknkqxy第三类边界条件：0xxyyuunknkhuuxy其中，q0是边界上热流的给定值,nx,ny分别为边界表面外法线在方向的方向余弦，h表示物体与周围介质的热交换系数，u0表示环境温度。传热问题的基本方程有限元分析针对二维问题，根据上面的瞬态热传导方程可得其积分形式为：(()()xyVVuuuckkQudVtxxyy其中，δu为温度的虚位移()()xyVxxyyVuuuuucukkdVtxxyyuuQudVnknkudxy将边界条件代入上式(注意，对于已知温度边界条件，虚位移δu为0，可得:0()qxyVVuuuuucukkdVQudVqudtxxyy工程背景一个长方形截面的冷空气通道，几何模型如下图所示。假设在垂直于纸面的方向上，通道内的初始温度为0℃。通道的导热系数为0.044W/m·℃，比热和密度的乘积为1J/(m3·℃)，内壁维持在0℃，外壁与流体发生对流交换，且与周围环境间的热换系数为10W/m·℃，环境温度为30℃，求3s后通道壁面中的温度和热流密度。几何模型瞬态热传导有限元分析工程建模1、点击“工程向导”进入公式库2、选择“热学”研究领域3、选择“坐标系”瞬态热传导ELAB1.0软件实现5、选择“问题类型”4、选择“单元类型”6、定义工程名和工程路径，完成工程设置定义材料参数点击工具栏“参数设置”→“材料参数”，如下图所示：材料参数对话框中设定相应的材料参数，如下图所示：a场体单元材料参数图a场边界单元材料参数b场单元材料参数几何建模：点击工具栏中“前处理”按钮进入GID。首先建立一个小的矩形面，利用gid中copy命令中的拉伸功能建立如下图所示的几何模型，详细步骤可以参考《有限元分析基础与应用》相关章节。有限元计算模型前处理注：进入GID后要进行ELAB1.0的数据转化data→problemtype→ELAB•施加材料属性：在condition窗口中为a场(温度)和b场(热流)分别施加材料属性和边界条件，该模型只有一种材料，材料赋值如下图所示：模型内壁保持0℃，外壁与外界发生对流交换(由边界条件文件来实现，在gid中通过赋边界材料来实现)，边界赋值如下图所示：a场面材料添加b场面材料添加•施加边界条件：a场固定温度a场对流交换边界•划分网格：为模型各条边设置尺寸0.1，划分如下图所示的有限元计算网格：结构化四边形网格(网格尺寸0.1)设置时间步长为1s，总计算时间为3s。点击工具栏中“求解计算”按钮，完成模型的求解计算。点击工具栏中的“后处理”按钮进入GID，查看计算结果，以下云图是3s时的计算结果。温度场u分布云图热流场x方向分布云图热流场y方向分布云图工程求解后处理有限元语言描述文件为生成该问题有限元计算的所有程序源代码，针对之前的ELAB有限元分析得到的微分方程弱形式，ELAB软件提供简洁的有限元语言描述文件，包括微分方程描述文件、多物理场描述文件以及求解命令流控制文件。针对该问题的有限元描述文件包括heatxy.fde(温度场fde文件)，hfxy.fde(热流场fde文件)，heat.mdi，heat.gcn在heatxy.fde给出单元的待求未知量，涉及到的材料参数，单元的形函数表达式，刚度矩阵表达式和载荷表达式，以及为描述刚度矩阵和载荷向量而自定义的函数。以下给出微分方程描述文件中与微分方程弱形式对应的部分(详细的解析见《有限元分析基础和应用》中相关章节)：微分方程描述文件heatxy.fde(温度场fde文件)DISPu未知变量定义微分方程弱形式中的变量uMATEekecq1.01.00.0材料参数行对应微分方程弱形式中的变量(考虑各向同性材料，各向热传导系数相同即kx=ky=ek)kx(ky)ρcq微分方程弱形式：未知变量：材料参数：0()qxyVVuuuuucukkdVQudVqudtxxyydist=+[gu_i;gu_i]*ek*vol(其中gu是一向量，其分量为vectguguxguygu的表达式在该fde中对应：@lgrad.xyffe@wgufe也就是未知量对x和y的导数。)单元刚度矩阵对应微分方程弱形式中的左端第二项单元刚度矩阵：yuuuukxxyymass%1ec*vol单元质量项对应微分方程弱形式中的左端第一项，其中的ec表示密度ρ与比热容c的乘积单元质量矩阵：ucut多物理场描述文件heat.mdi2dxy#a11ufdeheatxyq2fbcheatxyl2#b02uxuyfdehfxyq2#坐标系(二维直角坐标系)a场1个初值1个自由度u方程描述文件+单元类型和积分方法边界描述文件+单元类型和积分方法b场0个初值2个自由度ux，uy方程描述文件+单元类型和积分方法结束标志load=+[u]*q*vol单元载荷向量对应微分方程弱形式中的右端项单元载荷向量：VQudV求解命令流控制文件heat.gcnDEFIaparbbstraSTARTal1:BFTSOLVCaSOLVSTRbagidres(coor0);if(stop==0)gotol1;a场+算法b场+算法+耦合场(空一行)初始化a场时间循环开始标志时间和边值更新直接法求解a场最小二乘求解b场输出gid格式的结果文件循环结束标志