The interactions between incompressible fluid flows and immersed structures are nonlinear multi-physics phenomena that have applications to a wide range of scientific and engineering disciplines. In this article, we review representative numerical methods based on conforming and non-conforming meshes that are currently available for computing fluid-structure interaction problems, with an emphasis on some of the recent developments in the field. A goal is to categorize the selected methods and assess their accuracy and efficiency. We discuss challenges faced by researchers in this field, and we emphasize the importance of interdisciplinary effort for advancing the study in fluid-structure interactions.

In this paper, we review the recent development of phase-field models and their numerical methods for multi-component fluid flows with interfacial phenomena. The models consist of a Navier-Stokes system coupled with a multi-component Cahn-Hilliard system through a phase-field dependent surface tension force, variable density and viscosity, and the advection term. The classical infinitely thin boundary of separation between two immiscible fluids is replaced by a transition region of a small but finite width, across which the composition of the mixture changes continuously. A constant level set of the phase-field is used to capture the interface between two immiscible fluids. Phase-field methods are capable of computing topological changes such as splitting and merging, and thus have been applied successfully to multi-component fluid flows involving large interface deformations. Practical applications are provided to illustrate the usefulness of using a phase-field method. Computational results of various experiments show the accuracy and effectiveness of phase-field models.

The computational work and storage of numerically solving the time fractional PDEs are generally huge for the traditional direct methods since they require total memory and work, where NT and NS represent the total number of time steps and grid points in space, respectively. To overcome this difficulty, we present an efficient algorithm for the evaluation of the Caputo fractional derivative of order α∈(0,1). The algorithm is based on an efficient sum-of-exponentials (SOE) approximation for the kernel t–1–α on the interval [Δt, T] with a uniform absolute error ε. We give the theoretical analysis to show that the number of exponentials Nexp needed is of order for T≫1 or for TH1 for fixed accuracy ε. The resulting algorithm requires only storage and work when numerically solving the time fractional PDEs. Furthermore, we also give the stability and error analysis of the new scheme, and present several numerical examples to demonstrate the performance of our scheme.

Centered numerical fluxes can be constructed for compressible Euler equations which preserve kinetic energy in the semi-discrete finite volume scheme. The essential feature is that the momentum flux should be of the form where and are any consistent approximations to the pressure and the mass flux. This scheme thus leaves most terms in the numerical flux unspecified and various authors have used simple averaging. Here we enforce approximate or exact entropy consistency which leads to a unique choice of all the terms in the numerical fluxes. As a consequence novel entropy conservative flux that also preserves kinetic energy for the semi-discrete finite volume scheme has been proposed. These fluxes are centered and some dissipation has to be added if shocks are present or if the mesh is coarse. We construct scalar artificial dissipation terms which are kinetic energy stable and satisfy approximate/exact entropy condition. Secondly, we use entropy-variable based matrix dissipation flux which leads to kinetic energy and entropy stable schemes. These schemes are shown to be free of entropy violating solutions unlike the original Roe scheme. For hypersonic flows a blended scheme is proposed which gives carbuncle free solutions for blunt body flows. Numerical results for Euler and Navier-Stokes equations are presented to demonstrate the performance of the different schemes.

