Elastic properties identification of a bio-based material in tertiary packaging: Tools and methods development

Abstract

This study focuses on the use of bio-based materials for structural purposes in the packaging field, which requires the identification of their mechanical properties at a representative scale. The mechanical properties of bio-based materials are more variable than those of traditional composite materials. In a standard characterization approach using elemental coupons under uniaxial loading, the variability depends on the chosen representative elementary volume (REV), free edges, boundary conditions, etc. for elastic properties that are not identified for representative working conditions; this could lead to the ineligibility of these bio-based materials as structural materials. This paper contributes to the debate on how to study the response of bio-based materials within a structure, here a packaging structure as a logistic unit (LU) subjected to a compressive load simulating storage and stacking conditions. In the set of tools and methods for the design of packaging materials made of bio-based materials, an elastic nonlinear geometric finite element model (FEM) and an experimental approach are presented. The FEM allows the numerical identification of zones of interest within the LU. Inevitably, the FEM classically requires input data which are elastic properties of the equivalent homogeneous material. The design of the FEM is based on a calculation-test approach using an existing reference LU and it can be summarized in two main steps. The first step concerns the development of a FEM able to restore the experimental conditions of vertical compression imposed by transport standards for packaging. The second step is based on updating the input properties of the FEM by reverse identification, to achieve the representative working condition properties, using experimental results obtained on the existing reference LU. For the reverse identification a multi-scale investigation is mandatory. For this purpose, the linear elastic part of the load/vertical displacement curves (at the LU stiffness scale) and the displacement and strain fields measured (at the local LU scale) by 3D digital image correlation (3D DIC) are evaluated. Then, FEM property updating is carried out by reducing the deviation of displacement/strain fields between FEM and experimentally measured results (3D DIC). Finally, we explain how FEM and 3D DIC help in decision-making by allowing the recognition of zones of interest in a phase of design of new LUs with the concept of Multi-Instrumented Technological Evaluator (MITE).

Keywords

Mechanical design bio-based packaging multi-scale analysis compressive loading reverse identification FEM updating 3D DIC

Introduction

Types of packaging

One of the main roles of packaging is to protect the product from mechanical loadings that may occur during distribution and transportation. Three types of industrial packaging are used¹: (i) primary packaging, which is the container that is in direct contact with the product, (ii) secondary packaging that is used to cover the primary packaging in order to group a number of items together to make a sales unit and (iii) tertiary packaging, which is used to unitize large quantities of products or packages and transport them safely. Industrially this last packaging type is called unitized load or a logistic unit (LU). In this study we will focus on tertiary packaging. According to Jukes,² the level of content protection required is determined by the rigidity of the packaging box or the unit used, as well as the severity of the risks to which it may be exposed during transportation.

Interest of bio-based materials for packaging

In an era where ecological impact is one of the major preoccupations in logistics, many companies are beginning to incorporate bio-based materials into their packaging strategy to counter some of the environmental impacts generated by traditional plastic packaging,^3–5 and reduce the use of non-renewable raw materials, particularly petroleum.⁶ This reduced dependence on hydrocarbons is one of the long-term economic benefits promoting the use of bio-sourced packaging, as it is derived in most cases from abundant renewable resources such as wood, bamboo or flax.⁷ These materials have mechanical, thermal and functional properties that make them candidates for bio-based packaging applications.^4,5,8–14 To date, bio-based packaging has been widely applied in the food and cosmetic field as primary or secondary packaging.^15–19

Bio-based material mechanical properties for structural applications - case of tertiary packaging

No published study was found on the use of bio-based materials for the development of more structural components such as tertiary packaging or logistic units (LUs). For these units, the knowledge of the mechanical characteristics of the materials as a function of their manufacturing process is essential for the dimensioning.

It is well known that composite materials exist as structures. For example, in aeronautics, composite materials for structural applications are developed using a test pyramid approach.²⁰ This approach involves five levels of mechanical tests at different scales. The first level consists of characterization of numerous standardized test coupons to determine key mechanical properties such as the Young’s Modulus and ultimate strength and failure modes. The second level concerns the evaluation of the load capacity of structural basic elements. The third level is defined as the assembly of two or more basic elements to create structural details. The fourth level is then devoted to gathering several details to assemble a component. The final level considers testing of the whole component or structure. As can be inferred from the description, the test pyramid approach is long and costly, thus, becoming economically prohibitive in the field of packaging that is extremely cost sensitive.

In this paper, the aim is to focus on the response of the material within the structure. The development of tools and methods adapted to the design of structures using composite materials with high variability^21,22 such as bio-composites is proposed. The development phase is based on an existing logistics unit, which is called the reference LU, presented in §2. One of the tools is the finite element model (FEM) of the reference LU. To develop the FEM presented in §4, it is proposed to study the response of the bio-based material employed in the reference LU. This study is based on a dialogue between FEM (geometrically nonlinear elastic model) and experimental tests using Digital image correlation for the case of compression loading presented in §3. An approach to update mechanical property values of a REV (from elementary coupons) by reverse identification^23–26 is performed until a reasonable correlation is obtained between strain fields measured on the reference logistic unit and those calculated from the FEM presented in §5. Then, the FEM with updated input data is used to enable the numerical localization of zones of interest within the logistic unit under the compressive loading presented in §6. Based on experimental findings, the zones of interest are those where the behavior of the material is assumed to be elastic without damage, according to FEM. These zones of interest have areas and boundary conditions that are more representative of the material within the structure and its working conditions compared to the REV of an elementary coupon under uniaxial loading far away from the working conditions.

