A continuum damage model for creasing and folding of paperboard

A continuum damage framework is proposed for modelling the delamination process occurring in paperboard during mechanical loading. The main application of interest is line creasing and subsequent line folding used in package forming. To adequately capture creasing and folding, a continuum damage framework that uses a single isotropic damage variable that evolves with the plastic strains associated with out‐of‐plane shearing is used. The damage evolution is calibrated against folding experiments for a specific reference mesh, and a simple scaling strategy is proposed to reduce the inherent mesh dependency. To highlight the potential of the proposed model, an illustrative 3D example is considered where a paperboard sheet is creased by two plates and folded to resemble the corner of a package.


| INTRODUCTION
Paperboard is a cellulose-based material, and it is one of the main components used in packages for food products.On sheet level, paperboard is anisotropic due the manufacturing process where fibres are continuously sprayed on a moving web.The majority of the fibres will align in the web direction which is denoted the machine direction (MD) while the in-plane orthogonal direction to MD is defined as the cross direction (CD).In addition, the stacking direction of the fibres is termed the out-of-plane z direction (ZD).The elasticity modulus along MD is typically two times the magnitude of CD while ZD can be around 100 times lower compared to MD.The degree of anisotropy between the in-plane directions can be reduced by controlling the fibre direction; however, the large anisotropy between the in-plane and out-of-plane direction remains and this is one of the prime characteristics of paperboard.
Paperboard models exist on different scales, and they can be divided into fibre network scale models, where individual fibres are modelled, and continuum scale models in which the fibre network is approximated as a homogeneous continuum; see, for example, the review article by Simon. 1 The individual fibres in a fibre network can be modelled using beam elements.This idea was presented in Heyden 2 for dry fibres and further developed in Kulachenko and Uesaka 3 where the bonding between fibres were modelled using contact mechanics.In Borodulina et al, 4 cohesive zones were used to more accurately account for fibre-fibre bonding.Another approach was presented in Persson and Isaksson 5 where the end points of fibres and the fibre bonds were modelled as particles using the discrete element method (DEM).In Lavrykov et al, 6 the finite element method was used to model individual fibres with solid elements, which requires more computational effort compared to the beam approach but is important when studying the response in the stacking direction of fibres.While network models can be useful in order to tailor the material properties during the production phase, they are not suitable for large-scale applications, such as converting processes, for example, creasing and forming, which require computationally efficient continuum models.
The application of interest in the current work is the converting process that involves creasing and folding.Creasing is a process in which the material is locally deformed along predefined lines in order to guide the upcoming forming and create stable folds; see Coffin and Nygårds. 7Creasing followed by folding has been simulated in the works of Huang et al 8,9 and Beex and Peerlings. 10,11In these works, the in-plane behaviour of paperboard was modelled with relatively simple constitutive models using small strain theory while the outof-plane deformation, associated with large deformation and delamination, was modelled by cohesive interfaces.A challenge in this approach is to identify the location of the cohesive interfaces, which is a complex and important task in, for instance, 3D-rotational creasing; see Borgqvist et al. 12 Another method for modelling creased lines was developed in Giampieri et al 13,14 where an interface element, placed between two shell elements, was tailored to capture the response during folding.
An alternative continuum scale approach is to model paperboard using large strain theory where the complex out-of-plane deformation is included in the continuum constitutive model.Such models are typically based on thermodynamics with an anisotropic yield surface; see, for example, Xia et al. 15 In this class of models, the evolution laws are able to account for irreversible deformation in both the in-plane and out-of-plane directions.One model within this category was developed by Li et al 16 and Li et al 17 which was used to simulate a cylindrical compression test.Another anisotropic large strain elasto-plastic model was developed in Borgqvist et al 12 where line creasing was considered.This model was later utilized in Borgqvist et al 18 to study the short-span compression test and folding of uncreased paperboard.In addition, the model was combined with a solid-shell element in Robertsson et al 19 where the forming sequence to produce a package was simulated.Furthermore, in Robertsson et al, 20 the model developed by Borgqvist was extended to include viscoelasticity and viscoplasticity to model the experimentally observed rate dependency of paperboard.For materials that share characteristics with paperboard, notable works include the anisotropic viscoelasto-viscoplastic model for pressboard by Tjahjanto et al 21 and the large strain elasto-plastic model by Harrysson and Ristinmaa 22 for corrugated board.While models accounting for large strains have considered creasing and folding separately, for example, 12,18 no continuum model of this type has been able to obtain a realistic response during creasing followed by folding.As this process is the basis for packaging forming, this will be addressed in current work.
In the current work, continuum damage modelling (CDM) is combined with the framework in Borgqvist et al 12,18 to model the softening response during creasing and subsequent folding.CDM was originally introduced by Kachanov 23 and later formalized by, for instance, Lemaitre. 24The continuum damage concept has been used to model the in-plane deformation of paperboard in Isaksson et al. 25,26 In our model, a single damage variable is introduced since shearing is assumed to be responsible for the delamination process and the local material failure.The dissipated energy due to outof-plane shearing drives the damage evolution.The proposed framework shares characteristics with the ductile damage model presented in Larsson et al. 27 The damage model is calibrated against creasing and folding experiments.A simple scaling method is proposed in order to mitigate the mesh dependency caused by the damage evolution.In the constitutive driver a staggered approach is utilized in the sense that damage evolution is applied after the elasto-plastic response has been calculated.
The article is structured as follows.In Section 2, the overall framework together with the specific modelling approach for paperboard is presented.In Section 3.1, the line creasing and subsequent line folding is simulated and compared to simulations without damage.
Furthermore, a mesh study is performed using the proposed scaling method.In Section 3.2, an illustrative 3D example is simulated where paperboard is creased by two plates and folded to resemble the edge of a package.

