BFSC plate element (`bfscplate2d.bfscplate2d`)#

BFSC plate element and the functions computing its matrices

The module has the element class BFSCPlate2D and the functions update_KC0(), update_KCNL(), update_KG(), update_KG_cte_N(), update_M() and update_KCfollower(), which write the contributions of one element in sparse row, column and value arrays, from the positions given by the init_k_* attributes of the element, update_fint(), which adds the internal forces of one element to a global internal force vector, and update_fext_follower(), which adds the force vector of a follower pressure on one element to a global force vector.

The element is a rectangle with four nodes, of length lex along \(x\) and ley along \(y\). The natural coordinates xi and eta go from -1 to 1 over the element, node 1 being at xi = eta = -1, node 2 at xi = 1, eta = -1, node 3 at xi = eta = 1 and node 4 at xi = -1, eta = 1. Each node has 10 degrees of freedom, in the order \(u\), \(u_{,x}\), \(u_{,y}\), \(v\), \(v_{,x}\), \(v_{,y}\), \(w\), \(w_{,x}\), \(w_{,y}\), \(w_{,xy}\): the displacements \(u\), \(v\), \(w\) along \(x\), \(y\), \(z\), with dimension of length, their first derivatives with respect to \(x\) and \(y\), dimensionless, and the second derivative of \(w\) with respect to \(x\) and \(y\), with dimension of one over length. The degrees of freedom of nodes 1, 2, 3 and 4 are at the positions c1 to c1 + 9, c2 to c2 + 9, c3 to c3 + 9 and c4 to c4 + 9 of the global vectors and matrices, c1 to c4 being attributes of the element.

The functions other than update_M() integrate numerically over the element with nint by nint points, the products of the one-dimensional points and weights given to them, e.g. from bfscplate2d.quadrature.get_points_weights(). update_M() computes its entries with explicit expressions, without numerical integration. The functions evaluate the interpolation functions themselves: they do not read or change the attributes Bm, Bb, Su, …, Sw_y of BFSCPlate2D, which are updated only by the methods of the element.

The module also defines:

KC0_SPARSE_SIZE = 1600, KCNL_SPARSE_SIZE = 1024, KG_SPARSE_SIZE = 256, M_SPARSE_SIZE = 544 and KCfollower_SPARSE_SIZE = 1056: number of entries that one element writes in the row, column and value arrays of KC0, KCNL, KG, M and KCfollower, starting at the position given by its attribute init_k_KC0, init_k_KCNL, init_k_KG, init_k_M or init_k_KCfollower. For n elements the arrays need n times these sizes;
DOUBLE = np.float64: NumPy type of the arrays of values.

Notes#

The element is the rectangular plate element of Bogner, Fox and Schmit (BFS) [1], whose deflection \(w\) is interpolated with products of one-dimensional cubic Hermite polynomials, with \(w\), \(w_{,x}\), \(w_{,y}\) and \(w_{,xy}\) as nodal degrees of freedom. In the BFSC element [2] the in-plane displacements \(u\) and \(v\) are interpolated with the same polynomials, with \(u\), \(u_{,x}\), \(u_{,y}\) and \(v\), \(v_{,x}\), \(v_{,y}\) as nodal degrees of freedom, without the cross derivatives \(u_{,xy}\) and \(v_{,xy}\).

Interpolation

The element is a rectangle with sides parallel to \(x\) and \(y\), such that \(x = x_1 + lex (\xi + 1)/2\) and \(y = y_1 + ley (\eta + 1)/2\), with \((x_1, y_1)\) the position of node 1, and

\[\frac{\partial}{\partial x} = \frac{2}{lex} \frac{\partial}{\partial \xi}, \qquad \frac{\partial}{\partial y} = \frac{2}{ley} \frac{\partial}{\partial \eta}, \qquad dA = dx \, dy = \frac{lex \, ley}{4} d\xi \, d\eta\]

For node \(i\) at \((\xi_i, \eta_i)\), with \(\xi_i, \eta_i = \pm 1\), the interpolation functions are [3]

\[ \begin{align}\begin{aligned}H_i = \frac{1}{16} (\xi + \xi_i)^2 (\xi \xi_i - 2) (\eta + \eta_i)^2 (\eta \eta_i - 2)\\H_{x,i} = -\frac{lex}{32} \xi_i (\xi + \xi_i)^2 (\xi \xi_i - 1) (\eta + \eta_i)^2 (\eta \eta_i - 2)\\H_{y,i} = -\frac{ley}{32} (\xi + \xi_i)^2 (\xi \xi_i - 2) \eta_i (\eta + \eta_i)^2 (\eta \eta_i - 1)\\H_{xy,i} = \frac{lex \, ley}{64} \xi_i (\xi + \xi_i)^2 (\xi \xi_i - 1) \eta_i (\eta + \eta_i)^2 (\eta \eta_i - 1)\end{aligned}\end{align} \]

and the displacements are

\[ \begin{align}\begin{aligned}u = \sum_{i=1}^{4} \left( H_i u_i + H_{x,i} u_{,x i} + H_{y,i} u_{,y i} \right), \qquad v = \sum_{i=1}^{4} \left( H_i v_i + H_{x,i} v_{,x i} + H_{y,i} v_{,y i} \right)\\w = \sum_{i=1}^{4} \left( H_i w_i + H_{x,i} w_{,x i} + H_{y,i} w_{,y i} + H_{xy,i} w_{,xy i} \right)\end{aligned}\end{align} \]

i.e. \(u = \boldsymbol{S}^u\boldsymbol{u}_e\), \(v = \boldsymbol{S}^v\boldsymbol{u}_e\) and \(w = \boldsymbol{S}^w\boldsymbol{u}_e\), with the rows Su, Sv and Sw of BFSCPlate2D and the 40 element degrees of freedom \(\boldsymbol{u}_e\), ue in the code. Each function gives the value 1 to the nodal degree of freedom it multiplies and 0 to the other nodal degrees of freedom of the element, e.g. \(\partial H_{x,i}/\partial x = 1\) at node \(i\). The functions are bicubic, of degree 3 in \(\xi\) and in \(\eta\). Without the functions \(H_{xy,i}\) in \(u\) and \(v\), the interpolation gives \(u_{,xy} = v_{,xy} = 0\) at the nodes.

Kinematics

The kinematics is that of the classical laminated plate theory, without transverse shear deformation, with the nonlinear strains of von Karman [4]: the displacements at the distance \(z\) from the reference surface are \(u - z w_{,x}\), \(v - z w_{,y}\) and \(w\), and, for small strains and moderate rotations, only the squares and products of the slopes \(w_{,x}\) and \(w_{,y}\) are kept in the membrane strains:

\[\begin{split}\boldsymbol{\epsilon} = \begin{Bmatrix} \epsilon_{xx} \\ \epsilon_{yy} \\ \gamma_{xy} \end{Bmatrix} = \begin{Bmatrix} u_{,x} + w_{,x}^2/2 \\ v_{,y} + w_{,y}^2/2 \\ u_{,y} + v_{,x} + w_{,x} w_{,y} \end{Bmatrix}, \qquad \boldsymbol{\kappa} = \begin{Bmatrix} \kappa_{xx} \\ \kappa_{yy} \\ \kappa_{xy} \end{Bmatrix} = \begin{Bmatrix} -w_{,xx} \\ -w_{,yy} \\ -2 w_{,xy} \end{Bmatrix}\end{split}\]

In terms of the element degrees of freedom

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

