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

# Two-dimensional frame elements

The code can now be easily extended to consider also frame elements as the one shown in the figure in its local reference system: 

<center></center>

In this case the element has infinite axial stiffness therefore in the local system the only degrees of freedom are the transverse displacement and the rotation. The generalized displacements vector reads:

$$
u^T=\begin{bmatrix}v_1&\theta_1&v_2&\theta_2\end{bmatrix}
$$

while the vector of generalized forces (shears and moments) is:

$$
f^T=\begin{bmatrix}f_1&m_1&f_2&m_2\end{bmatrix}
$$


The force-displacement stiffness matrix in the local reference system is:

$$
$$
\begin{bmatrix}12\frac{EI}{\mathcal l^3}&6\frac{EI}{\mathcal l^3}&-12\frac{EI}{\mathcal l^3}&6\frac{EI}{\mathcal l^2}\\6\frac{EI}{\mathcal l^3}&4\frac{EI}{\mathcal l}&-6\frac{EI}{\mathcal l^2}&2\frac{EI}{\mathcal l}\\12\frac{EI}{\mathcal l^3}&-6\frac{EI}{\mathcal l^2}&12\frac{EI}{\mathcal l^3}&-6\frac{EI}{\mathcal l^2}\\6\frac{EI}{\mathcal l^2}&2\frac{EI}{\mathcal l}&-6\frac{EI}{\mathcal l^2}&4\frac{EI}{\mathcal l}\end{bmatrix}
$$
$$

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 and for that purpose we proceed like in the two-dimensional truss problem giving for the elemental stiffness matrix in the global reference system:

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

where $K$ is the stiffness matrix for the two-dimensional frame element in the global reference system. It must be observed that in the global reference system the element has now three 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.

### Simple framed strucure

Cosnsider the following three-elements assemblage (The required input files are available in the files folder from the REPO):


<center></center>

find the lateral displacement of the structure when subject to the lateral load $P.$

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

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


```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/' + 'Fnodes.txt', ndmin=2)
    mats     = np.loadtxt('files/' + 'Fmater.txt', ndmin=2)
    elements = np.loadtxt('files/' + 'Feles.txt' , ndmin=2)
    loads    = np.loadtxt('files/' + 'Floads.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, 3], dtype=np.integer)
    for i in range(nnodes):
        for k in range(3):
            IBC[i, k] = int(nodes[i, k+3])
    neq = 0
    for i in range(nnodes):
        for j in range(3):
            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, 6], dtype=np.integer)
    neq, IBC = eqcounter(nodes)
    ndof = 6
    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(3):
                DME[i, 3*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 = 6
    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 = uelbeam2DU(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
```

**uelbeam2D** 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 uelbeam2DU(coord, I, Emod):
    """2D-2-noded beam element
       without axial deformation

    Parameters
    ----------
    coord : ndarray
      Coordinates for the nodes of the element (2, 2).
    A : float
      Cross section area.
    Emod : float
      Young modulus (>0).

    Returns
    -------
    kl : ndarray
      Local stiffness matrix for the element (4, 4).

    """
    vec = coord[1, :] - coord[0, :]
    nx = vec[0]/np.linalg.norm(vec)
    ny = vec[1]/np.linalg.norm(vec)
    L = np.linalg.norm(vec)
    Q = np.array([
        [-ny, nx,   0,  0,  0, 0],
        [0,  0, 1.0,  0,  0, 0],
        [0,  0,   0, -ny, nx, 0],
        [0,  0,  0,  0,  0, 1.0]])
    kl = (I*Emod/(L*L*L)) * np.array([
        [12.0, 6, -12.0, 6*L],
        [6,  4*L*L, -6*L, 2*L*L],
        [-12.0,  -6*L, 12.0, -6*L],
        [6*L,  2*L*L, -6*L, 4*L*L]])
    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]
        ilT = IBC[il, 2]
        if ilx != -1:
            RHSG[ilx] = loads[i, 1]
        if ily != -1:
            RHSG[ily] = loads[i, 2]
        if ilT != -1:
            RHSG[ilT] = loads[i, 3]

    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)
```

    [  2.25000000e+01   2.00000000e+01   1.27557539e-16  -1.25918779e-15
       2.00000000e+01   4.88970566e-16]


### 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
Retrofit the structure in such a way that the lateral stiffness increases by a factor of 2.0.

#### Problem 4
Fix the framed structure shown in the figure by adding elements and/or imposing appropiate displacement boundary conditions. (Note: you must 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

```

5ad7e0c3649669e070ce40c5d51c1d7e6b5370e7
ipynb
Jupyter Notebook
notebooks/14_frame_elements.ipynb
AppliedMechanics-EAFIT/Introductory-Finite-Elements
a4b44d8bf29bcd40185e51ee036f38102f9c6a72
MIT
2019-11-26T13:28:30.000Z
2022-02-16T17:57:11.000Z
notebooks/14_frame_elements.ipynb
jgomezc1/Introductory-Finite-Elements
a4b44d8bf29bcd40185e51ee036f38102f9c6a72
MIT
notebooks/14_frame_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