Problem Description

My model includes Two-Phase Flow, Level Set or the Two-Phase Flow, Phase Field interface and it is either not converging, taking a very long time to solve, or not converging to a reasonable result. This article contains guidance for solving such models.

Solution

Background

Fluid flow is governed by the Navier-Stokes equations, which solve for fluid velocity and pressure fields everywhere within the modeling domain. The fluid properties (viscosity and density) are either constant, or vary relatively smoothly in space as a function of temperature, pressure, shear rate, etc... If, however, we want to model two immiscible fluids then the fluid properties will vary significantly across the interface between the two fluids, and the surface tension effects can become important, as well as the effect of contact angles at wetted walls.

To model this in COMSOL Multiphysics, we can use either the Level Set or the Phase Field methods. Both of these methods introduce an additional scalar field (the level set or phase field function) to the modeling domain. These scalar fields vary smoothly between -1 and +1 everywhere, and are used to define the the fluid viscosity and density in the governing Navier-Stokes equations. The transition (where the fields vary from 0 to 1 for the level set implementation and -1 to +1 for phase field) is quite abrupt in space, thus giving a good resolution of the two phases. The governing equations for the Level Set and Phase Field methods are a type of Convection-Diffusion equation, with the advective term coming from the Navier-Stokes equations. Such equations are, however, quite numerically challenging to solve due to the combination of the significant advective term, abrupt transition in the field, and strong coupling to the Navier-Stokes equations.

The level set function varies smoothly and defines the interface between the fluids. The boundary is advected by the fluid flow (back arrows.)

Guidelines for Two-Phase Flow Modeling

If you want to model two-phase flow that does not involve any topology changes, that is, no droplet breakup, free surface formation, etc... and if the free surface does not undergo significant shape changes, you may be able to use the Laminar Two Phase Flow, Moving Mesh interface available with the Microfluidics Module only. This formulation, if it can be used, will require less computational resources.

Geometry

If at all possible, begin your modeling in 2D or 2D axisymmetry. Such models will be significantly faster to solve than 3D models, and serve as a useful test for trying out different model settings that will be applicable for 3D models.

If reasonable, avoid sharp edges and corners at any boundaries where you expect the fluid interface to traverse. Introduce a fillet in the geometry so that there is a smooth transition between boundaries.

Mesh

The mesh should have roughly equal mesh size throughout the modeling and the mesh elements should be isotropic, that is, the elements should not be stretched or compressed in one direction. The mesh size must be small enough to give a good resolution of the fluid interface. A good rule of thumb is to begin with a global element size that is one-tenth of the smallest expected droplet size, and to refine the mesh from there. Use Global Definitions > Parameters to define a parameter h_max that is used within the Mesh > Size setting to define the Maximum element size in all domains. Once you have arrived at a converged solution on a particular mesh, repeat the analysis with smaller mesh size. Once the solution does not change significantly with mesh refinement, the solution can be converged with respect to the mesh.

If it is known that one sub-domain of the entire modeling domain will only ever contain the same fluid phase, then a coarser mesh can be used in that sub-domain.

There is not a significant difference between structured and unstructured meshes.

Physics Interfaces

All fluid inlets should be in regions where there is a single phase entering the modeling domain, and the phase entering at that inlet should match the initial phase inside the adjacent domain. You may need to alter the geometry slightly to achieve this.

For both the Level Set and Phase Field approaches, there is a Two-Phase Flow, Level Set/Phase Field feature within the Multiphysics branch. Make sure to choose the appropriate material from the drop-down list for Fluid 1 and Fluid 2 in the settings.

Within the Laminar Flow physics interface, make sure to use the incompressible formulation. Make sure that all specified velocities or pressures at inlets ramp smoothly up from consistent initial conditions, as described here: Solving time dependent CFD simulations In case of problems involving gravity, make sure to enable the Include gravity option under the settings for Laminar Flow instead of using a Volume Force condition so that the hydrostatic pressure is automatically taken into account under Initial Values and under other pressure conditions.