To study LUs experimentally, Digital Image Correlation (DIC) is used. Indeed, the DIC technique has been successfully applied to monitor in situ displacement fields of homogeneous^27–31 and heterogenous^23,32–35 materials under mechanical loadings. DIC has also been used to validate FE models that simulate compressive loadings on cardboard packaging.^36–38 Fadiji et al.³⁶ used 3D DIC to visualize the full-field displacement of corrugated paperboard packaging subjected to compression loading. Also, strain field distributions were derived from displacements, where the technique proved to be efficient in providing reliable information on the displacement and strain fields at the surface panels of this type of package. Viguié et al.³⁷ used DIC to analyze the out-of-plane displacement behavior of cardboard boxes during their compression. This technique provided relevant information about outer liner displacement and strain fields; therefore, the stress field of the panel could also be estimated. Furthermore, the displacement field was smoothed because it couldn’t be used in its raw form to calculate the strain field given the heterogeneous behavior of the material. However, these latter studies were limited to studying the mechanical behavior of the packages with small dimensions belonging to the secondary packaging category. To date, there are no studies that develop a reverse properties identification approach based on the DIC method on tertiary packaging to identify material properties considering not only its geometrical aspect but also its associated boundary conditions.

This paper presents a set of original tools and methods for the identification of mechanical properties of bio-based materials for structural applications. The proposed work focuses on the mechanical characterization during the compressive loading of the packaging. This choice depends on the loading conditions during storage and stacking during transport. A FEM of a compression test on a reference LU, simulating stacking conditions, is developed. Vertical compression tests on the same unit, instrumented by 3D DIC, are carried out to allow a calculation-test dialogue for material mechanical property updating by reverse identification to reach a reasonable correlation.

Bio-based material of the reference logistic unit

A unit qualified for transporting goods, the plywood box RAJA CB1286,³⁹ was chosen as “reference LU” cf. Figure 1. It is a container with the following dimensions: 1200 mm length, 800 mm width and 600 mm height. The unit consists of a set of 4 birch plywood lateral panels 6 mm thick. The plywood panel is made of 5 birch plies with a stacking sequence of [90/0/90/0/90], where 90° corresponds to a ply with its fibers in the length direction. The stacking sequence provided by the supplier is not optimal in terms of fiber placement in the panels because the fibers are not predominantly in the direction of stress flow (i.e. vertical compression being height direction). To manufacture these panels, the plies are assembled by gluing and hot pressing. Then they are linked together by angle irons which allow storage or transport of the LU in a space-saving manner by folding them up. The bottom and the cover are also made with the same plywood material and process.

Figure 1.

Reference Logistic Unit³⁹ (on left) and plywood panel (on right) with the detail of a homogenized REV.

Table 1 presents mechanical properties of the bio-based material used for the manufacturing of the birch plywood reference LU. Table 2 summarizes the dimensions of the reference LU. Both the in-plane and out-of-plane properties are taken into account in the numerical calculations. The hypothesis of the classical laminated plate theory is considered because of the ratio thickness-to-side ratio: from 1/200 to 1/100, with the use of first order shear deformation theory (FSDT), which considers the effect of out-of-plane properties. The properties of the homogeneous material are presented in Table 1. Two sources of material property values are provided, one from literature^40–42 and the other one from characterization tests performed during the present study. The order of magnitude of the values obtained from these tensile tests is the same as those reported in the literature. According to our tests, the behavior is linear elastic up to strain values of about 4 × 10⁻³. The mechanical response of the laminated material is that of a REV studied using elementary coupons, whose boundary conditions are those imposed by the tensile set-up configuration. In this article, our interest lies in the possibility of conducting a calculation-test dialogue under conditions different from the specific response of the laminated material within an elementary coupon, but rather with the proper conditions of plywood panels loaded in a LU subjected to compressive loading. Elastic properties are assumed to be the same in tension and compression.⁴²

Table 1.

Mechanical properties of the homogenized material (five ply birch plywood) from elementary coupons (1, 2 and 3 directions defined in Figure 1).

Five ply birch plywood in-plane: Homogenized mechanical and physical properties	Literature^40–42	Characterization tests (DIN 52377)
Total thickness: 5 ply thickness (mm)	6.5	6.1 ± 0.3
Volume mass (kg/m³)	690	-
Tensile strength long. (1) and transverse (2) (MPa)	110	87 ± 19
Compressive strength long. (1) and transverse (2) (MPa)	52	-
In-plane shear strength (MPa)	12	10.5 ± 1
E₁ tensile longitudinal elastic modulus (MPa)	9844	11942 ± 396
E₂ tensile transverse elastic modulus (MPa)	7656	7246 ± 233
ν₁₂ tensile Poisson coefficient	0.034	-
G₁₂ in-plane shear modulus (MPa)	610	602 ± 94
G₁₃ out-of-plane shear modulus (MPa)	206	-
G₂₃ out-of-plane shear modulus (MPa)	186	-

Table 2.

Dimensions of the reference logistic unit.

Parts	Length (mm)	Width (mm)	Height (mm)
Crate	1200	800	600
Pallet	780	95	144

Experimental set-up for LU in compression and measurement method

Experimental set-up

As already discussed above, our objective is to put the bio-based material into representative working conditions within the structure of a chosen LU. To do this, we propose to recreate the compression loading applied on the LU produced under storage and stacking conditions. Compression tests were performed in compliance with the ISTA 3E transport standard⁴³ by means of a compression machine dedicated to the qualification of packaging containers; the Zwick Roell Z100 TEW compression machine with a 100 kN load cell, available on the ILIPACK platform⁴⁴ cf. Figure 2. To obtain representative conditions, a pallet is placed on the top of studied LU, cf. Figure 2. The LUs are preconditioned at the ambient temperature and humidity of the laboratory for 24 h before the test. Five tests were performed to evaluate the standard deviation of the plywood material and the structural behavior. A 12.7 mm/min (half-inch per minute) compression displacement rate of the was applied on the LU loaded up to 800 daN (approximately to the force applied by an 800 kg container placed on top of the LU) to stay within the material elastic domain.

