CHEMICAL ENGINEERING TRANSACTIONS  
 

VOL. 61, 2017 

A publication of 

 

The Italian Association 
of Chemical Engineering 
Online at www.aidic.it/cet 

Guest Editors: Petar S Varbanov, Rongxin Su, Hon Loong Lam, Xia Liu, Jiří J Klemeš 
Copyright © 2017, AIDIC Servizi S.r.l. 

ISBN 978-88-95608-51-8; ISSN 2283-9216 

Computational Fluid Dynamic Analysis of Cooling in a Mixed 

Vessel using a Non-Newtonian Medium 

Zsolt Harsfalvi*, Christian Jordan, Bahram Haddadi, Michael Harasek 

Technische Universität Wien, Institute of Chemical Engineering, Getreidemarkt 9/166, A-1060 Vienna, Austria 

zsolt.harsfalvi@tuwien.ac.at 

In the food industry, the efficient production of high quality products leads to increased investigation efforts in 

the production process. High shear rates, improper heat transfer during the production process can lead to low 

quality products. In many cases the produced food preparation contains additives e.g. pectin and starch, which 

are affecting the products flow behaviour especially at lower temperatures. High shear rates during the 

production can be observed rather during the cooling and mixing processes of these products. The goal of this 

paper is to develop a CFD solver in the open source software OpenFOAM® which will be able to recognise the 

weak spots of the production. For this purpose, a detailed analysis of the material properties was completed 

with the help of measurements. In the food industry, the average volume of a produced batch is up to 1,000 L. 

To speed up the development of the CFD code a model production process was set up in small scale using a 

10 L stainless steel vessel with a simple mixer blade and heating/cooling jacket. The temperature of the medium 

inside the vessel was measured at various points. The torque of the shaft of the rotor blade was also measured 

during the production time. Using these measurement results the CFD computation could be validated. For the 

computed medium a shear rate and temperature dependent viscosity and a temperature dependent heat 

capacity, thermal conductivity and density model was created. To calculate the mixing during the process, 

dynamic mesh using the AMI concept was applied. The industrial cooling processes takes up to 30 min, but 

considering the inefficient mixer blade and the weak cooling power the cooling time at small scale usually takes 

more time. Because of the high computational need and relatively small time steps the simulation was calculated 

only for 25 min real time. The data of the installed PT 100 temperature sensors and the measured torque showed 

a close agreement with the simulation results. 

1. Introduction 

The increasing social and economic importance of food production, in addition to more complex production 

technologies requires a more detailed investigation as further development in the quality of the food or in the 

production efficiency. Creating fruit preparation for yogurt production especially in large amounts can have lot 

of pitfalls, which can lead to either an inefficient, long production time or/and a lower level of food quality. 

Considering the fragile food materials and the high-quality demand of the final product the optimization of the 

production will need special tools. (Magerramov et al., 2006)  

In this paper, the development and validation of a CFD solver is presented for a pilot scale stirred cooling vessel. 

The model can later be applied for industrial cooling processes, e.g. equipped with complicated mixer blades 

and heat exchangers. The simulated process starts with a temperature of 90 °C which has to be cooled down 

to 35 °C within the shortest possible time to make it ready for further operations. An extremely challenging goal 

is to provide a suitable material property model for the CFD solver. The viscosity of the non-Newtonian medium 

(shear thinning) strongly increases with decreasing temperature. Properties as thermal conductivity and specific 

heat have also to be considered as a temperature function. The developed CFD solver can be used for the 

identification of zones of high shear rate and ineffective heat transfer. The cooling experiments for testing and 

validation reasons were accomplished in a cooled mixed vessel (macro-viscosimeter) developed at Technische 

Universität Wien (Pohn et al., 2010). 

                               
 
 

 

 
   

                                                  
DOI: 10.3303/CET1761092

 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 

Please cite this article as: Harsfalvi Z., Jordan C., Haddadi B., Harasek M., 2017, Computational fluid dynamic analysis of cooling in a mixed 
vessel using a non-newtonian medium, Chemical Engineering Transactions, 61, 565-570  DOI:10.3303/CET1761092  

565



2. Computational Fluid Dynamics (CFD) 

