# Barycentric Rational Collocation Method for the Incompressible Forchheimer Flow in Porous Media

Authors: Qingli Zhao; Yongling Cheng
Journal: Journal of Mathematics
 (2021)

---

## Abstract

Barycentric rational collocation method is introduced to solve the Forchheimer law modeling incompressible fluids in porous media. The unknown velocity and pressure are approximated by the barycentric rational function. The main advantages of this method are high precision and efficiency. At the same time, the algorithm and program can be expanded to other problems. The numerical stability can be guaranteed. The matrix form of the collocation method is obtained from the discrete numerical schemes. Numerical analysis and error estimates for velocity and pressure are established. Numerical experiments are carried out to validate the convergence rates and show the efficiency.


---

## Body

## 1. Introduction

Darcy flow in porous media is of great interest in many science and engineering fields such as oil recovery and groundwater pollution contamination. Darcy’s law,(1)μkux=−dp^xdx+ρgx,x∈a,b,mainly describes the linear relationship between Darcy velocity u and derivative of pressure p^. Here, symbols μ, k, ρ, and gx represent the viscosity coefficient, permeability, the density of the fluids, and the gravitational term, respectively. This model is widely used and suitable for low velocity, small porosity, and permeability fluids [1–4].If the porosity is nonuniform and velocity is higher, a second-order term is needed to be added, the non-Darcy relationship has been researched by Forchheimer [1]. For example, the high-speed Forchheimer flow of single-phase incompressible fluid in porous medium is presented as follows:(2)μkux+βρuxux=−dp^xdx+ρgx,x∈a,b.Note that when Forchheimer numberβ=0, nonlinear model (2) degenerates to linear Darcy’s law (1).Model (2) is also called Darcy–Forchheimer law [5–10]. In [6], a block-centered finite difference method has been introduced to solve the Darcy–Forchheimer law. Discrete numerical scheme and error estimates were given. Mixed finite element method (MFEM) for equation (2) was studied in [7, 8]. Using this method, velocity and pressure can be approximated simultaneously. Two-grid and multigrid block-centered finite difference method (FDM) for the Darcy–Forchheimer flow in porous media was researched in [10, 11], respectively. This method can improve the efficiency of dealing with nonlinear problems. The barycentric formula is obtained by the Lagrange interpolation formula [12–16] and has been used to solve Volterra equation and Volterra integro-differential equation [12, 17, 18]. Floater and Hormann [19] have proposed a rational interpolation scheme which has higher accuracy on equidistant and special distributed nodes. Wang et al. [20–22] successfully applied the barycentric rational collocation method (BRCM) to solve initial value problem, boundary value problem, plane elasticity problem, and some nonlinear problems. These research studies extended the application fields of barycentric rational collocation method. In recent papers, Li et al. [23–27] have used the barycentric rational collocation method to solve heat conduction equation, biharmonic problem, and second-order Volterra integro-differential equation.In this paper, barycentric rational collocation method is introduced to solve the incompressible Forchheimer flow. We demonstrate that barycentric rational collocation method is highly accurate for both velocity and pressure.Ohd error estimates for velocity and pressure are given. Numerical experiments [28–32] are carried out to show the convergence rates. The paper is organized as follows. In Section 2, notations and barycentric formula are given. In Section 3, convergence analysis of barycentric rational collocation method for Forchheimer law and error estimates of velocity and pressure are presented. In Section 4, numerical examples are carried out to verify the convergence rates and show the efficiency. Throughout this paper, C denotes a positive constant independent of h.


## 2. Notations and Barycentric Rational Algorithm

