A fast and accurate dynamic relaxation approach for form-finding and analysis of bending-active structures

Abstract

Active bending is an increasingly popular construction technique that uses elastically bent structural members to form complex curved shapes. The design and analysis of bending-active structures requires an accurate simulation of the bending process, which is often complicated by the occurrence of large displacements. In this article, we propose to combine a previously developed implicit dynamic relaxation method with co-rotational beam elements to obtain a fast and accurate method for form-finding and analysis of bending-active structures. This approach is applied to four test cases. Implicit dynamic relaxation is compared to the classic Newton–Raphson method and conventional dynamic relaxation. The results show that the proposed implicit dynamic relaxation approach can be stabilized intuitively by changing the time step and damping ratio, making it more stable than the classic Newton–Raphson method. Moreover, the proposed approach converges fast compared to the conventional dynamic relaxation: the total computation time is considerably lower, even though the computation time per iteration is higher. Finally, a high accuracy is achieved due to the use of co-rotational beam elements. The combination of high accuracy and low computation time makes this approach well-suited for both form-finding and analysis of bending-active structures.

Keywords

active bending dynamic relaxation fictitious mass 6-degree-of-freedom beam element

Introduction

Active bending is an increasingly popular construction technique to build complex curved shapes in a simple and economical way. Lienhard² gives the following definition: “Bending-active structures are structural systems that include curved beam or shell elements which base their geometry on the elastic deformation from an initially straight or planar configuration.” Designing such structures requires an accurate simulation of the bending process. This is a challenging problem, as bending-active structures usually undergo very large displacements, making it necessary to use nonlinear solution techniques.³ The simulation process is often regarded as a form-finding problem. As simulating the bending process requires the material and section properties of the laths, only form-finding techniques that take these parameters into account can be used. This limits the choice of form-finding methods to stiffness matrix methods, such as the Newton–Raphson method, or dynamic equilibrium methods, such as dynamic relaxation.⁴

Dynamic relaxation was first proposed by Day⁵ as an explicit solver for the static analysis of structures. The method quickly became popular for the form-finding of tension/compression structures such as (grid) shells, cable nets, membranes, and tensegrity structures, with important contributions from Barnes,⁶ Topping,⁷ Papadrakakis,⁸ and Wakefield.⁹ The method solves a static problem by reformulating it as a dynamic problem and solving the reformulated problem iteratively. In its initial position, the internal forces of a structure are not in equilibrium with the externally applied loads. In the dynamic relaxation solution scheme, the structure starts to oscillate because of the lack of equilibrium, until, under influence of damping, it finally comes to rest in the sought equilibrium position. The movement of the structure is traced time step by time step, where each time step corresponds to one iteration. The number of time steps needed to reach equilibrium depends on the inertia/mass of the structure, the damping, and the length of the time steps. These parameters only affect the structure’s dynamic behavior and can therefore be considered as tuning parameters, allowing control over the number of iterations, without affecting the final result. Their value can be chosen freely to minimize the number of iterations, as long as the numerical stability of the method is maintained. For problems, including only tension/compression elements, effective heuristic rules exist to choose the time step and fictitious masses in such a way that the stability is ensured.^10,11

The recent popularity of bending-active structures has led to an increasing interest to include the effects of bending and torsion in the dynamic relaxation process. The most straightforward approach is to use three-dimensional beam elements with 6 degrees of freedom (6 DoF) per node.^12–15 The inclusion of rotational DoF, however, leads to slow convergence when a lumped mass matrix is used. Or, as Adriaenssens and Barnes¹⁶ formulate it, “… it is often the coupling of these (rotational DoF) with axial stiffness and translational DoF, which cause conditioning problems in numerical explicit methods such as dynamic relaxation.” To avoid this issue, several authors have tried reducing the number of DoF needed in the calculation. The effects of bending can be taken into account using only 3 DoF per node by replacing bending moments with equivalent translational shear forces.^11,12,17 While this approach was initially limited to torsion-free rods with so-called isotropic cross section (i.e. cross sections with an identical second moment of area around the cross section’s two principal axes), Barnes et al.¹⁸ attempted to include torsion without increasing the number of DoF. Apart from other limitations, the biggest concern with this method was the need for a so-called torsion factor to correct for the lack of lateral bending stiffness. Adding a degree-of-freedom per node is therefore inevitable to fully account for torsion and anisotropic cross sections.^19–21 However, beam elements with 4 DoF per node still have some important drawbacks: the accuracy of a 4-DoF approach is significantly worse than a 6-DoF approach, and the approach requires more elements for a comparable accuracy, resulting in longer computation time.²⁰

Several authors have suggested ways to overcome the slow convergence associated with the 6-DoF approach. Papadrakakis¹⁰ proposed a method for the automatic evaluation of the dynamic relaxation parameters through an approximation of the highest and the lowest eigenfrequency. Similarly, Underwood²² described an adaptive dynamic relaxation approach. This approach was modified and improved by Zhang and colleagues.^23,24 Kadkhodayan et al.²⁵ proposed a way to update the time step by minimizing the residual force after each iteration. Rezaiee-Pajand and colleagues^26,27 introduced new relations for the fictitious mass and damping matrix, allowing the time step to change each iteration. Finally, Alamatian²⁸ presented new fictitious masses for a varying time step in combination with kinetic damping. However, all these methods employ lumped mass matrices, which limit their effectiveness in speeding up the convergence.

In this article, we propose a dynamic relaxation approach that can be efficiently used for form-finding and analysis of bending-active structures. For this purpose, the implicit dynamic relaxation method presented in a previous paper²⁹ is combined with 6-DoF co-rotational beam elements, as these elements are particularly well suited to handle large displacements.³⁰

The remainder of this article is organized as follows. The following section recapitulates the implicit dynamic relaxation method developed by Rombouts et al.²⁹ Next, co-rotational beam elements are briefly explained. Subsequently, the performance of the proposed approach is demonstrated for form-finding and analysis of bending-active structures through four representative test cases. Finally, a conclusion is given.

Implicit dynamic relaxation

In the dynamic relaxation approach, the static problem is reformulated as a dynamic problem. In the linear case, and assuming a finite element discretization in space, the dynamic behavior of a structure is governed by the following equation

M a + C v + K u = p

(1)

where M is the mass matrix, a is the acceleration vector, C is the damping matrix, v is the velocity vector, K is the stiffness matrix, u is the displacement vector, and p is the external load vector. Note that when the velocity and the acceleration vectors are zero, equation (1) describes the linear static equilibrium. The full dynamic equilibrium, however, is a set of coupled differential equations, as the velocities v are the time derivatives of the displacements u, and the accelerations a are in turn the time derivatives of the velocities v. The system of equations is solved for the displacements u by discretization in time using a central difference scheme. At a given moment in time, the velocity vector v is written in terms of the displacement vector u as follows

