Multiphysics Simulation of Infrared Signature of an Ice Cube

This poster presents the Multiphysics numerical methodologies used to simulate the Infrared (IR) signature of an ice cube when taken out of the cold environment (-28 ooC) and allowed to warm at room temperature conditions by means of natural convection. The ice cube is chosen to simplify the geometry. The cubical geometry can be discretized simply and hence allows easy application of Multiphysics simulation tools. The aim of this work is to validate the various numerical methodologies. In order to model the behaviour, a 3D transient heat equation is solved using three different methodologies. In the first part of the work, the finite difference method is used to discretize the heat equation and solved using an FTCS (Forward-Time Central-Space) method in MATLAB® software. Then the same problem is modelled using the spectral method where the domain is discretized non-linearly for the appropriate solution. In the final part of the work, the problem is modelled in ANSYS® Multiphysics software. The results obtained through all methodologies are found in close agreement. This proves that Multiphysics tools can be employed to model thermal behaviour and hence predict IR profile. Further work will be carried out to compare with the experimental results. Thermal Diffusion 3D Domain 2D Cross-sectional View (Temperature Profile) 3D Visualisation of Thermal Diffusion Iso-surface Mapping H. Khawaja, T. Rashid, O. Eiksund, E. Brodal, K. Edvardsen UIT The Arctic University of Norway, Tromsø, Norway Conclusion The results from three different numerical methodologies FTCS (Forward-Time Central-Space) method, Spectral method and ANSYS® Multiphysics are found in good agreement. This proves that Multiphysics tools can be employed to model thermal behaviour and hence predict IR profile. Contact H. Khawaja Assoc. Professor, University of Tromsø, Tromsø, Norway E-mail: hassan.a.khawaja@uit.no Heat Equation Recommendation It is recommended to conduct the experiment and compare with the given results. This will validate the presented numerical modelling methodologies and the underlying principles of heat transfer. MULTIPHYSICS 2015 Numerical Results Discretized Domain (15 X 15 X 15 cm) Surface Temp. Profile (FTCS & Spectral Method – MATLAB®) Surface Temp. Profile (ANSYS®) Warm up time of 1500s Warm up time of 1500s Heat Equation ρρρρ ∂∂∂∂ ∂∂∂∂ = ?̇?q + ∂∂ ∂∂∂∂ kk ∂∂∂∂ ∂∂∂∂ where ρρ is density of the medium KKKK mm33 , ρρ is specific heat capacity JJ kkKK.KK , ?̇?q is volumetric energy generation term WW mm33 , kk is coefficient of thermal conduction WW mm.KK , ∂∂ is temperature (KK), ∂∂ refers to spatial position (mm) and ∂∂ is time (ss). Spatial expansion of Heat Equation (no source term) ∂∂∂∂ ∂∂∂∂ = αα ∂∂2∂∂ ∂∂∂∂2 + ∂∂2∂∂ ∂∂yy2 + ∂∂2∂∂ ∂∂zz2 where ∂∂, yy and zz refer to spatial positions (mm) in three dimensions and αα is the thermal diffusivity term mm 33 ss as given below, αα = kk ρρρρ Convective Boundary Condition −kk ∂∂∂∂ss ∂∂∂∂ = h(∂∂∞ − ∂∂ss) where ∂∂ss is the surface temperature (KK), ∂∂∞ is the surrounding temperature (KK) and hh is convective heat transfer coefficient WW mm22.KK


INTRODUCTION
Infrared (IR) refers to a band of electromagnetic waves between the 1 mm (frequency of 300 GHz) to 0.7 μm wavelengths (frequency of 430 THz) and photon energy from 1.24 meV to 1.7 eV as shown in the electromagnetic spectrum in Figure 1.Most of the thermal radiation from objects around 273K is in the range of IR.Thermal radiation is generated due to the interatomic motion of the particles.It happens in any matter above absolute zero (zero degrees Kelvin).Thermal radiation is emitted regardless of the physical state (solid, liquid and gas) of the matter [1].
The surface temperature can be determined using the Stefan-Boltzmann Law [2] based on the amount of thermal energy emitted from the object's surface as given in Equation (1).
Where  is heat transfer per unit time (W),  is emissivity in comparison to black body (dimensionless),  is Stefan-Boltzmann constant (W/(m 2 .K 4 )), A is area of emitting surface (m 2 ),   is surface temperature (K) and  ∞ is the room (surrounding) temperature (K).
The radiative emissivity of an object varies with the wavelength.Figure 2 shows how the emissivity of pure ice made from distilled water varies with wavelength [3].It shows that the value of ice emissivity varies from 0.965 to 0.995 in the range of 4μm to 13μm wavelengths which means that ice has high radiative emittance in the thermal and the far IR ranges.IR detection devices such as IR cameras scan over a range of wavelengths and average over the results to calculate the IR signature [4].The output of an IR device is a temperature profile superimposed over geometrical features as shown in Figure 3.In this work, three different methodologies are used to solve the 3D transient heat equation to simulate the infrared signature of an ice cube.These methodologies are the finite difference method (FDM), the spectral method, and the ANSYS® Multiphysics software.The size of the cube is 15cm x 15cm x 15cm.A cubical geometry was chosen to simplify the problem.Also, the thermal signature is symmetric in a cubical geometry, which allows for easy modelling using simulation techniques.
The results show the variation in temperature when an ice cube initially at a constant temperature of -28°C is left in room temperature conditions (25°C) to warm up.The temperature profiles on the surface of the cube are compared in the time and space domains.

