Stability of direct circuits integrating the equations of motion in the simulation of the dynamics of destruction

In this paper, it is discussed the approximation of the nonlinear equations of motion of mechanical systems used on different operators achieved level displacements and displacement increments. For the numerical solution of the problem, we used the scheme of direct integration of the equations of motion, which in the case of a linear problem reduces to the known Newmark scheme. The stability parameter of the integration scheme is determined. Analytical calculations are confirmed on the test meeting of the model task.


Introduction
The main task of structural mechanics is the development of such methods of analysis of structures for strength, rigidity and stability which with sufficient economic efficiency would ensure the safety of structures during the entire period of their life.Today, the increased requirements are applied to the methods of calculations -the analysis of structures should provide a comprehensive outlook on their performance at all stages of the life cycle.To solve a similar problem, it is necessary to use refined calculation models in order to adequately represent all the features of the system response to external perturbations in a numerical experiment.Thus, developing of effective methods of modeling of the dynamic response of non-linear problems, including problems of fracture mechanics, is an urgent task of modern structural mechanics.
The classification of non-linear problems of mechanics today considers three factors that determine the nonlinearity of the resolving equations.The first type -physical nonlinearity.Here the refined nonlinear relations between deformations and stresses, or between their increments, are considered.Rejection of the Navier hypothesis of small displacements leads to a refinement of geometric constraints and calculation in geometrically nonlinear statement.Classification of nonlinear problems of K. Novozhilov is complemented by the concept of constructive nonlinearity (at first time this kind of nonlinearity was presented by S. Timoshenko), which takes into account the change in the construction scheme at certain stages of its existence.The important aspectfracture mechanics -can be considered as a special case of physically and structurally nonlinear creations.
A sufficiently comprehensive consideration of the features of geometry of the calculated structures, their mutual influence, physical and geometric dependencies could only be possible today when using numerical calculation methods.Solving of dynamic problems in a nonlinear design leads to the fact that in the method of principal coordinates it is necessary to measure many times the frequency and shape of free oscillations.Therefore, because of the high complexity of the process of solving the spectral problem, even in a truncated formulation, the complexity of the principal coordinate method increases catastrophically.
Recently, in order to find dynamic response, the direct integration schemes are being used, in which a series of problems is being solved sequentially at small time phases, where the resolving equations are algebraic.In this case, at one phase in time, the complexity of the solution corresponds to the statics problem.Solving more complex and timeconsuming spectral problem here is not necessary.
In problems of the linear dynamics of the rigidity of the elements of the systems, the equations of equilibrium of a small element and the continuity equation do not change for the entire time interval where the solving of a problem is taking place.However, in nonlinear dynamics problems, the parameters of the basic equations change on each time layer.First the parameters of the linearized (within a small step on the time) equations of motion must be defined.These parameters seem to be "frozen" and are considered constant within the temporal layer.After determining the components of the stressstrain state of the system at the end of the time phase, the parameters of the equations of motion corresponding to the next phase are being determined.That is, within one time phase, a quasilinear problem with current parameters corresponding to the stress-strain state at the beginning of the time phase is considered.The simple loading hypothesis, which is used in many versions of the theory of plasticity, is being implemented within one time phase, but, in general, the process of complex loading of the system is applied.
One of the options of the equations of motion for an ensemble of finite elements linearized in a time phase has the form (1): and onwards M -is a mass matrix, C -a deformation matrix, K, n -a "secant" rigidity matrix for phase n of the time for the ensemble of finite elements, q (t) -a vector-function of nodal displacements, P (t) -the load vector function given to the nodes of the system.Matrix K, n is formed during the first time phase in accordance with the parameters of the stress-strain state of the system achieved at a time t = t n .It is assumed that in the time segment [t n, t n +1] the overall rigidity of the system does not change.Thus, at the time phase the potential energy of deformation is approximated by the quadratic function (2), and therefore the equations of motion are represented in the form (1): The second form of the equations of motion obtained on the basis of the decomposition of functional of the total potential energy in the range achieved at time t n of the level of stress-deformed state into a Taylor line with retention of terms up to second order and is presented as following For the approach used in (3), the linearization is designed in advance in the form of variation, so the overall algorithm of the method of finite elements forming "secant" K, n and "tangent" rigid matrices of the system is sustained (4): The disadvantage of (1) in comparison with (3) is a lower quality of approximation of the system reaction at a time phase.Indeed, it is straightforward showing that in (1) the augmentation in the reaction of the system is determined by the "secant" rigidity, which is not utterly correct: However, according to the formula, within the time phase (1), it corresponds to the equations of motion of the linear system.Therefore, for any numerical integration (1) any sustainable method can be used, such as Newmark, Tett-Wilson, and others.
It should be noted that if at the dawn of the development of direct integration methods the question of improving the accuracy of approximation at the time phase was considered to be dominant, later the question of numerical stability of methods became critical.Moreover, many methods that use the best approximation of the displacement vector function have an instability area, and their application is limited.A conditionally stable method of central differences can serve as an example.
Thence, we consider absolutely stable method for integrating the nonlinear equations of motion (3), here presented in the form (5). Earlier it was demonstrated that equations (5) can be handled as an extension of one of the best Newmark's methods for nonlinear problems.
( ) In (5), the parameter β is assigned due to the absolute stability conditions scheme.That is, in contrast to the classical form of the Newmark method, a directed search for the stability parameter is performed in (5).Previously, with no major justification, it was recommended that parameter β is defined as a form a ratio of squares of circular pitch frequencies in the oscillations of linear problems with the "secant" and "tangent" rigidity (6): However, this result (the ratio of the squares of frequencies) was obtained for the oscillator, and the propagation of the formula for a system with several degrees of freedom is unreasonable.
This paper shows that for a system with several degrees of freedom of the stability the parameter β should be defined according to (7) [ ] [ ] ( ) In this paper, it is demonstrated that such an assignment of the stability parameter according to (7) corresponds to the Neumann absolute stability criterion: the spectral radius of the integration scheme does not exceed unity.
Indeed, the stability test is carried out by considering the free oscillations of a nonlinear system in the absence of damping.The system (5) in this case can be written in the form (8): The matrix R in (8) is symmetric positive definite.We perform the splitting (8) with respect to the eigenvectors of the matrix R -1 Kc, n .Let μi be the eigenvalues of the matrix R -1 Kc, n, we obtain for the ith splitting form (9): For each eigenvalue μi, the representation (10) holds for the eigenvectors Zi of the matrix R -1 Kc, n: For the purpose of the assignment of the stability parameter β, according to (7), the Rayleigh inequality holds , which implies the restriction (12) of the eigenvalues of the matrix R -1 Kc, n : Obviously, that when the restriction (12) is satisfied, the spectral radius of the matrix operator in (9) does not exceed unity, i.e.The Neumann absolute stability condition for the scheme of direct integration of the nonlinear equations of motion (5) is satisfied.
Numerically, for a system with a large number of degrees of freedom, the realization of formula (7) can be quickly carried out by the Lyusternik iterative method.
Note that the matrix [Kk,n] -1 Kc,n is not symmetric, it takes a large amount of memory and the laboriousness of its construction is great.In order not to form the matrix [Kk,n] -1 Kc,n, the following modification is suggested within the limits of one iteration of the Lyusternik algorithm (13):

