Optimization of space and time discretization during the finite element method application to multicomponent diffusion simulation = Otimização da discretização espaço-temporal do método de elementos finitos aplicado a simulação da difusão multicomponente

Optimization of space and time discretization during the finite element method application to multicomponent diffusion simulation 1*

3

1

Evandro Bona , Luiz Henry Monken e Silva , Dionísio Borsato , Rui Sérgio Ferreira da 2 2 1 Silva , Dayanne Aline de Souza Fidelis e Adriana Araújo 1

Departamento de Química, Universidade Estadual de Londrina, Londrina, Paraná, Brasil. 2Departamento de Tecnologia de Alimentos, Universidade Estadual de Londrina , Rodovia Celso Garcia Cid, Pr 445, Cx. Postal 6050, 86051 -990, Londrina, Paraná, Brasil. 3Faculdade de Apucarana, Apucarana, Paraná, Brasil. *Autor para correspondência. e -mail: [email protected]

ABSTRACT. The space and time discretization of the finite element method was optimized for following application in multicompon ent diffusion simulation during Prato cheese salting, a traditional and much consumed foodstuff in Brazil originated from the European Gouda cheese. It was ascertained that the correct choice of the time intervals and mesh is fundame ntal in applying the method. After optimization the simulated result s were in agreement with the experimental and calculated results by the analytical method, showing that the method is a promising tool for simulation of diffusive processes when two solutes are considered, and is also a much less restrictive technique than the analytical method. Key words: response surface methodology, simplex optimization, desirability functions, prato cheese salting.

RESUMO. Otimização da discretização espaço -temporal do método de elementos finitos aplicado a simulação da difusão multicom ponente. Neste trabalho foi realizada a otimização da discretização espaço -temporal do método de elementos finitos para sua posterior aplicação na simulação da difusão multicomponente durante a salda de queijo prato, um alimento tradicional e muito consumi do no Brasil e similar ao queijo Gouda. Foi verificado que a escolha correta dos intervalos de tempo e da malha é fundamental para a aplicação do método. Após a otimização os resultados simulados concordaram com os experimentais e estimados pelo método ana lítico. Mostrando que o método é uma ferramenta promissora para a simulação de processos difusivos quando dois solutos são considerados, além de ser uma técnica muito menos restritiva que o método analítico. Palavras-chave: metodologia de superfície de res postas, otimização simplex, funções de desejabilidade, salga de queijo prato.

Introduction Many people have avoided consuming Prato cheese to reduce sodium chloride ingestion, and consequently reduce problems related to arterial hypertension (Rapacci, 1989). Sodium consumption is one of the main factors that have been proven to cause increase in arterial pressure (He and MacGregor, 1999). Several substitutes for sodium chloride have been studied, especially potassium chloride, because it presents similar ph ysical properties (Lynch, 1987). Therefore, brine containing NaCl/KCl mixture in adequate proportions, to avoid sensorial problems, is used to produce a cheese with reduced sodium content (Rapacci, 1989; Zorrilla and Rubiolo, 1994; Katsiari et al., 1998). The NaCl/KCl/water system is classified as ternary Acta Sci. Technol.

(Nauman and Savoca, 2001) and the flow of potassium chloride and its influence on the flow of sodium chloride should be considered. Thus, a more extensive mathematical modeling is needed to simulate the diffusion process that occurs in this multicomponent system. Prato cheese is an example of a food with high nutritional value and common in the Brazilian diet, where the quantity and homogeneous distribution of sodium chloride are relevant for its final quality. Thus, studies involving sodium chloride diffusion, in the presence of a substitute (KCl) are fundamental in estimating parameters (such as salting and maturation time) essential for its industrial scale manufacture. The finite element method (FEM) has already been successfully used in simulating sodium chloride diffusion in Prato Maringá, v. 28, n. 2, p. 141-150, July/Dec., 2006

Bona et al.

142

cheese, using the Fick’s second law as a theoretical basis for the phenomenon (Silva et al., 1998; 1999). However, as far as it is known, there are no studies in the literature using FEM to simulate multicomponent diffusion in solid foods by the generalized Fick equation. According to Wang and Sun (2003), there are few publications in food area that have used FEM to simulate three-dimensional conditions as proposed in the present study. The objective of this study was to optimize the discretization used for three-dimensional simulation of multicomponent diffusion during mixed salting of Prato cheese in brine at rest. Material and methods Experimental procedure

