Численное решение двумерных краевых задач типа Стефана в криохирургии тема автореферата и диссертации по математике, 01.01.03 ВАК РФ

л0.

НАЦИОНАЛЬНАЯ АКАДЕМИЯ НАУК УКРАИНЫ ОРДЕНА ТРУДОВОГО. КРАСНОГО ЗНАМЕНИ ИНСТИТУТ МАТЕМАТИКИ

На правах рукописи

Буздов Беслан Каральбиевич

УДК 517.946.9

ЧИСЛЕННОЕ ФЕШЕНИЕ ДВУМЕРНЫХ КРАЕВЫХ ЗАДАЧ ТИПА СТЕФАНА В КРИОХИРУРГИИ

(01.01.03 - математическая физика)

Диссертация на соискание ученой степени кандидата физико-математических наук

■ - и д

НАУЧНЫЕ РУКОВОДИТЕЛИ:

, ^иум ВАК рОсс®^Й^31Ш0_математических

- даукД профессор

:;лдатъ Диплом ~ ^ "

рктф физико-математических Росащки профессор

I М.Х.ШХАНУКОВ

клндидл^

А.А.БЕРЕЗОВСКМИ

-Ш-Ще: у1:

Киев - 1994

Содержание.

Введение ......................................................3

Глава I. Одномерные задачи типа Стефана в криохирургии.........20

§1.0 задачах Стефана....................................20

§2. Метод квазилинеаризации для нелинейных краевых

задач ...............................................25

§3. Задача плоско-параллельной криодеструкции биоткани...29

Глава II. Двумерные задачи типа Стефана в криохирургии.........42

§1. Локально-одномерный метод...........................42

§2 Двумерные постановки задач типа Стефана с явным

выделением границы раздела фаз......................48

§3. Двумерная краевая задача типа Стефана в прямоугольной области....................................53

§4. Двумерная краевая задача типа Стефана в

полукольце..........................................63

§5. Априорные оценки для решения параболических уравнений, возникающих в криохирургии.........................71

Приложение 1. Тексты программ.................................78

Приложение 2. Графики.........................................90

Литература....................................................95

БВЕДЕНИЕ

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

-3°С), скорости понижения температуры, времени достижения

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

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

криозонда фронт замораживания распространяется по ткани до некото poro предельного положения. Соответствующее предельное положение изотермы криопоражения определяет максимальный размер деструкции биоткани. На динамику криопоражения влияет геометрия охлаждающей поверхности криозонда и ее температура, теплофизические характеристики замороженной и незамороженной ткани,метаболическая скорость теплообразования ткани, температура крови, скорость кровотока в ткани, условия теплообмена на поверхности ткани.

Решение задач типа Стефана сопряжено со значительными трудностями в связи с их существенной нелинейностью и требует привлечения новых идей,использования всего арсенала конструктивных методов нелинейного анализа и возможностей современной вычислительной техники.

В последнее время при решении задач Стефана применяются различные приближенные аналитические методы. К таким методам следует отнести методы сведения к системе обыкновенных дифференциальных уравнений [ 1,9 ], метод интегральных уравнений [5,55] эффективность которого известна лишь в линейном случае, методы теории потенциала [61,73], метод малого параметра [37], метод Роте [46], метод нелинейных вариационных параметров [72] и др..

Однако в связи с быстрым совершенствованием вычислительной техники все большее развитие получают эффективные численные методы решения таких задач. К ним относится метод прямых и метод сеток [66], проекционно сеточный метод [40], развитый в работах [49,66],метод декомпозиции И метод фиксированных областей [ 45 ] и

