第09讲数值方法学案.ppt

有限单元法 将连续的求解域离散为一组单元的组合体，用在每个单元内假设的近似函数来分片的表示求解域上待求的未知场函数，近似函数通常由未知场函数及其导数在单元各节点的数值插值函数来表达。从而使一个连续的无限自由度问题变成离散的有限自由度问题。 步骤1：剖分: 将待解区域进行分割，离散成有限个元素的集合。元素（单元）的形状原则上是任意的。二维问题一般采用三角形单元或矩形单元，三维空间可采用四面体或多面体等。每个单元的顶点称为节点（或结点）。 步骤2：单元分析: 进行分片插值，即将分割单元中任意点的未知函数用该分割单元中形状函数及离散网格点上的函数值展开，即建立一个线性插值函数。 步骤3：求解近似变分方程 用有限个单元将连续体离散化，通过对有限个单元作分片插值求解各种力学、物理问题的一种数值方法。 为了适用于多维问题计算以及处理带有相变区间的情况, 引入焓H 、熵S 两个物理量 在整个区域用统一的能量方程, 这样式(1)和式(2)可写成一个焓形式的方程 等效热导率和等效热容 利用加权余量法,可导出相变传热温度场的有限元列式 式中,M为热容阵, K 为热传导阵, T 为节点温度向量, F 为等效右端项. 本文采用θ差分法求解方程式, 差分格式为 * 当定义 式中的函数 用 这种形式来代替，在相变温度点的半温度间隔内认为是有限大的，但不等于零。可以得到光滑的热容方程 类似的得到光滑的导热系数 * 问题的最终公式表示如下 Bonacina 等认为在 间隔内有效热容和导热系数分别为 上式用显式有限元公式很容易求解。 虽然表观热容公式具有求解方程简单，能预测相变界面位置和相应的等温线，但在没有发生相变的区域出现温度曲线扭曲和频繁的非正常无物理意义的波动。 Bonacina, Katayama ,Hattori, 和 Goodrich等给出了一些改进的计算方法。 * Schneider 、Raw 和 Raw 、 Schneider 分别用修正过的焓法模型计算了一维和多维的包含相变的热传导问题。 修改包括使用非线性迭代，它基于两个简单的规则即强化了熔解前锋的合理物理现象和在新的时间步长内发现相的正确分布。而且，和先前的焓法模型相比，CPU的时间少两个数量级。 Schneider用该方法计算了矩形封闭腔体内充满PCM材料，在具有自然对流的条件下，相变的能量迁移问题。 * 坐标变换法 求解具有自然对流的相变导热问题 设高为H的管子垂直地穿过处于溶解温度Ts的固体区域。在t=0时刻壁温升高为TwTs，于是周围的固体开始熔化。确定熔化区域随时间的变化和温度分布。 * * 熔化区内任一点的密度 有效压力 控制方程的守恒型 * 假定： 相界面上所传递的热量与界面的移动之间，允许有一个时间步长的滞后，即熔化区内，在一个时间步长的间隔内的流动与热量传递是在固定的相界面下计算的。 假定相界面的倾角很小，可认为平行于x轴。 长为dx的一段微元距离内界面上的热平衡 以上方程构成完整控制方程 初始条件，t=0，液体层厚度为零，为避免出现分母为零，实际计算中从一个很小的均匀厚度开始。形成这一厚度的时间用Stefan精确解得到. * 为了保持液体区域位于0与1之间，作坐标变换： 通过变换，把物理平面上的带运动边界的圆柱坐标系中的计算区域，变成计算平面上不动边界的矩形域。 为了把计算时间控制在合理的范围内，而又不严重影响计算结果，假定x方向的导数均忽略不计。 * 定义无量纲量 * 方程组转换为 * 经过变换，在ξ-η坐标系内液体区的计算范围已经简化成一个单位宽度的矩形，但控制方程极为复杂。由于计算从某一初始厚度开始，避免出现Ri=1的情形。 固液界面上的两相的密度差别不考虑，则边界条件为： * 采用控制容积积分法来建立离散方程。 在每一个时间间隔内，Ri为定值因而在离散方程中是常数。 在每一个时层上，这一方程组的求解可采用SIMPLE系列的算法，而在时间坐标上采用隐式格式，又选定的初始时刻向前推进。 在每一个时层的终了获得的收敛值后，可计算由于这段时间间隔内界面上的换热而造成的熔化界区的扩大。 * 计算中取 在r方向上均分，在X方向上逐密布置，顶部最密。 * SIMPLE算法 计算步骤 假定一个速度分布，记为U，以此计算动量离散方程中的系数及常数项。 假定一个压力场。 以此求解两个动量方程. 求解压力修正值方程。 据此改进速度值。 利用改进后的速度场求解那些源项物性等与速度场耦合的变量。 利用改进后的速度场重新