Donnell-type cylinder element (`bfsccylinder.bfsccylinder`)#

Donnell-type cylinder element (bfsccylinder.bfsccylinder)

The module provides the single-curvature BFSC element with Donnell-type kinematics, the class BFSCCylinder, and the functions update_KC0(), update_KCNL(), update_KG(), update_KG_constant_stress() and update_M(), that write the contributions of one element to the linear constitutive stiffness matrix KC0, the nonlinear stiffness matrix KCNL, the sum of the nonlinear constitutive and geometric stiffness matrices \(\boldsymbol{K}_{CNL} + \boldsymbol{K}_{GNL}\), the geometric stiffness matrix KG and the mass matrix M in sparse row, column and value arrays, starting at the positions given by its init_k_* attributes, and update_fint(), that adds its internal forces to a global internal force vector. For a pressure that follows the deformation of the shell, update_fext_follower() adds the force vector of one element to a global external force vector and update_KCfollower() writes its load stiffness matrix KCfollower in sparse arrays, see the Notes. The row, column and value arrays of the elements can be gathered in a sparse matrix of coordinate format, e.g. with scipy.sparse.coo_matrix, which sums the values at repeated positions.

Each of the four nodes has DOF = 10 degrees of freedom, in the order \(u\), \(u_{,x}\), \(u_{,y}\), \(v\), \(v_{,x}\), \(v_{,y}\), \(w\), \(w_{,x}\), \(w_{,y}\), \(w_{,xy}\), where \(x\) is the axial and \(y\) the circumferential coordinate, and \(u\), \(v\) and \(w\) are the axial, circumferential and normal displacements. Element vectors of size 4*DOF = 40 have the degrees of freedom of nodes 1, 2, 3 and 4 in sequence.

The module also defines:

KC0_SPARSE_SIZE = 1600, KCNL_SPARSE_SIZE = 1024, KG_SPARSE_SIZE = 256, M_SPARSE_SIZE = 544 and KCfollower_SPARSE_SIZE = 1312: number of entries that one element writes in the row, column and value arrays of KC0, KCNL, KG, M and KCfollower. For n elements the arrays need n times these sizes;
DOUBLE = np.float64: NumPy type of the arrays of values.

The module bfsccylinder.bfsccylinder_sanders has the same contents for the element with Sanders-type kinematics.

Notes#

The element is a rectangle of the developed mid-surface of a cylinder of radius \(R\) = R, with \(x\) the axial coordinate and \(y = R \theta\) the arc length along the circumference, of sides \(l_{ex}\) = lex along \(x\) and \(l_{ey}\) = ley along \(y\). The natural coordinates \(\xi\) and \(\eta\), from \(-1\) to \(1\), give \(x = x_1 + l_{ex} (\xi + 1)/2\) and \(y = y_1 + l_{ey} (\eta + 1)/2\), with \((x_1, y_1)\) the position of node 1, such that \(\partial/\partial x = (2/l_{ex}) \partial/\partial \xi\), \(\partial/\partial y = (2/l_{ey}) \partial/\partial \eta\) and \(dA = (l_{ex} l_{ey}/4) d\xi d\eta\).

Interpolation

For the node \(i\) at the corner \((\xi_i, \eta_i)\), \((-1, -1)\), \((+1, -1)\), \((+1, +1)\) and \((-1, +1)\) for the nodes 1, 2, 3 and 4, the interpolation functions are products of cubic Hermite polynomials in \(\xi\) and \(\eta\) [2]:

\[H_i = \frac{1}{16} (\xi + \xi_i)^2 (\xi \xi_i - 2) (\eta + \eta_i)^2 (\eta \eta_i - 2)\]

\[H_{xi} = -\frac{l_{ex}}{32} \xi_i (\xi + \xi_i)^2 (\xi \xi_i - 1) (\eta + \eta_i)^2 (\eta \eta_i - 2)\]

\[H_{yi} = -\frac{l_{ey}}{32} (\xi + \xi_i)^2 (\xi \xi_i - 2) \eta_i (\eta + \eta_i)^2 (\eta \eta_i - 1)\]

\[H_{xyi} = \frac{l_{ex} l_{ey}}{64} \xi_i (\xi + \xi_i)^2 (\xi \xi_i - 1) \eta_i (\eta + \eta_i)^2 (\eta \eta_i - 1)\]

At the nodes, \(H_i\) is one at its node and zero at the others, and its derivatives \(\partial/\partial x\), \(\partial/\partial y\) and \(\partial^2/\partial x \partial y\) are zero at all nodes; \(H_{xi}\), \(H_{yi}\) and \(H_{xyi}\) are zero at all nodes, and so are these derivatives, except \(\partial H_{xi}/\partial x = \partial H_{yi}/\partial y = \partial^2 H_{xyi}/\partial x \partial y = 1\) at their own node. The factors \(l_{ex}/2\) and \(l_{ey}/2\) convert the derivatives in \(\xi\) and \(\eta\) into derivatives in \(x\) and \(y\). The normal displacement is interpolated as in the element of Bogner, Fox and Schmit [1],

\[w = \sum_{i=1}^{4} \left( H_i w_i + H_{xi} {w_{,x}}_i + H_{yi} {w_{,y}}_i + H_{xyi} {w_{,xy}}_i \right)\]

and the in-plane displacements with the same functions but without the cross derivative, which is not a degree of freedom of \(u\) and \(v\):

\[u = \sum_{i=1}^{4} \left( H_i u_i + H_{xi} {u_{,x}}_i + H_{yi} {u_{,y}}_i \right), \qquad v = \sum_{i=1}^{4} \left( H_i v_i + H_{xi} {v_{,x}}_i + H_{yi} {v_{,y}}_i \right)\]

In row form \(u = \boldsymbol{S}^u\boldsymbol{u}_e\), \(v = \boldsymbol{S}^v\boldsymbol{u}_e\) and \(w = \boldsymbol{S}^w\boldsymbol{u}_e\), with \(\boldsymbol{u}_e\) the element displacement vector u and the rows Su, Sv and Sw of BFSCCylinder.update_Su(), BFSCCylinder.update_Sv() and BFSCCylinder.update_Sw(); the rows of the derivatives, e.g. Sw_x, are their derivatives. Along an edge, \(u\), \(v\), \(w\) and their derivatives normal to the edge depend only on the degrees of freedom of the two nodes of that edge, such that \(u\), \(v\), \(w\) and their first derivatives are continuous between elements that share an edge.

This interpolation of \(u\) and \(v\) with the Hermite functions is that of the BFSC plate element of Castro and Jansen [4], of which this element is the single-curvature version, first presented by Castro and Jansen [5] with Donnell-type and Sanders-type kinematics and later used by Wang et al. [6] and by dos Santos and Castro [7], who give its interpolation functions and the order of its nodal degrees of freedom.

Kinematics

The kinematics are those of Donnell [3] for a thin cylindrical shell under the Kirchhoff-Love hypothesis, such that transverse shear deformation is neglected, with \(w\) positive outwards, away from the axis of the cylinder. The membrane strains \(\boldsymbol{\varepsilon} = \{\varepsilon_{xx}, \varepsilon_{yy}, \gamma_{xy}\}^\top\) are the sum of a linear part \(\boldsymbol{\varepsilon}_L\) and of a nonlinear part \(\boldsymbol{\varepsilon}_{NL}\) in the slopes of \(w\), for small strains and moderate rotations, and the curvatures \(\boldsymbol{\kappa} = \{\kappa_{xx}, \kappa_{yy}, \kappa_{xy}\}^\top\) are linear:

\[\begin{split}\boldsymbol{\varepsilon}_L = \begin{Bmatrix} u_{,x} \\ v_{,y} + w/R \\ u_{,y} + v_{,x} \end{Bmatrix} = \boldsymbol{B}_m\boldsymbol{u}_e, \qquad \boldsymbol{\kappa} = \begin{Bmatrix} -w_{,xx} \\ -w_{,yy} \\ -2 w_{,xy} \end{Bmatrix} = \boldsymbol{B}_b\boldsymbol{u}_e\end{split}\]

\[\begin{split}\boldsymbol{\varepsilon}_{NL} = \begin{Bmatrix} w_{,x}^2/2 \\ w_{,y}^2/2 \\ w_{,x} w_{,y} \end{Bmatrix} = \frac{1}{2} \boldsymbol{B}_{NL}\boldsymbol{u}_e, \qquad \boldsymbol{B}_{NL} = \begin{bmatrix} w_{,x} \boldsymbol{S}^w_{,x} \\ w_{,y} \boldsymbol{S}^w_{,y} \\ w_{,x} \boldsymbol{S}^w_{,y} + w_{,y} \boldsymbol{S}^w_{,x} \end{bmatrix}\end{split}\]

