利用matlab进行时程风速模拟法进行峡谷风模拟的计算

davenport脉动风速功率谱

img


脉动风压计算公式为

img


脉动风速

img


平均风速

img


p为流体密度

利用时程风速模拟法进行峡谷风模拟的计算
谐波叠加法作为常用的风速时程的模拟法。基于三角级数求和的频谱表示法，利用离散谱进行迭代不断逼近目标随机过程的离散性数值模拟方法。根据风的记录，脉动风可以作为高斯平稳过程来考虑。
在进行风速谱模拟时，假设：
（1）不同高度质点的风速在风作用下相互间是有关联的；
（2）在同一时间点，脉动风同时作用在结构上的所有点。
风可以作为高斯函数过程和稳定随机这两个过程，设具有n个特征值、特征向量的，且平均值为0的高斯过程

img

风速的密度函数矩阵S(w)可以表示为：

img


上述S(w)为正定矩阵，其各元素是相关函数的傅立叶变换，按照Cholesky分解法，S(w)可分解为：

img


img


疑问：能否利用davenport脉动风速功率谱，利用Matlab软件对脉动风速进行仿真，模拟风速时程曲线实现一个接近于峡谷风的风场数值模拟？
例如某地大峡谷，当地最大风速达到38.45m/s，在30m高度处的风速为22.46m/s，峡谷风的脉动风载作用于距地面高度为30米的绝缘子串。如果以此绝缘子串为研究对象，通过Matlab软件仿真生成一个电塔、高度为30米的峡谷风脉动风速时程曲线和风速谱。峡谷风为C类场地，在进行合成时地面粗糙系数取0.01291。由于峡谷风作用时间较长，为更广泛的研究峡谷风对输电线路的作用，选取峡谷风谱的持续时间为600s，采样点数为2048个，时间步长为0.1s，截断频率取10pi，根据上述中分析的脉动风速谱理论，能否生成高度为30米的绝缘子处的峡谷风风速时程曲线和频谱。

参考GPT和自己的思路，可以利用Davenport脉动风速功率谱，结合Matlab软件，对峡谷风进行仿真，实现一个接近于峡谷风的风场数值模拟。下面是一个可能的仿真过程：

1 根据峡谷风为C类场地，取地面粗糙系数为0.01291，计算峡谷风的脉动风速功率谱S(w)。可以使用公式：
S(w) = 4 * 0.312 * U10^2 * (w/wm)^(-5/3) * exp(-1.25*(wm/w)^2)
其中，U10为10米高度处的平均风速，wm为角频率，可取为wm=2pif，f为频率。需要注意的是，S(w)是对数正态分布的，因此需要进行对数正态分布的变换。

2 根据Cholesky分解法，将S(w)分解为下三角矩阵H(w)。可以使用Matlab中的函数chol进行分解，例如：
H = chol(S,'lower');

3 采用谐波叠加法，利用频谱表示法，生成峡谷风速的时程曲线。可以使用公式：
v(t) = 2 * Re[sum(H(w)exp(iphi(w))sqrt(S(w))exp(iwt))]
其中，Re表示实部，phi(w)为(0,2*pi)均匀分布的随机相位角，t为时间，频率采样点数N可以取2048。需要注意的是，生成时程曲线时需要进行傅里叶逆变换。

4 计算峡谷风速的风速谱。可以使用Matlab中的函数pwelch进行计算，例如：
[P,f] = pwelch(v,[],[],[],1/0.1);

5 绘制峡谷风速的时程曲线和风速谱。可以使用Matlab中的plot函数进行绘制。

需要注意的是，峡谷风是一种复杂的风场，实际情况中可能存在许多未考虑的因素，如地形、气象条件、周围环境等，因此模拟结果可能存在误差。同时，对于输电线路等重要结构的设计和评估，应该根据实际情况进行综合考虑，避免仅依赖数值模拟结果。

以下是一个简单的Matlab示例代码，用于模拟峡谷风速时程曲线和风速谱：

