伊利湖和安大略湖污染问题的数学模型.pdf

伊利湖和安大略湖污染问题的数学模型 裘丛欣 梁曼 2004 4 23 2004 4 25 内内内容容容摘摘摘要要要 研究建立伊利湖和安大略湖的污染物迁移模型 用差分法对简单的零维维湖泊模型进行了较详 细地分析和求解 计算了污染物浓度降到特定值的具体时间 并对一维湖泊污染物扩散模型进行了 分析 最后对其他一些研究手段稍作讨论 关键词 污染物 扩散 湖泊 1有有有关关关全全全局局局假假假设设设的的的一一一些些些说说说明明明 这道题目已经对模型作了如下假设 全部安大略湖的水来自伊利湖 往伊利湖里注入的湖 水未受污染 并且 题目还提供了一些数值上的数据 伊利湖每年有 38 的湖水被更新 而对应 安大略湖是每年 13 题目对初始条件的限定是 最初安大略湖湖水中有 2500 单位的污染物 而 一年后污染物总数是 3150 单位 仅有的污染源 铝厂 每年向安大略湖注入 25 单位的污染物 在以上假设的基础上 我们对题目添加如下假设 两湖的水量都是稳定不变的 且湖水的更 新是等速的 污染物的浓度非常小 其所导致水流的密度变化可以忽略 污染物溶于水后与水的 分压相同1 2零零零维维维湖湖湖模模模型型型 作为最简单的情况 我们先来讨论零维湖的模型 在这个模型中 某个湖的污染状况可以用一 个标量来表示 这说明 我们对于湖的模型已经做了如下加强假设 污染物的扩散过程是量子化 的 即在某一个特定时刻 要么它完全未扩散 要么它已经扩散均匀 这个假设也等同于 湖水中污 染物的扩散速度是无穷大 亦即 污染物可以在瞬间在整个湖体中混合均匀 在这个条件下 我们便可以列出湖水浓度变化的差分方程 a n 1 0 38 a n 1 b n 0 38 a n 1 1 0 13 b n 1 25 1 其中 a n 和 b n 分别代表伊利湖 A湖 和安大略湖 B湖 的总污染量 而 n 代表经过的 年数 这个模型的意思是 在每一年的年末 A湖一次性地把它 38 的水 亦即 38 的污染物 排 出 并注入干净的新水 而B湖首先把它 13 的水 亦即 13 的污染物 排出 随即接收A湖刚刚 1也许应该换一种对于液体更适合的描述 大体的意思就是 污染物溶于水后的质量 体积 压强等与水的同类数据满足某 种关系 以使菲克 Fick 扩散定律成立 1 排出的水以及来自工厂的 25 单位污染物 在以上过程结束后 A湖和B湖立即把自己的湖水混合均 匀 通过方程 1 我们很容易推出这个离散动力系统中 a n 和 b n 的二阶递推关系 a n 1 0 38 a n 1 b n 2 0 38 0 13 b n 1 1 0 13 1 0 38 b n 2 0 38 25 以及 a n 和 b n 的通解 a n 0 62 na 0 b n 0 87 nb 0 0 87n 1h1 0 87n 2h2 0 87hn 1 hn 其中 hn 0 62n 1 b 1 0 87b 0 0 62n 2 1 9 5 并且 题目已提供了关于 b n 的两个初始条件 即 b 0 2500 b 1 3150 把这两个数值代入对应 n 1 时关于 b n 的方程 即可得出 a 0 2500 随即 我们便可以通 过迭代关系求出对应任何 n 时 a n 和 b n 的数值 变化情况如见图 1 和图 2 5101520 500 1000 1500 2000 2500 3000 图 1 前二十年 a n 和 b n 的变化 在这两张图中 横坐标代表时间 年数 n 而纵坐标代表污染物的量 上方的曲线是 b n 而下方的曲线是 a n 图 1 对应的是前 20 年 a n 和 b n 的变化情况 而图 2 对应的是 前一百年的 在图中我们可以看出 a n 一直是单调下降的 并且其实 a n 的所有数据点全部 落在一条指数衰减的曲线上 而 b n 在第一年和第二年上升 随后便单调下降 我们还可以看 出 b n 的衰减远比 a n 滞后很多 而当年数 n 时 a n 0 而 a n 192 3 即方程 0 13 a n 25 的解 即 从长尺度的时间上来讲 A湖和B湖的污染量都稳定在确切的数值上 并 且 污染远比目前的状况好 并且 b n 在第 34 年降到了其初始值 2500 的 10 以下 那年的污染 量 b 34 245 953 2 20406080100 100 200 300 400 500 600 700 图 2 前一百年 a n 和 b n 的变化 3一一一维维维湖湖湖模模模型型型 3 1相相相关关关地地地理理理信信信息息息及及及零零零维维维湖湖湖模模模型型型的的的局局局限限限性性性 在对湖体模型进行修正以前 我们先来分析一下伊利湖和安大略湖的地理状况 伊利湖呈西西南 东东北向 东西长 388 千米 最宽 92 千米 面积 2 57 万平方千米 在五大 湖中居第四位 湖面海拔 174 米 比安大略湖高 99 米 平均水深 18 米 最大深度 64 米 在五大湖 中最浅 蓄水量 455 立方千米 西经底特律河 圣克莱尔湖 圣克莱尔河接纳苏必利尔湖 密歇根 湖 休伦湖的湖水 东经尼亚加拉河注入安大略湖 安大略湖是北美洲五大湖最东和最小的一个 大致成椭圆形 主轴线东西长 311 公里 最宽处 85 公里 水面约 19554 平方公里 平均深度 86 米 最深 244 米 蓄水量 1688 立方公里 湖水通过 西端著名的尼亚加拉大瀑布从伊利湖流入 通过东端的圣劳伦斯河流出 在零维湖的模型中 我们事实上已经做了如下假设 污染物可以在瞬间在整个湖体中混合均 匀 这个假设显然是很不合理的 因为伊利湖和安大略湖都是东西宽三四百公里的巨大水体 混合均 匀是相当困难的事情 因此 污染物的扩散过程在这个数学模型中是一个必须考虑的问题 3 2一一一维维维湖湖湖模模模型型型的的的必必必要要要假假假设设设以以以及及及对对对假假假设设设合合合理理理性性性的的的分分分析析析 污染物的扩散应该大体上有以下几种方式 对流 弥散 分子扩散 吸附 化学反 应 如果对应单一孔隙空间并在微观水平上考虑 便只有对流迁徙和分子扩散两个组成部分 在前面我们已经假设 菲克 Fick 扩散定律 J D dN dz 2 3 对于这个水体中污染物的输运过程是适用的2 在这个公式中 D 为 与 的互扩散系数3 根据 这个扩散定律 我们便可以列出有关时间相应的微分方程 N t D 52N 0 3 其中 为流体即水体的密度 在有些书上也把 D 定义为扩散系数或是扩散率 由于在这个方程 中只出现了 N 关于 t 的一阶偏微商 在确定初始条件 或者对应任何时刻能够描述下面时刻发展状 态的参量 时 我们只需要描述出 N 即污染物浓度 对应空间坐标的分布即可 在已有假设的基础上 我们对模型添加以下假设 湖水的流动为定常流动 steady fl ow 即速度 v 仅为坐标 x y z 的函数 与时间 t 无关 湖水的环量 circulation 或对应的涡 度vorticity 为零 在单一孔隙空间并在微观水平上考虑 不考虑污染物的沉淀 从假设 可以推出不考虑对流 而从假设 可以推出不考虑吸附和化学反应 因而 污染物的 扩散仅由弥散造成 从以上关于湖体地理情况的分析 我们注意到 两湖有一些共同的 但对于 湖 这个大的客 体很特殊的性质 两湖面积大体相当 并且都呈长条状 两湖长度均长度 400 公里左右 而宽度只 有 90 公里左右 并且湖水都是从其一各端点注入 从另一个端点流出 因此 我们定义有关湖体 的直角坐标系 X 轴沿湖体的长轴方向 Y 轴沿湖体的短轴方向 而 Z 轴对应垂直方向 由假设 定常流动和假设 环量为零 我们可以推出 湖水在 Z 方向是干净的层流 则 由假 设 不考虑污染物的沉淀 我们便可以推出 在这个模型中 我们可以不考虑污染物在 Z 方向的变 化 3 3进进进一一一步步步简简简化化化 下面 我们以伊利湖为例 对湖水流动的性质作一些数量级上的估计 并且希望由此 遂与模 型作一些进一步的简化 伊利湖最大的入口 我们假设其为唯一的入口 在温莎及底特律附近 而唯一的出口在布法罗 及尼亚加拉瀑布附近 两者相距四百公里 又由已知条件 伊利湖每年有 38 的湖水被更新 即 每一个水分子在伊利湖中平均停留 2 6 年 由此 我们可以计算出湖水定向流动的平均速率为 v 400km 2 6yr 4 8 10 3m s 由此 我们便可以推断这个体系的雷诺数 Reynolds number R 我们知道 当小雷诺数的情况 下 体系与泊肃叶公式符合得较好 而当 R 4000 时 体系将会出现湍流 把 v 的值代入公式 R d v 水分子的运动粘性 1 4 10 6m2 s 求得当 d 流体微元与边界距离 大于 1m 时 会出现湍流 有关流体微元与边界的关系 也可以通过复壁边界层的理论加以理解 边界层的厚度 p 1 2 10 3m 2这件事情的合理性可能是非常值得仔细考虑 因为菲克定律描述的是 气体 分子热运动引起的浓度趋同现象 并且 在 菲克定律中 质量的流动是重要的 而在我们的模型中 污染物的浓度引起水流的密度变化是可以忽略的 但是 这种正比于 梯度的变化在诸如热传导 扩散等很多物理现象上是普适的 并且 我们目前能够确定 用于处理海水中盐分的输运这类问 题 所以 如果我们把菲克定律看作是一种在相对宽泛一些的范围内的说法 这种假设应该是合理的 并且这种推断也是正确 的 3我们可以假设 如果把本体的水去掉 那么在有污染物的地方目前只剩下污染物 而没有污染物的地方 目前剩下的是与 污染物分压相同的很少量的水 当然这个说法很蹊跷 我们也不妨把下标 和 去掉 而把 D 理解为一种从更宽泛意义角度 上对应浓度梯度的系数 叫做扩散系数 4 综上 边界对于水体的影响只在离边界很近的一米以内的区域起作用 而对于 Y 方向来说 湖 体的宽度为 90 公里左右 远大于一米的边界影响宽度 而 Y 方向上虽然有湍流 但湍流是小范围的 对应海水 湍流的大小应为厘米到千米的量级 对于湖水 应该也可以作类似估计 而 Y 方向 上各点在大范围上来看是平权的 且无法分辨的 因而我们只能推断出忽略污染物在 Y 方向上变化 的影响 对于安大略湖 数据并没有量级上的差别 因此 至此 我们便完成了湖体的一维简化模型 在这个模型中 我们只需要考虑污染物在 X 方向上的分布 3 4非非非扩扩扩散散散极极极端端端简简简化化化模模模型型型 作为一种最简单的情形 也作为预先前完全扩散的零维湖模型的对比 我们可以考虑另一种极 端情况 即 污染物完全不扩散 即 其完全随波逐流 在一维的简化模型中 湖水代谢必须遵循先 进先出 后进后出的原则 因而A湖在不到三年的时间便会完全干净 而B湖的污染物代谢可以分解 成完全独立的三个过程 附加的全局假设 分别代表B湖原有水的随波平移 A湖来水的随波平 移 和铝厂污染物的稳定排出 这样 在十一年 n 1 13 38 13 之后 B湖的污染物代谢也进 入稳定状态 如图 1 图 3 非扩散极端简化模型 当然 实际情况是 污染物在以一个有限大的速度扩散 因而 在一维湖的假设下 湖水污染 的数学模型也应该介于这两种模型之间 即 污染物含量降低的速度也应处于两者之间 5 3 5对对对于于于扩扩扩散散散的的的量量量级级级估估估计计计 扩散系数 D 应该是一个实验量 它与若干因素有关 由于扩散从本质上说是由无规则的热运动 造成的 因此扩散系数与温度 T 有关 又 很显然 扩散与污染物微元的质量 m 有关 显然 若 m 而污染物又不沉淀 则其将会完全随波逐流 甚至 若我们判断污染物究竟是分子还是 尘埃颗粒 处理的方法也会完全不同 分子为本征的运动 可能可以通过定义温度的 v2或者通过 麦克斯韦分布推出的 vp和 v来描述 显然速率是与扩散更密切相关的 而尘埃颗粒进行的是布朗 运动 可能需要用爱因斯坦扩散方程 s2 2 kBT 6 a t 来描述 4 炼铝厂排放的污染物是多种多样的 质量的差别也很大 包括各种废水 废气 粉尘 石油类 等等 因此 我们很难找到一个对于炼铝厂排放污染物普遍适用的扩散系数 而对于那些原先在湖中 的污染物 它们的污染源千差万别 它们的总类也会更加繁多 因而我们更难估计它们的扩散系数 我们可以通过盐在海水中扩散对此作一个大略的估计 对于盐 D 1 5 10 8kg m s 则 D 1 5 10 11m2 s 即 只靠单个盐粒子的运动 它扩散一米大约需要 1011s 的时间 合大约三 千年 显然 这个模型是错误的 因为盐的扩散主要并不是靠分子热运动而是靠水体的湍流 湍流引 起的扩散可以用涡动扩散系数 交换系数 来描述 但因为湍流的尺度弥散非常大 从厘米两极到千 米量级 涡动系数的弥散也相当大 大约是分子扩散系数的 102到 1010倍 因而很难通过理论 的方法对此进行准确估计 我们可以试图把以上理论用于分析湖水污染的模型的分析 如果考虑湍流 盐分子 以及假设 移植于我们需要考虑的铝厂排放污染物和湖中其它污染物 扩散的速度将介于 10 9m s 到 10 1m s 之间 而如果湖水需要用一年的时间混合均匀 需要的扩散速度大约是 10 2m s 量级 400km 1yr 介于这两者之间 因此 这个理论很难作出定性的判断 由此 我们便知道以上的假设 湖水的流动为定常流动 是不合理的 因为在定常流动的湖水 中与在有湍流的湖水中 污染物的扩散速度要差好几个数量级 这个假设必须舍去 但是在 Z 方向 上 由于存在着密度 温度和其他条件的差别5 湍流的形成要比水平方向上难一些 因此 Z 方向上 的情况可能将是稳定的分层结构和比较慢的热扩散6 因而 我们的一维湖简化模型应该是这样的 X 方向主要为湍流引起的扩散 Y 方向忽略 Z 方向很慢的热扩散及稳定的分层结构 暂且不考虑对流 并且把 Z 方向省略 3 6一一一个个个尝尝尝试试试性性性的的的扩扩扩散散散模模模型型型 由于湖水的环量为零 湖水中污染物的运动必然可以分解成随湖水的运动和其本身的扩散 无 论是热运动还是湍流造成的 两部分 在一维湖的假设下 湖水的运动必然遵循先进先出 后进后出 的原则 又因为其更新为等速 全局补充条件 所以湖水的运动完全与 非扩散极端简化模型 中的情况相同 这样 下面单独讨论在静水中的扩散模型 由于浓度 N 在空间坐标上只与 X 方向有关 方程 3 可以简化为有关 N x t 的偏微分方 程 并且对于湖A 它可以被描述为用隔板隔开的被污染的水和净水在抽取隔板后的扩散现象7 因 4有关速率的三个公式 以及爱因斯坦扩散方程均涉及玻耳兹曼常数 kB 而 kB 是通过摩尔气体常数 R 推出的 因而 我 们不能确信这个关于气体的公式是否能够用于流体 5对于湖这样浅的体系是否如此可能还需要再议 6在 Z 方向上 对流应该是不能忽略的 知识在这一部分 我们暂且不考虑这个问题 7这里 我们的起始时间是 n 0 并且我们假设在这个时候 除了铝厂 所有的污染源刚刚全部停工时 湖水的污染是均 匀的 6 而 我们很容易定出其边界条件如下 N t D 2N x2 0 N x N0 N x 0 N x x 0 N x x 0 N t 0 N0 x 0 4 这个方程用分离变量法求解似乎很困难 我们没有解出 但是我们可以通过差分的方法对这 个过程作一个简要的分析 首先 我们只做一个 1 2n 的表格 并且把右边一半涂黑 表示被污 染 然后作反复的迭代 迭代过程是 把每一个方格分成两部分 并把它的某一部分和另一各方各 相邻的部分作算术平均 这个过程其实就是菲克扩散方程 2 的差分表示 而格子的大小和迭代的 次数的比例关系确定扩散系数 D 事实上 如果选择比较小的格子和比较多的迭代次数 我们就可 以得到比较接近于方程 4 的折线 当然 如果格子的尺度固定 比较多的迭代次数则表示比较长 的时间 如图 4 图 4 差分菲克扩散定律的迭代图示 我们用一个简单的C程序演示了这种迭代过程 下面分别是迭代 9 次 29 次 59 次 99 次的 浓度分布图 以及一个它中间部分的细微结构 但是 这种迭代算法有一个明显的缺陷 在这种算法中 浓度的传播有一个绝对的边界 零或 非零 而事实上 这种边界是不存在的 考虑传播速度 v vc 不考虑远处的那个边界 只是 到了非常远的地方就非常弱 而长时间和短时间的浓度分布是相似的 还好 当迭代次数非常多时 差分便于微分的结果非常类似 只是 对于差分结果 只有中间很小的一段的数值 它们偏离了 0 或 者 1 比较多 对于我们来说是有意义的 7 40 202040 0 2 0 4 0 6 0 8 1 图 5 迭代 9 次 29 次 59 次 99 次的浓度分布图 15 10 5510 0 2 0 4 0 6 0 8 1 图 6 上图中间部分的精细结构 如果考虑湖水的流动 那么仅仅是在以上模型上加一个沿 X 方向的平动 而平动速度主要由涡 动扩散系数决定 由于使用非试验的办法 我们所能确定的涡动扩散系数在数值上有八个数量级的差 别 而其中有的数量级与湖体的尺度是相当的 因而 在这种情况下 下一步的计算必须依赖其他数 据 3 7用用用该该该扩扩扩散散散模模模型型型解解解决决决安安安大大大略略略湖湖湖问问问题题题的的的大大大略略略构构构想想想 安大略湖的问题比伊利湖复杂很多 若不考虑铝厂的影响 由于上言补充的全局假设 每一 个污染物粒子都是相互独立的 污染物粒子之间的影响可以忽略 因而可以完全独立地考虑它们的扩 散 新从伊利湖注入安大略湖的水和安大略湖原有的被污染的水的确有一个边界 作为扩散起始的 边界 并且这个边界在与那种极端简化模型完全相同的方式匀速地向右推进 这个模型与上言的 模型的差别在于 初始条件 N t 0对应 x 0是两个不同的非零值 虽然初始污染物都是 8 2500 单位 但是由于安大略湖的蓄水量比伊利湖大很多 因而它们污染物的浓度有很大差别 但 这个并不关键 只要改变零点的选择和比例就可以了 其实不用 这样也很好解 另一个关键的问 题是 对应 N x 的边界条件异常复杂 因为另一个边界的扩散 并且由于尼亚加拉瀑布的单向 性 那个边界的确是单向地对这个边界造成影响 解决这个问题还需要进一步的讨论 而如果考虑 铝厂 铝厂的位置所起的作用也是至关重要的 问题将更加复杂 只有一种情况可以使的边界 N x 对这个边界的影响消失 那就是那个边界的弥散足够慢 在涡动扩散系数足够大的取值范围内 这个条件很可能被取到 以至于无法对其造成影响 但其 实 若这种情况 便不如采用那种非扩散的极端简化模型 因为这两种模型由于这样的条件变得没有 本质的差别了 顺便提一点 在建立模型的过程中我们发现AdobePhotoshop或许是一个挺有意思的软件 它虽 然适用于图像处理 但是若在其中使用灰度图像 利用软件的模糊滤镜 高斯模糊 可以很容易的模 仿不同浓度污染水体相互融合的过程 而且通过使用工具栏中的拾色器功能方便的获取图像中任一点 的灰度值 可以认为它代表了该点的污染物浓度 图 7 25 不同的高斯模糊程度产生的混合效果图示 3 8另另另外外外一一一个个个关关关于于于扩扩扩散散散的的的解解解决决决办办办法法法 我们也曾经希望通过考虑单个分子的热运动来讨论污染物粒子的扩散 粒子的速度并不是在 0 上平权的 因而粒子扩散造成的相对原始位置的偏离并不是高斯分布 我们假设粒子只能取一 个速度 比如可以取由麦克斯韦分布推出的 v q 8kBT M 这样 偏离便呈现二项分布 我们尝试 做了 25 个单位粒子 铝厂 经过八次碰撞的偏离 当粒子数或碰撞次数增大时 计算量迅速增大 这使得