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
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:
Numerically this is approximated as the gradient of \(\phi\) and, following Youngs’ approach, this is first computed at each cell vertex:
which is implemented as follows:
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.,
which is implemented as follows:
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:
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\):
as
Utilising the so-called mid-point rule (first-order Gauss-Legendre quadrature) for the given computational coordinate system yields
or
which is implemented here:
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:
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:
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:
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:
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
which is implemented as follows:
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}
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}