Single Contaminated Drops Falling through Stagnant Liquid at Low Reynolds Numbers

: Numerical simulations of contaminated spherical drops falling through a stagnant liquid at low Reynolds numbers are carried out using the ﬁnite difference method. The numerical results are used to describe the behavior of the surfactant concentrations and to understand the surfactant effects on the ﬂuid motions in detail. The predicted interfacial surfactant concentration, Γ , is almost zero for angles, θ , below a certain value (the stagnant-cap angle, θ cap ), whereas it steeply increases and reaches a large value for θ > θ cap (the stagnant-cap region). The increase in the initial surfactant concentration, C 0 , in the drop enhances the adsorption from the drop to the interface, which results in the increase in Γ and the decrease in θ cap . Peaks appear in the predicted Marangoni stresses around θ cap , which causes similar peaks in the pressure distribution. The high-pressure spots prevent the ﬂuid motion along the interface, which results in the formation of the stagnant-cap region and the attenuation of the tangential velocity in the continuous phase. The surfactant ﬂux from the bulk to the interface decreases C in the vicinity of the interface for θ < θ cap and the weak diffusion cannot compensate for the reduction in C by adsorption, which results in C at the interface smaller than C 0 . The pattern of the low C region is determined by the advection and does not smear out because of a small diffusive ﬂux.


