ENG  Vol.12 No.9 , September 2020
Analytical Modeling and Computer Simulation of Heat Transfer Phenomena during Hydrothermal Processing Using SOLIDWORKS®
Abstract: Mathematical modeling of unit operations especially in biomaterials processing plays an important part in understanding simultaneous heat and mass transfer and fluid flow phenomena. Hydrothermal/supercritical processes are one such process which utilizes high temperature and pressure to synthesize materials for varied applications. The present study constitutes development of a model incorporating all standard transport equations and transient conditions to predict the behavior of process and thus materials upon heating to high temperature in an enclosed composite (Steel and Teflon (PFA) vessel. Commercial software package SOLIDWORKS® is employed to simulate and output is presented as animations after post-processing. Results yield very useful data and information about process; materials of construction and materials being processed which helps in optimization.

1. Introduction

Hydrothermal treatment [1] is a thermochemical process [2], widely used in various manufacturing processes [3] [4] [5]. These treatments are performed above 100˚C in an aqueous solution [6]. It is an old technique that has been employed in treatment of sewage water, power industry, synthesis of chemicals [1] [7] [8], nanotechnology, and petroleum industry. Supercritical hydrothermal processes involve high temperature (Tc) and pressure (Pc), usually above the critical temperature of the water (Tc = 647 K & Pc = 221 bar). Precise calculation of the heat generation and transfer can give more control over the hydrothermal process. Any hydrothermal process involves two steps i.e. heat generation and heat transfer. Heat generation may be attained via electrical power [9] [10]. Heat transfer depends upon the process being carried out, type of insulation and configuration of the hydrothermal reactor. For a successful modelling, it is important to incorporate all temperature & pressure dependent material properties. External heating elements are placed for more precise control of temperature during the hydrothermal treatment. The heat generation (Q) in joules is in accordance with the Joule heating law (Q = I2Rt), where I is the current in amperes, R is the resistance of the heating elements and t is the duration in seconds. Heating elements may be constructed through high electrical resistance materials (stainless steels, tungsten alloys & nickel chromium alloys).

Joule heating is also known as Ohmic heating and it can be thought of the similar heating process used by microwaves but at different frequencies. The main advantage of the Joule heating is the uniformity and better control in a hydrothermal process. It depends upon the electrical resistance of the heating elements, electrical field strength, residence time and type of waveforms. For the safe and reliable operation, it is important to understand factors on which heat generation depends and modelling is efficient way to understand the process. Some advantages of the Joule heating are homogeneity, better control, less time, more efficient and environmentally friendly [9] [10]. Second step involves heat transfer process that depends upon the process being carried out inside the reactor and material properties [11] [12]. To avoid the loss of the heat and for better control over the process, usually such hydrothermal reactors are lined with some insulating layer such as Teflon. The heat transfer process is analyzed via standard transport equations [13] [14]. The main challenges involved during modelling of the hydrothermal operations are lack of material thermal properties at various temperatures and pressures, adiabatic increase or decrease in the temperature due to local pressure variation. Incorporation of heat transfer coefficient can result in inaccuracies in the model and measures have been adopted to counter this. The main purpose of the study is to develop the numerical simulation for the heat generation and transfer process during hydrothermal operation. It aims to develop a computational platform for the heat control (mainly for generation and transfer). The process is analyzed via mathematical modelling and verified through simulation using commercial code (SOLIDWORKS®). The heat generation is calculated via Joule heating law and heat transfer process is modelled via standard transport equations [13] [14]. Previous studies have incorporated only single value of the heat transfer coefficient for the heat transfer throughout the process, however in real operation, it is a transient in nature [14]. The major contribution of this work is the incorporation of the transient heat transfer coefficients, which are calculated in each iteration. It provides more accurate model and precise values of output temperature. A model is developed and implemented via iterative approach via constitutive equations for the calculations of heat transfer coefficients during each step, based upon the correlations developed for conduction, convection, and radiation [14]. Heat generation and heat transfer simulations are carried out in SOLIDWORKS® software package for verification purposes. A computer simulation-based approach is developed for the accurate prediction of heat transfer during the hydrothermal process.

2. Mathematical Formulation

To formulate a generalized mathematical model corresponding to heat transfer phenomena during hydrothermal processing in sealed vessels, the following assumptions were made:

1) Main source of heating is electrical resistance heating in a Box type furnace. Electric heater is at constant temperature.

