方腔驱动流的simple算法（附matlab与c++程序）

已于 2023-03-01 16:05:29 修改
阅读量1.3w 收藏 133
点赞数 40
文章标签： 算法 cfd c++ matlab
于 2021-06-24 15:37:59 首次发布
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/Fasydeer/article/details/118187286
版权

一、问题：方腔驱动流

边界条件：

二、建立交错网格

三、边界条件的处理

U与V的边界条件看成虚拟节点与内部一层节点之和的一半等于零，即虚拟节点与内部一层节点互为相反数。

P的边界条件在网格中表现为第一层与第二层（或最后一层与倒数第二层）相等

四、有限差分

将当前时刻假设的U与V代入差分方程（1）与（2），求得当前时刻的U*与V*。

（其中i是横坐标，j是纵坐标）

（1）

（2）


五、修正

将U*与V*代入修正方程（3）中，得到P’，程序中利用高斯-塞德尔迭代法来求得方程组

（3）

，

将P’代入P=P*+0.1×P’中，得到当前时刻准确的P值。

再将P’代入公式（4）得到当前时刻准确的U与V。

（4）

再将当前时刻准确的U与V以及修正完毕的P代入差分方程（1）就可以得到下一时刻的U与V

然后开始循环

六、MATLAB代码

%写在前面：matlab一般无法吃满cpu，因此运算时间较长，可以通过暂停程序查看参数c来看已经迭代次数

m=101;

n=101;%修改除个位数外的位数来修改网格数的数量

P=zeros(m,n);

PX=P;

PP=P;

U=zeros(m,n+1);

V=zeros(m+1,n);

U(1,:)=1;%U的顶盖

VS=V;

US=U;%US是U星

RE=400;%高雷诺数建议增加网格数量进行运算

x=1/(m-1);

y=1/(n-1);

t=0.005;%△xyt的大小，△t<=△x收敛的可能性大，若e=NAN需将△t调小，例如Re=5000时可以调成0.001

e=1;

e1=1;%e1为d的最大值

e2=1;%e2用于选取各个d的值

e3=0;%e3用来保护e1

c=0;%外部迭代次数

l=0;%内部迭代次数

S1=P;%将U从交错网格上回归到P所在的原网格上

S2=P; %将V从交错网格上回归到P所在的原网格上

q=0;%总的内部循环次数

while  c<30000 %%%重要部分！！重要部分！！重要部分！！重要部分！！重要部分！！重要部分！！c<30000是流场收敛的大体迭代次数，已经可以满足大部分情况，也可以自己修改迭代次数缩短运行时间。若对于精度有极高的要求，可以将c<30000改为e1>0.0001，注意，这将消耗大量的运算时间。

     for a=2:n-1

    for b=2:m

  US(a,b)=U(a,b)-t*((U(a,b+1)^2-U(a,b-1)^2)/(2*x)+(U(a+1,b)*(V(a+1,b-1)+V(a+1,b))-U(a-1,b)*(V(a,b-1)+V(a,b)))/(4*y))+(t/RE)*((U(a,b+1)-2*U(a,b)+U(a,b-1))/(x*x)+(U(a+1,b)-2*U(a,b)+U(a-1,b))/(y*y))-(P(a,b)-P(a,b-1))*(t/x);

  VS(b,a)=V(b,a)-t*((V(b+1,a)^2-V(b-1,a)^2)/(2*y)+(V(b,a+1)*(U(b-1,a+1)+U(b,a+1))-V(b,a-1)*(U(b-1,a)+U(b,a)))/(4*x))+(t/RE)*((V(b,a+1)-2*V(b,a)+V(b,a-1))/(x*x)+(V(b+1,a)-2*V(b,a)+V(b-1,a))/(y*y))-(P(b,a)-P(b-1,a))*(t/y); 

    end

     end   %公式1与公式2

   US(:,1)=-US(:,2);

   US(:,n+1)=-US(:,n);

   VS(1,:)=-VS(2,:);

   VS(m+1,:)=-VS(m,:);

    PP=zeros(m,n);%每次内部迭代都要设置一个新的初始迭代值，于是PP于此更新

    PX=PP;

    q=q+l;%总的内部迭代次数

    l=0;

           while 1  %双while语句将无法实现套用的功能，需换成break语句

       l=l+1; %记录单次内部迭代次数     

     for d=2:n-1

    for f=2:m-1   

   PX(d,f)=1/4*(PP(d+1,f)+PP(d,f+1)+PX(d-1,f)+PX(d,f-1)-x*(US(d,f+1)-US(d,f)+VS(d+1,f)-VS(d,f))/t);%等式右边有PX为高斯-塞德尔迭代

    end

     end  %公式6-103

PX(1,:)=PX(2,:);

PX(:,1)=PX(:,2);

PX(m,:)=PX(m-1,:);

