Numerical thermal simulation of cryoexposure using Ansys

The current advances and problems of medical cryoexposure thermal simulation are considered. Recommendations on improving geometric models, thermophysical properties, boundary conditions, and for parameterizing of the computational domain are proposed. These recommendations can be applied to most modern FEA software packages (tested in Ansys CFX 19.2). Examples of cryoablation and cryotherapy simulation are presented.


INTRODUCTION
Two directions in the development of medical cryoexposure thermal simulation computer programs have been formed. The first is the most detailed numerical simulation of cryoexposure and on its basis, the second is the development of express planning programs. The prospect of using modern software to solve such problems [1][2][3][4][5] is based on increasing of the calculation power, which allows to make more and more complex calculations using FEA, including the use of increasingly sophisticated geometric models (finite element meshes consist of millions of elements). FEA software are becoming more sophisticated every year: turbulence models are being refined, the speed and quality of mesh generators are being improved, parallelization processes are being optimized, that is especially important in the context of increasing of the number of CPU cores. Some previously used approaches and simplifications become irrelevant and require revision. Assumptions made earlier due to limited computing power now can be avoided. However, the accuracy of simulation is made up of two components: the accuracy of the calculation and the accuracy of the input data. Therefore, increasing only the quality of software does not lead to a significant qualitative leap in the quality of simulation. It is necessary to refine the input data -the geometric models, thermal properties, boundary conditions, and other physical processes simulation features that affect cryoexposure [6]. This paper presents the features of increasing the accuracy of input data in such simulation programs.

Geometrical model
Computer technologies for supporting surgical operations (CAS) allows the use of data obtained from medical imaging (CT, MRI, etc.) to prepare the operation [7,8]. As a standard format for medical images obtained from imaging equipment, DICOM (Digital Imaging and Communications in Medicine) is used, which is a grid of color or tone pixels in the cut plane. Using a set of sequential slices (images), it is recommended to obtain a 3D image, as well as segment target organs. The set of slices with the selected target area form a segmented part, which is the surface that limits this area ( Figure 1). Image quality depends on the resolution of the medical equipment. However, a spatial raster model obtained by computed tomography cannot be directly used for FEA. Modern systems of geometric modeling (CAD) use solid-state models (geometric primitives, i.e. the basic set of geometric shapes that underlie all graphical constructions), segmentation results are in a surface model (for example, STL format -flat triangles). Development of reliable algorithms for automatic or semi-automatic transition from STL to solid-state model is not a completely solved problem that limits the possibilities of using DICOM images without simplifying the forms of biological tissues. Today, the development of medical imaging methods is aimed, among other things, at creating tools for automatic segmentation of biological tissues [9], which provides prospects for the using of geometric models that are closest to the real forms of various biological tissues. а) b) Figure 1. Obtaining the surface of the geometrical model using the example of kidney segmentation: a) segmentation of the target area in a flat section; b) resulting volumetric surface of the kidney

Thermal properties
Generally, the properties of biological tissues (density, heat capacity, thermal conductivity, heat of perfusion and metabolism) are determined by constants or empirical piecewise linear functions [10][11][12], which is mainly effective above 0 °C. Also, the greatest difficulty is the modeling of phase transitions (melting and boiling) of water in biological tissues and cryoagents of cryoequioment (sources of exposure). In these cases, a jump in enthalpy (most often as effective heat capacity), temperatures of the beginning and end of the phase transition, and mass fractions of ice and bound water are added. Moreover, during a phase transition, the properties of biological tissues also depend on the rate of temperature change. For example, various forms of intracellular ice produced by fast and slow freezing significantly affect the properties of this moisture-containing material. To increase the modeling accuracy, it is necessary to use the most accurate data on the thermal properties of biological tissues in a wide temperature range and on the features of phase transitions [13]. It is recommended to use experimentally obtained temperaturedependent thermal properties adapted to simulate a specific cryoexposure [14,15]. Obtaining the properties of biological tissues, depending both on temperature and on the speed of freezing, is the next step in improving the quality of modeling. Also, since the properties of all types of biological tissues cannot be experimentally obtained, it is recommended to use predicting methods based on their composition [16].
When modeling a phase transition, it is necessary to control and take into account the features of the movement of its front. Thermal properties should not change too quickly with respect to the size of the mesh element. In the case of a sharp jump in thermal properties, the solver may make mistakes in calculating the phase transition region. To prevent errors, it is useful to set abrupt changes in the properties so that they account for at least 3 mesh elements in the longitudinal direction of the phase transition front, then usually the approximation error is within the calculation error.
Also, when modeling a phase transition, it is necessary to take into account the features of the cryoagent flow during boiling (liquid nitrogen, refrigerants). To calculate boiling, it is needed to create a fluid domain in which two liquids will be indicated -the vaporous and liquid phases of the substance. Setting up a domain, mass transfer between phases must be considered. In the interphase mass transfer options, it is recommended to set the temperature of saturated vapors (as a constant or function of pressure) and take into account the heat transfer between the liquid and the vapor bubbles. One of the most accurate equations for this type of heat transfer is the RPI Model, which needs to be activated on each surface of the domain on which boiling will occur. At the same time, for each wall it is possible to set boiling settings (frequency of bubble formation, their diameter, temperature of steam overheating, etc.), overriding similar settings common to the entire domain. Equally important is the correct setting of dynamic calculation. The calculation step should be extremely small, otherwise the solver will not be able to accurately calculate the mass components of the phases. In this case, the smaller the size of the mesh elements, the smaller the time calculation step is necessary for adequate modeling.

