From 0a380a543f6607991dc294f4e2314c7906c96949 Mon Sep 17 00:00:00 2001 From: ignis Date: Mon, 28 Apr 2014 09:06:25 +0900 Subject: [PATCH] new model spectrum for initial velocity field --- init_velocity.f90 | 117 ++++++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 114 insertions(+), 3 deletions(-) diff --git a/init_velocity.f90 b/init_velocity.f90 index 45df2b3..ff2f690 100644 --- a/init_velocity.f90 +++ b/init_velocity.f90 @@ -26,6 +26,18 @@ subroutine init_velocity real*8 :: wmag, wmag2, ratio, fac, fac2 +!! New model spectrum--------------------------------------------------- +! real*8, allocatable :: model_e_spec(:) +! real*8 :: A,k0 +!!---------------------------------------------------See m_parameters.f90 + real*8 :: e_spec_log00=0. + real*8 :: e_spec_log01=0. + real*8 :: e_spec_log02=0. + real*8 :: e_spec_log03=0. + real*8 :: model_tke,model_eps_v,model_eta + real*8 :: model_d_spec(kmax) + real*8 :: kol_ts,kol_vs,rms_u_prime,model_sctmp,model_x_length,model_uvar + real*8 :: model_lambda,model_tau_e !-------------------------------------------------------------------------------- ! First, if it's a Taylor-Green vortex, then initialize and quit @@ -65,7 +77,7 @@ subroutine init_velocity !!$ call flush(out) seed1 = RN1 - write(out,*) "RANDOM SEED FOR VELOCITIES = ", seed1 + write(out,'(a30,3x,i10)') "RANDOM SEED FOR VELOCITIES = ", seed1 call flush(out) rseed = real(seed1,8) fac = random(-rseed) @@ -190,6 +202,98 @@ subroutine init_velocity !------------------------------------------------------------------------------- ! first, define the desired spectrum + +!--------------------------------------------------------J. Kwon. 201403 +! Model spectrum setup +! Ref. Flow, Turbulence and Combustion 73: 133-167, 2004 +! by L. Guichard, J. Reveillon and R. Hauguel +! ------------------------------------------------------------------------ + allocate(model_e_spec(kmax)) ; model_e_spec=zip + +! open(unit=404,file="model_spectrum.dat") +! write(404,*) 'VARIABLES = "X","Energy(k)"' ! +! write(404,*) '"Kolmogorov","Exponential","VonKarman"' ! + +! write(502,'(a15,f11.5,a1)')'ZONE T= " time=',time,'"' +! write(502,*) 'I=',kmax,'F=POINT' + + if (isp_type.eq.4) then !ksj + k0=(SQRT(two*PI))/t_length + model_A=16.0*(SQRT(two/PI))*(t_u_prime**2.0)/k0 + do i=1,kmax + k=real(i,8) + model_e_spec(i)=model_A*((k/k0)**4.0)*exp(-2.0*((k/k0)**2.0)) +! e_spec_log00=log10(model_e_spec(i)) + + wmag=real(k, 8) + ratio=wmag / peak_wavenum + e_spec_log01=wmag**(-5.d0/3.d0) + e_spec_log02=ratio**3 / peak_wavenum * exp(-3.0D0*ratio) + + fac = two * PI * ratio + e_spec_log03=fac**4 / (one + fac**2)**3 + + if(model_e_spec(i).le.0.) e_spec_log00=0. +! write(404,'(i3,4e15.5)') i,model_e_spec(i),& +! e_spec_log01,e_spec_log02,e_spec_log03 +! write(502,'(i5,e15.5)') i,model_e_spec(i) + enddo + endif + +! close(404) + +! write target stats at t=0 on tp_stats.dat + ! getting the total energy + model_tke = sum(model_e_spec(1:kmax)) + + ! finding dissipation spectrum and total dissipation + do k = 1,kmax + model_d_spec(k) = model_e_spec(k) * real(k**2,8) * two * nu + end do + model_eps_v = sum(model_d_spec(1:kmax)) + + ! finding Kolmogorov scale + model_eta = (nu**3/model_eps_v)**0.25 ! Kolmogorov length scale + + ! Kolmogorov time and velocity scale, J. Kwon + kol_ts=sqrt(nu/model_eps_v) ! Kolmogorov time scale(kol_ts) + kol_vs=model_eta/kol_ts ! Kolmogorov velocity scale(kol_vs) + + model_uvar=two/three*model_tke + ! rms velocity fluctuation, J. Kwon + rms_u_prime=sqrt(model_uvar) + + ! integral length scale(x_length) + model_sctmp = zip + do k = 1, kmax + model_sctmp = model_sctmp + model_e_spec(k) / real(k,8) + end do + model_x_length = PI / two * model_sctmp / model_uvar + + ! Taylor microscale + model_lambda = sqrt(15.d0 * model_uvar * nu / model_eps_v) + + ! Eddy turnover time + model_tau_e = model_x_length / sqrt(model_uvar) + + + ! outputting all this in the stat1 file +! inquire(file='tp_stat.dat', exist=there, opened=there2) +! if (.not.there) then +! open(69,file='tp_stat.dat',form='formatted') +! write(69,*) 'VARIABLES = "itime","time","tke","rms_u`","int_LS","Eddy_TO_time","Kol_LS"' ! +! write(69,*) '"Kol_VS","Kol_TS","Taylor_LS","dissipation"' +! end if +! if(there.and..not.there2) then +! open(69,file='tp_stat.dat',position='append') +! end if +! write(69,"(i10,f15.10,9e15.6)") itime,time,model_tke,rms_u_prime,model_x_length, & +! model_tau_e,model_eta,kol_vs,kol_ts,model_lambda,& +! model_eps_v + + +!---------------------------------------------------------J. Kwon. 201403 + do k = 1,kmax wmag = real(k, 8) @@ -208,6 +312,11 @@ subroutine init_velocity fac = two * PI * ratio e_spec1(k) = fac**4 / (one + fac**2)**3 + !ksj******************************************************* + else if (isp_type.eq.4) then + e_spec1(k) = model_e_spec(k) + !********************************************************** + else write(out,*) "ERROR: WRONG INITIAL SPECTRUM TYPE: ",isp_type call flush(out) @@ -216,8 +325,10 @@ subroutine init_velocity end if end do -! normalize it so it has the unit total energy - e_spec1 = e_spec1 / sum(e_spec1(1:kmax)) +!caution!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!! normalize it so it has the unit total energy +! e_spec1 = e_spec1 / sum(e_spec1(1:kmax)) +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! now go over all Fourier shells and multiply the velocities in a shell by ! the sqrt of ratio of the resired to the current spectrum