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

Abstract

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

Description

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

技术领域

本发明属于湖泊水环境数值模拟技术领域，更具体地，涉及一种基于分裂算法的湖泊三维水动力-水温-水质模拟预测方法。本发明提出的湖泊三维水动力-水温-水质模拟预测方法，可用于河湖连通工程、引水调度、湖泊水体富营养化、湖泊水环境修复和管理等方面，为湖泊生态环境的治理和保护提供有效的技术手段和科学的决策分析。

背景技术

湖泊作为陆地水资源的重要载体，是地球水圈的关键组成单元，是人类赖以生存和可持续发展的自然依托，但目前我国湖泊却遭受着湖泊水资源萎缩、水污染和富营养化等水环境问题，已给我国社会经济发展和人民的生活环境带来了严重的影响。湖泊水环境保护和综合治理已经成为资源、环境、生态领域亟需解决的关键科学问题。本发明所提出的湖泊三维水动力-水温-水质模拟预测方法，对于促进我国水生态文明建设，打造天蓝、地绿、水美的“美丽中国”和实现“中国梦”具有重要的科学意义。

原型观测和数值模拟是湖泊水环境研究的主要技术手段，但由于湖流运动在时间和空间上呈现出高度复杂性，原型观测往往受条件限制，无法做到同步监测，而数值模拟，由于费用低、物理意义明确，并能较真实地反映湖泊水体的水动力、水温和水质状况和理化生过程，已经成为支持水环境研究的重要技术手段。合适的湖泊水环境模型和稳定高效的模型求解算法是湖泊水环境数值模拟技术的关键。目前，湖泊水环境数值模拟技术存在以下不足：

1.传统湖泊水环境模型求解方法通常采用单一的数值格式对模型进行离散求解，导致计算要求的稳定条件无法满足或计算工作量大、计算效率低等问题。例如，采用纯显格式对模型进行离散，可能因不满足CFL稳定判据而造成计算结果发散的现象；而采用纯隐格式对模型进行离散，虽然保证了计算稳定性，但是必须在整个求解区域上解大规模非线性代数方程组，导致计算效率低；

2.在实际问题中，模型的数学控制方程往往被简化成二维甚至是一维，虽然合理的简化易于求解，但是自然水体是三维的，简化的一、二维模型并不能真实地反映水体在纵向、横向、垂向上的变化过程；

3.湖泊水环境中物质输移、能量转换等物理过程联系紧密、相辅相成，湖泊水动力过程是水温热量传递过程和污染物质输移过程的驱动条件。而某些研究基于过于简单的假设，将水温、水质过程视为孤立的单一过程，并未考虑它们与水动力等过程的耦合效应，导致模拟的结果与真实的湖体情况存在较大差异。

随着科研人员对湖泊水体复杂过程更深刻的认识，国内外近年来发表的与本发明相关的主要研究成果有：

《水力发电学报》2005年第1期《非静压假定的σ坐标下垂向二维浅水模型的应用研究》对非静压假定的σ坐标下垂向二维浅水模型进行了应用研究。该文采用的是垂向二维的浅水模型，而本发明构建的是三维水动力-水温-水质模型。

《水动力学研究与进展A辑》2009年第24卷3期《基于非结构化网格上的三维微幅表面流动非静压数值模型》建立了基于非结构化网格求解三维微幅自由表面流动非静压数值模型，一种有限差分法和有限体积法相结合的方法被用来在非结构化网格上离散控制方程。该文建立的非静压三维模型是解决微幅自由表面的流动问题，而本发明解决的是湖泊水动力-水温-水质模拟预测问题。

《环境科学学报》2012年第32卷12期《湖泊水质模型SALMO在太湖梅梁湾的应用》对水质模型在太湖梅梁湾的运用做出了细致的研究。该文的湖泊水质模型SALMO并没有考虑水动力等过程对水质过程的影响，而本发明构建的水动力-水温-水质模型体现了水动力过程对水温及水质过程的影响。

