Network model for predicting structural properties of paper

Introduction

Papermaking has been around for millennia. Though paper is ancient, the material and process to make it never stopped evolving. A lot has happened over the years, from using wooden presses to industrial machinery. Generations of experimentation and development have led to different paper products with a broad range of applications. Well into the digital age, simulations are becoming a realistic approach. Simulation is preferable, as testing new approaches on the ever-advanced industrial-sized machines are expensive in terms of resources, energy, and time.

Analyzing paper is challenging because of its heterogeneity. This complexity arises from the complex fibrous structure paper has on the microscopic level. These individual paper fibers are connected through mechanical interlocking and microscopic forces and form the paper (Hubbe 2007, Hirn and Schennach 2015, Orgéas et al. 2021).

There is a rich history of analyzing paper’s properties using continuum models. Cox presented 1952 (Cox 1952) formulas for deriving the elastic constants for straight fiber networks with non-uniform fiber orientation. This theory was expanded on by, among others, Perkins et al. in 1978 (Perkins et al. 1978), which took into account density and bond properties. In 1969, Page introduced the Page equation (Page 1969), which relates the tensile strength of a sheet of paper to the zero-span tensile strength, density, bond, and pulp properties. More recently, continuum models for analyzing paperboard have been developed (Borgqvist et al. 2015, Robertsson et al. 2021), and finite element models have been used by Lindberg and Kulachenko (2022) to model paperboard packaging.

Continuum models are usually dependent on parameters based on the structure of the paper that are derived experimentally or by analyzing the structures of simulated paper geometries. Another approach to continuum models is to simulate the entire microstructure of paper by modeling the individual fibers. In Hamlen (1991), paper permeability and mechanical properties were evaluated using linear elastic models on two-dimensional and three-dimensional network models. These simulations were deemed too demanding at the time, and he proposed removing less important deformation modes for efficiency. Heyden (2000) used network models to analyze the mechanical properties of different cellulose structures and domains, where individual fibers are modeled as Bernoulli beams. The forming and pressing process of paper was simulated by modeling individual fibers in Lavrykov et al. (2012). In Tojaga et al. (2021), individual fibers were modeled using elastoplastic Timoshenko beam finite elements with a quasi-brittle fracture model to recreate the nonlinear stress-strain response of paper.

A paper model that resolves the individual fibers can evaluate various mechanical properties, as the model is not inherently dependent on the specific test. The paper model presented in this work aims to act as a tool for paper developers. This tool might be helpful for evaluations of non-standardized mechanical tests or evaluations of fiber composites. Ideally, a papermaker could create a virtual paper model using a composite of pulps from a database and then simulate various structural experiments and get representable results. However, for this tool to be viable, it needs to be computationally feasible, predictable, and consistent with paper theory.

The paper model used in this work was developed to be simple and efficient. The model is similar to that of Hamlen (1991). Both are linear network models, model bonds as beam elements, and have the same approach to strength simulations. The main differences between the two models are that Hamlen (1991) handles moments as separate variables, and in this work, bonds are modeled as multiple beam elements. While more advanced models have been developed to analyze the paper in detail (Borodulina et al. 2012, Tojaga et al. 2021), this work aims to evaluate the usefulness of a simplified and computationally lenient fiber-based approach.

The implementation and validation of the network model in this work were developed in the industrial collaborative project called ISOP (Innovative Simulation of Paper). ISOP is performed by a consortium consisting of Albany International, Stora Enso, and Fraunhofer-Chalmers Center. The project’s end goal is to develop tools that are actively used in the product development of paper-based products. Within this project, paper forming has been simulated and evaluated by Svenning et al. (2012). In Kettil et al. (2020), the constitutive model used in this article is presented in detail, and in this article, the physical properties of the model are validated.

