Estimation of the Compressive Strength of Corrugated Cardboard Boxes with Various Openings

: This paper presents mixed analytical/numerical method for estimating the static top-to-bottom compressive strength of corrugated packaging with different ventilation openings and holes, in which the torsional and shear stiffness of corrugated cardboard as well as the panel depth-to-width ratio are included. Analytical framework bases on Heimerls assumption with a modiﬁcation to a critical force, which is here computed by a numerical algorithm. The proposed method is compared herein with the successful McKee formula and is veriﬁed with the large number of experiment results of various packaging designs made of different qualities of corrugated cardboard. The results show that, for various hole dimensions or location of openings in no-ﬂap and ﬂap boxes, the estimation error may be reduced up to three times than in the simple analytical approach.


Introduction
The packaging industry plays an important role in the modern world. The leading material used in this field is a corrugated cardboard. Corrugated boxes are used in many branches of various industries, mainly to transport goods or products from one place to another. During the transportation and storage, the boxes are subject to considerable loading. In order to ensure appropriate strength of the boxes, it is necessary to use appropriate estimation method and/or to conduct a series of experiments to determine a behavior of the corrugated material and the whole construction of the box. The strength estimation becomes challenging when box contains, e.g., hand holes and/or openings. Correct characterization of material properties of the corrugated board and precise estimation technique of the box strength allows to design the carton boxes in line with the sustainable development strategy. In the material such as corrugated board, there are two characteristic in-plane directions of orthotropy related to fiber orientation. The first one is perpendicular to the main axis of the corrugation and it is called the machine direction (MD). It is parallel to the paperboard fiber alignment. The second one is parallel to it and called cross direction (CD).
One possibility for analysis of strength of the corrugated boxes is the execution of several iterative experiments, in the loop in which testing, validation and optimization of the box design take place. In the packaging industry the most commonly used are compressive, bending or bursting strength tests of corrugated board. Another important material property to examine is also a humidity, which has a big impact on the box strength. However, the most important and practical are the box compression test (BCT) and the edge crush test (ECT).
Another possibility is to use analytical formulas for prediction compressive strength of boxes made from corrugated boards available in the literature. The formulas should be as simple as possible to make it useful for practical applications. Over the years, many different approaches has been proposed. The parameters used in these formulas can be categorized intro three groups: paper, board and box parameters [1]. In the first group, one can distinguish two tests-ring crush test (RCT), Concora liner test (CLT), liner type, weights of liner and fluting, corrugation ratio and a constant related to fluting. The parameters of board include thickness, flexural stiffnesses in MD and CD, ECT and moisture content. The parameters of boxes are dimensions, perimeter, applied load ratio, stacking time, buckling ratio and printed ratio. In 1952, Kellicutt and Landt used the parameters of paper (RCT, flute constant) and the box (perimeter, box constant) to propose one of the first formula for prediction of compressive strength of boxes [2]. Maltenfort related the critical force in the BCT to the parameters of the paper (CLT, liner type) and box dimensions [3]. The most popular and widely used formula has been proposed by McKee, Gander and Wachuta in 1963 [4]. They took into account the parameters of the paperboard (ECT, flexural stiffnesses) and the box perimeter. The formula is applicable for some relatively simple containers. Over the years, many researchers tried to extend the applicability of the McKee's approach. Probably one of the first is the work of Allerby et al. [5] from 1985, in which the authors modified constant and exponents of the McKee's formula. Two years later, Schrampfer et al. extended applicability of the McKee relationship for a wider range of cutting methods and equipment [6]. Batelka et al. included in this relationship also the dimensions of the box [7], while Urbanik et al. took into account also the Poisson's ratio [8]. In the original paper of McKee et al. [4], an arbitrary chosen constant value was used, which makes their approach limited only for some typical boxes. This constant was later analyzed in more complicated cases by Garbowski et al. [9]. If one would like to include also the transverse shear, the additional tests are required and therefore updated formulation should be used, which was recently modified by Aviles et al. [10] and later by Garbowski et al. [11,12].
On the contrary, in the analysis of corrugated boxes numerical methods are very often used to compute the strength of the boxes. The most common method is the finite element method (FEM), which was successfully applied in many research on numerical analysis of boxes made from corrugated cardboards, including investigations of buckling phenomena [13] or transverse shear stiffness of corrugated cardboards [14,15].
In numerical simulations, the corrugated paperboards are commonly modelled as sandwich panels [16]. The main goal in the numerical analysis is to obtain results as accurate as possible with computational costs as small as possible. Thus, there is a need of using models, which are simplification of the real objects without losing adequacy to reality. In the numerical analysis of layered materials [17][18][19] and constructions, made from such materials, it is useful to perform a process called homogenization [20,21]. As a consequence, instead of structure of layers made from different materials one can obtain one layer characterized by effective properties of the composite. Application of the external loads to the homogenized plate should give similar effect like for the composite material with heterogeneous core. One of the homogenization methods, which was based on strain energy and can be applied for sandwich panels, was proposed by Hohe [22]. In this approach, an equivalence of a representative element of the heterogeneous and homogenized elements was assumed. A period homogenization method, in which an equivalent membrane and pure bending characteristic of period plates is obtained, has been also proposed by Buannic et al. [23]. They applied a modified version of this approach to take into account the transfer shear effect in the analysis of sandwich panels. Biancolini used the FEM for analysis of a micromechanical part of the considered plate [24]. Appling the energy equivalency between the model and the equivalent plate, he obtained the stiffness properties of the considered plate. Abbès and Guo decomposed the plates into two beams in directions of the plate to obtain the torsion rigidity of the orthotropic sandwich plates [25]. Two approaches to homogenization have been compared by Marek and Garbowski. These considerations were based on the classical laminated plate theory and the deformation energy equivalence method [26] or based on inverse analysis [27]. The strength of the boxes made from corrugated cardboards can be reduced by many factors, including moisture, prints, perforations, etc. However, one of the most common factor for loss of the strength of boxes is the presence of openings on the walls of the package, as most of the boxes contain ventilation or at least hand holes [28]. This issue was a subject of some recently published papers. Fadji et al. [29] investigated an effect of the geometrical features of the holes, including height, area, orientation, shape and number of holes on the strength of boxes. The influence of ventilation holes on the stability and strength of the boxes using FEM was examined by Fadiji et al. [30]. They considered certain type of boxes, which are commonly used in the fresh fruit industry, in which the ventilation holes play an important role.
To the best knowledge of the authors, there is no research in which the effect of losing strength caused by holes is taken into account in the analytical prediction of the box strength. In this paper, we extend the applicability of the approach presented in [9] for corrugated boxes with holes.

