Tria3R - Triangular element with reduced integration (`pyfe3d.tria3r`)#

Triangular element with reduced integration, where a single point at the centroid (\(N1=N2=N3=1/3\)) and weight \(weight=1\) is evaluated, preventing shear locking. The drilling stiffness is evaluated with full integration.

The transverse shear stiffnesses \(A_{44}\), \(A_{45}\) and \(A_{55}\) are stabilised using the method of Lyly, Stenberg and Vihinen, often called Stenberg’s method:

Lyly, M., Stenberg, R., & Vihinen, T. (1993). A stable bilinear element for the Reissner-Mindlin plate model. Computer Methods in Applied Mechanics and Engineering, 110(3-4), 343-357. https://doi.org/10.1016/0045-7825(93)90214-I

Bischoff, M., & Bletzinger, K.-U. (2004). Improving stability and accuracy of Reissner-Mindlin plate finite elements via algebraic subgrid scale stabilization. Computer Methods in Applied Mechanics and Engineering, 193(15-16), 1517-1528. https://doi.org/10.1016/j.cma.2003.12.036

in the form adopted by Castro et al.:

Castro, S. G. P., Donadon, M. V., & Guimaraes, T. A. M. (2019). ES-PIM applied to buckling of variable angle tow laminates. Composite Structures, 209, 67-78. https://doi.org/10.1016/j.compstruct.2018.10.058

where the transverse shear terms are changed as per Eq. (26) in Castro et al., here repeated for convenience:

\[\left[\begin{matrix}\hat{A}_{44}, \hat{A}_{45} \ \hat{A}_{45}, \hat{A}_{55}\end{matrix}\right] = \frac{1}{1 + factor}\left[\begin{matrix}A_{44}, A_{45} \ A_{45}, A_{55}\end{matrix}\right]\]

with \(factor\) defined as:

factor = frac{alpha ell^2}{h^2}

where \(\alpha\) is a positive constant parameters (see the alpha_shear_locking attribute), \(\ell\) is the longest edge of the corresponding triangle, and \(h\) the total thickness of the element.

Note that \(factor \rightarrow 0\) as the element size shrinks, so the true stiffness is recovered under mesh refinement and the scheme is consistent. The rate matters though: \(\ell < 0.12 h\) is needed for the stabilised stiffness to be within 1 per cent of the true one, which for a thin shell is not a mesh anyone would build. What makes the scheme legitimate all the same is that in a load-driven problem the error it introduces is \(factor \times f_s\), with \(f_s\) the shear fraction of the response, and for a plate or beam of span \(L\) discretised with \(n\) elements

\[factor \times f_s \sim \frac{\alpha \ell^2}{h^2} c \frac{h^2}{L^2} = \frac{\alpha c}{n^2}\]

so the thickness cancels and the error behaves as an ordinary \(O(1/n^2)\) discretisation error. This is measured in tests/test_transverse_shear_stiffness.py.

Warning

The default \(\alpha = 0.7\) is roughly seven times the value used in all three references above, which is near \(0.1\): Eq. (24) of Bischoff and Bletzinger gives \(\alpha = 0.17 (\ell_g/\ell_n)/(2(1-\nu))\), i.e. \(0.12\) for a square element with \(\nu = 0.3\), and Castro et al. investigate \(0.05\) to \(0.15\) and recommend \(0.07\) to \(0.09\), reporting that above \(0.09\) the linear buckling behaviour becomes overly soft and converges from below.

The difference is not a retune but a consequence of where the factor is applied. In all three references it sits on top of a discrete shear gap (DSG) formulation, which is already free of shear locking, so there the factor is a mild stabilisation that improves coarse-mesh accuracy. This element has no DSG: the transverse shear is taken from a single point at the centroid, which still locks, and the factor is therefore doing the unlocking. Reducing \(\alpha\) to the literature value makes this element stiffen sharply. On the plate of test_tria3r_natural_freq.py, with the consistent mass matrix, the error in the first natural frequency goes

\[+0.6\% \ (\alpha = 0.7) \quad +30.6\% \ (\alpha = 0.1) \quad +77.3\% \ (\alpha = 0.01) \quad +96.1\% \ (\alpha = 0)\]

so with the stabilisation removed the element is nearly twice as stiff as it should be, which is the locking itself. For scale, pyfe3d.tria3dsg.Tria3DSG on the identical mesh, with no parameter of any kind, gives \(+2.8\%\).

