Solving elasticity equations with interfaces is a challenging problem for most existing methods. Nonetheless, it has wide applications in engineering and science. An accurate and efficient method is desired. In this paper, an efficient non-traditional finite element method with non-body-fitting grids is proposed to solve elasticity equations with interfaces. The main idea is to choose the test function basis to be the standard finite element basis independent of the interface and to choose the solution basis to be piecewise linear satisfying the jump conditions across the interface. The resulting linear system of equations is shown to be positive definite under certain assumptions. Numerical experiments show that this method is second order accurate in the L∞ norm for piecewise smooth solutions. More than 1.5th order accuracy is observed for solution with singularity (second derivative blows up) on the sharp-edged interface corner.

In this article we present a new family of high order accurate Arbitrary Lagrangian-Eulerian one-step WENO finite volume schemes for the solution of stiff hyperbolic balance laws. High order accuracy in space is obtained with a standard WENO reconstruction algorithm and high order in time is obtained using the local space-time discontinuous Galerkin method recently proposed in [20]. In the Lagrangian framework considered here, the local space-time DG predictor is based on a weak formulation of the governing PDE on a moving space-time element. For the space-time basis and test functions we use Lagrange interpolation polynomials defined by tensor-product Gauss-Legendre quadrature points. The moving space-time elements are mapped to a reference element using an isoparametric approach, i.e. the space-time mapping is defined by the same basis functions as the weak solution of the PDE. We show some computational examples in one space-dimension for non-stiff and for stiff balance laws, in particular for the Euler equations of compressible gas dynamics, for the resistive relativistic MHD equations, and for the relativistic radiation hydrodynamics equations. Numerical convergence results are presented for the stiff case up to sixth order of accuracy in space and time and for the non-stiff case up to eighth order of accuracy in space and time.

We discuss the development, verification, and performance of a GPU accelerated discontinuous Galerkin method for the solutions of two dimensional nonlinear shallow water equations. The shallow water equations are hyperbolic partial differential equations and are widely used in the simulation of tsunami wave propagations. Our algorithms are tailored to take advantage of the single instruction multiple data (SIMD) architecture of graphic processing units. The time integration is accelerated by local time stepping based on a multi-rate Adams-Bashforthscheme. A total variational bounded limiter is adopted for nonlinear stability of the numerical scheme. This limiter is coupled with a mass and momentum conserving positivity preserving limiter for the special treatment of a dry or partially wet element in the triangulation. Accuracy, robustness and performance are demonstrated with the aid of test cases. Furthermore, we developed a unified multi-threading model OCCA. The kernels expressed in OCCA model can be cross-compiled with multi-threading models OpenCL, CUDA, and OpenMP. We compare the performance of the OCCA kernels when cross-compiled with these models.

The split-operator technique for wave packet propagation in quantum systems is expanded here to the case of propagatingwave functions describing Schrödinger particles, namely, charge carriers in semiconductor nanostructures within the effective mass approximation, in the presence of Zeeman effect, as well as of Rashba and Dresselhaus spin-orbit interactions. We also demonstrate that simple modifications to the expanded technique allow us to calculate the time evolution of wave packets describing Dirac particles, which are relevant for the study of transport properties in graphene.

In this paper, a diffuse-interface immersed boundary method (IBM) is proposed to treat three different thermal boundary conditions (Dirichlet, Neumann, Robin) in thermal flow problems. The novel IBM is implemented combining with the lattice Boltzmann method (LBM). The present algorithm enforces the three types of thermal boundary conditions at the boundary points. Concretely speaking, the IBM for the Dirichlet boundary condition is implemented using an iterative method, and its main feature is to accurately satisfy the given temperature on the boundary. The Neumann and Robin boundary conditions are implemented in IBM by distributing the jump of the heat flux on the boundary to surrounding Eulerian points, and the jump is obtained by applying the jump interface conditions in the normal and tangential directions. A simple analysis of the computational accuracy of IBM is developed. The analysis indicates that the Taylor-Green vortices problem which was used in many previous studies is not an appropriate accuracy test example. The capacity of the present thermal immersed boundary method is validated using four numerical experiments: (1) Natural convection in a cavity with a circular cylinder in the center; (2) Flows over a heated cylinder; (3) Natural convection in a concentric horizontal cylindrical annulus; (4) Sedimentation of a single isothermal cold particle in a vertical channel. The numerical results show good agreements with the data in the previous literatures.