v^{t + Δ t / 2} = \frac{u^{t + Δ t} - u^{t}}{Δ t}

(2)

where $Δ t$ is the time step, $v^{t + Δ t / 2}$ is the velocity vector at time $t + Δ t / 2$ , $u^{t + Δ t}$ is the displacement vector at a time $t + Δ t$ , and $u^{t}$ is the displacement vector at time t. Note that the velocities are defined at intermediate time intervals $(Δ t / 2)$ when compared to the displacements. The velocities at a given time step are then simply obtained by linear interpolation

v^{t} = \frac{v^{t + Δ t / 2} + v^{t - Δ t / 2}}{2}

(3)

Similarly, the acceleration vector is discretized

a^{t} = \frac{v^{t + Δ t / 2} - v^{t - Δ t / 2}}{Δ t}

(4)

Equation (1) is rewritten as

{Ma}^{t} + {Cv}^{t} = r^{t}

(5)

where $r^{t}$ is the residual force vector at time t, defined as

r^{t} = p - {Ku}^{t}

(6)

For nonlinear problems, this linear relation does not hold, and equation (6) is generalized as

r^{t} = p - f^{t}

(7)

where $f^{t}$ is the internal force vector at time t.

Given time step $Δ t$ , damping matrix C, and mass matrix M, the displacements u are traced time step by time step until equilibrium is found under the influence of damping. Initially, displacements u, velocities v, and accelerations a are set to zero or some other arbitrary value. By combining equations (3), (4), and (5), the following expression is found to update the velocities v

v^{t + Δ t / 2} = {(\frac{M}{Δ t} + \frac{C}{2})}^{- 1} (\frac{M}{Δ t} - \frac{C}{2}) v^{t - Δ t / 2} + {(\frac{M}{Δ t} + \frac{C}{2})}^{- 1} r^{t}

(8)

Equation (2) is used to update the displacements u, based on the velocities found using equation (8). When after a number of time steps, the residual forces r become sufficiently small, the system is said to have converged. Alternatively, convergence is assumed if the velocities are below a certain threshold.

In order to limit the calculation time, the time step $Δ t$ , damping matrix C, and mass matrix M should be chosen in such a way as to minimize the number of iterations needed for the system to converge. However, an incorrect choice will result in numerical instability and will consequently divergence.

The stability limit for the time step is determined by the highest eigenfrequency of the discretized structure (see further, equation (13)). A larger time step is possible when the highest eigenfrequency is lower. Because a larger time step implies that less iterations are needed to reach equilibrium, the highest eigenfrequency should be minimized. On the other hand, since lower eigenmodes attenuate more slowly, the lowest eigenfrequencies should be maximized. Therefore, optimal convergence can be expected when all eigenfrequencies coincide. The eigenfrequencies can be manipulated by the choice of the fictitious mass, which has no impact on the final equilibrium position, as pointed out before.

The fictitious mass that makes all eigenfrequencies coincide is determined as follows. If the dynamic behavior of the structure is governed by equation (1), the modes and eigenfrequencies are found by analyzing the free vibration of the system. This means that no load is applied, and damping is ignored

Ma + Ku = 0

(9)

Non-trivial solutions to equation (9) are found by solving the following generalized eigenvalue problem

K ϕ - ω^{2} M ϕ = 0

(10)

where $ϕ$ is an eigenvector and $ω$ is an eigenfrequency. Instead of finding the eigenfrequencies for given mass and stiffness matrices, a fictitious mass matrix is sought that makes all eigenfrequencies coincide. Equation (10) shows that this is the case if the mass matrix is proportional to the stiffness matrix. This can be demonstrated as follows. Assume that the following equality holds

M = \frac{1}{ω^{*^{2}}} K

(11)

with $ω^{*}$ being the target eigenfrequency for all modes. In this case, equation (10) can be written as

K ϕ - \frac{ω^{2}}{ω^{*^{2}}} K ϕ = 0

(12)

This equation has non-trivial solutions if and only if $ω = ω^{*}$ . Therefore, equation (11) is used to choose the fictitious mass matrix.

The time step for which the dynamic relaxation method remains numerically stable follows from a stability analysis, as described by Bathe³¹ or Rombouts et al.²⁹ If all eigenfrequencies coincide, the numerical stability limit is given by

Δ t \leq \frac{2}{ω^{*}}

(13)

Note that, for the derivation of the optimal fictitious mass, the structural response was assumed to be linear. In a nonlinear case, the stiffness matrix should be updated at each iteration. Moreover, the structure’s eigenfrequencies will slightly shift, but if the time step is chosen somewhat below the limit in equation (13), fast and stable convergence is still obtained.

Next, appropriate values for the damping parameters must be determined. In order to critically damp a structure using viscous damping, the eigenfrequencies of the structure should be known, which usually requires a trial run or modal analysis. Instead, an artificial damping technique, so-called kinetic damping, is often applied. This technique resets the velocities to zero whenever the structure encounters a kinetic energy peak. In the implicit dynamic relaxation approach, all eigenfrequencies are known in advance, and viscous damping can be applied without requiring any additional computations. It can be shown that in the case of coinciding eigenfrequencies and classical damping, the damping matrix is proportional to the mass matrix

C = 2 ξ^{*} ω^{*} M

(14)

where $ξ^{*}$ is the damping ratio, assumed to be identical for all modes.

Finally, an updating scheme for the displacements is obtained by combining equations (2), (8), (11), and (14)

u^{t + Δ t} = u^{t} + Δ t (\frac{1 - Δ t ξ^{*} ω^{*}}{1 + Δ t ξ^{*} ω^{*}} v^{t - Δ t / 2} + \frac{ω^{*^{2}} Δ t}{1 + Δ t ξ^{*} ω^{*}} K^{- 1} r^{t})

(15)

This expression relates the displacements of the next time step to displacements, velocities, and residual forces from the current time step. Assume that the time step and the damping ratio are chosen as follows

\begin{matrix} Δ t = \frac{\sqrt{2}}{ω^{*}} \\ ξ^{*} = \frac{1}{\sqrt{2}} \end{matrix}

(16)

In this case, equation (15) simplifies to the well-known Newton–Raphson iteration scheme. (Equation (17) corresponds to the classic (force-controlled) Newton–Raphson method. Related solution strategies such as displacement control or arc-length control, which are primarily used to trace full load–displacement paths, are not discussed in this article.)

