Regarding contact angle in Allen-Cahn phase field simulation

[Shyam Sunder Yadav]

Dear All

I have implemented a second order Runge-Kutta method based Basilisk formulation for conservative Allen-Cahn  equation.

I am simulating the equilibrium shape of a drop placed on bottom boundary under gravity.

I am implementing a contact angle to the drop at the bottom boundary but the angle is not enforced.

I am attaching the code, can someone please check the implementation in the phase_field_RK2_chem_pot.h file and let me know what I am missing.

The contact angle information is from a paper which is inside the zipped file.

Thank you. Regards

--

Dr. Shyam Sunder Yadav
Associate Professor
Mechanical Engineering
BITS Pilani
09902346342
http://www.bits-pilani.ac.in/pilani/ssyadav/Profile

The information contained in this electronic communication is intended solely for the individual(s) or entity to which it is addressed. It may contain proprietary, confidential and/or legally privileged information. Any review, retransmission, dissemination, printing, copying or other use of, or taking any action in reliance on the contents of this information by person(s) or entities other than the intended recipient is strictly prohibited and may be unlawful. If you have received this communication in error, please notify us by responding to this email or telephone and immediately and permanently delete all copies of this message and any attachments from your system(s). The contents of this message do not necessarily represent the views or policies of BITS Pilani.

[Shyam Sunder Yadav]

Dear All

I observed that if I implement the contact angle using the lines below then the angle is not enforced:

    foreach_boundary(bottom)
    {
     phi[0,-1] = phi[] + sqrt(2.0*my_beta/my_kappa)*cos(mytheta)*phi[]*(1.0-phi[])*Delta; 
    } 

Whereas if I use the below code, the angle is enforced (we are modifying the phi values in the interior cells adjacent to the bottom boundary!!!):

    foreach_boundary(bottom)
    {
     phi[] = phi[0,1] + sqrt(2.0*my_beta/my_kappa)*cos(mytheta)*phi[0,1]*(1.0-phi[0,1])*Delta; 

    } 

Can someone please figure out what's wrong with the correct, first approach?

Thanks

[Yuanpeng Zhang (Yannick)]

Dear  Dr. Shyam Sunder Yadav,

I think the idea of your first approach is correct, but the way it is implemented in Basilisk is incomplete.

In Basilisk, manually assigning a ghost value with foreach_boundary() does not register a boundary condition for the scalar. Once phi[] is modified, the field is marked as dirty. If a later stencil requires neighboring values, Basilisk may automatically call boundary_internal() and reconstruct the ghost cells using the boundary condition registered for phi. If no wetting condition has been registered, the default homogeneous Neumann condition may therefore overwrite the ghost value you assigned manually.

So I would suggest registering the nonlinear Neumann condition directly, for both the main field and the RK intermediate field:

 #define wetting_Lambda  \ 

(-sqrt(2.0*my_beta/my_kappa)*cos(mytheta)) 

 #define wetting_gradient(q) \ 

 (wetting_Lambda*(q)*(1.0 - (q))) 

 phi[bottom] = neumann (wetting_gradient(phi[])); 

phi_1[bottom] = neumann (wetting_gradient(phi_1[]));  

Then Basilisk's automatic boundary updates will regenerate the required wetting ghost values consistently.

Your second approach appears to work because it modifies the first interior cell rather than the ghost cell. Automatic ghost-cell updates cannot undo this modification. However, this is no longer the original Neumann boundary condition: it overwrites the interior PDE solution and may also spoil conservation.

If higher accuracy at the wall is required, one can instead use the wall value \phi_w \simeq \frac{\phi_P+\phi_G}{2} in the nonlinear wetting condition and solve locally for the ghost value, but I would first test the registered Neumann implementation above.

Best regards,

Peng

[Shyam Sunder Yadav]

Dear Peng

Thank you very much for the response and nice explanation...

I tried your suggestion but it is still not working...

The Basilisk in my system may be corrupted,  can you please try in your system?

Thanks again dear...

Regards

[Yuanpeng Zhang (Yannick)]

Dear  Dr. Shyam Sunder Yadav,

Thank you for your reply.

I tested the case on my system. I have attached the code together with the corresponding results in a ZIP file for your reference.

Best regards,
Peng

[Shyam Sunder Yadav]

Dear Yuanpeng

Thanks for the efforts but the contact angle is still not enforced...

Please check again the images you have sent...

The drop spreads under gravity while maintaining an angle of 90 degree which is due to the homogeneous  Neumann condition taken by the code...

Don't know why it's behaving like this...

Please check again.

Regards  

--
You received this message because you are subscribed to the Google Groups "basilisk-fr" group.
To unsubscribe from this group and stop receiving emails from it, send an email to basilisk-fr...@googlegroups.com.
To view this discussion visit https://groups.google.com/d/msgid/basilisk-fr/04d991b0-4cb9-4c7b-a699-bbbe98907e18n%40googlegroups.com.

[Yuanpeng Zhang (Yannick)]

Dear   Dr. Shyam Sunder Yadav,

I did a simple test, and the contact-angle boundary condition seems to work. However, whether this treatment is physically correct still needs to be checked.

The problem is not caused by the homogeneous Neumann condition. It comes from the definition

#define EPSILON 3.0*L0/256.0

which should be changed to

#define EPSILON (3.0*L0/256.0)

Without the parentheses, expressions such as 4.0/EPSILON are expanded incorrectly.

Best regards,
Peng

[Shyam Sunder Yadav]

Dear Yannick

Thank you very much!!!

It takes the correct angle now!!!

I never thought that things can go wrong 

this way...

Thanks again for the efforts...

How did you figure it out?

Regards 

To view this discussion visit https://groups.google.com/d/msgid/basilisk-fr/4e1f9180-7f17-4e34-8362-a958c0c65a85n%40googlegroups.com.