吴红华，米慧敏

(湖南大学 建筑安全与节能教育部重点试验室，湖南 长沙 410082)

非高斯脉动风压的分形模拟研究*

吴红华†，米慧敏

(湖南大学 建筑安全与节能教育部重点试验室，湖南 长沙 410082)

基于风洞试验中采集的厦门沿海某高层建筑模型的表面风压时程,计算各测点脉动风压的分形维数,偏度系数和峰度系数,找出脉动风压呈现较明显的非高斯特性的4个测点.考虑脉动风时程具有分形的特性,结合相关性变形法,采用具有分形特性的Weierstrass-Mandelbrot函数对非高斯特性明显的4个测点的脉动风压进行了模拟.结果表明,模拟得到的4个测点的非高斯脉动风压时程与风洞试验的非高斯脉动风压时程的概率分布曲线吻合较好,功率谱图形比较相似,分形维数,偏度系数和峰度系数均比较接近.说明本文所提出的用分形方法模拟具有非高斯特性的脉动风压是可行的,弥补了传统的线性滤波器法和谐波叠加法模拟结果均不具有分形特性的缺点,对非高斯脉动风压的数值模拟有一定的参考价值.

分形;Weierstrass-Mandelbrot函数;非高斯脉动风压;模拟

20世纪70年代，分形理论的提出，揭示了非线性系统中有序与无序、确定性与随机性的统一，分形维数成为定性或者定量分析不规则的几何体和复杂现象的工具[1].随着分形方法和理论的发展，Mandelbrot[2]认识到Weierstrass-Mandelbrot函数具有分形特性，且在湍流分析过程中，融入了分形理论.Pinnington[3]分别利用随机型WM函数和确定型WM函数建立模型模拟路面的振动，得到了随机型WM函数的模拟结果与实测结果更接近的结论.

刘磊,West等[4-5]指出，自然风速脉动在功率谱密度中高频部分呈现出幂指数为-5/3的幂函数形式，称为“-5/3”律.脉动风速时程具有自相似分形特征[6-7]，反映了时间序列中的关联特性[8]，是风速脉动的一个重要特征.自相似分形和脉动风速的“-5/3”律、间歇性特征[9]均具有紧密联系.通过计算全国80个相似气象站点测量的风速时程的风速时程的分形维数，钟莉等[10]发现80个站点测得的风速时程的分形维数都比较接近.

一般情况下，常假定脉动风压为平稳高斯随机过程，但是大量研究证明，这种假设在很多情况下是不成立的.气流的分离、再附和涡脱落，使得建筑物表面区域风压呈现较复杂的变化[11-12].Holmes[13]的研究表明风压场中有非高斯特性的存在.Tsutsumi 等[14]通过风洞试验得出矩形建筑结构表面脉动风压有明显非高斯特性的结论，并分析了影响该特性的因素.

静态转换法(Static Transformation Methods)对高斯随机过程进行非线性变换，生成满足一定要求的非高斯随机过程.静态转换法分为相关函数转换法和谱更新法.目前，很多学者对非高斯脉动风压的模拟进行研究，但是并未考虑分形特性对脉动风压的影响.

本文结合相关性变形法对非高斯脉动风压进行了分形盒维数分析，进而应用随机型函数对于非高斯脉动风压时程进行了分形模拟，并且从脉动风压的概率分布、分形盒维数、偏度、峰度和功率谱5个方面对于模拟脉动风压和风洞试验脉动风压进行了对比，发现两者比较接近，说明了采用分形方法模拟非高斯脉动风压具有可行性.

1 基本方法

1.1 Weierstrass-Mandelbrot函数

Mandelbrot最早发现Weierstrass-Mandelbrot函数(WM函数)具有分形特性，将其应用于湍流中粒子的运动轨迹研究中.

WM函数的表达式为：

(1)

式中：t为频域分析区间内的时间参数；ω为频率或波数，ω=bn，b为实数，n为求和计数量.

Mandelbrot[2]提出了不具有特殊长度标度的WM函数：

(2)

式中：φn为独立的位相因子，且每个φn对应一个函数u(t)；D为时间序列的盒维数.式(2)所定义的函数即是著名的Weierstrass-Mandelbrot 函数.实际应用中，根据式(2)中φn的取值方案，WM函数可分为确定型WM函数和随机型WM函数.

