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.

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.

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 C_D that influences each particle in the flow. This coefficient is given by C_D =(1−φ)⁻ⁿC^∗_D, where φ is the volume fraction of particles, C^∗_D is the drag coefficient for a single particle, and n ≃ 4.67 for creeping flow. We obtain an approximation for CD/C^∗_D 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.

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.

The dynamics of the Liesegang type pattern formation is investigated in a centrally symmetric two-dimensional setup. According to the observations in real experiments, the qualitative change of the dynamics is exhibited for slightly different initial conditions. Two kinds of chemical mechanisms are studied; in both cases the pattern formation is described using a phase separation model including the Cahn-Hilliard equations. For the numerical simulations we make use of an adaptive grid PDE method, which successfully deals with the computationally critical cases such as steep gradients in the concentration distribution and investigation of long time behavior. The numerical simulations show a good agreement with the real experiments.

The multi-mode approximation is presented to compute the interior wave function of Schrödinger equation. This idea is necessary to handle the multi-barrier and high dimensional resonant tunneling problems where multiple eigenvalues are considered. The accuracy and efficiency of this algorithm is demonstrated via several numerical examples.

A method for gyrokinetic simulation of low frequency (lower than the cyclotron frequency) magnetic compressional modes in general geometry is presented. The gyrokinetic-Maxwell system of equations is expressed fully in terms of the compressional component of the magnetic perturbation, δB_‖, with finite Larmor radius effects. This introduces a "gyro-surface" averaging of δB_‖ in the gyrocenter equations of motion, and similarly in the perpendicular Ampere's law, which takes the form of the perpendicular force balance equation. The resulting system can be numerically implemented by representing the gyro-surface averaging by a discrete sum in the configuration space. For the typical wavelength of interest (on the order of the gyroradius), the gyro-surface averaging can be reduced to averaging along an effective gyro-orbit. The phase space integration in the force balance equation can be approximated by summing over carefully chosen samples in the magnetic moment coordinate, allowing for an efficient numerical implementation.

We investigate the critical properties of the Ising S=1/2 and S=1 model on (3,4,6,4) and (3⁴,6) Archimedean lattices. The system is studied through the extensive Monte Carlo simulations. We calculate the critical temperature as well as the critical point exponents γ/ν, β/ν, and ν basing on finite size scaling analysis. The calculated values of the critical temperature for S=1 are k_BT_C/J=1.590(3), and k_BT_C/J=2.100(4) for (3,4,6,4) and (3⁴,6) Archimedean lattices, respectively. The critical exponents β/ν, γ/ν, and 1/ν, for S=1 are β/ν=0.180(20), γ/ν=1.46(8), and 1/ν=0.83(5), for (3,4,6,4) and 0.103(8), 1.44(8), and 0.94(5), for (3⁴,6) Archimedean lattices. Obtained results differ from the Ising S = 1/2 model on (3,4,6,4), (3⁴,6) and square lattice. The evaluated effective dimensionality of the system for S =1 are D_eff =1.82(4), for (3,4,6,4), and D_eff=1.64(5) for (3⁴,6).

A lattice Kinetic Monte Carlo (KMC) approach is considered to study the statistical properties of the diffusion of Pt atom clusters on a Pt(111) surface. The interatomic potential experienced by the diffusing atoms is calculated by the embedded atom method and the hopping barrier for the allowed atomic movements are calculated using the Nudged Elastic Band method. The diffusion coefficient is computed for various cluster sizes and system temperatures. The obtained results are in agreement with the ones obtained in previous experimental and theoretical works. A simple scaling argument is proposed for the size dependence of the diffusion coefficient's prefactor. A detailed statistical analysis of the event by event KMC dynamics reveals two important and co-existing mechanisms for the diffusion of the cluster's center of mass. At low temperatures (below T=400K) the dominating mechanism responsible for the displacement of the cluster's center of mass is the periphery (or edge) diffusion of the atoms. At high temperatures (above T=800K) the dissociation and recombination of the clusters becomes more and more important.

The aim of the present work is to develop a general formalism to derive staggered discretizations for Lagrangian hydrodynamics on two-dimensional unstructured grids. To this end, we make use of the compatible discretization that has been initially introduced by E. J. Caramana et al., in J. Comput. Phys., 146 (1998). Namely, momentum equation is discretized by means of subcell forces and specific internal energy equation is obtained using total energy conservation. The main contribution of this work lies in the fact that the subcell force is derived invoking Galilean invariance and thermodynamic consistency. That is, we deduce a general form of the sub-cell force so that a cell entropy inequality is satisfied. The subcell force writes as a pressure contribution plus a tensorial viscous contribution which is proportional to the difference between the nodal velocity and the cell-centered velocity. This cell-centered velocity is a supplementary degree of freedom that is solved by means of a cell-centered approximate Riemann solver. To satisfy the second law of thermodynamics, the local subcell tensor involved in the viscous part of the subcell force must be symmetric positive definite. This subcell tensor is the cornerstone of the scheme. One particular expression of this tensor is given. A high-order extension of this discretization is provided. Numerical tests are presented in order to assess the efficiency of this approach. The results obtained for various representative configurations of one- and two-dimensional compressible fluid flows show the robustness and the accuracy of this scheme.

