汪江，陈喆，朱虹

(1.安徽省电力科学研究院，合肥市，230601;2.东南大学土木工程学院，南京市，210096)

0 引言

当前输电塔线设计时风荷载一般按照静力风考虑，忽略脉动风的影响［1］。实际上塔线体系中各部分动力耦合作用明显［2］，忽略脉动风作用可能会导致计算结果存在较大误差。但如采用现场实测和风洞试验则所需周期长、花费高，因此利用计算机程序针对输电塔结构进行脉动风速时程的精确模拟具有重要意义。

1 风速时程模拟的AR法

AR法通过线性自回归过滤器，将人工产生的均值为零的白噪声随机数，输出为具有特定频谱特性的随机数序列［3］。

M个点空间相关脉动风速时程V(X，Y，Z，t)的AR模型可表示为

式中:Ψk为AR模型的自回归系数矩阵;N(t)为M维独立正态分布随机过程向量。为简便起见，将V(X，Y，Z，t)简写为 V(t)，则根据式(1)可得

对式(2)做数学期望计算，并结合自相关函数的定义和性质，可得

式中:rV、rN分别为V(t)、N(t)的协方差矩阵。矩阵形式为

式(5)即为 AR模型的正则方程，式中，Ψ={I，Ψ1，…，Ψp}T，为(p+1)M ×M 阶矩阵，I是 M 阶单位矩阵;Op为pM×M阶零矩阵;R为(p+1)M×(p+1)M阶自相关Toeplitz矩阵

式中:rV(mΔt)是M×M阶矩阵，即:

式中:m=1，…，p。

根据随机振动理论，功率谱密度函数与相关函数满足Wiener-Khintchine公式，即:

式中:i，j=1，…，M。当 i≠j时，Sij(n)为空间第 i，j点的互功率谱密度函数;当i=j时，Sij(n)为第i点自功率谱密度函数。

通过式(7)、(8)求得式(6)后，代入式(5)，即可求解得到AR模型的系数矩阵Ψk和N(t)的协方差矩阵rN。

求得N(t)的协方差矩阵 rN后，N(t)按下式确定:

求出AR模型的系数矩阵Ψk和独立随机过程向量N(t)后，并假定初始时刻之前的风速为0，即t≤0时，V(t)为零向量，则由式(1)可得

式中:jΔt=0，…，T;k=1，…，p。计算可以得到 M 个空间点的具有时间空间相关、时间步长为Δt的离散脉动风速时程向量。

2 混合编程及模拟程序实现

2.1 混合编程方法

Matlab语言执行时程序逐行解释，当包含较多循环语句时存在执行效率低的缺点［4］。而 Visual C++拥有强大的执行效率，可与Matlab互补以便充分发挥各自的优势，适合大规模数值运算。通过Matlab接口程序(application programming interface，API)，混合编程有多种实现方式［5］。

本文采用Visual C++编写具有一定格式的程序(*.c函数的形式)，经过编译生成的Matlab动态链接子程序(*.dll函数的形式)，即MEX文件，然后以子函数或子过程的形式提供给Matlab应用程序，可以实现直接调用［6］。

2.2 程序实现

Matlab强大的矩阵运算功能及丰富的库函数，可简化风速模拟程序的编制，节省代码并减小出错，提高程序开发效率。如伪随机数的生成，通过调用rand函数，理论上其样本空间为21492个［7］。其他如矩阵相乘、除，矩阵的Cholesky分解，快速傅里叶变换(fast fourier transform，FFT)运算，风速时程的功率谱估计及绘图等均可以通过Matlab相应的运算格式和调用相应的库函数来实现，库函数可根据问题自动采用最优算法，保证了计算精度和模拟求解的稳定性。

在编制采用AR法模拟风速的程序时，为了保证数值积分运算求解精度，须经过大量循环迭代运算。Romberg定积分算法可提高计算效率和求解精度。为测试编制好的 Visual C++语言 Mex源文件rombergwk.c的迭代计算是否正确，下面采用该程序计算积分

＞＞ rombergwk(1，1，1，2，10－4) ［注:前2 个入口参数对此积分不需要，故任意取为1。其中截断误差(计算精度)为 10－4］

与ln2=0.693147相比，该结果的精度较高，而且还可以通过最后1个入口参数(10－4)来控制精度，证明了该程序迭代计算的正确性。通过结合各类风速谱函数编制程序时应写成Mex文件，在Matlab中编译连接后方可实现其对Visual C++语言程序的直接调用。

