О приложении теории перестройки трабекулярной костной ткани Текст научной статьи по специальности «Математика»

УДК 531/534:[57+61]

О ПРИЛОЖЕНИИ ТЕОРИИ ПЕРЕСТРОЙКИ ТРАБЕКУЛЯРНОЙ КОСТНОЙ ТКАНИ

А.А. Киченко, В.М. Тверье, Ю.И. Няшин, М.А. Осипенко, В.А. Лохов

Кафедра теоретической механики Пермского национального исследовательского политехнического университета, Россия, 614990, Пермь, Комсомольский проспект, 29, e-mail: nyashin@inbox.ru

Аннотация. Приведены определяющие соотношения, позволяющие описать напряжённо-деформированное состояние губчатой костной ткани с учётом её структуры, и эволюционные соотношения, описывающие адаптационные процессы, происходящие в кости. С учётом ранее поставленной авторами начально-краевой задачи, включающей в себя определяющие и эволюционные соотношения (S.C. Cowin), разработан алгоритм решения и на ряде примеров показана эволюция трабекулярной костной ткани при изменении напряжённо-деформированного состояния. Для определения коэффициентов, входящих в эволюционные уравнения, проведён вычислительный эксперимент, в котором учитывалось, что адаптация костной ткани должна происходить за 160 дней. Результаты показывают различный характер влияния изменения нагрузки на процесс формирования структуры и не противоречат закону Ю. Вольфа.

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

Введение

В предыдущей работе авторами (см. данный выпуск журнала) был показан вывод определяющего соотношения, включающего в себя тензор структуры, т.е. величину, способную отразить внутреннее строение губчатой кости [5, 7, 9, 17]. Данное соотношение связывает напряжённо-деформированное состояние в губчатой костной ткани с её строением (трабекулярной микроструктурой) и имеет следующий вид:

~ = (g1 + g2e)(tr ~)E + (g3 + g4e)~ + g5 (~ • K + K • ~)+ g6 (tr(K • ~) E + (tr ~)K)), (1)

где K - девиатор тензора структуры H, нормированный таким образом, что tr K = 0 [8]; e - изменение доли твёрдого объёма кости относительно отсчётной величины v0; g1-g6 - константы [8, 19], имеющие размерность [ГПа]. Эти константы были определены в работе [19] после серии экспериментов на различных образцах губчатых костей человека и крупного рогатого скота.

Также на основе закона Вольфа [3, 7, 15, 16, 20] был описан предполагаемый механизм перестройки трабекулярной микроструктуры, происходящий вследствие адаптационных процессов, перманентно протекающих в живой костной ткани. © Киченко А.А., Тверье В.М., Няшин Ю.И., Осипенко М.А., Лохов В.А., 2012

Киченко Александр Александрович, старший преподаватель кафедры теоретической механики, Пермь Тверье Виктор Моисеевич, к.т.н., доцент кафедры теоретической механики, Пермь Няшин Юрий Иванович, д.т.н., профессор, завкафедрой теоретической механики, Пермь Осипенко Михаил Анатольевич, к.ф.-м.н., доцент кафедры теоретической механики, Пермь Лохов Валерий Александрович, к.ф.-м.н., доцент кафедры теоретической механики, Пермь

В соответствии с законом Вольфа для губчатой костной ткани [3, 7, 12] были получены эволюционные соотношения, описывающие изменения в строении губчатой кости под воздействием различных нагрузок. Данные кинетические уравнения способны описать изменение девиатора тензора структуры К и доли твёрдого объёма кости е. А именно:

ИКС — —

---= (Н + Н>е)(ё-ё0) + Н41х(ё-ё0) К +

г Л 3 ( (2)

+н ^(1г( к ■ (ё-ё0))) -) (■ (ё-ё0)+(ё-ёо). к )

и

е = (/ + /2е) (^Г ё - tr ё 0) + ./3(*(К ■ (ё-ё °))), (3)

где к1-к4 и /1—/3 — константы [8], имеющие размерность [сут-1], определяемые эмпирически так, чтобы перестройка костной ткани происходила за время, соответствующее природной действительности (в норме адаптация губчатой кости происходит примерно за 160 дней [8, 9, 15, 16]).

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

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

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

Общие предположения

В качестве примера применения соотношений (1)—(3) (или, что то же самое, соотношений (41)—(43) из предыдущей статьи авторов) рассмотрим локальную область трабекулярной костной ткани, находящуюся в состоянии гомеостаза в течение достаточно долгого промежутка времени. В течение этого времени ( < 0) напряжённо-

деформированное состояние структуры описывается как ( ~0, ~ 0), архитектура костной

ткани — (К0, е0 ). В момент времени ^ = 0 происходит однократное изменение условий нагружения, приводящее к перестройке трабекулярной микроструктуры. Напряжённое

состояние в рассматриваемой области при этом описывается как а1 и не изменяется в течение достаточно долгого промежутка времени (^ > 0).

Будем считать, что начальные напряжения известны, главные значения тензора

заданы в следующем виде:

О 0 0

N о а 0 а2 0

0 V 0 а3

(4)