where \(\boldsymbol{B}_m\) and \(\boldsymbol{B}_b\) are the rows Bm and Bb of BFSCCylinder.update_Bm() and BFSCCylinder.update_Bb(), and \(\boldsymbol{B}_{NL}\), which depends on the slopes w_x and w_y of the current displacements, has the rows BmL1, BmL2 and BmL3 of update_KCNL() and update_fint(), such that the variation of the nonlinear strains is \(\delta\boldsymbol{\varepsilon}_{NL} = \boldsymbol{B}_{NL}\delta \boldsymbol{u}_e\). The slopes are \(\{w_{,x}, w_{,y}\}^\top = \boldsymbol{G}\boldsymbol{u}_e\), with \(\boldsymbol{G}\) of the rows \(\boldsymbol{S}^w_{,x}\) and \(\boldsymbol{S}^w_{,y}\), G1 and G2 in the code.

The in-plane displacements do not enter the curvatures and the rotations, which are those of a flat plate, and the curvature of the shell enters only through the term \(w/R\) of \(\varepsilon_{yy}\). These approximations are accurate when the displacements vary over lengths that are small compared with \(R\), e.g. for many circumferential waves, and lose accuracy for long waves. For instance, the curvatures do not vanish for the small rigid-body translation \(w = c \cos(y/R)\), \(v = -c \sin(y/R)\), for which \(\kappa_{yy} = w/R^2\); the kinematics of bfsccylinder.bfsccylinder_sanders give zero strains and curvatures for it.

Constitutive relation

The membrane stress resultants \(\boldsymbol{N} = \{N_{xx}, N_{yy}, N_{xy}\}^\top\) and the moment resultants \(\boldsymbol{M} = \{M_{xx}, M_{yy}, M_{xy}\}^\top\) are those of the classical laminated plate theory,

\[\begin{split}\begin{Bmatrix} \boldsymbol{N} \\ \boldsymbol{M} \end{Bmatrix} = \begin{bmatrix} \boldsymbol{A} & \boldsymbol{B} \\ \boldsymbol{B} & \boldsymbol{D} \end{bmatrix} \begin{Bmatrix} \boldsymbol{\varepsilon}_L + \boldsymbol{\varepsilon}_{NL} \\ \boldsymbol{\kappa} \end{Bmatrix}, \qquad \boldsymbol{A} = \begin{bmatrix} A_{11} & A_{12} & A_{16} \\ A_{12} & A_{22} & A_{26} \\ A_{16} & A_{26} & A_{66} \end{bmatrix}\end{split}\]

and likewise \(\boldsymbol{B}\) and \(\boldsymbol{D}\), with the values of A11, ..., D66 of the element at each integration point. They are the stiffnesses of a flat laminate referred to the mid-surface, without terms in \(h/R\).

Element matrices and internal forces

The strain energy of the element is

\[U = \frac{1}{2} \int_A \left( \boldsymbol{\varepsilon}^\top \boldsymbol{A} \boldsymbol{\varepsilon} + 2 \boldsymbol{\varepsilon}^\top \boldsymbol{B} \boldsymbol{\kappa} + \boldsymbol{\kappa}^\top \boldsymbol{D} \boldsymbol{\kappa} \right) dA\]

with \(\boldsymbol{\varepsilon} = \boldsymbol{\varepsilon}_L + \boldsymbol{\varepsilon}_{NL}\). The internal force vector of update_fint() is its gradient \(\partial U/\partial \boldsymbol{u}_e\), and the tangent stiffness matrix, its Hessian, is split into

\[\boldsymbol{K}_T = \frac{\partial^2 U}{\partial \boldsymbol{u}_e^2} = \boldsymbol{K}_{C0} + \boldsymbol{K}_{CNL}(\boldsymbol{u}_e) + \boldsymbol{K}_{GNL}(\boldsymbol{u}_e) + \boldsymbol{K}_G(\boldsymbol{u}_e)\]

computed by update_KC0(), update_KCNL() and update_KG(), whose Notes give each term. The mass matrix of update_M() comes from the kinetic energy, see its Notes.

The integrals over the element are computed with the Gauss-Legendre quadrature of bfsccylinder.quadrature.get_points_weights() in each direction,

\[\int_A f \, dA \approx \frac{l_{ex} l_{ey}}{4} \sum_{i=0}^{nint-1} \sum_{j=0}^{nint-1} weights_i \, weights_j \, f(points_i, points_j)\]

and the mass matrix in closed form. Since KCNL and KG are the derivatives of the same quadrature sums as the internal forces, \(\boldsymbol{K}_T\) is the exact Jacobian of the internal forces for any nint. With \(\boldsymbol{A}\), \(\boldsymbol{B}\) and \(\boldsymbol{D}\) constant over the element, the integrands are polynomials in \(\xi\) and \(\eta\): of degree up to 6 in each for KC0 and for update_KG_constant_stress(), integrated exactly with nint >= 4; up to 9 for KG, exact with nint >= 5; and up to 12 for KCNL and the internal forces, exact with nint >= 7.

Follower pressure

A pressure \(p\), per unit area of the undeformed mid-surface and a function of the undeformed coordinates, that stays normal to the deformed mid-surface and acts on its deformed area, a follower pressure, does the virtual work \(\delta W_p = \int_A p \boldsymbol{a}^\top \delta \boldsymbol{u} dA\), with \(\boldsymbol{a} = \boldsymbol{X}_{,x} \times \boldsymbol{X}_{,y}\) the area vector of the deformed mid-surface \(\boldsymbol{X} = \boldsymbol{X}_0 + \boldsymbol{u}\) per unit undeformed area, see the Notes of update_fext_follower(). With the Donnell-type kinematics the rotation of the normal about \(x\) is \(-w_{,y}\), the term \(v/R\) being neglected as in the strains. The area vector is split into its parts of order 0, 1 and 2 in the displacement gradients, \(\boldsymbol{a} = \boldsymbol{a}_0 + \boldsymbol{a}_1 + \boldsymbol{a}_2\), with \(\boldsymbol{a}_0 = \{0, 0, 1\}^\top\) the dead load and \(\boldsymbol{a}_1\) the rotation of the normal and the change of area. The truncation \(\boldsymbol{a}_0 + \boldsymbol{a}_1\), quadratic=0, is consistent with the moderate rotations of the strains, and \(\boldsymbol{a}_2\), added with quadratic=1, is of the order of the terms that the strains neglect. A positive \(p\) pushes along \(+z\), outwards, the direction of a positive \(w\), such that an external pressure is negative.

The force vector \(\boldsymbol{F}_p(\boldsymbol{u})\) of update_fext_follower() enters the residual \(\boldsymbol{R} = \boldsymbol{F}_{int}(\boldsymbol{u}) - \boldsymbol{F}_p(\boldsymbol{u})\), whose Jacobian is \(\boldsymbol{K}_{C0} + \boldsymbol{K}_{CNL} + \boldsymbol{K}_{GNL} + \boldsymbol{K}_G + \boldsymbol{K}_F\), with the load stiffness matrix \(\boldsymbol{K}_F = -\partial \boldsymbol{F}_p/\partial \boldsymbol{u}\) of update_KCfollower(). For quadratic=0, \(\boldsymbol{F}_p\) is affine in the displacements and \(\boldsymbol{K}_F\) is constant. The load stiffness matrix is unsymmetric in general. For quadratic=0, the skew part of the bilinear form \(P(\delta u, \Delta u) = \int_A p \delta \boldsymbol{u}^\top \boldsymbol{a}_1(\Delta \boldsymbol{u}) dA\) over the loaded region is, for both kinematics,

\[P(\delta u, \Delta u) - P(\Delta u, \delta u) = -\oint p \left( (\delta u \Delta w - \delta w \Delta u) n_x + (\delta v \Delta w - \delta w \Delta v) n_y \right) ds + \int_A \left( p_{,x} (\delta u \Delta w - \delta w \Delta u) + p_{,y} (\delta v \Delta w - \delta w \Delta v) \right) dA\]

with \(\{n_x, n_y\}^\top\) the outward normal of the contour of the loaded region. The curvature terms of the cylinder are symmetric and cancel in the skew part, which therefore holds for both kinematics. The load stiffness matrix is therefore symmetric, and the follower pressure conservative, if and only if:

\(p\) is uniform over the loaded region, such that the area integral vanishes; and
at each point of the contour of the loaded region, either \(w = 0\) or the in-plane displacement normal to the contour vanishes, \(u = 0\) on an edge \(x\) = constant and \(v = 0\) on an edge \(y\) = constant, such that the contour integral vanishes.

