脉动风速连续随机场的降维模拟

刘章军1,2,叶永友2,刘增辉2

(1.武汉工程大学土木工程与建筑学院,湖北,武汉 430074;2.三峡大学土木与建筑学院,湖北,宜昌 443002)

摘 要:基于标准正交随机变量的波数谱表示,通过定义标准正交随机变量集的随机函数形式,建立了连续时空随机场模拟的波数谱-随机函数方法。同时,引入快速傅里叶变换(FFT)的算法,极大地提高了波数谱-随机函数方法的模拟效率。在波数谱-随机函数模拟方法中,仅需两个基本随机变量即可在概率密度层次上描述时空随机场的概率特性,并利用数论方法选取基本随机变量的代表性点集,实现对连续时空随机场模拟的降维表达。数值算例表明,当模拟相同数量的样本时,综合考虑模拟的效率和精度两方面,该文方法与传统的波数谱表示方法不分伯仲,但该文方法所需的基本随机变量最少,生成的代表性样本数量少且构成一个完备的概率集,从而可结合概率密度演化理论实现结构随机动力反应及动力可靠度的精细化分析。最后,结合Kaimal风速谱及Davenport空间相干函数模型,模拟了水平向脉动风速连续随机场,验证了该文方法的有效性和优越性。

关键词:脉动风速;时空随机场;频率-波数谱;随机函数;FFT算法;降维

对于高耸建筑和大跨度桥梁,风荷载在结构设计中起控制作用。作用于结构上的风荷载包括平均风和脉动风,其中脉动风荷载是引起结构动力响应的主要因素[1],因此脉动风速随机场模拟是结构抗风分析的首要工作。

脉动风速随机场模拟最常用的两种方法是谱表示(SRM)[2―4]和本征正交分解(POD)[5―6]。谱表示最早由 Shinozuka等[3]提出并被应用于多维多变量随机场的模拟,Yang[7]在利用谱表示模拟多变量随机过程时,首次引入FFT算法,极大地提高了谱表示的模拟效率。随后,Deodatis[8]在Yang等的研究基础上,在谱表示中引入双索引频率,保证了模拟结果的各态历经性,这一改进使得谱表示成为随机风场模拟的经典方法。与谱表示相对应,Paolao[5―6]提出基于功率谱密度矩阵的本征正交分解,并将其用于模拟脉动风速随机场。与谱表示相比,本征正交分解物理意义明确,且可以采用模态截断技术,从而仅用少数几阶本征模态即可近似模拟脉动风速随机场[9]。上述谱表示和本征正交分解均是基于功率谱密度矩阵的谱分解来实现对脉动风速随机场或者随机向量过程的模拟,因而文献[10]对上述两种方法进行了统一表达。对于大跨度桥梁及超高层建筑等复杂工程结构,其脉动风速随机场的功率谱密度矩阵为高阶矩阵,对其进行分解时往往需要大量的计算。为此,Deodatis等[11]基于随机场的频率-波数谱,提出了连续时空随机场的波数谱表达,将时空随机场表达为连续的随机波形式;Benowitz等[12]进一步结合FFT算法,将该方法应用于脉动风速连续时空随机场的模拟。此外,Carassale等[13]基于脉动风速互功率谱密度函数,发展了一类脉动风速时空随机场模拟的连续 POD方法,并给出了特征问题的半解析解,该方法仅适用于特定的随机风场。

然而,在上述方法中,无论是谱表示和POD,还是波数谱,在模拟时空随机场时,均需要大量的随机变量才能保证模拟结果的精度,且模拟结果不能在概率密度函数层次上反映随机风场的概率特性。近年来,笔者提出了随机函数的降维思想,并与随机过程模拟的谱表示及 Karhunen-Loeve展开相结合[14―16],实现了仅用1~2个基本随机变量即可精细地模拟一维单变量随机过程,极大地减少了随机过程模拟的计算量。

为此,本文将在脉动风速时空随机场的波数谱表示基础上,引入随机函数的约束形式,实现脉动风速连续时空随机场的降维模拟,并利用FFT算法极大地提高模拟的计算效率,从而为脉动风速随机场模拟提供一种更为精细高效的方法。

