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.

1. 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.

Figure 1

Network models for coupled problems: (a) common approach in which the structural and transport network nodes are coincident. Both structural and transport elements are on the Delaunay edges; (b) simulated crack in structural network; and (c) improved approach in which transport elements are on the Voronoi edges and therefore aligned with potential cracks.

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.

2. 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.

Table 1

Nomenclature.

Symbol	(Units)	Definition
A_t	(m^2)	cross-sectional area of the tetrahedron face
a	(Pa)	parameter in van Genuchten model
B_1, B_2		matrices expressing rigid body kinematics
c	(s^2/m^2)	capacity of the material
C_e	(m s^2)	element capacity matrix
C	(m)	centroid of mid-cross-section
c_s		ratio of compressive and tensile strength
D_e	(Pa)	material stiffness matrix
d_min	(m)	minimum distance between nodes
e_p, e_q	(m)	eccentricities between the midpoint of the network element and the centroid C
E	(Pa)	Young’s modulus
f_s	(N)	acting structural forces
f_d		loading function
f_t	(Pa)	tensile strength
f_q	(Pa)	shear strength
f_c	(Pa)	compressive strength
f	(kg/m^2)	outward flux normal to the boundary
G_F	(J/m^2)	fracture energy
h	(m)	length of structural element
h_t	(m)	length of transport element
I		unity matrix
I_p	(m^4)	polar moment of area
I_1 and I_2	(m^4)	two principal second moments of area of the cross-section
K		element stiffness matrix
K_r		rotational stiffness at point C
l_c	(m)	crack length
L	(m)	length of specimen
m		parameter in van Genuchten model
n, p, q	(m)	coordinates of mid-cross-sections
P_c	(Pa)	capillary suction (tension positive)
P_w	(Pa)	pressure in the wetting fluid
P_d	(Pa)	pressure in the non-wetting fluid
q_s		ratio of shear and tensile strength
S		degree of saturation
t	(s)	time
u_n, u_p, u_q	(m)	displacement discontinuities
u_x, u_y, u_z	(m)	translational degrees of freedom
u_e		vector of degrees of freedom of structural element
u_t	(m)	vector of translational part of degrees of freedom
u_r		vector of rotational part of degrees of freedom
u_C	(m)	vector of displacement discontinuities
V	(m^3)	volume
V_avail	(m^3)	available volume to be filled
V_tot	(m^3)	total volume of the specimen
w_n, w_p and w_q	(m)	crack opening components
w_f	(m)	displacement threshold which determines the initial slope of the softening curve
w˜	(m)	equivalent crack opening
x, y, z	(m)	Cartesian coordinates
α_0	(s)	initial conductivity of the undamaged material
α_c	(s)	change of the conductivity due to fracture
α	(s)	conductivity
α_e		conductivity matrix
γ		input parameter, which controls Poisson’s ratio of the structural network
Γ_1, Γ_2		boundary segments
δ	(m)	load-point-displacement
ε		strain vector
ε_n, ε_p, ε_q		strain components
ε_0		strain threshold
θ	(kg/m^3)	moisture content
θ_r	(kg/m^3)	residual moisture content
θ_s	(kg/m^3)	saturated moisture content
κ	(m^2)	intrinsic permeability
κ_d		history variable in damage model
κ_r		relative permeability
μ	(Pa s)	dynamic (absolute) viscosity
ξ		tortuosity factor
ρ	(kg/m^3)	density of the fluid
σ^c	(Pa)	continuum stress
σ	(Pa)	stress vector
σ_n, σ_p, σ_q	(Pa)	stress components
ϕ_x, ϕ_y, ϕ_z		rotational degrees of freedom
ω		damage variable
∇		divergence operator

2.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 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).

Figure 2

Spatial arrangement of structural and transport elements of the 3D transport-structural network approach: (a) geometrical relationship between Delaunay and Voronoi tessellations; (b) structural element with cross-section defined by the associated Voronoi facet; and (c) transport element with cross-section defined by the associated Delaunay facet.

The discretisation of boundaries of the domain requires special attention. The procedure used in this work is illustrated in Figure 3.

Figure 3

Discretisation of domain boundaries: (a) Voronoi facet of Delaunay edge i–j located on the surface of the domain after initial tessellation; and (b) modified arrangement used for definition of transport nodes and elements.

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 which controls the average lengths of structural and transport elements.

2.2. 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 is the continuum stress. This equilibrium equation is approximated by a network of structural elements.

2.2.1. Structural Element

The discrete version of Equation (1) for the structural element shown in Figure 2b is where is the stiffness matrix, are the vector of degrees of freedom and 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 (, and ) and three rotational (, and ) degrees of freedom. The degrees of freedom of a structural element with nodes i and j are grouped in translational and rotational parts as , where and . These degrees of freedom and are used to determine displacement discontinuities at point C by rigid body kinematics [25] as where and are two matrices containing the rigid body information for the nodal translations and rotations, respectively, which are and where is a unity matrix. In matrix (5), and 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.

The displacement jump in Equation (3) is transformed into strains , where h is the length of the structural element. The strains are related to stresses by means of a material stiffness , where . 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 , 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, is a matrix containing the rotational stiffness at point C defined as where is the polar moment of area, and and are the two principal second moments of area of the cross-section. The factor in matrix (7) ensures that the rotational stiffness reduces to zero for a fully damaged cross-section (). 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.

2.2.2. 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 [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 . For pure shear and pure compression, the stress is limited by the shear strength and the compressive strength , respectively.

The damage function is determined by using an exponential stress-crack law in pure tension of the form where is the crack opening under monotonic tension and is the tensile strain. This crack opening is the first component of the crack opening vector , 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 by , 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 determines the initial slope of the softening curve and is related to the fracture energy as . The input parameters for the structural part of the model are the Young’s modulus E, the tensile strength , fracture energy , shear strength and compressive strength . These input parameters can be determined from inverse analysis of elementary structural tests of the specific geomeaterial of interest.

2.3. 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] subject to and where is the capillary suction, t is the time, α is the conductivity, f is the outward flux normal to the boundary (-direction) and is the position in the domain Ω. Furthermore, and are the boundary segments with prescribed suction and flux, respectively. The capillary suction in an unsaturated material is defined as , where is the pressure in the drying fluid and is the pressure in the wetting fluid. Here, is assumed to be zero, which is a common assumption for modelling the water retention in unsaturated materials subjected to ambient temperatures [29].

2.3.1. Transport Element

The discrete form of Equation (15) for a 1D transport element shown in Figure 2c is where and are the 1D element conductivity and capacity matrices, respectively, and are the external fluxes [20,30]. The degrees of freedom of the transport elements are the capillary suction . Within the context of a one-dimensional finite element formulation [30], Galerkin’s method is used to construct the elemental capacity matrix as where c is the capacity of the material, is the cross-sectional area of the tetrahedron face associated with the transport element (Figure 4), and is the length of the transport element.

Figure 4

Influence of cracking on transport.

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 is the initial conductivity of the undamaged material and is the change of the conductivity due to fracture.

2.3.2. 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 is defined as where ρ is the density of the fluid, μ is the dynamic (absolute) viscosity, κ is the intrinsic permeability and 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 and the saturated moisture content of the specific geomaterial [31]. Furthermore, the relative permeability 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 where and 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, 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), is assumed to act over (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 of the transport element cross-section.

The capacity c in Equation (19) is defined as . 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, and , 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.

3. 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.

3.1. Steady-State Potential Flow

For the first benchmark, a homogeneous material is discretised as shown in Figure 5a.

Figure 5

Steady-state simulation of potential flow: (a) Voronoi tessellation of domain; (b) conventional network solution; and (c) proposed network solution.

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.

Table 2

Network feature counts.

Network Type	Node Definition	Element Definition	Nodal Count *	Element Count *
Conventional	Delaunay vertex	Delaunay edge	330	1800
Proposed	Voronoi vertex	Voronoi edge	2880	5440

* rounded to nearest ten.

The material is subjected to a pressure difference between the x-faces of the domain: and . 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 . 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 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.

3.2. 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 for the constitutive model presented in Section 2.3.2 has been used. The initial condition at all nodes is . This assumption allows for a comparison with the analytical solution reported in [20].

Figure 6

Geometry and boundary conditions for the nonstationary transport benchmark.

Three transport networks with minimum distances between Delaunay vertices of , and are used. The coarse network with is shown in Figure 7.

Figure 7

Coarse network for the transport analysis.

The capillary suction distributions for four time steps and the three different network sizes are shown in Figure 8 for a symmetric portion of the model.

Figure 8

Capillary suction distributions due to non-stationary transport.

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.

3.3. 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 .

Figure 9

Geometry and loading setup of the structural (a) and transport (b) benchmark.

For the transport component of the analyses, three networks with minimum distances between Delaunay vertices of , and were used. The structural and transport networks with 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 GPa, MPa and N/m, which are representative of concrete materials. A notch of length is introduced by reducing the tensile strength of elements crossing the notch to 1% of the original value.

Figure 10

The coarse dual networks () for (a) structural and (b) transport analysis.

The load-displacement curves from the structural benchmark for the three networks are shown in Figure 11.

Figure 11

Load versus load-point-displacement results for three networks.

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 opening m are shown in Figure 12 for the three different networks at a load-point-displacement of mm in Figure 11.

Figure 12

Crack patterns for (a) coarse; (b) medium; and (c) fine network for a load-point-displacement of mm in Figure 11. The shaded polygons represent the mid-cross-sections of elements with 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 . Furthermore, the initial capillary suction of all other nodes is set to MPa, which for the chosen material parameters corresponds to an initial saturation of . Other input parameters for the transport problem are: m, , , MPa, and . The transport analysis is performed for crack patterns obtained at a displacement of 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 , where is the total volume of the specimen. The inflow is practically independent of the element size.

Figure 13

Influence of element size on the cumulative volume of inflow normalised by the domain volume.

Furthermore, contour plots of the capillary suction are shown for the three networks for the x–z plane (at m) and for the y–z plane (at 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.

Figure 14

Contour plots of capillary suction at 3.33 h for the (a) x–z plane at m and (b) y–z plane at m.

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 ). 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.

4. 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.