A conservative modification to the ghost fluid method (GFM) is developed for compressible multiphase flows. The motivation is to eliminate or reduce the conservation error of the GFM without affecting its performance. We track the conservative variables near the material interface and use this information to modify the numerical solution for an interfacing cell when the interface has passed the cell. The modification procedure can be used on the GFM with any base schemes. In this paper we use the fifth order finite difference WENO scheme for the spatial discretization and the third order TVD Runge-Kutta method for the time discretization. The level set method is used to capture the interface. Numerical experiments show that the method is at least mass and momentum conservative and is in general comparable in numerical resolution with the original GFM.

Efficient algorithms for the simulation of strained heteroepitaxial growth with intermixing in 2+1 dimensions are presented. The first of these algorithms is an extension of the energy localization method [T. P. Schulze and P. Smereka, An energy localization principle and its application to fast kinetic Monte Carlo simulation of heteroepitaxial growth, J. Mech. Phys. Sol., 3 (2009), 521-538] from 1+1 to 2+1 dimensions. Two approximations of this basic algorithm are then introduced, one of which treats adatoms in a more efficient manner, while the other makes use of an approximation of the change in elastic energy in terms of local elastic energy density. In both cases, it is demonstrated that a reasonable level of fidelity is achieved. Results are presented showing how the film morphology is affected by misfit and deposition rate. In addition, simulations of stacked quantum dots are also presented.

An all speed scheme for the Isentropic Euler equations is presented in this paper. When the Mach number tends to zero, the compressible Euler equations converge to their incompressible counterpart, in which the density becomes a constant. Increasing approximation errors and severe stability constraints are the main difficulty in the low Mach regime. The key idea of our all speed scheme is the special semi-implicit time discretization, in which the low Mach number stiff term is divided into two parts, one being treated explicitly and the other one implicitly. Moreover, the flux of the density equation is also treated implicitly and an elliptic type equation is derived to obtain the density. In this way, the correct limit can be captured without requesting the mesh size and time step to be smaller than the Mach number. Compared with previous semi-implicit methods, firstly, nonphysical oscillations can be suppressed by choosing proper parameter, besides, only a linear elliptic equation needs to be solved implicitly which reduces much computational cost. We develop this semi-implicit time discretization in the framework of a first order Local Lax-Friedrichs (or Rusanov) scheme and numerical tests are displayed to demonstrate its performances.

The important task of evaluating the impact of random parameters on the output of stochastic ordinary differential equations (SODE) can be computationally very demanding, in particular for problems with a high-dimensional parameter space. In this work we consider this problem in some detail and demonstrate that by combining several techniques one can dramatically reduce the overall cost without impacting the predictive accuracy of the output of interests. We discuss how the combination of ANOVA expansions, different sparse grid techniques, and the total sensitivity index (TSI) as a pre-selective mechanism enables the modeling of problems with hundred of parameters. We demonstrate the accuracy and efficiency of this approach on a number of challenging test cases drawn from engineering and science.

This article describes a number of velocity-based moving mesh numerical methods for multidimensional nonlinear time-dependent partial differential equations (PDEs). It consists of a short historical review followed by a detailed description of a recently developed multidimensional moving mesh finite element method based on conservation. Finite element algorithms are derived for both mass-conserving and non mass-conserving problems, and results shown for a number of multidimensional nonlinear test problems, including the second order porous medium equation and the fourth order thin film equation as well as a two-phase problem. Further applications and extensions are referenced.

In multi physics computations where a compressible fluid is coupled with a linearly elastic solid, it is standard to enforce continuity of the normal velocities and of the normal stresses at the interface between the fluid and the solid. In a numerical scheme, there are many ways that velocity- and stress-continuity can be enforced in the discrete approximation. This paper performs a normal mode stability analysis of the linearized problem to investigate the stability of different numerical interface conditions for a model problem approximated by upwind type finite difference schemes. The analysis shows that depending on the ratio of densities between the solid and the fluid, some numerical interface conditions are stable up to the maximal CFL-limit, while other numerical interface conditions suffer from a severe reduction of the stable CFL-limit. The paper also presents a new interface condition, obtained as a simplified characteristic boundary condition, that is proved to not suffer from any reduction of the stable CFL-limit. Numerical experiments in one space dimension show that the new interface condition is stable also for computations with the non-linear Euler equations of compressible fluid flow coupled with a linearly elastic solid.

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 new criterion is presented to detect global convergence to steady state, and to identify local transient characteristics, during rarefied gas flow simulations performed using the direct simulation Monte Carlo (DSMC) method. Unlike deterministic computational fluid dynamics (CFD) schemes, DSMC is generally subject to large statistical scatter in instantaneous flow property evaluations, which prevents the use of residual tracking procedures as are often employed in CFD simulations. However, reliable prediction of the time to reach steady state is necessary for initialization of DSMC sampling operations. Techniques currently used in DSMC to identify steady state convergence are usually insensitive to weak transient behavior in small regions of relatively low density or recirculating flow. The proposed convergence criterion is developed with the goal of properly identifying such weak transient behavior, while adding negligible computational expense and allowing simple implementation in any existing DSMC code. Benefits of the proposed technique over existing convergence detection methods are demonstrated for representative nozzle/plume expansion flow, hypersonic blunt body flow and driven cavity flow problems.

