книги / Пакеты прикладных программ

Рис. 5. Плоское тело

Требуется определить напряженно-деформированное состояние данной области.

Алгоритм метода конечных элементов для плоского сим- плекс-элемента реализуется в программной среде MATLAB.

1. Матрица жесткости.

function [K,N]=DBK(x,y,tri);

nu=0.3; t=1; E=2*1e6; D(1,1)=1; D(1,2)=nu; D(1,3)=0; D(2,1)=nu; D(2,2)=1; D(2,3)=0;

D(3,1)=0; D(3,2)=0; D(3,3)=(1-nu)/2; D=D*E/(1-nu^2);

N=inv(D)

K=zeros(2*length(x)); for i=1:length(tri)

l1=sqrt( (x(tri(i,2))-x(tri(i,1)))^2+(y(tri(i,2))- y(tri(i,1)))^2 );

l2=sqrt( (x(tri(i,3))-x(tri(i,1)))^2+(y(tri(i,3))- y(tri(i,1)))^2 );

l3=sqrt( (x(tri(i,3))-x(tri(i,2)))^2+(y(tri(i,3))- y(tri(i,2)))^2 );

p=(l1+l2+l3)/2;

S=sqrt(p*(p-l1)*(p-l2)*(p-l3)); B=zeros(3,6);

k1=1; k2=2; for j=1:3

k1=k1+1; k2=k2+1; if k1==4, k1=1; end

21

if k2==4, k2=1; end

B(1,2*j-1)=(y(tri(i,k1))-y(tri(i,k2)))/2/S; B(2,2*j)=(x(tri(i,k2))-x(tri(i,k1)))/2/S; B(3,2*j-1)=B(2,2*j);

B(3,2*j)=B(1,2*j-1);

end Ke=transpose(B)*D*B*S*t; for j=1:3

n(2*j-1)=2*tri(i,j)-1; n(2*j)=2*tri(i,j);

end

for j=1:6

for jj=1:6 K(n(j),n(jj))=K(n(j),n(jj))+Ke(j,jj);

end

end end

2. Граничные условия.

function

[K,F]=forces_and_gu(x,y,xz,yz,fzx,fzy,xp,yp,Px,Py,K,F); eps=1e-6;

for j=1:length(xz) for i=1:length(x)

if (abs(x(i)-xz(j))<eps)&(abs(y(i)-yz(j))<eps) K(2*i-1,2*i-1)=1e20;

K(2*i,2*i)=1e20; end

end end

for j=1:length(xp) for i=1:length(x)

if (x(i)==xp(j))&(y(i)==yp(j)) F(2*i-1,1)=Px(j); F(2*i,1)=Py(j);

end

end

end

3. Деформации сетки.

function [x,y]=DeformSetka(x,y,tri,u); f=figure; Ms=0.1*(max(x)-min(x))/(max(abs(u)));

22

for i=1:length(x) x(i)=x(i)+u(2*i-1)*Ms; y(i)=y(i)+u(2*i)*Ms;

end

plot(x,y,'.','color',[1,0,0]); for i=1:length(tri)

xx(1)=x(tri(i,1));

xx(2)=x(tri(i,2));

xx(3)=x(tri(i,3));

xx(4)=x(tri(i,1));

yy(1)=y(tri(i,1));

yy(2)=y(tri(i,2));

yy(3)=y(tri(i,3));

yy(4)=y(tri(i,1)); plot(xx,yy,'color',[1,0,0]); hold on

end

4. Поверхностные силы.

function F=V_S_Forces(r1,r2,x,y,tri,F,gx,gy,PE,PI,fi1,fi2,IE); t=1;

for i=1:length(tri)

L1=sqrt((x(tri(i,2))-x(tri(i,1)))^2+(y(tri(i,2))- y(tri(i,1)))^2);

L2=sqrt((x(tri(i,3))-x(tri(i,2)))^2+(y(tri(i,3))- y(tri(i,2)))^2);

