Computational Simulation of Heat Transfer in a Dip Shrink Tank Using Two Different Arrangements of Electrical Resistances

Abstract

Biohazard recontamination of food can occur in a meat processing plant during slicing, portioning, or racking. Subsequently, to protect them from external agents, not allow the loss of moisture contained in the product, and preserve its safety, they undergo a shrinking process; which consists of submerging in a tank with hot water at an approximate temperature of 87°C, for a certain time the food that has been wrapped with a heat shrink plastic, making it shrink. In this work, the behavior of heat transfer in a non-commercial shrink tank, built with two different arrangements of electrical resistances for water heating, is investigated. The study was carried out through numerical simulations with the implicit method of alternating directions (ADI). The results obtained from the heating times with their respective temperature distributions show that the arrangement with four resistances is the most efficient for the process of heating the water in the shrink tank, achieving a homogeneous temperature of 87°C, in times less than 9 minutes with a heat flux of q = 24.48 W. The validation of the simulations will be carried out in a subsequent work with experimental tests carried out in the shrink tank.

1. Introduction

The packaging of a food product consists of placing a material known as packaging, which completely covers the object, which can generally have two functions, the main and most important is to protect it and the second is the visual impact for commercial purposes.

Many foods require a container that guarantees their conservation, taking into account their different properties (liquid, solid, gel, pH, among others) and composition (proteins, lipids, vitamins, etc.). For this purpose, various materials have been used, which could have harmful properties for food, if the correct material and properties are not chosen.

The heat shrink process consists of immersing food previously covered in heat shrink plastic (in general, the material used for this process is polyethylene), in a tank with hot water, for a certain time with the intention of protecting it from external agents [1].

In many meat derivatives processing plants, a problem associated with product recontamination after its primary heat treatment has been detected, mainly associated with the use of a cooking sleeve that is changed for marketing the product, as well as the removal of that sleeve to distribute the product sliced or in portions and vacuum packed. Subsequently, these presentations are subjected to the shrinking process, preceded by some other procedure that allows control of recontamination, for example, the treatment by High Hydrostatic Pressures (HPP) [2].

For this reason, it is essential to obtain the normal operating conditions of temperature and time in the shrink tank, which allow increasing and sustaining these variables sufficiently to control the biological hazards that are identified.

In this work, the heat transfer was investigated in a shrink tank built in the facilities of the Faculty of Higher Studies Cuautitlán—UNAM, with two different arrangements of electrical resistances for water heating, in order to select an arrangement of resistances to ensure adequate heat transfer within the tank.

The behavior of heat transport that was studied is convection, which is considered an improved or modified form of conduction, in which a massive movement of the medium is also present [3]. There are two types of mechanisms for convection, free and forced. In this investigation we worked with free convection.

To model the shrink tub system, a rectangular cavity open at the top in two dimensions was considered, with two different configurations of electrical resistances, which serve as heat sources.

The resulting system of equations was discretized and a finite difference scheme was also used, which was solved using the Method of Alternate Directions Implicit (ADI), obtaining as a result the profiles of the stream function φ and vorticity ω.

2. Theory and models

2.1 Governing equations

Conduction, convection, and radiation are identified as ways to transfer heat. The mechanism of conductive heat transfer can be appreciated by heating a material with a heat source, as shown in Figure 1.

Figure 1.

Heating a material.

In the previous figure, there is a higher temperature on the left side of the material due to the heat source, and on the right side a lower temperature; Therefore, this temperature difference will result in heat transport by conduction in the material.

The heat flow is proportional to the area and the temperature difference, and can be quantified by means of Fourier’s first law.

$$\mathrm{Q}=\mathrm{kA}\left(-\frac{\mathrm{dT}}{\mathrm{dx}}\right)$$E1

In Eq. (1), k is the thermal conductivity and depends on the material, A is the cross-sectional area and (dT/dx) is known as the temperature gradient, and the negative sign indicates that the heat flux is in the opposite direction of the temperature gradient [4].

If one wanted to quantify the heat flux, Eq. (1) becomes

$$\mathrm{q}=\mathrm{k}\left(-\frac{\mathrm{dT}}{\mathrm{dx}}\right)$$E2

For the case of transient heat transport, Joseph Fourier developed a second law, which considers the variation of temperature with respect to time and is given by

$$\mathrm{\rho }{\mathrm{C}}_{\mathrm{p}}\frac{\mathrm{\partial T}}{\mathrm{\partial t}}=\mathrm{k}{\nabla }^2\mathrm{T}$$E3

In Eq. (3), ρ corresponds to the density and C_p to the heat capacity at constant pressure.

Regarding the mechanism of heat transport by convection, it can be seen by heating a metal container containing water, as shown in Figure 2.

Figure 2.

Heating of a container containing water.

In the previous figure, it can be seen that the fire produced by the combustion of the gas in the stove heats the container and this, in turn, heats the water. At a certain time, the fluid at the bottom of the container will have a lower density compared to the fluid near the surface, this is due to its thermal expansion. Therefore, the liquid at the bottom tends to rise and the liquid on the surface will go down, resulting in a heat transfer when mixed. Therefore, the phenomenon of convection occurs when mixing relatively hot portions of fluid with cold ones. The equation with which the heat flow by convection is quantified is Newton’s law of cooling.

$$\mathrm{Q}=-\mathrm{hA}∆\mathrm{T}$$E4

In Eq. (4), h is the convective coefficient, A is the area, and ΔT the temperature difference [5].

