hub.darcs.net :: basilisk -> basilisk -> files

Basilisk source code

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