Analytical Determination of the Bending Stiffness of a Five-Layer Corrugated Cardboard with Imperfections

Bending stiffness (BS) is one of the two most important mechanical parameters of corrugated board. The second is edge crush resistance (ECT). Both are used in many analytical formulas to assess the load capacity of corrugated cardboard packaging. Therefore, the correct determination of bending stiffness is crucial in the design of corrugated board structures. This paper focuses on the analytical determination of BS based on the known parameters of the constituent papers and the geometry of the corrugated layers. The work analyzes in detail the dependence of the bending stiffness of an asymmetric, five-layer corrugated cardboard on the sample arrangement. A specimen bent so that the layers on the lower wave side are compressed has approximately 10% higher stiffness value. This is due to imperfections, which are particularly important in the case of compression of very thin liners. The study showed that imperfection at the level of a few microns causes noticeable drops in bending stiffness. The method has also been validated by means of experimental data from the literature and simple numerical finite element model (FEM). The obtained compliance of the computational model with the experimental model is very satisfactory. The work also included a critical discussion of the already published data and observations of other scientists in the field.


Introduction
A sign of the present times is the constant pursuit of the purchase of various merchandise, and thus the need for their packaging and safe transport, both in traditional forms of sale and e-commerce. The foremost desirable characteristics of the packaging are naturally adequate strength in relation to light weight and, in the interest of the natural environment, reusable, recyclable and biodegradable. Corrugated cardboard packaging perfectly meets all the above-mentioned requirements. Going further, the popularity of this type of packages is associated with the intensive development of a separate branch of industry and research. In view of the laws governing the free market, manufacturers strive for the most cost-effective solutions while maintaining the appropriate load-bearing capacity of cardboard packaging. The scientists, who have been developing for many years new methods to determine the material properties of the corrugated cardboard [1,2] and are constantly trying to understand the nature of the packaging performance, through numerous studies, involving a variety of techniques, are here to help. The task is challenging mainly due to the layered structure of the corrugated cardboard with two characteristic in-plane directions of orthotropy associated with the mechanical strength of the paperboard-the machine direction (MD) perpendicular to the main axis of the fluting and parallel to the paperboard fiber alignment, and cross direction (CD) which is parallel to the fluting. Moreover, there are a number of factors that reduce the strength of a cardboard itself or corrugated cardboard packaging, the impact of which has been analyzed and is still is the subject of investigation, e.g., [3] in particular time and storing conditions [4,5], stacking load [6][7][8], respect to mechanical properties of liners and the flute geometric parameters was conducted to answer the question of which of the parameters have the greatest impact on BS. The main contribution of this study was the derivation of analytical relationships that explain the differences in the bending stiffness of asymmetrical corrugated boards when the layers on the B or E flute side are compressed.

The Four-Point Bending Test of a Sample with an Asymmetric Cross-Section
In the case of four-point bending of a sample with an asymmetric cross-section, especially when its cross-section consists of thin-walled faces (an example of such board is the corrugated cardboard), the dependence of the bending stiffness on the direction of the moment can be noticed. Using theoretical models as well as linear numerical models, this effect cannot be capture. This phenomenon belongs to the imperfection class of problems.
Since in the four-point bending test, the mid-segment is bent with a constant moment (see Figure 1) and all other section forces are not present, therefore, the problem is greatly simplified. In the case of asymmetrical sections in such test, the sample can be placed and, consequently, examined in two positions, which results in different values of the determined bending stiffness. Thin-walled structures, when they undergo bending (i.e., one part of the cross-section is compressed and other part is in tension), have higher bending stiffness if the "stronger" part of the cross-section is compressed (see Figure 2a). On the other hand, when the "weaker" part of the cross-section is compressed (see Figure 2b), the BS is lower-this is due to the preliminary buckling of the compressed fragments of the thinwalled cross-section. More information and a short discussion on this phenomenon can be found in the following subsections. In our case, where a five-layer corrugated cardboard sample is tested, the stronger part of the cross-section is on the E-wave side. Therefore, from now on, the following description is used to distinguish two cases: (a) EB-compression of a part of the cross-section on the B wave side (see Figure 2b) and (b) BE-compression of a part of the cross- Thin-walled structures, when they undergo bending (i.e., one part of the cross-section is compressed and other part is in tension), have higher bending stiffness if the "stronger" part of the cross-section is compressed (see Figure 2a). On the other hand, when the "weaker" part of the cross-section is compressed (see Figure 2b), the BS is lower-this is due to the preliminary buckling of the compressed fragments of the thin-walled cross-section. More information and a short discussion on this phenomenon can be found in the following subsections. respect to mechanical properties of liners and the flute geometric parameters was conducted to answer the question of which of the parameters have the greatest impact on BS.
The main contribution of this study was the derivation of analytical relationships that explain the differences in the bending stiffness of asymmetrical corrugated boards when the layers on the B or E flute side are compressed.

