A hybrid Taylor–Galerkin variational multi-scale stabilization method for the level set equation

Computers and Fluids 121 (2015) 192–205

Contents lists available at ScienceDirect

Computers and Fluids journal homepage: www.elsevier.com/locate/compfluid

A hybrid Taylor–Galerkin variational multi-scale stabilization method for the level set equation Ahmed Bakkar a,∗, Wagdi G. Habashi a,b, Marco Fossati a, Guido S. Baruzzi b a b

Department of Mechanical Engineering, CFD Laboratory, McGill University, 688 Sherbrooke Street West, Montreal, QC H3A 2S6, Canada Newmerical Technologies International, 680 Sherbrooke Street West, Montreal, QC H3A 2M7, Canada

a r t i c l e

i n f o

Article history: Received 13 June 2014 Revised 29 January 2015 Accepted 11 August 2015 Available online 18 August 2015 Keywords: Level set Variational multi-scale Taylor Galerkin Multiphase flow

a b s t r a c t A stabilized finite element formulation of the level set equation is proposed for the numerical simulation of water droplet dynamics for in-flight ice accretion problems. The variational multi-scale and Taylor–Galerkin approaches are coupled such that the temporal derivative in the weak Galerkin formulation is replaced with a Taylor series expansion improving the temporal accuracy of the scheme. The variational multi-scale approach is then applied to the semi-discrete equation, allowing the stabilization terms to appear naturally. Taylor series expansions up to the fourth order have been studied in terms of accuracy and convergence rates. A second order implicit expansion was found to provide a good trade-off between accuracy and computational cost when compared to a fourth order implicit expansion. Validation is done through a number of benchmark cases considering droplet stretching and high-speed advection. Results indicate good mass conservation characteristics compared to other methods available in the literature. © 2015 Elsevier Ltd. All rights reserved.

1. Introduction Understanding and modeling the impingement dynamics of supercooled large droplets (SLDs) is crucial to be able to address accurately the mechanism of in-flight icing on aircraft and rotorcraft. This has become critical after the introduction of the new Appendix O for certification [1]. SLD are droplets of diameters greater than 50 μm which are found in the atmosphere in liquid state at temperatures below the freezing point as a consequence of the surface tension [1]. While small droplets are generally assumed to retain a nearly spherical shape even under substantial aerodynamic stresses, this is not true for SLD, especially at speeds typical of aerospace applications. In these conditions, droplets undergo great deformations and the phenomenology of the impact with aircraft surfaces is fairly complex. Aerodynamic distortion, break-up before impact, splashing and bouncing are phenomena that have to be taken into consideration when studying SLD dynamics. These phenomena could lead to water deposition over areas that would have remained dry in the case of the impact of smaller droplets. Failing to take the SLD behavior into account when designing ice protection systems may lead to the formation of ice on unprotected surfaces—a safety disaster. The problem of SLD dynamics can be investigated either experimentally or numerically. Studying SLD experimentally in icing

∗

Corresponding author. Tel.: +15143981710. E-mail address: [email protected] (A. Bakkar).

http://dx.doi.org/10.1016/j.compfluid.2015.08.008 0045-7930/© 2015 Elsevier Ltd. All rights reserved.

tunnels is not straightforward due to technical limitations that can affect the efficacy of the study. Atmospheric droplet distributions can be very difficult to reproduce, particularly those proposed in the new certification regulations. Replicating the flow conditions experienced in an in-flight icing scenario is also a challenging and expensive task, especially when scaling is involved [2]. Accurate surface impingement data are also nearly impossible to obtain in a non-invasive manner. Numerical methods, on the other hand, are capable of reproducing the flow conditions without the need for scaling and therefore can be considered as a valid instrument for this type of studies. A major challenge to numerical modeling of SLD is the ability to address the evolution of the surface of the droplets due to the aerodynamic shear that makes it highly irregular and prone to folding, pinching, merging and eventually break up. Large droplets impacting a surface at high speed can shatter and produce sprays of smaller droplets that re-enter the flow and deposit further downstream, by-passing the ice protection system [32,33]. Many reports can be found in the literature that propose numerical approaches to address the dynamics of droplets in in-flight icing. Nevertheless the majority of these methods make use of heuristic correlations based on single or multiple droplet experiments, or simplified analytical approaches that do not consider the SLD regime explicitly, or extrapolate the SLD behavior from non-SLD/non-aeronautical conditions and eventually may fail to provide accurate predictions [3,4]. Better insight in SLD dynamics could be obtained by considering numerical approaches to consistently study the detailed surface dynamics and evolution of a single and/or multiple SLD to remove some of the

A. Bakkar et al. / Computers and Fluids 121 (2015) 192–205

empiricism currently adopted in macroscopic analyses of droplet impact. Single-droplet dynamics can be modeled with a Lagrangian, Eulerian, or mixed Eulerian–Lagrangian approach [5]. In a Lagrangian approach (also referred to as interface tracking), the grid coincides with the interface and tracking is explicitly accomplished by moving the grid points along with the interface. This method is generally accurate for cases involving small interface motion, but it is not effective when large deformations occur [5]. Lagrangian methods are usually computationally expensive due to the costs associated with a mesh that has to conform to the new interface location at every time step. In the mixed Eulerian–Lagrangian approach, the computations are done on a fixed Eulerian grid with marker points are used to define precisely the interface. These marker points are then used to guide and correct the interface location providing improved accuracy compared to a pure Eulerian approach. The major drawback of this approach is the difficulty associated with the tracking of the marker points when large deformations and topological changes occur, such as the merging of two droplets or the creation of non-simply connected regions [5]. The computational cost of mixed methods is reduced with respect to the Lagrangian approach, but still remains very high, especially when many markers are needed to accurately represent the interface. Lagrangian and mixed methods enjoy favorable properties when it comes to mass conservation since the tracking of the interface prevents any mass loss or gain, but their computational cost makes these methods unattractive. The Eulerian approach (also referred to as interface capturing) captures the interface on a fixed grid through a scalar function containing the interface information. The temporal and spatial evolution of the interface is updated by solving a transport equation. The Eulerian method offers: a greatly reduced computational cost and the possibility to define a consistent and unified numerical method that accounts simultaneously and accurately for the interaction between the different phases by resorting to mesh adaptation techniques. Despite the problems related to numerical stability and mass conservation, the Eulerian approach is recognized as a good candidate to address problems that involve the large deformation and topology changes of SLD dynamics. The two most commonly adopted Eulerian methods are the level set method (LSM) [6,7] and the volume of fluid (VOF) method [8,9]. Both methods use the same advection equation but differ in the scalar function adopted for the advection problem. The LSM uses a signed distance function while the VOF method uses the volume fraction. The LSM, introduced by Osher and Sethian [10], provides an efficient and simple approach for modeling interface motion. Unlike the VOF method where the interface is reconstructed using volume fractions, the LSM provides a continuous interface representation allowing for the accurate account of the surface tension forces. Computing the surface normal and curvature is also straightforward due to the use of a signed distance function. These features make the LSM an attractive choice for modeling high-speed droplet impingement. The LSM uses an implicit representation of the interface through the zero level set (LS) of a signed distance function. This function is evolved in time using a transport equation [6,11]. In the present work, a finite element (FE) implementation of the LSM was chosen due to its favorable numerical characteristics. These include the superior accuracy of the solution, the rigorous treatment of the boundary conditions, and the ability to handle anisotropic meshes with ease. The FE implementation of the LSM, however suffers from spurious oscillations for advection-dominated problems, which become even more noticeable in the case of the high speed flows typical of aeronautical applications. Stabilization techniques either introduce artificial diffusion, or an upwind discretization of the convection term. The addition of artificial diffusion offsets the negative diffusion introduced by the weak-Galerkin formulation, stabilizing the method, while the upwind formulation avoids the problems typical of the central difference approximation arising from the standard weak-Galerkin formulation