If modeling bubbles with surface tension, make sure to define the additional initial pressure inside of the bubble that balances the surface tension. A reasonable value for this initial pressure can be computed from the Young-Laplace equation.

When using either the Level Set or Phase Field interfaces, there is a Level Set/Phase Field Model feature which is used to define the Parameter controlling interface thickness. If this value is too small, it can lead to numerical instabilities and if it is too large, the interface movement is not captured correctly. A good value is half the maximum element size: h_max/2 where h_max is the parameter controlling mesh size described earlier.

When using the Level Set physics interface, the Level Set Model feature also defines the Reinitialization parameter. If the reinitialization parameter is too small, it can lead to numerical instabilities, on the other hand if it is too big the interface movement is not captured correctly. A good estimate for the Reinitialization parameter is the maximum expected velocity magnitude.

When using the Phase Field physics interface, the Phase Field Model feature also defines the Mobility tuning parameter. Again, too small a value leads to numerical instabilities, too large a value will not capture the interface movement correctly. A good initial estimate for the Mobility tuning parameter is:

2*u_max/(3*sqrt(2)*sigma)*(h_max/epsilon)

where u_max is the expected maximum velocity magnitude, sigma is the surface tension coefficient, h_max is the value of the parameter controlling maximum element size, and epsilon is the value of the Parameter controlling interface thickness, which is typically h_max/2, in which case the above expression can be reduced to:

(4/3/sqrt(2))*(u_max/sigma)

When using the Level Set approach, there is a Multiphysics > Wetted Wall feature, and when using the Phase Field approach, there is a Phase Field > Wetted Wall feature that specifies a Contact Angle. This contact angle should be in agreement with the actual initial angle of the two phases, as defined by the geometry and the Phase Field/Level Set > Initial Interface feature. If the contact angle as defined by the initial interface is different, then the specified Contact angle should be ramped up, in a way similar to the velocity ramping ( see: Solving time dependent CFD simulations )

Within the interface settings for the Laminar Flow, the Level Set and the Phase Field interface, there is a Discretization section where you can switch the element order. You will need to enable Discretization within the Show menu in the Model Builder to see this setting. The default discretization is P1+P1 for the Laminar Flow, and Linear for the Level Set and Phase Field. Higher order discretization will be significantly more computationally intensive and is not usually recommended, instead, refine the mesh.

Study Settings

The study always consists of two steps, first, a Phase initialization step initializes the level set or phase field variable such that it is smoothly varying everywhere. Next, a Time Dependent step simultaneously solves the Navier-Stokes equations and the Level Set or Phase Field equation.

Within the Time Dependent study step settings, there is a Tolerance: field that defaults to Physics controlled. This can be changed to User controlled and tighter solver relative tolerances should be investigated to confirm that the solution is converged with respect to relative tolerance. That is, once a solution is computed, repeat the solution with a finer solver tolerance. Once the solution does not change significantly with solver tolerance as well as mesh refinement, the solution can be considered converged. Do not study finer solver tolerances before addressing all of the above points.

Advanced Physics Settings

If, as you perform a mesh refinement study and a solver tolerance refinement study, you observe an issue with mass conservation (total mass entering and leaving the domain do not match, or the total mass in closed domain changes) then switch to the more computationally intensive conservative form of the governing equations. Within the interface settings for both the Level Set and Phase Field, there is an Advanced Settings section where you can switch between the Nonconservative form and Conservative form of the Convection term. You will need to enable Advanced Physics Options within the Showmenu in the Model Builder to see these settings.

Disclaimer

COMSOL makes every reasonable effort to verify the information you view on this page. Resources and documents are provided for your information only, and COMSOL makes no explicit or implied claims to their validity. COMSOL does not assume any legal liability for the accuracy of the data disclosed. Any trademarks referenced in this document are the property of their respective owners. Consult your product manuals for complete trademark details.