This article shows that the model can produce representable results for bending resistance, tensile stiffness, and tensile strength in both machine and cross direction using the proposed linear network model. The parameters of the model used for these validations are either from experimental data available at the paper product companies or by deducing reasonable values for the parameters using published results. Experimental results validate the models accuracy for a range of low-grammage sheets. This validation is strengthened by comparing the simulated results to predictions using continuum models. Moreover, the simulations are performed multiple times for different generated models to validate the approaches’ stability. The simplified nature of the model enables fast validation of macroscale properties on an off-the-shelf desktop computer. These resources and simulation times are accessible to a paper developer in the industry and may act as an intermediate tool between known papermaking identities.

Constitutive model

The edges of the spatial network model are modeled as beams using the equations presented in Kettil et al. (2020). These linear equations model axial stiffness (edge extension) and bending resistance (angular deviation) in the individual edges and provide a linear relationship between the displacement of the network structure and the forces resulting from the displacement. The resulting linear problem can be represented by the following linear system:

(1) K u = F , u = [ u 1 , u 2 , … , u n ] T , F = [ F 1 , F 2 , … , F n ] T ,

where n is the number of network nodes, K ∈ R 3 n × 3 n is the connectivity matrix, u i = [ u i , x , u i , y , u i , z ] are the three-directional nodal displacements, and F i = [ F i , x , F i , y , F i , z ] are the directional applied forces in the i:th network node.

Edge extension defines the forces arising in the change in the length of the edges. These changes are quantified by their strains ϵ = Δ L L , where L is the initial length of an edge and Δ L is its change. The strains used in the model are linearized by using orthogonal projections (see Figure 1). Using the projected strain in Hooke’s law, the magnitude of the force is approximated by:

(2) F = E · A · ϵ = E · A · Δ L L ,

where E is the axial stiffness (Young’s modulus), A is the cross-section area, and F is the quantified force arising from the displacement. The directions of these forces are chosen to be parallel to the edges’ initial axes, again for linearity. An illustration of the forces arising from a displaced edge by edge extension can be seen in Figure 1.

Figure 1

Edge extension forces that arise from a displacement. The change in length of the edge is calculated by projecting the displacement orthogonally with respect to the initial position.

Angular deviation provides a counterforce when the angle between two connected edges is changed. Two connected edges are called an edge pair. This angular change is evaluated in two planes, one plane containing the edge pair (in-plane) and the other orthogonal to the first (out of plane), and is handled similarly. For each of these planes, the deformation of an edge while fixing the other can be modeled using Euler-Bernoulli beam theory:

w ″ ( t ) = − M ( t ) E I = − F ( L − t ) E I ,

where w ( t ) and M ( t ) are displacement and the moment at the point t along the deformed edge, E I is Young’s modulus times the second moment of area, and L is the length of the edge. In the center-point of the edge pair, t = 0, this becomes:

(3) w ″ ( 0 ) = − F L E A I ⇔ F = − E A I · w ″ ( 0 ) L .

The displacement w in the selected plane in the network model is evaluated in the three nodes of the edge pair using projections (see Figure 2), and w ″ ( 0 ) is evaluated using the second-order central difference scheme on these three evaluations. With the second derivative in the node between the two edges approximated, the amplitude of the angular deviation forces are calculated by (3) and given directions in the plane, normal to the edges (see Figure 2).

Figure 2

The forces arising from the angular deviation equation given a displacement.

Experimental data

In the validation of the paper model, fourteen low-grammage sheets of paper are used as a reference. These sheets are single-ply, grammages in the range of 30–60 g/ m 2 , and with different wood fiber compositions. Over three days, these sheets were produced as a pilot project with Albany International and Stora Enso at The Packaging Greenhouse in 2012. The machine used was a pilot board machine with a standard configuration.

Experimental results on these low-grammage sheets conducted by Stora Enso are used for constructing the model and as the basis for the validation. Each paper was cut into smaller samples, and each sheet’s different structural and mechanical properties were analyzed ten times. Averages of these properties can be found in Table 1 and Figure 3. These properties are grammage, thickness, tensile stiffness, tensile strength, and bending resistance.

Table 1

The experimental data of the sheets produced at the Packaging Greenhouse used for the validation of the linear network model.