The price of the large \(\alpha\) is paid on the transverse shear itself. A field with \(w\) varying and all rotations zero has identically zero curvature, so it is resisted by transverse shear alone, and the stabilisation divides that resistance by \(1 + factor\). Measured on an 8 by 8 patch, this element’s stiffness in that subspace is 0.33, 0.019 and 0.0012 of the Quad4 and Quad4R value at \(\ell/h\) of 1.25, 6.25 and 25. A Donnell geometric stiffness matrix works on \(\partial w/\partial x\), which is exactly that subspace, so thin-shell buckling is where the scheme shows its cost, as measured for a cylinder in tests/test_quad4_linear_buckling_cylinder_displ.py. The failure there is not a uniformly wrong load but a change of critical mode: with \(factor\) at 1852 the element admits a mode at two elements per wavelength, the mesh Nyquist limit, carrying 36.6 per cent of its energy in transverse shear, and that mesh artefact undercuts the physical one. At the literature \(\alpha = 0.1\) the physical mode is critical again and this element agrees with the quadrilateral and with Tria3DSG to 3 per cent on it, while that same \(\alpha = 0.1\) makes the plate above 30.6 per cent too stiff. No single value serves both problem classes.

The proper fix is to give the element a locking-free transverse shear field, the DSG of Bletzinger, Bischoff and Ramm (2000) being the triangular counterpart of the MITC4 treatment used for quadrilaterals, after which no \(\alpha\) is needed at all. That element now exists as pyfe3d.tria3dsg, and it is the triangle to reach for by default; this one is kept for continuity with results obtained before it, and for the comparison itself. Where it is used, alpha_shear_locking is a parameter to be verified per problem class and not a constant of the element.

The transverse shear stiffnesses \(A_{44}\), \(A_{45}\) and \(A_{55}\) are read from the pyfe3d.shellprop.ShellProp object with the shear correction already applied, see pyfe3d.shellprop.ShellProp.calc_transverse_shear_stiffness(), and brought to the element coordinate system with pyfe3d.shellprop.ShellProp.calc_Ats_element(), before the stabilisation above is applied. As described in Castro et al., the shear correction is no longer a very relevant parameter when the stabilisation scheme presented above is used.

Drilling stiffness#

The FSDT kinematics contains no strain measure associated with \(r_z\), so the rows and columns of the drilling degree-of-freedom would be empty and a mesh of coplanar elements would give a singular global stiffness matrix. Two models are available, selected with the drilling_model attribute, and the formulation of each is the one documented for pyfe3d.quad4, to which the reader is referred for the derivation.

Physics-based, the default (drilling_model = 0). The in-plane displacement field is enriched with the hierarchical quadratic edge modes of Allman, so that the drilling rotations produce membrane strain energy, and the independently interpolated \(r_z\) is tied to the rotation of the membrane field by the regularised functional of Hughes and Brezzi:

Allman, D. J., 1984, “A compatible triangular element including vertex rotations for plane elasticity analysis,” Computers & Structures, 19(1-2), pp. 1-8. https://doi.org/10.1016/0045-7949(84)90197-4

Hughes, T. J. R., and Brezzi, F., 1989, “On drilling degrees of freedom,” Computer Methods in Applied Mechanics and Engineering, 72(1), pp. 105-121. https://doi.org/10.1016/0045-7825(89)90124-2

The amplitude of the mode of the edge \(k\) joining nodes \(i\) and \(j\) is \(a_k = \frac{\ell_k}{8}\left({r_z}_i - {r_z}_j\right)\), so the enrichment is the difference of the two drilling rotations already present at the ends of the edge and introduces no new degree-of-freedom. For the triangle the hierarchical bubble of that edge is

\[N_k = 4 S_i S_j\]

with \(S_i\) the area coordinates, equal to unity at the mid-point of the edge and zero at every vertex. The enrichment populates the drilling columns of the membrane operator, giving \(\pmb{\tilde B}_m\), and the drilling residual becomes \(\pmb{\tilde B}_{r_z} = \pmb{S}^{r_z} + \frac{1}{2}\pmb{\tilde S}^u_{,y} - \frac{1}{2}\pmb{\tilde S}^v_{,x}\), whose contribution to the element stiffness matrix is \(\gamma_{r_z} \int_A \pmb{\tilde B}_{r_z}^\top \pmb{\tilde B}_{r_z} dA\) with \(\gamma_{r_z} = A_{66}\), a modulus and not a user parameter. The curvature and transverse shear operators are untouched, the edge modes acting only on the in-plane translations, so the bending and the transverse shear response of the element are identical for the two drilling models.

