Simple modification techniques are proposed for making numerical fluxes amenable to unrealizable states (e.g., negative density) without degrading the design order of accuracy, so that a finite-volume solver never fails with unrealizable states arising in the solution reconstruction step and continues to run. The main idea is to evaluate quantities not affecting the order of accuracy but important for stabilization, e.g., a dissipation matrix, with low-order unreconstructed solutions. For the viscous flux, the viscosity is linearly extrapolated instead of being evaluated with linearly reconstructed temperatures to avoid a failure with a negative temperature. These ideas are quite general and may be applied to a wide range of numerical fluxes. In this paper, we illustrate them with the Roe flux and the alpha-damping viscous flux and demonstrate their effectiveness for cases, where a conventional technique encounters difficulties.

Conventional explicit electromagnetic particle-in-cell (PIC) algorithms do not conserve discrete energy exactly. Time-centered fully implicit PIC algorithms can conserve discrete energy exactly, but may introduce large dispersion errors in the light-wave modes. This can lead to intolerable simulation errors where accurate light propagation is needed (e.g. in laser-plasma interactions). In this study, we selectively combine the leap-frog and Crank-Nicolson methods to produce an exactly energy- and charge-conserving relativistic electromagnetic PIC algorithm. Specifically, we employ the leap-frog method for Maxwell's equations, and the Crank-Nicolson method for the particle equations. The semi-implicit formulation still features a timestep CFL, but facilitates exact global energy conservation, exact local charge conservation, and preserves the dispersion properties of the leap-frog method for the light wave. The algorithm employs a new particle pusher designed to maximize efficiency and minimize wall-clock-time impact vs. the explicit alternative. It has been implemented in a code named iVPIC, based on the Los Alamos National Laboratory VPIC code (https://github.com/losalamos/vpic). We present numerical results that demonstrate the properties of the scheme with sample test problems: relativistic two-stream instability, Weibel instability, and laser-plasma instabilities.

We present a simple numerical algorithm for solving elliptic equations where the diffusion coefficient, the source term, the solution and its flux are discontinuous across an irregular interface. The algorithm produces second-order accurate solutions and first-order accurate gradients in the L∞-norm on Cartesian grids. The condition number is bounded, regardless of the ratio of the diffusion constant and scales like that of the standard 5-point stencil approximation on a rectangular grid with no interface. Numerical examples are given in two and three spatial dimensions.

This work presents a numerical algorithm for using radial basis function-generated finite differences (RBF-FD) to solve partial differential equations (PDEs) on S2 using polyharmonic splines with added polynomials defined in a 2D plane (PHS+Poly). We introduce a novel method for calculating RBF-FD PHS+Poly differentiation weights on S2 using first a Householder reflection and then a projection onto the tangent plane. The new PHS+Poly RBF-FD method is implemented on two standard test cases: 1) solid body rotation on S2 and 2) 3D tracer transport within Earth's atmosphere. Compared to existing methods (including those in the RBF literature) at similar resolutions, this approach requires fewer degrees of freedom and is algorithmically much simpler. A MATLAB code to implement the method is included in the Appendix.

In this study, convection-pressure split Euler flux functions which contain weakly hyperbolic convective subsystems are analyzed. A system of first-order partial differential equations (PDEs) is said to be weakly hyperbolic if the corresponding flux Jacobian does not contain a complete set of linearly independent (LI) eigenvectors. Thus, the application of existing flux difference splitting (FDS) based schemes, which depend heavily on both eigenvalues and eigenvectors, are non-trivial to such systems. In the case of weakly hyperbolic systems, a required set of LI eigenvectors can be constructed through the addition of generalized eigenvectors by utilizing the theory of Jordan canonical forms. Once this is achieved for a weakly hyperbolic convective subsystem, an upwind solver can be constructed in the splitting framework. In the present work, the above approach is used for developing two new schemes. The first scheme is based on the Zha–Bilgen type splitting while the second is based on the Toro–Vázquez splitting. Both the schemes are tested on various benchmark problems in one-dimension (1-D) and two-dimensions (2-D). The concept of generalized eigenvectors based on Jordan forms is found to be useful in dealing with the weakly hyperbolic parts of the considered Euler systems.

We consider the initial boundary value problem for the time-fractional diffusion equation with a homogeneous Dirichlet boundary condition and an inhomogeneous initial data a(x)∈L2(D) in a bounded domain D⊂Rd with a sufficiently smooth boundary. We analyze the homogenized solution under the assumption that the diffusion coefficient κϵ(x) is smooth and periodic with the period ϵ>0 being sufficiently small. We derive that its first order approximation measured by both pointwise-in-time in L2(D) and Lp((θ,T);H1(D)) for p∈[1,∞) and θ∈(0,T) has a convergence rate of O(ϵ1/2) when the dimension d≤2 and O(ϵ1/6) when d=3. Several numerical tests are presented to demonstrate the performance of the first order approximation.

In this paper, the band structures of nanoscale phononic crystals based on the nonlocal elasticity theory are calculated by using a meshfree local radial basis function collocation method (LRBFCM). The direct method is applied to enhance the stability of the derivative calculations in the LRBFCM. A simple summation in the LRBFCM is proposed to deal with the integration related to the nonlocal stresses or tractions. The LRBFCM for the band structure calculations is validated by the results obtained with the first principle and the transfer matrix (TM) method for one-dimensional (1D) phononic crystals, as well as the comparison of the frequency responses of the two-dimensional (2D) periodical structures.

We propose a mesh refinement technique for solving elliptic difference equations on unbounded domains based on the fast lattice Green's function (FLGF) method. The FLGF method exploits the regularity of the Cartesian mesh and uses the fast multipole method in conjunction with fast Fourier transforms to yield linear complexity and decrease time-to-solution. We extend this method to a multi-resolution scheme and allow for locally refined Cartesian blocks embedded in the computational domain. Appropriately chosen interpolation and regularization operators retain consistency between the discrete Laplace operator and its inverse on the unbounded domain. Second-order accuracy and linear complexity are maintained, while significantly reducing the number of degrees of freedom and hence the computational cost.