Figure 2.

Vertical compression test representative of stacking and storage conditions with 3D DIC system targeting the large LU panel.

LUs show significant out-of-plane displacements at the center of the large LU panels compared to the LU small panels. The displacement field was found to be symmetrical for all large LU panels. For this reason, it was chosen to use only the displacement field from half of a large LU panel to perform a more detailed analysis. Among the zones associated with such displacement fields, the zones of interest were chosen. As mentioned before, these zones of interest have an extended area and boundary conditions that are more representative of the plywood material within the LU than those of the REV elementary coupon under uniaxial loading.

Determination of 3D digital image correlation parameters

The tests were monitored using the digital image correlation method (3D DIC) to identify the displacement fields on half of a large LU panel during the compression loading tests. Image processing using any image correlation software depends on several parameters to run in an optimal way and give acceptable results; these include the subset size, the step, and the filter size.⁴⁵ The subset size is the area of interest that is used to track the displacement between selected images. To obtain a credible correlation analysis in DIC, the size of a subset must be large enough for the intensity pattern contained in the subset to be sufficiently distinctive to differentiate it from other subsets, but not too large to prevent the loss of more localized resolution. The step size is a parameter that controls the density of points (subset centers) at which DIC data is computed, and it influences the spatial resolution of the measurements. Typically, a step size of one quarter to one third of the subset size is recommended,^45,46 so that neighboring subsets partially overlap, and also to avoid having repetitive data and long computational time. A step of 2/7 subset size is used in this work.

Since the displacement field obtained by the DIC method is noisy, its numerical derivation would lead to even noisier strain maps. In the VIC^® software (VIC-2D^®/VIC-3D^® from Correlated Solutions), for example, the strain provided to the user corresponds to an average of the strains calculated in a close vicinity of the point considered. The user specifies, in addition to the size of the DIC area of interest (subset) and the step, the value of a parameter called “filter”, that determines the number of points to be considered in the proximity to evaluate the strain.⁴⁵ In practice, the number of neighboring points used is mandatory. The size of the virtual gauge for strain calculation depends on this filter number and the step size. Thus, a parametric study of the VIC 3D^® software is essential in order to optimize these two parameters during the calculation of the displacement and strain fields of the structure.

Figure 3 shows the effect of subset size on the displacement error for the unit compression test for a particular lighting configuration. For small subsets, the error varies inversely with subset size. However, for larger subsets, the error approaches an asymptote value. Similar trends have been reported in literature studies.^47,48 The displacement errors obtained in the present study range from less than 5 × 10⁻³ to 2.2 × 10⁻² pixels. The subset size that offers the best compromise between spatial resolution and displacement accuracy lies just below the asymptote; larger subset sizes do not further reduce the displacement error, but they lower spatial local resolution. The optimum range of subset sizes that should be used in image analysis lies between 50 and 100 pixels for approximately 5 × 10⁻³ pixels of displacement uncertainty.

Figure 3.

3D DIC measurement error on displacements for different subset sizes, location of the chosen value (79 pixels).

Random strain errors are reduced with increasing subset size, due to greater confidence in image correlation when more details are available in larger subsets (smoothing effect), resulting in lower standard deviations (SD) cf. Figure 4, as reported previously.^45,49 The SD reflects the deviation of measured displacements from their mean value. Figure 4 shows the trend of strain standard deviation asymptotes towards a constant value, as seen with displacement errors from a range of 60 to 90 pixels. Based on these two graphs a subset size of 79 pixels was chosen for the image processing. Similarly for the size of the filter, several filter sizes were varied to examine their effect on the strain SD, cf. Figure 5. A 15 points filter was chosen to define the virtual strain gauge size for the calculation of the strain fields. The chosen parameters for displacement measuring are given in Table 3.

Figure 4.

Strain standard deviation variation at 780 daN of load as a function of subset size, location of the chosen value (79 pixels).

Figure 5.

Strain standard deviation variation at 780 daN of load as a function of filter size, location of the chosen value (15 points).

Table 3.

Image processing parameters on VIC-3D^®.

Correlation criterion	ZNSSD
Interpolation method	Optimized 8-tap
Subset weights	Gaussian
Subset size	79 pixels
Step size	23 pixels
Strain tensor	Lagrange
Filter size	15 points

Figures 6 and 7 show, respectively, the distribution field of the out-of-plane displacement U₃ and the transverse strain E₂₂ of the half of large LU panel face before and after changing the VIC-3D^® software parameters. A minor change can be seen in the distribution quality of the out-of-plane displacements, noting that the response is more averaged in the case of the 79 pixels subset. However, the strain field’s distribution is more homogenized in the case of the 15 points filter for a suitable comparison with the FEM, which is fed by homogeneous orthotropic material properties.

Figure 6.

Out-of-plane displacement U₃ field’s distribution with, (a) 39 pixels subset and (b) 79 pixels subset.

Figure 7.

Transverse strain E₂₂ field’s distribution along y direction for 79 pixels subset with, (a) 3-point filter and (b) 15-point filter.

Finite element model of the reference LU under compression

A geometrically nonlinear elastic Finite Element Model (FEM) was developed using ABAQUS^® software to simulate the compression loading on the studied LUs. The representation was made within the framework of Continuum Mechanics, using the elastic properties of an equivalent homogeneous material. The starting point of the calculation-test dialogue is the values of the elastic properties obtained from birch test coupons,^40–42 mentioned in Table 1. These calculation assumptions are obviously far removed from the reality of the bio-based material in the panel because the local bio-based material heterogeneity is not accessible^21,22 and our main objective is to develop operational tools for LU design.