with the operators of the code: \(\boldsymbol{B}_m\) (Bm), with the rows \(\boldsymbol{S}^u_{,x}\), \(\boldsymbol{S}^v_{,y}\) and \(\boldsymbol{S}^u_{,y} + \boldsymbol{S}^v_{,x}\); \(\boldsymbol{B}_b\) (Bb), with the rows \(-\boldsymbol{S}^w_{,xx}\), \(-\boldsymbol{S}^w_{,yy}\) and \(-2\boldsymbol{S}^w_{,xy}\); \(\boldsymbol{G}\) (G), with the rows \(\boldsymbol{S}^w_{,x}\) and \(\boldsymbol{S}^w_{,y}\), which give the slopes \(w_{,x} = \boldsymbol{S}^w_{,x}\boldsymbol{u}_e\) and \(w_{,y} = \boldsymbol{S}^w_{,y}\boldsymbol{u}_e\) (w_x and w_y); and \(\boldsymbol{B}_{NL}\) (BmL), with the rows \(w_{,x}\boldsymbol{S}^w_{,x}\), \(w_{,y}\boldsymbol{S}^w_{,y}\) and \(w_{,x}\boldsymbol{S}^w_{,y} + w_{,y}\boldsymbol{S}^w_{,x}\). In the functions of the module the nonzero entries of these operators are the local variables Bm1_01, …, BmL3_40, the digit after the name being the row, from 1, and the two last digits the column, from 1. Since the nonlinear strains are quadratic in \(\boldsymbol{u}_e\), \(\boldsymbol{B}_{NL}\boldsymbol{u}_e\) is twice the nonlinear strains, whereas the variation of the strains is

\[\delta\boldsymbol{\epsilon} = \left( \boldsymbol{B}_m + \boldsymbol{B}_{NL} \right) \delta\boldsymbol{u}_e\]

\(\boldsymbol{B}_{NL}\) and \(\boldsymbol{G}\) have nonzero entries only in the columns of \(w\), \(w_{,x}\), \(w_{,y}\) and \(w_{,xy}\).

Constitutive relation and strain energy

The membrane forces \(\boldsymbol{N} = \{Nxx, Nyy, Nxy\}^\top\) and the moments \(\boldsymbol{M} = \{Mxx, Myy, Mxy\}^\top\) are

\[\boldsymbol{N} = \boldsymbol{A}\boldsymbol{\epsilon} + \boldsymbol{B}\boldsymbol{\kappa}, \qquad \boldsymbol{M} = \boldsymbol{B}\boldsymbol{\epsilon} + \boldsymbol{D}\boldsymbol{\kappa}\]

with the symmetric blocks of ABD, \(\boldsymbol{A}\) of A11, A12, A16, A22, A26, A66, and likewise \(\boldsymbol{B}\) and \(\boldsymbol{D}\), constant over the element, which assumes a linear elastic laminate whose stiffness does not vary over the element. The strain energy of the element is

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

Element matrices and internal forces

The internal force vector is the gradient of \(U\), update_fint():

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

and its Jacobian, the tangent stiffness matrix, is split as

\[\frac{\partial \boldsymbol{F}_{int}}{\partial \boldsymbol{u}_e} = \boldsymbol{K}_{C0} + \boldsymbol{K}_{CNL} + \boldsymbol{K}_{GNL} + \boldsymbol{K}_G\]

with the linear constitutive stiffness matrix of update_KC0(), independent of \(\boldsymbol{u}_e\), the nonlinear constitutive stiffness matrix \(\boldsymbol{K}_{CNL}\) of the terms of \(\boldsymbol{B}_{NL}\) and the geometric stiffness matrix \(\boldsymbol{K}_{GNL}\) of the membrane forces of the nonlinear strains, whose sum is the matrix of update_KCNL(), null for null slopes, and the geometric stiffness matrix \(\boldsymbol{K}_G\) of the membrane forces of the linear strains, of update_KG(), linear in \(\boldsymbol{u}_e\). For small displacements \(\boldsymbol{F}_{int} \approx \boldsymbol{K}_{C0}\boldsymbol{u}_e\), the relation of a linear static analysis; update_KG() and update_KG_cte_N() give the geometric stiffness matrix of a linear buckling analysis. The mass matrix of update_M() follows from the kinetic energy of the same kinematics.

Follower pressure

A pressure \(p\) per unit undeformed area that stays normal to the deformed mid-surface and acts on its deformed area has the virtual work \(\int_A p \, \boldsymbol{a} \cdot \delta\boldsymbol{u} \, dA\), with the area vector \(\boldsymbol{a} = \boldsymbol{X}_{,x} \times \boldsymbol{X}_{,y}\) of the deformed position \(\boldsymbol{X}\) of the mid-surface [5] [6]. update_fext_follower() gives its force vector \(\boldsymbol{F}_p(\boldsymbol{u})\), with \(\boldsymbol{a}\) truncated to first order in the displacement gradients, the truncation consistent with the moderate rotations of the von Karman strains, or complete, and update_KCfollower() the load stiffness matrix \(\boldsymbol{K}_F = -\partial \boldsymbol{F}_p / \partial \boldsymbol{u}_e\), constant for the first-order area vector. A positive \(p\) pushes along \(+z\), and at \(\boldsymbol{u} = \boldsymbol{0}\) the force vector is the dead pressure load. The tangent stiffness matrix of the residual \(\boldsymbol{R} = \boldsymbol{F}_{int} - \lambda \boldsymbol{F}_p(\boldsymbol{u})\) is \(\boldsymbol{K}_{C0} + \boldsymbol{K}_{CNL} + \boldsymbol{K}_{GNL} + \boldsymbol{K}_G + \boldsymbol{K}_F\), with \(\boldsymbol{K}_F = \lambda \boldsymbol{K}_{F_{ref}}\) and \(\boldsymbol{K}_{F_{ref}}\) the load stiffness matrix of \(p\). The load stiffness matrix is unsymmetric in general, a non-conservative load. For the first-order area vector it is symmetric, and the follower pressure conservative, if and only if \(p\) is uniform over the loaded region and, at each point of the boundary of the loaded region, either \(w = 0\) or the in-plane displacement normal to the boundary is zero, see the Notes of update_KCfollower().

Numerical integration

The integrals over the element are computed as

\[\int_A f \, dA = \frac{lex \, ley}{4} \sum_{i=1}^{nint} \sum_{j=1}^{nint} weights_i \, weights_j \, f(points_i, points_j)\]

where 0.25*lex*ley*weight of the code is \(lex \, ley \, weights_i \, weights_j / 4\). With the Gauss-Legendre points of bfscplate2d.quadrature.get_points_weights(), which integrate exactly the polynomials of degree up to \(2 \, nint - 1\) in each of \(\xi\) and \(\eta\), and since the interpolation functions are bicubic, the integration is exact, for a constant ABD, with nint >= 4 for KC0 and for update_KG_cte_N(), whose integrands have a degree of at most 6 in \(\xi\) and in \(\eta\), with nint >= 5 for update_KG(), of degree at most 9, and with nint >= 7 for KCNL and the internal forces, of degree at most 12. With fewer points, KC0 + KCNL + KG is still the exact Jacobian of the internal forces computed with the same points. The force vector and the load stiffness matrix of a follower pressure constant over the element are integrated exactly with nint >= 4 for the first-order area vector, degree at most 6, and with nint >= 5 for the complete one, degree at most 8.

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]

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.

[3]

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

[4]

Reddy, J. N., Mechanics of Laminated Composite Plates and Shells: Theory and Analysis, 2nd ed., CRC Press, Boca Raton, 2003, doi:10.1201/b12409.

[5]

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.

[6]

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 bfscplate2d.bfscplate2d.BFSCPlate2D#

Four-node rectangular BFSC plate element

Nodal connectivity for plate element similar to Nastran’s CQUAD4:

 ^ y
 |

4 ________ 3
 |       |
 |       |   --> x
 |       |
 |_______|
1         2

The constructor takes no parameters. It sets n1, n2, n3, n4 and c1, c2, c3, c4 to -1, init_k_KC0, init_k_KCNL, init_k_KG, init_k_M and init_k_KCfollower to 0, lex, ley, h and rho to 0, and allocates ABD, Bm, Bb, Su to Sw_y and u filled with zeros. Before calling the functions of the module, the caller sets c1, c2, c3, c4, lex, ley, ABD and the init_k_* positions of the matrices to be computed, and h and rho for update_M().

