Numerical research of solidification dynamics with anisotropy and thermal fluctuations

The influence of thermal fluctuations and anisotropy on the solidification process of a pure component is analyzed. It helps to understand the unstable freezing process where complicated structures such as dendrites could be formed due to a supercooling event reached during the cooling of the biological material. The study consists of mathematical modeling, validation with an analytical solution, and shows the influence of thermal noises on phase front dynamics. The analysis has been modeled in the framework of the Phase Field Method with Cahn-Hilliard formulation of a free energy functional [1]. The phase front is described by the Eulerian approach of fluid fields and formulated as a Phase Field scalar variable (order parameter) with a predefined, diffused boundary thickness. The results describe an influence scale onto directional phase front propagation dynamics, and how significant are stochastic thermal noises in micro-scale freezing.


Introduction
Cell/tissue protection is becoming an integral part in many current medical fields like a transplantation, preservation of female fertility, or tissue engineering [1]. The principle of cryopreservation is based on lowering the metabolic activity. The biological materials chemical reactions, which sum is metabolism, could be reduced by decreasing it's temperature [1]. In the paper [2] the author estimated that for normothermic animals, metabolic activity decreases 1,5 to 2 times for every 10 K decrease in temperature. Thus, the lower the temperature is, the longer material could be stored.
It was shown in [1], when choosing the temperature below 143 K (−130 o C) material viability is remained over the years [3]. Unfortunately, cryopreservation methodology has to deal with several undesirable events, which could occur during the freezing and thawing process. This includes osmotic stresses, recrystallization of the ice, cryoprotectant agents toxicity or Intracellular Ice Formation (IIF). Among all of those listed events, IIF was found as a most dangerous for preserving cells and tissues [4]. It stands out as a formation of the ice inside a cell, and causes damage of organelles. The IIF occurs when a cell contains unbound (i.e. freezable) supercooled water, and the actual temperature of the cell interior is below the nucleation point [5]. To protect the cell interior before undesirable ice nucleation, the freezing process has to be performed in a way to avoid this event. One well known method of freezing biological material is the (so-called) slow freezing method, where the extracellular icefront is propagated relatively slowly to dehydrate the cell interior before reaching the nucleate temperature of water [6].
Many protocols exist concerning the slow freezing method. Some of them were completed after years of experiments. Alternatively, the research has been going in the direction of numerical simulations since 1984, where the first mathematical description of cell dehydration was proposed by Mazur [7]. In this paper, Mazur developed relative volume water content of the cell, taking as a parameter the cooling rate of freezing, permeability of the membrane and the actual temperature of the cell interior. By this simple method, it was possible to determine whether cooling avoids IIF or not. Recient mathematical descriptions of the cell freezing process deal with many variables aspects concerning freezing itself. The model presented by Jaeger et al. [8,9] is a two-dimensional finite element approach. The main emphasis was placed to the osmotic response of the cell membrane to an uneven flux of water during ice freezing propagation. This work presents an effect of the external ice propagation front, which changes the local concentration of the salt due to solidification phase movement. For modeling a complicated system of the phases an Arbitrary Lagrangian-Eulerian method (ALE) has been chosen, where an numerical mesh is adapted to the movment of the biological cell boundary. Kinetics of the water was modeled as a response of unequal salt concentration difference between extracellular and intracellular space. Later mathematical models were mainly focused on unstable freezing ice patterns due to the high subcooling of the solution. An example could also be cited, the two consecutive works of Udaykumar and Mao [10] and Yang and Udaykumar [11], in which numerical calculations of solidification itself were presented in a binary mixture solution (water-salt). The research was carried out on FEM fixed cartesian mesh by the immersed phase boundary method with the supply of level-set method. As the results show, the dendritic pattern of the freezing is shown due to instabilities caused by large subcooling.
Thermal instabilities together with the anisotropic solidification process is the main objective analized in the article. To investigate the influence of the dendritic pattern on the directional freezing front, the mathematical model based on Phase Field Method (PFM) has been proposed. This approach allows for the detection of the solidification pattern on the basis of primary free energy minimization in the numerical domain. This part is a background of the wider project concerning numerical research of solidification front influence on biological cell structure.