Discontinuous Galerkin (DG) methods have become mainstays in the accurate solution of hyperbolic systems, which suggests that they should also be important for computational electrodynamics (CED). Typically DG schemes are coupled with Runge-Kutta timestepping, resulting in RKDG schemes, which are also sometimes called DGTD schemes in the CED community. However, Maxwell's equations, which are solved in CED codes, have global mimetic constraints. In Balsara and Käppeli [von Neumann Stability Analysis of Globally Constraint-Preserving DGTD and PNPM Schemes for the Maxwell Equations using Multidimensional Riemann Solvers, Journal of Computational Physics, 376 (2019) 1108-1137] the authors presented globally constraint-preserving DGTD schemes for CED. The resulting schemes had excellent low dissipation and low dispersion properties. Their one deficiency was that the maximal permissible CFL of DGTD schemes decreased with increasing order of accuracy. The goal of this paper is to show how this deficiency is overcome. Because CED entails the propagation of electromagnetic waves, we would also like to obtain DG schemes for CED that minimize dissipation and dispersion errors even more than the prior generation of DGTD schemes. Two recent advances make this possible. The first advance, which has been reported elsewhere, is the development of a multidimensional Generalized Riemann Problem (GRP) solver. The second advance relates to the use of Two Derivative Runge Kutta (TDRK) timestepping. This timestepping uses not just the solution of the multidimensional Riemann problem, it also uses the solution of the multidimensional GRP. When these two advances are melded together, we arrive at DG(TD)2 schemes for CED, where the first “TD” stands for time-derivative and the second “TD” stands for the TDRK timestepping. The first goal of this paper is to show how DG(TD)2 schemes for CED can be formulated with the help of the multidimensional GRP and TDRK timestepping. The second goal of this paper is to utilize the free parameters in TDRK timestepping to arrive at DG(TD)2 schemes for CED that offer a uniformly large CFL with increasing order of accuracy while minimizing the dissipation and dispersion errors to exceptionally low values. The third goal of this paper is to document a von Neumann stability analysis of DG(TD)2 schemes so that their dissipation and dispersion properties can be quantified and studied in detail. At second order we find a DG(TD)2 scheme with CFL of 0.25 and improved dissipation and dispersion properties; for a second order scheme. At third order we present a novel DG(TD)2 scheme with CFL of 0.2571 and improved dissipation and dispersion properties; for a third order scheme. At fourth order we find a DG(TD)2 scheme with CFL of 0.2322 and improved dissipation and dispersion properties. As an extra benefit, the resulting DG(TD)2 schemes for CED require fewer synchronization steps on parallel supercomputers than comparable DGTD schemes for CED. We also document some test problems to show that the methods achieve their design accuracy.

This paper presents a topology optimization algorithm based on the lattice Boltzmann method coupled with a level-set method for increasing the efficiency of reactive fluid flows. The multi-relaxation time model is considered for the lattice Boltzmann collision operator, allowing higher Reynolds numbers flow simulations compared to the ordinary single-relaxation time model. The cost function gradient is obtained with the derivation of the adjoint-state formulation for the fully coupled problem. The proposed method is tested successfully on several numerical applications involving Reynolds numbers from 10 up to 1,000, as well as with different Damkohler and Peclet numbers. A limitation of the maximal pressure drop is also applied. The obtained results demonstrate that the proposed numerical method is robust and efficient for solving topology optimization problems of reactive fluid flows, in different operating conditions.

An arbitrary Lagrangian–Eulerian (ALE) finite element method for arbitrarily curved and deforming two-dimensional materials and interfaces is presented here. An ALE theory is developed by endowing the surface with a mesh whose in-plane velocity need not depend on the in-plane material velocity, and can be specified arbitrarily. A finite element implementation of the theory is formulated and applied to curved and deforming surfaces with in-plane incompressible flows. Numerical inf–sup instabilities associated with in-plane incompressibility are removed by locally projecting the surface tension onto a discontinuous space of piecewise linear functions. The general isoparametric finite element method, based on an arbitrary surface parametrization with curvilinear coordinates, is tested and validated against several numerical benchmarks. A new physical insight is obtained by applying the ALE developments to cylindrical fluid films, which are computationally and analytically found to be stable to non-axisymmetric perturbations, and unstable with respect to long-wavelength axisymmetric perturbations when their length exceeds their circumference. A Lagrangian scheme is attained as a special case of the ALE formulation. Though unable to model fluid films with sustained shear flows, the Lagrangian scheme is validated by reproducing the cylindrical instability. However, relative to the ALE results, the Lagrangian simulations are found to have spatially unresolved regions with few nodes, and thus larger errors.

In this paper, the efficiency in mesh updating (r-adaptivity) of the Transfinite Mean value Interpolation (TMI) and its generalization (k-TMI) are compared on three standardized problems to the well-known Inverse Distance Weighted interpolation (IDW) and Radial Basis Function interpolation (RBF) for unstructured data points and the new k-Transfinite Barycentric Interpolation (k-TBI) for structured data points such as, for instance, curves or surfaces in 3D. This is achieved by introducing a dynamical version of these interpolations via an ordinary differential equation that can be solved by standard ODE methods that are more economical than, for instance, solving vector partial differential equations as in the pseudo-solid method. A review of the very recent mathematical foundations of the k-TMI and k-TBI constructed from the function alone (standard) or from the function and its derivatives (enhanced) is provided in the first part of the paper.

We present a numerical method for interface-resolved simulations of evaporating two-fluid flows based on the volume-of-fluid (VoF) method. The method has been implemented in an efficient FFT-based two-fluid Navier-Stokes solver, using an algebraic VoF method for the interface representation, and extended with the transport equations of thermal energy and vaporized liquid mass for the single-component evaporating liquid in an inert gas. The conservation of vaporizing liquid and computation of the interfacial mass flux are performed with the aid of a reconstructed signed-distance field, which enables the use of well-established methods for phase change solvers based on level-set methods. The interface velocity is computed with a novel approach that ensures accurate mass conservation, by constructing a divergence-free extension of the liquid velocity field onto the entire domain. The resulting approach does not depend on the type of interface reconstruction (i.e. can be employed in both algebraic and geometrical VoF methods). We extensively verified and validated the overall method against several benchmark cases, and demonstrated its excellent mass conservation and good overall performance for simulating evaporating two-fluid flows in two and three dimensions.

This contribution is the numerically oriented companion article of the work [9]. We focus here on the numerical resolution of the embedded corrector problem introduced in [8], [9] in the context of homogenization of diffusion equations. Our approach consists in considering a corrector-type problem, posed on the whole space, but with a diffusion matrix which is constant outside some bounded domain. In [9], we have shown how to define three approximate homogenized diffusion coefficients on the basis of the embedded corrector problem. We have also proved that these approximations all converge to the exact homogenized coefficients when the size of the bounded domain increases. We show here that, under the assumption that the diffusion matrix is piecewise constant, the corrector problem to solve can be recast as an integral equation. In case of spherical inclusions with isotropic materials, we explain how to efficiently discretize this integral equation using spherical harmonics, and how to use the fast multipole method (FMM) to compute the resulting matrix-vector products at a cost which scales only linearly with respect to the number of inclusions. Numerical tests illustrate the performance of our approach in various settings.

In this work, a Lagrangian finite pointset model (FPM) is first developed to solve the 2D viscoelastic fluid governing equations, and then a corrected particle shifting technique (CPST) is introduced to eliminate the tensile instability in the long time simulations and added in the above FPM scheme (named as the FPMT). Subsequently, a coupled particle method (FPMT-SPH) is tentatively proposed to simulate the viscoelastic free surface flow, in which the FPMT method is used in the interior of fluid domain and the SPH method is adopted to treat the free surface near the boundary. The proposed FPMT method and FPMT-SPH method are motivated by: a) the spatial derivatives of the velocity and the viscoelastic stress are approximated and obtained by the weighted least squares method; b) the pressure is accurately solved by the application of projection method with a local iterative procedure; c) a corrected particle shifting technique is added to remedy the tensile instability (FPMT); d) an identification technique of free-surface particles is given in the FPMT-SPH method. The accuracy and the convergence of the proposed FPMT scheme for viscoelastic flow are first discussed by solving the planar flow based on the Oldroyd-B model, and compared with the analytical solutions. Secondly, the validity of the CPST is tested by several benchmarks, and compared with other numerical results. Thirdly, the viscoelastic lid-driven cavity flow at high Weissenberg number is simulated and compared with grid-based results, for further illustrating the robustness and the ability of the proposed FPMT method. Finally, the proposed coupled FPMT-SPH method is used to simulate the challenging free surface flow problem of a viscoelastic droplet impacting and spreading on rigid. All the numerical results show that the proposed particle method for the viscoelastic fluid or free surface flow is robust and reliable.

