颗粒物面源模式计算方法

摘要

对颗粒物面源模式做了深入研究，确定了颗粒物粒子面源开始沉降的位置。颗粒物粒子在面源内部不考虑沉降问题，面源源块内颗粒物浓度等于单位面积污染物源强与地面反射系数之积除以地面风速。面源源块外面开始考虑粒子沉降，采用虚点源后置法，参考颗粒物点源的部分反射倾斜烟云模式，全面给出了面源的颗粒物模式。面源边缘浓度采用线性内插方法处理，以免计算结果不连续。面源中心点位置应该在水平x，y方向平移一个微小距离，以消除计算点和面源中心点重合的可能。

正文

近年来颗粒物粒子对环境空气的影响逐渐显现出来, 对于颗粒物的控制和治理也提到各地环境保护部门的议事日程。与气质模式相比, 颗粒物模式计算方法虽然没有那么成熟, 但是也有一些可以应用的方法, 可是关于颗粒物的面源模式的计算方法却极少报道。而颗粒物却有可能是大面积的面源污染。究竟如何计算是一个新的研究课题。

颗粒物计算模式有许多种, 如源损耗模式、解析模式、倾斜烟云模式、部分反射的倾斜烟云模式等, 其中部分反射的倾斜烟云模式的实用性更好, 在大气环境模式计算里得到广泛应用^[1]。在颗粒物粒子直径小于10 μm 时采用气质模式;颗粒物粒子直径大于10 μm 时就应该采用部分反射的倾斜烟云模式。

在计算面源时一般采用虚点源后置法, 但是颗粒物面源计算有其特殊性, 因为颗粒物粒子要考虑沉降, 而从哪个位置开始沉降, 却无法确定。如果颗粒物到达面源网格中心点就开始沉降, 那么就有可能在污染物没有到达网格下风向边缘就已经沉降到地面, 其后的地面浓度应该为零, 这显然是不合理的。为此计算之前必须解决颗粒物面源从何处开始沉降的问题。谷清曾经提出过一个经验方法, 即颗粒物面源内部不考虑粒子沉降问题, 在面源内部地面浓度相同, 在面源外的一个区域内考虑粒子沉降, 从面源的边缘开始计算沉降距离和其影响^[1]。按照这样的假设, 颗粒物面源模式应该为:当x大于正方形面源边长的1/2时, 即在污染物传输、扩散到面源源块之外时, 采用虚点源后置法;当x小于等于正方形面源边长的1/2时, 面源源块内颗粒物浓度等于单位面积污染物源强乘上地面反射系数除以地面风速, 据此, 参考颗粒物粒子点源浓度的模式, 可以推导出颗粒物面源计算模式。

《1 面源内颗粒物浓度》

1 面源内颗粒物浓度

面源源块内部的颗粒物地面浓度计算公式为

$C = \frac{Q (1 + α)}{u_{d} R L^{2}}$ $C = \frac{Q (1 + α)}{u_{d} R L^{2}}$

式中:C—污染物地面浓度, mg/m³;

Q—污染物源强, mg/s;

α—地面反射系数;

u_d—平均风速, m/s;

RL—面源边长, m。

《2 面源外颗粒物浓度》

2 面源外颗粒物浓度

《2.1有风时 (u>1.0 m/s) 采用倾斜烟云模式》

2.1有风时 (u>1.0 m/s) 采用倾斜烟云模式

$\begin{array}{l} C = \frac{Q (1 + α)}{2 π U_{e} σ_{y m} σ_{z m}} \exp (\frac{- y^{2}}{2 σ_{y m}^{2}}) \times \\ \exp (\frac{- (Η - V_{g} \frac{X}{u})^{2}}{2 σ_{y m}^{2}}) \end{array}$ $\begin{array}{l} C = \frac{Q (1 + α)}{2 π U_{e} σ_{y m} σ_{z m}} \exp (\frac{- y^{2}}{2 σ_{y m}^{2}}) \times \\ \exp (\frac{- (Η - V_{g} \frac{X}{u})^{2}}{2 σ_{y m}^{2}}) \end{array}$

式中:C—污染物地面浓度, mg/m³;

u_d—平均风速, m/s;

x—面源中心起的下风向距离, m;

y—面源中心起的横向距离, m;

H—有效源高, m;

σ_ym—面源水平方向扩散参数,

$σ_{y m} = σ_{y} (x + x_{y}), m;$ $σ_{y m} = σ_{y} (x + x_{y}), m;$

σ_zm—面源垂直方向扩散参数,

$σ_{z m} = σ_{z} (x + x_{z}), m;$ $σ_{z m} = σ_{z} (x + x_{z}), m;$

x_y—根据反扩散参数求得的面源虚后移距离, m;

x_z—根据反扩散参数求得的面源虚后移距离, m;

V_g—颗粒物粒子沉降速度, m/s。

《2.2小风静风时 (u≤1.0 m/s) 采用沉降烟团模式》