Attributes:

n1, n2, n3, n4int: Node identification numbers. Not read by the functions of the 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, e.g. DOF times the index of the node. The 10 degrees of freedom of node 1 are at the positions c1 to c1 + 9, and likewise for the other nodes.
init_k_KC0, init_k_KCNL, init_k_KG, init_k_M, init_k_KCfollowerint: Position of the first entry written by the element in the sparse row, column and value arrays of KC0 (update_KC0()), KCNL (update_KCNL()), KG (update_KG() and update_KG_cte_N()), M (update_M()) and KCfollower (update_KCfollower()).
lex, leydouble: Lengths of the element sides along \(x\), from node 1 to node 2, and along \(y\), from node 1 to node 4. Dimension of length.
rhodouble: Density, mass per unit volume. Read only by update_M().
hdouble: Thickness, dimension of length. Read only by update_M().
ABDdouble[:, ::1]: Constitutive matrix [[A, B], [B, D]] of shape (6, 6), e.g. the ABD attribute of a composites.core.Laminate created with composites.isotropic_plate() or composites.laminated_plate(). A has dimension of force per unit length, B of force and D of force times length. Only the entries on and above the diagonal of the blocks ABD[:3, :3] (A11, A12, A16, A22, A26, A66), ABD[:3, 3:] (B11, B12, B16, B22, B26, B66) and ABD[3:, 3:] (D11, ..., D66) are read.
Bmdouble[:, ::1]: Array of shape (3, 40), updated by update_Bm(). The products of its rows with the element degrees of freedom give \(u_{,x}\), \(v_{,y}\) and \(u_{,y} + v_{,x}\), dimensionless.
Bbdouble[:, ::1]: Array of shape (3, 40), updated by update_Bb(). The products of its rows with the element degrees of freedom give \(-w_{,xx}\), \(-w_{,yy}\) and \(-2 w_{,xy}\), with dimension of one over length.
Su, Sv, Swdouble[::1]: Arrays of size 40, updated by update_Su(), update_Sv() and update_Sw(). Their products with the element degrees of freedom give \(u\), \(v\) and \(w\), dimension of length.
Su_x, Su_y, Sv_x, Sv_y, Sw_x, Sw_ydouble[::1]: Arrays of size 40, updated by the methods of the same name, update_Su_x() to update_Sw_y(). Their products with the element degrees of freedom give \(u_{,x}\), \(u_{,y}\), \(v_{,x}\), \(v_{,y}\), \(w_{,x}\) and \(w_{,y}\), dimensionless.
udouble[::1]: Array of size 40 with the element degrees of freedom, the 10 of node 1 followed by those of nodes 2, 3 and 4, with the dimensions given in the module documentation. Overwritten with values read from the global displacement vector by update_KCNL(), update_KG(), update_fint(), update_fext_follower() and update_KCfollower().

Notes

The rows Su to Sw_y are the interpolation functions \(H_i\), \(H_{x,i}\), \(H_{y,i}\), \(H_{xy,i}\) of the Notes of bfscplate2d.bfscplate2d, and their derivatives with respect to \(x\) and \(y\), at the positions of the degrees of freedom they multiply. Bm and Bb are the operators \(\boldsymbol{B}_m\) and \(\boldsymbol{B}_b\) of the same Notes, which give the linear membrane strains \(\{u_{,x}, v_{,y}, u_{,y} + v_{,x}\}^\top\) and the curvatures \(\boldsymbol{\kappa} = \{-w_{,xx}, -w_{,yy}, -2 w_{,xy}\}^\top\) of the classical laminated plate theory. ABD gives the membrane forces \(\boldsymbol{N} = \boldsymbol{A}\boldsymbol{\epsilon} + \boldsymbol{B}\boldsymbol{\kappa}\) and the moments \(\boldsymbol{M} = \boldsymbol{B}\boldsymbol{\epsilon} + \boldsymbol{D}\boldsymbol{\kappa}\) from the membrane strains \(\boldsymbol{\epsilon}\) and the curvatures, and h and rho the mass per unit area \(\rho h\) and the rotary inertia \(\rho h^3/12\) of update_M().

update_Bb(xi, eta)#

Update the bending rows Bb at a natural coordinate point

The entries of the three rows of Bb in the columns of the out-of-plane degrees of freedom \(w\), \(w_{,x}\), \(w_{,y}\), \(w_{,xy}\) of the four nodes are computed from lex and ley. The other entries are not changed (zero after construction).

Parameters:

xi, etadouble: Natural coordinates of the evaluation point. Dimensionless, from -1 to 1 over the element.

Notes

Bb is the operator \(\boldsymbol{B}_b\) of the curvatures \(\boldsymbol{\kappa} = \{-w_{,xx}, -w_{,yy}, -2 w_{,xy}\}^\top\) of the Notes of bfscplate2d.bfscplate2d, with the rows

\[-\frac{4}{lex^2} \frac{\partial^2 \boldsymbol{S}^w}{\partial \xi^2}, \qquad -\frac{4}{ley^2} \frac{\partial^2 \boldsymbol{S}^w}{\partial \eta^2}, \qquad -\frac{8}{lex \, ley} \frac{\partial^2 \boldsymbol{S}^w}{\partial \xi \partial \eta}\]

The curvatures are those of the rotations \(-w_{,x}\) and \(-w_{,y}\) of the normal to the reference surface.

update_Bm(xi, eta)#

Update the membrane rows Bm at a natural coordinate point

The entries of Bm in the columns of the in-plane degrees of freedom are computed from lex and ley: row 0 in the columns of \(u\), \(u_{,x}\), \(u_{,y}\), row 1 in those of \(v\), \(v_{,x}\), \(v_{,y}\) and row 2 in both, of the four nodes. The other entries are not changed (zero after construction).

Parameters:

xi, etadouble: Natural coordinates of the evaluation point. Dimensionless, from -1 to 1 over the element.

Notes

Bm is the operator \(\boldsymbol{B}_m\) of the linear membrane strains of the Notes of bfscplate2d.bfscplate2d, with the rows

\[\frac{2}{lex} \frac{\partial \boldsymbol{S}^u}{\partial \xi}, \qquad \frac{2}{ley} \frac{\partial \boldsymbol{S}^v}{\partial \eta}, \qquad \frac{2}{ley} \frac{\partial \boldsymbol{S}^u}{\partial \eta} + \frac{2}{lex} \frac{\partial \boldsymbol{S}^v}{\partial \xi}\]

i.e. Su_x, Sv_y and Su_y + Sv_x.

update_Su(xi, eta)#

Update the row Su at a natural coordinate point

The entries of Su at the positions of the degrees of freedom \(u\), \(u_{,x}\), \(u_{,y}\) of the four nodes are computed from lex and ley; the product of Su with the element degrees of freedom gives \(u\). The other entries are not changed (zero after construction).

Parameters:

xi, etadouble: Natural coordinates of the evaluation point. Dimensionless, from -1 to 1 over the element.

Notes

Su is the row \(\boldsymbol{S}^u\) of the interpolation of \(u\) of the Notes of bfscplate2d.bfscplate2d, with the functions \(H_i\), \(H_{x,i}\) and \(H_{y,i}\) of each node \(i\) at the positions of its degrees of freedom \(u\), \(u_{,x}\) and \(u_{,y}\).

update_Su_x(xi, eta)#

Update the row Su_x at a natural coordinate point

The entries of Su_x at the positions of the degrees of freedom \(u\), \(u_{,x}\), \(u_{,y}\) of the four nodes are computed from lex and ley; the product of Su_x with the element degrees of freedom gives \(u_{,x}\). The other entries are not changed (zero after construction).

Parameters:

xi, etadouble: Natural coordinates of the evaluation point. Dimensionless, from -1 to 1 over the element.

Notes

