下板不动, 上板匀速平板间流动（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)的流速
　　求解方程的程序代码：

1. 　　#include
   
2. 　　#include
   
3. 　　#include
   
4. 　　using namespace std;
   
5. 　　int main()
   
6. 　　{
   
7. 　　float u[100],u0[100];
   
8. 　　float b,t,dz,dt,dif,difmax,temp;
   
9. 　　int imax,imax1,iter,i,n;
   
10. 　　dz=0.1;
    
11. 　　dt=0.001;
    
12. 　　imax=11;
    
13. 　　for(i=1;i<=imax;i++)
    
14. 　　{
    
15. 　　u[i]=0;
    
16. 　　u0[i]=0;
    
17. 　　}
    
18. 　　imax1=imax-1;
    
19. 　　n=0;
    
20. 　　t=0;
    
21. 　　b=1.0/(1.0/dt+1.0/dz/dz);
    
22. 　　cout<
    
23. 　　cout<
    
24. 　　<
    
25. 　　<
    
26. 　　<
    
27. 　　<
    
28. 　　<
    
29. 　　<
    
30. 　　<
    
31. 　　<
    
32. 　　<
    
33. 　　<
    
34. 　　<
    
35. 　　<
    
36. 　　do
    
37. 　　{
    
38. 　　t+=dt;
    
39. 　　n+=1;
    
40. 　　if(t<0.1)
    
41. 　　u[imax]=t*100;
    
42. 　　else
    
43. 　　u[imax]=10;
    
44. 　　iter=0;
    
45. 　　do
    
46. 　　{
    
47. 　　difmax=0;
    
48. 　　iter+=1;
    
49. 　　if(iter>100)
    
50. 　　exit(1);
    
51. 　　for(i=2;i<=imax1;i++)
    
52. 　　{
    
53. 　　temp=u[i];
    
54. 　　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]);
    
55. 　　dif=fabs(temp-u[i]);
    
56. 　　if(dif>difmax)
    
57. 　　difmax=dif;
    
58. 　　}
    
59. 　　}while(difmax>0.00001);
    
60. 　　for(i=1;i<=imax;i++)
    
61. 　　u0[i]=u[i];
    
62. 　　if(n0==0)
    
63. 　　{
    
64. 　　cout<
    
65. 　　for(i=1;i<=imax;i++)
    
66. 　　cout<
    
67. 　　cout<
    
68. 　　}
    
69. 　　}while(n<1000);
    
70. 　　cout<
    
71. 　　return 0;
    
72. 　　}

复制代码

　　运行结果：






转自：http://blog.sina.com.cn/s/blog_14d64daa10102wkq1.html