The influence of structural variations on the constitutive response and strain variations in thin fibrous materials

The stochastic variations in the structural properties of thin fiber networks govern to a great extent their mechanical performance. To assess the influence of local structural variability on the local strain and mechanical response of the network, we propose a multiscale approach combining modeling, numerical simulation and experimental measurements. Based on micro-mechanical fiber network simulations, a continuum model describing the response at the mesoscale level is first developed. Experimentally measured spatial fields of thickness, density, fiber orientation and anisotropy are thereafter used as input to a macroscale finite-element model. The latter is used to simulate the impact of spatial variability of each of the studied structural properties. In addition, this work brings novelty by including the influence of the drying condition during the production process on the fiber properties. The proposed approach is experimentally validated by comparison to measured strain fields and uniaxial responses. The results suggest that the spatial variability in density presents the highest impact on the local strain field followed by thickness and fiber orientation. Meanwhile, for the mechanical response, the fiber orientation angle with respect to the drying restraints is the key influencer and its contribution to the anisotropy of the mechanical properties is greater than the contribution of the fiber anisotropy developed during the fiber sheet-making.


Introduction
Bio-based fibrous materials are broadly used in a number of industrial applications [1] of which paper and board-making are the largest by volume and total revenue. The usage of these materials has gradually advanced from conventional to more sophisticated applications where they are used as substrates for sensors, actuators and circuit boards [2] . The biodegradability, renewability, efficient production technique and low prices provide the preference of using those materials over synthetic materials like glass and plastic [3] . Furthermore, a key advantage of these materials is the ability to control the anisotropic properties and enhance the performance of the material in the loading direction [4] .
Paper-based materials, for which fibers are the main constitutive components, are characterized by randomness in their mechanical response. This randomness is originated in the production process, in which fibers, fiber segments, and other chemical that are correlated to the locus of failure initiation [18] . Moreover, in an experimental investigation, the influence of local structural properties on the local strain field has been quantified, showing that the thickness, basis weight, density, and local fiber orientation combined can justify 39% of the total variation in the local strain [19] .
Using full 3D fiber-network generation and simulations, the strain localization pattern was found to vary randomly as does the global mechanical response [20] . The complex inhomogeneous nature of paper materials requires the study of different scales in order to model the material behavior [21] . Three discrete scales can be recognized: microscale (fiber scale), mesoscale (network scale), and macroscale (sheet scale) [21] . The fiber scale tackles the mechanical response of a single fiber and its geometrical properties (see Section 3.1 ). The mesoscale is the size of the fiber network for which a continuum-based model with effective properties is proposed. In this study, it is defined by a network size of 4 × 4 mm 2 (see Section 3.2 ). The three different scales for the studied paperbased material is shown in Fig. 1 .
To bridge the gap between the scales of heterogeneous materials, multiscale and homogenization approaches are widely used for composites [22][23][24][25] and were recently reviewed in [26] . The heterogeneous fiber networks received less attention with some examples with collagen networks [27] and paper [28] . In [29] , the authors employed a stochastic multiscale approach consisting of microscale, mesoscale, and macroscale to model an isotropic fibrous sheet to capture the strain localization. The effect of drying was not accounted for. The data for the spatial strain and strength to failure fields were taken from the artificially created sheets rather than from experiments. Those fields have been shown to be the main predictors of the random strain localization pattern.
The aforementioned studies lack the ability to identify the individual impact of the structural properties on the local strain field as well as on the global mechanical response. In experimental studies, singling out the influence of one structural property by imposing the other properties to be constant is elusive. Also, using micromechanical models to investigate the heterogeneity in the material structure is computationally formidable at the product scale. Overcoming these limitations is crucial for proper reliability analysis [ 30 , 31 ] and product scale design [32] of fibrous materials considering local random variability.
In this work, the impact of spatial structural variability is assessed using a combined experimental, numerical and continuum model approach. The variability is first quantified through experimental local measurements of the thickness, density, fiber orientation and fiber anisotropy. Micromechanical fiber network simulations, with fiber properties derived from experimental observations, are thereafter employed to study the influence of the structural properties on the mesoscale mechanical response. A novelty in the micromechanical simulations is the ability to account for the impact of the restrained drying condition during the production process on the fiber properties. A continuum model of the mesoscale mechanical response is thereafter proposed with material parameters that are function of the local properties. The impact of variability is finally assessed by applying the measured spatial fields of thickness, density, fiber orientation and fiber anisotropy in a macromechanical Finite-Element (FE)-model where each element is governed by the proposed mesoscale continuum model. The accuracy of the proposed approach with regard to the prediction of the mechanical response and spatial strain field under uniaxial loading is evaluated by comparison to experimental measurements.