With discretized particle velocity space, a multi-scale unified gas-kinetic scheme for entire Knudsen number flows has been constructed based on the kinetic model in one-dimensional case [J. Comput. Phys., vol. 229 (2010), pp. 7747-7764]. For the kinetic equation, to extend a one-dimensional scheme to multidimensional flow is not so straightforward. The major factor is that addition of one dimension in physical space causes the distribution function to become two-dimensional, rather than axially symmetric, in velocity space. In this paper, a unified gas-kinetic scheme based on the Shakhov model in two-dimensional space will be presented. Instead of particle-based modeling for the rarefied flow, such as the direct simulation Monte Carlo (DSMC) method, the philosophical principal underlying the current study is a partial-differential-equation (PDE)-based modeling. Since the valid scale of the kinetic equation and the scale of mesh size and time step may be significantly different, the gas evolution in a discretized space is modeled with the help of kinetic equation, instead of directly solving the partial differential equation. Due to the use of both hydrodynamic and kinetic scales flow physics in a gas evolution model at the cell interface, the unified scheme can basically present accurate solution in all flow regimes from the free molecule to the Navier-Stokes solutions. In comparison with the DSMC and Navier-Stokes flow solvers, the current method is much more efficient than DSMC in low speed transition and continuum flow regimes, and it has better capability than NS solver in capturing of non-equilibrium flow physics in the transition and rarefied flow regimes. As a result, the current method can be useful in the flow simulation where both continuum and rarefied flow physics needs to be resolved in a single computation. This paper will extensively evaluate the performance of the unified scheme from free molecule to continuum NS solutions, and from low speed micro-flow to high speed non-equilibrium aerodynamics. The test cases clearly demonstrate that the unified scheme is a reliable method for the rarefied flow computations, and the scheme provides an important tool in the study of non-equilibrium flow.

We provide an overview of current techniques and typical applications of numerical bifurcation analysis in fluid dynamical problems. Many of these problems are characterized by high-dimensional dynamical systems which undergo transitions as parameters are changed. The computation of the critical conditions associated with these transitions, popularly referred to as ‘tipping points’, is important for understanding the transition mechanisms. We describe the two basic classes of methods of numerical bifurcation analysis, which differ in the explicit or implicit use of the Jacobian matrix of the dynamical system. The numerical challenges involved in both methods arementioned and possible solutions to current bottlenecks are given. To demonstrate that numerical bifurcation techniques are not restricted to relatively low-dimensional dynamical systems, we provide several examples of the application of the modern techniques to a diverse set of fluid mechanical problems.

Parallel computing has become an important subject in the field of computer science and has proven to be critical when researching high performance solutions. The evolution of computer architectures (multi-core and many-core) towards a higher number of cores can only confirm that parallelism is the method of choice for speeding up an algorithm. In the last decade, the graphics processing unit, or GPU, has gained an important place in the field of high performance computing (HPC) because of its low cost and massive parallel processing power. Super-computing has become, for the first time, available to anyone at the price of a desktop computer. In this paper, we survey the concept of parallel computing and especially GPU computing. Achieving efficient parallel algorithms for the GPU is not a trivial task, there are several technical restrictions that must be satisfied in order to achieve the expected performance. Some of these limitations are consequences of the underlying architecture of the GPU and the theoretical models behind it. Our goal is to present a set of theoretical and technical concepts that are often required to understand the GPU and its massive parallelism model. In particular, we show how this new technology can help the field of computational physics, especially when the problem is data-parallel. We present four examples of computational physics problems; n-body, collision detection, Potts model and cellular automata simulations. These examples well represent the kind of problems that are suitable for GPU computing. By understanding the GPU architecture and its massive parallelism programming model, one can overcome many of the technical limitations found along the way, design better GPU-based algorithms for computational physics problems and achieve speedups that can reach up to two orders of magnitude when compared to sequential implementations.

This paper is concerned with a new version of the Osher-Solomon Riemann solver and is based on a numerical integration of the path-dependent dissipation matrix. The resulting scheme is much simpler than the original one and is applicable to general hyperbolic conservation laws, while retaining the attractive features of the original solver: the method is entropy-satisfying, differentiable and complete in the sense that it attributes a different numerical viscosity to each characteristic field, in particular to the intermediate ones, since the full eigenstructure of the underlying hyperbolic system is used. To illustrate the potential of the proposed scheme we show applications to the following hyperbolic conservation laws: Euler equations of compressible gas-dynamics with ideal gas and real gas equation of state, classical and relativistic MHD equations as well as the equations of nonlinear elasticity. To the knowledge of the authors, apart from the Euler equations with ideal gas, an Osher-type scheme has never been devised before for any of these complicated PDE systems. Since our new general Riemann solver can be directly used as a building block of high order finite volume and discontinuous Galerkin schemes we also show the extension to higher order of accuracy and multiple space dimensions in the new framework of PNPM schemes on unstructured meshes recently proposed in [9].

