Influence of Imperfections on the Effective Stiffness of Multilayer Corrugated Board

There are many possible sources of potential geometrical inaccuracies in each layer of corrugated board during its manufacture. These include, among others, the processes of wetting the corrugated layers during profiling, the process of accelerated drying, the gluing process, and any mechanical impact of the pressure rollers on the cardboard. Work taking into account all the above effects in numerical modeling is not well described in the literature. Therefore, this article presents a simple and practical procedure that allows us to easily account for geometric imperfections in the calculation of the effective stiffness of corrugated board. As a main tool, the numerical homogenization based on the finite element method (FE) was used here. In the proposed procedure, a 3D model of a representative volumetric element (RVE) of a corrugated board is first built. The numerical model can include all kinds of geometrical imperfections and is used to calculate the equivalent tensile and bending stiffnesses. These imperfections were included in the 3D numerical model by appropriate modeling of individual layers, taking into account their distorted shape, which was obtained on the basis of a priori buckling analysis. This paper analyzes different types of buckling in order to find the most representative one. The proposed procedure is easy to implement and fully scalable.


Introduction
Currently, the growing consumption and the increase in the number of individual shipments make it necessary to ensure the safe storage and transport of goods. At the same time, more emphasis is placed on ecology, which is why companies quickly began to abandon plastic, which not only harms the environment, but also has a negative impact on the reputation of enterprises. For these reasons, corrugated cardboard has begun to conquer the packaging market. It owes its huge popularity to numerous advantages that largely meet the needs of today. Cardboard packaging is recyclable, easy to dispose of, biodegradable, durable under appropriate conditions, and takes up little space before folding. In addition to the positives for the environment, there are also useful properties for companies. Corrugated cardboard products can be easily shaped by adding holes, ventilation holes, or by printing a brand logo on them. It is also possible to create boxes with perforations, which are often used in the food industry. Such packages are opened by tearing off a part of the material along the previously designed perforations, and then immediately put on the shelf. Such a procedure allows for significant time savings, which brings profits for large companies.
In addition to all the above-mentioned advantages, packaging should, above all, effectively protect the goods. For this reason, the load-bearing capacity of packaging has been studied by scientists since the 1950s. The first methods were analytical methods determined for boxes with a rectangular base. In 1952, Kellicutt and Landt introduced an approach based on box perimeter, overall ring crush strength, and paper and box constants [1]. Four years later, Maltenfort presented the box compressive strength depending on the critical force, paper parameters, box dimensions, and empirical constants [2]. In 1963, McKee et al. proposed a relationship that has been commonly used in the packaging industry for many years to estimate the compressive strength of cardboard boxes [3]. The popular McKee formula is based on the edge crush resistance (ECT) value, flexural stiffnesses, box perimeter, and correction factors. Due to the simplicity of the formula, it is possible to obtain a result quickly, but it can only be used for simple standard boxes. In the following years, scientists tried to extend and modify the McKee formula. In 1985, Allerby et al. changed the constants and exponents [4] and, in 1987, Schrampfer et al. extended the approach to new cardboard cutting methods and equipment [5]. In 1993, Batelka and Smith included the dimensions of the packaging in the formula [6].
In practice, the compressive strength of the packaging depends on many factors related to the construction and storage conditions [7]. From the construction point of view, attention should be paid to aspects, such as openings, ventilation holes, perforations, and offsets [8][9][10][11][12], because they significantly affect the behavior of the box under pressure. Load-bearing capacity is also influenced by moisture content [13,14], stacking load [15], storage time and conditions [16], and many other factors.
With the development of technology, numerical computations have become popular. The well-known finite element method (FEM) is often used to determine the compressive strength of packaging. Urbanik and Frank compared the compressive strength of boxes calculated using FEM with the results of the McKee formula extended by Poisson's ratio [17]. Furthermore, FEM simulations were used by Nordstrand and Carlsson to compare analytical, experimental, and numerical values of transverse shear moduli [18,19]. The issue of transverse shear was studied by Aviles et al. [20] and Garbowski et al. [21,22], where the role of transversal shear stiffness in orthotropic sandwich material is presented. Urbanik and Saliklis used FEM to observe buckling and post-buckling phenomena in cardboard boxes [23]. Maneengam et al. used the FE model to study the vibration and damping characteristics of honeycomb sandwich structures reinforced with carbon nanotubes [24]. Corrugated board reinforced with metal liners has also been investigated by Gu et al., where uniaxial compression tests were conducted to analyze the damage and failure modes [25]. Sohrabpour and Hellström presented a review of analytical and numerical methods for estimating box compressive strength [26].
Due to the orthotropy of the paper and layered structure, corrugated board is quite difficult to analyze numerically. To facilitate the computations, the homogenization process is used, which consists of replacing a complicated cardboard cross-section with a homogeneous board with equivalent parameters. Homogenization methods can be divided into analytical and numerical. In analytical approaches, the equations of the classical theory of material strength and the classical theory of laminates are used [27]. Numerical homogenization is based on a finite element method framework. In this study, a method based on strain energy equivalence between a 3D cardboard model and a flat plate is used. This homogenization procedure was introduced by Biancolini [28] and later extended by Garbowski and Gajewski to include transversal shear stiffnesses [29]. In 2003, Hohe presented a representative element of heterogeneous and homogenized structures derived from the strain energy equations [30]. Homogenization was also used to determine substitute panel parameters, such as membrane and bending characteristics [31] and torsional stiffness [32]. In 2018, a multiscale asymptotic homogenization was used by Ramírez-Torres et al. for layered hierarchical structure analysis [33,34]. Suarez et al. used numerical homogenization to analyze seating made of multi-wall corrugated cardboard [35]. Garbowski et al. presented the influence of different creasing and perforation on cardboard parameters [36]. In 2022, Mrówczyński et al. presented non-local sensitivity analysis in the optimal design of three- [37] and five-layered [38] corrugated cardboard. A review of homogenization meth-ods and their accuracy was made by Nguyen-Minh et al., based on analysis of corrugated panels [39].
During the production process, the corrugated board may be deformed due to changes in temperature and humidity. Two types of emerging imperfections can be distinguished, namely global and local. A model describing systematic, large-scale deviation from the intended flat shape of cardboard was presented by Beck and Fischerauer [40]. However, more attention was paid to local imperfections. In 1995, Nordstrand presented the effect of the size of imperfections on the compressive strength of cardboard boxes [41], and in 2004 extended the nonlinear buckling analysis of Rhodes and Harvey orthotropic plates to include local imperfections [42]. Lu et al. investigated the effect of imperfections on the mechanical properties of cardboard during compression [43]. The problem of non-ideal shape during bending of double-walled corrugated cardboard was analyzed analytically by Garbowski and Knitter-Piątkowska [44]. Mrówczyński et al. presented a method of including the initial imperfections in the analysis of single-walled cardboard [45]. Recently Cillie and Coetzee presented an experimental and numerical study on the in-plane compression of corrugated paperboard panels with global and local imperfections [46].
In the article, the authors focused on the issue of local imperfections in double-walled corrugated board. The method presented here allows us to quickly and easily obtain the effective cardboard stiffnesses reduction due to the imperfect shape of its component layers. The described approach is an extension of the homogenization method proposed in the authors' previous works. Due to the specificity of the production process and very thin layers of paper, cardboard always has some imperfections, so it is important to be able to include this aspect in determining its mechanical parameters. The approach proposed in this paper extends the discussion carried out in our previous work [45], where the influence of imperfections on the stiffness of three-ply cardboard was investigated. The mechanics of both types differ significantly and, therefore, in this work other techniques, adequate to solve this problem, were used. The extension consists of taking into account the influence of the geometrical imperfections of individual layers on the calculated crosssection's effective stiffness. These imperfections were taken into account in the model by appropriate modeling of geometrical imperfections, which were determined in a priori buckling analysis and were appropriately included using the innovative technique based on numerical homogenization of double-walled corrugated board. Thanks to this approach, it is possible not only to quickly obtain homogenized stiffnesses data for complex crosssections, but also to easily take into account the effects of geometrical imperfections.

