1 Introduction

There is a large potential for wood-fiber based materials such as paper and paperboard in numerous engineering applications. Classical applications of these materials are for example in the packaging of products such as food packaging, which is an important and still growing industry [196]. However, because of sustainability reasons in terms of renewability and recyclability, paper and board are also gaining importance in other applications as a substitute for e.g. building materials.

In order to optimize the design of both the materials and the products made of these materials, computational methods have been utilized more and more in the last decades. The major challenges to be coped with are the multi-physics and the multi-scale nature of paper. Even though the effects and phenomena occuring on the different scales are well-known and well described, see e.g. the general overview on paper physics provided in [4], it is still a field of ongoing research to account for them in numerical models.

There exist recent reviews that deal with specific aspects of the modeling strategies for paper and paperboard. In particular, in [77] recent applications of the finite element method (FEM) to paper, paperboard, and corrugated paperboards used in food packaging are summarized. However, focus is mostly on basic findings on the sheet level and structural scale whereas the multi-scale and multi-physics issues are not discussed. Another exhaustive review has been provided in [198] with the scope on forming processes of paper. There, three-dimensional deformations and damage mechanisms were recapitulated especially for deep-drawing, press-forming, hydro-forming, and others. These processes and the corresponding effects are comprehensively discussed. Nevertheless, also in this work the multi-scale character of the material is not addressed.

Thus, the aim of the current work is to sum up recent trends in the modeling of paper and board on the different scales of interest which are fiber, network, sheet, and laminate scales. In addition, the purpose of this systematic review is to discuss both achievements as well as drawbacks and open questions.

2 Modeling Paper on the Fiber Scale

On the fiber scale, the key aspects in modeling deal with the single fiber behavior and the bonds between individual fibers.

2.1 Single Fiber Behavior

While it is well-known that the single cellulose fibers consist of several different layers, in most approaches to model their material response only the so-called S2-layer, which contributes most to the mechanical behavior, is considered. This layer is composed of a cellulose/hemicellulose/lignin matrix reinforced by microfibrils, see e.g. [52, 53].

The longitudinal material response in fiber direction strongly depends on the microfibril angle (MFA), which can vary significantly for different fibers [248]. Recent experimental methods to measure the mechanical behavior of fibers including their tensile strength are based on single fiber tensile tests, as shown in [159] for softwood and in [120] for hardwood fibers. In addition, the influence of humidity [121] and of refining as well as small-scale fiber deformations [132] on the single fiber strength have been investigated.

Even so, it is very challenging to measure the mechanical response of single fibers in the nonlinear regime. Thus, in most applications with consideration of the MFA, fibers are modelled linear elastic. Nevertheless, approaches to include the elasto-plastic behavior have been presented e.g. in [203], where the anisotropic yield criterion of Hill has been applied, and in [33], where nonlinear kinematic hardening with tension-compressoin asymetry has been taken into account for cyclic loading conditions. Further, the Tsai-Hill failure criterion has been used in [184] in an attempt to predict probable locations of failure initiation in wood fibers.

Another important issue is the geometrical description of the fibers including the fiber length, orientation, and thickness distributions. Furthermore, measuring the cross-sections accurately is important for a realistic consideration of single fibers. Nowadays, these geometrical data can be extracted from micro-CT images, see [160, 223] among many others.

2.2 Fiber-Fiber Bonds

In order to properly describe the mechanical behavior of fiber networks, accounting for the fiber-fiber-bonds is crucial [242]. Here, one challenge is determining the actual area of the bonds, which in many cases is much smaller than the analytical contact area [124,125,126]. Even more, it has been shown in [175] that only 30–40% of the measured overlap area is actually in bonding contact; according to [213], this value varies even in the range from 15 to 88%. However, nowadays, the bonding area is measured through computational image analysis of 3D X-ray microscopy [174, 247] or confocal laser scanning microscopy [150]. In addition, investigations have been made concerning the area in molecular contact in fiber bonds [105] and the corresponding energy contributions of bonding mechanisms [103, 104].

Another even more important aspect is the measurement of force-displacement curves for 2-fiber systems which are usually utilized for investigating the bonding strength. For this, two crossing fibers are considered, one of which is fixed on both ends while the other one is loaded. Different loading cases have been investigated in the literature in order to describe general loading scenarios (see Fig. 1):

  1. 1.

    peeling test (\(\approx\) mode I): one fiber is fixed at both ends, the other one is loaded perpendicular to the fiber direction at one end [168, 169, 220];

  2. 2.

    shearing test (\(\approx\) mode II): one fiber is fixed at both ends, the other one is loaded in fiber direction at one end [120, 168, 169, 217] (and [121] including the influence of relative humidity);

  3. 3.

    tearing test (\(\approx\) mode III): one fiber is fixed at one end, the other one is loaded in fiber direction at one end [81, 169];

  4. 4.

    z-directional test (\(\approx\) mode I): one fiber is fixed at both ends, the other one is loaded perpendicular to the fiber direction at both ends [139, 218].

Fig. 1
figure 1

Different loading cases applied in the fiber-cross experiment

Full size image

First finite element simulations of the fiber-cross can be found in [241], where both fibers have been modeled orthotropic elastic with simplified cross-sections without lumen. There, comparison between solid and beam elements has shown that the 3D solid elements can provide much more accurate results. Further, first qualitative results for normal traction and shear loadings have been shown in [166]. Based on these investigations, peeling and shearing loadings could be analyzed by using shell elements accounting for the curvature of the fibers [167]. A similar modeling approach has been applied in [71], where resultant forces and bending moments have been determined from the peak stresses in all three modes (peeling, shearing, and tearing). Finally, the influence of geometry changes—such as the fiber cross angle and the MFA—on the stiffness of fiber-fiber bonds has been investigated in [39].

In all the works mentioned above, the interface between the fibers has been modeled as perfectly bonded. In contrast, in recent works, the debonding of the fibers is modeled additionally by applying cohesive zone elements [85, 165]. Thereby, the damage progression within the interface can be evaluated until failure of the bond.

2.3 Open Problems on the Fiber Scale

