Basilisk source code (http://basilisk.fr/src/)

root / src / tension.h.page

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
/**
# Surface tension

Surface tension can be expressed as the interfacial force density
$$
\phi\nabla f
$$
with $f$ the volume fraction describing the interface and the potential
$$
\phi = \sigma\kappa
$$
with $\sigma$ the (constant) surface tension coefficient and $\kappa$
the interface mean curvature. */

#include "iforce.h"
#include "curvature.h"

/**
The surface tension coefficient is associated to each VOF tracer. */

attribute {
  double sigma;
}

/**
## Stability condition

The surface tension scheme is time-explicit so the maximum timestep is
the oscillation period of the smallest capillary wave.
$$
T = \sqrt{\frac{\rho_{m}\Delta_{min}^3}{\pi\sigma}}
$$
with $\rho_m=(\rho_1+\rho_2)/2.$ and $\rho_1$, $\rho_2$ the densities
on either side of the interface. */

event stability (i++) {

  /**
  We first compute the minimum and maximum values of $\alpha/f_m =
  1/\rho$, as well as $\Delta_{min}$. */

  double amin = HUGE, amax = -HUGE, dmin = HUGE;
  foreach_face (reduction(min:amin) reduction(max:amax) reduction(min:dmin)) {
    if (alpha.x[]/fm.x[] > amax) amax = alpha.x[]/fm.x[];
    if (alpha.x[]/fm.x[] < amin) amin = alpha.x[]/fm.x[];
    if (Delta < dmin) dmin = Delta;
  }
  double rhom = (1./amin + 1./amax)/2.;

  /**
  We then consider each VOF interface with an associated value of
  $\sigma$ different from zero and set the maximum timestep. */
  
  for (scalar c in interfaces)
    if (c.sigma) {
      double dt = sqrt (rhom*cube(dmin)/(pi*c.sigma));
      if (dt < dtmax)
	dtmax = dt;
    }
}

/**
## Definition of the potential

We overload the acceleration event to define the potential
$\phi=\sigma\kappa$. */

event acceleration (i++)
{
  
  /**
  We check for all VOF interfaces for which $\sigma$ is non-zero. */

  for (scalar f in interfaces)
    if (f.sigma) {
      
      /**
      If $\phi$ is already allocated, we add $\sigma\kappa$, otherwise
      we allocate a new field and set it to $\sigma\kappa$. */

      scalar phi = f.phi;
      if (phi.i)
	curvature (f, phi, f.sigma, add = true);
      else {
	phi = new scalar;
	curvature (f, phi, f.sigma, add = false);
	f.phi = phi;
      }
    }
}