The effect of bone tissue density on the stress-strain state near dental implants

Full Text

Введение

Эффективность дентальной имплантации помимо прочих факторов [1] зависит от плотности костной ткани, которая может изменяться в широком диапазоне [2, 3]. Ввиду этого изучение влияния плотности кортикальной и губчатой костных тканей челюсти на напряженно-деформированное состояние в зоне дентальной имплантации представляет большой научный и практический интерес [4, 5]. Изменение плотности костной ткани проводит к изменению ее модуля упругости [6]. Этот факт позволяет выполнять математическое моделирование влияния плотности костной ткани на напряженно-деформированное состояние при квазистатическом приложении нагрузки посредством соответствующего изменения модуля упругости кости [3, 7–9, 11].

Для моделирования напряженного состояния стоматологических имплантатов и костных тканей в последние десятилетия используются численные методы. Наибольшее распространение получил метод конечных элементов (МКЭ) [12, 13]. Mетод граничных элементов (МГЭ) в задачах стоматологической биомеханики используется не так широко. Одно из первых применений этого метода было выполнено в 1993 году для анализа напряженного состояния соединения цилиндрического имплантата и костной ткани [14]. МГЭ также применялся для исследования напряженного состояния костных тканей и винтовых стоматологических имплантатов, изготовленных из различных материалов (титана или керамики) [15]. Сравнительный анализ применения МКЭ и МГЭ для анализа напряженного состояния имплантатов и костной ткани был выполнен в [16]. При использовании МГЭ ввиду независимого моделирования перемещений и напряжений на поверхности тела, в зонах с высокими градиентами напряжений можно получить результаты, сопоставимые с результатами, полученными с помощью МКЭ, но без существенного измельчения сетки элементов. При решении задачи МГЭ выполняется дискретизация только поверхности (границы — в двумерном случае) расчетной модели, что является одним из основных преимуществ данного метода перед МКЭ, при использовании которого необходимо разбиение на элементы всего объема тела. В данной работе МГЭ используется для анализа влияния изменения плотности костных тканей на напряженное состояние вблизи дентальных имплантатов.

1. Постановка задачи и методика численного решения

Расчет напряженно-деформированного состояния выполнен в постановке плоской деформации для модели, состоящей из цилиндрического имплантата с коронкой и окружающих костных тканей при допущении, что костная ткань является изотропным и однородным упругим материалом. Значения модуля упругости кортикальной и губчатой костных тканей изменяются в достаточно широком диапазоне [17]. При выполнении данного исследования полагалось, что модуль упругости губчатой костной ткани изменяется в интервале ${E_s} = 0.5 \div 8$ ГПа, модуль упругости кортикальной костной ткани — ${E_c} = 6 \div 20$ ГПа.

Расчеты выполнялись при фиксированном модуле упругости костной ткани одного типа и вариации модуля упругости костной ткани другого типа: 

1. для каждого из трех значений модуля упругости кортикальной кости ${E_c} = \{6; 10; 18\}$ ГПа выполнялся расчет при различных значениях модуля упругости губчатой кости ${E_s} = \{0.5; 1; 3; 5; 8\}$ ГПа;
2. для каждого из трех значений модуля упругости губчатой кости ${E_s} = \{0.5; 1; 3\}$ ГПа выполнялся расчет при различных значениях модуля упругости кортикальной кости ${E_c} = \{6; 8; 10; 14; 18; 20\}$ ГПа.

Полагалось также, что коронка изготовлена из керамики, а остальные детали имплантата изготовлены из титана. Базовые механические свойства материалов имплантата и костных тканей приведены в табл. 1 по данным работ [10, 11, 17].

Таблица 1. Механические свойства имплантата и костных тканей
[Mechanical properties of the implant and bone tissues]

Material	Young's modulus $E$, GPa	Poisson's ratio	Yield strength (a), Tensile strength (b), MPa
Titanium	110	0.3	880 (a)
Ceramics	70	0.19	320 (b)
Cortical bone	18	0.3	130 (a)
Spongy bone	0.5	0.35	$10\div 20$ (b)