| Modelling base
The elasto-plastic description of the model herein closely follows the framework in Borgqvist et al, 12 and the reader is refereed to this source for additional details.The model is described in the spatial configuration with the dissipation inequality on the form where τ is the Kirchhoff stress tensor and d ¼ symð _ FF À1 Þ is the symmetric part of the spatial velocity gradient.A multiplicative split of the deformation gradient is assumed F ¼ F e F p where F e and F p are the elastic and plastic contributions, respectively.The Helmholtz free energy is postulated as where α models isotropic damage and Ψ e 0 and Ψ p 0 are the undamaged elastic and plastic contribution while κ ðνÞ with ν ¼ f1,…,12g being the internal hardening variables.The state variables that define Ψ e 0 are the elastic finger tensor b e ¼ F e F eT and the structural tensors m ðiÞ , introduced as m ð1Þ ¼ v ð1Þ v ð1Þ , m ð2Þ ¼ v ð2Þ v ð2Þ and m ð3Þ ¼ n ð3Þ n ð3Þ , ð3Þ where v ð1Þ and v ð2Þ are the in-plane characteristic directions of paperboard in the spatial configuration, that is, the MD and the CD, while n ð3Þ is the characteristic out-of-plane (ZD) direction.The n ð3Þ direction is defined to be perpendicular to the in-plane directions such that n ð3Þ ¼ v ð1Þ Â v ð2Þ holds.With this split, the dissipation (1) becomes From argument by Coleman and Gurtin, 28 we assume D e ¼ 0 such that τ ¼ ð1 À αÞτ 0 with which coincides with the Kirchhoff stress, τ 0 , in Borgqvist et al. 18 The remaining part of the dissipation in (4) reads where the conjugated forces are given by Associated plasticity is assumed which renders the evolution laws where f is the yield surface and _ λ is the plastic multiplier.The yield surface, which together with the Karush-Kuhn-Tucker (KKT) conditions, governs the evolution of plasticity.The undamaged part of the Kirchhoff stress is used to describe the yield surface.Similar to Borgqvist et al, 18 the yield surface is defined as where n ðνÞ s projects the stress tensor on the subspace ν while K ðνÞ 0 is the initial distance to the yield surface and k is a material parameter that controls the smoothness of the yield surface.Using the evolution equations ( 8), the yield surface (10) and assuming elasto-plastic loading, that is, f ¼ 0, the dissipation in (6), excluding the damaged part, can be written as