Properties	Values
Paper sheets analyzed	14
Samples per sheet analyzed	10
Sheet thickness (t)	28–54 µm
Sheet grammage ( g s )	41–67 g/ m 2
Sheet density, ρ s = g s / t	610–810 kg/ m 3
Elastic modulus, MD ( E MD )	5.4–8.5 GPa
Elastic modulus, CD ( E CD )	2.6–5.8 GPa
Stress at break, MD	41–68 MPa
Stress at break, CD	16–43 MPa
Bending resistance, MD	6.6–32 mN
Bending resistance, CD	2.8–32 mN

Figure 3

Mean structural properties of the indivdual sheets in the experiment presented in Table 1. The different sheets are color-coded to their wood composition, where SW stands for softwood and HW stands for hardwood.

Macroscale structural properties of the sheets analyzed are thickness and grammage. All of the sheet samples tested were weighed and measured in advance. Following the standard SCAN P:88-01, each sample sheet’s STFI-thickness was measured. Grammage is calculated by analyzing the weight of the different samples and provided for each sample.

The mechanical properties analyzed are tensile stiffness, tensile strength, and bending resistance. These experiments were performed under a controlled environment of 23 °C and 50 % relative humidity. Each mechanical experiment was performed ten times for different samples of each sheet in both machine and cross direction. The tensile experiments were performed on 15 mm wide paper samples, with a test span of 100 mm using an L&W Tensile Tester following the ISO 1924-3:2005 standard. For the bending resistance experiment, an L&W 15 ° machine was used following the ISO 2493-1:2010 standard on 38 mm wide samples, with the optional 10 mm bending length.

Model parameters

The macroscale parameters of the model are thickness, grammage, and fiber orientation. The thickness and grammage of the models are determined by the measurements performed on the sheets. The fiber orientation was not measured experimentally. In the model, the fiber orientation takes the 1-Cosine distribution (Cox 1952) with a 1 = 0.51, i. e.

Ψ ( γ ) = 1 π 1 + 0.51 cos ( 2 γ ) .

This value of a 1 was chosen such that the relation from Cox (1952):

E MD E CD = 6 + 4 a 1 6 − 4 a 1 ,

is consistent, on average with the experimental results E MD / E CD = 2.035, where E MD and E CD are the Young’s modulus in the machine direction and cross direction of the paper. This choice of a 1 is also within the range of experiments performed in Wahlström (2009).

The sheets analyzed are composed of different wood types. In Figure 3, the wood composition compared to the different structural properties of the sheets is presented. Slight variations of the structural properties are observable for the different wood types. These variations are out of scope for this initial validation, with only softwood fibers modeled.

Detailed experimental data govern the length and width distributions of the fibers in the model. Stora Enso collected this data in 2020 by analyzing a softwood pulp different from the structural experiments using L&W’s Fiber Tester Plus. The experimental data is presented in Table 2, which defines the geometric distribution for all the fibers in the simulation.

Table 2

Geometrical breakdown of a softwood pulp, provided by Stora Enso.

Properties	Values
Coarseness ( c f )	188 µg/m
Average length (L)	2 mm
Average width ( w f )	29 µm
Fiber length 0.2–0.5 mm	8 %
Average width	23 µm
Fiber length 0.5–1.5 mm	27 %
Average width	32 µm
Fiber length 1.5–3 mm	45 %
Average width	41 µm
Fiber length 3–4.5 mm	20 %
Average width	46 µm

The microscale parameters are the geometrical and structural properties of the paper fibers. Data from the pulp analysis (Table 2) defines the width and length distributions of the fibers, along with the coarseness. The cross-section areas of the fibers are not known from experiments. Instead, the cross-section area is based on interpolating the measurements performed by Horn (1974) by assuming a linear relationship between the cross-section area and coarseness. The cross-section area used in the model was chosen so that optimum results were achieved and, at the same time, stayed consistent with the other measurements. A graph illustrating this interpolation and the choice of cross-section area is provided in Figure 4. The known properties (coarseness, width, length) of the softwood pulp used in the paper model are indicative of earlywood (Levlin and Söderhjelm 1999 pp. 21–24), and thus the cell wall thickness of the fibers in the model was chosen to be 3 µm.