A new finite element level set method is developed to simulate the interface motion. The normal velocity of the moving interface can depend on both the local geometry, such as the curvature, and the external force such as that due to the flux from both sides of the interface of a material whose concentration is governed by a diffusion equation. The key idea of the method is to use an interface-fitted finite element mesh. Such an approximation of the interface allows an accurate calculation of the solution to the diffusion equation. The interface-fitted mesh is constructed from a base mesh, a uniform finite element mesh, at each time step to explicitly locate the interface and separate regions defined by the interface. Several new level set techniques are developed in the framework of finite element methods. These include a simple finite element method for approximating the curvature, a new method for the extension of normal velocity, and a finite element least-squares method for the reinitialization of level set functions. Application of the method to the classical solidification problem captures the dendrites. The method is also applied to the molecular solvation to determine optimal solute-solvent interfaces of solvation systems.

The equilibrium metric for minimizing a continuous congested traffic model is the solution of a variational problem involving geodesic distances. The continuous equilibrium metric and its associated variational problem are closely related to the classical discrete Wardrop’s equilibrium. We propose an adjoint state method to numerically approximate continuous traffic congestion equilibria through the continuous formulation. The method formally derives an adjoint state equation to compute the gradient descent direction so as to minimize a nonlinear functional involving the equilibrium metric and the resulting geodesic distances. The geodesic distance needed for the state equation is computed by solving a factored eikonal equation, and the adjoint state equation is solved by a fast sweeping method. Numerical examples demonstrate that the proposed adjoint state method produces desired equilibrium metrics and outperforms the subgradient marching method for computing such equilibrium metrics.

We propose a deterministic solver for the time-dependent multi-subband Boltzmann transport equation (MSBTE) for the two dimensional (2D) electron gas in double gate metal oxide semiconductor field effect transistors (MOSFETs) with flared out source/drain contacts. A realistic model with six-valleys of the conduction band of silicon and both intra-valley and inter-valley phonon-electron scattering is solved. We propose a second order finite volume method based on the positive and flux conservative (PFC) method to discretize the Boltzmann transport equations (BTEs). The transport part of the BTEs is split into two problems. One is a 1D transport problem in the position space, and the other is a 2D transport problem in the wavevector space. In order to reduce the splitting error, the 2D transport problem in the wavevector space is solved directly by using the PFC method instead of splitting into two 1D problems. The solver is applied to a nanoscale double gate MOSFET and the current-voltage characteristic is investigated. Comparison of the numerical results with ballistic solutions show that the scattering influence is not ignorable even when the size of a nanoscale semiconductor device goes to the scale of the electron mean free path.

A data-driven model reduction strategy is presented for the representation of random polycrystal microstructures. Given a set of microstructure snapshots that satisfy certain statistical constraints such as given low-order moments of the grain size distribution, using a non-linear manifold learning approach, we identify the intrinsic low-dimensionality of the microstructure manifold. In addition to grain size, a linear dimensionality reduction technique (Karhunun-Loéve Expansion) is used to reduce the texture representation. The space of viable microstructures is mapped to a low-dimensional region thus facilitating the analysis and design of polycrystal microstructures. This methodology allows us to sample microstructure features in the reduced-order space thus making it a highly efficient, low-dimensional surrogate for representing microstructures (grain size and texture). We demonstrate the model reduction approach by computing the variability of homogenized thermal properties using sparse grid collocation in the reduced-order space that describes the grain size and orientation variability.

A computational study on the enhancement of the second harmonic generation (SHG) in one-dimensional (1D) photonic crystals is presented. The mathematical model is derived from a nonlinear system of Maxwell’s equations, which partly overcomes the shortcoming of some existing models based on the undepleted pump approximation. We designed an iterative scheme coupled with the finite element method which can be applied to simulate the SHG in one dimensional nonlinear photonic band gap structures in our previous work. For the case that the nonlinearity is strong which is desirable to enhance the conversion efficiency, a continuation method is introduced to ensure the convergence of the iterative procedure. The convergence of our method is fast. Numerical experiments also indicate the conversion efficiency of SHG can be significantly enhanced when the frequencies of the fundamental and the second harmonic wave are tuned at the photonic band edges. The maximum total conversion efficiency available reaches more than 50% in all the cases studied.