др.. Вообще,все существующие в настоящее время методы численного решения задач типа Стефана можно условно разбить на два класса: с явным выделением фронта и сквозного счета. Численные алгоритмы, основанные на явном выделении фронта, см. [3,24]»обладают высокой точностью определения межфазной границы, но как правило весьма громоздки и для получения решения требуют больших затрат машинного времени. По этим причинам они зачастую неприемлимы для многомер- * ных и многофазных задач. Среди методов решения задач Стефана с явным выделением фронта раздела фаз следует отметить разностный метод с ловлей фронта в узел сетки, подвергнутый достаточно полному теоретическому исследованию в работах [22,25],однако он также мало пригоден, когда искомые функции зависят от нескольких пространственных координат. Также хотелось бы отметить работу [47], в которой предлагается новый подход к численному решению двухфазной задачи Стефана. В основу метода положена идея построения и использования адаптивной сетки. В предлагаемом подходе осуществляется автоматическое выделение фронта и перестановка расчетной сетки. Метод также учитывает кинетику фазовых превращений, что бывает важным во многих задачах.

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

фана с помощью экономичного алгоритма сквозного счета [23], эффективность которого особенно заметна в многомерных постановках,так как алгоритм не зависит ни от числа фазовых фронтов , ни от размерности задачи. Недостатком указанного метода является сравнительно низкая точность определения положения фазового фронта, а также чувствительность к выбору параметра сглаживания , определить значение которого априори в ряде случаев затруднительно.

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

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

Диссертация состоит из введения, двух глав, двух приложений и списка литературы.

В первой главе основное внимание уделено одномерной начально-краевой задаче типа Стефана и методу ее решения, который в дальнейшем В значительной степени используется при решении двумерных задач.

В первом параграфе дается краткая общая характеристика задач типа Стефана и вопросов, связанных с их решением.

Во втором параграфе рассматриваются вопросы сходимости метода квазилинеаризации [4] для следующих нелинейных краевых задач :

f Ы = - д(и)и = /(и), 0 < X < Ъ,

йх (Их

и(0) = и(Ь) = 0, к(х) > 0, д(х) > 0;

(1)

кроме задачи (1) рассмотрена задача :

ди =, сГи

at дз?

+ Г(и),

и(о,г) = и(1,г) = о,

„ и(х,0) = о.

(2)

Задача (2) решается в области £> = ((х^) :СКХ<1 . Соглас-

но методу квазилинеаризации для обеих задач строится итерационный процесс, который сходится при достаточно малых Т. При этом показывается , что скорость сходимости итерационного процесса является квадратичной. Задачи подобные (1),(2) возникают при решении двумерных задач криохирургии , что вызвало необходимость их рассмотрения .

В третьем параграфе решается следующая начально-краевая задача плоскопараллельной криодеструкции биоткани :

(\(и)—) - с(и)р(и)^ = -т(и) - Р ^*д(х-х*) дх дх дг сН

= и = соп^, х > О,

(3)

Х(и)— - <ш = -оси,., х = 0, í > О,

Ит и(х^ )■ = и, t > О х-*»

К написанным выше уравнениям добавляются следующие условия изотермично с ти :

и (хи(Х),г) = ип , и (х*(г)^) = и*,

х*а) = о, * ^ , япа; = 0, > Г4;

Здесь Я = Х(и), с = с (и), р = р(и) - коэффициенты теплопроводности, теплоемкости и плотности биоткани; ш = и>(г0 источники тепла; X*(t) и г = жп(I) - координаты изотермы замораживания и = и и изотермы криопоражения и = ц^; Р = р Лв; Лв - скрытая теплота кристаллизации воды ; р , рп - массивы внеклеточной и внутреклеточной воды в единице объема биологической ткани, 0(х) - дельта функция Дирака, й- начальная температура, ос - коэффициент теплообмена, иА= - температура охлаждающей

поверхности плоского криоинструмента; Ц Х2- моменты времени при которых поверхность биоткани х=0 охлаждается соответственно до температуры замараживания и(0,^)=и* и криопоражения и(0^2) = Цд. Определению подлежат функции и(х,Х), х*^)-, х^Х).

