白玉川1， 温志超1， 徐海珏1，廖世智2， 曹永港2， 夏华永3

(1. 天津大学 水利工程仿真与安全国家重点试验室， 天津 300072； 2. 国家海洋局 南海调查技术中心， 广州 510300； 3. 国家海洋局 南海预报中心， 广州 510310)

摘要： 三维水动力模型在准确模拟海洋物理特性中起着重要的作用，传统的σ坐标转换由于当时计算机能力所限，舍去了复杂的高阶项，在实际复杂地形(或水深变化)环境下，会带来一定的误差或计算失真等问题．由此，为了适应高精度计算结果的需求，对原有σ坐标三维水动力模型进行了重新修正．在改进后的模型中，综合考虑了经σ坐标变换引入的与流速、水位、地形相关的复杂高阶项，选用特定的插值函数，利用有限元和差分相结合的方法，进行求解σ坐标下的完整三维浅水模式方程．相比原模型，改进的模型对底坡、水深、潮汐振幅等变化适用范围更为广泛，能更好地模拟出复杂水深变化下的垂向流动分布特征，计算结果具有更高的精确度；改进的模型针对一些极端水位条件(潮汐振幅与水深比大于0.15)，其计算误差同样可保持在一个较低的范围内；同时，改进的模型只需更短的时间就可运行至稳定状态．

关 键 词： σ坐标转换； 三维水动力模型； 有限元法； 差分法； 水动力特征

感潮河段、大型湖泊及近海区域的水流运动，实际上均为三维运动，但由于其水深尺度相对于流场水平尺度较小，采用二维潮流数学模型也能有效模拟其运动状态．随着河口、海域地形及水深变化复杂程度加大，以及实际工程问题对潮流垂直结构及河口分层运动特性研究需求的提高，高精度的三维水动力数学模型研究与建立也变得愈加重要．

在三维模型中，垂向坐标的选取是决定其对模拟水域适应性的关键因素之一，因此，垂向坐标及坐标变换的选取也就成为了人们一直关注的热点问题．自1957年Philips[1]提出σ坐标变换起，σ坐标因具有能较好地拟合床面和自由水面起伏变化的优势，进而由此引出的三维σ坐标变换水流模式已在各海洋、水利学科中得到了广泛应用．如窦振兴等[2]以具有自由面的三维非线性N-S方程为基本方程，经σ坐标变换得到模式方程建立数值模型用以计算渤海的三维潮流运动；白玉川等[3-4]选用特殊的插值函数利用有限元和差分相结合的方法，建立了海域数值模型并完成了歧口至套儿河海区及广西廉州湾水域的三维潮流计算；陈虹等[5]采用σ坐标变换，对渤海湾湾顶典型潮流与泥沙运动进行了数值模拟，以较高的分辨率揭示了该海域潮流与泥沙运动的典型特征；Bai(白玉川)等[6]应用σ坐标变换，得到该坐标下泥沙运动方程，建立了三维水流泥沙运动模型，分析得到了陆源泥沙和海向泥沙的运动规律及海河河口冲淤特征；汪守东等[7]在算子分裂法思想的基础上，将两种高精度的离散格式推广应用于三维对流扩散方程，并对ADI格式的对流向进行了改进，将对流项空间精度提高至4阶；郝嘉凌等[8]通过σ坐标三维水动力模型计算结果对潮流流速分布模式进行研究，得出了近底流速分布的线性分布模式；Warner等[9]将ROM模型与SWAN模型中波浪模块进行耦合，经σ坐标变换后得到了三维模式方程用以计算河口、大陆架及近岸海域的水流泥沙运动；Liu(刘昭伟)等[10]探索推导出了σ坐标下的废水排出后的输移扩散方程，用于计算自然河流单侧废水排出后的输移扩散特征；于守兵等[11]采用水位积分平衡法并引入边壁滑移系数，解决了σ坐标系下有限体积离散时产生的静水压力项与底坡项的平衡问题和边壁阻力模拟问题，建立了非结构网格上的三维浅水模型，并通过计算Holtz丁坝试验，表明考虑边壁阻力能够模拟丁坝附近局部小回流区；田勇[12]建立了基于有限体积法的琥珀三维水动力水质耦合数值模型，模型采用无结构三角形网格和垂向σ坐标来剖分计算域，以武汉大东湖水系为研究对象，对大东湖水系六个主要湖泊的水动力、水温和富营养化进行了全面研究；吴毓儒[13]利用FOCOM模型建立了长江口深水航道三维水动力、盐度数值模型，垂向采用σ坐标转换，对深水航道整治工程不同阶段的水流、盐度运动规律进行了模拟研究并预测了来自于拐弯段内两侧浅滩的泥沙进入航道是造成长江口深水航道严重淤积的主要原因；赵旭东等[14]应用非结构化网格建立水动力模型，对不同计算核心分配方案的MPI并行计算，基于GPU的高性能计算技术，在CUDA开发平台下设计并行算法，建立非结构化网格的二维水动力模型；李大鸣等[15]建立三维局部流场模型，模拟分析了发电装置对试验场局部近场、大范围远场水动力的影响情况；黄牧涛等[16]建立了σ坐标系下的三维水动力数值模型，对武汉东湖的三维流场结构和变化规律进行了深入分析．很多学者对三维水动力模型中垂向σ坐标转换进行了详细研究和探讨，转换后得到σ坐标系下的连续方程和动量方程，部分计算中保留了两项．随着计算机技术的发展，大多数计算中保留项增加了四项，而坐标转换中得到的其他高阶项一般被忽略不计．

