Программный комплекс rtf для расчета нестационарных радиационных и тепловых полей, возбуждаемых лазерным облучением, в гетерогенных биологических тканях Текст научной статьи по специальности «Физика»

Вестник Челябинского государственного университета. 2015. № 7 (362).

Физика. Вып. 20. С. 64-71.

БИОФИЗИКА

А. С. Зарезина, А. В. Лаппа

ПРОГРАММНЫЙ КОМПЛЕКС RTF ДЛЯ РАСЧЕТА НЕСТАЦИОНАРНЫХ РАДИАЦИОННЫХ И ТЕПЛОВЫХ ПОЛЕЙ, ВОЗБУЖДАЕМЫХ ЛАЗЕРНЫМ ОБЛУЧЕНИЕМ, В ГЕТЕРОГЕННЫХ БИОЛОГИЧЕСКИХ ТКАНЯХ

Представлен программный комплекс для исследования радиационных и тепловых полей в биологических тканях, подвергаемых лазерному облучению. В программном комплексе реализованы: двумерная осесимметричная гетерогенная модель среды; кинетическая модель переноса излучения и метод Монте-Карло для моделирования радиационных полей; нестационарное «биотепловое» уравнение; метод конечных элементов и метод взвешенных невязок для моделирования тепловых полей.

Ключевые слова: лазер, биологическая ткань, излучение, температура, метод Монте-Карло, метод конечных элементов.

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

К настоящему времени в мире создано несколько программ для расчета радиационных и тепловых полей лазерного излучения в биологических тканях [1-4], но они практически недоступны для стороннего пользователя. Кроме того, предлагаемые в этих программах модели не дают полного решения поставленных проблем. Это связано, в первую очередь, со сложной структурой биологической ткани, а также с большим многообразием физических явлений, происходящих внутри ткани при взаимодействии с лазерным излучением: нагрев, тепломассоперенос, теплообмен с окружающей средой, фазовые изменения в структуре ткани и т. д.

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

поля, второй — моделирование тепловых полей с известными источниками тепла, вычисленными на предыдущем этапе.

Модель реализована в виде программного комплекса RTF (Radiation-Thermal Fields), снабженного интеллектуальным интерфейсом для ввода и хранения информации, а также для графического представления получаемых результатов. Комплекс адаптирован для работы на современных платформах (Windows 7, Windows 8).

Модель среды

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

Моделирование радиационных полей

Перенос лазерного излучения в биологических тканях рассматривается в рамках квазистационарной кинетической модели переноса излучения. Распространение излучения в ней трактуется как многократное рассеяние частиц на случайных неоднородностях среды. Традиционно в теории переноса эти частицы принято называть фотонами (этот термин будет использоваться и в настоящей работе).

г, см

Рис. 1. Модель гетерогенной биологической ткани

Основной характеристикой поля излучения является дифференциальная интенсивность излучения /(г,й). Предполагается, что она удовлетворяет уравнению переноса излучения:

ЙУ/(г,й) + ц(г)/(г,й)--ц,(г)|р(г, й'->й)/(г,йуй' = 5(г,П). (1)

Здесь 5"(г,й) — дифференциальная плотность мощности излучения (функция источника); p — индикатриса рассеяния (плотность распределения фотонов по углам рассеяния); ц = ц + ц ц ца — макроскопические сечения (см-1) взаимодействия, рассеяния и поглощения соответственно.

Индикатриса рассеянияp(r, й' ^ й) в приближении азимутальной симметрии лазерных пучков и изотропности среды является функцией только угла рассеяния 0 (cos9 = й'й):

p(r, й' ^ й) = p(r, й'й) = p(r,cos0).

В качестве нее используется традиционная для биооптики функция Хени — Гринштейна [5]:

p(cos 0) =

1 - g2

4 п(1 -2 g cos 0 + g2)

2Л3/2

(2)

единственный параметр которой g — средний косинус угла рассеяния

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

Определение характеристик радиационного поля лазерного излучения основывается на решении уравнения переноса излучения методом Монте-Карло [6], то есть моделируются траектории фотонов, начиная с их рождения источником и заканчивая поглощением средой. При этом используются действительные оптические параметры среды: ц^, ц и ^(соэ0). При пересечении фотоном границы раздела зон с разными показателями преломления п происходит отражение или преломление в соответствии с законами Френеля.

В настоящей работе рассчитываются следующие характеристики радиационных полей:

- плотность поглощенной мощности, равная количеству энергии, поглощенной единицей объема среды в единицу времени; эта величина является тепловым источником, генерируемым лазерным излучением, при моделировании тепловых полей:

б(г) = цв(г)|/(гЛ)Л =иа(г)-/(г).

4л

- плотность тока отраженного Q (r) и прошедшего Q (г) излучения:

б»(г) = \ /(г,й)^й, й(г)= | /(Г,й)с/й.

(Яп)<0 (£1п)>0

Эти величины используются при экспериментальном определении оптических параметров биотканей.

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

Моделирование тепловых полей

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

В качестве тепловой модели в настоящей работе принята модель «биотеплового» уравнения:

ф)р(г) = У(к (г )УТ (г, г)) -

