基于拉格朗日随机悬浮微粒模型的非高斯湍流模拟方法

摘要

本发明公开了一种基于拉格朗日随机悬浮微粒模型的非高斯湍流模拟方法，具体为：引入皮尔森IV概率密度函数，以及基于接受-拒绝方法适用于任意概率密度函数的随机数生成器，弥补原模型模拟非高斯湍流情况的不足，即控制其偏度和峰度。通过实测数据验证该方法能够在原模型基础上有效地控制其统计量中的偏度和峰度，使得该方案能够有效地模拟统计特征服从非高斯分布的湍流，提高模型的仿真度。

说明书

技术领域

本发明涉及基于拉格朗日随机悬浮微粒模型的非高斯湍流模拟方案，具体是一种关于近 地大气的非高斯湍流模拟方案，能够利用皮尔森随机数控制湍流模拟中的偏度和峰度统计量， 提高湍流仿真度和悬浮微粒的扩散预测能力，属于微气象学领域。

背景技术

拉格朗日随机模型(Lagrangian Stochastic Model简称LSM)，是一种悬浮微粒在三维空间 的运动轨道模拟模型，LSM利用数值模拟方式，替代使用动量方式模拟的欧拉微粒扩散模型， 克服欧拉模型的准确性以及在应用中的诸多限制。LSM是迄今为止世界上最广泛使用和流行 的一种微粒扩散模型，多应用于空气质量管理，自然资源管理(如农业和深林管理和防控)， 人体健康管理等。

LSM关注的悬浮微粒(Aerosol，简称微粒)具有以下性质：

1)所有模拟的微粒是单个存在的，即每个微粒的运动过程是独立的；

2)每个微粒是足够小的(一般小于10μm)，其大小小于湍流的最小长度规模(Kolmogorov 规模)，因此其扩散轨道与湍流扩散运动相同；

3)不考虑微粒之间的相互作用和微粒聚集。

微粒扩散(Aerosol Dispersion)：它是LSM的核心，其框架在一个运动维度上表示为：

$\mathrm{dx}=(u+\bar{u})\mathrm{dt}---\left(1\right)$

式中，dx是微粒运动时间dt下的运动距离；是平均速度；u是速度变化值；通常u的 均值为0，因此决定微粒速度概率分布的均值。

u＝qσ            (2)

式中，σ是微粒速度概率分布的方差。

q^t+dt＝αq^t+βr^t+dt         (3)

公式(3)是一个关于运动速度变化的马尔科夫链(Markov Chain)；其中α是加速度系数， 它是受历史状态所影响的，并且它是关于位置x，u，t的一个函数关系；而β是随机力(Random  Forcing，简称RF)r的系数；r是一个随机数，且模型假设RF服从高斯分布。

非高斯湍流(Non-Gaussian Turbulence)：LSM假设湍流是服从高斯分布，而实际测 量数据分布的均值和方差确与高斯分布相同，但实测数据分布具有较高的偏度和峰度。目前 在实际应用中，LSM仍假设湍流统计分布服从高斯分布，为了符合实测数据依据，更加准确 地模拟悬浮微粒运动状态，应修正其湍流分布假设为非高斯分布。

良混合条件(Well-Mixed Condition，简称WMC)：是业界认为最重要，严谨的一个评价模 型好坏的标准，已证明其它标准都可归纳为WMC。WMC要求框架公式(1)至公式(3)中 选择参数时应保证统计特征应与实际一致，且模拟出的统计特征应该具有稳定性和一致性。

目前针对非高斯湍流模拟的LSM主要采用以下所述几种措施：

1)非高斯RF：让原本服从高斯分布的RF服从一种非高斯分布，以控制其偏度和峰度。 虽然该方法在一定程度上能够模拟较为准确的偏度和峰度，但其稳定性差，不符合 WMC标准的要求。

2)两高斯联合分布(Two Gaussian Joint-PDFs)：结合两个高斯分布，并控制两分布的结 合系数，进而控制偏度和峰度。该方法可行，但实验结果不如原来的高斯LSM效果 好，其原因在于获得的目标概率分布函数不准确。

发明内容

本发的明目就是针对上述背景技术的问题与不足，提供一种基于拉格朗日随机悬浮微粒 模型的非高斯湍流模拟方案，引入皮尔森随机数，能够更准确，稳定地控制偏度和峰度统计 量，进而提高模型仿真度。

