（流体力学专业论文）哈尔滨三家子段二维流场的计算与应用.pdf

哈尔演t 稃大学硕十学位论文 摘要 河流数值模拟是以数值方法和计算机技术为手段，通过对河流的数值模 拟计算，解决河流工程所关心的问题。随着计算机技术的高速发展，为河流 数值模拟提供了更强的发展动力。如今，使用一维数学模型和二维数学模型 进行河流的数值模拟计算已经很普遍，且模拟技术也比较成熟，并在许多方 面都取得了题大的经济效益。三维模型虽然最接近实际河流的特点，但由于 问题太过复杂，计算量很大，因此发展较缓慢。总之，数值模拟已成为一种 研究水问题的重要手段。 应哈尔滨城市发展建设需求，需对松花江三家子段堤防进行改建。本文 主要做了如下的工作：采用有限体积法对松花江哈尔滨三家子流域进行二维 数值模拟计算。使焉? 羹占体坐标网格剖分求解域，羹占体坐标系下的水流控制 方程采用有限体积法进行离散，使用f o r t r a n 程序编程对离散方程进行求解， 逶过已有的水文资料对结栗进行水面线验证。验证模型正确后，运用该模型 对三家子堤线进行优化，得出优化方案。 关键词：数值模拟；数学模型；帖体坐标系；有限体积法 哈尔滨下程大等乏硕十学位论文 a b s t r a c t b y t h ec a l c u l a t i o no ft h en u m e r i c a lm e t h o da n dt h ec o m p u t e rt e c h n o l o g y , t h e r i v e rn u m e r i c a li m i t a t i o nn o ws o l v e sm a n yp r o b l e m sw ec a r ea b o u ti nt h ef i v e r e n g i n e e r i n g r a p i dd e v e l o p m e n to fc o m p u t e rt e c h n o l o g yp r o v i d e st h em o r e m o t i v i t i e sf o rt h er i v e rn u m e r i c a li m i t a t i o nd e v e l o p m e n t i ti su s u a lt ou s e o n e - d i m e n s i o nm a t h e m a t i c a lm o d e lo rt w o 。d i m e n s i o nm a t h e m a t i c a lm o d e lt o p r o c e s st h ec a l c u l a t i o n so ft h er i v e rn u m e r i c a li m i t a t i o n , a n dt h ei m i t a t i o n t e c h n o l o g yi sv e r ym a t u r e 。a l s o ，i tm a k e sl a r g ee c o n o m y b e n e f i t si nm a n yp l a c e s t h r e e d i m e n s i o nm o d e la p p r o a c h e sr e a l i t yw e l l b u tf o ri t sc o m p l e xc a l c u l a t i o n a n dl a r g et i m e - c o n s u m i n g ，i th a r d l yd e v e l o p s 。i naw o r d ，n u m e r i c a li m i t a t i o ni s b e c o m i n g a v e r yi m p o r t a n tm e a n s i nt h es t u d yo f t h er i v e rp r o b l e m s t h i st h e s i sc a l c u l a t e ds a n j i a z io fs o n g h u ar i v e ri n2 dn u m e r i c a li m i t a t i o n w i 囊f v m 。t h eb o d y - f i t t e dc o o r d i n a t e sg r i di su s e di nt h ef l o wa r e a a n dt h e g o v e r n i n ge q u a t i o n so ff l o wi nb o d y - f i t t e dc o o r d i n a t e sw e r ed i s p e r s e dw i t hf v m t h ee q u a t i o n sw e r ee v a l u a t e d 蛾搬f o r t r a np r o g r a m s b yv e r i f y i n gt h ew a t e r s u r f a c el i n e ，i tt e l l st h a tt h e2 dm o d e li m i t a t e dt h ef l o wf r e e l y a tl a s tt h ev e r i f i e d m o d e lw a su s e dt oo p t i m i z et h ee m b a n k m e n tl i n eo f s a n j i a z if o rt h ec a l l i n go ft h e d e v e l o p m e n to fh a r b i n k e yw o r d s ：n u m e r i c a li m i t a t i o n ：m a t h e m a t i c a lm o d e l ；b o d y - f i t t e dc o o r d i n a t e s ： f v m i l 哈尔滨工程大学 学位论文原创性声明 本人郑重声明：本论文的所有工作，是在导师的指导 下，由作者本人独立完成的。有关观点、方法、数据和文 献的引用已在文中指出。除文中已注明引用的内容外，本 论文不包含任何其他个人或集体已经公开发表的作品成 果。对本文的研究做出重要贡献的个人和集体，均己在文 中以明确方式标明。本人完全意识到本声明的法律结果由 本人承担。 作者( 签字) ：盔纽物 日期：刈年多月7e l 哈尔滨t 程大学硕士学位论文 第1 章绪论 1 1 概述 随着科技的发展，人们对河流的利用鹾趋多样，对河流的认识也日益深 入，因此对河流的模拟有着巨大的实践意义和理论意义。模拟河流的方法有 物理模拟和数值模拟两种类型。物理模拟即比尺模拟，以其直观性和物理概 念明确一直受到水利工程界的重视，但其模拟程度往往受模拟试验手段和院 尺关系的限制，而且费用比较高，试验时间、周期也相对较长。河流数值模 拟是以挟沙水流运动_ ； 羹河床演变控制方程为基础，以数值方法和计算机技术 为手段，通过对河流的数值模拟计算，解决河流工程所关心的问题。随着计 算枧科学的迅猛发展，使数值模拟有了强大的生命力，它以其省时、高效、 无比尺效应等优点，在众多河流工程问题上被广泛应用，取得了良好的效果。 如在水力学计算中，用一维数学模型或二维数学模型进行河流、河阈、河翻 的数值模拟已很普遍，模拟技术已较成熟，并且融在流域规划、水环境保护、 河网开发、河口整治等方面取得了巨大的经济效益，因此，计算机数值模拟 己成为一种水科学研究的重要手段。 1 2 数学模型的发展 1 2 1 一维数学模型简介 水流数学模型发展至今，其方法已十分丰富多样。按模拟的物理量运动 变化维数可分为一维、二维、准三维和三维模型。对于长距离的河道问题， 使用维模型采用有限差分法比较有效；对于局部问题，如河口、海港和海 口，采用二维或三维模型利用有限元或特征有限元较为合适；而对于剧毒性 的水质污染问题需建立较高精度的数学模型。方法的选取和采用与不同问题 的实质紧密相连。一维模型主要用于研究长时期内长河段的水流及河床变形 情况，可以得到水流及河床的总体变化情况，例如修建水利枢纽后的水库淤 积计算等。 1 8 7 1 年，法国人s a i n t - v e n a n t 哺3 就建立起了描述河流维水流运动的 s a i n t - v e n a n t 方程组。然恧由于该方程组是复杂的鼍# 线性偏微分方程组，数 哈尔滨下挥大学硕十学位论文 i i i l l l i l l li l i i l i ii f 7 i l l- - - i | 置黑 学上用解析法求解存在一定的难度，早期受限予计算能力的不足，人们只能 采取简化方法进行求解。这些方法受简化假定的限制，计算精度不高，适用 范圈也不广。后来随着计算机技术的发展，许多学者开始利雳计算机对完整 的s a i n t - v e n a n t 方程组进行离散求解，取得了很大的进展。1 9 5 3 年s t o c k e r 啪。 在对o h i o 潺的洪水计算中蓄次采用显式差分格式求解s a i n t v e n a n t 方程组， 1 9 7 5 年c u n g e m l 给出了数种显式差分格式的表达式，并对计算结果进行了分 析。由于显式差分格式的计算稳定性存在时间步长限制，因此不受时闻步长 影响无条件稳定的隐式差分格式受到了人们的重视，并得到迅速发展。目前 基于p r e i s s m a n n 四点偏心隐格式的一维s a i n t - v e n a n t 方程数值解法在实际 工程中得到了广泛的应用。结合张二骏融玛等提出的河网非恒定流三级联解算 法及王船海、李光炽提出的求解河网大型稀疏系数矩阵的矩阵标识法，对大 型涛网水流的进行数值模拟技术上也趋予成熟。2 0 世纪? o 年代屠期以来， 许多国家在发展冲积河流数学模型方面都作了大量的工作。美国、荷兰、法 国、日本、前苏联等函已研制出相当多的一维模型n3 ，且已经较力成熟。 1 2 2 二维数学模型简介 二维模型有两种，一种是平面二维模型，它以垂线平均的水流因素作为 研究对象，研究它们在平面上的变化情况。另一种是立面二维模型，它以在 水流中截取的纵剖瑟上的水流因素作为研究对象，主要研究垂向扩散问题。 目前较常用的是平面二维模型。二维模型不管是在建立模型的基本方程和辅 助方程上，还是针对某一具体实际闻题所建立的模型及其应用方面，都要比 一维模型复杂、困难。二维模型最早是应用在河口、海湾、潮泊等宽阔水域 的模拟，经过十多年的发展已经广泛应用到河流中。欧美等国家二维数学模 型发展较旱，数量较多，发展比较成熟，謦前已有许多比较实用的商业软件。 在我国，平面二维泥沙数学模型的研究最早可以追溯到2 0 世纪5 0 年代 初期，但限于当时的计算水平，不褥不对模型进行了很多篱化，因此有很多 不合理的地方，功能较少，使其应用受到很大限制。但由于所研究问题的复 杂程度高，在这之后二维模型的研究并未像一维模型那样迅速发展起来。近 些年来随着三峡工程中的二维泥沙数学模型的研究被列入“七五、“八五 国家重点科技攻关项圈和“黄河治理与水资源开发利用 的研究全面启动， 2 哈尔滨工程大学硕十学能论文 誓皇i i i i 蔫麓f qh ij 1 1 黼 二维泥沙数学模型强7 。取得了迅速发展，已有不少解决实际问题的成功先例， 模型水平已达到了国际先进。长江三峡水库变动回水区及其下游河段的模型 以南京水利科学研究院溺、中国水利水电科学研究院淄、武汉水利毫力学院酞3 、 河海大学长江航道局啼1 等模型为代表。黄河上具有代表性的泥沙二维模型主 要有：张红武嘲、韦直林矗3 窦医仁溯、韩其为嘲、李义天n 稍、程晓陶n 玎模型等。 二维数学模型近些年来在我国发展较快，且已经广泛用于实际计算。二 维数学模型能模拟更加复杂的问题，但无论是数值计算方法上还是物理处理 上都要比一维模型复杂得多，计算时间上也要比一维模型长。大多数二维数 学模型都是模拟顺直河道或弯道。对予分汊河型，可以采用现有二维模型模 拟分汊及合流的情况，但还无法反映成型堆积体的运动情况。至于游荡型河 流，这方面的研究更少。对于深水问题，或流速、含沙量等物理星在水深方 向上变诧较大的闯题，二维模型无法给出满意的答复，这些就有赖于三维模 型的开发、应用。 1 2 。3 三维数学模型简介 三维模型可以模拟比一维、二维更加复杂的问题，也更接近实际河流三 维的特点。露外以美国密嚣薅吃娴立大学w a n g 和a d e f f 建立的三维泥沙模 型为代表m 1 。在国内，相比予一、二维数学模型，三维模型发展比较缓慢， 耳前还很不成熟。基本方程离散求解常用的数值计算方法有垂向分层法、函 数展开法、破开算子法、边界拟合坐标法、二三维耦合法等。虽然也研制出 了一些三维模型但在实际中成功应用的并不多。 综上，由一、二、三维泥沙数学模型的发展情况来看，三者所研究闻题 由易到难，因此发展的速度也是由快到慢，模型大体都经过了一个从初步研 究到实际验证的过程。各种模型都是针对不同的问题丽提出的，对子具体的 实际问题，单靠一种模型往往不能满意的解决，这就需要依据不同模型各自 的特点联合解决。 1 3 数学模型计算方法 由于摇述水沙运动的微分方程组是j 线性的，想用解孝厅方法获得河流各 种复杂定解条件下的解析解几乎是不可能的。随着计算机的发展和广泛应用， 3 哈尔滨丁程大学硕十学位论文 运用数值方法求解水运动方程得到迅速发展。而数值解的收敛性、稳定性和 与实际水流的符合程度，很大程度上取决于数学方程的离散形式，即采用的 数值解法。因两数值计算中，计算方法选择得是否合理是数值计算成败的关 键。 1 3 1 有限差分法( f i n i t ed i f f e r e n c em e t h o d ) 有限差分法( 简写f d m ) 发展较早，数学概念直观，表达简单，计算方 便，用于水解流体力学闻题有较长的历史，其解在存在性、收敛性、稳定性 早已有研究成果，是较成熟的方法。有限差分法是最早提出的方程离散方法， 但真正应用于流体力学的计算是在六十年代后，最初多用正规网格和松弛解 法，1 9 6 8 年引入交替方向隐式差分法，七十年代末提出上游加权有限差分法， 至今已形成许多成熟麴算法格式和软件。目前的研究前沿是非线性问题的有 效算法及其理论分析。有限差分法是用有限个网格点代替求解域、用差分代 替微分的途径将微分方程组转化为差分方程组，从而形成一个以各节点函数 值为未知量的代数方程获得其近似解。有限差分法的缺点是对各类方程都是 用一条直线上的有限点的差分代替微分，若格式不当将引起解的不稳定甚至 出现伪物理现象，而且边界上往往只能取前差或器差，只有一点的信息可以 带入，使得精度只有一阶。因而有限差分法是不适于复杂的不规则边界。7 0 年代，t h o m p s o n 汹3 等提出的贴体坐标系( b o d y - f itt e dc o o r din a t e ) 在很大 程度上弥补了这一缺点。常采用的有限差分法有：特征线法、a d i 法陟1 和破 开算子法n 4 1 。 1 3 1 1 特征线法 特征线法，在数学上是通过坐标的转换寻求可以将偏微分方程转化为全 微分的一组特征线，然后沿特征线积分求解。由于特征线法的计算嘲格不援 则，需要同时贮存节点在( x ，t ) 平面上的位置坐标和节点上的水力参数值， 要求贮存量大。所以，在实际计算中一般均采用特征线差分法，它综合了特 征线法和一般差分法两种方法的优点，网格的节点位置仍按照一般差分法布 置，但格点之间的关系按照特征线法原理建立。特征线法把时闻离散和空闻 离散起处理，其优点是能反映问题中信息沿特征线传播的性质，算法符合 水流运动的物理机制，计算精度高，不受稳定性的限制。由于该方法是沿时 4 哈尔滨工程大学硕十学何论文 间的推进求解，故较适于双曲型和抛物型问题，对于求解周期短、变化急利 的问题比较适宜。如用于涌潮地区，这是其它方法所不能比拟的。因为涌潮 地区，潮头破裂，沿x 方向发生了不连续现象，这时用一般差分法不能求解， 但它的特征线仍然存在，因此可用特征法求解。而推广到二维问题时，由于 二维问题表现为一个特征锥面，嚣前一般是对特征锥面选用凡条特征母线， 因特征线法求解格式复杂，尤其对高维问题更是繁琐，因此目前很少直接用 于数值汁算。 1 3 1 2a d l 法 a d i 法也称方向交替隐式，是由p e a c e m a n - r a c h f o r d 和d o u g l a s 于1 9 5 5 年提出的一种用于专解二维问题的特殊分步法。这种方法是将一个时间步长 分为两个半步，前半步在x 、y 项分别用隐、显格式，后半步在x 、y 项分别 用显、隐格式，这样徽可以把原来较大的系数矩阵化为两个三对焦系数矩阵， 可用追赶法求解。可以证明，a d i 法的隐式格式具有垃，缸2 ) 阶精度，而且是 无条件稳定的。壹于在差分格式中弓| 入中间步长，在每半个步长中只需解一 个三对角矩阵的代数方程，所以计算工作量比较小。但是a d i 法由于坐标方 向的限制，水域边界不得不概化成锯齿形，加上网格划分必须是等间距的， 因而在计算中，由边界网格化所引起的失实比较严重，而且重点计算水域不 能局部加密网格。a d i 法解题范围窄，不适用于三维情况，即使在二维的情 况下，也不适合解具有二阶混合导数的微分方程，因此a d i 法的应用亦有一 定的局限性。 3 。 3 破秀算子法 破开算子法，也称分步法，是由苏联学者y a n a n k o 1 等人，于5 0 年代末 至6 0 年代初提出的，并在理论上逐渐加以完善。破开算子法是数值计算中的 种分裂解法。在数学上，它将一个微分算子分裂成几个简单算予的线性组 合，在各分数步长上，分别求出各个简单算子的解，从而得出整个步长的函 数值。破开算子法，不像a d i 法，比较灵活，可以使不同组的方程各自按其 最佳的方式求解，甚至可将欧拉空间与拉格朗日空间结合起来，以改善其稳 定性。焉且破开算子法可以用于多维问题翡数值计算。但由于方程组的分割， 也形成了一些新的矛盾，反而影响计算精度。 哈尔滨f 挥大学颂十学位论文 1 3 2 有限元法( f i n i t ee l e m e n tm e t h o d ) 有限单元法是1 9 5 0 年由一些飞机结构工程邶酋先提出的一种数值方法， 6 0 年代有限单元法开始应用予流体力学，最初是联解一些线性问题，到7 0 年代已推广应用于非线性问题的研究。有限元法的基础是极值原理和剖分插 值，它吸收了有限差分中离散处理思想，同时采用了变分计算中选择逼近涵 数及对任意形状( 三角形或四边形) 的许多微小单元进行积分处理的合理方 法，爱而具有很广泛的适应性，特别适合子足何、物理条件比较复杂的阀题。 该方法具有较强的适应性，计算精度较高，但存在计算格式复杂、计算量及 贮存量较大，大型系数矩阵较难求解等缺点。在多维数学模型计算中，由于 有限元方法贮存量比较大，大型系数矩阵求解较困难，直接影响着计算速度， 因藕在非恒定问题及其它对计算速度要求比较高的问题中应用不是很多。 有限元法的求解方法常见的有直接法、变分法、加权余量法和能量平衡 法等n 舶。其中变分法类的里兹法( 1 9 0 9 ) 和加权余量法类的伽辽金法 ( g a l e r k i n ，1 9 1 5 ) 、最小二乘法珏霹常属于明渠非恒定流数值计算。对于天 然河流，其边界一般都极不规则，有限元法是一种理想的数值计算方法。囡 蠹张二骏n 5 1 用加权余量法进行长江下游流场计算；谭维炎等、耿光铨、史金 松等用之进行二维明渠非恒定流计算。对于隐式的有限元法计算量和贮存量 较大的特点，在8 0 年代初，酲本学者k a w a h a r a 等提出了“选择参数质量集 中二步显式有限元法 ，较大的减少了计算量和贮存量。k a w a h a r a 等运用 此方法成功的解决了河口、海湾水流等方面的问题，k a w a h a r a 法的不足之处 是时间步长取得较短。韦直林对这种方法进行了改进，对时闻项采用“预报一 校正迭代”的格式，在黄河的水库泥沙计算中也取得了较满意的成果。 1 3 3 有限分析法( f i n i t ea n a l y t i cm e t h o d ) 有限分析法是由美国i o w a 大学的陈景仁n 卅教授等人首先提出，是解析法 与数值法相结合的一种近似方法，其实质是把数学物理方程的解耩方法结合 到偏微分方程的数值解中去，因此有限分析法又叫半解析法。它克服了有限 差分法在求解不可压缩粘性流体时在大雷诺数下麴各种困难，避免7 差分近 似中的各种数值效应，是求解大雷诺数的各种流体力学问题行之有效的方法。 目前，有限分析法在求解n - s 方程、地下水运动等方面有较为广泛的应用。 6 哈尔滨下程大学硕十学位论文 有限分析法的基本思想是：首先将待求问题的总体区域分为许多小的子 区域，在这些子区域中将方程局部线性化后求局部解析解；然后从局部解析 解导出一个代数方程，使子域上的内结点值与相邻的结煮值联系起来，利用 在求解区域中已经离散了的子域之间的相互搭接覆盖，使子域与予域之间光 滑过渡，接着把搿有的局部解析解汇集在一超，就得到所求问题的有限分析 数值解。有限分析法有较高的计算精度，并具有自动迎风特性，计算稳定性 好，收敛较快。但由予有限分析法系数中含有无穷级数，单元系数计算较复 杂，这给实际计算与理论分析都带来了些困难。在8 0 年代末、9 0 年代初， 李炜n 7 3 等提出混合有限分析法，这种方法是在有限分析方法思想的基础上， 把建立精确差分格式与非定常项的差分处理结合起来，建立一种混合有限分 析格式。这种格式保持了有限分析法的优点，同时避免了无穷级数带来的不 便，褥到了广泛的应用。 1 3 4 有限体积法( f i n i t ev o l u m em e t h o d ) 数值计算从二维渐渐发展到三维，特别是紊流模式的弓| 入，使三维数模 进入实用阶段，有限体积法渐渐引起人们的注意，目前，有限体积法在三维 水流计算中的应焉已囊益广泛。有限体积法也称有限控翻容积积分法，它是 把计算区域离散为若干点，以这些点为中心，把整个计算区域划分为若干巨 相连接但不重叠的控制容积，将控制方程对每一控制容积积分，得出一组以 计算节点上的物理量为未知量的代数方程组来求解。有限体积法交微分方程 在整个计算区域上加权积分，这与有限元法相似，但其离散沿坐标方向进行， 形成的离散网格和离散方程更象有限差分法。由于采用守恒型的微分并对每 一计算单元进行质量和动量守恒形式的离散，使得微分方程包含的守恒性质 在每一控制体积上都得到满足，若设法保持各单元赛面两侧相邻控制体的计 算输运通量相等，那么整个计算上都能保持守恒性。有限体积法兼具有有限 差分法物理概念清嚷翻有限元法适应复杂边雾形状，计算精度高的优点。有 限体积法最早用于解决热力学问题，如s v 帕坦卡早在本世纪七十年代就创立 了计算热力学及恒定不可压流问题的s i m p l e 算法n 副。谭维炎等n 9 3 用o s h e r 格式进行过潮流计算和溃坝洪水演进计算。 l 。4 三家子海段简介 7 哈尔滨下程大学硕士学位论文 1 4 1 流域概况 松花江流域三面环山，中部为广阔麴松嫩平原，松花江的出口地区是著 名的三江平原。松嫩平原及三江平原都是洪涝灾害频繁而严重的地区，而嫩 平原恰恰又是哈尔滨、长春等城市及耕地密集的地区，因此洪涝造成的损失 巨大，洪水灾害一直是松花江流域的心腹大忠。 松花江流域洪涝灾害频繁，局部地区较严重的洪涝灾害平均2 - 3 年发生 一次，大范围严重的洪涝灾害1 0 一- - 2 0 年发生一次，2 0 世纪以来的1 9 6 9 年、 1 9 8 8 年嫩江犬洪水，1 9 5 3 年、1 9 5 7 年、1 9 6 0 年、1 9 9 1 年和1 9 9 5 年第二松 花江大洪水，1 9 5 3 年、1 9 5 6 年、1 9 5 7 年、1 9 6 0 年和1 9 9 1 年松花弧大洪水以 及1 9 3 2 年全流域大洪水都造成严重的灾害n 幻。 哈尔滨市是黑龙涯省省会，全省政治、经济、文化、科技中心，是一座 门类齐全、经济发达的新兴工业城市，也是我国东北工农为产品的集散地之 一，铁路、航运、公路的重要交通枢纽。哈尔滨市位于松嫩平原东南部，座 落在松花江干流上游段的末端，主要市区位于汪遂右岸。哈尔滨市的防洪任 务艰巨。 。 三家子堤防n 3 3 位予松花江予流哈尔滨城区下游左岸，是在原有民堤基础 上逐渐加高培厚形成的，是哈尔滨松北堤防的一部分，南接东方红堤弯曲向 北，北端接滨北铁路，全长1 0 2 k m ，堤顶高程l1 9 1 5 l1 9 5 1 m 。三家子堤 处在松花江与呼兰河交汇的三角洲地带，属于松花江漫滩地，地势低平，l 霸 河汉、串沟较多，地面高程在1 1 4 0 0 - 11 7 o o m 之间。 松花江三家予堤段河道条件较复杂，上有阿什河，下有呼兰河支流流入松花 江。千流自滨洲桥至呼兰河口主槽处在“s 型弯道江段，下战备路处在“s 型弯道的中部。下战备路上游，在两岸堤防和水泥厂等固定节点的作用下， 江道基本稳定。下战备路下游，两岸堤距从3 8 4 k m 逐渐扩宽，到呼兰河口处 达到11 。o o k m ，过呼兰海日下游后又收缩到4 。5 0 k m 。松花汪于流的溺道比降 上缓下陡，平均为o 0 6 ，其中三家子堤段的c s 松4 0 一- - c s 松3 6 断面比降 为0 。1 3 。 1 4 2 地形资料 河道地形范重，松花江主河道从上游双城县林场、西姚家岗至下游 、家 哈尔滨t 程大学硕十学位论文 口村；呼兰河上游至团山子、岔拉火烧。本文采用的地形资料测量时间为2 0 0 2 年1 0 月 - - 2 0 0 3 年2 月，比例尺为l ：1 0 0 0 0 ，等高距为l m 。采用1 9 5 4 年北京 坐标系，1 9 5 6 年黄海高程系。 随着城市建设的发展，一些工程项目的兴建，河道地形也会有所变化。 2 0 0 5 年开工，2 0 0 7 年建成的新建啥尔滨绕城高速公路东段松花江大桥( 四环 桥) 位于松3 8 一l 断面东战备路上游6 0 0 m 处，从北岸东方红堤开始，桥长 2 2 7 5 m ，南岸为路堤。呼兰河从下到上有滨北铁路桥、哈呼公路桥和哈伊公路 桥，间距分别为5 5 k m 和3 5 k m 。哈呼公路2 0 0 2 年扩建后，主桥长2 9 4 m ，头 道滩桥长1 4 0 m ，二道滩桥长1 5 0 m ，总长5 8 4 m 。 流域地图见下页。 9 ：鎏士耋銮：翟老：垩兰苫 哈尔滨t 程大学硕十学位论文 。5 本文的主要研究内容 随着哈尔滨市的发展需要和行政区的重新划分，三家子堤段内将成为城 市的重要开发区，为保护铁路和形成完整的松北防洪体系，需改建为主堤。 三家子堤防是松干堤防建设目标中唯一未达标的堤段，直接威胁城市安全。 针对以上背景本文研究的主要内容包括： 1 对此段流域通过有限体积法建立数学模型； 2 把计算区域离散为以若干点为中心的控制容积，将控制方程对每一控 制容积积分，得到该点上物理量为未知量的代数方程组； 3 通过f o r t r a n 编程对这些方程组进行求解并验证，得出模型是否正确； 4 。运用验证后的模型对三家子段堤线改造的两套方案进行计算并分析比 较方案中的不足； 5 针对方案中的不足之处，对三家子段堤线进行优化。 哈尔滨t 程大学硕+ 学位论文 第2 章二维数学模型 2 1 二维曲线正交网格 2 。1 1 概述 实际工程问题往往具有复杂的几何边界，极其不规则。采用传统的方法 在笛卡尔坐标系嘲下将求解域用矩形网格划分对于模拟复杂几何边界会很困 难，该方法通常将边界概化为锯齿形状，这在模拟精度上存在一定的问题， 焉且关键区域的计算精度仍然得不到保证，此外需要占冀l 大量的内存来存储 边界数据。这就需要采用坐标变换，来建立其他坐标系下的网格。例如极坐 标系、正交曲线坐标系等。 坐标变换的基本思想是采用适当的计算格式，构造出一个规则的计算区 域，是该计算区域的坐标轴与曲线边界一一对应，生成帖体坐标晗钿系。其目 的就是减小数值误差，保证勰的稳定，节省计算时间，同时又能使边界条件 容易处理。 坐标变换时注意以下几个方面：首先，使诗算区域边界落在坐标线上。 这样能使边界条件容易处理，从而减小数值误差；其次，物理平面上的网格 点废与计算平面上的网格点一对应，同一簇中的蓝线不能相交，不同簇中 的两条曲线仅能相交次；再次，在物理平面上，计算隧域内部网格的疏密 程度要易于控制。在帖体坐标的边界上，网格线最好与边界线正交或者接近 正交，以便与边界条件的离散；最后，在物理量梯度大的区域，网格需密些， 以便边界条件的离散。 2 。1 2 生成帖体坐标系的基本原理 生成帖体坐标的方法主要有三种：复变函数法，代数变换法和微分变换 法。复变函数法和代数变换法要求已知曲线边界与坐标之间的函数关系。对 于复杂的几何外形，在多数情况下，这种条件是无法满足的。而微分变换法 则适合于复杂不规则边赛。通过微分方程的变换关系，褥到的坐标系麓理想 地拟合曲线边界。微分变换法通常有两类：l a p l a c e 方程变换和p o i s s o n 方 程变换。透过求解这两类方程变换中的种，即可得到所需要的坐标系。 1 2 哈尔滨t 稃大学硕十学位论文 i it t i 宣i i i i i 鼍 帖体坐标系的生成，实际上是在物理平面上求解l a p l a c e 方程和p o i s s o n 方程的边值问题。这种方法最早由w i n s l o w 提出。以后很多学者对该方法进 行了研究。比较系统全面的研究这一方法的主要有t h o m p s o n 口等人。他把这 种方法首先应用于二维物体绕流的数值计算。 描写边值问题最简单的方程是l a p l a c e 方程，对于二维问题而言，在物 理平面( 工，) ，) 上，只要给定合适的边界条件，通过求解两个l a p l a c e 方程，就 可以得到正交的边界拟合坐标。由于l a p l a c e 方程变换不能控制网格的疏密 程度，所以这种方法用得不是很广泛。在实际计算中，在物理量梯度较大的 区域，需要加密网格，以保证求解精度；对物理量变化较小的区域，需要稀 疏网格，以便节省计算时间。为达到这个目的，一般使用p o i s s o n 方程变换 跚3 生成帖体坐标系。 设( 石，y ) 为物理平面上的笛卡尔坐标系，( 善，7 7 ) 为相应计算平面上的变换 坐标系。采用以下的数学方法以可实现两个坐标系之间的变换，即： 孝= 孝( 工，y ) r = 刁( x ，y ) ( 2 1 ) 若采用p o i s s o n 方程进行变换，则有： v 2 善= p ( x ，y ) v 2 刀= q ( x ，y ) ( 2 - 2 ) 上式中的p 停，刁) 和q ( 善，刁) 为控制函数，用于控制网格的疏密程度。若 尸( 善，刁) 和q ( 善，7 7 ) 均为0 ，则为均匀网格，上式就变为l a p l a c e 方程： v 2 孝= 0 v 2 刀= 0 以上p ( 善，7 7 ) 和q ( 善，刁) 的控制函数有多种取法。 比较简单的取法如： p ( 善，7 ) = 一口，e x p ( - d ，i 善一喜i ) 1 3 ( 2 - 3 ) 哈尔滨t 稗大学硕士学何论文 iv q ( 善，町) = 一口e x p - d 加一r 川 ( 2 4 ) 其中口和d 为常数，毒，巩为某一给定的孝，刀线。p ，q 在磊，r 6 上达到最大值，向毒，r j 线的两边以由d 控制的速率衰减。这样的函数尸， q 将使善，7 线向毒，刁，线的一边集中，而向另一边散开。但如果采用变号 函数： 尸( 厶7 7 ) = - 口，s 咖( 孝一磊) e x p ( - d , l 孝一磊j ) q 皓，刁) = - a j s i g n ( 刁一r l a e x p ( - d s r 一，7 朋( 2 - 5 ) 则善，r 线将集中在，r j 线的两边。在这种方法中，“口 代表坐标 线的密集强度，“d 代表衰减因子。在一般的坐标变换中，口= 1 0 较弱， 口= 1 0 0 较为中等，口- - 1 0 0 0 则较强。d = 1 0 时密集的范围较小，d = 0 1 时包 括的范围就比较广。 t h o m p s o n 的方法中p 、q 函数表达式为： p = - - e a i s 劬( f 一量) e x p ( - c ，悟一当1 ) 一羔l s 忉 一身) e x p f 一嘭0 万= i f 玎而 j - - - i 、 7 q ：一窆叩劬仂一，7 ，) e x p ( - c 肋一仉i ) = 1 一扣劬h ) e x p ( - d j 厄丽) ( 2 _ 6 ) 上式中：a ，b 为密集强度；c ，d 为衰减因子；力为网格线数目；m 为 网格数。但t h o m p s o n 调节因子可能存在下列几个问题3 ：1 ) 可调整的参数过 多，当网格生成质量不理想时，很难确定调整哪一个参数才能改善网格的质 量：2 ) 网格生成迭代过程中不易收敛；3 ) 参数选取不当易引起双映射现象； 4 ) 不能保证生成的网格在边界处正交。 1 4 哈尔滨t 程大学硕十学付论文 葺i i i 置i i i i i i ii i i i i i i i i i i 宣i i i i i i i i i i 宣i i i i i i i 葺i 一 对于河道二维计算的曲线网格，在选择尸( 孝，刁) 和q ( 孝，刁) 的收缩因子时， 一方面需要考虑曲线网格的疏密分布，另一方面要使变换后的曲线网格正交 3 l 】 同时采用曲线网格的疏密分布与正交变换形式，将p o i s s o n 方程改写为 3 2 】 式中： k = ( 2 - 7 ) 其中，孑，六，印勤等分别为票，薯，妻，高等的简写。为了 便于计算又不失一般性，取孝= 刁= l ，则有 c f 丽。c ，= 丽( 2 - 8 ) 式中c ，、e 分别表示石、】，平面上曲线网格的长和宽，因此k 与曲线网格的 长宽比成比例。由流函数与势函数的正交性决定了根据以上方程变换后曲线 网格的正交。为了便于比较，可将上述方程改写为： a 2 善a 2 孝一a 善0 l n ( k ) a 善0 l n ( k ) 孤2 。o y 2 一孤玉砂砂 望+ 塑：一o r ：1 ln(1k)orlol n ( 1 k ) ( 2 - 9 ) 孤2 却2孤舐却巩 7 方程右边为收缩因子p 、q 。 可将以上方程进一步变换如下形式： c ；x 鹭+ c ；x i | l + jz b ；p + x q q l = 0 o = 一、nv一矿 翔也啼畦妊 a 一砂 a 一砂卜，=叫 矽瓦 却i 尼 1 一后 ，、 a h a 一孤 哈尔滨t 程大学硕十学何论文 c ；j ，嚣+ c ；y ，7 叩十j 2 眵f p + y 口e - - 0 ( 2 1 0 ) 上抻岛一吉掣，q = 吉掣。 文献 3 6 指出如果计算域四周边界上的网格点选择不当，则采用上式生 成的曲线网格的正交性精度不高。特别是对于河道，网格的长宽比相差较大， 单靠收缩因子p 、q 是很难使曲线网格达到足够精度的正交性。为了提高曲 线网格的正交精度，可以不把四个边界上的网格点都固定死，而是允许其中 一条岸边上的网格点，根据正交的要求，沿边界移动，这样不仅提高了正交 性的精度，且迭代收敛速度也加快。其基本思路简述如下：用迭代法求解方 程组，经刀次迭代后的解为x ( f ，) 、j ，( f ，) 。设7 7 = n 岸边线上的网格点位置 石( i ，) 、y ( f ，) 根据正交性要求调整，孝= f 线与岸边线r l = n 交点p 处的交 角口为：c o s 矽：f 兰墨掣1 ，因此，需要调整的角度为万：9 0 。一0 。网 l ， p 格点尸的位置应从工( f ，) ，y ( f ，) 沿边界向左边( 当万为负则向右边) 移动的 距离为d s = 砸c 玎1 8 0 。 2 1 3 坐标变换的数值计算 ( i ，j + d ( 卜- 1 ，j )( i ，j )( i + 1 j ) 、 、 ( i ，j 一1 ) 图2 1 离散节点示意图 对方程进行离散时需用到5 个节点，如图2 1 。一阶导数，在计算区域 内点采用中心差分，对边界点采用前差与后差。二阶导数采用中心差分格。 其二阶导数差分格式为： 嘞= 血铲嗝= 丛笋( 2 - 1 6 哈尔滨t 程大学硕十学位论文 y 劈= 业铲仙叩= 丛铲( 2 - 1 2 ) 求解以上两式有： 掣蜗q ) + 掣 5 1 虿1 虿一 虿+ 希 掣t j z 善p + y ，q ) + 掣 ”2 菘面r 一 莓+ 芹 ( 2 - 1 3 ) ( 2 - 1 4 ) 对上式可采用松弛迭代法进行计算，即有： 矽+ 1 ) = 口万( “1 ) + ( 1 一口切) ( 2 1 5 ) 上式中，口 l ，口 1 0当ip | l o 时扩散项可忽略不计。 c = 0 q 丸7 7 ，f i = 0 q ) 。a t ，c = p ql 善，c = ( v q 尢孝， 见2 。等朋。= 吲。是， 驴蚓。筹棚，= 蚓，筹。 式中：u 。，甜，。，咋分别为控制体垂直面上的流速，c ，l ，l ，l 为控制体面上紊动扩散系数，g ，q 布置在主网格点上，动量网格点上g ， g 采用相邻自然网格点线性插值，色，六，玑，玑为相邻节点距离。 由上述各式整理得结点p 上的离散化代数方程，经整理后可以得到节点 尸上的离散化代数方程： 口| p p2 口中e + a w 中矿+ 口j v 中+ a s s + 6 ( 2 - 3 0 ) 式中：= 囟。彳0 只i ) + m a x ( - c ，o ) k e ，唧= 【d 。彳0 r 1 ) + m a x ( r ，o 归矿 = 囟。彳0 只i ) + 瞰( - c ，o 肛，= 囟。彳0 e 1 ) + m a x ，o ) k s 口尸= 口e + 口矽+ 口5 + 口一s j p h 尸考7 7 ，b = s c 曰j p 辱，7 。 为了避免产生波状流速场和压力场，计算u 、 ，、h 时采用交错网格技 术，将速度的网格点布置在主网格点控制体积的交界面上 根据离散方程( 2 - 3 0 ) ，只有在水位分布给定以后才能求解，若给定的水 位不正确，则计算出的速度场不满足连续方程。本文在计算水位时采用 p a t a n k e r 水深校正法。令u = u + “，= ，+ v 7 ，h = h + h ，1 , 4 。，+ 及h 为求出的速度和水深值，甜7 ，y ，h 7 为速度和水深的校正值。即： 吼“：= 口：材二+ 以+ g 峨日l 似；一畦冷刁 ( 2 3 1 ) a w u * = 口二“二+ 6 。+ g 眩h ) 。似；一睇玲刁 ( 2 3 2 ) a n v * 。= 口二v 二+ 6 。+ g 忆日l 一碥冷f ( 2 3 3 ) a s v ；= a n b v n b + b s + g 睡日l - h ；h 孝 ( 2 3 4 ) 在离散的连续方程中代入u ，1 ，可得到水深校正方程。 口；办；= 口：办：= a w h + 口0 办：r + a s h s + 召 ( 2 3 5 ) 其中：口：= 龋，口= 畿 ，g ( c f 厶f ) ：，g ( ( j 月r 孝e 。瓦忑习习一s2 再忑巧习 a p 7 = g l e 7 + 口十口二，+ 曰= ( c h u ) 。一c h u + l j 刁+ h v i ( c c h v x 丛善 0 q 必7 7 一h 尸) r 一 址 h p = 酢十砟 ( 2 3 6 ) z，。=“：+端依；一：)，甜。=z，：+：g。(一h乙cl口)w二a爿rl二(h；一砟) ( 2 3 7 ) =y：+糍(；一。)，屹=v；-t口，g一(h乙c善口)二a爿孝厶(；一向；) ( 2 - 3 s ) 本文在求解离散方程时采用p a t a n k a r 和s p a l d i n g 啪1 给出的欠松弛技术， 块校正技术及逐行迭代法。 哈尔滨t 程大学硕十学位论文 i r1 i - 1 i 鼻| 2 2 5 流场计算步聚 整个流场计算步骤为： ( 1 ) 根据初始假定流场及经验确定初始水位场h ； ( 2 ) 求解动量离散方程( 2 - 3 0 ) ，得甜，+ ； ( 3 ) 计算( 2 - 3 5 ) 式中的口，检验u + ，1 ，+ 是否满足连续方程，若满足则执 行步骤( 7 ) ，否则继续下一步； ( 4 ) 求解水深校正方程( 2 - 3 5 ) ，得h 7 ； ( 5 ) 由式( 2 - 3 6 ) ，( 2 3 7 ) ，( 2 - 3 8 ) ，得到h ，甜和v ； ( 6 ) 利用改进后的速度，重复第( 2 ) 步到第( 6 ) 步： ( 7 ) 结束计算。 计算流程图示见图2 4 ： 是 图2 4 流场计算流程图 哈尔滨丁程大学硕十学位论文 2 2 6 流场计算控制及定解条件 作为恒定流计算，给出入麓的流量酝和出口的水位。根据选取断面的水 位及由上至下选取的测流控制断面的流量满足以下条件就认为计算稳定，并 结束计算。第一条件为各流量控制断面与q 0 的误差在3 以内，即 1 2 一酝| 姨3 ；第二条彳宰为选取断颈的计算流量满足第一条件，两次迭 代之间的最大水位的变化值s 小于0 0 0