渗流驱动问题的两层网格算法

陈艳萍;胡汉章

【摘 要】介绍了不可压缩的渗流驱动问题的特征有限元、特征混合有限元和特征扩张混合有限元两层网格算法;讨论可压缩的渗流驱动问题的特征有限元两层网格算法;通过数值例子验证了两层网格算法的有效性;讨论了渗流驱动问题的两层网格算法进一步的研究方向.

【期刊名称】《华南师范大学学报（自然科学版）》 【年(卷),期】2018(050)005 【总页数】8页(P98-105)

【关键词】非线性渗流耦合问题;特征有限元;混合有限元;两层网格算法 【作 者】陈艳萍;胡汉章

【作者单位】华南师范大学数学科学学院,广州510631;嘉应学院数学学院,梅州514015

【正文语种】中 文 【中图分类】O241.1

多孔介质中流体的运动是一个复杂而又普遍存在的物理现象,与实际生活当中的环境污染问题、海水入侵和油藏开采问题等密切相关,其数学模型相当复杂. 多孔介质中的流体运动通常由流体的输运(浓度)和流动(压力)两部分组成,宏观尺度下对其模拟是通过质量守恒方程及其Darcy 定律的变分形式实现的,CHEN和EWING[1]针

对该问题的数学模型和数值模拟做了全面的概述. 渗流驱动问题是这类流体问题中典型的一种,可分为：不可混溶、不可压缩油水两相驱动问题;可混溶、不可压缩油水两相驱动问题;可压缩相混溶的驱动问题. 其数学模型是由关于压力的流动方程(压力方程)和关于饱和度的对流扩散方程(浓度方程)耦合而成的非线性偏微分方程组.

过去的20余年,对渗流驱动问题数值方法的研究已有很多成果[2-10] . 对于压力方程,通常采用混合有限元逼近[4-5];对于浓度方程,主要离散方法有：特征有限差分[9]、特征有限元[7]、特征混合有限元[10]、间断有限元[8]和欧拉-拉格朗日局部伴随方法[11].

然而,对于非线性耦合渗流驱动问题,由上述离散方法导出的代数方程组是一个大系统的非线性耦合方程组. 两层网格算法作为求解非对称不定以及非线性问题的高效算法首次由XU[12-13]提出, 理论分析表明：两层网格算法在不降低求解精度的条件下,提高了计算效率. 几乎同时,CHEN和 HUANG[14]运用这种两层网格算法的思想,提出了求解非线性奇异两点边值问题的一种分层迭代校正方法,分析证明此方法可保证所有的高精度性质. 随着对这种高效的有限元两层网格算法研究的深入,国内外许多学者已经将其应用于各类不同的、具有实际应用背景的非线性或非线性耦合的偏微分方程问题[15-28]. 由于渗流驱动问题的非线性和耦合性,近几年来,我们针对不可压缩和可压缩的渗流驱动问题提出了特征有限元两层网格算法[29]、特征混合有限元两层网格算法[30].

本文介绍了不可压缩的渗流驱动问题的特征有限元两层网格算法、特征混合有限元两层网格算法,并在此基础上提出了特征扩张混合有限元两层网格算法;讨论可压缩的渗流驱动问题的特征有限元两层网格算法;利用文献[29-31]的数值例子,通过比较4个两层网格算法与相应的牛顿迭代法的收敛精度和运行时间说明两层网格算法的有效性.

1 不可压缩渗流驱动问题

1.1 不可压缩渗流驱动问题的数学模型

多孔介质中的油水两相渗流不可压缩、可混溶驱动问题的数学模型为： (1)

其中,xΩ,tJ=[0,T],p(x,t)、u(x,t)、c(x,t)分别是混合流体的压力、混合流体的Darcy速度和混合流体中某种物质的浓度,c0(x)是初始浓度,φ(x)是岩石的孔隙

度,a(c)=a(x,c)=k(x)/μ(c),μ(c)和k(x)分别是混合流体的黏度和绝对渗透度,q(x,t)是产量项,是非线性函数[4],(x,t)是注入流体的饱和度;D(u)为扩散和弥散张量,由分子扩散和机械弥散两部分组成: D(u)=φ[dmI+|u|(dlE(u)+dtE⊥(u))],

