Efficient finite element analysis using graph-theoretical force method; rectangular plane stress and plane strain Lagrange family elements

Applied Mathematics and Computation 266 (2015) 72–94

Contents lists available at ScienceDirect

Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc

Efficient finite element analysis using graph-theoretical force method; rectangular plane stress and plane strain Lagrange family elements A. Kaveh∗, M.S. Massoudi Centre of Excellence for Fundamental Studies in Structural Engineering, Iran University of Science and Technology, Narmak, Tehran-16, Iran

a r t i c l e

i n f o

Keywords: Finite elements force method Flexibility matrix Graph theory Plane stress Strain rectangular Lagrange elements

a b s t r a c t Formation of a suitable null basis for equilibrium matrix is the main part of finite elements analysis via force method. For an optimal analysis, the selected null basis matrices should be sparse and banded corresponding to sparse, banded and well-conditioned flexibility matrices. In this paper, an efficient method is developed for the formation of null bases of finite element models (FEMs) consisting of rectangular plane stress and plane strain Lagrange family elements, corresponding to highly sparse and banded flexibility matrices. This is achieved by associating special graphs with the FEM and selecting appropriate subgraphs and forming the self-equilibrating systems (SESs) on these subgraphs. The efficiency of the present method is illustrated through three examples. © 2015 Elsevier Inc. All rights reserved.

1. Introduction The force method of structural analysis in which the member forces are used as unknowns is appealing to engineers since the properties of members of a structure most often depend on the member forces rather than joint displacements. This method was used extensively until 1960. The advent of the digital computer and the amenability of the displacement method for computation attracted most researchers. As a result, the force method and some of the advantages it offers in non-linear analysis and optimization has been neglected. Five different approaches are adopted for the force method of structural analysis, classified as: 1. 2. 3. 4. 5.

Topological force methods. Graph theoretical methods. Algebraic force methods. Mixed algebraic-combinatorial force methods. Integrated force method.

Topological methods have been developed by Henderson [1], Maunder et al. [2] and Henderson and Maunder [3] for rigidjointed skeletal structures. Graph theoretical methods based on cycle bases of the graph models are due to Kaveh [4,5]. These methods are generalized to cover different types of skeletal structures such as rigid-jointed frames, pin-jointed planar trusses and ball-jointed space trusses in [6,7].

∗

Corresponding author. Tel.: +98 21 77240104; fax: +98 21 73223106. E-mail address: [email protected], [email protected] (A. Kaveh).

http://dx.doi.org/10.1016/j.amc.2015.05.064 0096-3003/© 2015 Elsevier Inc. All rights reserved.

A. Kaveh, M.S. Massoudi / Applied Mathematics and Computation 266 (2015) 72–94

73

Algebraic methods have been developed by Denke [8], Robinson [9], Topçu [10], Kaneko et al. [11], and Soyer and Topçu [12]. Mixed algebraic-topological methods have been used by Gilbert and Heath [13], Coleman and Pothen [14,15], Pothen [16], and Heath et al. [17]. The integrated force method has been developed by Patnaik et al. [18], Wei and Patnaik [19], in which the equilibrium and compatibility conditions are satisfied simultaneously in terms of the element force variables. Doiphode [20] utilized integrated force method for the analysis of frames and continuum structures. Other applications of the force method are due to Chi and Lee [21] who used the force methods for analysis of trusses with elastic boundary conditions. Konstantinos et al. [22] used a force-based approach for the limit analysis of three-dimensional frames. Nam-Il Kim et al. [23] employed the force method for damage identification of truss structures. Recently applications of the graph theory methods are extended to two classes of finite element models. The first class takes the element forces along the edges of the elements [21–25] and in the second class the element forces are concentrated at the mid-edge of the edges of the elements [26,27]. In this paper, an efficient method is developed for the formation of null bases for finite element models comprising of rectangular plane stress and plane strain Lagrange family elements leading to highly sparse and banded flexibility matrices, and can be used for optimal finite element analysis by the force method. This is achieved by associating a special graph to the finite element model and selecting subgraphs (known as γ -cycles [6]) for the formation of localized self-equilibrating stress systems (null vectors). Their numerical values are calculated by an algebraic process. The efficiency and accuracy of the present method is illustrated through simple examples. 2. Formulation of force method Consider a discrete or discretized structure which is statically indeterminate. The m-dimensional vector r contains independent element (member) forces, and an n-dimensional vector p denotes the nodal loads. The equilibrium equations of the structure can then be expressed as:

Ar = p

(1)

where A is an n×m equilibrium matrix. Assuming stability for the structure, the equilibrium matrix will have full rank, i.e. t = m − n > 0, rank(A ) = n. Also the member forces can be written as the sum of the particular and complementary solutions, where q is the tdimensional vector of the redundant forces.

r = B0 p + B1 q

(2)

B0 and B1 have m rows and n, and t columns, respectively. Pre-multiplying both sides of Eq. (2) by A and using Eq. (1) lead to

AB0 = I

(3)

AB1 = 0

(4)

Here, B0 and B1 are not unique for a structure and many of such matrices can be formed. B1 is called static basis or self-stress matrix. This basis is known as null basis in mathematics and each column of the null basis matrix is known as a null vector. The null space and null vectors are mathematical counterparts of the complementary solution space and self-equilibrating systems, respectively. Minimizing the complementary potential energy subjected to the constraint as in Eq. (1) requires r to minimize the quadratic form

minimize ( 12 rt F m r )

(5)

Here, Fm is a m×m block diagonal matrix known as the unassembled flexibility matrix containing the flexibility matrices of the elements of a structure in its block diagonal entries. Coupling Eqs. (5) and (2) results in

q = −(Bt1 F m B1 )−1 (Bt1 F m B0 ) p

(6)

According to Eq. (6) by solving a set of equations, redundant forces can be found. After the determination of the member forces, using the load–displacement relationship for each member, one can write member distortion as



[u] = [F m ][r] = [F m ] B0 B1

   p

(7)

q

Using virtual work, nodal displacements can be calculated as

 

t [v0 ] = B0 [u]

(8)

Combining Eqs. (7) and (8) leads to

v0 = Bt0 F m B0 p + Bt0 F m B1 q Substituting Eq. (6) in Eq. (9) and using Di j =



v0 = D00 −

D01 D−1 11 D10



