Time Step Of Numerical Integration Research Articles

An operator-splitting numerical scheme for relativistic magnetohydrodynamics

We describe a novel operator-splitting approach to numerical relativistic magnetohydrodynamics designed to expand its applicability to the domain of ultra-high magnetisation. In this approach, the electromagnetic field is split into the force-free component, governed by the equations of force-free degenerate electrodynamics (FFDE), and the perturbation component, governed by the perturbation equations derived from the full system of relativistic magnetohydrodynamics (RMHD). The combined system of the FFDE and perturbation equations is integrated simultaneously, for which various numerical techniques developed for hyperbolic conservation laws can be used. At the end of every time-step of numerical integration, the force-free and the perturbation components of the electromagnetic field are recombined and the result is regarded as the initial value of the force-free component for the next time-step, whereas the initial value of the perturbation component is set to zero. To explore the potential of this approach, we build a 3rd-order WENO code, which was used to carry out 1D and 2D test simulations. Their results show that this operator-splitting approach allows us to bypass the stiffness of RMHD in the ultra-high-magnetisation regime where the perturbation component becomes very small. At the same time, the code performs very well in problems with moderate and low magnetisation too.

Real-time hybrid simulation using analogue electronic computer technology

In the field of structural dynamics, real-time hybrid simulation (RTHS) continues to receive increased interest from researchers conducting component testing incorporating system-level behaviour. Increased digital computing power has allowed advances in RTHS. However, limitations in RTHS will persist due to the effects of discrete errors caused by quantisation and corresponding numerical integration time step limitations, which can limit the bandwidth and nonlinearities of the system studied. In this paper, an analogue electronic computer is constructed to conduct RTHS of a modelled base-isolated structure with a physically tested viscous damping device. The analogue computer models a two-degree-of-freedom (2DOF) structure and solves the equations of motion required in RTHS testing. Results of the RTHS tests using the analogue computer are compared to RTHS tests implemented with a digital computer in order to validate the proposed analogue RTHS method. Further applications of RTHS testing with analogue computing including high-frequency dynamic testing are discussed.

Real-time hybrid simulation using analogue electronic computer technology

In the field of structural dynamics, real-time hybrid simulation (RTHS) continues to receive increased interest from researchers conducting component testing incorporating system-level behaviour. Increased digital computing power has allowed advances in RTHS. However, limitations in RTHS will persist due to the effects of discrete errors caused by quantisation and corresponding numerical integration time step limitations, which can limit the bandwidth and nonlinearities of the system studied. In this paper, an analogue electronic computer is constructed to conduct RTHS of a modelled base-isolated structure with a physically tested viscous damping device. The analogue computer models a two-degree-of-freedom (2DOF) structure and solves the equations of motion required in RTHS testing. Results of the RTHS tests using the analogue computer are compared to RTHS tests implemented with a digital computer in order to validate the proposed analogue RTHS method. Further applications of RTHS testing with analogue computing including high-frequency dynamic testing are discussed.

Study on wheel polygonization of a metro vehicle based on polygonal wear simulation

A long-term iterative wear model is proposed to predict the development of polygonal wear of railway wheels. To improve the simulation accuracy, the time step for numerical integration, wheel profile updating strategy and curve fitting methods are analysed. Then a metro vehicle with 8th to 9th order polygonal wheel is presented as a case study where a vehicle resonance speed leading to a peculiar growth of wheel polygonization is identified by simulation and validated through field tests. Based on the simulated polygonal wear evolutions and analyses on vehicle vibrations, the mechanism of wheel polygonization is explained. Finally, two new strategies to investigate wheel polygonization are proposed.

From multiscale modeling to design of synchronization mechanisms in mesh antennas

