five.m

静态场有限元法的数值分析与仿真

clear
%这段程序处理区域为G1．且进行节点编号、单元编号。节点和单元编号都是按
%由下到上由左到右的规则。然后利用图论绘图方式绘制G1的网格图
xmin1=input('请插入区域左端的x值');
xmax1=input('请输入区域右端的x值');
ymin1=input('请输入区域下端的y值');
ymax1=input('请输入区域上端的y值');
n1=input('请输入区域划分的纵线数');
m1=input('请输入区域划分的横线数');
a1=zeros(m1*n1);b1=zeros(m1*n1);
xy1=zeros(m1*n1,2);x1=zeros(1,m1*n1);y1=zeros(1,m1*n1);
%计算节点坐标
for k=1:n1
    for i=1:m1
        x1(i+(k-1)*m1)=xmin1+(k-1)*(xmax1-xmin1)/(n1-1);
        y1(i+(k-1)*m1)=ymin1+(i-1)*(ymax1-ymin1)/(m1-1);
    end
end
clear k i
a1(m1,m1)=2;
a1((n1-1)*m1+1,(n1-1)*m1+1)=2;
a1(1,1)=3;
a1(n1*m1,n1*m1)=3;
for i=1:(n1-2)
    a1(i*m1+1,i*m1+1)=4;
    a1((i+1)*m1,(i+1)*m1)=4;
end
clear i
for j=2:(m1-1)
    a1(j,j)=4;
    a1((n1-1)*m1+j,(n1-1)*m1+j)=4;
end
clear j
for i=1:(n1-2)
    for j=2:(m1-1)
        a1(i*m1+j,i*m1+j)=6;
    end
end
clear i j
a1=sparse(a1);
clear i j
for i=1:n1*m1-1
    b1(i,i+1)=-1;
end
clear i
for i=1:(n1-1)*m1
    b1(i,i+m1)=-1;
end
clear i
for i=1:(n1-1)*m1-1
    b1(i,i+m1+1)=-1;
end
for k=1:(n1-1)
for h=1:(n1-1)
    b1(m1*k,m1*h+1)=0;
  end
end
clear k h
b1=sparse(b1);
b1=a1+b1'+b1;
clear a1
a1=sparse(b1);
clear b1
xy1=[x1',y1'];
subplot(2,2,1);
gplot(a1,xy1);%绘制网格图
e1=zeros(2,3);
u1=[1,1,1];
for k=1:(n1-2)
    for i=1:(m1-1)
        e1=[i+(k-1)*m1,i+k*m1,k*m1+1+i;...
        i+(k-1)*m1,k*m1+1+i,i+(k-1)*m1+1];
        u1=[u1;e1];
    end
end
clear k i e1
for i=1:(m1-1)
    e1=[i+(n1-2)*m1,i+(n1-1)*m1+(m1-1)/2,...
    (n1-1)*m1+1+i+(m1-1)/2;i+(n1-2)*m1,...
    (n1-1)*m1+1+i+(m1-1)/2,i+(n1-2)*m1+1];
    u1=[u1;e1];
end
u1(1,:)=[];
clear k i e1
%这段程序处理区域G2。功能同上。
xmin2=input('请输入区域左端的x值');
xmax2=input('请输入区域右端的x值');
ymin2=input('请输入区域下端的y值');
ymax2=input('请输入区域上端的y值');
n2=input('请输入区域划分的纵向线数');
m2=input('请输入区域划分的横向线数');
a2=zeros(m2*n2);b2=zeros(m2*n2);
xy2=zeros(m2*n2,2);x2=zeros(1,m2*n2);y2=zeros(1,m2*n2);
%calculation node axit
for k=1:n2
    for i=1:m2
        x2(i+(k-1)*m2)=xmin2+(k-1)*(xmax2-xmin2)/(n2-1);
        y2(i+(k-1)*m2)=ymin2+(i-1)*(ymax2-ymin2)/(m2-1);
    end
