CC BY
7
МЕХАНИКА ДЕФОРМИРУЕМОГО ТВЕРДОГО ТЕЛА
DEFORMATION IN SOLID MECHANICS
https://doi.org/10.21122/2227-1031-2023-22-5-376-386 УДК 517.97:615.47
О возможности применения программы COMSOL Multiphysics для топологической оптимизации пластин для остеосинтеза
Докт. техн. наук Д. А. Степаненко1),
асп. И. Мудинов1), студенты В. А. Охремчик1), А. А. Билейчик1)
^Белорусский национальный технический университет (Минск, Республика Беларусь)
© Белорусский национальный технический университет, 2023 Belarusian National Technical University, 2023
Реферат. В статье описана методика топологической оптимизации пластин для остеосинтеза, применяемых для внутренней фиксации переломов костей. Представленная методика основана на применении метода плотности и коммерчески доступного программного обеспечения COMSOL Multiphysics, предназначенного для моделирования с помощью метода конечных элементов. Проведен сравнительный анализ характеристик (продольной жесткости, объема и максимального напряжения, по Мизесу) базовой конструкции пластины и двух оптимизированных вариантов. Установлено, что оптимизированные варианты обеспечивают снижение объема пластины на 49-54 %. При этом продольная жесткость уменьшается на 43-53 %, что является положительным эффектом с точки зрения минимизации эффекта экранирования напряжений. Оптимизированные варианты конструкции имеют близкие значения продольной жесткости и максимального напряжения, по Мизесу, однако в одном из них возникает прогиб продольных сегментов, приводящий к повышению полной энергии деформации, которая используется в качестве целевой функции при оптимизации. В варианте 2 конструкции прогиб продольных сегментов пластины исключается за счет наличия поперечной перемычки между ними и полная энергия деформации принимает более низкое значение. Вариант конструкции без перемычки требует дополнительного исследования, так как сдвиговые напряжения, возникающие в результате контактного взаимодействия продольных сегментов пластины с костью, могут оказывать положительное влияние на регенерацию костной ткани.
Ключевые слова: пластины для остеосинтеза, метод конечных элементов, COMSOL Multiphysics, топологическая оптимизация, продольная жесткость, эффект экранирования напряжений
Для цитирования: О возможности применения программы COMSOL Multiphysics для топологической оптимизации пластин для остеосинтеза Д. А. Степаненко [и др.] // Наука и техника. 2023. Т. 22, № 5. С. 376-386. https://doi. org/10.21122/2227-1031-2023-22-5-376-386
On Possibility of Application of COMSOL Multiphysics Software for Topological Optimization of Osteosynthesis Plates
D. A. Stepanenko1), I. Mudinov1), V. A. Akhremchyk1), H. A. Bileichyk1)
^Belarusian National Technical University (Minsk, Republic of Belarus)
Abstract. The paper describes a technique of topological optimization of osteosynthesis plates used for internal fixation of bone fractures. The proposed technique is based on the application of the density method and the commercially available
Адрес для переписки
Степаненко Дмитрий Александрович Белорусский национальный технический университет ул. Я. Коласа, 22,
220013, г. Минск, Республика Беларусь Тел.: +375 17 293-91-01 [email protected]
Address for correspondence
Stepanenko Dmitry A.
Belarusian National Technical University
22, Ya. Kolasa str.,
220013, Minsk, Republic of Belarus
Tel.: +375 17 293-91-01
Наука
итехника. Т. 22, № 5 (2023)
COMSOL Multiphysics software intended for finite element modeling. A comparative analysis of the characteristics (axial stiffness, volume and maximum von Mises stress) is presented for initial design of the plate and two optimized variants of the design. It has been established that the optimized variants provide a reduction in the plate volume by 49-54 %. In this case, the axial stiffness decreases by 43-53 %, which is a positive effect in terms of minimizing the effect of stress shielding. The optimized variants of the design possess close values of axial stiffness and maximum von Mises stress, however, in one of them, deflection of the axial segments occurs, resulting in an increase in the total strain energy, which is used as an objective function during optimization. In the variant 2 of the design, the deflection of the longitudinal segments of the plate is eliminated due to the presence of a transverse bridge between them, and the total strain energy takes on a lower value. The variant of the design without a bridge should be additionally studied, since shear stresses resulting from the contact interaction of the longitudinal segments of the plate with the bone can have a positive effect on regeneration of the bone tissue.
Keywords: osteosynthesis plates, finite element method, Comsol Multiphysics, topological optimization, axial stiffness, stress shielding effect
For citation: Stepanenko D. A., Mudinov I., Akhremchyk V. A., Bileichyk H. A. (2023) On Possibility of Application of COMSOL Multiphysics Software for Topological Optimization of Osteosynthesis Plates. Science and Technique. 22 (5), 376-386. https://doi.org/10.21122/2227-1031-2023-22-5-376-386 (in Russian)
Общие сведения о методах остеосинтеза
Переломы костей являются одним из самых распространенных видов травм. Так, в 2019 г. общее число случаев переломов в мире составило около 178 млн [1]. Существующие методы лечения переломов можно разделить на консервативные, например наложение гипсовых повязок после закрытой репозиции костных отломков, и хирургические. Разновидность хирургических методов - накостный (экстрамедуллярный) остеосинтез, при котором отломки фиксируются в правильном взаимном положении, достигнутом путем репозиции, с помощью пластин-фиксаторов, соединяемых с костью шурупами или винтами. Простейший вариант пластин для остеосинтеза - пластины с круглыми отверстиями, однако они редко используются в настоящее время, так как не способны обеспечить компрессию отломков и по этой причине должны использоваться в сочетании с внешним компрессионным устройством. В динамических компрессионных пластинах (dynamic compression plate, DCP) используются отверстия особой формы, имеющие наклонный цилиндрический участок, при взаимодействии которого со сферической опорной поверхностью крепежного винта возникает горизонтальное перемещение пластины и связанного с ней отломка кости, приводящее к закрытию перелома и компрессии отломков [2]. Важным фактором, способствующим повышению эффективности экстрамедуллярного остеосинтеза, является сохранение целостности сосудистой системы надкостницы (периоста). Для решения этой проблемы разработаны динамические компрессионные пластины с ограниченным контактом (limited contact DCP, LC-DCP) и точечным контактом (point contact DCP, PC-DCP) [3].
В динамических компрессионных пластинах с ограниченным контактом на нижней поверхности пластины, контактирующей с периостом, создается профиль, способствующий снижению площади контакта, что, в свою очередь, приводит к стимуляции ангиогенеза и повышению плотности кровеносных сосудов [4]. В пластинах с точечным контактом профиль нижней поверхности обеспечивает формирование по обе стороны от крепежных отверстий заостренных выступов, обеспечивающих контакт пластины с периостом, близкий по характеру к точечному. В блокируемых компрессионных пластинах (locking compression plate, LCP) используются комбинированные отверстия, позволяющие применять как стандартные крепежные, так и блокирующие винты, головка которых имеет коническую резьбу, ввинчиваемую в резьбу крепежного отверстия, что обеспечивает стабильность углового положения пластины.
Эффект экранирования напряжений
и способы его минимизации
Применение пластин для хирургического лечения переломов может приводить к разгрузке (экранированию) кости от механических напряжений, действующих на нее у здорового человека; это явление называют эффектом экранирования напряжений (stress shielding) [5]. В результате экранирования напряжений возникает остеопения - снижение минеральной плотности и прочности костной ткани [6]. Данный негативный эффект объясняется законом Вольфа, согласно которому микроархитектура, свойства и общая морфология костной ткани адаптируются к изменениям действующих на нее внешних нагрузок [7]. В основу за-
Н Наука
итехника. Т. 22, № 5 (2023)
кона положено наблюдение немецкого хирурга Юлиуса Вольфа, установившего, что пучки трабекул в бедренной кости направлены вдоль траекторий главных напряжений [8]. Более глубокое понимание механизмов адаптации костной ткани к внешним нагрузкам дает так называемая ютская парадигма (Utah paradigm), согласно которой напряжения, возникающие в кости под действием внешней нагрузки, регистрируются механочувствительными клетками-остеоцитами и сравниваются ими с предписанными пороговыми значениями, нижнее из которых соответствует началу процесса разрушения (резорбции) костной ткани при недостаточном уровне нагрузок, а верхнее - началу процесса формирования дополнительной костной ткани при избыточном уровне нагрузок [9]. Остеоциты преобразуют механические напряжения в биологические сигналы, управляющие активностью остеобластов, участвующих в процессе формирования костной ткани путем преобразования остеобластов в остеоциты, и остеокластов, участвующих в процессе резорбции путем растворения минеральной составляющей и разрушения коллагена [10]. Эффект экранирования напряжений нежелателен на заключительных стадиях процесса срастания перелома, однако играет положительную роль на его начальном этапе, так как высокая жесткость пластины обеспечивает надежную фиксацию отломков кости. Идеален вариант, когда жесткость пластины и соответственно степень экранирования напряжений постепенно снижаются с течением времени по мере срастания перелома. Реализация такого сценария возможна при использовании пластин из биоде-градируемых материалов, например синтетических полимеров и магния [11], но в этом случае существует проблема согласования скорости деградации материала со скоростью срастания перелома. Преимуществом использования би-одеградируемых материалов является отсутствие необходимости в проведении повторного хирургического вмешательства для извлечения пластины. Идея изменения жесткости пластины по мере срастания перелома положена в основу концепции «динамизируемых» пластин [12], однако в отличие от пластин из биодеградиру-емых материалов процесс изменения жесткости динамизируемых пластин управляем и характеризуется скачкообразным снижением жесткости в определенный момент времени, выбор
которого зависит от скорости срастания перелома. Конструктивно динамизируемая пластина состоит из двух частей, возможность взаимного осевого перемещения которых определяется положением подвижного фиксатора: на начальных стадиях срастания перелома фиксатор находится в положении, исключающем взаимное перемещение частей пластины, что обеспечивает высокую жесткость, а на заключительной стадии переводится с помощью привода в положение, допускающее взаимное перемещение, соответствующее низкой жесткости. Использование магнитного привода позволяет управлять положением фиксатора неинвазив-ным образом. При использовании принципа селективного экранирования напряжений пластина проектируется таким образом, чтобы она обладала высокой жесткостью по отношению ко всем видам деформаций (кручение и изгиб), за исключением деформации растяжения-сжатия [13]. Такая пластина будет обеспечивать возможность осевых микроперемещений отломков кости и частичной передачи внешних нагрузок на кость в осевом направлении, исключая при этом взаимное смещение отломков по остальным степеням свободы. Технологически подобная пластина может быть реализована за счет использования анизотропных композиционных материалов, например ламинатов, состоящих из слоев «углеродное волокно/эпоксидная смола» и «льняное волокно/эпоксидная смола» [13].
Топологическая оптимизация
имплантатов и способы ее реализации
В настоящее время для создания оптимальных конструкций ортопедических имплантатов широко используется метод топологической оптимизации [14-18]. Топологическая оптимизация распределяет материал в изделии таким образом, чтобы обеспечить минимизацию некоторой целевой функции, характеризующей эффективность работы изделия, при заданных нагрузках, граничных условиях и ограничениях. В отличие от оптимизации формы, в которой возможные варианты конструкции формируются из исходного путем деформации существующих границ, топологическая оптимизация допускает формирование новых границ, например отверстий в изначально сплошном материале, то есть допускает изменение тополо-
Наука
итехника. Т. 22, № 5 (2023)
гии. Результатом топологической оптимизации обычно является создание облегченной по массе конструкции изделия, не уступающей по своим характеристикам базовой конструкции из сплошного материала. Классический пример использования топологической оптимизации -проектирование фермы Мессершмитта-Бёлкова-Блома (Messerschmitt-Bölkow-Blohm beam, MBB beam), обладающей минимальной изгиб-ной податливостью при минимальной массе. В этом случае топологическая оптимизация на сплошной прямоугольной области, соответствующей сплошной балке, приводит к созданию облегченной конструкции в виде решетчатой фермы. Минимизация податливости и минимизация массы являются взаимно конфликтными целями, так как снижение массы путем частичного удаления материала приводит к увеличению податливости, однако при некотором предписанном значении массы топологически оптимизированная конструкция будет обладать минимальной для данного значения массы податливостью (большей, чем податливость исходной конструкции). С точки зрения минимизации эффекта экранирования напряжений, повышение податливости будет положительным эффектом. Существует несколько подходов к решению задач топологической оптимизации: эволюционная структурная оптимизация, метод плотности и метод гомогенизации.
Метод эволюционной структурной оптимизации основан на постепенном удалении элементов из разбитой на дискретные элементы области оптимизации с использованием определенного критерия удаления последних, простейшим вариантом которого является критерий, основанный на уровне напряжений: элементы, для которых напряжение, по Мизесу, оказывается ниже порогового значения, заданного в процентах от максимального напряжения, считаются недостаточно эффективно используемыми и удаляются из конструкции [19]. Напряжения обычно рассчитываются с помощью метода конечных элементов. Процесс продолжается итеративным образом до тех пор, пока в конструкции не останется элементов, удовлетворяющих критерию удаления.
При использовании метода плотности область оптимизации разбивается на конечные элементы, каждому из которых в конечной оптимизированной топологии может соответство-
Н Наука
итехника. Т. 22, № 5 (2023)
вать нулевое значение безразмерной плотности p = 0 (отсутствие материала) или единичное значение плотности p = 1 (наличие материала). Недостатком такой дискретной формулировки является то, что она приводит к некорректно поставленной задаче. Для решения этой проблемы была предложена модель сплошного изотропного материала со штрафованием промежуточных плотностей (SIMP, Solid Isotropic Material with Penalization for intermediate densities), в которой плотность может принимать непрерывно изменяющиеся значения от pmin > 0 до pmax = 1 [20], а модуль упругости определяется выражением
E = Emm + (E0 - Emm)^,
где p >1 - целочисленный штрафной коэффициент; E0 - модуль упругости сплошного материала; Emin - значение модуля упругости при p = Pmin.
При использовании такой модели подавляющее большинство элементов в конечной топологии будет иметь плотность p = pmin или p = 1.
Использование метода гомогенизации рационально при изготовлении оптимизированного изделия из композиционных материалов, например ламинатов и ячеистых композитов (cellular materials), физические свойства которых (упругие свойства и плотность) могут регулироваться за счет изменения отношения х объемов фаз в композите. Развитие данного метода тесно связано с внедрением в производство аддитивных технологий, из которых наиболее важную роль в изготовлении имплан-татов играют селективная лазерная плавка (selective laser melting, SLM) и электронно-лучевая плавка (electron beam melting, EBM) [21]. Использование ячеистых композитов в постоянных имплантатах дает такие преимущества, как возможность врастания костной ткани и потенциальная возможность васкуляризации им-плантата. При использовании метода гомогенизации геометрия изделия остается неизменной, а определению подлежит оптимальное распределение свойств материала в нем. Для различных значений плотности p, соответствующих изменению величины отношения х объемов фаз в композите, предварительно рассчитываются эквивалентные значения упругих свойств (модулей продольной и сдвиговой упругости и коэффициента Пуассона), которые затем
используются для определения напряженно-деформированного состояния изделия в процессе оптимизации. Эквивалентные значения упругих свойств будут зависеть от реальных упругих свойств входящих в состав композита фаз, отношения x объемов фаз, а в случае ячеистых композитов - от геометрии элементарной ячейки.
Особенности топологической оптимизации, основанной на применении метода плотности. При решении задач топологической оптимизации методом плотности могут возникать численные неустойчивости, примерами которых являются паттерны типа «шахматная доска» (checkerboard patterns) [22] и сеточная зависимость решения. «Шахматная доска» представляет собой периодический двумерный паттерн в виде чередующихся областей высокой pmax и низкой pmin плотности, что соответствует среднему значению плотности р = (pmin + pmax)/2. Сеточная зависимость представляет собой зависимость оптимизированной топологии от параметров сетки конечных элементов, используемой при решении задачи. Для подавления численных неустойчивостей используется фильтрация плотности с помощью фильтров нижних пространственных частот [23]. Радиус фильтра rf (величина, обратно пропорциональная пространственной частоте среза) должен быть не меньше максимального размера элементов сетки. Радиус фильтра является параметром, определяющим минимальный возможный размер конструктивных элементов оптимизированного изделия, поэтому он должен выбираться с учетом разрешающей способности технологического оборудования, которое предполагается использовать для изготовления конструкции. Один из существующих подходов к фильтрации плотности - использование фильтра Гельмгольца, суть которого состоит в сведении задачи фильтрации к решению дифференциального уравнения типа Гельмгольца в частных производных [24]. Математически задача фильтрации описывается интегралом свертки между нефильтрованным полем плотности и импульсной характеристикой фильтра. В то же время, задача решения неоднородных дифференциальных уравнений в частных производных также сводится к интегралу свертки: решение уравнения представляется в виде свертки неоднородного члена с функцией Грина. Учитывая математическую эквивалентность
двух упомянутых задач, фильтрация может быть сведена к решению дифференциального уравнения с функцией Грина, вид которой допускает ее использование в качестве импульсной характеристики фильтра нижних частот. Удобство такого подхода состоит в том, что дифференциальные уравнения в частных производных эффективно решаются с помощью метода конечных элементов, являющегося основным инструментом топологической оптимизации.
Идеальное решение задачи топологической оптимизации должно содержать только черные (с плотностью ртт) и белые (с плотностью ртах) элементы, однако из-за непрерывного изменения плотности реальное решение будет содержать серые элементы с промежуточными значениями плотности, реализация которых при изготовлении изделия невозможна, так как они не могут быть интерпретированы в рамках бинарной классификации «отсутствие материала/наличие материала». Для решения этой проблемы может использоваться метод проецирования плотности, в котором на плотность дополнительно воздействуют проекционным фильтром, математически описываемым с помощью регуляризованной (сглаженной) функции Хевисайда [25], например гиперболическим тангенсным фильтром:
Рр = (Л(Р(рт - Ре)) + Л(Ррр))/(Л(Р(1 - Ре)) +
+ Л(РРр)), (1)
где рр - значение плотности после проецирования; рf - фильтрованное значение плотности; Рр - пороговое значение плотности (обычно принимается равным 0,5); Р - коэффициент фильтра (обычно принимает целочисленные значения от 1 до 8).
При использовании выражения (1) разрыв функции Хевисайда заменяется участком, на котором она плавно возрастает от нуля до единицы со скоростью, зависящей от коэффициента р. При увеличении коэффициента Р (1) приближается по форме к функции Хевисайда, которая ставит единицу в соответствие всем значениям, превышающим Рр, и нуль - всем значениям, меньшим Рр, то есть дает только черные и белые элементы. Таким образом, с увеличением коэффициента Р уменьшается количество серых элементов, однако при повышении Р ухудшается сходимость вычислений. Данное противоречие решается путем
Наука
итехника. Т. 22, № 5 (2023)
продолжения решения по параметру ß (beta-continuation): вначале ищется решение задачи оптимизации, соответствующее малому значению ß (нерезкому проецированию), например: ß =1, при котором вычисления имеют хорошую сходимость, но дают топологию с большим числом серых элементов, а затем полученное решение используется в качестве начального приближения для поиска уточненного решения с более высоким значением ß (более резким проецированием). Уточненное решение будет содержать меньшее число серых элементов, а сходимость процесса будет удовлетворительной благодаря использованию нетривиального начального приближения. При необходимости процесс повторяется для ряда значений ß.
Математически задача топологической оптимизации формулируется следующим образом:
min: c(p) = fTu; (2)
s t.: K(p)u = f; (3)
V = Epv < Vc: (4)
0 < Pmin < P < 1, (5)
то есть состоит в минимизации целевой функции (2), представляющей собой работу c(p) упругих сил, при ограничениях (3)-(5).
Минимизация работы упругих сил автоматически влечет за собой минимальное значение податливости. Работа упругих сил зависит от вектора u узловых перемещений и вектора f узловых сил. Ограничения представляют собой условие статического равновесия (уравнение (3)), где K - глобальная матрица жесткости, ограничение объема (уравнение (4)), где Vc - максимально допустимое значение объема V, vi - объем i-го конечного элемента, и ограничение диапазона допустимых значений плотности (уравнение (5)). В качестве исходных данных обычно указывается максимально допустимое значение относительного объема V/V0, где V0 = Zvi - начальный объем.
Реализация топологической оптимизации с помощью программы COMSOL Multiphysics. Для решения задачи топологической оптимизации пластин для остеосинтеза с помощью программы COMSOL Multiphysics использовалась двумерная геометрическая модель, представленная на рис. 1.
■ Наука
итехника. Т. 22, № 5 (2023)
0,045
0,035
0,025
0,015
0,005 0
-0,005 -0,015 -0,025
Рис. 1. Геометрическая модель базовой конструкции пластины
Fig. 1. Geometric model of the initial design of the plate
Учитывая геометрическую симметрию задачи, рассматривалась четверть пластины с наложением шарнирно-подвижных граничных условий (Roller Constraints) на границы 1-2, 2-3, 4-5 и 6-7. Геометрические параметры модели приведены в табл. 1.
Таблица 1
Геометрические параметры модели пластины Geometric parameters of plate model
Обозначение Описание Значение
параметра параметра параметра, мм
L Длина пластины 180
H Ширина пластины 30
d Толщина пластины 3
R Радиус отверстий 4,5
R1 Радиус прилегающих
к отверстиям кольцевых областей 6,75
L1 Координата центров
отверстий 22
L2 Координата центров
отверстий 44
В качестве материала был задан титановый сплав марки 21S - метастабильный Р-сплав титана с химическим составом Ti-15Mo-3N-3Al-0,2Si (в весовых процентах).
Так как расположение и форма крепежных отверстий в пластине должны оставаться неизменными на протяжении процесса оптимизации, отверстия были окружены прилегающими к ним кольцевыми областями, которые описаны в модели с помощью опции Prescribed Material, накладывающей на данные области ограничение p = 1. Максимальный размер конечных элементов задан равным 0,75 мм (1/20 от поперечного размера H/2 модели), а минимальный размер 0,375 мм. Параметры процесса топологической оптимизации приведены в табл. 2.
Таблица 2
Параметры процесса топологической оптимизации Topological optimization process parameters
Параметр Описание в COMSOL Значение
Радиус фильтра ^ Filter radius Rmin 1,5 мм
Коэффициент проекционного фильтра в Projection slope в {2, 4, 6, 8}
Пороговое значение плотности для проекционного фильтра Рр Projection point 0p 0,5
Штрафной коэффициент р SIMP exponent ^>SIMP {1, 2, 3, 4}
Минимальное значение плотности р^п Minimum penalized volume fraction 0min 10-3 (значение по умолчанию)
Максимально допустимое значение относительного объема У/У0 Upper bound of comp1.dtopo1.theta_avg 0,5
Переменная comp1.dtopo1.theta_avg представляет собой среднее значение плотности после оптимизации и численно совпадает с величиной V/V0.
В качестве целевой функции (2) в COMSOL Multiphysics используется глобальная переменная comp1.soild.Ws_tot, представляющая собой полную энергию деформации.
Для продолжения решения по параметру ß использовалась параметрическая прогонка (Parametric Sweep) с применением опции Reuse Solution from Previous Step (использование решения с предыдущего шага) с параметрами ß е {2, 4, 6, 8}, p е {1, 2, 3, 4}, то есть вначале решалась задача оптимизации без штрафования промежуточных плотностей (p = 1) с нерезким проецированием (ß = 2), а затем полученное решение использовалось при построении уточненных решений с постепенным усилением штрафования и резкости проецирования.
Фильтрация производилась с помощью фильтра Гельмгольца, а проецирование - с помощью гиперболического тангенсного фильтра.
При моделировании к наружным (по отношению к поперечной оси y симметрии пластины) частям контуров крепежных отверстий прикладывалась нагрузка Fx = -35 Н. В силу использования линейной модели упругости величина приложенной нагрузки не влияет на результаты топологической оптимизации.
В COMSOL Multiphysics используются четыре переменные плотности, которые могут использоваться для отображения результатов оптимизации:
1) 9c (comp 1 .dtopo1.theta_c, control material volume factor) - исходная плотность p (до фильтрации, проецирования и штрафования);
2) 9f (comp1.dtopo1.theta_f, filtered material volume factor) - фильтрованная плотность pf;
3) 9 (comp1.dtopo1.theta, output material volume factor) - плотность рр после проецирования;
4) 9p (comp1.dtopo1.theta_p, penalized material volume factor) - плотность после штрафования.
Результаты и их обсуждение
В результате оптимизации получено распределение плотности 9p, представленное на рис. 2.
Рис. 2. Расчетное распределение плотности 9p Fig. 2. Calculated distribution of density 9p
Как видно из рисунка, в оптимизированной топологии присутствуют серые элементы, что требует применения бинарной классификации для получения черно-белого дизайна. Классификация выполнялась с помощью программы Mathcad, для чего расчетные значения плотности 9p на регулярной сетке с шагом 0,1 мм экспортировались в текстовый файл, а затем импортировались в Mathcad в виде матрицы. Элементы матрицы трансформировались в целочисленный формат со значениями от 0 (сплошной материал) до 255 (отсутствие материала), что позволяло визуализировать матрицу в виде полутонового BMP-изображения. При этом учитывалось, что в точках сетки, лежащих в областях крепежных отверстий, плотность принимает значение NaN (нечисловое значение), так как данные области не входят в область оптимизации. В трансформированной матрице значения NaN заменялись значениями 255. Рис. 3 иллюстрирует влияние порога классификации на ее результаты.
Наука
итехника. Т. 22, № 5 (2023)
Рис. 3. Влияние порога классификации на ее результаты Fig. 3. Effect of threshold on the results of classification
Серая область соответствует оптимизированной топологии пластины при пороге классификации 50 (вариант 1), а черная область соответствует дополнительному объему, возникающему при повышении порога классификации до 200 (вариант 2). Как видно из рисунка, вариант 2 отличается от варианта 1 наличием тонкой перемычки, связывающей между собой верхний и нижний продольные сегменты пластины. Форма наружного контура оптимизированной пластины и формы возникающих при оптимизации дополнительных отверстий описывались таблицами координат принадлежащих им точек, которые извлекались из матрицы с помощью программы Mathcad. На основе полученных таблиц геометрия оптимизированных пластин представлялась в COMSOL Multiphysics интерполяционными кривыми (рис. 4).
Альтернативным вариантом построения геометрических моделей оптимизированных пла-
0,040 0,030 0,020 0,010 0
-0,010 -0,020
-0,08
-0,06
-0,04
-0,02
стин является преобразование классифицированного ВМР-изображения в векторный формат БХБ с помощью коммерчески доступных конвертеров, например программы Scan2CAD.
В целях сравнительного анализа характеристик базовой конструкции пластины и двух вариантов, полученных в результате топологической оптимизации, проводился статический анализ всех трех вариантов конструкции. Граничные условия, нагрузки, свойства материала и параметры сетки сохранялись такими же, как и при топологической оптимизации.
Результаты сравнительного анализа характеристик различных вариантов конструкции пластины представлены в табл. 3.
Как видно из таблицы, объем пластины уменьшается в результате оптимизации на 49-54 %, а продольная жесткость - на 43-53 %, что является положительным эффектом с точки зрения снижения степени экранирования напряжений. Отрицательный фактор, с точки зрения прочности пластины, - повышение максимального напряжения, по Мизесу, составляющее 19-27 % по отношению к базовой конструкции. Как видно, вариант 1 не является оптимальным в терминах полной энергии деформации, однако он представляет интерес с точки зрения характера возникающих деформаций (рис. 5).
- 0,040 m о
- 0,030 -
- 0,020 -
- 0,010 -
- 0 -
- -0,010 -
- -0.020 -
m
-0,08
-0,06
-0,04
-0,02
Рис. 4. Геометрические модели оптимизированных пластин Fig. 4. Geometric models of optimized plates
Характеристики различных вариантов конструкции пластины Characteristics of various plate designs
Таблица 3
0
0
Вариант конструкции Характеристика
Продольная жесткость, кН/мм Объем, см3 Максимальное напряжение, по Мизесу, МПа Полная энергия деформации, мкДж
Базовый 71,4 15,4 7,34 132,4
Вариант 1 (без перемычки) 33,8 7,1 9,32 240,4
Вариант 2 (с перемычкой) 41,0 7,9 8,72 198,8
H Наука
итехника. Т. 22, № 5 (2023)
т
0,050 0,040 0,030 0,020 0,010 0
-0,010 -0,020 -0.030 -0.040
-0,08 -0,06 -0,04 -0,02
т
0,050 0,040 0,030 0,020 0,010 0
-0,010 -0,020 -0.030 -0.040
-0,08 -0,06 -0,04
-0,02
Рис. 5. Деформированная форма оптимизированных пластин Fig. 5. Deformed shape of optimized plates
Как видно из рисунка, в варианте 1 конструкции возникает поперечная деформация (прогиб) продольных сегментов пластины, величина которой сопоставима с величиной продольной деформации. Это способствует повышению полной энергии деформации по сравнению с вариантом 2, в котором прогиб исключается за счет наличия перемычки. С точки зрения продольной жесткости и максимального напряжения, по Мизесу, варианты 1 и 2 конструкции сопоставимы: вариант 1 обладает несколько меньшей жесткостью, но демонстрирует более высокий уровень напряжений. В связи с этим представляет интерес исследование влияния прогибов пластины, возникающих в варианте 1 конструкции, на процесс срастания перелома: с учетом закона Вольфа сдвиговые напряжения, возникающие в результате контактного взаимодействия продольных сегментов пластины при их прогибе с костью, могут положительно влиять на регенерацию костной ткани, однако эта гипотеза требует экспериментальной проверки.
Качественно оптимизированная топология пластин согласуется с результатами, полученными другими исследователями [18, 26]. В частности, авторами работы [18] с помощью метода плотности спроектированы оптимизированные конструкции, характеризующиеся сужением пластины по мере удаления от поперечной плоскости симметрии и наличием дополнительных отверстий, разделяющих между собой крепежные отверстия. Аналогичный результат получен авторами работы [26], исследовавшими фасонки (gusset plates), используемые
при сборке металлоконструкций, с помощью метода эволюционной структурной оптимизации.
ВЫВОДЫ
1. Разработана методика топологической оптимизации пластин для остеосинтеза с применением программы С0М80Ь МиШрЬу8Ю8.
2. Проведен сравнительный анализ характеристик (продольной жесткости, объема и максимального напряжения, по Мизесу) базовой конструкции пластины и двух оптимизированных вариантов, полученных при различных значениях порога классификации расчетных значений плотности. Установлено, что оптимизированные варианты обеспечивают уменьшение объема пластины на 49-54 %, а продольной жесткости - на 43-53 %, что является положительным эффектом с точки зрения снижения степени экранирования напряжений.
3. Показано, что оптимизированные варианты конструкции пластины сопоставимы с точки зрения продольной жесткости и максимального напряжения, по Мизесу, однако в одном из вариантов возникает прогиб продольных сегментов, величина которого сопоставима со значением продольной деформации, что приводит к повышению полной энергии деформации, являющейся целевой функцией при оптимизации. Для варианта 2 конструкции полная энергия деформации принимает более низкое значение, так как за счет наличия перемычки между продольными сегментами пластины исключается их прогиб.
4. Вариант конструкции с более высокой полной энергией деформации требует дополнитель-
Наука
итехника. Т. 22, № 5 (2023)
m
m
ного исследования, так как сдвиговые напряжения, возникающие в результате контактного взаимодействия продольных сегментов пластины при их прогибе с костью, могут положительно влиять на регенерацию костной ткани.
ЛИТЕРАТУРА
1. Global, Regional and National Burden of Bone Fractures in 204 Countries and Territories, 1990-2019: a Systematic Analysis Form the Global Burden of Disease Study 2019 / GBD 2019 Mental Disorders Collaborators // The Lancet Healthy Longevity. 2022. Vol. 9, issue 2. P. 137-150. https://doi.org/10.1016/s2215-0366(21)00395-3.
2. Allgöver, M. A New Plate for Internal Fixation - the Dynamic Compression Plate (DCP) / M. Allgöver, S. Perren, P. Matter // Injury. 1970. Vol. 2, No 1. P. 40-47. https://doi.org/10.1016/s0020-1383(70)80111-5.
3. The Limited Contact Dynamic Compression Plate (LC-DCP) / S. M. Perren [et al.] // Archives of Orthopaedic and Trauma Surgery. 1990. Vol. 109, N 6. P. 304-310. https://doi.org/10.1007/bf00636166.
4. Antabak, A. Reducing Damage to the Periosteal Capillary Network Caused by Internal Fixation Plating: an Experimental Study / A. Antabak [et al.] // Injury. 2015. Vol. 46, N 6. P. S18-S20. https://doi.org/10.1016/). injury. 2015.10.037.
5. Dai, K. Rational Utilization of the Stress Shielding Effect of Implants / K. Dai // Biomechanics and Biomaterials in Orthopedics. London: Springer-Verlag London, 2004. P. 208-215. https://doi.org/10.1007/978-1-4471-3774-0_22.
6. Gilbert, J. A. Stress Protection Osteopenia Due to Rigid Plating / J. A. Gilbert // Clinical Biomechanics. 1988. Vol. 3, No 3. P. 179-186. https://doi.org/10.1016/0268-0033(88)90065-4.
7. Ahn, A. C. Relevance of Collagen Piezoelectricity to "Wolff's Law": a Critical Review / A. C. Ahn, A. J. Grodzinsky // Medical Engineering & Physics. 2009. Vol. 31, N 7. P. 733-741.
8. Boyle, C. Three-Dimensional Micro-Level Computational Study of Wolff's Law Via Trabecular Bone Remodeling in the Human Proximal Femur Using Design Space Topology Optimization / C. Boyle, I. Y. Kim // Journal of Biomechanics. 2011. Vol. 44, N 5. P. 935-942. https://doi. org/10.1016/j.jbiomech.2010.11.029.
9. Frost, H. M. From Wolff's law to the Utah Paradigm: Insights About Bone Physiology And its Clinical Applications / H. M. Frost // The Anatomical Record. 2001. Vol. 262, iss. 4. P. 398-419. https://doi.org/10.1002/ ar.1049.
10. Palumbo, C. The Osteocyte: From "Prisoner" to "Orche-strator" / C. Palumbo, M. Ferretti // Journal of Functional Morphology and Kinesiology. 2021. Vol. 6, N 1. https://doi.org/10.3390/jfmk6010028.
11. Sheikh, Z. Biodegradable Materials for Bone Repair and Tissue Engineering Applications / Z. Sheikh [et al.] // Materials. 2015. Vol. 8, N 9. P. 5744-5794. https://doi.org/ 10.3390/ma8095273.
12. Dichio, G. Engineering and Manufacturing of a Dyna-mizable Fracture Fixation Device System / G. Dichio [et al.] // Applied Sciences. 2020. Vol. 10, № 19. Article 6844. https://doi.org/10.3390/app10196844.
13. Samiezadeh, S. On Optimization of a Composite Bone Plate Using the Selective Stress Shielding Approach / S. Samiezadeh [et al.] // Journal of the Mechanical Behavior of Biomedical Materials. 2015. Vol. 42. P. 138-153. https://doi.org/10.1016/) .jmbbm.2014.11.015.
14. Kaymaz, I. A New Design for the Humerus Fixation Plate Using a Novel Reliability-Based Topology Optimization Approach to Mitigate the Stress Shielding Effect / I. Kaymaz [et al.] // Clinical Biomechanics. 2022. Vol. 99. Article 105768. https://doi.org/10.1016/j.clinbiomech.2022.105768.
15. Wang, M. Optimal Design and Biomechanical Analysis of a Biomimetic Lightweight Design Plate for distal Tibial Fractures: a Finite Element Analysis / M. Wang [et al.] // Frontiers in Bioengineering and Biotechnology. 2022. Vol. 10. Article 820921. https://doi.org/10.3389/fbioe.2022.820921.
16. Wu, N. The Advances of Topology Optimization Techniques in Orthopedic Implants: a Review / N. Wu [et al.] // Medical & Biological Engineering & Computing. 2021. Vol. 59, N 9. P. 1673-1689. https://doi.org/10.1007/ s11517-021-02361-7
17. Al-Tamimi, A. A. Novel bone fixation implants minimising stress shielding: PhD thesis / A. A. Al-Tamimi. University of Manchester, 2019. 253 p.
18. Gogarty, E. Hierarchical Topology Optimization for Bone Tissue Scaffold: Preliminary Results on the Design of a Fracture Fixation Plate / E. Gogarty, D. Pasini // Engineering and Applied Sciences Optimization. - Heidelberg: Springer, 2015. P. 311-340. https://doi.org/10.1007/978-3-319-18320-6_17.
19. Xie, Y. M. Evolutionary Structural Optimization / Y. M. Xie, G.P. Steven. London: Springer-Verlag London, 1997. 188 p. https://doi.org/10.1007/978-1-4471-0985-3.
20. Bendsee, M. P. Material Interpolation Schemes in Topology Optimization / M. P. Bendsee, O. Sigmund // Archive of Applied Mechanics. 1999. Vol. 69, № 9-10. P. 635-654. https://doi.org/10.1007/s004190050248.
21. Zadpoor, A. A. Additively Manufactured Porous Metallic Biomaterials / A. A. Zadpoor // Journal of Materials Chemistry B. 2019. Vol. 7, No 26. P. 4081-4226. https://doi. org/10.1039/c9tb00420c.
22. Diaz, A. Checkerboard Patterns in Layout Optimization / A. Diaz, O. Sigmund // Structural Optimization. 1995. Vol. 10, N 1. P. 40-45. https://doi.org/10.1007/bf01 743693.
23. Lazarov, B. S. Maximum Length Scale in Density Based Topology Optimization / B. S. Lazarov, F. Wang // Computer Methods in Applied Mechanics and Engineering. 2017. Vol. 318. P. 826-844. https://doi.org/10.1016/j.cma. 2017.02.018.
24. Lazarov, B. S. Filters in Topology Optimization Based on Helmholtz-Type Differential Equations / B. S. Lazarov, O. Sigmund // International Journal for Numerical Methods in Engineering. 2011. Vol. 86, № 6. P. 765-781. https://doi.org/10.1002/nme.3072.
25. Guest, J. K. Elimination Beta-Continuation from Heaviside Projection and Density Filter Algorithms / J. K. Guest, A. Asadpoure, S.-H. Ha // Structural and Multidisciplinary Optimization. 2011. Vol. 44, № 4. P. 443-453. https://doi. org/10.1007/s00158-011-0676-1.
26. Khalaf, A. A. Evolutionary Structural Optimization of Steel Gusset Plates / A. A. Khalaf, M. P. Saka // Journal of Constructional Steel Research. 2007. Vol. 63, № 1. P. 71-81. https://doi.org/10.1016/jjcsr.2006.03.002.
Поступила 22.05.2023 Подписана в печать 25.07.2023 Опубликована онлайн 29.09.2023
Н Наука
итехника. Т. 22, № 5 (2023)
REFERENCES
1. GBD 2019 Mental Disorders Collaborators (2022) Global, Regional and National Burden of Bone Fractures in 204 Countries and Territories, 1990-2019: a Systematic Analysis form the Global Burden of Disease Study 2019. The Lancet Healthy Longevity, 9 (1), 137-150. ttps://doi.org/ 10.1016/s2215-0366(21)00395-3
2. Allgöver M., Perren S., Matter P. (1970) A New Plate for Internal Fixation - the Dynamic Compression Plate (DCP). Injury, 2 (1), 40-47. https://doi.org/10.1016/ s0020-1383(70)80111-5.
3. Perren S. M., Mane K., Pohler O., Predieri M., Steinemann S., Gautier E. (1990) The Limited Contact Dynamic Compression Plate (LC-DCP). Archives of Orthopaedic and Trauma Surgery, 109 (6), 304-310. https://doi.org/10. 1007/bf00636166.
4. Antabak A., Papes D., Haluzan D., Seiwerth S., Fuchs N., Romic I., Davila S., Luetic T. (2015) Reducing Damage to the Periosteal Capillary Network Caused by Internal Fixation Plating: an Experimental Study. Injury, 46 (6), S18-S20. https://doi.org/10.1016/j.injury.2015.10.037.
5. Dai K. (2004) Rational Utilization of the Stress Shielding Effect of Implants. Biomechanics and Biomaterials in Orthopedics. London, Springer-Verlag London, 208-215. https://doi.org/10.1007/978-1-4471-3774-0_22.
6. Gilbert J. A. (1988) Stress Protection Osteopenia Due to Rigid Plating. Clinical Biomechanics, 3 (3), 179-186. https://doi.org/10.1016/0268-0033(88)90065-4.
7. Ahn A. C., Grodzinsky A. J. (2009) Relevance of Collagen Piezoelectricity to "Wolff's Law": a Critical Review. Medical Engineering & Physics, 31 (7), 733-741. https://doi. org/10.1016/j.medengphy.2009.02.006.
8. Boyle C., Kim I. Y. (2011) Three-Dimensional MicroLevel Computational Study of Wolff's Law Via Trabecu-lar Bone Remodeling in the Human Proximal Femur Using Design Space Topology Optimization. Journal of Biomechanics, 44 (5), 935-942. https://doi.org/10.1016/j. jbiomech.2010.11.029.
9. Frost H. M. (2001) From Wolff's Law to the Utah Paradigm: Insights About Bone Physiology and its Clinical Applications. The Anatomical Record, 262 (4), 398-419. https://doi.org/10.1002/ar.1049.
10. Palumbo C., Ferretti M. (2021) The Osteocyte: From "Prisoner" to "Orchestrator". Journal of Functional Morphology and Kinesiology, 6 (1). https://doi.org/10.3390/ jfmk6010028.
11. Sheikh Z., Najeeb S., Khurshid Z., Verma V., Rashid H., Glogauer M. (2015) Biodegradable Materials for Bone Repair and Tissue Engineering Applications. Materials, 8 (9), 5744-5794. https://doi.org/10.3390/ma8095273.
12. Dichio G., Call M., Terzini M., Putame G., Zanetti E. M., Costa P., Audenino A. L. (2020) Engineering and Manufacturing of a Dynamizable Fracture Fixation Device System. Applied Sciences, 10 (19), 6844. https://doi.org/10. 3390/app10196844.
13. Samiezadeh S., Avval P. T., Fawaz Z., Bougherara H. (2015) On Optimization of a Composite Bone Plate Using the Selective Stress Shielding Approach. Journal of the Mechanical Behavior of Biomedical Materials, 42, 138-153. https://doi.org/10.1016/jjmbbm.2014.11.015.
14. Kaymaz I., Murat F., Korkmaz I. H., Yavuz O. (2022) A New Design for the Humerus Fixation Plate Using a Novel Reliability-Based Topology Optimization Approach to Mitigate the Stress Shielding Effect. Clinical Biomechanics, 99, 105768. https://doi.org/10.1016/jxlin biomech.2022.105768.
15. Wang M., Deng Y., Xie P., Tan J., Yang Y., Ouyang H., Zhao D., Huang G., Huang W. (2022) Optimal Design and Biomechanical Analysis of a Biomimetic Lightweight Design Plate for Distal Tibial Fractures: a Finite Element Analysis. Frontiers in Bioengineering and Biotechnology, 10, 820921. https://doi.org/10.3389/fbioe.2022.820921.
16. Wu N., Li S., Zhang B., Wang C., Chen B., Han Q., Wang J. (2021) The Advances of Topology Optimization Techniques in Orthopedic Implants: a Review. Medical & Biological Engineering & Computing, 59 (9), 1673-1689. https://doi.org/10.1007/s11517-021-02361-7.
17. Al-Tamimi A. A. (2019) Novel Bone Fixation Implants Minimising Stress Shielding. PhD Thesis. University of Manchester. 253.
18. Gogarty E., Pasini D. (2015) Hierarchical Topology Optimization for Bone Tissue Scaffold: Preliminary Results on the Design of a Fracture Fixation Plate. Lagaros N., Pa-padrakakis M. (eds). Engineering and Applied Sciences Optimization. Heidelberg, Springer, 311-340. https://doi. org/10.1007/978-3-319-18320-6_17.
19. Xie Y. M., Steven G. P. (1997) Evolutionary Structural Optimization. London, Springer-Verlag London. 188. https://doi.org/10.1007/978-1-4471-0985-3.
20. Bendsee M. P., Sigmund O. (1999) Material Interpolation Schemes in Topology Optimization. Archive of Applied Mechanics, 69 (9-10), 635-654. https://doi.org/10.1007/ s004190050248.
21. Zadpoor A. A. (2019) Additively Manufactured Porous Metallic Biomaterials. Journal of Materials Chemistry B, 7 (26), 4081-4226. https://doi.org/10.1039/c9tb00420c.
22. Diaz A., Sigmund O. (1995) Checkerboard Patterns in Layout Optimization. Structural Optimization, 10 (1), 40-45. https://doi.org/10.1007/bf01743693.
23. Lazarov B. S., Wang F. (2017) Maximum Length Scale in Density Based Topology Optimization. Computer Methods in Applied Mechanics and Engineering, 318, 826-844. https://doi.org/10.1016/jxma.2017.02.018.
24. Lazarov B. S., Sigmund O. (2011) Filters in Topology Optimization Based on Helmholtz-Type Differential Equations. International Journal for Numerical Methods in Engineering, 86 (6), 765-781. https://doi.org/10.1002/nme.3072.
25. Guest J. K., Asadpoure A., Ha S.-H. (2011) Elimination Beta-Continuation from Heaviside Projection and Density Filter algorithms. Structural and Multidisciplinary Optimization, 44 (4), 443-453. https://doi.org/10.1007/s00158-011-0676-1.
26. Khalaf A. A., Saka M. P. (2007) Evolutionary Structural Optimization of Steel Gusset Plates. Journal of Constructional Steel Research, 63 (1), 71-81. https://doi.org/10. 1016/j.jcsr.2006.03.002.
Received: 22.05.2023 Accepted: 25.07.2023 Published online: 29.09.2023
Наука
итехника. Т. 22, № 5 (2023)