《昆明理工大学学报》杂志2013年38卷第一期《洱海湖泊及湖湾三维水动力模型构建及特征分析》针对洱海湖泊及湖湾建立了三维水动力模型，对其流场、水温进行了研究。虽然该文建立的模型适合描述湖泊三维水动力过程，但是该文的模型求解方法并未涉及本发明的分裂算法思想。

《InternationalJournalforNumericalMethodsinFluids》杂志2008年56卷第六期《Athree-dimensionalnon-hydrostaticverticalboundaryfittedmodelforfree-surfaceflows》建立了三维非静压有限体积模型，用来模拟在垂向边界适应网格上的自由水面流动，而本发明是为了解决湖泊水动力-水温-水质模拟预测问题。

中国专利公告号CN102156779A公告日是2011年8月17日名称为《地下水流仿真与预测分析方法》，通过对采区地下水数据进行动态观测与采集，利用数据引擎将地下水数据集成到图形工作站中，自动构建各水层动态水位的有限元网格模型，同时确定参数分区及其相应的参数值，从而实现水层水位动态模拟以及地下水运移仿真。该专利采用的是地下水模型，地下水模型的物理机制不适用于本发明研究的湖泊水体。

综上所述，从目前国内外研究的发展趋势看，如何考虑湖泊水环境各物理过程的耦合，如何解决湖泊水环境系统中各物理过程存在显著波频快慢差异的问题，针对湖泊水环境复杂模型求解方法如何兼顾稳定性、计算精度和计算效率的问题，这些都是亟待解决的技术难题。

发明内容

针对现有湖泊水环境模型亟需解决的上述技术问题和社会发展的需求，本发明提出一种基于分裂算法的湖泊三维水动力-水温-水质模拟预测方法，基于分裂算法的思想，本发明将湖泊水环境系统中的复合物理波动过程视作由若干不同波频的单一波动过程组成的集合，遵循地球物理问题一般将低频波动过程和高频波动过程分类处理的准则，将湖泊三维水动力-水温-水质模型中的各算子按照其物理波动过程的波频快慢特性，分为低频慢过程算子和高频快过程算子两类，对其中的低频慢过程算子采用显式处理，对高频快过程算子采用隐式处理，这种分裂算法解决了湖泊水环境系统中因各物理波动过程存在显著波频差异，而难以选用合适的模型求解方法的难题。本发明方法能对模型所有算子选用最适宜其物理特征的离散处理方法，较为准确地反映湖泊水体中的动量迁移、热量传递和污染物输移等复杂的物理过程，为湖泊三维水动力-水温-水质的数值模拟计算和预测提供了有效的技术支持。

本发明提供了一种基于分裂算法的湖泊三维水动力-水温-水质模拟预测方法，包括以下步骤：

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

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

其中，步骤(1)包括以下子步骤：

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

水动力控制方程组：

连续方程 ${\∂}_{x}Hu+{\∂}_{y}Hv+{\∂}_{\mathrm{\σ}}\mathrm{\ω}+{\∂}_t\mathrm{\η}=0---\left(1\right)$

水位方程 ${\∂}_t\mathrm{\η}+{\∂}_x\left(H{\∫}_0^{1}ud\mathrm{\σ}\right)+{\∂}_y\left(H{\∫}_0^{1}vd\mathrm{\σ}\right)=0---\left(2\right)$

x动量方程 ${\∂}_{t}u+u{\∂}_{x}u+v{\∂}_{y}u+\mathrm{\ω}{\∂}_{\mathrm{\σ}}u=fv-g{\∂}_x\mathrm{\η}+{\∂}_x\left({\mathrm{\ν}}_h{\∂}_{x}u\right)+\∂y\left({\mathrm{\ν}}_h{\∂}_{y}u\right)+{\∂}_{\mathrm{\σ}}\left(\frac{{\mathrm{\ν}}_v}{H^2}{\∂}_{\mathrm{\σ}}u\right)---\left(3\right)$

