CC BY
65
УДК 531/534:[57+61]
ПОСТАНОВКА НАЧАЛЬНО-КРАЕВОЙ ЗАДАЧИ О ПЕРЕСТРОЙКЕ ТРАБЕКУЛЯРНОЙ КОСТНОЙ ТКАНИ
А.А. Киченко, В.М. Тверье, Ю.И. Няшин, М.А. Осипенко, В.А. Лохов
Кафедра теоретической механики Пермского национального исследовательского политехнического университета, Россия, 614990, Пермь, Комсомольский проспект, 29, e-mail: nyashin@inbox.ru
Аннотация. Костная ткань является неоднородным анизотропным материалом; структурные особенности трабекулярной костной ткани могут быть описаны посредством тензора структуры. Задачи биомеханического моделирования требуют изучения истории формирования костных структур во времени как при физиологических, так и при патологических нагрузках. Это можно реализовать, имея как определяющее соотношение, связывающее тензор напряжений с тензорами структуры и деформации, так и кинетические уравнения, описывающие эволюцию тензора структуры и плотности костной ткани. В качестве таких уравнений были подробно проанализированы соотношения (S.C. Cowin) для трабекулярной костной ткани, в которых устранены некоторые неточности и недочеты. Осуществлена постановка начально-краевой задачи о перестройке трабекулярной костной ткани, решение которой позволяет проследить изменение напряженно-деформированного состояния при формировании трабекулярной структуры согласно закону Ю. Вольфа.
Ключевые слова: биомеханическое моделирование, начально-краевая задача, определяющее соотношение, эволюционное уравнение, структура костной ткани, трабекулярная (губчатая) костная ткань, тензор структуры, закон Вольфа, состояние гомеостаза (физиологическое равновесие), перестройка костной ткани.
Введение
В настоящее время существует множество неоднородных материалов, имеющих сложное внутреннее строение. То же самое относится и к объектам естественного биологического происхождения, таким как кость. Свойства таких материалов определяются, в том числе, их строением, структурой [5, 6, 11, 16, 17, 22, 32, 41, 52].
В частности, известно, что трабекулярная (губчатая) костная ткань является неоднородной пористой анизотропной структурой. Механические свойства губчатой костной ткани также анизотропны и в значительной мере определяются её внутренней архитектурой [17, 22, 26, 27, 29]. Многие задачи биомеханики зубочелюстной и опорнодвигательной систем человека требуют описания напряжённо-деформированного состояния губчатой костной ткани с учётом формирования её структуры во времени при изменении внешних нагрузок [5, 11-14, 17, 50, 51, 56].
Например, при постановке задач об определении напряжённо-деформированного состояния в нижней челюсти человека необходимо учитывать не только неоднородность свойств твёрдых и мягких тканей, но и их внутреннюю © Киченко А. А., Тверье В.М., Няшин Ю.И., Осипенко М. А., Лохов В.А., 2012
Киченко Александр Александрович, старший преподаватель кафедры теоретической механики, Пермь Тверье Виктор Моисеевич, к.т.н., доцент кафедры теоретической механики, Пермь Няшин Юрий Иванович, д.т.н., профессор, завкафедрой теоретической механики, Пермь Осипенко Михаил Анатольевич, к.ф.-м.н., доцент кафедры теоретической механики, Пермь Лохов Валерий Александрович, к.ф.-м.н., доцент кафедры теоретической механики, Пермь
структуру [11-14, 17, 22, 32]. Таким образом, необходимо иметь способ количественного описания формирующейся под воздействием изменяющегося биомеханического давления структуры костной ткани для различных отделов зубочелюстной системы [5, 6, 11-14, 17].
В связи с этим возникает необходимость введения величины, которая учитывала бы структурные особенности губчатой кости и могла бы быть легко встроена в зависимость строение - свойства материала, иначе говоря, в количественное описание микроструктуры кости [16, 52]. Для подобной величины необходимо сформулировать соотношение, способное описывать упругие свойства материала с учётом его строения.
Также данную величину следует включить в эволюционные уравнения, способные описывать адаптационные процессы и связанные с ними изменения строения костной ткани, происходящие при изменении условий нагружения в губчатой кости. Направленный характер перестройки трабекулярной архитектуры говорит о том, что векторная или тензорная характеристика больше подходит для объективизации процесса [52].
В настоящее время считается, что одним из наиболее удачных способов [5, 9, 17, 22, 32, 41, 52, 57, 60-62] описания локальной структуры многих пористых и композиционных материалов и, в частности, локальной структуры губчатой кости (в том числе степени её анизотропии) является симметричный, положительно определенный тензор второго ранга, названный тензором структуры {fabric tensor) и обозначенный как Н [22, 23, 26, 38, 39, 45, 46, 52].
Тензор структуры, построенный для губчатой костной ткани в соответствии с ранее описанной методикой [6, 17, 37, 57, 60-62], позволяет компактно в тензорной форме описать анизотропию костной структуры, причём его главные значения позволяют охарактеризовать распределение материала вдоль главных направлений [5, 17, 38, 52].
В настоящее время не существует единой формы записи соотношений, связывающих напряжённо-деформированное состояние материала (например, губчатой кости) с его строением (трабекулярной микроструктурой), хотя данной проблеме посвящён целый ряд работ [16, 18-27, 30-36, 38, 40-42, 46, 48, 49, 52, 54, 55, 58, 59, 64]. Однако имеется серия общепризнанных соотношений [21-27, 40, 52, 54, 55], включающих в себя тензор структуры и способных отразить внутреннее строение губчатой кости. То же касается и эволюционных уравнений, описывающих изменения в строении губчатой кости под воздействием различных нагрузок, например, под воздействием изменяющегося биомеханического давления [21, 22, 28, 31, 32, 41, 43, 49]. В данной работе соотношения [21] будут подробно проанализированы и уточнены. На их основе будет представлена постановка начально-краевой задачи о перестройке трабекулярной костной ткани, решение которой позволяет проследить изменение напряженно-деформированного состояния при формировании трабекулярной структуры согласно закону Ю. Вольфа [63].
Вывод определяющего соотношения
Известно, что плотность, пористость и ориентация трабекул в губчатой кости, так же как и упругие свойства костной ткани, изменяются в довольно широком диапазоне в зависимости от исследуемой области. Таким образом, определяющие соотношения должны описывать напряжённо-деформированное состояние губчатой костной ткани, если известны величины, определяющие ориентацию трабекул в рассматриваемой области, и некоторые упругие характеристики материала, которые зависят от плотности костного матрикса, но не зависят от ориентации трабекул (такие характеристики могут быть определены эмпирически после исследования ряда образцов костной ткани).
Экспериментально показано [21, 22, 41], что наблюдаемые в трабекулярной костной ткани деформации не превосходят 0,5%, т.е. они являются малыми [3, 7, 8]. В случае малых деформаций можно воспользоваться основными положениями линейной теории упругости.
Установлено [27], что механические свойства губчатой костной ткани в значительной мере определяются её внутренней архитектурой (геометрией). Костная ткань является анизотропной и неоднородной как в своём строении, так и в своих механических свойствах. При этом в первом приближении считается, что костный матрикс в губчатой кости (пористом упругом теле) изотропен [22, 27] и вся неоднородность губчатой кости связана с геометрией анизотропной микроструктуры трабекулярной костной ткани. В этом случае анизотропия губчатой кости описывается посредством тензора структуры, а тензор упругости и тензор структуры связаны некоторой функциональной зависимостью.
Механические свойства трабекулярной кости (такие как удельная упругость [22, 27]) также зависят от её пористости или связанной с пористостью величины - доли твёрдого объёма кости v, которая определяется как отношение объёма, занимаемого трабекулярной костной тканью в исследуемом образце Vbone, к объёму всего образца губчатой кости Vtotai [41], т.е.
V=VboneL (1)
V ,
total
Таким образом, тензор упругости должен связывать компоненты тензора напряжений с компонентами тензора деформации в упругом материале с наведённой анизотропией [22]. При этом упругие свойства материала зависят от пористости
губчатой кости (т.е. от v) и от ориентации трабекул (т.е. от тензора Н ). Математически [10, 27] это выразится в том, что в самом общем виде тензор напряжений а является
изотропной функцией тензора малых деформаций , тензора структуры Н и доли твёрдого объёма кости v:
5 = а(в, Н, v). (2)
Здесь мы имеем дело с тензорной функцией а двух тензорных аргументов и Н и также скалярного инварианта v. При этом должно выполняться условие
Q-5• QT = 5(Q• s-QT,Q Н• QT,v) (3)
для всех ортогональных преобразований Q [16, 44]. Дальнейшие преобразования связаны с поиском наиболее простой формы определяющего соотношения (2), не противоречащей теории определяющих соотношений и не уменьшающей их общности [10, 16].
Исходя из своего определения [17], тензор структуры является симметричным, положительно определённым тензором второго ранга. Тензор малых деформаций также является симметричным тензором второго ранга. Будем считать, что тензор напряжений также симметричный (т.е. предположим, что мы действуем в рамках теории симметричной упругости). В этом случае в соответствии с выводами, представленными в работах [2, 7, 10, 16] для описанного класса изотропных функций, можно конкретизировать соотношение (2) в виде следующего полиномиального разложения:
(5 = I0 Ё + /js +I2 s2 + /3Н + /4Н2 + I5(s -Н + Н -s) +
+I6(s-Н2 +Н2-s) + I7(s2-Н + Н-s2) + I8(s2 Н2 +Н2-s2), (4)
а б
Рис. 1. Тестовые изотропная (а) и структурированная (б) пористые микроструктуры
где Е - единичный тензор; /а - изотропные скалярные функции. Соотношение (4) будет получено в том случае [16], если функции /а (по а не суммировать!) будут представлены в виде полиномов от элементов целого рационального базиса для тензоров-аргументов ~ и Н [10]. То есть полиномиальные функции 1а должны зависеть как от инвариантов тензоров-аргументов, так и от скалярной величины V:
/а = /„(V, 2,3,1гН, 1гН2,1гН3,1г(8 • Н), 1г(8 • Н2), 1г(~2 • Н), ц(е2 • И2)). (5)
В дальнейшем костная ткань будет рассматриваться как линейно-упругий материал с наведённой анизотропией [22]. Известно, что для случая линейной упругости связь между тензором напряжений и тензором деформации выражается законом Гука [3], при этом тензор упругости полностью характеризует линейноупругое механическое поведение материала. В этом случае соотношение (4) может быть записано как
5 = /0Е + /1~ + /3Н + /4Н2 + /5(~Н + Н• ~) + /6(~Н2 +Н2 • ~), (6)
где полиномиальные функции /а можно представить в общем виде как
/а = /а (V, ^ ~, 1х Н, 1г Н2, 1г И3,1г(~ • Н), 1г(8 • И2)). (7)
Таким образом, соотношение (6) в общем виде принимает следующий вид:
~ = С (И, V) • •?, (8)
а упругие свойства материала (т.е. тензор упругости С) зависят от его удельной плотности, ориентации трабекул в пространстве и абсолютной степени анизотропии материала.
При описании структуры костной ткани посредством тензора структуры прежде всего интерес представляет относительная степень анизотропии губчатой кости [17] (за придание жёсткости материалу отвечает величина V [22, 27]). Ранее [17] авторами было показано, что компоненты тензора структуры главным образом отражают преимущественные направления распределения пор в трабекулярной костной ткани, их ориентацию в пространстве. Видно, что для тестовых микроструктур (см. рис. 1) компоненты тензора структуры в определённом диапазоне слабо зависят от размера пор или от их количества (рис. 2, 3 и [17]). Также в ряде работ [47, 55] было замечено, что упругие свойства пористого материала не зависят от размера пор. Итак, тензор структуры в дальнейшем можно нормировать таким образом, что 1г Н = 1.
0 5 10 15 20 R, 25 30 35 мм »0 45
а
Паа 00
мм 2000
[
1
с 2( К) 4< >0 Ы Ю 8( )0 ю 00 12
б
Рис. 2. Зависимость компоненты тензора структуры Паа для изотропного образца: а - от размера пор при их постоянном количестве; б - от числа пор при их постоянном размере
1 ♦ Чаа
\ — Прр 5
ZwU X
0 1 0 1 5 2 0 2 а, 5 3 4M 0 3 5 4 0 4
а
9000 8000 7000 6000 5000 ММ 4000 3000 2000 1000
10'Vr
о - *■
Лаа
*1РР
О 200 400 600 800 1000 1200
А/
б
Рис. 3. Зависимость компонент тензора структуры паа и Прр для анизотропного образца: а - от размера пор при их постоянном числе; б - от числа пор при их постоянном
размере
В этом случае от непосредственно тензора структуры Н удобно [21, 22] перейти к девиатору тензора структуры К, определяемому как
К = Н -1Е. (9)
3
Тензор К по-прежнему отражает распределение пор в губчатой кости, при этом будет выполняться условие 1г К = 0.
В этом случае с учётом введённых обозначений соотношение (6) может быть записано как
б = в0 Е + р18 + Р2К + Р3К2 + р4(в -К + К-в) + р5(в-К2 + К2-в), (10)
где полиномиальные функции ра могут быть представлены в общем виде как
Ра = Ра(у, *г ~^ К2, * К3,1х(в - К), 1х(в - КГ2)) . (11)
Раскроем полученное соотношение. Для этого подставим (11) в (10), при этом следует учитывать, что тензор ~ должен линейно входить в каждое слагаемое соотношения (10). То есть получаем, что
в0 = а11г в + а21х(в - К) + а31г (в - К2),
р1 = К
Р2 = с1 1Х 8 + с2 1х(£ - К) + С3 1Х (в - К2), (1 _)
(12)
вз = ^1 1г 8 + ^2 *г(б - К) + ^3 1г (¿5 - К ),
Р 4 = P0, в 5 = 4.0,
где аа, Ьа, са, ёа, ра и 4а - полиномиальные функции вида
/а= /а(у,1г£2,1ГК3) . (13)
В итоге получаем:
б = ( 1х в + а21х(в- К) + а31х(в- К2) ) + Ь0в + (с11х в + с21х(в- К) + с31х(в- К2)) 2 +
_ _ _ _ _ _ (14)
+ ( (гв + (1х(в-К) + ё31х(в-К2))К2 + р0 (в-К + (в) + ) (в-К2 + К2 - в).
Соотношение (14) аналогично определяющему соотношению из работы [27] (см. также [21, 22, 52]). Отличие от работы [27] заключается в том, что в качестве параметра структуры авторами используется девиатор тензора структуры.
Воспользуемся упрощающим предположением из работ [10, 27], а именно после подстановки в соотношение (13) полиномиальных функций (14) отбросим все слагаемые, чья общая степень превосходит п = 2. Линеаризация (п = 1) определяющего соотношения невозможна, поскольку в этом случае независимые компоненты тензора упругости станут функционально связанными [27]. В результате подобного упрощения было получено следующее соотношение:
б = (а10 +апу)1г в + а20 й(в- К) ) + (Ь00 + Ь01у)в + (с10 й в) К + р00 (в- К + К- в), (15)
где аш, а„, а20, Ь00, Ьш, сш, Р00 - константы.
Видно, что соотношение (15) с точностью до констант совпадает с соотношением, полученным в работе [21]. Осуществим переобозначения, приводящие соотношение (15) к принятому в литературе [22] виду. Для этого введём величину е - изменение доли твёрдого объёма кости относительно отсчётной величины v0, т.е.
e (16)
и заменим ею текущую долю твёрдого объёма кости v [21]. Также заменим константы a10, a11, a20, b00, b01, c10, p00 соответствующими работе [21] постоянными величинами (при этом необходимо полагать, что a20 = с10 [22]). В итоге получим:
6 = (g + ^2e)(tr ё)E + (g3 + g4e)g + g5 (ё •K + ■K • ё) + g6 (tr(K • ё) E + (tr ё) K), (17)
что эквивалентно соотношению (14) из работы [21]. Здесь g1-g6 - константы [21], имеющие размерность [ГПа], поскольку входящий в соотношение (17) тензор структуры был нормирован. Эти константы были определены в работе [55] после серии экспериментов на различных образцах губчатых костей человека и крупного рогатого скота. Тензор упругости, соответствующий соотношению (17), принимает следующий вид:
Cijkl = (g1 + g2e)5ij§kl + (g3 + g4e)5ik5 jl + g5 (5ikKlj + Kik5lj) + g6 (5ijKkl + Kij5kl )- (18)
Эволюционные уравнения закона Вольфа для губчатой костной ткани
Костная ткань живого человека является динамической сложной структурой, в которой происходят постоянный обмен веществ, анаболические и катаболические процессы, разрушение старых и создание новых костных трабекул. Заметим, что термин «анаболический» означает «относится к или способствует анаболизму» (процессу ассимиляции питательных веществ и превращению их в составные части клеток, сопровождаемому при этом энергетическими затратами). Термин «катаболический» означает «относится к катаболизму» (процессу расщепления в организме молекул на более простые, часто сопровождаемому высвобождением энергии).
Механические свойства костей подчиняются тем же принципам, что и несущие конструкции, сделанные людьми. Однако способность позвоночных приспосабливать строение своей костной ткани к приложенной нагрузке приводит к такой структуре, которая является очень сложной и, будучи здоровой, исключительно эффективной [17, 41].
Среди законов, описывающих поведение костной ткани под влиянием каких-либо факторов, например нагрузок, закон Вольфа (Wolff’s law) является наиболее известным [63], хотя его точная математическая запись до настоящего времени четко и однозначно не разработана. Данный вопрос до сих пор является краеугольным камнем современных исследований [5, 17]. Закон Вольфа отмечает изменение кости (или мягких тканей) вследствие функциональных требований. Каждое изменение в форме или функции сопровождается определенными изменениями во внутренней архитектуре и во внешней форме. Закон Вольфа применительно к живой костной ткани звучит следующим образом: кость приспосабливает свою внешнюю форму и внутреннюю структуру к тем механическим силам, которые она должна выдержать [5, 17, 41, 63].
При математическом описании адаптационных процессов, происходящих в кости с учётом её внутренней структуры (т.е. при построении кинетических уравнений, включающих в себя слагаемые, отображающие внутреннее строение кости), необходимо дать точную математическую форму записи для закона Вольфа. Для этого воспользуемся подходом, описанным в работах [21, 22, 28]. Закон Вольфа для костной ткани говорит о том, что трабекулярная архитектура губчатой кости 42 ISSN 1812-5123. Российский журнал биомеханики. 2012. Т. 16, № 4 (58): 36-52
в локальной области структурно приспосабливается к местному напряженному состоянию костной ткани. При этом структурная адаптация в живой губчатой кости носит направленный характер и трабекулы располагаются закономерно, сообразно тому, какие внешние нагрузки испытывает данная кость [21, 28]. В частности, установлено, что ориентация трабекул в рассматриваемой области губчатой кости совпадает с главными направлениями тензора напряжений в этой же области [22, 28].
Известно, что при изменении внешней нагрузки, оказывающей воздействие на кость, в соответствии с законом Вольфа происходит перестройка костной ткани. Отсутствие перестройки в кости (равновесное состояние, или гомеостаз) предполагает наличие определённого набора условий, при которых нет никаких изменений в геометрии трабекулярной микроструктуры, резорбции или роста кости. Г омеостатическое состояние обладает определённой специфической структурой
трабекулярной костной ткани, описываемой как (V0, К0) (или (V0, Н0)),
и соответствующим данной структуре специфическим напряженно-деформированным состоянием кости (<5°, ~0). Отметим, что напряжение <5° и деформация ~0 при гомеостазе фактически находятся в диапазоне изменения напряжений и деформации,
в пределах которого не происходит никакой перестройки кости. Величины б0 и ~0 являются здесь усреднёнными значениями напряжений и деформации среды за продолжительный промежуток времени, хотя можно использовать и другие способы
определения б0 и ~0.
Поскольку закон Вольфа отражает тот факт, что главные направления тензора напряжений в состоянии гомеостаза (при отсутствии перестройки) совпадают с ориентацией трабекул в губчатой костной ткани, то в равновесном состоянии главные
оси тензора напряжений б0 должны совпадать с главными осями девиатора тензора
структуры К0 (или тензора структуры Н0), см. [28]. Иначе говоря, тензоры ~0 и К0
должны быть соосными. Совпадение главных осей тензоров б0 и К0 имеет место в
случае, когда скалярное произведение тензора напряжений ~0 и девиатора тензора
структуры К0 коммутативно [4, 15, 28], т.е. если выполняется условие, что
50-К° = К0-50, (19)
то главные оси означенных тензоров совпадают [28].
В большинстве практических задач биомеханики опорно-двигательной системы необходимо описывать не гомеостатическое равновесное состояние, а перестройку костной структуры, происходящую с течением времени, её непрерывную направленную адаптацию к постоянно изменяющимся нагрузкам, например, биомеханическому давлению. При этом изменение нагрузки должно быть достаточным, чтобы вывести губчатую микроструктуру из состояния гомеостаза и запустить в ней адаптационную перестройку.
В частности, перестройка трабекулярной архитектуры губчатой кости может происходить при однократном достаточном изменении нагрузки. Данная ситуация может возникнуть, например, при вживлении в губчатую кость имплантата. Пусть
начальная равновесная губчатая микроструктура описывалась при помощи тензора К0 и соответствующего ему начального напряжённо-деформированного состояния
(~0, ~0). Далее произошло однократное изменение условий нагружения, запустившее процессы перестройки в губчатой кости и стремящееся привести трабекулярную архитектуру к новому гомеостатическому состоянию с характерной структурой К1 и соответствующим напряжённо-деформированным состоянием (б1, в1).
Гомеостаз,
¿<0
Новая нагрузка, і = 0
Перестройка,
* >0
Гомеостаз, t = U
Рис. 4. Схема перестройки трабекулярной костной ткани: а - начальное состояние гомеостаза; б - начало перестройки костной ткани; в - процесс перестройки костной ткани;
г - новое состояние равновесия
Описанный адаптационный процесс схематично показан на рис. 4, где тензоры напряжений 5, деформации ~ и структуры К представлены в виде эллипсоидов. Рис. 4, а отражает положение, существовавшее при I < 0 и соответствующее начальному состоянию гомеостаза. Видно, что три эллипсоида, представляющих
начальные напряжения ~0, деформацию ~0 и структуру К0, соосны, т.е. главные
направления тензоров а0, ~0 и К0 совпадают (см. формулу (19)).
Рис. 4, б отражает ситуацию начала перестройки костной ткани в момент времени I = 0, когда произошло изменение условий нагружения, соответственно, имеет
место новое напряжённое состояние костной ткани а1. Напряжённому состоянию а1
~0+ ~1 ~0+ соответствует новое деформированное состояние в , при этом эллипсоиды а и 8
пока что не соосны. В то же время эллипсоид структуры сохранил своё первоначальное
направление К0 , поскольку адаптационные процессы не могут протекать мгновенно.
Рис. 4, в отражает процесс перестройки костной ткани при I > 0. Видно, что
эллипсоиды а1, ~ (^) и К(^) не соосны, при этом главные направления тензоров ~ (I) и
К (^) изменяются таким образом, чтобы стать соосными главным направлениям нового
напряжённого состояния а1. Условием остановки перестройки может служить (19).
Рис. 4, г отражает новое состояние равновесия, которое было достигнуто по прошествии достаточно большого промежутка времени. Новое гомеостатическое
состояние может быть описано тензорами напряжений а1, деформации ~1 и девиатора
тензора структуры К1, причём представляющие их эллипсоиды вновь соосны. Таким образом, произошла адаптация трабекулярной костной ткани к новым условиям нагружения.
Из вышесказанного видно, что необходима математическая форма записи закона Вольфа, способная отразить не только отсутствие перестройки костной ткани, но и саму перестройку, т.е. изменение структуры, а следовательно, и изменение тензора К (или Н), а также изменение доли твёрдого объёма кости е с течением времени. Отметим, что в ряде работ [21, 22, 28, 31, 32, 49] были предложены уравнения, описывающие адаптацию костной ткани к изменяющимся нагрузкам и использующие тензор структуры. Также ряд авторов [22, 32, 43, 64] отмечал, что закон Вольфа может быть частью более общего закона, определяющего скорость изменения тензора структуры как функционал тензора напряжений, деформации, структуры, различных биомеханических факторов, возраста и времени.
При построении эволюционных соотношений необходимо учитывать изменение величин К и е, неизбежно происходящее при перестройке кости. Для построения кинетических уравнений, позволяющих описать перестройку костной ткани в соответствии с законом Вольфа, воспользуемся подходом, предложенным в работе [21]. Для этого введём величины, характеризующие скорости изменения К и е, т.е.
и е = [21, 22]. При этом скорости изменения К и е зависят от деформации
ёг ёг
~ , параметров структуры К и е.
По аналогии с ранее записанной тензорной функцией (2) скорости изменения К и е могут быть представлены в общем виде как функции двух тензорных аргументов ( ~ и К ) и скалярного инварианта е:
К = ¿(е, К е), (20)
е = /2(а, K, е) • (21)
Здесь К = 0 , т.е. отсутствует изменение ориентации трабекул, при ~ =~0. Изменение доли твёрдого объёма кости отсутствует (т.е. е = 0 ) при е = -V 0 (костный матричный материал полностью резорбирует), е = 1 - V 0 (костный матричный материал заполняет
весь объём исследуемого образца) или при ~ =~0 [21].
Предположим, что структура и изменение доли твёрдого объёма кости непосредственно не влияют друг на друга [21]. Тогда скорость изменения структуры зависит только от самой структуры и девиатора тензора деформации, а скорость изменения доли твёрдого объёма - только от объёмной плотности и объёмной деформации, тогда соотношения (20), (21) могут быть упрощены следующим образом:
К = ж5, К), (22)
е = /2 , е^ (23)
где 8 - девиатор тензора ~ , ~ = ~ - 3(*га) Е. При этом = 0 при ~ =~0, а е = 0 при
е =—^0 , е = 1 -V0 или при 8=80. Условие 8=80 означает, что у текущего ~ и начального ~0 тензоров деформации совпадают их главные значения, при этом тензор ~0 должен быть соосным тензору ~ .
Рассмотрим кинетические уравнения (20), (21). С учётом соотношения (10) уравнение (20) может быть записано как
К = У0Е + У18 +у2К + У3К2 + у4(К• ~ + ~• К) + У5(К~2 • ~ + ~• К2) (24)
и
tr K = 0, (25)
где K = 0 при ~ =~0.
Здесь полиномиальные функции уа могут быть определены аналогично (11) как
Ya = Ya (e, tr s, tr K2, tr K~3, tr(~ • K), tr(~ • K2)) . (26)
Эволюционное соотношение (21) для e можно записать в виде аналогичной полиномиальной функции, т.е.
e = /2(e,tr~, trK2, trK~3,tr(~- K), tr(~ • K2)) , (27)
где e = 0 при e = -Vo , e = 1 - v0 или при ~ =~0.
Далее подставим (26) в (24), при этом следует учитывать, что тензор ~ должен линейно входить в каждое слагаемое. Например,
Y0 = aa-1 tr s + aa2 tr( s • K) + aa3 tr (s • K ),
Y1 =
Y2 = CC1 tr~ + CC2 tr(s • K) + CC3 tr (~ • K2),
Y3 = dd1 tr ~ + dd2 tr(~ • K) + dd3 tr(~ • K2),
(28)
Y 4 = PP0,
Y 5 = Ш,
где aaa, bba, cca, dda, ppa и qqa - полиномиальные функции вида
Ia = Ia(e,tri~2,trK3) . (29)
В результате получим:
K = (aax trs + aa2 tr(s• K) + aa3 tr(s• K2))E + bb0s +
+ (c trs + cc2 tr(s• K) + cc3 tr(s• K2)) +
+ {dd1 trs + dd2 tr(s• K) + dd3 tr(s• K2))K2 +
+ pp0 ( • s + s •K ) + qq0 (2 • s + s •K2).
(30)
Применительно к соотношению (30) осуществим упрощение, при котором отбрасываются все слагаемые, чья общая степень превосходит члены второго порядка (п = 2). Тогда получим:
К = ((аа10 + аапе) 1х 8 + аа201х(8 •К)) Е + (йй00 + ¿¿01е)8 +
_ _ _ (31)
+ (сс101х8)К + рр00 (( + 8• К).
Далее необходимо учесть условие отсутствия перестройки трабекулярной кости,
т.е. что К = 0 при 8 =80. Тогда соотношение (31) преобразуется в
К = (а10 + аа11е) 1х (8 -80) + аа201х((8 -80) • К)) Е + (ЪЪ00 + ЪЪ01е) (8 -80) +
+ сс10№(8-80) К + рр00 (• (8-80) + (8 -в0) • К).
Учтём условие (25) и получим, что
(аа10 + аа11е) № (8 -80) + аа20№((8 -80) • К) =
= - 3 ((( + ЪЪ01е) М8 -80) + 2РР00^ (К •(8 -80))). Подставим (33) в (32):
К = ^(ЪЪ00 + ъКе) (8 -80) - 3 (ЪЪ00 + ЪЪ01е) М8 -80)Е^ + ссю ^ (8 - 80) К +
+РР00 Г к (8-80)+(8-8^ к - |*(К (8-8 0))
(32)
(33)
(34)
Видно, что первое слагаемое соотношения (34) включает в себя девиатор тензора ~ - ~0. Обозначим его как ~ - ~0 и приведём (34) к следующему виду:
К = (Ъ00 + ЪЪ01е) (8 -ё0)) + сс101х (ё - 80) К +
+Цl.(в-в0)+(в-в°).¿-ftr<¿.(в-в°))4 (35)
Соотношение (35) с точностью до констант совпадает с выражением (15) из работы [21]. Приведём уравнение (35) к принятому в литературе виду. Для этого
ёе/ ёе/ ёе/ ёе/ 3
введём следующие обозначения: ЪЪ00 = И1, ЪЪ01 = Н3, сс10 = ЪА, рр00 =-Ь2 [21].
В итоге получается соотношение
К = (^ + И3е)(:ё-ё0) + Л4№(8-80) К +
+и2 Г((-(8-в0))) - -3 ((8-в 0)+(8-8°). К)), (36)
где И1-И4 - константы [21], имеющие размерность [сут-1]. Данные константы
подбираются эмпирически таким образом, чтобы перестройка костной ткани происходила за время, соответствующее природной действительности (в норме адаптация губчатой кости происходит примерно за 160 дней [21, 22, 32]).
Аналогично из (27) с учётом обозначений для соотношения (16) из работы [21] можно определить, что
е = (/1 + /2е)0-г 5 - 0) + ./з(tr(i~•(а- 50))), (37)
где /1-/3 - константы [21]. Данные коэффициенты также имеют размерность [сут-1] и подбираются таким образом, чтобы перестройка костной ткани происходила за время, соответствующее природной действительности (т.е. аналогично постоянным Л1-Л4).
С учётом предположений (22) и (32) уравнения (36) и (37) могут быть упрощены следующим образом:
K = -~0) + h2^(tr(K~-(e- е0)))-)(~-(е-~0) + (е- е0)-K~)j, (38)
e = (fi + he)(tr ~ - tr~0) • (39)
Общая постановка начально-краевой задачи о перестройке
ТРАБЕКУЛЯРНОЙ КОСТНОЙ ТКАНИ
Для решения вопросов об исследовании напряжённо-деформированного состояния губчатой костной ткани и протекающих в ней адаптационных процессов будем использовать полученные соотношения (17), (36) и (37).
Рассмотрим ограниченную область V трёхмерного евклидова пространства Е3. Границу области обозначим как £ (£ = и 8и), тогда замыканию области
соответствует V ( V = V и £ ).
В этом случае полная постановка задачи о перестройке трабекулярной костной ткани должна включать в себя: уравнение равновесия
V • а = 0, х е V, ^ > 0, определяющее соотношение, соответствующее уравнению (17):
о = С(К, е) • • ~ , х е V, ^ > 0,
(40)
(41)
эволюционное уравнение, описывающее изменение ориентации трабекул в рассматриваемой области и соответствующее соотношению (36):
dK
dt
= fi(~, K, e), trK = 0, x eV, t > 0,
(42)
эволюционное уравнение, описывающее изменение плотности губчатой костной ткани в рассматриваемой области и соответствующее соотношению (37):
de
dt
= f2(s, K, e), X e V, t > 0,
геометрические соотношения Коши:
1
2
е = ^(Vu + wV), x eV, t > 0,
граничные условия:
П - а = P(t), X e Sa, t > 0,
u = U(t), X e Sw, t > 0,
начальные условия:
K = K0, e = e0, X eV, t = 0,
о
P = P0, X e S^, U = U0, X e S
0
t=0
(43)
(44)
(45)
(46)
(47)
(48)
Видно, что полная постановка задачи в декартовой системе координат в самом общем случае включает в себя 21 скалярное уравнение: три уравнения равновесия, шесть уравнений из определяющих соотношений, шесть кинетических уравнений и
шесть геометрических соотношений. В систему уравнений входит 21 неизвестная: шесть компонент тензора напряжения с j , шесть компонент тензора деформации s j,
пять компонент девиатора тензора структуры Kj, изменение доли твёрдого объёма
кости e и три компоненты вектора перемещений Uj. Таким образом, полученная
система будет статически определимой и можно определить все неизвестные величины.
Далее авторами будет рассмотрено применение соотношений (41)-(43) (или, что то же самое, соотношений (17), (36) и (37)) с учётом начальных условий (47), (48), но без учёта уравнений (40), (44)—(46). Это можно сделать, если предположить, что нагрузки, действующие на рассматриваемую область (а значит и напряжённое состояние рассматриваемой области), не изменяются непрерывно и не вызывают соответствующих перманентных адаптационных процессов в трабекулярной микроструктуре.
Такое возможно, например, при однократном изменении нагрузки и,
соответственно, изменении старого напряжённого состояния а0 на новое с1 . Тогда система уравнений (41)-(43) в прямоугольной декартовой системе координат в общем случае может быть преобразована в систему, включающую в себя двенадцать скалярных уравнений с двенадцатью скалярными неизвестными.
Заключение
На основе существующих подходов [17, 21, 22, 28] был показан вывод соотношений, способных описать перестройку трабекулярной костной ткани, происходящую вследствие изменяющихся условий нагружения. Осуществлена постановка начально-краевой задачи о перестройке трабекулярной костной ткани, включающая в себя эти соотношения. Далее авторами будет рассмотрена локальная область трабекулярной костной ткани, находящаяся в состоянии гомеостаза в течение достаточно долгого промежутка времени, причём соответствующие напряжённо-деформированное состояние структуры и её архитектура считаются известными. В начальный момент времени произойдёт однократное изменение условий нагружения, приводящее к перестройке трабекулярной микроструктуры. Новое напряжённое состояние в рассматриваемой области при этом также известно и не будет изменяться в течение достаточно долгого промежутка времени.
Благодарности
Работа выполнена при частичной поддержке РФФИ (грант 11—01—00910—а).
Список литературы
1. Билич Г.Л., Сапин М.Р. Анатомия человека. - М.: Издательский дом ОНИКС, 1998. - Кн. 1. - 236 с.
2. Грин А., Адкинс Дж. Большие упругие деформации и нелинейная механика сплошной среды. - М.: Мир, 1965. - 455 с.
3. Демидов С.П. Теория упругости. - М.: Высш. шк., 1979. - 432 с.
4. Димитриенко Ю.И. Тензорное исчисление. - М.: Высш. шк., 2001. - 575 с.
5. Киченко А.А., Тверье В.М., Няшин Ю.И., Симановская Е.Ю., Еловикова А.Н. Становление и развитие классической теории описания структуры костной ткани // Российский журнал биомеханики. - 2008. - Т. 12, № 1. - С. 68-88.
6. Киченко А.А., Тверье В.М., Няшин Ю.И., Заборских А.А. Экспериментальное определение тензора структуры трабекулярной костной ткани // Российский журнал биомеханики. - 2011. - Т. 15, № 4. -С. 78-93.
7. Лурье А.И. Нелинейная теория упругости. - М.: Наука, 1980. - 512 с.
8. Поздеев А.А., Трусов П.В., Няшин Ю.И. Большие упруго-пластические деформации. - М.: Наука, 1986. - 232 с.
ISSN 1812-5123. Российский журнал биомеханики. 2012. Т. 16, № 4 (58): 36-52 49
9. Салтыков С.А. Стереологическая металлография. - М.: Металлургия, 1958. - 122 с.
10. Спенсер Э. Теория инвариантов. - М.: Мир, 1974. - 160 с.
11. Тверье В.М., Симановская Е.Ю., Еловикова А.Н., Няшин Ю.И., Киченко А.А. Биомеханическое
описание структуры костных тканей зубочелюстной системы человека // Российский журнал
биомеханики. - 2007. - Т. 11, № 1. - С. 9-24.
12. Тверье В.М., Симановская Е.Ю., Няшин Ю.И. Атрофический синдром, связанный с изменениями биомеханического давления в зубочелюстной системе человека // Российский журнал биомеханики. - 2006. - Т. 10, № 1. - С. 9-13.
13. Тверье В.М., Симановская Е.Ю., Няшин Ю.И. Методы биомеханического моделирования зубочелюстной системы человека // Современные проблемы биомеханики. Биомеханика: достижения и перспективы: научный совет РАН по биомеханике. - М.: Изд-во МГУ, 2006. - Вып. 11. -С. 226-236.
14. Тверье В.М., Симановская Е.Ю., Няшин Ю.И., Киченко А.А. Биомеханический анализ развития и функционирования зубочелюстной системы человека // Российский журнал биомеханики. - 2007. -Т. 11, № 4. - С. 84-104.
15. Трусов П.В., Келлер И.Э. Тензорное исчисление. - Пермь: Изд-во ПГТУ, 2004. - 131 с.
16. Трусов П.В., Келлер И.Э. Теория определяющих соотношений. Ч. 1. Общая теория. - Пермь: Изд-во ПГТУ, 1997. - 98 с.
17. Экспериментальные методы в биомеханике / под ред. Ю.И. Няшина, Р.М. Подгайца. - Пермь: Изд-во ПГТУ, 2008. - 400 с.
18. Beaupre G.S., Orr T.E., Carter D.R. An approach for time-dependent bone modeling and remodeling -theoretical development // J. Orthop. Res. - 1990. - Vol. 8. - P. 651-661.
19. Carter D.R. Mechanical loading trabecular history and skeletal biology // J. Biomech. - 1987. - Vol. 20. -
P. 1095-1109.
20. Carter D.R., Fyhrie D.P., Whalen R.T. Trabecular bone density and loading history regulation of connective tissue biology by mechanical energy // J. Biomech. - 1987. - Vol. 20, No. 8. - P. 785-794.
21. Cowin S.C. An evolutionary Wolff’s law for trabecular architecture // J. Biomech. Engineering. - 1992. -Vol. 114. - P. 129-136.
22. Cowin S.C. Bone mechanics handbook. - Second ed. - New York: CRC Press, 2001. - 1136 p.
23. Cowin S.C. Fabric dependence of an anisotropic strength criterion // J. Mech. Materials. - 1986. - Vol. 5. -
P. 251-260.
24. Cowin S.C. Imposing thermodynamic restrictions on the elastic constant-fabric tensor relationship // J. Biomechanics. - 1998. - Vol. 31. - P. 759-762.
25. Cowin S.C. Mechanical modeling of the stress adaptation process bone // J. Calcified Tissue Int. - 1984. -Vol. 36. - P. S99-S104.
26. Cowin S.C. The mechanical and stress adaptive properties of bone // J. Annals of Biomechanical Engineering. - 1983. - Vol. 11. - P. 263-295.
27. Cowin S.C. The relationship between the elasticity tensor and the fabric tensor // J. Mech. Materials. -
1985. - Vol. 4. - P. 137-147.
28. Cowin S.C. Wolff’s law of trabecular architecture at remodeling equilibrium // J. Biomech. Engineering. -
1986. - Vol. 108. - P. 83-88.
29. Cowin S.C., Mehrabadi M.M. Identification of the elastic symmetry of bone and other materials // J. Biomechanics. - 1989. - Vol. 22. - P. 503-515.
30. Cowin S.C., Mehrabadi M.M. On the identification of material symmetry for anisotropic elastic materials // J. Mech. Appl. Math. - 1987. - Vol. 40. - P. 451-476.
31. Fritton S.P. Computational simulation of trabecular bone adaptation: Ph.D. dissertation. - New Orlean: Tulane University, Departament of Biomechanical Engeenering, 1994.
32. Fung Y.C. Biomechanics. - New York: Springer-Verlag, 1990. - 464 p.
33. Fyhrie D.P., Carter D.R. A unifying principle relating stress to trabecular bone morphology // J. Orthop.
Res. - 1986. - Vol. 4, No. 3. - P. 304-317.
34. Fyhrie D.P., Hollister S.J. Comparision of a trabecular tissue strain remodeling theory to experimental
results // J. Trans. Orthop. Res. Soc. - 1990. - Vol. 15. - P. 107.
35. Goulet R.W., Goldstein S.A., Ciarelli M.J., Kuhn J.K., Brown M.B., Feldkam L.A. The relationship between the structural and orthogonal compressive properties of trabecular bone // J. Biomechanics. -
1994. - Vol. 27. - P. 375-389.
36. Harrigan T.P., Hamilton J.J. Bone remodeling and structural optimization // J. Biomechanics. - 1994. -Vol. 27, No. 3. - P. 323-328.
37. Harrigan T.P., Mann R.W. Characterization of microstructural anisotropy in orthotropic materials using a second rank tensor // J. Mater. Sci. - 1984. - Vol. 19. - P. 761-767.
38. Kanatani K. Distribution of directional data and fabric tensor // Int. J. Eng. Sci. - 1984. - Vol. 22. -P. 149-164.
39. Kanatani K. Stereological determination of structural anisotropy // Int. J. Eng. Sci. - 1984. - Vol. 22. -P. 531.
40. Luo G.M., Cowin S.C., Sadegh A.M., Arramon Y.P. Implementation of strain rate as a bone remodeling stimulus // J. Biomech. Eng. - 1995. - Vol. 117, No. 3. - P. 329-338.
41. Martin R.B., Burr D.B., Sharkey N.A. Skeletal tissue mechanics. Second edition. - New York: SpringerVerlag, 1998. - 392 p.
42. Mullender M.G., Huiskes R., Weinans H. A physiological approach to the simulation of bone remodeling as a self-organizational control process // J. Biomech. - 1994. - Vol. 27, No. 11. - P. 1389-1394.
43. Mullender M.G., Huiskes R. Proposal for the regulatory mechanism of Wolff’s law // J. Orthop. Res. -
1995. - Vol. 13, No. 4. - P. 503-512.
44. Noll W.A. A mathematical theory of mechanical behavior of continuous media // J. Arch. Rational Mech. Anal. - 1958. - Vol. 2. - P. 197-226.
45. Oda M. Fabrics and their effects on the deformation behaviors of sand. - Saitama University: Dept. of Foundation Eng., 1976. - 92 p.
46. Odgaard A., Kabel J., van Rietbergen B., Dalstra M., Huiskes R. Fabric and elastic principal directions of cancellous bone are closely related // J. Biomechanics. - 1997. - Vol. 30. - P. 487-495.
47. Patel M.R. The deformation and fracture of rigid cellucal plactics under multiaxial stress: Ph.D. thesis. -Berkeley: University of California, 1969.
48. Sadegh A.M., Cowin S.C., Luo G.M. Inversion related to the stress-strain-fabric relationship // J. Mech. Mater. - 1991. - Vol. 11. - P. 323-336.
49. Sadegh A.M., Luo G.M., Cowin S.C. Bone ingrowth: an application of the boundary element method to bone remodeling at the implant interface // J. Biomech. - 1993. - Vol. 8, No. 20. - P. 785-794.
50. Simanovskaya E.Y., Bolotova M.Ph., Nyashin Y.I. Mechanical pressure as generator of growth, development and formation of the dentofacial system // Russian Journal of Biomechanics. - 2001. - Vol. 5, No. 3. - P. 14-17.
51. Simanovskaya E.Y., Bolotova M.Ph., Nyashin Y.I., Nyashin M.Y. Masticatory adaptation of the human dentofacial system // Russian Journal of Biomechanics. - 2002. - Vol. 6, No. 4. - P. 15-61.
52. Telega J.J., Jemiolo S. Fabric tensor in bone mechanics // J. Engineering Transactions. - 1998. - Vol. 46. -P. 3-26.
53. Truesdell C., Noll W. The nonlinear field theories of mechanics / ed. S. Flugge, C. Truesdell. - Berlin: Springer, 1965. - Bd. III/3. - 602 p.
54. Turner C.H., Cowen S.C. On the dependence of elastic constants of an anisotropic porous material upon porosity and fabric // J. Mater. Sci. - 1987. - Vol. 22. - P. 3178-3184.
55. Turner C.H., Cowen S.C., Rho J.Y., Ashman R.B., Rice J.C. The fabric dependence of the orthotropic elastic constants of cancellous bone // J. Biomechanics. - 1990. - Vol. 23. - P. 549-561.
56. Tverier V.M., Nyashin Y.I., Simanovskaya E.Y. Biomechanical description of the functional features of human masticatory apparatus in norm and in various pathological processes // J. Mechanika w Medycynie. -2008. - No. 9. - P. 147-160.
57. Underwood E. Quantitative stereology. - Mass.: Addision Wesley, 1970. - 370 p.
58. Weinans H., Huiskes R., Grootenboer H. Convergence and uniqueness of adaptive bone remodeling //
J. Trans. Orthop. Res. Soc. - 1989. - Vol. 14. - P. 310.
59. Weinans H., Huiskes R., Grootenboer H. Numerical comparisons of strain-adaptative bone remodeling theories // Trans. First World Congress of Biomechanics. - San Diego, 1990. - Vol. II. - P. 75.
60. Whitehouse W.J. Irregularities and asymmetries in trabecular bone in the innominate and elsewhere // J. Metab. Bone Dis. Rel. Res. - 1981. - Vol. 2. - P. 271-278.
61. Whitehouse W.J. The quantitative morphology of anisotropic trabecular bone // J. Microscopy. - 1974. -Vol. 101. - P. 153-168.
62. Whitehouse W.J., Dyson E.D. Scanning electron microscope studies of trabecular bone in the proximal end of the human femur // J. Anat. - 1974. - Vol. 118. - P. 417-444.
63. Wolff J. Das gesetz der transformation der knochen. - Berlin: Hirshwald, 1892. - 416 s.
64. Zysset P.K., Curnier A. An alternative model for anisotropic elasticity based on fabric tensor // J. Mech.
Mat. - 1995. - Vol. 21. - P. 243-250.
STATEMENT OF INITIAL BOUNDARY VALUE PROBLEM ON THE TRABECULAR BONE TISSUER REMODELLING
A.A. Kichenko, V.M. Tverier, Y.I. Nyashin, М^. Osipenko, V.A. Lokhov (Perm, Russia)
Bone tissue is a heterogeneous, anisotropic material; structural features of trabecular bone tissue can be described by means of the fabric tensor. Biomechanical modelling tasks demand to study the history of the formation of bone structures in time under both physiological and pathological loadings. This is possible to implement if there is both a constitutive relation which allows us to connect the stress tensor, the strain tensor, and the fabric tensor, and evolution equations which allows us to describe the evolution of the fabric tensor and bone density. These equations and relations of S.C. Cowin for trabecular bone are analyzed in detail; some inaccuracies and omissions are eliminated in these relations. An initial boundary value problem on the trabecular bone tissue remodelling is stated. The solution of this problem allows us to trace changes in the stress-strain state at the trabecular structure formation according to Wolff’s law.
Key words: biomechanical modelling, initial boundary value problem, constitutive relation, evolution equation, bone tissue structure, trabecular (cancellous) bone tissue, fabric tensor, Wolff’s law, state of homeostasis (physiological equilibrium), bone tissue remodelling.
Получено 01 ноября 2012