These conditions are met by simply supported and clamped edges, by symmetry planes, and by closed surfaces such as complete rings, cylinders and spheres, for which no boundary exists; for a cylinder closed around the circumference, the contour of the loaded region is that of its ends. In these cases the load potential is the pressure times the change of the enclosed volume. Conversely, the matrix is unsymmetric for free or partially restrained edges, for a pressure applied on a patch whose edges lie inside the shell, and for a pressure that varies over the surface. This is the result of Hibbitt [8] and of Schweizerhof and Ramm [9], who showed that the pressure load stiffness is symmetric for closed surfaces or for loaded surfaces with a fixed boundary. For the matrix of one element, the contour of the loaded region is the boundary of the element, where the second condition does not hold, and the element matrix is unsymmetric; in the assembled matrix the contour integrals along the edges shared by two elements cancel when \(p\) is continuous, since \(u\), \(v\) and \(w\) are continuous, such that the conditions apply to the contour of the loaded region of the shell. A pressure that depends on the current position, e.g. a hydrostatic head, is not modelled. The conditions are those of the first-order area vector; for quadratic=1 the matrix depends on \(\boldsymbol{u}\). With \(p\) uniform over the element, the integrands of the force vector and of the load stiffness matrix are polynomials in \(\xi\) and \(\eta\), of degree up to 6 in each for quadratic=0, integrated exactly with nint >= 4, and up to 9 for quadratic=1, exact with nint >= 5.

References#

[1]

Bogner, F. K., Fox, R. L., and Schmit, L. A., The generation of inter-element-compatible stiffness and mass matrices by the use of interpolation formulas, Proceedings of the Conference on Matrix Methods in Structural Mechanics, AFFDL-TR-66-80, Wright-Patterson Air Force Base, Dayton, Ohio, 26-28 October 1965, pp. 395-444.

[2]

Ochoa, O. O., and Reddy, J. N., Finite Element Analysis of Composite Laminates, Springer, Dordrecht, 1992.

[3]

Donnell, L. H., Stability of thin-walled tubes under torsion, NACA Report 479, 1933.

[4]

Castro, S. G. P., and Jansen, E. L., Displacement-based formulation of Koiter’s method: application to multi-modal post-buckling finite element analysis of plates, Thin-Walled Structures, 159, 107217, 2021, doi:10.1016/j.tws.2020.107217.

[5]

Castro, S. G. P., and Jansen, E. L., Displacement-based multi-modal formulation of Koiter’s method applied to cylindrical shells, AIAA SciTech 2022 Forum, AIAA 2022-0256, 2022, doi:10.2514/6.2022-0256.

[6]

Wang, Z., Almeida Jr, J. H. S., Ashok, A., Wang, Z., and Castro, S. G. P., Lightweight design of variable-angle filament-wound cylinders combining Kriging-based metamodels with particle swarm optimization, Structural and Multidisciplinary Optimization, 65, 140, 2022, doi:10.1007/s00158-022-03227-8.

[7]

dos Santos, R. R., and Castro, S. G. P., Lightweight design of variable-stiffness cylinders with reduced imperfection sensitivity enabled by continuous tow shearing and machine learning, Materials, 15(12), 4117, 2022, doi:10.3390/ma15124117.

[8]

Hibbitt, H. D., Some follower forces and load stiffness, International Journal for Numerical Methods in Engineering, 14(6), 937-941, 1979, doi:10.1002/nme.1620140613.

[9]

Schweizerhof, K., and Ramm, E., Displacement dependent pressure loads in nonlinear finite element analyses, Computers & Structures, 18(6), 1099-1114, 1984, doi:10.1016/0045-7949(84)90154-8.

class bfsccylinder.bfsccylinder.BFSCCylinder#

Single-curvature BFSC element with Donnell-type kinematics.

Nodal connectivity for cylindrical shell element similar to Nastran’s CQUAD4:

 ^ y (circumferential direction)
 |

4 ________ 3
 |       |
 |       |   --> x (axial direction)
 |       |
 |_______|
1         2

The constructor sets n1, n2, n3, n4 and c1, c2, c3, c4 to -1, lex, ley, rho, h and R to zero, and allocates every array filled with zeros; the init_k_* positions are zero. Before the functions of this module are called with the element, the caller sets c1, c2, c3, c4, the init_k_* positions of the matrices being computed, lex, ley and R, the stiffnesses A11, ..., D66, e.g. with bfsccylinder.utils.assign_constant_ABD(), and, for update_M(), rho and h.

Parameters:

nintint: Number of integration points in each direction: the size of the points and weights arrays later given to the functions of this module, which read the stiffness arrays A11, ..., D66 at the indices [i, j] from 0 to nint - 1.

Attributes:

n1, n2, n3, n4int: Node identification numbers. Not used by the functions of this module.
c1, c2, c3, c4int: Position of the first degree of freedom of nodes 1, 2, 3 and 4 in the global vectors and matrices. The DOF degrees of freedom of node 1 are at the positions c1 to c1 + 9, and likewise for the other nodes, e.g. c1 = DOF*i for the node of index i.
init_k_KC0, init_k_KCNL, init_k_KG, init_k_M, init_k_KCfollowerint: Position in the sparse arrays from which update_KC0(), update_KCNL(), update_KG() or update_KG_constant_stress(), update_M() and update_KCfollower() write the entries of the element, e.g. the index of the element times KC0_SPARSE_SIZE, KCNL_SPARSE_SIZE, KG_SPARSE_SIZE, M_SPARSE_SIZE and KCfollower_SPARSE_SIZE.
lex, leydouble: Length of the element along the axial direction \(x\) and along the circumferential direction \(y\), dimension of length.
rhodouble: Density, mass per unit volume. Used only by update_M().
hdouble: Thickness, dimension of length. Used only by update_M().
Rdouble: Radius of the cylinder, dimension of length. Terms of KC0, KCNL and KG, of the internal forces, of the follower pressure and of Bm are divided by R, so it must be set to a nonzero value before these are computed. Not used by update_KG_constant_stress() and update_M().
A11, A12, A16, A22, A26, A662-D array: Membrane stiffnesses, force per unit length. Shape (nint, nint), the entry [i, j] being used at the integration point xi = points[i], eta = points[j].
B11, B12, B16, B22, B26, B662-D array: Membrane-bending coupling stiffnesses, dimension of force. Same shape and use as A11.
D11, D12, D16, D22, D26, D662-D array: Bending stiffnesses, force times length. Same shape and use as A11.
Bm2-D array: Shape (3, 40), rows giving the linear membrane strains at the last point given to update_Bm(). Dimension of one over the dimension of each degree of freedom.
Bb2-D array: Shape (3, 40), rows giving the curvatures at the last point given to update_Bb(). Dimension of one over length times one over the dimension of each degree of freedom.
Su, Sv, Sw1-D array: Size 40, rows giving the displacements \(u\), \(v\) and \(w\) at the last point given to update_Su(), update_Sv() and update_Sw().
Su_x, Sv_x, Sw_x1-D array: Size 40, rows giving \(u_{,x}\), \(v_{,x}\) and \(w_{,x}\) at the last point given to update_Su_x(), update_Sv_x() and update_Sw_x().
Su_y, Sv_y, Sw_y1-D array: Size 40, rows giving \(u_{,y}\), \(v_{,y}\) and \(w_{,y}\) at the last point given to update_Su_y(), update_Sv_y() and update_Sw_y().
u1-D array: Size 40, element displacement vector, copied from the global displacement vector by update_KCNL(), update_KG() and update_fint(). Dimension of length for \(u\), \(v\), \(w\), dimensionless for the first derivatives and of one over length for \(w_{,xy}\).

Notes

The element and its theory are described in the Notes of bfsccylinder.bfsccylinder. The arrays A11, ..., D66 hold the laminate stiffness matrices \(\boldsymbol{A}\), \(\boldsymbol{B}\) and \(\boldsymbol{D}\) of the constitutive relation at each integration point, such that stiffnesses that vary over the element, e.g. of a laminate whose fiber angles vary, are integrated point by point. Bm and Bb are the operators \(\boldsymbol{B}_m\) and \(\boldsymbol{B}_b\) of the linear membrane strains and of the curvatures, and Su, …, Sw_y the rows of the interpolation of the displacements and of their derivatives; the update_* functions of the module evaluate the same operators at each integration point in local variables.

update_Bb(xi, eta)#

Update Bb at a natural coordinate point

The rows of Bb are computed at xi, eta from lex and ley. Their products with an element displacement vector, e.g. u, give at that point the curvatures \(\kappa_{xx}\), \(\kappa_{yy}\) and \(\kappa_{xy}\) for the rows 0, 1 and 2, dimension of one over length. Only the entries at the degrees of freedom \(w\), \(w_{,x}\), \(w_{,y}\), \(w_{,xy}\) of the nodes are written, the other entries keep the zeros of the constructor. The functions of this module do not read Bb.

Parameters:

xi, etadouble: Natural coordinates of the point, dimensionless, from -1 to 1 over the element, xi along \(x\) and eta along \(y\).

Notes

