（水文学及水资源专业论文）波流共同作用下长江口三维水流盐度泥沙数学模型.pdf

摘要 。在河口海岸地区，许多物理现象都是三维的，如盐水楔、泥沙异重流等等，深入研究这 些问题，就必须利用三维数学模型。本文在e c o m 模型和s w a n 模型的基础上，建立了一 套长江口波流共同作用下的三维水流盐度泥沙数学模型。 本文首先对原e c o m 模型中物质输运方程水平扩散项的计算方法进行改进。当水底地 形变化比较剧烈时，原e c o m 模式中水平扩散项的计算就可能产生虚假的流场和物质输运， 在本文中采用返回z 坐标系计算的方法，离散变量利用o r 层二次l a g r a n g e 插值得到。对于 改进方法，文中采用h u a n g & s p a u l d i n g 提出的数值试验进行验证，证明此方法大大地提高 了水平扩散项的计算精度。对于改进的e c o m ，利用长江口实测水文资料进行率定验证， 表明模型对水位流速流向的计算结果与实测值拟合较好，能比较合理地模拟出各测站的盐度 变化过程以及南北槽航道涨落急附近的盐水楔现象。 本文还利用s w a n 模型模拟了长江口不同风场作用下的波高分布，并利用j o s s o n 公式 计算了相应的底部最大流速和底部切应力。结果表明，底部最大流速和底部切应力的大小主 要取决于波高、水深和波浪频率，并随着波高的增大而增大，随着波频和水深的增大而减小。 底部最大流速和切应力对水深的变化更加敏感，而受波高的影响则较小，最大的底流速和切 应力不是分布在波高最大的地方，而是在浅滩上，说明波浪的作用是冲滩淤槽。文中还总结 和研究了波流相互作用机制，波流相互作用将导致两个显著效应：1 ) 引起了涡动粘性系数 的变化，进而导致水流结构的变化；2 ) 增加了水流紊动，引起底部切应力的增加。波流相 互作用的非线性效应取决于波周期的大小，周期越大，非线性效应越明显。 本文建立了波流共同作用下的三维悬沙数学模型，模型包括悬沙输移模块和波流相互作 用模块。其中波流相互作用从两方面来考虑波浪作用：一是涡动扩散系数的变化，采用v a n r i j n 三层模型；二是底部切应力的增加，采用g r a n t m a d s e n 模型和s o u l s b y 非线性模型。模 型把冲刷系数作为率定参数，充分考虑泥沙和河床的局地特性。通过无波情况下长江口实测 含沙量的率定验证，表明模型能基本模拟出各测点的含沙量变化过程和南北槽航道的最大浑 浊带。文中还计算了在波浪作用下的长江口泥沙场，充分表明了波浪对泥沙起悬的作用。 本文最后还研究分析了长江口北槽回淤强度的变化趋势，并利用正交函数法( e o f ) 研 究了回淤强度的时空变化规律。结果充分反映了回淤强度的年周期变化规律和年际变化趋 势，以及人类活动对回淤强度的影响。文中还利用通量分析法对北槽几个测站的流速和泥沙 通量进行深入分析，充分掌握了北槽流速和泥沙的短周期变化规律以及通量过程。 关键词：e c o m ，s w a n ，三维水流盐度泥沙模型，波流相互作用，g r a n t m a d s e n 模型，s o u l s b y 模型，长江口，止交函数法，同淤强度 a b s t r a c t t h r e e - d i m e n s i o n a ln u m e r i c a lm o d e l sc a nb eu s e dt oi n v e s t i g a t em a n yp h y s i c a lp h e n o m e n ai nt h e e s t u a f i n ea n dc o a s t a lr e g i o n s , f o ri n s t a n c e ，s a l i n ew e d g e ，s e d i m e n tg r a v i t y 加w a n ds oo n b a s e d o nt h ee c o mm o d e la n ds w a n m o d e l ，t h i sd i s s e r t a t i o np a y sag r e a tm u c ha t t e n t i o no ns e t t i n gu p at h r e e - d i m e n s i o n a lh y d r o d y n a m i cm o d e lt h a ti n c o r p o r a t e sf l o w , s a l t w a t e ri n t r u s i o n ，w a v e s i m u l a t i o n ，w a v e - c u r r e n ti n t e r a c t i o na n ds e d i m e n tt r a n s p o r tm o d u l e s f o ra s t e e pb o t t o ms l o p ，t h em e t h o dp r e s e n t e di ne c o mf o rc a l c u l a t i n gh o r i z o n t a ld i f f u s i o ni s t e n dt og e n e r a t es p u r i o u sv e l o c i t ya n dt r a n s p o r t h e n c e ，t h ep r e s e n tp a p e rf i r s t l ym o d i f i e st h e c a l c u l a t i o na p p r o a c ht oi m p r o v et h e a c c u r a c y t h eh o r i z o n t a l d i f f u s i o ni s c o m p u t e di n z c o o r d i n a t e s y s t e m ，t h ed i s c r e t i z e dv a r i a b l e sa r eo b t a i n e du s i n gs e c o n do r d e rl a g r a n g e i n t e r p o l a t i o ni ns i g m al a y e r s t h r o u g ht h ev e r i f i c a t i o no fan u m e r i c a lt e s tp r e s e n t e db yh u a n ga n d s p a u l d i n g ，i ti sp r o v e dt h a tt h ec o m p u t a t i o na c c u r a c yo ft h ep r e s e n tm e t h o di sg r e a t l ye n h a n c e d t h ei m p r o v e de c o mm o d e li sc a l i b r a t e da n dv e r i f i e db yt h ei n - s i t uh y d r o l o g i cd a t ao fy a n g t z e e s t u a r y , t h em o d e l e dr e s u l t sa l ei ng o o da g r e e m e n tw i t ht h em e a s u r e dd a t a ，t h et i m e - d e p e n d e n t s a l i n i t yo f e a c hs i t ea n ds a l i n ew e d g e a l eq u a l i t a t i v e l ym o d e l e d t h es w a nm o d e li su s e dt os i m u l a t et h ew a v eh e i g h td i s t r i b u t i o n su n d e rd i f f e r e n tw i n df o r c e si n y a n g t z er i v e rm o u t h t h e nj o s s o nf o r m u l ai sa p p l i e dt oc a l c u l a t et h em a x i i i l u mv e l o c i t ya n ds h e a r s t r e s sa tt h eb o t t o m i ti sf o u n dt h a tt h em a x i m u mv e l o c i t ya n ds h e a rs t r e s sa tt h eb o t t o mh a v e c l o s er e l a t i o n 稍t ht h ew a v eh e i g h t ，w a t e rd e p t ha n dw a v ef f e q u e n c y t h e yi n c r e a s ew i t h i n c r e a s i n gw a v eh e i g h ta n dd e c r e a s i n gf r e q u e n c ya n dw a t e rd e p t h ，b u tt h e y 铀屯m o r es e n s i t i v et o t h ew a t e rd e p t ha n dr a d i a nf r e q u e n c yt h a nt ot h ew a v eh e i g h t h e n c e ，t h em a x i m u ms h e a rs t r e s so f p u r ew a v ei sl o c a t e di nt h es h a l l o ww a t e rr e g i o n ，w h i c hr e s u l t si ns h o a le r o s i o na n ds e d i m e n t r e s u s p e n s i o n t h ew a v e - c u r r e n ti n t e r a c t i o na n dv e l o c i 哆p r o f i l e s 、) l ，i t h i na n do u t s i d et h eb o t t o m b o u n d a r yl a y e ra r ea l s oi n v e s t i g a t e da n ds u m m a r i z e di n t h i sp a p e r f o rs e d i m e n tt r a n s p o r t p r o b l e m s ，i n t e r a c t i o nb e t w e e nw a v e sa n dc u r r e n th a st w os i g n i f i c a n te f f e c t s ：1 ) i tr e s u l t si na c h a n g ei ne d d yv i s c o s i t yc o e f f i c i e n t ，a n d2 ) i te n h a n c e st h eb e dt u r b u l e n c ea n dc a u s e , sa l li n c r e a s e i nb o t t o ms h e a rs t r e s s m e a n w h i l e ，t h e r ei san o n l i n e a ri n t e r a c t i o nb e t w e e nt h ew a v ea n dc u r r e n t w h i c hi n c r e a s e sw i t ht h ei n c r e a s i n gw a v ep e r i o d t h i sp a p e ra l s oe s t a b l i s h e sat h r e ed i m e n s i o n a ls u s p e n d e ds e d i m e n tt r a n s p o r tm o d e lu n d e rt h e w a v e - c u r r e n ti n t e r a c t i o n 。i n c l u d i n gas e d i m e n tt r a n s p o r tm o d u l ea n daw a v e - c u r r e n ti n t e r a c t i o n m o d u l e t h ew a v e - c u r r e n ti n t e r a c t i o nm o d u l et a k e sw a v ee f f e c t si n t oa c c o u n tf r o mt w op o i n 埝o f v i e w ：1 ) ac h a n g eo ft h ee d d yd i f f u s i o nc o e f f i c i e n tu s i n gv a nr i j nt h r e e - l a y e rm o d e l ；2 ) a l l i n c r e a s eo ft h eb e ds h e a rs t r e s su s i n gg r a n t m a d s e nm o d e la n ds o u l s b ym o d e l t h es e d i m e n t m o d e lt a k e st h ee r o s i o nc o e f f i c i e n ta sap a r a m e t e rf o rc a l i b r a t i o ni no r d e rt of u l l yc o n s i d e rt h e l o c a lc h a r a c t e r i s t i c so fs e d i m e n ta n dr i v e rb e d 伽u g ht h ev e r i f i c a t i o no fi n - s i t us e d i m e n t c o n c e n t r a t i o nd a t ai ny a n g t z ee s t u a r yu n d e rt h ep u r ef l o ws i t u a t i o n i ti sf o u n dt h a tt h em o d e l 伽 r e a s o n a b l yr e f l e c tt h es e d i m e n td i s t r i b u t i o ni ny a n g t z ee s t u a r ya n dt h em a x i m u mt u r b i d i t yi n s o u t ha n dn o r t hp a s s a g e i na d d i t i o n t h i sd i s s e r t a t i o na l s ou s e sg r a n t m a d s e nm o d e la n ds o u l s b y 2 m o d e lt os h o wag e n e r a ls e d i m e n td i s t r i b u t i o ni n y a n g t z ee s t u a r yu n d e rt h ew a v e - c u r r e n t i n t e r a c t i o n t h i sp a p e rf i n a l l yi n v e s t i g a t e st h ec h a n g i n gt r e n do f b a c k - s i l t i n gs t r e n g t hi nt h en o r t hp a s s a g eo f y a n g t z ee s t u a r y , a n d 粕a l y z e si t ss p a t i a la n dt e m p o r a ld i s t r i b u t i o nu s i n ge m p i r i c a lo r t h o g o n a l f u n c t i o n ( e o f ) m e t h o d t h er e s u l t sf u u yr e f l e c tt h eb a c k - s i l t i n gs t r e n g t h sa n n u a la n di n t e r a n n u a l v a r i a t i o n s ，a n di n f l u e n c e so fh u m a na c t i v i t i e s af l u xa n a l y s i sm e t h o di sa l s oa p p l i e di nt h i sp a p e r t os t u d yt h es h o r t - t e r mv a r i a t i o no f f l o wa n ds e d i m e n tf l u xi nn o r t hp a s s a g e k e yw o r d s ：e c o m ，s w a n ，t h r e e - d i m e n s i o n a lh y d r o d y n a m i cm o d d ，w a v e - c u r r e n ti n t e r a c t i o n ， g r a n t - m a d s e nm o d e l , s o u l s b ym o d e l , y a n g t z ee s t u a r y , e m p i r i c a lo r t h o g o n a lf u n c t i o nm e t h o d ， b a c k - s i l t i n gs t r e n g t h 3 图1 1 图1 2 图2 1 图2 2 图2 3 图3 1 图3 2 图3 3 图3 4 图3 5 图3 6 图3 7 图3 8 图4 1 图4 2 图4 3 图4 4 图4 5 图4 6 图4 7 图4 8 图4 9 图4 1 0 图4 1 1 图辱1 2 图4 1 3 图4 1 4 图4 1 5 图4 1 6 图4 1 7 图4 1 8 图目录 长江河口分汉形势图 盐水入侵曲线 e c o m 模型水平和垂直方向上的变量配置 e u l e r - l a g r a n g e 方法示意图 模型的计算程序 坐标系统简图：口笛卡儿坐标系g ，z ) ；6 仃坐标系b ，盯) z 坐标下计算水平扩散项有限差分图 插值示意图 数值试验地形图 初始盐度场等值线图 原e c o m 中处理方法的计算结果 改进方法的计算结果 改进方法的计算结果( s = 0 0 5 ) 长江大通年径流总量与输沙总量变化过程 多年各向平均风速 模型计算网格 洪季潮位计算值与实测值比较图 洪季流速和流向计算值与实测值比较图 枯季潮位计算值与实测值比较图 枯季流速和流向计算值与实测值比较图 盐度验证站点位置图 长江口5 月份盐度验证 盐度计算值与实测值相关性图 北槽c s l 站和c s 5 站表底层盐度过程线 南槽涨落急附近等盐度线图( a ) 涨急( b ) 落急 北槽涨落急附近等盐度线图( a ) 涨急( b ) 落急 南北槽断面模型计算盐水通量过程( = 1 0 0 0 p p m ) 8 月份和2 月份表层盐度比较( 实线：8 月份，虚线：2 月份) 8 月份和2 月份底层盐度比较( 实线：8 月份，虚线：2 月份) 长江口2 月份表底层盐度分布( 实线：表层，虚线：底层) 长江口8 月份表底层盐度分布( 实线：表层，虚线：底层) 图4 1 9 不同流量下测站计算盐水过程线 图4 2 0 模型计算5 0 0 0 m 3 s 流量下涨急等盐度线 图4 2 1 模型计算1 5 0 0 0 m 3 s 流量下涨急等盐度线 图4 2 2 模型计算5 0 0 0 m 3 s 流量下落急等盐度线 图4 2 3 模型计算1 5 0 0 0 m 3 s 流量下落急等盐度线 图5 1错误边界条件导致计算失真区域示意图 图5 2长江口波浪场计算网格 图5 3 n n e 向1 0 m s 风场下的波高分布 图5 4n n e 向5 m s 风场作用下的波高分布 图5 5s e 向1 0 m s 风场下的波高分布 图5 6s e 向5 m s 风场下的波高分布 图5 7m 屹向1 0 m s 风场作用下波浪底部最大流速分布 图5 8n n e 向5 m s 风场作用下波浪底部最大流速分布 图5 9s e 向l o m s 风场作用下波浪底部最大流速分布 图5 1 0s e 向5 m s 风场作用下波浪底部最大流速分布 图5 1 1m 忸向1 0 m s 风场作用下波浪底部切应力分布 图5 1 2m 忸向5 m s 风场作用下波浪底部切应力分布 图5 1 3s e 向1 0 m s 风场作用下波浪底部切应力分布 图5 1 4s e 向5 m s 风场作用下波浪底部切应力分布 图5 1 5 波流作用下的涡粘性系数分布( f r e d s o ej & d e i g a a r dr ，1 9 9 2 ) 图5 1 6 涡粘性系数模型 图5 1 7 有波存在与无波情况下的平均流速分布( f r e d s o ej d e i g a a r dr ，1 9 9 2 ) 图5 1 8 波流床面应力的非线性作用( s o u l s b y , 1 9 9 7 ) 图5 1 9 波流共同作用下8 个模型计算的床面平均应力和最大应力比较( s o u l s b y , 1 9 9 7 ) 图5 2 0 波流共同作用下泥沙临界起动s h i e l d s 数( s o u l s b y , 1 9 9 7 ) 图6 1 河口地区泥沙动力机制( m a t h i j s v a nl e d d e n ，2 0 0 3 ) 图6 2 垂向泥沙输运过程( m a t h i j s v a nl e d d e n ，2 0 0 3 ) 图6 3 床面高程( e 表示e r o s i o n ，d 表示d e p o s i t o n ) 图6 4 长江口底沙中值粒径分布图 图6 5 模型含沙量率定图 图6 6 模型含沙量率定相关图 图6 7 模型含沙量验证图 图6 8 模型含沙量验证相关图 图6 9 南槽涨落急附近泥沙等值线图( a ) 涨急( b ) 落急 i l 图6 1 0 北槽涨落急附近泥沙等值线图( a ) 涨急( b ) 落急 图6 1 1 洪季涨急附近泥沙浓度分布 图6 1 2 洪季落急附近泥沙浓度分布 图6 1 3s o u l s b y 模型计算涨急附近的泥沙场 图6 1 4g r a n t m a d s e n 模型计算涨急附近的泥沙场 图7 1 北槽疏浚单元的位置图 图7 22 0 0 帖2 0 0 5 年度全航道平均淤积强度的分布 图7 32 0 0 帖2 0 0 5 年度w 2 - w 3 段航道平均淤积强度的分布 图7 42 0 0 0 。2 0 0 5 年度w 3 - w 4 段航道平均淤积强度的分布 图7 52 0 0 0 2 0 0 5 年间平均回淤强度空间分布 图7 6 e o f 计算特征值分布 图7 7 前8 个模式的p c ( 时间尺度) 图7 8前8 个e o f 模式的e o f s ( 空间分布) 图7 9通量分解各分量关系图 图7 1 02 0 0 4 年5 月c l - c 5 测站流速泥沙通量分析 图7 1 12 0 0 4 年5 月c b l - c s 5 测站流速泥沙通量分析 图7 1 2 北槽水文测点布置 1 1 i 表4 1 表4 2 表4 - 3 表4 4 表4 5 表7 1 表7 2 表目录 长江口潮差沿程变化 大通站月平均流量( 沈焕庭，2 0 0 1 ) 模型率定参数 洪季潮位计算值与实测值比较统计表 洪季流速计算值与实测值比较统计表 2 0 0 4 年5 月c 1 讲测站垂向积分通量 2 0 0 4 年5 月c b l - - e s 5 测站垂向积分通量 i v 学位论文版权使用授权书 本人完全了解同济大学关于收集、保存、使用学位论文的规定， 同意如下各项内容：按照学校要求提交学位论文的印刷本和电子版 本；学校有权保存学位论文的印刷本和电子版，并采用影印、缩印、 扫描、数字化或其它手段保存论文；学校有权提供目录检索以及提供 本学位论文全文或者部分的阅览服务；学校有权按有关规定向国家有 关部门或者机构送交论文的复印件和电子版；在不以赢利为目的的前 提下，学校可以适当复制论文的部分或全部内容用于学术活动。 学位论文作者签名：鸟润峰 秒年；月j 日 经指导教师同意，本学位论文属于保密，在年解密后适用 本授权书。 ， 指剥雠：刍i 丽兰学位论文作者虢马钢峥 萨而年岁月洲e l护年岁月f 日 同济大学学位论文原创性声明 本人郑重声明：所呈交的学位论文，是本人在导师指导下，进行 研究工作所取得的成果。除文中已经注明引用的内容外，本学位论文 的研究成果不包含任何他人创作的、已公开发表或者没有公开发表的 作品的内容。对本论文所涉及的研究工作做出贡献的其他个人和集 体，均已在文中以明确方式标明。本学位论文原创性声明的法律责任 由本人承担。 签名：马铜峰 炒z 年乡月1 日 第一章前言 第一章前言 1 1 概述 根据p r i t c h a r d ( 1 9 5 2 ) 的定义，河口是一个半封闭的沿岸水体，它与外海直接相连，在 这里海水被来自陆地的径流冲淡。在河口地区，潮流、径流、地形、风、波浪和科氏力相互 作用，直接导致其动力过程和地貌过程复杂多变，形成其独特的运动规律。深入了解这些运 动规律，是我们开发和利用河口资源的基础，也是维持河口生态系统可持续发展的必要条件。 长江河口为径流与潮流相互消长非常明显的多级分汉沙岛型中等潮汐河口( 图1 1 ) 。根 据其水动力和地貌形态综合特征，整个河口可分为：大通( 潮区界) 至江阴( 潮流界) 为河 口区的近口段，该河段的造床作用受单向径流控制，河道已发育成相对稳定的江心洲河型； 江阴至河口拦门沙浅滩( 外界大致为余山一牛皮礁一大戟山一线) ，为河口区的河口段，该 河段的造床作用受涨、落潮流往复运动控制，河床演变的基本特点为洪水塑造河床，潮量维 持水深，拦门沙河段为泥沙的汇集区，河口浅滩发育，滩槽之间水、沙交换频繁；长江口口 门以外( 余山一牛皮礁一大戟山一线以东) 至3 0 米等深线，为河口区的口外海滨段，该河 段为河流入海径流和泥沙的扩散沉积区，在地貌上形成水下三角洲，动力特征表现为流向旋 转多变，风浪掀沙及风成余流比较明显。 图1 1 长江河口分汉形势图 长江口的盐水入侵和泥沙运动主要受径流与潮流两股强劲动力的影响。在径流小，潮差 大时，北支有水、盐、沙向南支倒灌，盐水团影响的范围可到达南、北港的吴淞口及堡镇港， 严重影响南支水质和河势稳定。而在南支，盐水界( 盐度为0 5 ) 枯水期可达南港五号沟附 近，北港在六激港附近，崇明六激及南港马家港以下水域直接受外海盐水入侵上溯影响( 恽 第一章前言 才兴，2 0 0 3 ) 。另一方面，由于受到径流、潮流和盐水楔异重流的影响，长江河口四个入海 通道( 北支、北港、北槽、南槽) 都存在碍航拦门沙。其中北支由于分流量较小，几乎全受 潮流控制，所以拦门沙位置深居口内，南槽、北槽和北港的拦门沙都位于口门附近。长江口 拦门沙存在滩长、坡缓、变化复杂的特点，并有洪季淤、枯季冲，小潮淤、大潮冲的变化规 律( 匡翠萍，1 9 9 3 ) 。 本文以改进的e c o m 模型和s w a n 模型为基础，建立了一个考虑波流相互作用的长江 口三维潮流盐度泥沙数学模型，以用于研究长江口的盐水入侵和单独水流或波流相互作用下 的泥沙场，同时采用正交函数法( e o f ) 对长江a ：l l 槽航道近年的回淤强度进行深入地分析 研究。 1 2 国内外研究现状 河口海岸地区由于受潮流、径流、地形、风、波浪和科氏力等诸多因素的作用，其水动 力机制以及随之产生的泥沙和盐水入侵问题都相当复杂。为了模拟和预测这些物理现象的发 展态势，以及对周围生态环境的影响，数学模型是一种有效且可靠的方法。粗略地可以将这 些模型分为( v a nr i j n ，2 0 0 4 ) ： 行为分析模型( b e h a v i o r - o r i e n t e dm o d e l ) ：这类模型包括简单的工程定律、经验或半 经验公式、统计模型以及一些简单的对流扩散模型。这类模型大多对复杂系统作许 多重要的简化和假设，对边界条件只作粗糙的处理，模型具有较大的不确定性。 机制分析模型( p r o c e s s - o r i e n t e dm o d e l ) ：这类模型是以严格的物理机制为基础，建 立相应的、相对完善的数学模型，并用数值解法求得相应的物理量。这类模型也需 要对自然系统和物理现象作一些合理的假设，而这些假设就限制了这类模型适用的 时间和空间尺度。 1 2 1 三维数学模型研究 从l e e n d e r t s e ( 1 9 7 3 ) 建立第一个三维数学模型以来，国内外众多学者已经发展了大量 的浅水流动和输运三维数学模型。其中比较著名的有美国普林斯顿人学的p o m 和e c o m 模 式、佛罗里达大学的c h 3 d 模式、荷兰d e l f t 水力学研究所的d e l f t 3 d 模式、德国汉堡模式 等等。与国外模型相比，国内模型存在分散、缺少权威等缺点( 张修忠等，2 0 0 2 ) 。对于三 维数学模型的研究和问题主要有以下几个方面： 仃坐标变换 浅水流动和输运三维数学模型垂向空间离散通常在笛卡儿坐标系或d r 坐标系下进行。为 了适应不规则的床面地形，在笛 ：- j l 坐标系卜就必须对床面地形作阶梯化近似处理，这就造 成床面计算精度降低。同时，在笛卡儿坐标系下，在给定床面和水面的边界条件时会遇到很 2 第一章前言 大的困难。为此，p h i l l i p s ( 1 9 5 7 ) 提出了o 坐标变换来有效地拟合床面边界，即 仃。三旦 ( 1 - 1 1 )仃l 二1) 日+ ，7 。 式中：z 、，7 和h 分别为笛卡儿坐标下的垂向坐标、自由水面和水深。 从上式可以看出，通过o 变换以后，这个计算域均匀地处于0 和1 之间。这样，计算域 可以在垂向分为相同的层数，带来网格剖分和数值离散上的方便，同时还可以精确地给定床 面和水面的边界条件。对于床面和水面变化比较平缓的水域，仃坐标变换还能保证有较高 的分辨率。正是基于这些优点，许多浅水流动和输运模型都是采用o 坐标系，如e c o m ( b l u m b e r ga e & m e l l o rg l ，1 9 8 7 ) 、c h 3 d ( s h e n gy e ，1 9 8 7 ) 等。但是o 坐标系变换使 得控制方程形式复杂化，带来了很多数值处理上的麻烦和数值误差，许多学者为了减小这些 误差而作了很多的努力。 ( 1 ) 、仃坐标变换的垂向分辨精度 对式( 1 - 1 ) 微分可知， & i d a o t( 1 - 2 ) 由上式可知，在o r 坐标下，物理域的垂向网格尺寸只是水深的函数，所以当床面或水面 变化剧烈时，深水区的层高可能很大，致使垂向网格分布稀疏，难以精确地分辨床面或水面 处的高流速梯度和浓度梯度( 张修忠，2 0 0 2 ) 。 为了解决这个问题，h u a n gw r s p a u l d i n gm ( 1 9 9 5 ) 引进了一个正比于当地水深的 集中因子，使垂向网格的分布能够在床面或水面处集中，以此来改进其垂向分辨率。w a i q w h & l uq m ( 1 9 9 9 ) 考虑到改进的坐标网格分布不能随时间变化，不能精确反映输 运物质浓度梯度的变化，提出了梯度自适应o r 坐标变换，即引入求解变量的梯度及梯度的 变化率作为附加的控制因子，该坐标系能根据浓度的梯度变化自动地调节空间网格。a s a n t o se ta 1 ( 2 0 0 2 ) 采用双o 坐标变换系统，对计算区域垂向分成两个子域，对每个子域进 行仃变换求解，这样上层的子域就不受地形的影响，下层的子域在地形的变化程度上也可 以大大减小，很大程度地消除了地形变化太大带来垂向分辨率差的问题。 ( 2 ) 、仃坐标变换的水平扩散项 仃坐标系( p h i l l i p s ，1 9 5 7 ) 定义为 石i x ，y 。y ，仃。硼z ( 1 - 3 ) 以二维盐度场s 为例，在z 坐标下其控制方程的扩散项为 q 一百1 瓦0 ( 坳，) + i o q , ( 1 - 4 ) 其中，吼； i o s ，留：。k hi o s d xd z 3 第一章前言 经过坐标变换，可以得到 q = 吉 专小百oi a h 石a ( 脚叫+ 万1 【 a a q 了。】 ( 1 - 5 ) 其中，引一以( 等一百c ri o h 石0 s ，。) + 鲁( 0 ，芸) 式( 1 - 5 ) 的表达式中存在看非线性的导数i 贝，给数值离散带米了很多的小便，而且很咨 易产生较大的误差。为了解决这个问题，m e l l o r b l u m b e r g ( 1 9 8 5 ) 考虑床面坡降较缓的情 况，即o h 缸较小时，作如下简化 线。万1 斟0h a 一别 c l q h u a n g & s p a u l d i n g ( 1 9 9 6 ) 给出了另外一种水平扩散项的近似表达式 殴= 嘉( 4 等) 一2 面a 。i ( a h h o a 缸ha a 盯sj + i l 等( 豢) 石。叭( h 1a o s 孑) m 乃 从式( 1 - 6 ) 和( 1 - 7 ) 可以看出，式( 1 - 6 ) 在地形变化较大时就会产生很大的舍去误差，而 ( 1 - 7 ) 则仍然比较麻烦。h u a n g s p a u l d i n g ( 2 0 0 2 ) 采用一种新的方法来计算水平扩散项， 即水平梯度仍在g ，z ) 坐标系中处理，而离散变量的确定由b ，盯) 坐标下计算网格点的插值 确定。 敛t 鳖裔监“鳖铲 m 印 上式是扩散方程的精确描述，产生的只有离散误差。其中( f ，z 。) ，( f 1 ，) 表示在g ，z ) 下同 一水平线下值，s f + 1 嘲= s i + l ，s i 一1 = s i l ，。l ，而s + 1 ，s “吩则由仃坐标计算网格点 的插值确定，h u a n g 等采用了o r k _ 1 ，o r 七，o r 七+ 1 ) 三层上的变量值进行二阶l a g r a n g e 插值，如 此则必须满足“静压一致性”( h y d r o s t a t i cc o n s i s t e n c y ) 条件( h a n e y ，1 9 9 1 ) ： s1 ( 1 9 ) ( 3 ) 、仃坐标变换的斜压梯度力 h a n e y ( 1 9 9 1 ) 最早研究了o r 坐标变换在床面较陡时会带来较大的数值误差，并提出了 “静压一致性”的概念，此后，g lm e l l o r e ta 1 ( 1 9 9 4 ) 研究了o r 坐标下压力梯度项产生的 误差。在b ，y ，仃) 坐标下，压力梯度力为 v ：e f f i l v i 础= g v 叼+ 笥0 p p v 。p 一筹呈d 劝卜盯 m - 哆 p 弓【 d 仃 j 4 第一章前言 当地形变化剧烈、水平方向上密度变化较大时，v 。p 和竺要劝是两个大项，而其 代数和v 口p 一罢要劝可能远小于这两个项本身。由两个大项代数和求小项，很容易产 生与小项相当的计算误差( 朱建荣，2 0 0 3 ) 。g lm e l l o r ( 1 9 9 4 ) 研究并推导了这个误差 e ( 警) 。等警【( 罟) + 孚( 等) + 】【( 埘( 等) 2 】似 其中，6 ，h 表示相邻水平网格点的水深差，b 一p 。g p o ，p 为当地密度减去水平平均密 度，p 一p 一万，p o 为参考密度。m e l l o r x i a o h u a w a n g ( 1 9 9 6 ) 对斜压梯度力的物理机 制作了深入探讨，认为斜压梯度力主要是用来平衡表面的压力梯度力，即所谓的压力补偿。 m e l l o r 等( 1 9 9 8 ) 指出，通过口变换时会产生一类误差，这类误差只会稍微扭曲密度场， 而且随着计算的进行会自行消失，作者把它命名为第一类盯误差( s e f k ) ；对于三维模型， 还会产生另外一类误差，这类误差数值较小，但它不会随着计算而消失，即涡动误差( v o r t i c i t y e r r o r ) ，作者也把它命名为第二类误差( s e s k ) 。对于斜压梯度项产生的误差，e c o m ( b l u m b e r ge ta t 1 9 8 7 ) 采用扣除水平平均密度的方法来减小这种误差，史峰岩、朱建荣等 ( 2 0 0 0 ) 引进局域的概念，提出了扣除局域平均密度层结的思想，并采用返回z 坐标系下 计算斜压梯度力的方法。 系数，主要有以下几个参数：水平涡动粘性系数以、水平涡动扩散系数4 、垂直涡动粘性 系数以和垂直涡动扩散系数k 。一般来说，考虑到水平方向的速度梯度和垂直方向的速 在水平方向上，粘性系数以和扩散系数 可以取为常数( 匡翠萍，1 9 9 7 ：诸裕良等， 1 9 9 8 ；李蓓，唐士芳，2 0 0 0 ) 或采用s m a g o r i n s k y ( 1 9 6 3 ) 公式求解，即 一重。一( ? x 岁 ( 罢) 2 + ( 詈) 2 + 芝l 。( o 妙u ，+ 罢) 2 1 7 2 。1 1 2 ， 其中c 为常数，根据实际情况可取0 0 1 - 。0 5 ，实际应用中取以一a ，。 5 第一章前言 m u n k 型公式 k m k 。( 1 + 卢。r i ) k h a o ( 1 + f l h r i ) 口 ( 1 - 1 3 ) ( 1 - 1 4 ) ( 1 - 1 5 ) 其中，a 。，卢。，口 ，卢h , a o ) 。( - o 5 ，1 0 ，一1 5 ，3 3 3 1 0 4 ) ，彤为分层r i c h a r d s o n 数，有 指数型公式 其中凡，成为常数。 肌一万g 面o p 研l a zp 8 ufa z1 k m | k n e 蹦 k h 一4 , e 9 圆 ( 1 - 1 6 ) ( 1 - 1 7 ) ( 1 - 1 8 ) 双方程紊流模型是建立两个紊动量的微分输运方程，通过求解紊流闭合方程来计算垂向 涡动粘性扩散系数。对于这两个紊动量的选取，有一个共识就是紊动动能k ，而另外一个则 应与紊动混和长度，有关，但不一定是混和长度本身的紊动量( s yj a w & c j c h e n ，1 9 9 8 ) 。 目前根据后一个紊动量的不同，双方程紊流模型有很多模式：k k2 l 模式( k o l m o g o r o v , 1 9 4 2 ) 、k 模式b o ck 1 2j ( s p a l d i n g ，1 9 8 2 ) 、k k l 模式( m e i l o r y a m a d a ，1 9 7 4 ，1 9 8 2 ) 和七占模式k o ck 3 坨lj ( l a u n d e r & s p a l d i n g ，1 9 7 4 ) 。在这些模式中，最常用的是七一模 式( m u s l i mm u i n & m s p a u l d i n g ，1 9 9 7 ；w i h e la l t a ie ta 1 ，1 9 9 9 ；y o n gg l a ie ta 1 ，2 0 0 3 ) ，这 主要是由于占方程可以从n a v i e r - s t o k e s 方程中直接导出而不需要引入任何附加变量。另外 一类常用模式是k k l 模式，e c o m 和p o m 模型即是采用这种模式，在这种模式下，垂直 涡动枯性扩散系数可以通过紊动动能和紊动混合长度乘上一个稳定函数得到，但这种模式也 存在稳定性差的问题( h a n sb u r c h a r de ta 1 ，2 0 0 1 ) ，g a l p e r i nbe ta 1 ( 1 9 8 8 ) 专门研究了这个 问题，并对其中的稳定函数作了有效的改进，并指出紊动混合长度，必须受以下限制： zs z 。i 。0 5 f 3 q ( 1 - 1 9 ) 舯川一一誓) b r u n t - v a i s a h 瓣，g 煽锎叛 三维数模的应用 在河口海岸地区，由丁受到入海淡水径流和海水入侵双重影响，水动力、盐度、泥沙和 6 50 5 n 日 苫 值劫啦 卜 驯 ，il_i、 r z 帆柏 + 口 ，li，、i_-l 暑 4 置 k 第一章前言 污染物的垂向分布具有明显的不均匀性。为了更好地研究和模拟这种垂向分布结构，就必须 采用三维水动力模型，因此在河口海岸地区，三维数学模型得到了广泛的应用，下面是其中 部分学者所作的工作。a