计算传热学-传热基本原理和有限元应用

...wd......wd......wd...传热学的开展概述18世纪30年代首先从英国开场的工业革命促进了生产力的空前开展。生产力的开展为自然科学的开展成长开辟了广阔的道路。传热学这一门学科就是在这种大背景下开展成长起来的。导热和对流两种基本热量传递方式早为人们所认识，第三种热量传递方式则是在1803年发现了红外线才确认的，它就是热辐射方式。在批判“热素说〞确认热是一种运动的过程中，科学史上的两个著名实验起着关键作用。其一是1798年伦福特(B．T．Rumford)钻炮筒大量发热的实验，其二是1799年戴维(H．Davy)两块冰块摩擦生热化为水的实验。确认热来源于物体本身内部的运动开辟了探求导热规律的途径。1804年毕渥根据实验提出了一个公式，认为每单位时间通过每单位面积的导热热量正比例于两侧外表温差，反比例于壁厚，比例系数是材料的物理性质。傅里叶于1822年发表了他的著名论著“热的解析理论〞，成功地完成了创立导热理论的任务。他提出的导热定律正确概括了导热实验的结果，现称为傅里叶定律，奠定了导热理论的根基。他从傅里叶定律和能量守恒定律推出的导热微分方程是导热问题正确的数学描写，成为求解大多数工程导热问题的出发点。他所提出的采用无穷级数表示理论解的方法开辟了数学求解的新途径。傅里叶被公认为导热理论的奠基人。在傅里叶之后，导热理论求解的领域不断扩大。同样，自1823年M.Navier提出流动方程以来，通过1845年G.G.Stokes的改良，完成了流体流动基本方程的创立任务。流体流动理论是更加复杂的对流换热理论的必要前提，1909和1915年W.Nusselt开辟了在无量纲数原则关系正确指导下，通过实验研究对流换热问题的一种基本方法。1904年，L.Prandtl提出的对流边界层理论使流动微分方程得到了简化，1921年E.Pohlhausen基于流动边界层理论引进了热边界层的概念，为对流传热微分方程的理论求解建设了根基。在辐射传热研究方面，19世纪J.Stefan根据实验确定了黑体辐射力正比于它的绝对温度的四次方的规律，1900年M.Planck提出的量子假说奠定了热辐射传热理论根基。上述传热理论为传热分析解析、数值以及实验研究奠定了理论根基。还要特别提到的是，由于计算机的迅速开展，用数值方法对传热问题的分析研究取得了重大进展，在20世纪70年代已经形成一个新兴分支—数值传热学。近年来，数值传热学得到了蓬勃的开展[2-4]。传热分析计算理论热量传递主要有三种传递形式，分别是热传导、热对流和热辐射。热传导是指两个相互接触良好的物体之间的能量交换或一个物体由于其自身温度梯度而引起的内部能量的传递。其遵循傅里叶定律[5]：，其中是热导率，是温度梯度，是热流密度。热对流是指在物体与其周围介质之间发生的热量交换。热对流分为自然对流和强制对流，用牛顿冷却方程描述为，其中为外表传热系数，为物体外表的温度，为物体周围流体的温度。一个物体或两个物体之间通过电磁波形式进展的能量传递交换称为热辐射，通常由斯忒藩-波尔兹曼定律计算。就物体温度与时间的变化关系而言，热量的传递过程可以区分为稳态过程〔又称定常过程〕与非稳态过程〔又称非定常过程〕两类。但凡物体的各点温度不随时间而变化的热量传递过程都称为稳态热传递过程，反之温度随时间变化的热量传递过程则称为非稳态传热过程。2.1基本方程在进展传热分析时，主要用到的定律方程有能量守恒定律、动量守恒方程和质量守恒方程。能量守恒定律也是热力学第一定律，它是自然界基本的一个定律。它指出能量是不能消灭，也不能创造的，只能从一种能量形式转化为另一种能量形式，或者由一种物质传递到另一种物质，并且在这种能量转化和能量传递过程中其总量保持不变。同时，对流传热的描述还会用到动量守恒方程和质量守恒方程，动量守恒方程是描述粘性流体流动过程的控制方程。在数值模拟计算中，这些方程采用的是时均形式的微分方程。能量守恒方程质量守恒方程动量守恒方程式中：为流体压力；为流体温度；为流体所吸收的热量；为速度矢量；为流体的动力粘度；为作用在流体上的质量力，在重力场中；为导热系数；为流体的比热容；为能量耗散函数：其中为流体的变形张量，代表流体抑制粘性所消耗的机械能，他将不可逆转化为热而耗散掉；在充分开展的湍流区域，反映湍流脉动量对流场影响的湍流动能方程和湍流应力方程可以通过标准方程得到，其形式为：式中：为湍流粘度，；为湍流动能；为湍流动能耗散率。2.2基本控制方程求解的数值方法在利用数学方法进展热传递分析时，首先假定研究对象内各点的密度、温度、速度等都是空间坐标的连续函数。基本控制方程数值方法求解的基本思想是：把原来在时间和空间坐标中连续的物理场比方速度场、温度场等，用有限个一系列的离散点也就是节点上的值的集合来替代，再利用一定合理的原则建设这些有关离散点的表达变量值之间关系的代数方程即为离散方程，利用数学方法来求解所建设起来的这些代数方程并求得所求解变量的近似值。图2-1表示了基本控制方程的典型求解流程。图2-1控制方程数值求解流程计算流动传热常用到的数值方法主要包括：有限分析法、有限差分法、有限元法以及有限容积法。在有限元分析软件ANSYS中的有限容积法是指将计算区域划分成很多不互相重叠的网格，并且围绕每个网格节点都有一个控制体，再将每一个控制方程都在控制体上进展积分求解，可以得到包含一组节点计算变量值的离散化方程，可以保证具有守恒性，而且离散方程系数的物理意义明确，是目前流动与传热问题的数值计算中应用得最广的一种方法。有限元概述有限元分析方法是对真实的物理系统进展近似的数学模拟，用有限个单元去逼近无限未知量的过程。有限元的概念第一次提出是1943年Courant为研究St.Venant的扭转问题采用了三角形分片上的连续函数和运用最小势能的原理。有限元方法开展相当缓慢，直到1956年，Turner，Clough，Martin和Topp等人第一次真正通过运用直接刚度法来确定由弹性理论的方程求出三角单元特性解决平面应力问题，并且将其写入论文进展发表。由于计算机的出现，使得复杂的平面弹性问题求解更加容易，形成了新的研究方法。“有限元法〞这个名称，是1960年Clough发表的一篇平面弹性问题的论文中真正第一次出现。至此，工程师们开场注意到有限元法的作用，并把它进展广泛地应用。随着1970年代以后，计算机技术的飞速开展，也带来了有限元法的迅速开展进步，大量相关的学术论文相继发表，并且出现了更多相关专著，进入了有限元的全盛开展时期。迄今为止，有限元法主要被应用于流体力学、固体力学、电磁学、声学、热导学等各个领域；可以求解杆、梁、板、壳、块体等各类单元构成的弹性〔线性和非线性〕、弹塑性或塑性的问题；能计算温度场、电磁场、流体场等场分布问题的稳态和瞬态问题；还能求解水流管路、电路、润滑、噪声以及固体、流体、温度相互作用的问题[6]。有限元分析是建设真实的物理系统，包括几何条件和载荷工况，然后利用数学近似的方法进展数值模拟。有限元方法的主要基本思路是：“化整为零，积零为整〞。它的求解步骤包括：①将一个整体构造看作是由假设干个单个的构造元件构成，并且通过有限个连接点连接。单个的构造元件为“有限元〞或“单元〞，连接点为“节点〞。②在各单元上进展力学分析，并由相关力平衡条件建设相应的节点位移关系式及相关的节点力或节点位移的系列方程式。③求解得到的方程组，获得问题的解。如果形函数满足一定要求，解的精度会随着单元数目的增加而不断提高并且收敛于问题的准确解。如果无限制地增加单元的数目将会却增加计算机计算所消耗的时间。因此，在实际工程应用中，只要所得的数据能够满足工程需要就足够了。3.1有限元法的基本原理对于某些因为自身构造形状复杂或者具有非线性问题的工程实际问题，由于其边界值获得较为困难，解析解求解也很困难。这类问题主要可以由下面两种方法来解决，其一是把问题的连续体进展离散化处理，然后利用构造矩阵分析的方法进展处理，最后通过数值法进展求解；二是可以把实际问题进展相关简化处理，也就是进展简化问题的方程和边界条件，使其能够进展计算，进而求得它在简化状态下的解，但是过多的简化又可能造成求得的解不正确甚至是错误的。有限元法的基本思想是人为地把连续体的求解域划分成假设干单元，单元与单元之间只是通过节点相互连接，用构成一个单元的集合来替代本身的连续体。通过选定适当的插值函数求解单元内部各点的求解量。通过相关平衡关系或者能量关系来建设节点量相互之间的方程式，然后再将各个单元“集合〞在一块而形成总体的代数方程组，进展边界条件的求解。它是一种近似求解一般连续性问题的数值方法。3.2有限元法分析的步骤利用有限元法分析求解问题的基本步骤通常为：①明确问题，定义求解域：根据实际问题近似确定求解域的几何区域和物理性质；②有限元网络划分：将求解域近似划分为有限个具有不同形状和大小而且彼此相互连接的单元组成的离散域；③确定状态变量和控制方程：即用包含问题状态变量边界条件的微分方程来表示一个具体的物理问题，并且将其转化为等价的泛函形式进展有限元分析。④单元的推导：推导有限单元的列式指对单元构造一个适合的近似解，包括选择建设单元试函数和合理的单元坐标系，并且用某种适宜方法定义单元各状态变量的离散关系形成单元矩阵；⑤总矩阵方程求解：将单元总体组装成离散域的联合方程组，并且要满足一定的连续条件；⑥求解联立方程组：有限元法联立方程组的求解可用随机法、选代法和直接法。有限单元法分析从使用有限元程序的角度来分，又可以分成三大步骤，如图3－1所示。图3－1利用有限元程序进展分析的三个基本步骤前处理是指对研究对象进展网格划分并且形成计算模型的过程，主要包括选择计算单元类型、确定节点和单元网格及约束载荷的位移等。求解是指在得到总体刚度方程并进展约束处理后，联立线性方程组的求解，并且最后得到节点位移的总过程。后处理是指对计算结果的处理和有关数据的输出过程，包括各种温度、应变、应力或位移的整理，形成温度场分布图、应力图、变形图等[7]。3.3热传导问题的有限元描述从上述基本理论可以建设起热传导问题[8-9]的有限元描述方法。各向异性体传热问题能量方程为：(3.1)其中，，为单位面积热流向量分量，为单位面积内的热源率，为密度，为比热。对同一各向异性介质而言，傅立叶定律〔Fourier〕为：(3.2)其中为材料的导热系数张量。材料的密度、比热和导热系数均可随温度而改变。将Fourier定律代入能量方程中即可得到抛物线型传导方程式。结合边界条件和初始条件即可求解。初始条件设定零时刻的温度分布为：(3.3)考虑常用的边界条件：设定外表温度、设定外表热流及对流换热。上述边界条件为：在上在上在上其中，为可随时间变化的设定外表温度；，，为外表外法向的方向余弦；为单位面积的热流率；为对流换热系数，，分别为气流和内外表温度。将求解区域分成个单元而每个单元有个结点，单元内的温度及温度梯度可表示为：(3.4)以矩阵符号表示为：(3.5)其中为温度插值矩阵，为温度梯度插值矩阵。(3.6)Fourier定律写成矩阵形式，如下：(3.7)其中为导热系数矩阵。以结点温度表示温度梯度:(3.8)最后可推导出单元方程的矩阵形式:(3.9)其中(3.10)(3.11)(3.12)上面的式子中，为单元热容矩阵；，为与热传导以及对流有关的传导系数矩阵；，，分别为设定结点温度、内热源、外表热流和外表对流所产生的温度载荷向量。方程〔3.10〕为考虑传导和对流的各向异性介质的单元矩阵方程式，通过单元方程式组合可以很容易得到整体矩阵方程，引入边界条件和初始条件即可求解。如果系统的净热流率为0，即流入系统的热量加上系统自身产生的热量等于流出系统的热量，则系统处于热稳态。在稳态热分析中任一结点的温度不随时间变化，所以方程〔3.9〕中与时间有关的项都将消失，故稳态热分析的热传导方程简化为(3.13)式中：为传导矩阵，包含导热系数，对流系数及辐射率和形状系数；为结点温度向量；为结点热流率向量，包含热生成。本节内容为下一章ABAQUS有限元热传递分析计算奠定根基。ABAQUS软件[10]简介ABAQUS被广泛地认为是功能最强的有限元软件，可以分析复杂的固体力学和构造力学系统，特别是能够驾驭非常庞大复杂的问题和模拟高度非线性问题。在传热学领域，根据传热问题的分类和边界条件的不同，可以将热分析分成几种类型：与时间无关的稳态热分析和与时间有关的瞬态热分析；材料参数和边界条件不随温度变化的线性传热，材料和边界条件对温度敏感的非线性传热；包含温度影响的多场耦合问题。ABAQUS可以求解以下类型的传热问题：1〕非耦合传热分析。此类分析中，模型温度场不受应力应变场的影响。在ABAQUS/Standard中可以分析热传导、强制对流、边界辐射等传热问题，其分析类型可以是瞬态或稳态、线性或非线性。2〕顺序耦合热应力分析。此类分析中的应力应变场取决于温度场，但温度场不受应力应变影响。此类问题使用ABAQUS/Standard来求解，具体方法是首先分析传热问题，然后将所得到的温度场作为条件，进展热应力分析，得到应力应变场。分析传热问题所使用的网格和热应力分析的网格可以不一样。3〕完全耦合热应力分析。此类分析中的应力应变场和温度场之间有点强烈的相互作用，需要同时求解。可以使用ABAQUS/Standard或ABAQUS/Explicit来求解此类问题。4〕绝热分析；5〕热电耦合分析；6〕空腔辐射。本文通过使用ABAQUS建设铝合金厚板的热力耦合分析模型，其中铝合金的材料力学和热学材料属性都随着温度的变化而发生变化，分别使用顺序耦合分析和完全耦合分析，最后得到两组铝合金厚板的应力、应变、温度等的分布规律，对两组结果进展对比分析。计算实例分析如图1所示，一个两端固定的等截面铝合金厚板，铝合金厚板的弹性模量为E，线膨胀系数为α，在梁的上外表加随时间变化规律为Q(t)=20000+t(W/m2)的热流量和2×106(Pa)的均布载荷，求铝合金厚板上的应力、温度分布情况。铝合金厚板尺寸为：长0.8m，宽0.2m，高0.1m，密度为2810kg/m3，泊松比为0.33，选取材料为7075铝合金，7075铝合金材料属性如表1所示。表17075铝合金材料特性温度/导热系数/比热容/弹性模量/热膨胀系数/0155.083073.3322.625156.086072.0023.560158.387071.2024.0100161.090069.0724.9200175.097061.8728.4300185.0102053.8729.9400193.0112048.5331.4500197.0132044.5031.75.1顺序耦合热应力计算一、传热分析创立部件：翻开ABAQUS，进入Part模块，点击工具区“创立部件〞工具，输入部件名plate，选择三维、可变形、拉伸实体，模型空间大约尺寸设置为2，进入二维草图绘制界面。点击“创立线矩形〞，输入第一点坐标为(0,0)，第二点坐标为(0.8,0.1)，完成草图绘制，输入深度0.2，完成部件的创立，如图5-1所示。图5-1部件图图5-2材料属性定义图定义材料属性进入属性模块，点击“创立材料〞工具，输入材料名aluminumalloy，输入基本条件中的材料属性中的密度、导热系数、比热容、弹性模量和热膨胀系数，其中导热系数、比热容、弹性模量、热膨胀系数随温度变化发生变化，在输入数据时选择使用与温度相关数据，如图5-2所示。定义截面属性点击工具区中的创“建截面〞工具，输入截面名Plate-Section，类别选实体，类型选择连续，材料选取之前定义的aluminumalloy，完成截面创立操作，点击工具区中的分配截面，选择梁，将Plate-Section赋值给厚板。装配部件进入装配模块，点击工具区中的“将部件实例化〞工具，选择Plate，选择独立类型，完成装配定义。设置分析步进入分析步模块，点击工具区的“创立分析步〞工具，命名为Heat-flux,选择“热传递〞作为分析类型，基本信息为“瞬态〞响应，时间长度为600S，增量步选择“固定〞类型，最大增量步数100，增量步大小为10S，其他选项使用默认选项。可得增量步数应为60.定义载荷和边界条件进入载荷模块，点击“创立边界条件〞，分析步选择“Initial〞，类别选择力学中的“对称/反对称/完全固定〞，选择两个侧面，选中全部固定，加完力学边界条件。点击“创立预定义场〞，分析步选择“Initial〞，选择其他中的温度选项，选择整个厚板，输入温度大小为20℃，完初始温度设定。点击菜单栏中的工具选项，工具选项中的幅值，采用默认名称Amp-1，类型选择表，输入时间为0时刻值为20000，600时值为20600，完成对幅值的定义。点击“创立载荷〞，分析步选择“Heat-flux〞，选择热学中的外表热流，然后选中铝合金厚板上外表，大小为1，幅值选择刚刚定义的Amp-1，完成热流定义，结果如图3所示。图5-3载荷示意图图5-4网格示意图划分网格进入网格模块，可以看见构造显示为绿色，可知因构造形状很规则，可直接采用构造划网格。点击工具区的“为边布种〞，三条边布种数量为别为35，10，5，“网格控制属性〞选择六面体构造化网格，单元类型选择Standard，线性，热传递，得到八结点线性传热六面体单元DC3D8，然后点击为实例部件划分网格，网格划分完，点击检查网格，选择部件，进展检查，得到结果为：单元数:1850,分析错误:0(0%),分析警告:0(0%)，完成网格划分，结果如图4所示。运行计算进入作业模块，点击作业管理器，创立一个作业，各项均选择默认，创立完一个作业Job-1，提交作业，完成后点击结果，观察结果云图，保存模型FEM1.cae。二、力学分析翻开之前保存的模型FEM1.cae，仍使用之前所用模型，但是其中有几处需要修改：修改分析步(1)进入分析步模块，删除之前定义的分析步，重新建设一个静力学的分析步，名称为Load，类型选择“静力，通用〞，时间长度选择600S，增量中选择固定类型，最大增量步数输入100，增量步大小输入10，其他选择项默认选项。(2)点击场量输出管理器，删除自动生成的输出场量，新建一个场量，分析步选择选择默认的“Load〞，S，E，PE，U，NT为输出场量。修改载荷与预定义场(1)进入载荷模块，点击“创立载荷〞，分析步选择“Load〞，选择力学中的压强，然后选中铝合金厚板上外表，大小为2E6，完成压强定义。(2)点击“预定义场管理器〞，删除之前定义的预定义场，新建一个预定义场，分析步选择“Initial〞，选择其他中的温度选项，选择整个厚板，分布中选择“来自结果或输出数据的文件〞，选择之前计算完成的Job1.dob文件，分析步为1，增量为0，然后编辑“Load〞分析步中的预定义场，状态选择“已修改〞，开场分析步输入1，开场增量输入0，完毕分析步输入1，末尾增量输入60。完成对预定义场的修改。修改网格进入网格模块点击“指派单元类型〞，选择三维应力，其他选项使用默认选项，重新划分网格。运行计算进入作业模块，点击作业管理器，创立一个作业，各项均选择默认，创立完一个作业Job-2，提交作业，完成后点击结果，观察结果云图。保存模型FEM2.cae。5.2完全耦合热应力计算翻开之前保存的模型FEM1.cae，仍使用之前所用模型，但是其中有几处需要修改：修改分析步进入分析步模块，删除之前定义的分析步，重新建设一个静力学的分析步，名称为Temperature-displacement，类型选择“温度-位移耦合〞，时间长度选择600S，增量中选择固定类型，最大增量步数输入100，增量步大小输入10，其他选择项默认选项。修改载荷与预定义场(1)进入载荷模块，点击“创立载荷〞，分析步选“Temperature-displacement〞，选择力学中的压强，然后选中铝合金厚板上外表，大小为2E6，完成压强定义定义。再创立一个外表热流，然后选中铝合金厚板上外表，大小为1，幅值选择刚刚定义的Amp-1，完成热流定义。(2)点击“预定义场管理器〞，删除之前定义的预定义场，新建一个预定义场，分析步选择“Initial〞，选择其他中的温度选项，选择整个厚板，输入温度大小为20℃，完初始温度设定。修改网格进入网格模块点击“指派单元类型〞，选择温度-位移耦合选项，其他选项使用默认选项，重新划分网格。运行计算进入作业模块，点击作业管理器，创立一个作业，各项均选择默认，创立完一个作业Job-3，提交作业，完成后点击结果，观察结果云图。保存模型FEM3.cae。5.3计算结果分析通过ABAQUS计算，所得结果如以以下图，图5-5为顺序耦合热应力分析中第一步热传导分析所得温度云图，图5-6为在顺序耦合分析中，同时加上压力和温度时所得温度分布图，图5-7为完全耦合热应力分析中的温度云图，图5-8为顺序热传导中，同时加上压力和温度时的Mises应力云图，图5-9为完全耦合热应力分析中的Mises应力分布云图。通过对比图5-5和图5-6可知，在顺序热传导的分析过程中温度分布图一致，因此可以得知，在读取Job-29-1-1.dob过程中没有出错，可以先进展热传递分析，再进展应力分析。通过对比图5-6和图5-7可以得到在分析过程中顺序热力耦合分析和完全热力耦合分析所得到的温度分布图基本一致。通过对比图5-8和图5-9可以得到在分析过程中顺序热力耦合分析和完全热力耦合分析所得到的构造Mises应力分布图除了应力值大小外，应力集中所在区域基本一致，完全热力耦合分析中应