A unified gas-kinetic scheme for continuum and rarefied flows

With discretized particle velocity space, a multiscale unified gas-kinetic scheme for entire Knudsen number flows is constructed based on the BGK model. The current scheme couples closely the update of macroscopic conservative variables with the update of microscopic gas distribution function within a time step. In comparison with many existing kinetic schemes for the Boltzmann equation, the current method has no difficulty to get accurate Navier–Stokes (NS) solutions in the continuum flow regime with a time step being much larger than the particle collision time. At the same time, the rarefied flow solution, even in the free molecule limit, can be captured accurately. The unified scheme is an extension of the gas-kinetic BGK-NS scheme from the continuum flow to the rarefied regime with the discretization of particle velocity space. The success of the method is due to the un-splitting treatment of the particle transport and collision in the evaluation of local solution of the gas distribution function. For these methods which use operator splitting technique to solve the transport and collision separately, it is usually required that the time step is less than the particle collision time. This constraint basically makes these methods useless in the continuum flow regime, especially in the high Reynolds number flow simulations. Theoretically, once the physical process of particle transport and collision is modeled statistically by the kinetic Boltzmann equation, the transport and collision become continuous operators in space and time, and their numerical discretization should be done consistently. Due to its multiscale nature of the unified scheme, in the update of macroscopic flow variables, the corresponding heat flux can be modified according to any realistic Prandtl number. Subsequently, this modification effects the equilibrium state in the next time level and the update of microscopic distribution function. Therefore, instead of modifying the collision term of the BGK model, such as ES-BGK and BGK–Shakhov, the unified scheme can achieve the same goal on the numerical level directly. Many numerical tests will be used to validate the unified method.     

A unified incompressible lattice BGK model and its application to three-dimensional lid-driven cavity flow

A unified lattice Bhatnagar－Gross－Krook (ILBGK) model iDdQq for the incompressible Navier－Stokes equation is presented. To test its efficiency, the lid-driven cavity flow in three-dimensional space for Reynolds number Re=3200 and span aspect ratio SAR=1, 2 and 3 is simulated in detail on a 48×48×(48×SAR) uniform lattice using the model. The test results agree well with those in previous experiments and numerical works and show the efficiency and strong numerical stability of the proposed ILBGK model.     

Realizable high-order finite-volume schemes for quadrature-based moment methods

Dilute gas–particle flows can be described by a kinetic equation containing terms for spatial transport, gravity, fluid drag and particle–particle collisions. However, direct numerical solution of kinetic equations is often infeasible because of the large number of independent variables. An alternative is to reformulate the problem in terms of the moments of the velocity distribution. Recently, a quadrature-based moment method was derived for approximating solutions to kinetic equations. The success of the new method is based on a moment-inversion algorithm that is used to calculate non-negative weights and abscissas from the moments. The moment-inversion algorithm does not work if the moments are non-realizable, which might lead to negative weights. It has been recently shown [14] that realizability is guaranteed only with the 1st-order finite-volume scheme that has an inherent problem of excessive numerical diffusion. The use of high-order finite-volume schemes may lead to non-realizable moments. In the present work, realizability of the finite-volume schemes in both space and time is discussed for the 1st time. A generalized idea for developing realizable high-order finite-volume schemes for quadrature-based moment methods is presented. These finite-volume schemes give remarkable improvement in the solutions for a certain class of problems. It is also shown that the standard Runge–Kutta time-integration schemes do not guarantee realizability. However, realizability can be guaranteed if strong stability-preserving (SSP) Runge–Kutta schemes are used. Numerical results are presented on both Cartesian and triangular meshes.     

A Discontinuous Galerkin Method Based on a BGK Scheme for the Navier-Stokes Equations on Arbitrary Grids

A discontinuous Galerkin Method based on a Bhatnagar-Gross-Krook (BGK) formulation is presented for the solution of the compressible Navier-Stokes equations on arbitrary grids. The idea behind this approach is to combine the robustness of the BGK scheme with the accuracy of the DG methods in an effort to develop a more accurate, efficient, and robust method for numerical simulations of viscous flows in a wide range of flow regimes. Unlike the traditional discontinuous Galerkin methods, where a Local Discontinuous Galerkin (LDG) formulation is usually used to discretize the viscous fluxes in the Navier-Stokes equations, this DG method uses a BGK scheme to compute the fluxes which not only couples the convective and dissipative terms together, but also includes both discontinuous and continuous representation in the flux evaluation at a cell interface through a simple hybrid gas distribution function. The developed method is used to compute a variety of viscous flow problems on arbitrary grids. The numerical results obtained by this BGKDG method are extremely promising and encouraging in terms of both accuracy and robustness, indicating its ability and potential to become not just a competitive but simply a superior approach than the current available numerical methods.     

