On conservation of energy and kinematic compatibility in dynamics of nonlinear velocity-based three-dimensional beams

In this paper, we present an original energy-preserving numerical formulation for velocity-based geometrically exact three-dimensional beams. We employ the algebra of quaternions as a suitable tool to express the governing equations and relate rotations with their derivatives, while the finite-element discretization is based on interpolation of velocities in a fixed frame and angular velocities in a moving frame description. The proposed time discretization of governing equations directly relates the energy conservation constraint with the time-discrete kinematic compatibility equations. We show that a suitable choice of primary unknowns together with a convenient choice of the frame of reference for quantities and equations is beneficial for the conservation of energy and enables admissible approximations in a simple manner and without any additional effort. The result of this study is simple and efficient, yet accurate and robust numerical model.


Introduction
Geometrically nonlinear beam models are popular due to their suitability for modeling frame-like structures undergoing arbitrary large displacements and rotations but small strains with sufficient accuracy and relatively low computational costs. Despite the simplified kinematics, the governing equations remain challenging from the perspective of numerical solution methods. Computationally efficient, accurate, robust and stable numerical formulations demand firm understanding of the model, relationships between the quantities describing the beam and their dependence with respect to both space and time. Modern techniques in numerical analysis allow many possible approaches in handling the problem of such complexity. It is thus not a surprise that after more than three decades of intensive research this field is still a subject of interest for many researchers, as reflected by recent publications, see, e.g., [1][2][3][4][5][6][7][8][9][10][11][12].
Many challenges in beam formulations reported in the literature stem from the use of spatial rotations as members of the configuration space of the beam. Spatial rotations form a multiplicative group with interesting properties that are not easily preserved in numerical solution methods. Since the pioneering work of Simo and Vu-Quoc [13], beam finite elements have often been based on the interpolation of rotational degrees of freedom. In [13], a rotational vector is used for the parametrization of rotations, but this is only one among a variety of choices, see [14]. The chosen parametrization of rotations has an important impact on their dis-cretization and the behavior of solution methods. Thorough discussions and comparisons can be found, e.g., in [15][16][17][18][19]. Crisfield and Jelenić [20,21] pointed out that standard additive-type interpolation of rotational vectors could lead to nonobjective strain measures in many rotation-based approaches. Different interpolation strategies for rotational degrees of freedom were studied by Romero [22] and Bauchau and Han [23]. The use of only three parameters for the description of rotations results in lower storage demands for the primary unknown quantities. However, we can not easily avoid the calculation of rotation matrices and computationally more sensitive extraction of rotational parameters from rotation matrices [24]. This concept can be completely avoided by the introduction of rotational quaternions as a four-parameter representation of rotations. Not only that such parametrization is free from singularities, but also it is computationally more efficient regarding the number of calculations needed to rotate a vector. The use of quaternions seems to be standard for the efficient implementation of algorithms involving rotation matrices, see [54, p. 281] and the references therein. In beam formulations, quaternions were first employed by Bottasso [25] and Kehrbaum and Maddocks [26], while McRobie and Lasenby [27] used theoretically equivalent Clifford algebra. The interest in quaternion algebra as a suitable tool for beam formulations has increased in the last decade, see [2,[28][29][30][31][32][33]. Another alternative is the interpolation of all nine components of the rotation matrix by Betsch and Steinmann [34] and recently by Sonneville et al. [35].
The description of the relationship between the deformed configuration of the beam and the strain measures is an important aspect of beam theories. In the "geometrically exact" approach, see Reissner [36,37], strain resultants of the cross section are introduced and related to configuration variables in such a way that the kinematic equations are consistent with the virtual work principle. Several alternatives to this approach can be found, such as co-rotational method [38] and absolute nodal coordinate method [39]. In the present paper, we will stem from geometrically exact model as a suitable basis for the proposed numerical formulation.
In dynamics, a proper treatment of the configuration space with three-dimensional rotations might still result in a loss of numerical stability in the case of long-term calculations [40]. Several energy-preserving and energy-decaying algorithms were proposed as the answer to this phenomenon, e.g., [40][41][42][43][44][45][46]. These algo-rithms use similar strategy within different frameworks. For a more general approach, the reader is referred to the paper by Bauchau and Bottasso [47]. A common characteristic of these formulations is the special approximation of the kinematic equations imposed by preservation constraints. These approximations can be in contradiction with other properties of continuous systems, i.e., the orthogonality of rotation at the midtime is often violated. Despite this shortcoming, increased stability of long-term calculations using such methods was confirmed numerous times. Another important property of a continuous dynamic system, surprisingly rarely considered in beam formulations, is the direct relation between the strains and the velocities without any configuration variables being present, called the compatibility equations [48] or, alternatively, the intrinsic kinematic equations [49].
In the present paper, we will employ the rotational quaternions for the description of three-dimensional rotations. However, in contrast to our quaternion-based approach presented in [32], we will avoid the interpolation of rotational quaternions due to the abovementioned problems. The crucial idea exploited here assumes the spatial and temporal derivatives of configuration variables to be the natural quantities for the description of total mechanical energy of the system. The previously developed strain-based beam formulation by Zupan and Saje [50] for static analysis of three-dimensional beams was motivated by this idea. In dynamics, angular velocities have a similar role as the rotational strains in statics, which was exploited in the paper by Zupan and Zupan [51], where velocities and angular velocities were taken as the primary unknowns of both spatial and temporal discretization. The time discretization employed in [51] can be interpreted as a modification of the implicit Newmark scheme. Unfortunately, a more suitable choice of the primary variables does not automatically assure the energy preservation. For sparse meshes and larger time steps, we can observe in some cases a loss of convergence accompanied with unrealistic increase in the total mechanical energy [51]. To avoid this phenomenon, we here propose a novel energy-conserving method based on a more suitable form of discrete governing equations while preserving the advantages of the velocity-based approach. Three important computational benefits of our approach are: (i) the additivity of the primary unknowns enables the use of standard interpolation techniques in space and simplifies the iteration procedure; (ii) the kinematic compatibility equations are satisfied with the same level of accuracy as the governing equations; and (iii) the energy of conservative systems is preserved without the need to introduce any special approximation of rotation or strain field. The computational advantages of the quaternion representation of rotations are preserved; moreover, proper treatment of the primary unknowns and the governing equations provides us with a relatively simple, numerically stable and robust method. The accuracy and excellent numerical performance of the present approach will be demonstrated by several benchmark examples.

Governing equations
Among several possible approaches in modeling deformable bodies with one dimension more significant than the other two, we base our numerical formulation on Cosserat geometrically exact theory of rods [52]. The Cosserat rod model describes the beam with a reference curve (often, but not necessarily, the line of centroids), and rigid cross sections attached to each point of this curve. Such description is capable of representing an arbitrary deformation of a beam including bending, extension, shear and torsion. We introduce two coordinate systems: One is fixed in space and time, the other rigidly attached to the cross section. The corresponding base vectors will be denoted by { g 1 , g 2 , g 3 } for the fixed and { G 1 (x, t) , G 2 (x, t) , G 3 (x, t)} for the local frame of reference, see Fig. 1. In the following, t denotes the time and x ∈ [0, L] the arc-length parameter of the reference curve at the initial state of the beam with L being its initial length. An arbitrary configuration of such a beam can be described by the position vectors r (x, t) and the rotations between fixed and local bases. We will use the rotational quaternions as a suitable representation of rotations. Taking ϑ to be the angle of rotation and e denoting the unit vector on the axis of rotation, a rotational quaternion q is expressed as The first part of (1) is written in the polar form while the second part employs the quaternion exponential to get the same result. After introducing the conjugated From (2), it follows that rotational quaternions can be used for describing rotations in R 3 . This is admissible since R 3 is a linear subspace of the space of quaternions. While an arbitrary quaternion is usually represented as (s ∈ R is the scalar part and v ∈ R 3 is the vector part), the quaternions with zero scalar part can be identified with vectors: Quaternions with zero scalar part are called pure quaternions. In what follows, abstract vectors and abstract quaternions will be replaced by one-column representations and denoted by bold face symbols. The hat over the symbol will be omitted when the first component equals zero since the appropriate size of one-column representation will be evident from the context. Such notation will enable a clear distinction between pure and arbitrary quaternions. In component description, the rotational quaternion also serves for transformation between the fixed and local basis representation and vice versa: Here and henceforth, we use lowercase letters for the component representations with respect to the fixed basis and uppercase letters for the local basis representations. All equations in the sequel will be written in terms of quaternion algebra without presenting all the details. A more detailed derivation of equations of a three-dimensional beam in terms of quaternion algebra can be found in [2,25,29,31].
It should be noted that all vector quantities in the sequel take values on the current configuration of the beam. The initial length of the line of centroids is used only as a suitable parameter for computational purposes. The transformation between the spatial and material form of equations is straightforward since the pull-back and push-forward are performed by the quaternion rotation (3). For further details, see the discussion in [37,53].

Kinematic equations
In a geometrically exact beam theory, the resultant strains of the cross section are introduced and related to the configuration variables-positions of the reference curve and rotations of the cross sections. A thorough introduction, derivation and interpretation of these equations, crucial for the model was given by Reissner [36] and Simo [37]. We will only briefly present the kinematic equations expressed in terms of quaternion algebra. The translational strain vector, , consists of an extensional strain and two shear strain components. It is defined by the following equation where 0 is determined from known strains, position vectors and rotations at the initial configuration of the beam. When we have no initial strains and the cross sections are initially orthogonal to the axis of the beam, we get 0 = [−1, 0, 0] T . The prime denotes the derivative with respect to x and the coordinate transformation is carried out using rotational quaternions. Thus, translational strain describes the rate of change of the position vector along the length of the beam. The rotational strain, K, consists of a torsional and two bending strain components. In terms of quaternions, it is determined by Vector K describes the rate of change of the local basis with respect to x and is expressed with respect to the local basis. The analogy of the above formula to the standard definition using rotation matrices is evident, while the factor 2 that additionally appears arises from the definitions (1)- (2). In dynamics, the measures for the rate of change of configuration variables with respect to time are also needed. For the reference curve, we simply have v =ṙ, describing its velocity with respect to the fixed basis, while the angular velocity of the cross sections is defined as The dot above a symbol ( · ) denotes the time derivative. For the angular velocities, the description in local basis is more convenient, as we will show later. From the above definitions, it is evident that the strain measures and the velocities are mutually dependent. By comparing mixed partial derivatives of (4)- (7), we geṫ The significance of equations (8)-(9) was pointed out by Antman [48] and Hodges [49]. They are called compatibility equations and will play an important role in our numerical formulation.

Equations of motion
When expressed in terms of quaternion algebra, the balance equations of a three-dimensional beam read and the boundary conditions are: Here, n and m represent the stress-resultant force and moment vectors of the cross section with respect to the fixed basis;ñ andm are the external distributed force and moment vectors; f 0 , h 0 , f L and h L are the external point forces and moments at the two boundaries, x = 0 and x = L. ρ is mass per unit of the initial volume; A is the area of the cross section; J ρ is the mass-inertia matrix of the cross section. The kinetic energy of a beam consists of the translational part representing the inertia of motion of the reference axis and the rotational part representing the inertia of motion of cross sections: Note that the second term is written with respect to the local frame.

Constitutive equations
The stress resultants depend on the strain vectors and K. We will assume here that the stress-strain relationship is linear in the strain measures, leading to The matrix C = C C K C K C KK is called the crosssectional constitutive tangent matrix. It can be fully populated, but we will assume it is symmetric. N and M represent the stress-resultant force and moment vectors of the cross section in the local basis. For the sake of completeness, let us express the stress resultants in both frames using the coordinate transformation (3): Now, we define the strain energy of a linear elastic beam as An important property of the continuous system of equations is the conservation of the total mechanical energy, defined as W = W K + W D , for conservative external loads. As pointed out by many authors, this property is not automatically preserved after the discretization of equations, in spite of that being a highly desirable goal due its direct correlation to the increased long-term stability of calculations.

Numerical formulation
We are dealing with a system of nonlinear partial differential equations. The governing equations will thus be discretized with respect to time and space. In contrast to the widely used beam models where the configuration variables are usually interpolated, the velocities in fixed basis and the angular velocities in local basis are here employed as the primary unknowns. Despite its nonstandard nature, such an approach has firm theoretical grounding. Note that velocities and angular velocities directly appear in the description of kinetic energy (16) and are thus crucial for its evaluation. Another reason stems from noncommutativity and multiplicative nature of spatial rotations, which restrict and complicate the admissible treatment of rotational degrees of freedom. A multiplicative update or a special transformation of rotational increments that are (locally) additive is therefore often used in beam formulations. To avoid this, we here exploit the well-known property of additivity of the components of angular velocity vectors, when expressed in the local basis. We can then use standard interpolation functions in Euclidean vector spaces without any inconsistency introduced and the increments can be simply added to the current values in the iteration procedure. Instead of working on manifolds, most of the computational effort is done in vector spaces, which considerably reduces the complexity of calculations. Due to our choice of primary variables, it is convenient to express the angular momentum balance equation (11) in the local basis: After considering (4)- (5) and (20), we get its final form: The boundary conditions for internal moments need to be expressed in the local basis as well. By direct use of coordinate transformation, we have

Time discretization
After employing the midpoint rule [54], equations (10) and (23) The superscript denotes the time at which a particular quantity is evaluated, h = t n+1 −t n is the time step, and time t n+1/2 = t n + h/2 denotes the midtime between t n and t n+1 . The right-hand side of (26)- (27) consists of the values of the quantities of the beam at the midtime. The way they are evaluated is not yet fixed (in contrast to standard midpoint rule), leaving us enough freedom to satisfy additional constraints. We will especially focus on the newly introduced quantities n and M, which represent the approximations of stress resultants at the midtime. They will be introduced in Sect. 3.3. The kinematic compatibility equations (8) and are therefore satisfied with the same accuracy as the equations of motion. Observe that no accelerations appear in our discrete governing equations (26)- (29), which is an important computational advantage.
In accord with equations (26)- (27) and (28)- (29), we need to express all the quantities of the beam at the middle of the time step. Velocities and angular velocities at t n+1/2 are determined by averaging the values at the boundaries of time interval t n , t n+1 : For the configuration variables, we assume which is again a second-order approximation of exact relations (6)- (7). Additionally, multiplicative formula (31) is kinematically consistent as it preserves the unit norm of rotational quaternions. After applying (30)-(31) in (4)-(5) at times t n and t n+1/2 , we get Assuming that all the quantities of a beam are known at t n , equations can be expressed with average velocities and angular velocities, v and as the only unknown functions along the length of the beam.

Spatial discretization
Functions v(x, t) and (x, t) are at each time t replaced by a set of their space-discrete values v p (t) and p (t) and interpolated by a set of N interpolation functions The points x p ; p = 1, . . ., N , are chosen from the interval [0, L] with x 1 = 0 and x N = L. Thus, the sets of discrete vectors v p and p become the unknowns of the problem. In accord with the time discretization presented above, the only unknowns of the present formulation at each time step are the average velocities and angular velocities v p and p at interpolation points.
Equations (26)- (27) are, in accord with the method of weighted residuals, multiplied by test functions I p (x), p = 1, . . . , N , and integrated along the length of the beam. Some terms are integrated by parts and the boundary conditions are considered to get the final discrete set of governing equations: where It is evident that the governing equations are nonlinear and will be therefore solved iteratively. Before we present the solution procedure, we still need to find suitable formulae for n and M. They will be derived from the energy preservation constraint.

Energy conservation
The proposed method is general regarding the load conditions, but for the special case of conservative system the total mechanical energy can be algorithmically preserved. Let us start from the difference of kinetic energy between two successive times of the scheme Considering the chosen interpolation (34)-(35), we have When there are no external forces acting on a beam, the discrete equations (36)- (37) are simplified accordingly, leading to A rearrangement of terms gives and after considering the discrete kinematic compatibility equations (28)-(29), we get To preserve the energy of a conservative system within the time step, we should have Comparison of (38) and (39) leads to the admissible stress resultants. To preserve the total mechanical energy in conservative problems, the approxima-tive stress-resultant vectors need to be evaluated by the following formulae: For the simplest case where C is block-diagonal, we have In order to preserve the energy of conservative systems we only need to evaluate the internal stresses from the average strains. This concludes the derivation of the present model.
It is obvious that the use of above formulas for the evaluation of midtime stress resultants does not restrict the applicability of our approach for nonconservative problems. Let us stress that no special approximation of rotations or strains was needed in the present approach. Due to a suitable choice of primary variables and governing equations, we get the desired properties of the formulation, while the kinematic compatibility is satisfied with the same accuracy as the governing equations.

Computational aspects
At each discrete time, the nonlinear discrete governing equations (36)- (37) are solved for the nodal values of average velocities and angular velocities using the Newton-Raphson method. This requires the construction of the Jacobian matrix, which is in our case expressed analytically. The analytical Jacobian is advantageous compared to the approximative one due to the presence of rotational degrees of freedom. Its derivation is relatively simple as the primary variables were conveniently chosen.

Linearization of equations
The linearization of rotational quaternions demands the consideration of their multiplicative and noncommutative nature. In order to linearize q [n+1/2] with respect to the following directional derivative is expressed: After a straightforward derivation, we get where T 0 and T 1 are the 4 × 3 matrices . The remaining terms are varied in a standard manner. The details are therefore omitted. At each step of the iteration procedure, we solve a linear set of equations where J (i) is the Jacobian tangent matrix, f (i) the residual vector of discretized equations (36)- (37), both in iteration i, and δy a vector of corrections of all nodal unknowns Our choice of variables allows the iterative corrections to be directly added to the current values: at each nodal point of the structure.

Prediction of initial values
To begin the iteration procedure, we need to predict the initial values. A natural assumption that initial velocities are equal to the converged values of the previous time: is applied. The midtime configuration and strains are then evaluated in accord with this assumption and approximations (30)-(33).

Compatibility of boundary conditions
We wish to assemble at the structural level elements with different initial rotations that are connected at rigid joints as simple as possible. Let two elements have different initial orientations at a joint, described by rotational quaternions q I 0 and q II 0 , see Fig. 2. The total rotation is then a composition of initial and superimposed rotation where the superimposed rotations at the node must be equal:

Fig. 2 A rigid joint of two elements with different initial orientations
For the present approach, this yields v I = v II and I = II , as the angular velocities are expressed in different local frames: After considering (46) and (47), we have We can enforce the condition (48) using the method of Lagrangian multiplies, but it is computationally more efficient to recognize that the angular velocities can be expressed as By choosing an appropriate basis for the description of quantities and equations in the local coordinate system, we can thus simply use a standard Boolean identification of unknowns at the structural level.

Energy dissipation
It is often desired that a numerical scheme allows dissipation of the high-frequency response. The present scheme can be easily extended to have dissipative properties. Starting from the energy-preserving stressresultant approximations (40)-(41) dissipative terms can be added in the following manner We can interpret the above approximation as the decomposition of stress resultants to conservative and dissipative part. The dissipative part is proportional to the increment of strain measures and β ∈ [0, 0.5] is a parameter related to the magnitude of dissipation. For conservative problems, we now have which confirms that the contribution of additional terms to the total mechanical energy is always negative.

Numerical simulations
Several numerical examples will be presented to demonstrate the ability of proposed formulation to solve problems with large displacements and rotations on long time intervals. We will focus on welldocumented and widely studied problems from the literature with the main focus on computational performance of the formulation. The problems under consideration are classical benchmark examples with C K = C K = 0, E and G are the elastic and shear moduli of material. The remaining quantities represent the geometric properties of the cross section expressed with respect to its centroid: A 1 is the area; J 1 is the torsional inertial moment; A 2 and A 3 are the effective shear areas; J 2 and J 3 are the bending inertial moments. All problems involve deformable beam members. In first three examples, a part of the motion is conservative, while the last one is a nonconservative problem. The interpolation points in all of the examples were taken to be equidistant and standard Lagrange polynomials were used as shape functions. Integrals were evaluated numerically using the Gaussian quadrature rule. The number of integration points was taken to be equal to the number of interpolation points, N , for full integration and N − 1 for the reduced integration. Newton-Raphson iteration was terminated when the Euclidean norm of the vector of all unknowns at the structural level, δy 2 , was less than 10 −8 . A quadratic convergence of iteration procedure has been observed in all test problems at all time steps.

Free flight of a beam
This example, introduced by Simo and Vu-Quoc [13], shows the ability of formulation to consider very large displacements and rotations. This is also an excellent test for the long-term stability of time integrator. Simo et al. [40] shown a rapid growth of total mechanical energy causing a failure of the iteration process when using classical time-integration schemes.
Initially inclined beam is subjected at the lower end to triangular pulse load, consisting of force f X and moments h Y and h Z , see Fig. 3. The load vanishes at t = 5 leaving the beam to fly in free motion. The remaining data of the beam are: In our simulation, the finite-element mesh consisted of 10 quadratic elements, N = 3. To demonstrate the performance of the proposed approach, we study the long-term behavior of the beam. We show the results for constant time step h = 0.1 until the time t = 1000, but need to stress that our solver experienced no difficulties and we could continue with calculation. In Figs. 4 and 5, projection of the deformed shapes on the coordinate planes XY and X Z, respectively, is depicted at the beginning and at the end of the calculation period. The problem is nonlinear and nonperiodic. Still, we can observe a paddling-like pattern, which is preserved during the whole analysis. The energy remains constant after the load is removed as illustrated in Fig. 6.

Right-angle cantilever
This example, also presented by Simo and Vu-Quoc [13], studies a right-angle cantilever beam under triangular pulse load in the direction out of the plane of the beam, see Fig. 7. Fig. 7 The right-angle cantilever subjected to out-of-plane loading Both parts of the structure are modeled with only four cubic elements, N = 4. It has been reported by several authors, see, e.g., [21], that for this example numerical problems appear at approximately t = 50. We therefore used longer time interval [0, 100] and relatively large time step h = 0.2. For comparison reasons, we also present the results for a smaller time step h = 0.02. We can observe from Fig. 8 that the results for longer time step agree well with the ones obtained by smaller time increments. Moreover, the calculations remain stable during the whole time interval.
In Fig. 9, we compare our results with the ones presented by Simo and Vu-Quoc [13]. A good agreement  between the results can be observed. Moreover, we can see that the present approach gives excellent results even if a mesh of only two elements is used (one per each leg of the cantilever).

Circular beam
Another problem that reveals the energy growth when using classical trapezoidal or midpoint rule is a closed circular beam, [40]. This problem is interesting as such structures appear in many engineering problems. A ring of radius 5 is free in space and set in motion by two forces of equal magnitude in opposite directions at points A and B, see Fig. 10. Other data are as follows: [10,10,10] .
We used a mesh of sixteen straight elements of third order, N = 4. A solution was sought on a long time Fig. 11 The displacement of the node initially at the top of the circular beam interval [0, 500] with time step h = 0.1. The calculations remain stable for the whole simulation. Besides the constant energy after the load was removed, the long-term stability of the proposed formulation is evident from time histories of displacements at point A, see Fig. 11. Figure 11 also reveals large magnitudes of displacements and a complicated response of structure without any evident pattern of movement.

Four-bar mechanism
The four-bar mechanism with misaligned joint introduced by Bauchau [55] belongs to the well-documented benchmark examples [56,57], suitable for testing new formulations. The mechanism consists of three bars connected by revolute joints, as shown in Fig. 12. The axes of rotation of revolute joints A, B and D are orthogonal to the plane of mechanism. The revolute joint at point C is inclined with respect to the normal axis by the angle ϕ = 5 • simulating initial imperfection. Bars AB and BC share equal properties: A constant angular velocity = 0.6 is applied at support A. Initial imperfection results in three-dimensional motion of the system.
We modeled each beam with four elements of third order, N = 4. For the perfect joint, Boolean identification of degrees of freedom was used, while for the inclined joint the method of Lagrange multipliers was applied. In accord with the present model, we enforced the constraints at the level of angular velocities by satisfying the equations q * e • I − II • q e · 1 0 0 T = 0 q * e • I − II • q e · 0 1 0 T = 0, where I and II denote the angular velocities of the beams BC and C D at point C and q e is the rotational quaternion describing the initial imperfection of the joint. The simulation was carried out until t = 12. Following [57], we obtained the results presented here using numerical dumping and constant time step h = 0.004. The results for the displacements of joint C are presented in Fig. 13, while Fig. 14 shows the first Euler angle of the bar BC at point C and the relative rotations at the revolute joint D. The present results are in excellent agreement with other authors, see [56,57].

Conclusions
We have presented a novel energy-preserving scheme based on velocities and angular velocities, where we satisfy the kinematic compatibility equations with the same accuracy as the governing equations. The discrete kinematic compatibility equations enforce the admissible update of strain vectors. This update is in our approach completely harmonized with the energy conservation demands, which additionally provides the consistent approximation of stress resultants. The finite element proposed is based on the interpolation of velocities in a fixed basis and angular velocities in the local basis. The use of standard shape functions is therefore completely consistent with properties of the configuration space. The main advantage of the proposed solution method is in its simplicity, robustness and longterm stability. Its favorable behavior was demonstrated by several numerical examples.