Для решения задачи (3), (4) вводится функция теплосодержания :

Н (и) = X са)р(^т + Рг)(и-и*) + р(и-ип) после чего задача (3) , (4) переписывается в виде :

Ш(и), 0 < X < а? , г > 0,

В

0 ^ х ^ я ,

£5

X = 0, £ > 0, (5)

£ > 0,

где согласно [11]

, й X (и) <Ш

X = -1 £ -

в _ 0 _

/2 - и 1/2

(7 ш(и)Х(и)Ои) и

а (Ш)т) .щи)..

дх дх <Э£ и(х,0} = и = сопзг ,

\(и)— - оси = - оси,, дх А

и(хв,г) = и,

Для решения задачи (5) применяется метод [ 22 1, основанный на сглаживании функции теплосодержания Н(и) и коэффициента теплопроводности Х(и) . Сглаженная задача решается с использованием чисто неявной разностной схемы вида (1.3.15).

Нелинейное относительно уравнение решается итерацион-

ным методом Ньютона , где каждая из итераций находится методом прогонки .

Глава 2 настоящей работы посвящена решению двумерных краевых задач типа Стефана , возникающих в криохирургии.

В §1 описан локально-одномерный метод, применяемый при решении многомерных задач и используемый в дальнейшем. Метод предложен и обоснован в работах [63,64]. Он состоит в поэтапном решении по разным пространственным переменным одномерных уравнений теплопроводности, при помощи устойчивых неявных схем. Локально-одномерный метод пригоден для произвольных областей в случае краевых условий I рода,а в случае условий 3-го рода- для областей специального вида (см.[65]).

В §2 дается краткий обзор двумерных задач типа Стефана, рассмотренных другими авторами и имеющих непосредственное приложение в криохирургии.

В §3 рассматривается двумерная начально-краевая задача типа Стефана , возникающая в криохирургии в случае достаточно протяженного плоского криоинструмента . Задача решается в прямоугольной области, а соответствующие уравнения выглядят следующим образом:

—Гк(и)—] + -дщи)^] - с(и)р(и№ дх дх ду ду

-щи) + р 6(и _ и*} + р ди д(и _ }

at 0 дг ^

х € (-а,а), у € (0,Ъ), г > 0, (6)

и(х,у,0) = и = co7гsí,

Х(и- оси = оси., у = О, X € /"-Г_, Г > О,

ду А оо

Л-Ги;^ - Ти = Ти_, У = 0, х £ [~гп, г], í > О, ду 0 00

и(-а,у,г) = и , у € /"0,Ь7, t > О, и(а,у,г) = и , у € [0,Ы, г > О, и(х,Ъ,1) = и , я € £-а,а], t > О,

и(х,у*(х,1),г) = и*,

и(х,ул(х,г),г) = ип.

Здесь 2г0 - ширина плоского криоинструмента; Т коэффициент теплообмена с окружающей средой ; а, Ь известные параметры,

характеризующиеся тем , что вне рассматриваемого прямоугольника температура биоткани постоянна и равна и.

Определению подлежат функция температуры u = и(х,у,t), а

У

также пара изотермических поверхностей у (x,t), у (x,t), на

которых температура биоткани равна , соответственно , и и ип •Введя функцию теплосодержания , Н(и) и произведя процесс сглаживания , задача (6) заменяется аппроксимирующей

"сглаженной" задачей вида :

J? (Ш)ди } + _д (Ш)ди } _ Щи1 = _ w(u)t дх дх ду ду ■ dt

х € (-а,а), у € ГО,Ы, t > О, (7)