дг

-a(r)(T(r, t) - Tb) + Q(r, t), r e G;

(3)

б(г,0 = Цв(г)//(г,П)</а ге£; (4)

Т (г,0) = Г0(г); (5)

5иД (Т(г,г) - Т (г, г)) + 5я>2 (*(г) • (п • УТ(г,г)) --?(!■)) + (г) • (п • УТ (г, г)) + +А(г) • (Т(г, г) - Тз (г, г))) = 0, г е Г; (6) к а (г)п-УТА (г, г) --кв (г)пав -УТ (г, г) = 0, г еГ, Та (г,г)-Тв (г,г) = 0, г еГ'. (7)

Здесь:

G — множество в евклидовом пространстве М3, занятое средой, в котором рассматривается перенос тепла; множество G не включает границы ни внутренние (между смежными зонами), ни внешние (с окружающей средой), является открытым;

Г — внешняя, Г' — внутренняя границы; О и Г' и Г= О;

Т(г,Р) — температура среды в точке г в момент времени Р;

с Ть — удельная теплоемкость при постоянном давлении и температура крови.

а = 0с4иг>; — «объемный расход» кровотока (характеризует теплообмен «кровеносная система — среда» и равен массе крови, протекающей сквозь единицу объема среды в единицу времени);

0 — феноменологический параметр, характеризующий теплообмен между сосудом и тканью;

k — коэффициент теплопроводности среды;

с, р — удельная теплоемкость при постоянном давлении и плотность среды;

Q(г,t) — плотность мощности лазерного излучения, поглощаемой средой в точке г в момент времени Р;

ТА, kA, Тв, kB — температура и коэффициент теплопроводности в смежных зонах А и В (с общей границей Г');

п, пав — внешняя нормаль к Г и нормаль к Г' (направление от зоны А к зоне В);

т — признак граничного условия, заданного на некотором подмножестве границы

Надлежащим выбором параметра т возможен учет следующих физических условий на внешней границе Г:

- т = 1: фиксированная температура Т1(г) (граничное условие 1-го рода);

- т = 2: тепловой поток д(г) (граничное условие 2-го рода);

- т =3: теплообмен по закону Ньютона (граничное условие 3-го рода), ^г) — коэффициент

теплоотдачи, Т3(г) — температура окружающей среды вблизи границы

Условия (7) обеспечивают физические требования на внутренних границах: непрерывность температуры и плотности потока тепла через границу. Предполагается, что параметры с, р, к постоянны внутри отдельной зоны; сь, иь, Ть — постоянны во всей области G.

Принятая модель переноса тепла позволяет использовать для моделирования ту же модель среды (рис. 1): набор осесимметричных однородных кольцевых зон прямоугольного сечения со своими оптическими и теплофизическими параметрами.

Численное решение нестационарного уравнения теплопроводности осуществляется методом конечных элементов для пространственной части уравнения в сочетании с методом взвешенных невязок для решения временной части [7].

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

Результаты

Программный комплекс для расчета радиационных и тепловых полей лазерного излучения прошел ряд тестовых расчетов, включающих сравнение с известными аналитическими решениями, с результатами других программ, а также с результатами экспериментов [8].

Некоторые возможности нашего комплекса демонстрируются ниже.

Моделирование тепловых полей излучения в многослойной биоткани

На рис. 2 приведены результаты моделирования теплового воздействия лазерного пучка на 4-слой-ную модель кожи.

