Implementation of Complex Position Constraints in the Kinetic Dynamic Relaxation Method

(e traditional Kinetic Dynamic Relaxation (KDR) method can only deal with simple constraints such as the fixed joints, which restricts its practical applications. (is study has proposed a new algorithm to implement the complex position constraints in KDR. (e proposed algorithm is developed by transforming the position constraints to the acceleration form through combining time differentiation and Taylor expansion with dimensional analysis, and then solving the governing equation of the constraint forces with the Newton’s 2nd law. For the nonlinear constraints, projection technique is applied to avoid the drift-off phenomenon. Several tests have been performed to verify the proposed algorithm. (e numerical results show the algorithm is adaptive in both linear and nonlinear problems and works efficiently.


Introduction
e Dynamic Relaxation (DR) method is a way of finding the equilibrium solution for the static problem through the damping effect of the dynamic system. Since the method was proposed in the 1960s [1,2], it has been improved and applied in many fields [3][4][5], especially in structural mechanics such as plate analysis [6][7][8], frame and truss analysis [9], cable membrane structures [10][11][12], actively bent structures [13,14], and postbuckling analysis [15]. Recently, Rezaiee-Pajand et al. [16] made a good summary in the structural applications of DR. e most important feature of DR is that the time steps, the mass, and the damping terms are virtual [17], for only the final state of the dynamic system is concerned other than the motion. And also it is an explicit iterative technique because it does not solve linear equations or even assemble the stiff matrix [18], which are time and memory consuming. Hence, the method is suitable for highly nonlinear problems. Usually, the fictitious time step is set as a constant. However, some researchers also developed several ways to obtain better time steps, including the residual force or residual energy minimizing [19,20], zero damping [21], and error estimate [22]. For fictitious masses, Alamatian [23] proposed a general mass formulation for the Kinetic DR method. Based on whether there exists a damping matrix, there are two kinds of damping techniques used: Viscous Damping Relaxation and Kinetic Damping Relaxation (KDR). e viscous damping technique includes a fictitious damping matrix which is closer to the real dynamic system than the KDR. e most rapid convergence is obtained by damping the lowest mode of vibration. Several important methods [24][25][26] belong to this technique. e kinetic damping technique was firstly introduced by Cundall [3] for application to unstable rock mechanism. In the technique, only the time step and the fictitious nodal masses are required other than viscous damping term. When the peak of an undamped system kinetic energy is traced, the potential energy is minimized simultaneously, and then all velocity components are set to zeros. e process is restarted and repeated through further peaks until the most kinetic energy has been dissipated. High stability and rapid convergence make it widely used. It is reported that the KDR has good performance in the nonlinear analysis of frame and cable structures [9,27]. Also, KDR can be paralleled to accelerate the large problem in form finding [11,28]. e existing DR techniques always concentrate on the improvement of convergence speed, applications to various structures, or the key problems in the nonlinear analysis. However, researchers rarely focus on developing general techniques to deal with various complex position constraints, which are very common in the structure analysis. For instance, the traditional KDR method only implements simple constraints such as the fixed joints by assigning large masses or setting the zero residual forces [10,11], which is not enough in the structure modelling. is research mainly focuses on how to extend the KDR method to handle more complex constraints.
In the applications, many structures contain very stiff components so that not all the stiffness of degrees of freedom is in the same order. For example, the stiffness of stretching of the long thin beam is much larger than bending. According to the stability condition of KDR [17] the fictitious mass should be set in the same order with the largest stretching stiffness which led to extremely large masses. Apparently, this will slow down the convergence of KDR if the beam is under significant bending but neglectable stretching. If any deformations are not significant, position constraints should be introduced to accelerate the convergence. Common linear constraint is the Multipoint Constraint (MPC), which can be introduced for many situations such as the structures containing rigid components, the interface transition between high and low order elements, and implementation of cyclic symmetry boundary. Some inextensible problems such as the deformation of rods or fabrics can introduce nonlinear constraints [29]. All the mentioned constraints above are described as the complex position constraints. ey cannot be implemented in the traditional KDR for the missing corresponding constraint forces.
is study proposes a new algorithm dealing with the complex position constraints, which can be in the form of linear or nonlinear. As this method is based on the KDR, a brief method description is offered in Section 2. In Section 3, it shows how to obtain the constraint forces in detail and a modified KDR with constraints is followed. Section 4 verifies the proposed algorithm by several models including MPCs and nonlinear constraints. e final conclusions are made in Section 5.

A Short Summary of KDR
For a system with m-degree of freedom x � x 1 x 2 . . . x m T , the governing equation describing Newton's 2nd law is where f ext is the external force vector and f int is the internal force vector; and C � diag C 1 C 2 . . . C m is the viscous damping matrix and the elements of which are set to zeros in the KDR. e acceleration is assumed constant over the time step, which induces the following velocity iterative relations: and the position iterative relations: where the superscripts indicate the time point and Δt is time step.
Having the updated geometry, the new residual forces can be calculated. e total kinetic energy of the system is traced during iterations, which is calculated by the relationship When the current kinetic energy is less than the previous one at t − (Δt/2), an energy peak is detected. Topping and Ivanyi [11] assumed that the peak point of kinetic energy occurs at t − (Δt/2) and hence the new positions are set to e analysis can be restarted from * x 0 as the initial position. Meanwhile, all current velocities are reset to zero: u 0 � 0. For the first iteration or restarting the process after the peak, the velocities at time Δt/2 are set as which is obtained from the relation where R 0 should be evaluated from the position of (6). is process is iterated through several peaks until the convergence criterions are satisfied. In practice, the time step is set to 1 and the elements of the mass matrix are determined by equation (1) to ensure the stability.

Algorithm Implementing Complex Position Constraints in KDR
An m-degree freedom system contains n complex position constraints: e constraints can be linear or nonlinear. Taking time derivative of the position constraints (9) and according to the chain rule, the velocity constrains are obtained: in which G(x) is the Jacobian matrix with elements G ij � zg i /zx j . Based on Lagrange's equations of the first kind, the governing equations of a dynamic system of mdegree of freedom can be written as _ x � u, where λ � λ 1 λ 2 . . . λ n T is the corresponding constraint force vector.
e system is called Differential-Algebraic Equations (DAEs) of index 2 [30]. If the velocity constraints of (11) are replaced with position constraint (9), index 3 DAEs are obtained. From (9) to (10), it is the index reduction of DAEs by time differentiation.
Newton's 2nd law (the second equation in (11)) gives the acceleration at time t: However, the corresponding constraint force vector λ t remains unknown and needs to be solved. e inspiration is from the half-explicit method for DAEs, which means solving x and u in an explicit manner and the algebraic variables λ in an implicit manner [30]. By inserting the velocity and position iterative relations (3) and (4), the velocity constraint (10) at time t + Δt/2 gives en, the Taylor expansions of (13) at time t leads to the constraints with acceleration _ u t : where G x is the derivative of the Jacobian matrix with respect of x, with the components of (G x ) ijk � zG ij /zx k . Obviously, solving equation (14) combined with Newton's law (12) can obtain the acceleration and constraint forces. However, G x in (14) is too complicated for most cases. Dimension analysis technique is employed here to simplify the equation. e typical final deformation of a node l 0 is considered as the length scale, and the external force of a node F as the force scale. ere are two time scales which are the time step Δt and the typical time t 0 to reach equilibrium with the relation Δt � εt 0 , in which ε ≪ 1. e value of ε can be chosen by experience. Usually the larger the deformation l 0 , the smaller ε will be chosen to ensure more time steps are employed. e velocity scale can be set as u 0 � l 0 /t 0 . Supposing that the typical velocity u 0 can be reached in the time scale of several Δts from static status, based on the theorem of momentum, the mass scale can be written as which also indicates that the acceleration has the scale of _ u ∼ u 0 /(t 0 ε). If the system contains the elastic component, based on the balance equation in the equilibrium state, the stiffness scale is S ∼ F/l 0 . Accordingly, the right-hand side term of the stability condition has the scale of (Δt 2 /2) j |S ij | ∼ ε(FΔt/u 0 ) ∼ εM. It is obvious that the stability condition (1) is satisfied. Next, the scale of equation (14) is estimated. Supposing the scale of 2-norm of G is ‖G‖ ∼ G 0 , then other terms have the following scales: Zero order of ε for (14) gives out the approximate constraints for acceleration, which also have a simple form: Inserting Newton's 2nd law at time t (12) into (17) gives the equation for constraint forces: It is very efficient to get the inverse of diagonal mass matrices M − 1 in (18). By solving out the constraint force at time t, the new residual forces including the effect of constraints is and the new velocities iterative relation turns out to be When implementing the constraints at the first iteration or restarting the process after the peak, inserting (8) into constraints for acceleration (16) gives and by inserting (12) into (21), the constraint forces at time 0 can be solved out through us, the initial updated residual forces R 0 can be given out by (19). e velocity update relation is achieved by following the expression similar with (7): For most linear MPCs, G · x � b, G and b is independent of x, and (17) can describe the position constraints G · Δx � 0 exactly. However, for nonlinear problems, original position constraint (9) needs to be satisfied too. A position projection technique is employed after one or several position updates to prevent the drift-off phenomenon [30]. e Taylor expansion of g(x t + δx) � 0 at x t leads to where δx is the drift-off value and g(x t ) is the residual of original position constraints. Another equation from theorem of momentum read Mathematical Problems in Engineering Combining (24) and (25), μ can be solved out by e new position after projection will be By inserting (27) into g, the new residuals of original position constraints are calculated. e projection process is solved iteratively, until the norm of residuals satisfies the tolerance. LU factoring of GM − 1 G T is taken only in the first iteration to reduce the computational cost. e iterative procedure of the modified KDR is summarized as follows: (i) Set physical properties and number of nodes. (ii) Set the initial position of the structure. Set time step Δt to 1 and convergence criterion for the kinetic energy, and Tol � 1.0 × 10 − 12 is used for this study. (iii) Set initial fictitious velocities u and kinetic energy U k to 0. Set the first time step flag "ini" to 1. (iv) Time iterations: (1) Calculate the internal forces and the external forces related to x at time t to obtain the unbalanced forces R t � f ext − f int . If possible, the stiffness matrix is given here. (2) Estimate the mass matrix according to the stability condition from stiffness (1) or the mass scale (15). (3) Supply the constraint derivative matrix G and solve the constraint forces λ t according to (18) (if "ini" is 1, using (22) instead). (4) Calculate unbalanced forces R t according to (19). (5) Update the node speed using (20) (if "ini" is 1, using (23) instead and reset "ini" to 0). Update the node coordinate with (4). (6) If the kinetic energy (5) is less than the previous (t − (Δt/2)), calculate the initial position of restarted analysis by (6), set the previous kinetic energy to 0 and reset "ini" to 1. (v) Completion

Algorithm Verification and Discussion
Several static problems are chosen to verify the proposed algorithm. e first and the second are linear problems with MPCs. e third and the fourth ones are the geometry nonlinear problems with nonlinear constraints. Finally, a complex cyclically repeated structure with MPCs is solved for both linear and nonlinear cases.

Linear Problems with MPCs.
e first example is a linear problem with MPCs. A rigid bar of negligible mass with one end pinned is hanged by two rods with different physical properties, as shown in Figure 1 (this example is from [31]). And load P is applied on the other end. Two bar elements (number with circle) and five nodes are used to model the structure. Only the vertical displacement vector x � [x 1 , x 2 , x 3 , x 4 , x 5 ] T is considered to be unknown. is problem can be treated as linear because of the small displacements of nodes with two MPCs, including 3x 1 − x 5 � 0 and 6x 2 − 5x 5 � 0. It can be found that the fifth node has no elastic connection with other nodes except the constraints, which induced that the elements of the 5th row and column of the global stiffness matrix are all zeros. is will bring failure in estimation of the fifth nodes mass by equation (1). Alternatively, the estimations of all node masses are M i ∼ 1.5 × 10 8 kg from equation (15) Figure 2 shows the kinetic energy, relative error of displacement ‖x − x exact ‖/‖x exact ‖, and relative unbalanced force ‖R‖/P converge with the increase of the time step. It takes 49 time steps to reach kinetic energy tolerance. e dropping down of energy peaks, the exponential decrease of the displacement error and unbalanced forces can also be observed in the figure.
e second example is about uniformly stressed plate containing a circular hole. e geometry, the finite element model, and the loading and boundary conditions are shown in Figure 3. One-quarter of the square plate with the center hole is all needed to solve the problem for the symmetry of geometry and loading. Uniform pressure 25 Pa is acting parallel to the x-axis. e thick of the plate is 1 m. Young's modulus and Poisson's ratios are 20 GPa and 0.2, respectively. Hence, the problem is in plane stress and also a linear problem for the small displacement. Because of stress concentration, the stress gradient near the hole can be very high which needs high order elements such as eight node quadrilateral elements. But far away from the hole, stress distribution is flatter and four node quadrilateral elements are enough. One way of transition between two kinds of element is by using MPCs for maintaining compatibility. As shown in Figure 3, two linear MPCs are required between elements A and B for nodes i, j, and k, which are u i − 2u j + u k � 0 and v i − 2v j + v k � 0 (u and v are horizontal and vertical displacements). e mass is estimated by equation (1) because the stiffness matrix can be obtained by traditional finite element procedure. In the calculation by KDR with MPCs, it takes 451 time iterations to reach kinetic energy tolerance; meanwhile, the norm of relative unbalanced forces endure a continuous drop down and reach to ‖R‖/‖f ext ‖ � 8.09 × 10 − 5 , as shown in Figure 4(a). From Figure 4(b), the well-known distribution of stress σ x is obtained, and the maximum pressure is approximately 78 Pa, about 3 times of the load.

Geometry Nonlinear Problems with Nonlinear Position Constraints.
e third example involves the form finding of an inextensible suspension cable which is important to suspension bridges. Two forms of cable are studied here: parabola when the weightless cable subjected to a uniformly distributed deck weight load and catenary for a heavy cable with neglectable deck weight. e size and load condition is from [12]. e main span L is taken as 1000 m and the sag h is set to 90 m. e lengths of cable are calculated as S = 1021.1984 m for parabola and S = 1021.2831 m for catenary. e self-weight of the deck is taken as q = 72.4 kN/m and cable weight is mg = 26 kN/m. e half cable is divided into N bar elements with equal length Δl � S/N. Relation between elements and nodes is shown in Figure 5. e total external force for the ith element is qΔx i for parabola and mgΔl for the catenary case, where Δx i � x i+1 − x i is the x coordinate difference between two neighbouring nodes. Instead of using extremely large tension stiffness, zero stiffness is employed here which leads to zero internal force. Alternatively, the inextensible condition for the ith bar is guaranteed by the nonlinear constraints: where Δy i � y i+1 − y i . en, the nonzero elements of the sparse matrix G are setting as the elements of which belong to the ith row and the columns of 2i − 1 to 2i + 2. e node mass M i is estimated from equation (15), in which F ∼ ‖f ext ‖, ε ∼ 0.01, and l 0 ∼ Δl. N is set to 50 to simulate half span. Results are shown in Figure 6. It takes 752 and 603 steps to reach kinetic energy tolerance for parabola and catenary from a straight line, and the numerical results fit the exact curves very well, as shown in Figures 6(a)   Mathematical Problems in Engineering and 6(d) show both the kinetic energy and relative unbalanced force ‖R‖/‖f ext ‖ converge. But the relative displacement errors ‖x − x exact ‖/‖x exact ‖ do not go down after about 300 steps for both cases because these inevitable errors are from the space discretization.
e fourth example involves inextensible and incompressible long thin vertical cantilever beams subjected to moment M on the tip (Figure 7(a)) and vertical load P which causes large postbulking displacement (Figure 7(b)). e numerical model (Figure 7(c)) shows that the beam is  Mathematical Problems in Engineering divided into N bars with equal length Δl in the plane xy. X i � x i y i T presents the position vector of the ith node. e edge vector is defined as t i � (X i+1 − X i )/Δl, and the curvature at the ith node is [32] in which d 1 is along the z-axis and the normal of the plane. us, d i 2 � − t i 2 t i 1 T . d 1 and d 2 point along the principal direction of the section. Based on the geometry relation, the variation of the curvature can be given as in which I is the identity matrix. is relation corresponds to the continuous form of [33] δκ where t is the tangent vector of the curve. Around each node, an element is defined between the centers of the neighbouring bars. e internal virtual work δκ i EIκ i Δl of element i suggests the internal force vector of the element is in which E is Young's modulus and I is the moment of inertia. e incensement of the internal virtual work gives out the material stiffness of the element: is stiffness is only used to estimate the mass matrix. e nonlinear constraints (28) should be satisfied for the invariable length condition. And the nonzero elements of G have the same expressions with (29).
In both cases about beam, N is set to 50. e parameters for the cantilever beam subjected to moment M are the length L � 2, EI � 100, and M � 100π. As shown in Figure 8(a), the final calculated beam configuration is almost a perfect circle which is the analytical solution from the initial state of the straight line. e parameters for the cantilever beam subjected to vertical load P are the length L � 1, EI � 4/π 2 , and P � (4K 2 (k)/π 2 ) ≈ 4.0301. And the analytical solution is [34] x � 2k(cos ϕ − 1) in which k � sin(4π/9) and ϕ ∈ [0, 2π]. K(k) is the complete elliptic integral of the first kind and F(ϕ, k) and E(ϕ, k) are the elliptic integrals of the first and second kind separately. As shown in Figure 8  position constraints should be introduced for the node pairs at the cyclic symmetry boundary, which cannot be achieved by the traditional treatment of simple support. For example, a single layer circular truss dome is concerned, as shown in Figure 9(a) [23,35], which has 216 nodes and 600 elements for the whole truss. e structure is generated by cyclic replication of 24 similar pieces. Figure 9 e structure is solved by both linear small deflection and nonlinear large deformation (total Lagrangian bar elements) theories. Owing to the cyclic symmetry, just one piece of truss (17 nodes with 51 freedoms) is solved instead of solving the whole structure (216 nodes with 648 freedoms). Obviously, the cost of calculating one piece of truss with MPCs is much cheaper than the whole dome. However, the cyclic symmetry constraints for 8 node pairs must be introduced by 24 MPCs which have the form of where u, v, and w are displacements along x-, y-, and z-axis and i is the odd node number from 3 to 17 in Figure 9(b). For one piece, only node 1 is hinged, which leads to u(1) � v(1) � w(1) � 0. Figure 10

Discussion.
Similar with the existing KDRs, the algorithm proposed in this paper does not need fictitious damping matrix. And it also can deal with simple constraints such as the fixed joints, which led to the degeneration of the method to the existing KDR. is is implemented by (19) which equals to set zero residual forces. Numerical results show our method not only inherits good performance in stability and rapid convergence for both linear and nonlinear problem but also can reduce the computing costs dramatically by introducing the complex position constraints than the traditional one. Unlike the traditional full explicit KDR, the algorithm is half-explicit: solving linear equations is needed to obtain the constraint forces from (18)   in (18) is n × n, the numerical costs of KDR with constraints depend on the number of constraints n other than degrees of freedom m. Usually n ≪ m, which indicates the solving scale of the problem is small.
One limitation of the proposed algorithm is about estimating the mass matrix. For most cases, it can be determined based on the stiffness matrix. However, if it lacks stiffness for few degrees of freedom, it has to estimate the mass by using (15) which needs artificial intervention and experience to ensure less iterations. A better way to generate the proper mass matrix automatically is needed, which will be studied further.

Conclusion
In this study, the new algorithm is proposed to implement the complex position constraints in the traditional KDR methods. Based on the combination of the index reduction of DAEs and Taylor expansion with dimensional analysis, equation (18) is proposed to obtain the missing constraint forces. And drift-off avoiding projection method is further proposed for the nonlinear constraints. Stiff bar problem, stress concentration near the circular hole, inextensible hanging cables, inextensible beams with tip forces, and cyclically repeated structure problem are tested with linear MPCs or nonlinear position constraints. e results show that the newly proposed algorithm can extent the KDR method to the structures with complex position constraints and is of more modelling flexibilities.
Data Availability e data of the models and algorithms used to support the findings of this study are included within the article.

Conflicts of Interest
e authors declare that there are no conflicts of interest regarding the publication of this paper.