Level set based eXtended finite element modelling of the response of fibrous networks under hygroscopic swelling

06/23/2020 ∙ by P. Samantray, et al. ∙ TU Eindhoven 0

Materials like paper, consisting of a network of natural fibres, exposed to variations in moisture, undergo changes in geometrical and mechanical properties. This behaviour is particularly important for understanding the hygro-mechanical response of sheets of paper in applications like digital printing. A two-dimensional microstructural model of a fibrous network is therefore developed to upscale the hygro-expansion of individual fibres, through their interaction, to the resulting overall expansion of the network. The fibres are modelled with rectangular shapes and are assumed to be perfectly bonded where they overlap. For realistic networks the number of bonds is large and the network is geometrically so complex that discretizing it by conventional, geometry-conforming, finite elements is cumbersome. The combination of a level-set and XFEM formalism enables the use of regular, structured grids in order to model the complex microstructural geometry. In this approach, the fibres are described implicitly by a level-set function. In order to represent the fibre boundaries in the fibrous network, an XFEM discretization is used together with a Heaviside enrichment function. Numerical results demonstrate that the proposed approach successfully captures the hygro-expansive properties of the network with fewer degrees of freedom compared to classical FEM, preserving desired accuracy.



There are no comments yet.


page 3

page 4

page 14

page 15

page 16

page 19

page 21

page 24

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

At the micro-scale, 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, moisture-induced swelling takes place, which is called hygro-expansion (Larsson and Wagberg, 2008). In practical applications related to printing, this results in macro-scale 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 sheet-scale 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 hygro-expansive 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 elasto-plastic 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). Inter-fibre 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 two-dimensional 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.

(a) micrograph of the surface of a printing paper sheet
(b) an idealized bond between two fibres
Figure 3: Illustration of the microstructure of a paper sheet and the role of fibre bonds in hygro-expansion.

Several other studies addressed the hygro-expansivity and dimensional stability of paper. Some studies focused on the relation between hygro-expansion 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 macro-scale 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 hygro-expansion of paper fibres. However, most of these studies lack in modelling the hygro-mechanical behaviour of complex fibrous networks like paper subjected to moisture fields.

The aim of the present study is to model the hygro-mechanical behaviour of fibrous materials through a multi-scale analysis of the network using periodic homogenization (Guedes and Kikuchi, 1990; Boutin, 1996; Peerlings and Fleck, 2004). The fibres are modelled as two-dimensional ribbon-like elements in a network subjected to hygro-expansion, 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 non-conforming finite elements does not allow to accurately capture the geometry of the fibres and of the voids and bonded regions between them. Non-conforming 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 ribbon-shaped 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 non-conforming mesh and which is less computationally demanding than a conforming discretization strategy.

(a) complex fibre network
(b) magnified view
Figure 6: A model of a typical fibrous network represented as a periodic unit cell, corresponding to coverage  and the isotropic case (). Each rectangular outline represents a single fibre, where the colour is used to distinguish individual fibres more clearly.

Here, an advanced discretization scheme is combined with the hygro-mechanical model developed by (Bosco et al., 2017) to enable modelling of larger systems in two-dimensional configurations, with the potential to be extended towards three dimensions. The level-set 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 level-set methodology is a mathematical tool to describe complex and time-evolving boundaries (or, more generally, geometry) implicitly (Sethian, 1999). The boundary of a fibre is represented by the zero level-set of a higher dimensional function. This results in a versatile geometrical description that is de-coupled 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 level-set 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 level-set 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 hygro-mechanical 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 medium-complexity 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 fibre-level constitutive model formulated in (Bosco et al., 2015b) is adopted here. A two-dimensional 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.

Figure 7: The local and global coordinate axes.

The hygro-mechanical 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


where the hygroscopic strain  is given by


In these expressions, , , and  are the elastic constitutive tensor, the strain tensor and the hygro-expansivity tensor of the fibre respectively. In matrix notation, and  are represented as


