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

Authors: Samuel Jung; Tae-Yun Kim; Wan-Suk Yoo
Journal: Mathematical Problems in Engineering
 (2018)

---

## Abstract

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-size method 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.


---

## Body

## 1. 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–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 Section2. 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.


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


### 2.1. Algorithm of CKDR

A second-order dynamic system is represented by a set of nonlinear differential equations as follows:(1)at=ft,dt,vt,where t, d, v, and a denote the time, displacement, velocity, and acceleration vector, respectively. The aim of the CKDR method involves determining a static displacement vector dst that satisfies the static equilibrium condition as follows:(2)0=ft0,dst,0.According to [9], the CKDR algorithm is represented by the following implicit equations:(3)dk+1=dk+h2ak+1,(4)vk+1=0,(5)tk+1=t0,where h and k 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:(6)ψak+1=ft0,dk+1,0-ak+1=0.A nonlinear equations solver (e.g., the Newton–Raphson method) is used to determine the acceleration vector ak+1 that satisfies (6). Subsequently, dk+1 is calculated from (3), and the calculation of thek+1 step is completed. Furthermore, ak+1 and dk+1 converge gradually when this process is repeated. When ak+1 becomes smaller than a criterion, then dk+1 of the last step is output as a static displacement dst that satisfies (2), and the CKDR algorithm is terminated.


### 2.2. Characteristics of CKDR

To characterize CKDR, the equation of motion for a damped system with external force is given as follows:(7)a+2ζωnv+ωn2d=p,where ωn, ζ, and p denote natural frequency, damping ratio, and external force, respectively. The scalar forms of (3) and (4) of the CKDR method are as follows:(8)dk+1=dk+h2ak+1,(9)vk+1=0.According to [9], the displacement and acceleration convergence by the CKDR method can be expressed as the following difference equations:(10)dk+1-dst=ηdk-dst,(11)ak+1=ηak,where η denotes the amplification factor of the CKDR method: (12)η=11+hωn2.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 dk+1 and ak+1 from (10) and (11). Larger values of hωn 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 T:(13)hT=hωn2π.According to [9], the CKDR method has the following key features.(1)
It is absolutely stable (A-stable).(2)
The convergence rate is proportional toh2 (2-order method).(3)
No oscillation occurs in the process.(4)
It is an implicit method.Figure 1
Amplification factor of the CKDR.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.


## 2.1. Algorithm of CKDR

A second-order dynamic system is represented by a set of nonlinear differential equations as follows:(1)at=ft,dt,vt,where t, d, v, and a denote the time, displacement, velocity, and acceleration vector, respectively. The aim of the CKDR method involves determining a static displacement vector dst that satisfies the static equilibrium condition as follows:(2)0=ft0,dst,0.According to [9], the CKDR algorithm is represented by the following implicit equations:(3)dk+1=dk+h2ak+1,(4)vk+1=0,(5)tk+1=t0,where h and k 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:(6)ψak+1=ft0,dk+1,0-ak+1=0.A nonlinear equations solver (e.g., the Newton–Raphson method) is used to determine the acceleration vector ak+1 that satisfies (6). Subsequently, dk+1 is calculated from (3), and the calculation of thek+1 step is completed. Furthermore, ak+1 and dk+1 converge gradually when this process is repeated. When ak+1 becomes smaller than a criterion, then dk+1 of the last step is output as a static displacement dst that satisfies (2), and the CKDR algorithm is terminated.


## 2.2. Characteristics of CKDR

To characterize CKDR, the equation of motion for a damped system with external force is given as follows:(7)a+2ζωnv+ωn2d=p,where ωn, ζ, and p denote natural frequency, damping ratio, and external force, respectively. The scalar forms of (3) and (4) of the CKDR method are as follows:(8)dk+1=dk+h2ak+1,(9)vk+1=0.According to [9], the displacement and acceleration convergence by the CKDR method can be expressed as the following difference equations:(10)dk+1-dst=ηdk-dst,(11)ak+1=ηak,where η denotes the amplification factor of the CKDR method: (12)η=11+hωn2.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 dk+1 and ak+1 from (10) and (11). Larger values of hωn 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 T:(13)hT=hωn2π.According to [9], the CKDR method has the following key features.(1)
It is absolutely stable (A-stable).(2)
The convergence rate is proportional toh2 (2-order method).(3)
No oscillation occurs in the process.(4)
It is an implicit method.Figure 1
Amplification factor of the CKDR.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.


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


