Three-dimensional matched interface and boundary (MIB) method for treating geometric singularities

This paper reports the three-dimensional (3D) generalization of our previous 2D higher-order matched interface and boundary (MIB) method for solving elliptic equations with discontinuous coefficients and non-smooth interfaces. New MIB algorithms that make use of two sets of interface jump conditions are proposed to remove the critical acute angle constraint of our earlier MIB scheme for treating interfaces with sharp geometric singularities, such as sharp edges, sharp wedges and sharp tips. The resulting 3D MIB schemes are of second-order accuracy for arbitrarily complex interfaces with sharp geometric singularities, of fourth-order accuracy for complex interfaces with moderate geometric singularities, and of sixth-order accuracy for curved smooth interfaces. A systematical procedure is introduced to make the MIB matrix optimally symmetric and banded by appropriately choosing auxiliary grid points. Consequently, the new MIB linear algebraic equations can be solved with fewer number of iterations. The proposed MIB method makes use of Cartesian grids, standard finite difference schemes, lowest order interface jump conditions and fictitious values. The interface jump conditions are enforced at each intersecting point of the interface and mesh lines to overcome the staircase phenomena in finite difference approximation. While a pair of fictitious values are determined along a mesh at a time, an iterative procedure is proposed to determine all the required fictitious values for higher-order schemes by repeatedly using the lowest order jump conditions. A variety of MIB techniques are developed to overcome geometric constraints. The essential strategy of the MIB method is to locally reduce a 2D or a 3D interface problem into 1D-like ones. The proposed MIB method is extensively validated in terms of the order of accuracy, the speed of convergence, the number of iterations and CPU time. Numerical experiments are carried out to complex interfaces, including the molecular surfaces of a protein, a missile interface, and van der Waals surfaces of intersecting spheres.     

MIB method for elliptic equations with multi-material interfaces

Elliptic partial differential equations (PDEs) are widely used to model real-world problems. Due to the heterogeneous characteristics of many naturally occurring materials and man-made structures, devices, and equipments, one frequently needs to solve elliptic PDEs with discontinuous coefficients and singular sources. The development of high-order elliptic interface schemes has been an active research field for decades. However, challenges remain in the construction of high-order schemes and particularly, for nonsmooth interfaces, i.e., interfaces with geometric singularities. The challenge of geometric singularities is amplified when they are originated from two or more material interfaces joining together or crossing each other. High-order methods for elliptic equations with multi-material interfaces have not been reported in the literature to our knowledge. The present work develops matched interface and boundary (MIB) method based schemes for solving two-dimensional (2D) elliptic PDEs with geometric singularities of multi-material interfaces. A number of new MIB schemes are constructed to account for all possible topological variations due to two-material interfaces. The geometric singularities of three-material interfaces are significantly more difficult to handle. Three new MIB schemes are designed to handle a variety of geometric situations and topological variations, although not all of them. The performance of the proposed new MIB schemes is validated by numerical experiments with a wide range of coefficient contrasts, geometric singularities, and solution types. Extensive numerical studies confirm the designed second order accuracy of the MIB method for multi-material interfaces, including a case where the derivative of the solution diverges.     

Multiscale molecular dynamics using the matched interface and boundary method

The Poisson-Boltzmann (PB) equation is an established multiscale model for electrostatic analysis of biomolecules and other dielectric systems. PB based molecular dynamics (MD) approach has a potential to tackle large biological systems. Obstacles that hinder the current development of PB based MD methods are concerns in accuracy, stability, efficiency and reliability. The presence of complex solvent-solute interface, geometric singularities and charge singularities leads to challenges in the numerical solution of the PB equation and electrostatic force evaluation in PB based MD methods. Recently, the matched interface and boundary (MIB) method has been utilized to develop the first second order accurate PB solver that is numerically stable in dealing with discontinuous dielectric coefficients, complex geometric singularities and singular source charges. The present work develops the PB based MD approach using the MIB method. New formulation of electrostatic forces is derived to allow the use of sharp molecular surfaces. Accurate reaction field forces are obtained by directly differentiating the electrostatic potential. Dielectric boundary forces are evaluated at the solvent-solute interface using an accurate Cartesian-grid surface integration method. The electrostatic forces located at reentrant surfaces are appropriately assigned to related atoms. Extensive numerical tests are carried out to validate the accuracy and stability of the present electrostatic force calculation. The new PB based MD method is implemented in conjunction with the AMBER package. MIB based MD simulations of biomolecules are demonstrated via a few example systems.     

Fast finite difference solvers for singular solutions of the elliptic Monge–Ampère equation

The elliptic Monge–Ampère equation is a fully nonlinear Partial Differential Equation which originated in geometric surface theory, and has been applied in dynamic meteorology, elasticity, geometric optics, image processing and image registration. Solutions can be singular, in which case standard numerical approaches fail.In this article we build a finite difference solver for the Monge–Ampère equation, which converges even for singular solutions. Regularity results are used to select a priori between a stable, provably convergent monotone discretization and an accurate finite difference discretization in different regions of the computational domain. This allows singular solutions to be computed using a stable method, and regular solutions to be computed more accurately. The resulting nonlinear equations are then solved by Newton’s method.Computational results in two and three-dimensions validate the claims of accuracy and solution speed. A computational example is presented which demonstrates the necessity of the use of the monotone scheme near singularities.     

Time-Domain Numerical Solutions of Maxwell Interface Problems with Discontinuous Electromagnetic Waves

This paper is devoted to time domain numerical solutions of two-dimensional (2D) material interface problems governed by the transverse magnetic (TM) and transverse electric (TE) Maxwell's equations with discontinuous electromagnetic solutions. Due to the discontinuity in wave solutions across the interface, the usual numerical methods will converge slowly or even fail to converge. This calls for the development of advanced interface treatments for popular Maxwell solvers. We will investigate such interface treatments by considering two typical Maxwell solvers – one based on collocation formulation and the other based on Galerkin formulation. To restore the accuracy reduction of the collocation finite-difference time-domain (FDTD) algorithm near an interface, the physical jump conditions relating discontinuous wave solutions on both sides of the interface must be rigorously enforced. For this purpose, a novel matched interface and boundary (MIB) scheme is proposed in this work, in which new jump conditions are derived so that the discontinuous and staggered features of electric and magnetic field components can be accommodated. The resulting MIB time-domain (MIBTD) scheme satisfies the jump conditions locally and suppresses the staircase approximation errors completely over the Yee lattices. In the discontinuous Galerkin time-domain (DGTD) algorithm – a popular Galerkin Maxwell solver, a proper numerical flux can be designed to accurately capture the jumps in the electromagnetic waves across the interface and automatically preserves the discontinuity in the explicit time integration. The DGTD solution to Maxwell interface problems is explored in this work, by considering a nodal based high order discontinuous Galerkin method. In benchmark TM and TE tests with analytical solutions, both MIBTD and DGTD schemes achieve the second order of accuracy in solving circular interfaces. In comparison, the numerical convergence of the MIBTD method is slightly more uniform, while the DGTD method is more flexible and robust.     

Numerical solution of differential equations with colored noise

Using the general theory of numerical integration of stochastic differential equations, a constructive approach to numerical methods for a system with colored noise is proposed. Efficient methods up to the 5/2 strong order and up to the third weak order, including Runge-Kutta and implicit schemes, are presented. The algorithms are tested on the Kubo oscillator.     

A fully second order implicit/explicit time integration technique for hydrodynamics plus nonlinear heat conduction problems

