二维热传导有限元

二维热传导有限元分析学院：理学院姓名：邓高明学号：20091078日期：2010.6.27一、引言对于一般热传导问题，我们关心的是温度在介质内如何随位置而变化的，不管这种变化是由介质的边界条件引起的，还是由介质内部的热源产生的。我们希望确定系统中包括边界在内的不同点上的热流量。二、热力学基本概念：热量通过三种方式传递：热传导，热对流，热辐射。热传导是指由于介质内部存在温度梯度而引起的热交换，且能量通过分子运动从高温区传递到低温区。在二维直角坐标系下，热传递率遵循傅里叶定律：其中分别为X和Y方向的热传递率，k为介质的热传导系数，A是介质的横截面面积，为温度梯度。傅里叶定律也可表示为单位面积上的热传递率：其中，分别称为X，Y方向的热流量。当运动的流体接触到温度不同的表面时就会出现热传递。流体和表面间的总热传递率有牛顿冷却描述：，其中，h是热传递系数（可由试验得到），Ts是表面的温度，Tf是流动流体的温度。所有的物质都会发出热辐射。只要所讨论的物体的温度在有限的范围，这个规律就是正确的，表面发出的能量由以下方程确定：，其中，q//表示表面发出的单位面积上的热流率，ε为表面的辐射率，且0＜ε﹤，σ为斯特潘-波尔兹曼常数，约为5.67x10-8W/m2·K4。热辐射不同于热传导和热对流，它可以在真空中发生，并且所有的物质都会发出热辐射。稳态的二维热传导问题，在直角坐标系下，对系统应用能量守恒定律，得到热扩散方程：q表示单位体积产量的热量。热传导问题中常用的几种边界条件：1、忽略通过表面损失或得到的能量。一般是指绝热表面或绝缘性相当好的表面。在热传导问题中，对称线也表示绝热线。这类边界条件由下式表示：22220XYTTkkqXY(0,)0XYTX2、表面的热流量为常数，这类边界条件，可用以下方程表示：3、由对流引起的冷却或加热过程仅发生在物体表面，这类边界条件，可用以下方程表示：//(0)oXTkqX(0)[(0,)]fXTkhTyTX三、热力学有限单元公式推导：1、矩形单元：()[]ijeijmnmnTTTSSSSTT其中，形函数分别是：ijmnSSSS(1)(1)(1)(1)ijmnxyxySSlwlwxyyxSSlwwl对热扩散方程应用伽辽金方法，在局部坐标系x,y下，得到余差方程：TTiXYyxTnTmTjijmn矩形单元在有限元分析我们曾提到过，在所以单元组合完成后，应将余差设为零，将上述方程改写为：22()2222()2222()2222()22()()()()eiixyAejjxyAemmxyAennxyATTRSkkqdAxyTTRSkkqdAxyTTRSkkqdAxyTTRSkkqdAxy其中，ST为形函数的转置矩阵:2222()0TxyATTSkkqdAxyijTmnSSSSS对上式进行积分计算，可以得到双线性矩阵单元的传导矩阵：在矩形单元不同边上应用对流边界条件：()2211211222111221112212216611222112yexklkwKlw()()()()21000000120002100000012066000000000000200100000000002100006600121002ijjmeeeemnnihlhlKKhlhlKK单元的热荷载矩阵：()()()()10110122000100102211fijfjmeefmnfnieehTlhTlFFhTlhTlFF2、三角形单元：()[]ieijkjkTTSSSTT其中，形函数分别是：ijkSSS1()21()21()2iiiijjjjkkkkSXYASXYASXYAXTYTkTjTi(Xk,Yk)(Xi,Yi)(Xj,Yj)三角形单元A是单元面积，可以从以下方程计算的出：且:2()()()ijkjkikijAXYYXYYXYYijkkjijkikjjkiikjkijikkijjikijkjiXYXYYYXXXYXYYYXXXYXYYYXX应用伽辽金法，则三角形单元的三个余差方程的矩阵形式为：2222()0TxyATTSkkqdAXY其中，ST为形函数的转置矩阵:iTjkSSSS对上式积分得到三角形单元的传导矩阵为：22()222244iijikiijikeXYijjjkijjjkikjkkikjkkkkKAA三角形不同边上应用对流边界条件：()()()210000120021660000122010216102ijjkeeekihlhlKKhlK单元的热荷载矩阵：()()()101110222011fljfjkfkieeehTlhTlhTlFFF四、模型计算：混凝土结构的小工业烟筒，外正方形边长60cm，内正方形边长20cm，热传导率κ=1.4W/m·K，如下图所示，假设烟筒内表面的统一温度为100。C。外表面暴露在空气中，空气的温度为30。C，相应的自然热对流系数为h=20W/m2·K.。确定稳态情况下温度在混凝土内的分布。分析：利用问题的对称性，只分析烟筒界面的1/8的面，将所选择的界面部分分成9个节点，5个单元。单元（1），（2），（3）是矩形单元，而（4），（5）是三角形。T=100°CTsurr=30°Ch=20W/m2·K(1)(2)(3)(4)(5)123456789YX1、矩形单元由热传导引起的传导矩阵为：()2211211222111221112212216611222112ekwklKlw单元（1），（2），（3）的大小是相同的，因此：(1)(2)(3)2211211222111221(1.4)(0.1)(1.4)(0.1)112212216(0.1)6(0.1)11222112KKK为便于单元组合，在每个矩阵的上侧和右侧列出每个单元相应的节点编号：2（j）3（m）4（n）1（i）4321(1)0.9330.2330.4660.2330.2330.9330.2330.4660.4660.2330.9330.2330.2330.4660.2330.933K4（j）7（m）6（n）3（i）6743(2)0.9330.2330.4660.2330.2330.9330.2330.4660.4660.2330.9330.2330.2330.4660.2330.933K(3)0.9330.2330.4660.2330.2330.9330.2330.4660.4660.2330.9330.2330.2330.4660.2330.933K5（j）8（m）7（n）4（i）5（j）(4)0.700.700.70.70.70.71.4K4（k）2（i）三角形单元（4），（5）都有如下传导矩阵如前所述，对流边界条件对传导矩阵和荷载矩阵都有贡献。对流边界条件对单元（2），（3）传导矩阵的贡献为：(5)0.700.700.70.70.70.71.4K8（k）9（j）5（i）6(2)00000000000.6660.333000.3330.666K473沿单元（5）的j-k边也会有通过对流损失热量，因此：(3)00000000000.6660.333000.3330.666K5874(5)00000.6660.33300.3330.666K985对流边界条件沿m-n边对单元（2），（3）热荷载矩阵的贡献为：(2)003030F(3)003030F对流边界条件沿j-k边对单元（5）的热荷载矩阵的贡献为：(5)03030F将所有的单元矩阵组合起来，整体刚度矩阵为：()0.9330.2330.2330.466000000.2331.6330.4660.933000000.2330.4661.8660.46600.2330.466000.4660.9330.4664.1990.9330.4660.4660.46600000.9332.33300.4660.9330000.2330.46601.5990.100GK000.4660.4660.4660.13.1980.100000.4660.93300.13.6650.36700000000.3671.366123456789123456789在节点1和节点2应用常温边界条件，得到的总刚度矩阵：()1000000000100000000.2330.4661.8660.46600.2330.466000.4660.9330.4664.1990.9330.4660.4660.46600000.9332.33300.4660.9330000.2330.46601.5990.100000.4660.4660.4660.13.1980.100000.GK4660.93300.13.6650.36700000000.3671.366在节点1和节点2应用常温边界条件，将会得到如下最终热荷载矩阵：()10010000030606030GF最后的节点方程组为：1000000000100000000.2330.4661.8660.46600.2330.466000.4660.9330.4664.1990.9330.4660.4660.46600000.9332.33300.4660.9330000.2330.46601.5990.100000.4660.4660.4660.13.1980.100000.4660123456789.93300.13.6650.36700000000.3671.36610010000030606030TTTTTTTTT求解线性方程组，则会