Method

Using the second-order accurate finite-difference scheme or the first-order upwind scheme is unfavourable due to its unstable and too diffusive nature, respectively. To integrate the equation stably while preserving interfacial structures which are inherently sharp, we need a procedure called surface reconstruction, which is accomplished by means of the THINC/QQ scheme (Xie and Xiao, J. Comput. Phys. (349), 2017) in cylindrical coordinates.

Surface reconstruction

We consider \(\hat{H}\): the diffused representation of the indicator function \(H\) for each cell. For each cell, we consider a piecewise planar function

\[P \equiv \vec{n} \cdot \vec{x} + d.\]

Note that surface reconstructions are performed on the computational coordinate system rather than the original cylindrical coordinate system. Although this does not give planar representations in the physical coordinate systems, they give similar results to each other since both coordinates are orthogonal, and this greatly simplifies the whole implementation.

To start we focus on finding \(\vec{n}\): surface normal. By definition the direction of \(\vec{n}\) is given by the gradient of the indicator function:

\[\pder{H}{x_i} = \egr \frac{1}{\hvr} \pder{H}{\gvr} + \egt \frac{1}{\hvt} \pder{H}{\gvt} + \egz \frac{1}{\hvz} \pder{H}{\gvz}.\]

Numerically this is approximated as the gradient of \(\phi\) and, following Youngs’ approach, this is first computed at each cell vertex:

\[\vat{\dder{\phi}{x_i}}{\cpidx{i},\cpidx{j},\cpidx{k}} = \egr \frac{1}{\hvr} \dif{\phi}{\gvr} + \egt \frac{1}{\hvt} \dif{\phi}{\gvt} + \egz \frac{1}{\hvz} \dif{\phi}{\gvz},\]

which is implemented as follows:

src/interface/curvature_tensor.c
53#if NDIMS == 2
54    const double dvofdx = 1. / HXXF(i  ) * (
55        - VOF(i-1, j-1) + VOF(i  , j-1)
56        - VOF(i-1, j  ) + VOF(i  , j  )
57    );
58    const double dvofdy = 1. / HYXF(i  ) * (
59        - VOF(i-1, j-1) - VOF(i  , j-1)
60        + VOF(i-1, j  ) + VOF(i  , j  )
61    );
62    const double norm = sqrt(
63        + dvofdx * dvofdx
64        + dvofdy * dvofdy
65    );
66#else
67    const double dvofdx = 1. / HXXF(i  ) * (
68        - VOF(i-1, j-1, k-1) + VOF(i  , j-1, k-1)
69        - VOF(i-1, j  , k-1) + VOF(i  , j  , k-1)
70        - VOF(i-1, j-1, k  ) + VOF(i  , j-1, k  )
71        - VOF(i-1, j  , k  ) + VOF(i  , j  , k  )
72    );
73    const double dvofdy = 1. / HYXF(i  ) * (
74        - VOF(i-1, j-1, k-1) - VOF(i  , j-1, k-1)
75        + VOF(i-1, j  , k-1) + VOF(i  , j  , k-1)
76        - VOF(i-1, j-1, k  ) - VOF(i  , j-1, k  )
77        + VOF(i-1, j  , k  ) + VOF(i  , j  , k  )
78    );
79    const double dvofdz = 1. / hz * (
80        - VOF(i-1, j-1, k-1) - VOF(i  , j-1, k-1)
81        - VOF(i-1, j  , k-1) - VOF(i  , j  , k-1)
82        + VOF(i-1, j-1, k  ) + VOF(i  , j-1, k  )
83        + VOF(i-1, j  , k  ) + VOF(i  , j  , k  )
84    );

The normal vector \(\vec{n}\) is obtained by normalising this, i.e.,

\[\vat{\vec{n}}{\cpidx{i},\cpidx{j},\cpidx{k}} = \egr n_{\gvr} + \egt n_{\gvt} + \egz n_{\gvz},\]

which is implemented as follows:

src/interface/curvature_tensor.c
 92    const double norminv = 1. / fmax(norm, DBL_EPSILON);
 93#if NDIMS == 2
 94    DVOF(i, j)[0] = dvofdx * norminv;
 95    DVOF(i, j)[1] = dvofdy * norminv;
 96#else
 97    DVOF(i, j, k)[0] = dvofdx * norminv;
 98    DVOF(i, j, k)[1] = dvofdy * norminv;
 99    DVOF(i, j, k)[2] = dvofdz * norminv;
