Numerical Homogenization of Single-Walled Corrugated Board with Imperfections

Abstract

Numerical homogenization is an excellent tool for the quick simplification of complex structures with a model that is much simpler and, at the same time, accurately reflects the mechanical behavior of the original model. Corrugated cardboard modeling, with all geometrical nuances preserved, is a complicated and time-consuming process. The transfer of a full 3D model of corrugated board composed of two flat layers and a corrugated middle layer to one layer only, with substitute elastic parameters, greatly simplifies this process. Because the individual layers of corrugated cardboard are made of paper with a grammage in the range of 80–200 g/m², i.e., very thin plates, they are slightly buckled even in the initial configuration. These imperfections affect the equivalent parameters that are obtained in the homogenization process. This paper presents an approach of taking into account these imperfections when creating a simplified model. The numerical homogenization method based on the equivalence of elastic energy between a representative volumetric element (i.e., a part of a full 3D model) and an equivalent plate were applied. Different shapes of imperfections were analyzed in order to account for the buckling modes, notably for a specific unit deformation and curvature. Finally, one form of initial imperfections was proposed, which most accurately reflects the decrease in all plate stiffnesses.

1. Introduction

In recent years, there has been observed a constant growth in the popularity of corrugated board, which is widely used in the packaging industry. This was due to the great demand for packaging materials in the food, pharmaceutical and cosmetic industries. The COVID-19 global pandemic has contributed to the rapid growth of e-commerce, driven by a significant increase in individual shipments and home orders.

Nowadays, a strong argument in favor of the use of corrugated board is its eco-friendliness. Cardboard packaging is very easy to recycle or dispose of. Packaging can be made with any color prints, which is an attractive option for e-commerce companies as they can freely shape the appearance of the packaging.

The most important aspect of packaging is the appropriate protection of the goods during storage or transport. In order to correctly determine the load capacity of the packaging, a number of factors should be taken into account, such as the influence of humidity, the arrangement of perforations and holes, the type of the load or the location of the flaps. Manufacturers of cardboard packaging are searching for solutions that are economical and effective, which is why research on the durability and optimal design of corrugated cardboard products is constantly carried out [1,2]. Various tests are performed to assess the load capacity of boxes, including compressive, tensile and bending strength. The most common tests are the box compression test (BCT) and the edge crush test (ECT).

The load capacity of the packaging has been the subject of research conducted by many scientific groups for many years. As early as 1952, a formula for the strength of the packaging which takes into account the box perimeter, the overall ring crush strength and the box and corrugated layer constants was proposed [3]. Maltenfort presented the relationship between the critical force in BCT and the dimensions of the packaging based on the Concora liner test and empirical constants [4]. In 1963, McKee et al. proposed an approach that became one of the most popular methods of estimating the load capacity of packaging with a rectangular base [5]. The formula is based on the ECT value, box perimeter, bending stiffness of the walls and correction factors. The use of empirically determined correction factors reduces the versatility of the method because it is necessary to compute new values of these factors for each board quality and each packaging model.

In 1968, a formula for estimating the strength of the regular package was presented. As the McKee formula, it was based on the package perimeter, the ECT and plate bending stiffness, and, additionally, also took into account the buckling factor, stacking time and load ratio [6]. In 1985, Allerby et al. proposed other constants and exponents of the McKee formula [7], and in 1987, Schrampfer et al. modified the McKee formula for more cardboard cutting methods and equipment [8]. Kawanishi presented a formula for the load capacity of the packaging, which depends on several important factors, such as the grammage of the cardboard layers, cardboard thickness, take-up factor, average corrugation count, perimeter and type of the box, moisture of the walls and printed ratio [9]. In 1993, Batelka and Smith proposed an extension of the McKee formula to several other types of boxes, taking into account the width and depth of the package [10].

The finite element method (FEM) is often used for numerical analysis of packaging strength. Urbanik and Frank used FEM to compare the strength of the packaging with a formula based on, e.g., the Poisson’s ratio [11]. Nordstrand and Carlsson performed FEM simulations to compare the obtained effective cardboard transverse shear moduli with analytical predictions and experimental results [12,13]. Urbanik and Saliklis presented observations of buckling phenomena of corrugated cardboard boxes modeled with FEM [14]. In 2011, an overview of the analytical and numerical load capacity estimation of corrugated board packaging was presented by Sohrabpour and Hellström [15].

The layered structure of the cardboard and the anisotropy of the paper make numerical modeling difficult due to the necessity of knowing the material parameters of all corrugated board layers. The solution to this problem is a process called homogenization. Homogenization allows to reduce the complicated cross-section of the corrugated board to one layer, which significantly simplifies the numerical model and shortens the computation time. Using this method, the equivalent stiffnesses and the effective thickness of the model should be determined to ensure the same behavior before and after homogenization.

The analytical homogenization of corrugated board is based on the equations of the classical theory of material strength and the classical theory of laminates [16]. In numerical homogenization, the basic equations of the finite element method are also used. The authors in this work utilize the method proposed by Biancolini [17] and then extended by Garbowski and Gajewski [18], which is based on an energy balance between a representative volume element and an equivalent plate. Hohe presented a representative element of heterogeneous and homogenized elements, which is based on the strain energy [19]. By means of periodic homogenization, Buannic et al. obtained equivalent membrane and bending characteristics of the plates [20]. In 2010, Abbès and Guo determined the torsional stiffness of an orthotropic sandwich panel by splitting the panel into two beams [21]. In 2018, Ramírez-Torres et al. presented a multiscale asymptotic homogenization for a composite material with a layered hierarchical structure [22,23]. Gallo et al. proposed an approximation of the load capacity and deformation of a package while applying empirical observations [24]. The use of homogenization to analyze a thin-walled seating made of triple-wall corrugated cardboard is presented by Suarez et al. [25]. Garbowski et al. investigated the effect of creasing or perforation on the stiffness of corrugated board [26]. Mrówczyński et al. performed a sensitivity analysis in optimal design of single- and double-walled corrugated board packaging [1,2]. Nguyen-Minh et al. verified the accuracy and reliability of the homogenization methods in the static analysis of trapezoidally and sinusoidally corrugated panels [27].

The load capacity of the corrugated board depends on several factors. One of them is the imperfections of the cardboard layers resulting from the production process. In 1995, Nordstrand carried out parametric tests to investigate the sensitivity of the load capacity with the size of initial imperfections [28]. The influence of geometrical imperfections on the mechanical properties was presented by Lu et al. [29]. Nordstrand modified the nonlinear buckling analysis of Rhodes and Harvey orthotropic plates to account for initial cardboard imperfections [30]. Garbowski and Knitter-Piątkowska presented an analytical method of determining the bending stiffness of double-walled corrugated cardboard with imperfections [31].

In the present article, an approach to account for the geometric imperfections of corrugated cardboard in numerical homogenization of single-walled corrugated board is presented. The approach presented in this paper extends and generalizes the homogenization methods already presented in our previous works. This generalization concerns geometric imperfections built into the material. They are always present due to the very thin layers of paper used in the production of corrugated board and, therefore, cannot be ignored. The presented here method allows to easily take into account these geometric imperfections and thus correctly calculate the sought effective stiffnesses.

2. Materials and Methods

2.1. Material Parameters and Corrugated Cardboard Geometry

Corrugated cardboard is a fibrous material, and therefore it is highly orthotropic. The direction of the fibers in the individual cardboard layers determines the mechanical properties of the board components. For wood-based materials, two main directions can be distinguished: the machine direction (MD) and the cross direction (CD). The strength of paper and paperboard is more than twice as high in the MD than in the CD because the fibers are arranged along the MD in the production process (see Figure 1). In the cross direction, the material is more ductile, while in the machine direction, it is more resistant to tearing and crushing.

The linear behavior of an elastic orthotropic material can be represented by the strain–stress relationship:

$$\left[\begin{array}{c}{\epsilon }_{11}\\ {\epsilon }_{22}\\ 2{\epsilon }_{12}\\ 2{\epsilon }_{13}\\ 2{\epsilon }_{23}\end{array}\right]=\left[\begin{array}{ccccc}1/E_1 & -{\nu }_{21}/E_2& 0& 0& 0\\ -{\nu }_{12}/E_1& 1/E_2& 0& 0& 0\\ 0& 0& 1/G_{12}& 0& 0\\ 0& 0& 0& 1/G_{13}& 0\\ 0& 0& 0& 0& 1/G_{23}\end{array}\right]\left[\begin{array}{c}{\sigma }_{11}\\ {\sigma }_{22}\\ {\sigma }_{12}\\ {\sigma }_{13}\\ {\sigma }_{23}\end{array}\right],$$

(1)

where $E_1$ and $E_2$ are the Young’s moduli in the MD and the CD, respectively; ${\nu }_{12}$ and ${\nu }_{21}$ are the Poisson’s coefficients; $G_{12}$ is the Kirchhoff’s modulus; $G_{13}$ and $G_{23}$ are the transverse shear moduli. The symmetry of the material stiffness/compliance matrix allows us to determine the relationship between the Poisson’s coefficients and the Young’s moduli:

$$\frac{{\nu }_{12}}{E_1}=\frac{{\nu }_{21}}{E_2}.$$

(2)

The material orientation is always the same in all layers of corrugated cardboard. The machine direction coincides with the cardboard wave direction, which results from the technological conditions of the production process (see Figure 2). The lower strength in the cross direction is compensated by a take-up factor of the corrugated layer.

All corrugated board layers were modeled while using classical linear elastic orthotropy, see Equation (1). The material data were taken from the literature [17] and are presented in

Table 1. Material data used in orthotropic constitutive model.

Layers	${\mathit{E}}_1$	${\mathit{E}}_2$	${\mathit{\nu }}_{12}$	${\mathit{G}}_{12}$	${\mathit{G}}_{13}$	${\mathit{G}}_{23}$
	(MPa)	(MPa)	(-)	(MPa)	(MPa)	(MPa)
liners	3326	1694	0.34	859	429.5	429.5
fluting	2614	1532	0.32	724	362	362

. The thickness of both flat layers (liners) and the corrugated layer (fluting) was assumed to be 0.30 mm.

2.2. Homogenization Technique

The numerical homogenization method proposed by Biancolini [17], and then extended by Garbowski and Gajewski [18] was utilized in the present study. This method is based on the provision of strain energy equivalence between the full representative volume element (RVE) of the corrugated board and the simplified shell model. The RVE is a small, periodically repeating segment of the full corrugated board structure model that is then simplified to a single layer of flat shell elements in such a way as to ensure the same behavior of the full and simplified model. Only the most important assumptions of the described method are presented below. The entire theoretical derivation of the constitutive model can be found in [18]. However, it is important to note that for the evaluation of homogenized properties for macroscopic analysis, the concept of statistical RVE (SRVE) can be utilized [32,33,34].

The nodal displacements can be determined while involving the finite element method (FEM). The FEM equation for linear analysis is shown below:

$$K_e u_e=F_e,$$

(3)

where $K_e$ is a global stiffness matrix statically condensed to the external nodes of the RVE, $u_e$ is a displacement vector of the external nodes and $F_e$ is the vector of nodal forces applied to the external nodes. In Figure 3, the finite element mesh and nodes of the RVE are shown.

The global stiffness matrix was obtained by applying the static condensation process. It consists in removing selected unknown degrees of freedom (DOF) and determining the stiffness matrix for a smaller number of DOFs (called primary unknowns). In the discussed model, the internal nodes were removed, and for the external nodes, the stiffness matrix was determined from the following formula:

$$K_e=K_{ee}-K_{ei} {K_{ii}}^{-1} K_{ie},$$

(4)

where all the subarrays $K_{**}$ contain the stiffnesses for external (subscript $e$) and internal (subscript $i$) nodes, which in matrix form relate nodal displacements and nodal forces:

$$\left[\begin{array}{cc}{K}_{ee}& K_{ei}\\ K_{ie}& K_{ii}\end{array}\right]\left[\begin{array}{c}{u}_e\\ u_i\end{array}\right]=\left[\begin{array}{c}{F}_e\\ 0\end{array}\right].$$

(5)

As a result of static condensation, the total elastic strain energy was reduced to the work of external forces on the corresponding displacements:

$$E=\frac{1}{2} u_e^T F_e.$$

(6)

The balance of the total strain energy between the full 3D RVE model and the simplified shell model was ensured by appropriate determination of nodal displacements and by enabling the membrane and bending behavior. The details can be found in [18]. The relationship between generalized displacements and generalized strains can be presented by the formula:

$$u_i=H_i {\epsilon }_i,$$

(7)

where the $H_i$ matrix can be determined for each node ($x_i=x$, $y_i=y$, $z_i=z$):

$${\left[\begin{array}{c}{u}_x\\ u_y\\ u_z\\ {\theta }_x\\ {\theta }_y\end{array}\right]}_i={\left[\begin{array}{cccccccc}x& 0& y/2& xz& 0& yz/2& z/2& 0\\ 0& y& x/2& 0& yz& xz/2& 0& z/2\\ 0& 0& 0& -x^2/2& -y^2/2& -xy/2& x/2& y/2\\ 0& 0& 0& 0& -y& -x/2& 0& 0\\ 0& 0& 0& x& 0& y/2& 0& 0\end{array}\right]}_i{\left[\begin{array}{c}{\epsilon }_x\\ {\epsilon }_y\\ {\gamma }_{xy}\\ {\kappa }_x\\ {\kappa }_y\\ {\kappa }_{xy}\\ {\gamma }_{xz}\\ {\gamma }_{yz}\end{array}\right]}_i.$$

(8)

Using the definition of elastic strain energy for a discrete model:

$$E=\frac{1}{2}{u}_e^T K u_e=\frac{1}{2}{\epsilon }_e^T H_e^T K H_e {\epsilon }_e$$

(9)

and analyzing a finite element model as subjected to basic types of loads (i.e., bending, tension and transverse shear), the elastic internal energy can be represented as:

$$E=\frac{1}{2}{ϵ}_e^T H_{\mathrm{k}} {\epsilon }_e\left\{area\right\}.$$

(10)

The stiffness matrix for the homogenized composite can be determined from:

$$H_k=\frac{H_e^T K H_e}{area}.$$

(11)

The $H_k$ matrix consists of the stiffness matrices:

$$H_k=\left[\begin{array}{ccc}{A}_{3\times 3}& B_{3\times 3}& 0\\ B_{3\times 3}& D_{3\times 3}& 0\\ 0& 0& R_{2\times 2}\end{array}\right],$$

(12)

where the $A$ subarray contains tensile and shear stiffnesses, the $B$ represents tension-bending coupling stiffnesses, the $D$ contains bending and torsional stiffnesses, and in the $R$ subarray the transverse shear stiffnesses are gathered.

In the case of symmetrical cross sections, the $B$ matrix is the zero matrix. For asymmetrical RVE, non-zero matrix elements appear in the $B$ matrix, which affects the $D$ matrix values. This problem can be solved by minimizing the $B$ matrix by choosing the position of the neutral axis, or, alternatively, the following formula can be used:

$$D=D^{\prime }-B{A}^{-1}B,$$

(13)

where the $D^{\prime }$ matrix contains coupled bending and torsional stiffnesses for the non-zero $B$ matrix.

2.3. Numerical Model with Imperfections

