Modelling of complex interfaces for pendant drop experiments

Interfaces of fluid-fluid systems play an important role in the stability of foams and emulsions in chemistry, biology, consumer products, and foods. For most applications, surface active agents are added and adsorbed onto the interface to enhance stability, making the rheological behavior of the interface more complex. To understand the phenomena of these complex interfaces, various techniques are used to determine the interfacial properties. One of the most popular methods is the pendant drop technique. From the equilibrium state of the pendant drop, the inter-facial tension of a system can be obtained quite easily in the absence of surface active agents. But when complex viscoelastic interfacial characteristics are considered, in particular in oscillatory measurements, interfacial constitutive relations need to be defined. Interfaces containing proteins, particles or Langmuir monolayers formed by insoluble low weight surfactants appear to act like viscoelastic solid membranes. In this work, a two-dimensional axisymmetric finite element model is designed to study the behavior of complex interfaces in pendant drop experiments. The bulk fluid consists of a Newtonian fluid, while the interface behaves according to the Kelvin-Voigt model as elastic interfacial


Introduction
Complex fluid-fluid interfaces are found everywhere among us even though we are not always aware of them.One of the most important characteristics of an interface is the interfacial tension.The interfacial tension is defined as the interfacial free energy per unit area, which is the minimum amount of work required to create or expand the interface by a unit area.The interfacial tension, which can be seen as a material parameter of the fluid-fluid system, is dependent on temperature, but also on miscibility of the components, added particles, or surfactants.The recent review by Fuller and Vermant (2011) summarizes different examples and applications in daily life and industry and outlines the relevance of understanding the dynamics and physics of these complex interfacial transport phenomena.For example, interfacial stabilizers, surfactant molecules, proteins, or particles preferentially residing within the interfacial region are often used to lower the interfacial energy yielding, for example, stability to foams and emulsions, with applications for consumer care products, foods, paints, and other chemical products.Furthermore, these systems are found in living organisms as well.
Surface active agents, or "surfactants", are molecules that can adsorb onto the interface of a fluid-fluid system and alter the interfacial free energy.Rosen (2004) writes that they have a characteristic amphipathic structure, i.e., a molecular structure containing both a lyophilic group and a lyophobic group.The lyophilic group has strong attraction for the solvent, while the lyophobic group has little attraction for the solvent.At the interface, the surfactants orientate themselves such that the lyophobic groups have minimal contact with the solvent.Usually, this reduces the amount of work required to create or expand the interface and thus decreases the interfacial tension.
Besides classification of surface active agents based on the properties of the hydrophilic group, they can be characterized by the difference in measurable dynamic properties.Roughly, three groups can be distinguished: low molecular weight surfactants, proteins, and particles.The dilatational response dominates the weaker shear phenomena for all three groups.Bos and Van Vliet (2001) concluded that in dilatational rheology, the elasticity is lower if the molecules are present in the solvent as well, compared with insoluble surfactants, because the exchange of molecules between the surface layer and the solvent causes the surfactant layer to behave more viscous.Kotula and Anna (2015) found that the behavior of the interface in the presence of soluble surfactants depends on the surfactant isotherm parameters, transport parameters, and the geometry of the interface.In general, it holds that for more condensed, interacting, or solid-like interfaces, a more explicit viscoelastic response in shear is expected, according to Erni et al. (2012).In particular, the research within the field of interfacial rheology aims at quantifying the mechanical properties of complex fluidfluid interfaces.More specifically, the relation between deformation and induced stresses is examined.
A variety of experimental techniques are developed to measure several interfacial properties of systems with fluidfluid interfaces.Most procedures to measure the interfacial tension are based on a thermodynamical equilibrium state.An overview of these methods is given by Drelich et al. (2002).To measure other rheological characteristics of interfacial layers, the techniques are of a more dynamic nature.Derkach et al. (2009) wrote a review on methods of measuring rheological properties of interfacial layers.One technique used for both static and dynamic experiments is the pendant drop method and is widely used to measure the interfacial properties of a system.In particular, the radius of curvature at the apex is the key property of pendant drops together with the droplet volume that is then used to obtain the surface tension (Yeow et al. 2008).During static experiments, the equilibrium state between the surface tension forces and gravitational forces is considered.The dynamic interface characteristics are tested under small amplitude oscillatory conditions.By introducing surface active agents to the system, the interface properties alter, i.e., the behavior of the interface may become more viscoelastic.
To analyze most interfacial rheological experiments, twodimensional surface constitutive models are used, which are generalized bulk constitutive models.For purely viscous interfaces, the Boussinesq-Scriven model is the twodimensional analogue of the Newtonian fluid, described by Scriven (1960).On the other hand, for purely elastic interfaces, a generalization of Hooke's law is often used, as done by Barthès-Biesel and Rallison (1981).Since interfaces with surface active agents behave viscoelastic, constitutive relations like the generalized Maxwell model for interfaces, the Kelvin-Voigt model revised for interfaces, and more elaborate spring-dashpot models can be used to describe these characteristics even more effectively.An extensive review is written by Sagis (2011) on interfacial rheological experiments and the constitutive models used to analyze the behavior of the interfaces.
To computationally describe interfaces and transport on interfaces, different types of models can be used.In particular, the case of highly viscous drops with Newtonian components is studied extensively using the boundary integral model.Research groups of Pozrikis, Zinchenko, and Loewenberg have pioneered the field and demonstrated for a wide variety of flow conditions the complexity introduced by an interface even with a constant interfacial tension (Pozrikidis 1992;Zinchenko and Davis 2006;Zinchenko et al. 1997;Loewenberg and Hinch 1997).Stebe showed by boundary integral modeling that the addition of a surfactant molecule has a large effect on the modes of breakup in extensional flow (Eggleton et al. 2001).Bazhlekov et al. (2006) studied the effect of surfactants on true threedimensional drop deformation in shear flow, and the authors introduced a morphology diagram that links breakup modes to surfactant coverage.Although boundary integral modeling is very rich and powerful, it is limited to linear systems.Several groups have applied interface capturing techniques that in principle are more general, but they do not provide an explicit description of the interface.Examples are the front capturing work of James and Lowengrub (2004) who developed a surfactant-conserving volume-offluid method for interfacial flows with insoluble surfactant.Other approaches are by Voigt and co-workers, and these authors developed a Navier-Stokes-Cahn-Hilliard model for the macroscopic two-phase flow system that is combined with a surface phase-field-crystal model for the microscopic colloidal system along the interface (Aland et al. 2011).
In this work, an axisymmetric finite element interface tracking model is designed to investigate the behavior of complex viscoelastic interfaces for pendant drop experiments.