An Equation for the Dissipation Rate Correlation and Its Implications for the Intermittency Exponent μ in Turbulence

We derive an equation satisfied by the dissipation rate correlation function, for the homogeneous, isotropic state of fully-developed turbulence from the the Navier–Stokes equation. In the equal time limit we show that the equation leads directly to two intermittency exponents ₁=2– ₆ and ₂=z₄– ₄, where the 's are exponents of velocity structure functions and z₄ is a dynamical exponent characterizing the fourth order structure function. We discuss the contributions of the pressure terms to the equation and the consequences of hyperscaling.     

Efficient solutions to robust, semi-implicit discretizations of the immersed boundary method

The immersed boundary method is a versatile tool for the investigation of flow-structure interaction. In a large number of applications, the immersed boundaries or structures are very stiff and strong tangential forces on these interfaces induce a well-known, severe time-step restriction for explicit discretizations. This excessive stability constraint can be removed with fully implicit or suitable semi-implicit schemes but at a seemingly prohibitive computational cost. While economical alternatives have been proposed recently for some special cases, there is a practical need for a computationally efficient approach that can be applied more broadly. In this context, we revisit a robust semi-implicit discretization introduced by Peskin in the late 1970s which has received renewed attention recently. This discretization, in which the spreading and interpolation operators are lagged, leads to a linear system of equations for the interface configuration at the future time, when the interfacial force is linear. However, this linear system is large and dense and thus it is challenging to streamline its solution. Moreover, while the same linear system or one of similar structure could potentially be used in Newton-type iterations, nonlinear and highly stiff immersed structures pose additional challenges to iterative methods. In this work, we address these problems and propose cost-effective computational strategies for solving Peskin’s lagged-operators type of discretization. We do this by first constructing a sufficiently accurate approximation to the system’s matrix and we obtain a rigorous estimate for this approximation. This matrix is expeditiously computed by using a combination of pre-calculated values and interpolation. The availability of a matrix allows for more efficient matrix–vector products and facilitates the design of effective iterative schemes. We propose efficient iterative approaches to deal with both linear and nonlinear interfacial forces and simple or complex immersed structures with tethered or untethered points. One of these iterative approaches employs a splitting in which we first solve a linear problem for the interfacial force and then we use a nonlinear iteration to find the interface configuration corresponding to this force. We demonstrate that the proposed approach is several orders of magnitude more efficient than the standard explicit method. In addition to considering the standard elliptical drop test case, we show both the robustness and efficacy of the proposed methodology with a 2D model of a heart valve.     

Scaling Exponents for Active Scalars

We provide bounds for Dirichlet quotients and for generalized structure functions for 3D active scalars and Navier–Stokes equations. These bounds put constraints on the possible extent of anomalous scaling.     

High-order compact exponential finite difference methods for convection–diffusion type problems

A class of high-order compact (HOC) exponential finite difference (FD) methods is proposed for solving one- and two-dimensional steady-state convection–diffusion problems. The newly proposed HOC exponential FD schemes have nonoscillation property and yield high accuracy approximation solution as well as are suitable for convection-dominated problems. The O(h⁴) compact exponential FD schemes developed for the one-dimensional (1D) problems produce diagonally dominant tri-diagonal system of equations which can be solved by applying the tridiagonal Thomas algorithm. For the two-dimensional (2D) problems, O(h⁴ + k⁴) compact exponential FD schemes are formulated on the nine-point 2D stencil and the line iterative approach with alternating direction implicit (ADI) procedure enables us to deal with diagonally dominant tridiagonal matrix equations which can be solved by application of the one-dimensional tridiagonal Thomas algorithm with a considerable saving in computing time. To validate the present HOC exponential FD methods, three linear and nonlinear problems, mostly with boundary or internal layers where sharp gradients may appear due to high Peclet or Reynolds numbers, are numerically solved. Comparisons are made between analytical solutions and numerical results for the currently proposed HOC exponential FD methods and some previously published HOC methods. The present HOC exponential FD methods produce excellent results for all test problems. It is shown that, besides including the excellent performances in computational accuracy, efficiency and stability, the present method has the advantage of better scale resolution. The method developed in this article is easy to implement and has been applied to obtain the numerical solutions of the lid driven cavity flow problem governed by the 2D incompressible Navier–Stokes equations using the stream function-vorticity formulation.     

Statistical Estimates for the Navier–Stokes Equations and the Kraichnan Theory of 2-D Fully Developed Turbulence