With the rows \(\boldsymbol{S}^w\) of the interpolation of \(w\) of the Notes of bfsccylinder.bfsccylinder, the rows of the Donnell-type curvatures are

\[\begin{split}\boldsymbol{B}_b = -\begin{bmatrix} \boldsymbol{S}^w_{,xx} \\ \boldsymbol{S}^w_{,yy} \\ 2 \boldsymbol{S}^w_{,xy} \end{bmatrix}\end{split}\]

such that \(\boldsymbol{\kappa} = \{-w_{,xx}, -w_{,yy}, -2 w_{,xy}\}^\top = \boldsymbol{B}_b\boldsymbol{u}_e\). They do not depend on R: the curvatures are those of a flat plate.

update_Bm(xi, eta)#

Update Bm at a natural coordinate point

The rows of Bm are computed at xi, eta from lex, ley and R. Their products with an element displacement vector, e.g. u, give at that point the linear membrane strains \(\varepsilon_{xx}\), \(\varepsilon_{yy}\) and \(\gamma_{xy}\) for the rows 0, 1 and 2, dimensionless. Only the entries that are not identically zero are written, the other entries keep the zeros of the constructor. The functions of this module do not read Bm.

Parameters:

xi, etadouble: Natural coordinates of the point, dimensionless, from -1 to 1 over the element, xi along \(x\) and eta along \(y\).

Notes

With the rows of the interpolation of the Notes of bfsccylinder.bfsccylinder, the rows of the Donnell-type linear membrane strains are

\[\begin{split}\boldsymbol{B}_m = \begin{bmatrix} \boldsymbol{S}^u_{,x} \\ \boldsymbol{S}^v_{,y} + \boldsymbol{S}^w/R \\ \boldsymbol{S}^u_{,y} + \boldsymbol{S}^v_{,x} \end{bmatrix}\end{split}\]

such that \(\boldsymbol{\varepsilon}_L = \{u_{,x}, v_{,y} + w/R, u_{,y} + v_{,x}\}^\top = \boldsymbol{B}_m\boldsymbol{u}_e\).

update_Su(xi, eta)#

Update Su at a natural coordinate point

Su is computed at xi, eta from lex and ley. Its product with an element displacement vector, e.g. u, gives the displacement \(u\) at that point, dimension of length. Only the entries at the degrees of freedom \(u\), \(u_{,x}\), \(u_{,y}\) of the nodes, the positions 0, 1, 2 plus 0, 10, 20 and 30, are written, the other entries keep the zeros of the constructor.

Parameters:

xi, etadouble: Natural coordinates of the point, dimensionless, from -1 to 1 over the element, xi along \(x\) and eta along \(y\).

Notes

The entries of node i, from 1 to 4, at the positions 10*(i - 1) + 0, 10*(i - 1) + 1 and 10*(i - 1) + 2, are the Hermite functions \(H_i\), \(H_{xi}\) and \(H_{yi}\) of the Notes of bfsccylinder.bfsccylinder evaluated at xi, eta, such that \(u = \boldsymbol{S}^u\boldsymbol{u}_e\).

update_Su_x(xi, eta)#

Update Su_x at a natural coordinate point

Su_x is computed at xi, eta from lex and ley. Its product with an element displacement vector, e.g. u, gives the derivative \(u_{,x}\) at that point, dimensionless. Only the entries at the degrees of freedom \(u\), \(u_{,x}\), \(u_{,y}\) of the nodes, the positions 0, 1, 2 plus 0, 10, 20 and 30, are written, the other entries keep the zeros of the constructor.

Parameters:

xi, etadouble: Natural coordinates of the point, dimensionless, from -1 to 1 over the element, xi along \(x\) and eta along \(y\).

Notes

The entries are the derivatives \(\partial/\partial x = (2/l_{ex}) \partial/\partial \xi\) of those of BFSCCylinder.update_Su(), with the Hermite functions of the Notes of bfsccylinder.bfsccylinder.

update_Su_y(xi, eta)#

Update Su_y at a natural coordinate point

Su_y is computed at xi, eta from lex and ley. Its product with an element displacement vector, e.g. u, gives the derivative \(u_{,y}\) at that point, dimensionless. Only the entries at the degrees of freedom \(u\), \(u_{,x}\), \(u_{,y}\) of the nodes, the positions 0, 1, 2 plus 0, 10, 20 and 30, are written, the other entries keep the zeros of the constructor.

Parameters:

xi, etadouble: Natural coordinates of the point, dimensionless, from -1 to 1 over the element, xi along \(x\) and eta along \(y\).

Notes

The entries are the derivatives \(\partial/\partial y = (2/l_{ey}) \partial/\partial \eta\) of those of BFSCCylinder.update_Su(), with the Hermite functions of the Notes of bfsccylinder.bfsccylinder.

update_Sv(xi, eta)#

Update Sv at a natural coordinate point

Sv is computed at xi, eta from lex and ley. Its product with an element displacement vector, e.g. u, gives the displacement \(v\) at that point, dimension of length. Only the entries at the degrees of freedom \(v\), \(v_{,x}\), \(v_{,y}\) of the nodes, the positions 3, 4, 5 plus 0, 10, 20 and 30, are written, the other entries keep the zeros of the constructor.

Parameters:

xi, etadouble: Natural coordinates of the point, dimensionless, from -1 to 1 over the element, xi along \(x\) and eta along \(y\).

Notes

The entries of node i, from 1 to 4, at the positions 10*(i - 1) + 3, 10*(i - 1) + 4 and 10*(i - 1) + 5, are the Hermite functions \(H_i\), \(H_{xi}\) and \(H_{yi}\) of the Notes of bfsccylinder.bfsccylinder evaluated at xi, eta, such that \(v = \boldsymbol{S}^v\boldsymbol{u}_e\).

update_Sv_x(xi, eta)#

Update Sv_x at a natural coordinate point

Sv_x is computed at xi, eta from lex and ley. Its product with an element displacement vector, e.g. u, gives the derivative \(v_{,x}\) at that point, dimensionless. Only the entries at the degrees of freedom \(v\), \(v_{,x}\), \(v_{,y}\) of the nodes, the positions 3, 4, 5 plus 0, 10, 20 and 30, are written, the other entries keep the zeros of the constructor.

Parameters:

xi, etadouble: Natural coordinates of the point, dimensionless, from -1 to 1 over the element, xi along \(x\) and eta along \(y\).

Notes

The entries are the derivatives \(\partial/\partial x = (2/l_{ex}) \partial/\partial \xi\) of those of BFSCCylinder.update_Sv(), with the Hermite functions of the Notes of bfsccylinder.bfsccylinder.

update_Sv_y(xi, eta)#

Update Sv_y at a natural coordinate point

Sv_y is computed at xi, eta from lex and ley. Its product with an element displacement vector, e.g. u, gives the derivative \(v_{,y}\) at that point, dimensionless. Only the entries at the degrees of freedom \(v\), \(v_{,x}\), \(v_{,y}\) of the nodes, the positions 3, 4, 5 plus 0, 10, 20 and 30, are written, the other entries keep the zeros of the constructor.

Parameters:

xi, etadouble: Natural coordinates of the point, dimensionless, from -1 to 1 over the element, xi along \(x\) and eta along \(y\).

Notes

The entries are the derivatives \(\partial/\partial y = (2/l_{ey}) \partial/\partial \eta\) of those of BFSCCylinder.update_Sv(), with the Hermite functions of the Notes of bfsccylinder.bfsccylinder.

update_Sw(xi, eta)#

Update Sw at a natural coordinate point

Sw is computed at xi, eta from lex and ley. Its product with an element displacement vector, e.g. u, gives the displacement \(w\) at that point, dimension of length. Only the entries at the degrees of freedom \(w\), \(w_{,x}\), \(w_{,y}\), \(w_{,xy}\) of the nodes, the positions 6, 7, 8, 9 plus 0, 10, 20 and 30, are written, the other entries keep the zeros of the constructor.

Parameters:

xi, etadouble: Natural coordinates of the point, dimensionless, from -1 to 1 over the element, xi along \(x\) and eta along \(y\).

Notes

The entries of node i, from 1 to 4, at the positions 10*(i - 1) + 6 to 10*(i - 1) + 9, are the Hermite functions \(H_i\), \(H_{xi}\), \(H_{yi}\) and \(H_{xyi}\) of the Notes of bfsccylinder.bfsccylinder evaluated at xi, eta, such that \(w = \boldsymbol{S}^w\boldsymbol{u}_e\).

update_Sw_x(xi, eta)#

Update Sw_x at a natural coordinate point

Sw_x is computed at xi, eta from lex and ley. Its product with an element displacement vector, e.g. u, gives the derivative \(w_{,x}\) at that point, dimensionless. Only the entries at the degrees of freedom \(w\), \(w_{,x}\), \(w_{,y}\), \(w_{,xy}\) of the nodes, the positions 6, 7, 8, 9 plus 0, 10, 20 and 30, are written, the other entries keep the zeros of the constructor.

