Abstract: We report here an efficient implementation of the finite difference Poisson–Boltzmann solvent model based on the Modified Incomplete Cholsky Conjugate Gradient algorithm, which gives rather impressive performance for both static and dynamic systems. This is achieved by implementing the algorithm with Eisenstat's two optimizations, utilizing the electrostatic update in simulations, and applying prudent approximations, including: relaxing the convergence criterion, not updating Poisson–Boltzmann-related forces every step, and using electrostatic focusing. It is also possible to markedly accelerate the supporting routines that are used to set up the calculations and to obtain energies and forces. The resulting finite difference Poisson–Boltzmann method delivers efficiency comparable to the distance-dependent dielectric model for a system tested, HIV Protease, making it a strong candidate for solution-phase molecular dynamics simulations. Further, the finite difference method includes all intrasolute electrostatic interactions, whereas the distance dependent dielectric calculations use a 15-A cutoff. The speed of our numerical finite difference method is comparable to that of the pair-wise Generalized Born approximation to the Poisson–Boltzmann method. © 2002 Wiley Periodicals, Inc. J Comput Chem 23: 1244–1253, 2002

Abstract: This work investigates the application of a high-order finite difference method for compressible large-eddy simulations on stretched, curvilinear and dynamic meshes. The solver utilizes 4th and 6th-order compact-differencing schemes for the spatial discretization, coupled with both explicit and implicit time-marching methods. Up to 10th order, Pade-type low-pass spatial filter operators are also incorporated to eliminate the spurious high-frequency modes which inevitably arise due to the lack of inherent dissipation in the spatial scheme. The solution procedure is evaluated for the case of decaying compressible isotropic turbulence and turbulent channel flow. The compact/filtering approach is found to be superior to standard second and fourth-order centered, as well as third-order upwind-biased approximations. For the case of isotropic turbulence, better results are obtained with the compact/filtering method (without an added subgrid-scale model) than with the constant-coefficient and dynamic Smagorinsky models. This is attributed to the fact that the SGS models, unlike the optimized low-pass filter, exert dissipation over a wide range of wave numbers including on some of the resolved scales

Abstract: Holes in the skull may have a large influence on the EEG and ERP. Inverse source modeling techniques such as dipole fitting require an accurate volume conductor model. This model should incorporate holes if present, especially when either a neuronal generator or the electrodes are close to the hole, e.g., in case of a trephine hole in the upper part of the skull. The boundary element method (BEM) is at present the preferred method for inverse computations using a realistic head model, because of its efficiency and availability. Using a simulation approach, we have studied the accuracy of the BEM by comparing it to the analytical solution for a volume conductor without a hole, and to the finite difference method (FDM) for one with a hole. Furthermore, we have evaluated the influence of holes on the results of forward and inverse computations using the BEM. Without a hole and compared to the analytical model, a three-sphere BEM model was accurate up to 5-10%, while the corresponding FDM model had an error <0.5%. In the presence of a hole, the difference between the BEM and the FDM was, on average, 4% (1.3-11.4%). The FDM turned out to be very accurate if no hole is present. We believe that the difference between the BEM and the FDM represents the inaccuracy of the BEM. This inaccuracy in the BEM is very small compared to the effect that holes can have on the scalp potential (up to 450%). In regard to the large influence of holes on forward and inverse computations, we conclude that holes in the skull can be treated reliably by means of the BEM and should be incorporated in forward and inverse modeling.

Abstract: An efficient finite difference model of blood flow through the coronary vessels is developed and applied to a geometric model of the largest six generations of the coronary arterial network. By constraining the form of the velocity profile across the vessel radius, the three-dimensional Navier--Stokes equations are reduced to one-dimensional equations governing conservation of mass and momentum. These equations are coupled to a pressure-radius relationship characterizing the elasticity of the vessel wall to describe the transient blood flow through a vessel segment. The two step Lax--Wendroff finite difference method is used to numerically solve these equations. The flow through bifurcations, where three vessel segments join, is governed by the equations of conservation of mass and momentum. The solution to these simultaneous equations is calculated using the multidimensional Newton--Raphson method. Simulations of blood flow through a geometric model of the coronary network are presented demonstrating phy...

