CC BY
19
UDC 531.01:51-72: 624.074.5.06 Expansion Method for Continuating at Extended Solution at
Singular Points
V. V. Galishnikova
Department of Structural Engineering People's Friendship University of Russia 6 Mikluho-Maklaya str., Moscow 117198, Russian Federation
Stability of structures depends on the variation of the load factor on the continuation of their load path beyond a critical point. The computation of the continuation is known to be difficult due to special properties of the stiffness matrices of structures in the vicinity of singular points. A new expansion method is presented that leads to a robust and accurate continuation algorithm. The method is applied to the stability analysis of a spherical dome for symmetric and asymmetric snow loads.
Key words and phrases: geometrically nonlinear analysis, stability, singular points, continuation method.
1. Singular Points
Consider a structural model consisting of nodes, elastic finite elements, loads applied at the nodes and displacements prescribed at the nodes. The free displacements of the nodes in the direction of the applied loads and the reactions at the supports in the direction of the prescribed displacements are to be determined. The behaviour of the structure is described with a node displacement vector d and a node force vector f. The elements of d and f are ordered so that the free displacements precede the prescribed displacements. Vectors d and f are related by the total stiffness matrix K of the structure, which is decomposed into submatricesAiTOthat are compatible with the subdivision of d and f:
Kd
or
/ K11 K12 W di \ / fi \ ^ K21 K22 J - ^ f2 )
(1)
d1 — free displacements, f1 — applied loads, d2 — prescribed displacements, f2 — support reactions.
Behavior of the structure under load cannot be presented comprehensively with diagrams showing displacement coordinates as functions of load coordinates. The load path which shows the variation of the displacement norm with the force norm is more suited for this purpose. Each point on the load path in figure 1 corresponds to an equilibrium configration of the structure.
secondary branch C
primary branch B
' secondary branch D 'primary branch A
d
Figure 1. Load path of a structure
f
Received 29th December, 2010.
Consider a structure that is subjected to a load pattern and a displacement pattern so that the loading can be varied with a single parameter, which is called the load factor A:
fi = Afip, (2)
d2 = Ad2P, (3)
f1p — pattern load, d2p — prescribed displacement pattern.
The load path consists of a sequence of equilibrium configurations of the structure. As the load factor is increased from zero at the origin of the load path, the primary branch A in figure 1 is traversed. The displacement increment A d is related to the force increment A f by the tangent stiffness matrix Tk :
/ Tii T12 W A di \ f Afi \
A d =A f or (T2i TiJ (a dj ^AfJ (4)
Depending on the properties of the structure and its deformation under load, a critical load factor Ac with associated displacement dc and force fc may be reached. For this loading, an infinitesimal displacement increment 5d± is associated with a load factor increment A A = 0 so that 5f\ = 6d2 = 0.Denote the tangent stiffness matrix for displacement state dc by Tc. The following equation shows that determinant of submatrix T\ic must be null for Sdi = 0 :
/ Tiic Ti2C \ / ¿di \ f 0 \
T^d = f or (T2ic T22J (o ) = (fj (5)
The point on the load path corresponding to the critical load factor is called a singular point. The reliable detection and computation of singular points for complex structures is treated in paper [1].
2. Bifurcation Theory
The properties of load paths at singular points are studied by means of bifurcation theory [2]. Let the free node displacements be a function u(z) of a path parameter z, and the load factor a function A(z). Denote submatrix T\i by T and pattern load /ip by p. Equation (4) yields:
t £+„ ^=0. (6)
Matrix T is denoted by Tc at the singular point. It has at least one eigenvector x associated with eigenvalue null. If it has exactly one such eigenvector, the singular point is called simple.
Tc x = 0. (7)
Governing equation (6) is valid also at the singular point:
dwc dAc , ,
di" + ^^ = 0. (8)
If it is assumed that Tc is symmetric, multiplication of (8) from the left with xT and substitution of (7) leads to the condition for singularity of the stiffness matrix:
XTP d^ = 0. (9)
Equation (9) shows that there are two types of simple singular points which are illustrated in figure (1):
turning point Si : xTp =0 ^ = 0, (10)
bifurcation point S2 : xTp = 0. (11)
In order to determine a displacement increment on a load path beyond a turning point, the derivative of u is chosen proportional to the eigenvector x using a parameter di. Equation (12) satisfies the governing equation (8) due to result (10):
£ = »i«. (12)
In order to determine a displacement increment on a load path beyond a bifurcation point, the pattern load p is decomposed into components parallel and normal to the eigenvector x:
p = cx + q with xT q = 0 and xT x = 1, (13)
q = (I — xxT )p. (14)
The governing equation for the displacements due to q is obtained from (8) by multiplying from the left with I — xxT and using expressions (7) and (14):
(I — xxT )TC ^ = —(I — xxT )p ,
duc dAc
~dz = —^~dz . (15)
The derivative of u is chosen as a linear combination of the eigenvector x and an unknown vector w:
df = Oix + e2w, (16)
02TcW = — q . (17)
Since parameter d2 can be chosen freely, it is set equal to the derivative of A. The resulting equation cannot be solved for w because the tangent stiffness matrix is singular.
Tcw = —q. (18)
A classification of singular points by means of the bifurcation equation is given by Galishnikova [2, p. 227], Wagner [3, p. 110] and Crisfield [4, p. 346]. The classification is helpful in understanding structural behavior in the vicinity of singular points. In numerical analysis the classification is difficult to apply because higher order derivatives of the displacements are required and many cases must be distinguished.
The continuation method which is presented in this paper is applied without prior classification of the singular points. The classification topic is therefore not pursued further.
3. Continuation of the Load Path
A segment of the load path that is reached after a singular point has been passed is called a continuation of the load path. Equation (12) shows that there is a single continuation branch at a turning point such as point Si in figure 1. At a bifurcation point such as S2 in figure 1, there are three continuation branches: primary branch B as well as secondary branches C and D. If more than one eigenstate at a singular
point has eigenvalue null, linear combinations of the eigenvectors lead to an infinite set of continuations.
A structure does not necessarily become unstable when a singular point is reached on its load path. Stability depends on the variation of the load factor on the branch that is chosen for continuation. If the load factor increases on the branch, the structure is called stable for displacement states on this branch. If the load factor decreases on at least one branch the structure is called unstable at the singular point. In order to be able to decide whether or not a singular point is a point of instability of a structure it is necessary to continue the load path beyond the singular point.
The bifurcation theory in section 2 considers infinitesimal displacement and force increments. Numerical computations of the load path deal with finite displacement and force increments that take the structure from the equilibrium configuration at the singular point S to a neighboring equilibrium configuration C. The step from S to C will in general be associated with a load factor increment because the force increment is not null. This load factor increment is used to determine whether or not the structure is unstable on the branch.
Continuation is treated in monographs by Crisfield [4], Felippa [5], Galishnikova [2], Wagner [3] and Wriggers [6]. All methods of solution consist of branch switching followed by traversal of the selected branch. Branch switching leads from the primary branch on which the singular point was reached to another branch on which the analysis is continued. The traversal then follows this path. Wagner [3] in chapter 10 presents an engineering approach to branch switching. The analysis is continued on the branch which, at the simple singular point, is tangent to the eigenvector. If multiple eigenstates with eigenvalue null exist, the branch is tangent to a linear combination of the eigenvectors. In remark 10.1, Wagner observes (translation):
".... Since the load factor of the first equilibrium configuration on the secondary path can be less than at the bifurcation point, . . . . a curve following process must be used to permit load reduction. Too small a step size can lead to return to the primary path. Too large a step size can lead to divergence. Numerical experience shows that in the general case the transition to the secondary branch succeeds after 1 to 3 trials
Crisfield [4] in chapter 21 differentiates between simple branch switching (section 21.2), branch switching with higher order derivatives (section 21.3) and direct computation of the singular points (section 21.6) which is also presented by Wriggers [6] in section 7.2.1. In simple branch switching, a predictor step A d is taken in the direction described above for Wagner's method and the load factor increment is set to null. The trial displacement is corrected with a sequence of displacement increments orthogonal to A d. A Taylor series expansion of the equilibrium equations and the Newton-Raphson method are used for branch switching with higher order derivatives. On page 354 in [4] Crisfield emphasizes "that these techniques are not used in most current nonlinear finite element programs and, indeed, that it has yet to be demonstrated that they are economically viable."
The well-known difficulties encountered in branch switching are related by Galishnikova et. al [2] to the stiffness properties of structures in the vicinity of singular points. It is shown in chapter 7 of [2] with exact solutions for the load path of 2-bar trusses that the load path on the vicinity of the singular point is determined by the third order derivatives in the Taylor series expansion.
The traversal of the selected branch with an arc length method also proves difficult. Crisfield warns on page 355 in [4]: "We describe some arc length methods .... motivated by the problems with conven-tional arc length methods that have been reported by a number of workers". For tangent stiffness methods with negative eigenvalues, as encountered on branches where the structure is unstable, Crisfield concludes in section 21.7.3: "This is clearly an area where future work is required." This opinion has recently been confirmed by Galishnikova [7] who has determined the exact tangent stiffness matrix for the secondary branch of the load path of a steep regular tripod
under vertical load:
TK =
( 2w\ 2W\W2
2W1W2 2w2
2wi(1 + w3) 2W2 (1 + W3)
det TK = 0,
(19)
y2wi(l + W3) 2w2(1 + w3) 2(1 + W3)2 - m2 J
w\, w2 — normalized horizontal displacement coordinates of the apex, w3 — normalized vertical displacement coordinate of the apex, m — geometric ratio: base radius / height of tripod.
This tangent matrix is singular at all points of all continuations of the load path, so that arc length methods based on tangent stiffness matrices cannot be used to traverse these branches.
4. Expansion Method
The author has applied the following research strategy to overcome some of the difficulties that are encountered in the continuation of load paths as described in section 3:
— The singular point is computed with high accuracy using only properties of the primary branch.
— The continuation on the secondary branch is independent of the classification of the singular point.
— The change in load factor in the first step following the singular point is taken into account.
The strategy has lead to a new method of analysis which is based on an expansion of the governing equations with additional equations. The expansion differs from that presented by Crisfield [4], Wagner [3] and Wriggers [6] by including the finite increment of the load factor in the equations, as shown below. The method is based on the following hypotheses:
(a) The direction of the free displacement increment in the first load step beyond a singular point is approximately parallel to the eigenvector in the eigenstate with eigenvalue null, or to a freely chosen linear combination of multiple eigenvectors with eigenvalue null.
(b) The direction the free displacement increment in a load step s succeeding the first step of the continuation is approximately equal to the displacement increment in the preceding step s — 1.
The branch is traversed stepwise. The displacement vector at the beginning of the first step in the continuation equals that of the singular point. The following sequence of operations is required in each cycle of iteration of each load step of the continuation.
First trial state: In the first step s = 0 of the continuation the trial displacement is determined with hypothesis (a) using a specified displacement scale factor ^ :
( di \ /die + 'QXi \
d = dc + Vx or (^J = ^ d2c ) (20)
d — trial displacement at the end of the first cycle of iteration.
In the general step s > 0 of the continuation the trial solution is determined with hypothesis (b):
d = d(s) + (d(s) — d(s-2)) (21) ¿0) — displacement state at the start of load step s.
Force vector: The structure would be in equilibrium in the trial state if it were subjected to the force vector determined with equation (1):
Kd = f or
/ K11 K12 \ / di / fl \
^ KK21 K22 J ) = ^ f2 )
(22)
K — total stiffness matrix for the trial displacement state d.
The stiffness matrix is not assembled explicitly. Instead, the stiffness matrices for the elements are computed for the trial state, multiplied with the element displacement vectors and summed to obtain the trial force vector. The trial load is not proportional to the pattern load. It is decomposed into a component parallel to the pattern load and a component normal to the pattern load:
/1 = rf1P + e1 with
/1^1
= 0 and r =
fiPh fîpfi p
(23)
Termination in load step s: The trial state is assumed to be sufficiently accurate to terminate the iterative correction in load step s of the continuation if the norm of the normal component of the trial load is less than a specified fraction of the norm of the parallel component of the trial load:
V^i < (24)
Termination of the expansion method: The computation of the continuation with the expansion method is continued until the tangent stiffness matrix is sufficiently well-conditioned to permit analysis with the constant arc length method used on the initial primary path of the load path. This condition is tested with the decomposition D L DT of the tangent stiffness matrix of the trial state, where D is a diagonal matrix and L is a left triangular matrix with diagonal coefficients 1. The expansion method is terminated if the smallest absolute value ldmin|of the coefficients du in D exceeds a specified fraction vm of the smallest value deof the diagonal coefficients of De in the decomposition LeDeLJ of the linear elastic stiffness matrix of the structure:
Idminl ^ vmde ^ terminate the expansion method. (25)
Enforcement of the pattern load: The undesired normal load component could be removed by subtracting the displacements and reactions corresponding to this load component from the trial state. It is a special feature of the new expansion method presented in this paper that, instead, the displace-ments and reactions due to load ei + Qfip and prescribed displacement dd2p are subtracted from the trial state. The free parameter dis determined so that the displacement increment A ddue to loadings + 0/ip and 6d2p is normal to the displacement increment d — d(s).
T A d = e + 0fp+ A r,
(d, - d(s^ T A d = 0,
( T! 1
T21
Tx 2 T22
-hp \ 0
\(d1 - d^)T (da - 4S0
(s)\T
0
/
A d1 A d2 d
e1
A r-2 0
(26)
(27)
The unknowns A d2, A r2andd in equations (26) and (27) are determined as follows. Matrix T22 is decomposed and substituted into the first row of (28):
LDLT A di + 6Ti2d2p = ei + dfip. (29)
Equation (29) is multiplied from the left with D-2 L-2. The result is expressed with auxiliary variables:
z = LT A di = a + 6b, (30)
a = D-iL-2 ei, (31)
b = D-iL-i (fip — Ti2d2p). (32)
Condition (27) is satisfied by substituting the displacement increment from (30) into (27):
(di — d("))T L-T z + 6(d2 — 4S))T d2P = 0. (33)
Expression (33) is rewritten with auxiliary variables g and c as follows:
gT z + 9c = 0, (34)
g = L-i(di — d(]s))c = fa — d2])T d2P. Substitution of z from (30) into (34) yields the value of parameter 0 :
6 = — . (35)
g1 b + c ^ '
Improved trial state: The displacement increment is computed with expression (30). The reaction increment follows from the second row of equation (27):
Adi =L-i(a + db), (36)
A T2 = Adi + 6T22d2p. (37)
The improved trial displacement state is obtained by subtracting the incremental solution (36) from the trial displacement state as follows:
di =di— Adi, (38)
d2 = d2 — dd2P. (39)
The iteration in the load step is repeated starting with the computation of the force vector with expression (20). From the second cycle of iteration in the first load step onwards, the load factor increment should be computed so that it equals the value in the first cycle in the first load step.
5. Example of a Continuation
The expansion method described in section 4 has been implemented on the Java platform. The stability of a variety of space trusses such as shallow and steep tripods, columns, cantilevers, masts with rectangular and triangular sections, arches and several domes, some of these with arches at their edges, have been analyzed for a wide variety of load patterns. In each study, the existence of a singular point on the load path was investigated. If it existed, the critical load factor, the load paths including the continuation on a secondary path (if it existed) and the diagrams of the variation of the displacement of selected nodes with the load factor were studied. All analyses proved to be stable.
As an example, the stability analysis of the single layer spherical dome in figure 2 with spherical radius 29.0 m, base radius 19.5 m and height 7.531 m is presented in this section. All bars of the truss have an area of 0.002398 m2 and a modulus of elasticity of 0.7 x 105MPa. The dome has rigid supports and is subjected to the self weight of the bars and a sheet cover as well snow loads according to the Russian building code. The intensity of the symmetric snow load at the apex in load combination C1 is 1.57. The reference intensity of the eccentric snow load in combination C2 is 1.80.
Figure 2. Spherical lattice dome: (a) Plan, (b) Elevation
The critical load factor for buckling under load combination C1 is 2.194. The buckling zone is centered at the apex and the buckled shape is symmetric. The load path consisting of a primary branch and its continuation on a secondary branch is shown in figure 3(a). The singular point is a bifurcation point. The load path almost reverses on itself at the singular point.
1200
1000
800
600
400
200
1200
z
M
1000
0 40 80 120 160 200 240 calibrated displacement norm (m/(kN/m)
800
600
400
200
0 40 80 120 160 200 240 calibrated displacement normm/(kN/m)
Figure 3. Load paths of the dome in figure 2: (a) Load combination C1; (b) Load
combination C2
The critical load factor for buckling under load combination C2 is 1.694. Buckling is local on the fourth ring from the apex, in the direction of the load eccentricity. The
load path is shown in figure 3 (b). The singular point is a bifurcation point. The displacement versus load factor diagrams for the two nodes with the largest displacements are shown in figures 4 and 5.
1.6
-0.10 -0.08 displacement, m
ft
a X f b \ c
/ / zj I /" \ \ N \ j
/ \ \
/ \ \
J / w 1
- — - — < f \ \ \
1.6
0.8
0.0
Figure 4. Variation of load factor with displacement at node 42 in the buckling zone for C2 (a) vertical displacement component (b) (c) horizontal displacement components
-Ö 1.6
•S
Ü
o 1.2
0.8 0.4 0.0
a b c A
- / /
t (7 N
//
^»L f
1 \
§ H.
-0.020 -0.015
displacement, m
-0.010
-0.005
0.0
1.6
1.2
0.8 0.4
0.0
Figure 5. Variation of load factor with displacement at node 74 in the buckling zone for C2 (a) vertical displacement component (b) (c) horizontal displacement components
O 1.2
1.2
0.8
0.4
0.4
-0.06
0.04
-0.02
0.0
6. Conclusions
A new expansion method for the continuation of load paths beyond singular points has been developed and tested with stability analyses for space trusses. The example presented in this paper, and a set of examples which cannot be presented in the available space, show that the method is robust and highly accurate. The method can be applied to detect local buckling and to compute load paths which nearly double back on themselves near a singular point. Additional research is required to prove the suitabilty of the method for space frames, plates, folded plates and shells as well as structures with nonelastic material properties.
References
1. Galishnikova V. Stability Analysis of Space Trusses // International Journal for Computational Civil and Structural Engineering. — 2009. — Vol. 5, issue 1&2. — Pp. 35-44.
2. Galishnikova V., Dunaiski P., PahlP. J. Geometrically Nonlinear Analysis of Plane Trusses an d Frames. — Sun Press, Stellenbosch, 2009.
3. Wagner W. Zur Behandlung von Stabilitatsproblemen der Elastostatik mit der Methode der Finiten Elemente. — Bericht Nr. F91/1, Bereich Mechanik der Uni-versitat Hannover, 1991.
4. Crisfield M. A. Nonlinear Finite Element Analysis of Soliuds and Structures // Advanced Topics. — John Wiley & Sons, 1997. — Vol. 2.
5. Felippa C. Nonlinear Finite Element Methods. Course notes at the University of Colorado. — 2007.
6. Wriggers P. Nichtlineare Finite-Element-Methoden. — Springer Verlag, 2001.
7. Galishnikova V. V. Analytical Solution of the Nonlinear Stability Problem of a 3-D Truss // Journal of Volgograd State University for Architecture and Civil Engineering. Natural Science. — 2007. — Vol. 6 (18). — Pp. 53-64.
УДК 531.01:51-72: 624.074.5.06
Метод расширения для вычисления продолжения решения
в сингулярных точках В. В. Галишникова
Кафедра строительных конструкций и сооружений
Российский университет дружбы народов ул. Миклухо-Маклая, д.6, Москва, 117198, Россия
В работе излагается новый метод вычисления продолжения решения в сингулярных точках, основанный на расширении системы разрешающих уранений и обладающий хорошей сходимостью и точностью. Приводится пример расчёта сетчатого купола на устойчивость равновесия по разработанной автором программе, реализующей предложенный метод.
Ключевые слова: геометрическая нелинейность, устойчивость, сингулярные точки, продолжение решения.
