The quasilinear system of a composite type describing the spatial unsteady isothermal motion of a compressible fluid in a viscoelastic porous medium is considered. In this formulation, the influence of gravity is taken into account, the full equation of the balance of forces is considered, and the viscoelastic properties of the deformable porous skeleton are also taken into account. The problem is reduced to a single thirdorder equation for finding the porosity function in self-similar variables. The system is reduced to a second-order differential equation in the case of the predominance of viscous properties. A numerical study of this case by the method of determination is carried out.
Copyright ⃝c by the paper's authors. Use permitted under Creative Commons License Attribution 4.0 International (CC BY 4.0). relation to the problems of seismology, non-destructive testing, acousto-electronics, and a number of other areas [Bio56]. 2
Problem statement We consider the following quasilinear system of composite type describing the spatial nonstationary isothermal motion of a compressible fluid in a viscoelastic porous medium [Mor07],[Fow11]: ∂(1 − ϕ)ρs + div((1 − ϕ)ρsv⃗s) = 0, ∂t ∂(ρf ϕ) ∂t
+ div(ρf ϕv⃗f ) = 0, kϕn ϕ(⃗vf − ⃗vs) = − µ (∇pf − ρf⃗g), ∇ · ⃗vs = −a1(ϕ)pe − a2(ϕ)( ∂∂pte + ⃗vs · ∇pe),
( ∇ptot = ρtot⃗g + div (1 − ϕ)η ( ∂⃗vs + ( ∂⃗vs )∗)) ∂⃗x ∂⃗x , ρtot = ϕρf + (1 − ϕ)ρs, ptot = ϕpf + (1 − ϕ)ps, pe = (1 − ϕ)(ps − pf ). ρs, ρf = const, a1 = ϕηm , a2 = βϕϕb [Con98]: Here ϕ is the porosity; ρf , ρs, ⃗vs, ⃗vf are the true densities and velocities of the phases, respectively; pe is the effective pressure, ptot is the total pressure ρtot is the total density; ⃗g is the density of mass forces;βϕ is the compressibility coefficient of the solid skeleton, η is the dynamic viscosity of the solid phase, k is the permeability, µ is the dynamic viscosity of the liquid, σ is the total stress tensor. The true density of the solid phase ρs is assumed to be constant. The system is closed if ρs, ρf = const. The solvability of the self-similar problem for the original system of equations is established in [Tok16] in the case of an incomplete forces balance equation ∇ptot = −ρtot⃗g and ⃗g = 0. The solvability of the initial-boundary value problem for the equations of nonisothermal filtration in the case of the prevalence of the viscous properties of the skeleton was established in [Pap21]. Global in time solvability in the case of isothermal filtration was proved in [Pap19].
We arrive to a closed system of equations for ϕ, vs, vf , ps, pf in the one-dimensional case, under the condition ∂ϕ ∂t
+ ∂(1 − ϕ) + ∂t ∂ ∂x ∂ ∂x (ϕvf ) = 0, (vs(1 − ϕ)) = 0, k n ∂pf + ρf g), ϕ(vf − vs) = − µ ϕ ( ∂x ∂∂vxs = − ϕηm pe − βϕϕb( ∂∂pte + vs ∂∂pxe ), ∂ptot = −ρtotg + ∂x ∂ ( ∂x
∂vs ).
2η(1 − ϕ) ∂x t = t1t′, x = x1x, vf = v1vf′ , vs = v1vs′, pf = p1p′f , ptot = p1p′tot, pe = p1p′e, ps = p1p′s,
We pass in this system to dimensionless variables: (hereinafter, the primes are omitted), and also put: t1 = xv11 , α = µkvp121t1 , β = kµρvf1g , γ = p1ηt1 , λ = βf p1, ζ = 2η x1t1ρf g , ρ = ρs , κ = ρf
p1 . gx1ρf Then equations can now be written in the form ∂ϕ
+ ∂t ∂(1 − ϕ) + ∂t ∂(vf ϕ)
= 0, ∂x ∂(vs(1 − ϕ)) = 0,
∂x ϕ(vf − vs) = −ϕn(α ∂∂pxf + β ,
) ∂vs = −ϕmpeγ − ϕ λ ∂pe ),
b ( ∂pe + vs ∂x ∂x ∂t κ ∂ptot = ζ ∂ ( ∂x ∂x
∂vs ) (1 − ϕ) ∂x − ϕ − (1 − ϕ)ρ.
Next, we consider a self-similar solution of the ”traveling wave” type. Assuming that all the required functions depend only on the variable ξ = x − ct(ξ > 0, c is a constant parameter). After some transformations, we arrive to the following system of equations: dϕvf dξ
dϕ − c dξ
= 0, d dξ ((1 − ϕ)vs) − c d(1 − ϕ) = 0,
dξ ϕ(vf − vs) = −αϕn dpf
dξ − βϕn, ddvξs = −γϕmpe + λcϕb ddpξe − λϕbvs ddpξe , κ dptot = ζ d ( dξ dξ
dvs ) (1 − ϕ) dξ
− ϕ − ρ(1 − ϕ). vs(0) = vs0, vf = vf0, ϕ(0) = ϕ0, lim ϕ(ξ) = ϕ+,
ξ→∞ lim vf (ξ) = u+, lim ϕ(ξ) = ϕ+, ξ→∞ ξ→∞ c = A2 = ϕ+(1 − ϕ0)vs0 − ϕ0(1 − ϕ+)vf0 ϕ+ − ϕ0
, u+ = vf0ϕ0 + (1 − ϕ0)vs0, (1 − ϕ+)ϕ0(1 − ϕ0)(vf0 − vs ) ,
0 ϕ+ − ϕ0
ϕ+
A1 = 1 − ϕ+ A2. ϕ( A1 A2 ) = −ϕn(α dpf + β), ϕ − 1 − ϕ dξ The system is supplemented with boundary conditions: where vs0, vf0, ϕ0, ϕ+ are given constants satisfying the conditions
vs0 ̸= vf0, ϕ0 ̸= ϕ+.
From the equations (1) – (2) of the system, we obtain: Thus, the system is converted to the following form for finding functions ϕ, pf , ptot: (1) (2) (3) (4) (5) (6)
λA2ϕb dpe = −γϕm(ptot − pf ) − 1 − ϕ dξ , d ( κ dptot = ζA2 dξ dξ (1 − ϕ)
Next, we obtain the equation for the porosity function. We express from equations (6) and (8) the resulting system ddpξf and dpdtξot , respectively. Then we divide (7) by ϕm, differentiate it, and substitute the expressed derivatives. Thus, we obtain a third-order equation for finding the function ϕ: ( b −ϕm1+n− n + ζλA2 ϕb−m ( b − m κγ (1 − ϕ)3 ϕ + 3
1 ϕn−1(1 − ϕ)
1 1 − ϕ )( dϕ )3 dξ ) + + ( βκ − αρ ) ( b − m ακ ϕ
+ λζA2 ϕb−m ( b − m κγ (1 − ϕ)2 ϕ
In the case of the predominance of the viscous properties of the medium, only the first term in the right part will remain in the second equation of this system. Then, in the same way, we can obtain a second-order equation for finding the function ϕ: where +
1 ( ζ (1 − ϕ)2 κ + 2 1 m 1 γ ϕm(1 − ϕ) − γ ϕ1+m )( dϕ )2
dξ + A1 ϕ−n
A2 − 1κ−A2ρ ϕ − 1ϕ1−−nϕ +
β ρ αA2 − κA2
= 0.
d2ϕ A(ϕ) dξ2 + B(ϕ) dξ ( dϕ )2
+ C(ϕ) = 0, A(ϕ) = (
ϕ−m ζ γ(1 − ϕ)2 + κ(1 − ϕ) )
,
B(ϕ) = C(ϕ) = (
ζ 2ϕ−m m ϕ−m−1 ) κ(1 − ϕ)2 + γ(1 − ϕ)3 − γ (1 − ϕ)2 , A1 ϕ−n αA2
ϕ1−n − α(1 − ϕ) − ϕ + (1 − ϕ)ρ κA2 +
β αA2 . 3
The search for a solution to this equation is performed by the method of determination [Kha08]. The solution of a boundary value problem can be interpreted as an equilibrium state, which is approached by the solution of a non-stationary problem. Sometimes there is a situation when it is more convenient and more efficient, from a computational point of view, to solve such a unsteady problem than to directly search for a solution to the original boundary value problem. This problem can be solved by reducing semi-infinite interval [0; +∞) to the finite [0; ξ∗], where ξ∗ is found from the condition:
|ϕ(ξ∗) − ϕ+| ≤ ε, where ε – is the desired accuracy of the solution. The search for a solution is performed with the necessary accuracy using the condition |ϕin+1 − ϕin| ≤ ε and in the case under consideration, ε = 0.005.
In the region [0, ξ∗] × [0, 1] we construct a uniform grid ω¯hτ = ω¯h × ω¯τ : ω¯h = {ξi = ih, i = 0, 1, ...N, N h = ξ∗}, ω¯τ = {tn = nτ, n = 0, 1, ...M, M τ = 1}, h is the step in spatial coordinate, τ is the time step. Numerical solutions at grid nodes (xi, tn) are denoted by ϕin = ϕ(xi, tn). (8) (9) The iterative process is carried out using the following difference scheme ϕin+1 τ − ϕin = A(ϕin) [ ϕn+1 i−1 − 2ϕin+1 + ϕin++11 ] + B(ϕin) [ ϕin+1 − ϕin−1 h2 2h
To implement equation (10) by the sweep method, it is necessary to set the boundary conditions, and ξ∗ is determined in the course of numerical experiments according to condition (9).
Figures 1 – 2 show the dependence of the change in the porosity function on the self-similar variable.
Acknowledgements The work was carried out in accordance with the State Assignment of the Russian Ministry of Science and Higher Education entitled ‘Modern methods of hydrodynamics for environmental management, industrial systems and polar mechanics’ (Govt. contract code: FZMW-2020-0008, 24 January 2020). [Bea72] J. Bear. Dynamics of Fluids in Porous Media Elsevier, New York, 1972. 764 p. [Fri12]