Measurement of local strain
In this experimental study, four sack paper specimens were exposed to a uniaxial tensile load in the Machine Direction (MD) with a strain rate of 40 mm/min under a standard climate of 23 °°C and 50% relative humidity. Their local strain variation is recorded just prior to uniaxial tensile rupture which is estimated at 2.2% strain. The full size of the specimens is 40 × 65 mm 2 with a Region of Interest (RoI) of 32 × 56 mm 2 defined inside the deformed samples as seen in Fig. 2 . The local strain was measured with Digital Image Correlation (DIC) [33] with a resolution of 0.74 mm/pixel. The local strain map within RoI was resized and aligned with the measured structural property maps described in Section 2.2 with a resolution of 1 mm/pixel by utilizing a landmark based registration and a shape preserving coordinate transform [34] .

Measurement of thickness and density
After determining the local strain, the local thickness was measured with a twin laser profilometer [35] with a resolution of 100 μm/pixel within RoI. Local basis weight was determined with a β-radiographic transmission method [36] with a resolution of 50 μm/pixel within RoI. The measured thickness and basis weight maps were then aligned with each other with a resolution of 1 mm/pixel. Finally, the density map was calculated by dividing the aligned basis weight map with the aligned thickness map.

Measurement of fiber orientation and anisotropy
The local fiber orientation was measured with a sheet splitting method [37] . Each specimen was split in the thickness direction into 14-20 layers and subsequently scanning the layers by an optical scanner with a resolution of 13 μm/pixel. For each layer, the Fiber Orientation Distribution (FOD) is defined for subsets of 1 × 1 mm 2 . The fiber orientation map for all the layers in the local 1 × 1  mm ² subsets was then determined by averaging the FOD across the layers in the thickness direction and subsequently fitting ellipses to the averaged distributions.
The local Fiber Orientation Distribution (FOD) is defined by two parameters, orientation and anisotropy, derived from the direction and shape of the fitted ellipse as described in Section 3.2 , see where a and b are the major and minor axes of the local fiber orientation distribution ellipse, respectively.
The local variations in the measured structural properties as well as the strain fields are shown in Fig. 4 for two of the studied specimens. In the remainder of this paper, a methodology for quantifying the individual impact of local variability of each of the studied structural properties on both the strain field and mechanical response is presented.

Fiber network simulation
The sampling of the material data on the mesoscale from the experiments has fundamental difficulties such as a laborious way of analyzing the data by splitting, the need for testing in different directions beyond the yield stress, and performing well-controlled tests on relatively small samples. In order to avoid these limitations, we used a micromechanical simulation framework. In this framework, we reconstruct the fiber network with a controlled density, thickness, and fiber anisotropy. This enables performing the mesoscale (network scale) mechanical analyses in various directions, varying degrees of anisotropy, and yet maintaining other structural properties constant.

Microscale modeling
In order to study the influence of structural variability, we have employed micromechanical fiber network simulations. The geometries of the fiber including fiber length, width, shape factor, wall thickness and width to high ratio, are acquired using a characterization methodology outlined in [4] . These geometrical data are extracted based on measurements of wet kraft pulp using Fiber Morphology Analysis (FMA) [38] . The cross-sectional data are corrected by processing the scanned images of the kraft sheet with the Micro-Computed tomography ( μCT) [39] . The pulp used to produce the sheets was an unbleached (Kappa value is 43) softwood Kraft chemical pulp fibers being a mixture of spruce (~80%) and pine (~20%) wood and was characterized with the methods described in [4] . The relevant pulp data are listed in Table 1 .
The network generation is performed by randomly depositing the fibers [40] in a planar surface where the curvature of the fiber is constant and parallel to the deposition plane. During the deposition, the intersection between the fibers is determined and the intersection point is lifted up to avoid penetration, see Fig. 5 (a). The fiber geometry is thereafter smoothed in order to avoid kinking. This random deposition of fibers is continued until the specified grammage is reached. The density of the softwood pulp reported in the literature ranges between 1200 kg/m 3 [41] to 1500 kg/m 3 [42] and is dependent on the lignin content. The adopted value used during the deposition in this study is 1400 kg/m 3 . The code for generating the network along with the documentation is available as supplementary material in [29] .
The fibers are represented using a finite-element model by a series of 3-node 3D Timoshenko beam elements. The cross-section is circular with a solid or hollow cross-section. The pointwise contact beam to beam with traction and separation law is used to model the bond between fibers [40] . In this work, we did not consider debonding between the fibers as this phenomenon does not significantly influence the stress-strain curves in the considered strain interval before the ultimate failure and softening [43] . The softening behavior was not included in the constitutive model either. Not including them can, however, prohibit capturing accurately the strain field locally as the strain localization and the associated bond failures can occur prior to the global failure. The constitutive response of the fiber is modeled using bilinear plasticity with an isotropic hardening law. The flow stress is given by σ s, f + E tan,f ε e,pl where σ s,f is the initial yield stress of the fiber, E tan,f is the tangent modulus and ε e , pl is the equivalent plastic strain. The latter is incrementally computed from the plastic stain increments dε pl as dε e , pl = ( 2 3 d ε T pl d ε pl ) 1 / 2 which has three components (one normal and two transverse shear strains) in the used Timoshenko beam element. The material properties of the fibers are determined in Section 3.3 .  Table 1 The geometrical properties of fibers used in the micromechanical simulation tool.