In this paper we generalize a new type of compact Hermite weighted essentially non-oscillatory (HWENO) limiter for the Runge-Kutta discontinuous Galerkin (RKDG) method, which was recently developed in [38] for structured meshes, to two dimensional unstructured meshes. The main idea of this HWENO limiter is to reconstruct the new polynomial by the usage of the entire polynomials of the DG solution from the target cell and its neighboring cells in a least squares fashion [11] while maintaining the conservative property, then use the classical WENO methodology to form a convex combination of these reconstructed polynomials based on the smoothness indicators and associated nonlinear weights. The main advantage of this new HWENO limiter is the robustness for very strong shocks and simplicity in implementation especially for the unstructured meshes considered in this paper, since only information from the target cell and its immediate neighbors is needed. Numerical results for both scalar and system equations are provided to test and verify the good performance of this new limiter.

Several recently developed AUSM-family numerical flux functions (SLAU, SLAU2, AUSMM+-up2, and AUSMPW+) have been successfully extended to compute compressible multiphase flows, based on the stratified flow model concept, by following two previous works: one by M.-S. Liou, C.-H. Chang, L. Nguyen, and T.G. Theofanous [AIAA J. 46:2345-2356, 2008], in which AUSM+-up was used entirely, and the other by C.-H. Chang, and M.-S. Liou [J. Comput. Phys. 225:840-873, 2007], in which the exact Riemann solver was combined into AUSM+-up at the phase interface. Through an extensive survey by comparing flux functions, the following are found: (1) AUSM+-up with dissipation parameters of Kp and Ku equal to 0.5 or greater, AUSMPW+, SLAU2, AUSM+-up2, and SLAU can be used to solve benchmark problems, including a shock/water-droplet interaction; (2) SLAU shows oscillatory behaviors [though not as catastrophic as those of AUSM+ (a special case of AUSM+-up with Kp = Ku = 0)] due to insufficient dissipation arising from its ideal-gas-based dissipation term; and (3) when combined with the exact Riemann solver, AUSM+-up (Kp = Ku = 1), SLAU2, and AUSMPW+ are applicable to more challenging problems with high pressure ratios.

The conventional Poisson-Nernst-Planck equations do not account for the finite size of ions explicitly. This leads to solutions featuring unrealistically high ionic concentrations in the regions subject to external potentials, in particular, near highly charged surfaces. A modified form of the Poisson-Nernst-Planck equations accounts for steric effects and results in solutions with finite ion concentrations. Here, we evaluate numerical methods for solving the modified Poisson-Nernst-Planck equations by modeling electric field-driven transport of ions through a nanopore. We describe a novel, robust finite element solver that combines the applications of the Newton’s method to the nonlinear Galerkin form of the equations, augmented with stabilization terms to appropriately handle the drift-diffusion processes.

To make direct comparison with particle-based simulations possible, our method is specifically designed to produce solutions under periodic boundary conditions and to conserve the number of ions in the solution domain. We test our finite element solver on a set of challenging numerical experiments that include calculations of the ion distribution in a volume confined between two charged plates, calculations of the ionic current though a nanopore subject to an external electric field, and modeling the effect of a DNA molecule on the ion concentration and nanopore current.