This work combines the lattice Boltzmann equation (LBE) and the overset method to simulate moving boundary problems in Navier-Stokes flows in two dimensions (2D). The transformation of the velocity moments of the distribution functions between a moving frame of reference and the one at rest is analyzed. The flow past a cylinder moving with a prescribed motion is used to validate the proposed LBE-overset method. We show that the proposed LBE-overset method does reduce the numerical noise generated by the relative motion between a moving object and the underlying Eulerian mesh for flow fields by several orders of magnitude.

Recent implementations of the Electric-Field Integral Equation (EFIE) for the electromagnetic scattering analysis of perfectly conducting targets rely on the electric current expansion with the monopolar-RWG basis functions, discontinuous across mesh edges, and the field testing over volumetric subdomains attached to the surface boundary triangulation. As compared to the standard RWG-based EFIE-approaches, normally continuous across edges, these schemes exhibit enhanced versatility, allowing the analysis of geometrically non-conformal meshes, and improved accuracy, especially for subwavelength sharp-edged conductors. In this paper, we present a monopolar-RWG discretization by the Method of Moments (MoM) of the Combined-Field Integral Equation (CFIE) resulting from the addition of a volumetrically tested discretization of the EFIE and the Galerkin tested MFIE-implementation. We show for sharp-edged conductors the degree of improved accuracy in the computed RCS and the convergence properties in the iterative search of the solution. More importantly, as we show in the paper, these implementations become in practice advantageous because of their robustness to flaws in the grid generation or their agility in handling complex meshes arising from the interconnection of independently meshed domains. The hybrid RWG/monopolar-RWG discretization of the CFIE defines the RWG discretization over geometrically conformal and smoothly varying mesh regions and inserts the monopolar-RWG expansion strictly at sharp edges, for improved accuracy purposes, or over boundary lines between partitioning mesh domains, for the sake of enhanced versatility. These hybrid schemes offer similar accuracy as their fully monopolar-RWG counterparts but with fewer unknowns and allow naturally non-conformal mesh transitions without inserting additional inter-domain continuity conditions or new artificial currents.

We present a novel methodology to compute relaxed dislocations core configurations, and their energies in crystalline metallic materials using large-scale ab-intio simulations. The approach is based on MacroDFT, a coarse-grained density functional theory method that accurately computes the electronic structure with sub-linear scaling resulting in a tremendous reduction in cost. Due to its implementation in real-space, MacroDFT has the ability to harness petascale resources to study materials and alloys through accurate ab-initio calculations. Thus, the proposed methodology can be used to investigate dislocation cores and other defects where long range elastic defects play an important role, such as in dislocation cores, grain boundaries and near precipitates in crystalline materials. We demonstrate the method by computing the relaxed dislocation cores in prismatic dislocation loops and dislocation segments in magnesium (Mg). We also study the interaction energy with a line of Aluminum (Al) solutes. Our simulations elucidate the essential coupling between the quantum mechanical aspects of the dislocation core and the long range elastic fields that they generate. In particular, our quantum mechanical simulations are able to describe the logarithmic divergence of the energy in the far field as is known from classical elastic theory. In order to reach such scaling, the number of atoms in the simulation cell has to be exceedingly large, and cannot be achieved with the state-of-the-art density functional theory implementations.

We propose a multilevel approach for trace systems resulting from hybridized discontinuous Galerkin (HDG) methods. The key is to blend ideas from nested dissection, domain decomposition, and high-order characteristic of HDG discretizations. Specifically, we first create a coarse solver by eliminating and/or limiting the front growth in nested dissection. This is accomplished by projecting the trace data into a sequence of same or high-order polynomials on a set of increasingly h−coarser edges/faces. We then combine the coarse solver with a block-Jacobi fine scale solver to form a two-level solver/preconditioner. Numerical experiments indicate that the performance of the resulting two-level solver/preconditioner depends on the smoothness of the solution and can offer significant speedups and memory savings compared to the nested dissection direct solver. While the proposed algorithms are developed within the HDG framework, they are applicable to other hybrid(ized) high-order finite element methods. Moreover, we show that our multilevel algorithms can be interpreted as a multigrid method with specific intergrid transfer and smoothing operators. With several numerical examples from Poisson, pure transport, and convection-diffusion equations we demonstrate the robustness and scalability of the algorithms with respect to solution order. While scalability with mesh size in general is not guaranteed and depends on the smoothness of the solution and the type of equation, improving it is a part of future work.

A design tool was formulated for optimizing the efficiency of inorganic, thin-film, photovoltaic solar cells. The solar cell can have multiple semiconductor layers in addition to antireflection coatings, passivation layers, and buffer layers. The solar cell is backed by a metallic grating which is periodic along a fixed direction. The rigorous coupled-wave approach is used to calculate the electron-hole-pair generation rate. The hybridizable discontinuous Galerkin method is used to solve the drift-diffusion equations that govern charge-carrier transport in the semiconductor layers. The chief output is the solar-cell efficiency which is maximized using the differential evolution algorithm to determine the optimal dimensions and bandgaps of the semiconductor layers.

We construct a finite volume element scheme with a monotonicity correction for anisotropic diffusion problems on general quadrilateral meshes, where the strict convexity restriction on the meshes is removed. The main contributions of this paper include three aspects. Firstly, the classical finite volume element (C-FVE) method is extended to severely distorted quadrilateral meshes even with concave cells by virtue of a new overlapping dual partition and a new gradient approximation. In fact, the choice of this dual partition and gradient approximation is also conducive to the construction of a monotone scheme. Secondly, a new monotonicity correction is suggested, based on which we obtain a monotone finite volume element (M-FVE) method. The resulting M-FVE method still keeps the local conservation and is easy for implementation. Finally, we analyze theoretically the truncation error and the monotonicity for this scheme. Besides, the existence of a solution to this nonlinear scheme is proved by applying the Brouwer's fixed point theorem. Numerical results demonstrate that the M-FVE method has the approximate second-order accuracy and preserves well the positivity of the solution for both isotropic and anisotropic diffusion problems on severely distorted quadrilateral meshes.

The agglomeration of particles during the handling of powders results in caking, lumping or the local accumulation of electrostatic energy which represents a serious hazard to the operational safety of industrial facilities. In the case of dry powders the attraction in-between particles can be mainly attributed to van der Waals and electrostatic forces. Nonetheless, due to the challenges related to the small size and distance of relevant particles and the optical density of powder flows the detailed physical mechanisms of their interaction are so far little investigated. In this paper we present a novel numerical approach which is based on an algorithm developed by Erleben [1] in the field of computer graphics. This algorithm is extended to compute binary and multiple particle interaction with each other and solid surfaces. Therein, besides van der Waals and electrostatic forces also collisional forces and plastic particle deformation is accounted for. The herein presented results demonstrate that this algorithm allows to predict accurately and efficiently whether particles agglomerate or separate depending on their kinetic parameters. In particular, the imposed constrain forces prevent spurious velocity fluctuations and potential particle overlapping in statically overdetermined systems. Simulated test cases reveal how electrostatic and van der Waals forces lead to the growth of structures in case the particle restitution coefficient is sufficiently low.

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.