Su_x is \(\boldsymbol{S}^u_{,x} = (2/lex) \, \partial \boldsymbol{S}^u / \partial \xi\), with \(\boldsymbol{S}^u\) of update_Su().

update_Su_y(xi, eta)#

Update the row Su_y at a natural coordinate point

The entries of Su_y at the positions of the degrees of freedom \(u\), \(u_{,x}\), \(u_{,y}\) of the four nodes are computed from lex and ley; the product of Su_y with the element degrees of freedom gives \(u_{,y}\). The other entries are not changed (zero after construction).

Parameters:

xi, etadouble: Natural coordinates of the evaluation point. Dimensionless, from -1 to 1 over the element.

Notes

Su_y is \(\boldsymbol{S}^u_{,y} = (2/ley) \, \partial \boldsymbol{S}^u / \partial \eta\), with \(\boldsymbol{S}^u\) of update_Su().

update_Sv(xi, eta)#

Update the row Sv at a natural coordinate point

The entries of Sv at the positions of the degrees of freedom \(v\), \(v_{,x}\), \(v_{,y}\) of the four nodes are computed from lex and ley; the product of Sv with the element degrees of freedom gives \(v\). The other entries are not changed (zero after construction).

Parameters:

xi, etadouble: Natural coordinates of the evaluation point. Dimensionless, from -1 to 1 over the element.

Notes

Sv is the row \(\boldsymbol{S}^v\) of the interpolation of \(v\) of the Notes of bfscplate2d.bfscplate2d, with the functions \(H_i\), \(H_{x,i}\) and \(H_{y,i}\) of each node \(i\) at the positions of its degrees of freedom \(v\), \(v_{,x}\) and \(v_{,y}\).

update_Sv_x(xi, eta)#

Update the row Sv_x at a natural coordinate point

The entries of Sv_x at the positions of the degrees of freedom \(v\), \(v_{,x}\), \(v_{,y}\) of the four nodes are computed from lex and ley; the product of Sv_x with the element degrees of freedom gives \(v_{,x}\). The other entries are not changed (zero after construction).

Parameters:

xi, etadouble: Natural coordinates of the evaluation point. Dimensionless, from -1 to 1 over the element.

Notes

Sv_x is \(\boldsymbol{S}^v_{,x} = (2/lex) \, \partial \boldsymbol{S}^v / \partial \xi\), with \(\boldsymbol{S}^v\) of update_Sv().

update_Sv_y(xi, eta)#

Update the row Sv_y at a natural coordinate point

The entries of Sv_y at the positions of the degrees of freedom \(v\), \(v_{,x}\), \(v_{,y}\) of the four nodes are computed from lex and ley; the product of Sv_y with the element degrees of freedom gives \(v_{,y}\). The other entries are not changed (zero after construction).

Parameters:

xi, etadouble: Natural coordinates of the evaluation point. Dimensionless, from -1 to 1 over the element.

Notes

Sv_y is \(\boldsymbol{S}^v_{,y} = (2/ley) \, \partial \boldsymbol{S}^v / \partial \eta\), with \(\boldsymbol{S}^v\) of update_Sv().

update_Sw(xi, eta)#

Update the row Sw at a natural coordinate point

The entries of Sw at the positions of the degrees of freedom \(w\), \(w_{,x}\), \(w_{,y}\), \(w_{,xy}\) of the four nodes are computed from lex and ley; the product of Sw with the element degrees of freedom gives \(w\). The other entries are not changed (zero after construction).

Parameters:

xi, etadouble: Natural coordinates of the evaluation point. Dimensionless, from -1 to 1 over the element.

Notes

Sw is the row \(\boldsymbol{S}^w\) of the interpolation of \(w\) of the Notes of bfscplate2d.bfscplate2d, with the functions \(H_i\), \(H_{x,i}\), \(H_{y,i}\) and \(H_{xy,i}\) of each node \(i\) at the positions of its degrees of freedom \(w\), \(w_{,x}\), \(w_{,y}\) and \(w_{,xy}\), the interpolation of the BFS element.

update_Sw_x(xi, eta)#

Update the row Sw_x at a natural coordinate point

The entries of Sw_x at the positions of the degrees of freedom \(w\), \(w_{,x}\), \(w_{,y}\), \(w_{,xy}\) of the four nodes are computed from lex and ley; the product of Sw_x with the element degrees of freedom gives \(w_{,x}\). The other entries are not changed (zero after construction).

Parameters:

xi, etadouble: Natural coordinates of the evaluation point. Dimensionless, from -1 to 1 over the element.

Notes

Sw_x is \(\boldsymbol{S}^w_{,x} = (2/lex) \, \partial \boldsymbol{S}^w / \partial \xi\), with \(\boldsymbol{S}^w\) of update_Sw(), the first row of the operator \(\boldsymbol{G}\) of the Notes of bfscplate2d.bfscplate2d, which gives the slope \(w_{,x}\) of the nonlinear membrane strains.

update_Sw_y(xi, eta)#

Update the row Sw_y at a natural coordinate point

The entries of Sw_y at the positions of the degrees of freedom \(w\), \(w_{,x}\), \(w_{,y}\), \(w_{,xy}\) of the four nodes are computed from lex and ley; the product of Sw_y with the element degrees of freedom gives \(w_{,y}\). The other entries are not changed (zero after construction).

Parameters:

xi, etadouble: Natural coordinates of the evaluation point. Dimensionless, from -1 to 1 over the element.

Notes

Sw_y is \(\boldsymbol{S}^w_{,y} = (2/ley) \, \partial \boldsymbol{S}^w / \partial \eta\), with \(\boldsymbol{S}^w\) of update_Sw(), the second row of the operator \(\boldsymbol{G}\) of the Notes of bfscplate2d.bfscplate2d, which gives the slope \(w_{,y}\) of the nonlinear membrane strains.

bfscplate2d.bfscplate2d.update_KC0(plate, points, weights, KC0r, KC0c, KC0v)#

Update sparse vectors for linear constitutive stiffness matrix KC0

The element matrix, of size 40 by 40, is computed from lex, ley and the A, B and D blocks of ABD of plate, with nint by nint integration points, nint being the size of points. Its KC0_SPARSE_SIZE = 1600 entries, all the pairs of degrees of freedom of the element, are written from the position init_k_KC0 to init_k_KC0 + KC0_SPARSE_SIZE - 1: the global row and column indices to KC0r and KC0c, and the values added to KC0v. The attributes of plate are not changed.

Parameters:

plateBFSCPlate2D object: The plate element being added to KC0 Its attributes c1, c2, c3, c4, lex, ley, ABD and init_k_KC0 are read and must have been set.
pointsnp.array: Integration poins 1-D array of size nint with the points of the one-dimensional quadrature in natural coordinates, dimensionless, from -1 to 1, e.g. from bfscplate2d.quadrature.get_points_weights(). The same points are used along xi and eta.
weightsnp.array: Integration weights 1-D array with the weights of points, of size at least nint, dimensionless.
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. The entries at these positions are overwritten.
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 times length divided by the product of the dimensions of the row and column degrees of freedom, e.g. force per unit length between two of \(u\), \(v\), \(w\).

See also

update_KCNL, update_KG, update_M

Notes

The element matrix is the Hessian of the strain energy \(U\) of the Notes of bfscplate2d.bfscplate2d at \(\boldsymbol{u}_e = \boldsymbol{0}\), with the operators \(\boldsymbol{B}_m\) and \(\boldsymbol{B}_b\) of the linear membrane strains and of the curvatures:

\[\boldsymbol{K}_{C0} = \sum_{i=1}^{nint} \sum_{j=1}^{nint} \frac{lex \, ley}{4} \, weight \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)\]

