SUNTANS模型学习（4）——再学习cylinder算例（定解条件的设定）

在上一个博客中，我做了一个关于cylinder算例的学习笔记（SUNTANS模型学习（3）——学习cylinder算例）；主要涉及了网格和参数。除此之外，模型中还有一些极为重要的要素，如定解条件。
在SUNTANS模型中，定解条件一般包括初始条件与边界条件。上述条件都是通过模型代码给定的。对于Cylinder算例，这些条件包含在initialization.c和boundaries.c文件之中。此外，模型中还有些方程需要给定，比如水体的状态方程（确定水体的密度）；这部分的代码在state.c之中。

初始条件（initialization.c）

在initialization.c里，我们需要设定的条件有如下几项：

1. 设定垂向网格的间距：GetDZ
2. 设定静水深（或底面高程）：ReturnDepth
3. 设定初始水位：ReturnFreeSurface
4. 初始盐度、温度场：ReturnSalinity, ReturnTemperature
5. 初始速度场：ReturnHorizontalVelocity
6. 初始泥沙浓度场：ReturnSediment, ReturnBedSedimentRatio

垂向网格间距

int GetDZ(REAL *dz, REAL depth, REAL localdepth, int Nkmax, int myproc) {
  int k, status;
  REAL z=0, dz0, r = GetValue(DATAFILE,"rstretch",&status);

  if(dz!=NULL) {
    if(r==1) 
      for(k=0;k<Nkmax;k++)
		dz[k]=depth/Nkmax;
    else if(r>1 && r<=1.1) {    
      	dz[0] = depth*(r-1)/(pow(r,Nkmax)-1);
      	if(VERBOSE>2) printf("Minimum vertical grid spacing is %.2f\n",dz[0]);
      	for(k=1;k<Nkmax;k++) 
			dz[k]=r*dz[k-1];
    } else if(r>-1.1 && r<-1) {    
      r=fabs(r);
      dz[Nkmax-1] = depth*(r-1)/(pow(r,Nkmax)-1);
      if(VERBOSE>2) printf("Minimum vertical grid spacing is %.2f\n",dz[Nkmax-1]);
      for(k=Nkmax-2;k>=0;k--) 
		dz[k]=r*dz[k+1];
    } else {
      printf("Error in GetDZ when trying to create vertical grid:\n");
      printf("Absolute value of stretching parameter rstretch must  be in the range (1,1.1).\n");
      exit(1);
    }
  } else {
    r=fabs(r);
    if(r!=1)
      dz0 = depth*(r-1)/(pow(r,Nkmax)-1);
    else
      dz0 = depth/Nkmax;
    z = dz0;
    for(k=1;k<Nkmax;k++) {
      dz0*=r;
      z+=dz0;
      if(z>=localdepth) {
		return k;
      }
    }
  }
}

上述代码的主要功能是确定网格间距dz[ ]，其返回值为该水平空间位置点处最底网格的索引值。代码主体为两个部分：

1. 根据参数r（参数r在输入文件中给定），递推计算网格间距。[第4~23行]
2. 得到最底部网格的索引值，返回作函数值。[第24~35行]

静水深（或底面高程）和初始水位

在本算例中，数值水槽为平底水槽，初始静水深为0.25 m（或水槽底面高程为-0.25 m）。

REAL ReturnDepth(REAL x, REAL y) {
  return 0.25;
}

在初始时刻，全局的水位值均设置为0 m。

REAL ReturnFreeSurface(REAL x, REAL y, REAL d) {
  return 0;
}

初始盐度、温度场、泥沙浓度场

本算例不涉及盐度计算。故初始盐度场设置为零盐度。

REAL ReturnSalinity(REAL x, REAL y, REAL z) {
  return 0;
}

温度场按照如下代码设置：

REAL ReturnTemperature(REAL x, REAL y, REAL z, REAL depth) {
  if(y<.2)
    return 1;
  return 0;
}

