Two-dimensional Finite Difference Scheme Research Articles

Reprint of: Dirichlet-to-Neumann mappings and finite-differences for anisotropic diffusion

A general methodology, which consists in deriving two-dimensional finite-difference schemes which involve numerical fluxes based on Dirichlet-to-Neumann maps (or Steklov–Poincaré operators), is first recalled. Then, it is applied to several types of diffusion equations, some being weakly anisotropic, endowed with an external source. Standard finite-difference discretizations are systematically recovered, showing that in absence of any other mechanism, like e.g. convection and/or damping (which bring Bessel and/or Mathieu functions inside that type of numerical fluxes), these well-known schemes achieve a satisfying multi-dimensional character.

Dirichlet-to-Neumann mappings and finite-differences for anisotropic diffusion

A general methodology, which consists in deriving two-dimensional finite-difference schemes which involve numerical fluxes based on Dirichlet-to-Neumann maps (or Steklov–Poincaré operators), is first recalled. Then, it is applied to several types of diffusion equations, some being weakly anisotropic, endowed with an external source. Standard finite-difference discretizations are systematically recovered, showing that in absence of any other mechanism, like e.g. convection and/or damping (which bring Bessel and/or Mathieu functions inside that type of numerical fluxes), these well-known schemes achieve a satisfying multi-dimensional character.

Highly excited electronic image states of metallic nanorings.

We study electronic image states around a metallic nanoring and show that the interplay between the attractive polarization force and a repulsive centrifugal force gives rise to Rydberg-like image states trapped several nanometers away from the surface. The nanoring is modeled as a perfectly conducting isolated torus whose classical electrostatic image potential is derived analytically. The image states are computed via a two-dimensional finite-difference scheme as solutions of the effective Schrödinger equation describing the outer electron subject to this image potential. These findings demonstrate not only the existence of detached image states around nanorings but allow us also to provide general criteria on the ring geometry, i.e., the aspect ratio of the torus, that need to be fulfilled in order to support such states.

Dynamic Interface Fracture in a Stratified Geological Medium Accompanied by Primary Fault Rupture Propagation

Although the physical properties of fault rupture that causes a dip-slip earthquake (ordinary presumed as mode-II shear crack propagation) in “composite” or stratified media in Earth's crust may be significantly different from those in a monolithic one, fracture along interfaces is usually neglected in seismological analysis. That is, the primary fault rupture (main crack) is normally assumed to penetrate an interface and further propagate in the adjacent (similar or dissimilar) geological layer. However, as is well known in the field of fracture mechanics, in reality, the primary fault rupture may be deflected (bifurcated) at and propagated along the interface. In the static context, crack penetration and/or deflection in a bi-material system is often investigated for a main crack extending at an angle to (and until) the interface. Interfacial fracture ahead of the primary rupture front (tip) may also induce deflection, but the technical details of dynamic rupture deflection/penetration at or across an interface have not been fully clarified yet. The purpose of this study is to theoretically consider a possible seismic rupture model that takes into account dynamic interfacial deflection and penetration on (or ahead of, if any) a primary rupture front moving at an angle toward an interface in a heterogeneously stratified medium. Based on analytical dynamic fracture mechanics and utilizing a two-dimensional finite difference scheme, we illustrate the dynamic stress fields associated with the wave expansion from primary fault rupture and the wave interaction with geological features of discontinuity. The wave structures in the proximity of propagating fault rupture front in a linear elastic stratified medium show the possible effect of the Mach numbers (i.e. ratios between the rupture velocity (crack tip speed) and the velocities of the elastic waves in the surrounding medium) on the initiation and propagation of secondary cracks.

Methodology for Studying the Numerical Speed of Sound in Finite Differences Schemes

The space and time discretization inherent to all FDTD schemes introduce non-physical dispersion errors, i.e. deviations of the speed of sound from the theoretical value predicted by the governing Euler differential equations. A generalmethodology for computing this dispersion error via straightforward numerical simulations of the FDTD schemes is presented. The method is shown to provide remarkable accuracies of the order of 1/1000 in a wide variety of twodimensional finite difference schemes.

Development of least-square-based two-dimensional finite-difference schemes and their application to simulate natural convection in a cavity

In this paper, two-dimensional mesh-free finite-difference schemes for solving incompressible viscous flows are presented. The method is based on the use of a weighted least-square approximation procedure together with a Taylor series expansion of the unknown function. Discretization error for derivatives is investigated analytically on the uniform mesh and the convergence property of the method is numerically tested. The role of the weighting function playing in the method is studied. Neumann-type boundary condition is treated by applying locally orthogonal boundary grids. Application to a problem of natural convection in a cavity is demonstrated on three different types of point distribution.

