计算流体力学显示格式_计算流体力学_流体力学作业；显示格式_计算流体_差分计算流体_源码

load jiaohe2.mat KK = zeros(num_node*2,num_node*2); element=struct( 'node1' , [] , 'node2' , [] ,'node3',[],'k' , [] , 'stress' , [],'B',[]); i=1; for j=1:num_line-1 for k=1:num_list-1 l=k+1; m=j+1; element(i).eindex=i; %单元号 element(i).node1=[j,k]; element(i).node2=[j,l]; element(i).node3=[m,l]; element(i).B=planeB(node(j,k).x,node(j,k).y,node(j,l).x,node(j,l).y,node(m,l).x,node(m,l).y,a,b); element(i).k=a*b*t*element(i).B'*D*element(i).B/2; KK=Planezz(KK,element(i).k,(j-1)*num_list+k,(j-1)*num_list+l,(m-1)*num_list+l); i=i+1; element(i).eindex=i; %单元号 element(i).node1=[j,k]; element(i).node2=[m,k]; element(i).node3=[m,l]; element(i).B=planeB(node(j,k).x,node(j,k).y,node(m,k).x,node(m,k).y,node(m,l).x,node(m,l).y,a,b); element(i).k=a*b*t*element(i).B'*D*element(i).B/2; KK=Planezz(KK,element(i).k,(j-1)*num_list+k,(m-1)*num_list+k,(m-1)*num_list+l); i=i+1; end end index=[]; %index记录已知外力的节点的编号 p=[]; for j=1:num_line for k=1:num_list if node(j,k).Fx~=1e100 index=[index,2*((j-1)*num_list+k)-1]; p=[p,node(j,k).Fx]; end if node(j,k).Fy~=1e100 index=[index,2*((j-1)*num_list+k)]; p=[p,node(j,k).Fy]; end end end u=KK(index,index)\p'; for j=1:length(u) if mod(index(j),2)==0 if mod(index(j)/2,num_list)==0 node(floor((index(j)/2)/num_list),num_list).v=u(j); else node(floor((index(j)/2)/num_list+1),mod(index(j)/2,num_list)).v=u(j); end else if mod((index(j)+1)/2,num_list)==0 node(floor(((index(j)+1)/2)/num_list),num_list).u=u(j); else node(floor(((index(j)+1)/2)/num_list+1),mod((index(j)+1)/2,num_list)).u=u(j); end end end weiyi=[]; for i=1:num_line for j=1:num_list weiyi=[weiyi,node(i,j).u,node(i,j).v]; end end force=KK*weiyi'; for i=1:num_line for j=1:num_list node(i,j).Fx=force(2*((i-1)*num_list+j)-1); node(i,j).Fy=force(2*((i-1)*num_list+j)); end end S=struct('s',[],'stress',[]);%单元体形变和应力 for i=1:num_element S(i).s=Planewy(element(i).B,node(element(i).node1(1),element(i).node1(2)).u,node(element(i).node1(1),element(i).node1(2)).v,node(element(i).node2(1),element(i).node2(2)).u,node(element(i).node2(1),element(i).node2(2)).v,node(element(i).node3(1),element(i).node3(2)).u,node(element(i).node3(1),element(i).node3(2)).v); end for i=1:num_element S(i).stress=D*S(i).s; end