Membrane element

Membrane elements are special 3-node particles in DEM. They have their mass lumped into nodes and respond to plane-stress stretching (CST element) and bending (DKT element; optional, enabled by setting the bending flag). These two elements are superposed. The element has constant thickness. Each node has 6 DoFs:

1. two in-plane translations (gathered in woo.fem.Membrane.uXy for all nodes as a 6-vector), used by the CST element;

2. one out-of-plane translation – defines rigid motion of the whole element and causes no internal response;

3. two in-plane rotations (gathered in woo.fem.Membrane.phiXy for all nodes as a 6-vector), used by the DKT element;

4. one out-of-plane rotation (drilling DoF), which is ignored.

Mass and intertia of nodes are computed and lumped into nodes (this only requires geometry and material density), but it must be requested by calling woo.dem.DemData.setOriMassInertia on the node, once all its attached particles are in place. Details below.

Element coordinate system

When the element is first used, its configuration defines the reference configuration which is stored in woo.fem.Membrane.refPos and woo.fem.Membrane.refRot.

The element coordinate system (held in woo.fem.Membrane.node) is always such that

1. local origin coindices with centroid of the (deformed) triangle;

2. the local \(z\)-axis is normal to the plane defined by element nodes in their current positions;

3. local \(x\) and \(y\) axes are oriented in a way which minimizes displacements of nodes in the local plane relative to the reference positions, and tries to represent their mutual displacements as rigid body motion as much as possible.

   This point is trivially satisfied in the reference configuration (first time the element is used) when displacements are zero with arbitrary orientation of the local \(x\) and \(y\) axes. (The implementation currently chooses to rotate the \(x\)-axis so that it passes through the first node.)

   Otherwise, we are looking for a best-fit rotation; the algorithm is described in [Fel99b] in Appendix C. Nodes positions are projected onto (not yet properly rotated) local \(xy\) plane, and best-fit angle is computed (from reference positions \(X_i\), \(Y_i\) and current local positions \(x_i\), \(y_1\)) as

   \[\tan\theta_3=\frac{x_1 Y_1 + x_2 Y_2 + x_3 Y_3 - y_1 X1 - y_2 X_2 - y_3 X_3}{x_1 X_1 + x_2 X_2 + y_3 Y_3 + y_1 Y_1 + y_2 Y_2 + y_3 Y_3}\]

   where the numerator and denominator are computed separately in the implementation and \(\tan\phi_3\) is obtained vi via the atan2 function as to preserve quadrant information and handle zero denominator correctly. The \(\theta_3\) rotation is then added to the previous rotation of the \(xy\) plane, obtaining the best-fit coordinate system.

Once the local coordinate system is established, displacements are evaluated. Evaluating woo.fem.Membrane.uXy is done as trivial subtraction of reference nodal positions from the current ones in the local system. Rotations woo.fem.Membrane.phiXy are computed by rotation “subtraction” (conjugate product) using quaternions (woo.fem.Membrane.refRot holds the rotation necessary to obtain each node orientation from global element orientation). See the source code for details.

Constant Strain Triangle (CST) Element

The CST element is implemented following [Fel98], chapter 15. The strain-displacement matrix is computed as (15.17):

\[\begin{split}\mat{B}=\frac{1}{2A}\begin{pmatrix}y_{23} & 0 & y_{31} & 0 y_{12} & 0 \\ 0 & x_{23} & 0 & x_{12} & 0 & x_{21} \\ x_{32} & y_{23} & x_{12} & y_{31} & x_{21} & y_{12}\end{pmatrix}\begin{pmatrix} u_{x1} \\ u_{y1} \\ y_{x2} \\ u_{y2} \\ u_{x3} \\ u_{y3}\end{pmatrix}.\end{split}\]

The elastic stiffness matrix for plane-stress conditions reads

\[\begin{split}\mat{E}=\frac{E}{1-\nu^2}\begin{pmatrix}1 & \nu & 0 \\ \nu & 1 & 0 \\ 0 & 0 & \frac{1-\nu}{2}\end{pmatrix}\end{split}\]

where \(E\) is Young modulus and \(\nu\) is the Poisson ratio.

Since membrane thickness \(h\) is constant over the element, we may write the stiffness matrix as ([Fel98], (15.21)):

(35)\[\mat{K^{\mathrm{CST}}}=A h \mat{B}^T \mat E \mat B.\]

This matrix is stored as woo.fem.Membrane.KKcst.

Displacements \(\vec{u}\) (6-vector) are computed by subtracting reference positions from the current nodal positions in best-fit local coordinates. Nodal forces (6-vector) are

(36)\[\begin{split}\vec{F}^{\mathrm{CST}}=\begin{pmatrix}F_{x1}\\F_{y1}\\F_{x2}\\F_{y2}\\F_{x3}\\F_{y3}\end{pmatrix}=\mat{K}^{\mathrm{CST}}\vec{u}.\end{split}\]

This is an example of a CST-only mesh, no bending (examples/membrane1.py):

Stresses

Stresses don’t come up as such in the FEM, being hidden in eqs. (35) and (36):

\[\vec{F}^\mathrm{CST}=A h \mat{B}^T \underbrace{\mat{E}\mat{B}\vec{u}}_{\vec{\sigma}^\mathrm{CST}}.\]

Note that \(\vec{\sigma}^\mathrm{CST}\) are in local coordinates, that’s why they are then converted (by left-multiplication) to forces in global coordinates through \(Ah\mat{B}^T\).

In the implementation, the \(\mat{E}\mat{B}\) product is normally not stored in element, since it is not used in reaction computation. For that reason, it must be explicitly requested by setting Membrane.enableStress before stiffness matrix is evaluated. When stiffness matrix is then built, the \(\mat{E}\mat{B}\) product will be stored in Membrane.EBcst as side-effect. It can be then used to compute stresses either explicitly by multiplying EBcst × uXy, or using the shorthand Membrane.stressCst.

Note

Stresses returned are in element-local coordinates and they are constant over the entire element; the order is \((\sigma_x,\sigma_y,\sigma_{xy})\); to transform to global coordinates, call Membrane.stressCst(glob=True) (the result will be stress tensor represented as 3×3 matrix).

Discrete Krichhoff Triangle (DKT) Element

The DKT (bending) element is implemented following [BBH80]. The bending elasticity matrix for constant thickness is computed as

\[\begin{split}\mat{D}_b=\frac{E h^3}{12(1-\nu^2)}\begin{pmatrix}1 & \nu & 0 \\ \nu & 1 & 0 \\ 0 & 0 & \frac{1-\nu}{2}\end{pmatrix}\end{split}\]

where the thickness \(h\) may be different thickness than the one used for the CST element (different thicknesses are mostly not useful). The strain-displacement transformation matrix for bending reads ([BBH80], (30))

\[\begin{split}\mat{B}_{b}(\xi,\eta)=\frac{1}{2A}\begin{pmatrix}y_{31}\vec{H}_{x,\xi}^T+y_{12}\vec{H}_{x,\eta}^T \\ -x_{31}\vec{H}_{y,\xi}^T-x_{12}\vec{H}_{y,\eta}^T \\ -x_{31}\vec{H}_{x,\xi}^T-x_{12}\vec{H}_{x,\eta}^T+y_{31}\vec{H}_{y,\xi}^T+y_{12}\vec{H}_{y,\eta}^T \end{pmatrix}.\end{split}\]

where \(x_{ij}\equiv x_i-x_j\), \(y_{ij}\equiv y_i-y_j\) and \((\xi,\eta)\) are natural (≡barycentric) coordinates on the triangle. Refer to Appendix A of [BBH80] (or the source code of woo.fem.Membrane) for formulation of \(\vec{H}_{x,\xi}\), \(\vec{H}_{x,\eta}\), \(\vec{H}_{y,\xi}\), \(\vec{H}_{y,\eta}\), which are rather convolved.

The stiffness matrix is integrated over the triangle as ([BBH80], (31))

\[\mat{K}^{\mathrm{DKT}}=2A\int_0^1\int_0^{1-\eta}\mat{B}^T\mat{D}_b\mat{B}\d\xi\,\d\eta\]

and can be integrated numerically (see [Kan04], pg. 48) using the Gauss quadrature (Gauss points are mid-points of triangle sides, \(w\) is vector of point weights) over the unit triangle as

\begin{align*} \vec{\xi}&=\left(\frac{1}{2},\frac{1}{2},0\right),\\ \vec{\eta}&=\left(0,\frac{1}{2},\frac{1}{2}\right),\\ \vec{w}&=\left(\frac{1}{3},\frac{1}{2},\frac{1}{2}\right),\\ \mat{K}^{\mathrm{DKT}}&=2A\sum_{j=1}^{3}\sum_{i=1}^{3} w_j w_i \mat{B}_b^T(\xi_i,\eta_j) \mat{D}_b \mat{B}_b(\xi_i,\eta_j). \end{align*}

Generalized nodal forces are computed as

\[\begin{split}\vec{F}^{\mathrm{DKT}}=\begin{pmatrix} F_{z1} \\ T_{x1} \\ T_{y1} \\ F_{z2} \\ T_{x1} \\ T_{x2} \\ F_{z3} \\ T_{x3} \\ T_{y3}\end{pmatrix}=\mat{K}^{\mathrm{DKT}}\begin{pmatrix}u_{z1}\equiv0\\ \phi_{x1} \\ \phi_{y1} \\ u_{z2}\equiv0 \\ \phi_{x2} \\ \phi_{y2} \\ u_{z3}\equiv 0 \\ \phi_{x3} \\ \phi_{y3}\end{pmatrix},\end{split}\]

Since out-of-plane translations \(u_{zi}\) are always zero (they determine rigid body rotation of the element), we may condense those rows from \(\mat{K}_{\mathrm{DKT}}\) away, making it 9×6 rather than 9×9, and removing zero elements from the generalized displacement vector as well. Note that corresponding DoFs may nevertheless have non-zero forces \(F_{zi}\).