100#endif

In addition to the normal vector at cell vertices which are adopted to compute local curvature, we need to compute surface normal at cell centers to find \(P\), which are obtained by averaging the normal vector of the surrounding cell vertices, which are implemented as follows:

src/interface/curvature_tensor.c
138#if NDIMS == 2
139    const double lvof = VOF(i, j);
140#else
141    const double lvof = VOF(i, j, k);
142#endif
143    // for (almost) single-phase region,
144    //   surface reconstruction is not needed
145    if(lvof < vofmin || 1. - vofmin < lvof){
146      continue;
147    }
148#if NDIMS == 2
149    double nx = (
150        + DVOF(i  , j  )[0] + DVOF(i+1, j  )[0]
151        + DVOF(i  , j+1)[0] + DVOF(i+1, j+1)[0]
152    );
153    double ny = (
154        + DVOF(i  , j  )[1] + DVOF(i+1, j  )[1]
155        + DVOF(i  , j+1)[1] + DVOF(i+1, j+1)[1]
156    );
157    const double norm = sqrt(
158        + nx * nx
159        + ny * ny
160    );
161#else
162    double nx = (
163        + DVOF(i  , j  , k  )[0] + DVOF(i+1, j  , k  )[0]
164        + DVOF(i  , j+1, k  )[0] + DVOF(i+1, j+1, k  )[0]
165        + DVOF(i  , j  , k+1)[0] + DVOF(i+1, j  , k+1)[0]
166        + DVOF(i  , j+1, k+1)[0] + DVOF(i+1, j+1, k+1)[0]
167    );
168    double ny = (
169        + DVOF(i  , j  , k  )[1] + DVOF(i+1, j  , k  )[1]
170        + DVOF(i  , j+1, k  )[1] + DVOF(i+1, j+1, k  )[1]
171        + DVOF(i  , j  , k+1)[1] + DVOF(i+1, j  , k+1)[1]
172        + DVOF(i  , j+1, k+1)[1] + DVOF(i+1, j+1, k+1)[1]
173    );
174    double nz = (
175        + DVOF(i  , j  , k  )[2] + DVOF(i+1, j  , k  )[2]
176        + DVOF(i  , j+1, k  )[2] + DVOF(i+1, j+1, k  )[2]
177        + DVOF(i  , j  , k+1)[2] + DVOF(i+1, j  , k+1)[2]
178        + DVOF(i  , j+1, k+1)[2] + DVOF(i+1, j+1, k+1)[2]
179    );
180    const double norm = sqrt(
181        + nx * nx
182        + ny * ny
183        + nz * nz
184    );
185#endif
186    const double norminv = 1. / fmax(norm, DBL_EPSILON);
187    nx *= norminv;
188    ny *= norminv;
189#if NDIMS == 3
190    nz *= norminv;
191#endif

The position of the surface is found by approximating the definition of the volume-of-fluid \(\phi\):

\[\int_{V^\xi} \hat{H} \left( \xi^r, \xi^\theta, \xi^z \right) dV^\xi = \vat{J}{\ccidx{i},\ccidx{j},\ccidx{k}} \vat{\phi}{\ccidx{i},\ccidx{j},\ccidx{k}}\]

as

\[\sum_k \sum_j \sum_i w_i w_j w_k \hat{H} \left( \xi_i^r, \xi_j^\theta, \xi_k^z \right) = \vat{J}{\ccidx{i},\ccidx{j},\ccidx{k}} \vat{\phi}{\ccidx{i},\ccidx{j},\ccidx{k}}.\]

Utilising the so-called mid-point rule (first-order Gauss-Legendre quadrature) for the given computational coordinate system yields

\[\hat{H} \left( 0, 0, 0 \right) = \frac{ 1 }{ 1 + \exp \left( - 2 \beta d \right) } = \vat{\phi}{\ccidx{i},\ccidx{j},\ccidx{k}},\]

or

\[d = - \frac{1}{2 \beta} \log \left( \frac{1}{\vat{\phi}{\ccidx{i},\ccidx{j},\ccidx{k}}} - 1 \right),\]

which is implemented here:

src/interface/curvature_tensor.c
193const double d = - 0.5 / vofbeta * log(1. / lvof - 1.);

Advection