本文在已有研究的基础上，针对传统σ坐标转换舍去的复杂高阶项，在计算复杂地形与水深变化环境下的水流运动时，对引起的一定误差或计算失真等问题展开研究，进一步完善了σ坐标三维水动力模型．在白玉川等[3-4,6]模型的基础上，综合考虑补充了经σ坐标变换引入的与流速、水位、地形相关的复杂高阶项，形成σ坐标下的完整三维浅水模式方程．在改进后的模型中，采用特殊的单元插值函数并在水平方向利用三角形网格剖分计算域，采用集中质量有限元方法离散方程和加速计算，保留了原模型计算快速的优点；在垂直方向上利用差分方法进行求解．

在一般宽、浅水域中，垂向加速度远小于重力加速度，可忽略不计，N-S方程可化简为如下形式：

连续方程

(1)

运动方程

(2)

(3)

(4)

式中x，y和z分别为直角坐标系下的方向坐标，U，V和W分别为x，y和z方向上的速度，P表示压力，Tij是湍流应力项

下标i和j轮流以x和y代替， f为科氏力， f=2ω1sin ψ，ψ为当地纬度，ω1为地球自转角速度．

图1为直角坐标向σ坐标转换示意图，由图可见转换之后的地形和自由表面均变为平整面，并沿水深方向对研究区域进行分层．图1(a)中H为基准面到底床的水深，η为基准面到自由水面的水位，D为总水深，图1(b)中σk为σ坐标系中第k层的水深．

(a) 直角坐标系(b)σ坐标系(a) The Cartesian coordinate system(b) The σ coordinate system图1 坐标图Fig. 1 Coordinates

现有的海洋三维模型多为σ坐标转换，在文献[3-4,6]中分别考虑了经σ坐标转换后速度沿σ坐标方向的二阶项

此次改进后的模型，不仅考虑了由流速沿σ坐标方向的二阶项，同时加入了流速沿水平方向变化而引入的高阶项，以及地形和自由表面水平变化引入的高阶项．

设水体不可压缩且满足静压假定和Boussinesq近似，本文从直角坐标系下的连续方程和动量方程出发，推导出在σ坐标系下的连续方程和动量方程，即(x,y,z)→(α,β,σ)．将直角坐标系下方程转换到σ坐标系下：

