LES定义入口速度的问题（DSRFG方法）

sunss

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

霜染丹枫

@李东岳还没有，还在找原因@@

霜染丹枫

@sunss 感谢分享。我做的是槽道流，能谱是通过槽道流实验得到的（沿轴向的一维谱），通过DSRFG方法生成了入口速度场，得到的入口速度场能谱分布与实验能谱拟合的较好，但是在计算中会出现我上面提到下游脉动减弱的问题。下图是我还原Huang文章中的卡门谱结果，整体看起来方法的实现应该是没有问题。

xjwang

@霜染丹枫
我做的风场模拟，所以我检察风洞模拟是否正确是通过检察计算域内的风速时程的功率谱是否和目标谱一致，也就是是否和在入口处生成的风速满足同样的能量分布。
不知道你们槽流最关心的是哪一个数值特征。比如，你可以验证turbulence intensity也可以验证power spectra，但不能通过一个瞬态的流场来判断模拟的好坏。一直不清楚你通过离散能量谱得到的脉动速度为什么要通过观察瞬时速度场来判断模拟的好坏，就算有很强的紊流特性也不能保证这个速度场的能量分布就和你目标谱一样。
上文中@sunss 同学说的那个tinf，你可以看一下，其中有dfsem, dfm等方法，特别是对于槽流，openfoam的tutorial中有一个RE395的算例（也许把名字记错了），应用的是DFSEM方法。
如@sunss 同学所说，也许高斯谱更适合你的功率谱，这样的话你可以直接应用fluent自带的那个对应RFG（Smirnov et al,2001)的生成方法，具体名字不记得了，你可以自己去查一下。
在其他回复中我提到过，DSRFG方法虽然理论上来说是无散度的，但如果应用在数值模拟中，就不是严格的无散度了，需要考虑网格尺寸，所以我一开始就问你网格尺寸和紊流积分尺度的数值。我们完成了一篇这样的论文，希望早点可以和你们分享。

xjwang

@sunss
感谢分享，前年看到他们那个工具包的时候还很简单，最近看了看已经变得很完整了，还有一些算例。你也是做风场模拟的吗？

李东岳

@霜染丹枫我觉得你这个问题很难解决，本身这种矩形网格槽道流，FVM本身精度很低会把能量耗散掉。

xzdbjx

各位老师专家好，我目前在做ABL的LES湍流入口时也碰到一些问题想和大家交流学习。我之前一直用Fluent来生成LES的湍流入口，近半年开始对Openfoam感兴趣，于是开始尝试把Fluent的相关算法移植到Openfoam当中，我采用的算法是华南理工杨易老师课题组提出的NSRFG算法，参考文献：
[1] Yu Y, Yang Y, Xie Z. A new inflow turbulence generator for large eddy simulation evaluation of wind effects on a standard high-rise building[J]. Building and Environment, 2018, 138: 300-313.
先把Fluent和of的计算结果图贴两张，我主要对PSD和湍流强度与平均风剖面进行了验证。
计算域
计算域.png

Fluent流场图及验证图

P1 是入口处的测点，P2是模型处的测点
Openfoam流场图与验证图

这里貌似of的结果还行，但实际上我认为Fluent的结果更合理，因为数值耗散等因素的影响，湍动能一定会在计算过程中耗散，因此模型处的PSD一定会更加吻合目标值。所以我目前怀疑我的case设置有问题，目前我采用的LES是WALE模型，0文件夹下应该要设置k 和 nut两个文件，我的设置分别如下：

