A mathematical model is formulated and solved for the two-phase flow of a viscous drop or inviscid bubble in an immiscible, viscous surrounding fluid in the zero Reynold's number or Stokes flow limit. A surfactant that is present on the interface is also soluble in the exterior fluid, and the drop is deformed by an imposed linear flow. The geometry is two-dimensional and Cartesian.
The dissolved surfactant is considered in the physically realistic limit of large bulk Péclet number. That is, it convects and diffuses as a passive scalar in the bulk flow where the ratio of its convection to diffusion is large. The large bulk Péclet number limit presents a significant challenge for traditional numerical methods, since it implies that large gradients of bulk surfactant concentration can develop in a spatially narrow boundary or transition layer adjacent to the drop interface. The layer structure needs to be resolved accurately so that the bulk-interface surfactant exchange, surface surfactant concentration, and interfacial surface tension can be evaluated to determine the drop's dynamics and evolution.
To resolve this computational difficulty, the dynamics of the transition layer are modeled by a leading order singular perturbation reduction of the conservation law for dissolved surfactant that is derived in the large bulk Péclet number limit.
Two versions of the boundary integral equation for two-phase Stokes flow are presented and used as an underlying fluid solver in the absence of surfactant effects. To evaluate the influence of soluble surfactant the boundary integral solver is coupled to the evolution of surface surfactant concentration and the transition layer equation. The transition layer equation is first solved by a mesh-based numerical method that has a spectrally accurate equal arc length frame for discretization of the interface and a second order time-step. Results of numerical simulations are presented for a range of different physical parameters. An analytical solution of the transition layer equation by a Green's function representation is also derived, which leads to a second, mesh-free algorithm. Numerical results of the mesh-based and mesh-free methods are compared.