Mesoscale simulation
Direct simulations of 4 × 4 mm 2 fiber networks with grammage 70 g/m 2 are used to quantify the influence of fiber orientation and degree of anisotropy on the constitutive response. The chosen 4 × 4 mm 2 size is small enough to be representative of the mesoscale, i.e. to represent a discretization size that captures local variation in the macroscale. The relatively large size of the mesoscale ensures that 1) random realizations having the same thickness, density, fiber orientation and degree of anisotropy result in approximately the same constitutive response and 2) the influence of the boundary condition on the uniaxial response of the mesoscale fiber network is small.
The effect of fiber orientation θ on the constitutive response is studied using three uniaxial tests in θ = 0 , θ = 45 • and θ = 90 • with respect to the loading direction as seen in Fig. 5 (b). These The response of the uniaxial tests in Fig. 5 (b) for two different degrees of anisotropy, is shown in Fig. 6 (a). As can be seen, the difference in response between the x and y directions expectedly decreases with decreased degree of anisotropy, while the 45 °d irection remains almost unchanged. For the isotropic case, see Fig. 6 (b), a similar response is observed in all the three directions. These results do not take into account the effect of restrained drying in MD, which is modeled in Section 3.3 .

Effect of restrained drying condition during the production process
During the papermaking process, the fiber web, which is mainly constituted of fibers, is drained and the fiber network consolidation takes place. In the sack paper investigated in this study, the shrinkage is restrained in MD while free shrinkage is allowed in Cross Direction (CD) during the drying. This drying process is well known to influence the properties of the produced materials [ 45 , 46 ], by increasing the stiffness [47] in the restrained drying direction and reducing it in the free shrinkage direction. The physics behind this change is attributed to the impact on the microstructure of the fibers with micro-compressions and realignment of the microfibrils inside the fiber [ 47 , 48 ]. Consequently, in a sheet subject to drying, the fiber properties will change according to their orientation with respect to the drying direction [49] .
In the direct simulation tool, the effect of restrained drying in MD and free shrinkage in CD is implemented by modeling the fiber property as a function of its orientation angle θ f with respect to MD according to where E f is the fiber elastic modulus, E f MD and E f CD are fitting parameters representing the elastic modulus of a fiber oriented in MD and CD, respectively, see Fig. 7 . The selection of this function, although being empirical, is motivated by the analytical transformations presented in ref. [50] . The same relation is assumed to  Table 2 . These are found by matching the fiber network responses to the experimental uniaxial responses in MD, CD and 45 °directions. It is noted that Eq. (2) and the corresponding fitting parameters from Table 2 are only valid for the case of restraint drying in the MD and freely dried in the CD. This is typical for the commercially produced paper web with the exception of the web edges, where the effect of constraints may affect the shrinkage in the CD [51] .
From the experimentally studied specimens, it was found that the average anisotropy is λ= 0.34. Therefore, a fiber network with the same degree of anisotropy ( λ= 0.34) was generated using the micromechanical tool. The mechanical response of this fiber network was fitted to the stress-strain curves for the same material of the sack paper. Once the MD response was fitted, the CD and 45 °r esponses were computed. Fig. 8 (a) shows that with the measured anisotropy alone we were unable to achieve a reasonable agreement with the experiment conducted in all the directions. However, by adding the influence of drying a good fit was achieved, see Fig. 8 (b). This is explained as follows. As Fig. 7 shows, the selected function for the angular dependency of the fiber properties is not symmetric with respect to 45 °angle and is flatter toward the CD. This means the properties of the fibers change faster toward those in the CD in order to capture the experimental stress-strain curves in all directions as in Fig. 8 (b). Neglecting this angular dependency with respect to the drying direction and assuming the same fiber properties regardless of orientation results in the mismatch against the experiment shown in Fig. 8 (a). This result shows the impor-  tance of including the effect of drying as well as the impact it has compared to the anisotropy alone. It is concluded that the effect of the drying has a greater contribution to the anisotropy in the mechanical properties of the sheet compared to the effect of the fiber anisotropy alone.
In a large paper sample, the anisotropy may vary locally but this variation averages out on the continuum level and does not influence the global stress-strain response of the network as will be demonstrated later. Also, the network may show some small variability owing to its disordered nature even with the same structural parameters. This phenomenon may be size-dependent with smaller samples showing greater variability. To quantify it, we tested 10 networks with the same given set of anisotropy, thickness, and density but different random seeds used during the network generation. We also compared it with the variability observed in the experiment on a larger sample. The average stressstrain curves and the error bars representing one standard deviation are shown in Fig. 8 together with the experimental curves measured on a large sample. These results demonstrate that although the variation exists it is not significant and comparable to the one observed in the experiment. Therefore, the selected fitting parameters are still adequate. It is also important to note, that the greatest variability observed experimentally was in recorded strength values, which are not addressed in this study.

