一种下击暴流非平稳脉动风速的数值模拟方法

摘要

本发明提供了一种下击暴流非平稳脉动风速的数值模拟方法，通过下击暴流的时变平均风速模型和一个与时间频率有关的非均匀调制函数对功率谱进行调制来获得下击暴流非平稳脉动风速时变功率谱；将平稳随机过程的AR模型进一步考虑时变特征来建立非平稳随机过程的TAR时变模型；通过时变功率谱、TAR时变模型对下击暴流非平稳脉动风速进行有效模拟。本发明模拟的下击暴流非平稳脉动风速的相关性随着不同位置之间的距离增大而减弱，振幅随着时变平均风速的增大而增大，与实际风场特性相吻合；模拟的下击暴流非平稳脉动风速的功率谱具有明显的时变特征，且与目标时变谱的时变特征相吻合，而且样本统计平均功率谱、相关函数与目标也均吻合。

说明书

技术领域

本发明属于风暴模拟技术领域，尤其涉及一种下击暴流非平稳脉 动风速的数值模拟方法。

背景技术

在雷暴天气中，往往会形成局部强下沉气流猛烈冲击地面，产 生一种极具突发性和破坏性的强风即下击暴流。人们首次认识到下击 暴流的灾害性是在1975年美国一架波音727航班降落时发生的一次 坠机事件中。对该事故调查后，Fujita首次提出了下击暴流，即能在 地面产生17.9m/s以上辐散风的一种强烈的下沉气流^[1]。Chay和 Letchford^[2]从本质上详细的阐述了下击暴流风速的形成方式。在雷暴 环境中，尺度较小的微下击暴流发生的频率很高，在近地面产生的最 大风速可达75m/s^[3]。显然其破坏性特别大，属于近地面灾害性强风。 下击暴流在美国、日本、中国、澳大利亚等世界范围内均有发生，建 筑物主体及其围护结构在其作用下会产生严重的损伤、破坏，甚至垮 塌^[4,5]。

人们已逐渐认识到下击暴流研究的重要性和必要性，国内外研 究人员通过现场观测、实验模拟、理论分析和数值仿真等方法对下击 暴流展开了广泛的研究。通过对下击暴流的数次实测，国外学者发现 下击暴流风速场中最大风速出现在近地面，而距离地面较高处的风速 相对较小，这明显不同与传统大气边界层风速场分布^[3]。Holmes和 Oliver^[6]通过碰撞射流理论研究了下击暴流行进过程中的某一固定位 置所产生的平均风速，并提出了随时间变化平均风速的模型。Wood 等^[7]通过下沉气流的风洞实验数据，提出了与计算流体动力学分析模 型相吻合的风速场平均风速模型。Savory等^[8]运用Holmes建立的模 型研究了击暴流作用下的输电塔破坏情况。由于该模型没有考虑风速 场的随机波动成分，这可能会低估了下击暴流风荷载作用下结构的动 力响应。下击暴流风速具有很强的非平稳性。我国风工程领域专家瞿 伟廉教授^[9]将下击暴流分解为一个确定的时变平均风速和一个调制 非平稳的脉动风速，进行了数值模拟研究。李春祥等^[10]采用Deodatis 提出的均匀调制非平稳随机场模拟方法模拟了下击暴流非平稳脉动 风速。张文福等^[11]通过时变函数均匀调制基于AR模型的平稳脉动风 速来生成了具有空间相关性的下击暴流非平稳脉动风速。根据 Priestley提出的进化谱理论^[12]，非平稳脉动风速的模拟，应该通过一 个与时间频率有关的非均匀调制函数对功率谱进行调制获得非平稳 风速的时变功率谱即进化谱，然后通过非平稳随机过程的模拟方法来 进行模拟。然而，为了简化，通常将非均匀调制函数假定为仅与时间 有关的均匀调制函数，将非平稳随机过程的模拟转化为一均匀调制函 数对平稳随机过程的调制，回避了根据时变功率谱模拟非平稳随机过 程的难点。而且，目前对于非平稳脉动风速的模拟，通常将时变平均 风速函数假定为平稳脉动风速的均匀调制函数并没有相关理论依据。 为此，李锦华等^[13]将非平稳脉动风速离散成若干段在足够短时间Δt内 可以近似为平稳脉动风速的短时时间序列，从而严格推导出与时间和 频率有关的非均匀调制函数和相应的非平稳脉动风速的进化谱。

发明内容

本发明的目的在于提供一种下击暴流非平稳脉动风速的数值模 拟方法，旨在解决现有技术中回避了根据时变功率谱模拟非平稳随机 过程难点的问题。

本发明是这样实现的，一种下击暴流非平稳脉动风速的数值模拟 方法，包括以下步骤：

