Численное моделирование динамических процессов в глазу при лазерной экстракции катаракты Текст научной статьи по специальности «Математика»

УДК 617.741-004.1, 519.635.8

М. С. Богомолова, И. Б. Петров

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

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

The article presents a numeric approach to the problem of eye biomedium reaction during the procedure of surgeon cataract extraction by a low-powered laser radiation.

Введение

37

В настоящее время одним из самых актуальных и перспективных направлений в научных исследованиях является численное изучение физиологических, патологических и биологических процессов, происходящих в организме человека. На сегодняшний день медицина представляет собой экспериментальную науку, которая на основе констатации фактов клинических операций может рекомендовать различные средства для снижения патологических процессов. Поэтому построение математических моделей — как патологического, так и нормального функционирования тех или иных органов — необходимо для изучения процессов, протекающих в организме, также важное значение имеет моделирование последствий хирургических операций. Моделирование таких процессов помогает лучше изучить последствия хирургических вмешательств и минимизировать возможные осложнения. Представление о строении биосред и физиологии процессов строится в основном по книгам [6; 7].

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

Хрусталик — это прозрачная линза, расположенная непосредственно за зрачком. Заболевание, при котором возникает помутнение хрусталика, называется катарактой. Потеря прозрачности хрусталика препятствует проникновению лучей света внутрь глаза, зрение человека ухудшается. Если ничего не

предпринимать, катаракта может привести к и прозрачный хрусталик слепоте.

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

хрусталик хрусталик

Рис. 1. Мутный

Вестник РГУ им. И. Канта. 2007. Вып. 10. Физико-математические науки. С. 37—43.

38

воздействие, диабет, некоторые глазные заболевания или наследственные особенности организма человека. Чаще всего катаракта появляется в процессе старения организма, который приводит к увеличению плотности хрусталика и его помутнению. Единственным способом лечения катаракты является хирургическое удаление мутного хрусталика.

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

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

Для проведения операции в глаз вводится иголка с расположенным в ней оптоволокном. При диаметре прокола менее 3 мм можно обойтись без наложения швов. Но даже в этом случае могут возникнуть побочные явления и послеоперационные осложнения. Во время поглощения лазерного излучения выделяется тепло, воздействующее на роговицу. Кроме того, происходит быстрое расширение материала хрусталика, имеющее характер микроразрыва. Это, в свою очередь, приводит к появлению возмущений в биосреде глаза, распространяющихся к его границам, в частности к сетчатке, что делает ее зоной риска.

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

Рис. 4. Ввод гибкой линзы в капсулу хрусталика

Рис. 5. Установленная линза

Модель глаза

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

тело — прозрачный гель, состоящий из внеклеточной жидкости с коллагеном и гиалуроновой кислотой в коллоидном растворе. В задней части глаза его внутренняя поверхность выстлана сетчаткой. Промежуток между сетчаткой и плотной склерой, окружающей глазное яблоко, заполнен сетью кровеносных сосудов — сосудистой оболочкой. У заднего полюса глаза человека в сетчатке есть небольшое углубление — центральная ямка — место, где острота зрения при дневном освещении максимальна. Диаметр глаза — 2,5 см, толщина хрусталика — 4 мм, расстояние от роговицы до хрусталика — 6 мм. Поведение элементов глазного вещества моделировалось с использованием линейно-упругой либо вязко-упругой модели с выделением в ней четырех различных материалов: склеры, хрусталика, стекловидного тела и роговицы. На границах поставлены условия отсутствия нормальных напряжений.

Математическая модель деформируемой биосреды

Под углом 60° к оптической оси глаза вводится световод, по которому в хрусталик передается лазерное излучение. Интенсивность поглощаемого излучения рассчитывается по закону Ламберта — Бера:

Q(r, Н) = 1(г)е-н/Н, 1(г) = 10е-г 1 к)2 ,

где г — радиус облучаемой зоны, Н — вертикальная координата, Н — глубина поглощения (Н = 3,2 мм), 10 — интенсивность в центре световода, Я — диаметр световода (Я = 0,3 мм).

Полная энергия излучения считалась равной 250 мДж.

Р1& I - V ] Ст у = (1)

^ у = кк Q , ; (2)

рС т0

40

© = —8кк + ^ + -£, г', ] = 1,2. (3)

РС Цт0 РС

Уравнения (1) — (3) — уравнения динамики деформируемой среды, используются для математического описания поведения биосреды глаза, уравнение (1) — уравнение движения, уравнение (2) — реологические соотношения, уравнение (3) — уравнение энергии. Здесь

