Non-Local Sensitivity Analysis and Numerical Homogenization in Optimal Design of Single-Wall Corrugated Board Packaging

The optimal selection of the composition of corrugated cardboard dedicated to specific packaging structures is not an easy task. The use of lighter boards saves material, but at the same time increases the risk of not meeting the guaranteed load capacity. Therefore, the answer to the question “in which layer the basis weight of the paper should be increased?” is not simple or obvious. The method proposed here makes it easy to understand which components and to what extent they affect the load-bearing capacity of packages of various dimensions. The use of numerical homogenization allows for a quick transformation of a cardboard sample, i.e., a representative volume element (RVE) into a flat plate structure with effective parameters describing the membrane and bending stiffness. On the other hand, the use of non-local sensitivity analysis makes it possible to find the relationship between the parameters of the paper and the load capacity of the packaging. The analytical procedures presented in our previous studies were used here to determine (1) the edge crush resistance, (2) critical load, and (3) the load capacity of corrugated cardboard packaging. The method proposed here allows for obtaining a comprehensive and hierarchical list of the parameters that play the most important role in the process of optimal packaging design.


Introduction
The growing ecological awareness and concern for the natural environment with the simultaneous, inevitable increase in the needs for the purchase, production, transport, and storage of diverse products, have forced the pursuit for environmentally friendly solutions, including easy to dispose of and, more importantly, recyclable packaging. At this point, it should be emphasized that corrugated cardboard boxes fit perfectly into this trend. The merits of such a package include also the simplicity of shaping through appropriately designed creasing, easy formation of openings, ventilation holes and perforations, and convenient color printing, which attracts the vendee in the case of shelf-ready packaging (SRP) or retail-ready packaging (RRP).
The high demand on the market has caused the intensive development of a distinct branch of industry, i.e., packaging production, and thus the rapid progress of scientific research in this field. The issue of strength evaluation of corrugated cardboard products is the subject of continuous, extensive studies. The corrugated cardboard consists of layers, therefore their proper selection and combination determine its relevant load-bearing capacity. Two characteristic in-plane directions of orthotropy indissolubly connect to the mechanical strength of the paperboard, namely the machine direction (MD), which is (RVE) is created [37]. A multiple scales asymptotic homogenization approach can be found in [38], whereas the asymptotic homogenization technique is in [39]. Homogenization of the corrugated board can be carried out in two variants, i.e., homogenization to one layer or homogenization of fluting to the inner layer of the laminate. Such a method has been extensively employed over the last years [40][41][42][43][44][45][46][47][48][49] due to substantial savings in computation time while maintaining the accuracy of the results.
Measurement from an experiment, in other words physical testing, is very common in the paper industry for the estimation of the corrugated board load-bearing capacity. A number of typical tests have been developed to standardize the process of the characterization of corrugated cardboard mechanical properties. Apart from the above-mentioned ECT and BCT, the bending test (BNT), the shear stiffness test (SST), and the torsional stiffness test (TST) are applicable. Bursting and humidity testing are also common. In order to collect the data from the outer surface of the sample during testing, video extensometry can be applied. This technique is based on the measurement of the relative distances between pairs of points traced across images captured at different load values [50,51]. Such a procedure is analogous, but simpler, to the digital image correlation (DIC), which is a full-field noncontact optical measurement procedure. Due to the very high accuracy of data acquisition, it is gaining more recognition in the field of experimental mechanics [29,[52][53][54][55][56][57][58].
One cannot forget that there are many factors that diminish the compression strength of corrugated paperboard boxes [59], such as openings, ventilation holes and perforations or indentations [60][61][62][63][64][65], shifted creases on the flaps [66], time and conditions of storage [67,68], and stacking load [47,69,70]. The influence of the box geometry as well as the composition and arrangement of the corrugated board layers on the change of the buckling force, edge crushing (ECT), and the compressive box strength resistance (BCT) are the elements that need to be considered when assessing the load capacity of the box. Another very important factor that needs to be mentioned because of its high impact on the cardboard strength is the profile of a corrugated web, labeled by letters A, B, C, E, and F. The difference between them is in flute height (A type is the tallest, F is the lowest), wavelength, and take-up factor, which is a quantification of the fluting length per unit length of the board [25]. For the most common packages, the cardboard with B and C flutes is applied; for big boxes A flute is applied and for the smallest one, e.g., cosmetic packaging, E and F fluting is applied. Moreover, while speaking of the double-wall corrugated cardboard, different combinations of fluting are applied, e.g., BC, BE, AE, FE, or EB, and vice versa.
A detailed analysis of the sensitivity in the context of observing changes in the value of the global buckling force, ECT, and BCT resulting from minor perturbations of the parameters of the calculation model, allows for answering the questions of many manufacturerswhich parameters of the paper/corrugated board have the greatest impact on the loadbearing capacity of the box? In the paper, after in-depth research, the authors present a comprehensive and hierarchical list of parameters that play the most important role in the process of optimal packaging design.

Material Parameters and Corrugated Cardboard Geometry
Corrugated cardboard is a fibrous material, which results in a strong orthotropy. The mechanical properties of the components depend on the orientation of the fibers in the layers. Two main directions can be determined: the machine direction (MD) along the fibers, and the cross direction (CD) across the fibers, perpendicular to the MD. In machine direction, the material is two to three times stiffer in tensile/bending and almost two times the shear/torsion than in the cross direction. The MD is along the waves (see Figure 1), which compensates for the weaker material properties in the cross direction.
The mechanical properties of the layers, such as the moduli of elasticity E 1 (E MD ) and E 2 (E CD ), and the compression strength SCT CD , can be determined based on the grammage of the paper from the Mondi technical data [71].
An example of calculating the stiffness modulus E 1 and E 2 from the MONDI specifications is as follows (see Table 1): where grm is grammage (g/m 2 ), thk is paper thickness (mm), TS MD is the tensile stiffness index in MD (Nmm/g), and TS CD is the tensile stiffness index in CD. It is assumed that the paper with a grammage of 100 g/m 2 has a thickness equal to 160 µm.
Materials 2022, 15, 720 4 of 18 the shear/torsion than in the cross direction. The MD is along the waves (see Figure 1), which compensates for the weaker material properties in the cross direction. The mechanical properties of the layers, such as the moduli of elasticity and , and the compression strength , can be determined based on the grammage of the paper from the Mondi technical data [71].
An example of calculating the stiffness modulus and from the MONDI specifications is as follows (see Table 1): where is grammage (g/m ), ℎ is paper thickness (mm), is the tensile stiffness index in MD (Nmm/g), and is the tensile stiffness index in CD. It is assumed that the paper with a grammage of 100 g/m has a thickness equal to 160 m.
Poisson's ratio can be computed from the empirical formula [72]: The in-plane shear stiffness is approximated by the formula [72]: The transverse shear stiffnesses can be determined from [73] as: Wave height, period, and take-up factor are selected based on the wave type [25], as specified in Table 2.   Poisson's ratio ν 12 can be computed from the empirical formula [72]: The in-plane shear stiffness G 12 is approximated by the formula [72]: The transverse shear stiffnesses can be determined from [73] as: Wave height, period, and take-up factor are selected based on the wave type [25], as specified in Table 2.

Homogenization Technique
To calculate the BCT value of the packaging, it is necessary to know the stiffness of the corrugated cardboard, which can be obtained using the numerical homogenization method. In the present study, the method based on the elastic energy equivalence between the full RVE model and the simplified shell model was applied. The described method was proposed by Biancolini [37], and was then extended to include transversal shear stiffness by Garbowski and Gajewski [74]. The RVE (representative volume element) is a small, periodic section of the 3D corrugated cardboard structure. All theoretical issues related to the constitutive model can be found in [74]. Only basic information is presented below.
The finite element formulation for a linear analysis can be expressed as follows: where K e is a statically condensed global stiffness matrix of the RVE, u e is a displacement vector, and F e is a vector of the nodal forces. Subscript e means values for external nodes only. In Figure 2, the finite element mesh and mesh nodes of RVE are shown.

Homogenization Technique
To calculate the BCT value of the packaging, it is necessary to know the stiffness o the corrugated cardboard, which can be obtained using the numerical homogenization method. In the present study, the method based on the elastic energy equivalence between the full RVE model and the simplified shell model was applied. The described method was proposed by Biancolini [37], and was then extended to include transversal shear stiff ness by Garbowski and Gajewski [74]. The RVE (representative volume element) is a small, periodic section of the 3D corrugated cardboard structure. All theoretical issues related to the constitutive model can be found in [74]. Only basic information is presented below.
The finite element formulation for a linear analysis can be expressed as follows: where is a statically condensed global stiffness matrix of the RVE, is a displace ment vector, and is a vector of the nodal forces. Subscript means values for externa nodes only. In Figure 2, the finite element mesh and mesh nodes of RVE are shown. Static condensation is the process of removing unknown degrees of freedom (DOF and leaving behind selected degrees of freedom, called principal DOFs (or primary un knowns). In this case, internal nodes are removed and external nodes are the principa DOFs. The condensed stiffness matrix of the external nodes can be computed from the following formula: , (6 where the subarrays are related to the external (subscript e) and internal (subscript i nodes: .

(7
After static condensation, the total elastic strain energy can be presented as the work of external forces on the corresponding displacements: Static condensation is the process of removing unknown degrees of freedom (DOF) and leaving behind selected degrees of freedom, called principal DOFs (or primary unknowns). In this case, internal nodes are removed and external nodes are the principal DOFs. The condensed stiffness matrix of the external nodes can be computed from the following formula: where the subarrays are related to the external (subscript e) and internal (subscript i) nodes: After static condensation, the total elastic strain energy can be presented as the work of external forces on the corresponding displacements: The balance of the total energy between the full RVE model and the simplified shell model is ensured by appropriate definition of displacements of external nodes and taking into account membrane and bending behavior [74]. The generalized displacements are related to the generalized strains on the RVE edges: where the H i matrix can be determined for each node (x i = x, y i = y, z i = z): The transformation matrix, H i , presented above, links the generalized displacement of each node on the boundary with the generalized strain vector of the corrugated board RVE model. More details on the derivation of such a matrix within the Kirchhoff−Love assumption can be found in [37], while the derivation within the Reissner−Mindlin theory can be found in [74].
Using the definition of the elastic strain energy: and considering a finite element subjected to bending, tension and transverse shear, the elastic internal energy can be represented by the formula: The stiffness matrix for a homogenized composite can be determined as follows: The described homogenization method turns the full 3D model into a simplified shell model, which allows for shortening the duration of the computations while maintaining high accuracy of the results.
The matrix H k is composed of matrices A, B, D, and R according to the following equation: where A contains tensile and shear stiffnesses, B contains coupling of tensile and bending stiffnesses, D contains bending and torsional stiffnesses, and R contains transverse shear stiffness. In cases of symmetrical cross-sections, the matrix B is the zero matrix. However, if the cross-section is asymmetric, non-zero terms appear in matrix B, which results from the coupling between bending/torsional curvatures and tensile/shear forces, and affects the values in the matrix D. Traditionally, this problem was solved by minimizing matrix B with an appropriate selection of the neutral axis. The uncoupled matrix D can be computed using the following equation: where D contains bending and torsional stiffnesses for non-zero matrix B.

Edge Crush Test
For the analytical determination of the BCT value, it is necessary to know the ECT value, which can be obtained by summing the strength of all layers, including the takeup factor: where α i is the take-up factor (see Table 2) and p i max is the maximum load of the i-th layer. The value of this load can be the compressive strength SCT CD i or critical load P i cr , whichever occurs first (see Figure 3): als 2022, 15, 720 7 of 18

Edge Crush Test
For the analytical determination of the BCT value, it is necessary to know the ECT value, which can be obtained by summing the strength of all layers, including the take-up factor: where is the take-up factor (see Table 2) and is the maximum load of the -th layer. The value of this load can be the compressive strength or critical load , whichever occurs first (see Figure 3): The critical load can be computed in many ways. An overview of the formulae was presented by Garbowski et al. [14]. For the determination of ECT, the critical load for rectangular orthotropic panels was calculated from the following: where: , , where is the number of half-waves for which reaches the minimum, is the modulus of elasticity in MD, is the modulus of elasticity in CD, and are Poisson's coefficients in the plane, is the in-plane shear modulus, and is the thickness of the layer. The critical load can be computed in many ways. An overview of the formulae was presented by Garbowski et al. [14]. For the determination of ECT, the critical load for rectangular orthotropic panels was calculated from the following: where: where m is the number of half-waves for which P L cr reaches the minimum, E 1 is the modulus of elasticity in MD, E 2 is the modulus of elasticity in CD, ν 12 and ν 21 are Poisson's coefficients in the plane, G 12 is the in-plane shear modulus, and t is the thickness of the layer.

Box Compression Test
After determining the ECT value and the critical load, the BCT value can be computed. For a rectangular package, the compressive strength can be obtained from the formula [14]: where P L cr and P B cr are the critical loads of the packaging walls, and γ L and γ B are the reduction coefficients, which can be computed from the following: The critical loads can be evaluated from Equation (18), but can also be obtained from the formula in which the transverse shear stiffness is also included: where: This approach is crucial when the corrugated cardboard is relatively thick (especially for B and C flutes, and double-walled cardboards), and its transverse shear modulus is low, e.g., due to unintentional crushing, during printing, or the lamination process.

The Non-Local Sensitivity Analysis
The non-local sensitivity in this study was carried out for several of the abovedescribed quantities, namely for edge crush resistance (ECT), critical load (P cr ) and box strength to static crushing (BCT). In each case, the parameters of the model are the basis weights of the individual layers of the corrugated board collected in vector x. If by h(x) we denote the quantity whose sensitivity is determined, then through small perturbations, ∆x i of the i-th layer grammage, a change in the determined quantity h(x ± e i ∆x i ) can be computed. Here, e i is a unit vector of i-th grammage in the parameter space. Then, by determining a numerical gradient through, e.g., the central difference, the sensitivity at a specific point in the parameter space (weights of the component papers) can be obtained. Therefore, the sensitivity can be described by the following formula: Non-local means here that the sensitivity is checked at many points in the model parameter space in order to build information about the gradients of the studied quantities in the full range of the parameter, and not locally at a specific point in this space. In Figure 4, an algorithm for the determination of non-local sensitivity is shown. In the flowchart, i is the iteration number and n is the number of all perturbing parameters.
Materials 2022, 15, 720 9 of 1 quantities in the full range of the parameter, and not locally at a specific point in this space In Figure 4, an algorithm for the determination of non-local sensitivity is shown. In the flowchart, is the iteration number and is the number of all perturbing parameters.

Results
Before proceeding to the sensitivity analysis, which is the main goal of this work, the quality of the homogenization method used here was first checked and validated. For this purpose, in the first step, an analysis of the influence of the number of finite elements on the homogenization result was carried out. Such a check also allows for determining the influence of static condensation on the quality of the solution. The example uses a simple three-layer model with the middle layer in the form of rectilinear sections (zigzag shape see Figure 5) in order to avoid the effect of discretization of the undulating layer (in which the number of elements affects the exact representation of the waveform).

Results
Before proceeding to the sensitivity analysis, which is the main goal of this work, the quality of the homogenization method used here was first checked and validated. For this purpose, in the first step, an analysis of the influence of the number of finite elements on the homogenization result was carried out. Such a check also allows for determining the influence of static condensation on the quality of the solution. The example uses a simple three-layer model with the middle layer in the form of rectilinear sections (zigzag shape, see Figure 5) in order to avoid the effect of discretization of the undulating layer (in which the number of elements affects the exact representation of the waveform).
flowchart, is the iteration number and is the number of all perturbing param

Results
Before proceeding to the sensitivity analysis, which is the main goal of this w quality of the homogenization method used here was first checked and validated. purpose, in the first step, an analysis of the influence of the number of finite elem the homogenization result was carried out. Such a check also allows for determi influence of static condensation on the quality of the solution. The example uses three-layer model with the middle layer in the form of rectilinear sections (zigza see Figure 5) in order to avoid the effect of discretization of the undulating layer (i the number of elements affects the exact representation of the waveform).    Table 3 summarizes all elastic material parameters and thicknesses of the individual layers used in the calculation model.  Table 4 shows the main stiffnesses obtained by homogenization using various numerical models, ranging from a model composed of only a few finite elements to a model composed of several hundred elements. Table 4 also lists the analytically determined bending stiffnesses D 11 and D 22 , and the tensile/compressive stiffnesses A 11 and A 22 . In the case of the direction 11 (MD), only flat layers were included in the analytical calculation of both stiffnesses (assuming that the undulating layer has no influence on the result in this direction). This assumption is true for a sine-shaped corrugated layer, and to a lesser extent for the zig-zag shape of the fluting. However, it was made to simplify the analytical calculations, knowing a-priori that the calculated values will be slightly lower than the real ones (as can be seen in Table 4). The presented results indicate a good agreement between the analytically calculated stiffnesses and the stiffnesses obtained as a result of homogenization. In addition, the convergence of the numerical models along with the increase in the number of elements is clearly noticeable, which makes it possible to conclude that the homogenization method is correct, and the influence of the static condensation method does not adversely affect the obtained results.
Returning to the main thread of this work, which are sensitivities, all values presented in the following tables and graphs are computed by Equation (32), where h becomes ECT, P cr , or BCT. Table 5 presents the sensitivity of ECT computed by Equation (16). The ECT value depends on the SCT in CD, the stiffness in CD and MD (indirectly through a critical load). Thus, ECT becomes the quantity described as a function of the grammage, see, e.g., Equation (1).
The sensitivity of the ECT was computed for the four flutes (B, C, E, and F) and the combinations of the basis weight of the corrugated board layers. To create combinations, the following ranges were adopted: liner grammage from 100 every 20 to 200 g/m 2 and fluting grammage from 80 every 20 to 160 g/m 2 .  Table 5 shows the specific sensitivity values for the analyzed waves. The second and third columns show the minimum and maximum values that can be obtained with the adopted ranges of the basis weights. The fourth and fifth columns show the average values of sensitivity from all of the analyzed combinations with liner and fluting perturbation. Table 6 shows the sensitivity of ECT, P cr , and BCT depending on the basis weight of the corrugated boards layers (liners and fluting). The presented sensitivities are the averaged values from 120 of different boxes with various dimensions. The smallest dimension of the box base is 100 × 100 and the largest considered dimension is 500 × 300, while the box height varies from 50 to 500. The results presented here are limited to the B wave only and three indices of corrugated cardboard, i.e., three-ply with a two liners grammage of 100 g/m 2 and a fluting grammage of 160 g/m 2 , marked as 100-160-100, and two subsequent grades marked as 160-80-160 and 140-100-140. Table 7 shows the sensitivity of ECT, P cr , and BCT depending on the basis weight of the corrugated boards layers (liners and fluting). The presented sensitivities are the averaged values from nine different boxes, which were higher than 400 mm with a base dimension lower than 200.  Table 8 shows the sensitivity of ECT, P cr , and BCT depending on the basis weight of the corrugated boards layers (liners and fluting). The presented sensitivities are the averaged values from 36 boxes lower than 150 mm.  Figure 6 shows the sensitivity of ECT, P cr , and BCT depending on the basis weight of the corrugated boards layers (liners and fluting). The presented sensitivities are the averaged values from 120 different boxes with various dimensions.      Figure 7 shows the sensitivity of P cr and BCT for high boxes, while Figure 8 shows the sensitivity for stocky boxes. In all cases, just three selected grades were considered and presented, namely: (a) 100-160-100, (b) 160-80-160, and (c) 140-100-140.   Figure 6 shows the sensitivity of ECT, , and BCT depending on the basis weight of the corrugated boards layers (liners and fluting). The presented sensitivities are the averaged values from 120 different boxes with various dimensions.     . The presented sensitivities are the averaged values from 120 different boxes with various dimensions. Figure 10 shows the sensitivity of and BCT for high boxes, while Figure 11 shows the sensitivity for stocky boxes. In all cases, just three selected grades were considered and presented, namely: (a) 100-160-100, (b) 160-80-160, and (c) 140-100-140.   Figure 9 shows the sensitivity of P cr and BCT depending on the bending stiffnesses D 11 and D 22 . The presented sensitivities are the averaged values from 120 different boxes with various dimensions. Figure 10 shows the sensitivity of P cr and BCT for high boxes, while Figure 11 shows the sensitivity for stocky boxes. In all cases, just three selected grades were considered and presented, namely: (a) 100-160-100, (b) 160-80-160, and (c) 140-100-140.  . The presented sensitivities are the averaged values from 120 different boxes with various dimensions. Figure 10 shows the sensitivity of and BCT for high boxes, while Figure 11 shows the sensitivity for stocky boxes. In all cases, just three selected grades were considered and presented, namely: (a) 100-160-100, (b) 160-80-160, and (c) 140-100-140.   . The presented sensitivities are the averaged values from 120 different boxes with various dimensions. Figure 10 shows the sensitivity of and BCT for high boxes, while Figure 11 shows the sensitivity for stocky boxes. In all cases, just three selected grades were considered and presented, namely: (a) 100-160-100, (b) 160-80-160, and (c) 140-100-140.

Discussion
The conducted analyses allowed for obtaining a complete picture of the sensitivity of both the ECT and the critical load (the main components of the packaging load capacity to static loads), as well as the BCT itself, to small perturbations of the grammage in individual layers of the corrugated cardboard. Table 5 presents the values of the ECT sensitivity to grammage change in the liners and fluting. It is clearly seen that ECT is more sensitive to grammage perturbation for the C wave than for the E wave. At a 10% grammage increase, the minimum sensitivity of ECT is comparable for both wave types (1.9% and 2.0% for the C and E wave, respectively), but the maximum sensitivity of ECT is greater for the C wave (10.9% to 5.2% for the C and E wave, respectively). On average, the ECT sensitivity for the C wave is 53% higher for liner grammage perturbation and 51% higher for fluting grammage perturbation than for the E wave. The main reason is quite obvious-it is because of the distances between the fluting crests and the wave geometry (see Table 2). In the case of the C wave, the wave period and its height are 8 mm and 3.61 mm, respectively, while for the E wave, they are 2.5 mm and 0.76 mm, respectively. Therefore, in the case of the C wave, the loss of stability of both liners and fluting occurs much more often, so that the maximum value of in Equation (17) is the critical load and not the SCT value.
Further observations regarding the obtained results of the ECT sensitivity analysis are as follows: 1. As the liner grammage increases, the ECT sensitivity to liner grammage perturbation increases. At the same time, ECT sensitivity to the second liner and fluting grammage perturbation (the grammage of which does not change) decreases.
2. The increase in fluting grammage reduces the ECT sensitivity to liners perturbation.
3. The lowest sensitivity is achieved with low liner grammage perturbation, where the other liner and fluting grammages are high. The greatest sensitivity occurs to the perturbation of the liner with a high grammage with a low grammage for the other liner and fluting.
The results presented in Tables 6-8 allow for drawing the conclusions that the sensitivities of the critical load of the longer and shorter walls of the box are similar and range between 3.05% and 4.48% when the basis weight of liner is changed by 10%, while when the basis weight of fluting is changed by 10%, it ranges between 1.73% and 3.00%. The shorter boxes have a slightly higher sensitivity of regarding the grammage of fluting for higher boxes, and varies from 1.94% to 3.00% for high boxes and from 2.90% to 3.45% for stocky boxes.
The situation is slightly different in the case of BCT sensitivity to changes in the grammage of the liners and fluting. The sensitivity of BCT, as for that of the ECT, is dependent

Discussion
The conducted analyses allowed for obtaining a complete picture of the sensitivity of both the ECT and the critical load (the main components of the packaging load capacity to static loads), as well as the BCT itself, to small perturbations of the grammage in individual layers of the corrugated cardboard. Table 5 presents the values of the ECT sensitivity to grammage change in the liners and fluting. It is clearly seen that ECT is more sensitive to grammage perturbation for the C wave than for the E wave. At a 10% grammage increase, the minimum sensitivity of ECT is comparable for both wave types (1.9% and 2.0% for the C and E wave, respectively), but the maximum sensitivity of ECT is greater for the C wave (10.9% to 5.2% for the C and E wave, respectively). On average, the ECT sensitivity for the C wave is 53% higher for liner grammage perturbation and 51% higher for fluting grammage perturbation than for the E wave. The main reason is quite obvious-it is because of the distances between the fluting crests and the wave geometry (see Table 2). In the case of the C wave, the wave period and its height are 8 mm and 3.61 mm, respectively, while for the E wave, they are 2.5 mm and 0.76 mm, respectively. Therefore, in the case of the C wave, the loss of stability of both liners and fluting occurs much more often, so that the maximum value of p max in Equation (17) is the critical load and not the SCT value.
Further observations regarding the obtained results of the ECT sensitivity analysis are as follows: 1.
As the liner grammage increases, the ECT sensitivity to liner grammage perturbation increases. At the same time, ECT sensitivity to the second liner and fluting grammage perturbation (the grammage of which does not change) decreases.

2.
The increase in fluting grammage reduces the ECT sensitivity to liners perturbation.

3.
The lowest sensitivity is achieved with low liner grammage perturbation, where the other liner and fluting grammages are high. The greatest sensitivity occurs to the perturbation of the liner with a high grammage with a low grammage for the other liner and fluting.
The results presented in Tables 6-8 allow for drawing the conclusions that the sensitivities of the critical load of the longer and shorter walls of the box are similar and range between 3.05% and 4.48% when the basis weight of liner is changed by 10%, while when the basis weight of fluting is changed by 10%, it ranges between 1.73% and 3.00%. The shorter boxes have a slightly higher sensitivity of P cr regarding the grammage of fluting for higher boxes, and varies from 1.94% to 3.00% for high boxes and from 2.90% to 3.45% for stocky boxes.
The situation is slightly different in the case of BCT sensitivity to changes in the grammage of the liners and fluting. The sensitivity of BCT, as for that of the ECT, is dependent on the configuration of the papers in the corrugated board. The sensitivity of BCT to changes in the grammage of liners by 10% ranges from 2.98% to 5.22%, while a change in the grammage of fluting by 10% results in a change in BCT from 3.53% to 6.08%. The difference between the sensitivity of the BCT for low and high boxes is very small and reaches a maximum of 5%.
The results presented in Figures 6-8 are the graphical representation of the data compiled in Tables 6-8. The results in Figures 9-11 exhibit the sensitivity of the critical load and BCT to the change in stiffness. There is a clear trend that the stiffness D 11 generates higher sensitivities for all quantities than D 22 . The change in P cr resulted from a change in stiffness D 11 by 10% ranges from 3.48% to 4.49%, while the change in P cr resulted from a change in stiffness D 22 ranges from 1.72% to 2.32%. The sensitivity of BCT regarding stiffness D 11 reaches 1.1%, while D 22 reaches 0.5%. The change of P cr resulted from changes of D 11 and D 22 are almost the same in the case of low boxes, while in the case of high boxes, the sensitivity of P cr regarding D 22 is several times lower than for D 11 .

Conclusions
Nowadays, it is very important for lightweight material to be used in the production of various structures, including corrugated cardboard packaging. Therefore, understanding and checking the impact of changing the grammage of individual layers of corrugated cardboard on the changes in its mechanical properties is crucial in the optimization process. The paper presents the results of extensive numerical analyzes carried out in order to determine the sensitivity of various quantities for determining the mechanical properties of corrugated cardboard and the packaging that is made of it. The study examined the sensitivity of edge crush resistance (ECT), critical load (P cr ), and static crushing resistance of packaging (BCT) to changes in the grammage of individual layers of corrugated cardboard. Based on the numerical and computational analyzes carried out here, it is possible to make decisions about changing the composition of the three-layer corrugated cardboard in an easier and more conscious way.