金属凝固计算机模拟b.ppt

1 金属凝固计算机模拟 前言温度场数值模拟浓度场数值模拟 北京科技大学材料学院吴春京 尼加拉瓜瀑布2000 05 18摄 前言 材料加工过程铸造 锻造 焊接 挤 压 拉 拔 热处理 专业设置 大学 研究生 专业学位 2 计算机在材料加工中应用领域有五个方面 1 科学计算 数值模拟 优化计算如最低成本配料优化计算 其是应用线性规划解法 求满足化学成分要求范围内 在已知化学成分和价格的几种或几百种原料中 同时还可考虑某种原料最多 最少或必须用多少的条件 求出最低成本的配料单来 也可以用此方法求出最高成本的配料单 传统手算配料单 由于其不考虑成本 也无法考虑成本 因此 传统手算的配料成本落在最高和最低成本之间 3 4 5 6 7 8 9 前言 2 信息处理信息处理的特点是计算简单 数据量大 现在75 85 的计算机用于信息处理 数据库数据库 DATABASE 的概念由美国军方于1963年提出 生产 经营 人事 工资等管理都可应用数据库管理 材料加工的生产 经营 人事 工资等管理都可应用数据库管理 10 前言 专家系统第一个专家系统是1965年由美国斯坦福大学开发的 用质谱仪得到的数据来决定一个未知化合物分子结构的程序 专家系统是将专家的知识和经验按一定规则存在知识库或数据库里 用户则可通过专家系统 查寻检索到存在专家系统内的专家的知识和经验 材料加工的计算机专家系统正在研究之中 如 铸造缺陷分析专家系统 辊型设计的专家系统 11 12 前言 斯坦福大学 StanfordUniversity 创建于1885年 位于加州圣克拉拉县 总面积超过8000英亩 校园约5700英亩 住宅约500英亩 剩下大半个校园 除了一些缓冲地段之外 就是植物园 高尔夫球场和若干个科学实验场 校园构建誉为美国城市规划和庭院设计的典范 紧凑 集中 和谐 主要建筑 土黄色石墙环绕下的红屋顶建筑 拱廊相接 棕榈成行 最有名的建筑是斯坦福纪念教堂 校园内公路长达46英里 三个水坝和湖 1英亩 6 070市亩 1英里 1 609公里 13 前言 斯坦福的工业园是闻名遐迩的 硅谷 诞生地 在64名硅谷最重要的先驱中 有28人是斯坦福大学的校友 教授或研究人员 占42 7个学院 地球科学 工 文理 商 教育 法 医 专职教师1714人 从建校到现在已有25位诺贝尔奖得主 现有教师中 有17位诺贝尔奖得主 4位普利策奖得主 23位麦克阿瑟奖得主 21位国家科学奖得主 3位国家技术奖得主 220位美国文理科学院院士 128位国家科学院院士 83位国家工程院院士 25位国家教育科学院院士 42位美国哲学社会科学学会成员 6位Wolf基金奖获得者 6位Koret基金奖获得者 2位总统自由奖获得者 14 前言 3 计算机辅助设计和制造 计算机辅助设计 简称CAD 计算机辅助制造 简称CAM 近年来由于在计算机上进行设计 就连设计图的概念也有所改变 如在波音767飞机设计中 图纸定义为点列的集合 然后从美国把录有这些点列数据的磁盘寄到日本等国即可生产 输出设计图只是为了便于人的理解 而并非照图加工 机械加工时全部使用数控加工机械 材料加工工艺设计的CAD 材料加工工装制造的CAM 材料加工车间设计的CAD系统正在研究发展之中 由于材料加工的复杂性 材料加工车间的多样性和研究发展水平限制 现在材料加工的CAD CAM只局限于某一类型的CAD CAM 15 前言 现代英国铸造的CAD CAM技术发展到每周能够出图18000多幅 日本日立公司研究的设计金属压铸金属型CAD系统 能够在获得产品数据之后 启动画出型腔图 并对模型进行计算 输出用于CAM的数据磁盘 4 计算机测试和控制系统计算机和检测仪表 控制部件结合 可形成计算机测试 控制系统 材料加工生产中的测试和控制问题 一般地说 都可以采用计算机测试 控制系统 5 机械手和机器人 16 前言 数值模拟两个不同领域的现象 能用同一数学形式来描述 称这两个现象彼此是可模拟 模拟的方法是把一个领域内求解的问题过渡到另一个领域中去解决 随着计算机和计算数学的发展 用计算机数值计算法求解问题的计算精度已经达到很接近解析解 此法称为计算机数值模拟 仿真的定义 与模拟的差异 工艺优化计算机数值模拟的最终目的是解决工艺优化设计问题 一旦全面实现了材料加工过程的计算机数值模拟 材料的加工生产将会产生深刻的变革 17 前言 计算机数值模拟的步骤 给定研究对象几何条件 物理条件 初始条件 边界条件 离散化处理将过程所涉及的区域在空间和时间上进行离散化处理 建立数值方程建立内节点和边界节点的数值方程 选择计算方法选择适当的计算方法求解线性代数方程组 编制计算机程序编制计算机程序 由计算机算出结果 并用数据 曲线 图形输出 优化工艺分析计算机模拟的结果 如果结果不理想 调整工艺参数 再进行计算机模拟 直到模拟结果为最佳结果 18 前言 讲解的主要内容通过一 二个例子 让大家亲自经历计算机数值模拟的一般步骤 从中掌握其基本原理和方法 为更好应用计算机数值模拟 优化工艺参数打下基础 材料加工过程中 比较常用 典型的数值模拟有 温度场数值模拟 浓度场数值模拟 教学安排和考核 讲课18学时 自己上机练习4学时 考试 笔试开卷 2学时 成绩评定 笔试开卷 70 交上机的程序及运行结果 30 19 20 1温度场计算机数值模拟 1 1传热的基本知识1 1 1传热的基本方式 导热导热属于接触传热 是连续介质就地传递热量 没有各部部分物质之间宏观的相对位移 在不透明固体实体内部 由于各部分物质之间无法作宏观的相对位移 不透明无法传递辐射能 实体保证接触 所以只能依靠导热方式传递热量 导热的基本定律是傅立叶定律 即导热的热流密度q与温度梯度成正比 即 q gradT gradT W m2 1 1 1 n式中 q 传热方向上单位时间 单位面积上所通过的热量 J s m2 W m2 材料的导热系数 W m K 温度梯度 K m n负号表示导热的热流量向温度低的方向传递 即与温度梯度的方向相反 热流密度是个向量 即它有大小和方向 21 1 1 1传热的基本方式 如果热流密度的分量和 X Y Z 坐标系相联系 那么X Y Z方向的热流密度分量应是 qx qy qz 1 1 2 X Y Z热流密度q iqx jqy kqz1 1 3导热系数物理意义 沿导热方向的单位长度上 温度减低 物质所容许通过的热流量 方向性 大多数液体和固体属于各向同性的物质 各向异性材料的导热系数具有方向性 如石墨 温度函数 值还随温度而变化 大多数金属的导热系数随温度的升高而降低 大多数液体 水和甘油除外 其导热系数随温度的升高而降低 气体的导热系数随温度的升高而增加 22 1 1 1传热的基本方式 对流对流是流体 气体和液体 中温度不同的各部分相互混合的宏观运动引起热量传递的现象 工程上最具有实际意义的是 相对运动着的流体与所接触的固体壁面之间的热量交换过程 一般称为对流换热 工程上在研究固体壁面和流体之间的对流换热时 除了高度稀薄的气体外 人们不去注意流体的单个质点 而把流体看成是连续介质 实际的流体总有粘性 流动时 受粘性和壁面摩擦的影响 在靠近壁面附近的流体将降低流速 在壁面上完全被滞止不动 即X 0时 0 如图1 1所示 因此 热量从壁面传给贴壁的那部分流体 将依靠导热 T K V m s TwqwcTfVX m 图1 1邻近壁面的流体速度分布和温度分布 23 1 1 1传热的基本方式 流体和固体壁接触面上的 相 际热流密度为 qwc f x 01 1 4 X式中 qwc 热流密度 W m2 f 流体的导热系数 W m K T 流体的温度 K 液体的导热系数较小 气体的导热系数更小 所以受热时 在贴壁处的流体温度势必沿着 轴的反方向急剧升高 随着离壁面的距离 的增加 流体的流动将带走更多的热量 使 轴方向的温度梯度连续下降 温度分布趋向平坦化 正是通过这种导热和对流的共同作用 使热量在流体内部得到传播 越临近壁面 导热越起主导作用 图1 1所示 假如厚度为 的贴壁静止膜 膜内温度线性地从壁面温度TW降到远离壁面 尚未被加热的流体温度Tf 则Tf TWqwc f 1 1 5 无界对流时壁面与流体的换热 钢锭与周围空气的对流换热属于这种情况 24 1 1 1传热的基本方式 流体在管和槽道内部的流动 称为有界对流 这时不存在远离壁面 尚未被加热的流体温度 则采用截面平均温度作为流体的特征温度Tf 则qWC aC TW Tf W m21 1 6这就是所谓的 牛顿冷却定律 式中 aC为放热系数 W m2 K 其实 牛顿冷却定律 并不是表达对流换热现象本质的物理定律 凡能影响流体流动的各种因素 包括流体的种类和状态 运动的速度 流道的形状和大小 以及固体壁面粗糙度等 都会对aC值产生影响 式1 1 6只不过形式地把放热过程的一切复杂性和计算上的困难 都转移到并且集中在放热系数这样一个物理量上罢了 aC代表流体和所接触的固体表面之间温度每相差 该流体与表面之间 相 际热流量的大小 运用式1 1 6可以进行对流换热的计算 但由于对流换热的复杂性 该式中的放热系数aC需从相应的准则方程式求出 准则方程式是针对不同对流换热情况 在综合实验结果的基础之上 运用相似理论将表征某现象的物理量整理成一些相似准则 用因次分析法得到的各种类型的表达式 连铸钢坯二冷区 25 辐射只要温度高于绝对零度 任何物体都随时向周围空间辐射能量 辐射热流密度用斯蒂芬 玻耳兹曼定律表达 E w m21 1 7式中 物体的辐射率或黑度 0 1 斯蒂芬 玻耳兹曼常数或绝对黑度的辐射常数 5 67 10 8W m2K4 温度 K 实际上 辐射往往涉及二个物体间辐射热交换 如果二个物体的表面温度不同 中间被空气所隔开 T1 T2时 则相互辐射的结果 表面温度为T1的物体发射出去的辐射热超过了吸收来自表面温度为T2的物体辐射热 引起辐射换热的热流量则为 Q1 2 12 1 2 F1 12或q1 2 12 1 2 121 1 8式中 12 物体1与2综合黑度 F1 物体1的表面积 12 物体 的表面向外发射出去的辐射热量中 能投射到物体 表面上的份额 称为角系数 0 1 1 1 1传热的基本方式 26 轧辊喷淬时的工作环境 27 喷淬机的简易俯视图 28 红外测温仪 吹扫装置1 红外测温仪 2 固定支架 3 瞄准管 4 进气管 5 瞄准管前端 29 辊身不同深度处的冷却曲线 30 喷雾过程辊身处表面综合换热系数与时间的关系 31 1 1 1传热的基本方式 壁面在气体中冷却 存在辐射换热和对流换热 考虑到壁面与气体之间存在着辐射换热 其热流量 或热流密度 为Qr arF Tw Tf w或qwr ar Tw Tf w m21 1 9式中 qwr 单位时间 单位面积辐射的热量 辐射热流密度 w m2 ar 辐射放热系数 w m2 k Tw 辐射物体表面温度 k Tf 透明的气体介质的特征温度 k 考虑到壁面与气体之间还存在着对流换热 其热流密度为qwc ac Tw Tf w m21 1 10由壁面传走的总热流密度qw应是qwr和qwc二者之和qw ac Tw Tf ar Tw Tf 1 1 11令a0 ac ar则qw a0 Tw Tf w m21 1 12用辐射放热系数ar 可以形式地把辐射换热折合成对流换热 用总放热系 32 1 1 1传热的基本方式 数a0兼顾辐射换热的影响 从而有利于简化复杂传热的分析和计算 如远离表面的外界表面温度趋于环境温度Tf 并且 12 1时 由式1 1 8得qwr 12 w f 1 1 13由式1 1 9得 qwr ar w f 1 1 14由式1 1 13和式1 1 14得 12 w f ar w f 1 1 15因 w f w f w3 w2 f w f2 f3 设 m 1 2 w f w3 w2 f w f2 f3 4 m3ar 4 12 m3 w m2 k 1 1 16从此可见 w f 降为零时 ar并非零值 而是以4 12 m3 随着温差的扩大和平均温度 m的升高 ar值将迅速增加 由于ac随温差的变化较小 在高温范围和大温差情况下 ar有可能成为a0的主要组成部分 ar与气体的运动状况无关 而ac与气体速度的降低而减小 在气体自然对流的情况下 ac 5 w m2 k 即使 m只有300K ar就可能占a0的一半左右 33 1 1 2导热的偏微分方程式 钢锭一般属于各向同性的物质 其加热或冷却过程数学模拟计算依据的基本数学模型 是不稳定导热偏微分方程 下面讨论各向同性材料导热方程式的建立 直角坐标系导热偏微分方程导热偏微分方程的建立 是通过考察处于导热过程中的物质的微元体积 x y z 的能量平衡来建立 如图1 2所示 在时间 内 通过六个面的导热所获得的能量 加上微元体内产生的内热源热量 要等于微元体积内物质积蓄热量的改变 即温度升高或降低 x z y dQx x dQx dQz dQy dQz z dQy y 图1 2直角坐标系导热方程式的微元体 34 1 1 2导热的偏微分方程式 图1 2中 微元控制体尺寸 x y z 按照傅立叶导热定律 在X方向流入微元体左表面的热流可表示为 TdQx y z W1 1 17 X式中 导热系数 W m k y z X方向微元体表面积 m2 T X方向的温度梯度 k m X在X方向流出微元体右表面的热流可表示为 TdQx x y z W X在X方向导热的净热流为 dQx dQx 0 35 1 1 2导热的偏微分方程式 图1 2中 微元控制体尺寸 x y z 按照傅立叶导热定律 在X方向流入微元体左表面的热流可表示为 TdQx y z W1 1 17 X式中 导热系数 W m k y z X方向微元体表面积 m2 T X方向的温度梯度 k m X在X方向流出微元体右表面的热流 可以应用泰勒级数展开 保留级数的第一项和第二项而得到 TdQx x dQx dQx x dQx y z x X X X T dQx x y z1 1 18 X X在X方向导热的净热流为 TdQx dQx x y z1 1 19 X X 36 1 1 2导热的偏微分方程式 用同样的方法 可以得出Y Z方向与式1 19相类似的导热的净热流方程式 即 TdQy dQy y x y z1 1 20 Y Y TdQz dQz z x y z1 1 21 Z Z在三个坐标方向净热流的总和为 T T T x y z1 1 22 X X Y Y Z Z如果单位时间 单位空间所产生的热量为Q x y z 那么微元体单位时间的发热量 热流 为 Q x y z x y z1 1 23由于导热传进微元体的净热流式1 1 22和微元体内产生的热量式1 1 23一起用于增大微元体的内能 微元体的内能的增加表现在微元体能量存储随时间的变化率 即 T Cp x y z 1 1 24 37 1 1 2导热的偏微分方程式 式中 密度 kg m3 Cp 比热 J kg k 时间 s 对微元体进行能量平衡 使能量存储的时间变化率与由导热引起的流入微元体的净热流和微元体内产生的热流之和相等 得下式 T T T T Cp Q 1 1 25 X X Y Y Z Z 圆柱坐标系导热偏微分方程实际问题常常涉及柱面对称问题 且边界条件给定在一个表面上 因此表面具有一个坐标保持不变的性质 在这种情况下 采用圆柱坐标系是合适的 图1 3圆柱坐标系导热方程式的微元体 38 1 1 2导热的偏微分方程式 图1 3所示圆柱坐标系 直接按内外两个圆弧面和其它四个平面组成的微元体积 在导热过程中热量必须按收支平衡的原则导出 此时 微元体积为 dv rd dzdr1 1 26沿内圆弧面流入微元体积的热流 TdQr rd dz 1 1 27 r沿外圆弧面流出微元体积的热流 TdQr dr dQr dQr dr dQr rd dz dr r r r T1 T dQr r d dzdr dQr r dv1 1 28 r rr r r沿 平面流入微元体积的热流 TdQ drdz 1 1 29r 沿 d 平面流出微元体积的热流 TdQ d dQ dQ rd dQ drdz rd r r r 39 1 1 2导热的偏微分方程式 1 dQ dv1 1 30r2 沿z面流入微元体积的热流 TdQz rd dr 1 1 31 z沿z dz面流出微元体积的热流 TdQz dz dQz dQz dz dQz rd dr dz z z z T dQz dv1 1 32 z z根据直角坐标系导热微分方程推导的思路 可得到 1 T1 T T r Q Cp 1 1 33r r rr2 z z 40 1 1 2导热的偏微分方程式 球坐标系导热偏微分方程对于球坐标系 如图1 4所示 图1 4球坐标系导热方程式的微元体 由内 外两个球面 两个圆弧面和两个平面所组成的微元体积为 dv rd rSin d dr1 1 34 41 1 1 2导热的偏微分方程式 沿半径方向流入微元体积的热流 TdQr rSin d rd 1 1 35 r沿半径方向流出微元体积的热流 TdQr dr dQr dQr dr dQr rSin d rd dr r r r T1 T dQr r2 Sin d d dr dQr r2 dv1 1 36 r rr2 r r沿 方向流入微元体积的热流 TdQ rSin d dr 1 1 37r 沿 方向流出微元体积的热流 TdQ d dQ dQ rd dQ rSin d dr rd r r r T T dQ Sin rd dr rd dQ Sin dvr2 Sin r2 1 1 38 42 1 1 2导热的偏微分方程式 沿 方向流入微元体积的热流 TdQ dr rd 1 1 39 rSin 沿 方向流出微元体积的热流 dQ d dQ dQ rSin d rSin T dQ dr rd rSin d rSin rSin T dQ dv1 1 40 r2Sin2 同理可整理得 T T T r2 Sin Q Cp r2 r rr2Sin r2Sin 1 1 41 43 1 2导热方程的有限差分解法 求解不稳定导热偏微分方程的数值解法 主要有 有限差分解法 有限元法 边界元法三类 边界元法正在研究和完善之中 目前常用的是有限差分解法和有限元法 我们这门课专门讨论有限差分解法的数学基础 数值方程的建立 差分方程的稳定性和收敛性等问题 有限差分解法是用差分方程近似地代替微分方程 建立差分方程有直接法和能量平衡法两种 1 2 1直接代换法直接代换法是从微分形式出发 用差商直接代换微商的办法建立差分方程 1 2 1 1数学基础 微商和差商的定义若T x 是连续函数 则其导数为 dT x x T x lim lim 1 2 1dx x 0 x x 0 x 44 1 2 1 1数学基础 式1 2 1右边 x是有限的差商 x和 都不为零 式左边d dx是 T x当 x 0时的差商 称之为微商 在 x没有达到零之前 T x只是dT dx的近似 实际应用 x 0 如果把 T x趋于dT dx的过程认为是近似向精确过渡 那么 用 T x代替dT dx 则两者的差值 T x dT dx 表示差商代替微商的偏差 误差多大 需要做误差分析 才能大胆地应用 误差分析假设函数T x 在x时的值为T x 在x x时的值为T x x 如果函数T x 在x处的各阶导数存在 则按照泰勒级数展开 T x 与T x x 的关系如下式所示 dT x 2d2T x ndnTT x x T x x 1 2 2dx2 dx2n dxn整理后得 TT x x T x dT xd2T x n 1dnT 1 2 3 x xdx2 dx2n dxn从上式可知 T x 在x处的差商 T x等于函数T x 在x处的各阶导数的线性组合 只能是近似地等于差商 两者之间也必然有偏差 45 1 2 1 1数学基础 图1 2 1表示了一阶差商与一阶微商之间的关系 用不同方法得到的差商去代替微商 它们带来的误差是不同的 即向前差商 dT dx T x x T x x1 2 4向后差商 dT dx T x T x x x1 2 5中心差商 dT dx T x x T x x 2 x1 2 6 图1 2 1差商与微商 46 1 2 1 1数学基础 按照泰勒级数展开 T x 与T x x 的关系如下式所示 dT x 2d2T x ndnTT x x T x x 1 2 7dx2 dx2n dxn整理后得 T x x T x dT xd2T O x 1 2 8 xdx2 dx2即向前差商的偏差是截去了泰勒级数展开式中的高阶项而引起的 常称 截断误差 其截断误差为与 x同级的小量O x 同理dT x 2d2T x 3d3TT x x T x x 1 2 9dx2 dx23 dx3整理后得 T x T x x dT xd2T O x 1 2 10 xdx2 dx2即向后差商的截断误差为与 x同级的小量O x 47 1 2 1 1数学基础 由式1 2 7减式1 2 9 将2 xdT dx移至等式左边 两边再除以2 x 得 T x x T x x dT x 2d3T O x 2 1 2 112 xdx3 dx3即中心差商的截断误差为与 x 2同级的小量O x 2 当 x固定时 上述一阶差商一般仍是x的函数 对它们还可以求差商 这种一阶差商的差商称为二阶差商 它是二阶微商的近似 常用向前和向后差商来二阶微商 即d2TT x x T x T x T x x xdx2 x xT x x 2T x T x x 1 2 12 x 2由式1 2 7和式1 2 9相加 经简化后得 T x x 2T x T x x d2T O x 2 1 2 13 x 2dx2即二阶差商的截断误差为与 x 2同级的小量O x 2 48 1 2 1 2建立内节点差分方程 一维系统假定有一高宽无限 即高宽方向上无热流 厚度为L的平板 T f x 即温度是x方向位置和时间的函数 所谓一维系统是指几何空间为一维 初始时刻 0 T T0 为了简化 考虑无内热源 Cp均为常数 选取网格点间距 x和时间步长 将研究对象离散化 显式差分格式 T T一维不稳定导热方程为 Cp X X该方程在区域 0 0 x L内全部点都成立 如图1 2 2所示 将方程应用于内节点i可写成 Tp 2Tp Cp p 1 2 13 i X2i上式等号左端的微分式用温度对时间的一阶向前差商来近似 即 TpTp 1i Tpi O 1 2 14 i 上式等号右端的微分式用温度对空间的二阶差商来近似 即 TpTpi 1 2Tpi Tpi 1 O x 2 1 2 15 x2i x 2 49 1 2 1 2建立内节点差分方程 图1 2 2一维显式差分将式1 2 14和式1 2 15代入式1 2 13得 Tp 1i TpiTpi 1 2Tpi Tpi 1 Cp O O x 2 1 2 16 x 2若在上式中去掉O 和O x 2 整理得 Tp 1i Tpi F0 Tpi 1 2Tpi Tpi 1 1 2 17式中 F0 Cp x 2 导热速率 蓄热速率 称F0为傅立叶数 50 1 2 1 2建立内节点差分方程 F0的数值小意味着加热或冷却此物体所需要的时间长 反之 所需要的时间短 F0是一个无因次数字 当初始条件和边界条件已知时 用式1 2 17就可模拟区域内各节点随时间 增长的温度值Tpi i 2 3 N 1 p 1 2 3 隐式差分格式一维隐式差分如图1 2 3所示 将一维不稳定导热微分方程应用于内节点i 则 Tp 2Tp 1 Cp 1 2 18 i X2i式1 2 18和1 2 13相比 式的左边完全一样 温度对时间的一阶偏微商 仍用一阶向前差商来近似 而式1 2 18和1 2 13右边有所不同 式1 2 13中温度对距离的二阶偏微商是对应时刻p的 在用差商近似微商时 用p时刻的T值 而式1 2 18中 温度对距离的二阶偏微商是对应时刻p 1的 用差商近似微商时 用p 1时刻的T值 即式1 2 18相应的差分方程为 Tp 1i TpiTp 1i 1 2Tp 1i Tp 1i 1 Cp O O x 2 1 2 19 x 2若在上式中去掉O 和O x 2 整理得 Tp 1i Tpi F0 Tp 1i 1 2Tp 1i Tp 1i 1 1 2 20 51 1 2 1 2建立内节点差分方程 式中 F0 Cp x 2 比较式1 2 17和1 2 20 可以看出显式差分格式的突出优点是每个节点方程都可以独立求解 整个计算过程十分简便 但是 的选取要受到限制 有时为了满足差分格式稳定性条件 可能选得很小 使计算工作量加大 隐式差分格式的最大优点是 它对任意F0值都是稳定的 这种稳定是绝对的 即不受边界条件 x 热物参数的影响 于是 可以选的较大 计算速度加快 但是 对于节点i 只从式1 2 20不能独立求解 必须涉及Tp 1i 1 Tp 1i Tp 1i 1的联立线性代数方程组才能求解 也就是说 它含有三个未知数 时间步长每前进一步 从坐标0 x L 网格点1 i N整个区域的每个点 上述方程都要列出一次 见图1 2 3 因此 向一个新的时间步长每移动一步就必须解一个方程组 当按顺序列出这些方程时 除了要的第一点和最后一点的方程只有两个未知数外 其余每一个方程都含有三个未知数 于是方程是三对角的 对于这种情况 已有适用的求解程序 后面将讨论 关于差分方程稳定性将在后面讨论 52 1 2 1 2建立内节点差分方程 图1 2 3一维隐式差分 二维系统显式和隐式差分格式建立的方法和两种差分格式的特点在前面讨论过 对二维系统同样适用 为简略起见 在此只讨论二维系统显式差分方程的建立 为使问题简化 仍然假定热物性值为常数 考虑无内热源 首先把所讨论的区域离散化 如图1 2 4所示 53 1 2 1 2建立内节点差分方程 网线 用平行于X Y轴的直线 网线 进行空间离散化节点 网线与网线的交点步长 节点间的距离图1 2 4二维均质物体的分割 x与 y可以是不变的常量 即等步长 也可以是变量 即在区域内的不同处是不同的 即变步长 如果区域内各点处的温度梯度相差很大 则在温度变化剧烈处 网格布得密些 在温度变化不剧烈处 网格布得疏些 至于网格多少 步长取多少为宜 要根据计算精度与计算工作量等因素而定 对于0 x L1和0 y L2的矩形区域内 将二维不稳定导热方程式应用于节点 i j 可写成 Tp 2T 2Tp Cp 1 2 21 i j x2 y2i j 54 1 2 1 2建立内节点差分方程 式1 2 21左边一阶偏微商用温度对时间一阶向前差商近似 TpTp 1i j Tpi j O 1 2 22 i j 式1 2 21右边的二阶偏微商所对应的差商可表示为 TpTpi 1 j 2Tpi j Tpi 1 j O x 2 1 2 23 x2i j x 2 TpTpi j 1 2Tpi j Tpi j 1 O y 2 1 2 24 y2i j y 2当 x y较小时 忽略O O x 2 O y 2 项 当 x y时 即x y方向网格划分步长相等 将上三式代入式1 2 21则得到节点 i j 的差分方程 Tp 1i j Tpi j F0 Tpi 1 j Tpi 1 j Tpi j 1 Tpi j 1 4Tpi j 1 2 25式中 F0 Cp x 2 上述条件 隐式差分格式 则除Tpi j项外 其它项的Tp改为Tp 1 55 1 2 1 2建立内节点差分方程 三维系统和二维系统的假定相同 热物性值为常数 无内热源 设 x y z 即均匀网格间距 将三维不稳定导热方程式应用于节点 i j k 可以写成 Tp 2T 2T 2Tp Cp 1 2 26 i j k x2 y2 z2i j k即可得到节点 i j k 的近似式Tp 1i j k Tpi j kTpi 1 j k 2Tpi j k Tpi 1 j k Cp x 2Tpi j 1 k 2Tpi j k Tpi j 1 kTpi j k 1 2Tpi j k Tpi j k 1 1 2 27 y 2 z 2上式经简化可得节点 i j k 的近似式 Tp 1i j k Tpi j k F0 Tpi 1 j k Tpi 1 j k Tpi j 1 k Tpi j 1 k Tpi j k 1 Tpi j k 1 6Tpi j k 1 2 28式中 F0 Cp x 2 如果写成隐式差分格式 则除Tpi j k项外 其它项的Tp改为Tp 1 56 1 2 2能量平衡法 1 2 2能量平衡法能量平衡法是不再从微分方程出发用差商代替微商的办法建立差分方程 而是将导热的基本定律直接近似 局部地应用于围绕每个节点的各单元控制体积 用热量平衡法建立差分方程 由于它的基本思想是把能量守恒原则应用到每个单元体 因此常称为能量守恒原则的有限差分法 1 2 2 1建立内节点差分方程 一维系统高宽无限 即高宽方向上无热流 一定厚度的平板 T f x 初始时刻 0 T T0 为了简化 考虑无内热源 Cp均为常数 沿板厚方向 即热流方向 把均质物体分割为若干单元体 各单元体的端面积为F 几何步长为 x 时间步长为 用P上角标表示第几个时间步长 单元体中心为节点 节点温度代表该单元体的温度 二相邻节点间的导热面积是F 如图1 2 5所示 图1 2 5一维均质物体分割 57 1 2 2 1建立内节点差分方程 考虑内节点i 从时间等于 开始的 时间间隔内 通过截面F从 i 1 单元体流入到i单元体的热量为 Tpi 1 TpiQi 1 i F 1 2 29 x从 i 1 单元体流入到i单元体的热量为 Tpi 1 TpiQi 1 i F 1 2 30 x由于热量流入i单元体 引起i单元体的内能变化 在 时间内 i单元体的内能变化为 UTp 1i Tpi Cp F x 1 2 31 因为能量平衡或守恒 所以 Q U 即 Tp 1i TpiTpi 1 TpiTpi 1 Tpi Cp F x F F 1 2 32 x x整理可得 Tp 1i Tpi F0 Tpi 1 2Tpi Tpi 1 1 2 33式中 F0 Cp x 2 上式和直接代换法式1 2 17完全一样 58 1 2 2 1建立内节点差分方程 二维系统如果物体的形状使之只在X Y方向上有热流 Z方向上无热流 或相对于X Y方向可忽略不计 这种情况就可以当作二维系统处理 对于如图1 2 6所示的二维矩形区域 划分成有限个小单元体 图中虚线所示的 i j 单元体及周围四个相邻的单元体 用于推导内单元体差分方程 每个小单元体的几何步长为 x和 y 时间步长 为使问题简化 仍然假定热物性值为常数 考虑无内热源 图1 2 6二维系统的分割 59 1 2 2 1建立内节点差分方程 对于 i j 单元体 时间间隔内的内能变化为 UTp 1i j Tpi j Cp x y 1 1 2 34 在 时间内从周围四个相邻单元体流入 i j 单元体的热量分别为 Tpi 1 j Tpi jQi 1 j i j y 1 1 2 35 xTpi 1 j Tpi jQi 1 j i j y 1 1 2 36 xTpi j 1 Tpi jQi j 1 i j x 1 1 2 37 yTpi j 1 Tpi jQi j 1 i j x 1 1 2 38 y因为 Q U 若 x y 经整理后得 Tp 1i j Tpi j F0 Tpi 1 j Tpi 1 j Tpi j 1 Tpi j 1 4Tpi j 1 2 39式中 F0 Cp x 2 上式和直接代换法式1 2 25完全一样 三维系统同理可得和式 1 2 28 一样的节点差分方程 60 1 2 2 2能量平衡法与直接代换法的比较 1 2 2 2能量平衡法与直接代换法的比较尽管用能量平衡法与直接代换法所得到的差分方程在形式上完全一样 但它们所依赖的基础有所不同 用直接代换法得出差分方程时 要求区域内温度分布连续而光滑 因为按泰勒级数展开 要求温度T x 在x处的各阶导数都存在 连续而光滑即可导 而能量平衡法只要求每个单元体边界上热流可积 即可以算出每个单元体边界上的热量 它比温度分布连续光滑的要求显然放宽了 如果区域由导热系数截然不同的两种物质构成 并在交界面上存在接触热阻 或在交界面上存在间隙 那么在那些交界面上或间隙处 温度分布就不连续或不光滑 在这种情况下 可用能量平衡法建立差分方程 只需在单元划分时 使那些交界面正好落在单元的边界即可 如图1 2 7所示 图1 2 7异种物质接触时单元的划分图1 2 8界面温度分布示意图 61 1 2 2 2能量平衡法与直接代换法的比较 现假定二维矩形分割 交界面两侧的物体导热系数不同 讨论相邻两节点 i 1 j 节点流出 i j 节点的热量公式的建立 当异种物质接触时 除传导传热外 还要考虑界面热阻Rb 则Tpi 1 j Tp2Qi 1 j T2 i 1 j y 1 1 2 40 x 2Tp2 Tp1QT2 T1 y 1 1 2 41RbTp1 Tpi jQT1 i j i j y 1 1 2 42 x 2因为三个Q相似 令其相等 由式1 2 40和式1 2 42得 Q x 2 Tpi 1 j Tp21 2 43 y 1 i 1 jQ x 2 Tp1 Tpi j1 2 44 y 1 i jQ x 2 x 2 Tpi 1 j Tpi j Tp2 Tp1 1 2 45 y 1 i 1 j i j将式1 2 41代入上式得 62 1 2 2 2能量平衡法与直接代换法的比较 y 1Q Tpi 1 j Tpi j 1 2 46Rb x 2 i 1 j x 2 i j 与直接代换法相比 能量平衡法的主要优点是 在二维系统中能使用三角形 四角形 在三维系统中能使用四 五 六面体等多种分割单元 因此能解复杂形状的问题 对于每个节点单元能指定热物性值 易于求解由多种物质组成的系统的问题 能量平衡法的缺点是 输入数据量多 程序复杂 计算时间较长 所以 对于矩形或圆柱形等简单形状 能进行有规律分割的系统 采用直接代换法是有利的 可是 能有规律分割时 能量平衡法的程序与直接代换法的程序相同 一般说 能量平衡法更实用些 1 2 3边界节点差分方程的建立物体内部的温度场必然受到物体表面条件的影响 如表面散热条件差 表面温度高 物体内部的温度场平缓 反之 物体内部温度场的变化也影响着表面上的条件 如物体内部的导热系数大 物体内部温度场平缓 使表面温度升高 表面散热量也大 为了数值模拟 必须建立边界节点的差分方程 在此 先讨论一般换热条件的边界节点差分方程的建立 63 1 2 3边界节点差分方程的建立 图1 2 9所示了绝热 给定热流密度qr 对流 定温 辐射换热等五种边界表面的二维矩形区域网格 现用能量平衡法建立各种换热边界条件下的节点方程 图1 2 9各种换热边界表面图1 2 10绝热边界单元1 2 3 1绝热边界把图1 2 9中AB面上任一节点 i j 及其相邻节点取出分析 如图1 2 10所示 Tpi 1 j Tpi jTpi j 1 Tpi j y 1 x 2 1 x yTpi j 1 Tpi jTp 1i j Tpi j x 2 1 0 Cp x 2 y 1 1 2 47 y 64 1 2 3 1绝热边界 若 x y 经整理后得 Tp 1i j Tpi j F0 2Tpi 1 j Tpi j 1 Tpi j 1 4Tpi j 1 2 48式中 F0 Cp x 2 1 2 3 2给定热流密度qr边界把图1 2 9中BC面上任一节点 i j 及其相邻节点取出分析 如图1 2 11所示 用能量平衡法建立给定热流密度qr边界节点方程 Tpi 1 j Tpi j y 2 1 xTpi 1 j Tpi j y 2 1 xTpi j 1 Tpi j x 1 qr x 1 yTp 1i j Tpi j Cp x y 2 1 1 2 49 若 x y 经整理后得 图1 2 11给定热流密度边界单元Tp 1i j Tpi j F0 Tpi 1 j Tpi 1 j 2Tpi j 1 4Tpi j 2qr Cp x 1 2 50式中 F0 Cp x 2 65 1 2 3 3对流边界 1 2 3 3对流边界把图1 2 9中CD面上任一节点 i j 及其相邻节点取出分析 如图1 2 12所示 在CD面上 ac为对流换热系数 Tf为物质周围介质温度 用能量平衡法建立对流边界节点方程 Tpi 1 j Tpi j y 1 xTpi j 1 Tpi j x 2 1 yTpi j 1 Tpi j x 2 1 ac y 1 Tf Tpi j yTp 1i j Tpi j Cp x 2 y 1 1 2 51 若 x y 经整理后得 图1 2 12对流边界单元Tp 1i j Tpi j F0 2Tpi 1 j Tpi j 1 Tpi j 1 4Tpi j 2ac Tf Tpi j Cp x 1 2 52式中 F0 Cp x 2 1 2 3 4给定温度边界图1 2 9中AD面上给定温度Tw 把AD面上任一节点 i j 及其相邻节点取出分析 如图1 2 13所示 对这些单元直接写出结果 Tp 1i j Tw1 2 53 66 1 2 3 5辐射边界 把图1 2 9中AD面上受辐射热流的作用 把AD面上任一节点 i j 及其相邻节点取出分析 如图1 2 13所示 为黑度 1 为斯蒂芬 玻耳兹曼常数 Tf为物质周围介质温度 Tpi 1 j Tpi j y 2 1 xTpi 1 j Tpi j y 2 1 xTpi j 1 Tpi j x 1 y x 1 Tf4 Tpi j 4 Tp 1i j Tpi j Cp x y 2 1 1 2 54 若 x y 经整理后得 图1 2 13定温或辐射边界单元Tp 1i j Tpi j F0 Tpi 1 j Tpi 1 j 2Tpi j 1 4Tpi j 2 Tf4 Tpi j 4 Cp x 1 2 55式中 F0 Cp x 2 67 1 2 3 6混合边界 图1 2 9中四个角上的边界单元 除AD边定温边界条件时 A D两个拐角节点的温度为一定值 Tw 的情况以外 在其余换热边界条件下 它们的两个外边界从属于两种边界条件 如 C拐角节点 从图1 9可知 它的一个边界处于对流换热 另一个边界处于给定热流密度qr 在这种情况下 将 C拐角节点及其相邻节点取出分析 如图1 14所示 按能量平衡法可得 Tpi 1 j Tpi j y 2 1 xTpi j 1 Tpi j x 2 1 yqr x 2 1 ac y 2 1 Tf Tpi j Tp 1i j Tpi j Cp x 2 y 2 1 1 2 56图1 2 14混合边界单元 若 x y 经整理后得 Tp 1i j Tpi j F0 2Tpi 1 j 2Tpi j 1 4Tpi j 2qrF0 x 2acF0 x Tf Tpi j 式中 F0 Cp x 2 1 2 57 68 1 2 3 6混合边界 对图1 2 9中 A B D等拐角节点 在绝热 给定热流密度qr 对流和辐射换热等混合边界条件下 都可以按上述类似方法 分别建立它们的边界单元差分方程 对于物体出现内拐点 凹角 的情况 也可按能量守恒定律得出凹角节点的差分方程 如图1 2 15所示的凹角节点 i j 若凹角两边处于对流换热边界条件时 把 i j 节点及其相邻节点取出分析 则可得 Tpi 1 j Tpi j y 1 xTpi 1 j Tpi j y 2 1 xTpi j 1 Tpi j x 2 1 yTpi j 1 Tpi j x 1 yac x 2 y 2 1 Tf Tpi j Tp 1i j Tpi j Cp 3 4 x y 1 1 2 58图1 2 15凹角边界单元 69 1 2 4圆柱坐标系的有限差分方程 圆柱坐标系更容易表示圆柱体工件的形状 如果物体无内热源 热物性值为常数 圆柱坐标系的一般导热方程式如下 T1 T1 T Cp r r r rr2 z z1 r T1 T1 T r r r rr r rr2 z z1 T 2T1 2 T 1 2 59r r r2r2 2 z 上式的有限差分方程 可以用直接代换法和能量平衡法导出 直接代换法图1 2 16表示圆柱坐标的几何图形 z r 分别具有相等的增量 对于三维系统来说 直接代替偏微分方程之后 则可以得到节点 的显式和隐式的有限差分方程 Tp 10 Tp0Tp2 Tp4Tp2 2Tp0 Tp4Tp1 2Tp0 Tp3Tp5 2Tp0 Tp6 Cp 1 2 60 r0 r r 2r02 2 z 如果上式右端各项都以新的时间层p 1时刻计算 则得到隐式的差分方程 70 1 2 4圆柱坐标系的有限差分方程 图1 2 16圆柱坐标 能量平衡法在 时间内 从周围相邻各点流入 点的热量为 Tp1 Tp0Q1 0 z r 1 2 61r0 Tp3 Tp0Q3 0 z r 1 2 62r0 Tp2 Tp0Q2 0 r0 r 2 z 1 2 63 r 71 1 2 4圆柱坐标系的有限差分方程 Tp4 Tp0Q4 0 r0 r 2 z 1 2 64 rTp5 Tp0Q5 0 r0 r 1 2 65 zTp6 Tp0Q6 0 r0 r 1 2 66 z UTp 10 Tp0 Q Cp r0 r z 1 2 67 经整理后 可以得出节点0的显式有限差分方程 Tp 10 Tp0Tp2 Tp4Tp2 2Tp0 Tp4Tp1 2Tp0 Tp3Tp5 2Tp0 Tp6 Cp 1 2 68 r0 r r 2r02 2 z 式1 2 68和式1 2 60完全一样 除了对称轴上的一点外 该式对其它任一节点都是适用的 因为在轴上r0 0 分母为零 用能量平衡法导出对称轴上的点的有限差分方程 如图1 2 17所示 从 点传入到 点的热量Tp2 Tp0Q2 0 2 r 2 z N 1 2 69 r从 点传入到 的热量 72 1 2 4圆柱坐标系的有限差分方程 注 半径 r 2上下点为5 6点 节距 z 点周围温度为T2 或平均温度 设 角方向分成N份图1 2 17对称轴上点Tp5 Tp0Q5 0 r 2 2 N 1 2 70 zTp6 Tp0Q6 0 r 2 2 N 1 2 71 z UTp 10 Tp0 Q Cp r 2 2 z N 1 2 72 整理后得 Tp 10 Tp0Tp2 Tp0Tp5 2Tp0 Tp6 Cp 4 1 2 73 r 2 z 2 73 1 2 4圆柱坐标系的有限差分方程 T1 T T Cp r r r r z z1 r T1 T T r r r rr r r z z1 T 2T T r r r2 z 当r 0时 上式左侧第一项是不定的 所以需要采用变换微分方程的形式 根据罗毕塔法则 T1 T 2T T Cp r r r2 z 1 T 2T r2 2Tlim lim r 0r rr 0 r r r2 T 2T T Cp 2 r2 z 74 1 2 5差分方程的稳定性和收敛性 前面介绍了数值模拟中常用的两种差分格式 显式和隐式 并讨论了不同情况下差分方程的建立 这些差分方程从形式上看都是可以进行计算的 但还有两个问题尚需要研究 那就是 当空间步长 x 时和时