Three-Dimensional Network Model for Coupling of Fracture and Mass Transport in Quasi-Brittle Geomaterials

Dual three-dimensional networks of structural and transport elements were combined to model the effect of fracture on mass transport in quasi-brittle geomaterials. Element connectivity of the structural network, representing elasticity and fracture, was defined by the Delaunay tessellation of a random set of points. The connectivity of transport elements within the transport network was defined by the Voronoi tessellation of the same set of points. A new discretisation strategy for domain boundaries was developed to apply boundary conditions for the coupled analyses. The properties of transport elements were chosen to evolve with the crack opening values of neighbouring structural elements. Through benchmark comparisons involving non-stationary transport and fracture, the proposed dual network approach was shown to be objective with respect to element size and orientation.


Introduction
The influence of fracture and mass transport affects physical processes that govern many engineering applications, such as deterioration of construction materials and performance of waste barriers. In these applications, it is important to consider the influence of fracture induced pathways for the ingress of fluids. Modelling the discrete crack formations, and the mass transport along these cracks and through the surrounding uncracked material, is challenging. Models for this coupling are commonly based on continuum mechanics combined with a discrete representation of cracks [1][2][3][4][5]. Alternatively, discrete approaches, such as discrete element method, lattice and network models, have been proposed to model these processes [6][7][8][9][10][11][12][13][14][15].
One network approach, based on the Delaunay tessellation of a random set of points, has been shown to be suitable for modelling fracture [16][17][18][19] and mass transport [20], providing mesh insensitive results. In this approach, the physical processes are modelled by a multi-dimensional network of one-dimensional elements, which are placed on the Delaunay edges ( Figure 1a); the element properties are determined by the corresponding Voronoi tessellation. The nodes of the elements of structural and transport network models coincide, which is suitable for modelling the coupling of continuum fields. However, once cracks are formed, the transport elements in this approach are orientated perpendicular to the crack path, which is aligned with the mid-cross-section of the structural elements ( Figure 1b). This misalignment of the transport elements with the crack path complicates the modelling of crack-assisted transport and its dependence on crack opening. To resolve this deficiency, several researchers [8,10,13,21] have placed transport elements on the Voronoi edges, whereas the structural elements remain on the Delaunay edges ( Figure 1c). With this dual network approach, the influence of fracture on transport is more naturally represented, since the transport elements are aligned with the potential crack directions. So far, most of this work was either limited to 2D or did not provide coupling between fracture and transport.  This work proposes a three-dimensional dual network approach for modelling fracture and mass transport. Structural elements are placed on the edges of Delaunay tetrahedra and transport elements are placed on the edges of Voronoi polyhedra. Special attention is given to the discretisation of the dual networks at domain boundaries. Simple geometric relationships based on Voronoi and Delaunay tessellations are proposed for describing the change of permeability as a function of crack opening. By a series of benchmarks, it has been demonstrated that the present approach can describe fracture, transport, and the increase of permeability due to fracture mesh-independently. Fracture is modelled by means of a cohesive-frictional approach, which is suitable for geomaterials, such as concrete and rocks in which the size of the fracture process zone is large compared to the size of the structure. Transport is modelled by means of Darcy's flow equation. The proposed model is designed to describe the effect of cohesive fracture on conductivity.

Network Approach
The new network approach uses one-dimensional elements connected in a three-dimensional network to describe continuum fields as well as evolutions of discontinuities in the form of fracture process zones. In the present section, the discretisation and mechanical equations of the structural and transport parts are discussed. At the end of each section, the input parameters for the individual parts are presented. A list of symbols can be found in Table 1.