PX(:,n)=PX(:,n-1);

   e=max(max(abs(PP-PX)));%检验内部迭代收敛

   E(q+l,1)=e;%将每个误差存放于一个列向量中方便查看收敛趋势

   PP=PX;

   PX=PP;%高斯塞德尔迭代需要这一步来使第一步迭代准确

       if e<=0.0001

           break;%收敛的精度

       end

          if l>2000

           break;%收敛的迭代次数上限

          end

           end %求公式3的迭代

   P=P+0.1*PX;   %得到当前时刻的P值

P(1,:)=P(2,:);

P(:,1)=P(:,2);

P(m,:)=P(m-1,:);

P(:,n)=P(:,n-1); 

e3=0;

   for g=2:n-1

    for h=2:m-1   

   e2=abs(-US(g,h+1)/x+US(g,h)/x-VS(g+1,h)/x+VS(g,h)/x);

         if e2>e3

             e3=e2;

             e1=e2; %求d的最大值来与精度进行比较

         end

    end

   end  %公式3中d趋近于0

 for a=2:n-1

    for b=2:m   

  U(a,b)=US(a,b)-t/x*(PX(a,b)-PX(a,b-1));

  V(b,a)=VS(b,a)-t/y*(PX(b,a)-PX(b-1,a));

    end

 end %UV的修正方程即公式4

  U(:,1)=-U(:,2);

   U(:,n+1)=-U(:,n);

   V(1,:)=-V(2,:);

   V(m+1,:)=-V(m,:);

     for a=2:n-1

    for b=2:m

  US(a,b)=U(a,b)-t*((U(a,b+1)^2-U(a,b-1)^2)/(2*x)+(U(a+1,b)*(V(a+1,b-1)+V(a+1,b))-U(a-1,b)*(V(a,b-1)+V(a,b)))/(4*y))+(t/RE)*((U(a,b+1)-2*U(a,b)+U(a,b-1))/(x*x)+(U(a+1,b)-2*U(a,b)+U(a-1,b))/(y*y))-(P(a,b)-P(a,b-1))*(t/x);

  VS(b,a)=V(b,a)-t*((V(b+1,a)^2-V(b-1,a)^2)/(2*y)+(V(b,a+1)*(U(b-1,a+1)+U(b,a+1))-V(b,a-1)*(U(b-1,a)+U(b,a)))/(4*x))+(t/RE)*((V(b,a+1)-2*V(b,a)+V(b,a-1))/(x*x)+(V(b+1,a)-2*V(b,a)+V(b-1,a))/(y*y))-(P(b,a)-P(b-1,a))*(t/y); 

    end

     end  %将修正好的UVP再代入修正方程中求得下一时刻的UV

   US(:,1)=-US(:,2);

   US(:,n+1)=-US(:,n);

   VS(1,:)=-VS(2,:);

   VS(m+1,:)=-VS(m,:);

U=US;

V=VS;

c=c+1;%外部迭代次数

end

for a=1:n

    for b=1:m 

        S1(a,b)=(US(a,b)+US(a,b+1))/2;

        S2(b,a)=(VS(b,a)+VS(b+1,a))/2;%将UV从交错网格上重新放到P的网格中

    end

end %将UV化成101*101网格

xx=0:m-1;

yy=0:n-1;

[X,Y]=meshgrid(xx,yy);%生成坐标

figure
pcolor(X,Y,P);%生成压强图

shading interp%输出图像

figure

streamslice(S1,S2);%生成流线图

七、c++代码

#include<iostream>

#include <fstream>

#include<math.h>

using namespace std;

double A[2000][2000],B[2000][2000];

double A1[2000][2000],B1[2000][2000];

double P[2000][2000],P1[2000][2000];

double P2[2000][2000],P3[2000][2000];

double A2[2000][2000],B2[2000][2000];

double P4[2000][2000],P5[2000][2000];

void main()