Seven cheeses samples (Queijo Prato Lanche Di Carlo, Laticínios Campina Alta, Manoel Ribas , Estado do Paraná) were salted for 11 hours in a brine at rest containing 15.0 g NaCl/100 g NaCl + KCl + water and 5.6 g KCl/100 g NaCl + KCl + water at 10º 1ºC. Periodically two cylindrical samples were removed to estimate the NaCl and KCl concentrations (Figure 1) using a CELM FC-280 atomic emission photometer. The moisture, fat content and initial quantity of sodium chloride and potassium chloride were also determine d in a sample that was not salted (Bona et al., 2005).

A x is Y : 1 0c m

5 cm 1 .5 c m

4 cm

4 cm

A x is Z : 1 2 c m Figure 1. Dimensions of the Prato cheese used (axis X, hidden in the Figure, measured 4cm) and sampling adopted to study multicomponent diffusion.

 R3,   [-R1, R1]x[-R2, R2]x[-R3, R3], associated to a system of Cartesian coordinates x, y, z with origin located in geometric center of the Prato cheese. The system under study was considered not reactional, bearing in mind that the most significant chemical reactions occurred during the maturation period (Zorrilla and Rubiolo, 1998). The diffusion coefficient or mass diffusivity was considered constant in relation to the concentration (regardless of time and position). The external resistance was equal for the two solutes, because they are very similar ionic compounds (Zorrilla and Rubiolo, 1994). The sample contraction was considered negligible, because according to the literature the variation in volume is minimal during cheese salting (Silva et al., 1998; Gerla and Rubiolo, 2003). The process is done under very approximately isothermal conditions. Under these conditions it was proposed to analyze the cheese salting process considering it immersed in an unstirred aqueous solution containing NaCl and KCl in the mass proportion recommended by Rapacci (1989) to ensure good sensorial acceptability. The C1(x,y,z,t) and C2(x,y,z,t), concentration of the NaCl and KCl solutes, respectively, at a point P(x,y,z)   and at instant t, can be described by Onsager (1945) equations for the solute concentrations:

C1  D11 2C1  D12 2C2 t C2  D21 2C1  D22 2C2 t

(1)

where Dii are the main coefficients, Dij crossed coefficients, that combine the flows and 2(.) = .(.), is the Laplacian operator. In the salting process, the initial conditions are given by

C1 ( x, y, z ,0)  C1, 0 C 2 ( x, y, z ,0)  C 2, 0

x, y , z  

(2)

Three-dimensional modeling of the mult icomponent diffusion during salting by the finite element method

where C1,0 and C2,0 are known. The Cauchy boundary condition for no stirred brine were (Luna and Bressan, 1986):

The following considerations or simplifying hypotheses were used to formulate the simultaneous diffusion in finite elements: The diffusion of two solutes was modeled in a three-dimensional cheese that occupied the volume 

C1  R, t  hm C1  C1,s    m with x, y, z  , t0 C2  R, t  hm C2  C2,s    m

Acta Sci. Technol.

(3)

Maringá, v. 28, n. 2, p. 141-150, July/Dec., 2006

Optimization of finite element method applied to diffusion

where  is the set of points on the surface that wraps the Prato cheese; C1,s and C2,s are the solutes concentrations in brine; hm (g cm-2 h) is the mass transfer coefficient; λm (g cm-1 h) is the mass conductivity and ∂/∂η is the normal derivative operator. The hm and m coefficients are related to Biot mass-exchange number by

Bi 

hm  Ri ; i = 1,2,3. m

(4)

where Ri is the characteristic length (cm). The diffusion coefficients and the Biot massexchange number (Bona et al., 2005) were estimated (Table 1) from the experimental data obtained and by a one-dimensional analytical model for multicomponent diffusion in a brine at rest. Table 1. Fitted values for the parameters of the analytical solution (Bona et al., 2005). NaCl 0.225 (D 11) 0.027 (D 12)

Main coefficients (cm 2/day) Cross coefficients (cm 2/day) Biot mass-exchange number hm/λm (cm-1)

40.660* 20.330

KCl 0.240 (D 22) 0.045 (D 21)

