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

DOI: 10.15593/RZhBiomeh/2018.3.01 УДК 531/534: [57+61]

БИОМЕХАНИЧЕСКОЕ МОДЕЛИРОВАНИЕ ТРАБЕКУЛЯРНОЙ КОСТНОЙ ТКАНИ В СОСТОЯНИИ РАВНОВЕСИЯ

Т.Н. Чикова, А.А. Киченко, В.М. Тверье, Ю.И. Няшин

Пермский национальный исследовательский политехнический университет, Россия, 614990, Пермь, Комсомольский проспект, 29, e-mail: ChikovaTN@gmail.com, kichenko.alex@inbox.ru, tverier_55@perm.ru, nyashin@inbox.ru

Аннотация. В XXI веке J. Wolff отметил, что кость здорового человека или животного адаптируется к тем нагрузкам, которым подвергается. На уровне ткани в кости просматриваются слабопористые участки компактного вещества и сильнопористые -трабекулярного (губчатого). Известно, что в трабекулярной костной ткани механизм адаптации реализуется посредством выстраивания трабекул (костных балочек, образующих твердый матрикс) вдоль линии действия преобладающей нагрузки. При достижении трабекулярной тканью оптимальной структуры для существующего в локальной области нагружения кость переходит в состояние равновесия (гомеостаза). В своих работах S. Cowin предложил описывать положение трабекул в каждый дискретный момент времени главными направлениями тензора структуры, отыскиваемыми из решения системы кинетических уравнений. На основе предложенных соотношений авторами решены многие задачи динамики трабекулярной костной ткани, результатом которых являются функции некоторой величины, зависящей от времени. Например, рассмотрено изменение пористости, компоненты девиатора тензора структуры. В данной работе кинетические уравнения используются с целью определения структуры, напряженного состояния или упругих свойств костного образца, находящегося в состоянии равновесия. Математическая постановка представлена задачей теории упругости анизотропного тела. Все численные расчеты выполнены с использованием программного продукта ANSYS Mechanical APDL на примере растяжения прямоугольной пластины. Проведено сравнение результатов задачи, полученных аналитическим методом и методом конечных элементов. Предполагается, что полученные соотношения могут быть полезны для определения реализуемого в кости известной структуры напряженно-деформированного состояния или, наоборот, для предсказания оптимальной структуры кости при заданных условиях нагружения.

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

Введение

Рассмотрим стационарную задачу о плоской деформации образца из трабекулярной костной ткани как задачу линейной теории упругости анизотропного тела, дополненную кинетическими уравнениями, описывающими перестройку костной ткани. Математическая постановка задачи состоит из 21 скалярного уравнения [3]:

© Чикова Т.Н., Киченко А.А., Тверье В.М., Няшин Ю.И., 2018

Чикова Татьяна Николаевна, аспирант кафедры теоретической механики и биомеханики, Пермь

Киченко Александр Александрович, старший преподаватель кафедры теоретической механики и биомеханики,

Пермь

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

уравнения равновесия:

V-д = 0, х еУ, г >0, (1)

определяющего и геометрического соотношений: д = (Я + Я2е) (1г а)ё + (Яэ + Я4е) е + Яз (е - к + к - е) + Яб (* (к - а)ё + (1г а)к), (21)

(2.2)

|уи+иу |,

эволюционного уравнения, описывающего изменение ориентации трабекул в

е =1 (Уу и?+иУУ)

авнени

рассматриваемой области:

к = (h + vxe-e°) + htr (ё-Ё0 ) K +

+h[ (tr(k• (e-e0)))E-3(k• (e-e0)+(e-e0)• k)

(3)

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

е = 01 + /2е) (* е - 1г 8°) + / (1г (к - (е - е°))), (4)

граничных условий:

п - д = Р(г), х = ^, г > 0, (5)

и - и(г), х = яи, г > 0,

начальных условий:

к = к°, е = е0, х еУ, г = 0, (6)

Р = Р0, х е Я, и = и , х е Я, г = 0,

' а ^ 0 э и ? ?

