Complex Computational Domains Research Articles

NeuroBox: Computational Mathematics in Multiscale Neuroscience.

The brain is a complex organ operating on multiple scales. From molecular events that inform electrical and biochemical cellular responses, the brain interconnects processes all the way up to the massive network size of billions of brain cells. This strongly coupled, nonlinear, system has been subject to research that has turned increasingly multidisciplinary. The seminal work of Hodgkin and Huxley in the 1950s made use of experimental data to derive a coherent physical model of electrical signaling in neurons, which can be solved using mathematical and computational methods, thus bringing together neuroscience, physics, mathematics, and computer science. Over the last decades numerous projects have been dedicated to modeling and simulation of specific parts of molecular dynamics, neuronal signaling, and neural network behavior. Simulators have been developed around a specific objective and scale, in order to cope with the underlying computational complexity. Often times a dimension reduction approach allows larger scale simulations, this however has the inherent drawback of losing insight into structure-function interplay at the cellular level. This paper gives an overview of the project NeuroBox that has the objective of integrating multiple brain scales and associated physical models into one unified framework. NeuroBox hosts geometry and anatomical reconstruction methods, such that detailed three-dimensional domains can be integrated into numerical simulations of models based on partial differential equations. The project further focusses on deriving numerical methods for handling complex computational domains, and to couple multiple spatial dimensions. The latter allows the user to specify in which parts of the biological problem high-dimensional representations are necessary and where low-dimensional approximations are acceptable. NeuroBox offers workflow user interfaces that are automatically generated with VRL-Studio and can be controlled by non-experts. The project further uses uG4 as the numerical backend, and therefore accesses highly advanced discretization methods as well as hierarchical and scalable numerical solvers for very large neurobiological problems.

A simple diffuse interface approach on adaptive Cartesian grids for the linear elastic wave equations with complex topography

In most classical approaches of computational geophysics for seismic wave propagation problems, complex surface topography is either accounted for by boundary-fitted unstructured meshes, or, where possible, by mapping the complex computational domain from physical space to a topologically simple domain in a reference coordinate system. However, all these conventional approaches face problems if the geometry of the problem becomes sufficiently complex. They either need a mesh generator to create unstructured boundary-fitted grids, which can become quite difficult and may require a lot of manual user interactions in order to obtain a high quality mesh, or they need the explicit computation of an appropriate mapping function from physical to reference coordinates. For sufficiently complex geometries such mappings may either not exist or their Jacobian could become close to singular. Furthermore, in both conventional approaches low quality grids will always lead to very small time steps due to the Courant-Friedrichs-Lewy (CFL) condition for explicit schemes. In this paper, we propose a completely different strategy that follows the ideas of the successful family of high resolution shock-capturing schemes, where discontinuities can actually be resolved anywhere on the grid, without having to fit them exactly. We address the problem of geometrically complex free surface boundary conditions for seismic wave propagation problems with a novel diffuse interface method (DIM) on adaptive Cartesian meshes (AMR) that consists in the introduction of a characteristic function 0≤α≤1 which identifies the location of the solid medium and the surrounding air (or vacuum) and thus implicitly defines the location of the free surface boundary. Physically, α represents the volume fraction of the solid medium present in a control volume. Our new approach completely avoids the problem of mesh generation, since all that is needed for the definition of the complex surface topography is to set a scalar color function to unity inside the regions covered by the solid and to zero outside. The governing equations are derived from ideas typically used in the mathematical description of compressible multiphase flows. An analysis of the eigenvalues of the PDE system shows that the complexity of the geometry has no influence on the admissible time step size due to the CFL condition. The model reduces to the classical linear elasticity equations inside the solid medium where the gradients of α are zero, while in the diffuse interface zone at the free surface boundary the governing PDE system becomes nonlinear. We can prove that the solution of the Riemann problem with arbitrary data and a jump in α from unity to zero yields a Godunov-state at the interface that satisfies the free-surface boundary condition exactly, i.e. the normal stress components vanish. In the general case of an interface that is not aligned with the grid and which is not infinitely thin, a systematic study on the distribution of the volume fraction function inside the interface and the sensitivity with respect to the thickness of the diffuse interface layer has been carried out. In order to reduce numerical dissipation, we use high order discontinuous Galerkin (DG) finite element schemes on adaptive AMR grids together with a second order accurate high resolution shock capturing subcell finite volume (FV) limiter in the diffuse interface region. We furthermore employ a little dissipative HLLEM Riemann solver, which is able to resolve the steady contact discontinuity associated with the volume fraction function and the spatially variable material parameters exactly. While the method is locally high order accurate in the regions without limiter, the global order of accuracy of the scheme is at most two if the limiter is activated. It is locally of order one inside the diffuse interface region, which is typical for shock-capturing schemes at shocks and contact discontinuities. We show a large set of computational results in two and three space dimensions involving complex geometries where the physical interface is not aligned with the grid or where it is even moving. For all test cases we provide a quantitative comparison with classical approaches based on boundary-fitted unstructured meshes.

