new model spectrum for initial velocity field
This commit is contained in:
ignis
parent
89a77a9f32
commit
0a380a543f
1 changed files with 114 additions and 3 deletions
|
|
@ -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
|
||||
|
|
|
|||
Loading…
Add table
Reference in a new issue