A fast and accurate dynamic relaxation scheme

Abstract

Dynamic relaxation method (DRM) is one of the suitable numerical procedures for nonlinear structural analysis. Adding the fictitious inertia and damping forces to the static equation, and turning it to the dynamic system, are the basis of this technique. Proper selection of the DRM artificial factors leads to the better convergence rate and efficient solutions. This study aims to increase the numerical stability, and to decrease the analysis time. To fulfil this objective, the reduction rate of analysis error for consecutive iterations is minimized. Based on this formulation, a new time step is found for the viscous dynamic relaxation. After combining this novel relationship with the other DRM factors, various geometrical nonlinear structures, such as trusses, frames, and shells, are analyzed. The obtained results verify the efficiency of authors’ scheme.

Keywords

viscous dynamic relaxation / time step / displacement error / geometric nonlinear analysis

Introduction

Nowadays, the analysis of structure, which exhibits large deflections, is of great importance for engineers. Because of economy, optimal material application, and also lightweight use, the designers are forced to plan very flexible structures. For these bodies, the equilibrium condition in the non-deformed position is not acceptable, since their shapes are continuously changing. It should be added that the analytical solution is rarely available for merely simple structures under rationalization certain conditions. Moreover, utilizing the superposition of effects is not allowed for structure with the large deformations, as it is applicable only for the linear behavior. As a result, numerical techniques with the incremental and iterative procedures are the main tools for analyzing complicated structures exhibiting large deformations. One of these schemes is the dynamic relaxation method (DRM).

In 1965, DRM was first used by Day [¹]. Nonlinear analysis of structures with this process was utilized by Brew and Brotton [²]. They proposed that the mass of each degree of freedom is a multiple of diagonal stiffness matrix element. In 1972, Bunce [³] estimated the critical damping for this procedure. Cassell and Hobbs [⁴] utilized Gerschgorin’s theory for the numerical stability of DRM. For the first time, error analysis of this scheme was employed by Papadrakakis [⁵]. He suggested an automatic process to calculate the DRM fictitious factors. In 1983, one of the famous formulations for DRM, was proposed by Underwood [⁶]. To obtain minimum frequency of the artificial dynamic system, he employed Rayleigh’s principle. Furthermore, this investigator employed a constant time step. Based on Rayleigh’s quotient, Qiang [⁷] offered equations for damping and time step. Munjiza et al. [⁸,⁹] obtained the damping matrix as \mathrm{2}\mathit{M}{({\mathit{M}}^{-\mathrm{1}}\mathit{S})}^{\mathrm{0.5}}. In 1989, modified DRM (maDR) to reach the steady-state response was introduced by Zhang and Yu [¹⁰]. In another investigation, Zhang et al. [¹¹] applied a novel technique to find fictitious damping factor. They supposed a similar factor for all degrees of freedoms of each node.

Barnes [¹²] analyzed tension structures by utilizing kinetic dynamic relaxation. It should be added, in the kinetic DRM, damping factor is equal to zero. Passing from a peak point of the kinetic energy diagram indicates that the kinetic energy has a reducing trend. Using the result of this time leads to divergence process. To remedy this shortcoming, Topping and Ivanyi [¹³] suggested that the displacement and velocity must be calculated at the mid-point. Many researches have been done to improve the DRM, so far. Rezaiee-Pajand and Taghavian-Hakkak [¹⁴] proposed a new formulation for DRM by utilizing the first three terms of Taylor’s series. An estimation for the time step via minimizing residual force was obtained by Kadkhodayan et al. [¹⁵]. Rezaiee-Pajand and Alamatian [¹⁶] studied some structures with geometric nonlinear analysis. In another study [¹⁷], they proposed new mass and damping matrices to ameliorate the performance of this process by using Gerschgorin’s theory.

In most cases, Rayleigh’s principle has been employed to estimate the damping. This principle is an upper bound for the minimum eigenvalues [³]. Rezaiee-Pajand and Sarafrazi [¹⁸] suggested the power iterative procedure to improve the convergence rate of DRM. In another investigation [¹⁹], they propounded the zero-damping technique to introduce a new time step. The comparison of different well-known methods has been done by Rezaiee-Pajand et al. [²⁰]. Displacement-based method for calculating the buckling load and tracing the post-buckling regions with DRM was offered by Alamatian [²¹]. Rezaiee-Pajand and Rezaee [²²] proposed a new formulation for determination of the time step in kinetic dynamic relaxation, based on the displacement error. Lately, mixing DRMs with the load factor and displacement increments was studied by Rezaiee-Pajand and Estiri [²³]. In another event, these researchers performed a study on finding equilibrium paths using DRM [²⁴]. Rezaiee-Pajand et al. [²⁵,²⁶] conducted a review of researches about the DRM.

In this study, by utilizing the displacement error in the analysis, a new time step is proposed for viscous DRM. By a combination of the suggested novel relationship with the other DRM factors, several solution algorithms are obtained. To demonstrate these strategies’ abilities, some miscellaneous structures with geometrically nonlinear behavior are analyzed. The obtained results verify the substantial reduction of the analysis time.