We present a fully second order implicit/explicit time integration technique for solving hydrodynamics coupled with nonlinear heat conduction problems. The idea is to hybridize an implicit and an explicit discretization in such a way to achieve second order time convergent calculations. In this scope, the hydrodynamics equations are discretized explicitly making use of the capability of well-understood explicit schemes. On the other hand, the nonlinear heat conduction is solved implicitly. Such methods are often referred to as IMEX methods [2], [1], [3]. The Jacobian-Free Newton Krylov (JFNK) method (e.g. [10], [9]) is applied to the problem in such a way as to render a nonlinearly iterated IMEX method. We solve three test problems in order to validate the numerical order of the scheme. For each test, we established second order time convergence. We support these numerical results with a modified equation analysis (MEA) [21], [20]. The set of equations studied here constitute a base model for radiation hydrodynamics.     

Adaptively deformed mesh based interface method for elliptic equations with discontinuous coefficients

Mesh deformation methods are a versatile strategy for solving partial differential equations (PDEs) with a vast variety of practical applications. However, these methods break down for elliptic PDEs with discontinuous coefficients, namely, elliptic interface problems. For this class of problems, the additional interface jump conditions are required to maintain the well-posedness of the governing equation. Consequently, in order to achieve high accuracy and high order convergence, additional numerical algorithms are required to enforce the interface jump conditions in solving elliptic interface problems. The present work introduces an interface technique based adaptively deformed mesh strategy for resolving elliptic interface problems. We take the advantages of the high accuracy, flexibility and robustness of the matched interface and boundary (MIB) method to construct an adaptively deformed mesh based interface method for elliptic equations with discontinuous coefficients. The proposed method generates deformed meshes in the physical domain and solves the transformed governed equations in the computational domain, which maintains regular Cartesian meshes. The mesh deformation is realized by a mesh transformation PDE, which controls the mesh redistribution by a source term. The source term consists of a monitor function, which builds in mesh contraction rules. Both interface geometry based deformed meshes and solution gradient based deformed meshes are constructed to reduce the L(∞) and L(2) errors in solving elliptic interface problems. The proposed adaptively deformed mesh based interface method is extensively validated by many numerical experiments. Numerical results indicate that the adaptively deformed mesh based interface method outperforms the original MIB method for dealing with elliptic interface problems.     

Nonlinear development of acoustical instabilities in supersonic jets

In contrast with the roll-up of fluid interfaces through Kelvin-Helmholtz instability, recent numerical simulations with small amplitude perturbations of supersonic jets reveal another very different coherent mode of nonlinear acoustical instability of jets through the appearance of regular zig-zag shock patterns which traverse the interior of the jet and amplify as time evolves. In this paper, through a combination of appropriate ideas from linear and nonlinear high frequency geometric optics, the authors develop a quantitative theory which predicts the nonlinear development of zig-zag modes with a structure like those observed in the numerical simulations. The perturbation analysis is developed via a systematic application of nonlinear small amplitude high frequency geometric optics to the complex free surface problem defined by the perturbed jet; this procedure automatically yields simplified asymptotic equations which are analyzed explicitly and lead to the development of regular amplifying “zig-zag” shock structures in the jet. For a given streamwise period, Mach number, and jet width, the asymptotic theory gives explicit criteria for the number and structure of different regular zig-zag shock patterns which amplify with time. For Mach numbers M < 1, there are no amplifying acoustic zig-zig modes while for M > 1, there are a finite number of such modes depending on Mach number, jet width, and streamwise period. Explicit criteria to select the most destabilizing of these nonlinear eigenmodes are developed as well as several new quantitative predictions regarding the nonlinear development of acoustical instabilities in supersonic jets including the phenomenon of “super-resonance” for special values of the streamwise period.     

Application of Backward Differentiation Formulas to Neutron Kinetics Problems

The neutron kinetics equations belong to the class of stiff equations for numerical time integration schemes. In this work, the accuracy and speed of algorithms based on backward differentiation formulas (BDFs) are studied with regard to point and spatial neutron kinetics problems. Using a BDF algorithm with the automatic selection of time step and order, solutions to a number of model problems with both a positive and negative reactivity insertion are analyzed. Plots of numerical cost dependences on local and global errors are presented. The results indicate that algorithms based on BDFs are highly efficient and that their application to nuclear reactor simulation is justified.     

Sharp interface tracking using the phase-field equation