Discretisation
The dual network approach is based on the Delaunay and Voronoi tessellations of a set of points placed randomly within the domain. The points are placed sequentially while enforcing a minimum distance d min between all points; trial points that fail the minimum distance criterion are rejected. The Delaunay tessellation decomposes the domain into tetrahedra whose vertices coincide with the randomly placed points; the Voronoi tessellation divides the domain into polyhedra associated with the random points [22]. These geometrical arrangements of Delaunay and Voronoi tessellations are used to define the structural and transport elements. Figure 2a shows a Delaunay tetrahedron and the Voronoi facet associated with Delaunay edge i-j. The structural elements are placed on the Delaunay edges with their mid-cross-sections defined by the facets of the Voronoi polyhedra ( Figure 2b). Analogous to the structural network, the transport elements are placed on the edges of the Voronoi polyhedra, with their cross-sections formed by the facets of the Delaunay tetrahedra ( Figure 2c). The discretisation of boundaries of the domain requires special attention. The procedure used in this work is illustrated in Figure 3. Prior to the sequential, random filling of points within the domain, points are placed randomly on the domain surfaces. The minimum distance criterion is enforced during the placement of all of these points. Each interior point is then mirrored with respect to all surfaces of the domain, similar to the procedure of Yip et al. [23]. The tessellations for this set of random points results in Delaunay edges located on the domain surfaces, with their corresponding Voronoi facets traversing the domain boundaries as shown in Figure 3a. Here, Delaunay edge i-j lies on the surface of the domain. Furthermore, Voronoi vertices 1 and 5 are inside, and 2, 3 and 4 are outside the domain. In constructing the transport network, the Voronoi edges within the domain are retained. For edges that cross a surface, only the portion within the domain is kept. For example, edges 1-2 and 4-5 become edges 1-2 and 4 -5, respectively, where nodes 2 and 4 lie on the surface (Figure 3b). These truncated edges define transport elements that are perpendicular to the surface. The modified set of Voronoi edges defines the mid-cross-section of the structural element associated with nodes i and j.
Information exchange between the structural and transport networks is based on the geometrical relationship between neighbouring elements. Herein, a one-way coupling is considered, in which crack openings supplied by the structural network affect the conductivity of the associated transport elements. Details regarding this coupling are provided in Section 2.3. The input parameter for the discretisation is the minimum distance d min which controls the average lengths of structural and transport elements.

Structural Network Model
For the 3D structural analysis, the equilibrium equation for the quasi-static case without body forces [24] is where ∇ is the divergence operator and σ c is the continuum stress. This equilibrium equation is approximated by a network of structural elements.

Structural Element
The discrete version of Equation (1) for the structural element shown in Figure 2b is where K is the stiffness matrix, u e are the vector of degrees of freedom and f s are the acting forces. The formulation of the structural element is presented in the local coordinate system, i.e., the coordinate system (x, y and z) of the nodal degrees of freedom coincides with the coordinate system (n, p and q) of the quantities used for evaluating the constitutive response. Each node has three translational (u x , u y and u z ) and three rotational (φ x , φ y and φ z ) degrees of freedom. The degrees of freedom of a structural element with nodes i and j are grouped in translational and rotational parts as These degrees of freedom u t and u r are used to determine displacement discontinuities u C = u n , u p , u q T at point C by rigid body kinematics [25] as where B 1 and B 2 are two matrices containing the rigid body information for the nodal translations and rotations, respectively, which are and where I is a 3 × 3 unity matrix. In matrix (5), e p and e q are the eccentricities between the midpoint of the network element and the centroid C in the directions p and q of the local coordinate system, respectively (Figure 2b). The local coordinate system is defined by the direction n, which is parallel to the axis of the element, and p and q, which are chosen as the two principal axes of the mid-cross-section.
where h is the length of the structural element. The strains are related to stresses σ = σ n , σ p , σ q T by means of a material stiffness Here, E is the Young's modulus and ω is the damage variable, which is further discussed in Section 2.2.2. Furthermore, γ is an input parameter, which controls Poisson's ratio of the structural network. For γ = 1, Poisson's ratio equal to zero is obtained, which is used in this study. For this case, the structural network is elastically homogeneous under uniform modes of straining. For the case that the global coordinate system coincides with the local one, the element stiffness matrix is Here, K r is a matrix containing the rotational stiffness at point C defined as where I p is the polar moment of area, and I 1 and I 2 are the two principal second moments of area of the cross-section. The factor 1 − ω in matrix (7) ensures that the rotational stiffness reduces to zero for a fully damaged cross-section (ω = 1). For an elastic constitutive model, the present structural element is identical to the one described in Berton and Bolander [26]. All geometrical information of the network element is contained in the element formulation. In this way, the constitutive model relating stresses to strains depends only on properties of the material. This structure is preferred over one that incorporates geometrical information in the constitutive model, since it facilitates the adoption of constitutive modelling frameworks that are commonly used for continuum approaches.

