Stochastic Constitutive Model of Isotropic Thin Fiber Networks Based on Stochastic Volume Elements

Thin fiber networks are widely represented in nature and can be found in man-made materials such as paper and packaging. The strength of such materials is an intricate subject due to inherited randomness and size-dependencies. Direct fiber-level numerical simulations can provide insights into the role of the constitutive components of such networks, their morphology, and arrangements on the strength of the products made of them. However, direct mechanical simulation of randomly generated large and thin fiber networks is characterized by overwhelming computational costs. Herein, a stochastic constitutive model for predicting the random mechanical response of isotropic thin fiber networks of arbitrary size is presented. The model is based on stochastic volume elements (SVEs) with SVE size-specific deterministic and stochastic constitutive law parameters. The randomness in the network is described by the spatial fields of the uniaxial strain and strength to failure, formulated using multivariate kernel functions and approximate univariate probability density functions. The proposed stochastic continuum approach shows good agreement when compared to direct numerical simulation with respect to mechanical response. Furthermore, strain localization patterns matched the one observed in direct simulations, which suggests an accurate prediction of the failure location. This work demonstrates that the proposed stochastic constitutive model can be used to predict the response of random isotropic fiber networks of arbitrary size.


Introduction
Materials are characterized by certain degrees of random variations in their mechanical properties. This is true for all materials but especially pronounced in disordered materials such as thin fiber networks [1]. Such variations can be the cause of unexplained occasional failures that cannot be predicted by deterministic material models [2][3][4]. It is, therefore, crucial to develop a sound stochastic approach in studying mechanical failure of thin fiber networks of arbitrary size.
The rapid development of characterization tools enables the quantification of randomness at different scales and the construction of random realizations of a fiber network. In our previous works, the mechanical behavior of randomly generated networks was investigated using detailed direct micromechanical simulations [5][6][7][8] (see Section 2.1). However, although such simulations can capture the complicated mechanisms of failure, they cannot yet be employed for product development due to the overwhelming computational costs required to capture the relevant product sizes. For instance, simulation of the uniaxial mechanical response up to strain localization and failure of a 24 mm × 24 mm fiber network takes two days on a modern 28-core, 128 GB RAM supercomputer. Therefore, the straightforward use of direct micromechanical simulations for predicting the mechanical response of large structures of complex geometry will be impossible in the foreseeable future. In Figure 1, examples of paper-based products composed of random fiber networks are shown. A micro-tomography of a typical fiber network (Figure 1b) followed by numerical reconstruction (Figure 1c) is also shown. An example of random occurrence of breaks of paper-based products can be found in paper-making machines. Every break costs around 6000 € and some paper machines can experience up to six breaks a day. With the speed of the machine reaching 2000 m/min, the break can, therefore, occur once per 480 km of produced paper, which makes it an extremely stochastic event.
The aforementioned difficulties associated with predicting the mechanical response of larger structures based on random local material properties has been the topic of extensive research for a wide variety of materials including fiber systems [9][10][11][12][13][14][15][16][17]. In recent works, different methods based on stochastic volume elements (SVEs) [18] have been used in a variety of applications to alleviate the computational burden. In Reference [19], microstructure-property relations of heterogeneous materials were assessed using a hierarchical decomposition of statistically representative volume elements (RVEs) [20] into smaller SVEs. In Reference [21], analysis of particle reinforced viscoelastic polymer nanocomposites was also performed using SVEs. Similarly, a stochastic multiscale model for polycrystalline materials based on SVEs was developed in Reference [22]. In general, three different scales are involved in such analysis: (1) the micro-scale, which is a characteristic size of the microstructure; (2) the mesoscale, an intermediate scale such as the size of the volume element over which a homogenization can be performed; and (3) the macro-scale which is the size of the structural problem.
In this work, a stochastic constitutive model is proposed to address the aforementioned computational difficulties associated with modeling large random fiber networks. The model is based on SVEs and random field generation of local material properties. It includes three major components: a deterministic constitutive model, characterization of SVE size-dependent model parameters, and generation of two correlated non-Gaussian random spatial fields describing local material properties. The model is numerically validated by ensuring (1) accurate prediction of mechanical failure, i.e., strength and strain to failure as well as strain localization pattern; and (2) generation of continuum random realizations that are statistically equivalent to the ones generated using a direct simulation of fiber networks. supercomputer. Therefore, the straightforward use of direct micromechanical simulations for predicting the mechanical response of large structures of complex geometry will be impossible in the foreseeable future. In Figure 1, examples of paper-based products composed of random fiber networks are shown. A micro-tomography of a typical fiber network ( Figure 1b) followed by numerical reconstruction (Figure 1c) is also shown. An example of random occurrence of breaks of paper-based products can be found in paper-making machines. Every break costs around 6000€ and some paper machines can experience up to six breaks a day. With the speed of the machine reaching 2000 m/min, the break can, therefore, occur once per 480 km of produced paper, which makes it an extremely stochastic event.
The aforementioned difficulties associated with predicting the mechanical response of larger structures based on random local material properties has been the topic of extensive research for a wide variety of materials including fiber systems [9][10][11][12][13][14][15][16][17]. In recent works, different methods based on stochastic volume elements (SVEs) [18] have been used in a variety of applications to alleviate the computational burden. In Reference [19], microstructure-property relations of heterogeneous materials were assessed using a hierarchical decomposition of statistically representative volume elements (RVEs) [20] into smaller SVEs. In Reference [21], analysis of particle reinforced viscoelastic polymer nanocomposites was also performed using SVEs. Similarly, a stochastic multiscale model for polycrystalline materials based on SVEs was developed in Reference [22]. In general, three different scales are involved in such analysis: 1) the micro-scale, which is a characteristic size of the microstructure; 2) the mesoscale, an intermediate scale such as the size of the volume element over which a homogenization can be performed; and 3) the macro-scale which is the size of the structural problem.
In this work, a stochastic constitutive model is proposed to address the aforementioned computational difficulties associated with modeling large random fiber networks. The model is based on SVEs and random field generation of local material properties. It includes three major components: a deterministic constitutive model, characterization of SVE size-dependent model parameters, and generation of two correlated non-Gaussian random spatial fields describing local material properties. The model is numerically validated by ensuring 1) accurate prediction of mechanical failure, i.e., strength and strain to failure as well as strain localization pattern; and 2) generation of continuum random realizations that are statistically equivalent to the ones generated using a direct simulation of fiber networks.  The methodology proposed in this work is generic but demonstrated for isotropic fiber networks subjected to uniaxial loading until failure. The paper is organized as follows. An introduction to the direct fiber network simulation as well as important concepts in stochastic multiscale modeling and random spatial field modeling are presented in Section 2. The proposed stochastic constitutive model is presented in Section 3. In Section 4, a validation of the model is first performed for small specimens by comparing the results to direct fiber network simulations. The approach is thereafter demonstrated by random generation and mechanical simulation of large continuum models of fiber networks.