Two quadrature choices matter. Because the derivatives of the bubbles are linear in the area coordinates, the drilling columns of \(\pmb{\tilde B}_m\) vary linearly, and the terms quadratic in them are integrated with the three-point rule of Cowper, the same one the element already used for its drilling terms. A single point at the centroid would leave those terms unsampled and would give the element five zero eigenvalues over its nine in-plane degrees-of-freedom instead of three. The Hughes-Brezzi term itself is integrated with a single point at the centroid, following Ibrahimbegovic et al. (1990), which is what makes the element insensitive to \(\gamma_{r_z}\).

Unlike the penalty below, the added term is consistent rather than artificial: stationarity with respect to \(r_z\) gives \(r_z = \theta_z\) pointwise, so it contributes no energy at the exact solution for any positive \(\gamma_{r_z}\), and the nodal moments about the shell normal recovered in the internal force vector are physical. Note that the transverse shear stabilisation documented above scales the transverse shear stiffness only and does not interact with the drilling term, which draws its scale from \(A_{66}\) of the extensional stiffness matrix.

Fictitious penalty (drilling_model = 1), the default before version 0.10.0, following the approach adopted in MSC Nastran and Autodesk Nastran through their K6ROT parameter. It provides a small artificial stiffness whose only purpose is to remove the singularity, so the forces associated with it are spurious and any moment recovered about the shell normal is meaningless. The penalty energy is defined per element as:

\[U_{drill} = \frac{1}{2} K6ROT \cdot 10^{-6} \cdot \int_A A_{66} (r_z - \theta_z)^2 dA\]

where \(10^{-6}\) is a scaling factor suggested by MSC Nastran’s approach (CQUAD4) to make the artificial drilling stiffness sufficiently small. AUTODESK NASTRAN’s quick reference guide recommends \(K6ROT = 100\) for static analysis. For modal solutions, \(K6ROT = 10^4\) is suggested. MSC NASTRAN’s quick reference guide states that \(K6ROT > 100\) should not be used, thus contradicting AUTODESK NASTRAN. The rotation \(r_z\) represents the drilling degree-of-freedom in element’s coordinates, whereas \(\theta_z\) the in-plane rotation strain, defined as \(\theta_z = \frac{1}{2}\left(v_{,x} - u_{,y}\right)\), such that the penalty is built from the operator

\[B_{drill} = S^{r_z} + 1/2 S^u_{,y} - 1/2 S^v_{,x}\]

which is the same operator as \(\pmb{\tilde B}_{r_z}\) above, evaluated on the unenriched field. Being built from that operator and not from an addition on the diagonal terms is what keeps the penalty from stiffening a rigid rotation of the element about its normal, so both models represent all rigid-body motions and all constant-strain states exactly. \(A_{66}\) is assumed constant over the element. The approach herein presented is very similar to the one presented in Eq. 2.20 of:

Adam, F. M., Mohamed, A. E., and Hassaballa, A. E., 2013, u201cDegenerated Four Nodes Shell Element with Drilling Degree of Freedom,u201d IOSR J. Eng., 3(8), pp. 10u201320.

Choosing between them. The physics-based model is the default because it is the one that is correct when the drilling moment is part of the load path, when shells are connected to beams or stiffeners that must transmit in-plane moments, or when the mesh is too coarse for the unenriched membrane response to be trusted. It is markedly more accurate in in-plane bending: on Cook’s skew membrane with a four by four mesh of split quadrilaterals it gives 20.7 against the reference 23.9, where the penalty gives 11.3. The penalty remains available for reproducing results obtained before 0.10.0.

class pyfe3d.tria3r.Tria3R#

Nodal connectivity for the triangular element similar to Nastran’s CTRIA3:

3
|\
| \    positive normal in CCW
|  \
|___\
1    2

The element coordinate system is determined identically what is explained in Nastran’s quick reference guide for the CTRIA3 element, as illustrated below.

Attributes:

eid,int

Element identification number.

pid,int

Property identification number.

area,double

Element area.

alpha_shear_locking,double

Factor used to prevent shear locking, adopted from the DFG element, affecting the transverse shear stiffness terms A44, A45, A55, already in the element coordinate system and with the shear correction applied (see ShellProp.calc_Ats_element()):