Parameters:

xi, etadouble: Natural coordinates of the point, dimensionless, from -1 to 1 over the element, xi along \(x\) and eta along \(y\).

Notes

The entries are the derivatives \(\partial/\partial x = (2/l_{ex}) \partial/\partial \xi\) of those of BFSCCylinder.update_Sw(), with the Hermite functions of the Notes of bfsccylinder.bfsccylinder.

update_Sw_y(xi, eta)#

Update Sw_y at a natural coordinate point

Sw_y is computed at xi, eta from lex and ley. Its product with an element displacement vector, e.g. u, gives the derivative \(w_{,y}\) at that point, dimensionless. Only the entries at the degrees of freedom \(w\), \(w_{,x}\), \(w_{,y}\), \(w_{,xy}\) of the nodes, the positions 6, 7, 8, 9 plus 0, 10, 20 and 30, are written, the other entries keep the zeros of the constructor.

Parameters:

xi, etadouble: Natural coordinates of the point, dimensionless, from -1 to 1 over the element, xi along \(x\) and eta along \(y\).

Notes

The entries are the derivatives \(\partial/\partial y = (2/l_{ey}) \partial/\partial \eta\) of those of BFSCCylinder.update_Sw(), with the Hermite functions of the Notes of bfsccylinder.bfsccylinder.

bfsccylinder.bfsccylinder.update_KC0(shell, points, weights, KC0r, KC0c, KC0v)#

Update sparse vectors for linear constitutive stiffness matrix KC0

The element matrix, with 40*40 entries, is integrated with the integration points and weights given by points and weights, from lex, ley, R and the stiffnesses A11, ..., D66 of shell. All its entries are written: the global row and column indices to KC0r and KC0c, and the values added to KC0v, at the positions shell.init_k_KC0 to shell.init_k_KC0 + KC0_SPARSE_SIZE - 1. Before this function is called, c1, c2, c3, c4, init_k_KC0, lex, ley, R and the stiffnesses of shell must be set.

Parameters:

shellBFSCCylinder object: The cylindrical shell element being added to KC0
pointsnp.array: Integration poins 1-D array of DOUBLE of size nint, the nint given to the constructor of shell, with natural coordinates from -1 to 1, e.g. from bfsccylinder.quadrature.get_points_weights(). Used in both directions, for nint*nint integration points.
weightsnp.array: Integration weights 1-D array of DOUBLE of the same size as points.
KC0rnp.array: Array to store row positions of sparse values 1-D array of integers (INT) of size at least init_k_KC0 + KC0_SPARSE_SIZE.
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 DOUBLE, same size as KC0r. The values are added to the existing ones. Dimension of force per unit length between the displacements \(u\), \(v\), \(w\), multiplied by a length for each of the row and the column that is a first derivative, and by a length squared for each that is \(w_{,xy}\).

See also

update_KCNL, update_KG

Notes

KC0 is the Hessian of the strain energy of the Notes of bfsccylinder.bfsccylinder with the linear strains only, \(\boldsymbol{\varepsilon} = \boldsymbol{\varepsilon}_L\):

\[\boldsymbol{K}_{C0} = \int_A \left( \boldsymbol{B}_m^\top \boldsymbol{A} \boldsymbol{B}_m + \boldsymbol{B}_m^\top \boldsymbol{B} \boldsymbol{B}_b + \boldsymbol{B}_b^\top \boldsymbol{B} \boldsymbol{B}_m + \boldsymbol{B}_b^\top \boldsymbol{D} \boldsymbol{B}_b \right) dA\]

with the rows \(\boldsymbol{B}_m\) and \(\boldsymbol{B}_b\) of BFSCCylinder.update_Bm() and BFSCCylinder.update_Bb(), evaluated at each integration point as the local variables Bm1_*, Bm2_*, Bm3_*, Bb1_*, Bb2_* and Bb3_*, and the \(\boldsymbol{A}\), \(\boldsymbol{B}\) and \(\boldsymbol{D}\) of that point. KC0 does not depend on the displacements: it is the stiffness matrix of a linear static analysis and the tangent stiffness matrix at zero displacements.

bfsccylinder.bfsccylinder.update_KCNL(u, shell, points, weights, KCNLr, KCNLc, KCNLv)#

Update sparse vectors for nonlinear stiffness matrix KCNL

The element matrix is integrated with the integration points and weights given by points and weights, from lex, ley, R, the stiffnesses A11, ..., B66 of shell (the D stiffnesses are not used) and the displacements of the nodes of the element, copied from u to shell.u. The entries for which both the row and the column correspond to the in-plane degrees of freedom \(u\), \(u_{,x}\), \(u_{,y}\), \(v\), \(v_{,x}\), \(v_{,y}\) are not written. The global row and column indices of the other entries are written to KCNLr and KCNLc, and their values added to KCNLv, at the positions shell.init_k_KCNL to shell.init_k_KCNL + KCNL_SPARSE_SIZE - 1. Before this function is called, c1, c2, c3, c4, init_k_KCNL, lex, ley, R and the stiffnesses of shell must be set. KC0 + KCNL + KG, computed with update_KC0(), this function and update_KG() with the same u, points and weights, is the tangent stiffness matrix of the internal forces of update_fint().

Parameters:

unp.array: Global displacement vector 1-D array of DOUBLE with DOF entries per node of the model, in the order of the degrees of freedom of a node given in the module documentation. The entries of the nodes of shell are read at the positions c1+i, c2+i, c3+i and c4+i, for i from 0 to 9. Dimension of length for \(u\), \(v\), \(w\), dimensionless for the first derivatives and of one over length for \(w_{,xy}\).
shellBFSCCylinder object: The shell element being added to K
pointsnp.array: Integration poins 1-D array of DOUBLE of size nint, the nint given to the constructor of shell, with natural coordinates from -1 to 1, e.g. from bfsccylinder.quadrature.get_points_weights(). Used in both directions, for nint*nint integration points.
weightsnp.array: Integration weights 1-D array of DOUBLE of the same size as points.
KCNLrnp.array: Array to store row positions of sparse values 1-D array of integers (INT) of size at least init_k_KCNL + KCNL_SPARSE_SIZE.
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 DOUBLE, same size as KCNLr. The values are added to the existing ones, with the dimensions given for KC0v in update_KC0().

See also

update_KC0, update_KG, update_fint

Notes

Assuming that KCNL = KC0L + KCL0 + KCLL + KGNL, that is \(\boldsymbol{K}_{CNL} + \boldsymbol{K}_{GNL}\): the nonlinear constitutive stiffness matrix \(\boldsymbol{K}_{CNL} = \boldsymbol{K}_{C0L} + \boldsymbol{K}_{CL0} + \boldsymbol{K}_{CLL}\) of the terms of \(\boldsymbol{B}_{NL}\) and the geometric stiffness matrix \(\boldsymbol{K}_{GNL}\) of the membrane stress resultants of the nonlinear strains.

The first three groups are the constitutive terms coupling the linear and the nonlinear parts of the membrane strain. KGNL is the geometric term carrying the stress of the nonlinear membrane strain, A*BmL*u/2.

KGNL is geometric rather than constitutive, and would normally sit in update_KG. It is collected here so that update_KG can stay homogeneous of degree one in u, which is what a linear buckling analysis needs. With it here,

KT = KC0 + KCNL(u) + KG(u)

is the exact Jacobian of the vector assembled by update_fint, so a Newton-Raphson built on these three converges quadratically.

With the operators of the Notes of bfsccylinder.bfsccylinder, evaluated at each integration point with the slopes w_x and w_y of the element displacements shell.u:

\[\boldsymbol{K}_{C0L} = \int_A \left( \boldsymbol{B}_m^\top \boldsymbol{A} + \boldsymbol{B}_b^\top \boldsymbol{B} \right) \boldsymbol{B}_{NL} dA, \qquad \boldsymbol{K}_{CL0} = \int_A \boldsymbol{B}_{NL}^\top \left( \boldsymbol{A} \boldsymbol{B}_m + \boldsymbol{B} \boldsymbol{B}_b \right) dA\]

\[\begin{split}\boldsymbol{K}_{CLL} = \int_A \boldsymbol{B}_{NL}^\top \boldsymbol{A} \boldsymbol{B}_{NL} dA, \qquad \boldsymbol{K}_{GNL} = \int_A \boldsymbol{G}^\top \begin{bmatrix} N_{NLxx} & N_{NLxy} \\ N_{NLxy} & N_{NLyy} \end{bmatrix} \boldsymbol{G} dA\end{split}\]

