Конечно-элементное моделирование деформирования костной ткани на субмикроскопическом уровне Текст научной статьи по специальности «Физика»

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

УДК 539.32

А.С. Семенов, А.И. Грищенко, Б.Е. Мельников

Санкт-Петербургский государственный политехнический университет

КОНЕЧНО-ЭЛЕМЕНТНОЕ МОДЕЛИРОВАНИЕ ДЕФОРМИРОВАНИЯ КОСТНОЙ ТКАНИ НА СУБМИКРОСКОПИЧЕСКОМ УРОВНЕ

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

КОСТНАЯ ТКАНЬ, ЭЛЕМЕНТАРНЫЙ ПРЕДСТАВИТЕЛЬНЫЙ ОБЪЕМ, ГОМОГЕНИЗАЦИЯ, ГРАНИЧНОЕ УСЛОВИЕ ПЕРИОДИЧНОСТИ, РАЗОРИЕНТАЦИЯ КОНГЛОМЕРАТОВ, МЕЖКРИСТАЛЛИТНАЯ ПЕРЕМЫЧКА, УПРУГОСТЬ, ВЯЗКОУПРУГОСТЬ, МЕТОД КОНЕЧНЫХ ЭЛЕМЕНТОВ.

Введение

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

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

Под влиянием различных факторов, таких как оперативное вмешательство, хро-

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

Морфологическая модель костного матрикса

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

■ИЯВР^1

ЖЭШЙ

Мм, в

% 1) ЩИ

— 5

си:

— 6

Рис. 1 [2]. Схема распределения минеральных элементов в костном матриксе:

1 — коллагеновые фибриллы; 2 — кристаллические объединения на удалении от коллагеновых фибрилл; 3 — фибриллярная минеральная манжетка; 4 — кристаллические объединения коллагеновой фибриллы; 5 — молекулы коллагена; 6 — пласты кристаллов (их копланарные объединения)

вым фибриллам, формируют манжетки и контактируют с внутрифибриллярными объединениями. На отдалении от коллагеновых фибрилл кристаллиты образуют конгломераты, в которых они ориентированы в одном направлении. По отношению друг к другу эти конгломераты расположены под разными углами. Внутрифибрилляр-ные объединения кристаллитов образуют ряд параллельных, спирально закрученных пластов, которые ориентированы под углами 8—25° к осям фибрилл. Отдельные объединения минералов связывают соседние пласты, обеспечивая непрерывность минерального компонента фибрилл после удаления органической составляющей [3]. Наличие минеральных связей между элементами минерального матрикса подтверждено электронно-микроскопически; проведен также детальный анализ влияния минеральной связи между объединениями кристаллитов на механические свойства костного матрикса [3].

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

Элементарный представительный объем нанокомпозита кости

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

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

вающие наличие связи между кристаллами гидроксиапатита в виде твердотельных перемычек (рис. 2).

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

60 до 90 % (см. рис. 2). Изменение минерализации моделировалось пропорциональным изменением размеров конгломератов при сохранении габаритов ЭПО (фиксация местоположений центров конгломератов). Наряду с базовым вариантом рассматривались также и другие альтернативные сценарии роста минерализации, связанные с непропорциональным изменением размеров конгломератов.

Рассматривались также альтернативные варианты ЭПО (рис. 3), которые включали

Рис. 3. КЭ-модели элементарного представительного объема наноструктурных элементов костной ткани с различным числом конгломератов (к) и перемычек (п): а - модель П 0,5 х 0,5 (к = 0,5, п = 1, 304 КЭ); б - П 1 х 1 (к = 2, п = 4, 1216 КЭ); в - П 2 х 2 (к = 8, п = 16, 4864 КЭ); г - П 4 х 4 (к = 32, п = 64, 19456 КЭ); д - П 8 х 8 (к = 128, п = 256, 77824 КЭ)

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