2) Thermal properties of mold are not constant & do change with change of temperature.

3) Thermal conductivity of mold (steel) is very high and that of TEFLON is very low.

4) Rate of heat transfer by convection and radiation is negligible.

5) The furnace is very well insulated.

6) And most importantly, this model will consider the incorporation of the transient heat transfer coefficient during the modelling process.

Based upon these assumptions, following is considered

The steel and TEFLON (composite) vessel is considered identical to sealed pressurized reactor/vessel in which quantity of heat to and from system is zero (Adiabatic process).

3. Heat transfer for Hydrothermal Reactor

In order to formulate a mathematical model for heat transfer problem during heating, consider a stainless-steel engulfed TEFLON vessel placed inside an electric furnace (resistance heating).

There are three steps in which heat is transferred:

1) First step in which electrical power is used to heat the resistance element of furnace to desired temperature.

2) Second step in which heat energy generated in step one is used to heat the vessel (steel).

3) Third step in which heated steel vessel heats the TEFLON vessel.

Time of heating is determined separately during each interval using different relations and then individual times are summed up together to get total time taken for heating.

3.1. First Step of Heating—Heat Generation

When furnace is switched ON, electricity starts flowing from power (electrical) generation source into the elements. Material(s) of elements offer resistance to the flow of electricity. Consider a Voltage V flowing through resistor of resistance R in time t, then energy (heat) generated may be written as

Q t = V × I × time (1)

Q t = V 2 / R × time (2)


Qt = quantity of power consumed (Watt-hour) & heat energy produced (KJ).

V = Voltage (volts) (standard, different for different localities).

R = Resistance of heating element (ohms) (depend on cross-sectional area, length & electrical resistivity of material—changes with change of temperature) [9].

t = time of heating (hr).

All the heat generated in this step goes into overcoming material heat losses expressed in term of specific heat capacity of substance when it is heated across a temperature difference [14]. This may be written as.

Q t = m C p Δ T (3-1)

Q t = m C p ( T 2 T 1 ) (3-2)


m = mass of heating element (kg).

Cp = specific heat capacity of material (function of temperature) (KJ/Kg∙K).

T2 = final temperature (K).

T1 = initial temperate (K).

Comparing (2) & (3)

m C p ( T 2 T 1 ) = V 2 / R × time (4)

time = ( m C p ( T 2 T 1 ) × R ) / V 2 (5)

This is the time in which elements reach the set operating temperature (T2) from initial temperature (T1). Putting back the calculated time in (2) yields energy (heat) generated (hence available for transfer to raise the temperature of vessel).

Various factors affect time of heating during first step such as thermal conductivity of mold material [12] [15] [16], heating conditions, maximum power rating of furnace [9], safe operating limit, specific heat of metal, etc. All these should be taken into account while designing a heating system [10]. Soon after the element is heated to a designated temperature, second step of heating (heat transfer) begins.

3.2. Second Step of Heating—Heat Transfer

When resistance heating element is heated to its set operating temperature, the heat it produces (using electrical power) goes into raising the temperature of hydrothermal vessel [10] [12] [17]. Time to heat the vessel (bringing inside wall temperature to reaction temperature) is determined in that step using standard energy (heat) generation by electricity & heat transfer (conduction) equation(s) and is represented as time of reaction [14].

Energy (heat) generated once again may be written as

Q t = V 2 / R × time (6)

