下板不动, 上板匀速平板间流动（Crank-Nicolson格式）【转载】

摘自《FLUENT流体工程仿真计算实例与分析》，程序略有修改

两个间距为1cm水平平板，如下图所示：

下板不动, 上板匀速平板间流动（Crank-Nicolson格式）【转载】-LMLPHP 上板匀速平板间流动（Crank-Nicolson格式）【转载】">

充满着运动黏度系数υ=1cm/s的液体。上板做水平运动并在0.1s时间内，速度线性由0线性地增加到10cm/s，如下图所示：

下板不动, 上板匀速平板间流动（Crank-Nicolson格式）【转载】-LMLPHP 上板匀速平板间流动（Crank-Nicolson格式）【转载】" title="下板不动, 上板匀速平板间流动（Crank-Nicolson格式）【转载】">

通过对N—S方程的简化，可由下面的抛物线方程来描述

下板不动, 上板匀速平板间流动（Crank-Nicolson格式）【转载】-LMLPHP 上板匀速平板间流动（Crank-Nicolson格式）【转载】">

流动区域在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

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<<endl;

cout<<setiosflags(ios::left)

<<setw(8)<<"t(s)"

<<setw(8)<<"0"

<<setw(8)<<"0.1"

<<setw(8)<<"0.2"

<<setw(8)<<"0.3"

<<setw(8)<<"0.4"

<<setw(8)<<"0.5"

<<setw(8)<<"0.6"

<<setw(8)<<"0.7"

<<setw(8)<<"0.8"

<<setw(8)<<"0.9"

<<setw(8)<<"1"<<endl;

{

t+=dt;

n+=1;

if(t<0.1)

u[imax]=t*100;

else

u[imax]=10;

iter=0;

{

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<<endl<<setw(8)<<t;

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

cout<<setw(8)<<setiosflags(ios::fixed)<<setprecision(3)<<u0[i];

cout<<resetiosflags(ios::fixed)<<endl<<"----------------------------------------------------------------------------------------------";

}

}while(n<1000);

cout<<endl;

return
0;

}

运行结果：