渤海流场的数学模拟和应用

1、渤海流场的数学模拟和应用王万战 余欣 杨明（黄河水利科学研究院， 郑州 450003）摘要 利用渤海实测资料率定了渤海平面二维潮流模型，模拟了渤海潮流场。进行了水流参数灵敏度分析，发现涡粘系数变化对水位和流速影响不大，而糙率影响较大。数学模拟的成果表明，与河流水流相比，渤海潮流场的基本特点是，水流不仅存在“水向低处流”的过程、而且还存在“水向高处流”的过程，两个过程相互转化，其实质是动能和势能“此涨彼消”、相互转换的过程。与相对狭长海峡内的潮波、潮流比，渤海流场的特点是，渤海潮波、潮流是旋转型，而不是前进波型。这些特点也是黄河口实体模型水流控制技术的难点。把此数学模型应用于黄河口滨海区地形测验

2、中，发现如果使用数学模型给出各测点逐时水位，可减少误差±0.5m, 消除用常规方法造成的不合实际的“深海淤积现象”。水动力模拟成果表明，用把黄河口导向深海、或加大黄河口入海流量，都不能解决入海泥沙淤积口门附近的问题。黄河口治理的出路是，在于控制适当的入海水沙条件和用较大的容沙体积换取较长的行河时间。关键词 渤海；流场；M2分潮；二维潮流数学模型黄河口滨海区及渤海深海流场特性是影响黄河口水沙运动的基础。本文利用平面二维模型的水动力模块模拟渤海流场水位、流速大小和方向的动态变化；在此基础上，分析了黄河口实体模型的水流控制的难点所在，分析了黄河口滨海地形测验资料整编传统方法的缺点关键技术，

3、初步分析了用向深海延伸黄河口、治理黄河口的方法是否可行。1 基本方程基于水动力学和泥沙运动力学的河口模型能够计算在河流来水来沙、风、浪、潮汐、盐度、温度等因素共同作用下河口海域的水沙运动和地形演变。其基本方程为水流连续方程和动量方程。 X方向动量方程：Y方向动量方程：式中，为水位(m)，d为河底高程(m)，h为水深（d）（m），p、q 分别为 x,y方向的单宽流量(m3/s/m)，即p=hu,q=vh,u、v分别为流速在x、y方向的分量，C(x，y)谢才系数(m1/2/s)，与曼宁系数的关系为 c=n-1h1/6,g为重力加速度(m/s2)，f(v)为风摩擦系数，V、Vx、Vy 分别为风速及其

4、在x,y方向的分量(m/s)， 科氏力系数s-1， pa 为 大气压力(kg/m/s2),w水的密度(kg/m3),x，y为距离(m),t 为时间(s)，E为水流紊动粘滞系数，由Smagorinsky公式计算，由于本研究重点在于分析渤海潮流场的基本特性，所以暂时不考虑潮汐以外的其他因素如风、浪、温度、盐度和河流入汇的影响。2 模拟的基本条件和率定本研究模拟范围为黄河口以外、渤海海峡以西的整个渤海。黄河口附近海域地形采用黄河水利委员观测的1980年、1981年的滨海区地形，其他海域地形采用黄河水利委员设计院测绘总队和中国人民解放军4210工厂在1998年共同编汇的1：1百万的渤海地形图。计算网格

5、为边长5000m的正方型网格。边界条件：利用渤海海峡的北长山、北隍城、羊头洼三站潮汐调和常数，用调和分析模型预报出三站的潮汐过程，做空间线性内差作为模型的边界条件。模拟周期为1981年7月1日8月22日，模拟时间步长为30秒。模型的率定：需要率定的参数主要为Smagorinsky公式子涡扩散系数（Cs）和反映河床糙率的曼宁系数（n）。率定的结果为，曼宁系数为0.0125s/m1/3、Smagorinsky系数为0.5时，模拟结果，如M2分潮无潮点位置、渤海M2分潮潮差空间分布（图1）、渤海湾附近海域潮流逆时针方向旋转、K1、O1、S2等分潮潮差远小于M2分潮潮差等指标与实测资料相符。（a）模型