The main goal of this work was to implement geometric imperfections into the numerical model and to determine their influence on the stiffness reduction of corrugated board. Cardboard with a B flute (see

Table 2. Geometric parameters of waves.

Wave (Flute)	Wave Length (mm)	Height (mm)	Take-Up Factor (-)
B	6.5	2.46	1.32

) and material parameters contained in Table 1 were analyzed.

The application of imperfections in the numerical model of corrugated board consists in the appropriate definition of the buckling mode of the RVE. In the layers subjected to compression, the influence of imperfections is crucial. Therefore, due to the nature of the load, two cases should be considered: compression and bending. During compression, the entire cross section of the RVE is compressed, i.e., both the bottom and top liners. With regard to bending, one liner is compressed (assumed to be the bottom liner) and the other is tensioned. In the case of the liner, which is stretched, the imperfection does not affect the load capacity (in other words, it is not present); therefore, the tensile layers can be modeled as undeformed. During compression and bending, all or part of the middle cardboard layer is compressed. However, it was assumed that the influence of imperfections in this layer on stiffness is negligible due to the already existing corrugated shape of fluting.

The load direction of the corrugated board also affects the mode of buckling. During compression or bending in the machine direction, imperfections in the compressed layers are in the form of a half-sine wave across the entire depth of the RVE, where the support is the connection between liners and fluting (see Figure 4a,b). On the other hand, if the cardboard is loaded in the cross direction, then a half-wave appears in both directions, which causes dome-shaped imperfections. In this case, the support was provided by the external edges of the liners in the cross direction and the connection of the liners with fluting (see Figure 4c,d). Corrugated board may also contain imperfections that reflect in-plane shear and twisting of the RVE. For this reason, it was necessary to simulate these two load cases (see Figure 4e,f).

Based on the above description, six modes of imperfection were analyzed, and thus six numerical models, which are presented in Figure 4. In each case, the amount of imperfections, i.e., the maximum deviation from the ideal shape (red color), was 5% of the cardboard height, which was 0.123 mm. For finite element method computation, the four-node quadrilaterals shell elements with reduced integration (S4R elements) were used. In all seven numerical models (one model without imperfections and six models with imperfections), the same mesh size was assumed. The approximate global size of the element was 0.1625 mm, which gave 4920 nodes, 4800 elements and 29,520 degrees of freedom.

3. Results

The corrugated cardboard model without initial imperfections and six models with six buckling modes were analyzed.

Table 3. Constitutive stiffness matrix $H_k$ for the cardboard model without imperfections.

		A (N/mm) & B (N)			B (N) $\mathbf{\&\; D}(\mathbf{N}\cdot \mathbf{m}\mathbf{m})$			R $\left(\mathbf{N}/\mathbf{m}\mathbf{m}\right)$
		1	2	3	1	2	3	4	5
A (N/mm) & B (N)	1	2234	401	0	−18	−5	0
	2	401	1684	0	−5	−2	0
	3	0	0	687	0	0	0
B (N) & D $(\mathrm{N}\cdot \mathrm{mm})$	1	−18	−5	0	3291	574	0
	2	−5	−2	0	574	2047	0
	3	0	0	0	0	0	865
R $(\mathrm{N}/\mathrm{mm}$)	4							104	0
	5							0	95

shows an example of the $H_k$ matrix computed for the RVE without imperfections (see Equation (12)).

Table 4. Constitutive stiffness matrix $H_k$ for the cardboard model with imperfections in the case of compression in the MD.

		A (N/mm) & B (N)			B (N) & D$(\mathbf{N}\cdot \mathbf{m}\mathbf{m})$			R $\left(\mathbf{N}/\mathbf{m}\mathbf{m}\right)$
		1	2	3	1	2	3	4	5
A (N/mm) & B (N)	1	2144	386	0	140	22	0
	2	386	1683	0	22	43	0
	3	0	0	686	0	0	20