S1、根据进化谱理论，通过下击暴流的时变平均风速模型和一个 与时间频率有关的非均匀调制函数对功率谱进行调制来获得下击暴 流非平稳脉动风速时变功率谱；

S2、将平稳随机过程的AR模型进一步考虑时变特征来建立非平 稳随机过程的TAR时变模型；

S3、通过所述时变功率谱、TAR时变模型对下击暴流非平稳脉 动风速进行有效模拟。

优选地，在步骤S1中，所述非均匀调制函数为Kaimal调制函数：

$A(\omega ,t)=\sqrt{\frac{{\tilde{U}}_j\left(t\right)}{U_j}{\left[\frac{1+50\frac{\omega z_j}{2\pi U_j}}{1+50\frac{\omega z_j}{2\pi {\tilde{U}}_j\left(t\right)}}\right]}^{5/3}}$

式中:ω为圆频率；z_j为空间某点的垂直地面高度；为空间 某点处的时变平均风速；为空间某点处非平稳风速的 统计平均风速。

优选地，所述下击暴流非平稳脉动风速时变功率谱包括 Davenport谱、Harris谱、Simiu谱调制后的非均匀调制函数。

优选地，在步骤S1中，所述Davenport谱调制后的非均匀调制 函数为：

$A(\omega ,t)={\left(\frac{1+{\left(\frac{1200\omega }{2\pi U_{10}}\right)}^2}{1+{\left(\frac{1200\omega }{2\pi {\tilde{U}}_{10}\left(t\right)}\right)}^2}\right)}^{2/3};$

所述Harris谱调制后的非均匀调制函数为：

$A(\omega ,t)=\sqrt{\frac{{\tilde{U}}_{10}\left(t\right)}{U_{10}}{\left(\frac{2+{\left(\frac{1800\omega }{2\pi U_{10}}\right)}^2}{2+{\left(\frac{1800\omega }{2\pi {\tilde{U}}_{10}\left(t\right)}\right)}^2}\right)}^{5/6}}$

所述Simiu谱调制后的非均匀调制函数为：

$\left(\begin{array}{cc}A(\omega ,t)=\sqrt{\frac{{\tilde{U}}_j\left(t\right)}{U_j}{\left[\frac{1+50\frac{\omega z_j}{2\pi U_j}}{1+50\frac{\omega z_j}{2\pi {\tilde{U}}_j\left(t\right)}}\right]}^{5/3}},& \frac{\omega z_j}{2\pi {\tilde{U}}_j\left(t\right)}\le 0.2\left(1\right)\\ A\left(t\right)={\left(\frac{{\tilde{U}}_j\left(t\right)}{U_j}\right)}^{4/3},& \frac{\omega z_j}{2\pi {\tilde{U}}_j\left(t\right)}>0.2\left(2\right)\end{array}\right);$

以上各式中，ω为圆频率；z_j为空间某点的垂直地面高度；为空间某点处的时变平均风速；为空间某点处非平稳 风速的统计平均风速；为10米高处的时变平均风速； 为10米高处非平稳风速的统计平均风速。

优选地，在步骤S2中，所述AR模型为自回归模型：

$f\left(t\right)=\underset{i=1}{\overset{p}{\Sigma }}{A}_{i}f(t-\mathrm{i\Delta t})+\mathrm{Lw}\left(t\right);$

式中，A_i表示自回归模型系数，L表示待定的模型系数。

优选地，在步骤S2中，所述TAR时变模型用函数定义为：

$f\left(t\right)=\underset{i=1}{\overset{p}{\Sigma }}{A}_i\left(t\right)f(t-\mathrm{i\Delta t})+L\left(t\right)w\left(t\right);$

式中，时变自回归模型系数A_i(t)可用函数定义为：

$R_{\mathrm{ff}}(t,\mathrm{j\Delta t})=\underset{i=1}{\overset{p}{\Sigma }}{A}_i\left(t\right)R_{\mathrm{ff}}(t-\mathrm{i\Delta t},\mathrm{j\Delta t}-\mathrm{i\Delta t}),$

时变模型系数L(t)用函数定义为：

$L\left(t\right)=\sqrt{\underset{i=0}{\overset{p}{\Sigma }}{A}_i\left(t\right)R_{\mathrm{ff}}(t-\mathrm{i\Delta t},-\mathrm{i\Delta t})}$

上式中，t＝1,2,3,…；f(t)为零均值非平稳随机过程；R_ff(t,jΔt)为 f(t)的在t时刻的相关函数；A_i(t)为时变自回归模型系数；L(t)为待定 时变模型系数；w(t)为方差为1的白噪声序列，且满足相关函数 $R_{\mathrm{ww}}(t,\mathrm{j\Delta t})=E\left[w\left(t\right)w(t-\mathrm{j\Delta t})\right]=\left(\begin{array}{cc}1,& j=0\\ 0,& j\ne 0\end{array}\right)$，p为模型阶数。

