设计正压梯度力的计算范文

《海洋预报杂志》2016年第3期

摘要：

基于一套自主研制的无结构网格二维河口海洋数值模式A2D，在大圆湖理想模型下，通过与解析解进行比较分析，采用不同架构配置，改进设计正压梯度力计算方法。改进后的算法中引入了从算架构的配置，以配合主算架构，得到更佳的稳定性。通过水位场平面分布与单点过程线可以发现，三组试验的算法均获得了较好的精度和比原算法更好的稳定性，其中TSNS配置算法（中心点计算水位、边中点计算流速的主算架构，配合节点计算水位、边中点计算流速的从算架构）由于其主算架构更接近结构网格下的C网格，在守恒性、移动潮滩边界处理等方面具有一定优势和便利性，有利于在实际海洋中的计算。将TSNS配置算法在江浙沿海进行试算，水位验证结果与实测基本符合，与原A2D模式计算水位之间无显著差异。TSNS算法在稳定性方面的改进，有助于提升模式升级为三维后的稳定性，为今后模式成功升级为三维打下基础。

关键词：

无结构网格；海洋模式；解析解；正压梯度力；稳定性

1引言

无结构网格凭借其易拟合岸界、可局部加密等优势，正越来越多地被用于海洋水动力数值模式。当前国际上知名的模式多来自国外[1-12]，其中采用无结构网格的模式有FVCOM[10-11]、ADCIRC[12]等，这些模式在国内被广泛使用，而国内自主开发的模式[13-14]则非常缺乏。鉴于此，2006—2011年期间，笔者所在研究组自主研制了一套无结构网格二维河口海洋数值模式[15-16]，其中2011年最终版本后文称为A2D[16]。A2D模式采用无结构三角形网格，基于有限体积法求解。水位在网格中心（重心）求解，流速在网格边中点求解（同时求解x和y方向流速）。这种配置结合了各种已有的无结构网格模式的长处。一方面，它与国际上较为通行的C网格[17-18]配置相近，具有在网格中心求解水位之利于干湿判别的优点；另一方面，直接求解x-方向和y-方向流速使得离散简洁高效，无须水平坐标转换，且能更方便地利用有限体积法，提升守恒性。A2D模式在时间上显式求解，采用预估修正法[19-22]。通过理想试验、黄浦江和长江口等环境下的试验计算，有效地验证了A2D模式的精度[16]。但作为新生事物，A2D模式尚不成熟，存在一些不足，如模式还只是二维、在理想试验下有微小数值波动等。笔者所在研究组曾将A2D模式以相同架构和算法升级至三维版本（称之为A3D）[23]，但在对长江口海域的计算中发现A3D的水动力不如A2D稳定，所以三维版本暂未获得成功。A3D的不稳定性，很有可能源自于A2D的架构和算法，其中求解流速时正压梯度力项为主要项，是关注的重点。本文对A2D的正压梯度力算法进行分析与改进，提升模式在一个理想的大圆湖解析解模型下的稳定性表现，并投入实际海洋中试算。 

2大圆湖理想模型简介

Csanady提出了一个理想的大圆湖模型[24]，计算了该圆湖在风的作用下的表明波动。Birchfield指出了其中一些错误并给出了正确的解析解表达式[25]。该模型为一个半径67.5km，水深75m的平底大圆湖（见图1），初始时刻水位静止为0，施加以恒定西风3m/s，科氏力系数f取常数0.0001，忽略底摩擦。在风与科氏力的共同作用下，湖表将会产生水位波动，并且可以得到该波动的解析解表达式。通过这个模型，可以将数值模式的数值解与解析解进行比较，从而探讨模式的精度。由于这个大圆湖的空间尺度很大，根据量纲分析的结果，流速平流项和扩散项对水位和流速变化的贡献极其微小，可忽略不计[24-25]，所以影响模式计算精度的主要项为正压梯度力项，因此A3D在此模型下的稳定性问题，应是缘于正压梯度力项。重新设计模式的算法，使得正压梯度力项得到更合理地求解，并在大圆湖理想模型下取得较好的结果，是本文的主要目标。

