автореферат диссертации по информатике, вычислительной технике и управлению, 05.13.18, диссертация на тему:Численное моделирование механических факторов черепно-мозговой травмы

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

Агапов Павел Игореви

Численное моделирование

механических факторов черепно-мозговой травмы

Специальность 05.13.18 Математическое моделирование, численные методы и комплексы программ

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

АВТОРЕФЕРАТ

Москва - 2005

Работа выполнена на кафедре информатики Московского физико-технического института (государственного университета)

Научный руководитель:

доктор физико-математических наук, профессор Петров Игорь Борисович

Официальные оппоненты:

доктор физико-математических наук, профессор Лобанов Алексей Иванович,

кандидат физико-математических наук Колобов Андрей Владимирович

Ведущая организация:

Институт автоматизации проектирования РАН

Защита состоится « 2005 г. в часов на

заседании диссертационного совета К 212.156.02 при Московском физико-техническом институте (государственном университете) по адресу: 141700, Московская обл., г. Долгопрудный, Институтский пер., д. 9.

С диссертацией можно ознакомиться в библиотеке МФТИ.

Автореферат разослан «. » /^СЛ^ьЗ^- 2005 г.

Ученый секретарь диссертационного совета кандидат физико-математических наук — ____Федько О. С.

1 I

Общая характеристика работы

Актуальность темы

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

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

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

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

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

Цель работы

Целью работы является получение качествшнгА^щщтены распределения механических нагрузок на мозг человека т^Ф р^личМыхадд|йр1ть1х

"Р*

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

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

Научная новизна

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

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

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

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

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

Практическая ценность

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

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

Защищаемые положения

На защиту выносятся:

1. Математическая модель механической реакции системы череп-мозг на ударное воздействие для численного расчета пространственного распределения механических нагрузок и деформаций тканей мозга при черепно-мозговой травме.

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

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

4. Обобщение метода перестройки регулярных сеток для адаптации к деформирующимся границам на неструктурированные треугольные сетки.

Публикации

Научные результаты диссертации опубликованы в работах [1-11].

Апробация

Результаты работы были доложены, обсуждены и получили одобрение специалистов на следующих научных конференциях:

• международная конференция «Математика. Компьютер. Образование» (Пущино, 2005);

• научные конференции МФТИ «Современные проблемы фундаментальных и прикладных наук» (Долгопрудный, 2002 - 2004).

• научно-практическая конференция «Комбинированная и сочетан-ная патология: проблемы диагностики и лечения» (Москва, ГВКГ им. H. Н. Бурденко, 2003)

• четвертый Всероссийский научный симпозиум «Механизмы стресса в экстремальных условиях» (Москва, ГосНИИИ военной медицины Минобороны России, 2005)

Структура и объем диссертации

Диссертация состоит из введения, 3 глав и заключения. Общий объем диссертации составляет 88 страниц. Список литературы содержит ссылки на 36 публикаций.

Содержание работы Введение

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

Глава 1. Математическая модель

В главе 1 описана математическая постановка задачи моделирования механики черепно-мозговой травмы.

Двумерные модели черепно-мозгового отдела

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

Во всех случаях внешняя нагрузка моделируется в виде соударения системы череп-мозг, обладающей начальной скоростью =1 — 3м/с, с неподвижной абсолютно жесткой плоскостью.

Рис. 1: Использованные двумерные модели черепно-мозгового отдела: а - двух-компонентная модель, б - модель с желудочками, в — модель с желудочками и

мембраной (falx cerebri)

Уравнения движения

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

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

Здесь р - плотность среды, — компоненты скорости смещения, стч, £у - компоненты тензоров напряжений и деформаций, V} - ковариантная производная по ^'-й координате, Р13 - добавочная правая часть.

Вид компонент тензора 4-го порядка ц%3к1 определяется реологией среды. Для линейно-упругого тела

В этом соотпоптении, которое обобщает закон Гука, Л и ц параметры Ляме, 5г] - символ Кропекера. Плотность определяется из уравнения состояния

р = ро ехр (р/К),

где Р = - 5 <Укк ~ давление, К = А + - коэффициент всестороннего сжатия.

Уравпения (1) допускают запись в матричной форме:

рщ = Vj dij (уравнения движения),

cr'ij = qijki e'ki + Fij (реологические соотношения).

(1)

Qijki = Aijj 5ki + fi(5Zk53i + 5tiS3k).

