Tray forming operation of paperboard: A case study using implicit finite element analysis

The possibility to perform advanced forming operations of initially plane paperboard paves the way to making products like food trays, plates, cups and other containers as a part of shifting towards a circular bioeconomy. As a part of the ongoing efforts of expanding the product range using paperboard, we performed analyses of the forming operation using simulations. An implicit non‐linear finite element model is built to more accurately than previous studies simulate the tray forming process of paperboard. Two different commercial paperboards are investigated. The use of an implicit solver enabled the inclusion of the creasing pattern into the geometry of the paperboard blank resolving the formation of wrinkles during forming. The material data is extracted from tensile test curves of the investigated paperboards and was fitted accurately using Hill's plasticity with difference in tension and compression accounted for with subsequent failure evaluation. The results showed that the inclusion of the creases in the geometry is vital for getting a correct shape of the formed tray and important for decreasing the risk of failure. The results also showed that friction has a big impact on the formed shape and hence on the stress levels, and therefore supports the means of lowering friction between the blank holder and the blank during the tray forming operation. A stochastic approach is proposed to determine the probability of failure for the boards. The performed failure evaluation is consistent with the field observations. The developed approach enables more precise simulations of paperboard tray forming.


| INTRODUCTION
Today, the possibility to perform advanced forming operations of paperboard has a big interest in the industry to create products like food trays, plates, cups and other containers. [1][2][3] However, compared to plastics and sheet metal, paperboard has lower formability. Deep drawing of sheet metal is widely used in many industries, [4][5][6] and plastics are used to form a variety of trays and containers. 7,8 At the same time, paperboard is an important part of shifting towards a circular bioeconomy and the effort to increase formability and the processes enabling the forming operation is ongoing.
There are several types of forming operations for paperboard, such as hydroforming, press (tray) forming, deep drawing and molding.
The first three of them are all variants of pressing down a sheet of paper or paperboard, called the blank, into a die. Hydroforming presses the blank downwards utilizing air pressure. The press forming uses a punch, also called male die, to press the blank into the bottom of a die, while deep drawing uses the punch to press the blank into a bottomless die. The methods are described in e.g. Hagman et al. 3 Östlund and Söderberg 9 and Lowe et al. 10 Complex products such as egg packages are molded. 11 As a part of the product development process using paperboard, simulations of the forming operations are used. The behaviour of the paperboard can then be studied, and it can be evaluated if and how material properties should be altered to minimize the risk of failure of the paperboard during the forming operation. The present study focuses on the tray forming operation of paperboard. A numerical approach is presented where the simulation of the operation is further developed compared to what has been done previously.
Deep drawing and tray forming of paperboard has been simulated in earlier studies, [12][13][14] with different numerical approaches and objectives. In Wallmeier et al. 12 the effects of varying the blank holder force, the die temperature, and the thickness of the paperboard were investigated. They showed, amongst other things, that the friction of the blank holder and die have significant effects on the stress in the blank, implying that low-friction dies and blank holders can considerably reduce the failure probability. In Awais et al. 13 the effect of the number of creases on the strain levels was numerically investigated for paperboard. The creases were modelled with hinge connector elements, and not explicitly included in the geometry. It was seen that the number of creases aided in lowering the strain levels. It could also be seen that this effect was greater for the first paperboard having a higher Young's modulus compared to the second investigated paperboard, which had a larger strain to failure. The approach with hinge connector elements for the creases is investigated more in detail in Livill et al., 14 where the approach allows for spontaneous wrinkling, that is, the number of creases or the position of them do not need to be known by forehand. The approach is included in deep drawing simulations and shows that the number of creases is about the same as in experiments. The disadvantage of these approaches is the inability to assess the effect of creases on the frictional interaction between the paperboard and the forming unit as the wrinkles were not resolved physically. The inclusion of the creases in a numerical model is important not only to simulate the correct shape, friction interaction and strain levels, but also to investigate the possibilities to add a lid onto the tray to seal it. For the sake of successful sealing operations, the upper edge of the formed tray should be as smooth as possible.
The advances in the current study are as follows: 1. The use of an orthotropic material model with isotropic hardening according to Hill plasticity, fitted to actual tensile tests and compared with two sets of data from papers with different failure properties. The fitting accuracy was improved, which allowed a reliable comparison of different materials.
2. Accounting for the difference in tension and compression for yielding and failure.
3. The detailed resolution of the creases which are introduced through geometrical features rather than through the ad-hoc methods, and the considerably higher number of elements compared to earlier studies, allowing for more exact results. This enables to study the effect of creases on the frictional interaction.
4. The use of implicit time integration during the simulations, which has a direct impact on the accuracy of the test. 5. The analysis of the frictional conditions required to achieve the desired shape.
With all the aforementioned parts included at once in the model, the possibility to simulate advanced forming operations is enhanced.
The objective is to increase the precision of the results coming from a simulation model of tray forming, giving more exact results for stresses, strains and the risk of failure compared to previous models.
The verification is based on achieving the correct shape of the tray, something not fully attained in earlier studies of tray forming where the flange of the corners has been larger than what is seen in production, even for failed trays. The risk of failure is also determined and compared to reality.
The paperboard is commercially produced under the name Inverform by Iggesund paperboard, which is a part of the Holmen Group. The paperboard is produced on two board machines, but differences in converting performance have been observed and to understand this difference the present investigation was performed.
The paperboard coming from the first board machine is here called Board A, and the paperboard coming from the second machine is here called Board B.
In Figure 1, the failure mode of the studied tray using Board A is shown, where the material fails in the corner under the creases. With the new simulation approach, the tray forming operation is simulated using the different paperboards, and the difference in results is evaluated and compared to reality.