3原模式的稳定性分析

采用A2D和A3D分别对大圆湖进行模拟（网格见图1），时间步长均为5s，并在A点输出站位过程线作比较。模式A2D与解析解的过程线比较结果如图2所示。可以看到，模式的计算结果与解析解十分吻合，且总体比较光滑，没有出现明显的不稳定，说明模式的水动力达到了较高的精度。但在4.89d左右，水位的数值解（红色）存在微小的波动，即存在微小的不稳定迹象。而在相同架构和算法的三维模式A3D下，水位的不稳定性表现得更为明显（图3）。对比解析解（见图4）和A2D（见图5）的水位场分布也可以发现，在模式5d时，A2D的计算结果的等值线略显抖动，这也反映出模式在稳定性上略有欠缺。A2D在理想试验下的微小不稳定，在实际海洋中计算时并不会发生，因为无频散的流速平流项会使得模式保持稳定。但当三维开发版本A3D模式在实际海洋中计算时，由于不稳定性增大，所以难免出现个别区域流速或流向异常，影响计算结果的正确性。所以，本文改进A2D正压梯度力项算法，使得模式更为稳定，是模式升级为稳定的三维版本的基础。

4算法改进和探讨

4.1原模式控制方程组A2D模式包含水动力计算和盐度计算，而温度暂时不作计算，取为常数T=10°C。本文不牵涉温盐算法的改进，主要对水动力控制方程组的相关部分作简介。对流体不可压缩、Boussinesq和静力近似下的海洋动力学原始方程组作垂向积分，可得到垂向平均的二维控制方程组：式中：t为时间，ζ表示海表水位，D代表总水深（总水深D=H+ζ，H为固定不变的基准水深），U=1D∫-Hζudz和V=1D∫-Hζvdz分别表示水平x方向和y方向的垂向平均流速（其中u和v分别为空间中某点的水平x方向和y方向的局地流速），Fx和Fy为x方向和y方向的流速水平扩散，AM为水平湍流粘滞系数，f为柯氏力系数，g为重力加速度，<wu（0）>和<wv（0）>为海表应力，<wu（-1）>和<wv（-1）>为海底摩擦应力，ρ为密度，ρ0为参考密度。其中密度ρ=ρ（S,T,p）为盐度S、温度T（°C）和海压p（Pa）的函数，采用1980年国际海水状态方程（EOS80）计算。海表应力<wu（0）>和<wv（0）>为风应力在x和y方向上的分量：式中：ρa为空气密度；Ua和Va分别为x和y方向的风速，||Va=U2a+V2a为风矢量Va的绝对值大小；CD为海水对风的拖曳系数，它根据Large和Pond改进的稳定状态拖曳系数计算[26]：海底摩擦应力<wu（-1）>和<wv（-1）>根据谢才公式得到：式中：nb为曼宁系数，一般取值为0.015到0.018之间。

4.2原模式正压梯度力算法分析在已有的A2D模式中，正压梯度力项的计算方法如图6所示。构造至多4个控制体A1、A2、A3和A4（未必一定有4个，网格边缘和干点附近会有缺失），每个控制体的3个顶点均有计算得到的水位，于是可以通过格林公式计算控制体内的水位梯度。再将各控制体的水位梯度平均，得到边j上的水位梯度，从而计算得到边j上的正压梯度力结果。观察此算法，较为显著的问题为未能利用到节点水位进行正压梯度力计算，控制体远端的网格点距离较远，不利于数值稳定。而如果对水位进行插值，则由于水位梯度计算本来就对空间位置敏感，会导致精度不高，这在A2D研发阶段已做过尝试，其效果不如A2D最终方案理想[16]。在当前配置架构下暂时无法找到更好算法的情况下，尝试更多不同的配置架构是较好的途径。简便起见，将不同配置架构的试验采用其配置特征命名，分别用N、S、T3个字母代表三角形网格的节点、边中点、网格中心，试验的前两位字母代表水位计算点和流速计算点的位置。如原A2D模式的命名为：TS。有些试验采用一种配置架构作为主算架构，另一种配置架构作为从算架构，从算架构作为第2套同时计算的辅助配置架构为主算架构提升稳定性，这些试验的命名中第1—2位字母代表了主算架构的配置，第3—4位字母代表了从算架构的配置，如TSNS即代表网格中心水位、边中点流速的主算架构结合节点水位、边中点流速的从算架构的试验。在经过对多种不同配置算法的尝试后发现，部分配置架构下得到了较稳定的结果，其中包括TSNS。TSNS在节点计算辅助的水位，有助于提升模式的稳定性。这种配置的主算架构仍然是中心 位、边中点流速，但增加计算节点水位、边中点流速的从算架构，使得中心与节点同时具有水位，更有利于边中点位置的正压梯度力求解。

