Development of Embedded Element Technique for Permeability Analysis of Cracked Porous Media

The widely used approach of mesoscale finite element modeling for permeability analysis is to simulate the matrix and cracks with continuum elements (CE), whereas this process brings technical difficulties in generating a satisfying mesh conformity at the interface. In this work, an alternative method based on embedded element (EE) technique is developed for the prediction of water pressure field and effective permeability in the numerical simulation. Based on the mathematical similarity between elasticity and seepage problems, water pressure can derive from the corresponding displacement through “elastic analogy.” To assess the capability of the EE technique, different cases are simulated and compared with the CE model. The results show that there is a satisfactory agreement in water pressures and velocities between the CE and EE modeling. In the CE model, different factors, such as permeability contrast between matrix and cracks (Kcrack/Kmatrix) and mesh size, are considered. It is obvious to find that results will become stable when Kcrack/Kmatrix reaches 10, and the mesh size has little impact. The effective permeability of 3D porous media with random cracks is evaluated and the results show that the differential method is accurate for 3D permeability analysis when the crack density is not large.


Introduction
During the period of services, porous media like concrete and rocks are often subjected to various loads caused by mechanical, thermal, or physicochemical factors; the damage process generates new microcracks [1][2][3].If the microcracks continue clustering and connecting, the permeability of porous media will certainly increase, accelerating the deterioration process and threatening the safety of structures.The cracked porous media are usually regarded as a two-phase composite which includes a porous matrix and a large number of cracks for permeability problems [4,5].The X-ray CT scan images of porous media have shown that the crack distribution varies with locations, influencing mechanical and permeability properties [6][7][8].The geometrical information of individual cracks, such as orientation and connectivity, should be considered on cracked porous media.
The two main categories of models used in permeability analysis are the equivalent continuum model (ECM) and the discrete fracture model (DFM) [9,10].The ECM is widely used in field-scale studies, and it contains a limited number of subdomains in which permeability is assumed to be uniform [11].Under this assumption, the effective permeability is an average value over the subdomains that reflects the comprehensive effects of cracks and matrix [12,13].Although the ECM is simple to implement and requires low computational resources, it is difficult to present a reasonable representative volume elementary volume for permeability [14,15].Whereas, in small-scale studies, the influence of the individual cracks cannot be neglected, the DFM are quite suitable on this condition.In this approach cracks are represented by line elements or surface elements for 2D problems, and for 3D problems cracks are simulated with surface elements or body elements [5,[16][17][18].Although this method considers real crack geometry, it is computationally intensive for models containing large number of cracks.The most commonly used numerical approaches in three-dimensional DFM are conventional solid or continuum elements [19,20].

Composite
Host part Embedded part This approach not only describes the actual crack distribution but also considers water transportation between the matrix and the cracks.Unfortunately, this process brings technical difficulties in generating a satisfying mesh conformity at the interface; therefore software is not able to mesh a geometry with more than 100 fractures when both the matrix and fractures must be mesh [21].
Noncontinuous meshing techniques present a practical solution for the meshing problems in finite element models applied in DFM [22].Fish [23] introduced a mesh superposition technique to modeling of laminated composites and stress-stain fields.Takano et al. [24] proposed an enhanced mesh superposition method for multiscale finite element analysis of porous materials.They discretized the domain into three parts: global, local, and pores parts.Based on automatic image-based modeling and FE mesh superposition method, Kawagai et al. [25] study the multiscale modeling strategy for complex and heterogeneous microstructures of real materials.The embedded element of the commercial software ABAQUS can be applied in DFM as a mesh superposition technique [22].Under the EE technique, the host and embedded parts can be meshed separately and independently, resolving difficulties in meshing complex geometry of the matrix.
Some noncontinuous methods like s-version FE model are developed and implemented by some authors, which makes it difficult for other researchers to use and develop them further.In this case, this paper adopts the EE technique to simulate pressure fields of cracked porous media.The results of meso-FE simulations of different cases using the CE and EE models are investigated and discussed to verify the EE technique.Meanwhile, some factors such as permeability ratio between matrix and cracks and mesh size which could influence results are also taken into consideration.Finally, effective permeability of 3D porous media with random distributed cracks are evaluated.

Embedded Element Technique
In the EE model, the reinforced part is placed inside the matrix part.The host is the main part, which is considered an independent model from the point of view of translational degrees of freedom (DOFs).The simple form of embedded element technique is shown in Figure 1.The host part and the embedded part can be meshed separately.
The weight functions in the famous finite element analysis software ABAQUS are used to create the geometric relationship between the nodes of embedded and host elements [26], which is illustrated as follows: where U (Emb) and U () denote the DOF of the embedded and host nodes and   denotes the weight function.The value of the weight function is determined by the distance between the embedded node and the corresponding host node.A short distance will lead to a large weight function.If an embedded node is located within a host element, then the translational DOFs at the node are eliminated and the node becomes an "embedded node."The translational DOFs of the embedded node are constrained to the interpolated values of the corresponding DOF of the host element.The host and embedded strain are set to have the equivalent in the superposition regions, which is called the principal of strain equivalence.Therefore, the composite stiffness S  can be calculated as follows [27]: where S  and S  represent the stiffness for the host and embedded part separately.Finally, the stress in the embedded region   can be calculated as where   is strain matrix and D  is the elastic matrix.The relationships between stresses of different parts are shown in Figure 2 in a concise way.
Even with the possibility of applying the EE technique in permeability analysis for cracked porous media, the following limitations must be eliminated first.Firstly, the EE techniques are allowed to have rotational DOF, but these rotations are not constrained by embedding.Secondly, some DOFs such as rotational, temperature, pore pressure, acoustic pressure, and electrical potential at an embedded node are not constrained [28].Therefore, the direct application of the EE does not work.In this case, an indirect approach called "elastic analogy" is adopted as an auxiliary method combined with the EE technique.
The relationship between stresses of different parts.

Elastic Analogy
Heat conduction, mass diffusion, and seepage can be described by very similar mathematical equations [29].The mathematical similarity also exists between elasticity and seepage problems, which is illustrated in Table 1.For 3D problems, the displacement vector has three components, and the water head is merely a scalar.So a seepage field can be obtained through an analogy of the corresponding elastic field in the case of constrained DOF.Hence EE technique for elastic problems can be extended to the corresponding seepage problems by "elastic analogy" because of such correspondence.
In steady-state seepage analysis, the flow rate q and water head  obey the generalized Darcy law in Cartesian coordinate system as follows: where K is the hydraulic conductivity tensor, which is expressed as follows: And K is a second-order tensor with six independent elements.The equilibrium equation can be expressed as follows: where  denotes the water source.If the principal axes of the materials coincide with the coordinate axes, then ( 6) can be rewritten as follows: Strain  and stress  in the elasticity problems obey Hooke's law as follows: where flexibility tensor C is a fourth-order tensor with 21 independent elements.In an orthotropic material, C becomes where   and   (,  = 1, 2, 3) denote the elastic and shear moduli and ] denotes Poisson's ratio.In the case of ] = 0, C becomes a diagonal matrix.Furthermore, if the principal axes of the materials coincide with the coordinate axes but the displacements in the  and  directions are constrained to be zero, then the elastic equilibrium equation in the  direction can be simplified as follows: where  and   are the displacement and body force in the -direction, respectively.This equation is similar to (7).By equating  13 ,  23 , and  33 to   ,   , and   , the water head and flow velocity can be obtained from the values of displacement  and stress components  13 ,  23 , and  33 , which are determined by (10).This equation is used as the basis for a seepage field to be obtained by the "degradation" of a corresponding elasticity field.