u^{t + Δ t} = u^{t} + K^{- 1} r^{t}

(17)

Note that the damping ratio $ξ^{*}$ adopted in equation (16) is slightly lower than the critical damping ratio of 1, and the time step $Δ t$ is within the stability limits for linear problems (equation (13)). The choice for the time step and damping ratio directly affects the method’s convergence rate and numerical stability. A larger time step combined with a lower damping ratio generally results in faster, but less stable convergence, whereas a smaller time step combined with a higher damping ratio results in slower, but more stable convergence.

Because the proposed method requires the inversion of a global stiffness matrix, one iteration in the current scheme will require more calculation time compared to one iteration of classic dynamic relaxation methods. However, the considerable reduction in the number of iterations makes up for the increased iteration time. This will be discussed in more detail in the examples section.

Co-rotational formulation

Co-rotational beam elements are specifically developed to handle large displacements.³⁰ They are used in this article to describe the behavior of bending-active structures, which often undergo large displacements. Co-rotational elements are characterized by the decoupling of rigid body motion and local beam deformations. As a result, arbitrarily large translations and rotations are allowed, as long as the strains remain small. If this is not the case, a finer spatial discretization is required to maintain accuracy. Felippa and Haugen³² give a general explanation of the co-rotational formulation and compare it to the total and updated Lagrangian formulation. In the remainder of this section, the basis of the co-rotational beam element as described by Crisfield³³ is briefly recapitulated. A more extensive discussion can be found in Crisfield.^33,34 The considered beam element is based on the Euler–Bernoulli beam theory. Furthermore, the material is assumed to be homogeneous, isotropic, and linearly elastic.

To distinguish between rigid body motion and local beam deformations, a local co-rotated reference configuration is fitted to the actual deformed beam. Local deformations are traced by comparing the actual deformed beam to the co-rotated reference beam. The co-rotated reference beam is represented by a set of local element basis vectors $T_{e} = [e_{x} e_{y} e_{z}]$ , where e_x lies in the direction of the vector connecting the beam ends, representing the centerline of the reference beam, and e_y and e_z describe the direction of reference beam’s principal axis. On the other hand, the current configuration of the actual deformed beam is represented by nodal basis vectors $T_{a} = [a_{x} a_{y} a_{z}]$ and $T_{b} = [b_{x} b_{y} b_{z}]$ , which describe the orientation of the beam ends (Figure 1). Similar to the element-basis vectors, a_x and b_x point in the direction of the beam’s centerline, while a_y and a_z, and b_y and b_z point in the direction of the cross section’s principal axes. For an undeformed beam, both sets of nodal basis vectors coincide with the element basis vectors.

Figure 1.

Definition of the nodal basis vectors $(a_{x}, a_{y}, a_{z})$ and $(b_{x}, b_{y}, b_{z})$ as well as the element basis vectors $(e_{x}, e_{y}, e_{z})$ of a deformed beam element.

In each time step, the basis vectors are updated using Rodrigues’ rotation formula.^32,34 Equation (18) gives the relation between the current nodal basis vectors and the previous ones, formulated for beam end “a.” The calculation of the nodal basis vectors of beam end “b” is analogous

T_{a}^{t} = R (Δ λ_{a}) T_{a}^{t - Δ t}

(18)

Here, $T_{a}^{t}$ is the nodal basis matrix of beam end “a” at time t, with the nodal basis vectors as columns, $T_{a}^{t - Δ t}$ is the nodal basis matrix at time $t - Δ t$ , and $R (λ)$ is an orthogonal rotation matrix describing the rotation around a given pseudo-vector $λ = {α β γ}^{T}$ over a given angle $λ = ‖ λ ‖$ , where the double vertical lines denote the Euclidean norm

R (λ) = I + \frac{\sin (λ)}{λ} S (λ) + \frac{1 - \cos (λ)}{λ^{2}} S (λ) S (λ)

(19)

where I is a three-by-three identity matrix, and S is a skew-symmetric matrix defined as

S (λ) = [\begin{matrix} 0 & - γ & β \\ γ & 0 & - α \\ - β & α & 0 \end{matrix}]

(20)

In equation (18), the rotation vector $Δ λ_{a} = [Δ α_{a} Δ β_{a} Δ γ_{a}]^{T}$ describes the change in orientation of beam end “a” between this iteration and the previous one. The element basis vectors $T_{e}^{t}$ are calculated based on the current nodal basis vectors stored in $T_{a}^{t}$ and $T_{b}^{t}$ . The actual updating procedure of the element basis vectors can be found in the original paper.³³

Next, the deformation of the beam is described by seven independent local beam deformations, recorded at the beam ends. These deformations are referred to as “strains” in the original paper. However, because the concerning quantities cannot be regarded as strains in the strict mechanical sense, the term “deformations” is used here instead. An axial deformation $δ$ gives rise to a normal force N, and six rotational deformations $θ$ result in six corresponding bending moments M. These rotational deformations describe the beam’s twisting as well as bending around the cross section’s two principal axes at the beam ends.

The local rotational deformations $θ$ are expressed in terms of the nodal and element basis vectors as follows

\begin{matrix} 2 \sin (θ_{x, a}) = - {(a_{z})}^{T} e_{y} + {(e_{z})}^{T} a_{y} \\ 2 \sin (θ_{x, b}) = - {(b_{z})}^{T} e_{y} + {(e_{z})}^{T} b_{y} \\ 2 \sin (θ_{y, a}) = - {(a_{z})}^{T} e_{x} + {(e_{z})}^{T} a_{x} \\ 2 \sin (θ_{y, b}) = - {(b_{z})}^{T} e_{x} + {(e_{z})}^{T} b_{x} \\ 2 \sin (θ_{z, a}) = - {(a_{y})}^{T} e_{x} + {(e_{y})}^{T} a_{x} \\ 2 \sin (θ_{z, b}) = - {(b_{y})}^{T} e_{x} + {(e_{y})}^{T} b_{x} \end{matrix}

(21)

where $θ_{x}$ , $θ_{y}$ , and $θ_{z}$ are the deformation angles around the local x-, y-, and z-axes, respectively (Figure 2). The sine functions can be omitted for small rotational deformations. The derivation of equation (21) can be found in Crisfield.³⁴

Figure 2.

Orientation of the local beam deformations $(δ θ_{x, a} θ_{y, a} θ_{z, a} θ_{x, b} θ_{y, b} θ_{z, b})$ and corresponding internal forces $(N M_{x, a} M_{y, a} M_{z, a} M_{x, b} M_{y, b} M_{z, b})$ .