A general interface tracking method based on the phase-field equation is presented. The zero phase-field contour is used to implicitly track the sharp interface on a fixed grid. The phase-field propagation equation is derived from an interface advection equation by expressing the interface normal and curvature in terms of a hyperbolic tangent phase-field profile across the interface. In addition to normal interface motion driven by a given interface speed or by interface curvature, interface advection by an arbitrary external velocity field is also considered. In the absence of curvature-driven interface motion, a previously developed counter term is used in the phase-field equation to cancel out such motion. Various modifications of the phase-field equation, including nonlinear preconditioning, are also investigated. The accuracy of the present method is demonstrated in several numerical examples for a variety of interface motions and shapes that include singularities, such as sharp corners and topology changes. Good convergence with respect to the grid spacing is obtained. Mass conservation is achieved without the use of separate re-initialization schemes or Lagrangian marker particles. Similarities with and differences to other interface tracking approaches are emphasized.     

Efficient spectral-Galerkin methods for systems of coupled second-order equations and their applications

We construct in this paper two efficient spectral-Galerkin algorithms for solving systems of n coupled second-order equations. The computational complexity of these algorithms is essentially n times the cost of solving one second-order equation. We present numerical results which illustrate the accuracy and flexibility of these algorithms, as well as several interesting and challenging applications, including in particular a number of high-order nonlinear parabolic type equations.     

Differential geometry based solvation model I: Eulerian formulation

This paper presents a differential geometry based model for the analysis and computation of the equilibrium property of solvation. Differential geometry theory of surfaces is utilized to define and construct smooth interfaces with good stability and differentiability for use in characterizing the solvent-solute boundaries and in generating continuous dielectric functions across the computational domain. A total free energy functional is constructed to couple polar and nonpolar contributions to the salvation process. Geometric measure theory is employed to rigorously convert a Lagrangian formulation of the surface energy into an Eulerian formulation so as to bring all energy terms into an equal footing. By minimizing the total free energy functional, we derive coupled generalized Poisson-Boltzmann equation (GPBE) and generalized geometric flow equation (GGFE) for the electrostatic potential and the construction of realistic solvent-solute boundaries, respectively. By solving the coupled GPBE and GGFE, we obtain the electrostatic potential, the solvent-solute boundary profile, and the smooth dielectric function, and thereby improve the accuracy and stability of implicit solvation calculations. We also design efficient second order numerical schemes for the solution of the GPBE and GGFE. Matrix resulted from the discretization of the GPBE is accelerated with appropriate preconditioners. An alternative direct implicit (ADI) scheme is designed to improve the stability of solving the GGFE. Two iterative approaches are designed to solve the coupled system of nonlinear partial differential equations. Extensive numerical experiments are designed to validate the present theoretical model, test computational methods, and optimize numerical algorithms. Example solvation analysis of both small compounds and proteins are carried out to further demonstrate the accuracy, stability, efficiency and robustness of the present new model and numerical approaches. Comparison is given to both experimental and theoretical results in the literature.     

Invariantization of the Crank-Nicolson method for Burgers’ equation

In this article, a geometric technique to construct numerical schemes for partial differential equations (PDEs) that inherit Lie symmetries is proposed. The moving frame method enables one to adjust the numerical schemes in a geometric manner and systematically construct proper invariant versions of them. To illustrate the method, we study invariantization of the Crank-Nicolson scheme for Burgers’ equation. With careful choice of normalization equations, the invariantized schemes are shown to surpass the standard scheme, successfully removing numerical oscillation around sharp transition layers.     

Linear high-resolution schemes for hyperbolic conservation laws: TVB numerical evidence

The Osher–Chakrabarthy family of linear flux-modification schemes is considered. Improved lower bounds on the compression factors are provided, which suggest the viability of using the unlimited version. The LLF flux formula is combined with these schemes in order to obtain efficient finite-difference algorithms. The resulting schemes are applied to a battery of numerical tests, going from advection and Burgers equations to Euler and MHD equations, including the double Mach reflection and the Orszag–Tang 2D vortex problem. Total-variation-bounded (TVB) behavior is evident in all cases, even with time-independent upper bounds. The proposed schemes, however, do not deal properly with compound shocks, arising from non-convex fluxes, as shown by Buckley–Leverett test simulations.     

