Search

Inexpensive numerical methods are key to enable simulations of systems of a large number of particles of different shapes in Stokes flow. Several approximate methods have been introduced for this purpose. We study the accuracy of the multiblob method for solving the Stokes mobility problem in free space, where the 3D geometry of a particle surface is discretised with spherical blobs and the pair-wise interaction between blobs is described by the RPY-tensor. The paper aims to investigate and improve on the magnitude of the error in the solution velocities of the Stokes mobility problem using a combination of two different techniques: an optimally chosen grid of blobs and a pair-correction inspired by Stokesian dynamics. Optimisation strategies to determine a grid with a certain number of blobs are presented with the aim of matching the hydrodynamic response of a single accurately described ideal particle, alone in the fluid. Small errors in this self-interaction are essential as they determine the basic error level in a system of well-separated particles. With a good match, reasonable accuracy can be obtained even with coarse blob-resolutions of the particle surfaces. The error in the self-interaction is however sensitive to the exact choice of grid parameters and simply hand-picking a suitable blob geometry can lead to errors several orders of magnitude larger in size. The pair-correction is local and cheap to apply, and reduces on the error for more closely interacting particles. Two different types of geometries are considered: spheres and axisymmetric rods with smooth caps. The error in solutions to mobility problems is quantified for particles of varying inter-particle distances for systems containing a few particles, comparing to an accurate solution based on a second kind BIE-formulation where the quadrature error is controlled by employing quadrature by expansion (QBX)., 49 pages, 37 figures

In this work, we present a positivity-preserving high-order flux reconstruction method for the polyatomic Boltzmann--BGK equation augmented with a discrete velocity model that ensures the scheme is discretely conservative. Through modeling the internal degrees of freedom, the approach is further extended to polyatomic molecules and can encompass arbitrary constitutive laws. The approach is validated on a series of large-scale complex numerical experiments, ranging from shock-dominated flows computed on unstructured grids to direct numerical simulation of three-dimensional compressible turbulent flows, the latter of which is the first instance of such a flow computed by directly solving the Boltzmann equation. The results show the ability of the scheme to directly resolve shock structures without any ad hoc numerical shock capturing method and correctly approximate turbulent flow phenomena in a consistent manner with the hydrodynamic equations., 31 pages, 20 figures

Many multiphase fluid systems, such as those involving immiscible polymers or liquid-liquid systems with surfactants, have shown a breakdown of the no-slip condition at the material interface. This results in systems where the tangential velocity of the inner and outer fluid can differ, with the jump in velocity dependent not only the material properties of the interface but also the stresses applied by the surrounding fluid. In this work a numerical model is presented which is capable of investigating general multiphase fluid systems involving interfacial slip in both two- and three-dimensions. To make the system computationally feasible, a hybrid Navier-Stokes projection method is used, whereby the viscosity, density, and pressure are assumed to be continuous across the interface while the velocity field can experience a jump, which is handled via the Immersed Interface Method. The numerical model is compared to experimental results involving polymer-polymer mixtures and computational results for droplets in extensional flows, showing excellent agreement with both. It is then used to explore the influence of interfacial slip in a number of common multiphase fluid systems, including the shearing of a planar interface, droplet and filament relaxation, and droplets in shear flow, both unbounded and wall-bound.

We present a discrete exterior calculus (DEC) based discretization scheme for incompressible two-phase flows. Our physically-compatible exterior calculus discretization of single phase flow is extended to simulate immiscible two-phase flows with discontinuous changes in fluid properties such as density and viscosity across the interface. The two-phase incompressible Navier-Stokes equations and conservative phase field equation for interface capturing are first transformed into the exterior calculus framework. The discrete counter part of these smooth equations is obtained by substituting with discrete differential forms and discrete exterior operators. We prove the boundedness of the method for the first order Euler forward and predictor-corrector time integration schemes in the DEC framework. With a proper choice of two free parameters, the scheme remains phase field bounded without the requirement of any ad hoc mass redistribution. We verify the scheme against several standard test cases (for interface capturing) comprising not only the flat domains but also the curved domains, leveraging the advantage that DEC operators are independent of the coordinate system. The results show excellent properties of boundedness, mass conservation and convergence. Moreover, we demonstrate the ability of the scheme towards handling large density and viscosity ratios as well as surface tension in the simulation of various two phase flow physical phenomena on flat or curved surfaces., Comment: 27 pages, 25 figures

A simulation framework based on the level-set and the immersed boundary methods (LS-IBM) has been developed for reactive transport problems in porous media involving a moving solid-fluid interface. The interface movement due to surface reactions is tracked by the level-set method, while the immersed boundary method captures the momentum and mass transport at the interface. The proposed method is capable of accurately modeling transport near evolving boundaries in Cartesian grids. The framework formulation guarantees second order accuracy in space. Since the interface velocity is only defined at the moving boundary, an interface velocity propagation method is also proposed. The method can be applied to other moving interface problems of the ``Stefan'' type. Here, we validate the proposed LS-IBM both for flow and transport close to an immersed object with reactive boundaries as well as for crystal growth. The proposed method provides a powerful tool to model more realistic problems involving moving reactive interfaces in complex domains., 61 pages, 30 figures

The dynamics of cavitation bubbles are important in many flows, but their small sizes and high number densities often preclude direct numerical simulation. We present a computational model that averages their effect on the flow over larger spatiotemporal scales. The model is based on solving a generalized population balance equation (PBE) for nonlinear bubble dynamics and explicitly represents the evolving probability density of bubble radii and radial velocities. Conditional quadrature-based moment methods (QBMMs) are adapted to solve this PBE. A one-way-coupled bubble dynamics problem demonstrates the efficacy of different QBMMs for the evolving bubble statistics. Results show that enforcing hyperbolicity during moment inversion (CHyQMOM) provides comparable model-form accuracy to the traditional conditional method of moments and decreases computational costs by about ten times for a broad range of test cases. The CHyQMOM-based computational model is implemented in MFC, an open-source multi-phase and high-order-accurate flow solver. We assess the effect of the model and its parameters on a two-way coupled bubble screen flow problem., Comment: 19 pages, 9 figures, submitted to J. Comp. Phys

We present an interface and geometry preserving (IGP) method for the modeling of fully Eulerian fluid-structure interaction via phase-field formulation. While the hyperbolic tangent interface profile is preserved by the time-dependent mobility model, the proposed method maintains the geometry of the solid-fluid interface by reducing the volume-conserved mean curvature flow. To achieve the reduction in the curvature flow, we construct a gradient-minimizing velocity field (GMV) for the convection of the order parameter. The constructed velocity field enables the preservation of the solid velocity in the solid domain while extending the velocity in the normal direction throughout the diffuse interface region. With this treatment, the GMV reduces the normal velocity difference of the level sets of the order parameter which alleviates the undesired thickening or thinning of the diffuse interface region due to the convection. During this process, the time-dependent mobility coefficient is substantially reduced and there is a lesser curvature flow. The GMV ensures that the diffuse interface region moves with the solid bulk such that the fluid-solid interface conforms to the geometry of the solid. Using the unified momentum equation and the phase-dependent interpolation, we integrate the IGP method into a fully Eulerian variational FSI solver based on the incompressible viscous fluid and the neo-Hookean solid. We first demonstrate the ability of the phase-field-based IGP method for the convection of circular and square interfaces with a prescribed velocity field. The variational FSI framework with the IGP method is then examined for the flow passing a fixed deformable block in a channel domain. Finally, the vibration of a plate attached behind a stationary cylinder subjected to incoming flow is employed to assess the fully Eulerian framework for a large aspect ratio and sharp corners.

