下板不动, 上板匀速平板间流动（Crank-Nicolson格式）
摘自《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,u0;
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=0;
u0=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=t*100;
else
u=10;
iter=0;
do
{
difmax=0;
iter+=1;
if(iter>100)
exit(1);
for(i=2;i<=imax1;i++)
{
temp=u;
u=u0*b/dt+1.0*b/2/dz/dz*(u+u+u0+u0-2*u0);
dif=fabs(temp-u);
if(dif>difmax)
difmax=dif;
}
}while(difmax>0.00001);
for(i=1;i<=imax;i++)
u0=u;
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
页:
[1]