As for radiation, it refers to the radiant energy from a source to a receiver; in which part of the energy is absorbed by the receiver and part reflected by it. Boltzmann established an equation for the flow of heat by radiation

$$\mathrm{Q}=\mathrm{\epsilon \sigma A}{\mathrm{T}}^4$$E5

In Eq. (5), σ is the dimensionless Boltzmann constant and ε the emissivity [3].

For the case in which internal heat generation and convective heat transport are present, Eq. (3) takes the following form

$$\mathrm{\rho }{\mathrm{C}}_{\mathrm{p}}\frac{\mathrm{\partial T}}{\mathrm{\partial t}}=\mathrm{k}{\nabla }^2\mathrm{T}+\mathrm{G}-\nabla ·\mathrm{\rho }{\mathrm{C}}_{\mathrm{p}}\left(\mathrm{T}-{\mathrm{T}}_{\mathrm{ref}}\right)\overrightarrow{\mathrm{v}}$$E6

Eq. (6) contemplates the heat flow in a transitory state, the conduction heat mechanism (k∇²T), the internal heat generation (G) due to electrical resistances and convection $\nabla ·\mathrm{\rho }{\mathrm{C}}_{\mathrm{p}}\left(\mathrm{T}-{\mathrm{T}}_{\mathrm{ref}}\right)\overrightarrow{\mathrm{v}}$, being the velocity and T_ref, a reference temperature.

In the case of fluid mechanics, there are two equations that are essential to characterize the flow of fluids, these are the continuity equation and the Navier-Stokes equation, which are expressed in a vector formulation as follows [6]:

$$\frac{\mathrm{\partial \rho }}{\mathrm{\partial t}}=\nabla ·\mathrm{\rho }\overrightarrow{\mathrm{v}}$$E7

Eq. (7) corresponds to the continuity equation, in which ρ is the density and $\overrightarrow{\mathrm{v}}$, is the velocity.

$$\mathrm{\rho }\left[\frac{\partial \overrightarrow{\mathrm{v}}}{\partial t}+\overrightarrow{\mathrm{v}}·\nabla \overrightarrow{\mathrm{v}}\right]=\mathrm{\mu }{\nabla }^2\overrightarrow{\mathrm{v}}-\nabla \mathrm{P}+\mathrm{\rho }\overrightarrow{\mathrm{g}}$$E8

Eq. (8) is the well-known Navier-Stokes equation, the term μ denotes the viscosity, P is the pressure, and $\overrightarrow{\mathrm{g}}$ is the acceleration due to gravity.

2.2 Materials

The shrink tank used for the simulations can be seen in Figure 3, which was built with 304 stainless steel material for the parts that are in contact with the meat product. The dimensions of the tank are 0.728 m × 0.420 m.

Figure 3.

Shrink tank.

2.3 Modeling

For the numerical modeling, a two-dimensional rectangular cavity open at the top was considered, with the following configurations for the placement of electrical resistances as heat sources in red, as shown in Figure 4.

Figure 4.

(a) Arrangement of two resistors and (b) arrangement of four resistors for the rectangular cavity.

The cavity has a height H and a length L, a relationship between both magnitudes was considered as L = 3H. In the same way, the heat sources that are indicated with red color in Figure 4, were considered as 1/9 of L.

For the first configuration, the distance from the walls to the heat sources was estimated as 1/3 of L. For the second configuration, the distance from the wall to the heat sources is 1/9 of L and the distance between the different sources of heat is 1/9 of L.

All heat sources are considered equal and constantly emit heat.

For modeling, the temperatures (T₀) of the cavity walls and the upper zone were considered constant. The working fluid is water and no induced flow or pressure gradient is considered, in addition the upper zone is considered open to the environment, and it is assumed that there is no fluid exchange with the external environment.

It is considered that, due to temperature gradients, the fluid experiences density changes small enough to support the hypothesis of an incompressible fluid but large enough to produce a convection phenomenon, for which the Boussinesq approximation is taken.

According to the above, the conservation equations can be established as follows:

• The equation of conservation of mass for an incompressible fluid is obtained from Eq. (7) and remains as

  $$\nabla ·\overrightarrow{\mathrm{v}}=0$$E9

• The momentum conservation equation with the convection term is

  $$\frac{\partial \overrightarrow{\mathrm{v}}}{\mathrm{\partial t}}+\overrightarrow{\mathrm{v}}·\left(\nabla \overrightarrow{\mathrm{v}}\right)=-\frac{1}{{\mathrm{\rho }}_0}\nabla \left(\mathrm{P}-{\mathrm{\rho }}_0\overrightarrow{\mathrm{g}}\right)+{\nabla }^2\overrightarrow{\mathrm{v}}+\overline{\mathrm{g}}\mathrm{\beta }\left(\mathrm{T}-{\mathrm{T}}_0\right)$$E10

In the above equation, ρ is the initial density of the fluid, β is the coefficient of thermal expansion, $\overrightarrow{\mathrm{g}}$is the gravity vector (which only points in the y direction), t is time, $\mathrm{\nu }$is the kinematic viscosity, and T₀ is the external temperature.

• Equation of conservation of energy

$$\frac{\mathrm{\partial T}}{\mathrm{\partial t}}+\overrightarrow{\mathrm{v}}·\left(\nabla \mathrm{T}\right)=\mathrm{\alpha }{\nabla }^2\mathrm{T}$$E11