2.2小风静风时 (u≤1.0 m/s) 采用沉降烟团模式

$\begin{array}{l} C = \int_{0}^{Τ} \frac{Q (1 + α)}{(2 π)^{3 / 2} σ_{x m} σ_{y m} σ_{z m}} \exp (\frac{- R^{2}}{2 σ_{y m}^{2}}) \times \\ \exp (\frac{- (Η - V_{g} t) R^{2}}{2 σ_{z m}^{2}}) d t \end{array}$ $\begin{array}{l} C = \int_{0}^{Τ} \frac{Q (1 + α)}{(2 π)^{3 / 2} σ_{x m} σ_{y m} σ_{z m}} \exp (\frac{- R^{2}}{2 σ_{y m}^{2}}) \times \\ \exp (\frac{- (Η - V_{g} t) R^{2}}{2 σ_{z m}^{2}}) d t \end{array}$

式中:T—积分时间 (静风持续时间) , s;

R—水平距离, $= R \sqrt{X^{2} + Y^{2}}, m$ $= R \sqrt{X^{2} + Y^{2}}, m$ ;

σ_ym—面源水平方向扩散参数,

$σ_{y m} = σ (t + t_{0}) ‚ m ；$ $σ_{y m} = σ (t + t_{0}) ‚ m ；$

σ_zm—面源垂直方向扩散参数

$σ_{z m} = σ (t + t_{0}) ‚ m ；$ $σ_{z m} = σ (t + t_{0}) ‚ m ；$

σ_xm—顺风向扩散参数, σ_xm=σ_ym, m;

t₀—根据反扩散参数求得的面源虚拟扩散时间, s。

《2.3长期平均浓度模式》

2.3长期平均浓度模式

$\begin{array}{l} \begin{array}{l} C = \sum_{Ι = 1}^{6} \sum_{J = 1}^{16} \sum_{Ι = 1}^{n} F_{Ι, J, Κ} g_{3} + \sum_{Ι = 1}^{n} F_{Ι} g_{4} \\ g_{3} = \frac{16 (1 + α) Q}{(2 π)_{3 / 2} u σ_{z m} X} \times \end{array} \\ \exp (\frac{- (Η - V_{g} \frac{x}{u})^{2}}{2 σ_{z m}^{2}}) (\frac{0.414 x - | y |}{0.414 x}) \end{array}$ $\begin{array}{l} \begin{array}{l} C = \sum_{Ι = 1}^{6} \sum_{J = 1}^{16} \sum_{Ι = 1}^{n} F_{Ι, J, Κ} g_{3} + \sum_{Ι = 1}^{n} F_{Ι} g_{4} \\ g_{3} = \frac{16 (1 + α) Q}{(2 π)_{3 / 2} u σ_{z m} X} \times \end{array} \\ \exp (\frac{- (Η - V_{g} \frac{x}{u})^{2}}{2 σ_{z m}^{2}}) (\frac{0.414 x - | y |}{0.414 x}) \end{array}$

式中:F_{I, J, K}—I类稳定度、J风向、k级风速的发生频率;

F_I—I类稳定度静风发生频率;

F_In—风速等级;

F_Ig₄—即沉降烟团模式。

应用上式时, 应首先定出σ_y0和σ_z0, 然后根据扩散曲线和扩散参数公式计算从面源中心向上风向推移的虚源距离。例如设扩散参数采取如下形式:

$σ_{y} = γ_{1} x^{α 1} ‚ σ_{z} = γ_{2} x^{α 2}$ $σ_{y} = γ_{1} x^{α 1} ‚ σ_{z} = γ_{2} x^{α 2}$

由σ_y0和σ_z0定出虚源至面源中心的距离是

$x_{y} = (\frac{σ_{y 0}}{γ_{1}})^{\frac{1}{α_{1}}} ‚ x_{z} = (\frac{σ_{z 0}}{γ_{2}})^{\frac{1}{α_{2}}}$ $x_{y} = (\frac{σ_{y 0}}{γ_{1}})^{\frac{1}{α_{1}}} ‚ x_{z} = (\frac{σ_{z 0}}{γ_{2}})^{\frac{1}{α_{2}}}$

用反扩散参数可以表示为

$x_{y} = a r c f (σ_{y}) ‚ x_{z} = a r c φ (σ_{z})$ $x_{y} = a r c f (σ_{y}) ‚ x_{z} = a r c φ (σ_{z})$

只要给出扩散参数就可以求出下风距离, 更详细的反扩散参数的概念和方法见参考文献[1,2]。长期平均浓度模式中的最后一项是对面源横风向修正, 避免相临两块面源计算的浓度不连续, 具体方法也见参考文献[1,2]。

上面所介绍的只对某一级粒子而言, 实际计算时, 应该对不同直径的粒子采用不同的沉降速度和地面反射系数, 然后再把各级粒子的浓度相加, 得到总的颗粒物粒子浓度。

《3 面源中心点的安放》

