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

УДК 531/534:57+612.7

ОПТИМИЗАЦИЯ ДВУХЛИСТОВОГО УПРУГОГО ЭЛЕМЕНТА ПРОТЕЗА СТОПЫ С ИСПОЛЬЗОВАНИЕМ ЛИНЕЙНОЙ И НЕЛИНЕЙНОЙ ТЕОРИИ ИЗГИБА

О.И. Брынских, М.А. Осипенко, Ю.И. Няшин

Пермский государственный технический университет, кафедра теоретической механики, 614600, Россия, г. Пермь, Комсомольский пр., д. 29a, e-mail: oma@theormech.pstu.ac.ru

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

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

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

Введение

Многие современные конструкции протеза стопы содержат упругий элемент, который представляет собой консольно закрепленную многолистовую рессору с зазорами между листами или без них (см. рис. 1: a - серийный образец, изготовленный на экспериментальном заводе средств протезирования, г. Реутово, Московская область; б - опытный образец упругого элемента, изготовленный в Пермском НИИ композиционных материалов из углепластика КМУ-4Л; в - упругий элемент протеза стопы Carbon Copy, изготовленный фирмой Ohio Willow Wood Company, США [1]; г -упругий элемент протеза стопы Quantum Foot, изготовленный фирмой Hanger Inter Med, Великобритания). Основная проблема, возникающая при конструировании упругого элемента, - обеспечение заданных упругих характеристик и одновременно максимальной прочности.

Упругие элементы указанного вида сейчас активно изучаются, в том числе теоретически, как в России, так и за рубежом [1-9]. Однако полной и адекватной математической модели упругого элемента не имеется и в ближайшее время такая модель, вероятно, не будет создана ввиду значительной сложности конструкции. Разрабатываются упрощенные модели, отражающие поведение упругого элемента в некоторых частных случаях.

В данной работе за основу взята математическая модель, разработанная на кафедре теоретической механики Пермского государственного технического университета [7, 8]. В эту модель внесены новые элементы - учтены зазор между

© С.И. Брынских, М.А. Осипенко, Ю.И. Няшин, 2003

Рис. 4. Функция, описывающая форму изогнутого листа

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

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

Используемая модель упругого элемента показана на рис. 2. В этой модели упругий элемент представляет собой два плоских (в отсутствии нагрузки) листа прямоугольного сечения, отстоящих друг от друга на заданном расстоянии X; один край каждого листа защемлен, другой свободен. Листы изготовлены из материала с модулем Юнга Е и имеют одинаковые ширины н; длины 11 > 12 и толщины к1, к2 листов различны (нижний индекс 1 или 2 здесь и далее означает величины, относящиеся, соответственно, к нижнему и верхнему листам). К краю нижнего листа, перпендикулярно ему, приложена заданная следящая сила Е, равномерно распределенная по ширине листа (рис. 3). Эта сила моделирует реакцию поверхности при ходьбе на протезе. Под действием силы Е листы испытывают совместный изгиб (рис. 3; изгиб не предполагается малым). Геометрически толщины листов считаются равными нулю; реальные толщины входят только в изгибные податливости листов (см. ниже). Изгиб считается цилиндрическим, так что формы листов полностью описываются их профилями. Предполагается, что контакт профилей листов либо происходит только в одной точке (расположенной на краю верхнего листа), либо вообще отсутствует (при достаточно малой величине силы Е). Трением между листами пренебрегаем. Профиль изогнутого листа длины I (не меняющейся при изгибе) описывается функцией р(*), где 0 < * < I - длина дуги профиля листа от точки защемления до некоторой произвольной точки (криволинейная координата этой точки), р - угол между касательными к профилю в этих точках (рис.4).

Предполагается [10], что в точке с криволинейной координатой * на левую (относительно этой точки) часть листа со стороны правой части действует (сосредоточенный) момент

М = р'(*)/а, (1)

где а = 12/Енк3 - изгибная податливость листа, а “напряжение” (ст** на нижней поверхности листа) в этой точке

Ек

