首页 > 代码库 > 柯朗数(Courant number)研究

柯朗数(Courant number)研究

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

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 interface113                 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 cell123                 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 faces134       {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。

 

柯朗数(Courant number)研究