ab=((x(tri(i,3))-x(tri(i,2)))*(x(tri(i,2))- x(tri(i,1))))+((y(tri(i,3))-y(tri(i,2)))*(y(tri(i,2))- y(tri(i,1))));

cosalpha=(ab)/(L1*L2); sinalpha=sqrt(1-cosalpha^2); S=L1*L2*sinalpha/2;

for j=1:3 F(2*tri(i,j),1)=F(2*tri(i,j),1)+1/3*t*S*gy;

F(2*tri(i,j)-1,1)=F(2*tri(i,j)-1,1)+1/3*t*S*gx; end

end

eps=1e-10; j=1; jj=1; for i=1:length(tri(:,1))

usl1=sqrt(x(tri(i,1))^2+y(tri(i,1))^2)+eps;

usl2=sqrt(x(tri(i,2))^2+y(tri(i,2))^2)+eps;

usl3=sqrt(x(tri(i,3))^2+y(tri(i,3))^2)+eps; if (usl1>r2)&(usl2>r2)

23

ES(j,1)=i;

ES(j,2)=tri(i,1);

ES(j,3)=tri(i,2); j=j+1; end;

if (usl1>r2)&(usl3>r2) ES(j,1)=i; ES(j,2)=tri(i,1); ES(j,3)=tri(i,3); j=j+1; end;

if (usl2>r2)&(usl3>r2) ES(j,1)=i; ES(j,2)=tri(i,2); ES(j,3)=tri(i,3); j=j+1; end;

usl1=sqrt(x(tri(i,1))^2+y(tri(i,1))^2)-eps; usl2=sqrt(x(tri(i,2))^2+y(tri(i,2))^2)-eps; usl3=sqrt(x(tri(i,3))^2+y(tri(i,3))^2)-eps;

if (usl1<r1)&(usl2<r1) EI(jj,1)=i; EI(jj,2)=tri(i,1); EI(jj,3)=tri(i,2); jj=jj+1; end;

if (usl1<r1)&(usl3<r1) EI(jj,1)=i; EI(jj,2)=tri(i,1); EI(jj,3)=tri(i,3); jj=jj+1; end;

if (usl2<r1)&(usl3<r1) EI(jj,1)=i; EI(jj,2)=tri(i,2); EI(jj,3)=tri(i,3); jj=jj+1; end;

end l_1=0;

if (IE==1)|(IE==3)

for i=1:length(ES(:,1)) l=sqrt((x(ES(i,3))-x(ES(i,2)))^2+(y(ES(i,3))-

y(ES(i,2)))^2); l_1=l_1+l;

fi=atan(y(ES(i,2))/x(ES(i,2))); if x(ES(i,2))>=0

F(ES(i,2)*2-1)=F(ES(i,2)*2-1)+PE*cos(fi)*t*l/2; F(ES(i,2)*2)=F(ES(i,2)*2)+PE*sin(fi)*t*l/2;

24

else F(ES(i,2)*2-1)=F(ES(i,2)*2-1)-PE*cos(fi)*t*l/2; F(ES(i,2)*2)=F(ES(i,2)*2)-PE*sin(fi)*t*l/2; end;

fi=atan(y(ES(i,3))/x(ES(i,3))); if x(ES(i,3))>=0

F(ES(i,3)*2-1)=F(ES(i,3)*2-1)+PE*cos(fi)*t*l/2; F(ES(i,3)*2)=F(ES(i,3)*2)+PE*sin(fi)*t*l/2;

else F(ES(i,3)*2-1)=F(ES(i,3)*2-1)-PE*cos(fi)*t*l/2; F(ES(i,3)*2)=F(ES(i,3)*2)-PE*sin(fi)*t*l/2; end;

end

end l_1=0;

if (IE==2)|(IE==3)

for i=1:length(EI(:,1)) l=sqrt((x(EI(i,3))-x(EI(i,2)))^2+(y(EI(i,3))-

y(EI(i,2)))^2); l_1=l_1+l;

fi=atan(y(EI(i,2))/x(EI(i,2))); if x(EI(i,2))>=0

