槽道流添加扰动速度场MATLAB代码

%% write by rqli 23/9/11
%% create initial velocity fields for channel flow
%% 生成无量纲形式的速度场u*;u* = u/u_{\tau}
%% 无量纲长度L*;L* = L/h;h = Ly/2;
%% input1： Retau
%% input2:  Lx,Ly,Lz
%% input3:  nx,ny,nz
%% input4:  1-->'laminar' or 2-->'time_average'
%% output:  MATLAB_Retauxxx_Initial.bin（无量纲速度场）
%% 需要调用函数creatMesh.m
clear;clc

%% 1、输入摩擦雷诺数
Retau = 180;base_flow = 2;
%% 2、设置计算域Lx->streamwise;Ly->wall-norm;Lz->spanwise
Lx = 4;Ly = 2;Lz = 2;h = Ly/2;
%% 3、设置网格参数
nx = 401;ny = 201;nz = 201;

%% 生成初始网格点
[X,Y,Z] = creatMesh(nx,ny,nz,Lx,Ly,Lz);
U = zeros(nx*ny*nz,1);
V = zeros(nx*ny*nz,1);
W = zeros(nx*ny*nz,1); 

%% 层流速度场参数设置
nu = 1.58e-5;
utau = Retau*nu/h;
%% 首先生成base_flow
if base_flow == 1% 对应'laminar'
    Ubar = (utau^2)*h/3/nu;
    for i = 1:length(U) 
        y = min(Y(i),2*h-Y(i));
        U(i) = 3*Ubar*(y/h-0.5*(y/h)^2);
    end
elseif base_flow == 2% 对应'time_average'
    for i = 1:length(U) 
        y = min(Y(i),2*h-Y(i));
        yplus = y*Retau/h;
        if yplus<=10.8
            U(i) = utau*yplus;
        elseif yplus>10.8
            U(i) = utau*(log(yplus)/0.41+5);
        end
    end
    Ubar = mean(U);
end

%% 扰动场参数设置
duplus = Ubar*0.25/utau;
sigma = 0.00055;
%% spanwise wavenumber: spacing z+ = 200
betaPlus = 2*pi*(1/200);
epsilon = Ubar/200;
%% streamwise wave number: spacing x+ = 500
alphaPlus = 2*pi*(1/500);
%% 生成扰动速度场
for i = 1:length(U)
    deviation = 1+0.2*randn(1);%%添加扰动
    xplus = X(i)*Retau/h;
    y = min(Y(i),2*h-Y(i));
    yplus = y*Retau/h;
    zplus = Z(i)*Retau/h;
    U(i) = U(i)+(utau*duplus/2.0) * (yplus/40.0)...
                * exp(-sigma * (yplus^2) + 0.5)...
                * cos(betaPlus*zplus)*deviation;
    W(i) = epsilon*sin(alphaPlus*xplus)*yplus*exp(-sigma*(yplus^2))...
              * deviation;
end
%% 对速度进行无量纲化
U = U./utau;V = V./utau;W = W./utau;

%% 绘制云图
U3d = reshape(U,nx,ny,nz);
uxz = reshape(U3d(:,5,:),nx,nz);
pcolor(uxz)
shading interp

%% 将数据转换为二进制格式（AoXu）
data = [U,V,W];
writeXuBinaryfile(strcat('../4_binary_data/MATLAB_Retau',num2str(Retau),'_Initial.bin'),data)
disp('successful write!')

%% 写入数据函数
function [] = writeXuBinaryfile(file,data)
% 写只有U，V，W信息的二进制文件
fid = fopen(file,'w');
fwrite(fid,1,'int32');
fwrite(fid,data(:,1),'float64');%U
fwrite(fid,[1 1],'int32');
fwrite(fid,data(:,2),'float64');%V
fwrite(fid,[1 1],'int32');
fwrite(fid,data(:,3),'float64');%W
fwrite(fid,1,'int32');
fclose(fid);
end