Iterative methods for solving nonlinear dynamic problems of construction mechanics

. The paper considers direct methods for integrating nonlinear equations of dynamics. Two approaches have been identified for approximating the equations during a small provisional segment. In the first, a more approximate approach is considered, in which the displacements are approximated only with the cross rigidity of the system at the current moment. In the second, an approach in displacement increments is considered, when the change in the system response is differentiated by the achieved level of displacements and their increments at the time step. Different schemes of forth integration for the nonlinear equations in the first and second approaches are considered. The absolute stability of the new proposed schemes of forth integration is proved when using a more accurate approximation in displacements. Theoretical results are supported by the model examples solution for which the exact solution, or the solution in the series is known.


Introduction
The main task of structural mechanics is the development of methods for calculating buildings and structures for strength, rigidity and stability, which, with sufficient efficiency of structures, would ensure the safety of structures throughout their entire service life. Today, increased requirements are imposed on the calculation methods -the calculation of structures must give an exhaustive forecast of their work at all stages of the life cycle. To solve such a problem, it is necessary to use the refined computational models in order to adequately reflect all the features of the system response to external disturbances in a numerical experiment. Therefore, the development of effective methods for modeling the dynamic response of nonlinear problems, including the problems of fracture mechanics, is an urgent task of modern structural mechanics.
The classification of nonlinear problems in mechanics today considers three factors that determine the nonlinearity of the resolving equations. First typephysical nonlinearity. The refined non-linear relationships between deformations and stresses, or between their increments, are considered here. Rejection of the Navier hypothesis about the displacements' smallness leads to the refinement of the geometric dependencies and the calculation in geometrically nonlinear staging. Classification of the nonlinear problems of V.V. Novozhilov supplemented by the concept constructive nonlinearity (for the first time this type of nonlinearity was considered by S.P. Timoshenko), which takes into account the change in the design scheme of the structure at certain stages of its existence. An important aspectfracture mechanicscan be considered as a special case of physically and constructively nonlinear formulations.
A fairly complete consideration of the calculated structures' geometry features, their mutual influence, physical and geometric dependencies today is possible only with the use of numerical calculation methods. The solution of dynamic problems in a nonlinear formulation leads to the fact that in the method of principal coordinates it is necessary to repeatedly determine the frequencies and forms of free oscillations. Therefore, due to the high complexity of solving the spectral problem, even in a truncated formulation, the principal coordinates method complexity catastrophically increases.

Basic Provisions
Recently, to determine the dynamic response, forth integration schemes have been used, in which a series of problems are successively solved at small time steps, in which the solving equations are algebraic. In this case, at one time step, the complexity of the solution corresponds to the statics problem. No need to solve more time-consuming spectral problem.
In problems of the linear system elements' rigidity dynamics, the equilibrium equations of a small element and the continuity equations do not change for the entire time interval in which the solution is considered. However, in the nonlinear problems of dynamics, the parameters of the basic equations change at each time layer. First, the parameters of the linearized (within a fine pitch) the equation of motion are determined. These parameters, as it were, "freeze" and are considered constant within the time layer. After determining the components of the system stress-strain state at the end of the time step, the motion equations parameters corresponding to the next pitch are determined. That is, within onetime step, a quasilinear problem with current parameters corresponding to the stress-strain state at the beginning of the time step is considered. The simple loading hypothesis, used in many versions of the theory of plasticity, is realized within one-time step, but in general, the process of complex loading of the system is realized.
One of the variants of the equations of motion linearized at a time step for an ensemble of finite elements has the form (1): In (1) and further M is the mass matrix, С is a damping matrix, Кс, n defines "Cross" stiffness matrix at the time step n of the ensemble of finite elements, q(t) is the nodal displacement vector function, P(t) is the vector-function of the load reduced to the nodes of the system. The matrix Кс, n is formed at the beginning of the time step according to the achieved by the moment t=tn parameters of the system stress-strain state. It is assumed that for the time segment [tn, tn+1] the overall rigidity of the system does not change. Therefore, at a time step, the potential energy of deformation is approximated by a quadratic function (2), and therefore the motion equations are represented in the form (1): The second form of the motion equations is obtained on the basis of the total potential energy functional expansion in the vicinity of the achieved by the time instant tn the level of the stress-strain state in the Taylor series with retention of terms up to the second order inclusive and is written in the form (3): In the approach used in (3), the linearization is performed in advance in a variational form, therefore, the general algorithm of the finite element method for the formation of the "cross" Kс,n and "tangent" Kк,n system stiffness matrices (4): The disadvantage of (1) in comparison with (3) is the worse quality of the system response approximation at a time step. Indeed, it is easy to show that in (1) the increment of the system response is determined by the "cross" stiffness, which is not quite correct: However, in the form, within the time step (1) corresponds to the motion equations of a linear system. Therefore, for numerical integration (1), any stable method can be used, such as Newmark, theta-Wilson, etc.
It should be noted that if at the dawn of the forth integration methods' development, the central issue was the approximation accuracy improvement at the time step, then later the issue of the numerical stability of the methods became decisive. Moreover, many methods that use the best in accuracy approximation of the displacement vector function have a region of instability, and their application is limited. Example is the conditionally stable central difference method.

Stability of the forth integration configuration
Therefore, further we will consider an absolutely stable method of integration of nonlinear equations of motion (3), presented here in the form (5). It was shown earlier that the equation (5) can be interpreted as an extension of one of the best today Newmark's method for nonlinear problems.  n (5) the parameter β is assigned from the condition of the absolute stability of the configuration. That is, in contrast to the classical form of Newmark's method, a directed search for the stability parameter is performed in (5). Previously, not entirely justifiably, the parameter β was recommended to assign in the form of the fundamental tone circular frequencies squares ratio when oscillating the linear problems with "cross" and "tangent" rigidity (6): However, this result (the frequencies' squares ratio) was obtained for an oscillator, and the extension of the formula for a system with several degrees of freedom is unreasonable.
In this paper, it is shown that for the system with several degrees of freedom, the stability parameter β should be determined according to (7): It is shown in this work that such an assignment of the stability parameter according to (7) corresponds to the Neumann criterion of absolute stability: the spectral radius of the integration scheme does not exceed unity.
Indeed, the stability check of the integration scheme is performed by considering the free oscillations of a nonlinear system in the absence of damping. System (5) in this case can be written in the form (8):  Example 1. Let us consider an illustrative free vibrations example of a flat bracket made of a physically nonlinear material (Fig. 1). The number of dynamic degrees of freedom 2. In the equilibrium position, a static force is applied to the bracket Р, from which, at rest, stresses arise in the rods: σ1=1.463σТ, σ1=1.837σТ. Value "Tangent" module for both rods 0.1Е, "Cross" modulus is 0.26Е and 0.2Е for the first and second rods, respectively. Filling "Tangent" and "Cross" stiffness matrices, we get:   Table 1 the values of the stability parameter assigned by different methods are given. In Table 2 the resolving operator spectral radius values are given for different values of the time step.

Model examples
For (14), the relationship between the amplitude and the period of oscillations is known (15): In (15) А is an amplitude, ТА is an oscillation period, F() is a normal elliptic Legendre integral of the first kind.
The approximation of the motion equation (3) A series of test calculations was performed at different amplitudes of free oscillations and different pitches of integration over time. Table 3 shows a selection of the main results, in each cell two values are given -the absolute value and the relative error (in percent).