| Evolution of damage
The evolution law for the damage variable is inspired by Larsson et al 27 who postulated the damage evolution in a similar manner to plasticity theory, that is, where Φ is a potential function and _ μ a Lagrangian multiplier.The dissipated energy A T is defined as where P is a subset of the full summation in (11) and t being the current pseudo time.The potential Φ in ( 12) is chosen as where the material parameters c n control the damage evolution.The critical pseudo time t c for the initiation of damage is reached when the following criteria is fulfilled: where β n are material parameters.Similar to plasticity theory, the KKT conditions are assumed to control the onset of the damage evolution In conclusion, when t > t c , the criteria Φ ¼ 0 holds and the damage becomes The model above fulfils the dissipation inequality (6) since Ψ 0 ≥ 0 and _ α / P n P Dn which follows from (13).The model in ( 14) and ( 15) with the material parameters β n and c n allow for a flexible format when combined with the anisotropic yield surface in (10).As such, the contribution from the plastic dissipation for different modes of deformation can be combined in the expression for the damage evaluation.
For instance, the out-of-plane deformations associated with delamination can be influenced by, for example, debonding and fracture of fibres during in-plane tension; see Johansson et al. 29

| Specific elasto-plastic model
The undamaged part of the elastic energy, Ψ e 0 , is chosen in accordance with 18 as where J e ¼ detðF e Þ and H is the Heaviside function with the properties, H ¼ 1 if I 13 ≥ 1 otherwise zero.The invariants are defined as provides the undamaged part of the Kirchhoff stress, where the scalar functions are summarized in (A1).
The 12 projections n ðνÞ s that define the yield surface are given as where the normalized vectors v ð1Þ ¼ v ð1Þ =jv ð1Þ j, v ð2Þ ¼ v ð2Þ =jv ð2Þ j and v ð3Þ ¼ n ð3Þ =jn ð3Þ j were introduced.The parameter m in n ðνÞ s governs the plastic coupling between out-of-plane compression and shear while ν 12 and ν 21 are the in-plane Poisson numbers.Illustrative snapshots of projections for the yield function are shown in Figure 1.
Following Borgqvist et al, 18 the plastic spin tensor is chosen as which is the skew-symmetric counterpart to d p ¼ sym obtained from (8a), ( 9) and (10).Using the definition _ 8), ( 10) and ( 22) results in the elasto-plastic evolution equations, The specific format for the hardening K ðνÞ will be considered later on.The material model is implemented in the commercial finite element software Abaqus Standard as a Umat subroutine; see, for example, Dassault Systems. 30In the numerical implementation, the evolution equations in (23) are discretized using the Euler backward method and solved with the Newton-Raphson procedure to obtain the internal hardening variables κ ðνÞ and the plastic part of the deformation gradient F p .See Appendix B for additional details regarding the implementation of the model in Abaqus.