МЕТОД ЧИСЕЛЬНОГО МОДЕЛЮВАННЯ ГІДРОДИНАМІЧНИХ ПРОЦЕСІВ

This article presents a numerical simulation method that is used in the study of the dynamics of continuous viscous weakly compressed fluids based on the system of equations of continuity and Navier-Stokes. The proposed method uses a complex approach using the numerical solution of the continuity equation by the finite-volume method, and for solving the Navier-Stokes equation the splitting method by physical factors. The article shows that the finite volume method, which was used to describe the flow of both compressed and uncompressed liquids, has such important advantages as the presence of good conservative properties and assumptions of discretization of complex computational domains into simpler ones than isoparametric finite element formulation of the problem or the introduction of generalized coordinates. A component is introduced into the method of splitting according to physical factors, which takes into account the artificial compressibility of the test liquid, which allows you to first calculate the intermediate velocity field, which is then corrected taking into account the pressure gradients. The difference scheme of this method allows one to calculate the flow field without using the values ​​of the vortex and pressure on a solid surface. In the framework of the proposed approach, it is not necessary to calculate the value of the vortex on a solid surface. The latter can be found from the calculated velocity field using any of the difference representations of the expression for the vortex at the boundary points. To confirm the effectiveness of the proposed method, in the FlowVision CFD program, solutions were obtained for a number of problems in external hydrodynamics, using the example of a cylindrical surface optics, which confirmed the stability of the results obtained. This method allows one to compute the flow of flat, axisymmetric, and three-dimensional bodies of a complex configuration with a flow of a viscous weakly compressed fluid, as well as internal flows in a wide range of Reynolds numbers using a single algorithm.

Development of a three‐dimensional solver based on the local domain‐free discretization and immersed boundary method and its application for incompressible flow problems

SummaryRecently, the author and two other coauthors have proposed a two‐dimensional hybrid local domain‐free discretization and immersed boundary method (LDFD‐IBM), which can be used to solve the flow problem with complex geometries. In this paper, the LDFD‐IBM is extended to solve a three‐dimensional unsteady incompressible flow with the complex computational domain. The technical issues related to the implementation of the LDFD‐IBM in three‐dimensional problems are discussed in detail, particularly for the discretization of Navier‐Stokes equations, mesh strategies for a three‐dimensional flow, and the fast algorithm on the identification of the status of mesh nodes (ie, to identify if the mesh node is located in the solid domain, in the fluid domain, or near the immersed boundary). Numerical tests show that the LDFD‐IBM can accurately solve three‐dimensional incompressible problems with ease.

Boundedness-preserving implicit correction of mesh-induced errors for VOF based heat and mass transfer

Spatial discretisation of geometrically complex computational domains often entails unstructured meshes of general topology for Computational Fluid Dynamics (CFD). Mesh skewness is then typically encountered causing severe deterioration of the formal order of accuracy of the discretisation, or boundedness of the solution, or both. Particularly methods inherently relying on the accurate and bounded transport of sharp fields suffer from all types of mesh-induced skewness errors, namely both non-orthogonality and non-conjunctionality errors.This work is devoted to a boundedness-preserving strategy to correct for skewness errors arising from discretisation of advection and diffusion terms within the context of interfacial heat and mass transfer based on the Volume-of-Fluid methodology. The implementation has been accomplished using a second-order finite volume method with support for unstructured meshes of general topology. We examine and advance suitable corrections for the finite volume discretisation of a consistent single-field model, where both accurate and bounded transport due to diffusion and advection is crucial. In order to ensure consistency of both the volume fraction and the species concentration transport, i.e. to avoid artificial heat or species transfer, corrections are studied for both cases using distorted 2D meshes. The cross interfacial jump and adjacent sharp gradients of species concentration render the correction for skewness-induced diffusion and advection errors additionally demanding and has not so far been addressed in the literature.

