84 lines
2.5 KiB
Matlab
84 lines
2.5 KiB
Matlab
clear all;clc;close all;
|
|
global E u SD q0 Nx Ny
|
|
Nx=4;Ny=4;
|
|
E=2e6;%杨氏模量
|
|
u=0; %泊松比
|
|
SD=2;%二维问题
|
|
q0=-10;%每一步加载的荷载
|
|
%定义节点坐标%定义单元
|
|
Patch=1;%底部块体
|
|
[P1,E1,q0_node1,Boundary1] = Geometry(Patch, 0, 0);%返回节点单元信息、荷载施加节点、边界约束施加节点
|
|
dis_1=zeros(size(P1));
|
|
Patch=2;%上部块体
|
|
[P2,E2,q0_node2,Boundary2] = Geometry(Patch, 0, 0);
|
|
dis_2=zeros(size(P2));
|
|
%定义节点和单元
|
|
Node=[P1;P2];
|
|
Element=[E1;E2+size(P1,1)];
|
|
%绘制网格返回单元数节点数
|
|
[Nnod,Nele]=MeshPlot(Node,Element);
|
|
%定义材料刚度矩阵
|
|
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
|
|
%建立下部块体刚度阵
|
|
Patch=1;
|
|
[K_Patch_1, F_Patch_1] = Linear_Elasticity(P1, E1, q0_node1,Boundary1,Patch);
|
|
%建立上部块体刚度阵
|
|
Patch=2;
|
|
[K_Patch_2, F_Patch_2] = Linear_Elasticity(P2, E2,q0_node2,Boundary2, Patch);
|
|
%将两个全局patch的刚度矩阵和荷载矩阵组装为一体
|
|
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];
|
|
%考虑接触约束
|
|
Penalty=1e7;%惩罚因子
|
|
[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 绘制GIF动图
|
|
% [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
|
|
% 存储结点位移
|
|
gNTu = [dis_1;dis_2];
|
|
|
|
% 绘制变形图
|
|
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'); |