热弹性与有限元数值仿真(传热学2)

1、pvTTTTCqtxxyyzz 求解导热问题的关键是获得温度场求解导热问题的关键是获得温度场, ,而要获得温度场实而要获得温度场实质上归结为对如下导热微分方程式的求解。质上归结为对如下导热微分方程式的求解。 对上述偏微分方程：对上述偏微分方程： 对实际工程问题用纯数学的方法来解微分方程对实际工程问题用纯数学的方法来解微分方程非常困难；非常困难； 利用计算机来获得满足工程要求的数值解利用计算机来获得满足工程要求的数值解计算机数值仿真。计算机数值仿真。 一、常用的数值计算方法：一、常用的数值计算方法： 1 1。有限差分法、。有限差分法、2 2。有限单元法、有限单元法、3 3。边界元法等。边界元法等

2、二、二、有限元分析有限元分析软件的应用软件的应用： 目前，有限元理论及其应用已经很成熟，有许多目前，有限元理论及其应用已经很成熟，有许多商业软件可应用，如：商业软件可应用，如：ANSYSANSYS、PHOENICSPHOENICS、KIVA-2KIVA-2等。等。 导热问题数值解法的求解原理导热问题数值解法的求解原理 理论解 在规定的边界条件下积分，有很大局限性； 数值解 借助计算机，前景广阔。1.有限差分法原理（连续的问题 离散的问题） 以有限差分 无限微分 无限划分 实质 达到精度 以差分代数方程 微分方程 计算机帮助 （当离散点足够多时可以满足要求） xtdxdt 下面先对稳态导热问题中

3、位于计算区域内部的节点(简称内节点)介绍其离散方程的建立方法，而位于边界上的节点及非稳态导热中的非稳态项的离散将在以后讨论。 为讨论方便，把如图中的节点(m，n)及其邻点取出并放大，如图所示。图4-3 内节点离散方程的建立 下面以一个二维导热问题为例进行分析：把一个二维物体在X及Y方向上分别以 及 距离分割成矩形网格。则其中节点（m，n)的坐标为：X=m , Y=n ,其余节点类推。（举例） 三种基本差分格式：以节点(m，n)为例（1）向前差分：（2）向后差分：（3）中心差分：XYXYyttytxttxtnmnmnmnmnmnm,1, 1,yttytxttxtnmnmnmnmnmnm1, 1,