In this paper, we introduce a novel way to represent the interface for two-phase flows with phase change. We combine a level-set method with a Cartesian embedded boundary method and take advantage of both. This is part of an effort to obtain a numerical strategy relying on Cartesian grids allowing the simulation of complex boundaries with possible change of topology while retaining a high-order representation of the gradients on the interface and the capability of properly applying boundary conditions on the interface. This leads to a two-fluid conservative second-order numerical method. The ability of the method to correctly solve Stefan problems, onset dendrite growth with and without anisotropy is demonstrated through a variety of test cases. Finally, we take advantage of the two-fluid representation to model a Rayleigh--B\'enard instability with a melting boundary.

A range of optimization cases of two-dimensional Stefan problems, solved using a tracking-type cost-functional, is presented. A level set method is used to capture the interface between the liquid and solid phases and an immersed boundary (cut cell) method coupled with an implicit time-advancement scheme is employed to solve the heat equation. A conservative implicit-explicit scheme is then used for solving the level set transport equation. The resulting numerical framework is validated with respect to existing analytical solutions of the forward Stefan problem. An adjoint-based algorithm is then employed to efficiently compute the gradient used in the optimisation algorithm (L-BFGS). The algorithm follows a continuous adjoint framework, where adjoint equations are formally derived using shape calculus and transport theorems. A wide range of control objectives are presented, and the results show that using parameterised boundary actuation leads to effective control strategies in order to suppress interfacial instabilities or to maintain a desired crystal shape., Comment: 33 pages, 16 figures, preprint submitted to Journal of Computational Physics

The Navier-Stokes-Korteweg (NSK) system is a classical diffuse interface model which is based on van der Waals theory of capillarity. Diffuse interface methods have gained much interest to model two-phase flow in porous media. However, for the numerical solution of the NSK equations two major challenges have to be faced. First, an extended numerical stencil is required due to a third-order term in the linear momentum and the total energy equations. In addition, the dispersive contribution in the linear momentum equations prevents the straightforward use of contact angle boundary conditions. Secondly, any real gas equation of state is based on a non-convex Helmholtz free energy potential which may cause the eigenvalues of the Jacobian of the first-order fluxes to become imaginary numbers inside the spinodal region. In this work, a thermodynamically consistent relaxation model is presented which is used to approximate the NSK equations. The model is complimented by thermodynamically consistent non-equilibrium boundary conditions which take contact angle effects into account. Due to the relaxation approach, the contribution of the Korteweg tensor in the linear momentum and total energy equations can be reduced to second-order terms which enables a straightforward implementation of contact angle boundary conditions in a numerical scheme. Moreover, the definition of a modified pressure function enables to formulate first-order fluxes which remain strictly hyperbolic in the entire spinodal region. The present work is a generalization of a previously presented parabolic relaxation model for the isothermal NSK equations.

Fluid-Structure Interaction (FSI) is a crucial problem in ocean engineering. The smoothed particle hydrodynamics (SPH) method has been employed recently for FSI problems in light of its Lagrangian nature and its advantage in handling multi-physics problems. The efficiency of SPH can be greatly improved with the Adaptive Particle Refinement (APR) method, which refines particles in the regions of interest while deploying coarse particles in the left areas. In this study, the APR method is further improved by developing several new algorithms. Firstly, a new particle refinement strategy with the refinement scale ratio of 4 is employed for multi-level resolutions, and this dramatically decreases the computational costs compared to the standard APR method. Secondly, the regularized transition sub-zone is deployed to render an isotropic particle distribution, which makes the solutions between the refinement zone and the non-refinement zone smoother and consequently results in a more accurate prediction. Thirdly, for complex FSI problems with free surface, a new free-surface detection method based on the Voronoi diagram is proposed, and the performance is validated in comparison to the conventional method. The improved APR method is then applied to a set of challenging FSI cases. Numerical simulations demonstrate that the results from the refinement with scale ratio 4 are consistent with other studies and experimental data, and also agree well with those employing the refinement scale ratio 2. A significant reduction in the computational time is observed for all the considered cases. Overall, the improved APR method with a large refinement scale ratio and the new free-surface detection strategy shows great potential in simulating complex FSI problems efficiently and accurately., Comment: 47 pages, 26 figures, accepted to be published by Journal of Computational Physics

We present a meshless method for magnetohydrodynamics by evolving the vector potential of magnetic fields. A novel scheme and numerical techniques are developed to restrict the divergence of magnetic field, based on the Meshless Finite Mass/Volume with HLLD Riemann solver for conservative flux calculation. We found the magnetic field could be stabilized by a proper smoothing process and so the long-term evolution becomes available. To verify the new scheme, we perform the Brio-Wu shock tube, 2D and 3D Orszag-Tang vortex and Magnetorotational Instability test problems. Our results suggest that our method is robust and has better precision on central offset, amplitude and detailed pattern than an existing meshless code$-$GIZMO., 17 pages, 11 figures, accepted for publication in Journal of Computational Physics

Particles traveling through inertial microfluidic devices migrate to focusing streamlines. We present a numerical method that calculates migration velocities of particles in inertial microfluidic channels of arbitrary cross section by representing particles by singularities. Refinements to asymptotic analysis are given that improve the regularity of the singularity approximation, making finite element approximations of flow and pressure fields more effective. Sample results demonstrate that the method is computationally efficient and able to capture bifurcations in particle focusing positions due to changes in channel shape and Reynolds number.

In this work, we propose a novel phase-field model for the simulation of two-phase flows that is accurate, conservative, bounded, and robust. The proposed model conserves the mass of each of the phases, and results in bounded transport of the volume fraction. We present results from the canonical test cases of a drop advection and a drop in a shear flow, showing significant improvement in the accuracy over the commonly used conservative phase-field method. Moreover, the proposed model imposes a lesser restrictive Courant-Friedrichs-Lewy condition, and hence, is less expensive compared to other conservative phase-field models. We also propose improvements on computation of surface tension forces and show that the proposed improvement significantly reduces the magnitude of spurious velocities at the interface. We also derive a consistent and conservative momentum transport equation for the proposed phase-field model and show that the proposed model when coupled with the consistent momentum transport equation results in discrete conservation of kinetic energy, which is a sufficient condition for the numerical stability of incompressible flows, in the absence of dissipative mechanisms. To illustrate the robustness of the method in simulating high-density ratio turbulent two-phase flows, we present the numerical simulations of high-density ratio droplet-laden isotropic turbulence at finite and infinite Reynolds numbers., Comment: 29 pages, 12 figures

