(c) Juan Gomez 2019. Thanks to Universidad EAFIT for support. This material is part of the course Introduction to Finite Element Analysis

# Two-dimensional truss elements

The simple code implemented for the springs elements in a previous notebook can be easily adapted to consider an assemblage of two-dimensional truss elements. In fact the code can be converted into a general program to consider a variety of (finite) elements just by adding subroutines with the force-displacement (stiffness) equations for each element.

In this particular problem a typical element of axial stiffness $\frac{AE}{\mathcal l}$ is shown in its local reference system in the figure below: 

<center></center>

Note that since the element has only axial stiffness the force-displacement relationship in its local reference system is of the form: 

$$
\begin{Bmatrix}f_1\\f_2\end{Bmatrix}=\frac{AE}{\mathcal l}\begin{bmatrix}1&-1\\-1&1\end{bmatrix}\begin{Bmatrix}u_1\\u_2\end{Bmatrix}
$$

which is also analogous to the one in the spring element.


### Assemblage of truss elements.

Cosnsider the two-elements assemblage shown below.


<center></center>

To obtain the total stiffness from the structure it is now necessary to consider the stiffness contribution from all the elements in a common (**Global**) reference system.

Let:

* $U, F$ : Displacements degrees of freedom and forces in the global reference system.
* $u, f$ : Displacements degrees of freedom and forces in the local reference system.

and related by the rotational transformation matrix $\lambda$ like;

$$u=\lambda U.$$

Using the fact that the virtual energy of the forces upon imposition of virtual displacements is a reference-independent scalar quantity gives:

$$\delta U^TF=\delta u^Tf.$$

Using the first equation in the second yields;

$$\delta U^TF=\delta U^T\lambda^Tf$$

from which:

$$F=\lambda^Tf.$$

Now, conisdering the equilibrium relation for the element in the local system:

$$f=ku$$

where $k$ is the local stiffness matrix we can write:

$$\begin{array}{l}\lambda^Tf=\lambda^Tku\\\lambda^Tf=\lambda^Tk\lambda U\\F=KU\end{array}$$

from which:

$$K=\lambda^Tk\lambda$$

where $K$ is the stiffness matrix for the two-dimensional truss element in the global reference system. It must be observed that in the global reference system the element has two degrees of freedom per node.

In the actual implementation all the information required to compute $K$ is passed as input paramters to the elemental subroutine **uel()** as described below.

**Question: Find the rotational transformation matrix $\lambda$ required for the formulation of the stiffness matrix in the global reference system.**

### Example

The modifications that must be applied to the spring-elements based code are only related to the fact that now there are 2 degrees of freedom are each nodal point.

(The required input files containing the input data for the nodes, elements, materials parameters and loads are in the folder **files** of this REPO.)


```python
%matplotlib inline        
import matplotlib.pyplot as plt
import numpy as np
import sympy as sym
```

Read the input files from the **files** folder. 


```python
def readin():
    nodes    = np.loadtxt('files/' + 'Cnodes.txt', ndmin=2)
    mats     = np.loadtxt('files/' + 'Cmater.txt', ndmin=2)
    elements = np.loadtxt('files/' + 'Celes.txt' , ndmin=2)
    loads    = np.loadtxt('files/' + 'Cloads.txt', ndmin=2)

    return nodes, mats, elements, loads
```

**eqcounter** counts equations and generates the boundary conditions array in its second instance.


```python
def eqcounter(nodes):
    nnodes = nodes.shape[0]
    IBC = np.zeros([nnodes, 2], dtype=np.integer)
    for i in range(nnodes):
        for k in range(2):
            IBC[i , k] = int(nodes[i , k+3])
    neq = 0
    for i in range(nnodes):
        for j in range(2):
            if IBC[i, j] == 0:
                IBC[i, j] = neq
                neq = neq + 1

    return neq, IBC
```

**DME** computes the assembly operator.


