基于变分原理的二维热传导方程差分格式

© 1994-2010 China Academic Journal Electronic Publishing House. All rights reserved. http://www.cnki.net [文章编号 ] 　10012246X(2002) 0420299206 [收稿日期 ] 2001 - 01 - 03 ; [修回日期 ] 2001 - 04 - 19 [基金项目 ] 国家高技术 863 计划 项目进度计划表范例计划下载计划下载计划下载课程教学计划下载 惯性约束聚变委员会和中国工程物理研究院科学技术基金资助项目 [作者简介 ] 陈光南 (1937 - ) ,男 ,江苏常州 ,研究员 ,博士生导师 ,从事计算 数学 数学高考答题卡模板高考数学答题卡模板三年级数学混合运算测试卷数学作业设计案例新人教版八年级上数学教学计划 方面的研究 ,北京 8009 信箱 12 分箱. 基于变分原理的二维热传导方程差分格式 陈光南 , 　张永慧 (北京应用物理与计算数学研究所 ,北京　100088) [摘 　要 ] 　研究二维热传导方程的差分数值模拟. 用变分原理在不规则结构网格上建立热流通量形式的差分格 式.将热流通量作为未知函数求泛函极值 ,并与温度函数联立求解. 克服通常九点格式用插值方法计算网格边界上 的热传导系数和网格结点上的温度所引入的误差. [关键词 ] 　变分原理 ;二维热传导方程 ;不规则结构网格 [中图分类号 ] 　O24 [文献标识码 ] 　A 0 　引言 数值求解二维辐射流体力学方程组时 ,流体力 学运动使拉格朗日网格发生很大畸变. 因此 ,在不规 则网格上建立二维热传导方程的差分格式就成为至 关重要的问 题 快递公司问题件快递公司问题件货款处理关于圆的周长面积重点题型关于解方程组的题及答案关于南海问题 [1～5 ] . 尤其对许多应用性问题 ,如激光 惯性约束聚变研究中激光与靶耦合过程 ,内爆动力 学压缩过程 ,都会出现热传导方程系数和温度函数 在空间和时间上的强烈非均匀性 ,成为问题特有的 难点. 众所周知 ,数理方程的定解问题可以有不同的 数学提法[1 ] ,如微分方程提法 ,守恒原理提法 ,变分 原理提法. 相应地在网格的基础上进行差分离散化 时也有 3 种不同途径 ,如模拟微分方程及边界条件 时使用格点值的差商来代替微商 ;模拟守恒原理采 取在各网格单元及边界回路上建立离散守恒方程 , 即差分方程的解在某种离散意义下满足守恒律 ;在 网格分区累加的基础上模拟能量表达式 ,从而得到 有限自由度的二次函数极小问题. 这 3 种不同途径 各有其优缺点 :微分途径具有通用、简便等优点 ,但 物理、数学特性保持较差 ,容易导致解失真 ,对复杂 边界条件 ,奇异性、间断性的处理不容易统一等 ;而 变分原理途径是针对根本的物理规律进行模拟 ,容 易保证解的可靠性及保持问题的特性 ,对解、系数、 区域、网格等的不规则性也都比较容易适应. 虽然基 于变分原理所得到的差分格式在形式上可能很复 杂 ,但其构成原则却很简单且不难用程序实现. 在实 际应用中往往可以用不规则网格较少的结点来适应 解的不断变化 ,因此能够以较小的计算机存储量和 较少的计算时间取得较好结果. 在二维不规则网格上求解辐射热传导方程通常 选用微分方程提法 ,用积分插值法建立普通的九点 差分格式进行计算[5 ] . 与五点格式相比 ,尽管九点差 分格式克服了在非正交网格上计算热流通量漏失偏 大 ,以及温度分布变化明显依赖于网格剖分形状等 缺点 ,但是在计算公式中因涉及到网格结点上的温 度值 ,需要计算每个网格四周边界上的热传导系数 , 而这些值要借助插值公式求得[6 ] . 当网格发生强烈 畸变和不同物质之间出现系数间断时 ,这些插值公 式如果处理不当 ,将会导致计算结果失真 ,甚至出现 非物理解. 为了克服上述缺点 ,保证二维网格在大的 畸变情况下数值求解结果的可靠性 ,除了对上述插 值公式选取具有一定明确物理含义和计算精度的公 式外 ,另一种有效途径是本文将要讨论的基于变分 原理的方法. 先将热传导方程写成热流形式 ,使得热 流通量成为独立的未知函数 ,通过变分原理 ,研究热 流泛函的极值并与温度对流方程进行联立求解 ;与 此同时 ,解出每个网格内的温度值以及通过网格四 周边界上的热流通量. 1 　通量流形式下热传导方程变分公式 在单联通域 ( x , y) ∈Ω中讨论热传导方程5 u5 t = div( k grad u) , (1) 第 19 卷 第 4 期 2002 年 7 月 计 　　算 　　物 　　理 　CHINESE JOURNAL OF COMPUTATIONAL PHYSICS 　 Vol. 19 ,No. 4Jul . , 2002 　 © 1994-2010 China Academic Journal Electronic Publishing House. All rights reserved. http://www.cnki.net 其中 u 是温度 , k Ε c > 0 是热传导系数. 在初始时刻 给出温度分布 u ( x , y ,0) = u0 ( x , y) , 　( x , y) ∈Ω , (2) 在域Ω的边界Γ上给出边界条件 ( - k grad u) | Γ = φ, 　t > 0 , (3) 方程 (1)中散度所作用的向量 W = - k grad u 有明 显的物理意义 ,它表示热流通量. 数值求解热传导方 程 (1) ,可以将通量流函数作为独立未知函数求解. 这样方程 (1)又可写成方程组形式5 u5 t = - div W , (4) W = - k grad u , (5) 条件 (3)变为 ( W , n) | Γ = φ, 　t > 0 , (6) 由能量法 ,将边值问题化为等价的求泛函极值问题. 引理 1 　设 A 是正定算子 , u 是方程 Au = f 在 域Ω上的解的充分必要条件是 u 使泛函 J ( u) = ( Au , u) - 2 ( f , u)取极小值. 根据上述引理 ,可将方程 (1)在某一固定时刻看 作为 Au = - div( k grad u) , f = - 5 u5 t ,则有泛函 J ( u) = ( - div( k grad u) , u) - 2 ( - 5 u5 t , u) = ( k grad u ,grad u) + 2 (5 u5 t , u) , (7) 利用表达式 (4) 、(5) ,可将 (7)式改写成泛函 J ( W) = ( - W , - Wk ) + 2 ( - div W , u) =∫Ω | W | 2 k dΩ - 2∫Ωu div WdΩ , (8) 其中 u 和 k 认为是固定值. 这样将求解满足定解条 件 (2) 、(3)的热传导方程 (1) ,转化为在每一时刻求 满足边界条件 (6)的函数类中使泛函 (8)达到极小的 解 W ,再利用方程 (4)和条件 (2)求出未知函数 u . 2 　差分方程组 为了将方程离散化 ,先建立网格. 假设求解区域 Ω是由 4 条曲线围成的. 用两族依赖于参数 (ξ,η) 的曲线 (一般为折线) x = x (ξ,η) , 　 y = y (ξ,η) 将区域Ω剖分为四边形网格. 对于固定的ξ=ξi ( i = 0 ,1 , ⋯, I) , x = x (ξi ,η) , 　 y = y (ξi ,η) 表示一根 y 方向的网格线 , i = 0 和 i = I 分别表示 区域Ω的左、右边界. 同样 ,对于固定的η=ηj ( j = 0 ,1 , ⋯,J ) , x = x (ξ,ηj ) , 　 y = y (ξ,ηj ) 表示一根 x 方向的网格线 , j = 0 和 j = J 分别表示 区域Ω的上、下边界. 设ξ=ξi 和η=ηj 两条网格线 的交点为 Pi , j ( xi , j , yi , j ) , xi , j = x (ξi ,ηj ) , 　 yi , j = y (ξi ,ηj ) . 以 Pi , j , Pi + 1 , j , Pi + 1 , j + 1 , Pi , j + 1 为结点的四边形网格 用Ωi + 1Π2 , j + 1Π2表示 ,并用它表示这个网格的面积. 网 格 Ωi + 1Π2 , j + 1Π2 的 4 条边界依次用 si + 1Π2 , j , si + 1 , j + 1Π2 , si + 1Π2 , j + 1 , si , j + 1Π2表示 ,同样用它们表示相应边界的长 度. 网格 Ωi + 1Π2 , j + 1Π2 有 4 个内角 ,分别用φi + 1Π2 , j + 1Π2i , j , φi + 1Π2 , j + 1Π2i + 1 , j ,φi + 1Π2 , j + 1Π2i + 1 , j + 1 ,φi + 1Π2 , j + 1Π2i , j + 1 表示 ,这里下标对应 网格结点的编号 ,上标对应网格的编号 (图 1) . 图 1 　网格示意图 Fig11 　Grid sketch 建立离散化方程时 ,未知函数 u 取在网格中心 u n i + 1Π2 , j + 1Π2 ,其中上标 n 表示 t = tn = nτ时刻的量. 为 简化起见 ,下面凡是 n 时刻的量 ,都略去上标 n . 通 量流函数 W 取在网格边上. 在网格四周作法线单位 向量 , 其方向指向 i 和 j 增加的方向 : eη, i + 1Π2 , j , eξ, i + 1 , j + 1Π2 , eη, i + 1Π2 , j + 1 , eξ, i , j + 1Π2 . 通量流函数 W 在这 些法线向量上的投影记做 wη, i + 1Π2 , j , wξ, i + 1 , j + 1Π2 , wη, i + 1Π2 , j + 1 , Wξ, i , j + 1Π2 . 离散泛函 (8) 式时 ,需求出| W| 2i + 1Π2 , j + 1Π2的逼近 值. 借助上面定义的 4 条网格边界的法线方向的单 位向量构成 4 个局部坐标系 : ( eη, i + 1Π2 , j , eξ, i + 1 , j + 1Π2 ) , ( eξ, i + 1 , j + 1Π2 , eη, i + 1Π2 , j + 1 ) , ( eη, i + 1Π2 , j + 1 , eξ, i , j + 1Π2 ) , ( eξ, i , j + 1Π2 , eη, i + 1Π2 , j ) . 这样 ,利用向量 W 在任意坐标 系 ( eξ , eη)上的投影 ( wξ , wη)表示它的模 ,则有 | W | 2 = ( wξ) 2 + ( wη) 2 - 2 ( wξ) ( wη) cosψ sin2ψ , 　(9) 003 计 　　算 　　物 　　理 第 19 卷 　 © 1994-2010 China Academic Journal Electronic Publishing House. All rights reserved. http://www.cnki.net 其中ψ为基向量 eξ 和 eη 之间的夹角. 为求出| W| 2i + 1Π2 , j + 1Π2的值 ,需要先求出以上 4 个坐标系中 W 的模的表达式 ,然后再作算术平均 ,同时考 虑到夹角ψ和φ之间的关系 ,可得 4 | W | 2i +1Π2 , j+1Π2 = ( wξ, i , j +1Π2 ) 2 + ( wη, i +1Π2 , j ) 2 + 2 wξ, i , j+1Π2 wη, i +1Π2 , j cosφi +1Π2 , j +1Π2i , j(sinφi +1Π2 , j +1Π2i , j ) 2 + ( wη, i +1Π2 , j ) 2 + ( wξ, i +1 , j +1Π2 ) 2 - 2 wη, i +1Π2 , jwξ, i +1 , j+1Π2 cosφi +1Π2 , j+1Π2i +1 , j (sinφi +1Π2 , j +1Π2i +1 , j ) 2 + ( wξ, i +1 , j +1Π2 ) 2 + ( wη, i +1Π2 , j+1 ) 2 + 2 wξ, i +1 , j +1Π2 wη, i +1Π2 , j +1 cosφi +1Π2 , j+1Π2i +1 , j+1 (sinφi +1Π2 , j+1Π2i +1 , j+1 ) 2 + ( wη, i +1Π2 , j+1 ) 2 + ( wξ, i , j +1Π2 ) 2 - 2 wη, i +1Π2 , j+1 wξ, i , j+1Π2 cosφi +1Π2 , j+1Π2i , j +1 (sinφi +1Π2 , j+1Π2i , j +1 ) 2 , i = 0 ,1 , ⋯, I - 1 , 　 j = 0 ,1 , ⋯,J - 1. (10) 泛函 (8)中的另外一项在网格Ωi + 1Π2 , j + 1Π2上可近似为 ∫Ωi +1Π2 , j+1Π2 udiv WdΩ = ui +1Π2 , j+1Π2∮Γi+1Π2 , j+1Π2 W ·nd s = ui +1Π2 , j+1Π2 ( wξ, i +1 , j+1Π2 si +1 , j+1Π2 - wξ, i , j +1Π2 si , j+1Π2 + wη, i +1Π2 , j+1 si +1Π2 , j+1 - wη, i +1Π2 , jsi +1Π2 , j ) , 于是泛函 (8)离散化后 ,可写成 J h ( W) = ∑ J - 1 j = 0 ∑ I - 1 i = 0 | Wi +1Π2 , j+1Π2 | 2 ki +1Π2 , j+1Π2 Ωi +1Π2 , j+1Π2 - 2 ∑J - 1j = 0 ∑I - 1i = 0 ui +1Π2 , j+1Π2 ( wξ, i +1 , j+1Π2 si +1 , j +1Π2 - wξ, i , j +1Π2 si , j+1Π2 + wη, i +1Π2 , j+1 si +1Π2 , j+1 - wη, i +1Π2 , jsi +1Π2 , j ) , (11a) 其中| Wi + 1Π2 , j + 1Π2 | 2 用 (10)式代入. 这样 ,就化为求泛函 (8)式多元函数的极小值问题. 表达式 (11) 还可以写成 另一种近似计算公式 J h ( W) = ∑ J - 1 j = 0 ∑ I - 1 i = 0 1 4 { Ωi +1Π2 , j +1Π2i , j ki , j ( wξ, i , j +1Π2 ) 2 + ( wη, i +1Π2 , j ) 2 + 2 wξ, i , j+1Π2 wη, i +1Π2 , j cosφi +1Π2 , j +1Π2i , j (sinφi +1Π2 , j+1Π2i , j ) 2 + Ωi +1Π2 , j +1Π2i +1 , j ki +1 , j ( wη, i +1Π2 , j ) 2 + ( wξ, i +1 , j +1Π2 ) 2 - 2 wη, i +1Π2 , jwξ, i +1 , j+1Π2 cosφi +1Π2 , j+1Π2i +1 , j (sinφi +1Π2 , j +1Π2i +1 , j ) 2 + Ωi +1Π2 , j +1Π2i +1 , j+1 ki +1 , j+1 ( wξ, i +1 , j +1Π2 ) 2 + ( wη, i +1Π2 , j+1 ) 2 + 2 wξ, i +1 , j +1Π2 wη, i +1Π2 , j +1 cosφi +1Π2 , j+1Π2i +1 , j+1 (sinφi +1Π2 , j+1Π2i +1 , j+1 ) 2 + (11b) Ωi +1Π2 , j +1Π2i , j+1 ki , j+1 ( wη, i +1Π2 , j+1 ) 2 + ( wξ, i , j +1Π2 ) 2 - 2 wη, i +1Π2 , j+1 wξ, i , j+1Π2 cosφi +1Π2 , j +1Π2i , j +1 (sinφi +1Π2 , j+1Π2i , j +1 ) 2 - 2 ∑ J - 1 j = 0 ∑ I - 1 i = 0 ui +1Π2 , j+1Π2 ( wξ, i +1 , j+1Π2 si +1 , j+1Π2 - wξ, i , j+1Π2 si , j+1Π2 + wη, i +1Π2 , j+1 si +1Π2 , j +1 - wη, i +1Π2 , jsi +1Π2 , j ) , 其中Ωi + 1Π2 , j + 1Π2i , j ,Ωi + 1Π2 , j + 1Π2i + 1 , j ,Ωi + 1Π2 , j + 1Π2i + 1 , j + 1 ,Ωi + 1Π2 , j + 1Π2i , j + 1 分别是以结点 ( i , j) , ( i + 1 , j) ( i + 1 , j + 1) , ( i , j + 1) 为起 点 ,围绕 ( i + 1Π2 , j + 1Π2)网格求出面积的四分之一. ki , j , ki + 1 , j , ki + 1 , j + 1 , ki , j + 1分别是以 ( i , j) , ( i + 1 , j) , ( i + 1 , j + 1) , ( i , j + 1)为中心的相邻 4 个网格内 k 值的某种平均. 例如 kij = 1 4 ( ki +1Π2 , j +1Π2 + ki - 1Π2 , j+1Π2 + ki - 1Π2 , j - 1Π2 + ki +1Π2 , j - 1Π2 ) . 　　离散温度对流方程 (4) ,可将它在网格区域Ωi + 1Π2 , j + 1Π2上积分 ,得到差分格式 1 τ( ui +1Π2 , j+1Π2 - un - 1i +1Π2 , j+1Π2 )Ωi +1Π2 , j+1Π2 = - ( wξ, i +1 , j +1Π2 si +1 , j +1Π2 - wξ, i , j+1Π2 si , j +1Π2 + wη, i +1Π2 , j+1 si +1Π2 , j+1 - wη, i +1Π2 , jsi +1Π2 , j ) , i = 0 ,1 , ⋯, I - 1 , 　 j = 0 ,1 , ⋯,J - 1. 　　　　　(12) 这样 ,求解热传导方程 (1)的过程就成为先求出使泛函 (11)达到极小值的量 wξ, i , j + 1Π2 ( i = 1 ,2 , ⋯, I - 1 ; j = 0 , 1 , ⋯,J - 1)和 wη, i + 1Π2 , j ( i = 0 ,1 , ⋯, I - 1 ; j = 1 ,2 , ⋯, J - 1) ,然后从 (12) 式求出 ui + 1Π2 , j + 1Π2 ( i = 0 ,1 , ⋯, I - 1 ; 103　第 4 期 陈光南等 :基于变分原理的二维热传导方程差分格式 © 1994-2010 China Academic Journal Electronic Publishing House. All rights reserved. http://www.cnki.net j = 0 ,1 , ⋯,J - 1) . 使泛函 J h ( W)达到极值的 wξ, i , j + 1Π2和 wη, i + 1Π2 , j可以从 J h ( W) 对它们求导数并等于零的方程组中推出. J h ( W)对 wξ, i , j + 1Π2求导数并等于零的方程为 Ωi +1Π2 , j+1Π2 2 ki +1Π2 , j +1Π2 wξ, i , j+1Π2 + wη, i +1Π2 , j cosφi +1Π2 , j +1Π2i , j(sinφi +1Π2 , j +1Π2i , j ) 2 + wξ, i , j+1Π2 - wη, i +1Π2 , j +1 cosφi +1Π2 , j+1Π2i , j +1(sinφi +1Π2 , j +1Π2i , j +1 ) 2 + Ωi - 1Π2 , j+1Π2 2 ki - 1Π2 , j +1Π2 wξ, i , j+1Π2 + wη, i - 1Π2 , j+1 cosφi - 1Π2 , j+1Π2i , j +1(sinφi - 1Π2 , j +1Π2i , j+1 ) 2 + wξ, i , j +1Π2 - wη, i - 1Π2 , j cosφi - 1Π2 , j+1Π2i , j(sinφi - 1Π2 , j +1Π2i , j ) 2 + 2 ( ui +1Π2 , j+1Π2 - ui - 1Π2 , j+1Π2 ) si , j+1Π2 = 0 , 　　i = 1 ,2 , ⋯, I - 1 , 　j = 0 ,1 , ⋯,J - 1. (13) 同样 ,J h ( W)对 wη, i + 1Π2 , j求导数并等于零的方程为 Ωi +1Π2 , j+1Π2 2 ki +1Π2 , j +1Π2 wη, i +1Π2 , j + wξ, i , j +1Π2 cosφi +1Π2 , j +1Π2i , j(sinφi +1Π2 , j+1Π2i , j ) 2 + wη, i +1Π2 , j - wξ, i +1 , j+1Π2 cosφi +1Π2 , j +1Π2i +1 , j(sinφi +1Π2 , j +1Π2i +1 , j ) 2 + Ωi +1Π2 , j - 1Π2 2 ki +1Π2 , j - 1Π2 wη, i +1Π2 , j + wξ, i +1 , j - 1Π2 cosφi +1Π2 , j - 1Π2i +1 , j(sinφi +1Π2 , j - 1Π2i +1 , j ) 2 + wη, i +1Π2 , j - wξ, i , j - 1Π2 cosφi +1Π2 , j - 1Π2i , j(sinφi +1Π2 , j - 1Π2i , j ) 2 + 2 ( ui +1Π2 , j+1Π2 - ui +1Π2 , j - 1Π2 ) si +1Π2 , j = 0 , j = 1 ,2 , ⋯,J - 1 ; i = 0 ,1 , ⋯, I - 1. 　　　　(14) 方程组 (13) 、(14)中的 ui + 1Π2 , j + 1Π2 , ui - 1Π2 , j + 1Π2 , ui + 1Π2 , j - 1Π2等值用公式 (12)代入 ,由此可得到下列计算公式 - τsi +1 , j+1Π2 si , j +1Π2Ωi +1Π2 , j+1Π2 wξ, i +1 , j+1Π2 + Ωi +1Π2 , j+1Π24 ki +1Π2 , j+1Π2 1sin2φi +1Π2 , j+1Π2i , j + 1sin2φi +1Π2 , j+1Π2i , j+1 + Ωi - 1Π2 , j+1Π24 ki - 1Π2 , j +1Π2 1sin2φi - 1Π2 , j +1Π2i , j + 1 sin2φi - 1Π2 , j +1Π2i , j+1 +τs2i , j+1Π2 1Ωi +1Π2 , j+1Π2 + 1Ωi - 1Π2 , j+1Π2 wξ, i , j+1Π2 - τsi , j+1Π2 si - 1 , j+1Π2Ωi - 1Π2 , j +1Π2 wξ, i - 1 , j +1Π2 = Ωi +1Π2 , j+1Π2 4 ki +1Π2 , j +1Π2 cosφi +1Π2 , j+1Π2i , j+1sin2φi +1Π2 , j +1Π2i , j+1 + τsi , j+1Π2 si +1Π2 , j+1Ωi +1Π2 , j+1Π2 wη, i +1Π2 , j +1 - Ωi - 1Π2 , j+1Π24 ki - 1Π2 , j+1Π2 cosφi - 1Π2 , j +1Π2i , j +1sin2φi - 1Π2 , j+1Π2i , j+1 + τsi , j+1Π2 si - 1Π2 , j+1 Ωi - 1Π2 , j+1Π2 wη, i - 1Π2 , j+1 - Ωi +1Π2 , j+1Π24 ki +1Π2 , j+1Π2 cosφi +1Π2 , j +1Π2i , jsin2φi +1Π2 , j+1Π2i , j + τsi , j +1Π2 si +1Π2 , jΩi +1Π2 , j +1Π2 wη, i +1Π2 , j + Ωi - 1Π2 , j+1Π2 4 ki - 1Π2 , j +1Π2 cosφi - 1Π2 , j+1Π2i , jsin2φi - 1Π2 , j +1Π2i , j + τsi , j+1Π2 si - 1Π2 , jΩi - 1Π2 , j+1Π2 wη, i - 1Π2 , j - si , j +1Π2 ( un - 1i +1Π2 , j+1Π2 - un - 1i - 1Π2 , j+1Π2 ) , i = 1 ,2 , ⋯, I - 1 , 　j = 0 ,1 , ⋯,J - 1. 　　　(15) - τsi +1Π2 , j+1 si +1Π2 , jΩi +1Π2 , j +1Π2 wη, i +1Π2 , j +1 + Ωi +1Π2 , j +1Π24 ki +1Π2 , j+1Π2 1sin2φi +1Π2 , j+1Π2i , j + 1sin2φi +1Π2 , j+1Π2i +1 , j + Ωi +1Π2 , j - 1Π24 ki +1Π2 , j - 1Π2 1sin2φi +1Π2 , j - 1Π2i , j + 1 sin2φi +1Π2 , j - 1Π2i +1 , j +τs2i +1Π2 , j 1Ωi +1Π2 , j+1Π2 + 1Ωi +1Π2 , j - 1Π2 wη, i +1Π2 , j - τsi +1Π2 , jsi +1Π2 , j - 1Ωi +1Π2 , j - 1Π2 wη, i +1Π2 , j - 1 = Ωi +1Π2 , j+1Π2 4 ki +1Π2 , j+1Π2 cosφi +1Π2 , j +1Π2i +1 , jsin2φi +1Π2 , j+1Π2i +1 , j + τsi +1Π2 , jsi +1 , j +1Π2Ωi +1Π2 , j +1Π2 wξ, i +1 , j+1Π2 - Ωi +1Π2 , j - 1Π24 ki +1Π2 , j - 1Π2 cosφi +1Π2 , j - 1Π2i +1 , jsin2φi +1Π2 , j - 1Π2i +1 , j + τsi +1Π2 , jsi +1 , j - 1Π2 Ωi +1Π2 , j - 1Π2 wξ, i +1 , j - 1Π2 - Ωi +1Π2 , j+1Π24 ki +1Π2 , j+1Π2 cosφi +1Π2 , j +1Π2i , jsin2φi +1Π2 , j+1Π2i , j + τsi +1Π2 , jsi , j+1Π2Ωi +1Π2 , j +1Π2 wξ, i , j+1Π2 + (Ωi +1Π2 , j - 1Π24 ki +1Π2 , j - 1Π2 cosφi +1Π2 , j - 1Π2i , jsin2φi +1Π2 , j - 1Π2i , j + τsi +1Π2 , jsi , j - 1Π2Ωi +1Π2 , j - 1Π2 ) wξ, i , j - 1Π2 - si +1Π2 , j ( un - 1i +1Π2 , j +1Π2 - un - 1i +1Π2 , j - 1Π2 ) , j = 1 ,2 , ⋯,J - 1 , 　i = 0 ,1 , ⋯, I - 1. 　　　(16) 这是一组含有 ( I + 1) J + I ( J + 1)个未知数 wξ, i , j + 1Π2 ( i = 0 ,1 , ⋯, I ; j = 0 ,1 , ⋯, J - 1) 和 wη, i + 1Π2 , j ( i = 0 ,1 , ⋯, I - 1 ; j = 0 ,1 , ⋯,J )的方程组. 方程组 (15) 、(16)共有 ( I - 1) J + I ( J - 1) 个方程 ,还缺 2 ( I + J ) 个关系式. 这 正好由边界条件给出. 当边界满足绝热条件时 ,则有 wξ,0 , j +1Π2 = 0 , 　j = 0 ,1 , ⋯,J - 1 ; 　　　wξ, I , j+1Π2 = 0 , 　j = 0 ,1 , ⋯,J - 1 ; 　 (17) wη, i +1Π2 ,0 = 0 , 　i = 0 ,1 , ⋯, I - 1 ; 　　　wη, i +1Π2 , J = 0 , 　i = 0 ,1 , ⋯, I - 1. (18) 这样方程组 (15) ～ (18) 构成一个封闭的方程组 ,就可以用追赶迭代法求解 ,也可以用其它迭代方法进行 203 计 　　算 　　物 　　理 第 19 卷 　 © 1994-2010 China Academic Journal Electronic Publishing House. All rights reserved. http://www.cnki.net 求解. 3 　数值试验 为了检验本文提出的基于变分原理的二维热传 导方程差分格式的正确性和解题适应能力 ,在不均 匀、非正交网格剖分下 ,使用二维程序计算出一维对 称物理条件下的物理图像. 图 2 所示网格剖分为角度方向θ∈(0 ,πΠ4) 等 分为 12 片 ,θ∈(πΠ4 ,πΠ2)也等分为 12 片 ,二者剖分 方式一样 :径向共分 40 个点 ,在θ= 0°网格线上为均 匀等分 ,其步长为δr(0)1 =δr(0)2 = 10Π40 = 0125 ;在θ= 45°网格线上为分段均匀等分 ,步长分别为δr(πΠ4)1 = 715Π20 = 01375 ,δr(πΠ4)2 = (10 - 715)Π20 = 01125. 其它 角度网格线上的步长以等差级数变化剖分 ,即 δr( i)1 = δr( i - 1)1 + (δr(πΠ4)1 - δr(0)1 ) , δr( i)2 = δr( i - 1)2 + (δr(πΠ4)2 - δr(0)2 , i = (1 ,2 , ⋯,11) π48. 在外边界上各片加入相同的温度外源 T (θ, t ) = w1 exp[ w2 ( t - t0 ) 2 ] . 计算结果表明 ,尽管网格剖分 既不均匀又不正交 ,但温度分布基本上仍然是均匀 地向里传播 (图 3) . 这说明差分格式具有良好的对 称性和适应能力 ,它不受网格变形的强烈影响. 图 2 　二维网格的剖分 Fig. 2 　Grid on 2D space domain 另外 ,我们将本文提供的差分格式应用于激光 聚变研究中. 考察激光加源的不均匀性对激光直接 驱动玻璃壳充 DT 气靶和内爆压缩的二维影响. 根 据日本大阪大学 ILE设计的 # 3826 靶在 GEKKO XII 激光装置上的实验 ,我们进行了二维程序计算 ,并与 他们的 STA - 2D 程序计算作比较[7 ,8 ] . 图 3 　二维温度分布 Fig. 3 　Temperature distribution 　　3826 # 模型的基本参数 :靶丸直径 D = 1235μm , 玻璃壳厚度ΔR = 1131μm ,壳内 DT 气体的压强为 PDT = 612 atm. 激光功率为 13 kJ·ns - 1 ,高斯源 ,波长 为λL = 0153μm. 激光源的激光吸收功率 W (θ, t) = FPRT ( t) FPRF (θ)表示在 t 时刻靶球表面θ处单位 立体角的激光功率. FPRT ( t)取高斯源 , FPRF (θ) 为 空间分布因子. 表示式分别为 FPRT ( t) = 81915 ×105 exp [ - 27713( t - 0115) 2 ] , FPRF (θ) = (1 + A0 cos( nθ) ) 1 + 0125A0 ( - 1) n - 1 - 1 n - 1 - ( - 1) n+1 - 1 n + 1 , 其中 A0 是激光扰动振幅 , n 是谐波数 ,取 n = 6. 使用本文提出的热传导方程差分格式 ,编制了 LARED - D 激光内爆二维程序 ,并进行了数值模拟. 图 4 给出了 t = 119 ns 时刻网格界面的位置 ;图 5 给 出了玻璃壳内界面运动轨迹. 同时在表 1 中列出了 LARED - D 程序和 ILE 的 STA - 2D 程序的计算结 果. 对比表明一维与二维中子产额比较接近 ,符合相 同的物理规律. 图 4 　A0 = 015 , t = 119ns 时刻流线 Fig14 　Mesh position at t = 119ns ,A0 = 015 303　第 4 期 陈光南等 :基于变分原理的二维热传导方程差分格式 © 1994-2010 China Academic Journal Electronic Publishing House. All rights reserved. http://www.cnki.net 图 5 　A0 = 015 玻璃壳内界面运动轨迹 Fig15 　Moving locus of interface between SiO2 and DT 表 1 　LARED - D 程序和 STA - 2D 程序计算 中子产额结果对比 Table 1 　 Numerical results of neutron yield comparison between LARED - D and STA - 2D code Y1Dn Y2Dn Yexn 分片数 径向分点 STA - 2D 2 ×1013 1 ×1013 1 ×1013 24 20 + 80 LARED - D 3. 4 ×1013 116 ×1013 48 20 + 80 [参 　考 　文 　献 ] [ 1 ] 　冯康. 基于变分原理的差分格式 [J ] . 应用数学与计 算数学 ,1965 ,2 (4) . [ 2 ] 　Самарский А А. Численные методырешеиия многомерныхэадачмеханикиифиэики[J ] . ЖВМи МФТ,1980 ,20 (6) :1416 - 1464. [ 3 ] 　КорщияТК,ТищкинВФ, ФоворскийАП, Щашков М Ю. Вариационныйподходкпостр2оениюраэностных схем дляуравнения теп2лопроводности накри волинейныхсетках[J ] . ЖВМиМФТ,1980 ,20 (2) : 401 - 421. [ 4 ] 　ЛебоИГ. ПоповИВ, РоэановВБ, Ти2шкинВФ. Двумерноечисленноемоделированиенагреваисжатия лаэерныхмишеней[J ] . ИПМАНСССР,препринт2 , 1993. [ 5 ] 　李德元 ,陈光南. 抛物型方程差分方法引论 [M] . 北 京 :科学出版社 ,1998. [ 6 ] 　符尚武 ,付汉清 ,沈隆钧 ,等. 二维三温能量方程的九 点差分格式及其迭代解法 [J ] . 计算物理 , 1998 ,15 (4) :489 - 497. [ 7 ] 　Takabe H ,Yamanaka M ,Mima K,et al. Phys Fluids ,1988 , 31 :2284. [ 8 ] 　Takabe H.Lecture Note on“Theoretical base and numerical method for laser fusion research”[ R ] . Lectured at IAPCM Beijing ,1997 ,Jul. 28 - Aug. 2. DIFFERENCE SCHEME BASED ON VARIATIONAL PRINCIPLE FOR TWO DIMENSIONAL THERMAL CONDUCTION EQUATION CHEN Guang2nan , 　ZHANG Yong2hui ( Institute of Applied Physics and Computational Mathematics , Beijing 　100088 , China) [ Abstract ] 　Difference numerical simulation for two dimensional thermal conduction equation is studied. Using the variational principle ,the difference scheme with thermal2flux form in the irregular structured mesh is presented. The thermal flux as an unknown function is obtained by minimizing the functional and solved with the equation of temperature function simultaneously. This scheme overcomes the defects of general nine2point schemes ,which have to calculate the thermal conduction coefficients on the cell boundaries and the temperatures located in the cell vertices by interpolation. [ Key words] 　variational principle ;two dimensional thermal conduction equation ; irregular structured mesh [ Received date ] 2001 - 01 - 03 ; [ Revised date ] 2001 - 04 - 19 403 计 　　算 　　物 　　理 第 19 卷