As already mentioned above, the nonlinear single fiber behavior can be taken into account in numerical simulations, as shown e.g. in [33, 203] for elasto-plasticity and in [184] for failure analysis. Nonetheless, precise stress-strain data, that is needed to calibrate such models, is hardly available, yet. This holds even more when the dependence on loading rate, moisture, and temperature is to be considered in addition.

Moreover, obtaining experimental results on the fiber behavior in transverse direction is still an open topic. While the material response in fiber direction can be investigated via uniaxial tension tests, the transversal direction is much harder to test. In principle, the latter can be achieved for example by using AFM nanoindentation for elastic [82] and viscoelastic [62] properties. However, improvements of these methods are still necessary if results shall be used for accurately calibrating the numerical models.

Another critical aspect is the characterization of the cohesive zone models used to describe the debonding of the contact area between fibers. The numerical investigations mentioned above are able to provide an insight into the behavior of the interface and the influence of geometry variations. Nevertheless, the calibration of cohesive zone model parameters—in particular the determination of the critical fracture energies \(G_{Ic}\) and \(G_{I\!Ic}\)—is still problematic, since the fiber cross experiments are not pure mode tests [3].

3 Modeling Paper on the Fiber Network Scale

As shown in [248], the extensibility of paper relies mostly on three major factors: deformability of the single fibres, the ability of fibres to form strong and flexible bonds, and the three-dimensional structure of fibre network created during the manufacturing process. Hence, models are also needed to characterize the structure of fiber networks.

The first network model for paper has been introduced already in 1952 [60]. Since then, numerous different network models have been developed and applied to paper. These can be classified in several categories: (1) fiber networks that are reconstructed from micro-CT images, (2) structured networks, (3) synthetic networks with random distributions, (4) synthetic networks with statistical distributions.

3.1 Fiber Networks Reconstructed from Micro-CT

One way to generate fiber networks is the reconstruction of real network structures from micro-CT data. In particular, the fiber data as well as the network connectivity data can be analyzed from microtomography images as shown in [35, 251]. Alternatively, in order to avoid the separation of single fibers from the micro-CT images, local and global fiber orientation tensors can be calculated for porous fiber networks given in terms of binarized voxel images. This approach has been employed successfully for thermal conductivity simulations in [221]. Also, a method for determining the number of contacts in a fibrous network based on shortest-path measurements has been developed in [73]. Furthermore, 3D synchrotron X-ray microtomography has been applied to characterize the paper structure of z-structured paper by introducing micro nanofibrillated cellulose in [46]. Finally, the combination of results from synchrotron X-ray 3D microtomography with Mercury intrusion porosimetry (MIP) data has been shown to further improve the information about the considered structures [47].

Clearly, if the network structure is obtained from reconstructing real fiber structures, the results are very close to reality. On the other hand, only very specific structures can be studied without the possibility of varying parameters in a stochastic manner. Thus, it can also be advantageous to generate synthetic network structures for numerical investigations.

3.2 Structured Fiber Networks

One way to generate artificial fiber networks is to use simplified structural designs. For example, a 2D lattice model for simulating the failure of paper has been introduced in [157]. An extension of this lattice model accounting for bond elements and hence three-dimensional effects has been given in [158].

Two-dimensional lattice models based on a periodic (triangular) distribution of spring elements has also been extended already such that interfiber bond failure and subsequent frictional fiber sliding can be evaluated [252]. Such models can also be used as basis of multiscale investigations by simultaneously using the quasicontinuum method, as shown in [23].

Additionally, damage propagation within lattice structures has been investigated based on various criteria of the elimination of overstressed beams. Particularly, a comparison of damage patterns in triangular stretch-dominated and hexagonal bending-dominated lattices can be found in [50], whereas hybrid structures combining both triangular and hexagonal parts are given in [51].

Finally, the elastic-plastic response of 3D fiber networks has been investigated in [164]. While for the analysis of the effect of relative density on the uniaxial yield strength a random fiber network with periodic boundary conditions has been utilized, the dimensional analysis has been performed on idealized transversely isotropic models.

Structured fiber networks are advantageous because they are easy to generate and the computations are relatively cheap. On the other hand, structural and dimensional effects can only be approximated. Therefore, the evaluation of more complex structures can be necessary if structural effects shall be evaluated.

3.3 Random Fiber Networks

In order to allow for statistical evaluations, random fiber networks (RFN) have been generated in multiple ways such that the location of fibers in the network is random. Certain aspects of such random fiber networks can be analyzed by theoretical methods. For example, in [72], a theory has been developed, that describes the statistical properties of rod packings, while taking into account that the deposited rods cannot overlap and thus induce steric hindrances. Further, a new buckling theory including a statistical distribution of free-span lengths has been proposed in [130] and tested against experimental data.

Even so, for most investigations of practical relevance, numerical models need to be established to represent the mechanical behavior of RFN. Since normal paper sheets are mostly formed as highly oriented 2D networks (see e.g. [7]), several attempts have been made to model two-dimensional RFNs. For example, the influence of network deformation and fiber bond fracture on the macroscopic degradation and failure of paper materials has been investigated in [98, 117]. More recently, the mechanical behavior of cross-linked random fiber networks with inter-fiber adhesion have been investigated in tension and compression in [186], where the periodic networks have been created by performing Delaunay and Voronoi tessellations.

However, in order to cover the effects of fiber deformations and fiber bonding, 3D models have proven to be more accurate in many cases. Thus, 3D random fiber networks have been constructed already in [137] to simulate deformation and failure behavior of networks with dynamic bonding/debonding properties using the FEM. Here, the single fibers have been modeled via Tymoshenko beam elements. This approach has been adopted in several recent works. For instance, rotational constraints for 3D beams within random networks have been considered in [179]. Further, accounting for the debonding of fiber-fiber bonds has been performed by use of a cohesive zone model with linear traction-separation law in a penalty-based contact element [34]. The latter formulation has also been applied to examine the effect of changing paramters for the network structure on the short span compression strength [40]. Similarly, the effect of irreversible failure of fibres and of inter-fiber bonds with varying inter-fiber bond density and bond strength has been investigated in [86] for cellulose nanopapers. What is more, in [65], damage accumulation and failure in fiber networks have been examined based on 3D random networks defined on a Voronoi structure, where fibers have been modeled as Timoshenko beams of circular cross-section, while the bonds have been modeled as uncoupled springs with translational and rotational stiffness.