Figure 4

The cross-section area used (115 µ m 2 ) in the model (circle) compared to measurements in Horn (1974) (crosses) and the linear regression of the measurements (dashed line).

The structural parameter of the fibers: Young’s modulus and tensile strength, are chosen using published results. Young’s modulus is set to 35 GPa, in the range of Neagu et al. (2006) and Wang et al. (2011). For tensile strength, Wang et al. (2011) found tensile strengths for fir in the interval 0.5 GPa to 1.2 GPa, similarly to Groom et al. (2002), with 0.4 GPa to 1.3 GPa depending on age and position in the tree. Our model achieves optimum results for 0.6 GPa, which is consistent with those findings, and thus that value was chosen. For an overview of all the fiber-related parameters, see Table 2 and Table 3.

The last set of model parameters is related to bonding. Bonding parameters scale with respect to the bonding resolution in the model. The strength and stiffness of a bond should be relative to how bonded two fibers are but not dependent on how many bond edges are used to represent a bond. Therefore, the cross-section of a bond is scaled based on the number of edges representing a bond, more specifically, the width of the volume subsections used in the bonding process (bond delta). Beams with rectangular cross-sections represent the individual bond edges in the model, and they are 3 µm thick and wide as the length of the subsections. Bonds have angular stiffness in only one direction to not over-constrain the model. Young’s modulus of these bond edges is the same as the fiber edges, with the tensile strength being substantially lower at 45 MPa.

Simulation setup

In the tensile type simulations, two opposite sides of the model are clamped, where one is displaced to introduce stress (See Figure 5). The Dirichlet conditions are imposed explicitly by setting the displacement of all the network nodes (in all coordinate directions) in the clamped regions to the given displacement. For the tensile stiffness simulations, the models represent a 4.2 mm × 4 mm piece of paper, where two sides are clamped 0.1 mm, resulting in a 4 mm × 4 mm active domain. With the discretization parameters in the model (presented later), the 0.1 mm clamped domain guarantees that most fibers in those areas are fixed in at least two points. Adding additional length to the clamped domain results in little to no change for this model with these parameters.

Table 3

Parameters used in the linear network model that were deduced.

Properties	Values
Fiber cross-section area ( A f )	115 µm2
Fiber Young’s modulus ( E f )	35 GPa
Fiber strength	0.6 GPa
Fiber orientation, 1-Cosine ( a 1 )	0.51
Cell wall thickness	3 µm
Fiber density ( ρ f = c f / A f )	1635 kg/ m 3

Figure 5

Tensile simulation setup.

The model is strained to 0.5 % in the tensile stiffness simulation, and the forces arising are used to calculate the stiffness of the given paper model. In the tensile strength simulation, the strain is increased iteratively. The simulation strains the paper in 0.1 % increments, and the model is solved for each increment until no more fibers or bonds break (at most 15 iterations and continues if the stress is unchanged). A segment of fiber or bond breaks in the model if the projected axial stress of the component is higher than its strength. This iteration continues until the simulated force required to strain the paper falls below 90 % of the highest force observed in the previous iterations. At that point, the paper model is structurally compromised, and the strength is evaluated using the highest observed force in the iterations. The model can continue after this point, but it is stopped at that point for efficiency.

Bending resistance was experimentally evaluated by clamping one end of the sheet and bending the paper 15 degrees after 10 mm. The bending resistance is represented by the forces resulting from this displacement. The measurement/displacement probe is free to move along the paper in the experimental setup. This motion is assumed to be small and not considered for the simulation setup. A 10.2 mm × 4 mm model is created in the simulation. The first 0.1 mm length is clamped similarly to the tensile setup, followed by a 10 mm bending lever, and ends with an additional 0.1 mm after the displacement probe. This extra length provides some additional material after the displacement probe as in the experiment.