In the above equation, the relationship $\raisebox{1ex}{$\mathrm{k}$}\!\left/ \!\raisebox{-1ex}{$\mathrm{\rho }{\mathrm{C}}_{\mathrm{p}}$}\right.$ that appears in Eq. (3) represents the diffusion coefficient α.

For the system being studied, it is considered that the fluid is static and with external temperature when t = 0. Similarly, the system is considered only in two dimensions, therefore $\overrightarrow{\mathrm{v}}=\left(u,v,0\right)$, where “u” and “v” are the velocity components in the “x” and “y” direction respectively and are functions of time.

With the above, they were established as boundary conditions

$$\mathrm{u}\left(\mathrm{t},\mathrm{x},0\right)=0,\mathrm{u}\left(\mathrm{t},0,\mathrm{y}\right)=0,\mathrm{u}\left(\mathrm{t},\mathrm{L},\mathrm{y}\right)=0,\frac{\mathrm{\partial u}}{\mathrm{\partial y}}\left(\mathrm{t},\mathrm{x},\mathrm{H}\right)=0$$E12

$$\mathrm{v}\left(\mathrm{t},\mathrm{x},0\right)=0,\mathrm{v}\left(\mathrm{t},0,\mathrm{y}\right)=0,\mathrm{v}\left(\mathrm{t},\mathrm{L},\mathrm{y}\right)=0,\mathrm{v}\left(\mathrm{t},\mathrm{x},\mathrm{H}\right)=0$$E13

For the temperature we have

$$\mathrm{T}\left(\mathrm{t},\mathrm{x},0\right)={\mathrm{T}}_0,\mathrm{T}\left(\mathrm{t},0,\mathrm{y}\right)={\mathrm{T}}_0,\mathrm{T}\left(\mathrm{t},\mathrm{L},\mathrm{y}\right)={\mathrm{T}}_0,\mathrm{T}\left(\mathrm{t},\mathrm{x},\mathrm{H}\right)={\mathrm{T}}_0$$E14

Except for the points where the heat sources are located, as already mentioned above.

The previous system was expressed in dimensionless form considering the following characteristic variables

$${\mathrm{x}}^{\ast }=\frac{\mathrm{x}}{\mathrm{H}},{\mathrm{y}}^{\ast }=\frac{\mathrm{y}}{\mathrm{H}},\mathrm{T}=\frac{\mathrm{T}-{\mathrm{T}}_0}{{\mathrm{T}}_{\mathrm{H}}-{\mathrm{T}}_{\mathrm{o}}},{\mathrm{t}}^{\ast }=\frac{\mathrm{\alpha }}{{\mathrm{H}}^2}\mathrm{t}$$E15

$${\mathrm{p}}^{\ast }=\frac{{\mathrm{H}}^2}{{\mathrm{\rho }}_0{\mathrm{\alpha }}^2}\mathrm{p},{\mathrm{u}}^{\ast }=\frac{\mathrm{H}}{\mathrm{\alpha }}\mathrm{u},{\mathrm{v}}^2=\frac{\mathrm{H}}{\mathrm{\alpha }}\mathrm{v}$$E16

The constant T_H is defined from the heat source as

$${\mathrm{T}}_{\mathrm{H}}=\frac{\mathrm{qH}}{\mathrm{K}}+{\mathrm{T}}_0$$E17

By introducing the previous dimensionless variables to the system, we obtain

$$\frac{\partial {\mathrm{u}}^{\ast }}{\partial {\mathrm{x}}^{\ast }}+\frac{\partial {\mathrm{v}}^{\ast }}{\partial {\mathrm{y}}^{\ast }}=0$$E18

$$\frac{\mathrm{D}{\mathrm{u}}^{\ast }}{\mathrm{D}{\mathrm{t}}^{\ast }}=-\frac{\partial {\mathrm{P}}^{\ast }}{\partial {\mathrm{x}}^{\ast }}+Pr\left(\frac{{\partial }^2}{\partial {\mathrm{x}}^{\ast 2}}{\mathrm{u}}^{\ast }+\frac{{\partial }^2}{\partial {\mathrm{y}}^{\ast 2}}{\mathrm{v}}^{\ast }\right)$$E19

$$\frac{\mathrm{D}{\mathrm{v}}^{\ast }}{\mathrm{D}{\mathrm{t}}^{\ast }}=-\frac{\partial {\mathrm{P}}^{\ast }}{\partial {\mathrm{y}}^{\ast }}+Pr\left(\frac{{\partial }^2}{\partial {\mathrm{x}}^{\ast 2}}{\mathrm{u}}^{\ast }+\frac{{\partial }^2}{\partial {\mathrm{y}}^{\ast 2}}{\mathrm{v}}^{\ast }\right)+\text{RaPr}{\mathrm{T}}^{\ast }$$E20

$$\frac{\mathrm{D}{\mathrm{T}}^{\ast }}{\mathrm{D}{\mathrm{t}}^{\ast }}=\frac{{\partial }^2{\mathrm{T}}^{\ast }}{\partial {\mathrm{x}}^{\ast 2}}+\frac{{\partial }^2{\mathrm{T}}^{\ast }}{\partial {\mathrm{y}}^{\ast 2}}$$E21

where Pr is the Prandtl number and Ra is the Rayleigh number, which are defined as

$$Pr=\frac{\mathrm{\nu }}{\mathrm{\alpha }},\mathrm{Ra}=\frac{\mathrm{\beta g}{\mathrm{L}}^3\left({\mathrm{T}}_{\mathrm{H}}-{\mathrm{T}}_0\right)}{\mathrm{\nu \alpha }}$$E22

