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

Basilisk source code

## root / src / two-phase.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 91 92 93 94 95 96 97 98 99 100 101 102 103  /** # Two-phase interfacial flows This file helps setup simulations for flows of two fluids separated by an interface (i.e. immiscible fluids). It is typically used in combination with a [Navier--Stokes solver](navier-stokes/centered.h). The interface between the fluids is tracked with a Volume-Of-Fluid method. The volume fraction in fluid 1 is $f=1$ and $f=0$ in fluid 2. The densities and dynamic viscosities for fluid 1 and 2 are *rho1*, *mu1*, *rho2*, *mu2*, respectively. */ #include "vof.h" scalar f[], * interfaces = {f}; double rho1 = 1., mu1 = 0., rho2 = 1., mu2 = 0.; /** Auxilliary fields are necessary to define the (variable) specific volume $\alpha=1/\rho$ as well as the cell-centered density. */ face vector alphav[]; scalar rhov[]; event defaults (i = 0) { alpha = alphav; rho = rhov; /** If the viscosity is non-zero, we need to allocate the face-centered viscosity field. */ if (mu1 || mu2) mu = new face vector; } /** The density and viscosity are defined using arithmetic averages by default. The user can overload these definitions to use other types of averages (i.e. harmonic). */ #ifndef rho # define rho(f) (clamp(f,0.,1.)*(rho1 - rho2) + rho2) #endif #ifndef mu # define mu(f) (clamp(f,0.,1.)*(mu1 - mu2) + mu2) #endif /** We have the option of using some "smearing" of the density/viscosity jump. */ #ifdef FILTERED scalar sf[]; #else # define sf f #endif event properties (i++) { /** When using smearing of the density jump, we initialise *sf* with the vertex-average of *f*. */ #ifndef sf #if dimension <= 2 foreach() sf[] = (4.*f[] + 2.*(f[0,1] + f[0,-1] + f[1,0] + f[-1,0]) + f[-1,-1] + f[1,-1] + f[1,1] + f[-1,1])/16.; #else // dimension == 3 foreach() sf[] = (8.*f[] + 4.*(f[-1] + f[1] + f[0,1] + f[0,-1] + f[0,0,1] + f[0,0,-1]) + 2.*(f[-1,1] + f[-1,0,1] + f[-1,0,-1] + f[-1,-1] + f[0,1,1] + f[0,1,-1] + f[0,-1,1] + f[0,-1,-1] + f[1,1] + f[1,0,1] + f[1,-1] + f[1,0,-1]) + f[1,-1,1] + f[-1,1,1] + f[-1,1,-1] + f[1,1,1] + f[1,1,-1] + f[-1,-1,-1] + f[1,-1,-1] + f[-1,-1,1])/64.; #endif #endif #if TREE sf.prolongation = refine_bilinear; boundary ({sf}); #endif foreach_face() { double ff = (sf[] + sf[-1])/2.; alphav.x[] = fm.x[]/rho(ff); if (mu1 || mu2) { face vector muv = mu; muv.x[] = fm.x[]*mu(ff); } } foreach() rhov[] = cm[]*rho(sf[]); #if TREE sf.prolongation = fraction_refine; boundary ({sf}); #endif }