Dynamic relaxation

The following system of equation must be eventually solved for the structural static analysis:

(1)\mathit{S}\mathit{X}={\mathit{F}}_{\mathrm{r}}={\mathit{F}}_{\mathrm{r}},

where S, X, F_r and P_r are the stiffness matrix, displacement vector, internal and external load vectors in real condition, respectively. In DRM, artificial inertia and damping forces are added to Eq. (1). As a result, the static equilibrium converts to a fictitious dynamic system, as follows:

(2)\mathit{M}\ddot{\mathit{X}}+\mathit{C}\dot{\mathit{X}}+\mathit{S}\mathit{X}={\mathit{F}}_{\mathrm{f}},

where M, C, \dot{\mathit{X}} and \ddot{\mathit{X}} are fictitious mass matrix, damping matrix, velocity and acceleration vectors, respectively. Furthermore, F_f is the factitious internal load vector which was made by adding virtual inertia and damping forces (\mathit{M}\ddot{\mathit{X}} and \mathit{C}\dot{\mathit{X}}, respectively). Since DRM is an explicit scheme, it uses vector operator to reach the answer. To employ this numerical technique, pervious researchers have suggested different diagonal mass and damping matrices. It should be noted that the steady-state response of this artificial dynamic system can also satisfy the static equilibrium equation. To transfer the static governing equation to the dynamic space, inertia and damping forces should be assumed and added to static equation. In this way, the additional fictitious forces, which are known as the residual forces, are required to vanish. In other words, to have the static answer, the coming relation must be fulfilled for the dynamic system:

(3)\mathit{R}=\mathit{M}\ddot{\mathit{X}}+\mathit{C}\dot{\mathit{X}}={\mathit{P}}_{\mathrm{r}}-{\mathit{F}}_{\mathrm{f}},

where R is the residual force vector. Based on these assumptions and using the finite difference method, the fundamental iterative formulations for the DRM for the displacement and velocity are provided:

(4){\dot{x}}_i^{n+\frac{\mathrm{1}}{\mathrm{2}}}={\displaystyle \frac{\mathrm{2}-C_{ii}^n{h}^n}{\mathrm{2}+C_{ii}^n{h}^n}}{\dot{x}}_i^{n-\frac{\mathrm{1}}{\mathrm{2}}}-{\displaystyle \frac{\mathrm{2}{h}^n{r}_i^n}{m_{ii}^n(\mathrm{2}+C_{ii}^n{h}^n)}},\hspace{0.17em}i=\mathrm{1},\mathrm{2},\mathrm{...},ndof,

(5)x_i^{n+\mathrm{1}}=x_i^n+\Delta x_i^{n+\mathrm{1}}=x_i^n+h^{n+\mathrm{1}}{\dot{x}}_i^{n+\frac{\mathrm{1}}{\mathrm{2}}},\hspace{0.17em}i=\mathrm{1},\mathrm{2},\mathrm{...},ndof,

where superscripts n and i are number of the iteration and numerator of degrees of freedom, respectively. Furthermore, h and ndof represent the time step and the total number of degrees of freedom. According to the previous description, assuming appropriate values for fictitious factors of DRM is significantly effective on the solution speed, convergence rate and stability of the procedure. Some well-known assumptions to estimate these artificial factors are given in the following lines.

Fictitious factors

Different techniques have been proposed to calculate the fictitious parameters of DRM, so far. One of the most common approaches applied to the factitious diagonal mass matrix is suggested by Underwood [⁶]. This method which utilizes Gerschgorin’s circle theory has the next formula:

(6)m_{ii}={\displaystyle \frac{h^{\mathrm{2}}}{\mathrm{4}}}{\displaystyle \sum _{j=\mathrm{1}}^{ndof}\left|s_{ij}\right|}.

This study utilizes this equation in two ways, and they are called Method 1 and Method 2. Rezaiee-Pajand and Alamatian [¹⁷] introduced another well-known procedure for mass matrix. They minimized the displacement error in their process. The elements of this mass matrix are given as follows:

(7)m_{ii}=\max \left({\displaystyle \frac{h^{\mathrm{2}}}{\mathrm{2}}}{s}_{ii},{\displaystyle \frac{h^{\mathrm{2}}}{\mathrm{4}}}{\displaystyle \sum _{j=\mathrm{1}}^{ndof}\left|s_{ij}\right|}\right).

This matrix is applied in the third technique, and it is named Method 3. Using these mentioned two procedures, for estimating the factitious mass matrix, gives the better outcomes in the viscous dynamic relaxation, in comparison with other approaches [²⁰]. In this paper, the new proposed time step is combined with the abovementioned mass matrices. It is worth mentioning that the artificial dynamic system should be in critical damping state, for convergence to the answer in a shorter time. For this reason, many investigations have calculated the critical damping. In 1989, Zhang and Yu [¹⁰] suggested one of the simplest formulations for damping. This equation that is used in all mentioned three methods has the coming shape:

(8)C=\mathrm{2}\sqrt{{\displaystyle \frac{X^{T}F}{X^{T}MX}}}M.

Another fictitious parameter in DRM is the time step. Various studies have led to the appropriate value for this parameter. This paper attempts to improve the convergence rate and reduce the analysis duration by formulating a novel time step for the viscous DRM.

Initial displacement vector is another DRM parameter, which some values have assigned for it previously. After performing various numerical experiments, it was found that using these values for initial displacement is not very effective [²⁰]. As a result, the zero initial displacement vector is utilized for all three applied methods in this study.

Error analysis in viscous DRM

The displacement and velocity vectors have been shown by Eqs. (4) and (5), respectively. Utilizing these equations along with central finite difference, and substituting {\displaystyle \frac{x_i^n-x_i^{n-\mathrm{1}}}{h^n}} instead of {\dot{x}}_i^{n-\frac{1}{2}}, the following equation is obtained:

(9)x_i^{n+\mathrm{1}}=x_i^n+h^{n+\mathrm{1}}\left[{\displaystyle \frac{\mathrm{2}-C_{ii}^n{h}^n}{\mathrm{2}+C_{ii}^n{h}^n}}{\displaystyle \frac{x_i^n-x_i^{n-\mathrm{1}}}{h^n}}-{\displaystyle \frac{\mathrm{2}{h}^n{r}_i^n}{m_{ii}^n(\mathrm{2}+C_{ii}^n{h}^n)}}\right].

For each degree of freedom, Eq. (3) can be rewritten in the next form:

(10)r_i^n=p_i-f_i^n.

By substituting {\mathit{P}}_i with {\displaystyle \sum _{j=\mathrm{1}}^{ndof}{\mathit{S}}_{ij}^n}{\mathit{X}}_j^{*} and {\mathit{F}}_i^k with{\displaystyle \sum _{j=\mathrm{1}}^{ndof}{\mathit{S}}_{ij}^n}{\mathit{X}}_j^n, in which {\mathit{X}}_j^{\ast } is the exact displacement, in the last equation, the subsequent result can be found:

(11)r_i^n={\displaystyle \sum _{j=\mathrm{1}}^{ndof}{s}_{ij}^n}{x}_j^{*}-{\displaystyle \sum _{j=\mathrm{1}}^{ndof}{s}_{ij}^n}{x}_j^n.

Substituting Eq. (11) into Eq. (9) and simplifying, the following formula is achieved:

(12)x_i^{n+\mathrm{1}}=x_i^n+{\displaystyle \frac{h^{n+\mathrm{1}}}{h^n}}{\displaystyle \frac{\mathrm{2}-C_{ii}^n{h}^n}{\mathrm{2}+C_{ii}^n{h}^n}}(x_i^n-x_i^{n-\mathrm{1}})-{\displaystyle \frac{\mathrm{2}{h}^n{h}^{n+\mathrm{1}}}{m_{ii}^n(\mathrm{2}+m_{ii}^n{h}^n)}}\left({\displaystyle \sum _{j=\mathrm{1}}^{ndof}{s}_{ij}^n}{x}_j^{*}-{\displaystyle \sum _{j=\mathrm{1}}^{ndof}{s}_{ij}^n}{x}_j^n\right).

This equation can be rewritten in the below shape:

(13)x_i^{n+\mathrm{1}}-x_i^{*}=x_i^n-x_i^{*}+{\displaystyle \frac{h^{n+\mathrm{1}}}{h^n}}{\displaystyle \frac{\mathrm{2}-C_{ii}^n{h}^n}{\mathrm{2}+C_{ii}^n{h}^n}}(x_i^n-x_i^{*}-x_i^{n-\mathrm{1}}+x_i^{*})+{\displaystyle \frac{\mathrm{2}{h}^n{h}^{n+\mathrm{1}}}{m_{ii}^n(\mathrm{2}+C_{ii}^n{h}^n)}}{\displaystyle \sum _{j=\mathrm{1}}^{ndof}{s}_{ij}^n(x_j^n-x_j^{*})}.

By defining Eⁿ=Xⁿ−X*, as a vector showing the analysis error for the consecutive iterations, and assuming linear error reduction rate, r, for them, the following equations can be attained in vector space:

(14){\mathit{E}}^{n+\mathrm{1}}={\mathit{E}}^n+{\displaystyle \frac{h^{n+\mathrm{1}}}{h^n}}{\displaystyle \frac{\mathrm{2}-C_{ii}^n{h}^n}{\mathrm{2}+C_{ii}^n{h}^n}}({\mathit{E}}^n-{\mathit{E}}^{n-\mathrm{1}})+{\displaystyle \frac{\mathrm{2}{h}^n{h}^{n+\mathrm{1}}}{\mathrm{2}+C_{ii}^n{h}^n}}{\mathit{M}}^{-\mathrm{1}}{\displaystyle \sum _{j=\mathrm{1}}^{ndof}{s}_{ij}^n}({\mathit{E}}^n),

