Available online at [Link].

com

Applied Thermal Engineering 28 (2008) 450–459

[Link]/locate/apthermeng

Numerical simulation and experimental analysis of an industrial

glass melting furnace
A. Abbassi *, Kh. Khoshmanesh
Amirkabir University of Technology, Department of Mechanical Engineering, Tehran, Iran Hafez Ave., P.O. Box 15875-4413, Tehran, Iran

Received 12 December 2006; accepted 16 May 2007

Available online 2 June 2007

Abstract

This paper presents a three dimensional study of a gas-fired, regenerative, side-port glass melting furnace. The furnace is divided into
three regions: combustion space, glass tank and batch blanket. The combustion space is simulated separately, while the glass tank and the
batch blanket are coupled together as a whole. To couple the three regions, a new method is proposed that is more flexible and faster than
the existing conventional methods. The melting process of the batch blanket, the natural convection vortices of the glass tank and the
turbulence and the chemical reactions of the combustion space are all simulated. The simulation was preformed using Gambit-2 and
Fluent-6 software. The calculated temperatures are compared with actual values measured at different locations on the furnace surfaces.
Results indicate that the maximum relative error is less than 7.6%. This work can form the basis for the parametric analysis of the furnace
to reduce its fuel consumption.
 2007 Elsevier Ltd. All rights reserved.

Keywords: Glass furnace; Simulation; Combustion space; Glass tank; Coupling

1. Introduction the past three decades, remarkable progress has been made
in the use of computers to simulate the glass furnaces. This
Glass furnaces consist of a combustion space and a has enabled the numerical models to grow more rapidly in
glass tank. Fuel and air/oxygen are fed into the combus- comparison with other models. The simulations have
tion space and react with each other to release a large proved to be invaluable tools as replacements for expensive
amount of heat in the form of radiation and convection. and time-consuming testing. The introduction of advanced
Sand and cullet are fed into the glass tank and the processors and computational algorithms has enabled the
released heat makes them to participate in a series of three dimensional simulation of the glass furnace with
chemical reactions to be able to melt. The glass melt cir- more accurate treating of physical and chemical phenom-
culates in the glass tank to be homogenized until it leaves ena. The combustion reaction, the turbulent nature of com-
the furnace. bustion gases, the complicated melting process of batch
A lot of research work has been done to achieve the flow and the sophisticated behavior of gas bubbles and
reduction of fuel consumption and emissions and also solid particles remained in melt flow are of more interest
enhance the glass quality. Analytical, physical and numer- in the recent studies.
ical models have been employed as the main tools by McConnell and Goodson [1] presented a simplified
researchers to achieve the above-mentioned goals. Over model of whole furnace components. Three energy equa-
tions were solved for crown, batch and refractory temper-
ature. The radiative exchange in the combustion space
*
Corresponding author. Fax: +98 21 66419736. was calculated using the Hottel zonal method and the flow
E-mail address: abbassi@[Link] (A. Abbassi). patterns of glass tank were predefined. Mase and Oda [2]

1359-4311/$ - see front matter  2007 Elsevier Ltd. All rights reserved.
doi:10.1016/[Link].2007.05.011
 A. Abbassi, Kh. Khoshmanesh / Applied Thermal Engineering 28 (2008) 450–459 451

accomplished a two dimensional flow and temperature pat- 2. Mathematical simulation