This is the same as (2) above, with calculated time value from (5). This heat is responsible for raising the outside and then inside temperature of vessel. Its quantity keeps on changing (raising) with time as more and more time (which will become visible later) is needed to bring the inside temperature to reaction temperature. However, its average value remains constant (actually used in calculations), by automatic control of furnace. When vessel is placed in furnace, the outer (infinitesimal) skull of stainless steel body is instantaneously heated to high temperature (Reaction temperature, temperature of furnace chamber). Heat transfer also starts subsequently and temperature of reactor body as a whole starts increasing. As a result of this, temperature of furnace chamber starts decreasing. This is compensated by automatic switching ON and OFF of furnace elements. This in turn depends on the programming of furnace (not discussed here) such that every time chamber temperature raises a certain degrees above set operating temperature, Furnace shuts OFF, remain shut OFF, until the chamber temperature drops below another designated temperature, when it automatically powers ON and starts heating. In this way, the average temperature (thus heat produced) of furnace remains constant. The reliability of furnace to function satisfactorily between these limits depends on factors such as materials of construction of elements (type (KANTHAL, Nichrome, SiC and Stainless Steel), quality, geometry and amount), insulation (type, quality, geometry & amount), type & materials of construction of furnace, temperature reading, sensing & reporting devices (thermocouples, sensors, indicators, displays) i-e furnace will transfer the heat positively for longer period of time even in shut OFF condition (thus saves electricity) if aforementioned variables are of optimized value.

Heat transfer from high temperature T2 to lower temperature T1 may be represented by following equation

Q T = ( T 2 T 1 ) / R t (7)

where T2 and T1 are outside and inside temperatures (K) of stainless steel vessel at infinitesimal value of time t and Rt is the thermal resistance of materials (K/W) at the same time.

4. Determination of Thermal Resistance

This thermal resistance can further be expressed in the form [14]

R t = L / k A (8)


L = thickness of mold material(s) (m).

A = Surface area exposed/available for heat transfer (m2).

k = Thermal conductivity of material (W/m-K) (strong function of temperature).

For stainless steel (304) k may be expressed by general relation

k = 9.705 + 0.0176 T 1.60 × 10 6 T 2 (9)


T = temperature of operation (K) at infinitesimal time t.

This thermal conductivity is different for different materials and its temperature dependence is also different for different materials. A huge collection of literature is available to find it for different materials at different temperatures.

5. Results and Discussion

Results and discussions are presented in the form of readings taken both by simulation and by thermocouple. Effect of various output values is presented in the form of temperature as input function. Results are presented in animation form yielded by simulation software package to explain and co-relate the phenomena occurring during heat transfer. These observations primarily show time dependent behavior of heat transfer phenomena and pattern at various locations of vessel (inside, outside and across the wall). The readings were taken at small time intervals to take advantage of better and efficient heat transfer pattern across wall of vessel. The simulations results are performed at fine mesh and take into account various modes of heat transfer from open, exposed and enclosed vessel surfaces. All simulations, graphical and experimental results complement each other and exhibit similar patterns of heat transfer with the passage of time till a condition is reached where all inner surface of vessel gets heated to desired operating temperature (250˚C) uniformly. The simulation performed also takes into account small variation in geometry (contours, curvy surfaces, fillet radius, champers). The model also shows the effect of time on different properties of vessel material and describes heat transfer variation across mold wall with change of time and temperature. It also takes into account the dynamic change in heat transfer coefficients with time along with change in quantities of heat(s) released as affected by change of input and vessel material temperature.

5.1. Effect of Metal Temperature

