research/dns-hit3d-fdm
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

governing equation changed to consider variable diffusivity

Browse source
This commit is contained in:
ignis 2019-03-24 10:48:55 +09:00
parent 0776f3ef6e
commit f38b995827
1 changed files with 24 additions and 19 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

43
m_fdm_calc.f90
View file

@ -9,6 +9,7 @@ module m_fdm_calc
!variables
real*8, dimension(:,:,:), allocatable :: u_,v_,w_
real*8, dimension(:,:,:), allocatable :: dm
real*8, dimension(:,:,:,:), allocatable :: y1,y2,yf
real*8, dimension(:,:), allocatable :: fzzl, fzzu
@ -233,6 +234,10 @@ module m_fdm_calc
v_=0.0
w_=0.0
allocate(dm(nx,ny,nz))
dm = diff
! DQ initializing
fdmsavecount=1 !FDM save count
sum_wrate=0.
@ -529,7 +534,7 @@ module m_fdm_calc
DO j=1,yy
DO i=1,xx
DO k=1,zz
uz(1,k)=r1_(i,j,k,1) ! 1:rho
uz(1,k)=r1_(i,j,k,1)*dm(i,j,k) ! 1:rho*D
uz(2,k)=r1_(i,j,k,2)/r1_(i,j,k,1) ! 2:Y
uz(3,k)=uz(1,k)*ww_(i,j,k) ! 3:rho*w
uz(4,k)=uz(3,k)*uz(2,k) ! 4:rho*w*Y
@ -545,8 +550,8 @@ module m_fdm_calc
! -( d(rho*w*Yr)/dz ) + d(rho*D* d(Yr)/dz)/dz
! = -( d(rho*w*Yr)/dz )
! + D* (rho* d2(Yr)/dz2 + d(rho)/dz * d(Yr)/dz )
f_(i,j,k,2) = f_(i,j,k,2) - duz(4,k) + diff*(uz(1,k)*d2uz(k)+duz(1,k)*duz(2,k)) ! species conserv.
! + (rho*D * d2(Yr)/dz2 + d(rho*D)/dz * d(Yr)/dz )
f_(i,j,k,2) = f_(i,j,k,2) - duz(4,k) + (uz(1,k)*d2uz(k)+duz(1,k)*duz(2,k)) ! species conserv.
ENDDO
ENDDO
ENDDO
@ -563,7 +568,7 @@ module m_fdm_calc
! -( d(rho*w*Yr)/dz ) + d(rho*D* d(Yr)/dz)/dz
! = -( d(rho*w*Yr)/dz )
! + D* (rho* d2(Yr)/dz2 + d(rho)/dz * d(Yr)/dz )
! + (rho*D * d2(Yr)/dz2 + d(rho*D)/dz * d(Yr)/dz )
fbuf1(:,:,:) = r1_(:,:,:,2)*ww_(:,:,:) ! rho*w*Y
CALL pdfp (fbuf1, fzzl, fzzu, fbuf2, hy, xx*yy, zz, yy, 3)
@ -574,13 +579,13 @@ module m_fdm_calc
CALL pdfp (fbuf1, fzzl, fzzu, fbuf2, hy, xx*yy, zz, yy, 3)
CALL pd2fp(fbuf1, fzzl, fzzu, fbuf3, hy, xx*yy, zz, yy, 3)
fbuf1(:,:,:) = r1_(:,:,:,1) ! rho
fbuf1(:,:,:) = r1_(:,:,:,1) ! rho*D
CALL pdfp (fbuf1, fzzl, fzzu, fbuf4, hy, xx*yy, zz, yy, 3)
fbuf2 = fbuf2 * fbuf4
fbuf1 = fbuf1 * fbuf3 + fbuf2
fbuf2 = fbuf2 * fbuf4 ! dY/dz * d(rho*D)/dz
fbuf1 = fbuf1 * fbuf3 + fbuf2 ! rho*D * d2(rho*D)/dz2 + ...
f_(:,:,:,2) = f_(:,:,:,2) + diff*fbuf1(:,:,:) ! species conserv.
f_(:,:,:,2) = f_(:,:,:,2) + fbuf1(:,:,:) ! species conserv.
endif
@ -595,7 +600,7 @@ module m_fdm_calc
! -( d(rho*u*Yr)/dx ) + d(rho*D* d(Yr)/dx)/dx
! = -( d(rho*u*Yr)/dx )
! + D* (rho* d2(Yr)/dx2 + d(rho)/dx * d(Yr)/dx )
! + (rho*D * d2(Yr)/dx2 + d(rho*D)/dx * d(Yr)/dx )
CALL tp2mul (yxbuf1, r1_(:,:,k,2), uu_(:,:,k), yy, xx) ! rho*u*Y
CALL dfnonp(xx,hx,yxbuf1,yxbuf2,yy,1) ! d/dx(rho*u*Y)
@ -605,11 +610,11 @@ module m_fdm_calc
CALL dfnonp(xx,hx,yxbuf1,yxbuf2,yy,1) ! d/dx(Y)
CALL d2fnonp(xx,hx,yxbuf1,yxbuf3,yy,1) ! d2/dx2(Y)
CALL tp2 (yxbuf1, r1_(:,:,k,1), yy, xx) ! rho
CALL dfnonp(xx,hx,yxbuf1,yxbuf4,yy,1) ! d/dx(rho)
CALL tp2mul (yxbuf1, r1_(:,:,k,1), dm(:,:,k), yy, xx) ! rho*D
CALL dfnonp(xx,hx,yxbuf1,yxbuf4,yy,1) ! d/dx(rho*D)
yxbuf2 = yxbuf2 * yxbuf4 ! d/dx(Y) * d/dx(rho)
yxbuf1 = -diff*(yxbuf2 + yxbuf1 * yxbuf3) ! -D( ... + rho * d2/dx2(Y))
yxbuf2 = yxbuf2 * yxbuf4 ! d/dx(Y) * d/dx(rho*D)
yxbuf1 = -(yxbuf2 + yxbuf1 * yxbuf3) ! -( ... + rho*D * d2/dx2(Y))
CALL tp2subasgn (f_(:,:,k,2), yxbuf1, xx, yy) ! species conservation
@ -628,7 +633,7 @@ module m_fdm_calc
! -( d(rho*v*Yr)/dy ) + d(rho*D* d(Yr)/dy)/dy
! = -( d(rho*v*Yr)/dy )
! + D* (rho* d2(Yr)/dyy2 + d(rho)/dy * d(Yr)/dy )
! + (rho*D* d2(Yr)/dyy2 + d(rho*D)/dy * d(Yr)/dy )
xybuf1(:,:) = r1_(:,:,k,2)*vv_(:,:,k) ! rho*v*Y
CALL dfp(yy,hy,xybuf1,xybuf2,xx,2) ! d/dy(rho*v*Y)
@ -639,13 +644,13 @@ module m_fdm_calc
CALL dfp(yy,hy,xybuf1,xybuf2,xx,2) ! d/dy(Y)
CALL d2fp(yy,hy,xybuf1,xybuf3,xx,2) ! d2/dy2(Y)
xybuf1(:,:) = r1_(:,:,k,1) ! rho
CALL dfp(yy,hy,xybuf1,xybuf4,xx,2) ! d/dy(rho)
xybuf1(:,:) = r1_(:,:,k,1) * dm(:,:,k) ! rho*D
CALL dfp(yy,hy,xybuf1,xybuf4,xx,2) ! d/dy(rho*D)
xybuf2 = xybuf2 * xybuf4 ! d/dy(Y) * d/dy(rho)
xybuf1 = xybuf2 + xybuf1 * xybuf3 ! ... + rho * d2/dy2(Y)
xybuf2 = xybuf2 * xybuf4 ! d/dy(Y) * d/dy(rho*D)
xybuf1 = xybuf2 + xybuf1 * xybuf3 ! ... + rho*D * d2/dy2(Y)
f_(:,:,k,2) = f_(:,:,k,2) + diff*xybuf1(:,:) ! species conserv.
f_(:,:,k,2) = f_(:,:,k,2) + xybuf1(:,:) ! species conserv.
ENDDO