The Four-Point Bending Test of a Sample with an Asymmetric Cross-Section
In the case of four-point bending of a sample with an asymmetric cross-section, especially when its cross-section consists of thin-walled faces (an example of such board is the corrugated cardboard), the dependence of the bending stiffness on the direction of the moment can be noticed. Using theoretical models as well as linear numerical models, this effect cannot be capture. This phenomenon belongs to the imperfection class of problems.
Since in the four-point bending test, the mid-segment is bent with a constant moment (see Figure 1) and all other section forces are not present, therefore, the problem is greatly simplified. In the case of asymmetrical sections in such test, the sample can be placed and, consequently, examined in two positions, which results in different values of the determined bending stiffness. Thin-walled structures, when they undergo bending (i.e., one part of the cross-section is compressed and other part is in tension), have higher bending stiffness if the "stronger" part of the cross-section is compressed (see Figure 2a). On the other hand, when the "weaker" part of the cross-section is compressed (see Figure 2b), the BS is lower-this is due to the preliminary buckling of the compressed fragments of the thinwalled cross-section. More information and a short discussion on this phenomenon can be found in the following subsections. In our case, where a five-layer corrugated cardboard sample is tested, the stronger part of the cross-section is on the E-wave side. Therefore, from now on, the following description is used to distinguish two cases: (a) EB-compression of a part of the cross-section on the B wave side (see Figure 2b) and (b) BE-compression of a part of the cross- In our case, where a five-layer corrugated cardboard sample is tested, the stronger part of the cross-section is on the E-wave side. Therefore, from now on, the following description is used to distinguish two cases: (a) EB-compression of a part of the cross-section on the B wave side (see Figure 2b) and (b) BE-compression of a part of the cross-section on the E wave side (see Figure 2a). In the BE configuration, higher BS values are obtained.

Corrugated Cardboard-Samples
This study uses the results of the research presented in the work by Czechowski et al. [29]. The authors presented the mechanical parameters of the component papers (cor-rugated and flat layers; see Table 1), the geometric parameters of the five-layer corrugated cardboard (5EB; see Table 2) and the results of four-point bending tests for six boards made of various combinations of component papers. In this paper, the above experimental data were the starting point for in-depth studies on the cause of the difference in the bending stiffness of samples bent with a positive and negative moment. The geometric characteristics of the corrugated layers are also shown in Figure 3. It was assumed that the shape of wavy layers is described by a trigonometric function with the amplitude h i and the period 2π/p i . section on the E wave side (see Figure 2a). In the BE configuration, higher BS values are obtained.

Corrugated Cardboard-Samples
This study uses the results of the research presented in the work by Czechowski et al. [29]. The authors presented the mechanical parameters of the component papers (corrugated and flat layers; see Table 1), the geometric parameters of the five-layer corrugated cardboard (5EB; see Table 2) and the results of four-point bending tests for six boards made of various combinations of component papers. In this paper, the above experimental data were the starting point for in-depth studies on the cause of the difference in the bending stiffness of samples bent with a positive and negative moment. The geometric characteristics of the corrugated layers are also shown in Figure 3. It was assumed that the shape of wavy layers is described by a trigonometric function with the amplitude ℎ and the period 2 / . Since in the adopted calculation model (details will be discussed in the next subsection), only flat layers affect the machine direction (MD) bending stiffness, therefore only the MD stiffness moduli for liners only for all six boards are listed in Table 1. The corrugated layers play the role of keeping the liners at the right distance to ensure adequate bending stiffness. The geometry of the separated, undulating layers is presented in Table 2.  Since in the adopted calculation model (details will be discussed in the next subsection), only flat layers affect the machine direction (MD) bending stiffness, therefore only the MD stiffness moduli for liners only for all six boards are listed in Table 1.
The corrugated layers play the role of keeping the liners at the right distance to ensure adequate bending stiffness. The geometry of the separated, undulating layers is presented in Table 2.
In the geometrical description of the corrugated board, however, the thickness of the individual layers should also be taken into account, see Figure 4.
Therefore, the total height of the corrugated cardboard 5EB is: where h * i are the distances between the central axes of the liners. So the corrected E-flute height is: where ℎ * are the distances between the central axes of the liners. So the corrected E-flute height is: while the corrected B-flute height is: (3) Figure 4. Thickness of all corrugated board layers and axial distance between liners in the computational model. Table 3 summarizes the thicknesses of all corrugated cardboard layers for six boards and the calculated distances between liners (i.e., corrected heights of the corrugated layers), according to Equations (2) and (3). As already mentioned, the measured thicknesses of the constituent papers (see Table  3) and the geometry of the corrugated layers (see Table 2) were taken from [29].