In a previous work that investigated atomic collisional-radiative models, Boltzmann grouping was implemented to reduce and accelerate detailed-configuration-accounting or DCA-based CR simulations, exhibiting distinct advantages over other CR reduction techniques. However, the selection of level groups under this reduction technique was manually performed and required several iterations to construct the appropriate groups that properly showed the advantages associated with Boltzmann grouping. Therefore, a clustering technique was implemented to automatically group detailed (LS-coupled) atomic states with limited user interference. Clustering in conjunction with the Boltzmann grouping technique were applied to collisional and coronal simulations in this work to demonstrate the feasibility of automatic level grouping techniques to produce accurate reduced simulations under various plasma conditions.

Previous works have developed boundary conditions that lead to the L2-boundedness of solutions to the linearised moment equations. Here we present a spatial discretization that preserves the L2-stability by recovering integration-by-parts over the discretized domain and by imposing boundary conditions using a simultaneous-approximation-term (SAT). We develop three different forms of the SAT using: (i) characteristic splitting of moment equation's boundary conditions; (ii) decoupling of moments in moment equations; and (iii) characteristic splitting of Boltzmann equation's boundary conditions. We discuss how the first two forms differ in terms of their usage and implementation. We show that the third form is equivalent to using an upwind kinetic numerical flux along the boundary, and we argue that even though it provides stability, it prescribes the incorrect number of boundary conditions. Using benchmark problems, we compare the accuracy of moment solutions computed using different SATs. Our numerical experiments also provide new insights into the convergence of moment approximations to the Boltzmann equation's solution.

We describe an Arbitrary-Lagrangian-Eulerian (ALE) method for the compressible Euler system with capillary force. The algorithm is split in two steps. First, the Lagrangian step is based on cell-centered schemes [9], [20], [46]. The surface tension force is discretized in order to exactly verify the Laplace law at the discrete level. We also provide a second-order spatial extension and a low-Mach correction, which do not break the well-balanced property of the scheme. The Lagrangian scheme is assessed on several problems, particularly on a linear Richtmyer-Meshkov instability which is the targeted application. The second step is the rezoning and remapping done thanks to a swept-region method using exact intersections near the interface. We use a Volume Of Fluid (VOF) method to track the interface. We describe the treatment of mixed-cells, and in particular the thermodynamics closure and the curvature calculation. The new scheme is used to investigate the influence of surface tension on a non-linear Richtmyer-Meshkov instability.

This paper presents an accuracy-preserving p-weighted limiter for discontinuous Galerkin methods on one-dimensional and two-dimensional triangular grids. The p-weighted limiter is the extension of the second-order WENO limiter by Li et al. [W. Li, J. Pan and Y-X Ren, Journal of Computational Physics 364(2018)314-346] to high-order accuracy, with the following important improvements of the limiting procedure. First, the candidate polynomials of the p-weighted limiter are the p-hierarchical orthogonal polynomials of the current cell, and the linear polynomials constructed by minimizing the projection error on the face-neighboring cells. Second, the p-weighted procedure introduces a new smoothness indicator which has less numerical dissipation comparing with the classical WENO one. The smoothness indicator is efficiently computed through a quadrature-free approach that takes advantage of the orthogonal property of the basis functions. Third, the small positive number ϵ, which is introduced in the weights to avoid dividing by zero, is set as a function of the smoothness indicator to preserve accuracy near smooth extremas. Numerous benchmark problems are solved to test the p1, p3 and p5 discontinuous Galerkin schemes using the p-weighted limiter. Numerical results demonstrate that the p-weighted limiter is capable of capturing strong shocks while preserving accuracy in smooth regions.

We have found provably optimal algorithms for full-domain discrete-ordinate transport sweeps on a class of grids in 2D and 3D Cartesian geometry that are regular at a coarse level but arbitrary within the coarse blocks. We describe these algorithms and show that they always execute the full eight-octant (or four-quadrant if 2D) sweep in the minimum possible number of stages for a given Px×Py×Pz partitioning. Computational results confirm that our optimal scheduling algorithms execute sweeps in the minimum possible stage count. Observed parallel efficiencies agree well with our performance model. Our PDT transport code has achieved approximately 68% parallel efficiency with >1.5M parallel threads, relative to 8 threads, on a simple weak-scaling problem with only three energy groups, 10 direction per octant, and 4096 cells/thread. Our ARDRA code has achieved 71% efficiency with >1.5M cores, relative to 16 cores, with 36 directions per octant and 48 energy groups. We demonstrate similar efficiencies with PDT on a realistic set of nuclear-reactor test problems, with unstructured meshes that resolve fine geometric details. These results demonstrate that discrete-ordinates transport sweeps can be executed with high efficiency using more than 106 parallel processes.

The prediction of long-term dynamics of transitional environments, e.g., lagoon evolution, salt-marsh growth or river delta progradation, is an important issue to estimate the potential impacts of different scenarios on such vulnerable intertidal morphologies. The numerical simulation of the combined accretion and consolidation, i.e., the two main processes driving the dynamics of these environments, however, suffers from a significant geometric non-linearity, which may result in a pronounced grid distortion using standard grid-based discretization methods. The present work describes a novel numerical approach, based on the Virtual Element Method (VEM), for the long-term simulation of the vertical dynamics of transitional landforms. The VEM is a grid-based variational technique for the numerical discretization of Partial Differential Equations (PDEs) allowing for the use of very irregular meshes consisting of a free combination of different polyhedral elements. The model solves the groundwater flow equation, coupled to a geomechanical module based on Terzaghi's principle, in a large-deformation setting, taking into account both the geometric and the material non-linearity. The use of the VEM allows for a great flexibility in the element generation and management, avoiding the numerical issues connected with the adoption of strongly distorted meshes. The numerical model is developed, implemented and tested in real-world examples, showing an interesting potential for addressing complex environmental situations.

The objective of this work is to develop a procedure that allows for reconstructing three-dimensional flow fields from two-dimensional information contained in some representative planes, conveniently distributed throughout the domain. The reconstructing tool is based on two recent methods developed by two of the authors, namely the higher order dynamic mode decomposition and the spatio-temporal Koopman decomposition (STKD). The latter method decomposes a given spatio-temporal flow field as a series expansion in Fourier-like modes (including both wavenumbers/frequencies and spatial/temporal growth rates) in time and some distinguished longitudinal spatial coordinates, and spatial modes depending on the remaining transverse spatial coordinates. To obtain the (unknown) three-dimensional reconstruction, the STKD method is first applied to the (known) two-dimensional data in the considered planes. This application of STKD yields both the wavenumbers/frequencies and the spatial/temporal growth rates appearing in the three-dimensional reconstruction. Imposing that the three-dimensional reconstruction coincides with the two-dimensional data in the given set of planes, a system of linear equations results that permits computing the various ingredients appearing in the three-dimensional STKD expansion. The performance of the method is tested in the three-dimensional wake of a circular cylinder at Reynolds number (based on the cylinder diameter and the incoming free-stream velocity) equal to 280. The resulting flow is highly non-linear and quasi-periodic, with two fundamental temporal frequencies, associated with the well-known modes A and B. The method can be applied using experimental or numerical data, allowing to identify the full three-dimensional flow field and its main characteristics from limited information.