with weight = weights[i]*weights[j] and the operators evaluated at xi = points[i] and eta = points[j]. The terms are the membrane stiffness, the membrane-bending coupling, null when \(\boldsymbol{B} = \boldsymbol{0}\), e.g. for a symmetric laminate with the reference surface at mid-thickness, and the bending stiffness. The matrix is symmetric and independent of the displacements. The integration is exact with nint >= 4, see the Notes of bfscplate2d.bfscplate2d.

bfscplate2d.bfscplate2d.update_KCNL(u, plate, points, weights, KCNLr, KCNLc, KCNLv)#

Update sparse vectors for nonlinear stiffness matrix KCNL

The 40 degrees of freedom of the element are read from u and written to the attribute u of plate, which is overwritten. The element matrix is computed from them, with the slopes \(w_{,x}\) and \(w_{,y}\) they give at each integration point, and from lex, ley and the A and B blocks of ABD (D is not read), with nint by nint integration points, nint being the size of points. Its KCNL_SPARSE_SIZE = 1024 entries are the pairs of an in-plane degree of freedom (\(u\) to \(v_{,y}\) of a node) with an out-of-plane one (\(w\) to \(w_{,xy}\)), in both orders, and the pairs of two out-of-plane degrees of freedom; the pairs of two in-plane degrees of freedom are not written. They are written from the position init_k_KCNL to init_k_KCNL + KCNL_SPARSE_SIZE - 1: the global row and column indices to KCNLr and KCNLc, and the values added to KCNLv.

Parameters:

unp.array: Global displacement vector 1-D array of DOUBLE with the degrees of freedom of all nodes, read at the positions c1 to c1 + 9, …, c4 to c4 + 9, with the dimensions given in the module documentation.
plateBFSCPlate2D object: The plate element being added to K Its attributes c1, c2, c3, c4, lex, ley, ABD and init_k_KCNL are read and must have been set.
pointsnp.array: Integration poins 1-D array of size nint with the points of the one-dimensional quadrature in natural coordinates, dimensionless, from -1 to 1, e.g. from bfscplate2d.quadrature.get_points_weights(). The same points are used along xi and eta.
weightsnp.array: Integration weights 1-D array with the weights of points, of size at least nint, dimensionless.
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. The entries at these positions are overwritten.
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. Dimension of force times length divided by the product of the dimensions of the row and column degrees of freedom, e.g. force per unit length between two of \(u\), \(v\), \(w\).

See also

update_KC0, update_KG, update_fint

Notes

KCNL = KC0L + KCL0 + KCLL + KGLL is the part of the tangent stiffness matrix \(\partial \boldsymbol{F}_{int} / \partial \boldsymbol{u}_e\) of the Notes of bfscplate2d.bfscplate2d that is neither in KC0 nor in the KG of update_KG(), 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} = \boldsymbol{K}_{GLL}\) of the membrane forces of the nonlinear strains, with

\[ \begin{align}\begin{aligned}\boldsymbol{K}_{C0L} = \sum_{i=1}^{nint} \sum_{j=1}^{nint} \frac{lex \, ley}{4} \, weight \left( \boldsymbol{B}_m^\top \boldsymbol{A} \boldsymbol{B}_{NL} + \boldsymbol{B}_b^\top \boldsymbol{B} \boldsymbol{B}_{NL} \right)\\\boldsymbol{K}_{CL0} = \sum_{i=1}^{nint} \sum_{j=1}^{nint} \frac{lex \, ley}{4} \, weight \left( \boldsymbol{B}_{NL}^\top \boldsymbol{A} \boldsymbol{B}_m + \boldsymbol{B}_{NL}^\top \boldsymbol{B} \boldsymbol{B}_b \right)\\\boldsymbol{K}_{CLL} = \sum_{i=1}^{nint} \sum_{j=1}^{nint} \frac{lex \, ley}{4} \, weight \boldsymbol{B}_{NL}^\top \boldsymbol{A} \boldsymbol{B}_{NL}\\\begin{split}\boldsymbol{K}_{GLL} = \sum_{i=1}^{nint} \sum_{j=1}^{nint} \frac{lex \, ley}{4} \, weight \boldsymbol{G}^\top \begin{bmatrix} N_{NL,xx} & N_{NL,xy} \\ N_{NL,xy} & N_{NL,yy} \end{bmatrix} \boldsymbol{G}\end{split}\end{aligned}\end{align} \]

with weight = weights[i]*weights[j] and the operators evaluated at xi = points[i] and eta = points[j], \(\boldsymbol{B}_{NL}\) computed with the slopes \(w_{,x}\) and \(w_{,y}\) (w_x and w_y) of the element degrees of freedom at the point, and the membrane forces of the nonlinear strains

\[\begin{split}\begin{Bmatrix} N_{NL,xx} \\ N_{NL,yy} \\ N_{NL,xy} \end{Bmatrix} = \boldsymbol{A} \begin{Bmatrix} w_{,x}^2/2 \\ w_{,y}^2/2 \\ w_{,x} w_{,y} \end{Bmatrix}\end{split}\]

KC0L and KCL0 are linear in the slopes, KCLL and KGLL quadratic, and \(\boldsymbol{K}_{CL0} = \boldsymbol{K}_{C0L}^\top\). KGLL is the geometric term of the membrane forces of the nonlinear strains, whereas the geometric term of the membrane forces of the linear strains is the KG of update_KG(), such that KC0 + KCNL + KG is the exact Jacobian of the internal forces of update_fint() computed with the same points. Since \(\boldsymbol{B}_{NL}\) and \(\boldsymbol{G}\) have entries only in the columns of the out-of-plane degrees of freedom, KCNL has no entries between two in-plane degrees of freedom. The integration is exact with nint >= 7, see the Notes of bfscplate2d.bfscplate2d.

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

Update sparse vectors for the load stiffness matrix of a follower pressure

The 40 degrees of freedom of the element are read from u and written to the attribute u of plate, which is overwritten. With quadratic not 0, the displacement gradients \(u_{,x}\), \(u_{,y}\), \(v_{,x}\), \(v_{,y}\), \(w_{,x}\) and \(w_{,y}\) are computed from them at each of the nint by nint integration points, nint being the size of points; with quadratic=0 the matrix does not depend on u. The element matrix, minus the derivative of the force vector of update_fext_follower() with respect to the element degrees of freedom, is computed with lex, ley and the pressure p. It is unsymmetric in general and its KCfollower_SPARSE_SIZE = 1056 entries, written in full (not as an upper triangle), are the pairs of a row and a column degree of freedom of different displacements among \(u\), \(v\) and \(w\): the blocks \([u, v]\), \([u, w]\), \([v, u]\), \([v, w]\), \([w, u]\) and \([w, v]\), those of \(u\) and \(v\) including their nodal derivatives and those of \(w\) its nodal derivatives \(w_{,x}\), \(w_{,y}\), \(w_{,xy}\). With quadratic=0 the values of the blocks \([u, v]\) and \([v, u]\) are zero. They are written from the position init_k_KCfollower to init_k_KCfollower + KCfollower_SPARSE_SIZE - 1: the global row and column indices to KCfollowerr and KCfollowerc, and the values added to KCfollowerv.

Parameters:

unp.array: Global displacement vector 1-D array of DOUBLE with the degrees of freedom of all nodes, read at the positions c1 to c1 + 9, …, c4 to c4 + 9, with the dimensions given in the module documentation. Only used with quadratic not 0.
plateBFSCPlate2D object: The plate element being added to KCfollower Its attributes c1, c2, c3, c4, lex, ley and init_k_KCfollower are read and must have been set.
pointsnp.array: Integration poins 1-D array of size nint with the points of the one-dimensional quadrature in natural coordinates, dimensionless, from -1 to 1, e.g. from bfscplate2d.quadrature.get_points_weights(). The same points are used along xi and eta.
weightsnp.array: Integration weights 1-D array with the weights of points, of size at least nint, dimensionless.
pnp.array: Pressure at the integration points 2-D C-contiguous array of DOUBLE of shape at least (nint, nint), p[i, j] being the pressure at xi = points[i] and eta = points[j], the index i going with xi and j with eta. Force per unit area of the undeformed mid-surface, a function of the undeformed coordinates, already multiplied by the load factor if there is one. A positive value pushes along \(+z\).
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. The entries at these positions are overwritten.
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. Dimension of force times length divided by the product of the dimensions of the row and column degrees of freedom, e.g. force per unit length between two of \(u\), \(v\), \(w\).
quadraticint, optional: 0 for the area vector truncated to first order in the displacement gradients, any other value for the complete area vector, as in update_fext_follower(). Default 0.

