The immersed interface technique is incorporated into CIP method to solve one-dimensional hyperbolic equations with piecewise constant coefficients. The proposed method achieves the third order of accuracy in time and space in the vicinity of the interface where the coefficients have jump discontinuities, which is the same order of accuracy of the standard CIP scheme. Some numerical tests are given to verify the accuracy of the proposed method.

We propose a finite-difference ghost-point approach for the numerical solution of Cauchy-Navier equations in linear elasticity problems on arbitrary unbounded domains. The technique is based on a smooth coordinate transformation, which maps an unbounded domain into a unit square. Arbitrary geometries are defined by suitable level-set functions. The equations are discretized by classical nine-point stencil on interior points, while boundary conditions and high order reconstructions are used to define the field variables at ghost-points, which are grid nodes external to the domain with a neighbor inside the domain. The linear system arising from such discretization is solved by a multigrid strategy. The approach is then applied to solve elasticity problems in volcanology for computing the displacement caused by pressure sources. The method is suitable to treat problems in which the geometry of the source often changes (explore the effects of different scenarios, or solve inverse problems in which the geometry itself is part of the unknown), since it does not require complex re-meshing when the geometry is modified. Several numerical tests are successfully performed, which asses the effectiveness of the present approach.

Nonlinear filter problems arise in many applications such as communications and signal processing. Commonly used numerical simulation methods include Kalman filter method, particle filter method, etc. In this paper a novel numerical algorithm is constructed based on samples of the current state obtained by solving the state equation implicitly. Numerical experiments demonstrate that our algorithm is more accurate than the Kalman filter and more stable than the particle filter.

The unsteady flow separation of airfoil with a local flexible structure (LFS) is studied numerically in Lagrangian frames in detail, in order to investigate the nature of its high aerodynamic performance. For such aeroelastic system, the characteristic-based split (CBS) scheme combined with arbitrary Lagrangian-Eulerian (ALE) framework is developed firstly for the numerical analysis of unsteady flow, and Galerkin method is used to approach the flexible structure. The local flexible skin of airfoil, which can lead to self-induced oscillations, is considered as unsteady perturbation to the flow. Then, the ensuing high aerodynamic performances and complex unsteady flow separation at low Reynolds number are studied by Lagrangian coherent structures (LCSs). The results show that the LFS has a significant influence on the unsteady flow separation, which is the key point for the lift enhancement. Specifically, the oscillations of the LFS can induce the generations of moving separation and vortex, which can enhance the kinetic energy transport from main flow to the boundary layer. The results could give a deep understand of the dynamics in unsteady flow separation and flow control for the flow over airfoil.

We have developed efficient numerical algorithms for solving 3D steady-state Poisson-Nernst-Planck (PNP) equations with excess chemical potentials described by the classical density functional theory (cDFT). The coupled PNP equations are discretized by a finite difference scheme and solved iteratively using the Gummel method with relaxation. The Nernst-Planck equations are transformed into Laplace equations through the Slotboom transformation. Then, the algebraic multigrid method is applied to efficiently solve the Poisson equation and the transformed Nernst-Planck equations. A novel strategy for calculating excess chemical potentials through fast Fourier transforms is proposed, which reduces computational complexity from O(N2) to O(NlogN), where N is the number of grid points. Integrals involving the Dirac delta function are evaluated directly by coordinate transformation, which yields more accurate results compared to applying numerical quadrature to an approximated delta function. Numerical results for ion and electron transport in solid electrolyte for lithiumion (Li-ion) batteries are shown to be in good agreement with the experimental data and the results from previous studies.

The bubble packing method can generate high-quality node sets in simple and complex domains. However, its efficiency remains to be improved. This study is a part of an ongoing effort to introduce several acceleration schemes to reduce the cost of simulation. Firstly, allow the viscosity coefficient c in the bubble governing equations to change according the coordinate of the bubble which are defined separately as odd and normal bubbles, and meanwhile with the saw-shape relationship with time or iterations. Then, in order to relieve the over crowded initial bubble placement, two coefficients w1 and w2 are introduced to modify the insertion criterion. The range of those two coefficients are discussed to be w1 = 1, w2 ∈ [0.5,0.8]. Finally, a self-adaptive termination condition is logically set when the stable system equilibrium is achieved. Numerical examples illustrate that the computing cost can significantly decrease by roughly 80% via adopting various combination of proper schemes (except the uniform placement example), and the average qualities of corresponding Delaunay triangulation substantially exceed 0.9. It shows that those strategies are efficient and can generate a node set with high quality.

