# Numerical Modeling of Heat and Mass Transfer during Cryopreservation Using Interval Analysis
## 2 Methods  
### Interval Governing Equations  
Thermal processes proceeding in the axially symmetrical, heterogeneous articular cartilage sample can be described by the following interval energy equation:  
\[\overline{c}\frac{\partial\,\overline{T}(r,\,z,\,t)}{\partial\,t}=\left[\frac {1}{r}\frac{\partial}{\partial\,r}\left(\overline{\lambda}\,r\,\frac{\partial \,\overline{T}(r,\,z,\,t)}{\partial\,r}\right)\,+\,\frac{\partial}{\partial\, z}\left(\overline{\lambda}\,\,\frac{\partial\,\overline{T}(r,\,z,\,t)}{\partial\,z} \right)\right] \tag{1}\]  
where \(\overline{\lambda}=[\lambda^{-},\,\,\lambda^{+}]\) is the interval thermal conductivity, \(\overline{c}=[c^{-},\,\,c^{+}]\) is the interval volumetric specific heat, \(\overline{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:  
\[c\dot{T}(X,t)=\
abla(\lambda\
abla T(X,t)) \tag{2}\]  
where \(X\) refers to the coordinate system. The considered Equation (1) is supplemented by the boundary conditions of the 2 nd or 3 rd type (see Figure 1):  
\[\left\{ \begin{array}{ll}r=0,\,0\leq z\leq\frac{H}{2}:&\overline{\lambda} \frac{\partial\overline{T}(r,\,z,\,t)}{\partial\,r}=\overline{0}\\ 0\leq r\leq R,\,z=\frac{H}{2}:&-\,\,\overline{\lambda}\frac{\partial\overline{ T}(r,\,z,\,t)}{\partial\,z}=\overline{0}\\ r=R,\,0\leq z\leq\frac{H}{2}:&-\,\,\overline{\lambda}\frac{\partial\overline{ T}(r,\,z,\,t)}{\partial\,r}=\alpha\left(\overline{T}-T_{\text{bulk}}\right)\\ 0\leq r\leq R,\,z=0:&\overline{\lambda}\frac{\partial\overline{T}(r,\,z,\,t )}{\partial\,z}=\alpha\left(\overline{T}-T_{\text{bulk}}\right)\end{array} \right. \tag{3}\]  
and the initial condition defined as:  
\[t=0:\qquad\qquad\overline{T}(r,\,z,\,0)=T_{0} \tag{4}\]  
where \(T_{\text{bulk}}\) is the temperature of the bathing solution [57; 58]. The modeling of the cryoprotectant transport through the extracellular matrix could be written as follows:  
\[\frac{\partial\,\overline{\overline{C}_{d}}(r,\,z,\,\,t)}{\partial\,t}=\,\left[ \frac{1}{r}\frac{\partial\,}{\partial\,r}\left(\overline{D}_{d}\,r\,\frac{ \partial\,\overline{C}_{d}(r,\,z,\,\,t)}{\partial\,r}\right)+\frac{\partial}{ \partial\,z}\left(\overline{D}_{d}\,\frac{\partial\,\overline{\overline{C}_{d} }(r,\,z,\,\,t)}{\partial\,z}\right)\right] \tag{5}\]  
where \(\overline{\overline{C}}_{d}(r,\,z,\,\,t)\) is the interval cryoprotectant concentration in the extracellular matrix, while \(\overline{D}_{d}\,\) is the interval diffusion coefficient of the cryoprotectant in the extracellular matrix estimated by the Einstein-Stokes equation [46]:  
\[\overline{D}_{d}=\frac{k_{B}\,\overline{T}}{6\pi\eta\,r_{s}} \tag{6}\]  
where \(k_{B}\) is the Boltzmann constant, \(r_{s}\) is the radius of the spherical particle molecule, \(\eta\) is the dynamic viscosity. Equation (5) is based on [PERSON]'s second law:  
\[\dot{C}_{d}(X,t)=\,\
abla(D_{d}\
abla C_{d}(X,t)) \tag{7}\]  
The mathematical model (Equation (5)) should be supplemented by the boundary conditions (see Figure 1):  
\[\left\{ \begin{array}{l}r=0,\,0\leq z\leq\frac{H}{2}:\qquad\frac{\partial \overline{C}_{d}(r,\,z,\,t)}{\partial\,r}=\overline{0}\\ 0\leq r\leq R,\,z=\frac{H}{2}:\qquad-\frac{\partial\overline{C}_{d}(r,\,z,\,t )}{\partial z}=\overline{0}\\ r=R,\,0\leq z\leq\frac{H}{2}:\qquad\overline{C}_{d}(r,\,z,\,\,t)=0.9\,C_{ \text{bulk}}\\ 0\leq r\leq R,\,z=0:\qquad\overline{C}_{d}(r,\,z,\,\,t)=0.9\,C_{\text{bulk}} \end{array} \right. \tag{8}\]  
and the initial condition in the form:  
\[t=0:\qquad\overline{C}(r,\,z,\,\,0)=C_{0} \tag{9}\]  
where \(C_{\text{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.