See also

update_fext_follower: force vector of the follower pressure.

Notes

The matrix is the load stiffness of the follower pressure [1] [2],

\[\boldsymbol{K}_F = -\frac{\partial \boldsymbol{F}_p}{\partial \boldsymbol{u}_e} = -\sum_{i=1}^{nint} \sum_{j=1}^{nint} \frac{lex \, ley}{4} \, weight \, p_{ij} \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)\]

with \(\boldsymbol{F}_p\) and the area vector \(\boldsymbol{a}\) of update_fext_follower(), weight = weights[i]*weights[j] and \(p_{ij}\) = p[i, j], the rows evaluated at xi = points[i] and eta = points[j]. For the first-order area vector

\[\begin{split}\frac{\partial \boldsymbol{a}_1}{\partial \boldsymbol{u}_e} = \begin{bmatrix} -\boldsymbol{S}^w_{,x} \\ -\boldsymbol{S}^w_{,y} \\ \boldsymbol{S}^u_{,x} + \boldsymbol{S}^v_{,y} \end{bmatrix}\end{split}\]

is constant, which gives the blocks \(p (\boldsymbol{S}^u)^\top \boldsymbol{S}^w_{,x}\), \(p (\boldsymbol{S}^v)^\top \boldsymbol{S}^w_{,y}\), \(-p (\boldsymbol{S}^w)^\top \boldsymbol{S}^u_{,x}\) and \(-p (\boldsymbol{S}^w)^\top \boldsymbol{S}^v_{,y}\) integrated over the element, independent of the displacements. The complete area vector adds \(\partial \boldsymbol{a}_2 / \partial \boldsymbol{u}_e\), linear in the displacement gradients, to these blocks and to the blocks \([u, v]\) and \([v, u]\). In the tangent stiffness matrix of a nonlinear analysis with the load factor \(\lambda\) of the pressure, \(\boldsymbol{K}_{C0} + \boldsymbol{K}_{CNL} + \boldsymbol{K}_{GNL} + \boldsymbol{K}_G + \boldsymbol{K}_F\), with \(\boldsymbol{K}_F = \lambda \boldsymbol{K}_{F_{ref}}\) and \(\boldsymbol{K}_{F_{ref}}\) the matrix of p, is the Jacobian of the residual \(\boldsymbol{R} = \boldsymbol{F}_{int} - \lambda \boldsymbol{F}_p(\boldsymbol{u})\), see update_fext_follower(). In a linear buckling analysis of a follower pressure \(\lambda p\), \(\boldsymbol{K}_{F_{ref}}\) of p multiplies \(\lambda\) together with the \(\boldsymbol{K}_G\) of the linear static solution of p.

For the first-order area vector \(\boldsymbol{\delta}^\top \boldsymbol{K}_F \boldsymbol{\Delta} = -P(\delta, \Delta)\), with the bilinear form \(P(\delta, \Delta) = \int_A p \, \delta\boldsymbol{u} \cdot \boldsymbol{a}_1(\Delta\boldsymbol{u}) \, dA\) of two displacement fields \(\delta\) and \(\Delta\), and an integration by parts gives its skew part

\[P(\delta, \Delta) - P(\Delta, \delta) = - \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\}\) the outward normal of the boundary of the loaded region. 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 boundary of the loaded region, either \(w = 0\) or the in-plane displacement normal to the boundary vanishes, \(u = 0\) on the edges at constant \(x\) and \(v = 0\) on those at constant \(y\), such that the contour integral vanishes.

These conditions are met by simply supported and clamped edges and by symmetry planes, and the load potential is then 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 plate, and for a pressure that varies over the plate, the result of [1] and [2], 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 boundary 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 boundary of the loaded region of the plate. An unsymmetric matrix, a non-conservative load, must not be symmetrized and requires solvers of unsymmetric systems. The conditions are those of the first-order area vector; for quadratic=1 the matrix depends on \(u\). The integration is exact for the same nint as in update_fext_follower().

References

[1] (1,2)

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.

[2] (1,2)

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.

bfscplate2d.bfscplate2d.update_KG(u, plate, points, weights, KGr, KGc, KGv)#

Update sparse vectors for geometric stiffness matrix KG

The 40 degrees of freedom of the element are read from u and written to the attribute u of plate, which is overwritten. At each of the nint by nint integration points, nint being the size of points, the membrane forces Nxx, Nyy and Nxy are computed from them, linear in the degrees of freedom, with the A and B blocks of ABD (D is not read), and used with lex and ley to compute the element matrix. Its KG_SPARSE_SIZE = 256 entries are the pairs of out-of-plane degrees of freedom (\(w\), \(w_{,x}\), \(w_{,y}\), \(w_{,xy}\) of the four nodes). They are written from the position init_k_KG to init_k_KG + KG_SPARSE_SIZE - 1: the global row and column indices to KGr and KGc, and the values added to KGv.

Parameters:

unp.array: Global displacement vector 1-D array of DOUBLE with the degrees of freedom of all nodes, read at the positions c1 to c1 + 9, …, c4 to c4 + 9, with the dimensions given in the module documentation.
plateBFSCPlate2D object: The plate element being added to KG Its attributes c1, c2, c3, c4, lex, ley, ABD and init_k_KG are read and must have been set.
pointsnp.array: Integration poins 1-D array of size nint with the points of the one-dimensional quadrature in natural coordinates, dimensionless, from -1 to 1, e.g. from bfscplate2d.quadrature.get_points_weights(). The same points are used along xi and eta.
weightsnp.array: Integration weights 1-D array with the weights of points, of size at least nint, dimensionless.
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. The entries at these positions are overwritten.
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. Dimension of force times length divided by the product of the dimensions of the row and column degrees of freedom, e.g. force per unit length between two of \(u\), \(v\), \(w\).

See also

update_KG_cte_N: same matrix with constant membrane forces.

Notes

\[\begin{split}\boldsymbol{K}_G = \sum_{i=1}^{nint} \sum_{j=1}^{nint} \frac{lex \, ley}{4} \, weight \boldsymbol{G}^\top \begin{bmatrix} Nxx & Nxy \\ Nxy & Nyy \end{bmatrix} \boldsymbol{G}\end{split}\]

with weight = weights[i]*weights[j] and the operators evaluated at xi = points[i] and eta = points[j], \(\boldsymbol{G}\) of the Notes of bfscplate2d.bfscplate2d, and the membrane forces of the linear membrane strains and of the curvatures

\[\{Nxx, Nyy, Nxy\}^\top = \boldsymbol{A}\boldsymbol{B}_m\boldsymbol{u}_e + \boldsymbol{B}\boldsymbol{B}_b\boldsymbol{u}_e\]

It is the geometric term \(\partial (\boldsymbol{B}_{NL}^\top \boldsymbol{N}) / \partial \boldsymbol{u}_e = \boldsymbol{G}^\top \begin{bmatrix} Nxx & Nxy \\ Nxy & Nyy \end{bmatrix} \boldsymbol{G}\) of the tangent stiffness matrix for the part of the membrane forces \(\boldsymbol{N}\) that is linear in \(\boldsymbol{u}_e\); the part of the nonlinear strains is in the KGLL of update_KCNL(). Since KG is linear in the displacements, KG computed with the solution \(\boldsymbol{u}\) of the linear static problem \(\boldsymbol{K}_{C0}\boldsymbol{u} = \boldsymbol{F}_{ext}\) is proportional to the load, and the linear buckling problem

