搜索内容 搜索用户

# 5 麦片

050

1主题 0最佳答案

 % 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

# 319 麦片

50500

1主题 65最佳答案
• 关注者： 5

 你的error就没变过，确定不是死循环？

# 5 麦片

050

1主题 0最佳答案
楼主| 发表于 2019-11-10 19:14:58 | 显示全部楼层
 帕拉代斯 发表于 2019-11-8 17:54 你的error就没变过，确定不是死循环？ 您能说具体点吗？我也是初学者

 您需要登录后才可以回帖 登录 | 注册 本版积分规则 回帖后跳转到最后一页