* Related to X axis.

The system of partial differential equations formed by equation (1) and by the conditions (2) and (3), even with the simplifying hypotheses adopted in the

143

process, it is difficult to be resolved analytically. In this situation a transformation obtained with the eigenvalues and eigenvectors of the ma trix of the differential operator that appears on the right side of equation (1) was used (Bona et al., 2005). The disadvantages of this method are the restrictions in the geometry, in the process and further the fact that the transformed variables have no physical meaning. An alternative that does not have the se restrictions that can be applied in more general problems like this situation is the numerical solution of the system by the finite element method (FEM). Therefore the generalized Galerkin formulation was used (Huebner et al., 1995) to obtain the expression to be descritized by finite elements whose main steps are summarized as follows: For system (1) the expression of weighted residues formed is:  C1   0  t   0 D11 D12 2C1  (5) d        0   C2  0  D21 D222C2 d,    t 

where Φ and Ψ are weighting functions that are null in ∂Ω b) The right side of the equation (5) is integrated in parts, obtaining, after considering the boundary conditions (3), the following:

 C1   0   t   0    C2 d     t 

   x    0 

hm   m



0  x

 y 0

 0   D11 0    D    21

0  y

 z 0

 D11    D21 0   0     0  z   0   0

D12  C1  hm  ds   D22  C2  m

D12 D22 0 0 0 0

0 0 D11 D21 0 0

 0   D11 0    D    21

0 0 D12 D22 0 0

0 0 0 0 D11 D21

 C1   x   C  0  2  x 0   C   1 0   y   d 0   C2  D12   y    C  D22   1   Cz   2  z 

(6)

D12  C1,s   ds D22  C2,s 

where ∂Ω is the boundary of Ω and C1,s, C2,s are Acta Sci. Technol.

Maringá, v. 28, n. 2, p. 141-150, July/Dec., 2006

Bona et al.

144

known values.





f ( )d 

and





f ( s ) ds

are,

respectively, integrals of volume over the dominion and boundary. The discretization continued and the next steps were identified as follows (Singh, 1983): c) Partition of the domain () in subdomains called finite elements is adopted. d) In each finite element (with N nodes and two

A 2N  2N

 (7) A2N2NC 2N2N  B2N2NC2N F2N2NC2N G2N where,

0   1 0   1   2 0     0 2 0   N 0     0  2  1 d  0  1 0  2  0  N         0  N  0   N  

1 1 1  0 0 0   x  y  z    0 1 0 1 0 1   x y z D11 D12 0   2 2  2 0 0 0 D21 D22 0  x  0 0 D y z 11  B2N2N         2 2 2  0 0 0   0 0 D21 x y z   0 0 0         N N N  0 0 0 0 0 0    x  y  z    0 N 0 N 0 N   x y z  

C 2N

degrees of freedom per node) the variables of the problem were locally interpolated. Therefore, equation (6) can be rewritten, in matrix form, for each element as:

 C1  C   2     F   C  2N  2N  1 C2 

 1 0  2 hm    m  0     N 0 

0   1  0  D  2   11 D    21  0   N 

0

C

0

0 0 D12 0 D22 0 0 D11 0 D21

2N

 C1  t  C 2   t     C1  t  C 2   t

          

N 2 1  0   x 0 x 0  x  N  1 2 0  0 0  0  x x x  0 1  N 2 0 0  0  0  y  y  y d  N  1 2 0  0 0  0  y y y  D12   1  2 D22 0 0  N 0   z  z z  N  1 2 0  0 0  z z z  

D12  1 0 D 22   0  1

2 0

0 2

 

N 0

0  ds  N 

0   1 0    C1, s  1   C  2  0 2, s   0  0   0 D D   h    2 N 12   1  G 2N   m  0  2   11 ds   0    0   0  D D 1 2 N  21 22    m   C    1, s    0  C2, s   N 0   N  e) The global interpolation is formed by combining the elements, resulting in a system of Acta Sci. Technol.

Maringá, v. 28, n. 2, p. 141-150, July/Dec., 2006

Optimization of finite element method applied to diffusion

145

ordinary differential equations. f) The solution of this system is obtained by an implicit process of finite differences, supplying the values of the concentrations in the element nodes. The concentration profiles in any part (point) of the Prato cheese at any time could be determined by local interpolation.

