Robust Numerical Integration Research Articles

Numerical quadrature for the Prandtl—Meyer function at high temperature with application for air

When the stagnation temperature of the combustion chamber or ambient air increases, the specific heats and their ratio do not remain constant any more, and start to vary with this temperature. The gas remains perfect, except, it will be calorically imperfect and thermally perfect. A new generalized form of the Prandtl Meyer function is developed, by adding the effect of variation of this temperature, lower than the threshold of dissociation. The new relation is presented in the form of integral of a complex analytical function, having an infinite derivative at the critical temperature. A robust numerical integration quadrature is presented in this context. The classical form of the Prandtl Meyer function of a perfect gas becomes a particular case of the developed form. The comparison is made with the perfect gas model for aim to present a limit of its application. The application is for air.

The Ndynamics package—Numerical analysis of dynamical systems and the fractal dimension of boundaries

The Ndynamics package—Numerical analysis of dynamical systems and the fractal dimension of boundaries

A Reaction-Based River/Stream Water Quality Model: Reaction Network Decomposition and Model Application

This paper describes details of an automatic matrix decomposition approach for a reaction-based stream water quality model. The method yields a set of equilibrium equations, a set of kinetic-variable transport equations involving kinetic reactions only, and a set of component transport equations involving no reactions. Partial decomposition of the system of water quality constituent transport equations is performed via Gauss-Jordan column reduction of the reaction network by pivoting on equilibrium reactions to decouple equilibrium and kinetic reactions. This approach minimizes the number of partial differential advective-dispersive transport equations and enables robust numerical integration. Complete matrix decomposition by further pivoting on linearly independent kinetic reactions allows some rate equations to be formulated individually and explicitly enforces conservation of component species when component transport equations are solved. The methodology is demonstrated for a case study involving eutrophication reactions in the Des Moines River in Iowa, USA and for two hypothetical examples to illustrate the ability of the model to simulate sediment and chemical transport with both mobile and immobile water phases and with complex reaction networks involving both kinetic and equilibrium reactions.

Special Issue on Computational Multibody Dynamics

This issue of the Journal of Computational and Nonlinear Dynamics presents a snapshot of the state-of-the-art in Computational Multibody Dynamics (CMD). Since any attempt at approximating the time evolution of dynamical systems begins with the question of finding robust and efficient numerical integration methods, this special issue opens with several contributions that concentrate on this topic. Specifically, the first four papers cover both analytical and practical aspects related to numerical time integration, for smooth and non-smooth dynamics, for rigid and flexible body systems. Flexible body dynamics is further discussed in two contributions. The first outlines a procedure to reduce the computational burden that flexibility typically places on the numerical solution; the second illustrates how recent CMD advances are leveraged in the context of a challenging biomechanics application. One of the recurrent themes of the research effort within the CMD community pertains to the focus on system-level engineering applications. This usually requires leaving behind the convenience of dealing with simple systems, concentrating instead on complex problems such as active rotorcraft control, quantifying the uncertainty in the aerodynamic loads acting on ground vehicles, or designing an optimal set of experiments for parameter estimation for accurate vehicle dynamics simulation.Moving beyond the breath aspect discussed above, an attempt was made to illustrate the fact that the CMD community is robust and dynamic. Several contributions concentrate on topics in well established areas that provide the core competencies of this community: numerical integration, flexible body dynamics, nonlinear control, and frictional contact dynamics. Other contributions concentrate on application fields that, until recently, fell outside the scope of this community either because of a lack of dialogue with colleagues from other disciplines, or because the available analytical and computational tools were not mature enough. Uncertainty quantification, system decomposition for distributed computing, biomechanics, and molecular dynamics provide new areas of engineering and science where the CMD community can make fruitful contributions, as illustrated by several papers of this issue.Today, Computational Multibody Dynamics addresses real world, system-level applications leading to better engineering designs within a shorter development cycle and at a reduced cost. Looking forward, there are many reasons to be optimistic. Members of the community are making steady analytical advances that leverage ever increasing computational power. The emergence of commodity high performance computing, such as GPU computing, for instance, will surely be met in our community by an enthusiastic group of individuals ready to push the boundaries of Computational Multibody Dynamics to understand through simulation the dynamics of complex systems be it at the atomic (molecular dynamics), meso (granular dynamics), or macro (vehicle dynamics) scales.