Abstract: Simulation of very fast surge phenomena in a three-dimensional structure requires a method based on Maxwell's equations such as the finite difference time domain (FDTD) method or the method of moments (MoM), because circuit-equation-based methods cannot handle the phenomena. This paper presents a method of thin wire representation for the FDTD method that is suitable for the three-dimensional surge simulation. The thin wire representation is indispensable to simulate electromagnetic surges on wires or steel flames of which the radius is smaller than a discretized space step used in the FDTD simulation. Comparisons between calculated and laboratory-test results are presented to show the accuracy of the proposed thin wire representation, and the development of a general surge analysis program based on the FDTD method is also described in the present paper.

Abstract: The theoretical basis and the numerical implementation of free-vortex filament methods are reviewed for application to the prediction and analysis of helicopter rotor wakes. The governing equations for the problem are described, with a discussion of finite difference approximations to these equations and various numerical solution techniques. Both relaxation and time-marching wake solution techniques are reviewed. It is emphasized how the careful consideration of stability and convergence (grid-independent behavior) are important to ensure a physically correct wake solution. The implementation of viscous diffusion and filament straining effects are also discussed. The need for boundary condition corrections to compensate for the inevitable wake truncation are described. Algorithms to accelerate the wake solution using velocity field interpolation are shown to reduce computational costs without a loss of accuracy. Several challenging examples of the application of free-vortex filament methods to helicopter rotor problems are shown, including multirotor configurations, flight near the ground, maneuvering flight conditions, and descending flight through the vortex ring state

Abstract: The bidomain equations are the most complete description of cardiac electrical activity. Their numerical solution is, however, computationally demanding, especially in three dimensions, because of the fine temporal and spatial sampling required. This paper methodically examines computational performance when solving the bidomain equations. Several techniques to speed up this computation are examined in this paper. The first step was to recast the equations into a parabolic part and an elliptic part. The parabolic part was solved by either the finite-element method (FEM) or the interconnected cable model (ICCM). The elliptic equation was solved by FEM on a coarser grid than the parabolic problem and at a reduced frequency. The performance of iterative and direct linear equation system solvers was analyzed as well as the scalability and parallelizability of each method. Results indicate that the ICCM was twice as fast as the FEM for solving the parabolic problem, but when the total problem was considered, this resulted in only a 20% decrease in computation time. The elliptic problem could be solved on a coarser grid at one-quarter of the frequency at which the parabolic problem was solved and still maintain reasonable accuracy. Direct methods were faster than iterative methods by at least 50% when a good estimate of the extracellular potential was required. Parallelization over four processors was efficient only when the model comprised at least 500 000 nodes. Thus, it was possible to speed up solution of the bidomain equations by an order of magnitude with a slight decrease in accuracy.

Abstract: In this paper, a numerical model based on the finite difference method is presented to predict tool and chip temperature fields in continuous machining and time varying milling processes. Continuous or steady state machining operations like orthogonal cutting are studied by modeling the heat transfer between the tool and chip at the tool—rake face contact zone. The shear energy created in the primary zone, the friction energy produced at the rake face—chip contact zone and the heat balance between the moving chip and stationary tool are considered. The temperature distribution is solved using the finite difference method. Later, the model is extended to milling where the cutting is interrupted and the chip thickness varies with time. The time varying chip is digitized into small elements with differential cutter rotation angles which are defined by the product of spindle speed and discrete time intervals. The temperature field in each differential element is modeled as a first-order dynamic system, whose time constant is identified based on the thermal properties of the tool and work material, and the initial temperature at the previous chip segment. The transient temperature variation is evaluated by recursively solving the first order heat transfer problem at successive chip elements. The proposed model combines the steady-state temperature prediction in continuous machining with transient temperature evaluation in interrupted cutting operations where the chip and the process change in a discontinuous manner. The mathematical models and simulation results are in satisfactory agreement with experimental temperature measurements reported in the literature.

Abstract: We consider weak solutions of (hyperbolic or hyperbolic-elliptic) systems of conservation laws in one-space dimension and their approximation by finite difference schemes in conservative form. The systems under consideration are endowed with an entropy-entropy flux pair. We introduce a general approach to construct second and third order accurate, fully discrete (in both space and time) entropy conservative schemes. In general, these schemes are fully nonlinear implicit, but in some important cases can be explicit or linear implicit. Furthermore, semidiscrete entropy conservative schemes of arbitrary order are presented. The entropy conservative schemes are used to construct a numerical method for the computation of weak solutions containing nonclassical regularization-sensitive shock waves. Finally, specific examples are investigated and tested numerically. Our approach extends the results and techniques by Tadmor [in Numerical Methods for Compressible Flows---Finite Difference, Element and Volume Techniques, ASME, New York, 1986, pp. 149--158], LeFloch and Rohde [SIAM J. Numer. Anal., 37 (2000), pp. 2023--2060].