1 时空随机场的波数谱表达

f0(x,t)是一个具有零均值的1V-2D实值时空随机场,在时间t上平稳,在空间x上均匀,其单边波数谱为G0(κ,ω),则存在两个实的正交增量过程u(κ,ω)和v(κ,ω),将时空随机场f0(x,t)表示为如下积分形式[17]

式中:x/m为空间变量;κ/m-1为相应的波数;t/s为时间变量;ω/s-1为相应的圆频率。

定义在-∞<κ<∞ 和 0<ω<∞ 上的增量=满足如下基本条件:

式中:E[·]为数学期望;为 Kronecker符号。

于是,式(1)可近似写为有限级数形式:

式中:Δκ为波数步长;Δω为频率步长。一般地,ΔκΔω需足够小以保证式(3a)可近似代替式(1)。

在式(3a)中,定义为如下形式:

根据式(2)的基本条件,可知XnmYnm应满足如下基本条件:

显然,XnmYnm是零均值的标准正交随机变量。

将式(4)代入式(3a)中,并考虑到因此,式(3a)进一步写为:

式中,为零均值的标准正交随机变量,即:

式(6)即为1V-2D连续时空随机场的基于标准正交随机变量的波数谱表达。

对于式(6),其截断的相对误差为:

式中,为截断波数,ωu=为截断频率。一般地,相对误差本文取ε≤ 0 .05。

进一步地,令式(6)中的标准正交随机变量为:

式中,为区间[0,2π)上均匀分布的相互独立随机相位角。显然,式(9)满足式(7)的基本条件。

于是,将式(9)代入式(6)中,即可得到1V-2D连续时空随机场的基于随机相位角模拟公式[18]

式中,f1(x,t)为基于随机相位角的模拟时空随机场。

从上述推导中,可以发现基于随机相位角的波数谱模拟公式是基于标准正交随机变量的波数谱的一种特例,即通过在基于标准正交随机变量的波数谱中定义式(9)即可得到基于随机相位角的波数谱模拟式(10)。显然,只要能够满足基本条件式(7),还可以将标准正交随机变量定义为不同的随机函数形式,从而得到不同类型的波数谱模拟公式。

事实上,式(9)可以理解为是对标准正交随机变量施加的一种约束,从而使式(6)中随机变量的数量由4×N×M减少为式(10)中的2×N×M。可见,通过施加随机约束可以有效地减少时空随机场模拟的随机度(随机变量的数量)。

2 标准正交随机变量的随机函数表达

直接利用基于随机相位角的波数谱式(10)来模拟时空随机场f0(x,t)时,所需随机变量的数量巨大,高达2×N×M个,这对于脉动风速时空随机场的模拟及抗风动力可靠度分析带来极大困难。为此,本文采用随机函数的降维思想[10,16],在式(6)中,将标准正交随机变量定义为概率分布已知的基本随机向量Θ的正交随机函数形式,即显然,定义的随机函数形式必须满足基本条件式(7)。

为此,本文给出如下一类较常用的随机函数形式:

式中,基本随机变量Θ1Θ2相互独立,且在区间(0,2π)上服从均匀分布。可以证明,式(11)满足式(7)的基本条件。

一般地,还需要将式(11)定义的标准正交随机变量集按某种确定性方式一一映射到目标标准正交随机变量集其中分别为N×M阶矩阵的第个元素,记为同样地,分别为N×M阶矩阵的第个元素,记为X(p)[r]与Y(p)[r]。这种确定性的一一映射关系可以方便地采用 Matlab库函数rand(‘state’,0)和randperm(·)来实现。为此,可定义rand(‘state’,0)和temp=randperm(N*M), 并 令显然,rr之间的一一映射关系为于是,式(6)中所需的标准正交随机变量即可被唯一地确定。

将上述一一映射关系所确定的标准正交随机变量(p=1,2)代入式(6)中,得到连续时空随机场的波数谱-随机函数模拟公式:

式中,f2(x,t)为基于随机函数的模拟时空随机场;分别表示基本随机变量1Θ2Θ的函数,即:

式中,之间的一一映射关系即为之间的映射关系。