Здесь и = {г>1,г>2,<ти,¿т12,ст"22!^зз} ~ вектор искомых функций, / - вектор правых частей той же размерности, А, - матрицы 6 х б, х%, Х2 -независимые пространственные переменные, I - время.

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

Глава 2. Численные методы

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

Сеточно-характеристические методы на регулярных сетках

В первой части главы 2 описаны реализованные в программном комплексе конечно-разностные численные методы, сформулированные для гиперболической системы уравнений с одной пространственной переменной: схемы Лакса, Куранта-Изаксона-Риса, Лакса-Вендроффа, а также построенные на их основе гибридные и гибридизированные сеточно-характеристические схемы.

Принципы построения сеточно-характеристических разностных схем, использующих аппроксимацию производных вдоль характеристических кривых уравнения (2), рассмотрены на примере монотонной схемы Куранта-Изаксона-Риса первого порядка:

и^1 = - 0"П-1Л4 П«+1 - О - аП^Л-П« - (3)

Здесь П - матрица, строки которой являются левыми собственными векторами матрицы А, Л* - диагональные матрицы, содержащие соответствующие положительные и отрицательпые собственные числа.

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

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

= - <тА|(гС+т -

Здесь |А| обозначает диагональную матрипу, составленную из модулей собственных значений матрицы А. При а = 0 мы получаем схему первого порядка (3), при а — 1 - схему Лакса-Вендроффа второго порядка. Схема (4) называется гибридной, если коэффициент а выбирается в

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

~ т + ~ ~ "т-1)'

где коэффициент К = 0.5 подбирался экспериментально. В случае выполнения условия гладкости применялась схема Лакса-Вендроффа, иначе схема первого порядка (3).

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

Разностные схемы для неструктурированных треугольных сеток

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

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

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

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

Рис. 2: Реконструкция структурированного шаблона по треугольной сетке

Рис. 3: Результат расчета модельной задачи в круговой области: а - структурированная четырехугольная сетка, б - неструктурированная треугольная сетка

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

Расчет задач контактного разрыва в случае неоднородных тел

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

т

а и + (а ) и = о.

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

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

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

а б

Рис. 4- Примеры решения задач контактного разрыва- а - прохождение плоского импульса через область с неоднородностью; б косое соударение тел с

трением

Перестройка сетки для адаптации к деформирующимся грани-дам области интегрирования

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

В случае регулярной сетки оптимальное расположение узлов минимизирует функционал

1 = (5)

где (х, у) - декартовы координаты, а (£, 7?) - криволинейные координаты вдоль линий сетки. Дискретный аналог функционала (5) имеет вид

= + (6)

где суммирование проводится по всем треугольникам, прилежащим к углам четырехугольной ячейки. Здесь Д123 ~ площадь треугольника, а ¿12^23 ~ Длины его соответствующих сторон.

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

/г ^(ДтГ1^+ & + &). (7)

1 Иваненко С А , Чарахчъян А. А Криволинейные сетки из выпуклых четырехугольников // Ж вычисл матем и матем. физ - 1988. - Т28, N0.4. - С 503-514

Рис. 5' Расчет сверхзвукового соударения: а - вырождение сетки без перестройки, б - сетка в процессе расчета с перестройкой, в - изолинии давления

Каждое слагаемое в (7) представляет собой сумму квадратов сторон, деленную на площадь треугольника. Здесь суммирование производится по всем треугольникам сетки.

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

Глава 3. Результаты моделирования

Глава 3 содержит результаты численного исследования двумерных моделей черепно-мозговой травмы, сформулированных в первой главе.

Интегральная оценка механического воздействия

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

(Jrran — min min аг\ максимальное по времени растягивающее (поло-

te[o,T] «

жительное) напряжение гттих — max тах а%: максимальное по времени

te[o,T] г

касательное напряжепие ттах = тах ¿(тах аг — min <т,); максимальная

te[0,'i'] * г i

Рис. 6: Результаты расчета двухкомпонентной модели: а - максимальное сжатие;

б - максимальное растяжение; в - максимальные сдвиговые напряжения

по времени норма девиатора тензора напряжений Smax = max ^/Si7S%1,

te[o,r¡

также характеризующая сдвиговые напряжения.

Результаты расчетов и сравнение с клиническими данными

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

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

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

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

На рис. 8, а приведен пример КТ-снимка пациента с ушибом головного мозга тяжелой степени, полученным от удара слева при ДТП (данные предоставлены Главным военным госпиталем им. Н. Н. Бурденко). Стрелками на томограмме указаны места поражения мозгового вещества. На рис. 8, б приведен соответствующее распределение максимальных сдвиговых нагрузок.

Рис. 7' Распределение сдвиговых напряжений' а - полное слипание; б - скольжение; в - двухкомпонентная модель

Рис. 8: а - КТ-снимок пациента с ушибом головного мозга тяжелой степени, удар слева; б - соответствующее распределение максимальных сдвиговых наг

грузок

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

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

Заключение

Заключение содержит в развернутом виде формулировку основных результатов и выводов диссертации.

а

б

Основные результаты и выводы диссертации

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

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

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

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

Список публикаций по теме диссертации

1. Агапов П. И. Анализ результатов численного моделирования черепно-мозговой травмы // Процессы и методы обработки информации: Сб. ст. / Моск. физ.-тех. ин-т. - М.: МФТИ, 2005. - С. 186 - 193.

2. Агапов П. И. Численное моделирование патологических процессов при черепно-мозговых травмах // Тезисы 12-ой международной конференции «Математика. Компьютер. Образование».- Ижевск, 2005.-С. 71.

3. Агапов П. И., Петров И. Б. Численное моделирование патологических процессов при черепно-мозговых травмах // Труды XIV научной конференции МФТИ «Современные проблемы фундаментальных и прикладных наук». Часть VII «Прикладная математика и экономика». - М.: МФТИ, 2002. - С. 42.

4 Агапов П. И., Петров И. Б. Численное исследование механических факторов черепно-мозговой травмы // Моделирование и обработка

информации: Сб. ст. / Моск. физ.-тех. ин-т. — М.: МФТИ, 2003. С. 28 - 32.

5. Агапов П. И., Петров И. Б. Численное моделирование патологических процессов при черепно-мозговых травмах сеточно-характеристическими методами // Труды XLVII научной конференции МФТИ «Современные проблемы фундаментальных и прикладных наук». Часть VII «Управление и прикладная математика». — М.: МФТИ, 2004. - С. 44.

6. Агапов П. И., Петров И. Б. Расчет областей повреждения мозга при черепно-мозговой травме // Компьютер и мозг. Новые технологии: Сб. ст. / РАН. - М.: Наука, 2005. - С. 28 - 38.

7. Агапов П. И., Петров И. В., Челноков Ф. Б. Численное исследование задач механики деформируемого твердого тела в неоднородных областях интегрирования // Обработка информации и моделирование: Сб. ст. / Моск. физ.-тех. ин-т. — М., 2002. — С. 148 - 157.

8. Агапов П. И., Челноков Ф. Б. Сопоставление нескольких разностных схем для численного решения двумерных гиперболических уравнений в частных производных // Труды XLVI научной конференции МФТИ «Современные проблемы фундаментальных и прикладных наук». — М .-Долгопрудный: МФТИ, 2003. - С. 48.

9. Агапов П. И., Челноков Ф. Б. Сравнительный анализ разностных схем для численного решения двумерных задач механики деформируемого твердого тела // Моделирование и обработка информации: Сб. ст. / Моск. физ.-тех. ин-т. - М., 2003. - С. 19 - 27.

10. Компьютерное моделирование биомеханических процессов сеточно-характеристическим методом / П. И. Агапов, А. С. Обухов, И. Б. Петров, Ф. Б. Челноков // Управление и обработка информации: модели процессов: Сб. ст. / Моск. физ.-тех. ин-т. — М., 2001.— С. 95 - 114.

11. Численное решение динамических задач биомеханики сеточно-характеристическим методом / П. И. Агапов, А. С. Обухов, И. Б. Петров, Ф. Б. Челноков // Компьютерные модели и прогресс медицины: Сб. ст. / РАН. - М. -.Наука, 2001.- С. 275 - 300.

Тираж 80 экз. Подписано в печать 21.1О.20О5г. Типография ООО «Алфавит 2000» 2005 г.

--À-^l

РНБ Русский фонд

2007-4 5257

\

8 ФЕБ 70

Введение

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

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

Содержание работы.

Положения, выносимые на защиту.

Благодарности.

1 Математическая модель

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

1.1.1 Анатомическое строение головы человека.

1.1.2 Двумерные модели механической реакции головы человека на ударные нагрузки

1.2 Уравнения механики деформируемого твердого тела.

1.2.1 Уравнения движения и реологические соотношения

1.2.2 Матричная форма уравнений

1.2.3 Гиперболические свойства систем уравнений линейной упругости

2 Численные методы

2.1 Численные методы для решения системы уравнений механики деформируемого твердого тела.

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

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

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

2.1.4 Разностные схемы для неструктурированных треугольных сеток

2.2 Расчет граничных узлов и постановка граничных условий.

2.2.1 Расчет граничных узлов в одномерном случае.

2.2.2 Расчет граничных узлов в двумерном случае.

2.2.3 Примеры граничных условий

2.3 Решение задачи контактного разрыва.

2.3.1 Решение задачи контактного разрыва в одномерном случае

2.3.2 Решение задачи контактного разрыва в двумерном случае

2.3.3 Алгоритм обнаружения столкновений (collision detection)

2.3.4 Примеры контактных условий.

2.4 Перестройка расчетных сеток.

2.4.1 Структурированные четырехугольные сетки.

2.4.2 Неструктурированные треугольные сетки.

3 Результаты моделирования

3.1 Интегральная оценка механического воздействия.

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

3.2.1 Качественное описание механической реакции головы человека на ударную нагрузку

3.2.2 Сравнение с клиническими данными.

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

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

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

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

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

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

Научные результаты диссертации опубликованы в работах [1-11].

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

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

Научная новизна. Научная новизна работы заключается в следующем:

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

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

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

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

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

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

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

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

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

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

В работе [12] был предложен ряд двумерных моделей для трансверсально-го и сагиттального сечений для анализа влияния различий в механических свойствах белого и серого вещества на распределение максимальных сдвиговых напряжений. Автор показал, что различие в механических свойствах серого и белого веществ, а также включение в модель желудочков необходимо для соответствия вычисленного распределения максимальных сдвиговых напряжений наблюдаемым аксональным повреждениям в белом веществе. В данной модели желудочки описывались жидкостными элементами, а модуль Юнга белого вещества был на 60% больше, чем серого.

В исследовании [13] использовалась двумерная модель вертикального сечения головы человека. Контакт череп-мозг моделировался двумя способами: 1) жесткое сцепление, и 2) скольжение с возможностью отрыва. Кроме того, боль

Авторы Геометрия Реология Валидация Основные выводы

Ueno [15, 16] Сагиттальное сечение Недефор-мируемый череп, линейно-упругий мозг Давления [14] Сопоставление линейных ускорений с давлением, угловых ускорений со сдвиговым напряжением

Ruan [17] Фронтальное сечение Невязкая жидкость, линейная упругость Необходимость включения мембран

Willinger [18] Сагиттальное сечение Линейная упругость Затылочный удар приводит к лобовым ушибам

Chu [19] Сагиттальное сечение Линейная упругость Давления [14] Лучшая корреляция ушибов со сдвиговыми напряжениями, чем с давлениями

Kuijpers [13] Сагиттальное сечение Линейная упругость, вязкоупругость Давления [14] Давления больше зависят от контакта череп-мозг, чем от наличия большого отверстия

Zhou [12] Фронтальное сечение Линейная упругость, вязкоупругость Максимальные сдвиговые деформации соотносятся с местами поражения

Miller [20] Осесим-метрич-ная Недеформиру-емый череп, линейная упругость и гипер-вязкоупругость Контакт с трением лучше соответствует поражениям при дифф. ак. повреждении

Таблица 1. Сводная таблица двумерных конечно-элементных моделей головы человека шое отверстие также моделировалось тремя способами: 1) жесткое сцепление, 2) скольжение с возможностью отрыва, и 3) незакрепленное. Сравнение с экспериментом [14] показало наилучшее соответствие величины давления в месте удара в случае скользящего контакта, в то время как модели с жестким сцеплением показали результат, наиболее далекий от действительности. Ни одна из моделей не показала соответствия давлений в месте противоудара, что может объясняться неспособностью контактного алгоритма противостоять растягивающим нагрузкам. Авторы заключили, что величины давления как в месте удара, так и в месте противоудара намного более чувствительны к способу моделирования контакта череп-мозг, чем к наличию или отсутствию свободного большого отверстия.