B (N) & D $(\mathrm{N}\cdot \mathrm{mm})$	1	140	22	0	2975	520	0
	2	22	43	0	520	1949	0
	3	0	0	20	0	0	818
R $(\mathrm{N}/\mathrm{mm}$)	4							105	0
	5							0	95

and

Table 5. Constitutive stiffness matrix $H_k$ for the cardboard model with imperfections in the case of bending in the MD.

		A (N/mm) & B (N)			B (N) & D$(\mathbf{N}\cdot \mathbf{m}\mathbf{m})$			R $\left(\mathbf{N}/\mathbf{m}\mathbf{m}\right)$
		1	2	3	1	2	3	4	5
A (N/mm) & B (N)	1	2156	388	0	154	25	0
	2	388	1682	0	24	43	0
	3	0	0	686	0	0	20
B (N) & D $(\mathrm{N}\cdot \mathrm{mm})$	1	154	24	0	2990	522	0
	2	25	43	0	522	1945	0
	3	0	0	20	0	0	817
R $(\mathrm{N}/\mathrm{mm}$)	4							105	0
	5							0	95

show the stiffness matrices for imperfection modes in the case of compression and bending in the machine direction.

Table 6. Constitutive stiffness matrix $H_k$ for the cardboard model with imperfections in the case of compression in the CD.

		A (N/mm) & B (N)			B (N) & D$(\mathbf{N}\cdot \mathbf{m}\mathbf{m})$			R $\left(\mathbf{N}/\mathbf{m}\mathbf{m}\right)$
		1	2	3	1	2	3	4	5
A (N/mm) & B (N)	1	2140	349	0	101	41	0
	2	349	1653	0	41	51	0
	3	0	0	684	0	0	15
B (N) & D $(\mathrm{N}\cdot \mathrm{mm})$	1	101	41	0	3018	478	0
	2	41	51	0	478	1933	0
	3	0	0	15	0	0	826
R $(\mathrm{N}/\mathrm{mm}$)	4							105	0
	5							0	95

and

Table 7. Constitutive stiffness matrix $H_k$ for the cardboard model with imperfections in the case of bending in the CD.

		A (N/mm) & B (N)			B (N) & D$(\mathbf{N}\cdot \mathbf{m}\mathbf{m})$			R $\left(\mathbf{N}/\mathbf{m}\mathbf{m}\right)$
		1	2	3	1	2	3	4	5
A (N/mm) & B (N)	1	2161	359	0	127	54	0
	2	359	1659	0	53	58	0
	3	0	0	685	0	0	16
B (N) & D $(\mathrm{N}\cdot \mathrm{mm})$	1	127	53	0	3048	493	0
	2	54	58	0	493	1940	0
	3	0	0	16	0	0	828
R $(\mathrm{N}/\mathrm{mm}$)	4							104	0
	5							0	95

present the RVE stiffness matrices with imperfections in the case of compression and bending in the cross direction. In

Table 8. Constitutive stiffness matrix $H_k$ for the cardboard model with imperfections in the case of in-plane shear.

		A (N/mm) & B (N)			B (N) & D$(\mathbf{N}\cdot \mathbf{m}\mathbf{m})$			R $\left(\mathbf{N}/\mathbf{m}\mathbf{m}\right)$
		1	2	3	1	2	3	4	5
A (N/mm) & B (N)	1	2233	401	1	28	7	−1
	2	401	1684	0	8	7	0
	3	1	0	686	−1	0	2
B (N) & D $(\mathrm{N}\cdot \mathrm{mm})$	1	28	8	−1	3289	574	1
	2	7	7	0	574	2049	0
	3	−1	0	2	1	0	865
R $(\mathrm{N}/\mathrm{mm}$)	4							105	0
	5							0	95

and

Table 9. Constitutive stiffness matrix $H_k$ for the cardboard model with imperfections in the case of twisting.

		A (N/mm) & B (N)			B (N) & D$(\mathbf{N}\cdot \mathbf{m}\mathbf{m})$			R $\left(\mathbf{N}/\mathbf{m}\mathbf{m}\right)$
		1	2	3	1	2	3	4	5