Краевые задачи теории упругости для определения напряженно-деформированного состояния ЭПО решались методом конечных элементов. Последний основан на уравнении виртуальных работ, которое при консервативном внешнем воздействии приводит к минимизации потенциальной энергии системы. Введенные модели реализованы в КЭ программном комплексе РЛКТОСКЛТОК [4]; он позволяет автоматически генерировать дискретные модели ЭПО произвольной геометрии, получать решения краевых задач теории упругости, определять эффективные упругие модули, исследовать процессы разрушения ЭПО, анализировать распределения полей напряжений и деформаций для фрагмента морфологической структуры матрикса в наномасштабе.

Минерализация элементарного представительного объема определяется отношением плотностей компонент и геометрическими параметрами конгломератов и перемычек:

т =

1

1 + Г 2( а + 1)(Ь + ^ 11

аЬ + 1й Ра

(1)

где а, Ь — ширина и высота конгломерата; I, й — ширина и высота перемычки; g — вертикальный зазор между конгломератами; рс , рА — плотности коллагена и гидрокси-апатита соответственно.

В расчетах при задании различных уровней минерализации размеры конгломерата (объединения кристаллитов) варьировались в диапазонах: ширина а = 200 — 500 нм, высота Ь = 1800 — 3000 нм. Важно отметить, что указанные габариты минералов соответствуют размерам калькосферитов, описанных А. Бойдом (Воуёе) [5 — 6] в формирующихся участках кости. У взрослого человека в минерализующемся костном матриксе расстояния между фибриллами не превышают 20 — 30 нм, поэтому морфологический субстрат минерала включает в себя фрагменты минерализованных колла-геновых фибрилл (их диаметр равен 50 — 80 нм, длина — до 2000 нм) с окружающими их «манжетками», единичные матриксные везикулы (диаметр составляет 30 — 200 нм) и разнонаправленные плотные группы кристаллов, связывающие манжетки соседних фибрилл.

В расчетах принималось, что механические свойства отдельных компонент модели соответствуют изотропным материалам. Расчетные значения упругих модулей (табл. 1) выбирались как средние арифметические экспериментальных величин, приведенных в работах [7 — 10].

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

Формулировка адекватных моделей, которые бы позволяли описывать физико-механические свойства объектов, имею-

Таблица 1

Значения упругих модулей и плотностей отдельных компонент ЭПО, использованных в КЭ-расчетах

Параметр Обозначение Единица измерения Значение

Гидроксиапатит Коллаген Перемычка

Модуль Юнга Е ГПа 90 0,9 90

Коэффициент Пуассона V - 0,20 0,49 0,20

Плотность Р кг/м3 3160 1500 3160

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

В работе [11] показано, что соответствие механических характеристик нано-размерных объектов данным, полученным из макроскопических экспериментов, достигается лишь при количестве атомарных слоев не менее десяти; только тогда разница между микро- и макрохарактеристиками не превышает 11%.

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

Таким образом, для дальнейшего исследования свойств ЭПО нанокомпозита костной ткани необходимо определить порядок масштабных эффектов. Минимальный характерный размер элементов ЭПО составляет 30 нм — это ширина минеральной связи. Далее приведены размеры атомарных слоев основных веществ ЭПО [12, 13].

Вещество

Кристалл гидроксиапатита Молекула коллагена

Размер атомарного слоя,А

а = 9,41, с = 6,88 0 = 15

Таким образом, ширина минеральной связи, которая является наименьшим размером в рассматриваемых ЭПО, содержит в себе как минимум 30 атомарных слоев, что, согласно работе [11], позволяет пренебречь масштабными факторами. Однако необходимо учитывать, что ширина зазора между конгломератами не должна быть меньше 15 А -10 = 150 А = 15 нм, что накладывает дополнительные, но несущественные ограничения на исследуемый диапазон изменения минерализации.

Определение эффективных упругих свойств на основе КЭ-расчетов

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

£ = ^ • • С, (2)

где £ — тензор деформаций, с — тензор напряжений, 4 C — тензор упругих податли-востей 4-го ранга; двумя точками обозначена операция свертки (двойного скалярного умножения), которая в индексной форме имеет вид

8>у = а1к.

Черта над обозначениями введенных тензоров означает соответствие гомогенному (осредненному) материалу. Характеристики отдельных компонент гетерогенного материала будут обозначаться без верхней черты. Для ортотропного материала тензору 4 C в собственных осях анизотропии соответствует симметричная матрица упругих податливостей следующей структуры:

И =

-4Л -431 0 0

-Д32 0 0

1 V31

Д Е2 Е3

1 *32

Е1 Ё2 Е3

*23 1

Е1 Е2 Е3

0 0 0

0

0 0 0

1 О12 0 0

0 1 О23 0

0 0 1 О31

0 0

,(3)

где Е1, v¡J, Оу — модули Юнга, коэффициенты Пуассона и модули сдвига (/, у = 1,2,3).

В силу симметрии матрицы выполняются равенства

12 Е

г

Jl ^2 ^2 3

Упругие модули определяются при од-

у_31 __23_

у_32

ЕЕ

ноосном напряженном состоянии на основе соотношений

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

При решении задачи в двумерной постановке (плоское деформированное состояние) достаточно найти четыре упругие константы: Е1,Е2,у12,^12. Ряд упругих модулей в направлении оси 3 (ортогональна рассматриваемой плоскости) можно оценить без потери точности по правилу смесей. Разносопротивляемость костной ткани при растяжении и сжатии не учитывалась в рамках данной модели.

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

Кинематические граничные условия:

и = ¿'•г. (4)

Статические граничные условия:

п • о= п • 0*, (5)

Условия периодичности:

= 4, +«• (1*1-Г2). (6)

Чщ и

Здесь и — вектор перемещений; г — радиус-вектор; г* — заданный постоянный симметричный тензор, соответствующий различным деформируемым состояниям (осевые растяжения/сжатия и сдвиги); о* — заданный постоянный симметричный тензор, соответствующий различным напряженным состояниям (осевые растяжения/сжатия и сдвиги).

Условие периодичности (6) может быть переписано в виде

ul = ** • Г + w,

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

В этом случае вектор напряжений n • о принимает противоположные значения на противоположных сторонах ЭПО. Граничные условия (4) — (6) удовлетворяют условию макрооднородности Хилла

О •• Ё = О •• Ё,

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

Для определения компонент матрицы эффективных упругих податливостей [cj

(3) использовались осредненные по представительному объему значения тензоров деформаций и напряжений, полученных в КЭ-решении:

Ё = у-I* , = Ь dV (7)

ЭПО Уэпо ЭПО Кэпо

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

Оценка точности решений

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

(4) и напряжений (5) (рис. 4). Установлено,

что если задавать условия периодичности (6), то удовлетворительная (статистически существенная) точность, т. е. достаточная близость к асимптотическому (предельному) значению, достигается уже при использовании ЭПО, включающих одну ячейку периодичности (КЭ-модель П 1 * 1). Если же ставить граничные условия для перемещений и напряжений, то сходимость достигается только при использовании ЭПО, включающих 16 ячеек периодичности (П 4 * 4). Данный вывод подтверждается также анализом зависимости экстремальных значений напряжений в ЭПО.

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

Таким образом, корректно выбранные тип граничных условий и число ячеек периодичности, достаточное для обеспечения заданной точности, а также обоснованный выбор числа КЭ, обеспечивают проведение информативных вычислительных экспериментов с виртуальным варьированием

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

Анализ напряженно-деформированного состояния ЭПО при различных видах нагружения

Наличие в ЭПО гетерогенной микроструктуры приводит к реализации неоднородного многоосного напряженного состояния. Для сравнения уровней напряженного состояния для различных моделей и условий нагружения используется эквивалентная величина — интенсивность напряжений по Мизесу:

а,- = {(1/2)[( а^ -а )2 + (а - а )2 +

УУ

+(а -а )2] + 3[т2 +т2 +т2 ]}

V 77 ХХ/Л 1-ЛУ У1 7Х !>

I1/2

(8)

Она представляет собой эвклидову норму в пространстве компонент девиатора тензора напряжений.

Анализ распределения полей интенсивности напряжений по Мизесу о. при растяжении ЭПО в вертикальном (рис. 5, а, б) и в горизонтальном (рис. 5, в, г) направлениях показал, что в условиях вертикальной

Рис. 4. Зависимости эффективных значений модулей упругости в вертикальном (а) и в горизонтальном (б) направлениях от количества ячеек периодичности в ЭПО при различных граничных условиях. Каждое значение количества ячеек соответствует определенной КЭ-модели: П 0,25 * 0,25 (I), П 1 * 1 (II), П 2 * 2 (III), П 3 * 3 (IV), П 4 * 4 (V), П 5 * 5 (VI), П 8 * 8 (VII). Граничные условия: по формуле (4) — кривая 1; (5) — 3; (6) — 2

а) б) в) г)