连续方程为

(5)

运动方程为

(6)

(7)

式中νh为水平黏性系数，νz为垂向黏性系数，ω为σ坐标系下的垂向速度

(8)

其中

在σ坐标系下，沿整个水深(-1,0)积分连续方程，并考虑运动学边界条件，可得积深整体连续方程：

(9)

式中和为整体平均速度，表达式为[17]

(10)

动力学边界条件如下：

自由面

底平面

ω(x,y,-1,t)=0．

海岸边界，法向速度为0；在沿海与开路边界，由实测潮位或流速给出．

为方便起见，以(x,y,σ,t)代替(α,β,σ,tσ)．在σ坐标系下，取一系列与x和y轴平行的平面，利用三角形网格剖分所有平面，设Φi为线型三角元的坐标函数，则本文所取插值函数为

(11)

式中i取为1,2,3分别对应三角形单元中的三个节点，ui(σ,t)，vi(σ,t)，wi(σ,t)，Di(t)，ηi(t)对应的是三角形单元中节点i的值．

将式(11)代入式(5)～(9)中，即可得到有限元方程．与文献[3-4,6]比较，共增加了6个有限元方程系数矩阵，其形式如式(12)所示：

(12)

式中i， j，k取为1,2,3分别对应三角形单元中的三个节点．

2.1.1 设计模型试验

计算水域 平底，长L=3.5 km、宽B=1.6 km，初始水深见表1．

初始条件 全域初始水位为0.0 m，流速为0.0 m/s．

开边界条件 计算域的左侧边界为开边界， 潮位驱动为M2分潮， 振幅见表1， 潮周期为44 712 s； 右侧为自由出流； 其他两个边界为固边界， 固边界的法相流速为0 m/s．模拟总时长为864 000 s．

表1 潮位驱动模型试验数据

Table 1 The tidal driven model test data

experiment123456789depth D/m1010105202222tidal amplitude A/m10.80.60.840.60.50.40.3manning coefficient n0.0150.0150.0150.0150.0150.0150.0150.0150.015A/D0.10.080.060.160.20.30.250.20.15slope ξ000000000

设r=A/D(r=分潮振幅/水深)，试验1～3中r≤0.1，对应一般情况下分潮振幅与平均水深比值；试验4～9中，0.15≤r≤0.3，对应极端情况下分潮振幅与平均水深比值．

2.1.2 潮位驱动垂向平均值

根据Ippen[18]的理论，潮位驱动二维模型试验中，水位和垂向平均流速解析解为

η=Acos(ω2t),

(13)

(14)

式中η代表潮汐水位，m；A代表潮汐振幅，m；ω2代表潮汐角速率，代表垂向平均流速，m/s；D代表水深，m．

图2给出了水位、水深理论值与原模型及改进模型计算值对比的结果．由图可见，原模型及改进模型运算结果与理论值完全拟合，两种模型均能较好地模拟潮位变化过程和水深变化过程．

模型试验中将水体沿垂向平分为若干层，运用改进前后的两种数值模型对同一工况进行模拟，选取3个水平特征点：特征点1(1 000,800)、特征点2(1 700,800)、特征点3(2 400,800)．提取模型结果，根据式(10)计算得到各点垂向平均流速，并与式(14)计算得到的垂向平均流速理论值进行比较，结果见图3．

图2 水位、水深验证Fig. 2 The verification of the water level and the water depth

图3 特征点2的垂向平均流速Fig. 3 The vertical mean velocities of feature point 2

图3给出了试验1～9特征点2的垂向平均流速．从图中可以看出：试验1～3比值r较小，两种模型试验结果与垂向平均流速理论值拟合较好；试验4和6～9比值r较大而水深较小，可明显看出改进模型计算结果与垂向平均流速理论值更为接近；试验5比值r较大同时水深也较大，两种模型计算结果与该点垂向平均流速理论值拟合较好．因此，当比值r较小及比值r偏大，同时水深较大时，原模型与改进模型计算值均与垂向平均流速理论值有较高的吻合度；比值r较大同时水深较小时，改进模型试验结果与垂向平均流速理论值有较高的拟合度．