A mathematical formulation of the Kraichnan theory for 2-D fully developed turbulence is given in terms of ensemble averages of solutions to the Navier–Stokes equations. A simple condition is given for the enstrophy cascade to hold for wavenumbers just beyond the highest wavenumber of the force up to a fixed fraction of the dissipation wavenumber, up to a logarithmic correction. This is followed by partial rigorous support for Kraichnan's eddy breakup mechanism. A rigorous estimate for the total energy is found to be consistent with Kraichnan's theory. Finally, it is shown that under our conditions for fully developed turbulence the fractal dimension of the attractor obeys a sharper upper bound than in the general case.     

Navier–Stokes Limit for a Thermal Stochastic Lattice Gas

We study a stochastic particle system on the lattice whose particles move freely according to a simple exclusion process and change velocities during collisions preserving energy and momentum. In the hydrodynamic limit, under diffusive space-time scaling, the local velocity field u satisfies the incompressible Navier–Stokes equation, while the temperature field solves the heat equation with drift u. The results are also extended to include a suitably resealed external force.     

A semi-implicit lattice method for simulating flow

We propose a new semi-implicit lattice numerical method for modeling fluid flow that depends only on local primitive variable information (density, pressure, velocity) and not on relaxed upstream distribution function values. This method has the potential for reducing parallel processor communication and permitting larger time steps than the lattice-Boltzmann method. Several benchmark problems are solved to demonstrate the accuracy of the method.     

Simulating 2D Flows with Viscous Vortex Dynamics

Approximate solutions of the two-dimensional Navier–Stokes equation can be constructed as a superposition of viscous Lamb vortices. Requiring minimum deviation from the Navier–Stokes equation, one gets a set of ordinary differential equations for the positions, strength and width of the vortices. We calculate the deviation of the solution from the Navier–Stokes equation in the square norm. The time dependence of this error is determined and discussed.     

Transport Coefficients for Inelastic Maxwell Mixtures

The Boltzmann equation for inelastic Maxwell models (IMM) is used to determine the Navier–Stokes transport coefficients of a granular binary mixture in d-dimensions. The Chapman–Enskog method is applied to solve the Boltzmann equation for states near the (local) homogeneous cooling state. The mass, heat, and momentum fluxes are obtained to first order in the spatial gradients of the hydrodynamic fields, and the corresponding transport coefficients are identified. There are seven relevant transport coefficients: the mutual diffusion, the pressure diffusion, the thermal diffusion, the shear viscosity, the Dufour coefficient, the pressure energy coefficient, and the thermal conductivity. All these coefficients are exactly obtained in terms of the coefficients of restitution and the ratios of mass, concentration, and particle sizes. The results are compared with known transport coefficients of inelastic hard spheres (IHS) obtained analytically in the leading Sonine approximation and by means of Monte Carlo simulations. The comparison shows a reasonably good agreement between both interaction models for not too strong dissipation, especially in the case of the transport coefficients associated with the mass flux     

Navier–Stokes spectral solver in a sphere

The paper presents the first implementation of a primitive variable spectral method for calculating viscous flows inside a sphere. A variational formulation of the Navier–Stokes equations is adopted using a fractional-step time discretization with the classical second-order backward difference scheme combined with explicit extrapolation of the nonlinear term. The resulting scalar and vector elliptic equations are solved by means of the direct spectral solvers developed recently by the authors. The spectral matrices for radial operators are characterized by a minimal sparsity – diagonal stiffness and tridiagonal mass matrix. Closed-form expressions of their nonzero elements are provided here for the first time, showing that the condition number of the relevant matrices grows as the second power of the truncation order. A new spectral elliptic solver for the velocity unknown in spherical coordinates is also described that includes implicitly the Coriolis force in a rotating frame, but requires a minimal coupling between the modal velocity components in the Fourier space. The numerical tests confirm that the proposed method achieves spectral accuracy and ensures infinite differentiability to all orders of the numerical solution, by construction. These results indicate that the new primitive variable spectral solver is an effective alternative to the spectral method recently proposed by Kida and Nakayama, where the velocity field is represented in terms of poloidal and toroidal functions.     

A low-dissipation and time-accurate method for compressible multi-component flow with variable specific heat ratios

A low-dissipation method for calculating multi-component gas dynamics flows with variable specific heat ratio that is capable of accurately simulating flows which contain both high- and low-Mach number features is proposed. The technique combines features from the double-flux multi-component model, nonlinear error-controlled WENO, adaptive TVD slope limiters, rotated Riemann solvers, and adaptive mesh refinement to obtain a method that is both robust and accurate. Success of the technique is demonstrated using an extensive series of numerical experiments including premixed deflagrations, Chapman–Jouget detonations, re-shocked Richtmyer–Meshkov instability, shock-wave and hydrogen gas column interaction, and multi-dimensional detonations. This technique is relatively straight-forward to implement using an existing compressible Navier–Stokes solver based on Godunov’s method.     

Comparison of the generalized Riemann solver and the gas-kinetic scheme for inviscid compressible flow simulations