Computer based simulations allow a detailed analysis of fluid systems. This tool is known as Computational 

Fluid Dynamics or CFD and it is suitable for solving a wide range of industrial and non-industrial applications 

e.g. fluid flow, mass transfer, heat transfer, etc. (Fletcher, 1988). Using the 3D discretization of space, any 

computed phenomena can be investigated including even small details. Although CFD is a powerful tool, it is 

important, that the user has a proper knowledge regarding the investigated field. For validation reasons, it is 

helpful to accomplish lab scale measurements parallel to the computation to compare and develop the solver, 

to obtain reliable results. For the presented problem, the non-commercial open source software OpenFOAM® 

(version v3.0+) was applied. 

CFD is based on the Navier-Stokes and continuity equations Eq(1) and  Eq(2). These equations describe how 

the velocity, pressure, temperature and density of a moving fluid are related in time.  

𝜕𝒖

𝜕𝑡
+ (𝒖.∇)𝒖 = −

1

𝜌
∇𝑝 +

𝜇

𝜌
∇2𝒖 (1) 

𝜕𝜌

𝜕𝑡
+ ∇.(𝜌𝒖) = 0 (2) 

Where u [m/s] is the velocity, p [Pa] is the pressure, ρ is the density [kg/m³] and µ [Pas] is the dynamic viscosity. 
In the modelled system the energy storage and transport is calculated with the help of the energy equation 

Eq(3). 

𝜌 (
𝜕ℎ

𝜕𝑡
+ ∇.(ℎ𝒖)) = −

𝐷𝑝

𝐷𝑡
+ ∇.(𝐾∇T) + (�̅�.∇)𝒖 (3) 

h [J/kg] is enthalpy, K [W/(mK)] is the thermal conductivity of the fluid and τ is the shear stress. 
The investigated medium has a high Brix number [°Bx] because of the high pectin, starch and sugar content. 

These additives are affecting the viscosity. Therefore one of the most important part of Eq(1) is the diffusing 

viscous term and the viscosity itself. The energy transport during the cooling will be calculated using Eq(3), 

where the influence of the viscosity is also present in form of the shear stress. This term is rather small compared 

to the heat transfer and derives from the energy introduced by the rotation of the mixer blade.  

The vessel used in the cooling mixing process is usually not completely filled with the medium. The remaining 

volume of the vessel is filled with air. This will be taken into account in the simulation by using the VOF method 

and solving the transport equation. 

𝜕𝐶

𝜕𝑡
+ 𝒖.∇𝐶 = 0 (4) 

where C [-] is the fraction function with the interval 0 ≤ 𝐶 ≤ 1. The typical Reynolds number of the flow is smaller 
than 1, because of the high viscosity and low rotational speed. The flow characteristic is laminar, therefore no 

turbulence model was applied.  

2.1 Material properties model 
Due to the non-Newtonian behavior and temperature dependent thermophysical properties a complex property 

model was needed. The temperature dependence of heat capacity, thermal conductivity and density were 

modelled by the consideration of the Brix number (46,1 °Brix) of the medium using the experimental results from 

Carbal et al. (2007). 

  

Figure 1: Specific heat at 46.1 °Brix (a)  Calorimetry measurement, water/medium (b) 

2,700

2,800

2,900

3,000

3,100

3,200

270 290 310 330 350 370

H
e
a

t 
C

a
p

a
c
it
y
 [
J
/k

g
/K

]

Temperature [K]

Polyn. reg.

Carbal

Measurement
292

293

294

295

296

297

298

50 55 60 65 70 75 80

T
e

m
p

e
ra

tu
re

 [
K

]

Time [s]

a b 

566



  

Figure 2: (a) Thermal conductivity at 46.1 °Brix  (b) Density at 46.1 °Brix  

The experimental results for 46.1 °Brix were chosen and polynomial fitting functions for specific heat  

(Figure 1a), thermal conductivity (Figure 2a) and density were used (Figure 2b). The values of the specific heat 

were verified at selected points with the help of calorimetry measurements (Figure 1b) and in the case of density 

with pycnometer measurements. The verification measurement points showed acceptable results, hence the 

