02 材料科学与工程中典型物理场的数值模拟

第二章材料科学与工程中典型物理场的数值模拟,2.1温度场的计算,材料科学与工程领域存在大量传热问题金属材料的热加工高分子材料的成型陶瓷材料的成型、烧结,热量从一个物体传至另一个物体，或由同一物体的这一部分传至另一部分的过程称为传热或换热。只要两物体间或同一物体内部存在温度差，就会发生热量的传递，而且热量总是从高温部分向低温部分传递。,按物体温度是否随时间的变化而变化，热量传递过程可分为稳态过程（又称定常过程）与非稳态过程（又称非定常过程）两大类。凡是物体中各点温度不随时间改变而变化的热传递过程称为稳态热传递过程，反之称为非稳态热传递过程。,物体传热的三种基本形式导热对流辐射,物体各部分之间不发生相对位移，依靠分子、原子及自由电子等微观粒子的热运动进行的热量传递称为导热。导热现象的规律用傅立叶定律描述。,一.导热方程,1.傅立叶（Fourier）定律如图3-1所示的两个表面分别维持均匀恒定温度的平板，对于x方向上任意一个厚度的微元层来说，单位时间内通过该层的导热热量Q（W）与当地的温度变化率及平板面积A成正比，即:,图1通过平板的一维导热,2.热流密度（单位面积热流量）单位时间内通过单位面积的热量称为热流密度，记为q（W/m2）。,（3-2）,式中是比例系数，称为热导率(W/(mk)，又称导热系数，t/x是x方向上的温度梯度（K/m）；负号表示热量传递的方向与温度升高的方向相反。,3.对流,对流是指流体各部分之间发生相对位移，冷热流体相互掺混所引起的热量传递方式。对流仅发生在流体中，而且必然伴随着导热。,牛顿对流换热计算公式：对流换热所传递的热流量与流体和固体表面间的温度差以及两者的接触面积成正比。对流换热的基本计算式是牛顿冷却公式：流体被加热时：q（tw-tf）流体被冷却时：q（tf-tw）式中，称为对流换热系数，简称换热系数,4.热辐射,物体通过电磁波来传递能量的方式称为辐射。由于热的原因发出辐射能的现象称为热辐射。自然界中各个物体都不停地向空间发出热辐射，同时又不断地吸收其他物体发出的热辐射。辐射与吸收过程的综合就造成了以辐射方式进行的物体间能量传递辐射换热。,黑体单位时间内发射出的辐射热流量用斯蒂芬一玻尔兹曼定律表示：QF0T4T为黑体的绝对温度；0为黑体的辐射常数。实际物体辐射热流量的计算采用斯蒂芬玻尔兹曼定律的经验修正形式：QF0T4称为该物体的黑度（又称发射率）。,二、导热微分方程,对于一维导热问题，直接对傅立叶定律的表达式进行积分就可获得答案。多维导热问题的情况则较为复杂，虽然傅立叶定律仍然可以应用，但是还必须解决不同坐标方向间导热公式的互相联系问题。这时，导热问题的数学描写必须针对从物体中分割出来的微元平行六面体进行分析才能得到。这种数学描写称为导热微分方程式。它的建立，除了依靠傅立叶定律之外，还以能量守恒定律为基础。,导热的微分方程的建立,为了便于推导，把讨论对象局限于常物性（即物性，Cp都是常量）的各向同性的材料。对如下图所示的微元平行六面体进行导热分析，按照能量守恒定律，微元体的热平衡式可以表示为：（导入微元体的总热流量）（微元体内热源的生成热）（微元体内能的增量）（导出微元体的总热流量）,导入及导出微元体的总热流量可用傅立叶定律推出。任意方向的热流量总可以分解成为三个坐标轴方向的分热流量，上图中标出了微元的分热流量。根据傅立叶定律，通过x，y，z三个表面导入微元体的热量可直接写出：,同理，通过x+dx，y+dy，z+dz导出微元的热流量为：,微元体内的增量为：微元体内的生成热为：由上述的式子可以得到：,一般的三维导热微分方程（非稳态）,其中：材料的密度（kg/cm3）；c为材料的比热容（J/(kgK)）；t为时间（s）；x，y，z分别是材料沿x，y，z方向的热导率（W/(mK)）；q是物体内部的热源密度,等式左边：体元升温需要的热量；等式右边：前三项分别是x，y，z方向流入体元的热量；最后一项是体元内热源产生的热量微分方程的物理意义：体元升温所需的热量应等于流入体元的热量与体元内产生热量的总和。,二维热传导微分方程,三维稳态热传导微分方程,二维稳态热传导微分方程,定解条件是指使导热微分方程获得适合某一特定导热问题的求解的附加条件。,通过导热微分方程可知，求解导热问题，实际上就是对导热微分方程式的求解。预知某一导热问题的温度分布，必须给出表征该问题的附加条件。,三、初始条件与边界条件1、初始条件：所求问题的初始温度场，也就是零时刻的温度场。均匀温度场：T|t=0=T0温度场不均匀：T|t=0=T0(x,y,z),2、边界条件（1）第一类边界条件T=TwT=Tw(x,y,z，t),第一类边界条件,Tw,（2）第二类边界条件q=q=（3）第三类边界条件,l,n,T,=,qw,l,n,T,=,qw（x，y，z，t）,qw,第二类边界条件,-l,n,T,=,（TTf）,Tf，k,图34第三类边界条件,四、导热问题的数值解法,数值解法是以离散数学为基础，以计算机为工具的一种求解方法。常用的数值解法有有限差分法和有限元法。有限差分法是把原来求解物体内随空间和时间连续分布的温度问题，转化为求在时间领域内和空间领域内有限个离散点的温度值问题，再用这些离散点上的温度值去逼近连续的温度分布。它以差分代替微分，差商代替微商，建立以节点温度为未知量的线性代数方程组，然后求解得到各节点温度的近似值。,*泰勒级数展开法,根据泰勒级数展开式，用节点(m,n)的温度tm,n来表示节点(m+1,n)的温度tm+1,n用节点(m,n)的温度tm,n来表示节点(m-1,n)的温度tm-1,n,向前差商：向后差商：中心差商：,将上两式相加可得,将上式改写成的表达式，有,同样可得：,例1：以一维稳态导热为例来说明有限差分方法的应用。设有一炉墙，厚度为，已知内壁温t0500，外壁温度为t4100，求厚度方向上的温度分布。解：该问题可用下式求解：,t=c1x+c2x=0,t=c2=500 x=,t=c1+c2=100c1=100-500/=-400/即得：t=c1x+c2=（-400 x/）500因此在特定点上的值为：,采用有限差分求解：如图所示，把需要求解的空间区域0以间距x（称为步长）划分m等分，这些等分线称为网格线。,节点划分示意图,以每一网格线为中心，宽度为x组成一系列子区间，称为元体，如图中阴影所示。元体的中心点称为节点。以节点的温度表示元体的温度。,为了计算节点温度，需要建立节点的差分方程。将差分方程代入得：ti+1-2ti+ti-1=0,设m4，则：t0500t22t1+t00t32t2+t10t42t3+t20t40利用高斯赛德尔法可求解得到各节点温度值。,平面温度场的有限差分求解,x,y,Tf，k,qw,绝热,Tw,（1）划分网格根据求解区域的形状将连续的求解域离散为不连续的点，形成离散网格。,i,j+1,i-1,j,i,j-1,i+1,j,i,j,L1,L2,（2）差分方程的建立,对流传热边界条件：x0,0yL2,热流边界条件：y=0，0xL1,绝热边界条件:x=L1,0yL2,给定温度边界条件:y=L2,0xL1,T=Tw,热传导微分方程：,根据导热问题的控制方程(导热微分方程),若x=y则有,得,边界条件的差分格式：1）对流传热边界条件2）热流边界条件：3）绝热边界条件：4）给定温度边界条件：,qw,Ti,j-Ti1,j,=0,Ti,jTw,方程组：,例2：设有如图所示的温度区域，边界点温度已知，且区域内无内热源，求内部节点温度。已知T1T2T3T4100，T5860，T6400，T7300，T8200。,解：根据,可以得到：,对于节点A：T1T2TBT84TA=0对于节点B：TAT3TCT74TB=0对于节点C：TBT4T5T64TC=0,雅可比叠代法：,高斯赛德尔迭代法：,非稳态导热问题,从前面的介绍可以看出，有限差分数值解法的实质在于使一个连续体离散化，用一系列的代数方程式代替微分方程式，通过对一系列代数方程式的四则运算来获得温度场的近似数值解。非稳态导热的特点是温度不仅随空间坐标的变化而变化，而且还随时间的变化而变化，因此，离散化处理时必须同时将所研究的空间和时间范围各自等分成许多细小的间隔。,考察一个无内热源，导热系数为常数，宽、高无限，厚为L的平板的一维非稳态导热问题。当无内热源，cp均为常数时，导热微分方程为：,初始条件：T(x,0)=(x)(t=0,00,0xL),建立差分方程时，首先将空间和时间范围离散化。时间与空间坐标系如图所示，横坐标为空间距离，纵坐标为时间。用竖直网格线将空间区域按x划分成若干子区域，用水平网格线将时间步长按划分为若干时层。任一节点（i，n）上的温度值用tin。x称为距离步长，称为时间步长。建立导热问题的差分方程时，导热微分方程左边各项可以按上一时间层计算，也可以采用新的时间层计算，得到不同的差分格式。,显示差分格式：,当初始条件、边界条件已知时，就可以求各节点的温度。,隐式差分格式,Crank-Nicoson差分格式,有限元法求解,温度场有限元法计算的基本方程以二维问题为例，说明建立稳态热传导问题的有限元格式的过程。稳态二维温度场的微分方程为,（1）构造试探函数构造试探函数，并设满足s边界上的强制边界条件。,将试探函数代入上述方程式及相应边界条件式，因为的近似性，将产生余量，即有,式中，nx，ny是边界外法线的方向余弦。用加权余量法建立有限元格式的基本思想是使余量的加权计分为零，即,将相应余量代入上述方程并进行分部积分，可以得到：,（2）离散空间域将空间域离散为有限个单元体，在典型单元内各点的温度T可以近似用单元的节点温度Ti插值得到,式中，ne是每个单元的节点个数；是插值函数，它具有下述性质,（3）选择权函数用伽辽金法选择权函数；在边界上选择,因为已满足强制边界条件，因此在1边界上不再产生余量，可令w1在1边界上为零。综合上述各式可得：,上式写成矩阵形式有,上式是n个线性方程组成的方程组，求解即可以得到n个节点的温度Ti，简记为KT=P。式中K称为热传导矩阵；P是温度载荷矩阵，是节点温度列阵。,平面温度场的有限元法求解,（1）单元划分平面三节点三角形单元要求：三角形的三条边的边长尽量接近，一般最长边的长度不大于最短边的长度的三倍；各个三角形单元之间只能以顶点相交，两相邻三角形单元可以有一条公共边或一个公共顶点。（2）插值函数试探函数对于三角形单元，通常假设每个单元e上的温度分布T是坐标的线性函数。这一假设在单元面积足够小时可以得到满足，因此插值函数表示为：T=a1+a2x+a3y,式中，a1，a2，a3为待定系数。插值函数必须满足单元三个节点，即当x=xi，y=yi时，T=Ti；当x=xj，y=yj时，T=Tj；当x=xm，y=ym时，T=Tm；代入插值函数，得到方程组,写成矩阵形式,利用矩阵求逆的方法可以得到：,其中,A为三角形的面积，并且,将求得的a1，a2，a3代入插值函数得：,通常简记为,其中,（3）单元积分计算,热传导矩阵元素,三角形单元边界上的对流传热边界条件对单元热传导矩阵的影响,式中，l表示三角形单元三条边的长度,（4）单元的总体合成将单元热传导矩阵合成为总体热传导矩阵。总体节点载荷阵列合成时，要考虑对流传热边界条件对单元热载荷的影响：,（5）求解解方程组：KT=P,有限元分析实例,如图所示为一个工业用烟囱的截面图，烟囱所用材料为混凝土，其热导率=1.4W/（mK）。假设烟囱内表面的温度恒定为100，外表面暴露在温度为30的大气中。外表面与空气之间的对流传热系数k=20W/（m2K）。用有限元方法来计算烟囱壁中的温度分布。,解：（1）划分单元并给单元和结点编号：由于问题的对称性，选取求解域的1/8来计算，如图所示，将求解部分划分为8个三角形单元，具有9个结点。（2）单元热传导矩阵：三角形单元的热传导矩阵为：,单元，具有相同的形状和尺寸，因此它们的热传导矩阵是相同的。,由,得：且三角形单元的面积,因此,单元，具有相同的形状和尺寸，因此它们的热传导矩阵是相同的，根据上述公式有：,因此,（3）边界条件由于对流产生的散热产生在单元，的jm边，因此,对流边界条件对单元，的热载荷矩阵的贡献为,（4）总体热传导矩阵的合成通过上面的分析，可以得到各个单元的单元热传导矩阵为：,边界条件为：,按照单元矩阵在总体矩阵中的位置，将相应的元素放在总体热传导矩阵的相应位置，即得到总体热传导矩阵：,将结点1和2上的给定温度边界条件加到热传导矩阵中可得到：,同样可得到热载荷矩阵：,加上结点1和2的给定温度边界条件，最终热载荷矩阵为：,（5）求解方程组由以上可得结点方程为：,=,解此线性方程组，即得到各结点的温度值：,2.2浓度场计算,在固体中，扩散是物质传输的唯一方式。扩散与材料在生产和使用过程中的许多重要的物理化学过程密切相关。,扩散控制方程,1、Fick第一定律在稳态扩散的条件下，单位时间内通过垂直于扩散方向的单位截面积的扩散的物质的通量J与浓度梯度成正比，其数学表达式为：,2、Fick第二定律,数值解法,一维问题的微分方程为：,采用Crank-Nicoson差分格式,上述方程组可以表示成矩阵形式，解此方程组即可得出任一时刻渗碳层碳浓度分布。,2.3铸件充型过程的通用数学模型,连续性方程：,动量守恒方程：,体积函数方程：,能量守恒方程：,式中：u、v、w分别为x、y、z方向速度分量；为金属液密度；t为时间；P为金属液体内压力；为金属液动力粘度；gx、gy、gz为x、y、z方向重力加速度；F为体积函数，；c为金属液比定压热容；T为金属液