In this tutorial we present an example of surface-transport equations in a free-surface Navier–Stokes problem. This is a multi-domain, multi-physics problem because there is a coupling between the equations on the surface and those in the bulk. The coupling from the bulk to the surface transport arises through the surface velocity in the surface transport equation. The coupling from the surface transport to bulk arises in a more subtle manner because the surface concentration affects the surface tension. In this tutorial, we describe how to use the existing framework to create surface transport equations and how to include them in a free-surface problem.
The example problem
The problem to be solved is the evolution of an annular film of fluid on the inside of a solid cylinder in the presence of an insoluble surfactant on the interface: a modification of the classic Rayleigh–Plateau instability. For validation, we reproduce some of the results given in `A 2-D model of Rayleigh instability in capillary tubes — surfactant effects' by D. Campana, J. Di Paolo & F. A. Saita, Int. J. Multiphase Flow, vol 30 , pp 431–454, (2004). Our formulation, however, is different from their approach as described in detail in our free-surface theory document.
| The unsteady axisymmetric free-surface Navier–Stokes equations with insoluble surfactant . |
in the bulk fluid.
The governing equations are subject to the no slip boundary conditions
on the outer solid boundary ( ) and the symmetry boundary conditions
on the bottom ( ) and top ( ) boundaries.
We denote the position vector to the free surface by , which is subject to the kinematic condition
and the dynamic condition
where is the dimensionless surface tension relative to a reference value.
An insoluble surfactant of surface concentration is non-dimensionalised with respect to a reference value and obeys the surface transport equation
on the interface.
The surface tension is a function of the surfactant concentration and a linear equation of state is chosen
The symmetry boundary conditions on the bottom ( ) and top ( ) boundaries are
Initially, the system is at rest and . The free surface is moved into the position:
where is a small parameter and is the undeformed film thickness.
We choose parameters based on those used to compute Figures 8 and 9 in Campana et al; namely , , , , , , and . For these parameters, the system is unstable to the Rayleigh–Plateau instability and evolves towards a state in which the tube is completely occluded by the fluid at one end.
Time trace of the radius of the interface at z=0 showing dramatic collapse near t=80.
Evolution of the interface at times t=0, 10, 20, 30, 40, 50, 60, 70, 80 showing the developing lobe.
Evolution of the surfactant concentration on the interface at times t=0, 10, 20, 30, 40, 50, 60, 70, 80 accumulation in the developing lobe caused by the reduced surface area and advective flow into the lobe.
Global parameters and functions
The global parameters are simply the dimensionless parameters described above.
The driver code and problem class
The driver code and problem are very similar to those in the two-dimensional and axisymmetric interface-relaxation problems on which this driver was based. The main difference between this problem and standard free surface problems is that instead of oomph::SpineAxisymmetricFluidInterfaceElement, we use the custom oomph::SpineAxisymmetricMarangoniSurfactantFluidInterfaceElement. The symmetry boundary conditions for the surface concentration are the natural boundary conditions of our formulation, so we "do
nothing" for the additional field at the boundaries.
An additional member function
InterfaceProblem::compute_total_mass() is provided as a check on the implementation of the the surface transport equations. The surfactant cannot be removed from the surface, so its mass must be conserved. The function simply loops over the interface elements and sums their contribution to the total mass.
double mass = 0.0;
const unsigned n_interface_element = Interface_mesh_pt->nelement();
SpineAxisymmetricMarangoniSurfactantFluidInterfaceElement<ELEMENT>* el_pt =
mass += el_pt->integrate_c();
The SpineAxisymmetricMarangoniSurfactantFluidInterfaceElement class
This class is implemented in our driver code and inherits directly from oomph::SpineAxisymmetricFluidInterfaceElement. The class provides storage for the required additional dimensionless groups and the nodal index where the surface concentration will be stored.
class SpineAxisymmetricMarangoniSurfactantFluidInterfaceElement :
Most of the functionality is already provided by the underlying
FluidInterfaceElement and we need simply to overload a few functions. The constructor sets default values for the physical constants and adds the additional data value to nodes on the surface.
SpineAxisymmetricMarangoniSurfactantFluidInterfaceElement(FiniteElement* const &element_pt, const int &face_index) : SpineAxisymmetricFluidInterfaceElement<ELEMENT>
Beta_pt = &Default_Physical_Constant_Value;
Peclet_S_pt = &Default_Physical_Constant_Value;
Peclet_Strouhal_S_pt = &Default_Physical_Constant_Value;
unsigned n_node_face = this->nnode();
for(unsigned i=0;i<n_node_face;i++) additional_data_values[i] = 1;
C_index = this->node_pt(0)->nvalue()-1;
FluidInterfaceElement::sigma() is overloaded using the equation of state defined in the problem specification
double sigma(const Vector<double> &s)
const unsigned n_node = this->nnode();
C += this->nodal_value(l,C_index)*psi(l);
return (1.0 - Beta*(C-1.0));
The majority of the work is performed in
void add_additional_residual_contributions_interface(Vector<double> &residuals, DenseMatrix<double> &jacobian,
const unsigned &flag,const Shape &psif, const DShape &dpsifds,
const DShape &dpsifdS, const DShape &dpsifdS_div,
const Vector<double> &s,
const Vector<double> &interpolated_x, const Vector<double> &interpolated_n,
const double &W,const double &J)
which provides the additional surface transport equations. In the example code two formulations of the surface transport equations are provided the one used by Campana et al in which the curvature is computed explicitly and the formulation derived in our free-surface theory,
in which the curvature is not required. The version used is determined by an internal boolean
bool Integrated_curvature = true;
The remainder of the function adds the residuals associated with the surfactant transport equations which are described in the free-surface theory.
Note that an additional term arises due to the azimuthal curvature compared to the standard one-dimensional surface.
fill_in_contribution_to_jacobian(...) is also overloaded to that the effect of surfactant concentration on the bulk equations is computed by finite differences. This could be modified in the future so that the appropriate derivative terms are included in
Finally the elements contain a function
double integrate_c() const
unsigned n_node = this->nnode();
unsigned n_intpt = this->integral_pt()->nweight();
double mass = 0.0;
s = this->integral_pt()->knot(ipt,0);
double W = this->integral_pt()->weight(ipt);
double interpolated_c = 0.0;
interpolated_c += this->nodal_value(l,C_index)*psif(l);
interpolated_x[i] += this->nodal_position(l,i)*psif(l);
interpolated_t[i] += this->nodal_position(l,i)*dpsifds(l,0);
double r = interpolated_x;
double tlength = interpolated_t*interpolated_t +
double J = sqrt(tlength);
mass += interpolated_c*r*W*J;
double SpineAxisymmetricMarangoniSurfactantFluidInterfaceElement<ELEMENT>::Default_Physical_Constant_Value = 1.0;
template <class ELEMENT>
class MyHorizontalSingleLayerSpineMesh :
that computes the integral of the concentration over the elemental surface, representing the total mass of surfactant within the element.
- Investigate the difference between the solutions for the two formulations of the surfactant transport equations. Which conserves mass more accurately?
- Investigate the influence of variations in and try to reproduce the results found by Campana et al.
- Look at the three-dimensional (non-axisymmetric) version of the code found in the same directory. Confirm that the same results are produced. Is the instability stable to non-axisymmetric perturbations? Use the code to investigate what happens if you make the cross-sectional boundary slightly elliptical rather than circular?
Source files for this tutorial
A pdf version of this document is available.