（环境工程专业论文）插值法在二维水质模型数值计算中的应用.pdf

摘要 y 3 9 9 | ；7 ( ，在预测复杂的河流段面的b o d 、c 0 d 等污染物的数值时，利用有限 单元法进行数值计算处理复杂的二维水质模型得到了较广泛的应用。用 有限单元法把要研究的水质模型按流带和断面划分若干个网格，通过计 算每个网络中水质的变化来预测要求的断面的水质状况，这种方法，能 使复杂的水质模型变得简单化。 但是根据有限单元法原则，要得到准确的预测值必须尽可能多地选 取网格，因而造成了大量的数据需要处理，通常的作用法是通过读取数 据作出累积流量图，从而确定流带形状，再通过f o r t r a n 语言程序计算 出要求断面的水质预测值。这种作法不仅计算量大而且在作图读数据 时，存在系统内误差，往往使模型的水质预测值的精确性下降卉7 本文在有限单元法研究时，引入了特别针对矩阵的m a t l a b 语言处 理大量数据，从而简化了数据的计算，同时，用插值法取代了作图法， 避免了一定的系统内误差，给水质模型的预测增加了准确度，还节约了 大量的人力和物力。 在m a b t l a b 语言插值法中利用有限元来处理二维复杂河段的水质模 型，为有限单元法提供了一种新的工具，同时在二维水质模型的研究和 环境评价当中具有较高的应用价值。 关词二维水质模型，数值计算，插值法，m a t l a b 语言 a b s t r a c t n 坞f i n i t ee l e m e n tm e t h o dd e a l i n gw i t ht w o - d i m e n s i o n a lc o m p l e x w a t e r q u a l i t ym o d e l s i sw i d e l ya p p l i e dt op r e d i c tt h en u m e f i c a lv a l u e so f t h e w a t e rp o l l u t a n t ss u c ha sb o d ，c o do fc o m p l e xf i v e rs e c t i o n s a c c o r d i n g t ot h ef l o w i n gb e l t sa n dd i f f e r e n ts e c t i o n s ，t h ef i v e rc a nb ed i v i d e di n t o s e v e r a lg r i d sb yt h ef i n i t ee l e m e n tm e t h o d 1 1 1 ec o m p l i c a t e dw a t e rq u a l i t y m o d e lc o u l db es i m p l i f i e db ym e a n so fc a l c u l a t i n gt h es u b t l ec h a n g e si n e v e r ye l e m e n t ，t h e r e f o r e ，t h ew a t e rc o n d i t i o no f t h er e q u i r e ds e c t i o nc o u l d b ep r e d i c t e db yt h er e s u l t s c o n s i d e r i n gt h ep r i n c i p l eo f f i n i t ee l e m e mm e t h o d ，t h eg f i d sm u s t b es e l e c t e da sm a n ya sp o s s i b l ef o rt h es a k eo fa c c u r a t ep r e d i c t i o n a sa r e s u l t 1 0 r so f d a t as h o u l db ed e a l tw i t hs i m u l t m e o u s l y u s u a l l y t h ef l o w i n g b e l t s s h a p e c o u l db eo b t a i n e db yr e a d i n gt h ed a t ef r o mt h eg r a p ho f a c c u m u l a t i n gf l o wv o l u m e ，a n d t h e nt h er e s u l to f p r e d i c t i o nc o u l db eg a i n e d b yt h eu s eo ff o r t r a nl a n g u a g e h o w e v e r , t h et r a d i t i o n a lm e t h o do f t e n l c a d st oa1 0 wa c c u r a c yb e c a u s eo ft h ee n o r m o u sc a l c u l a t i o na n ds y s t e m a t i c e r r o r sd u r i n gt h ed a t aa c q u i r e m e n t t h e p a p e r u t i l i z e dm a r l a b l a n g u a g e w h i c hi se s p e c i a l l ys u i t a b l ef o r m a t r i xc a l c u l a t i o n ，t od e a lw i t hl a r g ed a t ai na p p l i c a t i o no ff i n i t ee l e m e n t m e t h o d i tc o u l dn o to n l ys i m p l i 母t h ec a l c u l a t i o n ，b u ta l s oa v o i ds o m e s y s t e m a t i cc 玎0 r sb y t h er e p l a c e m e n to fi n t e r p o l a t o rw i t hd ra 删g g r a p h c o n s e q u e n t l y , m o r ep r e c i s e r e s u l ta n de c o n o m i cc o s to f m a n p o w e r a n dm a t e r i a lr e s o u r c e sc o u l db ep r o d u c e d n l ea l g o f i t h i n w h c hu s e sm 舡l a bt os o l v et h et w o d i m e n s i o n a lf i v e r s i m u l a t i o np r o b l e m s ，p r o v i d e dan e wt o o lf o rt h ef i n i t ee l e m e n tm e t h o d ，a n d i ta l s oe x e r t sa h i g hp r a c t i c a lv a l u e i nt h er e s e a r c ho f t w o - d i m e n s i o n a lw a t e r q u a l i t ym o d e l s a n dt h ee n v i r o n m e n ti m p a c ta s s e s s m t , “n t s k e y w o r d st w o d i m e n s i o n a lw a t e rq u a l i t ym o d e l s ，n u m e f i c a lc a l c u l a t i o n ， i n t e r p o l a t o r , m 枷a bl a n g u a g e 第一章概述 1 1 河流水质模型 1 1 1 河流水质模型的介绍 河流水质模型是描述污染物在水体中随时间和空间迁移转化规律的数学 方程，它是目前环境系统数学模型中发展最早，应用最广的一种。如果从1 9 2 5 年s t r e e t e r p h e l p s 模型算起，水质模型的建立和应用已有五十多年的历史。 目前已出现了为数众多的各种水质模型。 从使用管理的角度来说，水质模型可分为江河模型、河口( 受潮汐影响) 模型、湖泊与水库模型、海洋模型等。一般河流模型比较能反映实际，湖、 海模型比较复杂、可靠性小。 从模型的水质组分来说，可以分为单组分、耦合组分、多组分和水生生态 模型等。应用得最多的是b o d - d o 耦合模型。这是因为它对水污染控制具有普 遍的重要性，而且它的模型已能较为真实地反映实际。对于水温，细菌，一 般碳、氮有机物，以及工厂排出的简单结构的有机化合物( 如酚等) 等这些 组分的模型已达到实用化的程度。对于重金属和许多复杂有机化合物和毒物， 由于对它们的各种反应过程的认识比较差，因此模拟还很困难。对于各种含 磷营养剂在湖泊中的非线性时变反应，虽然己作了很多研究，具备了模拟的 初步条件，但这些模型的建立要求很多的数据、时间和较高的技术。近来已 提出了比较复杂的多组分水生生态模型，它更为详细和综合地描述了水体中 存在的许多成分及其相互作用，需要更多的数据和计算时间。 从系统的状态来说，水质模型可区分为稳态和非稳态两类。两者的区别在 于水文情况和排污条件( 因而构成河流的水质状况) 是否随时间变化。我们 最关心的是河流的临界条件，即水质最不利的条件。一般在河流低流量时( 如 晚夏或早秋) 往往也是水温较高的时候，此时水质接近临界，水流状态则恰 恰属于稳定时期。因此一般河流和河口可以采用稳态的水质条件来进行水质 的模拟和规划。动态水质模型可用来估算暴雨径流和污染物事故泄漏等瞬时 变化的情况，适用于短期性的水质管理和控制，要求有较多的统计数据。同 时，由于输入河道中的点、面污染源的时变状况常被忽略，因此也就降低了 动态水质模型输出的可靠性。 从水质模型的空间维数来说，虽然所有的真空系统都是三维的结构，但在 使用上往往采用零维、一维或两维的水质模型已足够了。水质模型的空间维 数，主要取决于所研究的范围及其水体中污染物的混合情况。如果对区域性 水质进行粗略的规划估算，零维模型也就可以了。如果对一个较长的河段或 河流进行水质规划，一维模型已能得出较好的结果。如果要研究河流局部河 段范围内的水质情况、排污口附近的污染物分布以及污染带这类问题时，就 要考虑采用二维、甚至三维的模型。两维模型适用于某一坐标方向混合比较 均匀，例如竖向混合( 一般相对浅而宽的河流) 或横向混合( 如在分层的河 口) ，而其它二个坐标方向的浓度梯度较大的情况，本文第四章的实例就是一 个二维水质模型应用的例子。 从水质模型的确定性而言，虽然绝大多数数据本身要求随机或概率性的模 型。但随机性模型的识别，要求河流水质各种变量采用概率分布的数据来定 量，而不能采用它们的期望值和平均值，这是非常困难的，因此目前绝大多 数采用确定性模型来进行水污染控制的模拟和规划。 河流水质模型的建立，广泛地应用于环境影响评价的河流水质的模拟预测 中。一种较科学、优化的水质模型，能提高环境评价的准确度和简化计算过 程。 河流、河口、湖泊( 水库) 、地下水和海洋等天然水体虽有各自的水质模 式，但大同小异。污染物进入水体后，在水体中产生一系列的物理、化学和 生物化学等作用，根据其在水体中不同的稀释、扩散状况，可分别由空间不 同的维数推导出1 1 2 节中各种河流水质模型的基本方程。 1 1 2 河流水质模型基本形式 l 、零维水质模型 当一个水体呈完全混合状态可建立零维水质模型，比如：一段相当静止 的河流、一个湖泊、一个水库或一个局部水域等。 在某些特定的条件下，可以将所研究的水环境单元视作一个完全混合的 反应器。在这样一个反应器里，不存在环境质量的空间差异，即在任何一个 空间方向上都不存在环境质量的变化。图1 - 1 所示为一连续流完全混合反应 器，进入反应器的污染物能在瞬间内分散到反应器的空间各部位。 根据质量守恒原理，可以写出完全混合反应器的平衡方程，即零维模型： 2 y 百d c = q c o q c + s + r y( 1 1 ) 式中y 一反应器的容积；q 一流入流出反应器的物质流量； c o 一输入介质中的污染物浓度；c 一输出介质中的污染物浓度，即反 应器中的污染物浓度；r 一污染物的反应速度；s 一污染物的源与汇。 s q 图i - 1 零维模型图示 对于一个没有源、汇项的反应器( 如湖泊、水库) ，即当s = o 时，上式 可以写成： 矿箜d t = q c c o c m ( 1 2 ) 如果污染物的反应符合一级反应动力学的衰减规律，即，= 一k c ，式( 卜2 ) 可以写作： y 百d c = q ( c o - c ) 一k c 矿 ( 1 - 3 ) 式中k 一污染物衰减速度常数。( 1 - 3 ) 式是零维模型的基本形式，在湖泊、 水库等水质模型中广为采用。 2 、一维水质模型 当污染物在河流的横向、竖向均匀混合，比如污染源位于研究的河段上 游较远时，一般要进行一维河流的水质的模拟和预测。 第一章概述( i 页一8 页) 一维模型是通过一个微小的体积元的质量平衡推导的，在这个体积元中 只在一个方向( 设为x 轴向) 上存在浓度梯度。 图卜2 表示在x 方向上存在输入、输出的微小体积元，其边长分别为a x 、 缈和z 。 卜一纛缸 ，y 图1 2 体积元的质量平衡关系 石 卜( 一皿别驰 。， 叩+ - 警- a x + ( 一o x - 篆+ 丢( 一e 酬挑 m 。， 跹b姆位(1-6) 甏麓羲m 蛐h ，咿+ 警酬一眈争丢( 一所争缸】姚一黜姚“ 4 第一章概述( i 页一8 页) 鲁= 一警一去( 一砬篆) 一k c s ， 在均匀流场中，“，和皿都可以作为常数，则上式可以写作： 箜：d 。箕一“，丝一愆 ( 1 9 ) a t 4 a ra x 式中c 一污染物的浓度，它是时间t 和空间位置x 的函数；取一纵向弥散 系数；u 。一断面平均流速；k 一污染物的衰减速度常数 式( 1 - 9 ) 就是均匀流场中一维模型的基本形式，一维模型一般应用于河流 水质的模拟、预测。 3 、二维水质模型 如果需要模拟的河段较短，或宽度较大，污染物在宽度方向上的浓度梯 度较大，就要进行纵向和横向的模拟，必须建立二维水质模型，比如大型河 流、河口、海湾、浅湖等。 在沿河岸有排放口的河流水质中，污染物的排放在横向和纵向很不均匀， 是不断变化的。这时我们往往首先可认为在水深方向的污染物是均匀的，而 去研究水体中污染物在横向和纵向两个方向的二维水质数学模型。 二维模型是通过选取河流水体中一个微小体积元，在这个体积元中只存在 x ，y 方向的浓度梯度。利用质量平衡的原理推导出来的。 图卜3 表示在x ，y 方向上存在输入、输出的微小体积元，其边长分别为 ) 【，y ，厶z 。 盘 第二皇苎姿! ! 蔓二! 蔓! 一 单位时间内输入该体积元的污染量为： 阪c + ( _ d - 塞) a y a z 山p ( - d y 茜1 伽5 ( 1 - 1 0 ) 单位时间内输出的污染量为： 帆c + 警x + ( 一。x 襄) + 丧( 一。x 塞) x 】y z + m ，c + 巧a ”f y c y + 卜d y 考) + 昙( _ 。y 考y ) 拾【z ( 1 1 1 ) 若污染物在微小体积元内发生一级衰减反应，由衰减输出的量为： k c x v z( 1 1 2 ) 于是，单位时间内，该体积元的污染物输入和输出遵循质量平衡的原理， 有如下式： 害鸠& = 【u x c + c + ( 一m 塞) 】匈位+ 【u y c + c - 印考) 】x & ) 一 扣= c + 警缸+ c 一皿争+ 丢c m 萋刎匈垃 + - 等- ( a y + c 一印+ 茜卜缈考，她 。一。， 将上式简化，并令x o ，a y 专0 ，得： 鲁= 一警一去c m 塞，一警一专c 一跏争一k c。一，。， 在均匀流场中，肛x ，p ，d x 和d y 都可以认为是常数，因此式( 卜5 ) 可简化 害= 历等+ 印雾一虬雾一“，考一“，考一k c。一。， 式( 卜6 ) 即为在x 方向和y 方向上存在浓度梯度时二维基本模型。 整二主堡堕! ! 墨：! 蔓! 式中c 一污染物的浓度，它是时间t 和空间位置x 、y 的函数； d x - - 纵向弥散系数；“一断面x 方向的流速分量： “y - - - - 断面y 方向的流速分量；d y y 坐标方向的弥散系数； k - 一污染物的衰减速度常数。 4 、三维水质模型 如果在x ，y ，2 三个方向上都存在浓度梯度，可以用类似方法推导出三维水 质模型( 见公式卜1 0 ) 。 害= e 等+ b 等+ e 窘叫塞飞考q 老一k c ( 1 - 1 0 ) 式中 ”：一z 坐标方向的流速分量；e ，b ，e ：一x 、y 、z 方向的湍流扩 散系数。 在三维模型中，因为不采用断面平均值，所以不出现弥散系数。式( 卜1 0 ： 就是三维基本模型。三维模型大量应用在大气质量的模拟和预测中，在深海 排放污水时也可以用三维模型进行水质预测。 1 2 有限单元法 1 3 1 有限单元法的发展和应用 在科学领域内，对于许多环境问题、力学问题和物理问题，人们已经得到 了它们应遵循的基本方程( 常微分方程或偏微分方程) 和相应的定解条件。 但能用解析方法求出精确解的只是少数方程性质比较简单，且几何形状相当 规则的问题。对于大多数问题，由于方程的某些特征的非线性，或者由于求 解区域的几何形状比较复杂，则很难得到解析解。随着电子计算机的飞速发 展和广泛应用，数值分析方法已成为求解复杂问题的主要工具。而有限单元 法的出现，是数值分析方法研究域内重大突破性的进展。 现代有限单元法第一个成功的尝试，是将刚架位移法推广应用于弹性力学 平面问题，这是t u r n e r ，c l o u g h 等人在分析飞机结构时于1 9 5 6 年得到的成 果。他们第一次给出了用三角形单元求得平面应力问题的正确解答。三十多 7 第一章概述( 1 页一8 页) 年来，有限单元法的应用已由弹性力学平面问题扩展到空间问题、板壳问题， 由静力平衡问题扩展到稳定问题、动力问题和波动问题。分析对象从弹性材 料扩展到塑料性和复合材料等，从固体力学扩展到流体力学、传热学、环境 水污染控制等连续介质力学领域。在工程分析中的作用已从分析和校核扩展 到优化设计并和计算机辅助设计技术相结合。 1 3 2 基本思想 有限单元法的基本思想是将连续的求解区域离散为一组有限个、且按一定 方式相互联结在一起的单元的组合体。由于单元能按不同的联结方式进行组 合，且单元本身大可以有不同的形状，因此可以模型化几何形状复杂的求解 域。有限单元法作为数值分析方法的另一重要特点是利用每一个单元内假设 的近似函数来分片地表示全求解域上待求的未知场函数。单元内的近似函数 通常由未知场函数或及其导数在单元的各个结点的数值和其插值函数来表 达。这样一来，一个问题的有限元分析中，未知场函数或及其导数在各个结 点上的数值就成为新的未知量( 也即自由度) ，从而使一个连续的无限自由度 问题变成离散的有限自由度问题。一经求解出这些未知量，就可以通过插值 函数计算出各个单元内场函数的近似值，从而得到整个求解领域上的近似解。 虽然随着单元数目的增加，也即单元尺寸的缩小，或者随单元的自由度的增 加及插值函数精度的提高，解的近似程度将不断改进。如果单元是满足收敛 要求的，近似解最后将收敛于精确解。 第二章有限单元法在二维水质模型中的应用 ( 9 页一1 3 页) 第二章有限单元法在二维水质模型 中的应用 2 1 二维水质模型的有限单元法 2 1 1 问曩的提出 污染物进入河流后，一方面随着水流运动而产生推流迁移，同时因分子扩 散，湍流扩散和污染物在水中的弥散等作用而分散，另一方面由于污染物将 在运动过程发生化学反应、生物降解等作用还会产生衰减和转化。污染物在 迁移过程中将与河水在纵向( 河流方向) 、横向( 河宽方向) 和竖向( 河深方 向) 三个方向进行。一般情况下，河宽比河流水深大得多，竖向混合过程相 对而言是在较短时间内即完成，因此可采用二维水质数学模型来描述浅水河 流中污染物浓度场分布。第一章中式( 卜1 5 ) 描述了当河流处于稳定的流动 状态、污染源连续均匀排放的条件下的二维水质数学模型的解析解。而现实 中，条件往往比较复杂，因此，本节就讨论利用有限单元法建立二维水质模 型的数值解。应用有限单元法与流带划分的有关理论，在缺乏河流速度与流 量分布资料的情况下，建立利用河道水下地形图( 即河床的标高) 进行复杂 河流水质数值计算方法，并以大源渡水库对湘江衡阳地区段水质影响预测为 例，研究了河流弯曲、分叉、江心有洲的浅水河段水质的数值计算方法。 2 1 2 流带 在描述流体的运动时，为了反映流场中的流速，分析流场中的流动，常常形 象化的引用流线、流带的概念。流线是在流场中，某一时刻各点的切线方向 与通过该点的流体质点的流速方向重合的空间曲线，为了问题研究的方便， 根据流体的某些特点( 如流量、河床形状等) 划分一定的区域和形状，形成 直观的流带。 9 第二章有限单元法在二维水质模型中的应用 ( 9 页一1 3 页) 2 1 3 正交曲线坐标系统 在给定的河流中，沿水流方向将河宽分成m 个流带，同时垂直水流方向， 将河段分为n 个子河段，构成m xr 1 个有限单元的平面网格系统如图2 1 。对 每一个有限单元来说，水质变化的原因包括：由纵向或横向水流的携带作用 造成的输入与输出；纵向及横向弥散作用形成的输入与输出；污染物的转化 与衰减；系统外部输入。根据这些关系，可以针对每一个有限单元写出质量 平衡方程，然后联立求解m n 个方程，就可以获得二维系统中的污染物分布a 二维系统中的横向水流分量的确定是非常困难的。如果在划分流带时， 使得每条流带的流量保持恒定，就可以忽略横向的水流交换。为了保持流带 内的流量恒定，流带的宽度就必然要随河流的形状不断变化，因而每个网格 的形状也是不断变化的。 g h 图2 - i 正交曲线坐标系统 2 2 b o d 二维水质模型 线 在研究有限单元法预测污染物在水体中的分布场中，本文仅以b o d 的浓度 分布进行讨论，其余水质指标预测的原理相同。 2 2 1 b o d 二维水质模型 从图2 1 中取出任意一个有限单元( 设为i 、j 单元) 。通过该单元的质量 平衡关系，可以写出一个差微分方程。 i o 第二章有硪单元法在二维水质模型中的应用( 9 页一j 3 页) 由水流输入、输出该单元的b o d 总量为： q ( 上- 1 ，一上f ) 由纵向弥散作用输入、输出该单元的b o d 总量为： d ( i ，- 1 、l ( 工j d ，j 一 ) 一d f ( i + 1 j ) ( 上p 一l + 1 ，j ) 由横向弥散作用输入、输出该单元的b o d 总量为： d ( “一l w ( 厶j l 一工f ) 一d ；o + i ) ( 工f 一工“+ i ) 在i ，j 单元内b o d 衰减量为v ；| i k 耐l # 由系统外输入的b o d 总量为w 式中 ( 2 - 1 ) ( 2 - 2 ) ( 2 - 3 ) d l i - l ，。m = d i 。自，i j ：；g 兰盟( 2 - 4 ) ( i - i j ) 0 d 。# 。+ ，j ) = d 。 ；。j ) ：；! i 业 a i j , ( i + l j ) z ，t 。，一- x ”2 。t ，一- x “j 糍 。“2 川，歹a i i j , ( i , i j + 1 ) 式中：q j 第j 个流带中的流量： l t j 第i j 个单元中的b o d 浓度； v t ，第i j 个单元的容积； k j 第i j 个单元中b o d 的衰减速度常数 d i j , , 1 单元i j 和k 1 间的弥散系数。 ( 2 - 5 ) ( 2 - 6 ) ( 2 7 ) 第二章有限单元法在二维水质模型中的应用 ( 9 页一1 3 页) a i j k l 单元i j 和k l 间的界面面积； 一x l j , l d 单元游相邻单元间的平均距离； y 叫横向相邻单元之间的平均距离。 对i j 单元写出b o d 质量平衡关系： 屹鲁钆( l t _ i j - l u ) + d ( f _ 1 ) g ( 工f - l 一g ) 一d f ( 。1 ) ( 上口一工f + i 。j ) ( 2 8 ) + d ( 1 j i 埘( l i , j - l 一口) 一d f ( j ，一i ) ( l 一l s , j + 1 ) 一v q k du l u 七w ； 在研究稳态问题时，面d l = o ，则 一g ( 三f i ，一岛) 一d o ，卜i ( ，卜i l 口) + d f 。( f ，+ i ) ( 三一j “) 一d ( ，一l m ( 三f - l j 一三f ) ( 2 9 ) + d f ，( f + 1 j ) ( 工f i t + l j ) + k d f 工p = 嘭 如果将所有有限单元中的b o d 值写成n l n 维向量：l = ( l ，l 。：l 。l 。l l i j l n m ) ；将所有系统外输入出写成一个m n 维向量： w “= ( w 。lw 1 2 l 、纰w l 。w 。l w 毒w 二) 1 ，对整个研究的河段可以写出矩 阵方程： g = 谛。( 2 一1 0 ) g 是一个i l l n 阶矩阵，称为二维河流的b o d 变换矩阵。根据式( 2 - 9 ) ，可以 写出计算矩阵g 的各个元素g k e ( k = l ，m ，l = i ，n ) 的一般形式： 对l = k ，g k i = q j + d ，( i 十i ) i j + d i j ( ) ( i ，j + 1 ) + d ( i 1 巾，i j + d i j ，( i + 1 ，j ) + v l j k 刚( 2 - 1 1 ) 对l = k + l ，鼬= d i ，( i ，j + d ( 2 - 1 2 ) 对l = k 一1 ，g k l = d ( i ，j 1 ) ，i j( 2 - 1 3 ) 对l = k + m ，g k l = 一n i j ，( i + l ，j )( 2 - 1 4 ) 对l = k m ，g k l = q j d ( i 一1 ，j ) ，i j( 2 一1 5 ) 所有其余元素g k l - 0 。 第二章有限单元法在二维水质模型中的应用( 9 页一1 3 页) 由以上各式可以看出，矩阵g 的元素是河流流带的流量、弥散系数，各单 元的几何尺寸及b o d 衰减速度常数的函数，它们可以通过各种测量和计算方 法得到，下面的章节将会加以介绍。因此，当系统外部b o d 输入为已知时， 就可以用以下求得复杂河流的b o d 分布： = g “霄。工厂( 2 1 6 ) 式中，g 。1 称为二维河流的b o d 响应矩阵。 2 2 2 流带的划分l 假定河流的计算流量为q ，河宽为b ，横断面积为a ，断面形状如图2 2 所示。河流断面上任一单位宽度上的流量可以用以下式计算： 第三章二维水质模型的数值计算 3 1 数值分析 3 1 1 数值分析的对象与特点t 数值分析是研究各种数学问题求解的数值计算方法。在计算机成为数值计 算的主要工具以后，则要求研究适合于计算机的数值计算方法。为了更具体 地说明数值分析的的研究对象，我们考察用计算机解决科学计算问题时经历 的几个过程： 匡西i i 圊一匿i i i i 国一匿垂区圈一匿j i i i 圃一匡夏固唾重蚕i 圜 数值分析也称计算方法，但不应片面地理解为数值方法的简单罗列和堆 积，它是一种与计算机使用密切结合的实用性很强的计算方法。由实际问题 的提出到上机求得问题解答的整个过程都可以看作是应用数学的任务。如数 值分析就是研究用计算机解决数学问题的数值方法及其理论，它的内容包括 函数的数值逼近、数值微分与数值积分，非线性方程数值解、数值代数、常 微分和偏微分数值解等。 3 1 2 数值分析的特点 数值分析的特点概括起来主要有以下四点： 第一、应用计算机程序和语言。要根据计算机特点提供实际可行的有效 算法，即算法只能包括加、减、乘、除运算和逻辑运算，是计算机能直接处 理的。 第二、有可靠的理论分析。能注意接近并达到精度要求，对近似算法要 保证收敛性和数值稳定性，还要对误差进行分析，这都是建立在相应数字理 论的基础上。 第三、要有好的计算复杂性。时间复杂性好，能节省时间；空间复杂性 好能节省存储量，这也是建立算法要研究的问题，它关系到算法能否在计算 1 4 第三章二维水质模型的数值计算 ( 1 4 页一2 2 页) 机上实现。 第四、要有数值实验。即任何一个算法除了从理论上要满足以上条件之 外，还要能够解决实际中存在的问题，通过数值试验证明是行之有效的a 3 2 数值分析的主要方法 随着科学技术的发展，数值分析在科学领域中应用越来越广，目前主要 有插值法，函数逼近与计算，曲线拟合的最小二乘法，数值积分与数值微分， 解线性方程组的迭代法等等，本节重点介绍插值法。 3 2 1 插值法 许多实际问题都用函数y = f ( x ) 来表示某种内在规律的数量关系，其中相 当一部分函数是通过实验或观测得到的。虽然f ( x ) 在某个区间 a ，b 上是存在 的，有的还是连续的，但却只能给出 a ，b 上一系列点x - 的函数值 y i = f ( x ) ( i = o ，1 - n ) ，这只是一张函数表，有的函数虽有解析表达式，但由于 计算复杂，使用不方便，通常也造一个函数表，如大家熟悉的三角函数表、 对数表、平方根和立方根表等等，为了研究函数的变化规律，往往需要求出 不在表上的函数值。因此，我们希望根据给定的函数表做一个既能反映函数 e ( x ) 的特性，又便于计算的简单函数p ( x ) ，用p ( x ) 近似f ( x ) ，通常选一类较 简单的函数如代数多项式或分段代数多项式作为p ( x ) ，并使p ( x - ) = g ( x t ) 对 i = o ，1 n 成立，这样确定的p ( x ) 就是我们希望得到的插值函数。 设函数y = f ( x ) 在区间 a ，b 上有定义，且已知在点 口s z o 而 x 。 b 上的值y o ，y 。，y 。，若存在一简单函数p ( x ) ， 使 p ( x i ) = y 。( i = o ，l ，n ) ( 3 - 1 ) 成立，就称p ( x ) 为f ( x ) 的插值函数，点x 。，x i ，“。x n 称为插值节点，包含插值 节点的区间 a ，b 称为插值区间，求插值函数p ( x ) 的方法称为插值法，若p ( x ) 是次数不超过r l 的代数多项式，即 p ( x ) = a o + a i x + + a i x n 其中a 。为实数，就称p ( x ) 为插值多项式，相应的插值法称为多项式插值。若 p ( x ) 为分段的多项式，就是分段插值。若p ( x ) 为三角多项式，就称三角插值。 第三章二维水质模型的数值计算( 1 4 页 2 2 页) 本章只讨论多项式插值与分段插值。 从几何上看，插值法就是求曲线y = p ( x ) ，使其通过给定的n + 1 个点( x i ， y 。) ，i = o ，1 ，n ，并用它近似已知曲线= f ( x ) ，见图3 1 插值法是一种古老的数学方法，它来自生产实践，早在千多年前，我国 科学家在研究历法上就应用了线性插值与二次插值。但它的基本理论和结果 却是微积分产生以后才逐步完善的。其应用也日益增多，特别是在电子计算 机广泛使用以后，由于航空、造船、精密机械加工等实际问题的需要，使插 值法在实践上或理论上显得更加重要，并得到进一步发展，尤其是近十几年 发展起来的样条( s p l i n e ) 插值，l a g r a n g e 插值法更获得了广泛的应用。 ， 图3 - 1 插值法图示 3 2 2 拉格朗日播值 l 、线性插值与抛物插值 假定已知( x ) 的端点处的函数值y 。= f ( x 。) ，y 。= f ( x 。) ，要求线性插值 多项式l n ( x ) ，使它满足： l t ( x - ) = h ， l i ( x ) = y k + i 。 y = l t ( x ) 的几何意义就是通过两点( x 。y 。) 与( x 。，y 。) 的直线，如图( 五) 所示， l t ( x ) 的表达式可由几何意义直接给出： l i ( x ) = y i + 丝盥( x - x , ) ( 点斜式) 。 ( 3 2 ) 工i “一以 1 6 第三章二维水质模型的数值计算 ( 1 4 页- 2 2 页) l i ( 石) ：立丝二兰儿+ 三冬y ( 两点式) 。 ( 3 2 1 ) x i + 1 一x ix k + l x i c 由两占式看出，l i ( x ) 是由两个线性函数 ( x ) ：盥，厶+ l ( 工) = 竺 ( 3 2 2 ) x i x 女+ 1j k + l 一 i 的线性组合得到，其系数分别为y 。及y “，即： l t ( x ) = y k l 。( x ) + y i l h 。( x ) ( 3 3 ) 显然，l 。( x ) 及l 。( x ) 也是线性插值多项式，在节点x k 及x 。上满足条件： k ( x k ) = 1 ，l k ( xk + 1 ) = 0 ， l k + l ( x - ) = 0 ，l - ( x n , - i ) = 1 我们称函数l 。( x ) 及l 。( x ) 为线性插值基函数。 下面讨论n = 2 的情况，假定插值节点为x 。x 。，孙，要求二次插值多 项式l 2 ( x ) ，使它满足 l 。( x j ) = y ，( j ；k 一1 ，k ，k + 1 ) 我们知道y = l ：( x ) 在几何上就是通过三点( x y 。) 、( x 。y - ) 、( x 。y 。) 的 抛物线，为了求出l 2 ( x ) 的表达式，可采用基函数方法，此时基函数l 。( x ) ， l k ( x ) 及l 。( x ) 是二次函数，且在节点上满足条件： l k 1 ( x k ) = o( j = k ，k + 1 ) l k ( x k ) = 0( j = k l ，k + 1 ) ； ( 3 6 ) l k “( x k ) = 0( j = k 一1 ，k ) 满足条件( 3 6 ) 的插值函数是很容易求出的，例如求l 。( x ) ，因它有两个零 点x k 及x 。故可表为 l 。一。( x ) = a ( x x k ) ( x - - x 。) ， 其中a 为待定系数，可由条件l 。( x 。) = 1 定出： a ：j i 一， ( x k l x k ) ( x x k + 1 ) 张l k - l ( x ) 2 老岩善，【x k i x k 儿x k l x k “) 同理可得2 芒岩皋悬，【x k x k 一1j l x k x k + lj = l ) l = = p芯肌 【l 飞 一忆、 第三章二维水质模型的数值计算 ( 1 4 页也2 页) l “( x ) = 瓦( _ x - - 瓦x k _ 1 丽) ( x 而- - x k ) ， 二次插值基函数l 。( x ) ，l 。( x ) 分布在区间 x x m 上。 利用二次插值基函数l 。( x ) 、l 。( x ) 、k ( x ) , - e r 目p 得到二次插值多项式： l ：( x ) = y k ，l b 。( x ) + y k l k ( x ) + y h i l k + 。( x ) ( 3 - 7 ) 显然，它满足条件l ：( x 。) = y 。( j = k 一1 ，k ，k + 1 ) ，将上面求得的 l 。( x ) ，l 。( x ) ，l 。( x ) 代入( 3 - 7 ) ，得 芒鲁等 + v 堡二毡! ! 兰二塾! ! “( x k - - x ) ( 】【l ( 一x ) + y ( x xk - 1 ) ( x x k ) k + 1 ( x k + l x k 1 ) ( x k + l x k ) 拉格朗日插值多项式 上面对n = l 及n = 2 的情况，得到了一次与二次插值多项式l ( x ) 及l 2 ( x ) ， 它们分别用( 3 - 5 ) 与( 3 - 7 ) 表示，这种用插值基函数表示的方法容易推广到一 般情形，下面讨论通过n + 1 个节点x o ( x 的n 次插值多项式l 。( x ) ，假定 它满足条件 l ( x j ) = y j( j = o ，l ，n ) 。 ( 3 - 8 ) 为了构造1 ( x ) ，我们先定义n 次插值基函数。 定义l 若n 次多项式l j ( x ) ( j = o ，l ，n ) 在n + 1 个节点x o x x 。上满 足条件 l j ( 妒 甚嚣( j ，k - o ，1 ，n ) ( 3 - 9 ) 就称这n + 1 个n 次多项式l 0 ( x ) ，l ( x ) ，l m ( x ) 为节点x 。，x l ”，x n 上的n 次插 值基函数。 当n = 1 及n = 2 时的情况前面已经讨论，用类似的推导方法，可得到n 次 插值基函数为 l ；( x ) ：j 兰鱼堕竺d 坠垡出竺量l 1、7 ( x k x o ) ( x k x k 1 ) ( x k x k + i ) ( x k x n ) ( k = o ，l ，n )( 3 - 1 0 ) 显然它满足条件( 3 9 ) ，于是，满足条件( 3 8 ) 的插值多项式l n ( x ) 可表示 为 l 。( x ) = y k l t ( x ) k w o 由l 。( x ) 的定义，知 l 。( x ) = y k l k ( x ) 2 y j ( j 2 0 ，l ，n ) k 0 形如( 3 一1 1 ) 的插值多项式l 。( x ) 称为拉格朗日 ( 3 5 ) 与( 3 7 ) 是n = 1 和n = 2 的特殊情形。 若引入记号 w 。+ ( x ) ：( x x 。) ( x x ，) ( x x ) 容易求得 ( 3 1 1 ) ( l a g r a n g e ) 插值多项式，而 ( 3 1 2 ) w n + i ( x k ) = ( x k x 。) ( x i - x h ) ( x k - - xk + i ) ( x k - - x 。) 于是公式( 2 1 1 ) 可改写成 5 耋y k 高攀( 3 - 1 3 ) 注意，n 次插值多项式l ( x ) 通常是次数为1 1 的多项式，特殊情况次数可 能小于n ，例如，通过三点( x 。，y 0 ) 、( x 。y 。) 、( x 2 ，y ：) 的二次插值多项式l ：( x ) ， 如果三点共线，则y = l ，( x ) 就是一直线，而不是抛物线，这时l 2 ( x ) 是一次式。 3 2 3 三次样条插值( s p l i n e ) 上面讨论的分段低次插值函数都有一致收敛性，但光滑性较差，对于像 高速飞机的机翼，船体样头等值线往往要求有二阶光滑度。即有二阶连续导 数，早期工程师制图时，把富有弹性的细长木条( 所谓样条) 用压铁固定在 样点上，在其他地方让自由弯曲，然后画下长条的曲线，称为样条曲线，它 实际上是由分段三次曲线并接而成，在连接2 点即样点地要求二阶导数连续， 从数学上加以概括就得到数学样条一概念，最常用的是三次样条函数。其定 义如下： 若函数s ( x ) c 2 a ，b l ，且在每个小区间 x 。，x 。 上是三次多项式，其中 a = x 。 x x 。= b 是给定节点，则称s ( x ) 是节点x 。，x i ，“，x 上的三次样条函数。 若在节点x ，上给定函数值y l = f ( x j ) ( j = 0 ，l ，n ) ，并成立 s ( x j ) = y j ，( j = o ，= l ，n ) ，( 3 - 1 4 ) 则称( x ) 为三次样条插值函数。 对于推导三次样条插值函数比较复杂，但在m a t l a b 语言中，s p l i n e 作为 其内部函数命令已经定义。故使三次样条插值函数的应用更加方便和实用。 1 9 第兰章二维水质模型的数值计算 ( 1 4 页一2 2 页) 3 3b o d 二维水质模型的数值计算 在第二章中，根据有限单元法和b o d 物质的质量平衡原理已经推导出了 b o d 的二维水质模型为： 式( 2 1 0 )= g 。霄。 要准确预测出b o d 的浓度，我们就要求g 。才能得到二维河流的b o d 响应 矩阵6 - 1 。其中g 。是河流流带的流量、弥散系数、各单元的几何尺寸及b o d 的 衰减速度常数的函数。 3 3 1 已知量 l 、利用1 ：2 0 0 0 的河道地形图可以读出要研究的各断面的局部深度h t