The Dirichlet conditions for the displacement of the paper model in the bending resistance simulation occur in a thin domain perpendicular to the length. This Dirichlet domain is centered along the line where the probe in the bending resistance experiment would be. In the experiment, the probe only affects fibers on the bottom side of the paper. In the simulation, this is represented by only adding boundary conditions to the first quarter thickness of the model. The boundary conditions are added by discretizing the edges that pass through this thin domain and explicitly displacing them in the simulation. These Dirichlet conditions are only applied in the coordinate direction of the bend. For an illustration of the bending resistance simulation setup, see Figure 6. Bending resistance is evaluated by solving the equilibrium for this setup and analyzing the resulting forces.

Figure 6

Bending resistance simulation setup.

The model has two resolution parameters. These are element length, the length of the edges of the discretized fibers, and bond delta, the length of the further discretization around bonds. The resolution parameters are chosen by analyzing the tensile stiffness of a 30 g/ m 2 paper model with different element lengths and bond deltas. These results are presented in Figure 7. The discretization parameters used in the results of this article are 0.1 mm element lengths and 20 µm bond delta.

Figure 7

Study of the models resolution parameters (element length and bond delta). The parameters chosen in the model are an element length of 100 µm and 20 µm bond delta.

Figure 8

Tensile forces evaluated from models with different domains strained to 0.5%. For all simulations, the length of one side in the active domain is 4 mm.

Simulations of a 30 g/ m 2 paper for different domain sizes show the tensile forces in the stiffness simulations scale as one would expect of material with constant Young’s modulus. An illustration of the simulated results compared to the results assuming a constant Young’s modulus (fitted) is presented in Figure 8. These results show that 4 mm × 4 mm is a sufficiently large active domain for this specific model, with these parameters, to analyze tensile stiffness.

A similar domain study was performed for the tensile strength simulation as the tensile stiffness simulation. This study analyzes the tensile forces at the break, with the intrinsic strength (fitted) as a reference. The results from the study are presented in Figure 9. These results show that scaling is not as trivial as tensile stiffness, with some domain-related effects. However, considering the difference in strength for the domains analyzed, a 4 mm × 4 mm active domain should give a sufficient representation of the strength of the paper.

Figure 9

Tensile forces at the point at the break for strength simulations of models with different dimensions. In all active domains analyzed, one side is 4 mm.

In the bending resistance experiment, papers 38 mm wide were used and bent with a 10 mm lever. Similar to the tensile cases, a domain study was performed for the bending resistance simulations. Because bending resistance is not an intrinsic property, scaling the forces are not straightforward. However, the bending stiffness of the material can be calculated using the formula presented in (Mark et al. 2002, pp. 233–239):

(4) F = 3 · S · w l 2 · tan ( 15 ∘ ) ,

where l is the length of the lever, w is the width of the sample, and S is the bending stiffness of the material. This formula is derived from linear theory on isotropic sheets and will be used as the reference. The results of this study are presented in Figure 10 and compared to (4) with the bending stiffness estimated from the simulations by taking the average (excluding the smallest domain sizes). These results clearly show that the model scales according to linear theory and that a 4 mm × 4 mm active domain could be sufficient. However, for illustration purposes, the entire lever is simulated.

Figure 10

Forces from bending resistance simulations of models with different widths and lengths. The ideal curves are an equation on the form (4), and in all simulations, one side of the active domain is 4 mm.

The computer used for the simulations is an off-the-shelf desktop PC. It runs CentOS, a Linux distribution, with an Intel i7-7800x CPU and has 64GB of ram. Graphics card acceleration was not used for these results.

Results

Tensile stiffness, tensile strength, and bending resistance are simulated using the proposed paper model with the setups presented for each sample considered in the experiment in both machine direction and cross direction.