maxl = max(edge_12, edge_23, edge_31)
factor = alpha_shear_locking*maxl**2/thickness**2
A44 = 1 / (1 + factor) * A44
A45 = 1 / (1 + factor) * A45
A55 = 1 / (1 + factor) * A55

The adopted default is alpha_shear_locking = 0.7, based on a linear buckling analysis of a simply supported plate, such that the result approaches the one of the Quad4R element for an equivalent mesh (see the test case test_tria3r_linear_buckling_plate.py).

Warning

\(\alpha = 0.7\) is about seven times the value used in the references this scheme comes from, and it is the dominant error term on problems where transverse shear carries load. The reason, and the measurements, are in the section “The transverse shear stiffnesses” of the module documentation.

No single value serves every problem class, so this is a parameter to be verified per problem and not a constant of the element. Two measurements bracket it, and they pull in opposite directions:

on the plate of tests/test_tria3r_natural_freq.py with the consistent mass matrix, \(\alpha = 0.7\) gives a first natural frequency 0.6 per cent above the analytical value while the literature \(\alpha = 0.1\) gives one 30.6 per cent above it, so here the default is much the better of the two. With the stabilisation off the error is 96.1 per cent, which is how much of this element’s accuracy rests on \(\alpha\);
on the cylinder of tests/test_quad4_linear_buckling_cylinder_displ.py at ntheta = 60, where \(factor\) reaches 1852, the ordering reverses. At \(\alpha = 0.7\) the critical eigenvalue belongs to a mode at two elements per wavelength, the mesh Nyquist limit, carrying 36.6 per cent of its energy in transverse shear: a numerical mechanism rather than a physical mode. At \(\alpha = 0.1\) the physical eight-wave mode is critical instead, and there this element agrees with the quadrilateral and with the discrete shear gap triangle to 3 per cent, 2.520 against 2.453 and 2.534 times the reference load, all three overpredicting at so coarse a mesh.

So on a thin shell in buckling the default is not merely inaccurate, it can change which mode is critical, and a plausible-looking eigenvalue can belong to a mesh artefact. Where the answer matters, either sweep \(\alpha\) and confirm that the critical mode is physical and resolved by several elements per wavelength, or use Tria3DSG, whose discrete shear gap transverse shear field is locking-free and takes no such parameter.

drilling_model,int

Selects how the drilling degree-of-freedom \(r_z\) is given stiffness, see the module documentation. The default 0 is the physics-based stiffness of Allman (1984) and Hughes and Brezzi (1989), for which the drilling rotation is a kinematic variable that carries strain energy, the recovered nodal moments about the shell normal are physical, and no user parameter is involved. Any other value selects the fictitious penalty of MSC Nastran and Autodesk Nastran, which was the default up to version 0.9.0 and is controlled by K6ROT. Setting elem.drilling_model = 1 before calling update_KC0() is the way to reproduce results obtained before 0.10.0.

K6ROT,double

Dimensionless multiplier for the fictitious drilling stiffness, only read when drilling_model is not 0. It has no effect under the default physics-based model, which takes its regularisation parameter from the laminate stiffness instead, see gamma_rz. AUTODESK NASTRAN’s quick reference guide recommends K6ROT = 100. for static analysis. For modal solutions, K6ROT=1.e4 is suggested. MSC NASTRAN’s quick reference guide states that K6ROT > 100. should not be used, but this is contradicting AUTODESK NASTRAN.

gamma_rz,double

Regularisation parameter \(\gamma_{r_z}\) of the physics-based drilling stiffness, only read when drilling_model is 0. The default is a negative value, which means that \(A_{66}\) of the laminate extensional stiffness matrix is used, the value identified by Hughes and Brezzi (1989). This is a modulus and not a parameter that needs tuning: the element response has a broad plateau of insensitivity around it, and the attribute is exposed for the sensitivity study that the literature recommends rather than for normal use. Very large values over-constrain \(r_z = \theta_z\), and a zero value leaves the Allman enrichment rank-deficient by one.

r11, r12, r13, r21, r22, r23, r31, r32, r33double

Rotation matrix from local to global coordinates.

m11, m12, m21, m22double

Rotation matrix only for the constitutive relations. Used when a material direction is used instead of the element local coordinates.

c1, c2, c3: int

Position of each node in the global stiffness matrix.

n1, n2, n3: int

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.

