МОДЕЛИРОВАНИЕ БИОМЕХАНИЧЕСКИХ ПРОЦЕССОВ

 

Петров И.Б.

 

(Москва)

 

Представлено численное решение следующих задач биомеханики:

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

      исследование волновых явлений, происходящих при черепно-мозговых травмах с целью определения очагов поражения мозга;

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

 

NUMERICAL MODELING OF BIOMECHANICAL PROCESSES

 

Petrov I.B.

 

(Moscow)

 

Problem of biomedium behavior during the eye surgeon operation and cranial traumas are solving by numerical method. The behavior of such a solid medium was described by dynamic system of viscoelastic equations (system of partial differential equations). For the numerical solution of the equations grid-characteristic method was used.

 

1. Введение

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

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

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

2.Разработка численного метода, адаптированного к рассмотрению интересующих нас процессов (в данном случае - динамических).

3.Разработка численного алгоритма, реализующего метод для решения конкурентной задачи.

4.Реализация расчетной программы и систем визуализации расчетных данных.

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

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

2.1. Математическая модель глаза

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

Рис. 1. Схема горизонтального сечения правого глаза

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

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

2.2. Постановка задачи о черепно-мозговых травмах

В начале 80-х годов в работах, посвященных биомеханике черепно-мозговых травм, получила распространение гипотеза о противоударном механизме ушибов мозга при закрытых черепно-мозговых травмах. Гипотеза призвана объяснить появление видимых очагов деструкции, кровоизлияний и некрозов, наблюдаемых средствами компьютерной томографии, не только в месте приложения механической силы (область удара), но и на противоположной стороне мозга (области противоудара). Однако единой теории противоударного механизма, которая бы удовлетворительно объясняла большинство клинических случаев, на сегодняшний день нет. Основные продвижения в данной области связаны с механической и кавитационной теориями. Согласно первой, главным фактором, определяющим степень выраженности и распространенности, а также динамику развития патологического очага в области противоудара, являются 'первичные' повреждения, т.е. повреждения, обусловленные непосредственно механическим воздействием в момент удара; вторая отводит определенную роль мозгового кровообращения, вызывающим 'вторичные' повреждения в области противоудара, появление которых может растягиваться во времени на периоды 2-10 суток, а может происходить и в течение минут или часов, особенно если имелось повреждение многих сосудов вплоть до полного разрушения их стенки.

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

Кавитационная теория строится на предположении, что вследствие ротационного движения мозга или смещения по инерции его массы на противоположной стороне образуется вакуум. Отрицательное давление, действующее в течение 1/700 - 1/1000 с, вызывает в первую очередь в текущей жидкости (крови) появление пузырьков газа (кавитацию). Перемещаясь с потоком крови в область более высокого давления, кавитационный пузырек 'схлопывается', вызывая при этом гидродинамическую ударную волну, которая и разрушает стенку сосуда в той или иной степени. Кавитационная теория хорошо объясняет развитие во времени очагов поражения в области противоудара, отсутствие в них (в первые часы и сутки) крови и травматически поврежденных клеток. Однако ясного представления, насколько кавитационные эффекты действительно имеют место в крови артерий мозга, и каков масштаб этих эффектов, нет. За пределами теории остается также то, что локализация отрицательных давлений (положительных напряжений) происходит в основном в лобной и височной областях. При затылочных травмах противоударные очаги развиваются существенно сильнее, чем при лобовых.

а) Сагитальная проекция

б) Трансверсальная проекция

Рис. 2. Структура расчетных областей с выделенными материалами (по данным [5], перераб.): 1- кость, 2 - мозговое вещество

На коротких временах (порядка 1/1000 с) существенную роль играют ударно-волновые и акустические эффекты. Имеет смысл предложить, что именно сложная волновая картина, обусловленная как геометрическими особенностями анатомии черепа, так и различиями в упругих свойствах пластинчатых и губчатых костей черепа, определяет пространственное распределение областей положительных напряжений. Получение расчетным путем его качественной картины и сравнение с данными.

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

2.3. Строение и общая физиология сердца (сведения необходимые для построения математической модели)