An alternative approach to modeling RFN makes use of the discrete element method (DEM), where each end of fibers as well as each fiber bond are represented by discrete particles, as shown first in [204] for computing the dynamic fracture of thin network materials. Also using the DEM through a series of inter-connected particles in order to generate curly, kinky, and twisted fibers, the creping behavior of soft papers has been explored in [106]. Later, the same approach could be used also to research the uniaxial compression behavior [107]. Lastly, literature results concerning the scaling of elastic modulus as well as strength with increasing density could be reconfirmed in [26] through DEM simulations.

3.4 Statistically Oriented Fiber Networks

If fiber networks are created randomly, the actual distributions of geometrical characteristics are not accounted for. However, these distributions can be incorporated through statistical studies. For example, the fiber lenth distribution has been included in a statistical sense in [140].

Even more important is regarding statistical data on the fiber orientation distribution as shown e.g. in [227] for medium density fiberboards (MDF) made from wood fibers. Extensions of this methodology allowed for the inspection of multilayered structures [226] and of the viscoelastic material behavior of resin-bonded nonwovens [233]. Moreover, the influence of geometrical and spatial effects including the orientation of fibers in the network has been researched in [128], where focus is layed upon the realistic description of crossing fibers. Last but not least, orientation of fibres and positions of bond points in fabrics exposed to stretching in two main orthogonal directions (MD and CD) have been simulated and validated by experiments in [230].

3.5 Open Problems on Fiber Networks

As shown above, random fiber networks for paper have been extensively studied already (see Sect. 3.3). These works are mostly based on DEM or FEM with beam elements, which is computationally very efficient. However, it is not clear how realistic the results can be for out-of-plane loading scenarios. Therefore, comparisons to computations with solid elements as presented in in [140, 149] are advisable.

However, for a more realistic representation of the network behavior, statistical distributions should be considered in addition. Although there are several works dealing with the influence of fiber orientation distributions (see Sect. 3.4) on the overall network response, the majority of these works does not bear on paper but other fibrous networks such as MDF, resin-bonded nonwovens, and nonwovens composed of polypropylene fibers. In consequence, this method needs to be adopted also in the modeling of paper and paperboard with the corresponding statistical distributions and effects. One approach going into this direction has been published recently in [209], where a scheme to identify material parameter distributions using only a limited number of fibres has been developed, which is based on probability density functions (PDF) and Bayesian inference.

Further, it has been shown recently in [136], that local mass density variations play a crucial role, particularly for the fracture behavior of low density papers. Hence, these variations need to be included into the network models additionally.

In addition, the influence of moisture and temperature need to be implemented into the networks models. There exist models for describing the mechano-sorptive creep as result of moisture cycling [237, 238], dimensional instability of networks such as twist due to moisture change [222], as well as for hygro-expansion of the fibers [36,37,38, 180]. These models are included into random fiber networks. Nevertheless, these are based on idealized structural representations of the bonding between fibers, and thus more sophisticated strategies will be required in the long term.

Another critical issue is the homogenization of the network response to achieve effective properties for the macroscale. The fiber networks considered here cannot be periodized, present an infinite contrast of properties, and include an interconnected porous phase. As shown in [67], that leads to very large (’gigantic’) RVE-sizes if only one single RVE is to be investigated. A reduction of the RVE-size is possible if multiple realizations are considered. However, the smallest size of the volume element on the network level that needs to be considered such that the homogenized response is representative in a statistical manner and such that boundary conditions do not affect the result is still not known for paper. While such studies exist for other materials—for instance, for short fiber reinforced elastomers, where at least 4000 fibers should be considered in order achieve a converged RVE-size [48]—they are missing yet for paper and paperboard. First attempts to achieve a converged RVE-size based on numerous computations of multiple realizations have been presented in [8] and [149]. However, no real material parameters for the fibers’ nonlinear behavior or the fiber bonds have been included in these works. This is still an open topic for future research.

4 Modeling Paper on the Sheet Scale

On the microscale, paper is a composite that consists of fibers in a network and sourrounding air (see Sect. 3) or can be treated as single pores surrounded by cellulose [214]. In contrast, on the sheet scale, paper can be considered as a homogeneous, anisotropic material that is described in the framework of continuum mechanics. To account for the anisotropic nature, one usually defines three principal directions: machine direction (MD), cross direction (CD), and thickness direction (ZD), see Fig. 2.

Fig. 2
figure 2

Definition of principal directions in paper

Full size image

Nowadays, such anisotropy is often accounted for by use of the concept of structural tensors [231], which can be applied to problems with large deformations of nonlinear elasticity (see e.g. [210] among many others) and (hyper-)elastic-plasticity [96]. Alternatively, anisotropic 3D material models can be obtained from 1D material models accounting for orientation distribution functions by applying the unit sphere approach suggested in [178], which has already been adapted successfully to model composites with consideration of fiber misalignment [145]. Instead, stochastic continuum models can also be constructed by applying Gaussian fields to the 1D formulation [172]. However, independent of the method utilized, it is crucial to account for the anisotropy of paper on the sheet scale.

4.1 Elastic-Plastic Behavior of Paper Sheets

In most models, the assumption is made that the in-plane and out-of-plane behavior can be considered separately; this is justified by the finding that the out-of-plane Poisson’s ratio is close to zero, as shown in [234]. Therefore, most elasto-plastic material models for paper are suitable for either the in-plane or the out-of-plane behavior.

An overview of the models described in the following is provided in Table 1.

Table 1 Overview of elasto-plasticity models for paper on the sheet scale
Full size table

4.1.1 In-Plane Elasto-Plasticity of Paper Sheets