Defects in liquid crystals are of great practical importance and theoretical interest. Despite tremendous efforts, predicting the location and transition of defects under various topological constraint and external field remains to be a challenge. We investigate defect patterns of nematic liquid crystals confined in three-dimensional spherical droplet and two-dimensional disk under different boundary conditions, within the Landau-de Gennes model. We implement a spectral method that numerically solves the Landau-de Gennes model with high accuracy, which allows us to study the detailed static structure of defects. We observe five types of defect structures. Among them the 1/2-disclination lines are the most stable structure at low temperature. Inspired by numerical results, we obtain the profile of disclination lines analytically. Moreover, the connection and difference between defect patterns under the Landau-de Gennes model and the Oseen-Frank model are discussed. Finally, three conjectures are made to summarize some important characteristics of defects in the Landau-de Gennes theory. This work is a continuing effort to deepen our understanding on defect patterns in nematic liquid crystals.

A genuine finite volume method based on the lattice Boltzmann equation (LBE) for nearly incompressible flows is developed. The proposed finite volume lattice Boltzmann method (FV-LBM) is grid-transparent, i.e., it requires no knowledge of cell topology, thus it can be implemented on arbitrary unstructured meshes for effective and efficient treatment of complex geometries. Due to the linear advection term in the LBE, it is easy to construct multi-dimensional schemes. In addition, inviscid and viscous fluxes are computed in one step in the LBE, as opposed to in two separate steps for the traditional finite-volume discretization of the Navier-Stokes equations. Because of its conservation constraints, the collision term of the kinetic equation can be treated implicitly without linearization or any other approximation, thus the computational efficiency is enhanced. The collision with multiple-relaxation-time (MRT) model is used in the LBE. The developed FV-LBM is of second-order convergence. The proposed FV-LBM is validated with three test cases in two-dimensions: (a) the Poiseuille flow driven by a constant body force; (b) the Blasius boundary layer; and (c) the steady flow past a cylinder at the Reynolds numbers Re=10, 20, and 40. The results verify the designed accuracy and efficacy of the proposed FV-LBM.

In this paper, direct numerical simulation (DNS) is presented for spatially evolving turbulent boundary layer over an isothermal flat-plate at Ma∞ = 2.25,5,6,8. When Ma∞ = 8, two cases with the ratio of wall-to-reference temperature Tω/T∞ = 1.9 and 10.03 are considered respectively. The wall temperature approaches recovery temperatures for other cases. The characteristics of compressible turbulent boundary layer (CTBL) affected by freestream Mach number and wall temperature are investigated. It focuses on assessing compressibility effects and the validity of Morkovin's hypothesis through computing and analyzing the mean velocity profile, turbulent intensity, the strong Reynolds analogy (SRA) and possibility density function of dilatation term. The results show that, when the wall temperature approaches recovery temperature, the effects of Mach number on compressibility is insignificant. As a result, the compressibility effect is very weak and the Morkovin's hypothesis is still valid for Mach number even up to 8. However, when Mach number equal to 8, the wall temperature effect on the compressibility is sensitive. In this case, when Tω/T∞ = 1.9, the Morkovin's hypothesis is not fully valid. The validity of classical SRA depends on wall temperature directly. A new modified SRA is proposed to eliminate such negative factor in near wall region. Finally the effects of Mach number and wall temperature on streaks are also studied.

This paper is concerned with numerical solutions of the LDG method for 1D wave equations. Superconvergence and energy conserving properties have been studied. We first study the superconvergence phenomenon for linear problems when alternating fluxes are used. We prove that, under some proper initial discretization, the numerical trace of the LDG approximation at nodes, as well as the cell average, converge with an order 2k+1. In addition, we establish k+2-th order and k+1-th order superconvergence rates for the function value error and the derivative error at Radau points, respectively. As a byproduct, we prove that the LDG solution is superconvergent with an order k+2 towards the Radau projection of the exact solution. Numerical experiments demonstrate that in most cases, our error estimates are optimal, i.e., the error bounds are sharp. In the second part, we propose a fully discrete numerical scheme that conserves the discrete energy. Due to the energy conserving property, after long time integration, our method still stays accurate when applied to nonlinear Klein-Gordon and Sine-Gordon equations.

