《基于线性滤波法的脉动风速模拟及其MATLAB程序的实现》由会员分享，可在线阅读，更多相关《基于线性滤波法的脉动风速模拟及其MATLAB程序的实现（7页珍藏版）》请在金锄头文库上搜索。
1、第 23卷第 4期2007年 8月结 构 工 程 师Structural EngineersVo. l 23, No . 4Aug . 2007收稿日期: 2006- 12- 16基于线性滤波法的脉动风速模拟及其 MATLAB程序的实现袁 波 应惠清 徐佳炜(同济大学建筑工程系, 上海 200092)摘 要 详细介绍了自回归技术公式, 同时采用此技术模拟具有空间相关性的随机风荷载, 给出了用自回归技术模拟脉动风荷载的步骤以及具体 MATLAB程序, 并把所编程序用于江阴长江大桥风荷载模拟, 结果表明效果较好。关键词 自回归, 空间相关, 脉动风速, MATLAB程序Si mulation of
2、 TurbulentW ind Velocity Based on LinearFilterM ethod andMATLAB Program RealizationYUAN Bo YI NG Huiqing XU Jiawei(Department ofBuilding Engineering , TongjiUniversity , Shanghai200092, China)Abstract In this paper ,the technology of autogressivemodel is introduced . Stochasticw ind loadingw ith ti
3、meand spatial correlation is generated by the autogressive -type si mulationmethod ,giving autogressive -type si mu -lation process andMATLAB progra m.Turbulent w ind velocity of Jiangying Changjiang Bridge over YangtseRiver is generated w ith this program and the result indicates that the programme
4、 is good .Keywords autogressivemode,l spatial correlation , turbulentw ind , MATLAB program1 引 言目前的抗风计算主要分为在频域和时域计算。频域分析是风振分析的一种常用方法, 但在频域内只能对结构进行线性分析, 要假定瞬时风压与风力之间的关系是线性的, 结构的特性也假定为线性的, 等等。强风作用下, 细长的大跨结构在风荷载作用下往往会产生几何或材料非线性,所以频域分析有时并不能反应结构的真实特性。要进行较为精确的分析, 只能借助于时域分析方法。时域方法可以直接了解结构的特性, 求出力和位移随时间的变化规
5、律和最大值。尽管时域方法计算复杂且耗时, 但随着计算机的发展, 该问题已经能很好地解决, 况且时域方法对于实验验证和科学研究很有必要。目前, 国内外对风速时程的模拟方法主要是CAW S ( Constant Ampfitude W ave Superposition)法、 WAWA(W avesw ithW eighted Amplitude)法及线性回归滤波器法。研究表明 1, 3, 对大型工程结构而言, 其自由度是非常大的。特别是大型空间结构, 对风荷载的三维分布都比较敏感, 必须精确模拟各点的风谱。 CAWS法与 WAWA 法计算量巨大, 所产生的风速过程不能考虑时间相关性;Soarl提
6、出的线性回归滤波器法很容易求出模型参数, 但模型精度受风谱变化的影响, 风谱的差异越大, 精度越低; I watani提出的线性回归滤波器法具有较好的普适性, 但模型参数一般采用迭代、递推的方法求解, 容易产生并累积误差, 导致模型的精度不够。2 风的基本特性风对结构的作用可以看成由平均风作用和脉动风作用两部分组成, 其中平均风在一定的时间间隔内, 风的大小和方向不随时间变化; 经过实测风时程记录可知, 平均风剖面沿结构高度往往按指数或者对数规律变化; 而脉动风荷载是随机荷载, 是风力中的动力部分, 它使结构产生随机振动。国内外的研究表明 4 , 5, 一般可将脉动风近似看作高斯过程来考虑。脉
7、动风速可以用下式表示:v( t) = C u( t) ( 1)式中, u( t) 为互不相关的高斯过程; C为互相关矩阵, 可以由后面的公式求得。因而, 在某一时刻 t 、 高度 z处, 风速 V(z, t)可以看作高度 z处的平均风速 v)( z)与脉动风速 u(z, t)之和。V(z, t) = v)(z) + C u(z, t) ( 2)2 . 1 脉动风的自功率谱历来的观测结果表明, 脉动风速谱近似服从高斯分布。脉动风速谱的统计通常有两种方法:一种是将强风记录通过超低频滤波器, 直接测出脉动风速的功率谱曲线; 另一种是把强风记录经过相关性分析, 获得风的相关曲线, 然后通过傅氏变换求得
8、功率谱曲线 4。其中, 最著名、 最常用的就是 Davenport谱。Davenport根据近百次观测结果提出了脉动风速谱的经验公式 4:Sv(n) = 4 KV210x2n( 1+ x2)4/3( 3)式中 x= 1 200 n /V10;K) ) 地面粗糙度系数;V10) 标准高度处 (通常为 10 m)的平均风速;n) ) 脉动风速的频率 (H z)。观测结果表明, 实际上随着高度的增加, 脉动风速谱的峰值将有所降低, 然而 Davenport谱并没有考虑这一变化。2 . 2 脉动风的互功率谱阵风的特性除了用自相关性描述外, 还可以用空间相关性来表示。强风观测表明, 各点风速、风向并不是
9、完全同步的, 甚至可能是完全无关的。阵风的空间相关性包括侧向左右相关和竖向相关, 在必要时还包括前后相关。试验表明, 这种相关性是随着空间两点间的距离增大而近似地按一种指数形式衰减的, 一般可写成Qz(z1, z2) = exp -| z1- z2|luQx(x1, x2) = exp-| x1- x2|lu( 4)式中, Qz(z1, z2)表示上下相关系数;Qx(x1, x2)表示侧向相关系数。若需要同时考虑结构竖向和侧向相关性时,空间任意两点的相关系数一般采用 Davenport提出的表达式:Q(n, x1, x2, z1, z2) =exp- nC2x(x1- x2)2+ C2z(z1
10、- z2)212(Vz1+ Vz2)( 5)式中 n) ) 频率, n= X /2P ;Cx, Cz) 分别表示侧向和竖向的衰减系数, 其值的变化范围较大, 根据Em il的建议, 常取为 Cz= 10 ,Cx=16 5;x1, x2, z1, z2) ) 分别为垂直于来流风的水平、 竖向坐标;Vz1, Vz2) 相应点的风速。由以上各式可以求出脉动风的相干函数, 那么脉动风速的交叉谱就可以通过相干函数求出:SViVj(x1, x2, y1, y2, n) = Qij(x1, x2, y1, y2, n)#SVi(x1, x2, n)SVj(y1, y2, n)( 6) 从整个式子可以看出,
11、谱密度函数完全表征了脉动风随机过程的频率含量和能量性质; 相干函数以及相位角反映了脉动风的相关性和相位关系, 因此这些参数常常被看成脉动风速模拟过程的特征参数。3 线性滤波法3 . 1 回归系数的求解过程线性滤波法又名白噪声滤波法 (White NoiseF iltrationMethod), 将随机过程抽象为满足一定条件的白噪声, 然后经过某一假定系统进行适当变换而拟合出该过程的时域过程。目前, 最常见的#56#StructuralEngineersVo.l 23 , No. 4Earthquake andW ind Resistance 是采用自回归模型 AM ( Autogressive
12、 Model) 6-15法模拟多维风速随机过程, M 个相关的随机过程 u( x, y, z, t) = u1( x, y, z, t), , uM(x, y, z,t) T可由下式求算:u(x, y, z, t) =6pk= 1 Wk # u(x, y, z, t - k$t) + N ( t) ( 7)式中, u(x, y, z, t- k$t) = u1(x, y, z, t- k$t), uM(x, y, z, t- k$t) T;N ( t) = N1( t), ,NM( t) T, 其中 Ni( t)均值为 0 , 具有给定方差的正态随机过程, i= 1 , 2 , ,M; Wk为
13、 M M 阶自回归矩阵; k= 1 , 2 , ,p。对空间任意一点 i( i= 1 , 2 , , M )具有时间差的随机过程 ui( t)和 ui( t- k$t) 的协方差可以表示为Riu(x, y, z, k$t) = E ui(x, y, z, t - k$t) - E ui(x, y, z, t - k$t) # ui(x, y, z, t) - E ui(x, y, z, t) ( 8) 由于 ui(x, y, z, t)和 ui(x, y, z, t- k$t)为均值为 0的平稳随机过程, 其协方差的值仅为时间差的函数, 上式可改写成为Riu x, y, z, k$t = E
14、ui(x, y, z, t - k$t)ui(x, y, z, t) ( 9) 式 ( 7)两边同时乘 ui(x, y, z, t- k$t) = u1(x, y, z, t- k$t), , uM(x, y, z, t- k$t) , 并且同时对方程求期望, 考虑到 N ( t)的均值为 0且与随机过程 ui(x, y, z, t)独立, 以及协方差 Ru( k$t)为偶函数, 可得到协方差 Ru(k$t)与回归系数 7k之间的关系有:RMP M= R)MP MP WMP M( 10)式中 R MP M= Ru($t), , Ru(p$t) T; WMP M= WT1, , WTpT。R)M
15、 P MP的计算式如下:R)MP MP=R11(0)R12($t),R1(p-1) (p - 2)$tR1p (p - 1)$tR21($t)R22(2$t),R1(p-1) (p - 1)$tR2p(0)sssssR(p- 1)1 (p - 2)$tR(p-1)2 (p - 1)$t,R(p- 1)(p- 1) (p - 4)$tR(p- 1)p (p - 3)$tRp1 (p - 1)$tRp 2( 0),Rp (p-1) (p - 3)$tRpp (p - 2)$t式中Rik( j$t) =R11(j$t),R1 M( j$t)swsRM 1(j$t),RMM( j$t); ( i, k
16、, j = 1 , 2 , , p) Wj =W11j,W1MjswsWM 1j,WMMj, (j = 1 , 2 , , p)这里 Ru( j$t)可根据维纳 -辛欣 (W iener -Khintch -ine)公式求得Riku(j$t) Q0Siku(n) cos( 2Pj$t) dn( i, k = 1, 2 , ,M )( 11) 通过以上各式可以求出回归系数矩阵 W。由式 ( 10)确定了自回归系数以后, 同时对式( 7) 两 边进 行 右乘 u ( t) T= u1( t ), ,uM( t) , 同时对方程求期望, 可求得 RN为RN= R (0) -6pk = 1WkRu(k
17、$t)( 12) 对 RN 进行乔利斯基 ( Cholesky) 分解,RN = L LT, 则N ( t) = L n( t) ( 13)式中L =L110,0L12L22,0LM 1LM 2,LMMLij=Rij-6i- 1k= 1LikLjkLjj, Lii=Rii-6i- 1k= 1L2iki, j = 1 , 2 , , M 求得系数矩阵 W和 RN后, 可以求出 M 个相关的随机风过程。将式 ( 7) 按时间间隔 $t离散化并考虑 t0时, ui( t) = 0 , 于是可以得到求随机过程的表达式:#57# 抗震与抗风#结构工程师第 23卷第 4期 u1(j$t)suM(j$t)=
18、6pk = 1 Wku1 ( j- k) $tsuM ( j- k)$t+ N1(j$t)sNM(j$t), ( j$t = 0, , T;k j)( 14) 通过式 ( 14)就可以求出 M 个时间、 空间相关、 时间间隔离散的脉动风速过程向量。在计算中一般假定初始时刻之前的风速为 0 , 即 t 0时,v( t) = 0 。3 . 2 脉动风速系数矩阵的确定应用线性自回归器可以产生脉动风速向量u( t), 余下的问题就是如何把这 N 个统计无关的随机过程 uj( t)转化为 N 个具有特定相关特性的随机过程脉动风速 v( t)。可以通过如下关系进行转化。vj( t) =6Ni= 1Ciju
19、i( t)( 15) 把式 ( 15) 两边同乘 Vk( t), 然后对方程两边求期望, 可得Rvjvk= E vj( t)vk( t) = E 6Ni= 1Cijui( t)6Ni= 1Cikui( t) ( 16) 如果随机过程 u( t)互相关矩阵为 D, 那么脉动风速的互相关矩阵可表示成 12:R = CDCT( 17)式中, 矩阵 C= Cij是一个下三角矩阵:C =C110000C12C22000,w00Cj1Cj2sCjj0Cn1Cn2sCnn- 1Cnn 由于随机过程向量 ui( t)中的元素是统计无关的, 且其方差均为同一常数, 因此, 矩阵 D 可表示为D = R2uI(
20、18)式中, I为 n n阶单位矩阵。把式 (18)代入式 ( 17)中, 可得R = R2uCICT( 19) 由于 u( t)为均值, 为标准的高斯过程, 故 Ru为单位矩阵。从而矩阵 C 中的元素可以表示为如下的递推公式:Cij=Rij-6i- 1k= 1cikcjkcjjCii= Rii-6i- 1k= 1c2ik i, j = 1 , 2 , , M( 20)4 线性滤波法的 MATLAB程序实现线性滤波法的 MATLAB 程序实现如图 1所示。5 算 例江阴长江大桥曾是我国跨度最大的悬索桥,其桥梁主跨度 1 385 m, 单侧吊杆 85根, 吊杆间距为 16m。运用根据前面所编制的
21、通用程序, 模拟了桥面上沿跨度方向均匀分布的 100个点的水平脉动风速。模拟时所用到的参数见表 1 。 表 1模拟计算的主要参数参数取值跨度 /m1 385主梁离地面高度 /m60地面 粗糙 /m0 . 003主梁 处平均 风 速 /(m# s- 1)40 . 0模拟点数10模拟 间距 /m13 . 85参数取值初 始频率 /H z0. 001截 至频率 /H z20频率增量0. 000 1计算的 阶数4初 始时间 /s0. 1时 间增量 /s0. 1参数取值总的时间点数2 000模拟采样点数6 000采样时间长度 / s2 048 在模拟过程中, 所采用的功率谱函数为 Si m iu谱:S(
22、 X) =200u2*fn( 1+ 50f )5/3( 21)式中 f ) ) 相似坐标, f = nz /V (z);z ) 离开地面的高度;V(z) ) 高度 z处的平均风速;u*) ) 摩擦风速, u*=KV(z) / ln(z /z0),与地面粗糙度相关, K U 0. 4 , z0为地面粗糙长度。 根据本文推导出的风荷载模拟公式以及基于Matlab工具箱编制了程序, 把此程序运用到了江阴长江大桥实际工程中风脉动速模拟, 不同点模#58#StructuralEngineersVo.l 23 , No. 4Earthquake andW ind Resistance 拟出来的风波通过 F
23、ourier变换所获得的谱与原风谱曲线比较 (图 2) 图 7)可以看出二者吻合得很好。图 1 线性滤波法的 MATLAB程序实现#59# 抗震与抗风#结构工程师第 23卷第 4期 图 2 第 1点模拟风速图 3 第 2点模拟风速图 4 第 20点模拟风速图 5 第 40点模拟风速图 6 第 60点模拟风速图 7 第 80点模拟风速6 结 论通过所编制的程序对江阴长江大桥脉动风速模拟结果表明, 线性滤波法是一种很好的模拟方法, 它所占用的计算机内存较少, 计算速度快。本文编制的程序可供相关工程研究人员参考。但是在进行模拟过程中, 需要注意随机数产生要足够大, 否则要影响模拟的效果。参考文献 1
24、 I ANNUZZIA, SPI NELLI P. A rtificiatW ind Genera -tion and StructuralResponse J.JournalofStructur -alEngineering -ASCE, 1987 , ( 12): 113. 2 王修琼, 张相庭. 混台回归模型及其在高层建筑风响应时域分析中的应用 J. 振动与冲击, 2000, 19( 1): 5- 7 . 3 张相庭. 结构风压和风振计算 M . 上海: 同济大学出版社, 1985. 4 李桂青, 曹 宏, 李秋胜, 等. 结构动力可靠性理论及其应用 M . 北京: 地震出版社, 199
25、3. 5 李元齐, 董石麟. 大跨空间结构风荷载模拟技术研究及程序编制 J. 空间结构, 2001, (9): 3- 11. 6 赵 臣, 张小刚, 吕伟平. 具有空间相关性风场的#60#StructuralEngineersVo.l 23 , No. 4Earthquake andW ind Resistance 计算机模拟 J. 空间结构, 1996, ( 5): 21- 25 . 7 边建烽, 魏德敏. 大跨空间结构风速时程的数值模拟理论 J. 暨南大学学报, 2005, ( 2): 87- 90 . 8 丁泉顺. 大跨桥梁耦合颤抖振响应的精细分析D. 上海: 同济大学博士论文, 200
26、1. 9 曹映泓. 大跨度桥梁非线性颤振和抖振时程分析D. 上海: 同济大学, 1999. 1. 10 DAVENPOR W H. Probability and Random ProceasesM . New York: M c Graw-H il. l 1990. 11 王 栎. 空间立体桁架结构体系的抗风特性研究D. 西安: 西安建筑科技大学硕士论文, 2004. 12 Chen L izhong , ChrisW. Letchford, Si mulation ofmu- ltivariate stationaryGaussian stochastic processes: hy -b
27、rid spectral representation and POD approach J.Journal of Engineering M echanics ,2005 , 131 ( 8):801- 808 . 13 L izhong Chen ,ChrisW Letchford . A deter m inistic -stochastic hybrid modelof downbursts and its i mpacton a cantilevered structure J.Engineering Struc -tures , 2004, 26( 5) 619- 629. 14 王之宏. 风荷载的模拟研究 J. 建筑结构学报,1994 , 15( 1): 2.#61# 抗震与抗风#结构工程师第 23卷第 4期
