Numerical Modeling of Heat and Mass Transfer during Cryopreservation Using Interval Analysis

: In the paper, the numerical analysis of heat and mass transfer proceeding in an axially symmetrical articular cartilage sample subjected to the cryopreservation process is presented. In particular, a two-dimensional (axially symmetrical) model with imprecisely deﬁned parameters is considered. The base of the heat transfer model is given by the interval Fourier equation and supplemented by initial boundary conditions. The phenomenon of cryoprotectant transport (Me 2 SO) through the extracellular matrix is described by the interval mass transfer equation. The liquidus-tracking (LT) method is used to control the temperature, which avoids the formation of ice regardless of the cooling and warming rates. In the LT process, the temperature decreases/increases gradually during addition/removal of the cryoprotectant, and the articular cartilage remains on or above the liquidus line so that no ice forms, independent of the cooling/warming rate. The discussed problem is solved using the interval ﬁnite difference method with the rules of directed interval arithmetic. Examples of numerical computations are presented in the ﬁnal part of the paper. The obtained results of the numerical simulation are compared with the experimental results, realized for deterministically deﬁned parameters.


Introduction
Cryopreservation is a process which aims to preserve tissues or cells at low temperatures (between −80 and −196 • C). In this process, the biological activity of tissues is reduced or completely stopped and then the physiological temperature is restored again. Successful cryopreservation does not significantly impact the basic functions of cryopreserved biological tissues or cells, for example, mechanical properties [1,2].
Cryopreservation of cells or tissues has many medical applications, such as cryobanking human cells or organ transplantation. It allows the transport of organs between different medical centers on time. The long-term storage of stem cells, for example, is the initial step toward tissue engineering, which will enable the treatment of incurable diseases and the regeneration of soft tissues [3]. Cryopreservation of articular cartilage can be helpful in its complicated treatment, considering the damage to the articular cartilage, which can lead to osteoarthritis [4][5][6][7][8][9][10].
Heat and mass transfer problems are usually solved by assuming that the equations appearing in the mathematical model and all parameters in those equations are deterministic. Such assumption does not produce an exact image of the thermal and mass processes met in the cryopreservation process. It seems more reasonable to consider uncertainties related, for example, to the material parameters. Such analysis is particularly important in biomechanics because thermophysical parameters can change in a wide range (thermal conductivity, volumetric specific heat). This stems from the parameters depending on numerous individual traits such as age, sex, occupation, etc. [57,58].
In the present paper, a numerical analysis of the cryopreservation process, including heat and mass transfer proceeding, is presented. In the model, the LT method protocol for a cryopreserved articular cartilage sample was simulated. Furthermore, some thermophysical parameters cannot be determined deterministically; therefore, the values were approximated experimentally. This often causes the parameters to be inaccurate. Because of that, the interval arithmetic is used to properly interpret uncertain thermophysical parameters [57,58]. In this research, the thermal proceeding was modeled using the interval Fourier equation, whereas the cryoprotectant transport through the extracellular matrix was described by the interval mass transfer equation.

Interval Governing Equations
Thermal processes proceeding in the axially symmetrical, heterogeneous articular cartilage sample can be described by the following interval energy equation: where λ = [λ − , λ + ] is the interval thermal conductivity, c = [c − , c + ] is the interval volumetric specific heat, T(r, z, t) is the interval temperature, t is the time, r and z denote the cylindrical coordinates. Equation (1) is based on the Fourier equation: where X refers to the coordinate system. The considered Equation (1) is supplemented by the boundary conditions of the 2nd or 3rd type (see Figure 1): and the initial condition defined as: where T bulk is the temperature of the bathing solution [57,58].
and the initial condition in the form: where bulk C is the cryoprotectant concentration in the bathing solution [28]. The mathematical model does not consider the phase change phenomenon. This is due to the LT method being used to model heat and mass transfer. The LT protocol regulates the temperature and concentration so that the temperature of the sample is above or on the liquidus line, which eliminates the probability of ice crystallization in cells; see the calculated eutectic temperatures and the melting points of tissues in [21]. The modeling of the cryoprotectant transport through the extracellular matrix could be written as follows: where C d (r, z, t) is the interval cryoprotectant concentration in the extracellular matrix, while D d is the interval diffusion coefficient of the cryoprotectant in the extracellular matrix estimated by the Einstein-Stokes equation [46]: where k B is the Boltzmann constant, r s is the radius of the spherical particle molecule, η is the dynamic viscosity. Equation (5) is based on Fick's second law: The mathematical model (Equation (5)) should be supplemented by the boundary conditions (see Figure 1): (8) and the initial condition in the form: where C bulk is the cryoprotectant concentration in the bathing solution [28]. The mathematical model does not consider the phase change phenomenon. This is due to the LT method being used to model heat and mass transfer. The LT protocol regulates the temperature and concentration so that the temperature of the sample is above or on the liquidus line, which eliminates the probability of ice crystallization in cells; see the calculated eutectic temperatures and the melting points of tissues in [21].