Кровь может выполнять свои разнообразные функции, только находясь в постоянном движении. Это движение крови обеспечивается сердцем. Сердце можно рассматривать как два полых мышечных органа - 'левое' сердце и 'правое' сердце, каждое из которых состоит из предсердия и желудочка. Лишенная кислорода кровь от органов и тканей организма поступает к правому сердцу, выталкивающему ее к легким. В легких кровь насыщается кислородом, возвращается к левому сердцу и вновь поступает к органам. Таким образом, правое сердце перекачивает дезоксигенированную кровь, а левое - оксигенированную. Движение крови по сосудам легких от правого сердца к левому называется легочным кровообращением (малый круг). Кровоснабжение всех остальных органов (и отток крови от них) носит название системного кровообращения (большой круг). Разумеется, фактически оба этих отдела составляют единое кровеносное русло, в двух участках которого (правом и левом сердце) крови сообщается кинетическая энергия. Нагнетательная функция сердца основана на чередовании расслабления (диастолы) и сокращения (систолы) желудочков. Во время диастолы желудочки заполняются кровью, а во время систолы они выбрасывают ее в крупные артерии (аорту и легочный ствол). У выхода из желудочков расположены клапаны, препятствующие обратному поступлению крови из артерий в сердце. Перед тем как заполнить желудочки, кровь притекает по крупным венам (полым и легочным) в предсердия. Систола предсердий предшествует систоле желудочков; таким образом, предсердия служат как бы вспомогательными насосами, способствующими заполнению желудочков.

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

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

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

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

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

, , , i, j = 1, 2,

где - компоненты вектора скорости, U - плотность внутренней энергии, - компоненты тензоров деформаций и напряжений, - ковариантная производная по k-й координате, Q - плотность источников энергии, тензор четвертого ранга задается в соответствии с реологией среды. Для модели упруго-идеальнопластической среды

,

где - символ Кронекера, λ, μ - постоянные Ламе, k - предел текучести на единичный вектор, ρ - плотность, - девиатор тензора напряжения, - коэффициент линейного расширения, I определяется из условия текучести Мизеса:

Мы использовали матричное представление приведенных уравнений [2].

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

При использовании для аппроксимации системы уравнений явной схемы сеточно-характеристического метода, расчетными формулами во внутренних узловых точках , m = 1, 2,:, М - 1, l = 1, 2,:,L - 1, будут соотношения

, , , , j = 1, 2, i = 1, 2,:,7,

- собственные значения матриц k = 1, 2, определяемые из характеристических уравнений - неособенные матрицы, строками которых являются линейно независимые левые собственные векторы матриц , определяемые с точностью до их длины из совокупности линейных однородных систем уравнений , i = 1,:,7; - обратные к , матрицы; - транспонированные матрицы . Аналитическое выражение элементов матриц аналогично приведенному в [2].

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

Умножая полученные разностные соотношения скалярно на собственные векторы , получаем соотношения

l = 0,:, L,

аппроксимирующие с первым порядком точности условия совместности вдоль линий пересечения характеристических поверхностей системы и координатной плоскости (с уравнениями ):

, i = 1, 2,:,7.

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

, i = 1, 2 при ,

, i = 6, 7 при ,

причем необходимо, чтобы , , где, соответственно,

, .

Здесь , i = 1, 2,:,7, а , i = 1, 2,:,7 - собственные векторы матрицы . Для рассматриваемых ниже задач граничные условия выбирались полулинейными и после их аппроксимации имели вид:

, i = 1, 2 при , i = 6, 7 при .

Привлекая полученные разностные соотношения при i = 3,:, 7 для границы и при i = 3,:, 7 - для границы , получаем все необходимые расчетные формулы для точек, принадлежащих этим границам. Для расчета точек, принадлежащих границам , использовались обычные расчетные формулы с привлечением дополнительных 'лучей' (m = -1), (m = M+1), для которых компоненты искомого вектора u определялись по данным внутри области интегрирования с учетом соответствующей симметрии или периодичности решения, либо экстраполяцией (в зависимости от типа границы).

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

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

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

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

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

При a = 1 получим схему первого порядка аппроксимации, при a = 0 - второго; аналогично, положив 0 < a = const < l, получим некую промежуточную схему (также немонотонную, формально первого порядка аппроксимации), которая при хорошем экспериментальном выборе a позволит улучшить качество численного решения вблизи разрывов. Этот подход описан в работах [8, 9].

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

5.1. Моделирование офтальмологической операции

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

Склера:, ,

Хрусталик:,,