Implications of seismic scattering below the 8° discontinuity along PNE profile Kraton

Strong scattered arrivals are observed in the 800–1400-km offset range of the four Peaceful Nuclear Explosion (PNE) seismic sections along the 3500-km-long E–W-trending Russian profile Kraton. They give evidence of an ∼85-km-thick inhomogeneous scattering layer in the upper mantle, which has a slightly reduced (∼1%) average P-wave velocity compared to the layers above and below. The top of the anomalous zone is the 8° discontinuity at ∼100 km depth, which is constrained by a clear travel time delay of up to 2 s, consistently observed at ∼900–1000 km (∼8°) offset in the four PNE seismic sections. We model the inhomogeneous zone below the 8° discontinuity by random velocity fluctuations. The full waveform seismic response of a suite of stochastic velocity models is calculated using a two-dimensional finite-difference scheme. The calculated seismic sections resemble the PNE data for models where the random fluctuations have a horizontal correlation length of 5–10 km, a vertical correlation length not greater than 5 km and a standard deviation of about 2% of the average background velocity. Fluctuations with significantly larger correlation lengths give rise to coherent scattered arrivals and considerably shorter correlation lengths lead to significantly weaker scattering, none of which resemble the observed data. The preferred seismic velocity model is in agreement with an upper mantle model in which the low-velocity zone below the 8° discontinuity may contain small amounts of partially molten material. Alternatively, the combined effects of rheological weakening and small-scale anisotropy may also account for the low-velocity scattering layer.

Development of a sixth-order two-dimensional convection–diffusion scheme via Cole–Hopf transformation

In this paper, we develop a two-dimensional finite-difference scheme for solving the time-dependent convection–diffusion equation. The numerical method exploits Cole–Hopf equation to transform the nonlinear scalar transport equation into the linear heat conduction equation. Within the semi-discretization context, the time derivative term in the transformed parabolic equation is approximated by a second-order accurate time-stepping scheme, resulting in an inhomogeneous Helmholtz equation. We apply the alternating direction implicit scheme of Polezhaev to solve the Helmholtz equation. As the key to success in the present simulation, we develop a Helmholtz scheme with sixth-order spatial accuracy. As is standard practice, we validated the code against test problems which were amenable to exact solutions. Results show excellent agreement for the one-dimensional test problems and good agreement with the analytical solution for the two-dimensional problem.

An efficient Krylov-subspace-based algorithm to solve the dielectric-waveguide problem

An efficient scheme based on the bi-Lanczos algorithm has been developed for analysis of the dielectric-waveguide problem. A two-dimensional finite-difference scheme in the frequency domain is used to discretize the waveguide cross section. The resulting sparse eigenvalue problem is solved efficiently using the bi-Lanczos algorithm. Apart from solving the modes of the dielectric waveguide, a scheme to solve for the fields in the presence of a localized source is also described. Numerical results are also included to confirm the validity of the method.

A NOVEL TWO-DIMENSIONAL CONVECTION-DIFFUSION FINITE-DIFFERENCE SCHEME

In this article, we develop a two-dimensional finite-difference scheme for solving the convection-diffusion equation. The numerical method involves using transformation on the prototype scalar transport equation and transferring it to a Helmholtz equation. We apply the alternating-direction implicit scheme of Polezhaev to solve for the Helmholtz equation. As the key to success in simulating the convection-diffusion equation, we exploit the solution pertaining to the Helmholtz equation in the course of scheme development, thereby providing high-level accuracy to the prediction. Since this is a new method developed for solving the model equation, it is illuminating to conduct modified equation analysis on the discrete equation in order to make a full assessment of the proposed method. The results provide us with useful insights into the nature of the scheme. It is standard practice to validate the code by investigating test problems which are amenable to exact solutions to the working equation. Results show exact agreement for the one-dimensional test problem and good agreement with the analytic solutions for two-dimensional problems.

The delineation of drainage basins on the Greenland ice sheet for mass-balance analyses using a combined modelling and geographical information system approach

Balance fluxes have been calculated for the Greenland ice sheet using a two-dimensional finite difference scheme. They were calculated from a high-resolution digital elevation model (DEM) and mass-balance data for the Greenland ice sheet. The method applied assumes that there is a balance between the next accumulation and the ice flux at the margins and provides valuable information about the balance dynamics across the whole ice sheet. The balance fluxes, along with slope and slope direction information from the DEM were used in a geographical information system to delineate 10 major drainage basins. Each drainage basin can be considered as an independent flow unit and thus the individual mass balance of each basin can be calculated independently. Further, the fluxes across the grounding lines of several major outlet glaciers in northern Greenland have been calculated, resulting in a reassessment of the mass balance of northern Greenland. Copyright © 2000 John Wiley & Sons, Ltd.