The dynamics of compressible liquid-vapor flow depends sensitively on the microscale behavior at the phase boundary. We consider a sharp-interface approach, and propose a multiscale model to describe liquid-vapor flow accurately, without imposing ad-hoc closure relations on the continuum scale. The multiscale model combines the Euler equations on the continuum scale with molecular-scale particle simulations that govern the interface motion. We rely on an interface-preserving moving mesh finite volume method to discretize the continuum-scale sharp-interface flow in a conservative manner. Computational efficiency, while preserving physical properties, is achieved by a surrogate solver for the interface dynamics based on constraint-aware neural networks. The multiscale model is presented in its general form, and applied to regimes of temperature-dependent liquid-vapor flow which have not been accessible before.

A well-known drawback of the Volume-Of-Fluid (VOF) method is that the breakup of thin liquid films or filaments is mainly caused by numerical aspects rather than by physical ones. The rupture of thin films occurs when their thickness reaches the order of the grid size and by refining the grid the breakup events are delayed. When thin filaments rupture, many droplets are generated due to the mass conserving properties of VOF. Thus, the numerical character of the breakup does not allow obtaining the desired convergence of the droplet size distribution under grid refinement. In this work, we present a novel algorithm to detect and perforate thin structures. First, thin films or ligaments are identified by taking quadratic moments of an indicator obtained from the volume fraction. A multiscale approach allows us to choose the critical film thickness independently of the mesh resolution. Then, the breakup is induced by making holes in the films before their thickness reaches the grid size. We show that the method improves the convergence upon grid refinement of the droplets size distribution and of enstrophy.

A noteworthy aspect in blood flow modeling is the definition of the mechanical interaction between the fluid flow and the biological structure that contains it, namely the vessel wall. It has been demonstrated that the addition of a viscous contribution to the mechanical characterization of vessels brings positive results when compared to in-vivo measurements. In this context, the numerical implementation of boundary conditions able to keep memory of the viscoelastic contribution of vessel walls assumes an important role, especially when dealing with large circulatory systems. In this work, viscoelasticity is taken into account in entire networks via the Standard Linear Solid Model. The implementation of the viscoelastic contribution at boundaries (inlet, outlet and junctions), is carried out considering the hyperbolic nature of the mathematical model. A non-linear system is established based on the definition of the Riemann Problem at junctions, characterized by rarefaction waves separated by contact discontinuities, among which the mass and the total energy are conserved. Basic junction tests are analyzed, such as a trivial 2-vessels junction, for both a generic artery and vein, and a simple 3-vessels junction, considering an aortic bifurcation scenario. The chosen asymptotic-preserving IMEX Runge-Kutta Finite Volume scheme is demonstrated to be second-order accurate in the whole domain and well-balanced, even when including junctions. Two different benchmark models of the arterial network are implemented, differing in number of vessels and in viscoelastic parameters. Comparison of the results obtained in the two networks underlines the high sensitivity of the model to the chosen viscoelastic parameters. The conservation of the contribution provided by the viscoelastic characterization of vessel walls is assessed in the whole network, including junctions and boundary conditions., 42 pages, 19 figures

International audience; We present the construction of an original stochastic model for the instantaneous turbulent kinetic energy at a given point of a flow, and we validate estimator methods on this model with observational data examples. Motivated by the need for wind energy industry of acquiring relevant statistical information of air motion at a local place, we adopt the Lagrangian description of fluid flows to derive, from the 3D+time equations of the physics, a 0D+time-stochastic model for the time series of the instantaneous turbulent kinetic energy at a given position. Specifically, based on the Lagrangian stochastic description of a generic fluid-particles, we derive a family of mean-field dynamics featuring the square norm of the turbulent velocity. By approximating at equilibrium the characteristic nonlinear terms of the dynamics, we recover the so called Cox-Ingersoll-Ross process, which was previously suggested in the literature for modelling wind speed. We then propose a calibration procedure for the parameters employing both direct methods and Bayesian inference. In particular, we show the consistency of the estimators and validate the model through the quantification of uncertainty, with respect to the range of values given in the literature for some physical constants of turbulence modelling.

Permeability estimation of porous media from directly solving the Navier–Stokes equations has a wide spectrum of applications in petroleum industry. In this paper, we utilize a pressure-correction projection algorithm in conjunction with the interior penalty discontinuous Galerkin scheme for space discretization to build an incompressible Navier–Stokes simulator and to use this simulator to calculate permeability of real rock samples. The proposed method is accurate, numerically robust, and exhibits the potential for tackling realistic problems.

This paper presents a robust, adaptive numerical scheme for simulating high density ratio and high shear multiphase flows on locally refined Cartesian grids that adapt to the evolving interfaces and track regions of high vorticity. The algorithm combines the interface capturing level set method with a variable-coefficient incompressible Navier-Stokes solver that is demonstrated to stably resolve material contrast ratios of up to six orders of magnitude. The discretization approach ensures second-order pointwise accuracy for both velocity and pressure with several physical boundary treatments, including velocity and traction boundary conditions. The paper includes several test cases that demonstrate the order of accuracy and algorithmic scalability of the flow solver. To ensure the stability of the numerical scheme in the presence of high density and viscosity ratios, we employ a consistent treatment of mass and momentum transport in the conservative form of discrete equations. This consistency is achieved by solving an additional mass balance equation, which we approximate via a strong stability preserving Runga-Kutta time integrator and by employing the same mass flux (obtained from the mass equation) in the discrete momentum equation. The scheme uses higher-order total variation diminishing (TVD) and convection-boundedness criterion (CBC) satisfying limiter to avoid numerical fluctuations in the transported density field. The high-order bounded convective transport is done on a dimension-by-dimension basis, which makes the scheme simple to implement. We also demonstrate through several test cases that the lack of consistent mass and momentum transport in non-conservative formulations, which are commonly used in practice, or the use of non-CBC satisfying limiters can yield very large numerical error and very poor accuracy for convection-dominant high density ratio flows., Comment: Figures are compressed to comply with arXiv size requirements

We present a highly accurate numerical method based on a boundary integral formulation and the leaky dielectric model to study the dynamics of surfactant-covered drops in the presence of an applied electric field. The method can simulate interacting 3D drops (no axisymmetric simplification) in close proximity, can consider different viscosities, is adaptive in time and able to handle substantial drop deformation. For each drop global representations of the variables based on spherical harmonics expansions are used and the spectral accuracy is achieved by designing specific numerical tools: a specialized quadrature method for the singular and nearly singular integrals that appear in the formulation, a general preconditioner for the implicit treatment of the surfactant diffusion and a reparametrization procedure able to ensure a high-quality representation of the drops also under deformation. Our numerical method is validated against theoretical, numerical and experimental results available in the literature, as well as a new second-order theory developed for a surfactant-laden drop placed in a quadrupole electric field.