Direct Simulation of a Thin Fiber Network
A 3D network of fibers is generated using a deposition technique in which the fibers are sequentially deposited on a flat surface from two sides. The deposition algorithm can be outlined as follows: The fiber geometry is chosen from the fiber characterization data acquired with FiberLab (Valmet Fiber Image Analyzer), which is an apparatus for measuring fiber characteristics. It contains length, width, height, wall thickness, and curvature. The curvature is represented through an arc of constant curvature located in a single plane parallel to the deposition plane. The cross-sectional data is corrected using microtomography scans of a paper produced using the considered softwood kraft pulp. The details of the correction are described elsewhere [6].

2.
The fiber orientation is chosen randomly in this work although it can be controlled to match a specific distribution.

3.
The fiber position before deposition onto the domain is chosen randomly. 4.
The first fibers are deposited on the flat plane consecutively from above or below. For the subsequent fibers, we first find the intersection between them and the previously deposited fibers in the plane (Figure 2a). 5.
The found intersection points are lifted discretely to exclude penetration ( Figure 2b). The contact search diameter depicted in the figure corresponds to the height of the fiber, which is smaller or equal to the width of the fiber normally. 6.
The fiber geometry is smoothed to remove discontinuities caused by the previous step ( Figure 2c). During the smoothing, we control the maximum angle the fibers can form, and it was set to 5 degrees in this study. 7.
When the grammage (the basis weight or the weight per unit area) of the network has reached the prescribed value, the deposition procedure is stopped. The grammage used in this work was 28 g/m 2 , which is relatively low but corresponds to the set of handsheets used to calibrate the measurements in Reference [6]. 8.
The thickness of the network is evaluated and measured using the procedure described in Reference [6]. Thereafter, the thickness is brought to the target value by uniform scaling of the coordinates in the thickness direction with respect to the center plane of the network. The target value in this study was 68 micrometers, which corresponds to the measured value used in the calibration [6]. The scaling may result in interpenetration, which are zeroed out during the subsequent computations.
After generation of the geometry, the fibers are represented with a fine mesh of curved segments and imported into a custom finite element code. The fibers are converted into cubic splines and a full 3D network finite element model [5,6,23] is used, in which fibers are resolved as a chain of 3D quadratic Timoshenko/Reissner beam elements [24] with six degrees of freedom at each node. The implicit time-integration is used. During contact detection, cross-sections of the fibers are treated as circular and rigid, with the diameter equal to the mean of the values of fiber width and height. The mechanical Materials 2019, 12, 538 4 of 28 bonds behavior is described with traction-separation laws with a cohesive zone model based on the contact forces [6].  Table 1 summarizes the details of the fibers used in the analysis. The underlying data is identical to that used by Reference [6], who focused on the mean values of the tensile properties of the networks. The constitutive response of the fibers is described by a bilinear plasticity model with the material parameters as listed in Table 2. Table 1. Fiber geometrical data used in the network simulation, based on the direct measurement on wet pulp (FMA), on dry sheets (µCT), and numerical parameters in terms of length-weighted mean and standard deviation (SD) values.  1 The ratio between dry and wet measured radius and wall thickness, respectively. The contact conditions at the fiber bonds are governed by a bilinear cohesive tractionseparation law. It requires the definitions of a bond's stiffness, a bond's strength and separation, see Table 3. Since the beam's cross-section is rigid against the local normal forces and shear forces, the physical compliance of the fiber at the bonding sites is represented solely with the stiffness of the penalty-based contact element. The selection of the parameters listed in Table 2 and Table 3 are based on the experimental calibration in the tensile test described in Reference [6].

Mean SD
The software used for the fiber network reconstruction as well as the input data are uploaded as a Supplementary Material to the current paper.  Table 1 summarizes the details of the fibers used in the analysis. The underlying data is identical to that used by Reference [6], who focused on the mean values of the tensile properties of the networks. The constitutive response of the fibers is described by a bilinear plasticity model with the material parameters as listed in Table 2. Table 1. Fiber geometrical data used in the network simulation, based on the direct measurement on wet pulp (FMA), on dry sheets (µCT), and numerical parameters in terms of length-weighted mean and standard deviation (SD) values.  1 The ratio between dry and wet measured radius and wall thickness, respectively. The contact conditions at the fiber bonds are governed by a bilinear cohesive traction-separation law. It requires the definitions of a bond's stiffness, a bond's strength and separation, see Table 3. Since the beam's cross-section is rigid against the local normal forces and shear forces, the physical compliance of the fiber at the bonding sites is represented solely with the stiffness of the penalty-based contact element. The selection of the parameters listed in Tables 2 and 3 are based on the experimental calibration in the tensile test described in Reference [6]. The software used for the fiber network reconstruction as well as the input data are uploaded as a Supplementary Materials to the current paper.

Multi-Scale Modeling
Multi-scale mechanical modeling of thin fiber networks is a powerful technique aimed at predicting the mechanical response at the macro-mechanical level, based on properties at the meso-and microscale. The up-scaling methodology in multi-scale methods generally relies on either SVEs [18] or RVEs [25][26][27][28][29]. In the former, elements are characterized by a stochastic variation in their mechanical response. Put another way, SVEs taken from different locations within the macroscopic body have different material properties due to random variation at the micromechanical level. This is opposed to the concept of RVEs, for which the material properties or any chosen statistical descriptor remain constant regardless of the location of the RVEs within the macroscopic body [20]. An RVE can, therefore, be defined as the smallest volume over which statistical representation can be made for the considered material property [29]. It is used to determine the corresponding effective properties for a homogenized macroscopic model. According to the micro-meso-macro principle [30], the RVE should be sufficiently large in order to contain representative information about the microstructure, but its size should be small compared to the studied macroscopic body.
Multi-scale strategies aiming at reconstructing a macroscopic body using RVEs may not be successful, since an RVE may not exist [29]. An attempt to find an appropriate size of the RVE can be made by increasing the size of an SVE until its response becomes independent of its spatial location. In statistical terms, an appropriate RVE size is considered to be found when the standard deviation of the chosen material property vanishes, and its mean value converges to a constant value. It has been shown that, for instance, when considering a statistical descriptor related to the softening behavior of the material in question, its mean value does not converge with increasing element size [5]. In such cases, a multi-scale model based on the response of SVEs may be more appropriate.
As opposed to multi-scale applications where the focus is solely on the linear elastic response, prediction of strength to failure of thin fiber networks is the main target of the present work. In Figure 3, a multi-scale model for thin fiber networks based on SVEs is described. A fiber network is spatially discretized using a chosen SVE size, see Figure 3a. It should be noted that the discretization may involve overlapping SVEs. In Figure 3b, the response from direct numerical simulations on SVEs from different spatial locations is shown and the random variation is apparent. In order to correctly model the transition between different scales, the fiber network is created using a statistically equivalent microstructure, Figure 3c. However, direct numerical simulation on a large fiber network in Figure 3a is not possible due to formidable computational cost. The aim of this multi-scale strategy is therefore to predict the mechanical response of the fiber network based on the mechanical response of the SVEs in Figure 3b. made by increasing the size of an SVE until its response becomes independent of its spatial location. In statistical terms, an appropriate RVE size is considered to be found when the standard deviation of the chosen material property vanishes, and its mean value converges to a constant value. It has been shown that, for instance, when considering a statistical descriptor related to the softening behavior of the material in question, its mean value does not converge with increasing element size [5]. In such cases, a multi-scale model based on the response of SVEs may be more appropriate.

