（环境科学专业论文）厦门海域潮流及网格嵌套边界条件分析.pdf

国家海洋局第三海洋研究所 学位论文版权使用授权书 yiiiii11i1111911t1110111119111113111111811li1 9 i l ly 19 0 9 3 8 9 本学位论文作者完全了解国家海洋局第三海洋研究所有关保 留、使用学位论文的规定，同意本所保留并向国家有关部门或机构 送交论文的复印件和电子版，允许论文被查阅和借阅。本人授权国 家海洋局第三海洋研究所可以将本学位论文的全部或部分内容编入 有关数据库进行检索，可以采用影印、缩印或扫描等复制手段保存 和汇编本学位论文。 保密口，在年解密后适用本授权书。 本学位论文属于： 不保密团。 ( 请在以上方框内打“刀) 学位论文作者签名：- 丰秘骑坞 指导老师签名：嗨哆凝 日期：矿哆年，x 月k 多日 日期：d 7 年修月哆e l 3 1 流体动力学基本方程0 0 0 0 0 0 0 0 0000 000 0000 0 0 18 3 2 网格设置1 9 3 3 方程的差分格式2 0 3 4 动边界处理2 2 3 5 计算区域和参数取值2 3 3 6 3 7 边界条件2 7 大小网格的嵌套2 8 第四章数值模式的验证及讨论2 8 4 1 数值模式的验证2 8 4 2 边界条件对潮流场模拟精度的影响4 5 结语与展望5 0 致谢5 0 参考文献5 1 3 方法求解流体动力学方程 方法的优点在于它引入欧 易求解的优点和隐式良好 进行大小网格嵌套，提高 技术，模拟了大面积浅滩 验证，结果表明数值计算 s 2 两个分潮和包括m 2 、 s 2 等十个主要分潮进行模拟，讨论了分潮的选取对计算结果的影响。 接着在嵌套网格边界上分别使用线性插值、拉格朗日插值方法进行 计算，并对比实测结果探讨边界条件的影响。 关键词：潮流场数值模拟半隐式有限差分边界条件 4 n u m e r i c a ls i m u l a t i o no ft i d a lc u r r e n tf i e l df o rx i a m e n e r $ a n d n e s t e dg r i db o u n d a r ya n a l y s i s a b s t r a c t 1 1 1 ew a t e r sn e a rx i a m e na r ew i d ea n ds h a l l o w at w o d i m e n s i o n a l s h a l l o wf l o we q u a t i o ni su s e dt od e s c r i b et h et i d a l c u r r e n t s ，a n da e u l e r i a n - l a g r a n g i a ns e m i - i m p l i c i tf i n i t ed i f f e r e n c em e t h o dw h i c hw a s i n t r o d u c e db yvc a s s u l l ii sa p p l i e dt os i m u l a t et h ec u r r e n tf i e l di nt h i s a r e a w i t ht h e e u l e r i a n l a g r a n g i a ni n t e r p o l a t i o nf o r m u l a t i o n ，t h i s m e t h o dh a st h ea d v a n t a g e so f e a s ys o l v i n ga n dg o o ds t a b i l i t y t oi m p r o v et h ep r e c i s i o ni ns i m u l a t i o n ，an e s t e d 酊di su s e di n t o n g a nw a t e ra r e aa n dx i a m e nw e s tw a t e r s d u et ot h es i g n i f i c a n ta r e a o ff l o o da n dd e s i c c a t i o n ，av a r i a b l eb o u n d a r yt e c h n o l o g yi su s e d t h e c o m p u t e dr e s u l t sm a t c hw e l lw i t ht h eo b s e r v a t i o n s t h e nt w oi n s t a n c e s a r es i m u l a t e ds e p a r a t e l y , o n ew i t ht h em 2 、s 2c o n s t i t u e n t o n l y , t h eo t h e r w i t ht e nm a j o rc o n s t i t u e n t si n c l u d i n gm 2 、s 2 a n dl i n e a ri n t e r p o l a t i o n a n dl a g r a n g ei n t e r p o l a t i o na r eu s e di nt h ec o m p u t a t i o n st oa n a l y s i st h e i n f l u e n c eo ft h en e s t e dg r i da n db o u n d a r yc o n d i t i o n k e yw o r d s ：t w o - d i m e n s i o n a lt i d a lc u r r e n tf i e l d ；n u m e r i c a lm o d e l i n g ； s e m i - i m p l i c i tf i n i t ed i f f e r e n c em e t h o d ；b o u n d a r yc o n d i t i o n 5 第一章前言 1 1 研究背景 二十一世纪是海洋的世纪，海洋能源、资源的开发与利用，海 洋与全球变化、海洋环境与生态的研究是人类维持自身的生存与发 展，拓展生存空间，充分利用地球上这块最后的资源丰富的宝地的 最为切实可行的途径。我国是一个海洋大国，拥有辽阔的海域面积， 因此科学有效地开发利用海洋资源和能源，保护海洋环境是沿海经 济可持续发展的重要途径，也是成为二十一世纪强国的重要任务。 在有关海洋问题的研究中，水动力起着极其重要的作用，只有 了解水动力的过程，才有可能揭示水质等环境问题的发生、扩展机 制，从而进一步解决生态学问题。水动力是物质输运的基础，流场 对温度、盐度、污染物、泥沙、营养盐、浮游植物和浮游动物等的 输运起着决定的作用，它们相互作用，相互影响，相互关联。目前， 水动力数值模拟技术作为研究水动力学的重要方法和手段己得到广 泛应用，在此基础上对环境问题的进一步研究正在逐步展开。积极 开展对海洋水动力学的研究，并在深入分析海岸水动力学和环境承 载能力的基础上，进行海洋资源的合理开发利用，使生态环境得到 保护并促进区域经济的可持续发展，具有十分重要的现实意义。 1 2 水动力数值模拟技术发展现状与趋势 水动力数值模拟技术被广泛地用于解决各种海岸、河口动力问 题，它最早始于2 0 世纪6 0 年代，是随着计算流体力学( c o m p u t a t i o n a l f l u i dd y n a m i c ，c f d ) 技术的发展而兴起的- - n 综合性模拟技术。到 目前为止，水动力数值模拟历经了较为简单的一维数值模型和可以 求解多项水力参量( 如流速、温度、含盐量、含沙量等) 的二维数值模 型，直到今天大力发展的三维数值模型。在数值模拟不断发展的过 程中，国内外学者提出并完善了众多的数值解法，建立了多种多样 适合各种情况的数学模型，并在实际生产当中得到广泛应用。在科 学研究上，为揭示海洋、河流及海口变化的内在原因，单独和组合 考虑各种动力因素，通过数值试验揭示其变化的动力机制。数值模 拟技术还可以用来揭示许多海洋、近岸、河口物理现象，通过潮流、 泥沙数值模拟海洋环流、河口盐水入侵、河口入海物质输运、甚至 全球变化对河口、海洋的影响等。在社会发展和经济建设上，数值 模拟技术可在工程开工之前，定量预测它们对河口海洋的影响，减 少工程建成后的负面影响和工程的失误。 随着计算机技术的迅速发展和水动力学及相关理论的不断完 善，数值模拟在理论和实践上都发展甚快。目前水动力模型己较为 成熟，在此基础上联合污染物、泥沙、营养盐、浮游植物和浮游动 物等方程，采用三维高分辨率网格，设计质量和能量守恒的差分格 式，建立一个整体模型，是数值模拟技术发展的趋势。 有关潮流、泥沙数值模拟的研究工作，国内外学者从不同方面 做了不同程度的研究。国内有冯士笮、袁业立、方国洪、孙文心、 苏纪兰等海洋学家对物理海洋数值模拟计算进行了许多深入探讨， 李孟国、曹祖德等近年来的学科进展和成果进行了回顾和总结。国 b o t t 和d a v i e s 等学者对近年来三维潮 面的概括。下面简要回顾一下近年来 成果，评价分析数学模型的多种数值 型水库等广阔水域地区，水平尺度远 速、水深等) 在垂向方向上变化要远小 可用沿水深的平均流动量来表示，可 技术。而在另外一些水域，如窄深潮 汐通道、窄深河口地区，有关参量( 如流速、温度、含盐量、含沙量 等) 的垂向变化远大于水平方向的变化，这时可采用垂向二维数值模 拟技术。在解决二维潮流数值模拟的过程中，有多种数值解法可供 选择。下面给出的是几种常用且有效的计算方法，并分析其优缺点。 有限差分法( f i n i t ed i f f e r e n c em e t h o d f d m ) ： 有限差分法是工程中常用的一种离散方法，它将计算域分成有 限个网格，通过差分法求网格结点处微分方程的近似值。由于概念 明确，计算简单，有限差分法在计算流体力学中得到了广泛应用。 针对不同问题和不同精度的要求，人们提出了多种不同的差分格式， 比如显式差分，隐式差分，a d i 法，m a d i 法，破开算子法，准分 析法和特征线法。 a d i ( a l t e m a t i n gd i r e c t i o ni m p l i c i tm e t h o d ) 法( 即交替方向半隐半 显式求解法) 是由p e a c e m a n ，r a c h f o r d 和d o u g l a s 于19 5 5 年提出，该 方法兼备了显式和隐式格式的优点。和隐式格式相比，它虽然在时 间步长上有了限制，但不用每一时间步长都要计算一组大型线性方 程组，可用迭代法求解，节省了计算机内存。和显式格式相比，其 稳定性增强，不会出现数值波动现象。l e e n d e r t s e 最早将a d i 法用 于工程实践。m a d i ( m o v i n ga l t e r n a t i n gd i r e c t i o n a li m p l i c i tm e t h o d ) 法的特点是将水位、水深、流速等变量均布置在同一网格点上，由 此将基本运动方程离散而形成新的差分代数方程组。这种方法既具 有传统a d i 法的优点，又易于编程实现。破开算子法( s p l i t t i n g o p e r a t o rm e t h o d ) ，实质上是一种空间概念上的分步法，即将多维问 题分解成一系列的低维问题，通常是化为一系列的一维问题来求解。 但这种方法也会引入分步误差，并且误差与时间步长成正比。一般 当物理量变化剧烈时采用破开算子法时应使用较小的时间步长。准 分析法是将复杂的方程分解为一系列较简单的子方程，以子方程的 解析解为依据，综合而得最终结果。特征线法是从特征理论出发， 将原始方程进行处理，把偏微分方程转化为沿特征方向的常微分方 程，然后推导出特征差分格式求解，计算过程繁琐，主要适用于涌 潮等不连续现象的计算。 有限元法( f i n i t ee l e m e n tm e t h o d f e m ) ： 有限元法是近十多年来发展起来的一种计算方法，它将一个连 9 的控制体，并使每个网格周围有一个控制体，将待解的微分方程对 每一个控制体积分，得到一组离散方程，其直接的好处是即使在粗 网格下，也保证准确的积分守恒，提高了计算的稳定性。该方法近 年来在水流计算中得到广泛应用。 1 2 2 三维潮流数值模拟 2 0 世纪8 0 年代以前，由于受到计算机处理能力和计算方法的限 制，海岸河口潮流数学模型一般采用二维模型，三维模型还较少用 来解决实际生产问题。后来随着计算机技术和数值模拟技术的发展， 三维模型的应用逐渐成为可能。从n s 方程出发，建立三维水动力 学数学模型，根据静压假定，用水位变量z ( x ，y ，t ) 代替压力变量， 1 0 得到三维浅水方程，使得需求解的数学问题得到简化。l e e n d e r t s e 较早地在实际工程中应用三维模型，随后众多学者对三维模型进行 了大量的研究，取得了丰硕的成果。 国内方面，温生辉等建立了一个三维斜压潮流场模型，加入简 化的紊流闭合模型，并成功运用在厦门西海域整治和开发项目中。 李提来、窦希萍建立了边界拟合网格的三维潮流数学模型，模拟了 长江口复杂边界和地形下的三维潮流变化。刘桦等建立了涉及河口 密度分层效应的三维潮流、盐度数学模型，首次复演了完整的长江 口三维潮流场，并对三维盐度场进行了初步模拟。闰菊等运用一个 三维正压湍流封闭数值模型，成功模拟了胶州湾m 2 分潮的潮汐与 潮流分布。张越美、孙英兰基于e c o m 模型，引入干湿网格法模拟 潮滩涨落的改进，建立了渤海湾三维变动边界潮流模型。杨陇慧等 应用三维高分辨率非正交曲线网格河口海洋模型，模拟了以长江河 口、杭州湾及邻近海区作为整体的4 个主要分潮。 国际上一些研究机构利用资源优势提出了一些技术上较为成 熟，影响力较大的模型，如p o m ( 普林斯顿海洋模型) ，e c o m ( e s t u a r y , c o a s ta n do c e a nm o d e l ) ，荷兰的d e l t t 模型等。我国目前缺乏自行开 发的、国际上影响较大的三维水动力数值模型，一般使用较多的是 p o m 和e c o m ，但是一般的使用者缺乏对模型内在机理的了解，难 以充分利用模型的优势或对模型提出突破性的改进。 1 2 3 动边界处理 在浅水潮流场数值计算中，感潮水域因潮位涨落引起的水陆边 界变化以及波浪爬高等问题，如浅滩涨落潮、洪水波在干河床上传 播，都会涉及到动边界的处理。近2 0 年来，一些学者提出了各种处 理动边界的方法。目前，国内外最常用的方法干湿网格法。其基本 思路是判断计算节点的水深，若水深为负值，则认为该点露出水面， 并令该点的流速为零，或者令该点的糙率足够大，以迫使该节点的 流速为零。窄缝法是处理动边界的另外一种较好的方法，它原理明 确，计算简便、灵活。何少苓等将此方法运用到了杭州湾的潮流计 算中，取得了较好的效果。窄缝法在三维模型中的应用则始于近几 年。江毓武等将窄缝法引入三维数学模型，对厦门前埔围海工程进 行潮流数值模拟，取得了良好的效果。许金电等应用窄缝法对变边 界进行模拟，由此建立三维盯坐标潮流模型和三维盯坐标污染物扩散 模型，并将此模型应用于厦门海域，得出良好结果。温生辉的三维 潮流场模型中，在质量守恒的前提下，在无任何约束条件下，自动 判断干湿边和干湿单元计算，得到一个通用的简明的动边界处理方 法，成功应用于海洋工程数值计算。 1 2 4 潮流数值模拟发展趋势 综上所述，国内外学者对于潮流数学模型的研究虽已取得了丰 富的成果，但是仍有许多问题有待解决，主要表现在以下方面： ( 1 ) 数值模拟的软件化 将水动力数值模拟的全过程设计成软件包，配合良好的前、后 处理工具，降低技术门槛，将使其通用性大大加强，并且通过图形 界面操作，可大大减轻工作量。目前一些通用c f d 软件，如f l u e n t 、 s t a r - c d 、c f x 等，也可以用来进行水动力数值模拟，而它们的附 属程序包，如g a m b i t ，可以方便地根据复杂边界生成非结构网格。 ( 2 ) 开发基于并行计算的高效率高精度模型 过去由于受计算机技术的制约，在实际问题的应用中，只能对 方程进行大量简化，得到一些线性或拟线性方程，然后求得其数值 解，但实际的流动问题大都是复杂的强非线形问题。随着计算机技 术的迅速发展，人们试图在某些问题上放弃假设，直接求解和探索 流动机理问题，如抛开基于统计理论的湍流简化模型直接模拟湍流 问题，在这些问题中，网格数目异常惊人( 达1 0 9 量级) ，因此，开 发大规模高效率的并行算法成为必要。在我国，张理论、宋君强等、 郭照立、王能超等对并行算法进行了一系列探讨，并应用到实际工 程当中。目前，我国己建立了多个“国家高性能计算中心”，如上海 超算中心( s s c ) 等，随着计算机技术的不断发展，并行算法也将突飞 猛进。 ( 3 ) 与新兴学科相结合 数值模拟技术与c a d ( 计算机辅助设计) 、g i s ( 地理信息系统) 等新兴技术相结合也是未来发展的一个重要趋势。g i s 技术的发展 使潮流模拟的可视化大为加强，比如网格能够在g i s 上直接生成， 计算结果能够通过g i s 直接在海图上展现，从而提高了整体效率。 1 2 5 小结 海洋动力学数值模拟是研究海洋潮流、泥沙运动及水质污染扩 散的重要手段和方法，通过此技术并借助于地理信息系统、计算机 辅助设计、以及多媒体可视化等理论和技术的扩展，可以很好地满 足海岸海洋科学的研究需要。 到目前为止，水动力数值模拟历经了较为简单的一维数值模型 和可以求解多项水力参量的二维数值模型，直到今天大力发展的三 维数值模型。在数值模拟不断发展的过程中，国内外学者提出并完 善了众多的数值解法，建立了多种多样适合各种情况的数学模型， 并在实际生产当中得到广泛应用。水动力数值模拟发展的重点有： 网格的生成( 前处理) ，主要研究在具体几何特性，流动特性及流场求 解算法下，如何自动生成网格并使网格的光滑性、正交性、网格单 元的正则性以及网格的疏密等特性达到最优；数学模型的求解( 求解 器) ，主要研究的是如何快速，准确地求解离散方程，得到所需的数 值解；输出结果的可视化( 后处理) ，数值模拟可视化研究内容主要是 如何利用计算机的图像处理技术，将数值模拟结果在终端上显示出 来，使计算结果明显、直观、生动。 随着水动力学理论以及计算机技术的不断完善，潮流、泥沙数 值模拟在理论上和计算技术上也不断向前发展。可以预见，今后潮 流、泥沙数值模拟将朝着高效、高精度、可视化、软件化等方向发 展。 1 3 本文的主要工作 本文的主要工作包括： ( 1 ) 选取二维浅水波方程组为控制方程，采用vc a s s u l l i 半隐 式有限差分方法离散潮波方程，建立了二维潮流数学模型。并用此 模型对十个主要分潮( m e ，s 2 ，k 2 ，k l ，o l ，p l ，l v h ，m s 4 ，m 6 ， n 2 ) 作用下厦门海域潮波运动进行了数值模拟。 ( 2 ) 对九龙江河口湾、厦门西海域和同安湾海域进行4 0 0 m 和 5 7 m 细网格计算，提高模拟精度。 ( 3 ) 对比了选取m 2 、s 2 两个分潮和包括m 2 、s 2 等十个主要 分潮情况下计算结果的不同。边界条件分别采用线性插值、拉格朗 日插值处理，对比其对流场精度的影响。 第二章厦门海域水动力环境概述 为了充分了解厦门海域水动力环境，我们选用的外侧开边界取 为金门岛南北两段，北段开边界为围头至金门岛的北碇头断面，南 段为金门旧城岸边至流会断面。内侧开边界取为九龙江两条干流南 港和中港的上游( 在石码附近) 会合处附近，北港上游为闭边界。 本文讨论的海域即外侧和内侧开边界之间的全部海域( 见第三章图 2 1 ) ，总面积约为1 0 0 0 k m 2 ，包括九龙江下游入海河段( 北港、中港、 南港) ，九龙江河口湾( 九龙江入海口至鸡屿南北向断面) 、厦门外 港( 鸡屿断面以东、大担至岛美断面以西，鼓浪屿和厦门岛南岸经 白石至大担以南) 、厦门西港( 厦门岛以西，嵩鼓水道和厦鼓海峡以 北，高集海堤以南) 、同安湾( 五通至沃头以西) 、厦门东侧海区( 包 括厦门东侧水道和金门水道) 、厦门湾口区( 大担至岛美断面和金门 至流会断面之间) 、围头湾( 包括围头湾和大嶝南侧海区) 。 太平洋潮波进入东海和南海以后，从台湾海峡南北两端传入台 湾海峡。海峡西侧的福建沿岸，潮波自北向南传播，经过上述两段 外侧开边界传入本海域。南支从金门岛东南侧至流会之间传入，后 分为两支，一支通过厦门东侧海域进入同安湾，另一支经过厦f - j # l - 港分别传入厦门西港和九龙江河口湾，并溯河而上至北溪的江东桥 和西溪的榜山桥闸。北支经金门岛北侧进入围头湾，涨潮时在大嶝 南侧海域与进入同安湾的分支会合。 由于各局部海区所处位置不同，地形各异，因而它们的水动力 环境存在着较大的差异。厦门海洋站和石码站均有多年潮汐观测记 录，国家海洋局第三海洋研究所在本海域也设置过多处短期潮位站。 所有短期测站均与厦门海洋站进行了对比分析。这些观测比较全面 地反映了潮波传入本海域后的运动情况。外侧开边界处，潮时北早 南迟，进入本海域后越向湾底潮时逐渐推迟，潮差则逐渐增大；进 入九龙江下游河段后，潮时迅速推迟，河口附近与上游区潮时相差 约4 5 分钟，潮差则迅速减小并且落潮历时大于涨潮历时的情况越向 上游越显著。 厦门海域的潮汐性质为正规半日潮，厦门海洋站多年平均潮差 为3 9 8 m ，平均大潮差为4 9 5 m ，平均小潮差为2 8 5 m ，涨落潮历时 几乎相等。 由潮汐振动形成的潮流，是厦门海域的主要水流动力。潮流的 性质为正规半日潮流，形式以往复流为主( 潮流椭圆半径不大或很 小) ，主要流向多与岸线或水下地形的走向一致。转流时刻多数在高、 低潮时刻附近，停流和憩流时间不长( 一般不超过十几分钟) ，流速 最大出现时刻常在高、低潮后3 小时左右。表、底层转流时刻通常 略有差异。流速的垂向分布，以次表层或表层最大，向下逐渐减小， 接近底层则迅速减小；浅水区流速的垂向变化不显著。国家海洋局 第三海洋研究所在厦门海域进行过多次较大规模的全潮周日海流定 点观测，这些观测表明，本海域的强流区位于邻近外海的围头湾和 厦门湾口区的深水部分，潮流通道区( 厦门东侧水道、厦门外港) 和狭窄水道( 如东渡航道、厦鼓海峡、海门岛南侧) 以及九龙江河 口湾的浅航道也属强流区。弱流区出现在内湾( 宝珠屿海区、同安 湾) ，沿岸浅水区、大嶝南侧海区。潮间带浅滩为漫滩水流，流速一 般不大，仅在潮间带的退潮沟内观测到退潮时可以出现约一小时的 较强退潮流。大潮期间，强流流区流速可以超过l m s ( 最大可达 1 8 m s ) ，弱流区最大流速约0 3 - - - 0 5 m s ，两者差别甚大。 综上所述，厦门海域虽然整体上属于潮差较大的强潮区，但各 局部海区的地理位置不同，岸形和水下地形各异，造成各个局部海 区的水动力环境有明显的差别。 第三章数值模式的建立 3 1 流体动力学基本方程 潮波运动基本方程包括运动方程和连续性方程。前者是牛顿第 二定律在潮波运动中的特定形式；后者是在假定海水为不可压缩流 体的前提下反映海水过程所应遵循的质量守恒定律。这些基本方程 表述了潮波运动过程中，在引潮力、摩擦力、和地转等因子作用下 的波动形态和潮位变化与潮流变化之间的联系。在直角坐标系下， 采用右手坐标系，坐标原点取在平均海面，z 取垂直于等重力势面向 上为正。 z u x 厦门附近海域辽阔，水浅且湾口较开阔，潮差较大，潮流流速 u 、v 水平变化远大于垂向变化，可以采用二维浅水波方程组描述其 运动： r 宴+ “丝+ v 塑= 一g 里+ 和一，“ i 万栅瓦w 瓦= 一g 瓦+ 声吖“ j 加 加抛如 ， 弋瓦栅瓦+ v 面= 一g 面一声一w ( 1 ) i 拿+ 皇坠型+ 业型= 0 西锄加 舯，- 繇 x 、y 为平面直角坐标，u 、v 为x 、y 方向流速分量，t 为时间，h 为 平均水深，z 为潮位即瞬时水位，h - - h + z 为水深，f 为柯氏系数，r 为底摩擦系数，g 为重力加速度，e 为谢才系数，e = 日“6 n ，n 为 曼宁系数，本次研究n 取0 0 2 。 流体动力学方程组可用两种形式表达，即e u l e r 形式和l a g r a n g e 形式，两种形式都用来描述同一物理过程，它们实际上是等价的。 但是在数值计算上两者却有很大差异。由于拉氏坐标是建立在流体 质团上，跟随质团运动：随体导数丢说明沿流线相应物理量随时间 的变化率。拉格朗日时间导数鲁和欧拉时间偏导数箬之间存在全 微分和偏微分关系： 警= 署+ ( 云鲫d m 其中云= 甜f + 巧，删= 丢f + 昙歹 用拉格朗日时间导数表达方程组，不出现输运项( 云g r a d ) 五( 或；) ， 原方程组整理为： r 害= 竿叫妄咿删 j 害：, g l t + l - 一v i t ：一g 妻一f u 一州 ( 2 ) 石2 丁一g 万一州 u j l 妄= 竺= 一z r 罢+ 万0 v ，一r 警+ 等， 3 2 网格设置 uv z 设定网格位置如图1 所示： 按照示意图所标i + 1 2 , ji , j + l 2 分别表示当前网格u ，v 速度所在位 置，水位z i ，i 位于网格中心。u ，1 ，1 分别是( k + 1 ) at 时间层速 度和水位，u k ，v l 【，z k 分别是k a t 时间层的速度和水位，是流线上p 点在k a t 时间层流点所在流速和水位。 3 3 方程的差分格式 本文采用vc a s s u l l i 半隐式有限差分方法。首先用 e u l e r - l a g r a n g e 方法改写方程左边项为全微分项( 引用e l 插值公 式) ，可以减少方程中迁移项的差分离散过程人工数值扩散。另一方 面在方程右边项中，对动量方程水位梯度和连续方程的速度采用隐 式差分，使整个浅水潮波方程中质量方程在每个时间步长内求解过 程变成只求一个未知数( 水位值) 的代数方程组，这个对角线占优 势的线性方程组的解存在、唯一，且稳定。 根据网格配置，微分方程变为差分方程： ，= 嘭 ，一g 对a t , , k + l 一右1 ) +- r 抄u k ， g 考( 兹- 一菇；1 ) + f ( 一 z * * l = f z j k j 。i j l x x t m + - e 厂印一乇每 一尝, 巾k + l _ q 印一号( 幺小照吨v k + l 圭) 相邻四点u 、v 的平均值。 插值公示求得 及硗为显性非线性有限差分因子，由欧拉- 拉格朗e t 川，_ ，= o - a ) w ：j + 彳皇l 。_ ，】( 1 一b ) + ( 1 一彳) w 0 一l + 彳畦l ，一l 】b 其中彳 m 。可以用z j ，+ l ，2 _ ，v s + l ，2替换。 方程组( 3 ) 由于在动量方程中的水位梯度和连续方程中的速度采 用隐式求解，各个方程中的迁移项采用欧拉拉格朗日插值公式的显 式求解，经过整理变为只含一个未知数( 水位) 的数值方程 一g 等棼c 一，一 a t 2 一g 矿 z j 锄u 卜皂 l + + h f 一2 ， ( 右一如i ) 】 心k 川+ l 带鲁( 铲k + 1 1 ) 】l 锄一岩) 一稿( 岩气o 】 2 1 、- ， 一2 + 矿u 一2 i r ， 一 一2 卜 一 l 一2 凡 = 一2 p矿 l 一2 护 一 l 一2 “p 沪l 竺妙卜2 矿u = 8 竺血 一2 矿“ = 出 一2 七 一r ， +i 出 l 一2 + i r = 一尝c 巧z ! j + 争hc 嘭扣+ 够，一专z k + 争h t c 凡三扣+ 够， 一每c礴g；：j1hc q ，一等c 蜕p 吲， ( 4 ) 这个代数方程组是一个五对角线性方程组，只要保证： 式，+ q 广。 就是一个对角线占优势五对角线性方程组，它只有唯一解，而且它 的稳定性不取决于传播速度硒，由于差分方程是五对角格式，因 此求解z i ，j 值仅与它上下左右邻近四个网格的水位值有关，于是计算 时可以简单地以结点i ，j 为中心组成一个五点格式逐点迭代。全流 场各个网格上水位迭代过程的结束，由计算精度加以控制： m a x l z i l ；1 一吒 o ， 称( f ，) 单元此时刻为湿单元。 x 向干湿边的定义： h i n 2 吩+ i ，2 。+ 弓“，2 ，兰j 骂 o 且耳+ l j o ， + 脚+ 弓，当 o 觚+ - o ，( 5 ) 囊+ l ，2 + 刁“，当e 0 且县“， 0 ， 0 ，当县 0 ，称单元界面( i + 1 2 ，) 此时刻为湿边。y 向干湿 边的定义类推。在无任何特殊条件约束下，可自动判断干湿边和干 湿单元，并保证质量守恒。 3 5 计算区域和参数取值 图2 1 2 3 为计算区域图。坐标的y 方向为正北方向顺时针旋 转4 7 0 。九龙江南港和中港的上游( 在石码附近) 取为开边界，北港 上游为闭边界。同安湾湾项河流开边界为石浔支流和浦声支流；外 海流体开边界在围头一金门及金门一流会断面。整个海域采用8 0 0 m 8 0 0 m 等间距网格，时间步长a t = 3 6 秒，共有2 5 8 8 个网格点。嵌 套九龙江河口湾、厦门西海域和同安湾4 0 0 m 网格，西海域和同安 湾5 7 m 网格。 图2 18 0 0 m 网格计算区域图 ( 注：虚线位置为4 0 0 m 网格边界) 图2 24 0 0 m 网格计算区域图 ( 注：虚线位置为5 7 m 网格边界) 图2 35 7 m 网格计算区域图 3 6 边界条件 外海流体开边界采用强制水位，利用围头，流会测站一个月实 测潮位资料，求得十个分潮( m 2 ，s 2 ，k 2 ，k l ，o l ，p l ，m 4 ，m s 4 ， m 6 ，n 2 ) 的调和常数，各分潮组成协振潮水位。水位控制为时间的 已知函数，其形式为： 1 0 z = e , hi c o s ( o - f t + v o i + 1 i l i - 9 3 z 为水位，、e 分别为分潮的调和常数，q 为分潮的角速率， 1 l d ，为分潮格林威治天文初相角，、彳为分潮的交点订正角和交点 因子。为使模拟潮流场与实测资料同步，在水位表达式中，代入每 个分潮2 0 0 2 年5 月1 日的交点因子z 和格林威治天文初相角+ ， 预报出5 月1 日至5 月3 0 日围头，流会的潮位，并进行插值，作为 开边界强制水位。这是模式中潮流运动的驱动力。 外海开边界处的流速边界条件为： 0 u 2o 锄 却 2o 勿 九龙江河口流体开边界水位为：z - 风+ 彳，e , o s ( 2 刀r t ，一g o ) ，式 中风为附加水头，g o 为潮迟角，彳蔚为潮振幅，以此代表九龙江上游 河段的实际潮涨、潮落。在潮振动上再附加一个水头用来模拟平均 迳流量，附加水头按平均年迳流量3 2 1 m 3 s 模拟。同安湾湾项流体 开边界按平均年迳流量1 3 m 3 s 模拟。 河口流体开边界处流速边界条件为： 鱼+ 里盟：o ，y ：o 西缸 2 7 网格交界处的水位和流速数据由上一级网格的数据经插值方法得 到，并分别进行潮流场模拟。 第四章数值模式的验证及讨论 4 1 数值模式的验证 为了检验潮流场的可靠性，在5 7 m 网格下对计算结果进行了潮 位、潮流的验证。经过与实测资料的对比验证，可以看出数值计算 的结果是比较可靠的。 4 1 1 潮位比较 采用了国家海洋局第三海洋研究所2 0 0 2 年5 月在高集海堤附近 和澳头附近设立的临时潮位观测站所做的观测数据，测站位置如图2 所示。实测和计算的半个月潮位对比曲线见图3 ，由潮位曲线可以看 出两站的实测与计算结果吻合良好。 图3 潮位实测值与计算值曲线 4 1 2 潮流比较 潮位实测期间在厦门西海域设立m 1 、m 2 、m 3 、m 4 和m 5 测 流站，在同安湾设立n 1 、n 2 、n 3 、n 4 和n 5 测流站，位置如图2 所示，记录了2 0 0 2 年5 月1 4 日1 3 ：0 0 - - - 5 月1 5 日1 4 ：o o 大潮期 间每个站的海流数据，直接取逐时实测潮流与计算结果加以比较， 见图4 。由潮流曲线可以看出实测结果与计算结果流速和流向比较一 致，大部分时间流速大小相近，转流时流向存在一定误差。以同安 湾为例，n 1 站涨潮流向约为3 3 0 。，落潮流向约为1 5 0 。，n 5 站涨 潮流向约为2 7 0 。，落潮流向约为7 0 。n 3 站计算的最大垂向平均 流速达0 7 7 r n s ，为5 个站中最大，n 1 站达o 7 2m s ，n 2 站达o 7 0m s ， n 4 站达o 4 7r n s ，n 5 站达0 3 9m s ，基本接近实测的最大流速。总 体上说，潮流计算的结果是可信的。 t 向 噬) 舡 ( c m s ) t 向 暖】 1 l g 01 4 ：0 0t 1 1 目01 8 ：0 02 0 e 0 0 2 2 1 0 0哪 2 ：0 0枷t 9 08 c 0 01 0 ：0 0 枷1 4 ：0 0t 2 ：t 4 j 0 0t 1 1 h 0 0t l l c a 02 0 1 0 0 2 2 0 0哪 2 ：o o 4 o o e ：0 0螂1 0 ，0 01 2 | 0 01 4 , 0 0 t 向 $ 2 0 01 4 ：q 0i i i b 0 1 0 螂，啪2 2 0 0g 0 0枷枷蜘蛐 锄1 2 ：a 01 4 0 瑚1 4 ：0 01 4 ：1 0 2 0 1 0 0z 嘲哪枷枷螂嘲 旧1 2 a 01 4 ：9 0 艋 ( c a s ) t 2 0 0i g o1 8 ：a 02 0 1 0 0 嗣娜嘲硼4 ：0 0螂枷1 q a o1 2 0 01 4 4 0 01 2 ：0 0 t t 0 0 ：3 d o o 叠：0 ：o o枷枷蛐嘲1 q 0 01 2 ：0 0 堋 髓 刮s ) l 向 暖) 1 2 1 0 0t 4 ：0 0t l u 0 0t l l ：0 0 捆哪z 啪 q c o 2 口0枷螂枷1 0 ：0 0 枷t 4 ：0 0 l 世 向 镀) 髓 i 州1 ) i ：u 0 01 4 】0 0l 螂5 1 1 o3 0 ：a 0z 啪咖珊4 ：0 0枷蜘 1 0 ( o o 1 2 1 0 01 4 d 0 0 9 0 ， 嘲1 7 ：0 0l l p ：o o翻蕾i n n枷：7 9 ：0 01 舯 1 c o1 1 0 0 螂1 7 , 0 0 荆：d o 蕾i c o螂枷y ：0 0螂1 1 ：0 0 荫啪 t 向 幢) 1 1 ：0 0 枷1 6 0 01 7 ：0 0 螂钠姗1 ：0 0埘细7 ：0 00 c c o1 1 n n 螂 1 0 0l c 0 0 t 1 7 ：0 0t l i e q 0甜埔硼i 埔& g o7 ：0 0螂 1 1 】0 0 住一 瞄 酗 图4m i m s 站和n 1 n 5 站大潮实测垂线平均流速与计算流速比较 4 1 3 潮流场分析 图5 图7 分别为8 0 0 m 、4 0 0 m 、5 7 m 网格厦门海域大潮潮流场， 图中阴影表示浅滩的干出。以下对整个潮流场的流态做简要说明。 从流场图中可以看出流场分布情况与第二章厦门海域水动力环境概 述中实测潮流平面分布趋势很相似： ( 1 ) 外海潮波通过流会一大金门断面和北碇岛一围头断面分别 从南、北传入厦门海域，南面的潮波分成两支，一支通过厦门东侧 海区进入同安湾，另一支经过厦f - j 多, t - 港后又分为两支分别进入厦门 西港区和河口湾直至九龙江下游感潮河段。从澳头一五通断面进入同 安湾的外海潮波，涨潮时，主流沿着同安主航道上溯，而后在鳄鱼 屿以南附近水域分成三股涨潮流，一股沿集美次航道到达高集海堤， 一股沿洪塘头流向潘涂滩涂前沿海区，另一股沿着刘五店海岸方向 流向琼头方向；落潮时，湾内落潮流沿上述三条路径往湾外流出， 在鳄鱼屿附近水域汇成一股落潮流流向湾外。进入西港的厦f - j # l , 港 潮波，涨潮时由厦鼓、嵩鼓水道进入西港，厦鼓水道和嵩鼓水道两 股涨潮流汇合后，沿西港主航道进入宝珠屿海域，一股潮流北上形 成以宝珠屿为中心的逆时针环流，另一股潮流向西沿马銮水道逐渐 地进入马銮湾；落潮时，滩地水流向深槽归聚，通过排头断面进入 火烧屿水道，后经厦鼓水道和嵩鼓水道向外港运动。 北面的潮流经过围头湾，进入大嶝南侧大片浅水区，在大嶝南 端一金门古宁头断面附近海区，在涨潮过程与南面涨潮流会合；落 潮过程从该海区分为南北二支落潮流，逆涨潮路径各朝南，东北向 外海流出。 ( 2 ) 本海域强流区位于领近外海的围头湾和厦门湾口区的深水 部分、潮流通道区( 厦门东侧水道，厦f - l # 、港) 以及九龙江河口湾 属次强流区。弱流区出现在内湾( 宝珠屿海区，同安湾) 、沿岸浅水 区、大嶝南侧海区，潮间带浅滩为漫滩水流，流速一般不大。整个 海域潮流呈现明显的往复流特征，符合一般半封闭型海湾的潮流性 质，高低潮时刻流速最小，转流时刻出现在高、低潮时刻附近，涨 落潮半潮位时刻( 涨急、落急) 流速最大。 经过上述许多测流点和潮位站的潮流、潮位验证，可以说本模 型模拟的潮流场基本上接近实际的潮位潮流特征，具有良好的重现 性。 图5 a8 0 0 m 网格厦门海域潮流场( 落急) 图5 b8 0 0 m 网格厦门海域潮流场( 低潮) 图5 c8 0 0 m 网格厦门海域潮流场( 涨急) 3 5 图5 d8 0 0 m 网格厦门海域潮流场( 高潮) 图6 a4 0 0 m 网格厦门海域潮流场( 落急) 图6 b4 0 0 m 网格厦门海域潮流场( 低潮) 图6 e4 0 0 m 网格厦门海域潮流场( 涨急) 图6 d4 0 0 m 网格厦门海域潮流场( 高潮) 图7 b5 7 m 网格厦门海域潮流场( 低潮) 图7 c5 7 m 网格厦门海域潮流场( 涨急) 图7 d5 7 m 网格厦门海域潮流场( 高潮) 4 2 边界条件对潮流场模拟精度的影响 4 2 1 两个分潮与十个分潮情况对比 图8 、图9 是分别选取m 2 、s 2 两个分潮与包括m 2 、s 2 等十个 主要分潮进行模拟得出的潮位和潮流曲线。可以看出，在大潮期间， 仅代入两个分潮的模拟结果与实测资料比较相近，在一个潮周内十 个分潮的潮位模拟结果好于两个分潮潮位模拟结果，特别体现在小 潮期间的潮差上。由此可以看出，在大潮期间，m 2 、s 2 分潮近似代 表了本海域的潮汐驱动力，而小潮期间m 2 、s 2 分潮对潮汐驱动力 的贡献相对弱些，计算结果偏差较大。事实上，在水动力数值模拟 早期，许多工程常常选用m 2 、s 2 分潮来模拟大潮期间的潮波运动， 就是这个道理。另外，本海域潮汐形态属于正规半日潮，但仍然存 在潮汐日不等现象，只是不特别明显。用十个分潮模拟的潮流曲线 上可以看出本海域的潮流日不等现象，即涨落潮历时不相等，高低 潮潮流差不相等。m 2 、s 2 两个分潮的模拟结果涨落潮历时几乎相等， 高低潮潮流差也几乎相等，潮流日不等现象不明显。 位( c m 4 0 0 3 0 0 2 0 0 1 0 0 0 1 0 0 - 2 0 0 - 3 0 0 - 4 0 0 高集海堤以东 02 44 8 7 29 6 1 2 0 1 4 4 1 6 81 9 22 1 62 4 02 6 42 8 83 1 2 5 r 1 4 b 点 位c m ) 一1 - 2 - 3 4 实测值计算值澳头 时匍m 3 3 6 5 月2 7 日2 3 o2 4 5 月1 4 r 直 位c m ) 4 0 0 3 0 0 2 0 0 1 0 0 0 1 0 0 2 0 0 一3 0 0 - 4 0 0 4 87 29 6 图8 1 1 2 01 4 41 6 81 9 22 1 62 4 02 6 42 8 83 1 2 两个分潮模拟结果与实测潮位曲线 高集海堤以东 时同m 3 3 6 5 月2 7 日2 3 点 02 44 87 29 61 2 0 1 4 41 6 81 9 22 1