METHODOLOGY
The underlying physics of heat transfer through conduction in a solid medium can be solved mathematically using the heat equation [5] as given in Equation (2).
Where  is density of the medium (kg/m 3 ),  is specific heat capacity (J/(kg K)), ̇ is the volumetric energy generation term (W/m 3 ), k is coefficient of thermal conductivity (W/ (m.K)),  is temperature (K),  refers to spatial position (m) and  is the time (s).
The extended form of the above equation in three spatial dimensions with no energy generation term [6] is given in Equation (3).
Where ,  and  refer to spatial positions (m) in three dimensions and  is the thermal diffusivity term (m 2 /s) as given in (4).
To solve Equation ( 3), the boundary, and the initial conditions are required.The convective boundary conditions [7] are applied on each external surface of the cubical geometry as given in Equation (5).
Where   is the surface temperature (K),  ∞ is the surrounding temperature (K) and ℎ is convective heat transfer coefficient (W/(m 2 .K)).
This boundary condition represents the existence of convection heating or cooling at a surface and is obtained through the energy balance at the surface.
The initial condition, in this case, is the uniform temperature throughout the cube that is set to -28°C (245K).This temperature varies as the ice cube starts to warm up.Density(), specific heat capacity () and coefficient of thermal conduction () vary with temperature [8] as shown in Figure 4.
The convective heat transfer coefficient may also vary between 4-10 W/(m 2 .K) depending on the surrounding conditions.
The range of interest in this problem is from -28°C to 0°C.At 0°C ice will start phase change from solid to liquid water with no variation in temperature.Phase change is not simulated in this study.
Physical and thermal properties of ice between -28°C to 0°C are given in Table 1.As found, the standard deviation is less than 5% of the average values hence these values can be considered as constants for setting up the simulations.The values of the constants set are given in Table 2.
The following assumptions are considered in this study, • Energy transfer from ice cube through the mode of radiation is minimal.
• Energy transfer from the ice cube to the surrounding is only through natural convection.
• Variation in physical and thermal properties with temperature are not significant.
Values of set constants are given in Table 2.

FINITE DIFFERENCE AND SPECTRAL METHODS (MATLAB®)
Finite difference method (FDM) is a numerical method for solving differential equations such as the heat equation given in Equation ( 3).This method approximates the differentials with differences by discretizing the dependent variables (temperature) in the independent variable domains (space and time) [9].Each discretized value of the dependent variable is referred to as a nodal value.In this case, heat equation given in Equation ( 3) is discretized using FDM forward-time central-space (FTCS).The discretized equation is given in Equation (6).Where superscript  and subscript ,, refer to time and position for a value of nodal temperature respectively.∆ is a timestep size (s) and ∆,∆,∆ are the differences in the spatial position of temperature nodes.
The convective boundary condition is also discretised using FDM and only applied to the outer surfaces as given in Equation (7).
It is vital for the stability and accuracy of FDM to choose the correct time step value.In this work, Courant-Friedrichs-Lewy (CFL) condition [9,10] is used to decide the time step size.CFL condition for the heat equation is given in Equation (8).
Equations ( 6) and ( 7) are solved and post-processed in MATLAB®.Results are discussed in sections 5 and 6.
The spectral method is also a numerical method for solving differential equations such as heat equation given in Equation ( 3).This method assumes the solution can be written as a sum of certain basic functions such as Fourier series or as, in this case, polynomials.The spectral method gives a lower numerical error for mathematically continuous solutions in comparison to other numerical methods [11].
In this method, space dimensions are needed to be discretized using Chebyshev points, which gives an uneven spatial grading and is better for fitting polynomials.An example of 2D Chebyshev-Lobatto points grid [12] is shown in Figure 5.This method can be extended to multi-dimensions such as heat equation as given in Equation ( 3) by adding the derivatives of polynomials.
Time stepping is essential for stability and accuracy of the numerical results.In this case, Runge-Kutta time stepping method (RK4) is used.Equation ( 3) is solved using the spectral method in MATLAB®.Results are discussed in sections 5 and 6.