Flexible multibody simulation is supposed to be a promising method for studying the deployment dynamics of large mesh antennas and assisting the design of the related crucial mechanisms. However, fulfilling a full-scale deployment simulation of a mesh antenna is so challenging that it can be rarely found in the literature. In our opinion, the bottleneck of modeling some mesh antennas is dealing with the multiscale physics mainly sourced from the cable-pulley systems that are frequently used to transfer forces for long distance range. A traditional model of a cable-pulley is based on contact method, which involves fine cable mesh, a large amount of contact detection, and small time step for numerical integration, resulting in low calculation efficiency. To cope with this difficulty, we proposed a novel and efficient multiscale method to handle the cable-pulley systems. The basic idea is avoiding the contact detection by separating a cable into two segments: non-contact and contact segments, according to which part of the cable is in contact with the pulley. The non-contact segment is meshed with variable-length elements based on the arbitrary Lagrangian–Eulerian (ALE) formulation. The contact segment is considered an invisible segment to constrain the variable-length cable element on the pulley; its border is dynamically located by the relative configuration between the pulley and the cable. The accuracy and efficiency of the proposed method were confirmed by a numerical example. Then, it was used to evaluate, compare and design the synchronization mechanism in a kind of mesh antenna. The obtained results can directly guide the actual design of mesh antennas, and the proposed multiscale method can be applied to simulating other similar mechanical systems.

Автоматизований розрахунок запасу статичної стійкості у перетині електроенергетичної системи

The approach adopted in the work involves focusing on the control of static stability of the UES on the criterion of existence of the regime, subject to weighing of certain sections. In practice, this means that on the trajectory of weighing it is necessary to find the mode which has the greatest power in the investigated section, for which the coefficient of the reserve of static stability at the active power in the intersection is calculated and normalized by the "Guideline" on the stability of the power systems. A modification of the program of calculating the stationary regime in a linearized formulation, which allows us to perform weighting at the angle of the vectors of electromotive forces (EMF), without the iterative approach to the boundary regime and to specify the stock of static stability at the intersection of the Electric Power System (EPS). The linearization of the EPS model is carried out by replacing the loads in the nodes of the conductivities, and the generating nodes — by constant electromotive forces, which are fixed on the reactive supports, the magnitude of which corresponds to the staticity of the voltage of the automatic regulators of excitation of the synchronous generators. Such a model is used in the programs for calculating the dynamic stability for the simulation of EPS at the time step of numerical integration of the differential equations of dynamics of the system. A block diagram of the algorithm for determining the maximum and unacceptable fluxes in controlled sections of the electric power system is proposed. A separate module for the automated calculation of the static stability reserve at the intersection of the EPS with a choice of certain weighting trajectories from the library was created. Trajectories weighing are prepared in a specially created editor, which provides a step-by-step change in the power of power stations and electrical loads in different parts of the mathematical model of the power system and the change of the angle of the EMF generators. Calculations were made to determine the boundary power at the cross sections of the UES of Ukraine for a scheme containing 614 nodes and 975 branches.

A combined approach for modeling particle behavior in granular impact damper using discrete element method and cellular automata

A particle impact damper is a vibration absorber type that consists of a container attached to a primary vibrating structure. It also contains many particles that are constrained to move inside the container, whereby the damping effect can be obtained by collision between particles and the container. The discrete element method (DEM) has been developed for modeling granular systems, where the kinematics of each particle are calculated numerically using the equations of motion. However, the computational time is significant since the algorithm checks for particle contacts for all possible particle combinations. The use of a cellular automata (CA) modeling technique may provide increased computational efficiency due to the local updating of variables, and the discrete treatment of time and space. In this study, we propose a new approach combining DEM with CA for modeling a granular damper under a forced excitation. We use DEM to describe the particle motion according to the equations of motion, while CA is introduced for the particle contact checks in discrete space. We also investigate the effect of simplification in the contact force model, which allows the unit time step of numerical integration to become larger than that used in the strict model. It is shown that the suggested particle contact scanning method and the force approximation model contribute to the reduction of the computational time, and neither degenerates the calculation accuracy nor causes the numerical instability.

Three-Dimensional Dynamic Model of TEHD Tilting-Pad Journal Bearing—Part II: Parametric Studies

This paper presents a new analysis method for a thermo-elasto-hydro-dynamic (TEHD) tilting pad journal bearing (TPJB) system to reach a static equilibrium condition adopting nonlinear transient dynamic solver, whereas earlier studies have used iteration schemes such as Newton–Raphson method. The theoretical TPJB model discussed in Part I of this research is combined into a newly developed algorithm to perform a bearing dynamic analysis and present dynamic coefficients. In the nonlinear transient dynamic solver, physical and modal coordinates coexist for computational efficiency, and transformation between modal and physical coordinate is performed at each numerical integration time step. Variable time step Runge–Kutta numerical integration scheme is adopted for a reliable and fast calculation. Nonlinear time transient dynamic analysis and steady thermal analysis are combined to find the static equilibrium condition of the TPJB system, where the singular matrix issue of flexible pad finite element (FE) model is resolved. The flexible pad TPJB model was verified by comparison with other numerical results. Simulation results corresponding with the theoretical model explained in Part I are presented and discussed. It explains how the TPJB dynamic behavior is influenced by a number of eigenvector of flexible pad FE model and pad thickness. Preload change under fluid and thermal load is examined.

A ‘frozen electric-field’ approach to simulate repetitively pulsed nanosecond plasma discharges and ignition of hydrogen–air mixtures

High-fidelity modelling of nanosecond repetitively pulsed discharges (NRPDs) is burdened by the multiple time and length scales and large chemistry mechanisms involved, which prohibit detailed analyses and parametric studies. In the present work, we propose a ‘frozen electric-field’ modelling approach to expedite the NRPD simulations without adverse effects on the solution accuracy. First, a burst of nanosecond voltage pulses is simulated self-consistently until the discharge reaches a stationary state. The calculated spatial distributions and temporal evolution of the electric field, electron density and electron energy during the last pulse are then stored in a library and the electrical characteristics of subsequent pulses are frozen at these values. This strategy allows the timestep for numerical integration to be increased by four orders of magnitude (from 10−13 to 10−9 s), thereby significantly improving the computational efficiency of the process. Reduced calculations of a burst of 50 discharge pulses show good agreement with the predictions from a complete plasma model (electrical characteristics calculated during each pulse). The error in species densities is less than 20% at the centre of the discharge volume and about 30% near the boundaries. The deviations in temperature, however, are much lower, at 5% in the entire domain. The model predictions are in excellent agreement with measured ignition delay times and temperatures in H2–air mixtures subject to dielectric barrier NRPD over a pressure range of 54–144 Torr with equivalence ratios of 0.7–1.2. The OH density increases with pressure and triggers low-temperature fuel oxidation, which leads to rapid temperature rise and ignition. The ignition delay decreases by a factor of 2, with an increase in pressure from 54 to 144 Torr. In contrast, an increase in the H2–air equivalence ratio from 0.7 to 1.2 marginally decreases the ignition delay by about 20%. This behaviour is attributed to the insensitivity of OH production rates to the variation in the equivalence ratio.

An efficient approximate method for solving the problem of the establishment of chemical equilibrium in the products of explosion of gas mixtures

A method for calculating the changing composition of the explosion products in the case where the chemical equilibrium is absent but the bimolecular reactions are in quasi-equilibrium is developed. At each time step of numerical integration, the change in the total number of molecules in the system is calculated using a single differential equation written based on the selected kinetic mechanism. Then, the mixture composition is calculated under the assumption that the entropy of the mixture reaches its maximum value at a given internal energy, mass density, and molar mass. It is shown that the proposed method can be used to calculate the characteristics of explosive transformation processes after the induction period.

Equivalent circuits for implicit Runge-Kutta methods in circuit simulators for nonlinear circuits

This paper describes a method of implementing RadauIIA and LobattoIIIA of implicit Runge-Kutta formulas into circuit simulators for nonlinear circuits as numerical integration. These implicit Runge-Kutta methods have high orders and are A-stable. Equivalent circuits at discrete time for linear and nonlinear elements are proposed. Circuits at times between past and present time are needed in addition to the equivalent circuit at present time. Solutions at intermediate and present times must be estimated simultaneously. So, the size of equivalent circuit becomes larger than the numerical integration in conventional circuit simulators. However, since the orders of these algorithms are high, this problem is solved by using larger time step for numerical integration compared to conventional methods to save calculation time. The implicit Runge-Kutta and conventional methods are compared in terms of accuracy and computational costs using example circuits.

