Two Methods for Modeling Free Surfaces in COMSOL Multiphysics®

May 15, 2018

There are four methods for modeling free liquid surfaces in the COMSOL Multiphysics® software: level set, phase field, moving mesh, and stationary free surface. In the first part of this blog series, we discuss the level set and phase field methods, which are field-based methods that describe almost any type of free liquid surface. In part two, we will compare the results from this post with those obtained using the Moving Mesh interface for solving free surface problems.

What Are the Level Set and Phase Field Methods?

The level set and phase field methods are both field-based methods in which a free fluid surface is represented as an isosurface of the level set or phase field functions. The free liquid surface corresponds to the phase boundary between the liquid and the gas and is represented on a fixed mesh.

The figure below shows the surface of two droplets in a channel, taken from the droplet breakup model in the Application Library of the add-on Microfluidics Module. We can see in this plot that although the surface of the droplet is very sharp, the elements around the droplet do not adhere to the droplet’s surface.

A droplet breakup model illustrating that the droplet surface does not coincide with the element's surface.
The droplet surface does not coincide with the surface of the elements, neither in the level set method nor in the phase field method.

The level set and phase field functions are both advected by the velocity vector computed by the Navier-Stokes equations. For both the level set and the phase field methods, this looks like:

(1)

\[\frac{{\partial \phi }}{{\partial t}} + \nabla \phi \cdot u = F\]

Note that we are using Φ for both the level set and phase field functions. The right-hand side of the equation, F, is where the two methods differ. In the original level set method, F = 0, which gives a pure advective transport equation. However, the numerical solution with F = 0 is unstable and of small practical use in most cases. Instead, terms with higher-order derivatives of Φ, designed to keep the interface compact, are introduced in F in the level set method.

In the phase field method, F represents a term that tries to minimize the free energy of the system. Also, higher-order derivatives of Φ are introduced in this term. In fact, the source term in the phase field equation includes fourth-order terms. This means that, for practical reasons, the equation is often broken up into two equations, where an auxiliary dependent variable is defined as a function of the second derivatives of Φ. This is also what is done in COMSOL Multiphysics.

Both methods introduce the source term into the Navier-Stokes equations that originate from the surface tension at the free liquid surface. In the level set method, the curvature of the level set isosurface that represents the free boundary is used to describe surface tension. In the phase field method, the contribution to the Navier-Stokes equations from the surface tension is calculated from the chemical potential, which yields the surface tension and the gradient of the phase field function close to the interface.

The properties of the fluid transition smoothly from liquid to gas over a given value of the level set or phase field functions; i.e., the value that represents the free surface.

The level set function, Φ, varies between 0 and 1 across the free surface and is constant at 0 or 1 in the bulk of the two fluids. The level set function can, for example, be 0 for the liquid and 1 for the gas phase. The free liquid interface; i.e., the phase boundary between liquid and gas; corresponds to the value of the level set function Φ = 0.5. The density is thus a function of the level set function according to:

(2)

\[\rho = {\rho _1} + \phi \left( {{\rho _2}-{\rho _1}} \right)\]

As well as the dynamic viscosity:

(3)

\[\mu = {\mu _1} + \phi \left( {{\mu _2}-{\mu _1}} \right)\]

We can see in the two equations above that when Φ = 0, we get the properties of fluid 1 (for example, the liquid) and when Φ = 1, we obtain the properties of fluid 2 (for example, the gas). We obtain a smooth transition of fluid properties across the free surface; i.e., as smooth as the level set function.

The phase field function, Φ, varies between -1 and 1, where the free liquid surface is the isosurface for Φ = 0. The computation of the viscosity and density is done in an analogous way to the level set method, but with a different expression, since the phase field function varies between -1 and 1.

Settings in the Level Set and Phase Field Interfaces

The Level Set and Phase Field interfaces contain a number of parameters that need tuning in order for the methods to perform optimally. The values of these parameters depend on the modeled system and the numerical discretization of the equations. The size of the geometry (length, width, and depth); properties of the fluids; wetting of the walls; initial conditions; operating conditions; and size of the elements in the created mesh all have an impact on the selection of the model settings.

The Level Set Interface

