Co-rotational kinematics, EICR (pyfe3d.corotational)

The element-independent co-rotational (EICR) formulation turns the linear elements Tria3R, Tria3DSG, Quad4, Quad4R, BeamC, BeamLR, Truss and Spring into elements for large displacements and large rotations with small strains, without changing them: the local linear stiffness matrix of each element is taken once in its base configuration, and the rigid motion of the element is removed from the nodal displacements and rotations at every evaluation. With kinematics='corotational_von_karman' the local response is that of the von Kármán strains of the element instead, evaluated by its own routines at the deformational displacements, see CorotationalShell; this changes the results of BeamC only, since the von Kármán terms of the other elements vanish in their co-rotational frame.

• CorotationalModel assembles a model whose elements are added with kinematics='corotational', 'corotational_von_karman', 'von_karman' or 'linear', and provides the callables calc_fext, calc_fint, calc_kC and calc_kG of pyfe3d.solver.Analysis, which is used unchanged.

• CorotationalShell, CorotationalBeam, CorotationalTruss and CorotationalSpring wrap one element each.

• The finite-rotation utilities spin(), axial(), expso3(), logso3(), Lambda(), dLambdaT_m(), Hmat() and Lmat() are used by the wrappers and by the applied-moment helpers of CorotationalModel.

The module is a Cython extension, as the elements: the kinematics, the projector, and the internal force and tangent stiffness matrix of each wrapper are computed by typed C kernels, of which the finite-rotation utilities are the Python interface. The classes are Python classes, whose attributes can be read and changed as usual.

Module constants:

pyfe3d.corotational.SERIES_ANGLE = 0.2

Below this angle, in radians, the coefficients of the rotation matrices are evaluated with their Taylor series.

pyfe3d.corotational.SERIES_ANGLE_ETA_MU = 0.7

Below this angle, in radians, the coefficients of Hmat() and Lmat() are evaluated with their Taylor series.

pyfe3d.corotational.WARP_RATIO_WARNING = 0.01

Warp ratio of a quadrilateral above which CorotationalShell warns with warping=None, see its attribute warp_ratio.

pyfe3d.corotational.WARP_RATIO_WARNING_OFFSETS = 0.07

Warp ratio of a quadrilateral above which CorotationalShell warns with warping='offsets'.

The rotational degrees-of-freedom of the nodes are total rotation vectors, \(\boldsymbol{\psi}_a\), whose nodal rotation matrix is \(\boldsymbol{R}_a = \exp(\mathrm{Spin}(\boldsymbol{\psi}_a))\), and they are added by the solvers like the translations. No state is kept between the calls. The map is singular at \(|\boldsymbol{\psi}_a| = 2 k \pi\), \(k \ge 1\), which limits a node to less than one full revolution.

Notes

Notation: \(\boldsymbol{T}_0\) and \(\boldsymbol{T}_R\) are the base and the current co-rotational (CR) frames of an element, as 3x3 matrices whose columns are the element axes, the convention of the r11..r33 of the elements, such that global components are \(\boldsymbol{T}\) times local components. Felippa and Haugen [1] write the transposes of these matrices; Nour-Omid and Rankin [2] write \(\boldsymbol{E}_0\) and \(\boldsymbol{E}\).

The spatial spin \(\delta \boldsymbol{\omega}\) of a node, \(\delta \boldsymbol{R} \boldsymbol{R}^T = \mathrm{Spin}(\delta \boldsymbol{\omega})\), is related to the variation of the rotation vector by \(\delta \boldsymbol{\omega} = \boldsymbol{\Lambda}(\boldsymbol{\psi}) \delta \boldsymbol{\psi}\), see Lambda(). This \(\boldsymbol{\Lambda}\) is the \(\boldsymbol{H}(\boldsymbol{\theta})\) of Crisfield [3], Eq. (16.90), and the \(\boldsymbol{H}^{-1}\) of [1], Eq. (100); here \(\boldsymbol{H}\) is always the matrix of [1], Eq. (101), see Hmat().

The internal force and the consistent tangent stiffness matrix of an element are those of the consistent symmetrizable equilibrated (CSE) formulation of [1], Section 6.3, Eq. (53), with the local linear stiffness matrix of the element, written in terms of the spins, and transformed to the additive rotation vectors with \(\boldsymbol{\Lambda}\) as in [3], Section 17.4.2, Eqs. (17.173)-(17.177), see CorotationalShell.

References

[1] (1,2,3,4)

Felippa, C. A. and Haugen, B., A unified formulation of small-strain corotational finite elements: I. Theory, Computer Methods in Applied Mechanics and Engineering, 194(21-24), 2285-2335, 2005, https://doi.org/10.1016/j.cma.2004.07.035.

[2]

Nour-Omid, B. and Rankin, C. C., Finite rotation analysis and consistent linearization using projectors, Computer Methods in Applied Mechanics and Engineering, 93(3), 353-384, 1991, https://doi.org/10.1016/0045-7825(91)90248-5.

[3] (1,2)

Crisfield, M. A., Non-linear Finite Element Analysis of Solids and Structures, Volume 2: Advanced Topics, John Wiley & Sons, Chichester, 1997.

class pyfe3d.corotational.CorotationalBeam(elem, prop, x, q2_min=0.01, exact_KGP=True, local='linear', **elem_kwargs)

Co-rotational (EICR) wrapper of a BeamC or BeamLR element

As CorotationalShell, with the CR frame of the mean nodal triad.

Usually created by CorotationalModel.add().

Parameters:

elemBeamC or BeamLR

The element, with c1, c2 set and update_rotation_matrix(vxyi, vxyj, vxyk, x) called with the same x.

propBeamProp

Property of the element.

xarray_like

Initial nodal coordinates of the model, size 3*num_nodes.

q2_minfloat, optional

Smallest value of \(\bar{q}_2\), see the Notes, below which the CR frame is taken as undefined: the forces and the stiffness are then NaN, such that the non-linear solvers reduce the increment, and a warning is issued. Default 1.e-2.

exact_KGPbool, optional

Whether the projection geometric stiffness includes the variation of \(\boldsymbol{G}\), which depends on the nodal rotations, see the Notes. False gives Felippa and Haugen’s Eq. (47) alone. Default True.

localstr, optional

'linear' or 'von_karman', see CorotationalShell and the Notes. Default 'linear'.

**elem_kwargs

Passed to update_KC0 of the element.

Attributes:

q0np.ndarray

Base \(y\) axis of the element, \(\boldsymbol{q}_0 = \boldsymbol{T}_0 \boldsymbol{e}_2\).

The other attributes are those of :class:`CorotationalShell`, without

``warp_ratio``.

Methods

evaluate(self, u[, tangent])	Internal force vector and tangent stiffness matrix of the element in the additive degrees-of-freedom
kinematics(self, u)	Co-rotational frame and deformational displacements at a configuration
update_probe(self, u)	Write the deformational displacements into the probe of the element, for stress recovery with its own routines

Raises:

TypeError

When elem is not a BeamC or a BeamLR.

ValueError

When r11..r33 of the element are not orthonormal or their \(x\) axis is not along the line from node 1 to node 2 of x, i.e. update_rotation_matrix was not called with x, or when local is not 'linear' or 'von_karman'.

Notes

The CR frame is the mean nodal triad, Crisfield [1], Eqs. (17.11)-(17.14), the “midpoint interpolation” of Nour-Omid and Rankin [2], Section 5:

\[\boldsymbol{e}_1 = \frac{\boldsymbol{x}^D_2 - \boldsymbol{x}^D_1}{l}, \qquad \boldsymbol{q} = \frac{1}{2} (\boldsymbol{R}_1 + \boldsymbol{R}_2) \boldsymbol{q}_0, \qquad \boldsymbol{e}_3 = \frac{\boldsymbol{e}_1 \times \boldsymbol{q}}{|\boldsymbol{e}_1 \times \boldsymbol{q}|}, \qquad \boldsymbol{e}_2 = \boldsymbol{e}_3 \times \boldsymbol{e}_1\]

With \(\bar{\boldsymbol{q}} = \boldsymbol{T}_R^T \boldsymbol{q} = (q_1, q_2, 0)\) and \(\bar{\boldsymbol{q}}_i = \boldsymbol{T}_R^T \boldsymbol{R}_i \boldsymbol{q}_0\), its spin-fitter matrix (derived here) is

\[\omega_1 = \frac{1}{2 q_2} \sum_i (\bar{q}_{i2} \omega_{i,x} - \bar{q}_{i1} \omega_{i,y}) - \frac{q_1}{q_2} \frac{w_2 - w_1}{l}, \qquad \omega_2 = -\frac{w_2 - w_1}{l}, \qquad \omega_3 = \frac{v_2 - v_1}{l}\]

with \(\omega_{i,x}, \omega_{i,y}\) the local spins of the nodes. The frame is undefined when \(\boldsymbol{q}\) becomes parallel to the chord, \(q_2 \to 0\), see q2_min; [2], Appendix B.3, redefine \(\boldsymbol{q}\) in that case, which a stateless formulation cannot do without a jump of the forces.

\(\boldsymbol{G}\) depends on the nodal rotations, so the condition (37) of Felippa and Haugen [3] does not hold and their Eq. (47) misses a term of the order of the unbalanced force ([3], Eq. (48)). With exact_KGP=True the projection geometric stiffness is the exact variation of \(\boldsymbol{P}^T \boldsymbol{p}\) (derived here),

