郝林南，杨小宸

（1.辽宁省防汛抗旱指挥部办公室，辽宁沈阳110003；2.沈阳农业大学，辽宁沈阳110044；3.辽宁省水利水电科学研究院，辽宁沈阳110003）

辽河口三维水动力数值模拟

郝林南1,2，杨小宸3

（1.辽宁省防汛抗旱指挥部办公室，辽宁沈阳110003；2.沈阳农业大学，辽宁沈阳110044；3.辽宁省水利水电科学研究院，辽宁沈阳110003）

文中基于有限体积海岸海洋模型FVCOM，建立了辽河口三维水动力数值模型，利用潮位、流速实测资料对模型进行验证并分析辽河口三维水动力运动特性，可为进一步研究辽河口泥沙运动及盐度、污染物扩散等提供基础。

模型验证；辽河口；三维水动力；数值模型

辽河口地区位于辽宁省的南部，是我国海岸线的最北端，面临渤海；区域内自东向西有大辽河口、辽河口（双台子河口）、大凌河口三个独立入海口门。辽河流域属温带季风气候，年径流量为89亿m3；年降水量约为350～1 000 mm，并呈自东南向西北递减的特点。流域夏季多暴雨，强度大、频率高、集流快，常使水位陡涨猛落，造成下游地区洪涝。此外，辽河的含沙量较高，仅次于黄河、海河，为中国第三位，年输沙量达2 098万t。辽河口地区为平原型淤泥质海岸，有着丰富的海陆相来沙，近百年来发展的结果已成为拦门沙发达的堆积型河口。不断淤积的河口使中水河槽过水面积相对减少，使纳潮量减少和泄洪能力降低。此外，辽河口地区还存在淡水资源严重不足，水环境污染严重。

辽河口近海水体中的泥沙运动及水环境污染等问题是水动力作用过程的直接表现，对辽河口水动力特点的合理描述是研究上述问题的基础。因此，下文建立辽河口三维水运动数学模型，分析河口区域水流运动规律，可为整治工程提供指导依据，可缓解河口地区淤积情况，同时可为河口附近的港口选址、油气开发平台、航道疏通工程和海滩旅游等建设带来的相关环境问题的分析提供技术支撑，从长远来看，对于辽河口的资源开发利用、河口治理、环境变化和经济发展都具有重要的现实指导意义。

1 模型建立

基于有限体积海岸海洋模型FVCOM建立辽河口三维水动力数值模型。FVCOM模型是由美国马萨诸塞大学海洋科技研究院与伍兹霍尔海洋研究所共同开发的基于无结构网格、有限体积法、三维原始方程的海洋模型。模型包括动量方程、连续性方程、温盐输运方程和状态方程，通过二阶紊流闭合模型封闭方程。水平方向使用非结构化三角形网格，能更好地描述近岸地区的复杂岸线和建筑物，垂向采用坐标变换。数值方法采用有限体积法，能够很好地将有限差分法离散结构简单、计算效率高的优点和有限元方法处理复杂边界的几何灵活性结合起来，可以更好地保证具有复杂地形和岸线的近岸、河口区域的质量、动量、盐度以及热量的守恒。

1.1 控制方程

FVCOM三维水动力模型基于自由表面的三维原始控制方程，由动量方程、连续性方程和状态方程构成，基于笛卡尔坐标系的三维原始方程可表示为：

其中：x，y和z为笛卡尔坐标系下的三维坐标，分别为东西方向、南北方向以及垂向的坐标，m；u，v和w分别为x，y和z轴的速度分量，m/s；g为重力加速度，m/s2；t为时间，s；ρ为密度，kg/m3；p为压力，Pa；f为科氏力参数；Km为垂向涡粘系数，m2/s；Fu和Fv为水平动量扩散项，m/s；T为温度，℃；S为盐度，‰。

垂向采用坐标变换：

其中，ζ为自由表面，m；H为海床相对于基准面的距离，m。因此，D=H+ζ为总水深，σ坐标的变化范围为[-1，0]。

1.2 数值离散

