Flow in Porous Media and Its Simulation

Petr Girg, University of West Bohemia, Czech Republic

ABSTRACT

The flow of fluid in porous medium consists of very complex 3D motions in geometricaly  complicated net of millions of micro-channels, which is difficult to handle using fundamental equations of fluid motion such as Navier-Stokes equations. Alternatively, several phenomenological approaches were proposed. In this talk, we focus on mathematical models based on quasilinear parabolic equations driven by the so called p-Laplacian (as the elliptic part of the equation) (Subscript[Δ, p])u=div(|grad u((|)^(p-2))grad u), p>1. 
These models cover situations when linear constitutive law (Darcy’s law) fails to hold due to prevailing effect of turbulence for high Reynolds’ numbers or due to interaction between walls of medium and liquid at molecular level for very low Reynolds’ numbers. The former case corresponds to 3/2<p<2 (singular diffusion) and the latter to 2<p<5 (degenerate diffusion), where the values of p are obtained experimentally. Both of these situations appear in important practical applications such as ground water modelling, petroleum and natural gas extraction. In the talk we will present several typical real-world situations, their formulations as initial-boundary value problems, and their numerical solutions. We use implicit semi-discretization in time which has much better stability then the explicit one, but requires to solve strongly nonlinear elliptic problems at each step. To solve the elliptic problem, we will use several approaches based on  finite differences, SparseArray and FindMinimum  and/or new features of NDSolve in v12 using finite elements for nonlinear elliptic problems. The use of the implicit semi-discretization in time is justified by our recent findings obtained jointly with Benedikt, Kotrla and Takáč.

PDE under consideration

Parabolic problems of the form

∂b(u)/∂t  -c div(|∇u((|)^(p-2))∇u)=f(x,t)        in Ω×(0,T)
                                                 u(x,t)=0                      on ∂Ω×(0,T)
                                                 u(x,0)=(Subscript[u, 0])(x)             in Ω

appear as mathematical models of motion of natural liquids in porous media such as underground water,  petroleum, and/or natural gas. Here porous media are e.g. clays, sands, gravels and/or various types of rocks.

Here

u: Ω×[0,T)-> ℝ, (x,t)↦u(x,t),

c>0, and  p>1.   b(u) is a continuous, increasing function. Typically, b(u)=(u^q), q>0, or similar functions.

Note that the diffusion coefficient    c|∇u((|)^(p-2)) gets 
a)   singular  for |∇u|≈0      and  1<p<2 
b)  degenerate for |∇u|≈0   and  2<p

Finite speed of propagation may happen see e.g. in:
[Barenblatt, G. I. On some unsteady motions of a liquid and gas in a porous medium. (Russian) Akad. Nauk SSSR. Prikl. Mat. Meh. 16, (1952). 67–78.]

Porous medium

Main difficulties of direct approach:

1) Complicated 3D net of channels.

2) Turbulence may occur.
3) Effects of interaction of natural fluids with rocks at mollecular level may occur.

Thus such problem is intractable using Navier-Stokes Equations.

Phenomenological approach

piezometric head at given location  h=p/(ρ g) + z

p hydrostatic preasure, ρ density of the fluid, z height from reference level

Governing equation : material balance (fluid) for elementary volume
change of mass in time = what flows in - what flows out + production of external sources

In differential form:

∂ρ/∂t = -div (ρ (OverVector[v])) + f(x,t)
 
 where ρ=ρ(x,t) is volumetric (or other) density and (OverVector[v])=(OverVector[v])(x,t) is filtration velocity.

Constitutive law:   relation between piezometric head H and filtration velocity
Darcy (France, 1850’s)                                             Δh/ ΔL = c  v
Smreker (Austria-Hungary, 1880’s)                 Δh/ ΔL = f(v) v,    Δh/ ΔL = c (v^(3/2))
Forchheimer (Austria-Hungary, 1890’s)       Δh/ ΔL = a v+b (v^2),   Δh/ ΔL = c (v^(3/2))
Izbash (Soviet Union, 1930’s)                                            v=(((Δh/ ΔL))^n),   1/2<n<1
Missbach (Czechoslovakia, 1930’s)                  Δh/ ΔL = c (v^m),                    1<m<2
and others ...

For incompressible fluid ρ is constant. 
Thus pluging vector valued Darcy law  (OverVector[v])=-c∇h into governing equation we get:

- c div (∇u) = f (x, t)

If turbulence occurs (Re=v d/ϑ>500; [Watanabe, Journal of Hydrology, 58 (1982) 389--397]), we must use nonlinear law, e.g. Missbach’s/Izbash’s formulas. Then we get

- c div (|∇u((|)^(p-2))∇u) = f (x, t)  where   p = 1 + 1/m. Thus  3/2<p<2.

Constitutive law :   Experiments

Experiments conducted by Daniel Veselý and Zdeněk Veselý, summer 2017

Data

oblazky=;
drt=;

S=Pi (125-2*3.2)^2/4/1000^2;

River8to16Free=Transpose[{
Reverse[{37.53,36.59,34.89,32.33,31.09,30.41,28.14,24.98,21.68,18.57,14.25,9.03}]/60000,
Reverse[{1.03,0.94,0.85,0.77,0.68,0.60,0.51,0.43,0.34,0.26,0.17,0.09}]
}];

Crushed8to16Free=Transpose[{
{60.10,57.83,57.10,49.44,48.16,47.15,43.87,39.95,35.10,29.35,23.69,15.29}/60000,
{1.03,0.94,0.85,0.77,0.68,0.60,0.51,0.43,0.34,0.26,0.17,0.09}
}];

