Problem with C-H +Stokes + viscosity coupling

Hi! I’m trying to simulate the mixing of two liquids of different viscosities using Cahn-Hilliard and Stokes equations. Due to the non-linear nature of C-H equation I’m using Newton iterations. If the viscosity ratio (lambda) is 1 the code works without any problems but for any other ratio after around 100 iterations I start getting these minimums and maximums that spread and take over the entire domain as shown in the picture. Any ideas what’s wrong with my code? Thanks!


real radius = 0.05;
real M=1;
real centerX = 0.35;
real centerX2 = 0.5;
real centerY = 0.5;
real Inside = 1.0;
real Outside = -1.0;
real Pe= 10000;
real Ch=0.01;

func real initialCondition(real x, real y)
{
  if (y<0.5)
    return Inside;
  else
    return Outside;
}

mesh Th=square(72,72);
plot(Th,wait=1);
fespace Vh(Th,P1);
fespace Vh2(Th,[P1,P1]);

real lambda=10;
real a=0.01;
real dt=0.01;
int i,k;  

func real df(real u) { return u^3-u; }
func real ddf(real u) { return 3*u^2-1; }

Vh2 [u,w],[v,z],[oldu,oldz],[auxv,auxz],[vv,zz];
Vh dfalpha;  //  f'(u)
Vh ddfalpha; // f"(u)

fespace Uh(Th, P1);
Uh U, V;
Uh UU, VV;

fespace Ph(Th, P1);
Ph p, pp;

varf vdJ([v,z],[phi,psi]) = int2d(Th)
((u-oldu)*phi+dt*(dx(u)*U+dy(u)*V)*phi-dt*1/Pe*(dx(oldz)*dx(phi)+dy(oldz)*dy(phi))+oldz*psi- 
Ch*Ch*(dx(oldu)*dx(psi)+dy(oldu)*dy(psi))-dfalpha*psi);

varf vhJ([v,z],[phi,psi]) = int2d(Th)
(-v*phi-dt*1/Pe*(dx(z)*dx(phi)+dy(z)*dy(phi))+z*psi-Ch*Ch*(dx(v)*dx(psi) 
+dy(v)*dy(psi))-ddfalpha*v*psi);

[u,w]=[initialCondition(x,y),0];


plot (u,wait=1,cmm="Cahn-Hilliard",value=true);

macro SGrad(u,v) [[dx(u),0.5*(dx(v)+dy(u))],[0.5*(dx(v)+dy(u)),dy(v)]] //
macro div(u1,u2) (dx(u1)+dy(u2)) //
macro mu(a) (0.5*(a+1)-1/lambda*0.5*(a-1)) //

    

    for (int i=0;i<150;i++)
      {
    oldu[ ]=u[ ];
    solve Stokes( [U,V,p],[UU,VV,pp]) =
      int2d(Th)( 2* mu(oldu)*(SGrad(U,V):SGrad(UU,VV)) - div(U,V)*pp)
    - int2d(Th)( div(UU,VV)*p)
    + on(3, 2, 4, U=0, V=0)
    + on(1, U=16*(x)^2*(1-x)^2, V=0)
    ;
       for (int k=0;k<10;k++) // boucle de Newton
	{
	dfalpha = df( oldu ) ;
	ddfalpha = ddf( oldu );
	auxv[]=vdJ(0,Vh2);
          real res= auxv[]'*auxv[];
          cout << i <<  " residu^2 = " <<  res  << endl;
          matrix H;
          if( res< 1e-12) break;
          H= vhJ(Vh2,Vh2,factorize=1,solver=LU);
          vv[]=H^-1*auxv[];
          v[] -= vv[];
          oldu[]=v[];
	}
      u[]=v[];

      plot (u,wait=1,fill=true, cmm="Cahn-Hilliard"+i,value=true);

      }

Your Peclet number is huge. Try with a smaller Pe value or dramatically increase the resolution of your mesh. Also, you should use an inf-sup stable finite element discretization like P2-P1.

Thank you so much for the reply! I’ll try out those suggestions. If I can ask one more thing, did I add the advection correctly considering I’m using Newton iterations? I don’t have much experience with it so I’m not sure if I should add anything to vhJ?

EDIT: P2-P1 discretization and bigger mesh resolution fixed the problem. Thank you!

1 Like