查看: 62|回复: 2|关注: 0

[已答复] while循环只进行了一次,就完成了计算,请大家帮我分析一下原因

[复制链接]

新手

5 麦片

财富积分


050


1

主题

2

帖子

0

最佳答案
发表于 7 天前 | 显示全部楼层 |阅读模式
% Finite Volume Method
clear;clc;close all;
%M = input('Number of grids in x-direction = ')
%N = input('Number of grids in y-direction = ')
M=20;N=30;
dx = 0.01;
dy = 0.01;
p=zeros(M+1,N+1);
u=zeros(M+1,N+2);
v=zeros(M+2,N+1);
pn=zeros(M+1,N+1);
un=zeros(M+1,N+2);
vn=zeros(M+2,N+1);
U=zeros(M+1,N+2);
V=zeros(M+2,N+1);
P=zeros(M+1,N+1);
density = 997.05;% density
miu=0.8934;% Viscosity
u(12:20,2:31)=0.15;
u(11:21,1)=0.4-u(11:21,2);
u(2:11,12:21)=0.15;
u(11,1:11)=0;
u(11,22:32)=0;
u(21,1:32)=0;
u(11:21,32)=u(11:21,31);
u(1:11,11)=-u(1:11,12);
u(1:11,22)=u(1:11,21);
u(1,11:22)=0;
v(11:22,1)=0;
v(11:22,31)=0;
v(22,1:31)=-v(21,1:31);
v(11,1:10)=-v(12,1:10);
v(11,22:31)=-v(12,22:31);
v(1:11,11)=0;
v(1:11,21)=0;
v(1,11:21)=-v(2,11:21);
error=1;
k=0;
while error>1e-8

    for i=12:20
        for j=2:31
            awu(i,j)=miu+0.25*(u(i,j-1)+u(i,j))*density*dy;
            aeu(i,j)=miu-0.25*(u(i,j+1)+u(i,j))*density*dy;
            asu(i,j)=miu+0.25*(v(i,j)+v(i,j-1))*density*dx;       
            anu(i,j)=miu-0.25*(v(i+1,j)+v(i+1,j-1))*density*dx;
            au(i,j)=4*miu+0.25*(u(i,j+1)-u(i,j-1))*density*dy+0.25*(v(i+1,j-1)+v(i+1,j)-v(i,j)-v(i,j-1))*density*dx;
        end;
    end;
    for i=2:11
        for j=12:21
            awu(i,j)=miu+0.25*(u(i,j-1)+u(i,j))*density*dy;
            aeu(i,j)=miu-0.25*(u(i,j+1)+u(i,j))*density*dy;
            asu(i,j)=miu+0.25*(v(i,j)+v(i,j-1))*density*dx;       
            anu(i,j)=miu-0.25*(v(i+1,j)+v(i+1,j-1))*density*dx;
            au(i,j)=4*miu+0.25*(u(i,j+1)-u(i,j-1))*density*dy+0.25*(v(i+1,j-1)+v(i+1,j)-v(i,j)-v(i,j-1))*density*dx;
        end;
    end;
    for i=12:21
        for j=2:30
            awv(i,j)=miu+0.25*(u(i,j)+u(i-1,j))*density*dy;
            aev(i,j)=miu-0.25*(u(i-1,j+1)+u(i,j+1))*density*dy;
            asv(i,j)=miu+0.25*(v(i-1,j)+v(i,j))*density*dx;       
            anv(i,j)=miu-0.25*(v(i+1,j)+v(i,j))*density*dx;
            av(i,j)=4*miu+0.25*(u(i-1,j+1)+u(i,j+1)-u(i,j)-u(i-1,j))*density*dy+0.25*(v(i+1,j)-v(i-1,j))*density*dx;
       end;
   end;
   for i=2:11
        for j=12:20
            awv(i,j)=miu+0.25*(u(i,j)+u(i-1,j))*density*dy;
            aev(i,j)=miu-0.25*(u(i-1,j+1)+u(i,j+1))*density*dy;
            asv(i,j)=miu+0.25*(v(i-1,j)+v(i,j))*density*dx;       
            anv(i,j)=miu-0.25*(v(i+1,j)+v(i,j))*density*dx;
            av(i,j)=4*miu+0.25*(u(i-1,j+1)+u(i,j+1)-u(i,j)-u(i-1,j))*density*dy+0.25*(v(i+1,j)-v(i-1,j))*density*dx;
       end;
   end;
   error_u=1;
   while error_u>1e-6
        un_0(12:20,2:31)=un(12:20,2:31);      
        %%%%% caculate u %%%%%
        for i=12:20
            for j=2:31
                un(i,j)=(awu(i,j)*u(i,j-1)+aeu(i,j)*u(i,j+1)+asu(i,j)*u(i-1,j)+anu(i,j)*u(i+1,j)+dy*(p(i,j-1)-p(i,j)))/au(i,j);
            end;
        end;
        error_u=abs(un_0(12:20,2:31)-un(12:20,2:31));
   end;
   error_u=1;
   while error_u>1e-6
        un_0(2:11,12:21)=un(2:11,12:21);      
        %%%%% caculate u %%%%%
        for i=2:11
            for j=12:21
                un(i,j)=(awu(i,j)*u(i,j-1)+aeu(i,j)*u(i,j+1)+asu(i,j)*u(i-1,j)+anu(i,j)*u(i+1,j)+dy*(p(i,j-1)-p(i,j)))/au(i,j);
            end;
        end;
        error_u=abs(un_0(2:11,12:21)-un(2:11,12:21));
   end;
   error_v=1;
   while error_v>1e-6
        vn_0(12:21,2:30)=vn(12:21,2:30);
        %%%%% caculate v %%%%%
        for i=12:21
            for j=2:30
                vn(i,j)=(awv(i,j)*v(i,j-1)+aev(i,j)*v(i,j+1)+asv(i,j)*v(i-1,j)+anv(i,j)*v(i+1,j)+dx*(p(i-1,j)-p(i,j)))/av(i,j);
            end;
        end;
        error_v=abs(vn_0(12:21,2:30)-vn(12:21,2:30));
   end;
    error_v=1;
   while error_v>1e-6
        vn_0(2:11,12:20)=vn(2:11,12:20);
        %%%%% caculate v %%%%%
        for i=2:11
            for j=12:20
                vn(i,j)=(awv(i,j)*v(i,j-1)+aev(i,j)*v(i,j+1)+asv(i,j)*v(i-1,j)+anv(i,j)*v(i+1,j)+dx*(p(i-1,j)-p(i,j)))/av(i,j);
            end;
        end;
        error_v=abs(vn_0(2:11,12:20)-vn(2:11,12:20));
   end;
   un(21:1:32)=0;
   un(11,1:11)=0;
   un(11,22:32)=0;
   un(11:21,1)=0.4-un(i,2);
   un(11:21,32)=un(11:21,31);
   un(1:11,11)=-un(1:11,12);
   un(1:11,22)=un(1:11,21);
   un(1,11:22)=0;
   vn(:,1)=0;
   vn(:,31)=0;
   vn(22,:)=-vn(21,:);
   vn(11,1:10)=-vn(11,1:10);
   vn(11,22:31)=-vn(11,22:31);
   vn(1:11,11)=0;
   vn(1:11,21)=0;
   vn(1,11:21)=-vn(2,11:21);
   for i=12:20
        for j=2:31
            au_n(i,j)=4*miu+0.25*(un(i,j+1)-un(i,j-1))*density*dy+0.25*(vn(i+1,j-1)+vn(i+1,j)-vn(i,j)-vn(i,j-1))*density*dx;
        end;
   end;
   for i=2:11
        for j=12:21
            au_n(i,j)=4*miu+0.25*(un(i,j+1)-un(i,j-1))*density*dy+0.25*(vn(i+1,j-1)+vn(i+1,j)-vn(i,j)-vn(i,j-1))*density*dx;
        end;
   end;
   for i=12:21
       for j=2:30
           av_n(i,j)=4*miu+0.25*(un(i-1,j+1)+un(i,j+1)-un(i,j)-un(i-1,j))*density*dy+0.25*(vn(i+1,j)-vn(i-1,j))*density*dx;
       end;
   end;
   for i=2:11
        for j=12:20
            av_n(i,j)=4*miu+0.25*(un(i-1,j+1)+un(i,j+1)-un(i,j)-un(i-1,j))*density*dy+0.25*(vn(i+1,j)-vn(i-1,j))*density*dx;
        end;
   end;
   for i=12:20
       for j=2:30
           awp(i,j-1)=density*dy*(dy/au_n(i,j));
           aep(i,j+1)=density*dy*(dy/au_n(i,j+1));
           asp(i-1,j)=density*dx*(dx/av_n(i,j));
           anp(i+1,j)=density*dx*(dx/av_n(i+1,j));
           app(i,j)=awp(i,j-1)+aep(i,j+1)+asp(i-1,j)+anp(i+1,j);
           bp(i,j)=density*dy*(un(i,j)-un(i,j+1))+density*dx*(vn(i,j)-vn(i+1,j));
       end;
   end;
   for i=2:11
       for j=12:20
           awp(i,j-1)=density*dy*(dy/au_n(i,j));
           aep(i,j+1)=density*dy*(dy/au_n(i,j+1));
           asp(i-1,j)=density*dx*(dx/av_n(i,j));
           anp(i+1,j)=density*dx*(dx/av_n(i+1,j));
           app(i,j)=awp(i,j-1)+aep(i,j+1)+asp(i-1,j)+anp(i+1,j);
           bp(i,j)=density*dy*(un(i,j)-un(i,j+1))+density*dx*(vn(i,j)-vn(i+1,j));
       end;
   end;
   error_p=1;
   while error_p>1e-6
        pn_0(12:20,2:30)=pn(12:20,2:30);
        %%%%  correct pressure  %%%%
        for i=12:20
            for j=2:30
                pn(i,j)=(awp(i,j-1)*pn(i,j-1)+aep(i,j+1)*pn(i,j+1)+asp(i-1,j)*pn(i-1,j)+anp(i+1,j)*pn(i+1,j)+bp(i,j))/app(i,j);
            end;
        end;
        error_p=abs(pn_0(12:20,2:30)-pn(12:20,2:30));
    end;
    error_p=1;
   while error_p>1e-6
        pn_0(2:11,12:20)=pn(2:11,12:20);
        %%%%  correct pressure  %%%%
        for i=2:11
            for j=12:20
                pn(i,j)=(awp(i,j-1)*pn(i,j-1)+aep(i,j+1)*pn(i,j+1)+asp(i-1,j)*pn(i-1,j)+anp(i+1,j)*pn(i+1,j)+bp(i,j))/app(i,j);
            end;
        end;
        error_p=abs(pn_0(2:11,12:20)-pn(2:11,12:20));
    end;
    pn(11:21,30)=pn(11:21,31);
    pn(11:21,1)=0;
    pn(21,:)=0;
    pn(11,1:11)=0;
    pn(11,21:31)=0;
    pn(1:11,11)=0;
    pn(1:11,21)=0;
    pn(1,11:21)=0;
    error=abs(bp(i,j));

    for i=12:20
        for j=2:30
            p(i,j)=p(i,j)+pn(i,j);
        end;
    end;
    for i=2:11
        for j=12:20
            p(i,j)=p(i,j)+pn(i,j);
        end;
    end;
    p(11:21,1)=p(11:21,2);
    p(11:21,30)=p(11:21,31);
    p(21,:)=p(20,:);
    p(11,1:11)=p(12,1:11);
    p(11,21:31)=p(12,21:31);
    p(1:11,11)=p(1:11,12);
    p(1:11,21)=p(1:11,20);
    p(1,11:21)=p(2,11:21);

    for i=12:20
        for j=2:31
            u(i,j)=un(i,j)+(dy/au_n(i,j))*(pn(i,j-1)-pn(i,j));
        end;
    end;
    for i=2:11
        for j=12:21
            u(i,j)=un(i,j)+(dy/au_n(i,j))*(pn(i,j-1)-pn(i,j));
        end;
    end;
    for i=12:21
        for j=2:30
            v(i,j)=vn(i,j)+(dx/av_n(i,j))*(pn(i-1,j)-pn(i,j));
        end;
    end;
    for i=2:11
        for j=12:20
          v(i,j)=vn(i,j)+(dx/av_n(i,j))*(pn(i-1,j)-pn(i,j));
        end;
    end;
    u(11:21,1)=0.4-u(11:21,2);
    u(11,1:11)=0;
    u(11,22:32)=0;
    u(21,1:32)=0;
    u(11:21,32)=u(11:21,31);
    u(1:11,11)=-u(1:11,12);
    u(1:11,22)=u(1:11,21);
    u(1,11:22)=0;
    v(11:22,1)=0;
    v(11:22,31)=0;
    v(22,1:31)=-v(21,1:31);
    v(11,1:10)=-v(12,1:10);
    v(11,22:31)=-v(12,22:31);
    v(1:11,11)=0;
    v(1:11,21)=0;
    v(1,11:21)=-v(2,11:21);
    k=k+1;
end;
contour(u)
k

MATLAB 基础讨论
版块优秀回答者

入门

124 麦片

财富积分


50500


1

主题

121

帖子

21

最佳答案
发表于 7 天前 | 显示全部楼层
你的error就没变过,确定不是死循环?

新手

5 麦片

财富积分


050


1

主题

2

帖子

0

最佳答案
 楼主| 发表于 5 天前 | 显示全部楼层
帕拉代斯 发表于 2019-11-8 17:54
你的error就没变过,确定不是死循环?

您能说具体点吗?我也是初学者
您需要登录后才可以回帖 登录 | 注册

本版积分规则

关闭

站长推荐上一条 /3 下一条

快速回复 返回顶部 返回列表