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

added solution variable T

Browse source
This commit is contained in:
ignis 2019-03-27 13:35:36 +09:00
parent a61e753b27
commit c9c16c3f42
1 changed files with 43 additions and 62 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

105
code/ysolve.f90
View file

@ -9,7 +9,7 @@
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,y1,y2,yf
REAL, DIMENSION(:), ALLOCATABLE :: uxt,duxt
REAL, DIMENSION(:), ALLOCATABLE :: dm
@ -78,17 +78,17 @@
pflold=pflame
pflame=0.
DO i=1,nx
pflame=pflame+y1(i,1,1)*hx
uxt(i)=y1(i,1,1)
pflame=pflame+y1(i,1)*hx
uxt(i)=y1(i,1)
ENDDO
CALL dfnonp(nx,hx,uxt,duxt,1,1)
delf=1./MAXVAL(ABS(duxt))
oldu=u(1,1,1)
oldu=u(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
' // Old_U:',oldu,' // New_U:',u(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(*,*)
@ -129,7 +129,7 @@
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
ux(1,i)=y1(i,1) ! Yr
IF (ux(1,i).gt.1.) ux(1,i)=1.
ENDDO
@ -164,7 +164,7 @@
sdr=0.; sdd=0.
ENDIF
sd=sdr+sdd
uu=u(1,1,1)
uu=u(1,1)
onelw=(-dux(1,i))/c
dd=-dm(i)*d2ux(1,i)
if (c.eq.0.) onelw=0.
@ -191,13 +191,13 @@
OPEN (305,FILE='sfield.bin',form='unformatted',status='unknown')
DO i=1,nx
WRITE (305) y1(i,1,1)
WRITE (305) y1(i,1)
ENDDO
CLOSE (305)
OPEN (305,FILE='sfield.dat')
DO i=1,nx
WRITE (305,'(e30.20)') y1(i,1,1)
WRITE (305,'(e30.20)') y1(i,1)
ENDDO
CLOSE (305)
@ -210,7 +210,7 @@
REAL, DIMENSION(1,nx) :: ux, dux
DO i=1,nx
yr=y1(i,1,1)
yr=y1(i,1)
c=1.-yr
ux(1,i)=yr
@ -239,7 +239,7 @@
ENDDO
del_f=1./maxdy
WRITE(*,'(a13,e14.8,a25,e14.8)') ' Grid size : ',hx,' / Laminar flame speed : ',u(1,1,1)
WRITE(*,'(a13,e14.8,a25,e14.8)') ' Grid size : ',hx,' / Laminar flame speed : ',u(1,1)
WRITE(*,'(a19,e14.8,a3,f9.5,a20)') &
' Flame thickness : ',del_f,' / ',del_f/hx,' grids in the flame.'
WRITE(*,*)
@ -247,7 +247,7 @@
SUBROUTINE SET_BC
y1(1,1,1)=1.-ctmp
y1(1,1)=1.-ctmp
END SUBROUTINE SET_BC
SUBROUTINE SET_IC
@ -268,12 +268,14 @@
ENDIF
! y1(i,1,1)=(1.-xi)*1.
y1(i,1,1)=(1.-xi) ! reactant mass fraction
y1(i,1)=(1.-xi) ! reactant mass fraction
ENDDO
y1(:,2) = 1.0 - y1(:,1)
pflame=0.
DO i=1,nx
pflame=pflame+y1(i,1,1)*hx
pflame=pflame+y1(i,1)*hx
ENDDO
END SUBROUTINE SET_IC
@ -361,10 +363,10 @@
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( u(nx,1),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.
ALLOCATE(uxt(nx),STAT=ierr) ; uxt=0.
ALLOCATE(duxt(nx),STAT=ierr) ; duxt=0.
ALLOCATE(dm(nx),STAT=ierr) ; dm=diff
@ -380,8 +382,8 @@
END SUBROUTINE RK4
!------------------------------------------------------------------------
SUBROUTINE substep(ri,r1,r2,f)
REAL, INTENT(INOUT),DIMENSION(:,:,:) :: ri,r1,r2
REAL, INTENT(OUT),DIMENSION(:,:,:) :: f
REAL, INTENT(INOUT),DIMENSION(:,:) :: ri,r1,r2
REAL, INTENT(OUT),DIMENSION(:,:) :: f
INTEGER :: i,j,k
REAL :: at,bt
@ -390,75 +392,54 @@
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
r1=r1+at*f
r2=r1+bt*f
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
r1=r1+bt*f
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, INTENT(IN),DIMENSION(:,:) :: r1
REAL, INTENT(OUT),DIMENSION(:,:) :: f
REAL, DIMENSION(2,nx) :: ux, dux, d2ux
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
ux(1,i)=r1(i,1) ! Y
ux(2,i)=r1(i,2) ! T
ENDDO
CALL dfnonp(nx,hx,ux(1:3,:),dux(1:3,:),3,1)
CALL d2fnonp(nx,hx,ux(3:3,:),d2ux(1,:),1,1)
CALL dfnonp(nx,hx,ux(:,:),dux(:,:),2,1)
CALL d2fnonp(nx,hx,ux(:,:),d2ux(:,:),2,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=pre*ux(1,i)*exp(-ac/(1.+bc*(1.-ux(1,i)))) !wrate
IF ((1.-ux(1,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))+ &
IF ((1.-ux(1,i)).gt.c_cut) wrate=((refwr-min_wr)*exp(prof_wr*(1.-ux(1,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) - wrate
f(i,1) = - ( u(i,1)*dux(1,i) ) + diff * d2ux(1,i) - wrate
IF (i.eq.nx) THEN
Ly=ux(2,nx)*dux(3,nx) ! Ly = u*dYr/dx
f(nx,1,1)=-wrate
Ly=u(nx,1)*dux(1,nx) ! Ly = u*dYr/dx
f(nx,1)=-wrate
ENDIF
ENDDO
ENDDO
ENDDO
! Boundary conditionS
f(1,1,1)=0.
f(1,1)=0.
Dy=Ly
f(nx,1,1)=f(nx,1,1)-Dy
f(nx,1)=f(nx,1)-Dy
END SUBROUTINE fns
!------------------------------------------------------------------------