Time discretization

Spatial discretization

 C1  C02N   = A 2Nx2N  2N Δt   1 0  C  C2N   - B 2Nx2N  2N 2  

The spatial domain was represented by a set of serendipity hexahedron elements of the C0, type. Each element had twenty nodes with two degrees of freedom (one for each solute), distributed at the edges and vertices of its external surface, totalizing forty degrees of freedom (Brebbia and Ferrante, 1975; Silva et al., 1998). The element, defined in standard Cartesian coordinates -1 ≤ ξ, η, ζ ≤ 1 (Chung, 1978), the arrangement of the nodes and their local enumeration are shown in Figure 2. The interpolation functions preserve the continuity between the elements and ar e deduced from an incomplete quadratic polynomial, following the same procedure to obtain the Lagrange functions (Chung, 1978; Akin, 1982). In the formulation proposed i = ψi was used. 19

8 20

7 18

17 6

5

15

16

3

12

1

10

9

2

(a)

ζ=±1 6

3 4

2

1

ξ=±1

η=±1 5

(b)

Figure 2. (a) Standard hexahedron element, with 20 nodal points. (b) Orientation of the standard hexahedron element and enumeration of the surfaces.

Acta Sci. Technol.

as follows

(8)

 C1  C02N    G 2N  F2N 2N  2N 2   where, Δt is a discrete increase in time, C 02N is the solute concentration vector, already known, of the element before the adding of Δt and C12N represents the concentration vector that will be calculated, after adding time (Δt). Organizing the terms of the equation (8) it can be shown that:

 A 2Nx2N B 2Nx2N  F2N  2N  1   C 2N = 2  t   A 2Nx2N F2N  2N  B 2Nx2N  0   C 2N  G 2N 2  t 

(9)

14

11

4

13

The Crank-Nicholson scheme or central difference method was used for the dominion (Ω) and boundary (∂Ω) (Bickford, 1990). In this problem the  and C , equation (7), were discretized vectors C 2N 2N

Resolution of the linear equation system

The LU decomposition method (Sperandio et al., 2003) was used to solve the proposed system of linear equations. The computational implementation of the method was carried out by a set of subroutines (DLSLRB and DLFCRB) based on the LINPACK package (Cline et al., 1979; Dongarra et al., 1979). The DLSLRB subroutine was used to solve the linear equation system stored in band without using interactive refinement. This procedure uses the subroutines: DLFCRB that estimates the matrix condition number and performs the LU decomposition; DLFSRB that solves the system formed by the matrixes L and U. Mesh

A mesh generating subroutine was used because in practical applications the finite element method requires many input data. The global enumeration and the node coordinates of the mesh were obtained Maringá, v. 28, n. 2, p. 141-150, July/Dec., 2006

Bona et al.

146

automatically for rectangular dom ains. This procedure facilitates the identification of each node position at the element enabling sampling at any point or region of the mesh (Silva et al., 1998).

produced the best combination between deviation and computer time was chosen. At this stage all the simulations were carried out using a mesh of seven elements on each axis (totaling 343 elements).

Integration and interpolation of the simulation results

Table 2. Assessment of the intervals tested for time discretization.

The experimental values obtained represented the mean NaCl and KCl concentration for the cylindrical sample (Figure 1). However, the finite element method supplies point concentrations, therefore, the simulated points had to be interpolated onto the imaginary axis located at the center of the sample. The mean was calculated by integrating these points according to equation (10).

C 

0.05 2 2

Time interval (h) 0.20 0.30 0.50 1.00 11 22 5 10 3 10 3 2

Percent deviation in Computer Relation to experimental data Time 2.00 NaCl KCl (minutes)* 21.47% 11.28% 12.70 17.26% 10.06% 26.79 17.23% 10.23% 19.04 17.29% 10.24% 18.45 4 17.30% 10.25% 14.52

*For a microcomputer with Pentium III 800MHz processor and 256 Mbytes RAM memory.

 C x, y, z, t dV ; i = 1.2  dV V

Number of loops in each

(10)

i

The computational program for simulation of multicomponent diffusion, called Simul 3.0, was developed in Fortran language. The minimal requirement for his use is a Pentium III processor and 256 Mb of RAM memory.