end
clear k i
a2(m2,m2)=2;a2((n2-1)*m2+1,(n2-1)*m2+1)=2;a2(1,1)=3;
a1(n2*m2,n2*m2)=3;
for i=1:(n2-2)
    a2(i*m2+1,i*m2+1)=4;
    a2((i+1)*m2,(i+1)*m2)=4;
end
clear i 
for j=2:(m2-1)
    a2(j,j)=4;
    a2((n2-1)*m2+j,(n2-1)*m2+j)=4;
end
clear j
for i=1:(n2-2)
    for j=2:(m2-1)
        a2(i*m2+j,i*m2+j)=6;
    end
end
clear i j
a2=sparse(a2);
clear i j
for i=1:n2*m2-1
    b2(i,i+1)=-1;
end
clear i
for i=1:(n2-1)*m2
    b2(i,i+m2)=-1;
end
clear i
for i=1:(n2-1)*m2-1
    b2(i,i+m2+1)=-1;
end
for k=1:(n2-1)
    for h=1:(n2-1)
        b2(m2*k,m2*h+1)=0;
    end
end
clear k h
b2=sparse(b2);
b2=a2+b2'+b2;
clear a2
a2=sparse(b2);
clear b2
xy2=[x2',y2'];
hold on;
subplot(2,2,1);
gplot(a2,xy2);
e2=zeros(2,3);
u2=[0,0,0];
for k=1:(n2-1)
    for i=1:(m2-1)
        e2=[i+(k-1)*m2+m1*(n1-1),i+k*m2+m1*(n1-1),...
        k*m2+1+i+m1*(n1-1);i+(k-1)*m2+m1*(n1-1),...
        k*m2+1+i+m1*(n1-1),i+(k-1)*m2+1+m1*(n1-1)];
        u2=[u2;e2];
    end
end
u2(1,:)=[];
clear k i e2
x=[x1(1,1:m1*(n1-1)),x2(1,1:m2*n2)];
y=[y1(1,1:m1*(n1-1)),y2(1,1:m2*n2)];
u=[u1;u2];
%计算刚度矩阵[K]、节点数和节点坐标
N=m2*(n1+n2-1)-(n1-1)*(m1-1)/2;
K=zeros(N);
%开始计算矩阵[K]
for h=1:2*(m1-1)*(n1-1)+2*(m2-1)*(n2-1)
    mi=1;mj=2;mk=3;
    i=u(h,1);j=u(h,2);k=u(h,3);
    sm=[1,x(i),y(i);1,x(j),y(j);1,x(k),y(k)];
    s(h)=det(sm)/2;
    %计算Kij
    Kii=((y(j)-y(k))^2+(x(k)-x(j))^2)/(4*s(h));
    K(i,i)=K(i,i)+Kii;
    Kjj=((y(k)-y(i))^2+(x(k)-x(i))^2)/(4*s(h));
    K(j,j)=K(j,j)+Kjj;
    Kkk=((y(j)-y(i))^2+(x(j)-x(i))^2)/(4*s(h));
    K(k,k)=K(k,k)+Kkk;
    Kij=((y(j)-y(k))*(y(k)-y(i))+(x(k)-x(j))*(x(i)-x(k)))/(4*s(h));
    K(i,j)=K(i,j)+Kij;
    K(j,i)=K(j,i)+Kij;
    Kjk=((y(k)-y(i))*(y(i)-y(j))+(x(i)-x(k))*(x(j)-x(i)))/(4*s(h));
    K(j,k)=K(j,k)+Kjk;
    K(k,j)=K(k,j)+Kjk;
    Kki=((y(i)-y(j))*(y(j)-y(k))+(x(j)-x(i))*(x(k)-x(j)))/(4*s(h));
    K(k,i)=K(k,i)+Kki;
    K(i,k)=K(i,k)+Kki;
end
v=sparse(K);
clear i j k h mi mj mk sm Kii Kjj Kkk Kij Kjk Kki K
%解线性方程组
for i=0:m2
    v(N-i,:)=[];
