CN112749520B - 大型浅水湖泊三维水动力数值模型建模方法 - Google Patents

大型浅水湖泊三维水动力数值模型建模方法 Download PDF

Info

Publication number
CN112749520B
CN112749520B CN202110074012.9A CN202110074012A CN112749520B CN 112749520 B CN112749520 B CN 112749520B CN 202110074012 A CN202110074012 A CN 202110074012A CN 112749520 B CN112749520 B CN 112749520B
Authority
CN
China
Prior art keywords
water
wind
lake
wave
lcm
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202110074012.9A
Other languages
English (en)
Other versions
CN112749520A (zh
Inventor
吴挺峰
秦伯强
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Nanjing Institute of Geography and Limnology of CAS
Original Assignee
Nanjing Institute of Geography and Limnology of CAS
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Nanjing Institute of Geography and Limnology of CAS filed Critical Nanjing Institute of Geography and Limnology of CAS
Priority to CN202110074012.9A priority Critical patent/CN112749520B/zh
Publication of CN112749520A publication Critical patent/CN112749520A/zh
Application granted granted Critical
Publication of CN112749520B publication Critical patent/CN112749520B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00Computer-aided design [CAD]
    • G06F30/20Design optimisation, verification or simulation
    • G06F30/28Design optimisation, verification or simulation using fluid dynamics, e.g. using Navier-Stokes equations or computational fluid dynamics [CFD]
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00Details relating to CAD techniques
    • G06F2111/10Numerical modelling
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/14Force analysis or force optimisation, e.g. static or dynamic forces
    • YGENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02TCLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
    • Y02T90/00Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Mathematical Physics (AREA)
  • Fluid Mechanics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Computing Systems (AREA)
  • Pure & Applied Mathematics (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • General Engineering & Computer Science (AREA)
  • Algebra (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明涉及大型浅水湖泊三维水动力数值模型建模方法,对建模区域建立三维湖流数值模型LCM和浅水风浪数值模型SWAN;其中LCM模型在σ坐标系中建立,控制方程为连续方程、动量守恒方程和热力学方程,动量守恒方程中包含浪致辐射应力项;为LCM施加边界条件,包括水气边界条件、水土边界条件、热通量,其中水气边界施加风应力,水土边界考虑湖床摩阻;最后基于LCM和SWAN建立风浪和湖流耦合的数值模型WHM‑3D。本发明的建模方法建立的WHM‑3D模型不仅具有模拟波特征要素时空变化的能力;还显著提升了对流速的模拟精度,能够更加准确地模拟盛行风场条件下的太湖风生水体运动过程。

Description

大型浅水湖泊三维水动力数值模型建模方法
技术领域
本发明涉及水环境技术领域,具体涉及考虑波流相互作用的大型浅水湖泊三维水动力数值模型建模方法。
背景技术
以水动力为基础的湖泊水动力-生态数值模型是从全湖尺度研究湖泊生态环境过程的重要手段。通常,湖泊水动力是指湖水受风、密度和出入湖流的影响而在湖盆中的连续运动。水体运动导致的对流扩散也使得水体中的污染物质发生时空迁移。作为模拟这种时空迁移的主要手段,水动力-生态模型被广泛应用到全球著名大湖水生态环境过程研究中,包括:Laurentian Great Lakes(Hoffman&Hittinger,2017),Lake Biwa(Koue et al.,2018),贝加尔湖(Tsydenov et al.,2016),Lake Victoria(Nyamweya et al.,2016).Hutter et al.(2011)认为现代物理湖沼学是在全空间三维数值模型上发展而来的。
不过,目前的湖泊水动力-生态模型以湖流数值模型为基础(Huang,et al.,2010;Mooij,et al.,2010;-Leite and Casenave,2019),而统筹考虑波浪和湖流的模型研究鲜有报道。通常,应用于湖泊水动力及其生态环境效应研究的数值模型大多数是基于海洋环境开发的(Reiss et al.,2020;Soontiens et al.,2019;Ye et al.,2019)。在水深滤波作用下,深水海洋的表面波动影响范围占海水深的比例小,其对水流和生态要素影响也较小。但是,近些年对浅水湖泊的研究发现,风浪主导了浅水湖泊垂向水动力扰动,是导致受侵蚀沉积物悬浮和内源营养盐释放的主要作用力,并对光照、浮游植物和水生植物等均具有重要影响(Qin et al.,2007;Hofmann et al.,2008;Wu et al.,2019)。这意味着过去单纯以湖流为基础的湖泊水动力-生态数值模型并不能反映浅水湖泊真实的水动力结构及受其影响的生态系统。
实际上,风浪在浅水湖泊风生水体运动中占据重要地位。通常,由水气界面输入的风能激发的风生水体运动主要有两种形式:风浪和风生流。两者相伴存在于表层湖水中,并进行相互作用进而影响能量交换(Kirby,1984;Sun,et al.,2006)。Wüest和Lurke(2003)曾将表层湖泊水体定义为wave affected surface layer(WASL)。WASL层厚度随着风力的增强而增大。在浅水湖泊中,由于风浪能量能够轻易抵达水土界面,因此浅水湖泊整个水深范围均可以定义为WASL层。在WASL层中,风浪对水体机械能量的贡献占主导地位(Babanin,2006;Sullivan and McWilliams,2002;Cavaleri et al.,2012),并进一步对风生流产生影响。
目前,较为统一的看法是表面波浪将通过改变水气界面风动量传递效率和水面下水流的应力平衡而作用于水流。作为水流数值模型中最敏感参数(Li et al.,2015),拖曳系数是表征水气界面风动量传递效率的重要参数。以往研究已经证实,表面波浪将改变水气界面的粗糙度,进而影响湖面拖曳系数(Chen,et al.,2020;Foken,2008;Wei,et al.,2016;Wüest&Lurke,2003)。海洋观测表明风拖曳系数是随着风速而线性上升的(Wu,1980)。这是目前数值模型中拖曳系数表达式的主要依据(Zhou et al.,2009)。不过,最近浅水湖泊观测表明风拖曳系数随风速变化呈现先下降后上升的趋势(Xiao et al.,2013)。他们指出造成这种变化趋势与表面波浪发展状态有关。可见,还需要更多的研究来确定浅水湖泊水流数值模型的风拖曳系数。此外,表面波浪还将对下方的水体施加横向切应力,称为浪致辐射应力(Longuet-Higgins and Stewart,1964)。浪致辐射应力将导致下方水体应力平衡的破坏和水流的产生(Noda,1974)。Sun等(2006)推导了适用于三维水流数值模型的浪致辐射应力表达式,但是未将其应用到数值模型中。到目前为止,在湖泊水动力数值模型中还未考虑浪致辐射应力。
发明内容
本发明的目的在于针对目前大型浅水湖泊水动力数值模型研究存在的问题提供一种考虑波流相互作用的大型浅水湖泊三维水动力数值模型建模方法。
为实现上述目的,本发明采用如下技术方案:
一种大型浅水湖泊三维水动力数值模型建模方法,包括:
对建模区域建立三维湖流数值模型LCM和浅水风浪数值模型SWAN;其中LCM模型在σ坐标系中建立,控制方程为连续方程、动量守恒方程和热力学方程,动量守恒方程中包含浪致辐射应力项;
为LCM施加边界条件,包括水气边界条件、水土边界条件、边界热通量,其中水气边界施加风应力,水土边界考虑湖床摩阻;
基于LCM和SWAN建立风浪和湖流耦合的数值模型WHM-3D;LCM为SWAN提供水位和流速输入,SWAN输出的波要素特征值用于计算浪致辐射应力对风生流的作用。
作为本发明的进一步改进,施加于水气边界的风应力基于下式计算:
式中,AV为垂向涡粘系数;H为水深;u,v分别为x,y方向上的流速分量;ρa为空气密度;ρ为水密度;uw,vw为水面以上某一高度处的风速在x和y方向的分量;Cs为风拖曳系数。
进一步的,以uw,vw为变量,分别建立x方向和y方向的风拖曳系数表达式。这样的表达式考虑了动量传递的方向性,使得WHM-3D模拟流速值与实测值相关系数显著增大。目前被广泛应用的风动量传递公式仅考虑了动量量值的传递,并没有考虑动量的方向性,这导致了不同方向的风在x轴和y轴方向上的所产生的动量传递效率是一致的,本发明分别采用风速在x轴方向和y轴方向的分量定义风拖曳系数,较为合理的解决了动量传递的方向性问题。
进一步的,所述风拖曳系数的计算方程基于实测值拟合获取。
进一步的,采用分段函数计算风拖曳系数,在低风速段采用定常值,高风速段采用逻辑曲线表达。水表的风应力系数在高风速和低风速条件下是不连续点,在气动粗糙状态下,风向水体传递动量的效率更高,因此本发明采用分段计算,在高风速段,采用拟合曲线能力更强的logistics曲线来刻画,低风速阶段采用定常值。
进一步的,以7.5m s-1为界,将风速划分为低风速和高风速段。
进一步的,采用分段函数计算风拖曳系数,在低风速段采用定常值,高风速段采用逻辑曲线表达。
作为本发明的进一步改进,求解耦合模型时,以控制SWAN循环计算的时间变量作为WHM-3D的外循环控制变量;在外循环中,先启动SWAN进行当前时间步IT的运算,输出每个计算网格的特征波高,波周期,波长和波向与地理东向夹角,并通过二维数组变量将数据传递给对应的LCM计算网格,用于浪致辐射应力计算;
以SWAN的计算时间步长作为内循环计算周期,在每个内循环周期内,LCM完成N次循环计算,输出每个计算网格x,y方向上的流速分量和水位,传递给SWAN中的二维数组变量,以此作为IT+1步SWAN计算流速和水位输入。
作为本发明的进一步改进,所述LCM模型求解时,采用内外模分离方法求解表面重力波和内重力波的传播。
进一步的,所述LCM模型的求解时,采用带稳定度函数混合长度法计算垂向涡粘系数和热量垂直湍流扩散系数。
作为本发明的进一步改进,SWAN的计算时间步长≥10min,计算频率域为0.04-4Hz。考虑风浪的随机性,风浪特征值往往采用大于10min的统计值表示,基于此确定SWAN的计算时间步长。
本发明分析了大型浅水湖泊中风、风浪和风生流的相互作用,基于浅水湖泊水动力的特点,构建适用大型浅水湖泊的三维水动力数值模型(WHM-3D)。本发明的方法(1)新研发了三维湖流数值模型,并实现湖流数值模型与波浪数值模型的耦合;(2)重新定义浅水湖泊风能输入的表达方式;(3)引入浪致辐射应力来刻画波流相互作用过程。模型的精度评定证实,WHM-3D不仅具有模拟波特征要素时空变化的能力;还显著提升了对流速的模拟精度,能够更加准确地模拟盛行风场条件下的太湖风生水体运动过程。
附图说明
图1是太湖湖泊生态系统研究站(TLLER),五站水位(WL1,WL2,WL3,WL4,WL5)和湖流与气象原位观测站(LCWS)分布图。
图2是笛卡尔坐标系中湖底高程(h),水位(ζ)和水深(H);sigma(σ)坐标系中网格和第i(x方向),j(y方向),k(σ方向)个计算网格的三个流速分量设置。
图3是由SWAN和LCM耦合的风浪和湖流耦合数值模型(WHM-3D)的结构和变量定义和传递。
图4是2015年8月8日0:00至12日0:00间太湖WL1,WL2,WL3,WL4,WL5站处的WHM-3D模拟水位与实测水位变化过程。
图5是2015年8月8日0:00至12日0:00间LCWS站处实表层、中层和底层实测流与WHM-3D模拟流速变化过程。
图6是2015年8月10日12:00时刻的WHM-3D模拟的太湖水位和表、中和底层流场与EFDC模拟的太湖水位和三层流场对比。
图7是WHM-3D和EFDC模拟的2015年8月10日12:00太湖表、中和底层流场流线对比。
图8是2018年12月26日0:00至31日0:00间太湖WL1,WL2,WL3,WL4,WL5站处的WHM-3D1模拟水位与实测水位变化过程。
图9是2018年12月26日0:00至31日0:00间LCWS站处表层、中层和底层实测流速与WHM-3D1模拟流速变化过程。
图10是2018年12月26日22:00时刻的WHM-3D模拟的太湖水位和表、中和底层流场与EFDC模拟的太湖水位和三层流场对比。
图11是WHM-3D和EFDC模拟的2018年12月26日22:00太湖表、中和底层流场流线对比。
具体实施方式
实施例以太湖为例,对本发明的建模方法做进一步阐述。
实施例1
本实施例具体说明本发明的建模方法。
位于中国长江三角洲核心区的太湖(30°55′40″-31°32′58″N,119°52′32″-120°36′10″E)是一个呈浅碟形的大型浅水湖泊。在多年平均水位3m条件下,太湖水面面积2339km2,平均水深1.9m,平均坡度19.7″。太湖位于东亚季风气候区,风场主要受东亚季风影响:多年平均风速3.4±0.19m s-1;4月-8月盛行东南风;其他月盛行西北风。此外,在每年的8月和9月,太湖流域易受台风的影响。尽管我们还没有彻底搞清太湖的水动力结构,但是已经取得广泛认同的观点是太湖全湖尺度的水动力结构决定于风场,出入湖流对全湖尺度水动力结构影响不大(Li et al.,2011;Wu et al.,2018;Zhao et al.,2012)。风场所驱动的风浪和风生流是太湖水动力的两种主要形式,对太湖生态系统具有重要影响,尤其是太湖富营养化和蓝藻水华(Qin et al.,2007;Stone,2011)。
为了满足模型率定和验证需要,分别在夏季(2015年8月1日0:00到12日0:00)和冬季(2018年12月19日0:00到2019年1月1日0:00)开展了两次原位观测。在两次原位观测过程中(图1),利用中国水利部设立的环太湖5个水位站记录水位数据,并在建设有水上观测平台的LCWS站处安装便携式气象站(WXT520;Vaisala Inc.,Finland)和声学多普勒流速剖面仪(ADP)同步记录气象水文数据。
5个水位站的位置如图1中WL1~WL5所示,可以60min时间间隔记录太湖水位时空变化。WXT520被安装在LCWS站处水面以上5m高度处,可以10min时间间隔记录湖面风速、风向、气温、气压和相对湿度变化。夏季湖流观测采用ADP(SonTek Inc.,USA);冬季湖流观测采用XR(SonTek Inc.,USA)。ADP或XR固定在一块设置于湖床上的方形板上。ADP或XR可以每隔30分钟自动记录三维速度剖面图,精度为±0.5cm/s。
以太湖为例,本发明的建模方法流程如下:
一、对建模区域建立三维湖流数值模型(LCM);
本发明在垂向压缩坐标系中建立了一个能够较准确模拟湖泊水位、三维流速、水温和水密度变化的全新的三维湖流数值模型。这个新模型被命名为LCM。LCM模型采用内外模分离技术求解表面重力波和内重力波的传播,并采用带稳定度函数混合长度法计算涡粘系数和热量扩散系数。此外,依据浅水湖泊风生流实测资料,给出了一种新的风拖曳系数表达式,且考虑了浪致辐射应力对湖流的影响。
(1)控制方程
笛卡尔坐标系中的湖流模型的控制方程由连续方程,动量守恒方程和温度守恒方程组成,为了消除湖底地形对湖流模拟精度的影响,在垂直方向上引入σ坐标(图2)。
遵循复合函数求导法则,利用式(A.1-5)将笛卡尔坐标系(x′,y′,z,t′)中的连续方程、动量守恒方程和热力学方程进行坐标变换;
式中ψ为σ坐标系中的u,v,w,ρ,T;w′为笛卡尔坐标系中的垂向流速,m/s。
求得σ坐标系(x,y,σ,t)的方程为:
式中:u,v,w分别为x,y,σ方向上的流速分量,m s-1,m s-1,s-1;h,ζ和H分别为底高程,水位和水深,m;f为科氏力参数,s-1ω为地球自转角速度,/>为地理纬度;Fx,Fy为x,y方向上的浪致辐射应力;ρ,ρ0分别为水密度和参照密度,kg m3;g为重力加速度;AH,AV为水平涡粘系数,垂向涡粘系数,m2 s-1;T为水温,℃;KH,KV为热量在水平和垂向的湍流扩散系数,m2 s-1;Sh,CP为外部进入系统的热量和热容量,J m3 s-1,4179.98J kg-1-1;εU,εv,εT为坐标转换中产生的小项,表达式分别为A.6,A.7和A.8。
已有研究表明使用混合长度法计算的垂向水运动与高阶湍流计算方案无显著性差异(Koue et al.2018)。为了提高计算效率,本模型采用混合长度法计算AV和KV
式中,l,Ri是普朗特长度和Richardson数:
(2)边界条件
当风应力施加在水气边界(σ=1,w=0)时:
式中,ρa为空气密度;uw,vw为水面以上10m高度处的风速在x和y方向的分量,m s-1;Cs为风拖曳系数:据下文模型率定,本模型采用分段函数予以表达,在低风速段采用定常值,高风速段采用逻辑曲线来刻画。
x方向Cs
y方向Cs
对于水温方程的边界热通量为:
式中,ρS为表层水密度,kg m-3;QSN,QLN,QH,QE为向下太阳辐射、大气长波辐射以及湖面感热和潜热。具体计算方法和参数取值参见Koue et al.(2018)。
考虑湖床摩阻的水土边界条件(σ=1,w=0):
此外,模型假设水土边界热通量为0,即:
(3)方程的求解
采用过程分离法将方程(1)-(3)变换为描述表面重力波传播的快过程和内重力波传播的慢过程。对方程(1)-(3)沿垂直方向积分,并定义x和y方向上的平均流速和/>那么方程(1)-(3)变为描述快过程的方程如下:
式中,BU和BV见式(A.9)和(A.10):
由式(2)-(17),(3)-(18),得到慢过程表达式:
式中,DU和DV见式(A.11)和(A.12):
对于快过程,采用隐式差分格式和交错网格(图2和图3)离散式(16),(17),得到计算下一时刻U的方程:
式中,α为格式权重系数:当α=1,式(21)和(22)为显示;否者为隐式。各变量在交错网格上的定义如图2和图3所示。
利用式(21)和(22)求得的下一时刻U来计算ζ和V:
式(21)和(22)在x方向上或者式(23)和(24)在y方向上会形成1个三对角矩阵,可以使用追逐算法求解此矩阵。
对于慢过程的求解同样采用隐式差分格式求解式(19)和(20),得到式(25)和(26):
式中,β为格式权重系数:当β=1,式(25)和(26)为显示,否者为隐式。采用追逐算法求解式(25)或(26)形成的三对角矩阵。
相似的,热力学传输方程(4)也采用隐式差分格式求解:
式中,γ为格式权重系数:当γ=1,式(25)和(26)为显示,否者为隐式;DU的表达式见(A.13)。采用追逐算法求解式(27)形成的三对角矩阵。
依据上述方程,采用Intel Visual Fortran(Intel Inc.USA)编写LCM程序。
二、对建模区域建立浅水风浪数值模型(SWAN)
鉴于风浪对浅水湖泊水动力及生态环境的重要性,选择已被证明适用太湖的Simulating WAves Nearshore模型(SWAN;Booij and Holthuijsen,1999)模拟该湖风浪时空变化(Wang et al.,2016;Wu et al.,2019;Xu et al.,2013):
式中:N为动谱密度函数;t、x和y分别为时间和水平坐标方向;σ1和θ分别为频率和波向;cx,cy,cσ和cθ表示x,y,σ1和θ上的波浪传播速度;S为波能源项,表示了激发、耗散和波波间非线性作用等过程而造成的波能收支。上述方程利用有限差分法中的一阶迎风格式在笛卡尔坐标系统求解。
三、基于LCM和SWAN建立风浪湖流耦合数值模型(WHM-3D);
波流相互作用是一个非常复杂的过程(Mellor,2008)。Longuet-Higgins和Stewart(1964)首先提出了浪致辐射应力的概念。Sun等(2006)推导了地转条件下的浪致辐射应力:
式中,HS为特征波高,m;T0为波周期,s;L为波长,m;θ1为波向与地理东向夹角。
在弄清SWAN程序(version 40.41)和文件结构前提下,通过改写SWAN和LCM程序和文件,建立风浪和湖流耦合的数值模型WHM-3D(图3)。WHM-3D的计算流程为:以控制SWAN的循环计算的时间变量(MTC)作为WHM-3D的外循环控制变量。在外循环中,先启动SWAN进行当前时间步(IT)的运算,输出每个计算网格的HS,T0,L和θ通过二维数组变量VOQ将这些数据传递给相对应的LCM计算网格中(图2和图3),用于浪致辐射应力计算。由于SWAN和LCM平面计算网格相同,这种传递是很容易实现的。其次,以SWAN的计算时间步长600s作为内循环计算周期。在每个内循环周期内,LCM完成60次循环计算,输出每个网格u,v和ζ传递给SWAN中的二维数组COMPDA中,以此作为IT+1步SWAN计算流速和水位输入。
实施例2
本实施例说明太湖具体案例中WHM-3D模型的求解过程及精度评估。
一、WHM-3D设置
LCM计算网格设置为:在水平方向上(x-o-y)上采用1km正方形网格离散计算域,形成72×72个计算网格;垂直方向上(σ)分5层;计算时间步长为30s。α、β和γ均取0.5。Li etal.(2011)和Zhao et al.(2012)均认为出入湖流对太湖流场影响不大,在湖流数值模型中可以忽略。因此,在本次短期强风过程影响下的太湖水动力数值模拟中,出入湖流不予考虑。水气界面输入数据包括气温、气压、云量、相对湿度、风速和风向。其中,水面以上5m高度处的实测风速通过Coastal Engineering Research Center(1984)给出的公式调整到10m高度处(Wu et al.,2018)。LCM重要参数有12个(Table 1),需率定的参数为AH,AV,KH,KV,z0,rs,Cs和CB。其中,AV和KV的值由z0和rs确定,CB值由z0确定。
表1 WHM-3D中LCM所需参数取值
SWAN模型水平网格设置与LCM一致。考虑风浪的随机性,风浪特征值往往采用大于10min的统计值表示。因此,SWAN的计算时间步长确定为600s。SWAN的计算频率域设置为0.04-4Hz,波向计算域为360°,方向步长为6°。源项计算采用第二代风浪模型,考虑了风能输入,水深引起的波破碎,底摩擦和波-波间的作用力等过程。风拖曳系数Cs取0.00133;底摩擦采用Madsen公式,kn为与湖床糙率相关参数,取0.005m。有关这些参数的率定和验证工作已经在另外文献中予以报道(Xu et al.,2013;Wang et al.,2016)。
2015年原位观测数据用于模型的率定,2018年的原位观测数据用于模型的验证。考虑风场变化特征和模型预热,2015年8月8日0:00之前的数据不参与WHM-3D率定,2018年12月26日0:00之前的数据不参与WHM-3D验证。此外,在模型验证阶段,根据式(11)和式(12)计算Cs和浪致辐射应力的WHM-3D输出定义为WHM-3D1,而根据式(37)计算Cs且无浪致辐射应力的WHM-3D输出定义为WHM-3D2。
二、模型精度评定
(1)模拟与实测时间序列比较
使用统计学参数相关系数(r)和平均绝对误差评判模拟水位(MAEWL)和流速量值(MAEUV)与实测值拟合程度(Koue et al.,2018)。
式中,M为模型模拟值;O为观测值;N为样本数。相关显著性水平分别为P<0.05和P<0.01。湖流速用平均值±标准差表示。
湖流流向的率定采用平均绝对误差(MAEUVD;Carvalho et al.,2012):
此外,ArcGIS10.2(Esri Inc.,USA)和Tecplot 360(Tecplot Inc.,USA)被用于处理太湖空间数据和绘制等水位线、流场和流线。
(2)与EFDC对比
Huang et al.(2015)通过比对3个数值模型的模拟结果来论证模型在LakeOntario的适用性。本文也采用相同的策略,通过与太湖EFDC的模拟结果来论证WHM-3D的适用性。作为在全世界范围内浅水湖泊应用最为广泛的模型之一(Chen et al.,2020),EFDC是一个通用的建模包,用于模拟地表水系统,包括河流、湖泊、河口、水库、湿地和沿海地区的三维流动、传输和生物地球化学过程。(Ji et al.,2001).EFDC也是在太湖应用最多的模型之一(Li et al.,2011;Li et al.,2015;Wang et al.,2013),因此本文选择EFDC为参比模型。
EFDC模型的水动力学模型部分是由Hamrick(1992)研发的,其控制方程与式(1)~(5)相同,并且EFDC也采用内外模分离技术求解σ坐标系中的连续方程和动量守恒方程,并采用2阶矩湍流闭合模型来计算垂向湍流粘性和扩散系数(Ji et al.,2001)。在EFDC中,水上风应力可通过如下公式计算(Hamrick,1992;Li et al.,2015):
为便于对比,EFDC模型网格设置与WHM-3D相同。Li et al.(2015)给出了EFDC模型的水平涡粘系数和湍流扩散系数的取值范围均为1-100m2 s-1,并指出这两个参数为非敏感性参数。经与作者交流,EFDC模型的水平涡粘系数和湍流扩散系数最优值为1m2 s-1,糙率高度为0.005m,ws取值0.7。
三、模型参数率定
2015年8月8日0:00至12日0:00之间,太湖上空平均风速为8.6m s-1,主导风向为E发生率37.3%;最大风速发生于10日13:00,风速为13.5m s-1,风向为107.5°。2015年原位观测期间,太湖经历了一次显著的偏东风过程。
在2015年8月8日0:00至12日0:00之间,实测五站太湖水位均值为3.64±0.01m。实测水位最高值4.04m于10日12:00被WL1站记录;最低值3.33m于10日13:00被WL4站记录。
通过对实测值拟合计算,本发明确定了一个全新的Cs的表达式,如方程(11)和(12)所示;AH,KH和z0取与EFDC相同的值;rs取0.01(表1)。在此参数设置情况下,WHM-3D模拟2015年8月8日0:00至12日0:00水位与实测水位极显著相关(P<0.01),相关性均值为0.86(表1),MAEWL五站均值为0.04m。同时,EFDC模拟水位与实测水位极显著相关(P<0.01),相关性均值为0.87(表2),MAEWL五站均值为0.04m。
表2 2015年8月8日0:00至12日0:00间WHM-3D或EFDC模拟五站太湖水位与实测五站太湖水位相关系数(r)和绝对误差均值(MAEWL)
*P<0.05,**P<0.01.
依据数值模型的垂向分层和ADP记录的流速剖面中的单元层对应关系,分别采用表层、中层和底层水体流速对模型开展率定。在2015年8月8日0:00至12日0:00间,LCWS站处实测的三层(表层、中层和底层)流速(图5)均值分别为5.0±3.0cm s-1,5.5±3.5cm s-1和5.4±3.6cm s-1。三层流速实测值与WHM-3D或EFDC模型模拟值显著相关(P<0.01;表3)。其中,WHM-3D模拟的三层流速(图5)均值分别为4.6±3.3cm s-1,4.3±3.0cm s-1和4.3±2.7cm s-1;EFDC模拟三层流速均值为4.4±3.1cm s-1,3.8±2.7cm s-1和3.8±2.5cm s-1。湖流的其他统计参数见表3。
表3 WHM-3D和EFDC模拟的2015年8月8日0:00至12日0:00间LCWS站处表层、中层和底层水体水平流速值与实测水平流速值的相关系数(r)、绝对误差均值(MAEUV)和流向误差均值(MAEUVD)
*P<0.05,**P<0.01.
在8月10日12:00,WHM-3D和EFDC模拟的等水位面延伸相近,均呈现自西北向东南降低的趋势(图6)。两个模型模拟的表层流场流速均具有自东南向西北流动的趋势(图6)。同时刻的流线图更清晰的表现了这种趋势(图7)。同时,尽管中层和底层流场在南部湖区的流态与表层湖流相近,但是中部和北部湖区的水流以自西南向东北流动为主。两个模型模拟的流场区别在于:WHM-3D模拟的流场流速要显著大于EFDC模拟的流场流速。WHM-3D模拟的流场在上风区出现了涡,如胥口湾(表层流场),西山岛西北部(中层流场)。而EFDC模拟流场易在下风区形成涡,如:竺山湾湾口、湾底和梅梁湾湾底(中层和底层流场)。
四、模型验证
2018年12月26日0:00至31日0:00之间,太湖上空平均风速为8.7m s-1,主导风向为NNW发生率30.4%;最大风速发生于26日22:00,风速为13.6m s-1,风向为26.3°。2018年原位观测期间,太湖经历了一次显著的主导风向为NNW风过程。
2018年12月26日0:00至31日0:00间,实测太湖五站水位均值为3.46±0.01m。实测水位最高值3.64m于28日9:00被WL2站记录;最低值3.23m于26日22:00被WL4站记录(图9)。WHM-3D1模拟水位与实测水位极显著相关(P<0.01),相关性均值为0.88(表3),MAEWL五站均值为0.02m。WHM-3D2模拟水位与实测水位极显著相关(P<0.01),相关性均值为0.88(表3),MAEWL五站均值为0.04m。EFDC模拟水位与实测水位极显著相关(P<0.01),相关系数均值为0.87(表4),MAEWL五站均值为0.03m。
表4 2018年12月26日0:00至31日0:00间WHM-3D1、WHM-3D2和EFDC模拟五站太湖水位与实测太湖水位相关系数(r)和绝对误差均值(MAEWL)
*P<0.05,**P<0.01.
在验证期间,LCWS站处实测三层流速(图9)均值分别为3.7±2.0cm s-1,3.5±2.0cm s-1和4.2±2.2cm s-1。三层流速实测值与WHM-3D1模拟值极显著相关(P<0.01;表5)。其中,WHM-3D1模拟的三层流速(图9)均值分别为3.6±2.3cm s-1,3.6±2.2cm s-1和3.8±2.1cm s-1;WHM-3D2模拟三层流速均值为3.9±2.3cm s-1,4.0±2.2cm s-1和4.2±2.2cm s-1.三层流速WHM-3D2模拟值与实测值显著相关(P<0.05;表5);EFDC模拟三层流速均值为3.2±1.8cm s-1,3.4±1.8cm s-1和3.5±1.8cm s-1.三层流速EFDC模拟值与实测值显著相关(P<0.05;表5)。湖流的其他统计参数见表5.
表5 WHM-3D1、WHM-3D2和EFDC模拟的2018年12月26日0:00至31日0:00间LCWS站处表层、中层和底层水体水平流速值与实测水平流速值的相关系数(r)、绝对误差均值(MAEUV)和流向误差均值(MAEUVD)
*P<0.05,**P<0.01.
在12月26日22:00,WHM-3D和EFDC模拟的等水位面延伸相近,均呈现自西南向东北降低的趋势(图10)。两个模型模拟的表层流场流速均具有自北向南流动的趋势(图10)。同时刻的流线图更清晰的表现了这种趋势(图11)。同时,中层和底层流场流速以西北向东南流速为主。两个模型模拟的流场区别在于:WHM-3D模拟的流场流速要显著大于EFDC模拟的流场流速。EFDC模拟的三层流场均在贡湖湾底部形成了一个顺时针涡流。WHM-3D模拟的中层流场在贡湖湾底有一个顺时针涡流。
WHM-3D和EFDC的模型采用相似的数值解法和参数组合,两者最大区别在于垂向涡粘系数、风拖曳系数和浪致辐射应力。由于EFDC采用了更为先进的2阶矩湍流闭合模型(Jiet al.,2001),且其参数是通过对大量样本优选确定的(Li et al.,2015),因此理论上说EFDC能够给出更好的模拟结果。但是事实并非如此。模型验证结果表明,在使用与EFDC相同的拖曳系数和不考虑浪致辐射应力的条件下,WHM-3D2模拟精度与EFDC相近(表4和表5)。这说明,在参数设置相同的情况下,WHM-3D与EFDC具有相近的模拟能力,高阶湍流计算方案并不能改善湖流的模拟结果(Koue et al.2018)。不过,在使用新的拖曳系数表达式和浪致辐射应力条件下,表层、中层和底层流速的WHM-3D1模拟值与实测值的相关系数均值分别较WHM-3D2和EFDC模拟值与实测值相关系数提高48.2%和58.5%,而且模拟的流速量值和流速方向也有改善(表4和表5),尽管水位的模拟精度改善有限。
风拖曳系数是决定风力主导的湖泊水动力模拟精度的关键。在其他参数设置均采用Li et al.(2015)提供的方案下,本研究中EFDC模型的风拖曳系数的计算采用了已被证明适用太湖的方程(36)和(37),WHM-3D则构造了新的风拖曳系数表达式,式(11)和(12)。事实证明,新表达式更加适用于太湖风生流动态的模拟。
此外,水表的风应力系数在高风速和低风速条件下是不连续点,根据现有的报道,当风速高于7.5m s-1时,大气表面为气动粗糙状态;低于3m s-1时,大气表面为气动光滑状态;中间风速为过渡状态。在气动粗糙状态下,风向水体传递动量的效率更高。基于上述实验和观测成果,我们以7.5m s-1为界,将风速划分为低风速和高风速段。在高风速段,我们采用了拟合曲线能力更强的logistics曲线来刻画,低风速阶段采用定常值。
同时,新的表达式考虑了动量传递的方向性,在率定和验证阶段,流速的实测值与WHM-3D的模拟值相关系数比其与EFDC模拟值的相关系数平均提升了44.3%和58.5%。除拖曳系数外,引入的浪致辐射应力对WHM-3D精度改善有重要贡献。以2015年为例,在考虑浪致辐射应力的条件下,表、中和底层流速的WHM-3D模拟值均值分别为4.6±3.4cm s-1,4.1±3.0cm s-1和4.1±2.7cm s-1。这些值要比不考虑浪致辐射应力条件下的WHM-3D模拟值更加接近实测值。此外,在不考虑浪致辐射应力条件下,WHM-3D模拟的三层流速流向的MAEUVD分别为61.2°,54.8°和57.0°(表6),显著大于表3中的MAEUVD。据2018年的计算结果(表6)也能推出相近的结论。可见,浪致辐射应力的引入显著改善了WHM-3D模拟的流速量值和方向。
表6在不考虑浪致辐射应力条件下WHM-3D模拟的2015年和2018年LCWS站处表层、中层和底层水体水平流速值与实测水平流速值的相关系数(r)、绝对误差均值(MAEUV)和流向误差均值(MAEUVD)
*P<0.05,**P<0.01.
两个模型模拟的流场差异也表明了浪致辐射应力的对模拟结果的重要影响。尽管两个模型模拟的水位和流场存在极大的相似性,但是涡的出现位置存在较大差异。以2015年为例,EFDC模拟的中层和底层流场均在竺山湾和梅梁湾存在逆时针涡,但是WHM-3D模拟的流场并不具有相同的现象。这是由于这些涡所在的下风区沿岸带是风浪扰动最为强烈的水域。风浪引起的横向切应力降低了涡形成的可能性。实际上,在不考虑浪致辐射应力的条件下,WHM-3D模拟的中层和底层流场也会在竺山湾和梅梁湾形成逆时针涡。
本实施例中采用WHM-3D模型计算结果,尽管水位模拟精度改善有限,WHM-3D模型模拟的流速与实测流速的相关系数较EFDC平均提升了58.5%,流场流态也更为合理,能够比EFDC模型更好的模拟太湖流场。

Claims (6)

1.一种大型浅水湖泊三维水动力数值模型建模方法,其特征在于,包括:
对建模区域建立三维湖流数值模型LCM和浅水风浪数值模型SWAN;其中LCM模型在σ坐标系中建立,控制方程为连续方程、动量守恒方程和热力学方程,动量守恒方程中包含浪致辐射应力项,得到σ坐标系中LCM模型的控制方程如下:
式中:u,v,w分别为x,y,σ方向上的流速分量;h,ζ和H分别为底高程,水位和水深;f为科氏力参数,ω为地球自转角速度,/>为地理纬度;Fx,Fy为x,y方向上的浪致辐射应力;ρ,ρ0分别为水密度和参照密度;g为重力加速度;AH,AV为水平涡粘系数,垂向涡粘系数;T为水温;KH,KV为热量在水平和垂向的湍流扩散系数;Sh,CP为外部进入系统的热量和热容量;AV和KV采用混合长度法计算;εU,εv,εT为坐标转换中产生的小项,表达式分别为A.6,A.7和A.8;
为LCM施加边界条件,包括水气边界条件、水土边界条件、边界热通量,其中水气边界施加风应力,水土边界考虑湖床摩阻;施加于水气边界的风应力基于下式计算:
式中,AV为垂向涡粘系数;H为水深;u,v分别为x,y方向上的流速分量;ρa为空气密度;ρ为水密度;uw,vw为水面以上某一高度处的风速在x和y方向的分量;Cs为风拖曳系数;风拖曳系数的确定方式为:以uw,vw为变量,分别建立x方向和y方向的风拖曳系数表达式,风拖曳系数的计算方程基于实测值拟合获取,采用分段函数计算风拖曳系数,在低风速段采用定常值,高风速段采用逻辑曲线表达,如下:
x方向Cs
y方向Cs
边界热通量为:
式中,ρS为表层水密度;QSN,QLN,QH,QE分别为向下太阳辐射、大气长波辐射以及湖面感热和潜热;
考虑湖床摩阻的水土边界条件为:
式中,CB为河床阻力系数;此外,模型假设水土边界热通量为0;
基于LCM和SWAN建立风浪和湖流耦合的数值模型WHM-3D;LCM为SWAN提供水位和流速输入,SWAN输出的波要素特征值用于计算浪致辐射应力对风生流的作用。
2.根据权利要求1所述的方法,其特征在于,以7.5m s-1为界,将风速划分为低风速和高风速段。
3.根据权利要求1所述的方法,其特征在于,求解耦合模型时,以控制SWAN循环计算的时间变量作为WHM-3D的外循环控制变量;在外循环中,先启动SWAN进行当前时间步IT的运算,输出每个计算网格的特征波高,波周期,波长和波向与地理东向夹角,并通过二维数组变量将数据传递给对应的LCM计算网格,用于浪致辐射应力计算;
以SWAN的计算时间步长作为内循环计算周期,在每个内循环周期内,LCM完成N次循环计算,输出每个计算网格x,y方向上的流速分量和水位,传递给SWAN中的二维数组变量,以此作为IT+1步SWAN计算流速和水位输入。
4.根据权利要求1所述的方法,其特征在于,所述LCM模型求解时,采用内外模分离方法求解表面重力波和内重力波的传播。
5.根据权利要求1或4所述的方法,其特征在于,所述LCM模型的求解时,采用带稳定度函数混合长度法计算垂向涡粘系数和热量垂直湍流扩散系数。
6.根据权利要求1所述的方法,其特征在于,SWAN的计算时间步长≥10min,计算频率域为0.04-4Hz。
CN202110074012.9A 2021-01-20 2021-01-20 大型浅水湖泊三维水动力数值模型建模方法 Active CN112749520B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110074012.9A CN112749520B (zh) 2021-01-20 2021-01-20 大型浅水湖泊三维水动力数值模型建模方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110074012.9A CN112749520B (zh) 2021-01-20 2021-01-20 大型浅水湖泊三维水动力数值模型建模方法

Publications (2)

Publication Number Publication Date
CN112749520A CN112749520A (zh) 2021-05-04
CN112749520B true CN112749520B (zh) 2024-06-14

Family

ID=75652645

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110074012.9A Active CN112749520B (zh) 2021-01-20 2021-01-20 大型浅水湖泊三维水动力数值模型建模方法

Country Status (1)

Country Link
CN (1) CN112749520B (zh)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114818543B (zh) * 2022-05-20 2024-03-29 中国科学院南京地理与湖泊研究所 一种浅水湖泊典型重现期特征风浪计算方法

Family Cites Families (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2009155300A1 (en) * 2008-06-17 2009-12-23 Lifeflow Technologies, Inc. Fluid monitor
CN104331599A (zh) * 2014-09-30 2015-02-04 江苏省交通规划设计院股份有限公司 一种非结构化网格嵌套波浪数值模拟方法
CN105160162B (zh) * 2015-08-18 2016-12-14 华中科技大学 基于分裂算法的湖泊三维水动力‑水温‑水质模拟预测方法
CN108108911B (zh) * 2018-01-09 2020-09-15 中国科学院南京地理与湖泊研究所 基于水生态系统健康的大型湖泊分区水质目标制定方法
CN108573356A (zh) * 2018-05-08 2018-09-25 镇江海物信息技术有限公司 一种基于三维波-流-沙耦合模型的选择适宜作业海区的方法
CN108764535A (zh) * 2018-05-08 2018-11-06 镇江海物信息技术有限公司 一种大风条件下的船只废弃物对海滨浴场环境影响的预报方法
CN109582996A (zh) * 2018-08-19 2019-04-05 珠江水利委员会珠江水利科学研究院 一种小尺度岸滩剖面与大尺度岸线变化的耦合模拟方法
CN109657418B (zh) * 2019-01-31 2021-11-23 湖北省水利水电规划勘测设计院 一种基于mike21的湖泊水环境容量计算方法
CN109815608B (zh) * 2019-01-31 2019-12-31 湖北省水利水电规划勘测设计院 一种浅水湖泊群水质水量水生态耦合调度分析方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
Reconsideration of winds,wind waves,and turbulence in simulating wind-driven currents of shallow lakes in the Wave-current Coupled Model(WCCM) version 1.0;Tingfeng Wu等;《Geoscientific Model Development》;1-41 *

Also Published As

Publication number Publication date
CN112749520A (zh) 2021-05-04

Similar Documents

Publication Publication Date Title
CN109657418A (zh) 一种基于mike21的湖泊水环境容量计算方法
CN109815608A (zh) 一种浅水湖泊群水质水量水生态耦合调度分析方法
Wang et al. Seasonal circulation and influence factors of the Bohai Sea: a numerical study based on Lagrangian particle tracking method
CN105894106A (zh) 一种海洋模式和气象模式的一体化耦合方法
Chunhua et al. An irregularly shaped warm eddy observed by Chinese underwater gliders
Twigt et al. Coupled 1D–3D hydrodynamic modelling, with application to the Pearl River Delta
CN112749520B (zh) 大型浅水湖泊三维水动力数值模型建模方法
CN101865689B (zh) 一种流域咸潮预测方法
CN110287504B (zh) 一种模拟近岸海水交换规律的耦合模型的建立方法
Zhang et al. Simulation and analysis of Back siltation in a navigation channel using MIKE 21
Clementi et al. Mediterranean sea analysis and forecast (CMEMS MED-currents 2016-2019)
Zhu et al. Comparison and validation of global and regional ocean forecasting systems for the South China Sea
Nayak et al. Tidal and Residual Circulation in the Gulf of Khambhat and its Surrounding on the West Coast of India
CN112434423B (zh) 一种同心圆网格与新型台风场模式结合的风暴潮模拟方法
Raubenheimer et al. Observations and predictions of summertime winds on the Skagit tidal flats, Washington
CN114970393A (zh) 一种高分辨率南海三维斜压环流模型及其构建方法与应用
Zhang et al. Variational estimation of wave-affected parameters in a two-equation turbulence model
Pringle et al. Coupled Ocean Wave-Atmosphere Models for Offshore Wind Energy Applications
Ma et al. Identification and simulation the response of storm-induced coastal erosion in the China Yellow sea
Cawley et al. Sensitivity of a 2-dimensional hydrodynamic model to boundary conditions
Tang et al. A comparative study on wind profiles and surface aerodynamic parameters of typhoons over coastland and coastal sea
Sharbaty 3-D simulation flow pattern in the Gorgan Bay in during summer
CN118364673A (zh) 优化的区域海洋模型的构建方法及海洋区域分析方法
Zhang et al. Targeted observation analysis of the tides and currents in a Coastal Marine Proving Ground
Li et al. The impact of physical processes on pollutant transport in Hangzhou Bay

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant