NICE—An explicit numerical scheme for efficient integration of nonlinear constitutive equations

The paper presents a simple but efficient new numerical scheme for the integration of nonlinear constitutive equations. Although it can be used for the integration of a system of algebraic and differential equations in general, the scheme is primarily developed for use with the direct solution methods for solving boundary value problems, e.g. explicit dynamic analysis in ABAQUS/Explicit. In the developed explicit scheme, where no iteration is required, the implementation simplicity of the forward-Euler scheme and the accuracy of the backward-Euler scheme are successfully combined. The properties of the proposed NICE scheme, which was also implemented into ABAQUS/Explicit via User Material Subroutine (VUMAT) interface platform, are compared with the properties of the classical forward-Euler scheme and backward-Euler scheme. For this purpose two highly nonlinear examples, with the von Mises and GTN material model considered, have been studied. The accuracy of the new scheme is demonstrated to be at least of the same level as experienced by the backward-Euler scheme, if we compare them on the condition of the same CPU time consumption. Besides, the simplicity of the NICE scheme, which is due to implementation similarity with the classical forward-Euler scheme, is its great Advantage.     

A second-order scheme for integration of one-dimensional dynamic analysis

This paper proposes a second-order scheme of precision integration for dynamic analysis with respect to long-term integration. Rather than transforming into first-order equations, a recursive scheme is presented in detail for direct solution of the homogeneous part of second-order algebraic and differential equations. The sine and cosine matrices involved in the scheme are calculated using the so-called 2^N algorithm. Numerical tests show that both the efficiency and the accuracy of homogeneous equations can be improved considerably with the second-order scheme. The corresponding particular solution is also presented, incorporated with the second-order scheme where the excitation vector is approximated by the truncated Taylor series.     

Large stencil viscous flux linearization for the simulation of 3D compressible turbulent flows with backward-Euler schemes

The purpose of this article is to study different approximate linearizations of the RANS equations viscous fluxes, for numerical simulations of compressible turbulent flows with backward-Euler schemes. The explicit convective flux under consideration is centred with artificial dissipation. The discrete viscous flux, calculated from cell-centred evaluation of the gradients, needs less computations and memory storage than other discretizations. But, in other respects, the balance of this numerical flux has a large stencil, which is not coherent with the 3-point per mesh direction stencil of classical implicit stages. Therefore 3-point and 5-point per mesh direction approximate linearizations are built from the thin layer flux formula. The stability condition of the corresponding backward-Euler schemes is given for a scalar linear equation (for the basic non-factored version of scheme and with LU-relaxation). Multigrid and monogrid computations of turbulent flow around two external configurations are performed with Wilcox’s k-ω turbulence model. The 5-point per mesh direction linearizations, coherent with the differential of the fluxes balance of thin layer approximation of explicit viscous fluxes, leads to the most efficient implicit stages.     

Energy consistent algorithms for dynamic finite deformation plasticity

This paper addresses the theoretical development and numerical implementation of energy consistent algorithms for dynamic elastoplasticity, emphasizing finite strain constitutive formulations so that unconditional stability of the algorithms is assured even in the fully nonlinear regime. The key concept behind energy consistency is the requirement that the discretized system obey an a priori stability estimate, which may be derived in general using the first and second laws of thermodynamics. This approach to computational dynamic plasticity differs from typical application of traditional algorithms (such as Newmark or Hilber–Hughes–Taylor-α methods), where local time integration schemes for plasticity laws are developed somewhat independently from the global time integration scheme for the equations of motion, without explicit consideration of thermodynamical restrictions. Two algorithms based on both additive and multiplicative finite deformation plasticity model are formulated within the energy consistent framework. Both algorithms possess the desirable feature of nonlinear stability of previous energy–momentum algorithms for elastodynamics.     

A three-step wavelet Galerkin method for parabolic and hyperbolic partial differential equations

A three-step wavelet Galerkin method based on Taylor series expansion in time is proposed. The scheme is third-order accurate in time and O(2^?jp ) accurate in space. Unlike Taylor–Galerkin methods, the present scheme does not contain any new higher-order derivatives which makes it suitable for solving non-linear problems. The compactly supported orthogonal wavelet bases D6 developed by Daubechies are used in the Galerkin scheme. The proposed scheme is tested with both parabolic and hyperbolic partial differential equations. The numerical results indicate the versatility and effectiveness of the proposed scheme.     

Tuning of PIλDμ controllers based on sensitivity constraint

A graphical tuning method for fractional-order PID (PI^λD^μ) controllers is proposed based on the sensitivity function constraint of the closed-loop, which provides the information on robustness to plant uncertainties. The stabilizing regions in integral-derivative plane of the controller are first identified using a graphical stability criterion applicable to fractional-delay systems. Then, via Leibniz Sector Formula, the stabilizing region is optimized with respect to the two fractional orders of the controller to expect bigger stabilizing regions. Finally, the sensitivity function constraint of the closed-loop is mapped into stabilizing region by means of the explicit algebraic equations which can be solved efficiently. Numerical examples of a second-order integrating delay process are followed in each design procedure to show the effectiveness of the method.     

