# Implementation of Complex Position Constraints in the Kinetic Dynamic Relaxation Method

Authors: Yan Liu; Zaijin You; Yanqing Han
Journal: Mathematical Problems in Engineering
 (2020)

---

## 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.


---

## Body

## 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](1)Mi≥Δt22∑jSij,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 Section2. 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 withm-degree of freedom x=x1x2…xmT, the governing equation describing Newton’s 2nd law is(2)Mu˙=R−Cuwithx˙=u,in which M=diagM1M2…Mm is the diagonal fictitious mass matrix; u=u1u2…umT is the velocity vector; R=R1R2…RmT is the residual force vector implying R=fext−fint, where fext is the external force vector and fint is the internal force vector; and C=diagC1C2…Cm 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:(3)ut+Δt/2=ut−Δt/2+u˙tΔt,and the position iterative relations:(4)xt+Δt=xt+ut+Δt/2Δt,where the superscripts indicate the time point and Δt 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(5)Ukt+Δt/2=ut+Δt/2TMut+Δt/22.When the current kinetic energy is less than the previous one att−Δ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 xt−Δt/2:(6)x0∗=xt−Δt/2=xt+Δt−3Δt⋅ut+Δt/22+Δt2M−1Rt2,The analysis can be restarted fromx0∗ as the initial position. Meanwhile, all current velocities are reset to zero: u0=0. For the first iteration or restarting the process after the peak, the velocities at time Δt/2 are set as(7)uΔt/2=M−1R0Δt2,which is obtained from the relation(8)u−Δt/2=−uΔt/2=−12u˙0Δt,where R0 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