with \(\{N_{NLxx}, N_{NLyy}, N_{NLxy}\}^\top = \boldsymbol{A}\boldsymbol{\varepsilon}_{NL} = \boldsymbol{A}\{w_{,x}^2/2, w_{,y}^2/2, w_{,x} w_{,y}\}^\top\). KC0L and KCL0 are linear in the displacements, KCLL and KGNL quadratic, and KCNL vanishes at zero displacements. Since \(\boldsymbol{B}_{NL}\) and \(\boldsymbol{G}\) have entries only at the degrees of freedom of \(w\), the block of the in-plane degrees of freedom is zero.

bfsccylinder.bfsccylinder.update_KCfollower(u, shell, points, weights, p, KCfollowerr, KCfollowerc, KCfollowerv, quadratic=0)#

Update sparse vectors for the load stiffness matrix KCfollower

The load stiffness matrix of a follower pressure, minus the derivative of the force vector of update_fext_follower() with respect to the element displacements, is integrated with the integration points and weights given by points and weights, from lex, ley and R of shell, the pressure p at the integration points and the displacements of the nodes of the element, copied from u to shell.u, which enter only with quadratic=1. The matrix is unsymmetric in general and is written in full, it must not be symmetrized. Only the entries that are not identically zero for quadratic=1 are written, for both values of quadratic: those for which the row corresponds to the degrees of freedom \(u\), \(u_{,x}\), \(u_{,y}\) and the column to those of \(v\) or \(w\), the row to those of \(v\) and the column to those of \(u\) or \(w\), and the row to those of \(w\), \(w_{,x}\), \(w_{,y}\), \(w_{,xy}\) and the column to any degree of freedom. The entries between the degrees of freedom of \(u\) and \(v\) are zero for quadratic=0. The global row and column indices are written to KCfollowerr and KCfollowerc, and the values added to KCfollowerv, at the positions shell.init_k_KCfollower to shell.init_k_KCfollower + KCfollower_SPARSE_SIZE - 1. Before this function is called, c1, c2, c3, c4, init_k_KCfollower, lex, ley and R of shell must be set; the stiffnesses are not used. With quadratic=0 the values do not depend on u.

Parameters:

unp.array: Global displacement vector 1-D array of DOUBLE with DOF entries per node of the model, in the order of the degrees of freedom of a node given in the module documentation. The entries of the nodes of shell are read at the positions c1+i, c2+i, c3+i and c4+i, for i from 0 to 9. Dimension of length for \(u\), \(v\), \(w\), dimensionless for the first derivatives and of one over length for \(w_{,xy}\).
shellBFSCCylinder object: The shell element being added to KCfollower
pointsnp.array: Integration poins 1-D array of DOUBLE of size nint, the nint given to the constructor of shell, with natural coordinates from -1 to 1, e.g. from bfsccylinder.quadrature.get_points_weights(). Used in both directions, for nint*nint integration points.
weightsnp.array: Integration weights 1-D array of DOUBLE of the same size as points.
pnp.array: Pressure at the integration points 2-D array of DOUBLE of shape (nint, nint), with the meaning, indices, units and sign given in update_fext_follower().
KCfollowerrnp.array: Array to store row positions of sparse values 1-D array of integers (INT) of size at least init_k_KCfollower + KCfollower_SPARSE_SIZE.
KCfollowercnp.array: Array to store column positions of sparse values Same type and size as KCfollowerr.
KCfollowervnp.array: Array to store sparse values 1-D array of DOUBLE, same size as KCfollowerr. The values are added to the existing ones, with the dimensions given for KC0v in update_KC0().
quadraticint, optional: 0, the default, or 1, the truncation of the area vector of update_fext_follower().

See also

update_fext_follower

Notes

With the force vector \(\boldsymbol{F}_p(\boldsymbol{u}_e)\) and the area vector \(\boldsymbol{a}\) of the Notes of update_fext_follower(), the load stiffness matrix is

\[\boldsymbol{K}_F = -\frac{\partial \boldsymbol{F}_p}{\partial \boldsymbol{u}_e} = -\int_A p \left( (\boldsymbol{S}^u)^\top \frac{\partial a_x}{\partial \boldsymbol{u}_e} + (\boldsymbol{S}^v)^\top \frac{\partial a_y}{\partial \boldsymbol{u}_e} + (\boldsymbol{S}^w)^\top \frac{\partial a_z}{\partial \boldsymbol{u}_e} \right) dA\]

where the derivatives of \(\boldsymbol{a}\) are those of its displacement gradients, e.g. \(\partial u_{,x}/\partial \boldsymbol{u}_e = \boldsymbol{S}^u_{,x}\), the local variables Su2_* in the code, with Su1_*, Su2_* and Su3_* the rows \(\boldsymbol{S}^u\), \(\boldsymbol{S}^u_{,x}\) and \(\boldsymbol{S}^u_{,y}\) at the integration point, and likewise Sv* and Sw*. The term \(\boldsymbol{a}_0\) does not contribute. With quadratic=0 the matrix is constant, the blocks being, per unit \(p\), \((\boldsymbol{S}^u)^\top \boldsymbol{S}^w_{,x}\) between \(u\) and \(w\), \((\boldsymbol{S}^v)^\top \boldsymbol{S}^w_{,y}\) between \(v\) and \(w\), and \(-(\boldsymbol{S}^w)^\top \boldsymbol{S}^u_{,x}\), \(-(\boldsymbol{S}^w)^\top \boldsymbol{S}^v_{,y}\), \(-(\boldsymbol{S}^w)^\top \boldsymbol{S}^w/R\) between \(w\) and \(u\), \(v\) and \(w\); with quadratic=1 it is linear in \(\boldsymbol{u}_e\).

The tangent stiffness matrix of the residual \(\boldsymbol{R} = \boldsymbol{F}_{int}(\boldsymbol{u}) - \boldsymbol{F}_p(\boldsymbol{u})\) is \(\boldsymbol{K}_{C0} + \boldsymbol{K}_{CNL} + \boldsymbol{K}_{GNL} + \boldsymbol{K}_G + \boldsymbol{K}_F\). For a linear buckling analysis under a pressure \(p = \lambda p_0\), with u the linear static solution of \(p_0\) and \(\boldsymbol{K}_G\) of update_KG(), the load multipliers and the buckling modes solve \((\boldsymbol{K}_{C0} + \lambda (\boldsymbol{K}_G + \boldsymbol{K}_{F_{ref}}))\boldsymbol{\phi} = \boldsymbol{0}\), with \(\boldsymbol{K}_{F_{ref}}\) the \(\boldsymbol{K}_F\) of \(p_0\) and quadratic=0, an eigenvalue problem that is unsymmetric in general. The conditions under which the assembled matrix is symmetric are given in the Notes of bfsccylinder.bfsccylinder; the matrix of one element is unsymmetric.

bfsccylinder.bfsccylinder.update_KG(u, shell, points, weights, KGr, KGc, KGv)#

Update sparse vectors for geometric stiffness matrix KG

The element matrix is integrated with the integration points and weights given by points and weights, from lex, ley and R of shell. At each integration point the membrane stress resultants Nxx, Nyy and Nxy are computed from the stiffnesses A11, ..., B66 of shell and from the displacements of the nodes of the element, copied from u to shell.u. Only the entries for which both the row and the column correspond to the degrees of freedom \(w\), \(w_{,x}\), \(w_{,y}\), \(w_{,xy}\) are written: the global row and column indices to KGr and KGc, and the values added to KGv, at the positions shell.init_k_KG to shell.init_k_KG + KG_SPARSE_SIZE - 1. Before this function is called, c1, c2, c3, c4, init_k_KG, lex, ley, R and the stiffnesses of shell must be set. Nxx, Nyy and Nxy are those of the linear strains, such that the values are linear in u; the corresponding term of the nonlinear strains is in update_KCNL().

Parameters:

unp.array: Global displacement vector 1-D array of DOUBLE with DOF entries per node of the model, in the order of the degrees of freedom of a node given in the module documentation. The entries of the nodes of shell are read at the positions c1+i, c2+i, c3+i and c4+i, for i from 0 to 9. Dimension of length for \(u\), \(v\), \(w\), dimensionless for the first derivatives and of one over length for \(w_{,xy}\).
shellBFSCCylinder object: The cylindrical shell element being added to KG
pointsnp.array: Integration poins 1-D array of DOUBLE of size nint, the nint given to the constructor of shell, with natural coordinates from -1 to 1, e.g. from bfsccylinder.quadrature.get_points_weights(). Used in both directions, for nint*nint integration points.
weightsnp.array: Integration weights 1-D array of DOUBLE of the same size as points.
KGrnp.array: Array to store row positions of sparse values 1-D array of integers (INT) of size at least init_k_KG + KG_SPARSE_SIZE.
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 DOUBLE, same size as KGr. The values are added to the existing ones, with the dimensions given for KC0v in update_KC0().

See also

update_KG_constant_stress, update_KCNL

Notes

Built from the membrane stress of the LINEAR strain alone, A*Bm*u + B*Bb*u, which keeps KG homogeneous of degree one in u. A linear buckling analysis needs exactly that, so that lambda*KG represents the stress state at load lambda*P.

The geometric term of the nonlinear membrane stress is therefore NOT here; update_KCNL carries it, so that KC0 + KCNL + KG is still the consistent tangent of update_fint.

With the operators of the Notes of bfsccylinder.bfsccylinder, at each integration point:

\[\begin{split}\boldsymbol{K}_G = \int_A \boldsymbol{G}^\top \begin{bmatrix} N_{xx} & N_{xy} \\ N_{xy} & N_{yy} \end{bmatrix} \boldsymbol{G} dA, \qquad \begin{Bmatrix} N_{xx} \\ N_{yy} \\ N_{xy} \end{Bmatrix} = \boldsymbol{A}\boldsymbol{B}_m\boldsymbol{u}_e + \boldsymbol{B}\boldsymbol{B}_b\boldsymbol{u}_e\end{split}\]

with \(\boldsymbol{G}\) the rows \(\boldsymbol{S}^w_{,x}\) and \(\boldsymbol{S}^w_{,y}\), G1 and G2, such that \(\boldsymbol{G}\boldsymbol{u}_e = \{w_{,x}, w_{,y}\}^\top\). It is the second derivative of \(\int_A \boldsymbol{N}^\top \boldsymbol{\varepsilon}_{NL} dA\) with respect to \(\boldsymbol{u}_e\) for \(\boldsymbol{N}\) held fixed. In a linear buckling analysis, with u the linear static solution for a reference load, the load multipliers \(\lambda\) and the buckling modes \(\boldsymbol{\phi}\) solve \((\boldsymbol{K}_{C0} + \lambda \boldsymbol{K}_G)\boldsymbol{\phi} = \boldsymbol{0}\), with \(\boldsymbol{K}_{C0}\) of update_KC0().

bfsccylinder.bfsccylinder.update_KG_constant_stress(shell, points, weights, KGr, KGc, KGv, Nxx, Nyy, Nxy)#

Update sparse vectors for geometric stiffness matrix KG

The same matrix as update_KG(), with the membrane stress resultants Nxx, Nyy and Nxy given, constant over the element, instead of computed from displacements. It is integrated with the integration points and weights given by points and weights, from lex and ley of shell; the stiffnesses and shell.u are not used. The entries are those of update_KG(): the global row and column indices are written to KGr and KGc, and the values added to KGv, at the positions shell.init_k_KG to shell.init_k_KG + KG_SPARSE_SIZE - 1. Before this function is called, c1, c2, c3, c4, init_k_KG, lex and ley of shell must be set.

Parameters:

shellBFSCCylinder object: The cylindrical shell element being added to KG
pointsnp.array: Integration poins 1-D array of DOUBLE of size nint, the nint given to the constructor of shell, with natural coordinates from -1 to 1, e.g. from bfsccylinder.quadrature.get_points_weights(). Used in both directions, for nint*nint integration points.
weightsnp.array: Integration weights 1-D array of DOUBLE of the same size as points.
KGrnp.array: Array to store row positions of sparse values 1-D array of integers (INT) of size at least init_k_KG + KG_SPARSE_SIZE.
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 DOUBLE, same size as KGr. The values are added to the existing ones, with the dimensions given for KC0v in update_KC0().
Nxx, Nyy, Nxyfloat: Constant membrane stress state Force per unit length, used at every integration point.

See also

update_KG

Notes

The matrix is the \(\boldsymbol{K}_G\) of the Notes of update_KG() with the given \(N_{xx}\), \(N_{yy}\) and \(N_{xy}\), and with the same rows \(\boldsymbol{G}\), which do not depend on \(R\) in the Donnell-type kinematics. Since it is linear in the given stress resultants, the load multipliers of a linear buckling analysis multiply the given membrane state, which is assumed instead of computed from a linear static solution, e.g. a uniform axial compression without the effects of the boundary conditions on the membrane state.

bfsccylinder.bfsccylinder.update_M(shell, Mr, Mc, Mv)#

Update sparse vectors for mass matrix M

The element mass matrix is computed from lex, ley, rho and h of shell with closed-form expressions, without integration points; R and the stiffnesses are not used. Only the entries that couple the degrees of freedom \(u\), \(u_{,x}\), \(u_{,y}\) with themselves, \(v\), \(v_{,x}\), \(v_{,y}\) with themselves, and \(w\), \(w_{,x}\), \(w_{,y}\), \(w_{,xy}\) with themselves are written: the global row and column indices to Mr and Mc, and the values added to Mv, at the positions shell.init_k_M to shell.init_k_M + M_SPARSE_SIZE - 1. Before this function is called, c1, c2, c3, c4, init_k_M, lex, ley, rho and h of shell must be set.

Parameters:

shellBFSCCylinder object: The cylindrical shell element being added to M
Mrnp.array: Array to store row positions of sparse values 1-D array of integers (INT) of size at least init_k_M + M_SPARSE_SIZE.
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 DOUBLE, same size as Mr. The values are added to the existing ones. Dimension of mass between the displacements \(u\), \(v\), \(w\), multiplied by a length for each of the row and the column that is a first derivative, and by a length squared for each that is \(w_{,xy}\).

Notes

The consistent mass matrix of the interpolation of the Notes of bfsccylinder.bfsccylinder, for a density \(\rho\) = rho uniform over the thickness \(h\) = h, with the Kirchhoff-Love displacements \(u - z w_{,x}\), \(v - z w_{,y}\) and \(w\) at the distance \(z\) from the mid-surface:

\[\boldsymbol{M} = \int_A \rho \left( h \left( (\boldsymbol{S}^u)^\top \boldsymbol{S}^u + (\boldsymbol{S}^v)^\top \boldsymbol{S}^v + (\boldsymbol{S}^w)^\top \boldsymbol{S}^w \right) + \frac{h^3}{12} \left( (\boldsymbol{S}^w_{,x})^\top \boldsymbol{S}^w_{,x} + (\boldsymbol{S}^w_{,y})^\top \boldsymbol{S}^w_{,y} \right) \right) dA\]

The terms in \(h^3/12\) are the rotary inertia of the rotations \(-w_{,x}\) and \(-w_{,y}\) of the Donnell-type kinematics. There is no coupling between the in-plane and the normal displacements, which corresponds to a reference surface at the middle of the thickness, and the variation of the area with \(z\) is neglected. The integrand is a polynomial and the integral is exact. The matrix of bfsccylinder.bfsccylinder_sanders.update_M() differs: its rotary inertia uses the Sanders-type rotation \(-w_{,y} + v/R\), which couples \(v\) and \(w\).

bfsccylinder.bfsccylinder.update_fext_follower(u, shell, points, weights, p, fext, quadratic=0, reference=1)#

Update global external force vector of a follower pressure

The force vector of a pressure that follows the deformation of the element, normal to its deformed mid-surface and acting on its deformed area, is integrated with the integration points and weights given by points and weights, from lex, ley and R of shell, the pressure p at the integration points and the displacements of the nodes of the element, copied from u to shell.u. It is added to fext at the positions c1+i, c2+i, c3+i and c4+i, for i from 0 to 9. Before this function is called, c1, c2, c3, c4, lex, ley and R of shell must be set; the stiffnesses are not used.

Parameters:

unp.array: Global displacement vector 1-D array of DOUBLE with DOF entries per node of the model, in the order of the degrees of freedom of a node given in the module documentation. The entries of the nodes of shell are read at the positions c1+i, c2+i, c3+i and c4+i, for i from 0 to 9. Dimension of length for \(u\), \(v\), \(w\), dimensionless for the first derivatives and of one over length for \(w_{,xy}\).
shellBFSCCylinder object: The shell element being added to fext
pointsnp.array: Integration poins 1-D array of DOUBLE of size nint, the nint given to the constructor of shell, with natural coordinates from -1 to 1, e.g. from bfsccylinder.quadrature.get_points_weights(). Used in both directions, for nint*nint integration points.
weightsnp.array: Integration weights 1-D array of DOUBLE of the same size as points.
pnp.array: Pressure at the integration points 2-D array of DOUBLE of shape (nint, nint), the entry [i, j] being used at the integration point xi = points[i], eta = points[j], as the entries of A11 of shell. Force per unit area of the undeformed mid-surface, already multiplied by the load factor if there is one. A positive pressure pushes along \(+z\), outwards, away from the axis of the cylinder, the direction of a positive \(w\), such that an external pressure is negative.
fextnp.array: Global external force vector 1-D array of DOUBLE of the same size as u. The values are added to the existing ones. Dimension of force for \(u\), \(v\), \(w\), multiplied by a length for the first derivatives and by a length squared for \(w_{,xy}\).
quadraticint, optional: 0, the default, to truncate the area vector of the deformed mid-surface to its terms of order 0 and 1 in the displacement gradients, such that the force vector is affine in u; 1 to also include its terms of order 2, see the Notes.
referenceint, optional: 1, the default, to include the term of order 0 of the area vector, which gives the dead pressure load of the undeformed configuration \(\boldsymbol{F}_p(\boldsymbol{0})\); 0 to leave it out, which gives \(\boldsymbol{F}_p(\boldsymbol{u}_e) - \boldsymbol{F}_p(\boldsymbol{0})\).