tern for a glass furnace. The batch velocity over the melt
was assumed constant and an energy balance equation 2.1. Combustion space model
was solved to determine the batch temperature. The com-
bustion space temperature was simulated using the Hottel The thermal behavior of flame as well as combustion
zonal method. products is very complex due to turbulence, chemical reac-
Carvalho [3–5] introduced a comprehensive model com- tions and radiative exchange. The models applied to treat
posed of combustion space, glass tank and batch blanket the above phenomena are explained briefly here.
sub-models. The sand grain dissolution, redox and fining
processes were also taken into account. Chmelar et al. [6] 2.1.1. Turbulence model
offered a complete glass furnace model including three The turbulent nature of the combustion gases was simu-
dimensional glass melt flow with bubbling, combustion lated by applying realizable k–e model proposed by Shih
space and batch blanket simulations with the ability to con- [12]. This model was intended to address the deficiencies
sider electrical boosting. Lankhorst et al. [7] developed a of traditional k–e model by adopting two modifications.
complete glass furnace model including combustion space First, a new eddy viscosity formula involving a variable
and glass tank sub-models. The coupling was performed Cm was applied as proposed by Reynolds. Second, a new
through the distribution of heat flux on the melt surface. equation for the dissipation rate, e was applied based on
Hoke and Marchiando [8] presented a model using Fluent the dynamic equation for the mean-square vorticity
software, which included three dimensional combustion fluctuation.
space and glass tank models. The model was used to opti- Instead of standard wall function, a two-layer based,
mize the position of burners and firing rate with respect to non-equilibrium wall function was applied as introduced
glass quality. by Kimand and Choudhury [13]. The key elements in this
Golchert et al. [9] presented a spectral radiation heat wall function are as follows. First, Launder and Spalding’s
transfer model and a reduced NOx kinetics model coupled log-law is sensitized to pressure gradient effects. Second,
with a flow simulation code. This model was used to inves- the wall neighboring cells are assumed to consist of a vis-
tigate the radiative exchange, pollutant formation and flow cous sub-layer and a fully turbulent layer. Different profiles
behavior in the combustion space of a glass furnace. The of turbulent shear stress, production rate and dissipation
mathematical simulations of glass furnaces are reviewed rate are regarded in each layer. The model relaxes the local
comprehensively by Loch [10], in which the evolution and equilibrium assumption (production = dissipation) that is
trend of simulations are discussed. Recently Chang et al. adopted by the standard wall function at neighboring cells.
[11] accomplished a comprehensive computational model
using the Eulerian approach to simulate multiphase flows 2.1.2. Combustion model
in a glass tank. It included all the phases, divides solid par- The combustion reaction was simulated by applying the
ticles or gas bubbles into various size groups, and treats eddy break-up model proposed by Magnussen and Hjer-
each group as a continuum. tager [14]. The main feature of this model is that the influ-
In spite of progress achieved in three dimensional mod- ence of the turbulence on the reaction rate is taken into
eling of the combustion space, the glass tank and the batch account. The rate of reaction is proportional to the ratio
blanket, they either have not been coupled into an overall of dissipation rate to the turbulent kinetic energy of the
furnace model or have not been carefully validated against reactants and the products of combustion, e/k. The inverse
the operating furnace data. of this ratio represents the time scale of turbulent eddies
In this work, the furnace is divided into the combustion following the eddy dissipation model of Spalding [15].
space, the glass tank and the batch blanket. The glass tank
and the batch blanket are coupled to each other, where no 2.1.3. Radiation model
assumptions are made in their interfaces. However the The radiation exchange of hot gases was modeled using
combustion space is simulated separately. The main contri- the discrete ordinates radiation model (DORM) proposed
bution of this work is the introduction of a new method by Lockwood and Shah [16]. In this model imaginary
that is used to couple the combustion space to the glass spheres are considered around each cell. Each sphere is
tank-batch blanket models. In this method the heat flux divided into discrete solid angles, each associated with a
distribution over the glass tank’s and the batch blanket’s vector direction. The radiative transfer equation (RTE) is
free surfaces is estimated applying the mass flow rate, the solved along each vector to calculate the radiation inten-
temperature and the composition of the raw material, the sity. The model is useful to solve problems ranging from
refractory heat losses and the variation of gas flow rate surface to surface radiation to participating radiation in
through fuel ports along the furnace. Also the geometry combustion problems.
of the batch blanket is predicted by a simple technique. The radiation is emitted from the flame and combustion
The calculations are validated through the actual measure- products in different wave-lengths. This fact was taken into
ment of the temperature in eighteen locations on the fur- account by applying the weighted sum of gray gases as pro-
nace surface. posed by Hottel and Sarofim [17]. In this model, the total
452 A. Abbassi, Kh. Khoshmanesh / Applied Thermal Engineering 28 (2008) 450–459

emissivity of gases is considered as the sum of several imag- sure terms were coupled using SIMPLER algorithm in all
inary gray gases. The ratios used in summation correlation regions [21].
depend on the gas composition and the temperature and
were obtained from Rhine and Tucker [18]. 4. Physical properties

2.2. Glass tank model 4.1. Combustion space

The conservative equations include mass, momentum The thermophysical properties of combustion gases are
and energy equations. The molten flow is dominated by summarized in Table 1. The combustion gases can be
strong free convection cells created by severe temperature regarded as a mixture. The density of the mixture is defined
differences on the melt surface. The flow in the glass tank using the ideal gas law. The viscosity, the thermal conduc-
is laminar due to the low velocity and high viscosity of melt tivity and the specific heat of mixture are defined as a mass
[2,3,5]. The effect of radiation is considered by an effective fraction average of pure species. The specific heat of each
thermal conductivity [2,3,5]. species varies as a polynomial function of temperature.

2.3. Batch blanket model 4.2. Glass tank