end
clear i
for i=0:(n2-3)
    v(N-(2+i)*m2,:)=[];
end
clear i
for i=1:(m1+1)/2
    v(N-(n2-1)*m2-(m1-2)-i,:)=[];
end
clear i
for i=0:(n1-2)
    v(N-n2*m2-i*m1,:)=[];
    v(N-n2*m2+1-(1+i)*m1,:)=[];
end
clear i
q=v(:,1);
for i=0:(n1-2)
    q1=v(:,1+m1*i);q2=v(:,m1*(1+i));
    q=[q,q1,q2];
end
q(:,1)=[];
clear i
for i=1:(m1+1)/2
    q1=v(:,m1*(n1-1)+i);q=[q,q1];
end
clear i
for i=0:(n2-2)
    q1=v(:,m1*(n1-1)+m2*(i+1));q=[q,q1];
end
clear i
for i=(N-m2+1):N
    q1=v(:,i);q=[q,q1];
end
for i=0:m2
    v(:,N-i)=[];
end
clear i
for i=0:(n2-3)
    v(:,N-(2+i)*m2)=[];
end
clear i
for i=1:(m1+1)/2
    v(:,N-(n2-1)*m2-(m1-2)-i)=[];
end
clear i
for i=0:(n1-2)
    v(:,N-n2*m2-i*m1)=[];
    v(:,N-n2*m2+1-(1+i)*m1)=[];
end
clear i
M1=N-2*n1-n2-m2-(m1+1)/2+3;
f=ones(M1,1);f=4*s(1)*f;
for i=1:(m1-2)
    f(i,1)=2*s(1);
end
clear i
for i=0:(n2-3)
    f((m1-2)*n1+i*(m2-1)+1,1)=2*s(1);
end
clear i
l=1;
for i=1:(n1-1)
    l=[l,0];l=[l,10];
end
l(1)=[];
clear i
for i=1:(m1+1)/2
    l=[l,0];
end
clear i
for i=1:(n2-1)+m2
    l=[l,10];
end
clear i
g=f-q*l';r=v\g;
%后处理工作
w3=zeros(m1-2,1);
for i=1:n1
    w1(i)=0;w2(i)=10;
    w3=[w3,r((i-1)*(m1-2)+1:i*(m1-2),1)];
end
w3(:,1)=[];
w=[w1;w3;w2];
clear i
w4=0;
for i=1:(m1+1)/2-1
    w4=[w4;0];
end
clear i
w4=[w4;r((m1-2)*(n1-1)+1:(m1-2)*n1,1)];
for i=0:(n2-3)
    w4=[w4,r((m1-2)*n1+(m2-1)*i+1:(m1-2)*n1+(m2-1)*(i+1),1)];
end
clear i
w5=10*ones((m2-1),1);
w6=[w4,w5];w7=10*ones(1,n2);w6=[w6;w7];
t1=xmin1:(xmax1-xmin1)/(n1-1):xmax1;
T1=ymin1:(ymax1-ymin1)/(m1-1):ymax1;
t2=xmin2:(xmax2-xmin2)/(n2-1):xmax2;
T2=ymin2:(ymax2-ymin2)/(m2-1):ymax2;
subplot(2,2,2);meshc(t1,T1,w);axis([0 2 0 3 0 10]);hold on;
meshc(t2,T2,w6);
[t1,T1]=meshgrid(t1,T1);
[t2,T2]=meshgrid(t2,T2);
subplot(2,2,3);surfc(t1,T1,w);axis([0 2 0 3 0 10]);
hold on;surfc(t2,T2,w6);view(60,20);axis([0 2 0 3 0 10]);
subplot(2,2,4);contour(t1,T1,w);axis([0 2 0 3]);
hold on;contour(t2,T2,w6);axis([0 2 0 3]);
view(0,90);
w1=flipud(w),w2=flipud(w6)