Background

The lin3DHexa20 class creates a linearized three-dimensional hexahedron element with twenty-nodes. Each node has three degrees-of-freedom, and each degree-of-freedom represents a translation in the X-, Y- and Z-directions. The figure provides a simple representation of this element in which the local coordinates are represented by the r-, s-, and t-axis, and (1), (2), (3), (4), \(\ldots\) (17), (18), (19) and (20) represent the nodes of the element respectively. The element mass, damping and stiffness matrices, as well as internal force vector are computed in global-coordinate. Numerical integration is performed to compute matrices and vectors.

The coordinates of any point \(x_1,x_2,x_3\) in the element is written in terms of the shape functions \(\textbf{N}_\textrm{i}(r, s, t)\) in the reference coordinates and the values of the i-th element nodal coordinates \(x^\textrm{(i)}_\textrm{j}\) as follows

\[ x_\textrm{j}^\textrm{e}(r, s, t) = \sum_\textrm{i = 1}^{\textrm{N}_\textrm{n}} x^\textrm{(i)}_\textrm{j} \, \textbf{N}_\textrm{i}^\textrm{e}(r, s, t) \,. \]

The same transformation is used to evaluate the displacements within the element–i.e., \(u_\textrm{j}^\textrm{e}(r, s ,t)\). Then, the displacement field is approximated as

\[ u_\textrm{j}^\textrm{e}(r, s, t) = \sum_\textrm{i = 1}^{\textrm{N}_\textrm{n}} u^\textrm{(i)}_\textrm{j} \, \textbf{N}_\textrm{i}^\textrm{e}(r, s, t) \,. \]

Since both the displacement field and the point element coordinates are approximated using the same shape function, this formulation is commonly known as isoparametric transformation. Computation of the strain-displacement matrix \(\textbf{B}^\textrm{e}(x_1,x_2,x_3)\) now requires to be mapped into the reference coordinate systems. Application of the chain rule for differentiation to the shape functions yields

\[ \begin{bmatrix} \displaystyle{\frac{\partial N_\textrm{i}}{\partial r,}} \\ \displaystyle{\frac{\partial N_\textrm{i}}{\partial s}} \\ \displaystyle{\frac{\partial N_\textrm{i}}{\partial t}} \end{bmatrix} = \begin{bmatrix} \displaystyle{\frac{\partial x}{\partial r,}} & \displaystyle{\frac{\partial y}{\partial r,}} & \displaystyle{\frac{\partial z}{\partial r,}} \\ \displaystyle{\frac{\partial x}{\partial s}} & \displaystyle{\frac{\partial y}{\partial s}} & \displaystyle{\frac{\partial z}{\partial s}} \\ \displaystyle{\frac{\partial x}{\partial t}} & \displaystyle{\frac{\partial y}{\partial t}} & \displaystyle{\frac{\partial z}{\partial t}} \end{bmatrix} \begin{bmatrix} \displaystyle{\frac{\partial N_\textrm{i}}{\partial x}} \\ \displaystyle{\frac{\partial N_\textrm{i}}{\partial y}} \\ \displaystyle{\frac{\partial N_\textrm{i}}{\partial z}} \end{bmatrix} = \mathbf{J}^\textrm{e} \begin{bmatrix} \displaystyle{\frac{\partial N_\textrm{i}}{\partial x}} \\ \displaystyle{\frac{\partial N_\textrm{i}}{\partial y}} \\ \displaystyle{\frac{\partial N_\textrm{i}}{\partial z}} \end{bmatrix} \,, \]

where \(\mathbf{J}^e: \mathbb{R}^3 \to \mathbb{R}^3\) is called the Jacobian matrix. Hence, derivatives in \(\textbf{B}^\textrm{e}(x_1,x_2,x_3)\) can be transformed into \(\textbf{B}^\textrm{e}(r,s,t)\) in the reference coordinates using the Jacobian matrix. The twenty shape functions \(\textbf{N}_\textrm{i}^\textrm{e}(r,s,t)\) required to evaluate the strain-displacement matrix \(\textbf{B}^\textrm{e}(r,s,t)\), shape function matrix \(\textbf{N}^\textrm{e}(r,s,t)\), and Jacobian matrix \(\textbf{J}^\textrm{e}(r,s,t)\) are:

\[ \textbf{N}_\textrm{1}^\textrm{e}(r, s ,t) = -\frac{1}{8} (1 - r) (1 - s) (1 - t) (2 + r + s + t) \,, \\ \textbf{N}_\textrm{2}^\textrm{e}(r, s ,t) = -\frac{1}{8} (1 + r) (1 - s) (1 - t) (2 - r + s + t) \,, \\ \textbf{N}_\textrm{3}^\textrm{e}(r, s ,t) = -\frac{1}{8} (1 + r) (1 + s) (1 - t) (2 - r - s + t) \,, \\ \textbf{N}_\textrm{4}^\textrm{e}(r, s ,t) = -\frac{1}{8} (1 - r) (1 + s) (1 - t) (2 + r - s + t) \,, \\ \textbf{N}_\textrm{5}^\textrm{e}(r, s ,t) = -\frac{1}{8} (1 - r) (1 - s) (1 + t) (2 + r + s - t) \,, \\ \textbf{N}_\textrm{6}^\textrm{e}(r, s ,t) = -\frac{1}{8} (1 + r) (1 - s) (1 + t) (2 - r + s - t) \,, \\ \textbf{N}_\textrm{7}^\textrm{e}(r, s ,t) = -\frac{1}{8} (1 + r) (1 + s) (1 + t) (2 - r - s - t) \,, \\ \textbf{N}_\textrm{8}^\textrm{e}(r, s ,t) = -\frac{1}{8} (1 - r) (1 + s) (1 + t) (2 + r - s - t) \,, \\ \textbf{N}_\textrm{9}^\textrm{e}(r, s ,t) = -\frac{1}{4} (1 - r^2) (1 - s) (1 - t) \,, \\ \textbf{N}_\textrm{10}^\textrm{e}(r, s ,t) = \frac{1}{4} (1 + r) (1 - s^2) (1 - t) \,, \\ \textbf{N}_\textrm{11}^\textrm{e}(r, s ,t) = \frac{1}{4} (1 - r^2) (1 + s) (1 - t) \,, \\ \textbf{N}_\textrm{12}^\textrm{e}(r, s ,t) = \frac{1}{4} (1 - r) (1 - s^2) (1 - t) \,, \\ \textbf{N}_\textrm{13}^\textrm{e}(r, s ,t) = \frac{1}{4} (1 - r^2) (1 - s) (1 + t) \,, \\ \textbf{N}_\textrm{14}^\textrm{e}(r, s ,t) = \frac{1}{4} (1 + r) (1 - s^2) (1 + t) \,, \\ \textbf{N}_\textrm{15}^\textrm{e}(r, s ,t) = \frac{1}{4} (1 - r^2) (1 + s) (1 + t) \,, \\ \textbf{N}_\textrm{16}^\textrm{e}(r, s ,t) = \frac{1}{4} (1 - r) (1 - s^2) (1 + t) \,, \\ \textbf{N}_\textrm{17}^\textrm{e}(r, s ,t) = \frac{1}{4} (1 - r) (1 - s) (1 - t^2) \,, \\ \textbf{N}_\textrm{18}^\textrm{e}(r, s ,t) = \frac{1}{4} (1 + r) (1 - s) (1 - t^2) \,, \\ \textbf{N}_\textrm{19}^\textrm{e}(r, s ,t) = \frac{1}{4} (1 + r) (1 + s) (1 - t^2) \,, \\ \textbf{N}_\textrm{20}^\textrm{e}(r, s ,t) = \frac{1}{4} (1 - r) (1 + s) (1 - t^2) \,. \]

The isoparametric transformation has the advantage that every element is treated in the same manner. In this regard, mass and stiffness matrices as well as force vector are evaluated using the same reference coordinates as follows

\[ \textbf{M}^\textrm{e} = \int_{-1}^1 \int_{-1}^1 \int_{-1}^1 {\textbf{N}^\textrm{e}} (r, s,t)^\top \, \rho^\textrm{e} \, \textbf{N}^\textrm{e}(r, s,t) \det \mathbf{J}^\textrm{e} \, dr \, ds \, dt \,, \]

\[ \textbf{K}^\textrm{e} = \int_{-1}^1 \int_{-1}^1 \int_{-1}^1 {\textbf{B}^\textrm{e}} (r, s,t)^\top \, \mathbb{C}^\textrm{e} \, \textbf{B}^\textrm{e}(r, s,t) \det \mathbf{J}^\textrm{e} \, dr \, ds \, dt \,, \]