In Eq. (3), and  denote the elastic moduli in the longitudinal and transverse direction with respect to the fibre axis, is the in-plane shear modulus, and  and  in-plane 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 hygro-expansion and elasticity of the fibres. This implies full displacement strain compatibility inside a bond for all the fibres interconnected in that bond.

2.2 Level-set formalism

The level-set 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 level-set of a higher dimensional level-set function  (see Fig. 8).

Figure 8: Sign convention used for signed distance function representing the rectangular fibres.

In most cases, including the current work, the level-set function gives the signed distance of a point  to the interface of the rectangular fibre denoted by  in Fig. 8. The level-set function  is defined as


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.

(a) fibre orientation distribution function
(b) periodic unit cell
Figure 11: Anisotropic orientation function and a random network of fibres.

The orientations of the randomly generated fibres in the network satisfy a probability density function based on 

(Cox, 1952):


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 level-set function, , , a periodicity condition is enforced on the unit cell. The level-set function that allows identifying the fibres and voids for the entire network can be described as


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 geometry-conforming 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 non-conforming 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 level-set function  is adopted (Sukumar et al., 2001) as an enrichment function, since the geometry is already described by the level-set function . The displacement interpolation then reads (Daux et al., 2000, Section 4)


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 level-set (to capture the geometry), referred to as LS-XFEM in what follows, is discussed in detail in the next section. A LS-XFEM formalism that makes use of the level-set information in this particular XFEM setting thus greatly simplifies the set up of the model under consideration.

3.2 Discretization with the level-set and XFEM formalism (LS-XFEM)

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


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


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 strain-displacement relationship for each finite element, , where  is the standard strain-displacement 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


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.


or, using Eq. (10),


This equation is of the form of a linear system


in which  collects all nodal displacements, and stiffness matrix and the hygroscopic load vector are given by


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),


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 level-set 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 sub-divided 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:

  1. 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.

  2. The subdivision continues for triangles intersected by the fibre’s boundaries Figs. (c)c(d)d If a triangle lies completely inside or outside the fibre, the subdivision is not performed.

  3. The sub-triangulation is terminated when a specified tolerance limit in terms of a minimum triangle area is achieved.

  4. 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 sub-triangulation 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 sub-triangulation 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.

Figure 16: Partitioning of a triangular finite element intersected by the boundaries of two (rectangular fibres) for the area integration. (a) The finite element  is partially covered by two fibres. (b) As a first step, fibre  is considered. Since the finite element lies partially in the fibre (the level-set values at the vertices are not all of the same sign), it is bisected along the longest edge, resulting in triangles  and . (c) Both of the newly formed triangles  and  are again lying partially in fibre . They are hence further bisected. (d) Triangles , , and  are further split; remains unaffected as it is completely outside the fibre. The partitioning further continues until a certain tolerance is reached.

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.

Figure 22: Illustration of the longest edge propagation path (LEPP) refinement in a triangular mesh. (a) The LEPP is given by  with  being the terminal triangle. (b)  is bisected along its longest edge and the new LEPP is , (c)  and , being the terminal triangles, are bisected along their longest edge. The LEPP is . (d)  and  are bisected along their longest edge. The LEPP is  and . (e)  and  are bisected. This concludes the refinement.
(a) refined mesh
(b) magnified view
Figure 25: (a) A refined mesh. (b) Magnified view depicting the non-conforming character of the mesh relative to fibre interfaces.

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 LS-XFEM 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 LS-XFEM formalism for fibrous networks, simplified illustrations are first presented below. Subsequently, the formalism is used for medium-complexity and complex realistic networks to study the effect of microstructural features on the sheet-scale 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 LS-XFEM formalism as compared to the non-conforming FEM. The main reason for using a non-conforming 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 non-conforming fixed grid by the standard FEM and the LS-XFEM 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 non-conforming FEM. In the LS-XFEM 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.

