Двумерная краевая задача типа Стефана для полукольца Текст научной статьи по специальности «Математика»

МЕХАНИКА

УДК 517.549.8

ДВУМЕРНАЯ КРАЕВАЯ ЗАДАЧА ТИПА СТЕФАНА ДЛЯ ПОЛУКОЛЬЦА © 2007 г Б.К. Буздов

The cylindrical-shape cryoinstrument cooled biological tissue cold spreading model is proposed in the paper. The model is the two-dimensional Stephan's type boundary problem. The algorithm of its sufficiently complete analysis has been worked out. The computer program implementing the proposed algorithm has been developed. The program can be used by medical men in performing cryosurgical operations.

В медицине как метод лечения известна гипотермия - понижение температуры биологических тканей. Широко применяются низкие температуры при консервации и хранении биоматериалов - крови, костного мозга, отдельных органов. В последние десятилетия возникла новая область медицины - криохирургия, использующая криогенные температуры -120 °С и ниже. Цель криохирургии состоит в уничтожении клеток в локальном, четко ограниченном объеме биоткани, занимаемом злокачественной опухолью. Гибель клеток достигается в результате разрывов мембран, образующихся при криогенном охлаждении, кристаллами льда внеклеточной и внутриклеточной воды, а также осмотического разбухания при оттаивании биоткани. К основным достоинствам криохирургического метода лечения следует отнести локализацию, препятствующую миграции злокачественных клеток из разрушаемого объема, что обычно имеет место с кровотоком при хирургическом вмешательстве с помощью скальпеля. Локальные криопораженные участки отмирают и отвергаются биотканью, и основная задача состоит в контроле за гибелью всех злокачественных клеток в данном объеме. Очевидно, что развитие и внедрение криохирургического метода в широкую медицинскую практику во многом зависит от достоверного описания теплового процесса замораживания биоткани, сопровождающегося фазовым переходом вода - лед. По динамике изотерм криопоражения (-20: - 50 °С) и замораживания (0: - 3 °С) скорости понижения температуры, времени достижения заданного или стационарного состояния (экспозиции) и другим параметрам криогенного воздействия на биоткань можно прогнозировать его результаты и одновременно получить необходимые данные для расчета различного криохирургического инструмента и оборудования. В связи с этим наряду с экспериментальными исследованиями весьма актуальным является математическое моделирование тепловых процессов в замораживаемой биоткани, требующее разработки эффективных методов решения задач Стефана, отличительной особенностью которых является пространственная локализация и существование стационарных решений.

Методы решения задач типа Стефана можно разбить на две группы - с явным выделением границы раздела фаз и метод сквозного счета. Методы, относящиеся к первой группе [1], как правило, весьма

громоздки и малоэффективны при решении многомерных и многофазных задач, хотя они и более точно находят границу раздела фаз. В настоящей статье применялся метод сквозного счета, идея которого восходит к работам [2, 3].

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

Уравнения, описывающие процесс распространения тепла в указанной выше области, будут следующими:

1 д ( „ , , ди 1 1 д

r dr

(r Mu) du ] + Г2 ik(u )-c(u )p(u =

= -w(u)+ P —— S(u - и *) + Pq — ö(u - un), dt dt

0 <ф<л, r0 < r < R,

u(r, (р,0) = — = const, r0 < r < R, 0 <ф<п,

l(u)— -au = au^, r = r0, 0 <q><n, t > 0, dr

(1)

t(r,p, t) = i

r = R, 0 <Ф<П, t > 0,

А(и)1 -ди -уи Го < г < Я, ф = 0, t > 0,

г дф

./ \1 ди

- А(и)---уи = и, г0 < г < Я, ф = п, t > 0,

г дф

и(г,ф* (г^),t) = и * 1

и(г,фп (гЛ),t) = иП'

Здесь А. (и), с (и), р(и) - коэффициенты теплопроводности, теплоемкости и плотности биоткани соответственно; ф*(г^), фп(г^) - изотермические поверхности, на которых температура биоткани равна и * и и П соответственно; г0 - радиус цилиндрического криозонда; Я - некоторая известная постоянная величина, характеризующаяся тем, что вне полукруга радиуса Я температура биоткани постоянна и равна и ; Р,Р0,и0- заданные постоянные; г, ф- полярные радиус и угол.

(2)

Определению подлежат функции и(г,р, t) , р* (г, t), фп (т^). Будем считать, что коэффициенты

с(и), р(и), А(и) могут иметь разрывы 1-го рода при

*

и = и и и = и п и что выполнено условие

с(и) > Сшп > 0, Р(и) > Ртп > 0, Л(и) > 2т;п > 0.

Для решения поставленной выше задачи предлагается следующий алгоритм. Введем функцию

и

Н(и) = I с(£)+ Рп(и - и*) + Р0п(и - ип. 0

Уравнение (1) примет вид

рассматривать функцию дискретных аргументов уЩ, значения которой вычисляются в узлах сетки Ж. Обозначим Жгр = Жг ■ Для сеточной функции уЩ

получаем следующую схему: -B

. х Vn+1/2 vn

Ur2)) ,J

(r, 9J )eWn

T

j vNj112 = u,

2 V1,2)