Structural Material
The inelastic structural response of the material during fracture is described by a scalar damage model [27] of the form The damage variable ω is a function of the history variable κ d [28], which is, in turn, determined by the loading function and the loading-unloading conditions The equivalent strain corresponds to an ellipsoidal envelope in the stress space. For pure tensile loading, the stress is limited by the tensile strength f t = Eε 0 . For pure shear and pure compression, the stress is limited by the shear strength f q = q s f t and the compressive strength f c = c s f t , respectively. The damage function is determined by using an exponential stress-crack law in pure tension of the form where w n = ωhε n is the crack opening under monotonic tension and ε n is the tensile strain. This crack opening is the first component of the crack opening vector w = ωhε, which is used for the coupling of the structural and mass transport model. The normal stress in Equation (12) is also expressed in terms of the stress-strain law in Equation (8) as Comparing the right-hand sides of Equations (12) and (13), and replacing ε n by κ d , since a monotonically increasing tensile strain is assumed, the nonlinear equation is obtained from which the damage parameter ω is determined iteratively using a Newton method. In Equation (12), parameter w f determines the initial slope of the softening curve and is related to the fracture energy as G F = f t w f . The input parameters for the structural part of the model are the Young's modulus E, the tensile strength f t , fracture energy G F , shear strength f q and compressive strength f c . These input parameters can be determined from inverse analysis of elementary structural tests of the specific geomeaterial of interest.

Transport Model
For the transport part of the model, a 3D network of 1D transport elements is used to discretise the nonstationary transport equation [29] ∂P c ∂t − div (αgradP c ) = 0, and where P c is the capillary suction, t is the time, α is the conductivity, f is the outward flux normal to the boundary (n-direction) and x is the position in the domain Ω. Furthermore, Γ 1 and Γ 2 are the boundary segments with prescribed suction and flux, respectively. The capillary suction P c in an unsaturated material is defined as P c = P d − P w , where P d is the pressure in the drying fluid and P w is the pressure in the wetting fluid. Here, P d is assumed to be zero, which is a common assumption for modelling the water retention in unsaturated materials subjected to ambient temperatures [29].

Transport Element
The discrete form of Equation (15) for a 1D transport element shown in Figure 2c is where α e and C e are the 1D element conductivity and capacity matrices, respectively, and f are the external fluxes [20,30]. The degrees of freedom of the transport elements are the capillary suction P c = (P c1 , P c2 ) T . Within the context of a one-dimensional finite element formulation [30], Galerkin's method is used to construct the elemental capacity matrix C e as where c is the capacity of the material, A t is the cross-sectional area of the tetrahedron face associated with the transport element (Figure 4), and h t is the length of the transport element. Likewise, based on Galerkin's method [30], the elemental conductivity matrix is defined as where α is the conductivity of the material, which is the sum of two components where α 0 is the initial conductivity of the undamaged material and α c is the change of the conductivity due to fracture.

Transport Material
In the present example, the network approach is applied to mass transport in a general unsaturated geomaterial using techniques introduced originally by van Genuchten for soils [31], but also applied to other geomaterials, such as cementitious materials [32]. According to van Genuchten in [31], the conductivity of the undamaged material α 0 is defined as where ρ is the density of the fluid, µ is the dynamic (absolute) viscosity, κ is the intrinsic permeability and κ r is the relative permeability as a function of the degree of saturation. This degree of saturation is defined as with the moisture content θ, the residual moisture content θ r and the saturated moisture content θ s of the specific geomaterial [31]. Furthermore, the relative permeability κ r is where m is a model parameter [31]. The saturation is related to the capillary suction as where a is another model parameter. Physical justification of parameters m and a in Equations (24) and (25) are given by van Genuchten [31]. The second term in Equation (21) describes the influence of fracture on conductivity. It is defined as wherew i and l ci are the equivalent crack openings and crack lengths (see Figure 4) of neighbouring structural elements, which are located on the edges of the cross-section, and ξ is a tortuosity factor. For mortars, crack tortuosity considered by ξ may reduce flow by a factor of 4 to 6, relative to that between smooth parallel plates [33]. Here,w = |w| is the magnitude of the crack opening w defined in Section 2.2.2. The relation in Equation (26) expresses the well known cubic law, which has shown to produce good results for transport in fractured geomaterials [34]. In Equation (26), w i is assumed to act over l ci (i.e., the equivalent crack opening is uniform over the element crack length). The approach adequately represents variations in opening along the crack trajectory, provided the mesh is sufficiently fine. The way that the crack openings in the structural elements influence the conductivity of a transport element is schematically shown in Figure 4. For instance, for the transport element o-p, three structural elements (i-k, k-j and i-j) bound the cross-section of the transport element. Thus, the conductivity will be influenced by these three elements according to Equation (26) in proportion to their equivalent crack widths and the crack lengths. This crack length (shown by blue double lines in Figure 4) is defined as the length from the midpoint of the structural element to the centroid C t of the transport element cross-section.
The capacity c in Equation (19) is defined as c = −ρ∂θ/∂P c . Using Equation (23), this expression can be written as It is assumed that c is independent of the cracking described by the structural part. The input parameters of the transport part are the density ρ and dynamic viscosity µ of the wetting fluid, the permeability of the saturated uncracked material κ, the saturated and residual wetting fluid content, θ s , and θ r , respectively. Furthermore, parameters m and a of the van Genuchten constitutive model, and the tortuosity parameter ξ are needed.
The structural network is adept at simulating fracture in multi-phase representations of concrete, in which the matrix, aggregates, and matrix-aggregate interfaces are explicitly represented [18,35]. Study of the influence of interface fracture on effective permeability is one potential application of the proposed dual-network approach.

