The Elasto-plastic Behaviour of Three-dimensional Stochastic Fibre Networks with Cross-linkers

Fibre network materials constitute a class of highly porous materials with low density, promising for functional and structural applications; however, very limited research has been conducted, especially on simulation and analytical models. In this paper, a continuum mechanics-based three-dimensional periodic beam-network model has been constructed to describe the stochastic fibre network materials. In this model, the density of the cross-linkers is directly related to the relative density of the fibre network materials, and the cross-linkers are represented by equivalent beam elements. The objective of this work was to delineate the elasto-plastic behaviour of the stochastic fibre network materials. Characteristic stress and strain derived from the total strain energy density have been adopted to reveal the yielding behaviour of the fibre networks. The results indicate that the stochastic fibre network materials are transversely isotropic. The in-plane stiffness and strength are much larger than those in the out-of-plane direction. For the fibre network materials with a small relative density, the relationship between the uniaxial yield strength and the relative density is a quadratic function in the x direction and is a cubic function in the z direction, which agree well with our dimensional analysis and are consistent with the relevant experimental results in literature. The yield surface depends strongly on the relative density and the connection between fibres.

3 produced by overlapping the randomly distributed fibres by air-laid web-forming technology before compressing and sintering into three-dimensional materials with a fibre network (Xi et al., 2011). Experiments (Zhao et al, 2013) show that the MFSSs are transversely isotropic and that the in-plane stiffness and strength are much higher than those in the out-of-plane direction, which is different from foams or honeycombs. Connectivity in the fibre network materials plays a significant role in determining the mechanical properties. It was found that the in-plane deformation of porous materials with high nodal connectivity is stretching dominated. In contrast, the deformation mode of low density two-dimensional hexagonal honeycomb is a combination of cell wall stretching and bending when the nodal connectivity is low (Symons and Fleck, 2008). Jin et al (2013) developed a two-dimensional random micromechanics beam model to investigate the in-plane elasto-plastic behaviour of MFSSs ). In their model, all the fibres are completely overlapped with each other to form a two-dimensional stochastic fibre network, which leads to very strong bonding connections and high nodal connectivity between fibres. Xi et al (2011) proposed that the mechanical properties of metal fibre porous materials are highly dependent upon the fibre-fibre joints and the number of metallurgy nodes. The mechanical properties are enhanced with increasing sintering contact points per unit volume and the bonding intensity. Some analytical models have been developed based on the assumption of affine deformation of the fibre network Markaki and Clyne, 2005;Picu, 2011;Tsarouchas and Markaki, 2011).
The relative density can significantly affect the strength and stiffness of porous materials (Gibson and Ashby, 1997;Chen et al, 1999;Jin et al, 2013;Won et al, 2013;Zhu et al, 2001Zhu et al, , 2000Zhu et al, , 1997 ; however, it is difficult to control the relative density in the manufacturing process. Finite element method (FEM) which was originally developed for solving solid mechanics problems offers a means to probe the mechanical properties of intricate stochastic fibrous materials by controlling the relative density and other key parameters in the model. In 4 addition, optimized design of complex porous materials can be realized by finite element method. In Finite Element (FE) modelling, beams are mostly utilized to represent fibres and the way to treat the connections between the beams/fibres is crucial. A comprehensive study on the modelling of stochastic fibrous materials by mathematical treatment, for instance, the probability and distribution, can be found in reference (Sampson, 2008); however, the connection between fibres was not taken into consideration. Sastry and co-workers (Sastry et al, 2001; proposed a technique for modelling fibrefibre joints in which connection realized by a torsion spring can be regarded as flexible; however, the mechanical properties of fibrous networks with flexible bonding were not given.
It is suggested that the connectivity between fibres cannot be adequately described by single connection points in the beam modelling .
In this study, we have incorporated the density of intersections (or connections) into our model. This is important because the density of intersections is directly related to the relative density in stochastic fibrous materials. In our fibre network model, there is no fibre entanglement which is commonly found in woven fabrics (Grishanov et al, 2012). It is generally recognised that the macroscopic stresses and strains can be determined by the microscopic stresses and strains over a representative volume element (RVE). Hill (1963) proposed that a representative 'cell unit' (or a representative volume element) (RVE) can be a full-scale model to significantly reduce the computation complexity.
Traditionally the yielding of a material is defined according to the von Mises equivalent stress versus strain curve; however, it cannot be used to describe the yielding of a porous material when it is subjected to hydrostatic loading as the von Mises criterion is only based on the distortional part of the stored strain energy density. Some researchers have put forward the characteristic stress and strain which combine the hydrostatic energy density and deviatoric energy density to probe the yielding of two-dimensional isotropic honeycombs (Chen et al,5 1999; Chen and Lu, 2000), two-dimensional anisotropic cellular materials (Alkhader and Vural, 2009) and three-dimensional transversely isotropic foams (Ayyagari and Vural, 2015;Zhao et al, 2013). The deduction of characteristic stress and strain is based on the total stored strain energy density and is different from those phenomenological yield criterions, for instance, a shape parameter needs to be given to describe the mean-effective stress (Deshpande and Fleck, 2000).
It has been suggested that there exists a significant difference between different types of cellular materials, such as foams and honeycombs, and fibre network materials (Markaki and Clyne, 2003;Qiao et al, 2009;Zhao et al, 2013). It is crucial to build quantitative mechanical models to reveal the mechanical behaviour of stochastic fibre network materials. Due to the computation limitation, we have used a RVE to represent a whole fibre network material. The RVE thus must be periodic. It has been shown that periodic boundary conditions are more suitable than either mixed boundary conditions or prescribed displacement boundary conditions (Chen et al. 1999;Zhu et al, 2000Zhu et al, , 2001. To begin with, a three-dimensional periodic beam network model for stochastic fibre network materials has been constructed. We have introduced a novel method to deal with the connections between the fibres, in which additional beams are inserted to represent the connections (i.e. flexible cross-linkers). The elastic properties of the stochastic fibre network materials are investigated before proceeding to reveal their plastic behaviour. A scalar measure of characteristic stress and strain is utilized based on the total stored strain energy density irrespective of loading path. The objective has been to investigate the yield surface of 3D stochastic fibre network materials under multi-axial loadings and the dependence of elastic-plastic behaviour on the density of cross-linkers (i.e. connections) and the relative density of the fibre network materials. 6 a b Fig.1. The architecture of the stochastic fibre network materials: (a) 3-D reconstruction from the X-ray computed tomography (Tsarouchas and Markaki, 2011); (b) 3-D periodic random beam-network model.

Construction of three-dimensional stochastic periodic beam-network model
X-ray tomography has been utilized to extract the exact architectural characteristics from the very complicated porous fibre-network material and the data are used to reconstruct the geometrical model for this material Tsarouchas and Markaki, 2011) as shown in Fig.1a; however, this technique is very computationally challenging. Finite element models are an attractive approach to constructing the complex stochastic fibre-network structure based on fundamental geometrical parameters extracted from the x-ray computed tomography. In this paper, we have developed a code to construct a three-dimensional periodic stochastic beam-network model, as shown in Fig.1b, which shows a close similarity with the geometry from X-ray tomography.
In our study, the model generated is a three-dimensional stochastic fibre network structure with periodicity. In addition, we have properly incorporated the deformation mechanisms of the 7 cross-linkers into the stochastic fibre-network by inserting additional beam elements between the intersected fibres as shown in Fig.2. The establishment of the three-dimensional beamnetwork model is deduced progressively with the two-dimensional coordinates being considered first. A representative volume element (RVE) is adopted due to the computational limitation. The RVE must be periodic to meet the continuity and equilibrium between any two neighbouring RVEs. Thus, to construct a periodic stochastic fibre-network model in a two dimensional (i.e. the x-y plane) square domain of a side length , the x and y coordinates of the fibre centre, the fibre length L , diameter d and orientation  are all specified by independent random numbers which are generated by the computer and ranged from 0 to 1.0.
Both of the x and y coordinates are controlled to be in the range from 0.0 to 1.0, the fibre length is from 0.8 to 1.2 (the mean length is 1.0 L  ), the diameter is from 0.0086 to 0.0136 (the mean diameter is 0.111 d  ), and the orientation is from 0 to  . As shown in Fig.3, the parts of a fibre outside the square domain w w  are shifted into the domain by translating w in the x or y direction so that the RVE model is periodic in these directions.  The three-dimensional model is based on the two-dimensional one, but it requires a more complex design in terms of the intersections between fibres. All the information of the intersections between any two fibres is carried by the two-dimensional model, but there is a processing sequence of the fibres when we consider the growth in the out-of-plane dimension.  to be used to model the real structure, which would lead to enormous computational complexity even when the number of fibres is only 10 in the model. Alternatively, polylines could be used to represent the curved lines. We examined 10000 segments in our polyline model, the mean angle of the inclined segments with the x-y plane,  , is around  5 . 8 , and the relative error with 95% confidence associated with this mean angle is below 5%. The mechanical behaviour of a  Table 1. Note that the simplification by using polylines to represent the curved fibres does not cause significant error in predicting the performance of the fibre network model, moreover, the innovative modelling method optimises the simulation with both efficiency and accuracy.
In the numerical model, we have introduced an innovative way to incorporate the deformation mechanisms of the cross-linkers into the computer simulation. As shown in Fig process, if two metal fibres are partially overlapped because of, for instance, high temperature, the middle part of the inserted beam is usually thinner than the two end parts. If two fibres intersect at a right angle, the minimum diameter of the intersected part (which is in the middle of the inserted beam) is . The dependence of plastic behaviour of the fibre network materials on the dimension of the cross-linker will be elucidated in section 5.
The way to achieve periodicity in the z direction for a fibre-network model with N complete fibres is shown in Fig.5. After having built up a stochastic three-dimensional fibre-network model with N complete fibres, another N fibres are continuously put up in the z direction, whose x and y coordinates are identical to those of the previous N fibres correspondingly. Ideally the shape of the (i+N)th fibre is the same as that of the ith fibre when N is sufficiently large.
When taking the top layer fibre-network model with N complete fibres for simulation, there are some isolated cross-linkers (i.e. inserted beam elements) indicated by the blue dotted lines.
These isolated cross-linkers are moved up and the out-of-plane periodic boundary conditions are applied on their nodes. The relative density is a key parameter to elucidate the mechanical behaviour of porous materials including fibre network materials Zhu et al, 2001Zhu et al, , 2000Zhu et al, , 1997. In our study, we have conducted dimensional analysis to reveal the relationship between the yield strength and relative density of the stochastic fibre network materials as detailed in section 4.
By calculating the average distance between the constrained top and bottom nodes in the z 13 direction, the thickness t of the fibre-network model is obtained. Based on the thickness, the relative density of the fibre network materials is specified by: To obtain more intersections or cross-linkers with other fibres, each fibre has to drop down so as to connect with more fibres below. This consequently reduces the thickness of the structure and thus increases the relative density. As shown in Fig.6, the relative density of the fibrenetwork materials is linearly related to the density of the cross-linkers by

Mesh and boundary conditions
The fibre-network model is partitioned into a large number of 3-node beam elements (ANSYS   Table 3, and the results are the mean values of twenty paired stochastic models. It is worth noting that when the number of complete fibres increases, the ratio of the reaction force of two layers to that of one layer approaches 2.0. When N=200, the standard deviations of the relative density, the Young's modulus in the x direction and the reaction force are all very small. Also, with the consideration of computational efficiency, the finite element results presented below are the mean results over twenty stochastic models having the number of complete fibres fixed at N=200. The relative error with 95% confidence interval is 0.75% for relative density and 2.9% for Young's modulus.

Elasto-plastic behaviour
One of the significant features incorporated into our three-dimensional beam-network model is the anisotropic elasticity of the fibre networks. It is shown that the mean values of Young's modulus and the Poisson's ratio for 20 models are almost identical in the x and y directions, which suggests that the stochastic fibre network materials are transversely isotropic (i.e. isotropic in the x-y plane).

Characteristic stress and strain
Chen and Lu (2000) introduced characteristic stress and strain for isotropic foam, based on which, a yield criterion has been developed by hypothesizing that the yielding is driven by the total stored strain energy density, including the hydrostatic and deviatoric parts. Uniaxial response of a transversely isotropic material ) ( 2 1 E E  can be schematically represented by the stress-strain graphs in Fig.7. To obtain a scalar measure of the strain energy density independent of the loading directions, the normalized strain energy density can be written as, X are the yield stresses and strains along the x and z directions for the transversely isotropic fibre-network materials subjected to uniaxial loadings, and 13 Y , 12 Y , 13 X , 12 X are the shear yield stresses and strains along the 1-3 and 1-2 directions as shown in Fig.7. 18 Fig.7. The elastic-plastic response of transversely isotropic materials under uniaxial loading and pure shearing. By scaling back all the normalized components of the stress and strain tensors, the characteristic stress and strain are independent of the loading path, and given as, which all the characteristic stress and strain curves collapse along a master line whose slope is Ê . Note that the characteristic stress and strain carry the information on stiffness and strength anisotropy, such as 11 Y ) is taken as the stress at the point on the stress-strain curve at which the tangent gradient is 50% of the initial gradient (i.e. the Young's modulus or shear modulus) of uniaxial loading or pure shearing (Johnson, 1918).

Yield strength and yield surface
The homogenized total strain energy density It is hypothesized that yielding of the stochastic fibre network structure occurs when the homogenized total strain energy density Wˆ reaches a critical value. Under the stress state of uniaxial loading in the x direction, yielding occurs when Therefore, the yield criterion of the stochastic fibre-network materials under any arbitrary state of stresses is given as 2 11 From Eq. (7) it can be seen that the scaling back direction can affect the size of the yield surface (i.e. the amplitude of the right hand side of Eq. (7)

Characteristic strain offset value
Characteristic stress-strain curves can be used to determine the yield point for stochastic fibre network materials in a unique manner even under different states of stresses as all the characteristic stress-strain curves fall on the same line before yielding. The traditional yield strain for solid mild steel is 0.2% in experimental studies (Chen et al., 1999;Zhao et al., 2013).
In this study, the yield strain is taken as the strain at the point on the stress-strain curve of uniaxial loading in the x direction at which the tangent gradient is 50% the initial gradient (i.e. the initial Young's modulus 1 E ). Thus, the yield strain is 0.136% for stochastic fibre-network materials with a relative density % 16 . 8   . Note that the characteristic yield strain or stress is independent of the state of the stresses, thus the value is the same as those obtained by using the yield points of either the uniaxial loading in the z direction or pure shearing along the 1-3 or 1-2 directions.

Multiaxial normal stress states
With the energy-based yield criterion, the yield function can be fully calibrated in terms of the uniaxial tension (or compression) response in the x or z direction, or by pure shearing along the 1-3 or 1-2 directions instead of complex multiaxial loading responses.
The characteristic stress-strain curves from the FE simulations for various biaxial strain paths in both the plane of isotropy and the plane of anisotropy are plotted in Fig.9. All the data points 23 on the graphs are from the same stochastic fibre network structure with the relative density % 16 . 8   , the density of cross-linkers / 10 c Ll , the fibre aspect ratio / 90 Ld and the number of fibres N=200. As can be seen from Fig.9, for biaxial loadings in both the plane of isotropy (x-y plane) and plane of anisotropy ( 1  , 2 : 1 , show significant similarities. It is mainly attributed to the much larger strength and stiffness in the x direction than that in the z direction.
As the properties of the structure are dominated by the axial loading in the x direction in the plane of anisotropy, a slight change in loading of the z direction hardly affects the plastic response of the stochastic fibre network structure. The yield surface of transversely isotropic stochastic fibre network materials obtained from FE simulations under multiaxial loadings is plotted in Fig.8 in the space of characteristic effective stress and mean stress. As the strength and stiffness in the z direction are smaller than those in the x or y direction, the sign of effective stress of the fibre network materials in the multiaxial loadings is always the same as that of the 24 stress in the x or y direction rather than the z direction, e.g., 10 : 1 : 1 : as shown in Fig.8.

Combined axial loading and shearing
It is noted that the presence of the shear stress component does not affect the form of the yield criterion given in Eq. (8). In the plane of isotropy (i.e. the x-y plane), we do not need to consider the effect of the shear stress on the yield surface (i.e. Fig.8) because the principal stresses in this plane can always be obtained. In the plane of anisotropy (i.e. the x-z or y-z plane), the presence of the shear stress component contributes to the amplitude of the normalised characteristic deviatoric stress in Eq. (8), thus, it reduces the size, but does not change the shape of the yield surface, as can be seen in Fig.10 where the unit of the stress is MPa. When the shear stress is absent, the size of the yield surface in the plane of anisotropy arrives at its maximum. With the increase of the shear stress, the size of the yield surface decreases.

Effects of relative density on the uniaxial yield strength
The effects of the density of cross-linkers / c Ll on the elastic constants and yield strengths of random fibre-network materials can be evaluated from the relationships between the relative density of the fibre network and its stiffness and strength. The relative density can be adjusted by changing the density of cross-linkers (i.e. the number of intersections of a fibre with those below it) as the relative density of a fibre network material has a linear correlation with the density of cross-linkers (see Fig. 6).

Simulation results
In practical applications, one of the most important properties about the fibre network materials is the relationship between their yield strength and their relative density. To obtain the uniaxial yield strength, we first obtain the uniaxial tensile stress strain relationship for stochastic fibre network materials with different relative densities. The yielding is supposed to occur at the point on the uniaxial tensile stress-strain curve at which the tangent gradient is 50% of the elastic constant (i.e. the initial Young's modulus) (Johnson, 1918). a b Fig.11. The relationship between the relative density and the uniaxial yield strength of stochastic fibre network materials in the x direction (a) and in the z direction (b).

26
The finite element simulation results indicate that the uniaxial yield strength in the x direction is a quadratic function of the relative density of the stochastic fibre network materials when the relative density is small (i.e. 0.2   ), as shown in Fig.11(a). In contrast, the uniaxial yield strength in the z direction is constantly a cubic function of the relative density, as can be seen in Fig.11 A shear deformation is enforced by imposing a pair of shear strains in transverse directions.
The shear yield strength is taken as the stress at the point on the stress-strain curve where the tangent gradient is 50% of the shear modulus of pure shearing (Johnson, 1918). The effects of relative density on the shear strength of stochastic fibre network materials are shown in Fig.12.
The in-plane shear strength 12 Y is a quadratic function of the relative density when 0.2   (see Fig.12(a)), and it gradually becomes a linear function when 22 . 0   (see Eq. (13)). The outof-plane shear strength 13 Y of the fibre network materials is much lower than their in-plane shear strength 12 Y . Experimental tests of MFSSs suggest that the in-plane and out-of-plane shear strengths have linear and cubic dependence upon the relative density, respectively (Zhao and Chen, 2014). The simulation of our 3-D random beam model is broadly consistant with the 27 experimental results, especially a cubic function in the transverse direction as can be seen in Fig.12(b).

Results of dimensional analysis
To perform dimensional analysis, one has to simplify and/or idealise the geometrical structure of the fibre network materials. As the stochastic fibre network materials are transversely isotropic with the properties in the x direction the same as those in the y direction, the geometrical model for dimensional analysis is simplified, as shown in Fig.13(a). This model is periodic in all the x, y and z directions.
where d is the mean diameter of the fibres, 1.9 cx P l d  and x  is the effective uniaxial tensile stress applied to the fibre network materials in the x direction. Note that the distance between the centroidal lines of two intersected fibres is taken as 0.95d in the model.
The dimensional analysis results from Eq. (13) are presented in Fig. 11(a) for comparison with the simulation results and experimental results (Zhao et al, 2016), showing excellent agreement.
Eq. (13) suggests that when the relative density is small (i.e. 0.2   ), the in-plane tensile strength 11 Y is a quadratic function of  ; while when the relative density is large (i.e. 0.22 the denominator is dominated by the term 2 20.89 and the in-plane tensile strength 11 Y gradually becomes a linear function of  .

Uniaxial yield strength in the z direction
In the random beam model developed in this paper, the cross-linkers (or intersections) are represented by inserted beam elements which are in parallel to the z direction. When the 30 simplified model is stretched by a uniaxial tensile stress z  in the z direction, the contour of von Mises stress in the idealised fibre-network materials is shown in Fig.13 As the full bending capacity of a uniform beam with a circular cross-section is 1.7 times the bending moment at initial yielding, the uniaxial tensile strength (or yield strength) of the fibre network materials in the z direction can be approximated as The dimensional analysis results from Eq. (16) are presented in Fig.11(b) for comparison with the simulation results and experimental results (Zhao et al, 2016), confirming that the difference is just a constant factor and fibre bending is indeed the dominant deformation mechanism.

31
The results in Fig.11(a) and Fig.11(b) demonstrate that the above dimensional analyses obtained from a simplified model have correctly captured the main deformation mechanisms when the stochastic fibre network materials undergo uniaxial in-plane or out-of-plane tension.
It is noted that the dimensional analyses assume that the relative density  is the same or uniform over the whole model of the fibre network materials, while in reality, the local density could be double or half the average relative density at different locations. When the fibre network materials undergo uniaxial tension in the x direction, the beam/fibre-network can be treated as a number of springs (with different stiffnesses) which are deformed partly in parallel and partly in series. As 11 Y is proportional to 2  when 0.2   and to  when 0.22   , it is less affected by the non-uniform distribution of the relative density. When the fibre network materials undergo uniaxial tension in the z direction, the beam/fibre-network can be treated as a number of springs (with different stiffnesses) in series. As 33 Y is proportional to 3  , the yield strength of the fibre-network materials is mainly dominated by the locations where the relative density is much lower than the average. That is why the dimensional analysis results are closer to the simulation results in the x direction, but much larger than those in the z direction. The strong agreement from the comparison, as shown in Fig.11, indicates that our fibre network models are supported by experimental evidences. Jin et al (2013) have simulated the in-plane yield behaviour of fibre network materials and suggested that the deformation mechanism of porous metal fibre sintered sheets (MFSSs) is stretching dominated. The difference between their results and our findings is attributed to the geometrical models. For a 2-D beam model, the density of cross-linkers (intersections) is certainly significantly overestimated and this tends to predict a larger in-plane yield strength 11 Y .