6、计算的M2分潮潮差(b)基于实测资料的M2分潮潮差图1 实测和计算的渤海M2分潮潮差分布3 模型参数的灵敏度分析进行模型参数灵敏度分析的目的是研究每个参数对水流的影响程度，本文的可调参数为子涡扩散系数（Cs）和河床糙率。由图2可知，在曼宁系数不变的情况下，随着Cs增大，潮位、流速降低，但是变幅很小。由图3可知，在涡粘系数一定时，随着曼宁系数的增加，潮差和流速振幅明显变小。图2 曼宁系数为0.0125s/m3、不同Cs条件下的水位图3 Cs=0.5、不同曼宁系数(n)条件下的水位、流速过程4渤海流场特征4.1 水既向低水位处流、也向高水位处流常见的河流水流特征是水向低水位处流，而在潮汐影响的水域

7、，水流既可向低水位处流，也可向高水位处流，简述如下。图4表示渤海湾中部某处的流速（流速方向以北向为基线，顺时针为正，单位为弧度），图5各图中等值线上的数字表示水位，箭头方向表示流向，箭头长短表示流速大小。图4 计算渤海湾中部某处的水位、流速、方向9：00时渤海湾处于较高水位（图4），由图5可知，此时渤海海峡水位较低，从渤海湾到渤海海峡水面比降向东，水向东流，呈现水向低水位处流的特征， 而此时莱州湾水流却是向高水位处流。其后，由于水流流出渤海湾，渤海湾水位逐渐降低、流速变大，与此同时，渤海海峡附近由于水流流入，潮位逐渐升高、流速变小。(a)(b)(c)(d）图5（a）(b) (c) (d) 计算

8、的渤海水位和流速12：20时渤海湾水位接近低潮位、但水流继续东流，呈现水向高处流的特征，流速接近最大，此时渤海海峡附近水位接近高潮位，流速接近于零；其后至至14：20，渤海湾水位继续降低、流速开始降低，及至14：20时渤海湾水流流速接近于零，渤海海峡附近水位较高，水向低水位处流（即向西流）；15：20时，渤海湾水位由低逐渐升高、水流向西，呈现水向低水位处流的特征。与渤海湾水流相似，莱州湾、辽东湾水流也具有水向低水位处流和水向高水位处流两种相互转化的过程。这是典型的水流动能和势能相互转化的过程。此特征表明，在某一时段内，渤海水面比降是水流的主动驱动力，而在其他时段，其是水流的阻力此特征是一般河道

9、水流所没有的。这是黄河口实体模型试验水流不同于其他单向重力流的特征之一。42 旋转流和往复流图4最大流速出现在半潮位时刻，而最小流速则出现在高、低潮位附近，表明渤海潮波为旋转型。不同的地点流速矢量旋转的方向不同（顺时针、或逆时针）例如，在渤海湾为逆时针旋转。越近岸边，潮流椭圆长轴相对较长，且平行于海岸线，表现为越为明显的往复流。5应用5.1 黄河口实体模型水流控制的难点渤海水流既朝低水位处流，也朝高水位处流，同时，渤海水流为旋转流，神仙沟口外存在M2无潮点等，都是黄河口实体模型水流控制的难点所在。5.2 黄河口滨海地形测验中潮位的确定测验断面某点的高程Hb，是通过水位Z和该点瞬时实测水深（h）

10、得到的，即Hb=Z-h (图6)。瞬时水位 Z 海底高程 Hb h高程基面 图6 海底高程的求法问题是，如何求出水位？传统的方法是，假定任意测点水位等于附近岸边测点水位。然而，这个假定常常造成黄河口滨海区深水区出现严重的淤积显然不符合实际。用数学模型检查这样的假定。先看黄河口外滨海测验常用的35个断面中的第12断面（图7）。点绘第12断面两端点瞬时水位，计算两端点水位差，如图8所示，可见测验断面两端点水位差不是像假定那样为零，而是-0.4m 到 +0.6m。图7渤海湾滨海区地形测验断面分布图8 第12断面两端点处水位差同理，可计算处其他34断面两端点的水位差（图9），可见35个断面两端点水位差

11、约为-0.4m - +0.5m. 由上述可见，有必要使用数学模型为黄河口滨海区地形资料整编提供所需的水位。图9 黄河口滨海区每个断面两端点水位差5.3 筛选黄河口治理战略和方案如何治理黄河口一直是大家关心的大事之一，不同的专家提出了各种不同的观点。数学模型是筛选治理方案的经济手段之一。 黄河口治理观点1研究：黄河口向深海延伸有专家提出，尽量向深海延伸黄河口，希冀突出的沙嘴造成的流速增大能够把泥沙输送更远。下面用数学模型来论证此观点是否可行。考虑两种情形：黄河口向东延伸、以及黄河口向西延伸。首先看黄河口向东延伸的情形。用数学模型计算出当黄河口延伸0 km、10km、20km、30km、40km、

12、50km、60km、70km、80km、和90km时各点的流速，取流速最大值，然后用延伸后的流速最大值减去延伸前的流速最大值。图10为当黄河口延伸30 km时的流速最大值减去延伸20 km时流速最大值所得的差值分布图。可见，黄河口延伸的确可以增加沙嘴附近的流速，但是流速明显增加的区域面积较小。. 这些结果表明，即使黄河口延伸造成口门局部流场增强，输送到黄河口的泥沙的较粗沙仍不会被水流带到较远的地方，仍将淤积在口门附近。模拟结果初步否定了此种治河观点。图10 黄河口延伸30 km、20km时流速最大值差值几十年来黄河口研究表明，尽管专家们对黄河口淤积延伸对黄河下游河道演变影响有不同的认识，但是对

13、黄河口淤积延伸必然反馈抬高其上游河道的看法还是一致的。因此，用把黄河口延伸到深海、希冀借助突出沙嘴造成流速增大把泥沙带向深海的方法是既不可行、如果实施此方法，会加快黄河口上游河道河床抬升的速率。黄河口治理观点2研究：加大入海流量，希望大流量把泥沙带到深海有专家提出，黄河口治理可以通过加大入海流量，希冀大流量把泥沙带到深海。用数学模型模拟了清水沟口门大流量（5000 m3/s）、 出流方向为0度(N)、45度(NE)、90度(E)、135度(SE)、180度(S)时黄河口口外的流场。这些流场与清水沟口门流量（0 m3/s）时的流场相比较，可看出当清水沟口门大流量时的影响范围（图11），其大致为口

14、门附近40km (南北向)×30km(向东)。这个范围还是比较小的，无法把泥沙输送更远。 黄河口治理的根本出路：用堆沙空间换行河时间黄河口来沙多时，河口延伸长，反之，来沙少时则短钓口河1976年后停止行河后，钓口河附近海岸大量蚀退。这些事实表明，入海沙量是决定黄河口地貌演变的关键因素，因此，如何控制黄河口入海水沙应是黄河口治理的关键。 图11 黄河口流路向东时不同出流角度对向东西向流速的影响 (河口流量5000 m3/s)从理想的角度，适量（既不太多、也不太少）的泥沙既能防止黄河口河道萎缩、也能使河口海岸处于冲淤平衡。但是，实现起来可能不太容易：需要为黄河口建立一个或多个分水分沙系统

15、,还需要动态调度黄河口水沙过程。黄河口治理的出路是，在尽量控制黄河口来沙量的基础上，先相对稳定地使用一入海流路，然后有计划地摆动到另一流路，由海洋动力蚀退非行河流路海岸，待蚀推达到一定程度后，伺机再使用。即用有限的容沙体积换取黄河口流路的行河时间。至于何时实施人工流路摆动等，需要大量的数学模拟工作。6结论（1）平面二维数学模型能够再现渤海流场的基本特性。（2）在平面二维潮流参数中，子涡扩散系数的变化引起的水流变化较小，而糙率的变化引起的水流变化较大。（3） 渤海流场的基本特性是：渤海潮波、潮流是旋转型的，水流既可由高水位向低水位流，也可由低水位向高水位流。这些特点决定了黄河口实体模型水流控制的

16、难点所在。（4）黄河口滨海区水位不是水平的，岸边与深海处水位相差约-0.4m-+0.5 m；至今黄河口濒海区地形整编时使用的假定（深水区水位等于岸边水位）造成深水区出现虚假的淤积现象，可由数学模型提供动态的水位加以消除。（5）渤海水动力数学模型模拟表明，藉把黄河口牵引到深水区或增加黄河口入海流量的方法都只能增加口门附近局部的流场强度，对此小范围之外的流场影响不大，因此，这些方法无法把黄河口入海泥沙的较粗部分输沙到深海。（6）黄河口治理的根本出路是，在控制黄河口的水沙量的基础上，同时（b）给黄河口流路留以摆动的空间，用大范围的容沙空间换取相对较长的行河时间。 Numerical Modeling

17、 of the Flow in Bo Sea and its ApplicationsWang Wan Zhan Yu Xin Yang Ming（ Yellow River Institute of Hydraulic Research, Zhengzhou, 450003,Abstract: With a 2D depth-integrated hydrodynamic model, we did parameters analysis and simulation of the flow in the Bo Sea. On the basis we have found some app

18、lications in defining what the challenging problems in the physical modeling, compiling nearshore bathymetric data , and exploring feasibility of each of the strategies for training the Yellow River estuary. Keywords:Bo Sea; Flow Field; M2 tidal constituent; 2D numerical model Hydrodymics is fundame

19、ntal in affecting flow and sediment transport in estuaries and rivers as well. Simulation of the flow in the Bo Sea has been done with a 2D depth-integrated numerical model，to have a more in-depth understanding of the interactions between water level, speed, and flow direction there. With the modeli

20、ng, we could deduce its implications to the physical modeling for the Yellow River estuary, compiling of nearshore bathymetry, and feasibility of the scenarios for the various training strategies. 1 Governing EquationsThe 2D depth integrated numerical hydrodynamic model we used is a numerical modeli

21、ng tool for modeling of the flow field in estuaries and open seas due to river feedings, wind, tides, waves, salinity, temperature and other factors. The fundamental equations are given as follows.Continuity equation: Momentum equation in X direction: Momentum equation in Y direction:Where is flow l

22、evel (m)，d bed elevation(m)，h flow depth（d）（m），p and q flow rates per unit width in X and directions respectively(m3/s/m)，where p=hu, q=vh, u and v velocity components in X and Y directions respectively, C(x，y) Chezy Coefficient(m1/2/s)，bearing the relationship with Manning Coefficient as c=n-1h1/6,

23、 g gravitational acceleration(m/s2)，f(v) wind friction coefficient，V、Vx、and Vy wind speed and its components in X and Y directions respectively (m/s)， Corilis Force factor s-1， pa air pressure(kg/m/s2), w water density (kg/m3), x and y distance in X and Y direction(m),t time(s)，E eddy viscosity，whic

24、h is calculated with Smagorinsky Formula:In order to investigate the basic flow features in the Bo Sea, weinclude only tides in the modeling without considering the other factors for the time being.2 Conditions specification and model parameters calibrationThe simulated area covers the whole Bo Sea

25、west of the Bo Sea Strait, excluding the rivers flowing in. The bathymetrical data for the initial bed topography is composed of two parts: one is the data set of the measurements in the nearshore area bordering the Yellow River Delta, which was made in 1980 and 1981. The other part the sea chart of

26、 a scale of 1:1,000,000. The modeling domain is cut into squares, each side 5000 m long.Water level hydrographs at the Bo Sea Strait are taken as boundary conditions.Simulation period starts from July 1,1981 to Aug.22 1981 with the time step 30s.The parameters to be calibrated are (Cs) in the Smagor

27、insky formula and Manning Coefficient (n). When Manning Coefficient is 0.0125s/m1/3 and Smagorinsky Coefficient is 0.5, the model results are found in a nice agreement with the measured ones in terms of the tidal constituent (M2) (Fig.1). Details of the calibration and verification are omitted for c

28、oncise purpose.（a） calculated tidal ranges of constituent M2(b) Measured tidal ranges of the constituent M2 Fig.1 (a)(b) calculated and measured tidal ranges of the constituent M23Sensitivity analysis of model parametersSensitivity analysis is done to evaluate how the parameters (Cs and n) affect th

29、e flow in the Bo Sea. It is found that the increase in Cs, given a constant Manning coefficient, leads to the decrease in flow level and velocity with, however, a small margin (Fig.2). In contrast, it is found that for a constant value of Cs given, amplitudes of tidal level and velocity are lowered

30、by a larger margin each with Manning coefficient increasing (Fig.3)Fig.2 Water level hydrographs for various Cs values and n of 0.0125 s/m3Fig.3 flow level and velocity hydrographs for various n and Cs of 0.54Flow features4.1 Tidal current going from low to high level and from high to low level Rive

31、r flow usually goes from high to low elevation while water in the tides-affected estuary is found able to run from high to low at one moment and run from low to high at another moment, which is described in a little more details as follows. Fig.4 shows the level, speed and direction of the flow at a

32、 site in the center of the Bo Sea Bay, where the flow direction is measured clockwise with respect to the north. Fig.5 shows flow level marked on the contour lines and velocity vectors.Fig.4 Flow level, speed and direction near the center of Bo Sea Bay.At nine oclock the site in Bo Sea Bay is at a h

33、igh water (Fig.4) while water level at the Bo Sea Strait is much lower (Fig.5), consequently forming a surface gradient in the eastern direction, much the same as the flow direction. Meanwhile, flow in the Laizhou Bay goes from low to high level.Later on, water flowing out of the Bo Sea Bay graduall

34、y results in the flow level decreasing and current speed increasing while the water level at the Bo Sea Strait is seen rising and the current speed decreasing. Around 12:20, the Bo Sea Bay is at a low water, and the water is continuously flowing east, which means the flow, which goes from low to hig

35、h level, is approaching the maximum. Meanwhile, the water at the Strait is approaching the high water with flow speed close to zero. Later on, the flow level in the Bo Sea Bay is continuously going lower and lower with flow speed turning into a decreasing process, which is going on until 14:20 when

36、the east-going flow in the Bo Sea Bay is approaching to zero and the water at the strait flows westwards from the high level at the strait to the low level in the bay. After 15:20 the flow in the Bo Sea Bay turns into a rising process in water level, going westwards from high to low level.In a word,

37、 the flow in the Bo Sea have the moments when the water flows from high to low and the moments when flows goes from low to high. The processes are the ones when kinetic energy and potential energy are transferred between each other. It means that the water surface slope/gradient is a force to accele

38、rate the flow in a certain period while it is a force to slow down the flow for another period. This is one of the major differences from the unidirectional river. The feature poses more challenging requirements for the techniques and equipments for flow control for the physical model for the Yellow

39、 River Estuary.(a) (b)(c)(d）Fig.5（a）, (b), (c) and (d) Calculated current velocity and water level4.2 Rotary flow featuring the flow in Bo SeaWhen Tidal waves come into the Bo Sea, it is characterized by maximum velocity occurring usually at halfway from high waters and low waters, and the high wate

40、rs occurring usually at current-slack moment (Fig.4). The features indicate the nature of rotary waves in the Bo Sea.The tidal waves rotate clockwise in the Bo Sea bay while it does counterclockwise in the Laizhou Bay. Closer to the shore, the tidal current demonstrates more of to-and fro current th

41、an the rotary current. 5. Applications of numerical modeling 5.1 to define nutshells in controlling the flow in the physical model for the estuary The flow features as mentioned above are, in fact, the nutshells in control the flow in the physical model for the estuary.5.2 to prediction water level

42、for nearshore bathymetric data compiling The elevation (Hb) of any site in the nearshore area bordering the Yellow River Delta is obtained by subtracting the value of the water depth (h) from the value of the water level (Z) (i.e., Hb=Z-h) while Z is assumed the same as that measured at the shore (F

43、ig.6). The easy principle, unfortunately, doesnt lead to satisfactory results, which often shows heavy sedimentation in the deep water, opposite to the truth.Sea bed , Hb Water level Z ele. DatumhFig.6 Relationship of elevation of seabed, water depth, and water levelLets find out the reasons for the

44、 false sedimentation. With the modeling results, we have found the differences between the values of water levels at the two ends of the cross shore profile No.12 (Fig.7) varying between -0.4 m-+0.6 m. In the same way, we have found the water level differences at the two ends of all the 35 profiles

45、(Fig.7) varying between -0.4 m and +0.5 m (Fig.8).In a word, it is necessary to use the hydrodynamic model to give the water levels at any sites when compiling the bathymetrical data for accurate results. Fig.7 35 cross shore profiles in the nearshore around the Yellow River DeltaFig.8 water level d

46、ifferences at the both ends of CS 12Fig.9 Maximum differences of water levels at the two ends of the each cross shore profile5.3: to investigate feasibility of proposals for the estuary trainingHow to train the Yellow River estuarine course has been a hot issue in China so far. Experts made various

47、proposals, which could be cost-efficiently evaluated by numerical modeling. Study on Proposal 1One of the proposals for the training is by extending the course seawards as long as long as possible as to make the best of the velocity increased at the river mouth to carry as much sediment load to the

48、deeper waters offshore as possible. Lets check whether or not it is justifiable by modeling the scenarios where the river mouth extends eastwards and northwards, respectively.Fig.10 shows the differences between the maximum velocity values at the river mouth when it extending eastwards 30 km and 20

49、km. It is found that the river mouth extending can increase the maximum velocity near the river mouth as expected, which, unfortunately, covers quite a small area, unexpectedly. Similar results have been found when the river mouth extends northwards. The results above mean sediment delivered out of

50、the river mouth will not be transported very far, therefore, deposition of coarse sediment fraction can be expected to remain around the river mouth.In addition, considering previous research results to the effect that long extension of the river mouth leads to the rising of the river bed upstream o

51、f the river mouth although there are great differences between the understandings of the effect in terms of the length affected, we, therefore, are confident to make the conclusion that the proposal for training the Yellow River estuarine course by extending seawards far enough so as to deliver more

52、 the sediment load into deep water is not feasible. Fig.10 Differences in maximum current speed values when river mouth extends (note: positive values in the figures means the increasing of maximum velocity values while negative values the decreasing). Study on proposal 2Another proposal for the tra

53、ining by other experts is by increasing the flow rate released out of the river mouth so large as to carry more sediment with it into the farther water. Lets resort to the model to simulate the flow field when the flow out of the river mouth maintains as large as 5000 m3/s with various injection ang

54、les ranging from 0, 45, 90, 135, 180 degrees from the north clockwise, respectively. It is found from the modeling results (Fig.11) that the injected flow rate, can increase the local velocity as expected, around the river mouth, which, however, covers a small area, 40 km×30km, unexpectedly. Th

55、e result isnt a good supporter for the training proposal. Fig.11 Area affected by the Flow released out of the Yellow River Mouth 5.3.3 A practical way out: through vast space for the estuarine course to run for long term Lets consider the facts that the Yellow River mouth always extends seawar

56、ds longer when a larger sediment load coming to the estuary, and that the shore near the river mouth is found under significant erosion when less sediment load coming, say, since the Dikouhe River was abandoned in 1976. With these facts in mind, we could deduce that sediment load is one of the signi

57、ficant factors affecting the estuarine morphology. So how to control the incoming flow and sediment load will be key factor in making strategies for training the estuary.Ideally, a certain magnitude of sediment coming to the estuary, not too much or not too less, can prevent the estuarine channel from withering in channel geometry and flood conveyance, and can lead to less effect of the deltaic shore extension. To realize the proposal, however, will demand one large hydraulic complex, or even more, to regulate water and sediment