下击暴流风速具有很强的非平稳性，不仅平均风速具有时变特 征，而且脉动风速的功率谱也具有时变性，表现出明显的非平稳特征。 根据Priestley提出的进化谱理论，非平稳脉动风速的模拟，应该通过 一个与时间频率有关的非均匀调制函数对功率谱进行调制获得非平 稳风速的时变功率谱即进化谱，然后基于该进化谱对非平稳随机过程 进行模拟。然而，为了计算的简化，目前对于非平稳脉动风速的模拟， 通常将调制函数假定为仅与时间有关的函数，使得非平稳脉动风速的 模拟转化为时间函数对平稳脉动风速的调制，这并没有相关理论依 据。为此，本发明将根据进化谱理论，首先通过下击暴流的时变平均 风速模型和一个与时间频率有关的非均匀调制函数对功率谱进行调 制来获得下击暴流非平稳脉动风速时变功率谱，然后将平稳随机过程 的AR模型进一步考虑时变特征来建立非平稳随机过程的TAR时变 模型，最后通过时变功率谱、TAR时变模型实现下击暴流非平稳脉 动风速的模拟，通过模拟的下击暴流非平稳脉动风速的相关性、功率 谱分析验证了模拟的有效性。

附图说明

图1是本发明下击暴流非平稳脉动风速的数值模拟方法的步骤 流程图；

图2是本发明实施例中下击暴流径向风速与平移速度的矢量合 成图；

图3是本发明实施例中下击暴流时变平均风速图；

图4是本发明实施例中下击暴流的Kaimal非均匀调制函数图；

图5是本发明实施例中下击暴流风速时程和非平稳脉动成分图；

图6是本发明实施例中下击暴流非平稳脉动风速的时变功率谱 图；

图7是本发明实施例中单样本下击暴流非平稳脉动风速统计平 均估计谱与目标谱的对比图；

图8是本发明实施例中单样本下击暴流非平稳脉动风速统计平 均自相关函数与目标自相关函数的对比图；

图9是本发明实施例中单样本下击暴流非平稳脉动风速统计平 均自相关函数。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合 附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描 述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

本发明提供了一种下击暴流非平稳脉动风速的数值模拟方法，如 图1所示，包括以下步骤：

S1、根据进化谱理论，通过下击暴流的时变平均风速模型和一个 与时间频率有关的非均匀调制函数对功率谱进行调制来获得下击暴 流非平稳脉动风速时变功率谱

在步骤S1中，下击暴流时变平均风速具体描述为：

参照大气边界层中风速的描述方法，同时考虑到下击暴流的非平 稳性，在某一位置z高度t时刻的下击暴流风速可表示为随时间变化 的平均风速与脉动风速u(z,t)之和。与大气边界层近地风不同的 是，下击暴流风速时程具有较强的非平稳特性，下击暴流平均风速是 一种瞬时平均的时变过程。具有时变特征的下击暴流平均风速，可表 示为一个竖向分布函数与一个时间函数的乘积，即：

$\tilde{U}(z,t)=V\left(z\right)\times f\left(t\right)$

(1)

式中：f(t)为时间函数，其最大值为1；V(z)是最大平均风速的 竖向分布函数，Vicroy、Wood、OsegueraandBowles分别提出了3种 最大平均风速竖向风剖面^[14]。

Vicroy竖向分布模型

$V\left(z\right)=1.22\times [e^{-0.15z/z_{\max }}-e^{-3.2175z/z_{\max }}]\times V_{\max }$

(2)

Wood竖向分布模型

$V\left(z\right)=1.55{\left(\frac{z}{\delta }\right)}^{1/6}[1-\mathrm{erf}\left(0.7\frac{z}{\delta }\right)]\times V_{\max }$

(3)

OsegueraandBowles竖向分布模型

$V\left(z\right)=\left(\frac{\lambda R^2}{2r}\right)[1-e^{-{(r/R)}^2}](e^{-z/z^{*}}-e^{-z/ϵ})$

(4)

以上公式中，V_max为竖向分布风速中的最大风速；z_max为最大风 速所处的高度位置；erf为误差函数；δ为最大风速一半所处的高度 位置；r为距离下击暴流风场中心的径向位置；R为下击暴流风场的 辐射半径；z^*为大气边界层以外的特征高度；ε为大气边界层内的特 征高度；λ为尺度因子。

在碰撞射流理论中，径向风速是轴对称的，则在某一高度处径向 风速可以表示为^[6]：