Domain decomposition algorithm based on the group explicit formula for the heat equation

A finite difference domain decomposition algorithm (DDA) for solving the heat equation in parallel is presented. In this procedure, interface values between subdomains are calculated by the group explicit formula, whereas interior values of subdomains are determined by the classical implicit scheme. The stability and convergence for this DDA are proved. The stability bound of the procedure is derived to be eight times that of the classical explicit scheme. Though the truncation error at the interface is O(τ?+?h), L ²-error is proved to be O(τ?+?h ²). Numerical examples confirm the second-order convergence and indicate that the stability condition is sharp. A comparison of the numerical errors of this procedure with other known methods is also included.     

Adaptive element-free Galerkin method applied to the limit analysis of plates

The implementation of an h-adaptive element-free Galerkin (EFG) method in the framework of limit analysis is described. The naturally conforming property of meshfree approximations (with no nodal connectivity required) facilitates the implementation of h-adaptivity. Nodes may be moved, discarded or introduced without the need for complex manipulation of the data structures involved. With the use of the Taylor expansion technique, the error in the computed displacement field and its derivatives can be estimated throughout the problem domain with high accuracy. A stabilized conforming nodal integration scheme is extended for use in error estimation and results in an efficient and truly meshfree adaptive method. To demonstrate its effectiveness the procedure is then applied to plates with various boundary conditions.     

Runge-Kutta interpolants based on values from two successive integration steps

New interpolants of the explicit Runge-Kutta method for the initial value problem are proposed. These interpolants are based on values of the solution and its derivative from two successive integration steps. In this paper, three interpolants withO(h ⁶) local error (l.e.), for the fifth order solution, of the methods Fehlberg 4(5) (RKF 4(5)), Dormand and Prince 5(4) (RKDP 5(4)) and Verner 5(6) (RKV 5(6)) without extra cost are derived. An interpolant withO(h ⁷) (l.e.) for the sixth order solution of the Verner's method with only one extra function evaluation per integration step is also constructed. The above advantages are obtained without any cost in the magnitude of the error.     

Numerical stability of incompletely integrated shape functions

Finite element stiffness matrices formed by “unstable” integration schemes regain positive definiteness with a sufficient number of imposed boundary conditions. Moreover, the feared decline in the condition of these matrices does not materialize; their spectral condition number being O(h^?2m) as for the strictly variational elements.     

A Taylor Collocation Method for the Solution of Linear Integro-Differential Equations

In this study, a matrix method called the Taylor collocation method is presented for numerically solving the linear integro-differential equations by a truncated Taylor series. Using the Taylor collocation points, this method transforms the integro-differential equation to a matrix equation which corresponds to a system of linear algebraic equations with unknown Taylor coefficients. Also the method can be used for linear differential and integral equations. To illustrate the method, it is applied to certain linear differential, integral, and integro-differential equations and the results are compared.     

The Nørsett time integration methodology for finite element transient analysis

This paper presents a set of methods for time integration of problems arising from finite element semidiscretizations. The purpose is to obtain computationally efficient methods which possess higher-order accuracy and controllable dissipation in the spurious high modes. The methods are developed and analysed by a general collocation methodology which leads to the class of Nørsett approximants. An algorithmic parameter is used to achieve an effective control over numerical dissipation. Moreover, a simple and efficient implementation scheme is presented. At each time step, algorithms based on p-order collocation polynomials require the solution of p sets of linear algebraic equations with the same coefficient matrix. In this way, a single factorization is needed and no transformations are required to recover the approximate solution at the end or within the time interval. To demonstrate the performance of the proposed algorithms, a wide experimental evaluation is carried out on typical test problems in finite element transient analysis.     

Singularity-free time integration of rotational quaternions using non-redundant ordinary differential equations

A novel ODE time stepping scheme for solving rotational kinematics in terms of unit quaternions is presented in the paper. This scheme inherently respects the unit-length condition without including it explicitly as a constraint equation, as it is common practice. In the standard algorithms, the unit-length condition is included as an additional equation leading to kinematical equations in the form of a system of differential-algebraic equations (DAEs). On the contrary, the proposed method is based on numerical integration of the kinematic relations in terms of the instantaneous rotation vector that form a system of ordinary differential equations (ODEs) on the Lie algebra \(\mathit{so}(3)\) of the rotation group \(\mathit{SO}(3)\). This rotation vector defines an incremental rotation (and thus the associated incremental unit quaternion), and the rotation update is determined by the exponential mapping on the quaternion group. Since the kinematic ODE on \(\mathit{so}(3)\) can be solved by using any standard (possibly higher-order) ODE integration scheme, the proposed method yields a non-redundant integration algorithm for the rotational kinematics in terms of unit quaternions, avoiding integration of DAE equations. Besides being ‘more elegant’—in the opinion of the authors—this integration procedure also exhibits numerical advantages in terms of better accuracy when longer integration steps are applied during simulation. As presented in the paper, the numerical integration of three non-linear ODEs in terms of the rotation vector as canonical coordinates achieves a higher accuracy compared to integrating the four (linear in ODE part) standard-quaternion DAE system. In summary, this paper solves the long-standing problem of the necessity of imposing the unit-length constraint equation during integration of quaternions, i.e. the need to deal with DAE’s in the context of such kinematical model, which has been a major drawback of using quaternions, and a numerical scheme is presented that also allows for longer integration steps during kinematic reconstruction of large three-dimensional rotations.     