Расчетная модель имплантата и костных тканей содержит семь подобластей (см. рис. 1). Первые три подобласти моделируют костные ткани: 1 — внешний слой кортикальной кости (средняя толщина слоя 1.7 мм); 2 — губчатая кость — основная часть костной ткани челюсти; 3 — слой кортикальной кости толщиной 1 мм, формирующий стенки нижнечелюстного канала, имеющего диаметр 3 мм. Остальные подобласти соответствуют элементам конструкции имплантата (длина цилиндрической части имплантата — 12 мм, диаметр — 4.2 мм, см. подобласть 4). Вся модель полностью с указанием вида нагрузки и граничных условий, а также увеличенное изображение верхней части модели представлены на рис. 2 (для удобства графического представления внешний слой кортикальной кости, цилиндрическая часть имплантата и коронка заштрихованы). Между всеми подобластями модели полагается выполнение условий идеального сцепления, включая зону соединения цилиндрической части имплантата и костных тканей (полная остеоинтеграция). Нагрузка величиной $p=6$ МПа, соответствующая усилию жевания [18], прикладывается к окклюзионной поверхности коронки в двух вариантах (в вертикальном направлении и под углом $45^\circ$ к вертикали). Полагая, что площадь контактной поверхности коронки $S\approx 0.25 \text{|см|}^2$, можно считать, что заданное усилие $p$ соответствует приложенной силе $F\approx 150 \text{|Н|}$. По внешней границе кортикальной кости полагаем нулевые перемещения по направлениям осей координат, тонкая заштрихованная зона на рис. 2 показывает область заданных перемещений по внешней границе модели (это обозначение только на рисунке, не является частью модели).

Интенсивность напряжений (эквивалентные напряжения, напряжения по Мизесу) для состояния плоской деформации определяется следующим образом:
\[ \begin{equation*}
\sigma _i = \sqrt {
\frac{(\sigma _1 - \sigma _2)^2 + (\sigma _2 - \sigma _3 )^2 + (\sigma _3 - \sigma _1)^2 }{2}},
\end{equation*} \]
где $\sigma _1$, $\sigma _2$, $\sigma _3$ — главные напряжения.

Рис. 1. Подобласти модели: 1, 3 — кортикальная кость, 2 — губчатая кость, 4 — имплантат, 5 — винт, 6 — абатмент, 7 — керамическая коронка под действием наклонной нагрузки
[Figure 1. The model subregions: 1, 3 — cortical bone, 2 — spongy bone, 4 — implant, 5 — screw, 6 — abatment, 7 — ceramics crown under inclined load]

Рис. 2. Дискретизация границ подобластей имплантата и костных тканей; действие наклонной нагрузки; 7 подобластей; 1106 узлов во всей модели
[Figure 2. Boundary elements discretization of the implant and surrounded bone tissues under inclined load application; 7 subregions; total 1106 nodes]

Для численного анализа напряженно-деформированного состояния в работе используется прямой вариант МГЭ[19]. Двумерная составная конструкция моделируется набором однородных подобластей. Такой подход позволяет рассматривать соединения различных материалов в конструкции, а также условия неидеального соединения
между подобластями. Граничное интегральное уравнение (ГИУ) теории упругости для каждой однородной подобласти конструкции без учета массовых сил и температурных воздействий имеет вид [19]
\[ \begin{equation}
c_{ij} ( p )u_i ( p ) =
\int _\Gamma \bigl[ 
G_{ij} ( p, q ) t_i ( q ) - F_{ij} ( p, q ) u_i ( q ) 
\bigr] d\Gamma ( q ),\quad i, j = 1, 2,
\end{equation} \tag{1} \]
где точки $p$ и $q$ принадлежат границе подобласти $\Gamma$; $u_i ( q )$, $t_i ( q )$ — перемещения и усилия на границе подобласти; $c_{ij} ( p )$ — функции, зависящие от локальной геометрии в точке $p$; для гладкой границы $c_{ij} ( p )= \delta_{ij}/2$. Функции $ G_{ij} (p, q )$ и $ F_{ij} ( p, q )$ являются фундаментальными решениями Кельвина для задачи теории упругости (плоская деформация) [19].

На участках идеального контакта подобластей полагаем, что выполняются условия непрерывности для перемещений и равновесия для усилий:
\[ \begin{equation}
u_i^k ( p ) = u_i^n ( p ),
\quad 
t_i^k(p) = - t_i^n (p) ,
\end{equation} \tag{2} \]
где $k$ и $n$ — номера соседних подобластей. 