上述代码表示全场初始温度分布为：
$$T=\{\begin{array}{ll}1& \text{, x < 0.2}\\ 0& \text{, x >= 0.2}\end{array}$$

除此之外，本算例不涉及泥沙输运模拟。故泥沙浓度均设置为0.

REAL ReturnSediment(REAL x, REAL y, REAL z, int sizeno) {
  return 0;
}
REAL ReturnBedSedimentRatio(REAL x, REAL y, int layer, int sizeno,int nsize) {
   return 0;
}

初始速度场

模型的初始速度场可以通过数据文件来导入（热启动）。对于本算例，初始水平速度通过函数ReturnHorizontalVelocity给出；它在phys.c中被调用，用于初始水平流速的生成。代码如下：

REAL ReturnHorizontalVelocity(REAL x, REAL y, REAL n1, REAL n2, REAL z) {
  return n1;
}

此处n1的值即 grid->n1。在模型中，表示某一网格边（edge）中心点的法向向量为(n1, n2)，n1即为法相向量的x方向分量。

边界条件（boundaries.c）

首先回顾该算例的边界设置。如下图所示，该平面图展示了水槽边界处所采用的条件及各边界在模型中的编号（mark）。这个编号mark指代了边界的类型，它将在网格边界的数据文件edges中体现；例如，若某一网格边恰为水槽的固壁边界，则其采用的边界条件是Closed，其edge的编号为1。

boundaries.c文件中包含了诸多函数。对于本算例，主要涉及的函数有：

1. OpenBoundaryFluxes：计算处于开边界处的网格边的通量；
2. BoundaryScalars：开边界处的温度、盐度边界条件；
3. BoundarySediment：开边界处的泥沙浓度，在本算例中未被设置和调用（此处略过）；
4. BoundaryVelocities：确定边界处的速度、水位；
5. SetUVWH：计算了开边界处网格中心的u, v, w，不过该函数未在主程序中被调用（此处略过）；
6. WindStress：风应力（在本算例中，风应力为0，此处略过）；
7. 其余函数，包括InitBoundaryData、AllocateBoundaryData，这些函数代码均位于主程序中（此处略过）。

开边处的通量计算（OpenBoundaryFluxes）

void OpenBoundaryFluxes(REAL **q, REAL **ub, REAL **ubn, gridT *grid, physT *phys, propT *prop) {
  int j, jptr, ib, k, forced;
  REAL *uboundary = phys->a, **u = phys->uc, **v = phys->vc, **uold = phys->uold, **vold = phys->vold;
  REAL z, c0, c1, C0, C1, dt=prop->dt, u0, u0new, uc0, vc0, uc0old, vc0old, ub0;

  for(jptr=grid->edgedist[2];jptr<grid->edgedist[5];jptr++) {
    j = grid->edgep[jptr];

    ib = grid->grad[2*j];

    for(k=grid->etop[j];k<grid->Nke[j];k++) 
      ub[j][k]=0;

    for(k=grid->etop[j];k<grid->Nke[j];k++) 
      ub[j][k]=phys->boundary_u[jptr-grid->edgedist[2]][k]*grid->n1[j]+
	phys->boundary_v[jptr-grid->edgedist[2]][k]*grid->n2[j];
  }
}

本函数的主体过程可以这么理解，它对mark = 2~4的所有网格边（体现在for循环的循环条件中）都进行了通量计算。我们知道，在SUNTANS模型中，速度的空间离散方式如下：

各速度均定义在网格的界面的中心处，而非网格中心。对于每个网格边，我们都定义了一个法向量 (n1, n2)，如上图中的各个箭头方向。
在上段代码的11 ~ 12行，我们首先指定了各个网格的边界速度ub=0，以进行初始化；在14 ~ 16行，我们通过边界水平速度矢量 (boundary_u, boundary_v) 与法向量 (n1, n2) 的点积求得网格界面的速度。
注意：boundary_u和boundary_v也是在boundaries.c中给定的，关于它的设定详见BoundaryVelocities函数。

开边处的速度、水位（BoundaryVelocities）

void BoundaryVelocities(gridT *grid, physT *phys, propT *prop, int myproc, MPI_Comm comm) {
  int jptr, j, ib, k;
  REAL z, utheta, r, u, v, x0, y0;

  for(jptr=grid->edgedist[2];jptr<grid->edgedist[3];jptr++) {
    j = grid->edgep[jptr];

    ib = grid->grad[2*j];

    if(grid->xv[ib]>1)
      for(k=grid->etop[j];k<grid->Nke[j];k++) {
	// Right boundary
		phys->boundary_u[jptr-grid->edgedist[2]][k]=1;
		phys->boundary_v[jptr-grid->edgedist[2]][k]=phys->vc[ib][k];
		phys->boundary_w[jptr-grid->edgedist[2]][k]=phys->w[ib][k];
      }
    else
      for(k=grid->etop[j];k<grid->Nke[j];k++) {
	// Left boundary
		phys->boundary_u[jptr-grid->edgedist[2]][k]=1;
		phys->boundary_v[jptr-grid->edgedist[2]][k]=0.0;
		phys->boundary_w[jptr-grid->edgedist[2]][k]=0;
      }
  }

  // Add circulation with utheta nonzero to simulate cylinder with lift
  r = 0.1;
  utheta = 0.0;
  x0 = 0.3;
  y0 = 0.15;
  for(jptr=grid->edgedist[4];jptr<grid->edgedist[5];jptr++) {
    j = grid->edgep[jptr];

    ib = grid->grad[2*j];
    u = utheta*(grid->ye[j]-y0)/r;
    v = -utheta*(grid->xe[j]-x0)/r;

    for(k=grid->etop[j];k<grid->Nke[j];k++) {
      phys->boundary_u[jptr-grid->edgedist[2]][k]=u;
      phys->boundary_v[jptr-grid->edgedist[2]][k]=v;
      phys->boundary_w[jptr-grid->edgedist[2]][k]=0;
    }
  }
}

本函数主体分为了两个部分，第一部分（10 ~ 24行）给定了左右边界的速度值，第二部分（27 ~ 43行）设定了柱体边界（mark = 4）的no-slip条件。
对于流出边界（右边界），模型设定x方向速度boundary_u为1；而boundary_v和boundary_w均采用边界网格中心v和w的计算值，这相当于零梯度边界条件。对于流入边界（左边界），boundary_u设为1，而其它方向速度均为0。
对于柱体边界（mark = 4），我们首先设定其边界点处的水平速度为utheta = 0 （第35 ~ 36行），再者将这些速度分别存入boundary_u和boundary_v中；最后，设定boundary_w也等于0。这样一来，柱体边界上的速度为0。

开边界处的温度、盐度（BoundaryScalars）

void BoundaryScalars(gridT *grid, physT *phys, propT *prop, int myproc, MPI_Comm comm) {
  int jptr, j, ib, k;
  REAL z, Tb;

  for(jptr=grid->edgedist[2];jptr<grid->edgedist[3];jptr++) {
      j=grid->edgep[jptr];
      ib=grid->grad[2*j];

      if(grid->yv[ib]<0.2) 
		Tb = 1.0;
      else
		Tb = 0.0;
      for(k=grid->ctop[ib];k<grid->Nk[ib];k++) {
		phys->boundary_T[jptr-grid->edgedist[2]][k]=Tb;
		phys->boundary_s[jptr-grid->edgedist[2]][k]=phys->s[ib][k];
      }
  }
}

本算例不涉及盐度计算，而仅有温度计算。上述代码表示，在流入边界，当y<0.2，边界温度boundary_T被设定为1.0；当y>0.2，boundary_T被设定为0.0。
在模拟中，温度可被视作流入水体的一个被动示踪剂。它不会影响水体密度，也不影响流场，但能对流入水体的运动进行失踪。其模拟结果将在本blog的最后展示。

状态方程（state.c）

在state.c中只有一个函数StateEquation。在主程序中，它被调用作求解水体密度的函数。在SUNTANS中，水体密度通常以相对密度的方式表示，即密度与参考密度（一般为1000 kg/m³）的差值。该函数的主要内容

REAL StateEquation(const propT *prop, const REAL s, const REAL T, const REAL p) {
  return 0;
}

这表明，在本案例中，水体密度为常数，相对密度为0。
此外在其它的一些算例中，这个StateEquation还有以下这样的形式：

REAL StateEquation(const propT *prop, const REAL s, const REAL T, const REAL p) {
  return prop->beta*s;
}

即相对密度和盐度s成正比，beta为比例系数。
所以，我们可以按实际需要，在state.c中给定水体密度的状态方程。

模拟结果中的标量场（温度场）

当 t = 0.4 s，

当 t = 2.0 s,