Mathematical modeling of surface deformation during vitrification.

Mathematical modeling of surface deformation during cryopreservation by vitrification is presented in this study. The specific problem under consideration is of a cryoprotective agent (CPA) solution vitrifying in a vial, following previously obtained cryomacroscopy observations. A multiphysics solution is proposed in this study, combining coupled effects associated with heat transfer, fluid mechanics, and solid mechanics. Consistent with previous investigations, this study demonstrates that surface deformation is the result of material flow, which is the combined outcome of temperature gradients developed during the inward cooling process, the tendency of the material to change its volume with temperature, and the exponential increase in material viscosity with the decreasing temperature. During this process, the behavior of the CPA changes from liquid to a solid-like amorphous material, where the arrested flow in the vitrified state results in mechanical stresses. Results of this study show a good qualitative agreement of surface deformation with previously obtained experimental data, and support prior investigations to explain fracture tendencies propagating from the deformed surface. Results of this study also highlight the effect of heat convection in the CPA at the early stage of cooling.

Introduction

There is a global unmet need for tissue and organ preservation for the benefit of organ banking and transplant medicine, with a broad spectrum of applications ranging from critical care to improving the quality of life [1]. Cryopreservation via vitrification is the only alternative for long-term preservation of large specimens, where vitreous means an ice-free, solid-like state (vitreous in Latin means glassy) [2,3]. The challenges associated with cryopreservation by vitrification increase dramatically with the specimen size. Cryopreservation by vitrification of some small-scale systems, in the range of μm to mm, has been established as an effective and repeatable process [4-6]. Despite the reported successes in the cryopreservation of mammalian organs, such as rabbit kidney [2], whole sheep ovary [7], and rat hindlimb [8], cryopreservation of larger specimens call for further development of cryoprotective materials, protocols, and instrumentation to overcome challenges associated with their sheer size and unique characteristics. These challenges can be conveniently classified as relating the toxicity of cryoprotective agents (CPA), chilling and ischemic injuries, ice formation, repair and revival protocols, and thermomechanical stress [9-15] – the subject matter of this study.

Thermomechanical stress is the mechanical stress that develops due to the tendency of all materials to change their volume with temperature (thermomechanical stress and thermal stress are used interchangeably in this study). When that change in volume is constrained, thermal stress will develop, which may exceed the strength of the material, leading to structural damage in the cryopreserved material. The specific property that describes that volume-change tendency is the thermal expansion coefficient, while the resulting stress may be affected by the bulk temperature change of the material, localized temperature gradients [9], crystallization [16], structural relaxation near glass transition [17-19], and even thermal expansion differences between the specimen and its surrounding media or the container which holds it [12,20]. Structural integrity may also be affected by the natural degradation of the building blocks of the tissue [4,21], which is beyond the scope of the current study.

A notable phenomenon that affects the developing thermomechanical stress is surface deformation [18,22-25]. Such deformations result from the gradual transition of the material to glassy, subject to a non-uniform thermal history [18,24]. Unlike the development of thermal stress in the solid state, thermal stress development in a vitrifying material is associated with a fifteen orders of magnitude increase in viscosity with the decreasing temperature, from a water-like magnitude at room temperature to such a high value at glass transition temperature that the material can be considered solid in any practical time scale. The combination of changes in thermal expansion and in viscosity above the glass transition temperature causes material flow, essentially leading to deformations which are far more dramatic than those occurring solely due to thermal expansion (or contraction) in the solid state. Essentially, those large deformations are associated with residual stresses, which is the stress locked into the material at cryogenic storage. Furthermore, abrupt changes in the vitrified material surface might lead to stress concentration, which is a localized stress considerably higher than the average. Eventually, even if not resulting in structural failure, stresses built into the system during cooling may intensify the thermal stress developed during rewarming, making the material more prone to structural damage at that later stage of the cryogenic protocol [13–15, 26].

The current study aims at explaining the free surface deformation that often develops at the center of vitrifying CPAs contained in vials and cuvettes, such as those displayed in Fig. 1 [23,25,27]. This deformation has been found to be thermal history-dependent, and has been observed to cause stress concentration and even trigger fractures [24]. For this purpose, the current study uses multiphysics simulations to solve the coupled problems of heat transfer, fluid mechanics, and solid mechanics using COMSOL Multiphysics. To the best of my knowledge, this coupled problem has never been solved before in the context of cryopreservation.

Fig. 1.