To assess the efficiency of the spatial discretization, which depends on the quantity of elements and their distribution, preliminary tests were carried out. This previous evaluation indicated that at least six elements per axis are necessary to simulate the proposed system. Otherwise, the results are far away from the experimental values. Further, due to the restrictions imposed by the microcomputer used, a t most ten elements per axis could be used. An experimental design (conveniently randomized) of the BoxBehnken type with a central point (Table 3) was applied to determine the quantity of elements to be used on each axis, to minimize the deviations and computer time.

Statistical test

Table 3. Results obtained using the Box-Behnken design.

V

The integral was calculated numerically by a Gauss formula (Sperandio et al., 2003) with the abscissas and weights obtained according to Davis and Polonsky (1965). The computer program

The concentrations calculated by FEM, the analytical method (Bona et al., 2005) and the experimental method were compared by percentage deviation (Heldman, 1974) to assess the quality of fit. 2

%deviation  100

 Ccalc  Cexp   1     C   N 1 i 1  exp i   N

(11)

where, C calc is the mean concentration estimated by the numerical solution; C exp is the experimental mean concentration and N is the number of observations considered. Optimization of the space and time discretization

The objective of optimization was to determine values for the discrete time and space intervals that permit to minimize the deviations without significant increase in computer time. Some tests were performed for time discretization (Table 2) with different time intervals. The set of time intervals that Acta Sci. Technol.

Factors or inputs Elements per axis* X Y 6 (-1) 6 (-1) 10 (1) 6 (-1) 6 (-1) 10 (1) 10 (1) 10 (1) 6 (-1) 8 (0) 10 (1) 8 (0) 6 (-1) 8 (0) 10 (1) 8 (0) 8 (0) 6 (-1) 8 (0) 10 (1) 8 (0) 6 (-1) 8 (0) 10 (1) 8 (0) 8 (0)

Z 8 (0) 8 (0) 8 (0) 8 (0) 6 (-1) 6 (-1) 10 (1) 10 (1) 6 (-1) 6 (-1) 10 (1) 10 (1) 8 (0)

Responses or outputs % Deviation relative to experimental data Computer NaCl (%) KCl (%) time (min)** 12.26 7.60 10.67 7.62 6.87 17.53 8.32 7.63 40.17 3.54 6.79 76.88 22.11 12.12 11.01 16.11 10.00 17.64 7.83 7.19 39.44 3.33 6.29 72.30 19.74 11.47 7.63 16.23 11.01 28.09 6.28 7.01 24.66 2.54 7.11 117.97 3.43 7.15 30.43

*The codified levels are in parenthesis. **For a microcomputer with Pentium III 800 MHz processor and 256 Mbytes RAM memory. A time discretization with two 0.05h, three 0.30h, two 1.00h and four 2.00h intervals were adopted.

The replications at the central point were not included, because the responses are computer calculation and thus would be the same. This type of design was chosen because only three levels of variation were needed for the independent variables and is possible to fit a quadratic model (Box and Draper, Maringá, v. 28, n. 2, p. 141-150, July/Dec., 2006

Optimization of finite element method applied to diffusion

1987). From the results obtained for the design, regression models were fitted relating the number of elements on each axis with the dependent variables : computer time and deviation among the calculated and experimental results (StatSoft, 2005). The models obtained were then submitted to simultaneous optimization by the Simplex method (Bona et al., 2000) using the desirability functions (Derringer and Suich, 1980) to determine the best quantity of elements on each axis of the mesh. Results and discussion Table 2 shows the most promising results selected from the different tests to assess the best time discretization. The deviation compared to the experimental values was greater than 10% in all the attempts made because a poor mesh was used (with 343 elements). However, at this stage, the objective was to assess only the influence of time discretization. Table 2 shows that shorter time intervals are necessary at the start of the simulation, but after the first step, the intervals can be increased (Lyra, 1993). An indication of the importance of using short time intervals at the start of the simulation can be obtained from the profile of s alt distribution along X axis (Figure 3).

147