The partition of intervalΩ=a,b is as follows:(3)a=x0<x1<⋯<xn−1<xn=b.Define(4)Ωi=xi−1,xi,1≤i≤n,hi=xi−xi−1,1≤i≤n,h=max1≤i≤nhi.For the functionux, the interpolation function rx (d=0,1,…,n) is given as(5)rx=∑i=0n−dλixpix∑i=0n−dλix.Symbolpix denotes the d-order interpolation polynomial such that pixk=uxk for k=i,i+1,…,i+d,(6)pix=∑k=ii+d∏j=i,j≠ki+dx−xjxk−xjuk,where uk=uxk and λix is a blending function(7)λix=−1ix−xi…x−xi+d.For the numerator term in (5), we deduce that(8)∑i=0n−dλixpix=∑i=0n−d−1i∑k=ii+d∏i,j≠ki+d1xk−xjuk=∑k=0nωkx−xkuk.Here,(9)ωk=∑i∈Jk∏j=i,j≠ki+d1xk−xj,and Jk=i∈I;k−d≤i≤k,I=0,1,…,n−d.Note that(10)∑k=ii+d∏j=i,j≠ki+dx−xjxk−xj=1,and for the denominator term in (5),(11)∑i=0n−dλix=∑k=0nωkx−xk.Through further deduction, we get(12)rx=∑j=0nωj/x−xjuj∑j=0nωj/x−xj=∑j=0nrjxuj.Here,ωk is described as (9). The basis function rjx of barycentric rational interpolation is(13)rjx=ωj/x−xj∑j=0nωj/x−xj,j=0,1,…,n.Then, we get the derivative formula at nodexi as(14)umxi≔uim=dmuxidxm=∑j=0nrjmxiuj=∑j=0nDijmuj,m=0,1,2,….Its matrix formulation can be given as(15)um=Dmu,where(16)um=u0m,u1m,…,unm,Dijm=rjmxi.The derivative formulation of the basis functionrjx at node xi is(17)rj′xi=ωj/ωixi−xj,j≠i,(18)ri′xi=−∑j≠irj′xi.According to induction (14)–(18), we obtain the recurrence formula of Dijm as(19)Dijm=mDiim−1Dij1−Dijm−1xi−xj,i≠j,Diim=∑j≠iDijm.


## 3. Convergence Rates and Error Estimates

Define the error betweenux and barycentric rational interpolation function rx as(20)ex=ux−rx.According to rational interpolation error theory, we know(21)ex=x−xi…x−xi+duxi,xi+1,…,xi+d,x.Combing (21) with (5), we see(22)ex=∑i=0n−dλixux−pix∑i=0n−dλix=AxBx,where(23)Ax=∑i=0n−d−1iuxi,xi+1,…,xi+d,x,Bx=∑i=0n−dλix.Define the error norm of e(x)as(24)ex=maxa≤x≤bex.The following lemma has been proved in [12].Lemma 1 (see [12]).
For the errorex defined as (20), we have(25)ex≤Chd+1,u∈Cd+2a,b,exx≤Chd,u∈Cd+3a,b.Now, we deal with the barycentric rational collocation schemes for the following Forchheimer equations:(26)Dp^≔dp^xdx=−μkux−βρuxux+ρgx,x∈a,b,Du≔duxdx=fx,x∈a,b,u0=u0,p^0=p^0.For the second equation of (26), the approximate formula is(27)∑j=0nrj′xuj=fx.Takingx=xi in (27), the numerical scheme is(28)∑j=0nrj′xiuj=fxi,i=0,1,…,n.For the first equation of (26), the approximate formula is as follows:(29)∑j=0nrj′xp^j=−μkux−βρuxux+ρgx.Then, the calculation scheme is(30)∑j=0nrj′xip^j=−μkui−βρuiui+ρgxi,i=0,1,…,n.Note that, in practical calculation, first step, we approximate the second equation of (26) and then the first equation of (26).Letuxn denote the numerical solution of ux, then we have(31)Duxn=fx,limn⟶∞Duxn=fx.Based on the above states, the next theorem gives the error analysis of Darcy velocity.Theorem 1.
Letuxn:Duxn=fx and fx∈Ca,b. If ux∈Cd+3a,b, then we have(32)ux−uxn≤Chd.Proof.
For the second equation of (26), using the notation of differential matrix, the discrete form of the collocation method is(33)L1u=D1u=f,where(34)L1=D001…D0n1⋮…⋮Dn01…Dnn1,u=u0,u1,⋅,unT,f=f0,f1,⋅,fnT.Furthermore, we have(35)Dux−Duxn=uxx−uxxn=∑j=0n−d−1juxxj,xj+1,…,xj+d,x∑j=0n−dλjx=exx=Ohd.The proof of this theorem is completed.Letp^xn denote the numerical solution of p^x, then we have(36)Dp^xn=−μkux−βρuxux+ρgx,limm⟶∞Dp^xn=−μkux−βρuxux+ρgx.The following theorem presents the error analysis of pressurep^.Theorem 2.
Letp^xn:Dp^xn=−μ/kux−βρuxux+ρgx and fx∈Ca,b. If p^x∈Cd+3a,b and ux∈Cd+3a,b, then we have(37)p^x−p^xn≤Chd.Proof.
For the first equation of (26), the discrete numerical scheme is(38)L2p^=D1p^=b,where(39)L2=D001…D0n1⋮…⋮Dn01…Dnn1,p^=p^0,p^1,…,p^nT,b=b0,b1,…,bnT,bj=−μkuj−βρujuj+ρgxj,0≤j≤n.Furthermore, we see(40)Dp^x−Dp^xn=−μkux−uxn−βρuxux−uxnuxn=E1+E2.AsE1, note that μ and k are positive constants, we have(41)E1=−μkux−uxn≤Cux−uxn=Ohd.Similarly, forE2, according to the monotonicity of the nonlinear term, we know(42)E2=−βρ|ux|ux−uxnuxn≤Cux−uxn=Ohd.Combing (40)–(42), the proof is finished.Remark 1.
In the above proof of Theorem2, coefficients μ,k,β, and ρ are supposed to be positive constants. If they are functions that depend on variable x and bounded, the proof is similar.