The axial deformation $δ$ is simply the difference of the Euclidian distance between beam ends “a” and “b” at the current time $l^{t}$ and the initial distance $l^{0}$

δ = l^{t} - l^{0}

(22)

Next, the local internal forces are expressed in terms of the previously defined deformations

{\bar{f}}^{e} = {\bar{K}}^{e} {\bar{u}}^{e}

(23)

where ${\bar{f}}^{e} = [N M_{x, a} M_{y, a} M_{z, a} M_{x, b} M_{y, b} M_{z, b}]^{T}$ is the vector collecting all local element forces, ${\bar{u}}^{e} = [δ θ_{x, a} θ_{y, a} θ_{z, a} θ_{x, b} θ_{y, b} θ_{z, b}]^{T}$ is the vector collecting all local deformations, and ${\bar{K}}^{e}$ is the local element stiffness matrix

{\bar{K}}^{e} = \frac{1}{l^{0}} [\begin{matrix} E A & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & G J & 0 & 0 & - G J & 0 & 0 \\ 0 & 0 & 4 E I_{y} & 0 & 0 & 2 E I_{y} & 0 \\ 0 & 0 & 0 & 4 E I_{z} & 0 & 0 & 2 E I_{z} \\ 0 & - G J & 0 & 0 & G J & 0 & 0 \\ 0 & 0 & 2 E I_{y} & 0 & 0 & 4 E I_{y} & 0 \\ 0 & 0 & 0 & 2 E I_{z} & 0 & 0 & 4 E I_{z} \end{matrix}]

(24)

The global forces f^e are related to the local forces ${\bar{f}}^{e}$ as follows

f^{e} = T {\bar{f}}^{e}

(25)

where $f^{e} = [f_{x, a} f_{y, a} f_{z, a} m_{x, a} m_{y, a} m_{z, a} f_{x, b} f_{y, b} f_{z, b}$ $m_{x, b} m_{y, b} m_{z, b}]^{T}$ is the vector containing forces and moments in the global x-, y-, and z-axes for beam end “a” and beam end “b”, and T is a transformation matrix defined as

T = {(\frac{\partial {\bar{u}}^{e}}{\partial u^{e}})}^{T}

(26)

where $u^{e} = [x_{a} y_{a} z_{a} α_{a} β_{a} γ_{a} x_{b} y_{b} z_{b} α_{b} β_{b} γ_{b}]^{T}$ is a vector collecting translations and rotations of beam end “a” and beam end “b” in the global x-, y-, and z-directions. For the actual calculation of the transformation matrix T, the reader is referred to the original paper.³³ The full internal force vector f is assembled from the element forces f^e.

Finally, the element tangent stiffness matrix is obtained from differentiation of the internal forces f^e

\begin{matrix} K^{e} = \frac{\partial f^{e}}{\partial u^{e}} = \frac{\partial f^{e}}{\partial {\bar{f}}^{e}} \frac{\partial {\bar{f}}^{e}}{\partial u^{e}} + \frac{\partial f^{e}}{\partial T} \frac{\partial T}{\partial u^{e}} \\ = T {\bar{K}}^{e} T^{T} + K_{g}^{e} = K_{e}^{e} + K_{g}^{e} \end{matrix}

(27)

The resulting stiffness matrix consists of an elastic part $K_{e}^{e}$ and a geometric part $K_{g}^{e}$ . Often, the stiffness matrix can be effectively approximated by neglecting the geometric stiffness, as will be shown in the examples section. This simplifies equation (27) to

K^{e} = T {\bar{K}}^{e} T^{T}

(28)

For the exact tangent stiffness matrix, the calculation of the geometric stiffness matrix is required, for which the procedure can be found in the original paper.³³ The full tangent stiffness matrix K is assembled from the element stiffness matrices K^e following a standard finite element procedure.

Examples

In this section, four cases are discussed to illustrate the performance of the proposed method for form-finding and analysis of bending-active structures. For all cases described later, the convergence of three solution techniques is compared. (1) The first one corresponds to the implicit dynamic relaxation method described earlier, for which the fictitious mass matrix is calculated using only the elastic part of the stiffness matrix (equation (28)) to reduce the computation time per iteration. To make sure this simplified implicit dynamic relaxation method remains numerically stable for all examples considered, we have chosen an arbitrary value of 1 rad/s for the target eigenfrequency $ω^{*}$ , 1 s for the time step $Δ t$ , and 1 for the damping ratio. (2) This approach is compared to the classic Newton–Raphson method (equation (17)), which is equivalent to the implicit dynamic relaxation method combined with a time step and damping ratio given by equation (16)). Moreover, the full tangential stiffness matrix is used (equation (27)). (3) The classic Newton–Raphson method is known to encounter convergence problems when nonlinear effects cause a significant change in stiffness. For cases where the classic Newton–Raphson method does not converge, the Newton–Raphson updating scheme (equation (17)) is stabilized by choosing the time step and/or damping ratio more conservative than the values given by equation (16). In the following, this method is referred to as stabilized implicit dynamic relaxation. (4) Finally, these solution procedures are compared to a dynamic relaxation method that employs a diagonal fictitious mass matrix, for which the diagonal elements are given by

M_{ii} = c \frac{Δ t^{2}}{2} \sum_{j = 1}^{n} | K_{ij} | i = 1, 2, . . ., n

(29)

where $M_{ii}$ is the ith element of the diagonal fictitious mass matrix, $| K_{ij} |$ is the absolute value of the element on row i and column j of the stiffness matrix, and c is a stability constant. This way to choose the fictitious masses is commonly regarded as the conventional approach.^28,35 Here, the stiffness matrix used in the conventional approach is also calculated using only the elastic part of the element stiffness matrix (equation (28)) to reduce the required computation time. Moreover, this approach is combined with kinetic damping instead of viscous damping, and a value of 1 s is chosen for the time step $Δ t$ . A few trial runs showed that $c = 1.1$ generally performed best. Therefore, c was set to 1.1 for all shown examples. All methods are implemented in MATLAB and run on a MacBook Pro (2015) with a 2.7 GHz Intel Core i5 processor. The convergence criterion is defined as: $∥ r ∥ \leq 10^{- 3} N$ , unless specified otherwise.

Initially curved cantilever beam

