A fast algorithm for simulating vesicle flows in three dimensions

Vesicles are locally-inextensible fluid membranes that can sustain bending. In this paper, we extend the study of Veerapaneni et al. [S.K. Veerapaneni, D. Gueyffier, G. Biros, D. Zorin, A numerical method for simulating the dynamics of 3D axisymmetric vesicles suspended in viscous flows, Journal of Computational Physics 228 (19) (2009) 7233–7249] to general non-axisymmetric vesicle flows in three dimensions.     

SHAPE EQUATION AND SHAPES OF AXISYMMETRIC VESICLES IN TOROIDAL TOPOLOGY

For toroidal topology we solve the shape equation for axisymmetric vesicles numerically. The phase diagram is found to be similar to that from another shape equation. This similarity is the result of the insensitive dependence of beading energy upon detailed shape of vesicles and the constant volume and area ensemble we are considering. We argue that the very small distance between two opposite cusps of the sickle-shaped vesicle makes it unstable and fuse into two encompassed spheres which might be observed in experiment. The spontaneous curvatures of the observed Clifford tori are also estimated.     

A boundary integral method for simulating the dynamics of inextensible vesicles suspended in a viscous fluid in 2D

We present a new method for the evolution of inextensible vesicles immersed in a Stokesian fluid. We use a boundary integral formulation for the fluid that results in a set of nonlinear integro-differential equations for the vesicle dynamics. The motion of the vesicles is determined by balancing the non-local hydrodynamic forces with the elastic forces due to bending and tension. Numerical simulations of such vesicle motions are quite challenging. On one hand, explicit time-stepping schemes suffer from a severe stability constraint due to the stiffness related to high-order spatial derivatives and a milder constraint due to a transport-like stability condition. On the other hand, an implicit scheme can be expensive because it requires the solution of a set of nonlinear equations at each time step. We present two semi-implicit schemes that circumvent the severe stability constraints on the time step and whose computational cost per time step is comparable to that of an explicit scheme. We discretize the equations by using a spectral method in space, and a multistep third-order accurate scheme in time. We use the fast multipole method (FMM) to efficiently compute vesicle–vesicle interaction forces in a suspension with a large number of vesicles. We report results from numerical experiments that demonstrate the convergence and algorithmic complexity properties of our scheme.     

Numerical investigation of the deformation mechanism of a bubble or a drop rising or falling in another fluid

A numerical method for simulating the motion and deformation of an axisymmetric bubble or drop rising or falling in another infinite and initially stationary fluid is developed based on the volume of fluid （VOF） method in the frame of two incompressible and immiscible viscous fluids under the action of gravity, taking into consideration of surface tension effects. A comparison of the numerical results by this method with those by other works indicates the validity of the method. In the frame of inviseid and incompressible fluids without taking into consideration of surface tension effects, the mechanisms of the generation of the liquid jet and the transition from spherical shape to toroidal shape during the bubble or drop deformation, the increase of the ring diameter of the toroidal bubble or drop and the decrease of its cross-section area during its motion, and the effects of the density ratio of the two fluids on the deformation of the bubble or drop are analysed both theoretically and numerically.     

Computation of multiphase systems with phase field models

Phase field models offer a systematic physical approach for investigating complex multiphase systems behaviors such as near-critical interfacial phenomena, phase separation under shear, and microstructure evolution during solidification. However, because interfaces are replaced by thin transition regions (diffuse interfaces), phase field simulations require resolution of very thin layers to capture the physics of the problems studied. This demands robust numerical methods that can efficiently achieve high resolution and accuracy, especially in three dimensions. We present here an accurate and efficient numerical method to solve the coupled Cahn–Hilliard/Navier–Stokes system, known as Model H, that constitutes a phase field model for density-matched binary fluids with variable mobility and viscosity. The numerical method is a time-split scheme that combines a novel semi-implicit discretization for the convective Cahn–Hilliard equation with an innovative application of high-resolution schemes employed for direct numerical simulations of turbulence. This new semi-implicit discretization is simple but effective since it removes the stability constraint due to the nonlinearity of the Cahn–Hilliard equation at the same cost as that of an explicit scheme. It is derived from a discretization used for diffusive problems that we further enhance to efficiently solve flow problems with variable mobility and viscosity. Moreover, we solve the Navier–Stokes equations with a robust time-discretization of the projection method that guarantees better stability properties than those for Crank–Nicolson-based projection methods. For channel geometries, the method uses a spectral discretization in the streamwise and spanwise directions and a combination of spectral and high order compact finite difference discretizations in the wall normal direction. The capabilities of the method are demonstrated with several examples including phase separation with, and without, shear in two and three dimensions. The method effectively resolves interfacial layers of as few as three mesh points. The numerical examples show agreement with analytical solutions and scaling laws, where available, and the 3D simulations, in the presence of shear, reveal rich and complex structures, including strings.     

XTOR-2F: A fully implicit Newton–Krylov solver applied to nonlinear 3D extended MHD in tokamaks

XTOR-2F solves a set of extended magnetohydrodynamic (MHD) equations in toroidal tokamak geometry. In the original XTOR code, the time stepping is handled by a semi-implicit method 1, 2 and 3. Moderate changes were necessary to transform it into a fully implicit one using the NITSOL library with Newton–Krylov methods of solution for nonlinear system of equations [4]. After addressing the sensitive issue of preconditioning and time step tuning, the performances of the semi-implicit and the implicit methods are compared for the nonlinear simulation of an internal kink mode test case within the framework of resistive MHD including anisotropic thermal transport. A convergence study comparing the semi-implicit and the implicit schemes is presented. Our main conclusion is that on one hand the Newton–Krylov implicit method, when applied to basic one fluid MHD is more computationally costly than the semi-implicit one by a factor 3 for a given numerical accuracy. But on the other hand, the implicit method allows to address challenging issues beyond MHD. By testing the Newton–Krylov method with diamagnetic modifications on the dynamics of the internal kink, some numerical issues, to be addressed further, are emphasized.     

Bifurcation to 3D helical magnetic equilibrium in an axisymmetric toroidal device

We report the first direct measurement of the internal magnetic field structure associated with a 3D helical equilibrium generated spontaneously in the core of an axisymmetric toroidal plasma containment device. Magnetohydrodynamic equilibrium bifurcation occurs in a reversed-field pinch when the innermost resonant magnetic perturbation grows to a large amplitude, reaching up to 8% of the mean field strength. Magnetic topology evolution is determined by measuring the Faraday effect, revealing that, as the perturbation grows, toroidal symmetry is broken and a helical equilibrium is established.     

Dynamic simulation of locally inextensible vesicles suspended in an arbitrary two-dimensional domain,a boundary integral method

We consider numerical algorithms for the simulation of hydrodynamics of two-dimensional vesicles suspended in a viscous Stokesian fluid. The motion of vesicles is governed by the interplay between hydrodynamic and elastic forces. Continuum models of vesicles use a two-phase fluid system with interfacial forces that include tension (to maintain local “surface” inextensibility) and bending. Vesicle flows are challenging to simulate. On the one hand, explicit time-stepping schemes suffer from a severe stability constraint due to the stiffness related to high-order spatial derivatives in the bending term. On the other hand, implicit time-stepping schemes can be expensive because they require the solution of a set of nonlinear equations at each time step.     

A plasma resistive diffusion model

A self-consistent study of the slow resistive evolution of an axisymmetric toroidal plasma gives rise to a set of transport equations involving one space variable which require input from the solution of a generalized differential equation obtained from the time-differentiated Grad-Shafranov equation. An iterative scheme is presented for the numerical solution of this generalized differential equation which overcomes the problems of the non-standard boundary conditions. As an illustration this method is used to compute the instantaneous diffusion velocity of a class of model toroidal equilibria. A more detailed study is presented of the time evolution of this model in the cylindrical limit in order to illustrate techniques which can be used in a more complete toroidal simulation.     

再入体头部粘性流场高分辨率格式数值模拟

TVD格式是目前数值研究以激波为主要特征之一的超声速、高超声速流场的最先进的算法之一。本文用二阶迎风TVD格式,对三种烧蚀外形的轴对称粘性流场和10°钝锥有攻角三维粘性流场进行了数值模拟,得到了高质量的头部脱体激波和与实验结果及直线推进法计算一致的物面压力分布,表明了TVD格式在再入体粘性绕流计算中的独特优势。     

A multislope MUSCL method on unstructured meshes applied to compressible Euler equations for axisymmetric swirling flows

A finite volume method for the numerical solution of axisymmetric inviscid swirling flows is presented. The governing equations of the flow are the axisymmetric compressible Euler equations including swirl (or tangential) velocity. A first-order scheme is introduced where the convective fluxes at cell interfaces are evaluated by the Rusanov or the HLLC numerical flux while the geometric source terms are discretizated to provide a well-balanced scheme i.e. the steady-state solutions with null velocity are preserved. Extension to the second-order space approximation using a multislope MUSCL method is then derived. To test the numerical scheme, a stationary solution of the fluid flow following the radial direction has been established with a zero and nonzero tangential velocity. Numerical and exact solutions are compared for classical Riemann problems where we employ different limiters and effectiveness of the multislope MUSCL scheme is demonstrated for strongly shocked axially symmetric flows like in spherical bubble compression problem. Two other tests with axisymmetric geometries are performed: the supersonic flow in a tube with a cone and the axisymmetric blunt body with a free stream.     

