PDF
Jet Instability — Moving Mesh1
Introduction
The Marangoni effect results in slip velocity in the tangential direction on a fluid-fluid interface due to gradients in the surface tension coefficient. When the surface tension coefficient is constant, a two-fluid system can exist in static equilibrium. This is because the surface tension force is exactly balanced by a jump in the pressure across the interface. The pressure is discontinuous across the interface, but the velocity field is zero everywhere. The presence of a gradient in the surface tension coefficient means that the flow must be nonstationary. This is due to the fact that the force arising from the variability of the surface tension coefficient acts only in the tangential direction on the interface. This force must be balanced by viscous forces which are only present in a moving fluid. In this example an infinitely long liquid jet breaks up due to a spatially varying surface tension coefficient. Such a situation arises in the fluid jets emitted by continuous inkjet printers as a result of density or thermal gradients. The jet is moving at constant velocity so the problem can be treated in the inertial reference frame in which the jet is initially stationary.
Model Definition
A cylindrical fluid domain with radius 20 microns and 60 microns high contains an initial cylinder of water with radius 5 microns. The surface tension coefficient varies in the axial direction:
where z is the z coordinate, l is the periodicity of the surface tension variation (60 μm) and σ0 is the reference surface tension coefficient (0.07 N/m). The geometry and surface tension coefficient are shown in Figure 1. The spatial variation of the surface tension results in a force in the radial direction at t = 0 s. The two fluids are air and water, whose physical properties are 1.225 kg/m3 and 1000 kg/m3 for the density and 1.79·10−5 Pa·s and 0.001 Pa·s for the dynamic viscosity.
Figure 1: Plot of model geometry (moving mesh) showing surface tension coefficient at t = 0 s.
First, the problem is solved using the Laminar Two-Phase Flow, Moving Mesh interface. In this case the interface is represented by a boundary and has no thickness. Then the problem is formulated on a fixed mesh where the interface is tracked by a level-set function. In this case the thickness of the interface is greatly exaggerated. The Laminar Two-Phase Flow, Level Set interface is used in this case. For practical mesh densities the interface is represented more accurately by the Laminar Two-Phase Flow, Moving Mesh interface, but it cannot handle topological changes and therefore it can only be used prior to the breakup of the droplets.
Problem Formulation — Moving Mesh
Domain Equations
In domains, the incompressible Navier–Stokes equations are solved:
(1)
and
(2)
Here um denotes the mesh velocity and arises from the definition of time derivatives in the coordinate system of the deformed mesh. Poisson’s equation is solved for the mesh displacement.
Boundary Conditions
In the absence of mass flow across the boundary the correct boundary condition on the fluid/fluid interface is
(3)
where the total stress tensor, T is defined for either fluid 1 or fluid 2 as
This boundary condition can be decomposed into a normal component
(4)
and a tangential component
(5)
The term on the right-hand side of Equation 4 is the force per unit area due to local curvature of the interface. The term on the right-hand side of Equation 5 is a tangential stress associated with gradients in the surface tension coefficient. Equation 5 reveals that whenever gradients in the surface tension coefficient exist, the flow must be nonstationary. This is because the pressure is continuous in the tangential direction so the gradient in the surface tension coefficient must be balanced by the tangential component of the viscous stress. A mesh velocity equal to the fluid velocity is imposed on the interface:
(6)
Equation 1 and Equation 2 along with the equation for the mesh displacement describe the evolution of the fluid at the domain level. Equation 3 and Equation 6 are suitable boundary conditions for the problem. These equations and boundary conditions are solved in the Two-Phase Flow, Moving Mesh interface.
Problem Formulation — Level Set Method
When the problem is transformed onto a fixed mesh, the interface is approximated from the spatial derivatives of a level-set function.
Domain Equations
The velocity field and pressure for the liquid phase are described by the Navier–Stokes equations
Here u is the fluid velocity and the last term inside the divergence operator on the right-hand side is the force due to surface tension. Implementing the surface tension force in this way is convenient because the Marangoni effect is taken into account naturally.
In the level set method the fluid-fluid interface is represented as the 0.5 contour of the level-set function. The level-set function, ϕ, represents the volume fraction of the ink (ϕ = 1 in the ink and ϕ = 0 in the air). The transport of the fluid interface separating the two phases is computed by solving the level-set equation:
where ε is a parameter that controls the thickness of the interface.
Boundary Conditions
Periodic boundary conditions are used on the top and bottom boundaries to mimic the effect of the jet being infinitely long in length:
Point Settings
In order to make the pressure unique, the pressure is constrained at a point.
Results and Discussion
Figure 2 shows the solution obtained from the level set method at 6 different time steps. The water (indicated by the red color) initially forms a perfectly axial column but the force due to the variation in surface tension coefficient perturbs the jet. Once this occurs the force due to surface curvature takes over and the jet breaks up into two main lobes and a satellite drop.
Figure 2: Liquid regions of the model (black), obtained from the level set method, shown at 6 different times:0, 6, 9, 10, 11, and 16.4 μs.
Figure 3: Comparison of solutions obtained using level-set (top row) and the moving mesh (bottom row) methods.
Figure 3 compares the solution using the moving mesh and level set methods. Because the moving mesh method cannot handle a topological change, the comparison is only provided during the initial onset of the instability. The agreement between the two methods is good which indicates that the smoothing of the interface over several mesh elements in the level-set method still results in an accurate solution for the evolution of the surface.
Notes About the COMSOL Implementation
The model is straightforward to set up using either the Laminar Two-Phase Flow, Moving Mesh or the Laminar Two-Phase Flow, Level Set interface. The moving mesh simulation can only be used for early time steps, because topological changes cannot be handled by the interface.
Application Library path: Microfluidics_Module/Two-Phase_Flow/jet_instability_mm
Modeling Instructions — Moving Mesh Method
From the File menu, choose New.
New
In the New window, click  Model Wizard.
Model Wizard
1
In the Model Wizard window, click  2D Axisymmetric.
2
In the Select Physics tree, select Fluid Flow > Multiphase Flow > Two-Phase Flow, Moving Mesh > Laminar Two-Phase Flow, Moving Mesh.
3
Click Add.
4
Click  Study.
5
In the Select Study tree, select General Studies > Time Dependent.
6
Click  Done.
Geometry 1
Since the viscosity of the air is negligible compared to that of the ink, only a domain for the ink is necessary.
1
In the Model Builder window, under Component 1 (comp1) click Geometry 1.
2
In the Settings window for Geometry, locate the Units section.
3
From the Length unit list, choose µm.
Rectangle 1 (r1)
1
In the Geometry toolbar, click  Rectangle.
2
In the Settings window for Rectangle, locate the Size and Shape section.
3
In the Width text field, type 5.
4
In the Height text field, type 60.
5
Locate the Position section. In the z text field, type -30.
6
Click  Build All Objects.
Define parameters for the ink properties.
Global Definitions
Parameters 1
1
In the Model Builder window, under Global Definitions click Parameters 1.
2
In the Settings window for Parameters, locate the Parameters section.
3
In the table, enter the following settings:
Name
Expression
Value
Description
rhoink
1e3[kg/m^3]
1000 kg/m³
Density, ink
etaink
1e-3[Pa*s]
0.001 Pa·s
Dynamic viscosity, ink
sigma0
0.07[N/m]
0.07 N/m
Reference surface tension coefficient
Define a spatially varying function for the surface tension coefficient. Making the surface tension dependent on Z rather than z ensures that it is defined in the material frame. This means that as the surface moves in the spatial frame, it retains its original surface tension coefficient.
Definitions
Variables 1
1
In the Definitions toolbar, click  Local Variables.
2
In the Settings window for Variables, locate the Variables section.
3
In the table, enter the following settings:
Name
Expression
Unit
Description
sigma
sigma0*(1-0.2*cos(2*pi*Z/60[um]))
N/m
Surface tension coefficient
Set the fluid viscosity and density to the previously defined parameters.
Laminar Flow (spf)
Fluid Properties 1
1
In the Model Builder window, under Component 1 (comp1) > Laminar Flow (spf) click Fluid Properties 1.
2
In the Settings window for Fluid Properties, locate the Fluid Properties section.
3
From the ρ list, choose User defined. In the associated text field, type rhoink.
4
From the μ list, choose User defined. In the associated text field, type etaink.
The viscosity of the air is neglected and only the air pressure is accounted for in the model. The Free Surface boundary condition can therefore be used. The surface tension is set to the previously defined function.
Free Surface 1
1
In the Physics toolbar, click  Boundaries and choose Free Surface.
2
Select Boundary 4 only.
3
In the Settings window for Free Surface, locate the Surface Tension section.
4
From the Surface tension coefficient list, choose User defined. In the σ text field, type sigma.
The Contact Angle boundary condition is used to apply contact forces to the fluid. The default contact angle of 90 degrees is appropriate for a symmetry boundary.
Contact Angle 1
The Navier Slip boundary condition must be used on all boundaries with a two phase contact. In this case the boundaries are symmetry boundaries.
Wall 1
1
In the Model Builder window, expand the Free Surface 1 node, then click Component 1 (comp1) > Laminar Flow (spf) > Wall 1.
2
In the Settings window for Wall, locate the Boundary Condition section.
3
From the Wall condition list, choose Navier slip.
Boundary conditions for the mesh need to be applied on all the nonfixed exterior boundaries (except for the Free Surface, which automatically includes the mesh constraint).
Component 1 (comp1)
Symmetry/Roller 1
1
In the Model Builder window, right-click Component 1 (comp1) and choose Symmetry/Roller.
2
Select Boundaries 2 and 3 only (top and bottom boundaries).
Mesh 1
Free Triangular 1
In the Mesh toolbar, click  Free Triangular.
Size
1
In the Model Builder window, click Size.
2
In the Settings window for Size, locate the Element Size section.
3
From the Calibrate for list, choose Fluid dynamics.
4
From the Predefined list, choose Extra coarse.
5
In the Model Builder window, right-click Mesh 1 and choose Build All.
Study 1
Step 1: Time Dependent
1
In the Model Builder window, under Study 1 click Step 1: Time Dependent.
2
In the Settings window for Time Dependent, locate the Study Settings section.
3
In the Output times text field, type range(0,0.5e-6,1e-5).
Solution 1 (sol1)
1
In the Study toolbar, click  Show Default Solver.
2
In the Model Builder window, expand the Solution 1 (sol1) node, then click Time-Dependent Solver 1.
3
In the Settings window for Time-Dependent Solver, click to expand the Absolute Tolerance section.
4
In the Variables list box, select Velocity Field (Spatial Frame) (comp1.u).
5
From the Method list, choose Unscaled.
6
Click  Run.
Results
2D Plot Group 5
1
In the Results toolbar, click  2D Plot Group.
2
In the Settings window for 2D Plot Group, locate the Plot Settings section.
3
From the Frame list, choose Spatial  (r, phi, z).
Surface 1
1
Right-click 2D Plot Group 5 and choose Surface.
2
In the Settings window for Surface, locate the Expression section.
3
In the Expression text field, type spf.rho.
4
Locate the Coloring and Style section. From the Coloring list, choose Uniform.
5
From the Color list, choose Black.
6
In the 2D Plot Group 5 toolbar, click  Plot.
Compare the resulting plot with Figure 3 at different time steps.
2D Plot Group 5
1
In the Model Builder window, click 2D Plot Group 5.
2
In the Settings window for 2D Plot Group, locate the Data section.
3
From the Time (s) list, choose 0.
4
In the 2D Plot Group 5 toolbar, click  Plot.
5
Click the  Zoom Extents button in the Graphics toolbar.
6
From the Time (s) list, choose 6E-6.
7
In the 2D Plot Group 5 toolbar, click  Plot.
8
From the Time (s) list, choose 9E-6.
9
In the 2D Plot Group 5 toolbar, click  Plot.
2D Plot Group 6
1
In the Results toolbar, click  2D Plot Group.
2
In the Settings window for 2D Plot Group, locate the Data section.
3
From the Time (s) list, choose 0.
Surface 1
1
Right-click 2D Plot Group 6 and choose Surface.
2
In the Settings window for Surface, locate the Expression section.
3
In the Expression text field, type sigma.
4
In the 2D Plot Group 6 toolbar, click  Plot.
1
This model is courtesy of Prof. E. Furlani, Depts. of Chemical/Biological & Electrical Engineering, SUNY Buffalo, USA (formerly at Eastman Kodak, Rochester, NY, USA).