Analyses
In the proposed coupled network approach, the transport elements, which describe both the transport through continuum and fractures, are placed on the edges of the Voronoi polyhedra. This differs from the commonly used approach in which the elements are located at the edges of the Delaunay tetrahedra [20]. The performance of this new approach is investigated by three benchmark tests. The numerical analyses are performed with OOFEM, an open-source object-oriented finite element program [36] extended by the present authors.

Steady-State Potential Flow
For the first benchmark, a homogeneous material is discretised as shown in Figure 5a. The Delaunay/Voronoi discretisation of the domain is based on a set of randomly inserted nodes. Table 2 compares the numbers of nodes/elements forming both network types depicted in Figure 1: the conventional approach (in which transport elements are on the Delaunay edges) and the proposed approach (in which transport elements are on the Voronoi edges). It is clear that the proposed approach is computationally more expensive. The material is subjected to a pressure difference between the x-faces of the domain: P c (x = 0) = 0 and P c (x = L) = 1. For this test, a special case of the constitutive model presented in Section 2.3.2 has been used by assuming the conductivity and capacity to be constant with values of α = c = 1. Both networks accurately represent the steady-state solution, as shown by the nodal potentials plotted in Figure 5. Pressure values are not plotted for the nodes associated with prescribed boundary conditions. The discrete error norms presented in the figures are: where r m = P c (x m ) − P ch (x m ) is the difference between the theoretical and numerical solutions, respectively, at the position of node m; and M is the number of unconstrained nodal points.

Nonstationary Transport Analysis
For the second benchmark, nonstationary mass transport through undamaged material was studied. The geometry and boundary conditions are shown in Figure 6. The two ends of the specimen are subjected to zero pressure whereas all other boundaries are considered to be sealed. For this test, again the special case of α = c = 1 for the constitutive model presented in Section 2.3.2 has been used. The initial condition at all nodes is P c (x, t) = P 0 sin πx L . This assumption allows for a comparison with the analytical solution P c = P 0 sin πx L exp − π 2 L 2 t reported in [20].  For comparing the network results with the analytical results, the vertices were divided into groups with respect to their x-coordinate. For each group of vertices the mean of the x-coordinate and capillary suction are presented. The Voronoi-edge based network agrees well with the analytical solution of the capillary suction distribution without exhibiting any dependence on the element size. Any differences between the numerical and analytical solution originate from the time discretisation, rather than the new spatial discretisation.