On the FEM modeling of mechanical systems controlled by relative motion of a member: A pendulum–mass interaction test case

A new approach for simulating the dynamics of mechanical systems with relative movements of members by the finite element method (FEM) is proposed. The moving member is modeled by a mass traveling along a guiding beam. The mass is incorporated into the beam element that is being traversed. The nodal inertia forces due to the relative mass’ motion are derived from the rate of changing the mass matrix of this composite element. These inertia forces, which are added to the regular nodal forces, reflect the Coriolis-type effects present in the system. Several ways of calculating and including them into the time-integrating procedure are presented. The approach allows treating the moving member as a controller; with a desired pattern of the member's motion imposed on the system as input to analyze its response as output. The proposed procedure is tested on the problem of controlling pendulum-like oscillations of a beam by a moving mass. The element size and the time-step of numerical integration to secure sufficient accuracy of the approach for different mass motion patterns are discussed in detail. In particular, it is demonstrated that the case of the mass motion pattern attenuating the pendulum oscillations with a constant damping ratio can be precisely simulated.

Simulation of Two-Phase Flow in Complex Systems

The paper presents integral methods for simulating two-phase flow transients in complex cooling systems, such as those in nuclear power plants. The methods are designed to simplify presently prevailing thermohydraulics simulation methods without a loss in simulation fidelity. The paper describes the inherent, but unnecessary, complexity of currently used simulation models, explains their inherent shortcomings, and foretells the impact of current code development trends on future capabilities to resolve safety issues in light of growing code complexities and inflexibility. The purpose of the paper is to present simpler alternatives.Integral methods described in the paper facilitate flexibility via computer-automated modularity and simplicity. They provide transparency through analytical methods. Integral methods replace partial by ordinary differential equations and thereby simplify the mathematical model formulation and achieve numerical integration with minimal numerical damping. The models connect important physical characteristic response times with the time step for numerical integration. The mixture model of nonhomogeneous, nonequilibrium two-phase flow, the integral of the volumetric flux divergence equation, and the integrals of the system of coupled loop momentum balances for interconnected loops in complex thermohydraulic systems play central modeling roles. A new and compact formulation of these equations facilitates a computerized system assembly of component models, each one being judiciously selected from a model library to impose the minimum necessary complexity. The method of assembly is based on linear algebra and accommodates any combination of phasic flow directions anywhere in the hydraulic system and at any time.

Scaling law for the transient behavior of type-II neuron models

We study the transient regime of type-II biophysical neuron models and determine the scaling behavior of relaxation times tau near but below the repetitive firing critical current, tau approximately or equal to C(I(c)-I)(-Delta). For both the Hodgkin-Huxley and Morris-Lecar models we find that the critical exponent is independent of the numerical integration time step and that both systems belong to the same universality class, with Delta=1/2. For appropriately chosen parameters, the FitzHugh-Nagumo model presents the same generic transient behavior, but the critical region is significantly smaller. We propose an experiment that may reveal nontrivial critical exponents in the squid axon.

Sensitivity testing for improved efficiency of a statistical oil-spill risk analysis model

The US Minerals Management Service performs oil-spill risk analyses using, in part, a statistical model of hypothetical oil-spill trajectories. The “Oil-spill Risk Analysis” (OSRA) model is driven by analyzed sea surface winds and model-generated ocean surface currents. The OSRA model calculates thousands of oil-spill trajectories over extended areas of US waters and tabulates the frequencies with which the simulated oil spills contact the geographic boundaries of designated natural resources within a specified number of days after the simulated spill events are initiated. The limits of computer mass storage and speed impose constraints on the number of oil-spill trajectory simulations and the time step of the numerical integration generating the spill trajectories. Also, the model winds and ocean currents have limited spatial and temporal resolution, producing sensible limits on the spatial density of the simulated oil spills and on the numerical integration time step. This investigation attempted to determine the smallest practical spatial interval between simulated oil spill release sites and the shortest integration time step beyond which the frequencies of oil-spill contact essentially do not change. In addition, we investigated the number of hypothetical oil spills initiated per day from a single location and two methods of numerical integration to generate the spill trajectories—simple forward time stepping and the fourth order Runge–Kutta method. Applying the OSRA model to the Gulf of Mexico using 9 years of winds and ocean current data, we found a smallest distance between simulated spills of 6 nautical miles and a shortest integration time step of 1 h using the fourth order Runge–Kutta integration method. A 20-min time step would be needed with simple forward time stepping. Initiating hypothetical spills eight times per day produced essentially the same contact probabilities as initiating spills once per day. These intervals were practical for OSRA model runs over the US portion of the Gulf of Mexico using commonly available desktop personal computers.