Abstract: A three-dimensional numerical model based on the full Navier–Stokes equations (NSE) in σ-coordinate is developed in this study. The σ-coordinate transformation is first introduced to map the irregular physical domain with the wavy free surface and uneven bottom to the regular computational domain with the shape of a rectangular prism. Using the chain rule of partial differentiation, a new set of governing equations is derived in the σ-coordinate from the original NSE defined in the Cartesian coordinate. The operator splitting method (Li and Yu, Int. J. Num. Meth. Fluids 1996; 23: 485–501), which splits the solution procedure into the advection, diffusion, and propagation steps, is used to solve the modified NSE. The model is first tested for mass and energy conservation as well as mesh convergence by using an example of water sloshing in a confined tank. Excellent agreements between numerical results and analytical solutions are obtained. The model is then used to simulate two- and three-dimensional solitary waves propagating in constant depth. Very good agreements between numerical results and analytical solutions are obtained for both free surface displacements and velocities. Finally, a more realistic case of periodic wave train passing through a submerged breakwater is simulated. Comparisons between numerical results and experimental data are promising. The model is proven to be an accurate tool for consequent studies of wave-structure interaction. Copyright © 2002 John Wiley & Sons, Ltd.

Abstract: This work is concerned with the development of a numerical method capable of simulating viscoelastic free surface flow of an Oldroyd-B fluid. The basic equations governing the flow of an Oldroyd-B fluid are considered. A novel formulation is developed for the computation of the non-Newtonian extra-stress components on rigid boundaries. The full free surface stress conditions are employed. The resulting governing equations are solved by a finite difference method on a staggered grid, influenced by the ideas of the marker-and-cell (MAC) method. Numerical results demonstrating the capabilities of this new technique are presented for a number of problems involving unsteady free surface flows.

Abstract: The point-source traveltime field has an upwind singularity at the source point. Consequently, all formally high-order, finite-difference eikonal solvers exhibit first-order convergence and relatively large errors. Adaptive upwind finite-difference methods based on high-order Weighted Essentially NonOscillatory (WENO) Runge-Kutta difference schemes for the paraxial eikonal equation overcome this difficulty. The method controls error by automatic grid refinement and coarsening based on a posteriori error estimation. It achieves prescribed accuracy at a far lower cost than does the fixed-grid method. Moreover, the achieved high accuracy of traveltimes yields reliable estimates of auxiliary quantities such as take-off angles and geometric spreading factors.

Abstract: In this paper we develop a Lax--Wendroff time discretization procedure for high order finite difference weighted essentially nonoscillatory schemes to solve hyperbolic conservation laws. This is an alternative method for time discretization to the popular TVD Runge--Kutta time discretizations. We explore the possibility in avoiding the local characteristic decompositions or even the nonlinear weights for part of the procedure, hence reducing the cost but still maintaining nonoscillatory properties for problems with strong shocks. As a result, the Lax--Wendroff time discretization procedure is more cost effective than the Runge--Kutta time discretizations for certain problems including two-dimensional Euler systems of compressible gas dynamics.

Abstract: New multigrid methods are developed for the maximum principle preserving immersed interface method applied to second order linear elliptic and parabolic PDEs that involve interfaces and discontinuities. For elliptic interface problems, the multigrid solver developed in this paper works while some other multigrid solvers do not. For linear parabolic equations, we have developed the second order maximum principle preserving finite difference scheme in this paper. We use the Crank--Nicolson scheme to deal with the diffusion part and an explicit scheme for the first order derivatives. Numerical examples are also presented.

