纯隐格式的混合有限分析法在广雅桥水流数值模拟中的应用

方神光,黄胜伟

(珠江水利委员会珠江水利科学研究院,广东广州 510611)

在河道涉水工程中,新建桥墩阻水最为常见。由于占用了河道部分过水面积,会引起河道局部壅水、近岸流速增大等,因此桥墩壅水及其对河道防洪的影响受到广泛关注。有关桥墩壅水的理论和应用研究成果较多[1-6]。目前对桩墩的模拟主要采用2种模式[6]:一是采用局部阻力修正;二是直接进行模拟。对于大范围的计算区域,为提高计算效率,常采用较粗的网格,而桥墩尺寸太小,因此,常采用局部网格节点阻力修正的方法来体现桥墩对水位壅高的影响,较常用的局部阻力修正方法如赵晓冬[5]提出的桩群阻力修正公式。当要了解桥墩前后的局部流态以及河势变化时,则需要采用精细网格直接进行模拟的方法,即将桥墩所占据的区域作为陆域进行模拟计算。

在桥墩局部水动力数值模拟计算中,采用纯隐格式的混合有限分析法[7]来离散和求解常用的曲线坐标系下的二维水动力数学模型。该数值方法常用于局部区域流场的精细模拟,保留了有限分析法的优点,同时避免了有限分析法中无穷级数带来的不便,其计算结果准确、有效[8-9]。笔者以柳州市部分柳江河道和该河道上拟建的广雅桥为例,分别建立大范围柳江河道的二维水动力数学模型和广雅桥局部二维水动力数学模型,研究探讨该桥梁工程建成后对局部河道壅水和水流流态的影响,并采用Fortran语言编制了计算软件平台。其中,大范围的二维水动力数学模型主要用于为桥墩局部水动力数学模型提供初始条件和边界条件。

1 数学模型

正交曲线坐标系下的二维水深平均水动力数学模型如下:

连续方程:

式中:ξ,η为曲线坐标系(无量纲);U,V分别为ξ,η方向上的垂线平均流速,m/s;ζ为水位(基准面到自由水面的距离),m;Q为源汇项,m3/s;H为总水深,m;,为曲线坐标到直角坐标的变换参数 ,m;f为柯氏系数,s-1;g为重力加速度,m/s2;νt为紊动黏性系数,m2/s;t为时间,s;C为谢才系数。

2 数值计算方法

二维潮流数学模型中的偏微分方程可以写成如下统一形式:

式中物理量符号含义见文献[7]。应用纯隐格式的混合有限分析法对该方程进行离散,更为详细的推导步骤参见相关文献[7]。此处采用C型网格结合SIMPLER算法进行速度和水位的耦合求解。大模型中的桥墩阻力概化采用南京水利科学研究院的桩群阻力试验研究成果[5]。

由于采用曲线规则网格,计算是沿行或列进行的,因此计算前首先需要对计算区域的复杂边界进行识别,其主要方法是将不参与计算的区域的地形数据设为一很大值(如10m或100m),通过编写边界识别程序,将每一计算行或列的水陆边界识别出来,并存入一边界数组中,计算时通过读取该边界数组中的边界数据来确定每一计算行或列参与计算的范围。计算区域存在浅滩,对浅滩的处理方式影响到计算的稳定、收敛以及计算结果的精度。此处将Leendertse[10]在推出ADI数值模式的同时给出的干湿网格判别准则与河道水流数值模拟中常用的冻结法结合起来,准确有效地处理了该计算区域的浅滩。其基本原理是:①某瞬时任一行(或列)计算完毕,发现某网格节点处总水深H小于设定的最小水深Hmin,则认为干出,将该处糙率设为极大值(如108),即将该干节点处的流速冻结,同时该点水位为该点底部高程加上最小水深;②在任一时层计算前,根据干网格节点周围网格节点水位与该干网格节点水位的比较来确定该干网格节点是否被淹没,若被淹没,则将其糙率恢复为正常值。

3 计算模型网格及成果验证

3.1 河道概况及网格布置

柳江流经柳州的河段全长74 km,相对较为平缓。由于经济和社会发展迅速,需要在河道上兴建新的桥梁,拟在柳州市区柳江河道上修建广雅桥。为研究分析该桥梁的建设对河道局部壅水和水流流态的影响,采用大小模型嵌套的方式,大模型选取柳州市凤山镇—西流村全长63.9km的河段作为计算河段,如图1所示。小模型计算区域长0.686km,广雅桥东岸接广雅路,通往市中心;西岸接河西路和磨滩路,与西环线相连,如图2所示。

图1 大模型计算河段示意图

图2 广雅桥局部计算区域

大模型计算结果一方面用作验证,另一方面为小模型提供初始和边界条件。柳州市区从维义村至西流村全长20.01km的河段地形采用最新实测资料,对不完整的河道地形基础资料进行插补。大模型沿河流方向布置741个网格点,节点平均间距约为73m;垂直于河道水流方向布置21个网格节点,节点平均间距约为24m。小模型沿河流方向布置209个网格节点,节点平均间距约为3.3m;垂直于河道水流方向布置321个网格节点,节点平均间距约为2m。

3.2 小模型的边界选取方法

如前所述,本研究采用大小模型嵌套的方式来模拟桥墩的壅水及其对桥墩附近河道水流流态的影响。大模型采用桥墩阻力概化方式[5],用来进行大范围的桥墩壅水计算;小模型将桥墩直接概化成陆域,用来模拟分析桥墩对局部河道壅水及水流流态的影响。因此,小模型中将桥墩概化成不过水的陆域更符合工程实际情况。由于模拟计算方法和技术手段等的限制,大模型中对桥墩采用阻力概化,用来模拟桥墩引起的壅水变化,但桥墩上下游一定范围内的流态处于失真状态,在确定小模型的上下游边界时,若选取的边界处于这一范围内,将导致小模型计算得到的桥墩附近水流流态失真。

因此,为避免以上情况的出现,小模型上下游开边界位置的选取主要考虑两方面因素:一是上下游边界选取确保处于大模型中桥墩上下游流态失真范围之外;二是考虑到计算效率和成本,尽可能控制小模型的计算范围。选取步骤如下:①得到大模型的计算结果后,在大模型中大致确定桥墩上下游开边界位置,提取各边界断面上各节点的水深、平均流速和水位;②按照一定方法,将提取到的桥墩上下边界的网格节点上的流速和水位值赋给小模型边界网格节点,并确保由小模型上边界各节点计算得到的流量与大模型上边界设计流量一致;③小模型进行计算时,上边界各节点的水位和流速始终采用大模型中的提取值,下边界各节点的水位(或流速)值也采用提取值,下边界采用二类开边界条件;④小模型计算完后,将得到的下边界各节点上的流速(或水位)与大模型中提取的流速(或水位)进行比较,主要比较下边界各节点上流速(或水位)的变化趋势及各节点上数值的差,若两者变化趋势一致且各节点上数值的差很小,则认为选取的上下边界合适,若变化趋势不一致或各节点上数值的差别较大,则重新调整上下游边界位置,重复①～④的步骤,直至符合要求。

3.3 数学模型的验证

柳江河道属山区性河流,河床质基本上为砾石和卵石,根据以往研究成果,柳江河道糙率n变化范围在0.025～0.050之间。另外,数学模型中,n除反映河道粗糙度外,还包括了其他因素的影响,诸如数学模型自身的误差以及离散计算时的截断误差等因素,因此不同的数学模型和计算方法得到的n值会有所差异。数学模型中 νt的选取一定程度上影响到模拟计算的精度,该参数的选取同样受到所采用计算方法和计算网格尺寸的影响,可参考以往柳江河道工程实例中的取值范围进行选取,在此基础上,采用试算并与实测水位和流速验证的方法,确定大小数学模型的νt取值范围在0.1～10之间。

为对模型进行验证,上边界采用2009年7月5日柳州市柳江所遭遇20年一遇洪水流量(26700m3/s),下边界西流村水位根据以往的流量～水位关系插值得到。采用该模型进行数值模拟计算,统计水位计算值与实测值的差值见表1,两者差距在0.01m以内。模型计算得到河道主槽流速一般为2.5～3.5m/s,边滩流速一般在0.5～1.0m/s左右,如图3所示(左图为右图方框中的流场放大图),这与柳江河道实际流速相符,故该二维模型可用于柳江河道工程二维水流计算。

表1 2009年7月5日洪水计算值与实测值的差值

4 桥墩局部区域壅水及流态变化分析