dataFitPlot[measData_]:=Plot[Evaluate[{Fit[#, {q},q],Fit[#, {q,q^2},q], Exp[Fit[#/.{a_, b_}:>{Log[a], Log[b]}, {1,t},t]] /.t->Log[q]}&[measData]], {q,0,1.1Max[Transpose[measData][[1]]]},PlotRange->{{0,1.1Max[Transpose[measData][[1]]]},{0, 1.2}}, Epilog->Prepend[Point/@measData, PointSize[0.015]], PlotStyle->{{Black,Thick, DotDashed}, {Black,Thick, Dashed}, {Black, Thick}}, PlotLabels->{Text[Style["Darcy",{Bold, Medium}]],Text[Style["Forchheimer",{Bold, Medium}]] ,Text[Style["Power Law",{Bold, Medium}]]}, 
AxesLabel->{Text[Style["q [(m^3)/s]", {Black,Bold, Medium}]],Text[Style["Δh/ΔL", {Black,Bold, Medium}]]},TicksStyle->Directive[Black,Medium],PlotLegends->Placed["Expressions",{0.25,0.75}],ImageSize->800]

River stones 8 - 16 mm

Grid[{{Show[oblazky, ImageSize->Medium],dataFitPlot[River8to16Free]}}]

{Fit[#, {q},q],(*Fit[#, {q,q^2},q],*)Exp[Fit[#/.{a_, b_}:>{Log[a], Log[b]}, {1,t},t]] /.t->Log[q]}&[River8to16Free]

Solve[l==#, q]&/@%

Flow through unit area (1(m^2))

q/S/. %

Crushed stones 8 - 16 mm

Grid[{{Show[drt, ImageSize->Medium],dataFitPlot[Crushed8to16Free]}}]

q/S /. Solve[l==#, q]&/@({Fit[#, {q},q],(*Fit[#, {q,q^2},q],*)Exp[Fit[#/.{a_, b_}:>{Log[a], Log[b]}, {1,t},t]] /.t->Log[q]}&[Crushed8to16Free])

Example - Two parallel ditches

We are interested in the flow of the water in the gravel between two parallel ditches after fast discharge of the water at some place between the ditches.

Phenomenological approach leads to BVP for PDE

∂|u((|)^((p-1)/p))/∂t  - c div(|∂u/∂x((|)^(p-2))∂u/∂x)=0  in Ω×(0,T)
                                                                  u(x,t)=H              on ∂Ω×(0,T)
                                                                  u(x,0)=(Subscript[u, 0])(x)    in Ω

For gravel made of riverstones 8-16 mm: c=0.0591711, p=1+0.582042

Substitution u=w+H, we get  zero Dirichlet conditions:

∂b(w)/∂t  - c div(|∂w/∂x((|)^(p-2))∂w/∂x)=0          in Ω×(0,T)
                                                                     w(x,t)=0           on ∂Ω×(0,T)
                                                                     w(x,0)=(Subscript[u, 0])(x)-H    in Ω

Here b(w)=|w+H((|)^((p-1)/p))-|H((|)^((p-1)/p)).

(Notation (Subscript[Δ, p])w=div(|∂w/∂x((|)^(p-2))∂w/∂x))

Theoretical result
[J. Benedikt; P.G.; L. Kotrla; P. Takáč: The strong maximum principle in parabolic problems with the p-Laplacian in a domain. Appl. Math. Lett. 63 (2017), 95–101.]
suggest to use of implicit discretization in time of the parabolic problem in time

((b((w^(n+1)))-b((w^n)))/(Δt))=(Subscript[Δ, p])(w^(n+1))

This requires solving the following elliptic problem at each step

-(Subscript[Δ, p])(w^(n+1))+((b((w^(n+1))))/(Δt))=((b((w^n)))/(Δt))

Discretization of the operator by finite differences

(A=SparseArray[{Band[{1,1}]->1,Band[{2,1}]->-1},{401,400}])

(p-1)/p  /.  p->1.58204248554277 // N

b[x_]=(x+0.5)^0.3679057236842047-0.5^0.3679057236842047

Integrate[(x+0.5)^0.3679057236842047-0.5^0.3679057236842047,x]

B[x_]=-0.7749065674050282` x+0.7310445322991136` (((0.5` +x))^(1.3679057236842047`))

h=10/401

Solve[-0.7749065674050282`+(((0.5` +x))^0.367905723684204715377440936)==1.,x]

initCond=SparseArray[Table[{i}->4.2564270105249395`, {i,140,160}], 400]

ListPlot[b[initCond], Joined->True, PlotRange->All]

goal[x_, y_]:=0.3679057236842047`*0.059171130638716214`*Norm[Dot[A,x],1.5820424855427713`]^(1.5820424855427713`)*h*(1/h)^(1.5820424855427713`)/1.5820424855427713`-b[y].x * h*10+Norm[B[x],1]*h*10

oneStepInit[data_]:=(x /.FindMinimum[goal[x, data],{x, ConstantArray[1., 400]}, MaxIterations->3000][[2]]);

oneStepInit[initCond]

ListPlot[b[%], Joined->True, PlotRange->All]

Integration of PDE in Time

oneStep[data_]:=(x /.FindMinimum[goal[x, data],{x, data}, MaxIterations->3000][[2]])

NestList[oneStep,oneStepInit[initCond],300];

Two parallel ditches:  river stones 8 - 16 mm (linear versus nonlinear model)

Model of a heavy rain on an island

We use the following assumption, the island is made of gravel of pebblestones 8-16mm in diameter.

The island is circular of the diameter 1000 m.

The depth of the sea around is 50 m.

We assume heavy rain 25 cm per day for 12 h and assume (our least realistic assumption) that the rain water leaks immediately to underground water.

Ω=Disk[{0,0},1000];

Set up the partial differential operator.

c=-0.05917113063871185`*Sqrt[(10^(-50)+Grad[u[x,y],{x,y}].Grad[u[x,y],{x,y}])]^(0.5820424855-1);
op=(InactiveDiv[({{c,0},{0,c}}.(InactiveGrad[u[x,y], {x,y}, Inactive])), {x,y}, Inactive]);

b[x_]=(x+50)^0.367905723684204715377-50^0.367905723684204715377

Heavy rain: 25 cm / day

rain=0.25/60/60/24

Implicit discretization of parabolic problem

Discretization of the parabolic problem in time

((b((u^(n+1)))-b((u^n)))/(Δt))=(Subscript[Δ, p])(u^(n+1))+f(x,t)

requires solving the following elliptic problem at each step

-(Subscript[Δ, p])(u^(n+1))+((b((u^(n+1))))/(Δt))=((b((u^n)))/(Δt))+f(x,t)

$HistoryLength=3;

bag={};

usol=NDSolveValue[{op+b[u[x,y]]/1800==rain+0,DirichletCondition[u[x,y]==0,True]},u,{x,y}∈Ω,InitialSeeding->{u[x,y]==(1000^2-((x^2)+(y^2)))/1000^2},Method->{"FiniteElement","MeshOptions"->{MaxCellMeasure->100.}}];

bag=Append[bag, b[usol[x,y]]]

progress=1;

Dynamic[progress]

Do[Block[{},usol=NDSolveValue[{op+b[u[x,y]]/1800==rain+b[usol[x,y]]/1800,DirichletCondition[u[x,y]==0,True]},u,{x,y}∈Ω,InitialSeeding->{u[x,y]==(1000^2-((x^2)+(y^2)))/1000^2},Method->{"FiniteElement","MeshOptions"->{MaxCellMeasure->90.}}];progress++;bag=Append[bag, b[usol[x,y]]]],23]

Join[{Plot3D[0,{x,y}∈Ω, PlotRange->{0,0.01}]}, Plot3D[#,{x,y}∈Ω, PlotRange->{0,0.01}]&/@bag];

ListAnimate[%]

linbag={};

usol=NDSolveValue[{-0.06845901805278921*Laplacian[u[x,y],{x,y}]+b2[u[x,y]]/1800==rain+0,DirichletCondition[u[x,y]==0,True]},u,{x,y}∈Ω,InitialSeeding->{u[x,y]==(1000^2-((x^2)+(y^2)))/1000^2},Method->{"FiniteElement","MeshOptions"->{MaxCellMeasure->100.}}];

linbag=Append[linbag, b2[usol[x,y]]]

progress=1;

Dynamic[progress]

Do[Block[{},usol=NDSolveValue[{-0.06845901805278921*Laplacian[u[x,y],{x,y}]+b2[u[x,y]]/1800==rain+b2[usol[x,y]]/1800,DirichletCondition[u[x,y]==0,True]},u,{x,y}∈Ω,InitialSeeding->{u[x,y]==(1000^2-((x^2)+(y^2)))/1000^2},Method->{"FiniteElement","MeshOptions"->{MaxCellMeasure->90.}}];progress++;linbag=Append[linbag, b2[usol[x,y]]]],23]

ListAnimate[Join[{Plot3D[0,{x,y}∈Ω, PlotRange->{0,0.01}]}, Plot3D[#,{x,y}∈Ω, PlotRange->{0,0.3}]&/@linbag]]

linusol=NDSolveValue[{-0.06845901805278921*Laplacian[u[x,y],{x,y}]==rain,DirichletCondition[u[x,y]==0,True]},u,{x,y}∈Ω,InitialSeeding->{u[x,y]==(1000^2-((x^2)+(y^2)))/1000^2},Method->{"FiniteElement","MeshOptions"->{MaxCellMeasure->90.}}]

Plot3D[b2[linusol[x,y]],{x,y}∈Ω]

usol=NDSolveValue[{op==rain,DirichletCondition[u[x,y]==0,True]},u,{x,y}∈Ω,InitialSeeding->{u[x,y]==(1000^2-((x^2)+(y^2)))/1000^2},Method->{"FiniteElement","MeshOptions"->{MaxCellMeasure->100.}}];

Plot3D[b[usol[x,y]],{x,y}∈Ω]