Models of cooling and freezing of living biological tissues with a flat ruler applicator

Abstract

In this work, we construct new mathematical models of cooling and freezing of living biological tissue with a flat, long ruler applicator located on its surface. The models are two-dimensional boundary-value problems (including Stefan-type problems) and have applications in cryosurgery. The method of numerical study of these problems is based on smoothing discontinuous functions and applying locally one-dimensional difference schemes to “smoothed” problems without explicitly identifying the boundaries of the influence of cold and the boundaries of the phase transition.

Full Text

1. Введение. Лечение холодом в медицине известно достаточно давно. Использование низких температур для охлаждения и замораживания биологической ткани применяется, в частности, при хранении и консервации биоматериалов, а также при проведении криохирургических операций. В криохирургии используются криогенные температуры в $-{150}^{\circ }$ C и более низкие. В отличие от консервации и хранения, цель криохирургии состоит в гибели клеток в локальном, четко ограниченном объеме биоткани, занимаемом злокачественной опухолью. Гибель клеток достигается в результате разрыва мембран, образующегося при криогенном охлаждении кристаллами льда внеклеточной и внутриклеточной воды, а также осмотического разбухания при оттаивании биоткани. Замерзание жидкости в биоткани всегда начинается во внеклеточном пространстве ( $\mathrm{0<mo>:</mo>}-3^{\circ }$ C), а при дальнейшем понижении температуры распространяется и на внутриклеточную жидкость ( $-\mathrm{20<mo>:</mo>}-{50}^{\circ }$ C). Основная задача криохирургии состоит в контроле за гибелью всех злокачественных клеток в заданном объеме биологической ткани. Очевидно, что успешное развитие криохирургического метода и прогнозирование его результатов во многом зависит от достоверного описания теплового процесса замораживания биоткани, сопровождающегося фазовым переходом вода-лед. Отличительной особенностью возникающих здесь задач типа Стефана, является пространственная локализация тепла и существование стационарных решений.

Существует большое количество научных публикаций (см., например, обзор литературы в [8]), в которых рассматриваются математические модели охлаждения и замораживания живой биологической ткани. Основой их, как правило, является уравнение Пеннеса (или его незначительные модификации; см. [11]), из которого виден линейный характер зависимости источников тепла биоткани от искомого температурного поля. Такой характер зависимости не позволяет описать реально наблюдаемую пространственную локализацию тепла. Кроме того, модель Пеннеса не учитывает того факта, что замерзание межклеточной жидкости происходит гораздо раньше, чем замерзание внутриклеточной жидкости, и соответствующее этим двум процессам тепло выделяется в разные моменты времени. Известно, что при охлаждении биологической ткани в криомедицине используются криоинструменты с различными формами охлаждающей поверхности. Криоинструменты могут располагаться как на поверхности биоткани (аппликаторы), так и полностью внедряться в нее (зонды). С понижением температуры охлаждающей поверхности в ткани возникает нестационарное температурное поле, зависящее в общем случае от трех пространственных координат и времени. Интерес представляют как распределение температурного поля в ткани, так и размеры зон криопоражения, замораживания и влияния холода, а также время выхода на стационар.

В предлагаемой работе построены новые математические модели охлаждения и замораживания живой биологической ткани плоским, достаточно протяженным линейчатым аппликатором, располагаемым на ее поверхности. Модели учитывает указанные выше особенности, представляют собой двумерные краевые задачи с нелинейными источниками тепла специального вида и имеют приложение в криохирургии. Предложен метод численного исследования поставленных задач, основанный на сглаживании разрывных функций и применении к <<сглаженным>> задачам локально-одномерных разностных схем без явного выделения границы влияния холода и границ фазового перехода. Приведены некоторые численные расчеты на ЭВМ.

2. Постановка задачи. Как было сказано выше, биоткань может охлаждаться криохирургическими инструментами различной формы (плоской, цилиндрической и т. д.). Далее рассматриваются математические модели охлаждения и замораживания живой биоткани плоским протяженным линейчатым аппликатором. Пренебрегая краевыми эффектами, возникающими вблизи концов, приходим к постановкам двумерных начально-краевых задач типа Стефана в прямоугольной области. Соответствующие уравнения и дополнительные условия будут выглядеть следующим образом:

$\frac{\partial }{\partial x}\left[\lambda \mathrm{<mo\; stretchy="false">(</mo>}u\mathrm{<mo\; stretchy="false">)</mo>}\frac{\partial u}{\partial x}\right]+\frac{\partial }{\partial y}\left[\lambda \mathrm{<mo\; stretchy="false">(</mo>}u\mathrm{<mo\; stretchy="false">)</mo>}\frac{\partial u}{\partial y}\right]+c\mathrm{<mo\; stretchy="false">(</mo>}u\mathrm{<mo\; stretchy="false">)</mo>}\rho \mathrm{<mo\; stretchy="false">(</mo>}u\mathrm{<mo\; stretchy="false">)</mo>}\frac{\partial u}{\partial t}\mathrm{<mo>=</mo>}$

