4.2 定解条件

4.2.1 初始条件

模型初始条件是数学模型中所有与计算有关的、需要在计算初始时刻明确赋值的量。

模型要求的初始条件可分为计算初始时刻的地形条件、水流条件与工程条件。地形条件包括河底高程与糙率。初始水流条件分为两个部分，有水流流动的计算网格可根据实测资料或经过资料分析等给出水位、流速，也可按进出口水位线性插值给出水位，并令流速为0；没有水流流动的部分，即干河床部位可直接赋水位为河床高程，流速为0。

如果模型已进行过模拟，且模拟初始时刻已有模拟结果，则模型的初始水位与流速为起始计算时刻的上一次计算结果。对于用户来讲，设置初始流速和初始水位时假设模型计算的启动方式为冷启动；如果进行模型计算时选择热启动，用户设置初始计算时间后，模型软件将从上一次的模型计算输出结果中选择时间最接近的输出结果，并自动将其流速和水位设置为模型计算的初始流速和初始水位，然后进行模型计算。

对于泥沙冲淤、温度场及浓度场的计算，可令初始时刻含沙量、温升、浓度为0。

4.2.2 进出口边界条件

数值边界条件的给定，一方面要符合流动的物理性态，且在数学上和控制方程相容；另一方面，数值边界条件的规定不应使所确定的输入特征变量和未规定的边界点流动要素不合理或有较大误差。

边界处理合适与否决定数值模拟的成败。若单元某边也是区域边界，该处给定了边界条件，则需求解边界黎曼问题。此时，q_L由域内已知解确定，q_R即边界处的解，可根据当地流态（缓流或急流）选择输出特征的相容关系，并利用给定的物理边界条件加以确定。这相当于应用边界特征差分格式，不需引入冗余的数值边界条件。

边界条件的形式主要有三种：①给定水位过程；②给定流量过程；③给定水位流量关系。

在理论上，根据具体问题选用水位或流量均可，但从数值计算角度最好是上下游给定不同的变量，且对缓流可能以上游给定流量、下游给定水位为最佳组合。如上下游同时给定水位（或流量），结果中另一因变量可能有明显的误差，甚至所得恒定流解可能不唯一。对于河道恒定流模拟，一般进口给定流量，出口给定水位；对非恒定流，一般进出口均给定潮位过程。在急流入流点唯一的选择是同时给定水位和流量，在出流控制断面（闸堰、跌槛等）有固定水位流量关系时，宜规定为边界条件。

当进口给定流量或含沙量过程时，需对进口断面的流速分布和含沙量分布进行划分。

进口断面流速分布按假定的断面流速分布公式给出，模型采用以下流速分布公式：

式中，k为流速分布系数，由进口流量反求。

进口断面含沙量分布，按挟沙能力公式分配：

式中：S_l为进口第l组泥沙含沙量；S_j，l为进口断面第j点第l组泥沙含沙量；Q为进口流量；Q_j为进口断面第j点的流量；H_j为第j点的水深；U_j为第j点的流速；ω_l为第l组泥沙沉速；n为进口断面上节点总数。

对于温排水的计算，排水口边界给出排水流量、温升值（考虑热回归影响），陆地边界按绝热条件考虑。

对于污染物浓度场的计算，排污口边界给出排水流量、污染物浓度值，陆地边界按偏导数为0考虑。

4.2.3 固壁边界

模型中令固壁边界法线方向流速为0，具体处理如下：

如图4.1所示，对于岸边界节点i，设N为岸边界法线方向，T为固壁边界切线方向，θ为点i相邻的前后两点连线与x轴的夹角，u₁、v₁为在迭代过程中的前一次流速计算值，则其法向流速分量U_N和法向流速分量U_T分别为

令U_N=0，则下一次迭代的流速分量u₂、v₂分别为

4.2.4 动边界处理

1.目前常用的处理方法

目前对动边界的处理方法很多，常用的有漫水法、水边界步进法、开挖法、冻结法、切削法、窄缝法和线边界法，各种方法各有其局限性。

