We formulate a one-dimensional time-dependent non-linear mathematical model for some types of physiological fluid flow in collapsible tubes with discontinuous material properties. The resulting 6 × 6 hyperbolic system is analysed and the associated Riemann problem is solved exactly. Although the solution algorithm deals with idealised cases, it is nonetheless uniquely well-suited for assessing the performance of numerical methods intended for simulating more general situations. Moreover, our model may be a useful starting point for numerical calculations of realistic flows involving rapid and discontinuous material property variations. One important example in mind is the simulation of blood flow in medium-to-large veins in humans. Finally, we also discuss some peculiarities of the model regarding the loss of strict hyperbolicity and uniqueness. In particular we show an example in which the solution of the Riemann problem is non unique.

The lattice Boltzmann equation (LBE) methods (both LBGK and MRT) and the discrete unified gas-kinetic scheme (DUGKS) are both derived from the Boltzmann equation, but with different consideration in their algorithm construction. With the same numerical discretization in the particle velocity space, the distinctive modeling of these methods in the update of gas distribution function may introduce differences in the computational results. In order to quantitatively evaluate the performance of these methods in terms of accuracy, stability, and efficiency, in this paper we test LBGK, MRT, and DUGKS in two-dimensional cavity flow and the flow over a square cylinder, respectively. The results for both cases are validated against benchmark solutions. The numerical comparison shows that, with sufficient mesh resolution, the LBE and DUGKS methods yield qualitatively similar results in both test cases. With identical mesh resolutions in both physical and particle velocity space, the LBE methods are more efficient than the DUGKS due to the additional particle collision modeling in DUGKS. But, the DUGKS is more robust and accurate than the LBE methods in most test conditions. Particularly, for the unsteady flow over a square cylinder at Reynolds number 300, with the same mesh resolution it is surprisingly observed that the DUGKS can capture the physical multi-frequency vortex shedding phenomena while the LBGK and MRT fail to get that. Furthermore, the DUGKS is a finite volume method and its computational efficiency can be much improved when a non-uniform mesh in the physical space is adopted. The comparison in this paper clearly demonstrates the progressive improvement of the lattice Boltzmann methods from LBGK, to MRT, up to the current DUGKS, along with the inclusion of more reliable physical process in their algorithm development. Besides presenting the Navier-Stokes solution, the DUGKS can capture the rarefied flow phenomena as well with the increasing of Knudsen number.

In this article we present a new class of high order accurate Arbitrary-Eulerian-Lagrangian (ALE) one-step WENO finite volume schemes for solving nonlinear hyperbolic systems of conservation laws on moving two dimensional unstructured triangular meshes. A WENO reconstruction algorithm is used to achieve high order accuracy in space and a high order one-step time discretization is achieved by using the local space-time Galerkin predictor proposed in. For that purpose, a new element-local weak formulation of the governing PDE is adopted on moving space-time elements. The space-time basis and test functions are obtained considering Lagrange interpolation polynomials passing through a predefined set of nodes. Moreover, a polynomial mapping defined by the same local space-time basis functions as the weak solution of the PDE is used to map the moving physical space-time element onto a space-time reference element. To maintain algorithmic simplicity, the final ALE one-step finite volume scheme uses moving triangular meshes with straight edges. This is possible in the ALE framework, which allows a local mesh velocity that is different from the local fluid velocity. We present numerical convergence rates for the schemes presented in this paper up to sixth order of accuracy in space and time and show some classical numerical test problems for the two-dimensional Euler equations of compressible gas dynamics.

The commonly used incompressible phase field models for non-reactive, binary fluids, in which the Cahn-Hilliard equation is used for the transport of phase variables (volume fractions), conserve the total volume of each phase as well as the material volume, but do not conserve the mass of the fluid mixture when densities of two components are different. In this paper, we formulate the phase field theory for mixtures of two incompressible fluids, consistent with the quasi-compressible theory [28], to ensure conservation of mass and momentum for the fluid mixture in addition to conservation of volume for each fluid phase. In this formulation, the mass-average velocity is no longer divergence-free (solenoidal) when densities of two components in the mixture are not equal, making it a compressible model subject to an internal con-straint. In one formulation of the compressible models with internal constraints (model 2), energy dissipation can be clearly established. An efficient numerical method is then devised to enforce this compressible internal constraint. Numerical simulations in confined geometries for both compressible and the incompressible models are carried out using spatially high order spectral methods to contrast the model predictions. Numerical comparisons show that (a) predictions by the two models agree qualitatively in the situation where the interfacial mixing layer is thin; and (b) predictions differ significantly in binary fluid mixtures undergoing mixing with a large mixing zone. The numerical study delineates the limitation of the commonly used incompressible phase field model using volume fractions and thereby cautions its predictive value in simulating well-mixed binary fluids.