1.2 分形盒维数D

分形维数作为分形理论中判断系统或事物的复杂程度、局部与整体的自相似性的度量工具，包括Hausdorff维数、相似维数、盒维数(又称计盒维数)、关联维数等.其中，盒维数因具有计算简单和适用范围广的优点，得到研究学者的广泛青睐[15].

设在平面R2内风速时程曲线为L，将R×R划分为尽可能小的网格kδ1×kδ2(k=1,2,3…，为网格放大倍数).与L相交的网格数量为Nkδ1(或Nkδ2)，则分形盒维数D的数学表达式为[1]：

(3)

D实际为无标度区内(-logkδi,logNkδi)(i=1或2)的双对数图的斜率.其中k1≤k≤k2，对应与网格尺寸范围为((k1δ1×k1δ2)～(k2δ1×k2δ2)).在该无标度区间内对应的(-logkδi，logNkδi) (i=1或2 )点对数量为k2-k1+1.依据上述分析及分形盒维数的定义，可知，在无标度区间内，-logkδi与logNkδi满足线性回归方程：

logNkδi=-Dlog(kδi)+b，i=1或2

(4)

盒维数D由-logkδi，logNkδi的关系唯一确定，用最小二乘法即可求得D：

k1≤k≤k2，i=1或2，1

(5)

1.3 偏度系数和峰度系数

通常采用多阶统计矩来描述具有非高斯特性信号的概率密度函数特征[16].

3阶统计量偏度系数(Wsk)的数学表达式为公式(6)，描述随机过程的分布偏离高斯分布的偏斜程度，高斯分布中，Wsk=0.当偏度系数小于0时，该随机过程为左偏态；反之，为右偏态.

4阶统计量峰度系数(Wku)的数学表达式为公式(7)，描述随机过程分布相比高斯分布的尖削或平坦程度，高斯分布中，Wku=3.当峰度系数小于3时，该随机过程概率密度函数曲线与高斯分布相比相对平坦，视为负的峰度；反之，为正峰态.所以，偏度系数、峰度系数从不同角度体现了随机过程的分布偏离高斯分布的程度.

(6)

(7)

1.4 脉动风压的功率谱

对于高层建筑，其侧风面脉动风压主要是由来流紊流和尾流涡脱构成，准定常假定不适用.Ohkuma与Kanaya[17]在考虑矩形高层建筑的来流紊流和涡激励对横风向激励贡献的基础上，给出了横风向风激励谱的表达式：

(8)

式中：η=fD/UH，为无量纲频率;St=0.135-0.069exp (-0.056H/D)，为斯托罗哈数;β=0.6exp (-0.3H/D)，为带宽系数.

金虎等[18]在Ohkuma与Kanaya给出的横风向风激励谱数学模型的基础上提出了横风向脉动风压谱的表达式(9)，并证明侧风面风压功率谱试验结果与式(9)吻合较好：

(9)

式中:A，B，St为待定参数；η=fD/UH，为无量纲频率;Sw(f)为横风向脉动风压谱.

采用Welch谱法[19]计算试验脉动风压的功率谱，具体计算方法如下.

将脉动风压时程x(n)的长度N分成L段，每段M个数据，则第i段的功率谱为：

(10)

对L分段周期进行平均，可得到整个信号xN(n)的功率谱:

(11)

1.5 脉动风压的互相关函数

x(t)是各态历经过程的一个样本函数，设X(iω)是x(t)的傅里叶变换，即：

(12)

(13)

用t+τ代替公式(13)中的t：

(14)

将自相关函数写成如下频率函数的形式：

(15)

公式(15)表明，自相关函数Rx(τ)描述了随机样本x(t)在任意两个时刻t1,t2(t1≠t2)的相关程度，且Rx(τ)是Sx(ω)的傅里叶逆变换.

根据式(15)，可以得到两个随机样本x(t)和y(t)的互相关函数Rxy(τ)：

(16)

同理，Rxy(τ)描述了随机样本x(t)和y(t)在任意两个时刻t1,t2(t1≠t2)的相关程度,Rxy(τ)是Sxy(ω)的傅里叶逆变换.

