diff --git a/input/input_example.dat b/input/input_example.dat
index affaac836788797f7aafdf4bced22277d0ff59d1..1d49f6eadc9a41fd073b8d903c18527caf486392 100755
--- a/input/input_example.dat
+++ b/input/input_example.dat
@@ -12,6 +12,8 @@ start time-step (-), time-steps to compute (-) =
1,10
interval of writing result files, restart files (-,-) =
10,10000
+turbulence model [-/Smagorinsky] =
+Smagorinsky
CFL number (-), max number interations (-), relative tolerance (-) =
0.4,200,5.E-3
in/initial flow (m/s), pressure gradient (N/m**3) =
diff --git a/src/bc.f90 b/src/bc.f90
index a727a9c1302302887f5ba72abfdd9ca6a94f6db4..4d98adab6c0040771b83b2577b6597711915389d 100755
--- a/src/bc.f90
+++ b/src/bc.f90
@@ -26,7 +26,7 @@
2 enddo
endif
elseif (bcx.eq.'p') then
- call sync(u); call sync(v); call sync(w) ! done through the sync routines
+ call sync(myu); call sync(myv); call sync(myw) ! done through the sync routines
elseif (bcx.eq.'i') then
! see inflow routine
if (myid.eq.nrprocs-1) then
@@ -85,7 +85,6 @@
!> @brief inflow boundary condition
subroutine inflow
use var
- use parallel
real(kind=pr),dimension(dimj*diml) :: random
real(kind=pr) :: alpha,dist
integer :: i,j,l,m,mm
diff --git a/src/electrostatics.f90 b/src/electrostatics.f90
index 7ed97e19fb7243161e3bb80a4c60dc00a8278104..dfc80284e600c7e78f8aa70761b77a97e7add82a 100755
--- a/src/electrostatics.f90
+++ b/src/electrostatics.f90
@@ -101,7 +101,7 @@
!####################################################################
!> @author Holger Grosshans
-!> @brief compute electric forces on particles from E-field
+!> @brief compute specific electric forces on particles from E-field
!> at seeding time-step E is not available, only Coulomb forces
subroutine forcesGauss
use var
@@ -126,7 +126,7 @@
!####################################################################
!> @author Holger Grosshans
-!> @brief compute electric forces on particles from Coulomb's law
+!> @brief compute specific electric forces on particles from Coulomb's law
subroutine forcesCoulomb
use var
use parallel
diff --git a/src/fluid.f90 b/src/fluid.f90
index 30cab7f97af1550f3f282afd98d8c098826c9afd..c30d1e7624e8fb78db946ebf590e767c7ba5d72f 100755
--- a/src/fluid.f90
+++ b/src/fluid.f90
@@ -5,14 +5,14 @@
subroutine solveFluid
use var
use parallel
- real(kind=pr),dimension(ii,jj,ll) :: du,dv,dw,dp,ddp,massRes
+ real(kind=pr),dimension(ii,jj,ll) :: du,dv,dw,dp,massRes
real(kind=pr) :: err(4),err0(4),errL2
integer itout,it,itmommax,itpCorr1max,itpCorr2max,velCorr1,velCorr2
!> itmax: outer iterations
itmommax=1; itpCorr1max=1; itpCorr2max=0 ! sequential SIMPLE
- !itmommax=4; itpCorr1max=4; itpCorr2max=0 ! SIMPLE
- !itmommax=4; itpCorr1max=4; itpCorr2max=4 ! PISO
+ !itmommax>1; itpCorr1max>1; itpCorr2max=0 ! SIMPLE
+ !itmommax>1; itpCorr1max>1; itpCorr2max>1 ! PISO
velCorr1=0; velCorr2=0 ! velocity correction yes/no (1/0)
if ((bcx.eq.'i').and.(myid.eq.0)) call inflow
@@ -20,11 +20,12 @@
do 1 itout=1,itmax
if ((mod(it,10).eq.0).and.(npTot.ne.0)) call momentumCoupling
+ if (turbModel.eq.'Smagorinsky') call turbSmagorinsky
do 2 it=1,itmommax ! solve momentum eq
call momx(du); call momy(dv); call momz(dw)
- call bcUVW(du,dv,dw)
- u=u+du*urfu; v=v+dv*urfv; w=w+dw*urfw
+ u=u+du*urfu; v=v+dv*urfv; w=w+dw*urfw
+ call bcUVW(u,v,w)
2 enddo
call mass(massRes,u,v,w) ! first pressure correction
@@ -35,10 +36,10 @@
3 enddo
p= p+dp*urfp
- if (velCorr2.eq.1) then ! first velocity correction
+ if (velCorr1.eq.1) then ! first velocity correction
call velocityCorrection(du,dv,dw,dp)
call bcUVW(du,dv,dw)
- u=u+du*urfu; v=v+dv*urfv; w=w+dw*urfw
+ u=u+du*urfu; v=v+dv*urfv; w=w+dw*urfw
endif
if (itpCorr2max.ge.1) then
@@ -48,12 +49,12 @@
call pressureCorrection(dp,massRes)
call bcNeumann(dp)
4 enddo
- p= p+dp*urfp
+ p=p+dp*urfp
if (velCorr2.eq.1) then ! second velocity correction
call velocityCorrection(du,dv,dw,dp)
call bcUVW(du,dv,dw)
- u=u+du*urfu; v=v+dv*urfv; w=w+dw*urfw
+ u=u+du*urfu; v=v+dv*urfv; w=w+dw*urfw
endif
endif
@@ -67,7 +68,7 @@
1 enddo
if (myid.eq.0) write(*,'(a,i3,3(a,es8.2e2))') 'fluid |L2 #',itout-1,' = ',errL2, &
- ' |mom. = ',sum(err(1:3)/err0(1:3))/3._pr,' |pc. = ',err(4)/err0(4)
+ ' |mom. = ',sum(err(1:3)/err0(1:3))/3._pr,' |pc. = ',err(4)/err0(4)
end
@@ -125,7 +126,7 @@
ux= (myu(i,j,l)-myu(i-1,j,l))/(xf(i)-xf(i-1)) ! 2nd order central
vy= (myv(i,j,l)-myv(i,j-1,l))/(yf(j)-yf(j-1))
wz= (myw(i,j,l)-myw(i,j,l-1))/(zf(l)-zf(l-1))
- massRes(i,j,l)=-(ux+vy+wz)
+ massRes(i,j,l)= -(ux+vy+wz)
enddo; enddo; enddo
end
@@ -136,8 +137,7 @@
subroutine momx(du)
use var
real(kind=pr),dimension(ii,jj,ll) :: tu,qu,du
- real(kind=pr) :: ua,va,wa,px,uux,vuy,wuz, &
- uxx,uyy,uzz,dudt,Cdudt,Cviscx,erru
+ real(kind=pr) :: ua,va,wa,px,uux,vuy,wuz,uxx,uyy,uzz,dudt,Cdudt,Cviscx
integer :: i,j,l
tu=0._pr; qu=0._pr; du=0._pr
@@ -150,7 +150,7 @@
va= 0.25_pr*(v(i,j,l)+v(i,j-1,l)+v(i+1,j,l)+v(i+1,j-1,l))
wa= 0.25_pr*(w(i,j,l)+w(i,j,l-1)+w(i+1,j,l)+w(i+1,j,l-1))
-! convective terns (2nd order central)
+! convective terms (2nd order central)
uux= ua*(u(i+1,j,l)-u(i-1,j,l))/(xf(i+1)-xf(i-1))
vuy= va*(u(i,j+1,l)-u(i,j-1,l))/(yc(j+1)-yc(j-1))
wuz= wa*(u(i,j,l+1)-u(i,j,l-1))/(zc(l+1)-zc(l-1))
@@ -171,13 +171,13 @@
(2._pr*tau01*dt)
! momentum eq -lhs+rhs (m/s**2)
- tu(i,j,l)= - dudt -(uux+vuy+wuz) -(px+dpdx)/rhof + nuf*(uxx+uyy+uzz) + Fsx(i,j,l)
+ tu(i,j,l)= - dudt -(uux+vuy+wuz) -(px+dpdx)/rhof + (nuf+nuf_ru(i,j,l))*(uxx+uyy+uzz) + Fsx(i,j,l)
! coefficients of u(i,j,l) of momentum eq
Cdudt= (1._pr+tau01)/(2._pr*tau01*dt)
- Cviscx= -2._pr*nuf * ( 1._pr/((xf(i+1)-xf(i))*(xf(i)-xf(i-1))) &
- + 1._pr/((yc(j+1)-yc(j))*(yc(j)-yc(j-1))) &
- + 1._pr/((zc(l+1)-zc(l))*(zc(l)-zc(l-1))) )
+ Cviscx= -2._pr*(nuf+nuf_ru(i,j,l)) * ( 1._pr/((xf(i+1)-xf(i))*(xf(i)-xf(i-1))) &
+ + 1._pr/((yc(j+1)-yc(j))*(yc(j)-yc(j-1))) &
+ + 1._pr/((zc(l+1)-zc(l))*(zc(l)-zc(l-1))) )
! coefficients lhs-rhs
qu(i,j,l)= Cdudt - Cviscx
@@ -194,8 +194,7 @@
subroutine momy(dv)
use var
real(kind=pr),dimension(ii,jj,ll) :: tv,qv,dv
- real(kind=pr) :: ua,va,wa,py,uvx,vvy,wvz, &
- vxx,vyy,vzz,dvdt,Cdvdt,Cviscy
+ real(kind=pr) :: ua,va,wa,py,uvx,vvy,wvz,vxx,vyy,vzz,dvdt,Cdvdt,Cviscy
integer :: i,j,l
tv=0._pr; qv=0._pr; dv=0._pr
@@ -208,7 +207,7 @@
va= v(i,j,l)
wa= 0.25_pr*(w(i,j,l)+w(i,j,l-1)+w(i,j+1,l)+w(i,j+1,l-1))
-! convective terns (2nd order central)
+! convective terms (2nd order central)
uvx= ua*(v(i+1,j,l)-v(i-1,j,l))/(xc(i+1)-xc(i-1))
vvy= va*(v(i,j+1,l)-v(i,j-1,l))/(yf(j+1)-yf(j-1))
wvz= wa*(v(i,j,l+1)-v(i,j,l-1))/(zc(l+1)-zc(l-1))
@@ -229,13 +228,13 @@
(2._pr*tau01*dt)
! momentum eq. -lhs+rhs
- tv(i,j,l)= - dvdt - (uvx+vvy+wvz) - py/rhof + nuf*(vxx+vyy+vzz) + Fsy(i,j,l)
+ tv(i,j,l)= - dvdt - (uvx+vvy+wvz) - py/rhof + (nuf+nuf_rv(i,j,l))*(vxx+vyy+vzz) + Fsy(i,j,l)
! coefficients of v(i,j,l) of momentum eq
Cdvdt= (1._pr+tau01)/(2._pr*tau01*dt)
- Cviscy= -2._pr*nuf * ( 1._pr/((xc(i+1)-xc(i))*(xc(i)-xc(i-1))) &
- + 1._pr/((yf(j+1)-yf(j))*(yf(j)-yf(j-1))) &
- + 1._pr/((zc(l+1)-zc(l))*(zc(l)-zc(l-1))) )
+ Cviscy= -2._pr*(nuf+nuf_rv(i,j,l)) * ( 1._pr/((xc(i+1)-xc(i))*(xc(i)-xc(i-1))) &
+ + 1._pr/((yf(j+1)-yf(j))*(yf(j)-yf(j-1))) &
+ + 1._pr/((zc(l+1)-zc(l))*(zc(l)-zc(l-1))) )
! coefficients lhs-rhs
qv(i,j,l)= Cdvdt - Cviscy
@@ -252,8 +251,7 @@
subroutine momz(dw)
use var
real(kind=pr),dimension(ii,jj,ll) :: tw,qw,dw
- real(kind=pr) :: ua,va,wa,pz,uwx,vwy,wwz, &
- wxx,wyy,wzz,dwdt,Cdwdt,Cviscz
+ real(kind=pr) :: ua,va,wa,pz,uwx,vwy,wwz,wxx,wyy,wzz,dwdt,Cdwdt,Cviscz
integer :: i,j,l
tw=0._pr; qw=0._pr; dw=0._pr
@@ -266,7 +264,7 @@
va= 0.25_pr*(v(i,j,l)+v(i,j-1,l)+v(i,j,l+1)+v(i,j-1,l+1))
wa= w(i,j,l)
-! convective terns (2nd order central)
+! convective terms (2nd order central)
uwx= ua*(w(i+1,j,l)-w(i-1,j,l))/(xc(i+1)-xc(i-1))
vwy= va*(w(i,j+1,l)-w(i,j-1,l))/(yc(j+1)-yc(j-1))
wwz= wa*(w(i,j,l+1)-w(i,j,l-1))/(zf(l+1)-zf(l-1))
@@ -287,13 +285,13 @@
(2._pr*tau01*dt)
! momentum eq. -lhs+rhs
- tw(i,j,l)= - dwdt - (uwx+vwy+wwz) - pz/rhof + nuf*(wxx+wyy+wzz) + Fsz(i,j,l)
+ tw(i,j,l)= - dwdt - (uwx+vwy+wwz) - pz/rhof + (nuf+nuf_rw(i,j,l))*(wxx+wyy+wzz) + Fsz(i,j,l)
! coefficients of u(i,j,l) of momentum eq
Cdwdt= (1._pr+tau01)/(2._pr*tau01*dt)
- Cviscz= -2._pr*nuf * ( 1._pr/((xc(i+1)-xc(i))*(xc(i)-xc(i-1))) &
- + 1._pr/((yc(j+1)-yc(j))*(yc(j)-yc(j-1))) &
- + 1._pr/((zf(l+1)-zf(l))*(zf(l)-zf(l-1))) )
+ Cviscz= -2._pr*(nuf+nuf_rw(i,j,l)) * ( 1._pr/((xc(i+1)-xc(i))*(xc(i)-xc(i-1))) &
+ + 1._pr/((yc(j+1)-yc(j))*(yc(j)-yc(j-1))) &
+ + 1._pr/((zf(l+1)-zf(l))*(zf(l)-zf(l-1))) )
! coefficients lhs-rhs
qw(i,j,l)= Cdwdt - Cviscz
@@ -303,3 +301,87 @@
enddo; enddo; enddo
end
+
+!####################################################################
+!> @author Holger Grosshans
+!> @brief calculate turbulent viscosity with static Smagorinsky model
+ subroutine turbSmagorinsky
+ use var
+ real(kind=pr) :: ux,uy,uz,vx,vy,vz,wx,wy,wz,l_s
+ integer :: i,j,l
+
+ do l=lmin,lmax; do j=jmin,jmax; do i=imin,imax
+
+ if (celltype(i+1,j,l).ne.wall) then
+ ux= (u(i+1,j,l)-u(i-1,j,l))/(xf(i+1)-xf(i-1)) ! (2nd order central derivatives)
+ uy= (u(i,j+1,l)-u(i,j-1,l))/(yc(j+1)-yc(j-1))
+ uz= (u(i,j,l+1)-u(i,j,l-1))/(zc(l+1)-zc(l-1))
+ vx= (v(i+1,j,l)+v(i+1,j-1,l)-v(i,j,l)-v(i,j-1,l))/(2._pr*(xc(i+1)-xc(i)))
+ vy= (v(i+1,j,l)+v(i,j,l)-v(i+1,j-1,l)-v(i,j-1,l))/(2._pr*(yf(j)-yf(j-1)))
+ vz= (v(i+1,j,l+1)+v(i,j,l+1)+v(i+1,j-1,l+1)+v(i,j-1,l+1)-v(i+1,j,l-1)-v(i,j,l-1)-v(i+1,j-1,l-1)-v(i,j-1,l)-1)/ &
+ (4._pr*(zc(l+1)-zc(l-1)))
+ wx= (w(i+1,j,l)+w(i+1,j,l-1)-w(i,j,l)-w(i,j,l-1))/(2._pr*(xc(i+1)-xc(i)))
+ wy= (w(i+1,j+1,l)+w(i,j+1,l)+w(i+1,j+1,l-1)+w(i,j+1,l-1)-w(i+1,j-1,l)-w(i,j-1,l)-w(i+1,j-1,l-1)-w(i,j-1,l-1))/ &
+ (4._pr*(yc(j+1)-yc(j-1)))
+ wz= (w(i+1,j,l)+w(i,j,l)-w(i+1,j,l-1)-w(i,j,l-1))/(2._pr*(zf(l)-zf(l-1)))
+
+ l_s= C_s*filteru(i,j,l)
+ nuf_ru(i,j,l)= l_s**2*S(ux,uy,uz,vx,vy,vz,wx,wy,wz) !< Pope Eq. (13.128)
+ endif
+
+ if (celltype(i,j+1,l).ne.wall) then
+ ux= (u(i,j+1,l)+u(i,j,l)-u(i-1,j+1,l)-u(i-1,j,l))/(2._pr*(xf(i)-xf(i-1)))
+ uy= (u(i,j+1,l)+u(i-1,j+1,l)-u(i,j,l)-u(i-1,j,l))/(2._pr*(yc(j+1)-yc(j)))
+ uz= (u(i-1,j+1,l+1)+u(i,j+1,l+1)+u(i-1,j,l+1)+u(i,j,l+1)-u(i-1,j+1,l-1)-u(i,j+1,l-1)-u(i-1,j,l-1)-u(i,j,l-1))/ &
+ (4._pr*(zc(l+1)-zc(l-1)))
+ vx= (v(i+1,j,l)-v(i-1,j,l))/(xc(i+1)-xc(i-1))
+ vy= (v(i,j+1,l)-v(i,j-1,l))/(yf(j+1)-yf(j-1))
+ vz= (v(i,j,l+1)-v(i,j,l-1))/(zc(l+1)-zc(l-1))
+ wx= (w(i+1,j,l)+w(i+1,j+1,l)+w(i+1,j,l-1)+w(i+1,j+1,l-1)-w(i-1,j,l)-w(i-1,j+1,l)-w(i-1,j,l-1)-w(i-1,j+1,l-1))/ &
+ (4._pr*(xc(i+1)-xc(i-1)))
+ wy= (w(i,j+1,l)+w(i,j+1,l-1)-w(i,j,l)-w(i,j,l-1))/(2._pr*(yc(j+1)-yc(j)))
+ wz= (w(i,j,l)+w(i,j+1,l)-w(i,j,l-1)-w(i,j+1,l-1))/(2._pr*(zf(l)-zf(l-1)))
+
+ l_s= C_s*filterv(i,j,l)
+ nuf_rv(i,j,l)= l_s**2*S(ux,uy,uz,vx,vy,vz,wx,wy,wz) !< Pope Eq. (13.128)
+ endif
+
+ if (celltype(i,j,l+1).ne.wall) then
+ ux= (u(i,j,l)+u(i,j,l+1)-u(i-1,j,l)-u(i-1,j,l+1))/(2._pr*(xf(i)-xf(i-1)))
+ uy= (u(i,j+1,l)+u(i,j+1,l+1)+u(i-1,j+1,l)+u(i-1,j+1,l+1)-u(i,j-1,l)-u(i,j-1,l+1)-u(i-1,j-1,l)-u(i-1,j-1,l+1))/ &
+ (4._pr*(yc(j+1)-yc(j-1)))
+ uz= (u(i,j,l+1)+u(i-1,j,l+1)-u(i,j,l)-u(i-1,j,l))/(2._pr*(zc(l+1)-zc(l)))
+ vx= (v(i+1,j,l)+v(i+1,j,l+1)+v(i+1,j-1,l)+v(i+1,j-1,l+1)-v(i-1,j,l)-v(i-1,j,l+1)-v(i-1,j-1,l)-v(i-1,j-1,l+1))/ &
+ (4._pr*(xc(i+1)-xc(i-1)))
+ vy= (v(i,j,l)+v(i,j,l+1)-v(i,j-1,l)-v(i,j-1,l+1))/(2._pr*(yf(j)-yf(j-1)))
+ vz= (v(i,j,l+1)+v(i,j-1,l+1)-v(i,j,l)-v(i,j-1,l))/(2._pr*(yc(i+1)-yc(i)))
+ wx= (w(i+1,j,l)-w(i-1,j,l))/(xc(i+1)-xc(i-1))
+ wy= (w(i,j+1,l)-w(i,j-1,l))/(yc(j+1)-yc(j-1))
+ wz= (w(i,j,l+1)-w(i,j,l-1))/(zf(l+1)-zf(l-1))
+
+ l_s= C_s*filterw(i,j,l)
+ nuf_rw(i,j,l)= l_s**2*S(ux,uy,uz,vx,vy,vz,wx,wy,wz) !< Pope Eq. (13.128)
+ endif
+
+ enddo; enddo; enddo
+
+ contains
+
+ real(kind=pr) function S(ux,uy,uz,vx,vy,vz,wx,wy,wz)
+ !> @brief calculate characteristic filtered rate of strain, Pope Eq. (13.74)
+ real(kind=pr) :: ux,uy,uz,vx,vy,vz,wx,wy,wz,Sij(3,3)
+
+ Sij(1,1)= ux !< filtered rate of strain tensor, Pope Eq. (2.70)
+ Sij(2,2)= vy
+ Sij(3,3)= wz
+ Sij(1,2)= 0.5_pr*(uy+vx)
+ Sij(2,1)= Sij(1,2)
+ Sij(1,3)= 0.5_pr*(uz+wx)
+ Sij(3,1)= Sij(1,3)
+ Sij(2,3)= 0.5_pr*(vz+wy)
+ Sij(3,2)= Sij(2,3)
+ S= sqrt(2._pr*sum(Sij**2))
+
+ end
+
+ end
diff --git a/src/parallel.f90 b/src/parallel.f90
index 404e4bf85371c4280eec00ba8c08cd45ca9f8887..cb2678619ada5e518e208d10be5478d5c3550b6d 100755
--- a/src/parallel.f90
+++ b/src/parallel.f90
@@ -212,7 +212,7 @@ contains
timenow=real(mpi_wtime(),kind=pr)
-! mpi_allreduce does not work with intel compiler
+! mpi_allreduce may not work with intel compiler
call mpi_allreduce(myvar,syncMax,1,mpi_pr,mpi_max,mpi_comm_world,mpierr)
timeend=real(mpi_wtime(),kind=pr)
@@ -220,6 +220,24 @@ contains
end
+!####################################################################
+!> @author Holger Grosshans
+!> @brief compute the min of a scalar over all processors
+ real(kind=pr) function syncMin(myvar)
+ use var
+ use mpi
+ real(kind=pr) :: myvar
+
+ timenow=real(mpi_wtime(),kind=pr)
+
+! mpi_allreduce may not work with intel compiler
+ call mpi_allreduce(myvar,syncMin,1,mpi_pr,mpi_min,mpi_comm_world,mpierr)
+
+ timeend=real(mpi_wtime(),kind=pr)
+ timecom(1)=timecom(1)+timeend-timenow
+
+ end
+
!####################################################################
!> @author Holger Grosshans
!> @brief compute the sum of a real scalar over all processors
diff --git a/src/particles.f90 b/src/particles.f90
index 2808a91ac7eb1596be6563f6222e9786fd01516a..fcad544580fe60cc9dd2ffd58c3827264b28198c 100755
--- a/src/particles.f90
+++ b/src/particles.f90
@@ -1,24 +1,3 @@
-!####################################################################
-!> @author Holger Grosshans
-!> @brief solve particulate phase
- subroutine solveParticles
- use var
-
- if (pnd.ne.0.or.npTot.ne.0) then
- if (((bcx.eq.'i').and.(nt.ge.ntseed)).or.(nt.eq.ntseed)) call particlesSeed
- call particlesDrag
- call particlesLift
- if ((elForceScheme.eq.1).or.(elForceScheme.eq.3)) call forcesGauss
- if ((elForceScheme.eq.2).or.(elForceScheme.eq.3)) call forcesCoulomb
- call particlesVelocityNext
- call particlesCollideNext
- call particlesTransportNext
- if (myid.eq.0) write(*,'(a,i8,2(a,es8.2e2))') &
- 'particles |transp. = ',npTot,' |dt/t_el = ',dtNext*rtau_el_max,' |dt/t_p = ',dtNext*rtau_p_max
- endif
-
- end
-
!####################################################################
!> @author Holger Grosshans
!> @brief compute weigths for interpolation from Lagrangian to
@@ -60,13 +39,13 @@
! volume on the Eulerian grid
if (direction.eq.0) then
- volE=(xf(iend)-xf(ibeg-1))*(yf(jend)-yf(jbeg-1))*(zf(lend)-zf(lbeg-1))
+ volE= sum(vol(ibeg:iend,jbeg:jend,lbeg:lend))
else if (direction.eq.1) then
- volE=(xc(iend+1)-xc(ibeg))*(yf(jend)-yf(jbeg-1))*(zf(lend)-zf(lbeg-1))
+ volE= sum(volu(ibeg:iend,jbeg:jend,lbeg:lend))
else if (direction.eq.2) then
- volE=(xf(iend)-xf(ibeg-1))*(yc(jend+1)-yc(jbeg))*(zf(lend)-zf(lbeg-1))
+ volE= sum(volv(ibeg:iend,jbeg:jend,lbeg:lend))
else if (direction.eq.3) then
- volE=(xf(iend)-xf(ibeg-1))*(yf(jend)-yf(jbeg-1))*(zc(lend+1)-zc(lbeg))
+ volE= sum(volw(ibeg:iend,jbeg:jend,lbeg:lend))
endif
deltaf=xc(iend)-xc(ibeg)
diff --git a/src/particlesTransport.f90 b/src/particlesTransport.f90
index b24d4064f79452f15c4dac62faaff89d736a9552..b8428d4d4dc7c84ddb5a6e4001f6f5de810cb744 100755
--- a/src/particlesTransport.f90
+++ b/src/particlesTransport.f90
@@ -1,3 +1,24 @@
+!####################################################################
+!> @author Holger Grosshans
+!> @brief solve particulate phase
+ subroutine solveParticles
+ use var
+
+ if (pnd.ne.0.or.npTot.ne.0) then
+ if (((bcx.eq.'i').and.(nt.ge.ntseed)).or.(nt.eq.ntseed)) call particlesSeed
+ call particlesDrag
+ call particlesLift
+ if ((elForceScheme.eq.1).or.(elForceScheme.eq.3)) call forcesGauss
+ if ((elForceScheme.eq.2).or.(elForceScheme.eq.3)) call forcesCoulomb
+ call particlesVelocityNext
+ call particlesCollideNext
+ call particlesTransportNext
+ if (myid.eq.0) write(*,'(a,i8,2(a,es8.2e2))') &
+ 'particles |transp. = ',npTot,' |dt/t_el = ',dtNext*rtau_el_max,' |dt/t_p = ',dtNext/tau_p_min
+ endif
+
+ end
+
!####################################################################
!> @author Holger Grosshans
!> @brief explicit first-order Euler forward integration of particle
@@ -25,31 +46,31 @@
!####################################################################
!> @author Holger Grosshans
-!> @brief particle drag force
+!> @brief specific particle drag force
subroutine particlesDrag
use var
use parallel
- real(kind=pr) :: rtau_p
+ real(kind=pr) :: tau_p
integer :: n
call fluidVelocity
- rtau_p_max=0._pr
+ tau_p_min=1.e10_pr
do 1 n=1,np
- call calcRtau_p(rtau_p,n)
- fx_d(n)= (uf(n)-up(n))*rtau_p
- fy_d(n)= (vf(n)-vp(n))*rtau_p
- fz_d(n)= (wf(n)-wp(n))*rtau_p
+ call calcTau_p(tau_p,n)
+ fx_d(n)= (uf(n)-up(n))/tau_p
+ fy_d(n)= (vf(n)-vp(n))/tau_p
+ fz_d(n)= (wf(n)-wp(n))/tau_p
1 enddo
- rtau_p_max=syncMax(rtau_p_max)
+ tau_p_min=syncMin(tau_p_min)
end
!####################################################################
!> @author Holger Grosshans
-!> @brief particle lift force (Saffman, 1965, 1968, correction by Mei, 1992)
+!> @brief specific particle lift force (Saffman, 1965, 1968, correction by Mei, 1992)
subroutine particlesLift
use var
real(kind=pr) :: urel,Reyp,Reyf,beta,Cls
@@ -78,25 +99,26 @@
!####################################################################
!> @author Holger Grosshans
-!> @brief reciprocal (to avoid NaN) particle response time (Putnam, 1961)
- subroutine calcRtau_p(rtau_p,n)
+!> @brief particle response time
+ subroutine calcTau_p(tau_p,n)
use var
- real(kind=pr) :: urel,Reyp,Cd,rtau_p
+ real(kind=pr) :: urel,Reyp,Cd,tau_p
integer :: n
urel= sqrt((uf(n)-up(n))**2+(vf(n)-vp(n))**2+(wf(n)-wp(n))**2)
- Reyp= max(1.e-10_pr,2._pr*radp(n)*urel/nuf)
+ urel= max(1.e-10_pr,urel) ! avoid NaN
+ Reyp= 2._pr*radp(n)*urel/nuf
- if (Reyp.lt.1000._pr) then
+ if (Reyp.lt.1000._pr) then ! Putnam (1961)
Cd= 24._pr/Reyp * (1._pr + 1._pr/6._pr*Reyp**(2._pr/3._pr))
else
Cd= 0.424_pr
endif
- rtau_p= 3._pr*rhof*urel*Cd/(8._pr*rhop*radp(n))
+ tau_p= 8._pr*rhop*radp(n)/(3._pr*rhof*urel*Cd)
- rtau_p_max=max(rtau_p_max,rtau_p)
- if (rtau_p.gt.(1._pr/dtNext)) rtau_p= 1._pr/dtNext ! stability criteria Euler forward
+ tau_p_min=min(tau_p_min,tau_p)
+ if (tau_p.lt.dtNext) tau_p= dtNext ! stability criteria Euler forward
end
@@ -399,13 +421,14 @@
npp=np
if (bcx.eq.'i') then
if (myid.eq.0) then
- npseed=int(dt*pnd)
+ npseed=int((t-tstart)*pnd)-npstart
+ npstart=npstart+npseed
else
npseed=0
endif
else
npseed=int(pnd*dimx*dimy*dimz)
- ! npseed=2
+! npseed=2
endif
npseedTot=syncSumI(npseed)
np=npp+npseed
diff --git a/src/post.f90 b/src/post.f90
index 9caaa17d0f2df7fe783e4f48e9a78627d145b938..cacd8cdb7ba71bc6b908bcd313926ca66aba1581 100755
--- a/src/post.f90
+++ b/src/post.f90
@@ -22,10 +22,10 @@
subroutine monitor
use var
use parallel
- real(kind=pr) :: ucl,avu,avv,avw,tke,gradu,graduy,graduz, &
+ real(kind=pr) :: ucl,avu,avv,avw,gradu,graduy,graduz, &
Rec,Reb,avvp,avwp,avyp,avzp,avvf,avwf,avxp,avqp, &
rmsv,rmsw,C05,A05
- integer :: i,j,l,m,N05
+ integer :: m,N05
logical :: file_ex
character(70) :: filename2
diff --git a/src/pre.f90 b/src/pre.f90
index d656408f1b6d6b62128e67b03816c9085b3145cf..44c5afeb60e1b044775709c41d7ca1db04411d50 100755
--- a/src/pre.f90
+++ b/src/pre.f90
@@ -45,6 +45,8 @@
read(10,*)
read(10,*) int_results,int_restart
read(10,*)
+ read(10,*) turbModel
+ read(10,*)
read(10,*) cfl,itmax,tol
read(10,*)
read(10,*) ubulk0,dpdx
@@ -117,6 +119,8 @@
'Grid in x,y,z-direction [(u)niform/(s)tretch] = ',gridx,gridy,gridz
write(20,'(a,4(i8))') &
'nt: start, compute, results, restart (-,-,-,-) = ',ntstart,ntime,int_results,int_restart
+ write(20,'(a,3x,a)') &
+ 'turbulence model: = ',turbModel
write(20,'(a,es11.2e2,i11,es11.2e2)') &
'CFL number (-), max number iterations (-), rel. tol. = ',cfl,itmax,tol
write(20,'(a,2(es11.2e2))') &
@@ -163,11 +167,12 @@
Ex_el(ii,jj,ll),Ey_el(ii,jj,ll),Ez_el(ii,jj,ll), &
u01(ii,jj,ll),v01(ii,jj,ll),w01(ii,jj,ll), &
u02(ii,jj,ll),v02(ii,jj,ll),w02(ii,jj,ll), &
- Fsx(ii,jj,ll),Fsy(ii,jj,ll),Fsz(ii,jj,ll))
+ Fsx(ii,jj,ll),Fsy(ii,jj,ll),Fsz(ii,jj,ll), &
+ nuf_ru(ii,jj,ll),nuf_rv(ii,jj,ll),nuf_rw(ii,jj,ll))
! particles
allocate(npic(ii,jj,ll))
- np=0; npp=0; nseedLoc=0;
+ np=0; npp=0; nseedLoc=0; npstart=0
call allocateParticleArrays
! init
@@ -175,9 +180,10 @@
u01=u; v01=v; w01= w
u02=u; v02=v; w02= w
Fsx=0._pr; Fsy=0._pr; Fsz=0._pr
+ nuf_ru=0._pr; nuf_rv=0._pr; nuf_rw=0._pr
celltype=active
t=0._pr; timecom=0._pr
- rtau_p_max=0._pr; Dup_max=0._pr
+ tau_p_min=0._pr; dup_max=0._pr
vol=0._pr; volu=0._pr; volv=0._pr; volw=0._pr;
Ca=0._pr; Cb1=0._pr; Cb2=0._pr; Cc1=0._pr; Cc2=0._pr; Cd1=0._pr; Cd2=0._pr;
@@ -217,7 +223,7 @@
subroutine gridGeneration
use var
real(kind=pr) :: hy_temp,hz_temp
- integer :: i,j,l,m
+ integer :: i,j,l
imin=1+gc; imax=ii-gc
jmin=1+gc; jmax=jj-gc
@@ -335,7 +341,7 @@
use var
integer :: i,j,l
- do i=imin,imax; do j=jmin,jmax; do l=lmin,lmax
+ do i=imin-1,imax+1; do j=jmin-1,jmax+1; do l=lmin-1,lmax+1
vol(i,j,l)= (xf(i)-xf(i-1))*(yf(j)-yf(j-1))*(zf(l)-zf(l-1))
volu(i,j,l)=(xc(i+1)-xc(i))*(yf(j)-yf(j-1))*(zf(l)-zf(l-1))
volv(i,j,l)=(xf(i)-xf(i-1))*(yc(j+1)-yc(j))*(zf(l)-zf(l-1))
@@ -343,10 +349,15 @@
enddo; enddo; enddo
volTot= sum(vol(imin:imax,jmin:jmax,lmin:lmax))
+! LES filter width
+ filteru= volu**(1._pr/3._pr)
+ filterv= volv**(1._pr/3._pr)
+ filterw= volw**(1._pr/3._pr)
+
! coefficients for PISO algorithm, not sure which version is better
do l=lmin,lmax; do j=jmin,jmax; do i=imin,imax
- Cb1(i,j,l)=0._pr; Cb2(i,j,l)=0._pr; Cc1(i,j,l)=0._pr
- Cc2(i,j,l)=0._pr; Cd1(i,j,l)=0._pr; Cd2(i,j,l)=0._pr
+ Cb1(i,j,l)=0._pr; Cc1(i,j,l)=0._pr; Cd1(i,j,l)=0._pr
+ Cb2(i,j,l)=0._pr; Cc2(i,j,l)=0._pr; Cd2(i,j,l)=0._pr
if (celltype(i+1,j,l).ne.wall) Cb1(i,j,l)= 1._pr/((xc(i+1)-xc(i))*(xf(i)-xf(i-1)))
if (celltype(i-1,j,l).ne.wall) Cb2(i,j,l)= 1._pr/((xc(i)-xc(i-1))*(xf(i)-xf(i-1)))
if (celltype(i,j+1,l).ne.wall) Cc1(i,j,l)= 1._pr/((yc(j+1)-yc(j))*(yf(j)-yf(j-1)))
@@ -355,15 +366,6 @@
if (celltype(i,j,l-1).ne.wall) Cd2(i,j,l)= 1._pr/((zc(l)-zc(l-1))*(zf(l)-zf(l-1)))
Ca(i,j,l)= Cb1(i,j,l)+Cb2(i,j,l)+Cc1(i,j,l)+Cc2(i,j,l)+Cd1(i,j,l)+Cd2(i,j,l)
enddo; enddo; enddo
- !do l=lmin,lmax; do j=jmin,jmax; do i=imin,imax
- ! Cb1(i,j,l)= 1._pr/((xc(i+1)-xc(i))*(xf(i)-xf(i-1)))
- ! Cb2(i,j,l)= 1._pr/((xc(i)-xc(i-1))*(xf(i)-xf(i-1)))
- ! Cc1(i,j,l)= 1._pr/((yc(j+1)-yc(j))*(yf(j)-yf(j-1)))
- ! Cc2(i,j,l)= 1._pr/((yc(j)-yc(j-1))*(yf(j)-yf(j-1)))
- ! Cd1(i,j,l)= 1._pr/((zc(l+1)-zc(l))*(zf(l)-zf(l-1)))
- ! Cd2(i,j,l)= 1._pr/((zc(l)-zc(l-1))*(zf(l)-zf(l-1)))
- ! Ca(i,j,l)= Cb1(i,j,l)+Cb2(i,j,l)+Cc1(i,j,l)+Cc2(i,j,l)+Cd1(i,j,l)+Cd2(i,j,l)
- !enddo; enddo; enddo
end
diff --git a/src/restart.f90 b/src/restart.f90
index e2778c0883dbeb6c43436ccb1605c872c63e2484..9130983f34d19cd98cd4245428e20b0d0c1dd389 100755
--- a/src/restart.f90
+++ b/src/restart.f90
@@ -99,5 +99,6 @@
ntstart=ntstart+1
ntend=ntend+1
+ tstart=t
end
diff --git a/src/var.f90 b/src/var.f90
index 33a4ac8d64bb439fe1b0a030bad6b9cddb99528e..b54c9a6da1bac691cb1589e5576d69c0eab5ede2 100755
--- a/src/var.f90
+++ b/src/var.f90
@@ -20,10 +20,11 @@
real(kind=pr), parameter :: &
eps_el= 8.85e-12_pr, & !< permittivity vacuum/air (F/m)
pi= 4._pr*atan(1._pr), &
- urfu= 0.25_pr, & !< under relaxation factor variable u
+ urfu= 0.25_pr, & !< under-relaxation factor variable u
urfv= urfu, &
urfw= urfv, &
- urfp= 4._pr
+ urfp= 4.0_pr, &
+ C_s= 0.17_pr !< Smagorinsky constant
real(kind=pr), parameter :: &
Ew=1.e11_pr, & !< duct Young's modulus (kg/s**2/m)
@@ -66,6 +67,9 @@
dimitot, & !< global number of cells in x-direction
ntseed !< time step particle seeding
+ character(11) :: &
+ turbModel !< turbulence model
+
! run
real(kind=pr) :: &
dpdx, & !< streamwise pressure gradient (N/m**3)
@@ -74,7 +78,8 @@
dt01, & !< time step size from (nt-2) to (nt-1) (s)
dt02, & !< time step size from (nt-3) to (nt-2) (s)
tau01,tau02, &
- t !< physical time (s)
+ t, & !< physical time (s)
+ tstart !< start time step
integer :: nt, & !< time step
ntend !< last time step of computation
@@ -122,11 +127,13 @@
p, & !< fluid pressure (N/m**2)
u01,v01,w01, & !< fluid velocity previous time-step
u02,v02,w02, & !< fluid velocity previous time-step
- Fsx,Fsy,Fsz !< source term particles -> fluid (m/s**2)
+ Fsx,Fsy,Fsz, & !< source term particles -> fluid (m/s**2)
+ nuf_ru,nuf_rv,nuf_rw, & !< turbulent viscosity (m**2/s)
+ filteru,filterv,filterw !< LES filter width (m)
! particles
real(kind=pr) :: &
- rtau_p_max, & !< max reciprocal particle time-scale in a time-step
+ tau_p_min, & !< min particle response time in a time-step
rtau_el_max, & !< max reciprocal Coulomb force time-scale in a time-step
dup_max !< max particle velocity change in a time-step
@@ -135,6 +142,7 @@
npTot, & !< total number of particles
npp, & !< local number of particles before some operation
maxnp, & !< maximum possible local number of particles
+ npstart, & !< number of particles seeded since tstart
nseedLoc !< local counter for seeded particles
real(kind=pr), allocatable, dimension(:) :: &