Для численного решения ГИУ (1) используются изопараметрические квадратичные граничные элементы. Дискретизация границ подобластей и аппроксимация переменных позволяет заменить вычисление интегралов по всей поверхности подобласти в ГИУ (1) суммой интегралов по совокупности граничных элементов, представляющих поверхность подобласти, и получить для каждой точки коллокации $p$ дискретное представление ГИУ. Сингулярное поведение фундаментальных решений в ГИУ (1) проявляется в основном при вычислении интегралов по элементам, содержащим точку коллокации $p$. Поэтому рассматривается два случая расположения точки $p$: $p \notin e$ и $p \in e$, где $e$ — граничный элемент, по которому выполняется интегрирование. В первом случае интегрирование выполняется по квадратурным формулам Гаусса – Лежандра с учетом сгущения точек интегрирования в направлении точки коллокации. При $p \in e$ фундаментальное решение плоской задачи теории упругости для перемещений $G_{ij} (p,q )$ содержит слабую логарифмическую особенность и интегралы в дискретном представлении ГИУ, содержащие эту функцию, вычисляются с использованием квадратурной формулы с логарифмической весовой функцией. Фундаментальное решение для поверхностных усилий $F_{ij} (p, q )$ содержит сильную особенность и вычисление суммы соответствующих интегралов и параметра $c_{ij} ( p )$ выполняется косвенным путем исходя из смещения тела как жесткого целого (см. описание алгоритма численного интегрирования для двумерных и пространственных задач в [20]). Обходя последовательно узловые точки границ всех подобластей, получаем дискретное представление ГИУ. Учет граничных условий на внешних границах подобластей и дополнительных условий на участках соединения подобластей (2) позволяет получить из системы дискретных ГИУ для всех подобластей конструкции систему линейных алгебраических уравнений для определения узловых неизвестных задачи. Методика численного решения ГИУ реализована в комплексе программ «МЕГРЭ» (Метод ГРаничных Элементов) [20–22], который использован для получения результатов, приведенных в данной статье.

2. Результаты расчетов при нормальной нагрузке на коронку

Первый этап расчета выполнен при сжимающей нагрузке, приложенной к окклюзионной поверхности коронки. Зависимости максимальных эквивалентных напряжений в костных тканях от модулей упругости кортикальной и губчатой костной ткани представлены на рис. 3–6.

При увеличении модуля упругости губчатой костной ткани $E_s$ (увеличении жесткости кости) возрастает часть нагрузки, передаваемая на эту костную ткань. Поэтому
максимальные напряжения в губчатой костной ткани возрастают, см. рис. 3. Напряжения в кортикальной костной ткани при увеличении модуля упругости губчатой кости $E_s$ снижаются за счет уменьшения нагрузки, передаваемой на эту часть кости, см. рис. 4. Отметим, что напряжения в губчатой кости снижаются при увеличении модуля упругости кортикальной кости $E_c$, см. рис. 3 и 5. С другой стороны, уровень напряжений в кортикальной кости возрастает при увеличении модуля упругости этой кости, см. рис. 4 и 6.

Рис. 3. Зависимость максимальных напряжений в губчатой кости от модуля упругости губчатой кости; сжатие
[Figure 3. The maximum stresses in the spongy bone vs the modulus of elasticity of the spongy bone; compression]

Рис. 4. Зависимость максимальных напряжений в кортикальной кости от модуля упругости губчатой кости; сжатие
[Figure 4. The maximum stresses in the cortical bone vs the modulus of elasticity of the spongy bone; compression]

Рис. 5. Зависимость максимальных напряжений в губчатой кости от модуля упругости кортикальной кости; сжатие
[Figure 5. The maximum stresses in the spongy bone vs the modulus of elasticity of the cortical bone, compression]

Рис. 6. Зависимость максимальных напряжений в кортикальной кости от модуля упругости кортикальной кости; сжатие
[Figure 6. The maximum stresses in the cortical bone vs the modulus of elasticity of the cortical bone; compression]

Рис. 7. Интенсивность напряжений $\sigma_i$ вдоль контура губчатой кости; сжатие, $E_c=18$ ГПа: a) $\sigma_{i,\max}=2.6$ МПа, $E_s=0.5$ ГПа; b) $\sigma_{i,\max}=3.9$ МПа, $E_s=5.0$ ГПа
[Figure 7. Stresses intensity $\sigma_i$ along the spongy bone contour; compression, $E_{c}=18$ GPa: a) $\sigma_{i,\max}=2.6$ MPa, $E_{s}=0.5$ GPa; b) $\sigma_{i,\max}=3.9$ MPa, $E_{s}=5.0$ GPa]

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

Рис. 8. Зависимость максимальных напряжений в губчатой кости от модуля упругости губчатой кости; наклонная нагрузка
[Figure 8. The maximum stresses in the spongy bone vs the modulus of elasticity of the spongy bone; inclined load]

Рис. 9. Зависимость максимальных напряжений в кортикальной кости от модуля упругости губчатой кости; наклонная нагрузка
[Figure 9. The maximum stresses in the cortical bone vs the modulus of elasticity of the spongy bone; inclined load]

Рис. 10. Зависимость максимальных напряжений в губчатой кости от модуля упругости кортикальной кости; наклонная нагрузка
[Figure 10. The maximum stresses in the spongy bone vs the modulus of elasticity of the cortical bone; inclined load]