Образец ткани представлял собой цилиндр бесконечного радиуса с теплоизолированным верхним торцом. Расположение нижнего торца соответствовало глубине, ниже которой изменение температуры, связанное с лазерным облучением, было пренебрежимо малым (или отсутствовало). Поэтому было принято, что температура на нижнем торце оставалась постоянной и равной 37 °С. Оптические и теплофизические параметры слоев ткани приведены в табл. 1 [9]. Начальная температура для всей ткани была принята одинаковой и равной 37 °С .

Д Z, см

Лазерный пучок

10 14 Расстояние от оси пучка, хЮ~2 см

(4) дерма

(3) кровь

(2) дерма I I

1 (1) эпидермис |

г, см 0,1 см

Лазерный пучок

Рис. 2. Изотермы в коже, облучаемой цилиндрическим мононаправленным пучком Не-Ые лазера в непрерывном режиме в момент времени 1,5 с после включения лазера. Лазерный пучок: диаметр равен 0,1 см, мощность равна 0,25 Вт. Начальная температура равна 37 °С . Параметры 4-слойной модели кожи приведены в табл. 1 для каждого слоя

Таблица 1

Оптические и теплофизические параметры слоев ткани

Слой Толщина, мм Ца, 1/см 1/см g n k, 103 Вт/см-K р, г/см3 с, Дж/(г-К)

1.Эпидермис 0,065 4,3 107,0 0,79 1,5 2,66 1,6 3,7

2. Дерма 0,565 2,7 187,0 0,82 1,4 4,98 1,0 3,2

3. Кровь 0,09 25,0 400,0 0,98 1,35 5,30 1,0 3,6

4. Дерма 0,565 2,7 187,0 0,82 1,4 4,98 1,0 3,2

Лазерный пучок He-Ne лазера (длина волны 633 нм) имел гауссов профиль распределения плотности мощности, диаметр 0,1 см (два стандартных отклонения), мощность 250 мВт. Режим облучения — непрерывный.

Моделирование такого рода позволяет не только найти предварительные параметры поля лазерного излучения, но и изучать, выявлять эффекты лазерного воздействия на различные биологические ткани.

Моделирование полей излучения в ткани мозга

На рис. 3-5 приведены распределение поглощенной мощности и стационарные изотермы в полубесконечной среде белого вещества головного мозга человека, облучаемой лазерными пучками различной длины волны (805 и 1064 нм).

Образец ткани представлял собой однородный цилиндр (рис. 3). Теплообмен на границах не учитывался (все границы — теплоизолированы). Оптические и теплофизические параметры ткани для соответствующих длин волн приведены в табл. 2 [9] . Начальная температура для всей ткани была принята одинаковой и равной 37 °C .

Лазерный пучок в обоих случаях имел равномерное распределение плотности мощности пада-

ющего излучения, диаметр 0,04 см, мощность 1 Вт. Режим облучения — непрерывный.

Результаты данного моделирования позволили объяснить различия в клиническом эффекте при проведении операции по лазерной деструкции корешка тройничного нерва лазерами с длиной волны 805 и 1064 нм [10].

t г 2 см w

\ Z

/ А

Т1 см

Рис. 3. Схема облучения

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

-г"^—I-г

1 1.5 2 3 Расстояние от оси Лазерный пучок пучка, х10см

Лазерный пучок

1 1,5 2 3 Расстояние от оси пучка, хЮ"1 см

Рис. 4. Пространственное распределение плотности поглощенной мощности в белом веществе мозга человека, облучаемого лазером с длиной волны X =805 нм (слева) и X = 1064 нм (справа) в полубесконечной однородной среде на различной глубине на оси пучка. Лазерный пучок с равномерным распределением интенсивности диаметром d = 0,04 см, мощностью 1 Вт. Параметры среды для соответствующих длин волн приведены в табл. 2

1 2 3

Расстояние от оси .Лазерный пучок

1 2 3 Расстояние от оси .азерный пучок ^^ х10-см

Рис. 5. Стационарные изотермы в белом веществе головного мозга человека, облучаемого лазерным пучком длиной волны (слева — X = 805 нм, справа — X = 1064 нм). Соответствующие теплофизические параметры приведены в табл. 2

Таблица 2

Оптические и теплофизические параметры ткани для моделирования полей излучения в ткани мозга

Белое вещество мозга человека и., 1/см г а и, 1/см г s, g п к, 10-3 Вт/см-К р, г/см3 с, Дж/(г-К)

1064 нм 3,2 469,0 0,87 1,45 6.48 1 4.18

805 нм 0,2 400,0 0,9 1,38

Параметры среды

Таблица 3

