BeamC - Consistent Timoshenko 3D beam element (`pyfe3d.beamc`)#

Defines the two-node beam element BeamC, with six degrees-of-freedom per node, its probe BeamCProbe and BeamCData, which gives the number of sparse entries per element of each structural matrix. The cross-section and material properties are given by a BeamProp object.

class pyfe3d.beamc.BeamC#

Timoshenko 3D beam element with consistent shape functions

Nodal connectivity for the beam element:

^ y axis
|
|
______   --> x axis
1    2

The constructor stores p in probe, sets eid, pid, n1, n2, c1 and c2 to -1, and sets the init_k_* positions, length, vxyi, vxyj, vxyk and the rotation matrix to 0. The methods use the attributes updated by the others: update_rotation_matrix() reads c1 and c2; update_probe_xe() reads the rotation matrix and updates length; the update_K*, update_M() and update_probe_finte() methods read length and the rotation matrix; update_KG(), update_KCNL() and update_probe_finte() also read probe.ue, updated by update_probe_ue().

Parameters:

pBeamCProbe object: Probe of the element, stored in probe. It can be shared with other elements.

Attributes:

eid,int: Element identification number.
pid,int: Property identification number.
length,double: Element length, dimension length. Updated by update_length().
r11, r12, r13, r21, r22, r23, r31, r32, r33double: Entries \(T_{ij}\) of the transformation matrix \(\boldsymbol{T}_{3\times 3}\) from element to global coordinates, whose columns are the element axes. Dimensionless.
vxyi, vxyj, vxykdouble: Components of a vector on the \(XY\) plane of the element coordinate system, defined using global coordinates. Only the direction of the vector is used.
c1, c2int: Position of each node in the global stiffness matrix, i.e. of the first of its 6 degrees-of-freedom. The coordinates of the node are read from the position c1//2 or c2//2 of the array of nodal coordinates.
n1, n2int: Node identification number.
init_k_KC0, init_k_KCNL, init_k_KG, init_k_Mint: Position in the arrays storing the sparse data for the structural matrices, where the entries of this element start, see BeamCData for their number.
probe,BeamCProbe object: Pointer to the probe.

Methods

