无结构网格二维海洋模式的正压梯度力算法改进
2016-03-18陈睿
陈睿
（国家海洋局东海预报中心,上海201206）
无结构网格二维海洋模式的正压梯度力算法改进
（国家海洋局东海预报中心,上海201206）
摘要：基于一套自主研制的无结构网格二维河口海洋数值模式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.5 km，水深75 m的平底大圆湖（见图1），初始时刻水位静止为0，施加以恒定西风3 m/s，科氏力系数f取常数0.0001，忽略底摩擦。在风与科氏力的共同作用下，湖表将会产生水位波动，并且可以得到该波动的解析解表达式。通过这个模型，可以将数值模式的数值解与解析解进行比较，从而探讨模式的精度。
图1 大圆湖网格及输出点A位置
由于这个大圆湖的空间尺度很大，根据量纲分析的结果，流速平流项和扩散项对水位和流速变化的贡献极其微小，可忽略不计[24-25]，所以影响模式计算精度的主要项为正压梯度力项，因此A3D在此模型下的稳定性问题，应是缘于正压梯度力项。重新设计模式的算法，使得正压梯度力项得到更合理地求解，并在大圆湖理想模型下取得较好的结果，是本文的主要目标。
图2 A2D模式的数值解（红色）与解析解（蓝色）的过程线比较（U是东向流速，D是水深）
图3 A3D模式的数值解（黑色）与解析解（红色）的水位过程线比较
图4 解析解在1 h、1 d、5 d时的水位场分布
图5 A2D模式在1 h、1 d、5 d时的水位场分布
3 原模式的稳定性分析
采用A2D和A3D分别对大圆湖进行模拟（网格见图1），时间步长均为5 s，并在A点输出站位过程线作比较。模式A2D与解析解的过程线比较结果如图2所示。可以看到，模式的计算结果与解析解十分吻合，且总体比较光滑，没有出现明显的不稳定，说明模式的水动力达到了较高的精度。但在4.89 d左右，水位的数值解（红色）存在微小的波动，即存在微小的不稳定迹象。而在相同架构和算法的三维模式A3D下，水位的不稳定性表现得更为明显（图3）。对比解析解（见图4）和A2D（见图5）的水位场分布也可以发现，在模式5d时，A2D的计算结果的等值线略显抖动，这也反映出模式在稳定性上略有欠缺。
A2D在理想试验下的微小不稳定，在实际海洋中计算时并不会发生，因为无频散的流速平流项会使得模式保持稳定。但当三维开发版本A3D模式在实际海洋中计算时，由于不稳定性增大，所以难免出现个别区域流速或流向异常，影响计算结果的正确性。所以，本文改进A2D正压梯度力项算法，使得模式更为稳定，是模式升级为稳定的三维版本的基础。
4 算法改进和探讨
4.1原模式控制方程组
A2D模式包含水动力计算和盐度计算，而温度暂时不作计算，取为常数T=10°C。本文不牵涉温盐算法的改进，主要对水动力控制方程组的相关部分作简介。对流体不可压缩、Boussinesq和静力近似下的海洋动力学原始方程组作垂向积分，可得到垂向平均的二维控制方程组：
式中：t为时间，ζ表示海表水位，D代表总水深（总水深D=H+ζ，H为固定不变的基准水深），分别表示水平x方向和y方向的垂向平均流速（其中u和v分别为空间中某点的水平x方向和y方向的局地流速），Fx和Fy为x方向和y方向的流速水平扩散，AM为水平湍流粘滞系数，f为柯氏力系数，g为重力加速度，
海表应力
式中：ρa为空气密度；Ua和Va分别为x和y方向的风速，为风矢量的绝对值大小；CD为海水对风的拖曳系数，它根据Large和Pond改进的稳定状态拖曳系数计算[26]：
海底摩擦应力
式中：nb为曼宁系数，一般取值为0.015到0.018之间。
4.2原模式正压梯度力算法分析
在已有的A2D模式中，正压梯度力项的计算方法如图6所示。构造至多4个控制体A1、A2、A3和A4（未必一定有4个，网格边缘和干点附近会有缺失），每个控制体的3个顶点均有计算得到的水位，于是可以通过格林公式计算控制体内的水位梯度。再将各控制体的水位梯度平均，得到边j上的水位梯度，从而计算得到边j上的正压梯度力结果。
观察此算法，较为显著的问题为未能利用到节点水位进行正压梯度力计算，控制体远端的网格点距离较远，不利于数值稳定。而如果对水位进行插值，则由于水位梯度计算本来就对空间位置敏感，会导致精度不高，这在A2D研发阶段已做过尝试，其效果不如A2D最终方案理想[16]。
图6 A2D正压梯度力算法示意图，浅蓝色区域为控制体A1，至多有4个控制体
在当前配置架构下暂时无法找到更好算法的情况下，尝试更多不同的配置架构是较好的途径。简便起见，将不同配置架构的试验采用其配置特征命名，分别用N、S、T 3个字母代表三角形网格的节点、边中点、网格中心，试验的前两位字母代表水位计算点和流速计算点的位置。如原A2D模式的命名为：TS。有些试验采用一种配置架构作为主算架构，另一种配置架构作为从算架构，从算架构作为第2套同时计算的辅助配置架构为主算架构提升稳定性，这些试验的命名中第1—2位字母代表了主算架构的配置，第3—4位字母代表了从算架构的配置，如TSNS即代表网格中心水位、边中点流速的主算架构结合节点水位、边中点流速的从算架构的试验。在经过对多种不同配置算法的尝试后发现，部分配置架构下得到了较稳定的结果，其中包括TSNS。TSNS在节点计算辅助的水位，有助于提升模式的稳定性。这种配置的主算架构仍然是中心水位、边中点流速，但增加计算节点水位、边中点流速的从算架构，使得中心与节点同时具有水位，更有利于边中点位置的正压梯度力求解。
图7 求解水位的控制体模型
4.3 TSNS配置算法和计算结果
在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的流速时，对正压梯度力项的算法作了修改。在大圆湖模型中，无斜压梯度力项，平流项也非常小可忽略不计，故对精度影响最大的项正是正压梯度力项。TSNS配置算法中，动量方程仅正压梯度力项较原A2D模式有改变，其它项均维持不变。TSNS配置算法在计算边中点j的流速时，通过格林公式
图8 求解正压梯度力项的控制体
由于在以上计算方法中，同时计算了网格中心与网格节点两套水位，在时间积分久后两套的数值解可能存在分离的现象，所以将从算架构中的节点水位按照一定速度向主算架构中的中心点水位回归。这里取每时间步0.05的回归比例进行趋近。
图9 TSNS配置结果分析（U是东向流速，D是水深，下同）
将上述TSNS算法在大圆湖模型中试验，时间步长仍取5 s，发现TSNS的精度与原A2D接近，稳定性更好（见图9），第4.89 d未出现不稳定抖动，水位场等值线也更平整。
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配置下进行。
图10 NSTS配置结果分析
图11 ST配置结果分析
5 江浙沿海海域试算
5.1模式网格和基本设置
基于改进了正压梯度力算法的TSNS配置模式和原A2D模式，对东海区范围内包含吕四测站、嵊山测站和定海测站的江浙沿海海域进行试算，以观察模拟效果，并进行水位验证，其中TSNS配置模式因配置变化对边界条件进行了部分完善。模式的网格范围见图12，基于54坐标，包括了整个长江口、杭州湾和邻近海区。东边至124.5°E附近，北边到34.3°N左右，南边到28.4°N左右，长江上游边界取在大通。长江口内、深水航道附近和岛屿附近的网格作了局部加密，最小网格分辨率可达100 m，而口外网格则被放大，最大超过10 000 m。模式中深水航道的导堤和丁坝涨潮时淹没、落潮时露出，是作为动边界处理的。
模式的时间步长统一取为1 s。外海开边界处的水位利用16个分潮（M2，S2，N2，K2，K1，O1，P1，Q1，MU2，NU2，T2，L2，2N2，J1，M1，OO1）的调和常数计算得出，而在长江口上游大通处则利用实测径流量资料给出通量边界条件。海表面的风场以每6 h为1组、分辨率为0.5°×0.5°经纬度的气象预报后处理结果给出（可从网址http ://dss.ucar.edu/datasets/ ds744.4/data/处下载）。底摩擦曼宁系数在模式中设为0.015。模式从2008年11月5日起计算，共计算30 d。通过搜集的吕四站、嵊山站和定海站的实测水位资料，对原A2D模式和改进后的TSNS配置模式进行比对验证。
表1 大圆湖模型下各模式在点A的水位误差统计表（单位：mm）
5.2水位验证结果
模式共运行30 d，输出第10—30 d的水位，与吕四站、嵊山站和定海站的实测资料进行比对（图13—15，图中上子图的蓝实线为原A2D计算水位，黑虚线为TSNS配置计算水位，红点为实测水位；下子图的黑实线为TSNS计算水位减去A2D计算水位的差，蓝实线为0轴）。可以发现，TSNS配置模式的水位计算结果与原A2D模式（TS配置）的水位计算结果非常接近，两者差异在3站均不到0.1 m，同时模式计算结果与实测基本符合，基本反映了3个测站的水位变化规律，其中嵊山站和定海站存在较明显的潮汐日不等现象，而定海站的振幅略偏大。
图12 江浙沿海海域网格
图13 吕四站TSNS配置与原A2D结果对比
由于此番改进主要目的是探讨是模式本身的性能，故在江浙沿海海域进行试算并验证时，底摩擦曼宁系数统一取了0.015，未针对不同海域进行分区设置，对验证效果略有影响。而TSNS配置模式与原A2D模式虽然在算法上存在差异，但由于案例设定的物理环境相同，两者受相同的边界条件驱动，故两者的水位计算结果与变化特征并没有出现显著差异。
图14 嵊山站TSNS配置与原A2D结果对比
图15 定海站TSNS配置与原A2D结果对比
6 结论
本文通过对模式配置架构的改变，重新设计正压梯度力项的算法，得到了3种比原二维A2D模式稳定性更高的算法，在大圆湖理想模型下模拟得到了较稳定的结果，同时精度与原A2D接近，均与解析解较为符合。其中TSNS配置算法中由于其主算架构TS配置更接近结构网格下的C网格，在守恒性、移动潮滩边界处理等方面较好，故在完善了相应的边界条件后，在江浙沿海海域进行了试算并与实测资料进行了对比验证。水位验证结果与实测基本符合，同时与原A2D模式计算水位之间无显著差异。
原A2D的算法虽然与TSNS算法一样也能在二维情形的实际海洋中稳定，且计算结果相近，但A2D的三维开发版本A3D在实际海洋计算中却遭遇稳定性不佳的问题。大圆湖理想模型下更稳定的表现，有助于提升模式升级为三维后的稳定性，所以TSNS算法在稳定性方面的改进，为今后模式成功升级为三维打下基础。经过改进后的TSNS配置模式可在风暴潮模拟等场合进行应用，同时由于其网格配置特性，可进一步优化边界条件设定，根据需要灵活设置堤坝等特殊情形，具有一定的发展前景。
参考文献：
[1] Blumberg A F. A Primer for ECOM-si[R]. Technical Report of HydroQual, Mahwah, New Jersey, 1994.
[2] Blumberg A F, Mellor G L. A Description of a Three-Dimensional Coastal Ocean Circulation Model[M]. Washington, DC: American Geophysical Union, 1987.
[3] Casulli V. Semi-implicit Finite Difference Methods for the Two-dimensional Shallow Water Equations[J]. Journal of Computational Physics, 1990, 86(1): 56-74.
[4] Casulli V. A Semi-implicit Finite Difference Method for Non-hydrostatic, Free-surface Flows[J]. International Journal for Numerical Methods in Fluids, 1999, 30(4): 425-440.
[5] Casulli V, Walters R A. An Unstructured Grid, Three-dimensional Model based on the Shallow Water Equations[J]. International Journal for Numerical Methods in Fluids, 2000, 32(3): 331-348.
[6] Ham D A, Kramer S C, Stelling G S, et al. The Symmetry and Stability of Unstructured Mesh C-grid Shallow Water Models Under the Influence of Coriolis[J]. Ocean Modelling, 2007, 16 (1-2): 47-60.
[7] Zhang Y L, Baptista A M, Myers III E P. A Cross-scale Model for 3D Baroclinic Circulation in Estuary-plume-shelf Systems: I. Formulation and Skill Assessment[J]. Continental Shelf Research, 2004, 24(18): 2187-2214.
[8] Fringer O B, Gerritsen M, Street R L. An Unstructured-grid, Finite-volume, Nonhydrostatic, Parallel Coastal Ocean Simulator [J]. Ocean Modelling, 2006, 14(3-4): 139-173.
[9] DHI. MIKE 21 & MIKE 3 Flow Model FM Hydrodynamic and Transport Module Scientific Documentation[R]. Copenhagen: DHI, 2011.
[10] Chen C S, Liu H D, Beardsley R C. An Unstructured Grid, Finite-volume, Three-dimensional, Primitive Equations Ocean Model: Application to Coastal Ocean and Estuaries[J]. Journal of Atmospheric and Oceanic Technology, 2003, 20(1): 159-186.
[11] Chen C S, Beardsley R C, Cowles G. An Unstructured Grid, Finite-volume Coastal Ocean Model: FVCOM User Manual[M]. New Bedford, Mass: SMAST/UMASSD, 2006.
[12] Luettich R, Westerink J. Formulation and Numerical Implementation of the 2D/3D ADCIRC Finite Element Model Version 44.XX[R]. 2004.
[13]顾峰峰,倪汉根,戚定满,等. AUSM格式在二维浅水方程求解中的应用[J].水科学进展, 2008, 19(5): 624-629.
[14]赵棣华,戚晨,庾维德,等.平面二维水流-水质有限体积法及黎曼近似解模型[J].水科学进展, 2000, 11(4): 368-374.
[15]陈昞睿,朱建荣,吴辉,等.无结构网格二维河口海岸水动力数值模式的建立及其应用[J].海洋学报, 2010, 32(2): 31-39.
[16]陈昞睿.无结构网格二维河口海洋数值模式的研制[D].上海:华东师范大学, 2011.
[17] Arakawa A, Lamb V R. Computational Design of the Basic Dynamical Processes of the UCLA General Circulation Model[J]. Methods of Computational Physics, 1977, 17: 173-265.
[18] Kleptsova O, Pietrzak J D, Stelling G S. On the Accurate and Stable Reconstruction of Tangential Velocities in C-grid Ocean Models[J]. Ocean Modelling, 2009, 28(1-3): 118-126.
[19]冯士筰,李凤岐,李少菁.海洋科学导论[M].北京:高等教育出版社, 1999.
[20]朱建荣.海洋数值计算方法和数值模式[M].北京:海洋出版社, 2003.
[21]朱建荣,杨陇慧,朱首贤.预估修正法对河口海岸海洋模式稳定性的提高[J].海洋与湖沼, 2002, 33(1): 15-22.
[22]朱建荣,朱首贤. ECOM模式的改进及在长江河口、杭州湾及邻近海区的应用[J].海洋与湖沼, 2003, 34(4): 364-374.
[23] Chen B R, Zhu J R, Li L. Accelerating 3D Ocean Model Development by using GPU Computing[M]//Deng W. Future Control and Automation: Proceedings of the 2nd International Conference on Future Control and Automation. Berlin Heidelberg: Springer, 2012, 1: 37-43.
[24] Csanady G T. Motions in a Model Great Lake Due to a Suddenly Imposed Wind[J]. Journal of Geophysical Research, 1968, 73(20): 6435-6447.
[25] Birchfield G E. Response of A Circular Model Great Lake to A Suddenly Imposed Wind Stress[J]. Journal of Geophysical Research, 1969, 74(23): 5547-5554.
[26] Large W G, Pond S. Open Ocean Momentum Flux Measurements in Moderate to Strong Winds [J]. Journal of Physical Oceanography, 1981, 11(3): 324-406.
Algorithm Improvement of Barotropic Force in an Unstructured Grid Two-dimensional Ocean Model
CHEN Bing-rui
（East China Sea Marine Forecasting Center, State Oceanic Administration, Shanghai 201206 China）
Abstract：Based on A2D, an independently developed unstructured grid two-dimensional ocean model, and by comparison with analytical solution under an ideal Model Great Lake, the algorithm of barotropic force was improved via employing different computational designs. In the improved algorithm, an assistant design was introduced to cooperate with the major design in order to acquire better stability. Results of elevation field and site time series showed that algorithms of three experiments got satisfactory accuracy and better stability. The TSNS algorithm was among the three, in which the major design solves elevation at centroid and velocity at mid-point of side, and the assistant design solves elevation at node and velocity at mid-point of side. Due to the similarity of its major design to C-grid design in a structured grid, the TSNS algorithm had advantages in conservation and movable tide-flat boundary treatments, which made the algorithm easier to apply for real ocean simulations. The TSNS algorithm was applied to the simulation in real sea near Jiangsu and Zhejiang. The simulated elevation had a good agreement with observed data, and was similar with results from the original A2D model. The improvement in stability will help TSNS algorithm get better stability in three-dimensional upgraded version, and will be the foundation of a successful three-dimensional version in future.
Key words：unstructured grid; ocean model; analytical solution; barotropic force; stability
作者简介：陈昞睿（1980-），男，工程师，博士，从事海洋数值预报研究。E-mail：cu238@163.com
基金项目：国家海洋局青年海洋科学基金（2013202）。
收稿日期：2015-07-20
中图分类号：P731
文献标识码：A
文章编号：1003-0239（2016）01-0027-10