0 < j < M,

= AV+1/2 +_ w 1 2

n+1/2 0 J

fe1/2)

ahr + A^ o j

j)

+v!+1/2 +

hrauA

ochr + A^ o j

j)'

1 f f rA(u )f0 + -L ff a(u )9 + w(u ) = ^.

r dr f dr J r 2 09 f 09 J dt

Далее проводим сглаживание функций H (и), А(и), после чего получаем последовательность ограниченных гладких функций H к (и), А (и), сходящихся при к и и Ф ип, и Ф и соответственно к H (и),

А(и) [1], причем Вк(и) = нк (и)> с = const > 0.

Задачу Стефана (1), (2) можно заменить аппрокси- 2 мирующей «сглаженной» задачей, для которой и бу-

0 < J < M,

Л vn+1/2 = 1 1 '1 "r

ri+1/2A

(n+1/2 ) V+1/2, jf-

n+1/2 n+1/2Л ,+1, J

-V,,

- r,

A(n+1/2 )*

1 -1/2AVi-1/2, j) *

n+1/2 ..n+1/2

,, J

-V

i-1,J

, 1 \,,n+1 ..n+1/2 , , v

2 + P^V— = Л Vj+1 + 2 w(+1) (,9j)

+ —wV 2

:W,

r9

дет строиться численный алгоритм:

1

r dr

r aA(uk)

du dr

к Л

r 2 89

(

Ак ()■

к Л

du 89

V n+1 V10

n +1

A^r1) , .. - I n+Л Vi1

fifty+Ау,0 ) riftv+A

rV,1+1 +-

h9Yucri

W),

a(vm+1)