(EI(i,2)*2-1)=F(EI(i,2)*2-1)+PI*cos(fi)*t*l/2; F(EI(i,2)*2)=F(EI(i,2)*2)+PI*sin(fi)*t*l/2;

else F(EI(i,2)*2-1)=F(EI(i,2)*2-1)-PI*cos(fi)*t*l/2; F(EI(i,2)*2)=F(EI(i,2)*2)-PI*sin(fi)*t*l/2; end;

fi=atan(y(EI(i,3))/x(EI(i,3))); if x(EI(i,3))>=0

F(EI(i,3)*2-1)=F(EI(i,3)*2-1)+PI*cos(fi)*t*l/2; F(EI(i,3)*2)=F(EI(i,3)*2)+PI*sin(fi)*t*l/2;

else F(EI(i,3)*2-1)=F(EI(i,3)*2-1)-PI*cos(fi)*t*l/2; F(EI(i,3)*2)=F(EI(i,3)*2)-PI*sin(fi)*t*l/2; end;

end

end

5. Определение напряжений и деформаций.

function

[Strain_x,Strain_y,Strain_xy,Stress_x,Stress_y,Stress_xy] =StressStrain(x,y,tri,u);

nu=0.5; t=1; E=2*1e6; D(1,1)=1; D(1,2)=nu; D(1,3)=0; D(2,1)=nu; D(2,2)=1; D(2,3)=0;

D(3,1)=0; D(3,2)=0; D(3,3)=(1-nu)/2; D=D*E/(1-nu^2); Strain_x=zeros(1,length(x));

Strain_y=Strain_x;

Strain_xy=Strain_x;

25

Stress_x=Strain_x;

Stress_y=Strain_x;

Stress_xy=Strain_x; Ssum=Strain_x;

for i=1:length(tri)

L1=sqrt((x(tri(i,2))-x(tri(i,1)))^2+(y(tri(i,2))- y(tri(i,1)))^2);

L2=sqrt((x(tri(i,3))-x(tri(i,2)))^2+(y(tri(i,3))- y(tri(i,2)))^2);

ab=((x(tri(i,3))-x(tri(i,2)))*(x(tri(i,2))- x(tri(i,1))))+((y(tri(i,3))-y(tri(i,2)))*(y(tri(i,2))- y(tri(i,1))));

cosalpha=(ab)/(L1*L2); sinalpha=sqrt(1-cosalpha^2); S=L1*L2*sinalpha/2; B=zeros(3,6);

k1=1; k2=2; for j=1:3

k1=k1+1; k2=k2+1; if k1==4, k1=1; end if k2==4, k2=1; end

B(1,2*j-1)=(y(tri(i,k1))-y(tri(i,k2)))/(2*S); B(2,2*j)=(x(tri(i,k2))-x(tri(i,k1)))/(2*S); B(3,2*j-1)=B(2,2*j);

B(3,2*j)=B(1,2*j-1);

end

for j=1:3

n(2*j-1)=2*tri(i,j)-1; n(2*j)=2*tri(i,j);

end

for j=1:6 u_e(j,1)=u(n(j)); end

Strain=B*u_e; Stress=D*Strain for j=1:3

Strain_x(tri(i,j))=Strain_x(tri(i,j))+Strain(1)*S; Strain_y(tri(i,j))=Strain_y(tri(i,j))+Strain(2)*S; Strain_xy(tri(i,j))=Strain_xy(tri(i,j))+Strain(3)*S; Stress_x(tri(i,j))=Stress_x(tri(i,j))+Stress(1)*S; Stress_y(tri(i,j))=Stress_y(tri(i,j))+Stress(2)*S; Stress_xy(tri(i,j))=Stress_xy(tri(i,j))+Stress(3)*S; Ssum(tri(i,j))=Ssum(tri(i,j))+S;

end

end

26

for i=1:length(x)
Strain_x(i)=Strain_x(i)/Ssum(i);
Strain_y(i)=Strain_y(i)/Ssum(i);
Strain_xy(i)=Strain_xy(i)/Ssum(i);
Stress_x(i)=Stress_x(i)/Ssum(i);
Stress_y(i)=Stress_y(i)/Ssum(i);
Stress_xy(i)=Stress_xy(i)/Ssum(i);
end
										-30.0125