Стекловидное тело: , , Роговица:, ,

Рис. 4. Поле ско-ростей, t=5*10-7c

Рис. 5. Поле ско-ростей, t=1.9*10-6с

Рис. 6. Поле ско-ростей, t=6*10-6c

Рис. 7. Зоны воз-можных повреж-дений

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

где  = 10⁶ эрг, Т₁ = 3ћ10^-6 с, Т₀ = 15ћ10^-6 с, P₀/P₁ = 10^-2, α = 5,2ћ10⁵ К^-1, V = 10^-3 см³.

(ε - энергия импульса, Т₁ - длительность импульса, Т₀ - период импульса, α - коэффициент теплового расширения, V - объем в котором действует излучение в начальный момент времени, n - количество импульсов, Р₀/Р₁ - характеризует уровень шума.

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

5.2. Моделирование деформации сердечной мышцы

В расчетах о деформации сердечной мышцы задавались изменение внешних деформаций во времени, упругие константы полагались равными: E = 1,04 10⁹ Па (модуль Юнга), γ = 0,41 - коэффициент Пуассона, 10³ кг/м³ (плотность). На границах области интегрирования задавался вектор скорости, в местах выхода крови - нормальные условия. В начальный момент времени построенная квазиортогональная расчетная сетка показана на рис. 8. Поля скоростей в разные моменты времени (при сжатии сердечной мышцы): t₁ = 2,53 10^-5 c, t₂ = 5,07 10^-5 c, t₃ = 5,3ћ10^-3c приведены на рис. 9-11. Хорошо виден процесс выталкивания крови из желудочков.

5.2. Моделирование последствий черепно-мозговых травм.

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

6,5ћ10⁹ н/м²; 6,67ћ10⁵ н/м²; коэффициенты Пуассона: 0,22; 0,44; плотности: 1,51 кг/м³, 1,02ћ10³ кг/м³, соответственно. На рис. 12, 14 изображены расчетные сетки для численного решения рассматриваемой задачи в двух проекциях: вертикальной и горизонтальной. Соответственно, ход волновых процессов изображен на рис. 13, 15. Видно, что возмущения распространяются значительно быстрее в костной ткани, что может привести к появлению областей поражения в местах, не прилегающих к зоне удара.

Рис. 12.

Рис. 13.

Рис. 14.

Рис. 15.

 

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

1.Работнов Ю.Н. Механика деформируемого твердого тела. - М.: Наука, 1988.

2.Петров И.Б., Холодов А.С. Численное исследование некоторых динамических задач механики деформируемого твердого тела сеточно-характеристическим методом. ЖВMиМФ. 1984. - Т. 24. - ? 5. - С. 722-739.

3.Магомедов К.М., Холодов А.С. Сеточно-характеристические численные методы. - М.: Наука, 1988.

4.Г. Тевс, Р. Шмидс. Физиология человека. М.: Мир, 1996. - Т. 1-3.

5.Петровский Б.В. Краткая медицинская энциклопедия. М.: Советская энциклопедия, 1989. - Т. 1-3.

6.Громов А.П. Биомеханика травмы. - М.: Наука, 1979.

7.Лебедев В.В., Крылов В.В. Замечания к патогенезу ушибов мозга, возникающих по противоударному механизму, в остром периоде их развития. // Нейрохирургия, 1998. - С. 22-25.

8.Петров И.Б., Холодов А.С. О регуляризации разрывных численных решений уравнений гиперболического типа. // ЖВMиМФ. 1984. - Т. 24. - ? 8. - С. 1172-1188.

9.Петров И.Б., Тормасов А.Г., Холодов А.С. Об использовании гибридизированных сеточно-характеристических схем для численного решения трехмерных задач динамики диформируемого твердого тела. // ЖВMиМФ. 1990 - Т. 30. - ? 8. - С. 1237-1244.

10.Иваненко С.А. Расчет течений в водоемах на криволинейных сетках. / препринт ВЦ АН СССР. М., 1991.

11.Бубнов А.В., Коротин П.Н., Холодова О.А. Математическое моделирование ультразвуковых колебаний факоэмульсификатора в среде хрусталика./ В сб. Информатика и медицина: М.: Наука, 1997, С. 144--155.

12.Белоцерковский О.М. Численное моделирование в механике сплошных сред. - М.: Физматлит, 1994.