CC BY
32
УДК 616.6
НЕКОТОРЫЕ АСПЕКТЫ КОМПЬЮТЕРНОГО МОДЕЛИРОВАНИЯ ПРОКСИМАЛЬНОГО ОТДЕЛА БЕДРЕННОЙ КОСТИ
© И.А. Норкин, О.Н. Ямщиков, Д.А. Марков,
Р.К. Абдулнасыров, Д.Н. Перегородов, А.Ю. Заигралов
Ключевые слова: компьютерное моделирование; бедренная кость; персональный компьютер; ЭЭ-модель; поверхность; механика.
В статье дается подробная информация о выполнении компьютерного моделирования при лечении повреждений и заболеваний бедренной кости.
1. НАЗНАЧЕНИЕ ПРОГРАММЫ
Программно-информационный комплекс (ПИК) предназначен для диагностики, подбора оптимальной оперативной тактики, проведения виртуальных хирургических и травматологических операций и предсказания возможных последствий.
Макет ПИК создан на основе математических методов компьютерного моделирования с использованием технологий параллельных вычислений на многопроцессорных вычислительных системах для прогнозирования последствий хирургического вмешательства при лечении переломов бедренной кости.
ПИК может быть использован в исследовательских центрах по разработке средств и методик диагностики для нужд клинической практики, а также в исследовательских центрах, специализирующихся на разработке средств защиты от травм разных видов.
1.1. ТЕХНИЧЕСКИЕ ХАРАКТЕРИСТИКИ
Макет ПИК построен на алгоритмах, использующих:
- биомеханические модели и численные методы, позволяющие проводить адекватное математическое моделирование травматологических и операционных процессов с учетом их специфики;
- сеточно-характеристические методы, позволяющие корректно учесть физику задачи и проводить адекватное моделирование динамических процессов;
- методы моделирования динамических биомеханических процессов в трехмерном случае.
Макет ПИК обеспечивает возможность интерактивной визуализации расчетных данных, позволяющую изучать протекающие процессы.
Макет ПИК состоит из следующих блоков и модулей:
- модуль импорта геометрических примитивов;
- модуль задания механических свойств модели, параметров внешнего воздействия;
- модуль подготовки модели для расчета на многопроцессорной системе;
- параллельный расчетный модуль;
- модуль управления расчетом, в т. ч. балансировки нагрузки и мониторинга системы;
- модуль визуализации результатов расчета и проведения виртуальных хирургических операций.
Макет ПИК реализован с использованием языков программирования Си, Си++ и средств параллельного программирования MPICH2 и функционирует под операционной системой Linux (Unix).
Интерфейс пользователя ПИК, в т. ч. средства визуализации результатов расчета и средства для проведения виртуальных хирургических операций, функционирует под управлением операционных систем Windows и Linux на клиентской персональной ЭВМ.
Макет ПИК функционирует на многопроцессорных вычислительных системах с распределенной памятью и общей памятью. Масштабируем его по числу процессоров до 8.
1.2. ОГРАНИЧЕНИЯ, НАКЛАДЫВАЕМЫЕ НА ОБЛАСТЬ ПРИМЕНЕНИЯ
Область применения макета в настоящий момент ограничена операциями по лечению перелома бедра. В качестве возможных вариантов фиксации рассматриваются самокомпрессирующая широкая пластина, аппарат Илизарова, универсальный канюлированный стержень [1-5].
На основе анатомической базы бедренных костей с учетом биомеханики повреждения по данным томографических исследований производится трехмерная визуализация поврежденной кости и выбирается оптимальный способ хирургической фиксации перелома с учетом особенностей кровоснабжения сегмента, прочностных характеристик кости, фиксаторов и биомеханики функционирования мышц [6-7].
2. УСЛОВИЯ ПРИМЕНЕНИЯ 2.1. ТРЕБОВАНИЯ К ТЕХНИЧЕСКИМ СРЕДСТВАМ
В состав системы входят:
1) сервер высокопроизводительных вычислений -HP ProLiant DL160 G5 с характеристиками не ниже:
- два четырехядерных процессора Intel Xeon E5405 2,0 ГГц;
- оперативная память 16 Гб Fully Buffered DIMM PC2-5300;
- четыре жестких диска SATA емкостью 1 Тб каждый;
- сетевая карта на два порта 10/100/1000 Ethernet;
2) сервер баз данных с характеристиками не ниже:
- процессор AMD Phenom II Х4 810;
- оперативная память DDR3 2048 Mb;
- видеосистема SVGA 896 Mb nVidia 260GT;
- жесткий диск SATA емкостью 500 Gb;
3) клиентские рабочие места с характеристиками не ниже:
- процессор AMD Phenom II Х4 810;
- оперативная память DDR3 2048 Mb;
- видеосистема SVGA 896 Mb nVidia 260GT;
- жесткий диск SATA емкостью 500 Gb.
2.2 .ТРЕБОВАНИЯ К ПРОГРАММНЫМ СРЕДСТВАМ
Должно быть установлено следующее программное обеспечение:
1) сервер высокопроизводительных вычислений:
- Linux Fedora 12 (операционная система с открытым исходным кодом, собранная и скомпилированная Исполнителем для восьмиядерного высокопроизводительного сервера);
- Elmer 6.0 (система конечноэлементных расчетов с открытым исходным кодом, собранная и скомпилированная Исполнителем для моделирования на восьмиядерном параллельном вычислительном комплексе травматологических и операционных процессов);
2) сервер базы данных:
- Microsoft Windows XP Professional (операционная система);
- MySQL 5.0.45 (система управления базами данных);
- PHP 5.2.4 (интерпретатор языка разработки вебприложений);
- Apache 2.2.4 (веб-сервер);
3) клиентское рабочее место:
- Microsoft Windows XP Professional (операционная система);
- Salome Version 5.1.3 (программа построения твердотельных моделей);
- Microsoft Internet Explorer 8.0 (браузер);
- PuTTY 0.6 (средство терминального доступа);
- WinSCP 4.2.5 (SFTP и FTP Windows-клиент с поддержкой протокола SCP);
- Xming 6.9.0.31 (клиент X Window System для операционной системы Microsoft Windows).
2.3. ОБЩАЯ ХАРАКТЕРИСТИКА ВХОДНОЙ И ВЫХОДНОЙ ИНФОРМАЦИИ
Входной информацией для ПИК служат файлы твердотельных моделей костей, металлофиксаторов и переломов в формате IGES, информация о механических свойствах костей и металлофиксаторов, а также информация об объемных и поверхностных силах, воздействующих на построенную модель [8].
Выходной информацией служат результаты визуализации, которые представляют собой трехмерную
модель, окрашенную в различные цвета в зависимости от возникших в модели смещений и напряжений как по осям, так и абсолютные или эквивалентные напряжения Фон Мизеса, соответственно.
Формат входных и выходных данных (файлов) создаваемой оптимизированной редакции решателя остался неизменным и совпадающим с форматом входных и выходных данных (файлов) исходной версии Elmer 5.5.5. Это обстоятельство позволяет использовать в качестве основы модулей препроцессинга и постпроцессинга свободно распространяемые компоненты пакета Elmer, в частности Elmer GUI для 32-разрядных ОС Windows. Кроме того, модули импорта-экспорта моделей могут работать с имеющейся версией пакета SALOME/Netgen для 32-разрядных ОС Windows. Заметим, что в 32/64-битных версиях OC Linux x86/x86_64 указанные модули ПИК могут быть созданы на основе свободно распространяемого исходного кода при помощи комплекта GNU (C/C++ - gcc/g++).
3. ОПИСАНИЕ ЗАДАЧИ
Программно-информационный комплекс (ПИК) предназначен для моделирования напряженно-деформированного состояния в системе «костная ткань-ме-таллофиксатор» в ходе проведения виртуальной операции остеосинтеза.
ПИК содержит три базовых компоненты:
1) препроцессор, в функции которого входит подготовка либо импорт ранее подготовленной (средствами CAD) трехмерной геометрической модели и генерация конечно-элементной сетки. Препроцессор генерирует данные, характеризующие граничные условия, физические свойства среды в различных пространственных областях и т. д. с помощью графического интерфейса пользователя;
2) решатель (Solver), получающий от препроцессора сгенерированную конечно-элементную сетку и выполняющий назначение свойств среды конкретным группам конечных элементов. В решателе выполняется дискретизация уравнений математической физики при помощи вариационного либо проекционного варианта метода Галеркина, после чего выполняется их численное решение;
3) постпроцессор, который принимает данные, полученные в результате работы решателя, и далее выполняет визуализацию расчетных данных (например, полей упругих смещений, деформаций и напряжений в контексте механики деформируемого твердого тела).
Между двумя «смежными» компонентами требуется лишь совместимость на уровне формата входных и выходных файлов.
При разработке ПИК использовалось свободно распространяемое программное обеспечение. CAE-система скомпилирована на основе созданного решателя (solvers). Вся структура программных модулей ПИК как CAE-системы и системное программное обеспечение адаптированы к решаемой задаче трехмерного моделирования переломов бедренной кости взрослого человека.
ПИК представляет собой распределенный программный комплекс, функционирующий в гетерогенной программно-аппаратной среде.
Препроцессор ПИК представляет собой стандартную для выполнения конечно-элементного моделиро-
вания свободно распространяемую систему трехмерного геометрического моделирования SALOME с поддержкой генерации трехмерных конечно-элементных сеток на основе свободно распространяемых библиотек Netgen, дополненную специально разработанными конвертерами STL-IGES.
Выполняющий генерацию конечно-элементных сеток пакет SALOME/Netgen получает в качестве входных данных трехмерные геометрические модели в формате IGES - «твердотельную модель». Генерируемый современными медицинскими томографами формат файлов DICOM несовместим с форматом входных файлов конечно-элементных препроцессоров. В связи с этим в ПИК реализован оригинальный эффективный алгоритм генерации трехмерных твердотельных объектов IGES по геометрическим моделям STL, учитывающий особенности структур данных, генерируемых медицинскими томографами.
Исходными данными для обратного инжиниринга в задаче получения SD-модели бедренной кости человека являлось медицинское изображение в стандарте DICOM, полученное с помощью КТ-сканера (томографа) TOCHIBA Aquilion. DICOM формат (.dcm файлы) представляет собой медицинские изображения, где в каждом файле содержится изображение среза сканируемого объекта (бедренной кости человека), сделанного с определенным максимально допустимым по уровню облучения для человека шагом [9-12].
В ПИК реализована следующая технология получения SD-модели бедренной кости человека [13-14].
1) используя программу 3D Slicer Version 3.4, создать на основе .dcm файлов формата DICOM, SD-модель формата mesh (создание SD-модели бедренной кости человека на основе сплайновой аппроксимации и треугольных граней);
2) преобразовать полученную SD-модель формата mesh из промежуточного формата .vtk (формат программы SD Slicer Version S.4) в формат .stl (формат программы MeshLab v1.23b) с использованием формата .vrml (формат программы «преобразователя» Kitware ParaViewe S.6.2);
S) используя программу MeshLab v 1.2.Sb, сегментировать и упростить полученную SD-модель в формате mesh (отделение бедренных костей человека от подложки, разделение левой бедренной кости от правой и т. д.);
4) используя разработанную программу Конвертор, преобразовать mesh формат .stl в формат поверхностей .igs и построить набор кривых и поверхностей в этой модели (преобразование формата SD изображения бедренной кости человека на основе граней в формат поверхностей);
5) используя программу Salome 5.13, выполнить обработку SD изображения бедренной кости человека (создание геометрической твердотельной модели).
На подготовленную таким образом адекватную трехмерную модель перелома бедренной кости взрослого человека накладывается расчетная сетка, задаются механические свойства металлофиксаторов и костных отломков, задаются нагрузки на бедренную кость, граничные условия и т. д., которые являются основой для дальнейшего моделирования на основе МКЭ.
В конечном итоге результатом моделирования может быть количественная оценка смещений, прочностные характеристики системы кость-фиксатор, стабиль-
ность (допустимая микроподвижность, не превышающая заданную травматологом). Весь дальнейший процесс моделирования выполняется с использованием разработанной программы «решателя», созданной на основе МКЭ и являющейся неотъемлемой частью CAE-системы.
За основу создаваемой кластерной редакции конечно-элементного решателя и других модулей ПИК был взят исходный код свободно распространяемого пакета конечно-элементного моделирования Elmer 5.5.5.
В оптимизированной кластерной версии конечноэлементного решателя (Solver^) использован проекционный вариант метода Галеркина в форме метода возможных перемещений. При этом моделирование полей упругих перемещений, деформаций и напряжений приводит к манипуляциям над т. н. матрицами жесткости и инерции (разреженными, симметричными и положительно определенными).
Движение сплошной среды, занимающей в пространстве в текущий момент времени t некоторый объем V, описывается т. н. уравнением импульсов:
pwk = pFk + У,акі, i, k = 1,2,3, (Д. 1)
где s = (стк/} , alk = ак - симметричный тензор напряжений, заданный своими контравариантными компонентами; p - плотность среды; F = {Fk}, ,
к = 1,2,3 - вектор массовых сил, действующих на индивидуальные частицы сплошной среды, заданный своими контравариантными компонентами; w = {wk},
к = 1,2,3 - ускорение индивидуальных точек сплошной среды, заданное своими контравариантными компонентами; V,, i = 1,2,3 - символ ковариантного дифференцирования в используемой системе координат (возможно, подвижной, криволинейной и деформируемой); по одноименному верхнему и нижнему индексу, как обычно в тензорном исчислении, выполняется суммирование. С точки зрения Лагранжа, свойственной механике твердого деформируемого тела, индивидуализация точек сплошной среды выполняется посредством указания их координат r с V0 в начальный момент
времени t = t0, когда среда заполняет объем V0 . При этом вектор перемещений и отсчитывается от некоторой инерциальной системы координат, и векторы скорости v и ускорения w индивидуальных частиц сплошной среды суть:
v = 3u / dt, w = d 2u / dt2. (Д. 2)
При необходимости в массовых силах F учитываются переносные и кориолиссовы силы инерции. Радиус-вектор r индивидуальной точки сплошной среды относительно используемой системы координат
задается набором криволинейных координат {xk}, к = 1,2,3 , а расстояние между двумя бесконечно близкими точками суть:
ds2 = dr2 = gk,dxkdx1, к,i = 1,2,3; gh■ = g,k . (Д. 3)
Ковариантные компоненты симметричного тензора деформаций е = {еы}, к, г = 1,2,3 связаны с ковари-
антными компонентами вектора перемещений и формулами Коши:
~^ки] ) ■
(Д. 4)
Здесь вектор перемещений задан своими проекциями на оси криволинейной системы координат, соответствующей «начальному» состоянию деформированной среды, а ковариантное дифференцирование выполняется в метрике «начального» состояния. С точки зрения Лагранжа, метрические тензоры в «начальном» gk0) и
«актуальном» gkj состоянии суть:
gkj - gíШ) + 2єкг .
(Д. 5)
Любой вектор а и тензор А могут быть заданы
как своими компонентами а
л ік
в метрике «на-
т.к
чального» состояния, так и своими компонентами а Агк в метрике актуального состояния, причем (3^ -символ Кронекера)
а2 - ~г(8к +Угик), Акг =
- А1 (3к. +У ик)(8І +Угиг).
(Д. 6)
С точки зрения Лагранжа, закон сохранения массы индивидуальной частицы среды принимает вид:
Ро\/ае1^к0)} =
т(0)
+ '2Єкі } =
( Д. 7)
где р0 , - плотность среды и компоненты метри-
ческого тензора в «начальном» состоянии; р - плотность сплошной среды в «актуальном» состоянии.
В механике деформируемого твердого тела связь тензоров напряжений (в метрике «актуального» состояния) и деформаций задается либо некоторым функционалом, учитывающим историю деформирования индивидуальной частицы сплошной среды, либо некоторыми функциями компонент тензора деформаций и компонент тензора скоростей деформаций
е = {,1} У, I = 1,2А е ,1 = & ,1 /д:
ст = ст(є, є) .
(Д. 8)
В математических моделях теории упругости кон-травариантные компоненты тензора напряжений связаны с производными объемной плотности свободной энергии деформирования (приведенной к единице не-деформированного индивидуального объема упругой среды - т. н. упругого потенциала) по ковариантным компонентам тензора деформаций
сткг -^^, и = и(е);
Ро дєк
(Д. 9)
причем предполагается, что свободная энергия деформирования (упругий потенциал) и зависит лишь от компонент тензора деформаций е и, возможно, температуры Т. Если полагать, что относительные деформа-
ции достаточно малы (порядка 10 ' формировании литых металлов), то
сткг =ди / дєкі, и -и (е).
при упругом де-
(Д. 10)
В частности, если свободная энергия деформирования и является квадратичной (положительно определенной) формой компонент тензора деформаций
и (е)-2 Скіі1 екі є а ,
СкЦЇ - с 11кг - с^к! - сШ1
(Д. 11)
то связь напряжений и деформаций становится линейной
- ск1\
1Г
(Д. 12)
Для уравнений импульса (Д. 1) требуется поставить начальные и граничные условия. Начальные условия требуются лишь в том случае, когда в левой части уравнения импульса (Д. 1) учитываются компоненты (относительного) ускорения те , и состоят в задании в начальный момент времени 10 координат и скоростей
индивидуальных точек сплошной среды, занимающей объем у. При задании граничных условий требуется
учитывать, что на некоторой части 5 (рис. Д. 1) граничной поверхности S = д¥ = 51 и S2 могут задаваться
смещения и(0) точек деформируемого твердого тела, а на другой части 52 граничной поверхности могут задаваться действующие на деформируемое твердое тело внешние поверхностные силы Р(0) .
Рис. Д. 1
То есть
(0)
ик - ик( ; (Ґ) при г є 51, к -1,2,3 :
стк пк - Р(к0) при г є 52, к -1,2,3
(Д. 13) (Д. 14)
к
и
краевые условия (Д. 14).
Вместе с тем в ряде случаев требуется отдельно рассмотреть деформирование каждой из частей ¥1 и ¥п объема ¥ = ¥: и ¥п (например, в силу существенной неоднородности среды). Полагаем: д¥ = и 53 ,
д¥п = 5П и 53, 5 = и 5П. В каждом из объемов ¥1 и ¥п должно выполняться уравнение импульсов (Д. 1), для которого на соответствующих фрагментах поверхности 51 п и 51 п 5П потребуется задать краевые условия типа (Д. 13), а на фрагментах поверхности 52 п и 52 п 5П - краевые условия типа
(Д. 14). На разделяющей объемы ¥1 и ¥п , лежащей внутри ¥ поверхности 53 должны выполняться условия:
и(1) - и(1) при г є Б3
ст(1)п і - стсщ«і при г є 53 .
(Д. 15) (Д. 16)
Далее, на поверхности контакта 5 твердых деформируемых тел I и II (рис. Д. 2) должны совпадать нормальные составляющие векторов перемещений и поверхностных сил
кі к гг **
и(і)пк- и(іі)пк приг є 5 ;
кг кг *
СТ(І)ПіПк - СТ(ІІ)ПіПк При г є 5
(Д. 17) (Д. 18)
Рис. Д. 2
При отсутствии трения на контактной поверхности на стороне I и на стороне II отсутствуют касательные составляющие векторов касательных сил, т. е. достаточно потребовать выполнения любых двух равенств из трех равенств (Д. 19) и любых двух равенств из трех равенств (Д. 20):
8(І)« - (СТуфП] )пк - 1 к -1,2,3 при г є 5**;
(Д. 19)
8(П)« - (ст(їі)піп])пк -1 к -1,2,3 при г є 5**.
(Д. 20)
В аналитической механике уравнения движения (равновесия) механических систем обычно получаются при помощи т. н. вариационного принципа возможных перемещений Эйлера-Лагранжа, согласно которому работа всех действующих на индивидуальные точки механической системы внешних и внутренних физических сил, а также сил инерции на любом виртуальном перемещении тождественно обнуляется. Под виртуальным (т. е. возможным) перемещением понимается любое бесконечно малое перемещение, совместное со связями, наложенными на механическую систему. При этом реакции связей по определению не совершают никакой работы на виртуальных перемещениях. Пусть, далее, 5и - векторное поле виртуальных перемещений точек твердого деформируемого тела. Выполнение краевых условий (Д. 13) приводит к следующему ограничению, накладываемому на поле виртуальных перемещений:
8и - 0 при г є 51.
(Д. 21)
Наличие поля виртуальных перемещений приводит к возникновению тензорного поля виртуальных деформаций, которое в метрике «актуального» состояния суть:
ЗЄкі-2(^кЧ +^8%).
(Д. 22)
Принцип Эйлера-Лагранжа применительно к механике деформируемого твердого тела принимает вид:
|[р(^к - wk)8ик -ст^в*^ +
V
+ ^ рк(0)8и кёБ - 0.
(Д. 23)
В силу инвариантности операций тензорного анализа, подынтегральные выражения в (Д. 23) могут вычисляться в любой подходящей криволинейной системе координат, однако элементы объема и площади соответствуют актуальному состоянию. Из (Д. 21) следует, что интеграл по фрагменту поверхности 51 в левой части (Д. 23) тождественно обнуляется. Если на некоторых фрагментах 5* поверхности 5 для некоторых значений индекса к = 1,2,3 ставятся краевые условия (Д. 13), а для остальных значений индекса к задаются краевые условия (Д. 14), то в поверхностном
о*
интеграле по 5 в левой части (Д. 23) выполняется суммирование по соответствующему подмножеству значений индекса к .
Из принципа Эйлера-Лагранжа (Д. 23) и произвольности поля виртуальных перемещений 8и следуют уравнение импульса (Д. 1), а также краевые условия (Д. 13) и (Д. 14). Далее, если выполняются кинематиче-
ские условия (Д. 15), из принципа Эйлера-Лагранжа (Д. 23) следует краевое условие (Д. 16). Исходя из (Д. 23), легко показать, что выполнение условия (Д. 17) влечет выполнение краевых условий (Д. 18)-(Д. 20).
Так как деформируемая среда является упругой либо гиперупругой, то из (Д. 7), (Д. 9) следует стк'5ек= = 5 и(í,)d¥0 , и уравнение (Д. 23) принимает вид:
Уравнение Эйлера-Лагранжа (Д. 23) должно выполняться при любых значениях вариаций обобщенных координат 5д(. . После группировки слагаемых, приравнивания нулю сомножителей при 5^к и преобразования соответствующих уравнений к метрике «начального» состояния находим:
|р(¥к - wk)5ukdV + 8|и(e)dV +
V.
+ | рк(0)5и^5 - 0, к -1,2,3
(Д. 24)
где ¥0 - объем, занимаемый средой в недеформиро-
ванном состоянии.
При численном решении (Д. 1) сплошную среду наделяют конечным числом степеней свободы, например, аппроксимируя поля упругих смещений отрезками Фурье-разложений по той или иной полной системе базисных функций ^к (г) = Тк (х1, х2, х3), удовлетворяющих краевым условиям (Д. 13):
ик(г, і) Чк3 (ікА (г)
1 -1
к - 1,2,3, г єV0.
(Д. 25)
При использовании метода конечных элементов индекс , в (Д. 25) сопоставляется с номером текущего узла конечно-элементной сетки, а базисные функции
^ (г) отличны от нуля лишь в пределах тех конеч-
к,
ных элементов, в которые входит текущий узел ,.
Особенностью метода конечных элементов является тот факт, что выполнение краевых условий типа (Д. 13) сводится к заданию значений обобщенных координат для узлов на граничной поверхности 5 . При этом
в (Д. 25) функции дк. (/), к = 1,2,3, , = Мк +1,...,N
становятся известными функциями времени и тем самым исключаются из числа искомых функций, а дальнейшему определению подлежат функции
(/), к = 1,2,3, , = 1,2,...,Nk . В качестве базисных
функций метода конечных элементов использовались базисные функции т. н. элементов сирендипова семейства либо базисные функции т. н. элементов Лагранже-ва семейства (где базисные функции в том или ином смысле представляют собой произведения интерполяционных полиномов Лагранжа локальных координат конечного элемента).
Поскольку варьируются лишь искомые обобщенные координаты, из (Д. 25) следует
5и
1 -1
(Д. 26)
к -1,2,3, г єV0.
(Я.к - Ъ-ик№&‘
$к](г)йУ0 +
див22 912
-ф(г )^50 +
811З22 - (зЪУ + |ро (Рк - ик)\^к.(,г)йУ0 = 0, к = 1,2,3, і = 1,2.......Ик.
(Д. 27)
Здесь величины и операции в метрике «актуального» состояния отмечены символом «~».
Из (Д. 25), (Д. 4), (Д. 8)-(Д. 11) следует, что (Д. 27) представляет собой записанную в неявной форме систему обыкновенных дифференциальных уравнений второго порядка (весьма большой размерности) относительно набора обобщенных координат
дк, к -1,2,3, а -1,2,...,Ык . Фактически (Д. 27)
представляет собой некоторый вариант проекционного метода Галеркина.
Элементы матрицы при обобщенных ускорениях С[к ,
к -1,2,3, а - 1,2,..., N к , называемой матрицей инерции,
задаются выражениями I Р0(г)gki(г)ук (г)уг (г^К, ■’V) 1 1
представляющими собой скалярные произведения по начальному объему Уд базисных функций с весовыми
коэффициентами Р0(г)gki (г) . Поскольку весовые коэффициенты образуют положительно определенную симметричную матрицу, матрица инерции является симметричной положительно определенной. Так как базисные функции отличны от нуля лишь в пределах конечных элементов, содержащих текущий узел, матрица инерции будет разреженной.
Предполагается, что сплошная среда является упругой либо гиперупругой, упругие перемещения и деформации относительно малы, и можно пренебречь (относительным) ускорением частиц сплошной среды. Следовательно,
єкі - киі +^ик):
(Д. 28)
и, с точностью до малых высшего порядка, положение упруго деформируемого тела в «актуальном» состоянии практически совпадает с положением в «начальном» состоянии, т. е. ¥ и ¥0 , 5 и 5(0) . Поскольку
компоненты вектора массовых Г - {^к } и поверхностных Р(0) - {Р(0)} сил не зависят от поля смещений
и, то (Д. 24) принимает вид:
5 Ф - 0, Ф - Ф[и] -1pFkukdV -
V (Д. 29)
- I и(e)dУ + Г рк(0)ukd5,
V 5 2
т. е. решение уравнений равновесия сплошной среды доставляет экстремальное (а фактически - минимальное) значение функционалу Ритца (Д. 30). Подстановка приближенного решения (Д. 26) в (Д. 28) и далее в (Д. 29) позволяет явно выразить величину функционала Ф через значения обобщенных координат д(0) - {дк },
к -1,2,3, 1 -1,2,...,N
3 Nk д Ф( )
Ф - Ф(q), т. е. 5Ф -УУ-т^д, . (Д. 30)
и д% 1
Из (Д. 29), (Д. 30) и произвольности вариаций обобщенных координат 5дк;- следует, что приближенному решению (Д. 25) уравнений равновесия сплошной среды соответствует выполнение требований:
= 0, к -1,2,3, 1 -1,2,...,N. . (Д. 31)
Так как деформируемая среда является линейноупругой, т. е. справедливы формулы (Д. 11), то функционал Ритца (Д. 30) будет квадратичной функцией обобщенных координат:
N 2
Ф -УУв^к -
У-1 к-1 (Д. 32)
1 NN 3 3 4 '
- "2 Аушк1дк„ Чгт ;
V-1 т-1 к-1 1 -1
Вук -1pFkVк„dУ + Iрк(0)^кvdS (Д. 33)
V 52
(суммирование по к отсутствует);
Avmn - I С^Чк^ • V1^ dV -Ата (Д. 34)
V
(суммирование по і и г отсутствует).
Матрица коэффициентов А%,тЫ квадратичной формні в правой части (Д. 32), называемая матрицей жесткости, очевидно, симметрична. Поскольку квадратичная форма в правой части (Д. 32) получается посредством интегрирования по объему V некоторой квадратичной формы, положительно определенной в каждой
точке объема V, матрица жесткости также будет положительно определенной. Поскольку базисные функции, типичные для метода конечных элементов, отличны от нуля лишь в пределах конечных элементов, содержащих текущий узел, матрица Avmkj будет разреженной. Решение исходной краевой задачи сводится к решению системы линейных уравнений:
3 N 3 N
Avmklqlm = Bvk ^ Avmklqlm . (Д 35)
l=1 m =1 l=1 m=Nl +1
с разреженной симметричной положительно определенной матрицей.
Результаты расчетов визуализируются в окне пользовательского интерфейса модуля визуализации результатов в виде модели, раскрашенной в зависимости от степени нагрузки.
ЛИТЕРАТУРА
1. Hadjicostas P.T., Thielemann F.W. The use of trochanteric slide osteotomy in the treatment of displaced acetabular fractures // Injury. 2008. V. 39 (8). P. 907-913.
2. Barr R.J., Santore R.F. Osteotomies about the hip-Adults // Chapman's Orthopaedic Surgery. 3-rd ed. Philadelphia, 2001. P. 2723-28.
3. Papagelopoulos P.J., Trousdale R.T., Lewallen D.G. Total hip arthroplasty with femoral osteotomy for proximal femoral deformity // Clin. Orthop. Relat. Res. 1996. V. 322. P. 151-162.
4. Bartonicek J., Skala-Rosenbaum J., Dousa P. Valgus intertrochanteric osteotomy for nonunion of trochanteric fractures // J. Orthop. Trauma. 2003. V. 17. P. 606-612.
5. Hammer A.J. Nonunion of subcapital femoral neck fractures // J. Orthop. Trauma. 1992. V. 6. P. 73-77.
6. Firth G.B., Robertson A.J., Schepers A., Fatti L. Developmental Dysplasia of the Hip: Open Reduction as a Risk Factor for Substantial Osteonecrosis // Clin. Orthop. Relat. Res. 2010. V. 468 (9). P. 2485-294.
7. Yang L., Jing Y., Hong D., Chong-Qi T. Valgus osteotomy combined with intramedullary nail for Shepherd's crook deformity in fibrous dysplasia: 14 femurs with a minimum of 4 years follow-up // Arch. Orthop. Trauma Surg 2010. V. 130 (4). P. 497-502.
8. Roshan A., Ram S. The neglected femoral neck fracture in young adults: review of a challenging problem // Clin. Med. Res. 2008. V. 6
(1). P. 33-39.
9. Barker K.L., Lamb S.E., Simpson H.R. Recovery of muscle strength and power after limb-lengthening surgery // Arch. Phys. Med. Rehabil. 2010. V. 91 (3). P. 384-388.
10. Paley D. Principles of Deformity Correction. New York; Berlin: Springer-Verlag; Heidelberg, 2002. P. 1-18.
11. Paliobeis C.P., Kanellopoulos A.D., Babis G.C.. Magnissalis E.A., Catling J.C., Papagelopoulos P.J., Soucacos P.N. Intrinsic passive stiffness of 2 constructs of varus proximal femoral osteotomy: external fixator or blade plate // J. Pediatr. Orthop. 2010. V. 30 (4). P. 351-356.
12. McGrory B.J., Estok DM. 2-nd, Harris W.H. Follow-up of intertrochanteric osteotomy of the hip during a 25-year period // Orthopedics. 1998. V. 21. P. 651-653.
13. Mehra A., Hemmady M.V., Hodgkinson J.P. Trochanteric non-union -does it influence the rate of revision following primary total hip replacement? A minimum of 15 years follow-up // Surgeon. 2008. V. 6
(2). P. 79-82.
14. Seki T., Hasegawa Y., Masui T., Yamaguchi J., Kanoh T., Ishiguro N., Kawabe K. Quality of life following femoral osteotomy and total hip arthroplasty for nontraumatic osteonecrosis of the femoral head // J. Orthop. Sci 2008. V. 13 (2). P. 116-121.
Поступила в редакцию 14 мая 2012 г.
Norkin I.A., Yamshchikov O.N., Markov D.A., Abdulnasi-rov R.K., Peregorodov D.N., Zaigralov A.Yu. SOME ASPECTS OF PROXIMAL FEMUR COMPUTER MODELING
The article provides detailed information on the implementation of computer simulation in the treatment of injuries and diseases of the femur.
Key words: computer modeling; femur; personal computer; 3D-model; surface; mechanics.