(15){\mathit{E}}^n=\rho {\mathit{E}}^{n-\mathrm{1}},

(16)\rho {\mathit{E}}^n=\left(\mathrm{1}+{\displaystyle \frac{h^{n+\mathrm{1}}}{h^n}}{\displaystyle \frac{\mathrm{2}-c_{ii}^n{h}^n}{\mathrm{2}+c_{ii}^n{h}^n}}\right){\mathit{E}}^n+{\displaystyle \frac{\mathrm{2}{h}^n{h}^{n+\mathrm{1}}{\mathit{M}}^{-\mathrm{1}}\mathit{S}}{\mathrm{2}+c_{ii}^n{h}^n}}{\mathit{E}}^n-{\displaystyle \frac{h^{n+\mathrm{1}}}{h^n}}{\displaystyle \frac{\mathrm{2}-c_{ii}^n{h}^n}{\mathrm{2}+c_{ii}^n{h}^n}}{\displaystyle \frac{{\mathit{E}}^n}{\rho }}.

These outcomes are simplified as follows:

(17)\alpha ={\displaystyle \frac{h^{n+\mathrm{1}}}{h^n}}{\displaystyle \frac{\mathrm{2}-C_{ii}^n{h}^n}{\mathrm{2}+C_{ii}^n{h}^n}},

(18)\beta =\mathrm{1}+\alpha ,

(19)\gamma ={\displaystyle \frac{\mathrm{2}{h}^n{h}^{n+\mathrm{1}}}{\mathrm{2}+C_{ii}^n{h}^n}},

(20)\mathit{B}={\mathit{M}}^{-\mathrm{1}}\mathit{S},

(21)\rho {\mathit{E}}^n=(\beta \mathit{I}+\gamma \mathit{B}){\mathit{E}}^n-{\displaystyle \frac{\alpha }{\rho }}{\mathit{E}}^n,

where I is identity matrix. Reordering Eq. (21) leads to the next equation:

(22)\left({\displaystyle \frac{{\rho }^{\mathrm{2}}-\beta \rho +\alpha }{\gamma \rho }}I-B\right)E^n=\mathrm{0.}

The eigenvalues of matrix B are determined by ({\lambda }^B\mathit{I}-\mathit{B}){\mathit{E}}^n=\mathrm{0}. Therefore, the next equation can be derived:

(23){\displaystyle \frac{{\rho }^{\mathrm{2}}-\beta \rho +\alpha }{\gamma \rho }}=-{\lambda }^B.

The last relation can be changed to the coming form:

(24){\rho }^{\mathrm{2}}-\left(\mathrm{1}+{\displaystyle \frac{h^{n+\mathrm{1}}}{h^n}}{\displaystyle \frac{\mathrm{2}-C_{ii}^n{h}^n}{\mathrm{2}+C_{ii}^n{h}^n}}-{\displaystyle \frac{\mathrm{2}{h}^n{h}^{n+\mathrm{1}}}{\mathrm{2}+C_{ii}^n{h}^n}}{\lambda }_i^{M^{-\mathrm{1}}S}\right)\rho +{\displaystyle \frac{h^{n+\mathrm{1}}}{h^n}}{\displaystyle \frac{\mathrm{2}-C_{ii}^n{h}^n}{\mathrm{2}+C_{ii}^n{h}^n}}=\mathrm{0.}

In this relation, {\lambda }_i^{M^{-\mathrm{1}}S} is the i-th eigenvalues of matrix {\mathit{M}}^{-\mathrm{1}}\mathit{S}. To achieve the maximum convergence, the ratio of error reduction must be minimal. Consequently, the discriminant of last equation must be set to zero.

(25)\mathrm{\Delta }=\mathrm{0}\Rightarrow {\left(\mathrm{1}+{\displaystyle \frac{h^{n+\mathrm{1}}}{h^n}}{\displaystyle \frac{\mathrm{2}-C_{ii}^n{h}^n}{\mathrm{2}+C_{ii}^n{h}^n}}-{\displaystyle \frac{\mathrm{2}{h}^n{h}^{n+\mathrm{1}}}{\mathrm{2}+C_{ii}^n{h}^n}}{\lambda }_i^{M^{-\mathrm{1}}S}\right)}^{\mathrm{2}}-\mathrm{4}{\displaystyle \frac{h^{n+\mathrm{1}}}{h^n}}{\displaystyle \frac{\mathrm{2}-C_{ii}^n{h}^n}{\mathrm{2}+C_{ii}^n{h}^n}}=\mathrm{0.}

To gain the time step in (n+1)-th iteration, the equation should be solved, which leads to the next formula:

(26)h^{n+\mathrm{1}}=-{\displaystyle \frac{(c{h}^n-\mathrm{2}{h}^{\mathrm{2}n}{\lambda }_i^{M^{-\mathrm{1}}S}-\mathrm{2}\pm \mathrm{2}\sqrt{\mathrm{2}{h}^{\mathrm{2}n}{\lambda }_i^{M^{-\mathrm{1}}S}(\mathrm{2}-c{h}^n)})h^n(\mathrm{2}+c{h}^n)}{{(\mathrm{2}{h}^{\mathrm{2}n}{\lambda }_i^{M^{-\mathrm{1}}S}-\mathrm{2}+c{h}^n)}^{\mathrm{2}}}}.

The i-th free oscillation frequency of the fictitious dynamic system is calculated from square root of {\lambda }_i^{M^{-\mathrm{1}}S}. Since the most effective response of the dynamic system is obtained from the minimum oscillation frequency, hence, {\lambda }_i^{M^{-\mathrm{1}}S} is used to achieve the structural response. Because the numerical instability is usually lower for smaller time steps, the small-time steps will be chosen from the following relationship:

(27)h^{k+\mathrm{1}}=-{\displaystyle \frac{(c{h}^k-\mathrm{2}{h}^{\mathrm{2}k}{\lambda }_{\mathrm{1}}^{M^{-\mathrm{1}}S}-\mathrm{2}\pm \mathrm{2}\sqrt{\mathrm{2}{h}^{\mathrm{2}k}{\lambda }_{\mathrm{1}}^{M^{-\mathrm{1}}S}(\mathrm{2}-c{h}^k)})h^k(\mathrm{2}+c{h}^k)}{{(\mathrm{2}{h}^{\mathrm{2}k}{\lambda }_{\mathrm{1}}^{M^{-\mathrm{1}}S}-\mathrm{2}+c{h}^k)}^{\mathrm{2}}}},

To determine the minimum frequency, the Zhang formula is used:

(28){\lambda }_{\mathrm{1}}^{{\mathit{M}}^{-\mathrm{1}}S}={\displaystyle \frac{{({\mathit{X}}^n)}^{\mathrm{T}}{\mathit{F}}^n}{{({\mathit{X}}^n)}^{\mathrm{T}}{\mathit{M}}^n{\mathit{X}}^n}}.

Numerical examples

To demonstrate the ability of the proposed procedure, a computer program in FORTRAN was written by the authors. Some different truss, frame and shell structures are solved by using this program, in this article. These analyses are based on the geometric nonlinear and elastic material behavior. All structures are solved by utilizing the three mentioned techniques. The first one is created by implementing the mass matrix and damping factor from Eqs. (6) and (4), respectively, together with constant time step. This technique is the benchmark method in this study, and it is also known as Zhang procedure. The method will be indicated as Method 1, hereinafter. In the second technique (Method 2), the time step of first method is replaced by the new suggested one, and it is called the proposed method, henceforth. The third scheme is obtained by substituting the mass matrix of Method 2 with the one which is calculated by Eq. (5). The name of this procedure is Method 3 and it is known as the combined technique, in this paper.

All of these three described methods are compared based on the analysis duration in each numerical example. It is found that the number of iteration for converging to the answers is the nearly same for all three procedures. Therefore, the comparison of this criterion is ignored. However, the number of iteration is reported in the results. It should be added that the response accuracies of all examples are the precisely same. To perform an impartial comparison, the residual force error is 10^‒6 for all trusses and frames, and 10^‒4 for shells. For all solutions, the forces are applied in ten steps. The load-displacement graph is illustrated for each structure. The convergence points in each step, are clearly demonstrated in all diagrams. To find the rank of each scheme, the analysis time criterion, E_T, is calculated from the next relation:

(29)E_T=\mathrm{100}{\displaystyle \frac{T_{\mathrm{1}}-T}{T_{\mathrm{1}}}},

where T is the analysis time to obtain the answers. Subscript 1 denotes the analysis time required to solve the benchmark method. Thus, E_T is equal to zero for this procedure. As a result, E_T is ignored for the benchmark technique. In the following lines, the numerical investigations are presented.

Arch truss

First numerical example is an arch truss. This two-dimensional structure, which is shown in Fig. 1, has 42 nodes and 81 members [²⁷]. Cross-section area and elastic modulus of the truss elements are 1 cm² and 205.8^{GPa (2.1×106 kg/cm2), respectively. The external Load P of 196 kN is applied at the highest structural node. The load-displacement curve of Node D is plotted in Fig. 2. Table 1 demonstrates the solution outcomes. Based on these results, Method 2 which is created by substituting constant time step in Zhang’s procedure with the new proposed one, is convergent to the near exact response, in a shorter time. Furthermore, Method 3 performs better than two other scheme, with regard to the analysis time. As mentioned previously, this approach is formed by the combination of the suggested time step, Zhang damping factor and mass matrix, which is presented by Rezaiee-Pajand and Alamatian [17]. It should be noted that number of iteration is not compared due to its proximity in all three procedures.}

33-member truss