/*--------------------------------*- C++ -*----------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  6
     \\/     M anipulation  |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    location    "0";
    object      k;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [0 2 -2 0 0 0 0];

internalField   uniform 0.0001;

boundaryField
{
    inlet
    {
        type            fixedValue;
        value           uniform 0;
    }
    ground
    {
		type	        kLowReWallFunction;		//For high and low Re
        value           $internalField;
    }
    top
    {
        type            symmetry;
    }
    side1
    {
        type            symmetry;
    }
    side2
    {
        type            symmetry;
    }
    outlet
    {
        type            zeroGradient;

    }
}


/*--------------------------------*- C++ -*----------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  6
     \\/     M anipulation  |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    location    "0";
    object      nut;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [0 2 -1 0 0 0 0];

internalField   uniform 1e-8;

boundaryField
{
    inlet
    {
        type            zeroGradient;
    }
    ground
    {
		type            nutUSpaldingWallFunction;
        value           uniform 1e-8;
    }
    top
    {
        type            symmetry;
    }
    side1
    {
        type            symmetry;
    }
    side2
    {
        type            symmetry;
    }
    outlet
    {
        type            zeroGradient;

    }
}


// ************************************************************************* //

这里我不太确定两个文件中internalField参数应该如何取值？另外，这两个参数对结果影响是否很大？请各位老师专家解答。谢谢！（P.S. of计算速度是真的快）

陈冠江

@xzdbjx
Hi，我也在用WALE算圆柱绕流的LES.

如果是LES的话你的kLowReWallFunction和kLowReWallFunction应该是要换成别的边界条件zeroGradient啥的，不需要用壁面函数。

internalField是一个初场，设的值不要太离谱能收敛就可以了，对最终结果影响不大。

另外WALE的输入参数似乎不需要k，你可以试试。或者就是我错了。

祝好.

xzdbjx

@陈冠江感谢您的回复，我修改下再看看结果。

xiaoyi0830

@coolhhh 对于第一个问题，我也是晕了，三维能谱，一维能谱，功率谱怎么在这儿就混用了。

xiaoyi0830

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

@coolhhh 对于第一个问题，我也是晕了，三维能谱，一维能谱，功率谱怎么在这儿就混用了。

另外第四个问题，就算满足了连续性条件，也不会满足动量方程，所以进入计算域的速度肯定和入口处是有差别的。

coolhhh

@xjwang
对第1个问题，其实从DSRFG到MDSRFG、CDRFG、NSRFG，在推导零散度条件过程中，都忽略了当生成inhomogeneous inflow时，系数k,p,q也是高度z的函数，将速度u的公式代入连续性方程，不仅仅需要满足kp=0，kq=0，还有类似于RFG方法推导零散度的额外项，是否近似等于0还需要要验证。

coolhhh

@xiaoyi0830
1、功率谱计算。时间序列$ x(t) $计算时间相关函数$ R_{xx}(\tau) $，再进行傅里叶变换得到功率谱$ S_{xx}(\omega) $。

(引用：Shin K, Hammond J. Fundamentals of signal processing for sound and vibration engineers[M]. England: John Wiley & Sons, 2008：245.)

2、风速功率谱。风速时程$ u(t) $对应时间序列$ x(t) $，求得的一维频谱（频域）就是对应此处的功率谱。一维频谱的自变量是频率或者角频率。

3、三维谱（波数域）。首先对三维空间相关函数进行三维的傅里叶变换得到的三维波谱，然后进行球面积分并乘以系数1/2，得到三维谱。因此要计算一个三维谱，需要整个三维空间的所有点的风速时程才能计算得到一个三维谱。

(引用：田内科斯 H., 兰姆利 J. L. 湍流初级教程[M]. 北京: 科学出版社, 2014：194-197.)

4、一维波谱（波数域）。一维波谱是由一维空间相关函数（即一个空间方向的相关函数），进行一维的傅里叶变换得到的。与一维频谱类似，只是一个是频域，一个是波数域。

5、三维von Karman谱（均匀流场）表达式如下。针对于非均匀流场的三维von Karman谱表达式，个人未找到相关文章阐述。

(引用：Durbin P A, Pettersson Reif B A. Statistical theory and modeling for turbulent flows[M].Wiley, 2011：29-261.)

6、三维von Karman谱与一维von Karman谱区别。

(引用：[1] Durbin P A, Pettersson Reif B A. Statistical theory and modeling for turbulent flows[M].Wiley, 2011：29-261.
[2] 陈铃伟. 基于傅里叶合成法的大气边界层脉动风场大涡模拟[D]. 黑龙江:哈尔滨工业大学,2017. DOI:10.7666/d.D01589559.)

7、从《Diffusion by a Random Velocity Field》(Kraichan, 1970)，到RFG方法，论文中都是未出现频谱的图，因为方法本身目标就是生成三维谱。DSRFG方法的理论推导也是实现三维谱，但是在应用方法时，用一维频谱表示三维谱，并最终的结果用频谱图表示，这与三维谱的原始定义是不同概念的，这两者不应该混淆。

8、CDRFG、NSRFG的理论推导就是为实现一维频谱。CDRFG方法既然目标实现频谱，可以简化一个循环，如下阐述

(引用：陈铃伟. 基于傅里叶合成法的大气边界层脉动风场大涡模拟[D]. 黑龙江:哈尔滨工业大学,2017. DOI:10.7666/d.D01589559.)

9、CDRFG方法的P、Q的原论文的公式有误，在原作者后来发的一篇文章中已对其进行了更正。

(引用：Melaku A, Bitsuamlak G, Elshaer A, et al. Synthetic inflow turbulence generation methods for LES study of tall building aerodynamics[J]. 2017.)

李东岳

@霜染丹枫流场衰变后来解决了么？

霜染丹枫

@李东岳李老师，问题算是解决了，正如楼上xjwang所说的那样，不应该通过监测瞬时速度云图的方式来判断，后来我通过数次调整积分尺度的值，并在下游检测不同位置点的能谱情况，发现是可以得到湍流场的，虽然能量有些衰减，但就目前来看，结果是可以接受的。

LiuSW

请问DSRFG或者相关的RFG方法，可以用于可压缩湍流大涡模拟吗？

sunss

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

coolhhh

@sunss 水平方向流场均一性问题，能否具体再说说吗，这个点不太理解具体表现是什么？

HITSC30

@sunss 看到我师兄做的DSRFG方法了，我也是做风场模拟的，可以加个联系方式探讨一下吗

HITSC30

@xzdbjx 朋友你好，我也是做这个方向的。可以加个联系方式聊一下吗

HITSC30

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

@sunss
应该是在constant/boundaryData/inlet文件夹下边有一个points的文件，然后每一个时间步都是一个文件夹，其中包含一个U文件。比如说，我在inlet文件夹下边在terminal中输入ls,会显示：

0 0.001 0.002 ...... points

而在每一个时间步文件夹下输入ls,比如在0.001文件夹下，会显示：

U

您好，我想请教一下这个数据的事情。在foam中运用并行运算，是在分块前面把inlet坐标输出然后运用算法求出inlet上每个点对应的风速，然后按照您的方法导入到constant/boundaryData/inlet中吗？还是在分块后，按照该方法进行呢？

AICFD

@hitsc30 decompose前

HITSC30

@aicfd 谢谢你的回复

BznW

@sunss 您好，感谢分享您的入口合成方法，我在使用您的程序生成的风速导入到空流场时似乎遇到了一些文章中提到的入口不满足divergence-free所导致的很大的虚拟压力波动（见后压力方差图），这会对放入建筑时建筑的表面脉动风载产生很大的影响，我看很多采用NSRFG的文章都没有提到这个事情，不知道您有没有遇到这种情况以及是如何解决的呢？

BznW

@xzdbjx 您好！我看您使用过NSRFG方法导入过OF中计算，请问您在瞬态计算中有遇到过入口附近有非常大的虚拟压力波动的问题吗？我想验证下出现这种情况是否是算法的问题。我看您对湍流强度平均风剖面以及功率谱进行了验证，我这里湍流强度、平均风剖面以及功率谱这些与风速相关的统计量都还正常，就是涉及到脉动压强就会出问题。

coolhhh

@BznW 压力脉动问题主要是因为入口质量通量不平衡导致的，NSRFG原论文中算出来的脉动压力系数也是非常大。简单的入口质量通量修正可以参考以下文章：
1.《Kim, Y., I.P. Castro and Z. Xie, Divergence-free turbulence inflow conditions for large-eddy simulations with incompressible flow solvers. Computers & Fluids, 2013. 84: p. 56-68.》
2.《Wang, Y. and X. Chen, Simulation of approaching boundary layer flow and wind loads on high-rise buildings by wall-modeled LES. Journal of Wind Engineering and Industrial Aerodynamics, 2020. 207: p. 104410.》

BznW

@coolhhh 好的好的，感谢回复！我去学习一下