Ф) = — р С*). (2)

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

= шахі тах ст1(5),шах ст2(£) I (3)

' 0<«<// 0<^</2

^тах

(максимальное напряжение в упругом элементе). При заданных Е, н , /1, Е величины 8 и сттах являются функциями /2, к1, к2:

5 = 8(/2, к1, М, (4)

СТтах = СТтах(/2, к1, к2). (5)

Как упоминалось во введении, размер /1 упругого элемента и его податливость (прогиб 80 при определенной силе Е) заранее заданы; требуется максимизировать

прочность (минимизировать сттах). Тогда математическая постановка задачи

оптимизации принимает следующий вид.

Задача. Для заданных Е, н, /1, Е, 80, найти /21, к*, к^, минимизирующие

СТтах (/2 , к1, к2) при условии

8(/2, к1, к2) = 80. (6)

Алгоритм численного решения задачи оптимизации

Сформулированная задача оптимизации решалась численно. Вначале был построен алгоритм вычисления функций (4) и (5). Опишем этот алгоритм.

Если контакт между листами отсутствует, то верхний лист остается плоским, а на нижний действует одна сила Е (рис. 5). Составляя с помощью (1) уравнение равновесия участка АВ нижнего листа (рис. 5), найдем

где

- pl(s)/a1 + FH (s) = 0,

H (s) = | cos(p1 (/1) - p1 (x))dx

(7)

(8)

- плечо силы Е относительно точки А. Подставляя (8) в (7), дифференцируя полученное соотношение по * и учитывая, что р1 (0) = 0 (условие защемления) и р1'(/1) = 0 (следует из (7) и (8)), получим краевую задачу для р1(*) :

Pl(s) = -a1 Fcos(p1 (/1) - р1 (s)), 0 < s < /1,

(9)

< P(0) = 0,

PM = 0

Задача (9) неудобна для численного решения, так как является краевой (а не задачей Коши) и, кроме того, содержит неизвестное заранее p1 (l1). Поэтому была введена вспомогательная функция

/1(s) = Pi (li) - Pi(li - s) (10)

(геометрический смысл этой функции - отсчитываемый от свободного края листа угол между касательными к профилю листа в точке на этом краю и в точке с криволинейной

координатой l1 - s ). Из (9) и (10) получим задачу Коши для /1(s) :

/1(s) = a1 F cos // (s), 0 < s < l1,

/1(0) = 0, (11)

/1 (0) = 0.

Решив задачу (11), найдем из (10)

Pl(s) = /l(/1)-/l(/1 - s).

(12)

/

/

/

Рис. 5. Схема изгиба листов при отсутствии контакта

Так как предполагалось, что верхний лист остается плоским, то

Piis) - 0. (13)

Это предположение справедливо, если (как следует из простых геометрических соображений)

s

|sinpj (x)dx < X

0

при всех s > 0, таких что

s

|cos (p1 (x)dx < l2 . (14)

0

Если условие (14) не выполнено, то между листами происходит контакт. Схема изгиба нижнего и верхнего листов (соответствующая сформулированному во введении предположению о характере контакта) приведена на рис. 6. На нижний лист действуют силы F и P, на верхний лист действует сила P' = -P , причем величина силы P, криволинейная координата lP точки контакта и угол а между силой P' и нормалью к профилю верхнего листа в точке контакта (он же угол между профилями листов в точке контакта) заранее неизвестны. Составляя уравнения равновесия участков листов, расположенных между свободным краем и точкой с криволинейной координатой s и вводя вспомогательную функцию

W2(s) = P2(l2) - P2(l2 - sЬ (15)

получим, аналогично (11), задачу Коши для у/1(s), у/2(s):

/1(s) =

a1 F cos/^s), 0 < s < lr - lP, аГF cos у/Г (s) - arP cos/ (s) - у/Г (lr - lP )), lr - lP < s < lr,

/ (s) = a2P cos(/2(s) + a), 0 < s < l2,

/i(0) = / 2 (0) = 0,

/1(0) = /2(0) = 0.

(16)