We present a numerical method for computing the single layer (Stokeslet) and double layer (stresslet) integrals in Stokes flow. The method applies to smooth, closed surfaces in three dimensions, and achieves high accuracy both on and near the surface. The singular Stokeslet and stresslet kernels are regularized and, for the nearly singular case, corrections are added to reduce the regularization error. These corrections are derived analytically for both the Stokeslet and the stresslet using local asymptotic analysis. For the case of evaluating the integrals on the surface, as needed when solving integral equations, we design high order regularizations for both kernels that do not require corrections. This approach is direct in that it does not require grid refinement or special quadrature near the singularity, and therefore does not increase the computational complexity of the overall algorithm. Numerical tests demonstrate the uniform convergence rates for several surfaces in both the singular and near singular cases, as well as the importance of corrections when two surfaces are close to each other.

Simulations involving free surfaces and fluid interfaces are important in many areas of engineering. There is, however, still a need for improved simulation methods. Recently, a new efficient geometric VOF method called isoAdvector for general polyhedral meshes was published. We investigate the interface reconstruction step of isoAdvector, and demonstrate that especially for unstructured meshes the applied isosurface based approach can lead to noisy interface orientations. We then introduce a novel computational interface reconstruction scheme based on calculation of a reconstructed distance function (RDF). By iterating over the RDF calculation and interface reconstruction, we obtain second order convergence of both the interface normal and position within cells even with a strict L ∞ error norm. In 2D this is verified with reconstruction of a circle on Cartesian meshes and on unstructured triangular and polygonal prism meshes. In 3D the second order convergence is verified with reconstruction of a sphere on Cartesian meshes and on unstructured tetrahedral and polyhedral meshes. The new scheme is combined with the interface advection step of the isoAdvector algorithm. Significantly reduced absolute advection errors are obtained, and for CFL number 0.2 and below we demonstrate second order convergence on all the mentioned mesh types in 2D and 3D. The implementation of the proposed interface reconstruction schemes is straightforward and the computational cost is significantly reduced compared to contemporary methods. The schemes are implemented as an extension to the Computational Fluid Dynamics (CFD) Open Source software package, OpenFOAM®. The extension module and all test cases presented in this paper are released as open source.

This paper designs and analyzes positivity-preserving well-balanced (WB) central discontinuous Galerkin (CDG) schemes for the Euler equations with gravity. A distinctive feature of these schemes is that they not only are WB for a general known stationary hydrostatic solution, but also can preserve the positivity of the fluid density and pressure. The standard CDG method does not possess this feature, while directly applying some existing WB techniques to the CDG framework may not accommodate the positivity and keep other important properties at the same time. In order to obtain the WB and positivity-preserving properties simultaneously while also maintaining the conservativeness and stability of the schemes, a novel spatial discretization is devised in the CDG framework based on suitable modifications to the numerical dissipation term and the source term approximation. The modifications are based on a crucial projection operator for the stationary hydrostatic solution, which is proposed for the first time in this work. This novel projection has the same order of accuracy as the standard $L^2$-projection, can be explicitly calculated, and is easy to implement without solving any optimization problems. More importantly, it ensures that the projected stationary solution has the same cell averages on both the primal and dual meshes, which is a key to achieve the desired properties of our schemes. Based on some convex decomposition techniques, rigorous positivity-preserving analyses for the resulting WB CDG schemes are carried out. Several one- and two-dimensional numerical examples are performed to illustrate the desired properties of these schemes, including the high-order accuracy, the WB property, the robustness for simulations involving the low pressure or density, high resolution for the discontinuous solutions and the small perturbations around the equilibrium state., Comment: 57 pages

Computational fluctuating hydrodynamics aims at understanding the impact of thermal fluctuations on fluid motions at small scales through numerical exploration. These fluctuations are modeled as stochastic flux terms and incorporated into the classical Navier-Stokes equations, which need to be solved numerically. In this paper, we present a novel projection-based method for solving the incompressible fluctuating hydrodynamics (FHD) equations. By analyzing the equilibrium structure factor spectrum of the velocity field for the linearized FHD equations, we investigate how the inherent splitting errors affect the numerical solution of the stochastic partial differential equations in the presence of non-periodic boundary conditions, and how iterative corrections can reduce these errors. Our computational examples demonstrate both the capability of our approach to reproduce correctly stochastic properties of fluids at small scales as well as its potential use in the simulations of multi-physics problems., Comment: To appear in J. Comput. Phys

Overset methods are commonly employed to enable the effective simulation of problems involving complex geometries and moving objects such as rotorcraft. This paper presents a novel overset domain connectivity algorithm based upon the direct cut approach suitable for use with GPU-accelerated solvers on high-order curved grids. In contrast to previous methods it is capable of exploiting the highly data-parallel nature of modern accelerators. Further, the approach is also substantially more efficient at handling the curved grids which arise within the context of high-order methods. An implementation of this new algorithm is presented and combined with a high-order fluid dynamics code. The algorithm is validated against several benchmark problems, including flow over a spinning golf ball at a Reynolds number of 150,000., Comment: Preprint accepted for publication in the Journal of Computational Physics. From previous version: Fixed typos Fixed incorrect plot of spinning golf ball forces; Reorganized discussion of code performance for clarity; Added additional code performance profiling results

It is well-known that by placing judiciously chosen image point forces and doublets to the Stokeslet above a flat wall, the no-slip boundary condition can be conveniently imposed on the wall [Blake, J. R. Math. Proc. Camb. Philos. Soc. 70(2), 1971: 303.]. However, to further impose periodic boundary conditions on directions parallel to the wall usually involves tedious derivations because single or double periodicity in Stokes flow may require the periodic unit to have no net force, which is not satisfied by the well-known image system. In this work we present a force-neutral image system. This neutrality allows us to represent the Stokes image system in a universal formulation for non-periodic, singly periodic and doubly periodic geometries. This formulation enables the black-box style usage of fast kernel summation methods. We demonstrate the efficiency and accuracy of this new image method with the periodic kernel independent fast multipole method in both non-periodic and doubly periodic geometries. We then extend this new image system to other widely used Stokes fundamental solutions, including the Laplacian of the Stokeslet and the Rotne-Prager-Yamakawa tensor., Comment: To appear in Journal of Computational Physics