The conservative equations include mass, momentum The thermophysical properties of the molten glass are
and energy equations. The batch blanket can be divided summarized in Table 2. As can be seen in the table, the den-
into three regions with respect to temperature. The region sity, the specific heat and the thermal expansion ratio of
with temperatures less than 1073 K is called unfused core, melt are constant, while the viscosity and the effective ther-
where the temperature of raw material is increased. The mal conductivity of melt strongly depend on temperature.
region with temperatures of 1073–1473 K is called chemical The molten flow is dominated by the natural convection
reaction zone, where reactions occur between the compo- force and is strongly susceptible for divergence from itera-
nents of raw material to form the glass melt. Finally the tive solution. To diminish this tendency the values of vis-
region with temperatures more than 1473 K is called mol- cosity and effective conductivity are kept constant for
ten layer, where a thin layer of molten glass flows until temperatures less than 1173 K [22].
sinking in the glass tank [19]. The equations are different
in each region. The energy equation of chemical reaction 4.3. Batch blanket
zone includes an enthalpy source term. But the energy
equation of molten layer includes a reduced effective ther- The thermophysical properties of the batch blanket are
mal conductivity to take into account the radiation summarized in Table 3. As can be seen in the table, the
exchange. The value of effective thermal conductivity is properties strongly depend on temperature as the batch
lower than that of glass tank due to the small thickness blanket can be divided into unfused core, chemical reaction
of molten layer [20]. zone and molten layer with respect to temperature. Gener-
ally the properties of unfused core are those of raw mate-
3. Method of solution rial, the properties of molten layer are those of molten

Gambit-2 was applied to generate the geometry and the Table 1

mesh generation of the furnace [21]. The combustion space Thermophysical properties of combustion gases [21]
was divided into six compartments according to the loca- Density (kg/m3) Incompressible-ideal gas
tion of ports. In each compartment, the combustion space Viscosity (kg/m s) Mass-weighted-mixing-law
was divided into 56 · 42 · 36 elements, while the air port Thermal conductivity (W/m K) Mass-weighted-mixing-law
Specific heat (J/kg K) Mixing-law
was divided into 30 · 20 · 24 elements. Around the fuel
ports, where chemical reaction occurs between fuel and
air, the concentration of the elements was increased by Table 2
applying unstructured elements. The glass tank and batch Thermophysical properties of molten glass [2,3,22,25]
blanket were divided into 40 · 40 · 206 and 40 · 10 · 100 Density (kg/m3) 2500
elements, respectively. The concentration of the elements
Viscosity (Pa s) 1813   T < 1073 K
was increased near the melt free surface and the batch 10 exp 2 þ T3600 T > 1073 K
573
boundary [22].
Thermal expansion 5 · 105
Fluent-6 was applied as the solver [21], in which the
coefficient (1/K)
finite volume method is used to descritize the governing dif-
ferential equations [23]. In the combustion space, the con- Specific heat (J/kg K) 1256
vection terms were modeled by QUICK method, while in Effective thermal 5.38 T < 1073 K
the glass tank and batch blanket they were modeled by conductivity 5.386  0.02168T + 0.00002058T2 T < 1073 K
(W/m K)
the second-order upwind method. The velocity and pres-
 A. Abbassi, Kh. Khoshmanesh / Applied Thermal Engineering 28 (2008) 450–459 453

Table 3 0.95
Thermophysical properties of batch blanket [3,19,22,25]

free surface emissivity

0.85
Density (kg/m3) 1400 Unfused core
Linear Chemical reaction zone 0.75
interpolation
2500 Molten layer 0.65

Viscosity (Pa s) 1813   Unfused core 0.55