Balance equations and constitutive model
In this work, an isothermal flow of an incompressible fluid is assumed and the system is described by the following set of equations for the momentum balance and mass balance, respectively, where ρ is the density, u the velocity, and f an external volume force acting on the fluid, which is assumed to be a gravitational force f = ρg applied in z-direction in this work.The Cauchy stress tensor σ is defined as with pressure p and extra stress tensor τ .It is assumed that the bulk fluid is Newtonian; therefore, the extra stress tensor τ is given by where D = 1/2 ∇u + (∇u) T is the symmetric part of the velocity gradient tensor and η is the viscosity.These equations hold in geometry , as shown in Fig. 1.

Boundary, symmetry, and interface conditions
At the inlet of the capillary in , two types of boundary conditions are applied: either a parabolic velocity profile is applied instantaneously by prescribing a certain flow rate Q or a certain pressure p is imposed A no-slip boundary condition is prescribed at the capillary wall wall u = 0 on wall .( To imply symmetry, the velocity in radial direction at the symmetry axis symmetry is set to zero The boundary condition at the interface is a Neumann boundary condition where ∇ s is the surface gradient operator, ς is the interfacial stress tensor, n is the outwardly directed unit normal vector, and the outside pressure p out = 0.The surface gradient operator is defined as ∇ s = I s • ∇ and I s = I − nn is the second order unit surface dyadic tensor.Depending on the rheology of the interface, different constitutive models for the interfacial stress tensor ς can be used.In the case of linear behavior, a constant interfacial stress can be assumed where γ is the static equilibrium value of the interfacial tension.The surface tension γ is kept constant over the interface in time. For complex interfaces with non-linear behavior, elastic and viscous terms need to be included into the interfacial stress tensor ς, for example by using a generalized Maxwell model for interfaces or a Kelvin-Voigt interface model.Insoluble low weight surfactants forming Langmuir monolayers, proteins, and particles typically have the response of a viscoelastic solid.The Kelvin-Voigt model is a suitable model for these systems with dominating elasticity and is a sum of the Boussinesq-Scriven constitutive model for purely viscous fluids (Scriven 1960) and the generalized Hooke's law for purely elastic solids (Barthès-Biesel and Rallison 1981).The model can be made quasi-linear by using a finite strain tensor which can describe the stress response for arbitrary large deformations.Verwijlen et al. (2014) used the Green-Lagrange finite strain tensor E s = 1/2(C s − I s ) for one formulation of the quasi-linear two-dimensional Kelvin-Voigt model and the Hencky strain tensor H s = 1/2 ln(C s ) for a second formulation.Herein, C s = F T s • F s is the Cauchy-Green deformation tensor of the interface with F s = ∇ s,0 x T is the surface deformation gradient tensor, expressing the deformation of a certain material point x on the interface with respect to the initial undeformed state, indicated by subscript 0. From application of these two formulations of the model to both pure dilatation and simple shear, the formulation using the Hencky strain tensor appears to be more appropriate for dilatational deformation, whereas the formulation using the Green-Lagrange strain tensor is more suitable for shear deformation.In both formulations, the shear and dilatational contributions are separated.Combining the dilatational contribution using the Hencky strain tensor with the shear contribution using the Green-Lagrange strain tensor, the quasi-linear twodimensional Kelvin-Voigt viscoelastic constitutive equation becomes Herein, κ is the surface dilatational viscosity, μ the surface shear viscosity, K the surface dilatational elasticity, G the surface shear elasticity, and J = det (C s ).Furthermore, J 1 2 represents the change in interfacial area A/A 0 .Particle laden interfaces may exhibit bending moduli that are not marginal, as shown by Yunker et al. (2012).These bending stresses arise when the thickness of the interface is nonzero.In this work, however, the influence of bending stresses is assumed to be negligible compared to the surface tension, dilatational, and shear stresses.Therefore, the interfaces are modeled with a zero thickness.
If a complex interface is assumed, the quasi-linear Kelvin-Voigt model of Eq. ( 11) is prescribed to the interface.The value of the surface tension γ may be different from the value of pure water.All parameters of the viscoelastic interface are constants.Furthermore, other interfacial dynamics like Marangoni flows at the interface and diffusion of material between the bulk and the interface are neglected.

Evolution of the pendant drop in time
The motion of the interface is determined using a FEMbased interface tracking method.A moving curvilinear coordinate system describes the interface according to where x is the function that maps the curvilinear coordinates ξ = (ξ 1 , ξ 2 ) onto the spatial coordinates x of the interface.The interface is tracked in a Lagrangian way.Therefore, the velocity of the interface is defined as Herein, u is the material velocity at the interface.

Surface strain
To quantify the deformation of the interface, the strain is calculated in tangential and circumferential direction, using the stretch ratio of the interface.The stretch ratio in the tangential direction is given by Herein, t is the unit tangential vector.The stretch ratio in circumferential direction is determined by where e θθ is the unit vector in circumferential direction, r is the current radius, and r 0 is the initial radius of the interface, respectively.
The stretch ratios are used to calculate the linear strain, defined as both in tangential and circumferential direction.

Weak form
As stated in the introduction, the finite element method is applied to solve the governing equations with appropriate boundary conditions as discussed in the previous section.In order to derive the weak form, the momentum balance and mass balance are multiplied with the test functions v and q in the appropriate function spaces Note that (•, •) defines the inner product on in the usual manner by for scalars a and c, vectors a and c, and tensors A and C. For a tensor A and vector a, the product rule holds and the surface divergence theorem, as described by Weatherburn (1955), is given by where b = t × n is the binormal and ∂ the boundary of the interface.In our geometry, the boundary of the interface ∂ is the connection of interface to the capillary wall wall .Using these two theorems, the boundary condition on the interface in the weak formulation can be written as The interfacial stress tensor ς is always tangential to the surface, as a result of which n • ς = 0 and the second term of the right hand side is zero.
The complete weak form including the boundary terms becomes: Find u and p such that for all possible test functions v and q and where D v = 1/2 ∇v + (∇v) T , (•, •) is an inner product defined on interface and (•, •) ∂ is an inner product defined on the boundary of the interface ∂ .Due to the kinematic boundary condition on ∂ at the adhesion of interface to the capillary wall, the test function v = 0, and due to the lack of area at the tip of the drop on the other side of ∂ , The solution of the instationary term of the momentum balance for the first time step is found using a backwards Euler integration scheme, which is of first order, using the solutions of the velocity u at time t n+1 and t n , respectively.The solution for all subsequent time steps will be of second order, using the second order backwards differencing method (BDF2) Herein, the solutions of the velocity u at time t n+1 , t n , and t n−1 are used.The convection term is solved iteratively.For the first iteration step, the Picard fix point iteration is used, and Newton's method is used for the following iteration steps.

Semi-implicit time integration of the surface tension
Assuming an equilibrium between the viscous and interfacial forces, a characteristic capillary time scale can be defined as t capillary = ηR/γ .Herein, η is the bulk viscosity, R is a characteristic length scale, and γ is the interfacial tension.On element level, the capillary time scale of the problem is t capillary = ηR/γ n elem = O(10 −7 s), where R = R o and n elem is the number of elements on the interface.The explicit numerical scheme is stable for time steps smaller than the capillary time scale.In order to be able to use larger time steps, the problem is written to a semiimplicit time integration scheme, as presented by Hysing (2006).Taking the contribution of the constant interfacial stress from Eq. ( 10) and using definitions from differential geometry, this equation can be rewritten to where x is the interface position.By writing the new interface position as a function of the old position, proposed in the work of Bänsch (2001), the time discretization can be made semi-implicit Herein, t is the time step and u n+1 is the velocity field at the new time.The interfacial stress tensor for the contribution of the surface tension becomes Thus, the second term of the right hand side of Eq. ( 29) is added to the left hand side of the system.

Interface mesh
The movement of the interface is Lagrange based where dx /dt is the velocity of the interface and u is the material velocity at the interface.The discretization is done using a first order scheme Herein, x n+1 and x n are the coordinates of the interface at times t n+1 and t n , respectively.

Bulk mesh
Since the problem has moving boundaries, the mesh is not stationary.To take the motion of the mesh into account, an Arbitrary Lagrangian Eulerian (ALE) method is applied.The material derivative of the momentum balance is changed by subtracting the grid velocity from the flow velocity in the convective term where δ(•)/δt is the grid derivative and u grid is the grid velocity.
The grid velocity u grid for time t n+1 is calculated using a time discretization based on an Euler backwards differencing scheme for the first time step and based on a BDF2 scheme for all subsequent time steps Herein, x n+1 , x n , and x n−1 are the mesh coordinates on time t n+1 , t n , and t n−1 , respectively.
The displacement of the mesh coordinates x is necessary to find the new mesh coordinates x n+1 = x n + x and is calculated by solving a Laplace equation on the mesh.The displacement of the interface x = x n+1 − x n is used as a boundary condition for interface , while the other boundary conditions are set to zero, as shown in Fig. 1.The following set of equations is solved x = 0 on in , wall , and symmetry (37) for both zand r-direction separately.Coefficient α is a constant coefficient per element with the value α = 1/A e , where A e is the area of the element (Hulsen 2015).

Remeshing, refining, and projection
When the geometry is highly deformed, the mesh may become too distorted to produce accurate solutions.To quantify the quality of the mesh, the deformation of each element is monitored, using the following equations where A e and A e 0 are the element areas of the current mesh and undeformed mesh and S e and S e 0 are the element aspect ratios of the current mesh and undeformed mesh, respectively.The aspect ratio of an element is defined as Herein, L e max is the maximum length of the sides of an element.Remeshing is performed if either f e 1 > 0.2 or f e 2 > 0.2 (Hulsen 2015).Furthermore, the elongation of the side of each element on the interface and symmetry axis is monitored by where L e and L e 0 are the length of the sides of an element at these two curves of the current mesh and undeformed mesh, respectively.An element is refined if f e 3 > 1.4.The refinement is performed by adding a new node on the curve, exactly in the center between the nodes in the reference element.Subsequently, the geometry of the bulk is remeshed.
The solutions u n and u n−1 on the new mesh are necessary to solve the implicit time integration of the Navier-Stokes equation.A projection problem is solved to acquire these solutions on the new mesh.For this projection problem, where the values of u old are found in the integration points defined on the new mesh.acquire meshes with an increasing number of nodes, which can be found in Table 1.The mesh is subject to remeshing and refining during the calculations.
The relative L 2 -errors are calculated on line , given by (z, r) = (s, 0.1), where s is a coordinate along , as shown in Fig. 2. The relative L 2 -error u and p for the velocities and pressures respectively are defined as where u h and p h are the solutions on one of the meshes given in Table 1 and u * and p * are reference solutions computed on mesh M6 * .The integrals of Eqs. ( 44) and ( 45) are computed by dividing line into 10,000 intervals and using a midpoint rule on each interval.The mesh-convergence plots are shown in Fig. 3 and show that both the velocity and pressure converge with at least order two.

Comparison with recent results from literature
To validate the results in the case of an interface with constant surface tension, the results of this study are compared with the results of Dieter-Kissling et al. ( 2014).In their study, the applicability limit of the pendant drop measuring method at high inflow rates is studied.For high flow rates, the results of the Profile Analysis Tensiometry instrument appear to be invalid.Using their geometry, as shown in  The volume of the drop, as depicted in Fig. 5 in gray, is calculated by The apex length is the z-coordinate of the apex of the drop measured from the origin.As shown in Fig. 4, the

Necking
The simulations as described in the previous section are performed up until necking and the detachment of the drop.The results for the three different flow rates are shown in Fig. 6.
From approximately V = 13 mm3 onward, the curves become more steep.This is where the necking occurs.In Fig. 7, the deformation of the mesh during the necking process is depicted.
For a better understanding of the necking process, the radius of the neck is plotted in time, as shown in Fig. 8, where t = 0 s is the onset of the necking process.
Detachment of the drop occurs theoretically according to the following force balance (Garandet et al. 1994) where p is the pressure, R is the wetted radius of the capillary, V is the volume of the drop, ρ is the density, g is  9 Mesh and velocity for the reference state system at time t = 0.01 s (a), t = 0.50 s (b), t = 1.00 s (c), t = 1.01 s (d), t = 1.04 s (e), t = 1.05 s (f), t = 2.00 s (g), and t = 2.01 s (h), the velocity is in mm/s the gravitational acceleration, and γ is the surface tension.The pressure at the end of the capillary is p = 68 Pa.From this follows that the theoretical volume for detachment of the drop is equal to V = 13.47 mm 3 , which is similar to the volume at which the necking occurs in the three simulations.

Viscoelastic interface
To investigate the influence of the five parameters of the viscoelastic constitutive model of Eq. ( 11), a parameter study is performed (Table 2).The ranges of the dilatational and shear viscosities and elasticities are taken from Cantat et al. (2013) and experiments performed by the group of Vermant (Verwijlen 2013).The dilatational viscosity is typically one or two orders of magnitude larger than the shear viscosity.
For an ideal elastic interface, the surface dilatational elasticity and surface shear elasticity are related to the surface Young's modulus Y according to where ν s is the surface Poisson's ratio (Sagis 2011).Theoretically, the surface Poisson's ratio can vary between −1 ≤ ν s < 1.However, for conventional materials, it holds that 0 ≤ ν s < 1.The ratio of the shear elasticity G and the dilatational elasticity K will therefore be 0 < G/K ≤ 1 .
Fig. 10 Velocity in the z-direction for the reference state system at time t = 0.01 s, t = 0.50 s, t = 1.00 s, t = 1.01 s, t = 1.04 s, t = 1.05 s, t = 2.00 s, and t = 2.01 s, the velocity is in mm/s The parameters in the gray boxes in Table 2 are chosen to be the reference state, where the ratios μ/κ = 1 and G/K = 1 are upper limits.For every simulation, one parameter is changed from its reference state.
For 0 ≤ t ≤ 1 s, a parabolic velocity field with flow rate Q = 2 mm 3 /s is given as input.The corresponding Reynolds number is Re = 5.66.Subsequently, for 1 < t ≤ 2 s, a flow rate Q = 0 mm 3 /s is given as input, enabling relaxation of the system.Finally, for t > 2 s, a pressure p = 0 Pa is set at the inlet of the capillary in , facilitating outflow of the fluid.For the calculations, the time step t = 1 • 10 −3 s and gravitational forces are taken into account.
The mesh and the flow field of the reference state at different times are shown in Fig. 9.
In pictures (a), (b), and (c) of Fig. 9, the inflow profile is shown.In the picture (d), a vortex flowing counterclockwise can be found.This is a result of the inertial forces in the fluid.In pictures (e) and (f), a flow along the interface towards the apex is shown, resulting in a second vortex close to the apex flowing clockwise.The first vortex vanishes and the second vortex dominates the flow.During the relaxation regime, the absolute velocity decreases rapidly.In picture (h), the outflow from in due to the pressure difference is visible.In Figs 10, 11, and 12, the velocity in zand r-direction, and the pressure is shown.
From simulations with the same volume using different inflow rates, the occurrence of the clockwise vortex around the apex during the relaxation regime turns out to depend on the initial flow rate; for smaller flow rates, the vortex flowing clockwise appears sooner.11 Velocity in r-direction for the reference state system at time t = 0.01 s, t = 0.50 s, t = 1.00 s, t = 1.01 s, t = 1.04 s, t = 1.05 s, t = 2.00 s, and t = 2.01 s, the velocity is in mm/s

Surface tension
From the Young-Laplace equation it is known that the shape of a drop depends on the equilibrium between the gravitational forces and the surface tension.The shape of a drop should be more spherical for higher values of the surface tension, i.e., the apex length should be smaller.This behavior is found for both the inflation and relaxation of the drop as shown in Fig. 13, where the apex length in time is given for different values of the surface tension γ .
For values of the surface tension γ ≥ 20 μN/mm, the shape of the drop remains constant during relaxation, while the drop with surface tension γ = 10 μN/mm becomes more elongated and the onset of necking of the drop is visible in the relaxation time from t = 1.80 s.According to Eq. ( 47), using a pressure p = 14.6 Pa at the end of the capillary, the drop will start to detach at a theoretical volume of V = 1.67 mm 3 .The volume of the drop during relaxation is 2.09 mm 3 .
Once the pressure at the inlet is set to p = 0 Pa, an outflow is generated due to the pressure difference.As shown in Fig. 13, the entire volume of the drop flows back into the capillary.For higher values of the surface tension γ , this occurs faster.The velocity of the decrease in apex length is approximately proportional to the surface tension γ .
The clockwise vortex around the apex, as described during the relaxation regime in the simulation of the reference state, occurs sooner for higher values of the surface tension γ .In Table 3, the maximum values of the surface strain in both tangential and circumferential direction are given at the end of the inflow regime t = 1 s and at the end of the relaxation regime t = 2 s for different values of the surface tension γ .For smaller values of the surface tension, the maximum tangential surface strain is higher.This is due to the elongation of the drop for these small surface tension values, driven by the gravitational forces.This is shown for both the end of the inflow time and the end of the relaxation time.During the relaxation regime, the maximum tangential surface strain increases and the maximum circumferential surface strain decreases.This is due to the vortex in the bulk, moving the interface towards the apex.The influence of the surface tension on the circumferential surface strain is small.
Looking at the distribution of the surface strain, the tangential surface strain is maximum at the attachment to the capillary, while the circumferential strain is maximum between the equator and the apex of the drop.This is shown in Fig. 14.
In these calculations, the parameters for the viscoelastic behavior are nonzero.When only the surface tension is taken into account, i.e., the parameters for the viscoelastic behavior are set to zero, holds for a surface tension γ = 50 μN/mm at t = 1 s a maximum tangential surface strain of ε t = 4.55 and a maximum circumferential surface strain of ε c = 2.13.At time t = 2 s, the maximum tangential surface strain is ε t = 4.52 and the maximum circumferential surface strain is ε c = 2.13.

Surface viscosity
Viscosity is a measure of resistance to deformation.For higher viscosity, the shape of the drop should be more spherical, but should become equal for all different values of the viscosity after total relaxation.In Fig. 15, where the apex length in time is depicted for different values of the surface dilatational viscosity κ, this behavior is shown for both regimes.From Fig. 16, where the apex length in time is depicted for different values of the surface shear viscosity μ, it can be seen that the difference in behavior for the used range of values of the shear viscosity μ is negligible.
During the relaxation time, for the dilatational viscosity κ ≥ 50 μNs/mm, the relaxation occurs fast, as can be seen in Fig. 15.
Once the pressure at the inlet is set to p = 0 Pa, an outflow is generated due to the pressure difference.For the dilatational viscosity κ ≥ 50 μNs/mm, the shape becomes more elongated although the volume decreases.This is due to a flow in negative radial direction in the region of the drop just outside the capillary and the resistance to deformation of the interface.This radial flow initializes necking and "collapsing" of the drop, while the interface withstands the change in area.The flow field for κ = 100 μNs/mm at the start of the outflow regime is shown in Fig. 17.For higher Apex length [mm] Apex length vs. time − κ 1 μNs/mm 10 μNs/mm 50 μNs/mm 100 μNs/mm shear viscosity μ, the decrease in volume is slower than for lower shear viscosity μ.
The clockwise vortex around the apex, as described during the relaxation regime in the simulation of the reference state, occurs later for larger values of the dilatational viscosity κ, the absolute value of the velocity is lower and originates at a larger distance from the interface.The different values for the shear viscosity μ have little influence on the flow field of the drop.
The maximum values of the surface strain in both tangential and circumferential direction are given at the end of the inflow regime t = 1 s and at the end of the relaxation regime t = 2 s for different values of the dilatational viscosity κ in Table 4 and for different values of the shear viscosity μ in Table 5.The influence of the dilatational viscosity on the tangential surface strain is small.For larger values of the dilatational viscosity, the maximum circumferential surface strain is larger.This is shown for both the end of the inflow time and the end of the relaxation time.For higher values of the shear viscosity, the tangential surface strain is smaller and the circumferential surface strain is larger.This is shown for both the end of the inflow time and the end of the relaxation time as well.During the relaxation regime, the vortex in the bulk moves the interface toward the apex, resulting in an increase of the maximum tangential surface strain and a decrease of the maximum circumferential surface strain.This is observed for both the different values of the dilatational viscosity κ and the different values of the shear viscosity μ.

Surface elasticity
Elasticity is the tendency of materials to return to their original shape, i.e., their state of zero stress.In this study, this state is given by the geometry t=0 of Fig. 1.For higher elasticity, the shape of the drop should be more spherical during inflation.During the relaxation time, the shape should not change anymore.This behavior is shown in Fig. 18, where the apex length in time for different values of Table 4 Maximum surface strain for different values of the dilatational viscosity κ at the end of the inflow regime and the relaxation regime the surface dilatational elasticity K is displayed.In Fig. 19, the apex length in time is shown for different values of the surface shear elasticity G.As depicted, the difference in behavior for the used range of values of the shear elasticity G is negligible.Since a larger pressure is needed for higher values of the dilatational elasticity K, the outflow velocity is higher once the pressure is set to p = 0 Pa at the inlet.This results in a faster decrease in volume and apex length.For values of the dilatational elasticity K ≤ 10 μN/mm, the surface tension seems to dominate the outflow behavior, since the volume Table 5 Maximum surface strain for different values of the shear viscosity μ at the end of the inflow regime and the relaxation regime of the drop decreases to zero.For dilatational elasticity K = 100 μN/mm, an equilibrium between the elasticity and the surface tension is reached before the volume of the drop is decreased to zero.As shown in Fig. 19, the difference in outflow behavior for the used range of values of the shear elasticity G is negligible.
The clockwise vortex around the apex, as described during the relaxation regime in the simulation of the reference state, occurs later for larger values of the dilatational elasticity K and the absolute value of the velocity is lower.For K = 100 μN/mm, the velocity is already zero before the vortex can occur.For larger values of the shear elasticity G, the clockwise vortex around the apex occurs sooner and the absolute value of the velocity is higher.
In Table 6, the maximum values of the surface strain in both tangential and circumferential directions are given at the end of the inflow regime t = 1 s and at the end of the relaxation regime t = 2 s for different values of the dilatational elasticity K and for different values of the shear elasticity G this is shown in Table 7.The influence of both the dilatational and shear elasticity on the tangential surface strain is small.For larger values of the dilatational elasticity, the maximum circumferential surface strain is larger, while for larger values of the shear elasticity, the maximum circumferential surface strain is smaller.During the relaxation regime, the vortex in the bulk moves the interface toward the apex, resulting in an increase of the maximum tangential surface strain and a decrease of the maximum circumferential surface strain.This is observed for both the different  6 Maximum surface strain for different values of the dilatational elasticity K at the end of the inflow regime and the relaxation regime

Elevated values of the interface parameters
The same flow regime as described in the previous section is carried out for elevated values of the surface dilatational and shear viscosity and elasticity.The surface tension is kept at γ = 50 μN/mm.For simulations with elevated surface viscosities, the dilatational viscosity is set to κ = 100 μNs/mm, and the shear viscosity is varied between 1 ≤ μ ≤ 100 μNs/mm, while the surface elasticities are kept at K = G = 1 μN/mm.Similarly, for simulations with elevated surface elasticities, the dilatational elasticity is set to K = 100 μN/mm, and the shear elasticity is varied between 1 ≤ G ≤ 100 μN/mm, while the surface viscosities are kept at κ = μ = 1 μNs/mm.For flow situations where the drop undergoes very large deformation, the results obtained are limited by the quasi-linear constitutive equation, and future work is focussed on developing more realistic interfacial constitutive relations.

Surface viscosity
As stated in "Surface viscosity", for higher surface viscosity, the shape of the drop should be more spherical during the inflow regime, but should become equal for all different values of the viscosity after total relaxation.In Fig. 20, Table 7 Maximum surface strain for different values of the shear elasticity G at the end of the inflow regime and the relaxation regime where the apex length in time is depicted for different values of the surface shear viscosity μ, this behavior is shown for the inflow regime.For values of μ ≥ 10 μNs/mm, it can be seen that a relaxation time of 1 s is not sufficient for the flow to fully relax.During the outflow regime, only for μ = 1 μNs/mm, the figure shows that the shape of the drop becomes more elongated although the volume decreases.
For more elevated values, the surface shear viscosity might prevent the "collapsing" of the drop.

Surface elasticity
As shown in Fig. 21, where the apex length in time is depicted for different values of the surface shear elasticity, the drop becomes more elongated during the inflow regime for increasing values of the elasticity.This is contradictory to the expectation that the shape of the drop should be more spherical during inflation, as stated in "Surface elasticity".For G = 100 μN/mm, a sharp slope transition can be observed.In this case, the resistance of interface deformation in r-direction is more profound than in z-direction, resulting initially in a flow in z-direction only, before flowing in r-direction.During the relaxation time, the shape of the drop does not change anymore, as expected.For shear elasticity G = 1 μN/mm, the equilibrium between the elasticity and the surface tension is reached before the volume of the drop is decreased to zero, as is seen in "Surface elasticity" as well.For values of G ≥ 10 μN/mm, the shape of the drop becomes more elongated during the outflow regime, and the simulations end before the prescribed time is over.This is due to the numerical boundary condition that states that the interface is fixed to the outer edge of the capillary; the fluid layer on the wetted tip of the capillary becomes too thin to be properly calculated.

Oscillatory measurements
To connect the presented model to experimental methods for testing the interface characteristics, a small amplitude oscillatory measurement is simulated with different interfacial properties.For 0 ≤ t ≤ 0.5 s, a parabolic velocity field with flow rate Q = 2 mm 3 /s is given as input at in .After that, a flow rate Q = 0 mm 3 /s is given as input for 0.5 < t ≤ 1 s, enabling relaxation of the system.For the next 12π seconds, a parabolic velocity field with oscillating flow rate Q = 0.1 cos(0.5t)mm 3 /s is prescribed.This simulation is performed for the case of only a surface tension of γ = 50 μN/mm, the reference state as described in "Viscoelastic interface", the reference state with κ = 100 μNs/mm as the only changed parameter, and the  22, where V is the volume of the drop, L apex is the apex length, p is the pressure at the inlet in , and R is the radius of curvature at the apex.The subscript eq indicates the value of the variables at t = 1 s.
Figure 22 shows that the responses of the apex length, the pressure, and the radius of curvature to an oscillating volume are nearly sinusoidal.For the cases (a), (b), and (d), the pressure is out of phase compared to the radius by approximately π radians.This is in accordance with the findings of Kotula and Anna (2015): for a purely elastic interface limit of κ = 0 μNs/mm, the phase shift goes to π radians.Furthermore, from Kotula and Anna (2015) follows that a phase shift of π/2 is the limiting case when the dilatational viscosity is much greater than the surface tension.In case (c), where the dilatational viscosity κ = 100 μNs/mm is higher than the surface tension γ = 50 μN/mm, the phase shift tends towards π/2 radians.

Conclusions
In this study, a two-dimensional axisymmetric model is developed to simulate pendant drop experiments using a finite element method.In this model, the dynamics of the interface are described by a Lagrange-based interface tracking model.Furthermore, the quality of the mesh is monitored, and mesh refinement is performed if necessary.A semi-implicit time integration scheme is used to calculate the surface tension to cope with the small capil- A quasi-linear two-dimensional Kelvin-Voigt viscoelastic constitutive model is implemented at the interface, following Verwijlen et al. (2014).The model is made quasi-linear by using a combination of two finite strain tensors; the Hencky strain tensor is used for the dilatational contribution and the Green-Lagrange strain tensor is used to describe the shear contribution.A computational tool is designed to calculate the tangential and circumferential surface strain to quantify the deformation of the interface.
A parameter study is performed to investigate the influence of the five parameters of the quasi-linear twodimensional Kelvin-Voigt model: the surface tension γ , the surface dilatational viscosity κ, the surface shear viscosity μ, the surface dilatational elasticity K, and surface shear elasticity G.During inflation, the apex length of the drop is significantly smaller for higher values of the surface tension and the dilatational elasticity, and is slightly smaller for higher values of the surface dilatational viscosity.Unless the force balance at the outlet of the capillary initiates the onset of detachment, the shape of the drop stayed at an equilibrium state during relaxation for different values of the surface tension and dilatational elasticity.For all different values of the dilatational viscosity, the shape of the drop became equal after total relaxation.For higher values of the surface tension, the entire volume of the drop flowed back into the capillary faster.The decrease of the apex length during the outflow regime was slower for higher values of the dilatational viscosity.For values κ ≥ 50 μNs/mm, the apex length increased at first due to a flow in negative radial direction in the region of the drop just outside the capillary and the resistance to deformation of the interface.The surface tension seemed to dominate the outflow behavior for values of the dilatational elasticity K ≤ 10 μN/mm.For higher values of the dilatational elasticity, an equilibrium between the elasticity and the surface tension is reached during the outflow regime.The influence of both the shear viscosity and the shear elasticity on the apex length is negligible for the used parameter range.
For elevated values of the surface viscosities and elasticities, the higher values of the surface shear viscosity delay the relaxation and prevent the drop to "collapse", while increasing surface shear elasticities give rise to a more elongated shape during both the inflow and outflow regime.
A small amplitude oscillatory measurement is simulated.The phase shift between the response of the radius of curvature and the pressure for more elastic interfaces is approximately π radians, while for dominant dilatational viscosity the phase shift tends towards π/2 radians.
The simulations should be compared with experimental results.A relevant example of surface active agent with viscoelastic behavior is the protein hydrophobin (HFBII) (Alexandrov et al. 2012;Danov et al. 2015;Knoche et al. 2013); it forms a film of high mechanical strength when adsorbed to the interface.Furthermore, the geometry from the numerical computations could be used as input for the experimental software, to compare the input parameters of the model with the solutions from the experimental setup.

Fig. 1
Fig. 1 Geometry t=0 of the pendant drop, where the outer radius R o = 0.35 mm, the inner radius R i = 0.225 mm, and the length of the capillary L = 0.8 mm

Fig. 2
Fig. 2 Mesh M1 of the mesh-convergence study, where line is given in red

Fig. 3 Fig. 5
Fig. 3 Mesh-convergence of u (a) and p (b) Fig. 1, and equal boundary conditions and parameters (as stated in the problem description), the same problems are solved to compare the results of the two computational methods.The only interface characteristic is the constant surface tension γ .Three different flow rates are examined: Q = 10 mm 3 /s, Q = 25 mm 3 /s, and Q = 42 mm 3 /s.The corresponding Reynolds numbers are Re = 28.29,Re = 70.73,and Re = 118.84,respectively.To compare the data of the work of Dieter-Kissling et al. (2014) with this study, the drop formation for different flow rates is investigated.The apex length is plotted versus the drop volume, as shown in Fig. 4.Both the experimental data and the data generated by the simulation of Dieter-Kissling et al. (2014) are displayed along the results of this study.

Fig. 8
Fig. 8 Necking radius in time for Q = 10 mm 3 /s, Q = 25 mm 3 /s, and Q = 42 mm 3 /s, respectively Fig. 11Velocity in r-direction for the reference state system at time t = 0.01 s, t = 0.50 s, t = 1.00 s, t = 1.01 s, t = 1.04 s, t = 1.05 s, t = 2.00 s, and t = 2.01 s, the velocity is in mm/s

Fig. 12
Fig. 12Pressure for the reference state system at time t = 0.01 s, t = 0.50 s, t = 1.00 s, t = 1.01 s, t = 1.04 s, t = 1.05 s, t = 2.00 s, and t = 2.01 s, the pressure is in μN/mm 2 Fig. 13 Apex length in time for different values of the surface tension γ

Fig. 14 Fig
Fig. 14 Strain distribution for the tangential surface strain and the circumferential surface strain for surface tension γ = 10 μN/mm (left) and γ = 50 μN/mm (right) at time t = 1 s

Fig. 16
Fig. 16 Apex length in time for different values of the surface shear viscosity μ Fig. 17Flow field for κ = 100 μNs/mm at time t = 2.01 s (a), a close up of the region at the end of the capillary (b), the velocity in z-(c) and r-direction (d), and the pressure (e)

Fig. 19
Fig. 19 Apex length in time for different values of the surface shear elasticity G

Fig. 20
Fig. 20 Apex length in time for different values of the surface shear viscosity μ at a surface dilatational viscosity κ = 100 μNs/mm

Fig. 21
Fig. 21 Apex length in time for different values of the surface shear elasticity G at a surface dilatational elasticity K = 100 μN/mm

Fig. 22
Fig. 22Simulated responses of the normalized volume V /V eq , the apex length L apex /L apex,eq , pressure p/p eq , and radius of curvature R/R eq , due to an oscillating flow rate Q with amplitude a = 0.1 mm 3 and angular frequency ω = 0.5 rad/s, for only a surface tension of lary time scale of the problem.The model is verified with the results as presented in a publication by Dieter-Kissling et al. (2014) for a Newtonian fluid in a Navier-Stokes flow, with gravitational forces and constant surface tension.

Table 1
Refined meshes used in the mesh-convergence study

Table 2
Parameter range where the filled values denote the reference values formation of this work matches very well with both the experimental results and computational data of Dieter-Kissling et al. (2014).Our results even exceed the range of deformation of their numerical simulation. drop

Table 3
Maximum surface strain for different values of the surface tension γ at the end of the inflow regime and the relaxation regime