474 lines
13 KiB
Fortran
474 lines
13 KiB
Fortran
MODULE ysolve
|
|
USE Compact
|
|
IMPLICIT NONE
|
|
PRIVATE
|
|
REAL, PARAMETER :: pi=3.14159265358979323846
|
|
REAL :: hx,dt,vis,sc,diff,pre,ac,bc,tf,t_now,t_uf,dt_uf
|
|
REAL :: rvis,cstar
|
|
REAL :: c_cut,c_ref,refwr,minf,tar_lo,u0,ctmp,l_0,lo_flm=0.
|
|
REAL :: er_lof=0., erdot=0.,min_wr,prof_wr
|
|
REAL :: pflame,pflold,oldu
|
|
INTEGER :: ncyc=0,int_pr,nx
|
|
REAL, DIMENSION(:,:,:), ALLOCATABLE :: u,y1,y2,yf
|
|
REAL, DIMENSION(:), ALLOCATABLE :: uxt,duxt
|
|
REAL, DIMENSION(:), ALLOCATABLE :: dm
|
|
|
|
INTEGER :: istage
|
|
REAL, DIMENSION(5) :: a=(/ 970286171893./4311952581923., &
|
|
6584761158862./12103376702013., &
|
|
2251764453980./15575788980749., &
|
|
26877169314380./34165994151039., &
|
|
0. /), &
|
|
b=(/ 1153189308089./22510343858157., &
|
|
1772645290293./4653164025191., &
|
|
-1672844663538./4480602732383., &
|
|
2114624349019./3568978502595., &
|
|
5198255086312./14908931495163. /)
|
|
|
|
PUBLIC :: solve
|
|
CONTAINS
|
|
!------------------------------------------------------------------------
|
|
SUBROUTINE solve
|
|
INTEGER :: i,j,k,savenum
|
|
REAL :: pflame,pflold,delf=0.
|
|
|
|
CALL READ_INTRO
|
|
|
|
CALL ludcmp(nx,5,5,1,0,0)
|
|
|
|
CALL SET_IC
|
|
|
|
t_now=0.; t_uf=0.
|
|
|
|
DO
|
|
CALL SET_BC
|
|
CALL RK4
|
|
IF(t_now.ge.tf) EXIT
|
|
IF(t_uf.ge.dt_uf) THEN
|
|
pflold=pflame
|
|
pflame=0.
|
|
DO i=1,nx
|
|
pflame=pflame+y1(i,1,1)*hx
|
|
uxt(i)=y1(i,1,1)
|
|
ENDDO
|
|
CALL dfnonp(nx,hx,uxt,duxt,1,1)
|
|
delf=1./MAXVAL(ABS(duxt))
|
|
oldu=u(1,1,1)
|
|
u=u+0.5*(hx*REAL(nx)*tar_lo-pflame)+0.5*(pflold-pflame)
|
|
t_uf=0.
|
|
WRITE(*,'(a3,f8.3,a10,f6.3,a10,f6.3,a10,f7.4,a10,f7.4,a8,f7.4)') &
|
|
' T:',t_now,' // Tar_L:',l_0*tar_lo,' // cur_L:',pflame/hx/REAL(nx)*l_0, &
|
|
' // Old_U:',oldu,' // New_U:',u(1,1,1),' // L_f:',delf
|
|
|
|
! WRITE(*,'(a7,f7.4,a7,f7.4,a10,f7.4)') ' cur_U:',oldu,' // dU:',u(1,1,1)-oldu,' // new_U:',u(1,1,1)
|
|
! WRITE(*,*)
|
|
ENDIF
|
|
ncyc=ncyc+1
|
|
t_uf=t_uf+dt
|
|
t_now=t_now+dt
|
|
IF (MOD(ncyc,int_pr).eq.0) THEN
|
|
! WRITE(*,'(a2,f8.3,a9,f10.7,a11,i6,a7,f9.5)') &
|
|
! 'T:',t_now,' // dT:',dt,' // NCYC:',ncyc,' // U:',u(1,1,1)
|
|
ENDIF
|
|
ENDDO
|
|
|
|
CALL write_sd
|
|
CALL write_pre
|
|
CALL save_final_field
|
|
|
|
WRITE(*,*)
|
|
WRITE(*,*) 'Fin.'
|
|
WRITE(*,*)
|
|
|
|
END SUBROUTINE solve
|
|
|
|
SUBROUTINE write_sd
|
|
REAL :: c,yr,wrate,dely,sdr,sdd,sd,uu,onelw,dd
|
|
INTEGER :: i,j,k,nd
|
|
REAL, DIMENSION(1,nx) :: ux,dux,d2ux
|
|
REAL, DIMENSION(1,nx) :: duy
|
|
REAL, DIMENSION(10,nx) :: sav
|
|
|
|
sav=0.
|
|
|
|
! refwr=pre*1.*exp(-ac/(1.+bc*c_cut))
|
|
! minf=exp((c_ref-c_cut)*prof_wr)
|
|
refwr=pre*1.*exp(-ac/(1.+bc*c_ref))
|
|
|
|
|
|
WRITE(500,*) 'VARIABLES = "X","Yr","C","U","Wrate","|DEL(Y)|","Sdr"'
|
|
WRITE(500,*) ' "Sdd","Sd","(1/C)/(dC/dx)","DIV(rho*Dmu*Gra(C))"'
|
|
|
|
DO i=1,nx
|
|
ux(1,i)=y1(i,1,1) ! Yr
|
|
IF (ux(1,i).gt.1.) ux(1,i)=1.
|
|
ENDDO
|
|
|
|
dm = diff * ((1.0 + bc * (1.0 - y1(:,1,1))) ** 0.76)
|
|
|
|
nd=1
|
|
CALL dfnonp(nx,hx,ux(1,:),dux(1,:),nd,1)
|
|
CALL dfnonp(nx,hx,dm,duy(1,:),nd,1)
|
|
nd=1
|
|
CALL d2fnonp(nx,hx,ux(1,:),d2ux(1,:),nd,1)
|
|
|
|
|
|
DO i=1,nx
|
|
yr=ux(1,i)
|
|
c=1.-yr
|
|
IF (c.lt.0.) c=0.
|
|
|
|
wrate=pre*yr*exp(-ac/(1.+bc*(1.-yr))) !wrate
|
|
|
|
! IF (c.le.c_cut) THEN
|
|
! wrate=min_wr
|
|
! IF (c.gt.c_ref) wrate= &
|
|
! ((exp((c-c_cut)*prof_wr)-minf)/(1.-minf)*(refwr-min_wr))+min_wr
|
|
! ENDIF
|
|
|
|
IF (c.le.c_ref) THEN
|
|
wrate=min_wr
|
|
IF (c.gt.c_cut) wrate=((refwr-min_wr)*exp(prof_wr*(c-c_ref))+ &
|
|
min_wr-refwr*exp(prof_wr*(c_cut-c_ref)))/(1.-exp(prof_wr*(c_cut-c_ref)))
|
|
ENDIF
|
|
|
|
dely=ABS(dux(1,i))
|
|
sdr=wrate/dely
|
|
sdd=-(dm(i)*d2ux(1,i)+duy(1,i)*dux(1,i))/dely
|
|
IF (dely.eq.0.) THEN
|
|
sdr=0.; sdd=0.
|
|
ENDIF
|
|
sd=sdr+sdd
|
|
uu=u(1,1,1)
|
|
onelw=(-dux(1,i))/c
|
|
dd=-(dm(i)*d2ux(1,i)+duy(1,i)*dux(1,i))
|
|
if (c.eq.0.) onelw=0.
|
|
|
|
sav(1,i)=sav(1,i)+yr
|
|
sav(2,i)=sav(2,i)+c
|
|
sav(3,i)=sav(3,i)+uu
|
|
sav(4,i)=sav(4,i)+wrate
|
|
sav(5,i)=sav(5,i)+dely
|
|
sav(6,i)=sav(6,i)+sdr
|
|
sav(7,i)=sav(7,i)+sdd
|
|
sav(8,i)=sav(8,i)+sd
|
|
sav(9,i)=sav(9,i)+onelw
|
|
sav(10,i)=sav(10,i)+dd
|
|
ENDDO
|
|
|
|
DO i=1,nx
|
|
WRITE(500,'(37e30.20)') (i-1)*hx,sav(1:10,i)
|
|
ENDDO
|
|
END SUBROUTINE write_sd
|
|
|
|
SUBROUTINE save_final_field
|
|
INTEGER :: i,j,k
|
|
|
|
OPEN (305,FILE='sfield.bin',form='unformatted',status='unknown')
|
|
DO i=1,nx
|
|
WRITE (305) y1(i,1,1)
|
|
ENDDO
|
|
CLOSE (305)
|
|
|
|
OPEN (305,FILE='sfield.dat')
|
|
DO i=1,nx
|
|
WRITE (305,'(e30.20)') y1(i,1,1)
|
|
ENDDO
|
|
CLOSE (305)
|
|
|
|
END SUBROUTINE save_final_field
|
|
|
|
SUBROUTINE write_pre
|
|
REAL :: yr,c,dy,maxdy=0.,del_f
|
|
REAL :: S_L=0.,wrate
|
|
INTEGER :: i
|
|
REAL, DIMENSION(1,nx) :: ux, dux
|
|
|
|
DO i=1,nx
|
|
yr=y1(i,1,1)
|
|
c=1.-yr
|
|
ux(1,i)=yr
|
|
|
|
wrate=pre*yr*exp(-ac/(1.+bc*(1.-yr))) !wrate
|
|
! IF (c.le.c_cut) THEN
|
|
! wrate=min_wr
|
|
! IF (c.gt.c_ref) wrate= &
|
|
! ((exp((c-c_cut)*prof_wr)-minf)/(1.-minf)*(refwr-min_wr))+min_wr
|
|
! ENDIF
|
|
IF (c.le.c_ref) THEN
|
|
wrate=min_wr
|
|
IF (c.gt.c_cut) wrate=((refwr-min_wr)*exp(prof_wr*(c-c_ref))+ &
|
|
min_wr-refwr*exp(prof_wr*(c_cut-c_ref)))/(1.-exp(prof_wr*(c_cut-c_ref)))
|
|
ENDIF
|
|
|
|
S_L=S_L+wrate*hx
|
|
ENDDO
|
|
|
|
WRITE(*,'(a31,e14.8)') ' INTEGRAL( Wrate x dx ) => Sc :',S_L
|
|
|
|
CALL dfnonp(nx,hx,ux(1,:),dux(1,:),1,1)
|
|
|
|
DO i=(nx/10),nx-(nx/10)
|
|
dy=dux(1,i)
|
|
maxdy=MAX(ABS(dy),maxdy)
|
|
ENDDO
|
|
del_f=1./maxdy
|
|
|
|
WRITE(*,'(a13,e14.8,a25,e14.8)') ' Grid size : ',hx,' / Laminar flame speed : ',u(1,1,1)
|
|
WRITE(*,'(a19,e14.8,a3,f9.5,a20)') &
|
|
' Flame thickness : ',del_f,' / ',del_f/hx,' grids in the flame.'
|
|
WRITE(*,*)
|
|
END SUBROUTINE write_pre
|
|
|
|
SUBROUTINE SET_BC
|
|
|
|
y1(1,1,1)=1.-ctmp
|
|
END SUBROUTINE SET_BC
|
|
|
|
SUBROUTINE SET_IC
|
|
INTEGER :: i, ifl, si
|
|
REAL :: xi
|
|
|
|
ifl=INT(2.0*pi*0.1/hx)
|
|
|
|
u=u0; si=INT(nx*(1.-tar_lo))
|
|
|
|
DO i=1,nx
|
|
IF(i< nx-(si+ifl/2)) THEN
|
|
xi=0.+ctmp
|
|
ELSE IF(i> nx-(si-ifl/2)) THEN
|
|
xi=1.
|
|
ELSE
|
|
xi=0.5+REAL(i-nx+si)/REAL(ifl)
|
|
ENDIF
|
|
|
|
! y1(i,1,1)=(1.-xi)*1.
|
|
y1(i,1,1)=(1.-xi) ! reactant mass fraction
|
|
ENDDO
|
|
|
|
pflame=0.
|
|
DO i=1,nx
|
|
pflame=pflame+y1(i,1,1)*hx
|
|
ENDDO
|
|
|
|
END SUBROUTINE SET_IC
|
|
|
|
SUBROUTINE READ_INTRO
|
|
CHARACTER(LEN=8) :: cdum
|
|
INTEGER :: itape=5,otape=301,ierr
|
|
|
|
|
|
! OPEN(itape,FILE='itape')
|
|
OPEN(otape,FILE='otape')
|
|
READ(itape,*) cdum,nx
|
|
WRITE(otape,*) cdum,nx
|
|
READ(itape,*) cdum,l_0
|
|
WRITE(otape,*) cdum,l_0
|
|
READ(itape,*) cdum,int_pr
|
|
WRITE(otape,*) cdum,int_pr
|
|
READ(itape,*) cdum,tar_lo
|
|
WRITE(otape,*) cdum,tar_lo
|
|
READ(itape,*) cdum,dt
|
|
WRITE(otape,*) cdum,dt
|
|
READ(itape,*) cdum,sc
|
|
WRITE(otape,*) cdum,sc
|
|
READ(itape,*) cdum,vis
|
|
WRITE(otape,*) cdum,vis
|
|
READ(itape,*) cdum,pre
|
|
WRITE(otape,*) cdum,pre
|
|
READ(itape,*) cdum,ac
|
|
WRITE(otape,*) cdum,ac
|
|
READ(itape,*) cdum,bc
|
|
WRITE(otape,*) cdum,bc
|
|
READ(itape,*) cdum,u0
|
|
WRITE(otape,*) cdum,u0
|
|
READ(itape,*) cdum,tf
|
|
WRITE(otape,*) cdum,tf
|
|
READ(itape,*) cdum,dt_uf
|
|
WRITE(otape,*) cdum,dt_uf
|
|
READ(itape,*) cdum,ctmp
|
|
WRITE(otape,*) cdum,ctmp
|
|
READ(itape,*) cdum,c_cut
|
|
WRITE(otape,*) cdum,c_cut
|
|
! READ(itape,*) cdum,cs
|
|
! WRITE(otape,*) cdum,cs
|
|
READ(itape,*) cdum,c_ref
|
|
WRITE(otape,*) cdum,c_ref
|
|
READ(itape,*) cdum,min_wr
|
|
WRITE(otape,*) cdum,min_wr
|
|
READ(itape,*) cdum,prof_wr
|
|
WRITE(otape,*) cdum,prof_wr
|
|
READ(itape,*) cdum,rvis
|
|
WRITE(otape,*) cdum,rvis
|
|
READ(itape,*) cdum,cstar
|
|
WRITE(otape,*) cdum,cstar
|
|
|
|
! CLOSE(itape)
|
|
|
|
! hx=l_0*pi/REAL(nx)
|
|
hx=l_0*pi/REAL(nx-1)
|
|
cdum='hx'
|
|
WRITE(otape,*) cdum,hx
|
|
cdum='Ta/Tu'
|
|
WRITE(otape,*) cdum,ac
|
|
cdum='Tb/Tu'
|
|
WRITE(otape,*) cdum,bc+1
|
|
diff=vis/sc
|
|
cdum='diff'
|
|
WRITE(otape,*) cdum,diff
|
|
|
|
! refwr=pre*1.*exp(-ac/(1.+bc*c_cut))
|
|
! minf=exp((c_ref-c_cut)*prof_wr)
|
|
|
|
refwr=pre*1.*exp(-ac/(1.+bc*c_ref))
|
|
l_0=l_0*pi
|
|
|
|
ALLOCATE(u(nx,1,1),STAT=ierr) ; u=0.
|
|
ALLOCATE(y1(nx,1,1),STAT=ierr) ; y1=0.
|
|
ALLOCATE(y2(nx,1,1),STAT=ierr) ; y2=0.
|
|
ALLOCATE(yf(nx,1,1),STAT=ierr) ; yf=0.
|
|
ALLOCATE(uxt(nx),STAT=ierr) ; uxt=0.
|
|
ALLOCATE(duxt(nx),STAT=ierr) ; duxt=0.
|
|
ALLOCATE(dm(nx),STAT=ierr) ; dm=diff
|
|
|
|
END SUBROUTINE READ_INTRO
|
|
|
|
SUBROUTINE RK4
|
|
istage=1; CALL substep(y1,y1,y2,yf)
|
|
istage=2; CALL substep(y1,y2,y1,yf)
|
|
istage=3; CALL substep(y2,y1,y2,yf)
|
|
istage=4; CALL substep(y1,y2,y1,yf)
|
|
istage=5; CALL substep(y2,y1,y2,yf)
|
|
END SUBROUTINE RK4
|
|
!------------------------------------------------------------------------
|
|
SUBROUTINE substep(ri,r1,r2,f)
|
|
REAL, INTENT(INOUT),DIMENSION(:,:,:) :: ri,r1,r2
|
|
REAL, INTENT(OUT),DIMENSION(:,:,:) :: f
|
|
|
|
INTEGER :: i,j,k
|
|
REAL :: at,bt
|
|
|
|
CALL fns(ri,f)
|
|
IF(istage<5) THEN
|
|
at=a(istage)*dt
|
|
bt=(b(istage)-a(istage))*dt
|
|
DO k=1,1 ! nz
|
|
DO j=1,1 ! ny
|
|
DO i=1,nx
|
|
r1(i,j,k)=r1(i,j,k)+at*f(i,j,k)
|
|
r2(i,j,k)=r1(i,j,k)+bt*f(i,j,k)
|
|
ENDDO
|
|
ENDDO
|
|
ENDDO
|
|
ELSE
|
|
bt=b(istage)*dt
|
|
DO k=1,1 ! nz
|
|
DO j=1,1 ! ny
|
|
DO i=1,nx
|
|
r1(i,j,k)=r1(i,j,k)+bt*f(i,j,k)
|
|
ENDDO
|
|
ENDDO
|
|
ENDDO
|
|
ENDIF
|
|
END SUBROUTINE substep
|
|
!------------------------------------------------------------------------
|
|
SUBROUTINE fns(r1,f)
|
|
REAL, INTENT(IN),DIMENSION(:,:,:) :: r1
|
|
REAL, INTENT(OUT),DIMENSION(:,:,:) :: f
|
|
REAL, DIMENSION(3,nx) :: ux,dux,d2ux
|
|
REAL, DIMENSION(1,nx) :: duy
|
|
INTEGER :: i,j,k
|
|
REAL :: wrate,Ly,Dy
|
|
|
|
DO k=1,1 !nz
|
|
! x-direction
|
|
DO j=1,1 !ny
|
|
DO i=1,nx
|
|
ux(1,i)=u(i,j,k)*r1(i,j,k) ! u*Y
|
|
ux(2,i)=u(i,j,k) ! u
|
|
ux(3,i)=r1(i,j,k) ! Y
|
|
ENDDO
|
|
|
|
dm = diff * ((1.0 + bc * (1.0 - r1(:,1,1))) ** 0.76)
|
|
|
|
CALL dfnonp(nx,hx,ux(1:3,:),dux(1:3,:),3,1)
|
|
CALL dfnonp(nx,hx,dm,duy(1,:),1,1)
|
|
CALL d2fnonp(nx,hx,ux(3:3,:),d2ux(1,:),1,1)
|
|
|
|
DO i=1,nx
|
|
wrate=pre*ux(3,i)*exp(-ac/(1.+bc*(1.-ux(3,i)))) !wrate
|
|
! IF ((1.-ux(3,i)).le.c_cut) THEN
|
|
! wrate=min_wr
|
|
! IF ((1.-ux(3,i)).gt.c_ref) wrate= &
|
|
! ((exp(((1.-ux(3,i))-c_cut)*prof_wr)-minf)/(1.-minf)*(refwr-min_wr))+min_wr
|
|
! ENDIF
|
|
IF ((1.-ux(3,i)).le.c_ref) THEN
|
|
wrate=min_wr
|
|
IF ((1.-ux(3,i)).gt.c_cut) wrate=((refwr-min_wr)*exp(prof_wr*(1.-ux(3,i)-c_ref))+ &
|
|
min_wr-refwr*exp(prof_wr*(c_cut-c_ref)))/(1.-exp(prof_wr*(c_cut-c_ref)))
|
|
ENDIF
|
|
|
|
! -0.5*( d(u*Y)/dx + u*dY/dx + Y*du/dx ) + D*d2Y/d2x
|
|
f(i,j,k)=-0.5*( dux(1,i) + ux(2,i)*dux(3,i) + &
|
|
ux(3,i)*dux(2,i) ) &
|
|
+ dm(i)*d2ux(1,i) + duy(1,i)*dux(3,i) - wrate
|
|
IF (i.eq.nx) THEN
|
|
Ly=ux(2,nx)*dux(3,nx) ! Ly = u*dYr/dx
|
|
f(nx,1,1)=-wrate
|
|
ENDIF
|
|
ENDDO
|
|
ENDDO
|
|
!! y-direction
|
|
! DO i=1,nx
|
|
! DO j=1,ny
|
|
! uy(1,j)=v(i,j,k)*r1(i,j,k) ! v*Y
|
|
! uy(2,j)=v(i,j,k) ! v
|
|
! uy(3,j)=r1(i,j,k) ! Y
|
|
! ENDDO
|
|
! CALL dfyz3(ny,hy,uy(1:3,:),duy)
|
|
! CALL d2fyz1(ny,hy,uy(3:3,:),d2uy)
|
|
! DO j=1,ny
|
|
! f(i,j,k)=f(i,j,k) + &
|
|
! -0.5*( duy(1,j) + uy(2,j)*duy(3,j) + uy(3,j)*duy(2,j) ) &
|
|
! + diff*d2uy(j)
|
|
! ENDDO
|
|
! ENDDO
|
|
ENDDO
|
|
! f(1,1,1)=1.0
|
|
! f(nx,1,1)=-0.5*( dux(1,nx) + ux(2,nx)*dux(3,nx) + ux(3,nx)*dux(2,nx) )
|
|
! f(nx,1,1)=0.0
|
|
!! z-direction
|
|
! DO j=1,ny
|
|
! DO i=1,nx
|
|
! DO k=1,nz
|
|
! uy(1,k)=w(i,j,k)*r1(i,j,k) ! w*Y
|
|
! uy(2,k)=w(i,j,k) ! w
|
|
! uy(3,k)=r1(i,j,k) ! Y
|
|
! ENDDO
|
|
! CALL dfyz3(ny,hy,uy(1:3,:),duy)
|
|
! CALL d2fyz1(ny,hy,uy(3:3,:),d2uy)
|
|
! DO k=1,nz
|
|
! f(i,j,k)=f(i,j,k) + &
|
|
! -0.5*( duy(1,k) + uy(2,k)*duy(3,k) + uy(3,k)*duy(2,k) ) &
|
|
! + diff*d2uy(k)
|
|
! ENDDO
|
|
! ENDDO
|
|
! ENDDO
|
|
|
|
! DO j=1,ny
|
|
! DO k=1,nz
|
|
! f(1,j,k)=0.0
|
|
! f(nx,j,k)=-0.5*(dux(1,nx)+ux(2,nx)*dux(3,nx)+ux(3,nx)*&
|
|
! dux(2,nx))
|
|
! ENDDO
|
|
! ENDDO
|
|
|
|
! Boundary conditionS
|
|
f(1,1,1)=0.
|
|
|
|
Dy=Ly
|
|
f(nx,1,1)=f(nx,1,1)-Dy
|
|
|
|
END SUBROUTINE fns
|
|
!------------------------------------------------------------------------
|
|
END MODULE ysolve
|