### 3.1. Estimation of Natural Frequency

#### 3.1.1. 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 p ofk andk+1 steps are approximately the same. This assumption allows the expression of (7) fork andk+1 steps as follows:(14)ak+2ζk+1ωnk+1vk+ωnk+12dk=pk+1,(15)ak+1+2ζk+1ωnk+1vk+1+ωnk+12dk+1=pk+1.The following equation is derived by subtracting (14) and (15) and applying (9) as follows:(16)ak+1-ak+ωnk+12dk+1-dk=0.From (16), the positive natural frequency of the system can be estimated as follows:(17)ω^nk+1=ak+1-akdk+1-dk.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:(18)ω^nk+1=ak+1-akdk+1-dk.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.2 and 3.1.3 describe the convergence characteristics and the upper bound of the estimated natural frequency.


#### 3.1.2. 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:(19)Ma+Cv+Kd=p,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:(20)ak+1-ak=-Mk+1-1Kk+1dk+1-dk.This equation can be expressed as the product of a matrix and a vector as follows:(21)Abk+1=Ak+1·bk+1,where (22)Abk+1=ak+1-ak,(23)Ak+1=-Mk+1-1Kk+1,(24)bk+1=dk+1-dk.The submultiplicative property [11] states that Abk+1, Ak+1, and bk+1 in (21) must satisfy the following conditional equality:(25)Abk+1bk+1≤Ak+1.Substituting (22) and (24) into the LHS of (25), it can be expressed with the estimated natural frequency ω^n of (18) as follows: (26)ω^nk+12=Ak+1bk+1bk+1=ak+1-akdk+1-dk,where Ak+1 is assumed to be a positive semidefinite matrix as both Mk+1 and Kk+1 can be a positive definite matrix or a positive semidefinite matrix [12]. Therefore, Ak+1 and bk+1 satisfy the following equations [11]: (27)Ak+1xk+1,i=ωnk+1,i2xk+1,i,(28)bk+1=∑i=1nck+1,ixk+1,i,where i is a modal index, ωn(k+1,i)<ωn(k+1,i+1) is the ith natural frequency, xk+1,i is the ith eigenvector, and ck+1,i is the ith linear combination coefficient. To express the estimated natural frequency ω^n(k+1) with ωn(k+1,i) and xk+1,i, substituting (27) and (28) into (26) yields the following: (29)ω^nk+12=Ak+1bk+1bk+1=Ak+1bk+1TAk+1bk+1bk+1Tbk+1=∑i=1nck+1,iωnk+1,i2xk+1,iT∑i=1nck+1,iωnk+1,i2xk+1,i∑i=1nck+1,ixk+1,iT∑i=1nck+1,ixk+1,i.According to the convergence characteristics of CKDR obtained from (10) and (12), the higher the natural frequency ωn(k+1,i) is, the faster the corresponding linear combination coefficient ck+1,i decreases. As the CKDR process is repeated, the coefficient ck+1,1 corresponding to the lowest natural frequency becomes much more dominant over all the other coefficients, that is, c∞,1≫c∞,{2,⋯,n}, which is reflected in (29) as follows:(30)limk→∞⁡ω^nk+1=c∞,12ωn∞,14x∞,1Tx∞,1c∞,12x∞,1Tx∞,11/4=ωn∞,1=defωnmin.From (30), as the CKDR process proceeds, the estimated natural frequency ω^nk+1 defined by (18) converges to the lowest natural frequency ωnmin of a dynamic model in static equilibrium.


#### 3.1.3. The Upper Bound of the Estimated Natural Frequency

