**2. Governing equations**

fixed in the Eulerian method and fluid regions change in shape and location on the mesh. It uses a fixed grid system which is not transformed during the solution procedure. The fluid is studied while passing this volume and continuously replaced in time. Therefore, this method is not appropriate for formulation of basic equations of fluid movement. The Eulerian method has some limitations. For example, when the portion of the perimeter to the area of a zone of fluid is large, the error of this method is increased. In the Eulerian method, it is not possible to decompose the equation on the boundaries with the same precision of inner region of fluid and accordingly, the finer mesh should be used near the boundaries. Therefore, when the free surface of a discontinuous region is modeled by this method, finer grid should be employed in order to achieve more precise results, specifically if this surface has large deflections. This is crucial when the portion of the area to the perimeter of a zone is low, for example on phase of a multiphase fluid. In this case, using finer mesh could increase the portion of the number of the inner elements to the boundary elements, which in turn, increases the precision of the numerical solution. The main superiority of the Eulerian description is the possibility of modeling of complicated surfaces. For example, the collapse of a column of a fluid could be modeled in the Eulerian grid

which is shown in **Figure 1**.

366 Numerical Simulation - From Brain Imaging to Turbulent Flows

**Figure 1.** Fluid column in Eulerian grid: (a) before collapse and (b) collapsing flow.

**Figure 2.** A sample of Lagrangian grid in vertical direction.

Governing equations for a compressible viscous fluid flow with no phase change are as follows:

$$\frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf{u}) = 0 \tag{1}$$

$$\frac{\partial \left(\rho \mathbf{u}\right)}{\partial t} + \nabla \cdot \left(\rho \mathbf{u} \mathbf{u}\right) + \nabla P = \rho \mathbf{g} + \mathbf{F}\_s + \nabla \cdot \left[\mu \left(\nabla \mathbf{u} + \nabla \mathbf{u}^\top\right)\right] \tag{2}$$

In these equations, *ρ*, **u**, *t*, *P*, *μ*, *g* =(*gx*, *gy*), and **F**s are density, velocity vector, time, total pressure, kinematic viscosity, gravity acceleration, and body forces, respectively. Body forces include forces due to surface tension in the interface. Here, properties of a fluid such as density and viscosity are included in the equations. However, it should be kept in mind that the information changes from one fluid to another. Thus for mesh‐based numerical methods, new properties based on fluid properties of both materials should be considered for Eq. (2) in the cell containing the free surface, and the governing equations should be rewritten in the following form:

$$\frac{\partial}{\partial t}(a\rho\_\mathbb{I}) + \nabla \cdot (a\rho\_\mathbb{I}\mathbf{u}\_\mathbb{I}) = 0 \tag{3}$$

$$\frac{\partial}{\partial t}((1-\alpha)\,\rho\_2) + \nabla \cdot \left((1-\alpha)\,\rho\_2 \mathbf{u}\_2\right) = 0\tag{4}$$

$$\frac{\partial \left(a\rho\_{\rm l}\mathbf{u}\_{\rm l}\right)}{\partial t} + \nabla \cdot \left(a\rho\_{\rm l}\mathbf{u}\_{\rm l}\mathbf{u}\_{\rm l}\right) + \nabla P\_{\rm l} = a\rho\_{\rm l}\mathbf{g} + \mathbf{F}\_{\rm s} + \nabla \cdot \left[a\mu\_{\rm l}\left(\nabla \mathbf{u}\_{\rm l} + \nabla \mathbf{u}\_{\rm l}^{\rm r}\right)\right] \tag{5}$$

$$\begin{aligned} \frac{\partial \left( (1 - \alpha) \, \rho\_2 \mathbf{u}\_2 \right)}{\partial t} + \nabla \cdot \left( (1 - \alpha) \, \rho\_2 \mathbf{u}\_2 \mathbf{u}\_2 \right) + \nabla P\_2 &= \left( 1 - \alpha \right) \rho\_2 \mathbf{g} + \mathbf{F}\_s + \\ \nabla \cdot \left[ \left( 1 - \alpha \right) \mu\_2 \left( \nabla \mathbf{u}\_2 + \nabla \mathbf{u}\_2 \right)^\top \right] &\end{aligned} \tag{6}$$

where indices 1 and 2 show first and second fluids properties, and *α* is a scalar phase indica‐ tor function which is defined as follows:

$$\alpha = \begin{cases} 1 & \text{Control volume is filled only with phase } l \\ 0 & \text{Control volume is filled only with phase } 2 \\ 0 < \alpha < l & \text{Interface present} \end{cases} \tag{7}$$

This phase indicator function is the fluid property or volume fraction, which moves with it and can be derived as follows:

$$\frac{\partial \alpha}{\partial t} + \frac{\partial U\_i \alpha}{\partial \mathbf{x}\_i} = \mathbf{0} \tag{8}$$

This function can be used to calculate the fluid properties in each phase as a weight function. In order to use a set of governing equations using the weight function, each fluid property should be calculated based on the volume occupied by this fluid in the surface cell as ex‐ pressed in Eqs. (9) and (10) [2]:

$$
\rho = a \cdot \rho\_1 + (1 - a) \cdot \rho\_2 \tag{9}
$$

$$
\mu = \alpha \cdot \mu\_1 + \left(1 - \alpha\right) \cdot \mu\_2 \tag{10}
$$

Free surfaces considered here are those on which discontinuities exist in one or more varia‐ bles. This has been the challenge for researchers to omit or reduce this problem as much as possible. The transient state as well as phenomena such as surface tension, changing of fluid phase and Kelvin‐Helmholtz instability makes numerical simulation of such problems cumbersome. It is expected that methods used to simulate interface of fluids have a number of characteristics. These include mass conservation, simulating the interface as thin as possible, being able to reproduce complicated topologies, generalization of expansion to 3D problems, and being able to model surface phenomena and be computationally efficient.