A (N/mm) & B (N)	1	2204	372	4	−35	−33	4
	2	372	1652	3	−32	−46	3
	3	4	3	685	4	3	−6
B (N) & D $(\mathrm{N}\cdot \mathrm{mm})$	1	−35	−32	4	3209	524	5
	2	−33	−46	3	524	1982	4
	3	4	3	−6	5	4	854
R $(\mathrm{N}/\mathrm{mm}$)	4							106	0
	5							0	95

, the constitutive stiffness matrices for imperfections in the case of in-plane shear and twisting are compiled.

It should be noted that in the $B$ subarrays non-zero components are formed, which is caused by the asymmetry of the cross-section of the corrugated cardboard.

Table 10. Selected stiffnesses for models without imperfections and with six imperfection modes.

	REF	M1	M2	M3	M4	M5	M6
$A_{11} \left(\mathrm{N}/\mathrm{mm}\right)$	2234	2144	2156	2140	2161	2233	2204
$A_{22} \left(\mathrm{N}/\mathrm{mm}\right)$	1684	1683	1682	1653	1659	1684	1652
$A_{33} \left(\mathrm{N}/\mathrm{mm}\right)$	687	686	686	684	685	686	685
$D_{11} \left(\mathrm{N}\cdot \mathrm{mm}\right)$	3290	2966	2979	3013	3040	3289	3208
$D_{22} \left(\mathrm{N}\cdot \mathrm{mm}\right)$	2047	1948	1944	1931	1937	2049	1980
$D_{33} \left(\mathrm{N}\cdot \mathrm{mm}\right)$	865	817	816	826	827	865	854

shows the stiffness values of the main diagonal matrices that are presented in Table 3, Table 4, Table 5, Table 6, Table 7, Table 8 and Table 9 for the model without imperfections (REF) and the six models with buckling modes (M1–M6). The imperfection modes were marked in the order presented in Figure 4, i.e., M1—compression in the MD, M2—bending in the MD, M3—compression in the CD, etc. Omitting the components ${(\ast )}_{12}$ does not introduce an error in the analysis of the results, because these components depend on the components ${(\ast )}_{11}$, ${(\ast )}_{22}$ and on the Poisson’s ratio. The $B$ subarray was also disregarded, but its influence was taken into account by determining the uncoupled $D$ matrix from Equation (13). The components of the $R$ matrix were also omitted because their values remained almost unchanged for each buckling mode. Figure 5, Figure 6, Figure 7, Figure 8, Figure 9 and Figure 10 show the reduction of the individual stiffnesses for the analyzed six imperfection modes compared to the reference value that was computed from the model without imperfection.

Table 11. Stiffness degradation for models with imperfections (%).

	M1	M2	M3	M4	M5	M6
$A_{11}$	−4.0	−3.5	−4.2	−3.3	−0.06	−1.4
$A_{22}$	−0.1	−0.1	−1.8	−1.5	−0.01	−1.9
$A_{33}$	−0.1	0.0	−0.5	−0.2	−0.03	−0.2
$D_{11}$	−9.9	−9.5	−8.4	−7.6	−0.05	−2.5
$D_{22}$	−4.8	−5.1	−5.7	−5.4	0.07	−3.3
$D_{33}$	−5.5	−5.6	−4.5	−4.4	0.00	−1.3

presents the differences between the stiffnesses for models with imperfections and the stiffnesses of the reference model (without imperfections).

4. Discussion

On the basis of the analyzed six buckling modes, reduced stiffness of the corrugated board was obtained due to the consideration of geometric imperfections (see Table 3, Table 4, Table 5, Table 6, Table 7, Table 8, Table 9 and Table 10). In

Table 12. Final stiffnesses of corrugated board with imperfections.

	Reference Stiffnesses	Final Stiffnesses	Stiffness Reduction (%)