式(12)在表达形式上与基于随机相位角的波数谱模拟式(10)相似,但二者存在本质的区别,基于随机相位角的波数谱模拟式(10)在本质上属于Monte Carlo模拟方法,模拟时,需要对2×N×M个随机相位角进行大量的Monte Carlo随机抽样以获取令人满意的模拟精度,不仅增加了随机模拟的计算量,也为结构动力响应计算和可靠度分析带来巨大的挑战。在式(12)中,随机变量是基本随机变量的函数,随机变量的数量由2×N×M降低为2,从而可以采用数论方法[10]对基本随机变量选取代表性点集,实现仅用较少的代表性时程即可获得令人满意的模拟精度。事实上,与式(10)类似,式(12)同样可以理解为基于正交随机变量的波数谱表达的一种特例,本文称为基于随机函数的波数谱降维表达。

3 FFT算法

式(12)与式(10)具有相同的表达形式,因而可采用 FFT算法加快模拟速率。考虑到于是,将式(12)改写为如下的形式:

式中:Re[·]表示实部;q2=mod(p2/M),其中 mod(·)表示取余。

在式(14)中,的表达式为:

其中:

显然,从式(15)中可以看出,分别为的双重傅里叶变换,即:

式中:C(1)C(2)分别为组成的矩阵;B(1)B(2)分别为组成的矩阵;IFFTκ(·)和IFFTω(·)分别表示对波数和频率进行快速傅里叶逆变换;FFTκ(·)表示对波数进行快速傅里叶变换。

对比式(15)和式(12)可以发现,引入 FFT算法后,式(12)中的频率求和项和波数求和项被FFT算法替代,从而极大地提高了随机模拟的计算效率。

4 波数谱与功率谱的关系

在利用波数谱模拟脉动风速时空随机场时,首先需要获取脉动风速时空随机场的波数谱密度函数。脉动风速时空随机场的波数谱密度函数可通过功率谱密度函数与空间相干函数的傅里叶积分得到[19]

式中:G0(κ,ω)为单边波数谱密度函数;G(ω)为单边自功率谱密度函数;γ(ξ,ω)为空间相干函数;ξ为空间两点之间的距离。

在本文中,空间相干函数γ(ξ,ω)采 用Davenport空间相干函数[20]

式中:λ是衰减因子;U(z)为高度z处的平均风速。

将式(19)代入式(18)中,可得:

其中:

式(21)的推导过程见附录A。

需要指出的是,当对波数谱密度函数G0(κ,ω)在波数上进行积分时,波数谱密度函数即可转化为功率谱密度函数。这说明利用式(18)构造的脉动风速时空随机场的波数谱G0(κ,ω)在本质上与所选取的自功率谱密度函数G(ω)具有一致性[12]。在经典的谱表示中,脉动风速随机场往往被离散为一维

多变量(1D-nV)的随机向量过程,其对随机场的模拟是通过功率谱密度矩阵的Cholesky分解来实现的,因此也称为脉动风速随机场表达的离散形式。然而,在式(12)中,脉动风速时空随机场被描述为一系列由基本随机变量调制的时空随机波的叠加,其关于时间t和空间坐标x是连续函数,因此可将式(12)视为脉动风速时空随机场表达的连续形式。

5 数值算例

为了说明本文方法的有效性,以水平向跨度L=1600 m的连续时空随机场进行分析,取左端为坐标原点,沿水平方向建立x坐标轴。模拟水平空间跨度L上 16385个离散点中的第 1026点即x1=1 00.10 m(位置1),第1076点即x2=104.98 m(位置 2),第 1176点即x3=1 14.75 m(位置 3)处的脉动风速时程,同时与Monte Carlo模拟方法进行对比以验证本文方法的有效性。

在本算例中,脉动风速自功率谱采用 Kaimal风速谱[21]

式中:z为距离地面的高度;u*为剪切速度;k0为Von Karman系数,可取k0=0.4;z0为地面粗糙度系数,可取z0=0 .03 m 。

根据式(20)和式(21)计算脉动风速时空随机场的波数谱如下:

本文算例的具体参数及取值如表1所示。

