Multi-Field Models of Fiber Reinforced Concrete for Structural Applications

Short fiber reinforcement is added to concrete materials to improve a variety of performance measures related to structural safety, serviceability, and health diagnosis. Mesoscale models have been developed to understand the individual and collective actions of these fibers on various material properties. Those modeling efforts have predominately focussed on the mechanical (i.e., stiffness and strength) contributions of fibers. This paper introduces computationally efficient, semi-discrete representations of fibers within coupled mechanical and transport processes in cement-based matrices. Basic simulations are done to study: the use of conductive fibers for self-sensing; and the influences of fibers on early-age plastic settlement. It is found that the models can account for directional bias on fiber orientation, as may occur during material casting. With respect to plastic settlement, fibers may play competing roles: mechanical restraint offered by the fibers reduces settlement, whereas enhanced hydraulic conductivity along the fiber–matrix interface may increase settlement by facilitating the bleeding process.


Introduction
Short-fiber reinforcement is often added to concrete to improve its mechanical properties [1], particularly fracture toughness. Short fibers are also used to modify other aspects of concrete, including volumetric stability at early ages [2][3][4], resistance to spalling at high temperatures [5,6], bond properties of repair materials [7], structural composite behavior with steel reinforcement [8][9][10], and various electrical properties [11][12][13]. Along with the amount of fibers, the efficacy of fiber additions often depends on the distribution of fibers and characteristics of the fiber-matrix interface.
Effective medium theories have been developed to understand the global significance of fiber additions on the properties of interest (e.g., stiffness, strength, conductivity) [14]. The usefulness of such theories is limited, however, due to assumptions made regarding the distribution of fibers within the material volume. Numerical models have been developed to understand the individual and/or collective actions of these fibers on various performance measures. These modeling efforts have predominately focussed on the mechanical (i.e., stiffness and strength) contributions of fibers in both pre-and post-cracking conditions [15][16][17][18][19].
This paper introduces computationally efficient, semi-discrete representations of individual fibers within cement-based matrices that participate in electro-mechanical or hygro-mechanical phenomena. Voronoi Cell Lattice Models are used to represent the composite materials [20]. The fibers are discretized within the computational domain, such that the effects of fiber spatial distribution and orientation can be explicitly represented [21]. Basic simulations are conducted to study: the use of electrically conductive fibers for self-sensing; and the influences of fibers on early-age plastic settlement. The methods are verified through benchmark comparisons with numerical and theoretical models. It is found that the models can account for unintended directional bias on fiber orientation, as may occur during material casting. With respect to plastic settlement, fibers may play competing roles: mechanical restraint offered by the fibers reduces settlement, whereas enhanced hydraulic conductivity along the fiber-matrix interface may increase settlement by facilitating the bleeding process. Although the simulation efforts are limited to small component or unit-cell models, the low computational cost of the proposed approaches suggests future applications to structural concrete components.

Voronoi-Cell Lattice Model
Two related forms of lattice models are used to represent the variables associated with the multi-field phenomena considered herein. The mechanical behavior of the fiberreinforced concrete is modeled using Rigid-Body Spring Networks, whose geometry is based on Voronoi/Delaunay tessellations of the nodal points. The potential field analyses utilize lattice models composed of conduit elements, based on the same connectivity defined by the Delaunay tessellation. In the present study, the collective approach is called the Voronoi-cell lattice model (VCLM) to emphasize the central role of the Voronoi diagram in establishing its desirable mechanical and potential flow qualities.
Homogeneous regions of the material body are typically represented by randomly placed nodes, subject only to a minimum distance constraint. To construct material features within the domain, subsets of nodes are sequentially placed to define the vertices, edges, and faces of the desired features, prior to random filling of the feature volumes [33]. The RBSN is an assembly of rigid-body-spring elements. Each element is composed of a zero-size spring set, located at the centroid C of the corresponding Voronoi facet, connected to the nodes via rigid-arm constraints ( Figure 1b). As mentioned above, element connectivity is determined by the Delaunay tessellation of the nodes; the dual Voronoi tessellation provides geometric measures that are used to determine element properties. With respect to the n − s − t coordinate axes, defined local to Voronoi facet ij (Figure 1b), the stiffnesses of the axial and rotational springs composing the spring set are where E is the elastic modulus of the material; A ij is the area of Voronoi facet; h ij is the distance between i and j; I s and I t are the second moments of the facet area about axes s and t, respectively; J p is the polar second moment of area of the facet; and β is a factor related to macroscopic Poisson ratio ν. Herein, β = 1 renders the lattice elastically uniform, albeit with ν = 0 [34]. Complete representation of the Poisson effect, both macroscopically and in an element-local sense, has recently been accomplished using an iterative procedure [35]. The spring coefficients of Equations (1) and (2) form the diagonal of the material matrix where 0 ≤ ω ≤ 1 is a scalar damage parameter used to model material fracture. Prior to fracture initiation, ω = 0. The element stiffness matrix (with respect to element local coordinates) is where B relates the spring generalized displacements to the nodal degrees of freedom of the element [22] B = −I B 12 I B 14 0 −I 0 −I (5) in which submatrices B 12 and B 14 are defined by Recently, this elemental formation has been extended to accommodate large displacement analyses [36,37]. After transforming K e to global coordinates, the direct stiffness approach is used to assemble element stiffness matrices and internal force contributions into the structural equation set in which ∆u is the increment in generalized nodal displacements; K is the system stiffness matrix assembled from the elemental components; R and F are the external and internal nodal force vectors, respectively. One advantage of this form of discrete model is in its ability to simulate tensile fracture with reduced mesh bias. In particular, implementations of the crack band model of Bažant et al. [38,39] within the rigid-body-spring model are objective with respect to mesh size and geometry [34]. Volumetric instabilities due to creep and thermal/hygral processes can be simulated by introducing rheological constructs and hygral/thermal units in the place of the axial springs of the rigid-body-spring elements [20]. Although volumetric instabilities and fracture are not considered in the simulations conducted herein, such modeling capabilities have been validated [18,34,40] and could be incorporated within extensions of this research.

