小波分析在反应堆物理中的应用

《核反应堆热工数值分析》课程报告1小波分析简介及其在反应堆物理中的应用小波分析是近二十年来迅速发展起来的全新数值分析方法，其基本思想是将函数用小波基函数来展开，计算其展开系数。由于小波基函数具有优良的紧支性，使得小波基函数可以聚集在研究的任意细节，被数学家和工程师们誉为“数学显微镜”[9]。即小波基函数可以很好的逼近各种剧烈变化的函数，尤其是局部震荡的函数，具有非常高的逼近精度。另外，小波分析具有多尺度、多分辨的特性[3]，能够提供多种基函数作为有限元插值函数，由此构造的小波基单元可以根据实际需要任意改变分析尺度，使在变化梯度小的求解域用大的分析尺度，而变化梯度大的求解域则采用小的分析尺度。这是一种优于传统单元网格加密和插值阶次升高的自适应有限元算法，这种变尺度算法数值稳定性好、运算速度快、求解精度高。本文简要介绍了小波基函数的基本性质、求解小波基函数的展开系数的数值方法、以及小波基函数在反应堆物理中的应用，包括用小波基函数离散中子通量密度中的角度变量、共振能量区间的能量变量，并结合有限元方法的思想，将所有的变量都用基函数展开，然后带入变分方程，得出一个代数方程组，从而可以求出展开系数。1.1理论基础1.1.1小波函数a.小波函数基本概念小波函数()x是一个在(−∞，+∞)区间积分为零的函数[9]，（容许条件）()0xdx(0-1)称()x母小波；“小”指的是函数的局部性，积分为零说明()x含有波动性，而且正负两个方向的波动是均等的；小波函数经过平移、伸缩可以得到一族小波函数，即：,1()(),,0abxbxabaaaR(0-2)其中a反应了特定函数的尺度（即伸缩情况），变量b指明了它沿着x轴的平移位置；当a和b取一系列离散值，如𝑎=2𝑗,𝑏=2𝑗𝑘(𝑗,𝑘∈𝒁)时，即可以得到一族小波函数，即：《核反应堆热工数值分析》课程报告22j,k()2(2),jjxxkjkZ(0-3)在一定条件下{𝜓𝑗,𝑘|𝑗,𝑘∈𝒁}可以构成𝐿2(𝑹)的基，用它可以表示𝑓(𝑥),j,k,,j,k()()((),())jkjkjkfxdxdfxx(0-4)b.Daubechies小波在数值计算中，通常使用Daubechies小波[9]，Daubechies小波具有紧支性和正交性；Daubechies小波没有显式的数学表达式，其尺度函数和小波函数由以下两尺度方程给出：210()()(2)NNNNkxpkxk(0-5)122()()(2)NNNkNxqkxk(0-6)(1)1kNNqkpk(0-7)式中N为Daubechies小波的阶数。由于Daubechies小波没有显示的数学表达式，其尺度函数和相应的小波函数通常用数值方法以数表和曲线方式给出；N=6时，Daubechies小波的尺度函数𝜙𝑁(𝑥)和小波函数𝜓𝑁(𝑥)如下：（a）尺度函数（b）小波函数c.Daubechies小波性质（1）紧支性若小波函数𝜓𝑁(𝑥)在区间[a,b]外恒为0，则称在区间[a,b]上紧支，[a,b]为其支集；N阶Daubechies小波（记为DN小波）的尺度函数𝜙𝑁(𝑥)和相应的小波函数𝜓𝑁(𝑥)的支集分别为Supp𝜙𝑁(𝑥)=[0,2N-1]，Supp𝜓𝑁(𝑥)=[1-N,N]。紧支性反应了尺度函数和小波函数在时域和空域的局部化能力，支集越小的小波，局部化能力越小。《核反应堆热工数值分析》课程报告3（2）正交性Daubechies小波是正交小波，其尺度函数𝜙𝑁(𝑥)和母小波𝜓𝑁(𝑥)满足如下正交条件：∫𝜙𝑁(𝑥−𝑗)𝜙𝑁(𝑥−𝑚)d𝑥+∞−∞=𝛿𝑗,𝑚𝑗,𝑚∈𝒁(0-8)∫𝜙𝑁(𝑥)𝜓𝑁(𝑥−𝑚)d𝑥+∞−∞=0𝑚∈𝒁(0-9)（3）完备性Dohon已经证明，Daubechies小波是平方可积函数的无条件基，即任意平方可积函数可以表示为：,j,k,,j,k()()((),())jkjkjkfxdxdfxx(0-10)Daubechies小波尺度函数可以精确的表征出不大于N-1的幂级数；Daubechies小波的紧支性保证了以小波尺度函数作为插值函数构造小波单元时能以最少的单元自由度最大限度的逼近求解函数；小波的正交性保证了小波有限元所形成的的刚度矩阵是稀疏的，这将大大减少奇异性问题有限元分析的运算量。1.1.2小波数值方法小波数值方法是将函数用小波基函数展开，然后用某些数值方法求解其展开系数，包括小波加权余量法、小波伽辽金方法、小波配置法、小波彼德洛夫伽辽金方法、小波有限元方法[3,8]，下面简要介绍小波加权余量法。数理问题通常是求解某个区域中某种或几种场函数的分布，记为u，通常可以由微分方程及其边界条件来描述：𝑨(𝒖)={𝐴1(𝒖)𝐴2(𝒖)⋮}=0𝒖∈𝛀(0-11)𝑩(𝒖)={𝐵1(𝒖)𝐵2(𝒖)⋮}=0𝒖∈Γ(0-12)其中𝛀和Γ分别为求解区域和边界。在区域𝛀内与微分方程组完全等价的积分形式为：1122 dΩ()dΩ0TvAvAvAuuu(0-13)若(0-13)积分方程对任意𝒗都成立，则微分方程组(0-11)一定满足。同理：1122 d()d0TvBvBvBuuu(0-14)若(0-14)积分方程对任意𝒗̅都成立，则微分方程组(0-12)一定满足。式(0-13)加上(0-14)得到：《核反应堆热工数值分析》课程报告4 dΩ d0TTvAuvBu(0-15)若(0-15)对任意𝒗、𝒗̅都成立，则(0-11)和(0-12)成立。假设未知场函数𝒖可以采用近似函数来表示，近似函数是一族带有待定参数的已知函数：1ˆniiiuuNuNu(0-16)其中𝒖̃𝒊为待定参数，𝑵𝒊为试探函数。通常n取有限项的情况下近似解不能精确满足微分方程(0-11)和边界条件(0-12)，它们将产生残差R以及R̅：R=𝑨(𝑵𝒖̃);R̅=𝑩(𝑵𝒖̃)(0-17)用n个规定的函数来代替任意函数𝒗及𝒗̅，即𝒗=𝑾𝒋;𝒗̅=𝑾̅̅̅𝒋(j=1,2,⋯,n)(0-18)得出近似的等效积分形式：∫𝑾𝒋𝑇𝑨(𝑵𝒖̃)dΩ+∫𝑾̅̅̅𝒋𝑇𝑩(𝑵𝒖̃)dΓ=0;j=1,2,⋯,n(0-19)通过选择待定系数𝒖̃𝒊，强迫余量在某种平均意义上等于0；𝑾𝒋和𝑾̅̅̅𝒋称为权函数，余量加权积分为0就得到一组求解方程，用以求解近似解的待定系数𝒖̃𝒊，从而得到原问题的近似解；任何独立的完全函数集都可以用来作权函数，按照对权函数的不同选择就得到了不同的加权余量计算方法；1.配点法：𝑾𝒋=𝜹(𝒙−𝒙𝒋)即强迫余量在域内n个点上为0；2.子域法：在n个子域Ωj内，𝑾𝒋=𝑰，在子域Ωj外𝑾𝒋=𝟎;即强迫余量在n个子域内积分为0；3.最小二乘法：当近似解取𝒖̂=∑𝑵𝒊𝒖̃𝒊𝒏𝒊=1时，权函数𝑾𝒋=𝝏𝝏𝒖𝒋𝑨(𝑵𝒖̃);此方法实质是使得函数∫𝑨𝟐(𝑵𝒖̃)dΩ取得最小值；4.力矩法强迫余量的各次矩为0，以一维为例，𝑾𝒋=1,𝑥,𝑥2⋯∫𝑨(𝑵𝒖̃)dΩ=0;∫𝑥𝑨(𝑵𝒖̃)dΩ=0;∫𝑥2𝑨(𝑵𝒖̃)dΩ=0,⋯5.伽辽金法《核反应堆热工数值分析》课程报告5取𝑾𝒋=𝑵𝒋,边界上𝑾̅̅̅𝒋=−𝑾𝒋=−𝑵𝒋；即简单的利用近似近似解的试探函数序列作为权函数；∫𝑵𝒋𝑇𝑨(𝑵𝒖̃)dΩ−∫𝑵𝒋𝑇𝑩(𝑵𝒖̃)dΓ=0;j=1,2,⋯,n小波加权余量法中1ˆniiiuuNuNu，其中试探函数𝑵𝒊为小波函数，即场函数利用小波基函数来展开，利用加权余量法来求其展开系数。利用其它方法也可以求得其展开系数，比如有限元方法等，其基本思想与加权余量法类似，具体过程参考文献[3,8]。1.2小波方法在反应堆物理中的应用反应堆物理中要求解中子输运方程，中子输运方程是一个与角度、能量、空间、时间相关的7个自变量的微分积分输运方程；一阶稳态中子输运方程形式如下：t''0,,Σ,,,Σ(,,,),,,,sEEEdEEEEdQEΩrΩrrΩrΩΩrΩΩrΩ(0-20)其中𝜴=(cosθ,sinθcosφ,sinθsinφ),是单位球面上的点(θ,φ)，代表一个给定的方向，𝑑𝜴是球面上的面微分，𝑑𝜴=sin(𝜃)d𝜃d𝜑；Σt(𝒓,𝐸),Σ𝑠(𝒓,𝐸′,𝜴′→𝐸,𝜴),𝑄(𝒓,𝐸,𝜴)分别为宏观总截面，概率散射函数和中子源；中子角通量密度𝜙(𝒓,𝐸,𝜴)是需要求解的量，是一个关于空间𝒓=(𝑥,𝑦,𝑧)，能量𝐸,方向𝜴=(θ,φ)六个自变量的函数。1.2.1中子输运方程Daubechies小波角度离散当中子通量密度随角度的分布呈现出剧烈的变化时，比如目前新型反应堆中广泛应用的MOX燃料中，传统的角度离散方法很难对𝜙(𝒓,𝐸,𝜴)进行很好的逼近，而由于Daubechies小波正交紧支的特点，可以很好的逼近各种剧烈变化的函数，所以可以利用Daubechies小波离散中子输运方程中的角度变量[1,4,5]。采用Daubechies小波对𝜙(𝒓,𝐸,𝜴)的角度变量进行离散就是将𝜙(𝒓,𝐸,𝜴)用Daubechies小波展开,,,()EEnnnrΩrwΩ(0-21)其中()nwΩ为小波基函数，,Enr为展开系数，通常还要对能量和空间进行离散，若利用分群理论对能量进行离散，有限元理论对空间进行离散，则《核反应堆热工数值分析》课程报告6)),((,gkkEg,n,knnNwrrΩΩ(0-22)其中,,gErΩ为第g群的中子通量密度，()kNr为空间基函数，有限元理论中为有限元空间中的基函数，()nwΩ为角度基函数，若用小波展开法，()nwΩ为小波基函数，g,n,k为展开系数。通过一些处理，式(0-20)可以写为算子形式[10]：gLqg(0-23)右端项可以用基函数展开)),((,gkkqEqg,n,knnNwrrΩΩ(0-24)通过加权余量法或变分方法可以得到式(0-23)的变分方程（暂不考虑边界）：,)((,)gqLg(0-25)其中g属于解空间V，属于与解空间相关的空间W，当V与W相同时为伽辽金方法，当V与W不相同时，为彼德洛夫伽辽金方法[8]。设空间V的基函数为(())knwNΩr，空间W的基函数为'''()()knwLNΩr，其中'L为一个算子，不同方法中'L不同，在伽辽金方法中'1L，在最小二乘方法中'LL，其它方法[]也可以求出相应的'L，将(0-22)带入到变分方程(0-25)中得到''''''(()()((),())((),()))()kkkkkkqg,n,kng,n,knnnnnLwΩLwΩwΩNrNrNrNrLwΩ(0-26)令空间响应矩阵''1(),()()kkNrNLrSL，''2(),()()kkNrNrSL，角度相应矩阵'(),())(nnwAw，所以12SASAq(0-27)其中代表Kronecker积。在求解过程中矩阵1S、2S、A向量q，都是已知的，在每个能群中求解该线性方程组(0-27)可以得出展开系数g,n,k，从而得出中子角通量密度,,ErΩ。《核反应堆热工数值分析》课程报告71.2.2中子输运方程Daubechies小波能量离散由于共振核素在共振区的总截面有剧烈的震荡，造成中子通量密度在共振区也有剧烈的震荡，所以可以利用Daubechies小波尺度函数对共振能量段中子角通量密度的