- Email: [email protected]
Derivative of Heaviside step function vs. delta function in continuum surface force (CSF) models
ARTICLE IN PRESS
JID: IJMF
[m5G;December 14, 2017;5:23]
International Journal of Multiphase Flow 0 0 0 (2017) 1–4
Contents lists available at ScienceDirect
International Journal of Multiphase Flow journal homepage: www.elsevier.com/locate/ijmulflow
Short communication
Derivative of Heaviside step function vs. delta function in continuum surface force (CSF) models Hanif Montazeri a,b,1,∗, Seyed Hadi Zandavi b,1 a b
Depatment of Mechanical & Industrial Engineering, University of Toronto, 5 King’s College Road, Toronto, ON M5S 3G8, Canada NuPhysics Consulting Ltd., 112 College Street, Toronto, ON M5G 1L6, Canada
a r t i c l e
i n f o
a b s t r a c t
Article history: Received 28 June 2017 Revised 15 November 2017 Accepted 6 December 2017 Available online xxx
In Continuum Surface Force (CSF) model for implementing surface tension forces in multiphase flows, singular delta function is used to merge two continuous sets of flow equations. Discretizing the delta function has attracted a great deal of attention and has led to developing many different approaches. It has been shown numerically that two identical mathematical operators: Dirac delta function (δ ) and the gradient of Heaviside function (∇ H), lead to two different numerical solutions, and surprisingly the latter is the recommended method. Mathematically speaking, one would expect using the exact delta function prone to less numerical error. We explain the mathematical and mechanical reasons on why using the exact delta function leads to the reported numerical errors. © 2017 Elsevier Ltd. All rights reserved.
Keywords: Sharp interface Two-phase flow Dirac delta function Heaviside step function
Among all the complex phenomena associated with the numerical simulation of the multiphase flow, analyzing the interface has always been a great challenge as still numerous contradictions exist between experimental observations and the available theories and numerical analysis (Malladi et al., 1995; Scardovelli and Zaleski, 1999). In continuum-mechanics based formulations the interface between two phases is represented by a Dirac delta function (δ ). The singular delta function is used to merge two continuous sets of flow equations (i.e. Navier-Stokes equations). A suitable discretization of the delta function is the key in obtaining physical solutions. Discretizing the delta function has attracted a great deal of attention. The most commonly used methods for discretizing the Dirac delta function are the Continuum Surface Force (CSF) methods (Brackbill et al., 1992). Mathematically, the delta function (δ ) is defined as the gradient of the Heaviside function (∇ H):
δ = ∇H Therefore if we smear-out Heaviside function using the local values of level set function ϕ :
⎧ ⎪ ⎨
0 1 ϕ 1 πφ H (ϕ ) = + + sin ε ⎪ ⎩ 2 2ε 2π 1
ϕ < −ε −ε ≤ ϕ ≤ ε
ϕ>ε
∗
Corresponding author. E-mail addresses: [email protected] (H. Montazeri), [email protected] (S.H. Zandavi). 1 The authors contributed equally.
Then by taking derivative of H(ϕ ), we can obtain a smeared-out delta function:
⎧ ⎪ ⎨
0 1 1 πφ δ (ϕ ) = + sin ε ⎪ ⎩ 2ε 2ε 0
ϕ < −ε −ε ≤ ϕ ≤ ε
ϕ>ε
where ɛ determines the size of the bandwidth of numerical smearing. Since δ (ϕ ) = ∇ H (ϕ ), one would expect that ∇ H(ϕ ) and δ (ϕ ) can be used interchangeably. However, numerical tests show that solutions using δ (ϕ ) are not equivalent to numerical solutions using ∇ H(ϕ ) (Herrmann, 2008). The goal of this communication is to explain why these two identical operators could lead to two different numerical solutions. In Continuum Surface Force (CSF) models (Brackbill et al., 1992), surface forces are smeared over a few cells around the interface; therefore, a smeared version of the delta or Heaviside function must be used. Herrmann (2008) and Montazeri et al. (2012), Montazeri and Ward (2014) mentioned that if instead of using a delta function in two-phase flow equations:
D (ρ u ) = −∇ p + μ∇ 2 u + σ κδ (ϕ ) Dt
(1)
one uses its mathematical equivalent - the gradient of the Heaviside step function, ∇ H(ϕ ):
D (ρ u ) = −∇ p + μ∇ 2 u + σ κ∇ H (ϕ ) Dt
(2)
https://doi.org/10.1016/j.ijmultiphaseflow.2017.12.006 0301-9322/© 2017 Elsevier Ltd. All rights reserved.
Please cite this article as: H. Montazeri, S.H. Zandavi, Derivative of Heaviside step function vs. delta function in continuum surface force (CSF) models, International Journal of Multiphase Flow (2017), https://doi.org/10.1016/j.ijmultiphaseflow.2017.12.006
ARTICLE IN PRESS
JID: IJMF 2
[m5G;December 14, 2017;5:23]
H. Montazeri, S.H. Zandavi / International Journal of Multiphase Flow 000 (2017) 1–4
(a)
8
(b) 36.3 35 30 25 20 15 10 5 0
6
4
2
0
0
2
4
6
8
Fig. 1. Pressure contours for a static bubble which is simulated using a) Direct delta function (δ ) where spurious velocities as large as 2.34 × 10−6 mm/s occur around the interface, b) the gradient of the Heaviside step function, ∇ H, where no spurious velocity occur. The exact curvature is used in both simulations. Side lengths are in millimeter.
Areal method conversion
(a)
Stair-step Conversion
(b) Fig. 2. Comparing how smeared Delta function (δ ) and smeared Heaviside ∇ H numerically represent an interface: (a) smeared delta function (δ ) would covert the physical interface more of inclined small lines within each control volume, just like in Areal method, (b) Smeared Heaviside ∇ H converts an interface into stair-step-like interfaces, just like in stair-step method.
then no spurious velocities would occur around an interface (provided that the exact curvature is known)! Here, u, p, ρ and μ are the velocity, pressure, density and viscosity of the flow; and σ and κ are the surface tension and curvature of an interface. However, there is no explanation on why using the exact values of the delta function results in erroneous spurious velocities. Note that, one expects more accurate solutions from directly using the smeared delta function than using the smeared Heaviside function. This is because when using the smeared Heaviside function, a numerical derivative of the Heaviside function is required to estimate ∇ H. Therefore, the numerical estimation of ∇ H involves truncation errors; whereas, there are no truncation errors involved in calculating the smeared delta function. Surprisingly, using ∇ H results in more accurate solutions compared to using δ in CSF models!
An example of an inaccurate solution obtained from Eq. (1) can be seen in a virtual simulation of a static bubble with a radius R = 2 × 10−3 m which is positioned in a cubical domain having a side length of eight millimeters (Fig. 1). The surface tension is Kg N taken to be 72 × 10−3 m , and the density inside the bubble is 1 m 3,
Kg while the outside density is 10 0 0 m 3 (Francois et al., 2006). Spurious velocities occur around the entire interface (Fig. 1a), and consequently an, inaccurate pressure jump of 36.31 pa is predicted when Eq. (1) has been solved. As shown by Herrmann (2008) and Montazeri et al. (2012), if Eq. (2) is used, the exact pressure jump is obtained without any spurious velocities (Fig. 1b). Montazeri et al. (2017) studied two different methods for the sharp implementation of interfaces: stair-step and Areal method. They presented a detailed study that shows in a Cartesian grid that stair-step methods can result in consistent numerical solu-
Please cite this article as: H. Montazeri, S.H. Zandavi, Derivative of Heaviside step function vs. delta function in continuum surface force (CSF) models, International Journal of Multiphase Flow (2017), https://doi.org/10.1016/j.ijmultiphaseflow.2017.12.006
JID: IJMF
ARTICLE IN PRESS
[m5G;December 14, 2017;5:23]
H. Montazeri, S.H. Zandavi / International Journal of Multiphase Flow 000 (2017) 1–4
3
ues for the delta function (which are exact up to machine errors):
fe = σ κδes fn = σ κδns
Fig. 3. Schematic of a control volume and an interface. The solid line is the interface and the dashed line is an iso-surface of its delta function.
tions. Methods such as the Areal method, which tends to preserve the original shape of an interface, encounter misplacing surface forces and therefore produce inconsistent solutions. Considering the numerical experiments and mechanical explanation in Montazeri et al. (2017), it is now possible to explain why using δ results in erroneous solutions. Due to the accuracy of the delta function (similar to the Areal Method in their paper), the smeared delta function smears out the interface parallel to the actual interface (Fig. 2a). Therefore the physical shape of the interface is preserved (as in the Areal method). In contrast, the derivative of the Heaviside step function treats the interface more like that of a stair-step-like interface (as in the stair-step method, Fig. 2b); therefore the physical interface would be converted into stairsteps like interfaces. The reader is advised to study Section 5 of Montazeri et al. (2017) paper. The numerical mechanisms of these two different methods are visually represented in Figs. 3 and 4. In Fig. 3, the dashed line is one of the iso-surfaces of a smeared delta function. Since the normal distance from the iso-surface to the interface is the same for all the points along the interface, the value of the smeared delta function is the same for all the points along this iso-surface. As seen in Fig. 3, this iso-surface does not pass through point n, as the normal distance of point n is not equivalent to the normal distance of point e. Consequently, when interfacial forces are estimated at points e and n, each force is estimated by different non-zero val-
Here f is the surface force and the superscript s denotes the smeared quantity, i.e. smeared delta function. Details on how to smear out the delta function are given by Osher and Fedkiw (2004). However, when the derivative of the step function is used, to estimate the surface forces at e and n, the values of the smeared Heaviside step function of points E, P and N are used:
fe = σ κ
HEs − HPs
x
fn = σ κ
HNs − HPs
y
Since a first order-finite-difference is used, the accuracy is only first order accurate. This less accurate estimation maps the interface very similar to stair-step method. For instance, for the interface shown in Fig. 3, if a non-smeared H is used, fn becomes zero, and fe takes the entire value of the surface force which is σ κ , which is identical to the stair-step method. More details about the stair-step methods can be found in Montazeri et al. (2017). Fig. 4 summarizes the above discussion in a more graphical way. Using the delta function, the interfacial force is perpendicular to the converted interface i.e. the straight line in the control volume. In Newtonian mechanics, a force can be moved along its direction; therefore, we can transfer the interfacial force until it intersects with the control volume face. Hence, the X-component of the force is retained as the force applied on this face. As shown in Fig. 4, for this particular scenario, the interfacial force does not pass through the center of the face. However, in discretizing the momentum equation, one must consider the X-component of the interfacial force that is passing through the center of the staggered control volume. The displacement from another location to the face center is permitted only if the corresponding moment is also included in the flow equations. This, however, is not the case when delta function is used. As a result of ignoring these moments, the spurious velocities occur when this method is used.
Fig. 4. Estimating the interfacial force for a face of a control volume, shown on top, using two different methods: Gradient of Heaviside function (left) and Direct delta function (right).
Please cite this article as: H. Montazeri, S.H. Zandavi, Derivative of Heaviside step function vs. delta function in continuum surface force (CSF) models, International Journal of Multiphase Flow (2017), https://doi.org/10.1016/j.ijmultiphaseflow.2017.12.006
ARTICLE IN PRESS
JID: IJMF 4
[m5G;December 14, 2017;5:23]
H. Montazeri, S.H. Zandavi / International Journal of Multiphase Flow 000 (2017) 1–4
4
4
2
2
0
0
2
4
0 (a)
2
0
0
2
4
0
2
4
2
0
2
4
0 (b)
Fig. 5. Application of smeared delta function in the formulation of CSF. (a) A 45° interface resolved on uniform mesh; (b) A 30° interface is resolved on uniform mesh.
The above explanation can be numerically tested in an idealized numerical experiment which was used by Montazeri et al. (2017). A flat interface with a constant surface tension pressure of σ κ = 10 pa was simulated using the delta function, Eq. (1). We considered two angles for the flat interface: 45 and 30° (Fig. 5). For the 45-degree interface, surprisingly no spurious velocities occurred around the interface, as shown in case (a) of Fig. 5. This can be explained using Fig. 3, where for the 45° interface the very same iso-surface which passes thorough e, would pass through point n; and therefore the value of δ is the same for both n and e. However, this is not the case for any other angle such as a 30° flat interface, as shown in Fig. 5b, where spurious velocities occurred around the interface (similar to the results of Areal method). The above in-depth analysis explains why the numerical observation by Herrmann (2008) and Montazeri et al. (2012) makes mathematical sense: the direct discretization of smeared δ results in inaccurate solutions when the staggered formulation of Francois et al. (2006) is used. Acknowledgments
References Brackbill, J.U., Kothe, D.B., Zemach, C., 1992. A continuum method for modeling surface tension. J. Comput. Phys. 100, 335–354. doi:10.1016/0021- 9991(92)90240- Y. Francois, M.M., Cummins, S.J., Dendy, E.D., Kothe, D.B., Sicilian, J.M., Williams, M.W., 2006. A balanced-force algorithm for continuous and sharp interfacial surface tension models within a volume tracking framework. J. Comput. Phys. 213, 141– 173. doi:10.1016/j.jcp.20 05.08.0 04. Herrmann, M., 2008. A balanced force refined level set grid method for twophase flows on unstructured flow solver grids. J. Comput. Phys. 227, 2674–2706. doi:10.1016/j.jcp.20 07.11.0 02. Malladi, R., Sethian, J.A., Vemuri, B.C.others, 1995. Shape modeling with front propagation: a level set approach. IEEE Trans. Pattern Anal. Mach. Intell. 17, 158–175. doi:10.1109/34.368173. Montazeri, H., Bussmann, M., Mostaghimi, J., 2012. Accurate implementation of forcing terms for two-phase flows into SIMPLE algorithm. Int. J. Multiph. Flow 45, 40–52. doi:10.1016/j.ijmultiphaseflow.2012.05.003. Montazeri, H., Ward, C.A., 2014. A balanced-force algorithm for two-phase flows. J. Comput. Phys. 257, 645–669. doi:10.1016/j.jcp.2013.09.054. Montazeri, H., Zandavi, S.H., Bazylak, A., 2017. Sharp interface models for two-phase flows: insights towards new approaches. Comput. Methods Appl. Mech. Eng. 322, 238–261. doi:10.1016/j.cma.2017.04.022. Osher, S., Fedkiw, R., 2004. Level set methods and dynamic implicit surfaces. Appl. Mech. Rev. doi:10.1115/1.1760520. Scardovelli, R., Zaleski, S., 1999. Direct numerical simulation of free-surface and interfacial flow. Annu. Rev. Fluid Mech. 31, 567–603. doi:10.1146/annurev.fluid.31. 1.567.
The authors would like to acknowledge the financial contributions to this work from Grand Challenges Canada (#0501-01-10 and R-ST-POC-1707-07082).
Please cite this article as: H. Montazeri, S.H. Zandavi, Derivative of Heaviside step function vs. delta function in continuum surface force (CSF) models, International Journal of Multiphase Flow (2017), https://doi.org/10.1016/j.ijmultiphaseflow.2017.12.006
JID: IJMF
[m5G;December 14, 2017;5:23]
International Journal of Multiphase Flow 0 0 0 (2017) 1–4
Contents lists available at ScienceDirect
International Journal of Multiphase Flow journal homepage: www.elsevier.com/locate/ijmulflow
Short communication
Derivative of Heaviside step function vs. delta function in continuum surface force (CSF) models Hanif Montazeri a,b,1,∗, Seyed Hadi Zandavi b,1 a b
Depatment of Mechanical & Industrial Engineering, University of Toronto, 5 King’s College Road, Toronto, ON M5S 3G8, Canada NuPhysics Consulting Ltd., 112 College Street, Toronto, ON M5G 1L6, Canada
a r t i c l e
i n f o
a b s t r a c t
Article history: Received 28 June 2017 Revised 15 November 2017 Accepted 6 December 2017 Available online xxx
In Continuum Surface Force (CSF) model for implementing surface tension forces in multiphase flows, singular delta function is used to merge two continuous sets of flow equations. Discretizing the delta function has attracted a great deal of attention and has led to developing many different approaches. It has been shown numerically that two identical mathematical operators: Dirac delta function (δ ) and the gradient of Heaviside function (∇ H), lead to two different numerical solutions, and surprisingly the latter is the recommended method. Mathematically speaking, one would expect using the exact delta function prone to less numerical error. We explain the mathematical and mechanical reasons on why using the exact delta function leads to the reported numerical errors. © 2017 Elsevier Ltd. All rights reserved.
Keywords: Sharp interface Two-phase flow Dirac delta function Heaviside step function
Among all the complex phenomena associated with the numerical simulation of the multiphase flow, analyzing the interface has always been a great challenge as still numerous contradictions exist between experimental observations and the available theories and numerical analysis (Malladi et al., 1995; Scardovelli and Zaleski, 1999). In continuum-mechanics based formulations the interface between two phases is represented by a Dirac delta function (δ ). The singular delta function is used to merge two continuous sets of flow equations (i.e. Navier-Stokes equations). A suitable discretization of the delta function is the key in obtaining physical solutions. Discretizing the delta function has attracted a great deal of attention. The most commonly used methods for discretizing the Dirac delta function are the Continuum Surface Force (CSF) methods (Brackbill et al., 1992). Mathematically, the delta function (δ ) is defined as the gradient of the Heaviside function (∇ H):
δ = ∇H Therefore if we smear-out Heaviside function using the local values of level set function ϕ :
⎧ ⎪ ⎨
0 1 ϕ 1 πφ H (ϕ ) = + + sin ε ⎪ ⎩ 2 2ε 2π 1
ϕ < −ε −ε ≤ ϕ ≤ ε
ϕ>ε
∗
Corresponding author. E-mail addresses: [email protected] (H. Montazeri), [email protected] (S.H. Zandavi). 1 The authors contributed equally.
Then by taking derivative of H(ϕ ), we can obtain a smeared-out delta function:
⎧ ⎪ ⎨
0 1 1 πφ δ (ϕ ) = + sin ε ⎪ ⎩ 2ε 2ε 0
ϕ < −ε −ε ≤ ϕ ≤ ε
ϕ>ε
where ɛ determines the size of the bandwidth of numerical smearing. Since δ (ϕ ) = ∇ H (ϕ ), one would expect that ∇ H(ϕ ) and δ (ϕ ) can be used interchangeably. However, numerical tests show that solutions using δ (ϕ ) are not equivalent to numerical solutions using ∇ H(ϕ ) (Herrmann, 2008). The goal of this communication is to explain why these two identical operators could lead to two different numerical solutions. In Continuum Surface Force (CSF) models (Brackbill et al., 1992), surface forces are smeared over a few cells around the interface; therefore, a smeared version of the delta or Heaviside function must be used. Herrmann (2008) and Montazeri et al. (2012), Montazeri and Ward (2014) mentioned that if instead of using a delta function in two-phase flow equations:
D (ρ u ) = −∇ p + μ∇ 2 u + σ κδ (ϕ ) Dt
(1)
one uses its mathematical equivalent - the gradient of the Heaviside step function, ∇ H(ϕ ):
D (ρ u ) = −∇ p + μ∇ 2 u + σ κ∇ H (ϕ ) Dt
(2)
https://doi.org/10.1016/j.ijmultiphaseflow.2017.12.006 0301-9322/© 2017 Elsevier Ltd. All rights reserved.
Please cite this article as: H. Montazeri, S.H. Zandavi, Derivative of Heaviside step function vs. delta function in continuum surface force (CSF) models, International Journal of Multiphase Flow (2017), https://doi.org/10.1016/j.ijmultiphaseflow.2017.12.006
ARTICLE IN PRESS
JID: IJMF 2
[m5G;December 14, 2017;5:23]
H. Montazeri, S.H. Zandavi / International Journal of Multiphase Flow 000 (2017) 1–4
(a)
8
(b) 36.3 35 30 25 20 15 10 5 0
6
4
2
0
0
2
4
6
8
Fig. 1. Pressure contours for a static bubble which is simulated using a) Direct delta function (δ ) where spurious velocities as large as 2.34 × 10−6 mm/s occur around the interface, b) the gradient of the Heaviside step function, ∇ H, where no spurious velocity occur. The exact curvature is used in both simulations. Side lengths are in millimeter.
Areal method conversion
(a)
Stair-step Conversion
(b) Fig. 2. Comparing how smeared Delta function (δ ) and smeared Heaviside ∇ H numerically represent an interface: (a) smeared delta function (δ ) would covert the physical interface more of inclined small lines within each control volume, just like in Areal method, (b) Smeared Heaviside ∇ H converts an interface into stair-step-like interfaces, just like in stair-step method.
then no spurious velocities would occur around an interface (provided that the exact curvature is known)! Here, u, p, ρ and μ are the velocity, pressure, density and viscosity of the flow; and σ and κ are the surface tension and curvature of an interface. However, there is no explanation on why using the exact values of the delta function results in erroneous spurious velocities. Note that, one expects more accurate solutions from directly using the smeared delta function than using the smeared Heaviside function. This is because when using the smeared Heaviside function, a numerical derivative of the Heaviside function is required to estimate ∇ H. Therefore, the numerical estimation of ∇ H involves truncation errors; whereas, there are no truncation errors involved in calculating the smeared delta function. Surprisingly, using ∇ H results in more accurate solutions compared to using δ in CSF models!
An example of an inaccurate solution obtained from Eq. (1) can be seen in a virtual simulation of a static bubble with a radius R = 2 × 10−3 m which is positioned in a cubical domain having a side length of eight millimeters (Fig. 1). The surface tension is Kg N taken to be 72 × 10−3 m , and the density inside the bubble is 1 m 3,
Kg while the outside density is 10 0 0 m 3 (Francois et al., 2006). Spurious velocities occur around the entire interface (Fig. 1a), and consequently an, inaccurate pressure jump of 36.31 pa is predicted when Eq. (1) has been solved. As shown by Herrmann (2008) and Montazeri et al. (2012), if Eq. (2) is used, the exact pressure jump is obtained without any spurious velocities (Fig. 1b). Montazeri et al. (2017) studied two different methods for the sharp implementation of interfaces: stair-step and Areal method. They presented a detailed study that shows in a Cartesian grid that stair-step methods can result in consistent numerical solu-
Please cite this article as: H. Montazeri, S.H. Zandavi, Derivative of Heaviside step function vs. delta function in continuum surface force (CSF) models, International Journal of Multiphase Flow (2017), https://doi.org/10.1016/j.ijmultiphaseflow.2017.12.006
JID: IJMF
ARTICLE IN PRESS
[m5G;December 14, 2017;5:23]
H. Montazeri, S.H. Zandavi / International Journal of Multiphase Flow 000 (2017) 1–4
3
ues for the delta function (which are exact up to machine errors):
fe = σ κδes fn = σ κδns
Fig. 3. Schematic of a control volume and an interface. The solid line is the interface and the dashed line is an iso-surface of its delta function.
tions. Methods such as the Areal method, which tends to preserve the original shape of an interface, encounter misplacing surface forces and therefore produce inconsistent solutions. Considering the numerical experiments and mechanical explanation in Montazeri et al. (2017), it is now possible to explain why using δ results in erroneous solutions. Due to the accuracy of the delta function (similar to the Areal Method in their paper), the smeared delta function smears out the interface parallel to the actual interface (Fig. 2a). Therefore the physical shape of the interface is preserved (as in the Areal method). In contrast, the derivative of the Heaviside step function treats the interface more like that of a stair-step-like interface (as in the stair-step method, Fig. 2b); therefore the physical interface would be converted into stairsteps like interfaces. The reader is advised to study Section 5 of Montazeri et al. (2017) paper. The numerical mechanisms of these two different methods are visually represented in Figs. 3 and 4. In Fig. 3, the dashed line is one of the iso-surfaces of a smeared delta function. Since the normal distance from the iso-surface to the interface is the same for all the points along the interface, the value of the smeared delta function is the same for all the points along this iso-surface. As seen in Fig. 3, this iso-surface does not pass through point n, as the normal distance of point n is not equivalent to the normal distance of point e. Consequently, when interfacial forces are estimated at points e and n, each force is estimated by different non-zero val-
Here f is the surface force and the superscript s denotes the smeared quantity, i.e. smeared delta function. Details on how to smear out the delta function are given by Osher and Fedkiw (2004). However, when the derivative of the step function is used, to estimate the surface forces at e and n, the values of the smeared Heaviside step function of points E, P and N are used:
fe = σ κ
HEs − HPs
x
fn = σ κ
HNs − HPs
y
Since a first order-finite-difference is used, the accuracy is only first order accurate. This less accurate estimation maps the interface very similar to stair-step method. For instance, for the interface shown in Fig. 3, if a non-smeared H is used, fn becomes zero, and fe takes the entire value of the surface force which is σ κ , which is identical to the stair-step method. More details about the stair-step methods can be found in Montazeri et al. (2017). Fig. 4 summarizes the above discussion in a more graphical way. Using the delta function, the interfacial force is perpendicular to the converted interface i.e. the straight line in the control volume. In Newtonian mechanics, a force can be moved along its direction; therefore, we can transfer the interfacial force until it intersects with the control volume face. Hence, the X-component of the force is retained as the force applied on this face. As shown in Fig. 4, for this particular scenario, the interfacial force does not pass through the center of the face. However, in discretizing the momentum equation, one must consider the X-component of the interfacial force that is passing through the center of the staggered control volume. The displacement from another location to the face center is permitted only if the corresponding moment is also included in the flow equations. This, however, is not the case when delta function is used. As a result of ignoring these moments, the spurious velocities occur when this method is used.
Fig. 4. Estimating the interfacial force for a face of a control volume, shown on top, using two different methods: Gradient of Heaviside function (left) and Direct delta function (right).
Please cite this article as: H. Montazeri, S.H. Zandavi, Derivative of Heaviside step function vs. delta function in continuum surface force (CSF) models, International Journal of Multiphase Flow (2017), https://doi.org/10.1016/j.ijmultiphaseflow.2017.12.006
ARTICLE IN PRESS
JID: IJMF 4
[m5G;December 14, 2017;5:23]
H. Montazeri, S.H. Zandavi / International Journal of Multiphase Flow 000 (2017) 1–4
4
4
2
2
0
0
2
4
0 (a)
2
0
0
2
4
0
2
4
2
0
2
4
0 (b)
Fig. 5. Application of smeared delta function in the formulation of CSF. (a) A 45° interface resolved on uniform mesh; (b) A 30° interface is resolved on uniform mesh.
The above explanation can be numerically tested in an idealized numerical experiment which was used by Montazeri et al. (2017). A flat interface with a constant surface tension pressure of σ κ = 10 pa was simulated using the delta function, Eq. (1). We considered two angles for the flat interface: 45 and 30° (Fig. 5). For the 45-degree interface, surprisingly no spurious velocities occurred around the interface, as shown in case (a) of Fig. 5. This can be explained using Fig. 3, where for the 45° interface the very same iso-surface which passes thorough e, would pass through point n; and therefore the value of δ is the same for both n and e. However, this is not the case for any other angle such as a 30° flat interface, as shown in Fig. 5b, where spurious velocities occurred around the interface (similar to the results of Areal method). The above in-depth analysis explains why the numerical observation by Herrmann (2008) and Montazeri et al. (2012) makes mathematical sense: the direct discretization of smeared δ results in inaccurate solutions when the staggered formulation of Francois et al. (2006) is used. Acknowledgments
References Brackbill, J.U., Kothe, D.B., Zemach, C., 1992. A continuum method for modeling surface tension. J. Comput. Phys. 100, 335–354. doi:10.1016/0021- 9991(92)90240- Y. Francois, M.M., Cummins, S.J., Dendy, E.D., Kothe, D.B., Sicilian, J.M., Williams, M.W., 2006. A balanced-force algorithm for continuous and sharp interfacial surface tension models within a volume tracking framework. J. Comput. Phys. 213, 141– 173. doi:10.1016/j.jcp.20 05.08.0 04. Herrmann, M., 2008. A balanced force refined level set grid method for twophase flows on unstructured flow solver grids. J. Comput. Phys. 227, 2674–2706. doi:10.1016/j.jcp.20 07.11.0 02. Malladi, R., Sethian, J.A., Vemuri, B.C.others, 1995. Shape modeling with front propagation: a level set approach. IEEE Trans. Pattern Anal. Mach. Intell. 17, 158–175. doi:10.1109/34.368173. Montazeri, H., Bussmann, M., Mostaghimi, J., 2012. Accurate implementation of forcing terms for two-phase flows into SIMPLE algorithm. Int. J. Multiph. Flow 45, 40–52. doi:10.1016/j.ijmultiphaseflow.2012.05.003. Montazeri, H., Ward, C.A., 2014. A balanced-force algorithm for two-phase flows. J. Comput. Phys. 257, 645–669. doi:10.1016/j.jcp.2013.09.054. Montazeri, H., Zandavi, S.H., Bazylak, A., 2017. Sharp interface models for two-phase flows: insights towards new approaches. Comput. Methods Appl. Mech. Eng. 322, 238–261. doi:10.1016/j.cma.2017.04.022. Osher, S., Fedkiw, R., 2004. Level set methods and dynamic implicit surfaces. Appl. Mech. Rev. doi:10.1115/1.1760520. Scardovelli, R., Zaleski, S., 1999. Direct numerical simulation of free-surface and interfacial flow. Annu. Rev. Fluid Mech. 31, 567–603. doi:10.1146/annurev.fluid.31. 1.567.
The authors would like to acknowledge the financial contributions to this work from Grand Challenges Canada (#0501-01-10 and R-ST-POC-1707-07082).
Please cite this article as: H. Montazeri, S.H. Zandavi, Derivative of Heaviside step function vs. delta function in continuum surface force (CSF) models, International Journal of Multiphase Flow (2017), https://doi.org/10.1016/j.ijmultiphaseflow.2017.12.006