{int i,j,n,m;

double x=0.025,y=0.025,t=0.01;/*空间、时间步长*/

double uva,uvb,vua,vub;/*存储平均速度的值*/

double Re=300,e;

double a,b,c,d;/*表示求P的系数*/

a=2*(t/(x*x)+t/(y*y));

b=-t/(x*x);

c=-t/(y*y);

for(i=1;i<=41;i++)

{for(j=1;j<=41;j++)

{P[i][j]=0;}}/*给定初始P*的值*/

for(i=1;i<=40;i++)

{for(j=1;j<=40;j++)

{A[i][j]=0;

 B[i][j]=0;}}/*给定n时刻u*,v*的值*/

for(j=0;j<=41;j++)

{B[1][j]=0;

B[41][j]=0;}

for(i=0;i<=41;i++)

{A[i][1]=0;

A[i][41]=1;}

for(j=1;j<=41;j++)

{ A[0][j]=0-A[1][j];

  A[41][j]=0-A[40][j]; }

for(i=0;i<=41;i++)

{A[i][0]=0-A[i][2];

A[i][42]=2-A[i][40];}

for(i=1;i<=41;i++)

{B[i][0]=0-B[i][1];

 B[i][41]=0-B[i][40];}

for(j=0;j<=41;j++)

{ B[0][j]=0-B[2][j];

  B[42][j]=0-B[40][j];}/*处理u*,v*的边界值*/

int a1=0,a2=0;

double f,s,sum1=0,sum2=0;

for(n=1;n<=500;n++)

{for(i=1;i<=40;i++)

  {for(j=1;j<=41;j++)

   {uva=(A[i][j+1]+A[i][j])*(B[i+1][j]+B[i][j])/4;

    uvb=(A[i][j]+A[i][j-1])*(B[i+1][j-1]+B[i][j-1])/4;

    A1[i][j]=A[i][j]-(A[i+1][j]*A[i+1][j]-A[i][j]*A[i][j])/x*t-(uva-uvb)/y*t+

                       1/Re*((A[i+1][j]-2*A[i][j]+A[i-1][j])/x/x+

              (A[i][j+1]-2*A[i][j]+A[i][j-1])/y/y)*t-t/x*(P[i+1][j]-P[i][j]);}}

 for(i=1;i<=41;i++)

  {for(j=1;j<=40;j++)

   {vua=(A[i][j+1]+A[i][j])*(B[i+1][j]+B[i][j])/4;

    vub=(A[i-1][j+1]+A[i-1][j])*(B[i-1][j]+B[i][j])/4;

    B1[i][j]=B[i][j]-(vua-vub)/x*t-(B[i][j+1]*B[i][j+1]-B[i][j]*B[i][j])/y*t+

                       1/Re*((B[i+1][j]-2*B[i][j]+B[i-1][j])/x/x+

              (B[i][j+1]-2*B[i][j]+B[i][j-1])/y/y)*t-t/y*(P[i][j+1]-P[i][j]);}}

             /*求出un+1*,vn+1**/

for(i=0;i<=42;i++)

{for(j=0;j<=42;j++)

 {P1[i][j]=0;

  P2[i][j]=0;

  P3[i][j]=0;}}/*给定初始P’的值为了迭代*/

 d=1/x*(A1[4][5]-A1[3][5])+1/y*(B1[4][5]-B1[4][4]);

 P1[4][5]=-b/a*P1[5][5]-b/a*P1[3][5]-c/a*P1[4][6]-c/a*P1[4][4]-d/a;

for(i=2;i<=40;i++)

{for(j=2;j<=40;j++)

{   if(i==2&&j!=2&&j!=40)

 {d=1/x*(A1[i][j]-A1[i-1][j])+1/y*(B1[i][j]-B1[i][j-1]);

P4[i][j]=-b/(a+b)*P1[i+1][j]-c/(a+b)*P1[i][j+1]-c/(a+b)*P2[i][j-1]-d/(a+b);

P2[i][j]=P1[i][j]+0.2*(P4[i][j]-P1[i][j]);}

else if(i==40&&j!=2&&j!=40)

{d=1/x*(A1[i][j]-A1[i-1][j])+1/y*(B1[i][j]-B1[i][j-1]);

P4[i][j]=-b/(a+b)*P2[i-1][j]-c/(a+b)*P1[i][j+1]-c/(a+b)*P2[i][j-1]-d/(a+b);

P2[i][j]=P1[i][j]+0.2*(P4[i][j]-P1[i][j]);}

else if(j==2&&i!=2&&i!=40)

{d=1/x*(A1[i][j]-A1[i-1][j])+1/y*(B1[i][j]-B1[i][j-1]);

P4[i][j]=-b/(a+c)*P1[i+1][j]-b/(a+c)*P2[i-1][j]-c/(a+c)*P1[i][j+1]-d/(a+c);

P2[i][j]=P1[i][j]+0.2*(P4[i][j]-P1[i][j]);}

else if(j==40&&i!=2&&i!=40)

{d=1/x*(A1[i][j]-A1[i-1][j])+1/y*(B1[i][j]-B1[i][j-1]);

P4[i][j]=-b/(a+c)*P1[i+1][j]-b/(a+c)*P2[i-1][j]-c/(a+c)*P2[i][j-1]-d/(a+c);

P2[i][j]=P1[i][j]+0.2*(P4[i][j]-P1[i][j]);}

else if(i==2&&j==2)

{d=1/x*(A1[i][j]-A1[i-1][j])+1/y*(B1[i][j]-B1[i][j-1]);

P4[i][j]=-b/(a+b+c)*P1[i+1][j]-c/(a+b+c)*P1[i][j+1]-d/(a+b+c);

P2[i][j]=P1[i][j]+0.2*(P4[i][j]-P1[i][j]);}

else if(i==2&&j==40)

{d=1/x*(A1[i][j]-A1[i-1][j])+1/y*(B1[i][j]-B1[i][j-1]);

P4[i][j]=-b/(a+b+c)*P1[i+1][j]-c/(a+b+c)*P2[i][j-1]-d/(a+b+c);

P2[i][j]=P1[i][j]+0.2*(P4[i][j]-P1[i][j]);}

else if(i==40&&j==2)

{d=1/x*(A1[i][j]-A1[i-1][j])+1/y*(B1[i][j]-B1[i][j-1]);

P4[i][j]=-b/(a+b+c)*P2[i-1][j]-c/(a+b+c)*P1[i][j+1]-d/(a+b+c);

P2[i][j]=P1[i][j]+0.2*(P4[i][j]-P1[i][j]);}

else if(i==40&&j==40)

{d=1/x*(A1[i][j]-A1[i-1][j])+1/y*(B1[i][j]-B1[i][j-1]);

P4[i][j]=-b/(a+b+c)*P2[i-1][j]-c/(a+b+c)*P2[i][j-1]-d/(a+b+c);

P2[i][j]=P1[i][j]+0.2*(P4[i][j]-P1[i][j]);}

else{d=1/x*(A1[i][j]-A1[i-1][j])+1/y*(B1[i][j]-B1[i][j-1]);

P4[i][j]=-b/a*P1[i+1][j]-b/a*P2[i-1][j]-c/a*P1[i][j+1]-c/a*P2[i][j-1]-d/a;

P2[i][j]=P1[i][j]+0.2*(P4[i][j]-P1[i][j]);}}}

for(i=2;i<=40;i++)

{for(j=2;j<=40;j++)

{     if(i==2&&j!=2&&j!=40)

 {d=1/x*(A1[i][j]-A1[i-1][j])+1/y*(B1[i][j]-B1[i][j-1]);

P5[i][j]=-b/(a+b)*P2[i+1][j]-c/(a+b)*P2[i][j+1]-c/(a+b)*P3[i][j-1]-d/(a+b);

P3[i][j]=P2[i][j]+0.2*(P5[i][j]-P2[i][j]);}

else if(i==40&&j!=2&&j!=40)

{d=1/x*(A1[i][j]-A1[i-1][j])+1/y*(B1[i][j]-B1[i][j-1]);

P5[i][j]=-b/(a+b)*P3[i-1][j]-c/(a+b)*P2[i][j+1]-c/(a+b)*P3[i][j-1]-d/(a+b);

P3[i][j]=P2[i][j]+0.2*(P5[i][j]-P2[i][j]);}

else if(j==2&&i!=2&&i!=40)

{d=1/x*(A1[i][j]-A1[i-1][j])+1/y*(B1[i][j]-B1[i][j-1]);

P5[i][j]=-b/(a+c)*P2[i+1][j]-b/(a+c)*P3[i-1][j]-c/(a+c)*P2[i][j+1]-d/(a+c);

    P3[i][j]=P2[i][j]+0.2*(P5[i][j]-P2[i][j]);}

else if(j==40&&i!=2&&i!=40)

{d=1/x*(A1[i][j]-A1[i-1][j])+1/y*(B1[i][j]-B1[i][j-1]);

P5[i][j]=-b/(a+c)*P2[i+1][j]-b/(a+c)*P3[i-1][j]-c/(a+c)*P3[i][j-1]-d/(a+c);

P3[i][j]=P2[i][j]+0.2*(P5[i][j]-P2[i][j]);}

else if(i==2&&j==2)

{d=1/x*(A1[i][j]-A1[i-1][j])+1/y*(B1[i][j]-B1[i][j-1]);

P5[i][j]=-b/(a+b+c)*P2[i+1][j]-c/(a+b+c)*P2[i][j+1]-d/(a+b+c);

P3[i][j]=P2[i][j]+0.2*(P5[i][j]-P2[i][j]);}

else if(i==2&&j==40)

{d=1/x*(A1[i][j]-A1[i-1][j])+1/y*(B1[i][j]-B1[i][j-1]);

P5[i][j]=-b/(a+b+c)*P2[i+1][j]-c/(a+b+c)*P3[i][j-1]-d/(a+b+c);

P3[i][j]=P2[i][j]+0.2*(P5[i][j]-P2[i][j]);}

else if(i==40&&j==2)

{d=1/x*(A1[i][j]-A1[i-1][j])+1/y*(B1[i][j]-B1[i][j-1]);

P5[i][j]=-b/(a+b+c)*P3[i-1][j]-c/(a+b+c)*P2[i][j+1]-d/(a+b+c);

P3[i][j]=P2[i][j]+0.2*(P5[i][j]-P2[i][j]);}

else if(i==40&&j==40)

{d=1/x*(A1[i][j]-A1[i-1][j])+1/y*(B1[i][j]-B1[i][j-1]);

P5[i][j]=-b/(a+b+c)*P3[i-1][j]-c/(a+b+c)*P3[i][j-1]-d/(a+b+c);

P3[i][j]=P2[i][j]+0.2*(P5[i][j]-P2[i][j]);}

else{d=1/x*(A1[i][j]-A1[i-1][j])+1/y*(B1[i][j]-B1[i][j-1]);

P5[i][j]=-b/a*P2[i+1][j]-b/a*P3[i-1][j]-c/a*P2[i][j+1]-c/a*P3[i][j-1]-d/a;

P3[i][j]=P2[i][j]+0.2*(P5[i][j]-P2[i][j]);}}}

for(m=1;m<=200;m++)

{for(i=2;i<=40;i++)

{for(j=2;j<=40;j++)

{    if(i==2&&j!=2&&j!=40)

 {d=1/x*(A1[i][j]-A1[i-1][j])+1/y*(B1[i][j]-B1[i][j-1]);

P4[i][j]=-b/(a+b)*P3[i+1][j]-c/(a+b)*P3[i][j+1]-c/(a+b)*P2[i][j-1]-d/(a+b);

P2[i][j]=P3[i][j]+0.2*(P4[i][j]-P3[i][j]);}

else if(i==40&&j!=2&&j!=40)

{d=1/x*(A1[i][j]-A1[i-1][j])+1/y*(B1[i][j]-B1[i][j-1]);

P4[i][j]=-b/(a+b)*P2[i-1][j]-c/(a+b)*P3[i][j+1]-c/(a+b)*P2[i][j-1]-d/(a+b);

P2[i][j]=P3[i][j]+0.2*(P4[i][j]-P3[i][j]);}

else if(j==2&&i!=2&&i!=40)

{d=1/x*(A1[i][j]-A1[i-1][j])+1/y*(B1[i][j]-B1[i][j-1]);

P4[i][j]=-b/(a+c)*P3[i+1][j]-b/(a+c)*P2[i-1][j]-c/(a+c)*P3[i][j+1]-d/(a+c);

P2[i][j]=P3[i][j]+0.2*(P4[i][j]-P3[i][j]);}

else if(j==40&&i!=2&&i!=40)

{d=1/x*(A1[i][j]-A1[i-1][j])+1/y*(B1[i][j]-B1[i][j-1]);

P4[i][j]=-b/(a+c)*P3[i+1][j]-b/(a+c)*P2[i-1][j]-c/(a+c)*P2[i][j-1]-d/(a+c);

P2[i][j]=P3[i][j]+0.2*(P4[i][j]-P3[i][j]);}

else if(i==2&&j==2)

{d=1/x*(A1[i][j]-A1[i-1][j])+1/y*(B1[i][j]-B1[i][j-1]);

P4[i][j]=-b/(a+b+c)*P3[i+1][j]-c/(a+b+c)*P3[i][j+1]-d/(a+b+c);

P2[i][j]=P3[i][j]+0.2*(P4[i][j]-P3[i][j]);}

else if(i==2&&j==40)

{d=1/x*(A1[i][j]-A1[i-1][j])+1/y*(B1[i][j]-B1[i][j-1]);

P4[i][j]=-b/(a+b+c)*P3[i+1][j]-c/(a+b+c)*P2[i][j-1]-d/(a+b+c);

P2[i][j]=P3[i][j]+0.2*(P4[i][j]-P3[i][j]);}

else if(i==40&&j==2)

{d=1/x*(A1[i][j]-A1[i-1][j])+1/y*(B1[i][j]-B1[i][j-1]);

P4[i][j]=-b/(a+b+c)*P2[i-1][j]-c/(a+b+c)*P3[i][j+1]-d/(a+b+c);

P2[i][j]=P3[i][j]+0.2*(P4[i][j]-P3[i][j]);}

else if(i==40&&j==40)

{d=1/x*(A1[i][j]-A1[i-1][j])+1/y*(B1[i][j]-B1[i][j-1]);

P4[i][j]=-b/(a+b+c)*P2[i-1][j]-c/(a+b+c)*P2[i][j-1]-d/(a+b+c);

P2[i][j]=P3[i][j]+0.2*(P4[i][j]-P3[i][j]);}

else{d=1/x*(A1[i][j]-A1[i-1][j])+1/y*(B1[i][j]-B1[i][j-1]);

P4[i][j]=-b/a*P3[i+1][j]-b/a*P2[i-1][j]-c/a*P3[i][j+1]-c/a*P2[i][j-1]-d/a;

P2[i][j]=P3[i][j]+0.2*(P4[i][j]-P3[i][j]);}}}

for(i=2;i<=40;i++)

{for(j=2;j<=40;j++)

{    if(i==2&&j!=2&&j!=40)

 {d=1/x*(A1[i][j]-A1[i-1][j])+1/y*(B1[i][j]-B1[i][j-1]);

P5[i][j]=-b/(a+b)*P2[i+1][j]-c/(a+b)*P2[i][j+1]-c/(a+b)*P3[i][j-1]-d/(a+b);

P3[i][j]=P2[i][j]+0.2*(P5[i][j]-P2[i][j]);}

else if(i==40&&j!=2&&j!=40)

{d=1/x*(A1[i][j]-A1[i-1][j])+1/y*(B1[i][j]-B1[i][j-1]);

P5[i][j]=-b/(a+b)*P3[i-1][j]-c/(a+b)*P2[i][j+1]-c/(a+b)*P3[i][j-1]-d/(a+b);

P3[i][j]=P2[i][j]+0.2*(P5[i][j]-P2[i][j]);}

else if(j==2&&i!=2&&i!=40)

{d=1/x*(A1[i][j]-A1[i-1][j])+1/y*(B1[i][j]-B1[i][j-1]);

P5[i][j]=-b/(a+c)*P2[i+1][j]-b/(a+c)*P3[i-1][j]-c/(a+c)*P2[i][j+1]-d/(a+c);

    P3[i][j]=P2[i][j]+0.2*(P5[i][j]-P2[i][j]);}

else if(j==40&&i!=2&&i!=40)

{d=1/x*(A1[i][j]-A1[i-1][j])+1/y*(B1[i][j]-B1[i][j-1]);

P5[i][j]=-b/(a+c)*P2[i+1][j]-b/(a+c)*P3[i-1][j]-c/(a+c)*P3[i][j-1]-d/(a+c);

P3[i][j]=P2[i][j]+0.2*(P5[i][j]-P2[i][j]);}

else if(i==2&&j==2)

{d=1/x*(A1[i][j]-A1[i-1][j])+1/y*(B1[i][j]-B1[i][j-1]);

P5[i][j]=-b/(a+b+c)*P2[i+1][j]-c/(a+b+c)*P2[i][j+1]-d/(a+b+c);

P3[i][j]=P2[i][j]+0.2*(P5[i][j]-P2[i][j]);}

else if(i==2&&j==40)

{d=1/x*(A1[i][j]-A1[i-1][j])+1/y*(B1[i][j]-B1[i][j-1]);

P5[i][j]=-b/(a+b+c)*P2[i+1][j]-c/(a+b+c)*P3[i][j-1]-d/(a+b+c);

P3[i][j]=P2[i][j]+0.2*(P5[i][j]-P2[i][j]);}

else if(i==40&&j==2)

{d=1/x*(A1[i][j]-A1[i-1][j])+1/y*(B1[i][j]-B1[i][j-1]);

P5[i][j]=-b/(a+b+c)*P3[i-1][j]-c/(a+b+c)*P2[i][j+1]-d/(a+b+c);

P3[i][j]=P2[i][j]+0.2*(P5[i][j]-P2[i][j]);}

else if(i==40&&j==40)

{d=1/x*(A1[i][j]-A1[i-1][j])+1/y*(B1[i][j]-B1[i][j-1]);

P5[i][j]=-b/(a+b+c)*P3[i-1][j]-c/(a+b+c)*P3[i][j-1]-d/(a+b+c);

P3[i][j]=P2[i][j]+0.2*(P5[i][j]-P2[i][j]);}

else{d=1/x*(A1[i][j]-A1[i-1][j])+1/y*(B1[i][j]-B1[i][j-1]);

P5[i][j]=-b/a*P2[i+1][j]-b/a*P3[i-1][j]-c/a*P2[i][j+1]-c/a*P3[i][j-1]-d/a;

P3[i][j]=P2[i][j]+0.2*(P5[i][j]-P2[i][j]);}}}

for(i=2;i<=40;i++)

{for(j=2;j<=40;j++)

{f=abs(P3[i][j]-P2[i][j]);

sum1=sum1+f*f;}}

e=sqrt(sum1)/1521;

if(e<=0.000001&&m!=1)

{m=1001;}

sum1=0;

a1++;}/*运用松弛迭代法求出p’*/

for(i=1;i<=41;i++)

{P3[i][1]=P3[i][2];

P3[i][41]=P3[i][40];}

for(j=2;j<=40;j++)

{P3[1][j]=P3[2][j];

P3[41][j]=P3[40][j];}/*p'的边界条件处理*/

for(i=1;i<=40;i++)

{for(j=1;j<=41;j++)

{A2[i][j]=-t/x*(P3[i+1][j]-P3[i][j]);}}

for(i=1;i<=41;i++)

{for(j=1;j<=40;j++)

{B2[i][j]=-t/y*(P3[i][j+1]-P3[i][j]);}}/*求u',v'的值*/

for(i=1;i<42;i++)

{for(j=1;j<42;j++)

{P[i][j]=P[i][j]+P3[i][j];}}/*pn+1=P*+P'*/

for(i=1;i<=40;i++)

{for(j=1;j<=41;j++)

 {A[i][j]=A1[i][j]+A2[i][j];}}

for(i=1;i<=41;i++)

{for(j=1;j<=40;j++)

 {B[i][j]=B1[i][j]+B2[i][j];}}

for(j=0;j<=41;j++)

{B[1][j]=0;

 B[41][j]=0;}

for(i=0;i<=41;i++)

{A[i][1]=0;

 A[i][41]=1;}

for(j=1;j<=41;j++)

{ A[0][j]=0-A[1][j];

  A[41][j]=0-A[40][j]; }

for(i=0;i<=41;i++)

{ A[i][0]=0-A[i][2];

  A[i][42]=2-A[i][40];}

for(i=1;i<=41;i++)

{ B[i][0]=0-B[i][1];

  B[i][41]=0-B[i][40];}

for(j=0;j<=41;j++)

{ B[0][j]=0-B[2][j];

  B[42][j]=0-B[40][j];}/*处理u，v边界虚拟点*/

for(i=2;i<=40;i++)

{for(j=2;j<=40;j++)

{d=abs(1/x*(A1[i][j]-A1[i-1][j])+1/y*(B1[i][j]-B1[i][j-1]));

sum2=sum2+d*d;}}

s=sqrt(sum2)/1521;

if(s<=0.00001&&n!=1)

{n=1001;}

sum2=0;

a2++;}/*simple算法*/

cout<<e<<" "<<s<<" "<<a1<<" "<<a2;

int g=0,h=0;

string sstr= "C:\\Users\\12489\\Desktop\\re407.txt";

 ofstream fout(sstr);

for(j=1;j<=41;j++)

 {for(i=1;i<=41;i++)

{fout <<g*0.025<<" "<<h*0.025<<" "<<(A[i][j]+A[i-1][j])/2<<" " <<(B[i][j]+B[i][j-1])/2<<" "<<P[i][j]<<"\n";

  g++;

  if(g==41)

   {g=0;

     h++;}}}

fout << endl;

fout.close();} /*将数据写入文件*/

