diff --git a/code/ysolve.f90 b/code/ysolve.f90 index 0f5d671..6450166 100644 --- a/code/ysolve.f90 +++ b/code/ysolve.f90 @@ -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 !------------------------------------------------------------------------