10 exp 2 þ T3600
573 Chemical reaction zone and
Molten layer 0.45
0 5 10 15 20
Specific heat 1100 Unfused core
furnace length (m)
(J/kg K) Linear Chemical reaction zone
interpolation Fig. 1. The variation of emissivity over the melt and batch blanket free
1256 Molten layer surface along the furnace [3,22].
Effective thermal 1.5 Unfused core
conductivity Linear Chemical reaction zone assumed to be adiabatic. The emissivity of all ports and
(W/m K) interpolation surfaces was assumed as 0.6, except for the batch blanket
3.5–4.5 Molten layer
and the molten free surface whose emissivity varied along
Melting heat (J/kg) 523,000 Chemical reaction zone the furnace (see Fig. 1) as basically proposed by Carvalho
[3,22].
glass and the properties of reaction zone are calculated by
interpolation. The reactions are endothermic with a latent 5.2. Operating conditions of glass tank-batch blanket
heat of 523 kJ/kg. For temperatures less than 1073 K the
viscosity is kept constant [22]. The mass flow rate of raw material was 200 ton/day or
2.31 kg/s. The temperature of raw material was taken as
5. Furnace description and operating conditions 300 K and no cullet was mixed with raw material [22].
For bottom and side walls, the natural convection bound-
In this work, the furnace No. 2 of Abguineh glass melt- ary condition was applied. The overall heat transfer coeffi-
ing plant was simulated. The furnace is gas-fired, regenera- cient and surrounding temperature were taken as 10 W/
tive and side-port with a capacity of 200 ton/day container m2 K and 300 K, respectively. At glass tank exit, where
glass. The glass tank is 26.3 m long, 10.3 m wide and 1.5 m the melt flow leaves the glass tank and passes under the
deep. It is equipped with bubblers at its bottom but they shadow wall, the fully developed boundary condition was
were out of service during the measurements. The glass applied. For the beginning parts of the batch blanket that
tank extends 3 m outside the back wall where the raw mate- have not entered the combustion space yet, the natural con-
rial is fed. The combustion space has an arc shape crown vection boundary condition was applied. The heat transfer
with the minimum and maximum heights of 1.5 m and coefficient and surrounding temperature were taken as
1.8 m, respectively. There are six air ports on each side 10 W/m2 K and 300 K, respectively. The glass tank and
and each with three fuel ports. The direction of flame is the batch blanket control volumes were coupled to each
changed in every 20 min [22]. other. Therefore no additional assumptions were made
on their interface. The lower interface was considered as
a wall to prevent penetrating of solid particles to the glass
5.1. Operating conditions of combustion space tank [22].
According to our observations, the batch blanket pro-
At fuel ports, the gas velocity was given such that the ceeds to the end of port 2 under the operating conditions
flow rate of ports 1–6 was 120, 240, 460, 460, 440 and of furnace. The total length of the blanket is 11 m including
240 m3/h, respectively [22]. The temperature of fuel was the 3 m length outside the combustion space. The initial
taken 300 K. The turbulence intensity ratio and the height of blanket is 0.2 m. According to the densities of
hydraulic diameter of fuel eddies were taken 10% and raw material and the molten glass, 0.11 m of initial height
0.15 m, respectively. At air ports, the air velocity was given sinks in the glass tank and the rest floats on it. It is assumed
such that the excess air was 20% proportional to gas flow that the upper boundary of blanket remains flat until enter-
rate. The air was preheated to 1373 K by means of regen- ing the combustion space [22]. According to Ungan and
erator. The turbulence intensity ratio and the hydraulic Viskanta [19], the final height of the blanket is nearly one
diameter of air eddies were taken as 10% and 0.5 m, respec-
tively. At outlet ports the fully developed boundary condi-
tion was applied. At crown and side walls, the natural
convection boundary condition was applied. The overall
heat transfer coefficient and surrounding temperature were
taken as 10 W/m2 K and 300 K, respectively. Regenerator
extensions, air-fuel mixing space and shadow wall were Fig. 2. The geometry of the batch blanket [22].
454 A. Abbassi, Kh. Khoshmanesh / Applied Thermal Engineering 28 (2008) 450–459

fourth of the initial height, which is 0.05 m for our case. 30

Now, there is enough data to draw the geometry of the 25

weight function (%)

blanket (see Fig. 2). 20
15
10 initial
6. Coupling of combustion space to glass tank-batch blanket 5 modified
models 0
-5
-10
The conventional method to couple the combustion 0 2 4 6 8 10 12 14 16 18 20 22
space and the glass tank-batch blanket models was applied distance along furnace (m)
by Carvalho [3]. In this method the temperature of the
Fig. 3. The variations of initial and modified weight functions of heat flux
glass tank and the batch blanket free surfaces is estimated. distribution over the molten and batch blanket free surface, point 0 locates
Applying this boundary condition, the combustion space is on the back wall of the combustion space, points 1–6 located along the
simulated and the heat flux over the free surface is calcu- ports and point 7 locates on the front wall of the combustion space.
lated. Then using this heat flux, the glass tank and the
batch blanket are simulated and the temperature of free surface as Fig. 3. Some modifications should be considered
surface is calculated. The calculated and the guessed tem- on the initial weight function. A second order curve is
peratures are compared. If difference is significant the passed through the points 3, 4 and 6 to calculate the new va-
above operation is repeated until convergence is achieved. lue of point 5. The values of points 0 and 7 are calculated by
This coupling method has two problems. extrapolation. The point 0 is located at the back wall while
The middle parts of the free surface have the maximum the point 7 located at the front wall. The curve which passes
temperature and receive the maximum heat flux simulta- through the above points is the modified weight function of
neously. If the temperature boundary condition is applied heat flux distribution over the free surface.
on the combustion space the middle parts of the free sur- Third, the heat flux distribution is regarded as the prod-
face will receive small heat fluxes, because the combustion uct of the heat factor by modified weight function. The
space does not have any information about the glass tank. value of the heat factor is calculated dividing the total
In addition, coupling may need several iterations between absorbed heat by the integral of the modified weight func-
the combustion space and the glass tank-batch blanket tion and the furnace width. The modified weight function is
models. The simulations are three dimensional and each integrated along the free surface using numerical methods
iteration takes several days. To remedy the above problems like Trapezoid integration.
a new coupling method is proposed as follows [22].
q00absorbed ðzÞ ¼ weight functionðzÞ  heat factor ð6Þ
First, the total heat absorbed by free surface is esti- Z L
mated. This heat is consumed to melt the raw material, Qabsorbed ¼ q00absorbed ðzÞdz  width ð7Þ
above the melting temperature, to compensate the heat loss 0
through the refractory of the bottom and the side walls and Qabsorbed
also the heat loss by the free surface of batch blanket out- heat factor ¼ R L ð8Þ
0
weight functionðzÞdz  width
side the combustion space. The heat needed to melt the raw
material is calculated exactly through the furnace pull rate, Fourth, the heat flux is imposed on the combustion space
the temperature of raw material, the composition of raw and the glass tank-batch blanket models as the boundary
material and the percent of cullet as follows: condition. Simulating the glass tank and the batch blanket,
the portion of the heat lost through the refractory of the
Qmelting ¼ Qtemp rising þ Qreaction þ Qmelt heat ð1Þ

