# Numerical Modeling of Heat and Mass Transfer during Cryopreservation Using Interval Analysis
## 2 Methods
### 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. Figure 1: Domain considered and boundary conditions. At first, the time grid was introduced with a constant step \(\Delta t\):  
\[t^{0}<t^{1}<\ldots<t^{f-2}<t^{f-1}<t^{f}<\ldots<t^{F}<\infty \tag{10}\]  
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:  
\[\left(r\,\overline{\lambda}_{\frac{1}{\partial\tau}}^{\overline{q}\overline{ T}}\right)_{i,j+0.5}^{f-1}=r_{i,j+0.5}\,\overline{\lambda}_{i,j+0.5}^{f-1}\, \frac{\overline{T}_{i,j+1}^{f-1}-\overline{T}_{i,j}^{f-1}}{\overline{h}}= \left(r_{i,j}+\frac{h}{2}\right)\frac{\overline{T}_{i,j+1}^{f-1}-\overline{T }_{i,j}^{f-1}}{\overline{R}_{i,j+1}^{f-1}} \tag{11}\]  
\[\left(r\,\overline{\lambda}_{\frac{1}{\partial\tau}}^{\overline{q}\overline{ T}}\right)_{i,j-0.5}^{f-1}=r_{i,j-0.5}\,\overline{\lambda}_{i,j-0.5}^{f-1}\, \frac{\overline{T}_{i,j}^{f-1}-\overline{T}_{i,j-1}^{f-1}}{\overline{h}}= \left(r_{i,j}-\frac{h}{2}\right)\frac{\overline{T}_{i,j}^{f-1}-\overline{T}_{ i,j-1}^{f-1}}{\overline{R}_{i,j-1}^{f-1}} \tag{12}\]  
where the interval thermal conductivities were assumed in the form of harmonic means of values corresponding to the basic nodes:  
\[\overline{\lambda}_{i,j+0.5}^{f-1}=\frac{2\,\overline{\lambda}_{i,j}^{f-1} \,\overline{\lambda}_{i,j+1}^{f-1}}{\overline{\lambda}_{i,j+1}^{f-1}+\overline {\lambda}_{i,j+1}^{f-1}},\qquad\overline{\lambda}_{i,j-0.5}^{f-1}=\frac{2\, \overline{\lambda}_{i,j}^{f-1}\,\overline{\lambda}_{i,j-1}^{f-1}}{\overline{ \lambda}_{i,j}^{f-1}+\overline{\lambda}_{i,j-1}^{f-1}} \tag{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:  
\[\overline{R}_{i,j+1}^{f-1}=\frac{0.5\,h}{\overline{\lambda}_{i,j}^{f-1}}+ \frac{0.5\,h}{\overline{\lambda}_{i,j+1}^{f-1}},\qquad\qquad\overline{R}_{i,j -1}^{f-1}=\frac{0.5\,h}{\overline{\lambda}_{i,j}^{f-1}}+\frac{0.5\,h}{ \overline{\lambda}_{i,j-1}^{f-1}} \tag{14}\]  
Using the mean differential quotient once again, the following formula was obtained:  
\[\left[\frac{1}{r}\frac{\partial}{\partial r}\left(r\,\overline{\lambda}\frac{ \partial\overline{T}}{\partial r}\right)\right]_{i,j}^{f-1}=\frac{1}{r_{i,j} \,h}\left[\left(r_{i,j}+\frac{h}{2}\right)\frac{\overline{T}_{i,j+1}^{f-1}- \overline{T}_{i,j}^{f-1}}{\overline{R}_{i,j+1}^{f-1}}+\left(r_{i,j}-\frac{h}{2 }\right)\frac{\overline{T}_{i,j-1}^{f-1}-\overline{T}_{i,j}^{f-1}}{\overline{ R}_{i,j-1}^{f-1}}\right] \tag{15}\]  
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):  
\[\left(\overline{\lambda}\frac{\partial\overline{T}}{\partial z}\right)_{i+0.5, j}^{f-1}=\overline{\lambda}_{i+0.5,j}^{f-1}\,\frac{\overline{T}_{i+1,j}^{f-1}- \overline{T}_{i,j}^{f-1}}{k}=\frac{\overline{T}_{i+1,j}^{f-1}-\overline{T}_{i, j}^{f-1}}{\overline{R}_{i+1,j}^{f-1}} \tag{16}\]  
Figure 2: Five-point star.