$\mathrm{<mo>=</mo>}-w\mathrm{<mo\; stretchy="false">(</mo>}u\mathrm{<mo\; stretchy="false">)</mo>}+P_1\frac{\partial u}{\partial t}\delta \mathrm{<mo\; stretchy="false">(</mo>}u-u_1\mathrm{<mo\; stretchy="false">)</mo>}+P_2\frac{\partial u}{\partial t}\delta \mathrm{<mo\; stretchy="false">(</mo>}u-u_2\mathrm{<mo\; stretchy="false">)</mo>,}x\in \mathrm{<mo\; stretchy="false">(</mo>}-a\mathrm{,}a\mathrm{<mo\; stretchy="false">)</mo>,}y\in \mathrm{<mo\; stretchy="false">(</mo>0,}b\mathrm{<mo\; stretchy="false">)</mo>,}t\mathrm{<mo>></mo>0,}$ (1)

$u\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}y\mathrm{,}z\mathrm{,0<mo\; stretchy="false">)</mo><mo>=</mo>}\overline{u}\mathrm{<mo>=</mo>}const\mathrm{,}$

$\lambda \mathrm{<mo\; stretchy="false">(</mo>}u\mathrm{<mo\; stretchy="false">)</mo>}\frac{\partial u}{\partial z}-\alpha u\mathrm{<mo>=</mo>}-\alpha u_3\mathrm{,}y\mathrm{<mo>=</mo>0,}x\in \mathrm{<mo\; stretchy="false">[</mo>}-r_0\mathrm{,}{r}_0\mathrm{<mo\; stretchy="false">]</mo>,}t\mathrm{<mo>></mo>0,}$

$\lambda \mathrm{<mo\; stretchy="false">(</mo>}u\mathrm{<mo\; stretchy="false">)</mo>}\frac{\partial u}{\partial y}-\gamma u\mathrm{<mo>=</mo>}\gamma u_4\mathrm{,}y\mathrm{<mo>=</mo>0,}x\notin \mathrm{<mo\; stretchy="false">[</mo>}-r_0\mathrm{,}{r}_0\mathrm{<mo\; stretchy="false">]</mo>,}t\mathrm{<mo>></mo>0,}$

$u\mathrm{<mo\; stretchy="false">(</mo>}-a\mathrm{,}y\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}\overline{u}\mathrm{,}y\in \mathrm{<mo\; stretchy="false">[</mo>0,}b\mathrm{<mo\; stretchy="false">]</mo>,}t\mathrm{<mo>></mo>0,}$

$u\mathrm{<mo\; stretchy="false">(</mo>}a\mathrm{,}y\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}\overline{u}\mathrm{,}y\in \mathrm{<mo\; stretchy="false">[</mo>0,}b\mathrm{<mo\; stretchy="false">]</mo>,}t\mathrm{<mo>></mo>0,}$

$u\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}b\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}\overline{u}\mathrm{,}x\in \mathrm{<mo\; stretchy="false">[</mo>}-a\mathrm{,}a\mathrm{<mo\; stretchy="false">]</mo>,}t\mathrm{<mo>></mo>0,}$

$u\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}{y}_1\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}y\mathrm{<mo\; stretchy="false">)</mo>,}t\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{u}_1\mathrm{,}u\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}{y}_2\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}y\mathrm{<mo\; stretchy="false">)</mo>,}t\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{u}_2$ (2)

Здесь $\lambda \mathrm{<mo\; stretchy="false">(</mo>}u\mathrm{<mo\; stretchy="false">)</mo>}$, $c\mathrm{<mo\; stretchy="false">(</mo>}u\mathrm{<mo\; stretchy="false">)</mo>}$, $\rho \mathrm{<mo\; stretchy="false">(</mo>}u\mathrm{<mo\; stretchy="false">)</mo>}$ $—$ коэффициенты теплопроводности, теплоемкости и плотности биоткани соответственно; $2r_0$ $—$ ширина плоского криоинструмента, $u_3\mathrm{<mo>=</mo>}{u}_3\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}$ $—$ температура его охлаждающей поверхности; $u_4$ $—$ температура окружающей среды; $\alpha$ $—$ коэффициент теплообмена биоткани с аппликатором на участке $\mathrm{<mo\; stretchy="false">[</mo>}-r_0\mathrm{,}{r}_0\mathrm{<mo\; stretchy="false">]</mo>}$; $\gamma$ $—$ коэффициент теплообмена с окружающей средой; $\overline{u}$ $—$ температура биоткани, до которой еще не дошел холод; $u_1$ $—$ температура замораживания биоткани, $u_2$ $—$ температура криопоражения, $P_1$, $P_2$ $—$ известные постоянные (их физический смысл описан ниже); $\delta \mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{<mo\; stretchy="false">)</mo>}$ $—$ дельта-функция Дирака.

