import scipy.io as sio import numpy as np import time from angles import angles import matplotlib.pyplot as plt import math from mpl_toolkits.axes_grid1.inset_locator import inset_axes 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 wew1fct=2 visc=1./8000. qm=4 delta=0.5 # we assume that we have a boundary layer thickness of 0.5 sli=0.3 # make the length scaler one tenth of delta up=1. epsm=up**3/sli # #%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% # number of grid points in y, z nj=97 # number of grid point in y direction nk=33 # number of grid point in z direction # generate x-grid xp=np.zeros(2) xp[0]=0 xp[1]=1 # generate z-grid dz=1.6/(nk-1) zp=np.zeros(nk) zp[0]=0 for k in range(1,nk): zp[k]=zp[k-1]+dz # load RANS data created by rans.m (which can be downloaded) data = np.loadtxt("yp_u_k_om_vist_uv_yc_PDH_8000.dat") # viscosity nu=visc y=data[:,0] u=data[:,1] rk_turb=data[:,2] om=data[:,3] uv_rans=data[:,5] yc1=data[:,6] ed=0.09*rk_turb*om yp=yc1; # here yp means face coordinates # search min grid step dminy=min(np.diff(yp[0:-1-1])) dminz=min(np.diff(zp)) dxmin=min(dminy,dminz) # don't let the minimum distance be too small. This will create unnecessay high wavenumbers dxmin=max(dxmin,1.7e-3) # create a seed from time (must be negative) sec=int(time.time()) iseed=-sec # generate nstep realizations nstep=200 dt=np.zeros(nstep+1) for i in range(1,nstep): dt[i]=1e-3 ustar=(nu*u[1]/y[1])**0.5 dudy=np.gradient(u,y) nj_rans=len(rk_turb) # set synth. fluct at old time step to zero uprim_old=np.zeros((nj,nk)) vprim_old=np.zeros((nj,nk)) wprim_old=np.zeros((nj,nk)) uu=np.zeros(nj) vv=np.zeros(nj) ww=np.zeros(nj) uv=np.zeros(nj) # lowe Re EARSM by # S. Wallin and A. V. Johansson, #"A New Explicit Algebraic Reynolds Stress Model for Incompressible #and Compressible Turbulent Flows", JFM, #Vol. 403, pp. 89-132, 2000 for j in range (0,nj-1): diss=ed[j] rk=rk_turb[j] ttau=rk/diss om12=ttau*0.5*dudy[j] om21=-om12 om22=0. om11=0. s11=0. s12=ttau*0.5*dudy[j] s21=s12 s22=0. s33=0. vor=(-2.*om12**2) str1=(s11**2+s12**2+s21**2+s22**2) cpr1=9./4.*(1.8-1.) p1=(1./27.*cpr1**2+9./20.*str1-2./3.*vor)*cpr1 p2=p1**2-(1./9.*cpr1**2+9./10.*str1+2./3.*vor)**3 if p2 > 0: if p1-p2**0.5 >= 0: sigg=1. else: sigg=-1. un=cpr1/3.+(p1+p2**0.5)**(1./3.)+sigg*(abs(p1-p2**0.5))**(1./3.) else: un=cpr1/3.+2.*(p1**2-p2)**(1./6.)*np.cos(1./3.*np.arccos(p1/(np.sqrt(p1**2-p2)))) const=6./5. beta1=-const*un/(un**2-2.*vor) beta4=beta1/un uu[j]=2./3.*rk+rk*beta1*s11+rk*beta4*(s12*om21-om12*s21) vv[j]=2./3.*rk+rk*beta1*s22+rk*beta4*(s21*om12-om21*s12) ww[j]=2./3.*rk uv[j]=rk*(beta1*s12+beta4*(s11*om12-om12*s22)) # max(uv) at j=19 j=18 stress=[[uu[j],uv[j],0.], \ [uv[j],vv[j],0.], \ [0.,0.,ww[j]]] diag_sum=np.trace(stress)/3. stress=stress/diag_sum lambd,V =np.linalg.eig(stress) v1=[V[0,0],V[1,0], V[2,0]] v2=[V[0,1],V[1,1], V[2,1]] v3=[V[0,2],V[1,2], V[2,2]] lambda_1=lambd[0] lambda_2=lambd[1] lambda_3=lambd[2] # largest eigenvalue in x1*, smallest in x2*, i.e. OK v1_new=v1 v2_new=v2 v3_new=v3 lambda_1_new=lambda_1 lambda_2_new=lambda_2 lambda_3_new=lambda_3 # 1st eigenvector in 1st or 3rd quadrant # 2nd eigenvector in 2nd or 4th quadrant # switch sign on 12 and 21 to fix the above requirements v1_new[0]=-v1_new[0] v2_new[1]=-v2_new[1] r11=v1_new[0] r12=v1_new[1] r13=v1_new[2] r21=v2_new[0] r22=v2_new[1] r23=v2_new[2] r31=v3_new[0] r32=v3_new[1] r33=v3_new[2] a11=lambda_1_new a22=lambda_2_new a33=lambda_3_new; print('eigenvector 1= ',r11,r21,r31) print('eigenvector 2= ',r12,r22,r32) print('eigenvector 3= ',r13,r23,r33) print('eigenvalue 1= ',a11) print('eigenvalue 2= ',a22) print('eigenvalue 3= ',a33) # 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('total 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] # invert the eigenvalue matrix (anisotropic) a11i=1./a11 a22i=1./a22 a33i=1./a33 # 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]) xcp=r11*xc+r21*yc+r31*zc ycp=r12*xc+r22*yc+r32*zc zcp=r13*xc+r23*yc+r33*zc # initiate turbulent velocities to zero. utrp=0. vtrp=0. wtrp=0. # loop over all wavenumbers for m in range(1,nmodes+1): # r21, r31, r12, r32, r13, r23 = 0 in isotropic turbulence kxi=r11*kxio[m]+r21*kyio[m]+r31*kzio[m] kyi=r12*kxio[m]+r22*kyio[m]+r32*kzio[m] kzi=r13*kxio[m]+r23*kyio[m]+r33*kzio[m] sxi=r11*sxio[m]+r21*syio[m]+r31*szio[m] syi=r12*sxio[m]+r22*syio[m]+r32*szio[m] szi=r13*sxio[m]+r23*syio[m]+r33*szio[m] # a11, a22, a33=1 in isotropic turbulence sx=np.sqrt(a11)*sxi sy=np.sqrt(a22)*syi sz=np.sqrt(a33)*szi kx=np.sqrt(a11i)*kxi*wnr[m] ky=np.sqrt(a22i)*kyi*wnr[m] kz=np.sqrt(a33i)*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*xcp+ky*ycp+kz*zcp+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): utrpn=r11*utrp+r12*vtrp+r13*wtrp; vtrpn=r21*utrp+r22*vtrp+r23*wtrp; wtrpn=r31*utrp+r32*vtrp+r33*wtrp; # store the fluctuations u, v, w in 2D arrays u[j,k]=utrpn v[j,k]=vtrpn w[j,k]=wtrpn # 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[:],ut[:,1],'b-') plt.xlabel('z') plt.ylabel('u') plt.matplotlib.pyplot.show() plt.savefig('u_inst_python.ps') plt.show()