untitledzuixinuvjieyoulianggang.m

用来计算水流经过台阶的计算流体力学有限元程序

re=0.001;n=5;l=1/n;s=5;t1=0;
N=n*s*(3*n+2)+2*n+1+(2*n+1)*n*t1;
for i=1:N
    pw(i)=1;b1(i)=0;bw1(i)=0;t(i)=0;x(i)=0;yy(i)=0;u(i)=0;v(i)=0;tt(i)=0;
    for j=1:N
        K1(i,j)=0;W1(i,j)=0;
    end
end
q=(n+1)*(n*s+1);
min=1,min1=1,
ll=0;
while (min>0.0000001&min1>0.0000001)
     ll=ll+1
    if ll>=300
       break;
    end
    if min>=200
       break;
    end
    if min1>=200
       break;
    end
    for i=(2*(n+1)):(n+1):q    %B1
        W1(i,i)=1;
        bw1(i)=-(-3*u(i)+4*u(i-1)-u(i-2))/(2*l);
    end
    for i=(q+1):(q+n-1)        %B2
        W1(i,i)=1;
        bw1(i)=(-3*v(i)+4*v(i+2*n+1)-v(i+2*(2*n+1)))/(2*l);
    end
    for i=(q+n):(2*n+1):N     %B3
        W1(i,i)=1;
        bw1(i)=0;
    end
    for i=(N-2*n+1):N-1               %B4
        W1(i,i)=1;W1(i,i-2*n-1)=-1;
        bw1(i)=0;
    end
    for i=(n+2):(n+1):(q-n)      %B5
        W1(i,i)=1;bw1(i)=(-3*u(i)+4*u(i+1)-u(i+2))/(2*l);
    end
    for i=(q+n+1):(2*n+1):(N-2*n)
        W1(i,i)=1;bw1(i)=(-3*u(i)+4*u(i+1)-u(i+2))/(2*l);
    end
    for i=1:n+1    %B6
        W1(i,i)=1;y=(n+1-i)*l;
        bw1(i)=-4+8*y; 
    end
    for i=(n+3):(2*n+1)      %inside left
        for j=i:(n+1):(i+(n+1)*(n*s-2))
            r=l*u(j)/(re*2);z=l*v(j)/(2*re);
            W1(j,j)=4;W1(j,j-1)=-(1-z);W1(j,j+1)=-(1+z);W1(j,j+n+1)=-(1-r);W1(j,j-n-1)=-(1+r);
            bw1(j)=0;
        end
    end
    for j=q-n+1:q-1
        r=l*u(j)/(re*2);z=l*v(j)/(2*re);
        W1(j,j)=4;W1(j,j-1)=-(1-z);W1(j,j+1)=-(1+z);W1(j,j+2*n+1)=-(1-r);W1(j,j-n-1)=-(1+r);
        bw1(j)=0;   
    end    
    for i=(q+n+2):(q+3*n)
        for j=i:(2*n+1):(i+(2*n+1)*(n*s-2))
            r=l*u(j)/(re*2);z=l*v(j)/(2*re);
            W1(j,j)=4;W1(j,j-1)=-(1-z);W1(j,j+1)=-(1+z);W1(j,j+2*n+1)=-(1-r);W1(j,j-2*n-1)=-(1+r);
            bw1(j)=0;
        end
    end
    pw=W1\bw1';
    for i=(n+1):(n+1):q    %B1
        K1(i,i)=1;
        b1(i)=0;
    end
    for i=(q+1):(q+n)       %B2
        K1(i,i)=1;
        b1(i)=0;
    end
    for i=(q+3*n+1):(2*n+1):N     %B3
        K1(i,i)=1;
        b1(i)=0;
    end
    for i=(N-2*n+1):(N-1)          %B4
        K1(i,i)=1;K1(i,i-2*n-1)=-1;
        b1(i)=0;
    end
    for i=(n+2):(n+1):(q-n)         %B5
        K1(i,i)=1;b1(i)=2/3;
    end
    for i=(q+n+1):(2*n+1):(N-2*n)
        K1(i,i)=1;b1(i)=2/3;
    end
    for i=1:n                       %B6
        K1(i,i)=1;y=(n+1-i)*l;
        b1(i)=2*y*y-4*y*y*y/3;
    end
    for i=(n+3):(2*n+1)             %inside left
        for j=i:(n+1):(i+(n+1)*(n*s-2))
            K1(j,j)=4;K1(j,j-1)=-1;K1(j,j+1)=-1;K1(j,j-n-1)=-1;K1(j,j+n+1)=-1;
            b1(j)=(l*l)*pw(j);
        end
    end
    for j=q-n+1:q-1                  %inside middle
        K1(j,j)=4;K1(j,j-1)=-1;K1(j,j+1)=-1;K1(j,j+2*n+1)=-1;K1(j,j-n-1)=-1;
        b1(j)=(l*l)*pw(j); 
    end    
    for i=(q+n+2):(q+3*n)            %inside left
        for j=i:(2*n+1):(i+(2*n+1)*(n*s-2))
            K1(j,j)=4;K1(j,j-1)=-1;K1(j,j+1)=-1;K1(j,j-2*n-1)=-1;K1(j,j+2*n+1)=-1;
            b1(j)=(l*l)*pw(j);
        end
    end
    p=K1\b1';
    for i=(n+1):(n+1):q    %B1
        v(i)=0;u(i)=0;
    end
    for i=(q+1):(q+n)
        u(i)=0;v(i)=0;     %B2
    end
    for i=(q+3*n+1):(2*n+1):N    %B3
        u(i)=0;v(i)=0;
    end
    for i=(N-2*n+1):(N-1)       %B4
        v(i)=0;u(i)=(p(i-1)-p(i+1))/(2*l);
    end
    for i=(n+2):(n+1):(q-n)     %B5
        u(i)=0;v(i)=0;
    end
    for i=(q+n+1):(2*n+1):(N-2*n)
        u(i)=0;v(i)=0;
    end
    for i=1:n     %B6
        y=(n+1-i)*l;
        u(i)=4*y*(1-y);v(i)=0;
    end
    for i=(n+3):(2*n+1)       %inside left
        for j=i:(n+1):(i+(n+1)*(n*s-2))
            u(j)=(p(j-1)-p(j+1))/(2*l);v(j)=-(p(j+n+1)-p(j-n-1))/(2*l);
        end
    end
    for j=q-n+1:q-1            %inside middle
        u(j)=(p(j-1)-p(j+1))/(2*l);v(j)=-(p(j+2*n+1)-p(j-n-1))/(2*l);
    end    
    for i=(q+n+2):(q+3*n)      %inside left
        for j=i:(2*n+1):(i+(2*n+1)*(n*s-2))
            u(j)=(p(j-1)-p(j+1))/(2*l);v(j)=-(p(j+2*n+1)-p(j-2*n-1))/(2*l);
        end
    end
    min=abs(pw(1)-t(1)); min1=abs(p(1)-tt(1));
    for i=2:N
        if min<=(abs(pw(i)-t(i)))
           min=(abs(pw(i)-t(i)));
        end
    end
    for i=2:N
        if min1<=(abs(p(i)-tt(i)))
           min1=(abs(p(i)-tt(i)));
        end
    end
    min,min1
    for i=1:N
        t(i)=pw(i);
    end
    for i=1:N
        tt(i)=p(i);
    end
