Multiphase SIngle Component Model

Dear All,

I just completed modify the program code from Sukop’s book to model phase separation phenomenon.

But, when I plotted the results, the density field did not change too much from the initial value (200+rand())

I attach here the code.

Can you help me find the bug?

Thank you

  parameter(ly=200,lx=200)
  
  integer is_solid_node(ly,lx)
  real rho(ly,lx),f(ly,lx,9),ftemp(ly,lx,9),ex(9),ey(9),
 +     u_x(ly,lx),u_y(ly,lx),x(ly,lx),y(ly,lx),feq(ly,lx,9)
real psi(ly,lx),forcex(ly,lx),forcey(ly,lx)
real psi_x, psi_y

  tau = 1.

C Set solid nodes at walls on top and bottom

  is_solid_node=0

C Set initial density
do j=1,ly
do i=1,lx
rho(j,i)=200.+ rand()
end do
end do

  do j=1,ly
  do i=1,lx
  f(j,i,1) = (4./9. )*rho(j,i)
  f(j,i,2) = (1./9. )*rho(j,i)
  f(j,i,3) = (1./9. )*rho(j,i)
  f(j,i,4) = (1./9. )*rho(j,i)
  f(j,i,5) = (1./9. )*rho(j,i)
  f(j,i,6) = (1./36.)*rho(j,i)
  f(j,i,7) = (1./36.)*rho(j,i)
  f(j,i,8) = (1./36.)*rho(j,i)
  f(j,i,9) = (1./36.)*rho(j,i)
  enddo
  enddo

C Define lattice velocity vectors

  ex(1)= 0
  ey(1)= 0
  ex(2)= 1
  ey(2)= 0
  ex(3)= 0
  ey(3)= 1
  ex(4)=-1
  ey(4)= 0
  ex(5)= 0
  ey(5)=-1
  ex(6)= 1
  ey(6)= 1
  ex(7)=-1
  ey(7)= 1
  ex(8)=-1
  ey(8)=-1
  ex(9)= 1
  ey(9)=-1

C Time loop

  do ts=1,300

    write(*,*) ts

C Computing macroscopic density, rho, and velocity, u=(ux,uy).
do j=1,ly

    do i=1,lx
        
  u_x(j,i) = 0.0
  u_y(j,i) = 0.0
  rho(j,i) = 0.0

    
        if(is_solid_node(j,i).eq.0) then
            
            do k=1,9
                
                rho(j,i) = rho(j,i) + f(j,i,k)
                                   
                u_x(j,i) = u_x(j,i) + ex(k)*f(j,i,k)
                u_y(j,i) = u_y(j,i) + ey(k)*f(j,i,k)
                
            enddo

            u_x(j,i) = u_x(j,i)/rho(j,i)
            u_y(j,i) = u_y(j,i)/rho(j,i)

        endif
     
  enddo
  enddo

C Compute Force term
do j=1,ly
do i=1,lx
if(is_solid_node(j,i).eq.0) then
psi(j,i) = 4.0*exp(-200./rho(j,i))
end if
end do
end do

do j=1,ly

    if (j.gt.1) then 
        jn = j-1
    else
        jn = LY
    endif
    
    if (j.lt.ly) then
        jp = j+1
    else 
        jp = 1
    endif

    do i=1,lx

        if (i.gt.1) then
            in = i-1
        else 
            in = LX 
        endif
        if (i.lt.LX) then 
            ip = i+1
        else 
            ip = 1
        endif
  •      neighbor 2
            psi_x = (1./9. )*ex(2)*psi(j,ip)
            psi_y = (1./9. )*ey(2)*psi(j,ip)
    
  •      neighbor 3
            psi_x = psi_x + (1./9. )*ex(3)*psi(jp,i)
            psi_y = psi_y + (1./9. )*ey(3)*psi(jp,i)
    
  •      neighbor 4
            psi_x = psi_x + (1./9. )*ex(4)*psi(j,in)
            psi_y = psi_y + (1./9. )*ey(4)*psi(j,in)
    
  •     neighbor 5
            psi_x = psi_x + (1./9. )*ex(5)*psi(jn,i)
            psi_y = psi_y + (1./9. )*ey(5)*psi(jn,i)
    
  •     neighbor 6
            psi_x = psi_x + (1./36. )*ex(6)*psi(jp,ip)
            psi_y = psi_y + (1./36. )*ey(6)*psi(jp,ip)
    
  •     neighbor 7
            psi_x = psi_x + (1./36. )*ex(7)*psi(jp,in)
            psi_y = psi_y + (1./36. )*ey(7)*psi(jp,in)
    
  •     neighbor 8
            psi_x = psi_x + (1./36. )*ex(8)*psi(jn,in)
            psi_y = psi_y + (1./36. )*ey(8)*psi(jn,in)
    
  •     neighbor 9
            psi_x = psi_x + (1./36. )*ex(9)*psi(jn,ip)
            psi_y = psi_y + (1./36. )*ey(9)*psi(jn,ip)
    
         forcex(j,i) = -120.*psi(j,i)*psi_x
         forcey(j,i) = -120.*psi(j,i)*psi_y
      enddo
    enddo
    