Функции р1(^) и р2(^) выражаются через решение задачи (16) по формуле (12) и по аналогичной формуле

Р2^) = ^2(/2) ~¥2(12 - 5)- (17)

Эти функции зависят от параметров P, lP

a

(Pi = P(P, Ip ,a, s),

(p2 = p2(P,lP,a,s) . Параметры должны быть далее найдены из: а) условия совпадения значений координаты х точки контакта (рис. 6) при их вычислении с использованием сначала рГ (s), а затем р2 (s); б) аналогичного условия для координаты у; в) из упоминавшегося выше условия совпадения угла между профилями листов в точке контакта с углом a . Эти условия представляют собой систему трех алгебраических уравнений относительно P , lP , a :

Ip l2

I cos рГ (P, lP, a, s)ds -1 cos p2 (P, lP, a, s)ds = 0,

0 0

Ip l2

|sin рГ (P, lP, a, s)ds -1sin p2 (P, lP, a, s)ds - Я = 0,

(18)

рГ (P, lP, a, lP) - p2 (P, lP, a, l2) - a = 0.

Окончательно алгоритм вычисления функций (4) и (5) выглядит следующим образом. Сначала численно (методом Рунге - Кутта [11]) решается задача (11). Затем по формуле (12) находится р1(^) и проверяется условие (14). Если это условие выполнено, то для р2(^) справедливо (13). После этого находится 5

Т

5 =| sin рГ (s)ds

(19)

и из (12), (13), (2) и (3) находится <гшах. Если условие (14) не выполнено, то численно решается совместно задача (16) (методом Рунге - Кутта) и, с учетом (12) и (17), система

(18) (методом Ньютона с численным нахождением производных [11]). После этого из

(19) находится 5 и из (12), (17), (2) и (3) находится <гшах .

Далее численно находился минимум функции (5) при условии (6). Для этого сначала условие (6) численно разрешалось (методом Ньютона) относительно 12 :

12 = /2(К ^2). (20)

Подставляя (20) в (5), получим (вообще говоря, негладкую) функцию двух переменных (А1, ^2), точка безусловного минимума которой (, И22) находилась численно (методом Нелдера - Мида [12]). Далее из (20) получаем оптимальное значение /2= 12(А1, А2) и из (5) получаем минимальное значение отах:

(/2, а; , %).

^max ^max

а, кГ/мм2 а, кГ/мм2

А А

Рис. 7. Напряжения в листах оптимального (а) и неоптимального (б) упругих элементов

Приведем пример результата реализации описанного алгоритма численного решения задачи оптимизации. Все дальнейшие численные результаты соответствуют

Е - 9533 кГ/мм2 (углепластик КМУ-4Л), ^ - 60 мм , Х = 8 мм, 11 - 180 мм,

Г - 1,5 • 80 кГ. При 50 - 100 мм описанный алгоритм дает Л* = 3,08 мм, Л2 = 2,70 мм,

/2 - 89,5 мм, атах - 112 кГ/мм2. На рис. 7а показаны зависимости сг1(я) и а2(я) для упомянутых оптимальных значений параметров. Для сравнения, на рис. 7б показаны те же зависимости для неоптимальных (но соответствующих 5 - 50) значений параметров

к1 - 3,25 мм, Н2 - 2,00 мм , /2 - 120 мм; для этих значений атах - 145 кГ/мм2. Другие примеры приведены ниже (при сравнении результатов, полученных при использовании нелинейной и линейной теорий изгиба).

Аналитическое решение задачи оптимизации с использованием линейной теории

изгиба

Проведем линеаризацию соотношений (2), (3), (11), (12), (14), (16)-(19) по параметрам Г и X. Получающиеся приближенные варианты этих соотношений соответствуют линейной теории (слабого) изгиба [10]. Решения (11), (16) и (18) тогда находятся аналитически; после простых преобразований получаем (символом ~ обозначены величины, найденные в линейном приближении):

(1) если X > Хс, то

3(12, к1, к2) = < г/Буй;5,

С?! (5) = 6(/1 - 5)Г,