```python
def DME(nodes, elements):

    nels = elements.shape[0]
    IELCON = np.zeros([nels, 2], dtype=np.integer)
    DME    = np.zeros([nels, 4], dtype=np.integer)
    neq , IBC = eqcounter(nodes)
    ndof   = 4
    nnodes = 2
    for i in range(nels):
        for j in range(nnodes):
            IELCON[i, j] = elements[i, j+3]
            kk = IELCON[i, j]
            for l in range(2):
                DME[i, 2*j+l] = IBC[kk, l]
    return DME , IBC , neq
```

**assembly** uses the model and the **DME** operator to compute the global stiffness matrix.



```python
def assembly(elements, mats, nodes, neq, DME, uel=None):

    IELCON = np.zeros([2], dtype=np.integer)
    KG = np.zeros((neq, neq))
    nels = elements.shape[0]
    nnodes = 2
    ndof = 4
    for el in range(nels):
        elcoor = np.zeros([nnodes, 2])
        im     = np.int(elements[el , 2])
        par0= mats[im ,0]
        par1= mats[im , 1]
        for j in range(nnodes):
            IELCON[j] = elements[el , j+3]
            elcoor[j , 0] = nodes[IELCON[j], 1]
            elcoor[j , 1] = nodes[IELCON[j], 2]
        kloc = ueltruss2D(elcoor, par0, par1)
        dme = DME[el , :ndof]
        for row in range(ndof):
            glob_row = dme[row]
            if glob_row != -1:
                for col in range(ndof):
                    glob_col = dme[col]
                    if glob_col != -1:
                        KG[glob_row, glob_col] = KG[glob_row, glob_col] +\
                                                 kloc[row, col]

    return KG
```

**ueltruss2D** uses the nodal point coordinates and the material parameters to compute the local stiffness matrix transformed to the global reference system.

**Question: Add comments to explains the different steps in the following subroutine. In particular identify the computation of the rotational transformation matrix $\lambda$.**


```python
def ueltruss2D(coord, A, Emod):
    vec = coord[1, :] - coord[0, :]
    length = np.linalg.norm(vec)
    nx = vec[0]/length
    ny = vec[1]/length
    Q = np.array([
        [nx, ny , 0 , 0],
        [0,  0, nx , ny]])
    kl = (A*Emod/length) * np.array([
        [1, -1],
        [-1, 1]])
    kG = np.dot(np.dot(Q.T, kl), Q)
    return kG
```

**loadassem** forms the vector of nodal loads.


```python
def loadasem(loads, IBC, neq, nl):
    """Assembles the global Right Hand Side Vector RHSG

    Parameters
    ----------
    loads : ndarray
      Array with the loads imposed in the system.
    IBC : ndarray (int)
      Array that maps the nodes with number of equations.
    neq : int
      Number of equations in the system after removing the nodes
      with imposed displacements.
    nl : int
      Number of loads.

    Returns
    -------
    RHSG : ndarray
      Array with the right hand side vector.

    """
    RHSG = np.zeros([neq])
    for i in range(nl):
        il = int(loads[i, 0])
        ilx = IBC[il , 0]
        ily = IBC[il , 1]
        if ilx != -1:
            RHSG[ilx] = loads[i, 1]
        if ily != -1:
            RHSG[ily] = loads[i, 2]

    return RHSG
```

The main program still retains the same structure as follows:

* Reads the model
* Builds the DME() operator
* Assembles the global system of equations
* Solves for the global displacements $UG$


```python
nodes, mats, elements, loads = readin()
DME , IBC , neq = DME(nodes, elements)
KG   = assembly(elements, mats, nodes, neq, DME)
RHSG = loadasem(loads, IBC, neq, 1)
UG = np.linalg.solve(KG, RHSG)
print(UG)
```

    [ 0.53009278 -0.17789364]


### Proposed problems
#### Problem 1
Implement a subroutine to compute the nodal forces in each element and verify the equilibrium of the system.