Anm-degree freedom system contains n complex position constraints:(9)g=g1xg2x…gnxT=0.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:(10)Gx⋅u=0,in which Gx is the Jacobian matrix with elements Gij=∂gi/∂xj. Based on Lagrange’s equations of the first kind, the governing equations of a dynamic system of m-degree of freedom can be written as(11)x˙=u,M⋅u˙=R−GTλ,G⋅u=0,where λ=λ1λ2…λnT 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:(12)u˙t=M−1Rt−GtTλt.However, the corresponding constraint force vectorλt remains unknown and needs to be solved. The 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(13)Gxt+ut+Δt/2Δt2⋅ut−Δt/2+u˙tΔt=0.Then, the Taylor expansions of (13) at time t leads to the constraints with acceleration u˙t:(14)Gxt+Gx⋅ut−Δt/2Δt2+u˙tΔt22⋅ut−Δt/2+u˙tΔt=0,where Gx is the derivative of the Jacobian matrix with respect of x, with the components of Gxijk=∂Gij/∂xk.Obviously, solving equation (14) combined with Newton’s law (12) can obtain the acceleration and constraint forces. However, Gx 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 l0 is considered as the length scale, and the external force of a node F as the force scale. There are two time scales which are the time step Δt and the typical time t0 to reach equilibrium with the relation Δt=εt0, in which ε≪1. The value of ε can be chosen by experience. Usually the larger the deformation l0, the smaller ε will be chosen to ensure more time steps are employed. The velocity scale can be set as u0=l0/t0. Supposing that the typical velocity u0 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(15)M∼FΔtu0∼FΔt2εl0,which also indicates that the acceleration has the scale of u˙∼u0/t0ε. If the system contains the elastic component, based on the balance equation in the equilibrium state, the stiffness scale is S∼F/l0. Accordingly, the right-hand side term of the stability condition has the scale of Δt2/2∑jSij∼εFΔt/u0∼ε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∼G0, then other terms have the following scales:(16)Gx⋅ut−Δt/2Δt2∼Gx⋅u˙tΔt22∼G0ε,ut−Δt/2∼u˙tΔt∼u0.Zero order ofε for (14) gives out the approximate constraints for acceleration, which also have a simple form:(17)Gxtut+Δt/2=Gxtut−Δt/2+u˙tΔt=0.Inserting Newton’s 2nd law at timet (12) into (17) gives the equation for constraint forces:(18)GtM−1GtTλt=GtM−1Rt+1ΔtGtut−Δt/2.It is very efficient to get the inverse of diagonal mass matricesM−1 in (18). By solving out the constraint force at time t, the new residual forces including the effect of constraints is(19)R˜t=Rt−GtTλt,and the new velocities iterative relation turns out to be(20)ut+Δt/2=ut−Δt/2+M−1R˜tΔt.When implementing the constraints at the first iteration or restarting the process after the peak, inserting (8) into constraints for acceleration (16) gives(21)G0u˙0=0,and by inserting (12) into (21), the constraint forces at time 0 can be solved out through(22)G0M−1G0Tλ0=G0M−1R0.Thus, the initial updated residual forcesR˜0 can be given out by (19). The velocity update relation is achieved by following the expression similar with (7):(23)uΔt/2=M−1R˜0Δt2.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]. The Taylor expansion of gxt+δx=0 at xt leads to(24)Gxt⋅δx=−gxt,where δx is the drift-off value and gxt is the residual of original position constraints. Another equation from theorem of momentum read(25)M⋅δx+GxtTμ=0.Combining (24) and (25), μ can be solved out by(26)GM−1GTμ=gxt.The new position after projection will be(27)xnew=x+δx=x−M−1GTμ.By inserting (27) into g, 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 GM−1GT 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Δ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 velocitiesu and kinetic energy Uk to 0. Set the first time step flag “ini” to 1.(iv)
Time iterations:(1)
Calculate the internal forces and the external forces related tox at time t to obtain the unbalanced forces Rt=fext−fint. 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 matrixG and solve the constraint forces λt according to (18) (if “ini” is 1, using (22) instead).(4)
Calculate unbalanced forcesR˜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.(7)
For nonlinear constraints, projectx 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 Figure1 (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=x1,x2,x3,x4,x5T is considered to be unknown. This problem can be treated as linear because of the small displacements of nodes with two MPCs, including 3x1−x5=0 and 6x2−5x5=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. This will bring failure in estimation of the fifth nodes mass by equation (1). Alternatively, the estimations of all node masses are Mi∼1.5×108kg from equation (15), in which F∼30×103N, ε∼0.02, and l0∼1×10−2m. The exact solution of the problem is xexact=27,67.5,0,0,81T/55375m. Figure 2 shows the kinetic energy, relative error of displacement x−xexact/xexact, and relative unbalanced force R˜/P 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.Figure 1
Sketch of a structure with rigid bar.Figure 2
Validation of convergence for rigid bar problem.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 Figure3. 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 ui−2uj+uk=0 and vi−2vj+vk=0 (u and v 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 R˜/fext=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.Figure 3
Finite element modelling of stress concentration problem.Figure 4
Numerical result of stress concentration problem calculated by KDR with MPCs. (a) Validation of convergence. (b) The distribution of stressσx.
(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 Δl=S/N. Relation between elements and nodes is shown in Figure 5. The total external force for the ith element is qΔxi for parabola and mgΔl for the catenary case, where Δxi=xi+1−xi 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:(28)Δxi2+Δyi2=Δl,where Δyi=yi+1−yi. Then, the nonzero elements of the sparse matrix G are setting as(29)Gi,2i−1:2i+2=−Δxi−ΔyiΔxiΔyiΔl,the elements of which belong to the ith row and the columns of 2i−1 to 2i+2. The node mass Mi is estimated from equation (15), in which F∼fext, ε∼0.01, and l0∼Δ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) and 6(b). Figures 6(c) and 6(d) show both the kinetic energy and relative unbalanced force R˜/fext converge. But the relative displacement errors x−xexact/xexact do not go down after about 300 steps for both cases because these inevitable errors are from the space discretization.Figure 5
Numerical model of half-span cable.Figure 6
The numerical results of cables using KDR with constraints. (a) and (b) are the simulation shapes of half-span cables for parabola and catenary, respectively. (c) and (d) are the validations of convergence for parabola and catenary, respectively.
(a)
(b)
(c)
(d)The fourth example involves inextensible and incompressible long thin vertical cantilever beams subjected to momentM 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 Δl in the plane xy. Xi=xiyiT presents the position vector of the ith node. The edge vector is defined as ti=Xi+1−Xi/Δl, and the curvature at the ith node is [32](30)κi=2Δlti−1×ti⋅d11+ti−1⋅ti,in which d1 is along the z-axis and the normal of the plane. Thus, d2i=−t2it1iT. d1 and d2 point along the principal direction of the section. Based on the geometry relation, the variation of the curvature can be given as(31)δκi=BiδXieΔl,in which(32)Xie=xi−1yi−1xiyixi+1yi+1,Bi=−d2iTΔti+d2i−1TΔti−1,Δti=02×2−I2×2I2×2,Δti−1=−I2×2I2×202×2.I is the identity matrix. This relation corresponds to the continuous form of [33](33)δκ=t×δt′⋅I−t⊗t⋅d1=−δt⋅d2′,where t 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 δκiEIκiΔl of element i suggests the internal force vector of the element is(34)finte=BiTEIκi,in 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:(35)kie=EIBiTBiΔl.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 G have the same expressions with (29).Figure 7
Sketches and numerical model of the inextensible and incompressible beam. (a) A cantilever beam under tip moment. (b) A vertical cantilever beam under vertical load which causes large postbulking displacement. (c) The numerical model of a beam in the planexy.
(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 P=4K2k/π2≈4.0301. And the analytical solution is [34](36)x=2kcosϕ−1Kk,y=2Eϕ,k−Fϕ,kKk,in which k=sin4π/9 and ϕ∈0,2π. Kk 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(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 R˜/fext 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.Figure 8
The numerical results of cables using KDR with constraints. (a) and (b) are the simulation shapes of cantilever beams under tip moment and vertical force, respectively. (c) and (d) are the validations of convergence for beams under tip moment and vertical forces, respectively.
(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 Figure9(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.Figure 9
Circular truss dome. (a) The full dome. (b) Details of one repeated piece.
(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 of(37)ui−ui−1cos15∘−vi−1sin15∘=0,vi+ui−1sin15∘−vi−1cos15∘=0,wi−wi−1=0,where u,v,andw 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 u1=v1=w1=0. 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 R˜/fext vs. the time steps of linear analysis have been plotted in Figure10(b), which verified the convergence of our model.Figure 10
The numerical results of the circular truss dome. (a) The load-deflection curve at tip of the truss tower. Data shown by stars and crosses are taken from Ref. [23] while dashed line and continuous line with circle are our results. (b) Convergence of the linear case.
(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 GtM−1GtT 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.


## 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 Figure1 (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=x1,x2,x3,x4,x5T is considered to be unknown. This problem can be treated as linear because of the small displacements of nodes with two MPCs, including 3x1−x5=0 and 6x2−5x5=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. This will bring failure in estimation of the fifth nodes mass by equation (1). Alternatively, the estimations of all node masses are Mi∼1.5×108kg from equation (15), in which F∼30×103N, ε∼0.02, and l0∼1×10−2m. The exact solution of the problem is xexact=27,67.5,0,0,81T/55375m. Figure 2 shows the kinetic energy, relative error of displacement x−xexact/xexact, and relative unbalanced force R˜/P 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.Figure 1
Sketch of a structure with rigid bar.Figure 2
Validation of convergence for rigid bar problem.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 Figure3. 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 ui−2uj+uk=0 and vi−2vj+vk=0 (u and v 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 R˜/fext=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.Figure 3
Finite element modelling of stress concentration problem.Figure 4
Numerical result of stress concentration problem calculated by KDR with MPCs. (a) Validation of convergence. (b) The distribution of stressσx.
(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 Δl=S/N. Relation between elements and nodes is shown in Figure 5. The total external force for the ith element is qΔxi for parabola and mgΔl for the catenary case, where Δxi=xi+1−xi 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:(28)Δxi2+Δyi2=Δl,where Δyi=yi+1−yi. Then, the nonzero elements of the sparse matrix G are setting as(29)Gi,2i−1:2i+2=−Δxi−ΔyiΔxiΔyiΔl,the elements of which belong to the ith row and the columns of 2i−1 to 2i+2. The node mass Mi is estimated from equation (15), in which F∼fext, ε∼0.01, and l0∼Δ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) and 6(b). Figures 6(c) and 6(d) show both the kinetic energy and relative unbalanced force R˜/fext converge. But the relative displacement errors x−xexact/xexact do not go down after about 300 steps for both cases because these inevitable errors are from the space discretization.Figure 5
Numerical model of half-span cable.Figure 6
The numerical results of cables using KDR with constraints. (a) and (b) are the simulation shapes of half-span cables for parabola and catenary, respectively. (c) and (d) are the validations of convergence for parabola and catenary, respectively.
(a)
(b)
(c)
(d)The fourth example involves inextensible and incompressible long thin vertical cantilever beams subjected to momentM 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 Δl in the plane xy. Xi=xiyiT presents the position vector of the ith node. The edge vector is defined as ti=Xi+1−Xi/Δl, and the curvature at the ith node is [32](30)κi=2Δlti−1×ti⋅d11+ti−1⋅ti,in which d1 is along the z-axis and the normal of the plane. Thus, d2i=−t2it1iT. d1 and d2 point along the principal direction of the section. Based on the geometry relation, the variation of the curvature can be given as(31)δκi=BiδXieΔl,in which(32)Xie=xi−1yi−1xiyixi+1yi+1,Bi=−d2iTΔti+d2i−1TΔti−1,Δti=02×2−I2×2I2×2,Δti−1=−I2×2I2×202×2.I is the identity matrix. This relation corresponds to the continuous form of [33](33)δκ=t×δt′⋅I−t⊗t⋅d1=−δt⋅d2′,where t 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 δκiEIκiΔl of element i suggests the internal force vector of the element is(34)finte=BiTEIκi,in 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:(35)kie=EIBiTBiΔl.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 G have the same expressions with (29).Figure 7
Sketches and numerical model of the inextensible and incompressible beam. (a) A cantilever beam under tip moment. (b) A vertical cantilever beam under vertical load which causes large postbulking displacement. (c) The numerical model of a beam in the planexy.
(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 P=4K2k/π2≈4.0301. And the analytical solution is [34](36)x=2kcosϕ−1Kk,y=2Eϕ,k−Fϕ,kKk,in which k=sin4π/9 and ϕ∈0,2π. Kk 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(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 R˜/fext 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.Figure 8
The numerical results of cables using KDR with constraints. (a) and (b) are the simulation shapes of cantilever beams under tip moment and vertical force, respectively. (c) and (d) are the validations of convergence for beams under tip moment and vertical forces, respectively.
(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 Figure9(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.Figure 9
Circular truss dome. (a) The full dome. (b) Details of one repeated piece.
(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 of(37)ui−ui−1cos15∘−vi−1sin15∘=0,vi+ui−1sin15∘−vi−1cos15∘=0,wi−wi−1=0,where u,v,andw 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 u1=v1=w1=0. 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 R˜/fext vs. the time steps of linear analysis have been plotted in Figure10(b), which verified the convergence of our model.Figure 10
The numerical results of the circular truss dome. (a) The load-deflection curve at tip of the truss tower. Data shown by stars and crosses are taken from Ref. [23] while dashed line and continuous line with circle are our results. (b) Convergence of the linear case.
(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 GtM−1GtT 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.


## 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.



---

*Source: 9895464-2020-03-09.xml*