<~2( 5) = 0;

(и) если X < Хс, то

3(/2,к1,к2) = 2(2/3Г - /22(3/1 - /2)Р)/Е^к\ , (21)

с?(5) = /б((/1 - 5)Г - (/2 - 5)Р)/^2, 0 < 8 < /2, (22)

1 { 6(/1 - 5) Г/^2, / 2 < 5 < 11,

<г2(5) = 6(/2 - 5)Р/м>к^ , (23)

где

~ = (3/1 - /2)г/а3 - ХБУ2/2. (24)

2/2 (V к + V к23 ) ’

Хс = 2/22 (3/1 - /2)Г/Бук! ; (25)

50, мм 70 90 50, мм

Рис. 8. Зависимости параметров оптимального упругого элемента от заданного прогиба 30 Зависимости в линейном приближении отмечены символом ~

сгтах (/ 2, к1, к2) далее вычисляется по формуле (3).

Построение аналитического решения задачи оптимизации, даже для рассматриваемого линейного приближения, сопряжено с весьма громоздкими вычислениями [7]. Поэтому используем здесь предложенный в [7] эвристический метод решения (используемый только в линейном приближении; строгое обоснование метода имеется только для X = 0). Этот метод состоит в том, что условие минимизации <гтах (/2, к1, к2) заменяется на условие

<г 1 (5) = <г2 (0) для 0 < 5 < /2 (26)

(оптимальный упругий элемент должен быть равнонапряженным). Заметим, что сама возможность выполнения условия (26) существует только в линейном приближении и только в случае (и). В этом случае условие (26), с учетом (22) и (23), равносильно

равенствам Г = Р и <г1(0) = <г2 (0). Добавляя эти равенства к условию 3 (/2, к1, к2) = 30

и используя (21)-(24), получим систему алгебраических уравнений для /2*, к1*, к2 :

2 Г (2/3 - /22 (3/! - / 2) )/ Еу^3 = 30,

(3/1 - /2)Г/к/ - ХЕы/2/22

2/2 (V к.3 + V к2) ’

(/1 -/2)1 к2 = /г!к2 .

Несложное исследование этой системы приводит к выводу, что при 30 > X (это неравенство будем считать выполненным, так как обычно 30 >> X) она имеет единственное (положительное) решение, которое можно представить в виде:

Г=

~*= /1 (2 Г (2 - ц\3 - ц))/( Еу80) )

13

(27)

Ь*=уЬ*, (28)

'2

*

// = М, (29)

где

!Л = у2/(1 + у2), (30)

V - (положительный) корень уравнения