Boundary conditions
The most accurate way to model the sources of exposure is to consider their internal structure as a whole, or to apply the boundary conditions "fluid -wall", "coolant -biological tissue". However, this approach is not always possible (for example, when modeling argon cryoprobes). In such cases, the boundary conditions "biological tissue-cryoinstrument" are used, which can introduce a significant error in the simulation results [14]. In this case, it is recommended to use experimentally obtained boundary conditions in the form of a heat flux, depending on the type of cryoinstrument and its operating modes.
The external boundary of biological tissue may be set in two ways -by the adiabatic condition or by the temperature of the core of the human body. The first method is used when the dimensions of the biological tissue geometrical model are significantly larger than the area of cryoexposure. The second method is used when the temperature is known in a certain area -for example, the stable temperature of a certain area due to the influence of blood flow in large blood vessel.

Parameterization
The tasks of automating the search for optimal cryoexposure parameters have not been solved yet. Simulation is usually verification, while a physician based on recommendations and his own experience [2] determines cryoexposure planning. That is why an extremely important direction is multivariate calculations. The purpose of the calculations can be optimization of the location of exposure sources or exposure modes in order to increase the safety and effectiveness of exposure. At the same time, modern computer systems can not only reduce the calculation time, but also reduce the time taken to prepare the calculation for launch while carrying out a series of calculations in which individual values vary, or the calculation geometry change can be parameterized. It is possible to parameterize the geometric model and other solver settings that you plan to vary. This also applies to the dynamic settings of the calculation -the total time and the time step of the calculation (number of iterations), to the settings of the boundary conditions, thermal properties and other parameters. Software systems allows to set various properties of substances, parameters of boundary conditions in the form of various dependencies -on time, temperature, coordinate, or another variable. These conditions allow us to simulate the output of the source of exposure, the anisotropy of the thermal properties, and transient physical processes.

Multiprobe minimally invasive cryoablation of a malignant prostate tumor
The operating modes and arrangement of cryoprobes were corresponded to the previous cryooperation. A typical algorithm consists of two freezing cycles of 10 minutes each, between which there was a break of 5 minutes for defrost. As a result of the simulation, the complete freezing of the target tumor region to the temperature of cryonecrosis was demonstrated (Figure 2). The described recommendations were applied to create a geometric model, taking into account thermal properties (the experimentally obtained properties of a healthy and malignant prostate in a wide temperature range were used [15]. To clarify the boundary conditions, the average capacities of cryoprobes with one (IceSeed) and two (IceRod) capillary tubes were estimated depending on the operating mode of the cryoprobe (Table 1). This presented data should be used as a first approximation, while accurate modeling was carried out with boundary conditions depending on temperature [17][18][19]. It should be noted that in the practice of cryoablation, the first freezing cycle is usually carried out at maximum cooling power. The second and subsequent cycles may be carried out at a different modes. In this case, the first cycle has a significant load in the form of a phase transition, the cryoprobe entering the mode and cooling of the biological tissue from the initial temperature to the temperature at which the phase transition begins. This explains that the 80% mode provides more cooling power than the 100% mode. Figure 3 shows the parameterization diagram in the ANSYS DesignModeler of a warming catheter surrounded by cryoprobes. At the stage of a flat sketch, the center of the catheter was set by two coordinates, a separate parameter was responsible for the diameter of the catheter. Similarly, all cryoprobes were parameterized. Each side of the surface participating in heat transfer is assigned a separate name, which is further assigned to the corresponding boundary conditions. This made it possible to save all the settings of the mesh generator. Changing the value of any geometry parameters, the finite element mesh is easily rebuilt without re-setting.

Partial body cryotherapy
The most difficult task was the physical model of the patient. Testing the program showed that the use of the geometric model closest to the real shape of the human body complicates the calculation and causes algorithm errors, which increases the accuracy of the result. Therefore, a simplification of the shape of the cooling object has been applied. Also, this task did not include phase transition, therefore, the basic properties of biological tissues are taken in the form of constants. However, blood perfusion was taken into account as a temperaturedependent source term [20]. The boundary condition "fluid -biological tissue" simplified the statement of this problem. Using parameterization, a functional dependence of the velocity and temperature of the incoming air was set.