2										-29.7047
										-29.3627
1.5										-29.0207
										-28.6787
1										-28.3367
										-27.9946
0.5										-27.6526
										-27.3106
0										-26.9686
										-26.6266
-0.5										-26.2846
										-25.9426
-1										-25.6005
										-25.2585
-1.5										-24.9165
										-24.5745
-2										-24.2325
										-23.8905
-2.5										-23.5485
-2.5	-2	-1.5	-1	-0.5	0	0.5	1	1.5	2	2.5-23.2064
										-15.0063
2										-14.5885
										-14.1243
1.5										-13.6601
										-13.1959
1										-12.7317
										-12.2676
0.5										-11.8034
										-11.3392
0										-10.875
										-10.4108
-0.5										-9.9466
										-9.4824
-1										-9.0183
										-8.5541
-1.5										-8.0899
										-7.6257
-2										-7.1615
										-6.6973
-2.5										-6.2331
-2.5	-2	-1.5	-1	-0.5	0	0.5	1	1.5	2	2.5 -5.7689
	Рис. 6. Изо линии напряжений (см. с. 28)

27

											-4.154
											-3.8362
2											-3.4831
1.5											-3.13
											-2.7769
1											-2.4239
											-2.0708
0.5											-1.7177
											-1.3646
0											-1.0115
											-0.65844
-0.5											-0.30536
-1											0.047726
											0.40081
-1.5											0.75389
											1.107
-2											1.4601
											1.8131
-2.5											2.1662
-2.5	-2	-1.5	-1	-0.5	0	0.5	1	1.5	2	2.5	2.5193
											2.8724

Рис. 6. Окончание

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

28

ГЛАВА 3. ИССЛЕДОВАНИЕ НАПРЯЖЕННО-ДЕФОРМИРОВАННОГО СОСТОЯНИЯ В ЗАДАЧАХ ПРОТЕЗИРОВАНИЯ

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

Эндопротезирование тазобедренного сустава в настоящее время является широко распространенным способом лечения заболеваний опорно-двигательного аппарата. Замена больного сустава на искусственный позволяет устранить или значительно уменьшить болевой синдром, обеспечить опороспособность конечности, восстановить движение в суставе. Однако срок службы современных эндопротезов, изготавливаемых из легированных сталей и титана, ограничивается в среднем 10 годами. Проведение повторной операции эндопротезирования зачастую невозможно из-за большого количества противопоказаний и высокого риска развития послеоперационных осложнений. Потребность в эндопротезировании в России составляет до 100–300 тысяч операций в год, поэтому продление эксплуатационного ресурса эндопротеза является актуальной медицинской, технической и социальной проблемой.

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

29

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

Интраоперационные переломы бедренной кости представляют собой продольное растрескивание проксимального отдела бедра. Частота данного типа осложнения составляет 3,6 %, при этом перелом диафиза бедра встречается в 2,1 % случаев, большого вертела – в 1,0 % случаев, малого вертела – в 0,5 % случаев. Основные причины возникновения перелома – это допущенные ошибки на стадии предоперационного планирования (определение по ретгенограммам длины шейки бедренной кости и диаметра костномозгового канала), грубая техника проведения оперативного вмешательства при желании обеспечить более надежную начальную фиксацию бедренного компонента. Факторами риска являются предшествовавшие оперативные вмешательства (остеотомия бедренной кости, остеосинтез по поводу перелома шейки бедра и т.п.), врожденная дисплазия тазобедренного сустава и врожденный вывих бедра, узкий и искривленный костно-мозговой канал, тяжелая сгибательная или наружно-ротационная контрактура в тазобедренном суставе, остеопороз, избыточный вес больного, бесцементная фиксация имплантата, применение эндопротеза с прямой ножкой. Профилактикой данного типа осложнения являются тщательная предоперационная рентгенометрия формы костно-мозгового канала бедренной кости и соблюдение техники проведения операции.

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

30