大气环境数值模式论文

合集下载
  1. 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
  2. 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
  3. 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。

大气环境数值模式
课程论文
题目酸沉降之区域尺度的三维欧拉模式及其应用
学生姓名
学号
院系
专业
指导教师
二OO七年六月二十日
酸沉降之区域尺度的三维欧拉模式及其应用
摘要:大量资料和研究表明,我国酸雨现象已经非常明显。

酸沉降是一个区域性的污染问题,由于他对人类环境的巨大危害性,人们采用多种手段对这一问题进行了广泛深入地研究,其中数值模拟就是一种重要的理论研究方法。

本文主要介绍了硫沉降过程,对应的欧拉数值模拟方法及其应用举例。

关键词:酸沉降;干沉降;湿沉降;欧拉模式
1. 引言
大气中的酸性物质以干或湿的形式降落在土壤、水体(河、湖)、森林及植物表面,引起了生态环境系统的破坏,使酸沉降问题成了当前一个重要的环境问题。

在大气中形成酸的前体物是硫和氮的氧化物及某些有机成分,这些物质有很大一部分是由于人类的活动造成的。

因此,为了达到减轻酸雨污染的目的,首先要研究排放的污染物与酸雨形成的关系。

通过野外观测可以了解当前环境被污染的状况,建立较为完善的数学模式,可以分析各种条件的变化,尤其是排放源的改变或布局的调整引起环境的改善或恶化,并且了解各种物理或化学现象的相互影响及其对形成酸雨的相对重要性,从而确定最优的防治措
施[1]。

多年来, 人们对酸雨进行了多方面的观测分析和模拟研究, 取得了丰硕的成果。

这些研
究主要包括云和降水酸度及化学成分的测定分析[3,4 ]和酸雨的数值模拟研究。

随着计算技术
和计算机的发展, 酸沉降模式的研究受到了人们的重视。

近年来发展了一批欧拉酸沉降模式,
其中比较著名有: RADM , STEM , ADOM [4]。

2.硫沉降模式中的物理和化学过程
2.1干沉降
干沉降过程是影响污染物浓度时空分布和长距离输送的重要过程, 越来越多的污染物质量和传输模式考虑了干沉降的作用。

一些污染物输送和沉降模式中, 假定各种污染物具有相同的干沉降速度, 而且干沉降速度不随大气条件和下垫面变化而取为常数。

对气体和粒子干沉降速度的观测研究表明, 干沉降速度的大小随大气状态、污染物种类及下垫面类型的不同出现很大变化。

为解释这种不同, 发展了一系列干沉降模型, 多数都以阻力模型为基础。

其理论基础是假定污染物的干沉降速度与其在干沉降过程中受到的阻力成反比,
1(@)d a b s
V r r r -=++ (1) 其中, r a + r b 合称为大气阻力, 表示了大气状态对沉降物的作用, r a 为空气动力学阻力, 主要由大气边界层内表面层以上湍流混合作用引起, r b 称为附加边界层阻力(副层阻力) , 由表面层内的湍流混合和分子扩散共同作用引起。

r s 为表面阻力(覆盖层阻力) , 由沉降物及下垫面的化学、物理及生物特征共同决定。

我们对输送模式中干沉降模型的设计, 力求简洁而又包括主要的物理过程, 也选用了
阻力模型。

参照RADM 中对干沉降的处理及Voldner [3]的工作, SO 2气体和SO 42-气溶胶粒子干沉
降过程中受到的阻力如下。

2.1.1 SO 2气体干沉降过程中受到的阻力
0a (z/z )r ln c ku
ψ-= (陆地和冰冻地表) (2) 23c
b 5S r u = (非雪面、非水体) (3)
c a b ln(ku Z/D)-ψr +r =ku (水体) (4) b r =2.0 (雪盖、冬季的草地和耕地) (5) u*是摩擦速度, κ为卡曼常数, Sc 为施米特数, D 为分子扩散率, ψc 为污染物稳定度相关函数。

r s 为常数, 其值随季节、昼夜和下垫面类型不同。

2.1.2 SO 42-干沉降过程中受到的阻力
ra 同SO 2的r a 。

r b 为常数,其值随季节、昼夜和下垫面类型不同, 见文献[5] 。

r s = 0 。