The first example is a frequently used benchmark case for nonlinear beam elements which involves the coupling of torsion, bending, and axial forces. This example is presented to validate our implementation of Crisfield’s³³ co-rotational beam element. A cantilevered arc beam is clamped at its base and subjected to a lateral load. In its undeformed state, the beam spans an angle of 45° in the x–z plane (Figure 3). The radius of curvature is 100 in (2.54 m). The beam has a square 1 × 1 in² (6.45 cm²) cross section and material properties: $E = 10^{7} psi$ (68.95 GPa), $ν = 0$ . The load is applied in the y-direction and has a magnitude of 600 lb (2.67 kN). Following Crisfield,³³ the beam is discretized into eight elements. The convergence criterion for this example was expressed in pounds: $∥ r ∥ \leq 10^{- 3} lb$ ( $4.45 \times 10^{- 3} N$ ).

Figure 3.

An initially curved cantilever beam (thin line) is deformed (thick line) by a vertical load applied at the beam’s tip (example 1).

Table 1 compares the results to those obtained by Crisfield.³³ Our results agree well with Crisfield’s³³ results (Table 1), with a relative difference of less than 0.3%.

Table 1.

Tip coordinates of the deformed beam obtained with our own implementation of the co-rotational beam element and tip coordinates reported by Crisfield³³ (example 1).

	Tip coordinates
Beam model	x (m)	y (m)	z (m)
Own implementation	46.89	53.68	15.57
Crisfield³³	46.84	53.71	15.61

Table 2 compares the number of iterations and computation time which are required for the loaded cantilever beam to reach equilibrium. For this case, the classic Newton–Raphson method was incapable of finding a solution directly. Consequently, the load must be applied in several steps, and equilibrium has to be found successively for each load step. Alternatively, a different choice for the time step and/or damping ratio can be a more intuitive way to stabilize the calculation. This technique is applied in the stabilized implicit dynamic relaxation approach; observations have shown that lowering the time step from $\sqrt{2} / ω^{*}$ to $1.3 / ω^{*}$ makes the method sufficiently stable to reach convergence in one load increment. The resulting stabilized implicit dynamic relaxation approach required the lowest number of iterations. However, the calculation of the geometric part of the stiffness matrix (see equation (27)) takes up a significant amount of computation time. As a consequence, the computation time for the simplified implicit dynamic relaxation approach was lower, even though it needed more iterations to converge. For reference, the calculation and inversion of the simplified stiffness matrix required approximately 8% of the total calculation time for the simplified implicit dynamic relaxation method, whereas the calculation and inversion of the exact tangential stiffness matrix needed in the stabilized implicit dynamic relaxation method required approximately 50% of the total calculation time. The conventional dynamic relaxation method requires significantly more iterations, making it much slower compared to the other two methods, even though every iteration is computationally cheaper because of the diagonal mass matrix.

Table 2.

Required number of iterations and computation time for example 1.

Method	K ^e	$Δ t$	$ξ^{*}$	Iterations	Time (s)
Simplified implicit DR	$K_{e}^{e}$	$\frac{1}{ω^{*}}$	1	40	0.1
Stabilized implicit DR	$K_{e}^{e} + K_{g}^{e}$	$\frac{1.3}{ω^{*}}$	$\frac{1}{\sqrt{2}}$	26	0.3
Newton–Raphson method	$K_{e}^{e} + K_{g}^{e}$	$\frac{\sqrt{2}}{ω^{*}}$	$\frac{1}{\sqrt{2}}$	–	–
Conventional DR	$K_{e}^{e}$	1	0	7435	17.7

DR: dynamic relaxation.

Results are shown for the simplified implicit DR approach, the Newton–Raphson method, stabilized implicit DR, and conventional DR. For each solution procedure, the components of the adopted element stiffness matrix, and the choice for the time step and damping ratio are given.

Slender beam buckling

The second example shows the accuracy of the proposed approach in a two-dimensional case. A slender rod of 10 m is simply supported and axially loaded by a compression load to cause buckling (Figure 4). The rod has a circular cross section with a diameter of 3 cm and material properties: $E = 30 GPa$ , $ν = 0$ . Four compression loads p are considered: 1.015p_c, 1.152p_c, 1.518p_c, and 2.541p_c, where $p_{c}$ is the critical buckling load, given by

p_{c} = \frac{π^{2} E I_{z}}{l^{2}}

(30)

which is equal to 117.7 N in this case. To compare the influence of the discretization on the accuracy of the results, the rod is discretized into 8, 16, 32, and 64 elements. In order to initiate the buckling, an initial displacement was introduced, which corresponds to the deformation caused by a lateral point load of magnitude $p_{c}$ acting in the middle of the rod. It was observed that for smaller initial displacements, the classic Newton–Raphson method converged to an unbuckled solution for some cases; the rod straightened out. Closer investigation of this phenomenon showed that the buckling of the beam causes an overestimation of the beam’s deformation in the first Newton–Raphson iteration, which leads to oscillations in the internal forces, especially in the normal forces. When the initial deformation was taken too low, these destabilizing oscillations led to the convergence to an unbuckled beam. This phenomenon was not observed with the other solving strategies, due to their more conservative convergence behavior.

Figure 4.

A slender rod is axially loaded to a post-buckled state (example 2).

The results are compared to analytical results given by the elastic theory.³⁶ The shape of the buckled rod is characterized by the position of the middle of the rod and the orientation of the rod ends. Table 3 compares the resulting buckling shapes for different discretizations to the analytical results. The numerical results agree well with the analytical results, especially for finer rod discretizations. Part of the difference between the numerical result and the analytical result is caused by the axial deformation of the rod, which is neglected in the analytical results. This is certainly the case for the first load (1.015p_c), as the influence of the axial deformation is relatively large for this load. Increasing the length of the rods, and consequently making them more slender, reduced the error, which confirms the assumption that the axial deformations causes part of the inaccuracy.

Table 3.

The rotation angle $α$ of the rod end, and the x and y coordinates of the middle of the rod are compared to the analytical result for four considered loads and four rod discretizations (example 2).

Load p	$α (^{°})$	$x (m)$	$y (m)$
$p = 1.015 p_{c}$
Analytical	20	4.85	1.10
8 Elements	7	4.98	0.41
16 Elements	18	4.88	0.96
32 Elements	19	4.86	1.06
64 Elements	20	4.85	1.08
$p = 1.152 p_{c}$
Analytical	60	3.71	2.97
8 Elements	58	3.81	2.88
16 Elements	59	3.73	2.95
32 Elements	60	3.71	2.96
64 Elements	60	3.71	2.97
$p = 1.518 p_{c}$
Analytical	100	1.75	3.96
8 Elements	99	1.83	3.98
16 Elements	100	1.77	3.96
32 Elements	100	1.75	3.96
64 Elements	100	1.75	3.96
$p = 2.541 p_{c}$
Analytical	140	–0.54	3.75
8 Elements	140	–0.47	3.82
16 Elements	140	–0.52	3.77
32 Elements	140	–0.53	3.76
64 Elements	140	–0.53	3.75