Рис. 11. Зависимость максимальных напряжений в кортикальной кости от модуля упругости кортикальной кости; наклонная нагрузка
[Figure 11. The maximum stresses in the cortical bone vs the modulus of elasticity of the cortical bone; inclined load]

3. Результаты расчетов при наклонной нагрузке на коронку

Второй этап расчета на модели имплантата выполнен при наклонной нагрузке, приложенной к коронке (см. рис. 2). Зависимости максимальных эквивалентных напряжений в костных тканях от модулей упругости кортикальной и губчатой кости представлены на рис. 8–11. При наклонной нагрузке тенденции изменения максимальных эквивалентных напряжений аналогичны полученным ранее при сжимающей нагрузке, но уровень максимальных напряжений выше, что связано с действием сдвиговой компоненты нагрузки $p_x$, приводящей к возникновению изгибающих моментов $M_1\approx p_xh$ в слое кортикальной кости и $M_2\approx p_x(h+L)$ — у основания имплантата, где $h$ — высота коронки, $L$ — длина имплантата.

Увеличение модуля упругости губчатой костной ткани $E_s$ приводит к возрастанию максимальных напряжений в губчатой кости, причем рост напряжений происходит более интенсивно, чем при нормальной нагрузке, см. рис. 3 и 8. Отметим, что при наклонной нагрузке и малых значениях модуля упругости губчатой кости ($E_s<2$ МПа) изменение напряжений в кортикальной кости слабо зависит от модуля упругости кортикальной костной ткани (см. рис. 4 и 9). Напряжения в губчатой кости снижаются при увеличении модуля упругости кортикальной кости $E_c$, см. рис. 8 и 10. При этом значения максимальных напряжений значительно (примерно в три раза) превышают соответствующие значения при нормальной нагрузке (см. рис. 3, 5 и рис. 8, 10), что связано с действием изгибающего момента. Уровень напряжений в кортикальной кости возрастает при увеличении модуля упругости этой кости, (см. рис. 9 и 11). Для максимальных напряжений в кортикальной кости при наклонной нагрузке отличия от случая нормальной нагрузки невелики (см. рис. 4, 6 и рис. 9, 11), что связано с меньшей величиной механического момента $M_1$ в слое кортикальной кости.

Распределение интенсивности напряжений по границе губчатой кости при действии наклонной нагрузки представлено на рис. 12. При этом характер распределения напряжений в губчатой кости заметно изменяется: при $E_s=0.5$ ГПа максимальные напряжения наблюдаются в нижней части зоны контакта губчатой кости и основания имплантата, а при $E_s=5.0$ ГПа положение максимума напряжений смещается в зону шейки имплантата на границу с губчатой и кортикальной костными тканями.

Рис. 12. Интенсивность напряжений $\sigma_i$ вдоль контура губчатой кости; сжатие, $E_c=18$ ГПа: a) $\sigma_{i,\max}=4.3$ МПа, $E_s=0.5$ ГПа; b) $\sigma_{i,\max}=12.2$ МПа, $E_s=5.0$ ГПа
[Figure 12. Stresses intensity $\sigma_i$ along the spongy bone contour; compression, $E_{c}=18$ GPa: a) $\sigma_{i,\max}=4.3$ MPa, $E_{s}=0.5$ GPa; b) $\sigma_{i,\max}=12.2$ MPa, $E_{s}=5.0$ GPa]

Заключение

Исследование зависимости интенсивности напряжений от плотности кортикальной и губчатой костных тканей выполнено посредством параметрической вариации модулей упругости соответствующих костных тканей для модели соединения цилиндрического имплантата и костных тканей при действии нормальной и наклонной нагрузок. Полученные результаты качественно (ввиду различия геометрических параметров моделей и механических свойств) согласуются с результатами работ [3, 8, 11].

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

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

Конкурирующие интересы. Конкурирующих интересов не имею.
Авторская ответственность. Я несу полную ответственность за предоставление окончательной версии рукописи в печать. Окончательная версия рукописи мною одобрена.
Финансирование. Работа выполнена в рамках государственного задания (№ госрегистрации 123021700050-1).

About the authors

Mikhail N. Perelmuter

Ishlinsky Institute for Problems in Mechanics, Russian Academy of Sciences

Author for correspondence.
Email: perelm@ipmnet.ru
ORCID iD: 0000-0002-8430-5412
SPIN-code: 1057-0990
Scopus Author ID: 8156746000
ResearcherId: J-1283-2014

Dr. Phys. & Math. Sci.; Leading Researcher; Lab. of Mechanics of Strength and Fracture of Materials and Structures

Russian Federation, 119526, Moscow, pr. Vernadskogo, 101–1

References

Supplementary files