From the definition of a spectral norm of a general matrix [13], Ak+1 is equal to the largest singular value σk+1max of Ak+1 as follows:(31)Ak+1=maxi⁡σk+1,i=defσk+1max.To derive the upper bound of the estimated natural frequency ω^nk+1, substituting (26) and (31) into (25) yields the following: (32)ω^nk+1≤σk+1max.By applying the submultiplicative property to (27), the largest singular value σk+1max has the following conditional equality:(33)ωnk+1,i2xk+1,i=Ak+1xk+1,i≤Ak+1xk+1,i=σk+1maxxk+1,i.Because ωn(k+1,{1,⋯,n-1})<ωn(k+1,n)=defωn(k+1)max, the following conditional equality is satisfied from (33):(34)σk+1max≥ωnk+1max.According to (32) and (34), the estimated natural frequency ω^nk+1 is less than or equal to the square root of the largest singular value σk+1max, which is greater than or equal to the highest natural frequency ωn(k+1)max.Note that the equality of (34) is satisfied only when Ak+1 is a symmetric matrix according to [14]. If Mk+1 and Kk+1 are diagonal matrices, Ak+1 is a symmetric matrix. That is, for uncoupled dynamic systems, the equality is satisfied as follows: (35)σk+1max=ωnk+1max.This is a special case; otherwise the equivalence is not guaranteed.In Sections3.1.2 and 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.


### 3.2. 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 h. Applying this, the step-size adaptation can be defined as the problem of determining a next step size hk+1 that makes the amplification factor correspond to the desired value as follows: (36)ηd=11+hk+12ω^nk+12,where ηd is the desired amplification factor and the desired discrete natural Ωd is defined as follows:(37)Ωd=defhk+1ω^nk+1.By substituting (37) into (36), the former can be expressed as(38)ηd=11+Ωd2.From (37), the next step size hk+1, which maintains the convergence rate at a desired level, can be obtained as follows:(39)hk+1=Ωdω^nk+1.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 hk+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: (40)Ωd=2π.Note that this selection makes the next step size hn+1 equivalent to the estimated natural period T^k+1 of (13) as follows:(41)hn+1=T^k+1=2πω^nk+1.Moreover, substituting (40) into (36) indicates that the desired amplification factor ηd approximately corresponds to 0.0247. This means the next step size hk+1 is adapted to maintain η≈ηd for each step as follows:(42)ηd=11+2π2≈0.0247.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 hk+1 should be limited as follows [10]:(43)rminhk<hk+1<rmaxhk,where rmin and rmax denote a minimum and maximum step-size ratio, respectively.


### 3.3. Procedure of Step-Size Adaptation

Based on the proposed framework, the pseudocode of the CKDR with the adaptive step-size method is provided in Algorithm1. The parameters selected for the algorithm correspond to the desired discrete natural frequency of Ωd(0<Ωd≤2π), the minimum step size ratio rmax(0<rmax<1), the maximum step size ratio rmax(>1), and the acceleration tolerance ap(>0). It should be noted that any model parameters, such as mass and stiffness matrices, are not used.Algorithm 1:Pseudocode of CKDR with the adaptive step-size method.
Input:d0, t0, h0, ap, Ωd, rmin, rmax a 0 ← 0, k←0 While Do Findak+1 satisfying ψ(ak+1)=0 dk+1←dk+hk2ak+1 ω^n(k+1)←ak+1-akdk+1-dk hk+1←Ωdω^n(k+1) hk+1=min(max(hk+1,rminhk),rmaxhk) Ifak+1∞<ap, Then break while, End k←k+1 End Return:dk+1, ω^n(k+1)


## 3.1. Estimation of Natural Frequency

### 3.1.1. 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 p ofk andk+1 steps are approximately the same. This assumption allows the expression of (7) fork andk+1 steps as follows:(14)ak+2ζk+1ωnk+1vk+ωnk+12dk=pk+1,(15)ak+1+2ζk+1ωnk+1vk+1+ωnk+12dk+1=pk+1.The following equation is derived by subtracting (14) and (15) and applying (9) as follows:(16)ak+1-ak+ωnk+12dk+1-dk=0.From (16), the positive natural frequency of the system can be estimated as follows:(17)ω^nk+1=ak+1-akdk+1-dk.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:(18)ω^nk+1=ak+1-akdk+1-dk.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.2 and 3.1.3 describe the convergence characteristics and the upper bound of the estimated natural frequency.