To compare the convergence behavior, the rod was discretized into 32 elements. Again for this case, the classic Newton–Raphson method was incapable of finding a solution directly (Table 4). Moreover, subdivision into smaller load steps would not be straightforward for this example, as the initial displacement would simply disappear if the first load increment is too small to initiate the buckling. Here, the method was stabilized by slightly reducing the time step to $1.41 / ω^{*}$ , which is shown in Table 4 as the stabilized implicit dynamic relaxation approach. In comparison, the simplified implicit dynamic relaxation approach performed well for $p = 1.152 p_{c}$ and $p = 1.518 p_{c}$ . However, for a load of 1.015p_c, this approach overestimates the actual stiffness of the deformed rod, leading to relatively slow convergence. On the other hand, for a load of 2.541p_c, neglecting the geometric stiffness results in an underestimation of the stiffness, leading to divergence. For the intermediate loads, the influence of the geometric stiffness is lower. Although the conventional dynamic relaxation approach has no difficulty in handling the large displacements, the method required much more iterations compared to the proposed implicit dynamic relaxation approach. Consequently, the required computation time is also much higher.

Table 4.

Required number of iterations and computation time for example 2.

Method	K ^e	$Δ t$	$ξ^{*}$	Iterations	Time (s)
$p = 1.015 p_{c}$
Simplified implicit DR	$K_{e}^{e}$	$1 / ω^{*}$	1	395	2.8
Stabilized implicit DR	$K_{e}^{e} + K_{g}^{e}$	$1.41 / ω^{*}$	$1 / \sqrt{2}$	37	1.2
Newton–Raphson method	$K_{e}^{e} + K_{g}^{e}$	$\sqrt{2} / ω^{*}$	$1 / \sqrt{2}$	12	0.7
Conventional DR	$K_{e}^{e}$	1	0	104,141	860.6
$p = 1.152 p_{c}$
Simplified implicit DR	$K_{e}^{e}$	$1 / ω^{*}$	1	70	0.6
Stabilized implicit DR	$K_{e}^{e} + K_{g}^{e}$	$1.41 / ω^{*}$	$1 / \sqrt{2}$	34	1.1
Newton–Raphson method	$K_{e}^{e} + K_{g}^{e}$	$\sqrt{2} / ω^{*}$	$1 / \sqrt{2}$	–	–
Conventional DR	$K_{e}^{e}$	1	0	104,972	882.1
$p = 1.518 p_{c}$
Simplified implicit DR	$K_{e}^{e}$	$1 / ω^{*}$	1	38	0.3
Stabilized implicit DR	$K_{e}^{e} + K_{g}^{e}$	$1.41 / ω^{*}$	$1 / \sqrt{2}$	41	1.3
Newton–Raphson method	$K_{e}^{e} + K_{g}^{e}$	$\sqrt{2} / ω^{*}$	$1 / \sqrt{2}$	–	–
Conventional DR	$K_{e}^{e}$	1	0	70,312	585.0
$p = 2.541 p_{c}$
Simplified implicit DR	$K_{e}^{e}$	$1 / ω^{*}$	1	–	–
Stabilized implicit DR	$K_{e}^{e} + K_{g}^{e}$	$1.41 / ω^{*}$	$1 / \sqrt{2}$	39	1.2
Newton–Raphson method	$K_{e}^{e} + K_{g}^{e}$	$\sqrt{2} / ω^{*}$	$1 / \sqrt{2}$	–	–
Conventional DR	$K_{e}^{e}$	1	0	61,696	526.4

DR: dynamic relaxation.

Strained gridshell

The third example was chosen to demonstrate the performance of the proposed method in a case involving a higher number of DoF. A strained gridshell is bent into shape by applying vertical point loads to every node of the initially flat grid (Figure 5). The deformation of the gridshell is calculated for four grids composed of 11 × 11, 21 × 21 (Figure 5), 31 × 31, and 41 × 41 initially straight rods, loaded, respectively, with point loads of 50, 25, $16.67$ , and $12.5 N$ . The undeformed grid covers 10 × 10 m². The rods have a diameter of 3 cm and material properties: $E = 30 GPa$ , $ν = 0$ . A rod segment between two adjacent intersections is modeled with one beam element. For simplicity, all rods are rigidly connected at the nodes. The gridshell is supported by roller supports at the four corners, allowing horizontal movement of the supports during the erection.

Figure 5.

The erection of a strained gridshell (thick lines) is simulated by vertically loading each node of an initially flat horizontal grid of straight rods (thin lines) (example 3).

Without further adjustments, the implicit dynamic relaxation method, and consequently, the Newton–Raphson method would not converge as the chosen boundary conditions allows for rigid body motion, which causes the fictitious mass matrix to be singular. Two solutions to this problem are proposed. The first approach is to add supports to prevent rigid body motion of the structure: the central node is only allowed to move vertically, whereas the middle node of the outer rods can only move vertically and horizontally toward the middle. The second approach is to add fictitious mass to all nodes, consequently increasing the diagonal elements of the fictitious mass matrix. For this example, an identity matrix is added to the fictitious mass matrix. This corresponds to adding 1 N s²/m of lumped mass to each DoF. For reference, the mean value of the diagonal elements of the original fictitious mass matrix is of order of magnitude $10^{7}$ N s²/m. Therefore, by adding the identity matrix, the diagonal elements of the fictitious mass matrix are only slightly increased. The latter approach results in the separation of the eigenfrequencies. Consequently, if the incremental mass is chosen too large, a reduction of the convergence rate is observed. On the other hand, if the incremental mass is chosen too small, the effectiveness is low and the mass matrix remains ill-conditioned.

Table 5 gives the results for the added support approach, the results for the added mass approach are virtually identical. For this example, the classical Newton–Raphson method did not converge, even when rigid body motion was eliminated using additional supports. A more conservative choice of $1.1 / ω^{*}$ for the time step, and $0.9$ for the damping ratio was chosen to stabilize the method. For the implicit dynamic relaxation approach, the number of iterations is near to insensitive to the number of DoF (Table 5). This is an important difference with the conventional dynamic relaxation method, where the difference between the lowest and the highest eigenfrequency increases for a finer mesh discretization, causing the convergence rate to decrease.

Table 5.

Required number of iterations and computation time for example 3.