We describe our implementation of a parallel fast multipole method for evaluating potentials for discrete and continuous source distributions. The first requires summation over the source points and the second requiring integration over a continuous source density. Both problems require (N2) complexity when computed directly; however, can be accelerated to (N) time using FMM. In our PVFMM software library, we use kernel independent FMM and this allows us to compute potentials for a wide range of elliptic kernels. Our method is high order, adaptive and scalable. In this paper, we discuss several algorithmic improvements and performance optimizations including cache locality, vectorization, shared memory parallelism and use of coprocessors. Our distributed memory implementation uses space-filling curve for partitioning data and a hypercube communication scheme. We present convergence results for Laplace, Stokes and Helmholtz (low wavenumber) kernels for both particle and volume FMM. We measure efficiency of our method in terms of CPU cycles per unknown for different accuracies and different kernels. We also demonstrate scalability of our implementation up to several thousand processor cores on the Stampede platform at the Texas Advanced Computing Center.

The paper is devoted to the development of an efficient deterministic framework for modelling of three-dimensional rarefied gas flows on the basis of the numerical solution of the Boltzmann kinetic equation with the model collision integrals. The framework consists of a high-order accurate implicit advection scheme on arbitrary unstructured meshes, the conservative procedure for the calculation of the model collision integral and efficient implementation on parallel machines. The main application area of the suggested methods is micro-scale flows. Performance of the proposed approach is demonstrated on a rarefied gas flow through the finite-length circular pipe. The results show good accuracy of the proposed algorithm across all flow regimes and its high efficiency and excellent parallel scalability for up to 512 cores.

The simulation of multiphase flows is an outstanding challenge, due to the inherent complexity of the underlying physical phenomena and to the fact that multiphase flows are very diverse in nature, and so are the laws governing their dynamics. In the last two decades, a new class of mesoscopic methods, based on minimal lattice formulation of Boltzmann kinetic equation, has gained significant interest as an efficient alternative to continuum methods based on the discretisation of the NS equations for non ideal fluids. In this paper, three different multiphase models based on the lattice Boltzmann method (LBM) are discussed, in order to assess the capability of the method to deal with multiphase flows on a wide spectrum of operating conditions and multiphase phenomena. In particular, the range of application of each method is highlighted and its effectiveness is qualitatively assessed through comparison with numerical and experimental literature data.

As a continuation of developing multiscale method for the transport phenomena, a unified gas kinetic scheme (UGKS) for multi-scale and multi-component plasma simulation is constructed. The current scheme is a direct modeling method, where the time evolution solutions from the Vlasov-BGK equations of electron and ion and the Maxwell equations are used to construct a scale-dependent plasma simulation model. The modeling scale used in the UGKS is the mesh size scale, which can be comparable to or much larger than the local mean free path. As a result, with the variation of modeling scales in space and time through the so-called cell's Knudsen number and normalized Larmor radius, the discretized governing equations can recover a wide range of plasma evolution from the Vlasov equation in the kinetic scale to different-type of magnetohydrodynamic (MHD) equations in the hydrodynamic scale. The UGKS provides a general evolution model, which goes to the Vlasov equation in the kinetic scale and many types of MHD equations in the hydrodynamic scale, such as the two fluids model, the Hall, the resistive, and the ideal MHD equations. All current existing governing equations become the subsets of the UGKS, and the UGKS bridges these distinguishable governing equations seamlessly. The construction of UGKS is based on the implementation of physical conservation laws and the un-splitting treatment of particle collision, acceleration, and transport in the construction of a scale-dependent numerical flux across a cell interface. At the same time, the discretized plasma evolution equations are coupled with the Maxwell equations for electro-magnetic fields, which also cover a scale-dependent transition between the Ampére's law and the Ohm's law for the calculation of electric field. The time step of UGKS is not limited by the relaxation time, the cyclotron period, and the speed of light in the ideal-MHD regime. Our scheme is able to give a physically accurate solution for plasma simulation with a wide range of Knudsen number and normalized Larmor radius. It can be used to study the phenomena from the Vlasov limit to the scale of plasma skin depth for the capturing of two-fluid effect, and the phenomena in the plasma transition regime with a modest Knudsen number and Larmor radius. The UGKS is validated by numerical test cases, such as the Landau damping and two stream instability in the kinetic regime, and the Brio-Wu shock tube problem, and the Orszag-Tang MHD turbulence problem in the hydrodynamic regime. The scheme is also used to study the geospace environment modeling (GEM), such as the challenging magnetic reconnection problem in the transition regime. At the same time, the magnetic reconnection mechanism of the Sweet-Parker model and the Hall effect model can be connected smoothly through the variation of Larmor radius in the UGKS simulations. Overall, the UGKS is a physically reliable multi-scale plasma simulation method, and it provides a powerful and unified approach for the study of plasma physics.