A class of asymptotic-preserving schemes for kinetic equations and related problems with stiff sources

In this paper, we propose a general time-discrete framework to design asymptotic-preserving schemes for initial value problem of the Boltzmann kinetic and related equations. Numerically solving these equations are challenging due to the nonlinear stiff collision (source) terms induced by small mean free or relaxation time. We propose to penalize the nonlinear collision term by a BGK-type relaxation term, which can be solved explicitly even if discretized implicitly in time. Moreover, the BGK-type relaxation operator helps to drive the density distribution toward the local Maxwellian, thus naturally imposes an asymptotic-preserving scheme in the Euler limit. The scheme so designed does not need any nonlinear iterative solver or the use of Wild Sum. It is uniformly stable in terms of the (possibly small) Knudsen number, and can capture the macroscopic fluid dynamic (Euler) limit even if the small scale determined by the Knudsen number is not numerically resolved. It is also consistent to the compressible Navier–Stokes equations if the viscosity and heat conductivity are numerically resolved. The method is applicable to many other related problems, such as hyperbolic systems with stiff relaxation, and high order parabolic equations.     

High order matched interface and boundary methods for the Helmholtz equation in media with arbitrarily curved interfaces

This work overcomes the difficulty of the previous matched interface and boundary (MIB) method in dealing with interfaces with non-constant curvatures for optical waveguide analysis. This difficulty is essentially bypassed by avoiding the use of local cylindrical coordinates in the improved MIB method. Instead, novel jump conditions are derived along global Cartesian directions for the transverse magnetic field components. Effective interface treatments are proposed to rigorously impose jump conditions across arbitrarily curved interfaces based on a simple Cartesian grid. Even though each field component satisfies the scalar Helmholtz equation, the enforcement of jump conditions couples two transverse magnetic field components, so that the resulting MIB method is a full-vectorial approach for the modal analysis of optical waveguides. The numerical performance of the proposed MIB method is investigated by considering interface problems with both constant and general curvatures. The MIB method is shown to be able to deliver a fourth order of accuracy in all cases, even when a high frequency solution is involved.     

基于联合平行投影迭代的叠栅层析重建

描述了一种新的叠栅层析代数迭代模型,并针对传统代数迭代法中的垂直投影算法收敛速度较慢的问题,结合新模型提出了一种改进了的平行投影重建算法.用该算法进行了数值模拟重建,结果表明：与叠栅层析中的变换类算法相比,对非完全数据问题,新算法具有有效结合先验知识进行重建的能力;与传统的垂直投影算法相比,新算法能在保证重建准确度的前提下大幅度的提高收敛速度.     

界面不稳定性数值模拟中的虚拟流动方法

研究了流体界面不稳定性的一类数值模拟方法——虚拟流动方法(Ghost Fluid Method).在算法中直接针对多维问题设定虚拟区域的流动参数,在流体力学方程的计算中采用了非分裂型的高分辨SCB格式,最后利用该方法完成了R-M和R-T不稳定性问题的数值计算,得到了满意的计算结果.     

Analysis of nonlinear dynamic response for delaminated fiber–metal laminated beam under unsteady temperature field

The nonlinear dynamic response problems of fiber–metal laminated beams with delamination are studied in this paper. Basing on the Timoshenko beam theory, and considering geometric nonlinearity, transverse shear deformation, temperature effect and contact effect, the nonlinear governing equations of motion for fiber–metal laminated beams under unsteady temperature field are established, which are solved by the differential quadrature method, Nermark-β method and iterative method. In numerical examples, the effects of delamination length, delamination depth, temperature field, geometric nonlinearity and transverse shear deformation on the nonlinear dynamic response of the glass reinforced aluminum laminated beam with delamination are discussed in details.