4.3TSNS配置算法和计算结果在TSNS配置的算法中，通过连续方程，利用边中点的流速同时求出网格中心点和网格节点的水位，并通过动量方程，利用网格中心点和网格节点的水位求出边中点的流速。求网格中心点i的水位时，根据连续方程（1），仍然沿用原A2D的求法，在三角形控制体A（图7a）中根据格林公式求得：式中：l为绕A一周的正向曲线，即l的方向为逆时针，A始终在其左边。求网格节点m的水位时，也以格林公式（13）求解，不过控制体A变为包围节点m的多边形（图7b）。控制体A有可能完全包围节点m，也有可能如图7b一般缺失部分角度，无论哪种情形，l均为绕A一周的正向曲线。在图7b的情形下，仅需计算j4-i3-j3-i2-j2-i1-j1上的线积分，因为j1-m-j4上线积分等于0无通量进出。每一小段的流速为该小段所连接边中点的流速。求边中点j的流速时，对正压梯度力项的算法作了修改。在大圆湖模型中，无斜压梯度力项，平流项也非常小可忽略不计，故对精度影响最大的项正是正压梯度力项gD∂ζ∂x和gD∂ζ∂y。TSNS配置算法中，动量方程仅正压梯度力项较原A2D模式有改变，其它项均维持不变。TSNS配置算法在计算边中点j的流速时，通过格林公式解出正压梯度力项中的∂ζ∂x和∂ζ∂y，其中控制体A取为连接节点和中心点的菱形区域（见图8），l为包围A的正向曲线。由于节点和中心点的水位值模式均有计算，故沿i2-m2-i1-m1-i2的线积分易于求得。由于在以上计算方法中，同时计算了网格中心与网格节点两套水位，在时间积分久后两套的数值解可能存在分离的现象，所以将从算架构中的节点水位按照一定速度向主算架构中的中心点水位回归。这里取每时间步0.05的回归比例进行趋近。将上述TSNS算法在大圆湖模型中试验，时间步长仍取5s，发现TSNS的精度与原A2D接近，稳定性更好（见图9），第4.89d未出现不稳定抖动，水位场等值线也更平整。

4.4其它配置的计算结果除了TSNS配置外，本文还测试了一些其它配置，其中NSTS配置和ST配置这两组试验也得到了不错的结果。在NSTS配置中，主算与从算架构与TSNS对换，NS主算、TS从算，回归系数仍取0.05。其结果与TSNS略有不同（见图10），精度接近TSNS，稳定性较好。而ST配置下也得到了稳定性尚可的结果（见图11）。通过统计点A水位时间序列的平均误差和均方根误差来比较这3组试验与原A2D模式的精度，可以发现TSNS配置的精度在所有4种配置中最好（见表1），同时误差在模式4-5d时较之前时间段有所增大。由于模式最终应用时，以TS作为主算架构的TSNS配置更接近结构网格下的C网格，在守恒性、移动潮滩边界处理等方面具有一定优势和便利性，而NSTS配置和ST配置在各种边界条件设计中存在一定的困难，故NSTS配置和ST配置不作重点介绍，具体算法细节在此省略，仅展示理想试验下的结果，后续在实际海洋中的试算仅基于TSNS配置下进行。  