experimental data fitted with polynomial regressions were applied in the material models. The uncertainties 

showed in the figures are derived from the inaccuracy of the used measurement devices and the measurement 

difficulty due to the handling of the high viscose medium (e.g. small air bubbles during density measurements, 

slow mixing and heat exchange during the mixing calorimeter measurements). The uncertainty of the data 

presented by Carbal et al. (2007) is not known. 

The viscosity of the shear thinning medium is not only influenced by the temperature (Macías et al., 2013), but 

also by the shear rate, which was implemented in the material model using Eq(5). In the case of constant 

temperature this can be described by the Ostwald de Waele relationship (Goodwin and Hughes, 2000).  

𝜇 = (𝐾0 + 𝐾1𝑇 + 𝐾2𝑇
2 + 𝐾3𝑇

3) (
𝜕𝑢

𝜕𝑦
)
𝑛−1

 (5) 

In the equation 𝐾𝑛  is the flow consistency index, n the flow behavior index and 𝜕𝑢 𝜕𝑦⁄  the shear rate. For 
creating the viscosity model first measurements were accomplished with an Anton-Paar rotational viscometer 

for various probes. Few probes showed a slight difference between the measured viscosity field, typically at the 

bottom and at the top of the vessel. For the modelling of the field only a single measurement was chosen from 

the viscosity fields showing similar measurement results. The measurement was completed by changing the 

temperature at different, constant shear rates. By this method the viscosity field of the medium was known. 

Using polynomial regression, the proper values for 𝐾𝑛  and n were identified. The created viscosity model 
showed an 𝑛 = 0.46 flow behaviour index. Representative results of the rotational viscometer measurements 

and the polynomial approach are shown for few shear rate isolines in Figure 3. 

 

Figure 3: Temperature and shear rate dependent viscosity model 

0.42

0.44

0.46

0.48

0.5

0.52

270 290 310 330 350 370

T
h

e
rm

a
l 
c
o

n
d

. 
[W

/(
m

K
)]

Temperature [K]

Polyn. reg.

Carbal

1,120

1,140

1,160

1,180

1,200

1,220

1,240

270 290 310 330 350 370

D
e

n
s
it
y
 [
k
g

/m
³]

Temperature [K]

Polyn. reg.

Carbal

Measurement

0

5

10

15

20

25

293 303 313 323 333 343 353 363

V
is

c
o

s
it
y
 [
P

a
s
]

Temperature [K]

5 [1/s]

10 [1/s]

25 [1/s]

45 [1/s]

85 [1/s]

Polyn. approx. 85 [1/s]

Polyn. approx. 5 [1/s]

Polyn. approx. 25 [1/s]

a b 

567



2.2 Geometry of the vessel and the computational mesh 
The simulation was first tested on a small mixed vessel, called macro-viscosimeter (Pohn et al., 2011). This 

device was initially developed to measure viscosity of mediums with large particles and filaments (Figure 4). It 

has a 10 L volume, a heating/cooling jacket and a simple mixer blade in the centre. The torque on the shaft of 

the mixer blade can be measured with high accuracy. This was considered as an ideal device to test and later 

validate the CFD model.  

To be able to calculate the mixing process in the vessel, moving mesh, provided by the Arbitrary Mesh Interface 

(AMI) concept was applied. The geometry was split up into three different regions as shown in Figure 4.  

Region 3 is a hollow cylinder domain with no rotation, region 2 is the rotating domain including the mixer blade. 

Region 1 behaves like region 2 and was generated with the purpose to coarsen the mesh in this volume.  

Region 1 includes during the simulation almost only the air phase but the investigation concentrates mainly on 

the non-Newtonian medium. The second reason of the coarser mesh in this region is the possibility of the rising 

hot air at the start of the cooling. This can result in higher velocities which will cause smaller time steps. For the 

presented geometry, simulations with different boundary conditions were started. In this publication only the 

case for a closed vessel was presented, where the convection effects are not present. 

The introduced regions are connected by the interfaces where the transport of the physical quantities can take 

place. The mesh was generated using GAMBIT (version 2.4.) and has approximately 400,000 hexahedral cells. 

 