In this paper, we investigate the coupling of the Multi-dimensional Optimal Order Detection (MOOD) method and the Arbitrary high order DERivatives (ADER) approach in order to design a new high order accurate, robust and computationally efficient Finite Volume (FV) scheme dedicated to solve nonlinear systems of hyperbolic conservation laws on unstructured triangular and tetrahedral meshes in two and three space dimensions, respectively. The Multi-dimensional Optimal Order Detection (MOOD) method for 2D and 3D geometries has been introduced in a recent series of papers for mixed unstructured meshes. It is an arbitrary high-order accurate Finite Volume scheme in space, using polynomial reconstructions with a posteriori detection and polynomial degree decrementing processes to deal with shock waves and other discontinuities. In the following work, the time discretization is performed with an elegant and efficient one-step ADER procedure. Doing so, we retain the good properties of the MOOD scheme, that is to say the optimal high-order of accuracy is reached on smooth solutions, while spurious oscillations near singularities are prevented. The ADER technique permits not only to reduce the cost of the overall scheme as shown on a set of numerical tests in 2D and 3D, but it also increases the stability of the overall scheme. A systematic comparison between classical unstructured ADER-WENO schemes and the new ADER-MOOD approach has been carried out for high-order schemes in space and time in terms of cost, robustness, accuracy and efficiency. The main finding of this paper is that the combination of ADER with MOOD generally outperforms the one of ADER and WENO either because at given accuracy MOOD is less expensive (memory and/or CPU time), or because it is more accurate for a given grid resolution. A large suite of classical numerical test problems has been solved on unstructured meshes for three challenging multi-dimensional systems of conservation laws: the Euler equations of compressible gas dynamics, the classical equations of ideal magneto-Hydrodynamics (MHD) and finally the relativistic MHD equations (RMHD), which constitutes a particularly challenging nonlinear system of hyperbolic partial differential equation. All tests are run on genuinely unstructured grids composed of simplex elements.

In this article we detail the methodology developed to construct an efficient interface description technique — the robust conservative level set (RCLS) — to simulate multiphase flows on mixed-element unstructured meshes while conserving mass to machine accuracy. The approach is tailored specifically for industry as the three-dimensional unstructured approach allows for the treatment of very complex geometries. In addition, special care has been taken to optimise the trade-off between accuracy and computational cost while maintaining the robustness of the numerical method. This was achieved by solving the transport equations for the liquid volume fraction using a WENO scheme for polyhedral meshes and by adding a flux-limiter algorithm. The performance of the resulting method has been compared against established multiphase numerical methods and its ability to capture the physics of multiphase flows is demonstrated on a range of relevant test cases. Finally, the RCLS method has been applied to the simulation of the primary breakup of a flat liquid sheet of kerosene in co-flowing high-pressure gas. This quasi-DNS/LES computation was performed at relevant aero-engine conditions on a three-dimensional mixed-element unstructured mesh. The numerical results have been validated qualitatively against theoretical predictions and experimental data. In particular, the expected breakup regime was observed in the simulation results. Finally, the computation reproduced faithfully the breakup length predicted by a correlation based on experimental data. This constitutes a first step towards a quantitative validation.

A stochastic approach to conditional simulation of flow in randomly heterogeneous media is proposed with the combination of the Karhunen-Loeve expansion and the probabilistic collocation method (PCM). The conditional log hydraulic conductivity field is represented with the Karhunen-Loeve expansion, in terms of some deterministic functions and a set of independent Gaussian random variables. The propagation of uncertainty in the flow simulations is carried out through the PCM, which relies on the efficient polynomial chaos expansion used to represent the flow responses such as the hydraulic head. With the PCM, existing flow simulators can be employed for uncertainty quantification of flow in heterogeneous porous media when direct measurements of hydraulic conductivity are taken into consideration. With illustration of several numerical examples of groundwater flow, this study reveals that the proposed approach is able to accurately quantify uncertainty of the flow responses conditioning on hydraulic conductivity data, while the computational efforts are significantly reduced in comparison to the Monte Carlo simulations.

In this paper, the mechanisms of material removal in chemical mechanical polishing (CMP) processes are investigated in detail by the smoothed particle hydrodynamics (SPH) method. The feature-scale behaviours of slurry flow, rough pad, wafer defects, moving solid boundaries, slurry-abrasive interactions, and abrasive collisions are modelled and simulated. Compared with previous work on CMP simulations, our simulations incorporate more realistic physical aspects of the CMP process, especially the effect of abrasive concentration in the slurry flows. The preliminary results on slurry flow in CMP provide microscopic insights on the experimental data of the relation between the removal rate and abrasive concentration and demonstrate that SPH is a suitable method for the research of CMP processes.