When simulation of temperature of vessel and its variation across vessel wall thickness performed as a function of time, following pattern is obtained which is shown in Figure 1. As can be seen when outside wall of vessel is heated uniformly by switching ON the heating elements of furnace, the outmost skin of vessel quickly reaches temperature of heaters which is agreement with assumption 1. The vessel is heated longitudinally from outside which gives rise to this effect. It can be easily seen that heat transfer pattern at time t = 0 results in maximum heat being transferred by conduction. It also shows that conduction is appreciably affected by variation in geometry i-e the same pattern is not uniform at the base. This is because at the base, there is no source of heat. Heat is only being supplied from elements installed parallel to walls of cylinder which does not have appreciable effect towards base. This heat transfer from the elements to base via conduction is minimum there and also because of radiation from inner surface of base, heat transfer from base to ambient is maximum there that results in the coldest spot in whole heat transfer pattern. This cold spot has its own benefits and drawbacks. Benefits being, it results in minimum deterioration on the surface, resulting in increase in service life of cylinder. It also helps avoid unnecessary heating of material since maximum mass is available at base (i.e. both effects (a) cold spot at base and (b) incipient heating due to heated mass are

Figure 1. (a) Mesh size of cylinder used in current simulation study (Total elements = 26,346, Total Nodes = 49,250) (b) Steady state heat transfer pattern at outside wall temperature T1 = 250˚C, time (t) = 0 and Node = 1 (c) same heat transfer pattern observed from top (isometric view) (d) Top (non-isometric view) (e) side view.

compensated at the base resulting in overall improvement in efficiency). Drawbacks are, this cold spot, if not properly controlled, can result in no heating of material at the base especially in the earlier part of heat transfer when overall mass is cold as a whole which eventually will incur more effort on heating elements to supply more heat to compensate for that which in turn finally can be in direct conflict with over heat transfer calculations carried out. Care must be taken to tune the heating from heating elements in such a way that they supply maximum heat at the start of heating cycle and avoid cold spot and once base and maximum mass of material at base, reaches high enough temperature, the heating from elements may reduce in conjunction with heating due to incipient heated mass itself to account for, and cater the, remaining heat, and heat transfer requirements. This is the best way to optimize the process. An industrial example of this phenomena is microcontroller based back loop heating programs such as one supplied by Inductotherm UK®. In metallurgical hearth type furnaces this effect is taken care of by introduction of heated metal parts from within furnace (which have been heated to high temperature already) at the top mouth from time to time. This ensures that top does not suffer from cold spot. In Foundry practice, the same phenomena could be avoided by introduction of purpose made heating pads (a self-burning refractory material) which gets fire when coming in contact with liquid metal and remains at high temperature for a long enough time thus reduces the chances of getting cold spots especially in isolated regions of mold.

Any other form of heating (such as transverse (Figure 1(e))) will only result in concentration of heart in a region near the horizontal half plane of vessel and dissipation of heat from there outwards gradually (i.e. decreasing) towards edges of vessel. This form of heating pattern is not desired as this is inefficient heating and gives rise not only to non-uniform pattern of heating but waste of maximum amount heat to surroundings. Further, this will unnecessarily heat the vessel at its horizontal half plane deteriorating the properties of stainless steel itself from middle. Furthermore, the material inside the vessel will also not get heated to its desired value as most of its mass is below horizontal half plane and in this configuration, maximum heat will never reach there. The variation in the temperature with respect to various conditions is shown in Figures 1(a)-(e).

Clearly, it can be seen that heat transfer occurs not only along straight line drawn perpendicular to tangent to outside surface but it occurs at various angles along tangent which results in a 3D steady state heat transfer pattern which gives very useful information about properties of material (stainless steel) as a function of temperature. It clearly shows that inward heat transfer is maximum along horizontal center line (or central cross-section region) even when heating elements are applied longitudinally and imparting heat uniformly all along length of cylinder as shown in Figure 2. This is because, the effect of heat transfer occurs not only by conduction but also by convection and radiation to open atmosphere at places near the surface and bottom of cylinder due to their close proximity to ambient. Although, their effects are negligible and can be ignored but in actual practice and in simulations performed at very fine mesh size, their effect becomes prominent and contributes towards overall heat transfer thus resulting in variant heat transfer pattern from point to point. This effect also results in deviation from steady state nature of heat transfer to non-steady state heat transfer. This effect can be avoided by the use of wool or insulating pads at top end and below bottom of cylinder which results in blockage of heat transfer and contributes towards decrease in temperature drop at such appreciable value.

5.2. Effect of Time

As the time changes from time t = 0 to 1, 2, 3, …, all the thermodynamics properties and parameters in the casting change as well. Especially the phenomenon predicted is a strong function of changes in heat transfer coefficients (HTC). All three HTC from convection, conduction and radiation change consequently. Though minutes, the changes in HTC strongly affect the overall temperature and temperature distribution (profile) across various sections of casting including temperature variation across wall of the vessel as shown in Figures 3(a)-(h).