The different parts of the LU were modeled using quadrilateral linear shell elements, except for the pallets which were modeled using deformable volume elements to facilitate the constraint description (i.e. load transfer with the cover and the bottom by surface/surface general contact and stiffness pallet representation), cf. Figure 8. The moving compression plate was modeled using rigid shell elements. The material assigned for the lateral plywood LU panels, bottom and cover is birch plywood with a thickness of 6 mm, with property values homogenized from five plies, cf. Table 1. The material used for the pallet was poplar (wood species), considered to be as an isotropic material to lighten the model.

Figure 8.

Finite element model of vertical compression test on the reference logistic unit (reference LU).

Assumptions

are proposed to simplify the calculation and the determination of the mechanical properties of the studied structure:

• The plywood is modeled as a one ply composite with homogenized orthotropic properties using shell elements,

• The hooks which are used to ensure the connection between the lateral panels, the bottom and the cover panels are modeled by partitions on the faces of the panel with Tie interactions (ABAQUS^® constraint function) placed between these partitions. These partition dimensions are 50 × 30 mm². Tie interactions in this case do not affect the stiffness of the LU. Also, this allows the transfer of the moment from the cover to the lateral panels (because the cover leans on the lateral panels), otherwise it would not be possible, since the shell elements do not have mesh elements through the thickness (which is the case of lateral panels).

• The angled metal borders on the cover and bottom are not modeled here.

Four FEMs are studied to evaluate the computation time and the mesh configuration influence on the LU stiffness, maximal strains (Ε₁₁ and Ε₂₂) and out-of-plane displacement (U₃). The models were:

• Model I with 5 mm mesh size on the entire geometry of the logistic unit,

• Model II with 8 mm mesh size on the entire geometry of the logistic unit,

• Model III with 25 mm mesh size in the central zones of the panels and with 5 mm mesh size near singular contact zones (i.e. panel/cover, panel/angle iron),

• Model IV with 25 mm mesh size in the central zones of the panels and with 10 mm mesh size near singular contact zones (i.e. panel/cover, panel/angle iron).

To evaluate the influence of these mesh types and sizes, the LU stiffness, the maximal strain and out-of-plane displacement values are determined and compared for each type of mesh size, cf. Table 4 and Figures 9–10.

Figure 10(a) shows that the maximum value of the LU stiffness is obtained using Model I with the finest mesh elements (i.e. 5 mm). However, the computation time of Model I is nearly 330 min with the workstation used here. It would be interesting to reduce the computation time to facilitate the calculation-test correlation. Figure 10(a) shows that Model III has the closest LU stiffness with the shortest computation time (30 min vs 330 min) with respect to the finest mesh Model I, with a difference of 5 % in LU stiffness.

Almost the same results are found on relative deviation versus number of elements as a function of maximal strains (Ε₁₁ and Ε₂₂) and out-of-plane displacement (U₃) at 500 daN load, cf. Figure 10(b). Note that Model I is considered as the reference model and the deviation is calculated with respect to results from this model. Figure 10(b) shows that the relative deviations of these parameters using Model III are small and converge towards a null deviation value.

It can be concluded that Model III provides satisfactory results with a good compromise between accuracy and computation time ratio (studied duration/shorter duration), cf. Table 4 (11 times faster than Model I).

Table 4.

Computation time required for each mesh size for shell elements (in bold for the chosen mesh size).

FEM	Mesh size (mm) common panel zones	Mesh size (mm) contact zones	Number of elements	Computation time ratio
Model I	5	5	28320	11
Model II	8	8	11172	1.5
Model III	25	5	2375	1
Model IV	25	10	1627	1

The [bold] was used to show the optimal and the chosen model (Model III).

Figure 9.

Location of maximum value zone at 500 daN (on Model III) with, (a) out-of-displacement U_3, (b) transverse strain E₂₂ and (c) Longitudinal strain E₁₁.

Figure 10.

Mesh size configuration at 500 daN load with respect to; (a) stiffness variation (top) and (b) relative deviation (bottom), versus the chosen mesh size.

Mechanical properties reverse identification

Different phases are required for reverse identification with a multi-scale analysis.

Baseline data

Figure 11 illustrates the five experimental and one numerical “load/vertical displacement” responses during the compression tests. The five experimental “load/vertical displacement” curves show repeatability. Table 5 indicates the LU stiffness, the maximum strains and out-of-plane displacement (U₃) and standard deviation for the five tests. Only the linear elastic part of the Load-Vertical displacement curve is studied, which is within the normal operating limits of the LU. “FEM 1” corresponds to the first numerical model with initial property values. The difference between the experimental and numerical LU stiffnesses (the slope of the load vs displacement curve) is 115 %. To reduce this discrepancy, a reverse identification is proposed in the framework of a multiscale analysis.

Figure 11.

Load/vertical displacement curve of the experimental (cf. five dotted lines) and numerical (cf. One continuous line of FEM1 before properties updating by reverse identification) compression tests.

Table 5.

Experimental compression test results at different scales at 780 daN load.

Studied parameters	Average	Standard deviation
LU experimental stiffness (daN/mm)	581	14 (2.4 %)
Out plan displacement U₃ _max (mm)	7,5	0,1 (1.33 %)
Longitudinal strain E_{11 max}	2 × 10⁻⁴	3 × 10⁻⁵ (15 %)
Transverse (vertical) strain E₂₂ _max	6 × 10⁻⁴	1,5 × 10⁻⁵ (2.5 %)

First phase of the reverse identification: LU stiffness

An adjustment of the numerical response of the LU with respect to the experimental one is a proposed way to identify the mechanical properties of the studied structure. This was achieved by adjusting the numerical stiffness with respect to the experimental stiffness by varying the material’s longitudinal modulus of elasticity E₁ and transverse modulus of elasticity E₂. The ratio E₂/E₁ = 0.77 is kept constant, thus it is only necessary to vary E₁. The remaining elastic properties were also kept constant. Several models with different property values were computed until “FEM 2” where the updated LU stiffness aligns with the experimental one, cf. Figure 12. The elasticity moduli of the studied structure were found to be 3000 MPa and 2330 MPa for the longitudinal and transverse directions, respectively.