Numerical Model
The numerical model of thermal and mass processes proceeding in the domain of the articular cartilage sample is based on the finite difference method (FDM) in the version presented in [59,60]. In this paper, the model is extended to the case of interval parameters. At first, the time grid was introduced with a constant step ∆t: The sample domain was covered by a regular geometrical mesh and the five-point stars created by the central node (i, j) and the adjoining ones ( Figure 2). The boundary nodes were located at 0.5 h or 0.5 k with respect to the real boundary (h, k are the steps of the regular mesh in directions r and z, respectively). This approach provides a better approximation of the Neumann and Robin boundary conditions [59,60].

Numerical Model
The numerical model of thermal and mass processes proceeding in the domain of the articular cartilage sample is based on the finite difference method (FDM) in the version presented in [59,60]. In this paper, the model is extended to the case of interval parameters.
At first, the time grid was introduced with a constant step Δt: The sample domain was covered by a regular geometrical mesh and the five-point stars created by the central node (i, j) and the adjoining ones ( Figure 2). The boundary nodes were located at 0.5 h or 0.5 k with respect to the real boundary (h, k are the steps of the regular mesh in directions r and z, respectively). This approach provides a better approximation of the Neumann and Robin boundary conditions [59,60]. The approximation of differential operators of the interval heat transfer Equation (1) using the mean quotient was applied: (12) where the interval thermal conductivities were assumed in the form of harmonic means of values corresponding to the basic nodes: 22 , (13) and the thermal resistances between the central node and the adjoining ones (the rules concerning the mathematical operations defined for interval numbers should be considered, of course) were the following: Using the mean differential quotient once again, the following formula was obtained: where , ij r is the radial coordinate of the node (i, j). The same approach was used for the second term on the right-hand side of Equation (1): The approximation of differential operators of the interval heat transfer Equation (1) using the mean quotient was applied: where the interval thermal conductivities were assumed in the form of harmonic means of values corresponding to the basic nodes: and the thermal resistances between the central node and the adjoining ones (the rules concerning the mathematical operations defined for interval numbers should be considered, of course) were the following: Using the mean differential quotient once again, the following formula was obtained: where r i,j is the radial coordinate of the node (i, j).
The same approach was used for the second term on the right-hand side of Equation (1): and after appropriate transformations, the following formula was obtained: where: The approximate form of the interval energy Equation (1) for the internal nodes (i, j) and the transition t f−1 → t f was the following: where: are the shape functions of the differential mesh. As visible, the explicit scheme of the interval FDM was applied here.
In a similar way, the derivatives appearing in the interval mass Equation (5) were approximated. Thus: where the interval diffusion coefficients were assumed in the form of harmonic means of values corresponding to the basic nodes: and the diffusion resistances between the central node and the adjoining ones were the following: Using the mean differential quotient again, the following formula was obtained: The same approach was used for the second term on the right-hand side of Equation (5), and after transformations, we obtained: while: The final approximate form of the interval mass Equation (5) for the internal nodes (i, j) and the transition t f−1 → t f was the following: To summarize, for each transition t f−1 → t f , f = 1, 2, . . . , F, based on the Equation (20) the interval temperature T f i,j is calculated and then based on the Equation (29) the interval cryoprotectant concentration C d f i,j is determined. The method of attaching boundary conditions is discussed in detail in [60]. The system of Equations (20) and (29) was solved using the assumption of the stability condition for the explicit differential scheme [59,60].
All mathematical operations leading to the determination of the temperature field and the concentration field corresponding to time level t f were performed according to the rules of directed interval arithmetic (see Appendix A). For example, using Equation (29), interval subtraction and interval division were performed (see Equations (A8) and (A9)), while in Equation (20), an interval multiplication operation was applied additionally (see Equation (A6)).
The experiment using the LT approach was carried as follows: The cylindrical sample of articular cartilage was immersed in a bathing solution whose temperature and concentration were regulated using computer control. The cryoprotectant was delivered to the sample in such a manner that prevented ice crystals forming without causing toxicity. This is the sample temperature being close to the liquidus temperature during the process (the melting point of the sample was changed by the cryoprotectant's impact) [21,22].

Results
In this work, we used the LT protocol reported by Pegg et al. [21]. It consisted of eight steps in the cooling phase (the temperatures in steps 1 and 2 were the same; therefore, it can be considered that the cooling phase for the temperature consisted of seven steps) and seven steps in the heating phase. In each step, which lasted a certain amount of time, the temperature and concentration of the bathing solution were kept at the same level. For example, step 1 in the cooling phase lasted 10 min, and the temperature as well as concentration of the bathing solution were equal to 22 • C and 10 % (w/w), respectively. Similarly, in step 1 in the heating phase, the temperature and concentration of the bathing solution were equal to −48.5 • C and 63 % (w/w), respectively. Such conditions in this step were maintained for 30 min. The temperature of the bathing solution, T bulk , and the concentration of the bathing solution, C bulk , regulated in accordance with the LT method presented in [21,28] for the remaining steps are shown in Table 1 and in Figure 3.
In this work, we used the LT protocol reported by Pegg et al. [21]. It consisted of eight steps in the cooling phase (the temperatures in steps 1 and 2 were the same; therefore, it can be considered that the cooling phase for the temperature consisted of seven steps) and seven steps in the heating phase. In each step, which lasted a certain amount of time, the temperature and concentration of the bathing solution were kept at the same level. For example, step 1 in the cooling phase lasted 10 min, and the temperature as well as concentration of the bathing solution were equal to 22 °C and 10 % (w/w), respectively. Similarly, in step 1 in the heating phase, the temperature and concentration of the bathing solution were equal to −48.5 °C and 63 % (w/w), respectively. Such conditions in this step were maintained for 30 min. The temperature of the bathing solution, Tbulk, and the concentration of the bathing solution, Cbulk, regulated in accordance with the LT method presented in [21,28] for the remaining steps are shown in Table 1   Cooling In the mathematical model, the bathing solution parameters refer to boundary conditions. On the basis of the determined boundary conditions, it was possible to calculate In the mathematical model, the bathing solution parameters refer to boundary conditions. On the basis of the determined boundary conditions, it was possible to calculate the interval temperature and interval concentration in the given nodes inside the sample using the interval heat transfer equation and the interval mass transfer equation.
The computations were performed using the interval FDM, with the assumption that the time step ∆t = 0.001 s, mesh steps are h = 0.0001 m, and k = 0.00005 m. Table 1 also presents a comparison of temperatures and the Me 2 SO concentrations obtained in the simulation, with the temperature set according to the stepwise LT protocol of Pegg et al. [21,28]. The obtained results refer to the point with the coordinates r = 0.1 mm and z = 0.45 mm.
Firstly, the simulation without using interval arithmetic rules for constant values of the parameters c and λ was performed. The results of these computations are provided in          Figure 4 presents the history of the temperature over time for transition from 22 to −5 • C from step 2 to step 3 in the cooling phase ( Figure 4a) and transition from −5 to 22 • C from step 6 to step 7 in the heating phase (Figure 4b). The diagram shows that the required temperature is reached within 20 s. Figure 5 shows the history of the concentration over time for transition from 22 to −5 • C from step 2 to step 3 in the cooling phase ( Figure 5a) and transition from −5 to 22 • C from step 6 to step 7 in the heating phase (Figure 5b).
Afterward, calculations for a model using interval arithmetic rules were performed. Figure 6 illustrates the distribution of temperatures T − (Figure 6a) and T + (Figure 6b) as well as the concentrations C d − (Figure 6c) and C d + (Figure 6d) for step 3 after 10 s in the cooling phase.
Appl. Sci. 2021, 11, x FOR PEER REVIEW 10 of 16 Figure 5 shows the history of the concentration over time for transition from 22 to −5 °C from step 2 to step 3 in the cooling phase ( Figure 5a) and transition from −5 to 22 °C from step 6 to step 7 in the heating phase (Figure 5b).
Afterward, calculations for a model using interval arithmetic rules were performed. Figure 6 illustrates the distribution of temperatures T − (Figure 6a) and T + (Figure 6b Figure 7 shows the cooling curves for the given point at the moment of transition from 22 to −5 °C from step 2 to step 3 in the cooling phase (Figure 7a) and the heating curves at the moment of transition from −5 to 22 °C from step 6 to step 7 in the heating phase (Figure 7b).  The numerical simulation achieved very narrow temperature intervals. Additionally, the temperature value obtained from the experiment was within the calculated temperature interval ( Table 1).
The obtained concentration intervals increased over time in the case of cooling, while during heating, they narrowed again, which can be seen by analyzing the history of the concentration over time presented in Figure 8. Comparing the obtained concentration intervals with the results of experimental studies carried out by Pegg et al. [21], one can see that the values determined by means of numerical simulation did not coincide with the results measured by the experimental method (Table 1). In turn, the concentration values presented in [28] are more similar to the experimental results determined by Pegg et al. On this basis, we concluded that the resulting simulation error was caused by using the Einstein-Stokes model to estimate the diffusion coefficient (compared with the diffusion coefficient in [28]).

Discussion
In this paper, the numerical analysis of heat and mass transfer proceeding in an axially symmetrical articular cartilage sample subjected to the cryopreservation process was presented. The problem discussed was solved using the interval FDM with the rules of The numerical simulation achieved very narrow temperature intervals. Additionally, the temperature value obtained from the experiment was within the calculated temperature interval ( Table 1).
The obtained concentration intervals increased over time in the case of cooling, while during heating, they narrowed again, which can be seen by analyzing the history of the concentration over time presented in Figure 8. Comparing the obtained concentration intervals with the results of experimental studies carried out by Pegg et al. [21], one can see that the values determined by means of numerical simulation did not coincide with the results measured by the experimental method (Table 1). In turn, the concentration values presented in [28] are more similar to the experimental results determined by Pegg et al. On this basis, we concluded that the resulting simulation error was caused by using the Einstein-Stokes model to estimate the diffusion coefficient (compared with the diffusion coefficient in [28]). The numerical simulation achieved very narrow temperature intervals. Additionally, the temperature value obtained from the experiment was within the calculated temperature interval ( Table 1).
The obtained concentration intervals increased over time in the case of cooling, while during heating, they narrowed again, which can be seen by analyzing the history of the concentration over time presented in Figure 8. Comparing the obtained concentration intervals with the results of experimental studies carried out by Pegg et al. [21], one can see that the values determined by means of numerical simulation did not coincide with the results measured by the experimental method (Table 1). In turn, the concentration values presented in [28] are more similar to the experimental results determined by Pegg et al. On this basis, we concluded that the resulting simulation error was caused by using the Einstein-Stokes model to estimate the diffusion coefficient (compared with the diffusion coefficient in [28]).

Discussion
In this paper, the numerical analysis of heat and mass transfer proceeding in an axially symmetrical articular cartilage sample subjected to the cryopreservation process was presented. The problem discussed was solved using the interval FDM with the rules of

Discussion
In this paper, the numerical analysis of heat and mass transfer proceeding in an axially symmetrical articular cartilage sample subjected to the cryopreservation process was presented. The problem discussed was solved using the interval FDM with the rules of directed interval arithmetic. The mathematical model is based on Fourier and mass diffusion interval equations, assuming interval values of thermophysical parameters.
This approach allows the mathematical model to include imprecise values that are closer to reality, as these parameters are determined by experimental means. The use of interval values allows obtaining solutions in the form of intervals, which better reflect the analyzed process. The obtained interval temperatures and concentrations were compared with the results of the experiment and [21]. While satisfactory results were obtained for temperatures, different results were obtained for concentrations than in the experiment, probably because of estimation of the diffusion coefficient by the Einstein-Stokes model differently than in [28].

Conclusions
In this paper, numerical analysis of the cryopreservation process, including heat and mass transfer proceeding, was considered. The interval version of the FDM allows finding a numerical solution in interval form, and such information may be important especially for parameters that are difficult to estimate, for example, tissue parameters.
The obtained results of the numerical simulation were compared with the results of the experiment realized for deterministically defined parameters, such as thermal conductivity and volumetric specific heat. Considering imprecise parameter values in the mathematical model allows obtaining temperatures and cryoprotectant concentrations in the form of intervals, which, in turn, allows modeling the cryopreservation phenomenon closer to reality.
Numerical modeling of the cryopreservation process is an opportunity for further development in medicine as it allows a better understanding of the physicochemical phenomena occurring in tissues and cells.
The next stage of the work will be to extend the presented mass transfer model with the phenomenon of cryoprotectant transport (Me 2 SO) across the cell membrane.  The left or the right endpoint of interval a can be denoted as a s , s ∈ {+, −}, where s is a binary variable. This variable can be expressed as a product of two binary variables and is defined as: An interval is called proper if a − < a + , improper if a − > a + and degenerate if a − = a + . The set of all directed interval numbers can be described as D = P ∪ I, where P denotes a set of all directed proper intervals and I denotes a set of all improper intervals.
Additionally a subset Z = Z P ∪ Z I ∈ D should be defined, where Z P = {a ∈ P| a − ≤ 0 ≤ a + } and Z I = {a ∈ I| a + ≤ 0 ≤ a − }.
For directed interval numbers, two binary variables are defined. The first of them is the direction variable: and the other one is the sign variable: In the set of directed interval numbers, the basic mathematical operations for a = [a − , a + ] and b = [b − , b + ] belonging to D can be defined as follows: [min(a − · b + , a + · b − ), max(a − · b − , a + · b + )], a, b ∈ Z P [max(a − · b − , a + · b + ), min(a − · b + , a + · b − )], a, b ∈ Z I 0, a ∈ Z P , b ∈ Z I ∪ a ∈ Z I , b ∈ Z P which now makes it possible to obtain the number 0 by subtracting two of the same intervals and the number 1 by dividing the interval by itself, which is impossible using classic interval arithmetic [62,63].