本发明所涉及的基于拉格朗日随机悬浮微粒模型的非高斯湍流模拟方法，其实施步骤如 下：

1)使用风速仪获取实地湍流数据：顺风速度u，侧向风速v，垂直风速w，实时温度t；

2)分段处理并分析湍流数据，计算模型输入：平均风速及其方向偏度S_w，峰度K_w， 摩擦速度u^*，奥布霍夫长度L；

3)模型模拟非高斯湍流运动状态；

进一步，所述步骤3)在原模型基础上做改进，具体为：

模型框架中的公式(2)改为公式(4)并添加公式(5)和公式(6)

u＝(ep+(1-e²)^1/2q)σ          (4)

p＝p(M_I＝0,V_I＝1,S_I,K_I)          (5)

式中，p是皮尔森随机数；e是结合系数，保证p和q结合之后均值M_w和方差V_w不变， 分别为0和1；S_w,K_w分别是实测数据的偏度和峰度统计量；p的输入均值M_I和方差V_I为默 认值0和1；

{S_I,K_I,e}＝f(S_w,K_w)          (6)

A.公式(6)，p的输入偏度S_I和峰度K_I，以及e是关于S_w,K_w函数关系，具体如下：

a)根据K_w计算e,K_I：

e＝-25.9969K_w²+161.4162K_w-249.9756,K_I＝4,K_w∈(3.0200,3.0968]        (7)

e＝-0.6774K_w²+5.1138K_w-8.7387,K_I＝4,K_w∈(3.0968,3.6041]           (8)

e＝0.3137K_w-0.3275,K_I＝5,K_w∈(3.6041,3.9184]                 (9)

e＝0.2593K_w-0.1653,K_I＝6,K_w∈(3.9184，4.1089]             (10)

e＝0.2335K_w-0.0989,K_I＝7,K_w∈(4.1089,4.2782]              (11)

e＝0.9,K_I＝19.2881K_w-77.0168,K_w∈(4.2782,4.4155]           (12)

b)根据S_w计算S_I：

S_I＝f_s(S_w)＝1.7172S_w+0.0568                       (13)

B.公式(4)中p是由皮尔森IV密度函数和随机数生成器组成，具体如下：

a)皮尔森IV密度函数表示为：

$f\left(x\right)=k{[1+{\left(\frac{x-\lambda }{\alpha }\right)}^2]}^{-m}\exp [-{v\tan }^{-1}\left(\frac{x-\lambda }{\alpha }\right)],(m>1/2)---\left(14\right)$

公式(7)中k是密度函数的正规化系数，由公式(8-10)计算获得：

$k=\frac{\Gamma \left(m\right)}{\sqrt{\pi }\Gamma (m-1/2)}{\left|\frac{\Gamma (m+\mathrm{iv}/2)}{\Gamma \left(m\right)}\right|}^2---\left(15\right)$

${\left|\frac{\Gamma (x+\mathrm{iy})}{\Gamma \left(x\right)}\right|}^2=\underset{n=0}{\overset{\infty }{\Pi }}{[1+{\left(\frac{y}{x+n}\right)}^2]}^{-1}---\left(16\right)$

$\Gamma \left(x\right)=\underset{0}{\overset{\infty }{\int }}{t}^{t-1}{e}^{-x}\mathrm{dx}---\left(17\right)$

皮尔森IV密度函数中输入参数λ，α，m和ν与M_w，V_w，S_w和K_w的关系：

β₁＝S_w²,β₂＝K_w-3              (18)

r＝6(β₂-β₁-1)/2β₂-3β₁-6              (19)

m＝0.5r+1(20)

$v=-r(r-2)\sqrt{{\beta }_1}/\sqrt{16(r-1)-{\beta }_1{(r-2)}^2}---\left(21\right)$

$\alpha =0.25\sqrt{\sigma [16(r-1)-{\beta }_1{(r-2)}^2]}---\left(22\right)$

$\lambda =M_w-0.25(r-2)\sqrt{{\beta }_1}\sqrt{V_w}---\left(25\right)$

b)随机数生成器：采用接受-拒绝方法(Acceptance-Rejection Method)根据任意概率密

度函数生成稳定的随机数，该算法具体步骤如下：

x＝r(min,max)                  (24)

