import scipy.io as sio import numpy as np import time from angles import angles import matplotlib.pyplot as plt import math plt.rcParams.update({'font.size': 22}) plt.interactive(True) # This file can be downloaded at # http://www.tfd.chalmers.se/~lada/projects/proind.html # choose "Synthetic Inlet Boundary Conditions" on the left side # # # Litterature # # L. Davidson, "Using Isotropic Synthetic Fluctuations as Inlet Boundary Conditions for Unsteady # Simulations", "Advances and Applications in Fluid Mechanics", # Vol 1, No =1, pp. 1-35, 2007 # # L. Davidson, "HYBRID LES-RANS: Inlet Boundary Conditions for Flows With Recirculation", # "Advances in Hybrid RANS-LES Modelling", # Notes on Numerical Fluid Mechanics and Multidisciplinary Design, # Springer Verlag, pp. 55-66, Vol. 97, 2008. # # L. Davidson, "Hybrid LES-RANS: Inlet Boundary Conditions for Flows Including Recirculation", # 5th International Symposium on Turbulence and Shear Flow Phenomena, # Vol. 2, pp. 689-694, Munich, Germany, 2007. # # M. Billson and L.-E. Eriksson and L. Davidson, "Modeling of Synthetic Anisotropic # Turbulence and its Sound Emission", The 10th AIAA/CEAS Aeroacoustics Conference, # AIAA 2004-2857, Manchester, United Kindom, 2004. # # M. Billson, "Computational Techniques for Turbulence Generated Noise", # Dept. of Thermo and Fluid Dynamics, Chalmers University of Technology", # Göteborg, Sweden, 2004. # # L. Davidson and S.-H. Peng, "Embedded Large-Eddy Simulation Using the Partially # Averaged Navier-Stokes Model", AIAA J, vol. 51, number=5, pp=1066-1079, # doi="doi: 10.2514/1.J051864", 2013. #=========================== chapter 1 ============================================ #!!! number of modes = nmodes #!!! smallest wavenumber = dxmin #!!! ratio of ke and kmin (in wavenumber) = wew1fct #!!! turb. velocity scale = up #!!! turb. kin. energy = qm #!!! diss. rate. = epsm #!!! kinetic viscosity = visc #!!! length scale = sli amp=1.452762113 nmodes =150 #nmodes =10 wew1fct=2 visc=15.2e-6 qm=4 delta=0.5 # we assume that we have a boundary layer thickness of 0.5 sli=0.1*delta # make the length scaler one tenth of delta up=np.sqrt(2.*qm/3.) epsm=qm**1.5/sli # #%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% # number of grid points in y, z nj=2 # number of grid point in y direction nk=63 # number of grid point in z direction # generate x-grid xp=np.zeros(2) xp[0]=0 xp[1]=1 # generate z-grid dz=0.004 zp=np.zeros(nk) zp[0]=0 for k in range(1,nk): zp[k]=zp[k-1]+dz # generate y-grid yp=np.zeros(nj) yp[0]=0. yp[1]=1. # search min grid step dminy=min(np.diff(yp)) dminz=min(np.diff(zp)) dxmin=min(dminy,dminz) # create a seed from time (must be negative) sec=int(time.time()) iseed=-sec iseed=-3 # generate nstep realizations nstep=200 #nstep=1 # zero all arrays to zero wnr=np.zeros(nmodes+2) fi=np.zeros(nmodes+2) teta=np.zeros(nmodes+2) psi=np.zeros(nmodes+2) wnrf=np.zeros(nmodes+2) wnr=np.zeros(nmodes+2) dkn=np.zeros(nmodes+2) kxio=np.zeros(nmodes+2) kyio=np.zeros(nmodes+2) kzio=np.zeros(nmodes+2) sxio=np.zeros(nmodes+2) syio=np.zeros(nmodes+2) szio=np.zeros(nmodes+2) u=np.zeros((nj,nk)) v=np.zeros((nj,nk)) w=np.zeros((nj,nk)) ut=np.zeros((nk,nstep+1)) vt=np.zeros((nk,nstep+1)) wt=np.zeros((nk,nstep+1)) iy=0 iv=np.zeros(33) # #=========================== chapter 2 ============================================ # print('time step no: ',nstep) for it in range(1,nstep+1): # compute random angles fi,psi,alfa,teta,iy,iv,iseed = angles(nmodes,iseed,iy,iv) print('time step no: ',it) #%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% # highest wave number wnrn=2.*math.pi/dxmin # # k_e (related to peak energy wave number) wnre=9.*math.pi*amp/(55.*sli) # # wavenumber used in the viscous expression (high wavenumbers) in the von Karman spectrum wnreta=(epsm/visc**3)**0.25 # smallest wavenumber wnr1=wnre/wew1fct # wavenumber step dkl=(wnrn-wnr1)/nmodes # wavenumber at faces for m in range(1,nmodes+2): wnrf[m]=wnr1+dkl*(m-1) # wavenumber at cell centers for m in range(1,nmodes+1): wnr[m]=.5*(wnrf[m+1]+wnrf[m]) dkn[m]=wnrf[m+1]-wnrf[m] # wavenumber vector from random angles for m in range(1,nmodes+1): kxio[m]=np.sin(teta[m])*np.cos(fi[m]) kyio[m]=np.sin(teta[m])*np.sin(fi[m]) kzio[m]=np.cos(teta[m]) # # sigma (s=sigma) from random angles. sigma is the unit direction which gives the direction # of the synthetic velocity vector (u, v, w) sxio[m]=np.cos(fi[m])*np.cos(teta[m])*np.cos(alfa[m])-np.sin(fi[m])*np.sin(alfa[m]) syio[m]=np.sin(fi[m])*np.cos(teta[m])*np.cos(alfa[m])+np.cos(fi[m])*np.sin(alfa[m]) szio[m]=-np.sin(teta[m])*np.cos(alfa[m]) # #=========================== chapter 3 ============================================ # # loop through grid for k in range(0,nk-1): for j in range(0,nj-1): # cell center coordinates xc=0.5 yc=0.5*(yp[j]+yp[j+1]) zc=0.5*(zp[k]+zp[k+1]) # initiate turbulent velocities to zero. utrp=0. vtrp=0. wtrp=0. # loop over all wavenumbers for m in range(1,nmodes+1): kxi=kxio[m] kyi=kyio[m] kzi=kzio[m] sx=sxio[m] sy=syio[m] sz=szio[m] kx=kxi*wnr[m] ky=kyi*wnr[m] kz=kzi*wnr[m] rk=np.sqrt(kx**2+ky**2+kz**2) # if the wavenumber, rk, is smaller than the largest wavenumber, then create fluctuations if rk < wnrn: arg=kx*xc+ky*yc+kz*zc+psi[m] tfunk=np.cos(arg) # von Karman spectrum e=amp/wnre*(wnr[m]/wnre)**4/((1.+(wnr[m]/wnre)**2)**(17./6.))*np.exp(-2*(wnr[m]/wnreta)**2) utn=np.sqrt(e*up**2*dkn[m]) # synthetic velocity field utrp=utrp+2.*utn*tfunk*sx vtrp=vtrp+2.*utn*tfunk*sy wtrp=wtrp+2.*utn*tfunk*sz # end if rk < wnrn # end for m in range(0,nmodes-1): # store the fluctuations u, v, w in 2D arrays u[j,k]=utrp v[j,k]=vtrp w[j,k]=wtrp # store the fluctuations for each time step it (we have only one cell in y-dirrection j=1) # these will be used for plotting and evaluation below ut[k,it]=u[j,k] vt[k,it]=v[j,k] wt[k,it]=w[j,k] # end for k in range(0,nk-2) # ebd for j in range(0,nj-2) # end for it in range(1,nstep): print('synthetic fluctuations created') # #=========================== chapter 4 ============================================ # umean=np.mean(ut,axis=(0,1)); vmean=np.mean(vt,axis=(0,1)) wmean=np.mean(wt,axis=(0,1)) print('umean=',umean) print('vmean=',vmean) print('wmean=',wmean) print('program stop') #**************** time history ************************************* #% plot u of one realization (time step 5) # fig1 = plt.figure("Figure 1") plt.plot(zp[0:nk-1],ut[0:nk-1,1],'b-') plt.xlabel('z') plt.ylabel('u') plt.savefig('u_inst_python.eps') plt.show()