$V_r\left(r\right)=\left(\begin{array}{cc}\Pi \times V_{r,\max }\times (r/r_{\max })& (0\le r<r_{\max })\\ \Pi \times V_{r,\max }\times \exp (-{[(r-r_{\max })/R_r]}^2)& (r\ge r_{\max })\end{array}\right)$

(5)

式中：V_r,max是风速场中某一高度处的最大风速；r_max是最大风速 点与下击暴流中心之间的水平距离；r是观测点与下击暴流中心之间 的距离；R_r是一个径向长度比例系数；Π为强度系数，表明雷暴强度 随时间的变化。考虑下击暴流自身在移动，所以观测点与雷暴中心的 距离r_t是随时间变化的，观测点的平均风速可表示为径向风速V_r与下 击暴流平移速度V_o的矢量合成，即

V_p(t)＝V_r(r_t)+V_ο

(6)

则时间函数f(t)可定义为空间点任意时刻的平均风速与最大平 均风速的比值，表示为竖向风剖面随时间变化的函数，即:

f(t)＝V_p(t)/max|V_p(t)|

(7)

在步骤S1中，下击暴流非平稳脉动风速具体为：

与大气边界层风速不同，下击暴流的平均风速具有明显的时变特 征，其脉动风速具有较强的非平稳特征。基于非平稳随机过程的进化 谱理论，下击暴流的非平稳脉动风速的模拟，可首先通过一调制函数 对某一功率谱进行调制获得时变功率谱，然后通过非平稳随机过程模 拟算法进行模拟。该模拟过程不仅需要获得有效的调制函数，而且还 需建立非平稳随机过程的模拟算法。

在步骤S1中，所述非均匀调制函数具体为：

根据非平稳随机过程的数值模拟理论，非平稳脉动风速的数值模 拟关键是获得时变功率谱。时变功率谱的获取可以通过对现场实测非 平稳风速进行时变谱估计，或根据进化谱理论通过调制函数对功率谱 进行调制来获得时变功率谱即进化谱。现场实测对于大多数的数值模 拟研究者来说并不可行，而更加方便可行的方法是采用功率谱调制。 为了获得有效的调制函数，李锦华等^[13]将非平稳脉动风速离散成若干 段在足够短时间Δt内可以近似为平稳脉动风速的短时时间序列，基于 Kaimal谱推导出与时间和频率有关的Kaimal谱非均匀调制函数，简 称Kaimal调制函数：

$A(\omega ,t)=\sqrt{\frac{{\tilde{U}}_j\left(t\right)}{U_j}{\left[\frac{1+50\frac{\omega z_j}{2\pi U_j}}{1+50\frac{\omega z_j}{2\pi {\tilde{U}}_j\left(t\right)}}\right]}^{5/3}}$

(8)

式中:ω为圆频率；z_j为空间某点的垂直地面高度；为空间 某点处的时变平均风速；为空间某点处非平稳风速的 统计平均风速。通过该方法，可进一步推导获得Davenport、Harris、 Simiu谱对应的非均匀调制函数。其中，Davenport调制函数：

$A(\omega ,t)={\left(\frac{1+{\left(\frac{1200\omega }{2\pi U_{10}}\right)}^2}{1+{\left(\frac{1200\omega }{2\pi {\tilde{U}}_{10}\left(t\right)}\right)}^2}\right)}^{2/3}$

(9)

Harris调制函数：

$A(\omega ,t)=\sqrt{\frac{{\tilde{U}}_{10}\left(t\right)}{U_{10}}{\left(\frac{2+{\left(\frac{1800\omega }{2\pi U_{10}}\right)}^2}{2+{\left(\frac{1800\omega }{2\pi {\tilde{U}}_{10}\left(t\right)}\right)}^2}\right)}^{5/6}}$

(10)

Simiu调制函数：

$\left(\begin{array}{cc}A(\omega ,t)=\sqrt{\frac{{\tilde{U}}_j\left(t\right)}{U_j}{\left[\frac{1+50\frac{\omega z_j}{2\pi U_j}}{1+50\frac{\omega z_j}{2\pi {\tilde{U}}_j\left(t\right)}}\right]}^{5/3}},& \frac{\omega z_j}{2\pi {\tilde{U}}_j\left(t\right)}\le 0.2\left(1\right)\\ A\left(t\right)={\left(\frac{{\tilde{U}}_j\left(t\right)}{U_j}\right)}^{4/3},& \frac{\omega z_j}{2\pi {\tilde{U}}_j\left(t\right)}>0.2\left(2\right)\end{array}\right);$

(11)

在步骤S1中，所述进化谱具体为：

根据Pristley的“进化谱”理论^[12]，一零均值非平稳随机过程f(t) 可表示为：

$f\left(t\right)={\int }_{-\infty }^{+\infty }A(\omega ,t)e^{\mathrm{i\omega t}}\mathrm{dZ}\left(\omega \right)$

