下板不动, 上板匀速平板间流动（c++）【转载】

摘自《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 t,dz,dt,umax;

    int i,n,imax,imax1;

    imax=11;

    dz=0.1;

    dt=0.001;

    imax1=imax-1;

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

    {

        u[i]=0;

        u0[i]=0;

    }

    t=0;

    n=0;

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

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

            u[i]=u0[i]+dt/dz/dz*(u0[i+1]+u0[i-1]-2*u0[i]);

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

            u0[i]=u[i];

        if(n0==0)

        {

            cout<<endl<<setw(8)<<setiosflags(ios::fixed)<<setprecision(3)<<t;

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

                cout<<setw(8)<<setprecision(3)<<u0[i];

            cout<<endl<<"----------------------------------------------------------------------------------------------"<<endl;

        }

    }while(n<1000);

    return 0;

}

运行结果：

不同时刻的速度变化：

生成该图的matlab代码：

position=[0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1];

t_01=[0.000 0.086 0.202 0.386 0.684 1.162 1.903 3.017 4.633 6.905 10.000];

t_02=[0.000 0.543 1.128 1.794 2.573 3.490 4.554 5.763 7.097 8.524 10.000];

t_03=[0.000 0.826 1.668 2.544 3.463 4.435 5.463 6.543 7.668 8.825 10.000];

t_04=[0.000 0.935 1.876 2.829 3.799 4.789 5.799 6.829 7.876 8.935 10.000];

t_05=[0.000 0.976 1.954 2.936 3.925 4.921 5.925 6.936 7.954 8.976 10.000];

t_06=[0.000 0.991 1.983 2.976 3.972 4.970 5.972 6.976 7.983 8.991 10.000];

t_07=[0.000 0.997 1.994 2.991 3.990 4.989 5.990 6.991 7.994 8.997 10.000];

t_08=[0.000 0.999 1.998 2.997 3.996 4.996 5.996 6.997 7.998 8.999 10.000];

t_09=[0.000 1.000 1.999 2.999 3.999 4.998 5.999 6.999 7.999 9.000 10.000];

t_10=[0.000 1.000 2.000 3.000 3.999 4.999 5.999 7.000 8.000 9.000 10.000];

plot(position,t_01,position,t_02,position,t_03,position,t_04,position,t_05,position,t_06,position,t_07,position,t_08,position,t_09,position,t_10,'marker','.')

xlabel('板上位置(m)')

ylabel('速度(m/s)')

legend('t=0.1s','t=0.2s','t=0.3s','t=0.4s','t=0.5s','t=0.6s','t=0.7s','t=0.8s','t=0.9s','t=1s','Location','Best')

legend('boxoff')

box off

不同高度的速度变化图

生成该图的matlab代码：

t=[0 0.1    0.2    0.3    0.4    0.5    0.6    0.7    0.8    0.9    1]

position_0=[0 0    0    0    0    0    0    0    0    0    0]

position_1=[0 0.086    0.543    0.826    0.935    0.976    0.991    0.997    0.999    1    1]

position_2=[0 0.202    1.128    1.668    1.876    1.954    1.983    1.994    1.998    1.999    2]

position_3=[0 0.386    1.794    2.544    2.829    2.936    2.976    2.991    2.997    2.999    3]

position_4=[0 0.684    2.573    3.463    3.799    3.925    3.972    3.99    3.996    3.999    3.999]

position_5=[0 1.162    3.49    4.435    4.789    4.921    4.97    4.989    4.996    4.998    4.999]

position_6=[0 1.903    4.554    5.463    5.799    5.925    5.972    5.99    5.996    5.999    5.999]

position_7=[0 3.017    5.763    6.543    6.829    6.936    6.976    6.991    6.997    6.999    7]

position_8=[0 4.633    7.097    7.668    7.876    7.954    7.983    7.994    7.998    7.999    8]

position_9=[0 6.905    8.524    8.825    8.935    8.976    8.991    8.997    8.999    9    9]

position_10=[0 10    10    10    10    10    10    10    10    10    10]

plot(t,position_0,t,position_1,t,position_2,t,position_3,t,position_4,t,position_5,t,position_6,t,position_7,t,position_8,t,position_9,t,position_10,'marker','.')

ylim([0 11])

xlabel('时间(s)')

ylabel('速度(cm/s)')

legend('z=0cm','z=0.1cm','z=0.2cm','z=0.3cm','z=0.4cm','z=0.5cm','z=0.6cm','z=0.7cm','z=0.8cm','z=0.9cm','z=1cm','Location','BestOutside')

legend('boxoff')

box off