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

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

Continuous and discontinuous Galerkin methods for a scalable three-dimensional nonhydrostatic atmospheric model: Limited-area mode

This paper describes a unified, element based Galerkin (EBG) framework for a three-dimensional, nonhydrostatic model for the atmosphere. In general, EBG methods possess high-order accuracy, geometric flexibility, excellent dispersion properties and good scalability. Our nonhydrostatic model, based on the compressible Euler equations, is appropriate for both limited-area and global atmospheric simulations. Both a continuous Galerkin (CG), or spectral element, and discontinuous Galerkin (DG) model are considered using hexahedral elements. The formulation is suitable for both global and limited-area atmospheric modeling, although we restrict our attention to 3D limited-area phenomena in this study; global atmospheric simulations will be presented in a follow-up paper. Domain decomposition and communication algorithms used by both our CG and DG models are presented. The communication volume and exchange algorithms for CG and DG are compared and contrasted. Numerical verification of the model was performed using two test cases: flow past a 3D mountain and buoyant convection of a bubble in a neutral atmosphere; these tests indicate that both CG and DG can simulate the necessary physics of dry atmospheric dynamics. Scalability of both methods is shown up to 8192 CPU cores, with near ideal scaling for DG up to 32,768 cores.     

An interface treating technique for compressible multi-medium flow with Runge–Kutta discontinuous Galerkin method

The high-order accurate Runge–Kutta discontinuous Galerkin (RKDG) method is applied to the simulation of compressible multi-medium flow, generalizing the interface treating method given in Chertock et al. (2008) [9]. In mixed cells, where the interface is located, Riemann problems are solved to define the states on both sides of the interface. The input states to the Riemann problem are obtained by extrapolation to the cell boundary from solution polynomials in the neighbors of the mixed cell. The level set equation is solved by using a high-order accurate RKDG method for Hamilton–Jacobi equations, resulting in a unified DG solver for the coupled problem. The method is conservative if we include the states in the mixed cells, which are however not used in the updating of the numerical solution in other cells. The states in the mixed cells are plotted to better evaluate the conservation errors, manifested by overshoots/undershoots when compared with states in neighboring cells. These overshoots/undershoots in mixed cells are problem dependent and change with time. Numerical examples show that the results of our scheme compare well with other methods for one and two-dimensional problems. In particular, the algorithm can capture well complex flow features of the one-dimensional shock entropy wave interaction problem and two-dimensional shock–bubble interaction problem.     

Finite volume TVD formulation of lattice Boltzmann simulation on unstructured mesh

A numerical scheme is presented for accurate simulation of fluid flow using the lattice Boltzmann equation (LBE) on unstructured mesh. A finite volume approach is adopted to discretize the LBE on a cell-centered, arbitrary shaped, triangular tessellation. The formulation includes a formal, second order discretization using a Total Variation Diminishing (TVD) scheme for the terms representing advection of the distribution function in physical space, due to microscopic particle motion. The advantage of the LBE approach is exploited by implementing the scheme in a new computer code to run on a parallel computing system. Performance of the new formulation is systematically investigated by simulating four benchmark flows of increasing complexity, namely (1) flow in a plane channel, (2) unsteady Couette flow, (3) flow caused by a moving lid over a 2D square cavity and (4) flow over a circular cylinder. For each of these flows, the present scheme is validated with the results from Navier–Stokes computations as well as lattice Boltzmann simulations on regular mesh. It is shown that the scheme is robust and accurate for the different test problems studied.     

Positivity preserving semi-Lagrangian discontinuous Galerkin formulation: Theoretical analysis and application to the Vlasov–Poisson system

Semi-Lagrangian (SL) methods have been very popular in the Vlasov simulation community , , , , , ,  and . In this paper, we propose a new Strang split SL discontinuous Galerkin (DG) method for solving the Vlasov equation. Specifically, we apply the Strang splitting for the Vlasov equation [6], as a way to decouple the nonlinear Vlasov system into a sequence of 1-D advection equations, each of which has an advection velocity that only depends on coordinates that are transverse to the direction of propagation. To evolve the decoupled linear equations, we propose to couple the SL framework with the semi-discrete DG formulation. The proposed SL DG method is free of time step restriction compared with the Runge–Kutta DG method, which is known to suffer from numerical time step limitation with relatively small CFL numbers according to linear stability analysis. We apply the recently developed positivity preserving (PP) limiter [37], which is a low-cost black box procedure, to our scheme to ensure the positivity of the unknown probability density function without affecting the high order accuracy of the base SL DG scheme. We analyze the stability and accuracy properties of the SL DG scheme by establishing its connection with the direct and weak formulations of the characteristics/Lagrangian Galerkin method [23]. The quality of the proposed method is demonstrated via basic test problems, such as linear advection and rigid body rotation, and via classical plasma problems, such as Landau damping and the two stream instability.     