Let’s start with the level set method’s settings. The reinitialization parameter, γ in the user interface, makes sure that the gradients in the level set function are concentrated to the free surface over time; i.e., to the interface between water and air in our example below. It does not set the thickness of the surface, but it makes sure that the variations of the level set function are contained within this thickness. If we set this value too low, variations in the level set function can be entrapped in the bulk of one of the fluids, creating liquid-gas interfaces out of nowhere. A value that is too high tends to give small time steps and large computation times.

A screenshot of the Level Set interface settings in the COMSOL software.
The Settings window for the Level Set interface.

The interface thickness parameter, denoted ε in the user interface, is more intuitive to understand. It simply controls the thickness of the region where the variations of the level set function are obtained. A large thickness smears out the shape of the free surface. The problem becomes easier to solve, but accuracy in the shape of the free surface is lost. A small value of ε yields a compact representation of the surface but requires a corresponding resolution of the mesh. There is no point in setting a value that is much smaller than the mesh element size, since the interface then becomes ragged and unresolved surface tension forces might cause spurious velocities. The value of ε should be set approximately to the same value of the largest elements close to the surface.

The Phase Field Interface

If we now look at the phase field method, the interface thickness parameter, ε, is almost identical to the level set method above and the same reasoning can be applied to the selection of its value.

A screenshot of the Phase Field interface settings in the COMSOL software.
The settings for the Phase Field interface.

The mobility tuning parameter for the phase field method, χ in the user interface, determines in a certain way the diffusivity of the phase field. This value has to be large enough for the phase field equation to be stable but small enough to give a sharp interface. The proper value is proportional to the speed at the surface and inversely proportional to the coefficient of surface tension.

(4)

\[\chi \propto \frac{{\left| {\mathbf{u}} \right|}}{{\sigma}}\]

This means that once we have found a proper value of χ for a certain set of operating conditions, we can use the relation above to set χ for a new set of conditions.

Modeling a Free Liquid Surface

The figures below show the example problem that we are going to use to investigate the level set and phase field methods for modeling free surfaces in COMSOL Multiphysics.

A solid bar is half immersed in water in a small channel. The bar is put into motion back and forth, tangential to the water surface, in order to generate a gentle wave on the water surface. The dimensions of the channel and bar are small enough to keep the flow laminar.

The model geometry of a rectangular bar moving in a liquid with a free surface.
Geometry and definition of the example problem.

Note that the moving mesh functionality is used to prescribe the movement of the little rectangular bar back and forth on the surface. However, in the level set and phase field methods, the mesh is not moved with the liquid surface.

Here are a few remarks about the model implementation in COMSOL Multiphysics:

  • We add gravity as a source in the momentum equations
  • We add a reference and a point constraint for the pressure field
  • In order to compare results, we select Navier slip conditions for the walls and use a slip length equal to the size of the element length

Simulation Results for the Level Set and Phase Field Methods

The results for the flow field and the shape of the water surface are shown below for the following times: 0.07 s, 0.57 s, and 1.0 s. The flow streamlines and the shape of the water surface are in decent agreement for both methods at the respective times. There is a slight difference in the maximum velocity and in the height of the water surface, which may be explained by the fact that the two methods treat surface tension differently.

There are also some differences in the streamlines after longer times, for example, in the recirculation zones at t = 1.0 s. The phase field method tends to give a slightly calmer behavior of the surface. The level set method computes the surface tension force from the curvature of the surface, which is obtained from the gradient of the level set function. This gives a more “spiky” surface force compared to the smoother force obtained with the phase field method.

Simulation results for the level set method after 0.07 seconds.
Simulation results for the phase field method after 0.07 seconds.
Simulation results for the level set method after 0.57 seconds.
Simulation results for the phase field method after 0.57 seconds.
Simulation results for the level set method after 1 second.
Simulation results for the phase field method after 1 second.

Results from the level set (left) and phase field (right) methods after 0.07 s (top), 0.57 s (center), and 1.0 s (bottom).

We can also visualize the full results from the simulation in an animation for the period from 0 s to 2 s. We can see that although the movement of the rectangular bar causes vigorous stirring, the free surface never breaks.

 

An animation created from the solution of the problem using the phase field method.