We describe in this work a discontinuous-Galerkin Finite-Element method to approximate the solutions of a new family of 1d Green-Naghdi models. These new models are shown to be more computationally efficient, while being asymptotically equivalent to the initial formulation with regard to the shallowness parameter. Using the free surface instead of the water height as a conservative variable, the models are recasted under a pre-balanced formulation and discretized using a nodal expansion basis. Independently from the polynomial degree in the approximation space, the preservation of the motionless steady-states is automatically ensured, and the water height positivity is enforced. A simple numerical procedure devoted to stabilize the computations in the vicinity of broken waves is also described. The validity of the resulting model is assessed through extensive numerical validations.

We develop a super-grid modeling technique for solving the elastic wave equation in semi-bounded two- and three-dimensional spatial domains. In this method, waves are slowed down and dissipated in sponge layers near the far-field boundaries. Mathematically, this is equivalent to a coordinate mapping that transforms a very large physical domain to a significantly smaller computational domain, where the elastic wave equation is solved numerically on a regular grid. To damp out waves that become poorly resolved because of the coordinate mapping, a high order artificial dissipation operator is added in layers near the boundaries of the computational domain. We prove by energy estimates that the super-grid modeling leads to a stable numerical method with decreasing energy, which is valid for heterogeneous material properties and a free surface boundary condition on one side of the domain. Our spatial discretization is based on a fourth order accurate finite difference method, which satisfies the principle of summation by parts. We show that the discrete energy estimate holds also when a centered finite difference stencil is combined with homogeneous Dirichlet conditions at several ghost points outside of the far-field boundaries. Therefore, the coefficients in the finite difference stencils need only be boundary modified near the free surface. This allows for improved computational efficiency and significant simplifications of the implementation of the proposed method in multi-dimensional domains. Numerical experiments in three space dimensions show that the modeling error from truncating the domain can be made very small by choosing a sufficiently wide super-grid damping layer. The numerical accuracy is first evaluated against analytical solutions of Lamb’s problem, where fourth order accuracy is observed with a sixth order artificial dissipation. We then use successive grid refinements to study the numerical accuracy in the more complicated motion due to a point moment tensor source in a regularized layered material.

In this paper, we introduce a high-order accurate constrained transport type finite volume method to solve ideal magnetohydrodynamic equations on two-dimensional triangular meshes. A new divergence-free WENO-based reconstruction method is developed to maintain exactly divergence-free evolution of the numerical magnetic field. In this formulation, the normal component of the magnetic field at each face of a triangle is reconstructed uniquely and with the desired order of accuracy. Additionally, a new weighted flux interpolation approach is also developed to compute the z-component of the electric field at vertices of grid cells. We also present numerical examples to demonstrate the accuracy and robustness of the proposed scheme.

This article is to discuss the bilinear and linear immersed finite element (IFE) solutions generated from the algebraic multigrid solver for both stationary and moving interface problems. For the numerical methods based on finite difference formulation and a structured mesh independent of the interface, the stiffness matrix of the linear system is usually not symmetric positive-definite, which demands extra efforts to design efficient multigrid methods. On the other hand, the stiffness matrix arising from the IFE methods are naturally symmetric positive-definite. Hence the IFE-AMG algorithm is proposed to solve the linear systems of the bilinear and linear IFE methods for both stationary and moving interface problems. The numerical examples demonstrate the features of the proposed algorithms, including the optimal convergence in both L2 and semi-H1 norms of the IFE-AMG solutions, the high efficiency with proper choice of the components and parameters of AMG, the influence of the tolerance and the smoother type of AMG on the convergence of the IFE solutions for the interface problems, and the relationship between the cost and the moving interface location.

