Different results for droplet impact simulation: domain boundary vs. embed surface

[Ivan Vozhakov]

Dear colleagues,

I am facing an issue while simulating a droplet impact on a flat surface in Basilisk. The core of the problem is as follows:

• Case 1: The surface is the bottom face of the computational domain (i.e., a standard boundary condition).

• Case 2: The surface is defined using the embed boundaries method.

The results for these two cases are significantly different (e.g., in the spreading shape, rebound time, etc.).

In this regard, I have two main questions:

1. Can such a large discrepancy be a known limitation of the embed method for this kind of configuration?

2. Or am I most likely configuring something incorrectly in the embed case?

I would be very grateful for any advice or comments on this matter.

[Wojciech Aniszewski]

Hi Ivan,
I tried to reproduce these curves but I'm travelling and these 3D simulations are too battery-costly if you know what I mean.
But here are some thoughts I had straight off the bat looking at your source files.

- I'd try to lift the f[bottom] = neumann(0.);
- I'd try to see how sensitive the D_max writing procedure is to the condition here:
,----
| - if (y > -L0/2 && y < -L0/2 + 3*Delta && f[] >= 0.5 && f[] <= 1 && x > x_max)
| + if (y > y_min && y < y_min + 3*Delta && f[] >= 0.5 && f[] <= 1 && x > x_max) ,
`----
namely to the values of vof. Perhaps a '&& f[]>=0' condition would yield a different result... same goes for your y_min computation in the precursor (cs[] <= 0.9) of the embed simulation;

-for the embed simulation, I'd agressively refine on the cs[] (so have adapt_wavelet{{f,u,cs}...) or whatamacallit the solid array );
- look up fractions_cleanup() as well.

Now, there could also be an entire discussion here about the contact angle but I'll just blatantly ignore that and finish here.
cheers
Voitek

On Wed, Sep 24, 2025 at 01:51:00AM -0700, Ivan Vozhakov wrote:
>
>
> Dear colleagues,
>
> I am facing an issue while simulating a droplet impact on a flat surface in
> Basilisk. The core of the problem is as follows:
>

> -

>
> Case 1: The surface is the bottom face of the computational domain
> (i.e., a standard boundary condition).

> -

>
> Case 2: The surface is defined using the embed boundaries method.
>
> The results for these two cases are significantly different (e.g., in the
> spreading shape, rebound time, etc.).
>
> In this regard, I have two main questions:
>

> 1.

>
> Can such a large discrepancy be a known limitation of the embed method
> for this kind of configuration?

> 2.

>
> Or am I most likely configuring something incorrectly in the embed case?
>
> I would be very grateful for any advice or comments on this matter.
>

> --
> 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/375613a9-f6a8-4686-98a0-17f9e7ee5faan%40googlegroups.com.


> #include "grid/octree.h"
> #include "navier-stokes/centered.h"
> #include "two-phase.h"
> #include "tension.h"
> #include "navier-stokes/perfs.h"
> #include "view.h"
>
> #define scale 0.001
>
> #define We 10
> #define D0 2*scale
> #define RHO1 1000
> #define RHO2 1.2
> #define MU1 1.47e-3
> #define MU2 1.8e-5
> #define SIGMA 0.072
> #define U_down sqrt(We * SIGMA / (RHO1 * D0))
>
> int MAXLEV = 8;
> int MINLEV = 8;
>
> double tEnd = 0.06;
>
> u.n[top] = neumann(0.);
> u.t[top] = neumann(0.);
> u.r[top] = neumann(0.);
> f[top] = neumann(0.);
> p[top] = dirichlet(0.);
>
> u.n[bottom] = dirichlet(0.);
> u.t[bottom] = dirichlet(0.);
> u.r[bottom] = dirichlet(0.);
> f[bottom] = neumann(0.);
>
> double shift = 2*scale;
>
> event init (t = 0) {
> printf("run event: init \n");
> if (!restore (file = "restart")) {
> fraction(f, -sq(x)-sq(y+shift)-sq(z) + sq(D0 / 2));
> printf("%f\n", U_down);
> foreach() {
> u.y[] = -U_down*f[];
> }
> boundary(all);
> }
> else {
> printf("restarting ! \n");
> }
> printf("finished event: init \n");
> }
>
> int poweroftwo(int p)
> {
> int res = 1;
> for (int i=0;i<p;i++) {
> res*=2;
> }
> return res;
> }
>
> int main()
> {
> init_grid (poweroftwo(MINLEV));
> size (10*scale);
> origin (-L0/2, -L0/2, -L0/2);
>
> rho1 = RHO1, rho2 = RHO2;
> mu1 = MU1, mu2 = MU2;
> f.sigma = SIGMA;
>
> printf("run \n");
> run();
> }
>
> event adapt (i++) {
> double uemax = 0.01;
> astats s = adapt_wavelet ({f,u}, (double[]){0.0001, uemax, uemax, uemax}, maxlevel = MAXLEV, minlevel = MINLEV);
>
> fprintf (stderr, "# refined %d cells, coarsened %d cells\n", s.nf, s.nc);
> }
>
> event logfile (i++)
> fprintf (stderr, "%d %g %d %d\n", i, t, mgp.i, mgu.i);
>
> event maxDiameter(t = 0;t += 0.00001; t <= tEnd) {
>
> double x_max = 0;
>
> foreach(reduction(max:x_max)) {
> if (y > -L0/2 && y < -L0/2 + 3*Delta && f[] >= 0.5 && f[] <= 1 && x > x_max)
> x_max = x;
> }
>
> int id = 0;
>
> #if _MPI
> MPI_Comm_rank(MPI_COMM_WORLD, &id);
> #endif
>
> if(id == 0) {
> FILE *parall_contr;
>
> parall_contr = fopen("D(t).dat", "a");
>
> fprintf(parall_contr, "%f %f\n", t, 2 * x_max);
>
> fclose(parall_contr);
> }
> }
>
> event movie (t = 0; t += 5.e-5; t <= tEnd)
> {
> clear();
> draw_vof ("f");
> box();
> save ("movie.mp4");
> }

> #include "grid/octree.h"
> #include "embed.h"
> #include "navier-stokes/centered.h"
> #include "two-phase.h"
> #include "tension.h"
> #include "navier-stokes/perfs.h"
> #include "view.h"
>
> #define scale 0.001
>
> #define We 10
> #define D0 2*scale
> #define RHO1 1000
> #define RHO2 1.2
> #define MU1 1.47e-3
> #define MU2 1.8e-5
> #define SIGMA 0.072
> #define U_down sqrt(We * SIGMA / (RHO1 * D0))
>
> int MAXLEV = 8;
> int MINLEV = 8;
>
> double tEnd = 0.06;
>
> u.n[top] = neumann(0.);
> u.t[top] = neumann(0.);
> u.r[top] = neumann(0.);
> f[top] = neumann(0.);
> p[top] = dirichlet(0.);
>
> u.n[embed] = dirichlet(0.);
> u.t[embed] = dirichlet(0.);
> u.r[embed] = dirichlet(0.);
> f[embed] = neumann(0.);
>
> double shift = 2*scale;
>
> event init (t = 0) {
> printf("run event: init \n");
> if (!restore (file = "restart")) {
> fraction(f, -sq(x)-sq(y)-sq(z) + sq(D0 / 2));
> solid (cs, fs, y + shift);
> printf("%f\n", U_down);
> foreach() {
> u.y[] = -U_down*f[];
> }
> boundary(all);
> }
> else {
> printf("restarting ! \n");
> }
> printf("finished event: init \n");
> }
>
> int poweroftwo(int p)
> {
> int res = 1;
> for (int i=0;i<p;i++) {
> res*=2;
> }
> return res;
> }
>
> int main()
> {
> init_grid (poweroftwo(MINLEV));
> size (10*scale);
> origin (-L0/2, -L0/2, -L0/2);
>
> rho1 = RHO1, rho2 = RHO2;
> mu1 = MU1, mu2 = MU2;
> f.sigma = SIGMA;
>
> printf("run \n");
> run();
> }
>
> event adapt (i++) {
> double uemax = 0.01;
> astats s = adapt_wavelet ({f,u}, (double[]){0.0001, uemax, uemax, uemax}, maxlevel = MAXLEV, minlevel = MINLEV);
>
> fprintf (stderr, "# refined %d cells, coarsened %d cells\n", s.nf, s.nc);
> }
>
> event logfile (i++)
> fprintf (stderr, "%d %g %d %d\n", i, t, mgp.i, mgu.i);
>
> event maxDiameter(t = 0;t += 0.00001; t <= tEnd) {
> double y_min = -1;
>
> foreach(reduction(max:y_min)) {
> if (y > y_min && cs[] <= 0.9)
> y_min = y;
> }
>
> double x_max = 0;
>
> foreach(reduction(max:x_max)) {
> if (y > y_min && y < y_min + 3*Delta && f[] >= 0.5 && f[] <= 1 && x > x_max)
> x_max = x;
> }
>
> int id = 0;
>
> #if _MPI
> MPI_Comm_rank(MPI_COMM_WORLD, &id);
> #endif
>
> if(id == 0) {
> FILE *parall_contr;
>
> parall_contr = fopen("D(t).dat", "a");
>
> fprintf(parall_contr, "%f %f\n", t, 2 * x_max);
>
> fclose(parall_contr);
> }
> }
>
> event movie (t = 0; t += 5.e-5; t <= tEnd)
> {
> clear();
> draw_vof ("f_rep");
> draw_vof("cs");
> box();
> save ("movie.mp4");
> }


--
/^..^\
( (••) ) Wojciech (Vôitek) ANISZEWSKI
(|)_._(|)~ Research Assistant @ Univ-Rouen
GPG ID : AC66485E
Twitter : @echo_dancers3
Mastodon: @w...@fediscience.org
BlueSky : @aniszewski.bsky.social
Scholar : https://tinyurl.com/y28b8gfp
OrcId : https://orcid.org/0000-0002-4248-1194
RG : https://www.researchgate.net/profile/Wojciech_Aniszewski

[Ivan Vozhakov]

Hi Voitek!

Thank you very much for your detailed reply! I sincerely appreciate that you took the time to analyze my question despite your busy schedule.

I completely agree with your remark regarding the conditions for determining the diameter. You are right, they certainly can affect the accuracy. However, I also expected that their variation should not exceed several cells, so the problem apparently goes deeper.

Thank you for the advice on adaptation and the fractions_cleanup() function. I am currently studying how to apply it correctly in my case.

You also mentioned the contact angle, which I agree is probably key here. I've noticed that even with f[embed]=neumann(0) (which should give a 90-degree contact angle), the contact angle looks much larger. It's clearly visible in the simulation (I've attached a video).

Best,
Ivan

четверг, 25 сентября 2025 г. в 20:27:29 UTC+7, Wojciech Aniszewski:

[Spencer Schwartz]

Dear Ivan,

I believe doing

f[embed] = neumann(0)

 

doesn't actually do anything, which is why it looks like the contact angle is much higher. Instead, you should try using https://basilisk.fr/sandbox/popinet/contact/contact-embed.h to set the contact angle w/embed.

You can also try doing f[bottom] = dirichlet(0) for the case w/o embed, like how Voitek recommended, since embed has a non-wetting condition w/o contact-embed.h.

Best,

Spencer