Numerical Quadrature for the Prandtl Meyer Function at High Temperature with Application for Air

When the stagnation temperature of the combustion chamber or ambient air increases, the specific heats and their ratio do not remain constant any more, and start to vary with this temperature. The gas remains perfect, except, it will be calorically imperfect and thermally perfect. A new generalized form of the Prandtl Meyer function is developed, by adding the effect of variation of this temperature, lower than the threshold of dissociation. The new relation is presented in the form of integral of a complex analytical function, having an infinite derivative at the critical temperature. A robust numerical integration quadrature is presented in this context. The classical form of the Prandtl Meyer function of a perfect gas becomes a particular case of the developed form. The comparison is made with the perfect gas model for aim to present a limit of its application. The application is for air.

Accurate solution of multi-region continuum biomolecule electrostatic problems using the linearized Poisson-Boltzmann equation with curved boundary elements.

We present a boundary-element method (BEM) implementation for accurately solving problems in biomolecular electrostatics using the linearized Poisson-Boltzmann equation. Motivating this implementation is the desire to create a solver capable of precisely describing the geometries and topologies prevalent in continuum models of biological molecules. This implementation is enabled by the synthesis of four technologies developed or implemented specifically for this work. First, molecular and accessible surfaces used to describe dielectric and ion-exclusion boundaries were discretized with curved boundary elements that faithfully reproduce molecular geometries. Second, we avoided explicitly forming the dense BEM matrices and instead solved the linear systems with a preconditioned iterative method (GMRES), using a matrix compression algorithm (FFTSVD) to accelerate matrix-vector multiplication. Third, robust numerical integration methods were employed to accurately evaluate singular and near-singular integrals over the curved boundary elements. Fourth, we present a general boundary-integral approach capable of modeling an arbitrary number of embedded homogeneous dielectric regions with differing dielectric constants, possible salt treatment, and point charges. A comparison of the presented BEM implementation and standard finite-difference techniques demonstrates that for certain classes of electrostatic calculations, such as determining absolute electrostatic solvation and rigid-binding free energies, the improved convergence properties of the BEM approach can have a significant impact on computed energetics. We also demonstrate that the improved accuracy offered by the curved-element BEM is important when more sophisticated techniques, such as nonrigid-binding models, are used to compute the relative electrostatic effects of molecular modifications. In addition, we show that electrostatic calculations requiring multiple solves using the same molecular geometry, such as charge optimization or component analysis, can be computed to high accuracy using the presented BEM approach, in compute times comparable to traditional finite-difference methods.

A reaction-based river/stream water quality model: Model development and numerical schemes

Summary This paper presents the conceptual and mathematical development of a numerical model of sediment and reactive chemical transport in rivers and streams. The distribution of mobile suspended sediments and immobile bed sediments is controlled by hydrologic transport as well as erosion and deposition processes. The fate and transport of water quality constituents involving a variety of chemical and physical processes is mathematically described by a system of reaction equations for immobile constituents and advective–dispersive–reactive transport equations for mobile constituents. To circumvent stiffness associated with equilibrium reactions, matrix decomposition is performed via Gauss–Jordan column reduction. After matrix decomposition, the system of water quality constituent reactive transport equations is transformed into a set of thermodynamic equations representing equilibrium reactions and a set of transport equations involving no equilibrium reactions. The decoupling of equilibrium and kinetic reactions enables robust numerical integration of the partial differential equations (PDEs) for non-equilibrium-variables. Solving non-equilibrium-variable transport equations instead of individual water quality constituent transport equations also reduces the number of PDEs. A variety of numerical methods are investigated for solving the mixed differential and algebraic equations. Two verification examples are compared with analytical solutions to demonstrate the correctness of the code and to illustrate the importance of employing application-dependent numerical methods to solve specific problems.

