期末作业 Final Assignment 2023年5月28日 数学建模课程 May 28th,2023 是排放点源(0,0，)提供的浓度： =8-（品+】 是全反射像源(0,0，-h)提供的浓度。 计算地面灰尘浓度即令之=0，得到高架连续点源扩散的地表浓度模型： Q =8.（+）】 2 若要计算地面最大浓度，即求： Q mg C(,0.0)-Udy.exp2) (2.9) 为简化计算不妨令0，/o:为常数S,在式(2.8)中令： dcc,0,0=0 (2.10) da 求解式2.9)得到： 20S Cnar-URe (2.11) S2.3拉格朗日扩散模型 拉格朗日大气扩散模型可以细分为两类：拉格朗日粒子模型与拉格朗日烟团模型。 2.3.1拉格朗日粒子模型 该模型是基于湍流统计理论得到的扩散模型，可以细致追踪到每一个粒子在大气中的运 动状态，计算其时间、空间概率分布情况，以估算出灰尘浓度变化同。 粒子模型将粒子速度视为风场平均速度与湍流速度之和： V=V+V (2.12) 而湍流速度矢量满足： V'(t+△)=V'()R(△)+p (2.13) 其中，t为时间：△t为时间步长；R为拉格朗日自相关系数：P为蒙特卡洛分量。 在实际计算中，粒子模型将带预测区域划分为若干个网络，通过记录第i个粒子在第 个网络中的停留时间T,算出第n个网络中灰尘的浓度： C.(..)-NATAUA (2.14) 其中，Q为源强：N为灰尘粒子总数：△工，△g,△2分别为,,z方向上的网格分辨率。 3
期末作业 2023 年 5 月 28 日 数学建模课程 Final Assignment May 28th, 2023 是排放点源 (0, 0, h) 提供的浓度； C2 = Q 2πUσyσz exp [ − ( y 2 2σ 2 y + (z + h) 2 2σ 2 z )] 是全反射像源 (0, 0, −h) 提供的浓度。 计算地面灰尘浓度即令 z = 0 ，得到高架连续点源扩散的地表浓度模型： C(x, y, 0) = Q 2πUσyσz exp [ − ( y 2 2σ 2 y + (−h) 2 2σ 2 z )] + Q 2πUσyσz exp [ − ( y 2 2σ 2 y + h 2 2σ 2 z )] = Q πUσyσz exp [ − ( y 2 2σ 2 y + h 2 2σ 2 z )] (2.8) 若要计算地面最大浓度，即求: max x C(x, 0, 0) = Q πUσyσz exp ( − h 2 2σ 2 z ) (2.9) 为简化计算不妨令 σy/σz 为常数 S ，在式 (2.8) 中令： dC(x, 0, 0) dx = 0 (2.10) 求解式 (2.9) 得到： Cmax = 2QS πUh2e (2.11) §2.3 拉格朗日扩散模型 拉格朗日大气扩散模型可以细分为两类 [4]：拉格朗日粒子模型与拉格朗日烟团模型。 2.3.1 拉格朗日粒子模型 该模型是基于湍流统计理论得到的扩散模型，可以细致追踪到每一个粒子在大气中的运 动状态，计算其时间、空间概率分布情况，以估算出灰尘浓度变化 [5]。 粒子模型将粒子速度视为风场平均速度与湍流速度之和： V = V + V ′ (2.12) 而湍流速度矢量满足： V ′ (t + ∆t) = V ′ (t)RL(∆t) + ρ (2.13) 其中，t 为时间；∆t 为时间步长；RL 为拉格朗日自相关系数；ρ 为蒙特卡洛分量。 在实际计算中，粒子模型将带预测区域划分为若干个网络，通过记录第 i 个粒子在第 n 个网络中的停留时间 Tni ，算出第 n 个网络中灰尘的浓度： Cn(x, y, z) = Q ∑N i=1 Tni N∆x∆y∆z (2.14) 其中，Q 为源强；N 为灰尘粒子总数；∆x, ∆y, ∆z 分别为 x, y, z 方向上的网格分辨率。 3
期末作业 Final assignment 2023年5月28日 数学建模课程 May 28th,2023 2.3.2拉格朗日烟团模型 该模型将同一种类的不同烟尘粒子合并为烟团，将烟尘视为有限个烟团的叠加。假设烟 团的浓度分布在水平与垂直方向上均满足高斯分布，利用22中的分析，第i个烟团对空间 内任一点(r,)的浓度作用： c=9a即(-)m(-) (2.15) 其中，Q:为改烟团的质量：？为距离点源的水平距离：z为距离点源的垂直距离。 注意到空间任一点会受多个烟团的影响，且空间内排放源是连续不断的，故该点的浓度 为所有的浓度贡献之和： C(r,2)=〉C(r,) (2.16) 对各个烟团的浓度分布进行模拟或经验求解后，通过简单叠加即可实现对整体的求解。 S2.4常见模型的分析 高斯模型是基于简化假设而建立的统计学模型，其优点在于其简单易用，适用于大气污 染物浓度分布较为均匀的情况。但是，该模型忽略了地形、建筑物等因素对大气扩散的影响 在某些场景下，如建筑密集区域、狭长山谷等地方，衍射和折射现象十分重要，而高斯模型 缺乏对这些现象的解释，因此得到的浓度分布精度很低：同时高斯模型通常不适用于时间尺 度较短、风速变化的情况，在这种情况下，必须使用微分方程模型来计算污染物扩散过程。 而拉格朗日粒子模型可模拟非稳态、复杂大气流场下污染物的扩散过程，并能有效考虑 污染物在大气中传输、湍流等因素对扩散的影响，对污染物浓度分布结果准确性更高。但是 由于其模拟了每个粒子的运动过程，计算中需耗费大量的计算资源和时间，参数比较繁多 使用不太方便，适用范围较窄，主要应用于小尺寸场景的污染扩散预测与评估。 拉格朗日烟尘模型将两者的优点结合，但仍然避免不了先验假设浓度在特定方向形成了 正态分布稳态，只要考虑地形因素，这个假设往往无法成立。 于是，我们希望能得到一个介于这两者间的模型，既不预先假设整体浓度分布，也不细 致到每个粒子，而是将浓度作为整体，考虑烟尘浓度的传输。 三、问题分析 本文考虑一个高度为h的烟肉单位时间排放质量Q的灰尘在风速为U的西风作用下 在周围地表的扩散过程。我们需要建立一个数学模型来描述灰尘的扩散过程，并且通过实验 验证模型的有效性。具体来说，除了己知条件以外，我们需要在求解浓度变化前先回答以下 问题：
期末作业 2023 年 5 月 28 日 数学建模课程 Final Assignment May 28th, 2023 2.3.2 拉格朗日烟团模型 该模型将同一种类的不同烟尘粒子合并为烟团，将烟尘视为有限个烟团的叠加。假设烟 团的浓度分布在水平与垂直方向上均满足高斯分布，利用 2.2 中的分析，第 i 个烟团对空间 内任一点 (r, z) 的浓度作用： Ci(r, z) = Qi (2π) 3/2σrσz exp ( − r 2 2σ 2 r ) exp ( − z 2 2σ 2 z ) (2.15) 其中，Qi 为改烟团的质量；r 为距离点源的水平距离；z 为距离点源的垂直距离。 注意到空间任一点会受多个烟团的影响，且空间内排放源是连续不断的，故该点的浓度 为所有的浓度贡献之和： C(r, z) = ∑ M i=1 Ci(r, z) (2.16) 对各个烟团的浓度分布进行模拟或经验求解后，通过简单叠加即可实现对整体的求解。 §2.4 常见模型的分析 高斯模型是基于简化假设而建立的统计学模型，其优点在于其简单易用，适用于大气污 染物浓度分布较为均匀的情况。但是，该模型忽略了地形、建筑物等因素对大气扩散的影响， 在某些场景下，如建筑密集区域、狭长山谷等地方，衍射和折射现象十分重要，而高斯模型 缺乏对这些现象的解释，因此得到的浓度分布精度很低；同时高斯模型通常不适用于时间尺 度较短、风速变化的情况，在这种情况下，必须使用微分方程模型来计算污染物扩散过程。 而拉格朗日粒子模型可模拟非稳态、复杂大气流场下污染物的扩散过程，并能有效考虑 污染物在大气中传输、湍流等因素对扩散的影响，对污染物浓度分布结果准确性更高。但是 由于其模拟了每个粒子的运动过程，计算中需耗费大量的计算资源和时间，参数比较繁多， 使用不太方便，适用范围较窄，主要应用于小尺寸场景的污染扩散预测与评估。 拉格朗日烟尘模型将两者的优点结合，但仍然避免不了先验假设浓度在特定方向形成了 正态分布稳态，只要考虑地形因素，这个假设往往无法成立。 于是，我们希望能得到一个介于这两者间的模型，既不预先假设整体浓度分布，也不细 致到每个粒子，而是将浓度作为整体，考虑烟尘浓度的传输。 三、 问题分析 本文考虑一个高度为 h 的烟囱单位时间排放质量 Q 的灰尘在风速为 U 的西风作用下 在周围地表的扩散过程。我们需要建立一个数学模型来描述灰尘的扩散过程，并且通过实验 验证模型的有效性。具体来说，除了已知条件以外，我们需要在求解浓度变化前先回答以下 问题： 4
期末作业 Final Assignment 2023年5月28日 数学建模课程 May 28th,2023 1.烟肉排放的灰尘有怎样的性质？ 2.烟囱排放的灰尘如何在风的作用下扩散？ 3.如何验证模型的有效性？ $3.1建模假设 3.1.1灰尘性质 现实中，工厂排放的烟尘可以看作很多种不同颗粒，每种具有各异的扩散系数、质量与 表面积。扩散系数会影响扩散的速度，质量与表面积则是会影响重力与风力的作用。 为方便计算，假设烟尘颗粒相对空气足够稀薄，它们的碰撞与相互作用可以忽略不计 于是，最终的结果也可以看作每种颗粒分别独立叠加，进而只需要对特定的某种粒子考察即 可。此后不妨假设灰尘颗粒是均匀的，且其需要满足在扩散过程中质量守恒的基本要求。 此外，由于无法确定工厂具体的排放情况，假定每个时刻排出的烟尘是完全一致的，且 排出时的速度可以忽略。 3.1.2迎风扩散 由已知条件，我们假定风速V恒定，地势对风场的影响较小。在忽略了灰尘颗粒的相 互作用后，灰尘在空气中的浓度分布由四个条件决定：沿x轴正方向的风力作用、z轴方向 的重力与空气浮力复合作用、颗粒在空气中的扩散作用、边界处的情况。 这其中，风力U是给定的常数，扩散系数D也可通过实验测量，由于粒子均匀，视为 常数。因此，后文对模型应用的讨论主要侧重于重力的影响和边界的不同情况。 3.1.3有效性验证 由于我们的模型只是对单一烟尘考虑，而实际情况有各种复杂的经验公式，我们无法直 接通过测量的数据进行结果验证。 不过，我们可以将方程模拟结果与前述的高斯模型结果进行对比，并考察不同地形对结 果影响的适用性，进而验证模型的有效性。 S3.2符号说明 由于实际问题中的物理量都是存在量纲的，为了方便数学求解，需要进行无量纲化。例 如，排放强度Q的大小差别仅在于每点浓度乘常数倍，不影响分布，可设为1或其他值：而 对z坐标进行乏=2/h的无量纲化后，则可设h=1(在对不同高度的效果进行对比时可以 在1附近调整)。其他量以各自单位基准进行无量纲化后，即可以进行数学上的建模。 为方便讨论，我们的模型中采用以下的符号： 5
期末作业 2023 年 5 月 28 日 数学建模课程 Final Assignment May 28th, 2023 1. 烟囱排放的灰尘有怎样的性质？ 2. 烟囱排放的灰尘如何在风的作用下扩散？ 3. 如何验证模型的有效性？ §3.1 建模假设 3.1.1 灰尘性质 现实中，工厂排放的烟尘可以看作很多种不同颗粒，每种具有各异的扩散系数、质量与 表面积。扩散系数会影响扩散的速度，质量与表面积则是会影响重力与风力的作用。 为方便计算，假设烟尘颗粒相对空气足够稀薄，它们的碰撞与相互作用可以忽略不计。 于是，最终的结果也可以看作每种颗粒分别独立叠加，进而只需要对特定的某种粒子考察即 可。此后不妨假设灰尘颗粒是均匀的，且其需要满足在扩散过程中质量守恒的基本要求。 此外，由于无法确定工厂具体的排放情况，假定每个时刻排出的烟尘是完全一致的，且 排出时的速度可以忽略。 3.1.2 迎风扩散 由已知条件，我们假定风速 U 恒定，地势对风场的影响较小。在忽略了灰尘颗粒的相 互作用后，灰尘在空气中的浓度分布由四个条件决定：沿 x 轴正方向的风力作用、z 轴方向 的重力与空气浮力复合作用、颗粒在空气中的扩散作用、边界处的情况。 这其中，风力 U 是给定的常数，扩散系数 D 也可通过实验测量，由于粒子均匀，视为 常数。因此，后文对模型应用的讨论主要侧重于重力的影响和边界的不同情况。 3.1.3 有效性验证 由于我们的模型只是对单一烟尘考虑，而实际情况有各种复杂的经验公式，我们无法直 接通过测量的数据进行结果验证。 不过，我们可以将方程模拟结果与前述的高斯模型结果进行对比，并考察不同地形对结 果影响的适用性，进而验证模型的有效性。 §3.2 符号说明 由于实际问题中的物理量都是存在量纲的，为了方便数学求解，需要进行无量纲化。例 如，排放强度 Q 的大小差别仅在于每点浓度乘常数倍，不影响分布，可设为 1 或其他值；而 对 z 坐标进行 z˜ = z/h 的无量纲化后，则可设 h = 1（在对不同高度的效果进行对比时可以 在 1 附近调整）。其他量以各自单位基准进行无量纲化后，即可以进行数学上的建模。 为方便讨论，我们的模型中采用以下的符号： 5
期末作业 Final assignment 2023年5月28日 数学建模课程 May 28th,2023 表1：符号说明 符号 说明 单位 h 烟囱（排放点）高度 m Q 单位时间内排放的灰尘质量 kg/s U 风速 m/s 刀 扩散系数 m2/s (红，)地表位置坐标，烟肉为原点(0,0)，西风吹向为x正方向 (m,m） 垂直地表高度 m t 扩散时间 C(A,t) 灰尘浓度，表示空间点A,时间为t时的灰尘浓度 kg/(m3.s) s(r,y) 地面高度，表示(x,)处地面的之坐标 m 四、模型建立 S4.1对流-扩散方程 由流体力学，扩散作用项的表达是D△C,其中 22 △=+亦+0 是拉普拉斯算子。 而由空气运动所产生的对流项的表达是1VC,其中u(红，，云，)为流速， -(品”) 为Nabla算子。 由于x方向风速恒定，其余只有之方向速度，可化简为得到以下偏微分方程： 架=架++p(祭+器+9)+6 (4.1) 式(4.1)是一个对流扩散方程，对其有多种不同的数值求解策略同。 根据流体力学推导，扩散系数D可以表示为红，其中k是玻尔兹曼常数，T是温度， 是空气粘度，r是灰尘粒子半径，这里可不妨设其为己知常数。 U与m的系数为正，意味着均指向对应轴的负方向。U指向x轴的负方向是为方便之 后作图，而心由于主要代表重力作用，也应指向之轴负方向。 这里Co表示新生成的浓度，由单位时间生成浓度为Q,C%应满足∫Codz dydz=Q, 无量纲化后可任意指定Q为方便计算的值。由于已经假定了每个时刻排出烟尘的一致，C 与t无关。 6
期末作业 2023 年 5 月 28 日 数学建模课程 Final Assignment May 28th, 2023 表 1: 符号说明 符号 说明 单位 h 烟囱 (排放点) 高度 m Q 单位时间内排放的灰尘质量 kg/s U 风速 m/s D 扩散系数 m2/s (x, y) 地表位置坐标，烟囱为原点 (0, 0), 西风吹向为 x 正方向 (m,m) z 垂直地表高度 m t 扩散时间 s C(A, t) 灰尘浓度，表示空间点 A，时间为 t 时的灰尘浓度 kg/(m3 · s) s(x, y) 地面高度，表示 (x, y) 处地面的 z 坐标 m 四、 模型建立 §4.1 对流-扩散方程 由流体力学，扩散作用项的表达是 D∆C，其中 ∆ = ∂ 2 ∂x2 + ∂ 2 ∂y2 + ∂ 2 ∂z2 是拉普拉斯算子。 而由空气运动所产生的对流项的表达是 u · ∇C，其中 u(x, y, z, t) 为流速， ∇ = ( ∂ ∂x, ∂ ∂y , ∂ ∂z) 为 Nabla 算子。 由于 x 方向风速恒定，其余只有 z 方向速度，可化简为得到以下偏微分方程： ∂C ∂t = U ∂C ∂x + w(z, t) ∂C ∂z + D ( ∂ 2C ∂x2 + ∂ 2C ∂y2 + ∂ 2C ∂z2 ) + C0(x, y, z) (4.1) 式 (4.1) 是一个对流-扩散方程，对其有多种不同的数值求解策略 [6]。 根据流体力学推导，扩散系数 D 可以表示为 kT 6πηr，其中 k 是玻尔兹曼常数，T 是温度， η 是空气粘度，r 是灰尘粒子半径，这里可不妨设其为已知常数。 U 与 w 的系数为正，意味着均指向对应轴的负方向。U 指向 x 轴的负方向是为方便之 后作图，而 w 由于主要代表重力作用，也应指向 z 轴负方向。 这里 C0 表示新生成的浓度，由单位时间生成浓度为 Q，C0 应满足 ∫∫∫ C0 dx dy dz = Q， 无量纲化后可任意指定 Q 为方便计算的值。由于已经假定了每个时刻排出烟尘的一致，C0 与 t 无关。 6
期末作业 Final Assignment 2023年5月28日 数学建模课程 May 28th,2023 此外，假定地面边界有形状s(红，)，可能具有一定的吸收与反射比例，将在之后边界条 件的部分进行详细讨论。 S4.2有限差分法 4.2.1基本介绍 为了求解上述偏微分方程，我们采用有限差分法进行数值模拟。具体地，我们将三维空 间离散化，得到一个网格状的空间，并在每个网格点上计算灰尘浓度，根据扩散方程进行更 新。 对偏微分方程的另一个常用的离散方法为控制容积法。由于选取一定的型线可以得到和 有限差分法相同的形式，且控制容积的写法不易于从数学上进行分析，我们采取较直接的差 分方式。 假设工，y,z的模拟范围分别是1，x,h,h,[a1,剑，选取的空间步长为九，时间步长 为，则问题变为，已知此刻的三维矩阵C经*XX,如何得到下一刻的矩阵C女x’其 中xs=2+1是x方向以h为步长离散出的格点数（含边界），其他方向同理。只要能 够正确迭代，通过取之=0即可得到每一刻的地表浓度分布。 接下来介绍两个常用的迭代求解方法，也是本文实现的做法。值得注意的是，本节中介 绍的方法均不考虑边界条件，也即迭代规律只对内部的点有效。 4.2.2FTCS格式 若Ck表示h:时刻(h,jh,kh)位置的浓度，=w(h,h:),Ck表示此点对a 求n阶导数的值，由泰勒展开可推得： r=f+f-+o 2h f=f+月-f@+O h g=+)+g-月-2@+0u 2 因此有 acC+o0h) h acA-c盛t-C达+o时 2h gC4-c出+C-2C+o的 h2 7
期末作业 2023 年 5 月 28 日 数学建模课程 Final Assignment May 28th, 2023 此外，假定地面边界有形状 s(x, y)，可能具有一定的吸收与反射比例，将在之后边界条 件的部分进行详细讨论。 §4.2 有限差分法 4.2.1 基本介绍 为了求解上述偏微分方程，我们采用有限差分法进行数值模拟。具体地，我们将三维空 间离散化，得到一个网格状的空间，并在每个网格点上计算灰尘浓度，根据扩散方程进行更 新。 对偏微分方程的另一个常用的离散方法为控制容积法。由于选取一定的型线可以得到和 有限差分法相同的形式，且控制容积的写法不易于从数学上进行分析，我们采取较直接的差 分方式。 假设 x, y, z 的模拟范围分别是 [x1, x2], [y1, y2], [z1, z2]，选取的空间步长为 h，时间步长 为 ht，则问题变为，已知此刻的三维矩阵 C i xs×ys×zs，如何得到下一刻的矩阵 C i+1 xs×ys×zs，其 中 xs = x2−x1 h + 1 是 x 方向以 h 为步长离散出的格点数（含边界），其他方向同理。只要能 够正确迭代，通过取 z = 0 即可得到每一刻的地表浓度分布。 接下来介绍两个常用的迭代求解方法，也是本文实现的做法。值得注意的是，本节中介 绍的方法均不考虑边界条件，也即迭代规律只对内部的点有效。 4.2.2 FTCS 格式 若 C l i,j,k 表示 lht 时刻 (ih, jh, kh) 位置的浓度，w l k = w(kh, lht)，∂ n a C l i,j,k 表示此点对 a 求 n 阶导数的值，由泰勒展开可推得： f ′ (x) = f(x + h) − f(x − h) 2h + O(h 2 ) f ′ (x) = f(x + h) − f(x) h + O(h) f ′′(x) = f(x + h) + f(x − h) − 2f(x) h 2 + O(h 2 ) 因此有 ∂tC l i,j,k = C l+1 i,j,k − C l i,j,k ht + O(ht) ∂xC l i,j,k = C l+1 i+1,j,k − C l i−1,j,k 2h + O(h 2 ) ∂ 2 xC l i,j,k = C l+1 i+1,j,k + C l i−1,j,k − 2C l i,j,k h 2 + O(h 2 ) 7