NSRFG湍流入口合成方法分享(我在验证结果的时候发现NSRFG方法和CDRFG方法均存在水平方向流场均一性问题，不知道对不对）。
%% NSRFG湍流合成方法
tic
clc,clear
%%生成0-fmax Hz范围的功率谱
fmax=25;df=0.001;fw=df:df:fmax;
N=length(fw);
fm=(2*(1:N)-1)/2*df;
Time_NSRFG=1000;%时间点数量
h=1/fmax/2;%时间步长
Fs=1/h;
%%入口尺寸和模拟点数量设置
x=1;%流场方向坐标
y=1;%水平方向坐标
z=[linspace(0.01,0.1,10),linspace(0.11,0.5,20),linspace(0.55,1,10)];%竖直方向坐标
N0_Y=length(y);N0_Z=length(z);N0_X=length(x);%X,Y,Z对应u,v,w，各方向上的点数
xyz=[repmat(x,length(z)*length(y),1),...
repmat(y',length(z),1),...
reshape(repmat(z,length(y),1),[],1)]';
%%切分Time_NSRFG成N0_T段，避免内存溢出
memory_capacity=0.4*10^9;%16G内存建议取值，内存越大计算速度越快
N0_t=min(fix(memory_capacity/(N0_Z*N0_Y)/N),Time_NSRFG);
N0_T=ceil(Time_NSRFG/N0_t);
%%相关性参数设置，
c_1=10;c_2=12;c_3=12;
Gammaix=1.8;Gammaiy=1.5;Gammaiz=1.8;
%%风场特性-风速
z_ref=0.4;%标准参考高度
alpha=0.25;%粗糙度指数
U_ref=10.0;%参考风速
I_10=0.23;
coefv=0.88;
coefw=0.55;
Uav=U_ref*(z/z_ref).^alpha;
Uav_map=(reshape(repmat(Uav,N0_Y,1),[],1))';
%%风场特性-湍流强度
Iu=0.1*(z/z_ref).^(-alpha-0.05);
Iv=coefv*Iu;
Iw=coefw*Iu;
I_map=[(reshape(repmat(Iu,N0_Y,1),[],1)),...
(reshape(repmat(Iv,N0_Y,1),[],1)),...
(reshape(repmat(Iw,N0_Y,1),[],1))]';=
%%风场特性-积分尺度
z_0=0.25;
Lu=linspace(0.25,1,N0_Z);
Lv=0.5*coefv^3*Lu;
Lw=0.5*coefw^3*Lu;
L_map=[(reshape(repmat(Lu,N0_Y,1),[],1)),...
(reshape(repmat(Lv,N0_Y,1),[],1)),...
(reshape(repmat(Lw,N0_Y,1),[],1))]';
%初始化参数
Su=zeros(N,N0_Z*N0_Y);Sv=zeros(N,N0_Z*N0_Y);Sw=zeros(N,N0_Z*N0_Y);
p_uvw=zeros(3,N);
p_u=zeros(N,N0_Z*N0_Y);p_v=zeros(N,N0_Z*N0_Y);p_w=zeros(N,N0_Z*N0_Y);
Lx=zeros(N,N0_Z*N0_Y);Ly=zeros(N,N0_Z*N0_Y);Lz=zeros(N,N0_Z*N0_Y);
q_1n=zeros(N,N0_Z*N0_Y);q_2n=zeros(N,N0_Z*N0_Y);q_3n=zeros(N,N0_Z*N0_Y);
A_n=zeros(N,N0_Z*N0_Y);B_n=zeros(N,N0_Z*N0_Y);
k_1n=zeros(N,N0_Z*N0_Y);k_2n=zeros(N,N0_Z*N0_Y);k_3n=zeros(N,N0_Z*N0_Y);
K_xyz=zeros(3,N0_Z*N0_Y,N);
x_1n=zeros(N,N0_Z*N0_Y);x_2n=zeros(N,N0_Z*N0_Y);x_3n=zeros(N,N0_Z*N0_Y);
X_xyz=zeros(3,N0_Z*N0_Y,N);
U_u0=zeros(N0_t,N0_Z*N0_Y,N);U_v0=zeros(N0_t,N0_Z*N0_Y,N);U_w0=zeros(N0_t,N0_Z*N0_Y,N);
U_u=zeros(1,N0_Z*N0_Y);U_v=zeros(1,N0_Z*N0_Y);U_w=zeros(1,N0_Z*N0_Y);
%球面随机采样
S = randn(N,3);
sign_n =bsxfun(@rdivide, S, sqrt(sum(S.*S, 2)));
plot3(sign_n(:,1),sign_n(:,2),sign_n(:,3),'.k', 'MarkerSize',10) %绘制原始散点数据
%%生成随机数\thea_n
thea_n=2*pi*rand(N+100,N0_Z*N0_Y);
fm=(2*(1:N)-1)/2*df;
%%生成参数矩阵
parfor n=1:N%频率范围
disp('n');disp(n);
Su(n,:)=4*(I_map(1,:).*Uav_map).^2.*(L_map(1,:)./Uav_map)./...
(1+70.8*(fm(n)*L_map(1,:)./Uav_map).^2).^(5/6);%Karman
Sv(n,:)=4*(I_map(2,:).*Uav_map).^2.*(L_map(2,:)./Uav_map).*...
(1+755.2*(fm(n)*L_map(2,:)./Uav_map).^2)./((1+283.2*(fm(n).*L_map(2,:)./Uav_map).^2).^(11/6));%Karman
Sw(n,:)=4*(I_map(3,:).*Uav_map).^2.*(L_map(3,:)./Uav_map).*...
(1+755.2*(fm(n)*L_map(3,:)./Uav_map).^2)./((1+283.2*(fm(n).*L_map(3,:)./Uav_map).^2).^(11/6));%Karman
p_u(n,:)=sign(sign_n(n,1))*sqrt(2*Su(n,:)*df);
p_v(n,:)=sign(sign_n(n,2))*sqrt(2*Sv(n,:)*df);
p_w(n,:)=sign(sign_n(n,3))*sqrt(2*Sw(n,:)*df);
Lx(n,:)=U_ref./((fm(n)))./c_1./Gammaix;
Ly(n,:)=U_ref./((fm(n)))./c_2./Gammaiy;
Lz(n,:)=U_ref./((fm(n)))./c_3./Gammaiz;
q_1n(n,:)=p_u(n,:)./Lx(n,:);
q_2n(n,:)=p_v(n,:)./Ly(n,:);
q_3n(n,:)=p_w(n,:)./Lz(n,:);
A_n(n,:)=sqrt((q_2n(n,:).^2+q_3n(n,:).^2).^2+(q_1n(n,:).*q_2n(n,:)).^2+(q_1n(n,:).*q_3n(n,:)).^2);
B_n(n,:)=sqrt(q_2n(n,:).^2+q_3n(n,:).^2);
k_1n(n,:)=-(q_2n(n,:).^2+q_3n(n,:).^2)./A_n(n,:).*sin(thea_n(n,1));
k_2n(n,:)=(q_1n(n,:).*q_2n(n,:))./A_n(n,:).*sin(thea_n(n,1))+q_3n(n,:)./B_n(n,:).*cos(thea_n(n,1));
k_3n(n,:)=(q_1n(n,:).*q_3n(n,:))./A_n(n,:).*sin(thea_n(n,1))-q_2n(n,:)./B_n(n,:).*cos(thea_n(n,1));
div1(n,:)=q_1n(n,:).* k_1n(n,:)+q_2n(n,:).* k_2n(n,:)+q_3n(n,:).* k_3n(n,:);
div2(n,:)= k_1n(n,:).* k_1n(n,:)+ k_2n(n,:).* k_2n(n,:)+ k_3n(n,:).* k_3n(n,:);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
K_xyz(:,:,n)=[k_1n(n,:);k_2n(n,:);k_3n(n,:)];%3x9
x_1n(n,:)=(xyz(1,:))./Lx(n,:);
x_2n(n,:)=(xyz(2,:))./Ly(n,:);
x_3n(n,:)=(xyz(3,:))./Lz(n,:);
X_xyz(:,:,n)=[x_1n(n,:);x_2n(n,:);x_3n(n,:)];%3x9
kjxj(n,:)=(diag(K_xyz(:,:,n)'*X_xyz(:,:,n)))';
end;disp('完成构建矩阵');
%%合成风速
leq_n=2*pi*rand(3,N);
t=(1:1:N0_t)*h;
ru=sign_n(:,1);rv=sign_n(:,2);rw=sign_n(:,3);
for it=1:N0_T
for n=1:N
disp(['N0_T--','it--','n','运行时间： ',num2str([N0_T,it,n,toc])]);
U_u0(:,:,n)=p_u(n,:).*sin(repmat(kjxj(n,:),N0_t,1)+2*pi*fm(n).*t'+repmat(leq_n(1,n),N0_t,1));
U_v0(:,:,n)=p_v(n,:).*sin(repmat(kjxj(n,:),N0_t,1)+2*pi*fm(n).*t'+repmat(leq_n(1,n),N0_t,1));
U_w0(:,:,n)=p_w(n,:).*sin(repmat(kjxj(n,:),N0_t,1)+2*pi*fm(n).*t'+repmat(leq_n(1,n),N0_t,1));
end;t=t+N0_t*h;
U_u=cat(1,U_u,sum(U_u0,3));
U_v=cat(1,U_v,sum(U_v0,3));
U_w=cat(1,U_w,sum(U_w0,3));
end
U=roundn(U_u(2:end,:)+Uav_map,-4)';
V=roundn(U_v(2:end,:),-4)';
W=roundn(U_w(2:end,:),-4)';
toc
disp(['运行时间: ',num2str(toc)]);
disp('=====accomplish======')
