Numerical Prediction of Equivalent Mechanical Properties of Corrugated Paperboard by 3D Finite Element Analysis

: This study is to obtain the equivalent (e ﬀ ective) mechanical properties of corrugated paperboard speciﬁcally that is widely applied in Korea for packaging of agricultural products. To analyze the equivalent properties of corrugated paperboard, ﬁnite element modeling, measurement of the reaction force, and superposition theory for the unit cell were applied. The stress-strain behavior obtained by applying the calculated equivalent mechanical properties to the simpliﬁed model of the corrugated paperboard is compared with experimental results. Nine equivalent mechanical properties governing the orthotropic behavior of corrugated paperboard were analyzed through ﬁnite element analysis (FEA). The stress-strain behavior of the corrugated paperboard experimentally showed elastic-plastic behavior, although the equivalent mechanical properties applied to the simpliﬁed model were elastic properties based on the theoretical approach, so that the ﬁnite element analysis results showed linearity. Therefore, when applying the calculated equivalent mechanical properties through FEA, the characteristics of the section where the strain only increases without the increase in the load due to the ﬂute should be taken into consideration.


Introduction
Corrugated paperboard is a sandwiched engineering structure, whose strength characteristic is in accordance with cross-sectional geometry of flute and material properties of linerboard and corrugating medium [1][2][3][4][5]. Corrugated paperboard is an environmental friendly packaging material with high stiffness and strength [2,[5][6][7][8]. In addition, it is efficient and easy to use to protect the contents of the package.
In Korea, linerboard and corrugating medium made by mixing imported old corrugated container (OCC) and Korea OCC at a certain ratio have been widely used for making corrugated paperboard. The corrugated paperboard with AB-flute double-wall (AB/F-DW) is widely used to maintain the strength of corrugated paperboard. Among corrugated paperboard produced in Korea in 2017, approximately 53.4% was the AB/F-DW corrugated paperboard, which was 1.15 times more than the used amount of SW (A/F and B/F) corrugated paperboard (KCCA, 2018). Recently, the application of corrugated paperboard in Korea has gradually expanded because a BB-flute double-wall (BB/F-DW) corrugated paperboard box having similar compressive strength to the AB/F-DW corrugated paperboard box has the advantage of reducing the package size. In addition, the use of corrugated paperboard has extended from the simple material of the box to materials of various cushion structures and pallets. The types of corrugated paperboard used in this study were DW-AB/F, DW-BB/F and SW-A/F. The board combinations and main specifications for the corrugated paperboard are summarized in Table 1.

Measurement of Material Properties of Linerboard and Corrugated Medium
The machine direction (MD) and cross direction (CD) tensile tests for linerboard and corrugated medium used for fabrication of the corrugated paperboard samples were conducted according to ISO 1924-2 [20] and 187 [21]. A universal testing machine (UTM) was used for the tensile test. To reduce the error due to pre-failure of the contact between the test sample and the grip during the tensile test, the shape of the tensile test samples was modified in accordance with ISO 6892-1 [22] (Figure 2).  Corrugated paperboard is an orthotropic material that has three symmetry planes for the elastic properties. It has different mechanical properties in each direction (machine direction (MD, x), the cross direction (CD, y), the thickness direction (TD, z) due to the shape of the flute. Alternatively, corrugated paperboard components such as linerboard and corrugated medium exhibit orthotropy due to the fiber orientation. The paper fibers are oriented in MD while forming the paper sheets of machine-made paper board.
The types of corrugated paperboard used in this study were DW-AB/F, DW-BB/F and SW-A/F. The board combinations and main specifications for the corrugated paperboard are summarized in Table 1.  due to the fiber orientation. The paper fibers are oriented in MD while forming the paper sheets of machine-made paper board. The types of corrugated paperboard used in this study were DW-AB/F, DW-BB/F and SW-A/F. The board combinations and main specifications for the corrugated paperboard are summarized in Table 1.

Measurement of Material Properties of Linerboard and Corrugated Medium
The machine direction (MD) and cross direction (CD) tensile tests for linerboard and corrugated medium used for fabrication of the corrugated paperboard samples were conducted according to ISO 1924-2 [20] and 187 [21]. A universal testing machine (UTM) was used for the tensile test. To reduce the error due to pre-failure of the contact between the test sample and the grip during the tensile test, the shape of the tensile test samples was modified in accordance with ISO 6892-1 [22] (Figure 2).

Measurement of Material Properties of Linerboard and Corrugated Medium
The machine direction (MD) and cross direction (CD) tensile tests for linerboard and corrugated medium used for fabrication of the corrugated paperboard samples were conducted according to ISO 1924-2 [20] and 187 [21]. A universal testing machine (UTM) was used for the tensile test. To reduce the error due to pre-failure of the contact between the test sample and the grip during the tensile test, the shape of the tensile test samples was modified in accordance with ISO 6892-1 [22] (Figure 2).
Appl. Sci. 2020, 10, x FOR PEER REVIEW 3 of 16 due to the fiber orientation. The paper fibers are oriented in MD while forming the paper sheets of machine-made paper board. The types of corrugated paperboard used in this study were DW-AB/F, DW-BB/F and SW-A/F. The board combinations and main specifications for the corrugated paperboard are summarized in Table 1.

Measurement of Material Properties of Linerboard and Corrugated Medium
The machine direction (MD) and cross direction (CD) tensile tests for linerboard and corrugated medium used for fabrication of the corrugated paperboard samples were conducted according to ISO 1924-2 [20] and 187 [21]. A universal testing machine (UTM) was used for the tensile test. To reduce the error due to pre-failure of the contact between the test sample and the grip during the tensile test, the shape of the tensile test samples was modified in accordance with ISO 6892-1 [22] (Figure 2).  Since the corrugated paperboard components (linerboard, corrugating medium) are also orthotropic, nine unknown parameters (E x , E y , E z , ν xy , ν xz , ν yz , G xy , G xz, and G yz ) should be determined.
However, some unknown parameters are not straightforward to measure because of the thin thickness of paper. The in-plane properties (E x and E y ,) can be easily obtained by standard tests, e.g., stress-strain curves (S-S curve). The approximate estimation method is commonly used to determine the rest of the unknown parameters. The elastic modulus in the out-of-plane direction and the shear modulus of each direction can be approximated by Equations (1) and (2) [23][24][25], respectively.

Theoretical Consideration for Modeling
The constitute equation of the stress-strain relationship, which represents the mechanical properties of orthotropic material, is given by: where ε 1 is the tensile strain in the x direction. ε 2 and ε 3 are the Poisson shrinkage strain in the y and z directions, respectively. E 1 is elastic modulus (Pa), G s are shear modulus (Pa), γ 12 are shear strain at a different direction (unit), ν ij is Poisson's ratio at i th number.
Since ν ij /E i = ν ji /E j and υ ij = − j /ε i , ν ij ν ji . Therefore, the number of engineering constants of constitute equation for describing the behavior of orthotropic material is 9. For the evaluation of the equivalent mechanical properties, 9 constants should be calculated from Equation (1).
To evaluate the equivalent mechanical properties of a structure, it is important to determine the unit cell representing the average response of the structure. The unit cell be selected as a repeating minimum basic unit in terms of the geometry and boundary conditions. Equivalent mechanical properties are evaluated through numerical experiments on the selected unit cell. In other words, the i-directional equivalent elastic modulus E i and the equivalent shear modulus G ij can be obtained by simulating the uniaxial tensile and in-plane shear load states that resemble the actual test, respectively.
To obtain the equivalent elastic modulus E 1 , it is assumed that a uniaxial tensile load is applied as shown in Figure 3. The uniaxial tensile load state is simulated by dividing into three sub-problems and then superimposing the results linearly. The solid line and dotted lines indicate the shape before and after deformation, respectively. Here, the tensile strain ε 1 in the x-direction and Poisson shrinkage strain ε 2 and ε 3 in the yand zdirection be obtained using the superposition method. Appl. Sci. 2020, 10, x FOR PEER REVIEW 5 of 16 In Sub-problem (1), the strain on the other directions is constrained to zero except for the tensile strain in the x-direction. For the Sub-problems (2) and (3), displacement boundary conditions are applied so that the deformation of the other side was 0, except for the tensile deformation in the yand z-directions. The reaction force at the interface can be obtained when FEA is performed on the sub-problems. To simulate the uniaxial tensile state, the obtained reaction forces are multiplied by weight coefficient α for Sub-problems (1), (2), and (3), respectively, and are combined as the following equation: where F is reaction force (N), α is weight (kg), P is applied load (kg) i indicates x, y, and z direction, j indicates sub-problems numbers. By using the determinant, the coefficient matrix [α]ji can be obtained, and the strain under an applied load [F]ij can be obtained using the matrix, as indicated in Equation (5).
The work Wx by the tensile load Px, and the internal energy Ux stored in the unit structure are estimated from the following equations: In Sub-problem (1), the strain on the other directions is constrained to zero except for the tensile strain in the x-direction. For the Sub-problems (2) and (3), displacement boundary conditions are applied so that the deformation of the other side was 0, except for the tensile deformation in the yand z-directions. The reaction force at the interface can be obtained when FEA is performed on the sub-problems. To simulate the uniaxial tensile state, the obtained reaction forces are multiplied by weight coefficient α for Sub-problems (1), (2), and (3), respectively, and are combined as the following equation: where F is reaction force (N), α is weight (kg), P is applied load (kg) i indicates x, y, and z direction, j indicates sub-problems numbers. By using the determinant, the coefficient matrix [α] ji can be obtained, and the strain under an applied load [F] ij can be obtained using the matrix, as indicated in Equation (5).
The work W x by the tensile load P x, and the internal energy U x stored in the unit structure are estimated from the following equations: Appl. Sci. 2020, 10, 7973 6 of 16 Since the work W x by the tensile load P x, is equal to the internal energy U x stored in the unit structure, the equivalent elastic modulus E 1 for Sub-problem (1) can be calculated from Equation (8). E 2 and E 3 for Sub-problem (2) and (3) are also obtained from Equations (9) and (10).
where ∆x, ∆y, and ∆z represent the dimensions of the unit cell. Equivalent Poisson's ratio is given by: The calculation of the equivalent shear coefficients for orthotropic materials should be conducted separately because the shear stress only produces shear deformation of the component. After performing FEA under the boundary condition in which the unit cell is in the state of pure shear deformation, the equivalent shear modulus can be calculated from Equation (12) derived from the relation between the work W by the shear load and the energy U stored in the unit cell ( Figure 4).
Appl. Sci. 2020, 10, x FOR PEER REVIEW 6 of 16 Since the work Wx by the tensile load Px, is equal to the internal energy Ux stored in the unit structure, the equivalent elastic modulus E1 for Sub-problem (1) can be calculated from Equation (8). E2 and E3 for Sub-problem (2) and (3) are also obtained from Equations (9) and (10).
where Δx, Δy, and Δz represent the dimensions of the unit cell. Equivalent Poisson's ratio is given by: The calculation of the equivalent shear coefficients for orthotropic materials should be conducted separately because the shear stress only produces shear deformation of the component. After performing FEA under the boundary condition in which the unit cell is in the state of pure shear deformation, the equivalent shear modulus can be calculated from Equation (12) derived from the relation between the work W by the shear load and the energy U stored in the unit cell ( Figure 4).
With the same calculation method for G12, the equivalent shear modulus G23 and G31 can be derived from Equations (15) and (16), respectively.  With the same calculation method for G 12 , the equivalent shear modulus G 23 and G 31 can be derived from Equations (15) and (16), respectively.

Model Development
Through the preliminary analysis of various repetition models, modeling for the unit cell in the MD of the corrugated paperboards used in this study was three times the minimum repetition level based on the integer multiples of the A/F and B/F wavelength. This accounted for the dimension that can overcome the excessive reaction force at the constraint points of both ends in the x-direction. The dimensions were 18 mm for AB/F-DW and BB/F-DW corrugated paperboard, and 27 mm for A/F-SW corrugated paperboard. However, the dimension of the unit cells in the CD was taken as unit lengths on all corrugated paperboard samples and modeled as solid state considering the thickness direction.
MiDAS NFX (2018R2, MiDAS IT) software was used for FE modeling [26]. The interface between the flute and liner of the corrugated paperboard was joined to simplify modeling. The geometrical shape of the flute was modeled as a sine function based on the data shown in Table 1. The FE modeling for the unit cells of the corrugated paperboard samples is indicated in Figure 5. The hexahedral mesh was used due to the shape of corrugated papers. The computational domains for AB/F-DW, BB/F-DW, and A/F-SW were discretized into 156,825 nodes and 102,400 mesh elements, 148,767 nodes and 97,200 mesh elements, and 144,228 nodes and 94,200 mesh elements, respectively.

Model Development
Through the preliminary analysis of various repetition models, modeling for the unit cell in the MD of the corrugated paperboards used in this study was three times the minimum repetition level based on the integer multiples of the A/F and B/F wavelength. This accounted for the dimension that can overcome the excessive reaction force at the constraint points of both ends in the x-direction. The dimensions were 18 mm for AB/F-DW and BB/F-DW corrugated paperboard, and 27 mm for A/F-SW corrugated paperboard. However, the dimension of the unit cells in the CD was taken as unit lengths on all corrugated paperboard samples and modeled as solid state considering the thickness direction.
MiDAS NFX (2018R2, MiDAS IT) software was used for FE modeling [26]. The interface between the flute and liner of the corrugated paperboard was joined to simplify modeling. The geometrical shape of the flute was modeled as a sine function based on the data shown in Table 1. The FE modeling for the unit cells of the corrugated paperboard samples is indicated in Figure 5 Table 1).
The boundary conditions applied when analyzing the three vertical displacements and shear displacements were applied according to the methods shown in Figures 3 and 4. As a representative example, the AB/F-DW case is shown in Figures 6 and 7.

Model Development
Through the preliminary analysis of various repetition models, modeling for the unit cell in the MD of the corrugated paperboards used in this study was three times the minimum repetition level based on the integer multiples of the A/F and B/F wavelength. This accounted for the dimension that can overcome the excessive reaction force at the constraint points of both ends in the x-direction. The dimensions were 18 mm for AB/F-DW and BB/F-DW corrugated paperboard, and 27 mm for A/F-SW corrugated paperboard. However, the dimension of the unit cells in the CD was taken as unit lengths on all corrugated paperboard samples and modeled as solid state considering the thickness direction.
MiDAS NFX (2018R2, MiDAS IT) software was used for FE modeling [26]. The interface between the flute and liner of the corrugated paperboard was joined to simplify modeling. The geometrical shape of the flute was modeled as a sine function based on the data shown in Table 1. The FE modeling for the unit cells of the corrugated paperboard samples is indicated in Figure 5

Model Suitability
The suitability of the analysis model was confirmed by calculating the Mesh Convergence Error using Equation (1) [27]: here, is the RMS (root mean square) value of the stress obtained over the entire area of the model with respect to the volume, and and represent the volume of the element and the analysis model, respectively.
≈ (19) where , and denote the stress, weight (1/N) and the number of adjacent elements, respectively.
Representative examples of mesh convergence errors in vertical and shear deformation analysis are shown in Figures 8 and 9, respectively, and the values for each analysis are summarized in Table 2.

Model Suitability
The suitability of the analysis model was confirmed by calculating the Mesh Convergence Error using Equation (1) [27]: here, σ ave is the RMS (root mean square) value of the stress obtained over the entire area of the model with respect to the volume, and ν i and ν t represent the volume of the element and the analysis model, respectively.
where σ i , w i and N denote the stress, weight (1/N) and the number of adjacent elements, respectively. Representative examples of mesh convergence errors in vertical and shear deformation analysis are shown in Figures 8 and 9, respectively, and the values for each analysis are summarized in Table 2.

Model Suitability
The suitability of the analysis model was confirmed by calculating the Mesh Convergence Error using Equation (1) [27]: here, is the RMS (root mean square) value of the stress obtained over the entire area of the model with respect to the volume, and and represent the volume of the element and the analysis model, respectively.

≈
where , and denote the stress, weight (1/N) and the number of adjacent elements, respectively.
Representative examples of mesh convergence errors in vertical and shear deformation analysis are shown in Figures 8 and 9, respectively, and the values for each analysis are summarized in Table 2.    Figure 10 shows the stress-strain(S-S) curves of the uniaxial tensile test on the linerboard and corrugating medium of corrugated paperboard samples. The high anisotropic and non-linear mechanical behavior was observed. In ISO 1924-2 [20], the elastic modulus of the paper is defined as the maximum slope on the S-S curve. However, when ISO 1924-2 method was applied to the test results of the corrugated paperboard components used in this study, very high elastic modulus was calculated. If the elastic modulus is applied to FEA, the result from FEA indicates a significant difference from the actual behavior of the corrugated paperboard. Thus, in this study, the elastic modulus in the MD and CD of the corrugated paperboard components was represented by the secant modulus, which was the slope of the line connecting the origin point and the rupture point.    Figure 10 shows the stress-strain(S-S) curves of the uniaxial tensile test on the linerboard and corrugating medium of corrugated paperboard samples. The high anisotropic and non-linear mechanical behavior was observed. In ISO 1924-2 [20], the elastic modulus of the paper is defined as the maximum slope on the S-S curve. However, when ISO 1924-2 method was applied to the test results of the corrugated paperboard components used in this study, very high elastic modulus was calculated. If the elastic modulus is applied to FEA, the result from FEA indicates a significant difference from the actual behavior of the corrugated paperboard. Thus, in this study, the elastic modulus in the MD and CD of the corrugated paperboard components was represented by the secant modulus, which was the slope of the line connecting the origin point and the rupture point.   Figure 10 shows the stress-strain(S-S) curves of the uniaxial tensile test on the linerboard and corrugating medium of corrugated paperboard samples. The high anisotropic and non-linear mechanical behavior was observed. In ISO 1924-2 [20], the elastic modulus of the paper is defined as the maximum slope on the S-S curve. However, when ISO 1924-2 method was applied to the test results of the corrugated paperboard components used in this study, very high elastic modulus was calculated. If the elastic modulus is applied to FEA, the result from FEA indicates a significant difference from the actual behavior of the corrugated paperboard. Thus, in this study, the elastic modulus in the MD and CD of the corrugated paperboard components was represented by the secant modulus, which was the slope of the line connecting the origin point and the rupture point.  The other physical properties (E z , G xy , G xz , G yz ) were calculated from Equations (1) and (2). The ν xy , ν xz, and ν yz values were approximated and determined according to the material values similar to corrugated paperboard samples (0.33 for the corrugating medium and 0.34 for the linerboard) used in previous studies [4,10]. The physical properties of corrugated paperboard samples are summarized in Table 3 were applied to FEA.

Calculation of Equivalent Material Properties of Corrugated Paperboards
In order to obtain the equivalent mechanical properties of the corrugated paperboard depending on flute type, the physical properties (Table 3) of linerboard and corrugating medium were applied to the FE model for the unit cells ( Figure 5) of corrugated paperboard samples. Since the contact surface conditions of FEA were the contact areas between the linerboard and the flute of corrugated paperboard, the number of contact surfaces in DW and SW corrugated paperboards was set to 4 and 2, respectively. As aforementioned, the interface between the flute and liner of the corrugated paperboard was considered to be joined for the convenience of FEA.
First, to obtain the equivalent elastic modulus, the reaction force was measured through FEA with each constraint given to the three sub-problems ( Figure 3). The reaction forces (F xi , F yi , and F zi ) were measured in the state in which forced displacement of 1 mm was applied in the x, y, and z-directions, respectively. In this case, all degrees of freedom were constrained in all directions except for the direction in which the forced displacement was applied: the y and z-directions when measuring F xi , the x, and z-directions when measuring F yi , and the x and y-directions when measuring F zi .
The measured directional reaction force for unit cells of corrugated paperboard samples through FEA are shown in Figures 11-13. In addition, the measurement results of the reaction force for the three sub-problems are provided in Table 4. The measured reaction forces through FEA were not significantly different from the measured reaction forces in tension and compression test. The corrugated paper board is mainly subjected to compressive force when it is applied to the product packaging and cushioning. Therefore, the reaction force measured under compressive load was applied to the analysis of the equivalent elastic modulus for the corrugated paperboard samples. The other physical properties (Ez, Gxy, Gxz, Gyz) were calculated from Equations (1) and (2). The νxy, νxz, and νyz values were approximated and determined according to the material values similar to corrugated paperboard samples (0.33 for the corrugating medium and 0.34 for the linerboard) used in previous studies [4,10]. The physical properties of corrugated paperboard samples are summarized in Table 3 were applied to FEA. Table 3. Material properties of corrugated paperboard components applied to FEA.

Calculation of Equivalent Material Properties of Corrugated Paperboards
In order to obtain the equivalent mechanical properties of the corrugated paperboard depending on flute type, the physical properties (Table 3) of linerboard and corrugating medium were applied to the FE model for the unit cells ( Figure 5) of corrugated paperboard samples. Since the contact surface conditions of FEA were the contact areas between the linerboard and the flute of corrugated paperboard, the number of contact surfaces in DW and SW corrugated paperboards was set to 4 and 2, respectively. As aforementioned, the interface between the flute and liner of the corrugated paperboard was considered to be joined for the convenience of FEA.
First, to obtain the equivalent elastic modulus, the reaction force was measured through FEA with each constraint given to the three sub-problems ( Figure 3). The reaction forces (Fxi, Fyi, and Fzi) were measured in the state in which forced displacement of 1 mm was applied in the x, y, and z-directions, respectively. In this case, all degrees of freedom were constrained in all directions except for the direction in which the forced displacement was applied: the y and z-directions when measuring Fxi, the x, and z-directions when measuring Fyi, and the x and y-directions when measuring Fzi.
The measured directional reaction force for unit cells of corrugated paperboard samples through FEA are shown in Figures 11-13. In addition, the measurement results of the reaction force for the three sub-problems are provided in Table 4. The measured reaction forces through FEA were not significantly different from the measured reaction forces in tension and compression test. The corrugated paper board is mainly subjected to compressive force when it is applied to the product packaging and cushioning. Therefore, the reaction force measured under compressive load was applied to the analysis of the equivalent elastic modulus for the corrugated paperboard samples.   The equivalent elastic modulus for the corrugated paperboard samples could be calculated by applying the reaction forces listed in Table 4 to Equations (2) and (5)- (7). The calculated results are summarized in Table 5. For example, the equivalent elastic modulus for AB/F-DW corrugated paperboard was calculated as follows. The coefficient ([α]ji, relative ratio) was calculated by substituting the reaction forces obtained from FEA for each sub-problem into Equation (2). The coefficient for each sub-problem was calculated by applying the cross-sectional area and the strain of AB/F-DW corrugated paperboard (Table 5) and arbitrary load ([P]ji = 100 N). The calculated coefficients are listed in Table 6.   The equivalent elastic modulus for the corrugated paperboard samples could be calculated by applying the reaction forces listed in Table 4 to Equations (2) and (5)- (7). The calculated results are summarized in Table 5. For example, the equivalent elastic modulus for AB/F-DW corrugated paperboard was calculated as follows. The coefficient ([α]ji, relative ratio) was calculated by substituting the reaction forces obtained from FEA for each sub-problem into Equation (2). The coefficient for each sub-problem was calculated by applying the cross-sectional area and the strain of AB/F-DW corrugated paperboard (Table 5) and arbitrary load ([P]ji = 100 N). The calculated coefficients are listed in Table 6.   The equivalent elastic modulus for the corrugated paperboard samples could be calculated by applying the reaction forces listed in Table 4 to Equations (2) and (5)- (7). The calculated results are summarized in Table 5. For example, the equivalent elastic modulus for AB/F-DW corrugated paperboard was calculated as follows. The coefficient ([α] ji , relative ratio) was calculated by substituting the reaction forces obtained from FEA for each sub-problem into Equation (2). The coefficient for each sub-problem was calculated by applying the cross-sectional area and the strain of AB/F-DW corrugated paperboard (Table 5) and arbitrary load ([P] ji = 100 N). The calculated coefficients are listed in Table 6.  The procedure for obtaining the equivalent shear modulus was similar to that for the equivalent elastic modulus. The constraint was given to three sub-problems, and the reaction force in each direction was measured through FEA. The reaction forces F xyi , F yzi, and F zxi were measured in the state in which forced displacement of 1 mm was applied in the x-direction of the upper face xz (fixing lower face xz with xy plane in front), the y-direction of the upper face xy (fix lower face xy with yz plane in front), and the x-direction of the upper face xy (fix lower face xy with xz plane in front), respectively. In these cases, all degrees of freedom were constrained in all directions except for the direction in which the forced displacement was applied: the y and z-directions when measuring F xyi , the x and z-directions when measuring F yzi , and the z and y-directions when measuring F zxi . Figures 14-16 show the directional reaction force measurement for the unit cell of corrugated paperboard samples, and the measurement results are provided in Table 7. The equivalent shear modulus was calculated by applying the reaction force listed in Table 7 to Equations (11)- (13). For example, the equivalent shear modulus for AB/F-DW corrugated paperboard was calculated as follows. The calculation of the equivalent shear modulus was much simpler than the calculation of the longitudinal modulus. Instead of combining nine reaction components by applying a load in each direction to each sub-problem, the equivalent shear modulus be easily obtained from the reaction force measured by applying a shearing force to each model. The cross-sectional area and shear deformation angle of AB/F-DW (Table 8) were substituted in Equations (11)-(13) to obtain the equivalent shear modulus.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 12 of 16 The procedure for obtaining the equivalent shear modulus was similar to that for the equivalent elastic modulus. The constraint was given to three sub-problems, and the reaction force in each direction was measured through FEA. The reaction forces Fxyi, Fyzi, and Fzxi were measured in the state in which forced displacement of 1 mm was applied in the x-direction of the upper face xz (fixing lower face xz with xy plane in front), the y-direction of the upper face xy (fix lower face xy with yz plane in front), and the x-direction of the upper face xy (fix lower face xy with xz plane in front), respectively. In these cases, all degrees of freedom were constrained in all directions except for the direction in which the forced displacement was applied: the y and z-directions when measuring Fxyi, the x and z-directions when measuring Fyzi, and the z and y-directions when measuring Fzxi. Figures 14-16 show the directional reaction force measurement for the unit cell of corrugated paperboard samples, and the measurement results are provided in Table 7. The equivalent shear modulus was calculated by applying the reaction force listed in Table 7 to Equations (11)- (13). For example, the equivalent shear modulus for AB/F-DW corrugated paperboard was calculated as follows. The calculation of the equivalent shear modulus was much simpler than the calculation of the longitudinal modulus. Instead of combining nine reaction components by applying a load in each direction to each sub-problem, the equivalent shear modulus be easily obtained from the reaction force measured by applying a shearing force to each model. The cross-sectional area and shear deformation angle of AB/F-DW (Table 8) were substituted in Equations (11)- (13) to obtain the equivalent shear modulus.    The procedure for obtaining the equivalent shear modulus was similar to that for the equivalent elastic modulus. The constraint was given to three sub-problems, and the reaction force in each direction was measured through FEA. The reaction forces Fxyi, Fyzi, and Fzxi were measured in the state in which forced displacement of 1 mm was applied in the x-direction of the upper face xz (fixing lower face xz with xy plane in front), the y-direction of the upper face xy (fix lower face xy with yz plane in front), and the x-direction of the upper face xy (fix lower face xy with xz plane in front), respectively. In these cases, all degrees of freedom were constrained in all directions except for the direction in which the forced displacement was applied: the y and z-directions when measuring Fxyi, the x and z-directions when measuring Fyzi, and the z and y-directions when measuring Fzxi. Figures 14-16 show the directional reaction force measurement for the unit cell of corrugated paperboard samples, and the measurement results are provided in Table 7. The equivalent shear modulus was calculated by applying the reaction force listed in Table 7 to Equations (11)- (13). For example, the equivalent shear modulus for AB/F-DW corrugated paperboard was calculated as follows. The calculation of the equivalent shear modulus was much simpler than the calculation of the longitudinal modulus. Instead of combining nine reaction components by applying a load in each direction to each sub-problem, the equivalent shear modulus be easily obtained from the reaction force measured by applying a shearing force to each model. The cross-sectional area and shear deformation angle of AB/F-DW (Table 8) were substituted in Equations (11)- (13) to obtain the equivalent shear modulus.   Equivalent Poisson's ratio for the target corrugated paperboards was calculated using Equation (8). The calculated elastic modulus, shear modulus, and Poisson's ratio in all directions are summarized in Table 9. Table 9. Equivalent material properties for the corrugated paperboard samples.

Stress-Strain Behavior Using Equivalent Materials Properties of Corrugated Paperboard
In order to perform FEA on the flat crush resistance, the calculated equivalent mechanical properties of the corrugated paperboard samples were applied to the simplified FE models ( Figure  17). The simplified model was discretized into 156,825 nodes and 102,400 mesh elements. The S-S behavior (flat crush resistance) of the simplified FE model was compared with experimental data of corrugated paperboard samples obtained by following ISO 3035 [27] (Figure 18). In case of the flat crush test of corrugated paperboard samples, the S-S curve was measured when a circular metal plate that has 90.6 mm in diameter and was attached the upper cross-head of tester exerted compressive force on corrugated paperboard sample (length × width = 30 × 30 cm). The loading rate was 12.5 ± 2.5 mm/min. Sandpaper was attached to the upper and lower compression plates of the tester in contact with the test piece to reduce the experimental error by the swelling of the flute in the test piece of the  Equivalent Poisson's ratio for the target corrugated paperboards was calculated using Equation (8). The calculated elastic modulus, shear modulus, and Poisson's ratio in all directions are summarized in Table 9.

Stress-Strain Behavior Using Equivalent Materials Properties of Corrugated Paperboard
In order to perform FEA on the flat crush resistance, the calculated equivalent mechanical properties of the corrugated paperboard samples were applied to the simplified FE models (Figure 17). The simplified model was discretized into 156,825 nodes and 102,400 mesh elements. The S-S behavior (flat crush resistance) of the simplified FE model was compared with experimental data of corrugated paperboard samples obtained by following ISO 3035 [27] (Figure 18). In case of the flat crush test of corrugated paperboard samples, the S-S curve was measured when a circular metal plate that has 90.6 mm in diameter and was attached the upper cross-head of tester exerted compressive force on corrugated paperboard sample (length × width = 30 × 30 cm). The loading rate was 12.5 ± 2.5 mm/min. Sandpaper was attached to the upper and lower compression plates of the tester in contact with the test piece to reduce the experimental error by the swelling of the flute in the test piece of the corrugated paperboard. However, the FEA for the simplified model of corrugated paperboard samples was reproduced by simulating the same conditions as the experiment.
The flat crush behavior of the corrugated paperboard samples showed elastic-plastic behavior when it exceeded a certain strain; however, the calculated equivalent value was the result of the theoretical approach assuming the elastic region.
As expected, the difference between the simulated and experimental values was observed in the region beyond the elastic region. In the actual behavior of the corrugated paperboard, the liner of the corrugated paperboard was deformed in the shear direction because of the change of the flute shape. Therefore, the experiment showed more significant deformation at smaller loads than ideal conditions such as FEA ( Figure 18). As a result, when the equivalent mechanical properties calculated by FEA were used to simulate S-S behavior of corrugated paperboard, it was necessary consider the characteristics of the region where the strain increased without an increase in load due to the flute of the corrugated paperboard.

Conclusions
Simulation techniques such as FEA must avoid the experimental factors in the design of various structures made of corrugated paperboard and to achieve the optimum design. However, due to the complexity of the corrugated paperboard itself in FEA for these structures, a significant time and effort are required for modeling and analysis. Therefore, it is necessary to analyze the simplified model and the equivalent mechanical properties for corrugated paperboard by flute type. when it exceeded a certain strain; however, the calculated equivalent value was the result of the theoretical approach assuming the elastic region.
As expected, the difference between the simulated and experimental values was observed in the region beyond the elastic region. In the actual behavior of the corrugated paperboard, the liner of the corrugated paperboard was deformed in the shear direction because of the change of the flute shape. Therefore, the experiment showed more significant deformation at smaller loads than ideal conditions such as FEA ( Figure 18). As a result, when the equivalent mechanical properties calculated by FEA were used to simulate S-S behavior of corrugated paperboard, it was necessary consider the characteristics of the region where the strain increased without an increase in load due to the flute of the corrugated paperboard.

Conclusions
Simulation techniques such as FEA must avoid the experimental factors in the design of various structures made of corrugated paperboard and to achieve the optimum design. However, due to the complexity of the corrugated paperboard itself in FEA for these structures, a significant time and effort are required for modeling and analysis. Therefore, it is necessary to analyze the simplified model and the equivalent mechanical properties for corrugated paperboard by flute type. When the upper surface of the model with a diameter of 90.6 mm was forcibly displaced at the same loading rate as the experiment in the z-direction, the reaction force transmitted to the lower surface of the model and S-S curves were measured and recorded. In FEA, the lower surface of the model was restricted to the translational and rotational motion in the x, y, and z-directions; however, the upper surface was restricted only to translational motion in x and y-directions.
The flat crush behavior of the corrugated paperboard samples showed elastic-plastic behavior when it exceeded a certain strain; however, the calculated equivalent value was the result of the theoretical approach assuming the elastic region.
As expected, the difference between the simulated and experimental values was observed in the region beyond the elastic region. In the actual behavior of the corrugated paperboard, the liner of the corrugated paperboard was deformed in the shear direction because of the change of the flute shape. Therefore, the experiment showed more significant deformation at smaller loads than ideal conditions such as FEA ( Figure 18). As a result, when the equivalent mechanical properties calculated by FEA were used to simulate S-S behavior of corrugated paperboard, it was necessary consider the characteristics of the region where the strain increased without an increase in load due to the flute of the corrugated paperboard.

Conclusions
Simulation techniques such as FEA must avoid the experimental factors in the design of various structures made of corrugated paperboard and to achieve the optimum design. However, due to the complexity of the corrugated paperboard itself in FEA for these structures, a significant time and effort are required for modeling and analysis. Therefore, it is necessary to analyze the simplified model and the equivalent mechanical properties for corrugated paperboard by flute type.
In this study, the equivalent mechanical properties for corrugated paperboard were analyzed by measuring the reaction force generated by applying a forced displacement to FE model for a unit cell of a corrugated paperboard with the same shape repeatedly. In addition, the analyzed equivalent mechanical properties were applied to a simplified model of the corrugated paperboard to analyze the S-S behavior through FEA. The results were compared with experiments.
The S-S behaviors of the corrugated paperboard analyzed by the FEA and experimental methods showed a significant difference above a certain amount of strain. The S-S behavior of the corrugated paperboard experimentally showed elastic-plastic behavior, although the equivalent mechanical properties applied to the simplified model of the corrugated paperboard were elastic based on the theoretical approach, so that the FEA results showed linearity. Therefore, when applying the calculated equivalent mechanical properties of the corrugated paperboard through FEA, consideration should be made about the characteristics of the section where the strain only increases without the increase of the load due to the flute of the corrugated paperboard. If the improved equivalent mechanical properties are applied to the FEA, it can be concluded that the FEA results close to the compressive behavior of corrugated paperboard can be obtained.
Although the nonlinearity was neglected in this study, the linearized and homogenized three-dimensional orthotropic equivalent mechanical properties be obtained, so that FEA in consideration of the effect of a certain section and the characteristics of the corrugated paperboard be easily performed.