Random Field Representation
Accurate models of random local spatial variations of material response, such as the one observed in Figure 3b, are necessary to correctly predict mechanical failure of the disordered fiber network. In order to describe these local variations, random spatial fields are needed. A random field can be represented as an infinite expansion of orthogonal basis functions and random expansion coefficients. Let h(r, ω) : D × Ω → R be such a random continuous spatial field defined over a spatial domain D, where r = x y z T is the spatial vector and ω ∈ Ω denotes an element of the sample space indicating that the involved quantity is random. The random field's covariance function is denoted as C(r, r ) = cov(h(r, ω), h(r , ω)). By the definition of a covariance function, it has the spectral decomposition [31] where λ i and η i (r) are the eigenvalues and eigenfunctions of the covariance function, respectively. That is, they are the solution to the Fredholm integral The Karhunen-Loeve (KL) expansion [31] of the random field h(r, ω) uses the eigenfunctions of the covariance function as expansion basis according to where h(r) is the mean function and ξ i (ω) are orthogonal random variables with zero mean and unit variance. Explicit expression for ξ i (ω) can be found by multiplying Equation (3) by η i (r) and integrating over the spatial domain D according to where use is made of the orthogonality property of the set {η i (r)} ∞ i=1 . The choice of basis function and expansion coefficients in the KL expansion is motivated by truncating the infinite expansion according to Equation (3) and minimizing the total truncation error. If the expansion is truncated at the M-th term, the truncation error is given by . It can be shown that the total mean square error D e M (r) 2 dr is minimized (subject to the orthogonality condition of {η i (r)}) if the Fredholm integral equation according to Equation (2) is satisfied. Therefore, of all possible choices of expansion basis, the Karhunen-Loeve (KL) expansion [31] is the best approximation of the original random field in the sense that it minimizes the total mean-square error resulting of its truncation.
The usefulness of the Karhunen-Loeve expansion hinges on the ability to determine the eigenvalues and eigenfunctions of the random field covariance function through the Fredholm integral equation. In this work, the continuous random field is discretized at N spatial points r 1 , r 2 , . . ., r N in the domain D. An N × N covariance matrix with elements C pq = C r p , r q can be defined, where r p and r q are two spatial locations. For this discrete and finite case, the Fredholm integral equation can be rewritten as are eigenvalues and eigenvectors of the covariance matrix C. That is, when the random field is discretized, operations on functions are transformed into operation on matrices. The KL expansion of the discrete random field [32][33][34][35] can be written as where denotes a random vector whose elements are the random field values at the N discrete points and The random coefficients are obtained from In general, the underlying covariance function C(r, r ) is not known. One crucial task is therefore the modelling of the covariance function based on sampling and characterization of the spatial random fields. It should be noted that the above description can be extended to represent multivariate random spatial fields. In this case, the covariance matrix C is composed of both auto-and cross-covariances, where the latter describes the correlation between different spatial fields, see Section 3.2.

Overview of Methodology
The aim of this work is to model computationally expensive random isotropic 3D-fiber network realizations by statistically equivalent random continuum realizations, see Figure 4. The model is based on the random spatial fields of strain to failure ε f (r, ω) and strength σ f (r, ω) which are found using direct mechanical simulation on SVEs. The latter is defined as the ratio of the applied force on the SVE and the average thickness of the whole fiber network. The characterization and reconstruction of the random spatial fields as well as the stochastic constitutive model used to predict the mechanical response of isotropic fiber networks of arbitrary size, is presented in Sections 3.2-3.5.
denotes a random vector whose elements are the random field values at the discrete points and = ℎ ( ) ℎ ( ) ⋯ ℎ ( ) . The random coefficients are obtained from In general, the underlying covariance function ( , ′) is not known. One crucial task is therefore the modelling of the covariance function based on sampling and characterization of the spatial random fields. It should be noted that the above description can be extended to represent multivariate random spatial fields. In this case, the covariance matrix is composed of both autoand cross-covariances, where the latter describes the correlation between different spatial fields, see Section 3.2.

Overview of Methodology
The aim of this work is to model computationally expensive random isotropic 3D-fiber network realizations by statistically equivalent random continuum realizations, see Figure 4. The model is

Characterization and Simulation of Random Spatial Fields of Strength and Strain to Failure Based on Stochastic Volume Elements (SVEs)
The proposed continuum model of isotropic 3D fiber networks is based on two continuous macro-mechanical spatial fields describing the material properties at each material point in a 2D spatial domain r = x y T . These are the uniaxial strength σ f (r, ω) and uniaxial strain to failure ε f (r, ω). The characterization methodology of the random spatial fields is detailed in Table 4.
The first step is the sampling of the spatial fields. Since the fiber network is isotropic, it is sufficient to determine the spatial fields in one coordinate direction, see Figure 5a. A specimen of length L is discretized using SVEs of size L SVE × W SVE . The discretization distance ∆L is the distance between the centers of two neighboring SVEs corresponding to material points in the continuum model. All SVEs are cut from the fiber network and a numerical uniaxial test is performed. The resulting one-dimensional spatial fields of the strain to failure ε f (x, ω) and strength σ f (x, ω) are schematically shown in Figure 5b. An important characteristic of these spatial fields is their mutual correlation. It is computed from the Pearson's correlation coefficient based on the one-dimensional spatial fields in Figure 5b according to In Equation (9), superscript "t" denotes that this is a target value that needs to be satisfied when simulating new realization of ε f (x, ω) and σ f (x, ω) or ε f (r, ω) and σ f (r, ω). Table 4. Characterization and simulation methodology.