图4中分别为试验1～9计算得到的涨潮、落潮时刻各特征点理论平均流速、原模型平均流速与改进模型平均流速．从图中可以看出：涨潮与落潮时刻三个特征点的变化趋势相同，随着比值r的增大，流速绝对值逐渐加大；同时，随着比值r的增加、水深降低，原模型计算得到的平均流速较理论平均流速偏小，而改进模型计算结果受比值r和水深变化的影响不大，具有更好的适用性．

(a) 特征点1，落潮 (b) 特征点1，涨潮(a) Point 1 ebb tide (b) Point 1 flood tide

(c) 特征点2，落潮 (d) 特征点2，涨潮(c) Point 2 ebb tide (d) Point 2 flood tide

(e) 特征点3，落潮 (f) 特征点3，涨潮(e) Point 3 ebb tide (f) Point 3 flood tide图4 试验1～9特征点平均流速比较Fig. 4 Feature points’ average velocities for experiment 1～9

(a) 模型分为10层(b) 模型分为5层(a) The model of 10 layers (b) The model of 5 layers图5 各特征点平均流速误差比较Fig. 5 Average velocity errors at the feature points

图5(a)为试验1～9沿水深方向平均分为10层时，计算得到的特征点1～3处原模型垂向平均流速值、改进模型垂向平均流速值与理论垂向平均流速间的误差．从图中可以看出：原模型计算得到的平均流速误差偏大，最大值可达到13.1%，除比值r较小的情况外，剩余试验计算得到的误差相对较大，均在5%～13%之间；改进模型计算得到的平均流速误差较小，最大误差为5.3%，大部分误差保持在3%以内，可达到理论计算与工程计算精度要求，即改进模型计算结果准确度更高．

图5(b)为试验1～9沿水深方向平均分为5层时，计算得到的特征点1～3处原模型垂向平均流速值、改进模型垂向平均流速值与理论垂向平均流速间的误差．与图5(a)比较，当垂向分为5层时，特征点1平均流速误差增加，其余两点平均流速误差相差不大．特征点1位于入口附近，更易受边界影响，增加垂向分层之后减弱了水位变化对流速的影响，计算误差相应减小．因此，综合考虑对结果精度和计算机性能的要求，在后续计算中垂向分层一般取在5～10层之间．

2.1.3 水流速度垂向分布

对于潮位驱动流速沿水深分布规律，可采用对数分布或1/7指数律两种形式．本文根据文献取以下形式：

1) 当0z/h≤0.2时，流速分布表示为对数形式为[19-20]

(15)

式中z0为速度零点的高度，对于水力光滑流，z0=ν/(9u*)；k=0.4为Karman常数，u*为摩阻流速，定义为为水底摩擦力，为摩阻系数，可由Manning系数n表示，CD=gn2/h1/3．

2) 当0.2z/h≤0.5时，流速分布表示为指数形式为[21]

(16)

3) 当0.5z/h≤1时，流速分布表示为[20]

(17)

选取特征点1处模型计算得到的垂向流速分布结果同式(15)～(17)计算的垂向流速分布理论值进行比较，见图6．

图6 特征点1垂向流速分布比较Fig. 6 Vertical velocity distributions at feature point 1