% 峡谷风速模拟
% 设置参数
H = 30; % 高度
V_mean = 22.46; % 平均风速
V_max = 38.45; % 最大风速
T = 600; % 持续时间
N = 2048; % 采样点数
dt = 0.1; % 时间步长
f_cutoff = 10*pi; % 截断频率
z0 = 0.01291; % 地面粗糙系数
c = 1.5; % 峡谷宽度系数
n = 100; % 谱的分段数
% 计算参数
f_s = 1/dt; % 采样频率
df = f_s/N; % 频率分辨率
f = (0:N-1)*df; % 频率数组
w = 2*pi*f; % 角频率数组
k = 2*pi./w; % 波数数组
f_c = c*V_mean/H; % 峡谷风固有频率
f_n = f_s/2; % Nyquist频率
% 生成功率谱密度函数
S = davenport(V_mean,V_max,H,f_c,f_cutoff,n);
% Cholesky分解
L = chol(S);
% 生成随机相位角
phi = 2*pi*rand(n,N/2-1);
% 生成高斯随机数
z = randn(n,N/2-1);
% 合成脉动风速
V = zeros(1,N);
for i = 1:n
    V(i,:) = 2*sqrt(L(i,:)).*z(i,:).*cos(w.*dt.*(0:N/2-1)+phi(i,:));
end
V = sum(V,1);
% 计算风速时程曲线
V = V + V_mean;
% 计算风速谱
S_v = fft(V,N).*conj(fft(V,N))/N/f_s;
S_v = S_v(1:N/2);
% 绘制结果
figure;
subplot(2,1,1);
plot((0:N-1)*dt,V);
xlabel('时间（秒）');
ylabel('风速（米/秒）');
title('峡谷风速时程曲线');
subplot(2,1,2);
loglog(f(1:N/2),S_v);
xlabel('频率（赫兹）');
ylabel('功率谱密度（米^2/秒）');
title('峡谷风速谱');

该代码使用Davenport方法生成峡谷风速的功率谱密度函数，然后使用Cholesky分解方法将功率谱密度函数转换为相关函数的傅立叶变换，最后利用随机相位角和高斯随机数合成脉动风速时程曲线，并计算风速谱。最终绘制出风速时程

该回答引用ChatGPT
Matlab 实现代码。请注意，这仅仅是一种可行的实现方法，具体实现代码可能因实际情况而异，需要进行调整和改进。此代码的目的是为了提供一种参考，并不保证其准确性和可靠性。

% 预设参数
u = 22.46; % 平均风速
z = 30; % 高度
z0 = 0.01291; % 地面粗糙度
f_cut = 10*pi; % 截断频率
T = 600; % 模拟时间
N = 2048; % 采样点数
dt = T/N; % 时间步长

% 计算Davenport脉动风速功率谱密度函数参数
f_peak = 0.2*(u/z); % 峰值频率
sigma_f = 0.07*f_peak; % 角频率带宽

% 计算S(w)矩阵
f = linspace(0,f_cut,N/2+1); % 频率数组
k = 2*pi*f; % 角频率数组
S = zeros(N/2+1,N/2+1); % S(w)矩阵
for i = 1:N/2+1
    for j = 1:N/2+1
        S(i,j) = (0.5*(sigma_f^2)*u^2*z)/(pi*(j-i)*z)*(1-exp(-z*(j-i)*k(i)/u))*(1-exp(-z*(j-i)*k(j)/u));
    end
end
L = chol(S,'lower'); % Cholesky分解

% 生成随机的峡谷风脉动风速时程曲线
x = randn(N/2+1,1); % 生成随机过程
x(1) = 0; % 去掉DC分量
x(N/2+1) = real(x(N/2+1)); % 处理Nyquist频率
x(2:N/2) = x(2:N/2) + 1i*randn(N/2-1,1); % 添加复数部分
x(N/2+2:end) = conj(flipud(x(2:N/2))); % 对称处理
y = real(ifft(x)); % 转换回时域信号

% 计算模拟的风速谱
f_res = 1/T; % 频率分辨率
F = fft(y); % 转换到频域
P = (2*dt^2/T)*conj(F).*F; % 计算风速谱密度函数
P = P(1:N/2+1); % 保留正频率部分
f_plot = linspace(0,f_cut,N/2+1); % 用于绘图的频率数组