В работе [21] исследуются сдвиговые характеристики мозгового вещества на основе анализа экспериментальных данных, полученных с помощью физической модели. Показано, что вязкоупругая реология, использованная для моделирования мозгового вещества в большинстве работ [12,22,23], справедлива лишь для деформаций, не превышающих 1%. Предложена нелинейная вязко-упругая модель на основе уравнений гиперупругости Муни-Ривлина. Нелинейная модель мозгового вещества применена для конечно-элементной модели из [22], и было показано улучшение соответствия линейных смещений мозгового вещества при фронтальном ударе с экспериментальными данными [14].

Смещение мозга относительно черепа описаны более пятидесяти лет назад. Это вызывало особый интерес к численному моделированию этого явления. Моделирование внутричерепной жидкости линейно-упругими элементами с малым сдвиговым модулем [12,24] показало, что этот подход вносит дополнительные касательные напряжения между черепом и мозгом, а также снижают устойчивость вычислений вследствие вырождения элементов. Другим подходом является введение контактных условий на границе череп-мозг. Конкретные условия различаются у разных авторов, от жесткого сцепления до скольжения с трением. Было проведено некоторое количество сравнительных исследований с использованием различных контактных условий [22]. Все они приводят к заключению, что эффект ударного воздействия на голову человека чувствителен к способу моделирования взаимодействия череп-мозг.

