Skip to content

Adding Phase Field

Jump to bottom
Matthew Hockley edited this page Nov 3, 2020 · 3 revisions

Phase flield (PF) models two immiscible fluids where flow (Navier-Stokes) and electrochemistry can be implemented alongside. The fluid properties vary significantly across the interface where contact angle and wetted walls affect PF evolution over time. An example in microfluidics is creation of oil in water droplets. Adding PF to a Navier-Stokes (NS) problem is relatively simple assuming Creating a CFD problem has been followed and support from PF examples such as flow focusing (NS/PF) and flow around a cylinder (NS/PF).

The difficulties with using PF come when implementing the constants interface_thickness, surface_tension and pf_mobility_coeff. All 3 of the parameters need to be controlled to produce realistic PF conditions. Often surface_tension is retrieved from the literature and interface_thickness is around half the element size of the mesh. surface_tension controls how flexible the interface between the fluids. The interface_thickness controls the size of the scalar, -1 and +1, between the two fluids. The remaining parameter, pf_mobility_coeff, acts a as stabilisation coefficient to keep the droplet sustained and stops the simulation becoming unstable. If the pf_mobility_coeff becomes too low, the phase may never move and if too high, the two fluids diffuse and become indistinguishable. More information is available from Two-Phase Flow Modelling - Comsol.

Parameter Increase Decrease
interface_thickness PF interface too large, increased diffusion Fails to create PF
surface_tension Increases PF front/shape Decreases PF front/shape
pf_mobility_coeff Increases PF mobility Stops or decreases PF mobility

If the parameters do not create a solution, the mesh may required higher element density to allow flexibility of the interface_size and adjustment of the pf_mobility_coeff. The trade off is increased simulation time and mesh size increasing RAM usage which the workstation or simulation may not support depending on the solver used. Ultimately, these parameters control a hyperbolic tangent, tanH, which defines the scalar between '-1' and '+1' as shown below.

problem()

As mentioned in Creating a CFD problem, the most important function is problem(). This function contains a dictionary called parameters which needs the following amended or added...

Parameter Description
enable_PF Enable Phase Field study
contact_angle Phase Field parameter to control the "wetted wall" boundary conditions
interface_thickness Phase Field parameter to control the thickness between both phases in metres
surface_tension Surface tension in Newtons per metre
pf_mobility_coeff Important for forming phase liquids. Acts similar to a stabilisation factor
density Density of the two fluids. In NS problems, both values are the same. In PF, the values are either fluid where the first value is '-1' and the second value is '+1'. In kilograms per metre cubed
viscosity Viscosity of the two fluids. In NS problems, both values are the same. In PF, the values are either fluid where the first value is '-1' and the second value is '+1'. In kilograms per metre a second

In the flow focusing (NS/PF) example, there are 2 calls to FaceLength() function. This is to retrieve both inlets on the X and Y axis. The both X axis inlets have the same H, boundary face length, and Hmin, boundary face minimum co-ordinate, so does not need to be called twice. See image below where X axis inlets, North and South (marked with pink triangle), have the same dimensions and on the same axis compared to Y axis inlet, West.

In the flow focusing (NS/PF) demo, we assign PF '-1' to be oil (left) and '+1' to be water (right).

#         -1          +1
density=[1000,      998.2], # Kg/m3
viscosity=[6.71e-3, 1.003e-3],# Kg/m.s

initialize()

initialize() function allows for assignment of a pre-defined assumption for the problems. Read more in initialize() main. In PF problems, the problem assumes an initial start which is '+1' indicating the second values defined in density and viscosity, see PF problem(). In the flow focusing (NS/PF) example, we assign the mesh to be '+1' indicating entirely water fluid where assignment of PF '-1' is made by the inlet. Additionally, Navier-Stokes uses a pre-determined assumption using velocity_init.

    w_init_field = dict()
    if not restart_folder:
        if enable_NS:
            try:
                subspace = field_to_subspace["u"].collapse()
            except:
                subspace = field_to_subspace["u"]
            #length inlet, water inflow,
            #   X (0) or Y (1) dir flow, 
            #   Positive/neg flow along axis (+1/-1),
            #   Hmin value
            u_init = velocity_init(H[0], inlet_velocity, 0, 1, Hmin[0])
            w_init_field["u"] = df.interpolate(u_init, subspace)
        # Phase field
        if enable_PF:
            w_init_field["phi"] = df.interpolate(
                df.Constant(1.),
                field_to_subspace["phi"].collapse())

create_bcs()

Boundaries follow a similar definition as in 'create_bcs()' main where enable_PF adds expression for PF inlet/starting point and boundary specific constants such as wetting. The PF inlet/starting point controls where the immiscible fluid, '-1', begins in the channel and mixes with the bulk fluid, '+1'. This is defined using 'tanH' which creates a hyperbolic tangent from '-1' to '+1'. The boundary conditions are important when studying microfluidic flow to understand the flow focusing and "pinching" dynamics. Example below is from flow focusing (NS/PF).

    # Phase field uses an expression `tanH` which defines PF drop off
    if enable_PF:
        phi_expr = df.Expression(
            "tanh((abs((x[1]-Hmin)-H/2)-H/16)/(sqrt(2)*eps))",
            H=H[0], Hmin = Hmin[0], eps=interface_thickness,
            degree=2)
        phi_inlet = Fixed(phi_expr)

        ## PF Fixed across boundary
        # as no water can enter oil inlet
        #   and vice-versa
        bcs["inlet"]["phi"] = Fixed(df.Constant(-1.))
        bcs["inletT"]["phi"] = Fixed(df.Constant(1.))
        bcs["inletB"]["phi"] = Fixed(df.Constant(1.))
        ## Add contact angle to NS No-Slip Boudnaries
        bcs["wallLR"]["phi"] = ContactAngle(contact_angle)
        bcs["wallRL"]["phi"] = ContactAngle(contact_angle)

initial_phasefield()

The initial_phasefield() controls the size and shape of the PF at the beginning of the problem. Often this defines the PF as a droplet in flow as PF entering via an inlet is defined solely by the boundary condition in create_bcs.

Parameter Description
x0 X axis position for centre of PF
y0 Y axis position for centre of PF
rad Radius of PF droplet
eps Thickness of interface. Also known as interface_thickness
function_space function space used to apply PF on to 
    expr_str = "tanh((x[0]-x0)/(sqrt(2)*eps))"
    phi_init_expr = df.Expression(expr_str, x0=x0, y0=y0, rad=rad,
                                  eps=eps, degree=2)
    phi_init = df.interpolate(phi_init_expr, function_space)
    return phi_init

pf_mobility()

The PF needs the 3 parameters interface_thickness, surface_tension and pf_mobility_coeff to create a stable state. The pf_mobility() allows for finer control of the pf_mobility_coeff using phi a trial function for PF and gamma which is pf_mobility_coeff to adapt PF mobility over time as the PF changes. Depending on evolution of the PF problem, the gamma may need to be adjusted to increase or decrease diffusion and mobility. In the flow focusing (NS/PF), gamma is returned as pf_mobility_coeff.

Parameter Description
phi Current trial function for PF from simulation
gamma Coefficient to stabilise PF. Also known as pf_mobility_coeff

Menu

Clone this wiki locally