Chapter 5 Results and Discussion
5.2 The VOF schemes test in an oblique velocity field
Two different hollow shapes are convected in an oblique velocity field and Figure 5.1 presents the initial condition of these two shapes. The following describe the information about these two shapes as:
1. A hollow square (Figure 5.1(a)) aligned with the coordinate axes of an outer side length 0.8 and inner side length 0.4.
2. A hollow circle (Figure 5.1(b)) with an outer diameter of length 0.4 and inner diameter of 0.8. Because the volume fraction of each cell can not be obtained easily, there is a special
treatment for the initial setting. In the transition region, the physical surface intersects the cell faces in two points. An approximated surface can be obtained by connecting these two points with a straight line. This line and cell faces build an approximated area (see the hitched area in Figure 5.2). Thus, the portion of an approximated area of a consider cell is the volume fraction value.
The computational domain is a 4× square, subdivided into 4 100×100 grids. All shapes are initially at (0.8, 0.8) with their exact positions centered at (2.8, 1.8) after 1.0 seconds.
Three different time steps, Δt=0.00333, 0.00667 and 0.01, yield Courant number of 0.25, 0.5, 0.75 which will be denoted by low, medium and high, over the uniform grid systems.
Contours are display for
0 . 05 ≤ α ≤ 0 . 95
with 10 levels.For the purpose of comparison the solution errors, the solution error defined as [36]:
∑ the initial condition.
(a)
(b)
Figure 5.1 The flow field domain and initial condition of (a) hollow square, (b) hollow circle
Figure 5.2 Method for catching the value of volume fraction.
The comparison of different schemes
There are many linear and non-linear schemes are considered for solving indicator equation. Linear schemes are upwind difference scheme, central difference scheme, downwind difference scheme and cubic upwind scheme (CUS in Table 3.1). These linear schemes are revised and satisfy the convective boundedness criteria (CBC). In other words, the normalized volume fraction value at a consider face α~ tracks the upwind difference line f in NVD while the donor cell value α~D less than zero or greater than unity. Some famous non-linear schemes are chosen in the test and these schemes (see Figure 5.3 and 5.4) are Bounded Downwind scheme (BD), SMART, MUSCL, OSHER, Van Leer, CHARM and SUPERBEE. The NVD equation and flux limiter function of these non-linear schemes are detailed in Table 3.2.
Among these schemes, upwind difference scheme is a diffusive scheme in each case.
Downwind difference scheme produces highly numerical dispersion and compress shapes, especially at high Courant number (CN=0.75). In hollow circle, downwind difference scheme transforms the original shape into a polygon. Bounded downwind (BD) scheme and SUPERBEE are more accurate but compress shapes in hollow circle case (see Figure 5.4 (c)).
These two schemes are taken as compressive schemes in composite scheme tests. Cubic upwind difference scheme, MUSCL, STOIC and Van Leer are treated as high resolution scheme in the blending strategy when the flow direction perpendicular to interface.
Scheme Upwind Central Difference Downwind Bounded Downwind
Error 1.25E+00 7.59E-01 6.49E-01 3.52E-01
Min: 0.000 Max: 0.472
Min:-0.040 Max:0.989
Min:-0.464 Max:1.467
Min:-0.013 Max:0.949
Scheme Cubic Upwind SMART MUSCL STOIC
Error 5.97E-01 6.50E-01 5.93E-01 6.40E-01
Min:-0.030 Max:1.002
Min:-0.021 Max:0.933
Min:-0.005 Max:0.844
Min:-0.022 Max:0.961
Scheme OSHER Van Leer CHARM SUPERBEE
Error 7.37E-01 6.87E-01 7.47E-01 4.47E-01
Min: 0.000 Max:0.751
Min: 0.000 Max: 0.783
Min: -0.003 Max: 0.798
Min:-0.009 Max0.915:
Figure 5.3 (a) Comparison of different schemes with grids 100*100 and Courant number=0.75 (Δt=0.01)
Scheme Upwind Central Difference Downwind Bounded Downwind
Error 1.25E+00 7.23E-01 4.39E-01 3.51E-01
Min:0.000 Max:0.471
Min:-0.028 Max:0.999
Min:-0.444 Max:1.344
Min:-0.007 Max:0.926
Scheme Cubic Upwind SMART MUSCL STOIC
Error 5.45E-01 6.11E-01 5.81E-01 5.95E-01
Min:-0.017 Max:0.980
Min:-0.002 Max:0.873
Min:0.000 Max:0.842
Min:-0.004 Max:0.899
Scheme OSHER Van Leer CHARM SUPERBEE
Error 7.55E-01 6.76E-01 7.24E-01 4.40E-01
Min:0.000 Max:0.751
Min:0.000 Max:0.774
Min:0.000 Max:0.787
Min:0.000 Max:0.901
Figure 5.3(b) Comparison of different schemes with grids 100*100 and Courant number=0.50 (Δt=0.00667)
Scheme Upwind Central Difference Downwind Bounded Downwind
Error 1.25E+00 6.90E-01 2.03E-01 3.40E-01
Min:0.000 Max:0.471
Min:-0.010 Max:0.954
Min:-0.206 Max:1.190
Min:-0.002 Max:0.943
Scheme Cubic Upwind SMART MUSCL STOIC
Error 5.15E-01 5.72E-01 5.80E-01 5.46E-01
Min:-0.013 Max:0.952
Min:0.000 Max:0.874
Min:0.000 Max:0.851
Min:0.000 Max:0.891
Scheme OSHER Van Leer CHARM SUPERBEE
Error 7.69E-01 6.74E-01 7.04E-01 4.27E-01
Min:0.000 Max:0.761
Min:0.000 Max:0.782
Min:0.000 Max:0.783
Min:0.000 Max:0.914
Figure 5.3(c) Comparison of different schemes with grids 100*100 and Courant number=0.25 (Δt=0.00333)
Scheme Upwind Central Difference Downwind Bounded Downwind
Error 1.23E+00 7.01E-01 7.11E-01 2.25E-01
Min:0.000 Max:0.452
Min:-0.060 Max:0.944
Min:-0.537 Max:1.454
Min:-0.062 Max:1.039
Scheme Cubic Upwind SMART MUSCL STOIC
Error 5.39E-01 5.85E-01 5.18E-01 5.70E-01
Min:-0.045 Max:0.983
Min:-0.026 Max:0.912
Min:-0.006 Max:0.852
Min:-0.03 Max:0.952
Scheme OSHER Van Leer CHARM SUPERBEE
Error 6.57E-01 6.23E-01 6.85E-01 3.17E-01
Min:0.000 Max:0.728
Min:0.000 Max:0.780
Min-0.007 Max:0.796
Min:-0.024 Max:0.990
Figure 5.4(a) Comparison of different schemes with grids 100*100 and Courant number=0.75 (Δt=0.01)
Scheme Upwind Central Difference Downwind Bounded Downwind
Error 1.23E+00 6.46E-01 5.45E-01 2.47E-01
Min:0.000 Max:0.452
Min:-0.034 Max:0.936
Min:-0.380 Max:1.371
Min:-0.037 Max:1.020
Scheme Cubic Upwind SMART MUSCL STOIC
Error 4.67E-01 5.31E-01 4.99E-01 5.08E-01
Min:-0.022 Max:0.969
Min:-0.005 Max:0.873
Min:-0.002 Max:0.823
Min:-0.010 Max:0.905
Scheme OSHER Van Leer CHARM SUPERBEE
Error 6.81E-01 6.07E-01 6.59E-01 3.22E-01
Min:0.000 Max:0.730
Min:0.000 Max:0.761
Min:0.000 Max:0.783
Min:-0.014 Max:0.948
Figure 5.4(b) Comparison of different schemes with grids 100*100 and Courant number=0.50 (Δt=0.00667)
Scheme Upwind Central Difference Downwind Bounded Downwind
Error 1.23E+00 6.03E-01 3.29E-01 2.36E-01
Min:0.000 Max:0.452
Min:-0.012 Max:0.909
Min:-0.211 Max:1.202
Min:-0.016 Max:0.997
Scheme Cubic Upwind SMART MUSCL STOIC
Error 4.22E-01 4.85E-01 4.95E-01 4.50E-01
Min:-0.014 Max:0.943
Min:0.000 Max:0.869
Min:0.000 Max:0.832
Min:0.000 Max:0.913
Scheme OSHER Van Leer CHARM SUPERBEE
Error 6.98E-01 6.03E-01 6.34E-01 3.12E-01
Min:0.000 Max:0.733
Min:0.000 Max:0.760
Min:0.000 Max:0.753
Min:-0.006 Max:0.942
Figure 5.4(c) Comparison of different schemes with grids 100*100 and Courant number=0.25 (Δt=0.00333)
The comparison of composite schemes
As mentioned in the section 3.6, the blending strategy and composite schemes are introduced. When the normal vector of interface is parallel to the flow direction, compressive scheme is used. On the other hand, in the case of the normal vector of interface is perpendicular to the flow direction, it is proper to adopt high resolution scheme. There are three well-known composite schemes, HRIC scheme of Muzaferija [28], CICSAM scheme of Ubbink [30] and STACS of Darwish [29]. The following Figure 5.5 and 5.7 show the hollow square and hollow circle using HRIC, CICSAM and STACS schemes.
HRIC scheme is base on a blending of the Bounded Downwind Scheme and Upwind Differencing Scheme. Because bounded downwind is compressive and introduces negative numerical diffusion, the original shapes hollow circle and hollow square are compressed.
CICSAM scheme is based on the blending the HYPER-C scheme and the ULTIMATE-QUICKEST (UQ) scheme. Similarly, Upwind Different Scheme plays an important role and numerical diffusion is produced at high Courant number. STACS scheme has strong oscillation which caused by Downwind Difference Scheme. These well-known schemes all well performed in the case of hollow square especially at lower Courant number.
Figure 5.6 and Figure 5.8 illustrate the numerical of composite schemes. The former takes as compressive scheme and the later treats high resolution scheme. In the above mentioned, it is clear that Bounded Downwind Scheme and SUPERBEE scheme is more compressive. The combinations of SUPERBEE or Bounded Downwind scheme, which are taken as compressive schemes, with other high resolution schemes (STOIC, MUSCL, Van Leer and cubic upwind difference scheme (CUS)) are tested. Because Bounded Downwind Scheme introduces too much numerical dispersion and compress shapes strongly, SUPERBEE scheme is the better choice. The SUPERBEE/Cubic Upwind scheme introduces lower error but violets the boundedness criteria and cause divergence in the computation process. In the hollow circle test, the SUPEREE/Van Leer and SUPERBEE/STOIC scheme introduce numerical diffusion in the flow direction. By the comparison of numerical results, the blending of SUPERBEE and MUSCL is a proper choice.
The results are performed well at lower Courant number but cost two much computational resource. On the other hand, Numerical diffusion or dispersion will increase with Courant number. In conclusion, to derive reasonable and acceptable numerical simulations, the Courant number has to be limited under a proper value. In the following, the Courant number in each case are described in detail.
Scheme╲Co. 0.25 0.5 0.75
HRIC
Error:3.72E-01 Min:0.000 Max:0.904
Error:1.92E-01 Min:-0.056 Max:1.051
Error:2.43E-01 Min:-0.018
Max:1.024
CICSAM
Error:1.22E-01 Min:0.000 Max:0.990
Error:3.899E-01 Min:0.000 Max:0.887
Error:1.00E+00 Min:0.000 Max:0.603
STACS
Error:1.14E-01 Min:-0.147 Max:1.092
Error:2.567E-01 Min:-0.245 Max:1.237
Error:4.26E-01 Min:-0.379
Max:1.324
Figure 5.5 Comparison of HRIC, CICSAM and STACS with grids 100*100.
Scheme╲Co. 0.25 0.5 0.75
Bounded Downwind
╱ Cubic Upwind
Error:2.92E-01
Scheme╲Co. 0.25 0.5 0.75
SUPERBEE
╱ Cubic Upwind
Error:4.06E-01
Figure 5.6 Comparison of different composite schemes with grids: 100*100
Scheme╲Co. 0.25 0.5 0.75
HRIC
Error:3.25E-01 Min:0.000 Max:0.898
Error:2.39E-01 Min:-0.096 Max:1.021
Error:2.276E-01 Min:-0.080
Max:1.088
CICSAM
Error:1.26E-01 Min:-0.003 Max:1.002
Error:3.96E-01 Min:0.000 Max:0.876
Error:9.99E-01 Min:0.000 Max:0.553
STACS
Error:1.50E-01 Min:-0.152 Max:1.110
Error:2.667E-01 Min:-0.228 Max:1.169
Error:4.13E-01 Min:-0.334
Max:1.259
Figure 5.7 Comparison of HRIC, CICSAM and STACS with grids: 100*100
Scheme╲Co. 0.25 0.5 0.75
Figure 5.8 (continue) Comparison of different composite schemes with grids: 100*100
Scheme╲Co. 0.25 0.5 0.75
Figure 5.8 Comparison of different composite schemes with grids: 100*100.
5.3 Surface tension
Surface tension operates at the interface and produces a force especially at sharp curvature. A numerical test illustrating this phenomenon is to apply it to an initially square bubble [32] and the surface tension force will make it change into circle. The computational domain is 1*1 and on a 30*30 uniform grids with Δx= Δy= 0.0.033cm. The time step size Δt is 0.001. A square bubble with density ρ1=1, viscosity μ1=1 and surface tension σ= 1 is dipped in the outer fluid which has the densityρ2=0.8 and viscosityμ2=0.5.
The time evolution of the changes is shown in Figure 5.9. The relationship between girds units and vector magnitude is 1.5. To see the interface clearly, contour line sets in α equal to 0.5. The largest velocity vectors occur at the corners of the square initially. High surface tension forces produces at these positions which the curvature is maximum. The velocities pull inward on the diagonal direction and pull outward at other parts of surface. To continue, the interface overshoots the stable shape (see Figure 5.9 at t=0.3s) and the velocities pull outward at each original corners. This bubble does a period motion and the amplitude of oscillation will decrease as time goes by. Finally, the bubble will achieve to a stable stage.
The velocity field surrounding the interface has a low amplitude, even in the steady-state.
Lafaurie et al. [37] refer to these velocities as inevitable parasite currents and presented a dimensional analysis of the viscous and surface forces reveal that:
σ
maxμ K V
v
= (5.2)
Where K has no physical meaning, but is useful for predicting the behavior of the parasite currents. Figure 5.10 illustrates different ratio of viscosity and surface tension. The upper left panel shows the final result in the Figure 5.9. The ratioμ1/σ=1 holds in the case and parasite currents are visible around the interface. The next panel uses the ratioμ1/σ=10-1 which the viscosity of fluid 1 is ten times higher surface tension coefficient. The velocity vectors are ten times size of the original vectors. The lower left panel presents the ratioμ1/σ=10. The velocity vectors for this case almost vanish. The last panel shows the results obtained with both the surface tension coefficient and viscosity 10 times than the original values, resulting
in ratio ofμ1/σ=1. The magnitude of the velocity vectors are the same with the first panel.
Hence, these parasite currents increase with decreasing viscosity and increasing surface tension.
t=0.1s t=0.5s
t=1.0s t=3.0s
Figure 5.9 Time evolution of the shape changes of a square subjected to surface tension forces.
μ1=1, σ= 1, ρ1/σ=1 μ1=1, σ= 10,ρ1/σ=10-1
μ1=1, σ=10-1 ,ρ1/σ=10 μ1=10, σ= 10, ρ1/σ=1
Figure 5.10 Parasite currents for different viscosity of fluid 1 and surface tension coefficient.