Adaptive Step-Size Control for Dynamic Relaxation Using Continuous Kinetic Damping

Dynamic relaxation (DR) is a widely used numerical method to determine the static equilibrium of a dynamic system. However, it is difficult to apply conventional DR methods to nonlinear models because they require estimation of a stiffness matrix of the model. To resolve the forementioned problem, a new dynamic relaxation method using continuous kinetic damping (CKDR) was proposed in previous research. The CKDR method does not require any model parameters including the stiffness matrix, and it possesses absolute stability and a second-order convergence rate. However, the convergence rate is proportional to square of the step size, and this may result in a low convergence rate if the selected step size is excessively small.This problem leads to difficulties in the practical use of CKDR.Thus, an adaptive step-sizemethod is proposed in this paper to control the convergence rate of CKDR. The proposed method estimates natural frequency of the model and determines adaptive step size. Static equilibrium simulations were performed for three different models to verify the method. The results revealed that the computational cost of CKDR with a variable step size was very efficient when compared to fixed step sizes and that the convergence rate was also controlled as intended. In addition, the lowest natural frequencies of models in static equilibrium were accurately estimated.


Introduction
Dynamic relaxation (DR) is the most widely used approach for static equilibrium analyses.The basic DR method was proposed in 1965 by A.S. Day [1], and, since then, there have been many applications in the field of structural dynamics [2], geomechanics [3], and biomechanics [4].Research studies continue to develop faster and more efficient DR methods, and the state-of-the-art methods and examples are reported by Rezaiee-Pajand et al. [5][6][7].The DR method basically compels the dynamic system to converge to a static equilibrium through the application of fictitious damping.The type of DR method depends on the type of fictitious damping, and many variants are introduced in [6].Among them, dynamic relaxation with kinetic damping (KDR) is a widely used method proposed in a study by Cundall [8].The kinetic damping in KDR is applied by resetting the velocity to zero when each kinetic energy peak occurs.However, this discontinuous procedure makes it difficult to predict the convergence rate and stability of KDR.Furthermore, these conventional DR methods require the estimation of a stiffness matrix of the model, and, therefore, it is difficult to apply these methods to nonlinear models.
To resolve the aforementioned problem, a new dynamic relaxation method using continuous kinetic damping (CKDR) is proposed in the previous research [9].The basic concept of CKDR involves resetting the velocity to zero at every time step irrespective of kinetic energy.Unlike the conventional DR method, CKDR does not require any model parameters including the stiffness matrix.In addition, it is possible to analyze the convergence rate by eigenvalue analysis.The eigenvalue analysis showed that the convergence rate of the CKDR depends on the step size.An excessively small selected step size slows down the convergence rate and results in an inefficient calculation.Conversely, the selection of an excessively large step size can result in the failure of convergence, since CKDR is an implicit method.Thus, proper step-size selection is important for practical use of the CKDR method due to the forementioned reasons.
This paper proposes an adaptive step-size method of CKDR for proper step-size selection.The proposed method allows maintenance of the convergence rate of CKDR at a desired level while preventing calculation failures.In addition, the lowest natural frequency in static equilibrium can be accurately estimated.The adaptive step-size method is derived from the characteristics of CKDR, and thus it is summarized in Section 2. Additionally, Section 3 shows the derivation of the adaptive step-size method in detail, and an algorithm is provided for its implementation.Section 4 describes the verification of the effectiveness of the adaptive step-size method for CKDR through static equilibrium simulations of three extremely different models.Section 5 provides the conclusions.

Summary of CKDR
As the adaptive step-size method is derived from the characteristics of CKDR, the algorithm of CKDR, which is proposed in [9], is briefly summarized in this section.

