Mechanical Behavior Modeling of Containers and Octabins Made of Corrugated Cardboard Subjected to Vertical Stacking Loads

The aim of this paper is to characterize the mechanical behavior of corrugated cardboard boxes using simple models that allow an approach to the load capacity and the deformation of the boxes. This is very interesting during a box design stage, in which the box does not exist yet. On the one hand, a mathematical model of strength and deformation of boxes with different geometry is obtained from experiments according to the Box Compression Test and Edge Crush Test standards. On the second hand, a finite element simulation is proposed in which only the material elastic modulus in the compression direction is needed. For that, corrugated cardboard sheets are glued to build billets for testing, and an equivalent elastic modulus is obtained. This idea arises from the fact that the collapse of the box is given by the local bucking of the corrugated cardboard panels, due to the slenderness itself, and the properties in the compression direction are predominant. As a result, the numerical models show satisfactory agreement with experiments, concluding that it is an adequate methodology to simulate in a simple and efficient way this type of boxes built with corrugated cardboard.


Introduction
Corrugated cardboard boxes are used in the transport of goods, housing loads of more than one ton of various products, such as frozen food in bags, bulk substances, rigid pieces and so on. In recent years, its use has increased with the emerge of e-commerce, making the analysis of these boxes a subject of interest [1,2]. According to [3,4], this type of product generated a market of USD 70 billion in the UK alone, and just in the first ten months of 2020. With an increase of nearly 10% of the global courier, express and parcel market from the previous year, goods transportation has become a strategic issue on the agendas of firms and administrations [5,6].
These boxes are usually mounted on pallets in order to be manipulated. Due to the need to optimize warehouses, the containers are stacked one on top of the other, so that the lower container can support a load of several tons, besides the specific actions that the pallet base does on the edge of the ring. Initially, a low range of standardized sizes was common in the industry, but nowadays, companies like Amazon have set what is called the "science of packaging", by increasing the diversity of sizes, and therefore minimizing the usage of raw materials and optimizing transport. Cardboard is key in this configuration, due to its versatility, price compared with other alternatives and environmental footprint [4,7]. These cardboard boxes are designed not only by considering how various holes (to handle it) might affect the compressive strength but also taking into account how the vibrations during transportation might affect the products delivered [1,[8][9][10][11]. It is always essential to define efficiently the product design based on the characteristics requested by the client and the technical infrastructure with the quality that guarantees its functionality and response for the quality it has been conceived [12].
The load capacity of the boxes is determined by a compression test called Box Compression Test (BCT), which is carried out with standard humidity and temperature conditions (50% and 23 • C). The maximum applied load is determined by dividing the result of this test, named the BCT load, by a safety factor. This safety factor is determined by taking into account several variables, such as the humidity, the type of pallet, the load eccentricity, the effect of punctual loads, the vibration that occurs during handling and transport, fatigue and operational life of the container. The safety factor can be obtained by tables in standards as ASTM D4169. Usually, safety factors greater than two are adopted.
Since corrugated cardboard is a highly deformable material, the limit of its use may be delimited by the deformation of the box, so that the upper pallet does not touch the content of the lower box. For design purposes, it is interesting to have a model that allows obtaining an approximation of the BCT load or the deformation that the box will undergo. Numerous empirical studies that propose a way to determine the BCT load have been published, [13][14][15][16] among others. The most commonly used model is McKee's model. This model allows obtaining the load capacity BCT from the Edge Crush Test (ECT) (that gives the load per unit of length ECT) of the cardboard box, perimeter p and wall thickness h as follows, where the parameter m result is 5.876 for the tested boxes. Regarding deformation, forcedeformation curves have been widely studied by several authors. Urbanik analyzed corrugated cardboard extensively, from plates made in the laboratory to box specimens [17]. His approach is based on testing individual samples. Corrugated cardboard is considered a homogeneous material. All the samples are single-wall cores made in a laboratory. The predicted results were 6-7% lower than the ones obtained by McKee's formula. However, the aforementioned formula is useful because of its simplicity and accuracy. For other authors [18][19][20][21][22][23], the corrugated cardboard is a structure that consists of several layers of paper in which a flat layer, usually called a liner, is joined to a wavy layer (from up to bottom). Certain ligatures are assumed at the junction points. In these models, the properties of the cardboard are determined by tensile and compression tests of the paper. In essence, they allow simulating the compressive behavior of the cardboard through the elastic properties of the combination of papers. This type of paperboard structurebased design has been used extensively to analyze the optimal paper combinations in the manufacture of corrugated cardboard. Nevertheless, its application results in difficulty in the design of packaging since the measures differ depending on the client's requirements. Therefore, an empirical mathematical model that allows a fast way to determine the stress and the deformation of a box is useful to select an adequate thickness and load capacity of the cardboard sheet.
Some studies [24][25][26] simplify the structure by transforming the fluting into a homogeneous material with equivalent elastic properties. In this way, the corrugated cardboard is transformed into a composite of three layers: the upper and the lower layer have the same mechanical properties as the paper, and the central layer has equivalent elastic properties and the same thickness as fluting has.
The ECT test, according to DIN EN ISO 3037 standard, is the most economical way to characterize the properties of a material. The ECT test is performed in a sample of 100 × 25 mm cut from the box material, in such a way that it can represent the entire box if a strength correction factor that is a function of some geometrical scale factor is applied. In this context, in [27], authors analyze corrugated cardboard as a homogeneous material with isotropic elastic properties. They perform diverse tests on cardboard beams in different directions, and they found that the elastic modulus is lower than 200 N/mm 2 . It is an analysis developed over thick cardboard, such as the ones analyzed in this article.
Concluding, in this paper, a characterization of the mechanical behavior of corrugated cardboard applied to containers and octabins subjected to vertical stacking loads is presented. More than thirty commercial-use packaging boxes are tested aimed at obtaining the force-displacement curve. From the results, McKee Equation (1) is applied in order to obtain the coefficient m for the tested samples, and a useful mathematical relationship to predict the deformation of the boxes is proposed. Finally, a methodology to simulate the behavior of the boxes by finite elements is proposed, in which an equivalent isotropic elastic modulus is considered for the cardboard, although the material is actually orthotropic. This is given as the assumption that the collapse of these kinds of containers occurs by local buckling of the vertical cardboard panels, so the material properties in the compression direction are predominant over the others. In conclusion, it should be pointed out that the mechanical behavior of this type of boxes built with corrugated cardboard can be easily and efficiently simulated from a practical and engineering point of view.

Materials and Methods
In this section, the experimental program is presented. Specifically, BCT, ECT and compression tests are performed for three different purposes. The first objective is to determine the coefficient m of McKee's model given by Equation (1) for three sets of boxes of different geometries. The second one is to obtain a model for the vertical deformation and for the deformation modulus of these boxes. The third objective is to obtain an equivalent Young modulus of corrugated cardboard blocks from compression tests performed on billets built of glued corrugated cardboard sheets, aimed at using this modulus in numerical simulations.

Description of the BCT and ECT Compression Tests
The BCT test was carried out in accordance with the UNE 131000 standard. This test is performed on box samples, and the maximum load before the sample collapses is measured, namely BCT. From this force, the BCT stress S y,BCT can be calculated by where p represents the perimeter of the box and h wall thickness. Besides, a forcedisplacement curve can be reached. For example, Figure 1 shows three curves corresponding to the test performed on three of the samples tested in this work (specifically, the Id. R-05 item of Table 1 that is presented later). In these curves, three different zones are distinguished. In zone A, there is high initial deformation due to the crushing of the lower flaps, or due to the crushing of the plate that forms the bottom. In zone B, the ring is fully loaded linearly and elastically. If the sample unloads in this area, the box returns to zero vertical deformation. The second zone ends with the point of maximum load BCT. In zone C, the material of the ring is under excessive stress, or the deformation is so high that cracks occur in the ring, weakening it. Once the BCT force is exceeded, the box deforms excessively until it collapses. From Figure 1, it can be highlighted that the three curves of three "identical" boxes are fairly different. This is an inherent problem in this type of testing on these materials because the manufacturing processes of both the cardboard sheets and the boxes are not repetitive. Deformation is considered from the point where the basis of the box sits (end of zones A and B), which is the point where the maximum stress is reached, as indicated as an example in the abscise axis of Figure 1. For example, in that case, is computed as 24−12 = 12 mm. Thus, from this test, the deformation modulus BCT can be computed as where z is the height of the sample. The deformation also allows us to obtain the deformation of the box in zone B caused by the compression of the vertical panels of the samples , given by The ECT standard is performed according to DIN EN ISO 3037. All specimens have dimensions of 100 × 25 mm, the material and thickness h are the same as the corresponding box. From the test, the load per unit of length ECT is obtained (the length is 100 mm, according to the standard), from which the stress ECT ,ECT is given by Taking into account Equations (1), (2) and (5), the m coefficient of McKee's model yields where m result is 5.876 for McKee's model, and f is the shape factor of the box, given by = .
All the tests were performed with the same hydraulic press. The laboratory environment was air-conditioned.

Id.
Type Deformation δ is considered from the point where the basis of the box sits (end of zones A and B), which is the point where the maximum stress is reached, as indicated as an example in the abscise axis of Figure 1. For example, in that case, δ is computed as 24−12 = 12 mm. Thus, from this test, the deformation modulus E BCT can be computed as where z is the height of the sample. The deformation δ also allows us to obtain the deformation of the box in zone B caused by the compression of the vertical panels of the samples δ u , given by The ECT standard is performed according to DIN EN ISO 3037. All specimens have dimensions of 100 × 25 mm, the material and thickness h are the same as the corresponding box. From the test, the load per unit of length ECT is obtained (the length is 100 mm, according to the standard), from which the stress ECT S y,ECT is given by Taking into account Equations (1), (2) and (5), the m coefficient of McKee's model yields where m result is 5.876 for McKee's model, and f is the shape factor of the box, given by All the tests were performed with the same hydraulic press. The laboratory environment was air-conditioned.

Description of the Box and ECT Samples
For the BCT test, industrial use packaging of rectangular and octagonal geometry is tested. The footprint perimeter of the packages ranges between 1855 and 5856 mm, and the thickness ranges between 6.6 and 20 mm. Boxes of thickness greater than 13 mm have been manufactured using the ECOBOX technology patented by "Cartonajes Lantegi S.L.", which allows lamination of single, double or triple layers of cardboard before shaping the box in a die-cutter adapted to the thickness (see Figure 2a). The height of the boxes is between 140 and 1950 mm, and the sides of the faces range approximately between 200 and 1500 mm. Table 1 shows, for all the samples, the identification Id., type of article (rectangular, octabin or ECOBOX) and the geometric data of the boxes: thickness h, height z, footprint perimeter p, larger side a, shorter side b, for octabins the length of the corner side x (see Figure 2b), shape factor f and the number of BCT and ECT samples n BCT and n ECT , respectively (the number of samples depends on the discarded results and the material availability at the time of testing). The samples for the ECT test are obtained by cutting 100 × 25 mm samples from lateral panels of the boxes.
All the samples were maintained in a climatic chamber for at least 72 h, in an environment at 20 • C and a relative humidity of 50%. However, the testing machine was not inside the climate chamber, given its enormous size. Nevertheless, less than 4 min elapsed between sample extraction and testing, so it is not assumed that there is a significant variation in humidity in the sample. This is important due to the significant variation in load capacity that containers experience when the humidity of the cardboard increases. and 1500 mm. Table 1 shows, for all the samples, the identification Id., type of article (rectangular, octabin or ECOBOX) and the geometric data of the boxes: thickness h, height z, footprint perimeter p, larger side a, shorter side b, for octabins the length of the corner side x (see Figure 2b), shape factor f and the number of BCT and ECT samples n BCT and n ECT, respectively (the number of samples depends on the discarded results and the material availability at the time of testing). The samples for the ECT test are obtained by cutting 100 × 25 mm samples from lateral panels of the boxes. All the samples were maintained in a climatic chamber for at least 72 h, in an environment at 20 °C and a relative humidity of 50%. However, the testing machine was not inside the climate chamber, given its enormous size. Nevertheless, less than 4 min elapsed between sample extraction and testing, so it is not assumed that there is a significant variation in humidity in the sample. This is important due to the significant variation in load capacity that containers experience when the humidity of the cardboard increases.

Description of the Corrugated Cardboard Block Samples
Simple compression tests are carried out on solid blocks formed by successive sheets of corrugated cardboard (see Figure 3) in order to obtain an equivalent modulus. The billets were tested with the same press that is used for the BCT and ECT tests, and the samples are climatically conditioned in the same way. The billets are manufactured by gluing with white glue, plates are cut with a cut plotter of the brand ESKO (model called Kongsberg); hence, the edges are not squashed. The plates have dimensions that vary from 300 × 300 mm 2 to 500 × 500 mm 2 . The plates are made of conventional double-wave cardboard (BC channel, made with Kraft paper of 180 to 300 g/m 2 and semi-chemical paper of 160 g/m 2 ). Between 6 and 30 of these plates are glued together, resulting in corrugated cardboard billets, with dimensions indicated in Table 2. The variables presented in Table 2 are the following: the identification Id., number of glued sheets num, thickness of a sheet H, total thickness HA num × H, width of the sheet B and the height Z.

Experimental Results
The results derived from the tests are: • For boxes, the m coefficient to obtain a model for the BCT-ECT relationship, a model for the vertical deformation δ as a function of the stress, and a model for the deformation modulus as a function of S y,ECT also. • For corrugated cardboard billet samples, the equivalent elastic modulus E as a function of the stress S y,B (B index refers to billets).

Results of the BCT and ECT Tests
The obtained results from the BCT and ECT tests performed on the samples are shown in Table 3. These results are, on the one hand, the experimental results of the maximum load BCT, the deformation δ, the stress S y,BCT given by Equation (2) and the deformation modulus E BCT given by Equation (3). On the other hand, the load per unit of length ECT and the stress S y,ECT computed by Equation (5).  As previously mentioned, it is very usual to find no repeatability in the tests due to the fabrication processes of materials and boxes. Hence, the results of the BCT load and the ECT load per unit of length of Table 3 are given as mean value plus/minus the mean absolute deviation (when only one sample was tested, N.A. is indicated instead of the deviation value). The other results are given only as mean values.
In order to illustrate the mechanical behavior of the boxes and the ECT samples, some force-displacement curves are shown. As examples, the case of an octabin and an ECOBOX are considered.
First, Figure 4a presents 10 BCT curves for the octabin O-05 item, and Figure 4b other 10 curves relative to the ECT test samples. The cyan curve of BCT test gave a very low value for BCT force, and it was discarded for the analysis. Then, the BCT results were comprised between 14.3 kN and 18.4 kN: the BCT mean value resulting from the remaining nine curves was 16.2 kN, and the mean absolute deviation 0.9 kN (5.6%). In relation to ECT results, the minimum value was 12.3 N/mm and the maximum one 14.2 N/mm. The average value was 13.2 N/mm, and the deviation was 0.5 N/mm, i.e., 3.4%.   Then, Figure 5a presents 10 BCT curves for the ECOBOX E-09 item, and Figure 4b the 10 curves of the ECT samples. The BCT load was comprised between 20.8 kN and 28.8 kN: the mean vale was 24.9 kN, and the deviation was 1.5 kN (6.0%). Concerning ECT results, the minimum value was 18.9 N/mm and the maximum one 26.4 N/mm. The mean result was 24.1 N/mm, and the deviation was 1.6 N/mm (6.6%).  In summary, for some cases, the deviation is lesser than 1% (e.g., E-12 and R-11), but in most cases the deviation is between 5-10%, as the previously illustrated examples. Although it should be noted that in some cases it reaches between 15-20% (e.g., O-01 and E-10). This dispersion is unavoidable when working with this type of material.

Mathematical Models Derived from the BCT and ECT Tests
The first objective is to obtain the coefficient m from Equation (6). In Figure 6, the relationship between S y,BCT and S y,ECT / f is represented, and a linear correlation for each type of box is performed. The result is m = 7.63 for octabins, m = 5.97 for boxes made of ECOBOX cardboard, and m = 5.61 for rectangular boxes made of double or triple corrugated cardboard. The mean value is m = 6.24. These values are not very far from the McKee's result, m = 5.876.
To characterize the stress-deformation relationship, in Figure 7 the deformation δ u calculated by Equation (4) is represented as a function of the stress S y,ECT . A correlation for the average value δ u is achieved, Therefore, according to Equation (4), a model for the deformation δ as a function of the stress S y,ECT yields From Figure 7, it can be noted that the higher the resistance, the lower the deformation, i.e., the higher the stiffness. The non-linearity present in this model is given because the deformations are relevant, and lateral deformations also occur given by buckling. This buckling appears progressively and locally in certain zones of the panels of the samples (this is visualized in the simulations of Section 4). The general mechanical behavior of the boxes is governed by this instability produced by the stress in the compression direction.
Therefore, according to Equation (4), a model for the deformation as a function of the stress ,ECT yields = 0.59 ⋅ ⋅ ,ECT . .     Finally, a relationship between the deformation modulus E BCT of the boxes and the stress S y,ECT is represented in Figure 8. The deformation modulus E BCT is computed by Equation (3) and gathered in Table 2 for the three types of boxes. This relationship between E BCT and S y,ECT of Figure 8 can be modelled by the next equations: E BCT = 9.74·S 3.23 y,ECT for octabins (in blue), E BCT = 7.89·S 3.32 y,ECT for ECOBOX (in red) and E BCT = 7.49·S 2.98 y,ECT for rectangular boxes (in discontinuous black). In order to obtain a unique equation to model the deformation modulus as a function of the stress as an average for the three types of boxes, the curve fitting procedure (in continuous black) results in tween BCT and ,ECT of Figure 8 can be modelled by the next equations: BCT = 9.74 ⋅ ,ECT .
for rectangular boxes (in discontinuous black). In order to obtain a unique equation to model the deformation modulus as a function of the stress as an average for the three types of boxes, the curve fitting procedure (in continuous black) results in BCT = 7.4 ⋅ ,ECT . .
(10) Figure 8. Relationship between deformation modulus and stress for the corrugated cardboard boxes.

Results of the Compression Test for Corrugated Cardboard Blok Samples
The results from the compression test performed on corrugated cardboard billets are shown in Table 4. The results directly obtained from the test are the maximum load R and the deformation . To illustrate the mechanical behavior of the billets, Figure 9 shows the load-displacement curves of the items B-01 and B-03 (see Table 2 for details of the samples). It can be remarked that the curve of the thickest sample, i.e., B-03, is nearly linear until failure, while the one that corresponds to the thinnest one (the B-01 sample), slightly loses linearity a bit before failure. However, from a practical engineering point of view, both curves can be considered linear.

Results of the Compression Test for Corrugated Cardboard Blok Samples
The results from the compression test performed on corrugated cardboard billets are shown in Table 4. The results directly obtained from the test are the maximum load R and the deformation δ. To illustrate the mechanical behavior of the billets, Figure 9 shows the load-displacement curves of the items B-01 and B-03 (see Table 2 for details of the samples). It can be remarked that the curve of the thickest sample, i.e., B-03, is nearly linear until failure, while the one that corresponds to the thinnest one (the B-01 sample), slightly loses linearity a bit before failure. However, from a practical engineering point of view, both curves can be considered linear.
Comparing the load capacity R of the specimens B-01 and B-03, they are 8.87 kN and 61.9 kN, respectively. The ratio between them is 7, while the ratio between the crosssectional area of both specimens is 5. This can signify that the load capacity does not depend only on the cross-sectional area, because the greater the thickness, the greater the stability.  Comparing the load capacity R of the specimens B-01 and B-03, they are 8.87 kN a 61.9 kN, respectively. The ratio between them is 7, while the ratio between the cross-s tional area of both specimens is 5. This can signify that the load capacity does not depe only on the cross-sectional area, because the greater the thickness, the greater the stabil Then, the stress , is computed by , = , ( Then, the stress S y,B is computed by and the equivalent elastic modulus E be obtained from where the ratio δ/Z represents average strain or unitary deformation in the vertical direction. The relationship between E and S y,B is represented in Figure 10. This figure also represents the deformation modulus E BCT as a function of obtained for boxes (see Table 3) and its corresponding correlation (12).
From these results, it can be pointed out that deformation modulus increases with stress. The higher the resistance, the higher the stiffness. The dispersion is very significant, but an average linear model can be fitted, yielding It can be concluded that the modulus for blocks E is higher than the one for boxes E BCT , because in the latter lateral deformations are given because the sheets of the boxes are slender.
stress. The higher the resistance, the higher the stiffness. The dispersion is very significant, but an average linear model can be fitted, yielding = 60.0 ⋅ . (13) It can be concluded that the modulus for blocks is higher than the one for boxes BCT , because in the latter lateral deformations are given because the sheets of the boxes are slender.

Discussion on Experimental Results
Generally, it should be pointed out that a significant dispersion has been found in experimental results. This is because the manufacturing process of these products does not provide significant repeatability, and this is the reason because in practical applications safety factors greater than two are employed if the load capacity is obtained from numerical models.
From the modeling point of view, Equations (1) and (9) are useful to define the range of use of packages. Therefore, if it must be determined the limit of boxes that can be stacked during transport or storage in a logistics center, the load would have to be divided by a safety factor. The deformation of the frame makes it possible to estimate whether the distance between the content and the cover plates is sufficient so that they do not come into contact, absorbing part of the load transmitted by the pallet overhead. These equations also allow us to obtain a simplified stress-deformation curve for engineering practical applications.

Discussion on Experimental Results
Generally, it should be pointed out that a significant dispersion has been found in experimental results. This is because the manufacturing process of these products does not provide significant repeatability, and this is the reason because in practical applications safety factors greater than two are employed if the load capacity is obtained from numerical models.
From the modeling point of view, Equations (1) and (9) are useful to define the range of use of packages. Therefore, if it must be determined the limit of boxes that can be stacked during transport or storage in a logistics center, the load would have to be divided by a safety factor. The deformation of the frame makes it possible to estimate whether the distance between the content and the cover plates is sufficient so that they do not come into contact, absorbing part of the load transmitted by the pallet overhead. These equations also allow us to obtain a simplified stress-deformation curve for engineering practical applications.
Equation (13) allows obtainment of a value for the modulus of elasticity of the material as if it were a homogeneous and isotropic solid. This simplifies the modeling of the boxes using a finite element model in industrial practical applications, since it is possible to consider the material as an elastic solid, as mentioned throughout the text. The lateral deformations that occur in boxes are considered by the numerical method itself, this is the reason for which the modulus E BCT of Equation (10) is not appropriate to model the material stiffness.

Finite Element Modeling
In this section, two CAE (Computer Aided Engineering) finite element models are developed aimed at verifying that simplifications about the elastic modulus of the material allows us to obtain simple models that provide results with enough accuracy. The results of the CAE models are compared with experimental results and the given by the mathematical models developed in Section 3.

Description of the Models
The models are developed by means of ABAQUS software. Two models are analyzed for the octabin O-04 and the rectangular R-04 items, respectively. The properties of these items are taken from Tables 1 and 3 and from the results of the Section 3. Specifically, Table 5 gathers ECT as a property of the material to be known from experimentation. With this, the stress S y,ECT is obtained by means of Equation (5). For the CAE, the modulus E is estimated by Equation (13). For the mathematical models, the ultimate stress S y,BCT is obtained from Equation (6) using f and m parameters, and the displacement at the maximum load δ max is determined by Equation (9). It can be assumed that the footprint perimeter p is much larger than the thickness h, hence S4R quadrilateral shell-type finite element is used. The O-04 model has 2200 elements and 2288 nodes. The R-04 model has 8232 elements and 8036 nodes.
The box is subjected to a vertical load F, which is uniformly distributed over its perimeter. The material model is assumed to be elastic and homogeneous with a uniform modulus of elasticity E. This is due to the hypothesis of this paper that the properties in the compression direction are the predominant ones. The Poisson ratio is taken 0.40 from [24].
The analysis is performed by means of the Riks method. Although the analysis could be carried out with the standard structural modulus of ABAQUS, the Riks method has been chosen because it performs arc-length control instead of force or displacement control, so it can even predict post-buckling behavior, such as the one shown in Section 4. The first 20 buckling modes are taken into consideration.
As a result, force-displacement curves are obtained for O-04 and R-04 models. To investigate the dependence of the mechanical behavior on the modulus E, and to take into account the dispersion of the results, three curves will be computed for each model: one for the reference value of Table 5, and the other two for the corresponding ±12%. Specifically, for O-04, the lowest and the higher values of the modulus are E = 104 N/mm 2 and E = 132 N/mm 2 , and those of the R-04 model are E = 95 N/mm 2 and E = 121 N/mm 2 . These curves are compared with the corresponding experimental ones and with the results of the mathematical models. Figure 11 shows a qualitative comparison of the deformed model with the equivalent tested box, previous to collapse. It is not possible to perform a quantitative comparison because in the test neither local stresses nor local displacements are measured. However, the tendencies are clearly represented by the model, in which the instability of vertical panels is evidenced. Figure 11a presents also in color the minimum principal stress field, corresponding to the compression direction. The zones with maximum lateral displacement due to buckling are subjected to the maximum compression stress. Thus, buckling is mainly governed by the stiffness properties in the compression direction.

Results for the Octabin-Type Box Models
The results of the maximum load and the corresponding vertical deformation are gathered in Table 6. From the CAE results, it can be concluded that the higher the modulus, the higher the resistance, and the higher the stiffness (i.e., the lower the deformation). This behavior can be also observed in the experimental curves. Specifically, ±12% variation of modulus implies approximately ±6% variation of the load capacity BCT and ∓10% variation of deformation. It can be also remarked that the mathematical method provides the lower BCT force and the lower deformation.   Figure 12 compares the force-displacement curves achieved by the CAE models with those of the mathematical model and the experimental results of the O-04 item. The experimental curves have been reproduced from those provided by the test machine software, represented under the legend. The mathematical model only gives the failure point, so the curve is represented as a straight line. Nevertheless, the CAE model takes into account the lateral deformation of the boxes, which implies that the stiffness is reduced with deformation, and the slope of the curve decreases with load. In addition, Riks method estimates post-buckling behavior also, as it can be appreciated in the three CAE curves. The results of the maximum load and the corresponding vertical deformation are gathered in Table 6. From the CAE results, it can be concluded that the higher the modulus, the higher the resistance, and the higher the stiffness (i.e., the lower the deformation). This behavior can be also observed in the experimental curves. Specifically, ±12% variation of modulus implies approximately ±6% variation of the load capacity BCT and ∓10% variation of deformation. It can be also remarked that the mathematical method provides the lower BCT force and the lower deformation.

Results for the Rectangular-Type Box Models
In the same way as in the previous section, a qualitative comparison of the CAE deformation model and the experimental one is presented in Figure 13. In the photography of Figure 13b, it can be observed that the failure of the box is given by the instability by local buckling of the shorter lateral panel, that of 785 mm (the side b of Table 1). This has been also revealed in the simulation: the shorter panel presents the maximum lateral deformation, and it is subjected to the maximum compression stress (in blue), inducing buckling. On the contrary, the larger panel is curved by bending, but it has not buckled due to the higher resistance to buckling, because it is less slender. From the simulation, it can be also concluded that the corners are subjected to the maximum compression stress, because they are local rigid zones.   Table 7. From a qualitative point of view, similar conclusions as for the O-04 item are driven. Quantitatively, CAE results indicate that ±12% variation of modulus implies approximately ±5% variation of the load capacity BCT and only ∓5% variation of deformation. It should be noted that the mathematical method gives the lowest deformation, and the BCT is similar to experimental and to the stiffest CAE ones.    Table 7. From a qualitative point of view, similar conclusions as for the O-04 item are driven. Quantitatively, CAE results indicate that ±12% variation of modulus implies approximately ±5% variation of the load capacity BCT and only ∓5% variation of deformation. It should be noted that the mathematical method gives the lowest deformation, and the BCT is similar to experimental and to the stiffest CAE ones.

Conclusions
In this paper, a characterization of the mechanical behavior of corrugated cardboard applied to containers and octabins subjected to vertical stacking loads has been presented. Experiments have been completed in order to obtain mathematical models for load capacity of the boxes, for vertical deformation and for deformation modulus. From the experimental results, it can be pointed out the significant existing dispersion, inherent to the characteristics of this kind of products. Thus, the obtained mathematical models are functional in practical engineering applications design stages, but safety factors have to be applied for the final application.
In order to employ simple finite element models to simulate the corrugated cardboard boxes behavior, an equivalent Young modulus has been drawn from experimental tests on corrugated cardboard billets. This allows modeling of boxes by means of shell finite elements with a homogeneous equivalent material. Simulations have been performed with different properties to take into account the dispersion. Simulation results have been compared with experimental and mathematical model results. Taking into account the existing dispersion in experimental results, it can be concluded that CAE models based on homogeneous and isotropic materials is an effective way to easily simulate the behavior of containers made of corrugated cardboard, although the material is actually orthotropic. This is due to the fact that these types of containers fail due to buckling of the vertical panels, so the properties in the compression direction are the most significant.
Author Contributions: J.G., E.A. and F.C. writing, conceptualization and investigation; A.G. writing, conceptualization and funding acquisition. All authors have read and agreed to the published version of the manuscript.

Conclusions
In this paper, a characterization of the mechanical behavior of corrugated cardboard applied to containers and octabins subjected to vertical stacking loads has been presented. Experiments have been completed in order to obtain mathematical models for load capacity of the boxes, for vertical deformation and for deformation modulus. From the experimental results, it can be pointed out the significant existing dispersion, inherent to the characteristics of this kind of products. Thus, the obtained mathematical models are functional in practical engineering applications design stages, but safety factors have to be applied for the final application.
In order to employ simple finite element models to simulate the corrugated cardboard boxes behavior, an equivalent Young modulus has been drawn from experimental tests on corrugated cardboard billets. This allows modeling of boxes by means of shell finite elements with a homogeneous equivalent material. Simulations have been performed with different properties to take into account the dispersion. Simulation results have been compared with experimental and mathematical model results. Taking into account the existing dispersion in experimental results, it can be concluded that CAE models based on homogeneous and isotropic materials is an effective way to easily simulate the behavior of containers made of corrugated cardboard, although the material is actually orthotropic. This is due to the fact that these types of containers fail due to buckling of the vertical panels, so the properties in the compression direction are the most significant.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.