$A_{11} \left(\mathrm{N}/\mathrm{mm}\right)$	2234	2144	−4.0
$A_{22} \left(\mathrm{N}/\mathrm{mm}\right)$	1684	1653	−1.8
$A_{33} \left(\mathrm{N}/\mathrm{mm}\right)$	687	686	−0.03
$D_{11} \left(\mathrm{N}\cdot \mathrm{mm}\right)$	3290	2979	−9.5
$D_{22} \left(\mathrm{N}\cdot \mathrm{mm}\right)$	2047	1937	−5.4
$D_{33} \left(\mathrm{N}\cdot \mathrm{mm}\right)$	865	854	−1.3

, the final stiffnesses of the corrugated cardboard with imperfections are shown. The individual stiffness matrix components were selected from the corresponding imperfection modes as follows: $A_{11}$—compression in the MD, $A_{22}$—compression in the CD, $A_{33}$—in-plane shear, $D_{11}$—bending in the MD, $D_{22}$—bending in the CD, $D_{33}$—twisting.

The data presented in Table 10 and Table 11 and Figure 5, Figure 6, Figure 7, Figure 8, Figure 9 and Figure 10 show that the stiffness $A_{33}$ did not change significantly for all analyzed buckling modes. The maximum decrease in the stiffness $A_{33}$ was 0.5%. The presented data, in particular the data summarized in Table 11, clearly show that taking into account the imperfection mode M5 causes negligible reduction of all stiffnesses. It can also be noticed that the greatest differences in stiffness (due to taking imperfections into account) occur for the stiffness $D_{11}$, and a slightly smaller reduction occurs for the stiffness $D_{22}$ (see Table 12). The obtained values show that the bending stiffnesses are the most sensitive to disturbances in the ideal shape of the corrugated cardboard.

In order to save computational time, it was possible to analyze only one buckling mode instead of all six. Therefore, it was necessary to choose the buckling mode that gives the most similar results to those determined from the six imperfection cases. To determine this mode, the final stiffnesses from Table 12 were compared with the stiffness values for each buckling mode from Table 10. The smallest total error was obtained for the buckling mode M3, i.e., the compression in the CD, and it amounted to 5.3%. Slightly worse results were obtained for the modes M4—6.5%, M1 and M2—7.1%, and M6—12.9%. The worst fit was obtained for the in-plane shear case (mode M5). For this imperfection mode, the total error was as high as 23.5%. Based on the above results, it can be concluded that the most reliable results are obtained by analyzing the imperfections in the case of compression in the CD (see Figure 4c). This means that by using the imperfection geometry from the M5 model, it was possible to obtain drops in the values of the main stiffnesses, which correspond to the drops of the same stiffnesses but while using models with dedicated imperfections.

5. Conclusions

Nowadays, great emphasis is placed on the most optimal use of natural resources, including the production of corrugated cardboard. For this purpose, many research groups are trying to understand the behavior of cardboard depending on various factors in order to correctly model it. Taking the imperfections into account in the process of determining the equivalent stiffnesses is crucial because the geometrical deformations of the individual paper layers in the corrugated board are always present and cannot be ignored. The article presents an approach to taking into account such geometric imperfections in the modeling of corrugated board. The homogenization method, consisting in simplifying the full 3D model of the corrugated board to a single layer with equivalent parameters, was used. The numerical homogenization process was based on the elastic strain energy equivalence between the representative volume element and the equivalent plate. Taking into account geometric imperfections in the homogenization process is not obvious. The material may have initial imperfections of various shapes or imperfections, which are shaped by the loading conditions. Therefore, finding a universal geometry that reflects possible imperfections is very important for the correct determination of the effective parameters of a homogenized plate with a periodic core. On the basis of the conducted analyses, several conclusions were formulated regarding the influence of imperfections on the stiffness of the corrugated board, which enable a better understanding of the behavior of the corrugated board and thus more conscious design of the packaging structure made of this material.