1.6 多节点非高斯脉动风压的模拟

0

(17)

吴红华等[19]针对公式存在的模拟周期无限长导致的模拟效率大大降低的缺点，同时考虑脉动风具有分形特性这一特点，对式(17)进行了改进，改进后的公式如下：

(18)

对多节点的脉动风压进行模拟时，互谱密度函数矩阵为：

(19)

第i个与第j个待模拟点的脉动风压互谱密度函数表达式如下：

(20)

(21)

coh(ω)=

(22)

对互谱密度函数矩阵Su(ω)进行Cholesky分解：

(23)

式中：H*T(ω)为H(ω)的复数共轭转置矩阵.

现就非高斯脉动风压的模拟将分形特性融入到Gurley提出的相关性变形法，以生成具有分形特性的非高斯脉动风压时间序列.

相关性变形法通过如下公式将高斯过程转换为非高斯过程：

(24)

(25)

式中：Rw(τ)为高斯脉动风压时程的相关函数；Rnw(τ)为非高斯脉动风压时程的相关函数；wm(t)为高斯脉动风压时程；wnm(t)为非高斯脉动风压时程；

μ4=γ4-3，γ4为峰度系数.

因此，考虑分形特性的相关性变形法模拟多点非高斯脉动风压时程的过程如下：

2 非高斯脉动风压时程的模拟

2.1 非高斯脉动风压时程的选取

本次试验在湖南大学风洞实验室的HD-3 大气边界层风洞中进行，试验段长为10m，风速为0～20m/s且连续可调，截面宽为2.5m， 高为3m，转盘直径为1.8m.实验模型为厦门沿海某高层建筑，测点沿模型高度分层布置，除3层布置22个测点外，其余层布置15个测点，测点布置如图1所示.

图1 测点布置图Fig.1 The arrangement of measuring points

本次试验采用A类风场.试验的几何比尺、风速比尺和时间比尺分别为1/200，1/8 和 1/25.风压数据采样频率为312.5Hz，试验中共进行80次重复独立采样，每个样本的数据均为10 000个，对应于实际采样时间 800s.对采集的风压时程样本取平均值得到平均风压，样本数据减掉平均风压即得到所需的脉动风压时程.

楼文娟[22]通过研究认为，建筑物非高斯风压区别于高斯风压的标准为Wsk>0.2且Wku>3.5，在0°风向角下，迎风面几乎都是高斯区；左、右侧面则存在高斯区与非高斯区，高斯区主要位于下风区，非高斯区主要位于前缘气流分离区和近地面区域；背风面大部分区域属于非高斯区域.

在分析大跨屋盖结构时，孙瑛[23]得到非高斯的标准为|Wsk|≥0.2且Wku≥3.7.在对具有非高斯特性的高层建筑的研究中，韩宁等[24]认为区分高斯非高斯的标准为|Wsk|≥0.25且Wku≥3.2或|Wsk|≥0.45或Wku≥4.0.

至此，通过MATLAB分析0°风向角下各测点脉动风压时程的偏度系数和峰度系数，结合对已有研究的总结，选取实际高度为97.750m的E1，E2和121.150m的L16,L17四个测点的脉动风压时程数据为模拟非高斯脉动风压的样本数据，其中，E1，E2测点位于建筑物的右侧面，L16，L17位于建筑物的左侧面.E1，E2，L16，L17的概率密度图如图2所示.

图2 试验脉动风压概率分布与高斯分布对比Fig.2 The comparison among probability distribution of fluctuating wind pressure by measuring and Gaussian distribution

根据1.4节所述方法计算试验脉动风压的功率谱，并且调用MATLAB的非线性拟合函数Lsqcurvefit，选择恰当的A，B，St将试验脉动风压功率谱与式(9)相拟合，为模拟非高斯脉动风压时程选择恰当的指定风压功率谱.4个测点的脉动风压功率谱与拟合结果的对比如图3所示.

由图3可以看出，通过选择不同的A，B，St，可使式(9)与4个试验脉动风压的功率谱曲线之间达到很好的拟合.本文分别将E1和E2测点A，B，St的平均值作为E1和E2的指定功率谱的A，B，St；将L16和L17测点A，B，St的平均值作为L16和L17的指定功率谱的A，B，St.

