1.2 河口区潮流数学模型
珠江三角洲河道网状结构复杂，河道狭长，适合采用一维河网数学模型进行模拟，河口湾水域开阔，水深不大，适合采用二维数学模型模拟，而这两部分水体又是相互贯通、整体不可分割的，开发适用于珠江河口的一维、二维联解数学模型尤显必要。
1.2.1 模型范围
模拟范围包括河网主要水道100多条，及珠江八大出海口门与伶仃洋、黄茅海两大河口湾。综合考虑河口水系特点、水文（潮位）站分布与研究范围需要，模型上边界取至各江河控制水文（位）站：西江-马口站，北江-三水站，东江-博罗站，流溪河-老鸦岗站，增江-麒麟嘴站，潭江-石嘴站；下游开边界南取至外海25m等深线附近，即担杆列岛一线，西至镇海湾，东至大鹏湾港口站；一维、二维联解点设在出海口门控制断面：虎门大虎站，蕉门-南沙站，洪奇门-冯马庙站，横门-横门站，磨刀门-灯笼山，鸡啼门-黄金，虎跳门-西炮台，崖门-官冲，见表1.2-1及图1.1-1。
表1.2-1 整体数学模型范围及外边界选定
1.2.2 一维潮流模型控制方程
1.方程组
连续方程：
动量方程：
式中：Z为断面平均水位，m；Q为断面流量，m3/s；A为过水面积，m2；B为水面宽度，m；β为动量修正系数；g为重力加速度，m/s2；Sf为摩阻坡降，采用曼宁公式计算，Sf=g/C2，C=h1/6/n；h为水深，m；q为旁侧入流，m3/s；x为距离，m；t为时间，s。
2.汊口连接条件
网河区内汊口点是相关支流汇入或流出点，汊口点水位、流量要满足下列连接条件：
流量连接条件：
水位连接条件：
式中：Qi为流入汊口节点第i条支流流量，流入为正，流出为负；Zi，j等表示汊口节点第i条支流输第j号断面的平均水位；m为流入汊口的支流数。
3.方程离散与求解
方程离散采用四点加权Preissmann固定网格隐式差分格式。沿河段布置变量，模型采用水位、流量同一网格布置法（即水位、流量点布置在同一点），如图1.2-1所示。
图1.2-1 河段变量布置
以S代表流量Q和水位Z，则S在河段Δx、时段Δt内的加权平均量及相应偏导数表示为
式中：θ为差分系数。
经推导，一维潮流连续方程式（1.2-1）离散为
式中：a1、b1、c1、d1、e1为差分方程的已知系数。
一维潮流动量方程式（1.2-2）离散为
式中：a2、b2、c2、d2、e2为差分方程的已知系数。
略去未知数的上标n+1，就式（1.2-6）和式（1.2-7）求解变量Zi与Qi得
其中
设全河段有i个计算断面，未知量共2i个。每一计算河段上均有如式（1.2-8）的一对方程，则共有2（i-1）个方程，加上两个边界方程Zi=Z（t），Qi=Q（t）或Z=f（Q），其中Z（t）、Q（t）、f（Q）为边界上的已知函数，构成有定解的方程组，用高斯消去法，使系数化为上三角矩阵，然后回代求出各变量。
1.2.3 二维潮流模型控制方程
1.方程组
直角坐标系下潮流连续方程
直角坐标系下u、v向潮流动量方程
其中
式中：x、y为直角坐标；t为时间；z为水位，m；h为水深，m；u、v分别为垂线平均流速在x、y方向的分量，m/s；c为谢才系数；νt为紊动黏性系数；u*为摩阻流速，m/s；f为摩阻系数，a为综合系数，与河道形态及水流条件等有关；g为重力加速度，m/s2。
2.曲线坐标转换技术
考虑模拟区域边界、地形的复杂性，采用直角坐标系，不但网格剖分困难，边界也难以精确拟合，计算工作量也大。为此，采用坐标变换技术，将直角坐标下的变量方程转换为贴体曲线坐标系下的方程，即将物理坐标内计算域D的网格变换到计算域D*的网格，如图1.2-2所示。对于域D，变换函数可写为
式中：ξ、η为独立变量；x、y为因变量。
图1.2-2 坐标变换示意
两种坐标系中变量之间的转换关系
其中
式（1.2-13）即为生成正交曲线网格的方程和控制函数。
经坐标变换为贴体曲线坐标系下的潮流控制方程组为
其中
式中：Cξ、Cη为拉梅系数。
为了进行工程结构等小尺度模拟，引进通度系数θ改进后二维模型连续方程
式中：θc对应离散单元的面通度系数，为网格中能够被流体通过的面积（网格面积减去网格中固体或障碍物的面积）与整个网格面积之比，定义在网格中心；θζ、θη分别为对应于离散单元的ζ、η方向线通度系数，为该方向上能够被流体通过的网格长度与该网格总长之比，定义在网格边界上；通度系数根据结构占据网格大小，取值在0与1之间，占满网格取0，不占网格取1。
3.离散方程
采用交替方向隐式法（ADI法），基本方法是（图1.2-3）：设Δt、Δx、Δy分别为时间步长和x、y方向空间步长，n、i、j分别为时层数和x、y的步长数，h为实际水深，z为水位；在x-y平面上采用交错网格，并给定各变量（h，u，v，z）的计算点；在时间上采用将Δt分成两个半步长，计算采用隐、显格式交替隐、显进行，在nΔt→（n+1/2）Δt半步长上用隐格式离散连续方程和x方向上的动量方程，并用追赶法求得（n+1/2）Δt时层上的z和u，对y方向上的动量方程则用显格式离散，并求得（n+1/2）Δt时层上的v，然后在（n+1/2）Δt→（n+1）Δt半步长上用隐格式离散连续方程和y方向上的动量方程，并用追赶法求得（n+1）Δt时层上的z和v，对x方向上的动量方程则用显格式离散，并求得（n+1）Δt时层上的u。
图1.2-3 ADI法网格分层与变量布置示意
经推导，前半个时间步长，二维潮流连续方程式（1.2-14）离散如下：
其中
二维潮流动量方程式（1.2-15）离散为
其中
经推导，后半个时间步长，二维潮流连续方程式（1.2-14）离散如下：
其中
二维潮流动量方程式（1.2-16）离散为
其中
将方程式（1.2-18）和式（1.2-19）联立可写成以下简单形式：
式中：i=1，2，…，n。
式（1.2-22）和式（1.2-23）可用追赶法求得un+1/2和zn+1/2。
同样，将方程式（1.2-20）和式（1.2-21）式联立可写成以下简单形式：
式中：i=1，2，…，n。
式（1.2-24）和式（1.2-25）可用追赶法求得vn+1和zn+1。
1.2.4 一维、二维模型联解条件
整体模型在一维、二维模型连接点上的水位、流量满足以下连接条件：
水位连接条件：
流量连接条件：
式中：Z1为一维模型在内边界断面上的水位，m；Z2为二维模型在内边界上各节点的平均水位，m；Q1为一维模型在一维、二维模型连接断面上的流量，m3/s；Uξ为二维模型在一维、二维模型连接断面法向上的流速，m/s。
一维、二维模型联解基本思想是一维模型以流量传递给二维模型，二维模型以水位传递给一维模型。联解条件推导如下：
假设一维模型下边界A断面的流量为Q1，过水面积为A1，水力半径为R1，谢才系数为C1，水力坡降为J1。按谢才公式，有
变换为如下形式：
A断面同时又是二维模型的上边界，是一维、二维模型的连接断面，A断面刚好是二维模型网格中的第i1行，一维模型和二维模型的网格搭接情况如图1.2-4所示。
图1.2-4 一维模型和二维模型的网格搭接示意
假设A断面上第j条垂线的流速为Ui1-1/2，j，谢才系数为Ci1，j，水力坡降为J2，水力半径为Ri1，j，按谢才公式，有
变换为如下形式：
对于同一断面，J1=J2，所以
则
对于A断面，又因
所以
故
又因珠江河口区水较深，水力半径R可近似用水深h来代替，则上式变为
式（1.2-28）即为一维、二维模型水力要素的传递公式。
根据一维、二维连接方程式（1.2-26）、式（1.2-27）通过消元可导出二维潮流模型的控制方程
左边界为流量边界条件时σI1，J=0，左边界为水位边界条件时λI1，J=0。
通过消元可导出二维潮流模型对一维潮流模型的控制方程：
将方程式（1.2-30）作为一维模型的边点方程，通过河网非恒定流三级联解即可解出一维、二维模型连接断面上的水位及流量，利用连接断面上的水位及流量，分别回代给一维、二维模型即可计算所有计算点上的物理量。
1.2.5 整体模型建立与验证
考虑模拟范围太大，建立与验证困难，先将模型分为一维、二维分别建立与验证，再实现联解。
1.建模地形
珠江三角洲河网主要采用1999年水利部门联合测量的河道地形资料。河口湾伶仃洋浅海区：采用2000年实测地形图 （1∶10000）及1998年实测地海图 （1∶25000）；黄茅海及鸡蹄门浅海区：采用1989年实测地形图 （1∶10000）；磨刀门及澳门浅海区：采用1999年实测地形图 （1∶10000）、2000年实测地形图 （1∶10000）、2001年实测海图（1∶10000）及1998年实测海图 （1∶25000）；香港水域及外海：现有最新海图。
2.一维模型验证
模型共布置226个汊口，97条水道，352个河段，3449个河道断面，31个边界（包括河网区水闸内边界）。水文资料主要采用珠江三角洲河网“1998·6”洪水、“1999·7”中洪水、“2001·2”枯水典型组合，三组均包括了河口大、中、小潮时段。
经反复调试计算，“1999·7”中洪水组合，在102个潮（水）位验证站点中，最高潮（水）位计算误差在0.1m以内的有94个站点，占整个验证站的92.2%，计算误差在0.2m范围内的测流断面百分比为98.1%；同样，在102个潮（水）位验证站点，最低潮（水）位计算误差在0.1m范围内的测流断面百分比为90.2%，计算误差在0.2m范围内的测流断面百分比为98.1%；符合《水利工程水利计算规范》（SL 104—95）感潮河段非恒定流计算的验证标准相位相对误差应小于1h，峰谷值误差应不大于0.1～0.2m，单站检验合格率应为85%以上的要求。各测流断面计算流量过程与实测过程基本吻合。
“2001·2”枯水组合，48个潮（水）位验证站点最高潮（水）位误差绝对值小于0.2m的有46个站点，占整个验证站的95.8%；最低潮（水）位误差绝对值小于0.2m的有44个站点，占整个验证站的91.7%；符合《水利工程水利计算规范》（SL 104—95）感潮河段非恒定流计算的验证标准相位相对误差应小于1h，峰谷值误差应不大于0.1～0.2m，单站检验合格率应为85%以上的要求。各测流断面计算流量过程与实测过程基本吻合。
“1998·6”洪水组合，53个水位站点洪峰水位计算与实测误差均小于0.1m，单站检验合格率达100%，符合《水利工程水利计算规范》（SL 104—95）感潮河段非恒定流计算的验证标准相位相对误差应小于1h，单站检验合格率应为85%以上的要求。
经率定河网各河道糙率在0.014～0.030之间，糙率值较大的水道是：海州水道、虎跳门水道的横坑段、潭江、容桂水道、陈村支涌，其糙率值为0.030，口门附近水道的糙率值在0.020左右。模型主要河道糙率值网河区上部河道糙率值较大，近河口则较小，符合河口区河道糙率沿程变化的自然变化规律。
3.二维模型验证
二维模型水域较大，分为伶仃洋、磨刀门和黄茅海等三个区域重点验证，模拟水域面积约124亿m2，共布置网格950×844个。选取“1991·12”枯水、“1992·7”中水、“1998·6”洪水等几组散点资料较多的组合，经模型调试计算，流速、流向验证结果良好，与实测过程相符，流场分布基本合理，计算误差也达到了有关规范规程要求。
4.整体联解模型验证
收集河网、河口湾有同步资料的水文组合“1998·6”洪水组合进行模拟计算，计算成果合理，联解成功，可基本反映河网、河口湾水动力的整体联动性。部分验证成果见图1.2-5～图1.2-7。
图1.2-5 珠江河口区整体模型一维潮（水）位验证成果（“1998·6”洪水组合）
图1.2-6 珠江河口区整体模型二维潮（水）位验证成果（“1998·6”洪水组合）
图1.2-7 珠江河口区整体模型二维流速、流向验证成果（“1998·6”洪水组合）