See also

update_KCfollower

Notes

With \(\boldsymbol{X} = \boldsymbol{X}_0 + \boldsymbol{u}\) the position of the deformed mid-surface, in the basis of the undeformed mid-surface \(\boldsymbol{e}_x\), \(\boldsymbol{e}_y\) and its normal \(\boldsymbol{e}_z\), the virtual work of the pressure is

\[\delta W_p = \int_A p \boldsymbol{a}^\top \delta \boldsymbol{u} dA, \qquad \boldsymbol{a} = \boldsymbol{X}_{,x} \times \boldsymbol{X}_{,y}\]

where \(\boldsymbol{a}\) is the area vector of the deformed mid-surface per unit undeformed area, \(\delta \boldsymbol{u} = \{\delta u, \delta v, \delta w\}^\top\) and, with \(\boldsymbol{e}_{y,y} = -\boldsymbol{e}_z/R\) and \(\boldsymbol{e}_{z,y} = \boldsymbol{e}_y/R\),

\[\begin{split}\boldsymbol{X}_{,x} = \begin{Bmatrix} 1 + u_{,x} \\ v_{,x} \\ w_{,x} \end{Bmatrix}, \qquad \boldsymbol{X}_{,y} = \begin{Bmatrix} u_{,y} \\ 1 + v_{,y} + w/R \\ w_{,y} - k_v v \end{Bmatrix}\end{split}\]

with \(k_v = 0\) for the Donnell-type kinematics, which neglect the term \(v/R\) of the rotation of the normal as their strains do, and \(k_v = 1/R\) for the Sanders-type kinematics, for which \(\boldsymbol{X}_{,x} \times \boldsymbol{X}_{,y}\) is the exact area vector. Its parts of order 0, 1 and 2 in the displacement gradients are

\[\begin{split}\boldsymbol{a}_0 = \begin{Bmatrix} 0 \\ 0 \\ 1 \end{Bmatrix}, \qquad \boldsymbol{a}_1 = \begin{Bmatrix} -w_{,x} \\ -w_{,y} + k_v v \\ u_{,x} + v_{,y} + w/R \end{Bmatrix}\end{split}\]

\[\begin{split}\boldsymbol{a}_2 = \begin{Bmatrix} v_{,x} (w_{,y} - k_v v) - w_{,x} (v_{,y} + w/R) \\ w_{,x} u_{,y} - u_{,x} (w_{,y} - k_v v) \\ u_{,x} (v_{,y} + w/R) - v_{,x} u_{,y} \end{Bmatrix}\end{split}\]

\(\boldsymbol{a}_0\) gives the dead load, \(\boldsymbol{a}_1\) the rotation of the normal and the change of area. The force vector added to fext is

\[\boldsymbol{F}_p(\boldsymbol{u}_e) = \int_A p \left( (\boldsymbol{S}^u)^\top a_x + (\boldsymbol{S}^v)^\top a_y + (\boldsymbol{S}^w)^\top a_z \right) dA, \qquad \boldsymbol{a} = r \boldsymbol{a}_0 + \boldsymbol{a}_1 + q \boldsymbol{a}_2\]

with \(r\) = reference and \(q\) = quadratic, \(\boldsymbol{a}\) being ax, ay and az in the code, evaluated at each integration point with the displacement gradients u_x, u_y, v_x, v_y, w, w_x and w_y of u. With the ordering of the moderate rotations of the strains of the Notes of bfsccylinder.bfsccylinder, \(u_{,x} \sim \epsilon\) and \(w_{,x} \sim \epsilon^{1/2}\), \(\boldsymbol{a}_0 + \boldsymbol{a}_1\) retains every term of \(p \boldsymbol{a}^\top \delta \boldsymbol{u}\) up to the order \(\epsilon p \delta w\) of the strains, and \(\boldsymbol{a}_2\) adds terms of order \(\epsilon^2 p \delta w\), of the order of those that the strains neglect. quadratic=0 is therefore the consistent truncation, for which \(\boldsymbol{F}_p\) is affine in \(\boldsymbol{u}_e\), and quadratic=1 gives the complete area vector of the kinematics, e.g. to check the sensitivity of a result to the truncation. At \(\boldsymbol{u}_e = \boldsymbol{0}\) the force vector is the dead load \(\int_A p (\boldsymbol{S}^w)^\top dA\) for both truncations.

In a nonlinear analysis the follower pressure enters the residual \(\boldsymbol{R} = \boldsymbol{F}_{int}(\boldsymbol{u}) - \boldsymbol{F}_p(\boldsymbol{u})\), with \(\boldsymbol{F}_{int}\) of update_fint(), whose Jacobian is the sum of the matrices of update_KC0(), update_KCNL(), update_KG() and update_KCfollower(). For a pressure \(p = \lambda p_0\) proportional to a load factor \(\lambda\), reference=0 gives the part \(\lambda (\boldsymbol{F}_p(\boldsymbol{u}) - \boldsymbol{F}_p(\boldsymbol{0}))\) that depends on the displacements, to be kept with the internal forces, while the dead part \(\lambda \boldsymbol{F}_p(\boldsymbol{0})\) can be given as the external force vector of the load factor.

bfsccylinder.bfsccylinder.update_fint(u, shell, points, weights, fint)#

Update global internal force vector

The internal forces of the element are integrated with the integration points and weights given by points and weights, from lex, ley, R, the stiffnesses A11, ..., D66 of shell and the displacements of the nodes of the element, copied from u to shell.u, including the terms nonlinear in the displacements. They are added to fint at the positions c1+i, c2+i, c3+i and c4+i, for i from 0 to 9. Before this function is called, c1, c2, c3, c4, lex, ley, R and the stiffnesses of shell must be set.

Parameters:

unp.array: Global displacement vector 1-D array of DOUBLE with DOF entries per node of the model, in the order of the degrees of freedom of a node given in the module documentation. The entries of the nodes of shell are read at the positions c1+i, c2+i, c3+i and c4+i, for i from 0 to 9. Dimension of length for \(u\), \(v\), \(w\), dimensionless for the first derivatives and of one over length for \(w_{,xy}\).
shellBFSCCylinder object: The shell element being added to fint
pointsnp.array: Integration poins 1-D array of DOUBLE of size nint, the nint given to the constructor of shell, with natural coordinates from -1 to 1, e.g. from bfsccylinder.quadrature.get_points_weights(). Used in both directions, for nint*nint integration points.
weightsnp.array: Integration weights 1-D array of DOUBLE of the same size as points.
fintnp.array: Global internal force array 1-D array of DOUBLE of the same size as u. The values are added to the existing ones. Dimension of force for \(u\), \(v\), \(w\), multiplied by a length for the first derivatives and by a length squared for \(w_{,xy}\).

See also

update_KCNL

Notes

The internal forces are the gradient of the strain energy of the Notes of bfsccylinder.bfsccylinder:

\[\boldsymbol{F}_{int} = \int_A \left( \left( \boldsymbol{B}_m + \boldsymbol{B}_{NL} \right)^\top \boldsymbol{N} + \boldsymbol{B}_b^\top \boldsymbol{M} \right) dA\]

with the stress resultants of the total membrane strains,

\[\boldsymbol{N} = \boldsymbol{A} \left( \boldsymbol{B}_m + \frac{1}{2} \boldsymbol{B}_{NL} \right) \boldsymbol{u}_e + \boldsymbol{B}\boldsymbol{B}_b\boldsymbol{u}_e, \qquad \boldsymbol{M} = \boldsymbol{B} \left( \boldsymbol{B}_m + \frac{1}{2} \boldsymbol{B}_{NL} \right) \boldsymbol{u}_e + \boldsymbol{D}\boldsymbol{B}_b\boldsymbol{u}_e\]

Nxx, Nyy, Nxy and Mxx, Myy, Mxy in the code, evaluated at each integration point with the slopes w_x and w_y. The tangent stiffness matrix \(\partial \boldsymbol{F}_{int}/\partial \boldsymbol{u}_e\) is the sum of the matrices of update_KC0(), update_KCNL() and update_KG() computed with the same u, points and weights. At equilibrium \(\boldsymbol{F}_{int}\) balances the external forces.

Donnell-type cylinder element (bfsccylinder.bfsccylinder)#

Notes#

References#

Donnell-type cylinder element (`bfsccylinder.bfsccylinder`)#