Ultimate Compressive Strength of a Single Panel
The most popular and the simplest approach used these days for computation of ultimate compressive load of corrugated cardboard boxes was proposed by McKee et al. [4]. The starting point in their study was the determination of the ultimate plate load P f of a rectangular panel with dimensions a × b, loaded vertically on the upper edge, in which the plate under consideration is a separated panel of a box (Figure 1a), despite the fact that the box is loaded as a three-dimensional structure during the compression test ( Figure 1b). All four edges of the each separated plate are pinned in the out-of-plane direction, and the bottom edge is additionally constrained along b in the vertical direction. The basis for the derivation given by McKee et al. is the empirical formula presented in Heimerl [31], where the load P f at failure [N/mm] is the result of a combination of buckling force and compressive strength, namely: where r is an exponent, r ∈ (0, 1), k is a certain constant, P cr is a critical load of a corrugated panel [N/mm] and ECT is the edge-loaded compressive strength of a corrugated board in cross direction [N/mm]. These considerations were based on the classical laminated plate theory and the deformation energy equivalence method [26] or based on inverse analysis [27]. The strength of the boxes made from corrugated cardboards can be reduced by many factors, including moisture, prints, perforations, etc. However, one of the most common factor for loss of the strength of boxes is the presence of openings on the walls of the package, as most of the boxes contain ventilation or at least hand holes [28]. This issue was a subject of some recently published papers. Fadji et al. [29] investigated an effect of the geometrical features of the holes, including height, area, orientation, shape and number of holes on the strength of boxes. The influence of ventilation holes on the stability and strength of the boxes using FEM was examined by Fadiji et al. [30]. They considered certain type of boxes, which are commonly used in the fresh fruit industry, in which the ventilation holes play an important role.
To the best knowledge of the authors, there is no research in which the effect of losing strength caused by holes is taken into account in the analytical prediction of the box strength. In this paper, we extend the applicability of the approach presented in [9] for corrugated boxes with holes.

