P. Wesseling, E. O˜nate and J. P´eriaux (Eds) c
TU Delft, The Netherlands, 2006
MODELING OF THE WAVE BREAKING WITH CICSAM
AND HRIC HIGH-RESOLUTION SCHEMES
Tomasz Wac lawczyk∗, Tadeusz Koronowicz†
∗†The Szewalski Institute of Fluid Flow Machinery
Polish Academy of Sciences, Centre for Mechanics of Liquid, ul. J. Fiszera 14, Gda´nsk 80-952, Poland,
e-mail: tomi@imp.gda.pl web page: http://www.imp.gda.pl
Key words: Free surface flow, Volume of Fluid method, High-Resolution schemes Abstract. The paper concerns modeling of two-phase flow with the volume of fluid method (VOF) and two high-resolution advection schemes based on the normalized variable dia-gram (NVD). Compressive Interface Capturing Scheme for Arbitrary Meshes (CICSAM) and High Resolution Interface Capturing scheme (HRIC). Both considered schemes are used to discretize convective term in the scalar equation for the transport of the volume fraction. High-resolution schemes are employed to minimize influence of the artificial numerical dissipation and to keep the shape of the step interface profile. Original contri-bution of this work is a detailed comparison of the two high-resolution schemes CICSAM and HRIC in the case of the breaking wave phenomenon. It is shown that using relatively simple to apply, high-resolution schemes it is possible to obtain good agreement with an experimental evidence and other authors results. However, some difficulties connected with optimal size of the local Courant number are adressed. Additionally, since HRIC scheme is used in commercial software (e.g. Comet, Fluent) it is important to understand its capabilities and deficiencies when compared to other schemes.
1 INTRODUCTION
Numerical modelling of the multiphase flows is challenging field of the computational fluid dynamics, which experiences rapid growth. Recently, several interesting numer-ical approaches to the problem were introduced and developed, among them: level-set methods1,2, particle finite-element method3,4 or smoothed particle hydrodynamics5,6. Moreover, the well established Volume of Fluid7 method gained a large number of
im-provements to the interface reconstruction techniques8,9and to the modelling of the surface
tension10,11.
position is captured using two high-resolution advection schemes CICSAM12 and HRIC13,
based on the normalised variable diagram NVD14. Unlike geometric interface reconstruc-tion methods, high-resolureconstruc-tion schemes do not introduce geometric representareconstruc-tion of the interface. Here to avoid smearing of the step interface profile and preserve boundedness criterion, properly designed discretization scheme is applied.
To model dynamics of the system of immiscible fluids, the solution of the equation for transport of the volume fraction is coupled with an in-house code solving Navier-Stokes equations. Implicit finite-volume second order accurate code employs SIMPLE procedure15,16 for pressure-velocity coupling.
2 EQUATIONS GOVERNING THE FLOW OF IMMISCIBLE FLUIDS
To model dynamics of the set of immiscible, incompressible and viscous fluids one needs to solve Navier-Stokes equation (1) and continuity equation (2):
∂ρui ∂t + ∂ρujui ∂xj =−∂p ∂xi + ∂ ∂xj " µ ∂ui ∂xj + ∂uj ∂xi !# + ρgi+ σκniδs, (1) ∂ui ∂xi = 0, (2)
where uidenotes i-th velocity component, ρ is density, p is pressure, µ is dynamic viscosity,
gi denotes i-th component of the gravitational acceleration. Additional term that appears
in Eq. (1) is inherited from one fluid formulation and describes influence of the surface tension. When surface tension coefficient σ is constant the resulting force acts in the direction normal to the interface ni. κ = ∇~n denotes mean Gaussian curvature, δs is
Dirac delta function that indicates local character of the surface tension term, since it is non-zero only at the interface.
In this paper so called one-fluid formulation is used i.e. one treats the whole system of immiscible fluids as a single continuum medium. This assumption leads to the sharing of the same velocity in one control volume, between all fluids that are modelled. This last fact is clearly a disadvantage of the one-fluid formulation since one is not able to model directly effects of the slip phenomenon.
To introduce material properties of the fluids that built aforementioned system, one needs to consider additional constitutive relations. Since in the case of the present work only dynamics of the set of two fluids (water and air) is modelled, it is enough to choose one fluid as the background (water) and to consider only its volume fraction φ. Taking this into account, one obtains:
ρ = (1− φ)ρ1+ φρ2, (3)
µ = (1− φ)µ1+ φµ2, (4)
where ρj, µj denotes density and dynamic viscosity of the j-th fluid respectively. Finally,
that substitution of the Eq. (3) to the conservative form of the continuity equation leads to the equation for transport of the volume fraction that writes:
∂φ ∂t + uj
∂φ ∂xj
= 0, (5)
where, as mentioned before, φ denotes volume fraction of the background fluid. In the case of the VOF framework, value of the volume fraction φ in the control volume indicates its presence φ = 1 or absence φ = 0. When 0 < φ < 1 volume fraction distribution carries information about position of the interface.
The crucial issue for the modelling of the multiphase flow is a proper solution of Eq. (5) i.e. the discretization of time derivative and nonlinear convective term. The discretization have to avoid smearing of the step interface profile and assures that the boundedness criterion is satisfied. Here, Eq. (5) is discretized using finite volume method, employing Crank-Nicolson method for integration in time. This discretization procedure gives the implicit, unsplit time discretization scheme12:
VP ∆tφ t+∆t P + 1 2 n X f =1 φt+∆tf Sfuti,f = VP ∆tφ t P − 1 2 n X f =1 φtfSfuti,f, (6)
where VP denotes the volume of the control volume P , n is the number of faces of the
control volume, Sf is the face area and finally ∆t denotes the time step size. Convective
flux (value of the volume fraction at the face of the control volume φf) is obtained from
two high-resolution schemes CICSAM12 and HRIC13 that are compared in this work.
3 HIGH-RESOLUTION SCHEMES
Discretization of the convective term in Eq. (5) is critical for proper capturing of the interface position. One can show17 that commonly used schemes e.g. upwind (UDS), central (CDS) differencing schemes, introduce effects of artificial diffusion or dispersion respectively to their order of accuracy. Also other higher order schemes result in local oscillations18 of the variable, here volume fraction, that is unacceptable when one wants to employ them for the discretization of convective term in Eq. (5). To avoid this artificial effects high-resolution schemes are built that assure both lack of the numerical diffusion and compressive character i.e. sharpening of the step interface profile12,13.
3.1 Normalised Variable Diagram NVD
Normalised Variable Diagram14 (NVD) provides the foundation for both CICSAM12
and HRIC13 schemes. It is based on the convective boundedness criterion12,13,14 (CBC) that states that the variable distribution between the centers of the neighbourhood control volumes e.g. D and A should remain smooth φD ≤ φf ≤ φA, see Fig. 1a. Using this
constraint, and information about value of the variable in the upwind control volume φU,
f A D U V A φU φD φ φf φ ~ D UD 1 0.5 1 DD 0.5 ~ φf a) b)
Figure 1: a) Boundedness criterion, U upwind, D donor, A acceptor cells, b) Normalized Variable Diagram NVD, shaded region indicates where CBC is satisfied, notice that UDS satisfies CBC criterion in the whole range of eφD e φf = φf − φU φA− φU , (7) e φD = φD− φU φA− φU . (8)
Equations (7), (8) allow to rewrite the boundedness criterion basing on normalised vari-ables φeD ≤ φef ≤ 1, what can be depicted on a diagram, cf. Fig. 1b, that shows region
where the boundedness criterion is satisfied by any differential scheme. One needs to notice that when using Eq. (7)-(8) it is possible to find the value of the volume fraction at the control volume face φf:
φf = 1−βef φD +βefφA, (9) e βf = e φf −φeD 1−φeD , (10)
which is used for the calculation of the volume fraction flux in Eq. (6).
3.2 Compressive Interface Capturing Scheme for Arbitrary Meshes CICSAM In the case of the CICSAM scheme, additional assumption about the dependence of the region where the CBC is satisfied on the CFL condition is used12. The local value of
the Courant number, defined at the face of the control volume Cf = ufSf∆t/V together
with CBC gives the following formula φeD ≤ φef ≤ min
1,φeD/Cf
, which can be plotted at NVD, see Fig. 2. One needs to notice, that for explicit schemes, if local value of the Courant number is equal to Cf = 1, only the UDS satisfies the CBC criterion. On the
other hand, the smaller the value of the local Cf the lager the domain where the CBC
Cf φ~D UD 1 0.5 1 DD 0.5 ~ φf
Figure 2: Dependence of the CBC region on local CFL condition
The combining of the donor-acceptor7 scheme (DAS) with NVD formulation
depen-dent on the CFL condition, results in the first component of the CICSAM known as the HYPER-C12 scheme: e φfCBC = e φf =φeD : 0 < φeD,φeD > 1 min 1,φeD Cf : 0 ≤ φeD ≤ 1 . (11)
The HYPER-C scheme Eq. (11) satisfies the CBC criterion and is shown to be compres-sive, which means that it changes any gradient to a step profile due to the downwind differencing scheme employed (DDS).
However, compressive character of the HYPER-C is not always desirable. When inter-face is tangential to the flow direction it is shown that aforementioned scheme tends to artificially deform its shape. For this reason it is found to be necessary to switch between scheme Eq. (11) and other less compressive formulation. In the case of the CICSAM, the ULTIMATE-QUICKEST14 (UQ) scheme is employed, that is bounded version of the
third order accurate QUICK:
e φfU Q = e φD : 0 < φeD,φeD > 1 min 8CfeφD+(1−Cf)(6eφD+3) 8 ,φefCBC : 0 ≤ φeD ≤ 1 . (12)
To switch smoothly between both schemes, linear blending is used, with blending factor 0≤ γf ≤ 1. Value of the γf depends on the angle θf Eq. (14) between unit vector normal
to the interface ~n = ∇φD/|∇φD| and unit vector parallel to the line between centers of
the donor D and acceptor A cells ~d = ~DA/| ~DA|, cf. Fig. 3a. When interface position is normal to the direction of the flow γf = 1 and Eq. (11) is used. In the case of tangential
orientation of the interface γf = 0 and Eq. (12) is employed. e
φf = γfφef
θf = arccos ~ d~n , (14) γf = min ( 1 + cos 2θf 2 , 1 ) . (15)
It should be noticed that above derivation of the high-resolution CICSAM scheme was carried out only in one-dimension. To extend above formulation to multiply dimensions cell Courant number is introduced CD = Pnf =1max(Cf, 0), where n denotes number of
the control volume faces.
3.3 High Resolution Interface Capturing scheme HRIC
To simplify above procedure an get rid of the explicit dependence on the CFL condition the HRIC13 scheme was introduced. As it was mentioned this scheme also relies on the
NVD and normalised variables. Application of the HRIC scheme can be divided in to
d n θf D A φf ~∗∗ φf ~* Cf 0.7 0.3 φ ~ D 1 a) b)
Figure 3: a) Definition of the vector normal to the interface ~n = ∇φ/|∇φ| and vector parallel to the line between centers of D, A control volumes ~d = ~DA/| ~DA|, b) Continuity of the HRIC scheme in time domain
three steps. Firstly, normalised cell face value will be estimated from a scheme that continuously connects upwind and downwind schemes on the NVD diagram, cf. Fig. 1b:
e φf = e φD : φeD < 0,φeD > 1 2φeD : 0 ≤ φeD < 0.5 1 : 0.5 ≤ φeD ≤ 1 . (16)
Secondly, since DDS can cause alignment of the interface with the mesh and its artificial deformation (as in the case of the HYPER-C scheme) one needs other scheme that satisfies the CBC. In the case of the HRIC, as the most straightforward choice the UDS scheme is employed. Again the blending factor γf connected with angle θf is introduced, see Eq.
(18), to switch smoothly between the schemes.
e
γf = q
| cos θ|. (18)
Blending of the UDS and the DDS schemes is dynamic and takes into account local distribution of the volume fraction. In the case when the CFL condition is not satisfied the dynamic character of this scheme can cause stability problems. Therefore,φe∗
f, see Eq.
(17), is corrected with respect to the local Courant number Cf, cf. Eq. (19). The goal of
this correction is to force continuous switching between schemes also in time domain, cf. Fig. 3b. e φ∗∗f = e φ∗f : Cf < 0.3 e φD : Cf > 0.7 e φD+ (φe∗ f −φeD)0.7−C0.7−0.3f : 0.3 ≤ Cf ≤ 0.7 . (19)
When using this scheme in multiple dimensions the local Courant number Cf is again
replaced by its cell definition CD.
One can notice that main difference between the CICSAM and the HRIC are the order of accuracy of the component schemes and differently included CFL condition. Next several tests will be performed to asses how those differences influence results obtained with two presented high-resolution schemes.
4 ADVECTION TEST CASES
It is known that when convection is the only transport phenomenon, distance between two points advected in the uniform flow field remains constant. Therefore, properties of the high-resolution schemes are assessed using known in literature advection test cases (only Eq. (5) is solved, velocity field is assumed to be known). Firstly, advection of the circular spot in uniform, oblique velocity field is investigated, cf. Fig. 4. During this test case, orders of accuracy of the CICSAM and the HRIC are estimated.
r = [0.15,0.15]mI r = [1.55,0.85]m u = [0.1, 0.05]m/s 2m 1m F φ=0 φ=1 outlet inlet
Moreover, since both formulations depend on the value of the cell Courant number its influence is checked.
The second solid body revolution test case, see Fig. 5, is performed to check properties of both schemes concerning conservation of the shape of the interface in more compli-cated rotational velocity field. Results obtained with the HRIC and the CICSAM high-resolution schemes are compared with two versions of constrained interpolation profile (CIP) semi-Lagrangian scheme18.
symmetry plane ω s=0.1 dy=1 R=0.2 dx=1
Figure 5: Revolution of the solid body, advection test case
4.1 Circle advection test
The first, circle advection test case is carried out to check shape preserving properties of both the CICSAM and the HRIC. Since both schemes are convolved with the local value of the cell Courant number, see Eq. (11)-(12), Eq. (19), advection test is carried out using four different cell Courant numbers CD, see Fig. 6. The values of CD are
chosen to check the characteristic points in the range 0 < CD < 1, cf. Fig. 3b. Figure 4,
depicts computational domain and applied boundary conditions: inlet at the south and west boundaries and outlet on the east and north boundaries, initial ~rI = [0.15, 0.15] and
final ~rF = [1.55, 0.85] position of the spot. Let us notice that in the case of the largest
value of the cell Courant number CD = 0.7, the position of the bubble is captured at
~rF = [0.73, 0.43]. Because it was the last position where the circular spot was observed
on the coarsest 64× 32 grid, compare Fig. 6a and Fig. 6d.
On the first grid 64× 32, results obtained with the CICSAM and the HRIC both suffer from an artificial cell alignment, cf. Fig. 6a-d, i.e. initial circular spot shape deform to the square one. However, one can notice that when CD = 0.1, shape reconstruction is
CICSAM scheme dependence on large CD values is much stronger than in the case of the
HRIC. This effect might be explained by implicit dependence of the CICSAM scheme on the CFL condition, while the HRIC includes it only in the one corrector step. In the case of the second 128× 64 grid the superiority of the CICSAM, in circle shape reconstruction over the HRIC is clearly visible, see CD = 0.1, cf. Fig. 6e.
64× 32, CD= 0.1 64× 32, CD= 0.3 64× 32, CD= 0.5 64× 32, CD= 0.7 a) b) c) d) 128× 64, CD= 0.1 128× 64, CD= 0.3 128× 64, CD= 0.5 128× 64, CD= 0.7 e) f) g) h) 256× 128, CD= 0.1 256× 128, CD= 0.3 256× 128, CD= 0.5 256× 128, CD= 0.7 i) j) k) l) 512× 256, CD= 0.1 512× 256, CD= 0.3 512× 256, CD= 0.5 512× 256, CD= 0.7 m) n) o) p)
Figure 6: Comparison of the CICSAM and the HRIC shape preserving properties on four gradually refined grids and using four values of the cell Courant number CD, figures depicts contours of the volume fraction φ = 0.5
The circle shape is much better preserved by the first aforementioned scheme, while the latter still suffers from mesh alignment problems. Again, with increasing CD shape
Solution obtained on the third grid 256× 128 shows that both high-resolution schemes start to fit the exact shape of the circular spot, when CD = 0.1 and CD = 0.3. It is
interesting to notice that for CD = 0.5, locally, mesh alignment phenomenon can be
observed in the case of both formulations, cf. Fig. 6k. For the CD = 0.7 this phenomenon
is not visible at Fig. 6l but one needs to keep in mind that this figure is made half way from the final position of the spot ~rF. The results obtained on the last grid 512× 256,
show that both schemes converge to the exact solution for both CD = 0.1 and CD = 0.3.
Additionally, one can notice that the error in reconstruction of the shape is almost the same, see Fig. 7.
CD = 0.1 a) CD = 0.3 b) CD = 0.5 c) CD = 0.7 d)
Figure 7: Order of accuracy of the CICSAM and the HRIC high-resolution schemes calculated from results on four gradually refined grids a) CD= 0.1, b) CD= 0.3, c) CD= 0.5, d) CD= 0.7
As it was mentioned, since exact solution of this problem is known, results of the test case where used to estimate order of accuracy of both high-resolution schemes. The error of the numerical solution was calculated from the following average error formula15:
where φexact
i , φi are exact and numerical solutions, respectively, N is a number of volume
fraction contour samples. From Fig. 7a-c it can be seen that both the CICSAM and the HRIC are first order accurate12. In the case of the CD = 0.7 one can notice that order of
accuracy of the method (approximately, since this error estimation is made for different final position of the spot) is smaller then one. Moreover, in the case of the largest CD,
average errors of the CICSAM scheme are larger then the average errors of the HRIC, cf. Fig. 7d, on each considered grid. When the mesh is gradually refined, the difference in error of the solution between both high-resolution schemes becomes smaller independently on the value of the Courant number CD. Finally for the finest grid the average error of
the solution is almost the same for all values of CD.
4.2 Revolution of the solid body
The second advection test case was performed to compare the shape preserving prop-erties of the CICSAM and the HRIC high-resolution schemes with the CIP18
semi-Lagrangian scheme in rotational velocity field.
a) b) c)
d) e)
Figure 8: Volume fraction distribution after one revolution of the solid body, CD ≈ 0.15 at the circle contour: a) initial condition, b) CICSAM, c) HRIC, d) rational CIP18, e) tangent-transformed CIP18
In this test case, circular solid body with a slot is revolved with constant angular velocity ω that is chosen to satisfy CD ≈ 0.15 at the circle edges. Final solution was
initial shape of the solid body and solutions obtained with aforementioned high-resolution schemes. Comparing results obtained with the CICSAM and the HRIC, one can notice, that the first one preserves the shape of the interface better then the latter, compare Fig. 8b and Fig. 8c. The interface calculated with the HRIC is more smeared, non-zero volume fraction values are visible in the slot after one revolution. Conclusion from this test case is that the CICSAM does better work, what is in agreement with the results obtained in the first circular spot advection test for small CD values.
Finally, the CICSAM and the HRIC results can be compared with results obtained using two CIP18 formulations. It is noticeable, that the CICSAM shape preserving prop-erties lie between capabilities of the rational CIP18 and the tangent-transformed CIP18
formulations, however, closer to the last mentioned scheme. Unlike the CICSAM results, the HRIC schemes shape preserving properties are closer to rational CIP formulation, however one needs to notice that interface deformation is smaller in the case of the HRIC scheme.
Since the constrained interpolation profile CIP18 formulation is a semi-Lagrangian one i.e. additionally three equations are solved to trace advection of each component of the gradient of the volume fraction, with step interface profile constrained by higher-order polynomial function. In opinion of the authors comparison of the CICSAM with the last very good result obtained with tangent-transformed CIP18 should be found satisfactory.
5 BREAKING OF THE WATER WAVE
Next simulation is performed to check how both considered schemes can manage with the real physical example of the multiphase flow. In this case Eq. (5) is coupled with Eqs.(1)-(4), thus the whole set of equations that describe flow of the two fluids is solved. Here, the water column with height H = 600mm and width 2H, see Fig. 9, collapses in the gravitational field. Results of the simulations are compared with an experimental data6 consisting of height of the water surface recorded at two positions H1, H2 and the pressure history that was measured at point P1, cf. Fig. 9. Moreover, flow field parameters obtained with the CICSAM and the HRIC schemes are compared with the SPH6 (Smoothed Particle Hydrodynamics) method and with results of the LS19 (Level-Set) method.
Initial volume fraction distribution is shown in Fig. 9, initial pressure field is set accord-ing to the hydrostatic pressure connected with density distribution, initial velocity field is set to zero. Material properties of fluids are for water: ρ1 = 998 kg/m3, µ1 = 0.99× 10−5
kg/(ms) and for air: ρ2 = 1.205 kg/m3, µ2 = 1.81× 10−3 kg/(ms). To trace behaviour
of the water-air interface a non-dimensional time T , based on initial water column height H is introduced:
T = tqg/H. (21)
slip wall no g=9.8 m/s2 H H2 H1 5.366H 3H 3.713H 4.541H φ=0 0.26H P1 2H φ=1
Figure 9: Computational domain for the breaking wave test case, H denotes initial height of the water column, H1, H2 positions of the water height probes, P1 position of the pressure probe
value of the cell Courant number is CD = 0.2.
5.1 Comparison with water height history from H1, H2 probes
Initially water column collapses in the gravitational field colliding with the opposite wall (see Fig. 12, time T1, T2). After collision, it climbs up across the west wall of the
tank (see Fig. 12, time T3, T4), afterwards, when its whole kinetic energy changes into
potential energy it starts to move downward. During descending motion the water wave is formed that finally overturns and breaks, (cf. Fig. 12, time T5).
Results obtained during the simulation compare well with experimental evidence and results of the SPH method6, when T < 6.3 in the case of height probe H1 and for times T < 5.7 for probe H2, cf. Fig. 10. One can notice that water height predicted with the CICSAM and the HRIC schemes fit experiment better in the case of H1 probe when water initially changes its height, see Fig. 10, (H1 probe about T = 1.5). The differences between experiment and numerical results obtained with high-resolution schemes is noticeable, but character of the flow is preserved. The results start to differ more after the water wave overturns and breaks. The sudden peak in pressure history might be observed, see Fig. 11 about T = 6 when broken wave hits the surface of the water. The same phenomenon is recorded at Fig. 10 in rapid increase of water level first noticed at H2 probe and then observed at H1 probe. The data recorded by the probe H1 indicate that both the CICSAM and the HRIC schemes predict correctly maximal height of the water wave, slightly better then the SPH method6. However, change of the water level when T > 6.7
is taken into account and compressibility effects are partially included. Despite the fact that in current implementation of the VOF method both fluids are incompressible, results obtained are closer to experimental evidence then SPH results6 where only dynamics of the water is taken in to account, compare Fig. 10. It proves that it is necessary include both the motion of the water and air, when wave braking phenomena are modelled. Additionally, it might be noticed that change of the water height reconstructed with
Figure 10: Water height obtained with the CICSAM and the HRIC schemes (solid line), compared with experiment (black dots) and SPH method6(dashed-line - model for water-air system with compressibility included, doted line - model for water)
the HRIC compares better with experimental data than this obtained with the CICSAM. Comparison of the results from the H2 probe shows that both high-resolution schemes and the SPH6 only qualitatively predict behaviour of the interface, compare 5.6 < T < 7.5.
5.2 Comparison with pressure history from P1 probe
Pressure history obtained with the CICSAM and the HRIC schemes compares well with experimental evidence and results obtained with the SPH method6, see Fig. 11. It
might be observed that the moment of the first pressure peak, about T = 6 is better predicted with high-resolution schemes. The magnitude of the pressure peak is almost the same as in the case of SPH results6. Fast oscillations of the pressure from the SPH,
see Fig. 11, are explained as an effect of the air oscillations inside the bubble entrapped by the breaking wave. In the case of the solutions obtained with the CICSAM or the HRIC change in the pressure of an entrapped air is also modelled, however resulting bubble “oscillations” have much longer period. Additionally, one can notice that pressure history obtained with the CICSAM, is closer to the experimental evidence and the SPH6 results than in the case of the HRIC, see Fig. 11, T > 10.
Figure 11: Pressure history obtained with the CICSAM and the HRIC schemes (solid line) compared with experiment (solid dots) and results of the SPH method6 (dashed line - model for water-air system compressibility effects included)
5.3 Comparison of the wave profiles
Finally, differences in wave profiles predicted with both high-resolution schemes may be compared. Figure 12 depicts evolution of the free surface in non-dimensional time T . The difference between following time moments is constant and equal Tn+1− Tn = 1.2. The
first noticeable differences in the shape of the wave profiles between both high-resolution schemes are visible at T5 = 5.65 at the wave front. Wave profiles for times larger then
T5 are compared with results obtained using the LS (Level-Set) methods6,19, see Fig.
Figure 12: Free surface evolution in time obtained with the CICSAM (solid line) and the HRIC (dashed line) schemes, first differences in solution are noticable at T5
Differences appear in the volume and shape of the air bubble entrapped by the breaking wave and localisation and shape of the water filament that reflects after collision with the free surface.
6 CONCLUSIONS
In present paper, the VOF method with high-resolution scheme methodology was tested in two advection test cases and was applied to model dynamics of the system of immiscible fluids in the case of the breaking wave phenomenon. Two high-resolution schemes the CICSAM and the HRIC were employed for discretization of the nonlinear convective term in equation for transport of the volume fraction, cf. Eq. (5). Both high-resolution schemes proved their ability to capture interface position and preserve shape of the step interface profile.
Results obtained with the CICSAM and the HRIC schemes in two advection test cases show that both formulations are first order accurate and that their shape preserving or interface capturing properties are dependent on the size of the cell Courant number. Since the HRIC dependence on the cell Courant number is explicit, cf. Eq. (18), it seems to be less sensitive on the value of CD then the CICSAM. However, in opinion of the authors for
calculations with both considered high-resolution schemes Courant number value should be chosen smaller then CD < 0.5.
Comparison of both high-resolution schemes with the CIP18 type schemes, shows that
interface capturing properties of the CICSAM are close to the tangent-transformed CIP18. Taking into consideration complexity of both approaches, the first one seems to be prefer-able. However, this issue needs further investigation.
Finally, solution in the case of the breaking wave problem shows that it is possible to obtain results that compare well with an experimental evidence and other author results using the CICSAM and the HRIC high-resolution schemes.
REFERENCES
[1] S. Osher, R.P. Fedkiw. Level-Set Methods: An Overview and Some Recent Results, J. Comput. Phys., 169, 463–502, (2001).
[2] S. Osher, R. Fedkiw. Level Set Methods and Dynamic Implicit Surfaces, Springer, (2003).
[3] E. O˜nate, S.R. Idelsohn, F. Del Pin, R. Aubry. The particle finite element method. An overview. I. J. Comput. Meth., 1, 267–307, (2004).
[5] J.J. Monaghan. SPH without a Tensile Instability, J. Comput. Phys, 159, 290–311, (2000).
[6] A. Colagrossi, M. Landrini. Numerical simulation of interfacial flows by smoothed particle hydrodynamics, J. Comput. Phys, 191,448–475, (2003).
[7] C.W. Hirt, B.D. Nicholls. Volume of fluid method for dynamics of free boundaries, J. Comput. Phys., 39, 201–221, (1981).
[8] J.E. Pillod., E.G. Puckett. Second order-accurate volume-of-fluid algotithms for tracking material interfaces, J. Comput. Phys., 199, 718–742, (2004).
[9] J. Lopez, J. Hernandez, P. Gourez, F. Faura. A volume of fluid method based on multidimensional advection and spline interface reconstruction, J. Comput. Phys., 195, 718–742, (2004).
[10] Y. Renardy, , M. Renardy, V. Cristini. A new volume-of-fluid formulation for surfac-tants and simulations of drop deformation under shear at a low viscosity ratio. Eur. J. Mech. B - Fluids, 21, 49-59, (2002).
[11] D. Gerlach, G. Tomar, G. Biswas, F. Durst. Comparison of volume-of-fluid methods for surface tension-dominant two-phase flows. I. J. Heat and Mass Trans., 49, 740– 754, (2006).
[12] O. Ubbink, R.I. Issa. Method for Capturing Sharp Fluid Interfaces on Arbitrary Meshes, J. Comput. Phys., 153, 26–50, (1999).
[13] S. Muzaferija, M. Peric, P. Sames, T. Schelin. A two-fluid Navier-Stokes solver to simulate water entry, Twenty-Second Symposium on Naval Hydrodynamics, (1998). [14] B.P. Leonard. The ULTIMATE conservative difference scheme applied to unsteady
one-dimensional advection, Comput. Methods Appl. Mech. and Engng., 88, 17–74, (1991).
[15] J.H. Ferziger, M. Peric. Computational Methods for Fluid Dynamics, Springer, (2002).
[16] H.K. Versteeg, W. Malalasekera. An introduction to Computational Fluid Dynamics, Longman, (1998).
[17] T. Wac lawczyk, T. Koronowicz. Modeling of the flow in systems of immiscible fluids using Volume of Fluid method with CICSAM scheme, Turbulence, 8, 267–276, (2005). [18] T. Yabe, F. Xiao, T. Utsumi. The Constrained Interpolation Profile Method for
