Micro-mechanical modeling of irreversible hygroscopic strain in paper sheets exposed to moisture cycles International Journal of Solids and Structures

Paper is a complex material consisting of a network of cellulose ﬁbres at the micro-level. During manu- facturing, the network is dried under restraint conditions due to tension in the paper web in machine direction. This gives rise to internal strains that are stored in the produced sheet. Upon exposure to a moisture cycle, these strains may be released. This results in permanent shrinkage that may cause insta-bilities such as curl or waviness of the sheet. The prime objective of this paper is to model this irreversible shrinkage and to link its magnitude to the properties of the ﬁbres and of the network. For this purpose, randomly generated ﬁbrous networks of different coverages (i.e. ratio of the area occupied by ﬁbres and that of the sheet) are modeled by means of a periodic representative volume element (RVE). Within such RVEs, a ﬁnite element method combined with a kinematic hardening plasticity model at the scale of the ﬁbres is used to capture the irreversible response. The computational results obtained demonstrate that the magnitude of the irreversible strains increases with coverage until a certain coverage and beyond that coverage decreases in magnitude. This phenomenon is explained by considering the area fraction of free-standing ﬁbre segments relative to bonded ﬁbre segments in the network. A structure–property depen- dency of irreversible strains at the sheet-level on the micro-structural parameters of the network is thereby established. (cid:1) 2021 The Authors. Published by Elsevier Ltd. ThisisanopenaccessarticleundertheCCBYlicense(http:// creativecommons.org/licenses/by/4.0/).


Introduction
Paper belongs to a category of materials that are composed of a network of fibres produced from wood pulp. At the micro-level, parts of the fibres are bonded with each other at some locations while they are free-standing elsewhere. During the manufacturing process, the paper fibres obtain a preferential orientation (in machine direction) which results in an anisotropic behavior of the paper. Upon exposure of a sheet of paper to a humid environment, or to liquid water, anisotropic moisture induced swelling, or hygro-expansion, takes place (Larsson and Wågberg, 2008). This swelling originates at the scale of fibres and gets propagated through the network by the inter-fibre bonds, where interactions between the hygroscopic and the mechanical properties of the fibres occur -see Fig. 1a. In applications such as printing, the consequences of these deformations may be observed at the macroscale through curl, waviness and cockling. Understanding these phenomena and their dependence on the properties of the fibres and the network is essential to predict the overall response of the material at the sheet-level.
The hygroscopic response of paper is dependent on the drying phase during its manufacturing process. In the case of restrained drying, an initial strain is generated in the paper fibres due to the fibre segment activation mechanism (Vainio and Paulapuro, 2007). According to this mechanism, restrained drying of the fibres leads to a decrease in the orientation angle of the micro-fibrils with respect to the axis of the fibres, as sketched in Fig. 1b. As a result, the fibres are stretched and develop an internal stress. This mechanism induces dried-in strains that remain in the paper after manufacturing. When the manufactured paper is exposed to a moisture cycle beyond a certain moisture content, the dried-in strains are released to some extent, resulting in irreversible shrinkage. This effect was measured for instance by Mäkelä (2010) and is illustrated schematically in Fig. 2b. It is generally interpreted as a permanent or irreversible hygroscopic strain. However, this dried-in strain is absent if the paper is dried freely in the manufacturing process. Therefore, a reversible hygroscopic response is observed upon exposure to a moisture cycle, as sketched in Fig. 2a. In order to fully understand this characteristic behavior of paper and assess https://doi.org/10.1016/j.ijsolstr.2021.03.011 0020-7683/Ó 2021 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). the micro-structural factors affecting it, a model of the underlying complex fibrous network is essential.
Published works on the relationship between the macroscopic behavior of fibrous networks and micro-structural parameters address the dependence of stress transfer on the correlation length in a network (Åstrom et al., 1994); fracture of the network dominated by fibres longer than the average length (Alava and Ritala, 1990); the mechanical response of fibrous networks depending on the curl ratio (Yi et al., 2004); the magnetic and mechanical response of bonded networks of metal fibres due to filling in voids (Clyne et al., 2005); the influence of the aspect ratio and fibre concentration on the effective stiffness of planar fibre networks (Wu and Dzenis, 2005). Early attempts at structure-property relations for fibrous networks specifically addressed the mechanical response. Cox (1952), van den Akker (1962) and Bronkhorst (2003) studied the stiffness of fibre networks by taking into account the axial forces and transverse properties of fibres. Some of the papers considered the bonds in the networks enabling models to transmit axial and shear forces along with torsional and bending moments (Schulgasser and Page, 1988;Starzewski and Stahl, 2000;Ramasubramanian and Perkins, 1988). Some recently developed numerical models for fibrous networks describe only the mechanical response of a random fibre network, based on strong micro-structural simplifications (Kulachenko and Uesaka, 2012;Shahsavari and Picu, 2013;Dirrenberger et al., 2014).
Among the available studies, some focused on the hygroexpansion and dimensional stability of the paper (Uesaka and Qi, 1994;Nordman, 1958;Salmén et al., 1987;Uesaka, 1994;Kulachenko and Motamedian, 2019). In them, the hygroexpansivity of networks is studied by establishing a relation between the dimensional stability and different parameters that affect it, like the fibre-fibre contact ratio, relative humidity, wet The reversible part of the hygro-expansion is believed to stem mainly from fibreto-fibre bonds. Fibres expand predominantly in lateral direction. The free-standing parts of a fibre can do this freely, without affecting the sheet-scale dimensional stability. In bonded areas, however, the lateral expansion of one fibre is translated into longitudinal elongation of the other fibre -which does affect the sheet-scale dimensions. (b) Due to the fibre activation mechanism, micro-fibrils in the fibre wall are more aligned mutually and with the applied load or restraint in the post-drying state (Y) compared with their relaxed state (X).

Fig. 2.
Sketch of the typical hygro-mechanical response observed in experiments on hand sheets dried under different degrees of restraint. (a) Freely dried sheets show a more or less linear dependence of hygroscopic strain on moisture content, which is furthermore largely reversible. (b) Restrained-dried sheets exhibit a non-linear, partially irreversible response which results in a net shrinkage of the sheet in the first wetting-drying cycle. Subsequent cycles are reversible, but with a lower expansivity than that of freely dried sheets.
pressing of sheets, density or exchange of heat and moisture between paper and surroundings. Bosco et al. (2015) specifically studied the micro-structural origin of irreversible shrinkage following restrained drying, on the basis of an idealised meso-scale network model. The underlying mechanism was postulated to be fibre activation (Fig. 1b). It was shown to be active mainly in the free-standing parts of the fibres, as the interaction between differently oriented fibres in a bond largely restricts the mechanism. However, these findings were based on a highly idealised, periodic network model, in which the fibres are infinitely long fibres and aligned perfectly with either the machine direction or cross direction.
The main objective of the present work is to study the influence of a more realistic, random network geometry on the predicted irreversible sheet-level hygro-mechanical response of a fibre network which was previously subjected to restrained drying conditions. For this purpose, a two-dimensional periodic random network model similar to that of Bosco et al. (2017) is used, exploiting the kinematic hardening plasticity model of Bosco et al. (2015). One of the main features of the constitutive model is the occurrence of irreversible shrinkage in the longitudinal direction of fibres based on the fibre segment activation mechanism. This particular modeling aspect of irreversible shrinkage of complex paper fiber networks was not addressed so far.
Following our earlier work, the fibre network is generated in a two-dimensional framework with fibres being modeled as ribbon-like rectangular domains (Bosco et al., 2017). In the present contribution, an isotropic orientation distribution of the fibres is used in the generation of the network, so that the resulting networks are akin to hand sheets. A linear kinematic hardening model is used to represent the irreversible shrinkage in longitudinal direction of the fibres associated with the activation mechanism (Bosco et al., 2015). The phenomenon of manufacturing induced dried-in strains is thus explicitly incorporated. The release of strains upon moisture exposure is assumed to occur instantaneously (quasi-statically) so that a rate-independent model is adopted. This is an adequate assumption for understanding the effects of network properties on the sheet-scale behavior.
This paper is organized as follows. In Section 2, the modeling methodology for capturing the irreversible behavior of fibres is described, along with the details of the fibre constitutive model. Results highlighting the irreversible shrinkage of the complex network are presented in Section 3. Interactions between the activation mechanism and micro-structural features (e.g. coverage) are also investigated in this section. Finally, in Section 4, conclusions and perspectives are provided.
Throughout this contribution, the following notations are used. Scalars, vectors and (Cartesian) tensors are represented as a;ã and A respectively; 4 th order tensors are represented as 4 A. For tensor and vector operations, the following equivalent index notations are used, with Einstein's summation convention on repeated indices: A : B ¼ A ij B ji , with i; j ¼ x; y; z for the global reference system and i; j ¼ l; t; z for the local, fibre-reference system. The dyadic product between vectors is represented asãã. The Voigt notation used to represent tensors and tensor operations in a matrix format is as follows: a and A denote a column matrix and a matrix of scalars respectively. Matrix multiplication is denoted as A a ¼ A ij a j .

Random network creation
In the current study, assuming the fibres to be ribbon shaped in nature, a set of rectangular shaped fibres are randomly generated.
The length of the fibres is denoted l f and the width w f ¼ l f =50. The centroids of the fibres are randomly generated in the domain ðx; yÞ 2 ð0; lÞ Â ð0; lÞ within a periodic unit cell of edge length l ¼ l f as shown in Fig. 3. The orientations of the fibres are randomly generated as well. All fibres have an identical rectangular shape and dimensions. Parts of fibres which cross the boundary of the periodic cell are clipped and copied to the opposite boundary segment, thus ensuring periodicity of the geometry. In regions of overlap between two (or multiple) fibres, the fibres are assumed to be perfectly bonded. In order to characterise the denseness of the network, the coverage c is employed. It is defined as the ratio of the total area occupied by all fibres and the area of the periodic unit cell, l 2 .

Fibre longitudinal and transverse constitutive behavior
In order to describe the irreversibility phenomenon in the longitudinal direction of a fibre only due to tension, a plane stress plasticity model with linear kinematic hardening is adopted. After the paper is manufactured under restrained conditions, it is subjected to free expansion due to a moisture cycle, representing infiltration of ink during printing or absorption of moisture from the environment. Accordingly, there is no externally applied stress responsible for the irreversible phenomenon, as observed in the experiments described earlier (Mäkelä, 2010). However, the decreasing moisture content and the internally stored stress during the restrained drying process determine when plasticity may occur in the re-wetting phase via the release of dried-in strain, driven by the internal stress. A kinematic hardening plasticity model is used for modeling this storage and release of strain in the absence of an external stress, but triggered by the internal stress only. We will assume that this phenomenon takes place over a sufficiently long period of time so that time-dependent deformations can be assumed to have saturated.
In this contribution, a two-dimensional generalization of the fibre constitutive model proposed by Bosco et al. (2015) is adopted. The general constitutive law for a fibre makes use of a plane stress assumption in the thickness (z) direction. It furthermore assumes the material to be transversely isotropic, with the isotropy plane perpendicular to the fibre axis.
The total strain tensor e f in a fibre consists of a hygroscopic, an elastic and a plastic part: where the superscript f distinguishes these (fibre) strains from the macroscopic, sheet-scale strain. The hygroscopic and elastic parts of the strain are governed by the usual constitutive laws. The hygroscopic strain h e f is given by where b f is the hygroexpansivity tensor and Dv ¼ v À v 0 is the change of moisture content relative to the reference value v 0 . The moisture content v is defined here as the mass of the water in the paper divided by its total mass. The stress in the fibre is given by where 4 D f is the elasticity tensor.
In Voigt matrix notation, 4 D f and b f are represented with respect to the local fibre frame as where E l and E t are the elastic moduli in the longitudinal and transverse direction with respect to the fibre axis, G lt is the in-plane shear modulus, m lt and m tl are the in-plane Poisson's ratios. b l and b t denote the coefficients of hygroscopic expansion in the longitudinal and transverse direction of the fibre axis respectively. These material parameters are assumed to be independent of the moisture content v. This will give a better picture of the sole qualitative effect of the micro-structural parameters like coverage, free-standing fibres etc. on the irreversible strains, while keeping the model relatively simple and the number of parameters manageable. Quantitatively, however, this assumption may be of influence, as a significant dependence of e.g. the elastic stiffness of fibres on wetting conditions has been observed in experiments. Jentzen (1964) for instance reports the Young's modulus of a dry fibre to be 2.5 times that of the same fibre in a (fully) wet state. The local reference system ðl; t; zÞ of a fibre m is oriented at an angle h ðmÞ with respect to the global reference system ðx; y; zÞ, see Fig. 4. Therefore, it is important to transform the above expressions from the local frame of the fibre to the global frame for each fibre in the network before assembling the mechanical network -see e.g. Roylance (1996). The plastic part of the constitutive model represents the fibre activation mechanism. We postulate that this mechanism affects only the axial strain of the fibre and that its occurrence also depends only on the stress component along the fibre axis. Defin-ing the unit vectorp which points along the fibre axis, i.e. in the direction of the local l coordinate, the yield function is thus defined in terms of the projected longitudinal effective stress ðr f À q f Þ :pp (indicated in Fig. 4): where q f is the back stress and r y the (current) yield stress. Assuming associative plastic flow, the plastic strain tensor, p e f , and kinematic strain hardening tensor conjugate to the back stress, n f , evolve according to where the flow direction N is given by, according to (5), and the plastic multiplier _ c is governed by the Kuhn-Tucker The above definitions imply that the plastic strain tensor and the hardening variable have solely a normal component in axial direction, whereas all other components (in the local fibre frame) vanish, i.e. they may be written as Furthermore notice that as a result the plastic flow is necessarily non-isochoric.
The back stress q f appearing in (5) is given by the linear with the fourth-order tensor 4 H f represented in Voigt matrix notation as where H is the kinematic hardening modulus of the fibre.
Finally, the yield stress, r y , is postulated as a piece-wise linear function of the moisture content v, such that it decreases with an increase of the moisture content: where r y0 is the maximum yield stress, corresponding to a vanishing moisture content, and v w is the moisture content at which the yield stress vanishes. The longitudinal behavior in compression and the transverse behavior of the fibres both in tension and compression are assumed to follow a hygro-elastic constitutive behavior. This is justified by the fact that the irreversibility is attributed to tensile stresses in the longitudinal direction of the fibres according to the fibre segment activation mechanism explained before.

Fibre interaction modelling and finite element implementation
The hygro-elasto-plastic model developed in the previous section applies to each of the fibres of the random network according to their individual orientations. In the bonded regions, full kinematic compatibility is assumed between the fibres involved in the bond, i.e. the displacement vector is solely a function of the two-dimensional coordinate vectorx (plus time) and applies iden- tically for all fibres to which any given positionx belongs. Accordingly, the total fibre strain e f is equal in all fibres atx, and it equals the global, network strain e. The response of each fibre inx is computed individually and the resulting stresses, multiplied by the respective fibre thicknesses, enter the conventional twodimensional equilibrium equations. Complemented by periodicity conditions on the displacement, this results in a closed boundary value problem, which may be solved to yield the overall, sheetscale, response of the network -as well as all relevant local quantities. For details we refer the reader to the analogous discussion on the hygro-elastic case by Bosco et al. (2017).
The periodic unit cell containing the fibrous network is discretized by a regular finite element mesh and solved with a nonlinear incremental-iterative finite element method using Newton-Raphson iterations. In order to update the stresses in the linear kinematic hardening model, a return mapping algorithm based on a trial state is adopted. For the solution procedure, expressions for the plastic multiplier and the consistent tangent operators are required; these are derived in Appendix A. Further details of the nonlinear solution procedure are given in Appendix B.

Results and discussion
In this section, we subject the model as detailed above to the different phases of the manufacturing process of paper, followed by its exposure to a wetting-drying cycle. This entails the irreversible shrinkage of the fibrous network, which is investigated as a function of the micro-structural properties of the fibres and the overall properties of the network.

Modeling the effect of manufacturing constraints
During the manufacturing process of paper, strain may be 'dried in' due to web tension. We consider in particular the case of hand sheets which are dried under restraint. In the model, the stressfree, wet network is first dried under the constraint of zero overall (total) strain. At some point during the drying process the latter restraint is removed and the subsequent, remaining drying is hence free from restraint, i.e. the network (RVE) does not experience any overall stress. Internally, stresses will remain in the network, but these are internally equilibrated and do not require an externally applied stress. They are also accompanied by irreversible (plastic) strain which models the dried-in strain.
The resulting, 'as-manufactured' hand sheet is subsequently subjected to a wetting-drying cycle during which no constraint is applied. The increase in moisture content allows the dried-in strain to be released, as the plastic strain which was induced during manufacturing is removed under the influence of the internal stresses -which simultaneously are relaxed. Upon drying, the sheet is hence expected to have contracted irreversibly relative to the as-manufactured state.
The entire process of manufacturing followed by a wetting-drying cycle is represented schematically in Fig. 5. The diagram shows for each phase in the process how the moisture content v is varied and under which mechanical restraining conditions this is done. The moisture content evolution is applied uniformly in space. In the diagram, v low is the lowest moisture content reached, in the as-manufactured sheet as well as at the end of the moisture cycle, v high is the highest moisture content at the start of fibre bonding and also during the moisture cycle, v cons is the moisture content to which the network is dried under restraint, which is subsequently removed as the moisture content further drops. The moisture content range is chosen in accordance with the work performed by Larsson and Wågberg (2008), with v low ¼ 0:038; v cons ¼ 0:08; v high ¼ 0:12 and v w ¼ 0:13.

Fibre and network parameters considered
Following the work performed by Bosco et al. (2015) and adopting its values, the longitudinal elastic stiffness of the fibre is taken as E l ¼ 5:176 Â 10 2 r y0 , whereas in the transverse direction, (Strömbro and Gudmundson, 2008). The shear modulus is taken as G lt ¼ E l =10. The hygro-expansive coefficients along the longitudinal and transverse directions are given by b l ¼ 0:03 and b t ¼ 20b l (Niskanen, 1998). The in-plane Poisson ratios are m lt ¼ 0:3 and m tl ¼ 0:05 (Schulgasser and Page, 1988). The kinematic hardening modulus adopted in the plasticity model is H ¼ 100r y0 (Bosco et al., 2015). All of the above quantities are assumed to be constant irrespective of the moisture content. The numerical simulations have been performed with linear triangular elements, with a uniform element size h ¼ l=200.
The number of fibres n in a network is characterised by the coverage c, which is defined as the ratio of the total area occupied by all fibres and the area of the RVE. In what follows we consider coverages in the range 0:25 6 c 6 5:0. The coverage is c related to the grammage g via where q f is the mass density of the fibres and t is their thickness. For example, considering a density of the fibres q f ¼ 1500 kg=m 3 (Niskanen, 1998) and a fibre thickness t ¼ 10 lm, a coverage of c ¼ 1:0 as used in this section corresponds to a grammage of 15 g=m 2 . As the coverage is changed, with density and thickness remaining constant, the grammage scales linearly with the coverage as indicated by Eq. (13).

Simplified meso-scale network
Before we consider random networks, it is instructive to first study a periodic network consisting of two perpendicular families of equidistant, infinitely long fibres. This is the same network model as proposed by Bosco et al. (2015), but now solved numerically. This case illustrates the mechanisms at play in the clearest fashion and serves as a reference and verification for the complex networks considered below. Fig. 6 shows the plastic strain in fibre direction, p e f ll , normalised by b l Dv, where Dv ¼ v high À v low , before and after the wetting -i.e. before and after Phase 4 in Fig. 5. The result shown is for a network of coverage c ¼ 0:5. The before-wetting state, shown in Fig. 6a, corresponds with the as-manufactured sheet. It shows a significant amount of plastic strain, which would be visible at the macroscopic, 'sheet' level as dried-in strain. After subjecting the network to an increase in moisture content, these plastic strains are largely released, as evidenced in Fig. 6b. Note that little plastic strain develops in (and is subsequently released from) the bonded regions. Rather, most of the plastic activity takes place in the free-standing parts of the fibres. Fig. 7a shows the evolution of the plastic strain p e f ll in the freestanding fibre parts throughout the five phases of paper manufacturing and wetting-drying as defined in Fig. 5. It shows that the plastic strain is induced at the end of restrained drying. Subsequently, the restraint is released at a constant moisture content and then the network is dried freely. In both of these phases, the plastic strain remains constant as there is no restraint to produce the stress level necessary to cause plastic yielding -which increases during drying by virtue of the moisture dependent yield stress given by Eq. (12). As the network is subjected to a moisture increase again (in the moisture cycle -Phase 4), the yield stress decreases and the back stresses introduced during Phases 1-3 cause yielding in reverse direction, leading to the release of the plastic strains. Finally, when the moisture is decreased, in Phase 5, the yield stress increases again and there is no plastic activity anymore.
At the macro-scale, the local evolution of plastic strain as discussed above results in a global release of (dried-in) strain during the wetting-drying cycle, as depicted in Fig. 7b for three values of the coverage. In the initial phase of the moisture cycle, the overall strain increases linearly with moisture content. At this stage, the plastic strains present in the network remain constant and the response is hence hygro-elastic. The slope of this part of the response increases with coverage as a consequence of the fact that most of the hygro-elastic expansion is generated in the bonds and the area occupied by bonds is larger in networks with a higher coverage (Bosco et al., 2015). At some point of each of the curves, a kink is observed. At this point, the plastic strains start to get released until the moisture content starts to decrease. This occurs for the meso-scale networks of all coverages. Also, we can notice that as the coverage increases, the irreversible shrinkage also increases in magnitude. This computed response is in agreement with the findings by Bosco et al. (2015).

Complex networks
We proceed by applying the numerical modelling to more complex two-dimensional fibrous networks. Random isotropic fibrous networks have been generated for coverages c ¼ 0:5; 1:0; 2:0 and 5:0, and subsequently been subjected to the moisture and restraint history of Fig. 5. Fig. 8 shows the computed local plastic strains in the networks at the end of the manufacturing process and just before the wetting cycle, i.e. just before Phase 4. Upon subsequently subjecting these networks to an increase of the moisture level from v low to v high in Phase 4, the plastic strains have been partly released, as can be observed in the post-wetting plastic strain distributions in Fig. 9. At the macroscopic scale, this release in the networks is observed as an overall irreversible shrinkage  during the wetting-drying cycle (Phases 4 and 5), as illustrated in the hygro-expansive strain plots in Fig. 10. Fig. 8 reveals that the area of plastically deformed fibre segments increases as the coverage is increased from c ¼ 0:5 to c ¼ 2:0, followed by a decrease for the higher coverage, i.e. c ¼ 5:0. At the macro-level, the irreversible shrinkage also first slightly increases, from c ¼ 0:5 to c ¼ 1:0, and subsequently decreases, as noticed in Fig. 10. This figure also shows that the hygro-elastic expansivity monotonically increases with coverage, consistently with the observation made by Bosco et al. (2017) for hygro-elastic random networks -see also below.
A further observation which may be made in Fig. 10 is that the release of dried-in strain during wetting is no longer observed as a kink as in the idealised network ( Fig. 7b and (Bosco et al., 2015)), but as a smoother transition. This is despite the fact that a linear hardening law was used, which in the highly idealised network    Fig. 8 (with, left to right, coverage, c ¼ 0:5; 1:0; 2:0 and 5:0) after the wetting part of the applied moisture cycle. A significant part of the plastic strains which were present in the as-manufactured networks (Fig. 8) have been released -particularly from the free-standing fibre segments. resulted in a bi-linear hygroscopic network response, the response of the present random network is non-linear, as is also observed in experiments. This is due to the heterogeneity of the network, as a result of which the plastic strain release is initiated at different moisture contents in different fibre segments.
To establish that the irreversible strain observed at the macroscale (in Fig. 10) is stemming from the free-standing fibre parts, the hygro-mechanical response of networks with a low coverage, c ¼ 0:5 and c ¼ 1:0, has also been computed by allowing for plasticity only in the free-standing parts of the fibres -and conversely only in the bonded parts of fibres. The resulting macro-scale hygroexpansive strains, highlighting the irreversible strains, are shown in Fig. 11. Clearly, at the sheet-level, the irreversible shrinkage depends on plasticity occurring in the free-standing parts of fibres in the network, as the response obtained with plasticity only in the free-standing parts is virtually identical to that of the reference network. On the other hand, when only the bonded parts of the fibres in the network are allowed to deform plastically, the irreversible shrinkage at the sheet-level becomes negligible. This further confirms the finding that the observed irreversibility phenomenon is essentially due to plastic strains and internal stress stored in the free-standing fibre regions.
So far, we showed results for single networks with coverage c ¼ 0:5; 1:0; 2:0 and 5:0. To explore the variability between different realisations, the simulations have been repeated for five realisations at each coverage, still using randomly generated quasiisotropic networks. The irreversible shrinkage is defined as the difference of the overall strain, e xx , between the beginning and the end of the moisture cycle. The magnitude of the irreversible shrinkage, e irr , normalised by ðb l DvÞ, observed at the sheet-level for each of these networks and the area fraction distribution of free-standing fibre parts are depicted in Fig. 12. It is observed in Fig. 12b that as the coverage increases, the irreversible shrinkage increases until a coverage c ¼ 1:0, with a subsequent decrease until coverage c ¼ 5:0. This dependency is closely correlated with the area fraction of the free-standing part of fibres, in shown in Fig. 12a, which peaks at c ¼ 1:0. This confirms that the freestanding fibre parts in the network are responsible for the irreversible behavior of the networks upon exposure to moisture cycles.
Finally, as observed earlier, as the fraction of bonded regions increases with the coverage, the hygro-elastic expansivity of the networks increases with coverage. This is evidenced by the macroscopic coefficient of hygroscopic expansion extracted from the five realisations per coverage as shown in Fig. 13 -fully in agreement with the earlier predictions of Bosco et al. (2017).

Conclusions
In the present study, the dependency of the (experimentally observed) irreversible shrinkage at the sheet-level on the microstructural parameters in a fibrous network of paper is studied by a micro-structural modelling approach. The irreversible shrinkage  is postulated to be associated with micro-structural plastic strains caused by the restrained drying during the manufacturing process. These plastic strains get subsequently released upon wetting. For this study, ribbon shaped fibres were randomly generated. Using a standard homogenisation approach, the fibrous network which they constitute is represented as a periodic RVE. Relying on the hypothesis that the fibre segment activation mechanism is the physical cause of the permanent deformations in paper, it is assumed that the plastic deformation occurs in the longitudinal direction of the fibres. A constitutive model incorporating kine-matic hardening plasticity is adopted to model the plastic strain in the longitudinal direction of the fibres in the network. This allows modeling the irreversible deformation in paper subjected to moisture change with no applied external stresses, reflecting the real conditions to which the paper is subjected during manufacturing with restraint.
The key findings of the current work are: 1. The magnitude of the irreversible shrinkage at the macro-scale in the fibrous network model is a function of the coverage, i.e. density, of the network. Initially it increases with coverage until c ¼ 1:0, but subsequently decreases to small values at higher coverages, e.g. c ¼ 5:0. 2. The degree of irreversible strain release strongly correlates with the area fraction of the free-standing parts of fibres in the network. A higher coverage causes a lower density of free-standing fibres, which results in a lower irreversible shrinkage. 3. The non-linear response of the network at the macro-scale during the moisture cycle is captured even with a linear kinematic hardening model. It may be explained solely based on the heterogeneity of the random network.
One of the limitations of the present work is the assumption of full kinematical bonding of the fibres in the bonded regions. This may influence the observed decrease of irreversible shrinkage predicted by the model at higher coverages. A useful extension would therefore involve relaxing the kinematic constraints in the bonded regions. Furthermore, in the current study, the release of the driedin strains is represented as quasi-static phenomenon, i.e. any timedependent effects have been neglected. As a next step, a suitable rate-dependent model may be adopted to study the influence of time on the irreversible shrinkage.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. tangent operator     In order to predict the irreversible response of the twodimensional fibrous network, a finite element scheme is adopted using a non-linear incremental iterative methodology based on Newton-Raphson iterations. The RVE is discretized using a regular grid with constant strain triangular finite elements. The finite elements with their centroid lying inside the fibres contribute to the stiffness and hygroscopic load in the assembly of the global system of equations, while the other elements lie in voids. Now, as we are incorporating the hygro-expansion model in the non-linear equilibrium equations, it is important to obtain the correct initial estimate of the plastic multiplier Dc in the first iteration of every time increment to achieve convergence. We therefore derive his estimate for the first iteration of each time increment below.
The increment of plastic strain may be written as D p e f ¼ Dc N ¼ Dc N trial ðB:1Þ The stress at time t+Dt is written as  which is our desired prediction of stress at ðt þ DtÞ based on data at t.

Appendix C. Supplementary data
Supplementary data associated with this article can be found, in the online version, at https://doi.org/10.1016/j.ijsolstr.2021.03. 011.