(12)

式中，A(ω,t)为非均匀调制函数；Z(ω)为一正交过程，且满足

E[dZ(ω)]＝0

(13)

E[dZ^*(ω)dZ(ω′)]＝J(ω)δ(ω-ω′)dωdω′

(14)

其中*表示共轭，δ为狄拉克函数。则非平稳随机过程的均值为

$E\left[f\left(t\right)\right]={\int }_{-\infty }^{+\infty }A(\omega ,t)e^{\mathrm{i\omega t}}E\left[\mathrm{dZ}\left(\omega \right)\right]=0$

(15)

相关函数为：

$\left(\begin{array}{c}{R}_{\mathrm{ff}}(t,\tau )=E\left[f^{*}\left(t\right)f(t+\tau )\right]\\ ={\int }_{-\infty }^{+\infty }{\int }_{-\infty }^{+\infty }{A}^{*(\omega ,t)A({\omega }^{\prime },t+\tau )e^{-\mathrm{i\omega t}}{e}^{i{\omega }^{\prime }(t+\tau )}E\left[d{Z}^{*}\left({\omega }^{\prime }\right)\right]}\\ ={\int }_{-\infty }^{+\infty }{A}^{*}(\omega ,t)A\left(\text{ω,t+τ}\right)e^{\mathrm{i\omega \tau }}J\left(\omega \right)\mathrm{d\omega }\end{array}\right)$

(16)

当τ＝0时，

$E\left[f^2\left(t\right)\right]={\int }_{-\infty }^{+\infty }{\left|A(\omega ,t)\right|}^{2}J\left(\omega \right)\mathrm{d\omega }$

(17)

因此，进化谱S(ω,t)可通过一时频函数A(ω,t)对功率谱的非均匀调 制来表示^[16]，即

S(ω,t)＝|A(ω,t)|²J(ω)

(18)

根据公式(17)和公式(18)，可建立非平稳随机过程的相关函 数与时变功率谱之间的关系：

$\left(\begin{array}{c}{R}_{\mathrm{ff}}(t,\tau )=E\left[f^{*}\left(t\right)f(t+\tau )\right]\\ ={\int }_{-\infty }^{+\infty }{\left[S^{*}(\omega ,t)\right]}^{1/2}{\left[S(\omega ,t+\tau )\right]}^{1/2}{e}^{\mathrm{i\omega \tau }}\mathrm{d\omega }\end{array}\right)$

(19)

S2、将平稳随机过程的AR模型进一步考虑时变特征来建立非平 稳随机过程的TAR时变模型

在步骤S2中，非平稳脉动风速TAR时变模型的建立过程具体如 下：

在平稳随机过程的模拟中，较为成熟的方法主要有谱表示和线性 滤波方法。为了进一步建立非平稳随机过程的模拟方法，Liang等^[16] 详细推导了考虑时变特征的非平稳随机过程的谱表示方法。在平稳随 机过程模拟中，谱表示方法的模拟精度较高，但模拟效率较低。对于 实际大型工程的风场模拟，往往需要在谱表示方法基础上改进模拟效 率。将谱表示进一步考虑时变特征模拟非平稳随机过程，显然其模拟 效率更加难以满足实际工程的应用。线性滤波法由于计算量小、速度 快，已广泛用于随机过程的模拟中，但主要针对平稳高斯随机过程的 模拟。李锦华和李春祥^[17]考虑了随机过程的非高斯特性，实现了基于 线形滤波法中AR、ARMA模型的非高斯随机过程模拟。本节将继续 对线性滤波法中的AR模型进一步考虑时变特征，建立非平稳随机过 程的TAR时变模型。

设TAR时变模型阶数为p，非平稳随机过程模拟样本点数为N， 采样时间间隔为Δt，则基于TAR(p)时变模型的非平稳随机过程模 拟公式可表示为：

$f\left(t\right)=\underset{i=1}{\overset{p}{\Sigma }}{A}_i\left(t\right)f(t-\mathrm{i\Delta t})+L\left(t\right)w\left(t\right)$

(20)

式中：f(t)为零均值非平稳随机过程；A_i(t)为时变自回归模型系 数；L(t)为待定时变模型系数；w(t)为方差为1的白噪声序列，且满 足$R_{\mathrm{ww}}(t,\mathrm{j\Delta t})=E\left[w\left(t\right)w(t-\mathrm{j\Delta t})\right]=\left(\begin{array}{cc}1,& j=0\\ 0,& j\ne 0\end{array}\right)$。TAR时变模型的建立关键是 确定时变模型系数A_i(t)和L(t)。

在某时刻t＝t′时，有

