<jupyter_start><jupyter_text>Finite element analysis of an elastic frame ______________________________*Date:* 2018-01-14*Author:* <NAME> ()*Course:* Gestion des Risques Structure (CentraleSupélec, option / MS "Aménagement et Construction Durables")______________________________ This document presents the principal steps in a finite element analysis program for a simple  structure.For any FE program, the overall procedure is the same.The three main steps are:- Pre-processing: the system is described by the user (inputs)- Processing: the FE mesh is generated and the problem is numerically solved- Post-processing: plotting and analysis of the outputs We first import the Python packages and modules we will need:<jupyter_code>import numpy as np
import numpy.matlib
from numpy import linalg
import mesGen
import strMod
import impGM
import modAna<jupyter_output><empty_output><jupyter_text>1. Description of the system User inputs here the description of the structure she/he wants to model:- Structure geometry- Element section geometry- Element material properties- Boundary conditions (loading)- Parameterization for the solution procedure- Output data  1.1. Structure geometry We consider a 2D frame structure.We first create an array with, for each structural element, the following list: [node I, node J, element type (1: beam/column, 2: rotational spring), element length, global to local rotation angle] (see figure with the description of the structure): Nodes coordinates (2D problem):<jupyter_code>Lf = 8      # bay length [m]
Hf = 4      # story height [m]
nodCoo = []
nodCoo.append([0,0])     # node 0 [X coordinate,Y coordinate]
nodCoo.append([Lf,0])    # node 1 [X,Y]
nodCoo.append([0,Hf])    # ...
nodCoo.append([Lf,Hf])
nodCoo.append([0,2*Hf])
nodCoo.append([Lf,2*Hf])<jupyter_output><empty_output><jupyter_text>Elements connect a node I to a node J; we consider different:- element types (0->beam/column; 1->rotational spring)- element sections<jupyter_code>strEle = []
strEle.append([0,2,0,0])  # elt 0 [node I, node J, element type, section type] 
strEle.append([2,4,0,0])  # elt 1
strEle.append([1,3,0,1])  # ...
strEle.append([3,5,0,1])
strEle.append([2,3,0,2])
strEle.append([4,5,0,3])
strEle.append([0,3,0,4])
strEle.append([1,2,0,4])<jupyter_output><empty_output><jupyter_text>1.2. Element properties (material and section geometry) The geometric properties of selected steel beam sections (see figure) can for instance be found here: ![Geometry of a steel beam section.](img/SteelBeamSection.jpg)*Remark:* We will consider that the beams can be described using Euler-Bernoulli kinematic assumptions. Other options are generally available in FE analysis programs such as Timoshenko kinematics.E [N/m2]: material Young's modulusA [m2]  : element cross-section areaIy [m4] : element cross-section moment of inertia wrt centroid<jupyter_code>eltPro = []
eltPro.append([2.1e11,33.4e-4,2772e-8])   # section type 0  [E, A, Iy] -> IPE 220
eltPro.append([2.1e11,33.4e-4,2772e-8])   # section type 1
eltPro.append([2.1e11,33.4e-4,2772e-8])   # ...
eltPro.append([2.1e11,33.4e-4,2772e-8])
eltPro.append([2.1e11,2.5e-3,0])<jupyter_output><empty_output><jupyter_text>1.4. Boundary conditions In 'fixDof' array, we indicate for every node whether the DOFs are free or not.The line number of the matrix corresponds to the node number.For each of the 3 DOFs $(u,v,\theta)$, there is a 1 if the DOF is free and a 0 otherwise:<jupyter_code>fixDof = []
fixDof.append([0,0,0])   # node 0
fixDof.append([0,0,0])   # node 1
fixDof.append([1,1,1])   # ...
fixDof.append([1,1,1])
fixDof.append([1,1,1])
fixDof.append([1,1,1])<jupyter_output><empty_output><jupyter_text>The frame is fixed in the ground at its base.  2. Mesh generation In this simple example, the only thing that has to be done for generating the mesh is computing element length and orientation ($cos(\theta)$ and $sin(\theta)$ where $\theta$ is the angle between local and global coordinate systems). We get these mesh data as:<jupyter_code># mesh data
mesDat = mesGen.eltL0(strEle,nodCoo)