The level set and phase field methods also compute the flow field in the air domain above the free liquid surface. We can see that the movement of the rectangular bar leads to a vigorous flow field pattern that is continuous over the phase interface.

Simulation results of the flow field in water and air after 0.57 seconds, computed with the phase field method.
The flow field in the water and air domains as computed with the phase field method after 0.57 s.

If the agitation of the surface is increased so that it becomes more vigorous, then the surface can break up and later coalesce, which is shown in the animation below. This is also one of the benefits with the level set and phase field methods: It is relatively easy to deal with topology changes in the free surface in any of these two methods.

 

The agitation with the rectangle-shaped bar is increased both in frequency and amplitude until a small wave breaks and generates entrapped air pockets in the water phase.

Although the level set and phase field methods are similar, the treatment of the surface tension has a large impact on the stability of the two methods, at least in the implementation in COMSOL Multiphysics. For problems involving a strong effect from surface tension, the phase field method shows better performance regarding computation time compared to the level set method. The reason for this difference is that the computation of the surface curvature in the level set method forces the time-dependent solver to take substantially smaller time steps compared to the phase field method.

In our example, the time steps are on average five times larger in the phase field method, which results in a five-times-larger computation time for the level set method. So, for free surface problems at a small scale and for laminar flows where surface tension has a large impact (for example, in microfluidics), the phase field method is in general a better option.

Adaptive Mesh Refinement and Free Surfaces

In the 2D example that we have investigated in this blog post, the mesh is dense enough over the whole region where the free surface is expected to be found for the level set and phase field method. However, for 3D cases, we cannot always afford this type of resolution. One approach is to use the functionality for adaptive mesh refinement, which automatically creates a denser mesh according to whatever function we choose.

For example, the animations below show the mesh adaption according to the location of the largest gradients in the velocity (left) and according to the largest gradient in the phase field function (right). The plot shows both the air and the water phases. Note that the adaptive meshing with respect to the shear rate gives a refined mesh in the air phase and not that much in the water phase. In this case, we are interested in the water phase so we have to change the error indicator for adaptive meshing by scaling with the phase field function.

 

 

Adaptive mesh refinement with velocity gradients as error indicator (left) and with the phase field gradient as error indicator (right).

Moving Mesh and Free Surfaces

This blog post discusses the two field-based methods for modeling free surfaces: the phase field and the level set methods. In the second blog post in this two-part series, we will compare the modeling of free surfaces using moving mesh with the field-based methods shown here. Stay tuned!

Update: Part two is now live and you can read it here.

Next Steps

Learn more about the specialized features for simulating fluid flow available in the CFD Module, an add-on to COMSOL Multiphysics, by clicking the button below.

Further Reading

Editor’s note: This blog post was updated on 5/25/2018 to include a link to the MPH-file of the example model, which has been added to the Application Gallery.


Comments (20)

Leave a Comment
Log In | Registration
Loading...
maohua lin
maohua lin
May 16, 2018

This is quite interesting case. Can I get this model?

Thanks,

Andrei Kozhevnikov
Andrei Kozhevnikov
May 17, 2018

Indeed, an .mph file would help a lot in understanding of the practical implementation of these models.

Fanny Littmarck
Fanny Littmarck
May 25, 2018

Hello Maohua and Andrei,

The MPH-file has now been added to the Application Gallery: https://www.comsol.com/model/free-surfaces-with-level-set-phase-field-and-moving-mesh-a-comparison-62551

Enjoy!

gaurav singh
gaurav singh
December 17, 2021

I don’t see the mph file there, is it deleted now?

Qiang Huang
Qiang Huang
December 23, 2021

The link you provided does not have the mph file. Is it still available somewhere?

maohua lin
maohua lin
May 25, 2018

Great. Thanks,

Karl-Heinz Leitz
Karl-Heinz Leitz
May 29, 2018

Would it also be possible to provide the model cases with the adaptive mesh refinement? Would it be practicable to efficiently use this approach for 3D multi-Phase CFD simulations?

Anh-Duc LE
Anh-Duc LE
July 5, 2018

This is really interesting. Thanks

Harri Wijaya
Harri Wijaya
August 15, 2018