где / — /, ^ — - константы [3, 4], сут-1; ^ — я6 - константы [8, 9], ГПа; К - девиатор тензора структуры (отвечает за направление трабекул); е - изменение доли твердого объема (отвечает за пористость).

Интересно посмотреть на поведение трабекулярной костной ткани под нагрузкой в момент времени, когда ее структура находится в состоянии равновесия. Для примера кость моделируется как прямоугольная пластина, нагруженная с одного торца растягивающей нагрузкой, а с другого конца пластина жестко закреплена (рис. 1).

Рис. 1. Схематичное представление модели

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

1 тт т> о Б0 0 " ~0 ~0

1. Даны К , Р , е , найти с ,8 .

о тт ~0 ~0 ~ Л- 0 0

2. Даны с ,8 , найти К , е .

3. Даны К° = 0, у = 1...3 (изотропия), с0,80, найти Е,V.

Решение стационарных задач

Определение напряженно-деформированного состояния

Решение первой задачи является ответом на вопрос: какое напряженно-деформированное состояние реализуется в кости определенной микроструктуры под заданной нагрузкой?

Определяющее соотношение (2.1) представляет собой форму записи обобщенного закона Гука с = С •• 8, в котором тензор упругих констант является симметричным тензором четвертого ранга, зависящим от структуры материала следующим образом (используется нотация Фойгта):

Си = (а + &е) + (& + §4е) + 2Ки (& + §), Си = (§ + ^2е) + § (Ки + К22),

