matlab小程序-平面应力有限元求解器

% A=csvread('W-M fractal surface_D1_G1.5.csv'); % x=A(1,:); % y1=A(2,:); % y2=cos(x); % % figure % [AX,H1,H2]=plotyy(x,y1,x,y2,'plot'); % set(get(AX(1),'Ylabel'),'string','sinx') % set(get(AX(2),'Ylabel'),'string','cosx') % set(H1,'linestyle',':') % set(H1,'Color','r') % yyaxis left % H1=plot(x,y1); % % xlim([0,2*pi]) % yyaxis right % H2=plot(x,y2); % yyaxis left % ylabel('sinx') % set(H2,'linestyle',':') % set(H2,'Color','g') % yyaxis right % ylabel('cosx') clc clear %读取平面三角单元 %%------------- %读取单元 e='ELIST.dat'; ELE=fopen(e); element=textscan(ELE,'%s %*s %*s %*s %*s %*s %s %s %s %*s ','HeaderLines',2); fclose(ELE); id=element{1}; n1=element{2}; n2=element{3}; n3=element{4}; a=length(id); ID=zeros(a,1); N1=zeros(a,1); N2=zeros(a,1); N3=zeros(a,1); for i=1:a ID(i,1)=str2double(id(i)); N1(i,1)=str2double(n1(i)); N2(i,1)=str2double(n2(i)); N3(i,1)=str2double(n3(i)); end ELEMENT=[rmmissing(ID),rmmissing(N1),rmmissing(N2),rmmissing(N3)]; %%-------------- %读取节点 N='NLIST.lis'; no=fopen(N); nodes=textscan(no,'%s %s %s %*s','HeaderLine',3); fclose(no); nid=nodes{1}; nx=nodes{2}; ny=nodes{3}; b=length(nid); NID=zeros(b,1); NX=zeros(b,1); NY=zeros(b,1); for j=1:b NID(j,1)=str2double(nid{j}); NX(j,1)=str2double(nx{j}); NY(j,1)=str2double(ny{j}); end NODES=[rmmissing(NID),rmmissing(NX),rmmissing(NY)]; %节点单元读取结束 %%-------------------------------------------------- %计算单元刚度矩阵 %材料参数 mu=0.3; E=1e9; t=0.1; %获取单元节点号 NID1=ELEMENT(:,2); NID2=ELEMENT(:,3); NID3=ELEMENT(:,4); %获取节点坐标并计算单元刚度矩阵 for i=1:length(NID1) index1=NODES(:,1) == NID1(i); index2=NODES(:,1) == NID2(i); index3=NODES(:,1) == NID3(i); NODE1=NODES(index1,2:end); NODE2=NODES(index2,2:end); NODE3=NODES(index3,2:end); A=stiffness(NODE1(1),NODE1(2),NODE2(1),NODE2(2),NODE3(1),NODE3(2),mu,E,t); K{1,i}=A; end %组装整体刚度矩阵 EID=NODES(:,1); DOF=2*length(EID); KK=zeros(DOF); for i=1:length(ELEMENT) B=Assembly(K{i},NID1(i),NID2(i),NID3(i),EID); KK=KK+B; end %组装完成 %%-------------------------------------------------- %施加边界条件 %位移边界 ADOF=[1 26 46 47 48];%全约束节点 for i=1:length(ADOF) KK(2*ADOF(i),:)=0; KK(:,2*ADOF(i))=0; KK(2*ADOF(i),2*ADOF(i))=1; KK(2*ADOF(i)-1,:)=0; KK(:,2*ADOF(i)-1)=0; KK(2*ADOF(i)-1,2*ADOF(i)-1)=1; end %应力边界条件 R=zeros(DOF,1); RDOF=22; R(RDOF*2,1)=-1e3; %求解位移 c=det(KK);%检测是否奇异 if c==0 answer='刚度矩阵奇异'; else answer='计算完毕'; delta=KK\R; ST={length(ELEMENT),2}; DisE=zeros(6,1); %计算应力应变 for i=1:length(NID1) index1=NODES(:,1) == NID1(i); index2=NODES(:,1) == NID2(i); index3=NODES(:,1) == NID3(i); NODE1=NODES(index1,2:end); NODE2=NODES(index2,2:end); NODE3=NODES(index3,2:end); ID1=2*NID1(i)-1; ID2=2*NID1(i); ID3=2*NID2(i)-1; ID4=2*NID2(i); ID5=2*NID3(i)-1; ID6=2*NID3(i); DisE=[delta(ID1);delta(ID2);delta(ID3);delta(ID4);delta(ID5);delta(ID6)]; [strain,stress]=strain_compu(NODE1(1),NODE1(2),NODE2(1),NODE2(2),NODE3(1),NODE3(2),mu,E,DisE); ST{i,1}=strain; ST{i,2}=stress; end STE=cell2mat(ST); end answer dd=min(delta(2:2:end)); figure Displacement=patch('Faces',ELEMENT(:,2:4),'Vertices',NODES(:,2:3),'Edgecolor','none','facevertexCdata',delta(2:2:end),'FaceColor','interp'); ylim([-4,6]); colorbar; figure Strain=patch('Faces',ELEMENT(:,2:4),'Vertices',NODES(:,2:3),'Edgecolor','none','FaceVertexCdata',STE(2:3:end,2),'FaceColor','flat'); ylim([-4,6]); colorbar;