For the in-plane elastic-plastic behavior of paper sheets, the assumption of small deformations is justified in many technical applications. Then, the basic assumption of an additive split of the total strain into an elastic and a plastic part, \(\varvec{\varepsilon }=\varvec{\varepsilon }_e+\varvec{\varepsilon }_p\), can be used. However, if finite deformations (that is finite strains and/or finite rotations) have to be taken into consideration, a multiplicative split of the deformation gradient into an elastic and a plastic part, \({\mathbf{F}} = {\mathbf{F}}_{e} {\mathbf{F}}_{p}\).

One way to account for the material anisotropy is applying the concept of an isotropic plasticity equivalent (IPE) material, where the IPE-material is a fictitious isotropic material, subjected to a stress state that equals the corresponding stress state in the actual anisotropic material. This concept, initially introduced in [127] for pressure independent materials, has been applied to paper e.g. in [170]. Recently, this model has been shown to be suitable for predicting results for the box compression test in [163]. In the IPE-model, as in most elasto-plasticity models that have been formulated for pressure independent materials such as metals (see e.g. [42] for an overview), the yield function is given in dependence of the deviatoric stress. Since this might no be realistic for paper, numerous alternative models have been invented specifically for paper.

Another way to include anisotropy into the elastic-plastic model is the use of an orthotropic yield function that accounts for different yield stress values in different preferential directions. The most common one is the formulation of Hill [102], that includes six yield stress values. In particular in forming processes such as creasing and folding, this model has been applied successfully e.g. in [90, 109, 110] amongst others (see also Sect. 7). Additionally, Hill’s yield function has been utilized for the simulation of the box compression test [254] and of lateral compression tests of layered paperboard tubes [61].

Further, a well-known model for the in-plane elastic-plastic material behavior of paper is the one presented by Xia et al. in [253]. Here, the assumption is made that the out-of-plane material response is purely elatic, such that the plasticity model is not affected by any stress component that refers to ZD. Then, the yield surface can be defined as the sum of six sub-surfaces, four of which describe the tension and compression in MD and CD, respectively, whereas the remaining two surfaces correspond to positive and negative shear, respectively. As shown in [192, 193], this model captures the in-plane response of paper in creasing and folding very well, while the out-of-plane behavior should be modified to also account for delamination in thickness direction. Noteworthy, the formulation of [253] and [193] has already been implemented into LS-DYNA as \(^*\)MAT_PAPER as material type 274 [219].

A modification of the model of Xia et al. [253] has been presented in [30]. There, the yield function has been taken over from the original model but the harding part has been modified such that a coupling between the hardening of the six different sub-yield surfaces is introduced, allowing for non-isotropic hardening. Based on this work, in [31], a similar formulation has been developed that included hardening in tension and shear, while ideal plasticity was assumed in compression. Furthermore, in order to consider also viscous effects in model, the formulation of Xia et al. has been extended in [240] accounting for strain-rate dependent material behavior. Even more, also the densification in thickness direction has been incorporated (see Sect. 4.1.2) and the hardening model has been extended to be anisotropic-kinematic without introducing additional parameters. Another variant of the Xia et al. model has been provided in [146], where—as in the latter one—a small strain formulation has been derived including both isotropic and kinematic hardening in order to account for the anisotropic nature of hardening behavior. In addition, an alternative small-strain version of the Xia et al. model has been suggested recently in [205], where focus is laid on introducing a non-associated flow rule.

Moreover, the yield criterion can be formulated based on the Tsai-Wu criterion, which consists of a linear combination of quadratic and linear function in stress. This criterion has been applied in [97] for modeling the in-plane behavior of paper successfully. Moreover, based on a previous investigation [96], the flow rule has been formulated in a non-associated manner.

Finally, it is also possible to account for the material’s anisotropy by defining several different yield functions in the framework of multi-surface plasticity. For example, in [250], four different yield functions have been defined depending on whether MD or CD loading is dominant and discriminating between tension and compression. The formulation is simple, only requires a relatively low number of parameters, and includes even thermal dependence of the mechanical parameters. Another model based on multi-surface plasticity has been recently derived in [20], where three distinct yield functions are applied to distinguish between normal loading cases in MD and CD, respectively, as well as shear loading. It has been pointed out that the results are similar to the ones obtain from a model with one single non-quadratic yield function.

4.1.2 Out-of-Plane Elasto-Plasticity of Paper Sheets

In many technical applications, the out-of-plane response of paper is relevant, too. Consequently, several models have been invented also for elasto-plasticity in thickness direction.

For instance, based on experimental findings [235], a small-strain model has been suggested in [236], in which the nonlinear elastic behavior has been treated by accounting for the porous nature of the material as well as the plastic response under high compressive loads in ZD direction through defining a bounding surface.

The different material behavior in compression and shear can also be tackled by the concept of multi-surface plasticity. In particular, two separate yield functions with corresponding separate hardening laws have been introduced in [188]. One of the two functions represents compressive stresses in ZD while the other refers to shear.

To this author’s best knowledge, there are only two publications available in literature, in which in-plane elasto-plasticity models have been extended such that also the out-of-plane direction is treated. The first one is [31], in which the yield function of Xia et al. is augmented by six additional sub-surfaces to include all ZD-loading components. The model is based on the assumption that the fiber-layer normal direction remains unchanged during out-of-plane shear deformations. Hardening is only active for ZD compression, whereas all other out-of-plane plastic deformations are modeled as ideal plastic. Similarly, also in the second work [240], the yield function of Xia et al. is expanded such that in total 14 sub-surfaces are considered. An important aspect in this model is accounting for the material through-thickness densification in both the elastic and the plastic parts of the model.

Finally, the densification effect has also been addressed in [147]. This formulation is thermodynamically consistent and valid for finite deformations. It includes three sub-surfaces in the yield function. As in the model of [31], hardening only occurs for ZD compression. Also, the internal friction effect is taken into consideration such that the plastic response in compression and shear is coupled.

4.2 Modeling Viscous Effects in Paper Sheets

It is a well-known fact that paper shows rate-dependent material behavior and viscous effects [95] because of the hygroscopic nature of the fibers (not the bonds [64]). Particularly the creep and relaxation behavior are of high significance and have thus been subject to extensive studies—see e.g. the overview provided in [55].