I run the MPH (phase_field_laminar_wave_II.mph) on COMSOL 5.3a (COMSOL Multiphysics 5.3.1.180) as it is (no modification; just Open MPH –> Study 1 –> Compute). However I encountered the following error:

– Feature: Time-Dependent Solver 1 (sol1/t1)
Nonlinear solver did not converge.
Maximum number of segregated iterations reached.
In Velocity u, Pressure p:
There was an error message from the linear solver.
The relative error (1.3e+002) is greater than the relative tolerance.
Time : 0.062734590029399975
Last time step is not converged.

What could be the cause?
Thanks.

Ed Fontes
Ed Fontes
August 15, 2018

Hi Harri,

I saw that you have sent the same question to our support department. There have been four updates since the version that you are running, we are at Update 4. We have done a few small changes along the way, which may explain why we do not see the problem that you describe in the latest update. Try updating to Update 4. If you cannot do that, for whatever reason that may be, I am sure that our support can help you achieving the same settings in your installed version compared to Update 4.

Best regards,
Ed

Harri Wijaya
Harri Wijaya
August 16, 2018

Thanks a lot, it’s work after updating to Update 4.

KINZAH NOOR
KINZAH NOOR
September 11, 2018

This file is not opening on version 4.4. Is there any solution for it ?

Ed Fontes
Ed Fontes
September 11, 2018

Hi Kinzah,

You have to upgrade to 5.3a, that is the solution.

Best regards,
Ed

Nebyat Weldearegay
Nebyat Weldearegay
November 14, 2018

i really enjoy on the tutorial. but I want to throw a question on modeling using comsol. I kindly request to give its solution.
am working to extract the frequencies of a cantilever beam of partially immersed in liquid however I can’t obtain a result at the final that is it display an error ”Cannot evaluate expression.”
”Undefined variable.
– Variable: comp1.u_fluid
– Geometry: geom1
– Domains: 2-3
Failed to evaluate expression.
– Expression: comp1.u_fluid
– Plot: arwv1 (Arrow volume)”
so how can I fix these errors.

Santanu Das
Santanu Das
December 29, 2018

Dear Ed,

Thank you for sharing the simulation details. It is quite interesting to a newbie like me in this area of multiphase flows. I am working in the area of droplet microfluidics. However, I have a small concern regarding choosing the value of ε and χ. I saw your .mph file where you have chosen values (default/3.2) and 5 respectively. I want to know how one can decide these values based upon the physics.
Thanks in advance.

Best regards.

Ed Fontes
Ed Fontes
January 7, 2019

Dear Santanu,

The value of ε should be set approximately to the same value of the largest elements close to the surface. The proper value of χ is proportional to the speed at the surface and inversely proportional to the coefficient of surface tension. Try with a proportionality factor between 1-50 and make sure that you get enough but not too large diffusivity.

Best regards,
Ed

Ed Fontes
Ed Fontes
January 7, 2019

Dear Nebyat,

Please contact support. It looks to me like you are trying to plot the velocity field in a domain that does not contain a fluid but I could not tell without looking at the model.

Best regards,
Ed

Ed
Ed
January 11, 2021

Very nice blog post! One question about the Phase Initialization step: in the Laminar Flow interface if I modify the acceleration of gravity in the vertical direction from -g_const to, say, -g_const*setp1(t), then the Phase Initialization step complains about the variable t not being defined. Do you understand why? Thanks.

Mamadou Lamine NIANE
Mamadou Lamine NIANE
February 4, 2022

Bonjour,
Merci pour ce article. je rencontre des problèmes pour la mise en œuvre de mon cas. Pourrais je avoir le fichier MPH SVP ?
Merci

Andrei Ludu
Andrei Ludu
May 25, 2023

Hi,
I wonder if anyone knows of an application consisting in a free floating elastic object (disk or prism) over surface water waves. The waves can be generate somewhere far from the object, as simple as possible, like linear, plane, monochromatic waves. Small amplitude, also. Depth of water is not critical (deep is the best). Shape of the basin is not critical either, just way larger than the size of the raft. The equations can be linear. I am interested in the motion and deformation of this floating object.

EXPLORE COMSOL BLOG