Modeling and simulation of anisotropic cross-linked cellulose fiber networks with an out-of-plane topography

Non-woven cellulose fiber networks of low areal density are widely used in many industrial applications and consumer products. A discrete element method (DEM) modeling framework is advanced to simulate the formation of strongly anisotropic cellulose fiber network sheets in the dilute limit with simplified hydrodynamic and hydroelastic interactions. Our modeling accounts for in-plane fiber orientation and viscous drag indirectly by using theories developed by Niskanen (2018 Fundamentals of Papermaking, Trans. 9th Pulp and Paper Fundamental Research Symp. Cambridge, 1989 (FRC) pp 275–308) and Cox (1970 J. Fluid Mech. 44 791–810) respectively. Networks formed on a patterned and flat substrate are simulated for different fiber types, and their tensile response is used to assess the influence of the out-of-plane topographical pattern, specifically, on their stiffness and strength. Sheets with the same grammage and thickness, but composed with a higher fraction of softwood fiber (longer fibers with large diameter), have higher strength and higher strain to failure compared to sheets made from hardwood fibers (short fibers with small diameter). However, varying the fiber fraction produces only an insignificant variation in the initial sheet stiffness. The above simulation predictions are confirmed experimentally for sheets comprised of fibers with different ratios of Eucalyptus kraft and Northern Bleached Softwood Kraft fibers. Sheets with out-of-plane topography show an unsymmetric mass distribution, lower tensile stiffness, and lower tensile strength compared to those formed on a flat substrate. The additional fiber deformation modes activated by the out-of-plane topography, such as bending and twisting, explain these differences in the sheet mechanical characteristics.