（1）漫水法。漫水法既考虑了处理问题的复杂性和适应性，又考虑了使数学模型计算方法上不至过于复杂，并满足流动中泥沙与水流物质守恒。在漫水法中，对于给定的网格单元，考虑水流脱离流动和恢复流动两种模式，对于当前处于流动状态的计算网格点，当计算水深小于某一临界值时，该网格脱离流动区域成为陆地边界不再参加计算，对于当前处于非流动状态的陆地网格点i，每一计算时段循环时，逐一检查周围网格点的计算状态，周围网格点上如果有流动计算点，则将其流动水面与当地网格点的河床床面比较，如该临近网格点水面高程高于当地网格点河床床面时，设其网格点序号为j，则网格i与j之间存在流动，这样，一个新的流动网格点参与整个流场计算及冲淤计算，该网格点的水位假定与周围所有流动网格点的平均水位相同，其水量来自于各参与流动的周围网格点、各自的贡献由相应的相对水深加权平均给出。由于模型各网格计算起算水深很小，因此一个网格从静到动时，对周围附近区域的流动产生的影响是很小的。

（2）水边界步进法。水边界步进法是处理动边界的一种常用方法。基本原理是：在水位变动过程中，以浅滩水边线最接近的网格节点作为边界点，因为水边线是活动的，作为边界的网格节点也是变化的。如图4.2所示，建立边界网格节点的判别式为：

当z_s≤（z_a+z_b）/2时，边界网格节点应选在较低的a点；

当z_s>（z_a+z_b）/2时，边界网格节点应选在较高的b点。

其中，z_s为水边线的竖向坐标值；z_a邻近水边线的较低网格节点的竖向坐标值；z_b邻近水边线的较高网格节点的竖向坐标值。

（3）开挖法。开挖法（周发毅，1995）将滩地开挖至可能出现的最低水位之下（见图4.3）。为使水量平衡，将岸边界向水边界移动，并增大开挖部分的糙率以求得动量上的平衡。这种方法一般只适用于潮滩问题，而且滩地面积占整个计算海区较小的情况，而对于水位变化小、坡度较缓的地形，这种处理则易失真。

（4）冻结法。冻结法（程文辉等，1988；魏文礼等，1994）根据水深判断网格单元是否露出水面，对露滩单元取糙率系数为一接近无穷大的正数，这时糙率和潮位都布置在网格中心点上，使单元四周的流速为一趋于零的无穷小量，这样处理的结果相当于使该单元的潮位在计算过程中被“冻结”不变（见图4.4）。这种方法适用于宽浅、坡度较平坦的露滩问题，而对潮滩相间的海岸河口水域，则因水量和动量的过分“冻结”而失真。方春明（1991）在计算河道平面二维恒定水流运动方程时采用冻结法处理岸边界。

（5）切削法。切削法（夏云峰，1993）对露滩单元并不“冻结”水位，而是引入一个富裕水深（相当于岸滩上存在很薄的水层）以保证计算过程的完整，相当于将原始地形切削降低，而一旦判断实际水深大于富裕水深时，即恢复原始地形。该法也有开挖法类似的局限性。辛文杰（1993）针对这种限制，提出了一种修正措施，有一定的效果，但这种修正较繁且带有一定的经验性。

（6）窄缝法。窄缝法（何少苓等，1986）假想在岸滩的每个网格上存在一条很窄的缝隙，它的深度与岸滩前的水深一致。根据水量平衡原理将窄缝内的水量平铺到岸滩上，相当于将岸滩前的水域延伸到岸滩内，从而可以把计算网格点设在岸滩上。这种方法只适用于岸边滩的露滩问题，对水域内浅滩的显露则不适用，而且连续方程的形式需作一定的修正。

（7）线边界法。线边界法（林珲等，2000）为一种新方法，也是一种较先进的动边界技术。该方法不像传统方法以网格（面）为考虑对象，而以网格的边线为对象进行露滩处理，这在物理上更为真实。而且还涵盖了网格陆露的现象和工程边界的处理（潜堤、潜坝时的边界问题），具有较广的适用性。