3(к -1>4 + 2у3 + 6ку2 + 2к = 0, (31)

к = . (32)

Для V можно найти явное выражение, которое здесь не приводим ввиду его громоздкости. Из (22) и (29) получаем минимальное значение с~тах :

( х2тг2гл .лр У/3

. (33)

~ * 6 ^гпах = -

А

52 Е 2(1 - ц) Е

4(2 + 2ц-ц2)2 у

Формулы (27)-(33) дают аналитическое решение задачи оптимизации с использованием линейной теории изгиба. Заметим, что при Л = 0 из (28)-(32) получается известный результат [7]: /2*Д = 4/13, А2*/А1* = 2/3 .

Сравнение результатов оптимизации, полученных с помощью нелинейной и

линейной теорий изгиба

На рис. 8 показаны зависимости оптимальных толщин листов, оптимальной длины верхнего листа и минимизированного максимального напряжения от заданного прогиба 30. Каждая зависимость приведена для нелинейной теории изгиба и в

линейном приближении. При 30 = 100 мм погрешности величин к1*, А2*, /2*, с~тах по

отношению к А* , А2* , /2*, ^ах равны, соответственно, 5,8, 16, 2,7 и 6,3%.

Следовательно, линейное приближение дает достаточную точность при вычислении /2*,

7 * * 7 *

удовлетворительную - при вычислении й1 и <Гтах и непригодно при вычислении Й2 . Следует, однако, заметить, что качество линейного приближения определяется не только приведенными погрешностями, но и другими факторами. Во-первых, линейное

приближение дает заниженное значение ^тах, что нежелательно вблизи предела

прочности материала (94 кГ/мм для углепластика КМУ-4Л). Во-вторых, если с помощью нелинейной теории подсчитать 3 для параметров, считающихся оптимальными согласно линейной теории, то эта величина будет заметно отличаться от 50 . Например, при 50 = 100 мм описанный выше численный алгоритм дает 5 = 86,0 мм

(для Й1 = А1 = 3,26 мм, Й2 = /г2* = 3,15 мм, /2 = ~2* = 87,1 мм ).

Выводы

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

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

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

1. Geil M.D., Parnianpour M., Berme N. Significance of nonsagittal power terms in analysis of a dynamic elastic response prosthetic foot // Journal of Biomechanical Engineering. 1999, Vol. 121. P. 521-524.

2. GeilM.D. An iterative method for viscoelastic modeling of prosthetic feet // Journal of Biomechanics. 2002. V. 35. P. 1405-1410.

3. Allard P., Trudeau F., Prince F., Danserau J. Labelle H., Duhaime M. Modelling and gait evaluation of asymmetrical-keel foot prosthesis // Med. & Biol. Eng. & Comput. 1995. Vol 33. P. 2-7.

4. Barr A.E., Siegel K.L., Danoff J.V., McGarvey C.L.,Tomasko A., Sable I., Stanhope S.J. Biomechanical comparison of the energy storing capabilities of SACH and Carbon Copy II prosthetic feet during the stance phase of gait in a person with below-knee amputation // Physical Therapy. 1992. V. 72. № 5. P. 344-354.

5. Брынских С.И., Осипенко М.А. Численное исследование и оптимизация двухлистового упругого элемента протеза стопы с зазором между листами при сильном изгибе // 11-я Всероссийская конференция молодых ученых “Математическое моделирование в естественных науках”. Пермь, 2002. C. 98-99.

6. Брынских С.И. Моделирование и оптимизация двухлистовой рессоры // Областная научная конференция молодых ученых, студентов и аспирантов “Молодежная наука Прикамья - 2002”. Пермь, 2002. C. 34.

7. Osipenko M.A., Nyashin Y.I., Rudakov R.N., Ostanin A.V., Kuleshova E.N., Zhuravleva T.N. Mathematical modelling of the foot prosthesis elastic element under bending // Russian Journal of Biomechanics. 2001. V. 5. № 2. P. 18-29.

8. Няшин Ю.И., Осипенко М.А., Рудаков Р.Н. К теории изгиба листовой рессоры // Механика твердого тела. 2002. № 6. C. 134-143.

9. Смольский Ю.И., Яцынин Н.А., Толмачев И.А., Суханов Д.В., Копыл Н.И., Песошников Е.М. Композиционные конструкции в протезно-ортопедических изделиях // XV Российская школа по проблемам проектирования неоднородных конструкций. Миасс, 1996. C. 14-16.

10. РаботновЮ.Н. Механика деформируемого твердого тела. М.: Наука, 1988.

11. КалиткинН.Н. Численные методы. М.: Наука, 1978.

12. Банди Б. Методы оптимизации. М.: Радио и связь, 1988.

THE OPTIMIZATION OF THE TWO-LEAF ELASTIC ELEMENT OF THE FOOT PROSTHESIS USING LINEAR AND NON-LINEAR BENDING THEORIES

S.I. Brynskikh, M.A. Osipenko, Y.I. Nyashin (Perm, Russia)

The mathematical models of linear and non-linear bending of the elastic element of the foot prosthesis that is a two-leaf spring with the gap between the leaves are considered. On the basis of these models the construction of the elastic element was optimized: thicknesses of the leaves and length of the upper leaf were found that enabled to minimize the stress of the elastic element with given size and flexibility. Comparison of the models was carried out.

Keywords: foot prosthesis, elastic element, two-leaf spring, linear bending, non-linear bending, optimization.

Получено 30 апреля 2003