© = Т - Т0 — разница между текущей и начальной температурой, р — плотность среды, о у, 8] — компоненты тензоров напряжений и деформаций, V г, г = 1, 2 — компоненты вектора скорости смещения, с — удельная теплоемкость, у = 3Ка /, К = Ц + 2ц / 3 — коэффициент всеобщего сжатия, Ц, ц — коэффициенты Ламэ, а / — коэффициент линейного расширения, Vу, ] = 1, 2 — ковариантная производная по у-й координате, Бу = оу - 8у окк /3 — компоненты девиатора тензора деформаций, 8 у — символ Кронекера, т0 — время релаксации, Q = Q0 + дгу(9^адТ), Q0 — плотность объемных источников энергии, 0 — коэффициент теПЛопровоДности, у = Ц8у8и - У°к/8у /(РС) + Ц(8й8у/ +8]к8Й ) .

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

Р^г - VгР = 0, Р = (Р-1 + Р2У2 /(РС®))8кк + У/Р/(Рсу)Ч,

где р = -окк/3 — гидродинамическое давление, у / = р -1(др / дТ)в,

Су — теплоемкость в изохорном процессе, 8 кк = ду к / дхк .

Для твердой фазы использовалось уравнение состояния вида Р = Р0ехр(р-у©/К), где р0 — начальная плотность, © = (-у/рТгкк + ^)/(рсв).

Использовалось матричное представление приведенных уравнений

дм . дм . дм с .

(1) “^7 + --= / , где м = (у 1, у 2 , 011,012 , 0 22,0 33, ©) — вектор-

дt дх 1 дх 2

столбец искомых функций, / — вектор-столбец правых частей, а А1, А2 — квадратные матрицы размера 7 * 7. Для ортогональной системы координат вид этих матриц приведен, например, в [8].

Метод численного решения

При аппроксимации системы уравнений явной схемой сеточно-характеристического метода расчетными формулами во внутренних узловых точках (п1т,П2/), т = 1, 2, ..., М — 1, / = 1, 2, ..., I — 1, будут соотношения

м^*!1 = мпт/ +т/т/ + ЬПт/ + ЬПт/,

где

