辐射扩散方程组的应用领域十分广泛,包括惯性约束聚变(ICF)、磁约束聚变(MCF)、天体物理、高超声速、燃烧等.当辐射与物质间没有达到平衡时,用非平衡辐射扩散来近似.此时,若物质(包括电子、离子)自身达到平衡,则为两温方程组,即非平衡辐射扩散方程组；若物质自身没有达到平衡时,则是三温方程组(即辐射三温扩散方程组)或多群方程组.因为这类方程组具有系数强间断性、强非线性、高度耦合性及存在突变区等特点,所以数值求解非常困难.有效地数值模拟这类方程组是一个有意义的间题,引起了国内外很多学者的注意[1-11].有限体积元法(FVEM)也被称为广义差分法,控制体积法,盒式格式.这类方法被广泛应用于科学和工程领域,比如电磁学,计算流体力学,石油油藏模拟等.文[12-17]详细讨论了解线性椭圆方程的线性元、二次元有限体积法及收敛性.本文将以上针对椭圆方程的有限体积元法推广应用到了数值模拟非平衡辐射扩散方程组和辐射三温能量方程组.辐射扩散方程组的解在计算区域内局部变化非常剧烈,如果选用均匀的网格剖分,需要的单元数很大,求解效率较低,解决这类间题的一个有效途径就是使用自适应计算.文[18]基于后验误差估计提出了解非线性抛物方程的自适应算法,我们将这种方法推广改进,与线性元有限体积法相结合,应用于二温、三温方程组的数值模拟.同时,本文设计了一种基于能量相对变化的时间自适应方法.本文的主要内容具体由以下四个部分组成：1.解非平衡辐射扩散方程组的低阶元有限体积法考虑非平衡辐射扩散方程组其中,E是辐射能量,T是材料温度,k=0.01T5/2是材料的热传导系数,D是辐射扩散系数,选用不带限流项和带限流项的两种形式.光子的吸收截面为σα=(z3)/(T3),其中的原子质量数z是与材料相关的一个函数.能量交换受σ。控制,大的z值和小的T值导致了两个方程间更紧密的耦合关系.系数D,K和σα使得方程组具有强非线性.<1>两种低阶元有限体积法对于时间离散,使用向后Euler法.对于空间离散,采用低阶元有限体积法,包括三角网上的线性元有限体积法和四边形网上的等参双线性元有限体积法.取试探函数空间Uh为相应于三角剖分Γh的线性有限元空间或相应于四边形剖分Γh的等参双线性有限元空间.对三角形网,用重心对偶剖分Γh*；对四边形网,用平均中心对偶剖分Γh*.检验函数空间Vh取相应于对偶剖分Γh*的分片常数空间.对求解非平衡辐射扩散方程组的全离散线性元有限体积法中定义出现的重积分,分别采用两种不同的处理方法,得到的格式Ⅰ利用质量集中法,格式Ⅱ采用精确积分法.<2>单调性及修补技术令其中n表示时间层,3是Γh的节点数.因此,构造的全离散低阶元有限体积格式可以写为用冻结系数法对(3)式进行线性化,得到命题1.如果U0≥0且矩阵AT(U(n))是一个M-矩阵,则Un+1≥0(n≥0).由命题1知,如果我们可以提供一些条件,使得AT(U(n))是一个M-矩阵,那么该格式单调.对于格式Ⅰ,当剖分满足一定条件时(即对三角形剖分,所有原始单元非钝角；对平行四边形网格剖分,所有原始单元的邻边比在√3/3与√3之间,且最小内角不小于60°),能保证AT(U(n))是M-矩阵,从而格式单调.对于格式Ⅱ,通常不满足M-矩阵的性质.针对格式Ⅱ这一缺陷,我们采用了保证能量守恒的两种修补技术来弥补,包括局部修补法和整体修补法.修补技术的基本思想是：当某些单元的能量E(或温度T)小于限定的最小值时,在保持整个区域上总能量εtotal不变的情况下,所需能量由其他单元补给.其中,局部修补法是按单元逐个修补E和T.而整体修补法是同时修补所有单元上的E和T.数值模拟结果表明两种格式最终得到的数值解的守恒误差都很小,且温度、能量在计算过程中均保正.2.解非平衡辐射扩散方程组的高阶元有限体积法对非平衡辐射扩散方程组(1)-(2),采用向后Euler法进行时间离散,并利用冻结系数法线性化.对于空间离散,采用矩形网上双二次元有限体积法.使用矩形网上双二次元有限体积法求解二温方程组时,网格剖分可以比较稀疏；而用低阶元时则要求网格比较稠密.双二次元在64×64矩形网上的自由度等同于双线性元在128×128矩形网上的自由度.双线性元有限体积法在矩形网上是九点格式,得到的系数矩阵稀疏.在质量集中法推导系数矩阵元素的前提下,形成的系数矩阵是M-矩阵,从而保证格式单调,温度和能量不出负.双二次元破坏了低阶元系数矩阵结构的简单性,使M-矩阵性质难以满足.双二次元在矩形单元的顶点处构成了21点格式,在单元的边中点处是15点格式,在单元中心点处是9点格式.从而双二次元格式的系数矩阵非零元增多,带宽加大.用双二次元有限体积法求解二温方程组时,用质量集中法导出的格式在求解过程中尚未出现温度或能量出负的现象.3.解非平衡辐射扩散方程组的自适应方法<1>空间自适应以非平衡辐射扩散方程组(1)-(2)为模型,对于空间离散采用三角形网上的线性元有限体积法.令V=H1((2)×H1((2).相应于V的对偶空间记为V*.设方程组(1)-(2)的弱解存在.则有u=(E,T)∈V,满足以下变分形式：在以上式子中,(·,·)是定义在空间V上的一种内积,(?)tu=((?)tE,(?)tT),F(t,u)是从[0,T]×V到V*的一个算子,具体有关于自适应算法的后验误差指示子,先引入局部残量RK和单元边界跳量Jτ,这里的RK=(RKE,RKT)和Jτ=(JτE,JτT)分别定义在单元K和单元边τ上.令uh=(Eh,Th)是用三角网上线性元有限体积法得到方程组(1)-(2)的数值解.则有单元K上的残量：单元边τ上的跳量：其中K,K1,K2∈τh,τ是三角单元K1和K2的公共边,ni(i=1,2)是τ上分别关于单元Ki的单位外法向量.定义整体残量命题2.假设μ,υ是两个正的常数,且仅相关于三角剖分Γh中单元的最小内角,则对于t∈[0,T],以下不等式成立其中KH是K的直径,Hτ是τ的长度.命题2中的ηk：=(ηKE,ηKT)称为后验误差指示子.Rh受误差指示子ηK的控制,因此,我们可以依靠ηK来加细或粗化网格.自适应算法具体叙述如下.先根据初边值条件,对空间求解区域Ω构造初始时刻的自适应网格剖分Γh0.假设己有时间tn层的自适应网格Γhn,接下来需要构造时间tn-1层的自适应网格Γhn+Step1:在网格Γhn上用线性元有限体积法求解非平衡辐射扩散方程组,得到数值解；Step2:利用数值解计算Γhn上每个单元的后验误差指示子ηKα(α=E,T)；Step3:设定一个阀值Toln,找出ηKα较大的三角单元标记,构成集合R,同时找出ηKα较小的三角单元标记,组成集合C；Step4:将集合R中的所有单元加密,将集合C中的所有单元粗化,由此得到了新的网Γhn+1.令n=n+1,再从Step2开始重复以上步骤.试验结果表明自适应网格可以正确追踪数值解的突变区,在保证解满足一定精确性的同时,节省了CPU计算时间.<2>时间自适应由于能量E与温度T的变化趋势一致,且E的变化更加明显,因而我们设计的时间步长控制算法基于能量E的相对变化,即△E/E.在迭代收敛的情况下,当迭代数大于一定界限时,取△tn+1=min(△tren|1,0.8Atn),而当迭代数小于一定界限时,取△tn+1=min(△tren+1,1.1△tn,△tmax).4.解辐射三温能量方程组的低阶元有限体法与自适应考虑二维三温能量方程组：其中,Te、Ti、Tr分别为电子、离子、光子的温度函数.ρ为材料的密度,在每种材料区域内部都是常数,在不同材料的交界处间断.ce、ci、crTr3是比热.Ke、Ki、Kr分别为电子、离子、光子的扩散系数,具体有Kα=AαTα5/2(α=e,i),Kr=ArTr3+β,而Kα▽Tα(α=e,i,r)为热流密度.电子和离子的能量交换系数ωei=ρAeiTe-2/3,电子和光子的能量交换系数ωer=ρAerTe-1/2.参数cα,Aα(α=e,i,r),β,Aei,Aer仅依赖于介质材料,在每个子区域上是常量,在不同子区域的交界面处间断.方程的求解区域为Ω={(x,y)∈Ωxy,0≤t≤tmax},其中的Ωxy=∪i=13Ωi是由三种不同介质区组成的一个互不重叠的连通计算区域.Ω1是以直角坐标系的原点为圆心、以90μm为半径的半圆；Ω2是内半径和外半径分别为90μm和95μm的半圆环；Ω3是内半径和外半径分别为95μm和132μm的半圆环.初值条件：Tα(x,y,0)=Tα0(x,y)=3.0×10,α=e,i,r.边值条件：·固壁上的轴对称条件：Kα▽Tα·n|Γ2=0,α=e,i,r,其中n为单位外法向量；·自由面边界条件：Kα▽Tα·n|Γ1=0,α=e,i,Tr=Tr(x,y,t)|Γ1=2.0.<1>低阶元有限体积法对于时间离散,使用向后Euler法,并用冻结系数法线性化.对于空间离散,采用三角形网上的线性元有限体积法,并利用质量集中的方法处理重积分,得到全离散的线性元有限体格式.<2>自适应算法我们设计了一种新的基于后验误差指示子的h型自适应方法.先引入局部残量RK和单元边界跳量Jτ,这里的RK=(RKe,RKi,RKr)和Jτ=(Jτe,Jτi,Jτr)分别定义在单元K和单元边τ上.令uh=(Te,h,Ti,h,Tr,h)是用线性元有限体积法求得(14)-(16)式的数值解.则有单元K上的残量：和单元边τ上的跳量：其中K,K1,K2∈Γh,τ是三角单元K1和K2的公共边,ni(i=1,2)是τ上分别关于单元Ki的单位外法向量.后验误差指示子为ηK：=(ηKe,ηKi,ηKr),其中μ,y是两个正的常数,且仅相关于三角剖分Γh中单元的最小内角,HK是K的直径,Hτ是τ的长度.我们用后验误差指示子ηK作为自适应计算中调整网格的依据,并具体给出了自适应算法.数值试验结果表明当不用网格自适应时,随着网格的加细,能量守恒误差变得越来越小.并且在网格一致加密时,能量守恒误差的收敛速度几乎是一阶的.使用网格自适应算法时,网格加密区与光子温度突变区的移动趋势相一致.自适应算法有效地降低了能量守恒误差,并缩短了计算时间.