Numerical Modelling in Geo engineering

CIE4366 Numerical Modelling in GEO-ENGINEERING

Coursework 3:

2‐Dimensional Finite Element Method Programming

Name: Student Nr:

CX. Azua-Gonzalez 4490894

INDEX 1. ABSTRACT............................................................................................................................... 3 2. INTRODUCTION...................................................................................................................... 3 3. THEORETICAL BACKGROUND............................................................................................... 4 a. Skyline Storage.................................................................................................................. 4 b. Cholesky decomposition................................................................................................... 4 c.

Isoparametric elements.................................................................................................... 6

d. Bending of Euler-Bernoulli beams................................................................................... 6 e. Equivalent nodal loading systems...................................................................................7 4. MATHEMATICAL FORMULATIONS......................................................................................... 8 a. Gauss-Legendre quadrature............................................................................................. 8 b. Shape functions and local derivatives..........................................................................10 c.

Jacobian matrices............................................................................................................. 13

d. Global derivatives............................................................................................................ 13 e. Constitutive relations...................................................................................................... 14 f.

[B] matrix.......................................................................................................................... 15

g. Stresses............................................................................................................................. 15 5. GENERAL TWO-DIMENSIONAL FEM FORMULATION........................................................16 6. VALIDATION OF RESULTS.................................................................................................... 23 a. Axial compression of a column(laterally restrained).................................................23 b. Axial compression of a column(laterally unrestrained)............................................24 c.

Bending of a cantilever beam under concentrated load...........................................27

7. ANALYSIS OF RESULTS........................................................................................................ 28 8. CONCLUSIONS...................................................................................................................... 29 9. REFERENCES......................................................................................................................... 30 2

10.

APPENDIX.......................................................................................................................... 30

3

1.

ABSTRACT

A two general FEM program was developed in Fortran 95 as part of the assignment 3A and 3B for the course CIE 4366 CIE Numerical Modelling in Geo Engineering. Different types of elements were incorporated including linear, a quadratic elements. This code was used to solve simple two- dimensional structural mechanics problems under plain strain formulation, and one- dimensional structural elements under compression and bending. A standard FEM procedure was proposed for the code especially for the calculation of the element stiffness matrices, and calculation of stresses and internal forces. The standardized procedure in combination with the elaboration of subroutines was proven to assist in writing a clearer code to the reader. The main program was modified to assist the user in generating automatic node numbering element numbering. Furthermore, validation of the results was done by means of simplified analytical solutions. Finally, the influence of the vicinity to restraints, type of element, and number of elements was studied.

Keywords: FEM, Gauss-Lengendre quadrature, skyline storage, solid mechanics

2.

INTRODUCTION

Two dimensional structural analyses of solid bodies can be performed by means of the Finite Element Method (FEM), and its complexity can be controlled by the “available” degrees of freedom in each formulation. It might be intuitive that the mathematical formulation of each element type could influence the results, especially if approximate integration procedures are used. In fact, the use of Gauss-Legendre quadrature appears to be common for solving approximate integrals under the domain of dimensionless local coordinates. For these purpose, the Integrals are evaluated in the so called integration points derived after the Legendre Polynomials [1]. The current report aims to explore the use of a Fortran 95 Code to perform a structural analysis of two-dimensional solid bodies using Gauss-Legendre quadrature and classical formulations of stresses, and strains in mechanics. The program solves the problems in displacements as primary variables; as a result, boundary conditions are regarded as fixed nodes. Nodal forces are applied in the framework of equivalent external work, which means all the forces are applied on the nodes in such a manner that the external work would be equivalent as that of the original loading system. Additionally, one- dimensional solid bodies are solved to make a comparison among the results achieved by a two- dimensional analysis and analytical solutions if available. 4

The goals of this assignment are: 1. Explain how finite element programs work. 2. Use basic scientific programming techniques, including limited modification of FEM programs. To achieve these goals, a general FEM code was written partially and delivered to the students of the course CIE 4366 CIE Numerical Modelling in Geo Engineering. The task of the students was to complete the missing parts of the code using available subroutines. Additional subroutines were asked to be elaborated where needed. Finally, it should be mentioned that the structure of the program delivered in this report aims to find an acceptable solution. In fact, the code does not represent the most efficient manner to solve a structural analysis of solid bodies in terms of computational resources. Nevertheless, the implementation of skyline storage procedures of element stiffness matrices into a global stiffness matrix was required to make the problems resolvable by means of a standard computer.

3.

THEORETICAL BACKGROUND a. Skyline Storage

Under the situation of extremely sparse matrices the procedure of skyline storage is necessary. This process implies the storage of the values which are taken into account in the resolution of the matrix linear system. The procedure aims to be an alternative to old versions of banded storage [2]. The procedure requires the storage of the values in a skyline matrix and the position of the main diagonal terms. An elucidation of the skyline storage procedure is depicted blow:

Figure 1. Elucidation of skyline storage procedure Furthermore, it can be seen that sometimes it is useful to start the position vector at zero (depending on the further calculation procedures adopted.)

b.

Cholesky decomposition 5

Cholesky decompositions use the symmetry and positive definiteness of the stiffness matrix as a starting point of the analysis. In fact, if the previous condition is true, this procedure requires one half of the floating point operations than Gauss elimination [3]. There are two varians of cholesky decomposition: 1.

K ¿ BT B

2.

K ¿ BT DB

In general, the second variant of Cholesky decomposition is preferred to avoid the calculation of square roots. For that reason, this procedure is explained below, after defining matrices as depicted:

Equation 1. Definition of matrices in Cholesky decomposition Step1. Decompose

K ¿ BT DB

A. For each i=1, j-1 in increments of 1 a. H=kij b. Bij=h/di c. For each k=i+1,j in increments of 1 kkj is replaced by kkj-hbik B. dj=ajj Step2. Forward substitution

f ¿ BT z , Db=z

for each j=1,n in increments of

A. ZJ=FJ B. For each i=1,j-1 in increments of 1 zj:=zj-bijzi C. bj=zj/dj Step3. Back substitution

Bu=b

for each j=n,1 in increments of -1

A. uj=bj B. for each i=j+1,n in increments of 1 uj:=uj-bjiui The determinant of K is calculated with

|K|=|BT||D||B|=|D|=( d 1 d 2 … . d n )

6

Figure 2. Elucidation of the second Cholesky decomposition variant [3]

c.

Isoparametric elements

The formulation of interpolation functions depends on the application of the element (degrees of freedom). Nevertheless, in two dimensional structural analysis of solid bodies the use of the shape functions, often called interpolation functions, can be extended to interpolate coordinates from a local system to a global system. The elements which posses this kind of formulation are called isoparametric elements. The convenience of using isoparametric elements lie in the unique formulation of the shape function and derivatives to elaborate the Strain-displacement matrix [B]. If the chain rule is applied to the partial derivatives of the shape functions, a relationship between derivative in the local system and in the global system can be found [1]. The term that relates these derivatives is called the jacobian. The mathematical relationship is depicted below:

Equation 2. Jacobian matrix formulation

d. Euler-Bernoulli beams

Bending of

Bending is produced by transverse loading in a beam. The relationship of the curvature k, the modulus of elasticity and the moment of inertia shown in figure below, can be used to back calculate internal bending moments in a FEM formulation, after the resolution of the problem in terms of the primary variables (transverse displacement and rotational slope theta). Stresses can be calculated as the multiplication of the curvature, the modulus of elastic and the distance with respect to the neutral axis.

7

N.B: The fortran code developed in this report calculates the multiplication Ek automatically, however the user must perform the calculation of EKy to calculate the stress at a desired distance.

Figure 3. Elucidation of stresses, curvature and bending moments in an Euler-Bernoulli Beam [4].

e.

Equivalent nodal loading systems

The formulation of element stiffness matrices, and loading matrices is based upon the simplification of loading systems. In fact, distributed loads in a surface or defined as volume loads (loads per unit area, or body loads applied by gravity). Different procedures can be applied to achieve this goal, variational approaches are used to find equivalent nodal loads which ensure the same external work as the original system. Nevertheless, this procedure might not be the most common due to its mathematical formulation. In contrast, structural engineers in their desire of getting good estimation of results tend to use the strip method. This procedure consists in dividing distributed loads over its area of influence for each node as it is shown below. The procedure tends to be called a node by node lumping [2].

8

Figure 4. Elucidation of the node by node direct lumping [2]. In that sense, when the internal forces are calculated in a system resolved by means of an equivalent system the real forces must be computed subtracting the equivalent forces produced over each element.

{ f internal }=[ k m ] { u } −{f equivalent } Equation 3. Computation of real internal forces

N.B: The program code developed in the current report calculates the product of [km] and {u}, and the user is in charge of subtracting the {fequivalent} in case that the nodal nodes were derived from another loading system. This scheme gives the users of the code the freedom of defining the loads via variational approach or other procedures by themselves. The examples shown in this report, are based on simplified direct nodal conditions in which {fequivalent} is a zero vector. Hence, internal forces shown in the appendix are final answers and do not need an extra manipulation.

4.

MATHEMATICAL FORMULATIONS

The summary of mathematical formulations, and equations to be solved are depicted in this section of the report. It is assumed that the reader has a basic background of linear algebra, and solid mechanics.

a.

Gauss-Legendre quadrature 9

An integral can be represented by its approximate summation of the function evaluated in certain coordinates if its domain goes from -1 to 1. This is an advantage for FEM analysis, due to the situation that dimensionless local coordinates are defined in that domain. The weightings com from a mathematical approximation of Legendre polynomials of the same order into the expression shown below. The derivation of these coordinates and coefficients is not considered under the scope of this report. Nevertheless, it is important to mention that there is a minimum number of terms required in the summation, this values apply for a uni-variable polynomial and can be extrapolated to multi-variable polynomials as a the one for bilinear rectangular elements. The coordinates and weighting coefficients are shown in the chart below.

Equation 4. Gauss-Legendre Quadrature

chart 1. Gauss-Legendre Quadrature coordinates and weighting coefficients [3]

10

Figure 5. Degree of the polynomial for typical quadrilateral elements [3] For the integration of the [B] matrix, the correct number of points to be incorporated in the Gauss-Legendre Quadrature is essential. For that reason, a preliminary analysis of the nature of the shape functions must be performed. Even without knowing the shape function formulas, but its general definition in powers the maximum order of the polynomial can be foreseen. The B matrix collects partial derivatives and decreases the order of the polynomial in every case, except for the quadratic quadrilateral element. Virtually, the product T

[ B ] [ B ] rises this exponentials to the square as it can be observed in the figure shown.

Figure 6. Degree of the polynomials in the process of integration [3]

b.

Shape functions and local derivatives 11

A summary of the shape functions and partial derivatives for the elements used in the Fortran code are shown below. Quadrilateral elements are regarded with two degrees of freedom per node, while the rod elements are treated with only one degree of freedom. Exceptionally one dimensional beams have degrees of freedom per node. 2

3 4

1

Figure 7. Distribution of local numbering of nodes in a Bilinear Quadrilateral element

[ ][

N1 0.25 ( 1−ξ ) (1−η ) N2 0.25 ( 1−ξ ) ( 1+ η ) = N3 0.25 ( 1+ξ ) (1+ η ) N4 0.25 ( 1+ξ )( 1−η )

]

.

Equation 5. Shape functions for a Bilinear Quadrilateral element

[ ] ∂N1 ∂ξ ∂N2 ∂ξ ∂ N3 ∂ξ ∂ N4 ∂ξ

∂ N1 ∂η ∂ N2 −0.25 ( 1−η ) −0.25 ( 1−ξ ) −0.25 ( 1+η ) 0.25 ( 1−ξ ) ∂η = ∂ N3 0.25 ( 1+ η ) 0.25 ( 1+ξ ) ∂η 0.25 ( 1−η ) −0.25 ( 1+ξ ) ∂ N4 ∂η

[

]

.

Equation 6. Local derivatives for a Bilinear Quadrilateral element 4 3 2

5 6

1

7 8

Figure 8. Distribution of local numbering of nodes in a Quadratic Quadrilateral element

12