p = Fp

(9) Bti F m B j

leads to (10)

74

A. Kaveh, M.S. Massoudi / Applied Mathematics and Computation 266 (2015) 72–94

Therefore the overall flexibility matrix of structure is obtained as

F = D00 − D01 D−1 11 D10

(11)

For free vibration of linear structure without damping we have





[K] − ω2 [M] [v0 ] = 0

(12)

Obviously K v0 = p and substituting Eq. (10) in Eq. (12) leads to





[I] − ω2 [M][F ] [p] = 0

(13)

Then the frequency equation of the system in the force method is obtained as

|[M][F ] − λ[I]| = 0 and λ = 1/ω2

(14) Bt1 F m B1

Efficiency of this analysis depends on the required time for the formation of the matrix G = and its characteristics, i.e. sparsity and bandedness together with its conditioning. For the formation of a well-structured matrix G, one should select a well-structured B1 matrix. Many algebraic procedures based on various matrix factorizations such as Gauss–Jordan elimination, LU, QR, LQ exist for the formation a null basis matrix B1 of an equilibrium matrix A. Basic concept of these methods is described briefly in the following. Let matrix A be partitioned using a column permutation matrix P as below:

AP = [A1 , A2 ]

(15)

Where A1 is a n×n non-singular matrix. Obviously matrix B1 can be written as



B1 = P

−A−1 1 A2 I



(16)

The complete description of algebraic force method can be found in Ref. [30] and for brevity is not repeated in here. 3. Independent element forces and flexibility matrix of elements For the generation of the equilibrium matrix A of an FEM, a set of independent forces system should be defined and also their relations with the element nodal forces should be established [28]. In displacement method we have two forces at each node of the element. For an element with N nodes, 2 × N nodal forces can be defined. Using three equilibrium equations, 2N − 3 independent forces will remain. In other words, there are 2N − 3 independent element forces in an element with N nodes. The nodal forces and element forces systems are shown in Fig. 1 for rectangular plane stress and plane strain Lagrange family elements with various numbers of boundary nodes. For the higher order elements, the element forces system can be obtained with the same procedure. These element forces can be related to the nodal forces for a rectangular element by a (2N ) × (2N − 3 ) transformation matrix using Eq. (17) as

S = TF

(17)

Transformation matrix can be formed simply as

(n1 , n2 ) = end nodes of element force Fj For i = 1 : N For j = 1 : 2N − 3 If i = n1 T (2i − 1, j ) = mn1 n2 and T (2i, j ) = nn1 n2 If i = n2 T (2i − 1, j ) = mn2 n1 and T (2i, j ) = nn2 n1 End End where xi and yi are the Cartesian coordinates of node i, mi j =

xi −x j li j

,

ni j =

yi −y j li j

are the direction cosines and li j is the length of

the line between nodes i and j. Formulation of a discrete element equivalent to the actual continuous structure is the first step in matrix structural analysis. For a linear system it can be assumed that the stresses σ are related to the forces F by linear equation as

¯ σ = cF

(18)

The matrix c¯ represents a statically equivalent stress system due to the unit force F. The flexibility matrix of an element can be written as



fm =

V

¯ c¯t ϕ cdV

(19)

The integration is taken over the area of the element, where ϕ is the matrix relating the stresses to strains ε = ϕσ in two dimensional problems [29]. The primary step in the formation of the flexibility matrix of an element is determining the matrix

A. Kaveh, M.S. Massoudi / Applied Mathematics and Computation 266 (2015) 72–94

Element

Nodes number

Size of T

Nodal forces

Element forces

Linear

N=4

8×5

S8×1 = [S1...S8 ]t

F5×1 = [ F1...F5 ]t

Quadratic

N=9

18× 15

S18×1 = [ S1...S18]t

F15×1 = [ F1...F15]t

Cubic

N = 16

32 × 29

S32×1 = [ S1...S32]t

F29×1 = [ F1...F29]t

Quartic

N = 25

50 × 47

S50×1 = [S1...S50]t

F47×1 = [ F1...F47]t

75

Element forces

Fig. 1. A set of rectangular Lagrange family elements.

c¯ . It is obvious that the ith column of c¯ represents the resultant stresses due to unit element force Fi in the force method and also stresses due to nodal forces S is equal to the ith column of T utilizing the displacement method. Hence, we can form matrix c¯ using the stiffness properties of the rectangular element using the displacement method. Now the flexibility matrix of the element in the force method is formed from Eq. (19) using Gauss numerical integration method with four Gauss points. 4. Graphs associated with finite element models 4.1. Basic graph theory definitions A graph S consists of a set of elements, N(S), called nodes and a set of elements, M(S), called members, together with a relation of incidence which associates two distinct nodes with each member, known as its ends. Two nodes of a graph are called adjacent if these nodes are the end nodes of a member. A member is considered incident with a node if it is an end node of the member. The degree of a node is the number of members incident with that node. A subgraph Si of a graph S is a graph for which N(Si ) ⊆ N(S) and M(Si ) ⊆ M(S), and each member of Si has the same ends as in S. A path of S is a finite sequence Pi = {n0 , m1 , n1 , ..., mp , np } whose terms are alternately distinct nodes ni and distinct members mi of S for 1 ≤ i ≤ p, and ni-1 and ni are the two ends of mi . Two nodes ni and nj are said to be connected in S if there exists a path between these nodes. A cycle is a path (n0 , m1 , n1 , …, mp , np ) for which n0 = np and p ≥ 3; i.e. a cycle is a closed path. The cycles of a graph form a vector space known as the cycle space. The dimension of this space for a connected graph S is known as the first Betti number, b1 (S) = M(S)−N(S)+1, of the graph, where M(S) and N(S) are the number of members and nodes of S, respectively. In order to transfer the topological property of a finite element model to the connectivity property of a graph model, ten different types of graphs are introduced in [30,31]. 4.2. An interface graph The interface graph of a finite element model denoted by IG (FEM) can easily be constructed for rectangular FEM using the following rules: 1. This graph contains all the nodes of the FEM.

76

A. Kaveh, M.S. Massoudi / Applied Mathematics and Computation 266 (2015) 72–94

Fig. 2. A quadratic, cubic and quartic rectangular FEM with their interface graphs.