A visualization of a strained model in a tensile stiffness simulation is presented in Figure 11. In the figure, each fiber is rendered to its cross-sectional dimensions. Moreover, each fiber is color-coded to illustrate the stresses applied to the fiber in the solution. In the figure, both strained and compressed fibers are observable, showing that the model has a positive Poisson effect, as expected for paper.

Figure 11

Tensile stresses in individual fiber segments from a tensile stiffness simulation.

Figure 12 presents the stresses in the model in one of the bending resistance simulations. Similarly to Figure 11, each fiber is rendered three-dimensionally using its cross-sectional information, and the individual stresses are colored. Three perspectives are presented in the figure: from above, in profile, and from below. The fibers are mostly compressed on the top side, and the fibers are stretched on the bottom as expected. The largest stresses occur where the paper is clamped, with the stresses in the fibers dissipating along the bend profile. Interestingly, it is possible to see some positive Poisson effects in this figure, with oppositely signed stresses present perpendicular to the lever.

Figure 12

Illustration of individual fiber stresses in a bending resistance simulation. The top figure is the model from above, the central figures are the model in profile, and the bottom picture is the model from below.

Stress-strain curves from the strength simulations of a sheet in the machine direction and cross direction are presented in Figure 13. In the figure, the stiffnesses of the sheet in the two directions are illustrated by dashed lines. Some plasticity is observed for both stress-strain curves from bonds breaking throughout the simulation. Paper fibers in these simulations break, but they mostly break around the breaking point of the model. The strains of the model do not represent the strains observed in the experiments, with the strain at break in the experiments being almost twice as large as the simulations.

Figure 13

Stress-strain curves from a strength simulation in both machine direction and cross direction for one of the paper models. The solid lines are the stresses evaluated at each strain, and the dashed lines are the stiffnesses.

Each sheet and property are simulated ten times with different seeds in the randomization, with the mean and standard deviation presented as the simulation result. These results from the simulations are compared to the experimental values in Figure 14. The results show that the model can accurately predict tensile stiffness and tensile strength, with the bending resistance simulations producing stiffer results than the experiments.

Figure 14

Tensile stiffness, tensile strength, and bending resistance validation of the network model. The gray markers are the experimental values, and the colored markers are the simulated results. Each point in the graph represents ten simulations and experiments, with the spread representing the standard deviation.

Figure 15

Numerical metrics for the tensile stiffness and tensile strength simulations.

The tensile stiffness and bending resistance simulations are fast to compute, taking a few minutes on a consumer-grade computer. Tensile strength simulations require a substantially longer time to simulate as it is an iterative process. Hover, these simulations still take less than an hour for the sheets considered. The computational metrics for simulating sheets of different grammage in the tensile stiffness and tensile strength simulations are presented in Figure 15.

Discussion

In the following discussion, the results from the experimental validation of the network model are evaluated and compared to known papermaking identities. First, the stability and scaling of the model are discussed. Then the simulated results for the three structural experiments are compared to constitutive models of the problem. Finally, the approach’s usefulness is discussed, along with potential developments that can be made.

The model is stable and predictable. Stability is observed in Figure 7 for the discretization parameters analyzed. In the domain studies presented in Figures 8 and 10 it is clear that the simulation results for tensile stiffness scales as a material with constant Young’s modulus and the bending resistance simulations scale according to linear theory. The forces at break in the tensile strength simulation get slightly stronger for wider models and slightly weaker for longer models.

The validation results in Figure 14 show that the model can predict tensile stiffness in both machine and cross direction. These simulated results can be compared to the values predicted by Perkins formula for tensile stiffness (Perkins et al. 1978). Here, the same simplification in Rigdahl et al. (1983) is used (assuming the mechanical anisotropy is determined solely by the fiber orientation):

where constant ϕ accounts for the finite fiber length and the coupling between fibers. In Perkins et al. (1978), the following formula for ϕ is presented:

where l is the average length between two bonds’ centers along the fibers and N c is the number of segments before the next segment has the same strain as the network in the direction of the fiber. For the sheets analyzed in this work, N c = 1 (Perkins et al. 1978), resulting in the first term disappearing.