[ ][ ] 0.25 (1−ξ )( 1−η ) (−ξ−η−1 ) N1 0.5 ( 1−ξ ) ( 1−η 2) N2 0.25 ( 1−ξ ) ( 1+ η ) (−ξ +η−1 ) N3 N4 0.5 ( 1−ξ 2) ( 1+η ) = N5 0.25 ( 1+ξ )( 1+η )( ξ +η−1 ) N6 0.5 ( 1+ξ ) ( 1−η2 ) N7 0.25 (1+ ξ ) ( 1−η )( ξ−η−1 ) N8 0.5 ( 1+ξ ) ( 1−η2 )

.

Equation 7. Shape functions for a Quadratic Quadrilateral element

[] ∂N1 ∂ξ ∂N2 ∂ξ ∂ N3 ∂ξ ∂ N4 ∂ξ ∂ N5 ∂ξ ∂ N6 ∂ξ ∂ N7 ∂ξ ∂ N8 ∂ξ

∂ N1 ∂η ∂ N2 ∂η 0.25 ( 1−η ) ( 2 ξ+ η ) 0.25 ( 1−ξ ) ( 2 η+ ξ ) ∂ N3 −0.5 ( 1−η2 ) −η ( 1−ξ ) ∂η 0.25 ( 2ξ−η )( 1+η ) 0.25 ( 1−ξ ) ( 2η−ξ ) ∂ N4 −ξ ( 1+η ) 0.5 ( 1−ξ2 ) ∂η = ∂ N5 0.25 ( 1+η ) ( 2 ξ+ η ) 0.25 ( 2 η+ξ )( 1+ξ ) ∂η 0.5 ( 1−η2 ) −η ( 1+ ξ ) ∂ N6 0.25 ( 2ξ−η )( 1−η ) 0.25 ( 1+ξ ) ( 2η−ξ ) ∂η −ξ ( 1−η ) −0.5 ( 1−ξ2 ) ∂ N7 ∂η ∂ N8 ∂η

[

]

.

Equation 8. Local derivatives for a Quadrative Quadrilateral element

1 2 Figure 9. Distribution of local numbering for a linear rod element and a one dimensional

beam

[ ][

N 1 0.5 ( 1−ξ ) = N2 0.5 ( 1+ ξ )

]

. 13

Equation 9. Shape functions for a linear rod or two-node beam element (spatially)

[ ][

∂ N1 ∂ ξ = −0.5 ∂ N2 0.5 ∂ξ

]

.

Equation 10. Local derivatives for a linear rod or a two-node beam element (spatially) 1 3

2

Figure 10. Distribution of local numbering for a quadratic rod

[ ][

N 1 −0.5 ξ (1−ξ ) N 2 = ( 1−ξ ) (1+ ξ ) N3 0.5 ξ ( 1+ξ )

]

.

Equation 11. Shape functions for a quadratic rod

[]

∂ N1 ∂ξ −0.5 ( 1−2 ξ ) ∂ N2 = −2 ξ ∂ξ 0.5 ( 1+2 ξ ) ∂ N3 ∂ξ

[

]

.

Equation 12. Local derivatives for a quadratic rod N.B.: One dimensional beams in the current Fortran code are treated as non isoparametric elements. Thus, the shape functions used for a rod element, were used only for the spatial interpolation of the one-dimensional beam element. The shape functions to interpolate the transverse displacement of a beam are depicted below. Note that there are two degrees of freedom in each node despite the one dimensional formulation: 1 2

Figure 11. Distribution of local numbering and degrees of freedom in a one dimensional beam element

[ ][

2 0.25 ( 1−ξ ) ( 2+ξ ) N 1v 2 N 1θ 0.125 l (1−ξ ) ( 1+ξ ) = 2 N 2v 0.25 ( 1+ξ ) (2−ξ ) N 2θ −0.125 l ( 1+ξ )2 (1−ξ )

]

.

Equation 13. Shape functions for a one dimensional beam(interpolation of {U}) 14

[] ∂2 N 1 v

∂ ξ2 ∂2 N 1 θ

[

1.5 ξ ∂ ξ2 0.25 l ( 3 ξ−1 ) = 2 ∂ N2v −1.5 ξ −0.25 l ( 3 ξ +1 ) ∂ ξ2

]

.

∂2 N 2 θ ∂ ξ2

Equation 14. Local second order derivatives N.B.: The curvature-displacement matrix [B] for a Beam is elaborated with global second order shape function derivatives which are derived by means of the Jacobian matrix from the local second order derivatives as shown in the just previous equation.

c.

Jacobian matrices

A summary of the Jacobian matrices for the elements used in the Fortran code are shown:

[

∂ N1 ∂ξ J= ∂ N1 ∂η

∂ N2 ∂ξ ∂ N2 ∂η

∂ N3 ∂ξ ∂ N3 ∂η

][ ]

∂ N 4 x1 ∂ ξ x2 ∂ N 4 x3 ∂η x 4

y1 y2 y3 y4

.

Equation 15. Jacobian matrix of a bilinear quadrilateral element

[

∂ N1 ∂ξ J= ∂ N1 ∂η

∂ N2 ∂ξ ∂ N2 ∂η

∂ N3 ∂ξ ∂ N3 ∂η

∂N4 ∂ξ ∂N4 ∂η

∂ N5 ∂ξ ∂ N5 ∂η

∂N6 ∂ξ ∂N6 ∂η

∂ N7 ∂ξ ∂ N7 ∂η

∂N8 ∂ξ ∂N8 ∂η

]

[] x1 x2 x3 x4 x5 x6 x7 x8

y1 y2 y3 y4 y5 y6 y7 y8

.

Equation 16. Jacobian matrix of a quadratic quadrilateral element J=

[

∂ N1 ∂ξ

][ ]

∂ N 2 x1 ∂ ξ x2

.

Equation 17. Jacobian matrix of a linear rod or one-dimensional beam 15

[

∂ N1 J= ∂ξ

∂ N2 ∂ξ

][

x ∂ N3 1 x2 ∂ξ x3

]

.

Equation 18. Jacobian matrix of a quadratic rod

d.

Global derivatives

The global derivatives are obtained by means of the inverse of the jacobian matrix. A especial case is presented for the elaboration of global second order derivatives for onedimensional formulation of an Euler-Bernoulli beam. In the later case, the chain inversion procedure is used twice. It can be seen straight forward that because the Jacobian of a two node element is filled with constants, the inverse of the Jacobian is filled with constants and gets out of the derivative in the second inversion step.

[

∂ N1 ∂x ∂ N1 ∂y

∂ N2 ∂x ∂ N2 ∂y

∂N3 ∂x ∂N3 ∂y

][

∂ N4 ∂x =J −1 ∂ N4 ∂y

∂ N1 ∂ξ ∂ N1 ∂η

∂ N2 ∂ξ ∂ N2 ∂η

∂N3 ∂ξ ∂N3 ∂η

∂ N4 ∂ξ ∂ N4 ∂η

]

.

Equation 19. Global derivatives for a bilinear quadrilateral element

[

∂ N1 ∂x ∂ N1 ∂y

∂ N2 ∂x ∂ N2 ∂y

∂ N3 ∂x ∂ N3 ∂y

∂ N4 ∂x ∂ N4 ∂y

∂ N5 ∂x ∂ N5 ∂y

∂N6 ∂x ∂N6 ∂y

][

∂ N7 ∂x ∂ N7 ∂y

∂N8 ∂x =J −1 ∂N8 ∂y

∂ N1 ∂ξ ∂ N1 ∂η

∂ N2 ∂ξ ∂ N2 ∂η

∂ N3 ∂ξ ∂ N3 ∂η

∂ N4 ∂ξ ∂ N4 ∂η

∂ N5 ∂ξ ∂ N5 ∂η

. Equation 20. Global derivatives for a quadratic quadrilateral element

[

∂ N1 ∂x

] [

∂ N2 ∂N1 = J −1 ∂x ∂ξ

∂ N2 ∂ξ

]

.

Equation 21. Global derivatives for a linear rod

[

∂ N1 ∂x

∂ N2 ∂x

] [

∂N3 ∂ N1 =J −1 ∂x ∂ξ

∂ N2 ∂ξ

∂ N3 ∂ξ

]

.

Equation 22. Global derivatives for a quadratic rod

[

∂2 N 1 v

∂2 N 1 θ

∂2 N 2 v

∂2 N 2 θ

∂ x2

∂ x2

∂ x2

∂ x2

] [ =( J −1 )

2

∂2 N 1 v

∂2 N 1 θ

∂2 N 2 v

∂2 N 2θ

∂ ξ2

∂ ξ2

∂ ξ2

∂ ξ2

]

.

Equation 23. Global second order derivatives for a one-dimensional beam 16

∂N6 ∂ξ ∂N6 ∂η

∂ N7 ∂ξ ∂ N7 ∂η

∂N8 ∂ξ ∂N8 ∂η

]

N.B: The square of the inverse of the Jacobian is conceptually valid for this particular case given the fact that the Jacobian for a linear rod or a two node one dimensional beam element is a scalar.

e.

Constitutive relations

In general, the constitutive relation matrix varies with its application. For the current general FEM program the formulation of hyperelastic materials was considered. In that, sense plasticity is not taken into account nor failure envelopes to lock the stress state of infinitesimally small particles in the continuum. Due to the addition of the terms of Area, and moment of Inertia during the integration of the element stiffness matrices, the constitutive relation [D] in the Fortran code can be seen to be merged with the Area, or moment of Inertia. This appears to be convenient in the standardization of the Finite element procedure. Nevertheless, the original formulation of [D] as in mechanics is considered for the computation of stresses.

[ D ]=

[

1

E(1−v ) v ( 1+ v )( 1−2 v ) 1−v 0

v 1−v

0

1

0

0

1−2 v 2 ( 1−v )

]

.

Equation 24. Constitutive Elasticity matrix for plain strain formulation of twodimensional solid bodies

[ D ] =E . Equation 25. Constitutive relation for one-dimensional rod and beam elements

f. [B] matrix The derivation of the [B] matrix for different solid mechanics problems is not considered as part of the scope of this report. The [B] matrix called strain-displacement matrix for mechanical problems of two- dimensional elements for 4(8)-node quadrilateral elements, and one dimensional rod elements is presented. Additionally the curvature-displacement matrix [B] for one dimensional bending problems is shown as well.

17

[ [

∂ N1 ∂ N2 0 0 ∂x ∂x ∂ N1 ∂ N2 [ B ]= 0 0 ∂y ∂y ∂ N 1 ∂ N1 ∂ N2 ∂ N 2 ∂y ∂x ∂y ∂x

∂ N3 ∂ N4 0 0 ∂x ∂x ∂ N3 ∂ N4 0 0 ∂y ∂y ∂ N3 ∂ N3 ∂ N4 ∂ N 4 ∂y ∂x ∂y ∂x

]

.

Equation 26. Strain-displacement matrix for a Bilinear Quadrilateral element. ∂ N1 ∂ N2 0 0 ∂x ∂x ∂ N1 ∂ N2 [ B ]= 0 0 ∂y ∂y ∂ N 1 ∂ N1 ∂ N2 ∂ N 2 ∂y ∂x ∂y ∂x

∂N3 ∂x

… .0

∂ N7 ∂y ∂ N7 …. ∂x

0 …. ∂ N3 ∂y

∂N8 0 ∂x ∂ N8 0 ∂y ∂ N 8 ∂ N8 ∂y ∂x

]

.

Equation 27. Strain –displacement matrix for a Quadratic Quadrilateral element

[

[ B ]=

∂ N1 ∂x

∂ N2 ∂x

].

Equation 28. Strain-displacement matrix for a linear rod element

[

[ B ]=

∂ N1 ∂x

∂ N2 ∂x

∂ N3 ∂x

].

Equation 29. Strain-displacement matrix for a quadratic rod element

[

[ B ]=

∂2 N 1 v ∂ x2

∂2 N 1θ ∂ x2

∂2 N 2 ∂ x2

∂2 N θ ∂ x2

]

.

Equation 30. Curvature displacement matrix for one-dimensional beam element