Bending Stiffness of Assymetric Corrugated Board with Imperfections
It is assumed in this study that only liners are involved in bending in the MD, which means that fluting has only the role of supporting the liners in the correct position. Therefore, their tensile/compression and bending stiffnesses are negligible. In order to derive the model, first, all liners were segmented between the supporting wave crests in the corrugated layer (see Figure 5). Preliminary geometric imperfections were included in the compressed segments, which reduced the stiffness of these elements. These assumptions allow to capture the difference in BS depending on the sign of a bending moment in the asymmetric boards.  Table 3 summarizes the thicknesses of all corrugated cardboard layers for six boards and the calculated distances between liners (i.e., corrected heights of the corrugated layers), according to Equations (2) and (3). As already mentioned, the measured thicknesses of the constituent papers (see Table 3) and the geometry of the corrugated layers (see Table 2) were taken from [29].

Bending Stiffness of Assymetric Corrugated Board with Imperfections
It is assumed in this study that only liners are involved in bending in the MD, which means that fluting has only the role of supporting the liners in the correct position. Therefore, their tensile/compression and bending stiffnesses are negligible. In order to derive the model, first, all liners were segmented between the supporting wave crests in the corrugated layer (see Figure 5). Preliminary geometric imperfections were included in the compressed segments, which reduced the stiffness of these elements. These assumptions allow to capture the difference in BS depending on the sign of a bending moment in the asymmetric boards.
Since the five-layer corrugated board consists of three liners, three different segment lengths can be distinguished: L 1 and L 3 correspond to the E and B wave periods, respectively. On the other hand, the length L 2 can reach the maximum value of p 1 (where p 1 is a lower wave period, see Figure 3). However, usually every second segment is divided into two sections by the crest of wave B, indicated by L 2 on Figure 5. Therefore, the average length equal to 2/3 p 1 was adopted in the middle liner for further analyses. Since the five-layer corrugated board consists of three liners, three different segmen lengths can be distinguished: and correspond to the E and B wave periods, respec tively. On the other hand, the length can reach the maximum value of (where is a lower wave period, see Figure 3). However, usually every second segment is divided into two sections by the crest of wave B, indicated by ' on Figure 5. Therefore, the av erage length equal to 2/3 was adopted in the middle liner for further analyses. For a compressed i-th segment (see Figure 6) with a geometric imperfection, the lon gitudinal shortening of i-th segment, can be computed by the classical equation: where: = is the normal force; = 1 is a virtual normal force; = is the cross section area (with -segment width and -segment thickness); = 12 ⁄ is th cross-section moment of inertia; = is the bending moment; and = 1 is a virtual bending moment. By inserting all the relationships described above into Equa tion (4) and taking all the constant values out of the integral, the longitudinal deflection takes the form: where the deflection function can be described by e.g., the sine function: with the maximum deflection in the middle of the element span (see Figure 6) which was assumed here as a small fraction of the element length: = ⋅ 10 whil the value of is a quantity assumed between 2 and 4. When measuring the bending stiffness of five-layer corrugated board one can plac a corrugated board sample with the E wave facing upwards or vice versa. As the result in one case, two liners are compressed on the E wave side (see Figure 7) while in the othe case, a single liner on the B wave side (see Figure 8) is compressed. In a case where two liners on the E-flute side are compressed, a higher value of bending stiffness is usually obtained. This is because larger number and shorter (therefore, less slender) segments ar For a compressed i-th segment (see Figure 6) with a geometric imperfection, the longitudinal shortening of i-th segment, δ i can be computed by the classical equation: where: N i = P i is the normal force; N i = 1 is a virtual normal force; A i = bt i is the crosssection area (with b-segment width and t-segment thickness); I i = bt 3 i /12 is the crosssection moment of inertia; M i = P i w i (x) is the bending moment; and M i = 1 w i (x) is a virtual bending moment. By inserting all the relationships described above into Equation (4) and taking all the constant values out of the integral, the longitudinal deflection takes the form: where the deflection function w i can be described by e.g., the sine function: with the maximum deflection f i in the middle of the element span L i (see Figure 6), which was assumed here as a small fraction of the element length: f i = L i · 10 −k while the value of k is a quantity assumed between 2 and 4.