2.2 气相化学
气相化学在成酸过程中的作用十分重要, 有人估计, 大气中的SO 2有百分之几是通过气相氧化过程而去除的, 同时, 许多在液相化学中起着重要作用的氧化剂, 如H 2O 2、O 3等都与气相化学过程密切相关。

显然, 一个完整的酸沉降输送模式应该包括气相化学过程, 而且气相化学模式的好坏直接影响着液相化学过程。

何东阳等[6]在Atkinson [7]模式(CBM IV) 的基础上简化得到一个光化学反应模式, 该模
式包含28个物种, 46个化学方程。

模式节省了计算工作量, 而模拟结果与CBM IV 结果的对比表明, 该模式能够模拟出大气中主要的化学过程及物质的浓度变化。

模式是为研究大气中
O 3的变化而设计的, SO 2并不是模式所关心的物质。

我们在何东阳等[6]工作的基础上加入了SO 2
在大气中的主要氧化机制, 并对模式中光化学反应及对光强的依赖性做了更明确的描述, 建立了一个新的气相化学机制用以研究大气中SO 2的气相氧化问题。

该机制包括52个化学反应方程, 31个物种。

下列反应是对流层中SO
2气相氧化的主要机制。

气相化学模式中主要物质的初始浓度如表1 。

表1 气相化学模式中主要物质的初始浓度
OH + SO2 →SO 42- k = 1.50 ×1013 T -4 (cm 3mol - 1s - 1) (6)
CH 3CHOO + SO 2 →CH 3CHO + SO 42- k = 3.0 ×104 T - 1 (cm 3mol - 1s - 1) (7)
CH 3O 2 + SO 2 →CH 3O + SO 3 k = 4.3 ×10 - 15 (cm 3mol - 1s - 1) (8)
HO2 + SO2→OH + SO3k = 1.2 (cm3mol - 1s - 1) (9) 2.3液相化学及SO2和SO42-的湿清除
当有云水和雨水存在时发生的液相化学反应是SO2转化成SO42-的有效途径, 与之相伴的SO2 和SO42-等污染物从大气中的去除过程称为湿清除。

湿清除过程是SO2 气体和SO42-气溶胶清除的重要机制。

云滴、雨滴吸收污染气体并发生液相化学反应, 这一过程包括溶解、离解和化学反应三个步骤。

溶解过程和电离过程是可逆过程, 在一定的温度下达到平衡, 存在一个溶解平衡常数H 或电离平衡常数K。

不可逆化学反应是SO2在液相中氧化成SO42-的实现过程。

SO2的液相氧化包括O3、H2O2和O2的催化氧化。

云水和雨水中的H2O2和O3是比O2重要得多的氧化剂, 而H2O2尤为重要, 虽然大气中H2O2的浓度很低, 平均只有10- 9左右, 但H2O2是高水溶性物质, 当达到平衡时, 液相中的H2O2浓度几乎比液相中O3浓度高6个数量级。

大量的研究表明, 在酸性条件下, H2O2可迅速氧化水中的S( IV) 生成S(VI) , 此反应可在云内发生也可在雨滴下落过程中进行。

模式中包括了SO2在液相中被H2O2和O3氧化的不可逆化学反应。

有关化学反应及常数[8]列于表2。

表2液相化学反应及其常数
污染气体和大气气溶胶的湿清除过程都可分为云内过程和云下过程, 两种过程分别洗脱大气中不同层次上的气体和气溶胶。

对气体而言, 其湿清除过程是从云滴形成时开始的。

云滴形成以后, 它们也吸收大气中的微量气体, 降水发生时, 这些微量气体被带到地面, 这个过程叫做微量气体的雨冲刷湿清除。

降水滴在下降过程中还将继续吸收云下的气体并把它们带到地面, 这个过程为微量气体的水冲刷湿清除过程。

由于气体的性质不同, 其被清除的量及过程也不相同, 如果被清除的气体不与云滴或雨滴内的其他物质发生化学反应或只发生快速平衡可逆反应, 则问题比较简单, 这种气体的清除量完全由它们在水中的溶解度决定。

但是, 在实际大气中有许多气体能与液滴中的其他物质发生相互作用, 其清除过程也就复杂化了。