s＝r(0,upper)                  (25)

y＝p(x)                        (26)

第一步，输入目标随机数的最小值min，最大值max和p最大概率值upper；第二步，计 算公式(24)和公式(25)，其中r是均值随机数生成器，p是目标密度函数；第三步，当s>y 时，重复公式(24)至公式(26)；否则，返回x作为生成的随机数；第四步，重复步骤直到 获得预定数量的随机数。

本发明加入了皮尔森随机数并向模型，根据实测数据对其峰度和偏度进行控制，更符合实 际的非高斯湍流现象，提高了湍流模拟仿真度；皮尔森随机数的引入方式与原模型的均值和 方差处理方式相同，都是从模型主体(马尔科夫链)中剥离出来，更利于分离物理意义和计 算处理；随机数生成器能够生成稳定的，符合目标密度函数的随机数，使得该方案满足良混 合条件标准。

附图说明

图1至图13是实测数据情况(Wind)，原模型(Gaussian)和改进后模型(Pearson)模拟 情况对比：其中包括在三个维度(即顺风方向u，侧风方向v，垂直方向w)，关于统计量均 值(Mean)，方差(Variance)，偏度(Skewness)和峰度(Kurtosis)，以及摩擦速度(u^*)的 对比情况：

图1至图3分别是u，v，w三维度的均值的对比情况；

图4至图6分别是u，v，w三维度的方差的对比情况；

图7至图9分别是u，v，w三维度的偏度的对比情况；

图10至图12分别是u，v，w三维度的峰度的对比情况；

图13是关于摩擦速度u^*的对比情况。

具体实施方式

下面结合附图和实施例对本发明进一步说明。应理解该实施例仅用于说明本发明而不用 于限制本发明的范围，本领域技术人员对本发明的各种等价形式的修改均落于本申请所附权 利要求所限定的范围。

实施例：

本实施例是基于拉格朗日随机悬浮微粒模型的非高斯湍流模拟方案的实例，具体如下：

1)利用10Hz 3D超声风速仪检测目标地区的湍流数据：顺风速度u，侧向风速v，垂直 风速w，实时温度t；

2)基于实测数据进行处理和分析，计算出以30分钟为区间的模型输入，其中包括：

A.计算各区间平均风速及其方向

B.计算出波动速度u',v',w'：

$u^{\prime }=u-\bar{u},v^{\prime }=v-\bar{v},w^{\prime }=w-\bar{w}---\left(27\right)$

C.计算各区间统计量均值M_w，方差V_w，偏度S_w，峰度K_w；

D.计算摩擦速度u^*：

$u^{*}={({\left(\overline{u^{\prime }{w}^{\prime }}\right)}^2+{\left(\overline{v^{\prime }{w}^{\prime }}\right)}^2)}^{0.25}---\left(28\right)$

式中，形式是指x与y的相关系数：

$\overline{x^{\prime }{y}^{\prime }}=E\left((x-E\left(x\right))(y-E\left(y\right))\right)=E\left(\mathrm{xy}\right)-E\left(x\right)E\left(y\right)---\left(29\right)$

其中E(x)是x的期望；

E.计算奥布霍夫长度L：

$L=-{\left(u^{*}\right)}^3\bar{\theta }\left(\mathrm{Kg}\left(\overline{w^{\prime }{\theta }^{\prime }}\right)\right)---\left(30\right)$

式中，K是Von Karman系数，设为0.4；g是重力加速度，设为9.81；是用平均潜在 温度代替θ，其具体计算方式如下：

$\theta =T^{\prime }{\left(P_{0}P\right)}^{(R/c_p)}---\left(31\right)$

式中，P₀是标准参考虚温度，设为1000hPa；P是测量地大气压，设为133.32hPa；T' 是绝对温度，等于t+273.15；R是气体常数，设为8.314；c_p是相关系数，设为1.005；

F.根据u^*和L计算模拟的均值

式中，K同公式(30)中系数；z₀是粗糙度，设为0.002；

G.根据u^*和L计算模拟的方差σ：

σ_u＝σ_v＝2.4u^*,σ_w＝1.25u^*,L≥0    (33)

σ_u＝σ_v＝u^*(4+0.6(-z_i/L)^2/3)^1/2,σ_w＝1.25u^*(1-3z/L)^1/3,L＜0     (34)