On some topics for the numerical simulation of ductile fracture

In this work, we analyze some aspects of the macroscopic Gurson–Tvergaard–Needleman (GTN) constitutive model when it is addressed to solve ductile fracture problems by means of numerical simulations: (i) The analytical solutions of the material discontinuous bifurcation problem is performed. Closed and exact formulas are obtained. The so determined critical conditions and the developed strain localization mode are afterward studied and compared in crack growth problems. Even when this methodology of analysis is rather standard at the present, the conclusions drawn from this study differ significantly from that obtained with a similar analysis in quasi-brittle fracture problems. (ii) A new very robust numerical integration method for the GTN model, namely the Impl–Ex Method, is proposed. It is a low computational cost algorithm, equivalent to a linear problem per each integration step, with a reasonable precision for engineering purposes. Its accuracy and convergence rate is assessed by means of an error study applied to a ductile fracture test simulation. (iii) A detailed analysis of a plane strain ductile crack growth problem is performed in a material containing two size-scale of voids. In the analysis, particular attention is given to the mesh size dependence and to the coalescence of the larger void.

A reaction-based paradigm to model reactive chemical transport in groundwater with general kinetic and equilibrium reactions

This paper presents a reaction-based water quality transport model in subsurface flow systems. Transport of chemical species with a variety of chemical and physical processes is mathematically described by M partial differential equations (PDEs). Decomposition via Gauss–Jordan column reduction of the reaction network transforms M species reactive transport equations into two sets of equations: a set of thermodynamic equilibrium equations representing N E equilibrium reactions and a set of reactive transport equations of M– N E kinetic-variables involving no equilibrium reactions (a kinetic-variable is a linear combination of species). The elimination of equilibrium reactions from reactive transport equations allows robust and efficient numerical integration. The model solves the PDEs of kinetic-variables rather than individual chemical species, which reduces the number of reactive transport equations and simplifies the reaction terms in the equations. A variety of numerical methods are investigated for solving the coupled transport and reaction equations. Simulation comparisons with exact solutions were performed to verify numerical accuracy and assess the effectiveness of various numerical strategies to deal with different application circumstances. Two validation examples involving simulations of uranium transport in soil columns are presented to evaluate the ability of the model to simulate reactive transport with complex reaction networks involving both kinetic and equilibrium reactions.

A practical computational framework for the multidimensional moment-constrained maximum entropy principle

The maximum entropy principle is a versatile tool for evaluating smooth approximations of probability density functions with a least bias beyond given constraints. In particular, the moment-based constraints are often a common prior information about a statistical state in various areas of science, including that of a forecast ensemble or a climate in atmospheric science. With that in mind, here we present a unified computational framework for an arbitrary number of phase space dimensions and moment constraints for both Shannon and relative entropies, together with a practical usable convex optimization algorithm based on the Newton method with an additional preconditioning and robust numerical integration routine. This optimization algorithm has already been used in three studies of predictability, and so far was found to be capable of producing reliable results in one- and two-dimensional phase spaces with moment constraints of up to order 4. The current work extensively references those earlier studies as practical examples of the applicability of the algorithm developed below.

A simple robust numerical integration algorithm for a power‐law visco‐plastic model under both high and low rate‐sensitivity

AbstractThis note describes a simple and extremely robust algorithm for numerical integration of the power‐law‐type elasto‐viscoplastic constitutive model discussed by Perić (Int. J. Num. Meth. Eng. 1993; 36: 1365–1393). As the rate‐independent limit is approached with increasing exponents, the evolution equations of power‐law‐type models are known to become stiff. Under such conditions, the solution of the implicitly discretized viscoplastic evolution equation cannot be easily obtained by standard root‐finding algorithms. Here, a procedure which proves to be remarkably robust under stiff conditions is obtained by means of a simple logarithmic mapping of the basic backward Euler time‐discrete equation for the incremental plastic multiplier. The logarithm‐transformed equation is solved by the standard Newton–Raphson scheme combined with a simple bisection procedure which ensures that the iterative guesses for the equation unknown (the incremental equivalent plastic strain) remain within the domain where the transformed equation makes sense. The resulting implementation can handle small and large (up to order 106) power‐law exponents equally. This allows its effective use under any situation of practical interest, ranging from high rate‐sensitivity to virtually rate‐independent conditions. The robustness of the proposed scheme is demonstrated by numerical examples. Copyright © 2003 John Wiley & Sons, Ltd.