# number of free DOFs
nDof = np.sum(fixDof)<jupyter_output><empty_output><jupyter_text>3. External loading  3.1. Static loading Uniformly distributed static loadings are applied on both beams; self-weight is neglected:<jupyter_code>q1 = 1.5e5   # [N] 2nd-floor beam
q2 = 1.5e5   # [N] 3rd-floor beam
H1 = 0.1*q1*Lf/2<jupyter_output><empty_output><jupyter_text>In 'staLoa' array, we input the static loading actions on the nodes of the FE mesh:<jupyter_code>Fy1 = q1*Lf/2
Mz1 = q1*Lf**2/12
Fy2 = q2*Lf/2
Mz2 = q2*Lf**2/12
staLoa = []
staLoa.append([0,0,0])         # node 0 [Dx, Dy, Rz] (displacements)
staLoa.append([0,0,0])         # node 1 [Dx, Dy, Rz] (displacements)
staLoa.append([H1,-Fy1,-Mz1])   # node 2 [Fx, Fy, Mz] (forces)
staLoa.append([0,-Fy1,Mz1])    # ...
staLoa.append([2*H1,-Fy2,-Mz2])
staLoa.append([0,-Fy2,Mz2])<jupyter_output><empty_output><jupyter_text>3.2. Seismic loading We also consider lumped masses at nodes 2, 3, 4 and 5 which are active for DOFs 1 (horizontal) only:<jupyter_code>gra = 9.81       # [m/s2]
mas1 = Fy1/gra   # [kg]
mas2 = Fy2/gra   # [kg]
inpMas = []
inpMas.append([2,mas1,mas1,mas1])   # node,DOF1,DOF2,DOF3
inpMas.append([3,mas1,mas1,mas1])
inpMas.append([4,mas2,mas2,mas2])
inpMas.append([5,mas2,mas2,mas2])

# mass active DOFs in the horizontal direction
r_a1 = [[1],[0],[0],[1],[0],[0],[1],[0],[0],[1],[0],[0]]
r_v = np.matrix(r_a1)<jupyter_output><empty_output><jupyter_text>Ground acceleration is recorded in file 'ps_fv.csv' [time, acceleration] and we import it in the 2D-array 'GM_a2':![Ground motion time history and response spectrum.](img/GroundMotion.pdf)<jupyter_code>GM_a2,timFin,timSteGM = impGM.impCsv('./GroundMotions/accePS-FV.csv')<jupyter_output><empty_output><jupyter_text>4. Resolution procedure  4.1. Static analysis<jupyter_code># external forces vector
fsta_v = strMod.extLoa(nDof,fixDof,staLoa)
    
# structural stiffness matrix K
K_m = strMod.stiMat(nDof,strEle,mesDat,eltPro,fixDof)

# compute static displacements at free DOFs
dis_v = linalg.solve(K_m,fsta_v)

print(dis_v)<jupyter_output>[[ 0.00131453]
 [-0.00602847]
 [-0.05424632]
 [ 0.00425264]
 [-0.00716835]
 [ 0.02642456]
 [ 0.13603423]
 [-0.00930051]
 [-0.11907751]
 [ 0.13202723]
 [-0.01073976]
 [ 0.0706236 ]]<jupyter_text>4.2. Modal analysis<jupyter_code># structural mass matrix M
M_m = strMod.masMat(nDof,inpMas,fixDof)

# eigen properties of the structure
ome_v,phi_m = modAna.eigPro(M_m,K_m)
print('ome_v =',ome_v)
print('phi1 =', phi_m[:,0])
# generalized masses
genMas_v = np.zeros((nDof,1))
for i in range(nDof):
    genMas_v[i,0] = np.transpose(phi_m[:,i])*M_m*phi_m[:,i]
    
# normalize eigen vectors wrt mass
phiNor_m = np.zeros((nDof,nDof))
for i in range(4,nDof):
    for j in range(nDof):
        phiNor_m[j,i] = phi_m[j,i]/np.sqrt(genMas_v[i,0])
    
#print(np.transpose(phiNor_m[:,6])*M_m*phiNor_m[:,6])
    