Abstract: Finite-difference (FD) techniques are widely used to model wave propagation through complex structures. Two main sources of error can be identified: (1) from numerical dispersion and numerical anisotropy and (2) by modeling the response of internal grid boundaries. Conventional discretization criteria to reduce the effects of numerical dispersion and numerical anisotropy have long been established (5-8 gridpoints per wavelength for a fourth-order accurate FD scheme). We analyze the second source of errors, comparing different staggered-grid FD solutions to the Cagniard-de Hoop solution in models with fluid-solid contacts. Our results confirm that it is sufficient to rely on conventional discretization criteria if the fluid-solid interface is aligned with the grid. If accurate modeling of the Scholte wave is required, then a new imaging method we propose should be used to allow for conventional sampling of the wavefield to minimize numerical dispersion. However, for an interface not aligned with the grid (irregular interfaces), a spatial sampling of at least 15 gridpoints per minimum wavelength is required to obtain acceptable results, particularly in seismic seabed applications where Scholte waves may need to be modeled more accurately.

Abstract: A compact two-dimensional (2-D) full-wave finite-difference frequency-domain method is proposed for the analysis of dispersion characteristics of a general guided wave structure. Because the longitudinal field components are eliminated in the proposed method, only four transverse field components are involved in the final resulting eigen equation. This feature considerably reduces the required CPU time as compared to the existing approaches by which six field components are comprised. Additionally, unlike other 2-D finite-difference schemes that determine the eigenfrequency for a given propagation constant, the new method finds the propagation constant /spl beta/ for a given k/sub o/ (frequency). The new method has been verified by examining the computed results of a number of typical guided wave structures with the published results. Very good agreement is achieved.

Abstract: Mode multiplexing and demultiplexing devices by using multimode interference couplers are proposed. Mode separating as well as superposing based on single and multiple self-imaging are utilised. Simulation results of wave propagation by the finite difference beam propagation method demonstrate the validity of the synthesis method.

Abstract: Existing methods for the analysis of transient flows in pipe networks are often geared towards certain types of flows such as gas flows vis-a-vis liquid flows or isothermal flows vis-a-vis non-isothermal flows. Also, simplifying assumptions are often made which introduce inaccuracies when the method is applied outside the domain for which it was originally intended. This paper describes an implicit finite difference method based on the simultaneous pressure correction approach which is valid for both liquid and gas flows, for both isothermal and non-isothermal flows and for both fast and slow transients. The problematic convective acceleration term in the momentum equation, often neglected in other methods, is retained but eliminated by casting the momentum equation in an alternative form. The accuracy and stability of the method, depending on a time-step weighing factor α, are illustrated by analysing fast transients in a pipeline and simple branching network. The proposed method compares very well with the second-order-method of characteristics and the two-step Lax–Wendroff method. The advantages of the present method is its speed over a range of problems including both fast and slow transients, its accuracy, its stability and its flexibility. Copyright © 2001 John Wiley & Sons, Ltd.

Abstract: For pt.1 see ibid., vol.44, no.1, p.134-42 (2002). Higher-order schemes for the finite-difference time-domain (FDTD) method - in particular, a second-order-in-time, fourth-order-in-space method, FDTD(2,4) - are applied to a number of problems. The problems include array analysis, cavity resonances, antenna coupling, and shielding effectiveness case studies. The latter includes a simplified model of a commercial airliner, with a personal electronic device operating in the vicinity of the aircraft. The FDTD computations are also compared to measured data for this case. Incorporating PEC and other types of material boundaries into higher-order FDTD is problematic; a hybrid approach using the standard FDTD method in the proximity of the boundary is proposed, and shown to perform well.

Abstract: A methodology is presented that allows the derivation of low-truncation-error finite-difference representations or the two-dimensional Helmholtz equation, specific to waveguide analysis. This methodology is derived from the formal infinite series solution involving Bessel functions and sines and cosines. The resulting finite-difference equations are valid everywhere except at dielectric corners, and are highly accurate (from fourth to sixth order, depending on the type of grid employed). None the less, they utilize only a nine-point stencil, and thus lead to only minor increases in numerical effort compared with the standard Crank-Nicolson equations.

Abstract: The auxiliary differential equation method used for treating Lorentz media with the FDTD algorithm is reinterpreted to reduce computer memory requirements while maintaining full time synchronism. This new formulation enables us to save up to 20% of computer memory when compared to the previous approach. The validity of our formulation is proved analytically by comparing the resulting equations with those from the previous method. We also employ this approach for Debye media and show that no disadvantage arises in terms of memory requirement.