$f\left(t^{\prime }\right)=\underset{i=1}{\overset{p}{\Sigma }}{A}_i\left(t^{\prime }\right)f(t^{\prime }-\mathrm{i\Delta t})+L\left(t^{\prime }\right)w\left(t^{\prime }\right)$

(21)

对于A_i(t′)的确定，将公式(21)两边同时右乘f(t′-jΔt)，并取数学 期望，有

$E\left[f\left(t^{\prime }\right)f(t^{\prime }-\mathrm{j\Delta t})\right]=E\left[\underset{i=1}{\overset{p}{\Sigma }}{A}_i\left(t^{\prime }\right)f(t^{\prime }-\mathrm{i\Delta t})f(t^{\prime }-\mathrm{j\Delta t})\right]+E\left[L\left(t^{\prime }\right)w\left(t^{\prime }\right)f(t^{\prime }-\mathrm{j\Delta t})\right]$

(22)

根据相关函数的定义可得，左项＝R_ff(t′,jΔt)；；右二项＝L(t′)R_wf(t′,jΔt)，即

$R_{\mathrm{ff}}(t^{\prime },\mathrm{j\Delta t})=\underset{i=1}{\overset{p}{\Sigma }}{A}_i\left(t^{\prime }\right)R_{\mathrm{ff}}(t^{\prime }-\mathrm{i\Delta t},\mathrm{j\Delta t}-\mathrm{i\Delta t})+L\left(t^{\prime }\right)R_{\mathrm{wf}}(t^{\prime },\mathrm{j\Delta t})$

(23)

对于右二项，R_wf(t′,jΔt)可理解为w(t′)与f(t′-jΔt)的互相关函数。 w(t)为该系统的输入，而f(t)为系统的输出。当前的输出只依赖于当 前和过去的输入，而与将来的输入无关，因此w(t′)与f(t′-jΔt)互相独 立，故

R_wf(t′,jΔt)＝0

(24)

将公式(24)代入公式(23)，可确定A_i(t′)

$R_{\mathrm{ff}}(t^{\prime },\mathrm{j\Delta t})=\underset{i=1}{\overset{p}{\Sigma }}{A}_i\left(t^{\prime }\right)R_{\mathrm{ff}}(t^{\prime }-\mathrm{i\Delta t},\mathrm{j\Delta t}-\mathrm{i\Delta t})$

(25)

其展开式为

对于L(t′)的确定，将公式(21)两边同时右乘w(t′)，并取数学期望，

有

$R_{\mathrm{fw}}(t^{\prime },0)=\underset{i=1}{\overset{p}{\Sigma }}{A}_i\left(t^{\prime }\right)R_{\mathrm{fw}}(t^{\prime }-\mathrm{i\Delta t},\mathrm{i\Delta t})+L\left(t^{\prime }\right)R_{\mathrm{ww}}(t^{\prime },0)$

(26)

当前输入w(t′)与过去输出f(t′-iΔt)互相独立，则

R_fw(t′-iΔt,iΔt)＝0

(27)

将公式(27)代入公式(26)，可得

R_fw(t′,0)＝L(t′)R_ww(t′,0)

(28)

根据白噪声的特性，有R_ww(t′,0)＝1，则

R_fw(t′,0)＝L(t′)

(29)

再一次将公式(21)两边同时右乘f(t′)，并取数学期望，有

$R_{\mathrm{ff}}(t^{\prime },0)=\underset{i=1}{\overset{p}{\Sigma }}{A}_i\left(t^{\prime }\right)R_{\mathrm{ff}}(t^{\prime }-\mathrm{i\Delta t},-\mathrm{i\Delta t})+L\left(t^{\prime }\right)R_{\mathrm{wf}}(t^{\prime },0)$

(30)

又因为

R_wf(t′,0)＝E[w(t′)f(t′)]＝E[f(t′)w(t′)]＝R_fw(t′,0)

(31)

将公式(29)和(31)代入公式(30)，有

$R_{\mathrm{ff}}(t^{\prime },0)=\underset{i=1}{\overset{p}{\Sigma }}{A}_i\left(t^{\prime }\right)R_{\mathrm{ff}}(t^{\prime }-\mathrm{i\Delta t},-\mathrm{i\Delta t})+L^2\left(t^{\prime }\right)$

(32)

令A₀(t′)＝-1，则

$L\left(t^{\prime }\right)=\sqrt{\underset{i=0}{\overset{p}{\Sigma }}{A}_i\left(t^{\prime }\right)R_{\mathrm{ff}}(t^{\prime }-\mathrm{i\Delta t},-\mathrm{i\Delta t})}$

(33)

因此，时变模型系数A_i(t′)和L(t′)，t′＝1,2,3,…，可分别根据公式 (25)和(33)来确定。

S3、通过所述时变功率谱、TAR时变模型对下击暴流非平稳脉 动风速进行有效模拟