Analysis and modelling
Considering only the moving boundary problem two basic approaches have to be pointed out: Eulerian and Lagrangian.
Lagrangian methods maintain the interface as discontinuity and explicitly follow the evolution of tracking mesh points. In the opposite extreme, Eulerian methods deals with fixed grid and implicit volume fraction tracking. Unquestionable advantage of using the Eulerian approach is related to the formulation of the problems involving free boundaries, where deformation of the free surface is so large that Lagrangian methods cannot be used [12]. From the Eulerian tracking methods, one, the PFM was choosen for further analysis. PFM formulation was introduced by Cahn and Hilliard in 1958 [13] and its principles were based on Ginzburg and Landau free energy functional [14]. The evolution of the interface is represented by implicitly given phase-field variable called order parameter (η) [14] and has diffusive characteristics across the phase boundary. The main idea of diffuse interface is associated with the free energy minimization [14]. Thanks to the diffusiveinterface approach, PFM could track unstable and complex solid-liquid topologies without the necessity of explicit implementation of boundary conditions on the phase border [15]. Thus, all variables are smooth through the two-phase medium and remain continuous over the numerical domain. The basic form of PFM governing equation is shown below: The order parameter usually represents microscopic variables, like magnetic spins, but in transition between solid and liquid, the order parameter has a macroscopic character and is associated with density difference. Where solid/liquid phase transition is considered, each phase in the system is represented by the finite value of the phasefield parameter, in that case; η = 1 corresponds to solid, and η = 0 to liquid phase. In PFM the free energy functional of the single component liquid/solid transition is described as: In literature, this representation is called a Cahn-Hilliard free energy functional. The RHS parts of the above presented formulation describe the actual local density of free energy a ( , ) in homogenous system. The related gradient term is responsible for the diffusive character of the interface. The local free energy density of a pure component is assumed to be: where two additional functions, the double-well function ( ) = (1 − ) , and interpolating function ℎ( ) = (6 − 15 + 10) are implemented directly. Double-well function of that form ensuring the stable minimum / = 0 for = 1 and 0, where the second, ℎ( ) function provides smooth interpolation of the local free energy difference between two phases. The final form of the PFM equation has been presented in the following form: The local free energy difference ∆ ( ) is represented by Gibbs free energy per unit volume [16].
As it can be seen, local free energy density depends on the latent heat of fusion and actual temperature level. This consideration allows for the tracking of the phase front implicitly without additional implementation of the melting temperature on the border. It also consistents with the Gibbs-Thomson equation at vanishing interface thickness -the sharp-interface limit condition. The temperature itself provides energy field potential which establishes the direction of the phase boundary force. For the temperature equal to the relation above is equal to zero, where for < local free energy difference is positive, and for > obtains a negative value. Temperature field is derived from the first law of thermodynamics, and is represented by: An additional term related to latent heat of fusion was implemented as a source term interpolated by function ℎ( ) which comes from PFM. It was assumed, that the Peclet number in cell freezing is Pe << 1. For such a low value of Pe number, motion of the fluid inside the described systems is transported mainly due to diffusion without any convection effects.
Morphologically complex structures begin to grow until constrained by surface tension and interface kinetic effects. The anisotropy of the surface free energy plays a significant role in the dynamics of dendritic growth [17]. The first approach to describe the instability phenomena during solidification was presented by Mullins and Sekerka [10,18]. Interfacial energy anisotropy in the model was introduced based on the assumption that scalar gradient coefficient depends on the interface orientation [19]. Function of that coefficient in 2-D Cartesian space has the following form [19,20]: where = tan ( / ) is the angle between the normal to the interface and reference axis. Subscripts in , represent the order parameter derivatives at x and y axis, and the is a isotropic constant coefficient. It was shown by Kobayashi [14], that increasing the anisotropy magnitude parameter ( ) significantly influences the solid structure morphology. Roughly speaking, the patters for higher values of anisotropy magnitude result on more individual dendrites formation, with weak side-branching. It should be noticed, that the mechanisms of side-branching in the crystal form are mainly governed by the thermal fluctuation occurring near the dendrite tip. These stochastic noise forces are added into the model, to consider microstructure evolution changes [15]. The stochastic noise is included into the main governing equation. Final form of the order-parameter is equal to: The random "r" number is simulated on C++ language by a pseudo random number generator. The described mathematical formulation has been implemented by numerical packages OpenFOAM 2.3.0 [21].