| Specific damage model
The motivation for introducing damage in the paperboard model is to capture the complex behaviour of paperboard during creasing and subsequent folding.As the purpose of the creasing process is to weaken the material in specific zones, damage is needed to capture the material softening.The driving force for the damage evolution, see (17), is taken as the undamaged plastic dissipation associated with the out-of-plane shear deformation since this mode is observed to have large influence on hinge formation during subsequent folding.
The theoretical framework allows for further damage coupling, for instance, damage evolution during out-of-plane tension such that a corresponding uniaxial test can be predicted.However, for simplicity, such couplings are omitted in the current study.The assumptions implies that the subset P ¼ 9,10,11,12 f ggoverns the damage evolution such that ( 15) and ( 17) become The parameters c n and β n associated with out-of-plane shearing along MD, n f9,10g, and CD, n f11,12g, are calibrated against line creasing and folding when the MD and CD, respectively, are aligned with the creasing setup using a specific element area A e r as reference discretization.The reference area is chosen such that the thickness of the elements are on the same order as the thickness of a fibre.It is well known that local damage models give rise to pathological mesh dependency since the damage variable tends to localize in few elements.A sophisticated way to mitigate the mesh dependency is to introduce non-local variables, using, for instance, gradient enhancement methods, where the damage is spatially averaged for a characteristic length; see, for example, Peerlings et al. 31 Since the material model is implemented in a commercial software, an alternative engineering approach to mitigate the mesh dependency is proposed where the length scale is based on the finite element discretization.
As such, from empirical observations, the following scaling of the damage parameters are suggested: where A e md and A e cd are projected areas along the MD-ZD and CD-ZD in the undeformed configuration.For an eight-node brick element, the areas are computed from where A e i are the six sides of the brick element with its corresponding normal directions n i .Notably, the scaling (25) is very simple to implement as it is explicit.

| Calibration procedures
The paperboard used in Borgqvist et al 18 is considered in the current work, that is, thickness and density are 0.39 mm and 788kg=m 3 , respectively.The paperboard is entirely composed of softwood fibres with representative length of 3:6mm and a thickness range of 10 -25μm.During the manufacturing process, the paperboard is designed as a single ply structure; hence, the material parameters are assumed to be uniform along the thickness direction.The elastic parameters, E i , in (18) are tabulated in Table C1.Following Borgqvist et al., 32  Similar to Borgqvist et al, 18 plastic hardening variables in (10) are modelled by where the hardening parameters a ν and b ν together with the initial yield distances K ðνÞ 0 are tabulated in Table C2.The hardening parameters for in-plane shear and tension, that is, ν ¼ f1,2,3,6g in (27), are calibrated against uniaxial tension along MD, CD and 45 while the hardening a 7 and K The exponent in the yield surface (10) is set to k ¼ 3 while the coupling parameter between out-of-plane shear and compression is estimated to be m ¼ 0:7 from Stenberg et al. 35 The damage parameters in (24) will be identified by comparing the model to creasing and folding experiments as discussed in next section.