| The tray
Paper and paperboard are anisotropic, heterogeneous, and hygroscopic materials. They are built up by a network of cellulose fibers from softwood or hardwood, or from a combination of the two. The fibers are randomly distributed over the sheet with, partly, random orientations.
The fibers are connected via fibre bonds which form spontaneously when water disappears from the web in the papermaking process. The anisotropic material behaviour must be considered in numerical simulations of paperboard. A common way is to model the paperboard as orthotropic, 3,[12][13][14][15][16] where material properties are specified in the paperboard machine-direction (MD), cross-direction (CD), and in the Z-direction (ZD), that is, through the thickness of the paperboard.
In Figure 2 the paperboard blank is shown as it is prepared by the tray manufacturer for the forming operation. The blank is laminated with a polymer that is extruded over the blank since the tray must withstand moisture during usage. The creasing pattern with 30 creases in each corner has been pressed into the paperboard so that it folds F I G U R E 1 The failure mode occurring for the studied tray using Board A. In (A) seen from the inside the tray, and in (B) seen from the outside F I G U R E 2 The paperboard blank before forming operation. (A) The full blank with dimensions and (B) close-up of the creases F I G U R E 3 The studied tray after a successful forming operation. In (A) seen slightly from above and in (B) seen from the side and shapes correctly. The grammage is 330 g/m 2 and the thickness of the paperboards, including the thin (30 μm) and compliant PET coating, is 0.5 mm. Figure 3 shows the studied tray after a successful forming operation, that is, without detectable failure. The linear dimensions of the formed tray are 185 Â 125 Â 25 mm.

| Material data
Paperboard is an anisotropic material, which may be approximated as an orthotropic material. The material shows different responses in tension and compression which may not always be captured by the materials models available in the standard libraries in the commercial finite element tools. The source for the input data was the physical tensile tests of the two paperboard types considered in this study.

| Elastic material properties
In the following theory, the principal material directions MD, CD and ZD are described with indices 1, 2 and 3, respectively.
The elastic part of the paperboard is described using Hooke's law ε ¼ Cσ: For an orthotropic material the full expression reads Equation 3 along with the symmetry condition of the compliance matrix, giving ν 21 =E 2 ¼ ν 12 =E 1 , give the in-plane elastic material parameters which are listed in Table 1.