在步骤S3中，下击暴流时变平均风速模拟具体为：

位于同一高度的6个下击暴流风速观测点P₁～P₆，如图2所示。 初始时刻，观测点与下击暴流中心距离均为r₀＝3500m，方位角分别为 θ₀₁＝0°、θ₀₂＝15°、θ₀₃＝30°、θ₀₄＝45°、θ₀₅＝60°、θ₀₆＝90°。考虑下击暴流自 身的移动，则某一高度观测点与雷暴中心的距离r_t随时间变化，观测 点的平均风速可表示为径向风速V_r与下击暴流平移速度V_o的矢量合 成：

$V_p\left(t\right)=\sqrt{V_r^2+V_o^2+2V_r·V_o\frac{r_0\cos \left({\theta }_0\right)-V_o·t}{\sqrt{{\left(r_0\sin \left({\theta }_0\right)\right)}^2+{(r_0\cos \left({\theta }_0\right)-V_o·t)}^2}}}$

(34)

下击暴流的平均风速竖向分布模型V(z)采用Vicroy模型，则下击 暴流时变平均风速可表示为：

$\tilde{U}(z,t)=V\left(z\right)·\frac{V_p\left(t\right)}{\max \left|V_p\left(t\right)\right|}$

(35)

考虑下击暴流竖向分布风速中的最大风速V_max＝80m/s，最大风速 所处的高度位置z_max＝67m；风速场中某一高度处径向最大风速V_r,max取 47m/s，与下击暴流中心之间的水平距离r_max取1000m，径向长度比例 系数R_r取700m，雷暴强度Π随时间的变化，可表示为：

$\Pi =\left(\begin{array}{cc}t/5& 0\le t\le 5\min \\ \exp [-(t-5)/11.542]& t>5\min \end{array}\right)$

(36)

下击暴流的平移速度取V_o＝8m/s。

各观测点的时变平均风速如图3所示，对于径向距离相同的各观 测点，时变平均风速的波动随着方位角的增大而减弱。从公式(8)～ (11)可知，非平稳脉动风速的功率谱非均匀调制函数与时变平均风 速有关，因此时变平均风速的波动性将进一步影响非平稳脉动风速功 率谱的时变性。考虑到观测点P1的时变平均风速具有较强的波动性， 本文将基于建立的TAR时变模型进一步对该位置处竖向不同点的非 平稳脉动风速进行模拟。

在步骤S3中，基于TAR时变模型的下击暴流非平稳脉动风速模 拟具体为：

现采用建立的TAR时变模型来模拟下击暴流行进路线上初始距 离雷暴中心3500米处的竖向(Z1＝15m，Z2＝25m，Z3＝45m)三个不 同点的非平稳脉动风速。模拟采用的上限截止圆频率取为4rad/s，模 拟的时间长度为800s，TAR时变模型阶数P＝15。下击暴流非平稳脉 动风速的时变功率谱采用Kaimal非均匀调制函数对Kaimal谱进行调 制。三个不同点处的下击暴流Kaimal非均匀调制函数，如图4所示。

不同点之间的相干函数采用Davenport相干函数模型：

$\mathrm{Coh}\left(\omega \right)=\exp (-\frac{\omega }{2\pi }\frac{C_z|z_i-z_j|}{\frac{1}{2}(U_i+U_j)})$

(37)

式中，z_i、z_j分别为第i、j点的竖向空间坐标，C_z为空间竖向两 点之间的衰减系数，U_i、U_j分别为第i、j点的平均风速，取为时变 平均风速的算术平均值。

基于TAR模型模拟的下击暴流风速时程和其脉动风速成分如图 5所示。从图5中可以发现脉动风速的振幅与时变平均风速大小有关， 即时变平均风速越大，脉动风速的振幅越大，这与实际风场特性相吻 合，也说明了非均匀调制函数的有效性。为了说明模拟的下击暴流脉 动风速具有目标非平稳特征，本文对模拟的100组下击暴流脉动风速 样本进行了时变功率谱估计如图6所示，样本估计的功率谱具有明显 的时变特征，且与目标时变谱的时变特征相吻合。

为了进一步说明模拟的有效性，图7、8分别展示了单样本下击 暴流非平稳脉动风速的统计平均功率谱、自相关函数与目标统计平均 功率谱${\bar{S}}_{\mathrm{ff}}\left(\omega \right)=\frac{1}{T}{\int }_0^T{S}_{\mathrm{ff}}(\omega ,t)\mathrm{dt}$、自相关函数${\bar{R}}_{\mathrm{ff}}\left(\tau \right)=\frac{1}{T}{\int }_0^T{R}_{\mathrm{ff}}(t,\tau )\mathrm{dt}$的对比， 单样本统计平均功率谱、自相关函数与目标均相互吻合。此外，在图 5中，不同高度的下击暴流非平稳脉动风速的波动具有明显的相关性， 且相关性随着不同位置之间的距离增大而减弱如图9所示，这与实际 风场特性相吻合。