\[\boldsymbol{K}_{GP} = -\frac{\partial (\boldsymbol{G}^T \boldsymbol{z})}{\partial \bar{\boldsymbol{v}}} - \boldsymbol{G}^T \boldsymbol{F}_n(\boldsymbol{p})^T \boldsymbol{P}, \qquad \boldsymbol{z} = \boldsymbol{S}^T \boldsymbol{p}\]

with \(\boldsymbol{p} = \bar{\boldsymbol{H}}^T \boldsymbol{K}_e \boldsymbol{d}\) the unprojected force, \(\delta l = \delta \bar{u}_2 - \delta \bar{u}_1\) and \(\delta \bar{\boldsymbol{q}}_i = \mathrm{Spin}(\bar{\boldsymbol{q}}_i)(\boldsymbol{G} \delta \bar{\boldsymbol{v}} - \delta \bar{\boldsymbol{\omega}}_i)\); it is checked by finite differences in tests/test_corotational_core.py.

With local='von_karman' (see CorotationalShell), the deformational translations of the nodes are along the chord, \(\bar{v} = \bar{w} = 0\). BeamLR interpolates \(v\) and \(w\) linearly, so its von Kármán terms vanish and the result is that of 'linear' to round-off. BeamC interpolates them with the nodal rotations, and its von Kármán axial strain \(\epsilon_{xx} = u_{,x} + \bar{\epsilon}_{NL}\), with \(\bar{\epsilon}_{NL}\) the element mean of \((v_{,x}^2 + w_{,x}^2)/2\) (see pyfe3d.beamc.BeamC.update_KCNL()), gives the local “higher-order terms” of [1], Section 17.1.5: with an infinite transverse shear stiffness, it is the strain of [1], Eqs. (17.56a)-(17.56b), for the initially straight element, cubic \(v\) and \(w\), which is checked in tests/test_corotational_core.py. Measured on the cantilever of Sze, Liu and Lo [4], Section 3.1, at the end force P_max, against the elastica: with 'linear', 8 and 16 BeamC give errors of \(-U\) and \(W\) of -0.07 %, +0.17 % and -0.02 %, +0.05 %; with 'von_karman', 4, 8 and 16 BeamC give -0.013 %, -0.010 %; +0.001 %, +0.006 % and +0.002 %, +0.007 %.

References

[1] (1,2,3)

Crisfield, M. A., Non-linear Finite Element Analysis of Solids and Structures, Volume 2: Advanced Topics, Wiley, 1997.

[2] (1,2)

Nour-Omid, B. and Rankin, C. C., Finite rotation analysis and consistent linearization using projectors, CMAME 93, 353-384, 1991, https://doi.org/10.1016/0045-7825(91)90248-5.

[3] (1,2)

Felippa, C. A. and Haugen, B., A unified formulation of small-strain corotational finite elements: I. Theory, CMAME 194, 2285-2335, 2005, https://doi.org/10.1016/j.cma.2004.07.035.

[4]

Sze, K. Y., Liu, X. H. and Lo, S. H., Popular benchmark problems for geometric nonlinear analysis of shells, Finite Elements in Analysis and Design, 40(11), 1551-1569, 2004, https://doi.org/10.1016/j.finel.2003.11.001.

evaluate(self, u, tangent=True)

Internal force vector and tangent stiffness matrix of the element in the additive degrees-of-freedom

Parameters:

unp.ndarray

Vector of all the degrees-of-freedom of the model, size N.

tangentbool, optional

Whether the tangent stiffness matrix is computed. Default True.

Returns:

fnp.ndarray

Internal force vector of the element, size 6*nn, at the degrees-of-freedom dofs: forces, and moments conjugate to the rotation vectors.

KM, KGnp.ndarray

Only with tangent=True. Material and geometric parts of the tangent stiffness matrix, shape (6*nn, 6*nn), whose sum is the Jacobian of f with respect to u[dofs].

All the arrays are filled with NaN when the CR frame is undefined,

see CorotationalBeam.

kinematics(self, u)

Co-rotational frame and deformational displacements at a configuration

Parameters:

unp.ndarray

Vector of all the degrees-of-freedom of the model, size N: nodal translations and total rotation vectors.

Returns:

TRnp.ndarray or None

Current CR frame \(\boldsymbol{T}_R\), shape (3, 3), columns = element axes; None when the frame is undefined, see CorotationalBeam.

dnp.ndarray or None

Deformational displacements and rotations in the CR frame, \([\bar{\boldsymbol{u}}_{d,a}, \bar{\boldsymbol{\theta}}_{d,a}]\) per node, size 6*nn, in the order of the degrees-of-freedom of the element.

update_probe(self, u)

Write the deformational displacements into the probe of the element, for stress recovery with its own routines

probe.xe receives the base local coordinates \(\tilde{\boldsymbol{x}}_a\) and probe.ue the deformational displacements d of kinematics(), both in element coordinates, such that e.g. elem.update_probe_finte(prop) gives the local forces \(\boldsymbol{K}_e \boldsymbol{d}\).

Parameters:

unp.ndarray

Vector of all the degrees-of-freedom of the model, size N.

Raises:

ValueError

When the CR frame is undefined at u.

class pyfe3d.corotational.CorotationalModel(x, N=None)

Model with co-rotational, von Kármán and linear elements, with the callables of pyfe3d.solver.Analysis

The elements are added with add(), each with its kinematics. The vectors and matrices returned by calc_fext(), calc_fint(), calc_kC() and calc_kG() are those of all the N degrees-of-freedom, to be used with the bu and uk of pyfe3d.solver.Analysis. The rotational degrees-of-freedom are total rotation vectors at the nodes of the co-rotational elements, see pyfe3d.corotational.

Parameters:

xarray_like

Initial nodal coordinates, size 3*num_nodes, as used by the elements.

Nint or None, optional

Number of degrees-of-freedom. None gives 6*num_nodes. Default None.

Attributes:

xnp.ndarray

Copy of the initial nodal coordinates.

Nint

Number of degrees-of-freedom.

fextnp.ndarray

Constant external load vector, size N, used by calc_fext(): forces, and moments conjugate to the rotation vectors, see add_moment(). It may be edited.

corotationallist

The wrappers (CorotationalShell, CorotationalBeam, CorotationalTruss, CorotationalSpring) of the elements added with kinematics='corotational' or 'corotational_von_karman', in order; their attribute local tells the two apart.

von_karman, linearlist

Tuples (elem, prop, elem_kwargs) of the elements added with kinematics='von_karman' and 'linear'.

Methods

add(self, elem, prop[, kinematics])	Add an element to the model
add_moment(self, node, M[, kind])	Apply a concentrated moment at a node
calc_fint(self, c[, inc, silent])	Internal force vector
calc_kC(self[, c, NLgeom, silent])	Constitutive (material) stiffness matrix
calc_kG(self[, c, NLgeom, inc, silent])	Geometric stiffness matrix, including the load stiffness
deformational_displacements(self, c)	Deformational displacements of the co-rotational elements
nodal_rotations(self, c)	Rotation matrices of the nodes

See also

pyfe3d.solver.analysis.Analysis

the non-linear static solvers, which take the callables of the model.

Notes

The residual seen by the solvers is

\[\boldsymbol{R}(\boldsymbol{c}, \lambda) = \boldsymbol{F}_{int}(\boldsymbol{c}) - \lambda \left(\boldsymbol{F}_0 + \sum_k \boldsymbol{F}_k(\boldsymbol{c})\right)\]

with \(\boldsymbol{F}_0\) the attribute fext and \(\boldsymbol{F}_k(\boldsymbol{c})\) the configuration-dependent moments of add_moment(). As described in pyfe3d.solver.callbacks, calc_fint() returns \(\boldsymbol{F}_{int} - \lambda \sum_k (\boldsymbol{F}_k(\boldsymbol{c}) - \boldsymbol{F}_k(\boldsymbol{0}))\), calc_kG() contains the load stiffness \(-\lambda \sum_k \partial \boldsymbol{F}_k/\partial \boldsymbol{c}\), and calc_fext(c=c) returns \(\boldsymbol{F}_0 + \sum_k \boldsymbol{F}_k(\boldsymbol{c}) = -\partial \boldsymbol{R}/\partial \lambda\).

The contributions of the kinematics are:

• 'corotational': CorotationalShell.evaluate() or CorotationalBeam.evaluate(), the material part in calc_kC() and the geometric part in calc_kG();

• 'corotational_von_karman': the same, with the von Kármán strains of the element in its CR frame (local='von_karman' of the wrappers, see CorotationalShell); the local K_e + KCNL(d) enters the material part and the local KG(d) the geometric part;

• 'von_karman': update_fint(..., nonlinear=1) of the element, KC0 + KCNL in calc_kC() and KG in calc_kG();

• 'linear': KC0 u and KC0 in calc_kC().

Limits of validity: 'von_karman' keeps the rotations of the elements small to moderate, the assumption of the von Kármán strains; the two co-rotational kinematics allow large displacements and rotations with small strains, each nodal rotation vector staying below one full revolution, \(|\boldsymbol{\psi}_a| < 2 \pi\), see pyfe3d.corotational. In the CR frame, the von Kármán terms of the shells and of BeamLR vanish, and 'corotational_von_karman' gives the results of 'corotational' for them, to round-off; it changes the results of BeamC only, see CorotationalBeam.

