基于分裂算法的湖泊三维水动力‑水温‑水质模拟预测方法

摘要

本发明公开了一种基于分裂算法的湖泊三维水动力‑水温‑水质模拟预测方法。构建湖泊三维水动力‑水温‑水质模型，将湖泊离散成若干网格单元，并对网格采用Arakawa C模式布置变量，基于分裂算法将湖泊三维水动力‑水温‑水质模型中各算子按照其物理波动过程的波频快慢特性进行分类，对低频慢过程算子采用显式处理，对高频快过程算子采用隐式处理，运用分裂算法分步离散求解模型，得到湖泊水域内不同位置和时间的三维流场、水温和水质指标浓度。本发明提出的基于分裂算法的湖泊三维水动力‑水温‑水质模拟预测方法能较准确地反映湖泊水体中的动量迁移、热量传递和污染物输移等复杂的物理过程，具有计算稳定性好、计算精度和计算效率高的特点。

主权项

一种基于分裂算法的湖泊三维水动力‑水温‑水质模拟预测方法，其特征在于，包括以下步骤：(1)构建湖泊三维水动力‑水温‑水质模型，将湖泊计算域离散成若干网格单元，并对网格采用Arakawa C模式布置变量，采用分裂算法将湖泊三维水动力‑水温‑水质模型水动力控制方程组中水平动量方程各算子按照其物理波动过程的波频快慢特性进行分类，对低频慢过程算子采用显式处理，对高频快过程算子采用隐式处理，通过分裂算法分步离散并求解水动力控制方程组，得到湖泊水域内不同位置和时间的三维流场u、v、w及其水深H，其中u、v、w分别是x、y、z方向的流速；具体包括以下子步骤：(1.1)构建湖泊三维水动力‑水温‑水质模型，模型由反映湖泊水动力、水温和水质复杂变化过程的一组数量物理控制方程组构成，包括：水动力控制方程组：连续方程水位方程x动量方程y动量方程水温控制方程：水质控制方程：式中，是σ坐标变换；是σ变换后的垂向流速；x，y为水平坐标，z为垂向坐标；h为自由水面到水底的水深值，η为静水位偏移位移，H＝h+η为总水深；u，v分别为x，y方向的流速；t为时间；T为水温；C为污染物浓度；Ω为地球自转角速度，为所处纬度；V_h和V_v分别为水平涡黏系数和垂向涡粘系数；V_b为垂向扩散系数；g为重力加速度；A_HT和A_VT为水平温度扩散系数和垂向温度扩散系数；A_HC和A_VC为水平水质扩散系数和垂向水质扩散系数；(1.2)将湖泊离散成若干网格单元，并对网格采用Arakawa C模式布置变量：对一个正方体网格单元，规定i、j和k分别为x、y和z方向的网格编号索引，变量ω设置在上下表面的中点上(i，j，k±1/2)，在每一个横截面中，变量u设置在网格左右两边的中心(i±1/2，j，k)，变量v设置在前后两边的中点上(i，j±1/2，k)，水深H、水质浓度C、水温T设置在网格中心(i，j，k)；(1.3)采用分裂算法将湖泊三维水动力‑水温‑水质模型水动力控制方程组中水平动量方程各算子按照其物理波动过程的波频快慢特性进行分类：对应于缓行的内重力波产生的作用力和表面重力长波产生的作用力的差异，将水平动量方程(3)、(4)的算子分为低频慢过程算子和高频快过程算子两类，其中，低频慢过程算子包括对流项、科氏力项和水平涡粘项，高频快过程算子包括重力梯度项和垂向涡粘项；(1.4)对对流项、科氏力项和水平涡粘项这些低频慢过程算子做显式处理，求得x方向中间流速u1和y方向中间流速v1，以x方向为例：$\begin{array}{c}u{1}_{i+1/2,j,k}=-\frac{\Delta t}{2\Delta x}{u}_{i+1/2,j,k}^n\left(u_{i+3/2,j,k}^n-u_{i-1/2,j,k}^n\right)-\frac{v_{i+1/2,j,k}^{n*}\Delta t}{2\Delta y}\left(u_{i+1/2,j+1,k}^n-u_{i+1/2,j-1,k}^n\right)\\ -\frac{{\omega }_{i+1/2,j,k}^{n*}\Delta t}{{\mathrm{\Delta \sigma }}_k+0.5\left({\mathrm{\Delta \sigma }}_{k+1}+{\mathrm{\Delta \sigma }}_{k-1}\right)}\left(u_{i+1/2,j,k+1}^n-u_{i+1/2,j,k-1}^n\right)+\frac{v_h\Delta t}{{\mathrm{\Delta x}}^2}\left(u_{i+3/2,j,k}^n+u_{i-1/2,j,k}^n\right)\\ +\frac{v_h\Delta t}{{\mathrm{\Delta y}}^2}\left(u_{i,j+1,k}^n-2u_{i,j,k}^n+u_{i,j-1,k}^n\right)+\left(1-\frac{2v_h\Delta t}{{\mathrm{\Delta x}}^2}\right)u_{i+1/2,j,k}^n+{\mathrm{fv}}_{i+1/2,j,k}^{n*}\end{array}---\left(7\right)$式中：$v_{i+1/2,j,k}^{n*}=\frac{1}{4}(v_{i,j-1/2,k}^n+v_{i,j+1/2,k}^n+v_{i+1,j-1/2,k}^n+v_{i+1,j+1/2,k}^n)---\left(8\right)$${\omega }_{i+1/2,j,k}^{n*}=\frac{1}{4}({\omega }_{i,j,k+1/2}^n+{\omega }_{i,j,k-1/2}^n+{\omega }_{i+1,j,k+1/2}^n+{\omega }_{i+1,j,k-1/2}^n)---\left(9\right)$(1.5)显式离散水位方程(2)，求得下一个时间节点的静水位偏移位移ηⁿ⁺¹，根据已经求到的中间流速u1,v1和静水位偏移位移ηⁿ⁺¹，对高频快过程算子重力梯度项做隐式处理，求得x方向中间流速u2和y方向中间流速v2，以x方向为例：$\frac{u{2}_{i+1/2,j,k}-u{1}_{i+1/2,j,k}}{\Delta t}=-g\frac{{\eta }_{i+1,j}^{n+1}-{\eta }_{i,j}^{n+1}}{\Delta x}---\left(10\right)$(1.6)根据已经求到的中间流速u2和v2，对高频快过程算子垂向涡粘项做隐式处理，整理后分别得到关于下一个时间节点的x方向流速uⁿ⁺¹和y方向流速vⁿ⁺¹的线性方程式，以x方向为例：$UT·u_{k+1}^{n+1}+UB·u_{k-1}^{n+1}+UC·u_k^{n+1}=UF---\left(11\right)$式中UT、UB、UC，UF是包含已知变量u2的参数，在所有网格上写出式(11)，得到一个涉及三维流场时空关系的线性矩阵方程组，采用托马斯法求解矩阵方程组，求得uⁿ⁺¹，并用同样的方法求出vⁿ⁺¹；(1.7)根据已经得到的x方向流速uⁿ⁺¹、y方向流速vⁿ⁺¹和连续方程(1)，求得垂向流速ωⁿ⁺¹，根据σ反坐标变换求得笛卡尔坐标系下垂向流速wⁿ⁺¹；(1.8)重复步骤(1.4)‑(1.7)，循环迭代计算，得到湖泊水域内不同位置和时间的三维流场u、v、w及其静水位偏移位移η，并根据静水位偏移位移η进一步求得水深H；(2)通过步骤(1)求得包括湖泊三维流场u、v、w和水深H在内的水动力参数，采用分裂算法分别将水温控制方程和水质控制方程中的算子分为低频慢过程算子和高频快过程算子两类，对其中的低频慢过程算子采用显式处理，对高频快过程算子采用隐式处理，通过分裂算法分步离散并求解水温控制方程和水质控制方程，得到湖泊水域内不同位置和时间的水温T和水质指标浓度C。