CFD中文网

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======')

@霜染丹枫 看前面的回复得知你要计算的是槽道流，能量谱应该不需要满足卡门谱，高斯谱可能更合适。TinF这个工具能完美的生成满足无散度条件和高斯谱的湍流（目前这个工具包还在开发中，还不能生成卡门谱的湍流），我有万嘉伟博士给我们讲课的视频，你可以看看。链接：https://pan.baidu.com/s/1iJy7CW6S3iwcDIvzZiAIZg
提取码：70xw

@李东岳 解决了，要感谢@xjwang大佬的帮助。

@xjwang 在 LES定义入口速度的问题（DSRFG方法） 中说：

@sunss
文件格式应该是：
points:
.......................................................
(
(1 1 0)
(0 1 0)
(1 0 0)
(0 0 0)
)
...........................................................
U:
..........................................................
(
(三维速度向量)
(三维速度向量)
(三维速度向量)
(三维速度向量)
)
.....................................................
所以你的设置有两个问题，一个是标出了4或者8这种数量，这个需要删掉，一个是需要把速度分别存储在对应的文件夹下边，不能在一个文件中。还有那些openfoam的星号可以删掉哈哈

多谢！已经成功导入风速数据。

@霜染丹枫 在 LES定义入口速度的问题（DSRFG方法） 中说：

@李东岳 如李老师所言：“像是这种速度波动越往后越少的，类似开始有湍流后来无湍流”。我有个疑惑，采用中心差分格式的扩散项可以将扰动向四周传递，迎风格式的对流项具有迁移性，即可以将扰动向下游传递。从上面来看，入口的扰动应该是可以持续向下游传递的。如果几何轴向长度非常长的话，在很远的下游是不是也会保持一定湍流状态？还是由于壁面摩擦力及流体粘性对扰动有削减作用，在下游一直衰减甚至脉动全无，就是这种入口给定方法本身就会出现的现象呢？

湍流减少的现象应该是湍流生成方法不满足无散度条件，而像满足无散度条件ATSM-L这种方法下游湍流强度不仅不会减少，还会稍微增强。推荐你们去了解一下万嘉伟开发的一款湍流入口生成工具。https://nheri-simcenter.github.io/TinF-Documentation/examples/examples.html

@xjwang 在 LES定义入口速度的问题（DSRFG方法） 中说：

@sunss
这个设置不是很难 你可以看一下这个视频 以及Eddylicious这个软件的document,他们里边也是用的这个边界条件
https://www.bilibili.com/video/BV1CV411C7Jr

你说的这个视频和工具包，我之前都看了，没能找到解决方法。我简化一下问题，假设inlet面上只有四个点，分别是（1 1 0）（1 2 0）（2 1 0）（2 2 0），我要导入两个时间步的速度信息，也就是总共有八个速度的信息，那timeVaryingMappedFixedValue边界条件的point文件和U文件写成下面的形式行不行。

你好同学，我最近再用NSRFG方法生成大气边界层风速入口，但不知道怎么将生成的风速数据导入openFoam中，看到你的这篇帖子后，尝试使用timeVaryingMappedFixedValue这个边界条件，但网上关于这个边界条件的用法几乎没有，而openFoam里有一个timeVaryingMappedFixedValue算例，但这个算例里U文件只有一步时间的数据，多步时间的格式不知道怎么写，请问你能告诉我U文件和point文件的格式吗，谢谢！我可以把我写的NSRFG方法分享出来，相互学习。