Robust numerical integration and pairwise independent random variables

A numerical integration method by means of random samples is called robust, if the variance of its error is as small as that of the i.i.d.-sampling for any integrand. To reduce the randomness of robust numerical integration, we use pairwise random samples instead of i.i.d. samples. Among others, we recommend the discrete random Weyl sampling for the quick generation of pairwise independent samples.

Implicit numerical integration for a kinematic hardening soil plasticity model

AbstractSoil models based on kinematic hardening together with elements of bounding surface plasticity, provide a means of introducing some memory of recent history and stiffness variation in the predicted response of soils. Such models provide an improvement on simple elasto‐plastic models in describing soil behaviour under non‐monotonic loading. Routine use of such models requires robust numerical integration schemes. Explicit integration of highly non‐linear models requires extremely small steps in order to guarantee convergence. Here, a fully implicit scheme is presented for a simple kinematic hardening extension of the Cam clay soil model. The algorithm is based on the operator split methodology and the implicit Euler backward integration scheme is proposed to integrate the rate form of the constitutive relations. This algorithm maintains a quadratic rate of asymptotic convergence when used with a Newton–Raphson iterative procedure. Various strain‐driven axisymmetric triaxial paths are simulated in order to demonstrate the efficiency and good performance of the proposed algorithm. Copyright © 2001 John Wiley & Sons, Ltd.

Numerical problems in the solution of oxidation and combustion models

Faced with large and detailed kinetic schemes of pyrolysis, partialoxidation and combustion of hydrocarbon mixtures, the complex processmodels usually constitute stiff systems of both differential and coupleddifferential-algebraic equations. This paper discusses different typicalapplications of a new robust and efficient numerical integration methodto solve these problems. The examples refer to gas solid catalyticreactors, partial oxidation of hydrocarbons and homogeneous andheterogeneous combustion processes. The results show the clear advantages of this numerical package whencompared with those commonly referred to and used in the scientificliterature. In our experience, the robustness of the new method allows theuser to overcome the failures of other codes and to reduce the criticalselection of relative and absolute tolerances. Moreover, the efficiency ofthe method, in terms of computing time, is significantly improved mainlydue to a better policy of Jacobian evaluations.

MASTER: macroscopic traffic simulation based on a gas-kinetic, non-local traffic model

We present a gas-kinetic (Boltzmann-like) traffic equation that is not only suited for low vehicle densities, but also for the high-density regime, as it takes into account the forwardly directed interactions, effects of vehicular space requirements like increased interaction rates, and effects of velocity correlations that reflect the bunching of cars, at least partially. From this gas-kinetic equation, we systematically derive the related macroscopic traffic equations. The corresponding partial differential equations for the vehicle density and average velocity are directly related to the quantities characterizing individual driver–vehicle behavior, and, as we show by calibration of the model, their optimal values have the expected order of magnitude. Therefore, the model allows to investigate the influences of varying street and weather conditions or freeway control measures. We point out that, because of the forwardly directed interactions, the macroscopic equations contain non-local instead of diffusion or viscosity terms. This resolves some of the inconsistencies found in previous models and allows for a fast and robust numerical integration, so that several thousand freeway kilometers can be simulated in real-time. It turns out that the model is in good agreement with the experimentally observed properties of freeway traffic flow. In particular, it reproduces the characteristic outflow and dissolution velocity of traffic jams, as well as the phase transition to “synchronized” congested traffic. We also reproduce the five different kinds of congested states that have been found close to on-ramps (or bottlenecks) and present a “phase diagram” of the different traffic states in dependence of the main flow and the ramp flow, showing that congested states are often induced by perturbations in the traffic flow. Finally, we introduce generalized macroscopic equations for multi-lane and multi-userclass traffic. With these, we investigate the differences between multi-lane simulations and simulations of the effective one-lane model.