Examples of surface deformation during vitrification observed (a) with polarized-light cryomacroscope in 50% DMSO at −124 °C [22], and (b) by visual inspection of 60% ethylene glycol at −130 °C [25].

Mathematical Formulation

Problem definition

Conceptually, this is a multiphysics problem, integrating concepts from the traditional engineering areas of heat transfer, fluid mechanics, and solid mechanics, as specified in more detail in the following subsections. The analysis in this study focuses on a cylindrical CPA container (a vial) of similar dimensions to the rectangular cuvette typically used in cryomacroscopy experiments [22-24,27], as illustrated in Fig. 2. The container under consideration is a cylindrical vial having an overall diameter of 17 mm, a wall thickness of 1 mm, a CPA height of 36 mm. The vial is subject to convective cooling in a controlled-rate cooler [22]: free convection under the vial cap and forced convection on all other boundaries.

Fig. 2.

Schematic presentation of the 2D and axisymmetric computational model and its boundary conditions.

Heat transfer model

The heat transfer problem is modeled by including a heat convection term in the classical heat conduction equation used in previous studies, to account for buoyancy effects due to the variation of density with temperature (the so-called natural convection effect): where ρ is the density, c is the specific heat, T is the temperature, t is the time, k is the thermal conductivity, and is the velocity field, which is the solution of the fluid mechanics problem (a bold symbol herein represents a vector field).

Convective heat transfer is assumed between the outer surface of the vial’s wall and the cooling chamber that contains the vial: where is the normal to the surface, h is the heat transfer coefficient by convection between the vial and the cooling chamber, and the index c denotes the chamber’s temperature (that is the thermal history in the controlled-rate cooler used for cryobiology applications). The boundary condition between the inner surface of the vial’s wall and the CPA is of continuity in temperature and heat flux.

Fluid mechanics model

The Navier-Stokes equation for creeping flow is used to describe the fluid mechanics problem in the CPA, where the coupling with the heat transfer problem comes about through the temperature-dependent density and viscosity: where p is the pressure, is the identity matrix, μ is the dynamic viscosity, and is the gravitational acceleration. The rate of change of density is taken into account through the continuity equation:

A zero normal stress boundary condition is assumed at the CPA-air interface, while neglecting the surface tension effect:

Finally, no-slip condition is assumed on all solid walls (i.e., no relative movement between the fluid and the solid walls at the common interface).

Note that the high surface deformation observed during the process calls for a special treatment in the finite element analysis (FEA), where elements must also be allowed to deform significantly. To address this difficulty while avoiding the need for incremental finite elements remeshing as the process progresses, the Arbitrary Lagrangian-Eulerian method (ALE) with Yeoh smoothing is applied [28]. This method uses the fluid velocity calculated from Eq. (3) to determine the top surface deformation, while the side and bottom boundaries are restricted in the radial and axial directions, respectively.

Solid mechanics model

During cooling on the path to vitrification, the CPA solution transforms from liquid to a solid-like amorphous material. During this process, the viscosity elevates exponentially by fifteen orders of magnitude, between its value at room temperature and the glass transition temperature [29]. Consistent with previous studies [11,14,15,18,19,24], the Maxwell fluid model is used with a single viscoelastic branch. The total stress rate is calculated as: where ε is the elastic strain tensor, ε is the creep strain tensor, and ε the thermal strain tensor. These strain rates are computed as: where E is Young’s modulus, ν is the Poisson ratio, σ is the stress rate, is the deviatoric stress tensor, and α is the coefficient of thermal expansion.

Finally, CPA adherence to the container walls is assume as the CPA vitrifies, which is consistent with prior experimental observations [12, 15,22,27]. While this boundary condition resembles the boundary condition of the fluid mechanics problem at the same location, it is redefined here for solving the multiphysics problem using the computation framework described below.

Note that brittle materials, such as the vitrified CPAs, tend to have a much lower strength in tension than in compression. Although never measured for glassy CPAs before to the best of my knowledge, other glassy materials, such as silica glass, exhibit at least one order of magnitude higher strength in compression than in tension [30]. In broad terms, the principal stresses are defined as the normal stresses calculated at an angle such that the shear stresses become zero. The first principal stress is the largest of the three normal stress tensor components, which is the stress displayed in the figures in this study, where a positive sign means tension.

Computation framework