Numerical Application
This section presents the numerical results obtained using EE methods described in the previous sections.First, two examples which contain a single crack are addressed as the reference problems to verify the validity of the EE method by comparing with the results obtained from CE method.In the first one, a single crack in a 2D infinite matrix is modelled.Similarly, a 3D example with a single rectangular crack in a porous domain is presented in the second case.Moreover, an example with eight intersecting cracks is also presented to illustrate the capacities of the EE method.Finally, the effective permeability of 3D cracked porous media with hundreds of random distributed cracks is evaluated and compared with various analytical effective medium predictions.

2D Example of a Single Crack.
In this part, an inclined crack located in a rectangular matrix is shown in Figure 3 [10].The matrix has an isotropic permeability, and a constant conductivity is supposed for the crack along its direction.Two constant values of pressure are applied on the bottom and up edges, and a zero normal flux is prescribed on the left and right edges.The CE and EE methods are used to generate meshes for this case (Figure 3).It is obvious to find that the EE modeling has independent meshes between crack and matrix, and the grid around the single crack seems intensive as expectation.The pressure fields determined by the CE method and EE method are compared as shown in Figure 4.The contours obtained from these two methods have similar general characteristics which vary rapidly near the crack tips.However, unlike the displacements, water head varies continuously through the crack.The pressure profiles along two vertical sections ( = 0 and  = 3) through the domain were also created to compare the pressure change (shown in Figure 5).The results show that there is an excellent agreement between the two models, which proves the validity of EE method.