Automated algebraic integration of products of interpolation functions

This paper describes a method and a computer program for the explicit evaluation of the algebraic form of the integral $\text{∝ ∝ NN}{}^{\mathrm{<mtext>t</mtext>}}\text{d}\text{A}$ for triangular finite elements with straight sides, where N may be any vector of interpolation functions in which alphameric as well as numeric coefficients may occur. The method may be extended to any similar integration provided an explicit integration formula is available for the integral of a general polynomial term over the element area.     

Finite difference method for time-space-fractional Schrödinger equation

In this paper, an implicit finite difference scheme for the nonlinear time-space-fractional Schrödinger equation is presented. It is shown that the implicit scheme is unconditionally stable with experimental convergence order of O(τ^2?α+h²), where τ and h are time and space stepsizes, respectively, and α (0<α<1) is the fractional-order in time. In order to reduce the computational cost, the explicit–implicit scheme is proposed such that the nonlinear term is easily treated. Meanwhile, the implicit finite difference scheme for the coupled time-space-fractional Schrödinger system is also presented, which is unconditionally stable too. Numerical examples are given to support the theoretical analysis.     

On the complexity of algorithms for the translation of polynomials

Three algorithms for computing the coefficients of translated polynomials are discussed and compared from the point of view of complexity. The two classical translation algorithms based on explicit application of the Taylor expansion theorem and the Ruffini-Horner method, respectively, have complexityO (n ²). A third algorithm based on the fast Fourier transform is shown to have complexityO (n logn). However, when the cost of arithmetic operations is explicitly taken into consideration, the Ruffini-Horner algorithm is one order of magnitude better than the one based on the Taylor expansion and competes quite well with the algorithm based on the fast Fourier transform.     

Order-(n+m) direct differentiation determination of design sensitivity for constrained multibody dynamic systems

With the complexity and large dimensionality of many modern multibody dynamic applications, the efficiency of the sensitivity evaluation methods used can greatly impact the overall computation cost and as such can greatly limit the usefulness of the sensitivity information. Most current direct differentiation approaches suffer from prohibitive computational cost, which may be as great as O(n⁴+n²m²+nm³) (for systems with n generalized coordinates and m algebraic constraints). This paper presents a concise and computationally efficient sensitivity analysis scheme to facilitate such sensitivity calculations. A unique feature of this scheme is its use of recursive procedures to directly embed the algebraic constraint relations in forming and simultaneously solving a minimal set of equations. This results in far fewer operations than more traditional, or so-called O(n), counterparts. The algorithm determines the derivatives of generalized accelerations in O(n+m) operations overall. The resulting equations are exact to integration accuracy and enforce constraints exactly at both the velocity and acceleration levels.     

The finite volume element method with quadratic basis functions

The paper presents a box scheme with quadratic basis functions for the discretisation of elliptic boundary value problems. The resulting discretisation matrix is non-symmetrical (and also not an M-matrix). The stability analysis is based on an elementwise estimation of the scalar product <A _h u _h ,u _h >. Sufficient conditions placed on the triangles of the triangulation lead to discrete ellipticity. Proof of anO(h ²) error estimate is given for these conditions.     

Stability and error analysis of mixed finite-volume methods for advection dominated problems

We consider a convection-diffusion-reaction problem, and we analyze a stabilized mixed finite-volume scheme introduced in [1]. The scheme is presented in the format of discontinuous Galerkin methods, and error bounds are given, proving O(h^½) convergence in the L²-norm for the scalar variable, which is approximated with piecewise constant elements.     

A numerical scheme for integrating the rate plasticity equations with an “a priori” error control

An alternative scheme for the integration of incremental plasticity constitutive equations is presented and discussed. The method is based on a linear inequality theory approach and exploits the idea of the “a posteriori” local linearization of the yield surface. Results of both the flow theory of plasticity and the deformation theory can be obtained with a predetermined maximum violation of both the yield condition and the flow rule. Numerical results are shown, and compared to analytical solutions as well as to the results of other methods for the simple case of the von Mises yield function.