Figure 4: Geometry and mesh structure of the macro-viscosimeter; left image taken from (Pohn et al., 2011) 

2.3 Boundary and initial conditions 
The cooling phase of the production process which is happening immediately after the pasteurization was 

investigated here. At the beginning of the cooling the medium has a nearly constant temperature of 360 K. This 

is the initial temperature of the medium in the simulation except of the top 2 cm layer which was set to a lower  

340 K temperature as well as the air at the top, because of the heat loss due to convection. During the 

pasteurization, the top of the vessel is open and in contact with ambient air. Right after the finish of the 

pasteurization an isolating cylindrical cap was placed on the top of the macro-viscosimeter to avoid further 

convection effects in the air phase. For this reason, the top of the vessel was closed in the simulation and 

considered as adiabatic. The temperature on the inner vessel wall was measured during the cooling process. 

Using the measurement results time dependent polynomial functions were created which are matching with the 

real wall temperatures for the cylinder and for the bottom of the vessel. This function is only time dependent and 

the boundary condition does not include the thermal conductivity in the wall. In the case of the mixer blades the 

effect of the conductivity at the steel blade is also neglected. Here the boundary condition is considered also as 

adiabatic.  

These conditions were chosen because the thermal conductivity of the medium is two orders of magnitude lower 

as the steel vessel or mixer. Another reason for the chosen boundary conditions is the low cooling power, which 

results in a small cooling water temperature difference at the inlet and at the outlet of the heating/cooling jacket. 

This allows the assumption of an isothermal wall at the vessel boundary. Since the bottom of the vessel differed 

from the wall temperature a modified different time dependent boundary condition was applied here.  

The rotational speed of the blade was set to 30 min-1. Inserting the parameters to the implemented polynomial 

boundary conditions the simulation case was fully defined. 

568



3. Measurement 

The goal of the measurement was to obtain as much information as possible about the cooling process to create 

a matching model for the CFD simulation. This happened by the installation of 8 pieces of PT 100 temperature 

sensors into the vessel to investigate the temperature field. The temperature sensors were installed at different 

wall positions and also in the medium at different heights and radius. An additional help to support the validation 

process was the torque measurement of the system during the cooling using a Viscopakt measurement device. 

4. Results and discussion 

Industrial cooling processes for similar mediums usually take up to 30 min. Considering the presented device, 

the power of the heat exchange and the mixing is less effective compared to industrial technology. (Scargiali et 

al., 2013) This means that the cooling process using the macro-viscosimeter will take more time. Therefore, 

only the first 25 min were simulated and validated with experimental results to limit computational demand. For 

the post processing, the open source program Paraview (version 5.2.0.) was used. Selected contour plots of 

the results are shown in Figure 5. The streamlines of the velocity field in the vessel show that the medium is 

sucked in at the top and at the bottom of the rotating blade and it is pushed out in radial direction at the middle 

of the blade. This makes it possible that the colder medium at the bottom and at the top reaches the core of the 

vessel. This progress is clearly visible in the temperature field. The heat exchange and the gradient near to the 

wall can be also seen as expected. In spite of the presented flow the mixing process is slow due to the high 

viscosity. The non-Newtonian rheology can be recognised in the viscosity field - the shear thinning in the vicinity 

of the mixer blades adversely affects the mixing behaviour. The viscosity field represents also the interface of 

the liquid and air phases at the sudden change of the viscosity values. 

The results of the measurements (M) and the comparison with the simulation (CFD) data are shown in  

Figure 6. The “R” and “H” letters are representing the position of the PT 100 temperature sensors in the vessel. 

“R” stands for the radius measured from the axis of the rotor, “H” stands for the height of the sensor measured 

from the bottom of the vessel. For validation purposes, the temperature and the measured torque on the rotor 

shaft was used. The torque derives completely from the resistance of the fluid which is transmitted to the rotor. 

The force and the total torque were calculated by integration on the blade surface using the shear rate 

distribution and the viscosity on the rotor blade.  

The simulation results for the temperature are in close agreement for the measured locations. The slight 

difference could derive from the imperfect time dependent boundary condition or/and the inhomogeneity of the 