就SO 2气体的湿清除而言, 不仅要考虑SO2溶入液滴中的质量传输, 也要考虑在液滴中的溶解、电离和氧化反应。

SO 2溶于水后一部分被氧化成S(VI) , 另一部分仍以S( IV) 的形式存在于液滴中。

这就是说SO 2一部分是以S(VI) 的形式被清除的, 另一部分则以S( IV) 的形式被清除。

气溶胶粒子(SO 42-) 的湿清除与气体有很大的不同, 气溶胶粒子的湿清除效率与其粒径
有关。

大气中气溶胶粒子尺度跨度非常大, 难于统一考虑, 故在模式中分三段处理: 爱根核(0.001μm ≤r ≤0.1μm) , 大核(0.1μm < r < 1.0μm) 和巨核(1.0μm <r ≤10.0μm)。

那么气溶胶粒子的质量守恒方程为 1b h1dM =-M -M dt
(10) 2h2p2dM =-M -M dt
(11) 3h3p3dM =-M -M dt
(12) M1 、M2 和M3 分别表示爱根核、大核和巨核的浓度, Mb 为扩并下沉项, Mh 为核化作用项, Mp 为碰并作用项。

假定云滴谱服从赫尔吉安- 马逊分布, 雨滴谱服从Mar2shall - Palmer 分布。

半径小于0.1μm 的爱根核段气溶胶粒子服从Gamma 分布, 其谱表示为
N(r)=Arexp(-40r) (13)
N( r) 表示半径为r 的气溶胶粒子的数密度, A 为常数。

大核和巨核谱段的谱分布假定为Junge 分布, 谱表示为 i i 3C N (r)=r (14)
i = 2 , 3 分别表示大核和巨核。

C 为无量纲常数。

在以上的谱假定下, 可以求出方程(10) ~ (12) 右端的具体表达式, 解方程得到由于湿清除引起的气溶胶浓度的变化。

2.4非均相化学
尽管关于在气溶胶表面上发生的非均相反应对NOx 和SO 2的氧化的贡献也进行了一些研究, 但是, 由于实验方法, 颗粒物试样及模拟方法的不同, 所得的结论差别较大, 而且有
些还是相互矛盾的。

关于大气中气溶胶对SO 2转化成SO 42-的模式研究的工作还不多, 我们把气
相化学反应、液相化学反应和气溶胶谱模型结合起来, 设计了一个气溶胶表面上SO 2氧化的非均相化学反应模型, 并把该模型引入长距离输送模式, 期望对气溶胶表面上的非均相反
应增加一些认识。

模式设计的思路是这样的, 大气中可溶性气溶胶在一定的温度和湿度条件下, 发生潮解, 即在气溶胶表面形成一层薄薄的水层, 环境大气中的SO 2气体和氧化剂O 3、H 2O 2等通过气溶胶吸附和分子扩散作用进入水中形成溶液。

在溶液中发生液相化学反应, 生成的S( IV) 被吸附在气溶胶表面, 这样在潮解的气溶胶表面上不断地有S( IV) 被氧化成S (VI)。

气溶胶表面上吸附气态物质的量的多少与气溶胶的物理化学特性、其表面的潮湿程度及气相中的物质浓度有关。

在气溶胶表面上单位时间内S(VI) 生成量的多少则取决于液相反应的机制和速率。

气溶胶尺度分布也假定为前文中的三段谱即Gamma - Junge 谱分布。

求取各段气溶胶的平均半径分别为r1 , r2 , r3 。

根据Kohler 公式和气溶胶谱, 可得到一定温度、湿度和气溶胶浓度时单位体积大气中凝结在气溶胶表面上的液水总量, 进而结合气相化学和液相化学各模式而得到了描述无降水发生时气溶胶表面非均相化学反应模型。

联立以上各物理和化学模块即可建立完整的硫沉降模式。

3.欧拉型模式概况
本文所描述的是区域尺度的三维欧拉模式,计算的水平和垂直范围可由实际问题的需要及计算机的性能来定。

模式要求的输入资料包括污染源的排放清单,地面各网格点的海拔高度,气象台站的观测资料(如风向、风速、温度、气压、云量、露点及降水量等),包括地面站和高空站。