In the last years, experimental studies on the creep behavior of paper have mostly focused on its influence on the compressive strength of corrugated paperboard packages. For example, in [78], box compression tests as well as creep compression tests have been performed. Amongst others, it has been found that refrigerated conditions increased the creep rate of the packages in comparison with the standard conditions. In addition, the mechanosorptive effect, which describes the combined action of mechanical creep loading and changing moisture content, has been included into testing procedures in [114] and [134, 135]. In these works, testing has been executed in climate chambers that allowed for controlled cyclic changes of relative humidity leading to accelerated creep. Concerning relaxation, the effects of strain rate on the tension-strain curve and relaxation of wet paper have been researched in [133]. Further, both the creep [183] and relaxation [182] behavior have been studied recently with the help of a folding experiment that has been performed on creased paperboards.

One way to model the mechanosorptive effect is to consider the stresses that are produced by moisture content changes in the constraint fibers of the underlying network [6, 237, 238]. For instance, in [6], from the relations derived on the network level, a rheological model with springs and dampers has been developed that allows the description of creep via a continuum model. However, no structural problem has been considered.

There are only few works in which the viscous effects have been tackled also on the structural level by means of finite element simulations. In [207], this is performed for corrugated paperboard in two steps: The transient analysis has been conducted first in order to calculate the relative humidity values as a function of time, and the second stage of the analysis has been calculating the fiberboard deformation response (hygroexpansion) due to the change in relative humidity values. Creep is incorporated through a simple creep equation that relates the strain to stress and time.

Besides, the effect of the through-thickness moisture content gradient on the moisture accelerated creep has been simulated by using an isotropic hygro-viscoelastic model in [229]. It has been shown that the rate of the moisture content change and the steepness of the gradient affect the creep rate greatly. In addition, based on these results it has been suggested that the internal stresses generated by hygroexpansion may increase the creep rate when material layers with different moisture content exist.

Alternatively, the relaxation behavior can also be included to the model by directly affecting the internal stress decrease. This has been done in [75, 142] in combination with a simple plasticity model with von Mises yield criterion for plane stress states to capture cyclic and irreversible straining phenomena.

Furthermore, the multi-surface plasticity of Xia et al. (2002) [253] (see Sect. 4.1) has been extended to viscoelastic and viscoplastic behavior in [240]. For the viscoelastic model, the generalized Maxwell model has been adopted, such that the stress is split additively into the steady-state and the transient parts, where the transient stresses are computed from simplified integral expression accounting for the elastic strain rate. For the viscoplastic model, in order to account for time- or rate-dependent plastic deformation, a Perzyna-type power-law kinetic relation is introduced for the viscoplastic strain evolution. Applications of this viscoelastic viscoplastic model to compressive creep and compressive relaxation in thickness direction can be found in [84].

4.3 Modeling Effects of Temperature and Moisture in Paper Sheets

Table 2 Overview of time-/rate-, moisture-, and temperature-dependent models for paper on the sheet implemented into mechanical structural analyses (FEA)
Full size table

In the previous section, the influence of moisture on the material behavior of paper has partly been described already, because it has a significant effect on the creep behavior. However, the dependence on moisture and moisture change has been studied as a separate topic in several works. An overview about the factors affecting the hygroexpansion of paper is provided in [151], where mostly the fiber and the network level are addressed. However, numerical models are not discussed in detail at all.

The moisture transport in cellulose-based materials can be modeled using the framework of porous materials. For instance, a theoretical water vapor transport model has been developed for cellulose-based materials in [19] based on the assumption of isotropic material properties. Further, a steady-state finite element model of moisture transport in corrugated fibreboard has been presented in [41], where values of moisture permeability have been determined experimentally. Finally, also on the sheet level the moisture adsorption can be described, as has been conducted in [27] using a computational fluid dynamics (CFD) model. However, in the latter models, no coupling to the mechanical behavior has been conducted.

One way to do account for the moisture dependence of the mechanical material response on the sheet level is through homogenization of network simulations [37, 38, 246]. The corresponding fiber networks have been generated randomly (see Sect. 3.3), whereas idealized geometries have been studied for the hygroexpansion. Another strategy to investigate the effect of moisture diffusion on the microstructure has been presented recently in [123]. The model established therein is based on 2D scanning electron microscope images of real microstructures for describing moisture transport in paper sheets in order to explain the time-dependent natural aging process of a paper sheet and book stacks.

As an alternative, the moisture-dependence can be implemented directly into the continuum model approach on the sheet scale in a phenomenological sense. For example, the moisture induced out-of-plane deformation of a paper sheet has been investigated in [143], where elastic behavior assumed and anisotropy is taken into account. Moreover, in a similar approach to the already mentioned works [75, 142] (see Sect. 4.2), an elastio-plasticity model has been adopted in [156] for the plane stress case, where the yield point as well as the hardening have been formulated as functions of both the anisotropy and the dry solids content.

The other very important factor influencing the mechanical behavior of paper is temperature. One thermo-mechanical simulation has been undertaken in [250] for the deep-drawing process of paperboard. In particular, one of the relevant process parameters to be varied has been the die temperature. To include the temperature dependence into the deep drawing simulation, the mechanical parameters have been treated as functions of temperature.

However, in many applications with high temperature variations, also the moisture change is not negligible. Then, the combined effects of moisture and temperature on the mechanical response of paper have to be captured [152]. This is true e.g. for printing processes of paper, which are therefore subject to several examinations. For example, the coupled heat and moisture transport problem has been tackled in [256], where the considered paper sheet has been moving over a warm print surface. Likewise, in [141], a mathematical model of moisture and heat transport has been derived and applied to a paper sheet that moves through a fusing nip in a printer.

Another process, in which both moisture and temperature play crucial roles, is hydroforming of paper. That process has been modeled in [154] using the simple elasto-plasticity model derived in [250] mentioned above. In this work, the material parameters are defined as functions of moisture as well as temperature. However, the problem has not been solved in a completely coupled manner, but instead some key assumptions have been made, which allow the drying effect to be captured within the framework of the material model without the need of an external solver for the moisture gradient.