We present a three-dimensional (3D) common-refinement method for non-matching meshes between discrete non-overlapping subdomains of incompressible fluid and nonlinear hyperelastic structure. To begin, we first investigate the accuracy of common-refinement method (CRM) to satisfy traction equilibrium condition along the fluid-elastic interface with non-matching meshes. We systematically assess the accuracy of CRM against the matching grid solution by varying grid mismatch between the fluid and solid meshes over a cylindrical tubular elastic body. We demonstrate second-order accuracy of CRM through uniform refinements of fluid and solid meshes along the interface. We then extend the error analysis to transient data transfer across non-matching meshes between fluid and solid solvers. We show that the common-refinement discretization across non-matching fluid-structure grids yields accurate transfer of the physical quantities across the fluid-solid interface. We next solve a 3D benchmark problem of a cantilevered hyperelastic plate behind a circular bluff body and verify the accuracy of coupled solutions with respect to the available solution in the literature. By varying the solid interface resolution, we generate various non-matching grid ratios and quantify the accuracy of CRM for the nonlinear structure interacting with a laminar flow. We illustrate that the CRM with the partitioned NIFC treatment is stable for low solid-to-fluid density ratio and non-matching meshes. Finally, we demonstrate the 3D parallel implementation of common-refinement with NIFC scheme for a realistic engineering problem of drilling riser undergoing complex vortex-induced vibration with strong added mass effects., 38 pages, 16 figures

We extend the recently introduced divergence-conforming immersed boundary (DCIB) method [1] to fluid-structure interaction (FSI) problems involving closed co-dimension one solids. We focus on capsules and vesicles, whose discretization is particularly challenging due to the higher-order derivatives that appear in their formulations. In two-dimensional settings, we employ cubic B-splines with periodic knot vectors to obtain discretizations of closed curves with C^2 inter-element continuity. In three-dimensional settings, we use analysis-suitable bi-cubic T-splines to obtain discretizations of closed surfaces with at least C^1 inter-element continuity. Large spurious changes of the fluid volume inside closed co-dimension one solids is a well-known issue for IB methods. The DCIB method results in volume changes orders of magnitude lower than conventional IB methods. This is a byproduct of discretizing the velocity-pressure pair with divergence-conforming B-splines, which lead to negligible incompressibility errors at the Eulerian level. The higher inter-element continuity of divergence-conforming B-splines is also crucial to avoid the quadrature/interpolation errors of IB methods becoming the dominant discretization error. Benchmark and application problems of vesicle and capsule dynamics are solved, including mesh-independence studies and comparisons with other numerical methods., For supplementary movies go to https://www.andrew.cmu.edu/user/hugocp/research.html

In this paper a fully Eulerian solver for the study of multiphase flows for simulating the propagation of surface gravity waves over submerged bodies is presented. We solve the incompressible Navier–Stokes equations coupled with the volume of fluid technique for the modeling of the liquid phases with the interface, an immersed body method for the solid objects and an iterative strong–coupling procedure for the fluid–structure interaction. The flow incompressibility is enforced via the solution of a Poisson equation which, owing to the density jump across the interfaces of the liquid phases, has to resort to the splitting procedure of Dodd and Ferrante [1] , [1] . The solver is validated through comparisons against classical test cases for fluid–structure interaction like migration of particles in a pressure–driven channel, multiphase flows, ‘water exit’ of a cylinder and a good agreement is found for all tests. Furthermore, we show the application of the solver to the case of a surface gravity wave propagating over a submerged reversed pendulum and verify that the solver can reproduce the energy exchange between the wave and the pendulum. Finally the three–dimensional spilling breaking of a wave induced by a submerged sphere is considered.

Volume conservation and shape preservation are two well-known issues related to the advection and remeshing in front tracking. To address these issues, this paper proposes a divergence-preserving velocity interpolation method and a parabolic fit vertex positioning method for remeshing operations for three-dimensional front tracking. Errors in preserving the divergence of the velocity field when interpolating the velocity from the fluid mesh to the vertices of the triangles of the front are a primary reason for volume conservation errors when advecting the front. The proposed interpolation method preserves the discrete divergence of the fluid velocity by construction and is compared in this work with other known interpolation methods in divergence-free and non-divergence-free test cases, with respect to volume conservation and shape preservation of the front. The presented interpolation method conserves the volume and shape up to an order of magnitude better than the conventionally used interpolation methods and is within the range of higher order interpolation methods at lower computational cost. Additionally, the parabolic fit vertex positioning method for remeshing operations locally approximates the front with a smooth polynomial surface, improving volume conservation and shape preservation by an order of magnitude compared to conventional remeshing algorithms.

High-order methods and hybrid turbulence models have independently shown promise as means of decreasing the computational cost of scale-resolving simulations. The objective of this work is to develop the combination of these methods and analyze the effects of high-order discretizations on hybrid turbulence models, particularly with respect the optimal model parameters and the relative accuracy benefits compared to approaches such as under-resolved direct numerical simulation (URDNS). We employ the Partially-Averaged Navier-Stokes (PANS) approach using the flux reconstruction scheme on the flow around a periodic hill and the wake flow of a circular cylinder at a Reynolds number of 3900, the latter of which we provide direct numerical simulation results and novel statistical analysis. By increasing the order of the discretization while fixing the total degrees of freedom, it was observed that larger improvements in the prediction of the statistics and flow physics were generally seen with PANS than URDNS. Furthermore, less sensitivity to the resolution-control parameter was observed with a high-order discretization, indicating that high-order discretizations may be an effective approach for increasing the accuracy and reliability of hybrid turbulence models for scale-resolving simulations without a significant increase in computational effort., Comment: 20 pages, 13 figures

Entropy-Stable (ES) schemes, specifically those built from [Tadmor \textit{Math. Comput.} 49 (1987) 91], have been gaining interest over the past decade, especially in the context of under-resolved simulations of compressible turbulent flows using high-order methods. These schemes are attractive because they can provide stability in a global and nonlinear sense (consistency with thermodynamics). However, fully realizing the potential of ES schemes requires a better grasp of their local behavior. Entropy-stability itself does not imply good local behavior [Gouasmi \textit{et al.} \textit{J. Sci. Comp.} 78 (2019) 971, Gouasmi \textit{et al.} \textit{Comput. Methd. Appl. M.} 363 (2020) 112912]. In this spirit, we studied ES schemes in problems where \textit{global stability is not the core issue}. In the present work, we consider the accuracy degradation issues typically encountered by upwind-type schemes in the low-Mach-number regime [Turkel \textit{Annu. Rev. Fluid Mech.} 31 (1999) 285] and their treatment using \textit{Flux-Preconditioning} [Turkel \textit{J. Comput. Phys.} 72 (1987) 277, Miczek \textit{et al.} \textit{A \& A} 576 (2015) A50]. ES schemes suffer from the same issues and Flux-Preconditioning can improve their behavior without interfering with entropy-stability. This is first demonstrated analytically: using similarity and congruence transforms we were able to establish conditions for a preconditioned flux to be ES, and introduce the ES variants of the Miczek's and Turkel's preconditioned fluxes. This is then demonstrated numerically through first-order simulations of two simple test problems representative of the incompressible and acoustic limits, the Gresho Vortex and a right-moving acoustic wave. The results are overall consistent with previous studies [...]