The generalized Riemann problem (GRP) scheme for the Euler equations and gas-kinetic scheme (GKS) for the Boltzmann equation are two high resolution shock capturing schemes for fluid simulations. The difference is that one is based on the characteristics of the inviscid Euler equations and their wave interactions, and the other is based on the particle transport and collisions. The similarity between them is that both methods can use identical MUSCL-type initial reconstructions around a cell interface, and the spatial slopes on both sides of a cell interface involve in the gas evolution process and the construction of a time-dependent flux function. Although both methods have been applied successfully to the inviscid compressible flow computations, their performances have never been compared. Since both methods use the same initial reconstruction, any difference is solely coming from different underlying mechanism in their flux evaluation. Therefore, such a comparison is important to help us to understand the correspondence between physical modeling and numerical performances. Since GRP is so faithfully solving the inviscid Euler equations, the comparison can be also used to show the validity of solving the Euler equations itself. The numerical comparison shows that the GRP exhibits a slightly better computational efficiency, and has comparable accuracy with GKS for the Euler solutions in 1D case, but the GKS is more robust than GRP. For the 2D high Mach number flow simulations, the GKS is absent from the shock instability and converges to the steady state solutions faster than the GRP. The GRP has carbuncle phenomena, likes a cloud hanging over exact Riemann solvers. The GRP and GKS use different physical processes to describe the flow motion starting from a discontinuity. One is based on the assumption of equilibrium state with infinite number of particle collisions, and the other starts from the non-equilibrium free transport process to evolve into an equilibrium one through particle collisions. The different mechanism in the flux evaluation deviates their numerical performance. Through this study, we may conclude scientifically that it may NOT be valid to use the Euler equations as governing equations to construct numerical fluxes in a discretized space with limited cell resolution. To adapt the Navier–Stokes (NS) equations is NOT valid either because the NS equations describe the flow behavior on the hydrodynamic scale and have no any corresponding physics starting from a discontinuity. This fact alludes to the consistency of the Euler and Navier–Stokes equations with the continuum assumption and the necessity of a direct modeling of the physical process in the discretized space in the construction of numerical scheme when modeling very high Mach number flows. The development of numerical algorithm is similar to the modeling process in deriving the governing equations, but the control volume here cannot be shrunk to zero.     

Modified finite-volume method based on a cell vertex scheme for the solution of radiative transfer problems in complex 3D geometries

A modified finite-volume method based on a cell vertex scheme was applied to solve radiative transfer problems within a participating medium of complex three-dimensional shaped domain. The computational spatial domain of interest was divided into four-node tetrahedron elements with unstructured meshes while the adopted formulation was combined with a closure relation based on an exponential scheme. The studied medium was assumed to be grey, non-scattering and was bounded by black surfaces. Our results were then compared with those found in other articles on the subject. The approach shows a very good level of performance for wall heat transfer evaluation. Accurate results were obtained on coarse computational meshes and solution errors were found to decrease with grid refinement.     

Artificial compressibility, characteristics-based schemes for variable density, incompressible, multi-species flows. Part I. Derivation of different formulations and constant density limit

The paper presents various formulations of characteristics-based schemes in the framework of the artificial-compressibility method for variable-density incompressible flows. In contrast to constant-density incompressible flows, where the characteristics-based variables reconstruction leads to a single formulation, in the case of variable density flows three different schemes can be obtained henceforth labeled as: transport, conservative and hybrid schemes. The conservative scheme results in pseudo-compressibility terms in the (multi-species) density reconstruction. It is shown that in the limit of constant density, the transport scheme becomes the (original) characteristics-based scheme for incompressible flows, but the conservative and hybrid schemes lead to a new characteristics-based variant for constant density flows. The characteristics-based schemes are combined with second and third-order interpolation for increasing the computational accuracy locally at the cell faces of the control volume. Numerical experiments for constant density flows reveal that all the characteristics-based schemes result in the same flow solution, but they exhibit different convergence behavior. The multigrid implementation and numerical studies for variable density flows are presented in Part II of this study.     

Analytical and numerical study of coupled atomistic-continuum methods for fluids

The stability and convergence rate of coupled atomistic-continuum methods are studied analytically and numerically. These methods couple a continuum model with molecular dynamics through the exchange of boundary conditions in the continuum-particle overlapping region. Different coupling schemes, including velocity–velocity, flux–velocity, velocity–flux and flux–flux, are studied. It is found that the velocity–velocity and flux–velocity schemes are stable. The flux–flux scheme is weakly unstable. The stability of the velocity–flux scheme depends on the parameter T_c which is the length of the time interval between successive exchange of boundary conditions. It is stable when T_c is small and unstable when T_c is large. For steady-state problems, the flux–velocity scheme converges faster than the other coupling schemes.     

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.