- Email: [email protected]
Numerical modeling for freezing and cryogenic preservation for viability of biological tissue
Available online at www.sciencedirect.com
ScienceDirect Materials Today: Proceedings 5 (2018) 18823–18832
www.materialstoday.com/proceedings
ICMPC_2018
Numerical modeling for freezing and cryogenic preservation for viability of biological tissue Pallab Kundua, S.Sukumarb,S.P.Karc,* a
M.Tech Scholar, KIIT ,Bhubaneshwar, Odisha,India P.hd Scholar,KIIT , Bhubaneshwar, Odisha, India c Associate Professor, KIIT , Bhubaneshwar, Odisha, India b
Abstract A numerical model has been developed to analyze the temperature distribution inside a biological tissue under cryogenic conditions. The energy equation is discretized using finite volume method (FVM). The left and right walls of the tissue have been subjected to cryogenic temperature of nitrogen while the top and bottom walls of the tissue have been insulated. Tri-Diagonal Matrix Algorithm (TDMA) has been used to solve the algebraic equation to find temperature distribution inside the tissue. The method is validated with experimental data from the existing literature. The enthalpy-porosity method has been employed to track the moving solid-liquid interface at different time inside the specimen and the position of solid-liquid interface has plotted at different time. The position of velocity vectors have been also plotted at different time. As the convective cooling proceeds, it has been seen that the velocity current shrinks before finally seizing to exist.
© 2018 Elsevier Ltd. All rights reserved. Selection and/or Peer-review under responsibility of Materials Processing and characterization. Keywords:Finite Volume Method; Cryogenic temperature; biological tissue; Enthalpy-porosity method.
1. Introduction Cryopreservation is a technique of achieving rapid cooling of a tissue for longer preservation. The metabolic activities are arrested at sub-zero temperatures such that when the subjected tissue is gradually retrieved, the functional properties are not hampered. The cooling rate for cryopreservation is of the order of 1000’s of ºC/min [5]. The process if not carried out at optimum cooling rate, causes either extracellular or intracellular ice formation. * Corresponding author. Tel.: +919437263907 E-mail address:[email protected]
2214-7853 © 2018 Elsevier Ltd. All rights reserved. Selection and/or Peer-review under responsibility of Materials Processing and characterization.
18824
Pallab Kundu et al./ Materials Today: Proceedings 5 (2018) 18823–18832
Nomenclature TL temperature of the left wall temperature of the right wall TR initial temperature Tint specific heat capacity Cp k thermal conductivity H total enthalpy sensible heat hsen enthalpy source term Sh L length of the tissue W width of the tissue p pressure t time U resultant velocity u component of velocity in x-direction v component of velocity in y-direction coefficient in discretized algebraic equation ap Greek Symbols ρ density α thermal diffusivity λ under-relaxation parameter β volumetric thermal expansion coefficient σ Stefan-Boltzmann constant size of time step t volume of control volume V latent heat H Subscripts f freezing point P control volume n neighbouring nodes Superscripts n iteration number 0 previous time step
which is found to be functionally improper for proper preservation. Therefore the requirement arises to obtain such a cooling rate that preservation of the tissue will be functionally efficient. The temperature distribution, freezing front propagation and velocity vector distribution inside the biological tissue have been analyzed to achieve an optimum cooling rate. Cryopreservation is not a recent technique and has been in research as early as the eighteenth century [1]. It has a significant impact on organ transplantation and tissue engineering for the artificial organ implants [2]. A detailed description of the techniques and modern day developments were studied by Jouhara et.al [3]. Chua et.al [4] investigated the accurate transient temperature distribution in a solid tumour inside a complex blood vessel network by using a numerical and nano-simulation method. Ahmedikia et.al [5] used a mathematical model for the nonFourier heat transfer during phase change of a biological tissue. Dombrovsky et.al [6] used a two-temperature model to understand the transient temperature distribution in a human tissue during repetitive freezing. Nakatsuka et.al [7] carried out experiments to sequence the freezing-thawing process for curing lung tumors using a cry probe at the tip. Rubinksy et.al [8] used an irreversible thermodynamic technique to create a mathematical model to analyze the heat transfer during freezing of biological cylindrical tissue and axial blood vessel. Kandra et.al [9] has studied the interactions of biological tissue exposed to cryogenic environment and laser irradiation. Ahmadikia et.al [10] studied Fourier and non-Fourier heat conduction inside a two-dimensional tissue. An enthalpy-porosity method has been developed by Brent et.al [11] to track the moving solid-liquid interface during melting of pure Gallium. Swaminathan et.al [12] used the enthalpy approach to model an implicit fixed grid method using enthalpy as the source term. Kumar et.al [13] incorporated an enthalpy porosity approach to a biological tissue containing blood
Pallab Kundu et al./ Materials Today: Proceedings 5 (2018) 18823–18832
18825
vessel. He also took into account the non-ideal tissue property and heat from metabolism. Xu et.al [14] studied the multi-scale methods of cryopreserving cells and tissues at the macroscopic level and also the changes in cell volume and the intracellular mass transfer at microscopic level. Li et.al [15] produced an analytical solution for three dimensional heat transfers inside a biological tissue subjected to a freezing environment. Wang et.al [16] studied one dimensional heat transfer using a finite difference model to register temperature distribution and freezing time in food. Deng et.al [17] used a non-Fourier model to extract the transient temperature and thermal stress in a cryogenically preserved skin tissue. Sanz et.al [18] produced a mathematical simulation model using finite element method to determine the cooling rates in food. Studholme et.al [19] designed a compartmental model to imitate the tissue with its solute and solvents at their osmotic pressure to study the heat and mass transfer in a biological tissue. Deng et.al [20] incorporated a numerical method on dual reciprocity boundary element method (DRBEM) to determine the solid-liquid moving interface of a biological system. Kuwayama et.al [21] studied vitrification methods for the preservation of human oocytes cryogenically. Rubinsky et.al [22] studied the damage tendency of frozen liver during cryosurgery to understand their optimum cooling rate. A critical review of methods of cryopreservation and the various principles of cryobiology has been given by Karlsson et.al [23]. Kar and Rath [24] devised a fixed-grid method for a repetitive pulsed laser in a melting process. The current work deals with the formulation of a two-dimensional numerical method to study the moving solidliquid interface and the temperature distribution inside beef tissue, the properties of which have been derived from the literature [16]. The validation of the model has been done with the already existing literature [11]. Graphs have been plotted for variation of temperature, position of moving solid-liquid interface and position of velocity vector at different time inside the tissue. 2. Physical Model A two-dimensional beef tissue has been taken under consideration for the study, the property of which have been given in Table 1. The geometrical representation of the tissue has been shown in Fig. 1(a). The horizontal and the vertical length have been taken as 0.05 m each and the top and bottom walls have been insulted. The left and right walls have been maintained at 77 K which has been selected in close accordance to the literature [5].
Fig. 1(a). Geometrical representation of the tissue with boundary conditions.
18826
Pallab Kundu et al./ Materials Today: Proceedings 5 (2018) 18823–18832
Table 1. Thermo-physical properties of the beef tissue [16] Properties
Value
Units
Cp
3600
J/kg K
ρ
1040
kg/m3
∆H
188000
J/kg
K
0.51
W/m K
L
0.05
m
W
0.05
m
Tf
272.5
K
2.1. Assumptions The following assumptions have been taken into consideration for simplifying the current study. The thermo-physical properties of beef tissue have been assumed to be constant. Fluid flow has been considered to be laminar and incompressible. Tissue has been considered to be static when being frozen. Governing Equations: Energy equation: +
(
)=
( .
ℎ. ℎ
)+
ℎ
1. (a)
1. (b)
Momentum Equations: (
)
(
)
+
(
)=
( .
ℎ. ) −
+ .
+
(
)=
( .
ℎ. ) −
+
+
1. (C)
Continuity Equation: +
(
)= 0
1. (d)
Boundary Condition:
Top wall:
−
( , )
=0
2. ( a)
Pallab Kundu et al./ Materials Today: Proceedings 5 (2018) 18823–18832
( , )
Bottom wall: −
=
Left wall :
2. (b)
= 0; 0 ≤
,
=
Right wall:
=0
=
,
18827
≤
0≤
≤
2. (c)
2. (d)
3. Numerical Model Equation 1 (a) is discretized using finite volume method as described by patankar [24]. The resulting algebraic equation is solved using Tri-Diagonal-Matrix-Algorithm (TDMA). The source term in equation can be discretized as follows: =
(
+
3. (a)
)=
3. (b)
=ℎ+∆ =
3. (c)
+∆
+
3. (d)
∆
=
= −
=
−
( ∆ )= −
3. (e)
( ∆ )
3. (f)
3. (g)
∆
The enthalpy update equation is given by
(∆ )n+1 = (∆ )n +
−
(4)
T n P is the temperature of control volume at current time step and is the under relaxation factor.
18828
Pallab Kundu et al./ Materials Today: Proceedings 5 (2018) 18823–18832
Fig.2. Flow chart of the numerical solution.
4. Results and discussion The studied numerical model is validated with the results of existing literature [11].as shown in above Fig 2. Melting of pure gallium has been carried out in a rectangular cavity. The thermo-physical properties of gallium have been given in Table 2. The position of the solid-liquid interface was tracked and validated with the results from the literature. The produced results appeared to be in good agreement with the results of the literature as shown in the Fig.3. Therefore, the current study has been done using this numerical model. Table.2. Thermo-physical properties of gallium [11]
Properties
Value
Unit
Density (liquid)
6093
kg/m3
Reference Density
6095
Kg/m3
Reference temperature
29.78
Volumetric thermal expansion co-efficient
1.2 x 10
Thermal conductivity
32
Melting point
29.78
Latent heat of fusion
80160
Specific heat capacity
381.5
0 -4
C
/0 C W/m K 0
C
J/kg J/kg -3
Dynamic viscosity
1.81 x 10
Prandtl Number
2.16 x 10-2
Kg/ m s
Pallab Kundu et al./ Materials Today: Proceedings 5 (2018) 18823–18832
2 min
4 min
6 min
18829
17 min
0.060
Width (m)
0.050 0.040 0.030 0.020 NUMERICAL 0.010 0.000 0.000
0.020
0.040
0.060
0.080
Length (m) Fig. 3. Comparison of position of solid- liquid interface during melting of gallium inside a rectangular cavity with the literature [11].
(a)
(b)
(c)
Fig.4. Variation of position of solid-liquid interface with time at (a) 60 s, (b) 120 s, ( c) 1020 s.
(a)
(b)
Fig. 5. Variation of Temperature with time inside the computational domain at (a) 360 s, (b) 1020 s.
18830
Pallab Kundu et al./ Materials Today: Proceedings 5 (2018) 18823–18832
(a)
(b)
(c)
Fig. 6. Position of Velocity vectors at different time (a) 360 s , (b) 480 s, (c) 600 s. 305
Temperature (K)
300 295 290 285 280 275 270 0
10
20
30
40
50
60
Time (s)
Fig. 7. Variation of temperature at the centre of the tissue with time.
4.1. Variation of Temperature inside the tissue with time The temperature distribution inside the tissue at different time has been shown in the Fig. 5(a), (b) and (c). The variation of temperature at the center of the tissue at different time has been shown in the Fig. 7. When the tissue is subjected to the cryogenic temperature of liquid nitrogen which is 77K, its temperature reduces drastically. The tissue is cooled until all of the moisture content inside the tissue becomes amorphous solid. Vitrification is the method of preservation of living tissues at such a temperature that all of the moisture content inside the tissue gets converted into amorphous solid. A slower cooling rate of the tissue can leave a portion of the moisture content inside the tissue unfrozen and occurrence of even a slightest amount of moisture content in the preserved tissue can cause damage to the tissue. Even in case of a very rapid cooling rate, the moisture inside the tissue gets rapidly converted into crystalline solid which is also harmful for tissue preservation. Therefore achieving an optimum cooling rate is necessary for proper preservation. In the current study, the tissue has been cooled from both the sides simultaneously to achieve a proper cooling rate which is neither very fast nor very slow.
Pallab Kundu et al./ Materials Today: Proceedings 5 (2018) 18823–18832
18831
4.2. Variation of position of the solid-liquid interface with time The movement of the solid-liquid interface at different time during solidification has been shown in the Fig.4 (a), (b) and (c). As achieving an optimum cooling rate is the main objective of the current study, the tissue is cooled from the left and right sides simultaneously. Because of cooling from both the sides, the freezing time of the tissue decreases and achievement of vitrification becomes possible. The solid-liquid interfaces move from both the sides towards the center of the tissue until whole of the tissue is vitrified. 4.3. Variation of velocity vector with time The velocity vectors have been shown in Fig. 6 (a), (b) and (c). The tissue is initially kept at 300K and it is brought in contact with liquid nitrogen at its cryogenic temperature of 77K for the cooling purpose. When the cooling begins, the strength of the convective current inside the tissue is the greatest because of a large existing temperature gradient. As the cooling is performed from both the sides, two convective currents are generated. The left one is anti-clockwise while the right one is clockwise moving. Both the convective currents help in quickening the rate of heat transfer. With the advancement in cooling time, the portions of the tissue towards the cooling walls get solidified and temperature of the still liquid portion also reduces further. This reduces the temperature gradient and the strength of the convective current with time and both the convective currents collapse when the tissue is completely vitrified. 5. Conclusion An optimum cooling rate is badly necessary to achieve the moisture contents inside the tissue in the form of amorphous ice in the extracellular as well as intracellular sites. A suitable cooling rate is further required to avoid crystallization of ice that causes damage to the functional properties of the tissue. The following conclusions have been derived from the present work. 1. A suitable cooling rate, which is neither too quick nor too slow, is required to achieve moisture content of the tissue in the form of amorphous ice in the extracellular and intracellular sites which is ideal for vitrification and cryopreservation. 2. Cooling the tissue from both the sides simultaneously, helps in achieving a suitable cooling rate for vitrification as compared to cooling from only one side. 3. A greater temperature gradient between the cold walls and the interior of the tissue at the beginning of cooling process, causes formation of convective current of higher strengths and with the advancement in cooling time, the strengths of the currents get diminished until the collapse of both when complete solidification is achieved. The current work suggests the ways to achieve a suitable cooling rate for vitrification and cryopreservation which can be helpful for preservation of the biological tissues, but the present work is based on parabolic heat transfer model. It has been found out that hyperbolic heat transfer occurs inside the biological tissues. So, the current study can be done by using the hyperbolic heat transfer model, which will further increase the level of accuracy of the results obtained. References [1] B.Luyet, G.Rapatz,A review of basic researches on the cryopreservation of red blood cells, Third Cryopreservation Conference during the Annual Meeting of the Society for Cryobiology. Washington. D. C, August 7 to 9. 1967, Cryobiology, 6, No, 5, 1970 [2] Robert M. Nerem, Ph.d. and Athanassios Sambanis ,Tissue Engineering: From Biology to Biological substitutes, Ph.d 1995. [3] Theodora K. Nannou, Jon Trembley, Jens Herrmann, Hussam Jouhara, Cryopreservation: Methods, Equipment and Critical Concerns, International IIR Conference of Cryogenics and Refrigeration Technology, Bucharest, 2016 [4] K.J.Chua, Fundamental experiments and numerical investigation of cryo-freezing incorporating vascular network with enhanced nanofreezing, International Journal of Thermal Sciences 70 (2013), 17-3. [5] H. Ahmadikia, A.Moradi, Non-Fourier phase change heat transfer in biological tissuesduring solidification, Heat Mass Transfer (2012) 48:1559–1568 [6] Leonid A. Dombrovsky , Natalia B. Nenarokomova , Dmitry I. Tsiganov , Yuri A. Zeigarnik, Modeling of repeating freezing of biological tissues and analysis of possible microwave monitoring of local regions of thawing, International Journal of Heat and Mass Transfer 89 (2015) 894–902.
18832
Pallab Kundu et al./ Materials Today: Proceedings 5 (2018) 18823–18832
[7] Seishi Nakatsuka, Hideki Yashiro, Masanori Inoue, Sachio Kuribayashi, Masafumi Kawamura, Yotaro Izumi, Norimasa Tsukada, Yoshikane Yamauchi , Kohei Hashimoto , Kansei Iwata , Taisuke Nagasawa, Yi-Shan Lin, On freeze-thaw sequence of vital organ of assuming the cryoablation for malignant lung tumors by using cryoprobe as heat source, Cryobiology 61 (2010) 317–326. [8] B. Rubinsky, D. E. Pegg, A mathematical model for the freezing process in biological tissue, Proc. R. Soc. Lond. B 234, 343-358 (1988). [9] Deepak Kandra, Tissue interactions with lasers and liquid nitrogen: a novel cryopreservation method, A Thesis Submitted to the Graduate Faculty of the Louisiana State University and Agricultural and Mechanical College In partial fulfillment of the requirements for the degree of Master of Science in Mechanical Engineering in The Department of Mechanical Engineering. [10] Hossein Askarizadeh, Hossein Ahmadikia, Analytical study on the transient heating of a two-dimensional skin tissue using parabolic and hyperbolic bioheat transfer equations, Applied Mathematical Modelling 39 (2015) 3704–3720. [11] A. D. Brent, V. R. Voller, K. J. Reid, Enthalpy-porosity Technique for modeling convection-diffusion phase change: Application to a melting of a pure metal, Numerical Heat Transfer, vol. 13, pp. 297-318, 1988. [12] C.R.Swaminathan, V.R.Voller, A general enthalpy method for modelling solidification processes, Metallurgical Transactions, (23 B). [13] Sushil Kumar, V.K.Katiyar, A Mathematical modelling in freezing and thawing process in tissues: a porous media approach, International Journal of Applied Mechanics, Vol. 2, No. 3 (2010) 617–633. [14] Feng Xu, Sangjun Moon, Xiaohui Zhang, Lei Shao, Young Seok Song and Utkan Demirci, Multi-scale heat and mass transfer modelling of cell and tissue cryopreservation, Philosophical transactions of the Royal society, 368, 561-583. [15] Fang-fang Li, Jing Liu, Kai Yue, Exact analytical solution to three-dimensional phase change heat transfer problems in biological tissues subject to freezing, Appl. Math. Mech. -Engl. Ed.30 (1), 63–72 (2009). [16] Zhengfu Wang, Han Wu, Guanghua Zhao, Xiaojun Liao, Fang Chen, Jihong Wu, Xiaosong Hu, One-dimensional finite-difference modeling on temperature history and freezing time of individual food, Journal of Food Engineering 79 (2007) 502–510. [17] Zhong-Shan Deng, Jing Liu, Non-Fourier heat conduction effect on Prediction of temperature transients and Thermal stresses in skin cryopreservation, Journal of Thermal Stresses, 26: 779–798, 2003. [18] P.D. Sanza, C. de Elvira, M. Martino, N. Zaritzky, L. Otero, J.A. Carrasco, Freezing rate simulation as an aid to reducing crystallization damage in foods, Meat Science 52 (1999) 275-278. [19] Christopher Val Studholme, Modeling Heat and Mass Transport In biological tissues during freezing, thesis. [20] Zhong-Shan Deng, Jing Liu, Modeling of multidimensional freezing problem during cryosurgery by the dual reciprocity boundary element method, Engineering Analysis with Boundary Elements 28 (2004) 97–108. [21] Masashige Kuwayama, Gábor Vajta , Osamu Kato , Stanley P Leibo, Highly efficient vitrification method for cryopreservation of human oocytes, Vol 11. No 3. 2005 300-308. [22] B. Rubinsky, C. Y. Lee, J. Bastacky, G. Onik, The Process of Freezing and the Mechanism of Damage during Hepatic Cryosurgery, Vol. 234, No.1276 (Aug. 23, 1988) [23] Jens O.M. Karlsson ,Mehmet Toner, Long-term storage of tissues by cryopreservation: critical issues, Biomaterials 1996, Vol. 17 No. 3. [24] Satya Prakash Kar, Prasanjeet Rath, A fixed-grid based mixture model for pulsed laser phase change process,Computational Thermal Sciences: An International Journal, Volume 6, 2014 Issue 1.
ScienceDirect Materials Today: Proceedings 5 (2018) 18823–18832
www.materialstoday.com/proceedings
ICMPC_2018
Numerical modeling for freezing and cryogenic preservation for viability of biological tissue Pallab Kundua, S.Sukumarb,S.P.Karc,* a
M.Tech Scholar, KIIT ,Bhubaneshwar, Odisha,India P.hd Scholar,KIIT , Bhubaneshwar, Odisha, India c Associate Professor, KIIT , Bhubaneshwar, Odisha, India b
Abstract A numerical model has been developed to analyze the temperature distribution inside a biological tissue under cryogenic conditions. The energy equation is discretized using finite volume method (FVM). The left and right walls of the tissue have been subjected to cryogenic temperature of nitrogen while the top and bottom walls of the tissue have been insulated. Tri-Diagonal Matrix Algorithm (TDMA) has been used to solve the algebraic equation to find temperature distribution inside the tissue. The method is validated with experimental data from the existing literature. The enthalpy-porosity method has been employed to track the moving solid-liquid interface at different time inside the specimen and the position of solid-liquid interface has plotted at different time. The position of velocity vectors have been also plotted at different time. As the convective cooling proceeds, it has been seen that the velocity current shrinks before finally seizing to exist.
© 2018 Elsevier Ltd. All rights reserved. Selection and/or Peer-review under responsibility of Materials Processing and characterization. Keywords:Finite Volume Method; Cryogenic temperature; biological tissue; Enthalpy-porosity method.
1. Introduction Cryopreservation is a technique of achieving rapid cooling of a tissue for longer preservation. The metabolic activities are arrested at sub-zero temperatures such that when the subjected tissue is gradually retrieved, the functional properties are not hampered. The cooling rate for cryopreservation is of the order of 1000’s of ºC/min [5]. The process if not carried out at optimum cooling rate, causes either extracellular or intracellular ice formation. * Corresponding author. Tel.: +919437263907 E-mail address:[email protected]
2214-7853 © 2018 Elsevier Ltd. All rights reserved. Selection and/or Peer-review under responsibility of Materials Processing and characterization.
18824
Pallab Kundu et al./ Materials Today: Proceedings 5 (2018) 18823–18832
Nomenclature TL temperature of the left wall temperature of the right wall TR initial temperature Tint specific heat capacity Cp k thermal conductivity H total enthalpy sensible heat hsen enthalpy source term Sh L length of the tissue W width of the tissue p pressure t time U resultant velocity u component of velocity in x-direction v component of velocity in y-direction coefficient in discretized algebraic equation ap Greek Symbols ρ density α thermal diffusivity λ under-relaxation parameter β volumetric thermal expansion coefficient σ Stefan-Boltzmann constant size of time step t volume of control volume V latent heat H Subscripts f freezing point P control volume n neighbouring nodes Superscripts n iteration number 0 previous time step
which is found to be functionally improper for proper preservation. Therefore the requirement arises to obtain such a cooling rate that preservation of the tissue will be functionally efficient. The temperature distribution, freezing front propagation and velocity vector distribution inside the biological tissue have been analyzed to achieve an optimum cooling rate. Cryopreservation is not a recent technique and has been in research as early as the eighteenth century [1]. It has a significant impact on organ transplantation and tissue engineering for the artificial organ implants [2]. A detailed description of the techniques and modern day developments were studied by Jouhara et.al [3]. Chua et.al [4] investigated the accurate transient temperature distribution in a solid tumour inside a complex blood vessel network by using a numerical and nano-simulation method. Ahmedikia et.al [5] used a mathematical model for the nonFourier heat transfer during phase change of a biological tissue. Dombrovsky et.al [6] used a two-temperature model to understand the transient temperature distribution in a human tissue during repetitive freezing. Nakatsuka et.al [7] carried out experiments to sequence the freezing-thawing process for curing lung tumors using a cry probe at the tip. Rubinksy et.al [8] used an irreversible thermodynamic technique to create a mathematical model to analyze the heat transfer during freezing of biological cylindrical tissue and axial blood vessel. Kandra et.al [9] has studied the interactions of biological tissue exposed to cryogenic environment and laser irradiation. Ahmadikia et.al [10] studied Fourier and non-Fourier heat conduction inside a two-dimensional tissue. An enthalpy-porosity method has been developed by Brent et.al [11] to track the moving solid-liquid interface during melting of pure Gallium. Swaminathan et.al [12] used the enthalpy approach to model an implicit fixed grid method using enthalpy as the source term. Kumar et.al [13] incorporated an enthalpy porosity approach to a biological tissue containing blood
Pallab Kundu et al./ Materials Today: Proceedings 5 (2018) 18823–18832
18825
vessel. He also took into account the non-ideal tissue property and heat from metabolism. Xu et.al [14] studied the multi-scale methods of cryopreserving cells and tissues at the macroscopic level and also the changes in cell volume and the intracellular mass transfer at microscopic level. Li et.al [15] produced an analytical solution for three dimensional heat transfers inside a biological tissue subjected to a freezing environment. Wang et.al [16] studied one dimensional heat transfer using a finite difference model to register temperature distribution and freezing time in food. Deng et.al [17] used a non-Fourier model to extract the transient temperature and thermal stress in a cryogenically preserved skin tissue. Sanz et.al [18] produced a mathematical simulation model using finite element method to determine the cooling rates in food. Studholme et.al [19] designed a compartmental model to imitate the tissue with its solute and solvents at their osmotic pressure to study the heat and mass transfer in a biological tissue. Deng et.al [20] incorporated a numerical method on dual reciprocity boundary element method (DRBEM) to determine the solid-liquid moving interface of a biological system. Kuwayama et.al [21] studied vitrification methods for the preservation of human oocytes cryogenically. Rubinsky et.al [22] studied the damage tendency of frozen liver during cryosurgery to understand their optimum cooling rate. A critical review of methods of cryopreservation and the various principles of cryobiology has been given by Karlsson et.al [23]. Kar and Rath [24] devised a fixed-grid method for a repetitive pulsed laser in a melting process. The current work deals with the formulation of a two-dimensional numerical method to study the moving solidliquid interface and the temperature distribution inside beef tissue, the properties of which have been derived from the literature [16]. The validation of the model has been done with the already existing literature [11]. Graphs have been plotted for variation of temperature, position of moving solid-liquid interface and position of velocity vector at different time inside the tissue. 2. Physical Model A two-dimensional beef tissue has been taken under consideration for the study, the property of which have been given in Table 1. The geometrical representation of the tissue has been shown in Fig. 1(a). The horizontal and the vertical length have been taken as 0.05 m each and the top and bottom walls have been insulted. The left and right walls have been maintained at 77 K which has been selected in close accordance to the literature [5].
Fig. 1(a). Geometrical representation of the tissue with boundary conditions.
18826
Pallab Kundu et al./ Materials Today: Proceedings 5 (2018) 18823–18832
Table 1. Thermo-physical properties of the beef tissue [16] Properties
Value
Units
Cp
3600
J/kg K
ρ
1040
kg/m3
∆H
188000
J/kg
K
0.51
W/m K
L
0.05
m
W
0.05
m
Tf
272.5
K
2.1. Assumptions The following assumptions have been taken into consideration for simplifying the current study. The thermo-physical properties of beef tissue have been assumed to be constant. Fluid flow has been considered to be laminar and incompressible. Tissue has been considered to be static when being frozen. Governing Equations: Energy equation: +
(
)=
( .
ℎ. ℎ
)+
ℎ
1. (a)
1. (b)
Momentum Equations: (
)
(
)
+
(
)=
( .
ℎ. ) −
+ .
+
(
)=
( .
ℎ. ) −
+
+
1. (C)
Continuity Equation: +
(
)= 0
1. (d)
Boundary Condition:
Top wall:
−
( , )
=0
2. ( a)
Pallab Kundu et al./ Materials Today: Proceedings 5 (2018) 18823–18832
( , )
Bottom wall: −
=
Left wall :
2. (b)
= 0; 0 ≤
,
=
Right wall:
=0
=
,
18827
≤
0≤
≤
2. (c)
2. (d)
3. Numerical Model Equation 1 (a) is discretized using finite volume method as described by patankar [24]. The resulting algebraic equation is solved using Tri-Diagonal-Matrix-Algorithm (TDMA). The source term in equation can be discretized as follows: =
(
+
3. (a)
)=
3. (b)
=ℎ+∆ =
3. (c)
+∆
+
3. (d)
∆
=
= −
=
−
( ∆ )= −
3. (e)
( ∆ )
3. (f)
3. (g)
∆
The enthalpy update equation is given by
(∆ )n+1 = (∆ )n +
−
(4)
T n P is the temperature of control volume at current time step and is the under relaxation factor.
18828
Pallab Kundu et al./ Materials Today: Proceedings 5 (2018) 18823–18832
Fig.2. Flow chart of the numerical solution.
4. Results and discussion The studied numerical model is validated with the results of existing literature [11].as shown in above Fig 2. Melting of pure gallium has been carried out in a rectangular cavity. The thermo-physical properties of gallium have been given in Table 2. The position of the solid-liquid interface was tracked and validated with the results from the literature. The produced results appeared to be in good agreement with the results of the literature as shown in the Fig.3. Therefore, the current study has been done using this numerical model. Table.2. Thermo-physical properties of gallium [11]
Properties
Value
Unit
Density (liquid)
6093
kg/m3
Reference Density
6095
Kg/m3
Reference temperature
29.78
Volumetric thermal expansion co-efficient
1.2 x 10
Thermal conductivity
32
Melting point
29.78
Latent heat of fusion
80160
Specific heat capacity
381.5
0 -4
C
/0 C W/m K 0
C
J/kg J/kg -3
Dynamic viscosity
1.81 x 10
Prandtl Number
2.16 x 10-2
Kg/ m s
Pallab Kundu et al./ Materials Today: Proceedings 5 (2018) 18823–18832
2 min
4 min
6 min
18829
17 min
0.060
Width (m)
0.050 0.040 0.030 0.020 NUMERICAL 0.010 0.000 0.000
0.020
0.040
0.060
0.080
Length (m) Fig. 3. Comparison of position of solid- liquid interface during melting of gallium inside a rectangular cavity with the literature [11].
(a)
(b)
(c)
Fig.4. Variation of position of solid-liquid interface with time at (a) 60 s, (b) 120 s, ( c) 1020 s.
(a)
(b)
Fig. 5. Variation of Temperature with time inside the computational domain at (a) 360 s, (b) 1020 s.
18830
Pallab Kundu et al./ Materials Today: Proceedings 5 (2018) 18823–18832
(a)
(b)
(c)
Fig. 6. Position of Velocity vectors at different time (a) 360 s , (b) 480 s, (c) 600 s. 305
Temperature (K)
300 295 290 285 280 275 270 0
10
20
30
40
50
60
Time (s)
Fig. 7. Variation of temperature at the centre of the tissue with time.
4.1. Variation of Temperature inside the tissue with time The temperature distribution inside the tissue at different time has been shown in the Fig. 5(a), (b) and (c). The variation of temperature at the center of the tissue at different time has been shown in the Fig. 7. When the tissue is subjected to the cryogenic temperature of liquid nitrogen which is 77K, its temperature reduces drastically. The tissue is cooled until all of the moisture content inside the tissue becomes amorphous solid. Vitrification is the method of preservation of living tissues at such a temperature that all of the moisture content inside the tissue gets converted into amorphous solid. A slower cooling rate of the tissue can leave a portion of the moisture content inside the tissue unfrozen and occurrence of even a slightest amount of moisture content in the preserved tissue can cause damage to the tissue. Even in case of a very rapid cooling rate, the moisture inside the tissue gets rapidly converted into crystalline solid which is also harmful for tissue preservation. Therefore achieving an optimum cooling rate is necessary for proper preservation. In the current study, the tissue has been cooled from both the sides simultaneously to achieve a proper cooling rate which is neither very fast nor very slow.
Pallab Kundu et al./ Materials Today: Proceedings 5 (2018) 18823–18832
18831
4.2. Variation of position of the solid-liquid interface with time The movement of the solid-liquid interface at different time during solidification has been shown in the Fig.4 (a), (b) and (c). As achieving an optimum cooling rate is the main objective of the current study, the tissue is cooled from the left and right sides simultaneously. Because of cooling from both the sides, the freezing time of the tissue decreases and achievement of vitrification becomes possible. The solid-liquid interfaces move from both the sides towards the center of the tissue until whole of the tissue is vitrified. 4.3. Variation of velocity vector with time The velocity vectors have been shown in Fig. 6 (a), (b) and (c). The tissue is initially kept at 300K and it is brought in contact with liquid nitrogen at its cryogenic temperature of 77K for the cooling purpose. When the cooling begins, the strength of the convective current inside the tissue is the greatest because of a large existing temperature gradient. As the cooling is performed from both the sides, two convective currents are generated. The left one is anti-clockwise while the right one is clockwise moving. Both the convective currents help in quickening the rate of heat transfer. With the advancement in cooling time, the portions of the tissue towards the cooling walls get solidified and temperature of the still liquid portion also reduces further. This reduces the temperature gradient and the strength of the convective current with time and both the convective currents collapse when the tissue is completely vitrified. 5. Conclusion An optimum cooling rate is badly necessary to achieve the moisture contents inside the tissue in the form of amorphous ice in the extracellular as well as intracellular sites. A suitable cooling rate is further required to avoid crystallization of ice that causes damage to the functional properties of the tissue. The following conclusions have been derived from the present work. 1. A suitable cooling rate, which is neither too quick nor too slow, is required to achieve moisture content of the tissue in the form of amorphous ice in the extracellular and intracellular sites which is ideal for vitrification and cryopreservation. 2. Cooling the tissue from both the sides simultaneously, helps in achieving a suitable cooling rate for vitrification as compared to cooling from only one side. 3. A greater temperature gradient between the cold walls and the interior of the tissue at the beginning of cooling process, causes formation of convective current of higher strengths and with the advancement in cooling time, the strengths of the currents get diminished until the collapse of both when complete solidification is achieved. The current work suggests the ways to achieve a suitable cooling rate for vitrification and cryopreservation which can be helpful for preservation of the biological tissues, but the present work is based on parabolic heat transfer model. It has been found out that hyperbolic heat transfer occurs inside the biological tissues. So, the current study can be done by using the hyperbolic heat transfer model, which will further increase the level of accuracy of the results obtained. References [1] B.Luyet, G.Rapatz,A review of basic researches on the cryopreservation of red blood cells, Third Cryopreservation Conference during the Annual Meeting of the Society for Cryobiology. Washington. D. C, August 7 to 9. 1967, Cryobiology, 6, No, 5, 1970 [2] Robert M. Nerem, Ph.d. and Athanassios Sambanis ,Tissue Engineering: From Biology to Biological substitutes, Ph.d 1995. [3] Theodora K. Nannou, Jon Trembley, Jens Herrmann, Hussam Jouhara, Cryopreservation: Methods, Equipment and Critical Concerns, International IIR Conference of Cryogenics and Refrigeration Technology, Bucharest, 2016 [4] K.J.Chua, Fundamental experiments and numerical investigation of cryo-freezing incorporating vascular network with enhanced nanofreezing, International Journal of Thermal Sciences 70 (2013), 17-3. [5] H. Ahmadikia, A.Moradi, Non-Fourier phase change heat transfer in biological tissuesduring solidification, Heat Mass Transfer (2012) 48:1559–1568 [6] Leonid A. Dombrovsky , Natalia B. Nenarokomova , Dmitry I. Tsiganov , Yuri A. Zeigarnik, Modeling of repeating freezing of biological tissues and analysis of possible microwave monitoring of local regions of thawing, International Journal of Heat and Mass Transfer 89 (2015) 894–902.
18832
Pallab Kundu et al./ Materials Today: Proceedings 5 (2018) 18823–18832
[7] Seishi Nakatsuka, Hideki Yashiro, Masanori Inoue, Sachio Kuribayashi, Masafumi Kawamura, Yotaro Izumi, Norimasa Tsukada, Yoshikane Yamauchi , Kohei Hashimoto , Kansei Iwata , Taisuke Nagasawa, Yi-Shan Lin, On freeze-thaw sequence of vital organ of assuming the cryoablation for malignant lung tumors by using cryoprobe as heat source, Cryobiology 61 (2010) 317–326. [8] B. Rubinsky, D. E. Pegg, A mathematical model for the freezing process in biological tissue, Proc. R. Soc. Lond. B 234, 343-358 (1988). [9] Deepak Kandra, Tissue interactions with lasers and liquid nitrogen: a novel cryopreservation method, A Thesis Submitted to the Graduate Faculty of the Louisiana State University and Agricultural and Mechanical College In partial fulfillment of the requirements for the degree of Master of Science in Mechanical Engineering in The Department of Mechanical Engineering. [10] Hossein Askarizadeh, Hossein Ahmadikia, Analytical study on the transient heating of a two-dimensional skin tissue using parabolic and hyperbolic bioheat transfer equations, Applied Mathematical Modelling 39 (2015) 3704–3720. [11] A. D. Brent, V. R. Voller, K. J. Reid, Enthalpy-porosity Technique for modeling convection-diffusion phase change: Application to a melting of a pure metal, Numerical Heat Transfer, vol. 13, pp. 297-318, 1988. [12] C.R.Swaminathan, V.R.Voller, A general enthalpy method for modelling solidification processes, Metallurgical Transactions, (23 B). [13] Sushil Kumar, V.K.Katiyar, A Mathematical modelling in freezing and thawing process in tissues: a porous media approach, International Journal of Applied Mechanics, Vol. 2, No. 3 (2010) 617–633. [14] Feng Xu, Sangjun Moon, Xiaohui Zhang, Lei Shao, Young Seok Song and Utkan Demirci, Multi-scale heat and mass transfer modelling of cell and tissue cryopreservation, Philosophical transactions of the Royal society, 368, 561-583. [15] Fang-fang Li, Jing Liu, Kai Yue, Exact analytical solution to three-dimensional phase change heat transfer problems in biological tissues subject to freezing, Appl. Math. Mech. -Engl. Ed.30 (1), 63–72 (2009). [16] Zhengfu Wang, Han Wu, Guanghua Zhao, Xiaojun Liao, Fang Chen, Jihong Wu, Xiaosong Hu, One-dimensional finite-difference modeling on temperature history and freezing time of individual food, Journal of Food Engineering 79 (2007) 502–510. [17] Zhong-Shan Deng, Jing Liu, Non-Fourier heat conduction effect on Prediction of temperature transients and Thermal stresses in skin cryopreservation, Journal of Thermal Stresses, 26: 779–798, 2003. [18] P.D. Sanza, C. de Elvira, M. Martino, N. Zaritzky, L. Otero, J.A. Carrasco, Freezing rate simulation as an aid to reducing crystallization damage in foods, Meat Science 52 (1999) 275-278. [19] Christopher Val Studholme, Modeling Heat and Mass Transport In biological tissues during freezing, thesis. [20] Zhong-Shan Deng, Jing Liu, Modeling of multidimensional freezing problem during cryosurgery by the dual reciprocity boundary element method, Engineering Analysis with Boundary Elements 28 (2004) 97–108. [21] Masashige Kuwayama, Gábor Vajta , Osamu Kato , Stanley P Leibo, Highly efficient vitrification method for cryopreservation of human oocytes, Vol 11. No 3. 2005 300-308. [22] B. Rubinsky, C. Y. Lee, J. Bastacky, G. Onik, The Process of Freezing and the Mechanism of Damage during Hepatic Cryosurgery, Vol. 234, No.1276 (Aug. 23, 1988) [23] Jens O.M. Karlsson ,Mehmet Toner, Long-term storage of tissues by cryopreservation: critical issues, Biomaterials 1996, Vol. 17 No. 3. [24] Satya Prakash Kar, Prasanjeet Rath, A fixed-grid based mixture model for pulsed laser phase change process,Computational Thermal Sciences: An International Journal, Volume 6, 2014 Issue 1.