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

working version

Browse source
This commit is contained in:
ignis 2019-09-05 17:41:29 +09:00
parent d5bb7d388a
commit 4e04053f17
4 changed files with 58 additions and 19 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

5
code/m_chemistry.f90
View file

@ -63,7 +63,8 @@ real function rate1_2step (ya, yx, theta)
if(theta.lt.0.) t_reduce=0.
if(theta.gt.1.) t_reduce=1.
rate1_2step = 0.
rate1_2step = lambda1 * y1 * y2 * &
exp (-(beta1*(1. - t_reduce))/(1. - hrp*(1. - t_reduce)))
end function rate1_2step
@ -84,7 +85,7 @@ real function rate2_2step (yx, theta)
if(theta.lt.0.) t_reduce=0.
if(theta.gt.1.) t_reduce=1.
rate2_2step = 0.
rate2_2step = lambda2 * yx * yx
end function rate2_2step

25
code/m_parameters.f90
View file

@ -16,6 +16,8 @@ module m_parameters
real :: pre,ac,bc,c_cut,c_ref
real :: lambda1, lambda2, beta1, hrp
integer :: syp,eyp,twod
real :: SL_u
@ -28,7 +30,7 @@ module m_parameters
REAL :: hx,dt,vis,sc,diff,tf,t_now,t_uf,dt_uf
REAL :: lewis
REAL :: lewis, le_a, le_x
REAL :: minf,tar_lo,u0,ctmp,lo_flm=0.
@ -44,6 +46,15 @@ module m_parameters
contains
SUBROUTINE SET_CHEMISTRY
lambda1 = 100000.
lambda2 = lambda1 * 63.6
beta1 = 9.1
hrp = 0.8136765
END SUBROUTINE SET_CHEMISTRY
SUBROUTINE READ_INTRO
CHARACTER(LEN=8) :: cdum
INTEGER :: itape=300,otape=301
@ -96,6 +107,18 @@ contains
WRITE(otape,*) cdum,prof_wr
READ(itape,*) cdum,lewis
WRITE(otape,*) cdum,lewis
READ(itape,*) cdum,lambda1
WRITE(otape,*) cdum,lambda1
READ(itape,*) cdum,lambda2
WRITE(otape,*) cdum,lambda2
READ(itape,*) cdum,beta1
WRITE(otape,*) cdum,beta1
READ(itape,*) cdum,hrp
WRITE(otape,*) cdum,hrp
READ(itape,*) cdum,le_a
WRITE(otape,*) cdum,le_a
READ(itape,*) cdum,le_x
WRITE(otape,*) cdum,le_x
IF (.not.read_stdin) THEN
CLOSE(itape)

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

@ -9,7 +9,7 @@
ac 26.7 !22.250 !22.25 !22.56 !22.3 !22.39
bc 3. !7.000 !7.0 !7.0 !7.0 !7.0
u0 0.21
tf 80.
tf 200.
dt_uf 0.4
ctmp 0. !1.0e-14
c_cut 0.001 !0.012 ! c < c_cut -> wrate = 0.
@ -17,4 +17,9 @@
min_wr 0. ! 5.0e-14
prof_wr 1.0
lewis 1.0
lambda1 100000.
lambda2 63.6E+5
beta1 9.1
hrp 0.8136765
le_a 1.
le_x 0.2

38
code/ysolve.f90
View file

@ -223,7 +223,7 @@
OPEN (305,FILE='sfield.dat')
DO i=1,nx
WRITE (305,'(e30.20)') y1(i,2)
WRITE (305,'(4e30.20)') i*hx, y1(i,1), y1(i,2), y1(i,3)
ENDDO
CLOSE (305)
@ -273,14 +273,15 @@
END SUBROUTINE write_pre
SUBROUTINE SET_BC
y1(1,1)=1.-ctmp
y1(1,2)=1.
y1(1,1)=1.
y1(1,2)=0.
y1(1,3)=0.
END SUBROUTINE SET_BC
SUBROUTINE SET_IC
INTEGER :: i, ifl, si
REAL :: xi
REAL :: max_ysum
! initial flame thickness (0.2 pi n grids)
ifl=INT(2.0*pi*0.1/hx)
@ -303,13 +304,20 @@
! y1(i,1,1)=(1.-xi)*1.
y1(i,1)=(1.-xi) ! reactant mass fraction
y1(i,2) = (1./2.) * (lambda1/lambda2) * y1(i,1) * &
exp (-(beta1*(1. - xi))/(1. - hrp*(1. - xi)))
y1(i,3) = xi
ENDDO
! initialize intermediate species field
y1(:,2) = y1(:,1)
! initialize temperature field
y1(:,3) = 1.0 - y1(:,1)
! initialize species field
! max_ysum = maxval(y1(:,1) + y1(:,2))
! y1(:,1) = y1(:,1) / maxval(y1(:,1) + y1(:,2))
! y1(:,2) = y1(:,2) / maxval(y1(:,1) + y1(;,2))
! y1(:,1) = y1(:,1) / (y1(1,1) + y1(1,2))
! y1(:,2) = y1(:,2) / (y1(1,1) + y1(1,2))
! calculate flame position
pflame=0.
@ -423,21 +431,23 @@
DO i=1,nx
wrate=rate_1step(ux(1,i), ux(3,i))
wrate1=rate1_2step(ux(1,i), ux(2,i), ux(3,i))
wrate2=rate2_2step(ux(2,i), ux(3,i))
! - u*dY/dx + D*d2Y/d2x
f(i,1) = - ( u(i)*dux(1,i) ) + diff * d2ux(1,i) - wrate
f(i,1) = - ( u(i)*dux(1,i) ) + (diff/le_a) * d2ux(1,i) - wrate1
! - u*dY/dx + D*d2Y/d2x
f(i,2) = - ( u(i)*dux(2,i) ) + diff * d2ux(2,i) - wrate
f(i,2) = - ( u(i)*dux(2,i) ) + (diff/le_x) * d2ux(2,i) + wrate1 - 2.*wrate2
! - u*dT/dx + D*d2T/d2x
f(i,3) = - ( u(i)*dux(3,i) ) + (lewis * diff) * d2ux(3,i) + wrate
f(i,3) = - ( u(i)*dux(3,i) ) + (diff) * d2ux(3,i) + 2.*wrate2
! Boundary conditions
IF (i.eq.nx) THEN
f(nx,1) = -wrate - u(nx)*dux(1,nx)
f(nx,2) = -wrate - u(nx)*dux(2,nx)
f(nx,3) = wrate - u(nx)*dux(3,nx)
f(nx,1) = -wrate1 - u(nx)*dux(1,nx)
f(nx,2) = wrate1 - 2.*wrate2 - u(nx)*dux(2,nx)
f(nx,3) = 2.*wrate2 - u(nx)*dux(3,nx)
ENDIF
ENDDO