`update_KC0`(self, long[, long[, double[, ...)	Update sparse vectors for linear constitutive stiffness matrix KC0
`update_KCNL`(self, long[, long[, double[, ...)	Update sparse vectors for the nonlinear stiffness matrix KCNL
`update_KG`(self, long[, long[, double[, ...)	Update sparse vectors for geometric stiffness matrix KG
`update_M`(self, long[, long[, double[, ...)	Update sparse vectors for mass matrix M
`update_fint`(self, double[, BeamProp prop, ...)	Update the internal force vector
`update_length`(self)	Update element length
`update_probe_finte`(self, BeamProp prop, ...)	Update the internal force vector of the probe
`update_probe_ue`(self, double[)	Update the local displacement vector of the probe of the element
`update_probe_xe`(self, double[)	Update the 3D coordinates of the probe of the element
`update_rotation_matrix`(self, double vxyi, ...)	Update the rotation matrix of the element

Notes

Formulation based on reference [1].

The cross section is rigid. In the element coordinate system, whose \(x\) axis is the beam axis, the displacements of the point \((y, z)\) of the cross section are

\[u_x = u - y r_z + z r_y, \qquad u_y = v - z r_x, \qquad u_z = w + y r_x\]

where \(u, v, w, r_x, r_y, r_z\) are the translations and rotations of the beam axis, interpolated from the 12 element displacements ue. The axial strain, with the curvatures \(\kappa_y\) and \(\kappa_z\), and the transverse shear strains of the beam axis are

\[\epsilon_{xx} = u_{,x} + y \kappa_y + z \kappa_z, \qquad \kappa_y = -{r_z}_{,x}, \qquad \kappa_z = {r_y}_{,x}, \qquad \gamma_{xy} = v_{,x} - r_z, \qquad \gamma_{xz} = w_{,x} + r_y\]

and, with the twist \(\kappa_x = {r_x}_{,x}\), the linear strain energy is

\[U = \frac{1}{2} \int_0^L \boldsymbol{\epsilon}^\top \boldsymbol{C} \boldsymbol{\epsilon} \, dx, \qquad \boldsymbol{\epsilon} = \{u_{,x}, \gamma_{xy}, \gamma_{xz}, \kappa_x, \kappa_y, \kappa_z\}^\top\]

where \(L\) is length and \(\boldsymbol{C}\) is the sectional stiffness matrix C of BeamProp, of entries \(C_{ij}\), \(i, j = 0, \ldots, 5\), in the order of \(\boldsymbol{\epsilon}\), symmetric, with 21 independent entries; see BeamProp for its relation to the stiffness matrices of pyfe3d.vabs. When C is not set, it is built from the cross-section integrals of BeamProp, and

\[\begin{split}U = \frac{1}{2} \int_0^L \Big[ E \left( A u_{,x}^2 + I_{zz} \kappa_y^2 + I_{yy} \kappa_z^2 + 2 A_y u_{,x} \kappa_y + 2 A_z u_{,x} \kappa_z + 2 I_{yz} \kappa_y \kappa_z \right) \\ + G_{12} \left( A \gamma_{xy}^2 - 2 A_z \gamma_{xy} {r_x}_{,x} \right) + G_{13} \left( A \gamma_{xz}^2 + 2 A_y \gamma_{xz} {r_x}_{,x} \right) + G J {r_x}_{,x}^2 \Big] dx\end{split}\]

where \(G_{12}\) = G12 and \(G_{13}\) = G13 are the transverse shear moduli of the planes \(xy\) and \(xz\), and \(G J\) = GJ is the torsional stiffness, see BeamProp. The torsion is of Saint-Venant type, without warping, with GJ in place of the torsional stiffness \(G_{12} I_{yy} + G_{13} I_{zz}\) of the rigid cross section, and no shear correction factor is applied. The axial force is the first entry of \(\boldsymbol{C} \boldsymbol{\epsilon}\), \(N = C_{00} u_{,x} + C_{01} \gamma_{xy} + C_{02} \gamma_{xz} + C_{03} \kappa_x + C_{04} \kappa_y + C_{05} \kappa_z\), i.e. \(N = E (A u_{,x} + A_y \kappa_y + A_z \kappa_z)\) without C.

When the only non-zero off-diagonal entry of \(\boldsymbol{C}\) is the extension-twist coupling \(C_{03}\), e.g. without C and with Ay, Az and Iyz equal to zero, \(u\) and \(r_x\) are interpolated linearly between the nodes, and \(v\) and \(r_z\) are interpolated with the consistent shape functions of [1], cubic for \(v\) and quadratic for \(r_z\), which depend on the parameters

\[\alpha_y = -\frac{12 C_{44}}{C_{11} L^2}, \qquad \alpha_z = -\frac{12 C_{55}}{C_{22} L^2}, \qquad \beta_y = \frac{1}{1 - \alpha_y}, \qquad \beta_z = \frac{1}{1 - \alpha_z}\]

named alphay, alphaz, betay and betaz in the code, taken from the diagonal of \(\boldsymbol{C}\); without C, \(\alpha_y = -12 E I_{zz}/(G_{12} A L^2)\) and \(\alpha_z = -12 E I_{yy}/(G_{13} A L^2)\). With \(\xi = x/L\) and the nodal values \(v_1, {r_z}_1, v_2, {r_z}_2\):

\[\begin{split}v = \left[1 + \beta_y (2\xi^3 - 3\xi^2 + \alpha_y \xi)\right] v_1 + L \beta_y \left(\xi^3 + \frac{\alpha_y - 4}{2} \xi^2 + \left(1 - \frac{\alpha_y}{2}\right) \xi \right) {r_z}_1 \\ + \beta_y (-2\xi^3 + 3\xi^2 - \alpha_y \xi) v_2 + L \beta_y \left(\xi^3 - \frac{\alpha_y + 2}{2} \xi^2 + \frac{\alpha_y}{2} \xi \right) {r_z}_2\end{split}\]

\[r_z = \frac{6 \beta_y}{L} (\xi^2 - \xi) (v_1 - v_2) + \left[1 + \beta_y (3\xi^2 + (\alpha_y - 4) \xi)\right] {r_z}_1 + \beta_y (3\xi^2 - (\alpha_y + 2) \xi) {r_z}_2\]

\(w\) and \(r_y\) use the same functions, with \(v\) replaced by \(w\), \(r_z\) by \(-r_y\), \(\alpha_y\) by \(\alpha_z\) and \(\beta_y\) by \(\beta_z\). The transverse shear strains are then constant along the element, e.g. \(\gamma_{xy} = \alpha_y \beta_y \left((v_1 - v_2)/L + ({r_z}_1 + {r_z}_2)/2\right)\), and they vanish with \(\alpha_y = 0\), i.e. with an infinite G12, where the interpolation of \(v\) becomes the cubic Hermite interpolation of the Euler-Bernoulli beam.

The parameters \(\alpha_y\) and \(\alpha_z\) are negative. With \(\phi_y = 12 C_{44}/(C_{11} L^2)\), \(\beta_y = 1/(1 + \phi_y)\) and the functions above are the general solution of the homogeneous equilibrium equations of the strain energy above with a diagonal \(\boldsymbol{C}\), \((C_{11} \gamma_{xy})_{,x} = 0\) and \(C_{11} \gamma_{xy} = C_{44} {\kappa_y}_{,x}\), i.e. \(G_{12} A \gamma_{xy} = E I_{zz} {\kappa_y}_{,x}\) (and likewise for \(w\), \(r_y\), with \(C_{22}\) and \(C_{55}\)), which these functions satisfy for any nodal values only with this sign. The same solution is Eq. (17) of [1] and Eq. (3.236) of [2], whose denominator \(1 + 3\beta\), with \(\beta = E I/(\kappa A G a^2)\), where \(\kappa G\) is \(G_{12}\), and \(a = L/2\), is \(1 + \phi_y\); Eq. (21) of [1] gives \(\alpha_y\) and \(\alpha_z\) with the positive sign, inconsistent with its Eqs. (17) and (20), and its Eq. (20) prints the functions of \(w\) and \(r_y\) with the signs of those of \(v\) and \(r_z\), which, with \(\alpha_z = 0\), give \(r_y = w_{,x}\) instead of the \(r_y = -w_{,x}\) of the kinematics above. \(\beta_y\), \(\beta_z\) are finite for any length. The mass and geometric stiffness matrices use the same shape functions.

Shape functions with couplings (derived here). The functions above solve the equilibrium equations for any nodal values only when the only coupling is \(C_{03}\), with which \(u_{,x}\) and \(\kappa_x\) are constant. With any other coupling, including those of Ay, Az and Iyz when C is not set, the shape functions are derived from the same equations with the full \(\boldsymbol{C}\). For a prismatic element loaded only at its nodes, the virtual work of the strain energy above gives the equilibrium \(N_{,x} = {V_y}_{,x} = {V_z}_{,x} = {M_x}_{,x} = 0\), \({M_4}_{,x} = V_y\), \({M_5}_{,x} = V_z\) of the sectional forces \(\boldsymbol{C} \boldsymbol{\epsilon} = \{N, V_y, V_z, M_x, M_4, M_5\}^\top\), hence

\[\boldsymbol{C} \boldsymbol{\epsilon} = (\boldsymbol{I} + x \boldsymbol{Z}) \boldsymbol{q}, \qquad \boldsymbol{\epsilon} = \boldsymbol{S} (\boldsymbol{I} + x \boldsymbol{Z}) \boldsymbol{q}, \qquad \boldsymbol{q} = \{N, V_y, V_z, M_x, M_4(0), M_5(0)\}^\top\]

with \(\boldsymbol{S} = \boldsymbol{C}^{-1}\) and the only non-zero entries \(Z_{41} = Z_{52} = 1\). Integrating the kinematic relations from node 1, \(u = \int \epsilon_0 dx\), \(r_x = \int \epsilon_3 dx\), \(r_z = -\int \epsilon_4 dx\), \(r_y = \int \epsilon_5 dx\), \(v = \int (\epsilon_1 + r_z) dx\) and \(w = \int (\epsilon_2 - r_y) dx\), and adding the rigid-body motion of node 1 gives \(\boldsymbol{d} = \boldsymbol{R}(x) \boldsymbol{u}_1 + \boldsymbol{P}(x) \boldsymbol{q}\), with \(\boldsymbol{d} = \{u, v, w, r_x, r_y, r_z\}^\top\), quadratic in \(x\) for \(u, r_x, r_y, r_z\) and cubic for \(v, w\), with coefficients that depend on all entries of \(\boldsymbol{S}\). The statics of the element clamped at node 1 give \(\boldsymbol{q} = \boldsymbol{T} \boldsymbol{P}_2\) for the loads \(\boldsymbol{P}_2\) at node 2, and \(\boldsymbol{F}_c = \boldsymbol{P}(L) \boldsymbol{T}\) is the flexibility of this cantilever, symmetric and positive definite. With \(\boldsymbol{u}_2 = \boldsymbol{R}(L) \boldsymbol{u}_1 + \boldsymbol{P}(L) \boldsymbol{q}\):

\[\boldsymbol{q} = \boldsymbol{H}_f \boldsymbol{u}_e, \qquad \boldsymbol{H}_f = \boldsymbol{T} \boldsymbol{F}_c^{-1} \boldsymbol{G}, \qquad \boldsymbol{d} = \left( [\boldsymbol{R}(x), \boldsymbol{0}] + \boldsymbol{P}(x) \boldsymbol{H}_f \right) \boldsymbol{u}_e, \qquad \boldsymbol{K}_{C0_e} = \boldsymbol{G}^\top \boldsymbol{F}_c^{-1} \boldsymbol{G}\]

where \(\boldsymbol{G} = [-\boldsymbol{R}(L), \boldsymbol{I}]\). \(\boldsymbol{K}_{C0_e}\) is the Hessian of the strain energy with these shape functions, the axial force \(N = q_0\) is uniform along the element, and the mass, geometric and nonlinear matrices are integrated analytically with the same polynomial shape functions. With \(C_{03}\) as the only coupling they are the functions of [1] above, asserted symbolically in theory/deriving_equations/beamc/derive_beamc_coupled_shape_functions.py, whose closed-form expressions are then used, unchanged. Otherwise \(\boldsymbol{S}\), \(\boldsymbol{F}_c^{-1}\), \(\boldsymbol{H}_f\) and the coefficients of the shape functions are evaluated numerically for each element and call, from C, which must then be positive definite, otherwise a ValueError is raised. This evaluation reproduces the arrays of the closed-form expressions to 1.3e-14 of their largest entries, see tests/test_beam_sectional_stiffness.py, and costs about 4 to 7 µs more per element and call than the 1.5 to 2 µs of update_KC0, update_M, update_KG and update_fint with the closed-form expressions (measured on 2026-10-31, Windows, calls from Python).

The element is therefore exact under nodal loads for any \(\boldsymbol{C}\): one element gives the tip displacements of a cantilever under any tip load, measured in tests/test_beam_sectional_stiffness.py to 2.7e-15 of the largest for each single coupling \(C_{ij} = 0.3 \sqrt{C_{ii} C_{jj}}\) of a rectangular section, 50 mm by 30 mm, 0.5 m long, to 2.0e-15 for a \(\boldsymbol{C}\) with 21 non-zero entries and to 5.4e-16 for a composite box beam with extension-twist and bending-shear couplings, against up to 8.9e-2 and 0.42 with the functions of [1] taken from the diagonal of \(\boldsymbol{C}\). BeamLR, whose interpolation does not depend on \(\boldsymbol{C}\), converges as the square of the element length, e.g. 8.3e-2, 5.2e-3 and 3.3e-4 with 1, 4 and 16 elements for the coupling \(C_{15}\) under a transverse tip load.

The sign is verified against the literature in tests/test_timoshenko_beam_literature.py: one element reproduces the cantilever tip deflection \(P L^3/(3 E I_{zz}) + P L/(G A)\) of Eq. (1.2) of [3] (0.053 cm) and, for a cantilever as long as deep, the same formula, where the shear deflection is larger than the bending one; one element reproduces the frequency coefficients 3.322 and 26.569 of Example 3.10 of [2], which has the same shape functions and consistent mass matrix with rotary inertia, three elements the errors of Table 3.6 of [2], and refined meshes its exact frequencies. With the positive sign these tests fail, e.g. by -44 % for the deep cantilever and with 3.952 instead of 3.322. In tests/test_beam_shear_moduli.py, with different G12 and G13, the same tip deflection and frequencies are obtained in each plane with its own modulus, \(G = G_{12}\) with \(I_{zz}\) for \(v\) and \(G = G_{13}\) with \(I_{yy}\) for \(w\).

With nonlinear=1 in update_fint(), the geometric nonlinearity is that of the von Kármán axial strain, with its nonlinear part \((v_{,x}^2 + w_{,x}^2)/2\) replaced by its element mean, Crisfield [4], Section 7.1.3, Eqs. (7.31)-(7.32), and Crisfield [5], Section 17.1.5, Eqs. (17.56a)-(17.56b):

\[\epsilon_{xx} = u_{,x} + \bar{\epsilon}_{NL} + y \kappa_y + z \kappa_z, \qquad \bar{\epsilon}_{NL} = \frac{1}{L} \int_0^L \frac{1}{2} \left( v_{,x}^2 + w_{,x}^2 \right) dx\]

since the constant \(u_{,x}\) of the linear \(u\) cannot balance the pointwise value, of degree four in \(x\) with the cubic \(v\) and \(w\), which gives membrane locking [4]. The tangent stiffness matrix KC0 + KCNL + KG is the exact Jacobian of the internal forces, see update_KCNL() and update_KG(), whose axial force is the element mean of \(N\).

The matrices and the internal forces are computed in element coordinates and transformed to global coordinates with the transformation matrix of the element, e.g. \(\boldsymbol{K} = \boldsymbol{T} \boldsymbol{K}_e \boldsymbol{T}^\top\), where \(\boldsymbol{T}\) is block diagonal with four copies of the matrix \(\boldsymbol{T}_{3\times 3}\) of the entries \(T_{ij}\) stored in r11, ..., r33, whose columns are the element axes, one for the translations and one for the rotations of each node. The element coordinate system is not updated with the displacements.

References

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

Luo, Y., 2008, “An Efficient 3D Timoshenko Beam Element with Consistent Shape Functions,” Adv. Theor. Appl. Mech., 1(3), pp. 95–106.

[2] (1,2,3)

Petyt, M., 2010, Introduction to Finite Element Vibration Analysis, 2nd ed., Cambridge University Press, Section 3.9. https://doi.org/10.1017/CBO9780511761195

[3]

Bathe, K. J., 2014, Finite Element Procedures, 2nd ed., K. J. Bathe, Watertown, MA, ISBN 978-0-9790049-5-7, Section 1.2, p. 5, Eq. (1.2), Fig. 1.2(b).

[4] (1,2)

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

[5]

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

c1#: c1: ‘int’

c2#: c2: ‘int’

eid#: eid: ‘int’

init_k_KC0#: init_k_KC0: ‘int’

init_k_KCNL#: init_k_KCNL: ‘int’

init_k_KG#: init_k_KG: ‘int’

init_k_M#: init_k_M: ‘int’

length#: length: ‘double’

n1#: n1: ‘int’

n2#: n2: ‘int’

pid#: pid: ‘int’

probe#: probe: pyfe3d.beamc.BeamCProbe

r11#: r11: ‘double’

r12#: r12: ‘double’

r13#: r13: ‘double’

r21#: r21: ‘double’

r22#: r22: ‘double’

r23#: r23: ‘double’

r31#: r31: ‘double’

r32#: r32: ‘double’

r33#: r33: ‘double’

update_KC0(self, long[: :1] KC0r, long[: :1] KC0c, double[: :1] KC0v, BeamProp prop, int update_KC0v_only=0) → void#

Update sparse vectors for linear constitutive stiffness matrix KC0

The 144 entries of the matrix of the element, in global coordinates, are written at the positions init_k_KC0 to init_k_KC0 + 143: the row and column indices, built from c1 and c2, are stored in KC0r and KC0c, and the values are added to those already in KC0v. The matrix is computed from length, the rotation matrix and the sectional stiffness matrix C of prop, see BeamProp. length and the rotation matrix must be updated before, see update_rotation_matrix() and update_probe_xe().

Parameters:

KC0rnp.array: Array to store row positions of sparse values. 1-D array of integers of the C type long, of size at least init_k_KC0 + 144.
KC0cnp.array: Array to store column positions of sparse values. Same type and size as KC0r.
KC0vnp.array: Array to store sparse values. 1-D array of floats of size at least init_k_KC0 + 144.
propBeamProp object: Beam property object from where the stiffness and mass attributes are read from.
update_KC0v_onlyint, optional: The default 0 means that the row and column indices KC0r and KC0c should also be updated. Any other value will only update the stiffness matrix values KC0v.

Notes

KC0 is the Hessian, with respect to the element displacements ue, of the linear strain energy given in the notes of BeamC, integrated analytically along the element with the consistent shape functions, and transformed to global coordinates as described there. With couplings in C other than C[0, 3] it is \(\boldsymbol{G}^\top \boldsymbol{F}_c^{-1} \boldsymbol{G}\), from the flexibility of the element clamped at node 1, evaluated numerically, see “Shape functions with couplings” in the notes of BeamC; C must then be positive definite, otherwise a ValueError is raised.

update_KCNL(self, long[: :1] KCNLr, long[: :1] KCNLc, double[: :1] KCNLv, BeamProp prop, int update_KCNLv_only=0) → void#

Update sparse vectors for the nonlinear stiffness matrix KCNL

The matrix, evaluated at the displacements probe.ue, is computed from length, the rotation matrix and the attributes of prop read by _update_probe_KCNLve(). The tangent stiffness matrix of the internal forces of update_fint() with nonlinear=1 is KC0 + KCNL + KG, see update_KC0() and update_KG().

Before this function is called, the probe BeamCProbe attribute of the BeamC object must be updated using update_probe_ue() with the current displacements; and update_probe_xe() with the node coordinates.

The 144 entries of the matrix of the element, in global coordinates, are written at the positions init_k_KCNL to init_k_KCNL + 143: the row and column indices, built from c1 and c2, are stored in KCNLr and KCNLc, and the values are added to those already in KCNLv. The matrix in element coordinates is also stored in probe.KCNLve, overwriting it.

Parameters:

KCNLrnp.array: Array to store row positions of sparse values. 1-D array of integers of the C type long, of size at least init_k_KCNL + 144.
KCNLcnp.array: Array to store column positions of sparse values. Same type and size as KCNLr.
KCNLvnp.array: Array to store sparse values. 1-D array of floats of size at least init_k_KCNL + 144.
propBeamProp object: Beam property object from where the stiffness and mass attributes are read from.
update_KCNLv_onlyint, optional: The default 0 means that the row and column indices KCNLr and KCNLc should also be updated. Any other value will only update the stiffness matrix values KCNLv.

Notes

KCNL = KC0L + KCL0 + KCLL + KGNL, that is \(\boldsymbol{K}_{CNL} + \boldsymbol{K}_{GNL}\), the nonlinear constitutive stiffness matrix \(\boldsymbol{K}_{CNL}\) of the terms of the nonlinear axial strain plus the geometric \(\boldsymbol{K}_{GNL}\) of the axial force of the nonlinear axial strain. Crisfield [1], Section 7.1.1, Eq. (7.24), keeps that axial force in the initial-stress matrix instead, that is in \(\boldsymbol{K}_G + \boldsymbol{K}_{GNL}\); the tangent stiffness matrix is the same. The axial strain is that of von Kármán with its nonlinear part replaced by the element mean, Crisfield [1], Section 7.1.3, Eqs. (7.31)-(7.32), and Crisfield [2], Section 17.1.5, Eqs. (17.56a)-(17.56b):

\[\epsilon_{xx} = u_{,x} + \bar{\epsilon}_{NL}, \qquad \bar{\epsilon}_{NL} = \frac{1}{L} \int_0^L \frac{1}{2} \left( v_{,x}^2 + w_{,x}^2 \right) dx = \frac{1}{2} \boldsymbol{u}_e^\top \boldsymbol{H}_m \boldsymbol{u}_e\]

with \(\boldsymbol{H}_m\), the gradient \(\bar{\boldsymbol{B}}_{NL} = \boldsymbol{u}_e^\top \boldsymbol{H}_m\) and the row \(\overline{\boldsymbol{EB}}_L\) of the mean axial force of the linear strains given in _update_means(). With \(\bar{N}_{NL} = C_{00} \bar{\epsilon}_{NL}\), \(C_{00}\) = C[0, 0] of BeamProp, \(E A\) when C is not set (own derivation, asserted in theory/deriving_equations/beamc/derive_beamc_KC0L_KCL0_KCLL_KGLL_sparse.py):

KC0L = \(L \, \overline{\boldsymbol{EB}}_L^\top \bar{\boldsymbol{B}}_{NL}\)
KCL0 = KC0L`^top`
KCLL = \(L \, C_{00} \, \bar{\boldsymbol{B}}_{NL}^\top \bar{\boldsymbol{B}}_{NL}\)
KGNL = \(L \, \bar{N}_{NL} \boldsymbol{H}_m\)

The first three groups are the constitutive terms coupling the linear and the nonlinear parts of the axial strain. KGNL is geometric, carrying the axial force of the nonlinear axial strain. It is collected here so that update_KG() stays homogeneous of degree one in the displacements, which is what a linear buckling analysis needs. With it here,

\[\boldsymbol{K}_T = \boldsymbol{K}_{C0} + \boldsymbol{K}_{CNL}(\boldsymbol{u}) + \boldsymbol{K}_{GNL}(\boldsymbol{u}) + \boldsymbol{K}_G(\boldsymbol{u})\]

is the exact, symmetric Jacobian of the internal forces of update_fint() with nonlinear=1, and a Newton-Raphson iteration built on them converges quadratically.

With \(u\) linear and \(v\), \(w\) cubic, the pointwise \((v_{,x}^2 + w_{,x}^2)/2\) is a polynomial of degree four in \(x\) that the constant \(u_{,x}\) cannot balance, so that the axial strain cannot vanish in inextensional bending and coarse meshes are too stiff (membrane locking, [1], Section 7.1.3). The element mean removes it; [1] justifies it with the Hu-Washizu principle and a constant membrane strain. With zero transverse nodal displacements, as in the local frame of a co-rotational element, and an infinite transverse shear stiffness, \(\bar{\epsilon}_{NL}\) is the term \(\boldsymbol{\theta}_l^\top \boldsymbol{X} \boldsymbol{\theta}_l/60\) of [2], Eq. (17.56a). The slopes \(v_{,x}\) and \(w_{,x}\) include the transverse shear strains, so that \(\alpha_y\), \(\beta_y\), \(\alpha_z\) and \(\beta_z\) enter \(\boldsymbol{H}_m\), or, with couplings in C other than C[0, 3], all entries of \(\boldsymbol{C}^{-1}\), see _update_means(). The von Kármán axial strain is valid for small strains and moderate rotations: the nonlinear axial strain is uniform over the cross section and along the element, the curvatures and the transverse shear strains are linear in the displacements, and the element coordinate system is not updated with the displacements.

References

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

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

[2] (1,2)

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

update_KG(self, long[: :1] KGr, long[: :1] KGc, double[: :1] KGv, BeamProp prop, int update_KGv_only=0) → void#

Update sparse vectors for geometric stiffness matrix KG

The 144 entries of the matrix of the element, in global coordinates, are written at the positions init_k_KG to init_k_KG + 143: the row and column indices, built from c1 and c2, are stored in KGr and KGc, and the values are added to those already in KGv. The matrix depends on the axial force computed from the displacements probe.ue, from length and from the sectional stiffness matrix C of prop, see BeamProp: its entries C[0, 0], C[0, 3], C[1, 1], C[2, 2], C[4, 4] and C[5, 5] without other couplings, all of them otherwise. probe.ue, length and the rotation matrix must be updated before, see update_probe_ue() and update_probe_xe().

Parameters:

KGrnp.array: Array to store row positions of sparse values. 1-D array of integers of the C type long, of size at least init_k_KG + 144.
KGcnp.array: Array to store column positions of sparse values. Same type and size as KGr.
KGvnp.array: Array to store sparse values. 1-D array of floats of size at least init_k_KG + 144.
propBeamProp object: Beam property object from where the stiffness and mass attributes are read from.
update_KGv_onlyint, optional: The default 0 means that the row and column indices KGr and KGc are also updated. Any other value will only update the stiffness matrix values KGv.

Notes

The geometric stiffness matrix is

\[\boldsymbol{K}_G = N_m \int_0^L \left[ \left(\frac{\partial \boldsymbol{S}^v}{\partial x}\right)^\top \frac{\partial \boldsymbol{S}^v}{\partial x} + \left(\frac{\partial \boldsymbol{S}^w}{\partial x}\right)^\top \frac{\partial \boldsymbol{S}^w}{\partial x} \right] dx = L N_m \boldsymbol{H}_m\]

the Hessian of \(L N_m \bar{\epsilon}_{NL}\) at a constant \(N_m\), where \(\partial \boldsymbol{S}^v / \partial x\) and \(\partial \boldsymbol{S}^w / \partial x\) give the slopes \(v_{,x}\) and \(w_{,x}\) with the consistent shape functions, \(\bar{\epsilon}_{NL}\) and \(\boldsymbol{H}_m\) are those of update_KCNL() and

\[N_m = \frac{1}{L} \int_0^L N dx\]

is the element mean of the axial force of the linear strains, \(N = \boldsymbol{C}_{0,:} \boldsymbol{\epsilon}\), at probe.ue, with the first row of the sectional stiffness matrix, see the notes of BeamC. Without couplings in C other than C[0, 3], \(N_m = (C_{00} (u_2 - u_1) + C_{03} ({r_x}_2 - {r_x}_1))/L\), i.e. \(E A (u_2 - u_1)/L\) when C is not set and Ay, Az and Iyz are zero; with other couplings, \(N_m = q_0\), the first entry of \(\boldsymbol{q} = \boldsymbol{H}_f \boldsymbol{u}_e\), see “Shape functions with couplings” in the notes of BeamC. The mean follows from the mean nonlinear axial strain of update_KCNL(), with which the energy of the axial strain contains \(\int_0^L N dx \, \bar{\epsilon}_{NL}\), such that KC0 + KCNL + KG is the exact Jacobian of the internal forces of update_fint() with nonlinear=1 (own derivation, asserted in theory/deriving_equations/beamc/derive_beamc_KC0L_KCL0_KCLL_KGLL_sparse.py); the geometric matrix of Crisfield [1], Eq. (7.95), also carries an axial force constant along the element. \(N_m\) is the axial force of the nodal internal forces, -finte[0] = finte[6]. Since the shape functions solve the equilibrium equations, \(N\) is uniform along the element for any C and \(N_m = N\).

\(\boldsymbol{K}_G\) is homogeneous of degree one in probe.ue. With probe.ue from the linear static solution of a reference load, the critical load multipliers \(\lambda\) of the linear buckling analysis solve \((\boldsymbol{K}_{C0} + \lambda \boldsymbol{K}_G) \boldsymbol{u} = \boldsymbol{0}\), see linear_buckling(). The geometric term of the nonlinear axial strain is in update_KCNL().

References

[1]

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

update_M(self, long[: :1] Mr, long[: :1] Mc, double[: :1] Mv, BeamProp prop, int mtype=0) → void#

Update sparse vectors for mass matrix M

The entries of the matrix of the element, in global coordinates, are written from the position init_k_M: 144 entries with mtype=0 and 36 entries with mtype=1. The row and column indices, built from c1 and c2, are stored in Mr and Mc, and the values are added to those already in Mv. The matrix is computed from length, the rotation matrix and the attributes intrho, intrhoy, intrhoz, intrhoy2, intrhoz2 and intrhoyz, and the diagonal entries C[1, 1], C[2, 2], C[4, 4] and C[5, 5] of the sectional stiffness matrix of prop, or all its entries with couplings other than C[0, 3] and mtype=0. length and the rotation matrix must be updated before, see update_rotation_matrix() and update_probe_xe().

Parameters:

Mrnp.array: Array to store row positions of sparse values. 1-D array of integers of the C type long, of size at least init_k_M + 144 (init_k_M + 36 with mtype=1).
Mcnp.array: Array to store column positions of sparse values. Same type and size as Mr.
Mvnp.array: Array to store sparse values. 1-D array of floats, of the same size as Mr.
propBeamProp object: Beam property object from where the inertial attributes are read.
mtypeint, optional: 0 for consistent mass matrix, 1 for lumped mass matrix. Default 0. The element has no reduced-integration mass matrix, 2.

Raises:

ValueError: If mtype is not 0 or 1, if C is not positive semi-definite, see BeamProp, or if, with mtype=0, C has couplings other than C[0, 3] and is not positive definite.

Notes

The mass matrices come from the kinetic energy of the displacement field of the notes of BeamC, integrated over the cross section with the inertial attributes of BeamProp:

\[T = \frac{1}{2} \int_0^L \Big[ m (\dot{u}^2 + \dot{v}^2 + \dot{w}^2) + m_{yy} \dot{r}_z^2 + m_{zz} \dot{r}_y^2 - 2 m_{yz} \dot{r}_y \dot{r}_z - 2 m_y \dot{u} \dot{r}_z + 2 m_z \dot{u} \dot{r}_y - 2 m_z \dot{v} \dot{r}_x + 2 m_y \dot{w} \dot{r}_x + (m_{yy} + m_{zz}) \dot{r}_x^2 \Big] dx\]

with \(m\) = intrho, \(m_y\) = intrhoy, \(m_z\) = intrhoz, \(m_{yy}\) = intrhoy2, \(m_{zz}\) = intrhoz2 and \(m_{yz}\) = intrhoyz.

mtype=0 gives the consistent mass matrix, the Hessian of \(T\) with respect to the nodal velocities, interpolated with the consistent shape functions of BeamC, which depend on alphay and alphaz and hence on C[1, 1], C[2, 2], C[4, 4] and C[5, 5], i.e. on E, G12, G13, A, Iyy and Izz when C is not set. With couplings in C other than C[0, 3], the shape functions depend on all entries of C, see “Shape functions with couplings” in the notes of BeamC, and C must be positive definite. The integral along the element is computed analytically (exactly), and the rotary inertia is included.

mtype=1 gives the lumped mass matrix, diagonal in element coordinates: at each node, \(L m/2\) for the three translations, \(L (m_{yy} + m_{zz})/2\) for \(r_x\), \(L m_{zz}/2\) for \(r_y\) and \(L m_{yy}/2\) for \(r_z\). It corresponds to the two-point Gauss-Lobatto quadrature of \(T\), whose points are the nodes, without the coupling terms of \(m_y\), \(m_z\) and \(m_{yz}\), which are not included.

The lumped mass matrix follows the lumped mass matrix by Lobatto integration of Brockman (1987) [1], p. 2348, adapted to the beam: the two-point Gauss-Lobatto quadrature of \(T\), whose points are the nodes, gives a diagonal matrix only without the first moments of the density, which are therefore dropped. The consistent mass matrix is the classical one, see Petyt (2010) [2]: Sec. 3.3, Eq. (3.59), for the axial terms, Sec. 3.4, Eq. (3.102), with \(m_{yy} + m_{zz}\) in place of \(\rho I_x\), for the torsional terms, and Sec. 3.9, Eqs. (3.234) and (3.237), for the bending terms with the shape functions of static equilibrium and the rotary inertia, which BeamC reproduces in each plane of bending when \(m_y\), \(m_z\) and \(m_{yz}\) are zero.

References

[1]

Brockman, R. A., 1987, “Dynamics of the Bilinear Mindlin Plate Element,” Int. J. Numer. Methods Eng., 24(12), pp. 2343–2356. https://onlinelibrary.wiley.com/doi/pdf/10.1002/nme.1620241208

[2]

Petyt, M., 2010, Introduction to Finite Element Vibration Analysis, 2nd ed., Cambridge University Press, Cambridge. https://doi.org/10.1017/CBO9780511761195

update_fint(self, double[: :1] fint, BeamProp prop, int nonlinear=0) → void#

Update the internal force vector

Calls update_probe_finte(prop, nonlinear), transforms the forces and moments of probe.finte to global coordinates with the rotation matrix and adds them to fint at the positions c1 to c1 + 5 and c2 to c2 + 5. The displacements probe.ue must be updated before, see update_probe_ue().

Parameters:

fintnp.array: Array that is updated in place with the internal forces. The internal forces stored in fint are calculated in global coordinates. Method update_probe_finte() is called to update the parameter finte of the BeamCProbe with the internal forces in local coordinates. 1-D array of size 6*M, for \(M\) nodes in the model.
propBeamProp object: Beam property object from where the stiffness and mass attributes are read from.
nonlinearint, optional: The default 0 gives the linear internal forces, KC0*u. Any other value adds the geometrically nonlinear terms, for which the tangent stiffness matrix of the internal forces is KC0 + KCNL + KG, see update_KCNL().

Notes

The nonlinear terms are those of the von Kármán strains, for which the exact Jacobian of the internal forces is KC0 + KCNL + KG, see update_probe_finte() and update_KCNL().

update_length(self) → void#

Update element length

The distance between the two nodes, computed from the coordinates probe.xe in the element coordinate system, is stored in the attribute length. Called by update_probe_xe().

update_probe_finte(self, BeamProp prop, int nonlinear=0) → void#

Update the internal force vector of the probe

The attribute finte is updated with the BeamCProbe the internal forces in local coordinates. While using this function, mind that the probe can be shared amongst more than one finite element, depending how you defined them, meaning that the probe will always safe the values from the last udpate.

Note

The finte attribute of object BeamCProbe is updated, accessible using .probe.finte.

The 12 values of finte are overwritten. They are computed from length and from the displacements probe.ue, which must be updated before, see update_probe_xe() and update_probe_ue().

Parameters:

propBeamProp object: Beam property object from where the stiffness and mass attributes are read from.
nonlinearint, optional: The default 0 gives the linear internal forces, KC0*u. Any other value adds the geometrically nonlinear terms, for which the tangent stiffness matrix of the internal forces is KC0 + KCNL + KG, see update_KCNL().

Notes

The linear internal forces are \(\boldsymbol{K}_{C0_e} \boldsymbol{u}_e\), with the matrix \(\boldsymbol{K}_{C0_e}\) of update_KC0() in element coordinates. With nonlinear different from 0, the terms of the von Kármán strains are added, see _update_probe_finte_nonlinear(), such that finte is the gradient of the strain energy with the von Kármán axial strain averaged over the element, see update_KCNL(), for which the exact Jacobian of the internal forces is KC0 + KCNL + KG, see update_KCNL().

update_probe_ue(self, double[: :1] u) → void#

Update the local displacement vector of the probe of the element

Note

The ue attribute of object BeamCProbe is updated, accessible using .probe.ue.

The translations and rotations of the two nodes, read from u at the positions c1 to c1 + 5 and c2 to c2 + 5, are transformed to the element coordinate system with the transpose \(\boldsymbol{T}_{3\times 3}^\top\) of the transformation matrix, which must be updated before, see update_rotation_matrix().

Parameters:

uarray-like: Array with global displacements, for a total of \(M\) nodes in the model, this array will be arranged as: \(u_1, v_1, w_1, {r_x}_1, {r_y}_1, {r_z}_1, u_2, v_2, w_2, {r_x}_2, {r_y}_2, {r_z}_2, ..., u_M, v_M, w_M, {r_x}_M, {r_y}_M, {r_z}_M\). Translations of dimension length and rotations, 1-D array of size 6*M.

update_probe_xe(self, double[: :1] x) → void#

Update the 3D coordinates of the probe of the element

Note

The xe attribute of object BeamCProbe is updated, accessible using .probe.xe.

The coordinates of the two nodes, read from x at the positions c1//2 and c2//2, are transformed to the element coordinate system with the transpose \(\boldsymbol{T}_{3\times 3}^\top\) of the transformation matrix, which must be updated before, see update_rotation_matrix(). Then update_length() is called, updating length.

Parameters:

xarray-like: Array with global nodal coordinates, for a total of \(M\) nodes in the model, this array will be arranged as: \(x_1, y_1, z_1, x_2, y_2, z_2, ..., x_M, y_M, z_M\). Dimension length, 1-D array of size 3*M.

update_rotation_matrix(self, double vxyi, double vxyj, double vxyk, double[: :1] x) → void#

Update the rotation matrix of the element

Attributes r11,r12,r13,r21,r22,r23,r31,r32,r33 are updated, with the entries \(T_{ij}\) of the transformation matrix \(\boldsymbol{T}_{3\times 3}\) from element to global coordinates, whose columns are the element axes.

The element attributes vxyi, vxyj and vxyk are also updated when this function is called.

The element coordinate system is determined, identifying the \(ijk\) components of each axis: \({x_e}_i, {x_e}_j, {x_e}_k\); \({y_e}_i, {y_e}_j, {y_e}_k\); \({z_e}_i, {z_e}_j, {z_e}_k\).

The element \(x\) axis goes from node 1 to node 2, whose coordinates are read from x at the positions c1//2 and c2//2. The \(z\) axis is normal to the \(x\) axis and to the vector (vxyi, vxyj, vxyk), the \(y\) axis completes the right-handed system, and the three axes are normalized.

Parameters:

vxyi, vxyj, vxykdouble: Components of a vector on the \(XY\) plane of the element coordinate system. Only the direction of the vector is used. A vector parallel to the element axis gives a null \(z\) axis and a rotation matrix with nan entries.
xarray-like: Array with global nodal coordinates, for a total of \(M\) nodes in the model, this array will be arranged as: \(x_1, y_1, z_1, x_2, y_2, z_2, ..., x_M, y_M, z_M\). Dimension length, 1-D array of size 3*M.

vxyi#: vxyi: ‘double’

vxyj#: vxyj: ‘double’

vxyk#: vxyk: ‘double’

class pyfe3d.beamc.BeamCData#

Used to allocate memory for the sparse matrices.

Each attribute is the number of entries that one BeamC element writes in the arrays of rows, columns and values of a sparse matrix, at the consecutive positions that start at its init_k_KC0, init_k_KCNL, init_k_KG or init_k_M. The constructor takes no argument and sets the values below.

Attributes:

KC0_SPARSE_SIZE,int: KC0_SPARSE_SIZE = 144, written by BeamC.update_KC0().
KCNL_SPARSE_SIZE,int: KCNL_SPARSE_SIZE = 144, written by BeamC.update_KCNL().
KG_SPARSE_SIZE,int: KG_SPARSE_SIZE = 144, written by BeamC.update_KG().
M_SPARSE_SIZE,int: M_SPARSE_SIZE = 144, written by BeamC.update_M() with mtype=0; with mtype=1 only the first 36 are written.

KC0_SPARSE_SIZE#: KC0_SPARSE_SIZE: ‘int’

KCNL_SPARSE_SIZE#: KCNL_SPARSE_SIZE: ‘int’

KG_SPARSE_SIZE#: KG_SPARSE_SIZE: ‘int’

M_SPARSE_SIZE#: M_SPARSE_SIZE: ‘int’

class pyfe3d.beamc.BeamCProbe#

Probe used for local coordinates, local displacements, local stresses etc

The idea behind using a probe is to avoid allocating larger memory buffers per finite element. The memory buffers are allocated per probe, and one probe can be shared amongst many finite elements, with the information being updated and retrieved on demand.

Note

The probe can be shared amongst more than one finite element, depending on how you defined them. Mind that the probe will always safe the values from the last udpate.

The constructor takes no argument and allocates the arrays below, filled with zeros.

Attributes:

xe,array-like: Array of size NUM_NODES*DOF//2=6 containing the nodal coordinates in the element coordinate system, in the following order \({x_e}_1, {y_e}_1, {z_e}_1, {x_e}_2, {y_e}_2, {z_e}_2\). Dimension length. Updated by BeamC.update_probe_xe().
ue,array-like: Array of size NUM_NODES*DOF=12 containing the element displacements in the following order \(u_{e_1}, v_{e_1}, w_{e_1}, {r_x}_{e_1}, {r_y}_{e_1}, {r_z}_{e_1}, u_{e_2}, v_{e_2}, w_{e_2}, {r_x}_{e_2}, {r_y}_{e_2}, {r_z}_{e_2}\). Translations of dimension length and rotations. Updated by BeamC.update_probe_ue().
finte,array-like: Array of size NUM_NODES*DOF=12 containing the element internal forces corresponding to the degrees-of-freedom described by ue, forces and moments. Updated by BeamC.update_probe_finte().
KCNLvearray-like: KCNLve: ‘double[::1]’

KCNLve#: KCNLve: ‘double[::1]’

finte#: finte: ‘double[::1]’

ue#: ue: ‘double[::1]’

xe#: xe: ‘double[::1]’

BeamC - Consistent Timoshenko 3D beam element (pyfe3d.beamc)#

BeamC - Consistent Timoshenko 3D beam element (`pyfe3d.beamc`)#