A robust SN-DG-approximation for radiation transport in optically thick and diffusive regimes

We introduce a new discontinuous Galerkin (DG) method with reduced upwind stabilization for the linear Boltzmann equation applied to particle transport. The asymptotic analysis demonstrates that the new formulation does not suffer from the limitations of standard upwind methods in the thick diffusive regime; in particular, the new method yields the correct diffusion limit for any approximation order, including piecewise constant discontinuous finite elements. Numerical tests on well-established benchmark problems demonstrate the superiority of the new method. The improvement is particularly significant when employing piecewise constant DG approximation for which standard upwinding is known to perform poorly in the thick diffusion limit.     

An accurate and efficient method for treating aerodynamic interactions of cloud droplets

Motivated by a need to improve the representation of short-range interaction forces in hybrid direct numerical simulation of interacting cloud droplets, an efficient method for treating the aerodynamic interaction of two spherical particles settling under gravity is developed. An effort is made to ensure the accuracy of our method for any inter-particle separation by considering three separation ranges. The first is the long-range interaction where a multipole method is applied. After a decomposition into six simple configurations, explicit formulae for drag forces and torques are derived from an approximate Force–Torque–Stresslet (FTS) formulation. The FTS formulation is found to be accurate when the separation distance normalized by the average radius is larger than 5. The second range concerns the short-range interaction where the interaction force could be very large. Leading-order lubrication expansions are employed for this range and are found to be accurate when the normalized separation is less than about 0.01. Finally, for the intermediate range where no simple method is available, a third-order polynomial fitting is proposed to bridge the treatments for long-range and short-range interactions. After optimizing the precise form of polynomial fitting and matching locations, the force representation is found to be highly accurate when compared with the exact solution for Stokes flows. Using this method, collision efficiencies of cloud droplets sedimenting under gravity have been calculated. It is shown that the results of collision efficiency are in excellent agreement with results based on the exact Stokes flow solution. Collision efficiency results are also compared to previous results to further illustrate the accuracy of our calculations. The effects of particle rotation and the attractive van der Waals force on the collision efficiency are also studied. The efficient force representation developed here is more general than the usual lubrication expansion and thus can serve as a better approach to correct unresolved short-range interactions in particle-resolved simulations.     

Finite element form of FDV for widely varying flowfields

We present the Flowfield Dependent Variation (FDV) method for physical applications that have widely varying spatial and temporal scales. Our motivation is to develop a versatile numerical method that is accurate and stable in simulations with complex geometries and with wide variations in space and time scales. The use of a finite element formulation adds capabilities such as flexible grid geometries and exact enforcement of Neumann boundary conditions. While finite element schemes are used extensively by researchers solving computational fluid dynamics in many engineering fields, their use in space physics, astrophysical fluids and laboratory magnetohydrodynamic simulations with shocks has been predominantly overlooked. The FDV method is unique in that numerical diffusion is derived from physical parameters rather than traditional artificial viscosity methods. Numerical instabilities account for most of the difficulties when capturing shocks in these regimes. The first part of this paper concentrates on the presentation of our numerical method formulation for Newtonian and relativistic hydrodynamics. In the second part we present several standard simulation examples that test the method’s limitations and verify the FDV method. We show that our finite element formulation is stable and accurate for a range of both Mach numbers and Lorentz factors in one-dimensional test problems. We also present the converging/diverging nozzle which contains both incompressible and compressible flow in the flowfield over a range of subsonic and supersonic regions. We demonstrate the stability of our method and the accuracy by comparison with the results of other methods including the finite difference Total Variation Diminishing method. We explore the use of FDV for both non-relativistic and relativistic fluids (hydrodynamics) with strong shocks in order to establish the effectiveness in future applications of this method in astrophysical and laboratory plasma environments.     

Calculation of Conversion Efficiency of a Doubly Resonant Optical Parametric Oscillator

The problem of calculating the conversion efficiency for a doubly resonant optical parametric oscillator with a Fabry-Perot type resonator configuration is reinvestigated, extending the usual theoretical treatment, which is restricted to highly reflecting resonator mirrors, to the case of arbitrary reflectivity R of the output mirror. For decreasing R both the signal and the idler wave acquire growing portions of travelling waves, and it is shown by numerical analysis that this effect leads to a remarkable enhancement of the conversion efficiency (as a function of the relative excitation). In particular, the maximum efficiency may considerably exceed the value 50% to be attained by means of high quality resonators (1 — R < 1).     

Efficient schemes to compute diffusive barrier crossing rates