3D Example of a Single
Crack in a Cubic Domain.The 3D example contains a single square crack with 20 m long inside the cubic matrix, and the dimension of the matrix is 100 m × 100 m × 100 m (shown in Figure 6).The center of the square crack is located at the position (60, 60, and 60), and its normal direction is  → n = (0.63, −0.45, 0.63).The matrix and crack are both isotropic and their relative permeability coefficients are 1 and 10,000, respectively.On faces  ( = 100) and ℎ ( = 0), the water heads are set to be 1.0 and 0.0, and other faces are set to be impermeable.
The values of water heads at eight reference points are listed in Table 2; meanwhile the flow velocities are summarized in Table 3.The results shown in Table 2 reveal that difference of water head between the two methods is quite small, and the maximal relative error is less than 3%.In Table 3, the flow velocities in the -direction of the two methods are nearly the same; in addition the velocity components in  and  directions are so small that they can be omitted.However, what needs to be emphasized is that the relative errors of flow velocity at different reference points are all rather small except the point , which is very close to the crack tip.The accuracy of water head filed is excellent on the whole domain, and the accuracy of flow velocity field is good enough except the limited area near the crack tip.The good agreement between the EE and CE modeling indicates the EE model is accurate.

Effect of the Relative Permeability Ratio.
In order to examine the effect of permeability ratio between crack and matrix ( crack / matrix ) upon the flow properties, the 3D single crack model is studied.The variation of water pressure (water head) is plotted in Figure 9.
With the increase of  crack / matrix , water heads at different locations converge to constant values.When the ratio is larger than 10 4 , the water heads at the reference points reach stability.This result is attributed to the fact that when the opening of the crack is sufficiently large, the crack becomes a superconductor and the head losses that occur almost reach     zero as water flows through.In general, for crack-matrix composites such as cracked porous media, cracks usually can be regarded as superpermeable compared with the matrix [5].

Mesh Sensitivity.
Different mesh sizes of EE were investigated in the single crack model of the same 3D case to evaluate mesh sensitivity in EE modeling.The element size of the matrix remains constant (  = 5), and the mesh size of crack has three different values (  = 10, 5, 1).Therefore, mesh sensitivity analysis is investigated under three different mesh size ratios ( matrix/crack = 0.5, 1, 5). Figure 10 shows the water head fields with different mesh size ratios, and Table 4 provides the values of water heads at the eight different reference points.The very small differences among different results conclude that the influence of mesh size upon water head can be ignored, which are similar to results presented by Tabatabaei et al. [22].

Multiple Intersecting Cracks.
The case of intersecting cracks might be more interesting to us.Here consider a 6 m long square domain with nine intersecting cracks [30].
The model boundary condition and its corresponding CE and EE meshes are shown in Figure 11.The permeability of matrix and cracks is 1 × 10 −15 m 2 and 1 × 10 −9 m 2 , respectively.Both the lower and the upper boundary are impermeable.The left boundary has a constant pressure of 100 kPa, while the right boundary has a constant pressure of 0 Pa. Figure 12 shows that the displacement filed from the EE method is similar to the pressure filed from the CE method.In the center part of the domain, both fields have a relatively flat change, and the contour lines vary rapidly around different crack tips.The pressure values at different intersecting points are also investigated and Table 5 presents the results obtained from these two methods.It is found that the EE model provides a little higher values than the results obtained from the CE model, and the maximal relative error is less than 2%.The results conclude the EE technique is also feasible for intersecting cracks.