Coupled Structural-Transport Benchmark
In the third benchmark, the structural and transport models are coupled. Firstly, a double cantilever beam is used to assess the capability of the structural model to describe fracture without any pathological network dependence. Then, fluid transport through the fractured specimen at an intermediate loading stage of the structural analysis is modelled for different networks with different element sizes. The geometry and loading setup for the structural and transport tests are shown in Figure 9a,b, respectively. For the structural analysis, the load is applied at x = 0.25L. For the transport component of the analyses, three networks with minimum distances between Delaunay vertices of d min /L = 0.06, 0.045 and 0.03 were used. The structural and transport networks with d min /L = 0.06 are shown in Figure 10a,b, respectively. As noted earlier, transport elements local to the boundaries are perpendicular to the specimen surfaces. The input parameters for the structural constitutive model are E = 30 GPa, f t = 3 MPa and G F = 120 N/m, which are representative of concrete materials. A notch of length 0.25L is introduced by reducing the tensile strength f t of elements crossing the notch to 1% of the original value. The load-displacement curves from the structural benchmark for the three networks are shown in Figure 11. There is little difference between the responses obtained with the three networks. Fracture is indicated by shading the mid-cross-sections of elements in which the equivalent crack opening has reached a threshold value. The mid-cross-sections of elements with damage corresponding to an equivalent crack openingw > 10 µm are shown in Figure 12 for the three different networks at a load-point-displacement of δ = 0.15 mm in Figure 11.   Figure 11. The shaded polygons represent the mid-cross-sections of elements with w > 10 µm.
The transport network uses the same geometry as in the nonstationary transport test in Section 3.2. However, the boundary and initial conditions, and the material input parameters, are changed so that the influence of fracture could be studied more effectively. On the left-hand side of the model, the boundary is subjected to P c = 0. Furthermore, the initial capillary suction of all other nodes is set to P c = 1.736 MPa, which for the chosen material parameters corresponds to an initial saturation of S init = 0.5. Other input parameters for the transport problem are: α 0 = 1 × 10 −17 m 2 , θ s = 0.1, θ r = 0, a = 1 MPa, m = 0.5 and ξ = 0.001. The transport analysis is performed for crack patterns obtained at a displacement of δ = 0.15 mm in Figure 11.
Results for the cumulative volume of inflow at the left side of the specimen normalised by the available volume to be filled, from the time of initial wetting, are presented in Figure 13. The available volume to be filled is V avail = (1 − S init ) θ s V tot , where V tot = L × 0.25L × 0.25L is the total volume of the specimen. The inflow is practically independent of the element size.
Furthermore, contour plots of the capillary suction P c are shown for the three networks for the x-z plane (at y = 0.125 m) and for the y-z plane (at x = 0.3 m) in Figure 14. Darker regions correspond to lower values of capillary suction, which indicate higher amounts of intruded water. Slight broadening of the intrusion zone, lateral to the crack direction, is expected for the coarser network design. Otherwise, the network model simulates the transport field element size independently. Whereas this example involves mode I fracture, the scalar damage model presented in Section 2.2.2 allows for damage development under more general loading patterns. Modification of conductivity to account for fracture, according to Equations (21) and (26), is appropriate when the crack is open (i.e., when w n > 0). In this sense, the proposed model should be applicable to cases of mixed-mode loading within the tension-shear regime. Residual influences after crack closure, or possible modification of conductivity due to damage in the compression-shear regime, require additional study.

Conclusions
A new three-dimensional network approach for modelling the effect of fracture on mass transport has been proposed. The Delaunay tessellation of an unstructured set of points defines the structural network, which represents material elasticity and fracture. The edges of the corresponding Voronoi diagram define the network of transport elements, which simulate mass transport. A distinctive feature of the dual network approach is the alignment of transport elements with potential pathways for crack propagation. Several benchmark comparisons have been presented involving non-stationary transport, fracture, and their coupling. The following conclusions and remarks can be made.

•
The network of structural elements, defined by the Delaunay edges, provides element geometry and size independent load-displacement curves, as demonstrated through cohesive fracture simulations of double cantilever beams. The traction free condition is approached without stress locking. Local deviations of the fracture path due to random network generation has very little influence on the load-displacement curves. • The network of transport elements, defined by the Voronoi edges, provides results for non-stationary transport which are in very good agreement with analytical solutions, and are independent of element geometry and size. The proposed discretisation scheme for the transport network facilitates the enforcement of boundary conditions. Local to a domain boundary, transport elements have one node on the boundary and are directed perpendicular to the boundary. • The proposed method for coupling the effect of crack opening, determined by the structural network, with transport properties of the transport network yields objective results with respect to element geometry and size. This dual network approach facilitates the simulation of transport along crack paths and from crack faces into the bulk material.
The proposed coupling is limited to the effect of fracture on transport. A two-way coupling of field quantities (i.e., including the dependence of structural behaviour on the transport field [13]) is a natural extension of this work.