Balance velocities and measured properties of the Antarctic ice sheet from a new compilation of gridded data for modelling

AbstractThis paper presents a new compilation of gridded datasets for three-dimensional modelling of the Antarctic ice sheet. These are for surface elevation, ice thickness, bedrock elevation and accumulation rate as interpolated on a 281 × 281 mesh with 20 km spacing, and encompass all the ice sheet and surrounding continental shelf. Data sources include the Bamber digital-elevation model from ERS-1 radar-altimeter data, a redigitization of available ice-thickness data, the Giovinetto accumulation data, recent ice-thickness data from British and German expeditions as well as accumulation data from German and Norwegian expeditions. In particular, new data were incorporated for the Filchner-Ronne Ice Shelf and for Dronning Maud Land, Antarctica, arising from the EPICA pre-site survey. Special attention was devoted to matching the various data sources carefully, both among themselves and across the grounding line and below the ice shelves, to enable ice-sheet expansion and retreat in dynamic situations. As an application, the balance flow is calculated over the entire ice sheet using a two-dimensional finite-difference scheme and compared with a previous assessment. This brought to light the existence of ice-streaming features extending well inland. A detailed zoom over Dronning Maud Land exhibits the general flow characteristics of interest for locating a future deep-drilling site. As a by-product, an updated value of 26.4 × 106km3 was obtained for the total volume of the ice sheet and ice shelves, or equivalent to 61.1 m of global sea-level rise after removal of the ice sheet and subsequent oceanic invasion and isostatic rebound. The total accumulation over the grounded ice sheet, including the Antarctic Peninsula, is 1924 Gta−1, or between 5 and 20% higher than earlier estimates. Including all the ice shelves, the value is 2344 Gt a−1.

A comparison of balance velocities, measured velocities and thermomechanically modelled velocities for the Greenland ice sheet

AbstractBalance velocities for the Greenland ice sheet have been calculated from a new digital-elevation model, accumulation-rates compilation and an existing ice-thickness grid, using a two-dimensional finite-difference scheme. The pattern of velocities over the ice sheet is presented and compared with velocities derived from synthetic-aperture-radar interferometry for part of northern Greenland and a limited number of global positioning system data. This comparison indicated that the balance-velocity scheme and boundary conditions used here provide a remarkably good representation of the dynamics of the ice sheet inland from the margins. It is suggested, therefore, that these balance-velocity data could provide a valuable method of constraining a numerical ice-sheet model. The balance velocities were compared with the diagnostic velocity field calculated from several different configurations of a numerical ice-sheet model. The general pattern of flow agrees well. The detail, however, is quite different. For example, the large (>300km) ice stream in the northeast is not generated by the numerical model and much of the detailed flow pattern is completely lost due to the limited model resolution and limitations in the model physics.

An analysis of balance velocities over the Greenland ice sheet and comparison with synthetic aperture radar interferometry

AbstractBalance velocities for the Greenland ice sheet have been calculated from a new digital elevation model (DEM), accumulation rates and an existing ice-thickness grid, using a fully two-dimensional finite-difference scheme. The pattern of velocities is compared with velocities derived from synthetic-aperture radar (SAR) interferometry for three different regions of the ice sheet. Differences between the two estimates of velocity highlight the respective strengths and weaknesses of the datasets and techniques used. A comparison with ten global positioning system-derived velocities indicates that the balance-velocity scheme and input datasets used here provide a remarkably good representation of the velocity distribution inland from the margins. These balance-velocity data, therefore, could help constrain numerical ice-sheet models. The balance velocities were found to be unreliable close to the ice-sheet margins due to larger errors in ice thickness, surface slope and ablation rate in this region. Comparison of the balance velocities with SAR interferometry in the region of the “Northeast Greenland Ice Stream” indicates the importance of the smoothing distance that must be applied to the DEM before calculating balance velocities. A smoothing distance of 20 times the ice thickness gave good agreement between the two measures of velocity.

Hybrid three-dimensional finite-difference and finite-element analysis of seismic wave induced fluid–structure interaction of a vertical cylinder