\[(\boldsymbol{K}_{C0} + \lambda \boldsymbol{K}_G) \boldsymbol{\phi} = \boldsymbol{0}\]

gives the load multipliers \(\lambda\) of the critical loads \(\lambda \boldsymbol{F}_{ext}\) and the buckling modes \(\boldsymbol{\phi}\). The integration is exact with nint >= 5, see the Notes of bfscplate2d.bfscplate2d.

bfscplate2d.bfscplate2d.update_KG_cte_N(Nxx, Nyy, Nxy, plate, points, weights, KGr, KGc, KGv)#

Update sparse vectors for geometric stiffness matrix KG

The element matrix is computed from the membrane forces Nxx, Nyy and Nxy, the same at all the nint by nint integration points, nint being the size of points, and from lex and ley; ABD and the displacements are not read. Its KG_SPARSE_SIZE = 256 entries, at the same pairs of degrees of freedom as in update_KG(), are written from the position init_k_KG to init_k_KG + KG_SPARSE_SIZE - 1: the global row and column indices to KGr and KGc, and the values added to KGv. The attributes of plate are not changed.

Parameters:

Nxx, Nyy, Nxyfloat: Constant membrane forces within element of the plate parameter Dimension of force per unit length.
plateBFSCPlate2D object: The plate element being added to KG Its attributes c1, c2, c3, c4, lex, ley and init_k_KG are read and must have been set.
pointsnp.array: Integration poins 1-D array of size nint with the points of the one-dimensional quadrature in natural coordinates, dimensionless, from -1 to 1, e.g. from bfscplate2d.quadrature.get_points_weights(). The same points are used along xi and eta.
weightsnp.array: Integration weights 1-D array with the weights of points, of size at least nint, dimensionless.
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. The entries at these positions are overwritten.
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. Dimension of force times length divided by the product of the dimensions of the row and column degrees of freedom, e.g. force per unit length between two of \(u\), \(v\), \(w\).

See also

update_KG: same matrix with the membrane forces from the displacements.

Notes

The matrix is the KG of update_KG() for prescribed membrane forces, uniform over the element:

with weight = weights[i]*weights[j] and the operators evaluated at xi = points[i] and eta = points[j] and \(\boldsymbol{G}\) of the Notes of bfscplate2d.bfscplate2d. The membrane forces are positive in tension, e.g. Nxx = -1 in all the elements for a unit compressive force per unit length along \(x\), for which the load multiplier \(\lambda\) of the linear buckling problem \((\boldsymbol{K}_{C0} + \lambda \boldsymbol{K}_G) \boldsymbol{\phi} = \boldsymbol{0}\) is the critical compressive force per unit length. The membrane forces of a pre-buckling state computed with the stiffness of the structure, e.g. non-uniform because of the boundary conditions, are given instead by update_KG(). The integration is exact with nint >= 4, see the Notes of bfscplate2d.bfscplate2d.

bfscplate2d.bfscplate2d.update_M(plate, Mr, Mc, Mv)#

Update sparse vectors for mass matrix M

The element matrix is computed with explicit expressions, without numerical integration, from lex, ley, h and rho of plate; ABD is not read. Its M_SPARSE_SIZE = 544 entries are, for each pair of nodes, the pairs of degrees of freedom among \(u\), \(u_{,x}\), \(u_{,y}\), the pairs among \(v\), \(v_{,x}\), \(v_{,y}\) and the pairs among \(w\), \(w_{,x}\), \(w_{,y}\), \(w_{,xy}\). They are written from the position init_k_M to init_k_M + M_SPARSE_SIZE - 1: the global row and column indices to Mr and Mc, and the values added to Mv. The attributes of plate are not changed.

Parameters:

plateBFSCPlate2D object: The plate element being added to M Its attributes c1, c2, c3, c4, lex, ley, h, rho and init_k_M are read and must have been set.
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. The entries at these positions are overwritten.
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 times length squared divided by the product of the dimensions of the row and column degrees of freedom, e.g. mass between two of \(u\), \(v\), \(w\).

See also

update_KC0

Notes

The element matrix is the consistent mass matrix of the kinetic energy of the kinematics of the Notes of bfscplate2d.bfscplate2d, with the displacements \(u - z w_{,x}\), \(v - z w_{,y}\) and \(w\) at the distance \(z\) from the reference surface, for the density rho constant through the thickness h and the reference surface at mid-thickness:

\[\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 term with \(h^3/12\) is the rotary inertia of the normal. The terms with \(\int \rho \, z \, dz\), which couple \(u\) and \(v\) with the slopes of \(w\), vanish with the reference surface at mid-thickness, such that M has no entries between in-plane and out-of-plane degrees of freedom; they are not included for a reference surface elsewhere, e.g. an offset one. The integrands are polynomials, integrated in closed form.

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

Update global force vector of a follower pressure

The 40 degrees of freedom of the element are read from u and written to the attribute u of plate, which is overwritten. At each of the nint by nint integration points, nint being the size of points, the displacement gradients \(u_{,x}\), \(u_{,y}\), \(v_{,x}\), \(v_{,y}\), \(w_{,x}\) and \(w_{,y}\) are computed from them with lex and ley, and from these the area vector of the deformed mid-surface, along which the pressure p acts. The element force vector of the pressure is added to fext at the positions c1 to c1 + 9, …, c4 to c4 + 9. ABD and the init_k_* attributes are not used.

Parameters:

unp.array: Global displacement vector 1-D array of DOUBLE with the degrees of freedom of all nodes, read at the positions c1 to c1 + 9, …, c4 to c4 + 9, with the dimensions given in the module documentation.
plateBFSCPlate2D object: The plate element being added to fext Its attributes c1, c2, c3, c4, lex and ley are read and must have been set.
pointsnp.array: Integration poins 1-D array of size nint with the points of the one-dimensional quadrature in natural coordinates, dimensionless, from -1 to 1, e.g. from bfscplate2d.quadrature.get_points_weights(). The same points are used along xi and eta.
weightsnp.array: Integration weights 1-D array with the weights of points, of size at least nint, dimensionless.
pnp.array: Pressure at the integration points 2-D C-contiguous array of DOUBLE of shape at least (nint, nint), p[i, j] being the pressure at xi = points[i] and eta = points[j], the index i going with xi and j with eta. Force per unit area of the undeformed mid-surface, a function of the undeformed coordinates, already multiplied by the load factor if there is one. A positive value pushes along \(+z\).
fextnp.array: Global external force vector 1-D array of DOUBLE with the same positions as u. The values are added to the existing ones. Dimension of force times length divided by the dimension of the degree of freedom, e.g. force for \(u\), \(v\), \(w\).
quadraticint, optional: 0 for the area vector truncated to first order in the displacement gradients, any other value for the complete area vector. Default 0.
referenceint, optional: Any value other than 0 includes the force of the undeformed state, the dead pressure load; 0 leaves it out, which gives the part \(\boldsymbol{F}_p(\boldsymbol{u}) - \boldsymbol{F}_p(\boldsymbol{0})\) that depends on the displacements. Default 1.

See also

update_KCfollower: load stiffness matrix, minus the derivative of this force vector with respect to the degrees of freedom.

Notes

The pressure \(p\), per unit area of the undeformed mid-surface, is normal to the deformed mid-surface and acts on its deformed area [1] [2]. Its virtual work is

\[\delta W_p = \int_A p \, \boldsymbol{a} \cdot \delta\boldsymbol{u} \, dx \, dy, \qquad \boldsymbol{a} = \boldsymbol{X}_{,x} \times \boldsymbol{X}_{,y}\]

with \(\boldsymbol{u} = \{u, v, w\}^\top\), \(\boldsymbol{X} = \boldsymbol{X}_0 + \boldsymbol{u}\) the deformed position of the mid-surface, \(\boldsymbol{X}_{,x} = \{1 + u_{,x}, v_{,x}, w_{,x}\}^\top\) and \(\boldsymbol{X}_{,y} = \{u_{,y}, 1 + v_{,y}, w_{,y}\}^\top\). \(\boldsymbol{a}\) is the area vector, normal to the deformed mid-surface, with the length of the ratio between the deformed and the undeformed areas. Split by orders in the displacement gradients, \(\boldsymbol{a} = \boldsymbol{a}_0 + \boldsymbol{a}_1 + \boldsymbol{a}_2\), with

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

The function adds

\[\boldsymbol{F}_p = \sum_{i=1}^{nint} \sum_{j=1}^{nint} \frac{lex \, ley}{4} \, weight \, p_{ij} \left( (\boldsymbol{S}^u)^\top a_x + (\boldsymbol{S}^v)^\top a_y + (\boldsymbol{S}^w)^\top a_z \right), \qquad \boldsymbol{a} = r \, \boldsymbol{a}_0 + \boldsymbol{a}_1 + q \, \boldsymbol{a}_2\]

with weight = weights[i]*weights[j], \(p_{ij}\) = p[i, j], the rows \(\boldsymbol{S}^u\), \(\boldsymbol{S}^v\), \(\boldsymbol{S}^w\) of the Notes of bfscplate2d.bfscplate2d and the area vector evaluated at xi = points[i] and eta = points[j], \(r = 1\) for reference not 0 and \(q = 1\) for quadratic not 0, otherwise \(0\). The term of \(\boldsymbol{a}_0\) is the dead pressure load \(\int_A p (\boldsymbol{S}^w)^\top dA\), the force of the undeformed state.

With the moderate rotations of the von Karman strains, \(u_{,x} \sim \epsilon\) and \(w_{,x} \sim \epsilon^{1/2}\), such that \(\delta u \sim \epsilon^{1/2} \delta w\), \(\boldsymbol{a}_0 + \boldsymbol{a}_1\) keeps every term of \(p \, \boldsymbol{a} \cdot \delta\boldsymbol{u}\) up to the order \(\epsilon \, p \, \delta w\), and \(\boldsymbol{a}_2\) adds terms of the order \(\epsilon^2 p \, \delta w\), of the order of the terms that the strains neglect. quadratic=0 is therefore the consistent truncation, for which \(\boldsymbol{F}_p\) is affine in the displacements, and quadratic=1 the complete area vector of the plate, for checking the sensitivity of a result to the truncation.

A positive \(p\) pushes along \(+z\) in the undeformed state, like a positive load along \(w\). In a nonlinear analysis with the load factor \(\lambda\) of the pressure, the residual of the equilibrium equations is \(\boldsymbol{R} = \boldsymbol{F}_{int} - \lambda \boldsymbol{F}_p(\boldsymbol{u})\), where \(\boldsymbol{F}_p(\boldsymbol{u}) = \boldsymbol{F}_p(\boldsymbol{0}) + (\boldsymbol{F}_p(\boldsymbol{u}) - \boldsymbol{F}_p(\boldsymbol{0}))\): the dead part may stay in the external force vector of the solver and the part of reference=0, which vanishes at \(\boldsymbol{u} = \boldsymbol{0}\), be subtracted from the internal forces. The Jacobian of the residual is \(\boldsymbol{K}_{C0} + \boldsymbol{K}_{CNL} + \boldsymbol{K}_{GNL} + \boldsymbol{K}_G + \boldsymbol{K}_F\), with \(\boldsymbol{K}_F = \lambda \boldsymbol{K}_{F_{ref}}\) and \(\boldsymbol{K}_{F_{ref}}\) the matrix of update_KCfollower() for the same p and quadratic.

For a pressure constant over the element, the integrands are polynomials of degree at most 6 in \(\xi\) and in \(\eta\) for quadratic=0 and 8 for quadratic=1, integrated exactly with the Gauss-Legendre points of nint >= 4 and nint >= 5.

References

[1]

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.

[2]

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.

bfscplate2d.bfscplate2d.update_fint(u, plate, points, weights, fint)#

Update global internal force vector

The 40 degrees of freedom of the element are read from u and written to the attribute u of plate, which is overwritten. At each of the nint by nint integration points, nint being the size of points, the membrane forces and the moments are computed from them with the A, B and D blocks of ABD, including the terms that depend on the slopes \(w_{,x}\) and \(w_{,y}\) (there is no option for the linear internal forces only). The element internal forces computed from them are added to fint at the positions c1 to c1 + 9, …, c4 to c4 + 9. The init_k_* attributes are not used.

Parameters:

unp.array: Global displacement vector 1-D array of DOUBLE with the degrees of freedom of all nodes, read at the positions c1 to c1 + 9, …, c4 to c4 + 9, with the dimensions given in the module documentation.
plateBFSCPlate2D object: The plate element being added to fint Its attributes c1, c2, c3, c4, lex, ley and ABD are read and must have been set.
pointsnp.array: Integration poins 1-D array of size nint with the points of the one-dimensional quadrature in natural coordinates, dimensionless, from -1 to 1, e.g. from bfscplate2d.quadrature.get_points_weights(). The same points are used along xi and eta.
weightsnp.array: Integration weights 1-D array with the weights of points, of size at least nint, dimensionless.
fintnp.array: Global internal force array 1-D array of DOUBLE with the same positions as u. The values are added to the existing ones. Dimension of force times length divided by the dimension of the degree of freedom, e.g. force for \(u\), \(v\), \(w\).

See also

update_KC0, update_KCNL, update_KG

Notes

The internal forces are the gradient of the strain energy \(U\) of the Notes of bfscplate2d.bfscplate2d with respect to the element degrees of freedom:

\[\boldsymbol{F}_{int} = \sum_{i=1}^{nint} \sum_{j=1}^{nint} \frac{lex \, ley}{4} \, weight \left( (\boldsymbol{B}_m + \boldsymbol{B}_{NL})^\top \boldsymbol{N} + \boldsymbol{B}_b^\top \boldsymbol{M} \right)\]

with weight = weights[i]*weights[j] and the operators evaluated at xi = points[i] and eta = points[j], \(\boldsymbol{B}_{NL}\) computed with the slopes w_x and w_y at the point, and the membrane forces and moments of the total strains

\[ \begin{align}\begin{aligned}\{Nxx, Nyy, Nxy\}^\top = \boldsymbol{A} \left( \boldsymbol{B}_m + \frac{1}{2} \boldsymbol{B}_{NL} \right) \boldsymbol{u}_e + \boldsymbol{B}\boldsymbol{B}_b\boldsymbol{u}_e\\\{Mxx, Myy, Mxy\}^\top = \boldsymbol{B} \left( \boldsymbol{B}_m + \frac{1}{2} \boldsymbol{B}_{NL} \right) \boldsymbol{u}_e + \boldsymbol{D}\boldsymbol{B}_b\boldsymbol{u}_e\end{aligned}\end{align} \]

The factor \(1/2\) of the strains and the full \(\boldsymbol{B}_{NL}\) of their variation come from the nonlinear strains being quadratic in \(\boldsymbol{u}_e\), see the Notes of bfscplate2d.bfscplate2d. KC0 + KCNL + KG of update_KC0(), update_KCNL() and update_KG() is the exact Jacobian of these internal forces computed with the same points, and for small displacements \(\boldsymbol{F}_{int} \approx \boldsymbol{K}_{C0}\boldsymbol{u}_e\). In a nonlinear static analysis the residual is the difference between the assembled internal force vector and the external force vector, \(\boldsymbol{R} = \boldsymbol{F}_{int} - \boldsymbol{F}_{ext}\). The integration is exact with nint >= 7, see the Notes of bfscplate2d.bfscplate2d.

BFSC plate element (bfscplate2d.bfscplate2d)#

Notes#

References#

BFSC plate element (`bfscplate2d.bfscplate2d`)#