3 面源中心点的安放

一般来说, 面源基本源块大小与划分的网格大小相同, 这样便于统计面源源块内的污染物源强, 也符合常规的思维逻辑。当然, 这样一来面源的中心点与基本网格的中心点也重合。但是在某些情况下, 面源尺度大小需要变化, 面源与网格不重合, 其中心点在网格内随意安置, 这时可以给面源一个大小尺度的比例系数, 这样一来面源的大小可以随意给定, 其中心点也可以在网格内随意安放。计算时可以根据基本网格长度和比例系数确定实际面源大小, 再根据中心点位置分别处理面源内外不同点的计算。

《4 面源边缘的处理方法》

4 面源边缘的处理方法

过去在计算一个城市或区域的污染物浓度的时候, 由于计算机硬件设备的限制, 大都采用大网格计算, 一般是计算基本网格中心点和网格节点。现在计算机计算速度大大提高, 有可能加密计算点位, 以便模拟计算出计算区域的浓度的细微变化。按照上面的假设计算时发现, 在相临面源边界处, 计算浓度会偏高, 因为面源边界处被重复计算, 个别点可能重复多次, 从而造成面源边缘和节点处的浓度是面源内部的几倍。如果不计算面源边缘, 那么就有可能在相临两个面源交接处形成一条地面浓度或者几条地面浓度为零的深沟, 这两种现象显然是不对的。修改方法是在距离面源中心一定距离处到面源外面一个距离对浓度做线性内插。以x方向为例, 选择0.4面源边长和0.6边长之间内插, 这样在面源边缘处的浓度正好是内部浓度的二分之一, 0.6边长处的浓度为零。y方向也应该这样处理。经过这样处理, 基本可以满足浓度场连续的要求。但是在面源的4个角, 也就是4个节点附近还会出现类似这样的小毛病。按照这样的思路也可以处理平滑, 本文不再赘述。

《5 面源中心点的平移》

5 面源中心点的平移

作为面源来说, 特别是大城市和区域内计算时, 面源中心点和计算中心点非常容易重合, 一旦发生重合, 尽管设计计算程序时设计者增加了许多修正和限制, 但计算时, 还会发生计算浓度为零, 或者计算机死机和溢出现象, 这种现象发生概率还是很高的。前者相当于浓度场中出现一个黑洞, 而后者干脆不能够计算。为了避免出现这样的情况, 确定面源后, 在编制计算机程序时, 把面源中心点在平面内平移一定距离即可。如将面源中心点平移面源边长的1 %。最好x, y方向都移动。这样处理以后, 再发生上述现象的概率极小。

《6 计算实例》

6 计算实例

我国北方某地计划在海上建设一个大型煤炭港口, 按照环保要求应该对场区及周围环境大气中的总悬浮颗粒物浓度进行估算。该场区内有装船机6台, 码头6个, 卸煤机若干, 取煤机若干, 火车线, 输煤皮带运输机和堆放场。最主要的影响是堆放场。图1是模拟计算的颗粒物浓度, 图中心部位是填海形成的港口陆地区域, 向东北方向延伸的是栈桥, 三个圆点是大气监测点, 周围颜色浅的是海区和河流, 颜色深的为潮间带。堆放场形成两个大面源, 两块面源东西相临, 南北方向差半个面源长。作为TSP日平均浓度分布的图1, 上图浓度为20, 50, 100, 150, 200和250 μg/m³。在两块面源交接处出现高浓度, 最大为250 μg/m³。下图是按照本文第4章介绍的方法对面源边缘进行内插处理后的结果, 图中有10, 50, 100, 125和150 μg/m³, 最大浓度为150 μg/m³, 并且浓度分布连续合理。图中浓度是二块面源和多块小面源以及多个点源共同作用的结果 (此算例的数据发表时做了修正) 。

《图1》

图1 TSP浓度对比

Fig.1 Rate of concentration of TSP

《7 结论》

7 结论

1) 颗粒物面源模式

颗粒物粒子在面源内部不考虑沉降问题, 面源源块内颗粒物浓度等于单位面积污染物源强乘上地面反射系数除以地面风速;面源源块外面开始考虑粒子沉降, 采用虚点源后置法, 参考颗粒物点源的部分反射倾斜烟云模式, 可以全面推导出颗粒物面源模式。这一方法在地面面源和近地面面源中都可以应用。

2) 几个计算修正方法

面源应该给定一个大小尺度的比例系数来确定面源的大小, 其中心点也可以在网格内随意安放。可以根据基本网格长度和比例系数确定实际面源大小, 再根据中心点位置分别处理面源内外不同点进行计算。面源边缘浓度采用线性内插的方法处理, 以免相邻面源计算结果不连续。面源中心点位置应该在水平x, y方向平移一个微小距离, 以消除计算点和面源中心点重合的可能。

这几个问题都是属于计算方法和计算技巧, 如果处理不好会影响计算效果, 甚至无法计算。

展示更多