| Hardening model
Plasticity in paperboard has been modelled in many different studies.
In Harrysson and Ristinmaa 20 a large strain orthotropic elasto-plastic model was developed with a yielding surface based on the Tsai-Wu failure surface, 21 which made it possible to directly introduce the difference in tensional and compressive yield behaviour for paperboard.
Several models based on the complex anisotropic yield surface introduced by Xia et al. 18 have been developed, such as the in-plane paperboard models established in Li et al. 22 and Tjahjanto et al. 23 The latter model is a viscoelastic-viscoplastic small strain approach developed to capture creep and relaxation for transient uniaxial loading. One of the latest publications on the subject is the one by Robertsson et al. 24 where the continuum model is based on previous models 15,25 using numerous sub-surfaces. In Robertsson et al., 24 results from simulations using solid continuum elements and shell elements are compared for some forming operations. They showed, amongst other things, that the shell elements had a better performance compared to the continuum elements. For the example simulating an actual forming operation from the industry, frictionless contacts, an explicit solver scheme and ideal plasticity were used.
The evolution of the plastic strains in the current study is described using Hill's plasticity, 26 which is suitable for composites and a common way to model plasticity for orthotropic composite such as paperboard. Hill's plasticity is defined as where F, G, H, L, M and N are defined as The Hill's parameters R ij in Equation 5 are defined as and determine the shape of the yield surface, which initial size is determined by the initial yield stresses σ ij y and the isotropic yield stress σ y . The material parameters in Equation 6 are found with the previously mentioned fitting procedure. The stress-strain curve measured in the CD is used as a master curve for the multilinear hardening and the rest of the parameters are fitted in Matlab using the fminconfunction to minimize the error between the measured and the calculated tensile test curves. The paperboard is modelled to yield at R p ¼ 0.0001, that is, the initial yield stress is in this study the stress giving a permanent deformation of 0.01% strain. Such a low value is required to get a good fit between the experimental and the numerical tensile test curves. The quality of the fit is shown in Figure 5. The curves on the compressive side are only from the numerical tests since no data is given from experiments. In compression, the two paperboards have a 30% reduction of the yield stress and a 50% reduction of the ultimate stress compared to the tensional side, which renders in the curves on the compressive side in Figure 5.
The shape of the Hill yield surfaces for the two boards is shown in Figure 6, plotted for zero shear stress, τ 12 ¼ 0. As seen in Figure  The fitting procedure resulted in different yield surfaces for the two paperboards (see Figure 6) Board B has a yield surface close to circular and an earlier yield point, compared to Board A. In Table 2, the complete set of parameters for the Hill's plasticity used in this study is presented. Note that the fit of R 33 is important even though plane stress is approximated, since R 33 influences the shape of the yield surface in the MD-CD plane.

| Failure evaluation
Failure is not included in the numerical model but will be evaluated as a part of the post-processing of the final results using the Tsai-Wu stress failure criterion. 21 For plane stress, it reads 27,28 where F i and F ij are defined as In Figure 7 the failure surface is shown for the two boards for zero shear stress. The shape of the two surfaces is very similar, but Board A fails earlier compared to Board B. It is obvious that the combined stress state allows for considerably higher stress levels than can be concluded if only the uni-axial tensile tests are studied.
Strain failure is also evaluated, here using maximum strain theory according to For zero shear strain, the failure envelop is a rectangle limited by the uni-axial tensile and compressive strains.
In Equation 10 the tensile shear strain must be estimated and is in this project estimated in the same way as the tensile shear stress as In Table 3 31 it was seen that the number of formed wrinkles for blank paperboard depend highly on the blank holder force, but also on the height of the formed product. In Awais et al. 13 an optimum blank holder force was found to avoid failure during the forming process, and in Tanninen et al. 32 the use of advanced force control was studied. The current study focuses on the difference between the two paperboards, and so the blank holder force will not be further evalu-

| RESULTS AND DISCUSSION
In the following, the trays are evaluated for stresses, strains, failure and the effect of creases and friction. All results are displayed at the end of the punch stroke step, that is, when the punch is pressing the blank towards the very bottom of the die. The results are displayed for the case of zero friction in the area shown in Figure 10 since this gave the best shape of the tray. The effect of altering the friction coefficient in the corner is shown later as a separate section.
In Figure 11 the shape of the simulated tray is shown and compared to the real tray. The shape of the simulated tray is to a high

| Total normal strains
The total normal strains in the MD and CD, sampled from the midplane, are shown in Figure 13. At the observed failure location, that is, the lower corner, the strain in MD is about 2% and 1.3% for Board A