In Eq. (22), the operator $\raisebox{1ex}{$\mathrm{D}$}\!\left/ \!\raisebox{-1ex}{$\mathrm{Dt}$}\right.$represents the material derivative that is defined as $\frac{\mathrm{DA}}{\mathrm{Dt}}=\frac{\mathrm{\partial A}}{\mathrm{\partial t}}+\overrightarrow{\mathrm{v}}·\left(\mathrm{\nabla }\mathrm{A}\right)$. The term P* is defined as ${\mathrm{P}}^{\ast }={\mathrm{p}}^{\ast }+\mathrm{g}{\mathrm{y}}^{\ast }{\mathrm{L}}^{\mathrm{3}}/{\mathrm{\alpha }}^{\mathrm{2}}$.

For simplicity, from now on the use of (*) to indicate dimensionless variables will be omitted, assuming that all equations are in dimensionless form.

In order to reduce the system, the stream function formulation, φ, and vorticity, ω, are applied, which are defined as [7]

$$\mathrm{u}=\frac{\mathrm{\partial \phi }}{\mathrm{\partial y}},\mathrm{v}=-\frac{\mathrm{\partial \phi }}{\mathrm{\partial x}},\mathrm{\omega }=-\frac{\mathrm{\partial u}}{\mathrm{\partial y}}+\frac{\mathrm{\partial v}}{\mathrm{\partial x}}$$E23

so the system is rewritten as follows

$$\frac{{\partial }^2\mathrm{\phi }}{\partial {\mathrm{x}}^2}+\frac{{\partial }^2\mathrm{\phi }}{\partial {\mathrm{y}}^2}=-\mathrm{\omega }$$E24

$$\frac{\mathrm{D\omega }}{\mathrm{Dt}}=Pr\left(\frac{{\partial }^2\mathrm{\omega }}{\partial {\mathrm{x}}^2}+\frac{{\partial }^2\mathrm{\omega }}{\partial {\mathrm{y}}^2}\right)+\text{RaPr}\frac{\mathrm{\partial T}}{\mathrm{\partial x}}$$E25

$$\frac{\mathrm{DT}}{\mathrm{Dt}}=\frac{{\partial }^2\mathrm{T}}{\partial {\mathrm{x}}^2}+\frac{{\partial }^2\mathrm{T}}{\partial {\mathrm{y}}^2}$$E26

With this new formulation, the new initial conditions are

$$\mathrm{\phi }\left(0,\mathrm{x},\mathrm{y}\right)=0,\mathrm{\omega }\left(0,\mathrm{x},\mathrm{y}\right)=0,\mathrm{T}\left(0,\mathrm{x},\mathrm{y}\right)=0$$E27

and the boundary conditions are as follows

$$\mathrm{\phi }\left(\mathrm{t},0,\mathrm{y}\right)=0,\mathrm{\phi }\left(\mathrm{t},1,\mathrm{y}\right)=0,\mathrm{\phi }\left(\mathrm{t},\mathrm{x},0\right)=0,\mathrm{\phi }\left(\mathrm{t},\mathrm{x},1\right)=0$$E28

$$\mathrm{\omega }\left(\mathrm{t},0,\mathrm{y}\right)=0,\mathrm{\omega }\left(\mathrm{t},1,\mathrm{y}\right)=0,\mathrm{\omega }\left(\mathrm{t},\mathrm{x},0\right)=0,\mathrm{\omega }\left(\mathrm{t},\mathrm{x},1\right)=0$$E29

$$\mathrm{T}\left(\mathrm{t},0,\mathrm{y}\right)=0,\mathrm{T}\left(\mathrm{t},1,\mathrm{y}\right)=0,\mathrm{T}\left(\mathrm{t},\mathrm{x},0\right)=0,\mathrm{T}\left(\mathrm{t},\mathrm{x},1\right)=0$$E30

and for hot zones [8, 9, 10, 11, 12, 13, 14]

$$\frac{\mathrm{\partial T}}{\mathrm{\partial y}}=1$$E31

2.4 Numerical method

To solve the system, the Method of Alternate Directions Implicit (ADI) was used [15]. Regarding the discretization of the equations, the finite difference scheme is used. For the time derivative that appears on the left-hand side of Eq. (7), we have

$$\frac{\mathrm{\partial T}}{\mathrm{\partial t}}=\frac{{\mathrm{T}}_{\mathrm{i},\mathrm{j}}-{\mathrm{T}}_{\mathrm{i},\mathrm{j}}^{\mathrm{n}}}{\mathrm{\Delta t}}$$E32

In Eq. (32), T_i,j is the value of T at point (i, j) at the present instant, while ${\mathrm{T}}_{\mathrm{i},\mathrm{j}}^{\mathrm{n}}$ is the value of T at point (i, j) at the previous time instant. For the case of the spatial derivatives found in the first term on the right-hand side of (Eq. (6)), we have

$$\frac{\mathrm{\partial T}}{\mathrm{\partial x}}=\frac{{\mathrm{T}}_{\mathrm{i}+1,\mathrm{j}}-{\mathrm{T}}_{\mathrm{i}-1,\mathrm{j}}}{2\mathrm{\Delta x}}$$E33

$$\frac{\mathrm{\partial T}}{\mathrm{\partial y}}=\frac{{\mathrm{T}}_{\mathrm{i},\mathrm{j}+1}-{\mathrm{T}}_{\mathrm{i},\mathrm{j}-1}}{2\mathrm{\Delta y}}$$E34