THERMAL SIMULATION (ANSYS® MULTIPHYSICS)
ANSYS® Multiphysics software offers to simulate various physical phenomena [13].The method of solution is based on finite element method (FEM) [14].In this work, ANSYS® Multiphysics thermal module is used to solve thermal signature of an ice cube.To do so a cubic geometry is built in ANSYS® Multiphysics Graphics User Interface (GUI) and meshed using thermal mass Solid Brick 8 noded 278 elements [15].This mesh is tested for space and time sensitivity.The mesh built in ANSYS® Multiphysics is shown in Figure 6.

Figure 6: ANSYS® Multiphysics Mesh of 3D space
The initial condition is constant temperature throughout the mesh.The convective boundary condition is applied on three boundary surfaces.Symmetry (zero heat flux) is applied to the rest of the three boundary surfaces to reduce the run time of simulation.A program chosen algorithm is used to control the time stepping.Results are discussed in sections 5 and 6.

RESULTS AND DISCUSSION
This section discusses the results obtained by FDM, spectral method and ANSYS® Multiphysics simulations.Surface temperature contour plots obtained through FDM and Spectral methods at various time intervals are given in Figure 7.It is notable that the temperature profile is establishing at 100s.The difference between max and min values of temperature is around 3°C. Around 500s, the temperature profile is more established, and the difference between max and min values has risen to around 5°C.At the 1500s, the temperature profile is fully developed, and the gap between the maximum and minimum values of temperature is about 6°C.At 4500s, the pattern appears to be similar as noted earlier however the difference between max and min has started to drop (about 4.5°C).This difference falls slightly as surface temperature reaches close to melting temperature (0°C) and temperature throughout the cube tends to stabilize to a constant value.Since simulation does not take into account the phase change condition, therefore, results after 0°C are not valid.From all above plots, it can also be observed that the max value of temperature is at the corners, and the minimum value is in the surface centre.Surface temperature contour plots at various time intervals obtained through thermal simulation in ANSYS® Multiphysics are given in Figure 8.All indicated values are in °C.ANSYS® Multiphysics also demonstrated same behaviour as earlier discussed in FDM results.

COMPARISON OF SIMULATIONS RESULTS
Variation in the temperature of a corner of the cube is plotted against time for all three methods for a comparison as given in Figure 9.The result indicates the total time taken by an ice cube to reach melting temperature, which is found to be approximately 5500s.The comparison also shows close agreement between the three methods.

CONCLUSION & FUTURE WORK
Three different methods, namely finite difference method (FDM), spectral method and ANSYS® Multiphysics software, are used to simulate the thermal image of an ice cube, when warming from -28°C to reach melting point under room temperature conditions.Results revealed that ice cube of dimensions (15 X 15 X 15 cm) takes approximately 5500s to reach melting temperature.During this time temperature profile develops within the ice cube with a temperature difference of 5-6°C.These results are important to understanding the thermal behaviour of ice.Future work is proposed to capture the thermal image of an ice cube through a thermal imaging device (e.g.IR camera) and compare with given simulation methods.This work will also help to build the physical relation between thermal imaging devices and underlying physics of heat transfer.

ACKNOWLEDGMENT
The work reported in this paper is funded by the Norges forskningsråd, project no.195153/160 in collaboration with Faroe Petroleum.We would also like to acknowledge the support given by Prof. James Mercer at the UiT The Arctic University of Norway.

Figure 1 :Figure 2 :
Figure 1: Position of IR in an electromagnetic spectrum (CC BY-SA 3.0)

Figure 3 :
Figure 3: IR image superimposed over photographic image.Colour contours shows the false image for temperature in degree Celsius.

Figure 4 :
Figure 4: Variation in physical and thermal properties of ice with temperature.[8]

Figure 9 :
Figure 9: Comparison of corner temperature from three methods FDM, Spectral Method, and ANSYS® Multiphysics thermal simulation