风速模拟程序运行环境为Windows XP SP3系统，AMD Turion(tm)64 X2 TL －581.90 GHz处理器，1.00 GB内存，Matlab 版本为7.0。采用混编技术计算用时1267.5 s，未采用混编技术计算用时2520.7 s，几乎是采用混编技术用时的2倍。由此可看出混合编程后，运算速度与执行效率大为增加。

3 算例分析

自立式输电塔属于空间杆件体系，可采用空间杆单元和梁单元建立输电塔的有限元模型［8］。鉴于输电塔自重较轻，自重二阶效应不明显，数值分析时可将输电塔视为线性结构进行计算［9］。图1为2E2-SZ2型220 kV输电塔三维模型。

图1 2E2-SZ2型220 kV输电塔三维模型Fig.1 220 kV 2E2-SZ2 transmission tower model

塔架总高 58.7 m，节点总数 624，单元总数1602，其中迎风面节点总数为161，迎风面单元总数为325。

利用自编AR脉动风模拟程序，得到了考虑空间相关的多点水平脉动风模拟结果，模拟过程中，主要参数由表1给出。

表1 风模拟程序主要参数Tab.1 Main parameters of simulation program

采用Welch法［10］得到功率谱。部分节点的水平风速模拟时程曲线及功率谱如图2所示。

对输电塔结构节点沿高度绘制功率谱后发现，高度位于35 m以下的节点，低频部分与高频部分模拟谱与目标谱之间逼近较好，个别节点中部分模拟谱与目标谱误差略大;高度位于35 m以上部分高低频部分模拟谱均与目标谱有一定误差，且随着高度增加有扩大的趋势。事实上，若采用Bartlett［11］法绘制功率谱，低频部分为平直型折线，高频部分为近似以目标谱为分界，呈上下震荡型折线，低频部分误差更为明显。因此综合考虑为减小瞬态效应，计算时可考虑忽略前面部分时程点数据。程序运行期间快速稳定，说明Visual C++与Matlab混合编程切实有效，图2中各点模拟功率谱与目标谱总体吻合良好，证明所编制风速时程程序的正确性。需要注意的是，模拟时间步长不宜小于0.1 s，时间步长过小时模拟精度大为降低，易产生错误。

图2 节点水平风速模拟时程曲线及功率谱Fig.2 The simulation time-history curves of horizontal wind speeds and power spectrum

4 结语

(1)三维空间风速时程的精确模拟是进行输电塔结构非线性风振响应分析的基础。本文针对输电塔的三维空间分析模型，采用混合编程机制编制的风速时程模拟程序具有良好的计算效率和稳定性。

(2)不同高度处节点的脉动风速时程趋势基本接近，但随着高度的增加，脉动风速波动范围略有减小，表明输电塔结构体系的脉动特性随着高度的增加略有降低。

(3)由于Davenport谱紊流尺度沿高不变，对于一般构筑物，如常见输电塔、预应力网壳结构、桥梁等结构均可用本风速时程模拟程序模拟风速，因而本程序有一定的通用性。

［1］DL/T 5092—1999P 110～500 kV架空送电线路设计技术规程［S］.北京:中国电力出版社，1999.

［2］ MomomuraY，MarukawaH，OkamuraT，etal. Full-scale measurements of wind-induced vibration of a transmission line system in a mountainous area［J］.Journal of Wind Engineering and Industrial Aerodynamics，1997，72(1):241-252.

［3］周云.结构风振控制的设计方法与应用［M］.北京:科学出版社，2009.

［4］刘保柱，苏彦华，张宏林.Matlab7.0从入门到精通［M］.北京:人民邮电出版社，2010.

［5］王素立.Matlab混合编程与工程应用［M］.北京:清华大学出版社，2008.

［6］董振海.Matlab编译程序和外部接口［M］.北京:国防工业出版社，2010.

［7］舒新琳，周岱.风速时程AR模型及其快速实现［J］.空间结构，2003，9(4):27-32.

［8］徐建国，石华军，郭勇，等.大跨越输电塔线体系有限元建模［J］.电力建设，2011，32(3):15-19.

［9］王肇民，Peil U.塔桅结构［M］.上海:同济大学出版社，1989:196-253.

［10］李志敏.预应力网壳结构的风振响应研究［D］.南京:东南大学，2007.

［11］邹鲲，袁俊泉，龚享铱.Matlab6.x信号处理［M］.北京:清华大学出版社，2002.