Figure 5. Compressed and stretched segments of a corrugated board cross-section during ben
Since the five-layer corrugated board consists of three liners, three different segm lengths can be distinguished: and correspond to the E and B wave periods, res tively. On the other hand, the length can reach the maximum value of (wher is a lower wave period, see Figure 3). However, usually every second segment is div into two sections by the crest of wave B, indicated by ' on Figure 5. Therefore, th erage length equal to 2/3 was adopted in the middle liner for further analyses. For a compressed i-th segment (see Figure 6) with a geometric imperfection, the gitudinal shortening of i-th segment, can be computed by the classical equation: where: = is the normal force; = 1 is a virtual normal force; = is the c section area (with -segment width and -segment thickness); = 12 ⁄ is cross-section moment of inertia; = is the bending moment; and = 1 is a virtual bending moment. By inserting all the relationships described above into E tion (4) and taking all the constant values out of the integral, the longitudinal defle takes the form: where the deflection function can be described by e.g., the sine function: with the maximum deflection in the middle of the element span (see Figur which was assumed here as a small fraction of the element length: = ⋅ 10 w the value of is a quantity assumed between 2 and 4. When measuring the bending stiffness of five-layer corrugated board one can p a corrugated board sample with the E wave facing upwards or vice versa. As the re in one case, two liners are compressed on the E wave side (see Figure 7) while in the o case, a single liner on the B wave side (see Figure 8) is compressed. In a case where liners on the E-flute side are compressed, a higher value of bending stiffness is usu obtained. This is because larger number and shorter (therefore, less slender) segment compressed and even when imperfections are present, the effect is less pronounced. When measuring the bending stiffness of five-layer corrugated board one can place a corrugated board sample with the E wave facing upwards or vice versa. As the result, in one case, two liners are compressed on the E wave side (see Figure 7) while in the other case, a single liner on the B wave side (see Figure 8) is compressed. In a case where two liners on the E-flute side are compressed, a higher value of bending stiffness is usually obtained. This is because larger number and shorter (therefore, less slender) segments are compressed and even when imperfections are present, the effect is less pronounced.  In a four-point bending test, only the bending moment occurs in the center of the specimen, so the model simplifies to pure bending. In this model, the moment is balanced by the normal forces acting in the liners on the arms (see Figure 7) with respect to the neutral axis : The starting point for determining the bending stiffness is the kinematic excitation in the form of rotation (see Figure 7), which causes elongation or contraction of liners (see Figure 8). By taking a small value (e.g., 10 mm) and using the known values of , the remaining values can be determined (see Figure 8) and finally the rotation angle can be obtained: By solving the integrals in Equation (5), it is possible to determine the longitudinal deflection in liners under compression or tension: For the known deflection , the compressive force in the i-th liner can be determined: while the tensile force (for = 0) is: The i-th bending moment is:  In a four-point bending test, only the bending moment occurs in the center of the specimen, so the model simplifies to pure bending. In this model, the moment is balanced by the normal forces acting in the liners on the arms (see Figure 7) with respect to the neutral axis : The starting point for determining the bending stiffness is the kinematic excitation in the form of rotation (see Figure 7), which causes elongation or contraction of liners (see Figure 8). By taking a small value (e.g., 10 mm) and using the known values of , the remaining values can be determined (see Figure 8) and finally the rotation angle can be obtained: By solving the integrals in Equation (5), it is possible to determine the longitudinal deflection in liners under compression or tension: For the known deflection , the compressive force in the i-th liner can be determined: while the tensile force (for = 0) is: The i-th bending moment is: In a four-point bending test, only the bending moment occurs in the center of the specimen, so the model simplifies to pure bending. In this model, the moment is balanced by the normal forces P i acting in the liners on the arms z i (see Figure 7) with respect to the neutral axis z 0 : The starting point for determining the bending stiffness is the kinematic excitation in the form of rotation φ (see Figure 7), which causes elongation or contraction δ i of liners (see Figure 8). By taking a small value δ 1 (e.g., 10 −2 mm) and using the known values of z i , the remaining values δ i can be determined (see Figure 8) and finally the rotation angle φ can be obtained: By solving the integrals in Equation (5), it is possible to determine the longitudinal deflection in liners under compression or tension: For the known deflection δ i , the compressive force P i in the i-th liner can be determined: while the tensile force P i (for f i = 0) is: The i-th bending moment is: and the bending stiffness can be calculated as the sum of the integrals from the formula: If, instead of bending, the corrugated board cross-section is compressed or stretched in MD, the equation for compression/tensile stiffness can also be derived: The theoretical bending stiffness (valid for a perfect model without imperfections) as the product of the stiffness modulus in MD and the moment of inertia of liners only is: In order to normalize the theoretical values and the results of four-point bending tests, both values are divided by the width of the sample b. Hence, ultimately, the bending stiffness is: The presented derivation allows to explain the differences between the bending stiffness obtained from testing of the corrugated board sample placed with the E wave upwards or vice versa-B wave.

