# A Numerical Method for Compressible Model of Contamination from Nuclear Waste in Porous Media

Authors: Zhifeng Wang
Journal: Mathematical Problems in Engineering
 (2021)

---

## Abstract

This paper studies and analyzes a model describing the flow of contaminated brines through the porous media under severe thermal conditions caused by the radioactive contaminants. The problem is approximated based on combining the mixed finite element method with the modified method of characteristics. In order to solve the resulting algebraic nonlinear equations efficiently, a two-grid method is presented and discussed in this paper. This approach includes a small nonlinear system on a coarse grid with sizeH and a linear system on a fine grid with size h. It follows from error estimates that asymptotically optimal accuracy can be obtained as long as the mesh sizes satisfy H=Oh1/3.


---

## Body

## 1. Introduction

A compressible nuclear waste disposal contamination problem in porous media is presented by the following coupled systems of partial differential equations. The physical processes can be concreted to be a high-level waste disposal buried in a salt dome, and next the salt dissolves to generate a brine, radioactive elements decay to generate heat, and finally the radionuclides are transported by the flow.Fluid:(1)ϕ1∂p∂t+∇⋅u=−q+Rs′c^,(2)u=−kμc^∇p=−ac^−1∇p,where p and u are the fluid pressure and Darcy velocity, respectively, ϕ1=ϕcw and ϕ is the porosity. q=qx,t is a production term, Rs′c^=csϕKsfs/1+cs1−c^ is a salt dissolution term, kx is the permeability of the rock, and μc^, the viscosity of the fluid, is dependent upon c^, the concentration of the brine in the fluid.Brine:(3)ϕ∂c^∂t+u⋅∇c^−∇⋅Ec∇c^=gc^,where Ec is the diffusion tensor including the molecular diffusion and mechanical diffusion and Ec=D+DmI, D=Di,j=αTuδi,j+αL−αTuiuj/u, and gc^=−c^csϕKsfs/1+cs1−c^−qc+Rs. Here, Dm is molecular diffusion. uianduj are two direction cosines of Darcy velocity. I is an identity matrix.Heat:(4)d1p∂p∂t+d2∂T∂t+cpu⋅∇T−∇⋅E˜H∇T=Qu,T,c^,p,where T is the temperature of the fluid, d1p=−ϕcwU0+p/ρ, d2=ϕcp+1−ϕρRcpR, E˜H=Dcpw+Km′I, Km′=km/ρ0, and Qu,T,c^,p=−∇U0−cp∇T0⋅u+U0+cpT−T0+p/ρ−q+Rs′−qL−qH−qH.Radionuclide (componenti):(5)ϕKi∂ci∂t+u⋅∇ci−∇⋅Ec∇ci+d3ci∂p∂t=fic^,c1,…,cN,i=1,…,N,where ci is the trace concentration of the i th radionuclide, d3ci=ϕcwciKi−1, and(6)fic^,c1,…,cN=ciq−csϕKsfs1+cs1−c^−qci−qci+q0i+∑j=1NkijλjKjϕcj−λiKiϕci.We assume the following:(1)
Zero Neumann boundary conditions for the equations(2)
The initial conditions are assumed given(3)
The mediumΩ is vertically homogeneous and take Ω∈R2(4)
The solutions are smooth andΩ periodic(5)
k,ϕ,ϕ1,d2,andKi are bounded below by positive constants, and Rs′c^,μc^,gc^,d3ci,fici,andQT are twice continuously differentiable with bounded partial derivatives about the variables in parentheses(6)
D=0Chou and Li [1], Ewing et al. [2], and Li et al. [3] have presented and studied several numerical methods for system (1)–(5) and its incompressible case. In this paper, we use the mixed finite element method to approximate the fluid problem and treat the brine, heat, and radionuclides by a modified method of characteristic finite element. It is well known that the full discrete approximation scheme is coupled and nonlinear. If simply lagging, the evaluation of the nonlinear items is used to obtain a linear system; it would be inevitable to introduce the constraint conditions about the mesh grid due to the stability requirement. Moreover, it would take an expensive cost to choose the implicit scheme to nonlinear solutions. An efficient method motivated by Xu [4] is considered in this paper. The method is used by Bi et al. [5], Chen et al. [6–10], and Liu et al. [11, 12] for solving some nonlinear problem. We shall relegate all of nonlinear iterations on a coarse grid of size H much coarser than the original fine grid of size h. According to the error estimates in the context, it obtains the asymptotically optimal accuracy to take H=Oh1/3.The remainder of the paper is organized as follows. Notations and approximation assumptions are given inSection 2. A two-grid method is defined and the convergence error estimates are derived in Section 3. In Section 4, we give some conclusions and extensions.