Simulation
Simulation using Karhunen-Loeve expansion of 2D correlated Gaussian fields: g 1 (r, ω) and g 2 (r, ω) Equation (21) Transformation to original random field: ε f (r, ω) and σ f (r, ω) Equation (22) The next step is to fit a probability density function for the random variables ε f (ω) and σ f (ω). In Figure 5c, histograms of m random outcomes of ε fi and σ fi are shown. These are generated by cutting SVEs from random spatial locations in a large fiber network. The probability density functions where ϕ is the standard normal density function and h ε f , h σ f > 0 are smoothing bandwidth. The next step is to fit a probability density function for the random variables ( ) and ( ). In Figure 5c, histograms of random outcomes of and are shown. These are generated by cutting SVEs from random spatial locations in a large fiber network. The probability density functions (pdf) ( ) and ( ) are approximated based on the kernel density estimators where is the standard normal density function and ℎ , ℎ > 0 are smoothing bandwidth. Thereafter, a transformation from the non-Gaussian fields ( ( , ) and ( , )) to correlated Gaussian random fields ( ( , ) and ( , )) with zero mean and unit variance is performed. It should be emphasized that this step is performed to facilitate the modeling of the covariance function and that the Gaussian fields are later transformed back to the original spatial fields. The transformation to Gaussian random fields is performed using the relation [37,38] ( , ) = ( , ) where ( ) and ( ) are the cumulative density functions (CDFs) found by integration of the kernel density estimators according to Equation (10), and is the standard normal CDF. In Figure  5d), these Gaussian random fields are schematically shown. An important characteristic of each of these one-dimensional Gaussian fields is the average distance between two zero level up-crossings denoted by [36]. These are computed based on the number of up-crossings in the interval [0, ] corresponding to the specimen length in Figure 5a, defined as Thereafter, a transformation from the non-Gaussian fields (ε f (x, ω) and σ f (x, ω)) to correlated Gaussian random fields (g 1 (x, ω) and g 2 (x, ω)) with zero mean and unit variance is performed. It should be emphasized that this step is performed to facilitate the modeling of the covariance function and that the Gaussian fields are later transformed back to the original spatial fields. The transformation to Gaussian random fields is performed using the relation [36,37] where F ε f (ε f ) and F σ f (σ f ) are the cumulative density functions (CDFs) found by integration of the kernel density estimators according to Equation (10), and φ is the standard normal CDF. In Figure 5d), these Gaussian random fields are schematically shown. An important characteristic of each of these one-dimensional Gaussian fields g i is the average distance between two zero level up-crossings denoted by µ t upcri [38]. These are computed based on the number of up-crossings in the interval [0, L] corresponding to the specimen length in Figure 5a, defined as The average distance between two zero level up-crossings for g 1 and g 2 is therefore given by From Equation (13) it is seen that the required specimen length L has to be chosen sufficiently large so that both µ t upcr1 and µ t upcr2 converges to a constant value, see Section 4.1 and Figure 11c,d. These values are target values in the sense that they need to be satisfied when simulating new realization of ε f (x, ω) and σ f (x, ω) followed by the transformation according to Equation (11).
A crucial step in the characterization of random spatial fields is to choose an appropriate model for the underlying covariance function. This is easiest performed in the transformed Gaussian random variable space. The choice of model is important in that the properties of the Gaussian random field, such as smoothness or differentiability, are determined by the eigenvalues and eigenfunctions of the covariance function. In this work, the squared exponential kernel is used to model the auto-covariance for g 1 (r, ω) and g 2 (r, ω), denoted by K 11 (r, r ) and K 22 (r, r ), respectively, according to This choice of covariance function satisfies the condition of high smoothness and is infinitely differentiable. From Equation (14), it can be seen that the auto-covariance is only a function of the distance between two spatial locations |r − r | in relation to a characteristic length scale ( 1 and 2), see Figure 5d. The choice of length scale therefore affects how far away two spatial field values need to be in r-space before they can be significantly different. Furthermore, the spatial fields g 1 (r, ω) and g 2 (r, ω), representing the normalized strain to failure and strength, are assumed spatially correlated. Therefore, the cross-covariance function cov(g 1 (r, ω), g 2 (r , ω)) needs to be modeled to describe this spatial correlation. The cross-covariance is modeled as i.e., using an exponential kernel function with a correlation coefficient ρ 12 and characteristic length 12 . A total of four unknowns 1 , 2 , ρ 12 and 12 , need to be determined based on the one dimensional sample spatial field data in Figure 5d. For a discretized random field, the multivariate covariance matrix is given by where K 11,pq = K 11 r p , r q , K 22,pq = K 22 r p , r q , and K 12,pq = K 12 r p , r q . Using the one-dimensional spatial field data, see Figure 5d, these matrices become N × N covariance matrices. The sample data from both random fields can be collected into the vector where This set of observations g(ω) can be regarded as samples from a multivariate Gaussian distribution, i.e., g(ω) ∼ N (0, K). The log-marginal likelihood function [38] (MLII) is the natural logarithm of the distribution of g(ω), which can be written as A widely used method for the determination of unknown parameters of the Gaussian random fields is the maximization of the log-marginal likelihood function. In this work the unknown parameters 1 , 2 , ρ 12 , and 12 are determined by the maximization of the log-marginal likelihood function according to Equation (18) subject to three constraints. This can equivalently be written as a minimization problem according to The first constraint enforces that correlation coefficient between the strength and strain to failure ρ ε f σ f matches the one computed from the analysis of the L × W SVE specimen in Figure 5d, where ρ t ε f σ f is given by Equation (9). The second and third constraints enforces that the average distance between two zero levels up-crossings for the normalized spatial fields g 1 (x, ω) and g 2 (x, ω), respectively, matches the target ones, µ t upcr1 and µ t upcr2 , computed in Equation (13). Successful solution of the optimization problem according to Equation (19) is simplified if good initial guesses for the unknown parameters are used. The following initial values are used in the first iteration: . (20) As is seen, the initial guesses for the characteristic length scales 2 are proportional to the target average distance between zero-level up-crossings. The proportionality constant 2π is based on the analysis of univariate one-dimensional spatial fields with an underlying exponential Gaussian covariance function [38]. Furthermore, the initial guess ρ (0) 12 is given by ρ t ε f σ f . It should be observed that ρ 12 = ρ ε f σ f due to the non-linear transformation that relates the normalized spatial fields and the original spatial fields. The initial guess (0) 12 is chosen arbitrarily. In order to generate new realizations of 2D-spatial fields of strain to failure ε f (r, ω) and strength σ f (r, ω), simulation of Gaussian random fields g 1 (r, ω) and g 2 (r, ω) is first performed. The KL-expansion is used to simulate discrete and correlated Gaussian random fields g(ω) = g 1 (ω) T g 2 (ω) T T , see Section 2.3. It can be written as are eigenvalues and eigenvectors of the covariance matrix K in Equation (16). Since g(ω) ∼ N (0, K), the random variables ξ i (ω) become independent standard normal variables. A transformation is thereafter performed according to