и(х,у,0) = и = const, х € [-а,а], у е [0,Ы,

Х(и)~ - du = du., у = О, х € t-rn, rj, t > О, ду о о

Х(и)— - Ти = Ти„, у = 0, х € i-af-r.mi г ,al, t > О, ду с оо

u(-a,y,t) = и , у € [0,ЪJ, t > О, u(a,y,t) = и , у € [0,Ы, t > О, u(x,b,t) = и , х € [-a,a]f t > 0.

Для задачи (7) уже можно строить локально-одномерную разностную схему сквозного счета , так как фазовый фронт в такой постановке явно не выделяется.

В этом же параграфе для задачи (7) предлагается способ нахождения граничного условия на отрезках у = О,х € [-a,-rQ]¡J[ r0,aJ , когда коэффициент теплообмена Т является неизвестным. Если считать, что распространение тепла вдоль поверхности биоткани происходит с большей скоростью, чем по глубине, то представляется возможным неизвестное граничное условие ( для определенности рассматривается отрезок х € [rQ, aJ) находить как решение следующей одномерной начально - краевой задачи :

(Х(и) - _ W(U), х € fr_, a), t > О,

дх дх dt °

и(х,0) = и = const, х € írQf а],

%(и)— - <ш = - oíuA, х = rn, t > О, (8)

дх А 0

u(a,t) = и , t > 0.

Отметим, что задачи типа (8) необходимо решать на каждом шаге времени в исходной двумерной задаче. При этом граничное условие

3-го рода в исходной постановке заменяется на условие 7-го рода , что ускоряет процесс решения на ЭВМ . Скорость счета существенно увеличивается и за счет более быстрой сходимости итерационного процесса при решении нелинейного разностного уравнения методом Ньютона . Это связано с тем , что данные на границе становятся более "согласованными" с неизвестными значениями искомого решения в исследуемой области.

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

Предлагается следующая постановка задачи :

— — (г \(и)~) + — (Г Ш)—)~ С(и)р(и0 г дг дг г бф <Эф дг

= - ю(и) + Р — 6(и - и*) + Рп — д(и - и),

дг 0 дг ■ п

г0< г<Я,0<(р<%,г>0, (9)

и(Г,<р,0) = и = С0П8г, Г0 < Г ^ К , О^ф^ТС,

А- оси = оси., г = гп , О ^ Ц) ^ % , г > О,

дг А °

и(г, = и, г = й,0^ф<%,£>0,

^ - ти = тип, гп г ^ я , ф = о , г > о,

г д Ф со

~Х(и)~ — - Ти = Ти_, г * г й , ф = тс , £ > О, г (Эф со

и(г(г),г) = и*, игг,фпгг,г;д; = ип.

Здесь г - полярный радиус , ф - полярный угол ; г0 - радиус цилиндрического криозонда ; й - некоторая известная постоянная , характеризующаяся тем , что вне полукруга радиуса К температура биоткани постоянна и равна и; ф, ф -

пара изотермических поверхностей , на которых температура би-

* ^

откани постоянна и равна и и ип, соответственно . В остальном обозначения те же , что и ранее . Определению подлежат функции и(Г,ф,£) , <р*(г^) ,'ф . Метод решения задачи

(9) тот же , что и для задачи в прямоугольной области . После ввода функции теплосодержания и сглаживания получается следующая "сглаженная" задача :

Л (Г Х(и+ -I (г Х(и)^)~ ^МШ = - ю(и) , г дг дг г дф аф д1

г < г < Я , 0<ф<тс,1;>0,

(ю)

u(r,<p,0) = и = const, г € L'r ,R] , О < ф < ти , t > О,

Яц

k(U)— - СЙА = сСИА , Г =

дг А 0

ufr,9,tj = и, г = Д,0<ф<1и, t > О,

X(U)— ^ - Тu = гл < г < Д , Ф = О , t > О,

г 9ф со

-Х(и)— — - Ти = г « г < Я , ф = тс , i > О.

г 5ф с 0

Для задачи (10) далее применяется локально-одномерная разностная схема сквозного счета .

Краевые условия 3-го рода при ф = 0 ( и ф = 7U ) здесь также можно заменить на условие 1-го рода, которое находится реш