g.

Stresses

Stresses within this Fortran code are calculated at the integration points. Conceptually, if the displacements are calculated by means of the resolution of the system of equations, the strain field can be calculated at any point inside the domain of the elements. Nevertheless, it can be foreseen that the mathematical formulation of the shape functions influence greatly in the estimation of stresses since some shape functions are 18

more “rigid” than others. In that sense, quadrilateral elements are preferred over one dimensional elements due to the fact that the derivatives of the shape functions that define the gradient change over the domain, which is not the case for rod elements or even linear triangular elements. Note that stresses in a beam element on the distance from the neutral axis at which the stresses are calculated. The program computes de contribution of the Elasticity matrix [D] and curvature {k}. This result must be multiplied by the distance “y” at which the user of the program requires. σ =[ D ] { ε }=[ D ] [B] { U } . Equation 31. Stresses for a two dimensional plain strain formulation, or one dimensional rod elements σ =[ D ] { k } y= [ D ] [B ] {U } y . Equation 32. Stresses for one dimensional beam elements

5. GENERAL TWO-DIMENSIONAL FEM FORMULATION The general equation to be solved has the form as depicted below:

[ A ] T [ D ][ A ][ S ] { U }=−{ f }

.Equation 33

Due to the Galerkin’s procedure to reduce the residuals the previous equation lies as depicted below for a two-dimensional solid mechanics problem:

∫ [ S ] T [ A ] T [ D ][ A ][ S ] dx dy { U } =∫ [ S ]T dx dy { f }

. Equation 34

The product of the diferential elements in [A] and the shape functions in [S] lies in the generation of the [B] matrix. This matrix is called the strain-displacement matrix for solid problems, and curvature-displacement matrix for bending of beams. The general expression shown above is modified for the case of one-dimensional prismatic rod elements as depicted below:

∫ A [ S ] T [ A ]T [ D ][ A ][ S ] ds { U }=∫ A [ S ] T ds { f }

. Equation 35

In the case that a prismatic beam element is analyzed the expression lies as depicted below:

∫ Iner [ S ] T [ A ] T [ D ][ A ][ S ] ds { U }=∫ Iner [ S ] T ds { f }

. Equation 36 19

The expressions show clearly the standardized procedure to be taken. In fact, the major difference in the procedure for different problems is the decision in how the right hand side of the equations is considered, under traditional variational approach or simplified methodologies to find the equivalent nodal loads. Using the Gauss-Legendre Quadrature the element stiffness matrices lie as depicted below: nip

T

[ k m ]=∑ [ B ] [ D ] [ B ] det ( J ) W i . i=1 Equation 37. Element stiffness matrix for a general two dimensional FEM problem nip

T [ k m ]=∑ A [ B ] [ D ] [ B ] det ( J ) W i i=1

.

Equation 38. Element stiffness matrix for one-dimensional prismatic FEM problem nip

T [ k m ]=∑ Iner [ B ] [ D ] [ B ] det ( J ) W i i=1

.

Equation 39. Element stiffness matrix for one-dimensional prismatic bending problem

a. Conceptual Flow chart of the resolution in Fortran PROGRAM FE -Use library -Initialization The first step lies in the inititalization of variables. Dynamic arrays are preferred to set up the correct size and save computational resources. Here the variables are read from the input data file, and more important the number of equations and size of the [loads] vector is determined in combination with the position of diagonal elements in kdiag.

-Calculate Stiffness matrix The general procedure to calculate the [kv] stiffness matrix in skyline storage is shown in a descriptive manner: Step1. Set [kvacum] to zero. Which is a temporary variable for accumulation of [kv] constributions is=n skyline storage

Loop the elements from iel=1,nels 20

Step2. Set the [km] matrix to zero Step3. The coordinates, length (in the case of a beam) of the elements are obtained Subroutine: fcoord, flel Step4. The constitutive relation matrix [D] is calculated. Storage for further use. Subroutine: fconst

Loop the integration points DO inip=1,nip Step5.The local derivatives of the shape functions are obtained. If a one dimensional beam problem is being solved, the shape functions of a linear rod work for the jacobian of the beam. Nevertheless, the shape functions for interpolating {U} is different. If the reader requires the revision of the formulas. Subroutine: bending problem.

flocderive, and flocsecderiv in the case of a

Step6. The multiplication of the local derivatives and the coordinates is made to compute the Jacobian. Subroutine: fjac N.B.: The Jacobian matrix of a linear rod can be used for a onedimensional Euler Bernoulli beam. This comes from the fact that both elements are defined geometrically by two nodes. Nevertheless, the shape functions to interpolate the transverse are different, thus, the beam element can no longer be regarded as an isoparametric element. Step7. The inverse of the Jacobian matrix is calculated to change the reference system, from the local to the global system. If a bending problem is taken into account, the double inversion procedure must be performed. Here the function of inverse(matrix) from the library is used. Subroutine: fglobderiv, fglosecderiv in the case of a bending problem is activated automatically. Step8. The determinant of the Jacobian matrix is determined. The function determinant(matrix) is used inside the subroutine. Subroutine: fdet 21

Step9.The [B] matrix is calculated. The input values into this subroutine must be the description of the element and the global derivatives (second order in case of bending problem). Subroutine: fstrainsdisp Step10. The contribution of the inip-th integration point into the element stiffness is accumulated in [km]. The subroutine retrives the prebious value of [km], and the components for integration ([B], [D],det(J),Wi and A or Iner if required). Each [B] is stored for further use. Step 11. Output partial results END DO (for integration points)

Step 12. Calculate steering vector of the element, and store for further use. Store [km] for further use. Step13. Set to Zero [kv] to avoid undesirable previous values. Compute new contribution to [kv]. Subroutine: fsparv collects g, km and kdiag (computed in the intitialization stage) Step 14. Accumulate new [kv] in [kvacum](temporary name for the variable [kv] to avoid problems with the use of the subroutine sparv from the library) END DO (for elements) Step 15. Store [kacum] in [kv] for further usage in the code.

-Loading conditions Step 16. Retrieve the loads from the input data file. As soon as possible, close the space (10) in which the data file was being read. Nodes without force components are left with zeros.

-Equation solution Step 17. System of equations is solved and solutions are stored in [loads]. Subroutines: Sparin and spabac

-Output results Step 18. Displacements are re-stored in a global displacement matrix including zeros for fixed freedoms 22

Subroutine: fglobdisp

Loop the elements DO iel=1,nels Step 19. Call back [km] and retrieve the coordinates of the element. Afterwards the coordinates are used to compute the location of the integration points in the global system of reference Subroutine: fcoord Step 20. If a one dimensional rod or bending problem is being analyzed [D] matrix is restored as in mechanics that is [D]=AE/A or [D]=E(Iner)/(Iner). This is done due to the fact that [D] and the Area or moment of Inertia were merged in [D] to gain a more standard procedure for any kind of problem.

Loop the integration points DO inip=1,nip Step 21. Call back [B] to be used for calculation of stresses at Gaussian points Step 22. Retrieve the displacements at the nodes of the element from global array Subroutine: feldisp Step 23. Calculate the stresses as [D][B]{U}, in the case of bending the results will be normalized by the distance to the neutral axis. The user of the program must be aware of evaluating this expression at the neutral axis. It should be recalled that at y=0 from the neutral axis of a beam no stresses are developed. Step 24. The shape functions are evaluated at the local coordinates corresponding to the location of the inip-th integration point. Step 25. The global coordinates ate which the stresses were calculated are computed as the matrix multiplication of the shape functions and the global coordinates. Subroutine: fconv Step 26. Output of the stresses, and location at which the inip-th integration point at which calculations were performed END DO (for the integration points)

Step 27. Calculate internal forces without subtraction of equivalent nodal forces. The user is in charge of identifying if equivalent nodal forces were 23

used. If the later is true, the real internal forces are obtained by subtracting the equivalent nodal forces from the result of the program. The internal forces are obtained as the matrix multiplication of [Km] and {U}. N.B: If only nodal forces were considered in the original loading system, the result of the program is definitive. Revise equation 3 for more details. Step 28. Output of the internal forces END DO (for elements)

-Visualization of results Step 29. Generate the mesh using the global coordinates array and the global numbering array. The file generated will be named as “input data file.msh” Subroutine: mesh Step 30. Generate the deformed mesh with the name “input data file.dis”. The [loads] vector with solutions, and array with freedoms nf are used. Subroutine: dismesh Step 31. Generate the displacement vector file for visualization. The [loads] vector with solutions, and array with freedoms nf are used. This file will be named as “input data file.vec” Subroutine: vecmsh

STOP Contains SUBROUTINES

END PROGRAM FE

b. Automatic node and element numbering The previous Fortran code was slightly modified during the initialization stage. Complex schemes of meshing found in literature as in [5], in which Delauney triangulation is used, suggest that triangular elements are useful to force meshes into complex geometries. Nevertheless, due to the limited time to perform this assignment and regarding the lower rigidity of the shape functions in quadrilateral elements a simple scheme of automatic node and element numbering. 24

Figure 12. Elucidation of forced meshing techniques by means of Delauney triangulation

-Initialization Step1.The first step lies in the inititalization of variables. The number of columns and rows of elements (nx,ny) are read form the input data file. The number of node columns and node rows are calculated (nnx, nny). Additionally, the vectors of guiding coordinates xcoord and ycoord are allocated the correct size. nnx=nx +1 . Equation 40

nny=ny +1 . Equation 41 nn=nnx∗nny . Equation 42 Step 2. Check the axis in which the smallest number of nodes is required. If equal number of nodes is required, numbering will be performed along the x- direction, which means the row is kept constant until another row is involved. N.B.: A triple loop is proposed having a dummy variable k that increases the node numbering while the loops keep going. A similar procedure is applied to number the corresponding elements. Here the repetitive pattern along the moveable axis is taken as an advantage to relate the position of the local nodes (1,2,3,4) to the global numbering as shown below.