193

[12]. Various FE implementations of the LSM have been published in the literature and some of the most relevant work will be outlined. Nagrath et al. [13] introduced a streamline-upwind/Petrov–Galerkin (SUPG) stabilized FE LSM for incompressible two-phase flows, which was used to model a bubble rising in a liquid. Valance et al. [14] developed a Galerkin least-squares stabilized Bubnov–Galerkin formulation for the LSM for applications involving irregular domains and discontinuities. Cho et al. [15] developed a direct re-initialization scheme for a Taylor–Galerkin (TG) stabilized FE-LSM for incompressible two-phase flows. The method was subsequently refined to include flows with surface tension [16]. Farthing and Kees [17] introduced two different stabilization methods for the LS equation. The first employed a Runge–Kutta Discontinuous Galerkin method. It was shown that, in some circumstances, this lead to entropy-violating solutions which are overcome by the addition of shock-capturing diffusion and viscous stabilization. The second approach employed a variational multi-scale (VMS) continuous Galerkin formulation, which was also found to provide inaccurate results for some test cases [17] and required the introduction of isentropic shock capturing terms. Kees et al. [18] later improved the mass conservation of the schemes by coupling them with the volume fraction equation. In this paper, stabilization of the FE-LSM through VMS is further investigated. A hybrid TG-VMS approach is proposed in the attempt to enhance stabilization through the introduction of a Taylor series expansion for the temporal term. Temporal accuracy will also be improved with a higher-order temporal discretization, while the VMS approach will introduce the stabilization terms naturally. The outline of the paper is as follows: Section 2 introduces the equations to be modeled, Section 3 introduces the hybrid approach applied to the LS equation, Section 4 presents the numerical results and Section 5 contains the conclusions. 2. Mathematical model In the LSM, the interface is represented by the zero LS of a scalar function ϕ . The signed distance function is the most commonly used and it is defined as the minimum distance between a grid node and the interface. The scalar distance function ϕ is then advected according to the following partial differential equation (PDE):

ϕt + uˆ · ∇ϕ = 0

(1)

where uˆ is the interface advection velocity, and ϕ is the signed distance function. By virtue of Eq. (1), the second order temporal derivative of the LS equation can be written as:

ϕtt = −ut · ∇ϕ + (uˆ · ∇(uˆ · ∇ϕ))

(2)

where uˆ and ut are, respectively, the average advection velocity and its first order temporal derivative, which are given by:

uˆ =

un+1 + un 2

and ut =

un+1 − un t

(3)

where t is the time step used. The surface normal n and curvature k of the interface can be calculated using:

n=

∇ϕ ∇ϕ

k=∇ ·n=∇ ·

(4)

∇ϕ ∇ϕ

(5)

When using a signed distance function, the L2-norm of the gradient of the scalar function ϕ is equal to unity, simplifying Eqs. (4) and (5). However as Eq. (1) advects the interface with time, the function ϕ will drift from being a signed distance function and may become irregular [7]. This is because the velocity field will not necessarily transport all LSs with the same velocity [19]. This can lead to incorrect motion of the interface and may result in additional mass loss [20].

194

A. Bakkar et al. / Computers and Fluids 121 (2015) 192–205

To avoid such problems, the LS function is returned to a signed distance function through a re-initialization process. Re-initialization is essential when the velocity field depends on the thickness of the interface, such as in the case of modeling droplet dynamics, the curvature or surface normal, as well as when using a narrow-band LS approach. Re-initialization can be achieved either by explicitly defining the interface after a number of time steps and computing the new values of the scalar function, or by solving a re-initialization PDE in pseudo-time until a steady-state solution is reached. The reinitialization equation is written as:

ϕt  + s(ϕ)(∇ϕ − 1) = 0

Table 1 Values of the constants used in Eqs. (16) and (18).

s(ϕ) =



β

γ

δ

1/2

1

−1/3

1/12

where ϕ is the function, k is the desired order of expansion, and t is time. fi−1 (ϕ) expanding (12) up to the fourth derivative:

(6)

(7)

ϕ 2 + ∇ϕ2 (x)2

t 3 ∂ 3 ϕ 3! ∂ t 3

n+1

t

ϕ n = ϕ n+1 −

1!

ϕtn+1 +

t 2 2!

ϕttn+1 −

t 3 3!

n+1 ϕttt +

t 4 4!

n+1 ϕtttt + ···

(13)

ϕt  + c · ∇ϕ = s(ϕ)

(8)

where c represents the sign function pointing in the normal direction to the surface and is given by:

∇ϕ c = s(ϕ)n = s(ϕ) ∇ϕ

(9)

Eq. (13) can be re-written as: - First order

ϕtn+1 =

ϕ n+1 − ϕ n t

(14)

- Second order

ϕtn+1 =

Since the sign function is constant, the time derivative of Eq. (9) can be written as:

∂ (s(ϕ) − c · ∇ϕ) ϕt  t  = c · ∇(c · ∇ϕ) = (∇ · ∇)ϕ ∂t

−

or

ϕ

where x is the element characteristic length. The re-initialization equation can also be re-written as:

ϕt  t  =

n+1

t ∂ϕ n+1 t 2 ∂ 2 ϕ + 1! ∂ t 2! ∂ t 2 4 4 n+1 t ∂ ϕ + + ··· 4! ∂ t 4

ϕ n = ϕ n+1 −

is the scalar signed distance function, t is the pseudo-time,

where ϕ and s(ϕ ) is the sign function given by [6]:

α

ϕ n+1 − ϕ n t n+1 + ϕ t 2 tt

(15)

- Fourth order (10)

ϕtn+1 =

ϕ n+1 − ϕ n + α t ϕ˜tt t

(16)

3. Hybrid Taylor–Galerkin variational multi-scale formulation

where ϕ˜ is an intermediate value of the function ϕ given by:

Modeling high speed SLD impingement requires a scheme capable of handling accurately and in a stable fashion high velocities and severe topological changes at a reasonable computational cost. The VMS approach [21–23] assumes that the solution variable consists of two superimposed parts: a coarse scale solution which is computed on the available mesh and a fine scale part which is accounted for using analytical treatment. A TG [24] approach allows for a solution that is accurate in time. The approach proposed here attempts to combine the favorable characteristics from both schemes, thus creating a hybrid approach that is accurate in both space and time and effective for high-speed SLD dynamics. An implicit discretization is adopted.

