%風電塔項目—風荷載時程模擬-諧波疊加法-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