式中，z_i设为1000；

3)根据实测模型输入进行非高斯湍流模拟：

A.模型运行在一个圆形空间中，且微粒都是以中心点出发按照顺风方向移动，计算空间 从(x,y,z)转换为(r,θ,z)保存微粒的位置状态；其中r是离中心起点的半径，θ是与 以北方为0°且逆时针方向增加的夹角，θ₁是原始角度，z是高度：

$\left(\begin{array}{c}x=r\cos (q-q_1)\\ y=r\sin (q-q_1)\\ z=z\end{array}\right)---\left(35\right)$

$\left(\begin{array}{c}r=\sqrt{x^2+y^2}\\ q=\left(\begin{array}{c}{q}_1+{\tan }^{-1}(y/x),(x>0)\\ q_1+{\tan }^{-1}(y/x)+p,(x<0)\end{array}\right)\\ z=z\end{array}\right)---\left(36\right)$

B.计算马尔科夫链：

$q_u^{t+\mathrm{dt}}={\mathrm{aq}}_u^t+b(c_u{r}_u^{t+\mathrm{dt}}+c_w{r}_w^{t+\mathrm{dt}})---\left(37\right)$

$q_v^{t+\mathrm{dt}}={\mathrm{aq}}_v^t+{\mathrm{br}}_v^{t+\mathrm{dt}}---\left(38\right)$

$q_w^{t+\mathrm{dt}}={\mathrm{aq}}_w^t+{\mathrm{br}}_w^{t+\mathrm{dt}}+{\mathrm{gt}}_L\frac{{\mathrm{ds}}_w}{\mathrm{dz}}---\left(39\right)$

式中，α，β，γ，c_u和c_w参数表示为：

$\alpha =1-\mathrm{dt}/{\tau }_L,\beta =\sqrt{1-{\alpha }^2},\gamma =1-\alpha ---\left(40\right)$

dt＝0.025τ_L,τ_L＝l/σ_w              

$c_w=-u^{*}/{\sigma }_w,c_u=\sqrt{1-{c_w}^2}---\left(42\right)$

其中，尺度l表示为：

l＝0.5z(1+5z/L)^-1,L≥0           (43)

l＝0.5z(1-6z/L)^1/4,L＜0           (44)

C.计算波动速度：

$u=(e_u(c_u{R}_u+c_w{R}_w)+\sqrt{1-{e_u}^2}{q}_u^t){\sigma }_u---\left(45\right)$

$v=(e_v{R}_v+\sqrt{1-{e_v}^2}){\sigma }_v---\left(46\right)$

$w=(e_w{R}_w+\sqrt{1-{e_w}^2}){\sigma }_w---\left(47\right)$

式中e_u,e_v,e_w是结合系数；R_u,R_v,R_w是皮尔森随机数，在上述发明内容中详细阐述，由 皮尔森IV密度函数和随机数生成器组成：

$R_u=R_u(S_{I_u},K_{I_u}),R_v=R_v(S_{I_v},K_{I_v}),R_w=R_w(S_{I_w},K_{I_w})---\left(48\right)$

式中，R_u,R_v,R_w的各输入计算如下：

a)式中，是关于K_w的一个函数：

{e,K_I}＝f_k(K_w)             (49)

具体根据公式(7)至公式(12)计算；

b)式中，计算如下：

S_I＝f_s(S_w)                 (50)

具体根据公式(13)计算；

D.微粒移动距离计算

$\mathrm{dx}=(\bar{u}\left(z\right)+u)\mathrm{dt}---\left(51\right)$

dy＝vdt                 (52)

dz＝(w-v_s)dt              (53)

式中v_s是平均沉降速度，设为-4.7303e-04；

E.各个区间模拟的出的湍流速度为：

$U=\bar{u}\left(z\right)+u,V=v,W=w-v_s---\left(54\right)$

F.计算出各个区间湍流U,V,W的摩擦速度u^*和统计量：均值，方差，偏度，峰度

G.图1至图13比较了改进后模型模拟情况(Pearson)，与实际测量数据情况 (Wind)，和原模型模拟情况(Gaussian)的u^*以及在三个方向u,v,w的统计量；改 进后，偏度和峰度模拟情况明显改进，与实际测量情况很接近，仿真度高；均值， 方差，摩擦速度与原模型模拟情况差不多。