Three-dimensional MHD flow over a shrinking sheet: Analytical solution and stability analysis

The magnetohydrodynamic(MHD) steady and unsteady axisymmetric flows of a viscous fluid over a two-dimensional shrinking sheet are addressed. The mathematical analysis is carried out in the presence of a large magnetic field. The steady state problem results in a singular perturbation problem having an infinite domain singularity. The secular term appearing in the solution is removed and a two-term uniformly valid solution is derived using the Lindstedt–Poincaré technique. This asymptotic solution is validated by comparing it with the numerical solution. The solution for the unsteady problem is also presented analytically in the asymptotic limit of large magnetic field. The results of velocity profile and skin friction are shown graphically to explore the physical features of the flow field. The stability analysis of the unsteady flow is made to validate the asymptotic solution.     

Adaptive solution techniques for simulating underwater explosions and implosions

Adaptive solution techniques are presented for simulating underwater explosions and implosions. The liquid is assumed to be an adiabatic fluid and the solution in the gas is assumed to be uniform in space. The solution in water is integrated in time using a semi-implicit time discretization of the adiabatic Euler equations. Results are presented either using a non-conservative semi-implicit algorithm or a conservative semi-implicit algorithm. A semi-implicit algorithm allows one to compute with relatively large time steps compared to an explicit method. The interface solver is based on the coupled level set and volume-of-fluid method (CLSVOF) [M. Sussman, A second order coupled level set and volume-of-fluid method for computing growth and collapse of vapor bubbles, J. Comput. Phys. 187 (2003) 110–136; M. Sussman, E.G. Puckett, A coupled level set and volume-of-fluid method for computing 3D and axisymmetric incompressible two-phase flows, J. Comput. Phys. 162 (2000) 301–337]. Several underwater explosion and implosion test cases are presented to show the performances of our proposed techniques.     

Shape instabilities in vesicles: A phase-field model

A phase field model for dealing with shape instabilities in fluid membrane vesicles is presented. This model takes into account the Canham-Helfrich bending energy with spontaneous curvature. A dynamic equation for the phase-field is also derived. With this model it is possible to see the vesicle shape deformation dynamically, when some external agent instabilizes the membrane, for instance, inducing an inhomogeneous spontaneous curvature. The numerical scheme used is detailed and some stationary shapes are shown together with a shape diagram for vesicles of spherical topology and no spontaneous curvature, in agreement with known results.     

Numerical simulation of 3D bubbles rising in viscous liquids using a front tracking method

The rise of bubbles in viscous liquids is not only a very common process in many industrial applications, but also an important fundamental problem in fluid physics. An improved numerical algorithm based on the front tracking method, originally proposed by Tryggvason and his co-workers, has been validated against experiments over a wide range of intermediate Reynolds and Bond numbers using an axisymmetric model [J. Hua, J. Lou, Numerical simulation of bubble rising in viscous liquid, J. Comput. Phys. 22 (2007) 769–795]. In the current paper, this numerical algorithm is further extended to simulate 3D bubbles rising in viscous liquids with high Reynolds and Bond numbers and with large density and viscosity ratios representative of the common air–water two-phase flow system. To facilitate the 3D front tracking simulation, mesh adaptation is implemented for both the front mesh on the bubble surface and the background mesh. On the latter mesh, the governing Navier–Stokes equations for incompressible, Newtonian flow are solved in a moving reference frame attached to the rising bubble. Specifically, the equations are solved using a finite volume scheme based on the Semi-Implicit Method for Pressure-Linked Equations (SIMPLE) algorithm, and it appears to be robust even for high Reynolds numbers and high density and viscosity ratios. The 3D bubble surface is tracked explicitly using an adaptive, unstructured triangular mesh. The numerical model is integrated with the software package PARAMESH, a block-based adaptive mesh refinement (AMR) tool developed for parallel computing. PARAMESH allows background mesh adaptation as well as the solution of the governing equations in parallel on a supercomputer. Further, Peskin distribution function is applied to interpolate the variable values between the front and the background meshes. Detailed sensitivity analysis about the numerical modeling algorithm has been performed. The current model has also been applied to simulate a number of cases of 3D gas bubbles rising in viscous liquids, e.g. air bubbles rising in water. Simulation results are compared with experimental observations both in aspect of terminal bubble shapes and terminal bubble velocities. In addition, we applied this model to simulate the interaction between two bubbles rising in a liquid, which illustrated the model’s capability in predicting the interaction dynamics of rising bubbles.     

II. Recirculation flow induced by a bubble stream rising in a viscous liquid