Mesoscale continuum model
In this Section, we will present a continuum model for the mesoscale response taking into account the fiber anisotropy and density. The fiber orientation and thickness will be accounted for directly in the macro-mechanical finite-element analysis in Section 5 .

Plane stress state with hill's yield criterion and isotropic power law hardening
The 2D plane stress model can be written in rate form as ˙ σ = D tan ˙ ε where σ = [ σ x , σ y , τ xy ] T is the stress vector, ε = [ ε x , ε y , γ xy ] T is the strain vector and D tan = D el -D pl is the local tangent stiffness matrix computed from the elastic and elastic-plastic stiffness matrix D el and D pl , respectively. These are defined in a local coordinate system where the x -axis points in the fiber orientation direction following the definition in Fig. 3 and Fig. 5 (c). The elastic stiffness matrix with a plane stress condition can be written as where the elastic modulus in x and y direction, E x and E y , respectively, as well as the Poisson's ratio ν xy and shear modulus G xy are functions of the local properties. The Poisson's ratio is given by the empirical expression [52] ν xy = 0 . 293 E x / E y , where the constant 0.293 is the Poisson's ratio for the isotropic case. The shear modulus is computed from E x , E y and E 45 , where the latter is the elastic modulus in the 45 °direction with respect to the x -axis, according to [53] The plastic behavior of the thin fiber network is assumed to follow Hill's yield criterion [54] , which is expressed based on the ratio R ij of the yield stress in the direction ij with respect to a reference direction. Choosing the x -axis as the reference direction, i.e. R xx = 1, and setting R zz = 1, results in the 2D Hills criterion where σ f is the flow stress and H is the Hills orthotropic coefficient matrix given by The flow stress σ f is given by a power-law hardening where c and d are power-hardening constants, σ s is the initial yield stress and ε e , pl = where the derivatives of the yield surface and flow stress are computed from Eq. (6) and Eq. (8) as ∂ f

Material properties as a function of local anisotropy
The local tangent stiffness matrix derived in Section 4.1 depends on the local anisotropy λ. This dependency is modeled based on the uniaxial mesoscale fiber network responses in the local x , y and 45 °directions. In Fig. 9 (a) The shear modulus computed using Eq. (5) is also shown in Fig. 9 (b). In Eq. (10) , k 1 , k 2 and k 3 are fitting constants such as k 1 is the parameter value for λ = 0, k 2 is a dimensionless scale factor and k 3 is a dimensionless growth factor. The fitted values are presented in Table 3 for each material parameter. The sign of k 2 determines if the material property is increasing ( k 2 > 0) or decreasing ( k 2 < 0) with increased degree of anisotropy λ. The physical meaning of k 2 < 0 is that, as the degree of anisotropy increases, less fiber will be in oriented in the CD resulting in a decrease in E y , R yy and R xy . From Table 3 it is noted that the hardening exponent d is assumed to be constant, since k 3 = 0. Although d is not completely constant, its change with respect to λ is small (below 4%) and therefore neglected. It allowed reducing the number of varied parameters while retaining the good quality of the fit across all values of λ. It is also noted that, during the fitting procedure, the parameters are treated as independent. However, for physical reasons, a certain correlation between the parameters exists. For instance, as the degree of anisotropy increases, more fibers are oriented in MD and fewer fibers in CD. This results in both an increase of E x and a decrease of E y . Due to the restrained drying condition in MD, the increase in E x is not equal to the decrease in E y and the parameter fitting is therefore performed for each individual parameter independently.