| Line creasing and folding
To evaluate the response of the proposed model, we simulate the line-creasing setup in Borgqvist et al 12 followed by line folding.The geometry and boundary conditions for the creasing and folding steps are schematically outlined in Figure 2. The paperboard sheet is 110 mm long with a width of 38 mm while having a thickness of 0.39 mm.Since the width of the sheet is large compared to the localized area of deformation, plane strain is assumed.
Prior to creasing, the displacement of the paperboard edges are prescribed u x rendering an initial in-plane force of 57N that mimics webtension.The paperboard is then creased by displacing a male die, u z ¼ 0:2mm, into a female grove; see Figure 2A.After the webtension and dies are removed, the irreversible deformation of the material creates two localized shear zones that acts like hinges.
The line folding following creasing is initiated by two clamps compressing the paperboard sample such that the contact pressure becomes 0:2MPa.Rotation of the paperboard strip is prevented by a rigid block which is placed at a distance of 10mm from the clamps.
The paperboard is then folded by prescribing the rotation of the top and bottom clamps as illustrated in Figure 2B.
As previously mentioned, the simulation is performed using the commercial finite element software Abaqus Standard with the material model implemented as a Umat subroutine.The paperboard is discretized with linear continuum brick elements where the reference element size, at the critical region close to the creased zone, is 0:07mm in-plane and 0:0156mm through the thickness, that is, A e r ¼ 0:0011mm 2 .Finite deformation and quasi-static conditions are The calibrated damage parameters are tabulated in Table C3 where a reference element discretization of A e r ¼ 0:0011mm 2 was used.
The results during line creasing and line folding are presented in Figures 3 and 4.
As seen in Figure 3, the difference between the simulations with and without damage during creasing is minuscule.During folding, see Figure 4, the damage model predicts a plateau in the mechanical response at approximately 20 which is consistent with the measurements.This contrasts the model without damage which is unable to accurately predict the folding response as the force grows with increasing folding angle and does not reach the plateau that is experimentally observed.
The damage distribution after creasing and after creasing and folding is shown in Figure 5.It is observed that although the damage generated during creasing is relatively large, see Figure 3A, the influence on the mechanical response is negligible.Furthermore, in Figure 5B, the damage variable during folding has reached its maximum level for several elements which is expected for the creation of hinges.The increase in damage at the two hinges allows for a flexible response where the middle part of the crease deforms plastically from out-of-plane tension, which for simplicity is modelled as ideal plastic.
In Figure 6, the simulated deformation after folding with and without damage is compared to a representative experimental deformation pattern.Close to the hinges, the simulation with damage evolution included predicts a more flexible response than the simulation Reaction force during creasing when the experimental setup is aligned with (A) MD and (B) CD.
F I G U R E 4 Folding curve after creasing when the experimental setup is aligned with (A) MD and (B) CD.
without damage.This is consistent with the experimental finding in Figure 6C.In addition, the total plastic out-of-plane shear, defined as where it is concluded that the shear strain is concentrated at the hinges in both simulations.
Ductile damage models are, due to the softening, known to be plagued by inherent mesh dependency.Based on numerical experiments, the damage parameters are proposed to be scaled according to (25) where a rectangular reference area of A e r ¼ 0:0011mm 2 was used close to the critical region.A mesh size sensitivity study is performed in Figure 7 where the MD is aligned with the experimental setup and the ratio between the element thickness and length is kept fixed.The considered range of mesh sizes A e =A e r are separated into two figures for visibility.As seen in Figure 7, the refinement of the mesh has a large impact on the response if the scaling is omitted.In contrast, the results using the proposed scaling strategy (25) are in close agreement with the simulation of the reference mesh for the considered range of A e =A e r .In Figure 8, the deformation is shown at the folding angle 75 for three mesh sizes using the scaling approach.The deformation and damage distribution appears similar with a concentration of damage at the hinges and delamination at the middle part of the crease.While the deformation is similar, the amount of delamination is larger for the simulation with the finer mesh.
To further evaluate the model, the simulated response for two additional crease depths, when the setup is aligned along MD, is shown in Figure 9 and compared with measurements.As seen in Figure 9B, the damage model is accurately predicting the increase in bending force for the two additional crease depths.In addition, the model prediction with and without damage is simulated for uncreased F I G U R E 7 Mesh size sensitivity study.The element size ratio A e =A e r for constant aspect ratio is plotted against the folding response.male die is displaced towards the female grove until a crease depth of u z ¼ 0:2mm is reached, that is, the same crease depth as for the line-creasing example.The scaling in ( 25) and ( 26) is used, and the element size along the creased lines is 0:035mm in thickness and 0:07mm in width.In total, the paperboard is discretized with approximately 100 000 continuum eight-node brick elements.
The most critical region for the paperboard is at the intersection between the crease lines where the mesh density is increased along the x-y direction as shown in Figure 11B.The MD of the paperboard is aligned with the x direction; see Figure 11A, and before creasing is initiated, an initial displacement along the sheet edges is applied such that the initial webtension in the x direction becomes 57N.
The solution procedure is similar to the line-creasing example, that is, the quasi-static balance of linear momentum is solved using the implicit Abaqus Standard solver.The friction coefficient between the paperboard and the tools is μ ¼ 0:3 while the maximum amount of damage is α max ¼ 0:7.
After the creasing step, the folding phase that consists of three folding sequences is initiated.The procedure is shown in Figure 12 where the first folding sequence consists of restricting the motion of the top left (not shown) and top right quadrant of the paperboard, by engaging clamping tools that compress the paperboard with a pressure of 0:2MPa.In the first sequence, the rotation shown in The deformation after creasing and folding is shown in Figure 13A.The critical part of the package is at the vicinity of the corner.To highlight the deformation pattern, Figure 13 shows the damage variable when cutting out the corner of the package.
In Figure 14, the cut-out corner in Figure 13D is tracked throughout the deformation history.It is observed that the deformed pattern after the second folding step is highly complex and coincides with the creation of the flap in Figure 12C.In addition, the amount of damage is localized around the creased regions which is expected for successful folding.

| CONCLUSIONS
A ductile damage model for paperboard has been developed to enhance the capabilities of a large strain elasto-plastic continuum model with high degree of anisotropy.The main application was creasing and subsequent folding.For this endeavour, the plastic dissipation from the out-of-plane shear deformation was chosen as the driving force for damage evolution.An isotropic damage variable was used since shearing is assumed to be responsible for the delamination concept.The large-scale problem was numerically stable which allowed an implicit quasi-static solver to be used.At the central junction of the creased pattern, the deformation was shown to follow a complex history which could be studied to achieve insight into the forming process.Overall, the deformation followed the creased pattern as expected which lead to the formation of a realistic package.
To further validate the model, additional paperboard grades will be considered in future work.Specifically, paperboards produced as multi-ply structures will be addressed which requires a methodology to extract the nonuniform material parameters throughout the thickness.In addition, more sophisticated methods to mitigate the mesh dependency such as introducing non-local variables will be studied in future work.The scalar functions for the undamaged Kirchhoff stress in (20) are presented in (A1).The constitutive evolution equations (23) are discretized using the Euler backwards method which results in the residual equations with notation ð • Þ pre referring to the previous state of equilibrium.A compact notation for the residual system is obtained by introducing To achieve quadratic convergence when solving the balance of linear momentum with the implicit finite element method, the algorithmic tangent matrix is needed.In Abaqus, see Dassault Systems, 30 the format for the algorithmic tangent reads By use of τ ¼ ð1 À αÞτ 0 and the chain rule, we obtain since τ 0 is a function of F e , see (5), and the damage variable is an implicit function of the chosen internal variables, κ ðnÞ with n P; see (17).To calculate the partial derivatives ∂F p ∂F and ∂κ ðnÞ ∂F in (B5), RðYð FÞ, FÞ ¼ 0 is used such that where ∂F p ∂F and ∂κ ðnÞ ∂F can be extracted from

∂Ψ e 0 ∂b e b e þ ∂Ψ e 0 ∂m ð1Þ m ð1Þ þ ∂Ψ e 0 ∂m
ð2Þ m ð2Þ À m ð3Þ ∂Ψ e 0 ∂m ð3Þ þ ∂Ψ e 0 ∂m ð3Þ : m ð3Þ I , j being the Euclidean norm.The invariants I 11 , I 12 and I 13 are associated with the elastic stretch along MD, CD and ZD, respectively, while I 23 is associated with the area spanned by the in-plane characteristic directions, v ð1Þ and v ð2Þ .Use of the elastic potential(18) in(5)

E 1 ,
E 2 , E 4 and E 5 are obtained from uniaxial tension experiments along the MD, CD and 45 directions by comparing the initial stiffness tensor resulting from (20) with the initial orthotropic stiffness tensor.The in-plane shear modulus is computed from Lekhnitskii, 33 and ffiffiffiffiffiffiffiffiffiffiffiffiffi ffi ν 12 ν 21 p ¼ 0:30 is used for the in-plane Poisson ratios that were obtained by Baum et al 34 for a range of paperboards.With this, we obtain ν 12 ¼ 0:46 and ν 21 ¼ 0:19.The elastic parameter E 3 is estimated from the initial stiffness during out-of-plane tension while E 6 and E 7 are obtained from unloading during out-of-plane compression.

