柯朗数(Courant number)研究

Posted 2020-08-12 历久弥坚0820

在数值计算过程中，对于计算结果的准确性和效率有很高的要求，但是这两者之间往往互相矛盾；而使用柯朗数可用于平衡两者。

1、柯朗数的定义：

 C = sqrt(gh)*t/s

其中，t是时间步长，s是网格在水平方向的间距。

柯朗数的意义在于表示了在单位时间步长中，有多少个网格的信息发生了移动。经过正确的调整，可以更好地加速收敛和增强解的稳定性。

 

2、C语言实现柯朗数计算：

依据上述方程，在实际计算中采用C语言实现计算固液界面上的柯朗数，结果如下：

  1 void localCourantNumber()
  2 {
  3 
  4 
  5     double rhoe,rhon,rhot;    
  6 
  7     for(i=2;i<=nxm-1;i++)    //Calculation of local Courant number only at internal faces
  8       {
  9         ieast     = i + 1;
 10         dxpe = xc[ieast] - xc[i];
 11         fxe = (xf[i]-xc[i])/dxpe;        
 12         fxp = 1.0 - fxe;
 13 
 14             for(j=2;j<=nym;j++)
 15         {
 16             jnorth    = j + 1;
 17             dypn = yc[jnorth] - yc[j];
 18             fyn = (yf[j] - yc[j])/dypn;    
 19                 fyp = 1.0 - fyn;
 20     
 21               for(k=2;k<=nzm;k++)
 22             {
 23                 ktop      = k + 1;
 24                 dzpt = zc[ktop]-zc[k];
 25                 fzt = (zf[k] - zc[k])/dzpt;
 26                 fzp = 1.0 - fzt;
 27 
 28             
 29                 
 30                 //Calculating density at cell interface
 31                 rhoe = fxp * rho[i][j][k] +   fxe * rho[ieast][j][k];
 32                 
 33 /*                rhoe = 2.0 * rho[i][j][k] * rho[ieast][j][k]/( rho[i][j][k] + rho[ieast][j][k]);*/
 34 
 35                 s = (yf[j]-yf[j-1])*(zf[k]-zf[k-1]);
 36                       vole = dxpe * s;
 37 
 38 
 39                  //Sum of courant numbers of outflow faces of donor cell
 40                 Ce[i][j][k] = fabs(Fe[i][j][k]/(rhoe*vole))*dt;
 41 
 42 /*                printf("Ce=%e\n",Ce[i][j][k]);*/
 43             }
 44         }
 45     }
 46     
 47     
 48         for(i=2;i<=nxm;i++)    //Calculation of local Courant number only at internal faces
 49       {
 50         ieast     = i + 1;
 51         dxpe = xc[ieast] - xc[i];
 52         fxe = (xf[i]-xc[i])/dxpe;        
 53         fxp = 1.0 - fxe;
 54 
 55             for(j=2;j<=nym-1;j++)
 56         {
 57             jnorth    = j + 1;
 58             dypn = yc[jnorth] - yc[j];
 59             fyn = (yf[j] - yc[j])/dypn;    
 60                 fyp = 1.0 - fyn;
 61     
 62               for(k=2;k<=nzm;k++)
 63             {
 64                 ktop      = k + 1;
 65                 dzpt = zc[ktop]-zc[k];
 66                 fzt = (zf[k] - zc[k])/dzpt;
 67                 fzp = 1.0 - fzt;
 68 
 69             
 70                 
 71                 //Calculating density at cell interface
 72                 rhon = fyp * rho[i][j][k] +   fyn * rho[i][jnorth][k];
 73                 
 74 /*                rhon = 2.0 * rho[i][j][k] * rho[i][jnorth][k]/( rho[i][j][k] + rho[i][jnorth][k]);*/
 75 
 76 
 77                 s = (xf[i]-xf[i-1])*(zf[k]-zf[k-1]);
 78                 voln = s * dypn;
 79 
 80 
 81                  //Sum of courant numbers of outflow faces of donor cell
 82                 Cn[i][j][k] = fabs(Fn[i][j][k]/(rhon*voln))*dt;
 83                 
 84 /*                printf("Ce=%e\n",Ce[i][j][k]);*/
 85             }
 86         }
 87     }
 88     
 89     for(i=2;i<=nxm;i++)    //Calculation of local Courant number only at internal faces
 90       {
 91         ieast     = i + 1;
 92         dxpe = xc[ieast] - xc[i];
 93         fxe = (xf[i]-xc[i])/dxpe;        
 94         fxp = 1.0 - fxe;
 95 
 96             for(j=2;j<=nym;j++)
 97         {
 98             jnorth    = j + 1;
 99             dypn = yc[jnorth] - yc[j];
100             fyn = (yf[j] - yc[j])/dypn;    
101                 fyp = 1.0 - fyn;
102     
103               for(k=2;k<=nzm-1;k++)
104             {
105                 ktop      = k + 1;
106                 dzpt = zc[ktop]-zc[k];
107                 fzt = (zf[k] - zc[k])/dzpt;
108                 fzp = 1.0 - fzt;
109 
110             
111                 
112                 //Calculating density at cell interface
113                 rhot = fzp * rho[i][j][k] +   fzt * rho[i][j][ktop];
114                 
115 /*                rhot = 2.0 * rho[i][j][k] * rho[i][j][ktop]/( rho[i][j][k] + rho[i][j][ktop]);*/
116 
117 
118                 s = (xf[i]-xf[i-1])*(yf[j]-yf[j-1]);
119                 volt = s * dzpt;
120 
121 
122                  //Sum of courant numbers of outflow faces of donor cell
123                 Ct[i][j][k] = fabs(Ft[i][j][k]/(rhot*volt))*dt;
124 
125                 
126                 
127 /*                printf("Ce=%e\n",Ce[i][j][k]);*/
128             }
129         }
130     }
131     
132     
133     for(i=2;i<=nxm;i++)    //Calculation of local Courant number only at internal faces
134       {
135 
136             for(j=2;j<=nym;j++)
137         {
138     
139               for(k=2;k<=nzm;k++)
140             {
141             
142                 COutD[i][j][k] = Ce[i][j][k] + Cn[i][j][k] + Ct[i][j][k];
143 /*                printf("COutD=%lf\n",COutD[i][j][k]);*/
144 /*                printf("Ce=%e\n",Ce[i][j][k]);*/
145 /*                printf("Cn=%e\n",Cn[i][j][k]);*/
146 /*                printf("Ct=%e\n",Ct[i][j][k]);*/
147             }
148         }
149     }
150     
151     
152 } 

3、柯朗数使用的注意事项：

在fluent中，用courant number 来调节计算的稳定性与收敛性。一般来说，随着courantnumber 的从小到大的变化，收敛速度逐渐加快，但是稳定性逐渐降低。所以具体的问题，在计算的过程中，最好是把Courant number 从小开始设置，看看迭代残差的收敛情况，如果收敛速度较慢而且比较稳定的话，可以适当的增加courant number 的大小，根据自己具体的问题，找出一个比较合适的courant number，让收敛速度能够足够的快，而且能够保持它的稳定性。

 

Generally, in the explicit schemes of differential method, Courant number is an important number which should be less than 1 in order to assure the stability. However, if the Courant number is too small, much computation time will be lost. So the Courant number could be one of those important parameters affecting computation cost and stability. we could use Courant number to control the time step in the transient simulation in CFD codes. Here is some configuration parameters which could be used in simulation with OpenFOAM。