We are concerned with the approximation of solutions to a compressible two-phase flow model in porous media thanks to an enhanced control volume finite element discretization. The originality of the methodology consists in treating the case where the densities are depending on their own pressures without any major restriction neither on the permeability tensor nor on the mesh. Contrary to the ideas of [23], the point of the current scheme relies on a phase-by-phase “sub”-unpwinding approach so that we can recover the coercivity-like property. It allows on a second place for the preservation of the physical bounds on the discrete saturation. The convergence of the numerical scheme is therefore performed using classical compactness arguments. Numerical experiments are presented to exhibit the efficiency and illustrate the qualitative behavior of the implemented method.

The development of a novel Adaptive Mesh Enlargement (AME) method with ordinary Weight Essentially Non-Oscillatory (WENO) scheme for large-scale explosion problems is presented. The novel AME method adaptively adjusts the size of computational domains to accommodate the development of an explosion event which initially starts in a small domain. The grid resolution is verified by solving a simple one-dimensional (1-D) partial differential equation. The accuracy of AME methods is evaluated by calculations of 1-D Burgers equations and shock tube tests. It is found that the accuracy of AME methods depends on the number of enlargement operations. Numerical solutions with a reasonable number of enlargements are consistent with those of the fixed mesh method. At the same time, AME methods with multiple enlargements reduce the computational cost by several orders of magnitude. A three-dimensional (3-D) symmetric model with the AME method operating five enlargements is further developed to simulate the large-scale near-field explosion and validated against experiments. Propagation of the rarefaction and shock waves as well as the reflected wave are analyzed. The predicted peak incident pressures yield a variation of ∼3% compared to experiments indicating that the AME method is an effective construction method of computational domains for the problem of large-scale fluid dynamics.

We consider the interface advection problem by a prescribed velocity field in the special case when the interface intersects the domain boundary, i.e. in the presence of a contact line. This problem emerges from the discretization of continuum models for dynamic wetting. The kinematic evolution equation for the dynamic contact angle (Fricke et al., 2019) expresses the fundamental relationship between the rate of change of the contact angle and the structure of the transporting velocity field. The goal of the present work is to develop an interface advection method that is consistent with the fundamental kinematics and transports the contact angle correctly with respect to a prescribed velocity field. In order to verify the advection method, the kinematic evolution equation is solved numerically and analytically (for special cases). We employ the geometrical Volume-of-Fluid (VOF) method on a structured Cartesian grid to solve the hyperbolic transport equation for the interface in two spatial dimensions. We introduce generalizations of the Youngs and ELVIRA methods to reconstruct the interface close to the domain boundary. Both methods deliver first-order convergent results for the motion of the contact line. However, the Boundary Youngs method shows strong oscillations in the numerical contact angle that do not converge with mesh refinement. In contrast to that, the Boundary ELVIRA method provides linear convergence of the numerical contact angle transport.

This paper presents a coupling method of the level set and volume of fluid methods based on a simple local-gradient based re-initialisation approach that evaluates the distance function depending on the computational cell location. If a cell belongs to the interface, the signed distance is updated based on a search in the neighbouring cells and an interpolation procedure is applied depending on the local curvature or the sign of the level set function following [D. Hartmann, M. Meinke, W. Schröder, Differential equation based constrained reinitialisation method for level set methods, J. Comput. Phys. 227 (2008) 6821-6845]. The search algorithm does not distinguish between the upwind and downwind directions and hence it is able to be used for cells with an arbitrary number of faces increasing the robustness of the method. The coupling with the volume of fluid method is achieved by mapping the volume fraction field which is advected from the isoface evolution at a subgrid level. Consequently, the coupling with the level set approach is utilised without solving the level set equation. This coupled method provides better accuracy than the volume of fluid method alone and is capable of capturing sharp interfaces in all the classical numerical tests that are presented here.

We study the interior transmission eigenvalue problem for the elastic wave scattering in this paper. We aim to show the distribution of positive eigenvalues by efficient numerical algorithms. Here the elastic waves are scattered by the perturbations of medium parameters, which include the elasticity tensor C and the density ρ. Let us denote (C0,ρ0) and (C1,ρ1) the background and the perturbed medium parameters, respectively. We consider two cases of perturbations, C0=C1,ρ1≠ρ0 (case 1) and C0≠C1,ρ1=ρ0 (case 2). After discretizing the associated PDEs by FEM, we are facing the computation of generalized eigenvalues problems (GEP) with matrices of large size. These GEPs contain huge number of nonphysical zeros (for case 1) or nonphysical infinities (for case 2). In order to locate several hundred positive eigenvalues effectively, we then convert GEPs to suitable quadratic eigenvalues problems (QEP). We then implement a quadratic Jacobi-Davidson method combining with partial locking or partial deflation techniques to compute 500 positive eigenvalues.

Proper orthogonal decomposition (POD) and dynamic mode decomposition (DMD) are two complementary singular-value decomposition (SVD) techniques that are widely used to construct reduced-order models (ROMs) in a variety of fields of science and engineering. Despite their popularity, both DMD and POD struggle to formulate accurate ROMs for advection-dominated problems because of the nature of SVD-based methods. We investigate this shortcoming of conventional POD and DMD methods formulated within the Eulerian framework. Then we propose a Lagrangian-based DMD method to overcome this so-called translational problem. Our approach is consistent with the spirit of physics-aware DMD since it accounts for the evolution of characteristic lines. Several numerical tests are presented to demonstrate the accuracy and efficiency of the proposed Lagrangian DMD method.

This work introduces a new type of constrained algebraic stabilization for continuous piecewise-linear finite element approximations to the equations of ideal magnetohydrodynamics (MHD). At the first step of the proposed flux-corrected transport (FCT) algorithm, the Galerkin element matrices are modified by adding graph viscosity proportional to the fastest characteristic wave speed. At the second step, limited antidiffusive corrections are applied and divergence cleaning is performed for the magnetic field. The limiting procedure developed for this stage is designed to enforce local maximum principles, as well as positivity preservation for the density and thermodynamic pressure. Additionally, it adjusts the magnetic field in a way which penalizes divergence errors without violating conservation laws or positivity constraints. Numerical studies for 2D test problems are performed to demonstrate the ability of the proposed algorithms to accomplish this task in applications to ideal MHD benchmarks.

The radiation diffusion problem is a kind of time-dependent nonlinear equations. For solving the radiation diffusion equations, many linear systems are obtained in the nonlinear iterations at each time step. The cost of linear equations dominates the numerical simulation of radiation diffusion applications, such as inertial confinement fusion, etc. Usually, iterative methods are used to solve the linear systems in a real application. Moreover, the solution of the previous nonlinear iteration or the solution of the previous time step is typically used as the initial guess for solving the current linear equations. Because of the strong local character in ICF, with the advancing of nonlinear iteration and time step, the solution of the linear system changes dramatically in some local domain, and changes mildly or even has no change in the rest domain. In this paper, a local character-based method is proposed to solve the linear systems of radiation diffusion problems. The proposed method consists of three steps: firstly, a local domain (algebraic domain) is constructed; secondly, the subsystem on the local domain is solved; and lastly, the whole system will be solved. Two methods are given to construct the local domain. One is based on the spatial gradient, and the other is based on the residual. Numerical tests for a two-dimensional heat conduction model problem, and two real application models, the multi-group radiation diffusion equations and the three temperature energy equations, are conducted. The test results show that the solution time for solving the linear system can be reduced dramatically by using the local character-based method.