Определению подлежит функция температуры $u\mathrm{<mo>=</mo>}u\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}y\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}$, пара изотермических поверхностей $y_1\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}$, $y_2\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}$ на которых температура биоткани равна, соответственно, $u_1$, $u_2$; а также постоянные $a$, $b$, характеризующиеся тем, что вне рассматриваемого прямоугольника температура биоткани постоянна и равна $\overline{u}$. Их приближенные значения можно найти из решения одномерной стационарной задачи (см. [1]).

Будем предполагать, что коэффициенты $c\mathrm{<mo\; stretchy="false">(</mo>}u\mathrm{<mo\; stretchy="false">)</mo>}$, $\rho \mathrm{<mo\; stretchy="false">(</mo>}u\mathrm{<mo\; stretchy="false">)</mo>}$, $\lambda \mathrm{<mo\; stretchy="false">(</mo>}u\mathrm{<mo\; stretchy="false">)</mo>}$ могут иметь разрывы типа скачка при $u\mathrm{<mo>=</mo>}{u}_1$ и $u\mathrm{<mo>=</mo>}{u}_2$ так что

$c\mathrm{<mo\; stretchy="false">(</mo>}u\mathrm{<mo\; stretchy="false">)</mo>}\ge c_{\min }\mathrm{<mo>></mo>0,}\rho \mathrm{<mo\; stretchy="false">(</mo>}u\mathrm{<mo\; stretchy="false">)</mo>}\ge {\rho }_{\min }\mathrm{<mo>></mo>0,}\lambda \mathrm{<mo\; stretchy="false">(</mo>}u\mathrm{<mo\; stretchy="false">)</mo>}\ge {\lambda }_{\min }\mathrm{<mo>></mo>0.}$

В случае, когда происходит <<охлаждение>> биоткани ( $u_3\mathrm{<mo>></mo>0}$ ), уравнение (1) будет иметь вид

$\frac{\partial }{\partial x}\left[\lambda \mathrm{<mo\; stretchy="false">(</mo>}u\mathrm{<mo\; stretchy="false">)</mo>}\frac{\partial u}{\partial x}\right]+\frac{\partial }{\partial y}\left[\lambda \mathrm{<mo\; stretchy="false">(</mo>}u\mathrm{<mo\; stretchy="false">)</mo>}\frac{\partial u}{\partial y}\right]+c\mathrm{<mo\; stretchy="false">(</mo>}u\mathrm{<mo\; stretchy="false">)</mo>}\rho \mathrm{<mo\; stretchy="false">(</mo>}u\mathrm{<mo\; stretchy="false">)</mo>}\frac{\partial u}{\partial t}\mathrm{<mo>=</mo>}-w\mathrm{<mo\; stretchy="false">(</mo>}u\mathrm{<mo\; stretchy="false">)</mo>,}x\in \mathrm{<mo\; stretchy="false">(</mo>}-a\mathrm{,}a\mathrm{<mo\; stretchy="false">)</mo>,}y\in \mathrm{<mo\; stretchy="false">(</mo>0,}b\mathrm{<mo\; stretchy="false">)</mo>,}t\mathrm{<mo>></mo>0,}$ (3)

т.е. слагаемых с <<дельта-функциями>> не будет. В этом случае в постановке задачи также будут отсутствовать условия (2) на изотермических поверхностях. В остальном модель останется прежней.

3. Метод решения. Далее рассматривается задача для уравнения (1). Для определения приближенных значений температуры $u\mathrm{<mo>=</mo>}u\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}y\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}$ и положения изотермических поверхностей, отвечающих значениям $u\mathrm{<mo>=</mo>}u\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}y\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{u}_1$, $u\mathrm{<mo>=</mo>}u\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}y\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{u}_2$ предлагается следующий алгоритм.

Введем функцию