# participation factors
parFac_v = np.zeros((nDof,1))
for i in range(nDof):
    if genMas_v[i,0] == 0:
        parFac_v[i,0] = 0
    else:
        parFac_v[i,0] = np.transpose(phi_m[:,i])*M_m*r_v/genMas_v[i,0]

print(parFac_v)<jupyter_output>[[  1.85038500e+00]
 [  7.21726461e-01]
 [ -2.41641374e-16]
 [  1.81701563e-01]
 [ -5.72504177e-16]
 [ -1.48897639e-01]
 [ -1.18961907e-16]
 [ -9.51695256e-16]
 [ -5.46166933e-04]
 [ -6.90118872e-17]
 [ -1.40394832e-05]
 [ -9.59815643e-17]]<jupyter_text>4.3. Seismic analysis with response spectrum  4.4. Seismic analysis with time-stepping algorithm This is a dynamic nonlinear problem. We use Newmark time-stepping algorithm and Newton-Raphson iterative procedure to satisfy the equilibrium equation at each time step.The equilibrium equation reads:\begin{equation} \begin{pmatrix}  \mathbf{M}_{bb} & \mathbf{0} \\  \mathbf{0} & \mathbf{M} \end{pmatrix} \begin{pmatrix}  \mathbf{0} \\  \mathbf{a}(t) \end{pmatrix} + \begin{pmatrix}  \mathbf{C}_{bb} & \mathbf{C}_{b} \\  \mathbf{C}_{b} & \mathbf{C} \end{pmatrix} + \begin{pmatrix}  \mathbf{0} \\  \mathbf{v}(t) \end{pmatrix} + \begin{pmatrix}    \mathbf{f}_b^r(t) \\  \mathbf{f}^r(t) \end{pmatrix} = - \begin{pmatrix}  \mathbf{M}_{bb} & \mathbf{0} \\  \mathbf{0} & \mathbf{M} \end{pmatrix} \begin{pmatrix}  \mathbf{\Delta}_b \\  \mathbf{\Delta} \end{pmatrix} a^g(t) + \begin{pmatrix}  \mathbf{f}_b^s(t) \\  \mathbf{f}^s(t) \end{pmatrix}\end{equation}where the first lines (with subscript 'b') correspond to the $B$ DOFs that are fixed on the ground (base of the structure); the last $N$ lines correspond to the $N$ free DOFs. Also:- $\mathbf{M}$, $\mathbf{C}$ are the mass and stiffness matrices, and $\mathbf{f}^r$ is the resisting forces vector- $\mathbf{a}$, $\mathbf{v}$, and $\mathbf{d}$ are the vectors of the accelerations, velocities and displacements at the structural DOFs- $\mathbf{\Delta}a^g(t)$ is the seismic ground motion time history vector ($\Delta$ vector components are 1 for the horizontal DOFs and 0 at the other DOFs; $a^g(t)$ is the recorded seimic ground motion)- $\mathbf{f}^s$ is the vector of the external static loading (as opposed to the seismic loading); the components of the vector $\mathbf{f}_b^s$ are the reactions at the fixed DOFs.*Remark:* A lumped mass matrix is considered here (diagonal matrix).*Remark:* In what follows, we only solve the $N$ equations corresponding to the $N$ free DOFs. Then the $B$ other equations can be straightforwardly solved to compute the reactions as:\begin{equation} \mathbf{f}_b^s = \mathbf{C}_b \, \mathbf{v} + \mathbf{f}_b^r + \mathbf{M}_{bb} \mathbf{\Delta}_b a^g\end{equation}We know recast the problem introducing the residual\begin{equation} \mathbf{0} := \mathbf{r}(\mathbf{u},\mathbf{v},\mathbf{a}) = \mathbf{f}^e - \left( \mathbf{M} \, \mathbf{a} + \mathbf{C} \, \mathbf{v} + \mathbf{f}^r(\mathbf{u}) \right)\end{equation}where the external action is denoted as:\begin{equation} \mathbf{f}^e := -\mathbf{M}\mathbf{\Delta}a^g + \mathbf{f}^s\end{equation}Newmark time-stepping method:\begin{align} & \mathbf{u}_{n+1} = \mathbf{u}_n + \Delta t \, \mathbf{v}_n + \Delta t^2 \left( (\frac{1}{2}-\beta) \mathbf{a}_n + \beta \mathbf{a}_{n+1} \right) \\ & \mathbf{v}_{n+1} = \mathbf{v}_n + \Delta t \left( (1-\gamma) \mathbf{a}_n + \gamma \mathbf{a}_{n+1} \right)\end{align}Newton-Raphson iterative procedure to solve the nonlinear equation $\mathbf{r} = \mathbf{0}$:\begin{align} \mathbf{0} & := \mathbf{r}_{n+1}^{(i+1)} \approx \mathbf{r}_{n+1}^{(i)} + d\mathbf{r}_{n+1}^{(i+1)} \\  & = \mathbf{r}_{n+1}^{(i)} + \frac{\partial \mathbf{r}}{\partial \mathbf{u}}\Big|_{n+1}^{(i)} \, d\mathbf{u}_{n+1}^{(i+1)} + \frac{\partial \mathbf{r}}{\partial \mathbf{v}}\Big|_{n+1}^{(i)} \, d\mathbf{v}_{n+1}^{(i+1)} + \frac{\partial \mathbf{r}}{\partial \mathbf{a}}\Big|_{n+1}^{(i)} \, d\mathbf{a}_{n+1}^{(i+1)} \\  & = \mathbf{r}_{n+1}^{(i)} + \left( \frac{\partial \mathbf{r}}{\partial \mathbf{u}}\Big|_{n+1}^{(i)} + \frac{\partial \mathbf{r}}{\partial \mathbf{v}}\Big|_{n+1}^{(i)} \frac{d\mathbf{v}_{n+1}^{(i+1)}}{d\mathbf{a}_{n+1}^{(i+1)}} \frac{d\mathbf{a}_{n+1}^{(i+1)}}{d\mathbf{u}_{n+1}^{(i+1)}} + \frac{\partial \mathbf{r}}{\partial \mathbf{a}}\Big|_{n+1}^{(i)} \frac{d\mathbf{a}_{n+1}^{(i+1)}}{d\mathbf{u}_{n+1}^{(i+1)}} \right) d\mathbf{u}_{n+1}^{(i+1)} \\  & = \mathbf{r}_{n+1}^{(i)} - \left( \mathbf{Kt}_{n+1}^{(i)} + \frac{\gamma}{\beta \, \Delta t} \mathbf{C} + \frac{1}{\beta \, \Delta t^2} \mathbf{M} \right) d\mathbf{u}_{n+1}^{(i+1)}\end{align}where $\mathbf{Kt} := \partial \mathbf{f}^r / \partial \mathbf{u}$ is the so-called tangent stiffnes matrix.From the relation above, we can compute the updated displacements:\begin{equation} \mathbf{u}_{n+1}^{(i+1)} = \mathbf{u}_{n+1}^{(i)} + \left( \mathbf{Kt}^{\star (i)}_{n+1} \right)^{-1} \mathbf{r}_{n+1}^{(i)}\end{equation}where $\mathbf{Kt}^{\star} = \mathbf{Kt}_{n+1}^{(i)} + c_1 \mathbf{C} + c_2 \mathbf{M}$ with:\begin{equation} c_1 = \frac{\gamma}{\beta \, \Delta t} \qquad \textrm{and} \qquad c_2 = \frac{1}{\beta \, \Delta t^2}\end{equation}Then, the velocity and acceleration vectors can be updated as:\begin{align} & \mathbf{a}_{n+1}^{(i+1)} = \mathbf{a}_{n+1}^{(i)} + \frac{1}{\beta \, \Delta t^2} d\mathbf{u}_{n+1}^{(i+1)} \\ & \mathbf{v}_{n+1}^{(i+1)} = \mathbf{v}_{n+1}^{(i)} + \frac{\gamma}{\beta \Delta t} d\mathbf{u}_{n+1}^{(i+1)}\end{align}*Remark:* Here, the damping matrix $\mathbf{C}$ is considered as time independent, but it is time dependent for some damping models.  4.4.1. Parameterization<jupyter_code># structural damping ratio
