含强间断导热系数的三维扩散方程数值计算

第27卷第9期2015年9月H I G H P OW E R L A S E R A N D P A R T I C L E B E AM S 强激光与粒子束V o l . 27, N o . 9 , S e . 2015 p

*含强间断导热系数的三维扩散方程数值计算

郭少冬, 章明宇, 周海兵, 张树道

() 北京应用物理与计算数学研究所, 北京100094

借鉴各向异性导热的思想, 摘 要: 针对惯性约束聚变领域需要求解的含强间断导热系数的扩散方程,

给出了一种导热系数定义在网格节点的支撑算子改进格式㊂ 改进后的方法保留了支撑算子格式的算子相容性

适合于辐射流体力学问题的模拟㊂

关键词: 惯性约束聚变; 扩散方程; 支撑算子格式; 强间断; 强非线性

:/ 中图分类号: O532 文献标志码: A d o i 10. 11884H P L P B 201527. 092014优点㊂ 数值算例表明, 新方法可以很好地处理含强间断和强非线性的热扩散问题, 计算结果保持较高的精度,

]1-2惯性约束聚变的数值模拟中, 扩散方程的离散求解[是重要的组成部分㊂ 在这一领域 在辐射流体力学㊁

中, 扩散方程的求解有其特殊性和难点:一是由于要与流体力学的拉氏计算相耦合, 因此经常需要在大变形不

2-4]; 规则的网格上求解扩散方程[二是由于等离子体的辐射特性, 扩散方程通常具有强非线性和强间断的特5-6]㊂ 征[

[]-9针对在大变形网格上求解扩散方程的要求, 使用支撑算子方法来推导扩散方程的离散格S h a s h k o v 等7

式㊂ 其主要优点是, 格式能够保证离散后的梯度和散度算子满足连续方程的数学特性, 且算法对网格光滑性的]依赖较小, 即不论网格是否光滑, 算法均可保持二阶精度㊂ 文献[进一步将该算法改进, 使其能很好地处理10强扭曲网格上的求解㊂ 因此在处理大变形扭曲网格方面, 支撑算子格式具有较好的精度和健壮性㊂

在辐射流体力学中, 能量传输的一个特点是电子和辐射热传导率强烈地依赖于温度, 即热传导率在热等离子体中比在冷等离子体中大得多, 此时能量传输通常以陡峭的热波形式出现, 导致扩散方程呈现很强的非线性和间断性, 采用传统的支撑算子格式求解将会导致 数值热障 的出现㊂ 本文将借鉴各向异性扩散的思想, 对支撑算子格式进行改进, 使之能够处理强间断导热系数下的扩散问题㊂

1 支撑算子格式

非定常抛物型的扩散方程可以写为

T Ñ, c T =s x , z ɪ V -㊃ K Ñy , ρt 式中:c 为比热容; T 为温度; K 为导热系数; s 为源项㊂ ρ为密度; () 1

D ㊃ W +ΩT =S ; W =-K ÑT

式中:D 为散度; W 为通量; Ω与具体的时间离散格式有关; S 为离散后的源项㊂

选择散度作为基本算子, 对其进行离散㊂ 根据高斯散度定理, 散度D 的离散表达式可以写为将扩散方程写为一阶系统的形式() 2

1() D ㊃ W ʈ ð (W F ㊃ n 3F )A F V F ) 式(对任意一个三维网格单元都是成立的, 式中:3F 为构成三维网格单元的面; A F 为该面的面积; n F 为该面的

单位外法向量; W F 为该面的通量; V 为整个三维单元的体积㊂

对通量W 进行离散㊂ 在连续意义下, 对于任意的矢量W 和标量T 存在下述恒等式

㊂ 通量与散度算子的离散表达式需要满足积分恒等式(结合散度算子的离散 支撑算子方法的思想是:4)

*收稿日期:2015-05-04; 修订日期:2015-08-10㊃ W d T ÑV +W ㊃ ÑT d V -V =0ʏ ʏ ʏ Ñ㊃ (T W )d V V V () 4; 中国工程物理研究院科学技术发展基金项目(基金项目:国家自然科学基金项目(11205016, 11372050, 11472060) 2014B 0202031, ) 2013A 0101004 ) 作者简介:郭少冬(男, 博士, 副研究员, 从事辐射流体力学数值模拟研究; 1982, u o s d @163. c o m ㊂ g

092014-1

强激光与粒子束

) , ]) , 表达式(经过一系列推导(可参考文献[最终得到通量W 的表达式38W =G -1T

式中

T --T G =ð P n J n 1K -1J V n n P n () 5() 6) 式(的求和是针对单元i 内的所有节点n 进行的㊂ 式中:6P n 为节点n 相邻面的局部编号到全局编号的转换矩

阵; 这里由于导热系数定义在J K 为导热系数, n 为节点n 相邻的3个面的单位外法向矢量组成的雅克比矩阵;

单元中心, 故每个节点体积内的扩散系数均为单元中心值; V n 为节点的角体积㊂

和时间离散格式, 即可给出最终的离散方程, 这里不再赘述㊂ ) , 至此散度和通量的离散已经完成, 接下来只需将其代入扩散表达式(并结合单元界面的通量连续条件2n

数值热障 现象与格式改进2

在辐射流体力学中, 能量传输的特点是介质的导热系数强烈的依赖于自身的温度, 当能量从高温介质传向低温介质时, 高温介质的导热系数要比低温介质的导热系数大得多, 二者的差别可以达到十几个数量级, 这种情况下, 支撑算子格式会产生 数值热障 现象㊂

为了方便说明, 构造如下的非线性扩散问题

∂T Ñc T =㊃ K (T )Ñρ∂t ) i . c . T (r , 0δ(r ) =T 0其中:K (T )=T α㊂ () 7

0㊂ T 0为原点的初始温度㊂ 整个区域的总能量为

V 该问题是一个能量源的非线性扩散问题㊂ 初始时刻, 整个能量集中在坐标原点, 其余空间位置上温度为]d V =[r ) V =c T ρʏ (ρc T )d ʏ ρc T δ(V 00

由于导热系数K 是温度T 的α次方, 导致波前介质和波后介 坐标原点的能量以热波的形式向周围传播,

质的导热系数差别非常大, 存在极强的间断㊂

) 交均匀网格下, 单元界面的通量表达式(退化为5在支撑算子格式中, 导热系数定义在单元中心, 每个网格单元内的导热系数是均匀的, 没有方向性㊂ 在正() 8

() W =K 1(T c L =K 2(T f -T c L 91-T f )/2)/式中:L 为单元的边长; K 1, K 2和T c T c T f 为界面上的温度㊂ 当界面1, 2分别为界面两侧单元的导热系数和温度;

两侧的导热系数差别很大时, 如K 1≫K 2, 且K 2~0时, 等效导热系数为/(2K 1K 2K 1+K 2) ~2K 2~0这就使得热波不能有效的向前传播, 在数值上产生了热障现象㊂

) 上, 如图1中的A 点, 则式(写为5() 10为了克服支撑算子方法处理这类问题的缺陷, 借鉴各向异性导热系数的思想,

将导热系数定义在网格节点

如图1所示, K n 为网格节点A 上定义在笛卡尔坐标系下的

导热系数矩阵, 节点A 的通量W n 可以表示为(ð P J n T -1n n 1-T -K n J n P n V n )W =T () 11

(的面中心通量W f 其A D G F , A B C D , A B E F ) W f W f 1, 2, 3来表示,

中n 为面的单位外法向量() W n =-K n ÑT 12将节点A 的通量W n 用与其相邻的三个面 利用几何关系,

T T T T T () J K n ÑT =-[n W f n W f n W f 13n f 11, f 22, f 33]) 右端的面通量, 引入界面导热系数K f , 即 为了计算式(13[n W f n W f n W f 1, 2, 3]T f 1T f 2T f 3T i a K f K f =-d g (K f 1, 2, 3)J T n ÑT F i . 1 Ti c a l d i s c r e t e e l e m e n t g y p 图1典型离散单元

092014-2 () 14