广雅桥拟采用海鸥式双孔中承钢箱拱桥方案,桥梁全长966m。主桥采用海鸥式双孔中承钢箱,共有6跨,跨度分别为36m,60m,210m,210m,60m,36m,全长612m。下部结构:拱座采用C40钢筋混凝土,实心墩,承台尺寸为 12.40m×12.40m×6m,基础采用20根∅2.8m钻孔桩。广雅桥设计防洪标准为100年一遇,桥址平均阻水比为7.04%。

图3 2009年7月5日柳江局部河段洪水流场

从以上大模型河道计算结果中提取小模型中上下游边界条件和初始条件,以此计算工程前后广雅桥局部河道壅水和流态变化情况。图4给出了广雅桥上下游局部范围内水位和流速变化等值线,图5给出了工程前后流场。由图4可见,受桥墩阻水影响,桥墩上游一定范围内形成了壅水,墩后一定范围内形成了降水。桥墩河道断面平均壅水高度约为0.08m,紧邻墩前局部范围内壅水高度则更大。由流速和流态变化可见,广雅桥的建设对所处位置局部河道的主流向影响在3°以内。离岸和桥墩越近,变化值越大,岸边和桥墩局部区域流向变化可达到10°以上;桥墩主槽间的流速增大幅度在1m/s左右,岸边流速减小幅度在0.65m/s以内。以流速增大或减小0.5m/s为界,桥墩对下游流速的影响距离明显大于对上游的影响距离,桥上游的影响范围在50m以内。

图4 工程前后河道局部水位和流速变化等值线

图5 工程前后流场

表2给出了10年一遇洪水至100年一遇洪水下,计算和统计得到的广雅桥主墩和主槽的流速影响范围,左墩纵向和横向的影响距离分别在303m和121m以内,该范围内流速减小0.5m/s以上;左槽纵向和横向的影响距离分别在380m和124m以内,该范围内流速增大0.5m/s以上;中墩纵向和横向的影响距离分别在295m和61m以内,该范围流速减小0.5m/s以上;右槽纵向和横向的影响距离分别在400m和120m以内,该范围流速增大0.5m/s以上;右墩纵向和横向的影响距离分别在280m和70m以内,该范围内流速减小0.5m/s以上。由此可见,由于岸边和桥墩前后流速较大幅度的减小,广雅桥的建设会导致桥墩前后及两岸出现一定范围的淤积,主槽出现一定范围的冲刷,其影响范围主要集中在桥下游400m范围内。因此,广雅桥的建设对河道主槽的过流略有影响,可能引起局部河道地形和岸线出现一定程度的变化。

表2 广雅桥主墩和主槽的流速影响范围

5 结 论

a.将纯隐格式的混合有限分析法推广到工程实际应用中,采用Fortran语言编写了相应的计算求解软件平台,通过工程实例验证了该方法的适用性和有效性。

b.采用干湿网格判别准则与河道水流数值模拟中常用的冻结法结合起来,可以准确有效地处理计算区域中的浅滩。

c.计算结果表明,广雅桥建成后会对桥墩前后水位和流态形成一定的影响,桥墩上游一定范围内水位壅高约0.08m,墩间河槽流速明显增大,桥墩前后流速显著减小,影响范围主要集中在广雅桥下游400m范围内。

:

[1]陆浩,高冬光.桥梁水力学[M].北京:人民交通出版社,1991.

[2]SARKERM D.Flow measurement around scouredbridge piers usingacoustic-dopplervelocimeter(ADV)[J].Flow Measurement and Instrumentation,1998(9):217-227.

[3]张玮,解鸣晓.桩墩壅水数值计算方法研究[J].水利水电科技进展,2008,28(5):8-12.

[4]BALL D J.Simulation of piers in hydraulic models[J].Journal of the Waterways,Harbors and Coastal Engineering Division,ASCE,1974,100(1):18-27.

[5]赵晓冬.桩群阻力研究及模型码头桩群计算[R].南京:南京水利科学研究院,1996.

[6]解鸣晓,张玮,谢慧姣.桩群数值模拟中的概化方法研究[J].水动力学研究与进展,2008,23(4):464-471.

[7]李炜.黏性流体的混合有限分析解法[M].北京:科学出版社,2000.

[8]方神光.有障碍浮力射流特性的研究[D].武汉:武汉大学,2005.

[9]方神光,槐文信.静水中正方形孔口浮力射流特性的研究[J].应用基础与工程科学学报,2004,12(4):361-369.

[10]LEENDERTSE J J.A water quality simulation model for wellmixed estuaries and coastal seas:principle of computation[M].New York:the Rand Corporation,1970.