и , 1/см г а и, 1/см г s, g п к, 10-3 Вт/см-К р, г/см3 с, Дж/(г-К)

40,0 282,0 0,92 1,45 3,5 1,5 3,2

Параметры облучения выбирались согласно условиям проводимого эксперимента [11]. Лазерный пучок представлял собой сферический облучатель, позволяющий равномерно облучать всю поверхность полости. Облучатель имел равномерное распределение плотности мощности, диаметр 0,06 см, мощность 20 Вт. Режим облучения — непрерывный.

80

±1

z = 2 см >

Т|

0,06 см

5

7

\

J-

= 2 см

Z = 4 см -►

Рис. 6. Условия облучения

На рис. 7-8 представлены результаты моделирования распределения плотности поглощенной мощности и распределение температуры при облучении. Видно, что получено хорошее согласие температурных результатов с экспериментальными данными.

-2,0 -1,5 -1,0 -0,5 0,0 0,5 1,0 1,5 2,0

Расстояние от оси пучка, см

Лазерный пучок

Рис. 7. Распределение плотности поглощенной мощности в ткани головного мозга, облучаемой сферическим облучателем диаметром 0,06 см, мощностью 20 Вт. Параметры среды приведены в табл. 3

Полученные результаты позволили оценить площадь поверхности, подвергаемой лазерному облучению, что в свою очередь дало возможность скорректировать параметры и режимы облучения при проведении операции по удалению гематомы головного мозга и последующего облучения ее полости [11].

70

U

60

з

п. g 50

•

•

чч

•У • •V

излучение

40

30

0 10 20 30 40 50 60 70 80 90 Время, с

Рис. 8. Зависимость температуры от времени в точке на расстоянии г = 0,2 см от оси пучка на глубине г = 2 см. Сплошная линия — результаты экспериментального измерения, точки — результаты моделирования. Лазерный пучок представляет сферический облучатель

диаметром 0,06 см, мощностью 20 Вт, длительность облучения 30 с. Все поверхности образца теплоизолированы (нет теплообмена). Начальная температура равна 33 °С. Параметры среды представлены в табл. 3

Моделирование полей излучения в ткани печени

С помощью данного комплекса было проведено подробное исследование облучения ткани печени крысы Nd:YAG лазером (1064 нм) в условиях, соответствующих проведенному эксперименту [12]. Образец ткани представлял собой однородный цилиндр (рис. 9), теплообмен с окружающей средой не учитывался (торцы и боковая поверхность теплоизолированы). Оптические и теплофизические параметры ткани приведены в табл. 4. Оптические параметры ткани были определены экспериментально с помощью собственной установки [13]. Начальная температура для всей ткани была принята одинаковой и равной 37 °С .

i

2 см

0,04 см

IF

:.-

1 СМ

^ 1 см"

Рис. 9. Условия облучения для моделирования полей излучения в ткани печени

На рис. 10 приведено распределение плотности поглощенной мощности лазерного излучения. Видно, что излучение мало проникает внутрь ткани: уже на расстоянии 2 мм от лазерного пучка плотность поглощенной мощности уменьшается более чем в 100 раз. Следовательно, существенные изменения температуры будут локализованы вблизи источника, что позволяет для тепловой задачи принять границы образца теплоизолированными.

0,2 0,4 0,6 Расстояние от оси — пучка, см

Лазерный пучок

Рис. 10. Пространственное распределение плотности

поглощенной мощности лазерного излучения в однородной среде печени крысы. Лазерный пучок равномерной интенсивности, диаметром й = 0,4см, мощностью 10 Вт. Параметры среды приведены в табл. 4

2 3 4 5

Время облучения, с

Рис. 11. Зависимость температуры от времени в точке

на оси пучка (г = 0) на глубине z = 0,1 см для импульсного и непрерывного режимов облучения. Диаметр пучка 0,04 см, средняя мощность — 10 Вт для обоих режимов, длительность импульса 0,8 с, частота 1 Гц. Теплообмен с окружающей средой не учитывается. Начальная температура равна 37 °С. Параметры среды представлены в табл. 4

1,20 111111 1,20

1,10 а 1,10

1,00 1,00

0,90 0,90

0,80 0,80

0,70 / / / ч\ Б ^ 0,70

о )

0,60 0,60

0,50 \ 'Ъ ^ -_- 0,50

0,40 0,40

0,30 0,30

0,20 0,20