The average length between two bonds, l, is found numerically by analyzing the network models geometry for one seed in each sheet. This analysis of l in a network model is performed by generating a network model that represents the fibers and bonds with only one edge (unlike the structural simulation). These edges are spanned between the two closest points of each two bonded fibers. The average length between bonds, l, is then calculated by taking the average length of the edges representing fibers in the network. This analysis is performed for each sheet for one seed and is used to calculate ϕ for each sheet. On average the l is around 48 µm and ϕ around 0.95. These values are consistent with the sheets’ higher densities in the experiments (Perkins et al. 1978, Rigdahl et al. 1983). The tensile stiffness predictions using Perkins (5) are presented in Figure 16. Comparing the predictions from Perkins with the simulated results in Figure 14 it is clear that the spatial network model is consistent with Perkins’ formula for tensile stiffness for these sheets and parameters.

From Figure 14 it is clear that the network model can predict the tensile strength of the paper modeled. Strength simulations on fibrous materials are challenging. The irregularity of the network structure leads to non-affine deformations (Shahsavari and Picu 2013), with the structure re-organizing at large strains (Onck et al. 2005). However, as paper fails at low strains, these deformations are small, and rearrangement is limited (Deogekar and Picu 2018). Moreover, the domain size required for representable strength simulations varies between publications. Experimental results on paperboards show that the strength of the material is seemingly independent of sample sizes ((5–100) mm × (2.5–100) mm) (Hagman and Nygårds 2012). In contrast, the models that Heyden (2000) presented show a decline in strength for long slim domains, with constant strengths for different widths, similar to the network model presented in this work. Heyden argues that Weibull’s theory could explain this. Vincent et al. (2010) were able to reproduce experimental tensile strengths by analyzing 2.5 mm × 2.5 mm generated paper models geometric properties using a continuum model, which is smaller than the domains used to produce the results here.

Figure 16

Mean tensile stiffness results for each sheet in the experiment (gray) and Perkins’ formula for tensile stiffness (5) in color given the parameters of the model.

The stress-strain curve, Figure 13, shows that the model captures some plasticity. As mentioned, the experimental strains at break are almost twice as large as the strains at break observed in the simulations. This discrepancy could be due to insufficient simulation domains or the models linearity. Other research (Räisänen et al. 1996, Borodulina et al. 2012) has added plasticity to the fiber model. The linear network model could implement plasticity in the fibers by adding an additional threshold to the break condition where the fiber’s stiffness is changed. That said, the model has some plasticity in the stress-strain curve, unlike what is expected from linear models.

The tensile strength results using the spatial network model can be compared to the strengths predicted by the Page equation (Page 1969):

where T is tensile strength, Z is zero-span tensile strength, b shear bond strength, P perimeter of the cross-section, and R . B . A. the relative bonded area of the sheet.

Out of these constants, Z, R . B . A., and b need to be deduced. For b, one of the bond strengths for spruce sulfate fibers presented in Page (1969): 2.9 × 10 7 dynes / cm 2 is used. The relative bonded area, R . B . A., is defined using the following formula proposed by Batchelor and He (2005):

where f is the area of the fiber’s cross-section divided by the area of the smallest bounding box encapsulating the fiber.

The zero span tensile strength, Z, is evaluated using the spatial network model. These results are attained through the tensile strength simulation with a free length of 1 µm in both machine direction and cross direction for each sheet in the validation. Each setup is simulated ten times, with ten different seeds, and the averages determine the Zs for each sheet. With the Zs determined, tensile strengths predicted by Page (6) are calculated and presented along with the tensile strengths from the experiments in Figure 17.

Figure 17

Mean tensile strength for each sheet in the experiment (gray), compared to the tensile strength derived through Page (6) (color).

The results in Figure 17 show that the strengths predicted by Page (6) are similar but have less spread in the machine direction and the cross directions compared to the experiments. This spread is captured in the simulated results using the spatial network model.

