c **** ATTENTION : THIS PROGRAM USES FORTRAN NAME UP TO 10 CHARACTERS LONG *** C c c c c c c c c o..............o c : : c a(i-1,j) : c u(i,j) : c *------>------* S11, S22(i,j), pression, divergence c | : | u**2, v**2, cc(i,j) c | : | : c | S12(i,j)| v(i,j): c ^ o..... ^.......o c | | c(i,j-1) (c IS NOT cc) [a, c \simeq 1/\rho] c | | c | | c *------>------* c c c SKETCH OF AXISYMMETRIC BOUNDARY CONDITIONS c c !!!AXIS!!! c ! c : c -----*---:--*------*------* j = ny p(i,ny) c | : | | | c .... .....^......^......^......^.... v(i,ny) = 0 c | : | | | c ---------:------------>--- c | : | | | j = ny-1 c | : | | ^ c ---------:----------------| c | : | | | c c c | : | | | c ---------:---------------- c | : | | | c | : | | u(n-1,2) c -----*--->--*--->--*-->-- * p(n-1,2) | Boundary conditions c | : | | | | u(i,1) = u(i,2) c .... ............^......^......^... y = 0 v(i,2) | p(i,1) = p(i,2) c | : |p(2,1)| | | v(i,2) = FIXED(i) c -->--*--->--*--->--*-->-- * j=1 S12(i,2) = 0 c u(1,1)|u(2,1)| | | p(n-1,1) S12(i,ny) = 0 c ^ ^ ^ ^ v(n-1,1) c v(1,1) v(2,1) c c Sites on or outside the dotted lines are used to force boundary conditions c There is also one dummy column on each side (not represented). c c c u,v(i,j) the physical velocities *tau/h c c SUBROUTINE momentum(dcxt,dcyt,u,v,cc,p,a,c,ei, * S11,S22,S12,w1,w2,tau,mugas,muliq, * rhogas,g,sigma,nx,ny,outdiv,t,t0) include 'undefined.h' INTEGER nx,ny,outdiv,t,ifetch,t0 DOUBLE PRECISION tau, h, mugas, muliq, rhogas DOUBLE PRECISION u(nx,ny), v(nx,ny), cc(nx,ny), p(nx,ny) DOUBLE PRECISION S11(nx,ny),S12(nx,ny),S22(nx,ny) DOUBLE PRECISION w1(nx,ny),w2(nx,ny) DOUBLE PRECISION a(nx,ny), c(nx,ny), ei(nx,ny) INTEGER i,j,imax,jmax,cycles,ncycle,npre,npost,horgrav,relax_times * ,printres,relax_times_old,ncycle_old,nres DOUBLE PRECISION mucst1ctr, mucst2ctr, mucst1, mucst2 DOUBLE PRECISION rhocst1, rhocst1ov2,norm_old DOUBLE PRECISION hi2, sumfield, fmodmax1 DOUBLE PRECISION tauh2,g,tau2h,sigma,sigmad,normmean DOUBLE PRECISION s1,s2,t1,t2 DOUBLE PRECISION findninf,mag,mag2 DOUBLE PRECISION dcxt(nx,ny),dcyt(nx,ny),dcx,dcy DOUBLE PRECISION epsilon,mineps,maxeps,xstdgrav,ystdgrav COMMON /cycling/ mineps,maxeps,cycles,ncycle,npre,npost,nres, *relax_times,norm_old,printres,relax_times_old,ncycle_old c horgrav=0 i = ifetch("horgrav",horgrav) c c Numerical constants c h=1.d0/(nx-2) tauh2 = tau/(h*h) tau2h = tau**2/h hi2 = 1.d0/(h*h) c c Constants used to calculated local mu c Mu constant 1, centered c sigmad=(tau**2)*sigma/h**3 mucst1ctr = tauh2*mugas c c Mu constant 2, centered c mucst2ctr = 0.25d0*tauh2*( muliq - mugas) c c Mu constants 1 and 2, not centered c mucst1 = 2.d0*tauh2*mugas mucst2 = 2.d0*tauh2*(muliq - mugas) c c Constants used to calculate local rho c rhocst1 = 1.d0 - rhogas rhocst1ov2 = .5d0*( 1.d0 - rhogas) c write(6,*) "rhocst1ov2 =",rhocst1ov2 c c Gravity constants c xstdgrav = g*tau2h*horgrav ystdgrav = g*tau2h*(1-horgrav) do j=2,ny do i=2,nx S12(i,j) = - S12(i,j) * +(mucst2ctr*(cc(i,j)+cc(i-1,j)+cc(i,j-1)+cc(i-1,j-1)) * + mucst1ctr) * * (u(i,j) - u(i,j-1) + v(i,j) - v(i-1,j)) enddo enddo do j=2,ny-1 do i=2,nx-1 S11(i,j) = - S11(i,j) * + (mucst2*cc(i,j) + mucst1)*(u(i+1,j) - u(i,j)) S22(i,j) = - S22(i,j) * + (mucst2*cc(i,j) + mucst1)*(v(i,j+1) - v(i,j)) enddo enddo c call makebcsf(S11,S22,S12,nx,ny) c call afromc(a,c,ei,cc,rhocst1ov2,rhogas,nx,ny) c D.G. 4 nov 96 , we use T = - sigma*|dc|*(I - n@n), we calculate its divergence in c cartesian coordinates and we add the new axisymetric term -sigma/r nr |dc| n , cc is a pointer c on filtered (if nfilter > 0) or unfiltered (if nfilter = 0) volume fraction. do j=2,ny do i=2,nx dcxt(i,j)= 0.5* * (cc(i,j) + cc(i,j-1) - cc(i-1,j) - cc(i-1,j-1)) dcyt(i,j)= 0.5* * (cc(i,j) + cc(i-1,j) - cc(i,j-1) - cc(i-1,j-1)) enddo enddo do j=2,ny-1 do i=3,nx-1 dcx=0.5d0*(dcxt(i,j)+dcxt(i,j+1)) dcy=0.5d0*(dcyt(i,j)+dcyt(i,j+1)) mag = dsqrt(dcx**2 + dcy**2 + 1.d-50) c Now w1 and cc are radial and axial velocities. c cc is temporary the filtered volume fraction w1(i,j) = u(i,j) + a(i-1,j)*(S12(i,j+1) - S12(i,j) * + S11(i,j) - S11(i-1,j) * + (0.5d0/(i-(2.d0)))*( S11(i,j) + S11(i-1,j) * - 2*sigmad*mag * - (mucst2*0.5d0*(cc(i,j) + cc(i-1,j)) + mucst1) * *(2.d0/(i-2.d0))*u(i,j) )) w2(i,j) = v(i,j) + c(i,j-1)*(S12(i+1,j) -S12(i,j) * + S22(i,j) - S22( i,j-1) + *(0.5d0/(i-1.5d0))*(S12(i+1,j) + S12(i,j)) ) enddo enddo do j=2,ny-1 w2(2,j) = v(2,j) + c(2,j-1)*(S12(3,j) * + S22(2,j) - S22(2,j-1) + S12(3,j) ) w1(2,j) =0.d0 enddo do j=2,ny-1 do i=2,nx-1 w2(i,j) = w2(i,j) + ystdgrav enddo enddo call makebcf(w1,w2,nx,ny) C tauh4 = 0.25d0*tau/h do j=2,ny-1 do i=3,nx-1 u(i,j) = w1(i,j) + 0.25d0*( * (w1(i,j) + w1(i,j-1))*(w2(i,j) + w2(i-1,j)) * - (w1(i,j+1) + w1(i,j))*(w2(i,j+1) + w2(i-1,j+1)) * + (w1(i-1,j)+w1(i,j))**2 - (w1(i,j)+w1(i+1,j))**2 * - (1.d0/(i-2.d0))*((2.d0*w1(i,j))**2) ) enddo enddo do j=2,ny-1 do i=2,nx-1 v(i,j) = w2(i,j) + 0.25d0*( * (w1(i,j) + w1(i,j-1))*(w2(i,j) + w2(i-1,j)) * - (w1(i+1,j) + w1(i+1,j-1))*(w2(i+1,j) + w2(i,j)) * + (w2(i,j-1)+w2(i,j))**2 - (w2(i,j)+w2(i,j+1))**2 * - (1.d0/(i-1.5d0)) * *(w1(i,j)+w1(i+1,j)+w1(i+1,j-1)+w1(i,j-1))*w2(i,j) ) enddo enddo c call makebcf(u,v,nx,ny) c c Now w2=divergence c w1=pressure c do j=2,ny-1 do i=2,nx-1 w2(i,j) = ((i-1)*u(i+1,j)-(i-2)*u(i,j)+(i-1.5d0)*(v(i,j+1) * - v(i,j)))*hi2 enddo enddo if(t.le.(t0+1)) then write(6,*) "divergence sum=",sumfield(w2,nx,ny) endif if(t.le.(t0+1)) then ncycle=3*cycles else ncycle=cycles endif epsilon=fmodmax1(w2,nx,ny,imax,jmax) if(mod(t,outdiv).eq.0) then write(6,*) "t,ncycle, div before=",t,ncycle, $ normmean(w2,nx,ny) * ," div2 = ",epsilon endif call mglin(rhocst1ov2,rhogas,nx,ny) do j=2,ny-1 do i=2,nx-1 c print*,'i,j,p',i,j,p(i,j),u(i,j),v(i,j) c write(*,*) 'i,j,p',i,j,p(i,j) u(i,j) = u(i,j) - (p(i,j) - p(i-1,j))*a(i-1,j) v(i,j) = v(i,j) - (p(i,j) - p(i,j-1))*c(i,j-1) c write(6,*) 'a(',i,',',j,')=',a(i,j) enddo enddo call makebcf(u,v,nx,ny) do j=2,ny-1 do i=2,nx-1 w2(i,j) = ((i-1)*u(i+1,j)-(i-2)*u(i,j)+(i-1.5d0)*(v(i,j+1) * - v(i,j)))*hi2 enddo enddo epsilon=fmodmax1(w2,nx,ny,imax,jmax) if(mod(t,outdiv).eq.0) then write(6,*) "t,ncycles made=,div=",t,cycles, * normmean(w2,nx,ny) * ," div2 = ",epsilon, "imax, jmax",imax,jmax endif return c 100 format(g10.7) end