We propose an a-posteriori error/smoothness indicator for standard semi-discrete finite volume schemes for systems of conservation laws, based on the numerical production of entropy. This idea extends previous work by the first author limited to central finite volume schemes on staggered grids. We prove that the indicator converges to zero with the same rate of the error of the underlying numerical scheme on smooth flows under grid refinement. We construct and test an adaptive scheme for systems of equations in which the mesh is driven by the entropy indicator. The adaptive scheme uses a single nonuniform grid with a variable timestep. We show how to implement a second order scheme on such a space-time non uniform grid, preserving accuracy and conservation properties. We also give an example of a p-adaptive strategy.

This paper presents an extension of the lumped particle model in [1] to include the effects of particle collisions. The lumped particle model is a flexible framework for the modeling of particle laden flows, that takes into account fundamental features, including advection, diffusion and dispersion of the particles. In this paper, we transform a binary collision model and concepts from kinetic theory into a collision procedure for the lumped particle framework. We apply this new collision procedure to investigate numerically the role of particle collisions in the hindered settling effect. The hindered settling effect is characterized by an increase in the effective drag coefficient CD that influences each particle in the flow. This coefficient is given by , where ϕ is the volume fraction of particles, is the drag coefficient for a single particle, and n ≃ 4.67 for creeping flow. We obtain an approximation for CD/CD by calculating the effective work done by collisions, and comparing that to the work done by the drag force. In our numerical experiments, we observe a minimal value of n = 3.0. Moreover, by allowing high energy dissipation, an approximation for the classical value for creeping flow, n = 4.7, is reproduced. We also obtain high values for n, up to n = 6.5, which is consistent with recent physical experiments on the sedimentation of sand grains.

We propose a multigrid method to solve the molecular mechanics model (molecular dynamics at zero temperature). The Cauchy-Born elasticity model is employed as the coarse grid operator and the elastically deformed state as the initial guess of the molecular mechanics model. The efficiency of the algorithm is demonstrated by three examples with homogeneous deformation, namely, one dimensional chain under tensile deformation and aluminum under tension and shear deformations. The method exhibits linear-scaling computational complexity, and is insensitive to parameters arising from iterative solvers. In addition, we study two examples with inhomogeneous deformation: vacancy and nanoindentation of aluminum. The results are still satisfactory while the linear-scaling property is lost for the latter example.

This paper is concerned with a new version of the Osher-Solomon Riemann solver and is based on a numerical integration of the path-dependent dissipation matrix. The resulting scheme is much simpler than the original one and is applicable to general hyperbolic conservation laws, while retaining the attractive features of the original solver: the method is entropy-satisfying, differentiable and complete in the sense that it attributes a different numerical viscosity to each characteristic field, in particular to the intermediate ones, since the full eigenstructure of the underlying hyperbolic system is used. To illustrate the potential of the proposed scheme we show applications to the following hyperbolic conservation laws: Euler equations of compressible gas-dynamics with ideal gas and real gas equation of state, classical and relativistic MHD equations as well as the equations of nonlinear elasticity. To the knowledge of the authors, apart from the Euler equations with ideal gas, an Osher-type scheme has never been devised before for any of these complicated PDE systems. Since our new general Riemann solver can be directly used as a building block of high order finite volume and discontinuous Galerkin schemes we also show the extension to higher order of accuracy and multiple space dimensions in the new framework of PNPM schemes on unstructured meshes recently proposed in [9].

For a new nonlinear iterative method named as Picard-Newton (P-N) iterative method for the solution of the time-dependent reaction-diffusion systems, which arise in non-equilibrium radiation diffusion applications, two time step control methods are investigated and a study of temporal accuracy of a first order time integration is presented. The non-equilibrium radiation diffusion problems with flux limiter are considered, which appends pesky complexity and nonlinearity to the diffusion coefficient. Numerical results are presented to demonstrate that compared with Picard method, for a desired accuracy, significant increase in solution efficiency can be obtained by Picard-Newton method with the suitable time step size selection.

We present a novel adaptive finite element method (AFEM) for elliptic equations which is based upon the Centroidal Voronoi Tessellation (CVT) and superconvergent gradient recovery. The constructions of CVT and its dual Centroidal Voronoi Delaunay Triangulation (CVDT) are facilitated by a localized Lloyd iteration to produce almost equilateral two dimensional meshes. Working with finite element solutions on such high quality triangulations, superconvergent recovery methods become particularly effective so that asymptotically exact a posteriori error estimations can be obtained. Through a seamless integration of these techniques, a convergent adaptive procedure is developed. As demonstrated by the numerical examples, the new AFEM is capable of solving a variety of model problems and has great potential in practical applications.

In this paper the numerical solution of the two-dimensional sine-Gordon equation is studied. Split local artificial boundary conditions are obtained by the operator splitting method. Then the original problem is reduced to an initial boundary value problem on a bounded computational domain, which can be solved by the finite difference method. Several numerical examples are provided to demonstrate the effectiveness and accuracy of the proposed method, and some interesting propagation and collision behaviors of the solitary wave solutions are observed.