Fig. 3. An FEM with its natural associate graph (blue lines). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 4. Subgraph corresponding to SESs of Type II for two adjacent cubic elements.

2. With the all edges of an element of FEM, some graph members are associated. Therefore, in the interface of two elements, 2-multiple members are presented. 3. For each element with N nodes, 2N−3 members should be considered in the interface graph. Thus, ((2N−3) −edge members) internal members should be added. This graph for a quadratic, cubic and quartic FEM is shown in Fig. 2. The member numbering of the interface graph should be performed according to the numbering of the FEM, taking into account the primary nodal numbering of a considered element in the model. Thus, for each rectangular element 2N−3 members of the interface graph will be numbered sequentially according to the patterns which were illustrated in Fig. 1. 4.3. Natural associate graph The natural associate graph represented by NAG(FEM) is constructed by the following rules: 1. Nodes of the NAG(FEM) correspond to the elements of FEM. 2. For each pair of elements in FEM having 2, 3, 4 and 5 common nodes for linear, quadratic, cubic and quartic elements respectively, one member is added between the corresponding two nodes in NAG(FEM). NAG can be constructed using the following procedure: One of the preliminary steps in FEM is defining the elements with their connected nodes. In this way the element connectivity matrix is constructed which contains the element-node incidence relationships. In the process of constructing the element connectivity matrix, another matrix which contains node-element incidence properties can be formed. This matrix is named the node connectivity matrix. Now using the element connectivity and the node connectivity matrices leads to an algorithm with complexity O(n) for an efficient generation of NAG.

A. Kaveh, M.S. Massoudi / Applied Mathematics and Computation 266 (2015) 72–94

77

Fig. 5. The SES of Type III corresponding to the common node of four linear elements. Table 1 Generators of Type II SESs in two directions. Element

Horizontal direction (i < j)

Vertical direction (i < j)

Quadratic Cubic Quartic

15 × (i−1)+11 29 × (i−1)+[15,17] 47 × (i−1)+[19,22,24]

15 × (i−1)+13 29 × (i−1)+[21,23] 47 × (i−1)+[31,36,38]

Table 2 Generators of Type III SESs. Element

Linear Quadratic Cubic Quartic

i < { j, m, n}

5 × (i−1)+5 15 × (i−1)+9 29 × (i−1)+25 47 × (i−1)+42

In order to recognize the adjacent elements to the nth element which have specified common nodes or one common face, first the connected nodes to the nth element are identified from the element connectivity matrix. In the subsequent step using the node connectivity matrix, elements which have at least one common node with the nth element are identified. Now it is convenient to seek for the adjacent elements in this reduced search space. An FEM and its corresponding NAG are illustrated in Fig. 3. 5. Pattern corresponding to the self-equilibrating systems Considering Fig. 1, in order to find the patterns corresponding to the self-equilibrating systems, a rectangular element is simulated as a planar truss formed as the 1-skeleton of the rectangular element together with some diagonal members. This is possible since the independent element forces applied at in the nodes and are along the edges of the rectangular element. The statical indeterminacy of a planar truss with m members and n nodes is given as γ (S ) = m − 2n + 3; therefore, the degree of statical indeterminacy (DSI) of the entire model supported in a statically determinate manner can be calculated with the same relationship as

DSI = (2N − 3 ) × M − 2n + 3

(20)

where M is the total number of elements, N is the number of nodes within an element and n is the total number of nodes of the FEM. With the above simulation, the patterns of the self-equilibrating systems can be identified as follows: 5.1. Type I self-equilibrating systems Each multiple member of the interface graph is a subgraph on which one self-equilibrating system can be generated. In other words, on a 2-multiple member numbered as (i, j) with the condition (i < j), one self-equilibrating system can be constructed (extracted).

78

A. Kaveh, M.S. Massoudi / Applied Mathematics and Computation 266 (2015) 72–94

Fig. 6. Finite element model of Example 1. Table 3 Topological properties for the FEMs of Example 1. Element

N

Linear Quadratic Cubic Quartic

169 169 169 169

M

DSI

SESs Type I

SESs Type II

SESs Type III

144 36 16 9

385 205 129 88

264 120 72 48

– 60 48 36

121 25 9 4

Each pair such as (i, j) for which (i < j) corresponds to a null vector with their non-zero entries being located in rows i and j, and their numeric values are (−1, 1), respectively. The member with bigger member number (j) is called the generator. These pairs are called Type I self-equilibrating systems.

 For an FEM we have

Linear N−0 =1 4 =2 Quadratic N−1 4 × M Cubic N−4 =3 4 =4 Quartic N−9 4

self-equilibrating systems of Type I, where M is the number of members of the

associate graph of the model. 5.2. Type II self-equilibrating systems There are other types of self-equilibrating systems which are extracted from two adjacent elements of FEM. In other words, for two adjacent elements with N nodes, the DSI can be calculated as:

DSI = (2N − 3 ) × M − 2n + 3

⎛ ⎜

⇒ DSI = (2N − 3 ) × 2 − 2 × ⎝2N −



Linear N−0 4  Quadratic N−1 4 N−4 Cubic 4 N−9 Quartic 4

systems is

Type II =

⎧ ⎪ ⎨ Linear

N−0 4 Quadratic N−1 4 Cubic N−4 4 Quartic N−9 4

⎪ ⎩

⎫⎞

⎧

N−0 + 1⎪ ⎪ 2 ⎬ ⎨ Linear N−1 +1 ⎟ Quadratic 2 + 3 = ⎠ N−4 + 1⎪ ⎪ 2 ⎭ ⎩ Cubic N−9 +1 Quartic 2

⎫

− 1⎪ ⎬ −1 − 1⎪ ⎭ −1

(21)

self-equilibrating systems were generated as Type I systems. Thus the number of remaining self-equilibrating

⎧ ⎪ ⎨ Linear

N−0 2 Quadratic N−1 2 N−4 Cubic 2 Quartic N−9 2

⎪ ⎩ 



DSI

⎫ ⎧

⎫

⎧

N−0 − 1⎪ ⎪ Linear N−0 4 ⎪ 4 ⎬ ⎨ ⎬ ⎪ ⎨ Linear N−1 −1 Quadratic N−1 Quadratic 4 4 − = N−4 N−4 − 1⎪ ⎪ Cubic 4 ⎪ 4 ⎭ ⎩ ⎭ ⎪ ⎩ Cubic N−9 −1 Quartic N−9 Quartic 4 4

 



