非结构网格的LUSGS隐式算法在浅水方程组中的应用（可编辑）

非结构网格的LUSGS隐式算法在浅水方程组中的应用（可编辑） 非结构网格的LUSGS隐式算法在浅水方程组中的应用 河海大学 硕士学位 非结构网格的LU-SGS隐式算法在浅水方程组中的应用 姓名:杨彬 申请学位级别:硕士 专业:流体力学 指导教师:汪德 20080301 河海大学硕士学位论文 非结构网格的LU-SGS隐式算法在浅水方程组中的应用 摘要 鉴于隐式格式的高效率和较好的稳定性，以及近年来非结构网格技术的发 展，本文在前人的研究基础上，借鉴空气动力学领域中的相关算法，将用于结构 Gauss SymmetricSeidel 隐式算法应用于非结构 网格的LU-SGS Lower-upper 网格，并考唐浅水数值模拟的特点，推导并建立了在非结构网格上求解浅水方程 的LU―SGS隐式算法。本文引入斜底模型的概念，即将地形高程布置在网格节点， 而守恒变量定义在网格中心，并在本隐式算法的推导过程中对底坡项进行了修 正。本文在斜底模型概念的基础上，引入淹没节点法来处理动边界问。通过对 本文所选择的具有解析解或实测验证数据的算例的模拟，以及与显式Roe格式模 拟结果的比较。明本文所建立的LU―SGS隐式算法模型不但具有较好的稳定性 和快速的收敛效果，而且也能够有效的处理地形的起伏，比较好的捕捉激波。最 后对长江口实际潮位的模拟验证表明本算法模型在非结构网格上求解实际浅水 问题是切实可行的。 关键词:非结构网格，浅水方程，LU-SGS算法，复杂地形，移 动边界。 河海大学硕士学位论文 LU-SGSSchemeforshallow-water on equations Unstructured Grids Abstract Themethod缸two-dimensionshallow-water used equation compumfion’which LU―SGSschemebasedontmstructured in discussedthis implicit grids(was paper(As the willaffect and the ofthemathematical topograpbic订estmentstabilRyrationality model，a modelwas tothetreatmentof sloping applied topography(The used todealwiththe submerged-nodemethod SNM was problems concerning WerCcalculatedandtheresultswere boundary(Two moving examples compared、析tll the Roescheme(ItWasshowedthatthe and of explicit efficiency this convergence methodWer0 The resultof withthe goocL computed YangtzeEsmaryagreed observationwell(TheresultsindicatedthatthenCWschemewasaccurateand 踊cient( water equation,LU-SGS Keywords:unstructtwed鲥ds，shallow scheme,complex topography,movingboundary( (2( 学位论文独创性声明: 本人所呈交的学位论文是我个人在导师指导下进行的研究工作 及取得的研究成果。尽我所知，除了文中特别加以标注和致谢的地方 外，论文中不包含其他人已经发表或撰写过的研究成果。与我一同工 作的同事对本研究所做的任何贡献均已在论文中作了明确的说明并 表示了谢意。如不实，本人负全部责任。 论文作者 签名 : 年 月 日 学位论文使用授权说明: 河海大学、中国科学技术信息研究所、国家图书馆、中国学术期 刊 光盘版 电子杂志社有权保留本人所送交学位论文的复印件或电 子文档，可以采用影印、缩印或其他复制手段保存论文。本人电子文 档的和纸质论文的内容相一致。除在保密期内的保密论文外，允 许论文被查阅和借阅。论文全部或部分内容的公布I包括刊登 授权河 海大学研究生院办理。 论文作者 签名 : 年 月 日 河海大学硕士学位论文 第一章绪论 1(1概述 几十年来，国内外很多学者对浅水动力学的数值模拟采用各种计算方法作了 大量的研究n1。由于浅水方程的非线性，一般情况下无法获得解析解，主要通过 数值离散方法求解，目前，浅水方程的数值模拟技术已经比较成熟。许多工程中 的流动问题常常可以用浅水流动模型来描述，如近海、湖泊、水库及河口中的流 动。由于浅水方程在数学上属双曲方程组，在一定的条件下会产生如水跃、溃坝 波、涌浪等强间断解，给方程的数值求解带来一些困难，普通的有限差分格式不 能有效处理间断解，因而大大限制了所编制程序的应用范围12】。20世纪六七十 Volume 年代F(H(Harlow等01提出了有限体积法 FiniteMethods，简记为Fv岫 的基本思路，并应用于空气动力学领域的研究中(80年代至今，非结构网格技术 的发展和对数值通量算法的改进，使得FVM在方法设计和构造、处理激波等复杂流 体动力学，以及方法精度和收敛性的理论等方面均有了实质性的进展【4l，并结合 其他一些数值方法，如有限元法，Runge-Kutta方法，还有新发展起来的MUscL哪， RoeIs]，TVDLT]，EN0【g】等方法，产生了一系列新颖有效的高性能算法。近年来，很 多学者‘唧利用齐次浅水方程和可压缩Euler方程在数学形式上的相似性，借鉴空 气动力学中的这些高性能算法，研究具有守恒性、迎风性、单调性和大梯度分辨 率的浅水方程数值解法m”】，并成功应用于实际工程当中，取得很大的成功。 1(2问题的提出 到目前为止，浅水动力学的数值模拟方法已经相当的成熟，优秀的数值计算 格式层出不穷。在实际计算不规则天然地形上的浅水流动时，非 结构网格尺寸常 相差较大，我们在利用显式格式时，时间步长往往受到所有网格中稳定性限制中 最小者的控制，极大地降低了处理效率，特别是在计算历时较长的数值模拟中， 计算所需要的实际时间往往让人难以接受，从而使数值模拟受到限制。但是在非 结构网格上采用隐式算法比较困难，因为很多快速收敛的高效算法都是基于有结 构网格，比较典型的隐式算法如交替方向隐格式 ADI 法及近似因子分解法 四 河海大学硕士学位论文 等主要适用于有结构的矩形或曲线网格【9】。非结构网格舍去了网格节点的结构性 限制，使很多结构网格计算中的高效算法无法适用。目前在非结构网格的计算中， 已经有一些隐式方法用来加速收敛?„。在非结构网格上解线性方程，应用最多 的是迭代法和近似因式分解法。比较成熟的迭代方法主要有采用前处理器的 Gauss-Sedial及Krylovsubspace方法IJ删，但这些方法往往需要复杂的矩阵计 算，大量的存储空间，对于大型的数值计算问题，可能无法完成计算。 Symmetric Gauss Seidel 隐式算法最近有学者成功【26。1的将其应用于非结构网格求解齐次 欧拉方程，LU-SGS方法采用对隐式算子近似因式分解来完全消除对矩阵的存储， 使构造的L、U算子具有最大程度的对角占优，解的收敛速度和程度有很大的提 高。因此借鉴齐次欧拉方程的Lu―SGS隐式解法，建立适合在非结构网格上求解 浅水方程的LU-SGS隐式方法，不仅可以有效的加快收敛速度，而且可以增加格 式的稳定性，将有利于恒定流或解的时间变率小、流动过程历时长等一些实际浅 水问题的数值模拟。 由于齐次浅水方程与欧拉方程在数学形式上的相似性，因此数值通量的求解 方法是类似的。但是鉴于浅水方程与欧拉方程的不同，在非结构 网格上实现浅水 方程的LU-SGS隐式格式还有许多问题需要解决。 1(3浅水数值模型的研究与发展 1(3(1数学模型的主要数值计算方法 目前较为流行的数值计算离散方法主要有:有限差分法 FDM 、有限单元法 等。这几种方法都有优缺点以及各自不同的适用范围( 1 有限差分法 有限差分方法 FDM 是计算机数值模拟最早采用的方法，至今仍被广泛运 用。该方法将求解域划分为差分网格，用有限个网格节点代替连续的求解域。有 限差分法以Taylor级数展开等方法，把控制方程中的导数用网 格节点上的函数 值的差商代替进行离散，从而建立以网格节点上的值为未知数的代数方程组。该 方法是一种直接将微分问题变为代数问题的近似数值解法，数学概念直观，表达 2 河海大学硕士学位论文 简单，是发展较早且比较成熟的数值方法，其特点是采用矩形网格概化计算域， 各计算变量布置在网格节点上。对于有限差分格式，从格式的精度来划分，有一 阶格式、二阶格式和高阶格式。从差分的空间形式来考虑，可分为中心格式和逆 风格式。考虑时间因子的影响，差分格式还可以分为显格式、隐 格式、显隐交替 格式等。目前常见的差分格式，主要是上述几种形式的组合，不同的组合构成不 同的差分格式。差分方法主要适用于结构网格，网格的步长一般根据实际地形的 情况和柯朗稳定条件来决定。在三维计算中，往往对空间的不同方向分别采用不 同的差分格式，其中交替方向隐式差分法 ADI 是目前较为常用的一种方法。 有限差分法的优点是求解直观迅速，数学概念清晰，编程及基础资料的整理比较 简单，占用内存少，计算的收敛性和稳定性理论比较成熟;其最大的缺点是对复 杂的地形适应性差，网格布置不灵活，难以概化不规则的边界。为了克服其局限 性，许多学者都致力于不规则边界处理问题的研究，如美国的 ThompsonJ(F(等 FittedCoordinate 人提出的边界拟合坐标系 Boundary System 方法【2明，L( 谢利曼所提出的任意网格差分法唧，华祖林口ol等也对边界拟 合中的调节因子作过 有益的探索。利用这些方法在原则上可以把任意复杂的几何边界 变成规则的几何 边界求解，但计算域的规则化是以控制方程的复杂化为代价的， 而且当边界存在 尖角时，会出现局部奇异现象，使计算不能收敛。 2 有限单元法 有限单元法 FEM 是一种常用的数值分析方法，又称“有限元素法”。有限 单元法最初是在20世纪50年代作为处理固体力学问题的方法出现的，从20世 纪60年代后期开始，进一步利用加权余量法，主要是Galerkin法，来确定单元 特性和建立有限单元求解方程，使之应用于已知问题的微分方程和边界条件、但 变分的泛函尚未找到或者根本不存在的情况，进一步扩大了有限元法的应用领 域，从70年代起开始用于计算水力学中。有限单元法也是一种近似解法，是一 种用积分求解微分方程的方法。有限单元法以变分原理和剖分插值为基础，将连 续的结构离散为有限个单元的组合体，单元与单元之间仅靠节点 相连。有限元法 是应用局部的近似解来建立整个定义域的解的一种方法。先把注意力集中在单个 单元上，进行上述所谓的单元分析。基本前提是每一单元要尽可能小，以致其边 界值在整个边界上的变化也是小的。这样，边界条件就能取某一在节点间插值的 河海大学硕士学位论文 光滑函数来近似，在单元内也容易建立简单的近似解。这使有限元法具有明显的 优越性，由于其采用分块近似，只需对一个单元选择一个近似位移函数，且不必 考虑位移边界条件，只须考虑单元之间位移的连续性，因此对于具有复杂几何边 界的情况具有较高的精度，很好的克服了最初的有限差分法的规则网格划分不适 应复杂边界描述和不能满足计算域内局部加密的要求等缺点(有限单元法的基本 原理就是分单元对解逼近，使微分方程空间积分的加权残差极小化。有限单元法 的优点主要是划分网格灵活，拟合复杂边界和地形容易;其缺点是占用内存大， 程序编制比较复杂，对基础资料的整理比较繁琐，一定程度上限制了有限单元法 的推广应用。另外随着有限差分法在任意三角形和四边形结构网格的推广，有限 单元法的优势已不明显。 3 有限体积法 有限体积法 FvM 舢又称为控制体积法，其基本思路是:将计算区域划分 为一系列不重复的控制体积，并使每个网格点周围有一个控制体积，将待解的微 分方程对每,个控制体积积分，便得出一组离散方程，其中的未知数是网格点上 的因变量中的数值。为了求出控制体积的积分，必须假定中值在网格点之间的变 化规律，即设定中值的分段的分布剖面。从积分区域的选取方法看来，有限体积 法属于加权余量法中的子区域法;从未知解的近似方法来看，有限体积法属于采 用局部近似的离散方法。简言之，子区域法加离散，就是有限体积法的基本方法。 有限体积法离散方程的物理意义，就是因变量巾在有限大小的控 制体积中的守恒 原理，如同微分方程表示因变量在无限小的控制体积中的守恒原理一样。就离 散方法而言，有限体积法可视为有限单元法和有限差分法的中问物。有限单元法 必须假定值在网格点之问的变化规律 即插值函数 。并将其作为近似解。有限 差分法只考虑网格点上中的数值而不考虑中的数值在网格点之问如何变化。有限 体积法只寻求中的节点值，与有限差分法相类似;但有限体积法在寻求控制体积 的积分时，必须假定m值在网格点之间的分布，与有限单元法相类似。在有限体 积法中，插值函数只用于计算控制体积的积分(得出离散方程后，如果需要的话 可以对微分方程中的不同的项采取不同的插值函数。有限体积法的优点是离散 方程的解表示变量 如质量、动量等 在任意一组控制体积上的总体积分是守恒 的，对整个计算区域自然也得到守恒;对于有限差分法，仅当网格极其细密时， 4 河海大学硕士学位论文 离散方程才满足积分守恒:而有限体积法不仅能在精细的计算网格上获得较好的 结果，而且在粗网格的情况下，也显示出准确的积分守恒。在网格布置时。往往 采用交错网格，将不同的物理量布置在不同的节点上，此外，由于三维正交网格 的生成存在诸多困难，非正交网格下的有限体积法以及无结构网格下的控制体积 分也得到一定的应用。因此。有限体积法在流体力学的数值计算中得到了广泛的 应用，本论文就是利用有限体积法对方程进行离散的。 4 特征线法 特征线法 MOC 【321的最初思路是在x-t平面上绘制特征线，在其交点上确 定因变量来依次求解，后来在特征线理论的基础上发展了特征线法。该方法把时 问离散和空间离散一起处理，其优点是能反映问题中信息沿特征传播的性质，算 法符合水流运动的物理机制，稳定性好，计算精度高。由于该方法是沿时问推进 求解，故较适合于双曲型和抛物型问题，对于求解周期短、变化急剧的问题 如: 涌潮 比较适宜。推广到二维问题，由于二维问题中对应于一维问题的特征线是 两组特征曲面，表现为一个特征锥面，目前一般是对特征锥面选用几条母线，沿 对应的特征关系式积分来近似求解特征量。因特征线法求解格式复杂，尤其是对 高维问题更为繁琐，因此目前很少直接用于数值计算。但是，特征线法的原理仍 是很重要的，经常作为了解其它数值方法的基础。 《5 边界单元法 边界元法 BEM 是继有限元法之后发展起来的一种新数值方法，最初被称 为边界积分方程法，在流体力学中又称为有限奇点法或者有限基本解法。边界单 元法的实质是数学物理方程的Green函数法。通过格林公式给出解的积分表达 式，然后利用定解条件建立边界积分方程，所得到的边界积分方程通常不能解析 求解，因此可以利用有限元离散，将其转化为边界节点上未知量表示的代数方程 组，求解该代数方程组，即可求得边界节点所有的未知值。由此可见，边界单元 法可以看作格林函数法与有限单元法结合的产物。但与有限元法在连续体域内划 分单元的基本思想不同，边界单元法是在定义域的边界上划分单元，用满足控制 方程的函数去逼近边界条件。边界单元法可以将三维积分转化成二维面积分，再 通过格林公式将面积分变换成线积分，从而将空间三维问题简化为一维问题处 理，有效的减少了计算工作量和所需内存，对于三维水流计算中自由表面的处理 5 河海大学硕士学位论文 较为简单。边界元法与有限元相比具有单元未知数少，数据准备简单，精度高等 优点(但用边界元法解非线性问题时，遇到同非线性项相对应的区域积分，这种 积分在奇异点附近有强烈的奇异性，使求解遇到困难。边界单元法是计算椭圆问 题的有效方法，但由于需要控制方程的基本解，所以对于复杂的问题如完整的 N-S方程尚未得到广泛应用。 6 分步法 分步法 F跚 又称破开算子法，其基本思想是将待解的偏微分方程分解为 相对简单的若干部分，对于每一时间步长，分别求解偏微分方程的不同部分，逐 步达到偏微分方程的解。对于不同部分采用不同的离散方法，建立不同的离散方 程，用数值方法分别求解各部分方程，从而得出偏微分方程的解【，，1(分步法分为 空问概念上的分布 将多维问题分解成一系列一维或低维问题来求解 、物理概 念上的分布 将原始的包含若干个物理性质不同的复合过程分解成若干性质单一 的物理过程来求解 和解析上的分布 将原始解析方程进行分解 三种类型四。 分步法是求解二维水流的一种简单易行的数值计算方法，在处理上具有很大的灵 活性。理论上只要在分步之后，寻找出一个尽量正确、合理的数值离散方法，就 能取得良好的计算结果，分步法的主要缺点在于处理边界条件上。在此基础上， 湖泊水动力学数值模拟借鉴许多大气或海洋学中的先进方法，发展了不规则网格 的有限差分数值模拟技术、有限分析法、边界分析法、嵌套网格技术等，使数值 模拟方法能更加贴近实际，渗透到了湖泊学研究的各个领域。 1(3(2计算区域网格剖分 网格剖分与算法、边界处理等对数值模拟实际流动同样重要。提高精度和分 辨细部常常要依靠网格的合理布置和适当加密，随着计算机技术的发展，计算量 越来越不成为限制。网格建立的基本要求是:符合流动特点，易于建立，建立网 格要比较光滑、规则，满足精度和计算稳定性要求，便于组成节约，高效的数据 结构，必要时可随时根据解的梯度作适应性调整。网格分为有结构网格和无结构 网格p例。 1 有结构网格 有结构网格的网格主要包括矩形网格及贴体曲线网格，后者可以通过坐标变 换映射成矩形网格(有结构网格最早使用直角坐标系中的矩形网格，容易确定单 6 河海大学硕士学位论文 元格之间的邻接关系，有利于用差分逼近导数。但是每个格子的边长与相邻格子 边长比要满足一定的限制，以保证精度(在计算边界层时可以使用拉长的格子， 但是需要进行坐标变换处理。矩形网格便于组织数据结构，程序设计简单，适于 各种算法，处理效率较高(主要缺点是把计算区域概化成锯齿形边界，陆地边界 附近出现虚假的曲折水流，难以处理边界上的奇点，边界附近的误差较大。70 年代使用的贴体曲线网格，往往只能对几何形状简单的计算区域建立网格，形状 复杂时可以用若干割线将区域剖分为简单的子域后分别建网，但是需要解决分块 之间解的协调。有结构网格的主要优点是:生成技术比较成熟，网格的生成和安 捧是有序的和按照一定结构的，在平面情形下节点编号及节点变量可用二维数组 来表示，适合于求解区域有几何规律和求模型比较均衡的问题，正交网格比较适 合于差分计算的要求。 2 无结构网格 无结构网格是用任意三角形或四边形构成的不规则网格，最早是有限元方法 采用这类网格。目前多用凸四边形，是因为:用于显格式计算时时间步长比三角 形大;节点数相同时三角形网格的格子和边的数目为四边形网格的两倍或多倍， 计算量大;拉长的三角形上一阶方法的精度和稳定性较差;四边形网格的解较好， 二阶粘性项较容易处理。两者可以混合使用，以四边形为主，以三角形为补充， 后者用在局部地形巨变、租细网格过渡及曲折边界处。无结构网格的优点是:与 边界及水下地形拟合比较好，利于边界条件的实现;便于控制网格密度，易于修 改和适应性调整;建网比曲线网格容易，大型三角形网格可用程 序自动生成。缺 点是:格子排列不规则，需建立适当的数据结构以检索格子间的邻接关系，占用 内存多;间距寻址时，解的精度较低;隐格式求解效率低，粘性项处理比较麻烦， 数值解后处理工作量大。由于有限元求解效率低，无结构网格长期未能广泛应用， 直到80年代与有限体积法结合，经过多年的研究，才得到广泛的应用。 无结构网格是当前数值计算中一个受人关注的重要发展方向，其网格的生成 和安捧可以是无序，可以根据需要进行加密或减疏，易于构造自适应网格，这是 有结构网格难以做到的，无结构网格能适应复杂的求解区域，其与边界及水下地 形拟合较好，利于边界条件的实现，这一点也比有结构网格优越，原则上讲，无 结构网格可以具有任意结构，但为了实现和计算上的方便，多采用同一类型的网 7 河海大学硕士学位论文 格，对二维问题多为三角形，三维问题多为四面体。 1(3(3实际地形的处理 对二维浅水方程而言，复杂地形的模拟精度对计算结果有着重要的影响，这 是因为在推导二维浅水方程时隐含了地形变化的影响，忽略这种因素将影响到数 值解的合理性。若不考虑河床的起伏和摩阻影响，浅水方程成为 齐次方程，在数 学形式上类似于气体动力学中的Euler方程，源项的加入增加了数值求解的难 度。对其中的摩阻项既可采用显格式直接求解，也可以采用显隐格式求解[40(4“， 无论采用哪种方法对计算结果都影响不大，而其中的底坡项，即地形变化影响， 却是模型中必须考虑的问题。因为河床地形起伏不平，如忽略这种起伏变化，产 生的最明显的问题就是计算结果的不和谐。 为了解决这一难题，很多学者对此做了大量研究，并取得了丰硕的成果。最 早的地形处理方法采用中心差分格式，这在计算中容易产生伪流速。胡四一和赵 橡华[42删等学者在计算单元内采用平底假设，河道地形概化为阶梯状，从而省略 了底坡项，将单元界面处流动转换为阶梯流动，并针对特殊单元采用堰流公式加 以校正，这种处理在物理上把明渠流概化成类似多级跌水形式，在数学上则是用 阶梯状平台去逼近实际地形，这种方法比较适合地形变化平缓的水域水流模拟。 Gradient 周建国M镯提出“水面梯度法” Surface Method，SGii ，该方法用 水位重构代替传统的水深重构，对地形变化剧烈处 如台阶 也能取得较好的结 果，这种方法比较适合应用于结构网格下水流计算。许为厚和潘存鸿l删提出了“水 Level 面方程法” Water Formtion，1rLF ，采用分步法求解方程。谭维炎【9J 提出斜底模型的概念，即单元内的地形是倾斜的，假定河底高程在控制体内线性 变化，用斜面代替实际地面，相邻单元地面高程一般连续，但坡度不同。与平底 模型相比，相当满意地逼近水下地形，使得模拟出的流场更合理、精确。计算时 只用到各边中点和单元中心的高程，并不涉及单元内的地面形状。将底坡项作为 对流项中的静水压力项的修正。即取消底坡项，并将实际静水压力代之以计算到 该单元的公共基准面的换算静水压。 鉴于复杂地形的处理对浅水模型建立的重要性，本文有必要在吸取前人的经 验成果之上，建立一种适合本隐式算法模型的新的地形处理方法。 河海大学硕士学位论文 1(3(4动边界问题 对河口海岸地区的数值模拟中，伴随着水位的变化，水域的边界将发生改变， 这就是动边界问题。动边界处理方法有很多，如开挖法、冻结法、切削法、干湿 法、阻塞函数法和窄缝法等，其中目前常用的方法主要有冻结法、切削法、干湿 法、窄缝法和开挖法[9(47AS]( 1 冻结法 冻结法根据网格单元中心处水深判断网格单元是否露出水面，对露滩单元取 糙率系数为一接近无穷大的正数，使单元四周的流速都为趋于零的微小量，使该 单元水位在计算时被冻结不变。该方法适用于宽浅、坡度较平缓的岸滩，而对于 槽滩相间的海岸河口水域会由于水量和动量的过分冻结而失真。 2 切削法 切削法与冻结法一样判断单元是否露出水面，但并不冻结水位，而是引入一 个富裕水深来保证计算过程的完整和稳定，相当于将原始地形切割降低，而一旦 判断实际水深大于富裕水深时，即恢复原始地形。该方法解决了 冻结法存在的冻 结附近水位梯度过大的问题，但对差分方法要求比较严格，需要用滤波形式( 3 干湿法 干湿法又称水位判别法，是根据计算点的水深及其相邻半网格距处的水深和 水位值判别该计算点的干湿，而令露滩点流速沿此方向为零。这种方法具有清晰 的物理概念、简单的实现过程和良好的计算效果，应用比较广泛。但是临界水深 的判别标准比较难确定，临界水深过大则动边界模拟精度较差，临界水深过小， 则可能导致计算失稳。 4 窄缝法 窄缝法假想在岸滩的每个网格上存在一条很窄的缝隙，它的深度和岸滩前的 水深一致，其宽度是水位相对于单元中一13,高程距离的函数。根据水量平衡原理将 窄缝内的水量平铺到岸滩上，成为化引水深，相当于把岸滩前的水域延伸到岸滩 内，从而可以把计算网格点布置在岸滩上。这种方法只适用于岸边滩的露滩问题， 而且连续方程的形式需要修正。 4 开挖法 开挖法在程序设计上比较简单，在漫滩占计算区域较小的情况下还可以得到 9 河海大学硬士学位论文 较满意的计算精度。而对水位变化小，坡度较缓的地形，则比较容易失去真。 动边界问题也是浅水数值模拟的一个关键性难题，为了能够准确的模拟水域 边界的变化以及涨潮漫滩，落潮归槽的潮流运动等一些实际问题，结合本模型算 法引入适当的动边界处理技术是十分必要的。 1(4网格重新编号 所谓网格重新编号就是将已经生成好的网格单元编号进行重新排序，对于 uJ(sGS隐式算法，排序的结果要求每个网格单元的邻居单元编号既要有比其大 的(也要有比其小的。网格重新编号一般既可以通过编程实现 [49(50l，也可以用已 有软件来实现。 目前重新编号的方法主要针对于网格节点，其目的是缩短有限元方法的总体 刚度矩阵的带宽，从而提高线性方程组的求解效率。主要方法有波前法和矩形法。 波前法的基本思想是先寻找一个边角节点，使该节点的编号为o，新编号节点不断 地沿着已编号节点的前沿向外扩展，直到所有的节点都已编号为止;矩形法是模 拟形状简单的矩形区域按行 列 编号的一种方法，它将节点的编号转化为按节点 空间位置进行排序。 通过对比分析，网格单元重新编号采用类似的方法也可以达 到本文算法的要 求。目前有很多专业的网格生成和前处理工具，如GAMBIT、TGrid、GeoUesh、 preBFC、ICEMCFD等。一般生成的非结构网格在单元编号上具有无序性，但通过 相应软件的优化处理后，还是能达到提高网格质量的目的。本文针对LU―SGS隐式 算法的要求，提出了一种简单的网格编号重新排序的方法。 1(5本文主要研究工作 本文在前人的研究基础之上，借鉴气动力学中相关隐式算法，建立了在非结 构网格上求解浅水方程的叫(SOS隐式算法模型。并考虑浅水方程与齐次欧拉方 程的不同，在模型的建立过程中对地形和动边界问题进行了处 理。并利用所建立 的隐式算法模型对几个标准算例和实际工程资料进行模拟验证。本文主要在以下 几个方面做了尝试性研究(并取得了一定的成果。 1 根据齐次浅水方程与欧拉方程的相似性，对于数值通量的求解采用基于 0 塑墨查兰塑主兰垒堡奎 黎曼解的二阶Roe格式。 2 参照欧拉方程的LU-SGS隐式算法的建立，推导并建立求解浅水方程组 的LU(SGS隐式算法模型。 3 鉴于浅水数值模拟的特点，在斜底模型概念的基础上推 导本文的地形处 理方法。 4 结合本文采用的地形处理方法，引入淹没节点法解决移动边界问题( 通过对计算网格的重新编号后，利用本文所建立的LU-SGS隐式算法模型对 典型算例进行模拟验证，并与同条件下建立的显式Roe格式模型的模拟结果进行 比较，最后对长江口的实际潮位进行验证，结果表明本文所建立的LU-SGS隐式 算法可以有效的增加格式的稳定性，加快收敛速度，节省模拟计算时间，并且可 以用于实际浅水问题的数值模拟。 河海大学硕士学位论文 第二章 数值模型的建立 2(1控制方程 本模型采用二维浅水方程组，二维浅水方程组的连续方程和动量方程可以表 示如下洲: I丝+型+型:0 |毋盘 勿―― 2-1 其中:h一水深，乙一底高程，P一水的密度，v一运动粘性系数，g一重力加 速度，雄一曼宁糙率系数，“、v一分别为x和Y方向的垂向平均流速，名和瓦分 别表示科氏力项。 为了说明问题的需要，在本模型的控制方程中忽略了扩散项 和科氏力项。方 程组进而写成如下形式: 丝+旦盟+塑:o af 知05, 掣+扣+却+挈叫鼍+学协z， 掣+掣+耖+却 叫苦+吲浮 hu ，F E U i+G U j， hu2+ 1曲2 叫弘卟 二 hm, 12 河海大学硕士学位论文 加 O ^w s ? l―g^ ,，将方程组用矩阵形式来表示: G c， 一, 舻+三‖ 一曲 岛一, 2―3 丝+V(F ?:o_5_u+-扭--; -?--+_OG U :s ?、7 、’ 出 西 四 al， 其中:屯、,表衣和y方向的底坡项，, 鲁，岛 鲁??,、,表 示埔y方撇摩阻觚,:型h4噩13，,:呜霉。 2(2空间离散 2。2(1有限体积法离散原理 有限体积法的基本原理是在被离散化了的计算区域上，计算出通过每个控制 体边界沿法向输,X(或输出的流量和动量后，对每个控制体分别进行水量和动量的 平衡计算，最终得到计算时段末各控制体的平均水深和流速。 y X 图2(1单元有限体积离散示意图 本文对任意计算区域采用三角形网格，将各个变量定义在网格中心，即采用 网格中心格式 Cell(Centered，co 研，控制体即是网格单元。对方程在控制体积 塑堡查兰堡主兰篁!垒兰 即:f【警dQ+见 差+筹洒 儿姗 z―s 其中:Q为控制体平面区域 面积为彳 ，dQ表示面积分的微元?? 对E和G的面积分利用格林公式，将对控制体的面积分化为对其周边界的线 其中:s为控制体周边 逆时针方向 ，ds为线积分微元( 方程左边第二项在数值上等于被积函数在控制体各边上的法向值与该边长 度的乘积，并假定水力要素在各控制体内均匀分布，方程可以写成如下离散形式: A百OU+善暑?,， zs? 2-7 其中:A为控制体的面积，册为控制体的边数，J为控制体单元的第(，条边， 与为控制体单元第(，条边的边长?? fhucos目+hvsin口1 J 2―8 【huvcos0+ ^'，2+gh2,2 sin0J 利用欧拉方程的旋转不变性嘟1，可以将二维问题转换为一 系列一维问题求 解。州定义旋转矩阵H曰 和其逆矩阵r。’ 力: 叫i 0小sinO叫撼0 0] 由关系r 功E 【_r F “口 U 晌，所以‘ u r’ 力F 回。 旧c眇眩0菲]I[兰揣] 14 塑堂查兰堡主塑笙壅 疗 [兰箩]，Ft西 [兰毫，+朗2，z]，所以方程可以进一步化为: 4511:I二盟+羔1 1一 占妒 西岛 彳s ? 2_9 将非线性方程组掣+爿掣:o局部线性化成为警+j罢:o，其中j为由相邻 讲 斑 讲 啦 单元边两边局部变量叽和,求出的Roe平均状态，j j 叽，, (在 j冻结为 河海大学硕士学位论文 j 雩 陆2鲫最 隧叫心隧叫 其中:疗 及厅与前面的定义相同。 j矩阵的特征值分别为: A u+c，五 【，一c，五 u 对应的右特征向量为; ^ 1，厅+a乃7，乃 1，0一己乃7，乃 o'o’1 7 以^，乃，乃为基向量， 玩一玩 的特征分解为: ? 一U 一 lI 一, 一叽 坼圪 ，?? 舯嘶 业d等盟出，吃 限枷一q 龟 伪。,一,圪 一矿 致一k 所以法向数值通量可以写成如下形式: ， 西 三[最+丘卜i1刍3,I以I靠 2-11 但是R潞式不能正确识别含特征值零点的膨胀波，将其误认为问断。因此 可能将膨胀波作为膨胀间断来处理，从而违反了熵条件例。为了解决这个问题， 需要进行熵修正。处理方法是: 1 当五 o时，在上、下游分别增、减同样的通 量值。 2 当I五l 占时，取陋|_m双 o，五一五。，五^(-,b。对压缩间断，其中零 起作用而不作校正。 3 当|互I 占时，取l盂lI学，这样校正后的五连续可 微。 通过熵修正的Roe格式数值通量如下: 16 河海大学硕士学位论文 厶 二^Il 鼎以， 瞥拳一嘶“以州国 变量进行线性重构，以使空间离散达到二阶精度格式。设蚝和?分别为相邻控 制体单元公共边两边需要重新构造的守恒变量值，叱和?的重构是有限体积法 法，为了介绍方法，将控制体单元巧的中心的守恒变量设为,，设与控制体单元 有公共边乇的控制元圪。的中心守恒变量为坼+1( 度的数值格式。线性逼近引进参考量,(、Ul+:，所以这种方法以控制元为单一 制体的中心，设,。、“』+2分别为巧、K。的与公共边乇相对的角点上的守恒变量， ? ,。一兰烈?, ‘2-13 其中，一 2?一坼，? 2q一?，‘2哆+;7,’?5? 7砌(+;’ 17 河海大学硕士学位论文 伊是限制器函数(这是二阶精度方法，计算量中等，计算结果 比较好( Ii 图2(2MUSCLTj法示意图 2 Upwind 迎风 型格式1划 对于给定的中心点为口的控制体矿，对任意伍y Ey，将;f 似力在O点的作 Taylor 襄开，舍去二阶以上的高阶项，则有: u x,y ,+,??r 2一15 这里的，为点。到点 五， 的向量，,为“在。点的梯度，其近似值可以通过 以下公式计算: 其中积分路径翩是点。周围已知其越值的连线组成，』为动所围区域的有 向面积。分别在巧，,(上运用前面的方法进行就算，就可以得出叱，???为了确 保重构的单调性，实际计算中需要在公式中 2―15 加入限制器矿，即: u x，y Uo+伊 ， -,??r 2一17 这个是二阶精度方法，而且是标准的中心型有限体积法，计算结果要比 MUSCL格式要好，但是计算量比较大。 3 ENO或WENO重构05硎 在均匀矩形网格上可以构造出高阶逼近的，同时式中线积分的离散采用相应的高 翌塑茎兰塑主兰篁笙塞 阶格式，能够得到很好的计算结果，但计算量非常大，在非结构网格上的计算则 更加困难。 采用不同形式的数值通量函数和不同的吼，us重构方法，可以得到各种高阶 精度数值遁量的求解方法。根据各种因素的考虑，本文采用MUSCL方法来实现 二阶精度数值通量的求解。 2(2(4关于限制器 限制器 即限制函数 是对格式进行数值耗散的自适应调节，也就是能够自适 应地进行格式数值效应 耗散、色散 的调节，保持格式的单调性，以达到真正对 解的间断和大梯度的地方给出尖锐的图像。如果未知数的离散近似解有很大的梯 度或者振幅趋势，它可以适时的加强数值耗散;反之，在比较平缓时，可以降低 数值耗散，甚至局部适应地引入人工压缩，来达到数值解的高分辨率。妒 ， 中r 是前后相邻空间格点函数数值变差的比值四。通常，需要在上述二阶精度格式中 引入限制器。 在此之前已经有不少人提出了多种限制器，本文介绍几种建 立二阶格式【?‘I 常用的限制器( 1 Van Leer’monotionic 1974年 州带 藩:弗 cz一四 氟一 O，r出n 2，2r，0(5 1+r 2一i9 fo，当，? 烈， minmod[0，rain0，， 】 r，当10 r0 l 2-20 11，当ls， 4 Roe’superbee 1985年 河海大学硕士学位论文 o’ 当rs0 玩当o r?i1 烈， [O,min 2r，1 ，n血咿，2 】 Q。2D 1，当丢 ， l ，(当1 r?2 2’ 当2 r 5 Chakravarthy-Osher 1983年 fo， 当rso 识， [O,rain r,‖ 】 ，，当o rs妒 2―22 【,当妒 ， 其中P取l 缈?2 许多理论分析和数值实验已经证明TSupe:r'oee限制器在数值耗散效应上是 最弱的，VatLeer次之，Minmod而是最强的。但是Superbee数值耗散性是最强的， Minmod数值色散性最弱。因此，通常s叩erb?给出一个较陡的剖面图。Minmod 产生一个相对平坦的剖面图，VatLeer的居中[62(63】。使用时我们可以根据不同的 需要来选择他们。 一般来说，对于限制器的要求为: f o'当，?o 烈r o'当， o 2-23 【1，当， 0 为了保证格式保持单调性的要求为: 10茎烈， ?2r o?M， ’玛?2 ?24’ L ， 同时为了对向前、向后梯度等处理，常常要求“r 具有对称 性，即满足: 里盟:烈三 2―25 恿之，当烈， 满足式Q-23 ，且， O时，它的图形落在图 中的阴影区域内， 就可以保证格式是二阶级TVD格式。在允许区域内 即图中的阴影区域 选取不 同的非线性限制函数，就可以建立不同性能的格式。 河海大学硕士学位论文 图2(3限制器p ， 的要求 对于非线性守恒方程常采用双参数“，’，，+ 的限制器‘删: 咿，r+’畦刊么，匀’ 2-26 Au Au I 3 舯么2《’匀2毒’这里的卅一’，+燠种心对称的形式。双参 数限制器可使自适应耗散效果调节起来更灵活，并带有开关的作用，具有复杂灵 敏的形式，例如: 2-27 Yee’minmod型:烈，一，，+ rainmod 1，，一 +minrood 1，r+ 一1 2―28 2(2(5本文二阶精度格式建立 Van 是基于融锄蛐问题的近似解，对于捕捉激波非常有效，同时采用一种保持光 滑解高阶精度的通量限制器，抑制间断附近的数值波动，又不降 低光滑解区域的 精度。 对于空间离散现引入MUSCL方法，VanLeer提出的MUSCL方 法将相邻网格 2 塑塑查兰堡圭兰垡塑 节点出的雅插值到网格单元边界面处’定义单元边界面处中间 状态量为嚏、 嚏’采用以下线性插值得到’如下式( 钮2咯一一三 缸h 2-29 ?2吃+l su :!a 2-30 砌 乙 坼(厂坼川 认, 其中:北(， 是通量函数限制器，为了能在激波附近得到光滑间断解，防止 解在激波附近出现非物理数值振荡，本文中采用Roe’Superbee限制器，它具有单 调递增和对称的特性(形式如下: 嘶』 [O，rain 2,，1 ，fnin ,，2 】 2―31 其中:‘，:―Ul,j+|―--Ul,j 。 U1(J一吩(j-i 将重新构造的变量《一、屹毒分别替代一阶精度格式通量中的,J、坼?+-， 数值通量由原来的一阶精度变成二阶精度，所以其差分格式也为二阶精度，而且 解具有保持单调的特性。 2(2(6常规边界条件处理 实际问题的计算一般都是针对有限区域的，因此在区域的边界上需给定边界 条件，边界条件不能随意给定，它要求在数学上满足适定性，在物理上具有明显 的意义。在实际计算中，边界的处理很重要。 有限体积法具有三种不同形态的边界:开边界、内边界和陆地边界。因浅水 流计算中常用因变量为水深和流量，故边界条件主要有三种形式:给定水位、给 定流量、给定水位流量关系。 一、开边界 开边界可分为缓流开边界和急流开边界。对于急流开边界，其处理情况很简 单，因为其下游水力要素完全由上游决定，故不论是水位、流量或水位流量关系， 必须同时给定在入流边界上。以讨论缓流开边界的情况为例，未知状态的确定方 塑塑奎兰堡主兰垡丝兰 法如下: 1 给定水位 hz 1,R―z，r,R，z己知，根据黎曼不变量的关系?+2?万 叱+2?瓦得 到? 蚝+2?虱?瓦一0飘h vL; 2 定单宽流量 吼为已知单宽媳由方程组譬272届压一丙得孙I?5?十纠gt、,,一、,?J k2【 蚝+2?历- 一薏】2,49，给定?初始值，采用割线法迭代求解，求出,以 后，代入? ql,k，咋 屹。采用这种方法相当于平均分配流速，只对边界附 近有影响。 3 给定水位流量关系 己知边界处的水位流量关系靠 , , ，与给定流量关系情况相类似，可以 求解方程组: 譬,?。压厕，便可删相应的枷海„a 二、内边界 计算区域内有时候会存在水工建筑物，如堰、闸、堤、桥墩和涵洞等。当单 元边界与之相重合时，单元边上的法向数值通量的计算就属于内边界问题。对物 理边界上法向通量的计算要根据各种物理边界类型的不同区别对待。 三、陆地边界 单元体的交界面为陆地边界，表示边界处无水流通过。对于这种情况，设想 边界外面存在一个对称的虚拟控制体，令? ―叱，k ^工，缺点就是将陆边界作 为内部边界处理，利用静压假定但未作修正，适用于边界单元的流速近似与固壁 平行的情况( g,2+?I?I? 相容关系:圭鲥+吒kI吒 j1 h O,uR 0 河海大学硕士学位论文 陆地边界状态可表示为:? 吨，h 屹，k 吒 当局部弗劳德数C u,?历较大时，浅水方程组的静水压力假设在固壁处 不再成立，应考虑法向动量平衡的动水压力公式。 需要注意的是，以上处理中的越，v均表示局部坐标系中的法向和切向流速( 2(3LU-SGS隐式离散 2(2(7网格重新编号 一般来说，非结构网格的编号是随机的，这使得非结构网格在灵活性和网格 生成的贴体性方面有独特的优势。但是，在运用LIJ―SGS方法时，要求每个网格周 围编号比其大和小的网格都有，因此必须重排网格以满足这样的特定要求。有很 多文章介绍了针对非结构网格格点的重排方法M明，本文介绍一种基于非结构网 格格心的重排方法。 对任意一组非结构网格进行重排，在经过重捧、分组后。非结构网格应满足: 1 每组网格中的绝大多数都与其周围的网格相关，且周围网格的编号比其 大和小的均有: 2 每组中的每一网格都不与同组的其他网格相邻。 基于二维格心格式的三角形网格的重排步骤为 格点格式乃至更复杂的三维 问题同样适用 : a 将紧贴物面的一圈网格或者任意一端面的连续网格归于层I，任取其中 的一个网格编号为l，作为网格重排的起点。然后沿顺时针或逆时针方向将层I 中剩余的网格依次编号; b 把与层I相邻的各网格归于层?。根据网格的相邻性，依次把网格分别归 类到层?、?、V等。 c 检查是否所有的网格都已归属不同的层。 d 为在同一层面的网格中避免出现相邻的网格，需进行更细的层内“着色” 处理。 e 如果同一层内有相邻网格，就将其归入不同的“色层”，每一层内的“色 层”数由网格的相邻情况决定。各网格最终归属的不同层及各“色层”的情况 各 河海大学硕士学位论文 层的基础色层为0层，一旦同层内两网格相邻，就将色层数增加l，把其中一个归 为新的色层 。 f 根据每个网格所属的层号和。色层”号，对所有的网格 进行总体编号 从 层?开始 。为避免在重新编号后出现空穴，即某一网格湘邻网格的编号都大于 或都小于‘与从任意一个网格开始编号不同，本文从紧贴物面的一圈网格或者任 意一端面的连续网格开始编号，这样就增加了排序的鲁棒性。 在计算前只需进行一次网格重排，且不再需要矩阵的转置。这样在网格重捧 时，只需以较小的额外内存来存储排序后有关网格节点和边的信息，而因LU-SGS 法的计算过程和通量计算两者相互独立，所以内存可以共用。 为了方便理解，下面以四边形网格为例来简单说明本文的方法: 1 普通网格重排 a 将任意一端面的连续网格归于层I并进行顺序编号 图2(3 a 。 歉《_,鼍 涤:。;#蟹每w缈黼 l j6 „一灞 2 麓翻 嘲i一。滢 3 ?* 。’， 黧譬4灞 秽‘。;女也勰灞 4 „:， 璃 5 t1 ?_，蘧 “ ‘;’ i 《硪糖溺 图2(3 a b 把与层J相邻的各网格归于层?，根据网格的相邻性，依次把网格分别归 类到层?、?、V等 图2(3 b ( 。 鞲‘，饼甲“ |舅12”警爹’ 1 ? “? ? 弼 : 2 ? ’“? ’? 了 穗 3 ? ? ;? ; V通 4 ? 。? ? f” V萋 5 ? 。? 。? v蠹 棼4，?。 ??虢,,貅、 缀滋溯 图2(3 b 25 河海大学硕士学位论文 c 层内“着色”处理 图2(3 C ( 黪鬟。+鬻9《i 四?鬻”4 1 *。?墼。 》?1 ?辍?黔*零g骖势q 2 j?2、?2 魏; 誉鬻 3 菇H3一m3籍黼露澜 4 l?4，#?4。鬻融:，5 器瓣i 驾 5 ?5 。”m5 黔?旷 9诣?妻 驻戳翘僦鏊菇 盘莲濯 《鏊鑫rg城8?。一 图2(3 c d 根据每个网格所属的层号和“色层”号，对所有的网格进 行总体编号 图 2(3 d 。 1 6 ll 16 21 2 7 12 17 22 3 8 13 18 23 4 9 14 19 24 5 lO 15 20 25 图2(3 d 2 贴体网格重排 a 将紧贴物面的一圈网格归于层I并进行顺序编号 图2(4 a 。 棼#麟群i,#” ? ?辩《攀?#“ 嘲罐p?#‖帮掣@i??; 戮P"“者。 ((麝，V ‖。q脚磁g 《d„。?《，， 拼+、 《拶+ 器弼 I 16 1 2 3 4 l ; 嘉 g 。 一《 鍪! j 15 5 》 震 14 6 誓一 ”罐 器‘ !i j 一道 爹 j 13 7 萋’ 鸳 4 ， 舞 爹 12 11 10 9 8 《 _j ，。 l《?籀2瓣添、僦二漱翳狳+黢ti??褥1g教毋 蠡鑫蠢，城??臻氛 幽 图2(4 a 26 河海大学硕士学位论文 b 把与层I相邻的各网格归于层?。根据网格的相邻性，依 次把网格分别归 类到层?、?、V等 图2(4 b 。 虿 雕4_嚣瑟薹8蘸“曩聪鄹g 鬻阙 16 l 2 3 4 i 器，II一 ， , II鼋 静一 15 5” 誉6啊毳 14 鍪漫麓 ， 。。 6 鎏毋镬 13 7 誊，壁; i 。口4 菱4甄耋 擎口“ 12 11 10 9 8 纠I蔼 谶勰恩漱四里鏖纛陵磊点k。战。搬湖 图2(4 b c 层内“着色”处理 图2(4 c 。 ?5确麓酽嚣礴 ?毡譬 霉”宵爹 攀面9 慝, “,9 嘲 鍪。u19 16 l 2 3 4 ;iLl5蓬 15 5 赣 Iil8+ ;一I]6囊 罄1117(一 14 6 ;“117 l , 13 鏊囊-si ;， 7: 筹?s霪 篓?15。j 12 11 10 9 8 ;119礴 。々 i *、“一 ?14 II’讫 ，?IIL(1110 麴醯 ，以 :艇害‘ (j一 鞘 幽2(4 c d 根据每个网格所属的层号和“色层”号，对所有的网格进 行总体编号 图 』_ ， 37 36 17 18 19 20 38 ” 35 15 l 2 3 4 21 34 15