Numerical Study of a Self-Similar Problem of Fluid Filtration in a Viscoelastic Medium in the Field of Gravity Margarita A. Tokareva Rudolf A. Virts Victoria N. Larionova Altai State University, Altai State University, Altai State University, Institute of Mathematics Institute of Mathematics Institute of Mathematics and Information and Information and Information Technologies Technologies Technologies Russia, Barnaul, Russia, Barnaul, Russia, Barnaul, Lenin Avenue 61, 656049 Lenin Avenue 61, 656049 Lenin Avenue 61, 656049 tma25@gmail.ru virtsrudolf@gmail.com lazylazo801@gmail.com. Abstract The quasilinear system of a composite type describing the spatial un- steady 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 consid- ered, and the viscoelastic properties of the deformable porous skeleton are also taken into account. The problem is reduced to a single third- order 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. 1 Introduction The relevance of the theoretical study of filtration problems in porous media is associated with their wide application in solving important practical problems. Examples are filtration near river dams, reservoirs and other hydraulic structures [Bea72]; irrigation and drainage of agricultural fields; oil and gas production, in particular, the dynamics of hydraulic fracturing cracks; problems of degassing coal and shale deposits to extract methane [Fow99]; movement of physiological fluids in tissues; tumor growth processes [Fri12], [Ast07]. The construction of mathematical models of such processes is complicated by the fact that the fluid flow is often considered in a mobile inhomogeneous medium, which is characterized by the presence of variable porosity. A special feature of the model of liquid filtration in a porous medium considered in this paper is the consideration of the mobility of the solid skeleton and its poroelastic properties. Interest in this problem also arises in connection with the widespread use of surface waves that occur in viscoelastic medium, when three waves propagate independently of each other: fast and slow longitudinal, as well as transverse [Con98]. Surface waves are studied in detail in Copyright ⃝ c by the paper’s authors. Use permitted under Creative Commons License Attribution 4.0 International (CC BY 4.0). In: Sergei S. Goncharov, Yuri G. Evtushenko (eds.): Proceedings of the Workshop on Applied Mathematics and Fundamental Computer Science 2021 , Omsk, Russia, 24-04-2021, published at http://ceur-ws.org 1 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 ∂(ρf ϕ) + div((1 − ϕ)ρs v⃗s ) = 0, + div(ρf ϕv⃗f ) = 0, ∂t ∂t kϕn ϕ(⃗vf − ⃗vs ) = − (∇pf − ρf ⃗g ), µ ∂pe ∇ · ⃗vs = −a1 (ϕ)pe − a2 (ϕ)( + ⃗vs · ∇pe ), ∂t ( ( )) ∂⃗vs ∂⃗vs ∗ ∇ptot = ρtot⃗g + div (1 − ϕ)η +( ) , ∂⃗x ∂⃗x ρtot = ϕρf + (1 − ϕ)ρs , ptot = ϕpf + (1 − ϕ)ps , pe = (1 − ϕ)(ps − pf ). 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 m ρs , ρf = const, a1 = ϕη , a2 = βϕ ϕb [Con98]: ∂ϕ ∂ + (ϕvf ) = 0, ∂t ∂x ∂(1 − ϕ) ∂ + (vs (1 − ϕ)) = 0, ∂t ∂x k ∂pf ϕ(vf − vs ) = − ϕn ( + ρf g), µ ∂x ∂vs ϕm ∂pe ∂pe =− pe − βϕ ϕb ( + vs ), ∂x η ∂t ∂x ( ) ∂ptot ∂ ∂vs = −ρtot g + 2η(1 − ϕ) . ∂x ∂x ∂x We pass in this system to dimensionless variables: t = t1 t′ , x = x1 x, vf = v1 vf′ , vs = v1 vs′ , pf = p1 p′f , ptot = p1 p′tot , pe = p1 p′e , ps = p1 p′s , (hereinafter, the primes are omitted), and also put: x1 kp1 kρf g p1 t1 2η ρs p1 t1 = ,α = ,β = ,γ = , λ = βf p1 , ζ = ,ρ = ,κ = . v1 µv12 t1 µv1 η x1 t1 ρf g ρf gx1 ρf Then equations can now be written in the form 2 ∂ϕ ∂(vf ϕ) + = 0, ∂t ∂x ∂(1 − ϕ) ∂(vs (1 − ϕ)) + = 0, ∂t ∂x ( ) ∂pf ϕ(vf − vs ) = −ϕn α +β , ∂x ( ) ∂vs ∂pe ∂pe = −ϕ pe γ − ϕ λ m b + vs , ∂x ∂t ∂x ( ) ∂ptot ∂ ∂vs κ =ζ (1 − ϕ) − ϕ − (1 − ϕ)ρ. ∂x ∂x ∂x 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ϕ −c = 0, (1) dξ dξ d d(1 − ϕ) ((1 − ϕ)vs ) − c = 0, (2) dξ dξ dpf ϕ(vf − vs ) = −αϕn − βϕn , (3) dξ dvs dpe dpe = −γϕm pe + λcϕb − λϕb vs , (4) dξ dξ dξ ( ) dptot d dvs κ =ζ (1 − ϕ) − ϕ − ρ(1 − ϕ). (5) dξ dξ dξ The system is supplemented with boundary conditions: vs (0) = vs0 , vf = vf0 , ϕ(0) = ϕ0 , lim ϕ(ξ) = ϕ+ , ξ→∞ lim vf (ξ) = u+ , lim ϕ(ξ) = ϕ+ , ξ→∞ ξ→∞ where vs0 , vf0 , ϕ0 , ϕ+ are given constants satisfying the conditions vs0 ̸= vf0 , ϕ0 ̸= ϕ+ . From the equations (1) – (2) of the system, we obtain: ϕ+ (1 − ϕ0 )vs0 − ϕ0 (1 − ϕ+ )vf0 c= , ϕ+ − ϕ0 u+ = vf0 ϕ0 + (1 − ϕ0 )vs0 , (1 − ϕ+ )ϕ0 (1 − ϕ0 )(vf0 − vs0 ) A2 = , ϕ+ − ϕ0 ϕ+ A1 = A2 . 1 − ϕ+ Thus, the system is converted to the following form for finding functions ϕ, pf , ptot : A1 A2 dpf ϕ( − ) = −ϕn (α + β), (6) ϕ 1−ϕ dξ 3 ( ) d 1 λA2 ϕb dpe A2 = −γϕm (ptot − pf ) − , (7) dξ 1−ϕ 1 − ϕ dξ ( ) dptot d d 1 κ = ζA2 (1 − ϕ) ( ) − ϕ − ρ(1 − ϕ). (8) dξ dξ dξ 1 − ϕ Next, we obtain the equation for the porosity function. We express from equations (6) and (8) the resulting dp system dξf and dpdξtot , 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 ϕ: ( ) ζλA2 ϕb−m d3 ϕ 1 1 1 ζ d2 ϕ + + + κγ (1 − ϕ)2 dξ 3 1 − ϕ γ ϕm (1 − ϕ) κ dξ 2 ( ( ) ( ) λ ϕb−m b−m−n 1 1−ρ 1 + A1 + n − b−m+1+ − γ 1−ϕ ϕ1+n ϕ (1 − ϕ) κ ϕ(1 − ϕ) ( ) ( )( )) 1−n+b−m 1 βκ − αρ b−m 1 dϕ −A2 + + + + ϕn ϕn−1 (1 − ϕ) ακ ϕ 1−ϕ dξ ( )( )3 ( ) 2 ζλA2 ϕb−m b−m 1 dϕ λζA2 ϕb−m b−m 4 d ϕ dϕ + +3 + + + κγ (1 − ϕ)3 ϕ 1−ϕ dξ κγ (1 − ϕ)2 ϕ 1 − ϕ dξ 2 dξ ( )( )2 1 ζ 2 1 m 1 dϕ + + − + (1 − ϕ) κ γ ϕ (1 − ϕ) 2 m γ ϕ 1+m dξ A1 −n 1 − ρ ϕ1−n β ρ + ϕ − ϕ− + − = 0. A2 κA2 1 − ϕ αA2 κA2 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 ϕ: ( )2 d2 ϕ dϕ A(ϕ) 2 + B(ϕ) + C(ϕ) = 0, dξ dξ where ( ) ϕ−m ζ A(ϕ) = + , γ(1 − ϕ)2 κ(1 − ϕ) ( ) ζ 2ϕ−m m ϕ−m−1 B(ϕ) = + − , κ(1 − ϕ)2 γ(1 − ϕ)3 γ (1 − ϕ)2 A1 −n ϕ1−n ϕ + (1 − ϕ)ρ β C(ϕ) = ϕ − − + . αA2 α(1 − ϕ) κA2 αA2 3 Numerical study 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: |ϕ(ξ∗ ) − ϕ+ | ≤ ε, (9) where ε – is the desired accuracy of the solution. The search for a solution is performed with the necessary accuracy using the condition |ϕn+1 i − ϕni | ≤ ε 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 ϕni = ϕ(xi , tn ). 4 The iterative process is carried out using the following difference scheme [ ] [ n ] ϕn+1 − ϕn i n ϕ n+1 i−1 − 2ϕn+1 i + ϕ n+1 i+1 n ϕi+1 − ϕni−1 2 i = A(ϕi ) + B(ϕi ) + C(ϕni ) = 0. (10) τ h2 2h The equation (10) is supplemented with the following conditions: ϕ0 = 0.5, ϕ+ = 0.75. The initial condition can be selected in two ways: ϕ+ − ϕ0 ϕ(x, 0) = 1/2, ϕ(x, 0) = ξ + ϕ0 . ξ∗ 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. 0.85 0.8 0.75 0.7 φ 0.65 0.6 0.55 0.5 0 2 4 ξ 6 8 10 Figure 1: Dependence of the porosity function on the self-similar variable in the case of the following values of the parameters that determine the constants g = 9.8m/s2 ; µ = 0.009P a · s; η = 108 P a · s; ρf = 1000kg/m3 ; ρs = 916kg/m3 ; βϕ = 10−8 P a−1 ; k = 10−8 m2 ; p1 = 105 P a; t1 = 1000s; v1 = 0, 01m/s; x1 = 10m. Conclusion A self-similar problem of filtration of a viscous fluid in a viscoelastic porous skeleton is considered. The original system is reduced to a third-order differential equation for the porosity function in the case of a viscoelastic medium. If the viscous properties of the skeleton prevail over the elastic ones, the original system of governing equations is reduced to a second-order differential equation for the porosity function. A numerical study of the second case is carried out. In the future, it is planned to study the equation for a medium with both viscous and elastic properties. 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). References [Bea72] J. Bear. Dynamics of Fluids in Porous Media Elsevier, New York, 1972. 764 p. 5 0.85 0.8 0.75 0.7 φ 0.65 0.6 0.55 0.5 0 2 4 6 8 10 ξ Figure 2: Dependence of the porosity function on the self-similar variable in the case t1 = α = β = γ = λ = κ = ζ = ρ = 1. [Fow99] A. C. Fowler and X. Yang Pressure solution and viscous compaction in sedimentary basins J. Geophys. Res. 1999. Vol. 104. P. 12,98912,997. [Fri12] A. Friedman Cancer as Multifaceted Disease Math. Model. Nat. Phenom. - 2012. - Vol. 7, no. 1. - P. 3-28. [Ast07] S. Astanin, A. Tosin Mathematical model of tumour cord growth along the source of nutrient Math. Model. Nat. Phenom. - 2007. - Vol. 2, no. 3. - P. 3-28. [Con98] J. A. D. Connolly and Yu. Yu. Podladchikov. Compaction-driven fluid flow in viscoelastic rock. Geodin. Acta. 1998. Vol. 11. P. 5584. [Bio56] M. A. Biot. Theory of propagation of elastic waves in fluid-saturated porous solid. J. Acoust. Soc. Amer. 1956. Vol. 28, no. 2. P. 168191. [Mor07] C. Morency. A numerical model for coupled fluid flow and matrix deformation with applications to disequilibrium compaction and delta stability. Journal of Geophysical Research. 2007. [Fow11] A. Fowler. Mathematical Geoscience. Springer-Verlag London Limited, 2011. [Tok16] M. A. Tokareva and R. A. Virts Analytical and numerical study of the filtration problem in a poroelastic medium MAK P.75–79. – 2016. (in Russian) [Pap21] A. A. Papin, M. A. Tokareva and R. A. Virts. Filtration of liquid in a non-isothermal viscous porous medium. Journal of Siberian Federal Universit. Mathematics and Physics. - 2020. - Vol. 13, no. 6. - P. 763-773. [Pap19] A. A. Papin and M. A. Tokareva. Global Solvability of a System of Equations of One-Dimensional Motion of a Viscous Fluid in a Deformable Viscous Porous Medium. Journal of Applied and Industrial Mathematics. - 2019. - Vol. 13, no. 2. - P. 350-362. [Kha08] G. S. Khakimzyanov and S. G. Chernyi. Computational Methods. Novosibirsk, 2008. (in Russian) 6