The formulation of the classical barrier-crossing problem is reviewed in the context of numerical simulations, with the focus on barrier crossing problems where the reaction coordinate depends in a non-trivial way on the Cartesian coordinates of many particles. Often it is convenient to measure the barrier height using constrained dynamics. Such a calculation requires a knowledge of the Jacobian for the coordinate transformation between Cartesian and generalized (‘reaction’) coordinates, and it is shown that the calculation of this Jacobian can be simplified. The conventional expression for the crossing rate is found to become computationally inefficient when the barrier crossing is diffusive. An alternative formulation of the barrier-crossing rate is given that leads to much better statistical accuracy in the computed crossing rates.     

The calculation of electrical parameters of AT-cut quartz crystal resonators with the consideration of material viscosity

Electrical parameters like resistance and quality factor of a quartz crystal resonator cannot be determined through vibration analysis without considering the presence of material dissipation. In this study, we use the first-order Mindlin plate equations of piezoelectric plates for thickness-shear vibrations of a simple resonator model with partial electrodes. We derive the expressions of electrical parameters with emphasis on the resistance that is related to the imaginary part of complex elastic constants, or the viscosity, of quartz crystal. Since all electrical parameters are frequency dependent, this procedure provides the chance to study the frequency behavior of crystal resonators with a direct formulation. We understand that the electrical parameters are strongly affected by the manufacturing process, with the plating techniques in particular, but the theoretical approach we presented here will be the first step for the precise estimation of such parameters and their further applications in the analysis of nonlinear behavior of resonators. We calculated the parameters from our simple resonator model of AT-cut quartz crystal with the first-order Mindlin plate theory to demonstrate the procedure and show that the numerical results are consistent with earlier measurements.     

A High-Order Discontinuous Galerkin Method for the Two-Dimensional Time-Domain Maxwell's Equations on Curved Mesh

In this paper, a DG (Discontinuous Galerkin) method which has been widely employed in CFD (Computational Fluid Dynamics) is used to solve the two-dimensional time-domain Maxwell's equations for complex geometries on unstructured mesh. The element interfaces on solid boundary are treated in both curved way and straight way. Numerical tests are performed for both benchmark problems and complex cases with varying orders on a series of grids, where the high-order convergence in accuracy can be observed. Both the curved and the straight solid boundary implementation can give accurate RCS (Radar Cross-Section) results with sufficiently small mesh size, but the curved solid boundary implementation can significantly improve the accuracy when using relatively large mesh size. More importantly, this CFD-based high-order DG method for the Maxwell's equations is very suitable for complex geometries.     

Efficient prediction of the quality factors of micromechanical resonators

A high quality factor (Q-factor) is one of the major requirements of high-performance resonators. An understanding of the dissipation mechanism is crucial for maximizing the quality factor by reducing the energy loss. Thermoelastic damping has been well-known as the important intrinsic dissipation that affects the quality factor of micro-resonators. In this study, a finite element formulation based on the weak form of fully coupled thermoelastic problems is suggested. The coupled thermoelastic equation usually leads to a large-size complex eigenvalue problem, which is very massive and time-consuming to solve. Therefore, we also applied the model order reduction (MOR) scheme to this coupled multiphysical problem in order to achieve computational efficiency. The present approach is validated by comparing the numerical results and analytical solutions.     

The validation of the parallel three-dimensional solver for analysis of optical plasmonic bi-periodic multilayer nanostructures

Fundamentals of the three-dimensional spatial harmonic analysis (SHA) approach are reviewed, and the advantages of a fast-converging formulation versus the initial SHA formulation are emphasized with examples using periodic plasmonic nanostructures. First, two independent parallel versions of both formulations are implemented using the scattering matrix algorithm for multilayer cascading. Then, by comparing the results from both formulations, it is shown that choosing an advanced fast-converging scheme could be essential for accurate and efficient modeling of plasmonic structures. Important obstacles to the fast parallel implementation of this approach are also revealed. The results of test simulations are validated using the data obtained from a commercial finite-element method (FEM) simulations and from the experimental characterization of fabricated samples.     

Signal displacement in spiral-in acquisitions: simulations and implications for imaging in SFG regions