$$\frac{{\partial }^2\mathrm{T}}{\partial {\mathrm{x}}^2}=\frac{{\mathrm{T}}_{\mathrm{i}+1,\mathrm{j}}-2{\mathrm{T}}_{\mathrm{i},\mathrm{j}}+{\mathrm{T}}_{\mathrm{i}-1,\mathrm{j}}}{\mathrm{\Delta }{\mathrm{x}}^2}$$E35

$$\frac{{\partial }^2\mathrm{T}}{\partial {\mathrm{y}}^2}=\frac{{\mathrm{T}}_{\mathrm{i},\mathrm{j}+1}-2{\mathrm{T}}_{\mathrm{i},\mathrm{j}}+{\mathrm{T}}_{\mathrm{i},\mathrm{j}-1}}{\mathrm{\Delta }{\mathrm{y}}^2}$$E36

The numerical model was solved using a 200 × 100 uniform mesh and a time step $\mathrm{\Delta t}=1\times {10}^{-5}$. As a convergence criterion, it was established that $\mathrm{\Delta \phi },\mathrm{\Delta \omega },\mathrm{\Delta }\mathrm{T}\le 1\times {10}^{-6}$. The code was developed in the FORTRAN 90 programming language. It was compiled and executed using the commercial package Microsoft Visual Studio, using a 64-bit HP Pavilion 15-cw1xxx laptop, with an AMD Ryzen 53,500 U processor with a Radeon Vega video card. Mobile Gfx. [16, 17, 18, 19].

The general algorithm is expressed as follows (Figure 5):

Figure 5.

General algorithm for the Method of Alternate Directions Implicit.

The subroutines used to solve Eqs. (28)-(30) are given in Tables 1–3.

Subroutine for $\mathit{\phi }\mathit{:}$ !----------------------------------; SUBROUTINE funciondecorriente (y,dx,dy,dx2,dy2,phi,dphi,omega,u,v); USE Constantes; IMPLICIT NONE; REAL (kind=8), DIMENSION(mi) :: dx,dx2; REAL (kind=8), DIMENSION(nj) :: y,dy,dy2; REAL (kind=8), DIMENSION(mi) :: ai,bi,ci,ri,ui; REAL (kind=8), DIMENSION(nj) :: aj,bj,cj,rj,uj; REAL (kind=8), DIMENSION(mi,nj) :: v,u,res,omega,chi; REAL (kind=8), DIMENSION(mi,nj) :: phi,dphi,phiant; DO i= 2,mi-1; DO j= 2,nj-1; u(i,j) = (phi(i,j+1)-phi(i,j-1))/2.d0*dy(j); v(i,j) = -(phi(i+1,j)-phi(i-1,j))/2.d0*dx(i); res(i,j)= dti*((phi(i+1,j)-2.d0*phi(i,j)+phi(i-1,j))/dx2(i) &; & + ((phi(i,j+1)-2.d0*phi(i,j)+phi(i,j-1))/dy2(j)) &; & + omega(i,j)); END DO; END DO; DO j= 1,nj

bi(1) = 1.d0; ci(1) = 0.d0; ri(1) = 0.d0 - phi(1,j); ai(mi) = 0.d0; bi(mi) = 1.d0; ri(mi) = 0.d0 - phi(mi,j); DO i= 2,mi-1; ai(i) = -dti/dx2(i); bi(i) = 1.d0 + (2.d0*dti/dx2(i)); ci(i) = -dti/dx2(i); ri(i) = res(i,j); END DO; CALL tri(ai,bi,ci,ri,ui,mi); DO i= 1,mi; chi(i,j)=ui(i); END DO; END DO; DO i= 1,mi; bj(1) = 1.d0; cj(1) = 0.d0

rj(1) = 0.d0 - phi(i,1); aj(nj) = 0.d0; bj(nj) = 1.d0; rj(nj) = 0.d0 - phi(i,nj); DO j= 2,nj-1; aj(j) = -dti/dy2(j); bj(j) = 1.d0 + (2.d0*dti/dy2(j)); cj(j) = -dti/dy2(j); rj(j) = chi(i,j); END DO; CALL tri(aj,bj,cj,rj,uj,nj); DO j= 1,nj; dphi(i,j) = uj(j); END DO; END DO; DO i= 1,mi; DO j= 1,nj; phi(i,j)= phi(i,j) + dphi(i,j); END DO; END DO; END SUBROUTINE funciondecorriente; !--------------------

Table 1.

Subroutine for φ.

Subroutine for $\mathit{\omega }$: !----------------------------------; SUBROUTINE vorticidad (dx,dy,dx2,dy2,phi,u,v,theta,omega,domega); USE Constantes; IMPLICIT NONE; REAL (kind=8), DIMENSION(mi) :: dx,dx2; REAL (kind=8), DIMENSION(nj) :: dy,dy2; REAL (kind=8), DIMENSION(mi) :: ai,bi,ci,ri,ui; REAL (kind=8), DIMENSION(nj) :: aj,bj,cj,rj,uj; REAL (kind=8), DIMENSION(mi,nj) :: v,u,res,phi,chi1,theta; REAL (kind=8), DIMENSION(mi,nj) :: omega,domega,omegant; DO i= 2,mi-1; DO j= 2,nj-1; res(i,j) = -(dti*u(i,j)*((omega (i+1,j)-omega (i+1,j))/(2.d0*dx(i)))) &; & -(dti*v(i,j)*((omega (i,j+1)-omega (i,j-1))/(2.d0*dy(j)))) &; & +((dti*Pr)*((omega(i+1,j)-2.d0*omega(i,j)+omega(i-1,j))/(dx2(i)))) &; & +((dti*Pr)*((omega(i,j+1)-2.d0*omega(i,j)+omega(i,j-1))/(dy2(j))))&; & -(dti*Ra*Pr*((theta(i+1,j)-theta(i-1,j))/(2.d0*dx(i)))); END DO; END DO