$H\mathrm{<mo\; stretchy="false">(</mo>}u\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{\displaystyle \underset{0}{\overset{u}{\int }}}c\mathrm{<mo\; stretchy="false">(</mo>}\xi \mathrm{<mo\; stretchy="false">)</mo>}\rho \mathrm{<mo\; stretchy="false">(</mo>}\xi \mathrm{<mo\; stretchy="false">)</mo>}d\xi +P_1\eta \mathrm{<mo\; stretchy="false">(</mo>}u-u_1\mathrm{<mo\; stretchy="false">)</mo>}+P_2\eta \mathrm{<mo\; stretchy="false">(</mo>}u-u_2\mathrm{<mo\; stretchy="false">)</mo>,}$

которая при $u\mathrm{<mo>=</mo>}{u}_1$ и $u\mathrm{<mo>=</mo>}{u}_2$ имеет скачки соответственно $P_1$ и $P_2$; $\eta$ $—$ функция Хевисайда. Уравнение (1) теперь можно переписать в виде

$\frac{\partial H\mathrm{<mo\; stretchy="false">(</mo>}u\mathrm{<mo\; stretchy="false">)</mo>}}{\partial x}\mathrm{<mo>=</mo>}\frac{\partial }{\partial x}\left[\lambda \mathrm{<mo\; stretchy="false">(</mo>}u\mathrm{<mo\; stretchy="false">)</mo>}\frac{\partial u}{\partial x}\right]+\frac{\partial }{\partial y}\left[\lambda \mathrm{<mo\; stretchy="false">(</mo>}u\mathrm{<mo\; stretchy="false">)</mo>}\frac{\partial u}{\partial y}\right]+w\mathrm{<mo\; stretchy="false">(</mo>}u\mathrm{<mo\; stretchy="false">)</mo>.}$ (4)

Далее в (4) проводим сглаживание $H\mathrm{<mo\; stretchy="false">(</mo>}u\mathrm{<mo\; stretchy="false">)</mo>}$, $\lambda \mathrm{<mo\; stretchy="false">(</mo>}u\mathrm{<mo\; stretchy="false">)</mo>}$ по на $u$ интервалах

$u_1-{\overline{\overline{\delta }}}_k\le u\le u_1+{\overline{\overline{\delta }}}_k\mathrm{,}{u}_2-{\overline{\delta }}_k\le u\le u_2+{\overline{\delta }}_k\mathrm{,}$

где ${\overline{\delta }}_k\to 0$ и ${\overline{\overline{\delta }}}_k\to 0$ при $k\to \infty$ (см. [4]). В результате получаем последовательности ограниченных гладких функций $H^k\mathrm{<mo\; stretchy="false">(</mo>}u\mathrm{<mo\; stretchy="false">)</mo>}$, ${\lambda }^k\mathrm{<mo\; stretchy="false">(</mo>}u\mathrm{<mo\; stretchy="false">)</mo>}$, сходящихся при $k\to +\infty$ и $u\ne u_1$, $u\ne u_2$ соответственно к $H\mathrm{<mo\; stretchy="false">(</mo>}u\mathrm{<mo\; stretchy="false">)</mo>}$, $\lambda \mathrm{<mo\; stretchy="false">(</mo>}u\mathrm{<mo\; stretchy="false">)</mo>}$, причем $B^k\mathrm{<mo\; stretchy="false">(</mo>}u\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{H}_u^k\mathrm{<mo\; stretchy="false">(</mo>}u\mathrm{<mo\; stretchy="false">)</mo>}\ge D\mathrm{<mo>=</mo>}const\mathrm{<mo>></mo>0}$.

Задачу Стефана (1), (2) заменяем аппроксимирующей <<сглаженной>> задачей:

$\frac{\partial }{\partial x}\left[{\lambda }^k\left(u^k\right)\frac{\partial u^k}{\partial x}\right]+\frac{\partial }{\partial y}\left[{\lambda }^k\left(u^k\right)\frac{\partial u^k}{\partial y}\right]-B^k\left(u^k\right)\frac{\partial u}{\partial x}\mathrm{<mo>=</mo>}-w\mathrm{<mo\; stretchy="false">(</mo>}u\mathrm{<mo\; stretchy="false">)</mo>,}x\in \mathrm{<mo\; stretchy="false">(</mo>}-a\mathrm{,}a\mathrm{<mo\; stretchy="false">)</mo>,}y\in \mathrm{<mo\; stretchy="false">(</mo>0,}b\mathrm{<mo\; stretchy="false">)</mo>,}t\mathrm{<mo>></mo>0<mo>;</mo>}$ (5)

$u^k\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}y\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}\overline{u}\mathrm{<mo>=</mo>}const\mathrm{,}x\in \mathrm{<mo\; stretchy="false">[</mo>}-a\mathrm{,}a\mathrm{<mo\; stretchy="false">]</mo>,}y\in \mathrm{<mo\; stretchy="false">[</mo>0,}b\mathrm{<mo\; stretchy="false">]</mo>,}$

${\lambda }^k\left(u^k\right)\frac{\partial u^k}{\partial y}-\alpha u^k\mathrm{<mo>=</mo>}\alpha u_3\mathrm{,}y\mathrm{<mo>=</mo>0,}x\in \mathrm{<mo\; stretchy="false">[</mo>}-r_0\mathrm{,}{r}_0\mathrm{<mo\; stretchy="false">]</mo>,}y\mathrm{<mo>></mo>0,}$

${\lambda }^k\left(u^k\right)\frac{\partial u^k}{\partial y}-\gamma u^k\mathrm{<mo>=</mo>}\gamma u_4\mathrm{,}y\mathrm{<mo>=</mo>0,}x\in \mathrm{<mo\; stretchy="false">[</mo>}-\alpha \mathrm{,}{r}_0\mathrm{<mo\; stretchy="false">]</mo>}\cup \mathrm{<mo\; stretchy="false">[</mo>}{r}_0\mathrm{,}a\mathrm{<mo\; stretchy="false">]</mo>,}t\mathrm{<mo>></mo>0,}$

