import scipy.io as sio import numpy as np import matplotlib.pyplot as plt from mpl_toolkits.axes_grid1.inset_locator import inset_axes #from IPython import display plt.rcParams.update({'font.size': 22}) # makes sure figures are updated when using ipython #display.clear_output(wait=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 # # # exemple of 1d Channel flow with a PDH k-omegaa model [1]. Re=u_tau*h/nu=8000 (h=half # channel height). # # Discretization described in detail in # http://www.tfd.chalmers.se/~lada/comp_fluid_dynamics/ # # [1} S-H Peng and L. Davidson and S. Holmberg", # "A Modified Low-{R}eynolds-Number k-omega Model for Recirculating Flows", # Journal of Fluids Engng, vol. 119, pp. 867-875, 1997 # max number of iterations niter=30000 # friction velocity u_*=1 # half channel width=1 # # create the grid nj=30 # coarse grid #nj=98 # fine grid njm1=nj-1 #yfac=1.6 # coarse grid yfac=1.15 # fine grid dy=0.1 yc=np.zeros(nj) delta_y=np.zeros(nj) yc[0]=0. for j in range(1,int(nj/2)): yc[j]=yc[j-1]+dy dy=yfac*dy ymax=yc[int(nj/2)-1] # cell faces for j in range(0,int(nj/2)): yc[j]=yc[j]/ymax yc[nj-1-j]=2.-yc[j-1] yc[nj-1]=2. # cell centres yp=np.zeros(nj) for j in range(1,nj-1): yp[j]=0.5*(yc[j]+yc[j-1]) yp[nj-1]=yc[nj-1] # viscosity viscos=1./8000. # under-relaxation urf=0.5 # plot k for each iteration at node jmon jmon=8 # turbulent constants sigma_k=0.8 sigma_om=1.35 betas=0.09 gama0 = 0.42 beta = 0.075 cr1=0.75 small=1.e-10 great=1.e10 # initialaze u=np.zeros(nj) k=np.ones(nj)*1.e-4 y=np.zeros(nj) om=np.ones(nj)*1. vist=np.ones(nj)*100.*viscos dn=np.zeros(nj) ds=np.zeros(nj) dy_s=np.zeros(nj) fy=np.zeros(nj) tau_w=np.zeros(niter) k_iter=np.zeros(niter) om_iter=np.zeros(niter) # do a loop over all nodes (except the boundary nodes) for j in range(1,nj-1): # compute dy_s dy_s[j]=yp[j]-yp[j-1] # compute dy_n dy_n=yp[j+1]-yp[j] # compute deltay delta_y[j]=yc[j]-yc[j-1] dn[j]=1./dy_n ds[j]=1./dy_s[j] # interpolation factor del1=yc[j]-yp[j] del2=yp[j+1]-yc[j] fy[j]=del1/(del1+del2) u[1]=0. vist[0]=0. vist[nj-1]=0. k[0]=0. k[nj-1]=0. su=np.zeros(nj) sp=np.zeros(nj) an=np.zeros(nj) as1=np.zeros(nj) ap=np.zeros(nj) # do max. niter iterations for n in range(1,niter): for j in range(1,nj-1): # compute turbulent viscosity vist_old=vist[j] ret=k[j]/(viscos*om[j]) ret=max(ret,1.e-5) fmu=0.025+(1.-np.exp(-(ret/10.0)**0.75))*(0.975+(1.e-3/ret)*np.exp(-(ret/200.)**2)) vist[j]=urf*fmu*k[j]/om[j]+(1.-urf)*vist_old # solve u for j in range(1,nj-1): # driving pressure gradient su[j]=delta_y[j] sp[j]=0. # interpolate turbulent viscosity to faces vist_n=fy[j]*vist[j+1]+(1.-fy[j])*vist[j] vist_s=fy[j-1]*vist[j]+(1.-fy[j-1])*vist[j-1] # compute an & as an[j]=(vist_n+viscos)*dn[j] as1[j]=(vist_s+viscos)*ds[j] # boundary conditions for u u[0]=0. u[nj-1]=0. for j in range(1,nj-1): # compute ap ap[j]=an[j]+as1[j]-sp[j] # under-relaxate ap[j]= ap[j]/urf su[j]= su[j]+(1.0-urf)*ap[j]*u[j] # use Gauss-Seidel u[j]=(an[j]*u[j+1]+as1[j]*u[j-1]+su[j])/ap[j] # use TDMA # u=tdma(ap,an,as,su,u,nj) # monitor the development of u_tau in node jmon tau_w[n]=viscos*u[1]/yp[1] # print iteration info print('Iteration number: ',n) # check for convergence (when converged, the wall shear stress must be one) if abs(tau_w[n]-1.) < 0.001: # do at least 1000 iter if n > 1000: print('Converged!') break # solve k dudy=np.gradient(u,yp) dudy2=dudy**2 for j in range(1,nj-1): # production term su[j]=vist[j]*dudy2[j]*delta_y[j] # dissipation term ret=k[j]/(viscos*om[j]) ret=max(ret,1.e-5) fk=1.-0.722*np.exp(-(ret/10.)**4) sp[j]=-fk*betas*om[j]*delta_y[j] # compute an & as vist_n=fy[j]*vist[j+1]+(1.-fy[j])*vist[j] an[j]=(vist_n/sigma_k+viscos)*dn[j] vist_s=fy[j-1]*vist[j]+(1.-fy[j-1])*vist[j-1] as1[j]=(vist_s/sigma_k+viscos)*ds[j] # boundary conditions for k k[0]=0. k[nj-1]=0. for j in range(1,nj-1): # compute ap ap[j]=an[j]+as1[j]-sp[j] # under-relaxate ap[j]= ap[j]/urf su[j]= su[j]+(1.-urf)*ap[j]*k[j] # use Gauss-Seidel k[j]=(an[j]*k[j+1]+as1[j]*k[j-1]+su[j])/ap[j] # use TDMA # k=tdma(ap,an,as,su,k,nj) # monitor the development of k in node jmon k_iter[n]=k[jmon] #****** solve om-eq. dkdy=np.gradient(k,yp) domdy=np.gradient(om,yp) for j in range(1,nj-1): # compute an & as vist_n=fy[j]*vist[j+1]+(1.-fy[j])*vist[j] an[j]=(vist_n/sigma_om+viscos)*dn[j] vist_s=fy[j-1]*vist[j]+(1.-fy[j-1])*vist[j-1] as1[j]=(vist_s/sigma_om+viscos)*ds[j] ret=k[j]/(viscos*om[j]) ret=max(ret,1.e-5) fmu=0.025+(1.-np.exp(-(ret/10.0)**0.75))*(0.975+(1.e-3/ret)*np.exp(-(ret/200.)**2)) fomega=1.+4.3*np.exp(-(ret/1.5)**0.5) gama = gama0*fomega # production term su[j]=gama*fmu*dudy2[j]*delta_y[j] # dissipation term sp[j]=-beta*om[j]*delta_y[j] # cross diffusion term crosv=dkdy[j]*domdy[j] crost=cr1*vist[j]*crosv/k[j]*delta_y[j] su[j]=su[j]+max(crost,0,) sp[j]=sp[j]+min(crost,0,)/om[j] # b.c. south wall dy=yp[1] omega=6.*viscos/0.075/dy**2 sp[1]=-great su[1]=great*omega # b.c. north wall dy=yp[nj-1]-yp[nj-2] omega=6.*viscos/0.075/dy**2 sp[nj-2]=-great su[nj-2]=great*omega for j in range(1,nj-1): # compute ap ap[j]=an[j]+as1[j]-sp[j] # under-relaxate ap[j]= ap[j]/urf su[j]= su[j]+(1.-urf)*ap[j]*om[j] # use Gauss-Seidel om[j]=(an[j]*om[j+1]+as1[j]*om[j-1]+su[j])/ap[j] # use TDMA # om=tdma(ap,an,as,su,om,nj) om_iter[n]=om[jmon] # compute shear stress uv=-vist*dudy #./yp_u_k_om_vist_uv_yc_PDH_8000.dat # plot u fig1 = plt.figure("Figure 1") plt.clf() #clear the figure plt.plot(u,yp,'b-') plt.xlabel('u') plt.ylabel('y') plt.savefig('u_8000_python.eps') # plot k fig2 = plt.figure("Figure 2") plt.clf() #clear the figure plt.plot(k,yp,'b-') plt.xlabel('k') plt.ylabel('y') plt.savefig('k_8000_python.eps') # plot tau_w versus iteration number fig3 = plt.figure("Figure 3") plt.clf() #clear the figure plt.plot(tau_w,'b-') plt.title('wall shear stress') plt.xlabel('Iteration number') plt.ylabel('tauw') plt.savefig('tauw_python.eps') # plot k(jmon) versus iteration number fig4 = plt.figure("Figure 4") plt.clf() #clear the figure plt.plot(k_iter,'b-') plt.title('k in node jmon') plt.xlabel('Iteration number') plt.ylabel('k') plt.savefig('k_iter_python.eps') # plot om(jmon) versus iteration number fig5 = plt.figure("Figure 5") plt.clf() #clear the figure plt.plot(om_iter,'b-') plt.title('omega in node jmon') plt.xlabel('Iteration number') plt.ylabel('omega') plt.savefig('om_iter_python.eps') # save data data=np.zeros((nj,7)) data[:,0]=yp data[:,1]=u data[:,2]=k data[:,3]=om data[:,4]=vist data[:,5]=uv data[:,6]=yc np.savetxt('yp_u_k_om_vist_uv_yc_PDH_8000_python.dat', data)