CC BY
108
DOI: 10.15593/RZhBiomeh/2018.2.04 УДК 616.314.26-06:616.742.7-008.1-037
к Российский
Журнал / Биомеханики
www.biomech.ru
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ УДАРНОГО ВОЗДЕЙСТВИЯ НА ГОЛОВУ ПРИ ЧЕРЕПНО-МОЗГОВЫХ ТРАВМАХ
Институт механики Уральского отделения наук Российской академии наук, Россия, 426067, Ижевск, ул. Т. Барамзиной, 34, e-mail: karavaev-alexander@yandex.ru
Аннотация. Проведен обзор известных математических моделей головы человека, используемых в задачах контактной биомеханики. Построена конечно-элементная модель из объемных данных, состоящая из трех основных материалов: мягкая ткань, череп и мозг, созданных на основе компьютерной томографии с использованием разработанного авторами алгоритма. Для решения динамической задачи предложен итерационный метод декомпозиции Шварца для двух контактирующих трехмерных тел, позволяющий рассматривать решение задачи напряженно-деформируемого состояния в стандартной постановке для каждого из контактирующих тел. В соответствии с предложенным алгоритмом на каждой итерации выполняется два шага, где поочередно удовлетворяются условия сопряжения по перемещениям и напряжениям для узлов, расположенных на контактной поверхности. Проведен ряд вычислительных экспериментов по установлению критических величин контактной силы, приводящих к тяжелым последствиям. При значениях силы удара в диапазонах 4000-6000 Н можно говорить о получении травмы средней и серьезной тяжести, а величины силы выше 7000 Н расцениваются как удары, влекущие тяжелые последствия с длительной потерей сознания; при силе, равной 10000 Н, с большой долей вероятности наступит смерть потерпевшего.
Ключевые слова: математическое моделирование, контактное взаимодействие, биомеханика, черепно-мозговая травма, метод конечных элементов.
Травмы головы происходят относительно часто и составляют от 3 до 6% всех несчастных случаев в развитых странах. Нередко они бывают тяжелыми и приводят к значительным потерям рабочего времени по нетрудоспособности, коме и даже летальному исходу. Причиной травм являются дорожно-транспортные происшествия, спортивные, бытовые и иные несчастные случаи, аварии на производстве, удары головой во время падения или о какие-либо движущиеся предметы и т.д.
Травма головы - широкий медицинский термин, который относится к таким повреждениям, как:
- перфорация черепа при действии значительной силы на небольшой участок головы, например при ударе предметом с заточенным острием или острыми краями;
- перелом костей черепа или шейных позвонков при обширном сильном ударе, если его сила превысила предел упругой деформации черепа или сжатых шейных позвонков;
А.С. Караваев, С.П. Копысов
Введение
© Караваев А.С., Копысов С.П., 2018
Караваев Александр Сергеевич, научный сотрудник, Ижевск
Копысов Сергей Петрович, д.ф.-м.н., профессор, главный научный сотрудник, Ижевск
- поражение головного мозга без перелома костей черепа в результате смещения мозга в черепной коробке, которое приводит к контузии, сотрясению, внутреннему кровотечению и нарушению кровообращения мозга.
Выявление физических причин возникновения той или иной травмы весьма затруднительно, и современные исследования по этому вопросу не дают однозначного ответа. Некоторые специалисты считают основным травмирующим фактором силу, другие - энергию или количество движения. Возникновение травмы обусловливают также ускорение, разность ускорений или сила удара. На практике, по-видимому, каждый фактор в той или иной степени влияет на возникновение травмирующего воздействия. Во многих случаях данных о полученном ударе при травмах головы недостаточно и они противоречивы.
Сила удара, при которой может произойти травма головы, определяется экспериментально на трупах или на животных, однако полученные данные трудно экстраполировать на человека.
Простейшие сведения из теории ударного контактного взаимодействия показывают, что даже при рассмотрении случаев соударения однородных тел без их разрушения определение силы удара является весьма сложной задачей. При статической нагрузке свод черепа может выдержать без разрушения давление от 11 760 до 17 640 Н и выше [9]. При динамическом ударном воздействии дополнительно рассматривается работа удара, кинетическая энергия травмирующего предмета, удельная работа удара, импульс силы и прочее [8].
Моделирование процессов соударения является важнейшей областью исследования и находит широкое применение при оценке прочностных свойств защитных головных шлемов, систем безопасности автомобильного транспорта, определении в тех или иных случаях предельно допустимых норм ударных нагрузок и т.д. Решение такой задачи подразумевает вычисление деформаций биологических тканей, подвергнутых внешнему воздействию, с последующим определением возможных травматических последствий.
Применение для исследований математических моделей исключает необходимость проведения большого числа натурных экспериментов и позволяет рассматривать сложные виды контактного взаимодействия, воспроизведение которых на практике зачастую является очень трудоемким процессом. Кроме того, в данном случае существует возможность вычисления параметров нагрузок и деформаций, которые не могут быть определены в реальных экспериментах.
Математические модели головы человека
На раннем этапе исследований двумерные и трехмерные конечно-элементные модели головы человека упрощались представлением лишь двух тканей, черепа и мозга, построенных в виде эллипсоидных и сферических тел.
В начале 90-х гг. с развитием вычислительных возможностей компьютеров стали появляться более сложные и реалистичные 3.0-модели. Руан и др. [33] построили первый вариант хорошо известной модели для изучения травм головного мозга (ЖШВ/М). Расчетная сетка Ж5ЦВ/М содержала 6080 узлов и 7351 элементов. Череп, мозг и спинномозговая жидкость были построены из шестигранных элементов, а скальп, твердая мозговая оболочка и серп мозга описывались четырехугольными поверхностными ячейками. Модель неоднократно дорабатывалась путем более детального моделирования различных участков мозга и головы [40], а также введением контактного взаимодействия между черепом и мозговой оболочкой [22].
В 2002 г. Харди с соавт. [24] разработали модель головы взрослого мужчины, включающую дополнительно одиннадцать пар мостовых вен и упрощенную модель шеи.
Хорган и Гилкрист построили модель для рассмотрения простых аварий (UCDBTM) с пешеходами, включающую скальп, трехслойную структуру черепа, состоящую из различных частей мозга, а также слоя спинномозговой жидкости [18].
В работах И.Б. Петрова с соавт. [1] построена двумерная математическая модель механической реакции головы человека на ударные воздействия, описывающая пространственное распределение механических нагрузок на мозг. В дальнейшем авторы усовершенствовали разработанную модель, изменяя реологические свойства биоматериалов [3].
В целом разрабатывалось достаточное количество прототипов с различной степенью детализации границ и свойств биологических тканей. Как правило, построение таких моделей достаточно трудоемко ввиду необходимости ручного либо полуавтоматического проектирования контуров биологических органов на основании реально существующих образцов. Во многих случаях выделение дополнительных отделов мозга осуществлялось заданием разных свойств ячеек сетки без реального восстановления границы между материалами.
Отдельную группу составляют модели, созданные на основе объемных данных. В 1995 г. по данным магнитно-резонансной томографии Виллингер и др. [36] создали трехмерную модель для построения альтернативной манекен-модели головы. Бандак и др. [10] использовали данные компьютерной томографии для построения моделей с различной степенью разрешения и моделирования ударной нагрузки, которая показала, что максимальные напряжения возникают на внутренней поверхности черепа в точке удара.
В 2005 г. Белингарди с соавт. [11] предложили численную модель, построенную по данным компьютерной томографии, и установили, что включение в модель тенториальных мембран уменьшает расчетное пиковое внутричерепное давление в лобной и затылочной областях на 17 и 18% соответственно.
В работе [30] с использованием пакета MIMICS и исходных томографических данных построена трехмерная модель, состоящая из мозга, черепа для моделирования теста Наума.
В 2014 г. на основе поверхностной треугольной сетки, полученной непосредственно из двумерных топографичеких срезов, Янг и др. [38] разработали сложный прототип, содержащий одиннадцать различных типов тканей. Важной особенностью этой работы являлось рассмотрение спинномозговой жидкости как гидростатической жидкости, заполняющей свободное пространство между черепом и мозгом, что отличается от стандартного представления ее линейно-упругим материалом с коэффициентом Пуассона, близким к 0,5, что привело к увеличению возникающих в мозге максимальных напряжений и деформаций приблизительно на 18%.
Важно отметить следующее: несмотря на то что человеческие ткани считаются сложной, многофазной, гетерогенной и анизотропной структурой, в большинстве конечно-элементных моделей они рассматриваются как изотропные, поскольку на данный момент практически невозможно численно оценить анизотропные свойства того или иного биологического органа.
Как правило, все материалы, за исключением мозга, задаются линейно-упругими, а мозг, в свою очередь, имеет линейно-упругое или линейно-вязкоупругое поведение. Несмотря на существование по этому поводу различных подходов, следует отметить, что отличие в механических свойствах между линейной и вязкоупругой моделями увеличивается пропорционально времени проведения эксперимента.
В задачах кратковременного ударного взаимодействия применение упругих моделей не приводит к существенным погрешностям, однако не следует пренебрегать вязкоупругими свойствами в приложениях с небольшой скоростью возрастания нагрузок [26].
Кроме того, линейно-упругие свойства тканей головного мозга, получаемые в экспериментальных исследованиях, значительно отличаются. Так, в обзоре [12] модуль Юнга мозговой ткани изменяется в диапазонах от 6,67-102 до 68 МПа, коэффициент Пуассона варьируется в пределах 0,48-0,4996, плотность остается в среднем постоянной и составляет 1000-1040 кг/м3, также отмечаются модели более высокой плотности р = 1220 кг/м3.
Таким образом, мозг представляет собой сложную неоднородную биологическую структуру, механические свойства которой до конца не изучены, и представление ее в виде линейно-упругой или вязкоупругой модели является лишь некоторым приближением реальных физических свойств [12, 43]. Так, эмпирически установлено, что белое вещество на 60% жестче серого вещества, что объясняется наличием в нем аксональных волокон, отдельные части головного мозга (серп, желудочки, мозжечок и т.д.) также имеют различные механические свойства, во многих случаях сильно отличающиеся друг от друга [26]. Выбор тех или иных значений линейно-упругих свойств во многом связан типом проводимого вычислительного эксперимента. В некоторых работах [41, 42] модули Юнга тканей головного мозга отличались по значениям в десятки раз, при этом сама конечно-элементная модель головы не изменялась. В одном случае моделировалось ударное воздействие, в другом - угловое ускорение.
Немаловажным фактором [23], значительно влияющим на результаты вычислительного эксперимента, является способ задания взаимодействия между черепом и мозгом в результате полученной ударной нагрузки, приводящей к их относительному смещению.
Наиболее простая постановка задачи заключается в задании общей границы, состоящей из узлов, прилегающих к расчетным ячейкам как первого, так и второго материала [10, 19]. Иногда в моделях строится разделяющий слой спинномозговой жидкости, обладающий, как правило, линейно-упругими свойствами, близкими к несжимаемому материалу, и малым значением модуля Юнга [24, 30, 33, 36, 38]. Другой подход заключается в использовании различных типов контактных алгоритмов между поверхностями материалов, которые могут исключать взаимное движение контактных узлов либо допускать скольжение [1, 3, 12].
Использование слоя спинномозговой жидкости ведет к возникновению высоких напряжений сдвига и деформации в области контакта. Подобные модели наиболее адекватно описывают изменение внутричерепного давления при воздействии внешнего удара. Биомеханические модели, использующие контактные взаимодействия, предсказывают хорошее соотношение значений внутричерепного давления в лобной области, где задается ударная нагрузка, однако недостаточно хорошо воспроизводят экспериментальные результаты в задней части головы, где возникают отрицательные значения внутричерепного давления. Как показывают результаты представленного исследования, сходным поведением также обладают модели с типом взаимодействия по общей границе.
Конечно-элементная модель из объемных данных
Моделирование процессов биомеханики связано с рассмотрением объектов сложной формы и необходимостью построения границ, описывающих внутреннюю
структуру биологических тканей. Для повышения реалистичности проводимого вычислительного эксперимента предпочтительным является использование алгоритмов построения конечно-элементных моделей непосредственно по данным компьютерной или магнитно-резонансной томографии.
В данной работе конечно-элементная модель строилась на основании объемных данных и состояла из трех основных материалов: мягкая ткань, череп и мозг, созданных на основе компьютерной томографии с использованием алгоритма [5, 7] (рис. 1). Метод построения включал следующие основные этапы:
• построение уровней октодерева, используемого для хранения точек в пространстве, и уровней детализации объекта;
• генерация однородной шестигранной сетки на основе алгоритма двойственного контурирования, позволяющего сохранять острые углы и грани объекта;
• улучшение качества построенных ячеек с помощью различных процедур сглаживания узлов сетки, а также вставки дополнительного слоя граничных ячеек.
Расчетная неструктурированная шестигранная сетка содержит 1 363 369 узлов и 1 308 341 ячейку и характеризуется следующими мерами качества сетки: минимальный нормированный якобиан ячейки равняется Q = 0,2, что является допустимым в
конечно-элементных расчетах; максимальное аспектное соотношение - Q = 15; угол
искривленности четырехугольной грани 0w = 36,6°; метрика Одди Q = 178 [42].
Отметим, что наряду с использованием стандартных процедур оптимизации качества ячеек сетки, таких как вставка буферного слоя и сглаживания, применялись специальные функции корректировки формы области материала и удаления «зашумленных» участков реальной физической модели: заполнение внутренних полостей, поиск «висячих» ячеек, не связанных с основной моделью и т.д.
Свойства тканей лица, слоя скальпа и черепа, рассматриваемого как единое целое, приняты линейно-упругими на основании значений из работ [20, 35, 37] (табл. 1).
Введение материала лицевой ткани позволило скорректировать массу головы и привести ее к диапазону средних значений 4,5-5 кг, составляющих около 7-8% массы тела взрослого человека [29]. Масса построенной биомеханической модели составила 4,9 кг, в то время как мозг весил 1 ,4 кг, что также находится в пределах средних показателей (1,3-1,4 кг) [25].
В данной работе константы линейно-упругого материала для ткани мозга заимствованы из [19]. Предполагалось, что взаимодействие между черепом и мозговой оболочкой происходит через общую границу, что позволяет провести сопоставление результатов тестирования.
Таблица 1
Линейно-упругие свойства материалов биомеханической модели
Ткань Е, МПа V р, кг/м3 Источник
Слой скальпа 3,45 0,4 1200 Халил и др. [26]
Мозг 0,25 0,49 1000 Ханг и др. [21]
Череп (как целое) 6000 0,21 2100 Виллингер и др. [25]
Ткани лица 5000 0,23 2500 Виллингер и др. [24]
Рис. 1. Неструктурированные сетки с шестигранными ячейками для областей: мягкой ткани (а), мозга (б), черепа (в), конечно-элементная модель головы для трех тканей с увеличенным фрагментом сетки (г)
Конечно-элементная валидация по данным физических экспериментов
Проверка по данным реального физического эксперимента является немаловажным этапом разработки любой конечно-элементной модели, имитирующей свойства реального биомеханического объекта. В настоящее время в натурных экспериментах разработаны и используются несколько различных моделей черепно-мозговых травм, среди которых жидкостно-перкуссионая травма мозга, модель контролируемого коркового повреждения, травма инерционного ускорения и т.д.
Для оценки механических параметров, определяющих тяжесть повреждения, используют компьютеризированное измерение прилагаемой нагрузки на основе следующих физико-механических параметров: градиент давления воздействующей на мозг жидкости, скорость ударника, скорость прилагаемой нагрузки, вызывающей ускорение-замедление головы [2].
В большинстве исследований черепно-мозговых травм [12, 19, 23, 26, 30, 33, 38, 40, 43], не рассматривающих перелом костей черепа, тестировали модели, построенные с использованием данных внутричерепного давления, которое долгое время оставалось
единственным параметром, измеренным на практике в натурных экспериментах Наума [28] и Троселли [34] с человеческим трупом.
Многие исследователи считают более естественным сопоставление данных по возникающим деформациям. В этом случае диффузное аксональное повреждение тканей головного мозга в большей степени зависит от деформаций, чем от внутричерепного давления [26]. Такие данные стали доступны в работе [14], содержащей экспериментальные результаты относительного смещения между черепом и мозгом.
В тесте Наума приводятся значения массы и начальной скорости ударника, при этом остальные его физическо-механические характеристики не представлены. При моделировании линейно-упругие свойства ударника подбирались так, чтобы величина возникающей ударной нагрузки совпадала с экспериментальными данными. Несмотря на недостаточно подробное описание проведенного эксперимента Наума, его моделирование выполняется практически во всех исследованиях. В эксперименте ударная нагрузка является кратковременной, и, согласно работам [32, 37], влияние шеи на возникающее в мозге внутричерепное давление незначительно. Следовательно, данный тест также подходит для апробации моделей головы человека, разработанных без учета шейного отдела позвоночника.
Таким образом, в вычислительном эксперименте моделирование выполнялось при свободных граничных условиях на голове и задании только начальной скорости ударника. Физический же эксперимент состоял в следующем. Сидящий труп, удерживаемый от падения закрепленным в верхней части торсом, подвергался удару в область лобной кости вдоль сагиттальной плоскости головы цилиндрическим ударником массой 5,59 кг с постоянной скоростью 9,94 м/с. Голова была повернута по анатомической плоскости Франкфурта вперед на 45 градусов (рис. 2, а). Изменения давления во время удара измерялись и регистрировались с использованием датчиков, расположенных в слое спинномозговой жидкости для пяти областей головного мозга: лобной доли, где происходило воздействие; теменной; правой и левой затылочных долей; а также области задней ямки, состоящей из мозжечка и спинного мозга (см. рис. 2, б). На рис. 2, в показаны точки замера контрольных значений давления для каждой из областей в построенной расчетной модели. Так как слой спинномозговой жидкости не рассматривался, были выбраны точки непосредственно внутри мозга.
. ... ,
■ 5
1 .
.2, 3
а б в
Рис. 2. Схема эксперимента Наума: а - физический эксперимент; б - структура головного мозга (1 - лобная доля, 2 - теменная доля, 3 - затылочная доля, 4 - мозжечок, 5 - спинной мозг, 6 - височная доля); в - контрольные точки замера давления в вычислительном эксперименте (1 - лобная доля, 2, 3 - затылочные,
4 - область задней ямки, 5 - теменная)
1
3
4
4
Критерии оценки травматизма биологических тканей
К настоящему времени накоплен большой экспериментальный материал, описывающий характер, количество и распространенность повреждений при ряде произвольно выбранных условий удара [4]. Вероятность травматизма биологических тканей (в том числе сотрясения мозга) можно с достаточной степенью достоверности оценить использованием различных критериев черепно-мозговых травм.
Наиболее часто в тех или иных случаях используются критерии травматизма, вычисленные на основе ускорения центра тяжести головы: HIC (Head Injury Creterion), GSI (Gadd Severity Index), оценка максимального линейного и вращательного ускорения и т.д.
Оценки, основанные на значениях напряжений, деформаций и внутричерепного давления, применяются реже.
В работе [39] исследована возможность возникновения черепно-мозговой травмы от внутричерепного давления в области удара в диапазоне 66-114 КПа -наличие травмы, при 44-78 КПа - отсутствие травмы. При значениях на противоположенной стороне черепной коробки в направлении ударной нагрузки от -101 до -51 КПа - травма, при -59 до -23 КПа - без травмы. Возникновение внутричерепного давления более 200 КПа является индикатором контузии и гематомы мозговых тканей [30]. Кроме того, для обнаружения некоторых средних и тяжелых травм использования только параметра ускорения головы недостаточно и необходимо измерение пороговых значений внутричерепного давления [27].
Необходимо учитывать некоторые рекомендуемые значения в терминах возникающих напряжений и деформаций: максимальное напряжение сдвига 8-16 КПа - серьезная травма мозговых тканей [39]; значение интенсивности напряжения по Мизесу более 27 КПа - аксональное повреждение [13], максимальная деформация более0,19 - возникает 50%-ная вероятность получения травмы средней тяжести [39].
Силу удара можно классифицировать следующим образом [9]: до 160 Н -небольшая сила удара; от 160 до 1960 Н - значительная сила удара; от 1960 до 4900 Н -большая сила удара; более 4900 Н - очень большая сила удара.
Одним из первых критериев оценки потенциальной опасности возникновения
T
травмы был индекс SI, определяемый как SI = ja (t)2'5 dt, где T - время соударения,
0
a(t) - ускорение в момент времени t, значение степени 2,5 было выбрано опытным
путем [16]. Индекс SI порой показывает необоснованно высокие значения для медленных и длительных по времени соударений и таким образом предсказывает более высокую вероятность травмирования в сравнении с результатами, полученными в натурных экспериментах.
Критерий HIC [15], включающий эффекты как ускорения головы, так и продолжительности ускорения, позволяет избегать завышения вероятности получения травмы и является показателем вероятности появления повреждений мозговых тканей при возникновении ударной нагрузки
\ 2,5 ^
I 2 . .
HIC = max
у h
1
— ja (t) dt (t2 - ^ )
где a(t) - ускорение в момент времени t, величины ^ и t2 - начальное и конечное время интервала, в течение которого HIC достигает максимального значения.
Продолжительность интервала \tx, t2 ] ограничена диапазонами значений от 3 до 36 мс,
как правило, выбирается значение 15 мс [15].
Установленные соотношения между значениями HIC и тяжестью полученной черепно-мозговой травмы широко используются при оценке степени безопасности автомобильного транспорта, средств индивидуальной защиты и спортивного инвентаря. Как правило, параметр ускорения измеряется акселерометром, установленным в центре тяжести головы манекенов, используемых для краш-тестов. Считается, что значения HIC больше 1000 являются опасными для жизни. При HIC « 1000 вероятность тяжелой черепно-мозговой травмы составляет 18%, серьезной - 55%, травмы средней тяжести - 90% [16]. Однако некоторые автомобильные краш-тесты показывают критические значения HIC меньше 142 (пример - тестирование «Ауди 8» с подушкой безопасности в 1995 г.) [15, 16].
Для унификации результатов проводимых экспериментов также используется кодированная шкала оценки степени тяжести полученной травмы (AIS). Диапазон кодов в шкале изменяется от 0 (нет повреждений) до 6 (летальный исход) [30, 31]. В табл. 2 представлены соотношения шкалы AIS с полученным максимальным линейным ускорением головы и критерием HIC.
Таблица 2
Соотношение шкалы AIS с ускорением и критерием HIC
Шкала AIS Максимальное ускорение HIC Оценка степени травмы
0 50g 134 Нет травмы
1 50-100g" 135-519 Легкая - головная боль или головокружение
2 100-150g 520-899 Средняя - потеря сознания менее часа
3 150-200g 900-1254 Серьезная - потеря сознания 1-6 часов
4 200-250g 1255-1574 Тяжелая - потеря сознания в течение 6-24 часов
5 250-300g 1575-1859 Критическая - потеря сознания в течение более 25 часов, большая гематома
- - 1860 Летальный исход
Вычислительный алгоритм
В динамике контактного взаимодействия обычно применяется пространственно-временная дискретизация, включающая в себя конечно-элементную аппроксимацию в пространстве и пошаговую схему интегрирования по времени. В рассматриваемых задачах аппроксимации динамического контактного взаимодействия добавляется схема учета контактных сил. Эффективный совместный выбор всех этих схем затруднен. В проведенном выше обзоре не приводится каких-либо особенностей вычислительных алгоритмов решения динамических контактных задач при взаимодействии однородного ударника с многослойной средой с существенно различающимися физико-механическими характеристиками. Во многих случаях расчеты выполняются
в стандартных пакетах, в которых реализованы классические схемы интегрирования, и порой неподготовленными пользователями, а представляемые неполные данные и результаты порой вызывают обоснованное недоверие.
Применение классических схем пошагового интегрирования во времени в случае контактных задач осложняется тем, что точное выполнение условий на контактной границе при неявной схеме Ньюмарка для интегрирования по времени приводит к нефизичным колебаниям при существенном отклонении полной энергии в системе. Во многом возникновение таких колебаний связано с тем, что при стандартной аппроксимации матрицы масс всем узлам, в том числе расположенным на контактной поверхности, приписывается некоторая собственная масса, которая вносит вклад в полную энергию системы. В момент времени, соответствующий контактному взаимодействию, скорость узлов на поверности взаимодействия становится равной нулю, как и их кинетическая энергия, вызывая перераспределение вкладов энергии в полную энергию системы, что приводит к значительным нефизичным численным колебаниям. Консервативные свойства схем сохраняются только для задачи линейной упругости с постоянной поверхностью контакта, в то время как для случая изменяющейся контактной границы появление нового узла в области контакта уменьшает энергию дискретной системы, а выход узла из контакта её увеличивает.
Эффективные методы решения трехмерных динамических контактных задач характеризуются выбранной схемой численного интегрирования во времени, которая во многом определяет наличие паразитных колебаний, обусловленных тем, что нельзя достаточно точно воспроизвести вклад высокочастотных форм колебаний в решение динамической контактной задачи, сохранение полной энергии системы контактирующих тел, устойчивость решения.
Одним из решений по исключению паразитных осцилляций в области контакта является применение диссипативных схем. Как показывает расчетная практика, эффективным является использование обобщенных неявных методов на основе схемы Ньюмарка, таких как HHT (Шlber-Hughes-Taylor) - схема (а-метод), которая безусловно устойчива, имеет второй порядок точности и диссипативна для высоких частот [17]. Однако применение этой схемы для рассматриваемой математической модели ударного взаимодействия показало, что она теряет устойчивость.
Альтернативным подходом предотвращения контактных колебаний также является метод перераспределения массы на границе контактной области. Суть метода состоит в обнулении массы узлов в области контакта, в результате чего у них исключаются инерциальные характеристики. Данный способ значительно стабилизирует неизвестную контактную границу и может быть использован в комбинации с любыми численными схемами [21].
В данной работе в качестве контактного алгоритма был использован итерационный метод декомпозиции Шварца для двух контактирующих трехмерных тел, позволяющий рассматривать решение задачи напряженно-деформированного состояния в стандартной постановке для каждого из контактирующих тел и допускающий высокоуровневое распараллеливание процесса решения [6]. В соответствии с предложенным итерационным алгоритмом на каждой итерации выполняется два шага, где поочередно удовлетворяются условия сопряжения по перемещениям и напряжениям для узлов, расположенных на контактной поверхности.
Особенностью реализации динамической контактной задачи являлось совмещение итерационных процедур уточнения решения контактной задачи по алгоритму Ньютона-Рафсона непосредственно с итерациями метода декомпозиции Шварца. Как показали результаты тестовых расчетов, применение контактного алгоритма Шварца и HHT (Шlber-Hughes-Taylor) - схемы в комбинации с методом
коррекции массы узлов (SCA-HHT-MR) на контактной границе - является надежным, устойчивым и достаточно точным способом численного интегрирования контактных задач при ударном взаимодействии.
Результаты вычислительных экспериментов
Прежде чем перейти к сопоставлению значений H/C-критерия с данными линейного ускорения и внутричерепного давления в зависимости от значений приложенной контактной силы, рассмотрим тестирование расчетной модели и динамического контактного алгоритма.
Расчетные данные сравнивались с данными реального физического эксперимента Наума [28].
Как видно из рис. 3, а, максимальное сжимающее контактное напряжение возникает на 4 мс эксперимента и составляет около 3 МПа, при этом исходя из количества точек в контакте (рис. 3, б) этот момент также соответствует максимальной площади соприкосновения соударяемых тел.
Отклонение значений полной энергии в сравнении с исходной ее кинетической частью 276 Дж составило не более 2%, что можно считать удовлетворительным результатом (рис. 4).
На рис. 5, а наблюдается достаточно равномерное распределение расчетного поля давления на поверхности мозга, начиная от минимальных растягивающих отрицательных значений в затылочной области и заканчивая максимальными сжимающими положительными величинами в области удара.
Запаздывающие перемещения тканей мозга в сравнении с остальными частями модели можно видеть на рис. 5, б. Данный факт обусловлен более низкой скоростью распространения продольных волн внутри мозга.
В результате численного моделирования эксперимента Наума максимальная контактная сила, равная 6882 Н, хорошо согласуется с экспериментальными данными (6900 Н), при этом получено более высокое значение линейного ускорения 2154 м/с2 в сравнении с экспериментальным значением 2040 м/с (рис. 6).
Далее рассмотрим результаты распределения величин давления, замеренных в контрольных точках, изображенных на рис. 2, в.
P, Па
2,0E+06
5,0E+05
0,000 0,002 0,004 0,006 0,008 0,010 0,012
а б
Рис. 3. Зависимость от времени в алгоритме SCA-HHT-MR: а - максимального контактного давления в точке контакта; б - числа контактных точек
n
3.5E+06
3.0E+06
2,5E+06
1,5E+06
1,0E+06
0
0,0E+00
0,000
0,002
0,004
0,006
0,008
0,010
0,012
E, Дж
300 250 200 150 100 50 0
0,000
Рис. 4 Изменение энергии в алгоритме SCA-HHT-MR: полная - сплошная кривая; кинетическая - штриховая кривая, внутренняя - штрихпунктирная кривая
& X
P, Па 117034 100000
0
¡-100000 -133272
u, м
0,004421 0,004
0,003
10,002
0,001 0,000822
Рис. 5. Расчетные поля давления на поверхности мозга (а) и модуля перемещений в сечении на 0,004 с (б)
б
а
а б
Рис. 6. Сравнение экспериментальных (штриховая кривая - [30]) и численных результатов (сплошная кривая - SCA-HHT-MR): a - распределение модуля
контактной силы; б - ускорения
P, Па
160 000 140 000 120 000 100 000 80 000 60 000 40 000 20 000 0
-20 000
P,
80000
70000 60000 50000 40000 30000 20000 10000 0
10000
Па
/ \
1 \
/V 1 / А 1
/'Л; А\
/.' V/ \\
1 1 V /V Ii
Ji t, c
0,002 0,004 0,006 0,008 0,01 0,012
б
а
P, Па
P, Па
10000 0
-10000 -20000 -30000 -40000 -50000 -60000 -70000
t, c
в
г
д
Рис. 7. Сравнение экспериментальных (штриховая кривая - [28]) и численных результатов (сплошная кривая - SCA-HHT-MR): a - в лобной доле; б - теменной доле; в - левой затылочной доле; г - правой затылочной доле; д - в области задней
ямки
Расчетное значение в лобной доле имеет пиковую величину 107 КПа на 4 мс, что близко к значению 114 КПа, показанному в проверочном тесте из работы [19] (рис. 7, а). В натурном эксперименте максимальная величина составила 146 КПа, что примерно на 26% выше. Вместе с тем результаты, полученные в данной работе, так же, как и результаты в [19], можно признать удовлетворительными.
В работе [19] отсутствуют данные по отрицательным значениям давления в задней части головы. Однако для этих областей в рассматриваемой модели, как и в случае прототипов на основе контактного интерфейса взаимодействия между мозгом и черепом, наблюдаются завышенные по модулю величины внутричерепного давления мозга [12]. В левой затылочной области минимальный показатель давления отличается более чем на 10 КПа и составляет -60 КПа при тестовом значение около -50 КПа. Экстремальные значения давления для правой затылочной области примерно совпадают и равняются -50 КПа. Наиболее сильное отклонение давления зафиксировано в области задней ямки -104 КПа, что в 1,7 раза меньше расчетного -60 КПа. Пиковое значение давления в теменной области на 37% ниже, кроме того наблюдаются более сильные осцилляции данного показателя в сравнении с результатами натурного эксперимента (рис. 7, б-д).
В работе [12] рассматривают разные причины отличий в полученных показателей давления и ускорения.
Во-первых, разница в используемых объектах воздействия. Общая геометрия головы в числовой модели не совпадает с головами трупов, используемых в эксперименте. Кроме того, разработанный прототип не содержит представлений для различных структур черепа и биологических тканей, составляющих внутричерепное содержимое. Взаимодействие головы и шеи также не принимается во внимание. Представление границы черепа и мозга в виде общих узлов может сформировать слишком сильное отрицательное давление в задней области, вызванное растягивающими силами, оказываемыми черепом на более медленный в движении мозг, что создает значения намного ниже экспериментальных. Это также может являться причиной более низкого сжимающего давления во фронтальной области, так как жесткое закрепление мозга с черепом ограничивает движение узлов в граничной зоне. Более реалистичное давление в задней части головы наблюдается при использовании слоя спинномозговой жидкости или обеспечением возможности относительного движения черепа и мозга.
Во-вторых, различие в свойствах материала. Измерения проводились на трупах, для которых механическое поведение различных тканей неизвестно. Свойства материала, использованные для моделировании, взяты из [43].
Кроме того, в работах [1, 3] при моделировании бокового удара в височную область головы показано, что результаты распределения нагрузок внутри мозга сильно зависят от внутреннего строения его элементов. В частности, учет наличия желудочков практически не оказывает влияния на распределение областей максимального сжатия и растяжения, но существенно влияет на распределение сдвиговых нагрузок. Наличие мембраны более существенно для локализации областей сжатия растяжения при боковых ударах.
Также отметим, что измеренные в реальном эксперименте давления могут не располагаться в тех же местах, что и в численной модели. Более того, модель не содержит представления спинномозговой жидкости, в слое которой были расположены датчики, вместо этого давление было рассчитано в соответствующих частях внутри мозга. Относительно наблюдаемого различия в значениях ускорения можно сказать, что результаты других экспериментов [28] с близкими условиями удара показали значительно более высокое пиковое ускорение.
Были также выполнены численные эксперименты, в которых исследованы изменения начальной скорости и ее влияние на величину контактной силы в тестах Наума.
В табл. 3 представлены максимальные значения контактной силы, линейного ускорения и сжимающего давления в лобной доле мозга (см. рис. 2) для каждой из скоростей. На основании полученного линейного ускорения определялся
соответствующий индекс шкалы AIS и вычислялся критерий HIC. Для индексов 5 и 6 шкалы AIS найдено хорошее соотношение величин ускорения и критерия HIC в соответствии с диапазонами значений, приведенными в табл. 2. Вместе с тем для других индексов значение критерия HIC является несколько заниженным, что
объясняется вычислением критерия по интервалу , ^ ], равному времени удара, это не
во всех случаях обеспечивает определение максимального значения величины HIC.
Таблица 3
Зависимость ускорения и внутричерепного давления от приложенной контактной силы
Шкала AIS Скорость ударника, м/с Контактная сила, Н Максимальное ускорение Внутричерепное давление, КПА HIC
0 2,5 1438 45g 21 25
1 5 3181 99g 48 161
2 6,7 4442 139g 67 347
3 8,3 5606 175g 87 548
4 9,98 6882 215g 107 965
5 12,5 8906 278g 138 1746
6 15 10925 339g 170 2793
При значениях силы удара в диапазоне 4000-6000 Н можно говорить о получении травмы средней и серьезной тяжести, а величины выше 7000 Н расцениваются как удары, влекущие тяжелые последствия с длительной потерей сознания, при силе, равной 10 000 Н, с большой долей вероятности наступит смерть потерпевшего.
Согласно данным из работы [39], о возникновении черепно-мозговой травмы в рассматриваемом случае можно говорить при значении внутричерепного давления, равном 67 КПа, что соответствует травме средней тяжести.
Заключение
В настоящей работе проведена валидация трехмерной конечно-элементной модели головы человека, построенной из объемных данных, и вычислительного алгоритма, включающего решения контактной задачи методом декомпозиции Шварца и неявной обобщенной диссипативной схемы Ньюмарка с перераспределением масс на границе контакта.
Установлено, что распределение внутричерепного давления в лобной доле на 25% ниже экспериментального для отрицательных значений внутричерепного давления в задней части головы, как и в случае моделей на основе контактного взаимодействия между мозгом и черепом, наблюдаются завышенные по модулю величины.
Проведено сравнение значений H/C-критерия с данными линейного ускорения и внутричерепного давления при различных значениях приложенной контактной силы. Установлено, что значения силы удара более 7000 Н влекут тяжелые последствия для человека.
В целом полученные результаты можно признать удовлетворительными, что говорит о возможности применения разработанной модели в различных задачах контактной биомеханики для случаев кратковременных ударных нагрузок.
Благодарности
Работа выполнена при финансовой поддержке РФФИ грант № 17-01-00402.
Список литературы
1. Агапов П.И., Белоцерковский О.М., Петров И.Б. Численное моделирование последствий механического воздействия на мозг человека // Вычислительная матемематика и математическая физика. - 2006. - Т. 46, № 9. - С. 1711-1720.
2. Белошицкий В.В. Принципы моделирования черепно-мозговой травмы в эксперименте // Украинский нейрохирургический журнал. - 2008. - № 4. - С. 9-15
3. Васюков А.В., Петров И.Б. Моделирование механических факторов черепно-мозговых травм сеточно-характеристическим численным методом // Вестник Российского государственного университета им. И. Канта. - 2010. - Вып. 10. - С. 42-51.
4. Громов А.П. Биомеханика травмы. - М.: Медицина. - 1979. - 270 с.
5. Караваев А.С., Копысов С.П. Метод построения неструктурированных шестигранных сеток из объемных данных // Компьютерные исследования и моделирование. - 2013. - Т. 5. - С. 11-24.
6. Караваев А.С., Копысов С.П. Метод Шварца декомпозиции области в контактных задачах // Проблемы механики и материаловедения: тр. ИМ УрО РАН. - 2017. - С. 6-20.
7. Караваев А.С., Копысов С.П., Кузьмин И.М. Построение биомеханических конечно-элементных моделей // Проблемы механики и материаловедения: тр. ИМ УрО РАН. - 2016. - С. 7-13.
8. Корсаков С.А. Биомеханика ударного взаимодействия травмирующего предмета и человека // Медицинская экспертиза и право. - 2012. - Т. 12, № 3. - С. 10-12.
9. Недугов Г.В. Определение и экспертная оценка силы удара // Проблемы экспертизы в медицине. -2012. - Т. 12, вып. 45-46-1-2. - С. 14-17.
10. Bandak F.A., Vander Vorst M.J., Stuhmiller L.M., Mlakar P.F., Chilton W.E., Stuhmiller J.H. An imaging-based computational and experimental study of skull fracture: Finite element model development // Journal of Neurotrauma. - 1995. - Vol. 12, № 4. - P. 679-688.
11. Belingardi G., Chiandussi G., Gaviglio I. Development and validation of a new finite element model of human head // Proceedings of 19th International Technical Conference on the Enhanced Safety of Vehicles. - Washington, 2005. - P. 1-9.
12. Claessens M., Sauren F., Wismans J. Modelling of the human head under impact conditions: A parametric study // Transactions of SAE. - 1997, № 973338. - P. 3829-3848.
13. Deck C., Willinger R. Improved head injury criteria based on head FE model // International Journal of Crashworthiness. - 2008. - Vol. 13, № 6. - P. 667-678.
14. Hardy W.N., Foster C., Mason M., Yang K., King A., Tashman S. Investigation of head injury mechanisms using neutral density technology and high-speed biplanar X-ray // Proceedings of 45th Stapp Car Crash Conference. - San Antonio, 2001. - P. 337-368.
15. Head injury criterion [Электронный ресурс] // Wikipedia. - URL: https://en.wikipedia.org/wiki/ Head_injury_criterion (дата обращения: 02.02.2018).
16. Henn H. Crash tests and the head injury criterion // Teaching Mathematics and its Applications. - 1998. -Vol. 17, № 4.
17. Hilber H.M., Hughes T.J.R., Taylor R.L. Improved numerical dissipation for time integration algorithms in structural dynamics // Earthquake Engrg. Struct. Dynam. - 1977. - Vol. 5. - P. 283-292.
18. Horgan T. J., Gilchrist M. The creation of three dimensional finite element models for simulating head impact biomechanics // International Journal of Crashworthiness. - 2003. - Vol. 8, № 4. - P. 353-366.
19. Huang, H.M. Finite element analysis of brain contusion: an indirect impact study // Medical Biological Engineering Computing. - 2000. - Vol. 38, № 3. - P. 253-259.
20. Khalil T.B., Hubbard R.P. Parametric study of head response by finite element modeling // Journal of Biomechanics. - 1977. - Vol. 10, № 2. - P. 119-132.
21. Khenous H.B., Laborde P., Renard Y. Mass redistribution method for finite element contact problems in elastodynamics // European Journal of Mechanics A/Solids. - 2008. - Vol. 27, iss. 5. - P. 918-932.
22. King A.I., Yang K.H., Zhang L., Hardy W, Viano D. Is head injury caused by linear or angular acceleration // Proceedings of International IRCOBI Conference on the Biomechanics of Impacts. - Lisbon, 2003. - P. 1-12.
23. Kleiven S. Finite element modeling of the human head // Department of Aeronautics Royal Institute of Technology. - Stockholm, 2002.
24. Kleiven S., Hardy W.N. Correlation of an FE model of the human head with local brain motion-consequences for injury prediction // Proceedings of 46th Stapp Car Crash Conference. - Ponte Vedra, 2003, № 2002-22-0007. - P. 123-144.
25. Kristin L.B., Ahmad R.H., Daniel R.W. Neuroimaging genetics: principles and practices // Oxford University Press, 2015. - P. 157.
26. Kwong M.T. A review of head injury and finite element head models // American Journal of Engineering. -2014. - Vol. 1, № 5. - P. 28-52.
27. Marjoux D., Baumgartner D., Deck C., Willinger R. Head injury prediction capability of the HIC, HIP, SIMon and ULP criteria -New injury criteria for the head // Accident Analysis and Prevention. - 2008. -Vol. 40, № 3. - P. 1135-1148.
28. Nahum A.M., Smith R., Ward C.C. Intracranial pressure dynamics during head impact // Proceedings of 21st Stapp Car Crash Conference. - San Diego, 1977. - P. 339-366.
29. Paolo D.L. Adjustments to zatsiorsky-seluyanovs segment inertia parameters // Journal of Biomechanics. -1996. - Vol. 29, № 9. - P. 1223-1230.
30. Patel A., Goswami T. Comparison of intracranial pressure by lateral and frontal impacts - validation of computational model engineering biomedical engineering injury and skeletal biomechanics. - IntechOPEN, London, 2012.
31. Rezaei A., Karami G., Ziejewski M. Examination of brain injury thresholds in terms of the severity of head motion and the brain stresses // International brain injury association. - URL: http://www.internationalbrain.org/examination-of-bi-thresholds-in-terms-of-the-severity-of-head-motion-and-the-brain-stresses/ (дата обращения: 02.02.2018).
32. Ruan J.S., Khalil T., King A.I. Dynamic response of the human head to impact by three-dimensional finite element analysis // Journal of Biomechanical Engineering. - 1994. - Vol. 116, № 1. - P. 44-50.
33. Ruan J.S., Khatil T.B, King A.I. Finite element modeling of direct head impact // In Proceedings of 37th Stapp Car Crash Conference. - San Antonio, USA, 1993. - № 933114. - P. 69-81.
34. Trosseille X., Tarriere C, Lavaste F. Development of a FEM of the human head according to a specific test protocol // Proceedings of 30th Stapp Car Crash Conference. - Warrendale, 1992. - P. 235-253.
35. Willinger R., Baumgartner D. Numerical modeling of the human head under impact: new injury mechanisms and tolerance limits // IUTAM Symposium on Impact Biomechanics: From Fundamental Insights to Applications. - Netherlands, Springer, 2005. - P. 195-203.
36. Willinger R., Taleb L., Pradoura P. From the finite element model to the physical model // Proceedings of the International IRCOBI Conference on the Biomechanics of Impacts. - 1995. - P. 245-259.
37. Willinger R., Taleb L., Kopp C.M. Modal and temporal analysis of head mathematical models // Journal of Neurotrauma. - 1995. - Vol. 12, № 4. - P. 743-754.
38. Yang B. Development of a finite element head model for the study of impact head injury // Hindawi Publishing Corporation BioMed Research International. - 2014. - Vol. 2014.
39. Zhang L., Yang K.H., King A.I. A proposed injury threshold for mild traumatic brain injury // Journal of Biomechanical Engineering. - 2004. - Vol. 126, № 2. - P. 226-234.
40. Zhang L., Yang K.H., Dwarampudi R., Omori K., Li T., Chang K., Hardy W.N., Khalil T.B, King A.I. Recent advances in brain injury research: a new human head model development and validation // Proceedings of 45th Stapp Car Crash Conference. - San Antonio, 2001. - № 2001-22-0017. - P. 369-394.
41. Zhou C., Khalil T.B., King A.I. A new model comparing impact responses of the homogeneous and inhomogeneous human brain // Proceedings 39th Stapp Car Crash Conference. - 1995. - № 952714.
42. Zhou C., Khalil T.B., King A.I. Viscoelastic response of the human brain to sagittal and lateral rotational acceleration by finite element analysis // Proceedings of the International IRCOBI Conference on the Biomechanics of Impacts. - 1996. - P. 35-48.
43. Zong Z., Lee H., Lu C. A three-dimensional human head finite element model and power flow in a human head subject to impact loading // Journal of Biomech. - 2006. - Vol. 39, № 2. - P. 284-292.
MATHEMATICAL MODELLING OF HEAD IMPACT WITH CRANIOCEREBRAL INJURY
A.S. Karavaev, S.P. Kopusov (Izhevsk, Russia)
In the present paper, we review known mathematical models of a human head, which are applied in the contact problems of biomechanics. The finite element model was constructed from voxel data. This model involves of three main materials: soft tissue, skull and brain, defined on the basis of computed tomography with using the algorithm developed by the authors. In the paper, the iterative method of Schwarz domain decomposition was used
for two contacting 3D bodies which allowed us to consider the solving of the stress-strain state problem in the usual formulation for the each contacting body. According to the offered iteration algorithm, two steps are performed in each iteration, where the conjugation conditions for the displacements and stresses are fulfilled in turn for the nodes on the contact surface. The validation of the developed finite element model has been conducted for the case of a short duration impact. A number of computational experiments have been carried out for establishing critical values of a contact force leading to severe injuries. At the impact force values in the range of 4000-6000 N, moderate and serious injuries can occur; the impact force values above 7000 N are considered as resulting in severe injuries with loss of consciousness for a long period; the force of 10000 N most probably can cause the death of an injured.
Key words: mathematical modelling, contact interaction, biomechanics, craniocerebral injury, finite element method.
Получено 5 февраля 2018