Parametric finite elements with bijective mappings

The discretization of the computational domain plays a central role in the finite element method. In the standard discretization the domain is triangulated with a mesh and its boundary is approximated by a polygon. The boundary approximation induces a geometry-related error which influences the accuracy of the solution. To control this geometry-related error, iso-parametric finite elements and iso-geometric analysis allow for high order approximation of smooth boundary features. We present an alternative approach which combines parametric finite elements with smooth bijective mappings leaving the choice of approximation spaces free. Our approach allows to represent arbitrarily complex geometries on coarse meshes with curved edges, regardless of the domain boundary complexity. The main idea is to use a bijective mapping for automatically warping the volume of a simple parameterization domain to the complex computational domain, thus creating a curved mesh of the latter. Numerical examples provide evidence that our method has lower approximation error for domains with complex shapes than the standard finite element method, because we are able to solve the problem directly on the exact domain without having to approximate it.

Saturated-unsaturated groundwater modeling using 3D Richards equation with a coordinate transform of nonorthogonal grids

• Coordinate transform of nonorthogonal grids for complex computational domain. • Derivation of Richards equation in arbitrary coordinates and boundary conditions. • Saturated-unsaturated groundwater modeling using 3D Richards equation • Four test bench simulations: hillslope, coupled flow, tilted V-Catchment and slab. Saturated-unsaturated flow under a complex terrain is usually solved using the Richards equation. Finite difference or finite volume methods are commonly employed for discretization of Richards equation because of simplicity of coding. Complex subsurface boundary geometries lead to nonorthogonal grids in curvilinear grid systems, which leads to difficulty in discretization and mesh generation. This paper develops a vertical coordinate transform, enabling a computational domain regular in the vertical direction. As a result, the grid of curvilinear surfaces can be successfully transformed to a computational grid that allows solution of the Richards equation with efficient computation and simpler coding. The anisotropic Richards equation in the Cartesian coordinate system is transformed to the equation in the arbitrary coordinate system and further expressed as a form appropriate to the orthogonal coordinate system. The generalized third boundary condition is transformed to a form suited to the orthogonal coordinate system. The finite volume method is used to solve the Richards equation in the orthogonal coordinate system. Four cases are used to validate the present orthogonal coordinate system. The computational results from the orthogonal coordinate system are in good agreement with the results from Ansys Fluent solved in a Cartesian coordinate system for the subsurface flow case. For the coupled case of hill slopes, a good agreement between the computational results and the experimental data is obtained. The present results for V-titled catchment and slab case accord well with the results obtained from HydroGeoSphere and PAWS. The present algorithm can improve grid generation for solution of Richards equation in a hydrological model for a complex domain.

Generating Neuron Geometries for Detailed Three-Dimensional Simulations Using AnaMorph.

Generating realistic and complex computational domains for numerical simulations is often a challenging task. In neuroscientific research, more and more one-dimensional morphology data is becoming publicly available through databases. This data, however, only contains point and diameter information not suitable for detailed three-dimensional simulations. In this paper, we present a novel framework, AnaMorph, that automatically generates water-tight surface meshes from one-dimensional point-diameter files. These surface triangulations can be used to simulate the electrical and biochemical behavior of the underlying cell. In addition to morphology generation, AnaMorph also performs quality control of the semi-automatically reconstructed cells coming from anatomical reconstructions. This toolset allows an extension from the classical dimension-reduced modeling and simulation of cellular processes to a full three-dimensional and morphology-including method, leading to novel structure-function interplay studies in the medical field. The developed numerical methods can further be employed in other areas where complex geometries are an essential component of numerical simulations.

Phenomena of Nonlinear Diffusion in Complex 3D Media

Unbounded solutions (critical blow-up regimes) simulated by the 3D nonlinear diffusion equation in a spherical shell are studied. The coordinate splitting of the differential operator coupled with two spherical coordinate maps makes it possible to use periodic boundary conditions in the latitudinal and longitudinal directions and employ the computationally efficient Sherman-Morrison formula and Thomas algorithm. The resulting finite difference method is direct, with implicit and unconditionally stable schemes of second-order approximation in all the variables. Numerical tests demonstrate that it allows simulating different blow-up regimes in complex computational domains.

