Dear friends
I write a program to solve free energy two phase (wetting and spreading), but i am facing problem how to include external force in my source code (fortran).
I got satisfied results for simulating droplet at various contact angles, but when I applied gravitational force, something went wrong.
I paste here the simulation algorithm and subroutine for external force. Can anybody advice me?
thank you
DO ITER = 1, 200000
CALL COLLIDE(lx,ly,cd,f,feq,fnew,TAU)
CALL MOVE(lx,ly,cd,f,fnew)
CALL BOUNDARY(lx,ly,cd,f)
CALL MACRO(lx,ly,cd,rho,f,u,v)
CALL EQUILIBRIUM(LX,LY,CD,U,V,CX,CY,RHO,KAPPA,NU,A,B,T,FEQ)
IF(MOD(ITER,100) .EQ.0)THEN
WRITE(*,*)ITER
CALL fileoutput(lx,ly,u,v,rho,ITER)
END IF
IF(ITER .GT. 2000) THEN
CALL FORCE(LX,LY,CD,NG,CX,CY,F,RHO,FEQ,V,NL)
END IF
END DO
SUBROUTINE FORCE(LX,LY,CD,NG,CX,CY,F,RHO,FEQ,V,NL)
REAL*8 F(1:LX,1:LY,0:CD),FEQ(1:LX,1:LY,0:CD)
REAL*8 RHO(1:LX,1:LY),V(1:LX,1:LY)
REAL*8 CX(0:CD),CY(0:CD)
REAL*8 G,W,NG,WI,NL
G= 0.00012D0
DO K =0,CD
DO I = 1,LX
DO J = 1,LY
IF (K .EQ. 0) WI = 4.0D0/9.0D0
IF (K .GE. 1 .LE. 4) WI = 1.0D0/9.0D0
IF (K .GE. 5 .LE. 8) WI = 1.0D0/36.0D0
F(I,J,K) = F(I,J,K) - 3.0D0*(CY(K))*WI*(1.0D0-NG/RHO(I,J))*G
END DO
END DO
END DO
RETURN
END