Results
In the first step, the theoretical assumption, in which only liners affect the stiffness of the entire section, was validated. For this purpose, two simple numerical models of a five-layer corrugated cardboard in a plane state (i.e., a beam model) were built (see Figure 9). Both models consist of classic Bernoulli 2-node beam elements and were implemented in Matlab software (Mathworks Inc., Natick, MA, USA) [65]. Small rotation φ was applied in both ends and the corresponding reaction moments M were determined in order to calculate BS from Equations (13) and (16). In all cases, displacements resulting from φ rotation wrt neutral axis were applied on both ends of the model (in external nodes on the left and right sides of the model).

=
, (12) and the bending stiffness can be calculated as the sum of the integrals from the formula: If, instead of bending, the corrugated board cross-section is compressed or stretched in MD, the equation for compression/tensile stiffness can also be derived: The theoretical bending stiffness (valid for a perfect model without imperfections) as the product of the stiffness modulus in MD and the moment of inertia of liners only is: In order to normalize the theoretical values and the results of four-point bending tests, both values are divided by the width of the sample . Hence, ultimately, the bending stiffness is: The presented derivation allows to explain the differences between the bending stiffness obtained from testing of the corrugated board sample placed with the E wave upwards or vice versa-B wave.

Results
In the first step, the theoretical assumption, in which only liners affect the stiffness of the entire section, was validated. For this purpose, two simple numerical models of a fivelayer corrugated cardboard in a plane state (i.e., a beam model) were built (see Figure 9). Both models consist of classic Bernoulli 2-node beam elements and were implemented in Matlab software (Mathworks Inc.) [65]. Small rotation was applied in both ends and the corresponding reaction moments were determined in order to calculate BS from Equations (13) and (16). In all cases, displacements resulting from rotation wrt neutral axis were applied on both ends of the model (in external nodes on the left and right sides of the model). In the first model all layers were modeled according to their geometry and mechanical parameters, while in the second model, the stiffness of the corrugated layers was significantly reduced (by 100 times) to mimic a situation where only liners are active. The results are shown in Figure 10. Naturally, this assumption is not valid if one would like to derive the BS in CD, where all liners as well as both corrugated layers are equally important. In the first model all layers were modeled according to their geometry and mechanical parameters, while in the second model, the stiffness of the corrugated layers was significantly reduced (by 100 times) to mimic a situation where only liners are active. The results are shown in Figure 10. Naturally, this assumption is not valid if one would like to derive the BS in CD, where all liners as well as both corrugated layers are equally important.
In order to eliminate a possible error related to the discretization of numerical models, the influence of the number of finite elements and the number of waves in the model was also checked. The results are summarized in Table 4. All FE models consist of 2-node linear beam elements with a seed equal to 0.1 mm, which generated the following number of nodes and elements in four models: In order to eliminate a possible error related to the discretization of numerical models, the influence of the number of finite elements and the number of waves in the model was also checked. The results are summarized in Table 4. All FE models consist of 2-node linear beam elements with a seed equal to 0.1 mm, which generated the following number of nodes and elements in four models:   In the next step, the influence of imperfection amount on the bending stiffness in the analytical model was analyzed. The results for the parameter ranging from 2 to 4 are shown in Figure 11. The selected value of = 2.3 is marked on all graphs along with corresponding BS values for both case EB and BE. The selected value of gives the best agreement between the results obtained with the proposed model and the available experimental data.   In the next step, the influence of imperfection amount on the bending stiffness in the analytical model was analyzed. The results for the parameter k ranging from 2 to 4 are shown in Figure 11. The selected value of k = 2.3 is marked on all graphs along with corresponding BS values for both case EB and BE. The selected value of k gives the best agreement between the results obtained with the proposed model and the available experimental data.
Because the presented analytical model takes into account the initial imperfections of compressed segments in the corrugated board, thus allows to distinguish between the bending stiffness of the corrugated board whether the E wave or the B wave is compressed. The bending stiffness not only decreases with the increase of the initial imperfection, but also the BS difference between the EB and BE increases as the imperfections increase (see Figure 12). In other words, as the initial imperfections increase, the bending stiffness of the sample in the EB configuration (compression on the B wave side) decreases faster than the bending stiffness of the sample in the BE configuration. Because the presented analytical model takes into account the initial imperfections of compressed segments in the corrugated board, thus allows to distinguish between the bending stiffness of the corrugated board whether the E wave or the B wave is compressed. The bending stiffness not only decreases with the increase of the initial imperfection, but also the BS difference between the EB and BE increases as the imperfections increase (see Figure 12). In other words, as the initial imperfections increase, the bending stiffness of the sample in the EB configuration (compression on the B wave side) decreases faster than the bending stiffness of the sample in the BE configuration. The difference in bending stiffness between the EB and BE configurations can be as high as 25%. However, this applies to cross-sections in which the imperfection amounts to 1% of the initial length of the compressed segment, . In our case, the initial imperfection, for the selected value of the k-factor, is 0.5% of . The practically zero difference between EB and BE case can be observed for the coefficient equals to 3, i.e., initial imperfection equals to 0.1% of . Table 5 gathers all bending stiffness values for all 6 boards determined from experimental data [29], theoretical model (no imperfections), simplified FEM model (2D beam model-no imperfections), full 3D shell FE model [29]-no imperfections, and proposed here analytical model with imperfections. The difference in bending stiffness between the EB and BE configurations can be as high as 25%. However, this applies to cross-sections in which the imperfection amounts to 1% of the initial length of the compressed segment, L i . In our case, the initial imperfection, for the selected value of the k-factor, is 0.5% of L i . The practically zero difference between EB and BE case can be observed for the coefficient k equals to 3, i.e., initial imperfection equals to 0.1% of L i . Table 5 gathers all bending stiffness values for all 6 boards determined from experimental data [29], theoretical model (no imperfections), simplified FEM model (2D beam model-no imperfections), full 3D shell FE model [29]-no imperfections, and proposed here analytical model with imperfections. As the differences in the results summarized in Table 5, especially the differences between the experimental measurements and all computational models, suggest some errors in the experimental data, the sensitivity analysis was performed in the last step. This analysis was to show which of the parameters have the greatest impact on BS and therefore to point out which measurements require careful re-checking in order to find possible inaccuracies in experimental data presented in [29]. Figure 13 presents all sensitivities of BS in two configurations: EB and BE with respect to mechanical properties of corrugated board and the flute geometric parameters.
All graphs in Figure 13 show the BS sensitivity to 10% perturbations of (a) thickness of all corrugated cardboard layers, (b) liners stiffness moduli as well as (c) E and B wave heights. All other parameters do not affect the bending stiffness in both wave orientation (E wave up or B wave up). Certainly, the shape of the corrugated layer (apart from the amplitude) has no effect on BS because, as already proved in this paper, the flute itself contributes less than 1% to overall bending stiffness of corrugated cardboard.
Materials 2022, 15, x FOR PEER REVIEW 12 of 17 As the differences in the results summarized in Table 5, especially the differences between the experimental measurements and all computational models, suggest some errors in the experimental data, the sensitivity analysis was performed in the last step. This analysis was to show which of the parameters have the greatest impact on BS and therefore to point out which measurements require careful re-checking in order to find possible inaccuracies in experimental data presented in [29]. Figure 13 presents all sensitivities of BS in two configurations: EB and BE with respect to mechanical properties of corrugated board and the flute geometric parameters.

