[发明专利]基于分裂算法的湖泊三维水动力-水温-水质模拟预测方法

权利要求书

1.一种基于分裂算法的湖泊三维水动力-水温-水质模拟预测方法，其特征在于，包括以下步骤：

(1)构建湖泊三维水动力-水温-水质模型，模型由湖泊水动力控制方程组、水温控制方程和水质控制方程组成，将湖泊计算域离散成若干网格单元，并对网格采用ArakawaC模式布置变量，采用分裂算法首先将湖泊三维水动力-水温-水质模型中的水动力控制方程组的各算子按照其物理波动过程的波频快慢特性进行分类，对低频慢过程算子采用显式处理，对高频快过程算子采用隐式处理，通过分裂算法分步离散并求解水动力控制方程组，得到湖泊水域内不同位置和时间的三维流场u、v、w及其水深H，其中u、v、w分别表示x、y、z方向的流速；

(2)通过步骤(1)求得湖泊三维流场u、v、w和水深H等水动力参数的基础上，采用分裂算法分别将水温控制方程和水质控制方程中的算子分为低频慢过程算子和高频快过程算子两类，对其中的低频慢过程算子采用显式处理，对高频快过程算子采用隐式处理，通过分裂算法分步离散并求解水温控制方程和水质控制方程，得到湖泊水域内不同位置和时间的水温T和水质指标浓度C。

2.根据权利要求1所述的预测方法，其特征在于，步骤(1)包括以下子步骤：

(1.1)构建湖泊三维水动力-水温-水质模型，模型由反映湖泊水动力、水温和水质复杂变化过程的一组数学物理控制方程组构成，包括：

水动力控制方程组：

连续方程∂xHu+∂yHv+∂σω+∂tη=0---(1)]]>

水位方程∂tη+∂x(H∫01udσ)+∂y(H∫01vdσ)=0---(2)]]>

x动量方程∂tu+u∂xu+v∂yu+ω∂σu=fv-g∂xη+∂x(νh∂xu)+∂y(νh∂yu)+∂σ(νvH2∂σu)---(3)]]>

y动量方程∂tv+u∂xv+v∂yv+ω∂σv=-fu-g∂yη+∂x(νh∂xv)+∂y(νh∂yv)+∂σ(νvH2∂σv)---(4)]]>

水温控制方程：∂tT+u∂xT+v∂yT+ω∂σT=∂x(AHT∂xT)+∂y(AHT∂yT)+∂σ(AVTH2∂σT)---(5)]]>

水质控制方程：∂tC+u∂xC+v∂yC+ω∂σC=∂x(AHC∂xC)+∂y(AVC∂yC)+∂σ(vbH2∂σC)---(6)]]>

式中，是σ坐标变换；ω=wH-σH∂H∂t+uH(∂h∂x-σ∂H∂x)+vH(∂h∂y-σ∂H∂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)将湖泊离散成若干网格单元，并对网格采用ArakawaC模式布置变量：对一个正方体网格单元，规定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方向为例：

u1i+1/2,j,k=-Δt2Δxui+1/2,j,kn(ui+3/2,j,kn-ui-1/2,j,kn)-vi+1/2,j,kn*Δt2Δy(ui+1/2,j+1,kn-ui+1/2,j-1,kn)-ωi+1/2,j,kn*ΔtΔσk+0.5(Δσk+1+Δσk-1)(ui+1/2,j,k+1n-ui+1/2,j,k-1n)+vhΔtΔx2(ui+3/2,j,kn-ui-1/2,j,kn)+vhΔtΔy2(ui,j+1,kn-2ui,j,kn+ui,j-1,kn)+(1-2vhΔtΔx2)ui+1/2,j,kn-fvi+1/2,j,kn*---(7)]]>

式中：

vi+1/2,j,kn*=14(vi,j-1/2,kn+vi,j+1/2,kn+vi+1,j-1/2,kn+vi+1,j+1/2,kn)---(8)]]>

ωi+1/2,j,kn*=14(ωi,j,k+1/2n+ωi,j,k-1/2n+ωi+1,j,k+1/2n+ωi+1,j,k-1/2n)---(9)]]>

(1.5)由显式离散水位方程(2)，求得下一个时间节点的静水位偏移位移ηⁿ⁺¹，根据已经求到的中间流速u1、v1和静水位偏移位移ηⁿ⁺¹，对高频快过程算子重力梯度项做隐式处理，求得x方向中间流速u2和y方向中间流速v2，以x方向为例：

u2i+1/2,j,k-u1i+1/2,j,kΔt=-gηi+1,jn+1-ηi,jn+1Δx---(10)]]>

(1.6)根据已经求到的中间流速u2和v2，对高频快过程算子垂向涡粘项做隐式处理，整理后分别得到关于下一个时间节点的x方向流速uⁿ⁺¹和y方向流速vⁿ⁺¹的线性方程式，以x方向为例：

UT·uk+1n+1+UB·uk-1n+1+UC·ukn+1=UF---(11)]]>

式中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。

3.根据权利要求书2所述的模拟预测方法，其特征在于，步骤(2)包括以下子步骤：

(2.1)根据步骤(1.3)类似的原理，将水温控制方程(5)和水质控制方程(6)的算子分为低频慢过程算子和高频快过程算子两类，其中，低频慢过程算子包括对流项、水平扩散项，高频快过程算子包括垂向扩散项；

(2.2)对低频慢过程算子对流项和水平扩散项做显式处理，求得中间温度T1和中间浓度C1；

(2.3)根据求得的T1、C1，对高频快过程算子垂向扩散项做隐式处理，求得下一个时间节点的温度Tⁿ⁺¹、浓度Cⁿ⁺¹，整理得到一个线性方程式，以温度T为例：

MTi,j,kn+1+STi,j,k+1n+1+PTi,j,k-1n+1=R---(12)]]>

M、S、P、R是包含已知变量T1的参数，在各个网格写出式(12)，分别得到一个涉及水温时空关系和水质时空关系的线性矩阵方程组，采用托马斯法求得Tⁿ⁺¹，并用同样的方法求得Cⁿ⁺¹；

(2.4)重复步骤(2.2)-(2.3)，循环迭代计算，得到湖泊水域内不同位置和时间的水温T和水质指标浓度C。