图1为采用本文方法生成的脉动风速代表性样本(位置1、位置2和位置3处的第60条),从图中可知,代表性样本的时程曲线具有脉动风的基本特征。图2为采用Monte Carlo方法生成的脉动风速样本(位置1、位置2和位置3处的第60条)。从图1和图2可知,本文方法与Monte Carlo方法所生成的时程曲线相类似。

表1 计算参数及取值
Table 1 Values of calculated parameters

图1 本文方法生成的脉动风速代表性样本
Fig.1 Representative samples of wind velocity fluctuation generated by the random function method

图2 Monte Carlo方法生成的脉动风速样本
Fig.2 Sample functions of wind velocity fluctuation generated by the Monte Carlo method

表2为采用本文方法和Monte Carlo方法分别生成x1x2x3三点处脉动风速样本集合的均值相对误差平均值、标准差相对误差平均值与自功率谱相对误差平均值,其中,均值、标准差、自功率谱的相对误差定义见文献[22]。表3为采用本文方法和Monte Carlo方法分别生成x1x2x3三点处脉动风速样本集合所需平均时间。从表2和表3可知,采用本文方法生成233条代表性样本的相对误差均小于Monte Carlo方法生成相同数量样本的相对误差,虽然采用Monte Carlo方法生成1000条样本的相对误差略小于本文方法生成233条代表性样本的相对误差,但其需要的模拟时间是本文方法的2倍。这是因为本文方法在基于正交随机变量的波数谱表达基础上,引入了随机函数的约束,将原时空随机场中2×N×M=16,777,216个随机变量降低为2个基本随机变量,同时利用数论方法对基本随机变量选取代表性点集,从而仅需较少的代表性样本即可获得令人满意的模拟精度。此外,本文方法生成的脉动风速代表性样本构成一个完备的概率集,这为基于概率密度演化的结构随机动力响应与动力可靠度精细化分析提供基础。

表2 模拟精度的比较
Table 2 Comparison of simulation accuracy

表3 模拟效率的比较
Table 3 Comparison of simulation efficiency

图3和图4分别为本文方法生成的233条脉动风速代表性样本的均值、标准差与目标值的比较。从图3和图4可知,3个位置点处的233条代表性样本集合的均值、标准差与目标值均拟合一致,初步验证了本文方法的有效性。

图5为本文方法生成的233条脉动风速代表性样本及Monte Carlo方法生成的233条和1000条样本的自功率谱与目标谱比较,由于x1x2x3三点在同一高度处,因此本文仅给出x1点处自功率谱的比较图。从图5中可知,除低频部分少数点外,两种方法生成的脉动风速样本集合的自功率谱与目标谱均拟合一致,且本文方法生成的233条脉动风速代表性样本自功率谱与目标谱的拟合效果略好于Monte Carlo方法生成的233条样本,但略差于Monte Carlo方法生成的1000条样本,这与表2中自功率谱相对误差的计算结果相一致。

图3 均值的比较
Fig.3 Comparison between simulated mean values and target value

图4 标准差的比较
Fig.4 Comparison between simulated standard deviations and target value

图5 自功率谱的比较
Fig.5 Comparison between simulated auto-power spectrums and target value

图6为本文方法生成的233条脉动风速代表性样本的空间相干函数与目标Davenport相干函数比较,其中图6(a)、图6(b)和图6(c)分别为位置1与位置2、位置2与位置3及位置1与位置3处的模拟相干函数与目标相干函数比较。从图中可知,模拟的相干函数符合两点间隔越大则衰减速度越快的特征,且模拟值与目标值均拟合一致。图5和图6的模拟结果进一步说明了本文方法的有效性。

图6 相干函数的比较
Fig.6 Comparison between simulated coherent functions and target values

6 结论

在基于标准正交随机变量的波数谱方法基础上,通过引入标准正交随机变量集的随机函数表达形式,建立了连续时空随机场模拟的波数谱-随机函数方法。结合Kaimal谱及Davenport空间相干函数模型,进行了脉动风速时空随机场模拟的数值分析。研究表明,本文方法具有如下特点:

(1)脉动风速时空随机场模拟的波数谱-随机函数方法是一种连续表达形式,可以模拟建筑物上任一点的脉动风速代表性样本,为复杂结构的脉动风速时空随机场提供了一种有效的模拟方法。