The rotational degrees-of-freedom of a node shared by a co-rotational element and an element of another kinematics are interpreted by the latter as small rotations, valid only while the rotations of that node remain small; a warning is issued at the first evaluation.

property KC0

Linear stiffness matrix of all the elements, sparse, shape (N, N), from their update_KC0

add(self, elem, prop, kinematics='corotational', **elem_kwargs)

Add an element to the model

Parameters:

elemTria3R, Quad4, Quad4R, BeamC, Tria3DSG, BeamLR, Spring or Truss

The element, with its node positions c1, c2, ... set and update_rotation_matrix called with the coordinates x of the model.

propShellProp, BeamProp or None

Property of the element; None for a Spring, whose stiffnesses are attributes of the element.

kinematicsstr, optional

'corotational' (EICR, large rotations, small strains), 'corotational_von_karman' (EICR with the von Kármán strains of the element in its CR frame, see CorotationalShell), 'von_karman' (the non-linear strains of the element, moderate rotations) or 'linear'. Default 'corotational'. Spring and Truss, which have no von Kármán internal force vector of their own, do not accept 'von_karman', and 'corotational_von_karman' is treated as 'corotational' for them, without warning.

**elem_kwargs

Passed to update_KC0 and, with 'von_karman' and 'corotational_von_karman', to update_fint or update_probe_finte of the element, e.g. the hgfactor_* of Quad4R; with 'corotational' and 'corotational_von_karman', the keyword arguments of CorotationalShell (warping), CorotationalBeam (q2_min, exact_KGP), CorotationalTruss (frame, q2_min, exact_KGP) and CorotationalSpring (frame, q2_min) are also accepted.

Returns:

wrapperobject or None

The wrapper of a co-rotational element, a CorotationalShell, CorotationalBeam, CorotationalTruss or CorotationalSpring, also appended to corotational; None for the other kinematics.

Raises:

ValueError

When kinematics is not one of the four, or is not accepted by the element type; with 'corotational' and 'corotational_von_karman', also the errors of the wrapper classes, e.g. when update_rotation_matrix of the element was not called with the coordinates x of the model.

TypeError

When the element type is not supported.

add_moment(self, node, M, kind='rotation_vector')

Apply a concentrated moment at a node

Parameters:

nodeint

Node index, whose rotational degrees-of-freedom are 6*node + 3 to 6*node + 5.

Marray_like

Moment vector, size 3, in global components at the base configuration, for a unit load factor.

kindstr, optional

• 'rotation_vector': added to fext, conjugate to the rotation vector \(\boldsymbol{\psi}\) of the node, constant;

• 'fixed_axes': moment about fixed global axes, configuration-dependent;

• 'follower': moment that rotates with the node, \(\boldsymbol{R}(\boldsymbol{\psi}) \boldsymbol{M}\), configuration-dependent.

Default 'rotation_vector'.

Raises:

ValueError

When kind is not one of the three, or M is not of size 3.

Notes

The virtual work of a spatial moment \(\boldsymbol{m}\) is \(\boldsymbol{m} \cdot \delta \boldsymbol{\omega} = (\boldsymbol{\Lambda}^T \boldsymbol{m}) \cdot \delta \boldsymbol{\psi}\), see Lambda(). Hence:

• 'rotation_vector': \(\boldsymbol{F} = \boldsymbol{M}\); it is the moment about fixed axes when the rotation vector keeps a fixed direction, since \(\boldsymbol{\Lambda}^T \boldsymbol{n} = \boldsymbol{n}\) for \(\boldsymbol{\psi} \parallel \boldsymbol{n}\), e.g. in plane problems;

• 'fixed_axes': \(\boldsymbol{F} = \boldsymbol{\Lambda}(\boldsymbol{\psi})^T \boldsymbol{M}\), with the load stiffness \(-\partial \boldsymbol{F}/\partial \boldsymbol{\psi}\), see dLambdaT_m(). Moments about fixed axes are non-conservative (Crisfield [1], Section 17.5.2, Fig. 17.4), so the tangent stiffness matrix is unsymmetric;

• 'follower': \(\boldsymbol{m} = \boldsymbol{R} \boldsymbol{M}\) and \(\boldsymbol{F} = \boldsymbol{\Lambda}^T \boldsymbol{R} \boldsymbol{M} = \boldsymbol{\Lambda} \boldsymbol{M}\), since \(\boldsymbol{\Lambda}^T \boldsymbol{R} = \boldsymbol{\Lambda}\) (derived here and checked in the tests); the load stiffness is \(-\partial (\boldsymbol{\Lambda} \boldsymbol{M})/\partial \boldsymbol{\psi} = \partial (\boldsymbol{\Lambda}^T \boldsymbol{M})/\partial \boldsymbol{\psi}\) evaluated at \(-\boldsymbol{\psi}\), since \(\boldsymbol{\Lambda}(\boldsymbol{\psi}) = \boldsymbol{\Lambda}(-\boldsymbol{\psi})^T\). This is the moment counterpart of the follower force of [1], Eqs. (17.211)-(17.212), following the nodal triad; the follower moments of [1], Section 17.5.2, Eqs. (17.216)-(17.218), follow the element axes instead.

The configuration-dependent moments enter calc_fint(), calc_kG() and calc_fext() as described in CorotationalModel; use calc_fext() as the calc_fext of pyfe3d.solver.Analysis.

References

[1] (1,2,3)

Crisfield, M. A., Non-linear Finite Element Analysis of Solids and Structures, Volume 2: Advanced Topics, Wiley, 1997.

property calc_fext

External load vector, the calc_fext of pyfe3d.solver.Analysis

A callable calc_fext(inc=1., silent=True, c=None), with:

incfloat, optional

Load factor that multiplies the returned vector, dimensionless. Default 1..

silentbool, optional

Unused, for the signature of pyfe3d.solver.Analysis.

cnp.ndarray or None, optional

Configuration of all the N degrees-of-freedom, for the configuration-dependent moments of add_moment(). None gives the reference vector, with these moments at the base configuration. Default None.

It returns inc*(fext + sum_k F_k(c)), size N, see CorotationalModel.

Notes

pyfe3d.solver.Analysis treats the load as configuration-dependent when calc_fext declares c in its signature, see pyfe3d.solver.callbacks.accepts_kwarg(), and a linear static analysis then iterates on the load. The returned callable declares c only when configuration-dependent moments (kind='fixed_axes' or 'follower' of add_moment()) are present; otherwise c is accepted as a keyword argument and ignored. The attribute is read when it is passed to pyfe3d.solver.Analysis, so add the moments before.

calc_fint(self, c, inc=1., silent=True)

Internal force vector

Parameters:

cnp.ndarray

Configuration of all the N degrees-of-freedom: translations and rotation vectors.

incfloat, optional

Load factor, dimensionless, which scales the configuration-dependent part of the moments of add_moment(). Default 1..

silentbool, optional

Unused, for the signature of pyfe3d.solver.Analysis.

Returns:

fintnp.ndarray

Size N: the internal forces of all the elements, minus inc times the configuration-dependent part of the moments, see CorotationalModel.

calc_kC(self, c=None, NLgeom=False, silent=True)

Constitutive (material) stiffness matrix

Parameters:

cnp.ndarray or None, optional

Configuration of all the N degrees-of-freedom. Default None.

NLgeombool, optional

False or c=None give the linear stiffness matrix KC0. Default False.

silentbool, optional

Unused, for the signature of pyfe3d.solver.Analysis.

Returns:

kCcsr_matrix

Shape (N, N): KC0, or at c the linear KC0 of the linear elements, KC0 + KCNL of the von Kármán elements and the material part of the co-rotational tangent stiffness matrix.

calc_kG(self, c=None, NLgeom=False, inc=1., silent=True)

Geometric stiffness matrix, including the load stiffness

Parameters:

cnp.ndarray or None, optional

Configuration of all the N degrees-of-freedom. None gives a null matrix. Default None.

NLgeombool, optional

False gives a null matrix. Default False.

incfloat, optional

Load factor, dimensionless, which scales the load stiffness of the configuration-dependent moments. Default 1..

silentbool, optional

Unused, for the signature of pyfe3d.solver.Analysis.

Returns:

kGcsr_matrix

Shape (N, N): KG of the von Kármán elements, the geometric part of the co-rotational tangent stiffness matrix and the load stiffness -inc*dF_k/dc of the configuration-dependent moments. With calc_kC() it is the exact Jacobian of calc_fint().

deformational_displacements(self, c)

Deformational displacements of the co-rotational elements

Parameters:

cnp.ndarray

Configuration of all the N degrees-of-freedom.

Returns:

dslist

For each element of corotational, in order, the deformational displacements and rotations in its CR frame, size 6*nn, or None when its frame is undefined, see CorotationalShell.kinematics().

nodal_rotations(self, c)

Rotation matrices of the nodes

Parameters:

cnp.ndarray

Configuration of all the N degrees-of-freedom.

Returns:

Rnp.ndarray