Constitutive Model
An empirical relation between the equivalent stress σ eq and the equivalent strain ε eq of the form is assumed to model the response of each SVE in the network. In Equation (23), κ > 0 is a parameter controlling the elastic and elastoplastic modulus, n > 0 is a damage exponent and f 1−2 > 0. The relation is a variant of a Bammann-Chiesa-Johnson (BCJ) model previously used in stochastic constitutive modeling [39]. Enforcing that σ eq ε eq = ε f = σ f and ∂σ eq ∂ε eq ε eq =ε f = 0 , yields the following relations The elastic modulus can be written as Since f 1−2 are functions of n and κ, it is sufficient to study the effect of these two parameters in Equation (23). As can be seen from Figure 6, increasing the damage exponent n or decreasing the parameter κ has a similar effect on the stress-strain curve. In the proposed model, one value of the damage exponent n is assumed suitable to model the response of SVEs of same size, while κ is fitted to each SVE realization. The residual sum of squares resulting from fitting Equation (23) to a uniaxial test performed on the j th SVE is denoted by e j n, κ j . The optimal parameter n is found by minimizing the total sum of these errors for all SVEs, i.e., where N SVE is the total number of SVEs and the parameters κ * j N SVE j=1 are found by minimizing the error for each uniaxial test, as is seen from the constraints in Equation (26). The procedure is outlined in Figure 7. A large number of SVEs of the same size are randomly cut from a large fiber network, see Figure 7a. This assures that the SVEs are uncorrelated. A numerical uniaxial tensile test is performed on all extracted SVEs, see Figure 7b. All curves are scaled by their respective strength and strain to failure, see Figure 7c, and the relation in Equations (23) and (24) is fitted to the responses according to Equation (26). From Figure 6b, it is apparent that increasing the parameter κ, for a constant n (i.e., constant SVE size), also corresponds to an increasing elastic and plastic work in a uniaxial test. It is therefore clear that larger κ values are typical for responses characterized by larger strain to failure values. This leads to the assumption κ ∝ ε f , see Figure 7d and Section 4.1. A stochastic relation of the form is assumed, where c 1 and c 2 are found by least squares fit for SVEs of the same size and the residual error R ∼ N (0, s R ) is assumed to follow a normal distribution with a standard deviation s R computed from the residual fitting errors R j N SVE j=1 , see Figure 7d. The last term in Equation (27) ensures that the elastic modulus according to Equation (25) is not solely dependent on σ f and ε f but involves an additional random component. The uncertainty can be seen as an epistemic model error, resulting from the assumption according to Equation (27).
where is the total number of SVEs and the parameters are found by minimizing the error for each uniaxial test, as is seen from the constraints in Equation (26). The procedure is outlined in Figure 7. A large number of SVEs of the same size are randomly cut from a large fiber network, see Figure 7a. This assures that the SVEs are uncorrelated. A numerical uniaxial tensile test is performed on all extracted SVEs, see Figure 7b. All curves are scaled by their respective strength and strain to failure, see Figure 7c, and the relation in Equations (23) and (24) is fitted to the responses according to Equation (26).  From Figure 6b, it is apparent that increasing the parameter , for a constant (i.e., constant SVE size), also corresponds to an increasing elastic and plastic work in a uniaxial test. It is therefore clear that larger values are typical for responses characterized by larger strain to failure values. This leads to the assumption ∝ , see Figure 7d and Section 4.1. A stochastic relation of the form is assumed, where and are found by least squares fit for SVEs of the same size and the residual error ~(0, ) is assumed to follow a normal distribution with a standard deviation The softening part of the SVE response is modeled as where ξ > 0 is a parameter controlling the softening rate of the SVE response, see Figure 6. A small softening parameter ξ = 0.15 is chosen in order to increase the compliance of the system and assure numerical stability, but large enough for strain localization to onset in a large sample composed of SVEs. An isotropic plane stress state is assumed to prevail, and a non-linear hardening model is used to model the constitutive response of the thin fiber network. The relation between the increment in stress vector σ = σ x σ y τ xy T and strain vector ε = ε x ε y γ xy T is given by where D tan = D − D pl is the tangent stiffness matrix. The elastic stiffness matrix is given by with Poisson's ratio ν = 0.3 and elastic modulus given from Equation (25). The elastic-plastic tangent stiffness matrix is expressed as [40] where σ eq is the equivalent stress according to Equation (23) computed at the effective plastic strain ε e,pl = 2 3 ∑ i ∑ j ε ij,pl ε ij,pl 1/2 , ε ij,pl is the plastic strain component, s ij is the deviatoric stress component, s * xx = s xx + νs yy , s * yy = s yy + νs xx , and A is given by In Equation (32), H is the plastic modulus computed as H = ∂σ eq ∂ε eq ε eq =ε e,pl using Equation (23).
The consistency condition with a von Mises plasticity assumption states that the plastic strain increment is given by The plane stress condition implies that s zz = σ x + σ y /3 and therefore ∆ε z,pl = 0.

Finite-Element Implementation
The determined parameters n, c 1 , c 2 , and s R and the spatial fields ε f (x, y) and σ f (x, y) are used to specify the material properties in a Finite-Element (FE) model. The randomly generated spatial fields are applied on the integration points and extrapolated to the nodes. This ensures a continuous material property spatial field.

Summary of the Proposed Stochastic Constitutive Model
A flowchart of the proposed method is presented in Figure 8. The parameters 1 , 2 , ρ 12 , and 12 are used to construct the random spatial fields of strength and strain to failure σ f (r, ω) and ε f (r, ω). Given these spatial fields as input to the isotropic constitutive model in Section 3.3, as well as the parameters n, c 1 , c 2 , and s R , the random tangent stiffness matrix D tan (r, ω) = D(r, ω) − D pl (r, ω) can be determined. The randomness of the elastic stiffness matrix D given in Equation (30) at any spatial location r follows from the randomness of the elastic modulus. The latter can be expressed using Equations (25) and (27) according to In Equation (32), is the plastic modulus computed as = , using Equation (23).
The consistency condition with a von Mises plasticity assumption states that the plastic strain increment is given by The plane stress condition implies that = + 3 ⁄ and therefore Δ , ≠ 0.

Finite-Element Implementation
The determined parameters , , , and and the spatial fields ( , ) and ( , ) are used to specify the material properties in a Finite-Element (FE) model. The randomly generated spatial fields are applied on the integration points and extrapolated to the nodes. This ensures a continuous material property spatial field.

Summary of the Proposed Stochastic Constitutive Model
A flowchart of the proposed method is presented in Figure 8. The parameters ℓ , ℓ , , and ℓ are used to construct the random spatial fields of strength and strain to failure ( , ) and ( , ). Given these spatial fields as input to the isotropic constitutive model in Section 3.3, as well As can be seen the parameter s R contributes to the random response through the normal random variable R ∼ N (0, s R ). The same is true for the elastic-plastic tangent stiffness matrix. The parameters 1 , 2 , ρ 12 , 12 as well as n, c 1 , c 2 , and s R are all dependent on the chosen SVE size so that the constitutive response becomes independent of the latter. Poisson's ratio ν = 0.3 is assumed constant and independent of SVE size. The softening behavior at each material point is assumed constant and negligible with ξ = 0.15.

