Нелинейная система масс-с-пружинками для моделирования больших деформаций мягких тканей Текст научной статьи по специальности «Электротехника, электронная техника, информационные технологии»

5

УДК (517.958+539.3):004.45

НЕЛИНЕЙНАЯ СИСТЕМА МАСС-С-ПРУЖИНКАМИ ДЛЯ МОДЕЛИРОВАНИЯ БОЛЬШИХ ДЕФОРМАЦИЙ МЯГКИХ ТКАНЕЙ С.Н. Николаев

Описывается операция помещения имплантата под мягкие ткани. Показывается, что в такой операции ткани могут достигать растяжений, при которых проявляются нелинейные свойства. Разрабатывается модификация модели масс-с-пружинками, позволяющая моделировать нелинейную деформацию растяжения. Описывается подход к построению модуля упругости с помощью сплайнов. Моделирование коэффициента Пуассона производится за счет различных же-сткостей для различных типов пружин в кубической решетке. Для определения значений жесткостей решается система уравнений, описывающая растяжение материала. Модель проверяется с помощью эксперимента по растяжению квадратного образца материала. Тестирование показывает, что величина растяжения под действием различных внешних сил совпадает с задаваемым нелинейным модулем упругости. Также коэффициент Пуассона моделируется с погрешностью не более 35%, что лучше результатов для имеющегося способа моделирования коэффициента.

Ключевые слова: модель масс-с-пружинками, растяжение материала, нелинейный модуль упругости, коэффициент Пуассона, графики мягких тканей.

Введение

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

Модель конечных элементов является наиболее полной и качественной моделью для решения биомеханических задач. Однако врачам необходимо как можно быстрее увидеть результат заданной операции, поэтому одним из ключевых моментов также является скорость моделирования. Конечные элементы из-за использования сложных вычислений работают медленно. В связи с этим сейчас наряду с ними широкое распространение получает новая модель масс-с-пружинками [1], основным критерием которой является скорость вычислений. Эта модель уже была успешно использована для моделирования различных задач биомеханики сердца [2, 3], легких [3] и мышц [4].

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

51 = п х Я,

где Я - радиус имплантата. Площадь после растяжения будет равна площади полусферы:

52 = 2 х п х Я .

Таким образом, относительное удлинение тканей может достигать величины 1,0. Известно, что при таких растяжениях мягкие ткани проявляют нелинейные свойства, которыми нельзя пренебрегать [6].

В обычной модели масс-с-пружинками используются линейные формулы пружин [1]. В ряде работ описано улучшение этой модели, в которой пружины представлены в виде полиномов второго и третьего порядков [2] или как кусочная линейная функция [7]. Однако формы нелинейного поведения тканей значительно сложнее описаний полиномов, а линейные функции приближают графики тканей лишь на границах секций. Также при растяжении ткани происходит ее сжатие в плоскости, перпендикулярной направлению растяжения, которое описывается коэффициентом Пуассона. Его моделируют с использованием генетического алгоритма [8] или аналитического подхода [9]. Для генетического алгоритма необходимо иметь различные измерения опытного образца, что не всегда бывает возможно. В аналитическом подходе для диагональных пружин используется формула сдвига, которая правомерна лишь для малых деформаций. Таким образом, необходим новый подход к моделированию нелинейной упругости и коэффициента Пуассона. Целью настоящей работы является разработка модификации модели масс-с-пружинками для моделирования растяжения тканей при величинах относительного удлинения от 0,0 до 1,0.

МЕХАНИКА И МЕХАТРОНИКА

Моделирование нелинейности

В модели масс-с-пружинками объект представляется как набор точек, описывающих массу и положение объекта, соединенных пружинами (рис. 1).

Рис. 1. Модель масс-с-пружинками

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