数值离散采用有限体积法，将计算区域划分为若干互不重叠且相互邻接的三角形网格单元。模型采用二维外模和三维内模交替计算的方法，由外模态计算垂向平均流速和水位，由内模态计算三维流速、泥沙运动和温盐输运等物理量。二维外模将动量方程和连续性方程在控制体内积分，通过改进的四阶Runge-Kutta时间积分求解。三维内模采用隐式与显式结合的方法，对流项采用显式求解，垂向扩散项采用隐式求解。具体的数值离散过程见文献。

2 辽河口三维水动力数值模拟

2.1 模型设置

模型采用非结构化的三角形网格，通过大、小模型嵌套计算，由大模型为小模型提供边界潮位条件。大模型计算范围为122.5°E以西海域，包括整个渤海及渤海海峡。大模型区域共包含12 856个节点，划分为24 969个三角形单元，最小网格步长约为1 800 m。垂向平均分为5层，大模型开边界水位由中国海域潮汐预报软件Chinatide提供，外模计算时间步长为1.0 s。小模型计算范围包括西至葫芦岛、东至鲅鱼圈连线北侧海域，同时小模型北端向辽河口上游延伸，其中双台子河上游延伸至盘锦市盘锦县的盘山闸，大辽河延伸至鞍山海城市的三岔河，见图1。小模型区域共包含18 118个节点，划分为34 152个三角形单元，在河口及河道内加密，最小网格步长约为30 m。垂向平均分为10层，河口地形取自1∶500 00海图，河道地形资料采用2011年大断面资料。小模型开边界潮位由大模型计算结果提供，河流流量边界根据水文站实测资料确定，双台子河近似采用六间房测站数据，大辽河近似采用邢家窝棚和唐马寨测站之和，大凌河根据凌海测站资料显示流量相对较小，暂不考虑。模型外模时间步长取为0.1 s，内模时间步长为1.0 s。

图1 小模型计算范围及地形图

2.2 模型验证

收集到2013年9月20日10:00至 2013年9月21日16:00大潮期间在计算海域的水文全潮观测结果，用于模型验证。测站位置见图1，其中包括潮位测站T1，T2及分层流速、流向测站C1～C4。计算时间选为2013年9月18日10:00至2013年9月21日16：00，给水动力计算充分的稳定时间。通过大、小模型嵌套计算，由大模型为小模型提供开边界潮位条件，河流流量根据水文站数据确定为双台子河191 m3/s，大辽河111 m3/s。计算结果与实测数据吻合较好，能够较好地反映辽河口海域潮流变化情况。

计算结果表明：落潮流速略大于涨潮流速，表层和底层流向基本一致，表层流速大于底层流速。图2给出了涨急时刻和落急时刻计算区域表层的流场图，由图2中可以看出，辽河口海域潮流基本呈往复流运动。涨潮期间，外海潮流自西南向东北方向近岸传播，盖州滩东侧深槽（即测站附近）为主流区，流速大、流势强，水流主要通过这股涨潮主流分别流向双台子河口和大辽河口。在双台子河口外约12 km处，河流流量和外海潮流形成顶托，并向西北延伸至近岸区域，该现象可能会导致泥沙在此处滞留并沉积。落潮期间，流向与河流流向一致，双台子河口水流在盖州滩处分为两股流向外海，盖州滩东侧深槽依然是落潮流的主流区。

3 结 语

基于FVCOM模型建立的辽河口三维水动力数值模型，采用非结构化网格建立大、小模型并模拟了大潮期间辽河口海域的水动力运动，研究结论主要有：

1）通过与实测资料对比，模型计算结果与实测数据吻合较好，此文所建模型能够合理描述辽河口三维水动力特性，可为进一步研究辽河口泥沙运动及盐度、污染物扩散等提供基础。

2）辽河口海域潮流运动呈往复流，盖州滩东侧深槽为涨落潮的主流区。落潮流速略大于涨潮流速，表层流速大于底层流速。涨潮期间，双台子河口外潮水与河流的顶托作用明显，可能对泥沙产生滞留沉积作用。

［1］Chen C，Liu H，Beardsley R C.An unstructured,finitevolume,threedimensional，primitive equation ocean model：application to coastal ocean and estuaries［J］.Journal of Atmospheric and Oceanic Technology,2003，20（1）：159-186.

TV131

A

1002－0624（2016）12－0044－02