IntroductionCracks are common defects in quasibrittle materials such as rocks and concretes (Adler et al. 2013). These cracks can have different origins and result from various actions, e.g., the residual stress arisen during manufacturing processes (Li et al. 2020), fatigue or thermal loadings during service (Nguyen et al. 2020; Liu et al. 2021), or the aging of materials under internal and/or external processes (Torrijos et al. 2010). In porous solids, cracks facilitate the flow of fluids (liquids and gas) by creating additional transport channels and enhancing the connectivity of the pore structure in the solid matrix (Adler et al. 2013). Individual cracks can propagate and intersect each other, forming a crack network that could further promote the fluid flow (Adler and Thovert 1999; Li and Li 2015). With certain degree of clustering, the cracks may serve as preferential flow paths for aggressive agents, accelerate the aging processes of materials, and, consequently, affect the service lives of structures using these materials (Li 2016). In other cases such as extracting geothermal energy from rocks, the maximized overall flow capacity is sought as a desirable attribute by generating artificial crack networks in porous rocks (Hu et al. 2021). Therefore, it is a primary issue is to describe the impact of cracks on transport properties of materials and further on the target performance of these materials for engineering applications (Timothy and Meschke 2017b).The flow of fluids with low Reynolds number in a fully saturated porous solid can be described by Darcy’s law, which can also be derived by Navier-Stokes equation applied to steady flow on the scale of pores (Whitaker 1986). Introducing cracks into the porous matrix can have two consequences: physically, new channels, characterized by crack geometry, are added to pore channels for fluid flow (Wang et al. 2018), and chemically, a new interface (crack surface) is created, making possible the mass exchange between crack surface and flowing fluids such as the dissolution of solid phases or precipitation on the crack surface (Li et al. 2011). In this study, the chemical effect of a crack surface is neglected and only the physical effect is taken into account(Vinci et al. 2014; Mustapha et al. 2020). Thus, a cracked porous solid can be idealized in two phases: the cracks and the (homogeneous) porous matrix, and estimating its effective permeability requires averaging properly the permeabilities of the two phases, simultaneously taking into consideration the geometry of cracks (Guàuen et al. 1997). Two theories can help the estimation: effective medium theory (EMT) and percolation theory (PT).The effective medium theory develops various schemes to estimate the effective permeability of the equivalent homogeneous medium by evaluating the perturbation induced by inclusions. From the view of EMT, cracks could be taken as inclusions with high local permeability embedded in a low-permeable matrix (Eshelby 1957). Through appropriate assumptions of the far-field boundary conditions, EMT can further consider the interactions among multiple inclusions with multiple schemes available (Guàuen et al. 1997). The dilute scheme takes cracks to be nonoverlapping and positioned in homogeneous matrix without any interaction among cracks (Shafiro and Kachanov 2000). The self-consistent (SC) approximation idealizes a single crack embedded in the sought homogenized medium (Lemarchand et al. 2010; Jeannin and Dormieux 2017), and the SC scheme has been successfully applied to deal with stress–permeability coupling (Dormieux et al. 2011; Wang et al. 2018).The differential scheme assumes that cracks are positioned into the homogeneous matrix one by one, and the effective permeability obtained at each step is considered as the host matrix for the subsequent step (Norris 1985). The interaction direct derivative (IDD) scheme (Zheng and Du 2001) provides an explicit and simple structure for the prediction of permeability of cracked solids with low crack extent (Zhou et al. 2011; Li and Li 2019). The cascade continuum micromechanics (CCM) model considers the complexity of the embedded cracks through the hierarchization and the interaction of the cracks (Timothy 2016; Timothy and Meschke 2016). The multipole expansion method gives a complete solution to the boundary-value problem of the matrix–inclusion composites, comprising the continuous matrix and nonintersecting crack inclusions (Kushch 2020). The permeability upscaling scheme adapts Eshelby’s inhomogeneity problem for linear elasticity to the Darcy’s law (Dormieux and Kondo 2005), and the crack interaction is estimated by the SC scheme, which has been applied to cementitious materials (Pivonka et al. 2004; Damrongwiriyanupap et al. 2017) and bone tissue (Abdalrahman et al. 2015; Estermann and Scheiner 2018).Effective medium theory estimates the effective permeability from averaging; thus, it has its limit, i.e., the EMT models apply to situations where the perturbance of cracks to the fluid flow is moderate. As cracks perturb the fluid flow to a larger extent, e.g., by higher crack density or important clustering, the EMT begins to lose its validity (Guàuen et al. 1997). This is the very starting point for percolation theory to intervene, taking the cracks no longer as homogeneous inclusions but as networks. A central task of the percolation theory is to predict the percolation threshold (limit) of fluid flow and estimate the overall permeability when the crack perturbance is large (Hunt and Ewing 2009).The classical PT is established through lattice methods and has been applied to less regular or random networks through equivalent lattices (Langlois et al. 2018; Li et al. 2017). The continuum PT, more adapted to solids with random networks, has been used to evaluate the connectivity degree (Berkowitz 1995) and predict the effective permeability (Li and Li 2015). The percolation threshold and scaling exponent depend on the finite domain size and fractal morphology of random network (Liu and Regenauer-Lieb 2021). The correlation between the crack geometry and the permeability from percolation theory, lattice or continuum-based, has been reviewed by Li and Li (2015). PT intends to describe the flow behavior near the percolation threshold and cannot be applied to cases far below the threshold (Guàuen et al. 1997). Accordingly, the EMT and PT have their respective merits, suitable to address two characteristic ranges: the fluid flow is dominated by the solid matrix or the crack network.On the basis of the aforementioned state of the art, this paper focuses on the comprehensive modeling of permeability of microcracked porous solids for the entire range of crack density, employing the EMT and PT models to estimate the overall permeability in their respective valid ranges. In the range of moderate crack perturbance, the IDD model from EMT was used and extended to consider the finite crack connectivity and crack opening. As crack pertubance becomes strong with the appearance of dominant/spanning clusters of cracks, a percolation model was established to estimate the effective permeability considering the permeability of porous matrix and the geometry parameters of crack clusters. Following this logic line, the paper is presented as follows: the geometry characterization of random crack networks is presented for two-dimensional (2D) cases; then the IDD model and percolation scaling law for overall permeability are established; after that, numerical analysis and validation are provided; then, further analysis of crack opening and experimental validation are given; and the concluding remarks are summarized finally.Geometry of Random Crack NetworksThe geometry of random crack networks includes single cracks and crack clusters (Zhou et al. 2012; Li and Li 2015). The single crack geometry refers to the length, opening aperture, and tortuosity of the crack surface, and the network geometry includes the crack extent (density), orientation, and the crack cluster size (Adler et al. 2013). These geometry parameters have different influences on the permeability evaluation (Akhavan et al. 2012; Li and Li 2018). Three local parameters, namely opening aperture, clustering extent, and connectivity degree, and two averaging parameters, namely crack density and orientation, are retained to describe the crack networks. The geometry parameters are defined as follows and they are applied to illustrate impact of the crack network on the permeability.Crack DensityCrack density describes the damage extent of porous solids, defined as follows: (1) where li = ith crack length; N = total crack number; and As = cracked area. This area-based density definition is widely used in micromechanical models (Kachanov 1992) and also applied in effective properties estimate of cracked solids (Hasselman 1978). Fig. 1 illustrates three random crack networks with crack densities ρ=0.01, 0.1, 1.0, respectively.OrientationSupposing N cracks with length li(i=1,2,…,N) in a finite domain, the orientation degree ω can be determined through an orientation factor ι(2) ω=ι(θ)minι(θ)max,withι(θ)=1As∑i=1N(li2)2cos2θiwhere θi = angle between crack li and the flow direction; and ι(θ)min and ι(θ)max = minimum and maximum values of the orientation factor. The orientation degree ω=1 stands for a perfectly isotropic crack network, (Fig. 1) and ω=0 as all cracks are parallel to the transport direction.Fractal DimensionThe fractal dimension is a measure of the spatial organization of the cracks, which is independent on the crack length distribution and limited by the number of topological dimensions of the domain (three for a volume and two for a surface) (Davy et al. 1990). A two-point correlation function was applied to account for the fractal characteristics (Darcel et al. 2003). The correlation function C(r) scales with r as follows: (3) where NT = total crack number; N(r) = number of cracks for which the distance between crack centers is within r; and Dc = fractal dimension that can be regressed by C(r)−r on a logarithm scale. Fig. 2 illustrates three random crack networks with the same crack density ρ=0.5 and different fractal dimensions Dc=1.99, 1.75, 1.62.ConnectivityThe connectivity measures the cracking extent to which the flow paths are formed in random crack networks, a critical parameter for permeability estimation. A connectivity factor, f, based on percolation concepts of Li and Li (2015) considering local clustering, was used here (4) f=ξLP(3−Dc),withf∈[0,1]where ξ = cluster correlation length; L = domain size; Dc = fractal dimension; and P = likelihood of a crack being connected to the spanning clusters. If all cracks connect with each other and form into one spanning cluster across the domain, i.e., ξ=L and P=1, the connectivity factor f achieves the maximum value of 1. The connectivity of crack networks in Fig. 2 is 5×10−5 [Fig. 2(a)], 0.022 [Fig. 2(b)], and 0.089 [Fig. 2(c)]. When L→∞, the connectivity f near percolation threshold obeys the following scaling law: (5) f∝|ρ−ρc|−ν+(3−Dc)βwhere the percolation threshold, ρc, depends on the geometry characterization of crack clusters; beyond the threshold, a percolating crack flow path appears in statistical sense. The exponents β and ν are universal constants related to the topological dimension of the domain, i.e., β=5/36 and ν=4/3 for the two-dimensional case and β=0.41 and ν=0.875 for the three-dimensional case (Aharony and Stauffer 2003).Permeability of Cracked SolidsIDD Model for Nonoverlapping CracksFor the matrix-inclusion morphology, the crack inclusions are separately distributed in a continuous matrix phase. Permeability is taken as the linear physical property to be studied. So the permeability fluctuation Hi of crack inclusion i can be expressed (6) where the permeability tensors Ki,m stand for crack inclusion and matrix phase, respectively.The dilute estimate for the permeability increment Hid of an ellipsoidal matrix containing one ellipsoidal crack inclusion can be written (Shafiro and Kachanov 2000; Zhou et al. 2011) (7) Hid=ci(Hi−1+Ωi0)−1=ci(Hi−1+Km−1J)−1where ci = volumetric fraction of the ith crack inclusion; Ωi0 is the eigenstiffness tensor characterizing the ith inclusion shape and orientation; superscript 0 = matrix; susperscript d = dilute solution; and J is the geometry tensor related to the size of ellipsoidal inclusion. Its components Jj=1,2,3 are written (8) J1=a1a2a32∫0∞du(a12+u)(a12+u)(a22+u)(a32+u)with J2 and J3 obtained by cyclic permutations of all indices 1, 2, and 3. The terms a1,2,3 are semiaxes of the inclusion parallel to coordinate axes x1,2,3, and e1,2,3 are unit vectors along them. For some special cases, e.g., two-dimensional analysis, Jj can be simplified as follows: (9) J1=0,J2=a3a2+a3,J3=a2a2+a3The IDD estimate derives from the effective self-consistent (ESC) scheme; the root of them is the inclusion-matrix-composite model. The IDD estimate considers the inclusion spatial distributions as well as the interaction between inclusions and the immediate surrounding matrix (Zheng and Du 2001). This estimate has an explicit and simple structure for various inclusion shapes, inclusion volume fraction, as well as inclusion and matrix anisotropies. Explicit solutions with second-order accuracy can be provided by the IDD estimate in a statistical sense. The IDD scheme gives the permeability increment Hidd of a solid containing N groups of crack inclusions as follows (Zheng and Du 2001; Du and Zheng 2002): (10) Hidd=(I−∑i=1NHid·ΩDi0)−1·∑i=1NHidwhere ΩDi0 is the eigenstiffness tensor of the matrix atmosphere for inclusion i. The mathematical formulations for the classical IDD model are provided in the Appendix, and in the following, we develop further this model to account for crack density and opening in permeability estimation.The corresponding matrix atmosphere was selected as ellipsoidal for the three-dimensional case and elliptical for the two-dimensional case, with their principal axes aligned on the crack plane (Fig. 3). Considering the plane e2e3 and letting a1→∞ for two-dimensional analysis, the terms 2a20 and 2a30 denote the axis lengths of the elliptical matrix atmosphere, 2a2 is the crack length, 2a3 is the crack opening aperture, and a20=a2+a30. The crack density ρ in solid can be evaluated as follows: (11) where N = number of cracks (inclusions); and Aa = domain area. To maintain the same density in the crack-atmosphere cell equal to that in the whole domain, the geometrical parameters can be evaluated as follows: (12) a22π(a2+a30)a30=ρ⇒a30a20=a30a2+a30=1+4πρ−11+4πρ+1Using the dilute estimate in Eq. (7), the permeability increment Hid in the plane e2e3 can be calculated as follows: (13) Hid=ci[Hi−1+Km−1(a3a2+a3e2e2+a2a2+a3e3e3)i]−1By using the definition in Eq. (6), Hid can be rearranged as follows: (14) Hid=πa2a3Aa[(Ki−Km)·Km(a2+a3)Kma2+Kia3e2e2+(Ki−Km)·Km(a2+a3)Kma3+Kia2e3e3]iThe two-dimensional (2D) eigenstiffness tensor of the elliptical matrix atmosphere is (15) ΩDi0=Km−1·(a30a20+a30e2e2+a20a20+a30e3e3)iCombining Eqs. (14) and (15) with Eq. (10), the IDD estimate of permeability can be obtained.Now consider a crack with opening aperture w′ embedded in a solid matrix of low permeability. The classical intrinsic permeability of the crack between two parallel impermeable surfaces is expressed as Ki=w′2/12 (Timothy and Meschke 2017b). Serving as a reference, the intrinsic permeability of a cement paste (matrix) is around Km=10−17  m2 (Picandet et al. 2009) with crack permeability Ki=10−12  m2 with opening w′=5  μm, the lower limit for such cracks. Thus, it is rational to assume that Ki≫Km for physically meaningful cracks in engineering porous materials like concretes and rocks (Li et al. 2011; Dormieux et al. 2011). Using this assumption and considering a single crack family characterized by the same aspect ratio w=a3/a2 (Lemarchand et al. 2005), one can derive the relative permeability along the crack’s orientation for solids with parallelly arranged microcracks as follows: (16) KxKm=1+πρ(1+w)(1−γD)1−πρ(1+w)γDand the effective permeability for solids with randomly oriented microcracks as follows: (17) KeffKm=1+π2ρ(1+w)[w+(1+w)(1−γD)]1−π2ρ(1+w)[−w+(1+w)γD]where γD = geometrical parameter of an elliptical matrix atmosphere, expressed by (18) γD=a30a20+a30=1+4πρ−121+4πρTo investigate the influence of opening aperture w and crack density ρ in Eq. (17), the effective permeability Keff/Km was calculated for various crack opening aperture ratios (w=0.001 to 0.03 and crack densities (ρ=0 to 1) for nonoverlapping cracks (Fig. 4). The crack density contributes to the effective overall permeability in a monotone manner, and the crack opening does not change notably the effective permeability. It is rather reasonable according to the assumptions of the IDD estimate: the cracks with very high local permeability at large w are assumed positioned in a low permeable matrix, i.e., the cracks are separated and the matrix is continuous, serving as barriers for the flow.Extended IDD Model for Low and Medium Crack ClusteringFor crack density far below the percolation threshold, there is not yet a backbone cluster of cracks dominating the permeability, but local clustering can occur to increase the flow perturbance of flow. In this section, the classical IDD model is extended to account for this local crack clustering. Fig. 5 shows a random network with constant length l, crack density ρ=0.2, and domain size ratio L/l=5. Assume that the projected length of the largest cluster on the transport direction equals to ξ. Supposing the crack network contains Nξ clusters whose sizes are in the interval between 0.5ξ and ξ, the projected length of the ith cluster on the transport direction is denoted (19) ξj=ψjξ,ψj∈[0.5,1.0]As two cracks interconnect with each other, it seems reasonable to arrange a parallel crack with length ξj on the transport direction to evaluate the connectivity impact on the effective permeability (Fig. 5). The equivalent density ρh is (20) ρh=∑j=1Nξξj24A=ρ(∑j=1Nξψj2·ξ2)/(4·∑i=1N(a2)i2)which can further consider the connectivity definition in Eq. (4) and be rearranged as follows: (21) ρh=∑j=1Nξψj2[4f·P−(3−Dc)]2As shown in Fig. 5, a two-step scheme can be adopted to extend the IDD model to account for the local clustering. First, nonoverlapping cracks are introduced into the matrix so that one can evaluate the overall permeability using classical IDD model. The effective permeability K1 obtained at this step is attributed to the host matrix for the second step. In this step, a parallel crack is added on transport direction at the same position to consider the amplification effect of flow by the local crack connectivity (clustering). The relative permeability in this step can be estimated by Eq. (16) as follows: (22) KxK1=1+πρh(1+w)(1−γDx)1−πρh(1+w)γDxIn total, the classical IDD model for effective permeability can be extended to crack networks with connected cracks as follows: (23) KeffKm=1+πρh(1+w)(1−γDx)1−πρh(1+w)γDx·1+π2ρ(1+w)[w+(1+w)(1−γD)]1−π2ρ(1+w)[−w+(1+w)γD]where ρh = equivalent density related to the connectivity f; and coefficients γDx and γD can be derived from crack densities ρh and ρ by Eq. (18). The quantification of ρh will be given in the “Numerical Analysis” section.Percolation Model for High Crack ClusteringAs clustering effects become dominant and the spanning cluster of cracks is about to occur, the permeability of a cracked solid is to increase significantly, and the primary controlling parameter in the permeability of a crack network is connectivity. The effective permeability of cracked solids can be expressed (24) where Km = permeability of the homogeneous porous matrix; f = percolation-based connectivity in Eq. (5); τ = tortuosity of the backbone cluster; and ρ and w = crack density and crack opening aperture ratio, respectively. The continuum percolation theory allows the following scaling law for effective permeability: (25) Keff=K0(Km,w)|ρ−ρc|μwithμ=(3−Dc)β−(2−Dτ)νwhere exponent μ = connectivity and the tortuosity of dominating cluster; Dτ = fractal dimensions of the dominating cluster; and the term K0 = impact of crack opening and matrix permeability. Some further adaptations for a finite flow domain have been discussed by Li and Li (2015).Model SummaryThe permeability of porous solids incorporating random crack networks is expressed through two approaches: the IDD estimate and continuum percolation theory (Table 1). The established model covers three characteristic ranges of cracks: (1) for nonoverlapping cracks or cracks with low perturbance for fluid flow, the classical IDD model is employed, (2) for low and medium crack clustering, the extended IDD model is derived, and (3) for high clustering of cracks, the percolation scaling law is employed. The established model is clustering-sensitive and will be validated through numerical simulations for different clustering cases.Table 1. Effective permeability models applied for the whole range of crack density and clusteringTable 1. Effective permeability models applied for the whole range of crack density and clusteringCrack clustering extentρ<ρcρ→ρcρ>ρcNo clustering (nonoverlapped)IDD modelIDD model—Low and medium clusteringExtended IDD model——High clustering—Percolation scaling law—Percolated crack networks——Percolated flowNumerical AnalysisNumerical Generation of Random Crack NetworksA finite domain L×L with random distributed cracks is established. The crack length follows a lognormal distribution from the geometry characterization on cracks of concrete damage under compressive loading (Zhou et al. 2012) (Fig. 6). The statistical analysis showed that the cracking opening aperture has a clear correlation with crack length. So the ratio of opening aperture to length, w, was taken as a basic variable ranging from 0.001 to 0.028 from the results in Fig. 6, when the mean crack length lmean equals to 3.606 mm and the opening aperture ranges from 5 to 100 μm (Li and Li 2015). A length scale χ=5, 10, 15, 20, 25, 30 was retained as the ratio between the finite domain size L and the mean crack length lmean. The constant length distribution serves as comparative case in numerical analysis.A Monte-Carlo algorithm was designed to generate different crack networks. The crack density, crack opening aspect ratio, crack length distribution, and the average number of interaction per crack were preset before the simulation. According to the preset statistical parameters, one can generate random cracks and locate them one by one into the finite area. Once a random crack is placed, the total intersection number of all cracks is calculated; if the average interaction for each crack does not surpass the preset parameter, this crack will be adopted; otherwise, it will be deleted. Meanwhile, the isolated cracks and connected crack clusters are characterized and recorded on this step. These simulation steps are repeated until the preset parameters, i.e., crack density and average intersection number, are accomplished.The geometry characterization for the isolated cracks and crack clusters in the generated crack network, including the fractal dimension Dc and the connectivity factor f, was derived. Controlling the intersection ratio per crack produces correlated networks with different fractal dimensions Dc=1.43 to 1.99; releasing the control of the intersection ratio results in an uncorrelated network with Dc=2.0. Three correlated networks with the same density and domain size (ρ=0.5 and χ=20) are generated in Fig. 2. The intersection ratio per crack was controlled as 0.0 (Dc=1.99), 2.0 (Dc=1.75), and 3.0 (Dc=1.62), respectively. Fig. 7 shows two uncorrelated networks with the same crack density ρ=0.5, domain size χ=20, fractal dimension Dc=2.0, and different length distributions.Permeability SimulationThe generated domain containing crack networks was retained for permeability analysis. The assumptions of incompressible and Newtonian fluid were adopted in the permeability simulation. The effective permeability Keff can be calculated through Darcy’s law, Keff=−ηq/∇p, where q is the flow rate, η is the fluid viscosity, and ∇p is the pressure gradient. In the numerical solution, along the flow direction, a pressure gradient was applied and perpendicular to the flow direction, null flux conditions were exerted. The finite-element method (FEM) was used to solve the pressure field p(x,y): triangle elements were applied to mesh the matrix; line elements were prescribed to simulate the cracks; and the two nodes of line element as well as the intersection point of cracks were considered as shared nodes.When the intrinsic permeability of the porous matrix takes Km=10−17  m2, the local crack permeability Kcrack is eight magnitudes higher than the matrix through Kcrack=w′2/12 with crack opening w′ of 50 μm.Model ValidationIDD Model for Nonoverlapping CracksCrack networks of size χ=20 were constructed with randomly oriented and located cracks, which followed a constant length distribution. There was no intersection between generated cracks, shown in Fig. 2(a). The crack opening ratio w=0.0139 and the crack density ranged between 0 and 1 with an interval of 0.1. For each crack density level, 40 crack networks were generated, and the corresponding effective permeability Keff was solved by FEM. The average values of Keff are given in Fig. 8(a) together with the IDD estimate, dilute scheme, the Mori-Tanaka (MT) scheme (Timothy and Meschke 2018), and differential scheme (DS).The IDD estimate showed good agreement with the numerical results for the whole range of crack densities, whereas the MT scheme and the dilute scheme lost their prediction capacity at crack densities higher than 0.15, and the DS overestimated the effective permeability at crack densities higher than 0.35. Fig. 8(b) presents the effective permeability in terms of the crack opening ratio w ranging from 0.0014 to 0.0277 with ρ=1.0 for both numerical simulations and IDD estimates. One can see that the overall permeability shows very little change for the simulated values of crack opening.Extended IDD Model for Low and Medium Crack ClusteringTo illustrate the extended IDD model and obtain the equivalent horizontal density ρh, numerical simulations were performed for two different crack patterns with constant length, Dc=2.0, 1.75, shown in Figs. 2(b) and 7(a). The study domain size χ=20, and crack opening ratio w was 0.0139. The crack density ranges from 0 to 1.4 (1.6 for Dc=1.75), with an interval of 0.1, and 40 crack networks were generated for each crack density level. The equivalent horizontal density ρh can be obtained for each crack network. Figs. 9(a) and 10(a) show the average value of ρh in terms of crack density ρ for uncorrelated (Dc=2.0) and correlated crack networks (Dc=1.75), respectively.From Fig. 9(a), it can be seen that the equivalent horizontal density ρh increased slowly with ρ until ρ=0.4 but almost linearly afterward. That is because the connectivity factor f is ρ-dependent and the f values are very small at a low crack density, e.g., f=0.002 at ρ=0.2.The relative permeability was then evaluated by extended IDD model through Eq. (23), shown in Fig. 9(b). For ρ<0.4, the relative permeability obtained by IDD estimation for nonoverlapping cracks also showed good agreement with numerical results. From Fig. 10(a), the equivalent horizontal density ρh increased almost linearly with ρ due to the generation procedure of correlated networks. By using the obtained ρh, the effective permeability was evaluated through Eq. (23) and is illustrated in Fig. 10(b).Within the performed simulation ranges, namely ρ<1.1, ρh<0.335, and f<0.067 for Dc=2.0 and ρ<1.3, ρh<0.356, and f<0.168 for Dc=1.75, the extended IDD estimate was validated and showed good accuracy for crack networks with low and medium clustering. For higher crack density and notable crack clustering extents, the IDD estimate diverged from the numerical results because the clustering effect becomes dominating, and this effect is beyond the prediction capacity of the extended IDD model.Percolation Model for High Crack ClusteringThe high crack clustering effect on permeability was quantified through the percolation scaling law. In a forgoing study (Li and Li 2015), for uncorrelated networks, Dc=2.0, the percolation thresholds ρc were 1.43 and 1.67 for constant and lognormal length distributions. For correlated networks, ρc=1.68 and 1.78 for constant lengths with Dc=1.75 and 1.62, respectively. As crack density ρ→ρc, the scaling law in Eq. (25) is valid, i.e., |ρ−ρc|≤0.33 or 0.38 for uncorrelated and correlated networks. The scaling exponent μ of permeability is calculated as −0.922 (constant length) and −0.919 (lognormal length) for Dc=2.0, −1.122 (constant length) for Dc=1.75, and −1.257 (constant length) for Dc=1.62 by linear regression between the effective permeability and |ρ−ρc| on a logarithm scale. The regressed parameters are presented in Table 2.Table 2. Scaling parameters of permeability for different DcTable 2. Scaling parameters of permeability for different DcCrack length distributionDcυλμConstant2.0011.570.199−0.922Constant1.7556.260.477−1.122Constant1.6292.850.554−1.257Lognormal2.0013.790.223−0.919The effective permeability from the percolation scaling law was compared with results from the cascade continuum micromechanics (CCM) model proposed by Timothy and Meschke (2017a). The CCM model predicts the effective permeability through the Eshelby mean-field homogenization method. At the self-similar limit level [described as n→∞ by Timothy and Meschke (2017a)], the CCM model degenerates to the self-consistent scheme. For low matrix permeability, Km=10−18  m2, local crack permeability Kcrack is eight magnitudes higher than Km. The CCM model displays a sharp phase transition as the crack density increases. The order of magnitude of the permeability increases beyond a critical threshold of the crack density parameter, depending on the value of Kcrack/Km (Dormieux and Kondo 2004).An interesting connection between CCM model and the percolation model can be found. Although these two models were established based on two totally different assumptions, the effective permeability (Keff/Km) of uncorrelated crack networks with random distributed cracks, derived by the percolation scaling law in this study, is in good agreement with the effective permeability calculated by the CCM model (short-range interaction).The CCM model can predict the effect of the crack opening aspect ratio on the features of the conversion from a nonconnected solid to a completely connected solid, i.e., the percolation threshold depends on the crack opening aspect ratio (Leonhart et al. 2015, 2017). The percolation scaling law deals with the crack clustering effect and takes the cracks no longer as homogeneous inclusions but as networks, i.e., the percolation threshold relies on the local crack clustering effect and the connectivity factor.Further AnalysisImpact of Crack OpeningFig. 11 presents the effective permeability evolution Keff as the crack density ranges between [0,ρc] for cgs (w=0.001 to 0.028) for uncorrelated crack networks (Dc=2.0) with domain size χ=20. The influence of opening aperture w on Keff is less important as ρ≪ρc and becomes remarkable as ρ→ρc. Approaching the percolation threshold, a flow path is almost formed, and the flow rate is directly influenced by the crack opening aperture.To obtain the impact of crack opening w on Keff in Eq. (25) as ρ→ρc, the interested density range is ρ∈[ρc−b,ρc), with b=0.33 or 0.38 for uncorrelated and correlated crack networks, respectively. Figs. 12(a–d) illustrate the effective permeability Keff in terms of crack openings w for correlated/uncorrelated networks with constant and lognormal crack length distributions. The effective permeability Keff always increases with w=0.001 to 0.004 at a fixed crack density but keeps rather steady as w scales from 0.007 to 0.028. It is reasonable because for large w values, the local crack permeability becomes very high.Figs. 12(a–h) present the logarithmic plots of K0/Km to w for various ρ and Dc. On the logarithmic scale K0/Km is a linear function of w for the range 0.001≤w≤0.004, i.e., as follows: (26) K0∝wλfor  0.001≤w≤0.004The λ exponents are, respectively, 0.199, 0.477, 0.554, and 0.223 for Dc=2.0 (constant length), Dc=1.75 (constant length), Dc=1.62 (constant length), and Dc=2.0 (lognormal length). The exponent λ versus Dc is given in the following equation: (27) Substituting Eq. (26) into Eq. (25) gives (28) KeffKm=K0Km|ρ−ρc|μ={υw−0.937Dc+2.09|ρ−ρc|μ0.001≤w≤0.0044|ρ−ρc|μw>0.004where υ = fitting parameter, given in Table 2. Due to the amplification of local flow caused by clustering effects, υ increases for correlated crack networks (Dc<2.0). The expression in Eq. (28) reveals that Keff depends on the crack opening w via a power law through the term K0.With the specific parameters retained, the permeability of solids containing percolated crack networks presents an increase of several orders of magnitude compared with the unpercolated cases. When cracks intersect to form a percolated path across the boundaries of the domain, the effective permeability will be significantly altered (Leonhart et al. 2017; Renshaw et al. 2020; Benkemoun et al. 2018). After onset of percolation, the effective permeability is generally stable [Fig. 13(a)], and more likely to be determined by the total length of the percolated path. Although the permeability becomes higher when the crack density is larger, the overall impact of crack density on the permeability of percolated networks is limited. However, for the percolated crack networks, the effective permeability Keff exhibits a power law with respect to the crack opening w. In Fig. 13(b), the exponent of power law is regressed as 2.968, which is in good agreement with the results reported in literature (Li et al. 2019; Reinhardt and Jooss 2003; Picandet et al. 2009).Experimental DataThe correlation between a two-dimensional (2D) crack network and the effective permeability of cement-based materials was investigated by Li and Li (2018). The actual crack networks were extracted, and the geometry parameters including crack density, crack orientation, connectivity factor, and crack opening aperture were evaluated. The effective permeability was obtained indirectly using moisture transport modeling and drying tests. The cracked specimens had density in the range of 0.044–0.399 and relative permeability Keff/Km of 1–4.3. The experimental values Keff/Km and ρ were employed to compare with the permeability evolution derived from the model in Fig. 14. From the extracted actual crack geometry, i.e., ρ<0.4, the established model was vibrated against the experimental results. The missing clear comparison for ρ>0.4 is due the limited crack density extent.Appendix. Basis for IDD SchemeThis section is derived from the following studies: Zheng and Du (2001), Du and Zheng (2002), Shafiro and Kachanov (2000), and Zhou et al. (2011). The linear permeation problem in cracked solids can be expressed (29) where F = liquid flux; P = pressure gradient; and K is the permeability tensor. In the following, the permeability tensors Km, Ki, and K stand for the matrix, the ith inclusion, and the effective medium, respectively. The tensors Hi and H refer to the permeability increment of the ith inclusion and effective medium compared with the matrix, described as follows: (30) For the inclusion-matrix morphology, the IDD scheme takes into account the interaction between inclusions and their surrounding matrix, giving explicit solutions for linear physical properties (Zheng and Du 2001). The dilute estimation of the permeability increment Hd for a solid containing N groups of inclusions is expressed (31) Hd=∑i=1NHid=∑i=1Nci(Hi−1+Ωi0)−1where Hid and ci = dilute estimation and volumetric fraction of inclusion i, respectively; Ωi0 is the eigenstiffness tensor, which characterizes the shape and orientation of inclusion i; superscript 0 = dependency on the matrix; and d = dilute.By considering the interaction between inclusions, the permeability H obtained through IDD method can be written (Du and Zheng 2002) (32) Hidd=(I−∑i=1NHid·ΩDi0)−1Hdwhere ΩDi0 is the eigenstiffness tensor of matrix atmosphere; and Hidd is directly calculated by Hid, ΩDi0, and Hd. These quantities all have definite physical meanings and can be quantitatively expressed. If all inclusion and matrix atmospheres have same shape and orientation, ΩDi0=ΩD0, and Eq. (32) degenerates to (33) Hidd=((Hd)−1−ΩD0)−1Following the development of Shafiro and Kachanov (2000) and Zhou et al. (2011), considering an ellipsoidal matrix atmosphere with permeability Km containing one ellipsoidal inclusion, the effective permeability estimated by dilute scheme is described (34) Hd=∑i=1NHid,Hid=ci(Hi−1+Km−1J)−1where J is the geometry tensor related to the size of ellipsoidal inclusion. Assuming that a1,2,3 are semi axes of the inclusion and a1>a2>a3, the components Jj=1,2,3 are given (35) J1=a1a2a3(a12−a22)a12−a32[F(θ,κ)−E(θ,κ)]J3=a1a2a3(a22−a32)a12−a32[a2a12−a32a1a3−E(θ,κ)]J2=1−J1−J3with (36) E(θ,κ)=∫0θ1−κ2sin2wdwF(θ,κ)=∫0θ1/1−κ2sin2wdwand (37) θ=arcsin1−(a3/a1)2,κ=a12−a22a12−a32By equating Eqs. (31) and (34), the eigenstiffness tensor of inclusion i can be written (38) The eigenstiffness tensor of matrix atmosphere ΩDi0 can also be obtained by Eq. (38) with the known geometry parameters.References Abdalrahman, T., S. Scheiner, and C. Hellmich. 2015. “Is trabecular bone permeability governed by molecular ordering-induced fluid viscosity gain? Arguments from re-evaluation of experimental data in the framework of homogenization theory.” J. Theor. Biol. 365 (Jan): 433–444. Adler, P. M., and J. F. Thovert. 1999. Fracture and fracture networks. Dordrecht, Netherlands: Kluwer. Adler, P. M., J. F. Thovert, and V. V. Mourzenko. 2013. Fractured porous media. London: Oxford University Press. Aharony, A., and D. Stauffer. 2003. Introduction to percolation theory. New York: Taylor & Francis. Akhavan, A., S. M. H. Shafaatian, and F. Rajabipour. 2012. “Quantifying the effects of crack width, tortuosity, and roughness on water permeability of cracked mortars.” Cem. Concr. Res. 42 (2): 313–320. Benkemoun, N., H. Alkhazraji, P. Poullain, M. Choinska, and A. Khelidj. 2018. “3-D mesoscale simulation of crack-permeability coupling in the Brazilian splitting test.” Int. J. Numer. Anal. Methods Geomech. 42 (3): 449–468. Damrongwiriyanupap, N., S. Scheiner, B. Pichler, and C. Hellmich. 2017. “Self-consistent channel approach for upscaling chloride diffusivity in cement pastes.” Transp. Porous Media 118 (3): 495–518. Darcel, C., O. Bour, P. Davy, and J. R. de Dreuzy. 2003. “Connectivity properties of two-dimensional fracture networks with stochastic fractal correlation.” Water Resour. Res. 39 (10): 1272. Davy, P., A. Sornette, and D. Sornette. 1990. “Some consequences of a proposed fractal nature of continental faulting.” Nature 348 (6296): 56–58. Dormieux, L., and D. Kondo. 2004. “Coupling between permeability and damage: A micromechanical approach.” In Proc., XXI Int. Congress of Theoretical and Applied Mechanics. New York: Springer. Dormieux, L., and D. Kondo. 2005. “Diffusive transport in disordered media: Application to the determination of the tortuosity and the permeability of cracked materials.” In Applied micromechanics of porous materials, edited by L. Dormieux and F.-J. Ulm, 83–106. New York: Springer. Du, D. X., and Q. S. Zheng. 2002. “A further exploration of the interaction direct derivative (IDD) estimate for the effective properties of multiphase composites taking into account inclusion distribution.” Acta Mech. 157 (1): 61–80. Eshelby, J. D. 1957. “The determination of the elastic field of an ellipsoidal inclusion, and related problems.” Proc. R. Soc. London, Ser. A 241 (1226): 376–396. Estermann, S., and S. Scheiner. 2018. “Multiscale modeling provides differentiated insights to fluid flow-driven stimulation of bone cellular activities.” Front. Phys. 6 (Sep): 76. Hu, X., X. Song, K. Wu, Y. Shi, and G. Li. 2021. “Effect of proppant treatment on heat extraction performance in enhanced geothermal system.” J. Pet. Sci. Eng. 207 (Dec): 109094. Hunt, A., and R. Ewing. 2009. Percolation theory for flow in porous media. Cham, Switzerland: Springer. Kachanov, M. 1992. “Effective elastic properties of cracked solids: Critical review of some basic concepts.” Appl. Mech. Rev. 45 (8): 304–335. Kushch, V. 2020. Micromechanics of composites: Multipole expansion approach. Oxford, UK: Butterworth-Heinemann. Langlois, V., V. H. Trinh, C. Lusso, C. Perrot, X. Chateau, Y. Khidas, and O. Pitois. 2018. “Permeability of solid foam: Effect of pore connections.” Phys. Rev. E 97 (5): 053111. Lemarchand, E., C. A. Davy, L. Dormieux, and F. Skoczylas. 2010. “Tortuosity effects in coupled advective transport and mechanical properties of fractured geomaterials.” Transp. Porous Media 84 (1): 1–19. Lemarchand, E., L. Dormieux, and F.-J. Ulm. 2005. “Micromechanics investigation of expansive reactions in chemoelastic concrete.” Philos. Trans. R. Soc. London, Ser. A 363 (1836): 2581. Leonhart, D., J. J. Timothy, and G. Meschke. 2015. “Determination of effective transport properties of fractured rocks using the extended finite element method and micromechanics.” In Proc., ISRM Regional Symp.-EUROCK 2015. Austin, TX: OnePetro. Leonhart, D., J. J. Timothy, and G. Meschke. 2017. “Cascade continuum micromechanics model for the effective permeability of solids with distributed microcracks: Comparison with numerical homogenization.” Mech. Mater. 115 (Dec): 64–75. Li, K. 2016. Durability design of concrete structure: Phenomena, modelling and practice. New York: Wiley. Li, L., and K. Li. 2015. “Permeability of microcracked solids with random crack networks: Role of connectivity and opening aperture.” Transp. Porous Media 109 (1): 217–237. Li, L., and K. Li. 2018. “Experimental investigation on transport properties of cement-based materials incorporating 2d crack networks.” Transp. Porous Media 122 (3): 647–671. Li, L., W. Liu, Q. You, M. Chen, and Q. Zeng. 2020. “Waste ceramic powder as a pozzolanic supplementary filler of cement for developing sustainable building materials.” J. Cleaner Prod. 259 (Jun): 120853. Li, M., Y. B. Tang, Y. Bernabé, J. Z. Zhao, X. F. Li, and T. Li. 2017. “Percolation connectivity, pore size, and gas apparent permeability: Network simulations and comparison to experimental data.” J. Geophys. Res.: Solid Earth 122 (7): 4918–4930. Liu, J., J. Huang, K. Liu, and K. Regenauer-Lieb. 2021. “Identification, segregation, and characterization of individual cracks in three dimensions.” Int. J. Rock Mech. Min. Sci. 138 (Feb): 104615. Mustapha, H., K. Makromallis, and A. Cominelli. 2020. “An efficient hybrid-grid crossflow equilibrium model for field-scale fractured reservoir simulation.” Comput. Geosci. 24 (2): 477–492. Nguyen, T. T. N., M. N. Vu, N. H. Tran, N. H. Dao, and D. T. Pham. 2020. “Stress induced permeability changes in brittle fractured porous rock.” Int. J. Rock Mech. Min. Sci. 127 (Mar): 104224. Renshaw, C. E., E. M. Schulson, D. Iliescu, and A. Murzda. 2020. “Increased fractured rock permeability after percolation despite limited crack growth.” J. Geophys. Res.: Solid Earth 125 (8): e2019JB019240. Shafiro, B., and M. Kachanov. 2000. “Anisotropic effective conductivity of materials with nonrandomly oriented inclusions of diverse ellipsoidal shapes.” J. Appl. Phys. 87 (12): 8561–8569. Timothy, J. J. 2016. “Analytical and computational models for the effective properties of disordered microcracked porous materials.” Ph.D. thesis, Faculty of Civil and Environmental Engineering, Ruhr-Univ. Bochum. Timothy, J. J., and G. Meschke. 2017a. “Cascade continuum micromechanics model for the effective permeability of solids with distributed microcracks: Self-similar mean-field homogenization and image analysis.” Mech. Mater. 104 (Jan): 60–72. Timothy, J. J., and G. Meschke. 2017b. “The intrinsic permeability of microcracks in porous solids: Analytical models and Lattice Boltzmann simulations.” Int. J. Numer. Anal. Methods Geomech. 41 (8): 1138–1154. Timothy, J. J., and G. Meschke. 2018. “Effective diffusivity of porous materials with microcracks: Self-similar mean-field homogenization and pixel finite element simulations.” Transp. Porous Media 125 (3): 413–434. Torrijos, M. C., G. Giaccio, and R. Zerbino. 2010. “Internal cracking and transport properties in damaged concretes.” Supplement, Mater. Struct. 43 (S1): 109–121. Vinci, C., J. Renner, and H. Steeb. 2014. “A hybrid-dimensional approach for an efficient numerical modeling of the hydromechanics of fractures.” Water Resour. Res. 50 (2): 1616–1635. Wang, Y., L. Jeannin, F. Agostini, L. Dormieux, and E. Portier. 2018. “Experimental study and micromechanical interpretation of the poroelastic behaviour and permeability of a tight sandstone.” Int. J. Rock Mech. Min. Sci. 103 (Mar): 89–95. Zheng, Q. S., and D. X. Du. 2001. “An explicit and universally applicable estimate for the effective properties of multiphase composites which accounts for inclusion distribution.” J. Mech. Phys. Solids 49 (11): 2765–2788.

Source link

Leave a Reply

Your email address will not be published. Required fields are marked *