Finite-element simulation with a homogenization model and experimental study of free drop tests of corrugated cardboard packaging

– This paper presents experimental and numerical studies of drop tests of corrugated cardboard packaging containing diﬀerent foam cushions. An eﬃcient homogenization model for the corrugated card-board has been developed. In our homogenized model, the corrugated cardboard is represented by a 2D plate. Instead of using a local constitutive law (relating the strains to the stresses) at each material point, the homogenization gives the global rigidities (relating the generalized strains to the resultant forces) for the equivalent homogeneous plate. This model was implemented into the FE software ABAQUS. The foam behaviour was experimentally determined and modelled using a crushable foam model of ABAQUS. The packages are tested in free fall from a given height on a rigid ﬂoor. The deceleration of the packed product was recorded using a triaxial accelerometer. The numerical results obtained using the FE simulation with the homogenized model agree well with the experimental results. We have also shown that the contribution of the corrugated cardboard box to the shock response could not be neglected in the design of cushioning package.


Young modulus along direction i νij
Poisson ratio Gij Shear modulus {N } Internal membrane force vector {T } Internal transverse shear force vector {M } Internal moment vector [A] Membrane stiffness matrix [D] Bending-torsion stiffness matrix [F ] Transverse shear stiffness matrix [B] Membrane, bending-torsion coupling matrix {ε} Membrane strain vector {κ} Bending-torsion curvature vector {γc} Transverse During storage and transportation, products can accidentally fall onto the floor causing some damages on the products.The corrugated cardboard packages and the foam cushions are designed to protect the product from the shock it may undergo.As the product design always tends towards light, a high impact performance becomes extremely important for product design issues [1,2].Conventionally, a product reliability test to prevent impactinduced damage is carried out by a procedure of "design -prototype -test -redesign" which is highly cost and time consuming.A numerical modelling of the product and its packaging provides an efficient methodology to predict the structural strength during impact.Moreover, the finite-element simulation allows to avoid numerous experimental tests and to predict possible failures during the design stage.
The conventional method of packaging design uses cushion curves, and assumes that the foam elements of the package are rectangular blocks, loaded in uniaxial compression [3].For the design of complex shaped foam cushions the finite-element simulation was used in different works [4][5][6][7].For the corrugated cardboard, the mechanical properties under various solicitations were investigated by some authors [8][9][10].By using some FE models and commercial FE software the mechanical behaviour of corrugated cardboard was also studied [11][12][13].However, there is no published research on the design of packaging systems with both foam cushions and corrugated cardboard boxes.
This article deals with the drop impact problem of packaging systems using numerical and experimental methods.Different packaging systems, using corrugated cardboard boxes and different configurations of foam cushion, are studied in terms of shock resistance.
Firstly, we present the materials and experimental setup.Then, we describe the homogenization model for the corrugated cardboard, its implementation in the finite element code ABAQUS and the foam behaviour model.Finally, the decelerations measured during the shock by using the data acquisition system of the drop test platform are compared with the results obtained by ABAQUS simulations using the homogenization model.

Corrugated cardboard box
Several corrugated cardboard boxes with their inside dimensions of 150 × 150 × 90 mm 3 and 4 mm thickness are used for the experimental drop tests (Fig. 1).The product contained in the box is a thick steel plate of 87 × 87 × 30 mm 3 .The corrugated cardboard is composed of outer Kraftliner, inner test-liner and a flute recycled-liner.The Kraftliner is mainly made with virgin softwood and hardwood fibres but some recycled fibres are also added (e.g.clippings from corrugated board manufacture or from old corrugated containers).The test-liner mainly consists of selected recycled fibres.The recycled-liner is made with recycled fibres and its quality is lower than that of the test-liner.These constituents are separated according to the French norm NFQ03-043 in order to determine the material parameters of each one.The tensile test samples are cut using the ELCEDE MFT20 cutting table.The tensile tests are performed on the two liners and the fluting along the Machine Direction (MD-0 • ), Cross Direction (CD-90 • ) and 45 • Direction.The tensile strain-stress curves for the liners and fluting in the three directions are shown in Figure 2. It is seen that the strain-stress curves clearly depict the anisotropic elastoplastic hardening behaviour.For the Kraft-liner (Fig. 2a) and the test-liner (Fig. 2b), the Young's modulus and the tensile strength decrease when the angle increases from 0 • to 90 • ; because the fibres tend to be aligned parallel to the rolling direction (0 • ) during the manufacturing process.However, the fluting is made by running the steam-heated paper between two toothed rollers (like two wide gears) to form a sinusoidal shape.This process breaks some fibres along the machine direction so decreases its MD Young's modulus (Fig. 2c).
Many paper grades have an almost linear elastic behaviour and then brittle behaviour in the MD direction, but they can have a notable ductility in the CD direction (Fig. 2a).The rate of recycled fibres used in the composition of the cardboard changes this behaviour (Figs.2b  and c).In this study, the above paperboards are considered as elastic and orthotropic, having the following constitutive relations:  where the subscripts x and y refer to the MD and CD directions respectively.The Young's modulii are obtained by calculating the slopes of the tensile curves in the elastic zones.The empirical relations established by Baum et al. [14] are used to determine the Poisson ratios and the in-plane shear modulus as follows: These relations are used to obtain the elastic properties of the two liners and the fluting given in Table 1.

Foam cushions
The foam cushions are cut in Ethafoam 400 extruded polyethylene (PE) foam having a density of 58 kg.m −3 and supplied by Dow Chemical Company.In this study, two configurations are considered (Fig. 3): the foam cushion surrounding the entire product (Fig. 3a) and the foam cushion limited to the four corners of the product (Fig. 3b).
Uniaxial and hydrostatic compression tests are performed to characterize the foam.The nominal stresses versus nominal strains are plotted in Figure 4.These curves are used to estimate the elastic modulus and the hardening law of the foam.

Drop test system
The drop test is used to determine the package ability to retain and protect its contents under a shock due to a free fall.The drop test is a simple method and corresponds well to the reality.The package is placed on a fork at a given height, then the fork is withdrawn very fast and the package falls to the ground.The impact velocity depends on the drop height.This velocity is a crucial parameter, leading to great forces during a short time.
The drop test equipment used here is an AccruDrop Tester from L.A.B. Equipment Inc. (Fig. 5).The drop height is fixed to H drop = 1.05 m.To record the acceleration of the product, a triaxial accelerometer (Andevco Isotron 2258A-10) is glued to its horizontal top face.The analysis and the visualization of the measurements are done with Lansmont Test Partner 3 software.

Material behaviour models
The material models used in the FE simulations are presented in this section.The 3D corrugated cardboard panels are homogenized to 2D orthotropic elastic plates.This H-model is then implemented into ABAQUS using user's subroutine UGENS [15].The crushable foam model of ABAQUS is used for the foam which is considered as an isotropic material; with a hardening law in function of the volume change.

Corrugated cardboard homogenized model
A 3D geometrical modelling of the liners and the fluting of the corrugated cardboard is a very tedious and time consuming task.In our homogenization model, a corrugated cardboard panel is replaced by a 2D plate.Instead of using a local constitutive law (relating the strains to the stresses) at each material point, the homogenization leads to global rigidities (relating the generalized strains to the resultant forces) for the equivalent homogeneous plate.
Using the theory of laminated plates and integrating the stresses through the thickness (z) and then along the MD direction (x), we obtain the following generalised constitutive law [12]: where {N }, {T } and {M } are the internal forces and moments; [A], [D], [B] and [F ] are the stiffness matrices related to the membrane forces, the bending-torsion moments, the bending-torsion-membrane coupling effects and the transverse shear forces respectively; {ε} is the membrane strain vector, {κ} is the curvature vector and {γ c } is the transverse shear strain vector.The corrugated cardboard is more complex than a laminated plate because of the fluting core and the cavities between the two liners.Consequently some global effective stiffnesses in the matrix (3) obtained by the theory of laminated plates should be modified [12,16].

Traction and bending stiffnesses
On the MD section, if we apply a traction force (N x along x) or a bending moment (M x around y), the flute behaves like a thin membrane so it can almost not withstand the normal stress along x.Our theoretical and numerical studies [12] indicate that the flute has no contribution to the extension and bending rigidities related to N x and M x .The following changes are required in the matrices [A], [B] and [D]: unfold the flute to a flat plate, the shear stress τ xy can be easily calculated: where G f xy is the shear modulus of the fluting, e f its thickness and l its curvilinear length defined by: l = 4 with where h is the corrugation height (see Fig. 1).
Thus we obtain the in-plane shear behaviour of the fluting: Finally, the shear stiffness of the corrugated cardboard is obtained by adding the contribution of the two liners: where e o , e f and e i are the thicknesses of the outer liner, the fluting and the inner liner respectively; G o xy , G f xy and G i xy are their shear modulii.

Transverse shear stiffness on CD section
It is difficult to directly determine the rigidity of transverse shear on the CD section because of the coupling of bending and transverse shear.According to the reciprocity theorem, Nordstrand et al. [17] proposed to replace the transverse shear under the shear force T y (along z) by a shear on the thickness under a force along y.
The shear stiffness on the total thickness for corrugated board is obtained as follows:

Transverse shear stiffness on MD section
The transverse shear under the shearing force T x on the MD section (along z) is replaced by a shear on the thickness under a force along x.Thus this shearing problem becomes a bending problem of the fluting and liners which was solved by using the curved beam theory [17].Finally, the transverse shear stiffness relative to T x is obtained as follows: where G * is the apparent shear modulus calculated using the analytical formulas developed by Nordstrand et al. [17] and implemented into our H-model.

Torsion stiffness on CD and MD sections
In a previous work [16], we have theoretically demonstrated that the torsion of an orthotropic plate could be decomposed into two beam torsion problems, thus the torsion stiffness of a corrugated cardboard has been obtained: where GJ MD and GJ CD are the beam torsion stiffnesses on MD and CD sections respectively (GJ MD GJ CD for corrugated cardboards), L is the length along x and B is the width along y (see Fig. 1).
The beam torsion stiffness GJ CD can be analytically obtained using the Bredt's formula for a closed thin wall beam section:

Implementation of the homogenized model in a nonlinear FE analysis
For this dynamic non linear problem, the corrugated cardboard is modelled by shell elements, the internal forces and moments are calculated incrementally.Their increments {dN } and {dM } have the following relationship with the increments of the membrane strains and bending-torsion curvature: where the stiffness matrix [H] is determined by our homogenization model which is implemented into ABAQUS by programming a user subroutine UGENS.
In ABAQUS, the 4-nodal shell element S4R has 6 degrees of freedom per node (3 displacements and 3 rotations).In the UGENS subroutine, the generalised forces and strains are passed in using the arrays FORCE (6) and STRAN ( 6), and will be updated at the end of the subroutine; the known generalised strain and curvature increments {dε} and {dκ} are passed in using the array DSTRAN (6).The main purpose of the UGENS is to calculate the stiffness matrix [H] in Equation ( 15) which is stored in the array DDNDDE (6,6) at the end of the subroutine.
This array will be used by ABAQUS to define the element stiffness matrix in the FE solution.The transverse shear stiffness is directly introduced in the input file of ABAQUS.

Crushable foam model
Foam materials have a cellular structure and hence behave in a complex manner, especially under crush conditions.A typical stress-strain curve of polymeric foam has four distinct parts, which are caused by different mechanisms.A linear elastic region, followed by an almost horizontal plateau (Fig. 4) resulting from the cells collapsing; they buckle elastically and recover completely when the foam is unloaded.When the cells collapse completely, the stress rises steeply (densification) [18][19][20].For this part, a phenomenological non associated plasticity model is proposed for isotropic crushable foam materials.It consists of a modified plasticity model that accounts for a different behaviour in compression or tension, which incorporates uniaxial and hydrostatic compressive hardening laws.
This model with volumetric hardening uses a yield surface with an elliptical dependence of the deviator stress tensor [s] on the pressure p in the meridional plane.The elastic behaviour of the foam can be considered as isotropic and linear.The yield surface and the flow potential for the crushable foam model are defined in terms of the hydrostatic pressure stress p = − 1 3 trace [σ] and von Mises equivalent stress σ eq = 3 2 [s] : [s].The yield surface is defined as follows: where the parameter α represents the shape of the yield ellipse in the p-σ eq stress plane, p 0 is the centre of the yield ellipse, σ is the length of the vertical axis (σ eq ) of the yield ellipse; α can be calculated from the initial yield strength in an uniaxial compression (σ 0 c taken as a positive value), the initial yield strength in hydrostatic compression (p 0 c ) and the yield strength in hydrostatic tension (p t ): The plastic flow potential g for non associated plasticity is given by: The parameters of this model were identified from the experimental results of Figure 4, giving k = 1.1 and k t = 0.

Finite-element simulations
The product is considered as a rigid body.The corrugated cardboard box is modelled with 4-node S4R/ABAQUS shell elements and is considered as a homogenized material.The foam is modelled with 4-node C3D4/ABAQUS tetrahedral elements and is considered as elastoplastic material.The ground is considered as a fixed flat rigid body.A frictionless contact is assumed.The package is positioned 1 mm from the ground and has an initial velocity V = 2gH drop , with the drop height H drop = 1.049 m.In Figures 6 and 7, we have plotted the von Mises stresses for the box and for the two considered configurations (configuration 1: foam cushion surrounding the product, configuration 2: foam cushion at the corners of the product).We can see that the corrugated cardboard box undergoes higher stresses than the foam.This indicates that the box contributes also to the shock absorption.The highest stresses on the box are located on the central region of the topside.In fact, the cardboard under a very great deceleration is not supported by the foam cushion in this region.