The recirculation flow induced by the rising motion of a bubble stream in a viscous fluid within an open-top rectangular enclosure is studied. The three-dimensional volume averaged conservation equations are solved by a control-volume method using a hybrid finite differencing scheme to describe the liquid phase hydrodynamics. The momentum exhange between the bubbles and the liquid phase is modeled with a source term equals to the volumetric buoyancy force acting on the gas in the bubble stream. The volumetric buoyancy force accounts for in line interactions between bubbles through the average gas volume fraction in the gas liquid column which depends on the size and the rising velocity of bubbles. The fluid flow within an open-top rectangular enclosure is further investigated by particle image velocimetry for a bubble stream rising in a water-glycerol solution. The measured fluid velocities in a vertical plane are compared with the predictions of the numerical model over a wide range of fluid viscosity (43 mPa s-800 mPa s) and gas flow rates. Finally, the recirculation flows resulting from the interaction of two neighbouring vertical bubble streams are studied. Received: 23 July 1997 / Revised: 19 December 1997 / Accepted: 11 May 1998     

An unconditionally stable rotational velocity-correction scheme for incompressible flows

We present an unconditionally stable splitting scheme for incompressible Navier–Stokes equations based on the rotational velocity-correction formulation. The main advantages of the scheme are: (i) it allows the use of time step sizes considerably larger than the widely-used semi-implicit type schemes: the time step size is only constrained by accuracy; (ii) it does not require the velocity and pressure approximation spaces to satisfy the usual inf–sup condition: in particular, the equal-order finite element/spectral element approximation spaces can be used; (iii) it only requires solving a pressure Poisson equation and a linear convection–diffusion equation at each time step. Numerical tests indicate that the computational cost of the new scheme for each time step, under identical time step sizes, is even less expensive than the semi-implicit scheme with low element orders. Therefore, the total computational cost of the new scheme can be significantly less than the usual semi-implicit scheme.     

Fluid membranes in hydrodynamic flow fields: Formalism and an application to fluctuating quasispherical vesicles in shear flow

The dynamics of a single fluid bilayer membrane in an external hydrodynamic flow field is considered. The deterministic equation of motion for the configuration is derived taking into account both viscous dissipation in the surrounding liquid and local incompressibility of the membrane. For quasi-spherical vesicles in shear flow, thermal fluctuations can be incorporated in a Langevin-type equation of motion for the deformation amplitudes. The solution to this equation shows an overdamped oscillatory approach to a stationary tanktreading shape. Inclination angle and ellipticity of the contour are determined as a function of excess area and shear rate. Comparisons to numerical results and experiments are discussed. Received 20 August 1998     

An improved high-order scheme for DNS of low Mach number turbulent reacting flows based on stiff chemistry solver

We present an improved numerical scheme for numerical simulations of low Mach number turbulent reacting flows with detailed chemistry and transport. The method is based on a semi-implicit operator-splitting scheme with a stiff solver for integration of the chemical kinetic rates, developed by Knio et al. [O.M. Knio, H.N. Najm, P.S. Wyckoff, A semi-implicit numerical scheme for reacting flow II. Stiff, operator-split formulation, Journal of Computational Physics 154 (2) (1999) 428–467]. Using the material derivative form of continuity equation, we enhance the scheme to allow for large density ratio in the flow field. The scheme is developed for direct numerical simulation of turbulent reacting flow by employing high-order discretization for the spatial terms. The accuracy of the scheme in space and time is verified by examining the grid/time-step dependency on one-dimensional benchmark cases: a freely propagating premixed flame in an open environment and in an enclosure related to spark-ignition engines. The scheme is then examined in simulations of a two-dimensional laminar flame/vortex-pair interaction. Furthermore, we apply the scheme to direct numerical simulation of a homogeneous charge compression ignition (HCCI) process in an enclosure studied previously in the literature. Satisfactory agreement is found in terms of the overall ignition behavior, local reaction zone structures and statistical quantities. Finally, the scheme is used to study the development of intrinsic flame instabilities in a lean H₂/air premixed flame, where it is shown that the spatial and temporary accuracies of numerical schemes can have great impact on the prediction of the sensitive nonlinear evolution process of flame instability.     

粘性不可压流体流动问题用直角坐标网格的贴体解法

研究一种新的全贴体的求解粘性不可压流体流动问题的非结构化直角坐标网格方法.该方法在于利用直角坐标网格但通过在边界附近保留不规则控制体,使得算法是完全贴体的.这有别于目前流行的各种非结构化直角坐标网格方法.通过对两个典型流动问题的计算对该数值方法进行验证.对比结果表明,本方法计算的结果与精确解和STAR-CD的结果在一定Re数和网格数时是很接近的,可以满足一定的精度要求,说明该数值计算方法是可行的.还对二维钝头体周围的流场进行了计算,计算的流场与STAR-CD的结果相当吻和,说明该算法还可计算较复杂的流动现象.