Determination of SVE Size-Dependent Stochastic Constitutive Model Parameters
The eight parameters in the stochastic constitutive model, n, c 1 , c 2 , s R , 1 , 2 , ρ 12 , and 12 are dependent on the choice of SVE size, see Table 5. In this section, the determination of these parameters for an SVE size of 6 × 6 mm 2 and 4 × 4 mm 2 is demonstrated. The parameters n, c 1 , c 2 , and s R in the relation according to Equation (23) and Equation (27) are first determined. A total of 70 SVEs of the same size are cut from random locations in a large fiber network, and a uniaxial test is performed on each SVE. The responses for all SVEs are shown in Figure 9a,b. The response of three SVEs together with the corresponding fitted model according to Equation (23) is shown in Figure 9c,d. Based on the procedure outline in Section 3.3, a value of n = 40.7 was found to minimize the sum of all least square errors from fitting Equation (23) to all 6 × 6 mm 2 numerical tests. The corresponding optimal κ j values, one for each SVE test, resulted in c 1 = 0.5058 and c 2 = 0.7702 after performing the linear regression κ j = c 1 ε f j + c 2 , see Figure 9e. The standard deviation of the residual fitting error was s R = 0.0671 and was used as a deterministic parameter in the stochastic relation according to Equation (27). The random spatial field distributions of ( , ) and ( , ) were characterized by the parameters ℓ , ℓ , , and ℓ , see Table 5, which were determined by the optimization problem according Equation (19). In the following, the determination of the target values for the three constraints as well as initial guesses for the optimization parameters are detailed for the 6 mm × 6 mm SVE size, see Table 6.

Constraints
Initial Guess The random spatial field distributions of ε f (r, ω) and σ f (r, ω) were characterized by the parameters 1 , 2 , ρ 12 , and 12 , see Table 5, which were determined by the optimization problem according Equation (19). In the following, the determination of the target values for the three constraints as well as initial guesses for the optimization parameters are detailed for the 6 mm × 6 mm SVE size, see Table 6.  Following the methodology outlined in Figure 5, a realization of the spatial fields ε f (x, ω) and σ f (x, ω) computed by discretizing a 360 × 6 mm 2 specimen using 6 × 6 mm 2 SVEs is shown in Figure 10a,b. A total of N = 360 points or SVE uniaxial tests are performed with a discretization distance of 1 mm. The target correlation ρ t ε f σ f was determined by computing the correlation between ε f (x, ω) and σ f (x, ω), resulting in ρ t ε f σ f = 0.55. The initial guess for the parameter ρ 12 was given by ρ Figure 10c,d, the probability density functions of ε f and σ f , computed from the 70 SVE tests used in the determination of n, c 1 , c 2 and s R are shown together with the fitted kernel distributions. Following the methodology outlined in Figure 5, a realization of the spatial fields ( , ) and ( , ) computed by discretizing a 360 × 6 mm 2 specimen using 6 × 6 mm 2 SVEs is shown in Figure 10a,b. A total of = 360 points or SVE uniaxial tests are performed with a discretization distance of 1 mm. The target correlation was determined by computing the correlation between ( , ) and ( , ), resulting in = 0.55. The initial guess for the parameter was given by ( ) = = 0.55 . In Figure 10c,d, the probability density functions of and , computed from the 70 SVE tests used in the determination of , , and are shown together with the fitted kernel distributions. The resulting normalized spatial fields ( , ) = ( , ) and ( , ) = ( , ) , computed using the kernel approximation of the PDFs, are shown in Figure 11a,b.

Constraints Initial Guess
In Figure 11c,d, the average distance between two zero level up-crossings was computed for an increasing number of points along the length of the specimen. As can be seen, a specimen length of = 250 mm was necessary for convergence to = 6.3 mm and = 11.9 mm. The initial  The resulting normalized spatial fields , ω)]}, computed using the kernel approximation of the PDFs, are shown in Figure 11a,b. In Figure 11c,d, the average distance between two zero level up-crossings was computed for an increasing number of points along the length of the specimen. As can be seen, a specimen length of L = 250 mm was necessary for convergence to µ t upcr1 = 6.3 mm and µ t upcr2 = 11.9 mm. The initial guesses (0)

Continuum Reconstruction of A 18 × 18 mm 2 Sample and Choice of SVE Size
Validation of the methodology was performed by comparing results from uniaxial tensile tests on the fiber network and the SVE-based continuum model. The boundary conditions applied on the fiber network is shown in Figure 12a. For the proposed continuum model, the same boundary conditions would result in stress concentration and in turn erroneous strain localization initiation. Therefore, the clamped ends in the fiber network were modeled in the continuum model using a contact boundary condition, see Figure 12b. Very low contact penalty stiffness was used to allow free contraction at the boundaries. In order to avoid uneven deformation at the ends due to the low contact stiffness, the boundary nodes were coupled with respect to their axial deformation .
A uniaxial tensile test was performed on a specimen of size 18 × 18 mm 2 using direct simulation and the continuum model based on both a 4 × 4 mm 2 and a 6 × 6 mm 2 SVE size. The constitutive model parameters used are according to Table 5. The discretization distance (see Figure 5a), i.e., the distance between the centers of the SVEs, is half of the SVE size in both the length and width direction. The spatial distribution of stress and strain to failure for both SVE sizes is presented in Figure 13. As can be seen, the weakest material point using the 4 × 4 mm 2 SVE is at the lower right part of the specimen and has a value of 13.90 MPa. The corresponding value using 6 × 6 mm 2 SVEs is 14.86 MPa at approximately the same location. For the smaller SVE size, a higher variation in strength and strain to failure values can be observed within the specimen. (c) (d) Figure 11. Transformed Gaussian random field of (a) ε f and (b) σ f . Average distance between two zero-level up-crossings for (c) g 1 and (d) g 2 as a function of specimen length.

