CN111046542B - 一种评估二十面体的le网格计算特性方法及离散方法 - Google Patents

Abstract

本发明公开了一种评估二十面体的LE网格计算特性方法及离散方法，包含：(1)针对用LE跳点网格离散，定义二十面体三角形LE跳点网格中的要素；(2)基于LE网格进行离散，并给出离散算子，针对三角形网格分别按照散度、旋度、导数和切向导数的定义，给出其任一矢量场的散度表达式、其任一矢量场F的旋度表达式、其任一标量场g的法向导数表达式、其切向导数；(3)将动量方程分解为分量形式，并将步骤(2)中的离散算子代入分量形式中，得到相应的离散形式；同样将步骤(2)给出的离散算子代入连续方程，得到其空间离散形式。本发明的方法的离散格式比现在常用的二十面体六边形C网格上的格式更加保形，误差更小，质量和总能量守恒型更好。

Description

一种评估二十面体的LE网格计算特性方法及离散方法

技术领域

本发明涉及一种离散方法，具体涉及一种评估二十面体的LE网格计算特性方法及离散方法。

背景技术

球面网格用于全球大气数值模拟的动力框架设计，目前全球普遍采用普通经纬网格。但是，普通经纬网格存在极点奇异和经线在高纬度聚集的问题，采用有限差分方案会受到CFL稳定性条件限制，使得数值积分时间步长缩小、计算费时。因此，为了避免常规经纬网格存在的数值计算问题，众多科学家尝试构建准均匀球面网格来设计全球大气数值模式。

常见的准均匀球面网格包括：Kurihara(1965)最早提出的“精简网格”， Sadourny等(1968)提出的“正二十面体网格”，Sadourny(1972)提出的“立方体球面网格”，以及Kageyama等(2004)提出一种重叠网格系统(阴阳网格)，以下对各个准均匀球面网格进行介绍：

(1)精简网格，虽然从赤道向两极方向每一个纬度都减少一定的格点，可以节约经纬网格两极附近的计算量，并且使网格距较为均匀。但是，网格递减率会受到计算方案和极点的限制，纬向守恒计算比较困难，同时也会产生位相误差，造成波动的经向坡度加大，从而产生虚假的经向输送；

(2)立方体球面网格，其虽没有极点奇异，但球面的八个角是奇异点，并且每个面的边界会出现噪音，难以保证内部的精度，另外在立方体网格的 6个表面和12条棱上其矢量风场的转换比较困难。

(3)阴阳网格，其设计简单并且是正交的，同时也是准均匀、无奇异的结构化网格，但重叠区域的守恒条件难以满足、精度较低，并且重叠区构造的守恒格式计算量大，会引起计算负载的不平衡。

(4)二十面体网格，其没有奇点和没有分辨率聚集或膨胀问题，并且最大网格距和最小网格距之比接近于1，是一种比较理想的球面网格。从目前正在开发和业务化的全球大气格点模式来看，二十面体球面网格是最理想的选择。

但现在这些模式大部分是在六边形的C网格上进行离散的，后来公开了一种LE网格(一种新水平网格的Rossby波计算特性，地球物理学报，2006， 49(3):650～661)，将位势放在格点上，u，v风场同时放在x和y方向格点中间，可以有效避免科氏力项的空间平均，提高计算精度。在正方形网格上针对惯性重力波和Rossby波对C跳点网格和LE跳点网格进行过对比，认为 LE网格性能最好。从自由度数角度来分析，若动量和质量间的自由度数比为2：1，则不会产生虚假的计算模态。对于二十面体的六边形网格来说，其比率接近3：1，即动量节点太多，会产生虚假的计算模态。二十面体网格既可以用六边形网格划分，也可以用对等的三角形网格划分，那么如何评估三角形划分的LE网格上得到的计算频散性能与六边形C网格的，确定那个性能更佳便于日后模式设计者参考。