Figure 12.

Load/vertical displacement curves of the experimental (5 dotted lines) and numerical (2 continuous lines of FEM1 and FEM2) compression tests.

Moreover, FEM 2, fed by the new property values, can also be used to compare the numerical out-of-plane displacement and strain fields with those obtained by the 3D DIC. Figure 13 shows a comparison between the numerical and experimental out-of-plane displacements for three load stages. For each load stage, a general similarity between the numerical and experimental distributions of out-of-plane displacement fields within the panel face can be seen. However, there is a significant difference in the order of magnitude of displacements, nearly 60 % difference in the case of low loads (here 400 daN), up to 80 % in the case of a 780 daN load.

Figure 13.

Out-of-plane displacement U₃ fields during the compression test with, (a) experimental measurements and (b) numerical model respectively, at 400 daN (top), 600 daN (middle) and 780 daN (bottom).

Transverse strain Ε₂₂ fields are also compared, cf. Figure 14, and there is little similarity between measured and predicted values. The results of DIC and FEM analysis are not comparable in the case of 400 and 600 daN loads, where an out-of-plane displacement U₃ begins to form. On the other hand, a slight similarity can be found in the case of a 780 daN load, when the panel takes on a regular out-of-plane displacement shape, but a significant difference between DIC measurements and FEM values is observed. This investigation helps to conclude that fitting the model considering only load/vertical displacement response is not sufficient.

Figure 14.

Transverse strain Ε₂₂ fields during the compression test with, (a) experimental measurements and (b) numerical model respectively, at 400 daN (top), 600 daN (middle) and 780 daN (bottom).

Following this first attempt, a second phase of adjustment of the numerical input data was performed, based on minimizing the difference between experimental and numerical displacement-strain fields.

Second phase of the reverse identification: using the 3D DIC results

The reverse identification minimizes the difference between the maximal values of both experimental and FEM results at a local scale (in a chosen zone). Indeed, these input data (mechanical properties) were previously updated by reverse identification using the linear part of the load/vertical displacement curves.

However, it should be noted that the methods of representing the results field with VIC 3D^® and ABAQUS^® are different. They do not use the same data coordinate system, nor the same point distribution. It is therefore necessary to match the positioning of the measurement points, i.e., to place them in the same coordinate system, to be able to compare the values of the experimental fields (analyzed by VIC 3D^®) and the numerical ones (obtained by ABAQUS^®) at the same points of the panel. To achieve this, an interpolation algorithm is developed using MATLAB^® software, which involves the following requirements:

• choosing an identical analysis zone of the panel to be studied at 780 daN loading to compare experimental and numerical fields, cf. Figure 15 (500 × 300 mm², dashed line),

• positioning the 3D DIC experimental points (in an irregular surface) in accordance with the perfect virtual plan used in the FEM,

• interpolating the experimental points with reference to the model point’s coordinate system to match the coordinates of the points in the same reference.

Figure 15.

Selected area (dashed line) for the second phase of property adjustment with respect to out-of-plane displacement U₃ fields with respectively, (a) experimental measurements and (b) numerical model.

Figure 16 shows the out-of-plane displacement U₃ fields resulting from this interpolation. The similarity of the profiles is observed. However, the values of the experimental and numerical data are not so close, for the same values of module E₁ = 3000 MPa and E₂ = 2330 MPa (reverse identification first phase).

Figure 16.

Out-of-plane displacement U₃ field at 780 daN with, (a) experimental measurements and (b) numerical model before second phase (i.e. first phase) of property adjustment with respect to DIC field’s distribution on MATLAB^®.

The second phase of property updating of the numerical model was carried out. To do this, the two elastic moduli of the material are updated until a match is obtained between the “FEM 3” numerical and experimental out-of-plane displacement U₃ field values, cf. Figure 16. The identified longitudinal and transverse elastic modulus values are 4540 MPa for E₁ (vs 3000 MPa) and 3500 MPa for E₂ (vs 2330 MPa). Figures 17 and 18 display the updated numerical and experimental displacement U₃ and strain Ε₂₂ fields. An acceptable compromise is found between the numerical and experimental fields and values which validate the approach. Thus, these identified properties provide a description of the material’s response in the structure and under this type of loading, see Table 6.

Figure 17.

Out-of-plane displacement U₃ field at 780 daN with, (a) experimental measurements and (b) numerical model after the second phase of property adjustment with respect to DIC field’s distribution on MATLAB^®.

Figure 18.

Transverse strain E₂₂ field at 780 daN with, (a) experimental measurements and (b) numerical model after the second phase of property adjustment with respect to DIC field’s distribution on MATLAB^®.

Table 6.

Identified mechanical properties of the birch plywood panel in the LU subjected to compression loading.

	E₁ (MPa)	E₂ (MPa)	ν₁₂	G₁₂ (MPa)
Elementary coupon properties (i.e. before reverse identification)	9844	7656	0.034	610
Final properties by reverse identification with FEM3	4540	3500	0.034	610

The load/vertical displacement curve of the new model “FEM 3” is plotted, cf. Figure 19, showing where the panel stiffness of the final model is positioned compared to experimental panel stiffness. A 34 % difference between the “FEM 3” stiffness and the experimental data is noted (vs 115% for “FEM 1”). This difference can be explained because the finite element model does not account for the variability of boundary conditions associated with experimental tests, such as the influence of pallet sole geometric defects. These can affect the LU loading conditions and therefore the response of the panel. According to Table 6, it is also clear that the updated property values differ significantly from the property values provided by elementary coupons, by 54%. Note that a FEM, built with input data of tensile mechanical properties derived from elementary tensile tests, leads to an unsatisfactory calculation-test correlation. A satisfactory calculation-test correlation is achieved on the scale of the structure by updating values of the elementary coupon’s properties.