其中,E(u)是投影矩阵,(E(u))ij=uiuj/|u|2,E⊥(u)=I-E(u)是相补矩阵,dm 是分子扩散系数,dl、dt分别是流动方向、垂直方向的弥散系数. 问题(1)满足的边界条件和初始条件为: u·n=0 (x∂Ω,tJ), (2)

D(u)c·n=0 (x∂Ω,tJ), (3)

c(x,0)=c0(x) (xΩ), (4)

其中，n表示指向∂Ω的单位外法向量. 为了表述清楚,首先给出本文将要用到的函数空间. 记散度空间：H(div;Ω)={v(L2(Ω))2,·vL2(Ω)},定义空间和相应的范数为： V=H(div;Ω)∩{v·n=0 on

∂Ω},W={wL2(Ω),(w,1)=0},‖v‖V=‖v‖H(div;Ω)=(‖v‖2+‖·

‖w‖W=‖w‖.

设hc、hp表示Ω上四边形或者三角形拟一致剖分,其对应的剖分步长分别为hc、hp,记h=(hc,hp). 对浓度方程,逼近空间Mh⊂H1(Ω)是阶为l(≥1)且步长为hc的标准有限元空间[29],或者是逼近空间Mh×Vh⊂L(Ω)×V是阶为l(≥1)且步长为hc的混合有限元空间[30]；对压力方程,选择k阶Rarviart-Thomas混合有限元空间Vh×Wh⊂V×W,有·Vh=Wh.

1980年,EWINGT和WHEELER[2]对渗流驱动问题用Galerkin有限元法逼近;1983年,DOUGLAS[3]对不可压缩的两相渗流驱动问题提出了有限差分格式. 对于压力方程,由于浓度方程中只显含Darcy速度u而不是压力p,且混合有限元能同时逼近p和u,与标准有限元相比,在计算u时提高一阶精度,从而大大提高整体计算精度. 因此,DOUGLAS等[4-5]对压力方程采用混合有限元逼近. 由于浓度方程中对流占优的作用使抛物方程具有双曲性,为了克服经典方法可能出现的数值解振荡和失真,EWING等[7]利用特征线法与有限元法结合对浓度方程采用修正的特征有限元逼近. 随后,学者们对浓度方程提出了特征有限差分[9]、欧拉-拉格朗日局部伴随方法[11]、特征混合有限元[10]和间断有限元[8].

由于渗流驱动问题的非线性和耦合性,上述离散方法导出的代数方程组是一个大系统的非线性耦合代数方程组,研究求解这种非线性方程组的高效率、高精度算法有着重要的意义.

1.2 特征有限元两层网格算法

对于非线性耦合渗流驱动方程组,利用特征有限元-混合元获得的代数方程组是一个大系统的非线性耦合代数方程组,因此,HU等[29]考虑用全离散的特征有限元-混合元格式的两层网格算法(算法1)求解该非线性渗流耦合方程组. 这种格式的主要要素是定义在Ω上拟一致三角形或四边形剖分的粗网格有限元空间MH×VH×WH⊂Mh×Vh×Wh(h算法1 特征有限元-混合元两层网格算法

步骤1:在粗网格H上求解,,MH×VH×WH,使得它们满足如下的非线性耦合系统:

其中，D=φdmI为分子扩散张量.

步骤2:在细网格h上求解,,Mh×Vh×Wh,使得它们满足如下解耦的线性系统:

步骤3:在细网格h上求解,,Mh×Vh×Wh,使得它们满足如下解耦的线性系统：

此算法的程序是先在粗网格H上求解原非线性耦合问题,然后以步骤1的结果作为初值在细网格h上求解步骤2的线性解耦问题,最后以步骤2的结果作为初值在更细网格h上解步骤3的线性解耦问题. 文献[29]的理论分析和数值实验表明:算法1能保持非线性迭代法的误差精度且提高了计算效率. 但由算法1可知该算法没有考虑浓度的通量函数且不能保持单元上的质量守恒. 1.3 特征混合有限元两层网格算法

由于浓度方程实际上是对流占优的扩散方程,特征方法可以保证在流体峰线前沿逼近的高稳定性,消除数值弥散现象,并可以得到较小的时间截断误差,扩散项采用混合有限元离散,可以同时逼近浓度函数及其通量函数,而且可以保持单元上的质量守恒,因此,对浓度方程采用特征混合元逼近[30];注意到浓度方程中只显含u,对压力方程用混合有限元离散[4-5]. 用特征混合元-混合元获得的代数方程组也是一个大系统的非线性耦合代数方程组,因此,CHEN和HU[30]给出基于一次牛顿迭代的两层网格算法(算法2)，以求解该非线性渗流耦合方程组. 这种格式的主要要素是定义在Ω上拟一致三角形或四边形剖分的粗网格有限元空间MH×VH×VH×WH⊂Mh×Vh×Vh×Wh(h步骤1:在粗网格H上求解,,,MH×VH×VH×WH,使得它们满足如下的非线性耦合系统:

步骤2:在细网格h上求解,,,Mh×Vh×Vh×Wh,使得它们满足如下解耦的线性系统：

与算法1相比,算法2的优点是求出了浓度函数及其通量函数,而且可以保持单元上的质量守恒. 文献[30]的理论分析和数值实验表明:算法2同样能保持非线性迭代法的误差精度且提高了计算效率. 但扩散项很小时,它的逆在标准混合有限元中是不适用的[17]，这是算法2需改进的地方. 1.4 特征扩张混合有限元两层网格算法

由于实际问题中浓度方程的扩散项非常小,对浓度方程用特征扩张混合有限元离散会更加合理,基于1.2节和1.3节的分析,我们提出了特征扩张混合有限元-混合有限元两层网格算法(算法3).

算法3 特征扩张混合元-混合元两层网格算法

步骤1:在粗网格H上求解,,,,MH×VH×VH×VH×WH,使得它们满足如下的非线性耦合系统: (5)

其中，D=φdmI为分子扩散张量.

步骤2:在细网格h上求解,,,,Mh×Vh×Vh×Vh×Wh,使得它们满足如下解耦的线性系统:

与算法2相比,算法3的工作量并没有增加多少,但更加合理. 由于在粗网格上的方程组(5)是非线性耦合的方程组,为使两层网格算法的解误差阶达到渐近最佳逼近,在

分析两层网格解,,,,与真解(cn,,zn,un,pn)的误差时,需要用到、的Lq误差估计. 因此,必须先得到、的Lq误差估计.

为了得到与的Lq误差估计,首先把真解与数值解的Lq误差估计分解为真解与投影解和投影解与数值解的Lq误差估计. 接着用对偶论证和混合元投影算子的性质获得真解与投影解的Lq误差估计. 至于投影解与数值解的Lq误差估计,主要运用投影算子性质、方程组耦合性和超收敛性获得.

在两层网格算法的粗网格上得到对c、u误差的Lq估计基础上,再利用浓度沿特征方向导数的误差估计、Gronwall不等式获得两层网格算法的先验误差估计. 类似于文献[30]的理论分析,可以推出如下结果: ≤ Δt), (6) ≤ Δt). (7)

由式(6)、(7)可知:当粗细网格尺度满足H=O(h1/2)时,算法3可以获得方程组解的误差阶的最优阶逼近. 算法3的具体研究和分析将另文发表. 2 可压缩渗流驱动问题

2.1 可压缩渗流驱动问题的数学模型

多孔介质中的油水二相渗流可压缩、可混溶驱动问题的数学模型为: (8)

其中,xΩ,tJ=[0,T],ci(i=1,2)表示混合流体第i个分量的饱和度,c=c1=1-c2,b(c)=φ(x)c1{z1-∑zjcj},zj是压缩常数因子第j个分

量,d(c)=φ(x)∑,a(c)=a(x,c)=k(x)/μ(c),c(x,t)为相对饱和度,p(x,t)为混合流体的压力,k(x)是岩石的绝对渗透率,μ(c)是混合流体的黏度,u(x,t)是混合流体的Darcy速度,r(c)、d(x)分别是重力系数和垂直坐标,D(u)为扩散和弥散张量.

正如上节所述的不可压缩、可混溶油水两相渗流驱动问题的数学模型,DOUGLAS和ROBERTS[6]对可压缩渗流驱动问题的数学模型提出有限元-混合元和混合元-混合元离散格式;袁益让对该数学模型提出特征有限元[32]、特征差分[33]和迎风差分[34]等格式.

2.2 特征有限元两层网格算法

相比于不可压缩、可混溶油水两相渗流驱动问题的数学模型,可压缩、可混溶油水两相渗流驱动问题的数学模型的非线性和耦合性更强. 由有限元-混合元、混合元-混合元、特征有限元-混合元、特征混合元-混合元导出的离散代数方程组都是更大的非线性耦合方程组,其求解更加复杂. 因此,ZENG等[31]给出了基于一次牛顿迭代的特征有限元-混合元两层网格算法. 算法4 特征有限元-混合元两层网格算法

步骤1:在粗网格H上求解,,MH×VH×WH,使得它们满足如下的非线性耦合系统：

其中，D=φdmI为分子扩散张量.

步骤2:在细网格h上求解,,Mh×Vh×Wh,使得它们满足如下解耦的线性系统:

其中,B1=′

文献[31]的理论分析和数值实验表明:由于可压缩渗流驱动问题的非线性和耦合性更强,利用两层网格算法求解明显提高了计算效率. 3 数值实验

为了说明两层网格算法在不可压缩和可压缩渗流驱动问题中数值模拟的效果,利用

文献[29-31]的数值例子,分别用牛顿迭代算法和两层网格算法进行求解,通过对比相同网格尺度下的收敛精度和运行时间说明两层网格算法的效果. 浓度c、压力p属于分片常数空间,混合流速u、浓度的梯度和浓度的通量z属于最低阶的Raviart-Thomas元空间.

例1[29-30] 考虑如下不可压缩可混溶的渗流驱动问题: (9)

其中,Ω=[0,1]×[0,1],t[0,T],f(c)选取满足真解c=sin2(πx)sin2(πy)e-t,D=0.1×I,α(c)=c+2,

为了简单起见,假定Hp = Hc = H,hc = hp = h. 其中，算法1及相应的牛顿迭代法中假定时间步长为=1.0e-5,总时长T=1.0e-3；算法2、算法3及相应的牛顿迭代法中假定时间步长为=1.0e-5,总时长T=2.0e-4.

在不可压缩可混溶的渗流驱动问题中,由表1～表6可知：在相同的细网格尺度下,两层网格算法(算法1～算法3)与相应的牛顿迭代法的误差几乎相同,但随细网格空间步长的变小,算法1～算法3的计算效率明显提高. 即两层网格算法中,粗网格的网格数比细网格的网格数小得多,这样可以将大规模的非线性计算问题转化成小规模的非线性问题和大规模的线性问题进行求解,从而提高计算效率且不会影响细网格上有限元方法解的精度.

当两层网格算法是三步法时,由表1、表2可知：当粗细网格尺度满足H=O(h1/4)时,算法1可以获得方程组解的误差阶的最优阶逼近;当两层网格算法是二步法时,由表4～表6可以看出：当粗细网格尺度满足H=O(h1/2)时,算法2、算法3同样可以获得方程组解的误差阶的最优阶逼近.

表1 牛顿迭代法(特征有限元-混合有限元)的误差和运行时间Table 1 The error and CPU time of Newton iterative method (characteristic finite element-

mixed finite element)Hh误差误差阶‖cn-cnh‖ ‖un-unh‖‖pn-pnh‖‖cn-cnh‖ ‖un-unh‖‖pn-pnh‖运行时间/s1/41/165.617e-32.197e-16.873e-2———6.221/41/3.756e-45.503e-21.725e-2211105.431/41/2562.346e-51.415e-24.454e-32111 756.43 注：,,是该牛顿迭代法的解.

表2 两层网格算法(算法1)的误差和运行时间Table 2 The error and CPU time of two-grid algorithm(Algorithm 1)Hh误差误差阶‖cn-cnh‖ ‖un-unh‖‖pn-pnh‖‖cn-cnh‖ ‖un-unh‖‖pn-pnh‖运行时间/s1/41/166.003e-32.197e-16.874e-2———5.931/41/ 3.754e-45.500e-21.724e-2 2 1160.271/41/2562.346e-51.415e-24.456e-3 211853.46

表3 牛顿迭代法(特征混合有限元-混合有限元)的误差和运行时间Table 3 The error and CPU time of Newton iterative method(characteristic mixed finite element-mixed finite element)Hh误差误差阶运行时间/s‖cn-cnh‖‖zn-znh‖‖un-unh‖‖pn-pnh‖1/21/4 1.084e-1 8.100e-1 8.730e-12.856e-1—0.211/41/16 2.840e-2 2.177e-1 2.179e-1 6.880e-211.60 1/8 1/ 7.500e-3 5.0e-2 5.450e-2 1.720e-2132.071/161/1284.300e-3 3.080e-2 2.730e-2 8.600e-31168.15

注：,,,是该牛顿迭代法的解.

表4 两层网格算法(算法2)的误差和运行时间Table 4 The error and CPU time of two-grid algorithm(Algorithm 2)Hh误差误差阶运行时间/s‖cn-Cnh‖‖zn-Znh‖‖un-Unh‖‖pn-Pnh‖1/21/4 1.084e-1 8.100e-1 8.730e-1 2.856e-1 —0.371/41/162.840e-22.177e-12.179e-16.880e-210.781/81/7.500e-35.0e-25.450e-21.720e-2111.721/161/1284.300e-33.080e-22.730e-28.600e-3162.53

表5 牛顿迭代法(特征扩张混合有限元-混合有限元)的误差和运行时间Table 5 The error and CPU time of Newton iterative method(characteristic expanded mixed finite element-mixed finite element)Hh误差误差阶运行时间/s‖cn-cnh‖‖zn-znh‖‖zn-znh‖‖un-unh‖‖pn-pnh‖1/2 1/4 1.084e-1 8.092e-1 8.095e-2 8.633e-1 2.726e-1 —0.121/4 1/16 2.836e-2 2.173e-1 2.174e-2 2.365e-1 6.927e-2 1 0.911/8 1/ 7.510e-3 5.638e-2 5.638e-3 8.757e-2 1.7e-2 1 16.121/16 1/128 4.334e-3 3.085e-2 3.086e-3 4.403e-2 8.746e-3 1 1.08 注：,,,,是该牛顿迭代法的解.

表6 两层网格算法(算法3)的误差和运行时间Table 6 The error and CPU time of two-grid algorithm(Algorithm 3)Hh误差误差阶运行时间/s‖cn-Cnh‖‖zn-Znh‖‖zn-Znh‖‖un-Unh‖‖pn-Pnh‖1/2 1/4 1.084e-18.092e-18.092e-28.632e-12.726e-1—0.191/4 1/16 2.836e-22.173e-12.173e-22.366e-16.928e-21 0.831/8 1/ 7.510e-35.638e-25.638e-38.758e-21.7e-21 14.451/16 1/1284.334e-33.085e-23.085e-34.401e-28.744e-31 84.44 例2[31] 考虑如下可压缩可混溶的渗流驱动问题: (10)

其中,Ω=[0,1]×[0,1],t[0,T], f(c)选取满足真解

c=sin2(πx)sin2(πy)t,D=0.01×I,α(c)=1,b(c)=c(1-c),d(c)=c,p=-c+t/4. 为了简单起见,假定时间步长为=1.0e-5,总时长T=1.0e-3,Hp=Hc=H,hc=hp=h. 在可压缩可混溶的渗流驱动问题中,由表7、表8可以看出2种算法有相同的误差,但两层网格算法的运行时间少得多. 进一步可以看到：当粗细网格尺度满足H=O(h1/2)时,算法4同样可以获得解的误差阶的最优阶逼近.

表7 牛顿迭代法(特征有限元-混合有限元)的误差和运行时间Table 7 The error and CPU time of Newton iterative method(characteristic finite element-mixed finite element)Hh误差误差阶‖cn-cnh‖ ‖un-unh‖‖pn-pnh‖‖cn-cnh‖ ‖un-unh‖‖pn-pnh‖运行时间/s1/2 1/4 1.206e-5 1.719e-4 2.283e-5 ———1.321/4 1/16 4.921e-7 4.497e-5 6.099e-6 2 1 1 8.84 1/8 1/ 2.814e-8 1.234e-5 1.594e-6 2 1 1 96. 1/16 1/256 7.805e-9 6.509e-6 7.411e-7 21 1 267.28

注：,,是该牛顿迭代法的解.

表8 两层网格算法(算法4)的误差和运行时间Table 8 The error and CPU time of two-grid algorithm(Algorithm 4)Hh误差误差阶‖cn-Cnh‖ ‖un-Unh‖‖pn-Pnh‖‖cn-Cnh‖ ‖un-Unh‖‖pn-Pnh‖运行时间/s1/21/41.206e-51.719e-42.283e-5———1.231/41/1.921e-74.497e-56.099e-62114.021/81/2.814e-81.234e-51.594e-621132.571/161/2567.805e-96.509e-67.411e-7211103.73 4 小结

正如不可压缩渗流驱动问题讨论, 对可压缩渗流驱动问题同样讨论相应的特征有限元、特征混合元、特征扩张混合元及其他离散方法的两层网格算法. 但对于高维问题, 目前通常用分数步法分别与特征有限元、特征差分、迎风差分相结合的离散方法[35]. 对于其他离散方法的两层网格算法和高维问题的两层网格算法的文献还很少见,类似于特征有限元、特征混合有限元和特征扩张混合有限元等离散方法的两层网格算法可以进一步推广到其他离散方法的两层网格算法,但理论分析和数值算法实现将会更加复杂和困难,这将是下一步的研究方向. 参考文献：

【相关文献】

[1] CHEN Z X,EWING R E. Mathmatical analysis for reservior models[J]. SIAM Journal on Mathematic al Analy-sis,1999,30(2):431-453.

[2] EWING R E,WHEELER M F. Galerkin methods for miscible displacement problems in porous media[J]. SIAM Journal on Numerical Analysis,1980,17:351-365.

[3] DOUGLAS J J. Finite difference methods for two-phase incompressible flow in porous media[J]. SIAM Journal on Numerical Analysis,1983,20:681-696.

[4] DOUGLAS J J,EWING R E,WHEELER M F. A time-discretization procedure for a mixed finite elememt app-roximation of miscible displacement in porous media[J]. RAIRO:Analyse Numérique,1983,17:249-265.

[5] DOUGLAS J J,EWING R E,WHEELER M F. The app-roximation of the pressure by a mixed method in the simu-lation of miscible displacement[J]. RAIRO:Analyse Numé-rique,1983,17:17-33.

[6] DOUGLAS J J,ROBERTS J E. Numerical methods for a model for compressible miscible displacement in porous media[J]. Mathematics of Computation,1983,41:441-459. [7] EWING R E,RUSSELL T F,WHEELER M F. Convergence analysis of an approximation of miscible displacement in porous media by mixed finite elements and a modified method of characteristics[J]. Computer Methods in Applied Mechanics and Engineering,1984,47:73-92.

[8] KOU J S,SUN S Y. Analysis of a combined mixed finite element and discontinuous galerkin method for incompre-ssible two-phase flow in porous media[J]. Mathematical Methods in the Applied Sciences,2014,37(7):962-982.

[9] YUAN Y R. Characteristic finite difference methods for positive semidefinite problem of two phase miscible flow in porous media[J]. Systems Science and Mathematical Sciences,1999,12(4):299-306.

[10]SUN T J,YUAN Y R. An approximation of incompressible miscible displacement in porous media by mixed finite ele-ment method and characteristic-mixed finite element method[J]. Journal of Computational and Applied Mathematics,2009,228:391-411. [11]WANG H,EWING R E,QIN G,et al. A family of Eule-rian-Lagrangian localized adjoint methods for multi-dimensional advection-reaction equations[J]. Journal of Computational Physics,1999,152:120-163.

[12]XU J C. A novel two-grid method for semilinear equations[J]. SIAM Journal on Scientific Computing,1994,15:231-237.

[13]XU J C. Two-grid discretization techniques for linear and non-linear PDEs[J]. SIAM Journal on Numerical Analysis,1996,33:1759-1777.

[14]CHEN Y P,HUANG Y Q. A multilevel iterate correction method for solving nonlinear singular problems[J]. Natural Science Journal of Xiangtan University,1994,16:23-26. [15]CAI M,MU M,XU J C. Numerical solution to a mixed navier-stokes/darcy model by the two-grid approach[J]. SIAM Journal on Numerieal Analysis,2009,47:3325-3338.

[16]CHEN L P,CHEN Y P. Two-grid method for nonlinear reaction-diffusion equations by mixed finite element me-thods[J]. Journal of Scientific Computing,2011,49(3):383-401. [17]CHEN Y P,HUANG Y Q,YU D H. A two-grid method for expanded mixed finite-element solution of semilinear reaction-diffusion equations[J]. International Journal for Numerical Methods in Engineering,2003,57:193-209.

[18]CHEN Y P,LIU H W,LIU S. Analysis of two-grid methods for reaction-diffusion equations by expanded mixed finite element methods[J]. International Journal for Numerical Methods in Engineering,2007,69:408-422.

[19]CHEN Y P,LUAN P,LU Z L. Analysis of two-grid me-thods for nonlinear parabolic equations by expanded mixed finite element methods[J]. Advances in Applied Mathematics and Mechanics,2009,1(6):830-844.

[20]CHIEN C S,HUANG H T,JENG B W,et al. Two-grid discretization schemes for nonlinear schrödinger equations[J]. Journal of Computational and Applied Mathematics,2008,214(2):549-571.

[21]CHIEN C S,JENG B W. A two-grid discretization scheme for semilinear elliptic eigenvalue problems[J]. SIAM Journal on Scientific Computing,2006,27(4):1287-1304. [22]DAWSON C N,WHEELER M F. Two-grid methods for mixed finite element approximations of nonlinear parabolic equations[J]. Contemporary Mathematics,1994,180:191-203.

[23]HE Y N. Two-level method based on finite element and crank-nicolson extrapolation for the time-dependent na-vier-stokes equations[J]. SIAM Journal on Numerical Analysis,2003,41(4):1263-1285.

[24]JIN J C,SHU S,XU J C. A two-grid discretization method for decoupling systems of partial differential equations[J]. Mathematics of Computation,2006,75:1617-1626. [25]MU M,XU J C. A two-grid method of a mixed stokes-darcy model for coupling fluid flow with porous media flow[J]. SIAM Journal on Numerical Analysis,2007,45:1801-1813. [26]WU L. Two-grid mixed finite-element methods for nonlinear schrödinger equations[J]. Numerical Methods Partial Differential Equations,2012,28:63-73.

[27]WU L,ALLEN M B. A two-grid method for mixed finite-element solution of reaction diffusion equations[J]. Numerical Methods Partial Differential Equations,1999,15:317-332. [28]XU J C,ZHOU A H. A two-grid discretization scheme for eigenvalue problems[J].

Mathematics of Computation,2001,70:17-25.

[29]HU H Z,CHEN Y P,ZHOU J. Two-grid method for miscible displacement problem by mixed finite element me-thods and finite element method of characteristics[J]. Computers and Mathematics with Applications,2016,72(11):2694-2715.

[30]CHEN Y P,HU H Z. Two-grid method for miscible displacement problem by mixed finite element methods and mixed finite element method of characteristics[J]. Communications in Computational Physics,2016,19(5):1503-1528.

[31]ZENG J Y,CHEN Y P,HU H Z. Two-grid method for compressible miscible displacement problem by CFEM-MFEM[J]. Journal of Computational and Applied Mathematics,2018,337:175-1.

[32]袁益让. 多孔介质中可压缩可混溶驱动问题的特征有限元方法[J]. 计算数学,1992(4):385-406. YUAN Y R. Time stepping along characteristics for the finite element approximation of compressible miscible displacement in porous media[J]. Mathematica Numerica Sinica,1992(4):385-406.

[33]袁益让. 多孔介质中可压缩可混溶驱动问题的特征有限差分方法[J]. 计算数学,1993(1):16-28. YUAN Y R. Finite difference methods for a compressible miscible displacement problem in porous media[J]. Mathematica Numerica Sinica,1993(1):16-28.

[34]袁益让. 可压缩二相驱动问题的迎风差分格式及其理论分析[J]. 应用数学学报,2002,25(3):484-496.

YUAN Y R. The upwind finite difference methods for compressible two-phase displacement problem[J]. Acta Mathematicae Applicatae Sinica,2002,25(3):484-496. [35]袁益让. 可压缩两相驱动问题的分数步长特征差分格式[J]. 中国科学:A辑,1998,28(10):3-902.