The heat balance equation, Eq. (1), and the momentum equation, Eq. (3) are bidirectionally coupled and are solved simultaneously. The coupling effects are significant, as is evident from cryomacroscopy experiments [12,22,27]. Since mechanical stress at relatively low strain rates is not expected to generate a significant amount of heat compared with the external cooling and rewarming powers imposed on the system during the vitrification protocol, the thermomechanical problem is assumed to be unidirectionally coupled. Consistently, the coupled thermal-fluid problem is solved first and then its solution is used as an input to the solid mechanics problem.

The thermal expansion of the container affects the various coupled phenomena in different ways. When the material is free to flow as liquid, the container shrinkage bears minimal effect on the field of flow, compared with the free convection currents. However, thermal shrinkage of the container for example can impose compressive stress on the vitrifying CPA and, therefore, affects the solid mechanics solution. Hence, the thermal expansion of the container is taken into account in solving the solid mechanics problem but not in solving the fluids mechanics problem. Note that it is possible to use a highly compliant container, such as a cryobag, in order to minimize the thermal expansion effect of the container on the solid mechanics solution.

The computation framework has been reduced to practice using COMSOL Multiphysics software package (v5.5), with an Intel Core i7-9700 processor (8-Core, 12 MB Cache, up to 4.7 GHz). The typical computation run time for the conditions outlined below is in the range of 48 min and 5.2 h.

Material properties and thermal conditions

Consistent with prior experimental studies, the vial illustrated in Fig. 2 is made of polystyrene [22,23], with the relevant physical properties listed in Table 1. The reference CPA under consideration is 7.05 M DMSO, which has solid mechanics properties representative of high concentration CPA cocktails, such as DP6 and VS55 [31], with relevant physical properties listed in Table 1. Consistent with the previous cryomacroscopy experiments (using the Kryo 16, Planner, Inc.), the vial is placed in the chamber of a controlled-rate cooler. The air flow inside the cooling chamber is highly turbulent, where a heat transfer coefficient of 346 W/m2K from the container was previously measured [18,23].

Table 1

Material properties of 7.05 M DMSO used in the current study as a reference CPA solution and of a polystyrene vial used in the previous experimental studies (temperatures is in °C).

Parameter	Material	Value/function	Ref.
Density, kg/m3	CPA	1058 − 0.41 T	[32]
	Vial	1055 − 0.26 T	[33]
Young’s Modulus, GPa	CPA	0.8	[19]
	Vial	3.45	[17]
Poisson’s Ratio	CPA	0.25	[17]
	Vial	0.35	[17]
Thermal Conductivity, W/m-°C	CPA	0.356 + 7.42 × 10−4 T − 1.29 × 10−6 T2 − 6.87 × 10−8 T3 − 2.95 × 10−10 T4	[16]
	Vial	0.14 + 1.3 × 10−4 T	[34]
Thermal Expansion Coefficient, 1/°C	CPA	1.1 × 10−5	[32]
	Vial	8 × 10−5	[33]
Heat Capacity, W/kg-°C	CPA	2804 + 4.205 T − 0.054 T2 − 4.902 × 10−5 T3	[18]
	Vial	1121 + 3.94 T	[35]

The representative thermal protocol selected for the current study is based on a constant chamber cooling rate from an initial uniform temperature of 20 °C down to cryogenic storage at −140 °C, which is 8 °C below the glass transition temperature of the reference solution [17,24]. This study includes a range of cooling rates, between 5 °C/min, which surpasses the critical cooling rate of 7.05 M DMSO to achieve vitrification, and 50 °C/min, which is the maximum cooling rate achievable with the Kryo 16 chamber. While more elaborated thermal protocols may be conceived to avoid toxicity effects in higher temperatures and excessive residual stress at cryogenic storage, the effects of surface deformation take place at intermediate temperatures, and hence this protocol is considered adequate for the purpose of the current study.

The viscosity of the vitrifying CPA is highly temperature dependent. The current study uses compiled literature data, based on piecewise cubic interpolation [29,36]. The viscosity values range overall from 3.4 × 10−3 Pa s at 20 °C to 4.5 × 1012 Pa s at −140 °C. To facilitate convergence of the solid mechanics simulations, the viscosity above −120 °C was kept constant at a value of 1.1 × 107 Pa s, with a negligible effect on the developing stress above this temperature. By contrast, the viscosity values were taken with no modifications for the solution of the fluid mechanics problem.

Results and discussion