(2)脉动风速时空随机场的波数谱-随机函数方法是一类降维的模拟方法,所需基本随机变量的数量为 2,可以利用数论方法选取基本随机变量的代表性点集,实现仅用数百条代表性样本即可在概率密度层次上描述脉动风场的概率特征。与Monte Carlo方法相比,本文方法所需的随机变量最少,生成的代表性样本数量少,且模拟精度高。同时,通过引入FFT算法,可极大地提高模拟效率。

(3)本文方法在模拟脉动风速时空随机场时,不仅模拟所需的代表性样本数量少,而且每条代表性样本具有给定的赋得概率,因此生成的脉动风速代表性样本集合是一个完备的概率集,从而可结合概率密度演化理论,实现对复杂工程结构的风致动力响应及抗风可靠度的精细化分析。

附录A:

为方便起见,记则有:

其中:

于是,解得:

将式(A3)代入式(A1)中,可得:

参考文献:

[1]黄本才,汪丛军.结构抗风分析原理及应用[M].上海:同济大学出版社,2008.Huang Bencai,Wang Congjun.Principle and application of structural wind resistance analysis [M].Shanghai:Tongji University Press,2008.(in Chinese)

[2]Shinozuka M.Monte Carlo solution of structural dynamics [J].Computers and Structures,1972,2(5/6):855―874.

[3]Shinozuka M,Jan C M.Digital simulation of random processes and its applications [J].Journal of Sound and Vibration,1972,25(1): 111―128.

[4]苏延文,黄国庆,彭留留.与反应谱相容的多点完全非平稳地震动随机过程的快速模拟[J].工程力学,2015,32(8): 141―148.Su Yanwen,Huang Guoqing,Peng Liuliu.Simulation of stochastic processes of fully non-stationary and response-spectrum-compatible multivariate ground motions [J].Engineering Mechanics,2015,32(8): 141―148.(in Chinese)

[5]Paola M D,Gullo I.Digital generation of multivariate wind field processes [J].Probabilistic Engineering Mechanics,2001,16(1): 1―10.

[6]Paola M D.Digital simulation of wind field velocity [J].Journal of Wind Engineering and Industrial Aerodynamics,1998,74/75/76: 91―109.

[7]Yang J N.Simulation of random envelope processes [J].Journal of Sound and Vibration,1972,21(1): 73―85.

[8]Deodatis G.Simulation of ergodic multivariate stochastic processes [J].Journal of Engineering Mechanics,1996,122(8): 778―787.

[9]Chen L,Letchford C W.Simulation of multivariate stationary Gaussian stochastic processes: hybrid spectral representation and proper orthogonal decomposition approach [J].Journal of Engineering Mechanics,2005,131(8): 801―808.

[10]Liu Zhangjun,Liu Zenghui,Peng Yongbo.Simulation of multivariate stationary stochastic processes using dimension-reduction representation methods [J].Journal of Sound and Vibration,2018,418: 144―162.

[11]Deodatis G,Shinozuka M.Simulation of seismic ground motion using stochastic waves [J].Journal of Engineering Mechanics,1989,115(12): 2723―2737.

[12]Benowitz B A,Deodatis G.Simulation of wind velocities on long span structures: a novel stochastic wave based model [J].Journal of Wind Engineering and Industrial Aerodynamics,2015,147: 154―163.

[13]Carassale L,Solari G.Wind modes for structural dynamics: a continuous approach [J].Probabilistic Engineering Mechanics,2002,17(2): 157―166.

[14]刘章军,王磊,黄帅.非平稳随机地震作用的结构整体可靠度分析[J].工程力学,2015,32(12): 225―232.Liu Zhangjun,Wang Lei,Huang Shuai.Global reliability analysis of structures under non-stationary random earthquake excitations [J].Engineering Mechanics,2015,32(12): 225―232.(in Chinese)

[15]Liu Zhangjun,Liu Wei,Peng Yongbo.Random function based spectral representation of stationary and non-stationary stochastic processes [J].Probabilistic Engineering Mechanics,2016,45: 115―126.

[16]Liu Zhangjun,Liu Zixin,Peng Yongbo.Dimension reduction of karhunen-loeve expansion for simulation of stochastic processes [J].Journal of Sound and Vibration,2017,408: 168―189.