(a)E1测点 (b)E2测点

(c)L16测点 (d)L17测点图3 4个测点脉动风压功率谱Fig.3 The auto-power spectrum of fluctuating wind pressure

由4个测点脉动风压时程的偏度、峰度系数及图3可以看出，其概率分布具有显著的非高斯特性，可以将4个测点作为非高斯脉动风压的模拟点.E1,E2,L16和L17的试验脉动风压时程分别如图4(a),4(c),4(e)和4(g)所示.

2.2 脉动风压的分形特征分析

根据选取的风压时程数据，按照公式(5)计算得出各个测点的脉动风压时程的盒维数，计算结果如表1所示.

表1 试验脉动风压D，Wsk和Wku

2.3 非高斯脉动风压时程的模拟

根据前述非高斯脉动风压的模拟过程，对E1，E2，L16和L17四个测点的脉动风压分别进行分形模拟，模拟得出的非高斯脉动风压时程分别如图4(b),4(d),4(f),4(h)所示.

(a) E1测点试验脉动风压时程 (b) E1测点模拟脉动风压时程 (c) E2测点试验脉动风压时程

(d) E2测点模拟脉动风压时程 (e) L16测点试验脉动风压时程 (f) L16测点模拟脉动风压时程

(g) L17测点试验脉动风压时程 (h) L17测点模拟脉动风压时程图4 脉动风压时程Fig.4 The fluctuating wind pressure

3 试验与模拟脉动风压时程的对比

3.1 模拟与试验脉动风压时程的统计特性对比

分别将E1，E2，L16和L17四个测点模拟脉动风压与试验脉动风压的概率分布及高斯曲线进行对比.对于高斯分布曲线，由于采用试验脉动风压时程进行统计计算，其标准差采用试验脉动风压的标准差，如图5所示.

由图5可看出，采用本文方法模拟的脉动风压的概率分布与试验脉动风压的概率分布比较接近.

3.2 模拟脉动风压功率谱与指定谱的对比

4个模拟点的模拟脉动风压功率谱与指定谱相比较的结果如图6所示.

从图6可以看出，模拟的脉动风压功率谱与指定谱比较接近，因此可以认为模拟的脉动风压功率谱满足指定功率谱的要求.

(a)E1测点

(b) E2测点

(c) L16测点

(d) L17测点图5 实验脉动风压与模拟脉动风压概率分布、高斯分布对比Fig.5 The comparison among the probability distribution of fluctuating wind pressure by measuring and simulating and Gauss distribution

(a)E1测点

(b) E2测点

(c) L16测点

(d) L17测点图6 脉动风压功率谱与指定谱对比Fig.6 The comparison between auto-power spectrum of simulated fluctuating wind pressure and specified power spectrum

3.3 模拟与试验脉动风压时程的分形盒维数、偏度系数和峰度系数的对比

对通过分形模拟得出的各个测点的脉动风压时程的分形盒维数、偏度系数和峰度系数进行计算，具体结果如表2所示.

表2 模拟脉动风压D,Wsk和Wku

对比表2与表1对应测点的分形盒维数、偏度系数和峰度系数，发现它们的数值都较接近，证明采用具有分形特性的相关性变形法模拟非高斯脉动风压是可行的.

4 结 论

对试验得到的样本进行非高斯特性分析，考虑分形特性，运用相关性变形方法进行非高斯脉动风压时程的模拟.通过分析研究得到以下结论：

1)在风洞试验脉动风压时程的基础上，利用随机型函数模拟出4个测点具有分形特性的非高斯脉动风压时程.弥补了传统的线性滤波器法和谐波叠加法模拟结果均不具有自相似分形特性的缺点.

2)计算了4个测点的风洞试验非高斯脉动风压时程和模拟非高斯脉动风压时程的分形盒维数，通过比较，可以看出数值很接近，说明非高斯脉动风压时程的分形模拟是可行的.

3)通过文中方法模拟得到的非高斯脉动风压时程的概率分布、偏度系数、峰度系数和脉动风压功率谱与目标概率分布、偏度系数、峰度系数和脉动风压功率谱比较接近，说明相关性变形法可以对非高斯脉动风压进行模拟.