DO j= 1,nj; bi(1) = 1.d0; ci(1) = 0.d0; ri(1) = 2*(phi(1,j)-phi(2,j))/dx2(2) - omega(1,j); ai(mi) = 0.d0; bi(mi) = 1.d0; ri(mi) = (2.d0*(phi(mi,j)-phi(mi-1,j))/dx2(mi)) - omega(mi,j); DO i= 2,mi-1; ai(i) = dti*(u(i,j)/(2.d0*dx(i))) - dti*Pr/(dx2(i)); bi(i) = 1.d0 + (dti*2.d0*Pr)/(dx2(i)); ci(i) = -dti*(u(i,j)/(2.d0*dx(i))) - dti*Pr/(dx2(i)); ri(i) = res(i,j); END DO; CALL tri(ai,bi,ci,ri,ui,mi); DO i= 1,mi; chi1(i,j)=ui(i); END DO; END DO; DO i= 1,mi; bj(1) = 1.d0; cj(1) = 0.d0; rj(1) = 2.d0*(phi(i,1)-phi(i,2))/dy2(2) - omega(i,1); aj(nj) = 0.d0; bj(nj) = 1.d0; rj(nj) = 0 - omega(i,nj); DO j= 2,nj-1

aj(j) = dti*(v(i,j)/(2.d0*dy(j))) &; & -dti*Pr/(dy2(j)); bj(j) = 1.d0 + ((dti*2.d0*Pr)/(dy2(j))); cj(j) = -dti*(v(i,j)/(2.d0*dy(j))) &; & -dti*Pr/(dy2(j)); rj(j) = chi1(i,j); END DO; CALL tri(aj,bj,cj,rj,uj,nj); DO j= 1,nj; domega(i,j) = uj(j); END DO; END DO; DO i= 1,mi; DO j= 1,nj; omega(i,j)= omega(i,j) + domega(i,j); END DO; END DO; END SUBROUTINE vorticidad; !------------------

Table 2.

Subroutine for ω.

Subroutine for $\mathit{T}$: !---------------------------; SUBROUTINE energia (dx,dy,dx2,dy2,u,v,theta,dtheta); USE constantes; IMPLICIT NONE; REAL (kind=8), DIMENSION(mi) :: dx,dx2; REAL (kind=8), DIMENSION(nj) :: dy,dy2; REAL (kind=8), DIMENSION(mi) :: ai,bi,ci,ri,ui; REAL (kind=8), DIMENSION(nj) :: aj,bj,cj,rj,uj; REAL (kind=8), DIMENSION(mi,nj) :: v,u,res,chi2; REAL (kind=8), DIMENSION(mi,nj) :: theta,dtheta,thetant; DO i= 2,mi-1; DO j= 2,nj-1; res(i,j) = -dti*(u(i,j)*((theta(i+1,j)-theta(i-1,j))/(2.d0*dx(i)))) &; & -dti*(v(i,j)*((theta(i,j+1)-theta(i,j-1))/(2.d0*dy(j)))) &; & +(dti*(((theta(i+1,j)-2.d0*theta(i,j)+theta(i-1,j))/(dx2(i))) &; & +((theta(i,j+1)-2.d0*theta(i,j)+theta(i,j-1))/(dy2(j))))); END DO; END DO; DO j= 1,nj; bi(1) = 1.d0; ci(1) = 0.d0; ri(1) = 0.d0 - theta(1,j); ai(mi) = 0.d0; bi(mi) = 1.d0; ri(mi) = 0.d0 - theta(mi,j); DO i= 2,mi-1; ai(i) = dti*(u(i,j)/(2.d0*dx(i)))-(dti/(dx2(i))); bi(i) = 1+((dti*2.d0)/(dx2(i))); ci(i) = -dti*(u(i,j)/(2.d0*dx(i)))-(dti/(dx2(i))); ri(i) = res(i,j); END DO; CALL tri(ai,bi,ci,ri,ui,mi); DO i= 1,mi; chi2(i,j)=ui(i); END DO; END DO; DO i= 1,mi

IF(configura==1) THEN; IF(i<(mi/3)) THEN; bj(1) = 1.d0; cj(1) = 0.d0; rj(1) = 0.d0 - theta(i,1); ELSE IF(i>=(mi/3).and.i<=(4*mi/9)) THEN; bj(1) = 1.d0; cj(1) = -1.d0; rj(1) = 1.d0 - (theta(i,1)-theta(i,2)); ELSE IF (i>4*mi/9.and.i<(5*mi/9))THEN; bj(1) = 1.d0; cj(1) = 0.d0; rj(1) = 0.d0 - theta(i,1); ELSE IF (i>=5*mi/9.and.i<=(2*mi/3))THEN; bj(1) = 1.d0; cj(1) = -1.d0; rj(1) = 1.d0 - (theta(i,1)-theta(i,2)); ELSE IF (i>2*mi/3)THEN; bj(1) = 1.d0; cj(1) = 0.d0; rj(1) = 0.d0 - theta(i,1); END IF; END IF; IF(configura==2) THEN; IF(i<(mi/9)) THEN; bj(1) = 1.d0; cj(1) = 0.d0; rj(1) = 0.d0 - theta(i,1); ELSE IF(i>=(mi/9).and.i<=(2*mi/9)) THEN; bj(1) = 1.d0; cj(1) = -1.d0; rj(1) = 1.d0 - (theta(i,1)-theta(i,2)); ELSE IF (i>2*mi/9.and.i<(3*mi/9))THEN; bj(1) = 1.d0; cj(1) = 0.d0; rj(1) = 0.d0 - theta(i,1); ELSE IF (i>=mi/3.and.i<=(4*mi/9))THEN; bj(1) = 1.d0; cj(1) = -1.d0