A modified direct method for void fraction calculation in CFD–DEM simulations

The void fraction of computational cells in numerical simulations of particulate flows using computational fluid dynamics–discrete element method (CFD–DEM) is often directly (or crudely) calculated assuming that the entire body of a particle lies in the cell at which the particle centroid resides. This direct method is most inexpensive but inaccurate and may lead to simulation instabilities. In this study, a modified version of the direct method has been proposed. In this method, referred to as the particle meshing method (PMM), the particle is meshed and the solid volume in a fluid cell is calculated by adding up the particle mesh volume with the basic working principle being the same as that of the direct method. As a result, the PMM inherits the simplicity and hence the computational advantage from the direct method, whilst allowing for duplicating the particle shape and accurate accounting of particle volume in each fluid cell. The numerical simulation characteristics of PMM including numerical stability, minimum particle grid number, prediction accuracy, and computational efficiency have been examined. The results showed that for a specific cell-to-particle size ratio, there was a minimum particle grid number required to reach the stable simulation. A formula of estimating the minimum particle grid number was derived and discussed. Typically, a particle grid number of about 5 times the minimum number was suggested to achieve the best computational efficiency, which was comparable or even higher than that of simulations using the analytical approach. PMM also exhibited the potential to be applied for complex computational domain geometries and irregular shaped particles.

High-order finite volume shallow water model on the cubed-sphere: 1D reconstruction scheme

A central-upwind finite-volume (CUFV) scheme for shallow-water model on a nonorthogonal equiangular cubed-sphere grid is developed, consequently extending the 1D reconstruction CUFV transport scheme developed by us. High-order spatial discretization based on weighted essentially non-oscillatory (WENO) is considered for this effort. The CUFV method combines the alluring features of classical upwind and central schemes. This approach is particularly useful for complex computational domain such as the cubed-sphere. The continuous flux-form spherical shallow water equations in nonorthogonal curvilinear coordinates are utilized. Fluxes along the element boundaries are approximated by a Kurganov–Noelle–Petrova scheme. A fourth-order strong stability preserving Runge–Kutta time stepping scheme for time integration is employed in the present work. The numerical scheme is evaluated with standard shallow water test suite, which accentuate accuracy and conservation properties. In addition, an efficient yet inexpensive bound preserving filter with an optional positivity filter is used to remove spurious oscillations and to achieve strictly positive definite numerical solution. To tackle the discontinuities that arise at the edges of the cubed-sphere grid, we utilize a high-order 1D interpolation procedure combining cubic and quadratic interpolations.The results with the high-order scheme is compared with the results for the same tests for various schemes available in literature. Since, the scheme presented here uses local-cell information, it is expected to be scalable to high number of processors count in a distributed node high-performance computer.

An interface‐enriched generalized finite element analysis for electromagnetic problems with non‐conformal discretizations

SummaryAn interface‐enriched generalized finite element method is presented for analyzing electromagnetic problems involving highly inhomogeneous materials. To avoid creating conformal meshes within a complex computational domain and preparing multiple meshes during optimization, enriched vector basis functions are introduced over the finite elements that intersect the material interfaces to capture the normal derivative discontinuity of the tangential field component. These enrichment functions are directly constructed from a linear combination of the vector basis functions of the sub‐elements. Several numerical examples are presented to verify the method with analytical solutions and demonstrate its h‐refinement convergence rate. The proposed interface‐enriched generalized finite element method is shown to achieve the same level of accuracy as the standard finite element method based on conformal meshes. Two examples, involving multiple microvascular channels and circular inclusions of different radii, are analyzed to illustrate the capability of the proposed approach in handling complicated inhomogeneous geometries. Copyright © 2015 John Wiley & Sons, Ltd.

Advection step in the split-operator approach applied to river modeling

The split-operator approach proposes to divide the numerical solution of the governing flow and constituent transport equations into a series of successive steps, one of them being the numerically most challenging advection step. This paper presents tests of characteristics and finite difference methods used in the numerical solution of one-dimensional linear and nonlinear advection equations that describe the transport and flow advection steps, thus demonstrating their accuracy, oscillations and damping features. The studied methods were incorporated into the split-operator approach and evaluated as a part of the complete solution, which led to the development of numerical models implementing an algorithm that can endure larger time steps without losing much accuracy. This algorithm was tested on both one- and two-dimensional examples, for complex computational domains that predominate in natural watercourses.