Discussion
The results presented in the study were obtained while using derived analytical or numerical models, in which the experimental data presented in [29] were utilized. All experimental data used in the work are summarized in Tables 1 and 2. The heights of the corrugated layers (E-flute and B-flute) have been corrected and are compiled in Table 3. To the best of our knowledge, there are no other studies in the literature on this subject, although several observations made by various scientific groups have already indicated this phenomenon, e.g., Östlund and Niskanen [66]. The proposed analytical model does not take into account the corrugated layers in the calculation of the bending stiffness of the five-layer corrugated cardboard. Two numerical models were therefore built to validate this assumption. In the first of them, both flat and corrugated layers were used to determine BS. In the second one, the stiffness of the wavy layers was reduced to emulate a situation in which corrugated layers are excluded from the computation. For the comparison, the theoretical model was additionally employed, in which also just flat layers are considered in BS computation. All three models gave almost the same results presented in Figure 10. The differences were between 0.72% and 1.75%.
Since the presented model takes into account the influence of initial imperfections in the compressed segment of the corrugated board, the first attention was focused on determining the initial imperfection values. Figures 11 and 12 present BS in two configurations: (a) EB-E wave upwards and (b) BE-B wave upwards, as a function of imperfection value.
In Figure 11 it can be clearly noticed that for the imperfection value at the level of 0.1% of the initial length of the compressed segments (which corresponds to k equal to 3), not only the difference between the EB and BE configurations is not noticeable, but also the difference between the BS values for both EB and EB do not differ from the reference value (dashed lines) computed while using the theoretical model (see Equation (15)). The difference between the bending stiffness in the case of EB and BE increases with the augmentation of the imperfection coefficient and for the value k = 2 (i.e., imperfection equals to L i · 10 −2 ) it is between 12% and about 22% (see Figure 12). In this work the assumed imperfection coefficient is 2.3, which corresponds to the initial imperfections in the compressed elements at the level of 0.5% of the initial length of these segments. Table 5 summarizes all calculated and literature values of bending stiffness for six examples of five-layer corrugated board. It is clearly seen than in just two cases the theoretical BS is higher than the experimentally measured BS. It can be evidently noticed that only in two cases (Board 2 and Board 6) the theoretical BS is higher than experimentally measured BS. This is an alarming observation, because in the case of real structures made of corrugated board, the cross-section is rarely ideal (usually the corrugated board is slightly crushed [67,68]), which means that the measured bending stiffness values should rather be lower than theoretical. Not only the theoretical values of BS are lower than those measured experimentally. Virtually all the results presented in Table 5 follow a similar trend, both the results obtained with the use of analytical and numerical models, including the results from the literature (column 6) [29]. Due to this observation, the results of experimental research presented in [29] may contain some errors or are incorrectly ordered. Despite these doubts the results obtained while using the analytical model are very good for Examples 2 and 6 (marked in Table 6), for other Examples the results are not as good but still better than results presented in [29] (see Table 6). The mean absolute error generated by the analytical model is 11.7% for all cases while the mean absolute error of the results presented in [29] is 16.4%. Due to the relatively large discrepancies between the calculated and measured values of bending stiffness, and due to the suspected measurement error or incorrect compilation of results in [29], the sensitivity analysis of the analytical model was also carried out in this study. The graphs shown in Figure 13 clearly indicate that both the EB and BE models have the greatest sensitivity to the change in the stiffness modulus and thickness of the flat inner and outer layers (i.e., Liner-1 and Liner-3). The sensitivity of BS to changes in the height of the corrugated layer B, h 2 , is similarly high. Thus, even a small change of these parameters (just a few percent), can dramatically change the computational value of bending stiffness of the corrugated board.

Conclusions
In this study, a detailed analysis of the effect of imperfections in thin-walled asymmetrical sections bent with a constant moment was carried out. The main contribution of this work was the derivation of analytical relationships that accurately describe the phenomenon of the difference in bending stiffness depending on the sign of the moment loading the asymmetric corrugated cardboard sample in machine direction. The paper showed that the applied analytical model satisfactorily reflects the real behavior of bent fivelayer corrugated cardboard. The adopted simplifications did not affect the quality of the proposed solution, which was proved by a simple numerical model. Finally, the developed model was compared with the results of experimental research available in the literature. The obtained results are much closer to the experimental results than the results generated by other models available in the literature. Additionally, proposed model is very easy to implement, which makes it possible to use it in practice by cardboard manufacturers. This study also includes the sensitivity analysis, which indicates the most important parameters directly affecting the BS and, therefore, can be very helpful in more conscious design of optimal corrugated board.