Transport Behavior
Flow behavior is simulated using the same lattice topology, based on the Delaunay tessellation of the nodal points. The lattice elements serve as conduits of one-dimensional flow between the element nodes [20]. The cross-section area of a given conduit element is equal to the area of the Voronoi facet associated with the element nodes (Figure 1c). From a nodal point perspective, the lattice accurately simulates uniform flow through homogeneous materials subjected to a potential difference (i.e., the random lattice design does not affect the condition of uniform flow). The lattice model accommodates both steady-state and transient solutions. In cases where crack-assisted flow is of interest, a dual lattice approach is effective [41][42][43][44]. The field of displacements (for mechanical behavior) is represented at the Voronoi generator points, whereas the potential field is represented at the Voronoi vertices. Cracking is not considered in this study, so a conventional approach (where the mechanical and flow network nodes are coincident with the Voronoi generator points [20]) is employed.
The case of moisture diffusion, which is used to describe the modeling approach, is analogous to the cases of electrical and hydraulic conductivity considered later. Established approaches are used for modeling moisture diffusion in cement composites, the convective boundary conditions at exposed surfaces, and the coupling of the diffusion and strain analyses [45]. For isothermal conditions, with no sinks or sources present, the governing equation for transient diffusion problems is where h is the pore relative humidity and D is a humidity dependent diffusion coefficient, which is a scalar quantity for isotropic conditions. The problems considered hereafter do not involve the dependence of D on the field quantity. The semi-discrete form of Equation (8) is where H is the vector of nodal relative humidities and the dot over H indicates time derivative. The system matrices, M and K, are assembled from the elemental capacity and diffusivity matrices, respectively: as would be done for the finite element analysis of potential flow problems using twonode elements [46]. In the elemental capacity matrix, d = 1.0, 2.0, and 3.0 for 1D, 2D, and 3D networks, respectively. The ensuing discussions pertain to electrical or hydraulic conductivity, rather than diffusivity, but the equations representing potential flow take the same general form.

Models of Fiber Reinforcement
Much effort has been devoted to modeling the strength and stiffness contributions of fibers to the behavior of hardened concrete. There has been interest in explicitly representing individual fibers, such that their local and collective effects can be simulated [21,[47][48][49][50][51][52]. Such models permit the study of fiber distribution non-uniformity and its influence on the desired properties of the composite material, including cracking behavior during strain hardening [18,53].
A prime consideration is the balance between computational efficiency and accuracy in representing each fiber. On one end of the spectrum, fibers can be modeled using volume elements. This allows for accurate representation of the elasticity of the fiber and the possibility for directly modeling the fiber-matrix interface. Interface features (e.g., the lugs on ordinary steel reinforcing bars) can be explicitly represented [54,55]. This approach is computationally expensive, however, which limits the number of fibers that can be considered.
The number of computational degrees of freedom can be drastically reduced by modeling fibers as series of simple structural elements (e.g., line elements such as truss or frame elements). These line elements interact with the matrix via bond-link elements that account for relevant behaviors of the fiber-matrix interface [23,56]. Going a step further, by constraining the actions of the fibers to the kinematics of a particle network representing the cement-based matrix, fiber additions do not inflate the number of computational degrees of freedom [21]. The research presented herein employs this semi-discrete approach, which enables the inclusion of larger numbers of fibers. To have benchmark results for comparison, volumetric representations of individual fibers are also considered. For both approaches, emphasis is placed on modeling the individual and collective effects of fibers on transport processes, which has received scarce attention within the cement-based composites community.
The positioning of the fibers needs to be specified. Fiber locations can be based on sectional analyses [57][58][59] or 3D imaging [17,[60][61][62], which have the advantage of capturing directional and concentration biases caused by the casting process. Alternatively, as done herein, a procedure is used to sequentially place fibers in the computational domain without positional or directional bias. The spatial coordinates of the centroid of a trial fiber are assigned using a uniform random number generator. The orientation angles of the fiber are then determined using the algorithm of Knop [63]. To avoid fiber intersections, a minimum allowable distance is maintained between fibers [64]. Trial fibers that meet the minimum distance constraint are retained within the material volume. For the unit cell simulations of transport phenomena presented hereafter, periodic boundary conditions are implemented for fibers intersecting the domain boundaries.

Volumetric Representation of Fibers
The material is modeled as a two-phase composite consisting of short-fiber inclusions embedded within a homogeneous matrix. The cross-section and length of fibers form part of the continuum, as shown in Figure 2. Individual fibers are discretized by placing nodal points within the domain using regular patterns. Herein, cylindrical fibers are approximately modeled (Figure 2b), but other cross-section shapes are also possible. The matrix phase, which forms the remainder of the material volume, is then populated with randomly placed nodes. The lattice model formed from this nodal point set accommodates analyses of both the scalar field(s) associated with potential flow and the vector field of displacements associated with mechanical behavior. For the potential flow analyses, periodic boundary conditions are realized by: (1) positioning fibers irrespective of the domain boundaries; (2) placing any external portion of the fiber on the opposite side of the domain; and (3) introducing a special conductivity element that connects the nodes on the opposing faces, such that continuity of flow is maintained over the fiber length. This is shown in Figure 2c for the two-dimensional case.
In Figure 2, the lighter colored fibers intersect one of the faces of the unit cell, whereas the darker colored fibers reside fully within the unit cell. The same discretization scheme is used for simulating mechanical behavior, but without the mechanical connection between corresponding nodes on opposing faces.

Semi-Discrete Representation of Fibers
In contrast to the volumetric representation of fibers described above, fibers can be freely positioned within the computational domain, independent of the background mesh representing the matrix material. A fiber element is produced wherever the fiber intersects a Voronoi facet, as shown in Figure 3. The fiber element introduces stiffness and/or conductivity, as the case may be, between nodes i and j associated with the Voronoi facet. In this way, the addition of fibers does not increase the number of computational degrees of freedom. For the illustrative examples presented in Sections 4 and 5, the fibers are placed uniformly at random within the composite material domain. Other arrangements of fibers are possible, such as those resulting from casting [61] or extrusion processes. The semidiscrete fiber models are generally applicable to different fiber types, in the sense that the program routines utilize fundamental properties of the fibers (e.g., length, diameter, stiffness and/or conductivity of the fibers.) For cases where slippage occurs along the fiber-matrix interface, bond properties must also be specified.

Mechanical Behavior
An axial spring is inserted at each point of intersection of a fiber with a Voronoi facet [21]. The spring is aligned with the fiber direction. The force-displacement actions of the spring are related to the local (i − j) nodal degrees of freedom through rigid-arm constraints, similar to those of the rigid-body-spring network. Prior to matrix fracture, the fibers are assumed perfectly bonded with the matrix. The axial stiffness of the spring associated with each fiber element is determined using shear-lag theory [65]. In this way, stress transfer between the matrix and fibers appropriately depends on fiber aspect ratio and modular ratio.
This semi-discrete approach can simulate the debonding and pullout of fibers that bridge matrix cracks [21], based on the micro-mechanical model of Naaman et al. [66]. During the debonding/pullout processes, load is transferred between the material phases along the embedded lengths of the fibers, independently of mesh size [18,67]. For the illustrative examples considered herein, however, slippage along the matrix-fiber interface is not considered. The fibers are assumed to be fully bonded to the surrounding matrix material.

Transport Behavior
Akin to the approach for modeling the stiffness contributions of fibers, fibers contribute to the diffusivity (or conductivity) of the composite material wherever they intersect a Voronoi facet ( Figure 3). The elemental form of the diffusivity contribution between nodes i and j is similar to that of Equation (11) where A f is the cross-sectional area of the fiber; D f is the diffusivity of the fiber material in its axial direction; and h f is the projected length of h ij in the direction of the fiber. For the planar illustration of the lattice network shown in Figure 3, h f = h ij cos θ. In this way, the total length of the fiber is preserved when crossing multiple Voronoi facets where N is the number of facets crossed by the fiber. In Equation (13), strict equality is achieved if the fiber end conditions are considered.

Transport Behavior: Model Verification
The discrete (volumetric) approach and proposed semi-discrete approach are compared for the case of a fiber(s) embedded in a plate of homogeneous material. Domain discretization for the discrete approach is shown by the Voronoi diagram in Figure 4a. A similar diagram is used for the semi-discrete approach, but without the structured arrangement of cells representing the fiber. For both cases, the fiber is centrally positioned and aligned with the direction of flow, which is induced by a potential difference ∆φ = 1 acting across the model domain. The top and bottom surfaces act as perfect insulators. For the purposes of illustration, the fiber has aspect ratio l f /d f = 50 and conductivity ratio σ f /σ m = 1000. By comparing potential values along the plane of symmetry (Figure 4b,c), close agreement between the model results is evident. The semi-discrete approach is versatile in that fibers of different length, orientation, etc., can be placed in the domain without modification of the background mesh representing the matrix phase, as shown by the results in Figure 4d. Notably, the random lattice discretization of the matrix phase of the material does not produce artifacts in the contour plots. This lack of artificial heterogeneity in the background lattice allows for more precise descriptions of the fiber contributions to composite behavior.

Background
With respect to the electrical properties of fiber concrete, conductivity (or resistivity) of the composite material is commonly of interest. Changes in composite resistance due to applied loading can be measured and used for self-sensing of structural concrete [12]. Variations in voltage along the periphery of fibrous sensing skins have been used to produce conductivity maps and characterize damage within the enclosed region of interest [68]. Such self-sensing capabilities are potentially useful for the condition assessment of structural concrete.
Three behavioral regimes have been used to characterize composite conductivity. At low fiber contents, the fibers are electrically isolated from one another within the lowly conductive cement-based matrix. In this condition, conductivity of the composite increases with increasing fiber content, but at a relatively low rate. Effective medium theories have been developed to describe the dependence of composite conductivity on the amount, aspect ratio, conductivity, and orientation distribution of the fibers [14]. By increasing the amount of fibers, the number of direct contacts (or near contacts) between fibers also increases. Within the second regime, dramatic improvements in composite conductivity result from the establishment of those longer conductive pathways and, ultimately, percolation clusters across the region of interest [11,12]. In the third regime, composite conductivity continues to increase with increasing fiber content, but at a lower rate after achieving a sufficient degree of percolation. The simulations presented herein mainly involve the first behavioral regime, due to the relatively low fiber dosages and low aspect ratios considered.

Solution Scheme
The electro-mechanical coupling is established using a staggered solution scheme. The base resistance of the numerical specimen is calculated for the case of no applied loading. Loading is then applied incrementally. For each load increment, the nodal positions are determined. The associated element strains are then used to calculate updated values of composite resistance. For the examples considered herein, planar arrays of nodes are placed along the negative and positive x-faces of the 3D domain ( Figure 5). These nodes serve as electrode surfaces or grids for imposing a potential difference in the x-direction. The planar positioning of the nodes also facilitates the calculation of electrical current between the opposing faces.
Fibers are placed uniformly at random within the cell. At this stage of model development, lower fiber contents are analyzed and therefore fiber-to-fiber contact is not a consideration. The aspect ratio of the fibers is l f /d f = 50, where l f is 0.5L. The conductivity ratio is set at σ f /σ m = 1000. Whereas this ratio is much smaller than that between conductive fibers and concrete, which can be on the order of 10 8 [13], it is sufficient for examining basic operations of the model. For the discrete (volumetric) modeling of fibers, the fibers are discretized as indicated in Figure 2b. For both the discrete and semi-discrete models, random filling of the matrix volume with nodal points employed a minimum distance constraint of L/33.
The influence of conductive fibers on potential flow is illustrated in Figure 6, which is based on earlier simulations with lesser amounts of fibers. Flow initiating on the left side of the figure moves within and between the conductive fibers. The color coding of the streamlines is not a quantitative indicator but was done to distinguish the different flow paths.

Verification of Conductivity Simulations
As primary conditions for simulating electrical behavior of fiber-reinforced composites, the model should account for the amount, aspect ratio, and orientation of the fibers, as well as the differing conductivities of the fiber and matrix materials. Effective medium theories have been useful in estimating fiber composite properties (including conductivity) when basic requirements are met regarding the orientation distribution of the fibers. In this study, the equivalent inclusion method (EIM) of Hatta and Taya [14] provides benchmark values for verifying the lattice models. Whereas the EIM was developed for heat conduc-tion problems, it is analogous to Eshelby's equivalent inclusion method in elasticity [69]. The EIM applies to potential flow problems, in general, including the electric conductivity simulations presented herein. The expression for determining composite conductivity is given in Appendix A. For the assumed model parameters, Figure 7 compares estimates of composite conductivity according to the EIM of Hatta and Taya with bounds of the Hashin-Shtrikman solution [70]. To verify the lattice models, a potential difference ∆φ = 1 is applied across the unit cell shown in Figure 5. As previously mentioned, the potential is assigned to nodal arrays along each surface normal to the x-direction, separated by distance L e , according to φ(x) = ∆φx/L. Fibers were not placed close to the boundaries, which facilitated domain discretization and the assignment of nodal potentials. The influence of this condition on the conductivity simulations, which appears to be secondary, was modeled using parallel and series constructions of the regions with and without fibers [71]. Conductivity of the composite material, σ, was determined using where J is the average current density and φ(L) − φ(0) is the potential difference between the faces of the unit cell. Average current density was found by summing nodal flux values over either of the x-faces, and dividing by the face area. The procedure for calculating nodal flux values is outlined in [20]. Results of the discrete (volume element) approach to modeling fibers are compared with the EIM values in Figure 8, for small volume fractions of fibers up to V f = 0.01. The limitation on fiber dosage is due, in part, to the computational expense of models with larger volume fractions. For higher fiber contents, contact between fibers is difficult to avoid, along with the eventual formation of percolation clusters. Seven numerical specimens, differing only in their random distributions of fibers, were produced for each fiber volume fraction considered.
Several observations can be made regarding the conductivity values plotted in Figure 8. The lattice model results exhibit the expected trend: conductivity of the composite material increases as the dosage of conductive fibers increases. The mean values given by the lattice model agree well with those of the EIM of Hatta and Taya. In a separate study, experimental measurements of AC composite conductivity of steel-fiber reinforced cementitious materials agreed well with the EIM, when some degree of preferential alignment of the fibers was considered [13]. The simulated conductivity values exhibit scatter about the respective mean values. The error bars in Figure 8a indicate one standard deviation from the mean value. In Figure 8b, the simulated conductivity values have been normalized by the corresponding EIM values, σ EI M , such that the entire set of results can be plotted together. These results indicate that there is correlation between the degree of fiber alignment in the x-direction and the conductivity values. The degree of alignment is expressed by the average of cos α, where α is the angle between each fiber axis and the x-direction. As expected, preferential alignment produces higher conductivity values. The differing random realizations of the fiber-matrix pathways for conductive flow may also contribute to the scatter in the conductivity values. The model capabilities are potentially useful in understanding the effects of alignment bias due to the production/casting processes and container boundaries. Non-uniform or poorly aligned fiber distributions may strongly influence the fracture properties of fiber-reinforced concrete [57,60,61,72]. The severity of alignment bias on composite conductivity has not been examined herein, since all fibers distributions are, in a nominal sense, uniformly at random.
The same set of parameters is used for simulating composite conductivity using the proposed semi-discrete modeling approach. Regardless of the fiber dosage, the number of computational nodes representing the composite material, n c , does not change and is equal to n b (≈16,000), which is the number of nodes defining the background lattice representing the matrix material. This is in contrast to the discrete approach for which computational cost grows with fiber dosage, as shown in Figure 9. For perspective, nodal counts for a case of a higher fiber aspect ratio are included in the figure. It is apparent that the full volumetric discretization of the fibers becomes cost prohibitive as the volume fraction or aspect ratio of the fibers increases. In other words, the proposed semi-discrete approach has a clear advantage with respect to its lower computational complexity. For the parameter set and maximum fiber dosages considered herein, n c /n b = 7.8, which is nearly an order of magnitude difference in the number of nodes.
As seen in Figure 10, similarly good agreement with the EIM curve is observed for lower V f values, but the model results tend to drift upward as V f grows. At lower V f values, the background mesh is sufficiently fine such that the fibers are electrically isolated from one another. In other words, fiber elemental contributions to diffusivity (as per Equation (12)) are introduced to different matrix nodes, such that any two fibers are effectively separated by matrix material. The discrete approach satisfies this condition by design of the fiber cross-section, as shown in Figure 2b. The semi-discrete fibers could remain electrically isolated for much larger V f provided the density of nodes representing the matrix was increased accordingly. However, the computational efficiency advantages of the semi-discrete approach would then be lessened. Nonetheless, these simulations verify accuracy of the semi-discrete approach for modeling the transport behavior of fibers. Notably, the semi-discrete approach also captures the influence of directional bias of the fibers on composite conductivity, as shown in Figure 10b. In a separate study, the mechanical analogue of this fiber model was used to study the potential influences of fiber distribution non-uniformity within steel-fiber reinforced concrete slabs [61]. The modeling approach captured the salient influences of fiber amount and orientation on the fracture behavior of beam samples extracted from the slabs. In yet another study, spatially correlated random fields were used to position fibers within the computational domain, such that instances of fiber clumping could be approximated [73]. Regions of lower fiber content, caused by clumping elsewhere, can act as defects that promote larger crack openings and lower resistance to fracture localization.

Analyses and Results
For the sake of example, a unit cell containing a fiber dosage of V f = 0.2% is subjected to uniaxial load in the x-direction. As before, a fiber aspect ratio of 50 and a conductivity ratio of σ f /σ m = 1000 are assumed. Each load step is followed by calculations of steadystate potential flow, induced by a potential difference of ∆φ = 1 acting across the unit cell. The corresponding length changes of the matrix and fiber elements, due to the mechanical loading, cause changes in conductivity or resistivity (i.e., the inverse of conductivity). Figure 11 shows the change in specimen resistance under sinusoidal mechanical loading, where t is pseudo-time used to represent the stage of loading; ∆R is change in resistance of the composite material within the unit cell; R 0 is resistance in the unloaded state; and ε max is the maximum strain measured over length L. In typical applications, service loads are considered and their application is transient (e.g., due to vehicle movement over a bridge deck), such that concrete creep and fracture are not primary factors in the resistance measurements. If the change in resistance is based simply on the length changes of individual elements, the composite gauge factor is slightly less than unity. As seen in Figure 11, however, a composite gauge factor of unity is obtained if strain in the loading direction is used to define length change. Whereas R 0 for the discrete and semi-discrete approaches differed by a small amount, the normalized quantities plotted in Figure 11 are essentially the same. In actual applications, the sensitivity of fiber resistance to strain depends not only on axial length change, but also on the influences of Poisson's ratio and any piezoresistivity of the fiber material. These factors can be roughly considered by assigning a gauge factor η to the fiber material. In other words, the dependence of fiber resistivity on straining was amplified by a factor of η within the simulations, which produces changes in the composite gage factor ( Figure 11). As expected for this non-percolated state, however, changes in fiber gauge factor produce only small changes in composite gauge factor. The ratio of elastic moduli was assumed to be E f /E m ≈ 7, but the simulated resistance values are not sensitive to this parameter. The axial strains of the fiber are essentially compatible with those of the matrix over the fiber length, except near the fiber ends where shear lag behavior occurs.
Depending on the application area, larger fiber-matrix conductivity ratios and larger fiber aspect ratios need to be considered. The simulation of conductivity enhancements associated with fiber-to-fiber contacts (or near contacts) is also necessary for comprehensive analyses of fiber composite designs.

Background
Cementitious materials exhibit semi-fluid behavior during the initial several hours after casting, prior to significant degrees of cement hydration. Sizable volumetric reduction, termed generalized plastic shrinkage, is caused by two sequential processes: plastic settlement and capillary shrinkage. The former process is driven by the gravitational consolidation of cementitious particles and aggregate, accompanied with upward drainage of free water to the top surface of the member, which is usually called bleeding [74]. The latter process is a consequence of the rise of capillary pressure within the interstices between particles [75,76]. Pore water menisci are formed among solid particles near the top surface after the surface layer of bleed water is consumed by evaporation, as well as by the early stages of hydration. This produces capillary pressure and the associated shrinkage, which may result in cracking of the exposed surface.
The extent of plastic shrinkage is highly dependent on mixture composition and evaporation rate [77,78]. Composition of the mixture mainly affects the mechanical and moisture transport properties of concrete, which govern the maximum deformation. Evaporation rate generally determines duration of the plastic settlement period together with the composition-determined bleeding rate and the magnitude of capillary shrinkage. Thus, the problems caused by plastic shrinkage could be mitigated by controlling those factors. However, concrete mixture design is often governed by various other performance requirements. Rate of evaporation is also difficult to control in some circumstances, particularly for drying environments. The addition of fibers has thus been proposed as a means for controlling plastic shrinkage without compromising mechanical performance of concrete [2,3,[79][80][81][82][83].
While the concrete is in a plastic state, fibers contribute to consistency of the mixture and provide stiffness against contractive deformation. In addition, fibers may significantly modify the transport properties of the matrix. Fibers are believed to work as bleeding channels that assist the transport of pore water to the concrete top surface, replenishing bleed water that has been consumed by evaporation [80][81][82]84,85]. These bleeding channels may result from greater permeability of the fiber-matrix interface, particularly for the case of hydrophobic fiber materials. Through fiber additions, the occurrence of capillary shrinkage can therefore be significantly postponed, which provides additional time for strength development of the matrix. Other research suggests that fibers may increase tortuosity of the flow paths, thus decreasing the bleeding rate and the settlement rate [79,80]. Synthetic fibers, such as polypropylene, polyvinyl alcohol, and glass fibers, have been preferred in controlling plastic shrinkage cracking due to their large interfacial surface area as well as, in some cases, faster development of bonding strength.
Neglecting early-age hydration reactions, plastic settlement of cementitious material is analogous to the consolidation of soil. Terzaghi [86] proposed a one-dimensional smallstrain consolidation model for soil, which was refined and extended to three dimensions by Biot [87]. The correctness and capabilities of the Biot model were confirmed using dualfield (displacement and pore pressure) coupled equations considering mass conservation of water [88]. Tan et al. [89] modeled settlement along gravitational direction specially for plastic concrete. The model was later modified by Kwak and Ha [90,91] featuring thermal coupling and evaporation calculation to determine the initial time of capillary shrinkage. This dual-field approach was later implemented within 2D finite element models [92], which were validated through comparisons with experimental data. Combrinck [77] introduced fracture to the planar model using the same dual-field finite element method.
In this paper, the effects of individual fibers on the plastic cementitious matrix are analyzed via the VCLM. The proposed model extends the analysis of plastic settlement to three dimensions using the aforementioned dual-field approach, while integrating the semi-discrete fiber model to represent both the mechanical and transport properties of individual fibers. Since fluid water is the controlling medium of moisture transport before air-entry, the proposed settlement model is adequate and no capillary model is required for the research aims of this paper. Property development due to cementitious materials hydration is not considered, although the modeling framework can accommodate such behavior [40]. The model is first verified with results from an existing plastic settlement model. Thereafter, the effects of fibers on mechanical and transport properties of the composite are examined both separately and in combination via parametric analyses.

Formulation of the Matrix Phase Model
The modeling of the matrix phase, implemented within the VCLM, relies on several simplifying assumptions that are listed below. In future work, it may be possible to include the effects of early-age solidification associated with cementitious materials hydration [93][94][95]. • The load-deformation behavior of the mixture follows small-strain elastic theory. • No hydration takes place during the time period covered by the analyses, such that heat generation, chemical shrinkage, or hardening do not occur. • The material is isothermal and isolated from ambient conditions. • The surface of the formwork is considered to be frictionless and no wall effect is considered.
As a continuum, every point in the solid body of plastic concrete satisfies the equilibrium condition ∇σ + ρg = 0 (15) where ρ is the density of the two phase porous medium, taken as ρ = nS w ρ w + (1 − n)ρ s , in which ρ w is the density of water and ρ s is the density of solid skeleton; n is the porosity; S w (=1, herein) is degree of saturation and g is the acceleration of gravity. Stress tensor σ can expressed in terms of effective stress σ and pore water pressure p under isothermal conditions where Biot's coefficient α b = 1 − K d /K s , in which K d and K s are bulk moduli of the drained porous material and solid portion of the matrix, respectively. Because the solid itself can be treated as almost incompressible when compared to the drained mixture, we can assume α b = 1 here. For the VCLM, m = [1 0 0 0 0 0] T . Neglecting other strain sources, such as temperature variation, effective stress σ = Dε, where D is the stiffness matrix of drained solid skeleton. Transport of water in the porous media is driven by gravity and can be described using Darcy's Law in rate format. For the applications considered herein, Darcy's law is written in terms of hydraulic conductivity where q is the volumetric flux density, hydraulic conductivity K h = κγ w /µ, in which κ is the permeability of the porous media, µ is the viscosity of water, and γ w is the specific weight of water. The total mass of water of any phase stays constant during the transport process. For the case of settlement, only the liquid phase of water is considered, and the mass conservation equation of liquid water and solid skeleton together considering the mass diversion of hydrates [96] can be deduced: whereu is the velocity of the solid skeleton andṁ hydr is the production rate of new hydrates, which is assumed to be zero herein.
Equations (16) and (18) can be transformed into the following matrix form by applying a Galerkin residual method: (19) or in a simplified format where the elemental stiffness matrix K uu and diffusivity matrix K pp are defined according to Equations (4) and (11), respectively (with K h /γ w replacing D in Equation (11)). The coupling matrix K up and compressibility matrix C pp are where Q is the Biot compressibility parameter and d = 1, 2 or 3 for 1D, 2D, or 3D networks, respectively. Matrix N takes the following form: Initial conditions are assigned to the domain of interest Ω and its boundary Γ at time t = 0 as u = u 0 p = p 0 and Dirichlet boundary conditions are u =û on Γ u p =p on Γ p as well as Cauchy type with boundary traction and pore water flux After their transformation to the global coordinate system, the elemental matrices are assembled to form the system equations, which are solved using the generalized trapezoidal method where parameter θ = 0.5, which corresponds to the Crank-Nicolson method.

Verification of the Matrix Phase Model
The settlement component of the proposed model is verified with the numerical results of Kwak and Ha [92], who simulated settlement of a mortar specimen using planar finite elements. The virtual prismatic specimen has dimensions b × l × h equal to 800 × 800 × 400 mm. The upper surface of the specimen is exposed to the environment and is free to deform, while the other surfaces are constrained by a rigid mold with frictionless walls. As previously noted, moisture exchange with the environment and hydration of the cementitious materials are not considered in this simulation. The parameter values used for the simulations are indicated in Table 1. Hydraulic conductivity K h is considered to be a function of the initial hydraulic conductivity and porosity: in which porosity is assumed to evolve with volume change due to settlement [97] n = 1 + ( where n 0 is initial porosity, which can be calculated from the mixture design; ε v is the volumetric strain that is approximated as the sum of the principle strains, i.e., ε v ≈ ε 1 + ε 2 + ε 3 . Since the same constitutive model is applied in both modeling frameworks, the settlement computed using the VCLM should be comparable to that of the finite element method when applying identical material parameters (as given by Kwak and Ha [92]). The proposed model is examined with one-, two-, and three-dimensional Voronoi cell discretizations. The latter two cases are shown in Figure 12. The settlement is measured as the top-center node displacement for all cases. As seen in Figure 13, the VCLM results are practically identical to those of the twodimensional FE model, regardless of the dimensionality of the VCLM. Since the interface between the specimen and the mold is frictionless, settlement and pore pressure should be a function only of the vertical coordinate of a material point, provided the mesh density is adequate. Figure 14 displays the vertical displacement and pore pressure distributions over the specimen height at different times. The proposed settlement model produces consistent results regardless of the dimensionality of the VCLM. The smooth, curve-like nature of the distributions of nodal point values suggests the random discretization of the specimen does not introduce significant bias within the simulation results. Considering the close agreement between the 2D and 3D simulation results, 2D versions of the VCLM are used for the parametric analyses in the next section.

Influences of Short-Fiber Reinforcement: Parametric Study
The effects of straight, short fibers on the mechanical and transport properties of plastic cementitious material are assessed independently and in combination. The fibers are assumed to be 12 mm in length and 0.04 mm in diameter, and homogeneously distributed within the computational domain using coordinates and angles supplied by a uniform random number generator. The fibers are assumed to be hydrophobic, leading to enhanced porosity of the fiber-matrix interface. The effect of this enhanced porosity is lumped into the transport properties of the fibers, as described later. The elastic modulus of the fibers is chosen to be 40 GPa. It is further assumed that no relative movement (i.e., slip) occurs between the fibers and surrounding matrix. Fiber dosages up to 2% are considered.
The stiffening effect of fibers on plastic settlement is investigated first. Figure 15a shows the influence of fiber dosage on controlling plastic settlement. By increasing fiber dosage, the stiffness contributions of the fibers increase. As expected, higher volume fractions of fibers lead to less plastic settlement. Allowance for slip along the fiber matrix interface would reduce the magnitude of this trend. The cases of high fiber dosage also exhibit larger dispersion in the values of plastic settlement for any given nodal point elevation, as shown in Figure 16. This is reasonable considering the local stiffening effects of individual fibers, which is a source of heterogeneity.  The potential influence of enhanced transport, due to porosity of the fiber-matrix interface, on plastic settlement is investigated next. As noted above, the effective transport properties of the interface are considered by modifying the transport properties of the fibers. For this case, the Voronoi cell sizes are large relative to the fiber length. Accordingly, to account for the local influences of each fiber on the flow field, effective conductivity of the fiber depends on the location of the facet intersection point along the fiber length [98,99]. Conductivity grows from the fiber ends, similar to the shear lag effect present in the mechanical representations of the fibers [21]. In an average sense, the conductivity of the fibers was assumed to be greater than that of the matrix by a factor of about 10 (i.e., σ f /σ m ≈ 10). The results shown in Figure 15b support some of the experimental observations that fiber additions may increase plastic settlement, along with concurrent increases in the amount of bleed water. Figures 17 and 18 show how the distribution of pore pressure over the specimen height is potentially affected by conductivity of the fiber-matrix interface. Within increasing solution time, pore pressure diminishes over the specimen height. With the addition of fibers, however, the rate of pore pressure decrease is higher and the shape of the pore pressure profile also differs. Here, too, the introduction of locally enhanced transport properties, due to the presence of fibers, causes greater dispersion of the pore pressure values at any given nodal elevation.  The combined effect of both stiffness and conductivity enhancements, associated with fiber additions, is considered next. The relative influence of each form of enhancement depends on time after casting, as seen in Figure 15c. Initially, the rate of settlement is influenced mainly by the changes in conductivity associated with fiber dosage. Higher fiber dosages lead to higher conductivity of the composite material and thus higher rates of settlement. Ultimately, however, the stiffening effects of the fibers limit the amount of plastic settlement. It can be concluded that fiber additions may reduce the severity of plastic shrinkage phenomena by providing higher amounts of bleeding water during the early development period and reducing the maximum potential settlement. The competing influences of the stiffening effect and conductivity modification associated with fiber additions might explain, at least in part, the conflicting results of some experimental studies of plastic settlement.
The constitutive law based on Darcian flow is applicable not only for the plastic settlement period, but also for the period of early capillary pressure build-up until air-entry takes place since the water phase within the material is still continuous. The proposed model is therefore potentially useful for studying the onset of plastic shrinkage cracking of cementitious materials, including those containing short-fiber reinforcement.

Conclusions
This research has produced semi-discrete models of fiber-assisted transport in cementbased materials. Voronoi cell lattice models (VCLM) serve as a framework for the proposed developments. The models are used, in combination with analogous models of the mechanical behavior of short-fiber reinforcement, to simulate multi-field phenomena. The proposed models are verified through comparisons with theory and results based on discrete (volumetric) discretizations of the fibers. Capabilities of the semi-discrete approach are demonstrated through electro-mechanical analyses, with applications to self-sensing, and hygro-mechanical analyses of plastic settlement. Computational efficiency of the semi-discrete approach allows for larger amounts of fibers, or larger spatial dimensions, to be considered, which is important for structural applications. The following conclusions can be made:

1.
At low fiber dosages, both the discrete and semi-discrete approaches to modeling conductive fibers offer precise estimates of conductivity of the composite material, in comparison with the equivalent inclusion method (EIM) of Hatta and Taya.

2.
The semi-discrete approach is computationally efficient in its representation of the conductivity contributions of individual fibers. Fiber additions do not increase the number of computational degrees of freedom of the model. For V f = 0.01, which is one of the fiber dosages considered herein, the number of computational nodes was nearly one order of magnitude smaller, relative to the number of nodes used for modeling the composite with discrete (volumetric) fibers. The computational savings would be even more dramatic when considering higher fiber dosages or aspect ratios.

3.
When comparing with the EIM results, the semi-discrete approach requires sufficient resolution of the background lattice such that the fibers are electrically isolated from one another by matrix material. Increasing V f , while keeping nodal density constant, eventually causes composite conductivity estimates of the semi-discrete model to be higher than the corresponding EIM values.

4.
For a given fiber content, conductivity of the composite varies with directionality of the fiber distribution. The discrete and semi-discrete models capture the expected result: conductivity of the composite material increases with preferential orientation of fibers with respect to the electrical field. This capability should enable the study of fiber distribution effects in actual members, where boundary effects and other sources of fiber distribution non-uniformity are present.

5.
The coupled electro-mechanical model captures the changes in composite resistance due to axial loading. An updating of nodal point positions leads to a composite gauge factor of unity. Other composite gauge factors can be simulated by introducing a parameter associated with the sensitivity of fiber conductivity to applied strain. 6.
With respect to the simulations of plastic settlement of an unreinforced matrix, the results are in agreement with previously published results based on finite element analysis. Moreover, the results are shown to be equivalent, regardless of the dimensionality of the VCLM. 7.
The proposed model enables study of the relative benefits of introducing short fibers to reduce plastic settlement of concrete. Both the mechanical and transport modifications of the material are considered. The simulation results support some experimental observations that higher fiber content leads to reduced settlement and increased amounts of bleed water due to the stiffness and conductivity enhancements of the fibers, respectively. 8.
Model extensions and validation efforts, through comparisons with physical test data, are needed to address a variety of practical simulation needs. With respect to the validation efforts, it is difficult to measure some of the model input parameters (e.g., the mechanical and hydraulic conductivity properties of individual fibers within fresh concrete). Such quantities need to be determined, either experimentally or through inverse calculations.
Author Contributions: J.K. and Y.P. conceptualized the models and led early program development; S.I. and J.E.B. continued program development; Y.P., J.K., S.I., and J.E.B. conducted the analyses; Y.P. and J.E.B. wrote the paper with contributions from J.K. and S.I. All authors have read and agreed to the published version of the manuscript.
Funding: This research received no external funding.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript:

EIM Equivalent inclusion method RBSN
Rigid-body-spring network VCLM Voronoi-cell lattice model

Appendix A
According to the equivalent inclusion method of Hatta and Taya [14], conductivity of a 3D composite material containing uniformly at random fiber inclusions is where σ m is conductivity of the matrix material, σ f is conductivity of the fiber material, and V f is the volume fraction of fibers. In addition, in which tensor S ij depends only on the shape of the ellipsoidal inhomogeneity representing the fibers.