An Open-Source Simulation for Compression Analysis of Corrugated Boards

1. Introduction

Corrugated boards have become highly important in the packaging industry due to their lightweight, cost-effective nature and excellent ability to provide protection and cushioning during transportation [1,2,3]. In 2021, the global corrugated board market was valued at USD 199.8 billion, with projections estimating a rise to USD 254.5 billion by 2026, driven by a Compound Annual Growth Rate (CAGR) of 5.0% [1,4,5,6]. The significance of improving the knowledge of the mechanical behavior of corrugated boards—which display complex features such as anisotropy, stochasticity, contact behavior, and many failure modes—is highlighted by this increasing need [7].

Finite Element Method (FEM) simulations are frequently employed to predict the mechanical performance of corrugated boards, with popular commercial software packages like ANSYS and ABAQUS dominating the field. Studies in the literature [3,8,9,10] have demonstrated the effectiveness of these tools in analyzing compressive strength, deformation, and other mechanical properties of corrugated structures.

Recent studies in FEM-based modeling have focused on improving the accuracy of simulating corrugated boards’ compressive behavior. For instance, Park et al. [8] investigated the edgewise compression behavior of corrugated paperboard using FEM simulations. Their study highlighted the potential of FEM as a reliable tool for predicting mechanical performance, despite challenges in capturing the material’s inherent anisotropy and heterogeneities. Nordstrand [9] focused on the finite element analysis of paperboard packages under compressional force, examining buckling behavior and the impact of different loading conditions. The research was relevant in the observation of the structural integrity of paperboard packages and the effectiveness of FEM in predicting failure modes. Additionally, Ma et al. [10] investigated the buckling behavior of corrugated boards under compressive loads using FEM. The research aimed to understand the role of adhesive in the structural performance of such structures. Zaheer [11] developed a method for predicting the top-to-bottom compression strength of corrugated board containers using FEM. The model accommodated up to triple-wall corrugated board and was verified through extensive box compression tests. Similarly, Di Russo et al. [12] evaluated various wave configurations in corrugated boards through nonlinear FEM analysis, focusing on large displacements and plasticity. The research emphasized the computational intensity of detailed 3D numerical models and discussed homogenization processes to represent corrugated boards as an equivalent homogeneous plate, thereby reducing computational time while maintaining accuracy.

However, these approaches often come with significant drawbacks: fixed geometries that limit adaptability, reliance on expensive proprietary software, and a lack of alignment with the principles of sustainable and accessible design [3].

This paper introduces an open-source, automated, Pythonic approach for analyzing the compressive behavior of corrugated boards. The proposed framework has the intention of addressing the limitations of traditional FEM software, offering a flexible and cost-effective alternative that aligns with the principles of sustainable engineering. The authors want to tackle the research gap by providing a software framework that everyone can use, and whose results can be reproducible.

The main research questions explored in this work are as follows:

• How can an open-source Python-based approach accurately predict the compressive behavior of corrugated boards compared to commercial FEM software?
• To what extent does this framework improve adaptability and reduce dependency on fixed geometries?
• How can this methodology promote sustainable and accessible design in packaging applications?

The remainder of this paper is organized as follows: Section 2 provides a detailed review of the literature and highlights the challenges of using commercial FEM software for corrugated boards. Section 3 presents the methodology, including the Python-based implementation and the process of geometry adaptation. Section 4 discusses the results, comparing the proposed approach to traditional software in terms of accuracy, cost, and adaptability. Finally, Section 5 concludes the paper, summarizing the contributions and outlining future research directions.

2. Materials and Methods

This section is organized as follows: the first topic is the theoretical framework (strong and weak formulation), serving as a basis for the formulation of the corrugated board; after it, the specific formulation for the corrugated board, including the discretization and post-processing of the calculated displacements, is given; an introduction to Code Aster then follows it; and, finally, the methodology used to produce the results is presented.

2.1. Strong Formulation

The strong form of the equilibrium equation for a given domain $\Omega$ with boundary $\Gamma$ is expressed as:

$$\nabla ·\sigma +\mathbf{f}=\mathbf{0}\mathrm{in}\Omega ,$$

where:

• $\sigma =\mathbb{C}:\epsilon \left(\mathbf{u}\right)$ is the stress tensor,
• $\mathbb{C}$ is the elasticity tensor of the material,
• $\epsilon \left(\mathbf{u}\right)={\displaystyle \frac{1}{2}}\left[\nabla \mathbf{u}+{(\nabla \mathbf{u})}^{\mathrm{T}}\right]$ is the strain tensor,
• $\mathbf{f}$ is the body force acting on the domain.

2.2. Weak Formulation

The governing equation for the mechanical response of the corrugated board is derived from the principle of virtual work. For a given domain $\Omega$, the weak form of the equilibrium equation is:

$${\int }_{\Omega }\mathbb{C}:\epsilon \left(\mathbf{u}\right):\epsilon \left(\delta \mathbf{u}\right)d\Omega ={\int }_{\Omega }\mathbf{f}·\delta \mathbf{u}d\Omega +{\int }_{{\Gamma }_{\mathrm{top}}}P\delta u_{z}d\Gamma ,$$

(1)

where P is the compression load applied on the top liner at ${\Gamma }_{\mathrm{top}}$.

The boundary conditions are:

• Fixed bottom liner:

  $$u_x=0,u_y=0,u_z=0\mathrm{on}{\Gamma }_{\mathrm{bottom}},$$

  (2)
• Constrained flute extremes (side edges):

  $$u_x=0\mathrm{on}{\Gamma }_{\mathrm{left}}\cup {\Gamma }_{\mathrm{right}},$$

  (3)
• Compression load applied on the top liner:

  $${\sigma }_{zz}=-P\mathrm{on}{\Gamma }_{\mathrm{top}}.$$

  (4)

2.3. Contact Mechanics and Penalty Method

To model the interaction between the flute and liners, contact mechanics are incorporated. Based on the Kuhn-Tucker conditions, the non-penetration constraint is given by:

$$g_n\ge 0,{\lambda }_n\ge 0,{\lambda }_n{g}_n=0,$$

(5)

where:

• $g_n=({\mathbf{u}}_{\mathrm{liner}}-{\mathbf{u}}_{\mathrm{flute}})·\mathbf{n}$ is the normal gap between the surfaces,
• ${\lambda }_n$ is the contact force in the normal direction.

Using the penalty method, the contact constraint is approximated as:

$${\lambda }_n\approx k_p{g}_n,$$

(6)

where $k_p$ is the penalty stiffness parameter. A larger value of $k_p$ enforces the constraint more strictly but may lead to numerical instabilities. Substituting the penalty approximation into the weak formulation, the final form becomes:

$${\int }_{\Omega }\mathbb{C}:\epsilon \left(\mathbf{u}\right):\epsilon \left(\delta \mathbf{u}\right)d\Omega +{\int }_{{\Gamma }_{\mathrm{contact}}}{k}_p{g}_n\delta g_{n}d\Gamma ={\int }_{\Omega }\mathbf{f}·\delta \mathbf{u}d\Omega +{\int }_{{\Gamma }_{\mathrm{top}}}\mathbf{t}·\delta \mathbf{u}d\Gamma .$$

(7)

Here, the term ${\int }_{{\Gamma }_{\mathrm{contact}}}{k}_p{g}_n\delta g_{n}d\Gamma$ represents the penalty contribution for enforcing contact. The boundary conditions remain as described in the previous subsection. The penalty stiffness ($k_p$) controls the accuracy of contact enforcement and should be chosen based on the problem scale; the gap function ($g_n$) describes the relative normal displacement between the flute and liners; and the contact interface ${\Gamma }_{\mathrm{contact}}$ represents the region where the flute interacts with the liners. This penalty-based approach ensures a numerically stable approximation of the contact mechanics while maintaining computational efficiency.

2.4. Corrugated Board with Sine Wave Flute

To model the corrugated board with a sine wave flute geometry, the governing equations must incorporate the effects of the geometry on the mechanical response and contact mechanics. The flute is defined parametrically using a sine wave, influencing the domain $\Omega$ and the contact interface ${\Gamma }_{\mathrm{contact}}$. Figure 1 represents the structure aimed at being represented, with the defined domain $\Omega$ and the different boundaries $\Gamma$.

The geometry of the flute is given by:

$$y\left(x\right)=Asin\left({\displaystyle \frac{2\pi x}{\lambda }}\right),$$

(8)

where:

• A is the amplitude of the sine wave,
• $\lambda$ is the wavelength of the sine wave.

The unit normal vector to the flute surface is:

$$\mathbf{n}={\displaystyle \frac{1}{\sqrt{1+{\left(\frac{\partial y}{\partial x}\right)}^2}}}\left[\begin{array}{c}-{\displaystyle \frac{\partial y}{\partial x}}\\ 1\end{array}\right],$$

(9)

where:

$${\displaystyle \frac{\partial y}{\partial x}}=A{\displaystyle \frac{2\pi }{\lambda }}cos\left({\displaystyle \frac{2\pi x}{\lambda }}\right).$$

(10)

For each finite element e, the displacement field $\mathbf{u}$ and its variation $\delta \mathbf{u}$ are approximated using shape functions $N_i$:

$${\mathbf{u}}^e=\mathbf{N}{\mathbf{d}}^e,\delta {\mathbf{u}}^e=\mathbf{N}\delta {\mathbf{d}}^e,$$

(11)

where:

• $\mathbf{N}$ is the matrix of shape functions,
• ${\mathbf{d}}^e$ is the vector of nodal displacements for element e.

The strain-displacement relationship is:

$$\epsilon \left({\mathbf{u}}^e\right)=\mathbf{B}{\mathbf{d}}^e,\mathbf{B}={\displaystyle \frac{\partial \mathbf{N}}{\partial \mathbf{x}}},$$

(12)

where $\mathbf{B}$ is the strain-displacement matrix. The stiffness matrix for element e is given by:

$${\mathbf{K}}_e={\int }_{{\Omega }^e}{\mathbf{B}}^T\mathbb{C}\mathbf{B}d\Omega ,$$

(13)

where $\mathbb{C}$ is the elasticity tensor of the material. Using Gaussian quadrature, the stiffness matrix becomes:

$${\mathbf{K}}_e=\sum _{q=1}^{n_q}{\mathbf{B}}^T\left({\xi }_q\right)\mathbb{C}\mathbf{B}\left({\xi }_q\right)w_q\left|J\right|,$$

(14)

where:

• ${\xi }_q$ are the quadrature points,
• $w_q$ are the quadrature weights,
• J is the Jacobian determinant of the element transformation.

The penalty method introduces a contact stiffness matrix:

$${\mathbf{K}}_{\mathrm{contact}}^e={\int }_{{\Gamma }_{\mathrm{contact}}^e}{k}_p{\mathbf{N}}^T\mathbf{N}d\Gamma .$$

(15)

For a sine wave flute:

$${\int }_{{\Gamma }_{\mathrm{contact}}^e}d\Gamma ={\int }_{x_1}^{x_2}\sqrt{1+{\left(A{\displaystyle \frac{2\pi }{\lambda }}cos\left({\displaystyle \frac{2\pi x}{\lambda }}\right)\right)}^2}dx.$$

(16)

Thus:

$${\mathbf{K}}_{\mathrm{contact}}^e=\sum _{q=1}^{n_q}{k}_p{\mathbf{N}}^T\left({\xi }_q\right)\mathbf{N}\left({\xi }_q\right)w_q\sqrt{1+{\left(A{\displaystyle \frac{2\pi }{\lambda }}cos\left({\displaystyle \frac{2\pi x}{\lambda }}\right)\right)}^2}\left|J\right|.$$

(17)

The penalty load vector for contact is:

$${\mathbf{F}}_{\mathrm{contact}}^e={\int }_{{\Gamma }_{\mathrm{contact}}^e}{k}_p{\mathbf{N}}^T{g}_{n}d\Gamma ,$$

(18)

where the gap function $g_n$ is:

$$g_n=({\mathbf{u}}_{\mathrm{liner}}-{\mathbf{u}}_{\mathrm{flute}})·\mathbf{n}.$$

(19)

Substituting $d\Gamma$, this becomes:

$${\mathbf{F}}_{\mathrm{contact}}^e=\sum _{q=1}^{n_q}{k}_p{\mathbf{N}}^T\left({\xi }_q\right)g_n\left({\xi }_q\right)w_q\sqrt{1+{\left(A{\displaystyle \frac{2\pi }{\lambda }}cos\left({\displaystyle \frac{2\pi x}{\lambda }}\right)\right)}^2}\left|J\right|.$$

(20)

The global stiffness matrix and load vector are assembled by summing contributions from all elements:

$$\mathbf{K}=\sum _{e=1}^{n_e}{\mathbf{K}}_e+\sum _{e=1}^{n_e}{\mathbf{K}}_{\mathrm{contact}}^e,$$

