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

minor mods ysolve

Browse source
This commit is contained in:
ignis 2019-08-31 07:33:04 +09:00
parent cf620e1668
commit 1725f1eaf1
1 changed files with 35 additions and 25 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

60
code/ysolve.f90
View file

@ -65,14 +65,15 @@
INTEGER :: i,j,k,savenum
REAL :: pflame,pflold,delf=0.
CALL ludcmp(nx,5,5,1,0,0)
CALL ludcmp(nx,100,100,1,0,0)
CALL SET_IC
t_now=0.; t_uf=0.
t_now=0.
t_uf=0.
DO
CALL SET_BC
DO
CALL SET_BC
CALL RK4(fns)
IF(t_now.ge.tf) EXIT
IF(t_uf.ge.dt_uf) THEN
@ -249,7 +250,6 @@
END SUBROUTINE write_pre
SUBROUTINE SET_BC
y1(1,1)=1.-ctmp
END SUBROUTINE SET_BC
@ -257,10 +257,16 @@
INTEGER :: i, ifl, si
REAL :: xi
! initial flame thickness (0.2 pi n grids)
ifl=INT(2.0*pi*0.1/hx)
u=u0; si=INT(nx*(1.-tar_lo))
! initialize velocity field
u=u0;
! initial flame center (grid index)
si=INT(nx*(1.-tar_lo))
! initialize Yr field
DO i=1,nx
IF(i< nx-(si+ifl/2)) THEN
xi=0.+ctmp
@ -274,8 +280,10 @@
y1(i,1)=(1.-xi) ! reactant mass fraction
ENDDO
! initialize temperature field
y1(:,2) = 1.0 - y1(:,1)
! calculate flame position
pflame=0.
DO i=1,nx
pflame=pflame+y1(i,1)*hx
@ -320,28 +328,29 @@
END SUBROUTINE rhs
END INTERFACE
REAL, DIMENSION(5), PARAMETER :: a=(/ 970286171893./4311952581923., &
6584761158862./12103376702013., &
2251764453980./15575788980749., &
26877169314380./34165994151039., &
0. /), &
b=(/ 1153189308089./22510343858157., &
1772645290293./4653164025191., &
-1672844663538./4480602732383., &
2114624349019./3568978502595., &
5198255086312./14908931495163. /)
INTEGER :: istage
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)
CALL substep(1,y1,y1,y2,yf)
CALL substep(2,y1,y2,y1,yf)
CALL substep(3,y2,y1,y2,yf)
CALL substep(4,y1,y2,y1,yf)
CALL substep(5,y2,y1,y2,yf)
CONTAINS
SUBROUTINE substep(ri,r1,r2,f)
SUBROUTINE substep(istage, ri,r1,r2,f)
REAL, DIMENSION(5), PARAMETER :: a=(/ 970286171893.d0/4311952581923., &
6584761158862.d0/12103376702013., &
2251764453980.d0/15575788980749., &
26877169314380.d0/34165994151039., &
0.d0 /), &
b=(/ 1153189308089.d0/22510343858157., &
1772645290293.d0/4653164025191., &
-1672844663538.d0/4480602732383., &
2114624349019.d0/3568978502595., &
5198255086312.d0/14908931495163. /)
INTEGER :: istage
REAL, INTENT(INOUT),DIMENSION(:,:) :: ri,r1,r2
REAL, INTENT(OUT),DIMENSION(:,:) :: f
@ -349,6 +358,7 @@
REAL :: at,bt
CALL rhs(ri,f)
IF(istage<5) THEN
at=a(istage)*dt
bt=(b(istage)-a(istage))*dt