# need help on free energy multiphase

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``````

Hello,

You probably need to implement the proper forcing as in the paper Guo et al. in their paper. Also introduction of gravity sometimes can cause some additional condensation or evaporation which can make scheme unstable.