The examples/membrane2.py script shows the combined response of CST and DKT elements:

Stresses

Warning

This is not functional yet in the code.

Moments can be computed as ([BBH80], (32))

\[\mat{M}(x,y)=\mat{D}_b\mat{B}(x,y)\vec{U}\]

Similar as with with CST, stress evaluation needs some extra data (namely, the \(\mat{D}_b\mat{B}(x,y)\) term) of which availability must be requested by the user through Membrane.enableStress; DBdkt will be stored when stiffness matrix is computed as side-effect; stress can be obtained by hand-multiplication with rotations (phiXy) or using the shorthand function Membrane.stressDkt.

Total nodal forces

Total nodal forces are superimposed from both elements and also from contact forces on the particle and from a special load type on membranes, hydrostatic pressure. We note force and torque contributions as \(\Delta \vec{F}_i\), \(\Delta \vec{T}_i\) for \(i\)-th node, \(i\in\{1,2,3\}\).

Hydrostatic pressure

Elements may be under surface load which is always perpendicular to the element plane, thus representing hydrostatic pressure. Such pressure \(p\) is distributed into nodes as

\[\begin{split}\Delta\vec{F}_i = \begin{pmatrix} 0 \\ 0 \\ p\frac{A^*}{3} \end{pmatrix}\end{split}\]

where \(A^*\) is the current (not reference) element surface area.

Contact forces

Contact forces defined by contact point \(\vec{c}\), force \(\vec{F}_c\) and torque \(\vec{T}_c\) are simply distributed onto all nodes, so that each node \(i\) (positioned at \(\vec{x}_i\)) receives for each contact the contribution

\begin{align*} \Delta\vec{F}_i&=\frac{\vec{F}_c}{3}, \\ \Delta\vec{T}_i&=(\vec{x}_i-\vec{c})\times\vec{F}_c+\vec{T}_c. \end{align*}

CST + DKT

Generalized forces from elements are distributed to nodes as follows:

\begin{align*} \Delta\vec{F}_i&=\begin{pmatrix} \vec{F}^{\mathrm{CST}}_{xi} \\ \vec{F}^{\mathrm{CST}}_{yi} \\ -\vec{F}^{\mathrm{DKT}}_{zi} \end{pmatrix}, \\ \Delta\vec{T}_i&=\begin{pmatrix} \vec{F}^{\mathrm{DKT}}_{\phi_{xi}} \\ \vec{F}^{\mathrm{DKT}}_{\phi_{yi}} \\ 0 \end{pmatrix}. \end{align*}

Note

The value of \(\vec{T}_{zi}\) (drilling torque) is always zero; therefore, drilling motion of nodes is unconstrained (though still governed by dynamics). This can lead to some problems in special cases, manifesting as wobbly rotation of nodes which does not go away.

This video demonstrates the cylindrical triaxial test which includes CST+DKT elements with hydrostatic pressure and interaction with particles inside the membrane (displacements are scaled 10×):

Lumped mass and inertia

Triangle mass and inertia are computed supposing the entire thickness in the triangle plane. The “sleeve” rounding the edges is ignored for the purposes of computing mass and inertia.

Inertia computation is described at Wikipedia and summarized here. For triangle with vertices \(\vec{a}\), \(\vec{b}\), \(\vec{c}\), its double-area is computed as

\[A=\left|(\vec{b}-\vec{a})\times(\vec{c}-\vec{a})\right|.\]

Using covariance of unit triangle

\[\begin{split}\mat{S}=\frac{1}{24}\begin{pmatrix}2&1&1\\1&2&1\\1&1&2\end{pmatrix}\end{split}\]

and transformation matrix

\[\begin{split}\mat{V}=\begin{pmatrix}\vec{a}^T \\ \vec{b}^T \\ \vec{c}^T\end{pmatrix},\end{split}\]

we compute the covariance as

\[\mat{C}=\frac{A}{2}\mat{V}^T\mat{S}\mat{V}.\]

Inertia tensor (with respect to global origin and axes) is then evaluated as usually,

\[\mat{J}=\operatorname{tr}(\mat{C})\mat{I}_3-\mat{C}.\]

When lumping mass and inertia, only the part adjacent to each node is considered according to Voronoi tesselation; triangle is partitioned into 6 sub-triangles. With midpoints \(\vec{m}_{ab}=\frac{\vec{a}+\vec{b}}{2}\) etc., the part adjacent to \(\vec{a}\) are triangles \(\{\vec{a},\vec{m}_{ab},\vec{m}_{ac}\}\) and \(\{\vec{m}_{ab},\vec{C},\vec{m}_{ac}\}\), as shown on the image:

Fig. 29 Partitioning triangle where 2 sub-triangles belong to each node for the purposes of lumping mass and inertia (greenish, reddish and yellowish colors to \(\vec{a}\), \(\vec{b}\), \(\vec{c}\) respectively)

Tip

Report issues or inclarities to github.