摘自《FLUENT流体工程仿真计算实例与分析》，程序略有修改
两个间距为1cm水平平板，如下图所示：
上板匀速平板间流动（c++）【转载】" title="下板不动, 上板匀速平板间流动（c++）【转载】">
充满着运动黏度系数υ=1cm/s的液体。上板做水平运动并在0.1s时间内，速度线性由0线性地增加到10cm/s，如下图所示：
上板匀速平板间流动（c++）【转载】" title="下板不动, 上板匀速平板间流动（c++）【转载】">
通过对N—S方程的简化，可由下面的抛物线方程来描述
上板匀速平板间流动（c++）【转载】">
流动区域在z=0和z=0.1cm之间，初始条件为u(z,0)=0cm/s，边界条件为：u|=0，u|=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;
}
运行结果：
上板匀速平板间流动（c++）【转载】" title="下板不动, 上板匀速平板间流动（c++）【转载】">
不同时刻的速度变化：
上板匀速平板间流动（c++）【转载】">
生成该图的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
不同高度的速度变化图
上板匀速平板间流动（c++）【转载】">
生成该图的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