冯立文,姚清河, 苏炜, 王生, 徐文兵
(中山大学航空航天学院,广东 广州 510275)
自20世纪80年代以来,广东珠江口附近进行着广泛的围海造陆工程。这虽缓解了经济发展和用地不足的问题,却在一定程度上改变了潮流的流速、流向和水文条件,人为加剧了航道的淤塞情况。可见,河口、海岸的泥沙输移问题是河口海岸工程研究的重要内容,泥沙输移问题与港口的选址规划布置、海岸线的淤积侵蚀、海洋环境保护都密切相关。因数学模型试验不存在比例尺的限制,其成为了解决海岸、河口工程泥沙问题所广泛采用的方法。如:辛文杰[1]进行了潮流波浪综合作用下河口二维悬沙数学模型的研究;窦国仁[2]对河口海岸全沙模型相似理论进行了研究;朱志夏等[3]进行了海岸悬沙运移数学模型的研究。以上成果为采用数学模型进行潮流作用下泥沙运动的研究提供了重要理论基础和经验。然而,泥沙输移的影响因素较多,如波浪、水流、紊流强度、水下地形、泥沙类型及泥沙的沉积等,且现场观测困难。关于泥沙运动的基本问题还没解决,亦无成熟的输沙公式。泥沙输移的数学模型虽在一定程度上可作为实际工程的依据,却仍有许多问题需要深入研究。
洪湾渔港工程地处珠江河口,位于磨刀门水道和洪湾水道交汇处的左岸,靠近入海口,所以其既受珠江流域径流的影响,又受珠江磨刀门外海的潮汐作用,水域水动力特性复杂。受来自上游磨刀门水道排沙、以及西南沿岸洋流所携带的伶仃洋浅滩水沙的影响,水域河道及堤岸不断变化;受径流携带的水沙、外海潮汐动力的影响,泥沙冲淤变化速度较快[4-5]。本文根据收集到的实测资料、历史地形资料,建立了一维、二维潮流泥沙数学模型,并对模型进行了潮位、潮流、含沙量等方面的验证;计算了工程前后其附近水域的潮流、泥沙运动特征,预测分析了工程建设后的泥沙回淤情况,并根据计算结果提出了优化措施。
洪湾渔港位于珠海市洪湾建材码头西北侧、挂定角山东南侧,洪湾水道与磨刀门水道交汇处,具体位置如图1所示。规划陆域面积为0.6 km2,海域面积为0.4 km2;码头岸线和渔用岸线长度均超过1 000 m;规划可停泊渔船1 000艘以上,年卸港量8万t以上。规划方案中,洪湾渔港口门位于港池西北方,向西面对磨刀门水道。口门宽度约150 m。港池内水域北宽南窄,内航道宽度为60 m。口门处航道略宽于港池内,达90 m,通过支航道与磨刀门航道相接,支航道长度约350 m。规划方案见图2。
图1 工程地理位置图Fig.1 Geographic location of this project
图2 规划方案总平面布置图Fig.2 The general layout of the construction plan
为了分析工程所在河段的水沙特性以及预测工程建设对水沙运动的影响,本文采用二维潮流泥沙模型的方法。选取特征水文组合作为计算条件,用工程河段的二维潮流泥沙模型,计算详细的水动力特性及泥沙冲淤变化情况。
一般而言,由于水流对冲流动所引发的泥沙输运和迁移,要远远大于泥沙的自然扩散,所以我们采用简化的二维平面悬沙运动方程[6]:
(1)
式中:h为水深,C为水体含沙量,t为时间,u和v分别为x和y方向的速度,k为泥沙恢复饱和系数,W为泥沙沉降速度,C0为水流挟沙能力。而,
(2)
式中:α为经验系数,取0.023。β为经验系数,取0.000 4。γs为泥沙容重,γ为水的容重,Cz为谢才系数,T为波浪周期,hw为波浪高。
水流连续方程:
(3)
式中:ζ为水面高程;p,q为x、y方向的单宽流量m3/(s·m-1) ;S为源汇项[7-8]。
水流运动方程:
(4)
式中:g为重力加速度,f为风摩擦系数。V、Vx、Vy分别为风速及其在x、y方向上的分量;Ω为柯氏力参数;Pa为大气压强;ρ为水的密度;S、Six、Siy分别为源汇项及其在x、y方向上的分量;τxx、τxy、τyy为有效剪切力分量[9-10]。
泥沙连续方程:
(5)
其中,n为河床孔隙率,z为河床高程,Sx为全沙或悬移质沿x方向输移量,Sy为全沙或悬移质沿y方向输移量,△S为河床冲淤量[11-12]。模型采用不平衡输沙模式,这种模式下对流扩散方程可以写成如下形式:
ΔS=Φ0(η0)ωs(c-cε)
(6)
其中,η0是流速为零的床面,Φ0为单位河底高程线函数,ωs为悬移质泥沙沉速,c为床沙沿深度方向平均值,cε为垂向平均的饱和含沙量。该式表示,如果水体中的含沙量大于饱和含沙量则发生淤积,如果水体中的含沙量小于饱和含沙量则发生冲刷[13-14]。
水动力边界条件分为开边界与闭边界两种[15]。开边界是指有水流出入的边界,如分洪口门等;闭边界是指不过水的边界,如堤防。处理洪水问题通常要引入干湿边界。当水深小时,把单元格动量置为零,仅仅计算质量。当水深小到可以忽略时,则该单元格不考虑在计算中[16-17]。水深在每个单元格的值是可以监测的,根据水深的大小,单元格被分为干、半干和湿单元格。如果单元格水深比特定水深hdry小,并且单元格没有湿边界,则这个单元格不再计算;如果单元格水深比hdry小,另一个面是湿边界,认为动量为零仅仅计算质量;如果单元格水深比hwet大,是湿单元格,则动量和质量都将被计算。hwet和hdry必须满足hdry 默认的值是hdry=0.005 m,hwet=0.1 m。如果把hwet值设置的很小,则可能出现偏离实际的大流速值从而引起不稳定的问题[18-19]。 二维潮流模型的长约60 km、宽约25 km,从上游的磨刀门水道的控制站-竹银至下游口门, 如图3所示。模型范围既包含了拟建工程的影响范围,又可以充分考虑河道水流的特征。研究采用三角形网格划分程序对二维模型计算区域剖分,并对港址所在水域局部网格进行了加密,共布置网格26 315个。最小网格尺寸为10 m,如图4所示。 研究水域内有江心洲和边滩,为了正确模拟这些江心洲和边滩在涨、落潮期间淹没及出露的不同状况,模型采用动边界技术对其进行处理。即,将落潮期间出露的区域转化为滩地,同时形成新边界;反之,将涨潮期间淹没的滩地转化成计算水域[20-22]。 图3 二维网格模型Fig.3 Two-dimensional mesh 图4 工程附近局部二维网格Fig.4 Local two-dimensional mesh near the project 2.2.1 潮位 从图5可知,在该水文条件下,模型与原型的潮位过程线吻合良好,模型的涨、落潮历时和相位与原型实测资料一致,潮位特征值验证误差都小于±0.10 m,符合《海岸与河口潮流泥沙模拟技术规程》规定的精度要求。 2.2.2 含沙量 各测站含沙过程见图6。从图中可以看出,模型计算所得的各垂线悬沙浓度变化过程与观测资料的趋势一致。 图5 灯笼左的水位实测值与计算值对比Fig.5 Comparison of measured and calculated water levels on the left side of Denglongzuo 图6 灯笼左含沙量实测值与计算值对比Fig.6 Comparison of measured and calculated sediment concentration at the left of Denglongzuo 以上结果表明:模型能够较好的模拟工程水域的潮流运动、悬沙分布等情况,其精度满足验证要求,说明模型概化和计算参数选择合理,可以用于工程方案的潮流泥沙计算分析。 根据已有资料及二维潮流模型的计算结果,可对工程前后的流场变化情况进行分析。流场变化分析中涉及工程附近水域整体流态、港池和航道等位置的流速特征,及工程特征点位置的潮流变化情况等等。为此,需布置若干采样点,采样点位置如图7所示。模型中共布设有52个流速采样点,1#-12#采样点布置在港池内,13#-15#采样点在口门中心,16#-27#采样点位于进港支航道上,28#-34#采样点位于磨刀门航道上,35#-43#与44#-52#采样点分别位于工程上、下游。 图7 流速采样点位置图Fig.7 Sample point locations of velocity “99·7”洪水水文条件下,工程前后大潮涨急时刻流场及流速变化,如图8所示。工程口门附近地形为浅滩,水深较浅,加上上游挂锭角水闸凸出岸线的挑流作用,口门及航道受该挑流作用流速不大。港池内,由于水域纳潮量不大,进出口门流量小,其流速更小。 图8 方案实施前后,洪季大潮涨急时刻的流场和流速Fig.8 The flow field (a) and velocity (b) when the tide rises sharply in flood season 洪季大潮涨急时刻,工程建设前各特征点的最大流速为0.34 m/s,出现在靠近主流的磨刀门航道上。工程建设后,新增加的港池水域内流速一般不超过0.02 m/s。进港支航道内流速均有所减少,最大减少值为0.08 m/s,出现在19#点。工程下游流速有所增大,最大增幅为0.03 m/s(口门正南侧靠岸位置,46#采样点)。流向的变化主要发生在进港支航道位置,由于进港支航道走向与水流方向基本垂直,随着航道的浚深,水深在小范围内出现较大幅度的变化,对水流有折射作用。 “99·7”洪水水文条件下,工程前后大潮落急时刻流场及流速变化,如图9所示。相比于涨急时刻,洪季大潮落急时刻,渔港口门及航道处水动力有所增强;工程建设前各特征点的最大流速为0.92 m/s,出现在靠近主流的磨刀门航道上;港池内流速量值变化不大,流速一般不超过0.03 m/s。进港支航道内流速均有所减少,最大减少值为0.22 m/s,出现在22#点。工程下游流速有所增大,最大增幅为0.06 m/s(口门正南侧浅滩,45#、46#采样点)。 流向的变化主要发生在进港支航道位置,主要表现在靠近主流一侧的采样点流向逆时针偏转,而靠近口门一侧的采样点流速顺时针偏转,后者与主流和港池经口门向外的落潮流交汇有关。 图10为“01·2”枯水水文条件下,工程前后涨急时刻流场及流速变化情况。 图9 方案实施前后,洪季大潮落急时刻的流场和流速Fig.9 The flow field (a) and velocity (b) when the tide ebbs in flood season 图10 方案实施前后,枯季大潮涨急时刻的流场和流速Fig.10 The flow field (a) and velocity (b) at the high tide in dry season 由计算结果可知,在枯季大潮涨急时刻,工程建设前各特征点的最大流速为0.59 m/s,出现在靠近航道下游的47#点。工程建设后,新增加的港池水域内流速一般不超过0.03 m/s。特征点流速值的变化规律与流速等值线变化趋势基本相同,进港支航道内流速均有所减少,平均减少值约0.1 m/s;最大减少值为0.18 m/s,出现在靠近口门的航道主槽(20#采样点)。工程上下游流速有所增大,最大增幅为0.06 m/s,出现在航道上游,口门西北方的浅滩(42#采样点)。水流流向变化与洪季大潮趋势一致。 图11为“01·2”枯水水文条件下,工程前后落急时刻流场及流速变化情况。 图11 方案实施前后,枯季大潮落急时刻的流场和流速Fig.11 The flow field (a) and velocity (b) of spring tide in the dry season 从计算结果可知,在枯季大潮落急时刻,工程建设前各特征点的最大流速为0.66 m/s,出现在靠近主流的磨刀门航道上。工程建设后,新增加的港池水域内流速一般不超过0.04 m/s。特征点流速值的变化规律与流速等值线变化趋势基本相同,进港支航道靠近主流的位置流速均有所减少,最大减少值为0.20 m/s,出现在中部航道的主槽内(22#采样点)。流向的变化主要发生在进港支航道位置,其流向多为顺时针偏转。 工程附近的冲淤平面分布,如图12所示。从图中可以看出,在正常年份的来水来沙条件下,工程方案实施后,浅滩仍然呈现淤积态势,淤积呈带状分布,平均淤积厚度约0.30 m。拟建工程的进港支航道靠近口门位置出现明显的淤积,淤积的分布比较集中,最大淤积厚度达0.92 m;口门内和港池淤积则逐渐减小。 图12 规划方案实施后,年淤积分布图Fig.12 Years siltation distribution after the implementation 本文主要研究洪湾渔港工程建设对周围水域潮流运动和泥沙淤积的影响,经计算可以得出以下结论: 1)工程对于水流流速的影响主要在港池、口门及航道附近,流速相比工程前有所减小,涨落急时刻减小值约0.05-0.20 m/s。 2)工程对于水流流向的影响主要发生在进港支航道位置,对主流的影响较小。 3)规划方案实施后,港池内年平均淤积厚度为0.21 m,淤积量为8.5万m3;进港支航道主槽部分淤积厚度较大,年平均淤积0.53 m,淤积量2.0万m3;进港支航道边坡的年平均淤积厚度为0.26 m,淤积量为0.8万m3;工程开挖水域内一年总淤积量为11.3万m3。最大淤积强度约0.9 m,出现在口门外侧。 4)工程实施后进港支航道的平均回淤强度较大,最大淤积部位在口门附近,约为0.92 m。淤积的原因一方面是因为口门靠近岸边,且上下游有不断发育的浅滩,需要通过疏浚才能达到设计水深。另一方面是因为口门附近流速缓慢,动力很弱,这些因素会导致口门附近出现较明显的淤积。 5)针对口门淤积的问题,可通过以下措施进行优化:a.出口大范围开挖。按照设计深度,在约2倍口门宽度范围开挖,按照1∶10斜坡放坡,向两侧均等延伸;b. 两侧开挖区,具有预留沉砂池的功能,避免上游淤泥直接进入航道;c.砌筑短丁坝。口门两侧上下各20 m处,设20 m长的丁坝,浅滩淤泥可能在风浪条件下汇入航道,短丁坝主要是拦截近岸淤泥。2 模型建立与验证
2.1 计算域及模型的动边界条件
2.2 模拟结果
3 结果分析
3.1 工程前后洪季大潮涨急时刻流场分析
3.2 工程前后洪季大潮落急时刻流场分析
3.3 工程前后枯季大潮涨急时刻流场分析
3.4 工程前后枯季大潮落急时刻流场分析
3.5 工程水域泥沙回淤计算与分析
4 结 论