## 4. Numerical Experiments

In this section, we carry out some numerical experiments using barycentric rational collocation method to solve the Forchheimer equations.Example 1.
Consider the following incompressible Forchheimer model withΩ=0,1:(43)ux+310uxux=−dp^xdx+gx,x∈Ω,duxdx=πcosπx,x∈Ω,u0=0,p^0=0.
The analysis solution is chosen to be(44)ux=2sinπx,p^x=x2−x13.
Gravitational termgx is determined according to the first equation of (43). Define absolute error and relative error as(45)e1u=ux−uxn,er1u=ux−uxnux,e2u=ux−uxn2,er2u=ux−uxn2ux2,e1p^=p^x−p^xn,er1p^=p^x−p^xnux,e2p^=p^x−p^xn2,er2p^=p^x−p^xn2p^x2.
Numerical results are listed in Tables1–4 . The corresponding approximate figures between analysis solution and numerical solution can be seen in Figures 1 and 2. We test the barycentric rational with the uniform nodes for the direct methods. Tables 3 and 4 show that the convergence rates of velocity and pressure are Ohd with d=1,2,3,4. The theoretical convergence rate Ohd is reflected.Table 1
Errors of velocityu with n=20 for Example 1.
de1uer1ue2uer2u12.1667e – 031.0834e – 036.0702e – 033.0351e – 0329.6414e – 054.8207e – 052.7904e – 041.3952e – 0439.6877e – 064.8438e – 063.1901e – 051.5951e – 0546.4868e – 073.2434e – 072.2867e – 061.1433e – 0651.0361e – 075.1806e – 083.7468e – 071.8734e – 0768.9755e – 094.4877e – 093.3227e – 081.6614e – 0871.3998e – 096.9989e – 105.2848e – 092.6424e – 0981.5167e – 107.5836e – 115.8167e – 102.9083e – 1092.0588e – 111.0294e – 117.9930e – 113.9965e – 11104.6402e – 122.3201e – 121.8986e – 119.4929e – 12111.6687e – 128.3433e – 136.4924e – 123.2462e – 12Table 2
Errors of pressurep^ with n=20 for Example 1.
de1p^er1p^e2p^er2p^11.5603e – 012.5935e – 016.0581e – 011.0070e + 0025.5500e – 029.2256e – 022.2664e – 013.7673e – 0132.0024e – 023.3284e – 028.3781e – 021.3926e – 0147.0633e – 031.1741e – 023.0018e – 024.9897e – 0252.3541e – 033.9130e – 031.0104e – 021.6796e – 0267.3023e – 041.2138e – 033.1559e – 035.2458e – 0372.0459e – 043.4008e – 048.8865e – 041.4772e – 0385.1300e – 058.5273e – 052.2367e – 043.7180e – 0491.0396e – 051.7281e – 054.5468e – 057.5580e – 05102.0073e – 063.3367e – 068.7996e – 061.4627e – 05111.3584e – 072.2580e – 075.9683e – 079.9209e – 07Table 3
Errors of velocityu and convergence rates for Example 1.
hd=1Rated=2Rated=3Rated=4Rate1/202.1667e – 03—9.6414e – 05—9.6877e – 06–6.4868e – 07—1/404.7993e – 042.179.7770e – 063.304.9514e – 074.291.2656e – 085.681/801.0884e – 042.141.0429e – 063.232.5616e – 084.272.8257e – 105.491/1602.5226e – 052.111.1568e – 073.171.3571e – 094.246.7115e – 125.401/3205.9548e – 062.081.3223e – 083.137.3897e – 114.201.7519e – 135.26Table 4
Errors of pressurep^ and convergence rates for Example 1.
hd=1Rated=2Rated=3Rated=4Rate1/201.5603e – 01—5.5500e – 02—2.0024e – 02—7.0633e – 03—1/405.3988e – 021.531.0784e – 022.362.1501e – 033.224.1075e – 044.101/801.8789e – 021.521.9925e – 032.442.0891e – 043.362.0755e – 054.311/1606.5786e – 031.513.5968e – 042.471.9340e – 053.439.7970e – 074.401/3202.3128e – 031.516.4217e – 052.491.7487e – 063.474.4727e – 084.45Figure 1
Solution of fluid velocityu for Example 1 (d=1,h=1/320).Figure 2
Solution of pressurep^ for Example 1 (d=1,h=1/320).Example 2.
Consider the following incompressible Forchheimer flow:(46)ux+15uxux=−dp^xdx+gx,x∈Ω=0,1,duxdx=exsin8πx+8πcos8πx,x∈Ω=0,1,u0=0,p^0=0.
The analytical solution is set to be(47)ux=exsin8πx,p^x=2x+cos7πxlnx+1.
Gravitational termgx is determined according to the analytical solution. Numerical results are listed in Tables 5 and 6. The corresponding approximate results between analysis solution and numerical solution can be seen in Figures 3 and 4. We test the barycentric rational with the uniform nodes for the direct methods. Tables 5 and 6 show that the convergence rates are Ohd with d=1,2,3,4.Table 5
Errors of velocityu and convergence rates for Example 2.
hd=1Rated=2Rated=3Rated=4Rate1/201.2945e + 00—4.7246e – 01—1.7784e+00—1.6443e + 00—1/404.1691e – 011.631.5098e – 011.651.6003e – 013.477.8148e – 037.711/801.4368e – 011.543.2078e – 022.231.2557e – 023.671.4539e – 032.431/1605.0438e – 021.516.0933e – 032.401.0066e – 033.649.8861e – 053.891/3201.7802e – 021.501.1117e – 032.458.3782e – 053.595.0523e – 064.291/6406.2931e – 031.501.9941e – 042.487.1647e – 063.552.3727e – 074.411/12802.2254e – 031.503.5495e – 052.496.2243e – 073.521.0781e – 084.46Table 6
Errors of pressurep^ and convergence rates for Example 2.
hd=1Rated=2Rated=3Rated=4Rate1/201.4816e + 00—8.3941e – 01—2.6370e + 00—2.1833e + 00—1/404.6193e – 011.682.3137e – 011.862.0731e – 013.671.0082e – 027.761/801.5886e – 011.544.7049e – 022.401.5931e – 023.702.3613e – 032.091/1605.5677e – 021.518.7653e – 032.421.2503e – 033.671.4833e – 043.991/3201.9639e – 021.501.5856e – 032.471.0288e – 043.607.4072e – 064.321/6406.9404e – 031.502.8326e – 042.488.7430e – 063.563.4459e – 074.431/12802.4540e – 031.505.0319e – 052.497.5707e – 073.531.5589e – 084.47Figure 3
Solution of velocityu for Example 2 (d=1,h=1/1280).Figure 4
Solution of pressurep^ for Example 2 (d=1,h=1/1280).Remark 2.
Numerical experiments using the barycentric rational collocation method (BRCM) for the Forchhimer equations show the consistency of the convergence rates with the theoretical analysis. The main advantages of this BRCM are high precision and efficiency. The algorithm and program can be expanded to similar initial value problem and boundary value problem. It can effectively avoid the oscillation of other interpolation collocation methods. The numerical stability is guaranteed. We demonstrate that the proposed numerical scheme isOhd accurate for both Darcy velocity and pressure. In practical simulation, if the parameter d is further increased, we can get more accurate results. In the future, we will research higher-dimensional Forchheimer law and compressible fluids Forchheimer problems.



---

*Source: 5514916-2021-07-30.xml*