Abstract
The traditional Kinetic Dynamic Relaxation (KDR) method can only deal with simple constraints such as the fixed joints, which restricts its practical applications. This study has proposed a new algorithm to implement the complex position constraints in KDR. The 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. The numerical results show the algorithm is adaptive in both linear and nonlinear problems and works efficiently.
1. Introduction
The 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–5], especially in structural mechanics such as plate analysis [6–8], frame and truss analysis [9], cable membrane structures [10–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.
The 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). The viscous damping technique includes a fictitious damping matrix which is closer to the real dynamic system than the KDR. The most rapid convergence is obtained by damping the lowest mode of vibration. Several important methods [24–26] belong to this technique. The 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. The 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].
The 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. This 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. They cannot be implemented in the traditional KDR for the missing corresponding constraint forces.
This 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. The final conclusions are made in Section 5.
2. A Short Summary of KDR
For a system with m-degree of freedom , the governing equation describing Newton’s 2nd law isin which is the diagonal fictitious mass matrix; is the velocity vector; is the residual force vector implying , where is the external force vector and is the internal force vector; and is the viscous damping matrix and the elements of which are set to zeros in the KDR.
The 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 is time step.
Having the updated geometry, the new residual forces can be calculated. The 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 , an energy peak is detected. Topping and Ivanyi [11] assumed that the peak point of kinetic energy occurs at and hence the new positions are set to :
The analysis can be restarted from as the initial position. Meanwhile, all current velocities are reset to zero: . For the first iteration or restarting the process after the peak, the velocities at time are set aswhich is obtained from the relationwhere should be evaluated from the position of (6). This 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.
3. Algorithm Implementing Complex Position Constraints in KDR
An m-degree freedom system contains n complex position constraints:
The 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 is the Jacobian matrix with elements . Based on Lagrange’s equations of the first kind, the governing equations of a dynamic system of m-degree of freedom can be written aswhere is the corresponding constraint force vector. The 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 remains unknown and needs to be solved. The inspiration is from the half-explicit method for DAEs, which means solving and 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 gives
Then, the Taylor expansions of (13) at time t leads to the constraints with acceleration :where is the derivative of the Jacobian matrix with respect of , with the components of .
Obviously, solving equation (14) combined with Newton’s law (12) can obtain the acceleration and constraint forces. However, in (14) is too complicated for most cases. Dimension analysis technique is employed here to simplify the equation. The typical final deformation of a node is considered as the length scale, and the external force of a node as the force scale. There are two time scales which are the time step and the typical time to reach equilibrium with the relation , in which . The value of can be chosen by experience. Usually the larger the deformation , the smaller will be chosen to ensure more time steps are employed. The velocity scale can be set as . Supposing that the typical velocity can be reached in the time scale of several from static status, based on the theorem of momentum, the mass scale can be written aswhich also indicates that the acceleration has the scale of . If the system contains the elastic component, based on the balance equation in the equilibrium state, the stiffness scale is . Accordingly, the right-hand side term of the stability condition has the scale of . 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 is , 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 in (18). By solving out the constraint force at time t, the new residual forces including the effect of constraints isand 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) givesand by inserting (12) into (21), the constraint forces at time 0 can be solved out through
Thus, the initial updated residual forces can be given out by (19). The velocity update relation is achieved by following the expression similar with (7):
For most linear MPCs, , and is independent of , and (17) can describe the position constraints 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]. The Taylor expansion of at leads towhere is the drift-off value and is the residual of original position constraints. Another equation from theorem of momentum read
Combining (24) and (25), can be solved out by
The new position after projection will be
By inserting (27) into , the new residuals of original position constraints are calculated. The projection process is solved iteratively, until the norm of residuals satisfies the tolerance. LU factoring of is taken only in the first iteration to reduce the computational cost.
The 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 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 and kinetic energy 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 at time t to obtain the unbalanced forces . 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 and solve the constraint forces according to (18) (if “ini” is 1, using (22) instead).(4)Calculate unbalanced forces 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 (), calculate the initial position of restarted analysis by (6), set the previous kinetic energy to 0 and reset “ini” to 1.(7)For nonlinear constraints, project onto the geometric constraint equations by iterations with equations (26) and (27).(8)If “ini” is 1 and the kinetic energy is less than tolerance, jump out of time iterations, otherwise continue the time loop.(v)Completion
4. Algorithm Verification and Discussion
Several static problems are chosen to verify the proposed algorithm. The first and the second are linear problems with MPCs. The 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.
4.1. Linear Problems with MPCs
The 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 is considered to be unknown. This problem can be treated as linear because of the small displacements of nodes with two MPCs, including and . 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. This will bring failure in estimation of the fifth nodes mass by equation (1). Alternatively, the estimations of all node masses are from equation (15), in which , , and . The exact solution of the problem is . Figure 2 shows the kinetic energy, relative error of displacement , and relative unbalanced force converge with the increase of the time step. It takes 49 time steps to reach kinetic energy tolerance. The dropping down of energy peaks, the exponential decrease of the displacement error and unbalanced forces can also be observed in the figure.