et al., 2002), was damping when short time intervals were used in the beginning of the simulation. This fact was justified by the great salt concentration gradient between the brine and the cheese in the beginning of the process. Also, when short time intervals were used, the increase in concentration is reduced for each loop of the simulation. In addition, the results obtained (Table 2) indicated that the increase in the number of discretization intervals did not ensure a continuous reduction in the deviation. In fact, a convergence to values of around 17% was observed for NaCl and 10% for KCl. A greater reduction in the deviations could only be obtained by modifying the mesh. Therefore, discretization time was chosen with refinement in the initial time intervals and less computer time (last line of Table 2) and, it was used to optimize the mesh (Table 3). The regression models were attained (Table 4) from the results obtained to the proposed design (Table 3) in function of the number of elements per axis Table 4. Quadratic regression models fitted to the data obtained by the Box-Behnken design. % Deviation from NaCl (Computer time )-1/2 b experimental data a Model Coefficients Standard error Coefficients Standard error Intercept 223.6300 5.5618 1.8680 0.0711 X -13.3000 0.6197 -0.0680 0.0077 Y -7.6694 0.5798 -0.1383 0.0077 Z -28.4931 0.6197 -0.1346 0.0077 2 X 0.7066 0.0362 0.0020 0.0004 2 Y 0.4197 0.0362 0.0037 0.0004 2 Z 1.5222 0.0362 0.0037 0.0004 X*Y 0.0015 0.0003 X*Z 0.0937 0.0273 0.0014 0.0003 Y*Z 0.0040 0.0003 2 R (adjusted) 0.9989 0.9950 a

The model in relation to deviations of potassium experiments were not considered because they were highly correlated with the sodium deviations. b The computer time response was presented as a Box Cox transformation to stabilize variance.

Figure 3. Influence of the time discretization on the simulation of the salt distribution profile, in function of distance. The solid curve (—) represents the simulation of the first hour by a single time interval. In the dotted curve ( ---) two 0.05h loops and three 0.30h loops were used for the same simulation.