probe,Tria3RProbe 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 constitutive stiffness matrix KCNL
`update_KG`(self, long[, long[, double[, ...)	Update sparse vectors for geometric stiffness matrix KG
`update_KG_given_stress`(self, double Nxx, ...)	Update sparse vectors for geometric stiffness matrix KG
`update_M`(self, long[, long[, double[, ...)	Update sparse vectors for mass matrix M
`update_area`(self)	Update element area
`update_fint`(self, double[, ShellProp prop, ...)	Update the internal force vector
`update_probe_finte`(self, ShellProp 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[, ...)	Update the rotation matrix of the element

K6ROT#: K6ROT: ‘double’

alpha_shear_locking#: alpha_shear_locking: ‘double’

area#: area: ‘double’

c1#: c1: ‘int’

c2#: c2: ‘int’

c3#: c3: ‘int’

drilling_model#: drilling_model: ‘int’

eid#: eid: ‘int’

gamma_rz#: gamma_rz: ‘double’

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’

m11#: m11: ‘double’

m12#: m12: ‘double’

m21#: m21: ‘double’

m22#: m22: ‘double’

n1#: n1: ‘int’

n2#: n2: ‘int’

n3#: n3: ‘int’

pid#: pid: ‘int’

probe#: probe: pyfe3d.tria3r.Tria3RProbe

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, ShellProp prop, int update_KC0v_only=0) → void#

Update sparse vectors for linear constitutive stiffness matrix KC0

Parameters:

KC0rnp.array: Array to store row positions of sparse values
KC0cnp.array: Array to store column positions of sparse values
KC0vnp.array: Array to store sparse values
propShellProp object: Shell property object from where the stiffness and mass attributes are read from.
update_KC0v_onlyint: 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.

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

Update sparse vectors for the nonlinear constitutive stiffness matrix KCNL

Assuming that KCNL = KC0L + KCL0 + KCLL + KGNL, built from the von Karman membrane strains

\[\epsilon_{xx} = u_{,x} + \frac{1}{2} w_{,x}^2, \quad \epsilon_{yy} = v_{,y} + \frac{1}{2} w_{,y}^2, \quad \gamma_{xy} = u_{,y} + v_{,x} + w_{,x} w_{,y}\]

whose nonlinear part is \(\{\epsilon_{NL}\} = \frac{1}{2} [B_{mL}] \{u_e\}\), with \([B_{mL}]\) its variation. With \([B_m]\) and \([B_b]\) the linear membrane and bending strain-displacement matrices, \([G]\) the gradient of \(w\), and \([A]\), \([B]\) the laminate matrices:

KC0L = \([B_m]^T [A] [B_{mL}] + [B_b]^T [B] [B_{mL}]\)
KCL0 = KC0L`^T`
KCLL = \([B_{mL}]^T [A] [B_{mL}]\)
KGNL = \([G]^T [N_{NL}] [G]\), with \(\{N_{NL}\} = [A] \{\epsilon_{NL}\}\)

The first three groups are the constitutive terms coupling the linear and the nonlinear parts of the membrane strain. KGNL is geometric, carrying the stress of the nonlinear membrane 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,

\[K_T = K_{C0} + K_{CNL}(u) + K_G(u)\]

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

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

Parameters:

KCNLrnp.array: Array to store row positions of sparse values
KCNLcnp.array: Array to store column positions of sparse values
KCNLvnp.array: Array to store sparse values
propShellProp object: Shell property object from where the stiffness and mass attributes are read from.
update_KCNLv_onlyint: 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.

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

Update sparse vectors for geometric stiffness matrix KG

Two-point Gauss-Legendre quadrature is used, which showed more accuracy for linear buckling load predictions.

Before this function is called, the probe Tria3RProbe attribute of the Tria3R object must be updated using update_probe_ue() with the correct pre-buckling displacements; and update_probe_xe() with the node coordinates.

Parameters:

KGrnp.array: Array to store row positions of sparse values
KGcnp.array: Array to store column positions of sparse values
KGvnp.array: Array to store sparse values
propShellProp object: Shell property object from where the stiffness and mass attributes are read from.
update_KGv_onlyint: The default \(0\) means that only \(KGv\) is updated. Any other value will lead to \(KGr\) and \(KGc\) also being updated.

update_KG_given_stress(self, double Nxx, double Nyy, double Nxy, long[: :1] KGr, long[: :1] KGc, double[: :1] KGv, int update_KGv_only=0) → void#

Update sparse vectors for geometric stiffness matrix KG

Note

A constant stress state is assumed within the element, according to the given values of \(N_{xx}, N_{yy}, N_{xy}\).

Two-point Gauss-Legendre quadrature is used, which showed more accuracy for linear buckling load predictions.

Before this function is called, the probe Tria3RProbe attribute of the Tria3R object must be updated using update_probe_xe() with the node coordinates.

Parameters:

KGrnp.array: Array to store row positions of sparse values
KGcnp.array: Array to store column positions of sparse values
KGvnp.array: Array to store sparse values
update_KGv_onlyint: The default \(0\) means that only \(KGv\) is updated. Any other value will lead to \(KGr\) and \(KGc\) also being updated.

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

Update sparse vectors for mass matrix M

Different integration schemes are available by means of the mtype parameter.

Parameters:

Mrnp.array: Array to store row positions of sparse values
Mcnp.array: Array to store column positions of sparse values
Mvnp.array: Array to store sparse values
mtypeint, optional: 0 for consistent mass matrix using method from Brockman 1987 1 for reduced integration mass matrix using method from Brockman 1987 2 for lumped mass matrix using method from Brockman 1987

update_area(self) → void#: Update element area

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

Update the internal force vector

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 Tria3RProbe with the internal forces in local coordinates.
propShellProp object: Shell property object from where the stiffness and mass attributes are read from.
nonlinearint: The default 0 gives the linear internal forces, KC0*u. Any other value adds the geometrically nonlinear terms of the von Karman strains, for which the exact Jacobian of the internal forces is KC0 + KCNL + KG, see update_KCNL().

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

Update the internal force vector of the probe

The attribute finte is updated with the Tria3RProbe 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 Tria3RProbe is updated, accessible using .probe.finte.

Parameters:

propShellProp object: Shell property object from where the stiffness and mass attributes are read from.
nonlinearint: The default 0 gives the linear internal forces, KC0*u. Any other value adds the geometrically nonlinear terms of the von Karman strains, 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 Tria3RProbe is updated, accessible using .probe.ue.

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\).

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

Update the 3D coordinates of the probe of the element

Note

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

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\).

update_rotation_matrix(self, double[: :1] x, double xmati=0., double xmatj=0., double xmatk=0.) → void#

Update the rotation matrix of the element

Attributes r11,r12,r13,r21,r22,r23,r31,r32,r33 are updated, corresponding to the rotation matrix from local to global coordinates.

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\).

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\).
xmati, xmatj, xmatk: array-like: Vector in global coordinates representing the material direction. This vector is projected onto the plate element, thus becoming the material direction. The \(ABD\) matrix defining the constitutive behavior of the element is rotated from the material direction to the element \(x\) axis while calculating the stiffness matrices.

class pyfe3d.tria3r.Tria3RData#

Used to allocate memory for the sparse matrices.

Attributes:

KC0_SPARSE_SIZE,int: KC0_SPARSE_SIZE = 324
KCNL_SPARSE_SIZE,int: KCNL_SPARSE_SIZE = 324
KG_SPARSE_SIZE,int: KG_SPARSE_SIZE = 81
M_SPARSE_SIZE,int: M_SPARSE_SIZE = 270

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.tria3r.Tria3RProbe#

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.

Attributes:

xe,array-like: Array of size NUM_NODES*DOF//2=9 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\), \({x_e}_3, {y_e}_3, {z_e}_3\).
ue,array-like: Array of size NUM_NODES*DOF=18 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\), \({u_e}_3, {v_e}_3, {w_e}_3, {{r_x}_e}_3, {{r_y}_e}_3, {{r_z}_e}_3\).
finte,array-like: Array of size NUM_NODES*DOF=18 containing the element internal forces corresponding to the degrees-of-freedom described by ue.
BLexx, BLeyy, BLgxy, BLkxx, BLkyy, BLkxyarray-like: Arrays of size NUM_NODES*DOF=18 with the rows of the linear strain-displacement matrix for the membrane strains and curvatures, at the last evaluated integration point.
Gwx, Gwyarray-like: Arrays of size NUM_NODES*DOF=18 with the rows giving \(w_{,x}\) and \(w_{,y}\), at the last evaluated integration point.
KCNLvearray-like: KCNLve: ‘double[::1]’

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

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

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

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

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

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

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

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

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

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

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

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

Tria3R - Triangular element with reduced integration (pyfe3d.tria3r)#

Drilling stiffness#

Tria3R - Triangular element with reduced integration (`pyfe3d.tria3r`)#