Figure 3 shows the geometric of this structure. The truss consists of 18 nodes, 33 members and 32 degrees of freedom. The cross-section area and elastic modulus of this truss are the same and equal to 300 mm² and 300 GPa respectively. Also, the value of Load P is 400 kN. The static path of the load point, D, is illustrated in Fig. 4. Table 2 reports the outcomes of the analysis. As the results indicate, that the analysis duration time is reduced by 22.65% and 30.60% when using the Methods 2 and 3, respectively, compared to the Benchmark technique. Howbeit, comparison of the number of iteration was ignored in this study, but the results are shown a significant reduction in this term for the proposed and combined methods.

Truss tower

In the following, three-dimensional structures are solved. The first one is the truss tower which is illustrated in Fig. 5. This structure has 24 nodes and 90 members [¹⁹]. The elastic modulus is 205.8 GPa (2.1×10⁶ kg/cm²) and cross-section area of each member is 1 cm². The value of Load P is 294 kN (30000 kg). Figure 6 demonstrates the load-displacement curve for Node D. Table 3 represents the numerical outcomes. Based on the presented results in Table 3, the proposed method and the combined one decrease the analysis time by 12.54% and 20.28%, respectively. However, the total numbers of iteration are almost close together in all three methods.

Shallow truss

This structure which was first analyzed by Rezaiee-Pajand and Alamatian [²⁸], has extremely nonlinear behavior. The truss consists of 73 nodes, 168 elements and 147 degrees of freedom. Other geometric characteristics are displayed in Fig. 7. The cross-section area and elastic modulus of all members are 1 mm² and 100 GPa, respectively. The tip node is subjected to two different vertically concentrated loads of 1000 and 1200 N. Also, the lowest nodes are set on the simple support. Figure 8 shows the static equilibrium path of the top node. The number of iteration and analysis time of 1000 N load are inserted in Table 4. According to Table 4, despite the equal number of iteration for all three procedures, the analysis duration is diminished by about 6.54% and 9.03% for the proposed and combined method, respectively.

Circular truss dome

This example has 216 nodes, 16 elements and 576 degrees of freedom [²⁹]. As it is shown in Fig. 9, the structure is formed by 24 similar pieces. According to Fig. 9, each segment covers 15° of the truss. The model of underside nodes is all pinged support. Besides, the load of 35 kN is applied to all other nodes in the opposite direction of Z-axis. The cross-sectional area of members is 2×10^‒3 m² and the modulus of elasticity is 200 GPa. The load-displacement curve of nodes on the highest level is depicted in Fig. 10. The number of iterations and analysis time of all processes are arranged in Table 5. As the outcomes represent, decrement of the analysis duration is 8.65% for the proposed method and 14.65% for combined technique.

Portal frame

The first frame example which is solved in this section is demonstrated in Fig. 11 [²⁰]. All structure members are divided by three fame element, with elastic modulus of 8.2740379 MPa cross-section area of 12.9032×10^‒4 m² and inertia moment of 27.7487×10^‒8 m⁴. Moreover, the Load P is equal to 0.44484 N. Figure 12 illustrates the static equilibrium path of horizontal displacement of Node D versus the Load P. After solving this example, the obtained results are inserted in Table 6. Table 6 shows that using the proposed and combined methods leads to lessen the analysis time. Once again, the number of iterations is nearly the same for all three procedures.

Five-story frame

Figure 13 depicts the geometry of this frame. The structure is constructed by the columns of W2150 and beams of W1835 [²⁰]. The elastic modulus of the material is 196 GPa. Figure 13 clearly demonstrates the structural applied load. The load-displacement curve of Node u is drawn in Fig. 14. Table 7 presents the performance index of all approaches. According to Table 7, the superiority of proposed and combined processes against the benchmark one is evident. The main advantage is the decrement in the analysis time. This criterion is reduced by 24.46% for proposed technique and 34.87% for combined procedures.

Toggle frame

The other example is that of a two members frame which was analyzed having nonlinear behavior by many other researchers [²⁸]. The structure is shown in Fig. 15. Each member is modeled by five frame elements. The elastic modulus, cross-sectional area and inertia moment of the members are: E = 71 GPa, A = 118.06 mm² and I = 374.61 mm⁴, respectively. The value of concentrated load is 44.28 N. The static equilibrium path of the Node D is plotted in Fig. 16. All comparison results are inserted in Table 8. Based on the analysis duration, compared with the benchmark procedure, it is evident in this table that the Methods 2 and 3 are more appropriate. According to the recorded outcomes, these two processes shorten the time of analysis as much as 13.88% and 13.21%, respectively.

Cylindrical roof

In the following, the abilities of the propounded approaches are examined in the shell analysis. The first example is a cylindrical roof, which is displayed in Fig. 17 [³⁰]. The shell structural properties are as follows: Elastic modulus= 3.10275 GPa, thickness= 6.35 mm, and Poisson’s ratio= 0.3. As it is shown in Fig. 17, both curved edges of the shell are free and the others are on the simple supports. To solve the structure, it was divided into four equal sections which each of them is modeled by 100 shell elements. The concentrated load of 750 N is exerted on the Node A. The vertical displacement of this node versus the applied load is illustrated in Fig. 18. As it can be seen in this figure, the displacement between Steps 7 and 8 was greatly increased. This is due to the existence of load limit points. To indicate these points, the load steps must be shortened between main Steps 7 and 8. The load limit points acquired in this way, are included in Fig. 19. Table 9 presents the outcomes. Table 9 reveals that the proposed and combined methods reduce the analysis time by 11.79% and 23.53%, respectively. However, the number of iterations is almost the same for all three techniques. Furthermore, the deformed shape of the structure under this loading condition is depicted in Fig. 20.