Additionally, the heat and moisture transport can be relevant during storage of paper rolls, which has been treated in [5]. In this work, a macroscopic three-phase model—consisting of solid fiber, bound water, and pore gas—has been derived to describe transport and mass exchange processes in paperboard. Governing equations have been obtained through simplification of the general field equations: mass conservation, linear momentum, and energy conservation, established in the hybrid mixture theory framework. Similarly, a three-phase model has also been introduced in [10], where a two-scale framework has been employed. The resulting model has been applied [12] to failure analysis of moist packaging material exposed to excessive heating and in [11]—in combination with an anisotropic elasto-plasticity model taken from [31]—to investigate the response of paperboard in conditions similar to those present during a local sealing of two sheets of paperboard.

Concluding, an overview of the models that have been implemented into mechanical structural analyses (FEA) is provided in Table 2.

4.4 Damage Progression in Paper Sheets

Predicting the failure of paper and paperboard is relevant for improving the material performance. Here, the effect of different loading scenarios is of particular interest. These can be examined experimentally for example by using inflation tests [29] or biaxial in-plane tensile testing [153]. In both of these works, the experimental results on the failure of paper have been transfered into mathematical models by defining failure surfaces, where in [153] both stress-based and strain-based surfaces have been analyzed. However, even before the final failure state is reached, significant damage progression can be observed in paper and paperboard. This evolving damage leads to reduction of the stiffness, as shown experimentally e.g. in [99].

Recent experimental investigations focus on local effects on the damage behavior. For example, in [138], local structural properties have been measured and put in relation to local tensile deformations, which due to accumulation may yield to initializing the failure of paper. Further, the effect of small-scale deformations caused by mechanical treatment (refining) of pulps on the tensile behavior of fibers and papers has been studied in [132].

From modeling point of view, it is also possible to describe the damage progression on the fiber or network level in order to account for localized effects. One way to do that is establishing a relation between the fraction of broken bonds on the network scale—which represents the degree of material degradation—and the reduction of the network’s macroscopic stiffness [98, 116]. More recently, computations of deformation and damage in random fiber networks (see also Sect. 3.3) can be found, for instance, in [230].

Nevertheless, it is more convenient to model damage evolution in a continuum sense on the sheet scale. For this, the classical approach consists of introducing a damage variable D, which controls the material degradation in a continuous fashion. Since this damage variable is treated as an internal variable, an evolution equation needs to be defined in order to describe how the damage evolves under the current loading path. Such evolution equations can be defined either phenomenologically—as done in [115], where the damage evolution law is defined in dependence of the damage energy release rate—or based on microstructural observations on the network level [99].

One important issue is the modeling of continuum damage in an anisotropic manner, since the material response of paper is strongly anisotropic. One way to achieve an anisotropic formulation is the introduction of more than one scalar-valued damage variables (see e.g. [225] for details concerning composites in general). Particularly for paper, three different damage variables can be introduced in order to discriminate the behavior in MD, CD, and ZD [118]. Further, in [49], several differents fiber sets have been considered on the network level, where for each fiber set a scalar damage variable has been defined. As an alternative, a second order damage tensor can be introduced such that a tensor-valued evoltion law needs to be defined [206].

Finally, it should be mentioned that continuum damage models generally suffer from pathelogical mesh dependence. In order to obtain mesh-independent results, viscous regularization of Perzyna-type can be applied, as shown in [206]. Alternatively, non-local models, in which non-local quantities are incorporated into constitutive equations can preserve mesh-independence [49, 99, 118].

4.5 Open Problems on Sheet Level

Concerning the elasto-plasticity models, one major challenge is the combination of in-plane and out-of-plane material behavior. Even though there are models that take both into account [2, 31, 240], not all physical effects are captured by these models, yet. Moreover, one particular challenge is the determination of material parameters. As can be seen from Table 1, the more effects the models can cover, the more parameters are usually involved that need to be calibrated. Even though several strategies for parameter identification exist, see e.g. [187], there is still a need for more sophisticated experimental methods. In particular, the transverse shear properties are still difficult to get by direct measurement [101, 190].

One additional aspect, that has not been accounted for in elasto-plasticity models for papers is the plastic volume change that accompanies the densification effect. As shown in [239], this plastic volume change should be incorporated into the definition of the energies through \(J_p\) and hence considered in the Clausius-Duhem inequality. While this issue has been accounted for in models for other elastic-plastic materials—such as geomaterials [24, 25]—it has not been considered for paper, yet.

What is more, the combination of viscoelastic-viscoplastic models with mechanosorptive creep—as developed already for other materials such as wood [211, 212]—still needs to be applied to paper and paperboard. Even more, the fully coupled solution of the time-, rate-, moisture-, and temperature-dependence has not been achieved, yet. This is already challenging and computationally expensive in the isotropic case. However, the moisture and temperature transport might even be necessary to be described in an anisotropic fashion.

Furthermore, the coupling of the aforementioned effects with the continuum damage modeling of paper has not yet been tackled. Besides, the existing damage models for paper are limited to small deformations. Hence, there is a strong need to fill this gap by extending existing modeling strategies for the finite deformation regime.

5 Modeling of Paperboard Laminates

Paperboard is often used in packages. Then, usually laminates of several layers with different properties are applied. In principle, either the whole laminate can be modeled in a smeared manner with just one effective material model—see e.g. [195] for a recent application to modeling the drop test and compression test of a gable top package—or each of these layers can be described by the models presented in the previous section. In addition, several different materials are often applied that influence the overall performance of the laminate. For example, the softening of the aluminum foil which prevents light and oxygen penetration in beverage packaging can be considered [28, 29, 255].