## 2. Notations and Approximation Results

To analyze the temporal discretization on a time interval0,T, let M be a positive integer number, △t=T/M, tn=n△t0≤n≤M, and φn=φ⋅,tn. Let LpJ;Wj,qΩ denote the usual set of functions with the norm(7)ψLpJ;Wj,qΩ=∫Jψ⋅,tWj,qΩpdt1/p,where if p=∞, the integral is replaced by the essential supremum. Let lpJ;Wj,qΩ denote the time discrete analogue to the space LpJ;Wj,qΩ with the norm(8)ψlpWj,qΩ=∑n=1NψnWj,qΩpΔt1/p.LetW=L2Ω and V=v∈Hdiv;Ω;v⋅γ=0. The weak form is presented as follows:(9)ϕ1∂p∂t,w+∇⋅u,w=−q+Rs′c^,w,∀w∈W,(10)ac^u,v−∇⋅v,p=0,∀v∈V,(11)ϕ∂c^∂t,z+u⋅∇c^,z+Ec∇c^,∇z=gc^,z,(12)d2∂T∂t,z+cpu⋅∇T,z+E˜H∇T,∇z+d1p∂p∂t,z=Qu,T,c^,p,z,(13)ϕKi∂ci∂t,z+u⋅∇ci,z+Ec∇ci,∇z+d3ci∂p∂t,z=fic^,c1,…,cN,z,for z∈H1Ω and i=1,…,N.Assume thatVh×Wh is the Raviart–Thomas space of index at least k associated with a quasitriangulation of Ω such that the elements have diameters bounded by hp. Let Mh=Mhc and Rh=RhT be finite-dimensional subspaces of W1,∞Ω for the approximation of concentrations and temperature, respectively, and we take Mh and Rh as the piecewise-polynomial space of degree at least l and r, respectively. As in [2, 11, 13], the approximation properties for Vh×Wh and Mh,Rh are given by(14)infvh∈Vhv−vh≤Cvk+1hpk+1,(15)infvh∈Vh∇⋅v−vh≤Cvk+1+∇⋅vk+1hpk+1,for v∈V∩Hk+1Ω2 and ∇⋅v∈Hk+1Ω, and(16)infwh∈Whw−wh≤Cwk+1hpk+1,w∈Hk+1Ω,(17)infzh∈Mhz−zh1,q≤Czl+1,qhcl,z∈Wl+1,qΩ,1≤q≤∞,(18)infzh∈Rhz−zh1,q≤Czr+1,qhTr,z∈Wr+1,qΩ,1≤q≤∞.If the initial solutionsph0,uh0,c^h0,Tˇh0,cih0∈Vh×Wh×Mh×Rh×MhN, the characteristics-Galerkin and mixed finite element approximation schemes are to find phn,uhn,c^hn,Tˇhn,cihn∈Vh×Wh×Mh×Rh×MhN satisfying(19)ϕ1phn−phn−1Δt,w+∇⋅uhn,w=−q+Rs′c^hn,w,∀w∈Wh,(20)ac^hnuhn,v−∇⋅v,phn=0,∀v∈Vh,(21)ϕc^hn−c^¯hn−1Δt,z+Ec∇c^hn,∇z=gc^hn,z,∀z∈Mh,(22)d2Tˇhn−Tˇ¯hn−1Δt,z+E˜h∇Tˇhn,∇z+d1phnphn−phn−1Δt,z=Quhn,Tˇhn,c^hn,phn,z,∀z∈Rh,(23)ϕKicihn−c¯ihn−1Δt,z+Ec∇cihn,∇z+d3cihnphn−phn−1Δt,z=fic^hn,c1hn,…,cNhn,z,∀z∈Mh,where i=1,…,N and(24)c^¯hn−1=c^hn−1x¯n−1,x¯n−1=x−uhnΔtϕ,c¯ihn−1=cihn−1x¯n−1,x¯n−1=x−uhnΔtϕKi,Tˇ¯hn−1=Tˇhn−1x¯n−1,x¯n−1=x−cpuhnΔtd2.Remark. If x¯n−1 is located outside Ω, we can join x¯n−1 with Y∈∂Ω so that x¯n−1−Y/x¯n−1−Y is the outer-normal direction to the boundary ∂Ω at Y. Take x∗∈Ω so that Y−x∗n−1=x¯n−1−Y, then we define x¯n−1=x∗n−1.In order to deduce the error estimates, we employ the elliptic projections by labeling them with tildes.(25)∇⋅u−U˜,w+p−P˜,w=0,w∈Wh,(26)ac^u−U˜,v+∇⋅v,p−P˜=0,v∈Vh,(27)Ec∇c^−C˜,∇z+λc^−C˜,z=0,z∈Mh,(28)EH∇T−T˜,∇z+λT−T˜,z=0,z∈Rh,(29)Ec∇ci−C˜i,∇z+λci−C˜i,z=0,z∈Mh,where U˜,P˜:J⟶Vh×Wh,C˜:J⟶Mh,T˜:J⟶Rh and C˜i:J⟶Mh for t∈J and introduce the following notations:(30)η=p−P˜,α=P˜−ph,ϱ=u−U˜,β=U˜−uh,ζ^=c^−C˜,χ^=C˜−c^h,ζi=ci−C˜i,χi=C˜i−cih,θ=T−T˜,ψ=T˜−Tˇh.Subtracting (19) from (9) and taking w=dtαn−1, we get the error equation about the pressure function as follows:(31)ϕ1dtαn−1,dtαn−1+∇⋅βn,dtαn−1=−ϕ1∂ηn∂t,dtαn−1+ϕ1dtP˜n−1−∂P˜n∂t,dtαn−1+∂Rs′∂c^δ1nc^n−c^hn,dtαn−1,where dtαn−1=αn−αn−1/Δt and δ1n is between c^n and c^hn.Next, combining (20) from (10) at t=tn with the test function βn,(32)ac^hnβn,βn−∇⋅βn,αn=−∂a∂c^δ2nU˜nc^n−c^hn,βn.Whent=tn−1, we apply the Taylor expansion and obtain that(33)ac^hn−1βn−1,βn−1−∇⋅βn,αn−1=ac^hn−1βn−1,βn−1−βn−∂a∂c^δ2n−1U˜n−1c^n−1−c^hn−1,βn,where δ2 is between c^ and c^h.Combining (32) with (33), we get(34)12dtac^hn−1βn−1,βn−1−∇⋅βn,dtαn−1=−12Δtac^hnβn,βn+12Δtac^hn−1βn,βn−12Δtac^hn−1βn,βn+1Δtac^hn−1βn−1,βn−12Δtac^hn−1βn−1,βn−1−dt∂a∂c^δ2n−1U˜n−1c^n−1−c^hn−1,βn.By (31) and (34),(35)ϕ1dtαn−1,dtαn−1+12dtac^hn−1βn−1,βn−1=−ϕ1∂ηn∂t,dtαn−1+ϕ1dtP˜n−1−∂P˜n∂t,dtαn−1+∂Rs′∂c^δ1nc^n−c^hn,dtαn−1−12Δtac^hn−ac^hn−1βn,βn−12Δtac^hn−1βn−βn−1,βn−βn−1−dt∂a∂c^δ2n−1U˜n−1c^n−1−c^hn−1,βn≤−ϕ1∂ηn∂t,dtαn−1+ϕ1dtP˜n−1−∂P˜n∂t,dtαn−1+∂Rs′∂c^δ1nc^n−c^hn,dtαn−1−12Δtac^hn−ac^hn−1βn,βn−dt∂a∂c^δ2n−1U˜n−1c^n−1−c^hn−1,βn.Using the deduction as [1, 2], we have(36)dtαn−12+12dtac^hn−1βn−1,βn−1≤βn∞2+εdtχ^n−12+Cχ^n−12+βn2+hp2k+2+hc2l+2+Δt2.After making the induction hypothesis thatsupβn∞⟶0, we multiply (36) by Δt and sum over 1≤n≤M to get(37)∑n=1Mdtαn−12Δt+ac^hM−1βM,βM≤βn∞2+ε∑n=1Mdtχ^n−12Δt+C∑n=1Mχ^n−112+βn2Δt+hp2k+2+hc2l+2+Δt2.Combining (11) and (21), we get the following equality, in which we choose the test function z=χ^n−χ^n−1=dtχ^n−1Δt and sum over 1≤n≤M:(38)∑n=1Mϕdtχ^n−1,dtχ^n−1Δt+12Ec∇χ^M,∇χ^M−12Ec∇χ^0,∇χ^0≤∑n=1Mϕ∂c^n∂t+uhn⋅∇c^n−ϕc^n−c^¯n−1Δt,dtχ^n−1Δt+∑n=1Mϕχ^¯n−1−χ^n−1Δt,dtχ^n−1Δt+∑n=1Mϕdtζ^n−1,dtχ^n−1Δt+∑n=1Mϕζ^n−1−ζ^¯n−1Δt,dtχ^n−1Δt+∑n=1Mλζ^n,dtχ^n−1Δt+∑n=1Mun−uhn⋅∇c^n,dtχ^n−1Δt+∑n=1M∂g∂c^c^hnζ^n+χ^n,dtχ^n−1Δt,where χ^¯n−1=χ^n−1x¯n−1,ζ^¯n−1=ζ^n−1x¯n−1, and x¯n−1=x−uhn/ϕxΔt.Note that(39)∑n=1Mϕζ^n−1−ζ^¯n−1Δt,dtχ^n−1Δt=ϕζ^M−1−ζ^¯M−1,χ^M−ϕζ^0−ζ^¯0,χ^0−∑n=1M−1ϕdtζ^n−ζ^¯n,χ^n=CΔtζ^M−12+ζ^02+∑n=1M−1ζ^n−12+dtζ^n2+εΔtχ^012.The reminder of the right side items in (38) is just as [2], that is,(40)∑n=1Mdtχ^n−12Δt+Ec∇χ^M,∇χ^M≤C∑n=1M∇χ^n−1∞2χ^n−112+βn2Δt+hp2k+2+hc2l+2+Δt2.It follows from the assumptionsup∇χ^n−1∞≤C and Gronwall lemma that(41)χ^L∞H1+dtχ^l2L2+βL∞L2+dtαl2L2≤CΔt+hcl+1+hpk+1.Then, by the inverse estimate and (40), we know that the induction hypotheses hold if(42)hpk+1=ohc,hcl+1=ohp,Δt=ohc,Δt=ohp.Finally, from the approximation properties,(43)c^−c^hL∞h1+dtc^−c^hl2L2+u−uhL∞L2+dtp−phl2L2≤cΔt+hcl+hpk+1.Similar to the above analysis, we can obtain the error estimates for the radionuclide equation and heat equation as follows:(44)∑i=1NχiL∞H1+∑i=1Ndtχil2L2≤CΔt+hcl+1+hpk+1,ψL∞H1+dtψl2L2≤CΔt+hcl+1+hTr+1+hpk+1,where hpk+1=ohT,hcl+1=ohT,andΔt=ohT are satisfied.Note that the time step is limited to beoh due to the theoretical proof.Theorem 1.
Definephn,uhn,c^hn,Tˇhn,cihn∈Vh×Wh×Mh×Rh×MhN for n≥1 by system (19)–(23) and assume that the approximation properties (25)–(29) hold. If(45)hpk+1=ohc,hcl+1=ohp,Δt=ohc,Δt=ohp,hpk+1=ohT,hcl+1=ohT,Δt=ohT,then there exists a positive constant C independent of h and △t, such that(46)c^−c^hL∞H1+dtc^−c^hl2L2+u−uhL∞L2+dtp−phl2L2+∑i=1Nci−cihL∞H1+∑i=1Ndtci−cihl2L2+T−TˇhL∞H1+dtT−T⌣hl2L2≤CΔt+hcl+hTr+hpk+1.
Similarly, we can get the error estimates of fine grid scheme inL2 norm.Theorem 2.
Definephn,uhn,c^hn,Tˇhn,cihn∈Vh×Wh×Mh×Rh×MhN for n≥1 by system (19)–(23) and assume that the approximation properties (25)–(29) hold. Then, there exists a positive constant C independent of h,H, and △t, such that(47)c^−c^hL∞L2+u−uhl2L2+p−phL∞L2+∑i=1Nci−cihL∞L2+T−TˇhL∞L2≤CΔt+hcl+1+hTr+1+hpk+1.