Figure 19.

Load-displacement curve of the experimental (5 dotted lines) and numerical (3 continuous lines) compression tests after properties updating with respect to DIC field’s distribution on the LU scale.

The flow chart on Figure 20 shows the different steps of the study leading to the zone of interest identification and the MITE^50,51 design. It is the subject of the following paragraph (cf. §6).

Figure 20.

Flow chart of the different steps of the study leading to the zone of interest identification and the MITE design.

Zone of interest identification and multi-instrumented technological evaluators

Zone of interest identification

As mentioned previously, the validation of the defined FEM allows us to identify the different areas of interest within the logistic unit. The Tsai-Hill criterion is used to identify these areas as the “most stressed zones” within the unit during compression, cf. Figure 21. This criterion qualitatively translates the risk of a probable appearance of damage using material strength properties. Figure 21 illustrates areas (zones of interest) not only located at the center of the panel (Z.1) but also in the vicinity of the pallet soles (Z.2).

Figure 21.

Tsai-Hill distribution on the large LU panel face at 780 daN of load with the location of the zones of interest Z.1 and Z.2.

As explained previously, the zone 1 (at the center of the panel) represents an out-of-plane elastic displacement U₃ field because of a vertical compressive stress, which is an expected response in the case of rectangular packaging.³⁶ In the zone 2, a crushing (excessive compression) phenomenon is observed experimentally at loads higher than 800 daN, cf. Figure 22; this is due to the use of top and bottom soles, as well as high shear stresses close to the contact with the iron angles.

Figure 22.

View of the experimental crushing of the LU in the vicinity of boundary conditions.

In a design approach, the logistic unit, that is the final structure, does not exist, since obviously it has yet to be designed. We can resort to the concept of the multi-instrumented technological evaluator (MITE) for performing property value updating or for studying specific zones of interest.

Added value of multi-instrumented technological evaluators

Zone 1 (Z.1) is the subject of a specific study using the concept of MITE, introduced briefly in this paper. This area allows representative properties identification of the studied structure, far from combined loading areas, where boundary conditions affect the compression response. The evaluators are parts with smaller dimensions than the final structure and whose boundary conditions and loading kinematics allow the restoration of a local state equivalent to that of the zone of interest (Z.1) of the LU as illustrated in Figure 23. The MITE approach aims to design new logistic units made from different architectures of bio-based materials with a large zone of interest such as materials with bamboo fiber reinforcement.⁵² The advantage of the MITE approach⁵⁰ is to use only the technological evaluators, manufactured in large numbers, before investing in a validation test on the final logistic unit. Figure 24 shows the framework of the design investigation with a technological evaluator mounted on a multi-axial experimental set-up⁵⁰ to provide a similar response of the zone of interest (Z.1) of the logistic unit at a smaller scale. The benefit of this approach is illustrated in Figure 25 in the case of the studied reference LU. A calculation-test dialogue (cf. Figure 25) allows several opportunities such as the identification of elastic mechanical properties (not presented in this paper) and essentially applying the similar type of loading, using the multi-axial testing machine. The testing machine can apply combined loads thanks to four electric actuators. These actuators are attached to a steel quadrangle, which transmits the load to the evaluator for a chosen kinematic, cf. Figure 24.

Figure 23.

Transverse strain E₂₂ field on Z.1 (continuous and dashed line) of the LU at 780 daN.

Figure 24.

Multi-axial testing machine with the LU MITE mounted, painted with a speckled pattern.

Figure 25.

Transverse strain E₂₂ distribution in the zone of interest (Z.1) after compression loading with, (a) an experimental response on LU MITE, (b) a FEM response on LU MITE and (c) a FEM response on the LU.

Conclusion

This paper presents a study of the elastic heterogeneous composite material properties of a structure, here a logistic unit (LU), based on bio-based materials, under normal stacking and storage conditions. The initial assumption based on properties obtained with elementary coupons characterization as input data resulted in an unsatisfactory FE calculation-experimental tests correlation. A difference of 115 % between the experimental and numerical LU stiffness values was noted. Therefore, we proposed a property updating approach, via a reverse identification in the framework of a multiscale approach, by varying the longitudinal and transverse elasticity modulus values E₁ and E₂, respectively.

In this paper, numerical (FEMs) and experimental (protocols and measurement techniques at different scales) tools are developed. The proposed methodology is based on an original calculation-test dialogue, to be capable of considering the representative working conditions of bio-based materials as structural materials. Using ABAQUS^® software, the proposed FEM with non-linear geometry elastic material hypotheses reproduces the experimental vertical compression conditions imposed by the ISTA 3E standard. In a first phase, the adjustment (updating by reverse identification on the structural scale) of the input data (mechanical properties) of the FEM is then performed using the linear zone of the load/vertical displacement responses (LU stiffness) of the tested units. An acceptable LU stiffness correlation was obtained but the results proved that this property updating on the LU scale is insufficient to provide a satisfactory correlation at the local scale for out-of-plane displacement and strain fields in the lateral LU panels. In a second phase, FEM property updating was carried out by reverse identification using the maximum values of the out-of-plane displacement/strain fields defined by the 3D DIC technique and FE calculation. This required a MATLAB^® algorithm to obtain the correspondence between 3D DIC and FEM points in the same coordinate system. A reasonable correlation between the numerical and experimental results was achieved at the global scale with a 34 % difference between the experimental and numerical LU stiffnesses and, at the local scale, between the different fields, both in terms of distribution and order of magnitude. This final updating led to a 53% reduction between the initial property values and the final ones. The proposed approach also allows other more sophisticated choices to update more elastic properties using optimization algorithms which are not shown here. The advantage of the method presented is that it shows that there is a need to adjust the values of the elastic properties, and that this adjustment must consider the response of the structure at different scales. The disadvantage of the method is to reveal that if we aim for a satisfactory correlation at the local scale in what has been defined as the zone of interest, then we cannot necessarily correlate global stiffnesses at the macroscopic scale. This is why the use of Multi Instrumented Technological Evaluator (MITE) is discussed at the end of the paper, as it enables a direct dialogue to be established between the loading (in this case, compression) and the zone of interest, without depending on the influence of panel boundary conditions within the logistics unit under study. The last part of the paper introduces the benefits of the MITE concept to design new logistic units, based on sub-scale parts and a multiaxial set-up, to be close to the LU work conditions.