0,10 0,10

0,00 -1-1-1-1-1-1- —1-1- ^ 0,00

0,0 0,1 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0,9 1,0 1,1 1,2

0,0 0,1 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0,9 1,0 1,1 1,2

Рис. 12. Изотермы в ткани печени крысы, облучаемой мононаправленным лазерным пучком равномерной интенсивности, через 5 с после облучения. Режим облучения — непрерывный (а), импульсный (б). Диаметр пучка 0,04 см, средняя мощность — 10 Вт для обоих режимов, длительность импульса 0,8 с, частота 1 Гц. Теплообмен с окружающей средой не учитывается. Начальная температура равна 37 °С. Параметры среды представлены в табл. 4

Таблица 4

Оптические и теплофизические параметры ткани (печень крысы)

Биологическая ткань и , см 1 г а ц см 1 g п к, Вт/смК с, Дж/гК р (г/см3)

Печень крысы (1064 нм) 8,0 255,0 0,92 1,5 4,98-10-3 3,2 1

б

На рис. 11-12 представлено сравнение двух режимов облучения (импульсно-периодического и непрерывного) с одинаковой средней мощностью. Видно, что при импульсно-периодическом режиме облучения получен локальный нагрев ткани до более высоких температур (выше 100 °С), причем объем ткани, нагретой, например, выше 50 °С, практически одинаков для обоих режимов облучения.

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

Список литературы

1. London, R. A. Laser-tissue interaction modeling with LATIS / R. A. London, M. E. Glinsky, G. B. Zimmerman [et al ] // Appl. Opt. 1997. Vol. 36 (34). P. 9068-74.

2. Rossacci, M. J. A 3-d model for laser heating of a heterogeneous turbid medium / M. J. Rossacci, C. A. Di Marzio, S. C. Lindberg [et al.] // Proc. SPIE.

1997. Vol. 2970. P. 144-157.

3. Yilbas, B. S. 3-dimensional laser heating model including a moving heat source consideration and phase change process / B. S. Yilbas // Heat and Mass Transfer.

1998. Vol. 33. P. 495-505.

4. Щербаков, Ю. Н. Моделирование тепловых процессов при взаимодействии некоагулирующе-го лазерного излучения с многослойной тканью / Ю. Н. Щербаков, А. Н. Якунин [и др.] // Оптика и спектроскопия. 1994. Т. 76, № 5. С. 845-850.

5. Henyey, L. G. Diffuse radiation in the galaxy / L. G. Henyey, J. L. Greenstein // Astrophys. J. 1941. Vol. 93. P. 70-83.

6. Ермаков, С. М. Метод Монте-Карло и смежные вопросы / С. М. Ермаков. М. : Наука, 1975. 472 с.

7. Сегерлинд, Л. Метод конечных элементов / Л. Сегерлинд. М. : Мир, 1979. 392 с.

8. Lappa, A. V. A new computer code for calculation of radiation and heat fields in laser-irradiated tissues /

A. V. Lappa, A. S. Anikina, V. A. Kamalov // Proc. SPIE. 2000. Vol. 3911. P. 502-511.

9. Тучин, В. В. Лазеры и волоконная оптика в биомедицинских исследованиях / В. В. Тучин. Саратов : Изд-во Сарат. ун-та, 1998. 383 с.

10. Исмагилова, С. Т. Лазерная хирургия невралгии тройничного нерва : дис. ... д-ра мед. наук / С. Т. Исмагилова. Челябинск, 2000. 211 с.

11. Сорвилов, В. П. Профилактика повторных кровоизлияний расфокусированным высокоинтенсивным лазерным излучением при хирургическом лечении внутримозговых гематом (экспериментальное исследование): дис. ... канд. мед. наук /

B. П. Сорвилов. Челябинск, 2004. 146 с.

12. Бондаревский, И. Я. Опыт применения высокоинтенсивного лазерного излучения при операциях на печени (экспериментальное исследование) / И. Я. Бондаревский, В. Н. Бордуновский, Л. В. Астахова // Лазерные технологии в медицине : сб. науч. работ. Челябинск, 1999. С. 114-121.

13. Рогальский, Ю. К. Метод определения оптических параметров биоткани по пространственному распределению прошедшего образец излучения / Ю. К. Рогальский, А. В. Лаппа, М. В. Евневич // Лазерные технологии в медицине : сб. науч. работ. Челябинск, 2001. С. 166-172.