(a) infinitely long fibres and loading considered
(b) discretized unit cell with elements non-conforming to the geometry
Figure 28: (a) An infinite number of parallel, equispaced, infinitely long fibres subjected to a uniaxial stress . (b) Corresponding periodic unit cell with a regular mesh non-conforming to the geometry.
Table 1: Computed stresses and areas obtained by the standard non-conforming FEM and LS-XFEM.

As expected, the non-conforming 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 LS-XFEM formalism yields an accurate stress distribution with the same mesh. As the mesh is refined, the non-conforming FEM still struggles to yield the expected accuracy, showing a slow convergence rate. However, the LS-XFEM 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 LS-XFEM 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 LS-XFEM method.

(a) infinitely long fibre network and loading considered
(b) discretized unit cell with  elements non-conforming to the geometry
Figure 31: (a) An infinite number of two orthogonal families of parallel, equispaced, infinitely long fibres subjected to a uniaxial stress . (b) Corresponding periodic unit cell with a regular mesh non-conforming to the geometry. The  profiles along the cross-sections A–A and B–B are depicted in Fig. 40.
(a) reference conforming solution ()
(b) LS-XFEM ()
Figure 34: Stress distributions, , in the horizontal fibre.
(a) reference conforming solution ()
(b) LS-XFEM ()
Figure 37: Stress distributions, , in the vertical fibre.
(a) horizontal fibre, cross-section A–A
(b) vertical fibre, cross-section B–B
Figure 40: Stress distributions, , in cross-sections A–A and B–B as defined in Fig. (b)b for conforming FEM and LS-XFEM.