Lagrangian data assimilation of complex nonlinear turbulent flows is an important but computationally challenging topic. In this article, an efficient data-driven statistically accurate reduced-order modeling algorithm is developed that significantly accelerates the computational efficiency of Lagrangian data assimilation. The algorithm starts with a Fourier transform of the high-dimensional flow field, which is followed by an effective model reduction that retains only a small subset of the Fourier coefficients corresponding to the energetic modes. Then a linear stochastic model is developed to approximate the nonlinear dynamics of each Fourier coefficient. Effective additive and multiplicative noise processes are incorporated to characterize the modes that exhibit Gaussian and non-Gaussian statistics, respectively. All the parameters in the reduced order system, including the multiplicative noise coefficients, are determined systematically via closed analytic formulae. These linear stochastic models succeed in forecasting the uncertainty and facilitate an extremely rapid data assimilation scheme. The new Lagrangian data assimilation is then applied to observations of sea ice floe trajectories that are driven by atmospheric winds and turbulent ocean currents. It is shown that observing only about $30$ non-interacting floes in a $200$km$\times200$km domain is sufficient to recover the key multi-scale features of the ocean currents. The additional observations of the floe angular displacements are found to be suitable supplements to the center-of-mass positions for improving the data assimilation skill. In addition, the observed large and small floes are more useful in recovering the large- and small-scale features of the ocean, respectively. The Fourier domain data assimilation also succeeds in recovering the ocean features in the areas where cloud cover obscures the observations.

We develop a rapid and accurate contour method for the solution of time-fractional PDEs. The method inverts the Laplace transform via an optimised stable quadrature rule, suitable for infinite-dimensional operators, whose error decreases like $\exp(-cN/\log(N))$ for $N$ quadrature points. The method is parallisable, avoids having to resolve singularities of the solution as $t\downarrow 0$, and avoids the large memory consumption that can be a challenge for time-stepping methods applied to time-fractional PDEs. The ODEs resulting from quadrature are solved using adaptive sparse spectral methods that converge exponentially with optimal linear complexity. These solutions of ODEs are reused for different times. We provide a complete analysis of our approach for fractional beam equations used to model small-amplitude vibration of viscoelastic materials with a fractional Kelvin-Voigt stress-strain relationship. We calculate the system's energy evolution over time and the surface deformation in cases of both constant and non-constant viscoelastic parameters. An infinite-dimensional ``solve-then-discretise'' approach considerably simplifies the analysis, which studies the generalisation of the numerical range of a quasi-linearisation of a suitable operator pencil. This allows us to build an efficient algorithm with explicit error control. The approach can be readily adapted to other time-fractional PDEs and is not constrained to fractional parameters in the range $0

One of the central problems in the study of rarefied gas dynamics is to find the steady-state solution of the Boltzmann equation quickly. When the Knudsen number is large, i.e. the system is highly rarefied, the conventional iterative scheme can lead to convergence within a few iterations. However, when the Knudsen number is small, i.e. the flow falls in the near-continuum regime, hundreds of thousands iterations are needed, and yet the “converged” solutions are prone to be contaminated by accumulated error and large numerical dissipation. Recently, based on the gas kinetic models, the implicit unified gas kinetic scheme (UGKS) and its variants have significantly reduced the number of iterations in the near-continuum flow regime, but still much higher than that of the highly rarefied gas flows. In this paper, we put forward a general synthetic iterative scheme (GSIS) to find the steady-state solutions of rarefied gas flows within dozens of iterations at any Knudsen number. The key ingredient of our scheme is that the macroscopic equations, which are solved together with the Boltzmann equation and help to adjust the velocity distribution function, not only asymptotically preserve the Navier-Stokes limit in the framework of Chapman-Enskog expansion, but also contain the Newton's law for stress and the Fourier's law for heat conduction explicitly. For this reason, like the implicit UGKS, the constraint that the spatial cell size should be smaller than the mean free path of gas molecules is removed, but we do not need the complex evaluation of numerical flux at cell interfaces. What's more, as the GSIS does not rely on the specific collision operator, it can be naturally extended to quickly find converged solutions for mixture flows and even flows involving chemical reactions. These two superior advantages are expected to accelerate the slow convergence in the simulation of near-continuum flows via the direct simulation Monte Carlo method and its low-variance version.

With the aim of better understanding the numerical properties of the lattice Boltzmann method (LBM), a general methodology is proposed to derive its hydrodynamic limits in the discrete setting. It relies on a Taylor expansion in the limit of low Knudsen numbers. With a single asymptotic analysis, two kinds of deviations with the Navier-Stokes (NS) equations are explicitly evidenced: consistency errors, inherited from the kinetic description of the LBM, and numerical errors attributed to its space and time discretization. The methodology is applied to the Bhatnagar-Gross-Krook (BGK), the regularized and the multiple relaxation time (MRT) collision models in the isothermal framework. Deviation terms are systematically confronted to linear analyses in order to validate their expressions, interpret them and provide explanations for their numerical properties. The low dissipation of the BGK model is then related to a particular pattern of its error terms in the Taylor expansion. Similarly, dissipation properties of the regularized and MRT models are explained by a phenomenon referred to as hyperviscous degeneracy. The latter consists in an unexpected resurgence of high-order Knudsen effects induced by a large numerical pre-factor. It is at the origin of over-dissipation and severe instabilities in the low-viscosity regime., Comment: 59 pages, 16 figures