Более точная трехмерная геометрическая модель черепа была использована в [23] для оценки существующих критериев повреждаемости мозга. Модель включает в себя мембраны твердой оболочки, сагиттальные синусы и эмиссар-ные вен. Численное исследование проводилось с помощью пакета нелинейного конечно-элементного анализа LS-DYNA3D. Автором предложены способы коррекции критерия HIC (Head Injury Criteria), применяющегося в машиностроении, в зависимости от размера головы, размера мозга и направления удара.

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

Содержание работы

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

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

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

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

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

Обсуждаются методы построения адаптивных расчетных сеток и их перестройки в процессе расчета.

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

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

Положения, выносимые на защиту

На защиту выносятся:

Основные результаты и выводы диссертации

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

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

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

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

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

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

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

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

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

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

По итогам анализа полученных расчетных данных сделан ряд выводов о характере пространственных распределений различных типов механических воздействий (сжатий, растяжений и сдвиговых нагрузок):

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

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

3. Характер распределения максимальных сжатий и растяжений относительно оси удара практически не изменяется для моделей без мембраны.

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

5. Для модели с мембраной сжатия и растяжения в противоударной области при боковом ударе уменьшаются примерно вдвое по сравнению с моделями без мембраны. Это можно объяснить экранирующим действием мембраны, отражающей продольные упругие волны и сдерживающей смещения мозгового вещества.