In this paper we consider scattering resonance computations in optics when the resonators consist of frequency dependent and lossy materials, such as metals at optical frequencies. The proposed computational approach combines a novel hp-FEM strategy, based on dispersion analysis for complex frequencies, with a fast implementation of the nonlinear eigenvalue solver NLEIGS. Numerical computations illustrate that the pre-asymptotic phase is significantly reduced compared to standard uniform h and p strategies. Moreover, the efficiency grows with the refractive index contrast, which makes the new strategy highly attractive for metal-dielectric structures. The hp-refinement strategy together with the efficient parallel code result in highly accurate approximations and short runtimes on multi processor platforms.

The immersed boundary–finite element method (IBFE) is an approach to describing the dynamics of an elastic structure immersed in an incompressible viscous fluid. In this formulation, there are discontinuities in the pressure and viscous stress at fluid–structure interfaces. The standard immersed boundary approach, which connects the Lagrangian and Eulerian variables via integral transforms with regularized Dirac delta function kernels, smooths out these discontinuities, which generally leads to low order accuracy. This paper describes an approach to accurately resolve pressure discontinuities for these types of formulations, in which the solid may undergo large deformations. Our strategy is to decompose the physical pressure field into a sum of two pressure–like fields, one defined on the entire computational domain, which includes both the fluid and solid subregions, and one defined only on the solid subregion. Each of these fields is continuous on its domain of definition, which enables high accuracy via standard discretization methods without sacrificing sharp resolution of the pressure discontinuity. Numerical tests demonstrate that this method improves rates of convergence for displacements, velocities, stresses, and pressures, as compared to the conventional IBFE method. Further, it produces much smaller errors at reasonable numbers of degrees of freedom. The performance of this method is tested on several cases with analytic solutions, a nontrivial benchmark problem of incompressible solid mechanics, and an example involving a thick, actively contracting torus.

In this work, we consider a mathematical model for describing flow in an unsaturated porous medium containing a fracture. Both the flow in the fracture as well as in the matrix blocks are governed by Richards' equation coupled by natural transmission conditions. Using formal asymptotics, we derive upscaled models as the limit of vanishing ε, the ratio of the width and length of the fracture. Our results show that the ratio of porosities and permeabilities in the fracture to matrix determine, to the leading order of approximation, the appropriate effective model. In these models the fracture is a lower dimensional object for which different transversally averaged models are derived depending on the ratio of the porosities and permeabilities of the fracture and respective matrix blocks. We obtain a catalogue of effective models which are validated by numerical computations.

A realistic inflow boundary condition is essential for successful simulation of the developing turbulent boundary layer or channel flows. In the present work, we applied generative adversarial networks (GANs), a representative of unsupervised learning, to generate an inlet boundary condition of turbulent channel flow. Upon learning the two-dimensional spatial structure of turbulence using data obtained from direct numerical simulation (DNS) of turbulent channel flow, the GAN could generate instantaneous flow fields that are statistically similar to those of DNS. After learning data at only three Reynolds numbers, the GAN could produce fields at various Reynolds numbers within a certain range without additional simulation. Eventually, through a combination of the GAN and a recurrent neural network (RNN), we developed a novel model (RNN-GAN) that could generate time-varying fully developed flow for a long time. The spatiotemporal correlations of the generated flow are in good agreement with those of the DNS. This proves the usefulness of unsupervised learning in the generation of synthetic turbulence fields.

This paper examines the Laplace equation with mixed boundary conditions, the Neumann and Steklov boundary conditions. This models a container with holes in it, like a pond filled with water but partly covered by immovable pieces on the surface. The main objective is to determine the right extent of the covering pieces, so that any shock inside the container yields a resonance. To this end, an algorithm is developed which uses asymptotic formulas concerning perturbations of the partitioning of the boundary pieces. Proofs for these formulas are established. Furthermore, this paper displays some results concerning bounds and examples with regards to the governing problem.

In this study we are aimed to simulate incompressible fluid flows with moving interface separating liquid and gas phases. To achieve volume/mass conservation and to capture interface, the interface-capturing volume of fluid (VOF) method, which is exercised in compliance with requirement of volume/mass conservation, will be coupled with the other interface-capturing level-set (LS) method, which is suitable to capture interface accurately. In our proposed advection algorithm, VOF is the building block that solves the volume fraction and the level-set function is solely used to assist an accurate calculation of some geometrically relevant quantities at the interface. A high order scheme developed within the optimized compact reconstruction WENO framework has been applied to solve the advection equation. The novelty of this purposed advection algorithm is attributed to its efficient implementation without sacrifice of computational accuracy.

The modeling of atmospheric processes in the context of weather and climate simulations is an important and computationally expensive challenge. The temporal integration of the underlying PDEs requires a very large number of time steps, even when the terms accounting for the propagation of fast atmospheric waves are treated implicitly. Therefore, the use of parallel-in-time integration schemes to reduce the time-to-solution is of increasing interest, particularly in the numerical weather forecasting field. We present a multi-level parallel-in-time integration method combining the Parallel Full Approximation Scheme in Space and Time (PFASST) with a spatial discretization based on Spherical Harmonics (SH). The iterative algorithm computes multiple time steps concurrently by interweaving parallel high-order fine corrections and serial corrections performed on a coarsened problem. To do that, we design a methodology relying on the spectral basis of the SH to coarsen and interpolate the problem in space. The methods are evaluated on the shallow-water equations on the sphere using a set of tests commonly used in the atmospheric flow community. We assess the convergence of PFASST-SH upon refinement in time. We also investigate the impact of the coarsening strategy on the accuracy of the scheme, and specifically on its ability to capture the high-frequency modes accumulating in the solution. Finally, we study the computational cost of PFASST-SH to demonstrate that our scheme resolves the main features of the solution multiple times faster than the serial schemes.

In this paper, we develop high-order bound-preserving (BP) finite difference (FD) methods for the coupled system of compressible miscible displacements. We consider the problem with multi-component fluid mixture and the (volumetric) concentration of the jth component, cj, should be between 0 and 1. It is well known that cj does not satisfy a maximum-principle. Hence most of the existing BP techniques cannot be applied directly. The main idea in this paper is to construct the positivity-preserving techniques to all cj′s and enforce ∑jcj=1 simultaneously to obtain physically relevant approximations. By doing so, we have to treat the time derivative of the pressure dp/dt as a source in the concentration equation and choose suitable “consistent” numerical fluxes in the pressure and concentration equations. Recently, the high-order BP discontinuous Galerkin (DG) methods for miscible displacements were introduced in [6]. However, the BP technique for DG methods is not straightforward extendable to high-order FD schemes. There are two main difficulties. Firstly, it is not easy to determine the time step size in the BP technique. In finite difference schemes, we need to choose suitable time step size first and then apply the flux limiter to the numerical fluxes. Subsequently, we can compute the source term in the concentration equation, leading to a new time step constraint that may not be satisfied by the time step size applied in the flux limiter. Therefore, it would be very difficult to determine how large the time step is. Secondly, the general treatment for the diffusion term, e.g. centered difference, in miscible displacements may require a stencil whose size is larger than that for the convection term. It would be better to construct a new spatial discretization for the diffusion term such that a smaller stencil can be used. In this paper, we will solve both problems. We first construct a special discretization of the convection term, which yields the desired approximations of the source. Then we can find out the time step size that suitable for the BP technique and apply the flux limiters. Moreover, we will also construct a special algorithm for the diffusion term whose stencil is the same as that used for the convection term. Numerical experiments will be given to demonstrate the high-order accuracy and good performance of the numerical technique.