/ = к х (| *12| - ^ х-^.

Ы

Здесь /- сила, действующая на одну из вершин, к которой пружина прикреплена (сила, действующая на вторую вершину, та же по величине, но противоположна по направлению); х12 - текущее расстояние между вершинами; /12 - начальная длина пружины; к - коэффициент жесткости пружины. Стоит отметить, что формулу для расчета пружин можно записать в другом виде:

/ = Ба х Е х ——, Ь

где Е - модуль Юнга; Бд - размер топологического элемента; Ь - начальная длина пружины; ДЬ - абсолютное удлинение пружины. Эта форма записи более предпочтительна, поскольку она справедлива для пружины любого размера. В нелинейном случае эта формула преобразуется в более общую:

/ = Бд х Е/ (ДЬ),

где Е/- это некоторая функция, связывающая напряжение и деформацию.

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

Моделирование коэффициента Пуассона

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

—реб.

Рис. 2. Растяжение квадратной решетки В поперечных направлениях произошло сжатие на величину

ДЬ = е

= —V х е„

'-Ьреб =— V хДЬр

и г — V ^ О ^ г — V ^ \ I г

реб. раст реб. реб.

Здесь V - коэффициент Пуассона. В работе [9] авторы определяют относительные коэффициенты пружин по формулам сдвига и растяжения. Однако стоит ограничиться только операцией растяжения, поскольку пружины не моделируют сдвиг. Построенная система уравнений будет выглядеть следующим образом:

Р АЬ б АЬ б АЬ

-— = Е/ = к х Е/ + С х Е/ (-=

реб. реб. реб. диаг.

)

Ьреб.+АЬреб.

Ь + АЬ

диаг. диаг.

К2 х В/ (-

- V х АЬ

'реб.

Ь

) + С х В/ (-

'реб.

АЬ Ь б - V х АЬ б

• ^^диаг. ) х ре°._Реб_ = 0

Ь™ Ь + АЬ „„ " '

(1)

где к\ - это неизвестный коэффициент при реберных пружинах, а С - при диагональных, к2 описывает случай, когда реберная пружина сжимается.

Выражая С из второго уравнения системы (1) и подставляя его в первое уравнение, получим

Е/ (■

АЬ б АЬ б - V х АЬ ,

реб. ) = к, х Е/(—- к2 х Е/(- реб.

Ьр

Ьр

Ьр

^реб.+А^реб.

Ьреб. V хАЬреб.

(2)

^реб. реб. реб.

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

./реб.

(

;

Е/ (-

АЬ,

Е/ (-

- V х АЬ,

реб.

Л

реб.

) +

Ь

реб.

реб.

(1 + V) х Ь

реб.

Ьреб.+АЬреб.

- V

' АЬреб. — 0;

(3)

АЬреб.

^ х Е/ (у*)

Ьреб.

, АЬре6.< 0.

АЬд

и Ьдиаг. описывают в системе (1) длины диагональных пружин. По теореме Пифагора они

определяются следующими формулами:

Ь = Ь б

диаг. реб.

<>/2 ;

(Ьдиаг. + АЬдиаг. ) =

(4)

(5)

АЬреб. можно выразить через известные длины диагональных пружин из формул (4), (5):

АЬреб =-

Ьд

где J определяется как

./ = Ьр,

V -1 + ,

:(1+V2))

J =-

Г Ь +АЬ ^

диаг. диаг.

Ьдиаг

- (1 + 2v + V2)

(1 + V2)

При подстановке АЬреб. во второе уравнение системы (1) получается следующее выражение (к2 уже задано равным 1,0):

АЬ -л/2 х Е/ (-V х J)х (Ь +АЬ )

С х Цр( диаг. ) =__/ ^ диаг._диаг./

1 Ьдиаг. (1-V х '

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

/ =

диаг.

х Е/ (-V х J )х (Ьдиаг. + АЬдиаг. )

(1 - V х•/ ))

АЬдиаг.

^ х Е/ (■ диаг.

Ь

)

, АЬ — 0;

диаг.

, АЬдиаг. < 0.

(6)

Полученные формулы подставляются вместо формул обычных линейных пружин. Обращаясь еще раз к формулам (3) и (6) для расчета сил, необходимо отметить, что в новой модели в качестве искомых данных задаются функция Е/, описывающая нелинейный модуль упругости, и численное значение коэффициента Пуассона.

Результаты моделирования операции растяжения

После построения модели были проведены эксперименты по растяжению различных материалов тканей. Для них брался квадратный образец с размером ребра 4 см, который выравнивался по осям и разбивался на решетку с размером ребра 1 см. Поскольку данная модель разрабатывалась в первую очередь для описания растяжения тканей над имплантатом, в качестве исходных данных были использованы модули упругости мягких тканей. В простейшем случае имплантат помещается под две ткани - кожу и жир, поэтому именно они были выбраны для эксперимента. Отдельно для тканей по задаваемой функции Ef для различных относительных удлинений была рассчитана величина напряжения. Результаты представлены на графиках 1 и 2 (skin Ef и adipose 'Ef) .

Эксперименты на растяжение проводились следующим образом. К образцу вдоль одной из осей прикладывались две силы растяжения, направленные в противоположные стороны и действовавшие на каждую граничную точку по выбранной оси. Силы строились так, чтобы на каждые две вершины каждого граничного квадратного элемента действовали силы одинаковых величин. Другими словами, вершины, находящиеся на границах различных квадратных элементов, растягивались с удвоенными силами. При таком подходе к построению сил среднее напряжение, действующее на образец, можно вычислить как уполовиненную величину задаваемой силы. После 15000 итераций деформации методом ограничивающей рамки измерялись размеры растянутого образца по двум осям. По начальным и конечным размерам вычислялись относительные удлинения в направлении растяжения и в поперечном направлении.

Эксперимент повторялся несколько раз. Для каждого рассматриваемого материала деформации растяжения начинались с воздействия малых сил. Затем в последующих опытах значения сил увеличивались. Опыты проводились до тех пор, пока относительное удлинение при растяжении не начинало превышать 1,0.

Первой была рассмотрена кожа. Ее нелинейный модуль упругости при относительных деформациях от 0,0 до 1,0 условно делится на три участка [10]. При небольших деформациях он постоянен. Затем после удлинения 0,4 он начинает плавно возрастать, пока снова не станет постоянным примерно при удлинении 0,7. Таким образом, функцию Ef было решено разбить на три участка. Первый и третий куски описываются функциями первого порядка. Средняя часть представляет собой кубический сплайн, соединяющий линейные фрагменты. На рис. 1 представлены значения функции Ef (skin Ef) и экспериментальные результаты зависимости напряжения и относительного растяжения для различных заданных коэффициентов Пуассона. На рис. 3 показаны экспериментальные результаты зависимости относительного сжатия материала и его растяжения при различных заданных коэффициентах Пуассона. При этом по горизонтальной оси указаны значения растяжений, а по вертикальной - значения сжатий.

skin 'Ef v=0,5 v=0,3 v=0,1

Рис. 1. Экспериментальные результаты измерения нелинейного модуля упругости кожи при различных значениях коэффициента Пуассона. skin 'Ef - значения задаваемого нелинейного модуля упругости

6400 3200 1600

g 800

С

400 200 100 50

0 0,2 0,4 0,6 0,8 1 е adipose 'Ef v=0,5 v=0,3 v=0,1

Рис. 2. Экспериментальные результаты измерения нелинейного модуля упругости жира при различных значениях коэффициента Пуассона. adipose 'Ef - значения задаваемого нелинейного модуля упругости

е растяжения

0 0,5 1

—„...¥=0,5 ¥=0,3 ¥=0,1

Рис. 3. Экспериментальные результаты измерения коэффициентов Пуассона для кожи

0

-0,05 -0,1 ! -0,15

I

3 -0,2

-0,25 -0,3 -0,35 -0,4

е растяжения 0,4 0,6 0,8 1

Чх

—— Я -V

¥=0,5 ¥=0,3 ¥=0,1

График 4. Экспериментальные результаты измерения коэффициентов Пуассона для жира

Свойства жира описаны в [11]. При небольших деформациях модуль упругости жировой ткани линеен. Затем он постепенно уменьшается, и при относительном удлинении, большем 0,6, напряжение становится постоянным и больше не зависит от удлинения материала. Отмасштабированная функция арктангенса наиболее точно приближает модуль упругости жира. На рис. 2 представлены результаты значений нелинейного модуля упругости (adipose Ef) и экспериментальных отношений напряжения и деформации для различных значений коэффициента Пуассона. На рис. 4 показаны экспериментальные зависимости относительного сжатия жира от его растяжения. При этом по горизонтальной оси указаны значения растяжений, а по вертикальной - значения сжатий.

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

Из рис. 3 и 4 видно, что в среднем отношение между деформациями сжатия и растяжения меньше задаваемого коэффициента Пуассона. При увеличении коэффициента Пуассона погрешность также увеличивается, однако относительная погрешность не превышает 35%.

Подводя итог, можно отметить, что экспериментальные данные для предложенного подхода повторяют формы нелинейных графиков упругости, в отличие от подходов из [2, 7]. Кроме того, моделирование коэффициента Пуассона происходит с меньшей погрешностью, чем в [9]. Поскольку большинство мягких тканей имеет модули упругости, аналогичные коже и жиру, и коэффициент Пуассона в пределах от 0,35 до 0,5, данный способ с хорошей точностью будет описывать растяжения этих тканей. Таким образом, данный подход можно использовать для моделирования растягиваемых тканей при помещении под них имплантатов или эспандеров.

Заключение

В результате проделанной работы была создана модификация модели масс-с-пружинками, которая позволяет моделировать нелинейный модуль упругости и коэффициент Пуассона для относительных деформаций в диапазоне от 0,0 до 1,0. Были проведены тесты для двух формул и различных значений коэффициентов Пуассона, которые показали удовлетворительные результаты.

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

В дальнейшем планируется рассмотреть случай трехмерного пространства для описания деформаций объемных моделей.

Литература

1. Kass M. An introduction to Continuum Dynamics for Computer Graphics // SIGGRAPH Course Notes: Physically-based Modelling. - ACM SIGGRAPH, 1995. - 13 p. [Электронный ресурс]. - Режим доступа: http://www.cs.cmu.edu/~baraff/pbm/continuators.pdf, свободный. Яз. англ. (дата обращения 28.07.2013).

2. Jarrouse O. Modified Mass-spring System for Physically Based Deformation Modeling. - Karlsruhe: KIT Scientific Publishing, 2012. - 240 p.

3. Wang Y., Xiong Y., Xu K. A Mass-Spring Model for Surface Mesh Deformation Based on Shape Matching // Proceedings of the 4th international conference on computer graphics and interactive techniques in Australia and Southeast Asia. - GRAPHITE06, 2006. - P. 375-380.

4. Nedel L., Thalmann D. Real Time Muscle Deformations Using Mass-Spring Systems // Proceedings of Computer Graphics International. - CGI, 1998. - P. 156-165.

5. Николаев С. Программный модуль для трехмерного моделирования хирургической операции по увеличению груди // Компьютерные инструменты в образовании, Центр Информатизации образования «КИО», 2012. - № 3. - С. 38-47.

6. Fung Y. Biomechanics: Mechanical Properties of Living Tissues. - New York: Springer-Verlag, 1993. -592 p.

7. Basafa E., Farahmand F., Vossoughi G. A Non-linear Mass-spring Model for More Realistic and Efficient Simulation of Soft Tissue Surgery // Studies in Health Technology and Informatics. - IOS press, 2008. -P. 23-25.

8. Bianchi G., Solenthaler B., Szekely G., Harders M. Simultaneous Topology and Stiffness Identification for Mass-Spring Models based on FEM Reference Deformations // Medical Image Computing and ComputerAssisted Intervention - MICCAI 2004. - Lecture notes in computer science. - Springer, 2004. - V. 3217. -P. 293-301.

ЗАВИСИМОСТЬ РЕАКТИВНОГО СОПРОТИВЛЕНИЯ ПЬЕЗОЭЛЕКТРИЧЕСКОГО.

9. Baudet V., Beuve M., Jaillet F., Shariat B., Zara F. Integrating Tensile Parameters in Hexahedral Massspring System for Simulation // Laboratoire d'InfoRmatique en Image et Systemes d'information. - Research report, 2007. - 8 p. [Электронный ресурс]. - Режим доступа: http://liris.cnrs.fr/Documents/Liris-4449.pdf, свободный. Яз. англ. (дата обращения 28.07.20i3).

10. Hendriks F. Mechanical behavior of human epidermal and dermal layers in vivo. - The Netherlands: Technische University Eindhoven, Eindhoven, 2005. - 119 p. [Электронный ресурс]. - Режим доступа: http://alexandria.tue.nl/extra2/2005i094i.pdf, свободный. Яз. англ. (дата обращения 28.07.2013).

11. Geerlings M. Skin layer mechanics. PhD thesis. - The Netherlands, Universiteitsdrukkerij TU Eindhoven, Eindhoven, 2009. - 122 p. [Электронный ресурс]. - Режим доступа: http://alexandria.tue.nl/extra2/657803.pdf, свободный. Яз. англ. (дата обращения 28.07.2013).

Николаев Сергей Николаевич - Санкт-Петербургский государственный университет, аспирант,

ser.niev@gmail.com

УДК 53.091

ЗАВИСИМОСТЬ РЕАКТИВНОГО СОПРОТИВЛЕНИЯ ПЬЕЗОЭЛЕКТРИЧЕСКОГО ПРЕОБРАЗОВАТЕЛЯ ОТ МЕХАНИЧЕСКИХ

ПАРАМЕТРОВ ЕГО НАГРУЗКИ И.П. Попов

Показано, что инертная нагрузка пьезоэлектрического преобразователя может быть представлена в виде индуктивного сопротивления в его цени питания, а упругая - в виде емкостного. Предложены модели колебательных систем с элементами различной физической природы, в которых могут возникать свободные гармонические колебания. Показано, что в инертно-емкостной (mC) колебательной системе происходит взаимное превращение энергии электрического ноля конденсатора в кинетическую энергию массивного элемента; в упруго-индуктивной (kL) колебательной системе происходит взаимное превращение энергии магнитного ноля соленоида в потенциальную энергию пружины.

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

Введение

При разработке мехатронных систем следует учитывать, что механическая нагрузка может быть представлена в виде реактивного сопротивления в цепях питания электромеханических преобразователей. В индуктивных преобразователях инертная нагрузка создает емкостное сопротивление, а упругая -индуктивное [1]. В настоящее время получили распространение пьезокерамические электромеханические преобразователи [2], особенно в робототехнике [3]. Актуальной задачей является выявление характера реактивного сопротивления (их цепей питания), в виде которого представлена механическая нагрузка. В литературе нет непосредственного решения этой задачи. Предпосылкой ее решения является одна из двух систем аналогий между электромагнитными и механическими величинами [4], в соответствии с которыми масса и упругость связаны дуальными соотношениями с индуктивностью и емкостью: m ^ L, k ^ C.

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

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

Пьезоэлектрический преобразователь с инертной нагрузкой в электрической цепи

На рис. 1 изображен пьезоэлектрический преобразователь с инертной нагрузкой массой m. Работа преобразователя основана на прямом и обратном пьезоэффектах. Прямой пьезоэффект проявляется в том, что на обкладках пьезоэлемента при его деформации x появляется электрический заряд q:

q = d1x, (i)

где di - пьезомодуль. При подаче на обкладки напряжения u пьезоэлемент деформируется и развивает усилие F:

F = d2u. (2)

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

Пусть на обкладки пьезоэлемента подается напряжение u. В соответствии с третьим и вторым законами Ньютона, а также с учетом (2)