We study the computation of ground states and time dependent solutions of the Schrödinger-Poisson system (SPS) on a bounded domain in 2D (i.e. in two space dimensions). On a disc-shaped domain, we derive exact artificial boundary conditions for the Poisson potential based on truncated Fourier series expansion in θ, and propose a second order finite difference scheme to solve the r-variable ODEs of the Fourier coefficients. The Poisson potential can be solved within (MNlogN) arithmetic operations where M,N are the number of grid points in r-direction and the Fourier bases. Combined with the Poisson solver, a backward Euler and a semi-implicit/leap-frog method are proposed to compute the ground state and dynamics respectively. Numerical results are shown to confirm the accuracy and efficiency. Also we make it clear that backward Euler sine pseudospectral (BESP) method in [33] can not be applied to 2D SPS simulation.

In this paper we present an efficient algorithm for the calculation of photonic crystal band structures and band structures of photonic crystal waveguides. Our method relies on the fact that the dispersion curves of the band structure are smooth functions of the quasi-momentum in the one-dimensional Brillouin zone. We show the derivation and computation of the group velocity, the group velocity dispersion, and any higher derivative of the dispersion curves. These derivatives are then employed in a Taylor expansion of the dispersion curves. We control the error of the Taylor expansion with the help of a residual estimate and introduce an adaptive scheme for the selection of nodes in the one-dimensional Brillouin zone at which we solve the underlying eigenvalue problem and compute the derivatives of the dispersion curves. These derivatives are then employed in a Taylor expansion of the dispersion curves. We control the error of the Taylor expansion with the help of a residual estimate and introduce an adaptive scheme for the selection of nodes in the one-dimensional Brillouin zone at which we solve the underlying eigenvalue problem and compute the derivatives of the dispersion curves. The proposed algorithm is not only advantageous as it decreases the computational effort to compute the band structure but also because it allows for the identification of crossings and anti-crossings of dispersion curves, respectively. This identification is not possible with the standard approach of solving the underlying eigenvalue problem at a discrete set of values of the quasi-momentum without taking the mode parity into account.

This paper presents a modified Loop’s subdivision algorithm for studying the deformation of a single capsule, the hydrodynamic interaction between two capsules and the hydrodynamic diffusion of a suspension of capsules in bounded shear flow. A subdivision thin-shell model is employed to compute the forces generated on the surface of the elastic capsule during deformation. The capsule surface is approximated using the modified Loop’s subdivision scheme which guarantees bounded curvature and C1 continuity everywhere on the limit surface. The present numerical technique has been validated by studying the deformation of a spherical capsule in shear flow. Computations are also performed for a biconcave capsule over a wide range of shear rates and viscosity ratios to investigate its dynamics. In addition, the hydrodynamic interaction between two elastic capsules in bounded shear flow is studied. Depending on the wall separation distance, the trajectory-bifurcation points that separate reversing and crossing motions for both spherical and biconcave capsules can be found. Compared to the spherical capsules, the biconcave capsules exhibit additional types of interaction such as rotation and head-on collision. The head-on collision results in a large trajectory shift which contribute to the hydrodynamic diffusion of a suspension. A suspension of a large number of biconcave capsules in shear flow is also simulated to show the ability of the modified scheme in running a large-scale simulation over a long period of time.

In this work, we are concerned with a time-splitting Fourier pseudospectral (TSFP) discretization for the Klein-Gordon (KG) equation, involving a dimensionless parameter ɛ ∊ (0,1]. In the nonrelativistic limit regime, the small ɛ produces high oscillations in exact solutions with wavelength of (ɛ2) in time. The key idea behind the TSFP is to apply a time-splitting integrator to an equivalent first-order system in time, with both the nonlinear and linear subproblems exactly integrable in time and, respectively, Fourier frequency spaces. The method is fully explicit and time reversible. Moreover, we establish rigorously the optimal error bounds of a second-order TSFP for fixed ɛ = (1), thanks to an observation that the scheme coincides with a type of trigonometric integrator. As the second task, numerical studies are carried out, with special efforts made to applying the TSFP in the nonrelativistic limit regime, which are geared towards understanding its temporal resolution capacity and meshing strategy for (ɛ2)-oscillatory solutions when 0 < ɛ « 1. It suggests that the method has uniform spectral accuracy in space, and an asymptotic (ɛ−2Δt2) temporal discretization error bound (Δt refers to time step). On the other hand, the temporal error bounds for most trigonometric integrators, such as the well-established Gautschi-type integrator in, are (ɛ−4Δt2). Thus, our method offers much better approximations than the Gautschi-type integrator in the highly oscillatory regime. These results, either rigorous or numerical, are valid for a splitting scheme applied to the classical relativistic NLS reformulation as well.