衡量离散方案是否合理的标准一般是分析离散方案是否满足微分情况下局地波动的传播特性和一些不变量(例如：质量和总能量等)的全球守恒性质。正压原始方程模式满足质量守恒、平流过程中全球大气动量守恒、全球大气总能量守恒、全球大气位涡和位涡拟能守恒、不计地球自转作用无辐散大气涡度拟能和水平尺度守恒(平均波数守恒)。王斌等(1990)设计了显式完全平方守恒差分格式，并在正压原始方程模式中进行了检验，效果非常理想，但遗憾的是：二十面体网格是非结构化网格，仅能采用有限体积法，难以应用显式完全平方守恒差分格式。

发明内容

本发明的目的是提供一种评估二十面体的LE网格计算特性方法及离散方法，解决了二十面体网格难以应用显式完全平方守恒差分格式的问题，通过设计离散方案，在正压模式中证明了离散形式下的质量守恒和总能量守恒性质。

为了达到上述目的，本发明提供了一种用于评估二十面体的LE网格计算特性的离散方法，该方法包含：

(1)针对用LE跳点网格离散，定义二十面体三角形LE跳点网格中的要素，LE跳点网格把质量点放在三角形的顶点上，速度点放在六边形边和三角形边的交点上，对于三角形网格速度点上的离散矢量场F表示为

其中F^⊥和F^||分别表示矢量场F的法向和切向分量；

(2)基于LE跳点网格，对方程组进行离散，并给出离散算子，针对三角形网格分别按照散度、旋度、导数和切向导数的定义，给出其任一矢量场F 的散度表达式、其任一矢量场F的旋度表达式、其任一标量场g的法向导数表达式、其任一标量场g的切向导数表达式；

其中，所述方程组为：

式(1)和(2)中，

为速度矢量，f为科氏参数，

为z方向的单位矢量，g为重力加速度，h为深度，t为时间；

所述任一矢量场F的散度表达式为：

所述任一矢量场F的旋度表达式为：

所述任一标量场g的法向导数表达式为：

所述任一标量场g的切向导数表达式为：

式中，(g)_ls表示在l_s上的标量场，g_S(ls,2)表示在l_s边上的第二个三角形单元中心的标量场，g_S(ls，1)表示在l_s边上的第一个三角形单元中心的标量场；

(3)速度的矢量形式可表示为

将步骤(2)中的动量方程(1)分解为分量形式，并将步骤(2)中的离散算子代入分量形式中，得到相应的离散形式；将步骤(2)给出的离散算子代入步骤(2)中的连续方程(2)，得到其空间离散形式；

所述动量方程(1)分解为分量形式为：

其中，l_S为三角形网格单元的边，

为与边l_S相切并指向

左侧的单位矢量，

为垂直于边l_S并指向外的单位矢量；

将步骤(2)中的离散算子代入动量方程(1)的分量形式中，得到相应的离散形式为：

其中，符号

表示利用S点上的平均值在

方向进行导数的计算，S点为三角形网格单元的中心；

所述空间离散形式为：

其中，

是三角形网格单元速度点上的高度场值，由S点上的值平均得到，为

De_L为相邻三角形网格单元中心之间的距离；E_s(L_s)为三角形网格单元的边l_S的集合；A_S为三角形网格单元的面积；h_s为三角形网格单元中心的深度；h_S1为三角形网格单元S相邻的第一个三角形网格单元的深度；h_S2为三角形网格单元S相邻的第二个三角形网格单元的深度；v^⊥为垂直于三角形网格单元边的速度。

优选地，

在S点上的平均值为：

式中，N表示法线方向单位矢量；v^⊥为垂直于三角形网格单元边的速度； T表示切向方向单位矢量；v^‖表示平行于三角形网格单元边的速度；下标1、 2、3分别表示三角形网格单元S周围的三个三角形网格单元。

优选地，

在L点上的平均值为：

本发明还提供了一种评估二十面体的LE网格计算特性的方法，该方法基于所述的离散方法，针对该离散方法进行质量守恒和能量守恒分析。