ϕ˜ = βϕ n+1 + γ t ϕtn+1 + δ t 2 ϕttn+1

3.1. Weak Galerkin form of the LS equation Assuming a computational domain  discretized into nonoverlapping elements e with boundaries  e , the weak-Galerkin formulation of Eq. (1) can be written as:

(w, ϕt + uˆ · ∇ϕ)e = 0

(11)  where (·, ·)e = e ·d. Before applying the VMS approach to Eq. (11), the temporal derivative is expanded through an implicit Taylor series expansion. The order of the expansion will determine the temporal accuracy of the discretized equation. This is further explained in the following section.

with the constants α , β , γ and δ given in Table 1. The values of the constants can be obtained by substituting Eq. (17) into Eq. (16) and solving for the constant to match the coefficients of Eq. (13). It should also be noted that an iterative process is required since Eqs. (16) and (17) are dependent. The iterative process starts by assuming that ϕ n+1 is equal to ϕ n ; the intermediate value ϕ˜ is then obtained from Eq. (17). This intermediate value is then used in (16) to obtain the updated value of ϕ n+1 which is then used in the new iteration. 3.3. VMS approach Starting from the weak form of the LS equation, the temporal derivative is substituted with the Taylor series expansion presented in the previous section. Eqs. (14), (15), or (16) will replace the temporal term in the weak form of the LS equation presented in (11) resulting in a first, second, or fourth order expansion of the temporal derivative respectively. Only the coupling of the second order TG with the VMS approach will be explained in detail, a similar procedure can be followed to obtain the formulations for the higher orders. Starting from the weak Galerkin form of the LS equation and introducing the second order expansion of the temporal derivative, Eq. (11) can be written as:



3.2. Taylor series expansion

ϕ =ϕ n

n+1

+

k  i=1



∂ iϕ tn − t i! ∂ti

n+1 n w, ϕ −t ϕ



For the current scheme, an implicit formulation of the temporal derivative is proposed. The implicit Taylor series expansion can be written as:

 n+1 i





e

− w, 2t ut · ∇ϕ n+1



+ w, 2t uˆ · ∇ uˆ · ∇ϕ n+1







e

e

+ w, uˆ · ∇ϕ n+1



(18) e

=0

From Eq. (3):



w, (12)

(17)