Q ¼ m_ Cp
temp rising ð1073  T raw material Þ
temp rising
ð2Þ
Qreaction ¼ m_ Cpjreaction ð1473  1073Þ þ m_ DH reaction ð3Þ
Qmelt heat ¼ m_ Cpjmelt ðT melt  1473Þ ð4Þ

The heat needed to overheat the molten glass above the

melting point is calculated assuming that the temperature
rise is 50 K. This heat is small and does not have a signif-
icant effect. The total heat loss through the refractory of
bottom and the side walls and the free surface of batch
blanket outside the combustion space are assumed as one
fourth of the sum of Qmelting.
Qabsorbed ¼ Qmelting  1:25 ð5Þ

Second, the flow rate of fuel ports 1–6 is considered as the Fig. 4. The points located on the side walls and the crown of the
initial weight function of heat flux distribution over the free combustion space, where temperature was measured [24].
 A. Abbassi, Kh. Khoshmanesh / Applied Thermal Engineering 28 (2008) 450–459 455

The most deviation is related to point 7 with an absolute

error of 126 K and a relative error of 7.6%. Next, there is
the point 6 with an absolute error of 112 K and a relative
error of 6.4%.
The accuracy of the model was compared against two
existing models. In the Webb model [26], the maximum dif-
ference between the measured and calculated temperatures
or simply the maximum error is 90 K for the points located
on the free surface of the glass melt and 92 K for the points
located on the crown surface. In the Chang model [11], the
maximum error is estimated to be 44 K for the points
located on the free surface of the glass melt. In our model,
on the other hand, the maximum error is 86 K for the
points located on the free surface of glass melt and 126 K
Fig. 5. The points located on the molten and the batch blanket free
for the points located on the crown surface. The accuracy
surface, where temperature was measured [24].
of this model is lower than that of Chang [11] but its main
contribution is to introduce a faster and more flexible
bottom and the side walls and the free surface of the batch method to couple the combustion space model to the glass
blanket outside the combustion space are calculated and tank-batch blanket model.
compared with the initial value of 0.25. If the difference To assess the speed of the proposed coupling method,
is significant the above steps are repeated until convergence two simulations were carried out. In the first simulation,
is achieved. the combustion space model was coupled to the glass
tank-batch blanket model by using the conventional
method. In the second simulation, the proposed coupling
7. Results and discussion method was applied. The first simulation converged after
126 h. Using the same computer, the second simulation
7.1. Experimental validation converged after 54 h. Reduction of computational time
by 72 h can be of great importance in glass plants if the
The calculated temperatures were compared with the model is used to manage the furnace in different situations.
ones measured at 18 points on furnace surfaces as shown It should be noted that Webb [26] has improved the
in Figs. 4 and 5. The points 1–6 are located on the side accuracy of his model by considering the soot formation
walls, while the points 11–16 on the free surface of the later. Also it should be pointed out that Chang [11] has pre-
batch blanket and the glass tank. The temperature of those sented the normalised form of computational errors in his
12 points was measured by an infrared pyrometer. During work. Therefore the reference temperature was assumed to
the measurements there should be no flame in the combus- be 1523 K to achieve the above-mentioned error.
tion space. Therefore, the measurements were performed in A portion of the error could be addressed due to the
the 10 second gap between the flame direction intervals. measurement errors including the deterioration of thermo-
The points 7–10 are located on the crown, while the points couples and the wrong emissivities of the side walls, the
17 and 18 are located on the glass tank bottom. The tem- crown and the melt free surface being set for pyrometer.
perature of those 6 points was measured by the furnace The other portion arises from the computational errors
thermocouples embedded in its walls. The measurements caused by assumptions in the turbulence, the combustion
were carried out by Khoshmanesh [24]. The calculated and the radiation models to simplify the problem.
and the measured temperatures are compared in Fig. 6. In order to enhance the accuracy of current model, the
complicated phenomena of furnace should be modeled
1900 more accurately. Effect of the soot formation in the com-
1800 15
14 13 3 4
2
bustion space and the bubble formation on the batch blan-
9 5
ket surface should be taken into account. Furthermore, the
calculated temp (K)

