From d3fe42f947a635c0f027a0c95ca5c92fc982f587 Mon Sep 17 00:00:00 2001 From: Paul Date: Tue, 16 Jan 2018 18:06:37 +0000 Subject: [PATCH] Adding Froude numbers i.e. gravity --- convdiff.f90 | 180 +++++++++++++++++++++++++++++++++++++++++++++++ incompact3d.prm | 3 + module_param.f90 | 2 +- parameters.f90 | 27 ++++--- 4 files changed, 200 insertions(+), 12 deletions(-) diff --git a/convdiff.f90 b/convdiff.f90 index 9cdb8bd..819a906 100644 --- a/convdiff.f90 +++ b/convdiff.f90 @@ -69,6 +69,8 @@ subroutine convdiff(ux1,uy1,uz1,rho1,ta1,tb1,tc1,td1,te1,tf1,tg1,th1,ti1,di1,& integer :: code real(mytype) :: x,y,z + real(mytype) :: invfr + real(mytype), parameter :: ONETHIRD = 1._mytype / 3._mytype nvect1=xsize(1)*xsize(2)*xsize(3) @@ -429,6 +431,19 @@ subroutine convdiff(ux1,uy1,uz1,rho1,ta1,tb1,tc1,td1,te1,tf1,tg1,th1,ti1,di1,& tc1(:,:,:) = tc1(:,:,:) + xnu * tf1(:,:,:) endif + !! Gravity + if (frx.gt.0._mytype) then + invfr = 1._mytype / frx + ta1(:,:,:) = ta1(:,:,:) - rho1(:,:,:) * invfr + endif + if (fry.gt.0._mytype) then + invfr = 1._mytype / fry + tb1(:,:,:) = tb1(:,:,:) - rho1(:,:,:) * invfr + endif + if (frz.gt.0._mytype) then + invfr = 1._mytype / frz + tc1(:,:,:) = tc1(:,:,:) - rho1(:,:,:) * invfr + endif ! !! MMS Source term ! call momentum_source_mms(ta1,tb1,tc1) @@ -1119,3 +1134,168 @@ SUBROUTINE momentum_source_mms(mmsx1, mmsy1, mmsz1) ENDDO ! End loop over k ENDSUBROUTINE momentum_source_mms + +!!-------------------------------------------------------------------- +!! SUBROUTINE: entrainment_bcz +!! DESCRIPTION: Implements entrainment boundary conditions on z +!! stencil, based on work of Ioannou Vasilis. +!!-------------------------------------------------------------------- +SUBROUTINE entrainment_bcz(ux3, uy3, uz3, clx3, cly3, clz3) + + USE decomp_2d + USE variables + USE param + + IMPLICIT NONE + + REAL(mytype), DIMENSION(zsize(1), zsize(2), zsize(3)), INTENT(IN) :: ux3, uy3, uz3 + REAL(mytype), DIMENSION(zsize(1), zsize(2), zsize(3)) :: clx3, cly3, clz3 + + REAL(mytype) :: uu1, uv1, uw1 + REAL(mytype) :: y, z, yc, zc, ya + REAL(mytype) :: x1, x2, y1, y2, r1, r2 + INTEGER :: i, j, k + + IF (ncly.EQ.2) THEN + yc = yly / 2._mytype + zc = zlz / 2._mytype + IF (zstart(2).EQ.1) THEN + j = 1 + y = -yc + DO k = 1, zsize(3) + z = dz * (k - 1) - zc + + x2 = y + y2 = z + x1 = y + dy + y1 = y2 * x1 / x2 + r1 = SQRT(x1**2 + y1**2) + r2 = SQRT(x2**2 + y2**2) + IF (r1.GT.r2) THEN + PRINT *, "Bug1 in entrainment_bcz" + STOP + ELSE + IF (k.EQ.1) THEN ! First z-point + DO i = 1, zsize(1) + clx3(i, j, k) = ux3(i, j + 1, k + 1) + cly3(i, j, k) = uy3(i, j + 1, k + 1) * r1 / r2 + clz3(i, j, k) = uz3(i, j + 1, k + 1) * r1 / r2 + ENDDO + ELSEIF (k.EQ.((nz - 1) / 2 + 1)) THEN ! Mid z-point + DO i = 1, zsize(1) + clx3(i, j, k) = ux3(i, j + 1, k) + cly3(i, j, k) = uy3(i, j + 1, k) * r1 / r2 + clz3(i, j, k) = uz3(i, j + 1, k) * r1 / r2 + ENDDO + ELSEIF (k.EQ.nz) THEN ! Final z-point + DO i = 1, zsize(1) + clx3(i, j, k) = ux3(i, j + 1, k - 1) + cly3(i, j, k) = uy3(i, j + 1, k - 1) * r1 / r2 + clz3(i, j, k) = uz3(i, j + 1, k - 1) * r1 / r2 + ENDDO + ELSE ! General z-point + IF (z.GT.0._mytype) THEN + ya = y2 - dz + DO i = 1, zsize(1) + uu1 = ux3(i, j + 1, k - 1) & + + (ux3(i, j + 1, k) - ux3(i, j + 1, k - 1)) * (y1 - ya) / (y2 - ya) + uv1 = uy3(i, j + 1, k - 1) & + + (uy3(i, j + 1, k) - uy3(i, j + 1, k - 1)) * (y1 - ya) / (y2 - ya) + uw1 = uz3(i, j + 1, k - 1) & + + (uz3(i, j + 1, k) - uz3(i, j + 1, k - 1)) * (y1 - ya) / (y2 - ya) + + clx3(i, j, k) = uu1 + cly3(i, j, k) = uv1 * r1 / r2 + clz3(i, j, k) = uw1 * r1 / r2 + ENDDO + ELSEIF (z.LT.0._mytype) THEN + ya = y2 + dz + DO i = 1, zsize(1) + uu1 = ux3(i, j + 1, k + 1) & + + (ux3(i, j + 1, k + 1) - ux3(i, j + 1, k)) * (y1 - ya) / (ya - y2) + uv1 = uy3(i, j + 1, k + 1) & + + (uy3(i, j + 1, k + 1) - uy3(i, j + 1, k)) * (y1 - ya) / (ya - y2) + uw1 = uz3(i, j + 1, k + 1) & + + (uz3(i, j + 1, k + 1) - uz3(i, j + 1, k)) * (y1 - ya) / (ya - y2) + + clx3(i, j, k) = uu1 + cly3(i, j, k) = uv1 * r1 / r2 + clz3(i, j, k) = uw1 * r1 / r2 + ENDDO + ELSE ! z = 0 + ENDIF + ENDIF ! End checking where in z we are + ENDIF + ENDDO !! End loop over k + ENDIF !! End if (ztstart(2).eq.1), i.e. first y boundary + + IF (zend(2).EQ.ny) THEN + j = zsize(2) + y = yc + DO k = 1, zsize(3) + z = dz * (k - 1) - zc + x2 = y + y2 = z + x1 = y - dy + y1 = y2 * x1 / x2 + r1 = SQRT(x1**2 + y1**2) + r2 = SQRT(x2**2 + y2**2) + IF (r1.GT.r2) THEN + PRINT *, "Bug2 in entrainment_bcz" + STOP + ELSE + IF (k.EQ.1) THEN ! First z-point + DO i = 1, zsize(1) + clx3(i, j, k) = ux3(i, j - 1, k + 1) + cly3(i, j, k) = uy3(i, j - 1, k + 1) * r1 / r2 + clz3(i, j, k) = uz3(i, j - 1, k + 1) * r1 / r2 + ENDDO + ELSEIF (k.EQ.(nz - 1) / 2 + 1) THEN ! Middle z-point + DO i = 1, zsize(1) + clx3(i, j, k) = ux3(i, j - 1, k) + cly3(i, j, k) = uy3(i, j - 1, k) * r1 / r2 + clz3(i, j, k) = uz3(i, j - 1, k) * r1 / r2 + ENDDO + ELSEIF (k.EQ.nz) THEN ! Last z-point + DO i = 1, zsize(1) + clx3(i, j, k) = ux3(i, j - 1, k - 1) + cly3(i, j, k) = uy3(i, j - 1, k - 1) * r1 / r2 + clz3(i, j, k) = uz3(i, j - 1, k - 1) * r1 / r2 + ENDDO + ELSE ! General z-point + IF (z.GT.0._mytype) THEN + ya = y2 - dz + DO i = 1, zsize(1) + uu1 = ux3(i, j - 1, k - 1) & + +(ux3(i, j - 1, k) - ux3(i, j - 1, k - 1)) * (y1 - ya) / (y2 - ya) + uv1 = uy3(i, j - 1, k - 1) & + + (uy3(i, j - 1, k) - uy3(i, j - 1, k - 1)) * (y1 - ya) / (y2 - ya) + uw1 = uz3(i, j - 1, k - 1) & + + (uz3(i, j - 1, k) - uz3(i, j - 1, k - 1)) * (y1 - ya) / (y2 - ya) + + clx3(i, j, k) = uu1 + cly3(i, j, k) = uv1 * r1 / r2 + clz3(i, j, k) = uw1 * r1 / r2 + ENDDO + ELSEIF (z.LT.0._mytype) THEN + ya = y2 + dz + DO i = 1, zsize(1) + uu1 = ux3(i, j - 1, k + 1) & + + (ux3(i, j - 1, k + 1) - ux3(i, j - 1, k)) * (y1 - ya) / (ya - y2) + uv1 = uy3(i, j - 1, k + 1) & + + (uy3(i, j - 1, k + 1) - uy3(i, j - 1, k)) * (y1 - ya) / (ya - y2) + uw1 = uz3(i, j - 1, k + 1) & + + (uz3(i, j - 1, k + 1) - uz3(i, j - 1, k)) * (y1 - ya) / (ya - y2) + + clx3(i, j, k) = uu1 + cly3(i, j, k) = uv1 * r1 / r2 + clz3(i, j, k) = uw1 * r1 / r2 + ENDDO + ELSE ! z = 0 + ENDIF + ENDIF !! End checking where we are in k + ENDIF !! End error check + ENDDO !! End loop over k + ENDIF !! End if (zend(2).eq.ny), i.e. last y boundary + ENDIF +ENDSUBROUTINE entrainment_bcz diff --git a/incompact3d.prm b/incompact3d.prm index 618dcea..22a6d22 100644 --- a/incompact3d.prm +++ b/incompact3d.prm @@ -7,6 +7,9 @@ 400. #re # Reynolds number 1. #pr # Prandtl number 1. #sc # Schmidt number (if passive scalar) +0. #frx # Froude number in x direction (frx <= 0 gives no gravity) +0. #fry # Froude number in y direction +0. #frz # Froude number in z direction 1.0 #u1 # u1 (max velocity) (for inflow condition) 1.0 #u2 # u2 (min velocity) (for inflow condition) 1.0 #dens1 # dens1 (max density) (for inflow condition) diff --git a/module_param.f90 b/module_param.f90 index ec68e1e..d070625 100644 --- a/module_param.f90 +++ b/module_param.f90 @@ -141,7 +141,7 @@ module param integer :: iscalar, ilmn integer :: nxboite, istat,iread,iadvance_time real(mytype) :: xlx,yly,zlz,dx,dy,dz,dx2,dy2,dz2 - real(mytype) :: dt,xnu,noise,noise1,pi,twopi,u1,u2,sc,pr,dens1,dens2 + real(mytype) :: dt,xnu,noise,noise1,pi,twopi,u1,u2,sc,pr,dens1,dens2,frx,fry,frz real(mytype) :: t,xxk1,xxk2 integer :: itr,itime character :: filesauve*80, filenoise*80, & diff --git a/parameters.f90 b/parameters.f90 index aae7025..d5e70e8 100644 --- a/parameters.f90 +++ b/parameters.f90 @@ -67,6 +67,9 @@ subroutine parameter() read (10,*) re read (10,*) pr read (10,*) sc +read (10,*) frx +read (10,*) fry +read (10,*) frz read (10,*) u1 read (10,*) u2 read (10,*) dens1 @@ -145,6 +148,7 @@ subroutine parameter() print *,'Passive scalar : on' write (*,1113) sc endif +write(*,1114) frx,fry,frz if (ivirt.eq.0) print *,'Immersed boundary : off' if (ivirt.eq.1) then print *,'Immersed boundary : on old school' @@ -170,17 +174,18 @@ subroutine parameter() STOP endif - 1101 format(' Spatial Resolution: (nx,ny,nz)=(',I4,',',I4,',',I4,')') - 1102 format(' Boundary condition: (nclx,ncly,nclz)=(',I1,',',I1,',',I1,')') - 1103 format(' Domain dimension : (lx,ly,lz)=(',F6.1,',',F6.1,',',F6.1,')') - 1104 format(' High and low speed: u1=',F6.2,' and u2=',F6.2) - 1105 format(' High and low density: dens1=',F6.2,' and dens2=',F6.2) - 1106 format(' Reynolds number Re: ',F15.8) - 1107 format(' Time step dt : ',F15.8) - 1108 format(' Object centred at : (',F6.2,',',F6.2,',',F6.2,')') - 1110 format(' Object length : ',F6.2) - 1112 format(' Prandtl number : ',F6.2) - 1113 format(' Schmidt number : ',F6.2) +1101 format(' Spatial Resolution: (nx,ny,nz)=(',I4,',',I4,',',I4,')') +1102 format(' Boundary condition: (nclx,ncly,nclz)=(',I1,',',I1,',',I1,')') +1103 format(' Domain dimension : (lx,ly,lz)=(',F6.1,',',F6.1,',',F6.1,')') +1104 format(' High and low speed: u1=',F6.2,' and u2=',F6.2) +1105 format(' High and low density: dens1=',F6.2,' and dens2=',F6.2) +1106 format(' Reynolds number Re: ',F15.8) +1107 format(' Time step dt : ',F15.8) +1108 format(' Object centred at : (',F6.2,',',F6.2,',',F6.2,')') +1110 format(' Object length : ',F6.2) +1112 format(' Prandtl number : ',F6.2) +1113 format(' Schmidt number : ',F6.2) +1114 format(' Froude number : (frx,fry,frz)=(',F6.2,',',F6.2,',',F6.2,')') endif xnu=1._mytype/re