Anyways, one additional aspect that should be given consideration is the interface behavior between the layers. This can be done by introducing cohesive zone models (CZM), that describe the softening behavior during delamination via the definition of a traction-separation law, in each interface between layers. In the simplest case, the layers can be modeled with a linear elastic continuum mechanical model, while the bonding between the layers is considered via a CZM [93]. There, an orthotropic elastic-plastic cohesive law has been incorporated through an additive split of the separation into an elastic and a plastic portion. A similar interface model has also been used in [188] in combination with an elastio-plasticity continuum model for the layers. Likewise, such a model has been implemented already into LS-DYNA as *MAT_COHESIVE_PAPER as material type 279 [129].

Classical cohesive zone relations for composites use bilinear traction-separation laws, see [43, 224, 245] to name only a few, but also more complex laws exist—such as tri-linear [1, 89], quadratic [74], square root [228], exponential [161, 185], and general power laws [108, 202]—for special applications, e.g. when fiber bridging occurs. In contrast, in paperboard applications, the majority of traction-separation relations are exponential [21, 94, 144, 171, 188, 244], and less works use bilinear ones [58, 59]. Alternatively, the tranction-separation relations can also be defined directly from experimental results from fracture mechanics tests by evaluating the critical J-integral [257, 258].

In most technical applications, the loading of the interface is not only pure peeling (mode I), pure shearing (mode II), or pure tearing (mode III). Instead, combinations of these loadings occur in the interface that delaminates. Consequently, one relevant aspect for cohesive zone modeling is taking such mixed-mode loading cases into account. This can be achieved by defining traction-separation laws based on effective quantities [21, 94, 188]. However, when mixed-mode scenarios are considered, it is essential to prove consistency of the chosen traction-separation laws for any possible loading path, because otherwise non-physical results can be obtained as shown in [66, 200]. To overcome these issues and to guarantee thermodynamic consistent results, cohesive model formulations should be based on potentials that couple the behaviors in different modes in a consistent manner. For example, the potential developed in the PPR-model [45, 201, 232] can be adapted to paperboard applications, as shown in [144]. Similarly, a new mixed-mode framework has been proposed in [58, 59] in which the unilateral effect is accounted for and interpenetration upon interface closure is prevented.

An additional challenge is the calibration of these models. The classical experimental tests that can be conducted to measure the delamination resistance of paper and board include the Scott bond test [80, 119], the z-directional testing based on lamination techniques [9, 79], the Y-peel test [94], and the short-span uniaxial tension tests [243]. The aforementioned procedures are useful for the out-of-plane tensile response. The corresponding shear response is even harder to investigate. For this, the notched shear test (NST) has been invented [191, 194]. In addition, the short-span compression test can be interpreted such that conclusions on the shear properties can be drawn [91, 101].

An alternative approach is adapting test setups, that have been initially developed for fiber-reinforced composites, to paper and paperboard. For example, the mode I behavior has been investigated based on the double cantilever beam (DCB) test in [59, 144, 148], while the mode II response has been analyzed based on the three-point (3ENF) [144, 148] and four-point (4ENF) [59] end-notched flexure tests. For mixed-mode loadings, the use of mixed-mode bending (MMB) tests has been suggested in [58], but actual tests have not been conducted. Alternatively, a non-proportional separation test has been performed in [144]. Even so, accurate and reliable measurements for the mixed-mode behavior of paperboard are still an unsolved problem.

6 Modeling of Corrugated Paperboard

Packaging boxes are often formed from corrugated paperboard, i.e. by gluing fluted sheets of paperboard as core material to flat face sheets of linerboard. When these materials are modeled, this specific structure needs to be addressed.

One important issue is the influence of creasing on the overall behavior of the corrugated paperboard. For example, in [97], it has been shown that the creases result in local reduction of the bending stiffness, which is desirable because it simmplifies the folding operation. Therein, creasing of a corrugated board panel has been modeled by letting a punch deform the liner and fluting elasto-plastically.

A most relevant property of boxes made from corrugated paperboard is their compression strength. One way to quantify the box compression strength is by performing box compression tests (BCT). Already in the 1960’s, such tests have already been performed and evaluated by analytical methods, such as the McKee equation proposed in [176]. Some more recent observations on this model, its implications, and possible improvements can be found in [56].

Nowadays, however, most predictive models for simulating the BCT are based on the finite element method. For instance, the critical buckling loads have been distinguished using an orthotropic elasto-plasticity model in [254], where different geometries and the influence of creases have been considered. In addition, the influence of material orientation has been analyzed in [163] by comparing boxes that were cut to have compression axes along MD and CD, respectively. What is more, the BCT has also been described numerically in [78]. There, also the compressive creep behavior of ventilated corrugated paperboard packages has been accounted for by changing loads as well as environmental conditions.

Another test to determine the compression strength is the edge crush test (ECT). This test has been investigated by FEA e.g. in [100], where plasticity was defined by Hill’s model and failure by the Tsai-Wu criterion. Also, corrugated paperboard samples have been researched in the ECT and boxes from the same material in the BCT in [76], where only linear elastic material behavior has been assumed. Finally, the ECT has also been examined in [134, 135] using a new test rig and comparing to simulations. The setup allowed to subject the specimen to cyclic climate change in a climate chamber, thereby triggering creep behavior.

Further, the influence of humidity and moisture on the material response of corrugated paperboard has been in focus of research. For example, the humidity effect on compressive deformation and failure of recycled and virgin layered corrugated paperboard structures has been investigated in [181]. Also, the steady-state moisture transport through corrugated fiberboard packaging has already been explored in [41]. Aside from that, the influence of moisture loading cycles in terms of accelerated creep has been surveyed in [114].

Until now, the question of sustainability of corrugated paperboard and the corresponding manufacturing processes has not catched a lot of attention, yet. This is suprising, since corrugated paperboard is an energy-efficient material which is reusable, fully recyclable, and simple to produce. In consequence, the carbon and water footprints of corrugated board are less than of any other packaging material [44]. Nonetheless, the supply chain and its environmental impact is not well investigated, yet. In fact, to the knowledge of this author, the only work in this field is [69]. In consequence, much more effort needs to be invested in this topic as is done for other packaging processes, see e.g. the review [177] and the references therein.