Shape (N//6, 3, 3), R[a] = expso3(c[6*a + 3:6*a + 6]), the rotation of node a from the base configuration, see expso3().

class pyfe3d.corotational.CorotationalShell(elem, prop, x, warping='offsets', local='linear', **elem_kwargs)

Co-rotational (EICR) wrapper of a Tria3R, Tria3DSG, Quad4 or Quad4R element

The local linear stiffness matrix \(\boldsymbol{K}_e\) of the element is computed once, in the base configuration, by its own update_KC0 with the rotation matrix set to the identity; the element is not modified, and the attributes changed during that call are restored. evaluate() gives the internal force vector and the consistent tangent stiffness matrix at any configuration, in the additive degrees-of-freedom of the model: translations and total rotation vectors.

Usually created by CorotationalModel.add().

Parameters:

elemTria3R, Tria3DSG, Quad4 or Quad4R

The element, with its node positions c1, c2, ... set and update_rotation_matrix(x, ...) called with the same x, which also defines its material direction.

propShellProp

Property of the element.

xarray_like

Initial nodal coordinates of the model, size 3*num_nodes, as used by the elements.

warpingstr or None, optional

Treatment of a warped quadrilateral, see the Notes: 'offsets', rigid offsets between the nodes and the flat element of K_e; None, the flat element at the warped nodes, with the projector alone. Ignored for triangles and flat quadrilaterals. Default 'offsets'.

localstr, optional

Strains of the element in its CR frame: 'linear', the local linear stiffness matrix \(\boldsymbol{K}_e\) (pure EICR), or 'von_karman', the von Kármán strains of the element, see the Notes. Default 'linear'.

**elem_kwargs

Passed to update_KC0 of the element, e.g. the hgfactor_* of Quad4R, and with local='von_karman' to its update_probe_finte.

Attributes:

nnint

Number of nodes, 3 or 4.

dofsnp.ndarray

Positions of the 6*nn degrees-of-freedom of the element in the vectors of the model.

T0np.ndarray

Base frame \(\boldsymbol{T}_0\), shape (3, 3), columns = element axes, equal to r11..r33 of the element.

xtildenp.ndarray

Base local coordinates \(\tilde{\boldsymbol{x}}_a = \boldsymbol{T}_0^T (\boldsymbol{x}_a - \boldsymbol{x}_{C0})\), shape (nn, 3), with \(\boldsymbol{x}_{C0}\) the mean of the nodal coordinates.

Kenp.ndarray

Local linear stiffness matrix, shape (6*nn, 6*nn).

warpingstr or None

The warping argument.

localstr

The local argument.

warp_ratiofloat

max|xtilde[:, 2]|/sqrt(A), with \(A\) the area of the flat projection; zero for triangles. A warning is issued above WARP_RATIO_WARNING with warping=None and above WARP_RATIO_WARNING_OFFSETS with warping='offsets'.

Methods

evaluate(self, u[, tangent])	Internal force vector and tangent stiffness matrix of the element in the additive degrees-of-freedom
kinematics(self, u)	Co-rotational frame and deformational displacements at a configuration
update_probe(self, u)	Write the deformational displacements into the probe of the element, for stress recovery with its own routines

Raises:

TypeError

When elem is not a Tria3R, Tria3DSG, Quad4 or Quad4R.

ValueError

When r11..r33 of the element do not match the frame of x, i.e. update_rotation_matrix was not called with x.

ValueError

When warping is not 'offsets' or None, or local is not 'linear' or 'von_karman'.

Notes

Base configuration, Felippa and Haugen [1], Table 3, step 1. At a configuration with the deformed nodal positions \(\boldsymbol{x}^D_a\) and the nodal rotation matrices \(\boldsymbol{R}_a = \exp(\mathrm{Spin}(\boldsymbol{\psi}_a))\), the deformational displacements in the CR frame \(\boldsymbol{T}_R\) are ([1], Section 3.8, Eq. (22), Table 3, steps 2-3; Nour-Omid and Rankin [2], Eqs. (11), (14), (19))

\[\bar{\boldsymbol{u}}_{d,a} = \boldsymbol{T}_R^T (\boldsymbol{x}^D_a - \boldsymbol{x}_C) - \tilde{\boldsymbol{x}}_a, \qquad \bar{\boldsymbol{\theta}}_{d,a} = \mathrm{axial}(\log( \boldsymbol{T}_R^T \boldsymbol{R}_a \boldsymbol{T}_0))\]

with \(\boldsymbol{x}_C\) the mean of \(\boldsymbol{x}^D_a\). The positions \(\boldsymbol{x}^D_a\) are computed relative to the base centroid of the element, \((\boldsymbol{X}_a - \boldsymbol{X}_C) + \boldsymbol{u}_a\), which leaves \(\bar{\boldsymbol{u}}_d\) and \(\boldsymbol{T}_R\) unchanged and makes their round-off independent of the origin of the coordinates: from \(\boldsymbol{X}_a + \boldsymbol{u}_a\) it is of the order of \(\varepsilon |\boldsymbol{X}_a|\), with \(\varepsilon\) the machine precision, and the Newton-Raphson iterations stagnate above the round-off floor of pyfe3d.solver.convergence on fine meshes far from the origin (derived here, measured with test_cantilever_end_shear_translated of tests/test_corotational_benchmarks.py). The CR frame follows the frame rule of the element, applied to \(\boldsymbol{x}^D_a\):

• Tria3R and Tria3DSG, side 1-2 alignment ([1], Appendix B.2.1, Eq. (115); Battini and Pacoste [3], Eqs. (2)-(5)). The frame, and thus the solution, depends on the numbering of the nodes ([1], B.2.1; [3], Section 1).

• Quad4 and Quad4R, \(z\) along \(\boldsymbol{v}_{42} \times \boldsymbol{v}_{13}\) and \(x\) along \((\boldsymbol{v}_{13} + \boldsymbol{v}_{42})/2\), invariant to a cyclic renumbering of the nodes.

The spin-fitter matrix \(\boldsymbol{G}\) of these frames, \(\delta \boldsymbol{\omega}_R = \boldsymbol{G} \delta \bar{\boldsymbol{v}}\), was derived here; with \(\bar{\boldsymbol{x}}_a\) the current local coordinates and \(w, v\) the local translations:

• triangle, \(\bar{\boldsymbol{x}}_2 - \bar{\boldsymbol{x}}_1 = (a_1, 0, 0)\), \(\bar{\boldsymbol{x}}_3 - \bar{\boldsymbol{x}}_1 = (a_2, b_2, 0)\), \(n = a_1 b_2\): \(\omega_1 = -(a_2 (w_2 - w_1) - a_1 (w_3 - w_1))/n\), \(\omega_2 = -(w_2 - w_1)/a_1\), \(\omega_3 = (v_2 - v_1)/a_1\), the last one as in [3], Eq. (5);

• quadrilateral, \(\boldsymbol{d}_1 = \bar{\boldsymbol{x}}_3 - \bar{\boldsymbol{x}}_1 = (a_1, b_1, 0)\), \(\boldsymbol{d}_2 = \bar{\boldsymbol{x}}_2 - \bar{\boldsymbol{x}}_4 = (a_2, b_2, 0)\), \(n_z = a_2 b_1 - a_1 b_2\), \(c = |\boldsymbol{d}_1 + \boldsymbol{d}_2|/2\): \(\omega_1 = -(a_1 (w_2 - w_4) - a_2 (w_3 - w_1))/n_z\), \(\omega_2 = (b_2 (w_3 - w_1) - b_1 (w_2 - w_4))/n_z\), \(\omega_3 = ((v_3 - v_1) + (v_2 - v_4))/(2 c)\).

Both satisfy \(\boldsymbol{G} \boldsymbol{S} = \boldsymbol{I}\) ([1], Eq. (17)) and the condition (37) of [1], such that the projection geometric stiffness of [1], Eq. (47), is exact, which is verified by finite differences in tests/test_corotational_core.py.

The internal force and the tangent stiffness matrix are those of the CSE formulation of [1], Section 6.3, Eq. (53), with the projector \(\boldsymbol{P} = \boldsymbol{P}_u - \boldsymbol{S} \boldsymbol{G}\) ([1], Eqs. (11), (12), (31)-(34)):

\[\bar{\boldsymbol{p}} = \boldsymbol{P}^T \bar{\boldsymbol{H}}^T \boldsymbol{K}_e \boldsymbol{d}, \qquad \boldsymbol{K}_R = \boldsymbol{P}^T \bar{\boldsymbol{H}}^T \boldsymbol{K}_e \bar{\boldsymbol{H}} \boldsymbol{P} + \boldsymbol{P}^T \bar{\boldsymbol{L}} \boldsymbol{P} - \boldsymbol{F}_{nm} \boldsymbol{G} - \boldsymbol{G}^T \boldsymbol{F}_n^T \boldsymbol{P}\]

with \(\bar{\boldsymbol{H}} = \mathrm{diag}(\boldsymbol{I}, \boldsymbol{H}(\bar{\boldsymbol{\theta}}_{d,a}))\) ([1], Eq. (20)), \(\bar{\boldsymbol{L}} = \mathrm{diag}(\boldsymbol{0}, \boldsymbol{L}(\bar{\boldsymbol{\theta}}_{d,a}, \bar{\boldsymbol{m}}_a))\) with the moments of \(\boldsymbol{K}_e \boldsymbol{d}\) ([1], Eqs. (21), (46), (102)), and \(\boldsymbol{F}_{nm}\), \(\boldsymbol{F}_n\) built with \(\bar{\boldsymbol{p}}\) ([1], Eq. (44)). The first term is the material stiffness ([1], Eq. (43)), the others the geometric stiffness ([1], Eqs. (45)-(47)). With \(\boldsymbol{\mathcal{T}} = \mathrm{diag}(\boldsymbol{T}_R)\) and \(\hat{\boldsymbol{\Lambda}} = \mathrm{diag}(\boldsymbol{I}, \boldsymbol{\Lambda}(\boldsymbol{\psi}_a))\), the force and the stiffness in the additive degrees-of-freedom are (Crisfield [4], Section 17.4.2, Eqs. (17.176)-(17.177))

\[\boldsymbol{f} = \hat{\boldsymbol{\Lambda}}^T \boldsymbol{\mathcal{T}} \bar{\boldsymbol{p}}, \qquad \boldsymbol{K} = \hat{\boldsymbol{\Lambda}}^T \boldsymbol{\mathcal{T}} \boldsymbol{K}_R \boldsymbol{\mathcal{T}}^T \hat{\boldsymbol{\Lambda}} + \mathrm{diag}\left(\boldsymbol{0}, \frac{\partial (\boldsymbol{\Lambda}^T \boldsymbol{m}_a)}{\partial \boldsymbol{\psi}_a}\right)\]

with \(\boldsymbol{m}_a\) the moments of \(\boldsymbol{\mathcal{T}} \bar{\boldsymbol{p}}\), see dLambdaT_m(). evaluate returns the material part \(\hat{\boldsymbol{\Lambda}}^T \boldsymbol{\mathcal{T}} \boldsymbol{K}_M \boldsymbol{\mathcal{T}}^T \hat{\boldsymbol{\Lambda}}\) and the geometric part, the rest. The matrix is unsymmetric away from equilibrium ([1], Section 6; [4], Section 17.4.2).

Warped quadrilaterals. \(\boldsymbol{K}_e\) and KC0 use the flat projection of the element (update_KC0 ignores the local \(z\)), whose nodes are at \(\tilde{z}_a\) from the warped nodes along the normal of the mean plane, with \(\tilde{z}_1 = \tilde{z}_3 = -\tilde{z}_2 = -\tilde{z}_4\) for the frame of the diagonals. They do not annihilate the rigid rotations of the warped element, which makes the linear element too stiff in warped meshes, the “warping sensitivity” of [2], Section 1. With warping=None the flat element is evaluated at the warped nodes, and the projector built on the warped geometry restores the invariance to rigid motions ([2], Section 1; [1], Sections 4.2 and 6.5): the tangent at the base configuration is \(\boldsymbol{P}^T \boldsymbol{K}_e \boldsymbol{P}\). With warping='offsets', the default, the flat element is linked to the nodes by rigid offsets that rotate with the nodal triads, the eccentric nodes of [1], Sections 5.1 and 6, whose formulation (their reference [37]) is not given there. The following was derived here. With \(\boldsymbol{e}_a = \tilde{z}_a \boldsymbol{T}_0 \boldsymbol{e}_3\) and \(\boldsymbol{r}_a = \boldsymbol{R}_a \boldsymbol{e}_a\), the flat element has the nodes \(\boldsymbol{x}^D_a - \boldsymbol{r}_a\), the rotations \(\boldsymbol{R}_a\) and the base local coordinates \((\tilde{x}_a, \tilde{y}_a, 0)\), from which \(\boldsymbol{d}\) is measured; its degrees-of-freedom vary as

\[\delta \boldsymbol{u}^f_a = \delta \boldsymbol{u}_a + \mathrm{Spin}(\boldsymbol{r}_a) \boldsymbol{\Lambda}_a \delta \boldsymbol{\psi}_a, \qquad \delta \boldsymbol{\psi}^f_a = \delta \boldsymbol{\psi}_a\]

i.e. \(\delta \boldsymbol{v}^f = \boldsymbol{J} \delta \boldsymbol{v}\). With \(\boldsymbol{f}^f\) and \(\boldsymbol{K}^f\) the force and the tangent of the flat element given above, and \(\boldsymbol{n}_a\) the nodal forces of \(\boldsymbol{f}^f\),

\[\boldsymbol{f} = \boldsymbol{J}^T \boldsymbol{f}^f, \qquad \boldsymbol{K} = \boldsymbol{J}^T \boldsymbol{K}^f \boldsymbol{J} + \mathrm{diag}\left(\boldsymbol{0}, \frac{\partial (\boldsymbol{\Lambda}_a^T (\boldsymbol{n}_a \times \boldsymbol{r}_a))}{\partial \boldsymbol{\psi}_a}\right)\]

The last term, the eccentricity geometric stiffness of [1], Section 5.1, is \(\partial(\boldsymbol{\Lambda}_a^T \boldsymbol{m})/\partial \boldsymbol{\psi}_a\) at \(\boldsymbol{m} = \boldsymbol{n}_a \times \boldsymbol{r}_a\) (dLambdaT_m()) minus \(\boldsymbol{\Lambda}_a^T \mathrm{Spin}(\boldsymbol{n}_a) \mathrm{Spin}(\boldsymbol{r}_a) \boldsymbol{\Lambda}_a\), from \(\delta \boldsymbol{r}_a = -\mathrm{Spin}(\boldsymbol{r}_a) \boldsymbol{\Lambda}_a \delta \boldsymbol{\psi}_a\). The material part is \(\boldsymbol{J}^T \boldsymbol{K}^f_M \boldsymbol{J}\). At the base configuration \(\boldsymbol{K}\) is KC0 transformed by \(\boldsymbol{J}\), the linear flat element with rigid links, which annihilates the rigid motions of the warped element; it differs from KC0 by the order of the warp ratio. Both treatments keep the non-zero eigenvalues of KC0 at the base configuration. The wrapper sets the attribute warp_offsets of a Quad4 or Quad4R to 0, which switches off the linear rigid offsets of the element itself: those of 'offsets' are the exact ones above, and None uses the flat element without offsets. The stress recovery after update_probe() is therefore that of the flat element. KC0 of a CorotationalModel is computed with the element’s own offsets, warp_offsets = 1, for 'offsets', giving \(\boldsymbol{J}^T \boldsymbol{K}_{C0} \boldsymbol{J}\), the tangent of 'offsets' at the base configuration, and with warp_offsets = 0 for None. In the linear analysis of a cantilever twisted by 90 degrees (length 12, width 1.1, thickness 0.32, tip loads), with 6x1 Quad4 elements (warp ratio 0.049), the tip displacements with 'offsets' are within 0.7 % of those of a 96x16 mesh, 8-11 % too small with None and 74-91 % too small with the KC0 of the flat element, warp_offsets = 0; a coarser case is checked in tests/test_corotational_core.py.

Von Kármán local strains. With local='von_karman' the local element computation is step 8 of the procedure of [4], Section 17.1.5: the existing element routines are entered with the local displacements and give the local internal force vector and the local tangent stiffness matrix, which [4] transforms as \(\boldsymbol{K}_t = \boldsymbol{F}^T \boldsymbol{K}_{tl} \boldsymbol{F} + \boldsymbol{K}_{t\sigma}\) (step 11). Here the element’s own update_probe_finte(prop, nonlinear=1), update_KCNL and update_KG are evaluated with probe.ue = d, by the identity trick of \(\boldsymbol{K}_e\), giving \(\bar{\boldsymbol{f}}(\boldsymbol{d})\) and the local tangent \(\boldsymbol{K}_e + \boldsymbol{K}_{CNL}(\boldsymbol{d}) + \boldsymbol{K}_G(\boldsymbol{d})\), the exact Jacobian of \(\bar{\boldsymbol{f}}\) (update_fint of the elements). \(\bar{\boldsymbol{f}}\) replaces \(\boldsymbol{K}_e \boldsymbol{d}\) in \(\bar{\boldsymbol{p}}\), \(\bar{\boldsymbol{L}}\), \(\boldsymbol{F}_n\) and \(\boldsymbol{F}_{nm}\) above, and the local tangent replaces \(\boldsymbol{K}_e\) in the material stiffness. [1] obtain the material stiffness from the variation of the local force, their Eq. (41), with \(\boldsymbol{K}_e = \partial \boldsymbol{p}_e/\partial \boldsymbol{v}_d\) (their Eq. (42)), and note that their Eq. (43) holds only if \(\boldsymbol{K}_e\) does not depend on the deformational freedoms; the chain rule in Eq. (41) then gives \(\boldsymbol{P}^T \bar{\boldsymbol{H}}^T (\partial \bar{\boldsymbol{f}}/\partial \boldsymbol{d}) \bar{\boldsymbol{H}} \boldsymbol{P}\) (derived here), whereas the geometric terms of their Eqs. (44)-(47), which depend on the force and not on its derivative, keep their form. The material part returned by evaluate() takes \(\boldsymbol{K}_e + \boldsymbol{K}_{CNL}\) and the geometric part takes \(\boldsymbol{K}_G\), both transformed as the material stiffness: the split of kinematics='von_karman' in CorotationalModel. KCNL of the elements includes the geometric term of the axial force of the non-linear strains, see their update_KCNL.

In the CR frames of the shells the transverse translations of \(\boldsymbol{d}\) are zero (triangles, whose frame contains the three nodes) or the hourglass pattern \(\pm \bar{z}\) (quadrilaterals, whose \(x y\) plane is parallel to both diagonals and contains the mean of the nodes), the warp of the element and not a deflection, and the von Kármán terms of Tria3R, Tria3DSG, Quad4 and Quad4R take the slopes of \(w\) from the nodal translations only, at the centroid for the quadrilaterals, where the hourglass pattern has no slope. Quad4R takes them there by construction. Quad4 takes them at the two-by-two points with its default shear_model = 1, see the Notes of Quad4.update_KG(), where the pattern has a slope; the wrapper therefore sets the attribute centroid_slopes of a Quad4 to 1 while it calls the element, in the same way as it sets warp_offsets to 0, and restores it afterwards, also when the call raises, such that an element used outside the wrapper keeps its own rule. These terms therefore vanish, and local='von_karman' gives the internal force and the tangent of 'linear' to round-off, for both values of shear_model of Quad4 (derived here, checked in tests/test_corotational_core.py, where the element alone at the same \(\boldsymbol{d}\) with the slopes at the two-by-two points gives von Kármán forces of 1.4 % of the linear ones); the same holds for BeamLR, see CorotationalBeam.

References

[1] (1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19)

Felippa, C. A. and Haugen, B., A unified formulation of small-strain corotational finite elements: I. Theory, CMAME 194, 2285-2335, 2005, https://doi.org/10.1016/j.cma.2004.07.035.

[2] (1,2,3)

Nour-Omid, B. and Rankin, C. C., Finite rotation analysis and consistent linearization using projectors, CMAME 93, 353-384, 1991, https://doi.org/10.1016/0045-7825(91)90248-5.

[3] (1,2,3)

Battini, J.-M. and Pacoste, C., On the choice of local element frame for corotational triangular shell elements, Communications in Numerical Methods in Engineering, 20, 819-825, 2004, https://doi.org/10.1002/cnm.708.

[4] (1,2,3,4)

Crisfield, M. A., Non-linear Finite Element Analysis of Solids and Structures, Volume 2: Advanced Topics, Wiley, 1997.

evaluate(self, u, tangent=True)

Internal force vector and tangent stiffness matrix of the element in the additive degrees-of-freedom

Parameters:

unp.ndarray

Vector of all the degrees-of-freedom of the model, size N.

tangentbool, optional

Whether the tangent stiffness matrix is computed. Default True.

Returns:

fnp.ndarray

Internal force vector of the element, size 6*nn, at the degrees-of-freedom dofs: forces, and moments conjugate to the rotation vectors.

KM, KGnp.ndarray

Only with tangent=True. Material and geometric parts of the tangent stiffness matrix, shape (6*nn, 6*nn), whose sum is the Jacobian of f with respect to u[dofs].

All the arrays are filled with NaN when the CR frame is undefined,

see CorotationalBeam.

kinematics(self, u)

Co-rotational frame and deformational displacements at a configuration

Parameters:

unp.ndarray

Vector of all the degrees-of-freedom of the model, size N: nodal translations and total rotation vectors.

Returns:

TRnp.ndarray or None

Current CR frame \(\boldsymbol{T}_R\), shape (3, 3), columns = element axes; None when the frame is undefined, see CorotationalBeam.

dnp.ndarray or None

Deformational displacements and rotations in the CR frame, \([\bar{\boldsymbol{u}}_{d,a}, \bar{\boldsymbol{\theta}}_{d,a}]\) per node, size 6*nn, in the order of the degrees-of-freedom of the element.

update_probe(self, u)

Write the deformational displacements into the probe of the element, for stress recovery with its own routines

probe.xe receives the base local coordinates \(\tilde{\boldsymbol{x}}_a\) and probe.ue the deformational displacements d of kinematics(), both in element coordinates, such that e.g. elem.update_probe_finte(prop) gives the local forces \(\boldsymbol{K}_e \boldsymbol{d}\).

Parameters:

unp.ndarray

Vector of all the degrees-of-freedom of the model, size N.

Raises:

ValueError

When the CR frame is undefined at u.

class pyfe3d.corotational.CorotationalSpring(elem, prop, x, frame='auto', q2_min=0.01, **elem_kwargs)

Co-rotational wrapper of a Spring element

Usually created by CorotationalModel.add().

Parameters:

elemSpring

The element, with c1, c2 and its stiffnesses set and update_rotation_matrix(xi, xj, xk, vxyi, vxyj, vxyk) called.

propNone

Unused; the stiffnesses are attributes of the element.

xarray_like

Initial nodal coordinates of the model, size 3*num_nodes.

framestr, optional

• 'axis': the frame of CorotationalBeam, \(x\) along the current axis from node 1 to node 2 and \(y\) from the mean nodal triad, for a spring whose \(x\) axis is the line through its nodes;

• 'triad': the mean nodal triad, for any orientation and for coincident nodes, see the Notes;

• 'auto': 'axis' when the nodes are distinct, the \(x\) axis of the element is parallel to the line through them and kye = kze = 0, 'triad' otherwise.

With 'axis' the \(x\) axis of the frame follows the nodes, so there is no transverse relative displacement and kye, kze must be zero.

Default 'auto'.

q2_minfloat, optional

Smallest value of \(|\boldsymbol{a}|\) and \(\bar{b}_2\) ('triad') or of \(\bar{q}_2\) ('axis'), below which the CR frame is taken as undefined, see CorotationalBeam. Default 1.e-2.

**elem_kwargs

Passed to update_KC0 of the element.

Attributes:

framestr

'axis' or 'triad'.

p0, q0np.ndarray

Base \(x\) and \(y\) axes of the element, the first two columns of \(\boldsymbol{T}_0\).

The other attributes are those of :class:`CorotationalBeam`.

Methods

evaluate(self, u[, tangent])	Internal force vector and tangent stiffness matrix of the element in the additive degrees-of-freedom
kinematics(self, u)	Co-rotational frame and deformational displacements at a configuration
update_probe(self, u)	Write the deformational displacements into the probe of the element, for stress recovery with its own routines

Raises:

TypeError

When elem is not a Spring.

ValueError

When frame is not one of the three, when r11..r33 of the element are not orthonormal, or with frame='axis' when the nodes coincide, the \(x\) axis of the element is not along the line through them, or kye or kze is not zero.

Notes

The local stiffness matrix \(\boldsymbol{K}_e\) is that of the element, by the identity trick of CorotationalShell: six uncoupled springs between the relative translations and rotations of the nodes, in the element axes \(\boldsymbol{T}_0\), i.e. the user orientation. For distinct nodes with transverse stiffnesses, \(\boldsymbol{K}_e\) does not annihilate the rigid rotations, since the linear element does not balance the moment of the transverse forces; the projector adds that moment (Felippa and Haugen [1], Sections 4.2 and 6.2), so the tangent at the base configuration differs from KC0.

The 'triad' frame is derived here. With \(\boldsymbol{t}_1, \boldsymbol{t}_2\) the first two columns of \(\boldsymbol{T}_0\) and \(\boldsymbol{R}_i\) the nodal rotations,

\[\boldsymbol{a} = \frac{1}{2} (\boldsymbol{R}_1 + \boldsymbol{R}_2) \boldsymbol{t}_1, \qquad \boldsymbol{b} = \frac{1}{2} (\boldsymbol{R}_1 + \boldsymbol{R}_2) \boldsymbol{t}_2, \qquad \boldsymbol{e}_1 = \frac{\boldsymbol{a}}{|\boldsymbol{a}|}, \qquad \boldsymbol{e}_3 = \frac{\boldsymbol{e}_1 \times \boldsymbol{b}}{|\boldsymbol{e}_1 \times \boldsymbol{b}|}, \qquad \boldsymbol{e}_2 = \boldsymbol{e}_3 \times \boldsymbol{e}_1\]

i.e. the frame of CorotationalBeam with the chord replaced by the mean of the nodal \(x\) axes. It depends on the nodal rotations only, so it is defined for coincident nodes. With the local vectors \(\bar{\boldsymbol{p}}_i = \boldsymbol{T}_R^T \boldsymbol{R}_i \boldsymbol{t}_1\) and \(\bar{\boldsymbol{q}}_i = \boldsymbol{T}_R^T \boldsymbol{R}_i \boldsymbol{t}_2\), their means \(\bar{\boldsymbol{a}} = (A, 0, 0)\) and \(\bar{\boldsymbol{b}} = (B_1, B_2, 0)\), and \(\bar{\boldsymbol{\omega}}_i\) the local spins of the nodes, the spin-fitter matrix follows from the conditions that \(\bar{\boldsymbol{a}}\) stays along \(x\) and \(\bar{\boldsymbol{b}}\) in the \(x y\) plane:

\[\omega_3 = -\frac{1}{2 A} \sum_i \mathrm{Spin}(\bar{\boldsymbol{p}}_i)_2 \bar{\boldsymbol{\omega}}_i, \qquad \omega_2 = \frac{1}{2 A} \sum_i \mathrm{Spin}(\bar{\boldsymbol{p}}_i)_3 \bar{\boldsymbol{\omega}}_i, \qquad \omega_1 = \frac{1}{B_2} \left(B_1 \omega_2 - \frac{1}{2} \sum_i \mathrm{Spin}(\bar{\boldsymbol{q}}_i)_3 \bar{\boldsymbol{\omega}}_i\right)\]

with \(\mathrm{Spin}(\cdot)_k\) the row \(k\) of the matrix. It satisfies \(\boldsymbol{G} \boldsymbol{S} = \boldsymbol{I}\) ([1], Eq. (17)). The projection geometric stiffness is the exact one of CorotationalBeam, with \(\partial(\boldsymbol{G}^T \boldsymbol{z})/\partial \bar{\boldsymbol{v}}\) obtained from \(\delta \bar{\boldsymbol{p}}_i = \mathrm{Spin}(\bar{\boldsymbol{p}}_i) (\boldsymbol{G} \delta \bar{\boldsymbol{v}} - \delta \bar{\boldsymbol{\omega}}_i)\) and the same for \(\bar{\boldsymbol{q}}_i\); the rest is the CSE formulation of CorotationalShell. The deformational rotations of the two nodes are opposite, \(\pm \boldsymbol{\phi}/2\), so the rotational springs act on \(\boldsymbol{\phi}\), close to the relative rotation of the nodes. Zero internal force under rigid motions and the consistency of the tangent with central finite differences are checked in tests/test_corotational_core.py. The frame is undefined when the relative rotation of the nodes approaches 180 degrees (\(A \to 0\) or \(B_2 \to 0\)), see q2_min.

References

[1] (1,2)

Felippa, C. A. and Haugen, B., A unified formulation of small-strain corotational finite elements: I. Theory, CMAME 194, 2285-2335, 2005, https://doi.org/10.1016/j.cma.2004.07.035.

evaluate(self, u, tangent=True)

Internal force vector and tangent stiffness matrix of the element in the additive degrees-of-freedom

Parameters:

unp.ndarray

Vector of all the degrees-of-freedom of the model, size N.

tangentbool, optional

Whether the tangent stiffness matrix is computed. Default True.

Returns:

fnp.ndarray

Internal force vector of the element, size 6*nn, at the degrees-of-freedom dofs: forces, and moments conjugate to the rotation vectors.

KM, KGnp.ndarray

Only with tangent=True. Material and geometric parts of the tangent stiffness matrix, shape (6*nn, 6*nn), whose sum is the Jacobian of f with respect to u[dofs].

All the arrays are filled with NaN when the CR frame is undefined,

see CorotationalBeam.

kinematics(self, u)

Co-rotational frame and deformational displacements at a configuration

Parameters:

unp.ndarray

Vector of all the degrees-of-freedom of the model, size N: nodal translations and total rotation vectors.

Returns:

TRnp.ndarray or None

Current CR frame \(\boldsymbol{T}_R\), shape (3, 3), columns = element axes; None when the frame is undefined, see CorotationalBeam.

dnp.ndarray or None

Deformational displacements and rotations in the CR frame, \([\bar{\boldsymbol{u}}_{d,a}, \bar{\boldsymbol{\theta}}_{d,a}]\) per node, size 6*nn, in the order of the degrees-of-freedom of the element.

update_probe(self, u)

Write the deformational displacements into the probe of the element, for stress recovery with its own routines

probe.xe receives the base local coordinates \(\tilde{\boldsymbol{x}}_a\) and probe.ue the deformational displacements d of kinematics(), both in element coordinates, such that e.g. elem.update_probe_finte(prop) gives the local forces \(\boldsymbol{K}_e \boldsymbol{d}\).

Parameters:

unp.ndarray

Vector of all the degrees-of-freedom of the model, size N.

Raises:

ValueError

When the CR frame is undefined at u.

class pyfe3d.corotational.CorotationalTruss(elem, prop, x, frame='auto', q2_min=0.01, exact_KGP=True, **elem_kwargs)

Co-rotational wrapper of a Truss element

Usually created by CorotationalModel.add().

Parameters:

elemTruss

The element, with c1, c2 set and update_rotation_matrix(x) called with the same x.

propBeamProp

Property of the element.

xarray_like

Initial nodal coordinates of the model, size 3*num_nodes.

framestr, optional

• 'chord': the axial force only, along the current axis of the element; the rotational degrees-of-freedom receive no force and no stiffness, as in the linear element;

• 'mean_triad': the frame of CorotationalBeam, \(x\) along the current axis and \(y\) from the mean nodal triad, which also gives the torsion moment;

• 'auto': 'chord' when the torsional stiffness \(G J\) of the element is zero, 'mean_triad' otherwise.

Default 'auto'.

q2_min, exact_KGPoptional

See CorotationalBeam, used with frame='mean_triad'.

**elem_kwargs

Passed to update_KC0 of the element.

Attributes:

framestr

'chord' or 'mean_triad'. With 'chord', kinematics() gives the base frame turned by the smallest rotation that aligns its \(x\) axis with the current axis, and zero deformational rotations.

l0float

Initial length.

The other attributes are those of :class:`CorotationalBeam`.

Methods

evaluate(self, u[, tangent])	Internal force vector and tangent stiffness matrix, see CorotationalShell.evaluate()
kinematics(self, u)	Co-rotational frame and deformational displacements at a configuration
update_probe(self, u)	Write the deformational displacements into the probe of the element, for stress recovery with its own routines

Raises:

TypeError

When elem is not a Truss.

ValueError

When frame is not one of the three, or is 'chord' with a torsional stiffness.

Notes

With frame='chord', the element is the co-rotational truss of Crisfield [1], Section 3.6, with the engineering strain, Eqs. (3.123)-(3.124), written for the space truss of [1], Section 3.7: with \(l_0\) and \(l\) the initial and current lengths, \(\boldsymbol{e}\) the unit vector from node 1 to node 2 and \(N = (E A/l_0) (l - l_0)\), the internal force and the tangent stiffness matrix are

\[\begin{split}\boldsymbol{f} = N \begin{bmatrix} -\boldsymbol{e} \\ \boldsymbol{e} \end{bmatrix}, \qquad \boldsymbol{K} = \begin{bmatrix} \boldsymbol{K}_t & -\boldsymbol{K}_t \\ -\boldsymbol{K}_t & \boldsymbol{K}_t \end{bmatrix}, \qquad \boldsymbol{K}_t = \frac{E A}{l_0} \boldsymbol{e} \boldsymbol{e}^T + \frac{N}{l} (\boldsymbol{I} - \boldsymbol{e} \boldsymbol{e}^T)\end{split}\]

on the translations. The first term of \(\boldsymbol{K}_t\) is the material part ([1], Eq. (3.131)) and the second one the geometric part ([1], Eqs. (3.132)-(3.138) for the planar truss; the space form was written here and is checked by finite differences in tests/test_corotational_core.py). \(E A/l_0\) is the axial entry of the local stiffness matrix of the element. With frame='mean_triad' the torsion moment requires the nodal triads, and the frame is undefined when the mean of the nodal \(y\) axes becomes parallel to the axis, see q2_min of CorotationalBeam; with the nodal rotations constrained, this happens after a rotation of the axis of about 90 degrees, which 'chord' avoids.

References

[1] (1,2,3,4)

Crisfield, M. A., Non-linear Finite Element Analysis of Solids and Structures, Volume 1: Essentials, Wiley, 1991.

evaluate(self, u, tangent=True)

Internal force vector and tangent stiffness matrix, see CorotationalShell.evaluate()

kinematics(self, u)

Co-rotational frame and deformational displacements at a configuration

Parameters:

unp.ndarray

Vector of all the degrees-of-freedom of the model, size N: nodal translations and total rotation vectors.

Returns:

TRnp.ndarray or None

Current CR frame \(\boldsymbol{T}_R\), shape (3, 3), columns = element axes; None when the frame is undefined, see CorotationalBeam.

dnp.ndarray or None

Deformational displacements and rotations in the CR frame, \([\bar{\boldsymbol{u}}_{d,a}, \bar{\boldsymbol{\theta}}_{d,a}]\) per node, size 6*nn, in the order of the degrees-of-freedom of the element.

update_probe(self, u)

Write the deformational displacements into the probe of the element, for stress recovery with its own routines

probe.xe receives the base local coordinates \(\tilde{\boldsymbol{x}}_a\) and probe.ue the deformational displacements d of kinematics(), both in element coordinates, such that e.g. elem.update_probe_finte(prop) gives the local forces \(\boldsymbol{K}_e \boldsymbol{d}\).

Parameters:

unp.ndarray

Vector of all the degrees-of-freedom of the model, size N.

Raises:

ValueError

When the CR frame is undefined at u.

pyfe3d.corotational.Hmat(theta)

Jacobian of the rotation vector with respect to the spin

Parameters:

thetaarray_like

Rotation vector \(\boldsymbol{\theta}\), size 3, in radians, with \(|\boldsymbol{\theta}| < 2 \pi\); in the EICR a deformational rotation, small.

Returns:

Hnp.ndarray

Shape (3, 3), \(\boldsymbol{H}(\boldsymbol{\theta})\) with \(\delta \boldsymbol{\theta} = \boldsymbol{H} \delta \boldsymbol{\omega}\), the inverse of Lambda().

Notes

Felippa and Haugen [1], Eq. (101); Nour-Omid and Rankin [2], Eqs. (A.14)-(A.15); Crisfield [3], Eq. (16.94), where it is \(\boldsymbol{H}^{-1}\):

\[\boldsymbol{H} = \boldsymbol{I} - \frac{1}{2} \mathrm{Spin}(\boldsymbol{\theta}) + \eta \, \mathrm{Spin}(\boldsymbol{\theta})^2, \qquad \eta = \frac{1 - (\theta/2) \cot(\theta/2)}{\theta^2}\]

with the Taylor series of \(\eta\) for small angles.

References

[1]

Felippa, C. A. and Haugen, B., CMAME 194, 2285-2335, 2005, https://doi.org/10.1016/j.cma.2004.07.035.

[2]

Nour-Omid, B. and Rankin, C. C., CMAME 93, 353-384, 1991, https://doi.org/10.1016/0045-7825(91)90248-5.

[3]

Crisfield, M. A., Non-linear Finite Element Analysis of Solids and Structures, Volume 2: Advanced Topics, Wiley, 1997.

pyfe3d.corotational.Lambda(psi)

Matrix that maps the variation of a rotation vector to the spatial spin

Parameters:

psiarray_like

Rotation vector \(\boldsymbol{\psi}\), size 3, in radians.

Returns:

Lamnp.ndarray

Shape (3, 3), \(\boldsymbol{\Lambda}(\boldsymbol{\psi})\) with \(\delta \boldsymbol{\omega} = \boldsymbol{\Lambda} \delta \boldsymbol{\psi}\), where \(\delta \boldsymbol{R} \boldsymbol{R}^T = \mathrm{Spin}(\delta \boldsymbol{\omega})\) for \(\boldsymbol{R} =\) expso3(psi).

Notes

Crisfield [1], Eq. (16.90), where it is called \(\boldsymbol{H}(\boldsymbol{\theta})\); Felippa and Haugen [2], Eq. (100), last form, where it is \(\boldsymbol{H}^{-1}\):

\[\boldsymbol{\Lambda} = \boldsymbol{I} + \frac{1 - \cos \psi}{\psi^2} \mathrm{Spin}(\boldsymbol{\psi}) + \frac{\psi - \sin \psi}{\psi^3} \mathrm{Spin}(\boldsymbol{\psi})^2\]

It is singular for \(\psi = 2 k \pi\), \(k \ge 1\), and \(\boldsymbol{\Lambda}(\boldsymbol{\psi})^T = \boldsymbol{\Lambda}(-\boldsymbol{\psi})\).

References

[1]

Crisfield, M. A., Non-linear Finite Element Analysis of Solids and Structures, Volume 2: Advanced Topics, Wiley, 1997.

[2]

Felippa, C. A. and Haugen, B., CMAME 194, 2285-2335, 2005, https://doi.org/10.1016/j.cma.2004.07.035.

pyfe3d.corotational.Lmat(theta, m)

Derivative of \(\boldsymbol{H}(\boldsymbol{\theta})^T \boldsymbol{m}\) with respect to the spin, for a fixed vector \(\boldsymbol{m}\)

Parameters:

thetaarray_like

Rotation vector \(\boldsymbol{\theta}\), size 3, in radians.

marray_like

Vector of size 3, the moment of a node in the local frame.

Returns:

Lnp.ndarray

Shape (3, 3), \(\boldsymbol{L} = \partial (\boldsymbol{H}^T \boldsymbol{m})/\partial \boldsymbol{\theta} \, \boldsymbol{H}\).

Notes

Felippa and Haugen [1], Eqs. (102)-(103); Nour-Omid and Rankin [2], Eqs. (A.16)-(A.18):

\[\frac{\partial (\boldsymbol{H}^T \boldsymbol{m})}{\partial \boldsymbol{\theta}} = \eta \left[ (\boldsymbol{\theta} \cdot \boldsymbol{m}) \boldsymbol{I} + \boldsymbol{\theta} \boldsymbol{m}^T - 2 \boldsymbol{m} \boldsymbol{\theta}^T \right] + \mu \, \mathrm{Spin}(\boldsymbol{\theta})^2 \boldsymbol{m} \boldsymbol{\theta}^T - \frac{1}{2} \mathrm{Spin}(\boldsymbol{m}), \qquad \mu = \frac{\theta^2 + 4 \cos \theta + \theta \sin \theta - 4}{4 \theta^4 \sin^2(\theta/2)}\]

with the Taylor series of \(\eta\) and \(\mu\) for small angles.

References

[1]

Felippa, C. A. and Haugen, B., CMAME 194, 2285-2335, 2005, https://doi.org/10.1016/j.cma.2004.07.035.

[2]

Nour-Omid, B. and Rankin, C. C., CMAME 93, 353-384, 1991, https://doi.org/10.1016/0045-7825(91)90248-5.

pyfe3d.corotational.axial(W)

Axial vector of a skew-symmetric 3x3 matrix, the inverse of spin()

Parameters:

Wnp.ndarray

Shape (3, 3); only its skew-symmetric part is used.

Returns:

vnp.ndarray

Vector of size 3, (W[2, 1] - W[1, 2])/2 etc.

pyfe3d.corotational.dLambdaT_m(psi, m)

Derivative of \(\boldsymbol{\Lambda}(\boldsymbol{\psi})^T \boldsymbol{m}\) with respect to the rotation vector, for a fixed vector \(\boldsymbol{m}\)

Parameters:

psiarray_like

Rotation vector \(\boldsymbol{\psi}\), size 3, in radians.

marray_like

Vector of size 3, e.g. a moment.

Returns:

Knp.ndarray

Shape (3, 3), \(\partial (\boldsymbol{\Lambda}^T \boldsymbol{m})/\partial \boldsymbol{\psi}\).

Notes

Derived here from Lambda() and checked by finite differences in tests/test_corotational_core.py. With \(a = (1 - \cos \psi)/\psi^2\), \(b = (\psi - \sin \psi)/\psi^3\) and \(\boldsymbol{\Lambda}^T \boldsymbol{m} = \boldsymbol{m} - a \boldsymbol{\psi} \times \boldsymbol{m} + b \boldsymbol{\psi} \times (\boldsymbol{\psi} \times \boldsymbol{m})\):

\[\frac{\partial (\boldsymbol{\Lambda}^T \boldsymbol{m})}{\partial \boldsymbol{\psi}} = a \mathrm{Spin}(\boldsymbol{m}) - \frac{a'}{\psi} (\boldsymbol{\psi} \times \boldsymbol{m}) \boldsymbol{\psi}^T + b \left[ \boldsymbol{\psi} \boldsymbol{m}^T + (\boldsymbol{\psi} \cdot \boldsymbol{m}) \boldsymbol{I} - 2 \boldsymbol{m} \boldsymbol{\psi}^T \right] + \frac{b'}{\psi} (\boldsymbol{\psi} \times (\boldsymbol{\psi} \times \boldsymbol{m})) \boldsymbol{\psi}^T\]

with the Taylor series of \(a'/\psi\) and \(b'/\psi\) for small angles. At \(\boldsymbol{\psi} = \boldsymbol{0}\) it is \(\mathrm{Spin}(\boldsymbol{m})/2\), in agreement with Crisfield [1], Eq. (16.96). This is the extra term of the tangent stiffness matrix in additive rotation vectors, [1], Eq. (17.177).

References

[1] (1,2)

Crisfield, M. A., Non-linear Finite Element Analysis of Solids and Structures, Volume 2: Advanced Topics, Wiley, 1997.

pyfe3d.corotational.expso3(psi)

Rotation matrix of a rotation vector, the exponential map

Parameters:

psiarray_like

Rotation vector \(\boldsymbol{\psi}\), size 3, in radians.

Returns:

Rnp.ndarray

Shape (3, 3), the proper orthogonal matrix of the rotation by the angle \(|\boldsymbol{\psi}|\) about the axis of \(\boldsymbol{\psi}\).

Notes

Rodrigues’ formula, Crisfield [1], Eqs. (16.22) and (16.31); Felippa and Haugen [2], Eq. (89):

\[\boldsymbol{R} = \boldsymbol{I} + \frac{\sin \psi}{\psi} \mathrm{Spin}(\boldsymbol{\psi}) + \frac{1 - \cos \psi}{\psi^2} \mathrm{Spin}(\boldsymbol{\psi})^2, \qquad \psi = |\boldsymbol{\psi}|\]

with the Taylor series of the coefficients for small angles.

References

[1]

Crisfield, M. A., Non-linear Finite Element Analysis of Solids and Structures, Volume 2: Advanced Topics, Wiley, 1997.

[2]

Felippa, C. A. and Haugen, B., CMAME 194, 2285-2335, 2005, https://doi.org/10.1016/j.cma.2004.07.035.

pyfe3d.corotational.logso3(R)

Rotation vector of a rotation matrix, the logarithmic map

Parameters:

Rnp.ndarray

Proper orthogonal matrix, shape (3, 3).

Returns:

thetanp.ndarray

Rotation vector of size 3 with \(|\boldsymbol{\theta}| \le \pi\), such that expso3(theta) is R.

Notes

The quaternion of R is extracted with Spurrier’s algorithm, Felippa and Haugen [1], Eqs. (96)-(97), and Crisfield [2], Section 16.10, Eqs. (16.74)-(16.78), which picks the largest of the trace and the diagonal entries and is accurate for every angle, including 0 and \(\pi\); the closed form \(\sin \theta = |\mathrm{axial}(\boldsymbol{R})|\) of [1], Eq. (98), is not used since it is inaccurate near 0 and \(\pi\). The quaternion is taken with \(q_0 \ge 0\), and the rotation vector is \(\boldsymbol{\theta} = 2 \arctan(|\boldsymbol{q}|/q_0) \boldsymbol{q}/|\boldsymbol{q}|\).

References

[1] (1,2)

Felippa, C. A. and Haugen, B., CMAME 194, 2285-2335, 2005, https://doi.org/10.1016/j.cma.2004.07.035.

[2]

Crisfield, M. A., Non-linear Finite Element Analysis of Solids and Structures, Volume 2: Advanced Topics, Wiley, 1997.

pyfe3d.corotational.spin(v)

Skew-symmetric (spin) matrix of a 3-vector

Parameters:

varray_like

Vector of size 3.

Returns:

Snp.ndarray

Shape (3, 3), such that S @ w == np.cross(v, w).