Ultimate Compressive Strength of a Single Panel
The most popular and the simplest approach used these days for computation of ultimate compressive load of corrugated cardboard boxes was proposed by McKee et al. [4]. The starting point in their study was the determination of the ultimate plate load of a rectangular panel with dimensions × , loaded vertically on the upper edge, in which the plate under consideration is a separated panel of a box (Figure 1a), despite the fact that the box is loaded as a three-dimensional structure during the compression test ( Figure 1b). All four edges of the each separated plate are pinned in the out-of-plane direction, and the bottom edge is additionally constrained along in the vertical direction. The basis for the derivation given by McKee et al. is the empirical formula presented in Heimerl [31], where the load at failure [N/mm] is the result of a combination of buckling force and compressive strength, namely: where is an exponent, ∈ (0,1), is a certain constant, is a critical load of a corru-   (b) cardboard packaging in press during box compression testing; (c) box compression press used in laboratory tests [32]. To include the nonlinear material effect when no elastic buckling occurs the above formula was later modified by Urbanik et al. [8,13]: where n = 1 − r. Integer u plays an activating role and if it is equal 1, it corresponds to elastic buckling (ECT > P cr ), while if it is equal 0, it corresponds to inelastic buckling, so the ultimate load P f is proportional to ECT, namely: In both cases, if the selected panel of a corrugated box does not have any openings, see Figure 1a, the analytical formula describing the critical load for a rectangular orthotropic plate [9,[33][34][35] can be utilized: where P b cr is the critical load of the panel with dimensions a × b (a-panel height, b-panel width [mm]), k cr is a dimensionless buckling coefficient that depends, among others, on the ratio a/b, boundary conditions applied to the plate edges, the buckling mode and a material characteristics, D 11 is the effective bending stiffness in the MD [N/mm] and D 22 is the effective bending stiffness in the CD [N/mm]. In this study, the first buckling mode was used for computations in all cases. In this buckling mode, the displacements of panel with and without hole usually represent a shape of the half-wave in both directions.
All fibrous materials, including corrugated cardboard, are characterized by orthotropy, which means that the mechanical properties in such materials change depending on the direction, here CD/MD. The buckling coefficient for such panels takes the following form: where m defines a buckling mode and counts the half-waves along the load direction, D 33 is an effective torsion stiffness [N/mm] and D 12 = D 11 ν 21 = D 22 ν 12 , ν 12 and ν 21 are in-plane Poisson's ratios.

Buckling-McKee's Formula
If the analyzed panel does not contain any holes, a simplified formula for critical load computation can be utilized. Such a formula was suggested by McKee et al. [4], in order to maximally simplify the complicated equations describing the critical load of compressed orthotropic plates so that the use of these equations could be common and practical. Therefore, the authors replaced Equation (5) with a simpler formula: where θ = D −1 11 D 22 and K is a plate parameter assumed to be equal to 0.5, which leads to: They also assumed that the plate has equal sides, i.e., b = c and therefore b = Z/4 (where Z is the box perimeter), which led to further simplification of Equation (5) to the form: Energies 2021, 14, 155 5 of 20

Buckling-Finite Element Method
However, if the analyzed panel contains any openings, the analytical approach becomes impractical due to the lack of an universal solution for a panel with different holes. In such a case, the numerical method works best, although it is more difficult to implement, it is very general and allows to determine the critical load of panels with any shape and boundary conditions, including holes with complex shapes. In this study, the buckling was computed for each separated panel in the 3D space having all external edges fixed in out of plane direction, additionally the bottom edge was fixed in the vertical direction, other BC's were inactive.
The critical load can be calculated from the formula: where λ is a critical loading multiplier, i.e., the smallest eigenvalue and q is a distributed load on a panel upper edge [N/mm].
To calculate the critical force multiplier, the following two-step algorithm should be used. In the first step, we need to: • discretize the panel with finite elements (FE), see Figure 2a, • compute the stiffness matrix K e of each element, • assemble the global stiffness matrix K 0 of a whole panel, • compute nodal forces representing initial loading configuration P * , i.e., for loading multiplier λ = 1 (one-parameter loading assumed P = λP * ), • take boundary conditions into account, • solve equation 0 d * = P * to obtain nodal displacements in pre-buckling state: • extract element displacements d (e) * (from displacements of the system d * ) and compute in each element the generalized stresses s (e) * .
In the second step, we need to: • generate initial geometrical stress matrices for each element K e σ (s * e ) and the whole panel K σ (s * ), • formulate initial buckling problem: • solve the eigenproblem to determine the pairs ( Energies 2021, 14, x FOR PEER REVIEW 5 of 20

Buckling-Finite Element Method
However, if the analyzed panel contains any openings, the analytical approach becomes impractical due to the lack of an universal solution for a panel with different holes. In such a case, the numerical method works best, although it is more difficult to implement, it is very general and allows to determine the critical load of panels with any shape and boundary conditions, including holes with complex shapes. In this study, the buckling was computed for each separated panel in the 3D space having all external edges fixed in out of plane direction, additionally the bottom edge was fixed in the vertical direction, other BC's were inactive.
The critical load can be calculated from the formula: where ̅ is a critical loading multiplier, i.e., the smallest eigenvalue and is a distributed load on a panel upper edge [N/mm].
To calculate the critical force multiplier, the following two-step algorithm should be used. In the first step, we need to: • discretize the panel with finite elements (FE), see Figure 2a, • compute the stiffness matrix of each element, • assemble the global stiffness matrix of a whole panel, • compute nodal forces representing initial loading configuration * , i.e., for loading multiplier = 1 (one-parameter loading assumed = * ), • take boundary conditions into account, • solve equation * = * to obtain nodal displacements in pre-buckling state: * = · * , • extract element displacements ( ) * (from displacements of the system * ) and compute in each element the generalized stresses ( ) * .
In the second step, we need to: • generate initial geometrical stress matrices for each element ( * ) and the whole panel ( * ), • formulate initial buckling problem: • solve the eigenproblem to determine the pairs ( , ), … , ( , ), where is a number of degrees of freedom, is th eigenvalue, = ∆ is th eigenvector (post-buckling deformation mode).  In order to compute stiffness matrix of the plate, the proper formulation of the FE must be used. In the classical plate FE, the field variables are approximated through shape functions according to the associated node values as follows: where N j (x, y) is a shape function associated with node j and N n is a number of nodes in the element: is the displacement vector in the j-th node. The membrane, bending and shear strains associated to the displacement in Equation (11) can be therefore obtained as follows: where: The energy functional of the plate can be finally obtained as follows: where D and D s are stiffnesses of the plates [9]. By using Lagrange's equations for the Equation (19), the characteristic equation of the system is obtained as follows: where: and: As the geometry of each panel does not have to be regular, especially in the case of panels with openings, it was assumed that the FE mesh will be made of shell triangular elements. Most Reissner-Mindlin triangular shell FE described in the literature were developed using the assumed transverse shear strain techniques. Here, the 6-noded quadratic triangle FE with assumed linear shear strains was utilized. Zienkiewicz et al. [36] developed such a 6-noded plate triangle with a standard quadratic interpolation for the rotations and the deflection, and the assumed linear transverse shear strain field in the following form: i.e.,: Figure 2b,c show the position of the six shear sampling points and the ξ i directions. Following the procedure presented in Appendix A, we find: where: Matrix B s is computed by Equation (A24) in Appendix A. A full 3-point Gauss quadrature is used here for the numerical integration of element transverse shear stiffness matrix: This element performs well for most of applications, although sometimes is too flexible [36,37]. Improvements can be introduced by imposing a linear variation of the normal rotation θ n along all three triangle sides as: where i = 4, 5, 6 are the mid-side nodes (Figure 2a). Equation (31) can be later used for eliminating the normal rotation at the midsize nodes [36,38].

Box Compression Strength-McKee's Formula
To compute failure load of a single panel we substitute the buckling force derived from Equation (7) into Equation (1), after [4] we have: which can be simplified to: where k 1 , under the assumption [4] that k 1−r cr is constant if a/Z > 1/7 and equals approximately 1.33, while r = 0.746 and k = 0.4215, is equal to: To calculate the failure load of the entire box, it is enough to multiply Equation (33) by the Z, which is the perimeter of the box. The original assumption [4] was that the width and length of a box are equal (i.e., b = c) thus the compressive strength of a box, known as the long McKee formula, reads: The empirical observations [4] that √ D 11 D 22 ∼ = 66.1 ECT h 2 , led to a further simplification and finally to the most known short form of McKee's equation: where k 2 , with previously assumed value of r = 0.746, equals: So, Equation (36) takes the final form: This is the most commonly used formula. However, its accuracy is acceptable only for the simplest boxes. The McKee formula should not be used if the length to width ratio or the height to length ratio of the box is too large [9]. The same problem appear when box contains any openings or ventilation holes. To overcome this limitations we have derived an analytical formulation for box strength from Equation (1), in which the critical force is computed using numerical methods instead of analytical formulas. Both will be shown in the following sections.

Box Compression Strength-General Case
In order to calculate the compressive strength of a single panel of a box, see Figure 1a, the force P f described by the formula (1) should be multiplied by the i-th panel width b or c. To get the failure load of an entire box (BCT), it is required to sum up the compressive strength of all panels. Thus, in the general case, we get the following: where (·) i indicates the i-th panel of width b or c (i = 1, . . . , 2). P b cr and P c cr are the critical forces of panels of width b and c, respectively (see Figure 1a). γ b p and γ c p are the reduction factors taking into account the ratio of the compressive strength of the plate with and without a hole for plate width b and c, respectively. γ b and γ c are the reduction coefficients due to in-plate aspect ratio of the box dimensions, defined as:

Practical Considerations
Equation (39) uses different definitions of the critical force, i.e., Equations (4) and (8), depending on whether the i-th panel contains any hole or not. In any case, the values of constants k and r in Equation (39) are required. As suggested by other authors [4,[39][40][41], k = 0.5 and r = 0.75 were assumed.
The compressive strength of rectangular boxes with various openings (with base dimensions b × c, as shown in Figure 1a) after taking the typical values of k and r in Equation (39) takes the form: where P b cr and P c cr are defined in: (a) Equation (4)-for orthotropic rectangular plates without holes (analytical), (b) Equation (8)-for orthotropic rectangular plates with various openings (numerical -FEM).
In all of the above equations describing the critical force P c cr , all definitions of width b should be replaced by c, e.g., Equation (4) should read:

Box Compression Strength-Experiment vs. Estimation
In the following section, experimental and computational examples regarding boxes are presented and the applicability of the analytical-numerical approach is validated. In the paper, the analyzed cases represent the simplest slotted-type boxes with holes varied in size, position at the panel side and location on the longer or shorter wall. Two typical designs of FEFCO boxes with an addition of different holes were considered, i.e., FEFCO F200-packaging without upper flaps, see Figure 3, and FEFCO F201-packaging with flaps, see Figure 4. Two box dimensions were considered, namely, 300 × 200 × 200 mm and 300 × 200 × 300 mm (length, width and height, respectively). Moreover, three positions of the opening at the panel side were assumed, i.e., in the upper corner (O1), in the middle of the package height at the edge of the panel (O2) and in the center of the panel (O3). Also, two cases were considered, in which the opening was at the shorter wall (B), and at the longer wall (L) of the box. In most cases, the openings were square, but rectangular openings were also considered.   For the purposes of this study, the selected boxes were tested in a box compression testing machine, namely BCT-19T10 from FEMat [32]; see Figure 1c. According to the specifications, the machine measures the force up to 10 kN with a resolution of 0.1 N. The displacement accuracy is 0.001 mm. In order to get reliable results, 4-6 samples of each packaging design were tested. In the study, for all cases the number of tested samples  For the purposes of this study, the selected boxes were tested in a box compression testing machine, namely BCT-19T10 from FEMat [32]; see Figure 1c. According to the specifications, the machine measures the force up to 10 kN with a resolution of 0.1 N. The displacement accuracy is 0.001 mm. In order to get reliable results, 4-6 samples of each packaging design were tested. In the study, for all cases the number of tested samples For the purposes of this study, the selected boxes were tested in a box compression testing machine, namely BCT-19T10 from FEMat [32]; see Figure 1c. According to the specifications, the machine measures the force up to 10 kN with a resolution of 0.1 N. The displacement accuracy is 0.001 mm. In order to get reliable results, 4-6 samples of each packaging design were tested. In the study, for all cases the number of tested samples exceeded 180. The boxes were manufactured from different three-layer corrugated cardboard, namely, 350, 380, 400, and 480 g/m 2 . All those cardboards had B or E flute, their thickness was 1.49, 1.52, 1.59, 2.80 and 2.82 mm. All cases tested are presented in Tables 1 and 2 for FEFCO type boxes of F200 and F201, respectively. In Tables 1 and 2 "opening type" columns contain the information about the opening, for instance O1B1S, means that the opening was located in the corner (O1) of the shorter wall (B). "Cardboard quality" columns contain the code of cardboard, for instance 3B400-1, means that three layer cardboard and B flute was used with 400 g/m 2 , 1st material production. During the tests, the boxes were folded by one operator manually, and the flaps of the packaging were taped at the bottom and top, if applicable. All packaging structures were cut on a plotter, so the material was not converted. Also, each box was inspected for visual damage of the cardboard before it was submitted to displacement control protocol tests. TAPPI standard T402 was applied, the samples were laboratory conditioned, i.e., temperature 23 • C ± 1 • C and relative humidity of 50% ± 2% was hold. Figure 5 demonstrates examples of laboratory measurements for selected designs of packaging with holes, namely F200-O1B1S, F201-O3B1S1c and F201-O3L1S1d. The curves of force versus displacement are shown; a good repeatability of the measurement results is observed. The samples were not inspected after the measurements in order to determine the deformation or damage mechanism. visual damage of the cardboard before it was submitted to displacement control protocol tests. TAPPI standard T402 was applied, the samples were laboratory conditioned, i.e., temperature 23°C ± 1°C and relative humidity of 50% ± 2% was hold. Figure 5 demonstrates examples of laboratory measurements for selected designs of packaging with holes, namely F200-O1B1S, F201-O3B1S1c and F201-O3L1S1d. The curves of force versus displacement are shown; a good repeatability of the measurement results is observed. The samples were not inspected after the measurements in order to determine the deformation or damage mechanism.
(a) (b) (c) Moreover, the measured strength of the boxes with openings selected in this study serve as a validation data for those obtained with analytical-numerical approach proposed here, see Equation (39). In order to use Equation (39) the critical buckling loads of box panels (with and without openings) were required. The material data used for computing critical loads by FEM are presented in Table 3. The cardboard qualities were also presented in Table 1 and Table 2. All input data for Equation (39) used for the computations were attached in Supplementary Materials, separately for F200 and F201. Table 3. Material data for the cardboard qualities used in FEM computations of critical loads [42,43]. In Figure 6, the data from the tests were presented graphically (black plot), also separately for F200 ( Figure 6a) and F201 (Figure 6b) boxes. Table 1, Table 2 and Figure 6 show the mean values of the ultimate load measured at the box failure (BCT). In Figure 6, also the values obtained from the analytical-numerical approach proposed were presented (blue plot), see Equation (39). In the analytical-numerical approach presented, and used are the typical values taken from the literature, namely = 0.5 and = 0.75, see Equation (42). Moreover, the measured strength of the boxes with openings selected in this study serve as a validation data for those obtained with analytical-numerical approach proposed here, see Equation (39). In order to use Equation (39) the critical buckling loads of box panels (with and without openings) were required. The material data used for computing critical loads by FEM are presented in Table 3. The cardboard qualities were also presented in Tables 1 and 2. All input data for Equation (39) used for the computations were attached in Supplementary Materials, separately for F200 and F201. Table 3. Material data for the cardboard qualities used in FEM computations of critical loads [42,43]. In Figure 6, the data from the tests were presented graphically (black plot), also separately for F200 ( Figure 6a) and F201 (Figure 6b) boxes. Table 1, Table 2 and Figure 6 show the mean values of the ultimate load measured at the box failure (BCT). In Figure 6, also the values obtained from the analytical-numerical approach proposed were presented (blue plot), see Equation (39). In the analytical-numerical approach presented, k and r used are the typical values taken from the literature, namely k = 0.5 and r = 0.75, see Equation (42). The error of the method proposed may be calculated by the comparing the computed values with the ultimate loads measured, namely:

No. Cardboard Quality
where is an experimental average value of the box compression strength for a bo with a particular geometry, cardboard quality used and the type of opening and is its counterpart obtained via analytical-numerical approach proposed here. The mean error obtained for approach proposed here with typical parameters of and wa to 14.9%.
Also, similar computations for all boxes were performed by using the McKee for mula. The results are presented in Figure 6 (dashed magenta plot), see Equation (38). In similar way, the mean error of McKee short formula was computed, which was equal t 15.5%.

Reduction of the Estimation Error-Optimal Parameters
To determine how much the estimation error of the proposed method depends on the parameters ( , ), the computations for different sets of those parameters were per formed for the experimental data obtained in the course of this research. The values of and , see Equation (39), were modified and the mean error was computed for particula pair of those parameters. Notice that the typical values of and from the literature ar 0.5 and 0.75, respectively, see Equation (42). The mean errors obtained by systematic com putations are presented by contour plot in Figure 7. The contour plot values were com puted by averaging the magnitudes from Equation (44) for all box designs. The interva of and was the same for both parameters, namely, from 0.3 to 1.0. The error of the method proposed may be calculated by the comparing the computed values with the ultimate loads measured, namely: where BCT exp is an experimental average value of the box compression strength for a box with a particular geometry, cardboard quality used and the type of opening and BCT est is its counterpart obtained via analytical-numerical approach proposed here. The mean error obtained for approach proposed here with typical parameters of k and r was equal to 14.9%. Also, similar computations for all boxes were performed by using the McKee formula. The results are presented in Figure 6 (dashed magenta plot), see Equation (38). In a similar way, the mean error of McKee short formula was computed, which was equal to 15.5%.

Reduction of the Estimation Error-Optimal Parameters
To determine how much the estimation error of the proposed method depends on the parameters (k, r), the computations for different sets of those parameters were performed for the experimental data obtained in the course of this research. The values of k and r, see Equation (39), were modified and the mean error was computed for particular pair of those parameters. Notice that the typical values of k and r from the literature are 0.5 and 0.75, respectively, see Equation (42).  Table 4 summarizes the mean errors computed for typical set of parameters ( = 0.5 and = 0.75) and for the optimal parameter set obtained via optimization ( = 0.755 and = 0.435). The optimal solution was obtained by the systematic search in the − space, see contour plot in Figure 7. Optimal set of and parameters used in the Equation (39) gives the lowest computed mean error. The mean estimation errors using typical and optimal set of parameters were equal 14.9% (red dot in Figure 7) and 6.5% (black cross in Figure 7), respectively. Note that in Figure 7 the surface has a characteristic valley, in which optimal value was found. The value of mean error obtained by the use of the optimal parameters in the proposed analytical-numerical method is about three times smaller than the error obtained by the use of typical values of and . Having the optimal values of parameters ( , ), see Table 4, the analytical solutions may be revalidated for considered designs of boxes. The results are presented in Figure 6, see red plot. As expected, the value obtained from optimization are closer to the values measured than those computed for typical and .

Discussion
Given these results, the conclusion can be drawn that the best performance of the computational approach was obtained by the proposed analytical-numerical method with the optimal parameters. If the typical parameter were used the error was much bigger. The simple McKee formula gave almost three times bigger estimation error than for the method proposed with optimal parameters. Using optimally selected parameters of and greatly reduces the error in estimating the strength of boxes with openings considered in the study.
Notice that for the method proposed the difference in modelling boxes of FEFCO F200 and F201 was not observed, i.e., the errors are moderately similar. The mean error for optimal parameters for F200 boxes was equal 6.8%, while for F201 boxes was equal 5.2%. Considering the box strengths, see Figures 6a and 6b, one can notice that the greater grammage of the cardboard, the obtained value of the box strength is greater. However,  Table 4 summarizes the mean errors computed for typical set of parameters (k = 0.5 and r = 0.75) and for the optimal parameter set obtained via optimization (k = 0.755 and r = 0.435). The optimal solution was obtained by the systematic search in the k − r space, see contour plot in Figure 7. Optimal set of k and r parameters used in the Equation (39) gives the lowest computed mean error. The mean estimation errors using typical and optimal set of parameters were equal 14.9% (red dot in Figure 7) and 6.5% (black cross in Figure 7), respectively. Note that in Figure 7 the surface has a characteristic valley, in which optimal value was found. The value of mean error obtained by the use of the optimal parameters in the proposed analytical-numerical method is about three times smaller than the error obtained by the use of typical values of k and r. Having the optimal values of parameters (k, r), see Table 4, the analytical solutions may be revalidated for considered designs of boxes. The results are presented in Figure 6, see red plot. As expected, the value obtained from optimization are closer to the values measured than those computed for typical k and r.

Discussion
Given these results, the conclusion can be drawn that the best performance of the computational approach was obtained by the proposed analytical-numerical method with the optimal parameters. If the typical parameter were used the error was much bigger. The simple McKee formula gave almost three times bigger estimation error than for the method proposed with optimal parameters. Using optimally selected parameters of k and r greatly reduces the error in estimating the strength of boxes with openings considered in the study.
Notice that for the method proposed the difference in modelling boxes of FEFCO F200 and F201 was not observed, i.e., the errors are moderately similar. The mean error for optimal parameters for F200 boxes was equal 6.8%, while for F201 boxes was equal 5.2%. Considering the box strengths, see Figure 6a,b, one can notice that the greater grammage of the cardboard, the obtained value of the box strength is greater. However, in the cases 24-31 (480 g/m 2 ) the calculated values are lower than for the cases 18-23 (400 g/m 2 ), this is caused by the different papers (predominance of recycled or virgin fibers) used for the production of these corrugated boards, which is noticeable through measured material properties, see Table 3. Comparing the parameters of the first three cardboard qualities: 3B400A1-1, 3B400A1-3 and 3B480-1, one can observe that the values are very close to each other and the ECT values, used in the original McKee formula, are even greater for the cardboards of lower grammage.
In general, the box strength is greater if the hole is smaller and its location is closer to the center of the panel. It can be noticed in Figure 6b for the cases 24-26, which are ordered according to the increasing size of hole for the same cardboard quality and location of the hole-on the shorter wall (B) in the center of panel (O3). The only exception can be observed comparing the cases 27-31, in which experimentally for the greater hole (case 29) the greater box strength was obtained than for smaller hole (case 28). This issue is suspected to be caused by a initial imperfections of the corrugated board used for the production of these boxes.
It should be underlined that for particular cardboard quality, the McKee formula gives constant estimation of box strength. McKee formula results are independent from the opening type (its position, shape or size). For instance, see the results for boxes from no. 24 up to 31 in Figure 6b, which are presented for 3B480-1 cardboard quality. In the study, not only various design of boxes with holes were tested, also the cardboard with different quality was considered and packaging characterized with different overall dimensions of width, length and height. Such systematic design of experiments allows to use those data as a part of the input for numerical training, e.g., the artificial network. Using artificial intelligence algorithms could be especially beneficial for optimal selection of various structure parameters (not only two parameters selected in this study), if other crucial effects for reliable strength estimation of boxes would be determined. Among others, the structural parameters may describe the position of the opening, its area or shape. etc. This could be the interest of further research.

Conclusions
The strength of regular flap boxes is easily determined with a use of a simple McKee formula, however, if the box dimensions are not in standard proportions, the McKee formula becomes inaccurate, as shown in our previous study. If the box contains any openings the use of pure analytical formulas become even more impractical, because the loss of accuracy is very noticeable. The article presents the analytical formula feed with numerically obtained data for estimation of the strength of boxes with openings of various dimensions and locations.
The proposed method was verified on the basis of experimental data, which covered various examples of cardboard packaging. In the paper, the tested cases represent the simplest boxes with holes varied in size, position at the panel side and position at the longer or shorter wall. Also, two design of boxes with different cardboard qualities were considered, namely, modified design of FEFCO F200 and FEFCO 201. Qualities of three layer cardboards with E and B flute were used.
Due to the use of the new calculating approach proposed here, the estimation error was on average reduced three times compared to the simplified McKee analytical method. Proposed approach seems to be very promising in considering more complicated cases of boxes with complex geometries or/and materials. Thus, future studies will be devoted to adapting the analytical-numerical approach proposed here to estimate the strength of boxes with perforations of different properties. The study may concern shelf-ready boxes, in which the perforation of different cutting patterns may be used. Properties of the perforations could be also selected in the artificial intelligence algorithms as mentioned structural parameters, which among other factors, determines/decreases the strength of the box. Especially that new innovative methods of perforating are currently being tested in the corrugated board industry.

Derivation of the Substitute Transverse Shear Strain Matrix
An isoperimetric Reissner-Mindlin plate element with n nodes is assumed here, for which the rotation, the vertical deflection and the also transverse shear strains are independently interpolated according to the following equations: ω = N ω w (e) , θ = N θ θ (e) , γ = N γ γ (e) . (A1) We also assume that this interpolation satisfies the conditions: In the first step, the transverse shear strains are interpolated in the natural coordinate system ξ, η as: γ = γ ξ γ η = 1 ξ η ξη · · · ξ p η q | 0 0 0 · · · 0 0 0 0 0 · · · 0 | 1 ξ η · · · ξ r η s where n γ specifies the number of polynomial expansion sampling points for γ ξ and γ η in the element. The transverse shear strains in the cartesian system are then obtained as: where J is the in-plane Jacobian matrix [37,44,45]: The transverse shear along a natural direction ξ i strain is defined as: where β i is the angle between the ξ i direction and the natural ξ axis. The shear strains γ ξ i are sampled along the ξ i directions at each of the n γ points. From Equation (A6) we obtain the following: where: where (·) i , i = 1, . . . , n γ etc. denote transverse shear values at each sampling point. From Equation (A3) it is found that:γ Substituting Equation (A10) into (A7) gives: where P = TÂ is a n γ × n γ matrix. Vector α is obtained as Combining Equations (A3), (A7) and (A12) gives: The relationship between the Cartesian and the natural transverse shear strains at each sampling point is:γ where J i is the in-plane Jacobian matrix at the ith sampling point. Thus: Substituting Equations (A13) and (A15) into (A4) gives: where: are the shape functions. The relationship between the transverse shear strains in cartesian coordinate system and the nodal displacements in each element can be written by the weighted residual procedure [37,46] as: where ∇ = ∂ ∂x , ∂ ∂y T and W k are arbitrarily chosen weighting functions, e.g., Galerkin weighting, i.e., W k = N γ k . Substituting Equations. (A16) and the equation above into (A18) [37] gives:  Equation (A19) is then used to obtain the shear strains at the sampling points for each element as:γ =B s a (e) , (A20) If the point collocation method is used in Equation (A18), i.e., W i = δ i where δ i is the Dirac delta at the ith sampling point (i = 1, . . . , n γ ), we obtain: where B i s is the value of the original transverse shear strain matrix at the ith sampling point, so finally: γ = J −1 A P −1 T CB s a (e) = B s a (e) , where B s is the substitute transverse shear strain matrix given by: