Finite Volume Element Methods For Radiation Diffusion Equations

Radiation diffusion equations are used widely in ICF, MCF, astrophysics, hyper-sonic, combustion, and many other fields. When radiation and material can not reachequilibrium, non-equilibrium radiation diffusion is used to simulate. If the material (in-cluding electron, ion) balances, the question is simulated by two-temperature equations(i.e. non-equilibrium radiation diffusion equations). If the material itself is unbalanced,the question is simulated by three-temperature equations (i.e. three-temperature radia-tion diffusion equations) or group equations. Because such equations have multiscalecoefficients, strong nonlinearity, highly coupling, mutation regions and other charac-teristics, the numerical solution is very difficult to get. Correct and efficient simulationof such equations is a meaningful problem, so many scholars have researched on suchan issue [1–11].The finite volume element methods (FVEMs) are also called the generalized d-ifference methods, the control volume methods, or the box schemes. They are usedwidely in sciences and engineering, such as electromagnetics, computational fluid me-chanics, petroleum reservoir simulations. The linear element finite volume method andquadratic element finite volume method for solving the standard linear second orderelliptic equation and the corresponding theoretical derivation were studied in [12–17].According to features of this kind of radiation diffusion problems, we have improvedthe original finite volume element methods.The solution of radiation diffusion equations changes very intensely in some localregions,sosolvingthiskindofequationsonuniformgridsisunreasonable,whichbringsextremely time consuming and sometimes brings a huge amount of computation evenfar beyond the ability of computer hardware. An effective way to solve this problemis to use adaptive methods. Paper [18] proposed an adaptive algorithm for nonlinearparabolic equations based on a posteriori error. We have improved this method andextended to simulate two-temperature and three-temperature equations. Meanwhile,we designed a time direction adaptive method based on relative change of energy. This paper is composed of the following four parts:1. Low order element finite volume method for nonequilibrium radiation diffusion equationsNonequilibrium radiation diffusion equations where E is the radiation energy, T is the material temperature, D is the radiation d-iffusion coefficient, k=0.01T5/2is the material conduction coefficient. The photon absorption cross section is modeled by σa=(?), where the atomic mass number z rep-resents the material dependence. The energy exchange is controlled by σa, so high values of z and low values of T will result in tighter coupling between two equations. Coefficients D, k and σa all provide a strong nonlinearity for nonequilibrium radiation diffusion equations.<1>Two low order element finite volume schemesFor the temporal discreization, we use the backward Euler method. We utilize low order element finite volume methods, including the linear element finite volume method on triangular meshes and the isoparametric bilinear element finite volume method on quadrilateral meshes, to get spatial discreizations.The trial function space Uh is chosen as the linear finite element space correspond-ing to triangular partition Th or the isoparametric bilinear finite element space corre-sponding to quadrilateral partition. The dual element K*p is on barycenter dual partition or averaging center dual partition. The test function space Vh is chosen as a piecewise constant space corresponding to dual partition Th*. Based on various kinds of numerical methods we construct two discrete schemes for nonequilibrium radiation diffusion e-quations. We utilize a simple numerical integral formula, that is the mass concentration method, to get approximation of surface integrals. Because of the discontinuous nonlin-ear energy exchange coefficient σa, energy exchange term is approximated by a linear weighted combination, thereby obtaining Scheme I. Surface integrals of Scheme II are calculated exactly, and the discretizations of line integrals are the same as Scheme I.<2>Monotonicity and repair techniques Let Un=[En1,Tn1,En2,Tn2,...,Ens,Tns), where n denotes at time tn,(?) is the number of nodes of Th. Thus, discrete schemes can be written as We use the Freezing coefficient method to linearize (3) to getProposition1. Suppose that U0≥0and the matrix AT(U(n)) is an M-matrix, then Un+1≥0holds for all n≥0.The above proposition demonstrates that if we can offer some conditions to make AT(U(n)) to be an M-matrix, the corresponding scheme is monotone. For Scheme Ⅰ, AT(U(n)) is an M-matrix, when the partition can satisfy certain conditions (for triangular meshes, the maximal interior angle of the primary partition is not greater than90°; for parallel quadrilateral meshes, the ratio of adjacent sides of every primary element is between (?)3/3and (?)3, and its minimal internal angle is not less than60°). For Scheme Ⅱ, the property of M-matrix is hard to be satisfied. AT(U(n)) is an M-matrix for any vector U(n)≥0only when the number of primary elements is more than108even on equilateral triangular meshes. The constraint for meshes is quite harsh. Indeed, it is hard to be realized in practice. Therefore, two repair techniques used to deal with Scheme Ⅱ.The basic idea of repairing techniques is:if energy E (or temperature T) on some element is less than the minimum (Emin or Tmin), the needed energy is supplied by other elements and require that the modified solution has the same total energy εtotal. Local repair technique repaired E and T element by element. Global repair technique repaired E and T of all elements at the same time.Results of numerical tests show two kinds of fully-discrete linear element finite volume schemes are both energy conservative and robust on distorted meshes. We know the condition about M-matrix is sufficient but unnecessary for the monotonicity. Thus, we prove AT(U(n)) arising from Scheme Ⅰ is an M-matrix under certain conditions, but numerical results indicate that the solutions of Scheme Ⅰ all preserve positivity for given meshes. The solutions of Scheme Ⅱ can not preserve positivity even on uniform rect-angular meshes, but two repair techniques are effective to overcome this disadvantage. 2. High order element finite volume method for nonequilibrium radiation diffu-sion equationsChoosing the nonequilibrium radiation diffusion equations (1)-(2) as the model, we use backward Euler method for time discretization, and linearize equations by using frozen coefficient method. For spatial discretization, we choose biquadratic element finite volume method on rectangular meshes. Using biquadratic element finite volume method on rectangular meshes for solving two-temperature equations, meshes can be sparse relatively; but using bilinear element finite volume method, meshes should be dense. Degrees of freedom of biquadratic element finite volume method on the rectan-gle meshes (64x64) is equivalent with degrees of freedom of bilinear element finite volume method on the rectangle meshes (128x128). Bilinear element finite volume method on rectangular meshes is nine-point scheme, and the coefficient matrix is sparse. And using the mass concentration method to derive the coefficient matrix, the matrix is M-matrix, thus ensuring scheme to be monotonous, thereby ensuring temperature and energy not to be negative. Biquadratic element destroys the structure simplicity of the matrix derived from low order element, and it is harder to satisfy properties of M-matrix. Biquadratic element finite volume method on vertexes of rectangular meshes is21point scheme, and it is15point scheme on edge midpoints of rectangular meshes. And it is9point scheme on barycenters of rectangular meshes, so nonzero coefficients of matrix increase as a whole. When we use biquadratic element finite volume method to solve two-temperature equations, negative temperature or energy did not yet appear in the whole solving process.3. Adaptive algorithms for nonequilibrium radiation diffusion equations<1>Spatial adaptivityConsidering the nonequilibrium radiation diffusion equations (1)-(2), we choose linear element finite volume method on triangular meshes for spatial discretization. Let V=H1(Ω) x H1(Ω). The dual space of V is denoted by V*. We assume that the weak solution of equations (1)-(2) exists. It means there exists u=(E, T) e V, satisfying the following abstract variational formulation: In the above equation,(·,·> is a kind of inner product defined on V, dtu=(dtE, dtT), and F(t, u) is the operator defined from [0, T] x V to V*, satisfyingIn order to describe the posteriori error estimates conveniently, we need to give the definitions of the local residual RK and the jump JT, here RK=(RE/K,RT/K) and JT=(JE/T,JT/T) are defined on the triangular element K and edge T, respectively. Let Uh=(Eh, Th) be the solution of equations (1)-(2), and then the local residual of element K is defined by And the jump of element edge T is where K, K1, K2∈Th, T is the common edge of K1and K2, ni is the unit outward normal vector relatively to Kt(i=1,2).Let the total residual beProposition2. Suppose that μ, v are two positive constants, only relevant to the s-mallest inner angle of the triangulation Th-Then it holds that for t∈[0,T] whereηk:=(ηE/k,ηT/k) in Proposition2is the posteriori error indicator. Rh can be con-trolled by the error indicator ηk. Hence, we can refine or coarsen meshes depending on ηK. We present the adaptive mesh strategy.Firstly, we need to set up the initial triangular partition of computation domain Q. The initial mesh To/h is got based on the initial boundary value conditions. Suppose the adaptive mesh Tn/h at time t=tn has been obtained, we need to generate the adaptive mesh Tn/h1+1at time tn+Step1:Using linear element finite volume method to solve the nonequilibrium radiation diffusion equations on the mesh Tn/h, we get the numerical solution;Step2:We use the numerical solution to calculate the posteriori error indicator ηa/K (a=E, T) of every element on Tn/h;Step3:Set a threshold Toln, and then find out triangular elements with larger η(?) to mark, which constitute the setR, while find out triangular elements with smaller η(?) to mark, which constitute the set C;Step4:Refine all elements of%, and coarsen all elements of C to get the new mesh (?). Set n=n+1, and then begin from Step2again.The final test results show that the moving mesh can track the mutation of nu-merical solution correctly, and can save CPU time while ensuring the correctness of solution.<2> Temporal adaptivityDue to the energy E in accordance with the change trend of temperature T, and the change of E is more obvious, so we design the time step length control algorith-m based on relative change of energy E, viz.△E/E. Under the condition of iter-ative convergence, when the iteration number exceeds a certain limit, let△tn+1min(△t(?)+1,0.8△tn). When the iteration number is smaller than a certain limit, let Atn+1=min(△t(?),1.1△tn,△tmax)4. Low order element finite volume method and adaptation for3-T energy equa-tionsTwo dimensional energy equations coupled with electron temperature Te, ion tem- perature Ti and photon temperature Tr are defined by where p denotes the density of the medium, which is a constant within each subdo-main and discontinuous across interfaces of subdomains. ce、ci、crT3r are heat ca-pacities, Ke、 K、 Kr are diffusion coefficients of electron, iron, photon, where Ka=AaTa5/2(a=e,i), Kr=ArTr3+β, and Ka(?)Ta(a=e,i,r) are heat flow densities. The energy exchange coefficient between electron and iron is (?)ei=pAeiT-2/3and the en-ergy exchange coefficient between electron and photon is (?)er=pAerTe-1/2. The above parameters ca,Aa(a=e, i, r), β, Aei, Aer are constant within each subdomain, but they are discontinuous across interfaces of subdomains.Computation domain is Ω=[(x,y)∈Ωxy,0≤t≤tmax}, where Ωxy=U3/i=Ωi is a non-overlapping connected region consisting of three different media. Ω1is a semicircle with coordinate origin as the center of the circle,90μm as the radius; Ω2is the half of an annulus, where inside radius is90μm and outside radius is95μm; Ω3is the half of an annulus, where inside radius is95μm and outside radius is132μm.Initial conditions:Tα(x,y,0)=Tα/o(x,y)=3.0x10-4, a=e, i, r. Boundary conditions:· Wall:Ka(?)Ta· n|r2=0, a=e, i, r, where n is the unit outward normal vector;· Free:Ka(?)Ta· nT1=0, a=e,i, Tr=Tr(x,y,t)T1=2.0.<1>Low order element finite volume methodWe use backward Euler method for time discretization, and linearize equations by using frozen coefficient method. For spatial discretization, we choose linear element finite volume method on triangular meshes.<2>Adaptive algorithmWe design a new h-type adaptive method based on a posteriori error indicator. In order to describe the posteriori error estimates conveniently, we need to give the definitions of the local residual Rk and the jumpJT, here Rk=(Rek, RiK, RrK) and JT (JeT, JiT, JrT) are defined on the triangular element K and edge T, respectively. Let uh (Te,h, Ti,h, Tr,h) be the solution of equations (14)-(16), and then the residual of element K is defined by and the jump of element edge T is where K, K1, K2∈Th, T is the common edge of K1and K2, ni is the unit outward normal vector relatively to Ki (i=1,2).The posteriori error estimator ηK:=(ηeK, ηiK,ηrK), where μ, v are two positive constants, only relevant to the smallest inner angle of the triangu-lation Th. We present an adaptive algorithm based on the posteriori error estimator ηK for3-T energy equations.Results of numerical tests show along with the grid refinement, the energy conser-vation error is becoming smaller and smaller. When the grid refines consistently (all elements are divided into four), the convergence rate of energy conservation error is al-most first order. When using grid adaptive algorithm, grid refinement area and photon, electron, ion temperature mutation areas are moving consistently. All the isotherms are symmetrical half circles, and the positions of mutation areas have continued to change. The adaptive algorithm can reduce the energy conservation error and shorten the calcu-lation time effectively.