4、yttytxttxtnmnmnmnmnmnm2/1,2/1, 2/1, 2/1,对无内热源、稳态、二阶导热微分方程，有： 用中心差分格式因为：所以：02222ytxtxttxtxttxtnmnmnmnmnmnm, 1, 2/1, 1, 2/1,yttytyttytnmnmnmnmnmnm1,2/1,1,2/1,2, 1, 1, 2/1, 2/1,22)(2xtttxxtxtxtnmnmnmnmnmnm21,1,2/1,2/1,22)(2ytttyytytytnmnmnmnmnmnm2222ytxt2, 1, 1)(2xtttnmnmnm21,1,)(2ytttnmnmnm 最终得： 如果取正方

5、形网格，即取 ，则上式为： tm+1,n+tm-1,n+tm,n+1+tm,n-1-4tm,n=0 上式说明：在导热系数为常量时，热量的转移可用温度差来表达； 在稳态下，流向任何节点的热量的总和必须为零。 对于每个节点写出上式，然后联立求解方程组，即可求解。 （如边界温度已知，可逐步递推求解）yx 位于平直边界上的节点 这时边界节点(m，n)代表半个元体。如图4-4中有阴影线的区域所示。设边界上有向该元体传递的热流密度qw，于是该元体的能量守恒定律可表为图4-4 平直边界上的节点x=y时有(2)外部角点 在如图4-5所示的二维墙角计算区域中，节点AE均为外部角点，其特点是每个节点仅代表四分之一

6、个以x、y为边长的元体。今以外部角点D为例，假设其边界上有向该元体传递的热流密度qw，则其热平衡式为图4-5 外部角点与内部角点x=y时有（4-5a）（4-5b）(3)内部角点 图4-5中的F点为内部角点，代表了四分之三个元体。在同样的假设条件下有 x=y时有（4-6b）（4-6a） 现在来讨论关于边界热流密度的三种情况。 (1)绝热边界 令式(4-4)(4-6)中的qw0即可。 (2) qw值不为零 以给定的qw值代入上述方程，但要注意上述三式中以传入计算区域的热量为正。 (3)对流边界 此时qwh(tf tm,n)，将此表达式代入式(4-4)(4-6)，并将此项中的tm,n与等号前的tm,

7、n合并。对于x=y的情形有： 平直边界 当计算区域中出现曲线边界或倾斜的边界时，常常用阶梯形的折线来模拟真实边界，然后再用上述方法建立起边界节点的离散方程。例如，如要用数值方法确定如图4-6a所示二维区域的形状因子，显然，根据对称性我们只要考虑四分之一的计算区域即可。图4-6a中的内圆边界可以来用图4-6b所示的阶梯形的折线边界来近似。只要网格取得足够密，这种近似处理方法仍能获得相当准确的结果。处理不规则边界的更好的方法要用到坐标变换，这里不做介绍。图4-6 不规则区域的处理 下面讨论代数方程的求解方法。 前已指出，代数方程组的求解方法分为直接解法及迭代法两大类。直接解法是指通过有限次运算获得

8、代数方程精确解的方法，像矩阵求逆、高斯消元法等均属于此种方法。这一方法的缺点是计算所需的计算机内存较大，当代数方程的数目较多时使用不便。另一类方法称迭代法。在迭代法中先对要计算的场作出假设(设定初场)，在迭代计算过程中不断予以改进，直到计算前的假定值与计算后的结果相差小于允许值为止，称为迭代计算已经收敛。本书中只介绍迭代法。 迭代法中应用较广的是高斯-赛德尔(Gauss-Seidel)迭代法，现以简单的三元方程组为例说明其实施步骤。 设有一个三元方程组，记为 其中ai,j(i1，3，j1，3)及bi(i1，3)是已知的系数(设均不为零)及常数。采用高斯赛德尔迭代法求解的步骤如下： (1)将式(

9、a)改写成关于t1、t2、t3的显式形式(迭代方程)，如 (2)假设一组解(即迭代初场)，记为t1(0)、t2(0)及t3(0)，由式(b)逐计算出改进值t1(1)、t2(1)及t3(1)。每次计算均用t的最新值代人。例如当由式(b)中的第三式计算t3(1)时代入的是t1(1)及t3(1)之值。 (3)以计算所得之值作为初场，重复上述计算，直到相邻两次迭代值之差小于允许值，此时称为已达到迭代收敛，迭代计算终止。 判断迭代是否收敛的淮则一般有以下三种： 其中上角标k及(k十1)表示迭代次数，tmax(k)为第k次迭代计算所得的计算区域中的最大值。当计算区域中有接近于零的t时，采用式(4-10c)

10、比较合适。允许的相对偏差之值常在10-310-6之间，视具体情况而定。43 非稳态导热问题的数值解法 非稳态导热与稳态导热的主要差别在于控制方程中多了一个非稳态项，而其中扩散项的离散方法与稳态导热是一样的。因此，本节讨论重点将放在非稳态项的离散以及扩散项离散时所取时间层的不同对计算带来的影响上。1.泰勒展开法 首先以一维非稳态导热为例讨论时间空间区域的离散化。如图4-8所示，x为空间坐标，我们将计算区域划分为(N-1)等份，得到N个空间节点；为时间坐标，我们将时间坐标上的计算区域划分为(I-1)等份，得到I个时间节点。从一个时间层到下一个时间层的间隔称为时间步长。空间网格线与时间网格线的交点，

11、如(n，i)，代表了时间空间区域中的一个节点的位置，相应的温度记为tn(i)。 非稳态项 的离散有三种不同的格式。如果将函数在节点(n，i+1)对点(n，i)作泰勒展开，可有 由式(b)可得在点(n，i)处一阶导数的一种差分表示式 , 的向前差分: 类似地，将t在点(n，i-1)对点(n，i)作泰勒展开，可得 的向后差分的表达式： 如果将t在点(n，i+1)及(n，i-1)处的展开式相加，则可得一阶导数的中心差分的表达式： 在非稳态导热问题的数值计算中，非稳态项的上述三种差分格式都有人采用，本书主要采用向前差分的格式，但也简单介绍了向后差分的格式。 如果把式(4-14a)中的扩散项也用(i+1

12、)时层上的值来表示，则有 式中已知的是i时层的值tn(i)，而未知量有3个，因此不能直接由上式立即算出tn(i+1)之值，而必须求解(i+1)时层的一个联立方程才能得出(i十1)时层各节点的温度，因而式(4-15)称为隐式差分格式。从时空坐标系中的节点(n，i+1)来看，式(4-15)的左端是非稳态项的一种向后差分。隐式格式的缺点是计算工作量大，但它对步长没有限制，不会出现解的振荡现象。 以上是将一维非稳态导热方程中的两个导数项用相应的差分表示式代替而建立差分方程的，这种方法称为泰勒展开法。2.热平衡法 这种方法不受网格是否均分及物性是否为常数等限制，是更为一般的方法。 图4-9示出了一无限大

13、平板的右面部，其右侧面受到周围流体的冷却，表面传热系数为h。此时边界节点N代表宽度为x/2的元体(图中有阴影线的部分)。对该元体应用能量守恒定律可得经整理得 式中 是以x特征长度的傅里叶数，称为网格傅里叶数， 一项可作如下变化： 式中Fo及Bi分别为网格傅里叶数及网格毕渥数。于是式(4-16b)又可改写为 至此，我们可以把第三类边界条件下、厚度为2的无限大平板的数值计算问题作一归纳。由于问题的对称性，只要求解一半厚度即可，其数学描写见式(3-10)(3-13)，此处不再重复。设将计算区域等分为N-1等份(N个节点，见图4-10)，节点1为绝热的对称面，节点N为对流边界，则与微分形式的数学描写相

14、对应的离散形式为 其中式(4-20)是绝热边界的一种离散方式，在确定t1(i+1)之值时需要用到t-1(i) 。根据对称性该值等于t2(i)。这样，从已知的初始分布t0出发，利用式(4-17)及(4-19)可以依次求得第二时层、第三时层直到I时层上的温度值(见图4-8)。至于空间步长x及时间步长的选取，原则上步长越小，计算结果越接近于精确解，但所需的计算机内存及计算时间则大大增加。此外，x及之间的关系还受到显式格式稳定性的影响。下面我们从离散方程的结构来分析，说明稳定性限制的物理意义。 式(4-17)的物理意义是很明确的。该式表明，点n上i+l时刻的温度是在该点i时刻温度的基础上计及了左右两邻

15、点温度的影响后得出的。假如两邻点的影响保持不变，合理的情况是：i时刻点n的温度越高，则其相继时刻的温度也较高；反之，i时刻点n的温度越低，则其相继时刻的温度也较低。在差分方程中要满足这种合理性是有条件的，即式(4-17)中tn(i)前的系数必须大于或等于零。如用判别式表示，则为必须保证 ANSYSANSYS软件在求解柴油机零部件温度场的应用软件在求解柴油机零部件温度场的应用 180活塞二维轴对称模型稳态温度场 180活塞三维轴对称模型稳态温度场 二维结构耦合系统循环瞬态温度场动画演示 三维结构耦合系统循环瞬态温度场动画演示第二部分第二部分 有限元法原理有限元法原理v有限差分法要对网格规则划分而

16、有限单元法可以进行不规则的划分v缺点：死板，不利于复杂问题求解v变分原理可以允许网格任意划分v变分法是有限元法的基础v如何把微分方程转变为有限元法的基本计算公式?v古典的偏微分方程的近似（无限逼近达到要求）计算：1.有限差分法 2.变分法有限单元法 泛函的基本定义v定义函数y=f(x) 则定义泛函J=Jy(x) vJ值依赖于自变函数y(x)的变化v函数y(x)又依赖于自变量x的变化v泛函是函数的函数v研究泛函的极值问题的方法变分法v研究函数的极值问题的方法微分法v采用变分法求泛函J=Jy(x)的极值问题（变分法）v最简单的一维泛函一般形式Jy(x)= v变分v变分运算与微分运算几乎相同（注意在

17、变分计算中自变量x做常数处理v因为微分方程求解通常都比较困难，而计算泛函变分求极值比较方便，所以常常把微分方程的求解问题转化为求相应泛函的极值问题21)(),(,1xxdxxyxyxF)(xyJJv举例：v求泛函Jy(x)= ，在边界条件y(0)=0，y(1)=1下达到极值的曲线。v解：这里 v则v 既v一维欧拉方程dxxyy1022)(xyyyyxF2)(),(211 2)(2,2yyFdxdyyFxyFxyyx 02212 21cxy21361cxcxy210)(xxydxyFdxdyFJ即极值条件是： 利用边界条件：y(0)=0，y(1)=1得到： 极值曲线0)(yFdxdyFxxcxc

18、xy6561613213泛函，变分举例分析 在铅垂平面xoy上有2点A，B（它们不在同一铅垂线上），如图所示：连接A，B两点的曲线有很多，要求找到一条曲线，使垂物W凭自重又A点沿此曲线滑到B点（忽略摩擦）所需的时间最短？ 解：两点间直线距离最短，但所需的时间不是最短 连接A，B两点的曲线很多，设每一条曲线的 函数为 而每一条曲线对应一个时间量T，那么T即函数y (x)的泛函，T=Ty (x)4 . 3 . 2 . 1)(ixfyivA(0.0) Xv v p(x,y) B(x,y) v VY要求最短时间，即求T=Ty (x)的极值（最小值）泛函的变分问题21)(),(,Jy(x)1xxdxxy

19、xyxF二维欧拉方程v二维泛函的一般形式 v其中v变分v其中dxdyuuuyxFDyx,y)(x,Ju yyxuuxyxuuyx),(,),(dxdyuuFuuFuuFJyyxdx)(),(uyyuuuxxuuyx最后可变为极值条件是即二维欧拉方程为其中udxdyFyFxuFJyxuud)()(0J0yFxFuFyxuuyuxuuFFuFFyx,瞬态温度场的变分原理公式边界条件为第二类其中 为换热系数,Tf为周围介质温度,求其温度函数 02222ytxt)(Ttntf),(yxtt 这个温度场的泛函为:应用二维欧拉方程可证明,为极值函数时,必须满足微分方程 和边界条件参考书:(机械结构分析的有

20、限元方法,杨莱柏主编) dstttdxdyytxtyxtJRf)21()()(2),(2202222ytxt)(TTntf函数和泛函极值的类比若函数y=f(x)在x有极值,则 ， 为极值点若泛函在曲线 上达到极值,则在此曲线上有: 为极值函数,怎么把变分和有限元法联系到一起?0)(0 xf0 x)(0 xy0)(0 xyJ)(0 xy里兹法设要求得满足微分方程以及边界条件 的函数在区间a,b连续 ,这个问题等价于求解相应的泛函 达到极值时的极值函数)(22xfdxdBbAa)(,)()(),(xfxbadxxxfdxdxJ)()()(21)(2)(x里兹求解法极值函数展开成 是满足边界条件的可

21、能解 是待定系数则上式可变为:这样,泛函可看成n个变量的多元函数,对于n元函数的极值条件: niiinnxxxxxx1332211)()()()()()()()(1xxnn1,)(21nJxJ0,0,0,213212211nnnJJJ有限元法的变分解释对二维问题:由边界条件其泛函为:把求解区域D划分成有限个子区域,分成n个三角形单元,单元记作De02222yuxu),(),(yxfyxucdxdyyuxuyxuJD)()(21),(22则整个区域D为各个单元的组合同样:求解区域D的泛函可写成个单元泛函的组合. 若为三角形单元,有三个节点,按照里兹法,将函数近似地展开成K为总节点,由泛函极限条件

22、: neeDD1neeeyxuJyxuJ1),(),(),(yxuemmjjiieuNuNuNyxu),(,),(21kuuuJyxuJ0J即得到线性方程组,求k个节点的值从而获得整个区域的解.k) 1,2=(i, 0iuJkuuu21,单元剖分和温度场离散温度插值函数 YO XYO X),(jjJyxT),(iiiyxTmTmSjS),(yxTiS 用有限元法将区域D剖分成 个单元,n个节点,把区域内连续的温度场离散到几个节点上去,最后只需要求节点上的温度. 对于三角形单元内的温度T,假设其为x和y坐标的线性函数(因为单元足够小,总可以做到),即: 是待定系数 同理3E321321,aaay

23、axaaTmmmjjjiiiyaxaaTyaxaaTyaxaaT321321321矩阵简化写成 即 mjimmjjiiTTTaaayxyxyx321, 1, 1, 1mjimjimjimjiTTTcccbbbaaaaaa,21321)()()(21mmmmjjjjiiiiTycxbaTycbaTycxbaTmmjjiiTNTNTNTkjimjiTTTNNNT,eeTNT导热微分方程线形方程组对于第二类边界条件二维稳态温度场的微分方程为: 相应泛函为:对泛函求极值,极值条件为:)(, 02222TTnTyTxTfdSTTTdxdyyTxTyxTJkfd)21()()(2),(222ekekTJT

24、T0其中 为单元泛函, 为求和(三角形单元k为i,j,m)我们对内部单元有而我们前面推导得出,对三角形单元有: 所以 同理: eJeeiiiedxdyyTTxTTTJ)()(2220)()(eiidxdyyTTyTxTTxTijkmmjjiiTNTNTNT)(21)()()(mmjjiimmjjiiTbTbTbTNxTNxTNxxT)(21TcTcTcyTmijii以及 ijk注意到所以有:同理:求得. iiicyTTbxTT21)(,21)(edxdy0)()()(422mmimijjijiiiiieTccbbTccbbTcbTJmejeTJTJ,得到 mjimmmjmjjjmimijiji

25、iimejeieTTTcbccbbcbccbbccbbcbTJTJTJ222222, 0 , 0, 0,40,eemjimmmjmijmjjjiimijiimejeieTHTTThhhhhhhhhTJTJTJ其中 为单元温度刚度矩 为单元节点温度向量同理:对边界单元可得到对总体 eHeT)3 , 2 , 1( , 0nkPTHeeennnnnnnnPPPTTThhhhhhhhh2121212222111211,PTHnkTJTJekek3 , 2 , 1, 0第三部分 热弹性力学v基本概念v应力(材料力学定义)vP分解得到垂直与截面的分量称为正应力，相切于截面的分量称为剪应力。单位为v变形:

26、其中 为变形量，沿x方向的应变v剪应变，角应变v胡克定律： 或 （横向应变 ）v 或APPAlim02mNxSxlim0SEEuGrGr广义胡克定律：其中E为拉压弹性模量G为剪切弹性模量 为横向收缩系数（泊松系数）考虑温度应变影响：GrGrGrEEEzxxzyzyzxyxyyxzzzzyyzyxx,)(1)(1)(1)1 (2EG有 其中 为物体线膨胀系数,则有 TTxxTzyx000TETETEyxzzxzyyzyxx)(1)(1)(1热应力和热弹性概念P71平面热弹性问题平面热弹性平衡方程 00YxyXyxxyyyxxyxxy平面热弹性运动方程 (求解动态热应力)边界条件 2222tYxytXyxxyyyxxlmQmlQxyyyyxxx平面应力问题xyxyxxyyyyyxxxrEEE)1 (2)()(1)()(1002002用AN