参考文献

[1] 李春祥,李锦华,于志强.输电塔线体系抗风设计理论与发展[J]. 振动与冲击,2009,28(10):15-25.

[2] Chay MT,Letchford CW.Pressure distributions on a cube in a  simulated thunderstorm downburst-Part A:Stationary downburst  observations[J].Journal of Wind Engineering and Industrial  Aerodynamics,2002,90:711-732.

[3] Letchford C W,Mans C,Chay MT.Thunderstorms--their  importance in wind engineering(a case for the next generation wind  tunnel)[J].Journal of Wind Engineering and Industrial  Aerodynamics,2002,90(12-15):1415-1433.

[4] Holmes JD.Modeling of extreme thunderstorm winds for wind  loading of structures and risk assessment[C].Proceedings of 10th  International Conference on Wind Engineering,Copenhagen, Denmark,21-24June1999:1409-1415.

[5] Anindya Sengupta,Fred L.Haan,Partha P.Sarkar,Vasanth  Balaramudu.Transient loads on buildings and tornado winds in  microburst[J].Journal of Wind Engineering and Industrial  Aerodynamics,2008,96(10-11):2173-2187.

[6] Holmes JD,Oliver SE.An empirical model of a downburst[J]. Engineering Structures,2000,22:1167-1172.

[7] Wood GS,Kwok KCS,Motteram NA,Flethcer DF.Physical and  numerical modelling of thunderstorm downdraft[J].Journal of Wind  Engineering and Industrial Aerodynamics,2001,89(6):535-552.

[8] Savory E,Gerard ARP,Zeinoddini M,Toy N,Disney P.Modelling  of tornado and microburst-induced wind loading and failure of a  lattice transmission tower[J].Engineering Structures,2001, 23:365-75.

[9] 瞿伟廉,王锦文.下击暴流风荷载的数值模拟[J].武汉理工大学 学报,2008,30(2):70-74.

[10] 李春祥,刘晨哲,申建红,李锦华.土木工程下击暴流风速数值 模拟的研究[J].振动与冲击,2010,29(12):49-54.

[11] 张文福,谢丹,刘迎春,计静.下击暴流空间相关性风场模拟[J]. 振动与冲击,2013,32(10):12-16.

[12] Priestley MB.Power spectral analysis of non-stationary random  processes[J].Journal of Sound and Vibration,1967,6:86-97.

[13] 李锦华,李春祥,申建红.非平稳脉动风速的数值模拟[J].振动 与冲击,2009,28(1):18-23.

[14] Chen L,Letchford CW.A deterministic-stochastic hybrid model of  downbursts and its impact on a cantilever structure[J].Engineering  Structures,2004,26:619-26.

[15] Vicroy DD.Assessment of microburst models for downdraft  estimation[J].Journal of Aircraft,1992,29:1043-1050.

[16] Liang J,Chaudhuri SR,Shinozuka M.Simulation ofnon-stationary  stochastic processes by spectral representation[J].Journal of  Engineering Mechanics ASCE,2007,133(6):616-627.

[17] Li Jinhua,Li Chunxiang.Simulation of Non-Gaussian Stochastic  Process with Target Power Spectral Density and Lower-order  Moments[J].Journal of Engineering Mechanics ASCE,2012, 138(5):391-404.

[18] Chay MT,Albermani F,Wilson R.Numerical and analytical  simulation of downburst wind loads[J].Engineering Structures, 2006,28(2):240-254.

相比于现有技术的缺点和不足，本发明具有以下有益效果：为了 简化，现有技术通常将非均匀调制函数假定为仅与时间有关的均匀调 制函数，将非平稳随机过程的模拟转化为一均匀调制函数对平稳随机 过程的调制，回避了根据时变功率谱模拟非平稳随机过程的难点。而 且，对于非平稳脉动风速的模拟，通常将时变平均风速函数假定为平 稳脉动风速的均匀调制函数并没有相关理论依据。本发明没有回避根 据时变功率谱模拟非平稳随机过程的难点问题，首先通过下击暴流的 时变平均风速模型和一个与时间频率有关的非均匀调制函数对功率 谱进行调制来获得下击暴流非平稳脉动风速时变功率谱；然后将用于 模拟平稳随机过程的AR模型进一步考虑时变特征来建立模拟非平稳 随机过程的TAR时变模型，从而解决根据下击暴流时变功率谱模拟 非平稳随机过程的难点问题，最终实现下击暴流非平稳脉动风速的有 效模拟。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明， 凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等， 均应包含在本发明的保护范围之内。