Efficient FDTD-PML Simulation of Gain Medium Based on Exponential Time Differencing Algorithm

An exponential time differencing (ETD) algorithm is introduced to incorporate the homogeneously broadened Lorentzian oscillator of gain medium in a four-level atomic system into the finite difference time domain (FDTD). Compared with the well known auxiliary differential equation (ADE) method, the proposed algorithm shows the same accuracy but can save one-third of the additional memory space for treating the Lorentzian oscillators, and has simpler formulations. The ETD implementation of the stretched coordinates perfectly matched layer (SC-PML) with the complex frequency shifted (CFS) stretching variable is applied to truncate computational domain with gain media. Simulations involving both TE and TM waves indicate that compared with the modified conventional PMLs and the convolution PML (CPML), the proposed absorbing boundary formulations can lead to a significant improvement of the absorbing performance with a similar memory requirement. Compared with the application of the material dependent PMLs to the gain media in previous studies, the proposed absorbing boundary condition has simpler formulas and can be applied to complex computational domains much more straightforwardly.

A 3D hybrid grid generation technique and a multigrid/parallel algorithm based on anisotropic agglomeration approach

A hybrid grid generation technique and a multigrid/parallel algorithm are presented in this paper for turbulence flow simulations over three-dimensional (3D) complex geometries. The hybrid grid generation technique is based on an agglomeration method of anisotropic tetrahedrons. Firstly, the complex computational domain is covered by pure tetrahedral grids, in which anisotropic tetrahedrons are adopted to discrete the boundary layer and isotropic tetrahedrons in the outer field. Then, the anisotropic tetrahedrons in the boundary layer are agglomerated to generate prismatic grids. The agglomeration method can improve the grid quality in boundary layer and reduce the grid quantity to enhance the numerical accuracy and efficiency. In order to accelerate the convergence history, a multigrid/parallel algorithm is developed also based on anisotropic agglomeration approach. The numerical results demonstrate the excellent accelerating capability of this multigrid method.

A new triangular spectral element method I: implementation and analysis on a triangle

This paper serves as our first effort to develop a new triangular spectral element method (TSEM) on unstructured meshes, using the rectangle---triangle mapping proposed in the conference note (Li et al. 2011). Here, we provide some new insights into the originality and distinctive features of the mapping, and show that this transform only induces a logarithmic singularity, which allows us to devise a fast, stable and accurate numerical algorithm for its removal. Consequently, any triangular element can be treated as efficiently as a quadrilateral element, which affords a great flexibility in handling complex computational domains. Benefited from the fact that the image of the mapping includes the polynomial space as a subset, we are able to obtain optimal L 2- and H 1-estimates of approximation by the proposed basis functions on triangle. The implementation details and some numerical examples are provided to validate the efficiency and accuracy of the proposed method. All these will pave the way for developing an unstructured TSEM based on, e.g., the hybridizable discontinuous Galerkin formulation.

Solution of the equation of radiative transfer using a Newton–Krylov approach and adaptive mesh refinement

The discrete ordinates method (DOM) and finite-volume method (FVM) are used extensively to solve the radiative transfer equation (RTE) in furnaces and combusting mixtures due to their balance between numerical efficiency and accuracy. These methods produce a system of coupled partial differential equations which are typically solved using space-marching techniques since they converge rapidly for constant coefficient spatial discretization schemes and non-scattering media. However, space-marching methods lose their effectiveness when applied to scattering media because the intensities in different directions become tightly coupled. When these methods are used in combination with high-resolution limited total-variation-diminishing (TVD) schemes, the additional non-linearities introduced by the flux limiting process can result in excessive iterations for most cases or even convergence failure for scattering media. Space-marching techniques may also not be quite as well-suited for the solution of problems involving complex three-dimensional geometries and/or for use in highly-scalable parallel algorithms. A novel pseudo-time marching algorithm is therefore proposed herein to solve the DOM or FVM equations on multi-block body-fitted meshes using a highly scalable parallel-implicit solution approach in conjunction with high-resolution TVD spatial discretization. Adaptive mesh refinement (AMR) is also employed to properly capture disparate solution scales with a reduced number of grid points. The scheme is assessed in terms of discontinuity-capturing capabilities, spatial and angular solution accuracy, scalability, and serial performance through comparisons to other commonly employed solution techniques. The proposed algorithm is shown to possess excellent parallel scaling characteristics and can be readily applied to problems involving complex geometries. In particular, greater than 85% parallel efficiency is demonstrated for a strong scaling problem on up to 256 processors. Furthermore, a speedup of a factor of at least two was observed over a standard space-marching algorithm using a limited scheme for optically thick scattering media. Although the time-marching approach is approximately four times slower for absorbing media, it vastly outperforms standard solvers when parallel speedup is taken into account. The latter is particularly true for geometrically complex computational domains.