Continuum Reconstruction of A 18 × 18 mm 2 Sample and Choice of SVE Size
Validation of the methodology was performed by comparing results from uniaxial tensile tests on the fiber network and the SVE-based continuum model. The boundary conditions applied on the fiber network is shown in Figure 12a. For the proposed continuum model, the same boundary conditions would result in stress concentration and in turn erroneous strain localization initiation. Therefore, the clamped ends in the fiber network were modeled in the continuum model using a contact boundary condition, see Figure 12b. Very low contact penalty stiffness was used to allow free contraction at the boundaries. In order to avoid uneven deformation at the ends due to the low contact stiffness, the boundary nodes were coupled with respect to their axial deformation u x .
A uniaxial tensile test was performed on a specimen of size 18 × 18 mm 2 using direct simulation and the continuum model based on both a 4 × 4 mm 2 and a 6 × 6 mm 2 SVE size. The constitutive model parameters used are according to Table 5. The discretization distance (see Figure 5a), i.e., the distance between the centers of the SVEs, is half of the SVE size in both the length and width direction. The spatial distribution of stress σ f and strain to failure ε f for both SVE sizes is presented in Figure 13. As can be seen, the weakest material point using the 4 × 4 mm 2 SVE is at the lower right part of the specimen and has a value of 13.90 MPa. The corresponding value using 6 × 6 mm 2 SVEs is 14.86 MPa at approximately the same location. For the smaller SVE size, a higher variation in strength and strain to failure values can be observed within the specimen.  In Figure 14a, the uniaxial response in the x-direction using the continuum model is shown together with the response from the direct simulation. The results for both SVE sizes are in good agreement with direct simulation, with an absolute relative error in strength to failure of 0.4% and -2% for the 4 × 4 mm 2 and 6 × 6 mm 2 SVE size, respectively. The corresponding errors in strain to  In Figure 14a, the uniaxial response in the x-direction using the continuum model is shown together with the response from the direct simulation. The results for both SVE sizes are in good agreement with direct simulation, with an absolute relative error in strength to failure of 0.4% and -2% for the 4 × 4 mm 2 and 6 × 6 mm 2 SVE size, respectively. The corresponding errors in strain to In Figure 14a, the uniaxial response in the x-direction using the continuum model is shown together with the response from the direct simulation. The results for both SVE sizes are in good agreement with direct simulation, with an absolute relative error in strength to failure of 0.4% and −2% for the 4 × 4 mm 2 and 6 × 6 mm 2 SVE size, respectively. The corresponding errors in strain to failure, 6% and −6%, are, however, larger. In Figure 14b, the responses are also compared for an applied load in the transverse y-direction, with an overall performance similar to the x-direction.
In Figure 15, the spatial distribution of uniaxial strain in the loading direction is shown using direct simulation on the fiber network as well as the continuum model using 4 × 4 mm 2 and 6 × 6 mm 2 SVEs. The spatial distributions are shown for an increasing uniaxial load applied in the x-direction. By comparing Figures 13 and 14, it is clear that strain localization was initiated at weak material points. It is also seen that the strain localization pattern from direct simulation matched the ones from the continuum models relatively well. failure, 6% and −6%, are, however, larger. In Figure 14b, the responses are also compared for an applied load in the transverse y-direction, with an overall performance similar to the x-direction.
In Figure 15, the spatial distribution of uniaxial strain in the loading direction is shown using direct simulation on the fiber network as well as the continuum model using 4 × 4 mm 2 and 6 × 6 mm 2 SVEs. The spatial distributions are shown for an increasing uniaxial load applied in the x-direction. By comparing Figure 13 and Figure 14, it is clear that strain localization was initiated at weak material points. It is also seen that the strain localization pattern from direct simulation matched the ones from the continuum models relatively well.

Random Continuum Simulation of 18 × 18 mm 2 and A4 Paper Sample
In this section, random continuum realizations of the fiber network were subjected to a uniaxial loading. The proposed stochastic continuum model presented in Section 3 is used for both stochastic realization (Section 3.2) and constitutive response (Section 3.3). For a specimen size of 18 mm × 18 mm, the results are presented in Figure 16. Two correlated random spatial fields of strength and strain to failure, see Figure 16a,b, are first constructed based on an SVE size of 6 mm × 6 mm, i.e., using ℓ , ℓ , ℓ , and according to Table 5. As can be seen, these randomly generated spatial fields are similar to the spatial fields obtained by cutting and testing direct fiber networks, see Figure  13c,d. The mechanical response is thereafter computed based on these spatial fields as well as the parameters , , , and , see Figure 16c. It should be emphasized that the randomness in the model is not only determined by and but also depends on , as can be seen from the expression of the stiffness according to Equation (34). In Figure 16d, the uniaxial strain defined as the spatial derivative of the displacement field, is computed for the maximum applied stress marked in red in Figure 16c. As can be seen, damage localization was observed consistent with the observation from the analysis of the direct fiber network in Figure 15.

Random Continuum Simulation of 18 × 18 mm 2 and A4 Paper Sample
In this section, random continuum realizations of the fiber network were subjected to a uniaxial loading. The proposed stochastic continuum model presented in Section 3 is used for both stochastic realization (Section 3.2) and constitutive response (Section 3.3). For a specimen size of 18 mm × 18 mm, the results are presented in Figure 16. Two correlated random spatial fields of strength and strain to failure, see Figure 16a,b, are first constructed based on an SVE size of 6 mm × 6 mm, i.e., using 1 , 2 , 12 , and ρ 12 according to Table 5. As can be seen, these randomly generated spatial fields are similar to the spatial fields obtained by cutting and testing direct fiber networks, see Figure 13c,d. The mechanical response is thereafter computed based on these spatial fields as well as the parameters n, c 1 , c 2 , and s R , see Figure 16c. It should be emphasized that the randomness in the model is not only determined by σ f and ε f but also depends on s R , as can be seen from the expression of the stiffness according to Equation (34). In Figure 16d, the uniaxial strain defined as the spatial derivative of the displacement field, is computed for the maximum applied stress marked in red in Figure 16c. As can be seen, damage localization was observed consistent with the observation from the analysis of the direct fiber network in Figure 15. The major potential of the stochastic continuum model is demonstrated by the construction of larger specimens for which results cannot be validated using direct numerical simulation. A 297 × 210 mm 2 specimen (A4 paper size) was constructed and analyzed using the proposed model based on a 6 mm × 6 mm SVE size. Similarly to the result presented in Figure 16, the constructed spatial fields of strength and strain to failure are presented in Figure 17a,b, respectively, and the global uniaxial response of the specimen is presented in Figure 17c. The spatial distribution of uniaxial strain computed at the maximum global uniaxial stress, marked by a red point in Figure 17c, is presented in Figure 17d. By comparing Figure 17c and Figure 16c, it is seen that the response is more brittle for the larger specimen. This is, for instance, seen by comparing the softening branch of the response in both figures. The major potential of the stochastic continuum model is demonstrated by the construction of larger specimens for which results cannot be validated using direct numerical simulation. A 297 × 210 mm 2 specimen (A4 paper size) was constructed and analyzed using the proposed model based on a 6 mm × 6 mm SVE size. Similarly to the result presented in Figure 16, the constructed spatial fields of strength and strain to failure are presented in Figure 17a,b, respectively, and the global uniaxial response of the specimen is presented in Figure 17c. The spatial distribution of uniaxial strain computed at the maximum global uniaxial stress, marked by a red point in Figure 17c, is presented in Figure 17d. By comparing Figures 17c and 16c, it is seen that the response is more brittle for the larger specimen. This is, for instance, seen by comparing the softening branch of the response in both figures.

