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

Basilisk source code

## root / src / iforce.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 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150  /** # Interfacial forces We assume that the interfacial acceleration can be expressed as $$\phi\mathbf{n}\delta_s/\rho$$ with $\mathbf{n}$ the interface normal, $\delta_s$ the interface Dirac function, $\rho$ the density and $\phi$ a generic scalar field. Using a CSF/Peskin-like approximation, this can be expressed as $$\phi\nabla f/\rho$$ with $f$ the volume fraction field describing the interface. The interfacial force potential $\phi$ is associated to each VOF tracer. This is done easily by adding the following [field attributes](/Basilisk C#field-attributes). */ attribute { scalar phi; } /** Interfacial forces are a source term in the right-hand-side of the evolution equation for the velocity of the [centered Navier--Stokes solver](navier-stokes/centered.h) i.e. it is an acceleration. If necessary, we allocate a new vector field to store it. */ event defaults (i = 0) { if (is_constant(a.x)) { a = new face vector; foreach_face() a.x[] = 0.; boundary ((scalar *){a}); } } /** The calculation of the acceleration is done by this event, overloaded from [its definition](navier-stokes/centered.h#acceleration-term) in the centered Navier--Stokes solver. */ event acceleration (i++) { /** We check for all VOF interfaces for which $\phi$ is allocated. The corresponding volume fraction fields will be stored in *list*. */ scalar * list = NULL; for (scalar f in interfaces) if (f.phi.i) { list = list_add (list, f); /** To avoid undeterminations due to round-off errors, we remove values of the volume fraction larger than one or smaller than zero. */ foreach() f[] = clamp (f[], 0., 1.); boundary ({f}); } /** On trees we need to make sure that the volume fraction gradient is computed exactly like the pressure gradient. This is necessary to ensure well-balancing of the pressure gradient and interfacial force term. To do so, we apply the same prolongation to the volume fraction field as applied to the pressure field. */ #if TREE for (scalar f in list) f.prolongation = p.prolongation; boundary (list); #endif /** Finally, for each interface for which $\phi$ is allocated, we compute the interfacial force acceleration $$\phi\mathbf{n}\delta_s/\rho \approx \alpha\phi\nabla f$$ */ face vector ia = a; foreach_face() for (scalar f in list) if (f[] != f[-1]) { /** We need to compute the potential *phif* on the face, using its values at the center of the cell. If both potentials are defined, we take the average, otherwise we take a single value. If all fails we set the potential to zero: this should happen only because of very pathological cases e.g. weird boundary conditions for the volume fraction. */ scalar phi = f.phi; double phif = (phi[] < nodata && phi[-1] < nodata) ? (phi[] + phi[-1])/2. : phi[] < nodata ? phi[] : phi[-1] < nodata ? phi[-1] : 0.; ia.x[] += alpha.x[]/fm.x[]*phif*(f[] - f[-1])/Delta; } /** On trees, we need to restore the prolongation values for the volume fraction field. */ #if TREE for (scalar f in list) f.prolongation = fraction_refine; boundary (list); #endif /** Finally we free the potential fields and the list of volume fractions. */ for (scalar f in list) { scalar phi = f.phi; delete ({phi}); f.phi.i = 0; } free (list); } /** ## References See Section 3, pages 8-9 of: ~~~bib @Article{popinet2018, author = {S. Popinet}, title = {Numerical models of surface tension}, journal = {Annual Review of Fluid Mechanics}, pages = {1--28}, volume = {50}, year = {2018}, doi = {10.1146/annurev-fluid-122316-045034}, url = {https://hal.archives-ouvertes.fr/hal-01528255} } ~~~ */