Material properties as a function of density
The effect of density ρ is modeled by linear scaling [55] of the elastic moduli E x , E y and E 45 , the yield stress σ s and the power Table 3 Fitting constants k 1 , k 2 and k 3 in Eq. (10) , determining the influence of anisotropy on the mesoscale material properties.  law hardening constant c . The corresponding fitting parameters k 1 in Table 3 is multiplied by the density coefficient where ρ is the average density of the sheet and ρ per ≈ 200 kg/m 3 [55] is the percolation point of the density. A density lower than ρ per corresponds to insufficient number of fibers for establishing the connectivity across the network.

Model validation and impact of variability
In this section, the spatial strain distribution and the stressstrain response predicted by macro-mechanical FE-simulations are first validated against the experimental results. Thereafter, the impact of local structural variability is evaluated.

Comparison of simulation and experiment
The local variations in the measured structural properties are averaged over the 4 × 4 mm 2 mesoscale size and inputted in the FE-analysis as can be seen in Fig. 10 for specimen A and B, respectively. It shows that the local averaging introduces discontinuities across the domain. Each 4 × 4 mm 2 is further meshed with the mesh size of 1 × 1 mm 2 . These spatial fields of thickness, density, fiber orientation and fiber anisotropy are used as input in a FE shell model with boundary conditions according to Fig. 11 . Each finite-element is governed by the mesoscale continuum model presented in Section 4 with elemental thickness input and local coordinate system pointing in the local fiber orientation direction. The orientation with respect to the drying direction is set to be constant throughout the sheet. Fig. 12 (a) shows the comparison between the 4 stress-strain curves from the characterized samples and the simulated results. Although the curves are relatively close, the simulated results are lower and there is less variability in them. The reason for this is that the examined samples did not have the average Fiber Orientation (FO) equal to zero, as was assumed during the fitting. Since the samples were extracted from arbitrary positions across the width, the average FO was affected by the fact that the forming section of the paper machine may not ensure that the FO is strictly in the MD. It deviates from zero toward the edges of the web. By observing the results, one can see that paper with FO closer to the MD is stiffer and those with the greatest FO are more compliant. The fitting of the micromechanical model was done against an independent set of measurements on the samples with unknown fiber orientation and the assumed FO was zero. With this data fitted, and the coordinate system rotated during computation, this yielded a softer average response compared to the original fitting. To improve the match, one should have accounted for the individual fiber orientation for each sample during the fitting while keeping the drying direction along the MD. In this case, however, one should have the curves in all three directions (MD, CD, and 45 °) for a given FO, and getting this data for the same specimen is impossible as the testing in each direction surpasses the yield limit. Nevertheless, this result shows the importance of knowing the fiber orientation for meaningful comparison and its impact on the stressstrain curves, which will be further explored in Section 5.2 . Fig. 12 (b-e) shows the comparison between the computed and measured strain fields. In addition to the visual comparison of the strain fields, relative error plots quantifying the difference between the measured and simulated strain fields are shown in Fig. 12 (fi). They show an overall good agreement in samples B and D. In samples A and C, the agreement degrades locally where the strain localization takes place, while it remains good outside the localizations. The inability to account for strain softening that may take place locally is a plausible reason for not capturing the localization accurately. The point-wise Pearson correlation coefficient r [56] of measured and simulated strain fields are computed to r = 0.62 for Specimen A, r = 0.60 for Specimen B, r = 0.40 for Specimen C and r = 0.80 for Specimen D. Therefore, the multiscale FE-model is capable to capture the regions with higher strains and the regions which are almost intact relatively well.