Footnotes

Acknowledgments

The authors would like to thank the Occitanie Region (France) and the IUT Paul Sabatier (France) for funding Mohamed Hichem Saihi’s thesis. Our gratitude also goes to the Ilipack platform at IUT Paul Sabatier in Castres city (France), for their financial support in carrying out this work. The authors thank Nathalie Rocher, Laurent Crouzeix from ICA for their valuable support. We also would like to thank Cobratex French company at Carbonne, in the Occitanie Region, for the fruitful discussions about bio-based tertiary packaging. The authors would like to thank Dr Peter Davies from IFREMER Brest (France) for his precious help in the proofreading of this article.

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was supported by the Mixed syndicate, castres-mazamet; Institut Clément Ader, Université de Toulouse; Ilipack Platform.

ORCID iDs

Mohamed Hichem Saihi

Sonia Sassi

Francis Collombet

Redouane Zitoune

Data Availability Statement

Data sharing not applicable to this article as no datasets were generated or analyzed during the current study.

References

Saghir

. The concept of packaging logistics. In: Proceedings of the Fifteenth Annual POMS Conference. Mexico: Cancun, 2004, pp. 1–31.

Jukes

. Rigid boxes. Paper and Paperboard Packaging Technology. Amsterdam: Wiley Online Library, 2005, pp. 249–261.

Zhang

Zhao

. Green packaging management of logistics enterprises. Phys Procedia 2012; 24: 900–905.

Duhovic

Peterson

Jayaraman

. Natural-fibre–biodegradable polymer composites for packaging. Properties and performance of natural-fibre composites. Amsterdam: Elsevier, 2008, pp. 301–329.

Sahu

Gupta

. A review on the properties of natural fibres and its bio-composites: effect of alkali treatment. Proc Inst Mech Eng Part L 2020; 234: 198–217.

Ramesh

Palanikumar

Reddy

. Plant fibre-based bio-composites: sustainable and renewable green materials. Renew Sustain Energy Rev 2017; 79: 558–584.

Le Duigou

Davies

Baley

. Environmental impact analysis of the production of flax fibres to be used as composite material reinforcement. J Biobased Mater Bioenergy 2011; 5: 153–165.

McCracken

Sadeghian

. Green sandwich composites fabricated from flax FRP facings and corrugated cardboard cores. In: 17th International Conference on Non-Conventional Materials and Technologies (XVII NOCMAT 2017), Merida, Yucatan, Mexico, November 26th – 29th.

Betts

Sadeghian

Fam

. Structural behaviour of sandwich panels constructed of foam cores and flax FRP facings. In: Canadian Society for Civil Engineering (CSCE) Annual Conference. Vancouver, Canada, 31 May – 03 June 2017.

10.

Bambach

. Durability of natural fibre epoxy composite structural columns: high cycle compression fatigue and moisture ingress. Composites Part C: Open Access 2020; 2: 100013.

11.

Bledzki

Gassan

. Composites reinforced with cellulose based fibres. Prog Polym Sci 1999; 24: 221–274.

12.

Faruk

Bledzki

Fink

H-P

, et al. Biocomposites reinforced with natural fibers: 2000–2010. Prog Polym Sci 2012; 37: 1552–1596.

13.

Bernstad Saraiva

Pacheco

Gomes

, et al. Comparative lifecycle assessment of mango packaging made from a polyethylene/natural fiber-composite and from cardboard material. J Clean Prod 2016; 139: 1168–1180.

14.

Sydow

Bieńczak

. The overview on the use of natural fibers reinforced composites for food packaging. J Nat Fibers 2018; 16: 1455621.

15.

Darwin

Prajna

Tamba

. Bamboo fiber-PLA composite materials for disposable food and beverages packaging tools: a brief review. IOP Conf Ser Earth Environ Sci 2021; 926: 012087.

16.

Pandit

Nadathur

Maiti

, et al. Functionality and properties of bio-based materials. In: Green and sustainable advanced packaging materials. Berlin: Springer, 2018, pp. 81–103.

17.

Petersen

Væggemose Nielsen

Bertelsen

, et al. Potential of biobased materials for food packaging. Trends Food Sci Technol 1999; 10: 52–68.

18.

Wang

Euring

Ostendorf

, et al. Biobased materials for food packaging. Journal of Bioresources and Bioproducts 2022; 7: 1–13.

19.

Tawakkal

Cran

Bigger

. Effect of poly (lactic acid)/kenaf composites incorporated with thymol on the antimicrobial activity of processed meat. J Food Process Preserv 2017; 41: e13145.

20.

Rouchon

. Certification of large airplane composite structures. In: ICAS Congress Proceedings, Recent Progress and New Trends in Compliance Philosophy. Stockholm, Sweden: ICA, 1990, pp. 1439–1447.

21.

Davila

Crouzeix

Douchin