Introduction
Fibrous materials are ubiquitous: several biological materials are composed of collagen fiber networks [1]; paper products are networks of cellulose fibers [2,3]; and plant primary cell walls have cellulose microfibrils acting as a load bearing component embedded in a hydrated matrix of hemicelluloses and pectins [4].Engineered materials like composites, gels, polymers and paper also exhibit fibrous structures.Fibrous network structures are classified in various ways.Based on the effect of thermal fluctuations, constituent fibers can be either 'athermal' or 'thermal'.The mechanical response of athermal fibers is not influenced by thermal fluctuations.Thermal fibers behave entropically and their mechanical behavior is affected by thermal fluctuations.The mechanical behavior of semi-flexible fibers lies between that of athermal and thermal fibers.Interfiber interactions can be classified as excluded volume interaction (no interfiber penetration), interaction through surface forces (e.g.adhesion in fiber bundles), and cross-linked interaction.Based on the fiber arrangement, networks can be woven or nonwoven.Woven networks are mostly artificial (e.g.textiles), while nonwovens can be either natural (e.g.collagen network) or engineered (e.g.paper).Finally, based on the constituent phases, networks can be classified as: networks without any matrix (e.g.fiberglass insulation), networks embedded in a matrix (e.g.biological networks), and networks which embed inclusions (e.g.cyto skeleton embedding various organelle).
Paper is a network of bonded fibers, with fibers derived primarily from wood.Wood fibers can be broadly classified as softwood and hardwood.Softwood fibers are long and broad and provide strength to the paper structure.Hardwood fibers enhance the surface smoothness and bulk softness.In paper, fibers can bond to each other through six different mechanisms: interdiffusion, mechanical interlocking, Coulomb forces, hydrogen bonding, and Van der Waals forces [5][6][7][8].While hydrogen bonding is often considered a major fiber-fiber bonding mechanism [8], an experimental study by Hirn and Schennach [6] demonstrates that Van der Waals interactions can also be important.Commercial paper (or, any fibrous material) can be classified based on its density (mass per unit volume) and areal density (mass per unit area, or, grammage).The density and grammage of paper and other fibrous products are largely independent.For instance, tissue papers involving hygiene applications have a low grammage and low density; cleaning cotton wipes have a high grammage but low density; cardboards used for packaging applications are required to have a high grammage and high density.Low grammage but high density nano papers have recently been proposed as an eco-friendly substitute for applications requiring high strength, and high electrical and thermal stability [9][10][11].
Non-woven cellulose fiber networks like paper are formed on a moving substrate.Different forming methods exist for commercial paper production, such as Fourdrinier and twin wire forming, as reviewed in [12,13].During the traditional Fourdrinier process (figure 1), a jet of fiber-water suspension is created by pumping a mixture of fibers and water (furnish) through a narrow rectangular slot (the 'slice').The fiber mix consists of fibers of varying sizes and properties, blended according to a weight ratio determined by the manufacturer.Additionally, particulate substances, known as fines and fillers, are often introduced.These fines and fillers A schematic showing different stages in papermaking involving Fourdrinier type forming, dewatering, pressing, and drying.The paper sheet is about 80% water by weight at the end of the forming section, about 50% water after the press section, and about 5% water after drying section.
are significantly smaller in size (< 76 µm) than the constituent fibers and are added to achieve specific paper properties [14,15].The free-surface jet downstream of the slot is directed onto a permeable forming fabric (the 'wire'), moving at high speed (∼100 km hr −1 ).The direction of fabric motion is the machine direction (MD).The ratio of jet speed to the surface speed of the wire usually varies between 0.99 and 1.04 [16].Due to the speed difference between the jet and wire, hydrodynamic shear acts on the fibers in the MD.As a result, the commercially formed sheets manifest a preferential fiber orientation in the MD [17].The fabric is a permeable surface with an associated weave pattern, designed to facilitate water drainage while retaining fibers, fines, fillers, and additives.The draining water exerts a viscous drag on the fibers, causing them to conform around the fabric filaments, resulting in an out-of-plane topographical pattern in the wet sheet.[13,18,19].A special class of paper products called creped tissue papers also exhibit an out-of-plane structure due to the presence of folded micro-structure, which further modify their properties [20][21][22][23][24].In summary, paper can be regarded as a bonded (or crosslinked) network without a matrix with an out-of-plane topography composed of preferentially oriented athermal fibers.
During papermaking on a forming fabric (wire), the formed sheet often manifests wiremark.The wiremark is often categorized as topographical and hydrodynamic.Topographical wiremark is the physical out-of-plane pattern of the fabric embedded in the sheet due to fiber conformation, as explained in the previous paragraph.Hydrodynamic wiremark is caused by the variable spatial distribution of flow through the fabric, which transports different amount of tiny particulate substances (fines [25] and filler [26]) to different parts of the sheet.Hydrodynamic wiremark may be detected as a spatial variation in sheet density or chemical composition, and is generally less visible than topographical wiremark.It manifests mainly on the printed surface; thus, it is only a concern for finer print-grade paper of high grammage [19].From an application standpoint, wiremark is undesirable for uncreped products due to its association with printing issues and unfavorable optical properties [19].In creped tissue paper, the wiremark pattern can influence the crepe-fold frequency spectrum [27,28], thus serving as a controlling factor.In fiber network literature, the effects of micro-architectural parameters such as areal density, fiber orientation, and fiber curl on the network properties have been fairly understood [2,[29][30][31][32][33][34].However, the effect of out-of-plane topography (due to wiremark in paper products) on the sheet tensile properties such as stiffness and strength, when compared to a sheet without any out-of-plane topography, remains an area where our understanding is limited.The objectives of this paper are: (a) to simulate low-grammage formed fiber networks with an out-of-plane topography with relevant physical parameters as inputs, (b) to evaluate the effect of fiber size on the topography and tensile properties of the formed sheets, and (c) to assess the influence of topographical wiremark on the tensile properties of such formed sheets.The effects of hydrodynamic wiremark are beyond the scope of the present study, since two way fluid-structure interaction computations are required.Section 2 briefly reviews the present network modeling techniques, followed by the details of the proposed network-forming model, highlighting the merits of the current model.In section 3, the effect of topographical wiremark on the 2D Fast Fourier Transform (FFT) of the formed sheets, their tensile properties, mass distribution in the thickness direction, and strain energy distribution are investigated.Concluding remarks and directions for further research are offered in section 4.

Discrete element modeling
Forming is a high-speed manufacturing process involving fiber-fiber, fiber-wire, fluid-wire and fluid-fiber interactions.The underlying fluid-structure interactions are complex due to geometrically nonlinear fiber deformations [2], fiber crowding [35], and flow turbulence [36].Modeling the forming process of fiber networks requires the consideration of large fiber deformations, fiber bond-creation during contacts, and changes in the material properties of fibers that result from pressing and drying.Thus, creating a comprehensive forming model is challenging.The finite element method (FEM) [31,37,38] would be prohibitively expensive to use for these simulations for multiple reasons.FEM is challenged to consider hydroelastic interactions involving large deformations.Contact handling and bond formation and disassociation also pose severe computational challenges in FEM.An alternate framework based on the discrete element method (DEM) is appealing due to the ease with which bonds can be created, and their properties changed during the simulation [32,33,39].Also, the DEM framework can be extended to include fluid structure interactions (FSI) [40,41] to model hydrodynamic effects more faithfully.A more comprehensive comparison among FEM, DEM, and FSI based techniques can be found in our previous paper [30].
Our forming model is DEM based with approximate effects of fluid drag and hydrodynamic shear.Figure 2 summarizes the steps in numerical forming with DEM: fiber generation, fiber deposition, pressing (consolidation), and bond creation.Fibers are generated using a beads on a string approach, in which, spherical particles are bonded in series.Fibers can be initially straight or have resting-state curl, kinks, and twist.For these simulations, the fibers are assumed to be straight and of circular cross-section.Once the fibers are generated, a Monte-Carlo algorithm is employed to place the fibers in a given simulation space, which represents the pulp slurry leaving the headbox.Similar techniques are discussed in [42,43].Once the fibers are in place in the simulation box, they represent wet fibers suspended in water, as shown in grey color boxes in figure 2. In a real forming process, fibers suspended in water are driven onto the forming fabric due to hydrodynamic forces (these forces can arise from a suction generated below the fabric or a positive pressure generated in the pulp slurry).The next step in the simulations is to deposit the fibers on a substrate; the substrate can be flat, or have a three-dimensional pattern that represents the forming fabric.The specifications of the substrate used are discussed later in section 2.2.After deposition on the substrate, a flat plate, simulating the press roll in an industrial paper machine, is used to press the deposited fibers Network formation method using DEM involving fiber generation, deposition, and pressing resulting in a network.Note that wet fiber modulus is used in silico till pressing (grey color boxes), which is overridden by dry modulus in later stages to obtain a dry network (white color box).against the substrate.Pressing results in a consolidated, cross-linked wet network.The dry network is then created by changing the wet fiber modulus to that of dry fibers (white color box in figure 2), followed by equilibration.The present DEM model is implemented in the LIGGGHTS®-flexible fibers [44][45][46] software.More details on the DEM formulations can be found in [30,44].The performance of the implemented technique is compared with other forming models [39,47] elsewhere [30].A detailed explanation of the fiber, and bond properties, fabric specifications, and underlying assumptions in the deposition model are discussed in the next subsections.

Fiber modeling and bond properties
In DEM, each fiber, represented as a chain of bonded spherical particles (intrafiber bonds), can deform axially, and in shear, bending, and torsion.Since paper network failure is generally a result of cross-link (interfiber bond) failure [3], each fiber is assumed to be unbreakable.Moreover, we assume fibers to be elastic, and stress-free when straight.The fiber and bond properties used in the simulations are obtained from the literature on dry wood fiber, specifically Eucalyptus Kraft (Euc) and Northern Bleached Softwood Kraft (NBSK), as summarized in table 1.It is important to note that the current framework has the capability to accommodate realistic variations in fiber size distribution by allowing fibers follow to a specific distribution for both length and diameter.However, for simplicity in our current study, we assume that all fibers of each type possess uniform properties, as mentioned in table 1.The bead to bead distance is set to the particle diameter to avoid any interfiber penetration.For this study we mix both types of fibers in two different ratios by weight.For the grammage and sheet area specified in table 1, the total number of fibers for Euc-NBSK ratio of 60:40 is 2139 (33 217 particles).Whereas, for the 40:60 Euc-NBSK ratio, the number of fibers and particles are 1571 and 22 706, respectively.The total fiber count in the sheet (same grammage) with a During the papermaking process, most fiber lumens collapse, resulting in a ribbon-like fiber cross-section, thus having two principal bending stiffnesses.Out of the two, the softer bending stiffness governs the tensile mechanics of low-grammage fibrous structures like paper.The typical physical properties for dry Euc and NBSK fibers are listed in table 1.Since the fibers are cylindrical in the present numerical framework, we chose the fiber diameter by bending stiffness equivalence for the DEM and real fiber.Thus, if E f I f is the value of the softer bending stiffness of a real fiber, the diameter for DEM fiber (d f ) can be found using f /64 for a given fiber modulus E f .Additionally, the bead-spring discretization can introduce artificial fiber roughness [48], potentially leading to increased friction and, in some cases, interlocking, effects of which are shown to be important in non-crosslinked (woven) fiber networks [49].However, the present study focuses on low-density crosslinked networks without entanglements.Additionally, the networks are thin (≈3-4 fiber diameters thick), minimizing the possibility of interlocking.To assess the significance of friction, we analyze the work done by contact and friction forces as a percentage of the total work done on the network during tensile testing.Our observation indicates that the work attributed to friction and contact stays less than 4% of the total work.Considering these factors, the impact of friction is anticipated to be nonsignificant.For the present study, a minor influence of increased friction due to bead-spring approximation may manifest as a slightly extended post-failure region, attributed to increased resistance to fiber pullout.
As highlighted in the introduction, various mechanisms are responsible for interfiber bond formation.The interfiber bonds are assumed to be linearly elastic in all deformation modes.For both the inter and intrafiber bonds, the bond forces and vectors are calculated based of the position and velocities of the bonded particles.For a bonded-pair (ij) between particles i and j, with position vectors r i and r j respectively, the vectors of normal force ( f (ij)  n (t)), shear force ( f , and twisting moment (m (ij) n (t)) in the bond (ij) at time t can be written as: Here r ij is the position of particle i relative to j, i.e. r ij (t) = r i (t) − r j (t).Also, R i and R j are the radii of particles i and j, respectively.The relative velocity vector can similarly be expressed as v ij (t) = v i (t) − v j (t).The relative velocity vector can be broken down into its normal and tangential components.The normal component represents the rate of bond extension, while the tangential component combines the rates of bond shear and rotation.The normal component can be found as v The tangential component can now be found as v A similar exercise can be used to find the rate of twist (Θ (ij)  n (t)) and bending (Θ Using the above calculated normal force, shear force, twisting moment, and bending moment vectors, maximum normal and shear stress in a bond (ij) at time t are calculated as: , and The bond diameter d b is taken as the minimum bonded fiber diameter.Thus, the bond diameter for Euc-Euc, NBSK-NBSK, and Euc-NBSK interfiber bonds are respectively 30 µm, 15 µm, and 15 µm.The failure criterion for interfiber bonds is assumed to be based on maximum strength, i.e. a bond (ij) is considered broken if σ where σ b and τ b are strength in normal and shear respectively.The normal and shear strength, as mentioned in table 1, are taken from the Atomic Force Microscopy tests on softwood fibers [31,50].
Particles can collide with other non-bonded particles, causing inter-particle friction.In the present simulations, the inter-particle coefficient of friction is assumed to be 0.1.In DEM, the contact detection is realized through the neighbor list.All particle pairs within a neighbor cutoff distance equal to their force cutoff (or particle radius) plus the skin distance (maximum diameter of the present particles) are stored in the list.The contact is searched for the particle pairs in the neighbor list at each time step.Contact between two particles i and k is inferred when the distance between them |δr ik | is less than the sum of their respective radii, and interparticle penetration (δ (ik) ) is calculated as between particles i and k is calculated using Hertz's theory of elastic contact [55] as: where E c is the contact modulus and

Forming fabric modeling
In this section we will briefly discuss about the details of the forming fabric used for fiber deposition in both the experiments and the computational model.Forming fabrics are woven structures.Fabric filaments running in the MD are referred to as 'warp' while those in CD are 'weft'.The weave or fabric construction is often expressed in terms of shed, mesh, and count [56].The shed is the temporary separation between upper and lower warp yarns through which the weft is woven.Thus, the surface of a woven fabric is comprised of crests (referred to as 'knuckles') and troughs.It will be shown later that the knuckle patterns play an important role in embedding the out-of-plane topographical structure in the formed sheet.The meaning of warp, weft, shed, and knuckle are schematically represented in figure 3(a).The mesh is defined as the number of MD filaments (warp) per unit length, whereas, count is the number of CD filaments (weft) per unit length.In other words, the reciprocal of mesh and count give the knuckle-spacing in CD and MD respectively.A fabric can be woven in single or multiple layers.For instance, a single layer fabric is composed of one MD filament layer woven with one layer of CD filaments; while in double layer fabric, one MD filament layer woven with two layers of CD filaments.
For the present study, a commercially available 4-shed single-layer fabric with 250 µm filament diameter is used for experimental and numerical forming.The fabric has a mesh of 36 cm −1 , and a count of 21 cm −1 .An MD-CD view of the wire is shown in figure 3(b), along with its 2D-FFT in figure 3(c).The bright points in the first quadrants of the 2D-FFT in figure 3(c) are the dominant spatial frequencies in MD and CD, which correspond to the adjacent filament spacing in either direction, as marked in figure 3(b).As can be seen in figure 3(b), the wire has approximately rectangular pores (the openings shown in white), the largest of which is 0.116 mm × 0.245 mm.Smaller fiber fragments or fillers added to the paper may drain through these pores, although for the large fibers studied here all of the fibers are retained.
The fabric model used in silico is constructed from micro-CT images of a physical forming fabric.Image processing of the micro-CT images yields an STL file of the filament surfaces.This STL file is imported into the DEM platform for numerical simulations.Fabric deformation, which from physical arguments can be shown to be very small, is not considered here.The rigid fabric can interact with fibers through contact and friction.In the present simulations, the fabric-fiber coefficient of friction is assumed to be 0.1.The normal contact force is calculated using Hertz's theory of elastic contact [55] (see equation ( 6)).The model results exhibit general insensitivity to the coefficient of friction.However, excessively low friction can result in fiber slippage during the deposition of the first layer, leading to filtration.The chosen friction coefficient is maintained at a level sufficient to prevent filtration.In reality, friction also contributes to wear and tear in the forming fabric, which can have a slight impact on topography and tensile properties [57].It is important to note that since the fabric is rigid, the current model does not consider the effects of fabric wear.

Fiber deposition model
To examine the effect of fiber size on tensile response, fiber networks are formed with either a 60:40 Euc-NBSK ratio or with a 40:60 ratio.After fiber generation, the fibers are initially suspended in a 3D simulation box.The MD-CD dimension of the simulation box is taken as specified in the table 1, and the height of the box is kept large enough to avert any artificial contact forces caused by initial interfiber contacts.As mentioned in the introduction, the preferred fiber orientation is a complex function of jet, wire, and drainage speeds.Niskanen [17] showed that the orientation distribution of the fibers in the suspension can be approximated as an elliptic function of jet-wire speed difference and drainage speeds.In our simulations, we initially orient the suspended fibers with the same orientation distribution function.A brief explanation of the elliptic orientation function is given as follows.In the MD-CD plane shown in figure 4(a), the probability of finding a fiber at an angle φ from MD can be expressed as: , Here, q is the orientation parameter depending on the jet-to-wire speed difference (u s = |u jet | − |u wire |) and the drainage speed magnitude (|u d |).The orientation parameter q is zero for a randomly oriented fiber distribution (u s = 0), and is otherwise between 0 and 1.The angle θ is the angle of inclination of the ellipse's major axis from MD; for θ = 0 the major axis coincides with MD, while it coincides with CD for θ = π/2.In industrial practice θ is within few degrees of zero.To choose realistic values of q and θ for the simulations, we measure their values from the experimentally formed sheets, as explained in the next paragraph.
Sheets comprised of two different Euc:NBSK fiber ratios (60:40 and 40:60) are formed experimentally, using the fabric mentioned in section 2.2.Both types of sheets have the same grammage (25 g m −2 ), and are formed using sheet former.Figure 4(b) shows a schematic of a formed sheet in the ZD (deposition direction), where t is the sheet thickness.The orientation distribution for layers at different heights z is estimated.The estimated orientation functions for the Euc and NBSK ratios 60:40 and 40:60 are respectively plotted in figures 4(c) and (d) for five z/t ratios.For reference, the circle of radius 1/π (isotropic case) is plotted in black dashes.Figures 4(c) and (d) imply that: (i) there is only a minor change in q and θ in the thickness direction (z); (ii) both sheets have approximately the same values of the orientation parameter (≈ 0.1), with average value of θ = 0. Thus, for the simulations, we choose a constant value q = 0.1 and a constant θ = 0. Note that the experimentally obtained orientation parameters are measured after the sheets are formed, whereas the simulations require the initial fiber distribution.During simulated forming the orientation parameter showed an insignificant variation between the initiation and completion of the forming process.Consequently, the assumption that the experiments can provide a useful initial condition for the fiber distribution is justified.
The different stages of the forming are shown in figure 5. First, the oriented suspension of wet fibers is created in silico using the above experimentally obtained parameters (q = 0.1, θ = 0, a number of fibers corresponding with a basis weight of 25 g m −2 ).The fibers are assumed to be initially at rest.The fabric is then moved upwards (figure 5(a)) through the suspension at a speed |u d | [40,58,59], where |u d | represents a typical drainage velocity of water through a fabric in the forming section of a paper machine.In the industrial paper machine, the bulk velocity of water through fabric while draining is 0.05-0.5 m s −1 .Also, the drainage velocity shows spatial variation in the machine direction (MD), cross-direction (CD), and Zdirection (ZD) due to the complex fluid dynamics involved [16,60,61].For simplicity in our current simulations, we assume a uniform drainage velocity with a magnitude of 0.1 m s −1 .Early in this process, the fibers close to the fabric come into contact with the fabric and begin moving with the fabric speed.When this occurs, these fibers will be acted on by viscous drag (figure 5(b)) from the surrounding fluid, in addition to contact forces with the fabric.
To estimate the drag forces acting on the fibers we use Cox's slender-body theory [16,62].Figure 5(c) shows a fiber at an arbitrary angle α from the MD-CD plane suspended in a uniform flow field of the draining fluid.From Cox's theory the fluid drag tangential and normal to the fiber can be expressed as: Here l f is fiber length, d f is the fiber diameter, µ is dynamic viscosity of water ( = 10 −3 Pa • s at 20 • C), u t is the tangential component of fluid velocity, and u n is the normal component of fluid velocity.In this study we assume that the fibers are constrained to the MD-CD plane (α = 0), which will be approximately true due to the effect of hydrodynamic shear in the paper machine headbox.Therefore, u d ≈ u n , or |u t | ≈ 0. The viscous force on each particle can thus be written as [62]: After deposition, the deposited wet fibers are pressed from the top using a flat plate.During the pressing phase, the formation of cross-link bonds is permitted.A cross-link bond between two particles (i and k) is established if the distance between them (∆ ik ) falls below a specified threshold distance.It is crucial to select the threshold distance judiciously; a higher value might result in artificial cross-linking, whereas too small of a value could lead to minimal or no crosslinking due to high compression requirements.In this study, the threshold bonding distance is chosen to be 0.1% higher than the sum of the radii of the particle pair.Therefore, the bond is formed if ∆ ik < 1.001(R i + R k ), where R i and R k represent the radii of the respective particles.Pressing results in a consolidated, cross-linked wet network.To simulate a wet network formation, we use bond count saturation as a criterion for pressing termination (typical pressing pressures are approximately 1 MPa, which are at the low end of that used industrially); the pressure is gradually removed.The dry network is then created by overriding the wet fiber modulus to that of dry fibers (see table 1), and then allowing the sheet to equilibrate.In the present study, we do not consider the effects of restrained drying, which in a real papermaking process can influence the final sheet properties [63].After equilibration, the networks are subjected to a tensile test in silico to assess the influence of out-of-plane topography on their stiffness and strength.
In summary, below are the assumptions of the present model, along with their potential effects and justifications: • Each fiber is approximated as a chain of spherical particles representing cylindrical fibers (bead-spring approximation).The potential effects and justification of this discretization method are discussed in section 2.1.• Fiber and cross-link bonds are assumed to be elastic in all deformation modes.Fiber bonds (intrafiber bonds) are unbreakable, while cross-links (interfiber bonds) can fail based on the maximum strength criteria (equation ( 5)).• The fiber orientation distribution follows an elliptic distribution function in the MD-CD plane (see equation ( 7)), with negligible variation in the thickness direction.The justification of this assumption is provided in figures 4(b) and (c).• Fibers deposited onto the non-deformable fabric experience a viscous drag due to dewatering flow.This complex flow field is approximated by that of slender body in an unbounded flow field, for which Cox's theory applies (see equation ( 8)).• Fibers interact with other non-bonded fibers and the fabric primarily through normal contact and friction.To compute the normal contact force, we apply Hertz's theory of elastic contact [55].• The networks are assumed to have no initial residual stresses, which is an ideal case.The residual stresses are known to distort the shape and tensile curve of the sheet [64].
Note that, the method presented above is a single-wire forming representation, as the fiber deposition takes place from one side.The present method can also be extended to the twinwire forming method in which two wires, with fibers between them, move towards each other.

Results
Dry networks with MD-CD dimensions 5.5 mm × 2.5 mm are formed using the above method.
To avoid the edge effects, a region of interest (ROI) is selected for calculations of grammage, thickness, and tensile force, as will be discussed later.The method of ROI selection is explained in appendix A. The simulation results are presented in two sub-sections.Section 3.1 discusses the effect of the fiber-mixing ratio on the topography and tensile responses of the formed sheets.Results for 40:60 and 60:40 ratios of Euc-NBSK (by weight) are presented.The numerical results are also compared with the experiments in this sub-section.Section 3.2 discusses the effect of out-of-plane topography on the ZD mass variation, tensile responses, and strain energy distribution for the numerically formed sheets.

Effect of fiber mixing ratio: model vs experiments
In this sub-section, we study the effect of two fiber mixing ratios on the topography and tensile response of the sheets, using experiments and modeling.The sheets are experimentally formed on the fabric; both of the formed sheets have the same grammage (25 g m −2 ) and network thickness (110 µm), and thus the same network density.The sheets with the same fiber ratios are also formed numerically in-silico following the steps mentioned in the previous section.topography (away from the viewer), while the red represents peaks in the sheet topography.As can be seen, the MD-CD views for the sheets with two mixing ratios are almost identical.The 2D-FFTs for both ratios are found to be identical, therefore the 2D-FFT for only 40:60 Euc:NBSK is shown here.Figures 6(e) and (f) respectively show the same for fabric and the opposite sides.It can also be observed that the 2D-FFT of the formed sheet closely resembles that of the fabric (figure 3(b)), especially along MD (abscissa) and CD (ordinate) axes.Thus, it can be inferred that the blue troughs (depressions) in figures 3(a)-(c) are caused by embedment of the fabric's knuckle patterns (see section 2.2) in the formed sheet.
From the MD-CD views and the 2D-FFTs, it can be observed that the topographical wiremark is visible on both the sides of the formed sheets.However, a comparison of the color gradients in the MD-CD views of fabric and opposite side indicates that the fabric side (figures 6(a)-(c)) shows a larger amplitude of thickness variation (blue to red), as compared to the opposite side (yellow to red color variation in figures 6(b)-(d)).Also, the intensity of bright points in the 2D-FFT in figure 6(f) show a reduction as compared to that in figure 6(e).The diminution in the amplitude of the sheet topographical variations from the fabric side to the non-fabric side is believed to be due to the smoothing effect of the (roughly three) fibers that comprise a typical sheet thickness.The surface image of the experimentally formed sheet and its 2D-FFT is shown in figures 6(g)-(i).The procedure for taking surface images and evaluating their 2D-FFTs is explained in [28,65].The experimental 2D-FFTs are in qualitative agreement with the simulation results, and also show the two-sidedness of the wiremark.Furthermore, the FFT peak amplitudes in figure 6(h) are higher than in figure 6(i), under the same lighting conditions, which is consistent with in silico observation of less pronounced wiremark on the non-fabric side of the sheet.In the limit of very high grammage sheets, wiremark is expected only on the fabric side.The results presented so far indicate that the fiber mixing ratio does not have a significant effect of out-of-plane topography of the formed sheets.However, it is not the case with their tensile responses, as will be shown in the subsequent paragraphs.
We now study the effect of two fiber mixing ratios on the tensile response of the sheets, using experiments and modeling.Tensile tests were performed (INSTRON-5969) following TAPPI-T494 standard [66].Multiple tensile tests (n = 7) are performed for a single fiber mixing ratio; the averaged response is reported here.The numerically formed sheets are subjected to a virtual uniaxial tensile test along the Machine Direction (MD).Similar to the tensile tests, one end of the network is rigidly clamped in a vise, while the network is pulled along MD at a constant rate (14% strain per minute).The cross-links fail during numerical tensile test following equation (5). Figure 6(j) summarizes the effect of fiber mixing ratios on the tensile response of the sheets.The statistical and size convergence of the numerical simulations in shown in appendix A. Note that the ordinate in figure 6(j) shows tensile force per unit sheet width (in CD).It is observed experimentally and numerically that increasing the fraction of NBSK fibers (longer and wider fibers) increases the sheet strength and failure strain for the same grammage.However, the change in the initial stiffness (slope till 0.05% strain) of the sheets is not significant.
The two simulated networks share the same areal density (grammage) of 25 g m −2 .The mean fiber length for the 60:40 (Euc:NBSK) network is 1.08 mm, while it is 1.22 mm for the 40:60 (Euc:NBSK) network.The observed increase in strength with higher NBSK content in figure 6(j) can be attributed to the corresponding increase in fiber length.Due to the fixed grammage, the number of fibers in the network decreases by 22.5%.These findings align with the analysis conducted by Komori et al [67].Interestingly, the two networks exhibit minimal variation in their bonded area per fiber (≈0.2%) and mean segment stiffness (≈0.4%).As bonds can be viewed as additional stiffness parallel to fiber segments, the fact that both types of fiber segments have the same stiffness and parallel stiffness due to bonds may result in a non-significant variation in the initial stiffness of the network.
The comparison of the experimental and numerical tensile responses show a good agreement given the myriad approximations of the model (simplified fluid dynamics, monodisperse fiber length and diameter, idealized cylindrical fibers, uniform bond strength, etc).Specifically, the initial slope of the simulated tensile response is within 10% of the measured response.Although, the simulation overpredicts the strain to failure of the experiments, the disagreement is less than 15%.One obvious feature of disagreement between the simulations and experiments is that the experiments show a softening tensile response with smaller strain stiffness as compared to simulations.Wood fibers are known to exhibit plasticity [2,31], which may be a reason for the yielding of experimental tensile response.In contrast, the simulations assume elastic fibers and interfiber bonds, for which the simulated response is almost linear.
Studies on random fiber networks [68] have revealed size effects in tensile properties.However, finite element (FE) studies on cellulose paper sheets by Brandberg et al [69] demonstrated less sensitivity to sample size.This reduced sensitivity is attributed to the treatment of interfiber bond regions as deformable contact elements, in contrast to the approach of Shahsavari et al [68], where these regions were treated as merged finite element nodes.The study by Brandberg et al [69] indicates that mean tensile properties converge effectively for networks larger than 2 × 2 fiber length.Similarly, the discrete element method (DEM) based networks presented in this study also exhibit reduced sensitivity to sample size, as illustrated in appendix A. This diminished sensitivity can be attributed to the modeling of interfiber bonds as breakable and deformable in axial, shear, bending, and twisting modes.

Effect of out-of-plane topography
Figure 7 compares the sheet formed on a fabric to that formed on a flat substrate.A flat substrate represents a fabric with an extremely fine weave structure (small wavelength and small filament diameter compared to fiber diameter).Although the comparison is shown for the sheets with 40:60 Euc:NBSK ratio, the conclusions described below also apply for the sheet with 60:40 Euc:NBSK ratio.The probability density histogram for finding a particle in the thickness direction (ZD) is shown in figure 7(a).In this figure the height in the ZD (z) is nondimensionalized by the diameter of a Euc fiber (d f = 15 µm); z = 0 represents the side in direct contact with the fabric.The histogram for the sheet formed on the flat substrate is relatively symmetrical and uniform (excluding edges).However, for the sheet formed on the fabric, the histogram first gradually increases until z/d f ≈ 2.5, then increases rapidly to a maximum value, followed by a sudden drop as the opposite side of the sheet is approached.For both the fabric and flat sheets, the histograms in the sudden drop region are similar, perhaps because the top of the sheets are pressed with a flat plate in both cases.The difference between the histogram on the fabric side and the non-fabric side for sheet formed on a fabric is further evidence of the reduced prominence of topographical variation in the thickness direction with increasing distance from the fabric, as discussed above.The network formed on the fabric has a higher thickness (110 µm) than the one formed on the flat substrate (65 µm) for the same grammage, and thus has a lower apparent density.In his work, Adanur [70] also found the sheet thickness to decrease with a finer weave.
The tensile response of a 40:60 (Euc:NBSK) sheet, formed on the fabric, is compared in figure 7(b) with a sheet formed on a flat substrate.As shown, the sheet formed on the fabric has lower stiffness and lower tensile strength, as compared to the sheet formed on a flat substrate.This suggests that as the weave of the fabric becomes finer (or out-of-plane topography is less pronounced), the stiffness and strength increase.Experimentally, the structure of forming fabric is often characterized by the fiber support index (FSI) [71] and drainage index (DI) [72].FSI is based on the mean fiber length and average number of supports per fiber; higher FSI corresponds with a finer fabric structure.The DI is based on fabric air permeability, CD Mesh count and FSI.Higher values of DI reflect better drainage.Adanur [70] found that an increase in DI and FSI resulted in an increase in sheet strength, which is in qualitative agreement with our results.The flat network shows 40% higher bonded area as compared to the one with outof-plane topography, resulting an increase in the tensile strength.Figure 7(c) presents the bond fraction, which is the ratio of the number of broken bonds to the initial total number of bonds, as a function of strain.In the flat network, bond rupture progresses more gradually, whereas in the network formed on fabric, it exhibits a more intermittent pattern.The macroscopic failure of the sheet takes place after approximately 3% of the bonds have ruptured in the case of the flat network, but it occurs after approximately 5.5% of the bonds have failed in the networks formed on fabric.The abrupt change in the tensile responses depicted in figure 7(b) can be linked to the sudden decrease in bonds within the network, occurring at 0.5% strain for the fabric and 0.65% strain for the flat configuration.The experimental tensile response for the 60:40 Euc:NBSK ratio in figure 6(j) also demonstrates a kink near 0.5% strain, however, with a much smaller magnitude.The simulations exhibit a larger kink magnitude, possibly attributed to the smaller size of the sheet compared to the experimental setup.The precise cause of the kink remains inconclusive at this point.
The evolution of strain energy fraction with macroscopic strain for the above sheets is shown in figure 7(d).The fibers store more strain energy (SE) than do the bonds.For the sheet formed on the flat substrate (red curves), the SE stored in fibers initially gets transferred to bonds as strain increases, resulting in a reduction in SE fraction in fibers with increasing macroscopic strain.However, after a particular strain value (0.6 % for the present network), bond rupture initiates, and the SE fraction carried by the bonds begins to reduce.Finally, the network loses its stiffness around 0.8%, and the sheet fails (red curve in figure 7(c)).In contrast, for the sheet formed on the wire, the SE fraction in fibers decreases (and vice versa for bonds) until a very small strain ∼ 0.1 %.At higher strains, the SE fraction for fibers shows only a mild increase with the macroscopic strain.The sheet fails at a macroscopic strain of around 1%, and a large fraction of SE transfers from bonds to fibers.The mild variation in SE fraction in the case of the sheet with an embedded pattern suggests that the rate of SE transfer from the bond to fiber is insignificant.One reason may be the bending and torsion deformation mode activation due to out-of-plane topography, as explained in appendix B. The partition of fiber strain energy for the low strain (0.1%) and high strain (0.7%) values is shown in figure 7(e).At the low strain values (0.1%), the bending mode of deformation dominates in the sheets with out-ofplane topography.However, more energy accumulates in the axial mode for the sheets formed on the flat substrate.The differences due to topography become more contrasting at higher strains (0.7%).Significant energy accumulates in bending and twisting modes for the wireformed sheets; in contrast, almost all energy is stored in the axial mode in networks formed on flat substrates.A simple graphical explanation for the SE storage in bending and twisting mode for the sheets formed on the wire is presented in appendix C. Note that, as suggested in appendix C, the twisting mode gets activated exclusively due to out-of-plane topography.The observations in figures 7(d) and (e) suggest that in the case of sheets formed on fabric, the SE in fibers and bonds is nearly segregated.One reason for the segregation may be the activation of softer deformation modes (twisting and bending) in fibers due to the out of plane topography.

Conclusion
In this paper, we simulate the formation of low-grammage networks with out-of-plane topography.We investigated the effect of fiber size and out-of-plane topography (topographical wiremark) on the tensile response of a formed sheet of wood fibers.We use a DEM-based network formation technique with approximated effects due to hydrodynamic shear and viscous drag.Although the forming method is approximate, the formed networks have simulated tensile curves that are in good agreement with experiments.Key conclusions from the study are summarized below: • Increasing the fraction of softwood fiber (longer fiber with larger diameter) in the sheet increases its strength and strain to failure; however, the initial stiffness does not change significantly.Also, changing the fiber fraction does not show a significant impact on the out-of-plane topography of the formed sheets • The sheets with topographical wiremark formed on a fabric have a higher thickness and an asymmetric mass distribution in the thickness direction, as compared to the sheets formed on a flat substrate (fabric with a finer weave) • The tensile response of a sheet with pronounced topographical wiremark is softer with a lower tensile strength • Axial mode of fiber deformation is always salient in the tensile tests, but in the sheets with topographical wiremark bending and twisting modes of fiber deformation are also significant • Due to the activation of other deformation modes in the fibers, the strain energy (SE) transfer between fibers and bonds reduces during the tensile test.Thus, unlike sheets made on a flat substrate (fabric with an extremely fine weave), SE fraction variation is non-significant in sheets with embedded patterns The effects of fiber plasticity and the spatial variation of viscous drag have not been considered in this study.Future work will consider the plasticity of fibers, but an accurate consideration of the viscous drag acting on fibers would require simulations that are not currently being contemplated.It is essential to acknowledge that this study does not consider the influence of residual stresses resulting from sheet drying.While the simulated dry sheets possess some residual stress due to the overriding of dry fiber properties, the equilibrated samples exhibit an average residual stress of less than 1% of the mean stress near the network failure.Hence, the numerically tested samples can be considered as approximately stress-free.However, in reality, drying is recognized to impact the properties of the final sheet, attributed to factors such as micro-compression and the realignment of microfibers within the fiber network [63].Consequently, alterations in the fiber network topography due to drying are also expected, influencing tensile properties.However, this aspect is considered beyond the scope of the current study.We expect the presented conclusions to hold qualitatively after consideration of drying effects.Furthermore, although our study assumed constant cross-sectional fibers, the methodology employed here can be extended to accurately model the lumen as well as kinks and crimps in the fibers.Finally, the wiremark is known to govern the web-Yankee adhesion in the dry creping process used to produce tissue paper [20,27,73].Our future work will include in silico creping of a formed sheet to understand how the fabric structure affects the properties of the creped sheet.From a practical standpoint, a potential extension of our current research is the identification of an optimal fabric structure tailored to specific paper properties.

Appendix B. An analytical explanation for the strain energy variation in flat network and network with out-of-plane topography
In this section, we aim to elucidate the phenomenon depicted in figure 7(d).When a network forms on a flat substrate, the majority of the fibers remain in a straight configuration.In figure B1(a), a schematic representation illustrates two bonded straight fibers intersecting at an angle denoted as ϕ.We can approximate this intersection as a series combination of two spring-like elements: k f , representing the stiffness associated with the oriented fiber, and k b , representing the stiffness due to the bond.It's important to note that the stiffness k f is dependent on both the properties of the fiber and the angle of intersection ϕ.For loading in the horizontal direction, k f reaches its maximum value when ϕ = 0, and it decreases towards a minimum when ϕ = π.We can express this variation using the equation k f (ϕ) = k fm cos (ϕ/2).The fiber-bond stiffness ratio for a network formed on a flat substrate can now be defined as: The fraction of energy stored in the fiber and bonds can now be written as Upon the formation of a network on a fabric (wire), the fibers acquire an out-of-plane topography, as schematically illustrated in figure B1(b).In this scenario, we can approximate the fiber bond as a triangular link connected by bending stiffness.Using the formulations found in [23,24], we can express the bending stiffness as k θ = k f (t/l) 2 /(6 cos 2 θ).The fiber-bond stiffness ratio for network formed on a wire can now be defined as The variation in energy fraction within the network formed on the wire substrate exhibits much less sensitivity to changes in ϕ (or macroscopic strain) compared to the network formed on the flat substrate.Additionally, please note that this analysis does not account for the effect of fiber twist, which introduces an additional mode for storing strain energy.Consequently, the milder variation in strain energy fraction observed in simulations may be attributed to fibers serving as a relatively softer deformation mode, as other deformation modes are also activated alongside axial deformation.

Appendix C. Deformation modes in an oriented fiber with out-of-plane topography
For a sheet formed on a flat substrate, the fibers can be idealized as straight fibers, lying (approximately) in the MD-CD plane.For a sheet under a tensile test along MD, the force will transfer to the fibers through the bonds.A fiber segment oriented at an angle ϕ from MD, under force f in MD, is shown in figure C1(a).The force can be resolved into axial and transverse components, as shown.Here we assume the twisting stiffness of the bonds balances the couple formed by the transverse component of force (shown by the curled arrow at the origin); just the force components are shown for clarity.The axial component (f cos ϕ) will result in axial SE, while the transverse component (f sin ϕ) will result in shear and bending SE.However, the same fiber is now bent (due to out-of-plane topography) in ZD as shown in figure C1(b).The fiber is now subjected to additional bending moment (∝ fz cos ϕ) and twisting moment (∝ fz sin ϕ), thus, increasing the SE in bending and twisting, respectively.Note that when fibers have in-plane curl, the contribution from twisting mode disappears since z = 0.

Figure 1 .
Figure 1.A schematic showing different stages in papermaking involving Fourdrinier type forming, dewatering, pressing, and drying.The paper sheet is about 80% water by weight at the end of the forming section, about 50% water after the press section, and about 5% water after drying section.

Figure 2 .
Figure 2.Network formation method using DEM involving fiber generation, deposition, and pressing resulting in a network.Note that wet fiber modulus is used in silico till pressing (grey color boxes), which is overridden by dry modulus in later stages to obtain a dry network (white color box).

Figure 3 .
Figure 3. (a) Schematic showing the definition of weft, warp, and shed for a forming fabric.The eye at the top represents the direction of observation for knuckle definition; (b) top-view of forming fabric used for the fiber deposition in this study, showing filament spacing in MD (A) and CD (B).Scale bar = 1 mm; (c) 2D-FFT of the top-view of the forming fabric showing fundamental spatial frequencies.Note that both MD and CD filament spacing are also reflected in 2D-FFT, marked as points A and B.

Figure 4 .
Figure 4. (a) General representation of an elliptic orientation function; (b) schematic of network in MD-ZD plane, t is the network thickness and z is the height from the bottom.Orientation distribution plots for different z/t for (c) 60% Eucalyptus, 40% NBSK fibers, and (d) 40% Eucalyptus, 60% NBSK fibers.The dashed circle represents the isotropic case (q = 0).

Figure 5 .
Figure 5. Moving wire method for network formation: (a) fibers are initially (t = 0) suspended freely in a viscous medium, with MD-CD orientation described by equation (7), while the wire is moves upwards at speed |u d |; (b) at time t > 0, oriented fibers come in contact with the wall and are deposited under viscous drag described by equation (8); the deposited fibers are shown in dashed brown color; (c) the normal (f n ) and tangential (f t ) components of drag force on a fiber suspended in MD-ZD plane in uniform flow field u d .

Figure 6
shows the topography of numerically and experimentally formed sheets with 40:60 and 60:40 Euc-NBSK ratios (by weight), their 2D-FFTs, and their tensile responses.The top views of the side in contact with the fabric and the opposite one are shown in figures 6(a) and (b) for a sheet with 40:60 Euc-NBSK ratio.Figures 6(c) and (d) show the same for the sheet with 60:40 Euc-NBSK ratio.The blue color in figures 6(a)-(d) represents the troughs in the

Figure 6 .
Figure 6.Topographical wiremark of simulated formed sheet with 40:60 Euc:NBSK ratio on (a) fabric side and (b) non-fabric side.Simulated topographical wiremark with 60:40 Euc:NBSK ratio on (c) fabric side and (d) non-fabric side.2D-FFT of 40:60 Euc:NBSK for (e) wire side and (f) opposite side; (g) Surface image of an experimentally formed sheet along with 2D-FFT of the fabric and non-fabric sides is shown in (h) and (i) respectively; (j) effect of Euc:NBSK ratio on the tensile response (experiment and simulation).

Figure 7 .
Figure 7.The effect of fabric geometry for 25 g m −2 sheets with 40:60 Euc:NBSK ratio.These sheets, which were formed on a fabric, are compared to those created on a completely flat substrate in terms of the following aspects: (a) mass distribution in ZD; (b) the tensile curve of the sheet; (c) cross-link count; (d) strain energy (SE) fraction in fiber and bonds, and (e) strain energy partition in fibers.

A1.
(a) Density variation for sheet with 40:60 Euc-NBSK ratio along MD and CD for ROI selection; (b) domain-size and statistical convergence of the numerical simulations.Six realizations are shown for each sheet size mentioned in the legend.

Figure B1 .
Figure B1.Fiber intersection and equivalent lumped system for (a) flat network and (b) network with out-of-plane topography; (c) the strain energy fraction stored in fibers and bond for flat networks and networks with out-of-plane topography.

4 )
The fraction of energy stored in the fiber and bonds can now be written as: U T is the total strain energy stored.Equation (B.5) is shown in figureB1(c) in solid red, while equation (B.6) is shown in dotted red.It's worth noting that the highlighted green region in figureB1(c) bears a resemblance to the pre-failure region in figure7(d).

Figure C1 .
Figure C1.(a) A straight fiber lying on MD-CD plane and (b) a fiber with out-of-plane topography in ZD, under a force in MD.