（粒子物理与原子核物理专业论文）离散纵标法求解含有各向异性散射的输运方程.pdf

摘要 摘要 各向异性散射问题的理论计算，是反应堆物理、天体物理、大气辐射输运的 重要研究课题。由于输运方程的复杂性，除特别简单的情况外很难用解析方法求 解。对于一般的实际问题，都是采用数值方法求解。离散纵标法对所有的自变量 都采用直接离散，数值过程比较简单，因此在输运闫题的理论研究中被广泛应 用。 离散纵标法特指方向变量的离散，首先选定一组离散方向，每个方向对应一 定的立体角，所有的立体角之和为4n 。在特定方向上求解输运方程，假定每个 立体角内的角通量是常数，采用求积代替积分。在有限个离散方向条件下，角通 量高阶成分的求积误差造成了射线效应和总通量随空间位置振荡。 物理分析和数值结果表明，含各向异性敏射的一维平几何、球几何、有限跃 圆柱模型内，角通蠡的方向分布函数是分区光滑连续的，在某一方向( 或几个方 向) 存在极大。角通量各向异性越严重，极大方向上对角通薰各次矩的贡献份额 越大。利用计算得到的不同方向的角通量，从整体上构造角通量高阶部分，修正 求积误差，将会显著地提高整体计算精度。 对于一维平板各向异性散射问题，根据其存在8 函数特解的性质，我们在计 算源项时，把角通量函数分解成一个有限宽度的6 函数( 高阶) 和一个剩余函数， 对高阶函数用积分、低阶函数用求积确定。数值结果表明，较少的离散方向，能 够达到较高的计算精度。 本文分为六章，在简要说明了输运理论的基本概念之后，给出了微分输运方 程和积分输运方程，推导了输运算符在不同正交曲线坐标系中的表达式。然后介 绍了输运方程的几种求解方法。第四章，我们着重分析了各种空间和角度差分格 式的计算精度和求积组的构造。数值计算结果说明，采用精度高的差分格式和求 积组后，对各向同性和低阶各向异性散射的计算精度有一定的提高，但是对于高 阶的各向异性散射问题的计算，结果仍然不理想。在第五章，针对一维平板各向 异性散射问题，讨论了解析离散纵标法的概念和应用。最后，讨论了后续工作的 方向。 关键字：输运方程；离散纵标法( s n 方法) ：备向异性敖射 摘要 a b s t r a c t a n i s e 订o p i cs c a t t e r i n gi smi m p o r t a n tr e s e a r c ha r e amr e a c t o rp h y s i c s 神o n o m i c a lp h y s i c s a n dr a d i a t i o nt r a n s f e r f o rt h ec o m p l e x i t yo ft r a v a p o r te q u a t i o n , i t sh a r dt ob ea m a y t i cs o l v e d e x c e p tf o rs p e c i a lc o n d i t i o n s d i s c r e t eo r d i n a t em e t h o di sas i m p l ep r o c e s sa n dw i d e l ya p p l i e di n n u m e r i c a ls t u d yo f t r a n s p o r ti np r a c t i c a lp r o b l e m s i nt h ed i s c r e t eo r d i n a t em e t h o d ( d o m ) t h et r a n s p o r tc i i 噙d o ni ss o l v e df o ras e to fm d i s c r e t ed i r e c t i o n se a c hd i r e c t i o n 氖i sa s s o c i a t e dw i t has o l i da n g l ei 1 1w h i c ht h em t e n s i t yi s a s s u m e dt ob ec o n s t a n ta l ls o l i dw 塔l e sm n o n - o v e r l a p p i n ga n ds p a n n l 】a gt h et o t a la n s l er a n g eo f 4 nt h ei n t e g r a l so v c rt h ed i r e c t i o na r er e p l a c e db yn u m e r i c a lq u a d r a t u r ew e i g h t s ”ma s s u m e d o v e re a c ho r d i n a t ei ne a c hd i 栅o dt h et r a n s p o r te q u a t i o nm a yb es o l v e da n a l y t i c a l l y “e w i t h o u tt n m c a n o ne f f o r sd u et os p a t m ld i f f e r e n c i n g ) ，a n dt h es o i u u o n sa r ee x a c ti nt h ed i r e c t i o n 5 。t h er a ye f f e c ta n o m a l i e sa r i s e ，h o w e v e r , n o tf r o mt h ei n a b i l i t yt oc a l c u l a t et h ea n g u l a rf l u x i nt h ed i s c r e t eo r d i n a t ed i r e c t i o ne x a c t l yr a t h e rt h ea p p e a rf r o mt h ei n a b i l i t yo ft h eq u a d r a t u r e f o r m u l a i n1 - dp l a n g e o m e t r y , s p h e x i c a lg e o m e t r yo rf i n i t ec y l i n d e rg e o m e t r yi n c l u d i n ga n i s o l z o p i c s c a t t e r i n g ，p h y s i c a la n a l y s i sa n dn u m e r i c a lr e s u l t ss h o wt h a tt h ed i s t r i b u t i o nf u n c t i o no fa n g u l a r f l u xi ss l l q 0 0 t ha n dc o n l j n u o u 3 w i t hm a x m a u mi no n eo rt w od i r e c t i o n s n l eh i g h e rw a st h eo r d e r o fa m s o t r o p i cs c a t t e r i n g ，t h em o r ec o n t r i b u t i o nt h e g u l a rf l u xi nt h ed i r e c t i o mo fp e a kd o n e m a k i n gg o o du s eo fa n s u l a rf l u xi nd i f f e r e n td i r e c t i o m ，s 加l c t i m n g 出el l i g ho r d e rp a r to f t h e a n 则a rf l u xi nw h o l et or e v i s et h ee r f o r smq u a d r a t u r e ，w ec a l lm 驴r o v et h ec u l c u i a t i n ga c c u r a c y o i lt h ew h o l e f o rt h ea m s c - 口o p i cs c a t t e r i n gp r o b l e mi n1 - ds l a b ，a c c o r d i n gt ot h es p e c i a ls o l u t i o no fd e l t a f u n c t i o n , w es p l i tt h ea a a g u l a rf l u xf u n c t i o ni n t oa 矗n 沁w i d t hd e l t af i m c t i o n 。1 曲o r d e r ) a n da r e m a i n e df u n c t i o nt h e nt h eh i g ho r d e rf u n c t i o ni si n t e g r a t e dm a dt h el o w e ro r d e rf o n c t i o n q u a d r a t u r e n u m e r i c a l r e s u l t ss h o w t h a t w e c a n g e ts a t i s f i e d a c c u r a c y w i t h l e s s d i s c r e t e d i r e c t i o n s t h i sp a p e rh a ss i xc h a p t e r s a f i e rs i m p l yi n t r o d u c e dt h ee s s e n t i a lc o n c e p t so f t r a n s p o r tt h e o r y , w eg i v et h ed i f f e r e n t i a lt r a n s p o r te q u a t i o na n di n t e g r a le q u a u o n , a r md e r i v a t et h ee x p r e s s i o n so f t r a n s f e ro p e r a r t o l si nd i f f e r e n to r t h o g o n a lo r d i n a t e st h e nw es h o ws e v e r a ls o l u t i o nm e t h o d so f t r a n s p o r te q u a t i o ni nc h a p4 ，ma n a l y z et h ea c c u r a c yo fd i f f e r e n ts p a t i a la n d 酬a r d i f f e r i a l f o r ma n d t h ec o n s t r u c u o n o f q u a d r a t u r e f o r m u l a a f t e r t a k i n g h i s h a c c u f a c y d i f f e r t i a l f o r m sa n d q u a d r a t u r ef o r m u l a , n u m e r i c a l r e s u l t ss h o wt h a tt h ec o m p u t i n ga c c u r a c yo f i s o t r o p i co rl o wo r d e r a n i s o t r o p i cs c a t t e n n 8h a ss o m ep r o g r e s s , h o w e v e r , i th a k e $ p o o rm s p r o v e m e l q ti nh i g ho r d e r a m i s o t r o p i cs c a t t e r i n 塔p r o b l e m i nc h a p5 ，w ed i s c u s st h ec o n c e p t sa n da p p l i c a t i o mo fa n a l y s i s d i s c r e t eo r d i n a t em e t h o d sb y1 - ds l a bw i t ha n i s o w e p i cs c a t t e r i n g i nt h e 毗t h en e x ts t e pi s d i r e c t e d k e yw o r d s ：t r a n s p o r tt h e o r y ；d i s c r e t eo r d i n a t em e t h o d s f i ，e ) 中( ；，西，e ，t ) d 西d e + g ( ，e ，孬，f ) ：从上 面得到积分输运方程的过程可知，微分形式的输运方程与积分形式的输运方程 是等价的。对于各向同性散射和各向同性源，并对所有方向积分方程( 2 5 ) 得到： 痧c i ，e ，r ，= 皇里笔翥譬兰等幸尘蜴矿f 【j c ；，e e 卜 ，：，。、 m f ，e , t - - 止二丛) d e + 9 6 t ，e ， 一止二丛) ( 2 6 ) 式中e x p - r ( e ，；斗；) 】4 丌l ；一；f 表示；处出现的中子到达；处而不受碰 撞的几率；r ( e ，；f _ ；) 称为光学程长。 第二章粒予输运方程 2 2 输运算符 输运方程是根据物理中的中子守恒条件在有限体元和运动方向立体角元上 建立的，当体元和立体角元趋向无穷小时便得出微分输运方程。在数值求解输 运方程时，体元和立体角元是有限大小的，因而处理的是差分方程，如果能在 有限体元上保持推导方程时建立的守恒关系 1 】【乳，即源项等于随时间变化项、碰 撞移出项与表面流项之和，这样的差分方程就是守恒形式的。输运方程中表面 流与所围体积内散度积分之间满足高斯定理 l 矗垂瘩= v 西m d 矿( 2 - 7 ) 这时v 矗中是酗的散度，在曲几何下，由于方向角度和空间变量有关， v f 洒和西v 中形式是不同的。但是无论在笛卡尔坐标还是在曲线坐标下均有 v f 洒= 磊v o ( 2 8 ) 实际上，因为v 硒= 矗v 西+ o v 西，而v ，磊= 0 。笛卡尔坐标下可以直 接由定义得出；曲线坐标系下，由于矗向量应与坐标系选择无关，即西需要用 固定的坐标系来表示，所以也有v 西= 0 。 曲线坐标系下的散度可以写为 v = l 细丁= - 一札姗( 2 - 9 ) y “id v o s 2 2 1 正交曲线坐标系中的输运算符 在任意正交曲线坐标系下，设曲线坐标和微分线元分别为1 1 】 q t = q l ( x , y , z ) 幽= q 由，t 峨= ( 毒) 2 + ( 嚣 2 + ( 毒 2 了： 9 2 = q 2 ( x ，y ，：) td s 2 = h 2 由2 ， 面jlj n3 旦电， + nj 盟电，q l o) 堕电 一一 i i 第二章粒子输运方程 q 3 = q 3 ( x , y , z ) , ，_ ( 笥+ + 删j 瑚， 则体积元为d v = i s - s , 由l 由2 c q 3 。 设西= 垂( g 。，g ：，q 3 , , u ，) ，按( 2 - 9 ) 式求矢量殴数的散度v 五，先作一微小六 面体，由q lq l + 妃，9 2 ，吼+ 妃，9 3 ，吼+ d 毛六个面围成，因为是正交曲线，小 六面体可近似看成平行六面体。现计算各面的方向流量。通过g l ，q l + 由，两面的 万j 可执重，寺于 ( a i h 2 h 3 d q ：d q 引。幽一( a ，h 2 h 3 d q ：由引龟= 掣出，由：由3 。( 2 - 1 1 ) 同样，通过g ：，q 2 + 砌2 两面的方向流量和通过吼，q ，+ 由。两面的方向流量为别为 丛掣妃锄和塑笔婴嘲幽：妃( 2 - 1 2 ) 砌2口磊 将上面三式相加，再除以体元d ，= h 。日：h ，由l 由2 d 吼，并令由。由：d g ，趋向于零， 最后得到 v 永未百i 掣十掣+ 掣1 ( 2 - 1 3 ) 且也马i 却ta 9 2， j 因此有 v 西垂：一i 塑幽+ 塑幽+ _ o ( h t h 2 q 3 c b ) l h t 日2 爆l 抛ta 晚 j ( 2 - 1 4 ) v 两：! _ 型2 塑型+ 型! 坚盟+ 型擎鲤 甄玛、 却l却2电 虬拼变 + 1j 丝曼l 丝2 丝鱼竺2 + o t 垒! 些竖l 竺2 1 2 + 堕o ( h i h 2 9 d g ) l ( 2 15 ) i 却l舡 国2舡 却3础 j + 鬻堂笔竽+ 塑o q 塑袭掣+ 嚣塑秀竽p 2l l a 国 a 国固i a 国 j ， 第二章粒子输运方程 v 弧彘瓤芈+ 警+ 掣p k 变 + 云 嚣皿以q + 薏坞口，q ：+ 嚣且h ：q a i 垂c 。e ， a r ( a o h 卿h 。+ 豢乜且q + 罢一螭即 旦f 望 ：旦塑：o1 和l 。q 。ja g 。弘 旦f 丝 ：旦丝：o i 丽【瓦j _ 瓦瓦= 0 i 因为通常掣：o ，娑；o ，前式还可以简化为 0 9 l瑚l 此处“为选定的第一角变量。第一角变量也可能为f 。 上面就是一般的守恒形式的输运算符。 守恒形式输运算符带有空间偏导数的项与小体元( 和立体角元) 的乘积代 表中子流与表面积的乘积( 全微分) 。对小体元和立体角元记分后直接得出流出 表面的中子流积分值( 差值) ，即泄露。 带有方向角偏导数的项与立体角元( 和体元) 的乘积代表方向角度变化对 中子流的修正量与表面积的乘积( 全微分) 。在对立体角元( 和体元) 积分后代 表相关方向角度变化对泄露的修正。 虽然西v 4 ) 和v 磊。是等价的，但在曲线坐标下二算符的具体形式不同， v 西由具有以上守恒形式，而磊v 4 ) 不具有以上守恒形式。 下面给出几个实例： 盯，妲 r 陵g m p 一30 虬塑屯 哟一 h | 已 铲咿 峨一 翌 塑 皇a + 一v 蛆雌 忑一 0 一m擎去 31f，t，j盟j 竺 匕吼 印一 缸一庇 蹲一锄。一玛 也一 一 爿持 = 0 一址珏嗉 第二章粒子输运方程 例1 一般圆柱几何( ，o , z ) 的输运算符守恒形式 iq i h lq ： 1rl r = 1 一20 0 $ o ) 2臼r f = 4 1 一卢2 s i n c o 3z 1 中= 中( ，0 ，互从国) ( 对于选定角变量( ，国) 来说有国= z 一口，z 为q 在x - y 平面内投影与x 轴之间夹角) v - 面：上l 旦+ 型+ 亟型i 1 ，1l 静a 目 出j 。酃变 + 瓦0 ( 7 1 历0 t , f + ；警) 中 + 亳 ( ；等f + 警 中 即 v 磊o ：里螋+ 王旦! + “塑一三型( 2 - 1 8 ) r 却ra 口玉，a 国 例2 三维球几何( r ，0 ，妒) 的输运算符守恒形式 iq|ht q t lr1a = c o s 占 20r r = 4 1 一口2c o s 3rs i n 0f = 4 1 一2s i n 0 9 中= ( r ，护，伊；从0 9 ) ( 对于选定的角变量( ，国) ，有，= ( 占，妒) ，出= 出( 口，p ) ， v = ，2s i n 8 d r d o d o ) v 硒：一i o ( r 2s i n 0 庙) + o ( r s i n o r c ) + 旦型0 ，2s i n e 务0 0 却m 挪变 + 旦o p 肥l k r 塑o o 即+ 未历嚣f 卜 + 旦a c 阿l l ! r 塑a o 彳+ 志嚣彳卜 由坐标变换关系可以得到 竺：玎，娑：f s i n 0 ，塑：一s l n o ) ，塑：一c o s 8 一s i n e ”c o s c o 8 e 。 a 母。 0 0 7 s i ns 8 s i l l8 v 矗由：阜曼! 型+ r 。却 堡一型呈壁2 。上地 r s i n o 0 0r s i n oa 口+ 土业塑一一c o t 0 0 ( 圆) ( 2 1 9 )ro “r0 国 第二章粒子输运方程 2 2 2 一般坐标系下输运算符的具体形式1 - 3 1 1 s l 下面给出一般坐标系下输运算符矗v 西和v 盎。的表达式 1 直角坐标系 圈2 2 ：亘角坐标系 中子位置矢量尹：x , y , z 中子位置矢量磊：口， 这里，= 磊瓦，乏是z 轴方向的单位矢量；0 9 是矢量磊和t 所形成的平面与 矢量乏和t 形成的平面之间的夹角：t 和乞分别为x 和y 轴方向的单位矢量。 易见 q ：= 西已= q ，= 壶t = 正6 0 s 0 9 = 卵 1 4 第二章粒子输运方程 q ，= 西弓= 正s i n t o = f 一维平板：r b ( z ， 磊v 由：“塑， 以 二维x - y 几何：由( j ，y ，麒国) f t v m ：，7 垫+ f 塑， 出洲 三维长方体：由( j ，y ，z ，i t ，6 9 ) j 施= 2 万f ，d 1 t f 施：f 1d f 2 d 国 j j l j 0 f t v m ：节娑+ f 罢+ 罂 0 w洲出 直角坐标系下叠v 中和v 壶西形式是相同的。 2 柱坐标系 x 肛= f 。d e ”d o ) 圈2 3 ：柱坐标系 ( 2 - 2 0 ) 一一 兰三兰墼至塑重壹塑 中子位置矢量尹：，目，2 中子位置矢量矗：掣， 其中，护是极角；= 西t 乏；印是矢量磊和乏所形成的平面与矢量乏和位置矢 量尹所形成的平面之间的夹角。可得 r = ( z 2 + y 2 ) 2 船，j = r c o s o ，y = r s i n o ，曰= a r c t a n ( y x ) q ；= 磊砭= 卢 q ，= n 弓= ( i - 2 ) ”c 0 5 0 9 = 叩 q 口= 矗毛= ( 1 一2 ) ”s i n m = f 一维圆柱：西( ，岸，国) 衲中= 玎娑o r 一f 娑o c o ， 归= d - ir 。如r j 0 u 守恒形式v 两= 1 1o ( r c d ) 一! 型 r o r，口功 二维，口几何：由( ，0 ，2 ，曲) 西v o = 玎詈+ 事嚣一手慧， f 施_ f l j - 1 抛 j p 0 妇 。a rra 8r8 3 ” 守恒形式v 西西= 里型+ 生翌一! 燮 ，o rrd ， 刁国 二维，r y e 何：o ( ，：，a ，妫 a v 中= 疗娑o r + 罢一生r 罢， j 诵= j = ，d r 5 d 国 0 2d o j + jj o 守恒形式甲硒= 翌型+ a 塑一三型 ， o r出，d 三维，口，z ：中p ，占，= ，珊) f i v 母= 呀詈+ 手嚼一杀坶+ 警，f 施= f ，d r 8 拍 守恒形式v 西中= 翌螋+ 兰塑+ a 塑一1 0 ( 4 0 ) ro rra 目 色r8 ( 2 - 2 1 ) 第二章粒子输运方程 3 球坐标系 中子位置矢量尹：，矽，p q 一 圈2 4 ：球坐标系 中子位置矢量矗：从 其中，臼是极角，妒是方位角；= 磊i ；国是矢量矗和i 所形成的平面与矢量弓 和e 所形成的平面之间的夹角。可得 q ，= 矗i = c o s = u q 口= 西弓= ( 1 一t 2 ) “2 c o s = 叩 = 西弓= ( 1 - , u 2 ) ”2s i n w - - - c 一维球对称：西( ，a ) 衲中= 娑o r + 丝r 鲁，j 在2 霄胁a “ 。“ 守恒形式v 硒：尝亟拿塑+ 三叠坠；垒塑 r o rr弛 第二章粒子输运方程 二维，0 几何：o ( r ，毋，f ，) 两v o = 娑0 7 + 罟詈+，口 j d 磊= j 二，d r 5 d 缈 守恒形式 v 西m ：- - i - - = - - - 笔o ( r 2 0 ) + r d r 三：维，0 ，妒：中( ，目，仍，国) 1 一“。锄 r 8 n f c o t 口a ， d 脚 !l曼!竺壁!。三堡坚二岂：壁)cot00 ( g ) rs i n 8a era 越ra 6 ) 西v 西：“丝+ 翌塑+ l 塑+ 丝塑一g c o t o 塑 。静ra 8r s i n 88 母r 8 h r0 j 施= f 。d r 8 d o ) 守恒形式 v 偷= 笋宰+ 叩 ，s m o o ( s i n o o ) f 锄1a 【( 1 一2 ) m 。一。1。1。一 a pr s i n o a 汐 c o t 臼a ( 喾函) rd ( 2 2 2 ) 1 8 第二章粒子输运方程 2 4 散射函数 总的来说，相函数有以下性质 各向同性介质中，散射相函数仅依赖于两个方向和的夹角口 厂_ 国) = f ( c o s = f ( a ) 粒子数守恒或能量守恒 去f 他斗c o ) d _ 1 c o s ( 一毋) = c o s 0 下面是辐射输运中一些散射相函数模型： ( 1 )( 各向同性散射) 常函数 f ( c o s 0 ) 。专 粒子尺寸远小于光波长。 ( 2 ) 最基本的各向异性散射 f ( c o s 0 ) = l + a c o s 8 向前散射多于向后散射。 ( 3 ) l a m b e r ts u r f a c e s ，( c o s 口) ：婴( s i n p + 仞一p ) c o s 目) 根据l a m b e r t sl a w 的球面反射。 ( 4 ) r a y l e i g hs c a t t e r i n g f ( c o s = 言( 1 十c o s 28 ) 衍射效果占优势，粒子半径小于波长五，r 2 n 的项时，叫p n 近似，将被截断的展开式带入输运方程，可以求得展 开系数所满足的一系列微分方程。在适当的边界条件下求解这微分方程组，便可 得出问题的近似解。 球谐函数法一般用来求输运方程的的近似解析解，但在一定条件下( 例如在 补救离散纵标法的射线效应时) 也可用作求解输运方程的数值方法的基础。 下面以一维平几何为例来说明球谐函数法，这时，由于对称，方位角不出现。 球谐函数k ( q ) 约化为勒让德多项式。 第三章求解输运方程的方法 塑笔半蝎( 恻训) 钱r 彬。缸硒。，f i ) o ( x , d ) 悄训) ( 3 - 1 ) 式中f ( f i 。，西) 是散射概率函数，若命风= 西丘，可展开概率函数： f ( f i 。，f i ) = 妻簪) ( 3 - 2 ) 式中岛= 2 万。d , u o f ( p o 谒( 1 2 0 ) ，b o = 1 ，b 。_ - 。 利用加法定理 异( 矗f i ) = 一讥4 n t r ，2 ( f i ) r ( f i ) ( 3 - 3 ) 删= 一1 - 1 撕m 舻【等。嵩渊m 南c o s 跏州 睁。) 斗_ 7 r “+ lj ! j d 壶( 西) ( 磊) = 气6 二 可简化( 3 - 1 ) 式中散射项： r 槲。咖f ( f i ，f i ) 吣) = 砉圭。咖啄肋砌x 2 ，+ 1 ) 州训i ) ( 3 5 ) 于是( 3 - 1 ) 可写为 罢+ ，西= 孚喜异( 麒2 z + 1 渤f 。咖啄们m ( z ) + s ( z 棚( 3 6 ) 将m ( x ，) 及s ( x ，) 分别对勒让德多项式展开 巾( t ) = 刍2 百l ；+ 面1 ，( x ) b ( ) ，其中q ( x ) = 2 石f 。d 声西( 肋毋( ) 脚) = 善号( 堋办其嘲加2 万i _ ! d a s 眦( 再利用恒等式只= 寺告与。+ i f _ 只一，就可以通过对( 3 6 ) 做运算，得出 去警+ 嘉警郴，恻垆鼢l 删3 ，( 3 - 7 ) 这组无限多个微分方程和原来的输运方程等价。 存p n 斤m j 】审彝素 第三章求解输运方程的方法 塑盥：o d x 这为n + 1 个未知量舢，巾l ，辔2 ，抵给出n + 1 个微分方程。 以p ，近似为例，有方程组： 华+ ( ，一z 。) 中。：s o ， ；拿+ ( 。一瓦。) 。：s j 戤 记中。= ，o ，= j 并假设各向同性源，因而s = o ，上述p 1 方程组就化为扩散 方程： 一丢c 击剐= 氐 其中，= ，一_ s ，z 。= z ，一。 为考察一般p n 方程组解的结构，置西，( j ) = a l e ，带入方程( 3 7 ) ，可得 出确定k 的特征方程： l - c 曙。k 00 芷3 2 k0 02 足5 + 3 k 可以证明当n 为奇数时，这方程有n + 1 个根士玛嘲。与此相应的有解 当从奇阶到偶阶展开时，根的数目并不增加，而且偶阶的p n 展开会引起边 界条件的一些问题。困此一般只做奇阶( p 1 , p 3 , p 5 ，e t c ) 展开。 2p n 法对各向异性的处理 在上面中的推导中展开了散射函数和角通量，这种做法实际上可以处理各向 异性，是在实际计算种常用的方法。 3p n 方程和离散纵标方程的等价3 。4 拼】： x q 勺 。瑚 = 、，0 m 第三章求解输运方程的方法 m a r k 边界条件，相当于把真空看成纯吸收介质： 蚤( o ，鸬) = 零0 ，一硒) = o ，j = l ，2 ，( + 1 ) 2 这里m 是p n ( 儿) = 0 的正根。 可以证明： s ，。离散纵标法+ 高斯求积法= p n 法+ m a r k 边界条件。 如果在离散纵标法中，对于心陆定义 吣= 鬈等讯删= 萎n 蜊u ) 就可用离散纵标法定出p n 近似解。 其中蕺( z ) = 幼既e ( 如) 西( x ，心) 。 在阶数增多时，尤其对于多维复杂几何情况，球谐函数法将变得非常复杂和 困难。以p 3 近似为例，对最箍单的一维情况，也需要联立求解四个微分方程， 每个方程都含有耦合的三个未知函数，而对二三维圆柱几何系统，方程数目将达 到1 6 个，每个方程将耦合有7 个未知函数。这些方程组的数值求解将是非常麻 绶的。同时，当阶数不同时，求解的方程和数值过程都不相同，因此不可能编制 适用于不同阶数n 的通用程序，这给工程计算应用带来诸多不便。这些缺点严 重她阻碍了高阶球谐近似方法的应用。1 5 1 3 2 3 离散节块法 输运方程的离散节块法5 1 与输运方程的离散纵标法在角度的处理上是完全 相同的。在空间处理上，离散节块法以节块为单位，尺寸较大，因此空间上对通 量和源项采用高阶展开( 一般取二阶以上) ，对于多维问题采用横向积分方法， 将多维问题化为通过泄漏项互相耦合的一维问题。一维问题则通过沿特征线方向 的宜接积分方法求解。 离散节块法的要点是：从守恒形式的输运方程出发，通过横向积分转化成几 个一维的常微分方程；几个一维方程之间通过泄漏项作为源项的一部分耦合起 来：节块内的通量、源项和泄漏项用多次式展开，求解的未知量归结为通量展开 系数。 第三章求解输运方程的方法 3 2 4 m c 法 m o n t ec a r l o 法不同于前面所介绍的确定性数值方法，它是用来解决数学和 物理问题的非确定性的( 概率统计的或随机的) 数值方法。m o n t ec a r l o 方法 ( m c m ) ，也称为统计试验方法，是理论物理学两大主要学科的合并：即随机过 程的概率统计理论( j n 于处理布朗运动或随机游动实验) 和位势理论，主要是研 究均匀介质的稳定状态。它是用一系列随机数来近似解决问题的一种方法，是通 过寻找一个概率统计的相似体并用实验取样过程来获得该相似体的近似解的处 理数学问题的一种手段。运用该近似方法所获得的问题的解更接近于物理实验结 果，而不是经典数值计算结果。 普遍认为我们当前所应用的m c 技术，其发展约可追溯至1 9 4 4 年，尽管在 早些时候仍有许多未解决的实例。m c m 的发展归功于核武器早期工作期间l o s a l a m o s ( 美国国家实验室中子散射研究中心) 的一批科学家。l o sa l a m o s 小组 的基础工作刺激了一次巨大的学科文化的迸发，并鼓励了m c m 在各种问题中的 应用。 m o n t ec a r l o ”的名称取自于m o n a c o ( 摩纳哥) 内以赌博娱乐而闻名的一 座城市。 m o n t ec a r l o 方法的应用有两种途径：仿真和取样。仿真是指提供实际随机 现象的数学上的模仿的方法。一个典型的例子就是对中子进入反应堆屏障的运动 进行仿真，用随机游动来模仿中子的锯齿形路径。取样是指通过研究少量的随机 的子集来演绎大量元素的特性的方法。例如，( z ) 在t 2 x 0 ， ( m 4 n 2 + 1 ，一n ) 中。( 爿，) = 0 ， 0 ， ( m 21 ，n 2 ) ( 4 - 3 ) 西为z ，1 2 的函数，对角度用勒让德多项式展开，系数为x 的函数。 蚍舻m 2 n ，。+ l p ，( ) 五( 功 ( 4 _ 4 ) 五( 毒) ；2 1 中( t 功名( ) d ( 4 - s ) n _ o 时，r ( ) = 1 ，磊( 工) 为x 处的总通量。n = l 时，暑( ) = ， 最( x ) 为x 处 沿x 方向的流的总通量。 著利用( 4 1 ) 的求积公式，则有 吣棚= 蓑专囊( 旃x ) 睁e ) 第四章s n 方法概述 r 五( 工) = 2 万中( 石，膳) 只( “) ( 4 7 ) i = l 将区间【o ，x 】离散化为一组离散 r x k + 1 ，2 ，o x t l 2 x 3 2 o ) 咆圳= 篙咻。圳+ s 蕊a ，( 对约 0 ) ( 4 - 1 3 ) 4 1 3 各种空晦差分格式及误差l 对于平几何输运方程。由于它不存在角通量方向的导数项，其离散化处理比 曲几何输运方程容易。对于它的离散s n 方法求解，己提出很多计算格式。如何 评价一个计算格式，用以下五条标准鉴鄹比较合适： ( 1 ) 精确度： ( 2 ) 简单性； ( 3 ) 解的非负性； ( 4 ) 推广到其他几何的可能性； ( 5 ) 与迭代加速技巧的相容性。 平几何输运方程不存在角通量对方向余弦的导数项，下面讨论的各种计算格 式都是对空间变量而言，分析的误差也是空间差分的误差。并假定求轵公式 f 。中( 工，y = 芝中h ) 是精确的。 根据上面的五条标准，介绍以下几种格式： 1 菱形格式： 假定角通量m ，( z ) 在网格( x k 1 1 2 ，x k + m ) 内线性变化，就得到菱形格式。 中 。= ( m m 垃。+ 垂，2 ，。) ( 4 m 1 4 ) 第四章s n 方法概述 对于。 0 ，即1 n ；n 2 + 1 ，n 1 2 + 2 ， ，n ，可解出 ( 4 1 5 ) ：矿薯每：一惫蚤 均 由辔l 垃。= 辔i 2 ，+ l 一。可依次递推求得蚤3 2 o ，垂“，2 ，。 瓯= 学粪啪“幢一。如) ( 4 - 1 7 ) 误差：e i = o ( a x 3 ) ，即菱形格式的局部网格边界通量的误差是三阶的，整体误差 是二阶的 2 步函数格式 这个格式简单，非负，也容易推广至其他几何情况，缺点是精度不高。般 不单独使用它，有时做菱形格式的负通量修补格式用。 假定角通量在网格( 碓地，x j c + i 彪) 内是个常数，就得到步函数关系式 巾搠= o 2 拼，掣卅 0( 4 - 1 8 ) 。矿一k + 1 2 , m 一惫妾，心 o # t ：矿州也冲k - l ：2 , m + 去瞄e x p 【- z k ( x k + l , , 2 - - x ) s ( 蚺( 4 用s ( x ) n 格平均值代入后积分，得 吼+ 1 ，2 。：e 砸一k ) m 。缇。“1 一e x p ( 一k ) 】拿，心 o ( 4 - 2 5 ) 从上式求出网格边界通量值吼“：。后，通过 中。= 导一士( ：，m + 幢。) ( 4 - 2 6 ) 一o 。b 1 计算网格平均通量值面。从函。又可计算下一次迭代源的网格平均值 最：! 毕艺w i 嚷( 小玑( 4 - 2 7 ) ( 4 - 2 5 ) ( 4 2 6 ) ( 4 2 7 ) 三式组成了对于肛， 0 的步特征线格式的计算公式。 4 1 第四章s n 方法概述 对。 0 ，代替( 4 - 2 5 ) 1 j 寸有 m 。一。，：，= e x p ( 一。k 。) 。，2 ，。+ 【1 一e x p ( 一二。) 】桌o ( 4 - 2 8 ) z - , k ( 4 - 2 8 ) ( 4 - 2 6 ) ( 4 - 2 7 ) g j t 成了对于卢， 0 得 吆m r 卅= e x p ( k ) 西k1 2 ，d - 去：c x p 【_ t ( x k + 1 1 2 - - x ) 心似鼻( 4 3 0 ) 然后形m 柳代八r 侍 1 1 ) k + 1 1 2 , m ：e x p ( 一k ) 。垃，“1 一e x p ( 一k ) 】导+ 垒磐【1 一( 百1 + 士一e x p ( 一k ) ) 】 - k- t - - h ( 4 - 3 1 ) 由上式解出网格边界角通量值后，再用h = 导一( o 叠，+ 母缇。) 诗 _ t1 j 锄 算角通量的网格平均值，用& ：_ ! 譬芝中。( 工) + 吼计算下一次迭代源的网格 平均值s k ，用瓦= 訾计算源在网格内的斜率t k ，其中 ：华羔噙蚴( 4 - 3 2 ) 以k 是线性特征线格式关于“ 0 的计算过程。对于“。 0 ，有类似的计 第四孽s n 方法概述 算过程。 误差：e 5 = o ( x 5 ) ，即线性特征线格式的局部网格边界通量的误差是五阶的，整 体误差是四阶的。它是这里介绍的格式中精度最高的一种。 4 1 4 高斯和双高斯求积公式 一维时，求积公式的选取满足以下三个基本原则： ( 1 ) 权重因子w m 必须是非负的。因为积分通量是恒大于等于零， m o ) = 圣) d o ，也要求中( z ) ：n 曲。) o 。 ( 2 ) 由x = o 处的反射边界条件要求，若“在 心 中，则一= 一“也应在 ) 中，且 4 = - a n 。，和嵋= w 。即( 心) 和 是关于= o 对称的。 ( 但是在已知通量是既不对称的特殊情况，使用非对称求积组是有利的，比如已 知通量峰接近于= + l 处，使离散方向密集于此附近将是有效的。) ( 3 ) 当巾( 也) 是的低阶多项式时，求积公式必须给出积分的精确值，即满足矩 条件 善nw f = 黔州= 鬲2 ，n 为偶数 nq “一：。肛一咖，：o ， n 为奇数 高斯求积组是满足以上要求的，在一维时它具有最高代表精度。我们计算中 的离散方向h 和求积权重w ，采用高斯求积组。考虑区间【- l ，1 】内有个 求积点的单高斯求积组，则求积公式对于f 2 n 一1 的多项式求积是严格精确的， 即具有2 n i 阶求积精度。 只( “) = 0 ，有n 个解，i = l ，2 ，3 ，n 。 2 1 一。_ - _ _ _ _ _ _ _ _ _ _ _ - _ _ - - _ _ _ _ - _ 一 ，z 0 一，( 鸬) e ( “) ( 4 - 3 3 ) 第四章s n 方法概述 为了解决在= 0 处被积函数( 中子角通量密度) 的不连续性所可能引起的误 差，y v o n 提出将的高斯求积区间 - 1 ，l 】分成【1 ，o 】和 0 ，1 两个半区间，而把高斯 一勒让德求积公式分别应用于上述两个半区间上。这就是所谓双高斯求积组1 3 ， 也称为双球谐( d p n ) 求