6. Сдвиговые нагрузки для модели с мембраной увеличиваются в 1.25-1.5 раз по сравнению с моделями без мембраны. Это можно объяснить более сложной картиной деформации мозгового вещества вокруг отростков мембраны при их «обтекании».

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

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

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

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

Заключение

1. Агапов П. И., Петров И. Б. Расчет областей повреждения мозга при черепно-мозговой травме // Компьютер и мозг. Новые технологии: Сб. ст. / РАН.— М.: Наука, 2005. - С. 28 - 38.

2. Агапов П. И. Анализ результатов численного моделирования черепно-мозговой травмы // Процессы и методы обработки информации: Сб. ст. / Моск. физ.-тех. ин-т,- М.: МФТИ, 2005.- С. 186 193.

3. Агапов П. И., Петров И. Б. Численное исследование механических факторов черепно-мозговой травмы // Моделирование и обработка информации: Сб. ст. / Моск. физ.-тех. ин-т. М.: МФТИ, 2003. - С. 28 - 32.

4. Агапов П. И., Челноков Ф. Б. Сравнительный анализ разностных схем для численного решения двумерных задач механики деформируемого твердого тела // Моделирование и обработка информации: Сб. ст. / Моск. физ.-тех. ин-т. М., 2003. - С. 19 - 27.

5. Агапов П. И., Петров И. Б., Челноков Ф. Б. Численное исследование задач механики деформируемого твердого тела в неоднородных областях интегрирования // Обработка информации и моделирование: Сб. ст. / Моск. физ.-тех. ин-т. М., 2002. - С. 148 - 157.

6. Численное решение динамических задач биомеханики сеточно-характеристическим методом / П. И. Агапов, А. С. Обухов, И. Б. Петров, Ф. Б. Челноков // Компьютерные модели и прогресс медицины: Сб. ст. / РАН. М.: Наука, 2001. - С. 275 - 300.

7. Компьютерное моделирование биомеханических процессов сеточно-характеристическим методом / П. И. Агапов, А. С. Обухов, И. Б. Петров,