图6(a)、(b)对应第一个涨潮和落潮时刻；图6(c)～(i)对应第二个涨潮和落潮时刻．从图6(a)、(b)可以看出，原模型计算结果与理论值相比偏差较大，即原模型需要更长时间才能达到稳定状态．图6(c)～(e)沿水深方向平均分为10层，底层水深较大，即在底层中同时包含黏性底层和完全湍流层．底层下表面黏滞力由床面剪切应力计算公式确定，上表面黏滞力由湍流Reynolds应力公式计算，与实际浅海潮汐水流运动情况相符．因此，两种模型计算得到的垂向流速分布特征与垂向理论流速分布特征相近．图6(f)～(i)在增大比值r的同时减小平均水深至2 m，随着比值r的增加，原模型计算结果较理论值相比误差增大．特别是近底层，壁面附近流体剪切较为剧烈， 引起的扰动较大， 改进后的模型在考虑流速垂向变化的同时考虑流速水平方向的变化及自由表面水位变化， 与实际水体流动条件更为接近， 计算结果与理论值拟合较好．

由图3～6可知，改进模型较原模型具有更高的精准度和更好的适用性．当水深较大而潮汐振幅相对较小时，可选用原模型或改进模型，两种模型计算结果均能达到较高的精准度，但原模型需更长的运行时间才能达到稳定状态．当水深较浅而潮汐振幅较大时，应选用改进模型，即在模型计算中应考虑水平方向流速变化及自由水面变化．

2.2.1 设计模型试验

根据文献[22]建立8个不同水深、底坡斜率数模试验，试验参数见表2．试验区间为40 m×0.3 m的矩形明渠，由表中可知8个试验的宽深比在3.50～11.11之间，最大底坡坡度为2.7‰，最小底坡坡度为0.09‰．沿水深方向平均分为10层，运行时间为10 h，待模型运行至稳定状态后取特征点(20,0.15)处的流速对数值模型进行验证，结果见图7．

图7中给出了各试验特征点(20,0.15)处的平均流速，原模型与改进模型计算得到的平均流速均与实测平均流速拟合较好，但原模型误差较改进模型误差偏大．改进模型误差均小于3%，而原模型误差最大值接近8%，即改进模型计算结果更为精准．

表2 数值模型试验数据

Table 2 Numerical model experiment data

experimentdepth D/mbreadth B/mbreadth-depth ratio B/Dslope ξfriction velocity v/(m/s)10.0820.3003.6590.000 1000.008 5020.0500.3006.0000.000 5500.014 2130.0720.3004.1670.000 8380.021 4440.0400.3007.5000.000 3140.010 4450.0270.30011.1110.000 0900.003 5960.0850.3003.5290.001 0000.024 7370.0800.3003.7500.002 7000.038 2080.0570.3005.2630.000 7200.017 20

图7 模型试验验证Fig. 7 Verification of the numerical models

2.2.2 试验结果分析

图8为两种模型垂向流速与实测垂向流速对比误差．由图中可以看出，除个别水层外，改进模型垂向流速分布大多数情况下更接近实测数据．在设置的8个试验中，原模型垂向流速分布最大误差均大于改进模型，改进模型垂向流速最大误差为8.21%，而原模型垂向流速最大误差达到了12.17%．可见，改进模型计算得到的垂向流速分布结果更能反映明渠实际流体垂向运动特征．

图9为原模型及改进模型垂向流速最大误差对比折线图．由图可见原模型垂向流速分布最大误差值均大于改进模型计算结果．原模型中除试验2，4，5小于8%以外，剩余试验最大误差均超过9%，而改进模型中除试验6误差超过8%以外，剩余的7个试验最大误差均在7.5%以内．

图8 垂向流速误差分布对比Fig. 8 Vertical velocity errors of the numerical models

图9 最大误差对比Fig. 9 Comparison the maximun velocity errors

图10中为u/u*与u*z/v的关系曲线图，定义无量纲参数U+=u/u*，Z+=u*z/ν．随着Z+的增大，U+呈现出先迅速增大，后趋于平稳的趋势．在黏性底层内，流速变化梯度较大，U+会在短距离内迅速增大；当流速增大到一定数值后趋于稳定，即在完全湍流层内流速变化梯度小．根据改进模型结果计算的U+值较原模型计算结果更接近实测结果，特别是在底坡较小的条件下，改进模型计算结果具有更高的准确度．