. . Algorithm of CKDR.
A second-order dynamic system is represented by a set of nonlinear differential equations as follows: where , d, k, and a denote the time, displacement, velocity, and acceleration vector, respectively.The aim of the CKDR method involves determining a static displacement vector d st that satisfies the static equilibrium condition as follows: According to [9], the CKDR algorithm is represented by the following implicit equations: where ℎ and  denote a time step and step index, respectively.Substituting (3), (4), and ( 5) into (1) results in the error vector function for implicit calculation as follows: A nonlinear equations solver (e.g., the Newton-Raphson method) is used to determine the acceleration vector a +1 that satisfies (6).Subsequently, d +1 is calculated from (3), and the calculation of the k+1 step is completed.Furthermore, a +1 and d +1 converge gradually when this process is repeated.When a +1 becomes smaller than a criterion, then d +1 of the last step is output as a static displacement d st that satisfies (2), and the CKDR algorithm is terminated.
. .Characteristics of CKDR.To characterize CKDR, the equation of motion for a damped system with external force is given as follows: where  n , , and  denote natural frequency, damping ratio, and external force, respectively.The scalar forms of (3) and (4) of the CKDR method are as follows: According to [9], the displacement and acceleration convergence by the CKDR method can be expressed as the following difference equations: where  denotes the amplification factor of the CKDR method: The amplification factor  is a measure of the convergence rate of the CKDR method.This is because smaller values of  lead to faster convergence of  +1 and  +1 from (10) and (11).Larger values of ℎ  are associated with smaller values of  from 1 to 0, and this feature is also confirmed in Figure 1.The x-axis of Figure 1 is replaced with the ratio of time step to the natural period : According to [9], the CKDR method has the following key features.
(3) No oscillation occurs in the process.
(4) It is an implicit method.
The second feature that explains the relationship between step size and convergence rate is the most important feature of CKDR.The second feature suggests that it is necessary to select a larger step size to achieve a higher convergence rate.
Unfortunately, an excessively large step size causes failure of nonlinear equation calculations because this is an implicit method [10].Conversely, an excessively small step size lowers the convergence rate.Thus, proper step-size selection is an essential requirement for CKDR due to the forementioned reasons.To maintain the convergence rate and prevent the failure of the algorithm, the adaptive step size method for CKDR is derived in the following section.

Adaptive Step-Size Control for CKDR
According to (12), the convergence rate of the CKDR method increases as the step size and the natural frequency of the system increase.That is, if the natural frequency is given, the step size that yields the desired convergence rate can be determined.However, the natural frequency of the dynamic system is not generally given, and changes occur depending on the configuration in the case of the nonlinear model.Furthermore, in the case of a multi-degree-of-freedom (multi-DOF) model, its natural frequencies correspond to multiple numbers, and thus it cannot be expressed as a single scalar value.To resolve these problems, an estimation method for the natural frequency is proposed in the following section.
. .Estimation of Natural Frequency . . .Definition of the Estimated Natural Frequency.To estimate the natural frequency, it is assumed that the natural frequency  n and the damping ratio  and the external force  of k and k+1 steps are approximately the same.This assumption allows the expression of (7) for k and k+1 steps as follows: The following equation is derived by subtracting ( 14) and (15) and applying (9) as follows: From ( 16), the positive natural frequency of the system can be estimated as follows: However, this equation is only for single-DOF models and cannot be applied to multi-DOF models.Therefore, the estimated natural frequency for multi-DOF models is defined by modifying (17) as follows: However, without the physical meaning of this definition, this is just one of many definitions.To analyze the characteristics of the estimated natural frequency defined by (18), matrix analysis was performed, and Sections 3.1.2and 3.1.3describe the convergence characteristics and the upper bound of the estimated natural frequency.

. . . Convergence Characteristics of the Estimated Natural
Frequency.In order to analyze the physical meaning of the definition of ωn , the vector form of ( 7) is given as follows: where M, C, and K are the mass, damping, and stiffness matrices.By repeating the derivation procedure for single-DOF case, the matrix form of ( 16) is derived as follows: This equation can be expressed as the product of a matrix and a vector as follows: where The submultiplicative property [11] states that Ab +1 , A +1 , and b +1 in (21) must satisfy the following conditional equality: Substituting ( 22) and (24) into the LHS of (25), it can be expressed with the estimated natural frequency ωn of (18) as follows: where A +1 is assumed to be a positive semidefinite matrix as both M +1 and K +1 can be a positive definite matrix or a positive semidefinite matrix [12].Therefore, A +1 and b +1 satisfy the following equations [11]: where  is a modal index,  n(+1,) (<  n(+1,+1) ) is the  th natural frequency, x +1, is the  th eigenvector, and  +1, is the  th linear combination coefficient.To express the estimated Mathematical Problems in Engineering natural frequency ωn(+1) with  n(+1,) and x +1, , substituting ( 27) and ( 28) into (26) yields the following: According to the convergence characteristics of CKDR obtained from ( 10) and ( 12), the higher the natural frequency  n(+1,) is, the faster the corresponding linear combination coefficient  +1, decreases.As the CKDR process is repeated, the coefficient  +1,1 corresponding to the lowest natural frequency becomes much more dominant over all the other coefficients, that is,  ∞,1 ≫  ∞,{2,⋅⋅⋅ ,} , which is reflected in (29) as follows: ) ) From (30), as the CKDR process proceeds, the estimated natural frequency ωn(+1) defined by (18) converges to the lowest natural frequency  min n of a dynamic model in static equilibrium.