## 3. An Efficient Method

We now use and analyze a two-grid method for iteratively solving the nonlinear problem. The method has two steps.Stage 1. On the coarse grid TH with a mesh size H∈0,1, solve a small nonlinear system for PHn,UHn,C^Hn,THn,CiHn∈VH×WH×MH×RH×MHN given by (19)–(23).Stage 2. On the fine grid Th with a mesh size h∈0,1h≪H, solve the following linear system for Phn,Uhn,C^hn,Thn,Cihn∈Vh×Wh×Mh×Rh×MhN:(48)ϕ1Phn−Phn−1Δt,w+∇⋅Uhn,w−∂Rs′∂c^C^HnC^hn,w=−q+Rs′C^Hn−∂Rs′∂c^C^HnC^Hn,w,∀w∈Wh,(49)aC^HnUhn+∂a∂c^C^HnUHnC^hn,v−∇⋅v,Phn=∂a∂c^C^HnUHnC^Hn,v,∀v∈Vh,(50)ϕC^hn−C^¯hn−1Δt,z+Uhn⋅∇C^Hn,z+Ec∇C^hn,∇z−∂g∂c^C^HnC^hn,z=UHn⋅∇C^Hn,z+gC^Hn−∂g∂c^C^HnC^Hn,z,∀z∈Mh,(51)d2Thn−T¯hn−1Δt,z+E˜H∇Thn,∇z−∂QH∂TThn,z=−d1PhnPhn−Phn−1Δt,z+QH,z−∂QH∂TTHn,z,∀z∈Rh,(52)ϕKiCihn−C¯ihn−1Δt,z+Ec∇CiHn,∇z+∂d3∂ciCiHnPhn−Phn−1ΔtCihn,z−∂fiH∂c1C1hn+⋯+∂fiH∂cNCNhn,z=−d3CiHnPhn−Phn−1Δt,z+∂d3∂ciCiHnPhn−Phn−1ΔtCiHn+fiH−∂fiH∂c1C1Hn−⋯−∂fiH∂cNCNHn,z,∀z∈Mh,where(53)QH=QUhn,THn,C^hn,Phn,fiH=fiC^hn,C1Hn,…,CNHn,C^¯hn−1=C^hn−1x¯n−1,x¯n−1=x−UHnΔtϕ,C¯ihn−1=Cihn−1x¯n−1,x¯n−1=x−UhnΔtϕKi,T¯hn−1=Thn−1x¯n−1,x¯n−1=x−cpUhnΔtd2,and the projection on the fine grid is based on the numerical solutions on coarse grid.The sequential solution processes are defined as follows. Firstly, we apply the Newton iteration to the coupled system on the coarse grid and obtainPHn,UHn,C^Hn,THn,CiHn. Next, combining (48)–(50), we get PhnandUhn with RT1 and C^hn with piecewise linear finite element using a coupled linear system. Finally, from (51) and (52), we can get Cihn and Thn by parallel computation.In order to analyze the linear scheme on the fine grid, we define(54)π=P˜−Ph,σ=U˜−Uh,ξ^=C˜−C^h,ξ=C˜i−Cih,ω=T˜−Th.According to Taylor expansion, there exists a positiveδ3 such that(55)Rs′c^n−Rs′C^Hn−∂Rs′∂c^C^HnC^hn−C^Hn=Rs′C^Hn+∂Rs′∂c^C^Hnc^n−C^Hn−Rs′C^Hn−∂Rs′∂c^C^HnC^hn−C^Hn+∂2Rs′2∂c^2δ3c^n−C^Hn2=∂Rs′∂c^C^Hnc^n−C^hn+∂2Rs′2∂c^2δ3c^n−C^Hn2.According to (48), (9), and (25),(56)ϕ1dtπn−1,dtπn−1+∇⋅σn,dtπn−1=−ϕ1∂ηn∂t,dtπn−1+ϕ1dtP˜n−1−∂P˜n∂t,dtπn−1+ηn,dtπn−1+∂Rs′∂c^C^Hnc^n−C^hn+∂2Rs′2∂c^2δ1c^n−C^Hn2,dtπn−1.Like the deduction of (34), we see that(57)12dtaC^Hn−1σn−1,σn−1−∇⋅σn,dtπn−1=−12ΔtaC^Hnσn,σn+12ΔtaC^Hn−1σn,σn−12ΔtaC^Hn−1σn,σn+1ΔtaC^Hn−1σn−1,σn−12ΔtaC^Hn−1σn−1,σn−1−dt∂a∂c^C^Hn−1U˜n−1c^n−1−C^hn−1+∂a∂c^C^Hn−1U˜n−1−UHn−1C^hn−1−C^Hn−1,σn−12Δt∂2a∂c^2δ4U˜nc^n−C^Hn2−∂2a∂c^2δ4U˜n−1c^n−1−C^Hn−12,σn,where δ4 is between c^ and C^H.Hence, (56) and (57) can give that(58)ϕ1dtπn−1,dtπn−1+12dtaC^Hn−1σn−1,σn−1=−ϕ1∂ηn∂t,dtπn−1+ϕ1dtP˜n−1−∂P˜n∂t,dtπn−1+ηn,dtπn−1+∂Rs′∂c^C^Hnc^n−C^hn+∂2Rs′2∂c^2δ3c^n−C^Hn2,dtπn−1−12ΔtaC^Hn−aC^Hn−1σn,σn−12ΔtaC^Hn−1σn−σn−1,σn−σn−1,−dt∂a∂c^C^Hn−1U˜n−1c^n−1−C^hn−1+∂a∂c^C^Hn−1U˜n−1−UHn−1C^hn−1−C^Hn−1,σn−12Δt∂2a∂c^2δ4U˜nc^n−C^Hn2−∂2a∂c^2δ4U˜n−1c^n−1−C^Hn−12,σn.The error equation aboutc^ shows that(59)ϕξ^n−ξ^n−1Δt,z+Ec∇ξ^n,∇z=−ϕ∂c^n∂t+uHn⋅∇c^n−ϕc^n−c^¯n−1Δt,z+ϕξ^¯n−1−ξ^n−1Δt,z−ϕζ^n−ζ^¯n−1Δt,z+λζ^n,z−un−UHn⋅∇c^n−∇C^Hn+un−Uhn⋅∇C^Hn,z+∂g∂c^C^Hnζ^n+ξ^n,z+∂2g2∂c^2δ5c^n−C^Hn2,z,where ξ^¯n−1=ξ^n−1X¯n−1,ζ^¯n−1=ζ^n−1X¯n−1, and X¯n−1=x−UHn/ϕxΔt.Taking the test functionz=ξ^n−ξ^n−1=dtξ^n−1Δt and summing over 1≤n≤M, we have(60)∑n=1Mϕdtξ^n−1,dtξ^n−1Δt+12Ec∇ξ^M,∇ξ^M−12Ec∇ξ^0,∇ξ^0≤∑n=1Mϕ∂c^n∂t+uHn⋅∇c^n−ϕc^n−c^¯n−1Δt,dtξ^n−1Δt+∑n=1Mϕξ^¯n−1−ξ^n−1Δt,dtξ^n−1Δt+∑n=1Mϕζ^n−ζ^¯n−1Δt,dtξ^n−1Δt+∑n=1Mλζ^n,dtξ^n−1Δt+∑n=1Mun−UHn⋅∇c^n−∇C^Hn+un−Uhn⋅∇C^Hn,dtξ^n−1Δt+∑n=1M∂g∂c^C^Hnζ^n+ξ^n+∂2g2∂c^2δ4c^n−C^Hn2,dtξ^n−1Δt.Since(61)c^−C^H2≤c^−C^HL∞c^−C^HL2≤c^−C˜L∞+C˜−C^HL∞c^−C^HL2≤ClnHcHcl+1+lnHcHcl+1+Hpk+1Hcl+1+Hpk+1≤CHc2l+2lnHc+Hp2k+2lnHc,(62)u−UHc^−C^H≤u−UHc^−C^HL∞≤CHcl+1+Hpk+1lnHcHcl+1+Hpk+1≤CHc2l+2lnHc+Hp2k+2lnHc,(63)u−UH∇c^−∇C^H≤u−UHc^−C^H1,∞≤CHcl+1+Hpk+1Hcl+Hc−1Hcl+1+Hpk+1≤CHc2l+1+Hp2k+2Hc−1.It follows from (61)–(63) that(64)c^−C^hL∞H1+dtc^−C^hl2L2+u−UHL∞L2+dtp−Phl2L2≤CΔt+hcl+hpk+1+Hc2l+1+Hp2k+2Hc−1.At last, we present the error results of the radionuclide equation and heat equation as follows:(65)∑i=1Nci−CihL∞H1+∑i=1Ndtci−Cihl2L2≤CΔt+hcl+hpk+1+Hc2l+2lnHc+Hp2k+2lnHc,T−ThL∞H1+dtT−Thl2L2≤CΔt+hcl+hTr+hpk+1+HT2r+2lnHT+Hc2l+2lnHT+Hp2k+2lnHT.If settingl=1,r=1,k=1, hc=hT=hp, and Hc=HT=Hp, which accords with the assumption (45), we can get the following theorem of fine grid scheme.Theorem 3.
DefinePhn,Uhn,C^hn,Thn,Cihn∈Vh×Wh×Mh×Rh×MhN for n≥1 by system (48)–(52) and assume that the approximation properties (14)–(18) hold. Then, there exists a positive constant C independent of h,H, and △t, such that(66)c^−C^hL∞H1+dtc^−C^hl2L2+u−UhL∞L2+dtp−Phl2L2+∑i=1Nci−CihL∞H1+∑i=1Ndtci−Cihl2L2+T−ThL∞H1+dtT−Thl2L2≤CΔt+h+H3.Theorem 4.
DefinePhn,Uhn,C^hn,Thn,Cihn∈Vh×Wh×Mh×Rh×MhN for n≥1 by system (48)–(52) and assume that the approximation properties (14)–(18) hold. Then, there exists a positive constant C independent of h,H, and △t, such that(67)c^−C^hL∞L2+u−Uhl2L2+p−PhL∞L2+∑i=1Nci−CihL∞L2+T−ThL∞L2≤CΔt+h2+H3.


## 4. Conclusions and Extensions

The two-grid method presented in this paper reduces the complexity of problem. It involves a small nonlinear system on a coarse grid of sizeH and a linear system on a fine grid of size h. It is shown that the coarse space can be extremely coarse and still achieve asymptotically optimal approximation as long as the mesh sizes satisfy H=Oh1/3 in H1 norm. Compared with the implicit scheme, the two-grid method reduces CPU time. Moreover, the method is suitable to make the large-scale computation and long time duration. The future work is to use the discontinuous finite volume element method [14], block-centered finite difference method [15], higher-order finite volume [16], and SAV method [17, 18] to consider this problem.



---

*Source: 8872845-2021-05-19.xml*