(21)

$$\mathbf{F}=\sum _{e=1}^{n_e}{\mathbf{F}}_e+\sum _{e=1}^{n_e}{\mathbf{F}}_{\mathrm{contact}}^e.$$

(22)

The system of equations is solved for nodal displacements $\mathbf{d}$:

$$\mathbf{K}\mathbf{d}=\mathbf{F}.$$

(23)

Once the displacements $\mathbf{d}$ are computed:

• Strain is obtained as:

  $$\epsilon =\mathbf{B}\mathbf{d}.$$

  (24)
• Stress is obtained as:

  $$\sigma =\mathbb{C}\epsilon .$$

  (25)

2.5. Code Aster

Code Aster is a versatile open-source finite element analysis (FEA) software developed by Électricité de France (EDF) [13]. It is designed for solving a wide range of engineering problems, including structural mechanics, thermal analysis, and fluid-structure interactions. One of the primary advantages of Code Aster is its ability to handle both linear and nonlinear analyses, making it suitable for academic and industrial research [14]. However, as already mentioned, a distinguishing feature of Code Aster is the fact that it is free to use. It possesses robust support for scripting and automation through Python. This allows users to automate workflows, including geometry and mesh generation, boundary condition specification, and post-processing [14]. As a result, it is particularly suited for studies requiring iterative simulations and parametric analyses.

The software has been extensively used in various applications. For instance, it has been utilized for dynamic mooring simulations of floating wind turbines, demonstrating its capability in offshore engineering applications [15]. Additionally, Code Aster has been integrated into cloud-based platforms like SimScale, enhancing accessibility for structural analysis tasks [14].

A key strength of Code Aster lies in its extensive library of material models, solvers, and post-processing tools, allowing for a high degree of customization. It also integrates seamlessly with pre- and post-processing tools such as Salome-Meca and ParaView, enabling a streamlined workflow [14].

Due to its open-source nature, Code Aster is continuously being developed. This ensures regular updates, a wide range of documentation, and flexibility for researchers to adapt it to their needs.

2.6. Methodology

The methodology followed in this study integrates Python scripting, Code Aster for finite element simulations, and ParaView for post-processing. The entire workflow is depicted in Figure 2.

The steps include:

• Geometry and Mesh Definition: The geometry and mesh parameters are specified using Python scripts. This step ensures flexibility for parametric studies and rapid modifications.
• Material and Boundary Condition Specification: The material properties and boundary conditions are defined in Python scripts, which are then fed into Code Aster.
• FEA: Code Aster processes the input files to solve the FEA problem. The solver’s ability to handle large-scale computations and nonlinear mechanics is leveraged in this work.
• Post-Processing: The results are post-processed using ParaView, which provides advanced visualization capabilities for analyzing stress, strain, and deformation patterns.

While the process could be conducted manually, it is not as efficient for numerical optimization tasks that require the process to be automatic. Therefore, a tool called YACS is used.

YACS is an automation and workflow management tool integrated into the Salome-Meca platform. It provides a graphical and scriptable environment for defining, managing, and automating complex engineering processes. Designed to streamline simulations, YACS allows users to integrate input data, solvers, and post-processing tools into a seamless pipeline. This modular and customizable approach makes it ideal for iterative and parametric studies.

In this study, YACS is employed to automate the simulation process workflow to analyze the mechanical behavior of corrugated boards. Figure 3 illustrates the YACS workflow, where the process begins with a set of input parameters related to the geometry, material properties, and simulation conditions. The geometric definition of the corrugated board is based on a sine wave structure characterized by control points, weights, and scaling parameters. These parameters are defined to construct the sine wave flute, which is then scaled to a specified amplitude and wavelength to match the desired structure.

The inputs, such as the number of flutes, flute thickness, liner thickness, and material properties (e.g., Young’s moduli and Poisson’s ratios), are compiled into an XML file. This file is parsed and interpreted by YACS, which orchestrates the execution of the finite element simulations. The automation ensures that all parameters are consistently applied and outputs, such as stress distributions and load-deformation curves, are efficiently processed and stored for subsequent analysis and optimization.

