谐波叠加法模拟风速时程

%风电塔项目—风荷载时程模拟-谐波叠加法-Davenport谱
clear;
clc;
%定义基本信息
nd=21; %模拟点数
Um=37.5; %基本风速，对于风电塔结构基本风速是其风机处的风速
H0=65; %基本高度，对于风电塔结构是风机处的高度
dt=0.1; %模拟风速时程的时间步长
wmax=2*pi*10; %上限角频率
T=6000; %模拟总点数
N=ceil（T/（2*nd））; %频率采样点数
M=2*N; %逆傅里叶变换长度，取采样点数的2倍
nT=M*nd; %一个模拟点的时间点数
t=dt*（1:1:nT）; %时间序列
dw=wmax/N; %角频率步长

%定义各类常数
r=0.22; %指数型风剖面，考虑地表粗糙度影响的无量纲幂指数，按中国规范取0.16（B）类，荷载规范中可找到
k=0.00464; %地面粗糙度系数k=0.00129（A）0.00215（B）0.00464（C）0.01291（D），此数据来源于文献
Cz=10; %竖直方向指数衰减系数

%定义模拟点竖向高度（不考虑横向相干性）
z（1）=5.880;z（2）=8.820;z（3）=11.760;z（4）=16.554;z（5）=19.494;
z（6）=22.434;z（7）=25.374;z（8）=28.314;z（9）=31.254;z（10）=34.194;
z（11）=37.109;z（12）=39.809;z（13）=42.514;z（14）=45.224;z（15）=47.939;
z（16）=50.659;z（17）=53.384;z（18）=56.114;z（19）=58.849;z（20）=61.569;
z（21）=63.209;        

%求平均风速
U0=Um*（10/H0）^r; %标准高度为10m的平均风速
for i=1:nd
    U（i）=U0*（z（i）/10）^r;
end %求指数型分布风剖面各模拟点风速

%定义风速谱（Davenport谱）
w=wmax/N:wmax/N:wmax;
n=w./（2*pi）;
x=1200*n./U0;
Sv=4*k*U0^2*x.^2./n./（1+x.^2）.^（4/3）; %Davenport谱

%谐波叠加法运算，利用FFT优化运算速度
V=zeros（ndnT）;
D=zeros（ndndN）;
for j=1:nd                            
    rand（‘state‘0）; %模拟随机数序列
    thet=2*pi*rand（jN）; %生成随机相位
    for i=1:N %生成频率序列
    w1（i）=（i-1）*dw+j/nd*dw;
    end
    n1=w1./（2*pi）; 
    x1=1200*n1./U0; %谱的参数
    Sv1=4*k*U0^2*x1.^2./n1./（1+x1.^2）.^（4/3）; %Davenport谱
    for j1=1:nd
        for j2=1:nd
            for i=1:N
                Coh（j1j2i）=（exp（（-n1（i）*sqrt（Cz^2*（abs（z（j1）-z（j2）））^2）*2）/（U（j1）+U（j2））））;
                S（j1j2i）=Sv1（i）*Coh（j1j2i）;
            end
        end
    end %生成谱矩阵
    for i=1:1:N
        H（::i）=chol（S（::i））;
        H（::i）=H（::i）‘;
    end %Ch