The results from the bending resistance simulations are stiffer than the experimental values in Figure 14. The stiffer results are not surprising from a linear model, as the experiment has non-linear components. This increased resistance could be due to non-linearity or the modeling approach. The non-linearity hypothesis can be evaluated by comparing the results to a linear continuum model of the problem.

The bending resistance experiment can be modeled using the formula presented in (Mark et al. 2002, pp. 233–239), where the paper is assumed to be isotropic with constant Young’s modulus. Using Young’s modulus and thickness from the experiments, the bending stiffness in the machine direction and cross direction can be formulated as follows:

where t is the thickness of the sheet. With (4), bending resistances can be calculated using these bending stiffnesses. These predictions are presented and compared to the spatial network models bending resistances in Figure 18.

Figure 18

Bending resistance calculated assuming the paper is an isotropic sheet ((4) and (7)) with Young’s modulus equal to the average experimental results (gray), compared to the simulated bending resistance (color).

In Figure 18, it is clear that the model acts almost identically to the linear continuum model. The bending stiffness of the paper model is not an explicit parameter of the model but is derived entirely from the microstructure. This means that the deviation in the validation in Figure 14 might be fixed if a non-linear version of the model was used.

The results produced in this paper are performed on consumer hardware in a reasonable time. In Figure 15 the numerical metrics are presented. Tensile stiffness simulations can be performed and analyzed in about a minute, and strength simulations in less than an hour. Bending resistance simulations can also be performed in minutes. The strength simulations can be optimized. In the strength simulations, each strain is evaluated multiple times with increments of 0.1 %. If tensile strength is the only result that is of interest, some intermediate strains can be skipped, with only the breaking point analyzed in detail. Moreover, due to the iterative nature of the simulation, a non-linear version of the model might have similar execution times. In Heyden (2000), it was concluded that domains 1.2 times the longest fiber is sufficient for tensile stiffness simulations. With smaller domains, far denser paper can be analyzed without running into numerical constraints.

The simulation times are suitable for product development, where iterations and new ideas might be evaluated through simulation for numerical validation. The simulation framework enables quick changes of detailed parameters such as individual fiber properties, pulp distribution, and bond mechanics, where resulting changes in sheet mechanics can be attained and analyzed. These kinds of product development tools are necessary for an industry that significantly needs to reduce material use. Moreover, due to the decline in biodiversity, returning to more heterogeneous forests will most likely be necessary, promoting hardwood trees and introducing fiber types with significantly different properties compared to fir and pine. Simulations will be the faster method of understanding the mechanical properties of sheets based on new recipes.

Moving forward, several additions to these network models should be analyzed. It was observed in Görtz et al. (2019) that the LOD Method, a multiscale method, can handle these network models optimally. The method could enable large-scale simulations with periodicity, say 4 mm × 4 mm periods. Comparing the model with a non-linearized version might further justify a linear model for tensile stiffness simulations, and the bending resistance discrepancy could be resolved. In the strength tests, adding plasticity to fibers and adding time dependence might produce accurate stress/strain curves. Straight fibers are a good starting point, but adding shapes and using shape factor measurements might give insights into their effect on the structural properties of paper. The ISOP project developed a framework for simulating paper forming (Svenning et al. 2012). Using these simulated models geometries in the modeling step would provide a coherent simulation methodology of the complete paper manufacturing process.

Conclusions

A network-based linear model was presented to simulate the mechanics of different papers. This model adds cross-section information analytically, and most of its parameters were motivated by experimental data and published results in the literature. The parameters not motivated by experimental data have to do with bonding, and those parameters were chosen to scale appropriately. The network model can predict tensile stiffness, tensile strength, and bending resistance for both machine direction and cross direction for the different low-grammage papers. These results are consistent with Perkins’ formula for tensile stiffness, the Page equation for tensile strength, and linear elastic theory for bending resistance. The model presented acts as a linear elastic material when bent. The bending resistance simulations produce stiffer results than the experiments. This discrepancy might be solved by adding non-linearity to the model.