Finite Element Tutorial in Python

Introduction to the Python Laboratory for Finite Element Analysis (pyLabFEA). This tutorial covers FEA with elements in one or two dimensions to evaluate the mechanical behavior of laminate geometries with linear elastic materials. See the online documentation of the pyLabFEA package for detailed information. NumPy (http://www.numpy.org) is used for mathematical operations on arrays and matplotlib (https://matplotlib.org/) for the visualization of numerical results.

Author: Alexander Hartmaier, ICAMS / Ruhr-Universität Bochum

March 2020

This work is licensed under a Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International License (CC-BY-NC-SA)

The pyLabFEA package comes with ABSOLUTELY NO WARRANTY. This is free software, and you are welcome to redistribute it under the conditions of the GNU General Public License (GPLv3)

Python Laboratory for Finite Element Analysis (pyLabFEA)

In the first step, the pyLabFEA package, containing the modules Model and Material, is imported.

import pylabfea as fea
print('pyLabeFEA version', fea.__version__)

pyLabeFEA version 4.2

Model

The class Model is invoked to create an objecte for the finite element model to be used. This object contains all attributes of the model and defines the methods. To use the model for FEA, the following steps are performed: (i) geometry definition, (ii) material definition, (iii) applying boundary conditions, (iv) meshing, (v) solving of linear equations, and (vi) calculation of element stresses and strains by calling the corresponding methods of the class Model.

Parameters:
dim (optional): dim=1 for 1-d model, dim=2 for 2-d model (default: dim=1)
planestress (optional, used only in 2-d models): True if plane stress conditions shall be considered (default: planestress=False)

Methods:
geom: define (laminate) geometry with sections
assign: assign material properties to sections
bcleft, bcright, bctop, bcbot: define boundary conditions
mesh: generate mesh (1-d elements with linear or quadratic shape functions or 2-d mesh with linear quadriliteratl elements)
solve: setup stiffness matrix and calculate displacements satisfying mechanical equilibrium, i.e. $f=0$ for internal nodes, for given boundary conditions
calc_global: calculate global stress and strain based on displacements and residual forces of boundary nodes
plot: show graphical representation of mechanical quantities mapped into deformed shape

fem = fea.Model()  # create element of class Model

Geometry

Define the geometry of the model with method geom. The geometry is subdivided into sections, representing a 1-d or 2-d laminate structure normal to the $x$-direction. The sections are given as a list with the absolute length of each section. The total length of the model is the sum of all section lengths. All lengths are given in units of 1 mm.

Parameters:
sect: list with absolute length of each section in x-direction
LY (optional): height of model in y-direction (default: LY=1)
LZ (optional): thickness of model in z-direction (default: LZ=1)

fem.geom([3, 2, 4])  # define sections in absolute lengths

Materials

Material properties are given either in terms of anisotropic elastic constants $C_{11}$, $C_{12}$, $C_{44}$, or – in case of an elastically isotropic material – as Young's modulus $E$ and Poisson's ratio $\nu$. The parameters must be given with their name. The material is created by applying the class Material, which defines a container for the attributes and methods describing material properties and behavior. Elastic properties are defined by calling the method elasticity with the appropriate set of parameters.

Class:
Material

Methods:
elasticity: define linear elastic material either by anisoptropic or by isotropic constants. Note: One complete set of elastic parameters must be provided.
Paramaters:
C11, C12, C44 (optional): anisotropic elastic constants
E, nu (optional): Young's modulus and Poisson's ratio (isotropic elasticity)
one of both parameter sets must be provided

Attributes:
sets internal variables for $C_{11}$, $C_{12}$, $C_{44}$, $E$ and $\nu$

mat1 = fea.Material()              # create element of class Material
mat1.elasticity(E=200.e3, nu=0.3)  # assign isotropic elastic properties
mat2 = fea.Material()
mat2.elasticity(E=100.e3, nu=0.3)
mat3 = fea.Material()
mat3.elasticity(E=500.e3, nu=0.3)
print(f'MAT1: C11={mat1.C11:9.3e}, C12={mat1.C12:9.3e}, C44={mat1.C44:9.3e}, E={mat1.E:9.3e}, nu={mat1.nu:4.2f}')
print(f'MAT2: C11={mat2.C11:9.3e}, C12={mat2.C12:9.3e}, C44={mat2.C44:9.3e}, E={mat2.E:9.3e}, nu={mat2.nu:4.2f}')
print(f'MAT3: C11={mat3.C11:9.3e}, C12={mat3.C12:9.3e}, C44={mat3.C44:9.3e}, E={mat3.E:9.3e}, nu={mat3.nu:4.2f}')

MAT1: C11=2.692e+05, C12=1.154e+05, C44=7.692e+04, E=2.000e+05, nu=0.30
MAT2: C11=1.346e+05, C12=5.769e+04, C44=3.846e+04, E=1.000e+05, nu=0.30
MAT3: C11=6.731e+05, C12=2.885e+05, C44=1.923e+05, E=5.000e+05, nu=0.30

Assigning materials to model sections

Each section has to be assigned a material, which is accomplished by calling the according method of the class Model with the list of materials, which must have the same shape as the list of sections in the geometry definition.

Class: Model
Method:
assign: assign a material to each section<br/assign: assign a material to each section
Paramaters:
mats: list of materials, each defined by class Material, with same shape as in method geom. Each section is assigned to the according material.
Attributes
fem.mat: list with pointer to material for each section

fem.assign([mat1, mat2, mat3])  # assign materials to sections of model
print(fem.mat)

[<pylabfea.material.Material object at 0x15b92a8e0>, <pylabfea.material.Material object at 0x15b92aca0>, <pylabfea.material.Material object at 0x15b92acd0>]

Boundary Conditions

For left-hand-side and bottom nodes a displacement is given, typically $u_{\rm left}=u_{\rm bot}=0$.

The type of boundary conditions for right-hand-side and top nodes is specified as either displacement (type='disp') or force (type='force'). Value and type must be specified as parameters in this sequence.

For 1-d models only bcleft and bcright can be called.

Method:
bcleft: set displacement for lhs node(s)
bcbot: set displacement for bottom node(s), has no function for 1-d models
Paramaters:
scalar value: Displacement for boundary node(s)

Method:
bcright: set displacement or force for rhs node(s)
bctop: set displaement or force for top node(s), has no function for 1-d models
Paramaters:
scalar value: Displacement of force for boundary nodes
bctype: string indicating type of boundary conditions: 'disp' for displacements, 'force' for forces

fem.bcleft(0.)  # define displacement boundary condition on lhs node (u_x given as parameter)
fem.bcright(0.1*fem.lenx, 'disp')  # rhs node is subject to 10% strain (disp = 0.1 * length of model)

Meshing

Create a structured mesh with 1-d elements or 2-d quadrilateral elements. Nodes are situated at the corners of the elements and in the middle of the edges for quadratic shape functions. The finite elements of the model are stored in a list with elements of the subclass element with its own attributes and methods, which are used to setup the element shape function.

Method (class Model):
mesh
Parameters:
NX (optional): number of elements in x-direction (default: NX=10)
NY (optional): number of elements in y-direction (default: NY=1), only effective for 2-d models
SF (optional): degree of shape functions: 1=linear (default), 2=quadratic

Subclass: `Element` (parent class: `Model`)
Methods:
calc_Bmat: Calculate B matrix at Gauss point
calc_Kel: calculate element stiffness matrix
Attributes:
nodes: list of nodes of this element
Lelx, Lely: x- and y-dimensions of element
ngp: number of Gauss points
gpx, gpy: np.arrays of x- and y-locations of Gauss points in element
Bmat: list of B-matrices at Gauss points
Vel: volume of element
Kel: element stiffness matrix
Jac: determinant of Jacobian
wght: weight factor for each Gauss point in integration
Sect: Section of model in which element is located
Mat: Pointer to class Material to be applied to this element
CV: Voigt stiffness matrix of element material
eps: average element (total) strain in Voigt notation
sig: average element stress in Voigt notation

fem.mesh(NX=9)  # create mesh
print(f'nodal positions: {fem.npos.round(decimals=2)}')
print(f'nodes belonging to element #0: {fem.element[0].nodes}')
print('stiffness matrix of element #0:')
print(fem.element[0].Kel.round(decimals=2))

nodal positions: [0. 1. 2. 3. 4. 5. 6. 7. 8. 9.]
nodes belonging to element #0: [0, 1]
stiffness matrix of element #0:
[[ 269230.77 -269230.77]
 [-269230.77  269230.77]]

Solver

The system of linear equation defined by the stiffness matrix is solved by invoking the Linear Algebra (linalg) subroutine solve of NumPy. The solution yields the nodal displacements fulfilling mechanical equilibrium for the given boundary conditions. The nodal forces of internal nodes are zero, forces on boundary nodes are either given as boundary conditions or are residual forces on fixed boundary nodes. If no boundary conditions are specified, the force on such boundary nodes will be zero.

After solution is completed, the following attributes of the class `Model` are defined.
Attributes:
u: list of nodal displacements
f: list of nodal forces
element[i].eps, element[i].sig: Voigt strain and stress tensor of element #i

fem.solve()  # solve linear system of equations
print('nodal displacements u = ', fem.u.round(decimals=2))
print('nodal forces f = ', fem.f.round(decimals=2))
print('element strain eps = ', [element.eps[0].round(decimals=2) for element in fem.element])
print('Voigt strain tensor for element #0 = ', fem.element[0].eps.round(decimals=2))

nodal displacements u =  [0.   0.1  0.21 0.31 0.52 0.73 0.77 0.82 0.86 0.9 ]
nodal forces f =  [-28175.31      0.        0.        0.        0.        0.        0.
      0.       -0.    28175.31]
element strain eps =  [0.1, 0.1, 0.1, 0.21, 0.21, 0.04, 0.04, 0.04, 0.04]
Voigt strain tensor for element #0 =  [0.1 0.  0.  0.  0.  0. ]

Post-processing

Calculate global stresses and strains based on displacements and residual forces on boundary nodes. Provide a graphical output of the specified field, the parameter can be either 'strain' or 'stress'. A rectangle is plotted for each element, with the color given by the specified field and a color map.

Class `Model`
Method:
calc_global: Calculate global stresses and strains based on displacements and residual forces of boundary nodes.
Attributes:
glob: Python dictionary with global strains and stresses, contains the elements:
'ebc1', 'ebc2', 'sbc1', 'sbc2' : global strain and stress calculated from reaction forces and displacements of boundary nodes (type: float)
'eps', 'epl', 'sig', : global strain, plastic strain, and stress tensors homogenized from element solutions (type: Voigt tensor)

Method:
plot: Calculate stresses and strains for each element and stores them in lists assigned to model.
Parameters:
fsel: field selector, string with value 'strain1' for $\epsilon_{11}$, 'strain2' for $\epsilon_{22}$, 'stress1' for $\sigma_{11}$, 'stress2' for $\sigma_{22}$.
mag (optional): scaling factor (magnification) for displacements (default: mag=10)
cdepth (optional): depth of colormap (default: cdepth=20)
showmesh (optional): Boolean variable to set/unset plotting of lines for element edges (default: showmesh=True)
shownodes (optional): Boolean variable to set/unset plotting of nodes (default: shownodes=True)

See the online documentation of the pyLabFEA package for a full description of the attributes and methods of the class Model.

fem.calc_global()  # calculate global stress and strain
print(f'global stress: sig_11 = {fem.glob["sig"][0].round(decimals=3)} MPa')
print(f'global strain: eps_11 = {(fem.glob["eps"][0]*100).round(decimals=2)} %')
print('\nGraphical output of model with magnification factor 10')
fem.plot('strain1')  # create plot

global stress: sig_11 = 28175.313 MPa
global strain: eps_11 = 10.0 %

Graphical output of model with magnification factor 10

2-dimensional model 1: iso-strain

A 2-d model with three different sections is created and subjected to a uniaxial stress parallel to the sections (y-direction), resulting in an iso-strain condition, because all sections are subject to the same total strain as the entire model. The effective eleastic properties are calculated numerically and compared to the results of the Voigt homogenization rule (iso-strain assumption).

# 2-d model, iso-strain
comp=fea.Model(dim=2, planestress=True)
comp.geom([3,2,4], LY=9)        # define sections in absolute lengths
comp.assign([mat1, mat2, mat3]) # assign a material to each section
comp.bcleft(0.)                 # define boundary conditions
comp.bcright(0., 'force') 
comp.bcbot(0.)
comp.bctop(0.1*comp.leny, 'disp') # apply eps_22 strain at top boundary
comp.mesh(NX=18, NY=18)         # create structured mesh with 18x18 elements
comp.solve()                    # solve system of equations for mechanical equilibrium
comp.calc_global()              # evaluate global properties
mod_stiff = comp.glob['sig'][1]/comp.glob['eps'][1]
voigt_stiff = 3.*mat1.E/9 + 2*mat2.E/9 + 4*mat3.E/9 # weighted average of Young's moduli wrt volume fractions
print('2-d Model 1: iso-strain, plane stress, eps_22=10%')
print('Global strain: ',comp.glob['eps'][0].round(decimals=3), comp.glob['eps'][1].round(decimals=3))
print('Element strain Section 1: ', comp.element[0].eps.round(decimals=3))
print('Element strain Section 2: ', comp.element[6*18].eps.round(decimals=3))
print('Element strain Section 3: ', comp.element[10*18].eps.round(decimals=3))
print('Global stress (MPa): ',comp.glob['sig'][0].round(decimals=3), comp.glob['sig'][1].round(decimals=3))
print('Element stress (MPa) Section 1: ', comp.element[0].sig.round(decimals=3))
print('Element stress (MPa) Section 2: ', comp.element[6*18].sig.round(decimals=3))
print('Element stress (MPa) Section 3: ', comp.element[10*18].sig.round(decimals=3))
print('Stiffness (MPa): ',mod_stiff)
print('Target (MPa): ', voigt_stiff)
print('Error: ', 1.-mod_stiff/voigt_stiff)
print('----------------------------------------')

comp.plot('strain1', mag=1, shownodes=False)
comp.plot('strain2', mag=1, shownodes=False)
comp.plot('stress1', mag=1, shownodes=False)
comp.plot('stress2', mag=1, shownodes=False)

2-d Model 1: iso-strain, plane stress, eps_22=10%
Global strain:  -0.03 0.1
Element strain Section 1:  [-0.03  0.1  -0.03  0.    0.    0.  ]
Element strain Section 2:  [-0.03  0.1  -0.03  0.    0.   -0.  ]
Element strain Section 3:  [-0.03  0.1  -0.03  0.    0.    0.  ]
Global stress (MPa):  0.0 31111.111
Element stress (MPa) Section 1:  [    0. 20000.     0.     0.     0.     0.]
Element stress (MPa) Section 2:  [    0. 10000.     0.     0.     0.    -0.]
Element stress (MPa) Section 3:  [    0. 50000.     0.     0.     0.     0.]
Stiffness (MPa):  311111.1111111126
Target (MPa):  311111.1111111111
Error:  -4.6629367034256575e-15
----------------------------------------

2-dimensional model 2: iso-stress

A 2-d model with three different sections is created and subjected to a uniaxial strain perpendicular to the sections (x-direction) resulting in an iso-stress condition, because all section are subject to the same stress in loading direction as the entire model. The effective eleastic properties are calculated numerically and compared to the results of the Reuss homogenization rule (iso-stress assumption).

# 2-d model, iso-stress
comp2=fea.Model(dim=2, planestress=False)
comp2.geom([3,2,4], LY=9)  # define sections in absolute lengths
comp2.assign([mat1, mat2, mat3])
comp2.bcleft(0.)
comp2.bcright(0.1*comp2.lenx, 'disp') 
comp2.bcbot(0.)
comp2.bctop(0., 'disp')
comp2.mesh(NX=18, NY=18)
comp2.solve()
comp2.calc_global()
mod_stiff = comp2.glob['sig'][0]/comp2.glob['eps'][0]
reuss_stiff = 1./(3/(9*mat1.C11) + 2/(9*mat2.C11) + 4/(9*mat3.C11)) # Reuss average of stiffness
print('2-d Model 2: isostress, eps_11=10%')
print('Global strain: ',comp2.glob['eps'][0].round(decimals=3), comp2.glob['eps'][1].round(decimals=3))
print('Element strain Section 1: ', comp2.element[0].eps.round(decimals=3))
print('Element strain Section 2: ', comp2.element[6*18].eps.round(decimals=3))
print('Element strain Section 3: ', comp2.element[10*18].eps.round(decimals=3))
print('Global stress: ',comp2.glob['sig'][0].round(decimals=3), comp2.glob['sig'][1].round(decimals=3))
print('Element stress (MPa) Section 1: ', comp2.element[0].sig.round(decimals=3))
print('Element stress (MPa) Section 2: ', comp2.element[6*18].sig.round(decimals=3))
print('Element stress (MPa) Section 3: ', comp2.element[10*18].sig.round(decimals=3))
print('Stiffness (MPa): ',mod_stiff)
print('Target (MPa): ', reuss_stiff)
print('Error: ', 1.-mod_stiff/reuss_stiff)
print('----------------------------------------')

comp2.plot('strain1', mag=1, shownodes=False)
comp2.plot('strain2', mag=1, shownodes=False)
comp2.plot('stress1', mag=1, shownodes=False)
comp2.plot('stress2', mag=1, shownodes=False)

2-d Model 2: isostress, eps_11=10%
Global strain:  0.1 0.0
Element strain Section 1:  [ 0.105 -0.     0.     0.     0.     0.   ]
Element strain Section 2:  [ 0.209 -0.     0.     0.     0.     0.   ]
Element strain Section 3:  [ 0.042  0.     0.     0.     0.    -0.   ]
Global stress:  28175.313 12075.134
Element stress (MPa) Section 1:  [28175.313 12075.134 12075.134     0.        0.        0.   ]
Element stress (MPa) Section 2:  [28175.313 12075.134 12075.134     0.        0.        0.   ]
Element stress (MPa) Section 3:  [28175.313 12075.134 12075.134     0.        0.       -0.   ]
Stiffness (MPa):  281753.13059033663
Target (MPa):  281753.1305903399
Error:  1.1546319456101628e-14
----------------------------------------

Summary

FE models with linear elastic materials have been created and applied to analyze the mechanical response of elastic structures to given boundary conditions. The properties of linear-elastic structures are analysed more closely in the tutorial Composites.