8. Ф. Б. Челноков // Управление и обработка информации: модели процессов: Сб. ст. / Моск. физ.-тех. ин-т. — М., 2001. — С. 95 114.

9. Агапов П. И. Численное моделирование патологических процессов при черепно-мозговых травмах // Тезисы 12-ой международной конференции «Математика. Компьютер. Образование». — Ижевск, 2005.— С. 71.

10. Zhou С., Khalil Т. В., King A. I. A new model comparing impact responses of the homogeneous and inhomogeneous human brain / / 39th Stapp Car Crash Conf. / Society of Automotive Engineers. — 1995. — Pp. 121 137.

11. Kuijpers A. H., Claessens M. H., Sauren A. A. The influence of different boundary conditions on the response of the head to impact: a two-dimensional finite element study // J. Neurotrauma. — Vol. 12, no. 4. — Pp. 715 724.

12. Nahum A. M., Smith R. W., Ward С. C. Intracranial pressure dynamics during head impact // 21th Stapp Car Crash Conf. / Society of Automotive Engineers. — 1977.

13. Two dimensional finite element analysis of human brain impact responses: Application of scaling law / K. Ueno, J. W. Melvin, E. Lundquist, M. C. Lee // BED vol. 13 ASME. 1989. - Pp. 123 - 124.

14. Ueno K., Melvin J. W.and Lee M. C., Lighthall J. W. Development of tissue level brain injury criteria by finite element analysis //J. Neurotrauma. — 1995. — Vol. 12, no. 4. Pp. 695 - 706.

15. Ruan J. S., Khalil Т., King A. I. Human head dynamic response to side impact by finite element modeling // J. Biomechanical Engineering. — Vol. 113, no. 3.— Pp. 276 283.

16. Willinger R. Modal analysis of a finite element model of the head // IRCOBI Conf. Verona: 1992. - Pp. 283 - 297.

17. Finite element analysis of cerebral contusion / C. Chu, M. Lin, H. M. Huang, M. C. Lee // J. Biomechanics. 1994. - Vol. 27. - Pp. 187-194.

18. Miller R. T. e. a. Finite element modeling approaches for predicting injury in an experimental model of severe diffuse axonal injury // 42nd Stapp Car Crash Conf. / Society of Automotive Engineers. — 1998. — Pp. 155 166.

19. Brands D. W. A. Predicting brain mechanics during closed head impact: numerical and constitutive aspects: Ph.D. thesis / Eindhoven University of Technology. — 2002.

20. Claessens M. H. A. Finite Element Modeling of the Human Head under Impact Conditions: Ph.D. thesis / Eindhoven University of Technology. — 1997.

21. Kleiven S. Finite Element Modeling of the Human Head: Ph.D. thesis / Department of Aeronautics, Royal Institute of Technology, Sweden. — 2002.

22. Zhang L., Yang К. H., King A. I. Comparison of brain responses between frontal and lateral impacts by finite element modeling //J. Neurotrauma. — Vol. 18, no. 1,- Pp. 21 30.

23. Регирер С. А. Лекции по биологической механике. — М.: МГУ, 1980.

24. Седов Л. И. Механика сплошной среды. — М.: Наука, 1970.

25. Новацшй В. К. Теория упругости. — М.: Мир, 1975.

26. Магомедов К. М., Холодов А. С. О построении разностных схем для уравнений гиперболического типа на основе характеристический соотношений // Ж. вычисл. машем, и матпем. физ. — 1969. — Т. 9, № 2. — С. 373 386.

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

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

29. Courant R., Isaacson Е., Rees М. On the solutions of nonlinear hyperbolic differential equations by finite differences // Comm. Pure and Appl. Math. — 1952. Vol. 5, no. 5. — Pp. 243 - 254.

30. Lax P. D., Wendroff B. System of conservation laws // Comm. Pure and Appl. Math. 1960. - Vol. 13, no. 2. - Pp. 217 - 237.

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

32. Shewchuk J. R. Delaunay refinement algorithms for triangular mesh generation / / Computational Geometry: Theory and Applications.— 2002.— Vol. 22(1-3), no. 5. — Pp. 21-74. http://www.cs.berkeley.edu/~jrs/papers/2dj.ps.

33. Иваненко С. А., Чарахчьян А. А. Криволинейные сетки из выпуклых четырехугольников // Ж. вычисл. матем. и матем. физ. — 1988. — Т. 28, X8 4. — С. 503 514.