Figure 2. (a)-(c) shows actual heat transfer profiles across different sections of cylindrical vessel (a) 200 mm depth from top (b) at 400 mm from top and (c) at 500 from top (clearly it can be seen that maximum heat transfer efficiency is at distance of 400 mm (almost middle of cylinder).

Figure 3. (a)-(h) shows heat transfer profiles across different sections of cylindrical vessel at time t = 0.

These figures show temperature profile of an element chosen randomly at a certain point from cross section of wall of cylinder which is subjected to change in HTCs. This change in HTCs; result in variation in temperature from outside maximum (250˚C) to inside maximum (250˚C). That is with the change of time, as HTC change, the temperature across thickness of wall becomes equal on both sides indicating maximum heat transfer and uniform attainment of desired operating temperature on the inside (which is the objective of present study). It is also worth to note that initial change in HTC is considerably large (i.e. change of whole numbers) but as the time passes, to attain more efficient temperature profile on the inside of wall, and to accurately predict the temperature, the change in HTC become extremely narrow and vary only in fourth decimal place. This change is very small (in present case 0.00000004578 W/m2/K) and must be noted by carefully noting the values in simulation package as they reach a point where temperature at both faces of wall becomes uniform (i.e. 250˚C). This is the point which is the objective of present study and value of time and HTCs gives exact measure of efficiency of heat transfer process. The value of time at which inside temperature becomes equal to 250˚C is the time needed for completion of process. The values of HTC before and after this changes very quickly and before this, only a range of change of temperature can be observed in Figures 4(a)-(p).

While after that only one temperature is observed (not shown here) (both are not the objective of heat transfer process and are representation of in efficient heat transfer process. Only maximum efficiency is possible at value of HTC where both surfaces reach equal temperature). Accurate measurement of HTC

Figure 4. Temperature profile of an element from cross-section of wall of cylinder with varying values of HTC.

also provide path to vary the process in such a way as to create measures in which time required to reach equilibrium could be reduced and kept constant for a longer period of time. One way of doing this is to insulate the outside wall of vessel with glass wool or efficient refractory to ensure that no or minimum heat loss is suffered from outside wall to ambient and maximum heat is utilized in raising the temperature of wall. Another method is to vary and tune the HTC in such a way to account for and accommodate iteration values back into original values with the shortest possible time interval (usually user defined in a simulation program (e.g. C++ [13] [14], MATLAB®)). This iteration (usually expressed as iteration coefficient) in turn can be increased to several hundred units which will reduce the time in which HTC across both faces of wall will become uniform (also called smoothing of curve/values). Yet another method is to control the radiation from inside of vessel to outside (black body effect). This will also result in minimum heat loss from within and allow it to concentrate to inside vicinity which in turn help in raising the inside temperature of wall and quickly reach a point where both faces reaches equilibrium in temperature. Last but not least, way to improve time efficiency involved in heat transfer in hydrothermal process is the proper and efficient selection of surfaces from which heat transfer could occur in practical scenarios. This involves champers, fillet radius corners, contours etc. This is explained in detail in preceding section. Although, small in values for small designs (such a present case) this could considerably improve or hamper the overall heat transfer process as one improperly selected surface could result in extremely bad heat transfer patterns which result in heat accumulation or heat sink across different faces of vessel thus resulting in localized hot or cold spot which decreases the life of vessel considerably and decreases the efficiency of process as well.

6. Defects

Six types of defects are identified during simulation which may give rise to poor results as shown in Figure 5.

Figure 5. Six types of defects are identified during simulation which may give rise to poor results. (a) Radial heating only; (b1, b2, b3) Improper surface selection; (c1, c2) Convection from outside of base; (d) Zero radiation from inside base; (e1, e2) Radiation plus Convection from base (too much heat transfer – dead zone); (f1, f2, f3). Base hating only.

1) Radial heating (along horizontal half plane).

2) Improper surface selection.

3) Convection from outside of base (improper base insulation/protection) (cold spot outside as well as inside).

4) No radiation from inside base (it results in accumulation of heat at the base—deteriorates properties).

5) Radiation + convection from outside base (too much heat transfer resulting into dead zones).

6) Base heating only.

Results obtained by simulations are presented as animations to check the validity of model and accuracy and consistency of results. These observations primarily show time dependent behavior of heat transfer phenomena and determine the conditions of casting at different intervals of time. The readings were taken at small intervals of time to get better representation of thermal heat transfer pattern across mold surface. Both simulation and experimental results (visual observations) complement each other and exhibit similar patterns of heat transfer with the passage of time till a condition is reached at which metal losses all its heat and completely solidifies. The model also shows the effect of time on different properties of mold material as well as gives heat transfer pattern across mold wall with change of time and temperature. It also takes into account the dynamic change in heat transfer coefficients with time along with change in quantities of heat(s) released as affected by change of metal and mold temperature.