[1] 张济忠. 分形[M].2版.北京:清华大学出版,2011: 30-45.

ZHANG Jizhong. Fractal [M].2nd ed. Beijing: Tsinghua University Press,2011: 30-45.(In Chinese)

[2] MANDELBROT B B. The fractal geometry of nature [M]. San Francisco:Macmillan Publishers,1982: 466-469.

[3] PINNINGTON R J. A particle-envelope surface model for road-tyre interaction International[J].Journal of Solids and Structures,2012,49(3/4): 546-555.

[4] 刘磊,胡非,李军,等.基于Weierstrass-Mandelbrot函数的分形风速脉动仿真[J].气候与环境研究,2013,18(1): 43-50.

LIU Lei,HU Fei,LI Jun,etal. On the use of Weierstrass-Mandelbrot function to simulate fractal wind fluctuations [J]. Climatic and Environmental Research,2013,18(1): 43-50. (In Chinese)

[5] WEST B J,SHLESINGER M. The noise in natural phenomena [J]. American Scientist,1990,78(1): 40-45.

[6] SCHMITT E,SCHERTZER D,LOVEJOY S,etal. Estima-tion of universal for atmospheric turbulent multifractal indices for atmospheric turbulent velocity fields [J]. Fractals,1993,l(3) :568-575.

[7] SYU C Y,KIRCHHOFF R H. The fractal dimension of the wind[J]. Journal of Solar Energy Engineering,1993,l15(3) :151-154.

[8] FEDER J. Fractals [M]. New York: Plenum,1988:149-162.

[9] MANDELBROT B B. Intermittent turbulence in self-similar cascades divergence of high moments and dimension of the carrier [J]. Journal of Fluid Mechanics,1974,62: 331-358.

[10]钟莉,李杰. 风速时程的分形特征分析[J]. 华中科技大学学报:城市科学版,2008,25(4):273-275.

ZHONG Li,LI Jie. Fractal characteristic analysis of wind speed time history [J]. Journal of Huazhong University of Science and Technology: Urban Science,2008,25(4): 273-275. (In Chinese)

[11]沈国辉,钱涛,何国军,等. 风帆体型建筑表面风压的分布特征研究[J]. 湖南大学学报:自然科学版,2013,40(9): 1-6.

SHEN Guohui,QIAN Tao,HE Guojun,etal. Investigation of wind pressure distribution on a sail-shaped building [J]. Journal of Hunan University: Natural Sciences,2013,40(9): 1-6. (In Chinese)

[12]李秋胜,李永贵,郅伦海. 典型高层住宅建筑风压分布特性的试验研究[J]. 湖南大学学报:自然科学版,2011,38(4): 14-19.

LI Qiusheng,LI Yonggui,ZHI Lunhai. Experimental investigation of the wind pressure distributions on a typical tall residential building[J]. Journal of Hunan University: Natural Sciences,2011,38(4): 14-19. (In Chinese)

[13]HOLMES J D. Non-Gaussian characteristics of wind pressure fluctuations[J].Journal of Wind Engineering Industrial Aerodynamics,1981,7: 103-108.

[14]TSUTSUMI J,KATAYAMA T,NISHIDA M. Wind tunnel tests of wind pressure on regularly aligned buildings[J]. Journal of Wind Engineering and Industrial Aerodynamics,1992,43(1/3):1799-1810.

[15]PEITGEN H O,JURGENS H,SAUPE D. Chaos and fractals: new frontiers of science [M]. 2nd ed. Beijing: National Defense Industry Press,2008: 192-209.

[16]邱天爽,张旭秀,李小兵,等. 统计信号处理—非高斯信号处理及其应用[M].北京:电子工业出版社,2004: 114-119.

QIU Tianshuang,ZHANG Xuxiu,LI Xiaobing,etal. Statistical signal processing—non-Gaussian signal processing and its application [M]. Beijing: Publishing House of Electronics Industry,2004: 114-119. (In Chinese)

[17]OHKUMA T,KANAYA A. On the correlation between the shape of rectangular cylinders and characteristics of fluctuating lifts on them [C] // Proceedings of 5th Symposium on Wind Effects on Structures. Tokyo,Japan,1978: 147-154.