Cantilever plate

The configuration of this structure is shown in Fig. 21 [³¹]. To analyze this shell, 96 elements are used. The modulus of elasticity, thickness, and Poisson’s ratio are 1.2 GPa, 4 mm, and 0.3, respectively. The Load P of 35 N is extorted on the Node A. Figure 22 illustrates the load-displacement diagram for the Node B. The required iterations, and the analysis time for each method are given in Table 10. The excellence of both Methods 2 and 3 is deduced from this table. These techniques reduce the time to obtain the answer as much as 11.45% and 16.10%, respectively.

Cook membrane

The other example which is studied in this section is the cook structural membrane [³²]. The geometry of this structure is depicted in Fig. 23. This membrane is modeled by 322 elements. The following data are used in this example: The modulus of elasticity= 1, thickness= 1, and Poisson’s ratio= 0.3. All of these values are dimensionless. Furthermore, Load P is equal to 1. The static equilibrium path of the marked point is shown in Fig. 24. Table 11 indicates the obtained results. According to the data, which is presented in Table 11, the analysis duration is decreased by 10.24% for the proposed method and 15.48% for the combined procedure, in comparison with the benchmark technique. Nevertheless, the number of iteration is again nearly equal for all three processes.

Shallow spherical cap

Figure 25 shows a shallow spherical cap shell. All boundaries are simply supported [³³]. Because of symmetry, each quarter of the structure is modeled by 200 elements. The elastic modulus, thickness, and Poisson’s ratio are 68.95 GPa, 9.95 cm, and 0.3, respectively. The value of Load P, which applied on the center of structure, is 22241 N. Figure 26 clearly demonstrates the relationship between the load and displacement of the tip node. The numbers of iteration and analysis durations of three procedures are inserted in Table 12. Based on the information contained in Table 12, Methods 2 and 3 reduce the analysis time by 11.46% and 14.98%, respectively. Once again, since the number of iteration is almost the same for all three techniques, the comparison of this criterion is ignored.

Conclusions

DRM is one of the efficient numerical techniques for nonlinear analysis. Being explicit, by using only the vector operations, and requiring low memory, are the rationales for this superiority. In this approach, the fictitious dynamic system is made by adding the artificial inertia and damping forces to the static equilibrium equations. The steady-state response of this fake dynamic space is the same as the static answer of the original problem. It is important to notice that the well estimation of values for each DRM dummy factors significantly affects the convergence rate, accuracy of the response, and stability of the solution. Based on these facts, a novel time step is suggested by utilizing the displacement error in this study. The goal is to reduce the analysis time and better numerical convergence. In this regard, combinations of the proposed relationship with the others' previous propounded factors lead to the creation of two procedures, which are named proposed and combined methods. To determine the reduction percentage in the analysis time, various structures, such as, trusses, frames and shells are solved. Since the number of iteration is almost identical in all approaches, the comparison of this criterion is ignored. However, the times of obtaining the responses are different. Findings indicate that the proposed, and the combined techniques decrease the analysis time on average of 13.36% and 19.16%, respectively.

References

[1] Day A S. An introduction to dynamic relaxation. Engineer, 1965, 219: 218–221

[2] Brew J S, Brotton D M. Non-linear structural analysis by dynamic relaxation. International Journal for Numerical Methods in Engineering, 1971, 3(4): 463–483

[3] Bunce J W. A note on the estimation of critical damping in dynamic relaxation. International Journal for Numerical Methods in Engineering, 1972, 4(2): 301–303

[4] Cassell A C, Hobbs R E. Numerical stability of dynamic relaxation analysis of non-linear structures. International Journal for Numerical Methods in Engineering, 1976, 10(6): 1407–1410

[5] Papadrakakis M. A method for the automatic evaluation of the dynamic relaxation parameters. Computer Methods in Applied Mechanics and Engineering, 1981, 25(1): 35–48

[6] Underwood P.Dynamic relaxation. Computational Methods for Transient Analysis, 1983, 245–265

[7] Qiang S. An adaptive dynamic relaxation method for nonlinear problems. Computers & Structures, 1988, 30(4): 855–859

[8] Munjiza A. A K_m proportional damping for dynamic relaxation. International Journal for Engineering Modelling, 1996, 9: 1–9

[9] Munjiza A, Owen D R J, Crook A J L, AnM(M⁻¹K)Mproportional damping in explicit integration of dynamic structural systems. International Journal for Numerical Methods in Engineering, 1998, 41(7): 1277–1296 doi:10.1002/(SICI)1097-0207(19980415)41:7<1277::AID-NME335>3.0.CO;2-9