$u^k\mathrm{<mo\; stretchy="false">(</mo>}-a\mathrm{,}y\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}\overline{u}\mathrm{,}y\in \mathrm{<mo\; stretchy="false">[</mo>0,}b\mathrm{<mo\; stretchy="false">]</mo>,}t\mathrm{<mo>></mo>0,}$

$u^k\mathrm{<mo\; stretchy="false">(</mo>}a\mathrm{,}y\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}\overline{u}\mathrm{,}y\in \mathrm{<mo\; stretchy="false">[</mo>0,}b\mathrm{<mo\; stretchy="false">]</mo>,}t\mathrm{<mo>></mo>0,}$

$u^k\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}b\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}\overline{u}\mathrm{,}x\in \mathrm{<mo\; stretchy="false">[</mo>}-a\mathrm{,}a\mathrm{<mo\; stretchy="false">]</mo>,}t\mathrm{<mo>></mo>0.}$

Отметим, что условия (2) на изотермических поверхностях теперь содержатся в уравнении (5). Предполагая, что решение задачи (5) существует и является достаточно гладким, применим к этой задаче локально-одномерный метод, описанный в [9]. Ради простоты приведем схему локально-одномерного метода, когда задача (5) решается в прямоугольной области, хотя его можно использовать и для произвольной области $G$, используя лишь одно предположение: пересечение с областью $G$ прямой $L$, проведенной через любую внутреннюю точку $G$ параллельно какой-либо из осей координат, состоит из конечного числа интервалов. Введем произвольную неравномерную сетку:

$W_h\mathrm{<mo>=</mo><mo>\{</mo><mo\; stretchy="false">(</mo>}{x}_i\mathrm{,}{y}_j\mathrm{<mo\; stretchy="false">)</mo>,}i\mathrm{<mo>=</mo>1,2,}\dots \mathrm{,}{N}_1\mathrm{,}j\mathrm{<mo>=</mo>1,2,}\dots \mathrm{,}{N}_2\mathrm{,}{x}_0\mathrm{<mo>=</mo>}-a\mathrm{,}{y}_0\mathrm{<mo>=</mo>0,}{x}_{N_1}\mathrm{<mo>=</mo>}a\mathrm{,}{y}_{N_2}\mathrm{<mo>=</mo>}b\mathrm{<mo>\}</mo>}$

и обозначим через ${\nu }_{i\mathrm{,}j}^n$ значение искомой функции в точке $\mathrm{<mo\; stretchy="false">(</mo>}{x}_i\mathrm{,}{y}_j\mathrm{<mo\; stretchy="false">)</mo>}$ в момент времени $t_n\mathrm{<mo>=</mo>}n\tau$. Пусть ${\nu }_i^{n+\mathrm{1<mo>/</mo>2}}$ $—$ значение искомой функции на полуслое. Будем считать, что отрезок $\mathrm{<mo\; stretchy="false">[</mo>}-r_0\mathrm{,}{r}_0\mathrm{<mo\; stretchy="false">]</mo>}$ разбит на $m$ узлов. Локально-одномерная схема для задачи (5) будет иметь следующий вид:

$\frac{1}{2}B\mathrm{<mo\; stretchy="false">(</mo>}{\nu }_{ij}^{n+\mathrm{1<mo>/</mo>2}}\mathrm{<mo\; stretchy="false">)</mo>}\frac{{\nu }_{ij}^{n+\mathrm{1<mo>/</mo>2}}-{\nu }_{ij}^n}{\tau }\mathrm{<mo>=</mo>}{\Lambda }_1{\nu }_{ij}^{n+\mathrm{1<mo>/</mo>2}}+\frac{1}{2}w\left({\nu }_{ij}^{n+\mathrm{1<mo>/</mo>2}}\right)\mathrm{,}\left(x_i\mathrm{,}{y}_j\right)\in W_h\mathrm{,}$

${\nu }_{ij}^{n+\mathrm{1<mo>/</mo>2}}\mathrm{<mo>=</mo>}\overline{u}$, при $i\mathrm{<mo>=</mo>0}$, $i\mathrm{<mo>=</mo>}{N}_1$, $\mathrm{0<mo><</mo>}j\mathrm{<mo><</mo>}{N}_2$;

${\Lambda }_1{\nu }_{ij}^{n+\mathrm{1<mo>/</mo>2}}\mathrm{<mo>=</mo>}\frac{1}{h}\left[\lambda \left({\nu }_{i+\mathrm{1<mo>/</mo>2}}^{n+\mathrm{1<mo>/</mo>2}}\right)\frac{{\nu }_{i+\mathrm{1<mo>/</mo>2}}^{n+\mathrm{1<mo>/</mo>2}}-{\nu }_{ij}^{n+\mathrm{1<mo>/</mo>2}}}{h}-\lambda \left({\nu }_{i-\mathrm{1<mo>/</mo>2,}j}^{n+\mathrm{1<mo>/</mo>2}}\right)\frac{{\nu }_{ij}^{n+\mathrm{1<mo>/</mo>2}}-{\nu }_{i-\mathrm{1,}j}^{n+\mathrm{1<mo>/</mo>2}}}{h}\right]\mathrm{,}$