In corrugated cardboard papers, the material exhibits three orthogonal directions: machine direction (MD), cross direction (CD), and through-thickness direction (ZD). The machine direction aligns with the manufacturing process, where the cardboard paper is stronger due to the orientation of the fibers in that direction. The cross direction runs perpendicular to the MD and is weaker since the fibers are not aligned. The through-thickness direction reflects the properties across the cardboard’s thickness and is the weakest direction. Table 1 lists the material properties of the paper used in this study.

The aluminum block in the simulation was modeled to closely represent real-life conditions. Its material properties were defined to be significantly stiffer than those of the corrugated board, listed in Table 2

Moreover, the flute exhibits a sine-wave design, which is characterized by periodicity and odd symmetry in the thickness direction. This property suggests that the stress distribution in the top and bottom liners should exhibit similar patterns along the length of the flute. However, as this study considers only two wave cycles to reduce computational time, the influence of geometric asymmetry at the flute ends is expected to be more pronounced, potentially affecting the overall stress distribution.

3. Results

3.1. Flute Profiles

Table 3 summarizes the key parameters for three standard flute profiles (A, C, and E). The thickness values are consistent with industry standards, where profile A exhibits the largest amplitude and wavelength, while profiles C and E have lower values of flute amplitude and wavelength.

3.2. Simulation results

Figure 4 to Figure 6 depict the deformation behavior of Geometry 1 at displacements of 0.625 mm, 1.25 mm, and 2.5 mm, respectively. At lower displacement (Figure 4), the stress is concentrated around the contact zones, showing limited deformation. As the displacement increases (Figure 5), the deformation spreads across the structure, indicating an increase in strain energy. By the time displacement reaches 2.5 mm (Figure 6), large-scale deformation is evident, with noticeable deformation with higher buckling risk.

Figure 7 and Figure 8 present the deformation of Geometry 2 and Geometry 3 at selected displacements. Geometry 2, with a displacement of 1.25 mm (Figure 7), shows stress concentration localized to the middle flute. In contrast, Geometry 3 at 0.6 mm displacement (Figure 8) demonstrates stress propagation predominantly at the ends of the contact areas. This suggests a simpler stress distribution for Geometry 3 compared to Geometry 2, likely due to the effect of the thickness, which is maintained constant for both geometries as it happens in real life.

Figure 9 and Figure 10 detail the pressure distribution in the top and bottom liners of Geometry 1. In the top liner (Figure 9), the pressure is more uniformly distributed, whereas the bottom liner (Figure 10) exhibits pronounced pressure peaks at the flute contacts. These variations are likely due to the anisotropic nature of the material and the interaction effects between the liners and the flute.

Comparatively, Figure 11 and Figure 12, which display the pressure distributions for Geometry 2, reveal a less uniform pressure spread in the top liner (Figure 11), with higher stress concentrations. The bottom liner (Figure 12) for Geometry 2 demonstrates a marginally broader pressure distribution, indicating different stress transmission pathways than Geometry 1.

In Figure 13 and Figure 14, the pressure distribution for Geometry 3 is depicted. The top liner (Figure 13) shows distributed pressure with fewer localized peaks compared to Geometry 2. However, the pressure sharing tends to be higher compared to Geometry 1. Similarly, the same happens with the bottom liner (Figure 14).

In general, the results presented highlight the quasi-symmetric nature of the pressure distributions observed across different flute geometries and displacements. While the sine-wave design of the flutes is about periodicity and balanced stress patterns, slight deviations in symmetry are evident due to geometric asymmetry at the flute ends. This asymmetry becomes more pronounced at higher displacements, where localized stress concentrations disrupt the otherwise regular distribution. The comparison between the top and bottom liners across all geometries indicates that, despite these deviations, the quasi-symmetric design remains effective in achieving predictable mechanical behavior, with the bottom liners generally experiencing higher stress peaks due to reaction mechanisms.

A key difference noted among the geometries is the stress distribution and deformation response. Geometry 1 tends to localize stress, making it susceptible to localized failure, particularly at high displacements. Geometry 2, while exhibiting better stress distribution than Geometry 1, still shows pronounced peaks, indicating areas prone to failure. Geometry 3, however, achieves a balanced pressure profile, likely due to the chosen geometric parameters.

While stresses are not able to be measured in the lab, especially for comparison purposes, the distributions presented in the results most likely reflect expectations: a thinner flute will cause more stress distribution in order to resist the high load, while thicker flutes cause less stress distribution. Moreover, higher stress values are located in the bottom liner, a region where, in the absence of geometric imperfections, they are known to be more likely to deform first due to reaction mechanisms.