The immersed boundary (IB) method is an approach to problems of fluid-structure interaction in which an elastic structure is immersed in a viscous incompressible fluid. The IB formulation of such problems uses a Lagrangian description of the structure and an Eulerian description of the fluid. It is well known that some versions of the IB method can suffer from poor volume conservation. Methods have been introduced to improve the volume-conservation properties of the IB method, but they either have been fairly specialized, or have used complex, nonstandard Eulerian finite-difference discretizations. In this paper, we use quasi-static and dynamic benchmark problems to investigate the effect of the choice of Eulerian discretization on the volume-conservation properties of a formally second-order accurate IB method. We consider both collocated and staggered-grid discretization methods. For the tests considered herein, the staggered-grid IB scheme generally yields at least a modest improvement in volume conservation when compared to cell-centered methods, and in many cases considered in this work, the spurious volume changes exhibited by the staggered-grid IB method are more than an order of magnitude smaller than those of the collocated schemes. We also compare the performance of cell-centered schemes that use either exact or approximate projection methods. We find that the volume-conservation properties of approximate projection IB methods depend strongly on the formulation of the projection method. When used with the IB method, we find that pressure-free approximate projection methods can yield extremely poor volume conservation, whereas pressure-increment approximate projection methods yield volume conservation that is nearly identical to that of a cell-centered exact projection method.

The aim of the present work is to develop a general formalism to derive staggered discretizations for Lagrangian hydrodynamics on two-dimensional unstructured grids. To this end, we make use of the compatible discretization that has been initially introduced by E. J. Caramana et al., in J. Comput. Phys., 146 (1998). Namely, momentum equation is discretized by means of subcell forces and specific internal energy equation is obtained using total energy conservation. The main contribution of this work lies in the fact that the subcell force is derived invoking Galilean invariance and thermodynamic consistency. That is, we deduce a general form of the sub-cell force so that a cell entropy inequality is satisfied. The subcell force writes as a pressure contribution plus a tensorial viscous contribution which is proportional to the difference between the nodal velocity and the cell-centered velocity. This cell-centered velocity is a supplementary degree of freedom that is solved by means of a cell-centered approximate Riemann solver. To satisfy the second law of thermodynamics, the local subcell tensor involved in the viscous part of the subcell force must be symmetric positive definite. This subcell tensor is the cornerstone of the scheme. One particular expression of this tensor is given. A high-order extension of this discretization is provided. Numerical tests are presented in order to assess the efficiency of this approach. The results obtained for various representative configurations of one and two-dimensional compressible fluid flows show the robustness and the accuracy of this scheme.

A numerical procedure for the evaluation of equivalent permeability tensor for fractured vuggy porous media is presented. At first we proposed a new conceptual model, i.e., discrete fracture-vug network model, to model the realistic fluid flow in fractured vuggy porous medium on fine scale. This new model consists of three systems: rock matrix system, fractures system, and vugs system. The fractures and vugs are embedded in porous rock, and the isolated vugs could be connected via discrete fracture network. The flow in porous rock and fractures follows Darcy’s law, and the vugs system is free fluid region. Based on two-scale homogenization theory, we obtained an equivalent macroscopic Darcy’s law on coarse scale from fine-scale discrete fracture-vug network model. A finite element numerical formulation for homogenization equations is developed. The method is verified through application to a periodic model problem and then is applied to the calculation of equivalent permeability tensor of porous media with complex fracture-vug networks. The applicability and validity of the method for these more general fractured vuggy systems are assessed through a simple test of the coarse-scale model.