if (nx>ny) then !Numbering along y axis !Numbering of nodes k=1 DO WHILE (k

g_coord(1,k)=xcoord(1,i) g_coord(2,k)=ycoord(1,j) k=k+1 END DO END DO END DO !Numbering of elements k=1 DO WHILE (k
26

Step 3. The previous Fortran code continues without changes.

6.

VALIDATION OF RESULTS a. Axial compression of a column(laterally restrained) 2 1

3

4

Figure 13. Graphical representation of the sample problem, fixed basement and laterally restrained The sample problem given to solve was related to solid body, laterally restrained in the nodes and in the base. The problem in the sample output was related to a unique element. The same task was performed with the Fortran code attached to this report. The results matched until the 4th digit after the decimal point. Nevertheless, the solution provided must have been proved to be right. As result, and analytical solution was derived under the theory of elasticity. This problem was treated as a Biaxial test in which the lateral displacement is imposed to be zero. The analytical solution for the stresses is depicted below. For the problem v=0.3 and

∆ σ y =¿

F/A=-0.5*2/1.0=-1.0. As a result the lateral stress

must be ∆ σ x=

v 0.3 ∆ σ y= (−1.0 )=−0.4286 1−v 1−0.3

Results given by the Program match with this analytical solution with a minor error: Solutio n Analyti cal Progra

∆σx 0.4286 0 -

∆σ y

Error % -

-1.0

Error % -

0.01%

-0.74

26% 27

m

0.4285 7 chart 2. Comparison of solutions for the sample problem The internal forces at the base nodes was another proved that was made to this problem. In fact, the internal forces of the element at the basement must be equal to the total load on top. The linear distributed forces on the basement are

F2 =F 4 =0.5 As it can be seen in the input data file, the internal forces of this unique element are in fact the reaction forces that equilibrate the system.

b. Axial compression of a column(laterally unrestrained)

Figure 14. Column under axial compression, Bilinear quadrilateral elements (more than one element) The same problem was solved by means of a 4 element column. Nevertheless, the boundary conditions were changed to explore the stability of the calculation schemes as shown in the input data file sample. The chosen element type was a bilinear quadrilateral (4 nodes per element). The problem can be solved by means of a simplified system in which only axial compression is taken into account, for that analytical solution the displacement on top the column was computed:

∆ y=∆

σ∗y −1.0∗1 = =−1.0E-6 E 1.0E6

Solutio n Analyti cal Progra m

σy -1.0

Error % -

-1.0

0.0%

∆y

Error % -

-1.0E6 -0.879 12%

28

chart 3. Comparison of solutions for a 4 bilinear element column laterally unrestrained The results are at least from the same order of magnitude regarding displacements. However, it should be mentioned that perhaps the assumption of only axial interaction was a bit biased and so the results of the program are still ok. If the visualization tool of the program is used, the displacement vectors in fact appear to be slightly tilted to the outer part of the column. This should have been the first assumption for the analytical solution because of the two-dimensional solution of the FEM scheme. Below the mesh and displacement vectors generated by the program are shown.

Figure 15. Undeformed and deformed mesh, displacement vectors for bilinear element column laterally unrestrained 3 2

5 6

1

7

Figure 16. Equivalent nodal forces imposed for quadratic quadrilateral elements (more than one element) A row of 4 quadratic elements was used to analyze the column problem (unrestrained laterally). The equivalent nodal forces were calculated based on the Node by Node lumping. As a result, the middle node gets half of the total load on top of the column, while the nodes on the corners receive one fourth of the total force each. The situation explained is elucidated in the Figure 16. -0.12326E-05

∆y ∆σ y Solutio Error% Error n % Analyti -1.0E-6 -1.0 cal Progra -1.23E- 23% -1.4 40% m 6 chart 4. Comparison between approximate analytical solution and program results for 4 quadrilateral elements unrestrained laterally. 29

Figure 17. Undeformed, deformed mesh, displacement vectors for quadratic rectangular elements laterally unrestrained (more than one element). Another important aspect to be taken into account is that the restrained basement must be playing a great influence in the displacement field of the column particles. To explore that matter, the modified version of the Program was used to create a denser mesh. The results show that the displacement field starts to reduce as the location approaches to the restrained basement. In that sense, the analytical assumption could be valid only for the top particles of the body. The preliminary conclusion that the restraints affect the distribution of stresses arises so far.

Solutio n Analyti cal Progra m

∆y -1.E-6

Error % -

-0.74E- 26% 6 -0.87E- 13% 6 -0.2E-6 80% 0.0E+0 chart 5. Comparison of solutions for an irregular element array column laterally unrestrained (more than one element)

Figure 18. Mesh and displacement vectors for an irregular element array column laterally unrestrained 30

As a final procedure to test the sensibility of the problem additionally to the vicinity to restrained boundaries and the distribution of elements, the type of element was changed. The same problem was analyzed by means of one dimensional rod elements to observe how much it deviates from the “refined solution” of a denser array. Surprisingly, the result approaches exactly as the simplified solution. It must be mentioned that the program works with a two point integration Gauss-Legendre quadrature for rods. Linear polynomials tend to be “rigid” over its domain in terms of strains due to the fact that derivatives of shape functions are constant along this domain. However, reality might be slightly different and the strain fields vary along the domain of the element.

Solution

∆y

Error%

∆σ y

Error % 0.0%

Analytical -1.0E-6 -1.0 Program linear -1.0E-6 0.0% -1.0 rod Program -1.0E-6 0.0% -1.0 0.0% quadratic rod chart 6. Comparison of solutions for linear and quadratic rod element array column laterally unrestrained (more than one element)

An important aspect to be mentioned, its that both the linear and quadratic formulation of the rods gave the same results and converged to the analytical solution. These calculations can be assumed to be exact due to the amount of integration points that were used. However, it should be noticed that not always having a higher number of integration terms will allow getting stable results. The results, in fact should be interpreted critically otherwise wrong outcomes might arrive. The analytical solutions present for the unrestrained column are a very simplified calculation and does not represent the reality. The assumption for the analytical solution lies in one-dimensional compression of a solid body. As a result, the analytical solution and the

c. Bending of a cantilever beam under concentrated load

31

Figure 19. Cantilever beam under a concentrated load at the free extreme The problem of a beam perfectly embedded in one extreme, under a concentrated load on the free extreme was analyzed by means of an equivalent plain strain formulation, which means that an infinite slab is treated as adjacent beam disregarding the friction between the portions of the slab. Exact solutions for the maximum slope and transverse displacement of the beam at the free extreme are available and were used to validate the results of the program. The equivalent distributed load of -200.0 over a width of 0.5 is equivalent to -100.0 in an Euler-Bernoulli beam in a one-dimensional formulation The analytical solution considering the static equilibrium shows that the bending moment should be equal to M=200. In that sense the stress at extreme attached to the wall might be calculated and compared to the result given in the program. One stress point (Gauss point) was picked randomly to check its value against an analytical solution as depicted below:

σ=

My 200∗(2−0.0528) = ∗( 0.25−0.0132 ) =8850.23 Iner 2

Solution

x

y

Analytical

0.052 8 0.052 8

0.0132

σy 8850.23

Error % -

Program Quadrilateral 0.0132 8733.00 1.4% element chart 7. Comparison of a cantilever two dimensional solid beam and analytical solutions.

Figure 20. Deformed mesh, and displacement vectors of a cantilever beam. 32

N.B: The beam was modeled as a two dimensional solid body, in that sense the freedom of rotation is somehow involved. Furthermore, both solutions match in good accuracy. Finally, to explore the effect of the shape functions, one dimensional beam elements were used instead of a two dimensional element. The maximum slope and transverse displacement was computed from existing solutions and compared with the results given by the program.

−100 ( 22 ) −P l 2 θmax = = =−0.19194e-3 2 EIner 2 ( 200e6 ) (5.21e-3 )

v max =

3 −100 ( 23 ) −P l = =−2.5592e-4 3 EIner 3 ( 200e6 )( 5.21e-3 )

Solution

θmax

v max

Error%

Analytical

−0.19194e-3

−2.5592e-4

-

−0.19194e-3 −2.5592e-4 Program Quadrilateral 0.0% element chart 8. Comparison of a cantilever one dimensional solid beam and analytical solutions. As it can be seen, results from a one dimensional approach match with the analytical solution. Nevertheless, the rigidity in terms of the gradient produced by the shape functions might differ from those gradients produced in reality.

7.

ANALYSIS OF RESULTS

In general quadratic shape functions tend to give smoother results in terms of the stress gradients produced by the shape functions inside the domain of the elements. Nevertheless, the use of larger number of nodes might be cumbersome while preparing the input data file. This issue was overcome with the generation of automatic node and element numbering inside the code for quadrilateral elements. Additionally, checking internal forces at extreme elements was proved to be a good tool of validation of the results. From the group of sample input files analyzed by means of the Fortran code, it could be seen that at least the order of magnitude of stresses, and displacements must match with analytical solutions, no matter how simplified this analytical solutions are assumed as long as engineering judgment is used to test the assumptions. Furthermore, the effect of the number of nodes could be observed. The higher the number of nodes, the better the visualization outputs tend to be. These graphical outputs are a good tool to check if the boundary conditions are correctly applied especially if symmetry is dominant. The node by node lumping appeared to be effective in assisting in the load matrix components of the loading matrix. However, its effectiveness was only tested in loading systems with nodal forces, and calculation of reaction forces might get more difficult for the user due to the fact that the program does not distinguish between equivalent and nodal forces. This issue might be 33

overcome with a slight modification of the code, making it to calculate equivalent nodal forces by itself. Moreover, the use of skyline storage schemes was proved to be helpful in making a standard computer solve Two-dimensional problems without the risk of exceeding its available memory. The computational resources required for producing more nodes is definitely imperceptible for the kind of problems that were solved in this report. In that sense, it is highly valuable to retrieve a more realistic solution under the solution of two dimensional FEM equations. Finally, it can be mentioned that the “standard” procedure used for solving FEM equations tends to assist the programmer into the generation of shorter and clearer codes. In addition, the generation of specific subroutines and functions helped to re-organized in a clean manner the main program structure

8.

CONCLUSIONS

a. The use of the Jacobian matrix as a conversion path from local to global coordinate system or the other way around, tend to give good results if the appropriate number of integration points is chosen. b. Quadratic elements tend to be more reliable in terms of accurate and realistic gradients of stresses and displacements. This might be acknowledged to the higher order of the polynomials in the shape functions, and thus in the stress gradients obtained from the integration of the element stiffness matrix c. The restraints appear to affect the integration points in its vicinity in terms of stresses. As a result, displacement vectors seem to be considerably reduced if located in the surroundings of a fixed boundary. For that reason, while analyzing a two- dimensional problem the boundary conditions must be sufficiently far away to ensure that no restrictions are being imposed due to this vicinity, exception lies on the case in which this restriction is imposed in the real system. d. The use of a standard FEM procedures helps in generating more compact and clear codes. This practice enables the continuation of further modification of the code. e. The automatic node and element numbering proposed in the current report helps in saving time to the user and avoids undesirable mistakes produced during the writing input. f. The validation of FEM results appear to be necessary in addition to the verification of variables. This procedure of validation can be performed by means of checking the reactions, internal forces and comparison with analytical solutions. Additionally, solving the same problem with equivalent loads, and different types of elements appeared to be a good technique to check at least the order of magnitudes of the results.

34

9.

REFERENCES

[1] Moaveni S. (1999).Finite Element analysis: Theory and application with ANSYS. Prentice Hall [2] Opensource: Introduction to Finite Element Methods. Fall 2015. Department of Aerospace engineering sciences. University of Colorado at Boulder. http://www.colorado.edu/engineering/CAS/courses.d/IFEM.d/ [3]Kress G. Structural Analysis with FEM. Composite material and adaptive structures lecture notes. May 2014 [4] Opensource: Bending moments and beam curvatures. University of Cambdrigde http://www.doitpoms.ac.uk/tlplib/beam_bending/bend_moments.php [5] Per-olof P, and Gilbert S. A simple mesh generator in MATLAB [6] I.M. Smith and D.V. Griffiths . Programming the finite element method, 4th edition, John Wiley & Sons Limited, 2004.

35

10.

10. APPENDIX A. B. C.

INPUT SAMPLES: PART 1(A-D), PART 2(A-B), PART 3(A-B) RESULTS: PART 1(A-D), PART 2(A-B), PART 3(A-B) FORTRAN CODE PART 1, PART 2

36

INPUT SAMPLES: PART 1-“LINEAR ELEMENTS INPUT DATA FILES”: PART 1A: -1 QUADRILATERAL ELEMENT UNDER AXIAL COMPRESSION (LATERALLY RESTRAINED) PART 1B: -4 QUADRILATERAL ELEMENTS UNDER AXIAL COMPRESSION (LATERALLY UNRESTRAINED)

PART 1C: -4 ROD ELEMENTS IN A COLUMN UNDER AXIAL COMPRESSION ON TOP PART 1D: -4 BEAM ELEMENTS IN A CANTILEVER BEAM UNDER POINT LOAD AT

PART 2-“QUADRATIC ELEMENTS INPUT DATA FILES”: PART 2A: -4 ELEMENT COLUMN (2D CUADRATIC QUADRILATERAL ELEMENT) PART 2B: -4 ELEMENT COLUMN (1D CUADRATIC ROD ELEMENTS)

PART 3-“AUTOMATIC MESHING INPUT DATA FILES”: PART 3A: - IRREGULAR ELEMENT COLUMN (2D BILINEAR QUADRILATERAL ELEMENT) PART 3B: - IRREGULAR ELEMENT CANTILEVER BEAM (2D BILINEAR QUADRILATERAL ELEMENTS) 37

PART 1-LINEAR ELEMENTS INPUT DATA FILES. N.B.: Bold messages stand for explanation matters only.

PART 1A: 1 QUADRILATERAL ELEMENT UNDER AXIAL COMPRESSION (LATERALLY RESTRAINED): element type, no. nodes 'quadrilateral' 4 no. elements, no. nodes (total), no. ips (per element), no of freedoms per node, no. stress terms, no of dimensions 1 4 4 2 3 2 np_types 1 Prop (E, v) 1.0e6 0.3 Global nodal coordinates 0.0 0.0 1.0 0.0 0.0 -1.0 1.0 -1.0 Global element node numbering (global node numbers as local nodes) 3 1 2 4 Global freedoms (nf, node number 1=free 0=fixed) 4 1 0 1 2 0 1 3 0 0 4 0 0 Loaded nodes (nn, node number, load x load y) 2 1 0.0 -0.5 2 0.0 -0.5

PART 1B: 4 QUADRILATERAL ELEMENTS UNDER AXIAL COMPRESSION (LATERALLY UNRESTRAINED): element type, no. nodes 'quadrilateral' 4 no. elements, no. nodes (total), no. ips (per element), no of freedoms per node, no. stress terms, no of dimensions 4 10 4 2 3 2 np_types 1 Prop (E, v) 1.0e6 0.3 Global nodal coordinates 0.0 0.0 1.0 0.0 0.0 -0.25 1.0 -0.25 0.0 -0.5 1.0 -0.5 0.0 -0.75 1.0 -0.75 0.0 -1.0 1.0 -1.0 Global element node numbering (global node numbers as local nodes) 3 1 2 4 5 3 4 6 7 5 6 8 9 7 8 10 Global freedoms (nf, node number 1=free 0=fixed) 2 9 0 0 10 0 0 Loaded nodes (nn, node number, load x load y)

38

PART 1C: 4 ROD ELEMENTS IN A COLUMN UNDER AXIAL COMPRESSION ON TOP: element type, no. nodes 'rod' 2 no. elements, no. nodes (total), no. ips (per element), no of freedoms per node, no. stress terms, no of dimensions 4 5 2 1 1 1 np_types 1 Prop (E, A) 1.0e6 1.0 Global nodal coordinates 0.0 -0.25 -0.5 -0.75 -1.0 Global element node numbering (global node numbers as local nodes) 2 1 3 2 4 3 5 4 Global freedoms (nf, node number 1=free 0=fixed) 1 5 0 Loaded nodes (nn, node number, load y) 1 1 -1.0

PART 1D: 4 BEAM ELEMENTS IN A CANTILEVER BEAM UNDER POINT LOAD AT THE EXTREME: element type, no. nodes 'beam' 2 no. elements, no. nodes (total), no. ips (per element), no of freedoms per node, no. stress terms, no of dimensions 4 5 2 2 1 1 np_types 1 Prop (E, Iner) 200.0e6 52.1e-4 Global nodal coordinates 0.0 0.50 1.0 1.5 2.0 Global element node numbering (global node numbers as local nodes) 1 2 2 3 3 4 4 5 Global freedoms (nf, node number 1=free 0=fixed) 1 1 0 0 Loaded nodes (nn, node number, load y) 1 5 -100.0 0.0

39

PART 2-QUADRATIC ELEMENTS INPUT DATA FILES. N.B.: Bold messages stand for explanation matters only.

PART 2A: 4 QUADRILATERAL ELEMENTS UNDER AXIAL COMPRESSION (LATERALLY UNRESTRAINED): element type, no. nodes 'quadrilateral' 8 no. elements, no. nodes (total), no. ips (per element), no of freedoms per node, no. stress terms, no of dimensions 4 23 9 2 3 2 np_types 1 Prop (E, v) 1.0e6 0.3 Global nodal coordinates 0.0 0.0 0.5 0.0 1.0 0.0 0.0 -0.125 1.0 -0.125 0.0 -0.25 0.5 -0.25 1.0 -0.25 0.0 -0.375 1.0 -0.375 0.0 -0.5 0.5 -0.5 1.0 -0.5 0.0 -0.625 1.0 -0.625 0.0 -0.75 0.5 -0.75 1.0 -0.75 0.0 -0.875 1.0 -0.875 0.0 -1.0 0.5 -1.0 1.0 -1.0 Global element node numbering (global node numbers as local nodes) 6 4 1 2 3 5 8 7 11 9 6 7 8 10 13 12 16 14 11 12 13 15 18 17 21 19 16 17 18 20 23 22 Global freedoms (nf, node number 1=free 0=fixed) 3 21 0 0 22 0 0 23 0 0 Loaded nodes (nn, node number, load y) 3

PART 2B: 4 ROD ELEMENTS IN A COLUMN UNDER AXIAL COMPRESSION ON TOP: element type, no. nodes 'rod' 3 no. elements, no. nodes (total), no. ips (per element), no of freedoms per node, no. stress terms, no of dimensions 4 9 2 1 1 1 np_types 1 Prop (E, A) 1.0e6 1.0 Global nodal coordinates 0.0 -0.125 -0.25 -0.375 -0.5 -0.625 -0.75 -0.875 -1.0 Global element node numbering (global node numbers as local nodes) 3 2 1 5 4 3 7 6 5 9 8 7 Global freedoms (nf, node number 1=free 0=fixed) 1 9 0 Loaded nodes (nn, node number, load y) 1 1 -1.0 40

PART 3- AUTOMATIC MESHING INPUT DATA FILES N.B.: Bold messages stand for explanation matters only.

PART 3A: IRREGULAR ELEMENT COLUMN (2D BILINEAR QUADRILATERAL ELEMENT) 0625element type, no. nodes 'quadrilateral' 4 no. elements in x-,y- direction, no. ips (per element), no of freedoms per node, no. stress terms, no of dimensions 8 16 4 2 3 2 np_types 1 Prop (E, v) 1.0e6 0.3 Global x- coordinates, Global y- coordinates to create mesh 0.0 0.125 0.25 0.375 0.5 0.625 0.75 0.875 1.0 -0.0 -0.0625 -0.125 -0.1875 -0.25 -0.3125 -0.375 -0.4375 -0.5 -0.5625 -0.625 -0.6875 -0.75 -0.8125 -0.875 -0.9375 -1.0 Global freedoms (nf, node number 1=free 0=fixed) 9 145 0 0 146 0 0 147 0 0 148 0 0 149 0 0 150 0 0 151 0 0 152 0 0 153 0 0 Loaded nodes (nn, node number, load y) 9 1 0.0 -0.0625 2 0.0 -0.125 3 0.0 -0.125 4 0.0 -0.125 5 0.0 -0.125 6 0.0 -0.125 7 0.0 -0.125 8 0.0 -0.125 9 0.0 -0.

PART 3B: IRREGULAR ELEMENT CANTILEVER BEAM (2D BILINEAR QUADRILATERAL ELEMENTS) element type, no. nodes 'quadrilateral' 4 no. elements in x-,y- direction, no. ips (per element), no of freedoms per node, no. stress terms, no of dimensions 8 8 4 2 3 2 np_types 1 Prop (E, v) 1.0e6 0.3 Global x- coordinates, Global y- coordinates to create mesh 0.0 0.25 0.5 0.75 1.0 1.25 1.50 1.75 2.0 -0.0 -0.0625 -0.125 -0.1875 -0.25 -0.3125 -0.375 -0.4375 -0.5 Global freedoms (nf, node number 1=free 0=fixed) 9 1 0 0 10 0 0 19 0 0 28 0 0 37 0 0 46 0 0 55 0 0 64 0 0 73 0 0 Loaded nodes (nn, node number, load y) 9 9 0.0 -12.5 18 0.0 -25.0 27 0.0 -25.0 36 0.0 -25.0 45 0.0 -25.0 54 0.0 -25.0 63 0.0 -25.0 72 0.0 -25.0 81 0.0 -12.5 41

RESULTS: PART 1-“LINEAR ELEMENTS OUTPUT DATA FILES: PART 1A: -1 QUADRILATERAL ELEMENT UNDER AXIAL COMPRESSION (LATERALLY RESTRAINED) PART 1B: -4 QUADRILATERAL ELEMENTS UNDER AXIAL COMPRESSION (LATERALLY UNRESTRAINED)

PART 1C: -4 ROD ELEMENTS IN A COLUMN UNDER AXIAL COMPRESSION ON TOP PART 1D: -4 BEAM ELEMENTS IN A CANTILEVER BEAM UNDER POINT LOAD AT

PART 2-“QUADRATIC ELEMENTS OUTPUT DATA FILES”: PART 2A: -4 ELEMENT COLUMN (2D CUADRATIC QUADRILATERAL ELEMENT) PART 2B: -4 ELEMENT COLUMN (1D CUADRATIC ROD ELEMENTS)

PART 3-“AUTOMATIC MESHING OUTPUT DATA FILES”: PART 3A: - IRREGULAR ELEMENT COLUMN (2D BILINEAR QUADRILATERAL ELEMENT) PART 3B: - IRREGULAR ELEMENT CANTILEVER BEAM (2D BILINEAR QUADRILATERAL ELEMENTS) 42

RESULTS: PART 1A 1 QUADRILATERAL ELEMENT UNDER AXIAL COMPRESSION (LATERALLY RESTRAINED)

43

1 QUADRILATERAL ELEMENT UNDER AXIAL COMPRESSION (LATERALLY RESTRAINED): There are

2 equations and the skyline storage is

3

Element stiffness matrix and assemblage of global stiffness matrix:

Element No. 1 integration point 1 shape function derivatives with respect to local coordinates: -0.39433756470680237 -0.39433756470680237 -0.10566243529319763 0.39433756470680237 0.10566243529319763 0.10566243529319763 0.39433756470680237 -0.10566243529319763 jacobian matrix: 0.5 0. 0. 0.5 determinant of the jacobian matrix: 0.25 shape function derivatives with respect to global coordinates: -0.7886751294136047 -0.7886751294136047 -0.21132487058639526 0.7886751294136047 0.21132487058639526 0.21132487058639526 0.7886751294136047 -0.21132487058639526 strain displacement matrix, B, ip: 1 -0.7886751294136047 0. -0.7886751294136047 0. -0.7886751294136047 -0.7886751294136047 -0.21132487058639526 0. 0.7886751294136047 0. 0.7886751294136047 -0.21132487058639526 0.21132487058639526 0. cumulative km matrix, km 269138.27585577377 149521.26436431878 -3718.7611492722344 -73687.11730531833 -72115.38590972821 -40064.10328318234 -193304.12879677335 -35770.04377581811 149521.26436431878 269138.27585577377 -35770.04377581811 -193304.12879677335 -40064.10328318234 -72115.38590972821 -73687.11730531833 -3718.7611492722344 -3718.7611492722344 -35770.04377581811 74837.70798123218 -40064.10328318234 996.4390777682769 9584.554640913797 -72115.38590972821 66249.59241808664 -73687.11730531833 -193304.12879677335 -40064.10328318234 213623.82789161906 19744.404188336644 51795.686814882516 94006.81640016403 -72115.38590972821 -72115.38590972821 -40064.10328318234 996.4390777682769 19744.404188336644 19323.26001707742 10735.1444539319 51795.686814882516 9584.554640913797 -40064.10328318234 -72115.38590972821 9584.554640913797 51795.686814882516 10735.1444539319 19323.26001707742 19744.404188336644 996.4390777682769 -193304.12879677335 -73687.11730531833 -72115.38590972821 94006.81640016403 51795.686814882516 19744.404188336644 213623.82789161906 -40064.10328318234 -35770.04377581811 -3718.7611492722344 66249.59241808664 -72115.38590972821 44

integration point 2 shape function derivatives with respect to local coordinates: -0.10566243529319763 -0.39433756470680237 0.39433756470680237 0.39433756470680237 0.10566243529319763 -0.10566243529319763

-0.39433756470680237 0.10566243529319763

jacobian matrix: 0.5 0. 0. 0.5 determinant of the jacobian matrix: 0.25 shape function derivatives with respect to global coordinates: -0.21132487058639526 -0.7886751294136047 -0.7886751294136047 0.7886751294136047 0.7886751294136047 0.21132487058639526 0.21132487058639526 -0.21132487058639526 strain displacement matrix, B, ip: 2 -0.21132487058639526 0. -0.7886751294136047 0. -0.7886751294136047 -0.21132487058639526 -0.7886751294136047 0. 0.7886751294136047 0. 0.7886751294136047 -0.7886751294136047 0.7886751294136047 0. 0.21132487058639526 0. 0.21132487058639526 0.7886751294136047

cumulative km matrix, km 343975.9838370059 189585.36764750112 -7437.52229854447 -37917.07352950022 -144230.77181945642 -106313.69570126898 -192307.68971900508 -45354.598416731904 189585.36764750112 482762.10374739283 37917.07352950022 -386608.2575935467 -134070.91968334638 -144230.77181945642 -93431.52149365497 48076.92566561028 -7437.52229854447 37917.07352950022 343975.983837006 -189585.36764750112 -192307.68971900508 45354.598416731904 -144230.77181945642 106313.69570126898 -37917.07352950022 -386608.2575935467 -189585.36764750112 482762.1037473929 93431.52149365497 48076.92566561028 134070.91968334638 -144230.77181945642 -144230.77181945642 -134070.91968334638 -192307.68971900508 93431.52149365497 232947.08790869647 50799.24773711424 103591.37362976503 -10159.84954742285 -106313.69570126898 -144230.77181945642 45354.598416731904 48076.92566561028 50799.24773711424 94160.96799830958 10159.849547422848 1992.8781555365535 -192307.68971900508 -93431.52149365497 -144230.77181945642 134070.91968334638 103591.37362976503 10159.849547422848 232947.0879086965 -50799.24773711424 -45354.598416731904 48076.92566561028 106313.69570126898 -144230.77181945642

45

integration point 3 shape function derivatives with respect to local coordinates: -0.10566243529319763 -0.10566243529319763 0.10566243529319763 0.39433756470680237 0.10566243529319763 -0.39433756470680237

-0.39433756470680237 0.39433756470680237

jacobian matrix: 0.5 0. 0. 0.5

determinant of the jacobian matrix: 0.25 shape function derivatives with respect to global coordinates: -0.21132487058639526 -0.21132487058639526 -0.7886751294136047 0.21132487058639526 0.7886751294136047 0.7886751294136047 0.21132487058639526 -0.7886751294136047 strain displacement matrix, B, ip: 3 -0.21132487058639526 0. -0.21132487058639526 0. -0.21132487058639526 -0.21132487058639526 -0.7886751294136047 0. 0.21132487058639526 0. 0.21132487058639526 -0.7886751294136047 0.7886751294136047 0. 0.7886751294136047 0. 0.7886751294136047 0.7886751294136047 0.21132487058639526 cumulative km matrix, km 363299.2438540833 200320.51210143304 44358.164516338045 -28332.518888586426 -216346.15772918463 -146377.79898445134 -191311.2506412368 -25610.19422839526 200320.51210143304 502085.3637644702 57661.47771783687 -385611.8185157784 -174135.02296652872 -216346.15772918463 -83846.96685274118 99872.6124804928 44358.164516338045 57661.47771783687 557599.8117286251 -229649.47093068346 -385611.81851577846 -28332.518888586426 -216346.15772918463 200320.51210143304 -28332.518888586426 -385611.8185157784 -229649.47093068346 557599.811728625 57661.47771783687 44358.164516338045 200320.51210143304 -216346.15772918463 -216346.15772918463 -174135.02296652872 -385611.81851577846 57661.47771783687 502085.36376447027 200320.512101433 99872.6124804928 -83846.96685274118 -146377.79898445134 -216346.15772918463 -28332.518888586426 44358.164516338045 200320.512101433 363299.2438540834 -25610.19422839526 -191311.2506412368 -191311.2506412368 -83846.96685274118 -216346.15772918463 200320.51210143304 99872.6124804928 -25610.19422839526 307784.79588992865 -90863.35102029658 -25610.19422839526 99872.6124804928 200320.51210143304 -216346.15772918463 -83846.96685274118 -191311.2506412368 -90863.35102029658 307784.79588992865

46

integration point 4 shape function derivatives with respect to local coordinates: -0.39433756470680237 -0.10566243529319763 -0.10566243529319763 0.10566243529319763 0.10566243529319763 0.39433756470680237 0.39433756470680237 -0.39433756470680237 jacobian matrix: 0.5 0. 0. 0.5 determinant of the jacobian matrix: 0.25 shape function derivatives with respect to global coordinates: -0.7886751294136047 -0.21132487058639526 -0.21132487058639526 0.21132487058639526 0.21132487058639526 0.7886751294136047 0.7886751294136047 -0.7886751294136047 strain displacement matrix, B, ip: 4 -0.7886751294136047 0. -0.21132487058639526 0. -0.21132487058639526 -0.7886751294136047 -0.21132487058639526 0. 0.21132487058639526 0. 0.21132487058639526 -0.21132487058639526 0.21132487058639526 0. 0.7886751294136047 0. 0.7886751294136047 0.21132487058639526 0.7886751294136047 0. -0.7886751294136047 0. -0.7886751294136047 0.7886751294136047

cumulative km matrix, km 576923.0717457024 240384.61538461538 96153.85133122056 -48076.92307692307 -288461.54363891284 -240384.61538461538 -384615.37943801016 48076.92307692307 240384.61538461538 576923.0717457023 48076.92307692307 -384615.3794380101 -240384.61538461538 -288461.54363891284 -48076.92307692308 96153.85133122056 96153.85133122056 48076.92307692307 576923.0717457025 -240384.61538461538 -384615.37943801016 -48076.92307692307 -288461.54363891284 240384.61538461538 -48076.92307692307 -384615.3794380101 -240384.61538461538 576923.0717457024 48076.92307692307 96153.85133122056 240384.61538461538 -288461.54363891284 -288461.54363891284 -240384.61538461538 -384615.37943801016 48076.92307692307 576923.0717457024 240384.61538461535 96153.85133122056 -48076.92307692308 -240384.61538461538 -288461.54363891284 -48076.92307692307 96153.85133122056 240384.61538461535 576923.0717457025 48076.92307692307 -384615.37943801016 -384615.37943801016 -48076.92307692308 -288461.54363891284 240384.61538461538 96153.85133122056 48076.92307692307 576923.0717457024 -240384.61538461538 48076.92307692307 96153.85133122056 240384.61538461538 -288461.54363891284 -48076.92307692308 -384615.37943801016 -240384.61538461538 576923.0717457024 47

Cumulative global stiffness matrix in skyline storage: 576923.0717457024

96153.85133122056

576923.0717457025

kdiag: [ 1 3 ]T

Results: Node

x-disp

y-disp

1

0.00000E+00

-0.74286E-06

2

0.00000E+00

-0.74286E-06

3

0.00000E+00

0.00000E+00

4

0.00000E+00

0.00000E+00

Stresses at Gaussian points and Internal Forces: Element no. 1 The integration point(nip= 4 ) stresses are:

Element

x-coord

y-coord

sig_x

sig_y

tau_xy

1

0.2113

-0.7887

-0.42857E+00

-0.10000E+01

-0.92462E-17

1

0.2113

-0.2113

-0.42857E+00

-0.10000E+01

-0.19976E-16

1

0.7887

-0.2113

-0.42857E+00

-0.10000E+01

-0.19976E-16

1

0.7887

-0.7887

-0.42857E+00

-0.10000E+01

-0.92462E-17

The node (nod= 4 ) internal forces are: Element 1 1 1 1

Local number 1 2 3 4

Global number 3

Int. Hor. Force

0.21429E+00

Int. Vert. Force 0.50000E+00

1

0.21429E+00

2

-0.21429E+00

-0.50000E+00

4

-0.21429E+00

0.50000E+00

48

-0.50000E+00

RESULTS: PART 1B 4 QUADRILATERAL ELEMENTS UNDER AXIAL COMPRESSION (LATERALLY UNRESTRAINED)

49

4 QUADRILATERAL ELEMENTS UNDER AXIAL COMPRESSION (LATERALLY UNRESTRAINED): There are 16 equations and the skyline storage is 88

Element stiffness matrix and assemblage of global stiffness matrix: Cumulative global stiffness matrix in skyline storage: 624999.9943911778 -240384.61538461538 1826923.0605280576 144230.77483959153 48076.92307692307 624999.9943911777 -48076.92307692307 865384.6317796343 240384.61538461535 1826923.0605280576 -456730.7636219469 -48076.92307692307 -312500.0056088223 -240384.61538461538 1249999.9887823556 48076.92307692307 -1778846.1374511349 -240384.61538461538 -913461.5548565574 0. 3653846.1210561153 -312500.0056088223 240384.61538461538 -456730.76362194685 48076.92307692307 288461.54967918305 -7.275957614183426E-12 1249999.988782355 240384.61538461538 -913461.5548565574 -48076.92307692308 -1778846.1374511349 0. 1730769.2635592686 -2.9103830456733704E-11 3653846.1210561153 -456730.7636219469 -48076.92307692307 -312500.0056088223 -240384.61538461538 1249999.9887823556 48076.92307692307 -1778846.1374511349 -240384.61538461538 -913461.5548565574 0. 3653846.1210561153 -312500.0056088223 240384.61538461538 -456730.76362194685 48076.92307692307 288461.54967918305 -7.275957614183426E-12 1249999.988782355 240384.61538461538 -913461.5548565574 -48076.92307692308 -1778846.1374511349 0. 1730769.2635592686 -2.9103830456733704E-11 3653846.1210561153 -456730.7636219469 -48076.92307692307 -312500.0056088223 -240384.61538461538 1249999.9887823556 48076.92307692307 -1778846.1374511349 -240384.61538461538 -913461.5548565574 0. 3653846.1210561153 -312500.0056088223 240384.61538461538 -456730.76362194685 48076.92307692307 288461.54967918305 -7.275957614183426E-12 1249999.988782355 240384.61538461538 -913461.5548565574 -48076.92307692308 Results: Node

x-disp

y-disp

1

-0.19434E-06

-0.87849E-06

2

0.19434E-06

-0.87849E-06

3

-0.19331E-06

-0.65124E-06

4

0.19331E-06

-0.65124E-06

5

-0.18698E-06

-0.42478E-06

6

0.18698E-06

-0.42478E-06

7

-0.15548E-06

-0.20237E-06

8

0.15548E-06

-0.20237E-06

9

0.00000E+00

10

0.00000E+00

50

0.00000E+00 0.00000E+00

Stresses at Gaussian points and Internal Forces: Element no. 1 The integration point(nip= 4 ) stresses are: Element tau_xy 1 -0.91643E-03

x-coord

y-coord

0.2113

sig_x

-0.1972

sig_y

-0.33812E-02

-0.10003E+01

1

0.2113

-0.0528

-0.17775E-02

-0.99966E+00

-0.91643E-03

1

0.7887

-0.0528

-0.17775E-02

-0.99966E+00

0.91643E-03

1

0.7887

-0.1972

-0.33812E-02

-0.10003E+01

0.91643E-03

The node (nod= 4 ) internal forces are: Element

Local number

1

Global number

1

Int. Hor. Force

Int. Vert. Force

3

0.64484E-03

0.50000E+00 -0.50000E+00

1

2

1

-0.28027E-16

1

3

2

0.22765E-16

-0.50000E+00

Element no. 2 The integration point(nip= 4 ) stresses are:

Element

x-coord

y-coord

sig_x

sig_y

tau_xy

2

0.2113

-0.4472

-0.15593E-01

-0.10021E+01

-0.56244E-02

2

0.2113

-0.3028

-0.57502E-02

-0.99789E+00

-0.56244E-02

2

0.7887

-0.3028

-0.57502E-02

-0.99789E+00

0.56244E-02

2

0.7887

-0.4472

-0.15593E-01

-0.10021E+01

0.56244E-02


Local number

Global number

Int. Hor. Force

Int. Vert. Force

2

1

5

0.33127E-02

0.50000E+00

2

2

3

-0.64484E-03

-0.50000E+00

2

3

4

0.64484E-03

-0.50000E+00

2

4

-0.33127E-02

0.50000E+00

6

51


Element tau_xy

x-coord

y-coord

sig_x

sig_y

3

0.2113

-0.6972

-0.76723E-01

-0.10105E+01

-0.27977E-01

3

0.2113

-0.5528

-0.27763E-01

-0.98951E+00

-0.27977E-01

3 0.27977E-01

0.7887

-0.5528

-0.27763E-01

-0.98951E+00

3 0.27977E-01

0.7887

-0.6972

-0.76723E-01

-0.10105E+01


Local number

3

1

3

2

3

Global number 7 5

3

6

Int. Hor. Force

Int. Vert. Force

0.16373E-01

0.50000E+00

-0.33127E-02

-0.50000E+00

0.33127E-02

-0.50000E+00


Element tau_xy 4 -0.13810E+00

x-coord 0.2113

y-coord -0.9472

sig_x

-0.37855E+00

sig_y -0.10518E+01

4

0.2113

-0.8028

-0.13687E+00

-0.94821E+00

-0.13810E+00

4

0.7887

-0.8028

-0.13687E+00

-0.94821E+00

0.13810E+00

4

0.7887

-0.9472

-0.37855E+00

-0.10518E+01

0.13810E+00


Local number

4

1

4

2

4 4

Global number 9

Int. Hor. Force

Int. Vert. Force

0.80802E-01

0.50000E+00

7

-0.16373E-01

-0.50000E+00

3

8

0.16373E-01

4

10

-0.50000E+00

-0.80802E-01 52

0.50000E+00

RESULTS: PART 1C 4 ROD ELEMENTS IN A COLUMN UNDER AXIAL COMPRESSION ON TOP

53

4 ROD ELEMENTS IN A COLUMN UNDER AXIAL COMPRESSION ON TOP: There are

4 equations and the skyline storage is

7


Element No. 1 integration point 1 shape function derivatives with respect to local coordinates: -0.5 0.5 jacobian matrix: 0.125 determinant of the jacobian matrix: 0.125 shape function derivatives with respect to global coordinates: -4. 4. strain displacement matrix, B, ip: 1 -4. 4.

cumulative km matrix, km 2000000. -2000000. -2000000. 2000000 integration point 2 shape function derivatives with respect to local coordinates: -0.5 0.5 jacobian matrix: 0.125 determinant of the jacobian matrix: 0.125 shape function derivatives with respect to global coordinates: -4. 4. strain displacement matrix, B, ip: 2 -4. 4.

54

cumulative km matrix, km 4000000. -4000000. -4000000. 4000000. Cumulative global stiffness matrix in skyline storage: 4000000. -4000000. 8000000. -4000000. 8000000. -4000000. 8000000. kdiag: [ 1 3 5 7 ]T

Results: Nodal Displacements: Node

displacement

1

-0.10000E-05

2

-0.75000E-06

3

-0.50000E-06

4

-0.25000E-06

5

0.00000E+00

Stresses at Gaussian points and Internal Forces:


Element

coord

sig

1

-0.1972

-1.0000

1

-0.0528

-1.0000

The node (nod= 2 ) internal forces are: Element 1 1

Local number 1 2

Global number

Int. Force

2

0.10000E+01

1

-0.10000E+01

55

Element no. 2 The integration point(nip= 2 ) stresses are: Element

coord

sig

2

-0.4472

-1.0000

2

-0.3028

-1.0000

The node (nod= 2 ) internal forces are: Element 2

Local number 1

Global number 3

Int. Force

0.10000E+01


Element

coord

sig

3

-0.6972

-1.0000

3

-0.5528

-1.0000


Local number

3

1

3

2

Global number 4 3

Int. Force

0.10000E+01 -0.10000E+01


Element

coord

sig

4

-0.9472

-1.0000

4

-0.8028

-1.0000


Local number

4 4

Global number

1 2

5 4

Int. Force

0.10000E+01 -0.10000E+01

56

RESULTS: PART 1 D 4 BEAM ELEMENTS IN A CANTILEVER BEAM UNDER POINT LOAD AT THE EXTREME

57

4 BEAM ELEMENTS IN A CANTILEVER BEAM UNDER POINT LOAD AT THE EXTREME: There are

8 equations and the skyline storage is 24


Cumulative global stiffness matrix in skyline storage: 2.000639928184101E+8 0. 1.6671999551150631E+7 -1.0003199640920505E+8 -2.5007999102301262E+7 2.000639928184101E+8 2.5007999102301262E+7 4167999.775575315 0. 1.6671999551150631E+7 -1.0003199640920505E+8 -2.5007999102301262E+7 2.000639928184101E+8 2.5007999102301262E+7 4167999.775575315 0. 1.6671999551150631E+7 -1.0003199640920505E+8 -2.5007999102301262E+7 1.0003199640920505E+8 2.5007999102301262E+7 4167999.775575315 -2.5007999102301262E+7 8335999.775575316 kdiag: [ 1 3 6 10 13 17 20 24 ]T Results: Nodal Displacements and rotation: Node 1

displacement 0.00000E+00

theta

0.00000E+00

2

-0.21993E-04

-0.83973E-04

3

-0.79974E-04

-0.14395E-03

4

-0.16195E-03

-0.17994E-03

5

-0.25592E-03

-0.19194E-03

58



sig=(sig/y)*h where h is the distance to the neutral axis Element

coord

sig/y

1

0.1057

-36359.6463

1

0.3943

-30818.8566

The node (nod= 2 ) internal forces are:

Element Moment 1 1

Local number 1 2

Global number

Int. Shear. For

Int. Bending

1

0.10000E+03

0.20000E+03

2

-0.10000E+03

-0.15000E+03



coord

sig/y

2

0.6057

-26762.7173

2

0.8943

-21221.9276

The node (nod= 2 ) internal forces are: Element Moment 2 2

Local number 1 2

Global number 2 3

Int. Shear. For

0.10000E+03 -0.10000E+03

59

Int. Bending 0.15000E+03

-0.10000E+03



coord

sig/y

3

1.1057

-17165.7883

3

1.3943

-11624.9987

The node (nod= 2 ) internal forces are: Element Moment 3

Local number 1

3

3

2

Global number

Int. Shear. For

0.10000E+03 4

Int. Bending

0.10000E+03

-0.10000E+03

-0.50000E+02



coord

sig/y

4

1.6057

-7568.8593

4

1.8943

-2028.0697

The node (nod= 2 ) internal forces are: Element Moment

Local number

Global number

Int. Shear. For

Int. Bending

4

1

4

0.10000E+03

0.50000E+02

4

2

5

-0.10000E+03

0.69289E-12

60

RESULTS: PART 2A 4 ELEMENT COLUMN (2D CUADRATIC QUADRILATERAL ELEMENT):

61

4 ELEMENT COLUMN (2D CUADRATIC QUADRILATERAL ELEMENT): There are 40 equations and the skyline storage is 400 Results: Node

x-disp

y-disp

1

-0.43147E-06

-0.12326E-05

2

-0.13329E-19

-0.71118E-06

3

0.43147E-06

-0.12326E-05

4

-0.25232E-06

-0.10062E-05

5

0.25232E-06

-0.10062E-05

6

-0.18228E-06

-0.79311E-06

7

-0.81795E-20

-0.56543E-06

8

0.18228E-06

-0.79311E-06

9

-0.16068E-06

-0.63654E-06

10

0.16068E-06

-0.63654E-06

11

-0.15556E-06

-0.48141E-06

12

-0.35908E-20

-0.39018E-06

13

0.15556E-06

-0.48141E-06

14

-0.15100E-06

-0.35968E-06

15

0.15100E-06

-0.35968E-06

16

-0.13143E-06

-0.23925E-06

17

-0.47335E-21

-0.18495E-06

18

0.13143E-06

-0.23925E-06

19

-0.88493E-07

-0.11610E-06

20

0.88493E-07

-0.11610E-06

21

0.00000E+00

0.00000E+00

22

0.00000E+00

0.00000E+00

23

0.00000E+00

0.00000E+00

62



Element

x-coord

y-coord

sig_x

sig_y

tau_xy

1

0.1127

-0.2218

-0.18778E+00

-0.14051E+01

1

0.1127

-0.1250

-0.63812E-01

1

0.1127

-0.0282

0.23643E+00

-0.14051E+01

0.83570E-01

1

0.5000

-0.0282

0.64320E+00

-0.45591E+00

-0.33721E-15

1

0.8873

-0.0282

0.23643E+00

-0.14051E+01

-0.83570E-01

1

0.8873

-0.1250

-0.63812E-01

1

0.8873

-0.2218

-0.18778E+00

-0.14051E+01

-0.21529E+00

1

0.5000

-0.2218

0.21899E+00

-0.45591E+00

-0.28630E-15

1

0.5000

-0.1250

0.34296E+00

-0.49369E+00

-0.22755E-15

-0.14428E+01

-0.14428E+01

0.21529E+00 0.14943E+00

-0.14943E+00


Local number

Global number

Int. Hor. Force

1

1

6

1

2

4

0.15951E-16

0.32324E-14

1

0.20219E-15

-0.25000E+00

1

3

-0.41485E-01

Int. Vert. Force

0.21785E+00

1

4

2

-0.34407E-15

1

5

3

0.33732E-16

-0.25000E+00

1

6

5

-0.61149E-16

0.16395E-14

1

7

8

0.41485E-01

0.21785E+00

1

8

7

0.27368E-15

0.56430E+00

63

-0.50000E+00


Element

x-coord

y-coord

sig_x

sig_y

tau_xy

2

0.1127

-0.4718

-0.17019E+00

-0.11934E+01

0.12562E+00

2

0.1127

-0.3750

-0.16076E+00

-0.11991E+01

0.15818E+00

2

0.1127

-0.2782

-0.12470E+00

-0.11934E+01

0.19074E+00

2

0.5000

-0.2782

0.64242E-01

2

0.8873

-0.2782

-0.12470E+00

-0.11934E+01

-0.19074E+00

2

0.8873

-0.3750

-0.16076E+00

-0.11991E+01

-0.15818E+00

2

0.8873

-0.4718

-0.17019E+00

-0.11934E+01

-0.12562E+00

2

0.5000

-0.4718

0.18751E-01

-0.75254E+00

-0.65351E-15

2

0.5000

-0.3750

0.28186E-01

-0.75824E+00

-0.49957E-15

-0.75254E+00

-0.51860E-15


Local number

Global number

Int. Hor. Force

2

1

11

2

2

9

2

3

6

2

4

2

5

8

2

6

10

-0.46024E-16

-0.26372E-15

2

7

13

0.24138E-01

0.18381E+00

2

8

12

0.37147E-15

0.63237E+00

7

-0.24138E-01

Int. Vert. Force

0.81857E-17

0.18381E+00

-0.22829E-15

0.41485E-01

-0.21785E+00

-0.32082E-15

-0.56430E+00

-0.41485E-01

-0.21785E+00

64


Element

x-coord

y-coord

sig_x

sig_y

tau_xy

3

0.1127

-0.7218

-0.15086E+00

-0.10553E+01

0.13175E-01

3

0.1127

-0.6250

-0.11818E+00

-0.10501E+01

0.57949E-01

3

0.1127

-0.5282

-0.10977E+00

-0.10553E+01

0.10272E+00

3

0.5000

-0.5282

-0.58632E-01

3

0.8873

-0.5282

-0.10977E+00

-0.10553E+01

-0.10272E+00

3

0.8873

-0.6250

-0.11818E+00

-0.10501E+01

-0.57949E-01

3

0.8873

-0.7218

-0.15086E+00

-0.10553E+01

-0.13175E-01

3

0.5000

-0.7218

-0.99717E-01

-0.93602E+00

-0.54088E-15

3

0.5000

-0.6250

-0.67044E-01

-0.93082E+00

-0.50613E-15

-0.93602E+00

-0.48873E-15

The node (nod= 8 ) internal forces are: Element 3 3 3

Local number 1 2 3

Global number

Int. Hor. Force

Int. Vert. Force

16

0.14107E-02

0.17134E+00

14

-0.31191E-16

0.16881E-15

11

0.24138E-01

-0.18381E+00

3

4

12

-0.29293E-15

-0.63237E+00

3

5

13

-0.24138E-01

-0.18381E+00

3

6

15

0.38842E-16

0.50807E-15

3

7

18

-0.14107E-02

3

8

17

0.36607E-15

65

0.17134E+00 0.65731E+00


Element

x-coord

y-coord

sig_x

sig_y

tau_xy

4

0.1127

-0.9718

-0.41241E+00

-0.10850E+01

-0.23344E+00

4

0.1127

-0.8750

-0.26374E+00

-0.10692E+01

-0.12427E+00

4

0.1127

-0.7782

-0.18866E+00

-0.10850E+01

-0.15105E-01

4

0.5000

-0.7782

-0.11349E+00

-0.90956E+00

0.70188E-15

4

0.8873

-0.7782

-0.18866E+00

-0.10850E+01

0.15105E-01

4

0.8873

-0.8750

-0.26374E+00

-0.10692E+01

0.12427E+00

4

0.8873

-0.9718

-0.41241E+00

-0.10850E+01

0.23344E+00

4

0.5000

-0.9718

-0.33723E+00

-0.90956E+00

0.68800E-15

4

0.5000

-0.8750

-0.18856E+00

-0.89379E+00

0.63171E-15


Local number

Global number

Int. Hor. Force

Int. Vert. Force

4

1

21

0.64103E-01

0.19808E+00

4

2

19

-0.35383E-16

-0.21066E-15

4

3

16

-0.14107E-02

-0.17134E+00

4

4

17

-0.32865E-15

-0.65731E+00

4

5

18

0.14107E-02

-0.17134E+00

4

6

20

-0.31369E-16

-0.39693E-15

4

7

23

-0.64103E-01

0.19808E+00

4

8

22

-0.88566E-15

0.60383E+00

66

RESULTS: PART 2B 4 ELEMENT COLUMN (1D CUADRATIC ROD ELEMENTS):

67

4 ELEMENT COLUMN (1D CUADRATIC ROD ELEMENTS): There are

8 equations and the skyline storage is 18

Results: Nodal Displacements: Node

displacement

1

-0.10000E-05

2

-0.87500E-06

3

-0.75000E-06

4

-0.62500E-06

5

-0.50000E-06

6

-0.37500E-06

7

-0.25000E-06

8

-0.12500E-06

9

0.00000E+00



Element

coord

sig

1

-0.1972

-1.0000

1

-0.0528

-1.0000


Local number

Global number

Int. Force

1

1

3

0.10000E+01

1

3

1

-0.10000E+01

68


Element

coord

sig

2

-0.4472

-1.0000

2

-0.3028

-1.0000


Local number

Global number

Int. Force

2

1

5

0.10000E+01

2

3

3

-0.10000E+01


Element 3 3

coord

sig

-0.6972 -0.5528

-1.0000 -1.0000


Local number

Global number

Int. Force

3

1

7

0.10000E+01

3

3

5

-0.10000E+01

69


Element

coord

sig

4

-0.9472

-1.0000

4

-0.8028

-1.0000


Local number

Global number

4

1

9

4

2

8

4

3

7

Int. Force

0.10000E+01 0.19971E-15 -0.10000E+01

70

RESULTS: PART 3A IRREGULAR ELEMENT COLUMN (2D BILINEAR QUADRILATERAL ELEMENT)

71

IRREGULAR ELEMENT COLUMN (2D BILINEAR QUADRILATERAL ELEMENT) There are 288 equations and the skyline storage is 5804

Results: Node

x-disp

y-disp

1 -0.20202E-06 -0.89280E-06 2 -0.15311E-06 -0.88971E-06 3 -0.10332E-06 -0.88601E-06 4 -0.52138E-07 -0.88292E-06 5 -0.54613E-20 -0.88173E-06 6 0.52138E-07 -0.88292E-06 7 0.10332E-06 -0.88601E-06 8 0.15311E-06 -0.88971E-06 9 0.20202E-06 -0.89280E-06 10 -0.20061E-06 -0.83593E-06 11 -0.15164E-06 -0.83276E-06 12 -0.10193E-06 -0.82879E-06 13 -0.51273E-07 -0.82545E-06 14 -0.49025E-20 -0.82416E-06 15 0.51273E-07 -0.82545E-06 72

Stresses at Gaussian points and Internal Forces: Element no. 1 The integration point(nip= 4 ) stresses are: Element

x-coord

y-coord

sig_x

sig_y

tau_xy

1 0.91682E-03

0.0264

-0.0493

0.20531E-02

-0.99939E+00

1 0.80037E-03

0.0264

-0.0132

0.17584E-02

-0.99952E+00

1 0.63197E-03

0.0986

-0.0132

0.14091E-02

-0.10003E+01

1 0.74842E-03

0.0986

-0.0493

0.17038E-02

-0.10002E+01


Local number

Global number

Int. Hor. Force

1

1

10

-0.10819E-03

1

2

1

0.73360E-16

Int. Vert. Force

0.62452E-01 -0.62500E-01


Element

x-coord

y-coord

sig_x

sig_y

tau_xy

2 0.31352E-02

0.1514

-0.0493

0.93670E-02

-0.99820E+00

2 0.26421E-02

0.1514

-0.0132

0.98193E-02

-0.99801E+00

2 0.29006E-02

0.2236

-0.0132

0.83401E-02

-0.10015E+01

2 0.33937E-02

0.2236

-0.0493

0.78878E-02

-0.10017E+01


Local number

Global number

2

1

11

2

2

2

Int. Hor. Force

-0.45655E-03 -0.96799E-04 73

Int. Vert. Force

0.62328E-01 -0.62517E-01

RESULTS: PART 3B IRREGULAR ELEMENTS CANTILEVER BEAM (2D BILINEAR QUADRILATERAL ELEMENTS)

74

IRREGULAR ELEMENT CANTILEVER BEAM (2D BILINEAR QUADRILATERAL ELEMENTS)

28 0.00000E+00 0.00000E+00

There are 144 equations and the skyline storage is 2516

Results: Node

x-disp

y-disp

1 0.00000E+00 0.00000E+00

29 02

0.40882E-03

-0.10414E-

30 02

0.83126E-03

-0.39415E-

31 02

0.11759E-02

-0.83230E-

32 01

0.14563E-02

-0.13956E-

2 02

0.18049E-02

-0.14553E-

33 01

0.16756E-02

-0.20590E-

3 02

0.34151E-02

-0.42382E-

34 01

0.18323E-02

-0.27975E-

4 02

0.47956E-02

-0.85743E-

35 01

0.19231E-02

-0.35864E-

5 01

0.59214E-02

-0.14157E-

36 01

0.19575E-02

-0.43991E-

6 01

0.67972E-02

-0.20741E-

37 0.00000E+00 0.00000E+00

7 01

0.74220E-02

-0.28076E-

38 02

-0.16906E-16

-0.10161E-

8 01

0.77994E-02

-0.35906E-

39 02

-0.40443E-16

-0.39204E-

9 01

0.79387E-02

-0.44031E-

40 02

-0.63983E-16

-0.83065E-

41 01

10 0.00000E+00 0.00000E+00

-0.86730E-16

-0.13943E-

11 02

0.12839E-02

-0.12539E-

42 01

-0.10949E-15

-0.20580E-

12 02

0.25373E-02

-0.41051E-

43 01

-0.13115E-15

-0.27968E-

13 02

0.35620E-02

-0.84553E-

44 01

-0.13608E-15

-0.35861E-

14 01

0.44080E-02

-0.14063E-

45 01

-0.13125E-15

-0.43986E-

46

75

0.00000E+00

57 02

-0.16742E-02

-0.40040E-

70 01

-0.55342E-02

-0.28029E-

58 02

-0.23603E-02

-0.83725E-

71 01

-0.58132E-02

-0.35885E-

59 01

-0.29224E-02

-0.13996E-

72 01

-0.59156E-02

-0.44018E-

60 01

-0.33606E-02

-0.20621E-

73 0.00000E+00 0.00000E+00

61 01

-0.36740E-02

-0.27995E-

74 02

-0.18049E-02

-0.14553E-

62 01

-0.38574E-02

-0.35871E-

75 02

-0.34151E-02

-0.42382E-

63 01

-0.39258E-02

-0.44002E-

76 02

-0.47956E-02

-0.85743E-

77 01

-0.59214E-02

-0.14157E-

64 0.00000E+00 0.00000E+00



Element

x-coord

y-coord

sig_x

sig_y

tau_xy

1

0.0528

-0.0493

0.71132E+04

0.23004E+04

-0.13170E+04

1

0.0528

-0.0132

0.87331E+04

0.29946E+04

-0.14958E+04

1

0.1972

-0.0132

0.76600E+04

0.49066E+03

0.35545E+03

1

0.1972

-0.0493

0.60401E+04

-0.20357E+03

0.53430E+03


Local number

Global number

Int. Hor. Force

1

1

10

-0.89320E+02

1

2

1

-0.37234E+03

1

3

2

0.25215E+03

1

4

11

0.20951E+03

76

Int. Vert. Force -0.25138E+03

0.28143E+03 0.67448E+02 -0.97496E+02


Element

x-coord

y-coord

sig_x

sig_y

tau_xy

2

0.3028

-0.0493

0.54296E+04

-0.95971E+03

-0.69382E+03

2

0.3028

-0.0132

0.65388E+04

-0.48437E+03

-0.63318E+03

2

0.4472

-0.0132

0.69026E+04

0.36457E+03

0.63438E+03

2

0.4472

-0.0493

0.57935E+04

-0.11077E+03

0.57375E+03


Local number

Global number

2

1

11

2

2

2

2

3

3

2

4

12

Int. Hor. Force

-0.13323E+03 -0.25215E+03

Int. Vert. Force 0.69306E+02

-0.67448E+02

0.24472E+03

-0.69445E+01

0.14066E+03

0.50871E+01


Element

x-coord

y-coord

sig_x

sig_y

tau_xy

63

1.5528

-0.4868

-0.15342E+04

0.77391E+02

-0.31015E+03

63

1.5528

-0.4507

-0.12281E+04

0.20856E+03

-0.33325E+03

63

1.6972

-0.4507

-0.13667E+04

-0.11486E+03

0.16546E+02

63

1.6972

-0.4868

-0.16728E+04

-0.24603E+03

0.39647E+02


Local number

Global number

Int. Hor. Force

Int. Vert. Force

63

1

79

0.79061E+02

-0.49496E+01

63

2

70

0.11594E+02

0.14125E+02

63

3

71

-0.48294E+02

-0.18808E+02

63

4

80

-0.42360E+02

0.96331E+01

77


Element

x-coord

y-coord

sig_x

sig_y

tau_xy

64

1.8028

-0.4868

-0.52808E+03

0.11652E+03

-0.23125E+03

64

1.8028

-0.4507

-0.41335E+03

0.16569E+03

-0.23876E+03

64

1.9472

-0.4507

-0.45841E+03

0.60548E+02

-0.10764E+03

64

1.9472

-0.4868

-0.57314E+03

0.11378E+02

-0.10013E+03


Local number

Global number

Int. Hor. Force

Int. Vert. Force

64

1

80

0.42360E+02

64

2

71

-0.11533E+02

0.20223E+02

64

3

72

-0.30828E+02

0.19099E+01

64

4

81

-0.19327E-11

78

-0.96331E+01

-0.12500E+02

FORTRAN CODE: “LINEAR ELEMENTS AND QUADRATIC ELEMENTS”: PART 1: -LINEAR AND QUADRATIC TWO- AND ONE- DIMENSIONAL ELEMENTS GENERAL FEM MANUAL NODE AND ELEMENT NUMBERING PART 2: -LINEAR AND QUADRATIC TWO- AND ONE- DIMENSIONAL ELEMENTS GENERAL FEM AUTOMATE NODE AND ELEMENT NUMBERING

79

FORTRAN CODE: PART 1: -LINEAR AND QUADRATIC TWO- AND ONE- DIMENSIONAL ELEMENTS GENERAL FEM, MANUAL NODE AND ELEMENT NUMBERING

80

FORTRAN CODE: PART 2: -LINEAR AND QUADRATIC TWO- AND ONE- DIMENSIONAL

ELEMENTS GENERAL FEM, AUTOMATE NODE AND ELEMENT NUMBERING

81

Numerical Modelling in Geo engineering

Recommend Documents