strXi = 0.02

# time stepping algorithm (Newmark method)
timSte = 1e-2   # [s] time step
gamma = 0.5
beta = 0.25

# build Rayleigh damping matrix
alp,bet = strMod.rayDamCoe(strXi,ome_v)
C_m = alp*M_m + bet*K_m

# get time stepping coefficients c0, c1, and c2 (Newmark method)
c0 = 1/(beta*timSte)
c1 = c0*gamma
c2 = c0/timSte

# time steps ratio
timSteRat = timSteGM/timSte

# generalized stiffness matrix
K2_m = c1*C_m + c2*M_m
Kstar_m = K_m + K2_m<jupyter_output><empty_output><jupyter_text>4.4.2. Initialization<jupyter_code># time
time = 0

# initial forces vector
fi_v = np.matlib.zeros((nDof,1))  # !!!change if not zeros

# initial displacements vector
dis_v = np.matlib.zeros((nDof,1))  # !!!change if not zeros

# initial velocities
vel_v = np.matlib.zeros((nDof,1))

# initial ext. seismic forces vector
fsei_v = np.matlib.zeros((nDof,1))

# initial accelerations
acc_v = np.matlib.zeros((nDof,1))
f_v = fsei_v-C_m*vel_v-K_m*dis_v
for i in range(0,nDof):
    if M_m[i,i] == 0:
        acc_v[i] = 0
    else:
        acc_v[i] = 1/M_m[i,i]*f_v[i]

