research/incomp-flame-1d
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

solve T

Browse source
This commit is contained in:
ignis 2019-03-28 15:58:09 +09:00
parent 358d1fea3c
commit 9f0647865c
2 changed files with 45 additions and 54 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

3
code/sample/4pi-IC1/itape
View file

@ -16,6 +16,5 @@
c_ref 0.01 !0.003 !
min_wr 0. ! 5.0e-14
prof_wr 1.0
rvis 1.0
cstar 1.0
lewis 1.0

96
code/ysolve.f90
View file

@ -4,12 +4,13 @@
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,sigw
REAL :: lewis
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 :: u
REAL, DIMENSION(:,:), ALLOCATABLE :: y1,y2,yf
REAL, DIMENSION(:), ALLOCATABLE :: uxt,duxt
REAL, DIMENSION(:), ALLOCATABLE :: dm
@ -83,12 +84,12 @@
ENDDO
CALL dfnonp(nx,hx,uxt,duxt,1,1)
delf=1./MAXVAL(ABS(duxt))
oldu=u(1,1)
oldu=u(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),' // L_f:',delf
' // Old_U:',oldu,' // New_U:',u(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(*,*)
@ -113,9 +114,9 @@
END SUBROUTINE solve
SUBROUTINE write_sd
REAL :: c,yr,wrate,dely,sdr,sdd,sd,uu,onelw,dd
REAL :: c,yr,theta,wrate,dely,sdr,sdd,sd,uu,onelw,dd
INTEGER :: i,j,k,nd
REAL, DIMENSION(1,nx) :: ux,dux,d2ux
REAL, DIMENSION(2,nx) :: ux,dux,d2ux
REAL, DIMENSION(10,nx) :: sav
sav=0.
@ -130,20 +131,21 @@
DO i=1,nx
ux(1,i)=y1(i,1) ! Yr
ux(2,i)=y1(i,2) ! T
IF (ux(1,i).gt.1.) ux(1,i)=1.
ENDDO
nd=1
CALL dfnonp(nx,hx,ux(1,:),dux(1,:),nd,1)
nd=1
CALL d2fnonp(nx,hx,ux(1,:),d2ux(1,:),nd,1)
nd=2
CALL dfnonp(nx,hx,ux(:,:),dux(:,:),nd,1)
CALL d2fnonp(nx,hx,ux(:,:),d2ux(:,:),nd,1)
DO i=1,nx
yr=ux(1,i)
theta=ux(2,i)
c=1.-yr
IF (c.lt.0.) c=0.
wrate=pre*yr*exp(-ac/(1.+bc*(1.-yr))) !wrate
wrate=pre*yr*exp(-ac/(1.+bc*(theta))) !wrate
! IF (c.le.c_cut) THEN
! wrate=min_wr
@ -151,9 +153,9 @@
! ((exp((c-c_cut)*prof_wr)-minf)/(1.-minf)*(refwr-min_wr))+min_wr
! ENDIF
IF (c.le.c_ref) THEN
IF (theta.le.c_ref) THEN
wrate=min_wr
IF (c.gt.c_cut) wrate=((refwr-min_wr)*exp(prof_wr*(c-c_ref))+ &
IF (theta.gt.c_cut) wrate=((refwr-min_wr)*exp(prof_wr*(theta-c_ref))+ &
min_wr-refwr*exp(prof_wr*(c_cut-c_ref)))/(1.-exp(prof_wr*(c_cut-c_ref)))
ENDIF
@ -164,7 +166,7 @@
sdr=0.; sdd=0.
ENDIF
sd=sdr+sdd
uu=u(1,1)
uu=u(1)
onelw=(-dux(1,i))/c
dd=-dm(i)*d2ux(1,i)
if (c.eq.0.) onelw=0.
@ -197,32 +199,33 @@
OPEN (305,FILE='sfield.dat')
DO i=1,nx
WRITE (305,'(e30.20)') y1(i,1)
WRITE (305,'(2e30.20)') y1(i,1), y1(i,2)
ENDDO
CLOSE (305)
END SUBROUTINE save_final_field
SUBROUTINE write_pre
REAL :: yr,c,dy,maxdy=0.,del_f
REAL :: theta,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)
theta=y1(i,2)
c=1.-yr
ux(1,i)=yr
wrate=pre*yr*exp(-ac/(1.+bc*(1.-yr))) !wrate
wrate=pre*yr*exp(-ac/(1.+bc*theta)) !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
IF (theta.le.c_ref) THEN
wrate=min_wr
IF (c.gt.c_cut) wrate=((refwr-min_wr)*exp(prof_wr*(c-c_ref))+ &
IF (theta.gt.c_cut) wrate=((refwr-min_wr)*exp(prof_wr*(theta-c_ref))+ &
min_wr-refwr*exp(prof_wr*(c_cut-c_ref)))/(1.-exp(prof_wr*(c_cut-c_ref)))
ENDIF
@ -239,7 +242,7 @@
ENDDO
del_f=1./maxdy
WRITE(*,'(a13,e14.8,a25,e14.8)') ' Grid size : ',hx,' / Laminar flame speed : ',u(1,1)
WRITE(*,'(a13,e14.8,a25,e14.8)') ' Grid size : ',hx,' / Laminar flame speed : ',u(1)
WRITE(*,'(a19,e14.8,a3,f9.5,a20)') &
' Flame thickness : ',del_f,' / ',del_f/hx,' grids in the flame.'
WRITE(*,*)
@ -330,16 +333,8 @@
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
READ(itape,*) cdum,sigw
WRITE(otape,*) cdum,sigw
DO i = 0,100
WRITE(*,*) diffusivity_T076(i / 100.)
END DO
READ(itape,*) cdum,lewis
WRITE(otape,*) cdum,lewis
IF (.not.read_stdin) THEN
CLOSE(itape)
@ -363,7 +358,7 @@
refwr=pre*1.*exp(-ac/(1.+bc*c_ref))
l_0=l_0*pi
ALLOCATE( u(nx,1),STAT=ierr) ; u=0.
ALLOCATE( u(nx),STAT=ierr) ; u=0.
ALLOCATE(y1(nx,2),STAT=ierr) ; y1=0.
ALLOCATE(y2(nx,2),STAT=ierr) ; y2=0.
ALLOCATE(yf(nx,2),STAT=ierr) ; yf=0.
@ -405,7 +400,7 @@
REAL, INTENT(OUT),DIMENSION(:,:) :: f
REAL, DIMENSION(2,nx) :: ux, dux, d2ux
INTEGER :: i,j,k
REAL :: wrate,Ly,Dy
REAL :: wrate,Ly,Dy,Lt,Dt
! x-direction
DO i=1,nx
@ -418,20 +413,26 @@
CALL d2fnonp(nx,hx,ux(:,:),d2ux(:,:),2,1)
DO i=1,nx
wrate=pre*ux(1,i)*exp(-ac/(1.+bc*(1.-ux(1,i)))) !wrate
wrate=pre*ux(1,i)*exp(-ac/(1.+bc*(ux(2,i)))) !wrate
IF ((1.-ux(1,i)).le.c_ref) THEN
IF ((ux(2,i)).le.c_ref) THEN
wrate=min_wr
IF ((1.-ux(1,i)).gt.c_cut) wrate=((refwr-min_wr)*exp(prof_wr*(1.-ux(1,i)-c_ref))+ &
IF ((ux(2,i)).gt.c_cut) wrate=((refwr-min_wr)*exp(prof_wr*(ux(2,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,1) = - ( u(i,1)*dux(1,i) ) + diff * d2ux(1,i) - wrate
f(i,1) = - ( u(i)*dux(1,i) ) + diff * d2ux(1,i) - wrate
! - u*dT/dx + D*d2T/d2x
f(i,2) = - ( u(i)*dux(2,i) ) + (lewis * diff) * d2ux(2,i) + wrate
IF (i.eq.nx) THEN
Ly=u(nx,1)*dux(1,nx) ! Ly = u*dYr/dx
Ly=u(nx)*dux(1,nx) ! Ly = u*dYr/dx
f(nx,1)=-wrate
Lt=u(nx)*dux(2,nx) ! Ly = u*dYr/dx
f(nx,2)=wrate
ENDIF
ENDDO
@ -441,21 +442,12 @@
Dy=Ly
f(nx,1)=f(nx,1)-Dy
f(1,2)=0.
Dt=Lt
f(nx,2)=f(nx,2)-Dt
END SUBROUTINE fns
!------------------------------------------------------------------------
REAL FUNCTION sigmoid(x)
REAL, INTENT(IN) :: x
sigmoid = 1. / (1 + exp(-x))
END FUNCTION sigmoid
!------------------------------------------------------------------------
REAL FUNCTION diffusivity_sigmoid(c)
REAL, INTENT(IN) :: c
diffusivity_sigmoid = 1. + (rvis - 1.) * sigmoid(sigw*(c-cstar))
END FUNCTION diffusivity_sigmoid
!------------------------------------------------------------------------
REAL FUNCTION diffusivity_T076(c)
REAL, INTENT(IN) :: c
diffusivity_T076 = ((1.0 + bc * c) ** 0.76)
END FUNCTION diffusivity_T076
!------------------------------------------------------------------------
END MODULE ysolve