综合图7～10数据结果可得看出，改进后的模型计算结果准确度更高，能更好地反映明渠恒定流垂向流动分布特征．改进模型中加入x方向与y方向的高阶项，提高了计算结果精度．因此，当对计算结果有高精度要求，且需要获得精确的垂向流动分布特征时，应采用改进后的模型进行计算，以保证分析结果的准确度．

图10 U+-Z+关系曲线Fig. 10 The U+-Z+ curves relationship

本文根据以上分析计算内容，总结出以下结论：

1) 重新推导σ坐标系下的连续方程与动量方程，得到式(5)～(8)，选用有限元方法取特定插值函数，得到σ坐标系下的有限元方程．推导出的有限元方程较文献[3-4,6]中对应方程增加了6个有限元方程矩阵．

2) 通过建立潮位驱动模型试验比较原模型与改进模型计算结果精准度．根据分析结果得到：改进模型较原模型计算结果精度更高；适用潮位振幅与水深比值范围及水深范围更广；改进模型运算至稳定状态所需要的预热时间更短．因此，在研究某些极端条件下潮位驱动水动力特征时，应采用改进模型．

3) 建立变坡明渠恒定流模型试验，分析原模型、改进模型计算结果与实测数据进行比较，可以确定：改进模型计算得到的垂向平均流速更接近实测值，即改进模型计算结果与实测值间的误差更小；改进模型得到的垂向流速分布更接近实测流速分布．因此，在对水动力特征及垂向流动分布特征结果有较高的精度要求时，应采用改进模型加以计算．

参考文献

[1] PHILIPS N A. A coordinate system having some special advantages for numerical forecasting[J]. Journal of Atmospheric Sciences, 1957, 14(2): 184-185.

[2] 窦振兴, 杨连武, OZER J. 渤海三维潮流数值模拟[J]. 海洋学报, 1993, 15(5): 1-15.(DOU Zhenxing, YANG Lianwu, OZER J. Numerical simulation of three dimensional tidal current in Bohai Sea[J]. Acta Oceanologica Sinica, 1993, 15(5): 1-15.(in Chinese))

[3] 白玉川. 海岸三维潮流数学模型的研究及应用[J]. 海洋学报, 1998, 20(6): 87-100.(BAI Yuchuan. Research on three-dimensional tide current mathematical model of coast and its application[J]. Acta Oceanologica Sinica, 1998, 20(6): 87-100.(in Chinese))

[4] 白玉川, 于天一. 分步分层拟三维水流数学模型及其在廉州湾潮流计算中的应用[J]. 海洋学报, 1998, 20(5): 126-135.(BAI Yuchuan, YU Tianyi. Split step and layered quasi-three-dimensional mathematical model of flow and its application to tide current calculation in the Lianzhou Bay[J]. Acta Oceanologica Sinica, 1998, 20(5): 126-135.(in Chinese))

[5] 陈虹, 李大鸣. 三维潮流泥沙运动的一种数值模拟[J]. 天津大学学报(自然科学与工程技术版), 1999, 32(5): 573-579.(CHEN Hong, LI Daming. 3D numerical simulation of tidal current and sediment transportation[J]. Journal of Tianjin University (Science and Technology), 1999, 32(5): 573-579.(in Chinese))

[6] BAI Yuchuan, WANG Zhaoyin, SHEN Huanting. Three-dimensional modelling of sediment transport and the effects of dredging in the Haihe Estuary[J]. Estuarine Coastal & Shelf Science, 2003, 56(1): 175-186.

[7] 汪守东, 沈永明. 三维对流扩散方程的三种高精度分裂格式[J]. 应用数学和力学, 2005, 26(8): 921-928.(WANG Shoudong, SHEN Yongming. Three high-order splitting schemes for the 3D transport equation[J]. Applied Mathematics and Mechanics, 2005, 26(8): 921-928.(in Chinese))