ϕ − ϕn t ˆ ˆ + un · ∇ϕ + ∇( ∇ϕ)) =0 u · u · ( t 2  e

(19)

A. Bakkar et al. / Computers and Fluids 121 (2015) 192–205

where the (n + 1) superscripts is dropped for clarity. The VMS approach assumes that the scalar function ϕ and the weight function w can be decomposed into a sum of a coarse scale resolved by the grid, and a fine scale solved analytically.

ϕ(x) = ϕ( ¯ x) + ϕ  (x)  

 

Coarse scale

Fine scale

(20)

 

Coarse scale

Fine scale

ϕ  = w = 0 on e

(22)

Introducing the VMS approach to Eq. (19):

  ϕ¯ + ϕ  + un · ∇ ϕ¯ + ϕ  + w¯ + w , t



ϕn = w¯ + w , t





  t uˆ · ∇ uˆ · ∇ ϕ¯ + ϕ  2

e

×

    t  ϕ¯ + ϕ  n   ¯ + u · ∇ ϕ¯ + ϕ + uˆ · ∇ uˆ · ∇ ϕ¯ + ϕ w, t 2 e

n ϕ ¯ = w, t  e

    t  ϕ¯ + ϕ   + un · ∇ ϕ¯ + ϕ  + w, uˆ · ∇ uˆ · ∇ ϕ¯ + ϕ  t 2 e

ϕn  = w, t 

τ = b1







b2 we ,









+ b2 w e ,



t  2













uˆ · ∇ uˆ · ∇ b1 ϕ e





ϕ¯  t  

 t − b2 w e , ¯ (uˆ · ∇(uˆ · ∇ ϕ)) =



b2 we ,

ϕn t



2



−









− b2 we , un · ∇(ϕ) ¯

b2 w e ,





− w,

− w,



τ t 2

τ t 2

−1

b2 ,

t b1 + un · ∇ b1 + (uˆ · ∇(uˆ · ∇ b1 )) t 2





¯ un · ∇((uˆ · ∇(uˆ · ∇ ϕ)))



e

¯ uˆ · ∇(uˆ · ∇(un · ∇ ϕ))

τ t 2 4

e



¯ uˆ · ∇(uˆ · ∇(uˆ · ∇(uˆ · ∇ ϕ))) 

e

 τ  τ  ϕn n u · ∇ϕ n = w, 1 − − w, t t  t e e  τ − w, uˆ · ∇(uˆ · ∇ϕ n )



2

e

(33)

If linear shape functions are used, the higher-order derivatives in Eq. (33) can be dropped. The diffusion term that is introduced here is expected to improve the stability of the numerical scheme but may result in additional diffusion of the solution. For simplicity, the above formulation will be referred to as H2TV. Similarly, Eqs. (14) or (16) are used to substitute for the temporal derivative in Eq. (11). The hybrid VMS formulation for the first (H1TV) and fourth (H4TV) order Taylor series expansion of the temporal term are given by, respectively:



¯ 1− w,



τ  ϕ¯ t t







=

(28)



¯ 1− w,





¯ 1− w,

e

τ  ϕn t t

τ  ϕ¯ t t



¯ 1− + w,

¯ τ (uˆ · ∇(uˆ · ∇)ϕ)) −(w, ¯ e





−1  b1 n · ∇ b + t (u ˆ · ∇(uˆ · ∇ b1 ))  b2 ,  + u 1 t 2 

(31)

e





Taking the constant coefficients out of the integral, (28) can be rearranged as follows:

ϕe = 

b2 d 

−(w, τ un · ∇(un · ∇(ϕ))) ¯ e

(25)



(30) 

    2τ τ  ϕ¯ w, 1 − + w, 1 − un · ∇ ϕ¯ t t  t e e

  2τ t + w, 1− ¯ (uˆ · ∇(uˆ · ∇ ϕ)) 2 t 











(27)

+ b2 we , un · ∇ b1 ϕ e

t

ϕ ϕ¯ t + + un · ∇(ϕ) ¯ + ¯ (uˆ · ∇(uˆ · ∇ ϕ)) t t 2





− w,







n

Introducing the fine scale ϕ  presented in Eq. (31) into the coarse scale sub-problem and re-arranging, Eq. (24) can be re-written as:

(26)

on e

b1 ϕ e



(32)

(24)

where b1 , b2 , ϕ e and we are the bubble shape functions and fine scale coefficients for the trial solution and weight function, respectively. Substituting into Eq. (25):



· ∇ b1 + 2t (uˆ · ∇(uˆ · ∇ b1 ))

e

The fine scale problem, defined where (·, ·) = e=1 e element-wise, is then solved and a fine-scale solution of the scalar function is obtained. Introducing this into the coarse scale problem eliminates the fine scale component while still accounting for its effects. Assuming that the fine scale can be represented by “bubble functions” [22], the fine-scale variables can be written as:



b2 , −

+

un

where τ is the stabilization parameter defined as:

(·) de .

ϕ  = b1 ϕ e on e

(29) 

n

ϕ ϕ¯ t + + un · ∇(ϕ) ϕ  = −τ − ¯ + ¯ (uˆ · ∇(uˆ · ∇ ϕ)) t t 2

e

nelm 



or





−1 b1 b2 ,  t

(23) e

Since the fine scale vanishes on the element boundaries and by virtue of the linearity of the weight function, Eq. (23) can be split into a coarse-scale sub-problem (24) and fine-scale sub-problem (25).

w  = b2 w e

ϕn ϕ¯ t + + un · ∇(ϕ) ¯ + ¯ (uˆ · ∇(uˆ · ∇ ϕ)) t t 2

Substituting (29) into (26), the fine scale of the scalar function can be written as

(21)

where the fine scale is defined over the element ʹ but vanishes at the element boundary  ʹ .



b2 , −

ϕ  = b1 

w(x) = w¯ (x) + w (x)

 

×

195







¯ − w, e





e



 τ uˆ · ∇(ϕ n ) t e

¯ 1− + w,

¯ τ (uˆ · ∇(uˆ · ∇)ϕ)) −(w, ¯ e



2τ uˆ · ∇ ϕ¯ t e





2τ uˆ · ∇ ϕ¯ t e

(34)

196

A. Bakkar et al. / Computers and Fluids 121 (2015) 192–205



   τ  ϕn τ  ¯ 1− α t ϕ˜tt − w, t t  t e  e  τ ¯ ¯ ατ t uˆ · ∇ ϕ˜tt )e − w, + (w, uˆ · ∇(ϕ n ) t e =



where τ is the stabilization parameter defined as:

¯ 1− w,

(35)

where ϕ˜tt is defined as in Eq. (2). The τ stabilization parameter associated with H1TV and H4TV are, respectively:

τ = b1





b2 d  e

and

τ = b1

b1 + uˆ · ∇ b1 b2 , t





b2 d 

b2 ,

e

b1 + uˆ · ∇ b1 t

−1

(36) 

−1 

3.4. Re-initialization equation Following a procedure similar to the one presented above, the hybrid approach for the re-initialization can be obtained. To avoid redundancy, only the second order Taylor series expansion will be introduced. Starting from Eq. (8), the weak Galerkin form can be written as:

(w, ϕt  + c · ∇ϕ)e = (w, s)e

ϕ n+1 − ϕ n t n+1 + ϕ   + c · ∇ϕ n+1 t  2 tt



e

= (w, s)e

(39)

From (10), Eq. (39) can be written as:

w,

ϕ n+1 t  + (∇ · ∇)ϕ n+1 + c · ∇ϕ n+1 t  2



=

w, s +

e

ϕn t 

e

(40) Introducing the VMS approach presented earlier, the coarse (41) and fine (42) scale equations can be written as:



ϕ¯ + ϕ  t   ) + c · ∇(ϕ¯ + ϕ  ) + · ∇)( ϕ ¯ + ϕ (∇ t  2 e

n ϕ = w, s + t  

(41)

e

Substituting the fine scale presented in Eq. (44) into the coarse scale Eq. (41), the hybrid second order re-initialization equation can be written as:



    2τ τ  ϕ¯ ¯ 1− ¯ + 1 − c · ∇ ϕ ¯ w, w, t  t   t  e

e

   t ¯ ¯ τ t  c · ∇(∇ · ∇)ϕ¯ − 2τ (∇ · ∇)ϕ¯ + w, − w, e 

2





¯ τ − w,



e

2

2



(∇ · ∇)(∇ · ∇)ϕ¯ e

τ   τ  ϕn τ = w, 1 − s+ 1− − c · ∇ϕ n  t t  t  t   τ − (∇ · ∇)ϕ n 2

e

(46)

The fourth and fifth terms on the LHS can be dropped if linear shape functions are used. 4. Numerical results In this section, the proposed hybrid method is applied to a number of numerical benchmark cases. The first part of this section will compare both the spatial and temporal convergence rates of the method. This will help in better understanding the effect of the hybridization on the convergence accuracy. The stretching of a circular fluid element in a single vortex flow field was chosen as the most relevant to test the proposed formulation. The best scheme will then be considered for the subsequent section, where it is compared to other methods available in the literature. The test cases presented are the reversed stretching of a circular droplet in a single vortex and highspeed advection as proof of concept. Throughout this section, the error ε will be defined as:

ε=

|A(t ) − A(to)| A(to)

A(t ) =



(42)

Introducing the bubble functions presented in Eqs. (26) and (27), the fine scale equation can be re-written as:

ϕn ϕ¯ t  × b2 , −s − + + c · ∇ ϕ¯ + (∇ · ∇)ϕ¯   t t 2

(43) 

or



ϕn ϕ¯ t  ϕ  = −τ −s −  +  + c · ∇ ϕ¯ + (∇ · ∇)ϕ¯ t t 2

H (ϕ)d

(48)



H (ϕ) =

e

−1   b2 , b1t  + c · ∇ b1 + 2t (∇ · ∇)b1  

(47)

H(ϕ ) is a smoothed Heaviside function defined as:

  ϕ¯ + ϕ  t   + c · ∇(ϕ¯ + ϕ  ) + · ∇) ϕ ¯ + ϕ (∇ t  2 e

n ϕ = w , s + t  



t 



w ,

ϕ  = b1 

(45) 

where A(t) is the area (or volume in 3D) at time t, and A(to ) is the area at the initial time calculated using:

¯ w,



b2 d 

(38)

Replacing the temporal derivative with the second order implicit Taylor series expansion presented in (15), Eq. (38) can be re-written as:

w,

−1

b1 t  b2 , + c · ∇ b1 + (∇ · ∇)b  t 2

e

(37)

Detailed derivation of the H4TV approach can be found in Appendix A.



τ = b1





(44)

⎧ 0, ⎪ ⎪ ⎨  1

ϕ 1 1+ + sin ⎪ 2 ς π ⎪ ⎩ 1,

 π ϕ  ς

i f ϕ < −ς , i f |ϕ| ≤ ς if

(49)

ϕ>ς

where ς is the smoothing range, generally taken as 2x with x representing the element characteristic length. Observing Eq. (18), it is noticed that it has a resemblance to the advection-diffusion-reaction equation. Also, examining the stabilization parameter presented in (32), the presence of an advection component (u · ∇ ), a diffusion component (uˆ · ∇(uˆ · ∇)), and a reaction component (1/t) is evident. The same can be noticed for the stabilization parameters in Eqs. (36) and (37) except that they lack a diffusion component. For simplicity, the stabilization parameter τ presented in Eqs. (32), (36), and (37) will be replaced with a modified parameter similar to that used in stabilizing the advection-diffusionreaction equation. The modified stabilization used is given by [25,26]:

A. Bakkar et al. / Computers and Fluids 121 (2015) 192–205

τ=

1

τ12

+

1

τ22

+

by a stream function:

−1/2

1

(50)

τ32

h , 2 u 

ψ=

1

π

2

sin

(π x)sin2 (π y)

(52)

Thus, the velocity field components can be written as:

where

τ1 =

197

τ2 =

t 2

,

and τ3 =

h2 4υ

(51)

where h is the length scale defined in [27], and υ is the diffusion coefficient. The parameter τ 3 is only taken into account when computing the stabilization parameter for Eq. (32). A similar parameter can be deduced for the re-initialization equation. The hybrid methods were also compared with the explicit fourth order TG approaches, noted as (4TE), and results for the second order TG approach obtained by Cho et al. [15]. It should be noted that the implicit fourth order TG approach (4TI) failed due to instabilities arising from the presence of an advection term which was not stabilized, thus only the (4TE) was introduced. All simulations were performed using linear shape functions. 4.1. Convergence study To determine the rate of convergence of the proposed methods, a circular fluid element of 0.15 radius subjected to a single vortex flow field was chosen. The circular fluid element is initially located at (0.5,0.75) in a 1 × 1 square domain, with the flow field described

∂ψ 2 = 2sin (π x) sin (π y) cos (π y) ∂y ∂ψ 2 = −2 sin (π x) cos (π x)sin (π y) uy = − ∂x ux =

(53)

The total study time of 3 was used with all simulations performed on structured quadrilateral meshes with linear elements. For the spatial convergence study, all runs were performed using a time step of 0.005, while the temporal convergence study was performed on a 100 × 100 element structured quadrilateral mesh (Figs. 1–2). No reinitialization was performed to avoid any effects it might have on the performance of the methods. Table 2 shows that higher-order schemes, namely H2TV, H4TV and 4TE, have comparable conservation characteristics. Note that in Table 2, superscripts are used to denote whether the error indicates a gain or loss over the initial area as defined in (47). All methods resulted in a loss of area on the coarse meshes with the H1TV having the greatest loss of ∼80%. This loss was reduced as the mesh was refined. In some cases, refinement resulted in a gain of area except for H1TV. The higher-order hybrid methods and the fourth order 4TE had comparable accuracy throughout, with H4TV having the best overall performance, followed by H2TV, then 4TE. The effect of the time step

Fig. 1. Stretching of a circular fluid element of radius 0.15 in a single vortex on a quadrilateral mesh of (a) 50 × 50, (b) 100 × 100, (c) 150 × 150, and (d) 200 × 200 elements for a total time of 3 and a time step of 0.005.

198

A. Bakkar et al. / Computers and Fluids 121 (2015) 192–205

Fig. 2. Stretching of a circular element of radius 0.15 in a single vortex on a quadrilateral mesh of 100 × 100 elements for a total time of 3 with a time step of: (a) 0.01, (b) 0.005, (c) 0.001, and (d) 0.0005.

Table 2 Spatial convergence error % at time T = 3. Number of elements

50 × 50 100 × 100 150 × 150 200 × 200

Table 4 Spatial convergence absolute error % at time T = 3 using H2TV on structured and unstructured mesh.

Method H1TV

H2TV

H4TV

4TE

-79.660 -44.025 -33.188 -29.278

-38.556 -5.535 -0.382 +1.070

-38.547 -5.572 +0.310 +0.985

-38.674 -5.457 -0.778 +1.738

Table 3 Temporal convergence error % at time T = 3. Time step

0.0005 0.001 0.005 0.01

Method H1TV

H2TV

H4TV

4TE

-6.002 -9.302 -44.025 -88.549

-3.358 -3.395 -5.535 -9.113

-3.362 -3.405 -5.572 -9.091

-3.365 -3.411 -5.457 -9.052

Mesh type

Structured

Unstructured

Element Number of nodes Absolute error %

Quadrilateral 22801 0.382

Triangle 22814 1.696

reduced, the higher-order hybrid schemes remain slightly better at smaller time steps. It is thus clear that the higher-order schemes are superior to the H1TV scheme. It is also noticed that the hybrid approaches are better than the 4TE approach in terms of conservation. Given the fact that the H4TV approach requires an iterative process due to the equation splitting. The choice was made to utilize the H2TV approach in the rest of this paper. 4.2. Unstructured meshes

was also studied and reported in Table 3. Again, the H1TV approach had the greatest loss in area with the large time step, losing ∼89%. The higher order schemes were comparable with the hybrid approaches being slightly more dissipative at larger time steps. This is due to the effect of the stabilization parameter as it is partially based on the selected time step. While this difference diminishes as the time step is

To test the behavior of the H2TV scheme on an unstructured mesh, a circular fluid element of radius 0.15 initially located at (0.5,0.75) in a square domain of size 1 × 1 was subjected to single vortex flow such as that prescribed by Eq. (52). The unstructured and structured mesh had a comparable number of nodes and used the same time step of 0.005. As can be seen from Table 4, the conservation characteristics of

A. Bakkar et al. / Computers and Fluids 121 (2015) 192–205

199

Table 5 Change in area absolute error % for a single circle element. Number of elements

Time step 0.01

0.005

64 × 64 128 × 128 256 × 256 512 × 512

5.25 -

3.70 2.14 1.3 1.12

Table 6 Change in area absolute error % for two interesting circular elements. Number of elements

Fig. 3. Stretching of a circular fluid element of radius 0.15 in a single vortex for a total time of 3 and a time step of 0.002 on a 22,801-node structured mesh (dotted line) and a 22,814 nodes unstructured mesh (solid line).

the method on a structured quadrilateral mesh is slightly better than that on a triangular unstructured mesh. Upon inspecting Figs. 3–4, it is evident that the quadrilateral elements mesh was better at capturing the stretched tail region where the thickness of the interface was minimal. The leading region however was fairly similar on both meshes. 4.3. Re-initialization The proposed re-initialization formulation is tested by applying it to an initially distorted LS function. Two test cases will be presented, one of a single circular shape element and another of two intersecting circular shapes. The exact solution for these cases is well known and thus the final solution can be compared against the exact one. The error, representing the change in area enclosed by the interface is also provided. Starting with the single circular shape, the exact solution for a signed distance function is given by:

ϕ(x, y) =



x2

+

y2

−r



(54)

where r is the radius of the circle. This results in an infinity number of concentric signed distance functions with the zero LS representing the interface. This signed distance function is initially distorted such that the gradients vary widely by modifying Eq. (54) into:

  ϕ  (x, y) = (x − r)2 + (y − r)2 + 0.1 ϕ(x, y) where ϕ  is the distorted LS function.

(55)

64 × 64 128 × 128 256 × 256

Time step 0.01

0.005

2.77 -

0.76 0.58 0.53

The radius is set to one in a 4 × 4 domain with the solution time step and spatial discretization varied. Similarly, the exact solution for a signed distance function of two intersecting circular elements is given by:



⎧ √ 2 ⎪ 2 + (y ± 2 − a2 ) min x r ⎪ ⎪ ⎪ ⎨ x+a a−x a a ϕ(x, y) = and  if  ≥ ≥ 2 2 r r ⎪ 2 2 ⎪ (x + a) + y  ( a − x) + y  ⎪ ⎪ ⎩min (x ± a)2 + y2 − r Otherwise

(56) where (a) determines the intersection between the two circles and was set to 0.7. The contours are again distorted by modifying Eq. (60) into:

  ϕ  (x, y) = (x − r)2 + (y − r)2 + 0.1 ϕ(x, y)

(57)

where ϕ  is the distorted LS function. It should be noted that this is a slightly more difficult problem due to the presence of a kink at the intersection of the two circles. As it can be seen in Figs. 5 and 6, the initial solution is far from being a signed distance function. As the solution progresses in time, the contours deform and converge to form a signed distance function. Tables 5 and 6 list the total change in area at the final converged solution. It is evident that the interface has slightly drifted from its

Fig. 4. Tail (left) and leading edge (right) regions of a circular fluid element of radius 0.15 in a single vortex for a total time of 3 and a time step of 0.002 on a 22,801-node structured mesh (dotted line) and a 22,814 nodes unstructured mesh (solid line).

200

A. Bakkar et al. / Computers and Fluids 121 (2015) 192–205

Fig. 5. Level set values under re-initialization at (a) t = 0, (b) 200 iterations, (c) 400 iteration, (d) converged solution (dotted lines represent exact solution). Mesh size is 64 × 64, circle element radius of 1, and re-initialization time step of 0.005. Table 7 Change of area (%) for a circular fluid element of radius 0.15 in single vortex deformation in a (1×1) quadrilateral domain. Total time t = T = 8. Grid resolution

TG direct approach [29] TG [15] LSM [15] Hybrid particle LSM [29] Particle LSM no re-initialization [31] Particle LSM with re-initialization [31] H2TV no re-initialization H2TV with re-initialization

64 × 64 12.79 12.74 100 1.8 44.1 17.9 3.54 9.44

initial location, this is due to the fact that the interface does not coincide with the grid points. The information of the interface is interpolated from the element nodes using information from both sides of the interface [28]. This effect is reduced as the mesh is refined and a smaller re-initialization time step is used (Fig. 7). The sharp change in contours at the intersection point between the two circles is also noticed to slightly smear, this again is reduced and the corner tends to get sharper as the mesh is further refined (Fig. 8). 4.4. Stretching of a circular fluid element To demonstrate the capabilities of the H2TV method in modeling fluid stretching, the single vortex flow previously introduced consid-

128 × 128 6.84 12.33 39.8 0.7 2.6 4.2 2.36 2.66

256 × 256 3.55 6.37 10.3 0.4 0.8 2.1 0.305 2.15

512 × 512 1.99 3.48 0.028 0.071

ered again. The flow field however is reversed at half the study time to allow for comparison with literature. The modified velocity field can be written as:

 

ux = 2sin (π x) sin (π y) cos (π y) cos πτtt   2 uy = −2 sin (π x) cos (π x)sin (π y) cos πτtt 2

(58)

where t is the instantaneous time, and τ t is the total time. The fluid element studied is a circle of radius 0.15 with the center initially located at (0.50, 0.75) in a 1 × 1 square domain. The error at t = T = 8 associated with various mesh sizes as well as those of some of the methods available in literature are shown in Table 7 where results with and without re-initialization are presented. All simulations were run using a time step of 0.002 except for the 512 × 512

A. Bakkar et al. / Computers and Fluids 121 (2015) 192–205

201

Fig. 6. Level set values under re-initialization at (a) t = 0, (b) 200 iterations, (c) 400 iteration, (d) converged solution (dotted lines represent exact solution). Mesh size is 64 × 64 and re-initialization time step of 0.005.

Fig. 7. Re-initialized level set values on 64 × 64, 128 × 128, and 256 × 256 meshes, circle element radius of 1, and re-initialization time step of 0.005. Exact solution is dotted, 64 × 64: dashed, 128 × 128: dash dot, and 256 × 256: dash double dot.

element mesh with re-initialization where a time step of 0.001 was used. Re-initialization was applied every t = 0.5 using a pseudotime step of 0.0001. As the mesh was refined, it was noticed that the tail and leading edge regions at half time were better resolved. When returned to the original location, the shape of the element remained

Fig. 8. Re-initialized level set values on 64 × 64, 128 × 128, and 256 × 256 meshes, circle element radius of 1, and re-initialization time step of 5e − 3. Exact solution is dotted, 64 × 64: dashed, 128 × 128: dash dot, and 256 × 256: dash double dot.

fairly circular, even when using coarse meshes (Figs. 9 and 10). The proposed method succeeded in presersving the geometric features of the original shape despite being exposed to a severe stretching velocity field. Table 7 also shows that the proposed approach has improved conservation characteristics as compared to other approaches available in the literature. The exception was the hybrid particle LSM

202

A. Bakkar et al. / Computers and Fluids 121 (2015) 192–205

Fig. 9. Stretching of a circular fluid element of radius 0.15 in a single vortex on a quadrilateral mesh at (a) total time t = T = 8, and (b) half time t = 4. Time step was set to 0.002 and no re-initialization was used.

introduced by Enright et al. [29,30] which utilizes particles located on each side of the interface to correct its location as it progresses in time. This difference however diminishes as the mesh is refined. The current method, avoids the additional complications introduced by seeding and tracking particles around the interface. Performing reinitialization has led to additional losses in area. This is due to the additional iterations involved in the re-initialization process. The introduction of the re-initialization process also results in a scheme that is not time accurate as the order cannot be guaranteed. The method still provides improved accuracy when compared to other methods available in the literature. Stretching of a spherical fluid element in a three-dimensional vortex field was also conducted. The velocity field is given by: 2

ux = 2sin

(π x) sin (2π y) cos 2(π z) cos

uy = − sin (2π x)sin

2

(π y) sin (2π z) cos

uz = − sin (2π x) sin (2π y)sin

2

(π z) cos

 πt  τ  πtt  τ  πtt  τt

(59)

creating two rotating vortices. The spherical fluid element is of radius 0.15 and is initially located at (0.35,0.35,0.35) in a 1 × 1 × 1 cubic domain, with the total study time being 3. Figs. 11 and 12 show the stretching of the spherical element on 40 × 40 × 40 and 60 × 60 × 60 elements structured quadrilateral element meshes, respectively. It is clear that as the mesh is refined, the thin interface is better resolved.

The error at t = T = 3 associated with various mesh sizes, as well as those of some of the methods available in literature, are shown in Table 8.

4.5. Advection of a fluid element To evaluate the performance of the method for pure advection problems, a circular fluid element is transported under a constant velocity. As mentioned earlier, this approach was developed with the intent of modeling SLD impingement characterized by relatively large velocities. Consequently, advection velocities of 100 and 200 were used. The element was a circle of radius 0.15 initially located at (0.50, 0.25) in a 2 × 0.5 domain, and the total study time was 0.01 and 0.005 for the velocity of 100 and 200, respectively. The simulations were performed on a 200 × 50 structured quadrilateral meshes with a time step of 0.00001 (Fig. 13). Results indicate that the method is capable of preserving the shape of the circular element and the total change in area was negligible (Table 9). It is thus evident that the proposed scheme is capable of handling large advection velocities. Three-dimensional simulations were also performed, using a sphere of radius 0.15, initially located at (0.50,0.25,0.25) in a 2 × 0.5 domain. The total study time was 0.01 and 0.005 for advection velocities of 100 and 200, respectively (Figs. 14 and 15). The mesh used is 100 × 50 × 50 elements. The total change in volume can be found in Table 10.

A. Bakkar et al. / Computers and Fluids 121 (2015) 192–205

203

Fig. 10. Stretching of a circular fluid element of radius 0.15 in a single vortex on a quadrilateral mesh at (a) total time t = T = 8, (b) half time t = 4, (c) leading region at time t = 4, and (d) trailing region at time t = 4. Time step was set to 0.002 and the solution was re-initialization every t = 0.5 with a pseudo time step of 0.0001.

Fig. 11. Spherical fluid element in a single vortex flow on a structured 40 × 40 × 40 mesh, total time is 3.

Fig. 12. Spherical fluid element in a single vortex flow on a structured 60 × 60 × 60, total time is 3s.

Fig. 13. Advection of a circular fluid element in a constant velocity field on a (200 × 50) quadrilateral mesh.

Table 8 Change of area (%) for a spherical fluid element of radius 0.15 in a (1×1×1) quadrilateral domain. Grid resolution 40 × 40 × 40 Time (total = 3) LSM [15] Hybrid particle LSM [29] H2TV

1.5 26.14

60 × 60 × 60 3 9.26

1.5 10.75

100 × 100 × 100 3 0.814

1.5 49 1.9 -

3 80 2.6 -

204

A. Bakkar et al. / Computers and Fluids 121 (2015) 192–205 Table 10 Change of area (%) for a spherical fluid element in an advection velocity field at total time in 3D.

Table 9 Change of area (%) for a circular fluid element in an advection velocity field at total time in 2D. Velocity (units/s)

Grid resolution

100 200

0.0038 0.0115

Velocity (units/s)

Grid resolution

100 200

0.00017 0.0478

200 × 50

Fig. 14. Final location of an advected sphere in a constant velocity field of 100, total time 0.01.

100 × 50 × 50

was decided to utilize the H2TV. The difference in terms of accuracy versus computational robustness also favored the H2TV. Since the stabilization parameter chosen depends on the time step, the effect of the chosen time step was studied. The hybrid approaches were found to be slightly dissipative at a large time step. This effect diminished by reducing the time step. To test the capabilities of the chosen H2TV scheme, benchmark cases introducing stretching and high-speed advection were used. When subject to a single vortex flow, the scheme demonstrated good conservation characteristics compared to other methods available in the literature. The accuracy was comparable to the hybrid particle LS method introduced in the literature without the additional complications of tracking marker particles. The approach was also capable of preserving the geometric characteristics of the original circle, even on a coarse mesh. The method was also capable of handling high speed advection with ease. The next step of the current research is to couple the method with a flow solver in order to address multi-phase droplet dynamics by considering the effect of shear stresses acting on the surface of the droplet. Such individual droplet simulations could then be used to gradually refine our understanding of the SLD impingement at high speeds. This will help remove some of the current empiricism adopted in the Eulerian, and also Lagrangian, approaches to droplet dynamics simulations in in-flight icing. Appendix A Hybrid fourth order TG VMS (H4TV) Starting from the weak-Galerkin form of the LS equation and introducing the fourth order expansion of the time derivative, Eq. (11) can be written as:



w,

ϕ n+1 + u · ∇ϕ t





= e

w,

ϕn t −α ϕ˜ t 2 tt



(60) e

The VMS approach assumes that the scalar function ϕ and the weight w function can be decomposed into an overlapping sum of a coarse scale, resolved by the grid, and a fine scale solved for analytically.

Fig. 15. Final location of an advected sphere in a constant velocity field of 200, total time 0.005.

5. Conclusions A stabilized FEM implementation for the LS equation has been proposed, showing the potential of accurate modeling of droplet dynamics in a scenario typical of in-flight icing involving high speeds, large shear deformation due to aerodynamic forces and topological changes due to break-up of the droplets. The approach combines the VMS and the TG methods, providing stabilization for the LS equation. A convergence study was conducted to evaluate various order expansions of the Taylor series. The different orders of the Taylor series expansion were also compared with a pure VMS approach and an explicit fourth order TG approach. Due to the iterative process involved with the H4TV approach, which impacted the computational cost, it

ϕ(x) = ϕ( ¯ x) + ϕ  ( x)     Coarse scale

(61)

Fine scale

w(x) = w¯ (x) + w (x)

 

 

Coarse scale

Fine scale

(62)

where fine scale is defined over the element interiors ʹ but vanishes on the element boundary  ʹ.

ϕ  = w = 0 on e

(63)

Introducing the VMS approach to Eq. (60):





  (ϕ¯ + ϕ  ) + u · ∇ ϕ¯ + ϕ  t

ϕn t −α = w¯ + w , ϕ˜tt t 2 

w¯ + w ,

e

e

(64)

A. Bakkar et al. / Computers and Fluids 121 (2015) 192–205

Since it is assumed that the fine scale vanishes on the element boundaries, and utilizing the linearity of the weight function, Eq. (23) can be split into a coarse-scale (65) and fine-scale (66) sub-problems.







  (ϕ¯ + ϕ  ) + u · ∇ ϕ¯ + ϕ  t

¯ w,

w ,

  (ϕ¯ + ϕ  ) + u · ∇ ϕ¯ + ϕ  t

where (·, ·) =

n elm e=1

= e

¯ w,





ϕn − α t ϕ˜tt t

w ,

= 



ϕn − α t ϕ˜tt t

(65) e



(66) 



e e (·)d . Assuming that the fine scale can be

represented by “bubble functions” [22], the fine-scale variables can be written as: 

ϕ  = b1 ϕ e on e  w  = b2 w e

(67)

on e

(68) 



where b1 , b2 , ϕ e and we are the bubble shape functions and fine scale coefficients for the trial solution and weight function, respectively. Substituting in Eq. (66):







b2 we ,





b2 we ,







+ u · ∇ b1 ϕ e

t

=



b1 ϕ e





ϕ t ϕ¯ ϕ˜ − −α − u · ∇ ϕ¯ t 2 tt t n

(69) 

Taking the constants coefficients out of the integral, (69) can be rearranged as:

−1



ϕe = 



b b2 , (1t) + u · ∇(b1 )





ϕn ϕ¯ + α t ϕ˜tt + + u · ∇ ϕ¯ × b2 , − t t

(70) 

Substituting (70) into (67), the fine scale of the scalar function can be written as:

ϕ  = b1 

−1

× b2 , − or



ϕ  = −τ



b b2 , (1t) + u · ∇(b1 )



ϕ ϕ¯ + α t ϕ˜tt + + u · ∇ ϕ¯ t t n

ϕn ϕ¯ + α t ϕ˜tt + + u · ∇ ϕ¯ − t t

(71) 

(72)

where τ is the stabilization parameter defined as:

τ = b1





b2 d 

b2 ,

e

b1 + u · ∇ b1 t

−1

(73) 

Introducing the fine scale ϕ  presented in Eq. (72) into the coarse scale sub-problem and re-arranging, Eq. (66) can be re-written as:





τ  ϕ¯ ¯ 1− w, t t





e

¯ τ (u · ∇(u · ∇)ϕ)) −(w, ¯ e



=



¯ 1− w,



¯ 1− + w,

τ  ϕn t t









¯ 1− − w, 



2τ u · ∇ ϕ¯ t e

 τ  α t ϕ˜tt t e

 e τ n ¯ ¯ ατ tu · ∇ ϕ˜tt )e u · ∇(ϕ ) − w, + (w, t e 

(74)

where ϕ˜tt is defined as in Eq. (2). If linear shape functions are used, the higher order derivatives can be dropped. The intermediate variable ϕ˜ is obtained by applying the weak Galerkin formulation to Eq. (17) where no additional stabilization is required.

205

References [1] F.A.A. Airplane and engine certification requirements in supercooled large drop, mixed phase, and ice crystal icing conditions; 2014. [2] Papadakis M, Rachman A, Wong s-C, Yeong H-W, Hung KE, Vu GT, Bidwell CS. Experimental study of supercooled large drop impingement on ice shapes. DOT/FAA/AR-08/13; 2008. [3] Bilodeau DR, Habashi WG, Baruzzi G, Fossati M. An Eulerian re-impingement model of splashing and bouncing supercooled large droplets. In: 5th AIAA atmospheric and space environments conference, Paper Number: AIAA 2013-3058; 2013. [4] Honsek R, Habashi WG, Aubé MS. Eulerian modeling of in-flight icing due to supercooled large droplets. J Aircraft 2008;45(4):1290–6. [5] Uzgoren E, Singh R, Sim J, Shyy W. Computational modeling for multiphase flows with spacecraft application. Prog Aerospace Sci 2007;43(4–6):138–92. [6] Osher S, Fedkiw RP. Level set methods and dynamic implicit surfaces. Applied Mathematical Science, Vol. 153. Springer; 2003. [7] Sethian JA, Smereka P. Level set methods for fluid interfaces. Ann Rev Fluid Mech 2003;35:341–72. [8] Hyman JM. Numerical methods for tracking interfaces. Phys D: Nonlinear Phen 1984;12(1–3):396–407. [9] Scardovelli R, Zaleski S. Direct numerical simulation of free-surface and interfacial flow. Ann Rev Fluid Mech 1999;31(1):567–603. [10] Osher S, Sethian JA. Fronts propagating with curvature dependent speed: Algorithms based on Hamilton-Jacobi formulation. J Comput Phys 1988;79:12–49. [11] Sethian JA. Level set methods and fast marching methods: evolving interfaces in computational geometry, fluid mechanics, computer vision, and materials science. Cambridge University Press; 1999. [12] Codina R. Comparison of some finite element methods for solving the diffusionconvection-reaction equation. Comput Method Appl Mech Eng 1998;156(1– 4):185–210. [13] Nagrath S, Jansen KE, RTL Jr. Computation of incompressible bubble dynamics with a stabilized finite element level set method. Comput Method Appl Mech Eng 2005;194(42–44):4565–87. [14] Valance S, Borst R de, Réthor J, Coret M. A partition-of-unity-based finite element method for level sets. Int J Numer Method Eng 2008;76(10):1513–27. [15] Cho MH, Choi HG, Yoo JY. A direct reinitialization approach of level-set/splitting finite element method for simulating incompressible two-phase flows. Int J Numer Method Fluid 2011;67(11):1637–54. [16] Cho MH, Choi HG, Choi SH, Yoo JY. A Q2Q1 finite element/level-set method for simulating two-phase flows with surface tension. Int J Numer Method Fluid 2012;70(4):468–92. [17] Farthing MW, Kees CE. Evaluating finite element methods for the level set equation; Technical Report, Coastal and Hydraulics Laboratory - Engineer Research and Development Center; 2009. [18] Kees CE, Akkerman I, Farthing MW, Bazilevs Y. A conservative level set method suitable for variable-order approximations and unstructured meshes. J Comput Phys 2011;230(12):4536–58. [19] Min C. On reinitializing level set functions. J Comput Phys 2010;229(8):2764–72. [20] Sussman M, Smereka P, Osher S. A level set approach for computing solutions to incompressible two-phase flow. J Comput Phys 1994;114(1):146–59. [21] Hughes TJR. Multiscale phenomena: Green’s functions, the Dirichlet-to-Neumann formulation, subgrid scale models, bubbles and the origins of stabilized methods. Comput Method Appl Mech Eng 1995;127(1–4):387–401. [22] Masud A, Khurram RA. A multiscale/stabilized finite element method for the advection–diffusion equation. Comput Method Appl Mech Eng 28 May, 2004;193(21–22):1997–2018. [23] Khurram RA, Zhang Y, Habashi WG. Multiscale finite element method applied to the Spalart–Allmaras turbulence model for 3D detached-eddy simulation. Comput Method Appl Mech Eng 2012;233–236:180–93. [24] Donea J. A Taylor–Galerkin method for convective transport problems. Int J Numer Method Eng 1984;20(1):101–19. [25] Tezduyar TE, Osawa Y. Finite element stabilization parameters computed from element matrices and vectors. Comput Method Appl Mech Eng 2000;190(3–4):411– 30. [26] Tezduyar T. Calculation of the stabilization parameters in finite element formulations of flow problems. Applications of computational mechanics in structures and fluids, Cimne, Barcelona; 2005. p. 1–19. [27] Tezduyar TE, Park YJ. Discontinuity-capturing finite element formulations for nonlinear convection-diffusion-reaction equations. Comput Method Appl Mech Eng 1986;59(3):307–25. [28] Russoa G, Smerekab P. A remark on computing distance functions. J Comput Phys 2000;163(1):51–67. [29] Enright D, Fedkiw R, Ferziger J, Mitchell I. A hybrid particle level set method for improved interface capturing. J Comput Phys 2002;183(1):83–116. [30] Yu C-C, Heinrich JC. Petrov-Galerkin method for multidimensional, timedependent, convective-diffusion equations. Int J Numer Method Eng 1987;24(11):2201–15. [31] Hieber SE, Koumoutsakos P. A Lagrangian particle level set method. J Comput Phys 2005;210(1):342–67. [32] Beaugendre H. A PDE-Based 3D Approach to In-Flight Ice Accretion. Mechanical Engineering Department. Montreal, McGill University; PhD; 2003. p. 336. [33] Beaugendre H, Morency F, Habashi WG. Development of a second generation inflight icing simulation code. J Fluids Eng 2006;128(2):378–87.