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

第27卷第9期强激光与粒子束V o l .27,N o .92015年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 SS e p.,2015 含强间断导热系数的三维扩散方程数值计算*郭少冬, 章明宇, 周海兵, 张树道(北京应用物理与计算数学研究所,北京100094) 摘 要: 针对惯性约束聚变领域需要求解的含强间断导热系数的扩散方程,借鉴各向异性导热的思想,给出了一种导热系数定义在网格节点的支撑算子改进格式㊂改进后的方法保留了支撑算子格式的算子相容性优点㊂数值算例表明,新方法可以很好地处理含强间断和强非线性的热扩散问题,计算结果保持较高的精度,适合于辐射流体力学问题的模拟㊂关键词: 惯性约束聚变; 扩散方程; 支撑算子格式; 强间断; 强非线性中图分类号: O 532 文献标志码: A d o i :10.11884/H P L P B 201527.092014在辐射流体力学㊁惯性约束聚变的数值模拟中,扩散方程的离散求解[1-2]是重要的组成部分㊂在这一领域中,扩散方程的求解有其特殊性和难点:一是由于要与流体力学的拉氏计算相耦合,因此经常需要在大变形不规则的网格上求解扩散方程[2-4];二是由于等离子体的辐射特性,扩散方程通常具有强非线性和强间断的特征[5-6]㊂针对在大变形网格上求解扩散方程的要求,S h a s h k o v 等[7-9]使用支撑算子方法来推导扩散方程的离散格式㊂其主要优点是,格式能够保证离散后的梯度和散度算子满足连续方程的数学特性,且算法对网格光滑性的依赖较小,即不论网格是否光滑,算法均可保持二阶精度㊂文献[10]进一步将该算法改进,使其能很好地处理强扭曲网格上的求解㊂因此在处理大变形扭曲网格方面,支撑算子格式具有较好的精度和健壮性㊂在辐射流体力学中,能量传输的一个特点是电子和辐射热传导率强烈地依赖于温度,即热传导率在热等离子体中比在冷等离子体中大得多,此时能量传输通常以陡峭的热波形式出现,导致扩散方程呈现很强的非线性和间断性,采用传统的支撑算子格式求解将会导致 数值热障 的出现㊂本文将借鉴各向异性扩散的思想,对支撑算子格式进行改进,使之能够处理强间断导热系数下的扩散问题㊂1 支撑算子格式非定常抛物型的扩散方程可以写为ρc ∂T ∂t -Ñ㊃K ÑT =s , x ,y ,z ɪV (1)式中:ρ为密度;c 为比热容;T 为温度;K 为导热系数;s 为源项㊂将扩散方程写为一阶系统的形式D ㊃W +ΩT =S ; W =-K ÑT (2)式中:D 为散度;W 为通量;Ω与具体的时间离散格式有关;S 为离散后的源项㊂选择散度作为基本算子,对其进行离散㊂根据高斯散度定理,散度D 的离散表达式可以写为D ㊃W ʈ1V ðFW F ㊃n ()F A F(3)式(3)对任意一个三维网格单元都是成立的,式中:F 为构成三维网格单元的面;A F 为该面的面积;n F 为该面的单位外法向量;W F 为该面的通量;V 为整个三维单元的体积㊂对通量W 进行离散㊂在连续意义下,对于任意的矢量W 和标量T 存在下述恒等式ʏVT Ñ㊃W d V +ʏVW ㊃ÑT d V -ʏVÑ㊃T ()W dV =0(4) 支撑算子方法的思想是:通量与散度算子的离散表达式需要满足积分恒等式(4)㊂结合散度算子的离散*收稿日期:2015-05-04; 修订日期:2015-08-10基金项目:国家自然科学基金项目(11205016,11372050,11472060);中国工程物理研究院科学技术发展基金项目(2014B 0202031,表达式(3),经过一系列推导(可参考文献[8]),最终得到通量W 的表达式W =G -1T(5)式中G =ðnP T n J -1n K -1J -Tn P nV n (6)式(6)的求和是针对单元i 内的所有节点n 进行的㊂式中:P n 为节点n 相邻面的局部编号到全局编号的转换矩阵;J n 为节点n 相邻的3个面的单位外法向矢量组成的雅克比矩阵;K 为导热系数,这里由于导热系数定义在单元中心,故每个节点体积内的扩散系数均为单元中心值;V n 为节点的角体积㊂至此散度和通量的离散已经完成,接下来只需将其代入扩散表达式(2),并结合单元界面的通量连续条件和时间离散格式,即可给出最终的离散方程,这里不再赘述㊂2数值热障 现象与格式改进在辐射流体力学中,能量传输的特点是介质的导热系数强烈的依赖于自身的温度,当能量从高温介质传向低温介质时,高温介质的导热系数要比低温介质的导热系数大得多,二者的差别可以达到十几个数量级,这种情况下,支撑算子格式会产生 数值热障 现象㊂为了方便说明,构造如下的非线性扩散问题ρc ∂T ∂t=Ñ㊃()K T ÑT i .c . T (r ,0)=T 0δ(r )(7)其中:()K T =T α㊂该问题是一个能量源的非线性扩散问题㊂初始时刻,整个能量集中在坐标原点,其余空间位置上温度为0㊂T 0为原点的初始温度㊂整个区域的总能量为ʏVρ()c T d V =ʏVρc T 0δ(r [])d V =ρc T 0(8) 坐标原点的能量以热波的形式向周围传播,由于导热系数K 是温度T 的α次方,导致波前介质和波后介质的导热系数差别非常大,存在极强的间断㊂在支撑算子格式中,导热系数定义在单元中心,每个网格单元内的导热系数是均匀的,没有方向性㊂在正交均匀网格下,单元界面的通量表达式(5)退化为W =K 1T c 1-T ()f /L =K 2T f -T ()c 2/L (9)式中:L 为单元的边长;K 1,K 2和T c 1,T c 2分别为界面两侧单元的导热系数和温度;T f 为界面上的温度㊂当界面两侧的导热系数差别很大时,如K 1≫K 2,且K 2~0时,等效导热系数为2K 1K 2/(K 1+K 2)~2K 2~0(10)这就使得热波不能有效的向前传播,在数值上产生了热障现象㊂为了克服支撑算子方法处理这类问题的缺陷,借鉴各向异性导热系数的思想,将导热系数定义在网格节点F i g .1 T y pi c a l d i s c r e t e e l e m e n t 图1典型离散单元上,如图1中的A 点,则式(5)写为ðnP T n J-1nK -1n J -T n P n V ()n W =T(11) 如图1所示,K n 为网格节点A 上定义在笛卡尔坐标系下的导热系数矩阵,节点A 的通量W n 可以表示为W n =-K n ÑT (12) 利用几何关系,将节点A 的通量W n 用与其相邻的三个面(A D G F ,A B C D ,A B E F )的面中心通量W f 1,W f 2,W f 3来表示,其中n 为面的单位外法向量J T n K n ÑT =-n T f 1W f 1,n T f 2W f 2,n T f 3W []f 3T(13) 为了计算式(13)右端的面通量,引入界面导热系数K f ,即n Tf 1W f 1,n Tf 2W f 2,n Tf 3W []f 3T=-d i a g K f 1,K f 2,K ()f 3Jn TÑT强激光与粒子束式中:d i a g 表示以括号中的变量为主对角元素而形成的对角矩阵㊂根据式(13)和(14),得到节点A 的导热系数K n =J n -Td i a g K f 1,K f 2,K ()f 3J nT (15)式中:面导热系数K f 利用该面中心温度计算得到㊂在传统局部支撑算子格式中,导热系数定义在单元中心,面通量的推导只是在本单元内满足离散算子的积分恒等式,没有相邻单元的信息,导致该单元的面通量仅由本单元的物理量和几何量决定㊂当本单元的导热系数与相邻单元的导热系数差别很大时,即会产生如式(10)所示的数值热障㊂改进后的处理方法,将导热系数定义在单元节点和单元面上,构造相应的节点通量和面通量,根据节点通量和面通量的关系,得到节点导热系数与面导热系数的关系㊂由于面导热系数取决于面中心温度,而面中心温度由相邻两个单元的温度线性插值得到,故最终节点导热系数将依赖于相邻两个单元的温度,如式(15)㊂此时将会避免出现式(10)的数值热障㊂当网格为均匀正交六面体时,式(14)中的雅可比矩阵J n 变为单位矩阵,节点导热系数K n 变为对角矩阵K n =d i a g K T ()f 1,K T ()f 2,K T ()()f 3(16)将其代入式(5),注意到六面体单元有8个节点,所有节点求和,结合雅可比矩阵J n 和转换矩阵P n 的特点,得到通量W =d i a g K T ()f 1T c -T ()f 1L 1,K T ()f 2T c -T ()f 2L 2, ,K T ()f 6T c -T ()f 6L æèçöø÷6(17)式中:T c 为六面体单元中心的温度;T f 1,T f 2, ,T f 6分别为单元六个面上的温度㊂每个面的通量与该单元面中心的导热系数有关㊂现在考察相邻的两个单元㊂设T c 1和T c 2分别是这两个单元的中心温度,根据式(17),这两个单元的公共面上的通量,可以写为W =K f 1T c 1-T f L =K f 2T f -T c 2L(18)式中:T f 为公共面的面中心温度;K f 1和K f 2是面导热系数㊂对于均匀正交六面体单元K f 1=K T ()f 1=KT c 1+T c 2æèçöø÷2=K f 2(19)因此,两个单元的界面(公共面)的等效导热系数为2K f 1K f 2K f 1+K f 2=K T c 1+T c 2æèçöø÷2(20)此时该等效导热系数不会趋向于0,因而不会产生数值热障㊂对于任意的三维网格,面通量矩阵虽然不完全等于式(17)中的对角矩阵,但每个面的通量仍与该面导热系数有关㊂由于面导热系数取决于面中心温度,面中心温度由相邻两个单元中心的温度线性插值得到,使得最终通量矩阵的计算依赖于相邻单元中心的温度插值信息,因而也不会产生数值热障㊂需要进一步说明的是,支撑算子最重要的特性是通量与散度的离散表达式能够满足相容性:即积分恒等式(4)㊂依据该恒等式,离散的通量必须具有式(6)的形式㊂本节所做的改进是将导热系数定义至单元节点上,从而可以引入周围单元的贡献㊂根据新的导热系数推导出的通量离散表达式仍然满足式(6)㊂因此新的格式并没有破坏算子间的相容关系式㊂此外,与传统支撑算子格式一样,新格式也具有局部的模板点,离散得到的系数矩阵仍然是稀疏的对称正定矩阵,可以采用许多高效的迭代方法来求解㊂3 数值算例3.1 含强间断导热系数的两介质扩散问题采用上述方法对含有强间断导热系数的两介质扩散问题进行求解,相应的方程和边界条件为∂T ∂t=Ñ㊃K ÑT +z 2郭少冬等:含强间断导热系数的三维扩散方程数值计算b .c .∂T ∂x =0 x =0,x =1;∂T ∂y=0 y =0,y =1;∂T ∂z=0 z =0;T +2K ∂T ∂z =0 z =ìîíïïïïïïïïïï1(21) 在本问题中含有两种导热系数相差很大的介质,且两种介质的导热系数相差6个数量级㊂如式(22)所示,在z =0.5处,导热系数非连续,并存在强的间断㊂K =K 1=13, 0<z ɤ0.5K 2=13ˑ106, 0.5<z <ìîíïïïï1(22) 本问题具有稳态解析解T =1+8K 212K 2+K 2-K 1192K 1K 2-112K 1z 4, 0<z ɤ0.51+8K 212K 2-112K 2z 4,0.5<z <ìîíïïïï1(23) 采用新算法在正交网格和K e r s h a w 网格(图2所示)进行计算㊂图3给出了在不同规模的正交网格和K e r s h a w 网格(4ˑ4ˑ4,6ˑ6ˑ6,8ˑ8ˑ8,16ˑ16ˑ16,32ˑ32ˑ32)上数值解与解析解的相对误差,并将误差ε以网格平均尺寸(取对数)h 为自变量作线性拟合㊂可以看出,在正交网格上,数值解的误差接近2阶(1.9565)的速度收敛;在K e r s h a w 网格上,收敛的精度低于2阶(1.5569),这是因为K e r s h a w 网格是典型的大变形扭曲网格,并且随着网格的加密,变形扭曲程度会更加严重,网格品质也急剧下降㊂从图中看出,对于介质导热系数相差6个数量级的强间断扩散问题,在强扭曲的网格上,新算法仍具有高于1阶的精度,说明新算法具有较强的适应性㊂F i g.2 K e r s h a w g r i d 图2 K e r s h a w网格F i g .3 F i t t e d c o n v e r g e n c e r a t e i no r t h o go n a l a n dK e r s h a w g r i d 图3 正交网格和K e r s h a w 网格计算的拟合收敛阶3.2 强非线性热传导问题点源三维扩散问题本节考察一个具有强非线性导热系数的问题,导热系数K 为温度T 的α次方㊂这里采用第2节的点源三维扩散问题为例,该问题的控制方程为式(7)㊂此种情况下,热波以球对称的形式向四周传播㊂温度的空间分布应满足能量守恒条件,设初始总能量为ε0(量纲:M L 2T -2),根据式(8),有ʏ+¥-¥T 4πr 2d r =T 0=ε0/ρ()c (24) 该问题存在精确解:热波位置r f 和中心温度T c 分别为r f =ξ1b T α0()t 1/(3α+2)(25)T c =ξ13αξ122(3α2éëêêùûúú)1/αT 0r f3(26)强激光与粒子束郭少冬等:含强间断导热系数的三维扩散方程数值计算式中:b=1/ρ()c;ξ1是与α有关的常数㊂温度分布由式(27)给出T=T c1-r2/r f()21/α(27)在计算中,取α=3(导热系数是温度的3次方),总能量ε0=1.0J,计算区域为[0,1.2m]ˑ[0,1.2m]ˑ[0,1.2m],采用48ˑ48ˑ48的正交六面体网格离散,网格尺寸Δr=0.025m㊂如第2节所述,在初始时刻温度T c≫1K,导热系数相差27个数量级的情况下,传统支撑算子格式会产生 数值热障 ,导致热波无法前传,计算结果失真,而新方法可以保证热波有效的传播㊂图4是用新方法计算得到的2个不同时刻的温度分布情况(t=0.0002s,t=0.3s),图中的横坐标r表示到中心(能量源)的距离㊂F i g.4 A n a l y t i c a n dn u m e r i c a l s o l u t i o n s o f t e m p e r a t u r e a t t w od i f f e r e n t t i m e p o i n t s图4在2个不同时刻温度的解析解与数值解计算结果显示,在0.0002s时中心温度的数值解为433.3K,精确解为428.5K,相对误差为1.12%;在0.3s时,中心温度的数值解为0.6687K,精确解为0.6646K,相对误差仅为0.628%;此时,中心温度已经下降了4个数量级,计算得到的温度和热波位置与精确解符合得依然很好㊂4结论本文给出了一种改进的支撑算子方法,用于数值模拟辐射流体力学中强间断和强非线性的热扩散问题㊂采用两个具有强间断㊁强非线性的三维扩散算例,对算法进行了考核㊂计算结果表明,对于强间断多介质扩散问题,新方法在正交网格上可以保持2阶精度,在大变形非光滑网格上可以保持高于1阶的精度;对于强非线性热波扩散问题,新方法可以有效地克服传统方法的 数值热障 ,计算结果与解析解符合很好㊂因此,新的改进算法适合于辐射流体力学领域的强间断扩散问题的模拟㊂参考文献:[1]袁光伟,杭旭登,盛志强,等.辐射扩散计算若干研究进展[J].计算物理,2009,26(4):475-500.(Y u a nG u a n g w e i,H a n g X u d e n g,S h e n gZ h i q i a n g,e t a l.P r o g r e s s i nn u m e r i c a lm e t h o d s f o r r a d i a t i o n d i f f u s i o n e q u a t i o n s.C h i n e s e J o u r n a l o f C o m p u t a t i o n a lP h y s i c s,2009,26(4): 475-500)[2]裴文兵,朱少平.激光聚变中的科学计算[J].物理,2009,38(8):559-568.(P e iW e n b i n g,Z h uS h a o p i n g.S c i e n t i f i c c o m p u t i n g f o r l a s e r f u-s i o n.P h y s i c s,2009,38(8):559-568)[3]勇珩,高耀明,齐进,等.数值模拟辐射驱动2维内爆压缩过程[J].强激光与粒子束,2010,22(10):2340-2344.(Y o n g H e n g,G a oY a o m i n g,Q i J i n,e t a l.N u m e r i c a l s i m u l a t i o n o f2Dr a d i a t i o n-d r i v e n i m p l o s i o n p r o c e s s.H i g hP o w e rL a s e r a n dP a r t i c l eB e a m s,2010,22(10):2340-2344)[4]郁晓瑾,叶文华,吴俊峰.直接驱动球内爆点火的数值模拟[J].强激光与粒子束,2006,18(8):1297-1301.(Y uX i a o j i n,Y eW e n h u a,W uJ u n f e n g.N u m e r i c a l s i m u l a t i o no f d i r e c t-d r i v e I C Fi g n i t i o n i ns p h e r i c a l g e o m e t r y.H i g hP o w e rL a s e ra n d P a r t i c l eB e a m s,2006,18(8): 1297-1301)[5]李双贵,翟传磊,杭旭登,等.激光腔靶耦合二维辐射输运数值模拟[J].强激光与粒子束,2014,26:082002.(L i S h u a n g g u i,Z h a iC h u a n-l e i,H a n g X u d e n g,e t a l.T w o-d i m e n s i o n a l r a d i a t i o n t r a n s f e r c a l c u l a t i o no f l a s e r-t a r g e t h o h l r a u m.H i g hP o w e rL a s e r a n dP a r t i c l eB e a m s, 2014,26:082002)[6] S h e s t a k o vAI,P r a s a d M K,M i l o v i c h JL,e t a l.T h e r a d i a t i o nh y d r o d y n a m i c I C F3Dc o d e[J].C o m p u tM e t h o dA p p lM,2000,187(1/2):181-200.[7]S h a s h k o vMJ,S t e i n b e r g S.S o l v i n g d i f f u s i o n e q u a t i o n sw i t h r o u g h c o e f f i c i e n t s i n r o u g h g r i d s[J].J C o m p u t P h y s,1996,129(2):383-405.[8] M o r e l JE,R o b e r t sR M,S h a s h k o vMJ.Al o c a l s u p p o r t-o p e r a t o r s d i f f u s i o nd i s c r e t i z a t i o n s c h e m e f o r q u a d r i l a t e r a l r-z m e s h e s[J].JC o m-强激光与粒子束p u tP h y s,1998,144(1):17-51.[9] L i p n i k o vK,M o r e l J,S h a s h k o vMJ.M i m e t i c f i n i t e d i f f e r e n c em e t h o d s f o r d i f f u s i o n e q u a t i o n s o nn o n-o r t h o g o n a l n o n-c o n f o r m a lm e s h e s[J].JC o m p u tP h y s,2004,199(2):589-597.[10]郭少冬,周海兵,张树道.适应强扭曲网格的扩散方程时空二阶精度计算[J].强激光与粒子束,2012,24(11):2618-2622.(G u o S h a o d o n g,Z h o uH a i b i n g,Z h a n g S h u d a o.C a l c u l a t i o no f s e c o n d-o r d e r s p a t i a l a n dt e m p o r a l a c c u r a c y f o rd i f f u s i o ne q u a t i o no ns t r o n g l y d i s t o r t e d m e-s h e s.H i g hP o w e rL a s e r a n dP a r t i c l eB e a m s,2012,24(11):2618-2622)I m p r o v e dn u m e r i c a lm e t h o d f o r t h r e e d i m e n s i o n a l d i f f u s i o n e q u a t i o nw i t h s t r o n g l y d i s c o n t i n u o u s c o e f f i c i e n t sG u oS h a o d o n g, Z h a n g M i n g y u, Z h o uH a i b i n g, Z h a n g S h u d a o(I n s t i t u t e o f A p p l i e dP h y s i c s a n dC o m p u t a t i o n a lM a t h e m a t i c s,B e i j i n g100094,C h i n a)A b s t r a c t: T h en u m e r i c a l h e a t b a r r i e r i s i n d u c e dw h e n t h e t r a d i t i o n a l s u p p o r t o p e r a t o rm e t h o d i su s e d t o s o l v e t h e d i f f u s i o n e q u a t i o nw i t hs t r o n g l y d i s c o n t i n u o u s c o e f f i c i e n t s.I n t h i s p a p e r,a n e ws u p p o r t o p e r a t o rm e t h o d i s p r o p o s e d.I n t h i sm e t h o d,a n i-s o t r o p y c o n d u c t i v i t y t e n s o r s a r e i n t r o d u c e d a n d d e f i n e d o n g r i d v e r t i c e s.T h em e t h o d i s a p p l i e d t o t w o d i f f u s i o n c a s e sw i t h s t r o n g-l y d i s c o n t i n u o u s c o n d u c t i v i t y c o e f f i c i e n t s.T h e n u m e r i c a l r e s u l t s a g r e ew e l lw i t h t h e a n a l y t i c a l v a l u e s,w h i c h p r o v e s t h e f e a s i b i l i t y o f t h e p r o p o s e d s c h e m e t o s o l v e t h e d i f f u s i o ne q u a t i o nw i t hs t r o n g l y d i s c o n t i n u o u s c o e f f i c i e n t s.K e y w o r d s:i n e r t i a l c o n f i n e m e n t f u s i o n; d i f f u s i o ne q u a t i o n; s u p p o r to p e r a t o rm e t h o d; s t r o n g l y d i s c o n t i n u o u sc o e f f i-c i e n t s;s t r o n g n o n l i n e a r i t yP A C S:52.57.-z;47.70.M c;44.05.+e;02.70.B f。