We propose a positivity preserving entropy decreasing finite volume scheme for nonlinear nonlocal equations with a gradient flow structure. These properties allow for accurate computations of stationary states and long-time asymptotics demonstrated by suitably chosen test cases in which these features of the scheme are essential. The proposed scheme is able to cope with non-smooth stationary states, different time scales including metastability, as well as concentrations and self-similar behavior induced by singular nonlocal kernels. We use the scheme to explore properties of these equations beyond their present theoretical knowledge.

The fractional derivatives include nonlocal information and thus their calculation requires huge storage and computational cost for long time simulations. We present an efficient and high-order accurate numerical formula to speed up the evaluation of the Caputo fractional derivative based on the L2-1σ formula proposed in [A. Alikhanov, J. Comput. Phys., 280 (2015), pp. 424-438], and employing the sum-of-exponentials approximation to the kernel function appeared in the Caputo fractional derivative. Both theoretically and numerically, we prove that while applied to solving time fractional diffusion equations, our scheme not only has unconditional stability and high accuracy but also reduces the storage and computational cost.

A code that models signals produced by charge-exchange reactions between fast ions and injected neutral beams in tokamak plasmas is described. With the fast-ion distribution function as input, the code predicts the efflux to a neutral particle analyzer (NPA) diagnostic and the photon radiance of Balmer-alpha light to a fast-ion Dα (FIDA) diagnostic. Reactions with both the primary injected neutrals and with the cloud of secondary “halo” neutrals that surround the beam are treated. Accurate calculation of the fraction of neutrals that occupy excited atomic states (the collisional-radiative transition equations) is an important element of the code. Comparison with TRANSP output and other tests verify the solutions. Judicious selection of grid size and other parameters facilitate efficient solutions. The output of the code has been validated by FIDA measurements on DIII-D but further tests are warranted.

An all speed scheme for the Isentropic Euler equations is presented in this paper. When the Mach number tends to zero, the compressible Euler equations converge to their incompressible counterpart, in which the density becomes a constant. Increasing approximation errors and severe stability constraints are the main difficulty in the low Mach regime. The key idea of our all speed scheme is the special semi-implicit time discretization, in which the low Mach number stiff term is divided into two parts, one being treated explicitly and the other one implicitly. Moreover, the flux of the density equation is also treated implicitly and an elliptic type equation is derived to obtain the density. In this way, the correct limit can be captured without requesting the mesh size and time step to be smaller than the Mach number. Compared with previous semi-implicit methods, firstly, nonphysical oscillations can be suppressed by choosing proper parameter, besides, only a linear elliptic equation needs to be solved implicitly which reduces much computational cost. We develop this semi-implicit time discretization in the framework of a first order Local Lax-Friedrichs (or Rusanov) scheme and numerical tests are displayed to demonstrate its performances.

In this paper we consider two commonly used classes of finite volume weighted essentially non-oscillatory (WENO) schemes in two dimensional Cartesian meshes. We compare them in terms of accuracy, performance for smooth and shocked solutions, and efficiency in CPU timing. For linear systems both schemes are high order accurate, however for nonlinear systems, analysis and numerical simulation results verify that one of them (Class A) is only second order accurate, while the other (Class B) is high order accurate. The WENO scheme in Class A is easier to implement and costs less than that in Class B. Numerical experiments indicate that the resolution for shocked problems is often comparable for schemes in both classes for the same building blocks and meshes, despite of the difference in their formal order of accuracy. The results in this paper may give some guidance in the application of high order finite volume schemes for simulating shocked flows.