Susceptibility field gradients (SFGs) cause problems for functional magnetic resonance imaging (fMRI) in regions like the orbital frontal lobes, leading to signal loss and image artifacts (signal displacement and "pile-up"). Pulse sequences with spiral-in k-space trajectories are often used when acquiring fMRI in SFG regions such as inferior/medial temporal cortex because it is believed that they have improved signal recovery and decreased signal displacement properties. Previously postulated theories explain differing reasons why spiral-in appears to perform better than spiral-out; however it is clear that multiple mechanisms are occurring in parallel. This study explores differences in spiral-in and spiral-out images using human and phantom empirical data, as well as simulations consistent with the phantom model. Using image simulations, the displacement of signal was characterized using point spread functions (PSFs) and target maps, the latter of which are conceptually inverse PSFs describing which spatial locations contribute signal to a particular voxel. The magnitude of both PSFs and target maps was found to be identical for spiral-out and spiral-in acquisitions, with signal in target maps being displaced from distant regions in both cases. However, differences in the phase of the signal displacement patterns that consequently lead to changes in the intervoxel phase coherence were found to be a significant mechanism explaining differences between the spiral sequences. The results demonstrate that spiral-in trajectories do preserve more total signal in SFG regions than spiral-out; however, spiral-in does not in fact exhibit decreased signal displacement. Given that this signal can be displaced by significant distances, its recovery may not be preferable for all fMRI applications.     

Bounds on natural frequencies of composite circular membranes: Integral equation methods

This paper is concerned with finding upper and lower bounds for the natural frequencies of vibration of a circular membrane with stepped radial density. Such problems, involving discontinuous coefficients in the differential equation, may be treated by using classical variational methods. However, it is shown here that eigenvalue estimation techniques based on an integral equation formulation are more effective. Integral equation methods provide more accurate upper bounds, for a comparable amount of effort, and supply improvable lower bounds as well.     

Acoustic inversions for measuring boundary layer suspended sediment processes

Although sound has been applied to the study of sediment transport processes for a number of years, it is acknowledged that there are still problems in using the backscattered signal to measure suspended sediment parameters. In particular, when the attenuation due to the suspension becomes significant, the uncertainty associated with the variability in the scattering characteristics of the sediments in suspension can lead to inversion errors which accumulate as the sound propagates through the suspension. To study this attenuation propagation problem, numerical simulations and laboratory experiments have been used to assess the impact unpredictability in the scattering properties of the suspension has on the acoustically derived suspended sediments parameters. The results clearly show the commonly applied iterative implicit inversion can lead to calculated sediment parameters, which become increasingly erroneous with range, as the sound propagates through the suspension. To address this problem an alternative approach to the iterative implicit formulation is investigated using a recently described dual frequency inversion. This approach is not subject to the accumulation of errors and has an explicit solution. Here the dual frequency inversion is assessed and calculated suspended sediment parameters are compared with those obtained from the iterative implicit inversion.     

Planar waveguide-coupled, high-index-contrast, high-Q resonators in chalcogenide glass for sensing

High-index-contrast compact microdisk resonators in thermally evaporated As2S3 and Ge17Sb12S71 chalcogenide glass films are designed and fabricated using standard UV lithography and characterized. Our pulley coupler configuration demonstrates coupling of the resonators to monolithically integrated photonic wire waveguides without resorting to demanding fine-line lithography. Microdisk resonators in As2S3 support whispering-gallery-mode with cavity quality factors (Q) exceeding 2 x 10(5), the highest Q value reported in resonator structures in chalcogenide glasses to the best of our knowledge. We have successfully demonstrated a lab-on-a-chip prototype sensor device with the integration of our resonator with planar microfluidic systems. The sensor shows a refractive index sensitivity of 182 nm/RIU (refractive index unit) and a wavelength resolution of 0.1 pm through a resonant peak fit. This corresponds to a refractive index detection limit of 8 x 10(-7) RIU at 1550 nm in wavelength, which could be further improved by shifting the operating wavelength to a region where water absorption is reduced.     

Optical fiber ring resonator characterization by optical time-domain reflectometry

A novel method for characterization of optical fiber resonators by an optical time-domain reflectometry (OTDR) technique is reported. This easy-to-use technique yields accurate results for cavity lengths ranging from a few meters to several kilometers. A simple relationship is established between the round-trip cavity loss and the position where the OTDR signal is maximum. The value obtained for the round-trip cavity loss turns out to be quite insensitive to uncertainties in the determination of the OTDR maximum position.     

A new multi-layer discrete ordinate approach to radiative transfer in vertically inhomogeneous atmospheres

A recently developed matrix formulation of the discrete ordinate method is extended for application to an inhomogeneous atmosphere. The solution yields fluxes, as well as the complete azimuthal dependence of the intensity at any level in the atmosphere. The numerical aspects of the solution are discussed and numerical verification is provided by comparing computed results with those obtained by other methods. In particular, it is shown that a simple scaling scheme, which removes the positive exponentials in the coefficient matrix when solving for the constants of integration, provides unconditionally stable solutions for arbitrary optical thicknesses. An assessment of the accuracy to be expected is also provided, and it is shown that low-order discrete ordinate approximations yield very accurate flux values.