Type I







⎫

− 1 = 0⎪ ⎬ −1=1 − 1 = 2⎪ ⎭ −1=3

(22)



Type II

In other words, these SESs should be extracted from two adjacent elements. For example, the remaining subgraphs for two adjacent cubic elements are shown in Fig. 4. After deleting the generators corresponding to Type I SESs, the null vectors should be calculated from the remaining subgraph. These null vectors can easily be generated on the corresponding sub-structure utilizing an algebraic method.

A. Kaveh, M.S. Massoudi / Applied Mathematics and Computation 266 (2015) 72–94

79

Fig. 7. (a) FEMs, (b) interface graph and (c) natural associated graph of Example 1.

In an FEM, the total number of Type II SESs can be calculated as:

TypeII = M ×

⎧ ⎫ ⎪ ⎨ Linear 0⎪ ⎬ Quadratic1

⎪ ⎩ Cubic 2⎪ ⎭ Quartic 3

where

M

is the number of members of the associate graph of the model.

(23)

80

A. Kaveh, M.S. Massoudi / Applied Mathematics and Computation 266 (2015) 72–94

(a)

(b)

(c)

(d)

Fig. 8. Equilibrium matrix for FEMs of Example 1, (a) linear, (b) quadratic, (c) cubic and (d) quartic.

The most important point in Type II self-equilibrating systems is to select an appropriate generator. In fact by eliminating these generators from graph S, the sub-structure of Type III SESs and the primary structure of the structure S must remain stable. 5.3. Type III self-equilibrating systems Sub-structures which are topologically identical to the minimal cycles of the natural associate graph of FEM contains some Type I, II and one Type III self-equilibrating systems. 5.3.1. Type I minimal cycles of NAG(S) These minimal cycles of the natural associate graph of the FEM pass through four elements which have one common node. Corresponding interface graph of these elements have n nodes and m edges for an FEM with N-node elements.

m = 4 × (2N − 3 )

n = 4N − 4 ×

⎧ ⎪ ⎨ Linear

N−0 4 Quadratic N−1 4 Cubic N−4 4 Quartic N−9 4

⎪ ⎩

(24)

⎫

⎧

⎫

+ 1⎪ ⎪ ⎬ ⎨ Linear 3N − 3⎪ ⎬ +1 Quadratic3N − 2 +1= + 1⎪ ⎪ ⎭ ⎩ Cubic 3N + 1⎪ ⎭ Quartic 3N + 6 +1

(25)

Subsequently, the DSI of the interface graph is

DSI = m − 2n + 3 ⇒ DSI = 4 × (2N − 3 ) − 2 ×

⎧ ⎫ ⎪ ⎨ Linear 3N − 3⎪ ⎬ Quadratic3N − 2

⎪ ⎩ Cubic 3N + 1⎪ ⎭ Quartic 3N + 6

+3=

⎧ ⎫ ⎪ ⎨ Linear 2N − 3 ⎪ ⎬ Quadratic 2N − 5

⎪ ⎩ Cubic 2N − 11⎪ ⎭ Quartic 2N − 21

(26)

A. Kaveh, M.S. Massoudi / Applied Mathematics and Computation 266 (2015) 72–94

81

(a)

(b) Fig. 9. Null basis matrix of Example 2 for linear, quadratic, cubic and quartic element, (a) present method, (b) LU method.

 The

Linear N−0 4  Quadratic N−1 4 ×4 Cubic N−4 4 Quartic N−9 4

 SESs are Type I and there are

Linear N−0 −1 4 −1 Quadratic N−1 4 ×4 Cubic N−4 −1 4 −1 Quartic N−9 4

SESs of Type II.