cooled medium due to the ineffective mixing during the pasteurization. In the case of the torque results, the 

slope of the functions is similar, but the CFD results over-predict the measured torque values. This is explained 

by the inhomogeneity in the medium. The mixing of starch and pectin during the cooking phase is not completely 

perfect, especially at the top and at the bottom of the vessel.  

The measured viscosity relations could be modelled for a single sample well as showed previously in Figure 3. 

The uncertainty originates rather from the pasteurized and mixed mediums slight inhomogeneity due to the 

ineffective mixer blades.  

 

Figure 5: Results of the cooling simulation at t = 1,500 s for the temperature, velocity field (using streamlines) 

and viscosity field 

569



 
 

Figure 6: Comparison of the CFD and measurement results 

In industrial cases, the mixer blade shape and the cooling process are more complex and usually more effective 

which results in a more homogeneous medium and constant starch and pectin concentrations during the 

process. Therefore, at the simulation of the industrial processes a better agreement between model and physical 

measurements is expected. 

5. Conclusions 

A cooling process for a non-Newtonian (shear thinning) medium was investigated using measurements and 

CFD in small scale with the goal to obtain a reliable CFD solver which can be used in the future to analyze more 

complex vessels. Besides the implementation of the solver, the creation of a suitable material model with the 

proper parameters was very challenging. The method can be used only in the case of a well described, accurate 

material property model. Therefore, if other media have to be included, new material models have to be 

generated and validated. Another important factor is the homogeneity of the mixed product. In case of industrial 

coolers the mixing intensity is high and the cooled product can be described as homogeneous. In case of the 

lab setup (makro-viscometer) due to less effective mixing a small concentration inhomogeneity can be detected 

between the top and the bottom layer in the vessel. This concentration gradient can explain a lower model 

prediction performance which resulted in deviations of the torque calculation.  

The CFD model has been applied to an industrial cooling process (food prep) in order to observe the shear rate 

and heat transfer within the cooled medium. Since some of these products include fruit pieces dispensed in the 

homogenous medium it is considered to develop the model in this direction. 

References 

Carbal R.A.F., Orrego-Alzate C.E., Gabas A.L., Telis-Romero J., 2007, Rheological and Thermophysical 

Properties of Blackberry Juice, Food Science Technology, 27, 589-596. 

Fletcher C.A.J., 1988, Computational Techniques for Fluid Dynamics 1, Berlin, Germany, ISBN: 978-3-642-

58229-5. 

Goodwin J., Hughes R.W., 2000, Rheology for Chemists, Dunfermline, UK, ISBN: 9780854046164 

Macías-Salinas R., Aquino-Olivos M.A., García-Sánchez F., 2013, Viscosity Modelling of Reservoir Fluids over 

Wide Temperature and Pressure Ranges, Chemical Engineering Transaction, 32, 1573-1579. 

Magerramov M.A., Abdulagatov A.I., Abdulgatov I.M., Azizov N.D., 2006, Thermal Conductivity of Peach, 

Raspberry, Cherry and Plum Juices as a Function of Temperature and Concentration, Journal of Food 

Process Engineering, 29, 304-326. 

Pohn S., Kamarád L., Kirchmayr R., Harasek M., 2010, Design Calibration and Numerical Investigation of a 

Macro Viscosimeter. Czech Society of Chemical Engineering (ČSCHI), 61, 1099–1100. 

Scargiali F., Busciglio A., Grisafi F., Brucato A., 2013, Influence of Viscosity on Mass Transfer Performance of 

Unbaffled Stirred Vessels, Chemical Engineering Transaction, 32, 1483-1489. 

320

330

340

350

360

370

0 300 600 900 1,200 1,500

T
e

m
p
e
ra

tu
re

 [
K

]

time [s]

Cylindrical_wall_M Cylindrical_wall_CFD
R79/H150_M R79/H150_CFD
R96/H150_M R96/H150_CFD
R30/H130_M R30/H130_CFD

0

0.01

0.02

0.03

0.04

0.05

0.06

0.07

0 400 800 1,200 1,600

T
o

rq
u
e
 [

N
c
m

]

time [s]

Viscopakt measurement
CFD

570