Validation of the model
Validation was performed in two steps. First of all, the stable directional solidification process was compared with the one-dimensional "classic" Stefan problem formulation. The analytical solution of that problem was solved for pure water using the von Neumann approach without any thermal fluctuations or initial subcooling. The results are shown in Fig. 1. Secondly, the mathemathical model for instability of the freezing pattern was compared with a pure component by the directional freezing of pure metal (nickel) cast [19]. The obtained PFM results are presented in Fig. 2. The comparison data pattern could be found in the above listed article. In both cases, the developed numerical description of a unnstable solidification front showed very good agreement with the analytical solution of the Stefan problem, or dendritic growth pattern for the nickel cast. Some deviation in the matter of one-dimensional Stefan solution colud be explained due to the diffusive character of the solidification front itself in PFM. In the Stefan problem, solidification front is assumed to be a sharp-edge.

Results and discussion
The morphology changes during the unstable solidification of pure water were checked by this method. In case of pure water, the 6 th -fold symmetry of anisotropy was assumed [14]. The nucleation seed of the radius 2.5 • 10 was inserted in the middle of bulk liquid water which represets a solid sphere of a pure substance immersed into an undercooled melt. The numerical domain was equal to 2 2 . The single cell grid size was assumed to be about Δ = 2.5 • 10 which results in 640 000 mesh cells in total. Initial water temperature inside the numerical domain was set to the uniform value of 245.84 K. It corresponds to undercooling of Δ = 27.3 . The pure water properties used in the numerical calculations are presented in Tab. 1. The result of the simulation is studied in Fig. 3.
In Fig. 3 the snowflake like dendritic structure can be observed. Very characteristic is the sharp, visible dendrites pattern which follows anisotropy directions, with additional secondary branching from the main direction of anisotropy. This is a characteristic breakdown of the solid-liquid interface due to the Mullins-Sekerka instability [10,22]. It is worth mentioning, that dendrite growth is only due to thermal fluctuations, which is an invaluable advantage of taking into consideration free energy based on PFM. No initial sinusoidal perturbation on the phase boundary was applied. The temperature field propagates together with the solid-liquid interface. The very visiblie relation between order-parameter and temperature iso-lines dependency is shown in Fig. 4. The melting point temperature for pure water is set exactly on the interface border, where order-parameter achieve a value close to 1 (see Fig. 3 for comparison), which corresponds to 100% solid representation. Fig.3. Pure water unstable solidification with anisotropy and thermal noise -order parameter, the simulation time is equal to 1e-7 s.
The value of order-parameter between 0 and 1 represents the region which is called "the mushy zone"a mixture of solid and liquid. Therefore, for that region the temperature of the melt is lower than the temperature for pure ice, taking into consideration present boundary condition. Depending on the direction of dendrites growth, through time, the temperature field forms a strong thermal gradient. This feature could have an impact to cell freezing protocols like a red blood cell, which size is approx. 6-8m in diameter [21]. To show thermal gradients of temperature profile, Fig. 5 represents temperature distribution profile over probe lines. Those proble lines were placed horizontaly -in the direction of main dendrite structure gowth; and diagonally -in the directional approx. 45deg. inclined to the previous line (in between the main dendrites directions). The zero point for both probes was set in the center of the snow-like frozen water. As a result, it can be seen, that there are differences in temperature propagation over time The fastest temperature diffusion is in the direction of main dendrite growth orientation (horizontal probe). The temperature stip gradient profile reaches the position of 0.6 m after 1e-7 s for diagonal position, where the same temperature jump in Fig. 5 is reached after 7e-8 s for the main dendrite tip (horizontaly positioned line).

Conclusions
The new numerical approach for the directional freezing of biological cells, by using the PFM, has been presented and validated. Our interest is in the computation of dendritic solidification used in cryopreservation of cells and tissue. There is strong evidence, that the fluctuations and thermal noises have an impact on the solidifiction process under subcooling conditions. The temperature profile, as well as the direction of the solid front depends on anisotropy freezing and thermal fluctuations. Moreover, the resultant temperature field shows a strong relationship with the direction of solidifying instabilities during the time of solidification. The primary goal of this paper is to present a numerical technique for capturing a diffusive solid interface during the growth process. The microsegregation of solid interface is a usually an ignored issue due to the simplification of the cryopreservation models. The present model is a primiary step to develop a more adequate mathematical description of slow freezing methodology.