Experimental part
We consider model illustrative examples for nonlinear systems with several degrees of freedom.The stability of the solutions is shown according to (7) and the instability according to (6).Example 1: Consider an illustrative example of free oscillations of a flat bracket made of a physically nonlinear material (Fig. 1).The number of dynamic degrees of freedom equal 2. In a state of equilibrium when working in rest mode, the static force P, from which in a state of rest in a pair of voltages: σ1 = 1.463σТ, σ1 = 1.837σТ.The value of the "tangent" module for both rods 0.1E, the "secant" module is 0.26E and 0.2E for the first and second rods, respectively.Filling the "tangent" and "secant" stiffness matrices:  [ ] [ ] Table 2 shows the values of the spectral radius of the resolving operator for different values of the time step.
Table 2. Values of the spectral radius.As can be seen, only the stability parameter, assigned according to the proposed procedure, ensures the fulfillment of the Neumann criterion over the entire range of the step of integration.Example 2: Oscillations of a system with nonlinear restoring force.We consider the numerical solution of the known model equation: For ( 14), the relation between the amplitude and the oscillation period is known: In (15), A is the amplitude, ТА is the oscillation period, and F() is the normal elliptic Legendre integral of the first kind.Approximation of the equation of motion (3) for this example has the form (16): A series of test calculations is performed for different amplitudes of free oscillations and a different step of integration over time.Table 3 shows the selection of the main results, in each cell two values are given -the absolute value and the relative error (in percent).

Conclusions
The circuit has both a relatively amplitude errors and phase error.The phase error exceeds the value of the error in the amplitudes.It is expressed in an increase in the period of free oscillations in the approximate solution.Amplitude accuracy reduces the calculated amplitude relatively accurate.Both errors decrease rapidly with a decrease in the integration step.With the improved stability parameter of the computational process, we obtain a stable numerical solution.
Parameter assignment algorithms used above stability is not guaranteed stability of the iterative process.

Table 1
shows the values of the stability parameter, assigned according to different methods.

Table 1 .
Values of the stability parameter.

Table 3 .
Results on free oscillations of a nonlinear system.