### 3.1.2. 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:(19)Ma+Cv+Kd=p,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:(20)ak+1-ak=-Mk+1-1Kk+1dk+1-dk.This equation can be expressed as the product of a matrix and a vector as follows:(21)Abk+1=Ak+1·bk+1,where (22)Abk+1=ak+1-ak,(23)Ak+1=-Mk+1-1Kk+1,(24)bk+1=dk+1-dk.The submultiplicative property [11] states that Abk+1, Ak+1, and bk+1 in (21) must satisfy the following conditional equality:(25)Abk+1bk+1≤Ak+1.Substituting (22) and (24) into the LHS of (25), it can be expressed with the estimated natural frequency ω^n of (18) as follows: (26)ω^nk+12=Ak+1bk+1bk+1=ak+1-akdk+1-dk,where Ak+1 is assumed to be a positive semidefinite matrix as both Mk+1 and Kk+1 can be a positive definite matrix or a positive semidefinite matrix [12]. Therefore, Ak+1 and bk+1 satisfy the following equations [11]: (27)Ak+1xk+1,i=ωnk+1,i2xk+1,i,(28)bk+1=∑i=1nck+1,ixk+1,i,where i is a modal index, ωn(k+1,i)<ωn(k+1,i+1) is the ith natural frequency, xk+1,i is the ith eigenvector, and ck+1,i is the ith linear combination coefficient. To express the estimated natural frequency ω^n(k+1) with ωn(k+1,i) and xk+1,i, substituting (27) and (28) into (26) yields the following: (29)ω^nk+12=Ak+1bk+1bk+1=Ak+1bk+1TAk+1bk+1bk+1Tbk+1=∑i=1nck+1,iωnk+1,i2xk+1,iT∑i=1nck+1,iωnk+1,i2xk+1,i∑i=1nck+1,ixk+1,iT∑i=1nck+1,ixk+1,i.According to the convergence characteristics of CKDR obtained from (10) and (12), the higher the natural frequency ωn(k+1,i) is, the faster the corresponding linear combination coefficient ck+1,i decreases. As the CKDR process is repeated, the coefficient ck+1,1 corresponding to the lowest natural frequency becomes much more dominant over all the other coefficients, that is, c∞,1≫c∞,{2,⋯,n}, which is reflected in (29) as follows:(30)limk→∞⁡ω^nk+1=c∞,12ωn∞,14x∞,1Tx∞,1c∞,12x∞,1Tx∞,11/4=ωn∞,1=defωnmin.From (30), as the CKDR process proceeds, the estimated natural frequency ω^nk+1 defined by (18) converges to the lowest natural frequency ωnmin of a dynamic model in static equilibrium.


### 3.1.3. The Upper Bound of the Estimated Natural Frequency

From the definition of a spectral norm of a general matrix [13], Ak+1 is equal to the largest singular value σk+1max of Ak+1 as follows:(31)Ak+1=maxi⁡σk+1,i=defσk+1max.To derive the upper bound of the estimated natural frequency ω^nk+1, substituting (26) and (31) into (25) yields the following: (32)ω^nk+1≤σk+1max.By applying the submultiplicative property to (27), the largest singular value σk+1max has the following conditional equality:(33)ωnk+1,i2xk+1,i=Ak+1xk+1,i≤Ak+1xk+1,i=σk+1maxxk+1,i.Because ωn(k+1,{1,⋯,n-1})<ωn(k+1,n)=defωn(k+1)max, the following conditional equality is satisfied from (33):(34)σk+1max≥ωnk+1max.According to (32) and (34), the estimated natural frequency ω^nk+1 is less than or equal to the square root of the largest singular value σk+1max, which is greater than or equal to the highest natural frequency ωn(k+1)max.Note that the equality of (34) is satisfied only when Ak+1 is a symmetric matrix according to [14]. If Mk+1 and Kk+1 are diagonal matrices, Ak+1 is a symmetric matrix. That is, for uncoupled dynamic systems, the equality is satisfied as follows: (35)σk+1max=ωnk+1max.This is a special case; otherwise the equivalence is not guaranteed.In Sections3.1.2 and 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.


## 3.1.1. 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 p ofk andk+1 steps are approximately the same. This assumption allows the expression of (7) fork andk+1 steps as follows:(14)ak+2ζk+1ωnk+1vk+ωnk+12dk=pk+1,(15)ak+1+2ζk+1ωnk+1vk+1+ωnk+12dk+1=pk+1.The following equation is derived by subtracting (14) and (15) and applying (9) as follows:(16)ak+1-ak+ωnk+12dk+1-dk=0.From (16), the positive natural frequency of the system can be estimated as follows:(17)ω^nk+1=ak+1-akdk+1-dk.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:(18)ω^nk+1=ak+1-akdk+1-dk.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.2 and 3.1.3 describe the convergence characteristics and the upper bound of the estimated natural frequency.


## 3.1.2. 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:(19)Ma+Cv+Kd=p,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:(20)ak+1-ak=-Mk+1-1Kk+1dk+1-dk.This equation can be expressed as the product of a matrix and a vector as follows:(21)Abk+1=Ak+1·bk+1,where (22)Abk+1=ak+1-ak,(23)Ak+1=-Mk+1-1Kk+1,(24)bk+1=dk+1-dk.The submultiplicative property [11] states that Abk+1, Ak+1, and bk+1 in (21) must satisfy the following conditional equality:(25)Abk+1bk+1≤Ak+1.Substituting (22) and (24) into the LHS of (25), it can be expressed with the estimated natural frequency ω^n of (18) as follows: (26)ω^nk+12=Ak+1bk+1bk+1=ak+1-akdk+1-dk,where Ak+1 is assumed to be a positive semidefinite matrix as both Mk+1 and Kk+1 can be a positive definite matrix or a positive semidefinite matrix [12]. Therefore, Ak+1 and bk+1 satisfy the following equations [11]: (27)Ak+1xk+1,i=ωnk+1,i2xk+1,i,(28)bk+1=∑i=1nck+1,ixk+1,i,where i is a modal index, ωn(k+1,i)<ωn(k+1,i+1) is the ith natural frequency, xk+1,i is the ith eigenvector, and ck+1,i is the ith linear combination coefficient. To express the estimated natural frequency ω^n(k+1) with ωn(k+1,i) and xk+1,i, substituting (27) and (28) into (26) yields the following: (29)ω^nk+12=Ak+1bk+1bk+1=Ak+1bk+1TAk+1bk+1bk+1Tbk+1=∑i=1nck+1,iωnk+1,i2xk+1,iT∑i=1nck+1,iωnk+1,i2xk+1,i∑i=1nck+1,ixk+1,iT∑i=1nck+1,ixk+1,i.According to the convergence characteristics of CKDR obtained from (10) and (12), the higher the natural frequency ωn(k+1,i) is, the faster the corresponding linear combination coefficient ck+1,i decreases. As the CKDR process is repeated, the coefficient ck+1,1 corresponding to the lowest natural frequency becomes much more dominant over all the other coefficients, that is, c∞,1≫c∞,{2,⋯,n}, which is reflected in (29) as follows:(30)limk→∞⁡ω^nk+1=c∞,12ωn∞,14x∞,1Tx∞,1c∞,12x∞,1Tx∞,11/4=ωn∞,1=defωnmin.From (30), as the CKDR process proceeds, the estimated natural frequency ω^nk+1 defined by (18) converges to the lowest natural frequency ωnmin of a dynamic model in static equilibrium.


## 3.1.3. The Upper Bound of the Estimated Natural Frequency

From the definition of a spectral norm of a general matrix [13], Ak+1 is equal to the largest singular value σk+1max of Ak+1 as follows:(31)Ak+1=maxi⁡σk+1,i=defσk+1max.To derive the upper bound of the estimated natural frequency ω^nk+1, substituting (26) and (31) into (25) yields the following: (32)ω^nk+1≤σk+1max.By applying the submultiplicative property to (27), the largest singular value σk+1max has the following conditional equality:(33)ωnk+1,i2xk+1,i=Ak+1xk+1,i≤Ak+1xk+1,i=σk+1maxxk+1,i.Because ωn(k+1,{1,⋯,n-1})<ωn(k+1,n)=defωn(k+1)max, the following conditional equality is satisfied from (33):(34)σk+1max≥ωnk+1max.According to (32) and (34), the estimated natural frequency ω^nk+1 is less than or equal to the square root of the largest singular value σk+1max, which is greater than or equal to the highest natural frequency ωn(k+1)max.Note that the equality of (34) is satisfied only when Ak+1 is a symmetric matrix according to [14]. If Mk+1 and Kk+1 are diagonal matrices, Ak+1 is a symmetric matrix. That is, for uncoupled dynamic systems, the equality is satisfied as follows: (35)σk+1max=ωnk+1max.This is a special case; otherwise the equivalence is not guaranteed.In Sections3.1.2 and 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.


## 3.2. 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 h. Applying this, the step-size adaptation can be defined as the problem of determining a next step size hk+1 that makes the amplification factor correspond to the desired value as follows: (36)ηd=11+hk+12ω^nk+12,where ηd is the desired amplification factor and the desired discrete natural Ωd is defined as follows:(37)Ωd=defhk+1ω^nk+1.By substituting (37) into (36), the former can be expressed as(38)ηd=11+Ωd2.From (37), the next step size hk+1, which maintains the convergence rate at a desired level, can be obtained as follows:(39)hk+1=Ωdω^nk+1.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 hk+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: (40)Ωd=2π.Note that this selection makes the next step size hn+1 equivalent to the estimated natural period T^k+1 of (13) as follows:(41)hn+1=T^k+1=2πω^nk+1.Moreover, substituting (40) into (36) indicates that the desired amplification factor ηd approximately corresponds to 0.0247. This means the next step size hk+1 is adapted to maintain η≈ηd for each step as follows:(42)ηd=11+2π2≈0.0247.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 hk+1 should be limited as follows [10]:(43)rminhk<hk+1<rmaxhk,where rmin and rmax denote a minimum and maximum step-size ratio, respectively.


## 3.3. Procedure of Step-Size Adaptation

Based on the proposed framework, the pseudocode of the CKDR with the adaptive step-size method is provided in Algorithm1. The parameters selected for the algorithm correspond to the desired discrete natural frequency of Ωd(0<Ωd≤2π), the minimum step size ratio rmax(0<rmax<1), the maximum step size ratio rmax(>1), and the acceleration tolerance ap(>0). It should be noted that any model parameters, such as mass and stiffness matrices, are not used.Algorithm 1:Pseudocode of CKDR with the adaptive step-size method.
Input:d0, t0, h0, ap, Ωd, rmin, rmax a 0 ← 0, k←0 While Do Findak+1 satisfying ψ(ak+1)=0 dk+1←dk+hk2ak+1 ω^n(k+1)←ak+1-akdk+1-dk hk+1←Ωdω^n(k+1) hk+1=min(max(hk+1,rminhk),rmaxhk) Ifak+1∞<ap, Then break while, End k←k+1 End Return:dk+1, ω^n(k+1)


## 4. 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 ofMatlabtm, “fsolve.m,” is used to solve nonlinear equations. The parameters for proposed method were set and corresponded to Ωd=2π, rmin=0.5, rmax=2.0, and ap=10-4. These conditions were the same for all simulation cases.


### 4.1. 3-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:(44)ai+ωn,i2di=0,i=1,2,3,where ωn,1,2,3=1,10,100rad/s, di(t0)=1m, and the initial step size h0=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 12th step, the end condition was satisfied, and the simulation was terminated.Figure 2
Displacement convergence of the simple 3-DOF model.In Figure3, 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).Figure 3
Estimated natural frequency of the simple 3-DOF model.In Figure4, 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π.Figure 4
Step-size change of the simple 3-DOF model.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.Figure 5
Amplification factor of the simple 3-DOF model.