3.1大气扩散方程
控制大气中污染物的行为的大气扩散方程为 i i i i i C (c )K c R S t
v ∂+∇⋅=∇⋅⋅∇++∂ (15) 其中Ci 是第i 种气相物质的浓度,是风速向量,K 是涡旋扩散率,Si 是源项,Ri 是化学转化项,为了用于非均匀的复杂地形,取一无量纲地形订正因子为 z -h(x,y)
σ=H(x,y,t)-h(x,y)
(16) 其中z 是高度,h(x,y)是地形的高度,H(x,y,t)是模式的顶。

在这一无量纲坐标系中作坐标变换后,扩散方程变为 i i i i i i x i x y i i ()()()()()()()Hc u Hc v Hc Wc t x y c c K c K H K H x x y y H R H S H
σ
σσ
∂∂∂∂+++∂∂∂∂∂∂∂∂∂∂=+++∂∂∂∂∂∂+ (17)
其中
(,,)(,)
()()H H x y t h x y h H h H H W w u v x x y y t σσσ∆=-∂∂∆∂∂∆∂∆=-+-+-∂∂∂∂∂ (18)
3.2数值计算方法及边条件处理
在做模拟计算时,模拟范围取作700×700km 2,水平网格距△x=△y=20km.垂直方向从地
面到3km,共分15层,分别为σ=0.00032,0.0013,0.0034…0.5,0.7,1。

对方程(3)做数值积分时,积分的时间步长为△t=300s.对离散化的差分方程使用Crank-Nichelson 方法及二阶中心差分做数值积分。

化学反应项的计算是根据化学反应方程式进行的。

将模式边界取在离源足够远的地方时,计算的边条件为 i x max K 0,(0,)C x H x x x x
∂∆===∂方向 (19) i y max y K 0,(y=0,y=y )C H y
∂∆=∂方向 (20) i z i g K z ,(0.00032)C v H σσ
∂==∆∂方向 (21) z i 0,(1)C H σσ
∂==∆∂K
其中V x i 是第i 种物质的干沉降速度,按照定义
Vg=F/c(g)
其中F 是物质的垂直通量,c(g)是地面的物质浓度
F=Kz ∂c /∂z
Kz 是垂直涡旋扩散率,可由气象资料及地面粗糙度计算。

大气稳定度分类是根据Pasquill 的分类方法,按照太阳辐射强度、云量和风速分成稳定、中性和不稳定三类。

对不同的大气稳定度,大气扩散参数表示为
不稳定 1/3z z (14/)()1 4.7/ku z H K L z L H
⋅-=-+ (22) 中性 1/3()z H z K ku z H
-=⋅ (23) 稳定 1/31/2z (11/)()(116/)z ku z
H K L H
z L -⋅-=-- (24)
其中k 为卡曼常数,z 是计算高度,u *是摩擦速度,L 是Monin-Obukhov 长度,H 是混合层高度。

对不同的稳定度类,u *的计算式为
不稳定
11ln 2ln 211zw p r u ku arctgp arctgp p p ⎧⎫⎛⎫⎛⎫--=+⋅-+⋅⎨⎬ ⎪ ⎪++⎝⎭⎝⎭⎩⎭
(25) 其中 p=(1-15zw/L)
1/4 中性 0/ln(/)zw u ku zw z = (26)
稳定 004.7/ln(/)()zw u
ku zw z zw z L ⎡⎤=+-⎢⎥⎣⎦ (27) 计算时zw 取10m ,u zw 为10m 的风速。

Z 0是地面粗糙度,对不同类型地表特征有不同的值。

L 用Pasquill 和Smith 给出的值。

水平扩散参数取作k z =k y =Ak z ,对不同稳定度取不同A 值。

4.个例模拟
为了检验模式的综合模拟性能, 以1991 年6 月12~14日淮河流域的一次降雨过程作为
实例, 对我国东半部及其相邻地区(2400×2340km 2, 中心纬度36o N , 经度117o E) 的酸沉降
状况进行了数值模拟。

气象场由中尺度模式MM4预报得到, 与实况资料比较一致[9,10]。

这次
降雨过程由中尺度对流云团引起, 雨区呈东西向带状分布, 24 h 累积降水最大达197 mm (图
1)。