\[ \textbf{F}^\textrm{e} = \int_{-1}^1 \int_{-1}^1 \int_{-1}^1 {\textbf{N}^\textrm{e}} (r, s,t)^\top \, \textbf{F}_b^\textrm{e} \det \mathbf{J}^\textrm{e} \, dr \, ds \, dt \,, \]

where \(\rho^\textrm{e}\) is the element density, \(\textbf{N}^\textrm{e}\) is the shape function matrix, \(\mathbb{C}^\textrm{e}\) is the material constitutive law, and \(\textbf{B}^\textrm{e}(r, s)\) is the well-known linear (material) strain-displacement matrix. Note that \(\textbf{M}^\textrm{e}, \textbf{K}^\textrm{e} \in \mathbb{R}^{\textrm{N}_\textrm{dof}^\textrm{e} \times \textrm{N}_\textrm{dof}^\textrm{e}}\) with \(\textrm{N}_\textrm{dof}^\textrm{e} = 60\) in this case.

REFERENCE:

Bathe, K. Jurgen, "Finite Element Procedures", Chapter 5: pages 342-348. Prentice-Hall, 1996.

Pre-Analysis

The python Pre-Analysis in the 01-Pre_Process/Method/Attach.py file provides with an interface to populate the Entities dictionary. This file contains several functions to populate specific fields. For example, to create a lin3DHexa20, using json format, use:

addElement(tag, name='lin3DHexa20', conn, attributes):
- tag : The identifier of this element, i.e., tag > -1
- name : Seismo-VLAB element class name
- conn : Connectivity array of this element
- attributes : Specific properties for the created element, for example
  - 'material' is material identifier.
  - 'rule' the integration rule: 'Gauss', 'Lobatto'.
  - 'np' Number of quadrature-points used for the integration.
Example

A LIN3DHEXA20 element can be defined using the python interface as follows:
SVL.addElement(tag=1, name='lin3DHexa20', conn=[1,2,3,4,5,6,7,8,9,10,11,12,13,14,14,16,17,18,19,20], attributes={'np': 27, 'rule': 'Gauss', 'material': 1})

Application
Please refer to the I03-ST_Lin_3DPointLoad_Elastic_Hexa20.py file located at 03-Validations/01-Debugging/ to see an example on how to define a lin3DHexa20 element.

On the contrary, the 01-Pre_Process/Method/Remove.py file provides with an interface to depopulate the Entities dictionary. For example, to remove an already define Element, use:

delElement(tag):
- tag : The identifier of the element to be removed, i.e., tag > -1

Run-Analysis

The C++ Run-Analysis in the 02-Run_Process/04-Elements/10-Hexahedron/lin3DHexa20.cpp file provides the class implementation. A lin3DHexa20 element is created using the built-in json parse-structure provided in the Driver.hpp. A lin3DHexa20 is defined inside the "Elements" json field indicating its "Tag" as follows,

{
    "Elements": {
        "Tag": {
            "name" : "LIN3DHEXA20",
            "conn" : [ ],
            "attributes": {
                "np": int,
                "rule": str,
                "material": int
            }
        }
    }
}

Variable	Description
`Tag`	Unique element object identifier.
`conn`	The element connectivity node array.
`np`	Number of quadrature-points used for the integration.
`rule`	The quadrature integration name `rule`=GAUSS, LOBATTO.
`material`	The biaxial material identifier.

Attention: This formulation only allows small deformations.; Material non-linearity can be used in this element.; The material is considered to be the same at each quadrature point.; Currently the number of integration points np= 1,8,27,64,125,216,343 are supported.; The domain reduction force method is compatible with this element.; For the element size \(\Delta_\textrm{h}\) in wave propagation problems, one can consider about 12 points in one wavelength (this is 6 elements per wavelength) (the one with the maximum frequency). Once the mesh size is fixed, satisfying CFL condition to determine \(\Delta_\textrm{t}\) (i.e, \(\Delta_\textrm{t} \leq 0.9\,\Delta_\textrm{h}/V_p\), with \(V_p\) the p-wave velocity) for time marching is required.

Example

A LIN3DHEXA20 with 27 integration points between nodes 1 to 20, made of material 0 is constructed as:
{ "Elements": { "1": { "name" : "LIN3DHEXA20", "conn" : [1,2,3,4,5,6,7,8,9,10,11,12,13,14,14,16,17,18,19,20], "attributes": { "np": 27, "rule": "Gauss", "material": 0 } } } }