优选地，对式(3)进行全球所有三角形网格单元求和，相邻三角形网格单元中

的方向相反，相互抵消，式(3)变为：

该离散方法满足质量守恒。

优选地，把动能放在三角形网格单元的中心，根据动量方程，得到：

式中，g表示重力加速度；K_s表示动能；该离散方法满足能量守恒。

本发明的评估二十面体的LE网格计算特性方法及离散方法，解决了二十面体网格难以应用显式完全平方守恒差分格式的问题，具有以下优点：

本发明的方法评估了二十面体的LE网格的计算特性，设计了相应的概念和一些变量，并给出了对应的离散算子，在正压模式中证明了离散形式下的质量守恒和总能量守恒性质。本发明的离散格式比现在常用的二十面体六边形C网格上的格式更加保形，误差更小，质量和总能量守恒型更好。

本发明的方法可应用到现在的静力或非静力大气数值模式中，是一个值得推广的方法，还可将此方法应用到现有的某一全球模式中。

附图说明

图1为本发明的二十面体的三角形网格示意图。

图2为本发明初始时刻余弦钟的位势高度场示意图。

图3为本发明在方案1上模拟的第6天的位势高度场示意图。

图4为本发明在方案1上模拟的第12天的位势高度场示意图。

图5为本发明在方案2上模拟的第6天的位势高度场示意图。

图6为本发明在方案2上模拟的第12天的位势高度场示意图。

图7为本发明在方案1上模拟的第12天的解析解位势高度差示意图。

图8为本发明在方案2上模拟的第12天的解析解位势高度差示意图。

图9为本发明方案1和方案2计算的质量场示意图。

图10为本发明方案1和方案2计算的总能量示意图。

图11为本发明方案1的模拟误差示意图。

图12为本发明方案2的模拟误差示意图。

具体实施方式

下面将对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

一种评估二十面体的LE网格计算特性的方法，包含：

(1)针对用LE跳点网格离散，定义二十面体三角形LE跳点网格中的要素，三角形的LE跳点网格把质量点放在三角形的顶点上，也就是六边形单元的中心，速度点放在六边形边和三角形边的交点上，也就是六边形边的中点，对于三角形网格速度点上的离散矢量场F，可以表示为

其中F^⊥和F^||分别表示矢量场F的法向和切向分量(上标⊥表示法向，‖表示切向)；

(2)基于LE网格进行离散，并给出离散算子，针对三角形网格分别按照散度、旋度、导数和切向导数的定义，给出其任一矢量场的散度表达式、其任一矢量场F的旋度表达式、其任一标量场g的法向导数表达式、其切向导数；

(3)速度的矢量形式可表示为

将动量方程分解为分量形式，并将步骤(2)中的离散算子代入分量形式中，得到相应的离散形式；同样将步骤(2)给出的离散算子代入连续方程，得到其空间离散形式；

(4)在正压模式中验证离散形式下的质量守恒和总能量守恒性质：若对空间离散形式进行全球所有三角形网格单元求和，由于

是三角形网格单元 S的边l_S上的单位矢量，如果S单元为正，则在和边l_S相邻的另一个单元中为负，这样对全球所有三角形网格单元求和时，

的符号一正一负同时出现，相互抵消，则

离散方案满足质量守恒；把动能也放在三角形网格单元的中心，得到

总能量也是守恒的。

为了对本发明提供的评估二十面体的LE网格计算特性的方法进行详细说明，以下通过实施例1进行具体阐述。

实施例1

1、定义二十面体三角形LE网格中的要素

针对二十面体的LE网格，该正二十面体的每一个面可以进一步细分，网格细化的级数用k表示，单元的面数用NF(k)表示，NF(k)＝2^2kNF(0)，其中 NF(0)表示最初的二十面体网格的面数，即NF(0)＝20。网格的边数用NE(k)表示，

网格的顶点数用NV(k)表示，

例如，当网格细化到第6级时，单元面数为81920，网格边数为122880，网格顶点数为40962，平均分辨率为112公里。