Figs. 34 and 37 show the normalized stress distribution  in the bonded regions for the horizontal and vertical fibres as computed by the LS-XFEM with a non-conforming mesh and the conforming FEM reference solution. The coarse triangular mesh used for the LS-XFEM 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 LS-XFEM differs from the reference solution with a deviation of 20, as observed in Fig. (a)a. This is because with the LS-XFEM formalism employing a geometry non-conforming 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 LS-XFEM. However, in the geometry-conforming 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 LS-XFEM formalism. Most importantly, as noticed in these plots, the normalized stress distributions predicted by the LS-XFEM formalism are qualitatively and quantitatively similar to the reference solutions at other regions, which is further ascertained in the cross-sectional 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 LS-XFEM 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 Medium-complexity 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 LS-XFEM 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 hygro-mechanical 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 LS-XFEM 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 LS-XFEM (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.

(a) initial network of coverage  and 
(b) refined mesh for LS-XFEM
(c) conforming mesh for FEM
Figure 44: Geometry of the complex network and the magnified view of meshes, coverage  and anisotropy parameter .
(a) reference conforming FEM with  ( nodes)
(b) LS-XFEM with refined mesh ( nodes)
Figure 47: Deformed network as computed by the conforming FEM and LS-XFEM. The displacements have been magnified by a factor of .
conf. FEM
Table 2: Computed effective hygro-expansive coefficients normalized with  for the anisotropic network ( and ).

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 LS-XFEM 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 hygro-expansivity of the network with an increasing number of elements in the mesh for the LS-XFEM 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 LS-XFEM agrees within  with the reference conforming FEM results, indicating an adequate approximation of the hygroscopic properties of the considered network.

Figure 48: Convergence of the effective hygroscopic coefficients normalized with  obtained by the LS-XFEM as a function of the characteristic element size.
conf. FEM
Table 3: Computed effective hygro-expansive coefficients normalized with  for the isotropic network ( and ).
(a) initial network
(b) deformed network, LS-XFEM
Figure 51: A high coverage network ( and ). The displacements have been magnified by a factor of .

4.2.3 Local behaviour

Local strain distributions in the medium-complexity fibrous networks obtained by the LS-XFEM 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 LS-XFEM ( nodes). The strain distribution obtained by the LS-XFEM 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.

(a) normalized strain distribution,
(b) magnified view of the dashed region
Figure 54: Reference conforming FEM solution for medium-complexity network with mesh size  ( nodes, , and ).
(a) normalized strain distribution,
(b) magnified view of the dashed region
Figure 57: LS-XFEM solution for medium-complexity network obtained with refined mesh ( nodes, , and ).
cross-section A–A
(a) cross-section A–A
(b) cross-section B–B
Figure 60: Normalized strain in the network at cross-sections A–A and B–B.
(a) conforming FEM ( nodes)
(b) LS-XFEM formalism ( nodes)
Figure 63: A magnified view of the discretized network at cross-section A–A. An interfacial LX-XFEM element located close to  is highlighted.

The normalized strains along the cross-section 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 LS-XFEM, whereas there is no value for the conforming FEM. This is because there is a void for the conforming FEM, whereas in the LS-XFEM 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 sub-triangulation stated earlier. Similar effects would occur also in bonded areas.

4.3 Complex realistic network

In the last example, the capabilities of the proposed LS-XFEM 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 state-of-the-art 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 LS-XFEM 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 hygro-expansive 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 LS-XFEM are shown in Fig. 66, reflecting a much more complicated spatial dependence as compared to the medium-complexity network of Fig. (a)a.

The presented results demonstrate the ability of the proposed LS-XFEM formalism to resolve the effect of the complexity of such a dense network on the hygro-mechanical behaviour at the sheet-scale. With the LS-XFEM 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.

Table 4: Computed effective hygro-expansive coefficients normalized with  for the realistic isotropic network shown in Fig. 6 ( and ).
(a) normalized strain distribution,
(b) normalized strain distribution,
Figure 66: LS-XFEM solution for the realistic network shown in Fig. 6, obtained with a regular mesh ( nodes, , and ).

5 Conclusions

This paper presented a computational methodology for analysis of the hygro-mechanical response of fibre networks using an XFEM framework with a level set based geometry description. A two-dimensional 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 level-set function. The coupling between level-sets and the XFEM was carried out by means of a Heaviside function enrichment based on the level-set function. The edges of the fibres were captured from the nodal values of the level-set function. In this manner, a connection was established between the finite element mesh and the internal fibre geometry using the level-set 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 hygro-mechanical 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 LS-XFEM results were obtained with a relatively small system size, whereas a comparably accurate non-conforming FEM simulation required a significantly higher number of elements and degrees-of-freedom. It has been demonstrated, on a realistic network, that the LS-XFEM approach is capable of predicting hygro-expansive behaviour of high-coverage networks, for which we were unable to obtain a conforming triangulation within a reasonable amount of computing time. The LS-XFEM formalism is thus a suitable tool for modelling of realistic and high-coverage fibrous networks while capturing the interfaces with an adequate accuracy.

Possible extensions of this work include the modelling of irreversible shrinkage behaviour in paper-like 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.


The first author would like to acknowledge the financial support granted by the European Commision, grant agreement № 2013-0043, as well as Materials Innovation Institute (M2i) and Canon Production Printing.


  • 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 multi-scale modelling study on the role of fibre activation and micro-compressions. Mechanics of Materials 91, 76–94.
  • Bosco et al. (2015b) Bosco, E., Peerlings, R. H. J., Geers, M. G. D., 2015b. Predicting hygro-elastic properties of paper sheets based on an idealized model of the underlying fibrous network. International Journal of Solids and Structures 56-57, 43–52.
  • Bosco et al. (2017) Bosco, E., Peerlings, R. H. J., Geers, M. G. D., 2017. Asymptotic homogenization of hygro-thermo-mechanical properties of fibrous networks. International Journal of Solids and Structures 115-116, 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 two-dimensional elastic-plastic 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 3-D finite element mesh generator with built-in pre- and post-processing 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 fibre-fibre 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 hygro-expansion 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. Springer-Verlag 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 elastic-plastic behavior of paper. ASME Journal of Engineering Materials and Technology 110, 117–123.
  • Rivara (1997) Rivara, M. C., 1997. New longest-edge 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 in-plane 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 finite-element 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 paper-effects of fibre to fibre bonding. Journal of Pulp and Paper Science 20, 175–179.