Simulating free surface flow problems using hybrid particle element free Galerkin method

A hybrid particle element free Galerkin (HPEFG) method is proposed for simulating free surface flow problems. The method combines the advantages of the element free Galerkin (EFG) method in the Eulerian description and particle method in the Lagrangian description, while overcoming the disadvantages of each. It makes the EFG method possible for solving problems of fluid flow with free surfaces and overcomes the difficulties in dealing with complex computational domains and mesh distortions. Two sets of nodes are independently placed in the computational domain for the HPEFG method. The nodes in one set used to represent fluid are called particles, and those in the other set used to compute flow variables are called EFG nodes. In order to test the validity and accuracy of the proposed meshfree method, four free surface flow problems are simulated and compared with other numerical results obtained by other numerical techniques, including the dam breaking, fluid–solid/fluid impacting and rotation of a square fluid patch. All the results show feasibility and ability of the proposed method to solve incompressible viscous fluid dynamic problems with moving free surface boundaries.

Numerical Simulation of Shallow-Water Flow Using a Modified Cartesian Cut-Cell Approach

The Cartesian cut-cell method can be used to represent irregular and complex computational domains with less computational efforts by cutting the grid cells on the boundary surfaces in a background uniform Cartesian mesh. In this study, a modified Cartesian cut-cell grid technique is proposed to better represent complex physical geometries. A point shifting treatment was employed to determine the start and end points of a line segment in cut-cell grids. This led to an improved representation of sharply-shaped corners in surface polygons. Numerical simulation to solve a set of shallow-water equations was performed by incorporating a finite volume approach into the Cartesian cut-cell mesh. The advective fluxes at intercells were first estimated by a Harten, Lax and van Leer for contact wave approximate Riemann solver. In order to improve the model accuracy to the second order, a total variation diminishing-weighted average flux method was applied to work adaptively with the cut-cell mesh. The numerical model was then employed to simulate dam-break flow propagation in a small channel with a rectangular obstacle or a 45° bend. The numerical results show good agreement with available laboratory measurements.

A finite element algorithm for particle/droplet trajectory tracking, tested in a liquid–liquid system in the presence of an external electric field

In this paper we are concerned with modeling a dispersion of electrically charged droplets in motion in a second immiscible liquid phase continuum in the presence of an external electric field. The system is highly relevant to solvent extraction processes and to liquid–liquid reactions such as phase transfer catalysis. Electrostatic enhancement of liquid–liquid contacting processes is well known and relies upon intensified drop breakup and drop acceleration due to the presence of significant electrical forces. The enhanced drop motion results in greatly increased rates of mass transfer. The scope of the current work is the mathematical modeling and experimental validation of the trajectories of multiple charged drops in three-dimensional motion in a liquid–liquid system in the presence of an externally applied field. A new particle tracking algorithm using a tetrahedral finite element mesh to solve the relevant system of differential equations is described in detail. Such elements are shown to be especially convenient for real word simulation because of their simplicity, their high stability with strictly linear interpolation, and their flexibility for dealing with complex computational domains. The calculations of trajectory tracking for several hundred particles taking in a mesh consisting of 50,000 tetrahedral elements was shown only to require a few seconds of computational time using a PC class computer. Using a double step stepping algorithm, with intermediate velocity interpolation, second-order convergence was proven in a rigorous benchmark test. Comparison of simulations of individual drop motion with experimental measurements proved to be very accurate. The full three-dimensional simulation, having greater requirements and involving more degrees of freedom than the earlier two-dimensional case was successfully demonstrated, with greater accuracy and economic use of computational time. Preliminary simulation of swarming motion of droplets using a cloud model is also presented. This shows significant promise for the evaluation of the mass distribution of very small droplets (∼1 to 100 μm) in a real contactor for optimal shape design. The model, when coupled with the models for mass transfer kinetics and interfacial reaction kinetics, should provide a very valuable tool for liquid–liquid reactor design.