The computation of compressible flows becomes more challenging when the Mach number has different orders of magnitude. When the Mach number is of order one, modern shock capturing methods are able to capture shocks and other complex structures with high numerical resolutions. However, if the Mach number is small, the acoustic waves lead to stiffness in time and excessively large numerical viscosity, thus demanding much smaller time step and mesh size than normally needed for incompressible flow simulation. In this paper, we develop an all-speed asymptotic preserving (AP) numerical scheme for the compressible isentropic Euler and Navier-Stokes equations that is uniformly stable and accurate for all Mach numbers. Our idea is to split the system into two parts: one involves a slow, nonlinear and conservative hyperbolic system adequate for the use of modern shock capturing methods and the other a linear hyperbolic system which contains the stiff acoustic dynamics, to be solved implicitly. This implicit part is reformulated into a standard pressure Poisson projection system and thus possesses sufficient structure for efficient fast Fourier transform solution techniques. In the zero Mach number limit, the scheme automatically becomes a projection method-like incompressible solver. We present numerical results in one and two dimensions in both compressible and incompressible regimes.

This paper studies the numerical simulations for the Cahn-Hilliard equation which describes a phase separation phenomenon. The numerical simulation of the Cahn-Hilliard model needs very long time to reach the steady state, and therefore large time-stepping methods become useful. The main objective of this work is to construct the unconditionally energy stable finite difference scheme so that the large time steps can be used in the numerical simulations. The equation is discretized by the central difference scheme in space and fully implicit second-order scheme in time. The proposed scheme is proved to be unconditionally energy stable and mass-conservative. An error estimate for the numerical solution is also obtained with second order in both space and time. By using this energy stable scheme, an adaptive time-stepping strategy is proposed, which selects time steps adaptively based on the variation of the free energy against time. The numerical experiments are presented to demonstrate the effectiveness of the adaptive time-stepping approach.

In this paper, we investigate the coupling of the Multi-dimensional Optimal Order Detection (MOOD) method and the Arbitrary high order DERivatives (ADER) approach in order to design a new high order accurate, robust and computationally efficient Finite Volume (FV) scheme dedicated to solve nonlinear systems of hyperbolic conservation laws on unstructured triangular and tetrahedral meshes in two and three space dimensions, respectively. The Multi-dimensional Optimal Order Detection (MOOD) method for 2D and 3D geometries has been introduced in a recent series of papers for mixed unstructured meshes. It is an arbitrary high-order accurate Finite Volume scheme in space, using polynomial reconstructions with a posteriori detection and polynomial degree decrementing processes to deal with shock waves and other discontinuities. In the following work, the time discretization is performed with an elegant and efficient one-step ADER procedure. Doing so, we retain the good properties of the MOOD scheme, that is to say the optimal high-order of accuracy is reached on smooth solutions, while spurious oscillations near singularities are prevented. The ADER technique permits not only to reduce the cost of the overall scheme as shown on a set of numerical tests in 2D and 3D, but it also increases the stability of the overall scheme. A systematic comparison between classical unstructured ADER-WENO schemes and the new ADER-MOOD approach has been carried out for high-order schemes in space and time in terms of cost, robustness, accuracy and efficiency. The main finding of this paper is that the combination of ADER with MOOD generally outperforms the one of ADER and WENO either because at given accuracy MOOD is less expensive (memory and/or CPU time), or because it is more accurate for a given grid resolution. A large suite of classical numerical test problems has been solved on unstructured meshes for three challenging multi-dimensional systems of conservation laws: the Euler equations of compressible gas dynamics, the classical equations of ideal magneto-Hydrodynamics (MHD) and finally the relativistic MHD equations (RMHD), which constitutes a particularly challenging nonlinear system of hyperbolic partial differential equation. All tests are run on genuinely unstructured grids composed of simplex elements.

This paper establishes relations between the stability and the high-order truncated corrections for modeling of the mass conservation equation with the two-relaxation-times (TRT) collision operator. First we propose a simple method to derive the truncation errors from the exact, central-difference type, recurrence equations of the TRT scheme. They also supply its equivalent three-time-level discretization form. Two different relationships of the two relaxation rates nullify the third (advection) and fourth (pure diffusion) truncation errors, for any linear equilibrium and any velocity set. However, the two relaxation times alone cannot remove the leading-order advection-d if fusion error, because of the intrinsic fourth-order numerical diffusion. The truncation analysis is carefully verified for the evolution of concentration waves with the anisotropic diffusion tensors. The anisotropic equilibrium functions are presented in a simple but general form, suitable for the minimal velocity sets and the d2Q9, d3Q13, d3Q15 and d3Q19 velocity sets. All anisotropic schemes are complemented by their exact necessary von Neumann stability conditions and equivalent finite-difference stencils. The sufficient stability conditions are proposed for the most stable (OTRT) family, which enables modeling at any Peclet numbers with the same velocity amplitude. The heuristic stability analysis of the fourth-order truncated corrections extends the optimal stability to larger relationships of the two relaxation rates, in agreement with the exact (one-dimensional) and numerical (multi-dimensional) stability analysis. A special attention is put on the choice of the equilibrium weights. By combining accuracy and stability predictions, several strategies for selecting the relaxation and free-tunable equilibrium parameters are suggested and applied to the evolution of the Gaussian hill.

A robust and efficient solver coupling computational fluid dynamics (CFD) with discrete element method (DEM) is developed to simulate particle-laden flows in various physical settings. An interpolation algorithm suitable for unstructured meshes is proposed to translate between mesh-based Eulerian fields and particle-based La-grangian quantities. The interpolation scheme reduces the mesh-dependence of the averaging and interpolation procedures. In addition, the fluid-particle interaction terms are treated semi-implicitly in this algorithm to improve stability and to maintain accuracy. Finally, it is demonstrated that sub-stepping is desirable for fluid-particle systems with small Stokes numbers. A momentum-conserving sub-stepping technique is introduced into the fluid-particle coupling procedure, so that problems with a wide range of time scales can be solved without resorting to excessively small time steps in the CFD solver. Several numerical examples are presented to demonstrate the capabilities of the solver and the merits of the algorithm.

We present a new Finite Volume Evolution Galerkin (FVEG) scheme for the solution of the shallow water equations (SWE) with the bottom topography as a source term. Our new scheme will be based on the FVEG methods presented in (Noelle and Kraft, J. Comp. Phys., 221 (2007)), but adds the possibility to handle dry boundaries. The most important aspect is to preserve the positivity of the water height. We present a general approach to ensure this for arbitrary finite volume schemes. The main idea is to limit the outgoing fluxes of a cell whenever they would create negative water height. Physically, this corresponds to the absence of fluxes in the presence of vacuum. Wellbalancing is then re-established by splitting gravitational and gravity driven parts of the flux. Moreover, a new entropy fix is introduced that improves the reproduction of sonic rarefaction waves.

Due to the rapid advances in micro-electro-mechanical systems (MEMS), the study of microflows becomes increasingly important. Currently, the molecular-based simulation techniques are the most reliable methods for rarefied flow computation, even though these methods face statistical scattering problem in the low speed limit. With discretized particle velocity space, a unified gas-kinetic scheme (UGKS) for entire Knudsen number flow has been constructed recently for flow computation. Contrary to the particle-based direct simulation Monte Carlo (DSMC) method, the unified scheme is a partial differential equation-based modeling method, where the statistical noise is totally removed. But, the common point between the DSMC and UGKS is that both methods are constructed through direct modeling in the discretized space. Due to the multiscale modeling in the unified method, i.e., the update of both macroscopic flow variables and microscopic gas distribution function, the conventional constraint of time step being less than the particle collision time in many direct Boltzmann solvers is released here. The numerical tests show that the unified scheme is more efficient than the particle-based methods in the low speed rarefied flow computation. The main purpose of the current study is to validate the accuracy of the unified scheme in the capturing of non-equilibrium flow phenomena. In the continuum and free molecular limits, the gas distribution function used in the unified scheme for the flux evaluation at a cell interface goes to the corresponding Navier-Stokes and free molecular solutions. In the transition regime, the DSMC solution will be used for the validation of UGKS results. This study shows that the unified scheme is indeed a reliable and accurate flow solver for low speed non-equilibrium flows. It not only recovers the DSMC results whenever available, but also provides high resolution results in cases where the DSMC can hardly afford the computational cost. In thermal creep flow simulation, surprising solution, such as the gas flowing from hot to cold regions along the wall surface, is observed for the first time by the unified scheme, which is confirmed later through intensive DSMC computation.