在长江以南地区, 整个对流层中下层为一致的偏西南气流, 850 hPa 和700 hPa 等压面上分别有风速大于12m/s 和16m/s 的西南低空急流(图2) , 地面气压图上, 我国南方呈现东
高西低的形势[9]。

计算时水平格距取为60km,格点数41×40, 铅直方向分10层, 对应的R 值是1.0,0. 99, 0.97, 0.93, 0.85, 0.75, 0.60, 0.45, 0.3, 0.15, 0.0,顶层气压100hPa 。

利用算子分裂
法求解输送物种的浓度变化方程, 平流项采用Smolarkiewicz 修正的上游差分方案[11] , 铅
直扩散项采用显式的Crank- Nicolson 方案, 采用Hybrid 方案及误差分配质量守恒技术求
解非线性气相化学动力学方程组, 液相化学动力学方程组采用QSSA 方法求解[12]。

积分时源排
放、平流、扩散和干沉积的时步取为10 m in,气相化学过程的时步为60s, 与云雨有关的过程每小时强迫一次, 等过程时步取30s 。

模拟结果如下。

4.1 SO 2和SO 42-的浓度
模拟的SO 2日均地面浓度分布特征表现为浓度呈现明显的不均匀性, 在江苏及山东半岛、华北平原及辽东半岛和四川盆地三个地区存在高值中心, 对应着这些地区有较强的工业污染源, 这说明我国SO 2污染严重地区主要是局地源形成的, 排放源对地面浓度分布起了决定性的作用。

四川的总排放量虽然大于江苏, 但就排放强度而言, 江苏要大于四川, 然而在江苏地区并未形成特别大的高值浓度, 原因在于该地区多为平原, 东部邻海, 大气的稀释能力比较强, 污染物易于向外输送和扩散。

另外, 城市所在地区的网格平均浓度在10~
40ug/m 3 范围, 最大值83 ug/m 3出现在重庆及其周围地区; 远离城市的地方浓度较小,一般低于5ug/m 3;淮河流域、华中和东北地区浓度都较小, 这是降水清除带来的直接效果。

图1 6月12~13日24h累积降水(mm)分布图2 6月12日14: 00 2 km 高度流场分布
SO2浓度的高空分布表现为高层有与源排放和地面浓度相对应的大值区, 且浓度随高度递减, 在边界层顶以上的高度, 浓度已经很小, 水平方向的分布也比较均匀, 污染物受高空西风气流引导输送的轨迹非常明显。

SO42-的日均地面浓度分布基本和SO2相似, 也存在三个浓度的相对大值区, 但其峰值要比SO2低1~2个量级, 最大值出现在西南, 约3 ug/m3。

总的来说, SO42-的分布比SO2均匀, 且输送扩散的距离比SO2更远, 这是因为模式考虑的SO42-全部由气相SO2转化而来, 而这种转化需要一个时间过程。

4.2 SO2和SO42-的干沉积
图3 SO2的日干沉积量图4 SO42-的日干沉积量
单位:mg/m2,等值线间距为2 单位:mg/m2,等值线间距为0.03
图3、图4分别是SO2和SO42-的日干沉积量分布, 与地面浓度分布状况相类似, 存在三个高值区, 这意味着大气污染物在邻近源区有主要的沉降。

在城市区, SO2的干沉积量多在4~10 mg/m2, 最大可达44 mg/m2, 在广大非城市区, 其值只有1 mg/m2。

SO42-的干沉积量在数值上
远远小于SO2,这是由于雨天SO2的气相转化率较小, SO42-浓度较低, 且SO42-的干沉积速度比SO2小几倍。

4.3 云水和地面雨水中SO2和SO42-的浓度
地面雨水中SO2的浓度极低, 而SO42-的浓度较高, 在10~ 400umol/L范围, 说明SO2的氧化是SO42-形成的主要途径。

SO42在云水中的浓度通常要小于在雨水中的浓度, 前者为后者的30%~ 80% , 表明在酸雨的形成过程中, 云中清除和氧化占据主要地位, 云下过程也起到一定的酸化作用。

SO42-浓度的高值区并不完全和降水中心相对应, 而是偏于雨区边缘, 因为在强降水的地方, 虽然SO2的清除量大, 但由于降水量多, 使得浓度变稀。