Method	K ^e	$Δ t$	$ξ^{*}$	Iterations	Time (s)
Grid size: 11 × 11 (718 DoF)
Simplified implicit DR	$K_{e}^{e}$	$1 / ω^{*}$	1	39	2
Stabilized implicit DR	$K_{e}^{e} + K_{g}^{e}$	$1.1 / ω^{*}$	0.9	48	11
Newton–Raphson method	$K_{e}^{e} + K_{g}^{e}$	$\sqrt{2} / ω^{*}$	$1 / \sqrt{2}$	–	–
Conventional DR	$K_{e}^{e}$	1	0	9544	560
Grid size: 21 × 21 (2638 DoF) (Figure 5)
Simplified implicit DR	$K_{e}^{e}$	$1 / ω^{*}$	1	40	9
Stabilized implicit DR	$K_{e}^{e} + K_{g}^{e}$	$1.1 / ω^{*}$	0.9	50	44
Newton–Raphson method	$K_{e}^{e} + K_{g}^{e}$	$\sqrt{2} / ω^{*}$	$1 / \sqrt{2}$	–	–
Conventional DR	$K_{e}^{e}$	1	0	14,922	3297
Grid size: 31 × 31 (5758 DoF)
Simplified implicit DR	$K_{e}^{e}$	$1 / ω^{*}$	1	40	23
Stabilized implicit DR	$K_{e}^{e} + K_{g}^{e}$	$1.1 / ω^{*}$	0.9	53	105
Newton–Raphson method	$K_{e}^{e} + K_{g}^{e}$	$\sqrt{2} / ω^{*}$	$1 / \sqrt{2}$	–	–
Conventional DR	$K_{e}^{e}$	1	0	21,751	11,124
Grid size: 41 × 41 (10,078 DoF)
Simplified implicit DR	$K_{e}^{e}$	$1 / ω^{*}$	1	41	51
Stabilized implicit DR	$K_{e}^{e} + K_{g}^{e}$	$1.1 / ω^{*}$	0.9	51	183
Newton–Raphson method	$K_{e}^{e} + K_{g}^{e}$	$\sqrt{2} / ω^{*}$	$1 / \sqrt{2}$	–	–
Conventional DR	$K_{e}^{e}$	1	0	34,802	30,424

DR: dynamic relaxation.

Bending-active cable net structure

Finally, the fourth example shows a combination of bending rods and cable elements with a predefined force per unit length, also known as force density. The case was first described by Van Mele et al.¹⁷ Four 11.5-m-long rods are placed on the corners of a 14.14 × 14.14 m², with a cable net wrapped around them (Figure 6). The rods have a diameter of 2 cm and material properties: $E = 30 GPa$ , $ν = 0$ for two opposite rods, and $E = 90 GPa$ , $ν = 0$ for the two other rods. A force density of 1 N/m was assigned to the inner cable elements, and 10 N/m to the cable elements situated at the boundaries. The rods are discretized into nine beam elements. Although 6-DoF beam elements allow the introduction of simple clamping supports, the rods are fixed to the ground using two adjacent pin supports, similar to the original case. Two convergence criteria were considered; a force-based convergence criterion, stating that the norm of the residual forces must be lower than $10^{- 4} N$ , and a velocity based convergence criterion, stating that the norm of the velocity must be lower than $10^{- 5} m / s$ .

Figure 6.

A cable net (thin lines) is wrapped around four fixed initially straight rods (thick lines) to form a tent-like structure (example 4).

For this last example, the solution procedures discussed so far are compared to the dynamic relaxation approach proposed by Van Mele et al.¹⁷ Both the results obtained by Van Mele et al., and results from our own implementation of their method are used. Our own implementation is run in MATLAB, and on the same computer as the other methods. Van Mele’s approach does not require the calculation or inversion of a global stiffness matrix. Moreover, the method uses beam elements with only 3 DoF per node. Therefore, it is referred to as 3-DoF DR in Table 6.

Table 6.

Required number of iterations and computation time for the bending-active cable net structure (example 4).

Method	K ^e	$Δ t$	$ξ^{*}$	Iterations	Time (s)
Force convergence criterion
Simplified implicit DR	$K_{e}^{e}$	$1 / ω^{*}$	1	42	5
Newton–Raphson method	$K_{e}^{e} + K_{g}^{e}$	$\sqrt{2} / ω^{*}$	$1 / \sqrt{2}$	12	5
Conventional DR	$K_{e}^{e}$	1	0	69,723	7923
3-DoF DR (Van Mele et al.¹⁷)	NA	1	0	158,400	478
3-DoF DR (own implementation)	NA	1	0	58,863	66
Velocity convergence criterion
Simplified implicit DR	$K_{e}^{e}$	$1 / ω^{*}$	1	26	3
Newton–Raphson method	$K_{e}^{e} + K_{g}^{e}$	$\sqrt{2} / ω^{*}$	$1 / \sqrt{2}$	12	5
Conventional DR	$K_{e}^{e}$	1	0	38,744	4280
3-DoF DR (Van Mele et al.¹⁷)	NA	1	0	21,000	62
3-DoF DR (own implementation)	NA	1	0	24,939	27

DR: dynamic relaxation; DoF: degrees of freedom.

Results are shown for the simplified implicit DR approach, the Newton–Raphson method, conventional DR, and the DR approach proposed by Van Mele et al.¹⁷ For each solution procedure, the components of the adopted element stiffness matrix, and the choice for the time step and damping ratio are given.

The results show that the approach described by Van Mele et al. results in relatively inexpensive iterations: approximately 0.003 s per iteration, and 0.001 s per iteration for our own implementation (Table 6). This is mainly because this approach does not require the calculation and inversion of a stiffness matrix. Moreover, the use of 3-DoF beam elements reduces the total number of DoF in the model. Although the computation time per iteration is higher for the implicit dynamic relaxation method (approximately 0.118 s per iteration), the total computation time is much lower because of the strong reduction of iterations. The Newton–Raphson method requires the lowest number of iterations. However, with a computational cost of approximately 0.417 s per iteration, its total computation time is very similar to that of the simplified implicit dynamic relaxation approach.

Conclusion

This article proposes a fast and accurate dynamic relaxation approach for form-finding and analysis of bending-active structures. The implicit dynamic relaxation method as proposed by Rombouts et al.²⁹ was adopted. This method uses a stiffness proportional mass matrix to make all eigenfrequencies coincide, which results in an optimal trade-off between the amount of iterations and the computation time per iteration. This method was combined with co-rotational beam elements as described by Crisfield³³ to maintain high accuracy for large displacements, which are typical for bending-active structures.