Effective Permeability of 3D Cracked Porous Media.
To study the effective permeability of 3D cracked porous media, cracks are supposed to be isodiametric discs and follow unclustered random distribution.With such assumptions, the permeability of the cracked porous media is evaluated for different crack densities.Since the crack radius is kept constant, the crack density or volume fraction only depends on the number of cracks.The permeability contrast of crack to matrix is set to be 10000 to represent the great permeability   gap between crack and matrix in the real situation.Constitutive equations of flow are given by Darcy's law, and this law is based on the establishment of a uniform water gradient along the domain as a result of applied external water flux.The side length of the cubic matrix and the radius of all cracks are taken as  = 100 and  = 10.In the modeling, two different water heads are given on two opposite sides of the matrix, and other sides are impermeable.The mesh schematic with the EE technique is displayed in Figure 13, where the matrix and the cracks are meshed separately and independently.Meanwhile, the crack density adopts definition of volume fraction of cracks with respect to the reference volume [5,21]: where the crack volume is denoted by the cubic of its radius, and the reference volume is the cubic of the side length.The predicted results obtained from effective medium models, such as the dilute schemes, the differential scheme, and selfconsistent method, are compared with numerical results to investigate the potential application of the EE technique in DFM (shown in Figure 14).
The results shows that different approximations fit numerical results very well for crack density  ⩽ 0.04.As crack density increases, the dilute and difference methods seem to deviate the numerical solutions.On the other hand, the self-consistent approach still provides a relatively accurate fit to the numerical data.The effective medium methods are based on the assumption of a single inclusion in an infinite domain, so the effect of intersecting could not be taken into consideration.It is generally accepted that selfconsistent approach overestimates the effective property for nonintersecting cracks.However, self-consistent approach might be more reasonable when cracks comes intersecting.

Conclusion and Discussion
In this paper, we have developed a novel numerical method to calculate permeability of cracked porous media.The EE technique, combined with elastic analogy, can be utilized in the permeability modeling for crack-matrix composites like concrete or rocks.Compared to pressure (water head) and flow velocity fields obtained from the conventional CE method, the reliability of the EE technique is evaluated with different cases, and the results prove its validity and accuracy.Furthermore, some factors that influence the EE modeling are also taken into account.Further examinations indicate that stable calculation results can be induced if the relative permeability ratio reaches 10 4 .What is more, the influence of mesh size upon water pressure is not significant.
For the case of 3D cracks densely and randomly distributed, the EE technique consequently reduces the computational cost since matrix and cracks can be meshed separately and independently.The numerical results show that the EE method can be applied to simulate the effective permeability of cracked porous media, and the self-consistent approach has a good analytical estimation when cracks come intersecting.

Figure 1 :
Figure 1: The relationship between host and embedded part.

Figure 3 :
Figure 3: Two-dimensional cracked domain with a single crack: boundary condition (a); mesh of CE method (b); mesh of EE method (c).

Figure 4 :Figure 5 :
Figure 4: Pressure and displacement fields obtained from CE (a) and EE method for a single crack domain (b).

Figure 6 :Figure 7 :
Figure 6: Three-dimensional cracked domain with a single crack.

Figure 8 :Figure 9 :
Figure 8: Water pressure fields calculated for the single crack model.

Figure 11 :
Figure 11: Two-dimensional cracked domain with multiple intersecting cracks: boundary condition (a); mesh of CE method (b); mesh of EE method (c).

Figure 12 :Figure 13 :
Figure 12: Pressure and displacement fields obtained from CE and EE method for multiple intersecting cracks domain.

Table 1 :
Correspondence between elastic and conductive problems.

Table 2 :
Water heads at different reference points.

Table 3 :
Flow velocities at different reference points.

Table 4 :
Water head results under different mesh sizes.

Table 5 :
Pressure at different intersecting points.
Figure 14: Effective permeability of cracked porous media with 3D random cracks. of China (2017YFC0804602), National Basic Research Program of China (2013CB035902), National Natural Science Foundation of China (51339003), and State Key Laboratory of Hydroscience and Engineering (2016-KY-05).