| Failure results
The stress and strain levels for the analyses are used to study failure of the trays. In Figure 14 Figure 14. In both models, there are many locations where the strain state is outside of the failure surface, that is, ε F > 1, also in the lower corner for Board B where failure is not seen in production.
In Figure 15

| Failure evaluation
The Maximum Strain Theory was shown to be too conservative for this study. The Tsai-Wu stress σ TW for Board A implies that the risk of failure in the lower corner is very high, as seen in Figure 15A, where σ TW reaches 1.7 [-] which is far beyond the failure limit. Board B, however, has an area where σ TW reaches 0.8 [-] in the lower corner. In the following, these values are further analysed to show the probability of failure in these locations. The analysis is performed as follows: 1. Run FE analysis. where the scale η and shape β must be determined. The results of the fit of Equation 12 to the tensile tests in Figure 4 are shown in Table 4.
Since the Weibull distribution returns 0 for a negative variable, the distributions for the compressive stresses were created by dividing the tensile test data by two (since the failure in compression is 50% of that in tension) for the MD and the CD respectively, which were then used for fitting of the Weibull function.
The analysis continues with Monte Carlo simulations of the σ TW in the studied locations. As seen in Equation 7, the current stress state (σ 11 , σ 22 , σ 12 ) in the studied location is a part of the σ TW expression. In  The parameters in Table 5 only apply for the specific stress state   It should be emphasized that at σ TW ¼ 1 failure is initiated, and is   parameters, and the drying constraints during the manufacturing of paperboard will further contribute to a variation, as numerically proved in Alzweighi et al. 39 To include a variation of constitutive parameters requires further testing but, if acquired, is then straightforward to include in the analysis.
The approach may also be more exact if testing of the materials size dependency is performed. As mentioned earlier, the failure stresses and strains for paperboard have been shown 36 to have a size dependency, that is, locally allowing for higher stresses and strains than seen in standard tensile tests. Data for material size dependency was not accessible for this study and is hence not included, and in that sense the approach is conservative.
The effect of the polymer extruded on the paperboard blank before the tray forming operation could also be investigated. This coating may affect the upper ply allowing for higher failure strain. 40 Finally, the use of shell elements leaves the contribution of the out-of-plane (ZD) stress outside of the analysis. The ZD stress will influence the stress-based failure surfaces, locally allowing for both lower and higher stress levels.

| Effect of friction and creases
In the following the numerical model is used for investigating the influence of friction and the importance of including the creases in the geometry. We used only Board A for this evaluation as the results are similar to Board B.
The effect on the shape of changing the friction coefficient in the corner (the area shown in Figure 10) is shown in Figure 18, where the friction coefficient in (a) is 0.3, in (b) 0.1 and in (c) 0. It is obvious that the friction makes the forming operation harder, and that it is an advantage of keeping the friction low in the corner. The shape closest to reality is given by the simulation using zero friction in the corner, shown in Figure 18C, which is also the model used for the above presented results.
In Figure 19, the effect of the friction in the corner on the Tsai-Wu stress σ TW is shown. Here it is seen that the area where σ TW > 1 is very large for the cases with a friction coefficient of 0.3 and 0.1. In Figure 19C the friction coefficient is zero, and this is the same results as shown in Figure 15A. The results strongly support means of lowering the friction, such as the use of lubricants or increased die tool temperature, in the forming operation to aid problems with failure and bad shapes.
In Figure 20  Theory shows small differences between the two paperboards but is deemed too conservative. The failure evaluation using Tsai-Wu theory shows that Board A has a very high risk of failure, and that the risk of failure using Board B is lower, something that is in full agreement with what has been seen in production. The results from the stochastic failure evaluation based on the numerical model more precisely suggest that the risk of failure in the lower corner using Board A is close to 100%, and for Board B the risk of failure is about 1%. No known problems with Board B are reported from production and the here calculated failure risk is probably conservative.
The model may be further used to estimate the most critical material properties for the tray forming such as elongation, strength, and friction. Among these, in this study, we only probed the impact of the friction in the corner. A stochastic approach is used as a part of the post processing to study the impact of strength. However, the stochastic nature of paper makes it difficult to include the elongation and strength directly in the model as these parameters can vary locally and therefore require characterization and modelling which includes statistical distributions of the paper properties on the relevant scales.
The latter is outside the scope of this work.