Lastly, also alternative paperboard sandwich structures could be evaluated. For example, crumpled papers have been investigated in [173], and it has been shown that these exhibit mechanical properties that are comparable to flexible and rigid commercially-available polymer or aluminum foams of similar densities. Thus, this could be an interesting material for substituting core materials of paperboard sandwich structures.

7 Modeling of Forming Processes with Paperboard

Several forming processes of paper and paperboard have been explored by numerical models. From these, a brief selection is described in the following. A more complementary review on these forming processes is provided in [198], while an overview with special focus on creasing and folding is given in [57].

7.1 Creasing and Folding

Two important converting processes in the forming of laminated paperboard to packages for products are creasing and folding. Creasing is the manufacture of fold lines such that the bending stiffness of the board is reduced preventing it from breaking during the following folding step. For this, delamination is introduced between the layers.

There exist numerous works in the literature on creasing and folding, in which the layers are modeled with relatively simple elastic-plastic material models. For example, Hill’s yield function has been applied for the paperboard layers of coated paperboard in [16, 17], whereas the coating layers have been modeled using the von Mises yield criterion for ideal plasticity.

As already mentioned above, delamination plays a crucial role in the creasing and folding processes and thus should be considered also in the numerical models. This can be achieved conveniently by cohesive zone formulations as shown e.g. in [21], where the orthotropic Hill model with isotropic strain hardening has been applied for the layers, while a cohesive zone model with normal and shear components has been introduced for the interfaces. The latter modeling approach has also allowed investigating the influence of the delamination description and different numbers of delamination surfaces [22].

In addition, the cohesive behavior can also account for differences in MD- and CD-direction by introducing two different shear components as shown in [110, 111]. With the help of this modeling strategy, numerical simulations could be performed with different shear strength profiles of the considered laminates by altering the ply and interface properties. Thereby, different production strategies could be mimiced. It was shown that the interface strengths mainly influenced the folding behavior, whereas the ply properties mainly affected the creasing force [109]. Further, the delamination during the folding of creased paperboard has been modeled in [83] through a special interface element with both, elasto-plastic as well as damage behavior.

In other works on creasing and folding of paperboard, more complex material models have been used. For instance, the anisotropic elasto-plasticity model of Xia et al. [253] has been adopted for the in-plane behavior of the paperboard layers in [192, 193], while the out-of-plane behavior has been treated through a cohesive interface model. Alternatively, the elastic-plastic model of [31] has been inserted in the simulation of the line folding process in [32]. Moreover, in [54] the orthotropic elastic-plastic material response has been described on the basis of the Ramberg-Osgood plasticity model [208], where the strain in MD oder CD, respectively, is defined as function of stress via a power law. Moreover, in order to account for the tension-compression asymetry, the isotropic plastically unsymmetric Drucker-Prager model [70] with extensions proposed in [199] has been adopted in [68] for creasing and folding simulations. Finally, a finite element analysis of the folding process of creased white-coated paperboard has been performed in [122] using a combined fluffing resistance and shear yield glue model.

Although a lot of work has been done concerning the modeling of creasing and folding of paperboard, there are still several open issues. In general, the material models used are relatively simple in most cases. Thus, more complex phenomena such as creep and relaxation have not been included yet into the simulations, although the effects have been shown experimentally to be influential [182, 183]. Further, deeper investigations concerning the optimization of the shear strength profile should be undertaken [189] and potential improvements to increase the peak bending angle in folding should be developed, e.g. by using a needled middle layer [131]. Last but not least, more sophisticated comparisons with results of alternative testing setups such as the Scott Bond Test could yield deeper insigth into the processes [63].

7.2 Forming and Deep-Drawing

The forming and deep-drawing of paperboard is a relevant field of research, because these processes differ significantly from applications in sheet metal forming [87]. Nevertheless, only few works are available in the literature, yet.

There are three main methods for producing advanced shapes of dry paper and paperboard [92]. While in press-forming a solid male and female die are used to convert the paper, in hydro-forming the male die is substituted by a pressurized membrane. In contrast, a male punch is pushed through a whole during deep-drawing.

One important converting process is press-forming, which is particularly suitable to produce deep trays. In order to evaluate the resulting creasing patterns that are produced during the converting, an orthotropic elastic-plastic model based on Hill’s yield criterion has been applied in [14]. An extension of the latter model to predict failure has been given in [15]. In addition, the delamination behavior has an important effect on the modeling of the forming process, as shown in [112].

The hydro-forming process has also been investigated via finite element analysis. For example, in [88] the isotropic elastic-plastic von Mises criterion has been applied. More complex material behavior has been implemented for instance in [154], where both the moisture- and the temperature-dependence have been accounted for. The comparison between press-forming and hydro-forming can be found in [92].

What is more, the deep-drawing process has been modeled with an explicite FEA in [250]. A statistical approach has been applied additionally in [249] in order to include the occurrence of rupture. Furthermore, wrinkle prediction as well as the post-wrinkle behavior have been included into the model in [155].

Lastly, the simulation of the packaging forming process is provided in [215]. There, anisotropic elastic-plastic model presented in [31, 32] has been applied successfully to the large-scale formin of a packaging product. The material model has been implemented into a solid-shell finite element formulation such that the thin, shell-like structure could be considered using the fully three-dimensional material model.

8 Conclusions

Paper and paperboard are characterized by their multi-scale and multi-physics nature, since several different physical effects can be observed on different scales. Further, the material behavior is nonlinear and anisotropic as well as dependent on time, loading history, strain rate, temperature, and moisture. Additionally, in most cases, large rotations and strains have to be considered. In consequence, all these aspects have been included into models for paper and paperboard. Nonetheless, it can be concluded that no model exists yet, that is capable of accounting for all these effects simultaneously.

Moreover, the forming processes that are applied to paperboard products include very complex structural steps with complicated loading states. Therefore, simplified models have been used frequently in order to minimize the computational effort. Although these simplified approaches are necessary for simulations of structures with practical relevance, the more complex multi-scale and multi-physics models have to be further developed to enable engineers to further optimize the material properties and the forming processes. Only then, the field of applications of paperboard can be extended to other branches such as building construction [13, 162] or electrical engineering [197].