end
k=0;
for i=1:n+1:((n*s)*(n+1)+1)
    k=k+1;
    for j=i:i+n
        x(j)=l*(k-1);
    end
end
for i=q:q+n
    x(i)=x(q);
end
for i=(q+n+1):(2*n+1):(N-2*n)
    k=k+1;
    for j=i:i+2*n
        x(j)=l*(k-1);
    end
end
k=0;
for i=1:n+1
    k=k+1;
    for j=i:(n+1):i+(n*s-1)*(n+1)
        yy(j)=2-l*(k-1);
    end
end
k=0;
for i=q-n:q+n
    k=k+1;
    for j=i:(2*n+1):i+(n*(s+t1))*(2*n+1)
        yy(j)=2-l*(k-1);
    end
end
pw';
y1=1;x1=0:0.01:s;y2=2;x2=0:0.01:2*s;x3=s;y3=0:0.01:1;x4=s:0.01:2*s;y4=0;x5=2*s;y5=0:0.01:2;x6=0;y6=1:0.01:2;
plot(x1,y1,'b',x2,y2,x3,y3,'b',x4,y4,x5,y5,x6,y6)
hold on
z=[p];
xlin=linspace(0,(2*s+t1),200);
ylin=linspace(0,2,200);
[XX,YY]=meshgrid(xlin,ylin);
ZZ=griddata(x,yy,z,XX,YY,'cubic');
[c,h]=contour(XX,YY,ZZ,20);
clabel(c,h);