Temporal discretisation¶
The mass conservation is solved first by means of the THINC method, which is followed by integrating the momentum balance.
To enforce the incompressibility constraint, the updated velocity field is corrected:
\[\newcommand{bef}{n-\frac{1}{2}}
\newcommand{aft}{n+\frac{1}{2}}
\frac{
u_i^{n+1}
-
u_i^*
}{\Delta t}
=
-
\frac{1}{\rho^{n+1}}
\frac{1}{h_{\gcs^i}}
\pder{\psi^{\aft}}{\gcs^i}\]
where \(\psi\) is a scalar potential obtained by solving the variable-coefficient Poisson equation:
\[\frac{1}{J}
\pder{}{\gcs^i}
\left(
\frac{J}{h_{\gcs^i}}
\frac{1}{\rho^{n+1}}
\frac{1}{h_{\gcs^i}}
\pder{\psi^{\aft}}{\gcs^i}
\right)
=
\frac{1}{\Delta t}
\frac{1}{J}
\pder{}{\gcs^i}
\left(
\frac{J}{h_{\gcs^i}}
u_i^*
\right),\]
or discretely:
\[\frac{1}{J}
\dif{}{\gcs^i}
\left(
\frac{J}{h_{\gcs^i}}
\frac{1}{\rho^{n+1}}
\frac{1}{h_{\gcs^i}}
\dif{\psi^{\aft}}{\gcs^i}
\right)
=
\frac{1}{\Delta t}
\frac{1}{J}
\dif{}{\gcs^i}
\left(
\frac{J}{h_{\gcs^i}}
u_i^*
\right).\]
Since this is a variable-coefficient Poisson equation, I adopt the approach proposed by Dodd and Ferrante, J. Comput. Phys. (273), 2014 to utilise the orthogonal decomposition. With this approximation and \(\rho_{ref} \equiv \min \left( 1, \hat{\rho} \right)\), the Poisson equation is modified as
\[\frac{1}{J}
\dif{}{\gcs^i}
\left(
\frac{J}{h_{\gcs^i}}
\frac{1}{h_{\gcs^i}}
\dif{\psi^{\aft}}{\gcs^i}
\right)
=
\frac{1}{J}
\dif{}{\gcs^i}
\left\{
\frac{J}{h_{\gcs^i}}
\left(
\frac{\rho_{ref}}{\rho^{n+1}}
-
1
\right)
\frac{1}{h_{\gcs^i}}
\dif{\psi^{\bef}}{\gcs^i}
\right\}
+
\frac{\rho_{ref}}{\Delta t}
\frac{1}{J}
\dif{}{\gcs^i}
\left(
\frac{J}{h_{\gcs^i}}
u_i^*
\right).\]
src/fluid/compute_potential.c¶
497const double den_xm = + 0.5 * DEN(i-1, j )
498 + 0.5 * DEN(i , j );
499const double den_xp = + 0.5 * DEN(i , j )
500 + 0.5 * DEN(i+1, j );
501const double den_ym = + 0.5 * DEN(i , j-1)
502 + 0.5 * DEN(i , j );
503const double den_yp = + 0.5 * DEN(i , j )
504 + 0.5 * DEN(i , j+1);
505const double dpsi_xm = - PSI(i-1, j )
506 + PSI(i , j );
507const double dpsi_xp = - PSI(i , j )
508 + PSI(i+1, j );
509const double dpsi_ym = - PSI(i , j-1)
510 + PSI(i , j );
511const double dpsi_yp = - PSI(i , j )
512 + PSI(i , j+1);
513const double gp_xm = (refden / den_xm - 1.) * dpsi_xm / hx_xm;
514const double gp_xp = (refden / den_xp - 1.) * dpsi_xp / hx_xp;
515const double gp_ym = (refden / den_ym - 1.) * dpsi_ym / hy;
516const double gp_yp = (refden / den_yp - 1.) * dpsi_yp / hy;
517const double add = 1. / jd_x0 * (
518 - jd_xm / hx_xm * gp_xm + jd_xp / hx_xp * gp_xp
519 - jd_x0 / hy * gp_ym + jd_x0 / hy * gp_yp
520);
src/fluid/compute_potential.c¶
549const double den_xm = + 0.5 * DEN(i-1, j , k )
550 + 0.5 * DEN(i , j , k );
551const double den_xp = + 0.5 * DEN(i , j , k )
552 + 0.5 * DEN(i+1, j , k );
553const double den_ym = + 0.5 * DEN(i , j-1, k )
554 + 0.5 * DEN(i , j , k );
555const double den_yp = + 0.5 * DEN(i , j , k )
556 + 0.5 * DEN(i , j+1, k );
557const double den_zm = + 0.5 * DEN(i , j , k-1)
558 + 0.5 * DEN(i , j , k );
559const double den_zp = + 0.5 * DEN(i , j , k )
560 + 0.5 * DEN(i , j , k+1);
561const double dpsi_xm = - PSI(i-1, j , k )
562 + PSI(i , j , k );
563const double dpsi_xp = - PSI(i , j , k )
564 + PSI(i+1, j , k );
565const double dpsi_ym = - PSI(i , j-1, k )
566 + PSI(i , j , k );
567const double dpsi_yp = - PSI(i , j , k )
568 + PSI(i , j+1, k );
569const double dpsi_zm = - PSI(i , j , k-1)
570 + PSI(i , j , k );
571const double dpsi_zp = - PSI(i , j , k )
572 + PSI(i , j , k+1);
573const double gp_xm = (refden / den_xm - 1.) * dpsi_xm / hx_xm;
574const double gp_xp = (refden / den_xp - 1.) * dpsi_xp / hx_xp;
575const double gp_ym = (refden / den_ym - 1.) * dpsi_ym / hy;
576const double gp_yp = (refden / den_yp - 1.) * dpsi_yp / hy;
577const double gp_zm = (refden / den_zm - 1.) * dpsi_zm / hz;
578const double gp_zp = (refden / den_zp - 1.) * dpsi_zp / hz;
579const double add = 1. / jd_x0 * (
580 - jd_xm / hx_xm * gp_xm + jd_xp / hx_xp * gp_xp
581 - jd_x0 / hy * gp_ym + jd_x0 / hy * gp_yp
582 - jd_x0 / hz * gp_zm + jd_x0 / hz * gp_zp
583);
The velocity corrections are modified as
\[\newcommand{\old}[1]{
+
\left(
\frac{1}{\rho^{n+1}}
-
\frac{1}{\rho_{ref}}
\right)
\frac{1}{h_{\gcs^#1}}
\pder{\psi^{\bef}}{\gcs^#1}
}
\newcommand{\new}[1]{
-
\frac{1}{\rho_{ref}}
\frac{1}{h_{\gcs^#1}}
\pder{\psi^{\aft}}{\gcs^#1}
}
\frac{
u_i^{n+1}
-
u_i^*
}{\Delta t}
=
\new{i}
\old{i}.\]
There are two contributions, which are from the new scalar potential and the old scalar potential.
The contributions of the new scalar potential:
\[\new{\vx}\]
src/fluid/correct_velocity/ux.c¶
36for(int j = 1; j <= jsize; j++){
37 for(int i = 2; i <= isize; i++){
38 const double hx = HXXF(i );
39 const double psi_xm = PSI(i-1, j );
40 const double psi_xp = PSI(i , j );
41 UX(i, j) -= dt_new / refden / hx * (
42 - psi_xm
43 + psi_xp
44 );
45 }
46}
src/fluid/correct_velocity/ux.c¶
49for(int k = 1; k <= ksize; k++){
50 for(int j = 1; j <= jsize; j++){
51 for(int i = 2; i <= isize; i++){
52 const double hx = HXXF(i );
53 const double psi_xm = PSI(i-1, j , k );
54 const double psi_xp = PSI(i , j , k );
55 UX(i, j, k) -= dt_new / refden / hx * (
56 - psi_xm
57 + psi_xp
58 );
59 }
60 }
61}
\[\new{\vy}\]
src/fluid/correct_velocity/uy.c¶
36for(int j = 1; j <= jsize; j++){
37 for(int i = 1; i <= isize; i++){
38 const double psi_ym = PSI(i , j-1);
39 const double psi_yp = PSI(i , j );
40 UY(i, j) -= dt_new / refden / hy * (
41 - psi_ym
42 + psi_yp
43 );
44 }
45}
src/fluid/correct_velocity/uy.c¶
48for(int k = 1; k <= ksize; k++){
49 for(int j = 1; j <= jsize; j++){
50 for(int i = 1; i <= isize; i++){
51 const double psi_ym = PSI(i , j-1, k );
52 const double psi_yp = PSI(i , j , k );
53 UY(i, j, k) -= dt_new / refden / hy * (
54 - psi_ym
55 + psi_yp
56 );
57 }
58 }
59}
\[\new{\vz}\]
src/fluid/correct_velocity/uz.c¶
33for(int k = 1; k <= ksize; k++){
34 for(int j = 1; j <= jsize; j++){
35 for(int i = 1; i <= isize; i++){
36 const double psi_zm = PSI(i , j , k-1);
37 const double psi_zp = PSI(i , j , k );
38 UZ(i, j, k) -= dt_new / refden / hz * (
39 - psi_zm
40 + psi_zp
41 );
42 }
43 }
44}
The contributions of the old scalar potential:
\[\old{\vx}\]
src/fluid/correct_velocity/ux.c¶
70for(int j = 1; j <= jsize; j++){
71 for(int i = 2; i <= isize; i++){
72 const double hx = HXXF(i );
73 const double psi_xm = PSI(i-1, j );
74 const double psi_xp = PSI(i , j );
75 const double den_x0 = + 0.5 * DEN(i-1, j )
76 + 0.5 * DEN(i , j );
77 UX(i, j) += dt_new * coef * (1. / den_x0 - 1. / refden) / hx * (
78 - psi_xm
79 + psi_xp
80 );
81 }
82}
src/fluid/correct_velocity/ux.c¶
85for(int k = 1; k <= ksize; k++){
86 for(int j = 1; j <= jsize; j++){
87 for(int i = 2; i <= isize; i++){
88 const double hx = HXXF(i );
89 const double psi_xm = PSI(i-1, j , k );
90 const double psi_xp = PSI(i , j , k );
91 const double den_x0 = + 0.5 * DEN(i-1, j , k )
92 + 0.5 * DEN(i , j , k );
93 UX(i, j, k) += dt_new * coef * (1. / den_x0 - 1. / refden) / hx * (
94 - psi_xm
95 + psi_xp
96 );
97 }
98 }
99}
\[\old{\vy}\]
src/fluid/correct_velocity/uy.c¶
68for(int j = 1; j <= jsize; j++){
69 for(int i = 1; i <= isize; i++){
70 const double psi_ym = PSI(i , j-1);
71 const double psi_yp = PSI(i , j );
72 const double den_y0 = + 0.5 * DEN(i , j-1)
73 + 0.5 * DEN(i , j );
74 UY(i, j) += dt_new * coef * (1. / den_y0 - 1. / refden) / hy * (
75 - psi_ym
76 + psi_yp
77 );
78 }
79}
src/fluid/correct_velocity/uy.c¶
82for(int k = 1; k <= ksize; k++){
83 for(int j = 1; j <= jsize; j++){
84 for(int i = 1; i <= isize; i++){
85 double psi_ym = PSI(i , j-1, k );
86 double psi_yp = PSI(i , j , k );
87 const double den_y0 = + 0.5 * DEN(i , j-1, k )
88 + 0.5 * DEN(i , j , k );
89 UY(i, j, k) += dt_new * coef * (1. / den_y0 - 1. / refden) / hy * (
90 - psi_ym
91 + psi_yp
92 );
93 }
94 }
95}
\[\old{\vz}\]
src/fluid/correct_velocity/uz.c¶
51for(int k = 1; k <= ksize; k++){
52 for(int j = 1; j <= jsize; j++){
53 for(int i = 1; i <= isize; i++){
54 const double psi_zm = PSI(i , j , k-1);
55 const double psi_zp = PSI(i , j , k );
56 const double den_z0 = + 0.5 * DEN(i , j , k-1)
57 + 0.5 * DEN(i , j , k );
58 UZ(i, j, k) += dt_new * coef * (1. / den_z0 - 1. / refden) / hz * (
59 - psi_zm
60 + psi_zp
61 );
62 }
63 }
64}
Since I treat the diffusive terms fully-explicitly, the scalar potential and the pressure are simply related by
\[p^{n+1}
-
p^n
=
\psi^{\aft}.\]
src/fluid/update_pressure.c¶
20for(int j = 1; j <= jsize; j++){
21 for(int i = 1; i <= isize; i++){
22 P(i, j) += PSI(i, j);
23 }
24}
src/fluid/update_pressure.c¶
27for(int k = 1; k <= ksize; k++){
28 for(int j = 1; j <= jsize; j++){
29 for(int i = 1; i <= isize; i++){
30 P(i, j, k) += PSI(i, j, k);
31 }
32 }
33}
This is followed by updating \(\psi\):
src/fluid/update_pressure.c¶
43double * tmp = psi0->data;
44psi0->data = psi1->data;
45psi1->data = tmp;