MATHEMATICAL MODEL AND NUMERICAL SIMULATION OF OXYGEN TRANSPORT IN MULTILAYER FOOD PACKAGING: A SEMI-ANALYTICAL APPROACH.

: The purpose of this work was to develop a numerical simulation methodology using a free software capable of predicting the behavior of oxygen in multilayer polymeric systems to optimize packaging configurations with a greater barrier to O2. The Laplace Transform method was used to solve partial differential equations, with numerical inversion by the Stehfest algorithm together with the Gauss-Sidel method. It was possible to simulate the oxygen concentration profile in packaging systems containing two or three materials in the configuration. For some systems, the increase in the simulation time caused a numerical fluctuation at the material interface, especially if oxygen diffused from a material with greater diffusivity to one with lower diffusivity. The numerically simulated oxygen concentration profile can be associated with results from the literature, experimentally validated. The results showed to be promising for application in diffusion problems of multilayer systems, and the methodology can be used to describe systems with more layers. Mathematical modeling and numerical simulation have been allies to predict the behavior of oxygen in different materials and layer configurations, and it is relevant to provide a better understanding of the oxygen transport fundamentals involved and shorten product development cycle time and cost, contributing to the economy and new product development.


INTRODUCTION
The main function of packaging is to preserve the quality of food, promoting physical protection (against shock and light), protection against biological agents, or chemical protection, such as the transference of small molecules of water vapor, gases, and aromatic compounds (RESTUCCIA et al., 2015;WANI;SINGH;LANGOWSKI, 2014), in addition to information and convenience in handling, ensuring a cost-benefit ratio, and all aspects of food security (ANUKIRUTHIKA et al., 2020). All of these barrier properties are difficult to obtain from just one material, which is why multilayer packaging has been developed in recent decades so that each passive layer provides the properties of each polymer used in the packaging structure (NORIEGA; ESTRADA; LÓPEZ, 2014). Multilayer systems can also contain substances or additives that have some interaction with the food or the atmosphere that permeates it, capturing or releasing molecules of interest, making this system active (DURY-BRUN et al., 2007). These solutions extend shelf life and improve food quality and safety (HAN et al., 2018). The multilayer packaging on the market consists of 3 to 12 layers, and in the case of food packaging, the most common is from 3 to 7 layers (ANUKIRUTHIKA et al., 2020).
The ability of the polymer that makes up the packaging to promote or limit gas and steam exchanges between food and the environment is an important object of study, since, for some specific foods, limiting this permeant reduces chemical and biological degradation, as well as microbiological growth, and promotes the preservation of organoleptic properties (NORIEGA; ESTRADA; LÓPEZ, 2014).
However, a polymer with a good oxygen barrier, for example, may have a low water barrier (CRIPPA; SYDENSTRICKER; AMICO, 2008). Thus, multilayer polymeric packaging offers additional functionality, combining functional properties, as high barriers against gases and water vapor, even with reduced thickness, making it possible to produce thinner, lighter, and more compact packaging (ANUKIRUTHIKA et al., 2020;CRIPPA et al., 2007). Most techniques for the development of packaging have been applied in recent decades, due to the significant role played by oxygen in the deterioration of food, such as lipid oxidation, enzymatic activity, degradation of vitamins, and growth of deteriorating microorganisms and pathogens ( VAN BREE et al., 2010).
The mechanism of gas diffusion in polymeric materials was described by Thomas Graham, who proposed that the permeation process involved the dissolution of the gas with subsequent passage of the dissolved species through the membrane (MCKEEN, 2012), and occurs in three stages: (1) adsorption and solubilization of the permeant, that is, of the gas molecule on the surface; (2)  GOSWAMI; MAHAJAN, 2009;MORRIS, 2017). This mechanism is described by the laws of Henry and Fick, where the expression that relates the permeation rate with the area and thickness of the film can be obtained. Thus, it is understood that the concept of permeability is associated with the evaluation of the barrier properties of a plastic material (SIRACUSA, 2012), the properties of the penetrant, and the environmental conditions in which they are (LECHEVALIER, 2016).
To determine the best structures which have the required barrier capacity for certain foods and at the lowest possible cost, several tests of different materials and thicknesses are necessary. Mathematical modeling and numerical simulation have been great allies to predict the behavior of oxygen in different materials and layer configurations, also including studies with oxygen scavengers incorporated in these types of packaging. Such models are helpful to provide a better understanding of the oxygen transport fundamentals involved and shorten product development cycle time and cost, which is vital to the economy and new product development.
The methodologies used to simulate the behavior of oxygen in these works rely on analytical methods for solving the model's differential equations, as well as specific software for data application and simulation of results (CARRANZA; PAUL; BONNECAZE, 2010BONNECAZE, , 2012DI MAIO et al., 2017;FERRARI et al., 2009).
The Laplace Transform is an established technique for solving differential equations. In various circumstances, the transforms of a function and their inverses can be tabulated. However, in many practical situations, the use of tables of inverse transformations and the use of inversion formulas are not adequate, requiring a numerical inversion (DAVIES; CRANN, 1999).
The purpose of this work was to predict and compare the oxygen concentration that permeates a multilayer film with different materials and thicknesses, using a numerical simulation tool to select suitable materials and configurations for food packaging. The application of the developed tool, using the Laplace Transform combined with a Stehfest code for numerical inversion (STEHFEST, 1970) and the Gauss-Siedel method, made it possible to observe the oxygen concentration through the simulated layers using a free Fortran 77 compiler. Applying it to different polymeric materials with different oxygen barriers, it is possible to choose the best material for each type of food, since each one has a different need.  Fick's second Law describes diffusion in a non-stationary state, where the concentration gradient of the diffused substance is a function of time. In the absence of any chemical reactions and considering the flow in only one direction (x), we have: is the concentration of the diffused substance (oxygen), ( ) is the time in which the oxygen diffusion process occurs, ( ) is the coordinated onedimensional space measured in the polymeric material section, and ( 2 ) ⁄ is the diffusion coefficient of the gas molecule in the polymer.
The diffusion equations for each layer of the polymeric film in the present study are given by:

INITIAL AND BOUNDARY CONDITIONS
As an initial condition, the absence of oxygen in all layers was considered (Coxy (x,0) = 0). Given the structure of the multilayer film and the configuration of the arrangement, the symmetry of the oxygen concentration profile was established in the intermediate section of layer 2, while on the outer surfaces of the film the following condition was considered: where is the oxygen concentration that will permeate the material (A or C) and s is the oxygen solubility coefficient. Respecting the continuity of the mass flow at the interfaces between the layers, we have the following conditions:

PARAMETERS AND CONSTANTS
The parameters and constants used in the present study (table 1)   The thickness of the three layers of the polymeric film was fixed at 13µm, 9µm, 13µm for layers 1, 2 and 3, respectively, totaling 35µm, a configuration also used in a study (DI MAIO et al., 2017), in which they simulated applications for minimally processed fruits and vegetables.

NUMERICAL SOLUTION
In the equations for each layer of the polymeric film (2, 3, and 4), as well as in the initial and boundary conditions (5 to 9), the analytical Laplace Transform was applied. Therefore, the partial differential equations (PDEs) were transformed into ordinary differential equations (ODEs), equations 10 to 17, where p is the Laplace parameter.
Equations 10, 11, and 12 were solved analytically, resulting in Applying conditions 13 to 17 in ODE (10 to 12) and making the necessary substitutions, we have where T is the simulation time and M is an even value. The expression for the numerical inversion from the Laplace Transform is given by where the matrix is given by and must be the smallest integer of the operation.
Stehfest's method applies to simpler functions. In the case of a boundary value problem applied to diffusion, as in this work, the association of the finite difference method with the Stehfest algorithm and the Gauss-Siedel method is appropriate, to avoid numerical instabilities (DAVIES; CRANN, 1999).
With the definition of the equations for Laplace inversion using the methods described above, a numerical code was written in Fortran 77 language, in a free compiler (gFortran ® ) on Linux ® , where it was possible to obtain the graphical results (figures 2 to 6) for the oxygen concentration profile along the thickness of the polymeric film composed of three layers, which will be discussed in the next session.