Representative results for the coupled thermal, fluid, and stress fields are displayed in Fig. 3, subject to a cooling rate of 20 °C/min in the cooling chamber and at three representative points in time. Furthermore, two vial materials are considered for stress analysis there, the regular polystyrene vial (Figs. 3(g)–3(i)) and a hypothetically highly compliant container (Fig. 3(j)–3(l)), which in practice has mechanical properties similar to those of the CPA. At the early stage of cooling, when the material is free to flow as a liquid, temperature gradients lead to a convection cell, which is reflected by the closed streamlines in the velocity field, in Fig. 3(d). In turn, this convection cell intensify vertical heat transfer and creates a warmer area at the top of the vial, Fig. 3(a). The artifacts of the convection cell have been observed in recorded videos of cryomacroscopy due to unintended contamination in the form of floaters [22,27].

Fig. 3.

Representative Results of the coupled temperature, velocity, and stress fields, subject to a chamber cooling rate of 20 °C/min from the initial temperature of 20 °C to a storage temperature of −140 °C: (a)-(c) temperature field at 100 s, 500 s, and 1500 s, respectively; (d)-(f) the corresponding velocity field at the same instances in time, respectively; (g)-(i) the first principle stress distribution in a vial made of polystyrene, respectively; and (j)-(l) the first principle stress distribution in a hypothetical highly compliant vial, respectively. Note that the tensile strength of the CPA is about 3.2 MPa [12] and hence the stress color scheme is presented up to this value for better clarity, while higher stress values remain red. Further note that T∞ is the cooling chamber temperature.

As the cooling process progresses, the viscosity increases and the CPA transforms into a solid-like material, as is evident by the disappearance of the convection cell in Fig. 3(e)–3(f). During this later period of the protocol, the streamlines terminate at the walls and the free surface, which indicates material contraction towards the wall like a solid, as opposed to flowing along walls in Fig. 3(d). During that later time, the warmest region changes to be at the center of the domain, which is dictated by the underlying principles of linear heat conduction, Fig. 3(b). Eventually, the material approaches thermal equilibrium, Fig. 3(c), and the rate of flow (contraction actually) drops by seven orders of magnitude, Fig. 3(f), compared with the flow velocity in the low viscous range, Fig. 3(d). Subject to significant temperature gradients during the cooling process, the material can more freely flow at the center of the vial where temperatures are higher. In the absence of flow constraints on the upper free surface, where the temperatures are higher and with the assistance of gravity as a body force, the accumulated contraction effect causes the formation of a surface cavity along the centerline of the vial–the primary focus of the current study.

The stress level during the early stage of cooling is practically zero, Fig. 3(g). Early indication of significant stress starts to develop where the free surface of the CPA connects to the wall of the vial, Fig. 3(h). Note that a tensile stress level of 3.2 MPa has been found as the fracture strength of the material in previous studies [12]. From that point on, the critical stress level affects a growing area, until the material reaches thermal equilibrium at the storage temperature, Fig. 3(i). If the material will not fail during cooling – in the form of fractures for example – that stress distribution will remain as a residual stress during cryogenic storage [14]. The residual stress is the stress locked in into the material during the process of vitrification and it is path dependent. This means that different cooling protocols will result in different stress distributions and intensities [9,11,13-15,17,24]. It is possible to alleviate these stresses by designing intermediate steps of stress relaxation [14,17], which is out of the scope of the current study. Notably, the highest stresses are observed at the free surface, Fig. 3(i), which is consistent with experimental observation on the initiation of fracturing during cryomacroscopy studies [18,22,24,27].

The container effect on stress

It is possible to reduce the overall stress level and to affect its distribution by selecting a highly compliant container, Figs. 3(j)–(l). While the container does not impose an additional mechanical constraint in such a case, the sheer volume of the container and the gradual process of fluid turning into a solid-like material under variable temperature gradients may still lead to significant stresses. The stress level in the polystyrene vial case is threefold higher than that in the case of a highly compliant container (Fig. 3(i) vs. Fig. 3(l), respectively); the stress concentration at the free surface may far exceed those levels. It is worth noting that there are several highly compliant products in the market for cryogenic storage, popularly known as cryobags, which are made of ethyl vinyl acetate (EVA) and fluorinated ethylene propylene (FEP).

To put the application of those highly compliant materials in the appropriate context, the term compliance here is used in the context of the solid mechanics problem, when stresses are significant. Despite being a highly compliant material in the solid mechanics sense, the same material serves well as a container wall material in higher temperatures, when the CPA is free to flow as a liquid and stresses are negligible in comparison.