The two-dimensional finite-difference scheme has been extended to three dimensions to solve nonlinear hydrodynamic pressures and structural responses of a deformable, vertical and circular surface-piercing offshore cylinder during earthquakes. A complete three-dimensional analysis has been made with both the three-dimensional equations of motion and the simultaneous action of three components of ground acceleration included in the analysis. Not only the magnitude but also the direction of the acting ground motion can be varied with time. The dynamic response of a cylinder is approximated by the displacements in the fundamental modes of vibration. A comparison of the dynamic displacement of the cylinder with and without surrounding sea water has been made. The flexibility of the offshore cylinder can significantly increase the hydrodynamic pressures acting on cylinder faces, that is, the fluid–structure interaction is necessary in offshore cylinder analysis. Although the hydrodynamic pressure induced by the vertical ground acceleration of the El Centro 1979 earthquake is significant, the calculated structural dynamic response of a cylinder is very small and the corresponding resultant hydrodynamic force is almost nil. The hydrodynamic force induced by two-horizontal ground acceleration is about the same as that by three simultaneous components of ground acceleration. For a solid and stubbier circular cylinder, the vertical component of ground acceleration may be neglected.

Numerical Study of the Mechanism of Wavy Interface Generation in Explosive Welding.

Wavy interface generation is an interesting characteristic of explosive welding. There have been many theoretical and experimental discussions of the mechanism of wavy interface generation, but there is no consensus of opinion as yet. We report on a numerical approach to the analysis of the mechanism of wavy interface generation in explosive welding. The wavy interface generation process is calculated using a two-dimensional finite difference scheme for elastic-plastic materials. In this simulation the equations for a symmetric collision between copper plates are solved. The concept behind this study is that the wavy interface and vortex street are caused by the velocity distribution of shear flow and periodic disturbance at the interface. The calculation results show the deformation process of the wavy interface and we obtain qualitative agreements between numerical and experimental results.

Numerical Study of the Mechanism of Wavy Interface Generation in Explosive Welding.

Wavy interface generation is an interesting characteristic of explosive welding. There have been many theoretical and experimental discussions of the mechanism of wavy interface generation, but there is no consensus of opinion as yet. We report a numerical approach to the analysis of the mechanism of interface wave generation in explosive welding. The phenomena of the interface wave generation are calculated using a two-dimensional finite difference scheme for elastic-plastic materials. In this simulation the equations for a symmetric collision between copper plates are solved. The concept behind this study is that the wavy interface and vortex streets are caused by the velocity distribution of shear flow and periodic disturbances at the interface. The calculation results show the deformation process of the wavy interface and we obtain qualitative agreements between numerical and experimental results.

Stress wave fracturing of biological stones with lithotripter pulses

Stress waves contributions to stone fragmentation during lithotripsy are investigated both theoretically and experimentally. A two-dimensional finite difference scheme is developed to analyze the time evolution of the strain fields inside irregularly shaped solids subjected to ultrasonic pulses that simulate lithotripter shock waves. The reflections and superposition of stress waves inside the stones are analyzed to better understand the effects of stone parameters and geometry on the induced internal strains and fragmentation during lithotripsy. Numerical results show the focusing effect of the concave backsurface of a spherical stone, with the subsequent formation of focal zones (caustics). The focusing is reduced when a section of the back surface of the stone is removed. Principal strain contours depict the time evolution of the stress waves as they refract and reflect at the stone boundaries. Locations of maximum stresses are calculated and compared to locations of crack initiation in experiments with stones of similar geometry. The calculated time evolution of strain at fixed points within a stone is compared to imbedded silicon strain gauge measurements. Fracture characteristics of synthetic stones show internal crack initiation and subsequent propagation to external stone surfaces, indicative of internal stress fragmentation mechanisms. [Work supported by NIH.]

Vectorial simulation of passive TE/TM mode converter devices on InP

Accurate calculations of polarization rotation at abrupt waveguide interfaces are performed by a vectorial two-dimensional finite difference scheme. It is shown that the small orthogonal component of the vector field in rib waveguides can play an important role in mode polarization rotation. Polarization rotations in mode converters fabricated on InGaAsP/InP by periodically arranged tilted waveguides or periodic loading of the rib are explained by these calculations, and design rules for passive mode converter devices are given. >

Visualization of natural convection in a horizontal cylindrical annulus by the liquid crystal particle suspension method.

The temperature distribution and flow pattern of natural convection in a horizontal cylindrical annulus were made visible using the liquid crystal particle suspension method. Results were obtained using three kinds of silicone oil having different high Prandtl numbers. The Rayleigh number based on the gap width varied from 6.1×104 to 3.8×105. Resulting flows were observed to be steady in all cases. Numerical investigation was also performed with a transient two-dimensional finite difference scheme which have been reported in the previous paper. Comparison between the experimental and numerical results shows good agreement.