[18]金虎,楼文娟,沈国辉. 复杂体型超高层建筑风压脉动特性[J]. 哈尔滨工业大学学报,2009,41(8): 111-117.

JIN Hu,LOU Wenjuan,SHEN Guohui. Fluctuating character investigation on wind pressure of high-rise building with complex shape [J]. Journal of Harbin Institute of Technology,2009,41(8): 111-117. (In Chinese)

[19]吴红华,徐阳,呼腊梅,等.近地面脉动风速时程的分形模拟研究[J].地震工程与工程震动,2015,35(4): 121-129.

WU Honghua,XU Yang,HU Lamei,etal. Research on fractal simulation of subsurface fluctuating wind speed time-history [J]. Earthquake Engineering and Engineering Vibration,2015,35(4): 121-129. (In Chinese)

[20]PATZOLD M,KILLAT U,LAUE F. Adeterministic digital simulation model for suzuki processes with application to a shadowed rayleigh land mobile radio channel[J]. IEEE Transactions on Vehicular Technology,1996,45(2): 318- 331.

[21]DAVENPORT A G. The relationship of wind structure to wind loading[C]//Proceedings of the Symposium on Wind Effects on Buildings and Structures. London: National Physical Laboratory,1965,1: 53-102.

[22]楼文娟. 超高层建筑脉动风压的非高斯特性[J]. 浙江大学学报:工学版,2011,45(4): 671-677.

LOU Wenjuan. Non-Gaussian feature of wind-induced pressure on super-tall building[J].Journal of Zhejiang University: Engineering Science,2011,45(4): 671-677. (In Chinese)

[23]孙瑛. 大跨屋盖结构风压脉动的非高斯特性[J]. 土木工程学报,2007,40(4): 1-5.

SUN Ying. Non-Gaussian features of fluctuating wind pressures on long span roofs [J]. China Civil Engineering Journal,2007,40 (4): 1-5. (In Chinese)

[24]韩宁,顾明.方形高层建筑风压脉动非高斯特性分析[J]. 同济大学学报:自然科学版,2012,40(7): 971-976.

HAN Ning,GU Ming. Analysis on non-Gaussian features of fluctuating wind pressures on square tall buildings [J]. Journal of Tongji University:Natural Sciences,2012,40(7): 971-976. (In Chinese)

Research on Fractal Simulation of Non-GaussianFluctuating Wind Pressure

WU Honghua†,MI Huimin

(Key Laboratory of Building Safety and Efficiency of the Ministry of Education,Hunan University,Changsha 410082,China)

Based on the time history of wind pressure collected from the wind tunnel tests on a high-rise building at Xiamen coast，four points with obvious non-Gaussian characteristics were found out by the calculation of fractal dimension，skewness and kurtosis.Considering the fractal characteristics of time history of fluctuating wind and the correlation distortion method，the non-Gaussian fluctuating wind pressure was simulated by adopting fractal Weierstrass-Mandelbrot function.Comparison of simulation results with experimental results in the aspects of probability distribution，power spectrum，fractal box dimension，skewness and kurtosis，showed better consistency.Therefore，it is feasible to simulate the fluctuating wind pressure with non-Gaussian characteristics by using fractal method.The method makes up for the shortcomings of traditional linear filter and superposition of harmonic，and provides a reference for the simulation of non-Gaussian fluctuating wind pressure.

fractals; Weierstrass-Mandelbrot function; non-Gaussian fluctuating wind pressure; simulation

1674-2974(2017)07-0059-10

10.16339/j.cnki.hdxbzkb.2017.07.008

2016-05-09

国家自然科学基金资助项目(51278190，51478179)，Natural Science Foundation of China(51278190，51478179)；浙江省自然科学基金资助项目(LY12E08010),Natural Science Foundation of Zhejiang Province(LY12E08010)

吴红华(1967—)，女，湖北武汉人，湖南大学副教授，工学硕士†通讯联系人，E-mail: wuliyy@126.com

TU312

A

Davenport提出的空间任意两点的相关函数[21]来表示脉动风压的空间相关性，相关函数可按照公式(22)计算：