Influence of local variations
The spatial variability of each of the four studied structural properties is quantified by their Coefficient of Variation (COV) de-  fined as the standard deviation of the spatially varying property divided by its mean value. The influence of this variability on the strain field is studied as follows. The spatial variations of thickness, density, fiber orientation and anisotropy are excluded one by one by setting these properties equal to their mean values t , ρ, θ and λ , respectively. The corresponding simulated strain fields resulting from uniaxial loading at 2.2% strain are shown in Fig. 13 (a) for Specimen A and Fig. 13 (b) for Specimen B. The Reference strain fields in Fig. 13 are the simulated strain fields in Fig. 12 (b) and (c), which are added for a simpler visual comparison. To quantitatively assess the impact of spatial variability on the strain field, the point-wise Pearson correlation coefficients r between the Reference strain field and each of the other strain fields are computed. The impact of variability, defined as 1 − r, is plotted in Fig. 14 (a) as a function of the experimentally measured COV of each of the four structural properties for the four studied specimens. A linear regression with corresponding 95% confidence interval reveals that the spatial variability in density presents the highest impact on the strain field followed by the thickness and fiber orienta-tion, while the impact of spatial variability in fiber anisotropy is small. The linear regression starts from the origin which represents the theoretical case where the structural property is spatially constant (COV = 0) and consequently implying that there is no influence of spatial variability. It is noted that, even though the thickness and density both linearly scale the elemental properties in the FE-implementation, the higher influence of the density is predominantly attributed to the percolation point introduced in Eq. (11) . It is also noted that the confidence interval is relatively large due to the stochastic nature of the studied material and the limited number of specimens due to the laborious nature of the experimental study.
The corresponding impact on the uniaxial response is presented in Fig. 14 (b) for Specimen B, where the spatial variations of thickness, density, fiber orientation and anisotropy are excluded one by one by setting these properties equal to their mean values. Similar stress-strain responses are obtained for the other studied specimens. It is noted that the spatial strain fields at 2.2% strain are given in Fig. 13 (b). As can be seen, only a small influence of spatial variability in fiber orientation is observed and the impact of variability of the other structural properties is minimal. The results in Fig. 14 also show that fiber orientation variation is somewhat less important to the strain field than to the stress-strain-behaviors. The importance of fiber orientation has already been emphasized in the discussion of the results in Fig. 12 (a).

Discussion
In this work, we developed a combined numerical, modeling and experimental approach to assess the impact of local variability of thin fiber networks. The developed process consists of four major steps. fields are used as input in a macromechanical FE-simulation where each element is governed by the Multiscale continuum model . The impact of spatial variability is assessed from the FE-simulations by setting each spatial field to be constant and equal to its average value.
The advantage of the above approach is the ability to decouple the individual impact of local structural variabilities on the mechanical behavior of the material. This isolation of impact is difficult to achieve using a pure experimental method. In addition, the proposed approach bridges the gap between detailed computationally expensive micro-mechanical simulation [20] and the continuum approach [57] which omits the heterogeneity of the materials.
Experimental measurements of the strain field upon loading is performed in order to validate the prediction accuracy of the proposed model. The comparison between the simulation results and the experiments shows an excellent agreement for the mechanical response and a good agreement for the strain field. The deviation between the simulated and measured strain field is explained by 1) the experimental measurement uncertainties, 2) the difference between the continuum finite-element and the fiber network in terms of plasticity and bond behavior and 3) the assumption that the effect of drying on the fiber properties, modeled by Eq. (2) , is independent of the degree of anisotropy. Although the latter point is admittedly one of the limitations of the present approach, there are no reports available to put this assumption to the test.

Conclusions
A multiscale methodology is proposed to quantify the influence of spatial variability of structural properties due to the disordered nature of fiber networks. The proposed method combines detailed micro-mechanical simulations, physical measurement of fiber-level variability and continuum modeling at the mesoscale level. The experimental validation shows an excellent agreement of the macro-mechanical uniaxial response and a good agreement of the predicted strain field.
By using the proposed methodology, we could, for the first time, separate the effect of the fiber and drying anisotropy. We found that the drying anisotropy plays a major role in contributing to the anisotropy of the mechanical properties.
We used the proposed method to identify the role of the local structural variations on the strain variability and mechanical response in the uniaxial tensile test. Among the tested spatial fields, we found that the variations of density followed by thickness and fiber orientation present the largest impact on the strain field, while the fiber anisotropy does not. For the mechanical response, fiber orientation is the principal influencer while other properties are of less importance.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.