以上方法在数值模拟中都取得了一定的计算效果，但在模拟溃坝后干河床洪水演进过程中，直接采用上述处理方法有一定的局限性。如“切削法”中，由于方程中存在水位的偏导项且初始时刻水位一般取为与河床高程相同，因此当河道地形起伏时，对于洪水还没有传播到的位置会因水位偏导项不为0而产生水流虚假流动；“冻结法”中，对于下游为干河床时，会引起水流的停滞不动；“窄缝法”中窄缝的宽度难以确定，在模拟干河床洪水演进中计算误差较大等。

2.本书提出的新方法

目前，关于干湿边界的处理，目前已提出许多的方法，如“漫水法”、“切削法”、“冻结法”、“最小水深法”等，但各种方法均有其局限性。为此，本书提出了一种适合于干湿边界处理的新方法，并在计算中取得了较好的模拟效果。其主要思路为：首先采用最小水深法区分干湿节点；然后对干节点周围进行查找，若有湿节点且湿节点水位大于该干节点河床高程，即认为该干节点为“潜在”湿节点，可参与运算，否则为实际干节点，不参与运算。

与传统的“最小水深法”不同，本方法设置两个“最小水深”，即h_dry和h_wet，若某节点水深h<h_dry，将周围节点的水位减去本节点的河底高程，得到周围各点相对本节点的水深h′_i。若h′_i均小于h_dry，则表示该节点周围均为干河床，则该节点干出，即h=h_dry，该节点作冻结处理；若至少有一个节点h′_i＞h_wet，则令h=h_wet，即该节点为潜在的湿节点，可参与运算。在通常情况，可取h_dry=0.005m，h_wet=0.05～0.1m，如图4.5所示。

以上处理方法的优点在于对于洪水流经区域及水陆交接干节点可直接参与计算，而对于洪水还没有流到的干河床，由于采用了“冻结”处理，在计算过程中不会因为地形起伏而产生水流虚假流动。

3.概化溃坝模型检验

（1）部分溃坝模型实验。实际溃坝水流实测资料一般难以得到，因此本书引用Fraccrollo等进行的溃坝物理模型实验成果对模型进行检验。

物理模型为1m×2m的矩形容器，在矩形容器的长边中央有一个宽0.4m缺口，并用一矩形挡板挡住。该挡板能够在极短时间内开启，构成一个对称型部分溃坝模型（见图4.6），从而瞬间形成向下游推进的正波和向上游传播的负波。

在挡板上下游分别布置若干测点，分别采用压力表和波高计测量各点水位过程。

（2）数学模型计算条件：

计算域及网格划分：选择挡板上下游长3m、宽2m的区域作为数学模型计算区域，划分为75×50个正方形网格，网格边长为0.04m。

初始条件：初始时上游水深为0.6m，根据经验下游水深不能为0，计算中取为0.001m，初始时整个流场流速为0。

边界条件：挡板下游三面为开边界，容器内除挡板处为开边界外，其余为闭边界，容器内水量一定。

主要参系数：计算时间步长为1ms，不计底坡和底部阻尼。

（3）数学模型计算结果分析。图4.7给出溃坝初期（t=0.1s）和溃坝一定时间（t=0.6s）水位三维演示图。从计算成果来看，溃坝发生后约0.6s时水流到达计算域左右开边界，1.2s左右到达计算域最下游开边界，与物理模型实验观测结果一致，由此从定性说明数学模型计算结果较好反映了实际溃坝水流演进情况。

限于篇幅，图4.8给出了部分典型测点计算水位过程与实验结果对比情况。从中可见，计算和实测水位过程除在口门流动梯度变化较大的区域（如测点0）有一定的误差外，其他区域吻合较好，说明该模型能较好模拟溃坝波的演进过程。从中还可看出，随着上游容器内水流体积减少，坝上游水位不断下降，坝下游水位先升高后降低，且离口门越远，水位波动越小。

从计算与实验结果相比较可知，在溃口初期，数值计算成果与模型实验结果存在一定误差，尤其是在口门附近。主要原因是二维模型代替三维模型时浅水方程的静水压力假定，计算时忽略实际存在的底部摩阻、采用小水深以及实验时挡板开启方式和时间等综合因素不可避免带来的差异，但总体而言该模型能较好模拟溃坝波的运动特性。

4.天然干河床洪水演进检验