# output quantities
filDis = open('disp.txt','w')
filDis.write('{:1.3f},{:+.6E},{:+.6E},{:+.6E},{:+.6E} \n'
    .format(time,dis_v[0,0],dis_v[1,0],dis_v[2,0],dis_v[3,0]))<jupyter_output><empty_output><jupyter_text>4.4.3. Step-by-step resolution<jupyter_code># step forward
time += timSte

i = 1
while time < timFin:
#    print('time =',time)
    
    # interpolate ground acceleration
    xxx = i/timSteRat
    j = np.floor(xxx)
    j = np.int(j)
    res = xxx-j
    if res > 1e-12:
        groAcc = (1-res)*GM_a2[j][1] + res*GM_a2[j+1][1]
    else:
        groAcc = (1-res)*GM_a2[j][1]
    
    # compute external loading vector
    fsei_v = -M_m*r_v*groAcc*9.81

    # initialize
    vel1_v = (1-gamma/beta)*vel_v
    vel2_v = timSte*(1-gamma/(2*beta))*acc_v
    acc1_v = -1/(beta*timSte)*vel_v
    acc2_v = (1-1/(2*beta))*acc_v
    vel_v = vel1_v+vel2_v
    acc_v = acc1_v+acc2_v

    conv = False
    k = 0

    while conv == False:
        
        # compute residual
        res_v = fsei_v - (M_m*acc_v + C_m*vel_v + fi_v)

        # update increments
        disInc_v = linalg.solve(Kstar_m,res_v) 
        dis_v += disInc_v
        vel_v += c1*disInc_v
        acc_v += c2*disInc_v
        
        # for linear problems
        fi_v += K_m*disInc_v
        
        # for nonlinear problems
#        K1_m = strMod.stiMat(nDof,strEle,mesDat,eltPro,fixDof)
#        Kstar_m = K_m + K2_m
#        fi_v += K_m*disInc_v

        # check convergence
        norm = res_v.T*disInc_v
        if norm < 1e-12:
            conv = True
        elif k > 50:
            print('No convergence')
            break
        else:
            k += 1
    
    # output quantities
    filDis.write('{:1.3f},{:+.6E},{:+.6E},{:+.6E},{:+.6E} \n'
        .format(time,dis_v[0,0],dis_v[1,0],dis_v[2,0],dis_v[3,0]))

    # step forward
    i += 1
    time += timSte

filDis.close()
print('End of the computation')<jupyter_output>End of the computation
{'fraction_non_alphanumeric': 0.17645540209113206, 'fraction_numerical': 0.0467923066993675, 'mean_word_length': 4.365304709141275, 'pattern_counts': {'":': 0, '<': 61, '<?xml version=': 0, '>': 64, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '15', 'dir': 'jupyter-structured-clean-dedup', 'id': '439085', 'license': 'MIT', 'n_tokens_mistral': 6353, 'n_tokens_neox': 5970, 'n_words': 1402, 'path': 'TD2_2dFrameFEAdyn.ipynb', 'repo_name': 'jehelp/StructuralReliability', 'size': '15500'}