, et al. Identification and modelling of the in-plane reinforcement orientation variations in a CFRP laminate produced by manual lay-up. Appl Compos Mater 2018; 25: 647–660.

22.

Davila

Crouzeix

Douchin

, et al. Spatial evolution of the thickness variations over a CFRP laminated structure. Appl Compos Mater 2017; 24: 1201–1215.

23.

Petureau

Doumalin

Bremand

. Identification of local elastic parameters in heterogeneous materials using a parallelized FEMU method. Int J Appl Mech Eng 2019; 24: 140–156.

24.

Avril

Pierron

Sutton

, et al. Identification of elasto-visco-plastic parameters and characterization of Lüders behavior using digital image correlation and the virtual fields method. Mech Mater 2008; 40: 729–742.

25.

Collombet

Mulle

Grunevald

Y-H

, et al. Contribution of embedded optical fiber with Bragg Grating in composite structures for tests-simulations dialogue. Mech Adv Mater Struct 2006; 13: 429–439.

26.

Passieux

J-C

Bugarin

David

, et al. Multiscale displacement field measurement using digital image correlation: application to the identification of elastic properties. Exp Mech 2015; 55: 121–137.

27.

Bermudo Gamboa

Martín-Béjar

Trujillo Vilches

, et al. 2D–3D digital image correlation comparative analysis for indentation process. Materials 2019; 12: 4156.

28.

Candau

Pradille

Bouvard

J-L

, et al. On the use of a four-cameras stereovision system to characterize large 3D deformation in elastomers. Polym Test 2016; 56: 314–320.

29.

Száva

Pál-Gálfi

Ambrus

, et al. Video image correlation-based method used for the study of the torsional vibrations of an adder gearbox. Energies 2021; 14: 6129.

30.

Ramos

Braga

Eslami

, et al. Comparison between finite element method simulation, digital image correlation and strain gauges measurements in a 3-point bending flexural test. Procedia Eng 2015; 114: 232–239.

31.

Tang

Z-Z

Liang

Xiao

Z-Z

. Three-dimensional digital image correlation system for deformation measurement in experimental mechanics. Opt Eng 2010; 49: 103601.

32.

Crouzeix

Périé

Collombet

, et al. An orthotropic variant of the equilibrium gap method applied to the analysis of a biaxial test on a composite material. Compos Appl Sci Manuf 2009; 40: 1732–1740.

33.

Torres Arellano

Crouzeix

Douchin

, et al. Strain field measurement of filament-wound composites at ±55° using digital image correlation: an approach for unit cells employing flat specimens. Compos Struct 2010; 92: 2457–2464.

34.

Mehdikhani

Aravand

Sabuncuoglu

, et al. Full-field strain measurements at the micro-scale in fiber-reinforced composites using digital image correlation. Compos Struct 2016; 140: 192–201.

35.

Zitoune

Crouzeix

Collombet

, et al. Behaviour of composite plates with drilled and moulded hole under tensile load. Compos Struct 2011; 93: 2384–2391.

36.

Fadiji

Coetzee

Opara

. Evaluating the displacement field of paperboard packages subjected to compression loading using digital image correlation (DIC). Food Bioprod Process 2020; 123: 60–71.

37.

Viguié

Dumont

Orgéas

, et al. Surface stress and strain fields on compressed panels of corrugated board boxes. An experimental analysis by using Digital Image Stereo-correlation. Compos Struct 2011; 93: 2861–2873.

38.

Kueh

Dahm

Emms

, et al. Digital image correlation analysis of vertical strain for corrugated fiberboard box panel in compression. Packag Technol Sci 2019; 32: 133–141.

39.

RAJA packaging. https://www.raja.fr/

40.

Wang

Crocetti

, et al. In-plane mechanical properties of birch plywood. Construct Build Mater 2022; 340: 127852.

41.

Mechanical properties of wood. https://www.atomer.fr/1/Proprietes-mecaniques-des-essences-de-bois.html

42.

Omer

Ćehić

Hodžić

. Analysis of properties of plywood for some areas of usage in wooden constructions. In: International conference “Trends in the Development of Machinery and Associated Technology, Hammamet, Tunisia, 05-09 September 2007.

43.

ISTA 3E 2017 . Similar packaged-products in unitized loads for truckload shipment.

44.

Packaging transport qualification platform ILIPACK (Castres, France). https://ilipack.fr/

45.

International Digital Image Correlation Society . In: Jones

EMC

Iadicola

(eds) A Good Practices. Guide for Digital Image Correlation, 2018, ed1.

46.

Wigger

Lupton

Tong

. A parametric study of DIC measurement uncertainties on cracked metals. Strain 2018; 54: e12291.

47.

Rajan

Rossol

Zok

. Optimization of digital image correlation for high-resolution strain mapping of ceramic composites. Exp Mech 2012; 52: 1407–1421.

48.

Robert

Nazaret

Cutard

, et al. Use of 3-D digital image correlation to characterize the mechanical behavior of a fiber reinforced refractory castable. Exp Mech 2007; 47: 761–773.

49.

Pan

Xie

Wang

, et al. Study on subset size selection in digital image correlation for speckle patterns. Opt Express 2008; 16: 7037–7048.

50.

Collombet

Grunevald

Y-H

Crouzeix

, et al. Repairing composites. In: Advances in Composites Manufacturing and Process Design: Amsterdam. Elsevier, 2015, pp. 197–227.

51.

Torres

Collombet

Douchin

, et al. Numerical and experimental value added of multi-instrumented technological evaluator for the analysis of thick monolithic composite structures with singularity details. Compos Struct 2015; 127: 41–50.

52.

Robin

Davies

Arhant

, et al. Mechanical performance of sandwich materials with reduced environmental impact for marine structures. Jnl of Sandwich Structures Materials 2024; 26: 99–113.