Величины, описывающие начальную равновесную архитектуру (К0, е0), могут быть определены экспериментально [3-7, 9, 16]. При этом в равновесном состоянии главные оси тензора напряжений <5° должны совпадать с главными осями девиатора тензора структуры К0 (т.е. скалярное произведение тензоров 5° и К0 должно быть

коммутативно). Будем считать, что главные значения тензора К0 заданы в следующем виде:

' K? 0

*° = 0 K20

0 0

0

0

Л

(5)

а изменение доли твёрдого объёма кости равно е .

Начальное деформированное состояние может быть определено из соотношения

(1), т.е. компоненты тензора ~0 также известны. В главных осях тензор ~0 может быть представлен как

а0 =

^0

S1

0

0

0

0

S3

(б)

причём а°, K0 и s0 должны быть соосны [2, 12].

Необходимо отметить, что исходные характеристики архитектуры губчатой

костной ткани (K0, е0) не всегда можно определить экспериментально. Также

величины K0 и е0 нецелесообразно определять при решении модельных примеров. В работе [8] описан подход, при котором задаётся начальное деформированное состояние исходя из теории lazy zone (некоторый диапазон, в пределах которого отсутствует перестройка трабекулярной микроструктуры). Известно [8], что для определённого диапазона деформаций в губчатой кости адаптационные процессы

отсутствуют. Если задать начальное деформированное состояние ~0 из данного диапазона, то начальная архитектура губчатой кости (K0, e0 ) может быть определена из соотношения (1) с учётом известного начального напряжённо-деформированного состояния (5°, ~0 ). При этом после завершения адаптационных процессов и остановки перестройки трабекулярной микроструктуры её новое деформированное состояние вновь должно находиться в диапазоне, определяемом границами lazy zone. Фактически

имеет место равенство ~0 = в1 [8, 9], т.е. в своих главных осях главные значения тензоров ~0 и в1 совпадают.

Зададим новое напряжённое состояние, причём главные значения тензора а1 представлены в следующем виде:

(7)

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

~1 1 ~1

обозначить как К и е , деформированное состояние в этом случае обозначим как в .

а1 0 0

а = 0 а2 0

0 V 0 а3

0

0

Величины в1, К и е1 будут определяться из решения системы уравнений (1)—(3), при

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

К(г) и е^) .

Введём обозначения: пусть Р(<~0) - прямоугольная декартова система

координат, где её оси совпадают с главными направлениями тензора а0, Р(а1) -прямоугольная декартова система координат, где её оси совпадают с главными направлениями тензора а1. Тогда новое напряжённое состояние а1 в независящих от времени системах координат Р(а0) и ^(а1) будет связано следующим образом:

а1 р(~0) = д .51 р(а1). ят

(8)

где V - некоторое ортогональное преобразование, переводящее из главных осей тензора а1 к главным осям тензора 50 .

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

Рассмотрим случай, при котором действующая нагрузка изменилась по величине, но не по направлению, т.е. главные направления тензоров <50 и <51 соосны, а системы координат Р(а0) и Р(а1) совпадают (т.е. V = Е ).

Тогда текущие компоненты тензора деформации є (ї) и девиатора тензора

структуры К (ї) в системе координат Р(а1) будут иметь следующий вид:

є (ї)

р (~~0) =

0

К (ї)

Р(<0) =

К

0

0

0

К

0

0

0 Єп

(9)

2

/

0

0

0 0 — Кі — К 2

V 1 2 У

(10)