The most severe limitation of the standard Lattice Boltzmann Method is the use of uniform Cartesian grids especially when there is a need for high resolutions near the body or the walls. Among the recent advances in lattice Boltzmann research to handle complex geometries, a particularly remarkable option is represented by changing the solution procedure from the original “stream and collide” to a finite volume technique. However, most of the presented schemes have stability problems. This paper presents a stable and accurate finite-volume lattice Boltzmann formulation based on a cell-centred scheme. To enhance stability, upwind second order pressure biasing factors are used as flux correctors on a D2Q9 lattice. The resulting model has been tested against a uniform flow past a cylinder and typical free shear flow problems at low and moderate Reynolds numbers: boundary layer, mixing layer and plane jet flows. The numerical results show a very good accuracy and agreement with the exact solution of the Navier-Stokes equation and previous numerical results and/or experimental data. Results in self-similar coordinates are also investigated and show that the time-averaged statistics for velocity and vorticity express self-similarity at low Reynolds numbers. Furthermore, the scheme is applied to simulate the flow around circular cylinder and the Reynolds number range is chosen in such a way that the flow is time dependent. The agreement of the numerical results with previous results is satisfactory.

The adaptive generalized Riemann problem (GRP) scheme for 2-D compressible fluid flows has been proposed in [J. Comput. Phys., 229 (2010), 1448-1466] and it displays the capability in overcoming difficulties such as the start-up error for a single shock, and the numerical instability of the almost stationary shock. In this paper, we will provide the accuracy study and particularly show the performance in simulating 2-D complex wave configurations formulated with the 2-D Riemann problems for compressible Euler equations. For this purpose, we will first review the GRP scheme briefly when combined with the adaptive moving mesh technique and consider the accuracy of the adaptive GRP scheme via the comparison with the explicit formulae of analytic solutions of planar rarefaction waves, planar shock waves, the collapse problem of a wedge-shaped dam and the spiral formation problem. Then we simulate the full set of wave configurations in the 2-D four-wave Riemann problems for compressible Euler equations [SIAM J. Math. Anal., 21 (1990), 593-630], including the interactions of strong shocks (shock reflections), vortex-vortex and shock-vortex etc. This study combines the theoretical results with the numerical simulations, and thus demonstrates what Ami Harten observed “for computational scientists there are two kinds of truth: the truth that you prove, and the truth you see when you compute” [J. Sci. Comput., 31 (2007), 185-193].

A mass-conserving lattice Boltzmann model based on the Bhatnagar-Gross-Krook (BGK) model is proposed for non-Newtonian fluid flows. The equilibrium distribution function includes the local shear rate related with the viscosity and a variable parameter changing with the shear rate. With the additional parameter, the relaxation time in the collision can be fixed invariable to the viscosity. Through the Chapman-Enskog analysis, the macroscopic equations can be recovered from the present mass-conserving model. Two flow problems are simulated to validate the present model with a local computing scheme for the shear rate, and good agreement with analytical solutions and/or other published results are obtained. The results also indicate that the present modified model is more applicable to practical non-Newtonian fluid flows owing to its better accuracy and more robustness than previous methods.

In this paper, we numerically study the ground and first excited states of the fractional Schrödinger equation in an infinite potential well. Due to the nonlocality of the fractional Laplacian, it is challenging to find the eigenvalues and eigenfunctions of the fractional Schrödinger equation analytically. We first introduce a normalized fractional gradient flow and then discretize it by a quadrature rule method in space and the semi-implicit Euler method in time. Our numerical results suggest that the eigenfunctions of the fractional Schrödinger equation in an infinite potential well differ from those of the standard (non-fractional) Schrödinger equation. We find that the strong nonlocal interactions represented by the fractional Laplacian can lead to a large scattering of particles inside of the potential well. Compared to the ground states, the scattering of particles in the first excited states is larger. Furthermore, boundary layers emerge in the ground states and additionally inner layers exist in the first excited states of the fractional nonlinear Schrödinger equation. Our simulated eigenvalues are consistent with the lower and upper bound estimates in the literature.