Simulation of Stochastic Failure in Paper-Making Machines
Consider a paper-making machine, Figure 18a  Using the proposed method, random realizations of 1 m × 8 m samples can be simulated. In Figure 19, three constructed spatial fields of strength and strain to failure are shown. It is seen that, for all three realizations, the difference between the strongest and weakest material points within the specimen is high due to the large considered size. From the uniaxial response also shown in the figure, a small variation in strength can be observed. The variation in corresponding strain to failure is, however, larger. This statistical variation in material response demonstrated using the proposed model, together with the variation in the applied load, can be seen as a major reason for the random occurrence of breaks in paper machines. From this example, it can be concluded that the disordered nature of the fiber network contributes to the statistical variation in mechanical response, even for such large specimens.

Simulation of Stochastic Failure in Paper-Making Machines
Consider a paper-making machine, Figure 18a, with a construction that results in breaks that mostly occur in an open draw of 1 m [18], see Figure 18b.

Simulation of Stochastic Failure in Paper-Making Machines
Consider a paper-making machine, Figure 18a  Using the proposed method, random realizations of 1 m × 8 m samples can be simulated. In Figure 19, three constructed spatial fields of strength and strain to failure are shown. It is seen that, for all three realizations, the difference between the strongest and weakest material points within the specimen is high due to the large considered size. From the uniaxial response also shown in the figure, a small variation in strength can be observed. The variation in corresponding strain to failure is, however, larger. This statistical variation in material response demonstrated using the proposed model, together with the variation in the applied load, can be seen as a major reason for the random occurrence of breaks in paper machines. From this example, it can be concluded that the disordered nature of the fiber network contributes to the statistical variation in mechanical response, even for such large specimens. Using the proposed method, random realizations of 1 m × 8 m samples can be simulated. In Figure 19, three constructed spatial fields of strength and strain to failure are shown. It is seen that, for all three realizations, the difference between the strongest and weakest material points within the specimen is high due to the large considered size. From the uniaxial response also shown in the figure, a small variation in strength can be observed. The variation in corresponding strain to failure is, however, larger. This statistical variation in material response demonstrated using the proposed model, together with the variation in the applied load, can be seen as a major reason for the random occurrence of breaks in paper machines. From this example, it can be concluded that the disordered nature of the fiber network contributes to the statistical variation in mechanical response, even for such large specimens.

Discussion
In this paper, a stochastic constitutive model for disordered fiber networks, derived through the mechanical response of stochastic volume element, is presented. The model is based on two correlated random spatial fields representing the spatial distribution of strength and strain to failure of the SVEs. These spatial fields are constructed using a multivariate kernel characterized by four parameters that are determined by analyzing realizations of direct fiber networks. The multiaxial non-linear mechanical response is modeled using an isotropic plasticity model with four additional parameters. The randomness in the material response is determined through the stochastic spatial fields ( , ) and ( , ) as well as the parameter in the constitutive model. All model parameters are dependent on the chosen SVE size. This assures that the global response is relatively independent of SVE size. The independence of SVE size is achieved as long as the ratio between the spatial correlation length and mesh size is large [22]. However, a large enough SVE should be chosen in order to limit the influence of boundary conditions during the numerical tests performed on SVEs cut from the direct fiber network. It should be noted that, the constitutive model obtained using SVEs face two sources of uncertainties: one contribution resulting from the applied boundary

Discussion
In this paper, a stochastic constitutive model for disordered fiber networks, derived through the mechanical response of stochastic volume element, is presented. The model is based on two correlated random spatial fields representing the spatial distribution of strength and strain to failure of the SVEs. These spatial fields are constructed using a multivariate kernel characterized by four parameters that are determined by analyzing realizations of direct fiber networks. The multiaxial non-linear mechanical response is modeled using an isotropic plasticity model with four additional parameters. The randomness in the material response is determined through the stochastic spatial fields ε f (x, y) and σ f (x, y) as well as the parameter s R in the constitutive model. All model parameters are dependent on the chosen SVE size. This assures that the global response is relatively independent of SVE size. The independence of SVE size is achieved as long as the ratio between the spatial correlation length and mesh size is large [22]. However, a large enough SVE should be chosen in order to limit the influence of boundary conditions during the numerical tests performed on SVEs cut from the direct fiber network. It should be noted that, the constitutive model obtained using SVEs face two sources Materials 2019, 12, 538 26 of 28 of uncertainties: one contribution resulting from the applied boundary conditions and the other one from the uncertainties in the micro-structure. However as is seen in Reference [22], the uncertainties resulting from the micro-structure randomness, is assumed more important than the ones resulting from the applied boundary conditions.
The proposed methodology based on SVEs is validated by comparing the mechanical response of 18 mm × 18 mm samples using direct numerical simulation on the fiber network and the proposed continuum model. The results are compared through uniaxial loading in both the longitudinal and transverse direction. It has been shown that both mechanical response and strain localization pattern using the continuum model matches the one from the direct simulation relatively well. Based on the demonstrated accuracy, it may be concluded that using two spatial fields ε f (x, y) and σ f (x, y) is enough to characterize the material behavior and that the error from neglecting other spatial field variations is small. An important aspect explaining the error observed in the global strain and strength to failure is the assumption of isotropy of the SVE in the constitutive model. In future work, the use of an anisotropic plasticity model to describe the mechanical response of the SVE may further increase the accuracy in the global response. It should furthermore be noted that the results are mesh independent, since failure is assumed to occur when the maximum global stress is reached, and the mesh-dependent global softening behavior is therefore not of interest.
The proposed model paves the way to quantifying uncertainties in the response of thin fiber networks using the stochastic constitutive model. The eight parameters in the model can be linked to microstructural features [18] of the fiber network which opens the possibility to designing the network such that the variability in the global response is limited. For instance, it has been noted that the spatial fields of strain to failure and strength have different characteristic length. This is directly linked to microstructural features such as fiber connectivity and fiber orientation, which may affect the spatial fields differently. In the future, it is intended to enrich the present description by performing reliability estimation [41,42] and reliability-based design optimization (RBDO) [43], where both material and load variability are incorporated.

Conclusions
The conclusions in the papers can be summarized as follows.
• A stochastic constitutive model was proposed for isotropic fiber networks.

•
The effect of fiber disorder for isotropic fiber networks can be quantified at the macromechanical level using the proposed method.

•
In order to implement the stochastic constitutive model, eight SVE size-specific material parameters are necessary.

•
The randomness in the model was simulated using random field realizations of strain to failure and strength. An epistemic uncertainty parameter also contributes to the randomness. • Different characteristic lengths for the spatial fields of strain to failure and strength are assumed, which is in agreement with what is observed in direct fiber network simulation. • The proposed model speeds up the simulation time of fiber networks compared to direct numerical simulations. A 24 mm × 24 mm simulation takes 2 min on a modern 128 GB RAM supercomputer while a corresponding direct fiber network simulation takes 2 days. • The constitutive model captures the transition from ductile behavior for small sample size to semi-brittle behavior for larger sample size (Section 4.3).

Conflicts of Interest:
There is no conflict of interest.