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

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

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

充满着运动黏度系数υ=1cm²/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<<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;

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

}

运行结果：