Problem with C-H +Stokes + viscosity coupling

General Discussion
1

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!

2023-08-02 17_23_11-Window
2023-08-02 17_23_11-Window945×658 28.1 KB


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);

      }
2

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.

3

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
4

//Mesh
mesh Th=square(100,100);
//border C(t=0,2pi){x=cos(t);y=sin(t);};
//mesh Th=buildmesh (C(50));
plot(Th,wait=1);
macro Grad(u)[dx(u),dy(u)]//
macro uex(t)(exp(-2t)(sin(pix)^2)(sin(piy)^2)) //

//Fespace
fespace Vh(Th,P1);
fespace Vh2(Th,[P1,P1]);
//Vh uex;

//Parameter
real a=0.01;
real dt=0.10;
int i,k;
//real t;

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

Vh2 [u,w],[v,z],[oldv,oldz],[auxv,auxz],[vv,zz];
Vh dfalpha;
Vh ddfalpha;
//Vh uex;

//macro uex(t) (exp(-2t)(sin(pix)^2)(sin(pi*y)^2)//
// Problem
varf vdJ([v,z],[phi,psi])=int2d(Th)
((oldv-u)phi+dt(dx(oldz)dx(phi)+dy(oldz)dy(phi))+oldzpsi-
a(dx(oldv)*dx(psi)+dy(oldv)dy(psi))-dfalphapsi);

varf vhJ([v,z],[phi,psi])=int2d(Th)
(vphi+dt(dx(z)dx(phi)+dy(z)dy(phi))+zpsi-a(dx(v)dx(psi)
+dy(v)dy(psi))-ddfalphavpsi);

//u=uex(t);

[u,w]=[(sin(pix)^2)(sin(piy)^2),-2a*(pi^2)(sin(pi^2(x^2+y^2))
-2*(pi^2)(sin(pix)^2)(sin(piy)^2)(x^2+y^2))+((sin(pix)^2)^3)((sin(piy)^2)^3)
-(sin(pix)^2)(sin(pi*y)^2)];
//u=uex;

plot(u,wait=1,cmm=“Cahn-Hilliard” ,value=true);
//u=uex(t);
//macro uex(t) (exp(-2t)(sin(pix)^2)(sin(piy)^2))//
// u=uex(t);
//uex(t)=exp(-2t)(sin(pix)^2)(sin(piy)^2);
//real t=0;// start from t=0
real t;
//t=0.1;
for (int i=0;i<1./dt; i++) // Time loop
{
t=t+dt;
//u=uex(t);
//=exp(-2t)((sin(pix)^2)(sin(piy)^2));
{
oldv[]=u[];
for (int k=0;k<10;k++)// Newton loop
{
dfalpha=df( oldv );
ddfalpha=ddf( oldv );
auxv[]=vdJ(0,Vh2);
real res= auxv[]'auxv[]; // Residual
cout << i << " residu^2 = " << res << endl;
matrix H;
if (res< 1e-12) break;
H=vhJ(Vh2,Vh2, factorize=1,solver=LU);
vv[]=H^-1auxv[]; // Newton
v[] -=vv[];
oldv[]=v[];
}
u[]=v[];
//u=uex(t);
//uex=expsin(pix)^2sin(piy)^2;
plot(u,wait=1,cmm=“Cahn-Hilliard”,value=true);
//uex(t)=exp(-2t)sin(pix)^2sin(piy)^2;
//plot(uex(t),wait=1,cmm=“Cahn-Hilliard exact solution”,value=true);

}
 cout << "t" << t << "L^2-error=" << sqrt(int2d(Th)((u-uex(t))^2))<<endl;
 //cout << "t"<< t << "L^2-error=" << sqrt(int2d(Th)
 //(abs(u-(exp(-2*t)*(sin(pi*x)^2)*(sin(pi*y)^2)))^2))<<endl;
}

Can you tell me Please what is wrong in this code. Boundary of the approximate solution not matching with exact solution. Please help. It is very urgent.

5

Can you tell me how i post a queries in this platform?