1700 16 1 8

1600 10 6
variation of the heat flux distribution across the furnace (x
direction) should be regarded. In this case, the heat flux dis-
1500 7
12 tribution on the free surface of the batch blanket and glass
1400 melt is defined as follows; where the variations of heat flux
1300 18 along and across the furnace are represented by two weight
11
1200 17 functions, respectively. This point could be of great impor-
1200 1300 1400 1500 1600 1700 1800 tance in the cross-fired furnaces.
measured temp (K)
q00absorbed ðx; zÞ ¼ weight functionðzÞ
Fig. 6. Comparison of calculated and measured temperatures of the
furnace [24].  weight functionðxÞ  heat factor ð9Þ
456 A. Abbassi, Kh. Khoshmanesh / Applied Thermal Engineering 28 (2008) 450–459

7.2. Combustion space

The temperature contours of combustion gases at the

middle planes of ports 1–6 are shown in Fig. 7. The maxi-
mum temperature is 2237 K and is related to the flames of
ports 3, 4 and 5. These ports receive the maximum gas flow
rate to ensure that the batch blanket sinks until the middle
of glass tank. Approaching the crown and the molten free
surface the temperature decreases to 1670 K. The path of
preheated air can be seen as a band of 1373 K, which pro-
ceeds to the middle parts of combustion space. Approach-
ing the outlet ports the temperature of the combustion
products decreases. The mean flue temperature is 1801 K,
which enters the regenerator to preheat the combustion
air. The path of fuel can be seen as a region of 300 K, Fig. 9. The contours of CO2 mass fraction at the middle planes of ports
1–6.
which is limited around the fuel ports.
The temperature contours of the molten and the batch
blanket free surface calculated from the combustion space the ports 1 and 2 the temperature remains nearly constant.
model is shown in Fig. 8. It can be seen that the maximum The direction of the flame changes in every 20 min and the
temperature is 1725 K and is related to the regions near the location of hot and cold regions changes. Approaching the
outlet port 5. Approaching the outlet ports of 3, 4, 5 and 6 furnace back and front walls the mean temperature
the temperature increases. For example along port 5 the decreases to 1590 K and 1695 K, respectively.
temperature increases from 1575 to 1725 K. But along The contours of the CO2 mass fraction at the middle
planes of ports 1–6 are shown in Fig. 9. Above the molten
free surface the concentration of CO2increases until it
reaches to 12.7% near the outlet ports. The path of the pre-
heated air can be seen as a band of zero CO2 concentration.
Above this band the CO2 concentration increases to 5.5–
6.9%. Multiplying the values of CO2 concentration by
two the concentration of H2O is calculated. In the fossil
fuel fired combustion processes CO2 and H2O are the most
important emitters of gaseous radiation [18]. More details
about the simulations can be found in Abbassi and
Khoshmanesh report [22].

7.3. Glass tank

The temperature contours of the glass tank at its longi-

tudinal middle plane are shown in Fig. 10. The maximum
temperature is 1825 K and is related to the free surface at
Fig. 7. The temperature contours of combustion gases at the middle
planes of ports 1–6. the middle regions of the glass tank. The mean temperature

Fig. 8. The temperature contours of molten and batch blanket free surface Fig. 10. The temperature contours of the glass tank along its middle
calculated from combustion space model. plane.
 A. Abbassi, Kh. Khoshmanesh / Applied Thermal Engineering 28 (2008) 450–459 457