Four examples involving active bending were shown. It is concluded that the proposed 6-DoF dynamic relaxation approach is accurate, even for large displacements and the coupling of torsion, bending moments, and axial forces. Moreover, a conservative choice for the time step and damping ratio makes the method more stable than the classical Newton–Raphson method, which means that even challenging problems involving strongly nonlinear behavior can be solved in a single load step. Furthermore, the method scales very well to problems with a large number of DoF. Due to the very low number of iterations, it is faster than the conventional dynamic relaxation method, even though the average calculation time per iteration is higher. Finally, a simplification of the stiffness matrix speeds up the calculation. The combination of high accuracy and low computation time makes this 6-DoF implicit dynamic relaxation method well suited for both form-finding and analysis of spatial structures undergoing large displacements and small strains, such as bending-active structures.

Footnotes

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was supported by the Research Foundation—Flanders (FWO; grant number G0C2315N).

ORCID iDs

Jef Rombouts

Lars De Laet

References

Rombouts

Lombaert

De Laet

, et al. On the equivalence of dynamic relaxation and the Newton-Raphson method: application to the design and analysis of bending-active structures. In: Proceedings of the IASS annual symposium (IASS 2017), Hamburg, 25–27 September 2017.

Lienhard

. Bending-active structures. PhD Thesis, Universität Stuttgart, Stuttgart, 2014.

Cuvilliers

Yang

Coar

, et al. A comparison of two algorithms for the simulation of bending-active structures. Int J Space Struct 2018; 33: 73–85.

Veenendaal

Block

. An overview and comparison of structural form finding methods for general networks. Int J Solids Struct 2012; 49: 3741–3753.

Day

. An introduction to dynamic relaxation. Engineer 1965; 219(5688): 218–221.

Barnes

. Form finding and analysis of tension space structures by dynamic relaxation. PhD Thesis, City University, London, 1977.

Topping

BHV

. The application of dynamic relaxation to the design of modular space structures. PhD Thesis, City University, London, 1978.

Papadrakakis

. Gradient and relaxation nonlinear techniques for the analysis of cable supported structures. PhD Thesis, City University, London, 1978.

Wakefield

. Dynamic relaxation analysis of pretensioned networks supported by compression arches. PhD Thesis, City University, London, 1980.

10.

Papadrakakis

. A method for the automatic evaluation of the dynamic relaxation parameters. Comput Method Appl M 1981; 25(1): 35–48.

11.

Barnes

. Form finding and analysis of tension structures by dynamic relaxation. Int J Space Struct 1999; 14(2): 89–104.

12.

Adriaenssens

. Stressed spline structures. PhD Thesis, University of Bath, Bath, 2000.

13.

Knippers

. Rotation formulations for dynamic relaxation with application in 3D framed structures with large displacements and rotations. In: International symposium of the IASS-APCS, Seoul, South Korea, May 21–24, 2012.

14.

Senatore

Piker

. Interactive real-time physics: An intuitive approach to form-finding and structural analysis for design and education. Comput Aided Design 2015; 61: 32–41.

15.

Bessini

Lázaro

Monleón

. A form-finding method based on the geometrically exact rod model for bending-active structures. Eng Struct 2017; 152: 549–558.

16.

Adriaenssens

Barnes

. Tensegrity spline beam and grid shell structures. Eng Struct 2001; 23: 29–36.

17.

Van Mele

De Laet

Veenendaal

, et al. Shaping tension structures with actively bent linear elements. Int J Space Struct 2013; 28(3–4): 127–135.

18.

Barnes

Adriaenssens

Krupka

. A novel torsion/bending element for dynamic relaxation modeling. Comput Struct 2013; 119: 60–67.

19.

Du Peloux

Tayeb

Lefevre

, et al. Formulation of a 4-DoF torsion/bending element for the formfinding of elastic gridshells. In: Proceedings of the international association for shell and spatial structures (IASS) symposium 2015: “future visions,” Amsterdam, 17–20 August 2015.

20.

D’Amico

Zhang

Kermani

. A finite-difference formulation of elastic rod for the design of actively bent structures. Eng Struct 2016; 117: 518–527.

21.

Lefevre

Tayeb

Peloux

, et al. A 4-degree-of-freedom kirchhoff beam model for the modeling of bending–torsion couplings in active-bending structures. Int J Space Struct 2017; 32(2): 69–83.

22.

Underwood

. Dynamic relaxation. In: Belytschko

Hughes

(eds) Computational methods for transient analysis. Amsterdam: New Holland, 1983, pp. 245–265.

23.

Zhang

. Modified adaptive dynamic relaxation method and its application to elastic-plastic bending and wrinkling of circular plates. Comput Struct 1989; 33(2): 609–614.

24.

Zhang

Kadkhodayan

Mai

. Development of the maDR method. Comput Struct 1994; 52(1): 1–8.

25.

Kadkhodayan

Alamatian

Turvey

. A new fictitious time for the dynamic relaxation (DXDR) method. Int J Numer Meth Eng 2007; 74(6): 996–1018.

26.

Rezaiee-Pajand

Sarafrazi

. Nonlinear structural analysis using dynamic relaxation method with improved convergence rate. Int J Comput Methods 2010; 7: 627–654.

27.

Rezaiee-Pajand

Kadkhodayan

Alamatian

, et al. A new method of fictitious viscous damping determination for the dynamic relaxation method. Comput Struct 2011; 89: 783–794.

28.

Alamatian

. A new formulation for fictitious mass of the dynamic relaxation method with kinetic damping. Comput Struct 2012; 90–91: 42–54.

29.

Rombouts

Lombaert

De Laet

, et al. On the equivalence of dynamic relaxation and the Newton-Raphson method. Int J Numer Meth Eng 2018; 113(9): 1531–1539.

30.

Belytschko

Schwer

Klein

. Large displacement, transient analysis of space frames. Int J Numer Meth Eng 1977; 11(1): 65–84.

31.

Bathe

. Finite element procedures. 2nd ed. Watertown, MA: Klaus-Jurgen Bathe, 2014.

32.

Felippa

Haugen

. A unified formulation of small-strain corotational finite elements: I. Theory. Comput Method Appl M 2005; 194: 2285–2335.

33.

Crisfield

. A consistent co-rotational formulation for non-linear, three-dimensional, beam-elements. Comput Method Appl M 1990; 81: 131–150.

34.

Crisfield

. Non-linear finite element analysis of solids and structures, vol. 2. Chichester: Wiley, 1997.

35.

Rezaiee-Pajand

Alamatian

. The dynamic relaxation method using new formulation for fictitious mass and damping. Struct Eng Mech 2010; 34(1): 109–133.

36.

Timoshenko

Gere

. Theory of elastic stability. Mineola, NY: Dover Publications, 2009.