Superficial transportation model using finite volume method

The work presents a new modeling technique designed for transporting selected property while simultaneously tracking particular value of certain field variable. The model addresses the problem of transportation of superficial properties according to fluid/fluid interface evolution. Such problem arises when some superficial variable needs to be attached to the evolving interface or front defined. Several multiphase models in CFD technique lack such ability. The model proposed is illustrated by the experimental two phase system of toluene/water with surfactant specie (sodium dodecylsulfate SDS). The selected superficial property that is important to be tracked according to interface evolution is Gibbs surface excess. The computation was performed using Ansys/Fluent software with proposed models created and attached using C programming language.


List of symbols
Weber number Greek letters ø General scalar quantity α Phase volume fraction η Viscosity (Pa s)

Introduction
Interfacial phenomena are the subject of investigation in many branches of the industry. The phenomena occurring on the interface are analyzed in this work by the model proposed using experimental conditions which are good test for proposed modeling method. The evolving and aging interface during the process of drop formation is connected with the effects of changing interface curvature along with its physical properties from which the most important are surface excess and interfacial tension. In addition, the process of drop formation which is analyzed in this work has a great significance in several fields. The typical areas are the purification extraction processes [1][2][3], metallurgic processes of gas metal arc welding [4], emulsification processes [5][6][7], jet print processes [8][9][10], the oxygenation of blood equipment design [11,12] or pool boiling bubbles spreading [13] and manufacturing liquid microcapsules by centrifugal extrusion or spray drying [14]. Besides the pendant drop, the CFD simulation of liquid drops in different conditions is treated in literature, e.g., drop flow, impact, spreading, sliding, solidification on solid surfaces [15][16][17][18][19][20][21], the process of a liquid drop formation in the coaxial capillars [22], the effects of local interfacial nonhomogeneity due to local surface active specie concentration differences [23]. The problems of liquid drop oscillations is analyzed in other works of the authors [24,25]. The analysis of water-octane system [26] using commercial software FLUENT in 3D domain showed the model deviation is below 5% of relative error in a drop life time and its volume. Because the problem of multiphase flows and accompanying superficial effect are still being investigated, expanding the modeling capabilities of existing models of fluid/fluid interface has not only scientific and explorative significance but also may extend the capabilities of a designer's tools. An appropriate description of the surface processes allows to achieve more accurate predictions and consequently leads to a better understanding and utilization of the surface phenomena. The drop shape measurement is a good test for interface models because any model discrepancy can easily be detected in the condition of growing and aging a pendant drop. Moreover, it provides the opportunity to capture interface phenomena with high accuracy. Several techniques in laboratory practice, e.g., the capillary rise, Wilhelmy plate method, Du Noüy ring method, drop weight method, pendant/sessile drop and bubble methods, spinning drop method, or the maximum bubble pressure method are all derived for interface examination and research.
This work presents the extension to computational fluid dynamics volume of fluid model which allow taking into account the surface phenomena on the fluid/fluid interface. The key innovation proposed is the algorithm of superficial property transport along with the interface dynamic evolution. The original, typical multiphase VOF transient model allows only to describe interface only at current time step by its primary property-surface tension. It is not possible to use the typical VOF approach when calculation of dynamic evolution of interface requires that previous time step and previous interface location to be known for given property calculation. In other words, it is not possible to attach any field variable to specified interface location locally and translocate it together with its movement. At given time step, the VOF model is not able to locate previous time step locations, and previous time step values of any interface variable. The previous time step interface position in space and coupled with this location superficial variable must be known to describe some processes properly, for example in the case of interfacial adsorption of surfactants. One of the fundamental thermodynamic variables for the cases when surfactants are involved, is the Gibbs surface excess Γ , a quantity describing superficial concentration of adsorbed, surface active specie [27]. The dynamic simulations in the CFD technique are based on the time stepping algorithms which are applied by various numerical schemes. The calculation of surface excess Γ for current location on the liquid interface and for current time step can only be done referring to previous time step and previous location values. Due to this requirement and taking into account that well-known VOF algorithm contains no methods nor techniques for attaching any physical property to the interface, the proposed model is a useful extension allowing to simulate the surface excess dynamic advancement.
The Eulerian frame of reference of the VOF model reasons the need of interface reconstruction on the basis of specific color function. A locally defined surface tension, Gibbs surface excess and interface mass flux vector need to evolve according to the interface evolution. Proposed approach introduces lagrangian reference frame to the description of the interface by the use of markers indicating the surface position and time evolution. The set of markers is then translated according to the velocity field and phases volume fractions field. The markers might be thought of as a lagrangian virtual particles that translocate with the fluid elements. In the present work, they are independent from the underlying grid and are used to carry local information about the interface. The marker and cell (MAC) method is known in literature and was investigated to keep track of the front shock waves (shock-capturing methods) or interfaces between fluid phases. There are algorithms which rely on the use of particles to track desired variable or flow property. Some of the examples are the MAC method [28] and its several variants: SMAC [29] (simplified marker and cell), GENSMAC [30] (general domains simplified marker and cell), GSMAC [31] (generalized simplified marker and cell).
Similar to the presented works above, there are hybrid methods that combine the VOF and MAC techniques to achieve specific aims. Aulisa [32] utilizes concept of conservation markers to construct the interface. The algorithm is focused on correcting the interface calculated by standard VOF method to preserve its area. Their other works [33,34] focus on algorithm for the reconstruction and advection of interfaces using marker and cell grid but limited to the two-dimensional space. Enright [35] puts two sets of markers called "positive" and "negative" at either side of the interface to rebuild the region for specific cases that otherwise stay underresolved giving unrealistic geometric interface topologies.Ďurikovič [36] presents technique of multifluid (more than two phases) VOF simulation discretized on MAC grid. The technique is used mainly for visualization and animation purposes. The most similar to presented work use of markers with VOF model is the corrective algorithm proposed by Lopez [37]. The markers are placed at the midpoints of the reconstructed cell interface segments and are used to resolve the flow characteristic smaller than grid size. Typical usage of this approach is to account for the small fluid filaments creation phenomena.
The model proposed makes use of the Lagrangian method of fluid dynamics [38]. This technique utilizes an collection of virtual Lagrangian particles of zero spatial size whose paths in space are governed by the velocity field. The velocity field which are utilized to move the particles originates from CFD core algorithm implemented in Ansys/Fluent. The tracks of virtual particles map the streamlines of the velocity field, corrected at every numerical time step by the location of selected isovalue of volume fraction. By following the flow of virtual particles, and assigning a specified superficial variable (and possibly other properties) to them during calculations, the problem of interfacial properties translocation and their time evolution is addressed and solved.
Lagrangian approaches are extensively utilized in engineering, including discrete element methods [39] and smoothed particle hydrodynamics [40]. The Lagrangian structure for studying pathways corresponds to the analysis of tracers. One of the significant dissimilarities is the computational cost. For each time step, translocation of a Lagrangian particle takes only one calculation route of relatively simple computations. On the other hand, considering the advection-diffusion of a selected tracer variable, it takes several steps of computation to estimate its trajectory.

Modeling the interface
The concept of Gibbs dividing plane [41] is essential when defining the surface excess Γ of surfactant. The importance of the shape of the interface was realized and analyzed by Gibbs [42], Guggenheim [43], Gurtin [44], Oversteegen [45] and others. The position of dividing plane impacts the perceived surface excess value, and its location is chosen to conform with its zero value for clean solvents [46]. Most of computational fluid dynamics multiphase models present the interface as a surface geometrical manifold of an infinitesimal thickness which is located in the area of sudden change of physical and chemical parameters. Presence of such interface in the simulated case is a good testing environment for the algorithm for superficial property time evolution.
The formulation of the CFD technique involves the classical Navier-Stokes set of equations, which originate from Newton's Law of Motion: The mass conservation (or continuity equation) is used for total mass balance: The dynamic evolution of the interface as a result of the force and, consequently, the velocity field influence can be simulated using several approaches in CFD. The front tracking methods [47], level set methods [48] and volume of fluid methods [49] are the most commonly used ones for interface modeling. For the purpose of interface capturing Eq. (3), which specifies the evolution of the volumetric amount of the phase contained by the cell, is introduced in addition to the above. This quantity is the volume fraction α of a given phase, and the sum of all phase volume fractions in a given volume must equal one.
The explicit spatial discretization of geometric reconstruction was applied to the equation (so called Geo-Reconstruct). It is the most accurate and precise scheme available in Fluent, which unfortunately is more computationally demanding than the other schemes. Geo-Reconstruct is the best scheme when using lowquality meshes. The time advancement scheme chosen for this work is an explicit formulation. It is noniterative for given time step and is time-dependent scheme. It reveals better numerical accuracy compared to the implicit formulation, which is burdened with negative numerical diffusion phenomena. Though the time step size might be severely limited by a Courant number, typical stability criterion used in VOF modeling. Solution to Eq. (3) contains a discontinuity which is expressed as sudden jump of computed quantity α in the location of the interface. Such nondiffusive multiphase continuum as a limiting case of differential diffusion transport equation is studied by Mellado et al. [50]. Several algorithms are utilized to reconstruct the interface, having solved the distribution of α in the space [51][52][53][54][55][56].
The force source F generated by the interfacial tension at the interface is important quantity during interface reconstruction. Brackbill et al. [57] proposed the continuum surface force method which is based on the theory that the pressure jump across the interface is proportional to the curvature of the interface κ and the interfacial tension, according to the Young-Laplace equation: The wall adhesion and, consequently, contact angle θ play also significant role in determining the fluid behavior in the location of the fluid-fluid-wall contact line. The contact angle is accounted for by the use of formula (5) for unit vectorn normal to the interface at the wall contact point.
n =n w cos θ w +n t sin θ w (5) n t is a normal vector to the contact line between the interface and the wall at the contact point andn w is the unit wall normal directed to the wall. Implementation of the contact angle into fluid dynamics simulations may pose some serious difficulties as presented by Linder et al. [58]. In this work, the two immiscible liquids system is experimentally analyzed and used to make verification of the model proposed. In the system examined, the liquid droplet of heavy phase is formed into the lighter liquid phase. Several studies are devoted to simulation of droplet formation and its behavior [59][60][61]. In this work, the experimentally examined toluene-water system with sodium dodecyl sulfate (SDS) as a surface active specie was a base to validate the superficial transportation model proposed. The fundamental motivation for studying interface surfactants activity at liquid-liquid interface is to describe the process of liquid-liquid phase transfer catalysis by means of CFD. This process is a key element in interface chemical reactions in many industrial processes that involve green methodology [62,63].

Bulk phase mass transfer
The description of component mass transfer in model proposed is divided into two different mechanisms and consequently two distinct calculation approaches. The bulk phase and near-interface mass transfer differ due to different nature of driving forces. The calculations performed in this work utilize the Fickian approach for bulk diffusion [64][65][66]. The bulk phase binary diffusion coefficient of sodium dodecyl sulfate through water was calculated using Siddiqi-Lucas correlation which was derived based on over 600 different water solutions [67]: The above correlation is based on metric units and expresses diffusion coefficient in cm 2 /s. The value of SDS molecular volume V i was taken from literature [68]. The component flux is then calculated using Fickian approach: In the case of multiphase flow, the equation of fluid motion (1) and continuity (2) are additionally accompanied by the species mass balance equation for component i and given phase j: The equation above is solved for spatial locations where the jth phase exists. In location of interface, the weight of the transport is minimized by the value of volume fractions α j . In that case, near the interface, adsorptive mass transfer is main transportation mechanism.

The near-interface adsorptive flux computation
The key part of model proposed, where the markers introduced are involved, is the way the mass flux due to adsorption is calculated. The near-interface location is understood in model proposed as a set of cells containing the calculated, reconstructed interface. Unlike the bulk mass transfer, the mass flux near the interface is induced mainly by adsorption of surface active component. Applying appropriate mass transfer approach, expressed in modification either in diffusion coefficients or model equations themselves, is typical to describe the mass transfer in any media, where concentration gradients do not generate dominant driving force (e.g., ions flux in presence of electric field) or in cases where molecular diffusion coefficients are not known and are substituted by effective diffusivities [69][70][71]. Calculating adsorptive near-interface mass transfer creates the need for effective diffusion coefficient estimation. First approach to dynamics of the interface adsorption was presented by Ward and Tordai [72] who formulated first quantitative model of diffusion controlled adsorption, then the mixed adsorption kinetic model proposed by Baret [73] and, further, short and long time approximations proposed by Feinermann [74]. The thermodynamic description of system containing interface in the form of free energy functional is widely discussed [75][76][77][78][79] and also the detailed presentation of the models used to describe the adsorption processes at the liquid-liquid interfaces can be found in literature [80]. Such formulations give the adsorption dynamics description in a form similar to original Ward and Tordai equation but without assumption that initial surface excess is equal zero [81]. The significant part of the model proposed in this work is the ability of computation surface active components interface mass flux without the use of diffusion coefficients. In order to compute the interfacial mass flux generated by adsorption forces, it is formulated according to specific interface description. Previously mentioned Fickian relationship for bulk phase relates diffusion coefficient with concentration gradient. In the vicinity of the interface, the adsorptive diffusion flux is given by: Which is then modified by applying Gibbs surface excess to give known relationship: The above formulation introduces the surface excess as a time-dependent function Γ (t). The mass flux that directs into the interface is given by the change of surface excess in corresponding time period so to be most accurate a proper definition of such a function is later in the text proposed. The values of Γ (t) at consecutive time steps must be determined in order to estimate the derivative dΓ (t)/dt. Such function can easily be estimated with the use of Szyszkowski [82] isotherm coefficients although other similar isotherms (like Frumkin) could also be applied. To make distinction between functional forms and experimental data of the same quantities in the text, the subscript exp is added wherever experimental ones are used. The experimental values of interfacial tension γ exp are related with corresponding bulk concentrations c exp by Szyszkowski isotherm: Based on the experimental measurements of the interfacial tension γ Sz exp for different bulk concentrations c exp , the values of A sz and B sz coefficients are determined. This is done using least-squares method so the Szyszkowski curve fits the experimental data [83]. Figure 1 presents typical data acquired for SDS and corresponding, fitted using Szyszkowski isotherm. The experimental points are obtained on the basis of 600 continuous measurements, and corresponding statistical error is so small (between 10 −4 and 10 −5 N/m) that it is not presented on the graph.
The experimental values of a surface excess are estimated from the equation: The bulk concentrations are calculated by the rearrangement of Szyszkowski equation:

Fig. 2 Dynamic interface tension measurements
The dynamic surface tension γ dyn exp (determined experimentally) values are used to determine coefficients of the function Γ (t), which is formulated using arbitrary form suitable in the sense of its monotonicity.
The selected function and its derivative are continuous which is important for further use. The function is typical black-box approach and does not explain the process of adsorption in the sense of transport variables like diffusion coefficients. The choice of this exponential function form was made because it accurately renders the time evolution of Gibbs surface excess and due to its mathematical simplicity and robustness. Very important property of the choice made above is also the function differentiability, which is required for further steps and can be done on symbolic, analytical basis. The constants s in the above formulation are determined for each dynamic interfacial tension measurement based on the least-squares technique (Fig. 2). Finally, the function J ads (c) is estimated based on Eqs. (10) and (14) for several experimental values of bulk concentrations c exp . This allows the adsorptive flux J ads to be related to the bulk concentration c. The knowledge of the adsorptive flux together with the interface area is essential to mass transfer computation in the model proposed. In many commercial tools like Fluent, Tecplot or CFD-Post, for every volumetric cell, having known the values of phase volume fractions, the interface shape is reconstructed using marching cubes method [84]. This algorithm is widely used computational geometry algorithm for construction of isosurfaces at given evaluated location in the space, but it is not used in the present work. Instead own algorithm was used which is described in the next section.

The surface variables propagation algorithm
The proposed markers algorithm works in a few following computational steps.

Markers preparation
Before the iterative process of integrating the Navier-Stokes equations can be started, a computational domain must be initialized. During the initialization stage, when boundary and internal domain conditions are applied, the markers are placed according to a phase locations. According to Eq. 3, when a phase occupies whole control volume, its volume fraction α equals 1. An interface (geometric isosurface) is located at position for which α equals 0.5 (designated by α iso ) and so the initial markers are located. The procedure is straightforward, having values of α at cell nodes, the position vector x of all markers at each cell edge is determined according to ith node distance weighted mean of α i :

Fig. 3 Construction of initial markers
The locations x of edge markers located at α iso are then averaged to give initial markers locations in the cell interior. The situation when α iso equals α i means that the interface is located exactly at specified node of the cell. This condition, if happens, is handled by the code by if-then clause. Figure 3 presents the edge markers created using Eq. (15) and initial markers • used for further computations. The edge markers are discarded after initial ones are created. Such preparation technique ensures that the initial markers are located at a cell interface superficial elements centroids. This does not prevent them from clustering in the case when interface is located near one of the cell edges. This algorithm does not address such behavior as it has no impact on surface variables propagation and most important, it ensures that every cell in which interface is located has its own marker. At the final stage, for every cell containing the marker, the adsorption model, presented before, is computed to give superficial values of component surface excess. This value is stored at marker for further processing.

Markers projection
The algorithm of markers relocation is designed to follow the stepwise approach of interface transient evolution. Based on the original/initial locations of markers evaluated at previous time step and current time step interface location, the advancement procedure is performed. The superficial variable specified for previous time step and recorded at each marker has to be moved to new location of the interface. For this reason, a new spatial coordinates of the markers are calculated to locate actual, new position of interface. Each marker has its corresponding, individual velocity vector specified, and such a vector indicates the direction of advancement and relates previous to actual superficial variable position. However, the momentum equation is most often solved using different scheme of discretization when compared to the discretization of mass balance for given phase. For this reason, the new coordinates found by velocity vector might not exactly be the same as the new position of the interface. On the range of one time step, this discrepancy is very small but may increase during successive and repeated time steps. At this stage, the new set of markers (second, auxiliary set) is created at the exact location of the interface using the same method as described in the previous section. The two markers sets are kept in memory in order to copy the values of superficial variable from first markers relocated on the basis of velocity vectors to the new markers created on the actual interface. At this moment, it is important to associate markers from the first set with the appropriate new markers created. Such logic link between each marker from first set and each actual interfacial marker is required to perform the advancement of superficial variable in correct and appropriate direction. The link between the initial marker and actual marker is created on the basis of distance matrix using nearest neighbor algorithm. The two sets of markers are in close vicinity, typically in the distance much smaller than size of a single cell, but this is not a requirement because markers are not connected to mesh in any way. After the values of superficial variable are copied to the second markers, the first markers are discarded, the memory is freed, and the second set becomes actual set of markers, located at new interface position, containing the translocated value.

Transported quantity
This is important to note that the algorithm presented is suitable only to transport an intensive quantity due to balance considerations for computational mesh cells which have different volumes. The component surface excess is intensive quantity which is updated at every time step by the adsorption model described earlier in the text. Having computed the flux J ads of surface active specie into the interface, the value of surface excess is updated accordingly by the time advancement formula: In the case when the interfacial mass transfer is involved, the corresponding number of component moles are removed from the local cell containing the interface. The ith specie mass source (sink) is defined: and expresses the mass that is transported to the interface area S int due to adsorptive flux in the cell volume V cell . The m i is the source term included into the species mass balance equation Eq. (8).
The Ansys/FLUENT software allows to link own model, being an extension to existing CFD code, at different stages of a solution process. The attachment of the model procedures is allowed in a variety of computational states, for example for already available phenomena models, numerical states or geometric locations and in many other situations. The code that is written in C programming language is compiled in Ansys environment and then attached to the CFD computational core, for example at geometric boundaries as boundary conditions, in fluid volume as source terms, or in some specific point of computation like at the end of a numerical iteration. Such user-defined functions (UDF) were used in this work to formulate the described model and is attached to execute at the end of every numerical iteration (execute_at_end UDF function). The procedure loops over every mesh cell in a calculation domain and locates the cells containing the interface. For those cells, the calculations of the algorithm of the adsorptive mass transfer model are performed. The schematics of the procedure proposed are presented in Fig. 4.
The markers are programming objects defined as structured type that aggregates a fixed set of different types objects into one single object. The marker object structure contains: its location defined as three-dimensional

Superficial mass balance and numerical diffusion
A test simulation to perform a verification of the superficial mass conservation and numerical diffusion was prepared. In a cuboid domain (0.1 × 0.1 × 0.1 m), a two phase air-water system was defined and slow upward motion was applied by the inlet boundary located at the bottom. The vertical fluid movement was introduced to exert a momentum forces on the surface and to observe the superficial variables evolution under dynamic conditions. Actually, the surface excess Γ is the superficial, transported variable in the model proposed but for the test a unity value (g = 1) was applied instead. The test values g = 1 were distributed onto half of the interface S 1 (see Fig. 5) while setting the rest of the space with zero value (g = 0). A fluid velocity 0.001 m/s and 10 s of simulation time moved the surface upward enough to allow observation of superficial distribution of values applied initially. The algorithm proposed is not supposed to apply any surface diffusion component, but artificial diffusional effect is possible. To estimate the quantitative numerical effect of artificial diffusion component, the surface integral mean over the interface of transported testing values was calculated at every time iteration: The integral (18) defines total value of g over the surface part S i located in given part of the volume (Fig. 5). Any change of the averaged integral value G i indicates artificial diffusion effect between the volumes. The total sum of G i for both parts of the interface should result a unity value; any deviation from this value will indicate mass imbalance. Actual results for the test performed are presented on graph in Fig. 6 showing surface average integrals of test value g i . The curves represent the actual, recalculated for percentage scale, surface average of value G i . The ideal and expected values (for the case of no artificial diffusion) should obtain 100% for S1 (absolute value of G 1 = 0.005 m −2 due to the size of the domain) and 0% for S 2 (absolute value of G 2 = 0.0 m −2 ). Any deviation from those two values shows the numerical imbalance which is caused by artificial diffusion. Figure 6 shows typical character of any finite element method applied to given problem. The higher is the number of finite elements, which in the case of the CFD technique are finite volumes, the higher is the accuracy of results obtained. The results of the test performed (Fig. 6) show that the numerical imbalance and artificial diffusion, which is measured as deviation from value of 100% for S 1 and 0% for S 2 , is smaller for higher number of markers. Higher loss of mass is observed for number less than 4000 markers. The artificial Fig. 6 The surface average integral of test value g for both part of the surface analyzed: S1 and S2, used in the artificial surface diffusion test diffusion is visible for the parts where S 1 curve minimizes with corresponding S 2 curve increase (from 2000 to 4000 of markers). This indicates mass being transported unintentionally from the S 1 to the S 2 area. The test illustrates that such negative numerical effects of the algorithm proposed can be minimized by using higher number of finite elements during the calculation.
The visible error of the balance is mainly related to the VOF model used in this work. The markers follow the interface in accurate and exact way, tracking its position in every time step. Thus, the main source of the imbalance lays in the VOF numerical instability, namely spurious currents which affects the interface stability. This negative effect is caused by the VOF model assumption to use single, common velocity field approach. This phenomena is known as disadvantage of VOF method which becomes more visible for low velocities. In the example presented, it is deliberate to render the situation which is difficult from the modeling point of view. The low velocity of the interface transition which was applied in the case presented is typical setup for which such a model can be tested by drop volume method. There is possibility to make a choice of any other multiphase approach, where the limitation of common velocity field is not imposed. But because of the popularity of VOF model, this work is focused on the possibility to extend this approach with all its benefits and disadvantages.

Experimental part
A system of two immiscible liquids (toluene and water + SDS) was chosen for the verification of the CFD simulation. The Tracker tensiometer was used to perform the drop shape experiments. Seven samples at different SDS concentrations were submitted into the tests to build Szyszkowski isotherm. The twice distilled water was used as the aqueous and twice distilled toluene as organic phase. Samples were shaken for 4 h to assure equilibrium state at the interface. The measurements were carried out at 25 • C. Measurements are in general accordance to other works [85,86]. The fit goodness, measured as a square error, for the isotherm equals 4.367 × 10 −6 which gives about 2 mN/m isotherm mean deviation from the experimental values. Taking into account that determination coefficient for that isotherm is at value 0.963 which means that the isotherm explains 96.3% of all measured data, such accuracy assures that the Szyszkowski isotherm is a good description of the system.
The A and B coefficients of Szyszkowski isotherm were found to be 7.736 × 10 −4 and 0.069 accordingly. Following previously described procedure the SDS mass fluxes were computed for different bulk concentrations ( Table 1).
The values of flux were then related to bulk concentrations by the use of modified Weibull [87] function (19) by means of least-squares approximation.   The values of adjusted coefficient used in flux J calculation are given Table 2.
It was mentioned earlier that adhesion forces at the domain boundaries have major impact on the interface formation process. Incorporating the contact angle between the fluids and a solid wall was done by defining the contact angle θ at the wall by formula (5). Experimental data shows tendency of exponential growth of contact angle by increasing surface tension (Fig. 7).
The contact angle was related with surface tension by exponential function fitted to the obtained experimental values: The calculated by the least-squares method values for the equation according to experimental data presented in Fig. 5 are given in Table 3.

Modeling details
A finite elements domain rendering the capillary and surrounding volume was created using tetrahedral elements in three-dimensional space. The mesh element size must be chosen in such a way to be able to reproduce the specifics of the flow; here, the drop shape and consequently surface adsorption process itself. To render the specifics of the experimental cuvette setup the cubic shape space was created. The dimensions of the cube were chosen such to be able to receive and hold an investigated drop geometry. The cross section of the volume is presented at Fig. 8, the detailed view presents the capillary tip geometry. It is visible that the circular shape To reduce amount of computing work required to solve the field variables, more elements of smaller size are located near the capillary tip by applying a sizing function. More coarse mesh located in the remaining region is satisfactory to render the specifics of the surrounding flow. Appropriate boundary conditions were applied to the volume bounding faces. Typical wall boundary condition with no slip hydrodynamic condition was applied to the side and bottom faces. Top face was set as outflow boundary condition that lets to keep constant pressure in whole domain by calculating outflow and eventually backflow. The circular boundary inlet face inside the capillary was defined using velocity vector boundary condition normal to the inlet face, pointing inward the volume.

Calculation time
The process of computation involving the spatially variable interfacial tension and Gibbs free excess translocation is time consuming task. The tested case for single core (Intel Xeon X5680 at 3.33 GHz) shows that computation rate is about 0.001 calculated seconds per real-time hour. This rate is influenced by many factors where most important is the size of the mesh. Although processor affinity was not set the impact of operating system is not significant in this case because the machine is equipped with 24 processor cores and the operating system dispatcher (a task scheduler algorithm) provides effective workload balance. The rule of the thumb, which applied also in that case, showed that mesh improvement, that is a decreasing a cell size by half increases calculation time about eight times (factor 2 3 ). A computation rate for actual numerical step is shown versus simulated time at Fig. 9. At about third second of a simulated time, the sudden decrease can be noticed in a rate of calculation. This effect corresponds to drop detachment and necessity of simulation a sudden and fast moving interface and then the additional laying drop interface. The visible steep peaks are caused by the numerical convergence difficulties during calculation condition of fast changing interface during drop falling. The calculation was performed up to 9th second of simulated time that allowed to observe three consecutive drops detachments. The liquid already detached formed a laying drop at the bottom of simulated space and the mass transfer calculation for this part was still continued. In fact, the computation is highly time consuming, but the algorithm is prepared for single-core calculation in this work. The issues arising during parallelization introduce much more difficulties and will be, however, addressed in next works.

Droplet geometry
Forming a droplet through injection flow is presented by applying calculated markers onto experimental images. Figure 10 presents the stages of droplet creation for experimental volumes from 3 to 13 µL. The exact volumes are presented in Table 4. The choice of appropriate stages for graphical presentation and tabular comparison ( Table 4) was done based on experimental and calculated volumes closeness. A comparison of the geometry of the simulated drop at several presented stages indicates that the calculated shape is a little more spherical than experimental one.
The figures with droplet creation stages is made by overlapping two graphs, one which contains photograph taken from experimental setup during drop creation and the second from simulation. Special care was taken when rescaling the sizes. Both graphs contain capillary tip which is reconstructed in the simulation. The capillary was the reference element during rescaling process. The simulated graph contains marker located at interface reconstructed by VOF model. The similarity between experimental (underneath) and simulated (above) shows that marker model extension performs well in this case. There is although noticeable discrepancy at the capillary interface attachment location. The characteristic wave shape reproduced by simulation is not real interface behavior in this location. Although it has low impact on the overall droplet geometry and creation times, this evident negative phenomena is created due to the well-known numerical effect of spurious currents appearance. This has no relation to the model of markers introduced here but originates from the VOF model.

Mass transfer
The illustrative case was created to render the specifics of experiment performed. The surfactant SDS was added to the water phase entering the toluene phase at the concentration level equal 0.347 mol/m 3 . A typical situation of calculated stream fluxes and located markers is presented at high magnification in Fig. 11 for four typical computational cells. The thin lines represent the computational, volumetric cells, while the grayed areas constitute the interface. The markers are located at the centroids of a plane figures (triangle or tetragon) that represents the interface in three-dimensional space. The vectors are defined at the cells' centers and are attached to the barycenter of every tetrahedral cell containing interface.
The Gibbs surface excess evolution profiles at the interface (Fig. 12) show the process of adsorption at several stages of droplet creation. The results presented show a form of interface aging, and its development approximately evenly in all superficial directions.
The exact comparison of transient changes of surface tension with experimental data (Fig. 2) is subject to considerable difficulty, namely the calculation time. The calculated period of time includes 9 s of process which is much smaller than experimental data time period. In Fig. 13, dynamics of mean (over whole interface) surface tension and Gibbs surface excess shows the typical behavior, surface tension is decreased, while the mass flux of SDS increases Gibbs surface excess. The surface tension is calculated as the mean over the interface surface A drop by the formula:  Table 4) The same surface integral is used for the mean Gibbs excess estimation at whole drop area:  The typical situation that happens when multiphase interfaces are joining and/or separate. The simulated results present a few situations where the interface break (drop detachment) and the situation when two separate interfaces join themselves (drop hitting laying drop). The coalescing/decoalescing phenomena are presented to show the limitations of VOF method and that this limitation is inherited by the markers model. The film drainage phenomena are significant for a cases when tracking in very short time scales (bubble bouncing times can vary in the range of milliseconds depending on the species) and in interface proximity is required. The film drainage is then the reason the bubble bounces several times without interface breakup with high frequency and then eventually breaks by coalescing with the interface [88].
A droplets interface contact is presented by three consecutive stages: film drainage, film rupture, and neck growth [89]. Although successful in realistic interface simulation, the volume of fluid method is not able to reproduce the coalescence with high accuracy. The film drainage is the stage that is almost impossible to reproduce due to VOF model fundamental assumption that all phases are represented by common velocity field. During simulations when the two interfaces approach themselves, the resultant simulated effect is immediate film rupture almost without any effects of film drainage. This departure from realism can be accepted if the droplets kinetic energy is enough to rupture the interface. This is typical for bigger sizes of droplets, and during the experiments, it was a rare case to observe droplet bouncing. In fact the criteria of bouncing, coalescence and separation can be described in terms of Weber number (23), which estimates the relation of the translational kinetic energy and the surface energy of the droplets.
For the case presented, a simulated droplet radius and its velocity at the impact was calculated as integral volumetric mean over droplet volume: Interfacial tension was calculated according to the formula (21). The value of simulated mean droplet velocity at impact was 0.083 m/s, its radius 0.012 m and mean interfacial tension 0.035 N/m. Corresponding to that data, Weber number value 4.6 indicates that the bouncing stage is not expected [90], but it is still possible to happen in presented experimental system, and such behavior is sometimes observed. Consequently and due to the sharing of the velocity field among both simulated phases, the directions of vectors attached to the interface at markers locations show (Fig. 14) that film drainage do not exist in the simulated process. The interface separation is tracked by the markers according to the interface evolution and consequently their number and new locations are calculated. Because, for every time step, the calculated locations of markers depend entirely on the phase volume fraction, they are located without any numerical problems during separation according to the presented algorithm. The surface value handover using proposed nearest neighbor scheme provides continuity of surface property translocation. Proposed markers algorithm assures that no orphaned markers appear during interface rupture (see Fig. 15).

Markers distribution
The precise inspection of the three-dimensional markers locations shows their temporary tendency to accumulate. Forming such a groups seems to provide poorer interface description than when the markers locations were equally spaced. The statistical analysis of the markers locations was performed to indicate the markers distribution. Because the markers are assumed to be placed at every cell in which the interface is located the clustering tendency occurs for the case of small interface area per cell. In such situation, the interface surface is located near particular vertex of the cell which causes that adjacent cells markers are also each other at a   15 Droplet detachment in three consecutive time stages. The algorithm tracks the interface correctly during interface rupture without creating any orphaned markers smaller distance. Two ways of creating even locations could be performed; removing markers from the clusters or adding markers in the sparse regions. This problem is not addressed in this work, but statistical cluster analysis is provided (Fig. 16).
The agglomerative hierarchical clustering [91] method is most suitable for this case as one need not to assume any initial number of groups existing in the domain. This method is based on using progressive extension of a criteria forming group (e.g., distance) and creating the tree or agglomeration process graph. Such criteria for the case of analyzed markers are in this work defined based on the Euclidean and Chebyshev metrics-distance formulation. Increasing the binding distance the smaller number of bigger clusters are created. At the first stage, the small binding radius causes every marker to belong to its own group. As the distance value is increased, the more markers are included in the same sets which consequently causes creating smaller number of them in the domain. For a comparison, the equally spaced, reference surface (Fig. 16) was created and was subject of the same analysis. Both the Euclidean and Chebyshev metrics was performed to stress eventually differences between actual markers and equally spaced reference markers surface. Analysis shows (Fig. 17) that in both metrics cases the arising differences are generally similar. It is visible that most statistical differences between actual markers and reference ones occur at the smallest distance levels. This is obvious after observing animated view of the markers evolution-most clusters are formed when markers approach the cell vertexes. This happens when the interface is located near any vertex. It is also visible that most clusters' sizes, which present most deviation from ideal surface, are much below the mean mesh size. The strongest deviation from the ideal, equally spaced surface occurs at the smallest distances, smaller than cell size  level; so agglomeration tendency can be regarded as a model imperfection. But this is typical for any discrete method of problem solving that any physical interpretation of solution below discretization size (here cell size) becomes unreliable. The most important assumption, that every cell containing the interface should also contain one marker, is conserved. The jumps of the dashed lines representing equally spaced reference markers are also easily explained by the sudden markers agglomeration during hierarchical clustering process due to their equidistant character. It is then clearly visible that for sizes above the mean cell size (after first dashed line jump), the character of both the ideal and actual markers become similar and statistically corresponding to each other. Consequently, it can be stated that actual markers locations are equally spaced at the cell size level.

Conclusions
A new calculation algorithm was proposed which is designed for selected property or variable transportation along with the spatial and time evolution of predefined isovalue variable. The isovalue should indicate the location of some specific discontinuity, e.g., fluid interface or shock wave. For the test of the algorithm proposed, the multiphase system was chosen in which interface is present.
The experimental representation of the problem in which proposed modeling technique can be used is the process of adsorption at the water-organic liquid interface. The computational fluid dynamics models does not yet account for such specific interface processes. The surface excess, thermodynamic property describing amount of surfactant adsorbed at liquid interface, is an analogue of concentration. Surface excess has direct impact on the experimentally measurable surface tension where the higher excess values correspond to the lower surface tension, in the case of amphiphilic substances. The model proposed enables to track the evolution of any superficial property but in the case of surface excess it becomes more important because it opens a way to simulate cases unresolvable by CFD so far.
The results presented quantitatively describe the evolution of adsorption process and how it affects the key superficial property-surface tension. The experimental measurements show that the VOF model markers extension is not only suitable to render the geometry of the interface, the droplet volume and its creation time, but it also enables to observe dynamic changes in the character of the interface. It is not possible to illustrate local surface tension variation by experimental droplet creation methods (pendant drop method, drop shape method) which are widely used in the research field, but imaging its local changes by means of calculations lead to better description of complexity of the adsorption process. The markers model proposed in this work enables such superficial description and can be useful in the research process. In Fig. 13, calculated symmetry of the surface properties distribution is reasonable due to the symmetry of the droplet geometry. The results obtained for Gibbs surface excess and surface tension are in good accordance with experimental measurements. The markers extension to the well-known VOF model is beneficial in multiphase simulation because it gives a reasonable balance between simplicity and accuracy.
The entire proposed model is then the composite of two parts: calculating the interfacial adsorptive mass flux and evolution of surface properties using markers. Presented here adsorptive extension of volume of fluid model includes computation of interface mass fluxes, modeling spatially varying surface excess and surface tension and dynamic contact angle handling. In fact, the obtained results are burdened with error, but the exact nature of the flow and the behavior of the interface is maintained. The computed values of droplet volume and its geometry at different stages agree with the presented experimental values.
The main difficulty of the model presented is its high computational cost that reveals in long calculation times. While such negative attribute of the model is troublesome, it does not affect the quality of solutions obtained and ability to take into account the interface mass transfer phenomena.