非静压模型垂向网格分离计算方法
【技术领域】
[0001] 本发明涉及一种非静压模型垂向网格分离计算方法，属于技术领域。
【背景技术】
[0002] 目前，静压假设的数学模型已在河口海岸水动力研究中得到广泛应用。但在地形 剧烈变化、密度梯度较大、水流急剧变化等具有垂向流速变化较大的区域，静压模型难以准 确模拟水动力的变化。故此，如何发展非静压模型，已成为近年河口海岸数值模拟研究的热 点方向。大多数非静压模型压力值的计算采用分步法，即单独计算压力的动压部分与静压 部分，而求解动压时需解泊松方程，此部分在非静压模型中耗时严重，影响了非静压模型的 计算效率及其应用范围。
【发明内容】
[0003] 为解决现有技术的不足，本发明的目的在于提供一种非静压模型垂向网格分离计 算方法，通过分离求解泊松方程和动量方程的计算网格，解决非静压模型计算过程中耗时 过多的问题。
[0004] 为了实现上述目标，本发明采用如下的技术方案：
[0005] -种非静压模型垂向网格分离计算方法，其特征是，在垂直方向采用不同网格求 解压力泊松方程和动量方程，在求解泊松方程时采用粗网格，求解动量方程时采用细网格， 具体包括如下步骤：
[0006] 1)采用细网格求解不含动压项的动量方程得到速度中间值；
[0007] 2)采用粗网格求解压力泊松方程，通过分段平均得到粗网格速度中间值，得到粗 网格上的动压值，并将此动压值插值到细网格中得到细网格的动压值；
[0008] 3)通过动量方程的动压项修正速度场，得到最终的速度值。
[0009] 前述的非静压模型垂向网格分离计算方法，其特征是，所述步骤2)具体包括如下 步骤：
[0010]2. 1)通过分段平均法求得粗网格速度中间值if，分段平均法表达式如下：
[0011] (1)
[0012] 式中i表示粗网格上的网格编号，队表示在粗网格第i个网格内细网格的个数，Ji 表示第i个粗网格中细网格编号的最小值；
[0013]2. 2)通过泊松方程（3)求解粗网格上的动压值p'，
[0014]
(2)
[0015] 式中D表示水深，P表示密度，x，y，0分别表示水平方向坐标和垂向0坐 标，X#，/表示笛卡尔坐标系水平方向坐标，u$，/，V表示速度中间值；
[0016] 2. 3)由粗网格的动压值通过插值得到细网格上的动压值p。
[0017] 前述的非静压模型垂向网格分离计算方法，其特征是，所述步骤2)中插值的求取 采用三次样条插值法，插值函数需满足如下条件：
[0018] r3)
[0019] (4)
[0020] (5)
[0021] 式中Z= -h表示位于水底处，h表示静水面以下水深，n表示水体表面水位；p(1) 表示粗网格第i个网格中插值函数，p(1)'表示插值函数p(1)的一阶导数，P(1)"表示插值函 数P(1)的二阶导数，0i表示粗网格第i网格底面。坐标，Pl表示粗网格第i网格动压值。
[0022] 本发明所达到的有益效果：通过减少计算泊松方程的计算网格，在计算精度没有 降低的条件下，显著减少非静压模型的计算时间，测试算例显示，在计算泊松方程的垂向网 格数最少可以仅为动量方程网格数的十分之一，采用网格分离方法后，计算时间最多可减 少约80%。
【附图说明】
[0023] 图1是非静压模型网格分离方法的网格布置图；
[0024]图2是非静压模型网格分离方法计算步骤图；
[0025] 图3是Lock-Exchange算例初始状态图；
[0026] 图4是不同网格设置密度计算结果。
【具体实施方式】
[0027] 下面结合附图对本发明作进一步描述。以下实施例仅用于更加清楚地说明本发明 的技术方案，而不能以此来限制本发明的保护范围。
[0028]本发明涉及的一种非静压模型垂向网格分离计算方法，在垂直方向采用不同网格 求解压力泊松方程和动量方程，减少非静压模型求解压力泊松方程的时间，从而提高计算 效率。
[0029] 求解泊松方程的网格数较少，称为粗网格；求解动量方程的网格数较多，称为细网 格，网格布置图如图1所示，实线表示粗网格，虚线表示细网格，P'，u'，w'分别表示粗网格 上动压值和水平、垂向速度值，P，u,w表示细网格上动压值和水平、垂向速度值。
[0030] 非静压模型目前主要的计算方法为分步法，即将压力值分为动压部分和静压部 分，第一步通过求解不含动压项的动量方程得到速度中间值，然后求解压力泊松方程，得到 动压值，最后通过动量方程的动压项修正速度场，得到最终的速度值。
[0031] 而本发明采用的非静压模型垂向网格分离计算方法是在非静压模型计算分步法 的基础上进行改进，分步法的第一步不变，第二步由于求解泊松方程采用粗网格，需要求解 粗网格的速度中间值，因此增加通过分段平均得到粗网格速度中间值的步骤，在粗网格上 求解压力泊松方程得到动压值后需要增加动压值插值到细网格的步骤，而后的第三步修正 速度场不变。
[0032] 第一步和第三步由于都不变，就不进行阐述了，进行改进的第二步中包括如下步 骤：
[0033] 1)在细网格中求解不含动压项的动量方程（1)，得到速度中间值U%
[0034]
[0035] 式中，分别为动量方程中底摩阻项、斜压项和紊流扩散项，上标n表示模 型计算的第n时间步，F,G,H表示通量项；
[0036]2)通过分段平均法求得粗网格速度中间值if，分段平均法表达式如下：
[0037] (2)
[0038] 式中i表示粗网格上的网格编号，队表示在粗网格第i个网格内细网格的个数，Ji 表示第i个粗网格中细网格编号的最小值；
[0039]3)通过泊松方程（3)求解粗网格上的动压值p'，
[0040] (3)；
[0041] 4)由粗网格的动压值通过插值得到细网格上的动压值p，插值方法采用三次样条 插值，插值函数需满足如下条件：
[0042] C4)
[0043] (5)
[0044] 边界条件：
（6)
[0045] 5)通过上一步得到的动压值p，利用方程（7)修正细网格上的速度场，并利用方程 (8)更新自由表面水位：
[0046]
C7) (8)
[0048] 公式（7)中U(1)龙格库塔方法第一阶段的速度值，Sf:为动压源项，At为时间步 长；公式（8)中D表示水深，u，v为水平流速；
[0049]6)此步骤根据具体算例确定，考虑紊流时，可求解紊流模型，包含k-G，大涡模拟 等；在包含盐度、泥沙的模拟算例中，可通过盐度、泥沙对流扩散方程（9)求解盐度、泥沙等 的分布和输移。
[0050]
C9)
[0051] 其中可表示为温度、盐度、泥沙浓度等，i，j，k为水