From 9f0647865cf5206078a514c866cbdba94c619acd Mon Sep 17 00:00:00 2001 From: ignis Date: Thu, 28 Mar 2019 15:58:09 +0900 Subject: [PATCH] solve T --- code/sample/4pi-IC1/itape | 3 +- code/ysolve.f90 | 96 ++++++++++++++++++--------------------- 2 files changed, 45 insertions(+), 54 deletions(-) diff --git a/code/sample/4pi-IC1/itape b/code/sample/4pi-IC1/itape index 993a20d..2c08701 100644 --- a/code/sample/4pi-IC1/itape +++ b/code/sample/4pi-IC1/itape @@ -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 diff --git a/code/ysolve.f90 b/code/ysolve.f90 index 6450166..5cf3da0 100644 --- a/code/ysolve.f90 +++ b/code/ysolve.f90 @@ -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