y动量方程 ${\∂}_{t}v+u{\∂}_{x}v+v{\∂}_{y}v+\mathrm{\ω}{\∂}_{\mathrm{\σ}}v=-fu-g{\∂}_y\mathrm{\η}+{\∂}_x\left({\mathrm{\ν}}_h{\∂}_{x}v\right)+{\∂}_y\left({\mathrm{\ν}}_h{\∂}_{y}v\right)+{\∂}_{\mathrm{\σ}}\left(\frac{{\mathrm{\ν}}_v}{H^2}{\∂}_{\mathrm{\σ}}v\right)---\left(4\right)$

水温控制方程： ${\∂}_{t}T+u{\∂}_{x}T+v{\∂}_{y}T+\mathrm{\ω}{\∂}_{\mathrm{\σ}}T={\∂}_x\left(A_{HT}{\∂}_{x}T\right)+{\∂}_y\left(A_{HT}{\∂}_{y}T\right)+{\∂}_{\mathrm{\σ}}\left(\frac{A_{VT}}{H^2}{\∂}_{\mathrm{\σ}}T\right)---\left(5\right)$

水质控制方程： ${\∂}_{t}C+u{\∂}_{x}C+v{\∂}_{y}C+\mathrm{\ω}{\∂}_{\mathrm{\σ}}C={\∂}_x\left(A_{HC}{\∂}_{x}C\right)+{\∂}_y\left(A_{VC}{\∂}_{y}C\right)+{\∂}_{\mathrm{\σ}}\left(\frac{v_b}{H^2}{\∂}_{\mathrm{\σ}}C\right)---\left(6\right)$

式中， $\mathrm{\σ}=\frac{z+h}{H}$ 是σ坐标变换； $\mathrm{\ω}=\frac{w}{H}-\frac{\mathrm{\σ}}{H}\frac{\∂H}{\∂t}+\frac{u}{H}(\frac{\∂h}{\∂x}-\mathrm{\σ}\frac{\∂H}{\∂x})+\frac{v}{H}(\frac{\∂h}{\∂y}-\mathrm{\σ}\frac{\∂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)采用有限差分法离散模型控制方程组，对湖泊地形数据进行预处理：利用小波包变换方法进行纹理特征提取，采用反距离加权平均插值法进行空间插值获得湖底地形散点三维坐标；采用主成份分析方法对湖泊遥感影像进行光谱特征提取，获取湖泊边界的平面二维坐标；

(1.3)根据得到的湖底地形散点三维坐标和湖泊边界的平面二维坐标，将湖泊离散成若干网格单元，并对网格采用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.4)采用分裂算法将水动力控制方程组中的水平动量方程各算子按照其物理波动过程的波频快慢特性进行分类，缓行内重力波为低频慢过程，表面重力长波为高频快过程，由此将水平动量方程(3)、(4)的算子分为低频慢过程算子和高频快过程算子两类，其中，低频慢过程算子包括对流项、科氏力项和水平涡粘项，高频快过程算子包括重力梯度项和垂向涡粘项；

(1.5)对低频慢过程算子对流项、科氏力项和水平涡粘项做显式处理，求得x方向中间流速u1和y方向中间流速v1，以x方向为例：