The second example is about uniformly stressed plate containing a circular hole. The 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. The 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 and (u and are horizontal and vertical displacements). The 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 , as shown in Figure 4(a). From Figure 4(b), the well-known distribution of stress is obtained, and the maximum pressure is approximately 78 Pa, about 3 times of the load.


(a)

(b)
4.2. Geometry Nonlinear Problems with Nonlinear Position Constraints
The 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. The size and load condition is from [12]. The main span L is taken as 1000 m and the sag h is set to 90 m. The lengths of cable are calculated as S = 1021.1984 m for parabola and S = 1021.2831 m for catenary. The self-weight of the deck is taken as q = 72.4 kN/m and cable weight is mg = 26 kN/m. The half cable is divided into N bar elements with equal length . Relation between elements and nodes is shown in Figure 5. The total external force for the ith element is for parabola and for the catenary case, where 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 . Then, the nonzero elements of the sparse matrix are setting asthe elements of which belong to the ith row and the columns of to . The node mass is estimated from equation (15), in which , , and . 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) and 6(b). Figures 6(c) and 6(d) show both the kinetic energy and relative unbalanced force converge. But the relative displacement errors do not go down after about 300 steps for both cases because these inevitable errors are from the space discretization.


(a)

(b)

(c)

(d)
The 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)). The numerical model (Figure 7(c)) shows that the beam is divided into N bars with equal length in the plane xy. presents the position vector of the ith node. The edge vector is defined as , and the curvature at the ith node is [32]in which is along the z-axis and the normal of the plane. Thus, . and point along the principal direction of the section. Based on the geometry relation, the variation of the curvature can be given asin which is the identity matrix. This relation corresponds to the continuous form of [33]where is the tangent vector of the curve. Around each node, an element is defined between the centers of the neighbouring bars. The internal virtual work of element i suggests the internal force vector of the element isin which E is Young’s modulus and I is the moment of inertia. The incensement of the internal virtual work gives out the material stiffness of the element:This stiffness is only used to estimate the mass matrix. The nonlinear constraints (28) should be satisfied for the invariable length condition. And the nonzero elements of have the same expressions with (29).

(a)

(b)

(c)
In both cases about beam, N is set to 50. The 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. The parameters for the cantilever beam subjected to vertical load P are the length L = 1, EI = 4/π2, and . And the analytical solution is [34]in which and . is the complete elliptic integral of the first kind and and are the elliptic integrals of the first and second kind separately. As shown in Figure 8(b), the final beam configuration is close to the analytical solution from the initial state of the straight line. Figures 8(c) and 8(d) show the kinetic energy and relative unbalanced force converge for both cases. The dashed lines in Figures 8(c) and 8(d) show the distance between two endpoints and the error of maximum deflection comparing the analytic solution, respectively. Both of them do not go down further after certain steps because of the errors from space discretization.

(a)

(b)

(c)

(d)
4.3. Linear and Nonlinear Analysis of Cyclically Repeated Structure with MPCs
For cyclically repeated space structures under the cyclical load, if only one repeated piece is modelled instead of the whole structure, the computing costs will be reduced drastically. In this situation, complex 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. The structure is generated by cyclic replication of 24 similar pieces. Figure 9(b) shows one piece with 25 elements (continuous lines) covering 15° of the truss. The details of nodal numbers and coordinates are also illustrated in the figure. The cross-sectional area and the modulus of elasticity of members are 2 × 10−3 m2 and 200 GPa, respectively. All the underside nodes (z = 0.0) are pinged. 35 kN load is applied to all other free nodes of the dome in negative z-axis.

(a)

(b)
The 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 ofwhere 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 . Figure 10(a) shows the numerical results (the dashed line and the continuous line with circle) of the one-piece model with which MPCs coincide with the results from the literature very well. Both the kinetic energy and relative unbalanced force vs. the time steps of linear analysis have been plotted in Figure10(b), which verified the convergence of our model.

(a)

(b)
4.4. 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. This 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) but explicitly with respect to velocity and positions. Because the dimension of in (18) is , the numerical costs of KDR with constraints depend on the number of constraints n other than degrees of freedom m. Usually , 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.
5. 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. The 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
The data of the models and algorithms used to support the findings of this study are included within the article.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
This work was supported by the National Natural Science Foundation of China (nos. 11602102 and 51809135), Major Research Grant from both the National Natural Science Foundation of China and the Provincial Natural Science Foundation of Shandong (no. U1806227), and PhD Programs Foundation of Ludong University (no. LY2012017).