ð7Þ0
are fitted to out-of-plane compression experiments.The in-plane compression parameters, K ð4Þ 0 and K ð5Þ 0 , are estimated from the short-span compression test while the out-of-plane tension strength K ð8Þ 0 is fitted to out-of-plane experiments.Due to lack of reliable experimental data, the out-of-plane shear strengths K to Borgqvist et al12 where they were estimated from the force-displacement response during creasing.
Schematics of the line-creasing setup and (B) the subsequent line folding procedure.The radius of the tools are R ¼ 0:1 mm, all dimensions in mm.assumed, and the problem is solved using implicit time integration.The tools are modelled as rigid bodies.Contact between tools and paperboard are modelled using a tangential penalty formulation with friction coefficient μ ¼ 0:4.The maximum amount of damage allowed in a Gauss point is α max ¼ 0:96, whereafter the damage evolution is terminated.The experiments and simulations are performed with the setup aligned in the MD and CD.During line folding, the bending force at the rigid block and the rotation angle is monitored.Similarly, during creasing the male die force and displacement is monitored.For both operations, the experimentally measured response is compared to the simulation using the proposed model with and without damage.During calibration of the damage given by (24), the parameters associated with Subsurfaces 9 and 10 are assumed to be identical as they correspond to negative and positive shearing.Similarly, the parameters associated with Subsurfaces 11 and 12 are assumed to be identical.In addition, since no plastic dissipation evolves for Subsurfaces 11 and 12 when the experimental setup is aligned in the MD, the damage parameters associated with Subsurfaces 9 and 10 are decoupled from the 11 and 12 parameters during calibration.Hence, two damage parameters are fitted to the response in the MD and CD, respectively.

F I G U R E 5
Damage distribution and local deformation during (A) creasing and subsequent (B) folding at 90 .F I G U R E 6 Shape of the hinges after the folding procedure with (A) the damage framework and (B) without.A representative image of the hinge (C) from the experiment.
paperboard and compared with measurements in Figure9B.The uncreased simulations are similar to the response obtained in Borgqvist et al,18 and a minuscule difference is observed with and without damage since in-plane tension and compression are the dominant modes of deformation.This observation indicates that a realistic approach for damage evolution is used.Furthermore, as the bending stiffness obtained in Figure9Bfor uncreased paperboard is in close agreement with the measurements, the uniform assumption for material parameters along the thickness is deemed sufficient for this class of paperboard.

3. 2 |
Figure 10A.The paperboard is placed between a female and male die that has a horizontal, vertical and diagonal crease pattern; see Figure 10B,C.

F I G U R E 8
Figure 11A,C, while the remaining mesh consists of two elements in the thickness direction.As seen in Figure11C, the width of the male die and the female grove, at a cross section along the crease lines, are identical to the previous line-creasing setup; see Figure2A.In addition, the

Figure 12A !
Figure 12A !B is applied where the active blocks are highlighted in blue.In the second folding step, shown in Figure 12B !C, the lower rigid block is rotated to obtain a paperboard flap along the diagonally creased pattern.In the final sequence, Figure 12C !D, a rigid block is rotated to fold the paperboard flap and complete the structure of the package corner.
process.The model without damage was shown to overestimate the force response during the folding process while the model including damage evolution could match the measured response and the deformed structure, after creasing and folding, resembled a representative image found in experiments.A simple scaling strategy of the damage parameters was proposed and shown to reduce the inherent mesh dependency of the damage formulation.To highlight the potential of the proposed model, the industrially important creasing and folding operations were simulated to create a package as a proof of F I G U R E 1 4 Snapshot of the cutout corner of the package after (A) creasing, (B) first folding, (C) second folding and (D) third folding.