为反映模型对天然河道干河床的适应性，选取一天然河道进行洪水演进模拟计算。计算河段长约18km，采用矩形网格进行划分，网格边长为100m，时间步长取为1s。进口给定流量和水位，假定下游为干河床。图4.9给出了不同时刻洪水演进图（dt=600s），从中可见，模型能较好模拟洪水在干河床上的演进过程，说明提出了干河床及动边界的处理方法，能较好适应洪水演进过程中干湿边界的调整。

4.2.5 内部边界

在模型中水工建筑物处的质量和动量计算，也可在有限体积法的框架之下进行的，但法矢量的计算方法并非利用一般的差分格式。

内部边界处的质量通量计算一般采用水工建筑物过水经验公式。一般的水工建筑物为闸、堰、桥孔、管道、堤等。边界为堰，则采用相应的堰流经验公式；边界为闸门式溢洪道和桥，则视流态的不同采用相应堰流或孔流经验公式。边界为堤，但水位低于堤顶时，作为固壁边界处理；如水位高于堤顶时，则采用堰流经验公式计算流量，此时堰的宽度为堤长。

涵闸的存在，使本来连续的天然水流改变了走势，也改变了天然河网的水力连接。过闸水流有三种情况：关闸、自由出流、淹没出流。对于关闸的情况，闸上、闸下河道不再联系，Q=0；该条件可作为这两条河道的边界条件计算。开闸时，通过以下的计算可求得过闸流量。因为闸没有容积，即流入和流出的流量相等。在进行水量平衡计算时，该流量即为与闸相连的断面之流量。对于闸上、闸下的河道，是通过闸联系起来的。淹没出流和自由出流在任何时刻均可以发生，亦可以相互转换，由该时刻的水力条件决定。

过闸流量的计算：

开闸时，根据闸的上、下游水位，可能有自由出流和淹没出流两种情况，由水力学知识，其判别条件为闸后水深与闸前水深的比值，大于0.8时为淹没出流，否则为自由出流。其自由出流和淹没出流过闸流量Q可按宽顶堰公式计算：

或

上二式中：Q为过闸流量，m³/s；m为自由出流系数，一般取0.32～0.385；φ为淹没出流系数，一般取1.0～1.18；B为闸门开启总宽度，m；H₀为闸上水深，m；z_u为闸上游水位，m；z_d为闸下游水位，m。

式（4.21）适用于自由出流，式（4.22）适用于淹没出流。为了计算方便及两种出流之间的衔接连续，将自由出流公式写成与淹没出流公式统一的形式：

当δ≤σ时

当δ≥σ时

其中

式中，为判断流态的标准，可根据自由出流系数m与淹没出流系数φ求得。设临界流时的下游槛顶以上水深为σH₀，临界流时用式（4.21）或式（4.22）计算的流量相等，故

内部边界处的动量通量，可用下式计算：

式中：h为水深；u_n、v_t分别为边界法向和切向方向的流速，其中u_n=ucosφ+vsinφ；v_t=vcosφ-usinφ，u、v分别为x和y方向的流速。

4.2.6 工程边界

在河道数值模拟中，有时对某些萎缩河道进行封堵处理，枯水为堤，大水为堰。另外，随着经济的发展，在河道中修建了大量的桥梁、码头工程和航道整治工程（如丁坝、顺坝、潜坝、护滩带等）。往往许多结构物的尺寸小于网格尺度，不能作为网格单元进行简单概化。

以拦河坝为例，在已往的模型研究中，不得不以拦河坝为界分成坝上和坝下两个模型进行研究。为在同一个模型中模拟结构物的影响，首先需根据结构物上下游的水位条件进行结构物的水力学计算，然后将下泄流量提供给结构物下游河段进行水动力学计算。

下面以宽顶堰为例，介绍结构物与河道水流的嵌套计算。宽顶堰布置在河道中影响到的单元界面如图4.10（a）所示，不同水位下堰顶的宽度如图4.10（b）所示。

对于受影响的单元，需补充堰流方程，采用线性化公式表示为

式中：φ为流速系数；B为堰宽；H₁、H₂分别为上下游水位与堰顶之差；z_u、z_d分别为上下游水位。