7. Conclusions

Inverse relationship is found between time and inside/outside temperature as shown in the following Table 1.

Table 1. Relationship between inside/outside temperature and time.

TEFLON/Steel composite has shown good thermal resistance and cooling channels in the walls can give us more precise control over the heating and cooling of the reactor. Temperature is found to be a function of geometry, heating location, distance between heating source and vessel, and heat transfer properties of air/medium and vessel material. Rate of heat transfer, thermal conductance (K) decreases with time whereas Rt increases and Cp has shown variation with the time. Overall a mathematical model for the hydrothermal reactor has been developed using SOLIDWORKS® and employing transient heat transfer coefficients.

Cite this paper: Rafique, M. , Shah, U. (2020) Analytical Modeling and Computer Simulation of Heat Transfer Phenomena during Hydrothermal Processing Using SOLIDWORKS®. Engineering, 12, 682-697. doi: 10.4236/eng.2020.129048.

[1]   Byrappa, K. and Adschiri, T. (2007) Hydrothermal Technology for Nanotechnology. Progress in Crystal Growth and Characterization of Materials, 53, 117-166.

[2]   Sadat-Shojai, M. (2009) Preparation of Hydroxyapatite Nanoparticles: Comparison between Hydrothermal and Solvo-Treatment Processes and Colloidal Stability of Produced Nanoparticles in a Dilute Experimental Dental Adhesive. Journal of the Iranian Chemical Society, 6, 386-392.

[3]   Yoshimura, M. and Byrappa, K. (2008) Hydrothermal Processing of Materials: Past, Present and Future. Journal of Materials Science, 43, 2085-2103.

[4]   Zhang, X. and Vecchio, K.S. (2007) Hydrothermal Synthesis of Hydroxyapatite Rods. Journal of Crystal Growth, 308, 133-140.

[5]   Nian, J.-N. and Teng, H. (2006) Hydrothermal Synthesis of Single-Crystalline Anatase TiO2 Nanorods with Nanotubes as the Precursor. The Journal of Physical Chemistry B, 110, 4193-4198.

[6]   Rafique, M.M.A. (2018) Hydrothermal Processing of Phase Pure and Doped Hydroxyapatite and Its Characterization. Journal of Encapsulation and Adsorption Sciences, 8, 19.

[7]   Ashok, M., et al. (2007) Growth and Characterization of Hydroxyapatite Crystals by Hydrothermal Method. Journal of Materials Science: Materials in Medicine, 18, 895-898.

[8]   Earl, J.S., Wood, D.J. and Milne, S.J. (2006) Hydrothermal Synthesis of Hydroxyapatite. Journal of Physics: Conference Series, 26, 268-271.

[9]   Edward, H. (2010) Hughes Electrical and Electronic Technology. Pearson Education, London.

[10]   Theraja, B.L. (1977) Textbook of Electrical Technology. S. Chand, New Delhi.

[11]   Bergman, T.L., et al. (2011) Fundamentals of Heat and Mass Transfer. Wiley, Hoboken.

[12]   Holman, J.P. (1968) Heat Transfer.

[13]   Rafique, M.M.A. (2015) Modeling and Simulation of Heat Transfer Phenomena. In: Kazi, S.N., Ed., Heat Transfer Studies and Applications, IntechOpen, London, Chapter 9, 225-257.

[14]   Rafique, M.M.A. and Iqbal, J. (2009) Modeling and Simulation of Heat Transfer Phenomena during Investment Casting. International Journal of Heat and Mass Transfer, 52, 2132-2139.

[15]   Çengel, Y.A. and Boles, M.A. (2006) Thermodynamics: An Engineering Approach. McGraw-Hill Higher Education, New York.

[16]   Heine, R.W., Loper, C.R. and Rosenthal, P.C. (1967) Principles of Metal Casting. McGraw-Hill, New York.

[17]   Taylor, H.F., Flemings, M.C. and Wulff, J. (1959) Foundry Engineering. Wiley, Hoboken.