### 4.2. 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 Figure6, and the parameters of the quarter-car model are given in Table 1.Table 1
Parameters of the quarter-car model.
Parameters Values Sprung mass (ms) 320 kg Unsprung mass (mu) 32 kg Stiffness of suspension (ks) 20000 N/m Stiffness of tire (ku) 200000 N/m Damping of suspension (cs) 1700 Ns/m Free height of sprung mass (Ls) 1.0 m Free height of unsprung mass (Lu) 0.3 mFigure 6
Quarter-car model.Figure7 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.Figure 7
Displacement convergence of the quarter car.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 Figure8. 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.Table 2
Iteration numbers of the quarter-car model.
Step size Iteration number CPU time (s) Ratio(%) Fixed h = 0.25 10 0.27 44 h = 0.15 19 0.32 52 h = 0.05 117 0.61 100  Variable h 0 = 0.05 8 0.26 42Figure 8
Step-size adaptation of the quarter-car model.The estimated natural frequencies for variable step size calculation are shown in Figure9. 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 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).Figure 9
Estimated natural frequency of the quarter-car model.Figure 10
Amplification factor of the quarter-car model.


### 4.3. 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 Figure11. 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 11
Deflection of the cable model.Figure12 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.Table 3
Iteration numbers of the simple cable model.
Step size Iteration number CPU time (s) Ratio(%) Fixed h = 0.5 8 3.7 28 h = 0.3 12 4.2 32 h = 0.1 64 12.8 100  Variable h 0 = 0.1 8 4.5 35Figure 12
Step-size adaptation of the simple cable model.Figure 13
Amplification factor of the simple cable model.


## 4.1. 3-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:(44)ai+ωn,i2di=0,i=1,2,3,where ωn,1,2,3=1,10,100rad/s, di(t0)=1m, and the initial step size h0=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 12th step, the end condition was satisfied, and the simulation was terminated.Figure 2
Displacement convergence of the simple 3-DOF model.In Figure3, 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).Figure 3
Estimated natural frequency of the simple 3-DOF model.In Figure4, 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π.Figure 4
Step-size change of the simple 3-DOF model.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.Figure 5
Amplification factor of the simple 3-DOF model.


## 4.2. 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 Figure6, and the parameters of the quarter-car model are given in Table 1.Table 1
Parameters of the quarter-car model.
Parameters Values Sprung mass (ms) 320 kg Unsprung mass (mu) 32 kg Stiffness of suspension (ks) 20000 N/m Stiffness of tire (ku) 200000 N/m Damping of suspension (cs) 1700 Ns/m Free height of sprung mass (Ls) 1.0 m Free height of unsprung mass (Lu) 0.3 mFigure 6
Quarter-car model.Figure7 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.Figure 7
Displacement convergence of the quarter car.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 Figure8. 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.Table 2
Iteration numbers of the quarter-car model.
Step size Iteration number CPU time (s) Ratio(%) Fixed h = 0.25 10 0.27 44 h = 0.15 19 0.32 52 h = 0.05 117 0.61 100  Variable h 0 = 0.05 8 0.26 42Figure 8
Step-size adaptation of the quarter-car model.The estimated natural frequencies for variable step size calculation are shown in Figure9. 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 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).Figure 9
Estimated natural frequency of the quarter-car model.Figure 10
Amplification factor of the quarter-car model.


## 4.3. 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 Figure11. 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 11
Deflection of the cable model.Figure12 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.Table 3
Iteration numbers of the simple cable model.
Step size Iteration number CPU time (s) Ratio(%) Fixed h = 0.5 8 3.7 28 h = 0.3 12 4.2 32 h = 0.1 64 12.8 100  Variable h 0 = 0.1 8 4.5 35Figure 12
Step-size adaptation of the simple cable model.Figure 13
Amplification factor of the simple cable model.


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



---

*Source: 4061069-2018-10-22.xml*