如图1所示，为本发明的二十面体的三角形网格示意图，当把每个正二十面体六边形的中心相连就构成球面三角形网格，并且六边形网格的边与三角形网格的边正交，六边形网格和三角形网格的作用是完全对等的，即球面六边形网格和球面三角形网格是完全等效的，二者均可以作为球面划分的方法，本发明基于以二十面体三角形网格作为球面划分的网格，采用LE网格离散，需定义的变量如表1所示：

表1三角形网格要素定义

三角形网格单元的网格点S就是指三角形网格单元的网格中心。假定

向外(指向三角形外)为正，向内(指向三角形内)为负。同样，给定

表示指向

左侧为正，指向右侧为负。这样二者之间的关系为

k_l也是向外为正，向内为负。

指向外为正，指向内为负。

指向

左侧地为正，指向

右侧的为负。因此，

三角形网格和六边形网格是完全对等的，则二者之间的单位矢量之间存在下列关系：

关于三角形的LE跳点网格，是把质量点放在三角形的顶点上，也就是六边形单元的中心，速度点放在六边形边和三角形边的交点上(即六边形边的中点，也是三角形边的中点)。通过六边形网格单元的边l_L相连的两个三角形网格单元分别用S(l_L,1)和S(l_L,2)表示，从S(l_L,1)到S(l_L,2)的方向为三角形网格单元S(l_L,1)法向矢量

的正方向。对于给定的边l_L，L(l_L,1)和L(l_L,2)分别表示通过边l_L相邻两个六边形网格单元，也是穿过该边l_L三角形网格单元的两个顶点。从L(l_L,1)到L(l_L,2)的方向为三角形网格单元S(l_L,1)矢量

的正方向。对于一般的球面三角形网格速度点上的离散矢量场F，可以表示为

其中F^⊥和F^||分别表示矢量场F的法向和切向分量(上标⊥表示法向，‖表示切向)。

2、离散算子

基于LE网格，对浅水方程组进行离散，方程组为：

其中，

为速度矢量，f为科氏参数，

为z方向的单位矢量(采用笛卡尔坐标系)，g为重力加速度，h为深度，t为时间。

对上述方程组进行离散，需要给出相应的离散散度、旋度、切向导数和法向导数。对于三角形网格而言，按照散度的定义，其任一矢量场F的散度表达式为：

根据旋度的定义，其任一矢量场F的旋度表达式为：

根据导数的定义，其任一标量场g的法向导数表达式为：

其切向导数为：

3、正压方程的空间离散

在上述空间离散算子的基础上，速度的矢量形式可表示为

其中，上标⊥表示垂直，上标‖表示平行。

动量方程(1)分解为分量形式可表示为：

将步骤2给出的离散算子带入上式，得到相应的离散形式为

其中符号

表示利用S点上的平均值在

方向进行导数的计算，则

在S点上的平均值通过下式得到：

而

在L点上的平均值通过下式得到：

同样将步骤2给出的离散算子，代入连续方程(2)，得到其空间离散形式为：

其中，

是三角形网格单元速度点上的高度场值，由S点上的值平均得到，即

4、质量守恒

若对(3)式进行全球所有三角形网格单元求和，由于

是三角形网格单元S的边l_S上的单位矢量，如果S单元为正，则在和边l_S相邻的另一个单元中为负，这样对全球所有三角形网格单元求和时，

的符号一正一负同时出现，相互抵消。因此(3)式最后变为：

即离散方案是满足质量守恒的，对应于微分形式的

5、总能量守恒

当把动能也放在三角形网格单元的中心，也同样可以得到：

式中，g表示重力加速度；K_s表示动能。总能量也是守恒的，对应于微分形式的

对本发明的离散方法进行数值测试，具体如下：

1、数值测试以验证离散方案的精度和合理性

以下的数值测试的时间积分方案采用Runge－Kutta的三阶方法。

取预报变量