Continuum kinetic simulations of plasmas, where particle distribution functions are directly discretized in phase-space, permit fully kinetic simulations without the statistical noise of particle-in-cell methods. Recent advances in numerical algorithms have made continuum kinetic simulations computationally competitive. This work presents a continuum kinetic description of high-fidelity wall boundary conditions that utilize the readily available particle distribution function without coupling to additional physical models. The boundary condition is realized through a reflection function that can capture a wide range of cases from simple specular reflection to more involved first principles models. While the framework is usable for various numerical methods and boundary conditions, this work focuses on the discontinuous Galerkin implementation of electron emission using a first-principles quantum-mechanical model. Presented results demonstrate effects of electron emission from a dielectric material on formation of a classical plasma sheath.

A new method is proposed to solve Ampere's equation in an arbitrary toroidal domain in which all currents are known, given proper boundary conditions for the magnetic vector potential. The novelty of the approach lies in the application of singular value decomposition (SVD) techniques to tackle the difficulties caused by the kernel associated by the curl operator. This kernel originates physically due to the magnetic field gauge. To increase the efficiency of the solver, the problem is represented by means of a dual finite difference-spectral scheme in arbitrary generalized toroidal coordinates, which permits to take advantage of the block structure exhibited by the matrices that describe the discretized problem. The result is a fast and efficient solver, up to three times faster than the double-curl method in some cases, that provides an accurate solution of the differential form of Ampere law while guaranteeing a zero divergence of the resulting magnetic field down to machine precision.

We propose a two-layer model with two different axes of integration and a well-balanced finite volume method. The purpose is to study submarine avalanches and generated tsunamis by a depth-averaged model with different averaged directions for the fluid and the granular layers. Two-layer shallow depth-averaged models usually consider either Cartesian or local coordinates for both layers. However, the motion characteristics of the granular layer and the water wave are different: the granular flow velocity is mainly oriented downslope while water motion related to tsunami wave propagation is mostly horizontal. As a result, the shallow approximation and depth-averaging have to be imposed (i) in the direction normal to the topography for the granular flow and (ii) in the vertical direction for the water layer. To deal with this problem, we define a reference plane related to topography variations and use the associated local coordinates to derive the granular layer equations whereas Cartesian coordinates are used for the fluid layer. Depth-averaging is done orthogonally to that reference plane for the granular layer equations and in the vertical direction for the fluid layer equations. Then, a finite volume method is defined based on an extension of the hydrostatic reconstruction. The proposed method is exactly well-balanced for two kind of stationary solutions: the classical one, when both water and granular masses are at rest; the second one, when only the granular mass is at rest. Several tests are presented to get insight into the sensitivity of the granular flow, deposit and generated water waves to the choice of the coordinate systems. Our results show that even for moderate slopes (up to 30∘), strong relative errors on the avalanche dynamics and deposit (up to 60%) and on the generated water waves (up to 120%) are made when using Cartesian coordinates for both layers instead of an appropriate local coordinate system as proposed here.

This paper describes an innovative technique that is capable of automatically detecting shock-wave shapes and shock-interaction patterns by post-processing the results obtained from a shock-capturing calculation using any of the available, state-of-the-art, shock-detection techniques. This new technique makes use of the Hough transform, the least-squares fit and the fuzzy logic to post-process the data supplied by classical shock-detection techniques; it not only allows to significantly enhance the quality of the extracted shock-fronts, but, even more importantly, allows to recognize the different types of shock-shock and shock-wall interactions that may occur within the flowfield.

In this paper we propose a space-time adaptive finite element method for the phase field model of pitting corrosion, which is a parabolic partial differential equation system consisting of a phase variable and a concentration variable. A major challenge in solving this phase field model is that the problem is very stiff, which makes the time step size extremely small for standard temporal discretizations. Another difficulty is that a high spatial resolution is required to capture the steep gradients within the diffused interface, which results in very large number of degrees of freedom for uniform meshes. To overcome the stiffness of this model, we combine the Rosenbrock–Euler exponential integrator with Crank–Nicolson scheme for the temporal discretization. Moreover, by exploiting the fact that the speed of the corroding interface decreases with time, we derive an adaptive time stepping formula. For the spatial approximation, we propose a simple and efficient strategy to generate adaptive meshes that reduce the computational cost significantly. Thus, the proposed method utilizes local adaptivity and mesh refinement for efficient simulation of the corrosive dissolution over long times in heterogeneous media with complex microstructures. We also present an extensive set of numerical experiments in both two and three dimensional spaces to demonstrate efficiency and robustness of the proposed method.

Many phenomena, ranging from biology to electronics, involve flow over rough or irregular surfaces. We treat such surfaces as random fields and use an immersed boundary method (IBM) with discrete (random) forcing to solve resulting stochastic flow problems. Our approach relies on the Uhlmann formulation of the fluid-solid interaction force; computational savings stem from the modification of the time advancement scheme that obviates the need to solve the Poisson equation for pressure at each sub-step. We start by testing the proposed algorithm on two classical benchmark problems. The first deals with the Wannier problem of Stokesian flow around a cylinder in the vicinity of a moving plate. The second problem considers steady-state and transient flows over a stationary circular cylinder. Our simulation results show that our algorithm achieves second-order temporal accuracy. It is faster than the original IBM, while yielding consistent estimates of such quantities of interest as the drag and lift coefficients, the length of a recirculation zone in a cylinder's wake, and the Strouhal number. Then we use the proposed IBM algorithm to model flow over cylinders whose surface is either (deterministically) corrugated or (randomly) rough.

A multi-physics numerical model for laser-induced optical breakdown and laser-plasma interaction in a non-equilibrium gas is presented, accounting for: production of priming electrons via multi-photon ionization, energy absorption, cascade ionization, induced hydrodynamic response, and shock formation and propagation. The gas is governed by the Navier-Stokes equations, with non-equilibrium effects taken into account by means of a two-temperature model. The space-time dependence of the laser beam is modeled with a flux-tube formulation for the Radiative Transfer Equation. The flow governing equations are discretized in space using a second-order finite volume method. The semi-discrete equations are marched in time using an implicit-explicit (IMEX) dual time-stepping strategy, where diffusion and chemistry are solved implicitly, whereas convection is explicit. Application to a 20 ns long 50 mJ pulse laser-induced breakdown in quiescent O2 shows the advantages of this temporal discretization during and just after the laser pulse, while a less-expensive symmetric Strang splitting (for implicit chemistry) is sufficient for the post-breakdown gas dynamics after ≃ 0.1textmu s. The integrated model is shown to reproduce key features of corresponding experiments.