С13 = (§1 + §2е) - §6 С14 = К23 ^ С15 = К13( §5 + С16 = К12(§5 + §6 X

С22 = (§1 + §2е) + (§3 + §4е) + 2К22 (§5 + §6 ), С23 = (§1 + §2е) - §6К11 , С24 = К23(§5 + С25 = К13 §6 , С26 = К12 (§5 + §6 X

С33 = (§1 + &2е) + (§3 + §4е) - 2(К11 + К22 )(§5 + §6 X С34 = К23 (§5 + ) ,

(7)

С35 = К23(§5 + С36 = К12^ С44 = I ((§3 + §4е) - §5K11),

С44 = "2 ((§3 ^ &4е) — §5К11), С45 = "2 К12§5 , С46 = "2 К13§5 , С46 = "2 К13§5 , С55 = 1 ((§3 + §4е) — §5К22 ), С56 = IК23§5 , С66 = 1 ((§3 + §4е) + §5 (К11 + К22 )).

В ЛЫ8У8 уравнения (7) записываются в таблицу свойств материала командой ТБПЛТЛ. Напряженно-деформированное состояние определяется в результате решения задачи методом конечных элементов с учетом внешних сил и ограничений на перемещения (рис. 2).

Для решения задачи в ЛЫ8У8 используем конечный элемент Р1апе182.

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

Допустим, что реализуется плоское напряженное состояние (а « тХ2 « т « 0; уХ2 «у « 0 ). Тогда благодаря принципу Сен-Венана можно заменить

влияние заделки неким напряжением ах, действующим на поперечное сечение А = Ь

(см. рис 1). Из уравнения равновесия для оси х (ахЬ — Р = 0) находится ах = Р / Ь.

Так как касательные и нормальные напряжения, действующие на площадках, перпендикулярных оси у, отсутствуют, то а « 0, т « 0.

Л PRNS0L Command X

File

NODE sx sv sz sxy W sxz Л

7 D.1DDD3E+0D6 7.3651 0.0000 -2.1388 0.0000 0.0000

21 D.1DDD3E+0D6 7.3651 0.0000 2.1388 0.0000 0.0000

11 99996. 9.3095 0.0000 -3.3817 0.0000 0.0000

12 99983. 11.613 0.0000 0.16555E-007 0.0000 0.0000

13 99996. 9.3095 0.0000 3.3817 0.0000 0.0000 V

Л PRNSQL Command X

File

NODE EPELK EPELV EPELZ epelxy epelyz epelxz Л

7 0.16108E-003-0.27029E-001-0.27013E-001-0.91731E-008 0.0000 0.0000

Il 0.16108E-003-0.27029E-001-0.27013E-001 0.91731E-008 0.0000 0.0000

Al 0.16103E-003-0.27017E-001-0.27035E-001-0.12720E-007 0.0000 0.0000

Al 0.16101E-003-0.27005E-001-0.27033E-001 0.16368E-015 0.0000 0.0000

13 0.16103E-003-0.27017E-001-0.27035E-001 0.12720E-007 0.0000 0.0000 V

Рис. 2. Компоненты тензора напряжений и деформации для сечения х ~ I /2 в ЛЫБУБ

Оставшиеся 4 неизвестных (ех, еу, е2, у^) находятся из четырех уравнений закона

Гука

о х = Сие х + С12е у + Спг 2, (8)

с = C 2 8 + C22 8 + C2, 8 ,

y 12 x 22 y 23 z5

0 = C 8 + C 8 + C 8

0 C138x ^ C328y ^ C338z , TХУ _ C668xy .

Тогда

8 z =

(C13C22 C13C12C23 )

—Гэ Г + 1Г2Г Г —Г Г2 Г —Г2г г +г г г г

С13С22 ^ 2С13С12С23 С11С23С13 С12С33С13 ^ С11С33С13С22

С С — С С С С (9)

е = ч2^33 ч^н е , е =— ч3. е — чз8 , 8 = 0.

С13С22 — С12С23 С13 С13

Тот же результат получается из четырех уравнений, записанных через матрицу податливости,

8х = 0х ■ С1Г1, 8у = °у ■ С2Г1, О ■ С3Г1, Тху = 0. (10)

С учетом (1) и (2.1), а также зная, что компоненты тензора К равны нулю при Р = 1 кН (ах = 10 000 Па), получаются следующие значения:

вх = 0,16104 • 10-4, в, = вг -0,27034 -10"5, yw = 0 . (И)

Аналитическое решение (11) совпадает с численным с точностью А = 10—8. Связь компонент тензора малой деформации с перемещениями определяется формулами Коши:

_ди _ду _ Эю

= дх' 8уу =дУ' 8~ = &" • (12)

д

Рис. 3. Осевые перемещения (а, в, д) и относительная погрешность (б, г, е) численного решения в сечении х «I /2 (м), К = 0, г, ] = 1. 3

б

а

в

г

е

Определяя перемещения путем интегрирования, получим и = ех х + с, V = е^ у + с2, ^ = е2 2 + с. Постоянные и, V, ^ определяются из условий закрепления

Нх=0 = ^у=0 = Ч*=0 = 0 ^ ^ С3 = 0 .

Удовлетворяя им, получим

и = ехх, V = е у, w = е2 г.

(13)

(14)

(15)

Рис. 4. Изолинии перемещений по оси х при: а - К = 0, г,у = 1...3; б - Кп = —0,09, К2 = —0,05, К33 = 0,14 ; в - Кп = —0,09, К22 = —0,05, К33 = 0,14, а = 64°

а

б

в

Рассмотрим перемещения в поперечном сечении с координатой х = I /2 (рис. 3). На всех трех графиках относительная погрешность меньше 1%.

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

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

Рассмотрим контурные графики перемещений при задании различных компонент девиатора тензора структуры, в частности при его повороте на угол а (см. рис. 4, в).

-6.932 -6,94 их -6,948 -6,956 -6,954

•10

-6,972

0,4

►а н о о X

0 0,32 и

о 0,24 С

■10"

-6,25 -3,75 -1,25 1.25 3,75 6,25 У

К

Л

Ц

и н

С о я

н О

0,16 0.08 0

•ю-

-6.25 -3,75 -1.25 1,25 3,75 6,25 У

б

д е

Рис. 5. Осевые перемещения (а, в, д) и относительная погрешность (б, г, е) численного решения в сечении х «I /2 (м), Кп = —0,09, К22 = —0,05, К3з = 0,14

а

Аналитическое решение (15) говорит о том, что в общем случае анизотропии при растяжении стержень не только удлиняется в направлении силы и сокращается в поперечных направлениях, но еще и испытывает сдвиги во всех плоскостях, параллельных координатным [5]. Однако используемые в [5] формулы для деформации и перемещений в случае растяжения анизотропного стержня нельзя применять в данной задаче из-за сделанных изначально предположений (в частности, из-за предположения о плоском напряженном состоянии).

Рассмотрим поперечные сечения при анизотропии (К ^ 0). В качестве аналитического решения используем соотношение (15).

Из графиков на рис. 5 видно, что перемещения ив значительной степени зависят от

значений компонент тензора К, поэтому в случае анизотропии соотношение V = е у применять неверно.

Отыскание параметров структуры

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

Исходные данные возьмем из решения предыдущей задачи, а именно значения компонент тензора напряжений и деформации. Требуется доказать, что К = 0 при

/, j = 1...3 и е = —0,0179. Для этого определяющее соотношение (2.1) преобразуется в систему

{в, Кц, K22, K23, Ki3, Ki2 f = A J 8x , 8y , 8 z , 1 8yz , 1 8 xz , 1 8xy

где [ A ] =

6-6

g2tr ë + ^48x g2tr ë + g48y

g2tr ë + g48z g48 yz g48 xz g48 xy

2( g 5 + g6 )8 x + g68 y

g6 (8x -8 z )

-2(g5 + g6 )8z-g68y

—g58 yz 0

g58

xy

g6 (8 y -8 z )

2(g5 + g6 )8y + g68x

-2(g5 + g6 )8 z-g68x 0

- g58 xz

g5 8 xy

(16)

2g68yz 2(g5 + g6) 8yz 2(g5 + g6) 8yz (8y + 8z )(g5 + g6 ) + g6 8x

g 5 8 xy

g5 8 xz

2( g 5 + g6 )e.

2 g6 8 xz

2(g5 + g6 )e.

g5 8 xz

(8x + 8z )(g5 + g6 ) + g5 8

g5 8 yz

5 xy

2(g5 + g6 )8xy 2(g5 +g6 )8xy 2 g6 8 xy

g5 8 xy g5 8 yz

(8x + 8y )(g5 + g6 ) + g5 8z

Тогда е = —0,0179; Кп = 0; К22 = —1,3878-10—1'; К13 = 0; К23 = —1,3878-10—1'; К12 = —1,3878-10 Точность решения А = 10—1б.

T

Отыскание механических свойств материала

Задание к в определенном виде определяет тип упругой симметрии, реализуемый в материале. Когда К = 0, используя физическое соотношение для изотропного тела, можно найти технические константы Е, v.

E (1 - v)

-е„„ +-

Ev

е * +

Ev

(1 + v)(1 - 2v) " (1 + v)(1 - 2v) ^ (1 + v)(1 - 2v)

е = g ,

zz xx?

Ev

-е„„ +■

E (1 - v)

Ev

е H--е = g ,

yy /1 , о,Л zz yy'

(1 + v)(1 - 2v) xx (1 + v)(1 - 2v) yy (1 + v)(1 - 2v)

E

(1Hv)

е = t .

xy xy

Тогда модуль Юнга и коэффициент Пуассона находятся по формулам

(17)

v =

g е - е т

xx xy xx xy

E = -21(1 + v).

2gxxеxy - тxy (еxx + еyy + еzz ) е~

xy

(18)

Заключение

Опираясь на кинетические уравнения феноменологической теории адаптивной упругости, предложенные S. Cowin для описания перестройки архитектуры трабекулярной костной ткани, вызванной выходом локальной деформации кости за пределы lazy zone (диапазон деформаций, внутри которого не происходит изменений в структуре), авторы предложили уравнения, позволяющие анализировать костную ткань (напряженно-деформированное состояние, структуру и пористость, упругие свойства материала) в гомеостазе. Сравнение аналитических решений с результатами расчета в ANSYS показало работоспособность математической модели.

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

Работа поддержана грантом Российского фонда фундаментальных исследований № 18-01-00589 а.

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

1. Гороженинова Т.Н., Киченко А.А. Моделирование изгиба анизотропной консольной балки в ANSYS Mechanical // Master's Journal. - 2017. - № 1. - C. 225-229.

2. Гороженинова Т.Н., Киченко А.А. Создание интерфейса между ANSYS и MATLAB на примере перестройки трабекулярной костной ткани // Master's Journal. - 2018. - № 1. - С. 225-229.

3. Киченко А.А., Тверье В.М., Няшин Ю.И., Осипенко М.А., Лохов В.А. Постановка начально-краевой задачи о перестройке трабекулярной костной ткани // Российский журнал биомеханики. - 2012. - Т. 16, № 4. - С. 36-52.

4. Киченко А.А., Тверье В.М., Няшин Ю.И., Осипенко М.А. О приложении теории перестройки трабекулярной костной ткани // Российский журнал биомеханики. - 2012. - Т. 16, № 4. - С. 53-72.

5. Лехницкий С.Г. Теория упругости анизотропного тела - М.: Главн. ред. физ.-мат. лит. изд-ва «Наука», 1977. - 416 с.

6. Faina A.S., Gerasimov O.V., Sachenkov O.A. The evolution of the trabecular bone at a constant combined // International Journal of Pharmacy&Technology. - 2016. - Vol. 8, № 4. - P. 24261-24271.

7. Gerasimov O., Shigapova F., Konoplev Y., Sachenkov O. The evolution of the bone in the half-plane under the influence of external pressure // Mesh methods for boundary-value problems and applications: 11th International Conference. - 2016. - Ser. 158.

8. Cowin S.C. Wolff's law of trabecular architecture at remodeling equilibrium // J. Biomech. Engineering. - 1986. -Vol. 108. - P. 83-88.

9. Turner C.H., Cowin S.C., Rho J.Y., Ashman R.B., Rice J.C. The fabric dependence of the orthotopic elastic constants of cancellous bone // Journal of Biomechanics. - 1990. - Vol. 23. - P. 549-561.

10. Tverier V., Kichenko A., Nyashin Y., Lokhov V. Experimental construction of the fabric tensor for trabecular bone tissue // Series on Biomechanics. - 2015. - Vol. 29, № 4. - P. 33-38.

11. Wolff J. The law of bone remodelling. - Berlin, Germany: Hirshwald, 1986. - 126 p.

BIOMECHANICAL MODELLING OF TRABECULAR BONE TISSUE IN REMODELLING EQUILIBRIUM

T.N. Chikova, A.A. Kichenko, V.M. Tverier, Y.I. Nyashin (Perm, Russia)

In the 19th century, J. Wolff noted that the bone of a healthy person or animal adapts to the loads under which it is placed. At the tissue level, weakly porous regions of compact tissue and strongly porous regions of the cancellous (spongy) tissue are seen in the bone. It is known that in the cancellous bone tissue the adaptation mechanism is realized by aligning the trabeculae (bone bunches forming a solid matrix) along the principal stress trajectories. When the cancellous tissue reaches the optimal structure for the existing loading in the local area, the bone passes into a state of equilibrium (homeostasis). In his works, S. Cowin proposed to describe the position of the trabeculae at each discrete instant of time by the principal trajectories of the fabric tensor calculated from the solution of the system of rate equations. On the basis of the proposed relationships, many problems of the dynamics of cancellous bone tissue are solved, the results of which are functions of a certain quantity that depends on time. For example, the change in porosity, the components of the deviator of the fabric tensor. In this work, rate equations are used to determine the structure, stress state or elastic properties of a bone sample in equilibrium. The mathematical statement is presented by the anisotropic elasticity problem. All numerical calculations were performed using the ANSYS Mechanical APDL software on example of a rectangular plate tension. The results of the problem obtained by the analytical approach and the finite element method are compared. It is assumed that the obtained relationships can be useful for predicting of the stress-strain state in the bone of known structure, or vice versa, to predict the optimal bone structure for the given loading conditions.

Key words: cancellous (spongy) tissue, fabric tensor, equilibrium (homeostasis), boundary -value problems of biomechanics.

Получено 1 августа 2018