[17]Shinozuka M,Deodatis G.Simulation of multidimensional Gaussian stochastic fields by spectral representation [J].Applied Mechanics Review,1996,49(1): 29―53.

[18]Chen X,Kareem A.Proper orthogonal decompositionbased modeling,analysis,and simulation of dynamic wind load effects on structures [J].Journal of Engineering Mechanics,2005,131(4): 325―339.

[19]Zerva A.Seismic ground motion simulations from a class of spatial variability models [J].Earthquake Engineering and Structural Dynamics,1992,21(4): 351―361.

[20]Davenport A.The spectrum of horizontal gustiness near the ground in high winds [J].Quarterly Journal of the Royal Meteorological Society,1961,87(372): 194―211.

[21]Kaimal J C,Wyngaard J C,Izumi Y,et al.Spectral characteristics of surface-layer turbulence [J].Quarterly Journal of the Royal Meteorological Society,1972,98(417): 563―589.

[22]Liu Zhangjun,Liu Zixin,Chen Denghong.Probability density evolution of a nonlinear concrete gravity dam subjected to non-stationary seismic ground motion [J].Journal of Engineering Mechanics,ASCE,2018,144(1):04017157.

SIMULATION OF FLUCTUATING WIND VELOCITY CONTINUOUS STOCHASTIC FIELD BY DIMENSION REDUCTION APPROACH

LIU Zhang-jun1,2,YE Yong-you2,LIU Zeng-hui2
(1.School of Civil Engineering and Architecture,Wuhan Institute of Technology,Wuhan,Hubei 430074,China;2.College of Civil Engineering & Architecture,China Three Gorges University,Yichang,Hubei 443002,China)

Abstract:Based on the frequency-wavenumber spectrum representation correlating with the standard orthogonal random variables,a hybrid approach of frequency-wavenumber spectrum and a random function for simulating the continuous spatio-temporal stochastic field is proposed by introducing the random function of standard orthogonal random variable sets.Meanwhile,the simulation efficiency of the proposed approach is greatly enhanced by employing Fast Fourier Transform(FFT)algorithm technique.Benefiting from the proposed approach,the probability characteristics of the spatio-temporal stochastic field can be described on the probability density level with only two elementary random variables.Therefore,the complete representative point sets with assigned probabilities of the elementary random variables can be obtained through the number-theoretical method.As a result,the dimension reduction representation of the continuous spatio-temporal stochastic field can be realized.Numerical examples indicate that when using the same number of samples and taking the efficiency and accuracy into consideration at the same time,the proposed approach have a similar simulation result to the conventional frequency-wavenumber spectrum representation.However,the smallest number of the elementary random variables is needed in the proposed approach,which leads to a smaller number of representative samples with a complete probability set.Consequently,it could naturally be combined with the probability density evolutionary method(PDEM)to carry out the accurate analysis of stochastic dynamic response and dynamic reliability assessment of engineering structures.Finally,combining the Kaimal fluctuating wind velocity spectrum with Davenport spatial coherence function,a numerical example of simulation for horizontal-fluctuating-wind velocity continuous stochastic field is presented to verify the accuracy and superiority of the proposed approach.

Key words:fluctuating wind velocity; spatio-temporal stochastic field; frequency-wavenumber spectrum;random function; FFT technique; dimension reduction

刘增辉(1991―),男,湖北人,硕士生,主要从事工程结构抗风研究(E-mail: liuzenghuictgu@163.com).

叶永友(1993―),男,湖北人,硕士生,主要从事工程结构抗风研究(E-mail: yeyongyou1993@aliyun.com);

作者简介:

通讯作者:刘章军(1973―),男,湖北人,教授,博士,博导,主要从事工程结构抗震抗风研究(E-mail: liuzhangjun73@aliyun.com).

基金项目:国家自然科学基金项目(51778343,51278282,50808113)

收稿日期:2017-07-24;修改日期:2018-01-05

文章编号:1000-4750(2018)11-0008-09

doi:10.6052/j.issn.1000-4750.2017.07.0570

文献标志码:A

中图分类号:TU312.1