模式方程取为Φ_t＝R(Φ)，Runge-Kutta积分采用三步形式进行：

Φ^t+Δt＝Φ^t+ΔtR(Φ^**)

这里的Δt为时间步长，上标表示时间层。下面通过余弦钟测试对上述离散方案进行检验，以验证离散方案的精度和合理性，验证时采用以下误差范数：

其中，

这里的F_n(S)表示在三角形网格单元S上计算得到的结果，F₀(S)表示三角形网格单元S上的解析解。采用网格顶点数为40962的分辨率进行测试，其平均分辨率为112km。

该测试例在每一时间步长用解析指定的平流风场替换预报的风场。因为风场是无辐散的，因此自由面高度方程退化为平流方程。

余弦钟定义为：

v＝-u₀ sinλsinβ

其中，h₀＝1000，r是

和中心

之间的大圆距离，

钟的半径R＝a/3，a为地球半径，赤道速度u₀＝2πa/(12days)，大约为40m/s，参数β是固体旋转轴和球坐标系极轴之间的夹角。这里取β＝0。从定义式可以看出：该测试例模拟12天后，环绕地球一周后回到原地，形状保持不变，如图2所示，为初始时刻余弦钟的位势高度场。

2、检验保形特性

为了检验在三角形LE网格上(简称方案1)和六边形C网格上(简称方案2)计算结果的保形情况，给出了这两种方案在第6天和第12天模拟的位势高度场(图3-6)，图中等值线间隔为100米。比较图3、图5和图2可以看出：模拟计算的第6天，方案2上计算的结果等值线(图5)和初始时刻(图2)有所差异，在东侧等值线分布变稀，呈现不对称分布(初始时刻是对称分布的)。而方案1上模拟得到的等值线(图3)和初始时刻(图2) 一样仍呈现对称分布。模拟计算的第12天结果差别更加明显，方案1(图4) 在第12天都回到了初始时刻(图2)的位置，但方案2的结果(图6)中心没有回到初始时刻的位置，而是略偏西，并且等值线分布也不对称，西部更密，东部更疏。

为了更加直接反映二者在第12小时的对称情况，将这两种方案与初始时刻的结果相减，绘制了图7和图8，图中等值线间隔为2米，虚线为负值，实线为正值。从图7和8可以看到，方案1上在东西方向是对称分布的，误差的最小值为-8m，最大值为11m。而方案2上明显呈现不对称分布，中心西侧的误差为正值，中心东侧的误差为负值，外围也有些正值和负值等值线，但也是不对称分布，并且误差最小值为-199m，最大值为162m，由此可见，在方案1得到的结果更加保形。

3、检验守恒性质

按照公式(4)和(5)绘制了这两种方案相对于初始时刻的质量场和总能量相对误差随天数的变化情况，分别见图9和图10。从图中可以看出：和方案2相比，方案1的质量和总能量相对误差都较小。但从量级来看：质量的相对误差量级为10^-15，总能量的相对误差在10^-6左右，可以说，尽管方案 1的质量和总能量误差更小，但方案1和方案2的质量和总能量误差均很小，二者均满足质量和总能量守恒，只是方案1更好些。

图9和10中，两种方案上计算的质量场(单位：m·m²)和总能量相对误差(单位：m·m²/s²·m²)随天数的变化，图中纵坐标变量括号中的数值是该变量的量级大小。

4、误差分析

为了定量给出这两种方案的模拟误差，按照公式(4)和(5)计算他们的范数误差(图11和12)，图中在方案1(图11)上和在方案2(图12) 上计算的误差范数L₂、L_∞随天数的变化。

从图中可以看出，方案1的误差量级比方案2小10^-1左右，其中方案1 的L_∞最大值为0.0671，而方案2的L_∞最大值为0.256。可见，方案1确实比较优越。

尽管本发明的内容已经通过上述优选实施例作了详细介绍，但应当认识到上述的描述不应被认为是对本发明的限制。在本领域技术人员阅读了上述内容后，对于本发明的多种修改和替代都将是显而易见的。因此，本发明的保护范围应由所附的权利要求来限定。

