Development of geometrically nonlinear finite element «2-D prismatic joint» using the augmented Lagrangian method

An algorithm for obtaining the stiffness matrix and the vector of elastic forces of a finite element (FE) of kinematic constraints using the example of the geometrically nonlinear element “Two-dimensional prismatic coupling” is proposed in the article. The algorithm is based on the application of the extended method of Lagrange multipliers. The proposed approach is quite general, a wide class of finite elements of kinematic constraints can be developed with it. Also, method for taking into account friction between joined planes for the “Two-dimensional prismatic joint” FE is proposed. Coulomb friction was considered as a model of friction, however, the proposed approach can be extended to other friction models without significant changes. The developed element with friction can be used, for example, for the dynamic skew modelling in rotors that occurs during operation when friction forces are exceeded.


1.
Introduction. The development of finite elements for kinematic joints modeling is an important aspect of creating complexes for modeling the dynamics of system of many bodies. The elements of such joints can often make movements with large displacements and rotations. From a mathematical point of view, such joints are algebraic constraints, which can be formulated explicitly in rare cases. In the case of an explicit formulation, the constraints can be satisfied by eliminating the degrees of freedom of the dependent nodes directly, but in most cases the constraints are formulated implicitly, and in such case methods for solving problems with constraints should be applied [1,2]. The extended Lagrange multiplier method [3] is used as a method for solving constrained problems in this paper. The tangential stiffness matrix and the vector of nodal forces are obtained in an analytical form using this method for the finite element "Two-dimensional prismatic joint". The Euler vector is used to describe the finite rotations. Also, a method of accounting Coulomb friction between joined planes is proposed. The proposed method of friction accounting is quite universal and allows to apply various models of friction without significant modifications. The element considered in this paper can be used, for example, for the dynamic misalignment simulation in rotors that can occur during operation if external forces exceed the friction forces. The model task demonstrates this effect.

The stiffness matrix and the vector of elastic forces of the element
This connection models the conjunction between two nodes «a» and «b», each of which is movable in the general case (Ошибка! Источник ссылки не найден.). Each node has 6 degrees of freedom (3 displacements -u, 3 Euler rotation vector projections -ϑ). Another 6 variables are the Lagrange multipliers (λ), which mean the forces and moments. It is convenient to combine all unknowns into one state vector: Each node has its own coordinate system. The axes of the nodes coordinate systems coincide at the beginning; therefore, the unit vectors of both coordinate systems will be denoted as . After some element movement in the space, the axis of the coordinate system will take new position (ex, ey, ez). In each node, the Euler vectors ϑa and ϑb, respectively, which describe finite rotations, are known: , where -rotation tensors of nodes «a» and «b», respectively.
Let's make the conditions of connection between nodes. The element has 2 internal degrees of freedom, therefore, it is necessary to obtain 4 connection equations.
Three connection equations can be written on the condition that relative rotation is prohibited. This restriction in accordance with [4,5] can be represented as: The relative displacement vector uba can be calculated from nodal displacements, because the nodes coincide in the undeformed position: q q q e e e , ,  (1) allows to write the 4th connection equation (the relative displacement vector projection on the normal to the connected planes): The element potential using the extended method of Lagrange multipliers will be represented by: (2) where ϑba, fba -movement and rotation residuals; l 1, l 2 -Lagrange multipliers; p1, p2 -penalty factors.
The analytical expressions for the vector of nodal forces and stiffness matrix can be obtained by varying equation (2). The increments of the components that make up the element's potential will look like: where B(ϑ) -P. A. Zhilin tensor, L(ϑ) -rotation tensor (see [4,5]).
An expression for the vector of nodal forces Pe can be obtained from the first variation of the potential: The expression for the tangent stiffness matrix [K] follows from the second variation of the potential: ; ;