A robust numerical methodology to predict equilibrium interfaces over arbitrary solid surfaces is developed. The kernel of the proposed method is the distance regularized level set equations (DRLSE) with techniques to incorporate the no-penetration and volume-conservation constraints. In this framework, we avoid reinitialization that is typically used in traditional level set methods. This allows for a more efficient algorithm since only one advection equation is solved, and avoids numerical error associated with the re-distancing step. A novel surface tension distribution, based on harmonic mean, is prescribed such that the zero level set has the correct liquid-solid surface tension value. This leads to a more accurate triple contact point location. The method uses second-order central difference schemes which facilitates easy parallel implementation, and is validated by comparing to traditional level set methods for canonical problems. The application of the method in the context of Gibbs free energy minimization, to obtain liquid-air interfaces is validated against existing analytical solutions. The capability of the methodology to predict equilibrium shapes over both structured and realistic rough surfaces is demonstrated.

In the previous article Abushaikha et al. (2017) [1], we presented a fully-implicit mixed hybrid finite element (MHFE) method for general-purpose compositional reservoir simulation. The present work extends the implementation for mimetic finite difference (MFD) discretization method. The new approach admits fully implicit solution on general polyhedral grids. The scheme couples the momentum and mass balance equations to assure conservation and applies a cubic equation-of-state for the fluid system. The flux conservativity is strongly imposed for the fully implicit approach and the Newton-Raphson method is used to linearize the system. We test the method through extensive numerical examples to demonstrate the convergence and accuracy on various shapes of polyhedral. We also compare the method to other discretization schemes for unstructured meshes and tensor permeability. Finally, we apply the method through applied computational cases to illustrate its robustness for full tensor anisotropic, highly heterogeneous and faulted reservoirs using unstructured grids.

In the present work, the Cahn-Hilliard Phase-Field model of incompressible two-phase flows is considered. Conditions needed for consistency of reduction, consistency of mass and momentum transport, and consistency of mass conservation are proposed. The mass flux in the Navier-Stokes equations is defined such that it satisfies the proposed consistency conditions. The analysis in both continuous and discrete levels shows that violation of the consistency conditions result in unphysical solutions and the inconsistent errors are proportional to the density contrast of the fluids. After considering the conservative form of the inertial term, a consistent and conservative scheme for momentum transport is developed. The balanced-force algorithm for the sharp interface model is extended to the surface force derived from the Cahn-Hilliard model. The proposed scheme is formally 2nd-order accurate in both time and space, satisfies the consistency conditions, conserves mass globally and momentum essentially, and is balanced-force, in the discrete level. Its convergence to the sharp interface solution is systematically discussed in cases including large density and viscosity ratios, surface tension, and gravity. Various two-phase flow problems with large density ratios are performed to validate and verify the proposed scheme and excellent agreements with published numerical and/or experimental results are achieved. The proposed scheme is a practical and accurate tool to study two-phase flows, especially for those including large density ratios.

In this paper, a family of stencil selection algorithms is presented for WENO schemes on unstructured meshes. The associated freedom of stencil selection for unstructured meshes, in the context of WENO schemes present a plethora of various stencil selection algorithms. The particular focus of this paper is to assess the performance of various stencil selection algorithm, investigate the parameters that dictate their robustness, accuracy and computational efficiency. Ultimately, efficient and robust stencils are pursued that can provide significant savings in computational performance, while retaining the non-oscillatory character of WENO schemes. This is achieved when making the stencil selection algorithms adaptive, based on the quality of the cells for unstructured meshes, that can in turn reduce the computational cost of WENO schemes. For assessing the performance of the developed algorithms well established test problems are employed. These include the least square approximation of polynomial functions, linear advection equation of smooth functions and solid body rotation test problem. Euler and Navier-Stokes equations test problems are also pursued such as the Shu-Osher test problem, the Double Mach Reflection, the supersonic Forward Facing step, the Kelvin-Helmholtz instability, the Taylor-Green Vortex, and the flow past a transonic circular cylinder.

We describe a new and simple strategy based on the Gauss divergence theorem for obtaining centroidal gradients on unstructured meshes. Unlike the standard Green–Gauss (SGG) reconstruction which requires face values of quantities whose gradients are sought, the proposed approach reconstructs the gradients using the normal derivative(s) at the faces. The new strategy, referred to as the Modified Green–Gauss (MGG) reconstruction results in consistent gradients which are at least first-order accurate on arbitrary polygonal meshes. We show that the MGG reconstruction is linearity preserving independent of the mesh topology and retains the consistent behaviour of gradients even on meshes with large curvature and high aspect ratios. The gradient accuracy in MGG reconstruction depends on the accuracy of discretisation of the normal derivatives at faces and this necessitates an iterative approach for gradient computation on non-orthogonal meshes. Numerical studies on different mesh topologies demonstrate that MGG reconstruction gives accurate and consistent gradients on non-orthogonal meshes, with the number of iterations proportional to the extent of non-orthogonality. The MGG reconstruction is found to be consistent even on meshes with large aspect ratio and curvature with the errors being lesser than those from linear least-squares reconstruction. A non-iterative strategy in conjunction with MGG reconstruction is proposed for gradient computations in finite volume simulations that achieves the accuracy and robustness of MGG reconstruction at a cost equivalent to that of SGG reconstruction. The efficacy of this strategy for fluid flow problems is demonstrated through numerical investigations in both incompressible and compressible regimes. The MGG reconstruction may, therefore, be viewed as a novel and promising blend of least-squares and Green–Gauss based approaches which can be implemented with little effort in open-source finite-volume solvers and legacy codes.

The first paper of this series presents a discretely entropy stable discontinuous Galerkin (DG) method for the resistive magnetohydrodynamics (MHD) equations on three-dimensional curvilinear unstructured hexahedral meshes. Compared to other fluid dynamics systems such as the shallow water equations or the compressible Navier–Stokes equations, the resistive MHD equations need special considerations because of the divergence-free constraint on the magnetic field. For instance, it is well known that for the symmetrization of the ideal MHD system as well as the continuous entropy analysis a non-conservative term proportional to the divergence of the magnetic field, typically referred to as the Powell term, must be included. As a consequence, the mimicry of the continuous entropy analysis in the discrete sense demands a suitable DG approximation of the non-conservative terms in addition to the ideal MHD terms. This paper focuses on the resistive MHD equations: Our first contribution is a proof that the resistive terms are symmetric and positive-definite when formulated in entropy space as gradients of the entropy variables, which enables us to show that the entropy inequality holds for the resistive MHD equations. This continuous analysis is the key for our DG discretization and guides the path for the construction of an approximation that discretely mimics the entropy inequality, typically termed entropy stability. Our second contribution is a detailed derivation and analysis of the discretization on three-dimensional curvilinear meshes. The discrete analysis relies on the summation-by-parts property, which is satisfied by the DG spectral element method (DGSEM) with Legendre–Gauss–Lobatto (LGL) nodes. Although the divergence-free constraint is included in the non-conservative terms, the resulting method has no particular treatment of the magnetic field divergence errors, which might pollute the solution quality. Our final contribution is the extension of the standard resistive MHD equations and our DG approximation with a divergence cleaning mechanism that is based on a generalized Lagrange multiplier (GLM). As a conclusion to the first part of this series, we provide detailed numerical validations of our DGSEM method that underline our theoretical derivations. In addition, we show a numerical example where the entropy stable DGSEM demonstrates increased robustness compared to the standard DGSEM.