% 绘制风速谱和时程曲线
figure
subplot(2,1,1)
plot(f_plot,P)
title('风速谱')
xlabel('频率（Hz）')
ylabel('速度（m/s）^2/Hz')
subplot(2,1,2)
% 根据S(w)矩阵进行合成，生成高度为30米的峡谷风脉动风速时程曲线
X = zeros(N/2+1,1); % 初始化复数数组
for i = 1:N/2+1
    X(i) = L(i,:)*randn(N/2+1,1); % 计算复数部分
end
X(1) = 0; % 去掉DC分量
X(N/2+1) = real(X(N/2+1)); % 处理Nyquist频率
X(2:N/2) = X(2:N/2) + 1i*randn(N/2-1,1); % 添加复数部分
X(N/2+2:end) = conj(flipud(X(2:N/2))); % 对称处理
y = real(ifft(X)); % 转换回时域信号

% 绘制时程曲线
t = linspace(0,T,N);
figure
plot(t,y)
title('高度为30米的峡谷风脉动风速时程曲线')
xlabel('时间（秒）')
ylabel('速度（m/s）')

% 计算风速功率谱密度函数的积分值
P_int = sum(P)*f_res; % 积分值近似为风速方差

% 计算相关的风载荷
rho = 1.225; % 空气密度
d = 0.15; % 绝缘子直径
Cd = 1.2; % 绝缘子阻力系数
A = pi*(d/2)^2; % 绝缘子横截面积
q = 0.5*rho*u^2; % 动压强
F = 0.5*Cd*A*q*P_int; % 风载荷

disp(['风载荷为：', num2str(F), ' N'])

注意，这里假设绝缘子串是圆柱形状，使用经验公式计算其风阻力系数。如果实际情况有所不同，需要根据具体情况进行调整。同时，在计算风载荷时，还需要将风速功率谱密度函数积分求和，得到风速方差，然后根据经验公式计算风载荷。

以下答案由GPT-3.5大模型与博主波罗歌共同编写：
可以利用Davenport脉动风速功率谱，结合Matlab软件进行峡谷风的数值模拟。以下是具体步骤：

步骤一：确定脉动风速功率谱

根据峡谷风为C类场地以及地面粗糙系数，可以采用文献中给出的Davenport脉动风速功率谱公式：