[10] Zhang L G, Yu T X. Modified adaptive dynamic relaxation method and its application to elastic-plastic bending and wrinkling of circular plates. Computers & Structures, 1989, 33(2): 609–614

[11] Zhang L C, Kadkhodayan M, Mai Y W. Development of the maDR method. Computers & Structure, 1994, 52(1): 1–8

[12] Barnes M. Form finding and analysis of tension structures by dynamic relaxation. International Journal of Space Structures, 1999, 14(2): 89–104

[13] Topping B H V, Ivanyi P. Computer Aided Design of Cable Membrane Structures.Stirling: Saxe-Coburg Publications, 2008

[14] Rezaiee-Pajand M, Taghavian Hakkak M. Nonlinear analysis of truss structures using dynamic relaxation. International Journal of Engineering, 2006, 19(1): 11–22

[15] Kadkhodayan M, Alamatian J, Turvey G J. A new fictitious time for the dynamic relaxation (DXDR) method. International Journal for Numerical Methods in Engineering, 2008, 74(6): 996–1018

[16] Rezaiee-Pajand M, Alamatian J. Nonlinear dynamic analysis by dynamic relaxation method. Structural Engineering and Mechanics, 2008, 28(5): 549–570

[17] Rezaiee-Pajand M, Alamatian J. The dynamic relaxation method using new formulation for fictitious mass and damping. Structural Engineering and Mechanics, 2010, 34(1): 109–133

[18] Rezaiee-Pajand M, Sarafrazi S R. Nonlinear structural analysis using dynamic relaxation method with improved convergence rate. International Journal of Computational Methods, 2010, 7(4): 627–654

[19] Rezaiee-Pajand M, Sarafrazi S R. Nonlinear dynamic structural analysis using dynamic relaxation with zero damping. Computers & Structures, 2011, 89(13–14): 1274–1285

[20] Rezaiee-Pajand M, Sarafrazi S R, Rezaiee H. Efficiency of dynamic relaxation methods in nonlinear analysis of truss and frame structures. Computers & Structures, 2012, 112‒113: 295–310

[21] Alamatian J. Displacement-based methods for calculating the buckling load and tracing the post-buckling regions with dynamic relaxation method. Computers & Structures, 2013, 114–115: 84–97

[22] Rezaiee-Pajand M, Rezaee H. Fictitious time step for the kinetic dynamic relaxation method. Mechanics of Advanced Materials and Structures, 2014, 21(8): 631–644

[23] Rezaiee-Pajand M, Estiri H. Mixing dynamic relaxation method with load factor and displacement increments. Computers & Structures, 2016, 168: 78–91

[24] Rezaiee-Pajand M, Estiri H. Finding equilibrium paths by minimizing external work in dynamic relaxation method. Applied Mathematical Modelling, 2016, 40(23–24): 10300–10322

[25] Rezaiee-Pajand M, Alamatian J, Rezaee H. The state of the art in dynamic relaxation methods for structural mechanics Part 1: Formulations. Iranian Journal of Numerical Analysis and Optimization, 2017, 7(2): 65–86

[26] Rezaiee-Pajand M, Alamatian J, Rezaee H. The state of the art in dynamic relaxation methods for structural mechanics Part 2: Applications. Iranian Journal of Numerical Analysis and Optimization, 2017, 7(2): 87–114

[27] Rezaee H. Nonlinear structural analysis using dynamic relaxation method. Dissertation for the Doctoral Degree. Mashhad: Ferdowsi University, 2012 (in Persian)

[28] Rezaiee-Pajand M, Alamatian J. Automatic DR structural analysis of snap-through and snap-back using optimized load increments. Journal of Structural Engineering, 2011, 137(1): 109–116

[29] Alamatian J. A new formulation for fictitious mass of the dynamic relaxation method with kinetic damping. Computers & Structures, 2012, 90‒91: 42–54

[30] Rezaiee-Pajand M, Estiri H. Computing the structural buckling limit load by using dynamic relaxation method. International Journal of Non-linear Mechanics, 2016, 81: 245–260

[31] Kuo-Mo H. Nonlinear analysis of general shell structures by flat triangular shell element. Computers & Structures, 1987, 25(5): 665–675

[32] Kuznetsov V V, Levyakov S V. Phenomenological invariant-based finite-element model for geometrically nonlinear analysis of thin shells. Computer Methods in Applied Mechanics and Engineering, 196(49): 4952–4964

[33] Boutagouga D, Gouasmia A, Djeghaba K. Geometrically nonlinear analysis of thin shell by a quadrilateral finite element with in-plane rotational degrees of freedom. European Journal of Computational Mechanics, 2012, 19(8): 707–724 doi:10.3166/ejcm.19.707-724

RIGHTS & PERMISSIONS

Higher Education Press and Springer-Verlag GmbH Germany, part of Springer Nature