A 2D self-organized percolation model for capillary impregnation
Anh Khoa Nguyen, C Trang, E. Blond, E. de Bilbao, T. Sayet, A. Batakis
1 CSMA 2019
A 2D self-organized percolation model for capillary impregnation
A. K. Nguyen , C. B. Trang , E. Blond , E. De Bilbao , T. Sayet and A. Batakis Univ. Orléans, LaMé (EA7494, Univ. Orléans, Univ. Tours, INSA CVL), France, {anh-khoa.nguyen, eric.blond, thomas.sayet, cong-bang.trang}@univ-orleans.fr Univ. Orléans, Insitut Denis Poisson (UMR CNRS, Univ. Orléans, Univ Tours), France, [email protected] Univ. Orléans, CEMHTI (UPR CNRS), France, [email protected]
Abstract —A two-dimension extension of the S elf-organized G radient P ercolation ( SGP ) method initially developed for the one-dimensional simulation is proposed. The initialization in the two directions is considered as the analytic solution of the 2D (homogeneous) diffusion equation. The evolution of the saturation front is assumed to be the evolution of both standard deviations in each direction. The validation of the implementation is done by comparisons between SGP and finite element results.
Keywords — Impregnation, porous media, gradient percolation, simulation.
1. Introduction
To simulate the unsaturated (non-reactive) impregnation in a porous medium, classical numerical methods (FEM, theta-method, etc.) have often been employed to solve the Richards’ equation [1], [2]. Yet, the CPU time and spurious oscillations become prohibitive, as for instance, when the impregnation is coupled with chemical reactions. That is why a drastically different, and more reliable approach, called S elf-organized G radient P ercolation ( SGP ) method, was proposed [3]. The initialization of the algorithm is driven by an analytic solution of the (homogeneous) diffusion equation. It is a convolution between a P robability D ensity F unction ( PDF ) and a smoothing function [4]. The evolution of the capillary pressure profiles with time is reproduced by the self-evolution of the standard deviation of the PDF. To test this model, the capillary pressure profiles and the mass gain curve have been confronted with those obtained by FEM and experimental measurement, respectively [3]. The main goal of the present work is to extend the SGP method to the 2D unsaturated impregnation only considering the non-reactive case. The two-dimensional SGP algorithm has the same structure than the one-dimensional one. On the mathematical side, the analytic form is a convolution between a PDF of two directions and a smoothing function. In addition, the evolution of the saturation front is proposed to be the evolution of the two standard deviations of the PDF in each direction. The visualizations and comparisons between the numerical results from SGP method and FEM are given.
2. Self-organized Gradient Percolation (SGP) method for impregnation
The gradient percolation method is a probabilistic approach to predict a spatial evolution. The idea for using this method here is to avoid the numerical difficulties of the resolution of the Richards’ equation derived from the mass balance combined with Darcy’s law [3]. Richards’ equation is a nonlinear P artial D ifferential E quation ( PDE ), which requires a small-time increment and an appropriate fine space discretization to make sure the accuracy of the results. However, a small-time increment leads to a high computational cost and, moreover, a very small time-increment can give rise to spurious oscillations that impact the accuracy of the results [5]. The main goal for the development of this method is to reduce the CPU time and the number of parameters actually used for classical models. The S elf-organized G radient P ercolation ( SGP ) method, which is based on the gradient percolation method, has been developed to reach these goals. Two fundamental types of percolation are broadly used to simulate physical phenomena, the bond percolation and the site percolation [6]. The site percolation can be applied to the simulation of the non-reactive impregnation process for the determination of the saturation, first considered as a state variable describing the state (“occupied” or “empty”) of each local pore of the media. Considering a random network defined by numerous cells/squares, function 𝑋𝑋 ( 𝑧𝑧 ) , where 𝑧𝑧 is the position of the cell, is a random value uniformly distributed on the interval [0,1]. A site declared to be “occupied” or to be “empty” by the liquid is defined in comparison with constant probability 𝑝𝑝 as: • site 𝑧𝑧 is “occupied” when 𝑋𝑋 ( 𝑧𝑧 ) ≥ 𝑝𝑝 ; • site 𝑧𝑧 is “empty” otherwise. In classical percolation model, probability 𝑝𝑝 , which indicates the state of a site, is independent of the site, it is the same value for the whole domain. However, physics implies that the state of a site is given by capillary pressure, which is not constant on the lattice. From physics point of view, the C apillary P ressure C urve ( CPC ) is at equilibrium only under the gravity field. The capillary pressure is the equilibrium pressure, which is determined by the difference between the non-wetting and the wetting phases. In the capillary rising test, the saturation profile along the vertical direction (Figure 1), considered as a function of the level (height) with respect to the liquid regarding Jurin’s law, is an image of the CPC at each time 𝑡𝑡 . That is why we propose to transform probability 𝑝𝑝 into a function of the location of the site. This can be done by adapting gradient percolation theory developed by Pierre Nolin [7]. Choosing 𝑝𝑝 as a function of the location allows reproducing the impregnation gradient, but the results will not respect the continuity of the liquid phase. Moreover, it is impossible to manage the boundary conditions. It is thus proposed to define a convolution product to satisfy the continuity and the possibility to take into account various boundary conditions. Finally, the local saturation is obtained through two steps: (1) local state 𝑋𝑋 is determined by the classical percolation and (2) employing the convolution product gives: 𝑆𝑆 ( 𝑧𝑧 ) = 𝑋𝑋 ( 𝑧𝑧 ) ∗ 𝛿𝛿 ( 𝑧𝑧 ) ( ) where 〈∗〉 is the convolution operator and 𝛿𝛿 is a smoothing function (spatial weighted average). Figure 1 – The saturation profile as a probability density function
The saturation profile at initial time 𝑡𝑡 is then proposed by the following expression [4]: 𝑋𝑋�𝑃𝑃 𝑐𝑐𝑐𝑐𝑐𝑐 , 𝑡𝑡 � = 𝑆𝑆 𝑟𝑟 + ( 𝑆𝑆 𝑠𝑠 − 𝑆𝑆 𝑟𝑟 ) 𝑒𝑒𝑒𝑒𝑝𝑝 �− �𝑃𝑃 𝑐𝑐𝑐𝑐𝑐𝑐 − 𝑃𝑃 𝑐𝑐𝑐𝑐𝑐𝑐 , 𝑡𝑡 𝑆𝑆 𝑚𝑚𝑚𝑚𝑚𝑚 � 𝑚𝑚 𝑚𝑚𝜎𝜎 𝑚𝑚 � ( ) where 𝑆𝑆 𝑟𝑟 and 𝑆𝑆 𝑠𝑠 mean the residual and maximum saturation, respectively, 𝑃𝑃 𝑐𝑐𝑐𝑐𝑐𝑐 is the local capillary pressure, 𝑃𝑃 𝑐𝑐𝑐𝑐𝑐𝑐 , 𝑡𝑡 𝑆𝑆 𝑚𝑚𝑚𝑚𝑚𝑚 is the minimum pore pressure to initiate the impregnation, 𝜎𝜎 is the standard deviation of the PDF, 𝑚𝑚 is the empirical parameter to switch between probability distributions to fit the CPC (e.g. 𝑚𝑚 = 1 and 𝑚𝑚 = 2 designate Laplace and Gaussian distributions, respectively). From the physics point of view, 𝑚𝑚 is a dimensionless parameter which depends on the porosity network (such as porosity, tortuosity, etc.). To reproduce the physics of impregnation, it is necessary that the model can evolve autonomously with time following the physical laws of the problem. In order to ensure this self-evolution along time, a relationship between the standard deviation and the capillary pressure is proposed [3]: 𝜎𝜎 𝑛𝑛 �𝑃𝑃 𝑐𝑐𝑐𝑐𝑐𝑐 � = 𝜎𝜎 𝑛𝑛−1 �𝑃𝑃 𝑐𝑐𝑐𝑐𝑐𝑐 � + 𝐴𝐴 �𝑃𝑃 𝑐𝑐𝑐𝑐𝑐𝑐 , 𝑡𝑡 𝑆𝑆 𝑃𝑃 ℎ , 𝑛𝑛 − � ( ) where 𝑛𝑛 designates the time increment and 𝐴𝐴 designates a kinetic constant ( 𝑚𝑚 ∙ 𝑠𝑠 −1 ). Finally, there are only three parameters required for the SGP algorithm: 𝐴𝐴 ( 𝑚𝑚 ∙ 𝑠𝑠 −1 ), 𝑃𝑃 𝑐𝑐𝑐𝑐𝑐𝑐 , 𝑡𝑡 𝑆𝑆 (Pa) and 𝑚𝑚 (dimensionless). To identify these three parameters, the inverse identification method [3] is employed by fitting the mass gain curve obtained from the SGP method with the experimental one and then is validated by comparison with FEM simulation if the capillary pressure profiles are known. Then, 𝐴𝐴 is closed to the intrinsic permeability (but not fully explained today) while 𝑚𝑚 is directly linked to the shape of the capillary pressure curve and can be taken as the inverse of the parameter of van Genuchten’s model. Therefore, only unknown 𝑃𝑃 𝑐𝑐𝑐𝑐𝑐𝑐 , 𝑡𝑡 𝑆𝑆 is the minimum pore pressure to initiate impregnation.
3. SGP method for 2D problem
Since the SGP algorithm presented above for the (quasi) one-dimensional case permits to reduce the CPU time and to avoid the spurious oscillations [3]. Its extension to the 2D case thus keeps working on and is straightforward to construct: the main parameter will no longer be the height but the total pressure. It is to note that the Richards’ equation, where the height parameter is replaced by the total pressure parameter since from the physics point of view in the general case, is written as [3]: ϕ 𝜕𝜕𝑆𝑆𝜕𝜕𝑡𝑡 = 𝑑𝑑𝑑𝑑𝑑𝑑 �𝐾𝐾 ( 𝑆𝑆 ) 𝜂𝜂 𝑔𝑔𝑔𝑔𝑔𝑔𝑑𝑑 ( 𝑃𝑃 𝑡𝑡𝑡𝑡𝑡𝑡 ) � ( ) where 𝜙𝜙 is the porosity (dimensionless), 𝑆𝑆 is the saturation (dimensionless), 𝐾𝐾 ( 𝑆𝑆 ) is the nonlinear function of saturation ( 𝑚𝑚 ) , 𝜂𝜂 is the dynamic viscosity ( 𝑃𝑃𝑔𝑔 ∙ 𝑠𝑠 ) , 𝜌𝜌 𝑊𝑊 is the mass density of the wetting liquid ( 𝑘𝑘𝑔𝑔 ∙ 𝑚𝑚 −3 ) ; 𝑔𝑔𝑔𝑔𝑔𝑔𝑑𝑑 ( 𝑃𝑃 𝑡𝑡𝑡𝑡𝑡𝑡 ) = 𝑔𝑔𝑔𝑔𝑔𝑔𝑑𝑑�𝑃𝑃 𝑐𝑐𝑐𝑐𝑐𝑐 + 𝜌𝜌 𝑊𝑊 𝑔𝑔𝑧𝑧� is not a spatial description, but rather depends on the direction of gravity, where 𝑃𝑃 𝑡𝑡𝑡𝑡𝑡𝑡 is the total pressure ( 𝑃𝑃𝑔𝑔 ) , 𝑔𝑔 is acceleration due to the gravity ( 𝑚𝑚 ∙ 𝑠𝑠 −2 ) and 𝑧𝑧 is the height along the gravity direction ( 𝑚𝑚 ) . Three main steps of the 2D SGP algorithm are presented as follows.
Step 1:
Based on the main idea of solving the 1D (homogeneous) diffusion equation [3], the analytic solution of the 2D case, even in the general case, is pointed out to be a convolution between a Probability Density Function and a smoothing function. To be physically meaningful, the initialization is thus taken to be the analytic solution of the (homogeneous) diffusion equation (only considering initial time 𝑡𝑡 ) in the vertical and horizontal directions ( 𝑧𝑧 , 𝑧𝑧 ) as the following: 𝑋𝑋 ( 𝑃𝑃 𝑡𝑡𝑡𝑡𝑡𝑡 ( 𝑧𝑧 ), 𝑃𝑃 𝑡𝑡𝑡𝑡𝑡𝑡 ( 𝑧𝑧 ), 𝑡𝑡 ) = 𝑆𝑆 𝑟𝑟 + ( 𝑆𝑆 𝑠𝑠 − 𝑆𝑆 𝑟𝑟 ) 𝑒𝑒𝑒𝑒𝑝𝑝 �− �𝑃𝑃 𝑡𝑡𝑡𝑡𝑡𝑡 ( 𝑧𝑧 ) − 𝑃𝑃 𝑡𝑡𝑡𝑡𝑡𝑡 , 𝑡𝑡 𝑆𝑆 𝑚𝑚𝑚𝑚𝑚𝑚 � 𝑚𝑚 𝑚𝑚𝜎𝜎 ( 𝑧𝑧 ) 𝑚𝑚 − �𝑃𝑃 𝑡𝑡𝑡𝑡𝑡𝑡 ( 𝑧𝑧 ) − 𝑃𝑃 𝑡𝑡𝑡𝑡𝑡𝑡 , 𝑡𝑡 𝑆𝑆 𝑚𝑚𝑚𝑚𝑚𝑚 � 𝑚𝑚 𝑚𝑚𝜎𝜎 ( 𝑧𝑧 ) 𝑚𝑚 � ( ) where 𝑃𝑃 𝑡𝑡𝑡𝑡𝑡𝑡 , 𝑃𝑃 𝑡𝑡𝑡𝑡𝑡𝑡 , 𝑡𝑡 𝑆𝑆 𝑚𝑚𝑔𝑔𝑒𝑒 mean the total pressure and the minimum pore pressure to start the impregnation assuming to be the same in each direction, respectively, 𝑚𝑚 is the empirical parameters defining the shape of the distribution in each direction and 𝜎𝜎 ( 𝑧𝑧 𝑖𝑖 ) is the standard deviation function along the 𝑑𝑑 -direction. To take into account Eq. (1) and Eq. (5), we employ the dedicate function “norminv” in Matlab and then a multistate gradient percolation model is then given. Step 2:
In the next time steps, we assumed that the evolution of the capillary front along each direction is the evolution of each standard deviation, respectively. In particular, standard deviation 𝜎𝜎 𝑛𝑛 ( 𝑧𝑧 𝑖𝑖 ) along 𝑑𝑑 -direction at timestep 𝑛𝑛 can be written as: 𝜎𝜎 𝑛𝑛 ( 𝑧𝑧 𝑖𝑖 ) = 𝜎𝜎 𝑛𝑛−1 ( 𝑧𝑧 𝑖𝑖 ) + 𝛼𝛼�𝑃𝑃 𝑡𝑡𝑡𝑡𝑡𝑡𝑛𝑛 ( 𝑧𝑧 𝑖𝑖 ) � ( ) where 𝛼𝛼�𝑃𝑃 𝑡𝑡𝑡𝑡𝑡𝑡𝑛𝑛 ( 𝑧𝑧 𝑖𝑖 ) � means the evolution of the standard deviation at timestep 𝑛𝑛 along 𝑑𝑑 -direction. The link with the Poiseuille’s equation has been pointed out in the 1D capillary rising case, which gives: 𝛼𝛼�𝑃𝑃 𝑡𝑡𝑡𝑡𝑡𝑡𝑛𝑛 ( 𝑧𝑧 𝑖𝑖 ) � = 𝛾𝛾𝐷𝐷 𝜂𝜂 𝑃𝑃 𝑡𝑡𝑡𝑡𝑡𝑡𝑛𝑛 ( 𝑧𝑧 𝑖𝑖 ) ( ) where 𝐴𝐴 = 𝛾𝛾𝐷𝐷 𝜂𝜂⁄ depends on capillary diameter 𝐷𝐷 ( 𝑚𝑚 ) and 𝛾𝛾 is the surface tension ( 𝑁𝑁 𝑚𝑚⁄ ) . Step 3:
In the final step, note that the motion of liquid here is due to driving force 𝑔𝑔𝑔𝑔𝑔𝑔𝑑𝑑�𝑃𝑃 𝑐𝑐𝑐𝑐𝑐𝑐 + 𝜌𝜌 𝑊𝑊 𝑔𝑔𝑧𝑧� . Obviously, only in the capillary rising case, the simulation will stop when there is no difference between capillary pressure and hydrostatic pressure [3]. Nevertheless, in the general case, the simulation will not be able to stop until there is no more liquid on the boundary of the source. For example, the horizontal impregnation is never stopped, since there is no effect of the gravity in this direction. In the numerical practice, we thus suggest comparing the difference between the mass of liquid in the container (denoted Π ) and the mass gain in time in the porous medium (denoted 𝑚𝑚 𝑡𝑡 ). It leads to that the SGP algorithm will stop when there is no difference in this comparison. 𝜕𝜕Ω 𝑖𝑖 corresponds to the three types of boundary conditions [3]. At the initial timestep, depending on the treated problem, it is considered that almost one surface of the porous sample is in contact with the surface of liquid. This leads to impose that concerning row (first row) of the matrix is fully filled with liquid. A specific boundary condition will be applied on each interface. The choice of the boundary conditions depends on the interpretation of the physical phenomena. Figure 2 – Flowchart of the SGP algorithm for the 2D problem.
4. Application
In order to numerically simulate the test using FEM and the SGP algorithm, we decided to build the two models implemented both in (2D) Abaqus standard (v. 6.14) and Matlab (v. R2019a) software, respectively. To compare the CPU time of the two methods, it is essential that both models have to be built with the same mesh size, time increment, dimensions and boundary conditions.
For the application of the SGP algorithm extended to the 2D simulation of the impregnation, considering a numerical example, a cylindrical porous sample is in (partial) contact with a liquid at the bottom surface. The height and width of the sample are 0.04 m and 0.035 m (in Figure 3) respectively. The porous material is Alumina 99% and the impregnated liquid is glycerine. resdidual =
Π − 𝑚𝑚 𝑡𝑡 L oop on ti m e s t e p 𝑡𝑡 residual ≠ residual = 0 DATA 𝑃𝑃 𝑡𝑡𝑡𝑡𝑡𝑡 , 𝑡𝑡 𝑆𝑆 𝑚𝑚𝑚𝑚𝑚𝑚 , 𝜎𝜎 , 𝛾𝛾 , 𝐷𝐷 , 𝑚𝑚 , , 𝜂𝜂 , Π 𝜎𝜎 𝑖𝑖𝑛𝑛 = 𝜎𝜎 𝑖𝑖𝑛𝑛−1 + 𝛼𝛼�𝑃𝑃 𝑡𝑡𝑡𝑡𝑡𝑡𝑛𝑛 ( 𝑑𝑑 ) � 𝑋𝑋�𝑃𝑃 𝑐𝑐𝑐𝑐𝑐𝑐 ( 𝑧𝑧 ), 𝑃𝑃 𝑐𝑐𝑐𝑐𝑐𝑐 ( 𝑧𝑧 ), 𝑡𝑡 � (detailed in Eq. ( ) ) CONVOLUTION OPERATOR 𝑆𝑆 ( 𝑡𝑡 , 𝑃𝑃 𝑡𝑡𝑡𝑡𝑡𝑡 ( 𝑧𝑧 )) = 𝑋𝑋 ( 𝑡𝑡 , 𝑃𝑃 𝑡𝑡𝑡𝑡𝑡𝑡 ( 𝑧𝑧 )) ∗ 𝛿𝛿 ( 𝑧𝑧 ) END
START
Figure 3 – The specific dimensions and boundary conditions of the numerical test are used both for the models using FEM and SGP method. The mesh size is 0.03 m (along surface 𝜕𝜕𝜕𝜕 ), 0.025 m and 0.0267 m (along surfaces 𝜕𝜕𝐴𝐴 in the horizontal and vertical directions, respectively).
Due to the lack of experiment data, the validation of the performance of the 2D SGP algorithm is made by comparisons of the local saturation of each square of the model using FEM with the one of each square with the same coordinates of the SGP model. Obviously, the number of the squares of each model is equal to the number of the squares in row multiplied with the one in the column. It leads to compare a very large number of squares in the case of very smooth mesh. To be simple in our ongoing study, we decided to consider the model with a coarse grid, but, of course, must respect the condition between time step and mesh size to avoid the spurious problems [5]. For the boundary conditions in both models, there is an imposed pressure on boundary surface 𝜕𝜕𝜕𝜕 and the free draining on boundary surfaces 𝜕𝜕𝐴𝐴 used in the simulation (Figure 3). In particular, for the first one, the liquid is in contact with surface 𝜕𝜕𝜕𝜕 , so it corresponds to an imposed total pressure; for the second one, to make sure that the liquid can be free to flow out of these surfaces, the drainage-only boundary condition is imposed on surface 𝜕𝜕𝐴𝐴 . The “expression” of these boundary conditions in the SGP model is done thanks to the convolution procedure [3].
Table 1- CPU time and time ratios between FEM and SGP method.
Durations (in seconds) CPU time (in seconds) Time ratios
FEM SGP method 142 15 0.075 200 4254 48.3 0.5367 90 𝜕𝜕𝐴𝐴 𝜕𝜕𝐴𝐴 𝜕𝜕𝐴𝐴 𝜕𝜕𝐴𝐴 𝜕𝜕𝐴𝐴 𝜕𝜕𝜕𝜕
Figure 4 –Simulations using FEM and SGP method for the evolution of the saturation fronts at different time from (a) to (d). The absolute errors (e) and (f), calculated by the absolute minus of the local saturation at a square (a) (b) (c) (d) (e) (f) Absolute error at 4254s of the model using FEM with the local saturation of the square with the same coordinates of the SGP model at different time, is resulted in the scale color bar.
The computational costs for the test for the SGP method and FEM are obtained by running the software on 2.60 GHz, x64-based PC with 32 GB of RAM. Obviously, the computational time of the SGP algorithm is significantly lower than the one of FEM (Table 1). This variation in the ratio is due to the difference in the evolution of time step between both methods [3].
5. Conclusion
In this work, our ongoing study of the extending 2D SGP method to the impregnation process is proposed. The preliminary visualizations from the FEM and SGP method (Figure 4 (a) to (d)) seems to indicate that the behaviour of the evolution of the saturation fronts from the two methods are the same at the first time steps (Figure 4 (e) and (f)). Hence, there is still a huge work to do in order to study the mesh sensitivity, the impact of the time step, the impact of each of the parameters required for the SGP model and their link with the physics (properties of the material). For the perspectives, the SGP algorithm is expected to adapt itself to taking into account the kinetics of the liquid motion in respect with the gravity direction. Then, to validate the results in this case, it is necessary to fit the kinetics of the liquid motion along the porous sample by comparing with the mass gain curves from FEM and/or experimental data at each separate zone regarding the gravity direction of the porous sample.
5. References [1] L. Bergamaschi and M. Putti, ‘Mixed finite elements and Newton-type linearizations for the solution of Richards’ equation’,
Int. J. Numer. Methods Eng. , vol. 45, no. 8, pp. 1025–1046, Jul. 1999. [2]
R. Kodešová, J. Šimůnek, A. Nikodem, and V. Jirků, ‘Estimation of the Dual -Permeability Model Parameters using Tension Disk Infiltrometer and Guelph Permeameter’,
Vadose Zone J. , vol. 9, no. 2, p. 213, 2010. [3] A. K. Nguyen, E. Blond, T. Sayet, A. Batakis, E. de Bilbao, and M. D. Duong, ‘Self-organized gradient percolation method for numerical simulation of impregnation in porous media’,
Comput. Methods Appl. Mech. Eng. , vol. 344, pp. 711–733, Feb. 2019. [4]
Computer Vision and Applications . Elsevier, 2000. [5] P. A. Vermeer and A. Verruijt, ‘An accuracy condition for consolidation by finite elements’,
Int. J. Numer. Anal. Methods Geomech. , vol. 5, no. 1, pp. 1–14, Jan. 1981. [6] G. Grimmett,
Percolation , vol. 321. Berlin, Heidelberg: Springer Berlin Heidelberg, 1999. [7] P. Nolin, ‘Critical exponents of planar gradient percolation’,