Claims (6)

1.一种用于评估二十面体的LE网格计算特性的离散方法，其特征在于，该方法包含：

(1)针对用LE跳点网格离散，定义二十面体三角形LE跳点网格中的要素，LE跳点网格把质量点放在三角形的顶点上，速度点放在六边形边和三角形边的交点上，对于三角形网格速度点上的离散矢量场F表示为

其中F^⊥和F^||分别表示矢量场F的法向和切向分量；

(2)基于LE跳点网格，对方程组进行离散，并给出离散算子，针对三角形网格分别按照散度、旋度、导数和切向导数的定义，给出其任一矢量场F的散度表达式、其任一矢量场F的旋度表达式、其任一标量场g的法向导数表达式、其任一标量场g的切向导数表达式；

其中，所述方程组为：

式(1)和(2)中，

为速度矢量，f为科氏参数，

为z方向的单位矢量，g为重力加速度，h为深度，t为时间；

所述任一矢量场F的散度表达式为：

所述任一矢量场F的旋度表达式为：

所述任一标量场g的法向导数表达式为：

式中，(g)_ls表示在l_s上的标量场，g_S(ls,2)表示在l_s边上的第二个三角形单元中心的标量场，g_S(ls，1)表示在l_s边上的第一个三角形单元中心的标量场；

所述任一标量场g的切向导数表达式为：

(3)速度的矢量形式可表示为

将步骤(2)中的动量方程(1)分解为分量形式，并将步骤(2)中的离散算子代入分量形式中，得到相应的离散形式；将步骤(2)给出的离散算子代入步骤(2)中的连续方程(2)，得到其空间离散形式；

所述动量方程(1)分解为分量形式为：

其中，l_S为三角形网格单元的边，

为与边l_S相切并指向

左侧的单位矢量，

为垂直于边l_S并指向外的单位矢量；

将步骤(2)中的离散算子代入动量方程(1)的分量形式中，得到相应的离散形式为：

其中，符号

表示利用S点上的平均值在

方向进行导数的计算，S点为三角形网格单元的中心；

所述空间离散形式为：

其中，

是三角形网格单元速度点上的高度场值，由S点上的值平均得到，为

De_L为相邻三角形网格单元中心之间的距离；E_s(L_s)为三角形网格单元的边l_S的集合；A_S为三角形网格单元的面积；h_s为三角形网格单元中心的深度；h_S1为三角形网格单元S相邻的第一个三角形网格单元的深度；h_S2为三角形网格单元S相邻的第二个三角形网格单元的深度；v^⊥为垂直于三角形网格单元边的速度。

2.根据权利要求1所述的用于评估二十面体的LE网格计算特性的离散方法，其特征在于，

在S点上的平均值为：

式中，N表示法线方向单位矢量；v^⊥为垂直于三角形网格单元边的速度；T表示切向方向单位矢量；v^‖表示平行于三角形网格单元边的速度；下标1、2、3分别表示三角形网格单元S周围的三个三角形网格单元。

3.根据权利要求1所述的用于评估二十面体的LE网格计算特性的离散方法，其特征在于，

在L点上的平均值为：

4.一种评估二十面体的LE网格计算特性的方法，其特征在于，该方法基于如权利要求1-3中任意一项所述的离散方法，针对该离散方法进行质量守恒和能量守恒分析。

5.根据权利要求4所述的评估二十面体的LE网格计算特性的方法，其特征在于，对式(3)进行全球所有三角形网格单元求和，相邻三角形网格单元中

的方向相反，相互抵消，式(3)变为：

该离散方法满足质量守恒。

6.根据权利要求4所述的评估二十面体的LE网格计算特性的方法，其特征在于，把动能放在三角形网格单元的中心，根据动量方程，得到：

式中，g表示重力加速度；K_s表示动能；

该离散方法满足能量守恒。

Images