г т У1

i 1 _А 1

Рис. 5. Распределения интенсивностей напряжений а. (а, в) и деформаций е. (б, г) по Мизесу в центральном фрагменте ЭПО (1/4 часть) при различных направлениях воздействия нагрузки: вертикальном (а, б, д) и горизонтальном (в, г, е). На графиках д, е дополнительно представлено деформированное состояние перемычки с полями

интенсивности напряжений

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

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

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

Влияние различных факторов на механические свойства ЭПО костной ткани

Наличие перемычек. Результаты КЭ-расчетов для модели с перемычками сравнивались с аналогичными результатами, полученными при отсутствии перемычек

Т аб л и ц а 2

Сравнение расчетных значений эффективных упругих модулей для ЭПО с перемычками и без перемычек

Элементарный представительный объем (ЭПО) Модуль упругости, ГПа

Е2 Е1

С перемычками 18,28 15,71

Без перемычек 16,63 13,32

Различие ЭПО 1,65 (9,9%) 2,39 (17,9%)

Обозначения: Е1, Е2 — модули упругости в горизонтальном и вертикальном направлениях нагрузки, соответственно.

(прочие условия одинаковы). Установлено, что наличие перемычек при растяжении в горизонтальном направлении повышает уровень экстремальных напряжений (относительная интенсивность напряжений ст(. / ст ,) на 52 %. При этом уровень интенсивности напряжений в угловой точке минеральной составляющей снижается на 19 %.

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

Перемычки увеличивают модуль упругости Е1 (горизонтальное направление) на 17 % (табл. 2) и модуль упругости Е2 (вертикальное направление) на 10 %. Следует отметить, что наблюдаемое увеличение жесткости вызвано изменением структуры (топологии) ЭПО, вносимым появлением перемычек, а не повышением минерализации, влияние которой на упругие модули на два порядка меньше.

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

Углы наклона минеральной составляющей. Чтобы оценить влияние углов разори-ентации минеральной компоненты в ЭПО на ее напряженно-деформированное состояние и эффективные упругие свойства, нами проведены КЭ-расчеты с тремя различными значениями указанного угла: 0о, 4о и 8о. Результаты проведенных вычислительных экспериментов показали, что отклонение единичного конгломерата от вертикали фактически не оказывает влияния назначения эффективных упругих модулей Е1, Е2 (отличия составляют менее 1 %), но такое отклонение влияет на уровень максимальной интенсивности напряжений (рис. 6).

тах(а,. / ст,.)

2 ч

у

Г

1

012345678 Угол поворота конгломерата, град

Рис. 6. Зависимости максимальной интенсивности напряжений от угла поворота центрального конгломерата в фрагменте ЭПО при различных направлениях нагружения: по горизонтали (1) и по вертикали (2)

1,2 -

О о

5 <—

!кГ 1,1 -

о о

Р 1.0-

0,9-

1 \ ,

-40 -20 0 20

Угол поворота перемычек, град

40

Рис. 7. Зависимости эффективных упругих модулей ЭПО нанокомпозита костной ткани от угла поворота перемычек; 1, 2 — величины

Еу/Еу

угол 0°

и Е2/Е2

2 угол 0°

соответственно

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

Угол наклона перемычек. Чтобы оценить влияние угла наклона перемычек на напряженно-деформированное состояние и эффективные упругие свойства ЭПО, нами была проведена серия КЭ-расчетов с различными углами отклонения от горизонтального уровня: на —45°, —30°, —15°, +15°, +30° и +45°. Результаты расчетов показали, что максимальное отклонение Е2 от значения при горизонтальном положении перемычки не превышает 23 %, а такое же отклонение Е1 не превышает 13 % (рис. 7). Установлено также, что поворот перемычек сказывается на степени максимальной интенсивности напряжений.

Согласно полученным результатам расчета, поворот перемычек в ЭПО при одних условиях нагружения увеличивает уровень экстремальных напряжений в ней, а при других — уменьшает. Так, при отклонении перемычки на 15° в направлении против

часовой стрелки растяжение в горизонтальном направлении характеризуется снижением максимальной относительной интенсивности напряжений на 2%, при 30° — на 12 %, при 45° — на 47 %, соответственно, по сравнению с идеализированным вариантом (без поворота). Этот эффект вызван тем, что при растяжении в горизонтальном направлении ЭПО сжимается в вертикальном, и поэтому соседние конгломераты смещаются по вертикали друг относительно друга; это приводит к повороту соединяющей их перемычки и возникновению в ней дополнительных усилий. В зависимости от направления начального отклонения перемычки от горизонтали вертикальные смещения конгломератов могут создать в них как дополнительное сжимающее напряжение (для перемычек, у которых поворот под нагрузкой компенсирует начальный поворот), так и дополнительное растягивающее напряжение (для перемычек, у которых поворот под нагрузкой увеличивает начальное отклонение). Появление дополнительных напряжений может как компенсировать, так и усиливать первоначальное напряжение, возникшее от растяжения в горизонтальном направлении. Данный сценарий реализуется и при нагружении в вертикальном направлении, для которого отклонение перемычки на 15° в направлении против часовой стрелки приводит к увеличению максимальной относительной интенсивности напряжений на 65 %, при 30о — на 80 %, при 45о — на 47 %, по сравнению с идеализированным вариантом (без поворота).

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

Минерализация костной ткани. Чтобы оценить влияние минерализации ЭПО на его напряженно-деформированное состояние и эффективные упругие свойства, нами проведены КЭ-расчеты с различной минерализацией костной ткани (%): 25, 30, 45, 50, 60, 70, 80, 90 и 95. Результаты проведенных

£„£„ГПа

30-

70

60

20-

у

х' У

У

\ ' <

У >* "V

у

л' х I

(¡^ 4

0,2

0.3

0,4

0.5

0,6

0,7

0.3

0,0

т

Рис. 8. Зависимости эффективных упругих модулей ЭПО нанокомпозита костной ткани от степени его минерализации.

1, 2 — Е1, Е2 соответственно (данные получены методом КЭ-гомогенизации); 3, 4 — аналитические оценки Фойгта и Рейса (осредненные жесткости и податливости компонент соответственно)

вычислительных экспериментов показали, что изменение минерализации существенно влияет на эффективные модули (рис. 8) и на осредненную по объему интенсивность напряжений в ЭПО и в его составляющих. Так, при изменении минерализации с 60 до 70 %, Е2 (в вертикальном направлении) из-

меняется на 55,7 %, а Е1 (в горизонтальном направлении) — на 47,7 %. Осредненная же по всему объему интенсивность напряжения, при том же изменении минерализации, возрастает в ЭПО на 54 %.

Степень вытянутости конгломератов. Важным морфологическим параметром является соотношение ширины и высоты конгломератов. Чтобы выяснить влияние указанного параметра на величины эффективных упругих модулей, была поставлена серия вычислительных экспериментов с пропорцией Ь/а, меняющейся в пределах от 5 до 10. КЭ-модели для отношений Ь/а, равных 5, 7 и 9, представлены на рис. 9. Установлено, что свойства ЭПО костной ткани чрезвычайно сильно зависят от степени вытянутости конгломератов. Так, при увеличении отношения Ь/а от 5 до 10 величина Е2 возрастает почти в три раза, а Е1 — в полтора (рис. 10).

Вязкие свойства органического матрик-са. Для анализа влияния вязких свойств органического матрикса рассматривалась ползучесть ЭПО в течение одного часа при напряжении 10 МПа (12-кратная весовая перегрузка). Результаты расчетов показали, что вязкие эффекты крайне незначительны (рис. 11). Так, интенсивность деформаций ползучести, осредненная по ЭПО, не пре-

Рис. 9. КЭ-модели ЭПО, полученные при разных отношениях Ь/а: 5 (а), 7 (б); 9 (в); а, Ь — ширина и высота конгломерата соответственно

Е., Е.. ГПа

40-

35 -

30-

25-

20-

10

2

у'

—1—

Ь/а

Рис. 10. Зависимости эффективных упругих

модулей ЭПО в горизонтальном (1) и вертикальном (2) направлениях от степени вытянутости конгломератов; 1, 2 — Е1, Е2 соответственно

Рис. 11. Зависимости осредненной по ЭПО интенсивности деформаций ползучести от времени для случаев наличия (1) и отсутствия (2) перемычек

восходит 0,00014 %. Также видно, что наличие перемычек практически не оказывает влияния на осредненную интенсивность деформаций ползучести.

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

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

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

Заключение

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

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

Таблица 3

Сравнение расчетных значений эффективных упругих модулей ЭПО костной ткани

с экспериментальными данными

КЭ-модель ЭПО либо эксперимент Модуль упругости, ГПа

e2 E

Модель с перемычками, m = 0,7 18,3 15,7

Модель без перемычек, m = 0,694 16,6 13,3

Эксперимент Bonfield W. [14] 18,5 9,5

Эксперимент Ashman R.B. [15] 20,0 13,5

Эксперимент Turner C.H. [16] 20,6 16,5

Обозначение: m — степень минерализации костной ткани.

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

Проведены многовариантные вычислительные эксперименты с варьированием минерализации костной ткани. Выполнен анализ влияния минерализации ЭПО на напряженное состояние представительного объема и его эффективные модули.

СПИСОК ЛИТЕРАТУРЫ

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

Работа выполнена при финансовой поддержке Министерства образования и науки РФ.

1. Денисов-Никольский Ю.И., Миронов С.П., Омельяненко Н.П., Матвейчук И.В. Актуальные проблемы теоретической и клинической остео-артрологии. М.: ОАО «Типография "Новости"», 2005. 336 с.

2. Докторов А.А., Денисов-Никольский Ю.И., Жилкин Б.А. Структурная организация костного минерала // Бюллетень экспериментальной биологии и медицины. 1996. Т. 122. № 12. С. 687-691.

3. Аврунин А.С., Семёнов А.С., Фёдоров И.В., Мельников Б.Е., Докторов А.А., Паршин Л.К. Влияние минеральной связи между объединениями кристаллитов на механические свойства костного матрикса. Моделирование методом конечных элементов // Травматология и ортопедия России. 2013. № 2 (67). С. 1-12.

4. Семёнов А.С. PANTOCRATOR - конечно-элементный программный комплекс, ориентированный на решение нелинейных задач механики // В кн.: Научно-технические проблемы прогнозирования надежности и долговечности конструкций и методы их решения. Труды V Международной конференции. СПб.: Изд-во Политехнического ун-та, 2003. С. 466-480.

5. Boyde A., Hobdell M.H. Scanning electron-microscopy of lamellar bone, Z. Zellforsch. 1969. No. 93. S. 213-231.

6. Boyde A. Scanning electron microscope studies of bone. The biochemistry and physiology of bone. Ed. G.H. Bourne. N.Y., London: Academic Press, 1972. No. 1, pp. 259-310.

7. Lawson А.^, Czernuszka J.T. Collagen-calcium phosphate composites, Proc. Instn. Mech. En-

grs. 1998. Vol. 212. Part H, pp. 413-425.

8. Gilmore R.S, Katz J.L. Elastic properties of apatites. Journal of Materials Science. 1982. Vol. 17. No. 4, pp. 1131-1141.

9. Katz J.L., Ukraincik K. On the anisotropic elastic properties of hydroxyapatite. J. Biomechanics. 1971. Vol. 4, pp. 221-227.

10. Cowin S.C. Bone mechanics. Boca Raton, CRC Press, FL, 1989.

11. Кривцов A.M., Морозов Н.Ф. О механических характеристиках наноразмерных объектов // Физика твердого тела. 2002. Т. 44. № 12. С. 2158 -2163.

12. Bhatnagar V.M. Refinement of the synthetic hydroxyapatite cell parameters, Contributions to Mineralogy and Petrology. 1969. Vol. 22. Iss. 4, pp. 375-378.

13. Barkaoui A., Hambli R. Finite element 3D modeling of mechanical behavior of mineralized collagen microfibrils. J. Appl. Biomater. Biomech. 2011. Sep-Dec. Vol. 9. No. 3, pp. 199-205.

14. Bonfield W., Grynpas M.D. Anisotropy of the Young's modulus of bone, Nature. 1977. Vol. 270. No. 5636, pp. 453-454.

15. Ashman R.B., Cowin S.C., Van Bus-kirk W.C., Rice J.C. A continuous wave technique for the measurement of the elastic properties of cortical bone, J. Biomechanics. 1984. Vol. 5, pp. 349-361.

16. Turner C.H., Rho J.-Y., Takano Y., Tsui T.Y., Pharr G.M. The elastic properties of trabecular and cortical bone tissue are similar: Results from two microscopic measurement techniques. Journal of Biomechanics. 1999. Vol. 32, pp. 437-441.

СВЕДЕНИЯ ОБ АВТОРАХ

СЕМЕНОВ Артём Семенович — кандидат физико-математических наук, доцент кафедры «Механика и процессы управления» Института прикладной математики и механики Санкт-Петербургского государственного политехнического университета.

195251, Россия, г. Санкт-Петербург, Политехническая ул., 29 Semenov.Artem@googlemail.com

ГРИЩЕНКО Алексей Иванович — студент Института прикладной математики и механики Санкт-Петербургского государственного политехнического университета. 195251, Россия, г. Санкт-Петербург, Политехническая ул., 29 gai-gr@yandex.ru

МЕЛЬНИКОВ Борис Евгеньевич — доктор технических. наук, заведующий кафедрой сопротивления материалов Санкт-Петербургского государственного политехнического университета. 195251, Россия, г. Санкт-Петербург, Политехническая ул., 29 kafedra@ksm.spbstu.ru

Semenov A.S., Grishchenko A.I., Melnikov B.E. FINITE ELEMENT MODELING OF BONE DEFORMATION AT THE SUBMICROSCOPIC SCALE.

By means of the direct finite element simulation and homogenization the analysis of variation influence in the morphological characteristics (hydroxyapatite crystals disorientation, sizes and orientation of mineral bridges, mineralization) on mechanical properties of the representative volume element of bone at the nanoscale (at the collagen fibrils level) is carried out. The morphological model of bone with an account of the mineral bridges between the associations of the hydroxyapatite crystals is used in the computations. The purpose of the paper is to analyze the influence of the nanostructure parameters of bone on its elastic and strength properties. Such studies are important for the creation of artificial bone-substitute materials. The analysis of the stress-strain state of the RVE of bone tissue has been performed in order to determine the location of the most critical points and deformation mechanisms of bridges. The most loaded elements are the corners of conglomerates and corners of bridges. Taking account of the bridges leads to the decrease of the von Mises stresses in the corner of the conglomerate and to the increase of the effective elastic moduli in the vertical and horizontal directions. The effects of orientation of conglomerates and bridges on the stress state of the representative volume were analyzed. The influence of the bone mineralization on the effective elastic moduli and stress state were investigated. The comparison of obtained results with experimental data was performed and discussed.

BONE TISSUE, REPRESENTATIVE VOLUME ELEMENT, HOMOGENIZATION, BOUNDARY CONDITION OF PERIODICITY, MISALIGNMENT OF CONGLOMERATES, INTERCRYSTALLINE BRIDGE, ELASTICITY, VISCOELASTI-CITY, FINITE ELEMENT METHOD.

REFERENCES

1. Denisov-Nikolskiy Yu.I., Mironov S.P., Omelyanenko N.P., Matveychuk I.V. Aktualnyye problemy teoreticheskoy i klinicheskoy osteoartrologii. Moscow, OAO «Tipografiya "Novosti"», 2005, 336 s.

2. Doktorov A.A., Denisov-Nikolskiy Yu.I., Zhilkin B.A. Strukturnaya organizatsiya kostnogo minerala. Byulleten eksperimentalnoy biologii i meditsiny, 1996, T. 122, № 12, pp. 687-691.

3. Avrunin A.S., Semenov A.S., Fedorov I.V., Melnikov B.Ye., Doktorov A.A., Parshin L.K. Vliyaniye mineralnoy svyazi mezhdu obyedineniyami kristallitov na mekhanicheskiye svoystva kostnogo matriksa. Modelirovaniye metodom konechnykh elementov. Travmatologiya i ortopediya Rossii, 2013, № 2 (67), pp. 1-12.

4. Semenov A.S. PANTOCRATOR -konechno-elementnyy programmnyy kompleks, oriyentirovannyy na resheniye nelineynykh zadach mekhaniki. V kn.: Nauchno-tekhnicheskiye problemy prognozirovaniya nadezhnosti i dolgovechnosti

konstruktsiy i metody ikh resheniya: tr. V mezhdunarodnoy konf. SPb.: Izd-vo SPbGPU, 2003, pp. 466-480.

5. Boyde A., Hobdell M.H. Scanning electronmicroscopy of lamellar bone. Z. Zellforsch. 1969, N. 93, S. 213-231.

6. Boyde A. Scanning electron microscope studies of bone. In: The biochemistry and physiology of bone. Ed.:G.H. Bourne. Academic Press., N.Y. and London, 1972. No. 1, pp. 259-310.

7. Lawson A.C., Czernuszka J.T. Collagen-calcium phosphate composites. Proc. Instn. Mech. Engrs. 1998, Vol. 212, Part H, pp. 413-425.

8. Gilmore R.S, Katz J.L. Elastic properties of apatites. Journal of Materials Science. 1982, Vol. 17, No. 4, pp. 1131-1141.

9. Katz J.L., Ukraincik K. On the anisotropic elastic properties of hydroxyapatite. J. Biomechanics, 1971, Vol. 4, pp. 221-227.

10. Cowin S.C. Bone Mechanics, CRC Press, Boca Raton, FL, 1989.

11. Krivtsov A.M., Morozov N.F. O

mekhanicheskikh kharakteristikakh nanorazmernykh obyektov. Fizika tverdogo tela. 2002. T. 44. № 12, pp. 2158-2163.

12. Bhatnagar V.M. Refinement of the synthetic hydroxyapatite cell parameters. Contributions to Mineralogy and Petrology, 1969, Vol. 22, Iss. 4, pp. 375-378.

13. Barkaoui A., Hambli R. Finite element 3D modeling of mechanical behavior of mineralized collagen microfibrils. J. Appl. Biomater. Biomech. 2011, Sep-Dec; Vol. 9, No. 3, pp. 199-205.

14. Bonfield W., Grynpas M.D. Anisotropy of the Young's modulus of bone. Nature, 1977; Vol. 270, No. 5636, pp. 453-454.

15. Ashman R.B., Cowin S.C., Van Buskirk W.C., Rice J.C. A continuous wave technique for the measurement of the elastic properties of cortical bone. J. Biomechanics. 1984, Vol. 5, pp. 349-361.

16. Turner C.H., Rho J.-Y., Takano Y., Tsui T.Y., Pharr G.M. The elastic properties of trabecular and cortical bone tissue are similar: Results from two microscopic measurement techniques. Journal of Biomechanics, 1999, Vol. 32, pp. 437-441.

THE AuTHORS

SEMENOV Artem S.

St. Petersburg State Polytechnical University

29, Politekhnicheskaya St., St. Petersburg, 195251, Russia

Semenov.Artem@googlemail.com

GRISHCHENKO Alexei I.

St. Petersburg State Polytechnical University

29, Politekhnicheskaya St., St. Petersburg, 195251, Russia

gai-gr@yandex.ru

MELNIKOV Boris E.

St. Petersburg State Polytechnical University

29, Politekhnicheskaya St., St. Petersburg, 195251, Russia

kafedra@ksm.spbstu.ru

© Санкт-Петербургский государственный политехнический университет, 2014