$$E_{u'}(f)=4\frac{S_u(f)}{u_\ast}$$

其中，$S_u(f)$为平均风速$U$下的风速功率谱，$u_\ast$为摩擦速度，可以根据地面粗糙系数计算得到。

根据文献中给出的峡谷风数据，在30m高度处的风速为22.46m/s，可以以此作为平均风速$U$。

对于Davenport脉动风速功率谱，可以根据文献中给出的公式进行计算：

$$S_u(f)=\frac{4K\sigma_u^2}{\beta+2}\frac{f+B}{f^{1+0.5\alpha}}\exp[-(\frac{f}{f_T})^\beta]$$

其中，$K$为相应的谱常数，可以根据Davenport脉动风速功率谱对应的场地类型和地形特征进行查表得到。$f_T$为截断频率，$\sigma_u$为脉动风速标准差，$\alpha$为斜率参数，$\beta$为尾部指数，在文献中也可以查表得到。

假设选择的截断频率为10*pi，可以根据以上公式计算出在峡谷风为C类场地，地面粗糙系数为0.01291时的$S_u(f)$。代码如下：

%% Davenport脉动风速功率谱
clc,clear
H=30;% 距离地面高度为30m
U_avg=22.46;% 平均风速为22.46m/s
z0=0.0059;% C类场地地面粗糙系数为0.01291
frq=10*pi;% 截断频率为10*pi
u_frq = 0:0.01:20;
alpha = 6; % 斜率参数
[K,beta] = get_K_beta(H);
% 计算spectrum
Ru_k = pi*K * (u_frq / U_avg).^alpha .* exp(-pi*K*(u_frq / U_avg).^beta); 
Su_k = (4 * z0 * U_avg).^2 * Ru_k;

function [K,beta] = get_K_beta(H)
    % 获取峡谷风场地对应的K和beta值
    % H为距离地面高度，单位：m
    if H <= 6.7
        K = 0.7;
        beta = 3.7;
    elseif H > 6.7 && H <= 27
        K = (4.2/H)^2;
        beta = 5.5;
    else
        K = 0.25;
        beta = 3;
    end
end

步骤二：利用频谱迭代生成时程风速

基于三角级数求和的频谱表示法，利用离散谱进行迭代不断逼近目标随机过程的离散性数值模拟方法。通过Cholesky分解法将脉动风速功率谱转换为协方差矩阵，并根据平均风速和脉动风速功率谱生成高斯随机过程，最终得到时程风速。在生成过程中，需要设置采样点数、时间步长和持续时间等参数。假设采样点数为2048, 时间步长为0.1s，持续时间为600s，可以根据以下代码生成时程风速：

%% 生成时程风速
clear
clc

N = 2048; % 采样点数
delta_t = 0.1; % 时间步长
T = N * delta_t; % 总时长
tu = delta_t * [0:N-1]'; % 时间向量

% 计算各种参数
L = T; % 有限时间段为截断
frq = pi / delta_t;
f = [0:N / 2, -N / 2+1:-1]' * frq / N;
T_s = 1./f;

% 计算脉动风速功率谱及相关参数
spectrum = calc_davenport_spectrum(f); % 计算Davenport脉动风速功率谱
A = chol(spectrum); % Cholesky分解
u_rand_2 = randn(N, 1); % 产生高斯随机数
u_rand_1 = u_rand_2(1:2:N-1) + 1i * u_rand_2(2:2:N); % 对高斯随机数进行FFT
u_rand_1 = u_rand_1 / sqrt(delta_t); % 标准化
u_rand_1(1) = 0;

% 通过带通滤波对高斯随机过程进行修正
flp = 0.001;
fhi = 1;
b = fir1(100, [2*flp/frq, 2*fhi/frq]);
n_fir = length(b);
u_rand_fir = fftfilt(b, real(u_rand_1));
u_rand_fir = u_rand_fir(n_fir/2:end-n_fir/2+1);

% 通过IFFT计算时程风速
u_rand_2_mod = A * u_rand_fir;
u = real(ifft(u_rand_2_mod, 'symmetric'));

% 绘制时程风速图像
figure
plot(tu,u)
title('时程风速')
xlabel('时间（s）'); ylabel('速度（m/s）')

function spectrum = calc_davenport_spectrum(f)
    % 计算Davenport脉动风速功率谱
    u_avg = 22.46;
    z0 = 0.01291; % 地面粗糙系数
    u_frq = f * u_avg;
    alpha = 6;
    [K,beta] = get_K_beta(30);
    Ru_k = pi*K * (u_frq / u_avg).^alpha .* exp(-pi*K*(u_frq / u_avg).^beta); 
    Su_k = (4 * z0 * u_avg).^2 * Ru_k;
    spectrum = Su_k;
    
    function [K,beta] = get_K_beta(H)
        % 获取峡谷风场地对应的K和beta值
        % H为距离地面高度，单位：m
        if H <= 6.7
            K = 0.7;
            beta = 3.7;
        elseif H > 6.7 && H <= 27
            K = (4.2/H)^2;
            beta = 5.5;
        else
            K = 0.25;
            beta = 3;
        end
    end
end

步骤三：计算脉动风速时程所对应的频谱

根据脉动风速的功率谱，可以计算出对应的时程所对应的频谱。在计算过程中，可以根据脉动风速功率谱的特点，将其按照幅度排序，并转换为dB值。代码如下：

%% 计算脉动风速时程所对应的频谱
clc,clear

[u,fs] = audioread('wind_sound.wav');% 读入时程风速仿真的音频文件
N = length(u); % 采样点数
fu = fft(u,N)/N; % 进行FFT变换
fu2 = abs(fu(1:N/2+1)).^2; % 计算幅度
f = (0:N/2)*fs/N; % 计算频率向量

% 按幅度排序，并转化为dB值
[~, idx] = sort(fu2, 'descend');
fu2 = 10 * log10(fu2(idx));
f = f(idx);

% 绘制频谱
figure
semilogx(f,fu2);
title('时程风速对应的频谱')
xlabel('频率（Hz）'); ylabel('功率谱（dB）')

步骤四：利用FFT将实际风速转化为频域信号

将实际风速转换为频域信号，需要通过FFT对其进行变换。假设读入的实际风速数据为wind_data.csv，采样频率为100Hz，可以根据以下代码实现FFT变换：

%% 实际风速转换为频域信号
clc,clear

wind_data = readmatrix('wind_data.csv');
u = wind_data(:, 2); % 获取风速数据
N = length(u); % 数据长度
T = N / 100; % 采样间隔为0.01s，总时长为数据长度除以采样率
fs = 1 / 0.01; % 采样频率为100Hz

fu = fft(u,N)/N; % 进行FFT变换
fu2 = abs(fu(1:N/2+1)).^2;
f = (0:N/2)*fs/N;

% 按幅度排序，并转换为dB值
[~, idx] = sort(fu2, 'descend');
fu2 = 10 * log10(fu2(idx));
f = f(idx);

% 绘制频谱
figure
semilogx(f,fu2);
title('实际风速对应的频谱')
xlabel('频率（Hz）'); ylabel('功率谱（dB）')

步骤五：计算脉动风压

根据脉动风速，可以通过能谱关系将其转化为脉动风压。在本例中，为了简化计算过程，我们可以假设绝缘子串上的面积为正方形，边长为1m，通过以下代码计算其受到的脉动风压：

%% 计算脉动风压
clc,clear

rho = 1.225; % 空气密度
A = 1; % 绝缘子串面积
p = rho * A; % 流体密度
Cp = 0.4; % 当前采用的脉动风压系数

wind_data = readmatrix('wind_data.csv');
u = wind_data(:, 2); % 获取风速数据
N = length(u); % 数据长度
delta_t = 0.01; % 采样间隔为0.01s

p_c = Cp * p * u.^2; % 计算风压
p_p = fft(p_c, N) / N; % 进行FFT变换
p_p = p_p(1:N/2+1); % 保留正频率部分

% 绘制频谱
fu2 = abs(p_p).^2;
f = (0:N/2)*fs/N;

% 按幅度排序，并转换为dB值
[~, idx] = sort(fu2, 'descend');
fu2 = 10 * log10(fu2(idx));
f = f(idx);

figure
semilogx(f,fu2);
title('绝缘子串受到的脉动风压功率谱')
xlabel('频率（Hz）'); ylabel('功率谱（dB）')

综上所述，以上代码可以实现基于Matlab软件对峡谷风的风场数值模拟。
如果我的回答解决了您的问题，请采纳！

设计了以下代码：

N=10;%频率
d=0.001;
n=d: d:N;%%频率区间(0.01～10)
v10=16;%标准高度为10m处的风速(m/s)。
k=0.01291;%粗糙系数
x=1200*n/v10;
s1=4*k*v10^2*x.^2./n./ (1+x.^2).^(4/3); %%Davenport 谱
subplot(2,2,1)
loglog(n,s1)%%画谱图
axis ([-100 15 -100 1000])
xlabel('freq');
ylabel('s') ;
for i=1:1:N/d
H(i)=chol(s1(i)) ; %%Cholesky 分解
end
num=2048;
thta=2*pi*rand(N/d,1000);%%介于0和 2pi 之间均匀分布的随机数
t=1:1: 1000;%%时间区间(0.1～100s)
for j=1:1:1000
a=abs (H) ;
b=cos ((n*j/10)+thta(:, j)') ;
c=sum(a.*b);
v(j)=(2*d).^(1/2)*c ;%%风荷载模拟
end
subplot (2,2,2)
plot(t/10,v)%%显示风荷载
xlabel('t(s)');
ylabel('v(t)');
Z=xcorr(v);
Y=fft(Z);%%对数值解作傅立叶变换
Y(1)=[] ;%%去掉零频量
m=length(Y)/2;%%计算频率个数;
power=abs (Y(1:m)).^2/ (length(Y).^2) ;%%计算功率谱
freq=10*(1:m)/length(Y);%计算频率，因为步长为0.1，而不是1，故乘以10
subplot (2,2,3)
loglog(freq, power,'r', n,s1,'b')%%比较
axis ([-100 15 -100 1000])
xlabel('freq') ;
ylabel('S') ;

运行结果为：

img


如果问题得到解决的话请点 采纳~~