rj(1) = 1.d0 - (theta(i,1)-theta(i,2)); ELSE IF (i>4*mi/9.and.i<(5*mi/9))THEN; bj(1) = 1.d0; cj(1) = 0.d0; rj(1) = 0.d0 - theta(i,1); ELSE IF (i>=5*mi/9.and.i<=(2*mi/3))THEN; bj(1) = 1.d0; cj(1) = -1.d0; rj(1) = 1.d0 - (theta(i,1)-theta(i,2)); ELSE IF (i>2*mi/3.and.i<(7*mi/9))THEN; bj(1) = 1.d0; cj(1) = 0.d0; rj(1) = 0.d0 - theta(i,1); ELSE IF (i>=7*mi/9.and.i<=(8*mi/9))THEN; bj(1) = 1.d0; cj(1) = -1.d0; rj(1) = 1.d0 - (theta(i,1)-theta(i,2)); ELSE IF (i>8*mi/9)THEN; bj(1) = 1.d0; cj(1) = 0.d0; rj(1) = 0.d0 - theta(i,1); END IF; END IF; aj(nj) = 0.d0; bj(nj) = 1.d0; rj(nj) = 0.d0 - theta(i,nj); DO j= 2,nj-1; aj(j) = dti*(v(i,j)/(2.d0*dy(j)))-(dti/(dy2(j))); bj(j) = 1+((dti*2.d0)/(dy2(j))); cj(j) = -dti*(v(i,j)/(2.d0*dy(j)))-(dti/(dy2(j))); rj(j) = chi2(i,j); END DO; CALL tri(aj,bj,cj,rj,uj,nj); DO j= 1,nj; dtheta(i,j) = uj(j); END DO; END DO; DO i= 1,mi; DO j= 1,nj; theta(i,j)= theta(i,j) + dtheta(i,j); END DO; END DO; END SUBROUTINE energia; !----------------------

Table 3.

Subroutine for T.

3. Simulations and results

Before carrying out the simulations, the code used was validated using the fields of temperature and stream function obtained by Davis G., 1983.

As shown in the Figures 6–8 corresponding to the profiles for both φ, ω, and T obtained, they replicate the results reported in the literature.

Figure 6.

Stream function φ for the proposed arrangements with different Rayleigh numbers. (a) Two heat sources and Ra = 1 × 10², (b) four heat sources and Ra = 1 × 10², (c) two heat sources and Ra = 1 × 10³, (d) four heat sources and Ra = 1 × 10³, (e) two heat sources and Ra = 1 × 10⁴, (f) four heat sources and Ra = 1 × 10⁴, (g) two heat sources and Ra = 1 × 10⁵, (h) four heat sources and Ra = 1 × 10⁵.

Figure 7.

Vorticity ω for the proposed arrangements with different Rayleigh numbers. (a) Two heat sources and Ra = 1 × 10², (b) four heat sources and Ra = 1 × 10², (c) two heat sources and Ra = 1 × 10³, (d) four heat sources and Ra = 1 × 10³, (e) two heat sources and Ra = 1 × 10⁴, (f) four heat sources and Ra = 1 × 10⁴, (g) two heat sources and Ra = 1 × 10⁵, (h) four heat sources and Ra = 1 × 10⁵.

Figure 8.

Temperature T for the proposed arrangements with different Rayleigh numbers. (a) Two heat sources and Ra = 1 × 10², (b) four heat sources and Ra = 1 × 10², (c) two heat sources and Ra = 1 × 10³, (d) four heat sources and Ra = 1 × 10³, (e) two heat sources and Ra = 1 × 10⁴, (f) four heat sources and Ra = 1 × 10⁴, (g) two heat sources and Ra = 1 × 10⁵, (h) four heat sources and Ra = 1 × 10⁵.

All the simulations carried out were obtained with a fixed value of Pr = 7, with the characterization of the fluid as water. The computational experiments were performed with the two heat source configurations, using four values of the Ra number Ra (1 × 10², 1 × 10³, 1 × 10⁴, 1 × 10⁵) obtaining the temperature distributions, as well as the formation of vorticity and the changes in the current lines.

Figures 6–8 show that, by keeping the fluid constant, with the greater number of Ra there is a greater T_H, i.e., the magnitude of q increases. The results obtained were found by establishing the following values for the constants $\mathrm{H}=0.2\mathrm{m},\mathrm{\beta }=207\times {10}^{-6}{\mathrm{K}}^{-1},\mathrm{k}=0.58\mathrm{W}/\mathrm{mK}\mathrm{\alpha }=1.3882\times {10}^{-4}{\mathrm{m}}^2/\mathrm{s},\mathrm{\nu }=1.004\times {10}^{-6}{\mathrm{m}}^2/\mathrm{s},{\mathrm{T}}_0=25°\mathrm{C}.$

