Formulating the Applied Element Method: Linear 2D (Part I)

I first came to know about the Applied Element Method (AEM) during the one day seminar on "Geotechnics and Geo hazards" 2012 held in Nepal. Ramesh Guragain had mentioned this method in his presentation on the study of collapse of masonry structures (and something to do with fragility functions I believe). The results he showed using AEM were really compelling. I wanted to learn more about this method.

Here I discuss the basic procedures that I have followed to develop an AEM program on Python that can solve Linear 2D structural problems. Poisson's ratio is not considered here. All of the procedures discussed below are based on Meguro and Tagel-Den (2000). The developed program is then used to solve a classic cantilever beam problem. The displacement values obtained using AEM are compared with corresponding theoretical values and values obtained from Finite Element Method (FEM) for plane stress condition. AEM performs better compared to FEM even when the number of elements is small.

Introduction

AEM is modeled by dividing the structure into rigid elements connected with pairs of normal and shear springs as shown in figure 1. The stiffness values of each pair of springs represent the material property of certain area of both the elements as shown in figure 1 b. The stiffness values are determined as shown below:

where E and G are Young's and shear modulus, d is the distance between springs, T is the thickness of the element and a is the length of representative area.

Modeling of elements connected with springs has major advantages that include simple model and programming, and high accuracy of the results with relatively short CPU time. This is mainly because the ultimate stiffness matrix depends on the number of elements rather than the number of springs. So increasing the number of springs wont significantly affect the computing time. This model can be easily extended to model highly nonlinear behaviors including crack initiation, propagation, separation of elements, collapse, etc. This can be achieved mainly because the stiffness of each springs can be changed independently.

Figure 1: Modelling of structure to AEM (adopted from Meguro and Tagel-Den 2000)

Stiffness Matrix

First we derive the stiffness matrix for a single pair of springs connecting two rigid elements shown in Figure 2. Each of the two elements have three degrees of freedom, so the spring has a total of 6 degrees of freedom. [also note that in AEM we drive stiffness matrix for the springs and not for the element]. Similar procedure used in Chandrupatla and Belegundu (2002) section 8.7 to derive stiffness matrix of a plane frame is used to derive the stiffness matrix of the spring here. First the stiffness matrix (K'^e) is computed for the displacements in the local system (represented by q'). Then the local to global transformation matrix (L) is used to compute the stiffness matrix in the global system (q).

Figure 2: Two elements connected with a pair of normal and shear springs. 

Here K^e  represents the stiffness matrix for the single pair of springs shown in Figure 2. The governing equation for this spring pair is:

Here q is displacement vector and f is the load vector of the elements. It can be noted that, "elements" in AEM are equivalent to "nodes" in FEM. Similarly, "springs" in AEM are equivalent to "elements" in FEM. Hence, in FEM "nodes" are connected by "elements", where as in AEM "elements" are connected by "springs".

Similar procedure is followed to derive stiffness matrices for other spring pairs present. The orientation of the edge on which the springs are attached should be considered while calculating the values of φ, b_n1, b_n2, a₁, a₂, etc. It should also be noted that the values of bn1 and bn2 can be negative.

After this, all of the spring stiffness matrices are assembled to form the global stiffness matrix (K). This procedure is similar to the procedure used for plane frame elements (see Chandrupatla and Belegundu 2002).  Each row or column of the matrix represents a degree of freedom of a particular element. Penalty approach was used to incorporate boundary conditions into the global stiffness matrix. The loads on the structure will be incorporated into the global load vector (F). By solving the following governing equation, the global displacement vector (Q) can be calculated:

KQ = F

AEM Program

The object-oriented programming capability of python was used to develop the AEM program. Class definitions include:

• class Element
• Attributes:
• ele_no = element number
• edg_len = {1:a,2:b,3:a,4:b}
• b = height of element (edge 2 and 4)
• a = length of element (edge 1 and 3)
• T = thickness of element
• E = Young's modulus
• G = Shear modulus
• x = x coordinate of center of element
• y = y coordinate of center of element
• nodes = [node_1, node_2, node_3, node_4]
• r = rotation angle in radian of element
• u = displacements
• Methods:
• getPlotMat : get a matrix of x and y coordinates of the edges of element
• class Elements
• Attributes:
• list = list of Element objects
• dic = dictionary of element no linking to particular element object
• ele2id = convert element no. of element id (element id starts from 0)
• Methods:
• addElement
• addNewElements = creates new Element objects and adds it to its list. Input: element combination matrix.
• global_stiffness_matrix = creates global stiffness matrix by assembling spring stiffness matrices
• class Node
• class Nodes
• class Spring
• Attributes:
• ele_1 = element one
• ele_2 = element two
• edg_1 = edge no. of element one
• edg_2 = edge no. of element two
• d = height of spring
• T = thickness of spring
• K_n = Normal spring constant
• K_s = Shear spring constant
• x_1 = x coordinate of first point of spring
• y_1 = y coordinate of first point of spring
• x_2 = x coordinate of second point of spring
• y_2 = y coordinate of second point of spring
• Methods:
• stiffness_matrix = creates spring stiffness matrix
• class Springs
• class Load

Validation of program

Figure 3

In order to validate the program, a simple cantilever beam bending problem shown in Figure 3 was used. Mechanical properties of the cantilever beam that were used are given below:

Young's Modulus (E) = 207 GPa
Shear Modulus (G) = 79.1 GPa
Thickness (T) = 0.15 m
Moment of Area (I) = 0.0001 m⁴

The theoretical deflection of the cantilever beam at point B was calculated as shown below (for derivation  see article 21 of Timoshenko and Goodier 1951):

Figure 4: Deflection of the cantilever beam using the AEM program that was developed on Python. The deflection is magnified 500 times.

The displacement at B calculated using AEM and FEM for the different models is presented in Table 1. For FEM modeling Abaqus 6.13 was used. For the AEM model, 10 springs per edge was used.  It can be seen that AEM performs better than FEM even when the number of elements is small.

Table 1: Displacement at B computed for the different models shown in Figure 4 using both FEM and AEM.

Model	Total no. of elements	No. of elements per section	using FEM		using AEM
			Displacement	Percentage Error	Displacement	Percentage Error
a	10	1	-0.0448225	99.71%	-0.000131577	1.13%
b	38	2	-0.000172269	24.49%	-0.00013108	0.76%
c	136	4	-0.000138103	5.80%	-0.00013087	0.60%
d	300	6	-0.000133241	2.37%	-0.000130843	0.58%

Acknowledgements

Special thanks to Dr. Shanker Dhakal for his valuable comments and suggestions that has helped improve the quality of this article.

Reference

• Chandrupatla TR and Belegundu AD (2002) Introduction to Finite Elements in Engineering, Third Edition. Prentice Hall.
• Meguro K and Tagel-din H (2000) Applied Element Method for Structural Analysis: Theory and Application for Linear Material. Structural Eng./Earthquake Eng., JSCE, Vol. 17, No. 1, 21s-35s
• Popov EP and Balan TA (2010) Engineering Mechanics of Solids. PHI Learning, New Delhi
• Timoshenko SP and Goodier JN (1951) Theory of elasticity. New York, McGraw-Hill