#### Problem 2
Find the lateral stiffness of the structure using the relation:

$$k=\frac P\delta$$

#### Problem 3
By introducing the necessary modifications to the input files, find the resulting displacement solution after adding a third bar connecting nodes $0$ and $2$. Comment on your results.

#### Problem 4
Fix the truss structure shown in the figure by adding elements and/or imposing appropiate displacement boundary conditions. (Create a new set of input files.)

<center></center>

### References

* Bathe, Klaus-Jürgen. (2006) Finite element procedures. Klaus-Jurgen Bathe. Prentice Hall International.

* Juan Gómez, Nicolás Guarín-Zapata (2018). SolidsPy: 2D-Finite Element Analysis with Python, <https://github.com/AppliedMechanics-EAFIT/SolidsPy>.


```python
from IPython.core.display import HTML
def css_styling():
    styles = open('./nb_style.css', 'r').read()
    return HTML(styles)
css_styling()
```





<link href='http://fonts.googleapis.com/css?family=Fenix' rel='stylesheet' type='text/css'>
<link href='http://fonts.googleapis.com/css?family=Alegreya+Sans:100,300,400,500,700,800,900,100italic,300italic,400italic,500italic,700italic,800italic,900italic' rel='stylesheet' type='text/css'>
<link href='http://fonts.googleapis.com/css?family=Source+Code+Pro:300,400' rel='stylesheet' type='text/css'>

<style>

/*
Template for Notebooks for Modelación computacional.

Based on Lorena Barba template available at:

    https://github.com/barbagroup/AeroPython/blob/master/styles/custom.css
*/

/* Fonts */
@font-face {
font-family: "Computer Modern";
src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunss.otf');
}

/* Text */
div.cell{
width:800px;
margin-left:16% !important;
margin-right:auto;
}
h1 {
font-family: 'Alegreya Sans', sans-serif;
}
h2 {
font-family: 'Fenix', serif;
}
h3{
font-family: 'Fenix', serif;
margin-top:12px;
margin-bottom: 3px;
}
h4{
font-family: 'Fenix', serif;
}
h5 {
font-family: 'Alegreya Sans', sans-serif;
}	
div.text_cell_render{
font-family: 'Alegreya Sans',Computer Modern, "Helvetica Neue", Arial, Helvetica, Geneva, sans-serif;
line-height: 135%;
font-size: 120%;
width:600px;
margin-left:auto;
margin-right:auto;
}
.CodeMirror{
font-family: "Source Code Pro";
font-size: 90%;
}
/* .prompt{
display: None;
}*/
.text_cell_render h1 {
font-weight: 200;
font-size: 50pt;
line-height: 100%;
color:#CD2305;
margin-bottom: 0.5em;
margin-top: 0.5em;
display: block;
}	
.text_cell_render h5 {
font-weight: 300;
font-size: 16pt;
color: #CD2305;
font-style: italic;
margin-bottom: .5em;
margin-top: 0.5em;
display: block;
}
.warning{
color: rgb( 240, 20, 20 )
}
</style>









```python

```

d6237c8fba83bc5a6a5c7d3dc2173e6777617486
ipynb
Jupyter Notebook
notebooks/13_truss_elements.ipynb
AppliedMechanics-EAFIT/Introductory-Finite-Elements
a4b44d8bf29bcd40185e51ee036f38102f9c6a72
MIT
2019-11-26T13:28:30.000Z
2022-02-16T17:57:11.000Z
notebooks/13_truss_elements.ipynb
jgomezc1/Introductory-Finite-Elements
a4b44d8bf29bcd40185e51ee036f38102f9c6a72
MIT
notebooks/13_truss_elements.ipynb
jgomezc1/Introductory-Finite-Elements
a4b44d8bf29bcd40185e51ee036f38102f9c6a72
MIT
2020-02-17T07:24:59.000Z
2022-03-02T07:54:28.000Z
Qwen/Qwen-72B
1. YES
2. YES