CC BY
82
УДК 519.6
А. В. Васюков, И. Б. Петров
МОДЕЛИРОВАНИЕ МЕХАНИЧЕСКИХ ФАКТОРОВ ЧЕРЕПНО-МОЗГОВЫХ ТРАВМ СЕТОЧНО-ХАРАКТЕРИСТИЧЕСКИМ ЧИСЛЕННЫМ МЕТОДОМ
Приведены некоторые результаты численного исследования с применением сеточнохарактеристических методов на неструктурированных сетках механической реакции системы «череп - мозг человека» на ударные воздействия. Сформулированы механико-математические модели покровов мозга. Получено пространственное распределение нагрузок для различных моделей.
Some results of numerical modelling of mechanical factors of cranial-cerebral injury caused be percussive stress are described. The usage of grid-characteristic method on unstructured grid is discussed. The authors introduce some mechanical-mathematical models of brain teguments and provide results obtained using different models.
Ключевые слова: математические модели, черепно-мозговая травма, сеточно-характеристические численные методы, неструктурированные сетки.
Key words: numerical models, cranial-cerebral injury, grid-characteristic method, unstructured grid.
Введение
Численное изучение физиологических и патологических процессов, происходящих в организме человека, позволяет получить новые качественные и количественные характеристики функционирования органов в различных условиях и протекающих в них нормальных и патологических процессов, что необходимо для прогнозирования их развития, предсказания последствий патологий, выдачи медицинских рекомендаций и разработки новых принципов диагностики на ранних стадиях различных заболеваний. В частности, сегодня медицина черепномозговой травмы — исключительно экспериментальная область, способная произвести констатацию фактов и рекомендовать операционные или медикаментозные средства для ослабления патологических процессов, но проблема построения моделей для понимания механизмов протекающих патологических процессов открыта.
Из нейрохирургической практики известно, что области поражения мозга не всегда совпадают с областями, прилежащими к месту удара. Пример — известный феномен «противоудара»: при ударе затылком область поражения мозга локализуется в лобной части головы человека. Объяснение явлениям, наблюдаемым при черепно-мозговой травме, может дать только изучение сложных волновых картин, образующихся в неоднородной биосреде (или биоконструкции), которую представляет из себя голова (рис. 1). Слоистая конструкция ослабляет действие продольных (пористый слой черепной коробки) и поперечных (слой ликвора под костной частью) упругих волн, напряжений, вызванных нормальным (твердая костная часть черепа), а также скользящим (волосяной покров и кожа) ударами.
а б
Вестник Российского государственного университета им. И. Канта. 2010. Вып. 10. С. 42 — 51.
Рис. 1. Срезы черепной коробки: а — кости и твердая оболочка, сагиттальное сечение; б — оболочки мозга (фронтальный срез):
1 — кожа; 2 — надкостница; 3 — кость черепа; 4 — продольный шов; 5 — твердая оболочка; 6 — паутинная оболочка; 7 —
сосудистая оболочка; 8 — субарахноидальное пространство;
9 — венозная впадина; 10 — арахноидальные грануляции; 11 — серп большого мозга;
12 — нижний сагиттальный синус (расстояния между оболочками преувеличены)
Механико-математические модели реакции головы на ударные нагрузки позволяют исследовать качественные зависимости влияния индивидуальных особенностей строения головы человека на степень риска при черепно-мозговых травмах и предсказывать возможные физиологические патологии, которые могут появиться после черепно-мозговой травмы, поскольку карты мозга (области, отвечающие за те или иные функции организма человека) хорошо известны нейрофизиологам.
1. Уравнения механики деформируемого твердого тела
Для описания поведения биологических тканей как сплошной среды под воздействием механических нагрузок ударного типа использовалась система уравнений линейной теории упругости [2], уравнения движения и реологические соотношения:
Р Ч = Vі , ст ї = еы + Ру, (1)
где р — плотность среды; VI — компоненты скорости смещения;
Сту, е^ — компоненты тензоров напряжений и скоростей деформаций;
V — ковариантная производная по ]-й координате;
Ру — правая часть.
Тензор 4-го порядка Цщ определяет реологию среды. В случае линейноупругого тела его компоненты выражаются через две независимые постоянные — константы Ламе X и ц:
Цщ = Х8І) 8к/ + ц(8ік ). (2)
Также использовалась модель вязкоупругого тела Максвелла.
Плотность определяется из уравнения состояния р = р0е(р/К), где р = -3^^кк — давление;
К = Х +2 ц — коэффициент всестороннего сжатия. Уравнения (1; 2) допускают запись в матричной форме:
Iй+А +А ^ = I (3)
где й = (VI, vг, сти, СТ12, а22, Ъзз)т — вектор искомых функций; / — вектор правых частей той же размерности; Аі — матрицы 6 х 6, явный вид которых приведен в работе [3]; х1, х2 — независимые пространственные переменные; £ — время. Если матрицы Аі в уравнении (3) имеют 6 вещественных собственных чисел, то такая система называется гиперболической, ее решения соответствуют процессам, которые обычно называются волновыми — распространение возмущений вдоль характеристических конусов в пространстве (х1, х2, і).
2. Механико-математическая модель системы «череп — мозг»
Нами было сформулировано несколько механико-математических моделей головы человека. Реологические свойства биоматериалов подвергались вариации. Так, реология мозгового вещества полагалась как линейноупругой, так и вязкоупругой. Поведение костного материала моделировалось как изотропной линейноупругой сплошной средой со средними свойствами пластинчатой и губчатой кости, так и многокомпонентными моделями, явно выделяющими костную ткань с разными механическими характеристиками.
Моделирование взаимодействия между черепом и мозгом —сложная задача, так как в действительности мозг имеет большое количество различных по механическим свойствам оболочек, складчатых структур, врастающих друг в друга, с полостями, заполненными жидкостью (ликвором). В данной работе применялся метод явного выделения контактного разрыва с контактными условиями, которые варьировались от полного слипания до скольжения с возможностью отслоения.
3. Сеточно-характеристические методы для систем одномерных гиперболических уравнений
Для численного исследования моделей был использован ряд известных конечно-разностных методов, относящихся к классу сеточно-характеристических. Для системы уравнений с одной пространственной переменной
Iй+А аХй = 1 (4)
решение ищется в виде сеточной функции йп, заданной в узлах сетки {хШ = шк, Ґ1 = ит), где к и т — шаги по пространству и по времени.
Монотонная схема. Эта схема строится на основе анализа поведения характеристик системы уравнений (4) и приводит к следующим формулам [3]:
“Ш+1 = йШ-ъО-1Л+О (йШ+1 -йШ)-аО-1Л-О (йШ -йпт_1), (5)
где О — матрица, строки которой являются левыми собственными векторами матрицы А; Л± —
диагональные матрицы, содержащие соответствующие собственные числа. Схема (5) имеет порядок аппроксимации 0(к, т), обладает свойством монотонности и минимальной аппроксима-ционной вязкостью среди монотонных схем первого порядка, что важно при расчете динамических процессов в неоднородных средах.
Схема 2-го порядка точности. Единственная центральная схема
2-го порядка аппроксимации на трехточечном шаблоне (то есть использующая значения в точках Ш - 1, Ш, Ш + 1)) — схема Лакса — Вендроффа:
С1 = йШ -ъА )г(йШ+1- йШ-1)+°2 а2 У2 (йШ+1- 2йШ + йШ-1). (6)
Схема (6) обладает минимальным размазыванием волнового фронта, но не является
монотонной. Это имеет вид нефизичных осцилляций вблизи разрывов точного решения.
Гибридные разностные схемы. Если составить линейную комбинацию двух предыдущих схем, то гибридную схему можно записать в виде
Ш = йШ -ъА 12 (йШ+1- йШш-1)+((1 - «)ъ°-1М О+аъ2 а2 ) 12 (йШ+1- 2йШ + йШш-1). (7)
Здесь |л| обозначает диагональную матрицу, составленную из модулей собственных значений
матрицы А. При а = 0 мы получаем схему первого порядка (5), при а = 1 — схему второго порядка (6). Схема (7) называется гибридной, если коэффициент а выбирается в соответствии с локальными свойствами решения [4], и гибридизированной, если он имеет фиксированное значение [5], подбираемое экспериментально. В данной работе локальная гладкость решения определялась из условия, предложенного Р. П. Федоренко [9]:
Уі(Йт+1 ~-Йт +Йт-і) ^ ^ Уг (Йт+1 ~ Йт-1 )• (Ю
В расчетах полагалось К = 0,5. В случае гладких решений применяется схема второго порядка (6), если решение имеет разрывный характер, то используется схема первого порядка (5).
4. Разностные схемы для двумерных гиперболических уравнений
В случае двумерных гиперболических уравнений разностные схемы на регулярной расчетной сетке {х/ = /к, уШ = Шк, їи = ит) могут быть получены из перечисленных одномерных разностных схем следующим образом. Каждую явную схему на трехточечном шаблоне можно представить в виде действия оператора
К1 = иШ + С (А, т к) {иЩ-1, иШ, иШ+1}. (9)
Тогда соответствующая двумерная схема на пятиточечном шаблоне будет иметь вид
1^/,. лг и и ^ 1
иП1 = и: +
2 С(2Ах, Т к {ш , иІШ, и”+1,Ш}+ } С(2Ау , Т к) иИШ , иИ,Ш+1} . (10)
Идея другого метода — расщепление по пространственным направлениям — заключается в замене двумерной системы (3) тремя одномерными системами
5 5 5 5 5
и + А1 — и = 0, - и + А 2 —и = 0, ^ и = ^ (11) — (13)
дЬ дх1 дЬ дх2 дЬ
Решение их производится с помощью некоторой одномерной разностной схемы поэтапно: сначала во всей области интегрирования решается система (11), затем результат используется для решения системы (12) и, наконец, конечное решение получается в результате прибавления правой части (13).
5. Расчет контактных границ
Расчеты, в которых область интегрирования распадается на участки различных реологий и механических свойств, возможно, перемещающиеся относительно друг друга, требуют аккуратного решения задачи контактного разрыва на поверхностях соприкосновения тел. Сеточно-характеристические схемы позволяют формулировать условия на контактной границе двух областей интегрирования в явном виде. Как показывают результаты работ [7] и [8], явное выделение контактных границ имеет ряд преимуществ над их сквозным расчетом.
Полный алгоритм расчета контактных границ для двумерного случая приведен в работах [7; 8]. Данный алгоритм позволяет явно выделять контактные границы между телами разной реологии и рассчитывать их с учетом различных контактных условий.
Использовались условия полного слипания, то есть равенство скоростей, нормальных напряжений и касательных напряжений:
ЭД = 0[2], а„ [1] = а„ [2], ат [1] = ат [2], (14)
а также условие свободного скольжения, то есть равенство нормальных скоростей и напряжений и
отсутствие касательных напряжений:
Vп [1] = Ъп [2] °п [1] = ап [2Ь ат [1] = °т [2] = ° (15)
6. Обобщение на неструктурированные треугольные сетки
Нами был реализован алгоритм, позволяющий применить описанные двумерные разностные схемы для расчетов на неструктурированных треугольных сетках. При этом значения в узлах пятиточечного шаблона определялись линейной интерполяцией внутри прилежащих треугольников.
Отметим, что уже имеются программные инструменты промышленного уровня для построения треугольных сеток в областях произвольной формы и связности. Например, использовалась программа Triangle (автор Шевчук). Треугольная сетка позволяет управлять размерами ячеек, добиваясь либо полной однородности, либо при необходимости сгущения или разряжения в определенных зонах. Таким образом, можно сделать вывод, что использование гибридных схем с расщеплением на неструктурированных треугольных сетках вполне оправдано, так как позволяет получать численные решения удовлетворительного качества в областях интегрирования произвольной формы. Примеры использованных расчетных сеток приведены на рисунке 2.
а
в
Рис. 2. Расчетные сетки для моделей головы человека: а — двухкомпонентная модель; б — модель с желудочками; в — модель с желудочками и мембраной
7. Результаты расчетов
По результатам сравнения расчетных данных, полученных ранее, с клиническими данными по 18 пациентам, а также с имеющимися в медицинской литературе качественными описаниями биомеханики черепно-мозговой травмы сделан вывод об удовлетворительной описательной способности модели. Пример сравнения дан на рисунках 3, 4.
Рис. 3. КТ-снимок пациента с ушибом головного мозга тяжелой степени, удар слева, стрелками указаны места поражения мозгового вещества
Рис. 4. Распределение максимальных сдвиговых нагрузок
В предыдущих расчетах также получено, что картина распространения возмущений в стенке черепа носит явно выраженный волновой характер, наблюдается многократное отражение волн от контактных границ между отдельными компонентами. Это говорит о том, что внутреннее строение рассматриваемой системы может оказывать сильное влияние на происходящие процессы. В связи с этим было выполнено моделирование с использованием различных моделей с целью сравнения их описательной способности.
Геометрия и реология костей черепа достаточно сложны, поэтому даже для многокомпонентной модели был сделан ряд упрощений. Череп был принят состоящим из 3 однородных слоев костной ткани, также в модель были включены внутричерепная жидкость (ликвор) и мозговое вещество (рис. 5).
Рис. 5. Вид расчетной области (размер по вертикали не соблюден):
1 — мозговое вещество; 2 — ликвор; 3 — внутренний слой компактной костной ткани;
4 — слой губчатой костной ткани; 5 - внешний слой компактной костной ткани
Для сравнения использовались две упрощенные модели. В первой из них слои 3 — 5 (рис. 5) были объединены в один слой, состоящий из компактной костной ткани. Во второй модели слои объединялись так же, но итоговый слой состоял из губчатой ткани. Итак, три модели (многокомпонентная и две упрощенные) позволили рассмотреть влияние как сложности модели (числа слоев), так и механических характеристик материала черепа. Моделировался удар шариком из жесткого пластика сантиметрового диаметра, налетающим со скоростью 3 м/ с перпендикулярно стенке черепа. Принципиальная картина волновых процессов не зависит от скорости шарика, ее можно выбрать малой для удобства расчетов.
На правой границе, соответствующей внешней стенке черепа (рис. 5), ставилось условие свободной границы. На трех остальных границах также было условие свободной границы, но при
этом геометрия области выбиралась таким образом, чтобы за рассматриваемое время сильное возмущение не успело дойти до этих границ. Расчет различных слоев костной ткани черепа, внутричерепной жидкости и мозгового вещества производился сквозным образом.
На рисунках 6 и 7 даны изолинии напряжений в различные моменты времени, полученные с использованием разнообразных моделей. Из рисунков видно, что на начальной стадии, пока картина еще носит ярко выраженный волновой характер, распределение напряжений в случае использования разных моделей покровов мозга сильно отличается.
б
в
д
Рис. 6. Многокомпонентная модель. Изолинии напряжений. Ось х направлена вправо; ось у — вверх; шарик налетает справа: а — охх при і = 7,5 мкс; б — охх при і = 12,5 мкс; в — о при і = 7,5 мкс; г — о при і = 12,5 мкс; д — о при і = 7,5 мкс; е — о при і = 12,5 мкс.
а
г
е
а б в г де
Рис. 7. Модель с однородным черепом из компактной костной ткани.
Изолинии напряжений. Ось х направлена вправо; ось у — вверх; шарик налетает справа: а — охх при і = 7,5 мкс; б — охх при і = 12,5 мкс; в — оху при і = 7,5 мкс; г — оху при і = 12,5 мкс; д — оуу при і = 7,5 мкс; е — оуу при і = 12,5 мкс.
Итак, результаты применения различных моделей отличаются друг от друга на начальных стадиях, пока протекающие процессы носят ярко выраженный волновой характер. Можно утверждать, что для получения достоверной картины распределения нагрузок в мозге важно учитывать внутреннее строение его покровов.
Список литературы
1. Белоцерковский О. М. Численное моделирование в механике сплошных сред. М., 1994.
2. Новацкий В. К. Теория упругости. М., 1975.
3. Магомедов К. М., Холодов А. С. Сеточно-характеристические численные методы. М., 1988.
4. Петров И. Б., Холодов А. С. О регуляризации разрывных численных решений уравнений гиперболического типа // Журн. вычислительной математики и математической физики. 1984. Т. 24, № 8. С. 1172—1188.
5. Петров И. Б., Тормасов А. Г., Холодов А. С. Об использовании гибридизированных сеточно-характеристических схем для численного решения трехмерных задач динамики деформируемого твердого тела // Там же. 1990. Т. 30, №8. С. 1237—1244.
6. Агапов П. И., Петров И. Б., Челноков Ф. Б. Численное исследование задач механики деформируемого твердого тела в неоднородных областях интегрирования // Обработка информации и моделирование: сб. ст. М., 2002. С. 148—157.
7. Белоцерковский О. М., Агапов П. И., Петров И. Б. Моделирование последствий черепно-мозговой травмы / / Медицина в зеркале информатики. М., 2008. С. 113—124.
8. Петров И. Б., Тормасов А. Г., Холодов А. С. О численном изучении нестационарных процессов в деформируемых средах многослойной структуры // Изв. АН СССР. Механика твердого тела. 1989. № 4. С. 89—95.
9. Федоренко Р. П. Введение в вычислительную физику. М., 1994.
Об авторах
Алексей Виктрович Васюков — асп., МФТИ, Москва, e-mail: vasyukov@gmail.com
Игорь Борисович Петров — д-р физ.-мат. наук, проф., МФТИ, Москва, e-mail: petrov@mipt.ru
Authors
Aleksey Vasyukov — PhD student, MlPT, Moscow, e-mail: vasyukov@ gmail.com Professor lgor Petrov — MlPT, Moscow, e-mail: petrov@mipt.ru