The cooling rate effect on surface deformation

The relationship between the cooling rate and the resulting stress has been investigated previously under a small deformations condition, which is the result of thermal expansion (or contraction) without taking into account gross material flow [9,11,13-15,17,20,21,24]. Considering the potentially devastating effects of surface deformation and stress concentration, Fig. 4 displays the variation of the free surface shape with time. It can be observed from Fig. 4 that the depth of cavity at the center of the free surface increases with the cooling rate, ranging from 3.5 mm at 5 °C/min to 11.8 mm at 50 °C/min. Notably, the overall depth of cavity approaches a maximum value with the increasing cooling rate. More generally, results suggest that lower cooling rates are not only favorable for reducing the overall mechanical stress, but also for reducing localized stress concentration effects.

Fig. 4.

Surface deformation history subject to various chamber cooling rates: (a) 5 °C/min, (b) 20 °C/min, (c) 35 °C/min, and (d) 50 °C/min.

The effect of fluid flow on stress development

The development of the cavity signifies how the physical properties of the liquid change to solid-like material in different locations at different times. Consistently, the local mechanical stress also changes with time, and so is the location of the maximum stress. For example, at the beginning of cavity formation, the stress at its center is practically zero, where the material there is free to flow as liquid, while CPA has already transformed to a solid-like material along the wall. However, Fig. 3(i) demonstrates that the stress at the same point during cryogenic storage is extremely high. In fact, the location of the maximum principal stress in the domain at any given point in time is impossible to evaluate a priori. Nonetheless, in order to get some insight on how the stress in the domain is affected by accounting for surface deformation, Fig. 5 displays the thermal and stress history at the geometric center of the container.

Fig. 5.

Thermal and maximum principal stress histories at the geometric center of the CPA domain subject to a cooling rate of 20 °C/min for two representative cases: with and without accounting for material flow.

Recall that the CPA is initially cooled as a liquid, where the material is stress-free despite the temperature gradients and thermal expansion. The material changes to solid-like as the viscosity increases, where effects of temperature gradients and thermal expansion give rise to mechanical stress. Due to the gradual change of the CPA to a solid-like material, significant mechanical stress remains in the material even when it gets to thermal equilibrium at cryogenic storage [14].

In order to demonstrate the effect of modeling fluid flow on the resulting stress, Fig. 5 displays principal stress results at the geometric center of the CPA domain, for a cooling rate of 20 °C/min. It can be seen from Fig. 5 that the center of the domain cools faster when material flow is considered. For example, the temperature difference between the two cases is 19 °C at 250 s, while the same difference becomes insignificant at about 500 s. It can be seen from Fig. 4(b) that the maximum rate of cavity deepening occurs at about the same time. It can also be seen from Fig. 5 that the residual stress at the geometric center is 0.96 MPa when surface deformation is taken into account, and 1.25 MPa when it is not. Recall from Fig. 3 that the maximum stress in the domain is close to the deformed surface, while the current discussion demonstrates the effect of surface deformation on the stress well below the surface, at the geometric center.

Summary and conclusion

Mathematical modeling of surface deformation of the CPA during cryopreservation by vitrification has been presented in this study. For this purpose, the current study uses multiphysics simulations to solve the coupled problems of heat transfer, fluid mechanics, and solid mechanics using COMSOL Multiphysics. To the best of my knowledge, this coupled problem has never been solved before in the context of cryopreservation.

In general, surface deformation is the result of material flow, which is the combined outcome of temperature gradients developed during inward cooling, thermal expansion, and the exponential increase in material viscosity with the decreasing temperature. This investigation is presented as a proof of concept for the multiphysics analysis of surface deformation during vitrification, which is demonstrated by a qualitative comparison with experimental observations obtained in cryomacroscopy studies.

Results of this study display cavity formation at the center of the CPA surface, which is similar to experimental observations. The cavity depth increases with the cooling rate, where its maximum value appears to be dependent not only on the cooling protocol, but also on the container dimensions and material properties. The results also highlight heat convection effects at higher temperature, where the material is free to flow. This effect leads to a more uniform temperature distribution at the early stage of the vitrification process. Finally, the areas of maximum stress have been identified near the deformed surface, which is consistent with experimental observations of polarized light cryomacroscopy and direct cryomacroscopy observations of fracture.