This paper aims to study the numerical features of a coupling scheme between the immersed boundary (IB) method and the lattice Boltzmann BGK (LBGK) model by four typical test problems: the relaxation of a circular membrane, the shearing flow induced by a moving fiber in the middle of a channel, the shearing flow near a non-slip rigid wall, and the circular Couette flow between two inversely rotating cylinders. The accuracy and robustness of the IB-LBGK coupling scheme, the performances of different discrete Dirac delta functions, the effect of iteration on the coupling scheme, the importance of the external forcing term treatment, the sensitivity of the coupling scheme to flow and boundary parameters, the velocity slip near non-slip rigid wall, and the origination of numerical instabilities are investigated in detail via the four test cases. It is found that the iteration in the coupling cycle can effectively improve stability, the introduction of a second-order forcing term in LBGK model is crucial, the discrete fiber segment length and the orientation of the fiber boundary obviously affect accuracy and stability, and the emergence of both temporal and spatial fluctuations of boundary parameters seems to be the indication of numerical instability. These elaborate results shed light on the nature of the coupling scheme and may benefit those who wish to use or improve the method.

The construction of symplectic numerical schemes for stochastic Hamiltonian systems is studied. An approach based on generating functions method is proposed to generate the stochastic symplectic integration of any desired order. In general the proposed symplectic schemes are fully implicit, and they become computationally expensive for mean square orders greater than two. However, for stochastic Hamiltonian systems preserving Hamiltonian functions, the high-order symplectic methods have simpler forms than the explicit Taylor expansion schemes. A theoretical analysis of the convergence and numerical simulations are reported for several symplectic integrators. The numerical case studies confirm that the symplectic methods are efficient computational tools for long-term simulations.

In this paper, we develop and study numerical methods for the two-mode shallow water equations recently proposed in [S. STECHMANN, A. MAJDA, and B. KHOUIDER, Theor. Comput. Fluid Dynamics, 22 (2008), pp. 407-432]. Designing a reliable numerical method for this system is a challenging task due to its conditional hyperbolicity and the presence of nonconservative terms. We present several numerical approaches—two operator splitting methods (based on either Roe-type upwind or central-upwind scheme), a central-upwind scheme and a path-conservative central-upwind scheme—and test their performance in a number of numerical experiments. The obtained results demonstrate that a careful numerical treatment of nonconservative terms is crucial for designing a robust and highly accurate numerical method.

This paper focuses on the development of an efficient, three-dimensional, thermo-mechanical, nonlinear-Stokes flow computational model for ice sheet simulation. The model is based on the parallel finite element model developed in [14] which features high-order accurate finite element discretizations on variable resolution grids. Here, we add an improved iterative solution method for treating the nonlinearity of the Stokes problem, a new high-order accurate finite element solver for the temperature equation, and a new conservative finite volume solver for handling mass conservation. The result is an accurate and efficient numerical model for thermo-mechanical glacier and ice-sheet simulations. We demonstrate the improved efficiency of the Stokes solver using the ISMIP-HOM Benchmark experiments and a realistic test case for the Greenland ice-sheet. We also apply our model to the EISMINT-II benchmark experiments and demonstrate stable thermo-mechanical ice sheet evolution on both structured and unstructured meshes. Notably, we find no evidence for the “cold spoke” instabilities observed for these same experiments when using finite difference, shallow-ice approximation models on structured grids.

We investigate the critical nuclei morphology in phase transformation by combining two effective ingredients, with the first being the phase field modeling of the relevant energetics which has been a popular approach for phase transitions and the second being shrinking dimer dynamics and its variants for computing saddle points and transition states. In particular, the newly formulated generalized shrinking dimer dynamics is proposed by adopting the Cahn-Hilliard dynamics for the generalized gradient system. As illustrations, a couple of typical cases are considered, including a generic system modeling heterogeneous nucleation and a specific material system modeling the precipitate nucleation in FeCr alloys. While the standard shrinking dimer dynamics can be applied to study the non-conserved case of generic heterogeneous nucleation directly, the generalized shrinking dimer dynamics is efficient to compute precipitate nucleation in FeCr alloys due to the conservation of concentration. Numerical simulations are provided to demonstrate both the complex morphology associated with nucleation events and the effectiveness of generalized shrinking dimer dynamics based on phase field models.

In Markov Chain Monte Carlo (MCMC) simulations, thermal equilibria quantities are estimated by ensemble average over a sample set containing a large number of correlated samples. These samples are selected in accordance with the probability distribution function, known from the partition function of equilibrium state. As the stochastic error of the simulation results is significant, it is desirable to understand the variance of the estimation by ensemble average, which depends on the sample size (i.e., the total number of samples in the set) and the sampling interval (i.e., cycle number between two consecutive samples). Although large sample sizes reduce the variance, they increase the computational cost of the simulation. For a given CPU time, the sample size can be reduced greatly by increasing the sampling interval, while having the corresponding increase in variance be negligible if the original sampling interval is very small. In this work, we report a few general rules that relate the variance with the sample size and the sampling interval. These results are observed and confirmed numerically. These variance rules are derived for the MCMC method but are also valid for the correlated samples obtained using other Monte Carlo methods. The main contribution of this work includes the theoretical proof of these numerical observations and the set of assumptions that lead to them.