$\frac{1}{2}B\left({\nu }_{ij}^{n+1}\right)\frac{{\nu }_{ij}^{n+1}-{\nu }_{ij}^{n+\mathrm{1<mo>/</mo>2}}}{\tau }\mathrm{<mo>=</mo>}{\Lambda }_2{\nu }_{ij}^{n+1}+\frac{1}{2}w\left({\nu }_{ij}^{n+1}\right)\mathrm{,}\left(x_i\mathrm{,}{y}_j\right)\in W_h\mathrm{,}$

${\nu }_{i0}^{n+1}\mathrm{<mo>=</mo>}\frac{\lambda \left({\nu }_{i0}^{n+1}\right)}{\gamma h+\lambda \left({\nu }_{i0}^{n+1}\right)}{\nu }_{i1}^{n+1}+\frac{\gamma h{u}_4}{\gamma h+\lambda \left({\nu }_{i0}^{n+1}\right)}\mathrm{,}\mathrm{0<mo><</mo>}i\mathrm{<mo><</mo>}\frac{N_1-m}{2}\mathrm{,}\frac{N_1+m}{2}\mathrm{<mo><</mo>}i\mathrm{<mo><</mo>}{N}_1\mathrm{,}$

${\nu }_{i0}^{n+1}\mathrm{<mo>=</mo>}\frac{\lambda \left({\nu }_{i0}^{n+1}\right)}{\alpha h+\lambda \left({\nu }_{i0}^{n+1}\right)}{\nu }_{i1}^{n+1}+\frac{\alpha h{u}_3}{\alpha h+\lambda \left({\nu }_{i0}^{n+1}\right)}\mathrm{,}\frac{N_1-m}{2}\le i\le \frac{N_1+m}{2}\mathrm{,}$

${\nu }_{i0}^{n+1}\mathrm{<mo>=</mo>}\overline{u}\mathrm{,}\mathrm{0<mo><</mo>}i\mathrm{<mo><</mo>}{N}_1\mathrm{<mo>;</mo>}$

${\Lambda }_2{\nu }_{ij}^{n+1}\mathrm{<mo>=</mo>}\frac{1}{h}\left[\lambda \left({\nu }_{i\mathrm{,}j+\mathrm{1<mo>/</mo>2}}^{n+1}\right)\frac{{\nu }_{i\mathrm{,}j+1}^{n+1}-{\nu }_{ij}^{n+1}}{h}-\lambda \left({\nu }_{i\mathrm{,}j-\mathrm{1<mo>/</mo>2}}^{n+1}\right)\frac{{\nu }_{ij}-{\nu }_{i\mathrm{,}j-1}^{n+1}}{h}\right]$

 ${\nu }_{ij}^0\mathrm{<mo>=</mo>}\overline{u}$ при $t\mathrm{<mo>=</mo>0}$, $0\le i\le N_1$, $0\le j\le N_2$. Здесь

${\nu }_{i+\mathrm{1<mo>/</mo>2,}j}\mathrm{<mo>=</mo>}\frac{1}{2}\left({\nu }_{ij}+{\nu }_{i+\mathrm{1,}j}\right)\mathrm{,}{\nu }_{i\mathrm{,}j+\mathrm{1<mo>/</mo>2}}\mathrm{<mo>=</mo>}\frac{1}{2}\left({\nu }_{ij}+{\nu }_{i\mathrm{,}j+1}\right)\mathrm{.}$

Из написанных выше уравнений видно, что для нахождения значения $\nu$ (при $t\mathrm{<mo>=</mo>}{t}_i+1$ ) на $\mathrm{<mo\; stretchy="false">(</mo>}n+\mathrm{1<mo\; stretchy="false">)</mo>}$ -м временном слое по известному значению на $j$ -м временном слое нужно последовательно решать две серии одномерных задач, соответственно, по координатам $x$ и $y$. Каждая такая задача представляет собой нелинейную алгебраическую систему с трехдиагональной матрицей, и для ее решения лучше пользоваться методом прогонки совместно с каким-либо итерационным методом. В настоящей работе использовался итерационный метод Ньютона. При определении ${\nu }^{n+1}$ методом итераций коэффициенты $c$, $\rho$, $\lambda$ можно брать на предыдущей итерации (см. [10]).

Отличительной особенностью примененного для счета метода является возможность использования большого шага по времени, что существенно сокращает время счета при решении многомерных задач. Анализ результатов показывает, что область определения сглаживающей функции лучше выбирать так, чтобы охватывалось 2-3 счетные точки. В то же время можно сделать вывод, что вид самой функции (т.е. порядок сглаживания) слабо влияет на результаты счета. В настоящей работе велся расчет для задачи с кусочно постоянными коэффициентами. Однако, этого требовать вовсе не обязательно, и никаких усложнений методики при этом не возникает.