ьПт/=а1[(о-1л^1)т/ (мт-1,/ - мт/) - (о-^ол (мт+1,/ - мт)], ЪПтГ =02[(а21Л+2^ 2)т/ (мт,/-1 - мI/) - (О-^ ^ (м"т/+1 - мт/)], Л±у = (Л у ±Л у |)/2, о у = Н., Л у = {Ц}, Л у = {| Ц |}, г = 1 * 7, у = 1,2.

Здесь ХУ — собственные значения матриц Л:, ] = 1, 2, определяемые

из

характеристических уравнений деЬ(Лу -Х:Е) = 0, Оу = {ю^} — неособенные матрицы, строками которых являются линейно независимые левые собственные векторы юі матриц Лу, определяемые с точностью до их длины из совокупности линейных однородных систем уравнений (ЛТ -Х^юі = 0, і = 1, ..., 7, О-1 — обратные к Оу матрицы, ЛТ — транспонированные матрицы Л у. Аналитическое выражение элементов матриц О-1 аналогично приведенному в [8].

При построении расчетных формул на границах прямоугольной (в координатах (£, п1, п 2)) области интегрирования рассмотрим только верхнюю (п2 = 1) и нижнюю (п2 = 0) границы, имея в вицу, что остальные границы (п1 = 0, Пг) часто являются плоскостью (или осью) симметрии, либо выбираются таким образом, чтобы за рассматриваемое время І ?; I возмущения от неоднородностей в начальных данных не достигали этих границ. Умножая полученные разностные соотношения скалярно на собственные векторы (ю2 , получаем соотношения

(ю2гт1ипт+1 = в2 = (ю2)пт1 «г +т/тг + ь п1т1) ±-^ (х2) пт1 (ю2)пт1 - <г),

Н 2

I = 0, 1, ..., I, аппроксимирующие с первым порядком точности условия совместности вдоль линий пересечения характеристических поверхностей системні и координатной плоскости п1 = П1т (с уравнениями

й-ц 2 = Х2й£): ю2 иг + Х2 ю2 и п = ю2(/ - Л1 п1), і = 1, 2,...,7.

Как известно, число граничных условий для гиперболической системы уравнений определяется числом отрицательных (положительных) собственных значений матрицы А2 на верхней (соответственно, нижней) границе области интегрирования. В рассматриваемых ниже задачах на верхней границе п 2 = 1 имеет место Х^ < Х^ < 0, на нижней границе п 2 = 0, соответственно, Х2 > Х2, > 0 и, следовательно, на каждой из этих границ требуется постановка двух граничных условий. Запишем их в виде:

Ф і (£, п1, и1,..., и7) = 0, і = 1, 2 при п2 = 0,

Фі (Ї, п 1, и1,..., и7) = 0, і = 6, 7 при п2 = 1, причем необходимо, чтобы деі О- ф 0, деі О+ ф 0, где

II_ ___ 2 21Т

^^— = ю1 ю2 ю з ... ю^ , ^^+ =

2 2 — —

ю -\ ... ю с ю ^ ю о

Здесь юі = {дфі /ди1,..., дфі /ди7}, і = 1, 2, 6, 7, а ю2, і = 1, 2, ..., 7 — собственные векторы матрицы А2. Для наших задач граничные условия выбирались полулинейными и после их аппроксимации имели вид:

Фі = юі (~П+1, п 1т )<+1 - $і (~^ П1т ), і = 1 2, П2 = 0, и І = 6, ^, П2 = 1.

41

Т

42

Привлекая полученные разностные соотношения при г = 3, 7 для

границы п2 = 0 и при г = 1, ...,5 — для границы п2 = 1, получаем все необходимые расчетные формулы для точек, принадлежащих этим границам. Для расчета точек на границах п1 = 0, п1 = П* использовались обычные формулы с привлечением дополнительных «лучей» П1 = -к 1, (т = —1), п1 = П* + ^1, (т = М + 1), для которых компоненты искомого вектора и определялись по данным внутри области интегрирования с учетом соответствующей симметрии или периодичности решения либо экстраполяцией (в зависимости от типа границы). Вид матриц О у и собственных векторов матриц Ау дан в [8]; там же приведены формулы для расчета функций на границах области.

Шаг по времени выбирался из условия устойчивости Куранта — Фридриха — Леви, имеющего вид

f 2

Т <

-1

^ hi 1 max/

V i=l

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

Результаты расчетов

Параметры материалов глаза.

Склера: E = 2,1 • 1011 дин| см2, v = 0,33, p = 1 г| см3.

Хрусталик: E = 1,4 • 1011 дин|см2, v = 0,2, p = 1 г|см3.

Стекловидное тело: E = 1011 дин|см2, v= 0,45, p = 1 г|см3.

Роговица: E = 1011 дин| см2, v = 0,33, p = 1 г| см3.

Лазерное воздействие моделируется путем задания в правых частях уравнений динамики деформируемой среды источника тепла Q(t):

Q(t) =

є

V • Tl

є

V • Tl

, Ton < t < (To + T1)n,

P

•-0,(T, + Tl)n < t < T,(n +1), Pl

где є = 106 эрг., Т1 = 3 10—6 с, Т0 = 15 10—6 с, Р0/Р1 = 10—2, а = 5,2 105 К—1,

V = 10—3 см3 (є — энергия импульса, Ті — длительность импульса, Т0 — период импульса, а — коэффициент теплового расширения,

V — объем, в котором действует излучение в начальный момент времени, п — количество импульсов, Р0/Р1 — характеризует уровень шума.

На рисунке 7—10 приведены поля скоростей в материалах глаза. Видно, что в течение некоторого времени векторы скоростей направлены из глаза (рис. 7, 8), затем внутрь глаза (рис. 9), после взаимодействия с его границами. На рисунке 10 показаны зоны возможных разрушений в глазу (с отрицательными значениями давления).

Рис. 7. Поле скоростей,

Рис. 8. Поле скоростей, і = 1,9-106с

і = 5-10"7е

/

\

43

Рис. 9. Поле скоростей,

Рис. 10. Зоны возможных

і = 6-10"6с

повреждений

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

1. Балановский Н. Н., Бубнов А. В., Обухов А. С., Петров И. Б. Расчет динамических процессов в глазу при лазерной экстракции катаракты // Мат. моделир. 2003. Т. 15. № 11. С. 37-44.

2. Петров И. Б. О численном моделировании биомеханических процессов в медицинской практике // Информ. технол. и вычисл. системы. 2003. Т. 1 — 2. С. 102 — 111.

3. Федоров С. Н., Егорова Э. В., Холодов А. С., Марченкова Т. Е., Бубнов А. В. О численном моделировании процессов ирригации и аспирации при экстракапсу-лярной экстракции катаракты // Вопросы кибернетики. Биомединформатика и её приложения. 1988. С. 99—114.

4. Марченкова Т. Е., Бубнов А. В., Холодов А. С. Математическое моделирование ирригационно-аспирационной техники факоэмульсификации // Офтальмохирургия. 1991. № 1. С. 11 — 15.

5. Самарский А. А., Гулин А. В. Численные методы математической физики. М.: Научный мир. 2003.

6. Физиология человека. М.: Мир, 1996. Т. 1 — 3.

7. Краткая медицинская энциклопедия. М.: Сов. энциклопедия, 1989. Т. 1.

8. Петров И. Б., Холодов А. С. Численное решение некоторых динамических задач механики деформируемого твердого тела сеточно-характеристическим методом // Журн. вычисл. мат. и мат. физ. 1984. Т. 24. № 5. С. 722 — 739.

М. С. Богомолова — ст. препод., РГУ им. И. Канта, maria_ishanova@mail.ru. И. Б. Петров — д-р физ.-мат. наук, проф., МФТИ.

Об авторах