4.4 云雨水的pH值
从雨水的pH值分布图上可以看出, 降水强度大的地方往往雨水的pH值高, 而雨区边缘pH值较低, 可能是那里降雨量少, H+浓度大, 且有足够SO2和氧化剂充足的缘故。

云水的pH 值分布与雨水比较类似, 且云水pH值略高于雨水pH值, 说明导致酸雨形成的过程主要在云内, 云下过程对雨水只有弱的酸化作用。

此外, 比较pH值和SO42-的浓度可以发现, 两者有着较好的负相关关系。

由于缺少基础资料, 模式没有考虑云水和雨水对其他碱性气溶胶粒子的清除作用, 所以模拟得到的云雨水pH值总体上偏低。

4.5SO2和SO42-的湿沉积
图5 SO2的日湿沉积量图6 SO42-的日湿沉积量
单位:mg/m2单位:单位:mg/m2
图5、图6 给出了SO2和SO42-的湿沉积量分布, 可以发现, SO42-湿沉积量与地面降水分布形式是一致的, 沉降量的极大值出现在降水中心附近, 因为在强降水中心, H2O2、O3等氧化剂和SO42-气溶胶粒子被有效清除而进入雨水中, 非常有利于SO2的转化, 造成较大的湿沉降量, 数值一般在100mg/m2左右。

SO2的湿沉降量要比SO42-小得多, 几乎可以忽略, 表明由于清除而溶解于水中的SO2绝大部分被氧化掉了。

4.6地区侧边通量
侧边通量指的是通过给定地区侧边界的污染物量, 以浓度与风速分量的乘积来表示, 规定输出为正, 输入为负, 它可以反映一个区域与其周围污染物输送和交换的特征。

我们选择了江苏和山东这两个SO2排放(4085t/d、5586t/d) 较大的沿海省份, 计算了两省的侧边通量。

从计算结果来看, 在贴地层, 两省在边界交换的数量并不大, 一般不超过区域排放量的10% , 2000m 以上交换量已经很小, 把四个边界的收支相抵消后, 净通量更小。

就本个例的天气条件下, 在东部边界, 江苏为整层输出, 山东是低层输入、高层输出。

上述个例模拟所得到的定量结果和定性结论与同类工作相呼应, 这说明本文建立的简化模式具有良好的模拟性能。

[13]
参考文献:
[1
[3] Hegg , D1 A1 and P1 V1 Hobbs and L1 F1 Radke , 1984 , Measurements of the scavenging of sulfate and nitrate in clouds , Atmos1 Envi ron1 , 18 , 1939-19461 [4] 黄美元等.重庆地区云水和雨水酸度及其化学组分的观测分析.大气科学.1988.12:
389-3951.
[5] Voldner , E1 C1 et al1 , 1986 , A literature review of dry deposition of oxides of sulfur and nitrogen with emphasis on long - range transport modeling in north Americ1 , Atmos1 Envi ron1 , 20 , 2101-21231.
[6] 何东阳等.一个适用于对流层光化学模拟的大气化学模式.环境科学学报.1992.12:
182-1921
[7] Atkinson , R1 , Lloyd , A1 C1 and Winges , L1 , 1982 , An updated chemical mechanism for hydrocarbons/ NOx/ SO2 photo - oxidation suitable for inclusion in atmospheric simulation model , Atmos1 Env rion1 ,16, 1341-13551
[8] 雷恒池黄美元.对流云云中酸化的数值模拟-模式建立和初步结果, 环境科学学
报.1992.12:415-4261.
[9] 潘益农等.一次暴雨过程的分析研究和数值模拟.南京大学学报.1993.29:118-122.
[10] 郑维忠等.中尺度数值模式MM 4 模拟系统介绍.南京大学出版社1992:2-16.
[11] Smo lark iew icz, P io tr K. , 1983, A simp le po sitive definte advection scheme with small imp lict diffusion, M on. W ea.R ev. , 111, 479- 486.
[12] Odman, Mehmet T. , 1992, A comparison of fast chem ical k inetic so lvers fo rair qualitymodeling, A tm os. E nv iron. ,26A (9) , 1783- 1789.
[13] 王体健李宗恺等.《区域酸沉降的数值研究Ⅱ个例模拟和试验》.Scientia Atmosphere。

相关文档
最新文档