Подстановка выражений (8)—(10) в (1)-(3) приведёт к системе, состоящей из шести уравнений с шестью неизвестными: в1(^), в2(0, в3(^), К1(^), К2(^) и е(^). С учётом

выражений (7)—(10) соотношение (1) может быть представлено в компонентной форме

в виде системы нелинейных алгебраических уравнений:

а1 = (#1 + ^2е)(в1 + в2 +в3) + (ё3 + ё4е )в1 +

+2^5в1К1 + #6(2в1К1 +В 2(К1 + К 2)-взК2Х

а2 = (ё1 + ё2е)(в1 + в2 + в3) + (#3 + ё4е) в 2 +

+2 ё5в2К 2 + ёб (2в2 К2 + в1 (К1 + К2) - в3 К1),

= (ё1 + ё2е) (в1 + в2 +в3) + (ё3 + ё4е)в3 --2ё5в3 (К1 + К2) - ёб (в1 К2 + в2К1 + 2в3 (К1 + К2)).

!88М 1812-5123. Российский журнал биомеханики. 2012. Т. 16, № 4 (58): 53—72

(11)

(12)

(13)

Б

1

В

2

0

Аналогично с учётом выражений (7)—(10) кинетические уравнения (2) и (3) могут быть представлены в компонентной форме в виде системы нелинейных дифференциальных уравнений:

—К И + И^е 0 0 0чч , ^ \

—— =------3---((2в1 -в2-в3)-(2в1 -в2-вз)) + И4К1((в1 +в2 +в3) -

Л 3 (14)

_(в<0 + в2 +в3)) + И2(К2(в2 -в2) -2К1(в1 -в0) -(К1 + К2)(в3 -в0)Х

-К2 И + Ие 0 0 0\\ 7 тг и \

—— =------3---((2в2-в1 -в3)-(2в2-в0 -взХ) + И4К2((в1 +в2 +в3)-

Л 3 (15)

(в0 +в2 +в0)) + И2(К1(в1 -в0) - 2К2(в2 -в2) - (К1 + К2)(в3 -в0)),

— = (11 + /2е)((в1 + в2 + в3) - (в0 + в2 + в3)) +

Л (16)

+ЖК1(в1 -в<0) + К2(в2 -в2)-(К1 + К2)(в3-в'0)).

Полученная система уравнений (11)—(16) не может быть линеаризована. Это связано с тем, что соотношения, управляющие эволюционными процессами,

обязательно должны быть нелинейными, в противном случае будет отсутствовать механизм обратной связи между адаптационными процессами, протекающими в трабекулярной костной ткани, и её напряженно-деформированным состоянием [8, 10, 11]. К тому же будет отсутствовать взаимовлияние между пористостью губчатой кости и распределением в ней трабекул. Полученная система уравнений (11)—(16) может быть решена численно.

Для упрощения полученной системы можно воспользоваться следующим

допущением. Пусть структура (т.е. девиатор тензора структуры К(?)) и пористость кости (т.е. изменение доли твёрдого объёма кости е(^)) непосредственно не влияют друг на друга [8, 9]. Тогда скорость изменения структуры будет зависеть только от самой структуры и девиатора тензора деформации, а скорость изменения доли твёрдого объёма — только от объёмной плотности и объёмной деформации. Соотношения (2), (3) могут быть упрощены следующим образом:

К = //(е, К), (17)

е = /2 в, е). (18)

Рассмотрим численный пример для системы (11)—(16). Для этого в полученной системе зададим материальные константы. Коэффициенты для определяющего соотношения (1) равны следующим величинам: ё1 = 154,9; #2 = 1147; ёз = 612,9; ё4 = 4536; ё5 = 2384; #6 = 510,8 (эти коэффициенты имеют размерность [ГПа]) — см. [8, 19].

Для определения коэффициентов, входящих в кинетические уравнения (2) и (3), был проведён вычислительный эксперимент, где учитывалось, что адаптация костной ткани должна происходить за 160 дней. Эксперимент показал, что коэффициенты для эволюционных соотношений (2) и (3) с учётом предположения о независимости структурной и объёмной перестройки (17) и (18) должны иметь следующий вид: И1 = —5; И2 = 10; И3 = 0, И4 = 0, /1 = —2,5; /2 = 5; /3 = 0 (эти коэффициенты имеют размерность [сут—1]). Отметим, что в работе [8], где также определялись приведённые коэффициенты, была допущена неточность при поиске коэффициентов И1 и /1.

Рассмотрим случай всестороннего сжатия (пример 1) трабекулярной

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

;0 P(S0) =

Г-1,4024 0 0

0

-1,4024 0

0 0

-1,4024

Л

У

Для иллюстративного примера зададим начальное деформированное состояние

~0 исходя из теории lazy zone [8, 9, 12, 16] (поскольку для определения параметров

архитектуры губчатой костной ткани K0 и е0 возникает необходимость рассмотрения реальной костной микроструктуры). А именно: известно, что перестройка

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

:0 P(~0) =

Г- 0,0015 0

0 - 0,0015

0

0

0 0

- 0,0015

Л

У

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

Система уравнений (11)—(13) позволяет численно определить начальные параметры архитектуры губчатой костной ткани К0 и е0, а именно е0 = -0,0179 ,

K

0Р(а0) =

Г0 0 0^ 000 000

У

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

Зададим новое напряжённое состояние, равномерно увеличив всестороннее сжатие:

а Да1) =

Г-1,5472 0 0

0 0 >

-1,5472 0

0 -1,5472у

и осуществим численное решение системы уравнений (11)—(16). Следует заметить, что в данном и последующих примерах тензоры а° и а1 являются соосными, т.е. системы отсчёта Р(а0) и ^(а1) совпадают.

Результаты, полученные по истечении 160 суток, показаны на рис. 1. Видно, что костная микроструктура вновь пришла в гомеостатическое равновесие, описываемое

следующими величинами: е1 = -0,0059 ,

а

60 80 100 120 140 160

/. сут

б

Рис. 1. Изменение структуры (компонент девиатора тензора структуры К1(() (а), К2(/) (б)) и плотности (доли твёрдого объёма кости е(/)) трабекулярной микроструктуры (в)

Г1 де1) _

'- 0,0015 0 0 Л

0 -0 0015 0

V 0 0 - 0,0015,

K P(g1) _ ' 0 0 0 >

0 0 0

V0 0 0 V

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

Отметим, что вследствие применения теории lazy zone полученное конечное

деформированное состояние в1 соответствует гомеостатическому, а главные значения

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

в

а

б

Рис. 2. Изменение структуры (компонент девиатора тензора структуры К1(/) (а), К2(/) (б)) и плотности (доли твёрдого объёма кости в(/)) трабекулярной микроструктуры (в)

Рассмотрим другой пример (пример 2). Предположим, что начальное напряжённо-деформированное состояние и начальные архитектурные параметры такие же, как и в предыдущем примере 1. Новое напряжённое состояние задаётся следующим образом:

а Да1) _

Г-1,4024 0 0

0

-1,4024 0

0 0

-1,5472

у

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

перестаёт быть однородным.

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

описываемое следующими параметрами губчатой микроструктуры: в1 _-0,0139 и

К

1 Да1) _

Г- 0,005 0 0

0

- 0,005 0

0 ^ 0

0,010

Также видно, что губчатая микроструктура вследствие неравномерного нагружения перестала быть изотропной, при этом также изменилась доля твёрдого объёма кости (а значит, и её пористость).

в

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

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

Изменение нагрузки трабекулярной микроструктуры

ПО ВЕЛИЧИНЕ И НАПРАВЛЕНИЮ

Далее в качестве примера рассмотрим ситуацию, частично разобранную в работе [8]. В начальный момент времени известны (или могут быть определены)

начальные тензоры 50, s0 и K0 в виде соотношений (4)—(6), а также начальная доля

твёрдого объёма кости е0 . Отметим, что главные значения тензоров КГ° и s0 заданы

в системе координат P(5°). Новое напряжённое состояние 51 будет определено следующим образом:

а1 P(50) = Q(t) -51 P(51) • Q(xf, (19)

где Q(t) — ортогональное преобразование от главных направлений тензора а1 к

~0 ~1 ~1

главным направлениям тензора 5 . Тензор 5 в системе координат Р(а ) может быть

представлен в виде (7), ортогональное преобразование Q(t) задаётся следующим образом:

(cos т 0 - sin тЛ

Q(T) =

О 1 О

у sin т О cos Ту

(2О)

где т - угол в плоскости 1-3 между главными осями (системами координат) ^(а1)

и Р(д°) - см. рис. 3. Отметим, что в предложенном примере авторы [8] изменяли не только направление приложения нагрузки, но и её величину, что нецелесообразно для простого модельного примера.

Рис. 3. Иллюстрация процесса перестройки структуры трабекулярной костной ткани из одного состояния физиологического равновесия в другое (иллюстрация трабекулярной структуры взята из работы [8]): а - исследуемая область; б - сечение костной ткани в плоскости 1-3 до изменения нагрузки; в - изменение нагрузки в сечении костной ткани

в плоскости 1-3

Рис. 4. Механизм изменения напряженно-деформированного состояния и перестройки трабекулярной костной ткани в главных осях тензоров 5, s, K (иллюстрация трабекулярной структуры взята из работы [8]): а — начальное состояние гомеостаза; б — начало перестройки костной ткани; в — процесс перестройки костной ткани; г — новое

состояние равновесия

В работе [8] тензор деформации s1, исходя из теории lazy zone, считается

~1 1

известной величиной. Характеристики архитектуры губчатой костной ткани Кие могут быть определены из соотношения (1) с учётом напряжённо-деформированного

состояния (a1, s'1).

Таким образом, рассматривается два состояния гомеостаза: первое состояние

при t < 0, описываемое напряжённо-деформированным состоянием (50, s0) и

архитектурой костной ткани (КГ0, е0). Новое состояние гомеостаза, наступившее по истечении достаточно большого промежутка времени с момента начала перестройки,

описывается напряжённо-деформированным состоянием (a1, s1) и архитектурой костной ткани (КГ1, е1). При этом тензор текущей деформации s (t) и текущие параметры структуры: девиатор тензора структуры K(t ) и доля твёрдого объёма кости e(t) описывают, соответственно, деформированное состояние и структуру в момент времени t (t > 0) между двумя состояниями физиологического равновесия.

Пусть P(s (t)) обозначает систему координат, связанную с главными осями

тензора s (t), а Р( К (t)) систему координат, связанную с главными осями тензора К (t) (рис. 4). Тогда если Z(t) — угол в плоскости 1—3 между системами координат P(s (t)) и

Р(50), то текущее деформированное состояние в момент времени t s (t) в неизменной с течением времени системе координат Р(о°) может быть задано как 62 ISSN 1812-5123. Российский журнал биомеханики. 2012. Т. 16, № 4 (58): 53—72

T

a (t )Р<5 >=ек(')) • г (t) Рда»'(е(С(')))!

где ортогональное преобразование Q(Z(t)) записано следующим образом:

С cos Z(t) 0 - sin Z(t) ^

(21)

Q(Z(t)) =

О 1 О sin Z(t) О cos Z(t)

(22)

и

є (t)

P(~(0) -

О є

2

О

О О є

(23)

Аналогично, если K(t) - угол в плоскости 1-3 между системами координат P(K(t)) и P(5°), то текущий девиатор тензора структуры в момент времени t K(t) в неизменной с течением времени системе координат P(o° ) может быть представлен как

K(t)P<a°> - б(к(і)) • K(t)nK(>» • ((£?(к(/)))T, (24)

где ортогональное преобразование Q(K(t)) задано как

С cos K(t) О - sin K(t) ^

(25)

и

Q(K(t)) - О 1 О

v sin K(t) О cos K(t) у

Г Ki О О Л

K(t)P(K(t)) - О K2 О

О V О - K1 - K2 ,

(26)

Опишем процесс перестройки, показанный на рис. 4. Рис. 4, а отражает начальное состояние гомеостаза, существовавшее при і < 0. На рисунке видно, что главные направления тензоров 5°, ~0 и К0 совпадают (т.е. их скалярные

произведения коммутативны).

На рис. 4, б показана ситуация, соответствующая началу перестройки костной

ткани (і = 0), когда начальное напряжённое состояние костной структуры 5° изменилось на а1. На рисунке показан неизменный с течением времени тензор а1, чьи главные оси Р (51) образуют угол т в плоскости 1-3 с главными осями Р (50) тензора

первоначального напряжённого состояния ~0 . Также показано новое значение тензора деформации ~(і), в то время как девиатор тензора структуры сохранил своё

первоначальное значение К0, поскольку трабекулярная микроструктура не успела отреагировать на новое напряжённое состояние.

Рис. 4, в отражает процесс перестройки костной ткани при V > 0. Угол между системами координат Р(~ (V)) и Р(а0) обозначен как £(/), угол между системами

координат Р(К(V)) и Р(а0) - как к(^). Главные направления тензоров ~(V) и К(V) перемещаются таким образом, чтобы стать соосными главным направлениям нового

напряжённого состояния а1.

Наконец, на рис. 4, г показано новое состояние равновесия, которое достигается по прошествии достаточно большого промежутка времени. Новое гомеостатическое

состояние может быть описано тензорами напряжений а1, деформации ~1 и девиатора тензора структуры К1 с главными осями Р(а1), Р(а1) и Р(К1) соответственно, а также новой долей твёрдого объёма кости V1.

Далее, с учётом введённых обозначений, тензорные соотношения (1), (2), (3) могут быть преобразованы в систему скалярных уравнений. Подстановка выражений (19), (21) и (24) в (1), (2), (3) приведёт к системе, состоящей из восьми уравнений с восемью неизвестными: 8^), 82(0, 8з(0, £(0, К^О, К2(0, к(0 и в(£).

В дальнейшем все расчеты будут проводиться в системе отсчёта Р(К(V)). Тогда из рис. 4 видно, что угол в плоскости 1-3 между главными направлениями тензоров К(V) и а1, т.е. между осями систем координат Р(К(V)) и Р(а1), равен (к-т).

Аналогично угол в плоскости 1-3 между главными направлениями тензоров К(V) и 8(V), т.е. между осями систем координат Р(К(V)) и Р(~(V)), равен (к-£). Таким образом, в системе координат Р(К(V)) имеем следующие компоненты тензоров а1, 8 (V) и К~(^) :

( 11 1 а1 + а3 + а1

-1Р (К (/)) _

а3

соб(2к - 2т)

а1 -а3

б1п(2к - 2т)

0

а2

0

а1

а3

б1п(2к - 2т) 0

1

1

+ а3^^1со8(2к-2т)

2 2

, (27)

Р(К(0) _

+ ^!^сов(2к-20 0

81 -83

б1п(2к - 20

82

0

б1п(2к- 20 0

+ ^з^сов(2к-2С)

(28)

Получим скалярные уравнения, соответствующие соотношению (1). Обозначим правую часть (1) в системе координат Р(К(V)) как тензор А, а именно:

1А 0 А

А(оР(К(()) _ 0 А22 0

V А13 0 А33 у

(29)

Запишем соотношение (1) в компонентной форме, причём правую часть будем представлять согласно введённым обозначениям (29). Тогда получим следующие соотношения:

а1 + °3 - А11 + А33>

(01 — 03) соб(2к — 2т) — Ап — А33.

а2 _ А22 ,

(О - а3)в1и(2к - 2т) _ 2А13.

Преобразуем второе и четвёртое уравнения системы (30) таким образом, что

О -03 _ 2А13 бш(2к - 2т) + (Ац - А33)соб(2к - 2т),

0 _ 2А13 соб(2к - 2т) - (Ац - А33) бш(2к - 2т).

(30)

(31)

Для дальнейшей работы с правой частью соотношения (1) введём некоторые обозначения. Пусть тензор ~(V) в системе координат Р(К(V)) имеет следующие компоненты:

/ а 0 ё

йе!

є (і) —

(32)

а компоненты тензора К(ї) в этой же системе координат таковы:

/а 0 0^| 0 в 0

йе!

К(ї) —

0 0

(33)

Тогда компоненты тензора А с учётом соотношений (1), (32), (33) можно записать в виде

А11 _ (§1 + §2е)(а + ь + с) + (§3 + §4е)а + 2§5аа + §6 (2аа + Ь(а + в) + с(а + У)), (34)

А22 _ (§1 + §2е)(а + Ь + с) + (§3 + §4е)Ь + 2§5ЬР + §6 (2Ьв + а(а + в) + с(Р + У)), (3 5)

А33 _ (§1 + §2е)(а + Ь + с) + (§3 + §4е)с + 2§5сУ + §6(2сУ + а(а + У) + Ь(в + У)), (36)

А13 — (83 + 84е)й + 85й (а + У)-

(37)

Подставим соотношения (34)-(37) в правые части систем (30) и (31) и получим следующие выражения:

_1 _1

а1 - а3 _ (2(§3 + §4е)ё + 2§5ё(а + у)) Бт(2к - 2т) +

+ ((§3 + §4е)(а - с) + 2 §5(аа - сУ) + §6(а- У)( а + ь + с)) со§(2к - 2т),

_ 2(§1 + §2е)(а + ь + с) + (§3 + §4е)(а + с) + 2§5(аа + сУ) +

+ §6 (а(3а + у) + Ь(а + 2в + у) + с(а + 3у)),

а2 _ (§1 + §2е)(а + ь + с) +(§3 + §4е)Ь + 2§5ЬР + §6 (2ЬР + а(а + в) + с(Р + У)) =

0 _ (2(§3 + §4е)ё + 2§5ё(а + у)) соб(2к - 2т) -- ((§3 + § 4е)(а - с) + 2 §5(аа - сУ) + §6(а - У)(а + ь + с)) зт(2к - 2т).

!88М 1812-5123. Российский журнал биомеханики. 2012. Т. 16, № 4 (58): 53-72

(38)

(39)

(40)

(41)

Далее соотнесём компоненты тензора ~ (ґ) (см. выражение (28))

с аналогичными, введёнными согласно выражению (32), компонентами. Также соотнесём компоненты тензора К (ґ) (см. выражение (26)) с аналогичными

компонентами тензора, введённого согласно выражению (33). Воспользуемся полученными результатами, выразим соответствующие слагаемые из соотношений (38)—(41) и преобразуем их. В результате получается следующая система нелинейных алгебраических уравнений:

°1 - = (£з + - g5K2 ) (S1 - S3 ) COs 2(T-Z) +

(42)

(43)

+ [ + gб)(S1 + S3) + gбS2 ](2K1 + A^^2^-*^

°1 + °3 = 2(gi +eg2)(Si +S2 +S3) + (g3 + eg4 (g 5 + 2g6) K2)(S1 +S3) +

+ (g5 + gб ) (S1 - S3 ) (2K1 + K2 ) COs 2(Z - к) + gбS2K2 ,

^2 = (gi + eg2 ) (S1 + S2 + S3) + (g3 + eg4 ) S 2 + 2g5S2K2 +

Г 1 ^ 1 (44)

+gбK2 I 2S2 + 2 (S1 + S3) I + 2 gб(S1 - S3) (2K1 + K2 ) COs 2(Z - к),

(g3 + eg4 - g5K2)(S1 - s3)sin 2(C - T) +

+((g5 + gб)(S1 + S3) + gбS2 ) (2K1 + K2)sin 2(к - T) = 0

(45)

Данная система практически полностью соответствует соотношениям (33), представленным в [B] (однако в работе [B] в соотношениях был обнаружен ряд неточностей).

Запишем кинетическое уравнение (2) в компонентной форме. При этом будем учитывать, что производится дифференцирование тензорного соотношения в подвижной системе координат [2] и необходимо вычислять коротационную производную Яумана. Если

K(t) P(s0) = В P(~0), (4б)

где В P(~~ ) - правая часть уравнения (2) в неподвижной системе координат P(50), то при переходе в подвижную систему координат P(K(t)) получаем, что

K(t)P(K(t)) = ВP(K(t)) - Q • K(t)P(K~(t)) + K(t)P(K~(t)) • Q, (47)

где ВP(K(t)) - правая часть уравнения (2) в подвижной системе координат P(K(t)), имеющая следующий вид:

B(t)P(K (t)) =

При этом

сц4 0 В1З

0 В22 0

V В1З 0 - Cq4 -

(4B)

В

P(K(t)) = (t))) • ВP(~~0) • (~(к(t)) . (49)

Здесь Q определяется как

Q = U • U

T

(50)

T

где и = (бОФ)))5 - тензор перехода от неподвижной системы координат к подвижной. С учётом выражения (25) определяем, что

Q =

00 к 0

0

0

(51)

Также с учётом соотношения (26) находим, что:

(

K (t)

P(K(t)) _

Ktf) 0

0 K 2(t)

00

0

0

- K:1(t) - k2(t)

(52)

Далее определим правую часть уравнения (2) в компонентной форме

в подвижной системе координат Р(К (ї)) (т.е. определим ВР(К (ї))) и вычислим правую часть соотношения (47). При этом, учитывая, что компоненты из (52) будут равны соответствующим компонентам в правой части соотношения (48), получим

К = Вп,

(53)

K 2 B22,

1C = -

B„

2 K1 + K 2

Если принять, что ~0 = ~1, т.е. считать, что начальное и новое деформированные состояния находятся в диапазоне, определяемом границами lazy zone, то тензор ~0 в системе отсчёта P( K (t)) примет следующий вид:

0 , „0 „0 „0

+ £l_-£lcos(2T- 2к) 0

0 P(K(t))

2

2

00 S1 -s3

0

sin(2x - 2к)

„2

„0 - 0 1 2 3 sin(2x - 2к)

0

0 0 0 0 81 + „3- + ^l^^cos(2T- 2к)

2

2

. (54)

Введём обозначение, аналогичное (32), а именно:

~0 def „ °(<) =

' a 0 0 d0'

0 b О 0

d О 0 с0

(55)

Тогда с учётом эволюционного уравнения (2) и введённых обозначений (32), (33) и (55) получим следующие компоненты тензора В. При этом следует учитывать, что у = -а - в (поскольку 1хК = 0 ).

0

2

В11 = 1(И1 + екъ)(2(а - а0) - (Ъ - Ъ0) - (с - с0)) + И4а((а - а0) + (Ъ - Ъ0) +

3 (56)

+ (с - с0)) - К, (2(а - а0)а - (Ъ - Ъ0)в - (с - с0)у),

В22 = 3(А+ еК)(2(Ъ -Ъ0) - (а-а0) - (с - с0)) + К4а((а-а0) + (Ъ - Ъ0) + (57)

+ (с - с0)) - И2(2(Ъ - Ъ0)Р - (а - а0)а - (с - с0)у),

В13 = 2(ё - &0)(2(к1 + е^3) - 3^2(а + У)).

(58)

Далее в соответствии с равенствами (26), (28) и (54) можно конкретизировать компоненты тензора В (эти вычисления аналогичны ранее представленным - см. (30)-(41)). Окончательно получим систему трёх нелинейных дифференциальных уравнений:

К = 1 & 2

1

+ —

2

\^{И1+ек3 +И2(К2 -К1))( ^ -83)со82(<^ -к)-(в^° -в!°)со8 2(т-к)) J-

^ К2(К2 + К1) - ~3(К + еК3) ^((28 2 -81 -В3) - (20 2 - £0 -80))

+ (К4 - Ю К1 (81 -80) + (82 -82) + (83 -80))

(59)

&

(3 (А + еК,)- 2ККг]( (( -8, -(28° -80 -ф ) +

+

к

~2^(2К1 + К2)( -83)сов2(^-к)- (8° -8°!)СОБ2(т- к)) +

+ (К4 - К2) К2 ( (81 +82 +83) - (810 +82 + 83) ),

(60)

=2((4( ^2 (- 83)81п2(с- к) - (80 - 80)^2(т - к))- (61)

Запишем кинетическое уравнение (3) в компонентной форме. Подробный вывод представлен не будет, поскольку все выкладки аналогичны предыдущим. Отметим только, что с учётом введённых обозначений (32), (33) и (55) кинетическое уравнение (3) можно записать в следующем виде:

е = С/1 + е/2)((а - а°) + (Ъ - Ъ°) + (с - с0)) + /3((а - а0)а + (Ъ - ъ0)Р + (с - с0)т). (62)

В итоге, учитывая равенства (26) (28) и (58), получим дифференциальное уравнение вида:

ёе

' (/ + /2е) ( (81 +82 +83) - (810 +8 2 +80) ) +

/ ( (^1 +К2 ( ((81 - 83 ) СОБ 2(С - к) - (8!0 - 80 ) СОБ 2(т - к)) ) -

(63)

+ /3 К2 ((282 -81 -83)-(2822 -80 -80) ).

Полученная система скалярных уравнений включает в себя четыре нелинейных алгебраических (42)-(45) и четыре нелинейных дифференциальных уравнения (59)-(61) и (63). Данная система описывает перестройку трабекулярной костной ткани из первоначального гомеостатического состояния к новому гомеостатическому состоянию и не может быть линеаризована по вышеназванным причинам. Отметим, что в работе [8] в аналогичных соотношениях был обнаружен ряд неточностей.

Далее будет показано численное решение системы уравнений (42)—(45), (59)-(61) и (63). Для этого воспользуемся примером, ранее частично описанным в работе [8] (пример 3). Рассматривается модель трабекулярной микроструктуры (см. рис. 3), находящаяся в состоянии физиологического равновесия в течение достаточно долгого

промежутка времени ( < 0). Начальное напряжённое состояние а0 задаётся в виде следующего тензора напряжений:

ї0 P(~0) =

I-1,4024 0 0

0 0 ^

-1,4024 0

0 -1,6859

ґ0

J

Начальное деформированное состояние s задаётся исходя из теории lazy zone. А именно, было известно, что перестройка трабекулярной костной ткани отсутствует, если начальное деформированное состояние задать в следующем виде:

г0 Р(а0) =

I- 0,0015 0 0

0 0 ^

- 0,0015 0

0 - 0,0015 J

Система уравнений (42)-(45) позволяет численно определить начальные параметры архитектуры губчатой костной ткани К 0 и е0 (следует отметить, что для нахождения К0 и е0 также можно воспользоваться системой (11)—(13), поскольку до изменения напряженного состояния тензоры ~0, ~0 и К0 были соосными).

В результате получим, что е0 = -0,01,

I- 0,01 0

K

0 P(а0) =

0

- 0,01 0

0 0,02

J

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

Новое напряжённое состояние а1 определяется исходя из соотношений (7) и (8). В этом случае угол между главными направлениями тензоров а0 и а1 в плоскости 1-3 т = 30°. Тензор а1 в своих главных осях В(а1) представляется в виде

А Р(а1) =

I-1,5472 0 0

0 0 >

-1,5472 0

0 - 2,1141J

При этом в начальный момент времени (ї = 0) главные оси тензора деформаций в должны мгновенно повернуться на угол (^° = 48° (см. систему уравнений (42)-(45)).

0

Кг

-0,009

-0,010

-0,011

-0,012

-0,013

0,014

-0,015

-0,016

-0,017

-0,018

-0,019

0,020

0 20 40

60 80 /, сут

а

100 120 140 160

80 I, сут б

К,

/, сут в

и сут г

Рис. 5. Изменение структуры (компонент девиатора тензора структуры К1(() (а), К2(/) (б)), плотности (доли твёрдого объёма кости в(/)) трабекулярной микроструктуры (в) и угла к(/)

между главными осями тензоров 5° и К(() (в плоскости 1-3) (г)

Результаты, полученные по истечении 160 суток, показаны на рис. 5. Видно, что новая костная микроструктура описывается следующими величинами: в1 = 0,01,

(- 0,0015

г1 Да1) _

0

К

1 Да1) _

0 - 0,0015

00 (- 0,02 0 0 - 0,02 00

0 0

- 0,0015 0 ^

0 0,04 1

Л

При этом видно, что т = к(^1), т.е. главные оси тензоров 51 и К1 вновь стали соосными. Отметим, что полученные результаты совпадают с аналогичными из [8].

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

Заключение

На основе ранее полученных определяющих и эволюционных соотношений, способных описать перестройку трабекулярной костной ткани, происходящую вследствие изменяющихся условий нагружения, рассмотрен ряд модельных примеров. Для определения коэффициентов, входящих в кинетические уравнения (2) и (3), был проведён вычислительный эксперимент, где учитывалось, что адаптация костной ткани должна происходить за 160 дней. Результаты показывают различный характер влияния изменения нагрузки на процесс формирования структуры, подчиняющийся закону Вольфа.

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

Работа выполнена при частичной поддержке РФФИ (грант 11-01-00910-а).

Список литературы

1. Демидов С.П. Теория упругости. - М.: Высш. шк., 1979. - 432 с.

2. Димитриенко Ю.И. Тензорное исчисление. - М.: Высш. шк., 2001. - 575 с.

3. Киченко А.А., Тверье В.М., Няшин Ю.И., Симановская Е.Ю., Еловикова А.Н. Становление и развитие классической теории описания структуры костной ткани // Российский журнал биомеханики. - 2008. - Т. 12, № 1. - С. 68-88.

4. Киченко А.А., Тверье В.М., Няшин Ю.И., Заборских А.А. Экспериментальное определение тензора структуры трабекулярной костной ткани // Российский журнал биомеханики. - 2011. - Т. 15, № 4. -С. 78-93.

5. Тверье В.М., Симановская Е.Ю., Еловикова А.Н., Няшин Ю.И., Киченко А.А. Биомеханическое описание структуры костных тканей зубочелюстной системы человека // Российский журнал биомеханики. - 2007. - Т. 11, № 1. - С. 9-24.

6. Тверье В.М., Симановская Е.Ю., Няшин Ю.И., Киченко А.А. Биомеханический анализ развития и функционирования зубочелюстной системы человека // Российский журнал биомеханики. - 2007. -Т. 11, № 4. - С. 84-104.

7. Экспериментальные методы в биомеханике / под ред. Ю.И. Няшина, Р.М. Подгайца. - Пермь: Изд-во ПГТУ, 2008. - 400 с.

8. Cowin S.C. An evolutionary Wolff’s law for trabecular architecture // J. Biomech. Engineering. - 1992. -Vol. 114. - P. 129-136.

9. Cowin S.C. Bone mechanics handbook. - Second ed. - New York: CRC Press, 2001. - 1136 p.

10. Cowin S.C. Fabric dependence of an anisotropic strength criterion // J. Mech. Materials. - 1986. - Vol. 5. -P. 251-260.

11. Cowin S.C. The relationship between the elasticity tensor and the fabric tensor // J. Mech. Materials. -

1985. - Vol. 4. - P. 137-147.

12. Cowin S.C. Wolff’s law of trabecular architecture at remodeling equilibrium // J. Biomech. Engineering. -

1986. - Vol. 108. - P. 83-88.

13. Co win S.C., Mehrabadi M.M. Identification of the elastic symmetry of bone and other materials // J. Biomechanics. - 1989. - Vol. 22. - P. 503-515.

14. Fritton S.P. Computational simulation of trabecular bone adaptation: Ph.D. dissertation. - New Orlean: Tulane University, Departament of Biomechanical Engeenering, 1994.

15. Fung Y.C. Biomechanics. - New York: Springer-Verlag, 1990. - 464 p.

16. Martin R.B., Burr D.B., Sharkey N.A. Skeletal tissue mechanics. Second edition. - New York: Springer-Verlag, 1998. - 392 p.

17. Telega J.J., Jemiolo S. Fabric tensor in bone mechanics // J. Engineering Transactions. - 1998. - Vol. 46. -P. 3-26.

18. 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.

19. 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.

20. Wolff J. Das gesetz der transformation der knochen. - Berlin: Hirshwald, 1892. - 416 s.

ON APPLICATION OF THE THEORY OF TRABECULAR BONE TISSUE REMODELLING

A.A. Kichenko, V.M. Tverier, Y.I. Nyashin, M.A. Osipenko, V.A. Lokhov (Perm, Russia)

The constitutive relations which allow us to describe the stress-strain state of the cancellous bone tissue taking into account its structure and evolution relationships describing the adaptation processes in the bone are shown. Based on authors’ previously stated initial boundary value problem which includes Cowin’s constitutive relations and evolution relationships, the solution algorithm is developed and trabecular bone tissue evolution is demonstrated on the set of model examples when the stress-strain state is changed. The computing experiment is carried out to determine coefficients in evolution equations; it was considered that bone tissue adaptation should take place for 160 days. The results demonstrate different character of influence of changes of loading conditions on process of structure formation which follows from the Wolff’s law.

Key words: biomechanical modelling, 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