High-accuracy statistical simulation of planetary accretion: I. Test of the accuracy by comparison with the solution to the stochastic coagulation equation

The object of this series of studies is to develop a highly accurate statistical code for describing the planetary accumulation process. In the present paper, as a first step, we check the validity of the method proposed by Wetherill and Stewart (1989) by comparing the results obtained by their method with the analytical solution to the stochastic coagulation equation (or to a well-evaluated numerical solution). As the collisional probability Ai j between bodies with masses of im1 and jm1 (m1 being the unit mass), we consider the two cases: one is Ai j∝ i × j and another is Ai j ∝ min(i, j )(i 1/3 + j1/3)(i + j ). In both cases, it is known that runaway growth occurs. The latter case corresponds to a simplified model of the planetesimal accumulation. We assumed that a collision of two bodies leads to their coalescence. Wetherill and Stewart’s method contains some parameters controlling the practical numerical computation. Among these, two parameters are important: the mass division parameter δ, which determines the mass ratio of the adjacent mass batches, and the time division parameter ∈, which controls the size of a time step in numerical integration. Through a number of numerical simulations for the case of Ai j = i × j, we find that when δ ≤ 1.6 and ∈ ≤ 0.03 the numerical simulation can reproduce the analytical solution within a certain level of accuracy independently of the size of the body system. For the case of the planetesimal accumulation, it is shown that the simulation with δ ≤ 1.3 and ∈ ≤ 0.04 can describe precisely runaway growth. Because the accumulation process is stochastic, in order to obtain reliable mean values it is necessary to take the ensemble mean of the numerical results obtained with different random number generators. It is also found that the number of simulations, Nc, demanded to obtain the reliable mean value is about 500 and does not strongly depend on the functional form of Ai j. From the viewpoint of the numerical handling, the above value of δ(≤ 1.3) and Nc(∼500) are reasonable and, hence, we conclude that the numerical method proposed by Wetherill and Stewart is a valid and useful method for describing the planetary accumulation process. The real planetary accumulation process is more complex since it is coupled with the velocity evolution of the planetesimals. In the subsequent paper, we will complete the high-accuracy statistical code which simulate the accumulation process coupled with the velocity evolution and test the accuracy of the code by comparing with the results of N -body simulation.

Implicit SPH Drag and Dusty Gas Dynamics

In a multiphase flow such as gas and dust, the time step for numerical integration is controlled by the Courant condition, by the drag between the fluids, and by gravity. If the dust particles are sufficiently small the time step is determined by the drag. Since the primary effect of very large drag is to keep the fluids moving together, and the details of how this is achieved are not important when the drag is very large, it is desirable to treat the drag terms implicitly. For standard grid-based methods this is straightforward, but for particle methods like SPH it appears difficult because the drag appears in the form of pair interactions between particles, and any given particle can interact with ∼40 neighbours. The idea exploited in this paper is to treat each pair interaction separately. The velocities of the two particles involved are updated implicitly, according to their pair interaction, and the initial velocities of these particles are then replaced by the new velocities. This process is repeated for each pair interaction. Two sweeps over the particles gives satisfactory convergence. In this paper the method is tested first by keeping just the drag terms, then including pressure and drag to test wave propagation, and finally including gravity to consider the fall of a layer of gas in an isothermal atmosphere. The results are in good agreement with theory. The basic idea of working with implicit pair interactions can be extended to other SPH problems.