RESULTS AND DISCUSSION
Two important parameters were tested to improve the stability of the The semi-analytical model was tested to predict the barrier capacity of multilayer films for application in food packaging. With the numerical solutions, it was possible to obtain the oxygen concentration profile along the film thickness for each selected configuration (ABA or ABC), represented in figures 2 to 6, with the vertical lines representing the interface between the layers. In figure 2, it can be seen that the oxygen concentration profiles in the outer layers showed the same behavior since it was the same material (ABA configuration), Polyethylene, which has the same diffusion coefficient (DPE: 1.0*10 -6 cm 2 /s), while the diffusion coefficient of the inner layer, composed of Ethylene Copolymer and Vinyl Alcohol, is 10 4 times lower (DEVOH: 7.2*10 -10 cm 2 /s). It is possible to observe a small fluctuation in the concentration profile at the interface between layers A and B of this configuration (figure 2) for curves with T = 0.50tmax, 0.75tmax, and tmax, probably due to limitations of the methodology used, which caused a numerical fluctuation in the concentration function when migrating material, because EVOH presents a greater oxygen barrier than PE. When the simulation time was 0.25tmax, the curve did not show oscillation.
The same symmetrical behavior for the concentration profile is observed in figure 3, which presents ABA configuration, but with the outer layers replaced by Polypropylene (DPP: 3.76*10 -7 cm 2 /s), while the inner layer remained with EVOH. There was a small change in the concentration profile, possibly due to the difference variation between the diffusion coefficients which, in this case, in the outer layer is at the range of 10 3 times greater than in the inner layer. In this configuration, the perceived oscillations in the layer's interfaces were softer, and in layer B the curves drawn for the different simulation times overlapped. The developed methodology was able to simulate the behavior of oxygen in the multilayer system without any variation in the concentration function.   In films where the configuration of the layers tested was ABC, an asymmetric graphic behavior was observed, since, in these cases, there are three materials with different diffusivities. Figure 6 represents this behavior, where PET, EVOH, and PE polymers have diffusion coefficients of 2.7*10 -9 cm 2 /s, 7.2*10 -10 cm 2 /s, and 1.0*10 -6 cm 2 /s, respectively. The same behavior was observed in a study using the commercial software based on the finite element method, Comsol ® , which presented PET and PE outer layers, as in the configuration of figure 6, however, the intermediate layer of the study included an oxygen scavenger, which reduced the oxygen concentration to zero in this layer (DI MAIO et al., 2017), while the present study contains an inert material, but with a lower diffusion coefficient than that of the material of the outer layers. The graphic behavior of the study can be associated with the curves obtained in the present study when compared in the shorter simulation times (0.25tmax and 0.50tmax).
In this last simulation (figure 6), it can be observed that there was also a small fluctuation in the first points of layer B (EVOH), where the oxygen concentration increases at the beginning of the layer and then decreases again. This behavior was observed for longer simulation times (0.75tmax and tmax), while to shorter simulation times the numerical instability is reduced, and it is possible to obtain oxygen concentration profiles without an oscillation in the function occurring when the diffusion is analyzed in multilayer materials. The concentration of oxygen at a point is calculated by the concentration given at a point before and a point afterward, by the Gauss-Siedel method. Therefore, at the layer interfaces, the concentration is calculated from the equation that corresponds to the layer of one material (A or B) and the concentration that corresponds to the layer of another material (B or C). At each iteration, the approximation found is tested to analyze if it can be considered the solution to the problem, by a numerical error analysis, which must be less than or equal to 10 -5 mol/m 3 . Also checking equations 21 to 23, obtained from the solutions in the Laplace space and used to calculate the oxygen concentration, there is the interaction between the data of each material in the neighboring material. Regarding this, we see the influence that one layer has on the other and this may be causing these oscillations, as observed in figures 2 and 6, as there are influences both from the analytical solutions obtained, which carry data from the nearby material, as well as from the Gauss-Siedel method. Reiterating that these oscillations occur for longer simulation times, therefore, the use of shorter numerical simulation times becomes more suitable for the convergence of results, generating oxygen concentration profiles without the numerical fluctuations in the concentration function.

CONCLUSIONS
The proposed methodology was able to simulate the oxygen concentration profile in selected polymeric materials using a high-level programming language, adopting a free compiler. This can still be improved since, for some configurations, numerical oscillations were detected when the simulation time was increased.
Therefore, a more detailed analysis must be carried out to reduce this method's instability. The Laplace Transform is widely used for solving differential equations and was used in the present work to solve the diffusion equation in a non-stationary state (Fick's 2nd Law). The use of a numerical simulation tool that can calculate its inverse extends its application to problems found in engineering.
The tool developed and used in this study has shown to be promising for application in the study of diffusion in multilayer polymeric systems and can be extended to describe other systems, with more layers than those manipulated in this