Although the temperature profile does not express large changes for different values of Ra, the magnitude of the temperatures present in the system do.

The results of the simulations to reach the steady state times in the water heating process in the shrink tank show similar trends with the two arrangements of heat sources and very similar time values. For both arrangements, as the Rayleigh number increases, the times increase, but not significantly, having dimensionless time values that oscillate between 0.8647 and 0.9101 for the two and four arrangements respectively for the case of the Rayleigh value of 1 × 10². Y of 1.5503 and 1.6171 for the two and four arrangements respectively for Rayleigh = 1 × 10⁵, as can be seen in Figure 9.

Figure 9.

Steady state time for the Rayleigh numbers of 1 × 10², 1 × 10³, 1 × 10⁴, and 1 × 10⁵ for the two heat sources.

Regarding the temperature of the surface of the heat source, in the same way, similar trends are observed with the two arrangements of heat sources and very similar temperature values. But for Rayleigh values above 1 × 10⁴, a difference is observed between both arrangements as can be seen in Figure 10. Obtaining values that go from 25 to 26°C for Rayleigh numbers = 1 × 10² and 1 × 10³ respectively for both arrangements. However, for values of Rayleigh = 1 × 10⁵ there are temperatures of 176.95 and 193.84°C for arrangements of two and four resistors respectively.

Figure 10.

Tank surface temperature for the Rayleigh numbers of 1 × 10², 1 × 10³, 1 × 10⁴, and 1 × 10⁵ for the two heat sources.

In the same way, average temperature values of the water inside the shrink tank were obtained. There are similar trends with the two arrangements of heat sources and very similar temperature values below Rayleigh numbers = 1 × 10³. But for values of Rayleigh ≥1 × 10⁵ there are more marked temperature differences, as can be seen in Figure 11, for both arrays, but higher temperature values for the case of four arrays.

Figure 11.

Average tank temperature for the Rayleigh numbers of 1 × 10², 1 × 10³, 1 × 10⁴, and 1 × 10⁵ for the two heat sources.

For the case of Ra = 1 × 10², which implies a heat source of q = 0.0024 W, there is a temperature near the heat sources of 25.1519°C on average, above ambient temperature. On the other hand, with a value of Ra = 1 × 10⁴ with a heat source of q = 24.48 W, it produces a temperature of 176.9534°C, which is above the ambient temperature, however, in the areas closest to the heat sources this temperature can easily be exceeded.

From the above, we can say that, if we want an adequate equation through this study, it can be established that the heat sources, considering isothermal walls, must produce at least 25 W of heat, with the consideration that said sources, are in direct contact with the fluid.

On the other hand, the need to place four heat sources is contemplated, since, with the arrangement of two, it is insufficient to achieve a homogenization of the temperature in the heat shrink tank, generating areas where the temperature gradients are very small.

4. Conclusions

In the present work, computer simulations were carried out to evaluate the heat transfer in a shrink tank built in the facilities of the Faculty of Higher Studies Cuautitlán—UNAM.

For the heating of the water in the tank, two different arrangements of electrical resistances were implemented as heat sources.

In order to quantify the heating times and temperature distributions, the Method of Alternate Directions Implicit was used in combination with the finite difference scheme, obtaining profiles of the stream function φ and vorticity ω, which helped with the selection of the resistance arrangement that guarantees a better heat transfer of the water in the tank.

From the simulations carried out in this work, it was observed that the formation of convective cells favors the homogenization of the temperature in the tank and that increasing the value of the Rayleigh number increases the vorticity but not the temperature field, which allows keep lower power.

It is confirmed that, when working with non-isothermal sources, a higher energy accumulation is obtained, but a less homogeneous temperature field than that produced by isothermal sources, deeper studies on this were carried out by Ostrach in 1988 [9].

The results of the simulations to reach the steady state times in the water heating process in the shrink tank, showed that the type of arrangement does not interfere directly in said time. Although it is observed that an increase in the Rayleigh number brings as a consequence an increase in the times to reach the steady state.

Likewise, regarding the temperature of the surface of the heat source, it can be seen that for values of the Rayleigh number above 1 × 10³, the temperatures increase and there is a considerable difference with Ra = 1 × 10⁵ for both arrangements, but higher for a four heat source arrangement.

Regarding the average temperature values of the water inside the shrink tank, for values of Rayleigh ≥1 × 10⁵ there are more marked temperature differences between both arrangements, but higher temperature values are presented for the case of four arrangements, favoring the conditions required in the heat shrink process.

From the results obtained from the heating times with their respective temperature distributions, we can conclude that the arrangement that best optimizes the heat transport process in the shrink tank is the one corresponding to the four resistance arrangements, achieving a homogeneous temperature of 87°C, in times less than 9 min with a heat flux of q = 24.48 W.

The implementation of solid walls is considered for future work to study the effects of different insulation, in order to better conserve heat within the system, as well as the inclusion of the evaluation of heat transport in the shrink tub with the four arrangements, but incorporating a piece of meat with its shrink wrap, in order to quantify the distribution of temperatures with their respective validation with microbiological methods that indicate the null contamination of the product with bacteria.

Acknowledgments

The present work was developed under the sponsorship of the Facultad de Estudios Superiores Cuautitlán—Universidad Nacional Autónoma de México.

Conflict of interest

The authors declare no conflicts of interest regarding the publication of this paper.