This article introduces a representation of dynamic meshes, adapted to some numerical simulations that require controlling the volume of objects with free boundaries, such as incompressible fluid simulation, some astrophysical simulations at cosmological scale, and shape/topology optimization. The algorithm decomposes the simulated object into a set of convex cells called a Laguerre diagram, parameterized by the position of $N$ points in 3D and $N$ additional parameters that control the volumes of the cells. These parameters are found as the (unique) solution of a convex optimization problem -- semi-discrete Monge-Amp\`ere equation -- stemming from optimal transport theory. In this article, this setting is extended to objects with free boundaries and arbitrary topology, evolving in a domain of arbitrary shape, by solving a partial optimal transport problem. The resulting Lagrangian scheme makes it possible to accurately control the volume of the object, while precisely tracking interfaces, interactions, collisions, and topology changes., Comment: 37 pages, 17 figures

In the present study, a consistent and conservative Phase-Field model is developed to study thermo-gas-liquid-solid flows with liquid-solid phase change. The proposed model is derived with the help of the consistency conditions and exactly reduces to the consistent and conservative Phase-Field method for incompressible two-phase flows, the fictitious domain Brinkman penalization (FD/BP) method for fluid-structure interactions, and the Phase-Field model of solidification of pure material. It honors the mass conservation, defines the volume fractions of individual phases unambiguously, and therefore captures the volume change due to phase change. The momentum is conserved when the solid phase is absent, but it changes when the solid phase appears due to the no-slip condition at the solid boundary. The proposed model also conserves the energy, preserves the temperature equilibrium, and is Galilean invariant. A novel continuous surface tension force to confine its contribution at the gas-liquid interface and a drag force modified from the Carman-Kozeny equation to reduce solid velocity to zero are proposed. The issue of initiating phase change in the original Phase-Field model of solidification is addressed by physically modifying the interpolation function. The corresponding consistent scheme is developed to solve the model, and the numerical results agree well with the analytical solutions and the existing experimental and numerical data. Two challenging problems having a wide range of material properties and complex dynamics are conducted to demonstrate the capability of the proposed model., This is an accepted manuscript

A new sub-grid-scale model is developed for studying influences of the Hall term on macroscopic aspects of magnetohydrodynamic turbulence. Although the Hall term makes numerical simulations extremely expensive by exciting high-wave-number coefficients and makes magnetohydrodynamic equations stiff, studying macroscopic aspects of magnetohydrodynamic turbulence together with the Hall term is meaningful since this term often influences not only sub-ion-scales but also macroscopic scales. Our new sub-ion-scale sub-grid-scale model for large eddy simulations of Hall magnetohydrodynamic turbulence is developed in order to overcome the difficulties. Large eddy simulations by the use of the new model successfully reproduce statistical signatures such as the energies and probability density functions of the vorticity and current density, keeping some signatures intrinsic to Hall magnetohydrodynamic turbulence. Our new sub-grid-scale model enables numerical simulations of homogeneous and isotropic Hall magnetohydrodynamic turbulence with a small computational cost, improving the effective resolution of an LES from that carried out with earlier models, and retaining the ion-electron separation effects by the Hall term in the grid scales.

Large eddy simulation (LES) has been increasingly used to tackle vortex-dominated turbulent flows. In LES, the quality of the simulation results hinges upon the quality of the numerical discretizations in both space and time. It is in this context we perform a Fourier analysis of several popular methods in LES including the discontinuous Galerkin (DG), finite difference (FD), and compact difference (CD) methods. We begin by reviewing the semi-discrete versions of all methods under consideration, followed by a fully-discrete analysis with explicit Runge–Kutta (RK) time integration schemes. In this regard, we are able to unravel the true dispersion/dissipation behavior of DG and Runge–Kutta DG (RKDG) schemes for the entire wavenumber range using a combined-mode analysis. In this approach, we take into account all eigenmodes in DG and RKDG schemes. The physical-mode is verified to be a good approximation for the asymptotic behavior of these DG schemes in the low wavenumber range. After that, we proceed to compare the DG, FD, and CD methods in dispersion and dissipation properties. Numerical tests are conducted using the linear advection equation to verify the analysis. In comparing different methods, it is found that the overall numerical dissipation strongly depends on the time step. Compact difference (CD) and central FD schemes, in some particular settings, can have more numerical dissipation than the DG scheme with an upwind flux. This claim is then verified through a numerical test using the Burgers' equation.

In this paper, a fourth-order compact gas-kinetic scheme (GKS) is developed for the compressible Euler and Navier–Stokes equations under the framework of two-stage fourth-order temporal discretization and Hermite WENO (HWENO) reconstruction. Due to the high-order gas evolution model, the GKS provides a time dependent gas distribution function at a cell interface. This time evolution solution can be used not only for the flux evaluation across a cell interface and its time derivative, but also time accurate flow variables at a cell interface. As a result, besides updating the conservative flow variables inside each control volume, the cell averaged slopes inside each control volume through the differences of flow variables at the cell interfaces can be updated as well in GKS. So, with the updated flow variables and their slopes inside each cell, the HWENO techniques can be naturally implemented for the compact high-order reconstruction at the beginning of each time step. Therefore, a compact higher-order GKS, such as the two-stage fourth-order compact scheme can be constructed. This fourth-order compact GKS has the same stencil and as robust as the second-order scheme. In comparison with the fourth-order DG method, they have the same stencil. In order to get a fourth-order temporal accuracy, the GKS uses two stages and DG needs four stages in the time stepping methods. The CFL number used in GKS is on the order of 0.5 instead of 0.11 in the DG. This research concludes that beyond the first-order Riemann solver the use of high-order time evolution model at a cell interface is extremely helpful in the design of robust, accurate, and efficient higher-order compact schemes for the compressible flow simulations.

We present a coupled Eulerian-Lagrangian method to simulate cloud cavitation in a compressible liquid. The method is designed to capture the strong, volumetric oscillations of each bubble and the bubble-scattered acoustics. The dynamics of the bubbly mixture is formulated using volume-averaged equations of motion. The continuous phase is discretized on an Eulerian grid and integrated using a high-order, finite-volume weighted essentially non-oscillatory (WENO) scheme, while the gas phase is modeled as spherical, Lagrangian point-bubbles at the sub-grid scale, each of whose radial evolution is tracked by solving the Keller-Miksis equation. The volume of bubbles is mapped onto the Eulerian grid as the void fraction by using a regularization (smearing) kernel. In the most general case, where the bubble distribution is arbitrary, three-dimensional Cartesian grids are used for spatial discretization. In order to reduce the computational cost for problems possessing translational or rotational homogeneities, we spatially average the governing equations along the direction of symmetry and discretize the continuous phase on two-dimensional or axi-symmetric grids, respectively. We specify a regularization kernel that maps the three-dimensional distribution of bubbles onto the field of an averaged two-dimensional or axi-symmetric void fraction. A closure is developed to model the pressure fluctuations at the sub-grid scale as synthetic noise. For the examples considered here, modeling the sub-grid pressure fluctuations as white noise agrees a priori with computed distributions from three-dimensional simulations, and suffices, a posteriori, to accurately reproduce the statistics of the bubble dynamics. The numerical method and its verification are described by considering test cases of the dynamics of a single bubble and cloud cavitaiton induced by ultrasound fields., Comment: 28 pages, 16 figures

A computational framework based on nonlinear direct-adjoint looping is presented for optimizing mixing strategies for binary fluid systems. The governing equations are the nonlinear Navier-Stokes equations, augmented by an evolution equation for a passive scalar, which are solved by a spectral Fourier-based method. The stirrers are embedded in the computational domain by a Brinkman-penalization technique, and shape and path gradients for the stirrers are computed from the adjoint solution. Four cases of increasing complexity are considered, which demonstrate the efficiency and effectiveness of the computational approach and algorithm. Significant improvements in mixing efficiency, within the externally imposed bounds, are achieved in all cases., 23 pages, 55 figures

Fluid-structure systems occur in a range of scientific and engineering applications. The immersed boundary(IB) method is a widely recognized and effective modeling paradigm for simulating fluid-structure interaction(FSI) in such systems, but a difficulty of the IB formulation is that the pressure and viscous stress are generally discontinuous at the interface. The conventional IB method regularizes these discontinuities, which typically yields low-order accuracy at these interfaces. The immersed interface method(IIM) is an IB-like approach to FSI that sharply imposes stress jump conditions, enabling higher-order accuracy, but prior applications of the IIM have been largely restricted to methods that rely on smooth representations of the interface geometry. This paper introduces an IIM that uses only a C0 representation of the interface,such as those provided by standard nodal Lagrangian FE methods. Verification examples for models with prescribed motion demonstrate that the method sharply resolves stress discontinuities along the IB while avoiding the need for analytic information of the interface geometry. We demonstrate that only the lowest-order jump conditions for the pressure and velocity gradient are required to realize global 2nd-order accuracy. Specifically,we show 2nd-order global convergence rate along with nearly 2nd-order local convergence in the Eulerian velocity, and between 1st-and 2nd-order global convergence rates along with 1st-order local convergence for the Eulerian pressure. We also show 2nd-order local convergence in the interfacial displacement and velocity along with 1st-order local convergence in the fluid traction. As a demonstration of the method's ability to tackle complex geometries,this approach is also used to simulate flow in an anatomical model of the inferior vena cava., Added a non-axisymmetric example (flow within eccentric rotating cylinder in Sec. 4.3) - Added a more in-depth analysis and comparison with a body-fitted approach for the application in Sec. 4.7

An exponential time-integrator scheme of second-order accuracy based on the predictor-corrector methodology, denoted PCEXP, is developed to solve multi-dimensional nonlinear partial differential equations pertaining to fluid dynamics. The effective and efficient implementation of PCEXP is realized by means of the Krylov method. The linear stability and truncation error are analyzed through a one-dimensional model equation. The proposed PCEXP scheme is applied to the Euler equations discretized with a discontinuous Galerkin method in both two and three dimensions. The effectiveness and efficiency of the PCEXP scheme are demonstrated for both steady and unsteady inviscid flows. The accuracy and efficiency of the PCEXP scheme are verified and validated through comparisons with the explicit third-order total variation diminishing Runge-Kutta scheme (TVDRK3), the implicit backward Euler (BE) and the implicit second-order backward difference formula (BDF2). For unsteady flows, the PCEXP scheme generates a temporal error much smaller than the BDF2 scheme does, while maintaining the expected acceleration at the same time. Moreover, the PCEXP scheme is also shown to achieve the computational efficiency comparable to the implicit schemes for steady flows., Comment: Exponential scheme for Computational Fluid Dynamics (CFD). 37 pages, 11 figures

We simulate complex fluids by means of an on-the-fly coupling of the bulk rheology to the underlying microstructure dynamics. In particular, a macroscopic continuum model of polymeric fluids is constructed without a pre-specified constitutive relation, but instead it is actively learned from mesoscopic simulations where the dynamics of polymer chains is explicitly computed. To couple the macroscopic rheology of polymeric fluids and the microscale dynamics of polymer chains, the continuum approach (based on the finite volume method) provides the transient flow field as inputs for the (mesoscopic) dissipative particle dynamics (DPD), and in turn DPD returns an effective constitutive relation to close the continuum equations. In this multiscale modeling procedure, we employ an active learning strategy based on Gaussian process regression (GPR) to minimize the number of expensive DPD simulations, where adaptively selected DPD simulations are performed only as necessary. Numerical experiments are carried out for flow past a circular cylinder of a non-Newtonian fluid, modeled at the mesoscopic level by bead-spring chains. The results show that only five DPD simulations are required to achieve an effective closure of the continuum equations at Reynolds number Re=10. Furthermore, when Re is increased to 100, only one additional DPD simulation is required for constructing an extended GPR-informed model closure. Compared to traditional message-passing multiscale approaches, applying an active learning scheme to multiscale modeling of non-Newtonian fluids can significantly increase the computational efficiency. Although the method demonstrated here obtains only a local viscosity from the mesoscopic model, it can be extended to other multiscale models of complex fluids whose macro-rheology is unknown., Comment: 18 pages, 9 figures

We present a reduction-consistent and thermodynamically consistent formulation and an associated numerical algorithm for simulating the dynamics of an isothermal mixture consisting of $N$ ($N\geqslant 2$) immiscible incompressible fluids with different physical properties (densities, viscosities, and pair-wise surface tensions). By reduction consistency we refer to the property that if only a set of $M$ ($1\leqslant M\leqslant N-1$) fluids are present in the system then the N-phase governing equations and boundary conditions will exactly reduce to those for the corresponding $M$-phase system. By theromdynamic consistency we refer to the property that the formulation honors the thermodynamic principles. Our N-phase formulation is developed based on a more general method that allows for the systematic construction of reduction-consistent formulations, and the method suggests the existence of many possible forms of reduction-consistent and thermodynamically consistent N-phase formulations. Extensive numerical experiments have been presented for flow problems involving multiple fluid components and large density ratios and large viscosity ratios, and the simulation results are compared with the physical theories or the available physical solutions. The comparisons demonstrate that our method produces physically accurate results for this class of problems., 51 pages, 15 figures, 8 tables

This paper introduces a new approach to Large-Eddy Simulation (LES) where subgrid-scale (SGS) dissipation is applied proportionally to the degree of local spectral broadening, hence mitigated or deactivated in regions dominated by large-scale and/or laminar vortical motion. The proposed coherent-vorticity preserving (CvP) LES methodology is based on the evaluation of the ratio of the test-filtered to resolved (or grid-filtered) enstrophy, σ. Values of σ close to 1 indicate low sub-test-filter turbulent activity, justifying local deactivation of the SGS dissipation. The intensity of the SGS dissipation is progressively increased for σ 1 which corresponds to a small-scale spectral broadening. The SGS dissipation is then fully activated in developed turbulence characterized by σ ≤ σ e q , where the value σ e q is derived assuming a Kolmogorov spectrum. The proposed approach can be applied to any eddy-viscosity model, is algorithmically simple and computationally inexpensive. LES of Taylor–Green vortex breakdown demonstrates that the CvP methodology improves the performance of traditional, non-dynamic dissipative SGS models, capturing the peak of total turbulent kinetic energy dissipation during transition. Similar accuracy is obtained by adopting Germano's dynamic procedure albeit at more than twice the computational overhead. A CvP-LES of a pair of unstable periodic helical vortices is shown to predict accurately the experimentally observed growth rate using coarse resolutions. The ability of the CvP methodology to dynamically sort the coherent, large-scale motion from the smaller, broadband scales during transition is demonstrated via flow visualizations. LES of compressible channel are carried out and show a good match with a reference DNS.

In this paper we consider discontinuous Galerkin (DG) methods for the incompressible Navier-Stokes equations in the framework of projection methods. In particular we employ symmetric interior penalty DG methods within the second-order rotational incremental pressure correction scheme. The major focus of the paper is threefold: i) We propose a modified upwind scheme based on the Vijayasundaram numerical flux that has favourable properties in the context of DG. ii) We present a novel postprocessing technique in the Helmholtz projection step based on $H(\text{div})$ reconstruction of the pressure correction that is computed locally, is a projection in the discrete setting and ensures that the projected velocity satisfies the discrete continuity equation exactly. As a consequence it also provides local mass conservation of the projected velocity. iii) Numerical results demonstrate the properties of the scheme for different polynomial degrees applied to two-dimensional problems with known solution as well as large-scale three-dimensional problems. In particular we address second-order convergence in time of the splitting scheme as well as its long-time stability., Comment: 24 pages, revised submission to Journal of Computational Physics