研究生《固体力学中的数值方法》课程论文

学号研究生课程论文课程名称《固体力学中的数值方法》题目温度梯度板单元的热模态理论分析学院专业班级姓名指导教师2015年6月日温度梯度板单元的热模态理论分析摘要：热变形和热模态分析是结构分析的重要内容。本文首先介绍了温度和热应力场的有限元，通过板单元面内温度梯度等效成热载荷，热变形问题转换成弹性问题；同时引入小变形条件下板单元的几何刚度矩阵，热刚度等于线性刚度和几何刚度矩阵的叠加；假设板单元厚度方向的温度呈线性分布，并拆分成对称和反对称两部分，分别采用平面单元和弯曲单元将它们等效成面内热载荷以及弯曲热载荷，将计算的结果叠加，得到温度梯度下板单元的变形，建立了温度梯度板单元的热模态计算理论。关键字：温度梯度，几何非线性，几何刚度矩阵，热模态Abstract：Thermaldeformationandthermalmodalanalysisareimportantpartsofstructuralanalysis.Firstly,introductionofFEMoftemperatureandthermalstressfield.Thermaldeformationproblemscanbetransformedintoelasticonesoncethrough-thicknesstemperatureisequivalenttothermalload.Then,introductiontogeometricstiffnessmatrixofbeamandplateundersmalldeformation.Thermalmatrixequalslinearmatrixandgeometricstiffnessmatrix.Assumethelineardistributionoftemperaturethroughthickness,andsplititintosymmetricalandanti-symmetricalparts.Andthenuseplanestressandblendingelementsrespectivelytocalculaterelativein-planeandblendingthermalloads.Byadditionofthem,derivesthedeformation.Keywords:Thermalgradient,geometricnonlinearity,geometricstiffnessmatrix，thermalmodes.1绪论1.1热模态研究的背景和意义模态分析是结构动力学的重要内容，热模态分析是模态分析中的重点。随着现代高速、超高速航空航天器的飞速发展，热模态分析的重要性更是得到加强。为了减轻结构的质量，增加内部空间，提高结构的效能，现代的航空航天器多采用薄壁结构。当飞行器高马赫飞行时，外表面由于暴露在大气下，不断受到大气摩擦，温度急剧升高，在某些条件下飞行器局部外表面的温度能够上升至600到2000℃。由于壁板不同位置受到的空气摩擦的具体情况不同，所以壁板面内会产生较大的温度梯度，从而产生较大的温度应力，使得壁板局部发生收缩和伸长。同时由于壁板外表面温度骤然升高，内壁温升明显慢于外壁温升，在这样的条件下，壁板的厚度方向也会产生较大的温度差，壁板的上下表面热膨胀有别，在此作用下，壁板会发生弯曲，并承受弯曲应力。由于气动加热的影响，不仅结构的材料参数发生变化，同时结构内部产生热应力，从而影响结构的整体刚度，在这样的高温条件下的模态，称之为热模态。由于高超音速飞机和导弹反应速度，突防能力和破坏力上的极大优势，纵观当代世界军事强国，无论是美国，英国，法国还是俄罗斯对高超音速飞机和导弹均投入极大的热情，中国也正加快研制步伐，在高超音速飞行器上加大研制进度。高超音速飞行器已经成为未来飞行器研制的重点。由于高超音速飞行器的飞行速度能够达到5马赫，如此高的马赫数必然引起极大的气动热效应。飞行器在气动载荷和热载荷的作用下，有可能导致飞行器的性能下降甚至破坏。作为制约超音速飞行器发展和研制的关键因素之一,国内外已经对气动热弹性问题进行了广泛而深入的研究，并且取得了大量的成果[1]。1.2本文研究内容随着现代计算机技术的飞速发展，结构的有限元计算精度与计算的资源消耗之间的矛盾已经得到极大程度的缓解，在当代的计算水平条件下，对合理的精度的要求比对资源消耗量的平衡显得更重要。本文主要研究了温度梯度板单元的热变形以及热模态问题，其主要分4析步骤可以表述为“热边界-热分布-热变形-热应力-几何刚度矩阵-热模态”，涉及到热传导，热应力，有限元，几何非线性等多个方面的理论知识。本文首先阐述了热传导的相关理论，建立了热传导的基本平衡方程和三大类热边界条件，在建立平衡方程和边界条件的基础上，由虚功原理推导出热传导的有限元计算公式，为热应力的计算做准备。热应力的分析通常基于杜阿梅尔-诺依曼提出的线性热应力理论，同时对于复杂结构的热应力分析计算通常基于有限元法，随后通过基于最小势能原理，通过势能的变分推导出了结构受热情况下等效热载荷，对于板的面内温度梯度，通过计算可以将其等效成面内拉压热载荷，而对于板的厚度方向的温度梯度，同样通过有限元计算，最后可以将其等效成弯曲热载荷。如果同时存在面内和厚度方向温度梯度的情况下，可以将温度分布拆分成上下均匀分布和上下反对称温度分布，然后将计算结果线性叠加即可。接着根据第二节的内容，在温度梯度或者热膨胀，热收缩受到约束的情况下，在结构的内部会产生热应力，在不考虑热弹耦合的情况下，热应力和结构受载荷作用等其他情况下产生的初应力并没有本质区别。初始应力或者应变的存在影响了弹性体的应力应变关系以及应变位移关系，为此一般通过格林应变公式增加非线性应变项来解决。格林应变比微小变形下的应变更精确，但是由于引入了非线性项，增加了计算的难度。为了提高计算效率，传统的非线性板单元对板单元的格林应变公式进行了有针对性简化，建立了在此基础上的非线性应变公式。本文中详细推导了初位移刚度矩阵以及初应力刚度矩阵的表达式，由于壁板在热载荷下的初始位移相对影响较小，所以本文中并不考虑初始位移的影响。然后在此基础上，从变分原理出发，较为详细地推导了传统的非线性板单元的初始位移刚度矩阵和初始应力刚度矩阵，从初应力刚度矩阵的表达式可以发现，只有中面内的初始应力对弯曲刚度产生影响，并不能对拉弯耦合，或者面内的刚度产生影响。最后则进一步，以格林应变等式为基础，从哈密顿原理出发，推导了体单元在初始应力下的初应力刚度矩阵，由于精度的提高，体单元的初应力刚度矩阵计算较复杂，从计算的结果可以看到，体单元的初应力刚度矩阵较为完备，初始应力之间相互产生着钢化或者软化作用。最后小节，从弹性力学的角度，推导了温度梯度下板单元的热弹性基本方程。从热弹性基本方程可以看出，对于面内的温度梯度可以等效成为面内的拉压热载荷来处理，而对于厚度方向的温度梯度，则可以等效成弯曲热载荷进行处理。为了计算温度梯度对几何刚度矩阵的影响，在第二小节中对普通板的几何刚度矩阵进行了修正，采用体单元和完整的格林应变公式计算出修正板的几何刚度矩阵，从而建立了温度梯度板单元的热模态理论。2温度场和热应力场的有限元对于简单的工程结构如板，圆筒，可以通过理论解析方法计算它们的温度场，热应力场，甚至是他们的热模态，但是对于那些形状不规则，边界条件复杂以及材料参数随着温度改变的热传导或者热应力问题，理论解答是不可能的，因此一般采用数值法解法求解。数值法以离散数学为基础，虽然不如解析法严密，精度没有解析法高，但是具有非常好的边界适用性。温度场的数值算法有多种，如差分法，有限元法和体积法等方法，但是有限元法在描述复杂边界时更具有明显优越性。本文通过有限元法计算结构的温度场和热应力场，然后以此为基础计算结构的热模态。2.1温度场有限元法2.1.1三维热传导微分方程的矩阵形式由热力学第一定律可知，物体在单位时间内获得的能量等于物体在单位时间内内能的增加和对外界做功的总和。QEA（2.1）对于固体，发生热膨胀的时候变性很小，所以由于体积膨胀或收缩对外界所做的功可以忽略。即物体单位时间获得的能量全部用于内能的增加。于是(2.1)可以表示为QE（2.2）在热传导中可以表示为123QQQ（2.3）其中，Q1和Q2分别表示单位时间物体通过表面和内部热源获得的热量，Q3则表示单位时间物体温度升高所需要的能量[2]。设有物体V，表面为S，在物体内取任一子域V1，表面积为S1，对于这个子域有如下等式11111ssQqndskgradTnds（2.4）同时物体内部体热源强度为W，则单位时间内子域内部产生的热量为211vQWdV（2.5）且子域温度升高所需要的能量可以表示为311vTQcpdVt（2.6）由式2.3可得11[()]0vTdivkgradTWcpdVt（2.7）在直角坐标中(2.7)的表达形式分别为222222()TTTTWatxyzcp（2.8）设微元体ΔV=ΔxΔyΔz，内部温度分布为θ=θ(x,y,z,t)，热源强度为g=g(x,y,z,t)，则微元体内任一点P在Δt时间内的热流量平衡方程为(qq)yzt(qq)zt(qq)tgtpxxyyxxxxxxxxzzxxxxcxyzxxyxyz（2.9）式中xq，yq，zq分别是以x,y,z轴为发现的平面上的热流分量。p和pc为介质的密度和比热容。等式两边同时除以ΔxΔyΔzΔt，可得到P点的热平衡偏微分方程：()gyxzpqqqcxyz（2.10）对于正交各向异性材料，gpxxyyzzckkkxxyyzz（2.11）如果介质的主轴是沿着x'，y'和z'方向，而不是沿着x，y和z方向，则在两个坐标系中的温度梯度或者热流矢量的有如下的关系111213212223313233xxnnnnnnyynnnzz（2.12）111213212223313233xxyyzzqnnnqqnnnqqnnnq（2.13）式中1jn，2jn和3jn分别表示x'，y'和z'坐标轴的方向余弦。设沿着坐标系Ox'y'z'中坐标轴的正交各向异性热传导率分别为xk，yk和zk，由此得出如下热流关系式111213212223313233xxyyzzqnnnqqnnnqqnnnq（2.14）由（2.12）（2.13）和（2.14）可以得到111213111213212223212223313233313233000000xxyyzzxqnnnknnnqnnnknnnyqnnnknnnz（2.15）因此得到任意坐标系下人传导率矩阵的表达式TDNDN（2.16）其中111213212223313233nnnNnnnnnn,000000xyzkDkk。2.1.2稳态热传导问题的有限元解法设在子域V内是任一连续函数，V的表面由S1，S2和S3构成。在S1上有10s（2.17）对于正交各向异性材料的稳态热传导，根据式(2.11)有如下表达式g0xxyyzzkkkxxyyzz（2.18）用乘以(2.18)两端，在体元V上积分，可得到如下等式()()xxyyzzvxxxyyyzzzSvkkkdVxxyyzzknknkndSgdVxyz