Comparison of drop test experimental results and FE simulations
Two types of test are performed for the configurations 1 and 2 to study the influence of the packaging on the product protection: the product is protected only by the foam, or by both the foam and the corrugated cardboard box. Figure 8 shows typical curves of deceleration vs. time corresponding to the shock on the ground.The deceleration is presented in terms of the vector resultant magnitude.We observe that in the case without the corrugated cardboard box, the maximum acceleration  (deceleration) is 20% higher than that in the case with the box.This can be explained by the fact that the box strengthens the foam package in the lateral directions and gives some damping effects on the underside of the foam.Five dropping tests have been done for each configuration; the mean and the standard deviation of the maximum deceleration and the shock duration are given in Table 2.
The comparisons of experimental and numerical curves of deceleration vs. time for the configurations 1 and 2 in the case of underside shocks are shown in Figures 9 and 10 respectively.Figure 11 presents the comparison of the experimental and numerical curves of deceleration vs. time for the configuration 2 in the case of lateral shocks.We can see that the finite element simulation results agree well with the experimental results.In the drop tests, the underside surface of the package must be well flat and horizontal in order to obtain the underside impact.Otherwise, some impacts may take place on the package edges.In Figures 9 and 11, only one peak is present on the experimental results.However, we observe two peaks on the experimental curve in Figure 10.The reason is that the underside is not entirely flat and horizontal, so the first shock on an edge gives the first peak and the second shock on the opposite edge immediately gives the second peak.It is evident that an impact test in two shocks gives a less maximal deceleration and more impact duration than that in one shock.
In Table 3, we present the maximum deceleration and the shock duration for different configurations.Between the experimental and numerical results, we have obtained an error of 2.87% for the maximum deceleration in the case with the foam around the product, -13.96% in the case with the foam on the corners of the product and -11.76% in the case with the foam on the corners but the shock on the lateral side.
It can be noted that the corrugated cardboard box contributes about 20% to the shock absorption.The comparison between configurations 1 and 2 shows that the surrounding foam offers a more rigid support to the product leading to a shorter shock duration (more than 50%) but a greater maximum deceleration (about 50%).
The comparison between the underside (Fig. 10) and lateral side (Fig. 11) shocks for configuration 2 shows that the shock duration is almost the same and the maximum deceleration is about 7% lower for the lateral shock.This means that the foam cushions at the corners of the product also provide a good protection of the product on the lateral side of the box.

Conclusion
In the design stage of packaging systems, we need to simulate the contributions of both the corrugated     ferent foam cushion configurations.The drop tests were instrumented, giving the curves of deceleration vs. time.The FE simulation results agree well with the experimental results.We have also shown that the corrugated cardboard box with the corner foam cushions gives a more damping effect to the shock response of the product.
In this study, we have not presented the results for the shock on an edge or corner of the box.In fact, in these two cases the corrugated cardboard box undergoes a large plastic deformation, damage and rupture.These phenomena are not yet included in our homogenization model.We intend to extend our work on elastoplastic homogenization model considering large deformations of the corrugated cardboard.
shear strain vector P Corrugation period of the cardboard L Length of a cardboard along MD B Width of a cardboard along CD l Curvilinear length on a corrugation period of a fluting h Height of the corrugation of a fluting eo, e f , ei Thicknesses of outer liner, fluting and inner liner respectively G * Plastic flow potential for non associated plasticity (k, kt) Parameters of the crushable foam model 1 Introduction

Fig. 6 .
Fig. 6.Von Mises stresses (in Pa) for configuration 1: (a) on the corrugated cardboard box and (b) on the foam cushion (section).

Fig. 7 .
Fig. 7. Von Mises stresses (in Pa) for configuration 2: (a) on the corrugated cardboard box and (b) on the foam cushion.

Fig. 9 .
Fig. 9. Comparison of experimental and numerical curves of resultant magnitude deceleration vs. time for configuration 1.

Fig. 10 .
Fig. 10.Comparison of experimental and numerical curves of resultant magnitude deceleration vs. time for configuration 2.
cardboard box and the foam cushions.According to the drop tests, it is evident that the contribution of the corrugated cardboard box to the protection cannot be neglected.An efficient homogenization model for the corrugated cardboard has been developed and implemented into the FE software ABAQUS through the user subroutine UGENS.This model has been used to simulate the drop tests of the corrugated cardboard box containing dif-

Fig. 11 .
Fig. 11.Comparison of experimental and numerical curves of resultant magnitude deceleration vs. time for configuration 2 with lateral shock.

Table 1 .
Material characteristics of the corrugated cardboard constituents.

Table 2 .
Maximum deceleration and shock duration for drop tests.

Table 3 .
Comparison of the results obtained by the drop test and FE simulation.