3.3. Limitations

Using open-source finite element analysis, this study offers insightful information about the mechanical behavior of corrugated boards; however, several limitations should be noted.

One significant limitation lies in the material modeling. While the material properties were realistically assigned to the corrugated board and aluminum block, the model did not fully account for the anisotropic and nonlinear characteristics of the board under varying environmental conditions, such as moisture, temperature, and long-term loading. These factors are known to significantly influence the mechanical performance of corrugated structures and could introduce discrepancies when comparing simulation results to real-world applications.

Another limitation is the simplification of boundary and loading conditions. Although the study aimed to replicate real-life flat crush tests as closely as possible, practical applications often involve more complex scenarios, such as non-uniform loading, dynamic forces, or multi-directional stresses. The uniformity of the applied pressure and the fixed nature of the boundary conditions in the model may not adequately capture these real-world complexities. Additionally, the study focused on specific geometries of the flute, which, while informative, might not represent the entire spectrum of geometric variations used in industrial applications. This limits the generalizability of the findings to configurations outside the tested range. Future work could address these limitations by incorporating more comprehensive environmental and operational factors as well as exploring a broader variety of geometric designs.

3.4. Answering to the Research Questions

How can an open-source Python-based approach accurately predict the compressive behavior of corrugated boards compared to commercial FEM software?

The open-source Python-based approach demonstrated in this study achieves comparable accuracy to commercial FEM software by leveraging numerical methods and scripting capabilities in Code_Aster. The integration of Python scripting allows for the precise definition of material properties, boundary conditions, and geometric parameters, ensuring that the simulations are accurate and reproducible. While validation is limited, the results are highly consistent with the discussion done in [16] on the physical testing of paper, namely regarding corrugated board out-of-plane compression.

To what extent does this framework improve adaptability and reduce dependency on fixed geometries?

The framework significantly enhances adaptability by employing parametric modeling and automated workflows through YACS. Unlike commercial FEM software, which often relies on fixed geometries, the Python-based approach allows for the dynamic generation and modification of flute geometries. This flexibility is particularly advantageous for exploring a wide range of designs and configurations, as evidenced by the various geometries tested in this study. The ability to quickly adapt the model to new specifications reduces the time and cost associated with redesigns, offering a practical advantage for iterative design processes in the packaging industry.

How can this methodology promote sustainable and accessible design in packaging applications?

By relying on open-source tools, this methodology aligns with the principles of sustainability and accessibility. The avoidance of expensive proprietary software reduces financial barriers, making advanced FEM analyses accessible to a broader range of researchers and industry practitioners. Additionally, the integration of automated workflows minimizes computational resources and time, supporting the efficient evaluation of sustainable design alternatives. Even though the method does not fully capture the real complexities of a corrugated board geometry, such as heterogeneity and geometric imperfections, it is a good starting point for optimization tasks, especially when expected geometries (e.g., without imperfections or with a minimum number of them) are used.

4. Conclusions

This study demonstrated the effectiveness of an open-source Python-based framework for analyzing the compressive behavior of corrugated boards. By employing Code_Aster and parametric modeling through Python scripting, the methodology provided accurate and adaptable simulations comparable to commercial FEM software while significantly reducing the dependency on fixed geometries. The results highlighted the influence of geometric variations on stress distribution and deformation patterns.

The framework’s adaptability and automation capabilities were particularly evident in its ability to evaluate multiple flute geometries efficiently, demonstrating its potential for iterative design processes. Furthermore, the use of open-source tools aligns with the principles of sustainable and accessible engineering by minimizing financial barriers and promoting collaboration.

Future work should focus on addressing the limitations identified in this study, such as incorporating environmental factors and exploring a broader range of geometric configurations. Expanding the framework to include more advanced material models and dynamic loading conditions will enhance its applicability and accuracy. Representation of different geometries, including non-sine-wave flute types (e.g., FOAM-based geometries or using a different flute waving design), highly increases the applicability and need of use of the software due to the increase in the number of geometries to test. Moreover, in certain applications, the symmetry of the corrugated boards can significantly simplify the proposed software model. In practical scenarios, corrugated boards are often considered to have a much larger length compared to their thickness, effectively approximating an infinite number of waves. Under such conditions, the influence of asymmetries arising from the flute ends becomes negligible, further reinforcing the utility of the quasi-symmetric assumptions in the model.