import scipy.io as sio import numpy as np import time from angles import angles import matplotlib.pyplot as plt import math plt.interactive(True) plt.rcParams.update({'font.size': 22}) # 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 # the Reynols stress tensor it taken DNS # % Authors: Moser, Kim & Mansour # Reference: DNS of Turbulent Channel Flow up to Re_tau=590, 1999, # Physics of Fluids, vol 11, 943-945. # # y y+ R_uu R_vv R_ww R_uv R_uw R_vw # #2.7060e-02 1.5890e+01 7.6684e+00 3.1974e-01 1.4997e+00 -6.6206e-01 1.7267e-03 8.7423e-04 stress=[[7.6684,-6.6206e-01,0.], \ [-6.6206e-01,3.1974e-01,0.], \ [0.,0.,1.4997]] 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('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[0:nk-1],ut[0:nk-1,1],'b-') plt.xlabel('z') plt.ylabel('u') plt.savefig('u_inst_python.eps') plt.show()