84 lines
2.5 KiB
Mathematica
84 lines
2.5 KiB
Mathematica
|
clear all;clc;close all;
|
|||
|
|
global E u SD q0 Nx Ny
|
|||
|
|
Nx=4;Ny=4;
|
|||
|
|
E=2e6;%<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ģ<EFBFBD><EFBFBD>
|
|||
|
|
u=0; %<EFBFBD><EFBFBD><EFBFBD>ɱ<EFBFBD>
|
|||
|
|
SD=2;%<EFBFBD><EFBFBD>ά<EFBFBD><EFBFBD><EFBFBD><EFBFBD>
|
|||
|
|
q0=-10;%ÿһ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>صĺ<EFBFBD><EFBFBD><EFBFBD>
|
|||
|
|
%<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ڵ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>%<EFBFBD><EFBFBD><EFBFBD>嵥Ԫ
|
|||
|
|
Patch=1;%<EFBFBD>ײ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
|
|||
|
|
[P1,E1,q0_node1,Boundary1] = Geometry(Patch, 0, 0);%<EFBFBD><EFBFBD><EFBFBD>ؽڵ㵥Ԫ<EFBFBD><EFBFBD>Ϣ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ʩ<EFBFBD>ӽڵ㡢<EFBFBD>߽<EFBFBD>Լ<EFBFBD><EFBFBD>ʩ<EFBFBD>ӽڵ<EFBFBD>
|
|||
|
|
dis_1=zeros(size(P1));
|
|||
|
|
Patch=2;%<EFBFBD>ϲ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
|
|||
|
|
[P2,E2,q0_node2,Boundary2] = Geometry(Patch, 0, 0);
|
|||
|
|
dis_2=zeros(size(P2));
|
|||
|
|
%<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ڵ<EFBFBD><EFBFBD>͵<EFBFBD>Ԫ
|
|||
|
|
Node=[P1;P2];
|
|||
|
|
Element=[E1;E2+size(P1,1)];
|
|||
|
|
%<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ص<EFBFBD>Ԫ<EFBFBD><EFBFBD><EFBFBD>ڵ<EFBFBD><EFBFBD><EFBFBD>
|
|||
|
|
[Nnod,Nele]=MeshPlot(Node,Element);
|
|||
|
|
%<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ϸնȾ<EFBFBD><EFBFBD><EFBFBD>
|
|||
|
|
D=[E/(1-u^2),u*E/(1-u^2),0;
|
|||
|
|
u*E/(1-u^2),E/(1-u^2),0;
|
|||
|
|
0,0,E/(2*(1+u))];
|
|||
|
|
time=1;
|
|||
|
|
d_final=zeros(2*size(Node,1),1);
|
|||
|
|
|
|||
|
|
% f=figure;
|
|||
|
|
% interval=5; %plot gif
|
|||
|
|
while(time <100)
|
|||
|
|
% clf(f);%plot gif
|
|||
|
|
%<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>²<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ն<EFBFBD><EFBFBD><EFBFBD>
|
|||
|
|
Patch=1;
|
|||
|
|
[K_Patch_1, F_Patch_1] = Linear_Elasticity(P1, E1, q0_node1,Boundary1,Patch);
|
|||
|
|
%<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ϲ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ն<EFBFBD><EFBFBD><EFBFBD>
|
|||
|
|
Patch=2;
|
|||
|
|
[K_Patch_2, F_Patch_2] = Linear_Elasticity(P2, E2,q0_node2,Boundary2, Patch);
|
|||
|
|
%<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ȫ<EFBFBD><EFBFBD>patch<EFBFBD>ĸնȾ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ͺ<EFBFBD><EFBFBD>ؾ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>װΪһ<EFBFBD><EFBFBD>
|
|||
|
|
K_Global=[K_Patch_1,zeros(size(K_Patch_1, 1),size(K_Patch_2,2)) ;zeros(size(K_Patch_2, 1),size(K_Patch_1,2)) ,K_Patch_2];
|
|||
|
|
F_Global=[F_Patch_1; F_Patch_2];
|
|||
|
|
%<EFBFBD><EFBFBD><EFBFBD>ǽӴ<EFBFBD>Լ<EFBFBD><EFBFBD>
|
|||
|
|
Penalty=1e7;%<EFBFBD>ͷ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
|
|||
|
|
[K_Contact RC_force] = Contact_Formualtion(P1, P2,E1,E2, Penalty, dis_1, dis_2);
|
|||
|
|
K_Global=K_Global+K_Contact;
|
|||
|
|
d=linsolve(K_Global, F_Global);
|
|||
|
|
d1=d(1:SD*size(P1,1));
|
|||
|
|
d2=d(SD*size(P1,1)+1:SD*(size(P1,1)+size(P2,1)));
|
|||
|
|
dis_1=dis_1+transpose(reshape(d1, size(dis_1,2), size(dis_1,1)));
|
|||
|
|
dis_2=dis_2+transpose(reshape(d2, size(dis_2,2), size(dis_2,1)));
|
|||
|
|
d_final=d_final+d;
|
|||
|
|
|
|||
|
|
time=time+1;
|
|||
|
|
|
|||
|
|
|
|||
|
|
|
|||
|
|
% %plot gif <EFBFBD><EFBFBD><EFBFBD><EFBFBD>GIF<EFBFBD><EFBFBD>ͼ
|
|||
|
|
% [NodeStrain,NodeStress]=CalculateStrainAndStress2(d_final,D,Node,Element);
|
|||
|
|
% PlotContour(Node,Element,d_final,NodeStress(2,:))%sigma_11
|
|||
|
|
% axis([0 1 -0.1 0.7])
|
|||
|
|
% caxis([-1591 123])
|
|||
|
|
% colorbar off
|
|||
|
|
% drawnow;
|
|||
|
|
% F=getframe(f);
|
|||
|
|
% I=frame2im(F);
|
|||
|
|
% [I,map]=rgb2ind(I,256);
|
|||
|
|
% if time == 2
|
|||
|
|
% imwrite(I,map,'sigma_22.gif','gif', 'Loopcount',inf,'DelayTime',0);
|
|||
|
|
% elseif mod(time,interval)==0
|
|||
|
|
% imwrite(I,map,'sigma_22.gif','gif','WriteMode','append','DelayTime',0);
|
|||
|
|
% end
|
|||
|
|
|
|||
|
|
end
|
|||
|
|
% <EFBFBD>洢<EFBFBD><EFBFBD><EFBFBD><EFBFBD>λ<EFBFBD><EFBFBD>
|
|||
|
|
gNTu = [dis_1;dis_2];
|
|||
|
|
|
|||
|
|
% <EFBFBD><EFBFBD><EFBFBD>Ʊ<EFBFBD><EFBFBD><EFBFBD>ͼ
|
|||
|
|
SF=1;
|
|||
|
|
PlotResults2(Node,Element,gNTu,100);
|
|||
|
|
[NodeStrain,NodeStress]=CalculateStrainAndStress2(d_final,D,Node,Element);
|
|||
|
|
figure
|
|||
|
|
PlotContour(Node,Element,d_final,NodeStress(1,:))%sigma_11
|
|||
|
|
title('\sigma11');
|
|||
|
|
figure
|
|||
|
|
PlotContour(Node,Element,d_final,NodeStress(2,:))%sigma_22
|
|||
|
|
title('\sigma22');
|