of the batch blanket lower boundary is 1000 K, which cre- tex, then passes through the clockwise vortex and finally
ates a density difference of 103 kg/m3 with middle regions leaves the glass tank through the channel.
and feeds a free convection vortex. The mean temperature The temperature contours of the molten and the batch
of the free surface near the shadow wall is 1600 K, which blanket free surface, which were calculated from glass
creates a density difference of 28 kg/m3 with the middle tank-batch blanket model, is shown in Fig. 12. The maxi-
regions and feeds another free convection vortex. mum temperature is 1825 K and is related to the middle
The streamlines of the glass tank at its longitudinal mid- regions, which begin after the batch blanket boundary.
dle plane are shown in Fig. 11. The streamlines are colored Approaching the back and the front walls the temperature
by the molten velocity. The maximum velocity is 19.8 mm/s decreases to 1063 K and 1520 K, respectively. The begin-
and related to the free surface just after the batch blanket ning parts of the batch blanket have not entered the com-
boundary. The molten flow is dictated by two free convec- bustion space and their temperature is limited to 300 K.
tion vortices. A counter clockwise extends from the middle The temperature is constant across the glass tank. If the
regions of the glass tank to the back wall. This vortex plays imposed free surface heat flux is multiplied by a weight
an important role of heating the batch blanket lower function, which varies across the furnace the temperature
boundary, although some of its heat is lost through the variation across the furnace can be taken into account.
bottom and the side refractory of the glass tank.
Another clockwise vortex extends from the middle 7.4. Batch blanket
regions of the glass tank to the front wall. This vortex heats
the end regions of the glass tank. Some part of this heat is The maximum height of the batch blanket is 0.2 m and
emitted from the free surface to combustion space; the its contours can hardly be presented. Therefore, the varia-
other portion heats the molten flow and leaves the glass tion of the molten percentage at the middle plane of the
tank through its channel. The rest of the heat is lost glass tank is shown in Fig. 13. The upper boundary receives
through the refractory. The path of the molten particles direct radiation exchange form the combustion gases and
is determined by the interaction of those two vortices. therefore is more heated. The first melting reactions of
The particle first passes through the counter clockwise vor- the upper boundary begins at z = 1 m and are completed
at z = 5.425 m. The lower boundary receives convection
exchange from a free convection vortex and therefore is less
heated. The melting reactions of the lower boundary begin
at z = 3.8 m and are completed at z = 7 m. The distribution
of the unfused core, the chemical reaction zone and the
molten layer of the batch blanket can be seen from this
graph. The unfused core spreads from z = 3 m to
z = 3.8 m, the chemical reaction zone spreads from
z = 1 m to z = 7 m and finally the molten layer lay between
z = 3.8 m and z = 8 m.

7.5. Optimization of glass furnace

For a new furnace, the simulation can be used to opti-

mize its thermal performance. The furnace designer can
Fig. 11. The streamlines of the glass tank along its middle plane. change the furnace span, the combustion space height
and the location and flow rate of different ports. The tem-
perature distribution of the melt free surface and the

1
0.9
0.8
melting fraction

0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
-3 -2 -1 0 1 2 3 4 5 6 7 8
batch blanket length (m)

Fig. 12. The temperature contours of molten and batch blanket free Fig. 13. The variation of molten percentage at the middle plane of the
surface calculated from glass tank-batch blanket model. glass tank.
458 A. Abbassi, Kh. Khoshmanesh / Applied Thermal Engineering 28 (2008) 450–459