For the simulations showed in Table 2, it was observed that, except for the attempt with eleven onehour loops, all the others presented a similar profile. One can see in Figure 3, which depict the simulation curves for the first and the last line showed in Table 2, that the initial oscillatory phenomenon in the results, that is inherent to the finite element method (Britz Acta Sci. Technol.

Deviation relative to experimental data for NaCl and KCl were highly correlated (r=0.98), and therefore the use of only one of these dependent variables is recommended in the simultaneous optimization (Peterson, 2004). Thus, when one is minimized the other one will consequently be minimized as well. As the deviation in relation to the experimental NaCl values resulted in a better model, this response was adopted to carry out the simultaneous optimization. Regarding the computer time, it was observed this model violated one of the assumptions of the variance analysis. A Box-Cox power transformation was used to stabilize the variance (Box and Draper, 1987). The quality of the fitted models can be assessed in Figure 4. It can also be observed by the fitted coefficients that the deviation for the experimental values underwent greater Maringá, v. 28, n. 2, p. 141-150, July/Dec., 2006

Bona et al.

148

reduction with the increase in the number of elements on the X and Z axis. This observation is in line with the characteristics of the mesh adopted. Since the concentrations calculated and used for comparison with the experimental values were taken along the X axis, a greater refinement on this was important to reduce the deviations. The number of elements on Z axis, due to the characteristics of the subroutine used to generate the mesh, is mainly responsible for the matrix dimension of the linear equation system. Thus, increasing the number of elements on this axis also increases both numbers of equations used as the precision of the method. Regarding computer time, it was observed that an increase in the number of elements, increased the computing time, especially for Y and Z axis. This performance was also explained by the characteristics of the mesh generator subroutine.

optimization of the space and time discretization. The oscillatory damping during the salting first hour simulation could be observed in Figure 5. Table 5. Results obtained for optimized mesh with 8 elements on X axis, 7 on Y axis and 9 on Z axis. Model Deviation from experimental NaCl concentration 2.94 2.86 2.80

Response Predict Observed Error (%)

Computer time (min) 28.61 28.31 1.06

Figure 5. Influence of the spatial discretization on the simulation of the salt distribution profile, in function of distance. The solid curve (—) represents the simulation of the salting first hour using a mesh with seven elements in each axis. In the dotted curve (---) was used optimized mesh with eight elements along X axis, seven along Y and nine along Z axis.

Acta Sci. Technol.

2,75 (water + NaCl + KCl)

As can be observed from the characteristics of the fitted equations, the system led to an interesting problem of optimization. An increase in the number of elements led to a reduction in deviation (desirable), but also implied an increase in computer time (undesirable). This problem was solved using the Derringer and Suich desirability functions (Derringer and Suich, 1980) , which were adjusted to a maximum computer time value (30 minut es) and to minimize the deviation compared to the experimental results of sodium chloride. The general desirability was then optimized by the Simplex method (Bona et al., 2000). As a result of the optimization, a mesh was obtained with eight elements on X axis, seven on Y axis and nine on Z axis, totaling 504 elements (Table 5). Comparison of the dotted curve in Figure 5 and the solid curve in Figure 3 (for NaCl and KCl) shows the evident damping of the oscillation obtained by the

A comparison between the results obtained by the analytical method (AM) and experimental values (Bona et al., 2005) and those from finite element method (FEM) are shown in Figure 6.

Salt Concentration, g/100g

Figure 4. Observed versus predictive values by the regression models.

2,50 2,25 2,00 1,75 1,50 1,25 1,00 0,75 0,50 1

2

3

4

5

6

7

8

9

10

11

Salting time (h)

Figure 6. Comparison of the experimental results (dots), analytical method (dotted line) and the finite element method (solid line) after optimization of the spatial and time discretization.

It was verified that after optimization, the FEM obtained a good agreement with the experimental results (deviations of 2.86% for NaCl and 7.09% for Maringá, v. 28, n. 2, p. 141-150, July/Dec., 2006

Optimization of finite element method applied to diffusion

KCl) and the simulated results by AM (deviations of 3.14% for NaCl and 2.74% for KCl). Comparing the deviations between the optimized mesh with 504 elements (see the deviations above) and the starting mesh with 343 elements (see the last line of Table 2) there was a substantial reduction. It makes clear the importance of the optimization performed. Conclusion The proper choice of the time intervals and the mesh permitted the optimization of the results obtained by the finite element method ( FEM) by damping the oscillation without significant increase in computer time. Once again the importance of these parameters was confirmed in the application of FEM. The use of short time intervals was important to damp the inherent oscillation in the first loops of the simulation. However, this reduction was restricted and should be accompanied by a correct fit of the mesh. After optimizing the results simulated by the FEM, these were very close to the experimental values or values estimated by the analytical method (AM). Therefore, the use of FEM is recommended to simulate multicomponent diffusion because it presents three-dimensional characteristics (that makes it more appropriate for the cheese maturation stage) and also permits the use of various boundary conditions without great changes in the basic formulation. Furthermore, it enables study with the other geometry and permits simultaneous simulation of other parameters (for example, heat transfer). Acknowledgements Authors thank Fundação Arauc ária for the financial support and Capes for scholarships. References AKIN, J.E. Application and implementation of finite element methods. London: Academic Press, 1982. BICKFORD, W.B. A first course in the finite element method. Boston: Irwin, 1990. BONA, E. et al. Aplicativo para otimização empregando o método simplex seqüencial. Acta Sci., Maringá, v. 22, n. 5, p. 1201-1206, 2000. BONA. E. et al. Difusão multicomponente durante a salga mista de queijo prato. Cienc. Tecnol. Alim., Campinas, v. 25, n. 2, p. 394-400, 2005. BOX, G.E.P.; DRAPER, N.R. Empirical model-building and response surfaces. New York: Wiley, 1987. BREBBIA, C.A.; FERRANTE, A.J. The finite element technique: an introduction for engineers . Porto Alegre: Ed. UFRGS, 1975 Acta Sci. Technol.

149

BRITZ, D. et al. Damping of Crank-Nicolson error oscillations. Comp. Biol. Chem., Luton, v. 27, p. 253-263, 2002. CHUNG, T.J. Finite element analysis in fluid dynamics . New York: McGraw Hill, 1978. CLINE, A.K. et al. An estimate for the condition number of a matrix. SIAM J. Num. Anal., Essex, v. 16, p. 368-375, 1979. DAVIS, P.J.; POLONSKY, I. Numerical interpolation, differentiation and integration. In: ABRAMOWITZ, M.; STEGUN, I. (Ed.). Handbook of mathematical functions , New York: Dover Publications, 1965 . cap. 25, p. 875-924. DERRINGER, G.; SUICH, R. Simultaneous optimization of several response variables . J. Qual. Technol., Milwaukee, v. 12, n. 4, p. 214-219, 1980 DONGARRA, J.J. et al. LINPACK users’ guide . Philadelphia: SIAM, 1979. GERLA, P.E.; RUBIOLO, A.C. A model for determination of multicomponent dif fusion coefficients in foods. J. Food Eng., Davis, v. 56, p. 401-410, 2003. HE, F.J.; MACGREGOR, G.A. Potassium intake and blood pressure. Am. J. Hypert., New York, v. 12, p. 849851, 1999. HELDMAN, D.R. Predicting the relationships between unfrozen water fraction and temperature during food freezing using freezing point depression. Trans. ASAE, St. Joseph, v. 17, p. 63, 1974. HUEBNER, K.H. et al. The finite element method for engineers. 3. ed. New York: John Wiley & Sons, 1995 . KATSIARI, M.C. Manufacture of kefalograviera cheese with less sodium by partial replacement of NaCl with KCl. Food Chem., WhitKnigts, v. 61, n. 1/2, p. 63-70, 1998. LUNA, J.A.; BRESSAN, J.A. Mass transfer during brining of cuartirolo argentino cheese. J. Food Sci., Chicago, v. 51, n. 3, p. 829-831, 1986. LYNCH, N. In search of the salty taste. Food Technol., Chicago, v. 41, n. 11, p. 82-86, 1987. LYRA, P.R.M. Finite elements analysis of parabolic problems: combined influence of adaptive mesh refinement and automatic time step control. Rev. Bras. Cienc. Mec., Rio de Janeiro, v. 15, n. 2, p. 172-198, 1993. NAUMAN, E.B.; SAVOCA, J. An engineering approach to an unsolved problem in multicomponent diffusion. AIChE J., Newarke, v. 47, n. 5, p. 1016-1021, 2001. ONSAGER, L. Theories and problem s of liquid diffusion. Ann. N.Y. Acad. Sci., New York, v. 46, p. 241, 1945. PETERSON, J.J. A posterior predictive approach to multiple response surface optimization. J. Qual. Technol., Milwauke, v. 36, n. 2, p. 139-153, 2004. RAPACCI, M. Efeito da substituição do cloreto de sódio por cloreto de potássio nas características do queijo prato . 1989. Dissertação (Mestrado)-Universidade Estadual de Londrina, Londrina, 1989. SILVA, R.S.S.F. et al. Difusão do cloreto de sódio no processo de salga de queijos: modela gem matemática com Maringá, v. 28, n. 2, p. 141-150, July/Dec., 2006

Bona et al.

150 o emprego do método de elementos finitos. Cienc. Tecnol. Alim., Campinas, v. 18, n. 1, p. 87-92, 1998. SILVA, L.H.M. et al. Finite element analysis of diffusion problem during cheese salting: combined influence of space and time discretization. Acta Sci., Maringá, v. 21, n. 4, p. 873-879, 1999. SINGH, R.P. Numerical Techniques. In: SAGUY, I. (Ed.). Computer-aided techniques in food technology . New York: Marcel Dekker, 1983. cap. 2, p. 23-47. SPERANDIO, D. et al. Cálculo numérico: Características matemáticas e computacionais dos métodos numéricos. São Paulo: Prentice Hall, 2003. STATSOFT, Inc. Statistica (data analysis software system), version 7.1. www.statsoft.com., 2005. ZORRILLA, S.E.; RUBIOLO, A.C. A model for using

Acta Sci. Technol.

the diffusion cel in the determination of multicomponent diffusion coefficients in gels or food. Chem. Eng. Sci., Berkeley, v. 49, n. 13, p. 2123-2128, 1994. ZORRILLA, S.E.; RUBIOLO, A.C. Modelling average concentrations of salt and salt substitute in partial or total volumes of semihard cylindrical cheeses. Int. J. Food Sci. Technol., London, v. 33, p. 501-508, 1998. WANG, L.; SUN, D. Recent developments in numerical modelling of heating or cooling processes in the food industry – a review. Trends Food Sci. Technol ., Norwich, v. 14, p. 408-423, 2003.

Received on February 02, 2006. Accepted on November 28, 2006.

Maringá, v. 28, n. 2, p. 141-150, July/Dec., 2006