Introduction
Surface-active agents (surfactants) have various effects on the motion of a fluid particle, e.g., the drag coefficient of a spherical bubble under a fully contaminated condition corresponds to that of a solid particle [1,2], capillary waves formed on the surface of a bubble are dampened by the presence of surfactant [3,4], and so on. Boussinesq [5] assumed that the interface contaminated with surfactant behaves similar to a viscous membrane and modeled the increase in drag caused by the presence of surfactant as surface viscosity [1,6]. The drag model for clean spherical fluid particles of low Reynolds numbers, i.e., the Hadamard-Rybczynski model [7,8], can be easily extended to that for contaminated fluid particles by introducing the surface viscosity. On the other hand, Levich [9] succeeded in explaining the surfactant effect on the drag force by accounting for the surface tension gradient, i.e., the Marangoni stress, caused by non-uniform distribution of the surface tension due to the adsorption of surfactant to the interface. Numerical simulations of contaminated fluid particles have been carried out so far by taking into account the Marangoni stress [2,[10][11][12][13][14][15][16][17][18][19][20], in contrast, measurements of the Marangoni stress and the interfacial surfactant concentration have rarely been carried out.
Reduction in the surface tension due to surfactant adsorption takes place for various surfactants such as Triton X-100, sodium dodecyl sulfate (SDS), alcohols such as 1-pentanol and so on. In the presence of interfacial flow, the distribution of surfactant at the interface is non-uniform due to advection. For example, the interfacial flow formed at the interface of a drop falling through a stagnant liquid due to gravity transfers surfactant molecules adsorbing to the interface from the nose side to the rear side of the drop, which results in a high surfactant concentration in the rear side. The non-uniform surfactant distribution also results in a non-uniform distribution of surface tension. The distribution of the surface tension can be calculated if the Marangoni stress is measured, and the interfacial surfactant concentration can also be computed from the surface tension by making use of the Langmuir-Szyszkowski isotherm [21]. The Marangoni stress can be evaluated if the velocity gradients at the interface and the viscosities of the two phases are known since the jump of the tangential viscous stress at the interface balances with the Marangoni stress. Hosokawa et al. [22] measured velocity fields about single contaminated drops of a glycerol-water solution falling through a silicone oil at low drop Reynolds numbers using SFV (spatiotemporal filter velocimetry) [23] to evaluate the distributions of the interfacial surfactant concentration ( Figure 1). However, other field variables, e.g., the surfactant concentration in the drop and the pressure fields, required for full understanding of the phenomena, have not been obtained in the experiment yet.
In this study, the motions of contaminated drops of low Reynolds numbers are numerically simulated using the finite difference method under the same conditions as those in Hosokawa et al. [22] to discuss the distributions of the bulk and interfacial surfactant concentrations and the surfactant effects on the interface velocity and pressure in detail. In the following sections, the numerical method is first introduced. The validity of the numerical method is examined for clean drops, spherical particles and contaminated bubbles [10]. Then, the falling motions of contaminated drops are simulated, and the surfactant effects are discussed.  [22], where d is the drop diameter, a is the drop radius, r and θ are the radial and azimuthal coordinates, U is the drop velocity, and g is the magnitude of the gravitational acceleration.

Governing Equations
The governing equations are written in the frame of reference moving with a single drop. The drop shape is assumed to be spherical, and the drop Reynolds number is given as a parameter from the experimental data [22]. The flow passing through the spherical drop at the prescribed Reynolds number is therefore obtained by solving the governing equations. The (r, θ) two-dimensional spherical coordinates are used since the flow has axial symmetry at low drop Reynolds numbers. The flow is assumed to be incompressible and isothermal. The Navier-Stokes equation for Newtonian fluids is therefore given by [24] ρ k ∂u r,k ∂t where u r and u θ are the velocity components in the r and θ directions, respectively, t is the time, p is the pressure, ρ is the density, µ is the viscosity, and k denotes either the dispersed phase (k = D) or the continuous phase (k = C). The differential operators, u · ∇ and ∇ 2 , are given in Appendix A. The following continuity equation holds in each phase: The normal component of the velocity is continuous at the interface when there is no interfacial mass transfer. Hence the velocity boundary condition at the drop interface is given by u r,D = u r,C = 0 (4) and the following no-slip condition is used for the tangential velocity component: where the subscript int denotes the interface. The tangential viscous stress, τ rθ , in each phase balances with the Marangoni stress at the interface: where [X] int represents the jump of the quantity, X, at the interface, σ is the surface tension, and a is the drop radius.
In the present study, surfactant is assumed to be soluble in the dispersed phase, so that it is present inside the drop at the molar concentration, C. The transport equation of C is given by where D D is the diffusion coefficient. Adsorption and desorption of surfactant molecules take place at the interface and the magnitude of the adsorption-desorption flux depends on C at the interface, C S , and the molar concentration, Γ, of surfactant adsorbing to the interface. The transport equation of Γ is given by [2,10,25] ∂Γ ∂t where D S is the diffusion coefficient of surfactant at the interface, ∇ S is the surface gradient operator, and u S is the interface velocity. TheṠ Γ is the net adsorption-desorption flux, which is evaluated by making use of the following Frumkin-Levich model [9,10,26]: where k a is the adsorption rate constant, k d is the desorption rate constant, and Γ max is the saturation value of Γ. The diffusive flux of C at the interface balances withṠ Γ , i.e., The surface tension, σ, of the contaminated interface depends on Γ as [21] where σ 0 is the surface tension of a clean interface, R G is the universal gas constant, and T is the temperature.
The following dimensionless forms of the governing equations are useful to make clear the relevant dimensionless groups: where u * = u/U, p * = p/ρU 2 , C * = C/C 0 , Γ * = Γ/Γ max , r * = r/d, t * = Ut/d, ∇ * = d∇, ∇ * S = d∇ S , and C 0 is the initial concentration. The Reynolds numbers, Re k , the Peclet number, Pe D , and the surface Peclet number, Pe S , are defined by Re k = ρ k Ud/µ k , Pe D = Ud/D D , and Pe S = Ud/D S , respectively, where d is the drop diameter, and U is the drop velocity. The momentum jump condition in the tangential direction is given by where the Marangoni number, Ma, and the dimensionless surface tension, σ * , are defined by Ma = R G TΓ max /µ C U and σ * = σ/R G TΓ max , respectively. The dimensionless radius is given by a * = a/d = 1/2. The interfacial surfactant flux in the dimensionless form is given byṠ * where C * S = C S /C 0 , La is the Langmuir number defined by La = k a C 0 /k d , and α is the ratio of the adsorption rate to the characteristic velocity, i.e., α = k a C 0 d/U. The boundary condition of C at the interface is where K is the dimensionless adsorption length defined by K = Γ max /C 0 d. The σ * is given by where σ * 0 = σ 0 /R G TΓ max .

Velocity-Pressure Coupling
The present numerical method for obtaining the velocity and pressure fields is based on the fractional step method [27,28]. The Adams-Bashforth scheme is used for the advection terms of the Navier-Stokes equations, while the Crank-Nicolson scheme is used for the diffusion terms. Therefore, where ∆t * is the time step size and The superscript n denotes the time step, A and D are the advection and diffusion terms, respectively, and S is the source term in the Crank-Nicolson solver. Equations (20) and (21) are solved for u m , where the superscript m is the step number for the iteration of the Crank-Nicolson solver. Since the former and latter include the terms of u θ and u r , respectively, in the diffusion term, these terms are treated as source terms in each iteration. The converged solutions of the velocity components are denoted as u † .
The pressure contribution is then subtracted from u † to obtain the intermediate velocity The pressure at n + 1 is obtained by solving the following Poisson equation: where The velocity components are then updated as The staggered variable arrangement is used, i.e., the velocity components are located at the faces of computational cells, whereas the scalar variables are located at the cell centers. The advection terms are discretized using the ENO (essentially non-oscillatory) scheme [29]. The second-order centered difference scheme is used for the diffusion terms.

Interface Velocity and Numerical Treatment at Pole
The u θ,int is updated at each time step after obtaining u n+1 k . The momentum jump condition in the tangential direction is given by where κ is the viscosity ratio defined by κ = µ D /µ C . Applying the continuity of the velocity at the interface and the second-order one-sided finite difference to the velocity gradients for the configuration in Figure  Discretization of some derivatives in the governing equations needs a special treatment at the pole (r = 0). In the evaluation of the divergence of u and the Laplacian of the velocity components and the pressure, no special treatments are required since the terms rφ and r∂φ/∂r vanish at the pole, where φ = u r , u θ or p. On the other hand, for instance, the finite difference approximation of the advection term u r ∂u θ /∂r requires u θ at r = 0. A numerical treatment similar to that in Mohseni and Colonius [30] and Sugioka and Komori [31] is used in this study to set u θ at r = 0, i.e., u θ (0, θ) = (u θ (∆r/2, θ) + u θ (∆r/2, θ + π))/2 ( Figure 3). Although 0 ≤ θ ≤ π, u θ (∆r/2, θ + π) can be set to the value of u θ (∆r/2, θ − π) because of the axial symmetry.

Surfactant Transport Equations
The transport equation of Γ * is solved using the second-order Runge-Kutta scheme: The Runge-Kutta scheme is also used for the transport equation of C * as follows: The advection and diffusion terms are discretized with the second-order ENO and secondorder centered difference schemes, respectively. The bulk concentration, C * S , at the interface is evaluated by applying the second-order one-sided finite difference to ∂C * /∂r * , and is used to evaluateṠ * Γ . TheṠ * Γ is then used to set the boundary condition, Equation (18). The u * θ,int is updated for ∇ * S σ * (Γ * n+1 ).

Solution Procedure
The solution procedure of the numerical method is summarized in the following: 1. Set the dimensionless groups (ρ D /ρ C , κ, Re C , La, Pe D , Pe S , Ma, α, K and σ * 0 ); generate computational grids inside and outside the drop; and set the initial conditions for each field variable.
Compute the surfactant concentration, C * n S , at the interface using Equation (40) with C * n . 9.

Clean Drop and Solid Sphere
The numerical method is validated for simple problems with well-accepted correlations, i.e., uniform flows about solid spheres and clean spherical drops. Figure 4 shows the computational domain and the grid system. The dimensionless particle diameter is unity, and the domain size is 40 in the radial direction. The outer boundary of the domain is split into either inflow or outflow. The zero-gradient condition is used for the velocity at the outflow boundary and the outlet pressure is set to zero. The numbers of the computational cells in the r and θ directions are 25 and 120, respectively. The cell size is uniform in the θ direction, whereas in the r direction the cell size decreases (increases) toward the drop interface (the outer boundary of the domain) as r increases. For clean drops, it is not necessary to use the very thin computational cells in the vicinity of the drop interface for capturing the velocity boundary layer at low and moderate Reynolds numbers. On the other hand, in the contaminated drop cases discussed in Section 4, the thin cells are required to capture concentration boundary layers of C * forming at the interface. Cuenot et al. [10] estimated the dimensionless thickness of the concentration boundary layer as δ ∼ π/(4Pe C ) [9] and generated computational grids so as to have at least three computational cells within the boundary layer. The present grid is prepared to meet the same requirement. This grid will also be used in the contaminated drop simulations. The mesh dependence of the predictions was confirmed to be small enough at the present spatial resolution. Figure 5 shows predicted u * θ,int of clean drops at several drop Reynolds numbers, Re C , the range of which (0.1 ≤ Re C ≤ 20) covers Re C in the experiment of the contaminated drops [22], i.e., 0.5 ≤ Re C ≤ 0.73. The viscosity ratio in the simulations is κ = 0.021. The predicted u * θ,int at Re C = 0.1 agrees well with the Hadamard-Rybczynski (HR) solution, u * θ,int = sin θ/(2(1 + κ)). The slight difference between the prediction and the analytical solution for the Stokes flow is attributed to the presence of a weak inertial effect even at the small Reynolds number. The increase in Re C up to Re C = 1 makes the inertial effect obvious, i.e., u θ,int for θ < 3π/4 is larger than the HR solution. The deviation from the Stokes flow is increased by further increase in Re C . Figure 6 shows predicted C D of clean spherical drops. The open circles are for κ = 0.021 and they are compared with the HR drag (the black dashed line) and the empirical correlation proposed by Myint et al. [32] (the black solid line)  The C D at Re C = 0.1 agrees with both drag correlations. The predictions become larger than the HR drag as Re C increases due to the inertial effect. The trend of the drag curve of Equation (42) agrees with that of the data, the values of which are however smaller than those of the correlation since the correlation accounts for the increase in drag by shape deformation. The green and blue symbols in the figure are predictions for κ = 1 and 5, respectively. The simulations reproduce well the increase in C D by the factor of (2 + 3κ)/(1 + κ). Simulations of flows about spherical particles are also carried out as the limiting case of κ → ∞ by directly applying the no-slip boundary condition to the particle surface. The predicted drag coefficients are shown in the figure by the red symbols and compared with the following empirical correlation proposed by Schiller and Naumann [33]: The data agree well with the correlation. Thus, the present numerical method gives good predictions of the motion of drops for wide ranges of Re C and κ.

Contaminated Bubbles
Cuenot et al. [10] carried out numerical simulations of contaminated spherical bubbles in stagnant liquid at Re C = 100, Pe C = 1 × 10 5 and K = 1. They changed Ma, α and La to deal with different situations, i.e., the stagnant-cap regime and the fully covered regime. Case 1 (Ma = 61, α = 0.001, La = 0.112), Case 3 (Ma = 61, α = 0.001, La = 0.0112) and Case 4 (Ma = 610, α = 0.001, La = 0.112) in their study are simulated with the present numerical method. The flows inside bubbles are not solved; instead, the tangential viscous stress in the gas phase is set to zero in Equation (6). Figure 7 shows predicted u * θ,int and Γ * , which are compared with the predictions by Cuenot et al. [10]. The u * θ,int increases with increasing θ, whereas it steeply decreases to zero at θ ∼ π/4 and 5π/8 in Cases 1 and 3, respectively. On the other hand, the gradients of Γ * are large at these θ. The steep decrease in u * θ,int is due to the Marangoni stress caused by the gradient of Γ * . In comparison with Case 3, Γ * in the contaminated region is larger in Case 1 and the attenuation of u * θ,int is more significant because of the larger La. In Case 4, u * θ,int is almost zero, in other words the interface is immobile, due to the large Ma, and Γ * ≈ 0.09 even at the bubble nose. Thus, being similar to the predictions of Cuenot et al. [10], the effects of La and Ma on u * θ,int and Γ * are clearly observed.

Simulations of Contaminated Drops
The main difference in the simulations of contaminated drops carried out in this section is the value of La, which ranges from 0 to 10 2 , leading to varying surface concentrations. The fluid properties of a drop and a bulk liquid are the same as those in the experiment [22], i.e., those for a glycerol-water solution of 54 wt.% and a silicone oil. They are summarized in Table 1. The drop diameter is 8.3 mm. Under this condition, Re C is less than unity. Triton X-100 is used for surfactant. The adsorption-desorption properties and the diffusion coefficients of Triton X-100 are shown in Table 2 (see Appendix B for evaluation methods). The C 0 in a drop is initially set either to 0 (clean liquid), 0.0020, 0.0050, 0.010 or 0.10 mol/m 3 . These conditions are referred to as Cases 1, 2, 3, 4 and 5, respectively. The concentrations are less than the critical micellar concentration (0.28 mol/m 3 ). The measurements were carried out at 250 mm below the nozzle tip. The dimensionless time when drops reached the measurement region was 30. The relevant dimensionless groups in each case are summarized in Table 3. These values of the dimensionless groups are used as the input parameters of the numerical simulations. See Hosokawa et al. [22] for further details of the experimental conditions.    Figure 9 shows Γ * at t * = 10, 20 and 30. In Case 2 (C 0 = 0.0020 mol/m 3 ), Γ * increases in the rear region of the drop as t * increases. The contaminated area also increases, but is not so large even at t * = 30. In Case 3 (C 0 = 0.0050 mol/m 3 ) and 4 (C 0 = 0.010 mol/m 3 ), the trend of the temporal change in Γ * is similar to that in Case 1, whereas the contaminated area becomes much wider as C 0 increases. The Γ * in Case 5 (C 0 = 0.10 mol/m 3 ) is much larger than the other cases, and the whole interface is contaminated. The presence of surfactant reduces the region of the internal circulation both in the experiment and simulation of Case 3 as shown in Figure 8b (the flow field of Case 2, whose flow characteristics are in-between those of Cases 1 and 3, was omitted for saving the space). The magnitude of the drop velocity in the vicinity of the drop rear is very small. The increase in C 0 from 0.0050 mol/m 3 to 0.010 mol/m 3 (Case 4) makes the circulation region smaller (Figure 8c). The interface velocity in the rear half, θ > π/2, is close to zero, so that the formation of velocity boundary layer along the interface can be clearly observed in the continuous phase. In Case 5 (Figure 8d), no internal circulation is formed inside the drop. In the continuous phase, the grayscale values near the interface in the measurement are very low (almost black), which implies that the velocities there are very small. The prediction agrees with the measurement; the interface velocity is almost zero at the whole interface; in other words, the interface is fully contaminated and immobile, so that the velocity field outside the drop seems to be similar to that about a solid sphere.  Figure 10a shows the distribution of C * in Case 2. The t * increases from the left to the right. The rightmost figures are magnified views of the regions labeled by (i) and (ii). The C * near the interface in the rear region decreases at t * = 5. The low C * region develops toward the centerline of the drop as t * increases due to the advection. However C * near the rear does not decrease so much. Because of the low diffusion coefficient, the pattern of the low C * region is determined by the advection and does not smear out. Figure 11 shows the time evolution of C * S in Case 2. The surfactant flux from the bulk to the interface decreases C * in the vicinity of the interface and the weak diffusion cannot compensate for the reduction in C * by adsorption. The C * S at t * = 5 is therefore smaller than the initial value, C * S = 1. The θ for the minimum C * S corresponds to that for steep velocity decrease; in other words, the stagnant-cap angle, θ cap . For θ < θ cap , Γ * is very small, so that the adsorption is dominant in the surfactant flux. The increase in C * S with increasing θ for θ > θ cap implies that the desorption increases C * S in the stagnant-cap region. In the stagnant-cap region, Γ * approaches a value for the local adsorption-desorption equilibrium, Γ * ∼ 1/(1 + (C * S ) −1 ). The u * θ,int is very small in that region, but is still non-zero. The weak interface compression effect (Γ∇ S · u S ) appears and slightly increases Γ * . However, some excess amount of surfactant is desorbed from the interface to the bulk. Thus, C * S has an increasing trend in the high θ region. At larger t * , C * S decreases moderately for θ < θ cap , whereas for θ > θ cap C * S increases and becomes larger than unity.   Figure 10b shows the time evolution of C * in Case 4. As shown in the magnified view of C * near the interface (iii), a low C * region is formed around the drop equator. The interface area is more than half covered with surfactant at this moment (Figure 9c). The formation of the low C * region is therefore the result of surfactant adsorption from the bulk to the interface. The thickness of the concentration boundary layer is, as expected, very thin because of the high Pe D . In Case 5, the whole interface is covered by surfactant as shown in Figure 9d, which results in the formation of the thin concentration boundary layer of low C * in the nose region (Figure 10c(v)). As it is similar to Case 2, the interface compression effect increases Γ * in the rear half of the drop and the desorption becomes dominant near the rear, resulting in the formation of boundary layer of C * larger than unity in the rear region (the magnified view (iv)).  Figure 12 shows the pressure distribution in Case 3, in which p * D and p * C are subtracted by the pressure, p * 0 , at the drop center and the far field pressure, p * ∞ = 0, respectively. The θ cap represents the cap angle evaluated from the velocity distribution and θ cap ∼ 0.6π rad. In the pressure distribution, high-pressure spots are formed in the vicinity of the interface at θ slightly smaller than θ cap . Comparing p * with the velocity field clearly shows that these high-pressure spots prevent the fluid motion along the interface, which results in the formation of the stagnant-cap region in the rear half of the drop and the attenuation of the tangential velocity in the continuous phase.  Figure 13a shows the interface velocity, u * θ,int , at t * = 30. The lines and symbols are the predictions and the measured data [22], respectively. The predicted velocity of the clean drop (Case 1) agrees fairly with the data, except for the rear region. The lower values at large θ in the data are due to slight contamination in the rear region, the cause of which would be contamination brought into the drop by seeding particles. The retardation effect of surfactant in u * θ,int is clearly seen in both data and predictions of Cases 2-5. In Case 2, the predicted velocity profile agrees with the data. The u * θ,int decreases from θ ∼ π/2 as θ increases and becomes very small for θ > 3π/4, i.e., θ cap ∼ 3π/4. The increase in C 0 decreases θ cap , i.e., θ cap are 0.6π and 0.45π rad in Case 3 and 4, respectively. Although the predicted velocities are somewhat larger than the data, the dependence of u * θ,int on C 0 is reproduced. In Case 5, u * θ,int ∼ 0 at any θ as already seen in Figure 8. The distributions of Γ * are shown in Figure 13b. The predicted Γ * of Case 2 is almost zero for θ < θ cap , whereas it steeply increases and reaches a large value for θ > θ cap , which is larger than La/(1 + La) = 0.67 for the absorption-desorption equilibrium for C 0 since C * S in this region is larger than unity ( Figure 11). The measured Γ * shows a similar profile, while the following differences are present: the measured Γ * is smoother around θ = θ cap than predicted and approaches 0.67. The increase in C 0 (Cases 3 and 4) increases Γ * , whereas the drop nose is still almost clean. On the other hand, Γ * takes a large value even at the nose in Case 5 since the adsorption is much more dominant than the desorption as represented by the large La.

Nose
As shown in Figure 13c, peaks appear in the predicted Marangoni stresses in Cases 2, 3 and 4. However, in Case 5, there are no peaks. This is due to the smaller Γ * gradient. Despite a smaller Γ * gradient, the Marangoni stress immobilizes the interface almost completely since the Marangoni stress acts on the whole interface. In this situation, the Marangoni stress agrees very well with the viscous stress acting on the solid sphere given by the Stokes solution (the black dotted line in the figure), i.e., ∂u * θ /∂r * r * =a * = 3 sin θ. The smaller peaks in the data are attributed to a smoother distribution of Γ * .
The peaked distribution of the Marangoni stress also causes similar peaks in the pressure distribution as shown in Figure 13d. The locations of the strong Marangoni stress correspond to those of the steep velocity reduction (the cap angles). The peaks of the pressure however locate at slightly smaller θ. This trend was already observed in Figure 12. The pressure distribution in Case 5 is smooth as in the clean drop case and agrees well with the Stokes solution given by p * C = 3 cos θ/(2Re C ). In the distributions of C * S shown in Figure 13d, θ for the minimum C * S in Cases 2-4 roughly correspond to θ cap . The C * S in these cases exhibit the increasing trend in the stagnant-cap region. On the other hand, on the completely immobilized interface in Case 5, the profile of C * S monotonically increases with increasing θ. The measured Γ * is much smoother than predicted, which suggests that in reality Pe S is much lower than that used in the simulation or that there is some other transport process much more effective for interface transport. A brief discussion on this point is given in Appendix C.  Figure 13. Predictions of (a) u * θ,int , (b) Γ * , (c) −Ma∂σ * /∂(a * θ), (d) p * C and (e) C * S at interface. The symbols and lines are measured [22] and predicted values, respectively. The black dotted line in ∂σ * /∂(a * θ) represents the dimensionless velocity gradient, ∂u * θ /∂r * r * =a * = 3 sin θ, given by the Stokes solution for solid sphere. The black dotted line in p * C is also the Stokes solution given by p * C = 3 cos θ/(2Re C ).

Conclusions
The motion of contaminated spherical drops falling through a stagnant liquid were numerically simulated using the finite difference method to investigate the effects of surfactant on the flow inside and outside the drops. The numerical method was first validated through simulations of spherical particles, clean drops and contaminated bubbles. The predicted drag coefficients of the spherical particles and clean drops agreed well with available drag correlations. The interface velocity and the interfacial surfactant concentration of contaminated bubbles also agreed with those predicted in the literature [10]. The numerical simulations of contaminated drops were then carried out. The numerical conditions were the same as those in the experiment carried out in Hosokawa et al. [22], i.e., the fluid properties for the drop and the bulk liquid were for a glycerol-water solution and a silicone oil, respectively, and the surfactant was Triton X-100. The drop Reynolds numbers were less than unity and the drop shape kept spherical. The surfactant concentration, C 0 , was varied below the critical micellar concentration; C 0 = 0.0020, 0.0050, 0.010 and 0.10 mol/m 3 in Cases 1, 2, 3, 4 and 5, respectively. The characteristics of the field variables can be summarized as follows: • The predicted interfacial surfactant concentration, Γ, is almost zero for the angle, θ, smaller than a certain value (θ cap ), whereas it steeply increases and reaches a large value for θ > θ cap . The interface velocity, u θ,int , decreases as θ increases and becomes very small for θ > θ cap . The increase in C 0 increases Γ and decreases θ cap . The presence of surfactant attenuates the internal circulation and the increase in C 0 makes the circulation region smaller. • The surfactant flux from the bulk to the interface decreases C in the vicinity of the interface and the weak diffusion cannot compensate for the reduction in C by adsorption. The bulk concentration, C S , at the interface therefore tends to be smaller than C 0 for θ < θ cap . The pattern of the low C region is determined by the advection and does not smear out because of a small diffusive flux. • Peaks appear in the predicted Marangoni stresses in Cases 2-4, while in Case 5 no peaks develop due to a smaller Γ gradient. The peaked distribution of the Marangoni stress also causes similar peaks in the pressure distribution. The locations of the strong Marangoni stress correspond to θ cap , whereas the peaks of the pressure appear at slightly smaller θ. The high-pressure spots prevent the fluid motion along the interface, which results in the formation of the stagnant-cap region in the rear half of the drop and the attenuation of the tangential velocity in the continuous phase.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Differential Operators
The differential operators used in the governing equations are defined by

Appendix B. Surfactant Property
The isotherm is given in terms of C 0 : where χ = k d /k a At large C 0 , this equation can be approximated as The Γ max is therefore calculated by fitting the expression to data of σ(C 0 ), which can be measured using a sessile drop method ( Figure A1). Substituting Γ max and the interfacial surfactant concentration in the adsorption-desorption equilibrium state into Equation (11) yields σ(C 0 ) = σ 0 + R G TΓ max ln χ χ + C 0 (A7) Fitting this functional form to the σ data determines χ.
By eliminating the advection and diffusion terms in Equation (8) for a static interface with a uniform Γ distribution, Equation (8) can be immediately integrated, i.e.
where t is the time measured from the instant at which a drop intrudes into the bulk liquid. Substituting thus obtained Γ(t) into Equation (A3) yields σ(t) = σ 0 + R G TΓ max ln χ + C 0 e −(χ+C 0 )k a t χ + C 0 (A9) where only k a is the unknown constant. The least-square fitting to the data of σ(t) gives k a . The change of σ in time at a stationary interface is given by [34] σ(t) = σ(Γ eq ) + R G TΓ 2 The D D can therefore be calculated by applying the sessile drop method. Since there is no knowledge of D S of Triton X-100 at the interface between the glycerol-water solution and silicone oil, it is assumed that D S = D D [2,10].

Appendix C. Contaminated Drops of Lower Peclet Numbers
The measured Γ * is much smoother than predicted as shown in Figure 13, which leads to a speculation that in reality the diffusion of surfactant is much larger; in other words, the Peclet number is much lower than that in the simulation. As described in Appendix B, D D was evaluated from the data of σ obtained by the sessile drop method (Equation (A10)). By considering the measurement uncertainties in D D we increased the value of D D in the simulation by three times. The Pe D for the increased D D is given in Table A1. The simulations of the contaminated drops were carried out with these Pe D and the assumption of Pe S = Pe D . The predicted Γ * profiles were still much sharper than those in the measurement although the reduction in the Peclet numbers slightly affected the results. By further decreasing Pe S down to the order of O(10) (Table A1), we obtained better agreements between the predictions and the data as shown in Figure A2. The peaks in the pressure and the Marangoni stress observed in Figure 13 are mitigated.
The simulation thus implies the possibility of lower Pe in reality. However, according to values of D D and D S measured for different contaminated systems [35,36], the assumption, D D = D S , is considered to be acceptable, and Pe S of O(10) seems too small. Further studies are required to make clearer the cause of the deviation from the measurement.   Table A1). The legends are the same as those in Figure 13.