streamlines of the glass melt that are calculated in each case 2006. The authors wish to appreciate the industrial sector
guide the designer to select an optimum combination of of IFCO for their cooperation and financial support and
above values. Mrs. Sara Baratchi for her precious helps.
For an existing furnace, on the other hand, the simula-
tion can be used to modify its operating conditions in
response to different situations. For example, if one port References
is plugged, the simulation helps the operators to change
the flow rate or excess air of the other ports to manage [1] R.R. McConnell, R.E. Goodson, Modeling of glass furnace design
the new conditions. for improved energy efficiency, Glass Technology 20 (3) (1979) 100–
106.
Furthermore, the simulation can be employed to deter- [2] H. Mase, K. Oda, Mathematical model of glass tank furnace with
mine the portion of different heat consumptions in the fur- batch melting process, Journal of Non-Crystalline Solids 38–39 (1980)
nace. For example, this simulation revealed that the mean 807–812.
temperatures of the combustion flue, the combustion space [3] M.G. Carvalho, Computer simulation of a glass furnace, PhD Thesis,
London University, 1983.
refractory and the glass tank refractory are 1801, 1602 and
[4] M.G. Carvalho, M. Nogueira, Modelling fluid flow and heat transfer
1056 K, respectively. Those values were applied in the in an industrial glass furnace, in: Annual ESPRIT Conference,
energy balance equation of the furnace (see Eq. (10)). It Kluwer, Dordrecht, Netherlands, 1990, pp. 530–543.
was concluded that 4.4 MW (20%) of the total energy is [5] M.G. Carvalho, M. Nogueira, J. Wang, Mathematical modeling of
consumed for glass melting and another 7.4 MW (32%) is glass melting industrial process, in: 17th International Congress on
recovered in the regenerator. The rest of energy is lost Glass, 6, 69–74, Beijing, China, 1995.
[6] J. Chmelar, P. Schill, A. Franek, Mathematical models of glass
through the combustion space refractory, the glass tank melting furnaces, in: 4th International Conference on Advances in
refractory and the stack flue as 4.6 MW (21%), 1.6 MW Fusion and Processing of Glass, Wurzburg, Germany, 1995, pp. 63–
(7%) and 4.4 MW (20%), respectively. 66.
[7] A.M. Lankhorst, G.P. Boerstoel, H.P. Muysenburg, Complete
Qfuel þ Qpreheating air  Qflue ¼ Qmelting þ Qglass tank loss simulation of the glass tank and combustion space of the former
ford Nashville float furnace, including melter, refiner and working
þ Qcombustion space loss ð10Þ end, in: 4th International Seminar on Mathematical Simulation in
Glass Melting, Czech Republic, 1997.
The above conclusions motivated the management of [8] B.C. Hoke, R.D. Marchiando, Using combustion modeling to
the Abginue Company to enhance the thermal efficiency improve glass melt model results, 4th International seminar on
of the furnace. The proposed solutions included further mathematical simulation in the glass melting, Horni Becva, Czech
preheating of the combustion air and raw material, better Republic, 1997, pp. 87–95.
insulation of the combustion space refractory as well as [9] B. Golchert, C.Q. Zhou, S.L. Chang, M. Petrick, Investigation of
spectral radiation heat transfer and NOx emission in a glass furnace,
application of porous burners. in: International Mechanical Congress and Exposition, Orlando,
Florida, 2000.
8. Conclusions [10] H. Loch, D. Krause, Mathematical simulation in glass technology,
Springer, 2002.
[11] S.L. Chang, C.Q. Zhou, B. Golchert, Eulerian approach for multi-
Results presented in this paper indicate the capabilities phase flow simulation in a glass melter, Applied Thermal Engineering
of the two models developed for three dimensional simula- 25 (2005) 3083–3103.
tion of a complete glass furnace. The glass tank and the [12] T. Shih, W. Liou, A. Shabbir, J. Zhu, A new k–e eddy-viscosity model
batch blanket are simulated in tandem to each other, while for high Reynolds number turbulent flows model development and
the combustion space is simulated separately. The main validation, Computers Fluids 24 (3) (1995) 227–238.
[13] S.E. Kimand, D. Choudhury, A near-wall treatment using wall
contribution of this work is the introduction of a new functions sensitized to pressure gradient, ASME FED, vol. 217, 1995.
method to couple the combustion space to the other [14] F. Magnussen, B.H. Hjertager, On mathematical models of turbulent
regions. In this method, the heat flux distribution over combustion with special emphasis on soot formation and combustion,
the glass tank and the batch blanket free surfaces is esti- in: 16th Symposium on Combustion, 1976, pp. 719–729.
mated through the mass flow rate, the temperature and [15] B. Spalding, Mixing and chemical reaction in steady confined
turbulent flames, in: 13th Symposium on Combustion, 1971, pp.
the composition of raw material, the refractory heat losses 649–657.
and the variation of gas flow rate through the fuel ports [16] F.C. Lockwood, N.G. Shah, A new radiation solution method for
along the furnace. The actual measurements of the temper- incorporation in general combustion prediction procedures, in: 18th
ature indicate that the maximum relative error of calcula- Symposium on Combustion, 2003, pp. 1405–1414.
tions is less than 7.6%. The results of the simulations can [17] H.C. Hottel, A.F. Sarofim, Radiation Heat Transfer, McGraw-Hill,
1967.
be used to optimize the thermal efficiency of the furnace. [18] J.M. Rhine, R.J. Tucker, Modeling of gas-fired furnaces and boilers,
British Gas, 1991.
Acknowledgements [19] A. Ungan, R. Viskanta, Melting behavior of continuously charged
loose batch blankets in glass melting furnace, Glastech. Ber. 59 (1986)
279–291.
This work was carried out under the contract 1205–1577 [20] R. Viskanta, X. Wu, Effect of radiation on the melting of glass batch,
between Amirkabir University of Technology and Iran Glastech. Ber. 56 (1983) 138–147.
Fuel Conservation Organization (IFCO) during 2005– [21] Fluent 6.2, User’s guide, Fluent Inc., 2002.
 A. Abbassi, Kh. Khoshmanesh / Applied Thermal Engineering 28 (2008) 450–459 459

[22] A. Abbassi, Kh. Khoshmanesh, H. Abbassi, Simulation of combus- [25] F. Simonis, H. De waal, R.C.G. Beerkens, Influence of furnace design
tion and melting in a glass melting furnace, some applied suggestions and operational parameters on the residence time distribution of glass
to reduce the fuel consumption, Amirkabir University of Technology, tanks predicted by 3-D computer simulation, in: 14th International
Contract No. 1205–1577, 2006. Congress on Glass, New Dehli, 1986, pp. 118–128.
[23] H.K. Versteeg, W. Malalaskera, An Introduction to Computational [26] B.W. Webb, M.Q. McQuay, Development, experimental validation,
Fluid Dynamics: The Finite Volume Method, Longman, 1995. and application of advanced combustion space models for glass
[24] Kh. Khoshmanesh, Mathematical simulation and experimental mea- melting furnaces, The Final Report, Prepared for The United States
surements of a glass melting furnace, MS Thesis, Amirkabir Univer- Department of Energy, Agreement DE-SCO2-95CE41122, 2001.
sity of Technology, Tehran, Iran, 2000.