In this paper, a novel approach for quantifying the parametric uncertainty associated with a stochastic problem output is presented. As with Monte-Carlo and stochastic collocation methods, only point-wise evaluations of the stochastic output response surface are required allowing the use of legacy deterministic codes and precluding the need for any dedicated stochastic code to solve the uncertain problem of interest. The new approach differs from these standard methods in that it is based on ideas directly linked to the recently developed compressed sensing theory. The technique allows the retrieval of the modes that contribute most significantly to the approximation of the solution using a minimal amount of information. The generation of this information, via many solver calls, is almost always the bottle-neck of an uncertainty quantification procedure. If the stochastic model output has a reasonably compressible representation in the retained approximation basis, the proposed method makes the best use of the available information and retrieves the dominant modes. Uncertainty quantification of the solution of both a 2-D and 8-D stochastic Shallow Water problem is used to demonstrate the significant performance improvement of the new method, requiring up to several orders of magnitude fewer solver calls than the usual sparse grid-based Polynomial Chaos (Smolyak scheme) to achieve comparable approximation accuracy.

We develop in this paper efficient and robust numerical methods for solving anisotropic Cahn-Hilliard systems. We construct energy stable schemes for the time discretization of the highly nonlinear anisotropic Cahn-Hilliard systems by using a stabilization technique. At each time step, these schemes lead to a sequence of linear coupled elliptic equations with constant coefficients that can be efficiently solved by using a spectral-Galerkin method. We present numerical results that are consistent with earlier work on this topic, and also carry out various simulations, such as the linear bi-Laplacian regularization and the nonlinear Willmore regularization, to demonstrate the efficiency and robustness of the new schemes.

The purpose of this article is to summarize our recent progress in high-order and high accurate CFD methods for flow problems with complex grids as well as to discuss the engineering prospects in using these methods. Despite the rapid development of high-order algorithms in CFD, the applications of high-order and high accurate methods on complex configurations are still limited. One of the main reasons which hinder the widely applications of these methods is the complexity of grids. Many aspects which can be neglected for low-order schemes must be treated carefully for high-order ones when the configurations are complex. In order to implement high-order finite difference schemes on complex multi-block grids, the geometric conservation law and block-interface conditions are discussed. A conservative metric method is applied to calculate the grid derivatives, and a characteristic-based interface condition is employed to fulfil high-order multi-block computing. The fifth-order WCNS-E-5 proposed by Deng is applied to simulate flows with complex grids, including a double-delta wing, a transonic airplane configuration, and a hypersonic X-38 configuration. The results in this paper and the references show pleasant prospects in engineering-oriented applications of high-order schemes.

We present an unconditionally energy stable and uniquely solvable finite difference scheme for the Cahn-Hilliard-Brinkman (CHB) system, which is comprised of a Cahn-Hilliard-type diffusion equation and a generalized Brinkman equation mod-eling fluid flow. The CHB system is a generalization of the Cahn-Hilliard-Stokes model and describes two phase very viscous flows in porous media. The scheme is based on a convex splitting of the discrete CH energy and is semi-implicit. The equations at the implicit time level are nonlinear, but we prove that they represent the gradient of a strictly convex functional and are therefore uniquely solvable, regardless of time step size. Owing to energy stability, we show that the scheme is stable in the time and space discrete and norms. We also present an efficient, practical nonlinear multigrid method . comprised of a standard FAS method for the Cahn-Hilliard part, and a method based on the Vanka smoothing strategy for the Brinkman part . for solving these equations. In particular, we provide evidence that the solver has nearly optimal complexity in typical situations. The solver is applied to simulate spinodal decomposition of a viscous fluid in a porous medium, as well as to the more general problems of buoyancy- and boundary-driven flows.

Vortex and vorticity are two correlated but fundamentally different concepts which have been the central issues in fluid mechanics research. Vorticity has rigorous mathematical definition (curl of velocity), but no clear physical meaning. Vortex has clear physical meaning (rotation) but no rigorous mathematical definition. For a long time, many people treat them as a same thing. However, based on our high-order direct numerical simulation (DNS), we found that first, “vortex” is not “vorticity tube” or “vortex tube” which is widely defined as a bundle of vorticity lines without any vorticity line leak. Actually, vortex is an open area for vorticity line penetration. Second, vortex is not necessarily congregation of vorticity lines, but dispersion in many 3-dimensional cases. Some textbooks say that vortex cannot end inside the flow field but must end on the solid wall (and/or boundaries). Our DNS observation and many other numerical results show almost all vortices are ended inside the flow field. Finally, a more theoretical study shows that neither vortex nor vorticity line can attach to the solid wall and they must be detached from the wall.

A reconstruction-based discontinuous Galerkin method is presented for the solution of the compressible Navier-Stokes equations on arbitrary grids. In this method, an in-cell reconstruction is used to obtain a higher-order polynomial representation of the underlying discontinuous Galerkin polynomial solution and an inter-cell reconstruction is used to obtain a continuous polynomial solution on the union of two neighboring, interface-sharing cells. The in-cell reconstruction is designed to enhance the accuracy of the discontinuous Galerkin method by increasing the order of the underlying polynomial solution. The inter-cell reconstruction is devised to remove an interface discontinuity of the solution and its derivatives and thus to provide a simple, accurate, consistent, and robust approximation to the viscous and heat fluxes in the Navier-Stokes equations. A parallel strategy is also devised for the resulting reconstruction discontinuous Galerkin method, which is based on domain partitioning and Single Program Multiple Data (SPMD) parallel programming model. The RDG method is used to compute a variety of compressible flow problems on arbitrary meshes to demonstrate its accuracy, efficiency, robustness, and versatility. The numerical results demonstrate that this RDG method is third-order accurate at a cost slightly higher than its underlying second-order DG method, at the same time providing a better performance than the third order DG method, in terms of both computing costs and storage requirements.

We study the numerical performance of a continuous data assimilation (downscaling) algorithm, based on ideas from feedback control theory, in the context of the two-dimensional incompressible Navier-Stokes equations. Our model problem is to recover an unknown reference solution, asymptotically in time, by using continuous-in-time coarse-mesh nodal-point observational measurements of the velocity field of this reference solution (subsampling), as might be measured by an array of weather vane anemometers. Our calculations show that the required nodal observation density is remarkably less than what is suggested by the analytical study; and is in fact comparable to the number of numerically determining Fourier modes, which was reported in an earlier computational study by the authors. Thus, this method is computationally efficient and performs far better than the analytical estimates suggest.

We propose a fast local level set method for the inverse problem of gravimetry. The theoretical foundation for our approach is based on the following uniqueness result: if an open set D is star-shaped or x3-convex with respect to its center of gravity, then its exterior potential uniquely determines the open set D. To achieve this purpose constructively, the first challenge is how to parametrize this open set D as its boundary may have a variety of possible shapes. To describe those different shapes we propose to use a level-set function to parametrize the unknown boundary of this open set. The second challenge is how to deal with the issue of partial data as gravimetric measurements are only made on a part of a given reference domain Ω. To overcome this difficulty we propose a linear numerical continuation approach based on the single layer representation to find potentials on the boundary of some artificial domain containing the unknown set D. The third challenge is how to speed up the level set inversion process. Based on some features of the underlying inverse gravimetry problem such as the potential density being constant inside the unknown domain, we propose a novel numerical approach which is able to take advantage of these features so that the computational speed is accelerated by an order of magnitude. We carry out numerical experiments for both two- and three-dimensional cases to demonstrate the effectiveness of the new algorithm.

Cavities and fractures significantly affect the flow paths in carbonate reservoirs and should be accurately accounted for in numerical models. Herein, we consider the problem of computing the effective permeability of rock samples based on high-resolution 3D CT scans containing millions of voxels. We use the Stokes-Brinkman equations in the entire domain, covering regions of free flow governed by the Stokes equations, porous Darcy flow, and transitions between them. The presence of different length scales and large (ten orders of magnitude) contrasts in permeability leads to highly ill-conditioned linear systems of equations, which are difficult to solve. To obtain a problem that is computationally tractable, we first analyze the relative importance of the Stokes and Darcy terms for a set of idealized 2D models. We find that, in terms of effective permeability, the Stokes-Brinkman equations are only applicable for a special parameter set where the effective free-flow permeability is less than four orders of magnitude different from the matrix permeability. All other cases can be accurately modeled with either the Stokes or the Darcy end-member flows, depending on if there do or do not exist percolating free-flow regions. The insights obtained are used to perform a direct computation of the effective permeability of a rock sample model with more than 8 million cells.