Random Weyl sampling for robust numerical integration of complicated functions

Random Weyl sampling for robust numerical integration of complicated functions

Mechanical and numerical modeling of a porous elastic–viscoplastic material with tensile failure

The objective of this paper is to develop simple but comprehensive constitutive equations that model a number of physical phenomena exhibited by dry porous geological materials and metals. For geological materials the equations model: porous compaction; porous dilation due to distortional deformation and tensile failure; shear enhanced compaction; pressure hardening of the yield strength; damage of the yield strength due to distortional deformation and porosity changes; and dependence of the yield strength on the Lode angle. For metals the equations model: hardening of the yield strength due to plastic deformation; pressure and temperature dependence of the yield strength, and damage due to nucleation of porosity during tensile failure. The equations are valid for large deformations and the elastic response is hyperelastic in the sense that the stress is related to a derivative of the Helmholtz free energy. Also, the equations are viscoplastic with rate dependence occurring in both the evolution equations of porosity and elastic distortional deformations. Moreover, formulas are presented for robust numerical integration of the evolution equations at the element level that can be easily implemented into standard computer programs for dynamic response of materials.

Derivation, properties, and simulation of a gas-kinetic-based, nonlocal traffic model

We derive macroscopic traffic equations from specific gas-kinetic equations, dropping some of the assumptions and approximations made in previous papers. The resulting partial differential equations for the vehicle density and average velocity contain a non-local interaction term which is very favorable for a fast and robust numerical integration, so that several thousand freeway kilometers can be simulated in real-time. The model parameters can be easily calibrated by means of empirical data. They are directly related to the quantities characterizing individual driver-vehicle behavior, and their optimal values have the expected order of magnitude. Therefore, they allow to investigate the influences of varying street and weather conditions or freeway control measures. Simulation results for realistic model parameters are in good agreement with the diverse non-linear dynamical phenomena observed in freeway traffic.

Estimating imprecision profiles in biochemical analysis

We describe a computer program IMPROFIL which determines an imprecision profile of an analytical method from replicated measurements of samples. It calculates the variance function, the coefficient of variation, the power of definition, the critical limit, the limit of detection and the lower limit of quantification. The primary property, the variance function, is determined by two alternative methods: the conventional maximum approximate conditional likelihood method and the newly developed weighted absolute deviation method. For all quantities, confidence intervals are obtained using the bootstrap procedure. The program combines the use of robust numerical techniques, user-friendliness and integration into a spreadsheet program for data pre- and post-processing. The algorithms used are described in detail. Tests with synthetic data sets are used to validate the method and to establish its powers and limitations. Finally, its application to a practical analytical task (tumor marker CA 15-3 in human sera) is reported. For the method to yield a reliable estimate of the variance function and the derived properties, certain minimum requirements on the raw data must be met: They have to be spread throughout the concentration range of interest, there should not be less than three replicates per specimen, and there must be at least of the order of 25 (better at 50) specimens.

Stability Boundary for Pseudo-Random Parametric Excitation of a Linear Oscillator

Stability bounds are outlined for the null solution of the equation describing the response of a linear damped oscillator excited through periodic coefficients, the excitation being a form of Rice noise comprising equal amplitude sinusoids with frequencies at equal intervals in the vicinity of twice the natural frequency of the system, but with pseudo-random initial phases. Stability was investigated by the monodromy matrix method, which is exact apart from errors due to numerical integration, and by the approximate method due to R. A. Struble, which replaces the dependent variable by its amplitude and a phase variable. Struble’s method gives the main features of the stability diagram and leads to faster and more robust numerical integration with potential advantages for nonlinear and several degree-of-freedom systems, but loses much of the detail. When the frequency spacing is relatively large, the stability diagram is closely related to that for Mathieu’s equation, but the detailed shape becomes very complicated as the frequency spacing decreases. Quantitative comparison with the corresponding boundary for Gaussian white noise excitation shows very approximate equivalence.