Corrugated Cardboard-Material and Geometry
Corrugated cardboard is a composite material made of several flat and corrugated layers of paper, called "liners" and "flutings", respectively. Due to the fibrous structure of the paper, the cardboard is characterized by a strong orthotropy. For this reason, the mechanical properties of the entire composite depend on the direction of fiber arrangement in individual layers. Two main directions can be distinguished, namely the machine direction (MD) and cross direction (CD) (see Figure 1). They result from the production process, in which layers of paper are rolled from multi-tone bales. The cardboard in the MD is stiffer but less ductile than in the CD (see Figure 2). The lower strength of the layers in the cross direction is partly compensated by the corrugation of the layers, which "adds" the material in this direction.  Linear elastic orthotropic material can be represented by the stress-strain relationship, as follows: where and are the strain and stress vector components, and are the Young's moduli in the machine and cross directions, respectively, and are the Poisson's coefficients, is the in-plane shear modulus, and are the transverse shear moduli. The compliance/stiffness matrix is symmetrical, so the relationship between Poisson's ratios can be written as follows: ( In this work, the paper layers were modeled using the linear elastic classical orthotropy described above. The material data were taken from the literature [28] and are presented in Table 1. The thickness of both flat and corrugated layers was assumed to be 0.30 mm.   Linear elastic orthotropic material can be represented by the stress-strain relationship, as follows: where and are the strain and stress vector components, and are the Young's moduli in the machine and cross directions, respectively, and are the Poisson's coefficients, is the in-plane shear modulus, and are the transverse shear moduli. The compliance/stiffness matrix is symmetrical, so the relationship between Poisson's ratios can be written as follows: ( In this work, the paper layers were modeled using the linear elastic classical orthotropy described above. The material data were taken from the literature [28] and are presented in Table 1. The thickness of both flat and corrugated layers was assumed to be 0.30 mm.  Linear elastic orthotropic material can be represented by the stress-strain relationship, as follows: where ε ij and σ ij are the strain and stress vector components, E 1 and E 2 are the Young's moduli in the machine and cross directions, respectively, ν 12 and ν 21 are the Poisson's coefficients, G 12 is the in-plane shear modulus, G 13 and G 23 are the transverse shear moduli. The compliance/stiffness matrix is symmetrical, so the relationship between Poisson's ratios can be written as follows: In this work, the paper layers were modeled using the linear elastic classical orthotropy described above. The material data were taken from the literature [28] and are presented in Table 1. The thickness of both flat and corrugated layers was assumed to be 0.30 mm.

Imperfections-Numerical Study
The main goal of the work is the numerical analysis of many cases of a corrugated board samples with imperfections and their impact on the stiffness values. In Figure 3, two imperfection shapes are presented; these were considered in this study. The first shape is the most common and corresponds to compression in the machine direction (see Figure 3a). It is formed as a result of thermal and moisture processes during the production of cardboard. The second imperfection mode (see Figure 3b) corresponds to the compression of the sample in the cross direction. This buckling shape was selected in [45] as the most representative for single-walled corrugated board. In both cases, it was assumed that only liners are involved in buckling analysis, and flutes (waves) only support the liners in the right position.

Imperfections-Numerical Study
The main goal of the work is the numerical analysis of many cases of a corrugated board samples with imperfections and their impact on the stiffness values. In Figure 3, two imperfection shapes are presented; these were considered in this study. The first shape is the most common and corresponds to compression in the machine direction (see Figure 3a). It is formed as a result of thermal and moisture processes during the production of cardboard. The second imperfection mode (see Figure 3b) corresponds to the compression of the sample in the cross direction. This buckling shape was selected in [45] as the most representative for single-walled corrugated board. In both cases, it was assumed that only liners are involved in buckling analysis, and flutes (waves) only support the liners in the right position. All calculations were made for double-walled cardboard composed of high B flute and low E flute (EB cardboard). In Table 2, the geometry of the corrugated layers is presented. The variants also include four imperfection levels, namely 0, 1, 2, and 3% for both imperfection shapes. The amount of the imperfection was taken as a certain part of the buckling length of the liners. Three different buckling lengths can be specified for the analyzed double-walled corrugated board (see Figure 4). The lengths and correspond to the periods of the B and E waves, respectively. The buckling length ′ can take different values depending on the relative position of the crests of E and B waves. The maximum possible value of is equal to the period of the lower wave, but most often every second segment is divided into two parts by the top of the B wave. For this reason, the buckling length of the middle liner was assumed to be 2/3 of the E wave period.  All calculations were made for double-walled cardboard composed of high B flute and low E flute (EB cardboard). In Table 2, the geometry of the corrugated layers is presented. The variants also include four imperfection levels, namely 0, 1, 2, and 3% for both imperfection shapes. The amount of the imperfection was taken as a certain part of the buckling length of the liners. Three different buckling lengths can be specified for the analyzed double-walled corrugated board (see Figure 4). The lengths L 1 and L 3 correspond to the periods of the B and E waves, respectively. The buckling length L 2 can take different values depending on the relative position of the crests of E and B waves. The maximum possible value of L 2 is equal to the period of the lower wave, but most often every second segment is divided into two parts by the top of the B wave. For this reason, the buckling length of the middle liner was assumed to be 2/3 of the E wave period.  Additionally, due to the irregular arrangement of the B and E flutes, the effect of the flute offset on the reduction in cardboard stiffness was verified. For this purpose, for each set of shape and level of imperfection, 10 cases of shifting the lower wave relative to the higher one, ranging from 0% to 90% in increments of 10%, were considered. In Figure 5, all simulated cross-sections are shown. Taking into account all the variants described above, each case can be marked with the symbol XX-Y-ZZ. The first two letters XX can be "MD" or "CD", which indicates the imperfection shape corresponding to compression in the MD and CD, respectively. The third sign Y represents the imperfection value and can be 0, 1, 2, or 3, and ZZ is the offset of the lower wave relative to the higher flute and can be 00, 10, 20, 30, 40, 50, 60, 70, 80, or 90.

Homogenization Procedure
The effect of imperfections on the effective stiffnesses of corrugated board was investigated with the numerical homogenization method proposed by Biancolini in 2005 and extended by Garbowski and Gajewski in 2021. It consists of ensuring the equivalence of the strain energy between the full 3D model of a representative volume element (RVE) and a simplified model of a flat 2D plate. The RVE is a small repeatable fragment of the full 3D model, which in this case is one period of the higher flute of the corrugated board. In the paper, only the main assumptions of the homogenization method used were presented. Detailed information and full derivation can be found in [29].
The homogenization procedure is based on the finite element method. Displacements from linear analysis can be represented by the following basic formula: Additionally, due to the irregular arrangement of the B and E flutes, the effect of the flute offset on the reduction in cardboard stiffness was verified. For this purpose, for each set of shape and level of imperfection, 10 cases of shifting the lower wave relative to the higher one, ranging from 0% to 90% in increments of 10%, were considered. In Figure 5, all simulated cross-sections are shown.  Additionally, due to the irregular arrangement of the B and E flutes, the effect of the flute offset on the reduction in cardboard stiffness was verified. For this purpose, for each set of shape and level of imperfection, 10 cases of shifting the lower wave relative to the higher one, ranging from 0% to 90% in increments of 10%, were considered. In Figure 5, all simulated cross-sections are shown. Taking into account all the variants described above, each case can be marked with the symbol XX-Y-ZZ. The first two letters XX can be "MD" or "CD", which indicates the imperfection shape corresponding to compression in the MD and CD, respectively. The third sign Y represents the imperfection value and can be 0, 1, 2, or 3, and ZZ is the offset of the lower wave relative to the higher flute and can be 00, 10, 20, 30, 40, 50, 60, 70, 80, or 90.

Homogenization Procedure
The effect of imperfections on the effective stiffnesses of corrugated board was investigated with the numerical homogenization method proposed by Biancolini in 2005 and extended by Garbowski and Gajewski in 2021. It consists of ensuring the equivalence of the strain energy between the full 3D model of a representative volume element (RVE) and a simplified model of a flat 2D plate. The RVE is a small repeatable fragment of the full 3D model, which in this case is one period of the higher flute of the corrugated board. In the paper, only the main assumptions of the homogenization method used were presented. Detailed information and full derivation can be found in [29].
The homogenization procedure is based on the finite element method. Displacements from linear analysis can be represented by the following basic formula: Taking into account all the variants described above, each case can be marked with the symbol XX-Y-ZZ. The first two letters XX can be "MD" or "CD", which indicates the imperfection shape corresponding to compression in the MD and CD, respectively. The third sign Y represents the imperfection value and can be 0, 1, 2, or 3, and ZZ is the offset of the lower wave relative to the higher flute and can be 00, 10, 20, 30, 40, 50, 60, 70, 80, or 90.

Homogenization Procedure
The effect of imperfections on the effective stiffnesses of corrugated board was investigated with the numerical homogenization method proposed by Biancolini in 2005 and extended by Garbowski and Gajewski in 2021. It consists of ensuring the equivalence of the strain energy between the full 3D model of a representative volume element (RVE) and a simplified model of a flat 2D plate. The RVE is a small repeatable fragment of the full 3D model, which in this case is one period of the higher flute of the corrugated board. In the paper, only the main assumptions of the homogenization method used were presented. Detailed information and full derivation can be found in [29].
The homogenization procedure is based on the finite element method. Displacements from linear analysis can be represented by the following basic formula: where K e is the stiffness matrix of the RVE after static condensation, u e is a displacement vector of the external nodes, and F e is the vector of external forces applied to the considered nodes. The finite element mesh and external nodes are presented in Figure 6. where is the stiffness matrix of the RVE after static condensation, is a displacement vector of the external nodes, and is the vector of external forces applied to the considered nodes. The finite element mesh and external nodes are presented in Figure 6. The static condensation procedure is based on the removal of unnecessary unknown degrees of freedom (DOF) and leaving the considered degrees of freedom (called principal DOFs or primary unknowns). In the case of cardboard, the primary unknowns are external RVE nodes. This results in bringing the stiffness matrix only to the external nodes of the model. The FE stiffness matrix after static condensation can be calculated from the following equation: where the components of the stiffness matrix are written for internal (subscript ) and external (subscript ) nodes, as follows: By the static condensation, the equation of the total elastic strain energy is simplified to the multiplication of the nodal displacements and the external forces acting on these nodes, as follows: Maintaining the appropriate properties of the simplified model is ensured by the balance of the total energy between the 3D model of corrugated board and the 2D plate model. The appropriate definition of displacements at the external edges of the RVE and allowing for bending and membrane behaviors works to ensure the energy balance. The relationship between generalized displacements and strains can be represented by the transformation matrix , as follows: and considering a single node ( = , = , = ) the above equation can be written in matrix form, as follows: The static condensation procedure is based on the removal of unnecessary unknown degrees of freedom (DOF) and leaving the considered degrees of freedom (called principal DOFs or primary unknowns). In the case of cardboard, the primary unknowns are external RVE nodes. This results in bringing the stiffness matrix only to the external nodes of the model. The FE stiffness matrix after static condensation can be calculated from the following equation: where the components of the stiffness matrix are written for internal (subscript i) and external (subscript e) nodes, as follows: By the static condensation, the equation of the total elastic strain energy is simplified to the multiplication of the nodal displacements and the external forces acting on these nodes, as follows: Maintaining the appropriate properties of the simplified model is ensured by the balance of the total energy between the 3D model of corrugated board and the 2D plate model. The appropriate definition of displacements at the external edges of the RVE and allowing for bending and membrane behaviors works to ensure the energy balance. The relationship between generalized displacements and strains can be represented by the transformation matrix H j , as follows: and considering a single node (x j = x, y j = y, z j = z) the above equation can be written in matrix form, as follows: Substituting the basic FEM equation and the relationship between generalized displacements and strains into Equation (6) results in the following: and considering a finite element subjected to standard load conditions, such as bending, stretching, shear, elastic internal energy, this can be represented by the following equation: where the matrix H k is the RVE stiffness matrix, which can be determined from the following: The stiffness matrix H k consists of several submatrices corresponding to the basic load conditions, as follows: where A represents the tensile and in-plane shear stiffnesses, D represents the bending and torsional stiffnesses, B is the coupling subarray of the tension and bending stiffnesses, and R is the transversal shear stiffnesses.
The stiffness values in subarray A do not depend on the position of a neutral axis. This is not true in the case of subarray B, in which the values change depending on the position of the neutral axis. In most symmetrical cases, the stiffness matrix B is a zero matrix, unlike for asymmetric cross-sections (e.g., multi-walled corrugated cardboard), for which some components of the matrix B are not zero, which affects the stiffness in the subarray D. This effect can be suppressed in two ways. The first is to choose the position of the neutral axis so as to minimize the values in the subarray B. An alternative solution is to calculate the uncoupled matrix D from the following formula: where D contains the original (coupled) bending and torsional stiffnesses for the non-zero matrix B. Numerical computations consisted in applying the selected buckling shape and imperfection value, and then calculating the stiffness matrix of the RVE. In all analyses, the four-node quadrilaterals shell elements with full integration (S4 elements) were used [47]. All models were built in FE commercial software, from where the stiffnesses in all nodes were generated. Next, the effective stiffness matrix was determined using the homogenization method described above. An approximate global size equal to 0.1625 mm was assumed. Due to the different variants of offset, the number of elements changed. For the MD-0-00 case, the model consisted of 9922 nodes, 9760 elements, and 59,532 degrees of freedom. The sensitivity study of the mesh size was already studied in our previous works; therefore, here this part of the research is omitted. It is worth noting that in the presented homogenization method, the use of the finite element method is limited only to building a global stiffness matrix, and no formal loads or boundary conditions are defined. They are directly implemented through the H e matrix in Equation (11).

Results
The first step was to properly create the corrugated cardboard RVE model in FE software. In the cases of samples with imperfections, the geometry of the liners was changed and replaced with the correct buckling shape (see Figure 3). Then, the material parameters contained in Table 1 were assigned to both liners and waves. Following the homogenization procedure described in Section 2.3, stiffness matrices were determined for all analyzed cases. In Table 3, an example of the stiffness matrix H k is presented for a model without imperfections and with 0% offset. No buckling mode is assigned to this case, because the imperfection value is equal to zero; therefore, it was marked with the XX-0-00 symbol.  As expected, non-zero components of the B matrix appeared in the stiffness matrix H k , which results from the asymmetry of the double-walled corrugated board sample. Tables 4-10 show the main diagonal components of the stiffness matrix for all analyzed cases. The components ( * ) 12 are omitted, but this does not create an error in the analysis because these elements are related to the components ( * ) 11 and ( * ) 12 in each stiffness matrix. The following tables do not list the elements of matrix B, but its effect on matrix D has been taken into account by applying Equation (13). The components of the R matrix were also omitted. The data contained in Tables 4-10 are also presented in the form of graphs in Figure 7.   (N·mm)  2153  2152  2146  2134  2124  2127  2140  2149  2153  2154  Table 5. The selected stiffnesses with 1% imperfections of compression in the MD shape and with different offsets.  Table 10. The selected stiffnesses with 3% imperfections of compression in the CD shape and with different offsets. For each shape and level of imperfection, the average stiffness values from all offset cases were calculated. Then, the obtained values for models with imperfections were compared to the average stiffnesses of perfectly shaped cardboard. In Figure 8, the tensile and bending stiffness reduction due to geometrical imperfections is shown.  For each shape and level of imperfection, the average stiffness values from all offset cases were calculated. Then, the obtained values for models with imperfections were compared to the average stiffnesses of perfectly shaped cardboard. In Figure 8, the tensile and bending stiffness reduction due to geometrical imperfections is shown. For each shape and level of imperfection, the average stiffness values from all offset cases were calculated. Then, the obtained values for models with imperfections were compared to the average stiffnesses of perfectly shaped cardboard. In Figure 8, the tensile and bending stiffness reduction due to geometrical imperfections is shown. As can be seen in Tables 4-10 and Figure 7, the stiffnesses vary depending on the offset. Considering each set of shape and level of imperfection, the difference between the maximum and minimum stiffness values was calculated. Then, the obtained difference As can be seen in Tables 4-10 and Figure 7, the stiffnesses vary depending on the offset. Considering each set of shape and level of imperfection, the difference between the maximum and minimum stiffness values was calculated. Then, the obtained difference was related to the average stiffness obtained from all 10 offset cases, obtaining a percentage relative spread. In Table 11, the relative spread for all analyzed sets of imperfection is presented.

Discussion
All crucial results have been shown in Tables 3-11 and Figures 7 and 8. In Table 3, an example of a stiffness matrix determined for the XX-0-00 model is presented, in which the submatrices A, B, D, and R are distinguished. The B matrix contains non-zero components because the XX-0-00 case cross-section, as in all other models, is asymmetric. For this reason, Equation (13) was used in all corrugated board cases to calculate the uncoupled matrix D, whose values are gathered in Tables 4-10. Comparing the stiffness values D 11 , D 22 , and D 33 from Tables 3 and 4, it can be seen that the values in Table 4 are smaller, which results from the uncoupling of the D matrix.
Considering Tables 4-10 and Figure 7, where the stiffness results for all models are presented, it can be seen that as the level of imperfection increases, the cardboard stiffness decreases, as expected. The line diagrams also show that depending on the applied wave offset, the stiffnesses change their value. This is shown in more detail in Table 11, which shows the relative spread of the extreme stiffnesses. The greatest fluctuations occur for the A 11 stiffnesses, are slightly smaller for all bending stiffnesses, and are negligibly small for the A 22 and A 33 stiffnesses. In most cases, it is clear that the relative spread increases as the amount of imperfection increases. In addition, the data in this table show that the stiffnesses of the models with the MD imperfection shape are more sensitive to offset changes than the cases with the CD buckling mode.
Based on Tables 4-10 and Figures 7 and 8, it can be concluded that geometric imperfections have a huge influence on the bending stiffness D 11 . Large decreases were also noted for the stiffnesses A 11 , D 22 , and D 33 . The effect of imperfections on the stiffness A 22 is small, and the effect is negligible on the A 33 stiffness. These results are very similar to those obtained from the analysis of the effect of imperfections on the stiffness of single-walled corrugated board conducted by Mrówczyński et al. [45], where the largest decreases were also obtained for D 11 , and the smallest were obtained for A 33 . The implementation of the MD buckling mode results in a greater average reduction in stiffnesses A 11 , D 11 , and D 33 compared to the CD imperfection shape, for which a higher average reduction occur for the stiffnesses of A 22 , A 33 , and D 22 . These differences are not that significant and amount to a maximum of 2.75 % (for the stiffness D 11 between the MD-3-ZZ and CD-3-ZZ models).
To simplify computations and save computational time, it is desirable to select one case that should be computed as representative instead of carrying out an entire numerical study. To determine this, several observations from the conducted analyses should be considered. The imperfection shape of the compression in the CD has smaller relative spread values than the buckling mode of compression in the MD, so the change in the offset does not affect the stiffnesses so significantly. In addition, considering the MD shape does not significantly reduce the stiffnesses A 22 and A 33 , which seems to be an unfavorable effect. Based on the calculations, it can be concluded that a reasonable solution is to choose the CD imperfection shape as a representative model. Due to the quite small relative spread, the choice of a representative offset value is not so important, so its value can be chosen arbitrarily. The conclusions drawn from the presented numerical study are consistent with the work of Mrówczyński et al. [45], in which the authors showed that the CD imperfection shape is also the most representative buckling mode for single-walled corrugated cardboard.
All the observations discussed above prove that there is a direct influence of imperfections on the stiffnesses A, B, and D. However, the basic question arises whether it is possible to verify the presented theoretical approach with experimental data. Although testing corrugated board is not very laborious or particularly expensive, nevertheless validating the theoretical analysis presented here is not an easy task, because, firstly, it is practically impossible to produce corrugated board with a given level of imperfection and, secondly, it is also very laborious to check the actual level of imperfection in the cardboard produced. Nevertheless, below, we present a comparison (validation) of the results obtained using the method presented here and the method presented in paper [44] in comparison with the results of laboratory tests conducted by Czechowski et al. and presented in paper [48].
All experimental details are described in [48], while analytical assumptions, in which imperfections are also included, are given in [44]. Here, only the results of the experimental campaign (namely, results of 4-point bending tests of various asymmetric double-walled boards) are shown together with the results obtained using numerical and analytical models. All these results are presented in Table 12, together with the results obtained using the methods presented here with the same amount of imperfection of flat layers as assumed in [44]. It is worth noting that only one of all stiffnesses is checked, namely the bending stiffness in MD. The obtained results do not differ by more than 5% from the results obtained using the analytical model, while the average absolute error (compared to the experimental data) was 7.5%.

Conclusions
The article presents an approach of numerical homogenization of double-walled corrugated board, taking into account the initial imperfections of its layers. The creation of many cardboard models in various configurations made it possible to study various aspects of the considered problem. In the study, corrugated cardboard was homogenized using the numerical method based on the strain energy equivalence between a representative 3D model and a simplified plate. The analysis of five-layered cardboard with imperfections is not such an easy task. The cross-sectional shape of the cardboard depends on the storage and load conditions. For this reason, it is important to indicate a representative case that adequately reflects the properties of the material.
Based on the conducted analyses, it can be seen that the bending stiffnesses and the tensile stiffness along the wave are the most sensitive to the imperfect shape of the cardboard. Calculations were made for a wide range of cases, which allowed us to choose a representative imperfection model. The conclusions drawn from the numerical analyses are consistent with the research that was carried out in the previous work for single-walled corrugated board.

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