c*********************************************************************** c SKETCH OF THE C-GRID in 3-D c c S11,S22,S33,divergence,u**2,v**2,cc are defined in the same c position as p, while (i,j,k) are the default indices c c u p c p(i-1) *------------------*------------------* c /. /| c / . S12 / | c v(i-1) * . * v * | c / . / | c / . u(j-1) p(j-1) / | c p(i-1,j-1) *------------------*------------------* | c | . S13 | | c | + w(i-1) + | * w c | . | | c | . | S23 | c S23(i-1) | + . | * | c | . | | c | . S13(j-1) | | c w(i-1,j-1) * . * w(j-1) * | c | . | | c | + . . . . . . . . .+ . . . . . .| . . * p(k-1) c | . p(i-1,k-1) u(k-1) | / c | . | / c | + v(i-1,k-1) + | * v(k-1) c | . S12(k-1) | / c |. |/ c *------------------*------------------* p(j-1,k-1) c p(i-1,j-1,k-1) u(j-1,k-1) c c c*** c c c SKETCH OF THE TOP-BOTTOM BOUNDARY CONDITIONS c c c 1 2 3 N-2 N-1 N c *-------*-------*-- --*-------*-------* k=N; p(i,j,N) c | | | | | | c ....^.......^.......^.. ..^.......^.......^ w(i,j,N)=0 c | | | | u(N-1,j,N-1) | c *-------*-------*-- --*--->---*-------* k=N-1 c | | | | | | c | | | | | | nx=ny=nz=N c | | | | | | c *-------*-------*-- --*-------*-------* c | | | | | | c c c | S12(2,j,3) | | | | c *---x---*-------*-- --*-------*-------* k=3 c | | | | | | c | | | | | | c | | | | | | c *-------*-------*-- --*-------*-------* p(N,j,2) c | S23(2,j,2) | | | | c ....^.......x.......^.. ..^.......^.......^ z=0; w(i,j,2)=0 c | | | | | | c *-------*-------*-- --*-------*-------* k=1; p(i,j,1) c | | | | | | c ^ x ^ ^ ^ ^ ^ c S13(2,j,1) w(N-1,j,1) c c c*********************************************************************** SUBROUTINE momentum(u,v,w,cc,p,S11,S22,S33,S12,S13,S23, w2, w3, c called as u,v,w,zcc,zu,S11,S22,S33,S12,S13,S23,zres,zrhs % a, b, c, ei,w1,t) c called as za,zb,zc,zei,zd,t) c*** INCLUDE 'undefined.h' INCLUDE 'dimension.h' INCLUDE 'physics.h' INCLUDE 'algorithms_controllers.h' INCLUDE 'flags.h' INTEGER i,j,k,t,imax,jmax,kmax DOUBLE PRECISION u(nx,ny,nz), v(nx,ny,nz), w(nx,ny,nz) DOUBLE PRECISION p(nx,ny,nz) DOUBLE PRECISION w1(nx,ny,nz), w2(nx,ny,nz), w3(nx,ny,nz) DOUBLE PRECISION S11(nx,ny,nz),S22(nx,ny,nz),S33(nx,ny,nz) DOUBLE PRECISION S12(nx,ny,nz),S13(nx,ny,nz),S23(nx,ny,nz) DOUBLE PRECISION a(nx,ny,nz), b(nx,ny,nz), c(nx,ny,nz) DOUBLE PRECISION cc(nx,ny,nz), ei(nx,ny,nz) DOUBLE PRECISION tauh2,tau2h,r1,r2,r10,r20,hi2,stdgrav,temp DOUBLE PRECISION horgrav, horgr1, volliq DOUBLE PRECISION epsdiv DOUBLE PRECISION SUMFIELD,FMODMAX EXTERNAL SUMFIELD,FMODMAX INTRINSIC MOD DOUBLE PRECISION FINDN2S EXTERNAL FINDN2S c SUBROUTINES: abcei, bc_tensor, bc_vector, mglin c*** c define constants, renormalize gravity c*** tauh2 = tau/(h*h) tau2h = tau**2/h r20 = tauh2*mu2 r10 = 0.25d0*tauh2*(mu1-mu2) r1 = 2.d0*tauh2*(mu1-mu2) r2 = 2.d0*tauh2*mu2 stdgrav= g * tau2h horgrav = g1 * tau2h c*** c total liquid volume fraction and lift force to compensate c gravity and keep droplets/bubble steady c*** volliq = SUMFIELD(cc,nx,ny,nz)/((nx-2)*(ny-2)*(nz-2)) horgr1 = horgrav*(rho1*volliq + rho2*(1.d0 - volliq)) c*** c add the viscous part to the stress tensor c*** do k=2,nz-1 do j=2,ny-1 do i=2,nx-1 temp = r2 + r1*cc(i,j,k) S11(i,j,k) = - S11(i,j,k) + % temp * (u(i+1,j,k) - u(i,j,k)) S22(i,j,k) = - S22(i,j,k) + % temp * (v(i,j+1,k) - v(i,j,k)) S33(i,j,k) = - S33(i,j,k) + % temp * (w(i,j,k+1) - w(i,j,k)) S12(i,j,k) = - S12(i,j,k) + (r20 + r10 * % (cc(i,j,k)+cc(i-1,j,k)+cc(i,j-1,k)+cc(i-1,j-1,k))) * % (u(i,j,k) - u(i,j-1,k) + v(i,j,k) - v(i-1,j,k)) S13(i,j,k) = - S13(i,j,k) + (r20 + r10 * % (cc(i,j,k)+cc(i-1,j,k)+cc(i,j,k-1)+cc(i-1,j,k-1))) * % (u(i,j,k) - u(i,j,k-1) + w(i,j,k) - w(i-1,j,k)) S23(i,j,k) = - S23(i,j,k) + (r20 + r10 * % (cc(i,j,k)+cc(i,j-1,k)+cc(i,j,k-1)+cc(i,j-1,k-1))) * % (v(i,j,k) - v(i,j,k-1) + w(i,j,k) - w(i,j-1,k)) enddo enddo enddo c*** call bc_tensor(S11,S22,S33,S12,S13,S23,nx,ny,nz) c*** r2 = rho2 r1 = .5d0*(rho1 - rho2) c*** c*** Compute the inverse of rho and put in a,b,c call abcei(a,b,c,ei,cc,r1,r2,nx,ny,nz) c*** c advance velocity field in time: viscous and gravity terms c*** do k=2,nz-1 do j=2,ny-1 do i=2,nx-1 w1(i,j,k) = u(i,j,k) + a(i-1,j,k) * (S11(i,j,k) - % S11(i-1,j,k) + S12(i,j+1,k) - S12(i,j,k) + % S13(i,j,k+1) - S13(i,j,k)) + horgrav % - horgr1/(r1*(cc(i,j,k)+cc(i,j,k-1)) + r2) w2(i,j,k) = v(i,j,k) + b(i,j-1,k) * (S12(i+1,j,k) - % S12(i,j,k) + S22(i,j,k) - S22(i,j-1,k) + % S23(i,j,k+1) - S23(i,j,k)) w3(i,j,k) = w(i,j,k) + c(i,j,k-1) * (S13(i+1,j,k) - % S13(i,j,k) + S23(i,j+1,k) - S23(i,j,k) + % S33(i,j,k) - S33(i,j,k-1)) + stdgrav enddo enddo enddo c*** call bc_vector(w1,w2,w3,nx,ny,nz) c*** c advance velocity field in time: advection term c*** do k=2,nz-1 do j=2,ny-1 do i=2,nx-1 u(i,j,k) = w1(i,j,k) + 0.25d0 * ( % (w1(i,j,k)+w1(i-1,j,k))**2 - (w1(i+1,j,k)+w1(i,j,k))**2 % + (w1(i,j-1,k)+w1(i,j,k))*(w2(i,j,k) +w2(i-1,j,k)) % - (w1(i,j+1,k)+w1(i,j,k))*(w2(i,j+1,k)+w2(i-1,j+1,k)) % + (w1(i,j,k-1)+w1(i,j,k))*(w3(i,j,k) +w3(i-1,j,k)) % - (w1(i,j,k+1)+w1(i,j,k))*(w3(i,j,k+1)+w3(i-1,j,k+1)) ) v(i,j,k) = w2(i,j,k) + 0.25d0 * ( % (w2(i,j,k)+w2(i,j-1,k))**2 - (w2(i,j+1,k)+w2(i,j,k))**2 % + (w2(i-1,j,k)+w2(i,j,k))*(w1(i,j,k) +w1(i,j-1,k)) % - (w2(i+1,j,k)+w2(i,j,k))*(w1(i+1,j,k)+w1(i+1,j-1,k)) % + (w2(i,j,k-1)+w2(i,j,k))*(w3(i,j,k) +w3(i,j-1,k)) % - (w2(i,j,k+1)+w2(i,j,k))*(w3(i,j,k+1)+w3(i,j-1,k+1)) ) w(i,j,k) = w3(i,j,k) + 0.25d0 * ( % (w3(i,j,k)+w3(i,j,k-1))**2 - (w3(i,j,k+1)+w3(i,j,k))**2 % + (w3(i-1,j,k)+w3(i,j,k))*(w1(i,j,k) +w1(i,j,k-1)) % - (w3(i+1,j,k)+w3(i,j,k))*(w1(i+1,j,k)+w1(i+1,j,k-1)) % + (w3(i,j-1,k)+w3(i,j,k))*(w2(i,j,k) +w2(i,j,k-1)) % - (w3(i,j+1,k)+w3(i,j,k))*(w2(i,j+1,k)+w2(i,j+1,k-1)) ) enddo enddo enddo c*** call bc_vector(u,v,w,nx,ny,nz) c*** c now let w3 = (div v) / h c*** hi2=1.d0/(h*h) do k=2,nz-1 do j=2,ny-1 do i=2,nx-1 w3(i,j,k) = (u(i+1,j,k) - u(i,j,k) + v(i,j+1,k) - % v(i,j,k) + w(i,j,k+1) - w(i,j,k)) * hi2 enddo enddo enddo c*** c check accuracy of MG c*** if (t .le. 2) then write (6,*) "t = ",t," divergence sum = ", SUMFIELD(w3,nx,ny,nz) endif if (MOD(t,tcheck) .eq. 0) then write (*,*) "Norm2 div before = ",h**2*FINDN2S(w3,nx,ny,nz) endif c*** c call the pressure solver c*** call mglin() c*** c make the velocity field about divergence free c*** do k=2,nz-1 do j=2,ny-1 do i=2,nx-1 u(i,j,k) = u(i,j,k) - (p(i,j,k)-p(i-1,j,k))*a(i-1,j,k) v(i,j,k) = v(i,j,k) - (p(i,j,k)-p(i,j-1,k))*b(i,j-1,k) w(i,j,k) = w(i,j,k) - (p(i,j,k)-p(i,j,k-1))*c(i,j,k-1) enddo enddo enddo c*** call bc_vector(u,v,w,nx,ny,nz) c*** c check calls c*** if (MOD(t,tcheck) .eq. 0) then do k=2,nz-1 do j=2,ny-1 do i=2,nx-1 w3(i,j,k) = (u(i+1,j,k) - u(i,j,k) + v(i,j+1,k) - % v(i,j,k) + w(i,j,k+1) - w(i,j,k))*hi2 enddo enddo enddo write (*,*) "Norm2 div after = ",h**2*FINDN2S(w3,nx,ny,nz), % " Ncycles = ", cycles," relax_times = ",relax_times endif c*** return end c*********************************************************************** SUBROUTINE abcei (a,b,c,ei,cc,r1,r2,nx,ny,nz) c*** INCLUDE 'undefined.h' INTEGER nx,ny,nz,i,j,k DOUBLE PRECISION a(nx,ny,nz),b(nx,ny,nz),c(nx,ny,nz), % ei(nx,ny,nz),cc(nx,ny,nz) DOUBLE PRECISION r1,r2 c SUBROUTINES: bc_coef c*** do k=2,nz-1 do j=2,ny-1 do i=2,nx-1 a(i,j,k) = 1.d0 / (r1*(cc(i,j,k)+cc(i+1,j,k)) + r2) b(i,j,k) = 1.d0 / (r1*(cc(i,j,k)+cc(i,j+1,k)) + r2) c(i,j,k) = 1.d0 / (r1*(cc(i,j,k)+cc(i,j,k+1)) + r2) enddo enddo enddo c*** c need indices i=1 for a, j=1 for b and k=1 for c to calculate ei c*** call bc_coef(a,b,c,cc,r1,r2,nx,ny,nz) c*** do k=2,nz-1 do j=2,ny-1 do i=2,nx-1 ei(i,j,k) = 1.d0/(a(i,j,k) + a(i-1,j,k) + b(i,j,k) + % b(i,j-1,k) + c(i,j,k) + c(i,j,k-1)) enddo enddo enddo c*** return end