Итерации, которые применяются для решения указанных выше задач, сходятся быстро (5-7 итераций).

4. Источники тепла. Как и в обычных задачах Стефана, замораживание биоткани сопровождается выделением тепла при кристаллизации сначала внеклеточной, а затем и внутриклеточной воды. Это моделируется разрывностью удельной внутренней тепловой энергии в точках $u\mathrm{<mo>=</mo>}{u}_1$, $u\mathrm{<mo>=</mo>}{u}_2$. Величины скачков пропорциональны плотностям внеклеточной ${\rho }_1$ и внутриклеточной ${\rho }_2$ воды: $P_1\mathrm{<mo>=</mo>}{\rho }_1\Lambda$, $P_2\mathrm{<mo>=</mo>}{\rho }_2\Lambda$, где $\Lambda$ $—$ удельная теплота кристаллизации воды. Сами источники сосредоточены на подлежащих определению изотермических поверхностях: $u\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}y\mathrm{,}z\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{u}_1$, $u\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}y\mathrm{,}z\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{u}_2$.

Далее, в связи с тем, что биоткань пронизана разветвленной сетью капилляров, снабжающей кровью охлажденную не замороженную и замороженную, но не криопораженную области биоткани, в ней возникают внутренние источники тепла. По теплофизическому смыслу функциональная зависимость этих источников тепла от температуры должна быть ограниченной, непрерывной и монотонно убывающей в интервале положительных температур, а в интервале отрицательных температур $—$ ограниченной и монотонно возрастающей. Для существования пространственной локализации необходимым является условие $w^{\prime }\mathrm{<mo\; stretchy="false">(</mo>}\overline{u}\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}-\infty$. Оно имеет следующий физический смысл: при сколь угодно малом возмущении начальной температуры биоткани в ней возникают сколь угодно малые источники, скорость нарастания которых неограниченна (см. [2]). Кроме того, необходимым условием является сходимость интеграла

${\displaystyle \underset{\underset{¯}{u}}{\overset{\overline{u}}{\int }}}\frac{\lambda \mathrm{<mo\; stretchy="false">(</mo>}u\mathrm{<mo\; stretchy="false">)</mo>}du}{{\left({\displaystyle \underset{u}{\overset{\overline{u}}{\int }}}w\mathrm{<mo\; stretchy="false">(</mo>}u\mathrm{<mo\; stretchy="false">)</mo>}\lambda \mathrm{<mo\; stretchy="false">(</mo>}u\mathrm{<mo\; stretchy="false">)</mo>}du\right)}^{\mathrm{1<mo>/</mo>2}}}\mathrm{<mo><</mo>}\infty \mathrm{.}$

Поскольку $\lambda \mathrm{<mo>=</mo>}\lambda \mathrm{<mo\; stretchy="false">(</mo>}u\mathrm{<mo\; stretchy="false">)</mo>}$ $—$ строго положительная функция, то достаточно, чтобы

${\displaystyle \underset{\underset{¯}{u}}{\overset{\overline{u}}{\int }}}\frac{du}{{\left({\displaystyle \underset{u}{\overset{\overline{u}}{\int }}}w\mathrm{<mo\; stretchy="false">(</mo>}u\mathrm{<mo\; stretchy="false">)</mo>}du\right)}^{\mathrm{1<mo>/</mo>2}}}\mathrm{<mo><</mo>}\infty \mathrm{,}$

где $\underset{¯}{u}$ $—$ минимальная температура биоткани.

При проведении численных расчетов на ЭВМ (см. [5–7]) хорошо зарекомендовали себя следующие функции источников тепла:

$w_1\mathrm{<mo\; stretchy="false">(</mo>}u\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}\left\{\begin{array}{ll}\hfill & w_0{\mathrm{<mo\; stretchy="false">(</mo>}\overline{u}-u\mathrm{<mo\; stretchy="false">)</mo>}}^{\beta }\mathrm{,}{u}_1\le u\mathrm{<mo><</mo>}\overline{u}\mathrm{,}\hfill \\ \hfill & 0\text{в противном случае;}\hfill \end{array}\right.w_2\mathrm{<mo\; stretchy="false">(</mo>}u\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}\left\{\begin{array}{ll}\hfill & w_{0}ln\mathrm{<mo\; stretchy="false">(</mo>}\overline{u}-u\mathrm{<mo\; stretchy="false">)</mo>,}{u}_1\le u\mathrm{<mo><</mo>}\overline{u}\mathrm{,}\hfill \\ \hfill & 0\text{в противном случае;}\hfill \end{array}\right.$

$w_3\mathrm{<mo\; stretchy="false">(</mo>}u\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}\left\{\begin{array}{llll}\hfill & w_0{\mathrm{<mo\; stretchy="false">(</mo>}\overline{u}-u\mathrm{<mo\; stretchy="false">)</mo>}}^{\beta }\mathrm{,}\hfill & \hfill & u_1\le u\mathrm{<mo><</mo>}\overline{u}\mathrm{,}\hfill \\ -\hfill & \frac{w_0{\overline{u}}^{\beta }}{u_A}\mathrm{<mo\; stretchy="false">(</mo>}u-u_A\mathrm{<mo\; stretchy="false">)</mo>,}\hfill & \hfill & u_A\le u\mathrm{<mo><</mo>}{u}_1\mathrm{,}\hfill \\ \hfill & 0\hfill & \hfill & u\mathrm{<mo><</mo>}{u}_A\mathrm{.}\hfill \end{array}\right.$

Содержащиеся в указанных зависимостях параметры $w_0$ и $\mathrm{0<mo><</mo>}\beta \mathrm{<mo><</mo>1}$ должны определяться экспериментально

5. Расчетные данные и результаты численных экспериментов. Расчеты проводились при следующих значениях теплофизических характеристик биоткани (см. [3, 5, 6]):

${\lambda }_1\mathrm{<mo>=</mo>0,56}\frac{\text{Вт}}{\text{м}{\cdot }^{\circ }\text{C}}\mathrm{,}{\lambda }_2\mathrm{<mo>=</mo>}{\lambda }_3\mathrm{<mo>=</mo>2,22}\frac{\text{Вт}}{\text{м}{\cdot }^{\circ }\text{C}}\mathrm{,}{w}_0\mathrm{<mo>=</mo>48,5}\cdot {10}^3\frac{\text{Вт}}{\text{м}{\cdot }^{\circ }\text{C}}\mathrm{,}$

$c_1\cdot {\rho }_1\mathrm{<mo>=</mo>3,6}\cdot {10}^6\frac{\text{Вт}\cdot \text{с}}{{\text{м}}^3{\cdot }^{\circ }\text{C}}\mathrm{,}{c}_2\cdot {\rho }_2\mathrm{<mo>=</mo>2,01}\cdot {10}^6\frac{\text{Вт}\cdot \text{с}}{{\text{м}}^3{\cdot }^{\circ }\text{C}}\mathrm{,}{c}_3\cdot {\rho }_3\mathrm{<mo>=</mo>1,08}\cdot {10}^6\frac{\text{Вт}\cdot \text{с}}{{\text{м}}^3{\cdot }^{\circ }\text{C}}\mathrm{,}$

$P_2\mathrm{<mo>=</mo>300}\cdot {10}^6\frac{\text{Вт}\cdot \text{с}}{{\text{м}}^3}\mathrm{,}{P}_1\mathrm{<mo>=</mo>0,3}\cdot P_2\mathrm{,}\alpha \mathrm{<mo>=</mo>100}÷2\cdot {10}^5\frac{\text{Вт}}{{\text{м}}^2{\cdot }^{\circ }\text{C}}\mathrm{,}$

$u_3\mathrm{<mo>=</mo>}-\mathrm{90,}\overline{u}\mathrm{<mo>=</mo>36,7.}$

Индексы 1, 2, 3 относятся соответственно к охлажденной, замороженной и криопораженной области биоткани. Некоторые результаты проведенных расчетов приведены ниже в таблице 1. Расчеты проведены с использованием функции $w_1\mathrm{<mo\; stretchy="false">(</mo>}u\mathrm{<mo\; stretchy="false">)</mo>}$, $\beta \mathrm{<mo>=</mo>0,5}$, $u_4\mathrm{<mo>=</mo>20}$.

Таблица 1. Результаты численных экспериментов

Время (сек)	Радиус зоны криопоражения (мм)	Радиус замороженой зоны (мм)	Радиус зоны влияния холода (мм)
50	9,8	10,5	13,0
100	11,3	13,4	21,1
150	12,5	15,2	26,3
200	13,3	16,1	30,2
250	13,8	16,7	32,7
300	14,0	17,0	34,0

6. Выводы.

1. Результаты проведенных численных расчетов позволили обнаружить наличие реально наблюдаемого на практике эффекта пространственной локализации теплового поля. Кроме того, рассматриваемые модели учитывают тот факт, что замерзание межклеточной и внутриклеточной жидкости происходит при существенно различных температурах и в разные моменты времени.
2. Время выхода решения на стационар приблизительно равно 5-6 минутам.
3. Использованная методика сквозного счета c предварительным сглаживанием разрывных функций применима без изменений для численного исследования задач типа Стефана любой размерности, возникающих в криохирургии.
4. В результате счета обнаружена связь между шагами по пространству и по времени. Этот факт косвенно указывает на то, что построенные разностные схемы являются условно устойчивыми. Кроме того, численные эксперименты показали возможность выбора <<большого>> шага по времени (примерно на два порядка больше шага по пространству), что особенно важно при решении многомерных задач.
5. Итерационный процесс метода Ньютона сходится достаточно быстро: 5-7 итераций при выборе функций $w_1$, $w_2$ и 3-5 итераций для функции $w_3$. Однако факт сходимости не всегда имеет место. В таких случаях итерационный процесс приходится останавливать после нескольких итераций.