C Compute the equilibrium distribution function, feq.

  f1=3.
  f2=9./2.
  f3=3./2.

  do j=1,ly   
    
    do i=1,lx
               
        if(is_solid_node(j,i).eq.0) then
            

            rt0 = (4./9. )*rho(j,i)
            rt1 = (1./9. )*rho(j,i)
            rt2 = (1./36.)*rho(j,i)

            ueqxij =  u_x(j,i)+tau*forcex(j,i)/rho(j,i)
            ueqyij =  u_y(j,i)+tau*forcey(j,i)/rho(j,i)
            uxsq   =  ueqxij * ueqxij
            uysq   =  ueqyij * ueqyij
            uxuy5  =  ueqxij +  ueqyij
            uxuy6  = -ueqxij +  ueqyij
            uxuy7  = -ueqxij -ueqyij
            uxuy8  =  ueqxij -ueqyij
            usq    =  uxsq + uysq

            
            feq(j,i,0+1) = rt0*(1.                      - f3*usq)
            feq(j,i,1+1) = rt1*(1.+ f1*ueqxij+f2*uxsq - f3*usq)
            feq(j,i,2+1) = rt1*(1.+ f1*ueqyij+f2*uysq - f3*usq)
            feq(j,i,3+1) = rt1*(1.- f1*ueqxij+f2*uxsq - f3*usq)
            feq(j,i,4+1) = rt1*(1.- f1*ueqyij+f2*uysq - f3*usq)
            feq(j,i,5+1) = rt2*(1.+ f1*uxuy5 +f2*uxuy5*uxuy5-f3*usq)
            feq(j,i,6+1) = rt2*(1.+ f1*uxuy6 +f2*uxuy6*uxuy6-f3*usq)
            feq(j,i,7+1) = rt2*(1.+ f1*uxuy7 +f2*uxuy7*uxuy7-f3*usq)
            feq(j,i,8+1) = rt2*(1.+ f1*uxuy8 +f2*uxuy8*uxuy8-f3*usq)

            
        endif
     enddo
     
  enddo

C Collision step.
do j=1,ly
do i=1,lx

        if(is_solid_node(j,i).eq.1) then               

C Standard bounceback

            temp= f(j,i,1+1)
            f(j,i,1+1) = f(j,i,3+1)
            f(j,i,3+1) = temp
            temp= f(j,i,2+1)
            f(j,i,2+1) = f(j,i,4+1)
            f(j,i,4+1) = temp
            temp= f(j,i,5+1)
            f(j,i,5+1) = f(j,i,7+1)
            f(j,i,7+1) = temp
            temp= f(j,i,6+1)
            f(j,i,6+1) = f(j,i,8+1)
            f(j,i,8+1) = temp
            
        else  

C Regular collision

            do k=1,9
                               
                f(j,i,k) = f(j,i,k)-( f(j,i,k) - feq(j,i,k))/tau
                
            enddo 
            
        endif
    enddo
    
  enddo

C Streaming step; subtle changes to periodicity here due to indexing

  do j=1,ly

    if (j.gt.1) then 
        jn = j-1
    else
        jn = LY
    endif
    
    if (j.lt.ly) then
        jp = j+1
    else 
        jp = 1
    endif

    do i=1,lx

        if (i.gt.1) then
            in = i-1
        else 
            in = LX 
        endif
        if (i.lt.LX) then 
            ip = i+1
        else 
            ip = 1
        endif

        ftemp(j,i,0+1)  = f(j,i,0+1)
        ftemp(j,ip,1+1) = f(j,i,1+1)
        ftemp(jp,i,2+1) = f(j,i,2+1)
        ftemp(j,in,3+1) = f(j,i,3+1)
        ftemp(jn,i ,4+1) = f(j,i,4+1)
        ftemp(jp,ip,5+1) = f(j,i,5+1)
        ftemp(jp,in,6+1) = f(j,i,6+1)
        ftemp(jn,in,7+1) = f(j,i,7+1)
        ftemp(jn,ip,8+1) = f(j,i,8+1)

  enddo
  enddo

  f=ftemp

C End time loop
enddo

  open(unit=20,file='rho.dat',status='unknown')

  do j=1,ly
  do i=1,lx
  
    write(20,*) rho(j,i)
  enddo
  enddo


  end