$\begin{array}{c}u{1}_{i+1/2,j,k}=-\frac{\mathrm{\Δ}t}{2\mathrm{\Δ}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*}\mathrm{\Δ}t}{2\mathrm{\Δ}y}\left(u_{i+1/2,j+1,k}^n-u_{i+1/2,j-1,k}^n\right)\\ -\frac{{\mathrm{\ω}}_{i+1/2,j,k}^{n*}\mathrm{\Δ}t}{{\mathrm{\Δ\σ}}_k+0.5\left({\mathrm{\Δ\σ}}_{k+1}+{\mathrm{\Δ\σ}}_{k-1}\right)}\left(u_{i+1/2,j,k+1}^n-u_{i+1/2,j,k-1}^n\right)+\frac{v_h\mathrm{\Δ}t}{{\mathrm{\Δx}}^2}\left(u_{i+3/2,j,k}^n-u_{i-1/2,j,k}^n\right)\\ +\frac{v_h\mathrm{\Δ}t}{{\mathrm{\Δ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\mathrm{\Δ}t}{{\mathrm{\Δ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)$

${\mathrm{\ω}}_{i+1/2,j,k}^{n*}=\frac{1}{4}({\mathrm{\ω}}_{i,j,k+1/2}^n+{\mathrm{\ω}}_{i,j,k-1/2}^n+{\mathrm{\ω}}_{i+1,j,k+1/2}^n+{\mathrm{\ω}}_{i+1,j,k-1/2}^n)---\left(9\right)$

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

$\frac{u{2}_{i+1/2,j,k}-u{1}_{i+1/2,j,k}}{\mathrm{\Δ}t}=-g\frac{{\mathrm{\η}}_{i+1,j}^{n+1}-{\mathrm{\η}}_{i,j}^{n+1}}{\mathrm{\Δ}x}---\left(10\right)$

(1.7)根据已经求到的中间流速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.8)根据已经得到的x方向流速uⁿ⁺¹、y方向流速vⁿ⁺¹和连续方程(1)，求得垂向流速ωⁿ⁺¹，根据σ反坐标变换求得笛卡尔坐标系下垂向流速wⁿ⁺¹；

(1.9)重复步骤(1.5)-(1.8)，循环迭代计算，得到湖泊水域内不同位置和时间的三维流场u、v、w及其静水位偏移位移η，并根据静水位偏移位移η进一步求得水深H。

其中，步骤(2)包括以下子步骤：

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

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

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

${\mathrm{MT}}_{i,j,k}^{n+1}+{\mathrm{ST}}_{i,j,k+1}^{n+1}+{\mathrm{PT}}_{i,j,k-1}^{n+1}=R---\left(12\right)$

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

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

总体而言，本发明与现有技术相比，具有以下特点：

(1)基于分裂算法的思想，本发明将湖泊水环境系统中的复合物理波动过程视作由若干不同波频的单一波动过程组成的集合，遵循地球物理问题一般将低频波动过程和高频波动过程分类处理的准则，将湖泊三维水动力-水温-水质模型的各算子按照其物理波动过程的波频快慢特性，分为低频慢过程算子和高频快过程算子两类，对其中的低频慢过程算子采用显式处理，对高频快过程算子采用隐式处理；

(2)能对模型所有算子选用最适宜其物理特征的离散处理方式，解决了湖泊水环境系统中因各物理波动过程存在显著波频差异，而难以选用合适的模型求解方法的技术难题；

(3)能较准确地反映湖泊水体中的动量迁移、热量传递和污染物输移等复杂的物理过程，具有计算稳定性好、计算精度和计算效率高的特点。

附图说明

图1是本发明基于分裂算法的湖泊三维水动力-水温-水质模拟预测方法的总流程图；

图2是湖泊网格化及前处理过程图；

图3是本发明实施例中沙湖边界线示意图；

图4是本发明实施例中沙湖网格划分的示意图；

图5是本发明实施例中ArakawaC变量布置示意图；

图6是本发明实施例中湖泊某时刻的三维流场矢量图；

图7是本发明实施例中某网格水温时间序列图(垂向分为4层)；

图8是本发明实施例中湖泊某时刻COD浓度分布渲染图(垂向分为4层)；

图9是本发明实例中静水位偏移位移η的实测值与模拟值；

图10是本发明实例中水温T的实测值与模拟值。

具体实施方式

本发明以武汉市沙湖为例。沙湖是武汉市重要的城市湖泊，分为内沙湖(0.134平方公里)和外沙湖(3.197平方公里)。根据遥感影像对沙湖进行特征提取，构建沙湖三维水动力-水温-水质模型。

图1所示为本发明一种基于分裂算法的湖泊三维水动力-水温-水质模拟预测方法的总流程图，其总体思路是对湖泊水动力-水温-水质模型中的各算子根据其波频特性进行分类离散处理。包括以下步骤：

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

具体而言，本步骤包括以下子步骤：

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

水动力控制方程组：

连续方程 ${\∂}_{x}Hu+{\∂}_{y}Hv+{\∂}_{\mathrm{\σ}}\mathrm{\ω}+{\∂}_t\mathrm{\η}=0---\left(1\right)$

水位方程 ${\∂}_t\mathrm{\η}+{\∂}_x\left(H{\∫}_0^{1}ud\mathrm{\σ}\right)+{\∂}_y\left(H{\∫}_0^{1}vd\mathrm{\σ}\right)=0---\left(2\right)$

x动量方程 ${\∂}_{t}u+u{\∂}_{x}u+v{\∂}_{y}u+\mathrm{\ω}{\∂}_{\mathrm{\σ}}u=fv-g{\∂}_x\mathrm{\η}+{\∂}_x\left({\mathrm{\ν}}_h{\∂}_{x}u\right)+\∂y\left({\mathrm{\ν}}_h{\∂}_{y}u\right)+{\∂}_{\mathrm{\σ}}\left(\frac{{\mathrm{\ν}}_v}{H^2}{\∂}_{\mathrm{\σ}}u\right)---\left(3\right)$

y动量方程 ${\∂}_{t}v+u{\∂}_{x}v+v{\∂}_{y}v+\mathrm{\ω}{\∂}_{\mathrm{\σ}}v=-fu-g{\∂}_y\mathrm{\η}+{\∂}_x\left({\mathrm{\ν}}_h{\∂}_{x}v\right)+{\∂}_y\left({\mathrm{\ν}}_h{\∂}_{y}v\right)+{\∂}_{\mathrm{\σ}}\left(\frac{{\mathrm{\ν}}_v}{H^2}{\∂}_{\mathrm{\σ}}v\right)---\left(4\right)$

水温控制方程： ${\∂}_{t}T+u{\∂}_{x}T+v{\∂}_{y}T+\mathrm{\ω}{\∂}_{\mathrm{\σ}}T={\∂}_x\left(A_{HT}{\∂}_{x}T\right)+{\∂}_y\left(A_{HT}{\∂}_{y}T\right)+{\∂}_{\mathrm{\σ}}\left(\frac{A_{VT}}{H^2}{\∂}_{\mathrm{\σ}}T\right)---\left(5\right)$

水质控制方程： ${\∂}_{t}C+u{\∂}_{x}C+v{\∂}_{y}C+\mathrm{\ω}{\∂}_{\mathrm{\σ}}C={\∂}_x\left(A_{HC}{\∂}_{x}C\right)+{\∂}_y\left(A_{VC}{\∂}_{y}C\right)+{\∂}_{\mathrm{\σ}}\left(\frac{v_b}{H^2}{\∂}_{\mathrm{\σ}}C\right)---\left(6\right)$

式中， $\mathrm{\σ}=\frac{z+h}{H}$ 是σ坐标变换； $\mathrm{\ω}=\frac{w}{H}-\frac{\mathrm{\σ}}{H}\frac{\∂H}{\∂t}+\frac{u}{H}(\frac{\∂h}{\∂x}-\mathrm{\σ}\frac{\∂H}{\∂x})+\frac{v}{H}(\frac{\∂h}{\∂y}-\mathrm{\σ}\frac{\∂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为水平水质扩散系数和垂向水质扩散系数；

在本实施方案中，沙湖处于北纬29°58′，东经114°33′；水平涡黏系数V_h和垂向涡粘系数V_v分别为1.5×10^-6和1.0×10^-4；水平水温分子扩散系数A_HT和垂向温度扩散系数A_VT皆为6×10^-6；水平水质分子扩散系数A_HC和垂向水质扩散系数A_VC皆为2×10^-5；距水面10m处风速为2-5m/s；湖底摩擦系数取0.025；重力加速度g＝10m/s²；水体比热容c_p＝4.2×10³J/(kg℃)；

(1.2)如图2所示，采用有限差分法离散模型控制方程组，对沙湖地形数据进行预处理：利用小波包变换方法进行纹理特征提取，采用反距离加权平均插值法进行空间插值获得沙湖湖底地形散点三维坐标，然后将沙湖垂向分层为四层；

(1.3)通过谷歌地图获取沙湖的遥感影像，采用主成份分析方法对遥感影像进行光谱特征提取，获取沙湖边界的平面二维坐标，生成湖泊边界线L(图3)；根据湖泊边界线L的4个极点坐标，生成湖泊边界线L的一个外接矩形，将外接矩形划分为若干个方形网格。以网格中心点为标识，将网格中心与湖泊边界线L进行空间拓扑分析，中心点在湖泊边界线L内的网格保留，反之则删除，最终的网格数为305个(图4)；

(1.4)在划分好的网格上采用ArakawaC模式布置模型变量(图5)：对一个正方体网格单元，规定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.5)采用分裂算法将水动力控制方程组中的水平动量方程各算子按照其物理波动过程的波频快慢特性进行分类，缓行内重力波为低频慢过程，表面重力长波为高频快过程，由此将水平动量方程(3)、(4)的算子分为低频慢过程算子和高频快过程算子两类，其中，低频慢过程算子包括对流项、科氏力项和水平涡粘项，高频快过程算子包括重力梯度项和垂向涡粘项；

(1.6)对低频慢过程算子对流项、科氏力项和水平涡粘项做显式处理，求得x方向中间流速u1和y方向中间流速v1，以x方向为例：

$\begin{array}{c}u{1}_{i+1/2,j,k}=-\frac{\mathrm{\Δ}t}{2\mathrm{\Δ}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*}\mathrm{\Δ}t}{2\mathrm{\Δ}y}\left(u_{i+1/2,j+1,k}^n-u_{i+1/2,j-1,k}^n\right)\\ -\frac{{\mathrm{\ω}}_{i+1/2,j,k}^{n*}\mathrm{\Δ}t}{{\mathrm{\Δ\σ}}_k+0.5\left({\mathrm{\Δ\σ}}_{k+1}+{\mathrm{\Δ\σ}}_{k-1}\right)}\left(u_{i+1/2,j,k+1}^n-u_{i+1/2,j,k-1}^n\right)+\frac{v_h\mathrm{\Δ}t}{{\mathrm{\Δx}}^2}\left(u_{i+3/2,j,k}^n-u_{i-1/2,j,k}^n\right)\\ +\frac{v_h\mathrm{\Δ}t}{{\mathrm{\Δ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\mathrm{\Δ}t}{{\mathrm{\Δ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)$

${\mathrm{\ω}}_{i+1/2,j,k}^{n*}=\frac{1}{4}({\mathrm{\ω}}_{i,j,k+1/2}^n+{\mathrm{\ω}}_{i,j,k-1/2}^n+{\mathrm{\ω}}_{i+1,j,k+1/2}^n+{\mathrm{\ω}}_{i+1,j,k-1/2}^n)---\left(9\right)$

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

$\frac{u{2}_{i+1/2,j,k}-u{1}_{i+1/2,j,k}}{\mathrm{\Δ}t}=-g\frac{{\mathrm{\η}}_{i+1,j}^{n+1}-{\mathrm{\η}}_{i,j}^{n+1}}{\mathrm{\Δ}x}---\left(10\right)$

(1.8)根据已经求到的中间流速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.9)根据已经得到的x方向流速uⁿ⁺¹、y方向流速vⁿ⁺¹和连续方程(1)，求得垂向流速ωⁿ⁺¹，根据σ反坐标变换求得笛卡尔坐标系下垂向流速wⁿ⁺¹；

(1.10)重复步骤(1.6)-(1.9)，循环迭代计算，得到湖泊水域内不同位置和时间的三维流场u、v、w及其静水位偏移位移η，并根据静水位偏移位移η进一步求得水深H。

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

具体而言，本步骤包括以下子步骤：

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

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

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

${\mathrm{MT}}_{i,j,k}^{n+1}+{\mathrm{ST}}_{i,j,k+1}^{n+1}+{\mathrm{PT}}_{i,j,k-1}^{n+1}=R---\left(12\right)$

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

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

(2.5)根据求得的三维流场u、v、w、水深H、水温T和水质指标浓度C进行可视化后处理，具体为：

(a)根据某时刻的湖区三维流场的计算结果生成流场矢量图(图6)；

(b)对某个具体的网格，以时间为横坐标，以湖泊水环境某因子(流场、水温、水质指标)的值为纵坐标，生成该因子的时间序列图(图7)；

(c)根据某时刻的湖区垂向每层的水温(或水质)计算结果生成三维水温(或水质)分布渲染图(图8)。

模型验证

计算值与实测值对比

图9、图10分别是静水位偏移位移η、水温T的模拟值与实际观测数据的对比结果；表1是水质指标(COD)浓度C的模拟值与实际观测数据的对比结果：

表1沙湖水质指标COD模拟计算结果

计算稳定性与计算效率

本案例所有算法程序均运行于同一台PC(CPU:Celeron(R)E33002.50GHz；Memory:2GB,DDR2)。IDE(编译器)为visualstudio2010(c#)。采用本发明方法与传统的纯显格式处理方法对模型进行求解，计算结果表明纯显格式不能满足稳定条件，出现了计算结果发散的现象，而本发明方法在整个计算过程中保持着良好的稳定性。在网格数(305)和时间步数(100)的情况下，分别采用本发明方法与传统的纯隐格式处理方法对模型进行求解，本发明方法和纯隐格式的CPU耗时分别为1.3和3.6个机时单位，本发明方法比纯隐格式的计算效率提高了64％。

本发明方法的模拟结果与实测值吻合良好，应用实例证明本发明方法具有良好的稳定性，且计算精度和计算效率较高。综上所述，本发明提出的一种基于分裂算法的湖泊三维水动力-水温-水质模拟预测方法能较为准确地反映湖泊水体中的动量迁移、热量传递和污染物输移等复杂的物理过程，具有计算稳定性好、计算精度和计算效率较高的特点。

Claims (3)

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

水动力控制方程组：

连续方程 ${\∂}_{x}Hu+{\∂}_{y}Hv+{\∂}_{\mathrm{\σ}}\mathrm{\ω}+{\∂}_t\mathrm{\η}=0---\left(1\right)$

水位方程 ${\∂}_t\mathrm{\η}+{\∂}_x\left(H{\∫}_0^{1}ud\mathrm{\σ}\right)+{\∂}_y\left(H{\∫}_0^{1}vd\mathrm{\σ}\right)=0---\left(2\right)$

x动量方程 ${\∂}_{t}u+u{\∂}_{x}u+v{\∂}_{y}u+\mathrm{\ω}{\∂}_{\mathrm{\σ}}u=fv-g{\∂}_x\mathrm{\η}+{\∂}_x\left({\mathrm{\ν}}_h{\∂}_{x}u\right)+\∂y\left({\mathrm{\ν}}_h{\∂}_{y}u\right)+{\∂}_{\mathrm{\σ}}\left(\frac{{\mathrm{\ν}}_v}{H^2}{\∂}_{\mathrm{\σ}}u\right)---\left(3\right)$

y动量方程 ${\∂}_{t}v+u{\∂}_{x}v+v{\∂}_{y}v+\mathrm{\ω}{\∂}_{\mathrm{\σ}}v=-fu-g{\∂}_y\mathrm{\η}+{\∂}_x\left({\mathrm{\ν}}_h{\∂}_{x}v\right)+{\∂}_y\left({\mathrm{\ν}}_h{\∂}_{y}v\right)+{\∂}_{\mathrm{\σ}}\left(\frac{{\mathrm{\ν}}_v}{H^2}{\∂}_{\mathrm{\σ}}v\right)---\left(4\right)$

水温控制方程： ${\∂}_{t}T+u{\∂}_{x}T+v{\∂}_{y}T+\mathrm{\ω}{\∂}_{\mathrm{\σ}}T={\∂}_x\left(A_{HT}{\∂}_{x}T\right)+{\∂}_y\left(A_{HT}{\∂}_{y}T\right)+{\∂}_{\mathrm{\σ}}\left(\frac{A_{VT}}{H^2}{\∂}_{\mathrm{\σ}}T\right)---\left(5\right)$

水质控制方程： ${\∂}_{t}C+u{\∂}_{x}C+v{\∂}_{y}C+\mathrm{\ω}{\∂}_{\mathrm{\σ}}C={\∂}_x\left(A_{HC}{\∂}_{x}C\right)+{\∂}_y\left(A_{VC}{\∂}_{y}C\right)+{\∂}_{\mathrm{\σ}}\left(\frac{v_b}{H^2}{\∂}_{\mathrm{\σ}}C\right)---\left(6\right)$

式中，是σ坐标变换； $\mathrm{\ω}=\frac{w}{H}-\frac{\mathrm{\σ}}{H}\frac{\∂H}{\∂t}+\frac{u}{H}(\frac{\∂h}{\∂x}-\mathrm{\σ}\frac{\∂H}{\∂x})+\frac{v}{H}(\frac{\∂h}{\∂y}-\mathrm{\σ}\frac{\∂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方向为例：

$\begin{array}{c}u{1}_{i+1/2,j,k}=-\frac{\mathrm{\Δ}t}{2\mathrm{\Δ}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*}\mathrm{\Δ}t}{2\mathrm{\Δ}y}\left(u_{i+1/2,j+1,k}^n-u_{i+1/2,j-1,k}^n\right)\\ -\frac{{\mathrm{\ω}}_{i+1/2,j,k}^{n*}\mathrm{\Δ}t}{{\mathrm{\Δ\σ}}_k+0.5\left({\mathrm{\Δ\σ}}_{k+1}+{\mathrm{\Δ\σ}}_{k-1}\right)}\left(u_{i+1/2,j,k+1}^n-u_{i+1/2,j,k-1}^n\right)+\frac{v_h\mathrm{\Δ}t}{{\mathrm{\Δx}}^2}\left(u_{i+3/2,j,k}^n-u_{i-1/2,j,k}^n\right)\\ +\frac{v_h\mathrm{\Δ}t}{{\mathrm{\Δ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\mathrm{\Δ}t}{{\mathrm{\Δ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)$

${\mathrm{\ω}}_{i+1/2,j,k}^{n*}=\frac{1}{4}({\mathrm{\ω}}_{i,j,k+1/2}^n+{\mathrm{\ω}}_{i,j,k-1/2}^n+{\mathrm{\ω}}_{i+1,j,k+1/2}^n+{\mathrm{\ω}}_{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}}{\mathrm{\Δ}t}=-g\frac{{\mathrm{\η}}_{i+1,j}^{n+1}-{\mathrm{\η}}_{i,j}^{n+1}}{\mathrm{\Δ}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。

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

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

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

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

${\mathrm{MT}}_{i,j,k}^{n+1}+{\mathrm{ST}}_{i,j,k+1}^{n+1}+{\mathrm{PT}}_{i,j,k-1}^{n+1}=R---\left(12\right)$

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

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