1 Introduction
At the microscale, a paper sheet consists of a network of fibres that is produced from wood pulp as shown in Fig. (a)a. The paper fibres have a preferential orientation (machine direction) due to the manufacturing process, which results in the observed anisotropic behaviour (Sabit, 1997). In the network, fibres are bonded with each other in certain regions and free standing elsewhere. Upon exposure of a sheet of paper to a humid environment or liquid water, moistureinduced swelling takes place, which is called hygroexpansion (Larsson and Wagberg, 2008). In practical applications related to printing, this results in macroscale effects such as curling, waviness, and cockling. At the microstructural level, the changes in each individual fibre due to the variation in moisture content are transmitted to the neighbouring fibres by the bonds in the fibrous network as sketched in Fig. (b)b. These changes include geometrical variations of the shape in each fibre and the accompanying stress concentrations induced by the bonded areas. Understanding these phenomena and their dependence on the properties of the individual fibres and the network geometry is essential for predictive modelling of sheetscale phenomena in paper as well as other networks of natural fibres.
Gaining insight in the behaviour of a fibrous network subjected to hygroscopic swelling is essential to unravel the influence of different properties of the fibres and network on the overall hygroexpansive response. Among the early attempts to model the mechanical response of paper, Cox (1952) studied the effect of orientation of fibres on stiffness and strength of paper by assuming that fibres carry only axial forces. The flexural stiffness of the fibres was also taken into account (Akker, 1962). The deformation of bonds and the elastoplastic behaviour of fibres was also studied (Ramasubramanian and Perkins, 1988). Assuming constant strain, the transverse properties of fibres were included (Schulgasser and Page, 1988)
. Later, a model incorporating shear forces, axial and bending and torsional moments between rigid bonds was developed for a network of fibres
(Starzewski and Stahl, 2000). Interfibre bonding was also considered in modelling the fibrous network (Bronkhorst, 2003). Some of the recently developed numerical models for fibrous networks mostly modelled the fibrous material with twodimensional assumptions considering the individual fibres as trusses or beams with only isotropic properties (Kulachenko and Uesaka, 2012; Shahsavari and Picu, 2013; Dirrenberger et al., 2014). Most of these works were carried out for the description of the mechanical behaviour of the fibrous network, without any coupling to the hygroscopic response.Several other studies addressed the hygroexpansivity and dimensional stability of paper. Some studies focused on the relation between hygroexpansion and moisture content of paper in order to understand the factors affecting them (Nordman, 1958; Salmen and Olsson, 2014; Uesaka and Qi, 1994). Also, studies on the dimensional stability of paper at the macroscale were carried out (Page and Tydemann, 1961; Uesaka, 1994; Niskanen, 1998; Niskanen et al., 1997; Larsson and Wagberg, 2008; Erkkila et al., 2015; Bosco et al., 2015b). Several other works were done on studying the hygroexpansion of paper fibres. However, most of these studies lack in modelling the hygromechanical behaviour of complex fibrous networks like paper subjected to moisture fields.
The aim of the present study is to model the hygromechanical behaviour of fibrous materials through a multiscale analysis of the network using periodic homogenization (Guedes and Kikuchi, 1990; Boutin, 1996; Peerlings and Fleck, 2004). The fibres are modelled as twodimensional ribbonlike elements in a network subjected to hygroexpansion, for which a finite element framework was used in (Bosco et al., 2017). This study used a regular grid of triangular finite elements to model the network which inevitably is not aligned with most of fibres. The use of such nonconforming finite elements does not allow to accurately capture the geometry of the fibres and of the voids and bonded regions between them. Nonconforming FEM considers any fibre lying inside the centroid of a finite element to contribute to its stiffness. This leads to a geometrical representation of a fibre in the network with jagged boundaries. To better resolve the fibre boundaries, a very fine discretization was therefore employed in (Bosco et al., 2017), which increases the computational effort. A conforming discretization would avoid these limitations. However, in networks of realistic complexity it will be very hard to generate a conforming mesh for all fibres in the densely bonded regions. An example of such a network is shown in Fig. (a)a, where a paper sheet is represented as a collection of ribbonshaped rectangular fibres with coverage , where the coverage is defined as the ratio of the total area occupied by the fibres in the network to the area of the microstructural unit cell. The density of fibre intersections and overlapping regions is more clearly visible in the magnified view in Fig. (b)b. The main reasons rendering a conforming triangulation infeasible are: (i) finding the geometrical intersections of individual fibres and their overlapping regions requires a large number of operations and hence will be too expensive (if possible at all); (ii) meshing the decomposed geometry will result in an excessively fine triangulation, rendering the subsequent analysis of the network model untractable. Hence, there is a need to develop a methodology that captures the fibre boundaries with a better accuracy than a nonconforming mesh and which is less computationally demanding than a conforming discretization strategy.
Here, an advanced discretization scheme is combined with the hygromechanical model developed by (Bosco et al., 2017) to enable modelling of larger systems in twodimensional configurations, with the potential to be extended towards three dimensions. The levelset formalism is used in combination with the Extended Finite Element Method (XFEM) to capture the geometrical complexity of the network—here in two dimensions—in an efficient and inexpensive manner, providing the decomposed geometry and mutual intersections of individual fibres. The levelset methodology is a mathematical tool to describe complex and timeevolving boundaries (or, more generally, geometry) implicitly (Sethian, 1999). The boundary of a fibre is represented by the zero levelset of a higher dimensional function. This results in a versatile geometrical description that is decoupled from the spatial discretization. XFEM (Daux et al., 2000)
allows to account for the effect of interfaces (boundaries of fibres) on the mechanical behaviour of the problem. In the bonded regions, the interpolation functions classically used in the elements are modified by discontinuous enrichment functions, so that geometrical discontinuities associated with the fibre interfaces can be resolved. This allows to capture the fibre’s boundaries in the bonded regions accurately using reasonably coarse meshes, by coupling of the levelset functions with the XFEM enrichment. For realistic networks, which are as complex as the one shown in Fig.
6, the level set technology thus provides a feasible numerical tool for the description of their geometry.This paper is organized as follows: In Section 2, the geometry of the fibre network model is discussed briefly, along with the levelset formalism used to represent it. In Section 3, the XFEM discretization used for the fibrous network is described, together with the required specific numerical integration scheme to accommodate it. Section 4 presents the simulation results of the hygromechanical behaviour of paper using the XFEM enrichment, first for the simplified models presented in (Bosco et al., 2015b), next for somewhat more complex networks with different coverages, and finally for a complex realistic network. The local and global responses are analysed for both mediumcomplexity and complex networks. Finally, Section 5 reports the conclusions and perspectives.
Throughout this contribution, the following notations are used for operations on the Cartesian tensors. Scalars, vectors, and tensors are denoted by
, , and respectively. The th order tensors are represented by . For tensor and vector operations, the following equivalent notations are used with Einstein’s summation convention on repeated indices: with ( for the global reference system and for the local reference system). The Voigt notation used to represent tensors and tensor operations in a matrix format is depicted as follows: and denote a column matrix and a matrix of scalars respectively. The matrix multiplication is denoted as .2 Fibre network model
2.1 Fibre model
The fibrelevel constitutive model formulated in (Bosco et al., 2015b) is adopted here. A twodimensional plane stress model is assumed and the fibres in consideration are oriented in the plane of the global reference frame, and subjected to a uniform moisture change. For the fibre constitutive model, a local reference frame is considered along the directions of the fibre as shown in Fig. 7.
The hygromechanical properties of the paper fibres are assumed to be transversely isotropic with respect to their longitudinal axis. The general elastic constitutive law for a fibre, assuming plane stress state in the direction, and subjected to a moisture change is expressed as
(1) 
where the hygroscopic strain is given by
(2) 
In these expressions, , , and are the elastic constitutive tensor, the strain tensor and the hygroexpansivity tensor of the fibre respectively. In matrix notation, and are represented as
(3) 
In Eq. (3), and denote the elastic moduli in the longitudinal and transverse direction with respect to the fibre axis, is the inplane shear modulus, and and inplane Poisson ratios. The coefficients of hygroscopic expansion, and , are different in the longitudinal and transverse directions of the fibre.
In the local reference system, a fibre is oriented at an angle with respect to the global reference system . Therefore, the relationships of Eqs. (1)–(3) need to be transformed from the local (fibre) frame to the global frame of the paper material for each fibre in the network (Roylance, 1996). The fibre bonds are important because of their influence on the overall behaviour of the network. In the 2D modelling, the fibres are assumed to be perfectly bonded to simulate the interplay between hygroexpansion and elasticity of the fibres. This implies full displacement strain compatibility inside a bond for all the fibres interconnected in that bond.
2.2 Levelset formalism
The levelset formalism (Sethian, 1999; Osher and Fedkiw, 2003) is here used to describe the geometry of fibres. In this method, the fibre is described as the zero levelset of a higher dimensional levelset function (see Fig. 8).
In most cases, including the current work, the levelset function gives the signed distance of a point to the interface of the rectangular fibre denoted by in Fig. 8. The levelset function is defined as
(4) 
where the sign is negative if is outside and positive if it is inside the contour defined by .
2.3 Random fibre network creation
In order to understand the hygroscopic behaviour of a complex network of fibres, a set of rectangular fibres having a length and width is randomly generated in a unit cell of length . Each of the fibres is generated with random coordinates for its centroid as shown in Fig. (b)b. Now, depending on the number of fibres generated in the unit cell, a coverage is defined as the ratio of the total area occupied by all fibres and the area of the unit cell.
The orientations of the randomly generated fibres in the network satisfy a probability density function based on
(Cox, 1952):(5) 
where is the angle between the fibre and the machine direction, with . Fig. (a)a shows the probability density function for the orientation for two values of the anisotropy parameter : for the isotropic case , and an arbitrary degree of anisotropy (Bosco et al., 2015b).
Once the fibres are generated and each of them is assigned a levelset function, , , a periodicity condition is enforced on the unit cell. The levelset function that allows identifying the fibres and voids for the entire network can be described as
(6) 
where in fibres and in voids.
3 Discretization strategy
3.1 XFEM methodology
The extended finite element method (XFEM) is a numerical discretization method developed to model material discontinuities, singularities, and moving boundaries independently of the underlying finite element mesh. With the microstructures considered in this work, it is not feasible to produce a geometryconforming mesh for the free standing parts of fibres and even more in the bonded regions where fibres partially overlap. Accurate local results cannot be obtained by the conventional finite element method (FEM) if the mesh does not conform to the boundaries of the fibres. XFEM can be used instead to capture the fibre boundaries both in free standing and bonded regions accurately even with a regular nonconforming mesh. In terms of accuracy, an adequate solution can be obtained by enriching the classical shape functions with special functions capturing the geometrical discontinuities in XFEM.
Depending on the nature of the problem under consideration, a suitable enrichment function needs to be adopted in XFEM. To model the geometrical discontinuity between solid (fibres) and voids in a fibrous network with a regular mesh, a Heaviside function defined in terms of the levelset function is adopted (Sukumar et al., 2001) as an enrichment function, since the geometry is already described by the levelset function . The displacement interpolation then reads (Daux et al., 2000, Section 4)
(7) 
in which is the number of all elements of the underlying mesh and is the nodal displacement column of the finite element , for inside fibres and for inside void regions.
The interconnection between XFEM (to model discontinuities) and levelset (to capture the geometry), referred to as LSXFEM in what follows, is discussed in detail in the next section. A LSXFEM formalism that makes use of the levelset information in this particular XFEM setting thus greatly simplifies the set up of the model under consideration.
3.2 Discretization with the levelset and XFEM formalism (LSXFEM)
In the current work, the behaviour of the fibres is assumed elastic and small strains and displacements are considered. The total potential energy in the unit cell consisting of network of fibres is given by
(8) 
where is the total number of fibres occupying the periodic cell area , is the signed distance function of the th fibre, the strain tensor associated with Eq. (7), the hygroscopic strain in fibre , and its thickness. This expression may be written in matrix format as
(9)  
in which and are the elastic constitutive matrix and hygroscopic strain matrix transformed to the global frame, for a particular fibre aligned at an angle with respect to the global frame as discussed earlier. Using the straindisplacement relationship for each finite element, , where is the standard straindisplacement matrix, the integration over the entire network of fibres in Eq. (9) is carried out element wise. Swapping the summations over all finite elements with areas and the fibres, the assembled potential energy is obtained as
(10)  
Of all the possible displacements that satisfy the boundary conditions of an elastic structural system, the ones corresponding to the equilibrium configuration minimize the total potential energy, i.e.
(11) 
or, using Eq. (10),
(12) 
This equation is of the form of a linear system
(13) 
in which collects all nodal displacements, and stiffness matrix and the hygroscopic load vector are given by
(14)  
(15) 
Linear triangular finite elements are used for the discretization. Three cases can be distinguished for mapping fibres onto a finite element. First, if the finite element lies entirely in a void, the Heaviside function in Eqs. (14)–(15) vanishes and the element does not contribute to the stiffness matrix and hygroscopic load vector. If the finite element lies entirely in one or more fibres, its full contribution is accounted for in the computation of the stiffness matrix and hygroscopic load vector. Finally, if the finite element is intersected by fibre boundaries, a specific integration technique needs to be used, as described next.
3.3 Numerical integration scheme
In this section, the integration scheme used for elements intersected by fibre boundaries is detailed. The integration over a particular finite element for all fibres passing in it is considered for the stiffness and force vector of that element, recall Eqs. (16) and (17),
(16)  
(17) 
In the computation of the contribution of each fibre to the finite element, it is therefore only required to determine the area occupied by each fibre, i.e. in this particular element, where . The algorithm used to compute this area of each fibre in the finite element is illustrated by considering two fibres lying in the region of a triangular element with vertices as shown in Fig. 16. Using the levelset function for a particular fibre , if one of the vertices of the finite element is identified to lie inside the fibre, then the element is considered to lie in the fibre. Therefore, the triangle is subdivided along its longest edge. The subtriangulation continues and stops only when a tolerance is met in terms of the smallest triangle area. This procedure is executed for all fibres lying partially inside the finite element. The integration scheme can be summarized as follows:

Check for the first fibre if all element vertices lie in it (i.e. ). If the particular fibre lies partially in the element, then the element is split into two new triangles as shown in Fig. (b)b.

The subtriangulation is terminated when a specified tolerance limit in terms of a minimum triangle area is achieved.

After the termination of subtriangulation, the area occupied by the fibre in the finite element is sum of the areas of triangles lying completely in the fibre and the ones whose centroid lies inside the fibre. Thereafter, the above steps are repeated again for the next fibre.
The element subtriangulation for numerical integration is adopted because of its ability to control the integration accuracy in relation to the computational cost. The main reason here is that a conforming integration subtriangulation for all fibres might become relatively expensive, and for complex systems as the one shown in Fig. 6 it may reach a complexity that is comparable to creating a fully conforming triangulation.
3.4 Mesh refinement
Considering a unit cell consisting of a network of fibres with complex topology, a coarse triangulation is naturally insufficient for an accurate description of the interfaces. Using a globally fine triangulation would however dramatically increase the computational cost. Therefore, a refinement strategy based on the backward longest edge bisection algorithm (Rivara, 1997) is adopted here for elements intersected by boundaries. This approach preserves adequate accuracy at a lower cost, still allowing for large triangles in regions with little strain fluctuation within fibres and voids. This mesh refinement should not be confused with the subtriangulation discussed in the previous section. The subtriangulation is used to accurately capture the geometry, i.e. the volume of an element occupied by a particular fibre. Mesh refinement, on the other hand, improves the accuracy of the kinematics of the problem.
After the generation of the geometry of a network of fibres, the finite elements located at the boundaries of the fibres are revisited for mesh refinement. For each such element, the corresponding fibre signed distance functions are evaluated at all of its vertices of considered triangular finite element. When all three vertices are contained inside all fibres passing through it, or when none of them belong to the fibres, the element is not categorized as a boundary element. In all other cases, the element is identified as boundary element to be refined, for which a mesh refinement algorithm is applied.
As defined by Rivara (1997), for a triangle of any triangulation , the longest edge propagation path (LEPP) is an ordered list of triangles , such that is neighbour to along its longest edge. The (LEPP) terminates with (a) one triangle with its longest edge along the external boundary of the mesh, or (b) a pair of triangles sharing the same longest edge. Now, based on the LEPP of triangle marked for refinement, a backward longest edge refinement algorithm is used in which the longest edge of is bisected in the former case and both triangles and are bisected along longest edge for the later case. The LEPP is updated and the procedure is repeated until the initial triangle is bisected as well. In Fig. 22, the longest edge bisection algorithm is illustrated for a simple triangulation.
In the current work, the elements at the interfaces between a free standing fibre segment and void, as well as at the boundaries of bonds are refined. The newly generated refined triangular elements at these locations also better capture multiple fibres passing through the finite elements. Therefore, the accuracy of the mesh improves in representing the fibre geometry and, more importantly, its kinematics in these complex networks. In addition, the boundary edges of the unit which is periodic in nature are also refined to maintain periodicity of the nodes.
However, the refined mesh, as seen in Fig. 25 will not be conforming to the fibre boundaries. As mentioned earlier, the high number of fibre interfaces particularly in the bonds, would render it difficult to capture with a conforming mesh at cheaper computational effort. This again emphasizes the motivation to use the proposed LSXFEM formalism, which allows to decouple the mesh and fibre boundaries with a structured mesh that still captures the fibre edges adequately.
4 Results and discussion
To demonstrate the benefits of the LSXFEM formalism for fibrous networks, simplified illustrations are first presented below. Subsequently, the formalism is used for mediumcomplexity and complex realistic networks to study the effect of microstructural features on the sheetscale properties due to moisture infiltration.
4.1 Simplified networks
First, the stress level in a single family of parallel fibres subjected to a tensile load is assessed, followed by the study of the stress concentrations in the bond regions of an elementary network consisting of two families of fibres subjected to a macroscopic tensile load.
4.1.1 Parallel fibres subjected to uniaxial tension
In this problem, an infinite number of parallel, equispaced, infinitely long fibres is considered, as shown in Fig. (a)a. The fibres are subjected to a horizontal stress, . Using a square periodic unit cell of length for this simple geometry, half of a fibre occurs at the top and the another half at the bottom of the cell, see Fig (b)b. The chief reason to use this unit cell in the example is that there are regions with fibre and voids, which may also be the case in the complex fibrous network considered later. Also, in the discretized unit cell, the mesh is not conforming to the geometry of the fibre, i.e. the same situation occurs in a complex fibrous network, and therefore we will see the advantages of the LSXFEM formalism as compared to the nonconforming FEM. The main reason for using a nonconforming mesh for standard FEM is to measure the error induced by it. Although a mesh conforming to the considered geometry could easily be obtained for this simple system, it is not possible for much more complex systems, as discussed in Section 4.3 below. The input parameters for the anisotropic fibres are an elastic modulus in the longitudinal direction, , and in the transverse direction, , shear modulus, , Poisson ratios and . The fibre has a width and thickness . The exact area of the fibre is .
The problem is formulated as a plane stress case and solved with a nonconforming fixed grid by the standard FEM and the LSXFEM approach. In the standard finite element setting, a finite element is considered to be part of the fibre if the centroid of the finite element is located inside the fibre. Accordingly, these fibres contribute to the finite element stiffness in the FEM solutions, and a slow convergence rate with respect to the chosen element size of the underlying discretization is thus expected for the nonconforming FEM. In the LSXFEM formalism, a finite element partially covering a fibre still represents that fibre, but the area of the fibre in that finite element is determined using the area integration method presented in Section 3.3 and it contributes to the element stiffness accordingly. The results obtained by these two approaches for different grid spacings (element edge lengths) are presented in Tab. 1 in terms of the (local) fibre stress component as computed, and the fibre area considered by the numerical model.
As expected, the nonconforming FEM does not yield a good estimate of the stresses with a coarse mesh of 10 elements along the cell edge. On the contrary, the LSXFEM formalism yields an accurate stress distribution with the same mesh. As the mesh is refined, the nonconforming FEM still struggles to yield the expected accuracy, showing a slow convergence rate. However, the LSXFEM formalism captures the geometry with a better accuracy, which also results in a more accurate prediction of the stress level.
4.1.2 Bonded elementary network under tension
The attention is next shifted towards the simplest fibre arrangement that involves bonds between two orthogonal families of parallel fibres subjected to a stress , as introduced in (Bosco et al., 2015b), see Fig. (a)a. Because of periodicity, the unit cell depicted in Fig. (b)b is used, with half fibres along the edges. The material parameters are identical to those of the previous example. The LSXFEM formalism, with a regular mesh with as shown in Fig. (b)b, is used to predict the mechanical response, especially focusing on the stress distribution in the bonded regions. This network is also modelled with a very fine FE mesh of size elements conforming to the geometry, serving as a reference solution. It is emphasized that such a conforming mesh is of course only achievable for simple configurations as considered in this example, and is used here mainly to quantify the accuracy of the LSXFEM method.
Figs. 34 and 37 show the normalized stress distribution in the bonded regions for the horizontal and vertical fibres as computed by the LSXFEM with a nonconforming mesh and the conforming FEM reference solution. The coarse triangular mesh used for the LSXFEM formalism is evident in the jagged stress distribution in the free standing and bonded regions in Fig. (b)b. Also, these are noticed near the bonded region of the vertical fibres in Fig. 37. Note that only at locations (), close to the edges of the bond, the stress distribution predicted by LSXFEM differs from the reference solution with a deviation of 20, as observed in Fig. (a)a. This is because with the LSXFEM formalism employing a geometry nonconforming coarse mesh, the finite element at this location lies in the bond as well as in the free standing fibre. Now, the stress distribution is high in the free standing fibre resulting in the stress jumps obtained by LSXFEM. However, in the geometryconforming mesh employed by the FEM solution, the finite element is lying inside a bond with low . Later on, the stress distribution by the FEM also attains the same value as predicted by LSXFEM formalism. Most importantly, as noticed in these plots, the normalized stress distributions predicted by the LSXFEM formalism are qualitatively and quantitatively similar to the reference solutions at other regions, which is further ascertained in the crosssectional plots through the bonds shown in Fig. 40.
In a fibrous network, the bonded regions are vital for an accurate prediction of the overall response, as well as for the proper reproduction of the local behaviour of the fibres (Bosco et al., 2015b). Therefore, the ability of the LSXFEM formalism to make adequate predictions of the mechanical stress state inside the bond at a lower computational cost than a fine conforming FEM discretization makes it a suitable tool for modelling fibrous networks.
4.2 Mediumcomplexity network
The attention is now focused on somewhat more complex networks, to illustrate the ability of the proposed formalism to recover information at both the microstructural and macroscopic level. The error committed by the LSXFEM is quantified by comparison with a reference, finely discretized conforming FEM. The complexity of the networks considered in this section is limited by our ability to generate (and simulate) such a conforming discretization—see Section 4.3 for a complex realistic network, for which this is no longer feasible (but where our approach still works).
4.2.1 Network parameters
Different network configurations are considered, with coverages of and . Recall that coverage is defined as the ratio of total area occupied by the fibres in the network to the area of the microstructural unit cell. The characteristic size of the finite elements used in the mesh to model the unit cell is chosen as . After application of the refinement strategy at the fibre edges, the smallest finite element size reduces to . For the anisotropic behaviour of the fibres, the material parameters used are identical to those of the previous examples. The coefficients of hygroscopic expansion are taken according to for all cases. A unit change in moisture content, , is adopted which is assumed uniform over the entire unit cell.
4.2.2 Average expansivity and deformed geometry
The initial anisotropic network for coverage and anisotropy parameter is shown in Fig. (a)a. The hygromechanical response of the unit cell is computed by solving the static equilibrium problem for a unit change in moisture content, . This generates a hygroscopic load causing deformation in the network, which is computed by means of the LSXFEM formalism on a discretization with four levels of refinements, see Fig. (b)b. The response of the same network is computed as well using a conforming FEM with the maximum element size , the discretization of which is shown in Fig. (c)c. For the conforming discretization, the Gmsh mesh generator, version 4.5.6, has been used (Geuzaine and Remacle, 2009). It can be observed in Fig. 47 that the deformed geometry obtained by the LSXFEM (with a relatively coarse mesh), in Fig. (b)b, is similar to the deformed network obtained by the conforming FEM (with a very fine mesh), in Fig. (a)a.
Approach  
conf. FEM  
LSXFEM 
This is further illustrated by the overall behaviour on the basis of the computed overall hygroscopic coefficients of the network, as listed in Tab. 2. The anisotropic network fabric () causes a pronounced overall anisotropy: the expansivity in the cross direction exceeds that of the machine direction, , by more than a factor of three. This is due to the fact that the fibres are on average more oriented in the machine direction, i.e. their expansion occurs predominantly in the cross direction. The values obtained by the LSXFEM have a relative deviation of less than for and for , when compared against the conforming FEM solution for the same network. In Fig. 48, convergence of the normalized effective hygroexpansivity of the network with an increasing number of elements in the mesh for the LSXFEM formalism is shown. There is no significant change in the effective coefficients of the network with an increase in the number of elements. Even with the coarsest discretization considered, i.e. , we obtain a response which is within of the response obtained at .
In Fig. 51, an anisotropic network with a higher coverage is considered (). For this case, the bonded area in the network is larger, resulting in a comparatively higher hygroscopic strain and overall deformation. As in the previous case, the anisotropic orientation of fibres in the network results in a higher expansion along the cross direction. Hence, both the anisotropy and higher coverage contribute to a higher expansion in the cross direction when compared with networks having low coverages.
Finally, an isotropic network () with low coverage of is considered. The values of the expansivity in machine direction () and cross direction () are listed in Tab. 3. Theoretically, if the network would be truly isotropic and representative, the listed values for the machine and cross directions should be identical (i.e. ) and the shear component should vanish (i.e. ). However, and differ by approximately and due to the sparsity and the statistically small size of the network used. The LSXFEM agrees within with the reference conforming FEM results, indicating an adequate approximation of the hygroscopic properties of the considered network.
Approach  
conf. FEM  
LSXFEM 
4.2.3 Local behaviour
Local strain distributions in the mediumcomplexity fibrous networks obtained by the LSXFEM formalism are analysed in this section by comparing them against the high resolution conforming FEM solutions. In Figs. (a)a and (a)a, the normalized strain distributions are plotted for, respectively, a conforming FEM mesh ( nodes) and the LSXFEM ( nodes). The strain distribution obtained by the LSXFEM is accurately computed in the overall network as well as in the bonds. For a better comparison, the magnified views of the normalized strain distribution are plotted in Figs. (b)b and (b)b.
The normalized strains along the crosssection A–A and B–B are further plotted in Fig. 60. At a few locations, e.g. close to in the dashed box in Fig. (a)a, a strain value is predicted by LSXFEM, whereas there is no value for the conforming FEM. This is because there is a void for the conforming FEM, whereas in the LSXFEM solution, the corresponding finite element is identified as an interfacial element, as is noticeable in Fig. (b)b. Although there is a strain value predicted inside the entire interfacial element, the area of fibre in that element is captured accurately through the subtriangulation stated earlier. Similar effects would occur also in bonded areas.
4.3 Complex realistic network
In the last example, the capabilities of the proposed LSXFEM methodology are demonstrated on the realistic fibre network of Fig. 6, for which we were unable to obtain a conforming triangulation within a reasonable computing time due to its complex geometry. Several thousands of densely intersecting and overlapping fibres are present within the employed periodic cell, which proved to pose a challenge even for stateoftheart mesh generators such as Gmsh. The individual fibres are of length and width , and the material parameters used for the anisotropic behaviour of fibres, their hygroscopic expansion, as well as the integration tolerance of the LSXFEM method are identical to those of the previous examples. A unit change in moisture content, , is adopted, which is again assumed uniform over the entire unit cell. Because of the high density of the network, a uniform mesh with nodes is used.
The resulting effective hygroexpansive coefficients are summarized in Tab. 4. The coefficients significantly more accurately satisfy the conditions imposed by isotropy (i.e. and ), as compared to Tab. 3. This is expected because the considered cell contains many more fibres, sampling their orientation distribution much more closely. The remaining discrepancies are due to the relatively small periodic cell size used, of twice the fibre length. Local distributions of strains obtained by the LSXFEM are shown in Fig. 66, reflecting a much more complicated spatial dependence as compared to the mediumcomplexity network of Fig. (a)a.
The presented results demonstrate the ability of the proposed LSXFEM formalism to resolve the effect of the complexity of such a dense network on the hygromechanical behaviour at the sheetscale. With the LSXFEM approach, it is thus possible to systematically study the effect of different microstructural parameters of realistic networks with high coverages at an affordable computational cost.
Approach  
LSXFEM 
5 Conclusions
This paper presented a computational methodology for analysis of the hygromechanical response of fibre networks using an XFEM framework with a level set based geometry description. A twodimensional unit cell consisting of a network of fibres has been considered to this end. It was subjected to free expansion by a uniform field of hygroscopic loading caused by a change in moisture. Rectangular fibres were generated randomly to produce periodic cells representing complex networks. The initial geometry of each fibre was described using a levelset function. The coupling between levelsets and the XFEM was carried out by means of a Heaviside function enrichment based on the levelset function. The edges of the fibres were captured from the nodal values of the levelset function. In this manner, a connection was established between the finite element mesh and the internal fibre geometry using the levelset formalism, which simplifies and improves the efficiency of computations for complex geometries in XFEM.
An accurate tracking of the boundaries of the fibres in the network was achieved, especially in the bonded regions, which matter for the prediction of the hygromechanical response of the network subjected to moisture infiltration. Local fields in the network (stress, strains and displacements) could be evaluated accurately as compared to a reference finite element solution with a fine conforming mesh. The LSXFEM results were obtained with a relatively small system size, whereas a comparably accurate nonconforming FEM simulation required a significantly higher number of elements and degreesoffreedom. It has been demonstrated, on a realistic network, that the LSXFEM approach is capable of predicting hygroexpansive behaviour of highcoverage networks, for which we were unable to obtain a conforming triangulation within a reasonable amount of computing time. The LSXFEM formalism is thus a suitable tool for modelling of realistic and highcoverage fibrous networks while capturing the interfaces with an adequate accuracy.
Possible extensions of this work include the modelling of irreversible shrinkage behaviour in paperlike materials. During the manufacturing process of paper, internal stresses are developed when the paper is dried under tension. They are released when the paper is subjected to a moisture cycle, resulting in irreversible shrinkage (Bosco et al., 2015a). This behaviour can be included by considering suitable constitutive models.
Acknowledgements
The first author would like to acknowledge the financial support granted by the European Commision, grant agreement № 20130043, as well as Materials Innovation Institute (M2i) and Canon Production Printing.
References
 Akker (1962) Akker, J. A. V., 1962. Some theoretical considerations on the mechanical properties of fibrous structures. The Formation and Structure of Paper. Technical Section of British Papers and Board Makers Association 1, 205–241.
 Bosco et al. (2015a) Bosco, E., Peerlings, R. H. J., Geers, M. G. D., 2015a. Explaining irreversible hygroscopic strains in paper: a multiscale modelling study on the role of fibre activation and microcompressions. Mechanics of Materials 91, 76–94.
 Bosco et al. (2015b) Bosco, E., Peerlings, R. H. J., Geers, M. G. D., 2015b. Predicting hygroelastic properties of paper sheets based on an idealized model of the underlying fibrous network. International Journal of Solids and Structures 5657, 43–52.
 Bosco et al. (2017) Bosco, E., Peerlings, R. H. J., Geers, M. G. D., 2017. Asymptotic homogenization of hygrothermomechanical properties of fibrous networks. International Journal of Solids and Structures 115116, 180–189.
 Boutin (1996) Boutin, C., 1996. Microstructural effects in elastic composites. International Journal of Solids and Structures 33, 1023–1051.
 Bronkhorst (2003) Bronkhorst, C. A., 2003. Modeling paper as a twodimensional elasticplastic stochastic network. International Journal of Solids and Structures 40, 5441–5454.
 Cox (1952) Cox, H., 1952. The elasticity and strength of paper and other fibrous materials. British Journal of Applied Physics 3, 72–79.
 Daux et al. (2000) Daux, C., Moes, N., Dolbow, J., Sukumar, N., Belytschko, T., 2000. Arbitrary branched and intersecting cracks with the extended finite element method. International Journal for Numerical Methods in Engineering 48, 1741–1760.
 Dirrenberger et al. (2014) Dirrenberger, J., Forest, S., Jeulin, D., 2014. Towards gigantic RVE sizes for 3D stochastic fibrous networks. International Journal of Solids and Structures 51, 359–376.
 Erkkila et al. (2015) Erkkila, A., Leppanen, T., Ora, M., Tuovinen, T., Puurtinen, A., 2015. Hygroexpansivity of anisotropic sheets. Composites Science and Technology 30, 325–334.
 Geuzaine and Remacle (2009) Geuzaine, C., Remacle, J.F., sep 2009. Gmsh: A 3D finite element mesh generator with builtin pre and postprocessing facilities. International Journal for Numerical Methods in Engineering 79 (11), 1309–1331.
 Guedes and Kikuchi (1990) Guedes, J., Kikuchi, N., 1990. Preprocessing and postprocessing for materials based on the homogenization method with adaptive finite element methods. Computer Methods in Applied Mechanics and Engineering 83, 143–198.
 Kulachenko and Uesaka (2012) Kulachenko, A., Uesaka, T., 2012. Direct simulations of fiber network deformation and failure. Mechanics of Materials 51, 1–14.
 Larsson and Wagberg (2008) Larsson, P. A., Wagberg, L., 2008. Influence of fibrefibre joint properties on the dimensional stability of paper. Cellulose 15, 515–525.
 Niskanen (1998) Niskanen, K., 1998. Paper Physics (Papermaking science and technology)., 2nd Edition. Fapet Oy.
 Niskanen et al. (1997) Niskanen, K., Kuskowski, S. J., Bronkhorst, C. A., 1997. Dynamic hygroexpansion of paperboards. Nordic Pulp Paper Research Journal 12, 103–110.
 Nordman (1958) Nordman, L., 1958. Laboratory investigations into the dimensional stability of paper. Technical Association of the Pulp and Paper Industry 41, 23–30.
 Osher and Fedkiw (2003) Osher, S., Fedkiw, R., 2003. Level set methods and dynamic implicit surfaces, 1st Edition. SpringerVerlag New York.
 Page and Tydemann (1961) Page, D. H., Tydemann, P. A., September 1961. A new theory of the shrinkage, structure and properties of paper. Transactions of the symposium held at Oxford, ed. F. Bolam. London: Technical Section, British Paper and Board Makers’ Association, 397–421.
 Peerlings and Fleck (2004) Peerlings, R. H. J., Fleck, N. A., 2004. Computational evaluation of strain gradient elasticity constants. International Journal for Multiscale Computational Engineering 2, 599–619.
 Ramasubramanian and Perkins (1988) Ramasubramanian, M. K., Perkins, R. W., 1988. Computer simulation of the uniaxial elasticplastic behavior of paper. ASME Journal of Engineering Materials and Technology 110, 117–123.
 Rivara (1997) Rivara, M. C., 1997. New longestedge algorithms for the refinement and/or improvement of unstructured triangulations. International Journal for Numerical Methods in Engineering 40, 3313–3324.
 Roylance (1996) Roylance, D., 1996. Mechanics of Materials., 1st Edition. Wiley and Sons.
 Sabit (1997) Sabit, A., 1997. Paper machine clothing: Key to paper making process. CRC Press 1st edition.
 Salmen and Olsson (2014) Salmen, L., Olsson, A. M., 2014. Mechanosorptive creep in pulp fibres and paper. Wood Science and Technology 48, 569–580.
 Schulgasser and Page (1988) Schulgasser, K., Page, D. H., 1988. The influence of transverse fibre properties on the inplane elastic behaviour of paper. Composite Science and Technology 32, 279–292.
 Sethian (1999) Sethian, J. A., 1999. Level set methods and fast marching methods, 2nd Edition. Cambridge University Press.
 Shahsavari and Picu (2013) Shahsavari, A., Picu, R., 2013. Size effect on mechanical behavior of random fiber networks. International Journal of Solids and Structures 50, 3332–3338.
 Starzewski and Stahl (2000) Starzewski, O. M., Stahl, D. C., 2000. Random fiber networks and special elastic orthotropy of paper. Journal of Elasticity 60, 131–149.
 Sukumar et al. (2001) Sukumar, N., Chopp, D. L., Moes, N., Belytschko, T., 2001. Modeling holes and inclusions by level sets in the extended finiteelement method. Computer Methods in Applied Mechanics and Engineering 190, 6183–6200.
 Uesaka (1994) Uesaka, T., 1994. General formula for hygroexansion of paper. Journal of Material Science 29, 2372–2377.
 Uesaka and Qi (1994) Uesaka, T., Qi, D., 1994. Hygroexpansivity of papereffects of fibre to fibre bonding. Journal of Pulp and Paper Science 20, 175–179.
Comments
There are no comments yet.