. . . e Upper Bound of the Estimated Natural Frequency.
From the definition of a spectral norm of a general matrix [13], ‖A +1 ‖ is equal to the largest singular value  max +1 of A +1 as follows: To derive the upper bound of the estimated natural frequency ωn(+1) , substituting (26) and (31) into (25) yields the following: By applying the submultiplicative property to (27), the largest singular value  max +1 has the following conditional equality: Because  n(+1,{1,⋅⋅⋅ ,−1}) <  n(+1,) def     max n(+1) , the following conditional equality is satisfied from (33): According to (32) and (34), the estimated natural frequency ωn(+1) is less than or equal to the square root of the largest singular value  max +1 , which is greater than or equal to the highest natural frequency  max n(+1) .
Note that the equality of (34) is satisfied only when A +1 is a symmetric matrix according to [14].If M +1 and K +1 are diagonal matrices, A +1 is a symmetric matrix.That is, for uncoupled dynamic systems, the equality is satisfied as follows: This is a special case; otherwise the equivalence is not guaranteed.In Sections 3.1.2and 3.1.3,the convergence characteristics and the upper bound of the estimated natural frequency were analyzed, and these are confirmed in the verification simulations.In the following section, the adaptive step-size method is derived by using the estimated natural frequency, and the detailed derivation is described.

. . Adaptive
Step-Size Calculation.According to (12), the amplification factor  is a function of the product of the natural frequency  n and the step size ℎ.Applying this, the step-size adaptation can be defined as the problem of determining a next step size ℎ +1 that makes the amplification factor correspond to the desired value as follows: where  d is the desired amplification factor and the desired discrete natural Ω d is defined as follows: By substituting (37) into (36), the former can be expressed as From (37), the next step size ℎ +1 , which maintains the convergence rate at a desired level, can be obtained as follows: According to (36), a very high convergence rate can be obtained if an extremely large Ω d is selected.Unfortunately, In this case, the next step size ℎ +1 of (39) is adapted with respect to an excessively large value, and this leads to failure in solving the nonlinear equations.To prevent failure while obtaining a sufficient convergence rate, we recommend 2 as the value of the desired discrete natural frequency Ω d : Note that this selection makes the next step size ℎ +1 equivalent to the estimated natural period T+1 of (13) as follows: to 0.0247.This means the next step size ℎ +1 is adapted to maintain  ≈  d for each step as follows: Additionally, the natural frequency of nonlinear models can change very quickly when the analysis proceeds.This causes a sudden change in step size, and this can result in failure to determine the solutions to nonlinear equations [10].In order to prevent this, the next step size ℎ +1 should be limited as follows [10]: where  min and  max denote a minimum and maximum stepsize ratio, respectively.
. .Procedure of Step-Size Adaptation.Based on the proposed framework, the pseudocode of the CKDR with the adaptive step-size method is provided in Algorithm 1.The parameters selected for the algorithm correspond to the desired discrete natural frequency of Ω d (0 < Ω d ≤ 2), the minimum step size ratio  max (0 <  max < 1), the maximum step size ratio  max (> 1), and the acceleration tolerance  p (> 0).It should be noted that any model parameters, such as mass and stiffness matrices, are not used.

Verification Simulations
The verification of the adaptive step-size method for CKDR involved the simulation of three models, namely, 3-DOF modal model, quarter-car model, and multiple spring-mass model.These simple but representative models were chosen to clearly demonstrate the features of the proposed step-size control method and to enable the reproducibility of the verification results.The built-in function of Matlab tm , "fsolve.m," is used to solve nonlinear equations.The parameters for proposed method were set and corresponded to Ω d = 2,  min = 0.5,  max = 2.0, and  p = 10 −4 .These conditions were the same for all simulation cases.Step k . .-DOF Modal Model.To verify the frequency estimation and convergence control of the proposed method, a static equilibrium analysis of the following undamped 3-DOF modal model is performed as follows: where  n,{1,2,3} = {1, 10, 100}rad/s,   ( 0 ) = 1m, and the initial step size ℎ 0 = 0.01s.The displacement convergence of the model is shown in Figure 2. The progress of the simulation first results in the convergence of the high-frequency response, and this is followed by the gradual convergence of the low-frequency response.This is due to the characteristic of CKDR in which an increase in the frequency increases the convergence rate.At the 12 th step, the end condition was satisfied, and the simulation was terminated.
In Figure 3, the natural frequency estimated from (18) converged to the lowest natural frequency (1 rad/s), which can be deduced from (30).The estimated natural frequencies of all the steps were less than 100 rad/s, which corresponds to the square root of the largest singular value of this model, according to (32).As this ideal model is uncoupled, the square root of the largest singular value is equal to the highest natural frequency according to (35).
In Figure 4, the step size increased exponentially from the initial value of 0.01 to converge to a specific value of 6.283.The step size was calculated by using (39), and the exponential increase is due to the restriction of (43) to prevent sudden changes in step size.The converged step size of 6.28 s corresponds to the longest natural period of the 3-DOF model from (41).Note that this is because the desired discrete natural frequency was chosen as 2.
From the estimated natural frequency and step size, the amplification factor was calculated by (12), and the result is shown in Figure 5.The amplification factor converged to 0.0247, and this corresponds to the desired value of (42).This shows that the convergence rate was controlled as well as intended.The fluctuation of the value is due to the limitation in the step-size change, which was recovered after the step size was sufficiently increased.
. .Quarter-Car Model.A ground vehicle is a typical model that requires static equilibrium analysis prior to dynamic analysis.A 2-DOF quarter car was modeled as shown in Figure 6, and the parameters of the quarter-car model are given in Table 1.
Figure 7 shows the results of the CKDR with the dynamic analysis.In the dynamic analysis, the quarter-car model began to fall to the ground owing to gravity and contacted the road surface at approximately 0.4 s.The model then bounced from the road surface because of the reactive force of the suspension and tire.In contrast, the results of the CKDR algorithm converged without bouncing.
To verify the effectiveness of the adaptive step-size method, the variable step is compared with the fixed steps of 0.05, 0.15, and 0.25 as shown in Figure 8.In the case of  the fixed step, a smaller step size requires a more exponential increase in the number of iterations.It should be noted that, despite the same initial step size of 0.05, the result of the variable step converged significantly faster than that of the fixed step.These results are summarized in Table 2, and the displacement convergence of variable and fixed steps is shown in Figure 7.
The estimated natural frequencies for variable step size calculation are shown in Figure 9. Initially, the natural frequency was estimated as zero because the acceleration did not change during the free fall.Zero frequency corresponds to the natural frequency of the rigid body mode.Following contact Step k  with the road surface, the value increased and converged to 7.535 rad/s, which is the same as the lowest natural frequency of the quarter car.In Figure 10, the amplification factor also converged to 0.0247, and this was predicted by (42).
. .Multi-DOF Mass-Spring Model.An important use of static equilibrium analysis corresponds to determining the static configuration of a structure, that is, the form-finding problem.For example, a simple cable model with both ends fixed is given as shown in Figure 11.The simple cable model consists of 100 nodes of 0.01 kg, and each node is connected with a spring of 250 N/m at intervals of 0.1 m.Additionally, a flat barrier exists at -2 m, and this poses an obstacle to the cable deflection.Figure 11 shows the shape change of the simple cable model during static equilibrium analysis.As expected, the deflection of the cable is initiated by gravity, and it equilibrated on the barrier at -2 m. Figure 12 shows the adaptation of variable step size when compared with the fixed step sizes, and the iteration numbers are summarized in Table 3.Despite the relatively high DOF, the iteration number of the cable model was similar to that of the 2-DOF quarter-car model.In Figure 13, the amplification factor of the cable model converged to a predicted value of 0.0247 as shown in the previous models.

Conclusion
In the previous research [9], a new dynamic relaxation method termed as CKDR was proposed for static equilibrium analysis.In contrast to conventional DR methods, the CKDR method does not require a stiffness matrix of a model, and it exhibits a second-order convergence rate.However, many iterations are required if an excessively small step size is selected, and this limits the practical use of CKDR.
Hence, an adaptive step-size method was proposed in this paper to control the convergence rate of the CKDR Step size (s) 0.  method.The proposed adaptive step-size method estimated the natural frequency and determined the step size such that the convergence rate corresponded to the desired level.The natural frequency was estimated using only the displacement and acceleration responses, and the approach was proven by the matrix algebra.
The static equilibrium simulations of the simple 3-DOF model, the quarter-car model, and the simple cable model were performed to verify the adaptive step-size method.The results indicated that static equilibrium simulation was terminated with less than 10 iterations in all the model cases, although the degrees of freedom and boundary conditions of the models were very different.Furthermore, the lowest natural frequencies of the models in static equilibrium were accurately estimated, and the amplification factors converged to the desired values.In conclusion, the proposed adaptive step-size method does not require any model parameters, and it makes the CKDR method efficient and suitable for nonlinear models.

Figure 3 :
Figure 3: Estimated natural frequency of the simple 3-DOF model.

Figure 8 :
Figure 8: Step-size adaptation of the quarter-car model.

Figure 9 :
Figure 9: Estimated natural frequency of the quarter-car model.

Figure 10 :Figure 11 :
Figure 10: Amplification factor of the quarter-car model.

Figure 12 :
Figure 12: Step-size adaptation of the simple cable model.

Figure 13 :
Figure 13: Amplification factor of the simple cable model.

Table 1 :
Parameters of the quarter-car model.

Table 2 :
Iteration numbers of the quarter-car model.

Table 3 :
Iteration numbers of the simple cable model.