5江浙沿海海域试算

5.1模式网格和基本设置基于改进了正压梯度力算法的TSNS配置模式和原A2D模式，对东海区范围内包含吕四测站、嵊山测站和定海测站的江浙沿海海域进行试算，以观察模拟效果，并进行水位验证，其中TSNS配置模式因配置变化对边界条件进行了部分完善。模式的网格范围见图12，基于54坐标，包括了整个长江口、杭州湾和邻近海区。东边至124.5°E附近，北边到34.3°N左右，南边到28.4°N左右，长江上游边界取在大通。长江口内、深水航道附近和岛屿附近的网格作了局部加密，最小网格分辨率可达100m，而口外网格则被放大，最大超过10000m。模式中深水航道的导堤和丁坝涨潮时淹没、落潮时露出，是作为动边界处理的。模式的时间步长统一取为1s。外海开边界处的水位利用16个分潮（M2，S2，N2，K2，K1，O1，P1，Q1，MU2，NU2，T2，L2，2N2，J1，M1，OO1）的调和常数计算得出，而在长江口上游大通处则利用实测径流量资料给出通量边界条件。海表面的风场以每6h为1组、分辨率为0.5°×0.5°经纬度的气象预报后处理结果给出底摩擦曼宁系数在模式中设为0.015。模式从2008年11月5日起计算，共计算30d。通过搜集的吕四站、嵊山站和定海站的实测水位资料，对原A2D模式和改进后的TSNS配置模式进行比对验证。

5.2水位验证结果模式共运行30d，输出第10—30d的水位，与吕四站、嵊山站和定海站的实测资料进行比对（图13—15，图中上子图的蓝实线为原A2D计算水位，黑虚线为TSNS配置计算水位，红点为实测水位；下子图的黑实线为TSNS计算水位减去A2D计算水位的差，蓝实线为0轴）。可以发现，TSNS配置模式的水位计算结果与原A2D模式（TS配置）的水位计算结果非常接近，两者差异在3站均不到0.1m，同时模式计算结果与实测基本符合，基本反映了3个测站的水位变化规律，其中嵊山站和定海站存在较明显的潮汐日不等现象，而定海站的振幅略偏大。由于此番改进主要目的是探讨是模式本身的性能，故在江浙沿海海域进行试算并验证时，底摩擦曼宁系数统一取了0.015，未针对不同海域进行分区设置，对验证效果略有影响。而TSNS配置模式与原A2D模式虽然在算法上存在差异，但由于案例设定的物理环境相同，两者受相同的边界条件驱动，故两者的水位计算结果与变化特征并没有出现显著差异。

6结论

本文通过对模式配置架构的改变，重新设计正压梯度力项的算法，得到了3种比原二维A2D模式稳定性更高的算法，在大圆湖理想模型下模拟得到了较稳定的结果，同时精度与原A2D接近，均与解析解较为符合。其中TSNS配置算法中由于其主算架构TS配置更接近结构网格下的C网格，在守恒性、移动潮滩边界处理等方面较好，故在完善了相应的边界条件后，在江浙沿海海域进行了试算并与实测资料进行了对比验证。水位验证结果与实测基本符合，同时与原A2D模式计算水位之间无显著 差异。原A2D的算法虽然与TSNS算法一样也能在二维情形的实际海洋中稳定，且计算结果相近，但A2D的三维开发版本A3D在实际海洋计算中却遭遇稳定性不佳的问题。大圆湖理想模型下更稳定的表现，有助于提升模式升级为三维后的稳定性，所以TSNS算法在稳定性方面的改进，为今后模式成功升级为三维打下基础。经过改进后的TSNS配置模式可在风暴潮模拟等场合进行应用，同时由于其网格配置特性，可进一步优化边界条件设定，根据需要灵活设置堤坝等特殊情形，具有一定的发展前景。

作者：陈昞睿 单位：国家海洋局东海预报中心