八、运算结果（上流线图下压强分布图）

RE=50    101*101    △t=0.005s

RE=400    101*101    △t=0.005s

RE=1000    201*201    △t=0.005s

RE=5000     501*501    △t=0.001s

感谢西南交通大学力学与工程学院陈浩老师对本程序的指导

卡萨甲
关注
  • 40
    点赞
  • 踩
  • 133
    收藏
    觉得还不错? 一键收藏
  • 50
    评论
Simple算法算例
06-02
Simple算法3个典型算例，详细讲解。二维、结构网格单相流SIMPLE 算法的C语言实现。
SIMPLE算法Matlab代码
01-25
%--------------------------------------------------------% % SIMPLE Algorithm - staggered mesh % Incompressible Lid-driven Flow in Cavity %--------------------------------------------------------%
方腔驱动流源代码(SIMPLE算法).rar_simple 方腔流_simple 算法_方腔simple_方腔流_驱动流
07-15
从方程离散到代码，到最终结果的后处理都有详细的描述
基于多网格的 CFD 求解盖驱动腔问题附matlab代码
m0_60703264的博客
06-15 65
实现了多种时间、空间数值格式以及多种湍流模式与嵌补网格技术的有机结合.开发了基于嵌补网格技术,适合复杂外型亚跨超音速流场模拟的计算软件.经多个典型算例的验证,证明软件及所使用的方法具有良好的精度和可靠性.对复杂流场具有很强的模拟能力.
cavityFlow2D.zip_simple算法_方腔 simple_方腔simple_方腔流 matlab_方腔驱动
07-14
用matlab计算方腔驱动流.并且使用了simple算法对压力进行修正
有限体积法计算方腔流(F)_flow_流场计算_流体力学_cavity_方腔_
10-02
张德良有限体积法计算方腔流场fortran及c源码
力学流体simple
a_zxswer的博客
07-03 2678
目录问题问题描述离散动量方程SIMPLE算法基本思想Simple算法的流程图： 问题 顶盖驱动流动 问题描述 假设 的方腔内充满粘性不可压缩流体，左、右、下壁固定，上壁以 运动， 方腔内初始时刻压力和密度为它周围壁面（左右壁面和底面）固定不动，上壁面以量纲为一的速度沿着上壁面方向自左向右运动。方腔如图1所示。 图1 方腔内流动示意图 基本方程组、初始条件和边界条件 设流体是黏性流体。二维方腔流动问题的动量方程为： 连续方程为：
有限体积法（12）——SIMPLE算法
薛定谔的alpha狗
09-29 6940
承接上篇《交错网格》一文，本文介绍流动方程在交错网格上的离散以及SIMPLE算法。 方程离散 二维稳态的动量方程和连续性方程如下： ∂∂x(ρuu)+∂∂y(ρvu)=∂∂x(μ∂u∂x)+∂∂y(μ∂u∂y)−∂p∂x+Su(1-a)\begin{aligned} \frac{\partial}{\partial x} (\rho u u) + \frac{\partial}{\partial y} (\rho vu) &= \frac{\partial}{\partial x} \left( \
2.2 SIMPLE系列算法 | 2.3 PISO算法（OpenFOAM理论笔记系列）
热门推荐
CloudBird07的博客
07-31 1万+
2.2 SIMPLE系列算法 2.2.1标准SIMPLE算法 SIMPLE算法(Semi-Implicit Method for PressureLinked Equations)1最初被设计用来求解稳态问题，即控制方程中不包含瞬态项的计算。按照1.3.3节的约定，我们假设计算开始的时候有初始的压力和速度值Po,Uo⃗P^o,\vec{U^o}Po,Uo,待求的真值为Pn,Un⃗P^n,\vec{U^n}Pn,Un。实际上，速度的初始值可以完全随意给定，因为我们完全首先可以将初始压力带入式2.1解出一个速度
有限体积法（13）——SIMPLE算法的一维算例
薛定谔的alpha狗
11-19 6304
模型是一个二维平面喷嘴，如图，假设流动是稳态且无粘的，流体不可压。 采用后向交错网格，如下图，压力节点5个，速度节点4个，入口给定滞止压强，出口给定静压
LBGK_LBGK_方腔流_
09-29
LBGK 常用代码，格子玻尔兹曼方法描述二维方腔顶盖驱动流，适合初学者
方腔顶盖驱动流动.doc
10-08
方腔顶盖驱动流动.doc
SIMPLE算法-MATLAB.zip_simple_simple matlab_simple 算法_simple-matlab
07-14
SIMPLE算法编程语言，使用Matlab对SIMPLE算法进行了阐述。
MATLAB传热基础算例.rar_matlab房间传热_传热 matlab_换热matlab算例_方腔_非稳态
07-15
给出了传热学上的几个基础算例，包括一维稳态/非稳态，二维稳态/非稳态，以及方腔驱动流的matlab算法
标准方腔流动.zip_不可压缩流_方腔_方腔流 matlab_方腔流动 matlab_流函数matlab
07-13
用来计算流体力学中非常典型的不可压缩方腔流动。包含MATLAB程序绘制流函数图。
SIMPLE算法计算二维方腔流动（C语言描述）
11-18
SIMPLE算法，全名为压力耦合方程组的半隐式方法（Semi-Implicit Method for Pressure Linked Equations），是计算流体力学中一种被广泛使用的求解流场的数值方法。二维方腔流动是各计算流体力学教材经典算例。基于...
Java | Leetcode Java题解之第35题搜索插入位置
m0_57195758的博客
04-18 445
Java | Leetcode Java题解之第35题搜索插入位置
DFS与回溯专题：路径总和问题
孤舟一叶～的博客
04-19 759
对二叉树进行dfs搜索，递归计算每条路径的节点值之和，当某个节点的左右子节点都为空时，说明已经搜索完成某一条路径，将它与目标值进行比较，若相等，则为true。路径恢复： 向上回溯前，需要将当前节点从路径 list 中删除，即执行list.remove(list.size() - 1)。为什么要写list_all.add(new ArrayList(list))，而不是list_all.add(list)？注意：解释里面的sum.add(path)就是list_all.add(list)
算法训练营day21
最新发布
wugong_true的博客
04-23 242
接下来，根据计数count与最大计数maxCount的比较，更新众数的结果数组resList和最大计数maxCount。如果count大于maxCount，则清空resList并将当前节点值添加进去，并更新maxCount为count；接着，根据resList中的众数，创建一个整型数组res，并将其返回作为结果。然后，通过与前一个节点的值进行比较，判断当前节点值出现的次数。如果当前节点与前一个节点的值不同，说明出现了一个新的值，将计数count重置为1；，计算相邻两个元素的差值，并将较小的差值更新到。
方腔驱动流的simple算法
06-12
方腔驱动流是指在一个矩形方腔中通过施加外力来驱动流体运动的问题。简单的算法如下： 1. 建立数值模型：将方腔划分为网格，采用有限体积法或有限元法等数值方法建立模型。 2. 设置边界条件：对于流场中的边界，需要设置边界条件，如速度、压力等，以模拟实际情况。 3. 计算流场：采用数值方法求解流场中的速度场和压力场，可以使用迭代法或直接解法等方法。 4. 分析结果：根据计算结果分析流场特征，如涡旋、流线等，并对流体力学现象进行深入研究。 需要注意的是，方腔驱动流问题是一个复杂的流体力学问题，需要结合实际情况进行综合分析，建立合适的数值模型和边界条件，以获得准确可靠的计算结果。

“相关推荐”对你有帮助么？

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 50
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值