摘自《FLUENT流体工程仿真计算实例与分析》，程序略有修改
两个间距为1cm水平平板，如下图所示：
充满着运动黏度系数υ=1cm2/s的液体。上板做水平运动并在0.1s时间内，速度线性由0线性地增加到10cm/s，如下图所示：
通过对N—S方程的简化，可由下面的抛物线方程来描述
流动区域在z=0和z=0.1cm之间，初始条件为u(z,0)=0cm/s，边界条件为：u|z=0=0，u|z=1=Ucm/s，U为上板移动速度
取Δz=0.1cm，Δt=0.001s，沿板的铅垂方向把空间离散为11等步长的节点，计算各点处的流速，边界条件可取节点i=1(z=0)和i=11(z=1cm)的流速
求解方程的程序代码：
- #include
- #include
- #include
- using namespace std;
- int main()
- {
- float u[100],u0[100];
- float b,t,dz,dt,dif,difmax,temp;
- int imax,imax1,iter,i,n;
- dz=0.1;
- dt=0.001;
- imax=11;
- for(i=1;i<=imax;i++)
- {
- u[i]=0;
- u0[i]=0;
- }
- imax1=imax-1;
- n=0;
- t=0;
- b=1.0/(1.0/dt+1.0/dz/dz);
- cout<
- cout<
- <
- <
- <
- <
- <
- <
- <
- <
- <
- <
- <
- <
- do
- {
- t+=dt;
- n+=1;
- if(t<0.1)
- u[imax]=t*100;
- else
- u[imax]=10;
- iter=0;
- do
- {
- difmax=0;
- iter+=1;
- if(iter>100)
- exit(1);
- for(i=2;i<=imax1;i++)
- {
- temp=u[i];
- u[i]=u0[i]*b/dt+1.0*b/2/dz/dz*(u[i+1]+u[i-1]+u0[i+1]+u0[i-1]-2*u0[i]);
- dif=fabs(temp-u[i]);
- if(dif>difmax)
- difmax=dif;
- }
- }while(difmax>0.00001);
- for(i=1;i<=imax;i++)
- u0[i]=u[i];
- if(n0==0)
- {
- cout<
- for(i=1;i<=imax;i++)
- cout<
- cout<
- }
- }while(n<1000);
- cout<
- return 0;
- }
运行结果：
转自：http://blog.sina.com.cn/s/blog_14d64daa10102wkq1.html
|