In this paper, we propose a new conservative semi-Lagrangian (SL) finite difference (FD) WENO scheme for linear advection equations, which can serve as a base scheme for the Vlasov equation by Strang splitting [4]. The reconstruction procedure in the proposed SL FD scheme is the same as the one used in the SL finite volume (FV) WENO scheme [3]. However, instead of inputting cell averages and approximate the integral form of the equation in a FV scheme, we input point values and approximate the differential form of equation in a FD spirit, yet retaining very high order (fifth order in our experiment) spatial accuracy. The advantage of using point values, rather than cell averages, is to avoid the second order spatial error, due to the shearing in velocity (v) and electrical field (E) over a cell when performing the Strang splitting to the Vlasov equation. As a result, the proposed scheme has very high spatial accuracy, compared with second order spatial accuracy for Strang split SL FV scheme for solving the Vlasov-Poisson (VP) system. We perform numerical experiments on linear advection, rigid body rotation problem; and on the Landau damping and two-stream instabilities by solving the VP system. For comparison, we also apply (1) the conservative SL FD WENO scheme, proposed in [22] for incompressible advection problem, (2) the conservative SL FD WENO scheme proposed in [21] and (3) the non-conservative version of the SL FD WENO scheme in [3] to the same test problems. The performances of different schemes are compared by the error table, solution resolution of sharp interface, and by tracking the conservation of physical norms, energies and entropies, which should be physically preserved.

In this report, a charge preserving numerical resolution of the 1D Vlasov-Ampère equation is achieved, with a forward Semi-Lagrangian method introduced in [10]. The Vlasov equation belongs to the kinetic way of simulating plasmas evolution, and is coupled with the Poisson's equation, or equivalently under charge conservation, the Ampère's one, which self-consistently rules the electric field evolution. In order to ensure having proper physical solutions, it is necessary that the scheme preserves charge numerically. B-spline deposition will be used for the interpolation step. The solving of the characteristics will be made with a Runge-Kutta 2 method and with a Cauchy-Kovalevsky procedure.

The classical discrete element approach (DEM) based on Newtonian dynamics can be divided into two major groups, event-driven methods (EDM) and time-driven methods (TDM). Generally speaking, TDM simulations are suited for cases with high volume fractions where there are collisions between multiple objects. EDM simulations are suited for cases with low volume fractions from the viewpoint of CPU time. A method combining EDM and TDM called Hybrid Algorithm of event-driven and time-driven methods (HAET) is presented in this paper. The HAET method employs TDM for the areas with high volume fractions and EDM for the remaining areas with low volume fractions. It can decrease the CPU time for simulating granular flows with strongly non-uniform volume fractions. In addition, a modified EDM algorithm using a constant time as the lower time step limit is presented. Finally, an example is presented to demonstrate the hybrid algorithm.

We propose a fast local level set method for the inverse problem of gravimetry. The theoretical foundation for our approach is based on the following uniqueness result: if an open set D is star-shaped or x₃-convex with respect to its center of gravity, then its exterior potential uniquely determines the open set D. To achieve this purpose constructively, the first challenge is how to parametrize this open set D as its boundary may have a variety of possible shapes. To describe those different shapes we propose to use a level-set function to parametrize the unknown boundary of this open set. The second challenge is how to deal with the issue of partial data as gravimetric measurements are only made on a part of a given reference domain Ω. To overcome this difficulty, we propose a linear numerical continuation approach based on the single layer representation to find potentials on the boundary of some artificial domain containing the unknown set D. The third challenge is how to speed up the level set inversion process. Based on some features of the underlying inverse gravimetry problem such as the potential density being constant inside the unknown domain, we propose a novel numerical approach which is able to take advantage of these features so that the computational speed is accelerated by an order of magnitude. We carry out numerical experiments for both two- and three-dimensional cases to demonstrate the effectiveness of the new algorithm.

The spatial domain of Molecular Dynamics simulations is usually a regular box that can be easily divided in subdomains for parallel processing. Recent efforts aimed at simulating complex biological systems, like the blood flow inside arteries, require the execution of Parallel Molecular Dynamics (PMD) in vessels that have, by nature, an irregular shape. In those cases, the geometry of the domain becomes an additional input parameter that directly influences the outcome of the simulation. In this paper we discuss the problems due to the parallelization of MD in complex geometries and show an efficient and general method to perform MD in irregular domains.