Friction accounting
It should be noted that the subsequent calculations are general, if it is necessary, an arbitrary law of friction may be introduced and a tangent damping matrix for it may be obtained similarly. In this paper, for simplicity, we restrict ourselves to Coulomb friction.
Let's assume, that Coulomb friction, which is described by a continuous function of the relative velocity, acts between the planes that are connected by the prismatic connection [6]: where "# -relative velocity vector, μ -friction coefficient, N -normal force, v % -friction model sensitivity parameter. Equation (3) of the friction force can be introduced into the vector of nodal forces as follows taking into account the kinematic restrictions imposed on the element: The tangent damping matrix is formed on the basis of the first variation equation (3) taken with the opposite sign. Omitting simple but cumbersome transformations, a variation of equation (3)

Test cases
The test task was solved (figure 2) to verify the correctness of the relations describing the friction force. Nodes 1 and 2 are connected by the developed element, friction is taken into account between the joint planes. The friction coefficient μ = 0.1, the normal force N = 100 N, the sensitivity parameter of the friction model was v0 = 10 -11 ms -1 . The Newmark method was used for motion equations integration. The law of F force changing is presented below: It follows from the F(t) time dependence that sliding should occur at time of 10 s, and sticking must occur at time of 12 s. In this test, the system mass tends to zero, so the setting time should be at 12 seconds, in a more inertial system, setting would take place later. Figure 3 shows the moments of sliding and sticking, they correspond to the expected results, which confirms the correctness of the introduced friction forces vector and the obtained tangent damping matrix.
, , As it was mentioned earlier, the kinematic joint considered in this paper can be used to model the dynamic misalignment occurring in rotors. To demonstrate this effect, the problem is considered (figure 4): 500 mm shaft length, R = 20 mm outer radius, r = 10 mm inner radius, 10 kg point mass, 40 gsm unbalance, μ = 0.1 friction coefficient, N = 100 N normal force, the sensitivity parameter of the friction model was v0 = 10 -9 ms -1 . The rotation speed ω changing law is presented below: Figure 4. The scheme of the test task №2 Figure 5 shows a graph of the "Two-dimensional prismatic joint with friction" nodes relative displacement. An analysis of Figure 5 shows that at a certain rotation speed (approximately 3500 rpm), the friction forces are exceeded, and the nodes begin to slide. In the interval from 1 to 2 seconds, when the rotation speed is constant, the rotor precesses in an offset position, then braking starts and at a speed of about 3400 rpm sticking occurs. It is important to note that there remains some distortion, which actually violates the initial balancing of the system, by the time the initial coaxial rotor stops completely. The subsequent acceleration of such a system will have a completely different picture of vibrations.   Figure 5. The time dependence of the "Two-dimensional prismatic joint with friction" element nodes relative displacement It should be noted that the considered problem is a model one, it is presented solely to demonstrate the dynamic misalignment and the effects that it can lead to. In a real task, the maximum value of the relative displacement is usually limited by the gaps. One may account these gaps by different ways: • direct introduction of the gap into the program functional in the form of an inequality type constraint; • parallel setting of the limiter finite element.
The second method is a bit simpler and allows to use elements independently if it is necessary.

Conclusion
This paper presents analytical relationships of the tangent stiffness matrix and the nodal force vector for the "Two-dimensional prismatic joint" FE. To describe the finite rotations of the nodes, the Euler vector was used, which allows to extend the elements for large displacements and rotations accounting. The algorithm for obtaining the stiffness matrix and force vector is based on the extended method of Lagrange multipliers. The proposed approach is quite general, which allows to use it in the development of a wide class of FE kinematic constraints ("Hinge", "Rigid rod", etc.). Also, a method for taking into account friction between planes connected by the considered coupling is proposed. Coulomb friction is considered as an example, however, relations for another friction law can be obtained in a similar way, taking into account, for example, the Stribeck effect [7,8]. The correctness of the relations describing friction is confirmed by a test example. The considered finite element with friction, can be used, for example, to simulate the rotor skew, which varies during operation (when the friction force is exceeded by external loads). This effect is demonstrated on the model task. Such effects in a real design can lead to a violation of the system initial balancing, which in its turn may lead to an increased vibration level.