Through the surface reconstruction, we obtain the piecewise linear function \(P\) and its diffused representation \(\hat{H}\) for each cell and are ready to evaluate the fluxes integrated on cell faces:

\[\int_{\partial V^\xi} u_j \hat{H} dS_{\xi^j} = u_j \int_{\partial V^\xi} \hat{H} dS_{\xi^j}\]

to update \(\phi\).

Note that the velocities on the corresponding cell surfaces are taken out of the integrals as they are assumed to be constant on each cell face in the framework of finite-difference methods. To evaluate the surface integrals, we adopt the two-point Gaussian quadratures. Note that, depending on the sign of the local velocity, \(\psi\) computed in the upwind cell should be used due to the stability requirement.

Finally, by normalising the integrated values by the corresponding cell surface areas, we obtain the fluxes as follows.

Radial fluxes:

\[\vat{\mathcal{F}_\vr}{ i \pm \frac{1}{2}, j, k } \equiv \vat{\ur}{i \pm \frac{1}{2}, j, k} \sum_{n = 0}^1 \sum_{m = 0}^1 w_{\gvt_m} w_{\gvz_n} \hat{H} \left( \gvr, \gvt_m, \gvz_n \right)\]
src/interface/update/fluxx.c
 91  BEGIN
 92#if NDIMS == 2
 93    const double lvel = UX(i, j);
 94    double * flux = &FLUXX(i, j);
 95#else
 96    const double lvel = UX(i, j, k);
 97    double * flux = &FLUXX(i, j, k);
 98#endif
 99    const int    ii = lvel < 0. ?     i : i - 1;
100    const double  x = lvel < 0. ? - 0.5 : + 0.5;
101#if NDIMS == 2
102    const double lvof = VOF(ii, j);
103#else
104    const double lvof = VOF(ii, j, k);
105#endif
106    if (lvof < vofmin || 1. - vofmin < lvof) {
107      *flux = lvof;
108    } else {
109      *flux = 0.;
110#if NDIMS == 2
111      for(int jj = 0; jj < NGAUSS; jj++){
112        const double w = gauss_ws[jj];
113        const double y = gauss_ps[jj];
114        *flux += w * indicator(&NORMAL(ii, j), &(const vector_t){x, y});
115      }
116#else
117      for(int kk = 0; kk < NGAUSS; kk++){
118        for(int jj = 0; jj < NGAUSS; jj++){
119          const double w = gauss_ws[jj] * gauss_ws[kk];
120          const double y = gauss_ps[jj];
121          const double z = gauss_ps[kk];
122          *flux += w * indicator(&NORMAL(ii, j, k), &(const vector_t){x, y, z});
123        }
124      }
125#endif
126    }
127    *flux *= lvel;
128  END

Azimuthal fluxes:

\[\vat{\mathcal{F}_\vt}{ i, j \pm \frac{1}{2}, k } \equiv \vat{\ut}{i, j \pm \frac{1}{2}, k} \sum_{n = 0}^1 \sum_{l = 0}^1 w_{\gvr_l} w_{\gvz_n} \hat{H} \left( \gvr_l, \gvt, \gvz_n \right).\]
src/interface/update/fluxy.c
 91  BEGIN
 92#if NDIMS == 2
 93    const double lvel = UY(i, j);
 94    double * flux = &FLUXY(i, j);
 95#else
 96    const double lvel = UY(i, j, k);
 97    double * flux = &FLUXY(i, j, k);
 98#endif
 99    const int    jj = lvel < 0. ?     j : j - 1;
100    const double  y = lvel < 0. ? - 0.5 : + 0.5;
101#if NDIMS == 2
102    const double lvof = VOF(i, jj);
103#else
104    const double lvof = VOF(i, jj, k);
105#endif
106    if (lvof < vofmin || 1. - vofmin < lvof) {
107      *flux = lvof;
108    } else {
109      *flux = 0.;
110#if NDIMS == 2
111      for(int ii = 0; ii < NGAUSS; ii++){
112        const double w = gauss_ws[ii];
113        const double x = gauss_ps[ii];
114        *flux += w * indicator(&NORMAL(i, jj), &(const vector_t){x, y});
115      }
116#else
117      for(int kk = 0; kk < NGAUSS; kk++){
118        for(int ii = 0; ii < NGAUSS; ii++){
119          const double w = gauss_ws[ii] * gauss_ws[kk];
120          const double x = gauss_ps[ii];
121          const double z = gauss_ps[kk];
122          *flux += w * indicator(&NORMAL(i, jj, k), &(const vector_t){x, y, z});
123        }
124      }
125#endif
126    }
127    *flux *= lvel;
128  END

Axial fluxes:

\[\vat{\mathcal{F}_\vz}{ i, j, k \pm \frac{1}{2} } \equiv \vat{\uz}{i, j, k \pm \frac{1}{2}} \sum_{m = 0}^1 \sum_{l = 0}^1 w_{\gvr_l} w_{\gvt_m} \hat{H} \left( \gvr_l, \gvt_m, \gvz \right).\]
src/interface/update/fluxz.c
67BEGIN
68  const double lvel = UZ(i, j, k);
69  double * flux = &FLUXZ(i, j, k);
70  const int    kk = lvel < 0. ?     k : k - 1;
71  const double  z = lvel < 0. ? - 0.5 : + 0.5;
72  const double lvof = VOF(i, j, kk);
73  if (lvof < vofmin || 1. - vofmin < lvof) {
74    *flux = lvof;
75  } else {
76    *flux = 0.;
77    for(int jj = 0; jj < NGAUSS; jj++){
78      for(int ii = 0; ii < NGAUSS; ii++){
79        const double w = gauss_ws[ii] * gauss_ws[jj];
80        const double x = gauss_ps[ii];
81        const double y = gauss_ps[jj];
82        *flux += w * indicator(&NORMAL(i, j, kk), &(const vector_t){x, y, z});
83      }
84    }
85  }
86  *flux *= lvel;
87END

These values are used to update \(\phi\) defined at each cell center \(\left( \ccidx{i},\ccidx{j},\ccidx{k} \right)\) following

\[\pder{\phi}{t} = - \frac{1}{J} \dif{}{\gvr} \left( \jhvr \mathcal{F}_\vr \right) - \frac{1}{J} \dif{}{\gvt} \left( \jhvt \mathcal{F}_\vt \right) - \frac{1}{J} \dif{}{\gvz} \left( \jhvz \mathcal{F}_\vz \right),\]

which is implemented as follows:

src/interface/update/main.c
 93for(int j = 1; j <= jsize; j++){
 94  for(int i = 1; i <= isize; i++){
 95    const double hx_xm = HXXF(i  );
 96    const double hx_xp = HXXF(i+1);
 97    const double hy    = HYXC(i  );
 98    const double jd_xm = JDXF(i  );
 99    const double jd_x0 = JDXC(i  );
100    const double jd_xp = JDXF(i+1);
101    const double flux_xm = FLUXX(i  , j  );
102    const double flux_xp = FLUXX(i+1, j  );
103    const double flux_ym = FLUXY(i  , j  );
104    const double flux_yp = FLUXY(i  , j+1);
105    SRC(i, j) = 1. / jd_x0 * (
106        + jd_xm / hx_xm * flux_xm
107        - jd_xp / hx_xp * flux_xp
108        + jd_x0 / hy    * flux_ym
109        - jd_x0 / hy    * flux_yp
110    );
111  }
112}
src/interface/update/main.c
115for(int k = 1; k <= ksize; k++){
116  for(int j = 1; j <= jsize; j++){
117    for(int i = 1; i <= isize; i++){
118      const double hx_xm = HXXF(i  );
119      const double hx_xp = HXXF(i+1);
120      const double hy    = HYXC(i  );
121      const double jd_xm = JDXF(i  );
122      const double jd_x0 = JDXC(i  );
123      const double jd_xp = JDXF(i+1);
124      const double flux_xm = FLUXX(i  , j  , k  );
125      const double flux_xp = FLUXX(i+1, j  , k  );
126      const double flux_ym = FLUXY(i  , j  , k  );
127      const double flux_yp = FLUXY(i  , j+1, k  );
128      const double flux_zm = FLUXZ(i  , j  , k  );
129      const double flux_zp = FLUXZ(i  , j  , k+1);
130      SRC(i, j, k) = 1. / jd_x0 * (
131          + jd_xm / hx_xm * flux_xm
132          - jd_xp / hx_xp * flux_xp
133          + jd_x0 / hy    * flux_ym
134          - jd_x0 / hy    * flux_yp
135          + jd_x0 / hz    * flux_zm
136          - jd_x0 / hz    * flux_zp
137      );
138    }
139  }
140}