B ( = -w(uk )

n +1

h9Yucri

r0 < r < R,

(3)

u

0 <9<n, t > 0,

к(r,9,0) = u = const, re[r0,R], 0 <9<n,

л vn+1 = r^h1

a( )- r,ft9V U m-1 a( )-r,ft,

I \ v" (n+1 \i, j+1/2/

0 <i< N,

0 <i< N,

v n+11 V1

к

Ак (uк — auл = auA, r = r0, 0 <9<n, t > 0,

к

dr

-A(+1 J-

9

vl+1 V-1 ^

' (r,9, t )= u, r = R, 0 <9<n, t > 0,

duk к

-yu =yuc, r0 < r < R, 9 = 0, t > 0,

v(l} = u, t = 0, 0 < i < N, 0 < j < M. Здесь

r d9

v,+1/2, j = 0,5(vi, j +v,+1

, J 1

i,j+1/2

= 0,5V,

J +Vi, J+1i

«к[ к\1 du к ^ ^u t^n г,+1/2 = 0,5(г, + г,+1), г,-1/2 = 0,5(г, - г,-1).

-A lu I---yu =yuc, r0 < r < R, 9 = л, t > 0. ,+1/2 ' vi ,+u' , 1/2 ' vi , 1

r 89

У функций Bk, ик, Ак для наглядности индекс к

Условия (2) на изотермических поверхностях при ^ „ ,

4 ' * /1\ тт опущен. Для нахождения значений сеточной функции

такой постановке содержатся в уравнении (3). Для решения данной задачи применялся локально одномерный метод [4].

V на (п +1) -м временном слое по известному значению на п -м необходимо последовательно решать две

Прежде чем выписывать разностную схему, вве- серии одномерных задач с°°тветственн° го к°°рди-дем в области пространственно-временную сетку. Для гатам г и р. Каждая такая задача представляет с°-простоты рассматриваем равномерные сетки по каж- бой нелинейную алгебраическую систему с трехдиа-

дому из направлений.

Обозначим: Wr ={ri = ihr, i = 0,1,...,N - 1,N},

где

rN = R, ro = ro; Wp = \pJ = jhv, j = 0, 1,..., M -1, M} где <p0 = 0, фм =n; Wt ={tn = пт, n = 0,1,...,}, где to = 0.

гональной матрицей. Для ее решения использовался метод прогонки совместно с итерационным методом

Ньютона. При этом при определении Vп+1 коэффициенты с, р, Я брались на предыдущей итерации [3]. Вопрос о сходимости итерационного процесса рассматривался ранее в [5].

Отличительными чертами примененного для счета

Тогда на множестве W = Wr - Wt

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

что область определения сглаживающей функции

h

r

h

r

ч>

<Р

u

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

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

Выводы:

1. Результаты счета вполне согласуются с результатами других двумерных задач [6, 7], а также с результатами, полученными при рассмотрении одномерных нестационарных задач типа Стефана, возникающих в криохирургии [8].

2. Анализ выполненных расчетов позволил проследить стабилизацию перехода как поля температуры, так и поля изотерм к предельному стационарному состоянию.

3. В результате счета обнаружена связь между шагами по времени и пространству, что говорит о том, что полученная разностная схема, по всей видимости, является условно-устойчивой.

4. При исходных данных г0 = 1мм (радиус цилиндрического криоинструмента) и а =-90 °С ; ^(и) =

= (и - и)1 ¡2 ; (значения исходных данных, характеризующих биоткань, можно посмотреть в [9]) получены следующие результаты: стабилизация к предельному стационарному состоянию достигается примерно через 5-6 мин; при этом радиус зоны поражённой биоткани колеблется около отметки 1,6 мм; замороженной биоткани - « 1,9 мм, охлаждённой биоткани - и 35 мм. На рис. 1, 2 можно проследить динамику распространения температурного поля в биоткани, а также границы зон криопоражения, замораживания и влияния холода в различные моменты времени.

Рис. 1. Динамика распространения температурного поля

0,10

0,05

9,00

50 с

-1-1- -1- -1-1-

- 0,10 -0,05 0,00 0,05 0,10

0,10 -

0,05

0,00

300 с

-0,10 -0,05 0,00 0,05 0,10

Рис. 2. Изотермы криопоражения замораживания и влияния холода в разные моменты времени: 1 - изотерма влияния холода; 2 - изотерма замораживания; 3 - криопоражения

Литература

1. Мажукин В.И., Самарский А.А., Чуйко М.М. // Докл. РАН. 1999. Т. 368. № 3. С. 307-310.

2. Будак Б.М., Васильев Ф.П., Успенский А.Б. // Численные методы в газовой динамике. Вып. IV. М., 1965. С. 139-183.

3. Самарский А.А., Моисеенко Б.Д. // ЖВМ и МФ. 1965. Т.5. № 5. С. 816-827.

4. Самарский А.А. // ЖВМ и МФ. 1963. Т. 3. № 3. С. 431-466.

Кабардино-Балкарский государственный университет

5. Буздов Б.К. // Нелинейные краевые задачи математической физики и их приложения. Киев, 1993. С. 33-35.

6. Буздов Б.К., Буздов А.К. // Системные проблемы качества, математического моделирования, информационных, электронных и лазерных технологий. Материалы Междунар. конф. и Рос. науч. школы: М., 2002. С. 1416.

7. Буздов А.К., Буздов Б.К. // Изв. КБНЦ РАН. 2005. Т. 2. № 14. С. 1-5

8. Березовский А.А., Жураев К.О., Юртин И.И. // Задачи Стефана со свободными границами. Киев, 1990.

9. Березовский А.А., Леонтьев Ю.В. // Криобиология. 1989. № 3. С. 7-13.

_29 мая 2006 г