⎧ ⎫ ⎧ ⎫ N−0 ) × 4 − ( N−0 − 1 ) × 4 ⎪ ⎪ Linear 1⎪ ⎪ 4 4 ⎨ Linear (2N − 3 ) − ( N−1 ⎬ ⎨ ⎬ N−1 Quadratic1 Quadratic (2N − 5 ) − ( 4 ) × 4 − ( 4 − 1 ) × 4 DSI − (TypeI&II ) = = N−4 N−4 ) × 4 − ( 4 − 1) × 4⎪ ⎪ 4 ⎩ Cubic (2N − 11 ) − ( N−9 ⎭ ⎪ ⎩ Cubic 1⎪ ⎭ Quartic 1 − 1) × 4 Quartic (2N − 21 ) − ( 4 ) × 4 − ( N−9 4

(27)

Therefore, one independent SES should be extracted. This SES with 12 members can be formed for linear types of rectangular elements around the common node as is indicated bold in Fig. 5.

82

A. Kaveh, M.S. Massoudi / Applied Mathematics and Computation 266 (2015) 72–94

(a)

(b)

Fig. 10. Flexibility matrix of Example 2 for linear, quadratic, cubic and quartic element, (a) present method, (b) LU method.

A. Kaveh, M.S. Massoudi / Applied Mathematics and Computation 266 (2015) 72–94

83

Table 4 Comparison of the optimality characteristics of the null basis matrices B1 and flexibility matrices G for the FEMs of Example 1. Element

DSI

Linear Quadratic Cubic Quartic

385 205 129 88

nnz(B1 )LU method nnz(B1 )Present method

TimeLU method

DKI

Time

338 338 338 338

Present method

63.11 37.48 17.29 11.75

3.89 2.68 1.86 1.50

Condition number λλmax

Norms max |A × B1|

Present method

LU method

Present method

1.13 × 104 1.96 × 105 2.56 × 104 1.53 × 105

8.57 × 104 3.18 × 107 2.66 × 105 2.36 × 105

3.33 × 10−16 1.06 × 10−14 5.22 × 10−15 1.15 × 10−15

min

LU method 1.88 × 10−14 6.25 × 10−14 7.10 × 10−15 7.11 × 10−15

Fig. 11. (a) Finite element model of Example 1, (b) element numbering and (c) natural associated graph. Table 5 Topological properties for the FEMs of Example 2. Element

N

M

DSI

SESs Type I

SESs Type II

SESs Type III

Linear Quadratic Cubic Quartic

64 221 468 805

45 45 45 45

100 236 372 508

68 136 204 272

– 68 136 204

32 32 32 32

5.3.2. Type II minimal cycles of NAG(S) Each minimal cycle that surrounds an opening is called the Type II minimal cycle. Such a cycle passes through M ,(M ≥ 8 ), finite elements and its corresponding interface graph has N M × (2N − 3 ) members. The DSI of subgraph is

DSI = M × (2N − 3 ) −

⇒ DSI = M ×

⎧ ⎪ ⎨ Linear

⎫ ⎪ ⎬

⎪ ⎩

⎪ ⎭

3N−4 4 Quadratic 3N−3 4 2× 3N Cubic 4 Quartic 3N+5 4 Linear N−2 2 Quadratic N−3 2 +3 Cubic N−6 2 N−11 Quartic 2

⎧ ⎪ ⎨

⎫ ⎪ ⎬

⎪ ⎩

× M

−



Linear N−0 +1 4 +1 Quadratic N−1 4 × M +1 Cubic N−4 4 N−9 Quartic 4 +1

=



Linear 3N−4 4  Quadratic 3N−3 4 × M 3N Cubic 4 3N+5 Quartic 4

nodes and

× M + 3 (28)

⎪ ⎭

that some SESs of Type I and some SESs of Type II can be extracted.

⎧ ⎫ ⎧  N−2      ⎫ Linear × M + 3 − N−0 × M − N−0 − 1 × M ⎪ ⎪ Linear 3⎪ 2 4 4 ⎪ ⎪ ⎪  N−3      ⎨ ⎬ ⎨ ⎬ Quadratic 2 × M + 3 − N−1 × M − N−1 − 1 × M Quadratic3 4  4  N−6     DSI − (TypeI&II ) = = Cubic × M + 3 − N−4 × M − N−4 − 1 × M ⎪ ⎪ 2  4  4 ⎩ Cubic 3⎪ ⎭ ⎪ ⎪ ⎪  N−11  N−9  N−9  ⎩   ⎭ Quartic

2

×M +3−

4

×M −

4

−1 ×M

(29)

Quartic 3

Therefore, each Type II minimal cycle corresponds to three null vectors which are calculated utilizing an algebraic method.

84

A. Kaveh, M.S. Massoudi / Applied Mathematics and Computation 266 (2015) 72–94

(a)

(b)

(c)

(d)

Fig. 12. Equilibrium matrix for FEMs of Example 2, (a) Linear, (b) quadratic, (c) cubic and (d) quartic.

6. Selection of generators for SESs of Type II and Type III The most important point in Type II and Type III self-equilibrating systems is to select appropriate generators. This is by eliminating these generators from graph S, the sub-structure of primary structure of the structure S must remain stable. Tables 1 and 2 illustrate generators for Type II and Type III SESs. 7. Algorithm Step 1: Generate the associate graph of the FEM and use an efficient method for its nodal numbering [32,33]. It is obvious that good numbering of this graph corresponds to good numbering of elements of the FEM. This numbering leads to a banded adjacency matrix of the graph and correspondingly to a banded flexibility matrix. Since the numbering of the members of the interface graphs corresponds to the element numbering of the finite elements, such a numbering is the only parameter for controlling the bandwidth of the flexibility matrix. Step 2: Set up the equilibrium matrix of the FEM. Step 3: Generate the interface graph and perform its numbering. The numbering of this graph should be performed according to the element numbering of the considered FEM. After this numbering the interface graph can easily be formed and its members can be numbered. Step 4: Find the Type I self-equilibrating systems. All multiple members of the interface graph are identified and the values −1 and 1 are assigned to appropriate rows (corresponding to the member numbers) and the corresponding minimal null vectors are created. Step 5: Find the Type II self-equilibrating systems. The SESs of Type II should be extracted from two adjacent elements. Step 6: Find the Type III self-equilibrating systems. For each minimal cycle of natural associate graph of FEM with four members (one common node), one SES and with eight or more members (Opening), three SESs should be extracted. Step 7: Order the null vectors. At this step the constructed null vectors should be ordered such that their last non-zero entries form a list with an ascending order.

A. Kaveh, M.S. Massoudi / Applied Mathematics and Computation 266 (2015) 72–94

85

(a)

(b) Fig. 13. Null basis matrix of Example 2 for linear, quadratic, cubic and quartic element, (a) present method, (b) LU method.

Table 6 Comparison of the optimality characteristics of the null basis matrices B1 and flexibility matrices G for the FEMs of Example 2. Element

Linear Quadratic Cubic Quartic

DSI

100 236 372 508

DKI

128 442 936 1610

TimeLU method Time

Present method

14.23 47.72 67.95 98.04

nnz(B1 )LU method nnz(B1 )Present method

1.88 2.91 2.62 2.82

Condition number λλmax

Norms max |A × B1|

Present Method

LU Method

Present Method

5.88 × 102 1.29 × 103 2.57 × 103 6.39 × 103

1.08 × 103 1.24 × 104 2.23 × 104 8.43 × 104

2.22 × 10−16 2.22 × 10−16 4.21 × 10−15 6.66 × 10−16

min

LU Method 5.88 × 10−16 5.33 × 10−15 5.77 × 10−15 4.44 × 10−15

8. Numerical examples In this section three FEMs are considered which supported in a determinate fashion. The effect of the presence of additional supports can separately be included for each special case with no difficulty as discussed in [34,35]. The equilibrium matrices are formed. Null bases and the flexibility matrices are constructed and the required computational times, and the condition

86

A. Kaveh, M.S. Massoudi / Applied Mathematics and Computation 266 (2015) 72–94

(a)

(b)

Fig. 14. Flexibility matrix of Example 2 for linear, quadratic, cubic and quartic element, (a) present method, (b) LU method.

Element

Method

Stress (kN/m2 )

Element 17

Element 22

Element 23

Element 24

Element 29

Linear

Force method

σxx σyy σxy σxx σyy σxy σxx σyy σxy σxx σyy σxy σxx σyy σxy σxx σyy σxy

−41.6571461340121000 −53.2599313598811000 0.0000000000023000 −41.6571461340096000 −53.2599313598838000 0.0000000000000000 −19.8150309049441000 −52.2465128705436000 0.0000000000099000 −19.8150309049497000 −52.2465128705591000 −0.0000000000006000 −18.3501392729491000 −56.8446415698558000 −0.0000000000353000 −18.3501392729569000 −56.8446415697606000 0.0000000000005000

−106.3472777142010000 −92.2285227154218000 −79.9150829662687000 −106.3472777141990000 −92.2285227154259000 −79.9150829662709000 −97.5847257952226000 −52.7367712926246000 −75.5589536263580000 −97.5847257952208000 −52.7367712926275000 −75.5589536263611000 −93.7191526853095000 −46.6463027264143000 −63.9402623456037000 −93.7191526854048000 −46.6463027264167000 −63.9402623456028000

−43.9998969902881000 −51.7259163161676000 0.0000000000004000 −43.9998969902877000 −51.7259163161712000 −0.0000000000008000 −63.0131113936506000 −88.5240852475773000 0.0000000000021000 −63.0131113936457000 −88.5240852475836000 0.0000000000014000 −67.4751163685376000 −93.0602667786896000 −0.0000000000052000 −67.4751163686122000 −93.0602667786342000 0.0000000000020000

−106.3472777142020000 −92.2285227154238000 79.9150829662720000 −106.3472777142000000 −92.2285227154246000 79.9150829662713000 −97.5847257952279000 −52.7367712926234000 75.5589536263739000 −97.5847257952195000 −52.7367712926226000 75.5589536263616000 −93.7191526853221000 −46.6463027264209000 63.9402623455400000 −93.7191526853969000 −46.6463027264185000 63.9402623456045000

−71.1277868302311000 −213.0900972924210000 0.0000000000013000 −71.1277868302316000 −213.0900972924250000 −0.0000000000002000 −53.3198066094249000 −204.0757629652120000 0.0000000000059000 −53.3198066094234000 −204.0757629652200000 −0.0000000000018000 −45.0521766981594000 −198.9964201698490000 −0.0000000000289000 −45.0521766981803000 −198.9964201698400000 0.0000000000000000

Displacement method

Quadratic

Force method

Displacement method

Cubic

Force method

Displacement method

A. Kaveh, M.S. Massoudi / Applied Mathematics and Computation 266 (2015) 72–94

Table 7 Comparison among stresses resulting by the displacement method and the present force method for some elements of Example 2.

87

88

A. Kaveh, M.S. Massoudi / Applied Mathematics and Computation 266 (2015) 72–94

Fig. 15. Finite element models of Example 3, (a) linear, (b) quadratic, (c) cubic and (d) quartic.

numbers are calculated. In all the following examples, nnz represents the number of non-zero entries and λmax /λmin is the ratio of the extreme eigenvalues taken as the condition number of a matrix. The comparison between present algorithm and algebraic force method is shown for all three examples. Finally, the present method is validated through comparison of resulting stresses and deflections using the present graph-theoretical force method and the displacement method.

A. Kaveh, M.S. Massoudi / Applied Mathematics and Computation 266 (2015) 72–94

89

Table 8 Topological properties for the FEMs of Example 3. Element

N

M

DSI

SESs Type I

SESs Type II

SESs Type III

Linear Quadratic Cubic Quartic

481 481 481 481

432 108 48 27

1201 661 433 310

816 384 240 168

– 192 160 126

385 85 33 16

Table 9 Comparison of the optimality characteristics of the null basis matrices B1 and flexibility matrices G for the FEMs of Example 3. Element

Linear Quadratic Cubic Quartic

DSI

1201 661 433 310

DKI

962 962 962 962

TimeLU method Time

Present method

98.44 40.28 19.35 13.80

nnz(B1 )LU method nnz(B1 )Present method

6.31 5.16 4.49 8.49

Condition number λλmax

Norms max |A × B1|

min

Present method

LU method

4.96 × 10 9.56 × 102 3.71 × 102 7.76 × 103

3.20 × 10 3.34 × 104 7.92 × 103 8.06 × 105

3

4

Present method −16

4.44 × 10 4.99 × 10−15 5.32 × 10−16 2.36 × 10−14

(a)

(b)

(c)

(d)

LU method 9.99 × 10−15 8.43 × 10−15 2.88 × 10−15 6.76 × 10−10

Fig. 16. Equilibrium matrix for FEMs of Example 3, (a) linear, (b) quadratic, (c) cubic and (d) quartic.

8.1. Example 1 Consider the following finite element model comprising of 169 nodes as illustrated in Fig. 6. The model is idealized using plain stress rectangular elements. As indicated in Fig. 7, four types of elements are generated using linear, quadratic, cubic and quartic elements. All models have the same nodes with the following mechanical properties: Thickness = 0.01 m, E = 2e + 8 kN/m2 , ν = 0.3. Topological properties of the models are collected in Table 3. The pattern of equilibrium matrices for various orders are shown in Fig. 8. Fig. 9 displays pattern of the null basis matrices employing the present method and LU factorization method for the model comprising of various orders elements. The flexibility matrices are illustrated in Fig. 10. In this relation Table 4 contains other

90

A. Kaveh, M.S. Massoudi / Applied Mathematics and Computation 266 (2015) 72–94

(a)

(b) Fig. 17. Null basis matrix of Example 3 for linear, quadratic, cubic and quartic element, (a) present method, (b) LU method.

optimality characteristics of the force method procedures. It is clear that the graph theoretical method forms the most wellstructured null basis in smallest computational time. 8.2. Example 2 Consider the following finite element model comprising of 45 elements and four openings as illustrated in Fig. 11. The model is idealized using plain stress rectangular elements. As indicated in Fig. 11(c), there are four cycles in natural associated graph of FEM surround four openings. The FEM is discretized as 45 elements using linear, quadratic, cubic and quartic elements, alternatively. All models have the same members with the following mechanical properties: E = 2e + 8 kN/m2 , ν = 0.3. Topological properties of the models are collected in Table 5. The pattern of equilibrium matrices for various orders are shown in Fig. 12. Fig. 13 displays pattern of the null basis matrices employing the present method and LU factorization method for the model comprising of various orders elements. The flexibility matrices are illustrated in Fig. 14. In this relation Table 6 contains other optimality characteristics of the force method procedures. It is clear that the graph theoretical method forms the most wellstructured null basis in smallest computational time. Comparison among stresses resulting by the displacement method and the present force method for some elements 17, 22, 23, 24 and 29 is shown in Table 7. 8.3. Example 3 Consider the following finite element model comprising of 481 nodes as illustrated in Fig. 15. The model is idealized using plain stress rectangular elements. As indicated in Fig. 15, four types of elements are generated using linear, quadratic, cubic and

A. Kaveh, M.S. Massoudi / Applied Mathematics and Computation 266 (2015) 72–94

(a)

(b)

Fig. 18. Flexibility matrix of Example 3 for linear, quadratic, cubic and quartic element, (a) present method, (b) LU method.

91

92

A. Kaveh, M.S. Massoudi / Applied Mathematics and Computation 266 (2015) 72–94 Table 10 Comparison among deflections of 36 bottom nodes of the beam for Example 3 resulting by the displacement method and the present force method for some elements. X (cm)

0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36

Y-deflection (1e – 6 m) Linear Force method

Displacement method

Quadratic Force method

Displacement method

Cubic Force method

Displacement method

0 −0.04313 −0.0589 −0.07065 −0.08005 −0.08829 −0.09571 −0.10249 −0.10872 −0.11441 −0.11956 −0.12418 −0.12823 −0.13169 −0.13456 −0.13681 −0.13842 −0.1394 −0.13972 −0.1394 −0.13842 −0.13681 −0.13456 −0.13169 −0.12823 −0.12418 −0.11956 −0.11441 −0.10872 −0.10249 −0.09571 −0.08829 −0.08005 −0.07065 −0.0589 −0.04313 0

0 −0.04313 −0.0589 −0.07065 −0.08005 −0.08829 −0.09571 −0.10249 −0.10872 −0.11441 −0.11956 −0.12418 −0.12823 −0.13169 −0.13456 −0.13681 −0.13842 −0.1394 −0.13972 −0.1394 −0.13842 −0.13681 −0.13456 −0.13169 −0.12823 −0.12418 −0.11956 −0.11441 −0.10872 −0.10249 −0.09571 −0.08829 −0.08005 −0.07065 −0.0589 −0.04313 0

0 −0.05041 −0.06607 −0.07906 −0.0884 −0.09677 −0.10421 −0.11103 −0.11728 −0.123 −0.12818 −0.13282 −0.13688 −0.14037 −0.14324 −0.1455 −0.14712 −0.1481 −0.14843 −0.1481 −0.14712 −0.1455 −0.14324 −0.14037 −0.13688 −0.13282 −0.12818 −0.123 −0.11728 −0.11103 −0.10421 −0.09677 −0.0884 −0.07906 −0.06607 −0.05041 0

0 −0.05041 −0.06607 −0.07906 −0.0884 −0.09677 −0.10421 −0.11103 −0.11728 −0.123 −0.12818 −0.13282 −0.13688 −0.14037 −0.14324 −0.1455 −0.14712 −0.1481 −0.14843 −0.1481 −0.14712 −0.1455 −0.14324 −0.14037 −0.13688 −0.13282 −0.12818 −0.123 −0.11728 −0.11103 −0.10421 −0.09677 −0.0884 −0.07906 −0.06607 −0.05041 0

0 −0.05654 −0.07239 −0.08499 −0.09347 −0.10179 −0.10917 −0.11604 −0.12229 −0.12801 −0.13319 −0.13782 −0.14189 −0.14537 −0.14825 −0.15051 −0.15213 −0.15311 −0.15344 −0.15311 −0.15213 −0.15051 −0.14825 −0.14537 −0.14189 −0.13782 −0.13319 −0.12801 −0.12229 −0.11604 −0.10917 −0.10179 −0.09347 −0.08499 −0.07239 −0.05654 0

0 −0.05654 −0.07239 −0.08499 −0.09347 −0.10179 −0.10917 −0.11604 −0.12229 −0.12801 −0.13319 −0.13782 −0.14189 −0.14537 −0.14825 −0.15051 −0.15213 −0.15311 −0.15344 −0.15311 −0.15213 −0.15051 −0.14825 −0.14537 −0.14189 −0.13782 −0.13319 −0.12801 −0.12229 −0.11604 −0.10917 −0.10179 −0.09347 −0.08499 −0.07239 −0.05654 0

quartic elements. All models have the same nodes with the following mechanical properties: Thickness = 0.01m, E = 2e + 8 kN/m2 , ν = 0.3. Topological properties of the models are collected in Table 8. The pattern of equilibrium matrices for various orders are shown in Fig. 16. Fig. 17 displays pattern of the null basis matrices employing the present method and LU factorization method for the model comprising of various orders elements. The flexibility matrices are illustrated in Fig. 18. In this relation Table 9 contains other optimality characteristics of the force method procedures. It is clear that the graph theoretical method forms the most well-structured null basis in smallest computational time. The results are verified by the standard displacement method in Table 10 and Fig. 19. The sparsity of the final null basis obtained by the present method is less than 20% of the LU method.

9. Concluding remarks The main conclusions of this paper are as follows: • Solution of many examples reveals that good accuracy can be achieved by the present algorithm. • Flexibility matrices obtained are highly sparse with narrowly banded. This is due to the use of regional cycles of the natural associate graphs and appropriate ordering of the selected self-equilibrating systems. • The conditioning of the flexibility matrices generated by the present algorithm is better than those formed by the LU method. • Because of a high reduction in the number of floating point operations, the resulted null basis has better accuracy in comparison to other methods. This is obvious, since nearly 60% of the null vectors are selected without numerical analysis and the remaining null vectors are obtained with working on small and limited lists.

A. Kaveh, M.S. Massoudi / Applied Mathematics and Computation 266 (2015) 72–94

(a)

93

(b)

Fig. 19. 36 bottom nodes deflection of beam for Example 3, (a) present method, (b) displacement method.

• The method developed in this paper can easily be extended to FEMs with higher-order two-dimensional rectangular elements. It should be noted that in the present method the most important point is selecting independent element forces system. • The use higher order elements in the force method leads to fewer unknowns than the displacement method. • The computational time required for the present method is considerably lower than those of the algebraic methods. Since the complexity of the LU method is O(n3 ), if the DSI of the model increases, then the time difference dramatically rises. • In the present method, numbering the nodes of a finite element model is less important and only a suitable ordering of the members of the natural associate graph is required for reducing the bandwidth of the flexibility matrices. Acknowledgements The first author is grateful to the Iran National Science Foundation for the support. Appendix A Index FEM

Finite element model

SES A B1 Fm N M IG NAG Type I self-equilibrating system Type II self-equilibrating system Type III self-equilibrating system DSI DKI M’

Self-equilibrating system Equilibrium matrix Self-stress matrix Unassembled flexibility matrix Number of nodes of FEM Number of elements of FEM (The number of members of the NATURAL associate graph) Interface graph Natural associate graph Self-equilibrium system which is constructed on a 2-multiple member Self-equilibrium system which is extracted from two adjacent elements of FEM. Self-equilibrium system which is extracted from cycle of the NAG. Degree of statical indeterminacy Degree of kinematical indeterminacy Number of members of the associate graph of the model

References [1] J.C. de C. Henderson, Topological aspects of structural analysis, Aircr. Eng. 32 (1960) 137–141. [2] E.A.W. Maunder, Topological and linear analysis of skeletal structures, (Ph.D. thesis), London University, Imperial College, 1971. [3] J.C. de C. Henderson, E.A.W. Maunder, A problem in applied topology: On the selection of cycles for the flexibility analysis of skeletal structures, J. Instit. Math. Appl. 5 (2) (1969) 254–269. [4] A. Kaveh, Application of topology and matroid theory to the analysis of structures, Ph.D. thesis, London University, Imperial College, 1974. [5] A. Kaveh, An efficient program for generating subminimal cycle bases for the flexibility analysis of structures, Commun. Appl. Numer. Method 2 (4) (1986) 339–344. [6] A. Kaveh, A combinatorial optimization problem optimal generalized cycle bases, Comput. Methods Appl. Mech. Eng. 20 (1979) 39–52. [7] AC. Cassell, An alternative method for finite element analysis; a combinatorial approach to the flexibility method, Proc. Lond. Royal Soc. A 352 (1976) 73–89. [8] P.H. Denke, A General Digital Computer Analysis of Statically Indeterminate Structures, NASA-TD-D-1666, 1962. [9] J. Robinson, Integrated Theory of Finite Element Methods, Wiley, New York, 1973. [10] A. Topcu, A contribution to the systematic analysis of finite element structures using the force method, (Doctoral dissertation), Essen University, Germany, 1979 (in German).

94 [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35]

A. Kaveh, M.S. Massoudi / Applied Mathematics and Computation 266 (2015) 72–94 I. Kaneko, M. Lawo, G. Thierauf, On computational procedures for the force methods, Int. J. Numer. Method Eng. 18 (1982) 1469–1495. E. Soyer, A. Topc-u, Sparse self-stress matrices for the finite element force method, Int. J. Numer. Method Eng. 50 (2001) 2175–2194. J.R. Gilbert, M.T. Heath, Computing a sparse basis for the nullspace, SIAM J. Alg. Disc. Method 8 (1987) 446–459. T.F. Coleman, A. Pothen, The null space problem I; complexity, SIAM J. Alg. Disc. Method 7 (1986) 527–537. T.F. Coleman, A. Pothen, The null space problem II; algorithms, SIAM J. Alg. Disc. Method 8 (1987) 544–561. A. Pothen, Sparse null basis computation in structural optimization, Numer. Math. 55 (1989) 501–519. M.T. Heath, R. Plemmons, J. Ward, Sparse orthogonal schemes for structural optimization using the force method, SIAM J. Sci. Stat. Comput. 5 (1984) 514–532. SN. Patnaik, L. Berke, RH. Gallagher, Integrated force method versus displacement method for finite element analysis, Comput. Struct. 38 (1991) 377–407. X.F. Wei, S.N. Patnaik , Application of stochastic sensitivity analysis to integrated force method, Int. J. Stoch. Anal 2012 (2012) Article ID 249201 14 pp. G.S. Doiphode, Development of integrated force method for the analysis of frames and continuum structures, (Ph.D. thesis), The Maharaja Sayajirao University of Brado, 2012. H. Chi, T.J. Lee, Force methods for trusses with elastic boundary conditions, Int. J. Mech. Sci. 66 (2013) 202–213. V. Konstantinos, V. Spiliopoulos, N.G. Dais, A powerful force-based approach for the limit analysis of three-dimensional frames, Arch. Appl. Mech. 83 (2013) 723–742. N-I Kim, S. Lee, N. Ahn, J. Lee, Damage identification of truss structures based on force method, Adv. Appl. Math. Mech. 7 (2015) 229–244. A. Kaveh, M.J. Tolou Kian, Efficient finite element analysis using graph-theoretical force method with brick elements, Finite Elem. Anal. Des. 54 (2012) 1–15. K. Koohestani, An orthogonal self-stress matrix for efficient analysis of cyclically symmetric space truss structures via force method, Int. J. Solid Struct. 48 (2011) 227–233. A. Kaveh, M.S. Massoudi, M.J. Massoudi, Efficient finite element analysis using graph-theoretical force method; rectangular plane stress and plane strain serendipity family elements, Period. Polytech. 58 (2013) 1–20. A. Kaveh, Computational Structural Analysis and Finite Element Methods, Springer International Publishing, Springer Verlag, Switzerland, 2014. J.S. Przemieniecki, Theory of Matrix Structural Analysis, McGraw-Hill, New York, 1968. M. Kurutz, Stability of equilibrium and compatibility of structures with uniaxial material behavior, Period. Polytech. Civil Eng. 38/4 (1994) 415–430. A. Kaveh, G.R. Roosta, Comparative study of finite element nodal ordering methods, Eng. J. 20 (1–2) (1998) 86–96. A. Kaveh, S. Beheshti, Weighted triangular and circular graph products for configuration processing, Period. Polytech. Civil Eng. 56/1 (2012) 1–9. A. Kaveh, Structural Mechanics: Graph and Matrix Methods, 3rd ed., Research Studies Press (Wiley), Somerset , UK, 2004. A. Kaveh, Rahimi Bondarabady, A multi-level finite element nodal ordering using algebraic graph theory, Finite Elem. Anal. Des. 38 (2002) 245–261. A. Kaveh, H. Fazli, Analysis of frames by structuring technique based on using algebraic and graph methods, Commun. Numer. Method Eng. 23 (2007) 921–943. A. Kaveh, K. Koohestani, N. Taghizadieh, Force method for finite element models with indeterminate support conditions, Asian J. Civil Eng. 8 (4) (2007) 403–417.