[8] 郝嘉凌, 宋志尧, 严以新, 等. 河口海岸潮流速分布模式研究[J]. 泥沙研究, 2007(4): 34-41.(HAO Jialing, SONG Zhiyao, YAN Yixin, et al. Study on the tidal velocity profile in the estuarine and coastal areas[J]. Journal of Sediment Research, 2007(4): 34-41.(in Chinese))

[9] WARNER J C, SHERWOOD C R, SIGNELL R P, et al. Development of a three-dimensional, regional, coupled wave, current, and sediment-transport model[J]. Computers & Geosciences, 2008, 34(10): 1284-1306.

[10] LIU Z W, CHEN Y C, LI L, et al. Sigma-coordinate numerical model for side-discharge into natural rivers[J]. Journal of Hydrodynamics, 2009, 21(3): 333-340.

[11] 于守兵, 王万战, 余欣. 基于非结构网格的三维浅水模型[J]. 河海大学学报(自然科学版), 2011, 39(2): 195-200.(YU Shoubing, WANG Wanzhan, YU Xin. 3-D shallow water model based on unstructured grids[J]. Journal of Hohai University (Natural Sciences), 2011, 39(2): 195-200.(in Chinese))

[12] 田勇. 湖泊三维水动力水质模型研究与应用[D]. 博士学位论文. 武汉: 华中科技大学, 2012.(TIAN Yong. Development and application of a three dimensional hydrodynamic and water quality lake model[D]. PhD Thesis. Wuhan: Huazhong University of Science and Technology, 2012.(in Chinese))

[13] 吴毓儒. 长江口深水航道三维水动力与盐度数值模拟研究[D]. 硕士学位论文. 天津: 天津大学, 2013.(WU Yuru. Three-dimensional numerical simulation of hydrodynamics and salinity in the Changjiang Estuary deepwater channel[D]. Master Thesis. Tianjin: Tianjin University, 2013.(in Chinese))

[14] 赵旭东, 梁书秀, 孙昭晨, 等. 基于GPU并行算法的水动力数学模型建立及其效率分析[J]. 大连理工大学学报, 2014, 54(2): 204-209.(ZHAO Xudong, LIANG Shuxiu, SUN Zhaochen, et al. Foundation and analysis of computational efficiency for hydrodynamic model based on GPU parallel algorithm[J]. Journal of Dalian University of Technology, 2014, 54(2): 204-209.(in Chinese))

[15] 李大鸣, 张弘强, 卜世龙, 等. 潮流发电试验场水动力特性数学模型研究[J]. 水力发电学报, 2017, 36(12): 105-120.(LI Daming, ZHANG Hongqiang, BU Shilong, et al. Mathematical modeling of hydrodynamic characteristics of tidal power generation test field[J]. Journal of Hydroelectric Engineering, 2017, 36(12): 105-120.(in Chinese))

[16] 黄牧涛, 田勇. 湖泊三维流场数值模拟及其在东湖的应用[J]. 水动力学研究与进展, 2014, 29(1): 114-124.(HUANG Mutao, TIAN Yong. Three-dimensional lake hydrodynamic numerical modeling and its application to Lake Donghu[J]. Chinese Journal of Hydrodynamics, 2014, 29(1): 114-124.(in Chinese))

[17] 白玉川, 顾元棪, 邢焕政. 水流泥沙水质数学模型理论及应用[M]. 天津: 天津大学出版社, 2005.(BAI Yuchuan, GU Yuanyan, XING Huanzheng. Theory and Application of Mathematical Model for Water Flow Sediment and Quality[M]. Tianjin: Tianjin University Press, 2005.(in Chinese))

[18] IPPEN A T. Estuary and Coastline Hydrodynamics[M]. New York: Mcgraw-Hill Book Co, 1966.

[19] 汪亚平, 高抒, 贾建军. 海底边界层水流结构及底移质搬运研究进展[J]. 海洋地质与第四纪地质, 2000, 20(3): 101-106.(WANG Yaping, GAO Shu, JIA Jianjun. Flow structure in the marine boundary layer and bedload transport: a review[J]. Marine Geology and Quaternary Geology, 2000, 20(3): 101-106.(in Chinese))

[20] 邹志利. 海岸动力学[M]. 北京: 人民交通出版社, 2009.(ZOU Zhili. Coastal Hydrodynamics[M]. Beijing: China Communication Publication, 2009.(in Chinese))

[21] SOULSBY R L. Dynamics of Marine Sands[M]. Thomas Thelford, 1998.

[22] KIRKGÖZ M S. Turbulent velocity profiles for smooth and rough open channel flow[J]. Journal of Hydraulic Engineering, 1989, 115(11): 1543-1561.

BAI Yuchuan1, WEN Zhichao1, XU Haijue1,LIAO Shizhi2, CAO Yonggang2, XIA Huayong3

(1. State Key Laboratory of Hydraulic Engineering Simulation and Safety, Tianjin University, Tianjin 300072, P.R.China; 2. South China Sea Marine Survey and Technology Center, State Oceanic Administration, Guangzhou 510300, P.R.China; 3. South China Sea Prediction Center, State Oceanic Administration,Guangzhou 510310, P.R.China)

Abstract: The 3D hydrodynamic model plays an important role in accurate simulation of the physical characteristics of the ocean. The complex high-order terms are discarded in the traditional σ coordinate transformation due to the limitation of the computer ability, which causes certain errors or calculation distortions for actual complex terrains (or water depth variations). Therefore, the existent σ coordinate 3D hydrodynamic model was modified in order to meet the needs for high-precision calculation results. In the improved model, the complex high-order terms related to the flow velocity, the water level and the terrain introduced through the σ coordinate transformation were comprehensively considered. The specific interpolation function and the combination of the FEM and the FDM were used to solve the complete 3D shallow water model equations in the σ coordinate system. Compared with the existent model, the improved model has a wider range of applications for changes in the bottom slope, the water depth and the tidal amplitude, which could improve the simulation of the vertical flow distribution characteristics under complex water depth changes and promote the accuracy of the calculation results. The calculation error can be kept in a small range in the improved model under some extreme water conditions (with a tidal amplitude to water depth ratio greater than 0.15), and the improved model can reach a steady state in a short time.

Key words: σ coordinate transformation; 3D hydrodynamic model; finite element method; finite difference method; hydrodynamic characteristic

收稿日期： 2018-12-05； 修订日期： 2019-01-05

基金项目： 国家自然科学基金(41576093)；国家重点研发计划(2018YFC0407505)

作者简介：

白玉川(1967—)，男，教授，博士，博士生导师(E-mail: [email protected]);

温志超(1990—)，女，博士(通讯作者. E-mail: [email protected]).

中图分类号： O352; TV131.2

文章编号：1000-0887(2019)08-0840-16

ⓒ 应用数学和力学编委会，ISSN 1000-0887

文献标志码： A

DOI: 10.21656/1000-0887.390344

Foundation item: The National Natural Science Foundation of China(41576093)；The National Key R&D Program of China(2018YFC0407505)

引用本文/Cite this paper: 白玉川, 温志超, 徐海珏, 廖世智, 曹永港, 夏华永. 三维海洋水动力计算σ坐标转换存在的问题及改进方法[J]. 应用数学和力学, 2019, 40(8): 840-855.BAI Yuchuan, WEN Zhichao, XU Haijue, LIAO Shizhi, CAO Yonggang, XIA Huayong. Problems and improved methods of 3D ocean hydrodynamic calculation with the σ coordinate transformation[J]. Applied Mathematics and Mechanics, 2019, 40(8): 840-855.