Tutorials#
The tutorials are Jupyter notebooks, stored in doc/source/tutorials/ with
their outputs, which compare models of stringer-stiffened cylinders with
published tests: the skin is modelled with the elements of bfsccylinder and
each blade stringer with the BFSC plate element of bfscplate2d. Both notebooks build their
models with the module stiffened_cylinder.py of the same folder, below.
Running them requires bfscplate2d 0.5.0 or newer, with the shape functions
named Su, …, Sw_y, composites, structsolve and Matplotlib, besides
bfsccylinder.
Model of the stiffened cylinders#
The skin and the stringers are assembled as in Assembling a model, with the degrees of freedom of the stringer nodes along the skin eliminated by the constraint basis \(\boldsymbol{Z}\) described in the module docstring:
r"""Stringer-stiffened cylinder of BFSC elements
The skin is modelled with the BFSC cylinder elements of ``bfsccylinder``,
Sanders-type or Donnell-type, and each blade stringer with BFSC plate
elements of ``bfscplate2d``, in the plane `(x, r)` of the stringer, from the
middle surface of the skin to the tip of the stringer. The plate of a
stringer has the local coordinates `X = x` and `Y = r - R`, positive
outwards, and its displacements `U`, `V` and `W` are along `x`, along the
outward radius and along the circumference, the direction of `v` of the
skin.
Along the line where a stringer meets the skin, `Y = 0`, the nodes of the
stringer are tied to the nodes of the skin, the translations and the
rotation about `x` being the same:
.. math::
U = u, \quad V = w, \quad W = v, \quad W_{,Y} = \beta_y = -w_{,y} +
\frac{v}{R}
and the derivatives of these relations along `x`: `U_{,X} = u_{,x}`,
`V_{,X} = w_{,x}`, `W_{,X} = v_{,x}`, `W_{,XY} = -w_{,xy} + v_{,x}/R`. The
rotation `\beta_y` is that of the Sanders-type kinematics, with which
rigid-body motions of the stiffened cylinder give no strains; it is also
used with the Donnell-type skin, whose kinematics neglect `v/R`. The
derivatives `U_{,Y}` and `V_{,Y}` of the stringer, its in-plane shear and
radial strains at the root, stay free. The 8 tied degrees of freedom of
each root node are eliminated with the constraint basis
`\boldsymbol{u} = \boldsymbol{Z}\boldsymbol{u}_r`, and the matrices of the reduced degrees of freedom
`\boldsymbol{u}_r` are `\boldsymbol{Z}^\top \boldsymbol{K} \boldsymbol{Z}`, `\boldsymbol{Z}^\top \boldsymbol{K}_G \boldsymbol{Z}` and `\boldsymbol{Z}^\top \boldsymbol{M} \boldsymbol{Z}`.
The model is the half of the cylinder `0 \le \theta \le \pi` with two
planes of symmetry midway between two stringers, or the whole cylinder.
"""
import numpy as np
from numpy import pi
from scipy.sparse import coo_matrix
from composites import isotropic_plate
import bfsccylinder
import bfsccylinder.sanders
from bfsccylinder.quadrature import get_points_weights
from bfsccylinder.utils import assign_constant_ABD
import bfscplate2d
from bfscplate2d import BFSCPlate2D
DOF = 10
# positions of the degrees of freedom in a node
U, U_X, U_Y, V, V_X, V_Y, W, W_X, W_Y, W_XY = range(10)
class StiffenedCylinder:
r"""Skin of BFSC cylinder elements with blade stringers of BFSC plate
elements
Parameters
----------
R, L, h : float
Radius of the middle surface, length and thickness of the skin.
E, nu, rho : float
Young's modulus, Poisson's ratio and density of the skin and of the
stringers.
nstringers : int
Number of stringers, equally spaced, even when ``half=True``.
bs, Hs : float
Thickness of a stringer and height of its plate, from the middle
surface of the skin to the tip.
side : str
``'external'`` or ``'internal'``.
nx : int
Number of nodes along `x`.
nbay : int
Number of skin elements between two stringers, even.
nh : int
Number of stringer elements along the height.
kinematics : str
``'sanders'`` or ``'donnell'``, of the skin elements.
half : bool
Half model with two planes of symmetry, or the whole cylinder.
free_axial : bool
With ``True``, the mean of `u` along the edge `x = 0` is tied to
zero, which removes the rigid-body axial translation of the boundary
conditions with `u` free at both ends.
"""
def __init__(self, R, L, h, E, nu, rho, nstringers, bs, Hs, side, nx,
nbay=2, nh=2, kinematics='sanders', half=True,
free_axial=False, nint=4):
assert side in ('external', 'internal')
assert kinematics in ('sanders', 'donnell')
assert nbay % 2 == 0
if half:
assert nstringers % 2 == 0
self.R, self.L, self.h = R, L, h
self.E, self.nu, self.rho = E, nu, rho
self.nstringers, self.bs, self.Hs, self.side = nstringers, bs, Hs, side
self.nx, self.nbay, self.nh = nx, nbay, nh
self.kinematics, self.half = kinematics, half
if kinematics == 'sanders':
self.mod = bfsccylinder.sanders
Element = bfsccylinder.sanders.BFSCCylinderSanders
else:
self.mod = bfsccylinder
Element = bfsccylinder.BFSCCylinder
self.points, self.weights = get_points_weights(nint=nint)
# skin nodes, x varying slowest
dtheta = 2*pi/nstringers/nbay
ny_el = nstringers*nbay//2 if half else nstringers*nbay
nyn = ny_el + 1 if half else ny_el
self.nyn = nyn
self.theta = np.arange(nyn)*dtheta
self.xlin = np.linspace(0, L, nx)
skin = np.arange(nx*nyn).reshape(nx, nyn)
self.skin_nodes = skin
num_skin_nodes = nx*nyn
self.x = np.repeat(self.xlin, nyn)
self.y = np.tile(self.theta, nx)*R
# stringers: positions midway between the planes of symmetry
nstr_model = nstringers//2 if half else nstringers
self.j_stringer = np.arange(nstr_model)*nbay + nbay//2
self.theta_stringer = self.theta[self.j_stringer]
sign = 1. if side == 'external' else -1.
self.Y = sign*np.arange(nh + 1)*Hs/nh
str_nodes = (num_skin_nodes
+ np.arange(nstr_model*nx*(nh + 1)).reshape(
nstr_model, nx, nh + 1))
self.str_nodes = str_nodes
num_nodes = num_skin_nodes + str_nodes.size
self.N = N = DOF*num_nodes
# elements
lex = L/(nx - 1)
self.skin_elements = []
for i in range(nx - 1):
for j in range(ny_el):
jp = (j + 1) % nyn if not half else j + 1
shell = Element(nint)
shell.n1 = skin[i, j]
shell.n2 = skin[i + 1, j]
shell.n3 = skin[i + 1, jp]
shell.n4 = skin[i, jp]
shell.c1 = DOF*shell.n1
shell.c2 = DOF*shell.n2
shell.c3 = DOF*shell.n3
shell.c4 = DOF*shell.n4
shell.R = R
shell.lex = lex
shell.ley = R*dtheta
shell.h = h
shell.rho = rho
self.skin_elements.append(shell)
self.skin_prop = isotropic_plate(thickness=h, E=E, nu=nu, rho=rho)
for shell in self.skin_elements:
assign_constant_ABD(shell, self.skin_prop)
self.str_prop = isotropic_plate(thickness=bs, E=E, nu=nu, rho=rho)
self.str_elements = []
for s in range(nstr_model):
for i in range(nx - 1):
for k in range(nh):
# node 1 at the smallest X and Y
k1, k4 = (k, k + 1) if sign > 0 else (k + 1, k)
plate = BFSCPlate2D()
plate.n1 = str_nodes[s, i, k1]
plate.n2 = str_nodes[s, i + 1, k1]
plate.n3 = str_nodes[s, i + 1, k4]
plate.n4 = str_nodes[s, i, k4]
plate.c1 = DOF*plate.n1
plate.c2 = DOF*plate.n2
plate.c3 = DOF*plate.n3
plate.c4 = DOF*plate.n4
plate.lex = lex
plate.ley = Hs/nh
plate.h = bs
plate.rho = rho
plate.ABD = self.str_prop.ABD
self.str_elements.append(plate)
# constraint basis, u = Z u_r
dependent = {}
for s in range(nstr_model):
for i in range(nx):
c = DOF*skin[i, self.j_stringer[s]]
b = DOF*str_nodes[s, i, 0]
dependent[b + U] = [(c + U, 1.)]
dependent[b + U_X] = [(c + U_X, 1.)]
dependent[b + V] = [(c + W, 1.)]
dependent[b + V_X] = [(c + W_X, 1.)]
dependent[b + W] = [(c + V, 1.)]
dependent[b + W_X] = [(c + V_X, 1.)]
dependent[b + W_Y] = [(c + W_Y, -1.), (c + V, 1/R)]
dependent[b + W_XY] = [(c + W_XY, -1.), (c + V_X, 1/R)]
if free_axial:
# trapezoidal weights of the edge x = 0
wts = np.ones(nyn)
if half:
wts[0] = wts[-1] = 0.5
c0 = DOF*skin[0, 0]
dependent[c0 + U] = [(DOF*skin[0, j] + U, -wts[j]/wts[0])
for j in range(1, nyn)]
independent = np.ones(N, dtype=bool)
independent[list(dependent.keys())] = False
self.Nr = Nr = independent.sum()
red = -np.ones(N, dtype=int)
red[independent] = np.arange(Nr)
self.red = red
rows = list(np.flatnonzero(independent))
cols = list(red[independent])
vals = [1.]*Nr
for d, terms in dependent.items():
for d2, coef in terms:
assert red[d2] >= 0
rows.append(d)
cols.append(red[d2])
vals.append(coef)
self.Z = coo_matrix((vals, (rows, cols)), shape=(N, Nr)).tocsr()
# assembly of the matrices of all the degrees of freedom
def _assemble(self, skin_func, skin_size, str_func, str_size):
nse = len(self.skin_elements)
nst = len(self.str_elements)
r = np.zeros(skin_size*nse + str_size*nst, dtype=self.mod.INT)
c = np.zeros_like(r)
v = np.zeros(r.shape[0], dtype=float)
init_k = 0
for shell in self.skin_elements:
skin_func(shell, init_k, r, c, v)
init_k += skin_size
for plate in self.str_elements:
str_func(plate, init_k, r, c, v)
init_k += str_size
return coo_matrix((v, (r, c)), shape=(self.N, self.N)).tocsr()
def reduce(self, A):
r"""Matrix of the reduced degrees of freedom, `\boldsymbol{Z}^\top \boldsymbol{A} \boldsymbol{Z}`"""
return (self.Z.T @ A @ self.Z).tocsr()
def KC0(self):
"""Linear stiffness matrix, all the degrees of freedom"""
m, p, w = self.mod, self.points, self.weights
def skin(shell, k, r, c, v):
shell.init_k_KC0 = k
m.update_KC0(shell, p, w, r, c, v)
def stringer(plate, k, r, c, v):
plate.init_k_KC0 = k
bfscplate2d.update_KC0(plate, p, w, r, c, v)
return self._assemble(skin, m.KC0_SPARSE_SIZE, stringer,
bfscplate2d.KC0_SPARSE_SIZE)
def M(self):
"""Mass matrix, all the degrees of freedom"""
m = self.mod
def skin(shell, k, r, c, v):
shell.init_k_M = k
m.update_M(shell, r, c, v)
def stringer(plate, k, r, c, v):
plate.init_k_M = k
bfscplate2d.update_M(plate, r, c, v)
return self._assemble(skin, m.M_SPARSE_SIZE, stringer,
bfscplate2d.M_SPARSE_SIZE)
def KG(self, u):
"""Geometric stiffness matrix of the displacements ``u`` of all the
degrees of freedom, with the in-plane terms of the stringers of
:meth:`KG_inplane`"""
m, p, w = self.mod, self.points, self.weights
def skin(shell, k, r, c, v):
shell.init_k_KG = k
m.update_KG(u, shell, p, w, r, c, v)
def stringer(plate, k, r, c, v):
plate.init_k_KG = k
bfscplate2d.update_KG(u, plate, p, w, r, c, v)
KG = self._assemble(skin, m.KG_SPARSE_SIZE, stringer,
bfscplate2d.KG_SPARSE_SIZE)
return KG + self.KG_inplane(self.stringer_membrane_forces(u))
def KG_membrane(self, eps):
r"""Geometric stiffness matrix of the uniform axial strain ``eps``
of the membrane pre-buckling state, `N_{xx} = E h \epsilon` in the
skin and `E b_s \epsilon` in the stringers, without hoop stress,
with the in-plane terms of the stringers of :meth:`KG_inplane`"""
m, p, w = self.mod, self.points, self.weights
Nskin = self.E*self.h*eps
Nstr = self.E*self.bs*eps
def skin(shell, k, r, c, v):
shell.init_k_KG = k
m.update_KG_constant_stress(shell, p, w, r, c, v, Nskin, 0., 0.)
def stringer(plate, k, r, c, v):
plate.init_k_KG = k
bfscplate2d.update_KG_cte_N(Nstr, 0., 0., plate, p, w, r, c, v)
KG = self._assemble(skin, m.KG_SPARSE_SIZE, stringer,
bfscplate2d.KG_SPARSE_SIZE)
nel, npts = len(self.str_elements), len(p)**2
N = np.zeros((nel, npts, 3))
N[:, :, 0] = Nstr
return KG + self.KG_inplane(N)
def _stringer_operators(self):
"""Interpolation rows of the stringer elements at the integration
points, the same for all the stringer elements"""
if hasattr(self, '_ops'):
return self._ops
plate = self.str_elements[0]
p, w = self.points, self.weights
Bm, Sv_x, Sv_y, wts = [], [], [], []
for i, xi in enumerate(p):
for j, eta in enumerate(p):
plate.update_Bm(xi, eta)
plate.update_Sv_x(xi, eta)
plate.update_Sv_y(xi, eta)
Bm.append(np.array(plate.Bm))
Sv_x.append(np.array(plate.Sv_x))
Sv_y.append(np.array(plate.Sv_y))
wts.append(w[i]*w[j]*plate.lex*plate.ley/4)
dofs = np.array([[pl.c1 + d for d in range(DOF)]
+ [pl.c2 + d for d in range(DOF)]
+ [pl.c3 + d for d in range(DOF)]
+ [pl.c4 + d for d in range(DOF)]
for pl in self.str_elements])
self._ops = (np.array(Bm), np.array(Sv_x), np.array(Sv_y),
np.array(wts), dofs)
return self._ops
def stringer_membrane_forces(self, u):
"""Membrane forces `N_{xx}, N_{yy}, N_{xy}` of the stringer elements
at their integration points, shape ``(num_elements, nint**2, 3)``,
from the displacements ``u`` of all the degrees of freedom"""
Bm, Sv_x, Sv_y, wts, dofs = self._stringer_operators()
A = self.str_prop.ABD[:3, :3]
eps = np.einsum('pij,ej->epi', Bm, u[dofs])
return eps @ A.T
def KG_inplane(self, N):
r"""Geometric stiffness matrix of the radial displacement `V` of the
stringers in their plane, for the membrane forces ``N`` of
:meth:`stringer_membrane_forces`:
.. math::
\int (N_{xx} V_{,X}^2 + 2 N_{xy} V_{,X} V_{,Y} + N_{yy}
V_{,Y}^2) \, dX \, dY
These terms, the counterpart of the terms of `w` of the skin, are
not in the geometric stiffness matrix of the BFSC plate element,
which has only the terms of its out-of-plane displacement `W`."""
Bm, Sv_x, Sv_y, wts, dofs = self._stringer_operators()
Nw = N*wts[None, :, None]
Ke = (np.einsum('ep,pi,pj->eij', Nw[:, :, 0], Sv_x, Sv_x)
+ np.einsum('ep,pi,pj->eij', Nw[:, :, 1], Sv_y, Sv_y)
+ np.einsum('ep,pi,pj->eij', Nw[:, :, 2], Sv_x, Sv_y)
+ np.einsum('ep,pi,pj->eij', Nw[:, :, 2], Sv_y, Sv_x))
r = np.repeat(dofs, dofs.shape[1], axis=1).ravel()
c = np.tile(dofs, (1, dofs.shape[1])).ravel()
return coo_matrix((Ke.ravel(), (r, c)),
shape=(self.N, self.N)).tocsr()
def axial_force_membrane(self, eps):
"""Axial force of the whole cylinder for the membrane state of the
uniform axial strain ``eps``"""
return self.E*eps*(2*pi*self.R*self.h
+ self.nstringers*self.bs*self.Hs)
# boundary conditions, as boolean arrays of the reduced dofs
def _mark(self, bk, nodes, dofs):
for n in np.atleast_1d(nodes).ravel():
for d in dofs:
rd = self.red[DOF*n + d]
if rd >= 0:
bk[rd] = True
def constrained(self, ends='SS'):
r"""Constrained reduced degrees of freedom
``ends='SS'``: `v = w = 0` at both ends, `u` and `w_{,x}` free, and
`V = W = W_{,Y} = 0` along the end edges of the stringers.
``ends='CC'``: `u = v = w = w_{,x} = 0` at both ends and the end
edges of the stringers clamped. The derivatives along the edges of
the constrained displacements are constrained as well. In the half
model, the planes `\theta = 0` and `\theta = \pi` are planes of
symmetry: `v = v_{,x} = u_{,y} = w_{,y} = w_{,xy} = 0`.
"""
bk = np.zeros(self.Nr, dtype=bool)
skin, strn = self.skin_nodes, self.str_nodes
ends_skin = np.concatenate((skin[0, :], skin[-1, :]))
ends_str = np.concatenate((strn[:, 0, :], strn[:, -1, :]), axis=None)
if ends == 'SS':
self._mark(bk, ends_skin, [V, V_Y, W, W_Y])
self._mark(bk, ends_str, [V, V_Y, W, W_Y])
elif ends == 'CC':
self._mark(bk, ends_skin, [U, U_Y, V, V_Y, W, W_Y, W_X, W_XY])
self._mark(bk, ends_str, [U, U_Y, V, V_Y, W, W_Y, W_X, W_XY])
else:
raise ValueError(ends)
if self.half:
sym = np.concatenate((skin[:, 0], skin[:, -1]))
self._mark(bk, sym, [V, V_X, U_Y, W_Y, W_XY])
return bk
def axial_dofs(self, i):
"""Reduced degrees of freedom of `u` and `U` of the cross-section of
index ``i`` along `x`, and of their derivatives along the edge"""
bk = np.zeros(self.Nr, dtype=bool)
nodes = np.concatenate((self.skin_nodes[i, :],
self.str_nodes[:, i, :].ravel()))
self._mark(bk, nodes, [U])
bkd = np.zeros(self.Nr, dtype=bool)
self._mark(bkd, nodes, [U_Y])
return bk, bkd
def skin_field(self, u, dof=W):
"""Field of one degree of freedom of the skin, shape ``(nx, nyn)``,
from the displacements ``u`` of all the degrees of freedom"""
return u[DOF*self.skin_nodes + dof]
def wave_numbers(model, u, mmax=12, nmax=40):
r"""Axial half-waves `m` and circumferential waves `n` of the normal
displacement `w` of the skin, from the largest coefficient of its
expansion in `\sin(m \pi x/L) \cos(n \theta)` (`\theta` from a plane of
symmetry), for simply supported ends, or in `\cos(n \theta)` and the
number of sign changes along `x` of the dominant term, for other ends"""
w = model.skin_field(u, W)
x = model.xlin
th = model.theta
if model.half:
# even about theta = 0: cosine terms
wy = np.ones_like(th)
wy[0] = wy[-1] = 0.5
ns = np.arange(nmax + 1)
Cn = (w*wy) @ np.cos(np.outer(th, ns))
else:
Cn = np.abs(np.fft.rfft(w, axis=1))[:, :nmax + 1]
n = int(np.argmax(np.sum(Cn**2, axis=0)))
wn = Cn[:, n]
wx = np.ones_like(x)
wx[0] = wx[-1] = 0.5
ms = np.arange(1, mmax + 1)
Sm = (wn*wx) @ np.sin(np.outer(x, ms)*pi/x[-1])
m = int(ms[np.argmax(np.abs(Sm))])
return m, n
Model of the post-buckling of a cylinder#
The tutorial Post-buckling of an axially compressed cylinder with imperfections traces with
the Riks method of structsolve 0.10.0 or newer the post-buckling path of
the axially compressed cylinder of Simo, Fox and Rifai (1990), an
L/10 x 15 deg section meshed with the Donnell-type element, with an
initial geometric imperfection added exactly. It requires structsolve,
composites and Matplotlib besides bfsccylinder, and builds its model with
the module simo1990_cylinder.py of the same folder, also used by
tests/test_simo1990_cylinder.py:
r"""Axially compressed cylinder of Simo, Fox and Rifai (1990), with imperfections
Model of the tutorial ``simo1990_cylinder_postbuckling.ipynb`` and of
``tests/test_simo1990_cylinder.py``: Section 6.3.4 of Simo, Fox and Rifai
(1990), "Post-buckling analysis of a simply-supported cylinder; effect of
initial imperfections", pp. 54-57, Fig. 6.3.4 on p. 56. A cylinder of radius
``R = 36``, length ``L = 100`` and thickness ``t = 0.125``, of an isotropic
material with ``E = 1e7`` and ``nu = 0.3``, under axial compression. As in
the paper, after Stanley (1985), an ``L/10 x 15 deg`` section is modelled,
with symmetry conditions on its interior edges and a simple support on the
loaded edge, such that the diamond mode of 5 axial half-waves and 12
circumferential waves is the lowest buckling mode.
The section is meshed with the Donnell-type element :class:`.BFSCCylinder`,
with ``x`` from ``0``, the loaded edge, to ``a = L/10``, and ``y`` from ``0``
to ``b = 15 deg x R``:
- ``x = 0``, simply supported and loaded: ``v = v_y = w = w_y = 0``, ``u``
free, under the compressive force per unit length ``lbd*N_ref``;
- ``x = a``, symmetry plane of the axial half-wave: ``u = u_y = v_x = w_x =
w_xy = 0``;
- ``y = 0`` and ``y = b``, symmetry planes of the circumferential waves:
``u_y = v = v_x = w_y = w_xy = 0``.
The elements take no initial geometric imperfection. For the Donnell-type
kinematics it is added exactly by :class:`ImperfectCylinder`, from the
vector ``u0`` of the degrees of freedom of the imperfection ``w0``. The
membrane strains of the imperfect shell are ``eps(u + u0) - eps(u0)`` and
its curvatures ``kappa(u)``, see the tutorial, which gives
``fint_imp(u) = fint(u + u0) - fint(u0) - KG(N0) u``
``KT_imp(u) = KT(u + u0) - KG(N0)``
with ``KG(N0)`` the geometric stiffness matrix of the membrane resultants
``N0`` of the whole non-linear strains of ``u0``. :func:`.update_KG` builds
its matrix from the resultants of the linear strains only, ``KG(u0)``, so
the non-linear part ``KGNL0`` is added by :meth:`SectionModel.calc_kG_N`,
with the rows of ``w_x`` and ``w_y`` of the Donnell-type kinematics. The model
is therefore restricted to the Donnell-type element: with the Sanders-type
element of :mod:`bfsccylinder.bfsccylinder_sanders` the rotation is ``-w_y +
v/R``, ``KG(N0)`` has rows in ``v``, and the construction, not verified here,
would need that matrix.
"""
import time
import numpy as np
from numpy import isclose
from scipy.sparse import coo_matrix, csc_matrix
from composites import isotropic_plate
from structsolve import Analysis, lb, solve
from structsolve.stopping import value_reached
from bfsccylinder import (BFSCCylinder, update_KC0, update_KCNL, update_KG,
update_fint, DOF, DOUBLE, INT, KC0_SPARSE_SIZE, KCNL_SPARSE_SIZE,
KG_SPARSE_SIZE)
from bfsccylinder.quadrature import get_points_weights
from bfsccylinder.utils import assign_constant_ABD
R = 36.
L = 100.
t = 0.125
E = 1.e7
nu = 0.3
#: axial length of the section
A_SEC = L/10.
#: arc length of the 15 deg section
B_SEC = R*np.pi/12.
#: classical buckling load per unit length of the circumference
N_CL = E*t**2/(R*np.sqrt(3*(1 - nu**2)))
# Values read from Fig. 6.3.4 of Simo, Fox and Rifai (1990), p. 56, on the
# page rendered at 300 dpi with a grid of 0.001 x 0.05 drawn over the axes
# (readings of the benchmark simo1990_cylinder.py of structsolve, checked
# against the figure). Reading precision about +-0.0003 in the axial
# deflection and +-0.01 in the load normalized by the linear buckling load.
# Pairs (maximum axial deflection, load).
#: Fig. (a), solid line, Simo et al., 10 % imperfection
FIG_A_SIMO_10 = [
(0.0125, 0.575), (0.0150, 0.675), (0.0168, 0.715), (0.0158, 0.642),
(0.0151, 0.580), (0.0146, 0.512), (0.0142, 0.445), (0.0143, 0.399),
(0.0148, 0.366), (0.0158, 0.343), (0.0184, 0.326), (0.0200, 0.330),
(0.0250, 0.344), (0.0300, 0.372), (0.0350, 0.406), (0.0400, 0.459)]
#: Fig. (a), solid line, Simo et al., 50 % imperfection
FIG_A_SIMO_50 = [
(0.0108, 0.411), (0.0116, 0.366), (0.0126, 0.335), (0.0141, 0.319),
(0.0150, 0.315), (0.0200, 0.315), (0.0250, 0.344), (0.0300, 0.370),
(0.0350, 0.401), (0.0400, 0.449)]
#: Fig. (a), triangles, 9-ANS element of Stanley (1985), 10 % imperfection
FIG_A_ANS_10 = [
(0.0071, 0.331), (0.0095, 0.445), (0.0125, 0.570), (0.0148, 0.666),
(0.0159, 0.717), (0.0166, 0.713), (0.0163, 0.693), (0.0157, 0.653),
(0.0152, 0.619), (0.0148, 0.580), (0.0144, 0.532), (0.0141, 0.487),
(0.0140, 0.430), (0.0143, 0.382), (0.0157, 0.346), (0.0184, 0.319),
(0.0218, 0.324), (0.0296, 0.349), (0.0353, 0.395)]
#: Fig. (a), diamonds, 9-ANS element of Stanley (1985), 50 % imperfection
FIG_A_ANS_50 = [
(0.0055, 0.229), (0.0066, 0.273), (0.0077, 0.310), (0.0084, 0.342),
(0.0089, 0.366), (0.0094, 0.383), (0.0098, 0.397), (0.0106, 0.402),
(0.0109, 0.385), (0.0112, 0.367), (0.0118, 0.354), (0.0122, 0.332),
(0.0132, 0.312), (0.0151, 0.300), (0.0175, 0.297), (0.0204, 0.300),
(0.0237, 0.318), (0.0265, 0.349), (0.0363, 0.420)]
#: Fig. (d), solid line, Simo et al., perfect cylinder, branch switching
FIG_D_SIMO = [
(0.0183, 0.833), (0.0182, 0.742), (0.0171, 0.695), (0.0160, 0.599),
(0.0151, 0.530), (0.0149, 0.457), (0.0160, 0.360), (0.0190, 0.335),
(0.0200, 0.334), (0.0250, 0.355), (0.0300, 0.397), (0.0350, 0.436),
(0.0400, 0.480)]
#: limit loads read from Fig. 6.3.4(a), (maximum axial deflection, load)
FIG_A_LIMIT = {0.1: (0.0168, 0.715), 0.5: (0.0108, 0.411)}
#: slope of the straight line of the figures: axial deflection at load 1
FIG_LINEAR_U1 = 0.0205
class SectionModel(object):
r"""Mesh, boundary conditions and reference load of the section
Parameters
----------
nx, ny : int
Number of elements along ``x`` (axial) and ``y`` (circumferential).
nint : int
Number of Gauss-Legendre points in each direction of an element.
load : float
Compressive force per unit length applied at ``x = 0``.
Attributes
----------
elements : list of :class:`.BFSCCylinder`
N : int
Number of degrees of freedom of the model.
bu : array of bool
The unknown (free) degrees of freedom.
KC0 : csc_matrix
Linear constitutive stiffness matrix, all degrees of freedom.
fext : array
Reference external force vector, all degrees of freedom.
edofs : array of int, shape ``(num_elements, 40)``
Global degrees of freedom of each element.
loaded : array of int
Node indices of the loaded edge ``x = 0``, in increasing ``y``.
"""
def __init__(self, nx=8, ny=8, nint=4, load=N_CL):
self.nx, self.ny, self.nint = nx, ny, nint
self.points, self.weights = get_points_weights(nint=nint)
self.prop = isotropic_plate(thickness=t, E=E, nu=nu)
xlin = np.linspace(0, A_SEC, nx + 1)
ylin = np.linspace(0, B_SEC, ny + 1)
xmesh, ymesh = np.meshgrid(xlin, ylin, indexing='ij')
self.x = x = xmesh.ravel()
self.y = y = ymesh.ravel()
nids = np.arange(x.size).reshape(nx + 1, ny + 1)
n1s = nids[:-1, :-1].ravel()
n2s = nids[1:, :-1].ravel()
n3s = nids[1:, 1:].ravel()
n4s = nids[:-1, 1:].ravel()
self.N = N = DOF*x.size
self.lex = A_SEC/nx
self.ley = B_SEC/ny
self.elements = []
for i, (n1, n2, n3, n4) in enumerate(zip(n1s, n2s, n3s, n4s)):
shell = BFSCCylinder(nint)
shell.n1, shell.n2, shell.n3, shell.n4 = n1, n2, n3, n4
shell.c1, shell.c2 = DOF*n1, DOF*n2
shell.c3, shell.c4 = DOF*n3, DOF*n4
shell.R = R
shell.lex = self.lex
shell.ley = self.ley
assign_constant_ABD(shell, self.prop)
shell.init_k_KC0 = i*KC0_SPARSE_SIZE
shell.init_k_KCNL = i*KCNL_SPARSE_SIZE
shell.init_k_KG = i*KG_SPARSE_SIZE
self.elements.append(shell)
ne = len(self.elements)
self.edofs = np.array([np.concatenate([DOF*n + np.arange(DOF)
for n in nodes])
for nodes in zip(n1s, n2s, n3s, n4s)])
KC0r = np.zeros(KC0_SPARSE_SIZE*ne, dtype=INT)
KC0c = np.zeros(KC0_SPARSE_SIZE*ne, dtype=INT)
KC0v = np.zeros(KC0_SPARSE_SIZE*ne, dtype=DOUBLE)
for shell in self.elements:
update_KC0(shell, self.points, self.weights, KC0r, KC0c, KC0v)
self.KC0 = coo_matrix((KC0v, (KC0r, KC0c)), shape=(N, N)).tocsc()
self.KCNLr = np.zeros(KCNL_SPARSE_SIZE*ne, dtype=INT)
self.KCNLc = np.zeros(KCNL_SPARSE_SIZE*ne, dtype=INT)
self.KCNLv = np.zeros(KCNL_SPARSE_SIZE*ne, dtype=DOUBLE)
self.KGr = np.zeros(KG_SPARSE_SIZE*ne, dtype=INT)
self.KGc = np.zeros(KG_SPARSE_SIZE*ne, dtype=INT)
self.KGv = np.zeros(KG_SPARSE_SIZE*ne, dtype=DOUBLE)
# boundary conditions, dofs u, u_x, u_y, v, v_x, v_y, w, w_x, w_y, w_xy
bk = np.zeros(N, dtype=bool)
x0 = isclose(x, 0)
xa = isclose(x, A_SEC)
ye = isclose(y, 0) | isclose(y, B_SEC)
for dof in (3, 5, 6, 8):
bk[dof::DOF] |= x0
for dof in (0, 2, 4, 7, 9):
bk[dof::DOF] |= xa
for dof in (2, 3, 4, 8, 9):
bk[dof::DOF] |= ye
self.bk = bk
self.bu = ~bk
self.loaded = np.flatnonzero(x0)
self.loaded = self.loaded[np.argsort(y[self.loaded])]
# consistent nodal forces of the uniform force per unit length
# ``load`` along the edge x = 0, compressive, pointing to +x: the
# integrals of the Hermite functions along an element side of length
# ley are ley/2 for the values and +-ley**2/12 for the derivatives
# u_y, which cancel between neighbours and fall on the constrained
# u_y of the corners
self.fext = np.zeros(N, dtype=DOUBLE)
for na, nb in zip(self.loaded[:-1], self.loaded[1:]):
self.fext[DOF*na + 0] += load*self.ley/2
self.fext[DOF*nb + 0] += load*self.ley/2
self.fext[DOF*na + 2] += load*self.ley**2/12
self.fext[DOF*nb + 2] -= load*self.ley**2/12
self.fext[self.bk] = 0
self.load = load
# rows of the interpolation at the integration points, the same for
# every element of this regular mesh: Bm, Bb, Sw_x and Sw_y of shape
# (num_points, ..., 40), and the weights times the Jacobian
shell = self.elements[0]
Bm, Bb, Gx, Gy, wp = [], [], [], [], []
for i, xi in enumerate(self.points):
for j, eta in enumerate(self.points):
shell.update_Bm(xi, eta)
shell.update_Bb(xi, eta)
shell.update_Sw_x(xi, eta)
shell.update_Sw_y(xi, eta)
Bm.append(np.array(shell.Bm))
Bb.append(np.array(shell.Bb))
Gx.append(np.array(shell.Sw_x))
Gy.append(np.array(shell.Sw_y))
wp.append(self.weights[i]*self.weights[j]
*self.lex*self.ley/4)
self.Bm, self.Bb = np.array(Bm), np.array(Bb)
self.Gx, self.Gy = np.array(Gx), np.array(Gy)
self.wp = np.array(wp)
p = self.prop
self.A = np.array([[p.A11, p.A12, p.A16], [p.A12, p.A22, p.A26],
[p.A16, p.A26, p.A66]])
self.B = np.array([[p.B11, p.B12, p.B16], [p.B12, p.B22, p.B26],
[p.B16, p.B26, p.B66]])
self.D = np.array([[p.D11, p.D12, p.D16], [p.D12, p.D22, p.D26],
[p.D16, p.D26, p.D66]])
def expand(self, c):
"""The displacements of all dofs from those of the unknown dofs"""
u = np.zeros(self.N, dtype=DOUBLE)
u[self.bu] = c
return u
def restrict(self, K):
"""``K`` on the unknown dofs"""
return csc_matrix(K)[self.bu, :][:, self.bu]
def calc_fint(self, u):
"""Internal forces of the displacements ``u``, all dofs"""
fint = np.zeros(self.N, dtype=DOUBLE)
for shell in self.elements:
update_fint(u, shell, self.points, self.weights, fint)
return fint
def calc_KCNL(self, u):
"""``KCNL`` of :func:`.update_KCNL` at ``u``, all dofs"""
self.KCNLv *= 0
for shell in self.elements:
update_KCNL(u, shell, self.points, self.weights, self.KCNLr,
self.KCNLc, self.KCNLv)
return coo_matrix((self.KCNLv, (self.KCNLr, self.KCNLc)),
shape=(self.N, self.N)).tocsc()
def calc_KG(self, u):
"""``KG`` of :func:`.update_KG` at ``u``, all dofs: the membrane
resultants of the linear strains only"""
self.KGv *= 0
for shell in self.elements:
update_KG(u, shell, self.points, self.weights, self.KGr,
self.KGc, self.KGv)
return coo_matrix((self.KGv, (self.KGr, self.KGc)),
shape=(self.N, self.N)).tocsc()
def generalized_strains(self, u):
"""Linear membrane strains, non-linear membrane strains and
curvatures of ``u`` at the integration points, each of shape
``(num_elements, num_points, 3)``"""
ue = u[self.edofs]
epsL = np.einsum('pkj,ej->epk', self.Bm, ue)
kappa = np.einsum('pkj,ej->epk', self.Bb, ue)
wx = np.einsum('pj,ej->ep', self.Gx, ue)
wy = np.einsum('pj,ej->ep', self.Gy, ue)
epsNL = np.stack((wx**2/2, wy**2/2, wx*wy), axis=-1)
return epsL, epsNL, kappa
def membrane_resultants(self, u, nonlinear=True):
"""``(Nxx, Nyy, Nxy)`` of ``u`` at the integration points, of shape
``(num_elements, num_points, 3)``, with or without the non-linear
strains"""
epsL, epsNL, kappa = self.generalized_strains(u)
eps = epsL + epsNL if nonlinear else epsL
return eps @ self.A.T + kappa @ self.B.T
def calc_kG_N(self, N):
r"""Geometric stiffness matrix of the membrane resultants ``N`` given
at the integration points, shape ``(num_elements, num_points, 3)``:
`\int_A \boldsymbol{G}^\top [N_{xx}\ N_{xy};\ N_{xy}\ N_{yy}]
\boldsymbol{G} dA`, all dofs, with the quadrature of the elements"""
Nw = N*self.wp[None, :, None]
Gx, Gy = self.Gx, self.Gy
ke = (np.einsum('ep,pi,pj->eij', Nw[..., 0], Gx, Gx)
+ np.einsum('ep,pi,pj->eij', Nw[..., 1], Gy, Gy)
+ np.einsum('ep,pi,pj->eij', Nw[..., 2], Gx, Gy)
+ np.einsum('ep,pi,pj->eij', Nw[..., 2], Gy, Gx))
rows = np.repeat(self.edofs[:, :, None], 40, axis=2)
cols = np.repeat(self.edofs[:, None, :], 40, axis=1)
return coo_matrix((ke.ravel(), (rows.ravel(), cols.ravel())),
shape=(self.N, self.N)).tocsc()
def strain_energy(self, u, u0=None):
"""Strain energy of the strains ``eps(u + u0) - eps(u0)`` and the
curvatures ``kappa(u)``, by the quadrature of the elements,
independently of the kernels of the elements"""
if u0 is None:
u0 = np.zeros_like(u)
eL1, eN1, k1 = self.generalized_strains(u + u0)
eL0, eN0, k0 = self.generalized_strains(u0)
e = (eL1 + eN1) - (eL0 + eN0)
k = k1 - k0
dens = (np.einsum('epi,ij,epj->ep', e, self.A, e)
+ 2*np.einsum('epi,ij,epj->ep', e, self.B, k)
+ np.einsum('epi,ij,epj->ep', k, self.D, k))
return 0.5*np.sum(dens*self.wp[None, :])
def shortening(self, u):
"""Maximum and average axial displacement of the loaded edge ``x =
0``, the shortening of the section: the nodal maximum and the
integral of the Hermite interpolation of ``u`` over the edge"""
n = self.loaded
ui = u[DOF*n + 0]
uyi = u[DOF*n + 2]
l = self.ley
integral = np.sum(l/2*(ui[:-1] + ui[1:])
+ l**2/12*(uyi[:-1] - uyi[1:]))
return ui.max(), integral/B_SEC
def w_max(self, u):
"""Largest nodal ``|w|``"""
return np.abs(u[6::DOF]).max()
class ImperfectCylinder(object):
r"""Callables of structsolve.Analysis for the section with the initial
imperfection ``u0``
``u0`` has all dofs of the model, with only the dofs of `w` non-zero,
``u0 = 0`` gives the perfect shell. The callables work on the unknown
dofs ``c = u[bu]``:
- ``calc_fint(c) = fint(u + u0) - fint(u0) - KG(u0) u - KGNL0 u``
- ``calc_kC(c, NLgeom=True) = KC0 + KCNL(u + u0) - KGNL0``
- ``calc_kG(c, NLgeom=True) = KG(u + u0) - KG(u0) = KG(u)``
with ``KG`` of :func:`.update_KG`, which is linear in ``u``, and
``KGNL0`` the geometric stiffness matrix of the membrane resultants of
the non-linear strains of ``u0``, such that ``KG(u0) + KGNL0 = KG(N0)``.
``calc_kC + calc_kG`` is the Jacobian of ``calc_fint``, and ``calc_kG``
keeps the meaning of the geometric stiffness matrix of the perfect
shell. ``calc_kC()``, without ``NLgeom``, gives the linear stiffness
matrix of the imperfect shell, ``KC0 + KCNL(u0) - KGNL0 + KG(u0)``.
"""
def __init__(self, model, u0=None):
self.model = m = model
self.u0 = np.zeros(m.N) if u0 is None else np.asarray(u0, float)
self.fint0 = m.calc_fint(self.u0)
N0NL = (m.membrane_resultants(self.u0, nonlinear=True)
- m.membrane_resultants(self.u0, nonlinear=False))
self.KGNL0 = m.calc_kG_N(N0NL)
self.KG0 = m.calc_KG(self.u0)
self.KC0uu = m.restrict(m.KC0)
self.K0uu = m.restrict(m.KC0 + m.calc_KCNL(self.u0) - self.KGNL0
+ self.KG0)
self.fextuu = m.fext[m.bu]
def calc_fext(self, inc=1., silent=True):
return inc*self.fextuu
def calc_fint(self, c, silent=True):
m = self.model
u = m.expand(c)
f = (m.calc_fint(u + self.u0) - self.fint0 - self.KG0 @ u
- self.KGNL0 @ u)
return f[m.bu]
def calc_kC(self, c=None, NLgeom=False, silent=True):
if c is None or not NLgeom:
return self.K0uu
m = self.model
u = m.expand(c)
return m.restrict(m.KC0 + m.calc_KCNL(u + self.u0) - self.KGNL0)
def calc_kG(self, c=None, NLgeom=False, silent=True):
m = self.model
u = m.expand(c)
return m.restrict(m.calc_KG(u))
def linear_buckling(model, num_eigvalues=4):
"""Linear buckling load factors of the reference load of ``model`` and
the modes, all dofs, from the linear static solution"""
KC0uu = model.restrict(model.KC0)
clin = solve(KC0uu, model.fext[model.bu], silent=True)
KGuu = model.restrict(model.calc_KG(model.expand(clin)))
eigvals, eigvecs = lb(KC0uu, KGuu, num_eigvalues=num_eigvalues,
silent=True)
modes = np.zeros((model.N, eigvecs.shape[1]))
modes[model.bu] = np.real(eigvecs)
return np.real(eigvals), modes, model.expand(clin)
def imperfection(model, mode, xi):
"""The dofs of `w` of ``mode``, scaled to a largest nodal ``|w|`` of
``xi*t``"""
u0 = np.zeros_like(mode)
for dof in (6, 7, 8, 9):
u0[dof::DOF] = mode[dof::DOF]
w = u0[6::DOF]
return u0*xi*t/w[np.argmax(np.abs(w))]
def normalized_model(nx, ny, nint=4):
"""Model whose reference load is its linear buckling load, the lowest
buckling mode and the linear buckling load in units of ``N_CL``"""
model = SectionModel(nx, ny, nint)
eigvals, modes, _ = linear_buckling(model)
lbd_cr = eigvals[0]
model = SectionModel(nx, ny, nint, load=lbd_cr*N_CL)
return model, modes[:, 0], lbd_cr
def run_riks(xi, nx=8, ny=8, nint=4, umax=0.04, branch_switching=0,
silent=True, **kwargs):
r"""Riks analysis of the section until the maximum shortening ``umax``
``xi > 0``: imperfection of amplitude ``xi*t`` in the shape of the lowest
linear buckling mode; ``xi = 0``: the perfect shell, with ``inertia``,
SuperLU, and ``branch_switching``. ``kwargs`` are attributes of the
:class:`structsolve.Analysis`. Returns a dict with the load factors
``lbds`` (in units of the linear buckling load), the maximum and
average shortenings ``umax`` and ``uavg``, the largest ``|w|``
``wmax``, starting from the unloaded state, the linear buckling load
``lbd_cr`` in units of ``N_CL``, the run time ``elapsed`` and the
analysis ``an``.
"""
model, mode, lbd_cr = normalized_model(nx, ny, nint)
u0 = imperfection(model, mode, xi) if xi > 0 else np.zeros(model.N)
imp = ImperfectCylinder(model, u0)
an = Analysis(imp.calc_fext, imp.calc_fint, imp.calc_kC, imp.calc_kG)
an.NL_method = 'arc_length_riks'
an.final_load_factor = None
an.initialInc = 0.1
an.maxInc = 0.1
an.minInc = 1e-5
an.maxArcLength = 1e4
if xi == 0:
an.inertia = True
an.branch_switching = branch_switching
an.linear_solver = 'splu'
an.initialInc = an.maxInc = 0.03
an.stop_criterion = value_reached(
lambda c: model.shortening(model.expand(c))[0], umax)
for k, v in kwargs.items():
setattr(an, k, v)
t0 = time.perf_counter()
an.static(NLgeom=True, silent=silent)
elapsed = time.perf_counter() - t0
us = [model.expand(c) for c in an.cs]
sh = np.array([[0., 0.]] + [model.shortening(u) for u in us])
return dict(lbds=np.concatenate(([0.], an.increments)),
umax=sh[:, 0], uavg=sh[:, 1],
wmax=np.array([0.] + [model.w_max(u) for u in us]),
lbd_cr=lbd_cr, elapsed=elapsed, an=an, model=model, u0=u0)
def limit_point(lbds, u):
"""Load and shortening at the first limit point, the first decrease of
the load factor, ``None`` if there is none"""
drops = np.flatnonzero(np.diff(lbds) < 0)
if drops.size == 0:
return None
i = drops[0]
return lbds[i], u[i]
Model of the cantilever cylinder under a follower pressure#
The tutorial Cantilever cylinder under a follower pressure compares
the linear buckling factors of the cantilever cylinders of Schweizerhof and
Ramm (1984) under a dead and a follower external pressure, Table 4, and
traces with the Riks method of structsolve 0.10.0 or newer the path of the
imperfect shell 1 of Fig. 13, with the full, unsymmetric, load stiffness
matrix of update_KCfollower() and the count of the negative
eigenvalues of the unsymmetric tangent stiffness matrix. The sectors of the
paper are meshed with the Sanders-type element. It requires structsolve,
composites and Matplotlib besides bfsccylinder, and builds its model with
the module schweizerhof1984_cylinder.py of the same folder, also used
by tests/test_schweizerhof1984_cylinder.py:
r"""Cantilever cylinder under a follower pressure, Schweizerhof and Ramm (1984)
Model of the tutorial ``schweizerhof1984_cylinder_follower.ipynb`` and of
``tests/test_schweizerhof1984_cylinder.py``: Section 6.4 of Schweizerhof and
Ramm (1984), "Cylindrical cantilever shell under uniform pressure", Fig. 12
on p. 1110, Table 4 on p. 1111 and Fig. 13 on p. 1112. A cylinder of radius
``R``, length ``L`` and thickness ``t``, of an isotropic material with ``E =
2.1e6`` and ``mu = 0.3``, upper edge free, lower edge hinged, under a
uniform external pressure, either of constant direction (dead, index ``c``
of the paper) or normal to the deformed surface (follower, index ``f``).
The three shells of the paper are given by :data:`SHELLS`, with the number
``n`` of circumferential waves of their buckling mode. The pressure is
normalized by eq. (60) of the paper,
``p = k E/(12 (1 - mu**2)) (t/R)**3``
and the reference pressure of the models is that of ``k = 1``, such that the
load multipliers are the factors ``k`` of Table 4 and the load ``K`` of
Fig. 13. Here ``R = 100``, the factors do not depend on it.
As in the paper, a sector of the angle ``180 deg/n`` (60, 45 and 30 deg for
the shells 1, 2 and 3) is modelled, with ``x`` from ``0``, the lower edge, to
``L``, the free upper edge, and ``y`` from ``0`` to ``b = pi R/n``, with
symmetry planes at ``y = 0`` and ``y = b``, which keep the modes with
``n``, ``2n``, ... circumferential waves; ``full=True`` gives the closed
cylinder instead, with every ``n``. The degrees of freedom of a node are
``u, u_x, u_y, v, v_x, v_y, w, w_x, w_y, w_xy``:
- ``x = 0``, B1, "warping completely restrained", ``v = w = u = m_x = 0``:
``u = u_y = v = v_y = w = w_y = 0``;
- ``x = 0``, B2, "warping partially restrained", ``v = w = u0 = m_x = 0``:
``v = v_y = w = w_y = 0`` and ``u = 0`` at the node ``x = y = 0`` only,
which removes the axial translation of the sector (sector only);
- ``y = 0`` and ``y = b``, symmetry planes: ``u_y = v = v_x = w_y = w_xy =
0``;
- ``x = L``: free.
The element is the Sanders-type :class:`.BFSCCylinderSanders` by default,
``kinematics='donnell'`` gives the Donnell-type :class:`.BFSCCylinder`. With
``n = 3`` circumferential waves over a length of six radii the Donnell-type
kinematics are not accurate, see the tutorial. The follower pressure is that
of :func:`.update_fext_follower` with the first-order area vector
(``quadratic=0``), with its full, unsymmetric, load stiffness matrix of
:func:`.update_KCfollower`: on the free edge the axial displacement is not
restrained, which makes the pressure nonconservative.
The initial geometric imperfection ``u0``, of the dofs of `w` only, is added
exactly as in ``simo1990_cylinder.py``: the membrane strains of the
imperfect shell are ``eps(u + u0) - eps(u0)`` and its curvatures
``kappa(u)``, which gives
``fint_imp(u) = fint(u + u0) - fint(u0) - KG(N0) u``
with ``KG(N0)`` the geometric stiffness matrix of the membrane resultants
``N0`` of the strains of ``u0``. The derivation holds for both kinematics,
whose non-linear membrane strains are quadratic forms of the rotations
``G u``, with ``G`` of the rows ``Sw_x`` and ``Sw_y`` (Donnell) or ``Sw_x``
and ``Sw_y - Sv/R`` (Sanders), see :meth:`CantileverModel.calc_kG_N`. The
pressure acts on the imperfect surface: the follower force vector is that of
the displacements ``u + u0``.
"""
import time
import numpy as np
from numpy import isclose, pi
from scipy.sparse import coo_matrix, csc_matrix
from composites import isotropic_plate
from structsolve import Analysis, lb, solve
import bfsccylinder
import bfsccylinder.sanders
from bfsccylinder.quadrature import get_points_weights
from bfsccylinder.utils import assign_constant_ABD
DOF = 10
E = 2.1e6
MU = 0.3
R = 100.
#: shells 1, 2 and 3 of Section 6.4: ``(L/R, R/t, n)``
SHELLS = {1: (6., 100., 3), 2: (3., 100., 4), 3: (3., 500., 6)}
KINEMATICS = {
'donnell': (bfsccylinder, 'BFSCCylinder'),
'sanders': (bfsccylinder.sanders, 'BFSCCylinderSanders'),
}
# Table 4 of Schweizerhof and Ramm (1984), p. 1111, read on the page image.
# The k_f of B2 are those of the load stiffness matrix symmetrized by eqs.
# (63), (64) and (65); the paper gives a single k_f for B1 (none for shell 3).
#: ``TABLE4[(shell, bc)] = (k_c, k_f)``, ``k_f`` a tuple of (63), (64), (65)
TABLE4 = {
(1, 'B2'): (9.11, (8.10, 8.09, 8.12)),
(1, 'B1'): (10.88, (9.70,)),
(2, 'B2'): (16.43, (15.42, 15.38, 15.46)),
(2, 'B1'): (20.85, (19.58,)),
(3, 'B2'): (36.48, (35.48, 35.46, 35.50)),
(3, 'B1'): (45.83, ()),
}
#: Table 4, clamped-free shells, k_f of refs [36] and [15] of the paper
TABLE4_CLAMPED = {1: (9.85, 10.51), 2: (19.97, 21.01), 3: (None, 46.98)}
# Fig. 13 of Schweizerhof and Ramm (1984), p. 1112, imperfect shell 1,
# follower pressure, imperfection of 1 % of t in the shape of the linear
# buckling mode. Read on the page rendered at 300 dpi, from the centre line
# of the thick curve, with the axes calibrated on their ticks: 87.0 px per
# unit of K and 268 px per unit of w/t. The curve is about 6 px thick, the
# reading precision is about +-0.05 in K and +-0.01 in w/t (+-0.003 below
# w/t = 0.2, where the curve is steep). Pairs (w/t, K).
#: Fig. 13, imperfect shell, maximum displacement against K
FIG13_IMPERFECT = [
(0.058, 6.72), (0.095, 7.19), (0.133, 7.41), (0.170, 7.54),
(0.245, 7.71), (0.319, 7.77), (0.506, 7.85), (0.693, 7.91),
(1.002, 7.93), (1.439, 7.95), (2.007, 7.99), (2.559, 8.00),
(2.970, 8.01)]
#: Fig. 13, horizontal line of the perfect shell, K
FIG13_PERFECT = 8.06
#: Fig. 13, dashed curve, current stiffness parameter of Bergan, pairs
#: (CS, K), read on the dashes to about +-0.02 in CS and +-0.05 in K
FIG13_CS = [
(0.97, 0.98), (0.948, 1.44), (0.93, 1.90), (0.88, 2.36), (0.81, 2.82),
(0.747, 3.28), (0.667, 3.74), (0.57, 4.20), (0.475, 4.66),
(0.37, 5.11), (0.26, 5.57), (0.18, 6.03), (0.10, 6.49), (0.047, 6.95),
(0.02, 7.41)]
def shell_geometry(shell):
"""``(L, t, n)`` of the shell ``1``, ``2`` or ``3``, with ``R = 100``"""
LR, Rt, n = SHELLS[shell]
return LR*R, R/Rt, n
def p_ref(t):
"""Pressure of ``k = 1`` in eq. (60), positive"""
return E/(12*(1 - MU**2))*(t/R)**3
class CantileverModel(object):
r"""Mesh, boundary conditions and matrices of the cantilever cylinder
Parameters
----------
shell : int
``1``, ``2`` or ``3``, see :data:`SHELLS`.
bc : str
``'B1'`` or ``'B2'``, the boundary condition of the lower edge.
nx, ny : int
Number of elements along ``x`` (axial) and ``y`` (circumferential).
kinematics : str
``'sanders'`` or ``'donnell'``.
n : int or None
Number of circumferential waves of the sector of ``180 deg/n``,
``None`` for that of :data:`SHELLS`.
full : bool
The closed cylinder instead of the sector, only with ``'B1'``.
nint : int
Number of Gauss-Legendre points in each direction of an element.
Attributes
----------
L, t, b : float
Length, thickness and arc length of the model.
N : int
Number of degrees of freedom.
bu : array of bool
The unknown (free) degrees of freedom.
KC0, KF : csc_matrix
Linear constitutive stiffness matrix and load stiffness matrix of
the reference pressure (external, ``k = 1``), all dofs; ``KF`` is
constant with the first-order area vector.
fp0 : array
Force vector of the reference pressure on the undeformed shell, the
dead load, all dofs.
"""
def __init__(self, shell=1, bc='B1', nx=16, ny=4, kinematics='sanders',
n=None, full=False, nint=4):
self.shell, self.bc, self.kinematics = shell, bc, kinematics
self.nx, self.ny, self.nint, self.full = nx, ny, nint, full
self.L, self.t, n0 = shell_geometry(shell)
self.n = n = n0 if n is None else n
self.mod = mod = KINEMATICS[kinematics][0]
cls = getattr(mod, KINEMATICS[kinematics][1])
self.points, self.weights = get_points_weights(nint=nint)
self.prop = isotropic_plate(thickness=self.t, E=E, nu=MU)
self.b = b = 2*pi*R if full else pi*R/n
nyn = ny if full else ny + 1
xlin = np.linspace(0, self.L, nx + 1)
ylin = np.arange(nyn)*b/ny
xmesh, ymesh = np.meshgrid(xlin, ylin, indexing='ij')
self.x = x = xmesh.ravel()
self.y = y = ymesh.ravel()
nids = np.arange(x.size).reshape(nx + 1, nyn)
if full:
nids = np.hstack((nids, nids[:, :1]))
n1s = nids[:-1, :-1].ravel()
n2s = nids[1:, :-1].ravel()
n3s = nids[1:, 1:].ravel()
n4s = nids[:-1, 1:].ravel()
self.N = N = DOF*x.size
self.lex, self.ley = self.L/nx, b/ny
self.elements = []
for i, (n1, n2, n3, n4) in enumerate(zip(n1s, n2s, n3s, n4s)):
elem = cls(nint)
elem.n1, elem.n2, elem.n3, elem.n4 = n1, n2, n3, n4
elem.c1, elem.c2 = DOF*n1, DOF*n2
elem.c3, elem.c4 = DOF*n3, DOF*n4
elem.R = R
elem.lex = self.lex
elem.ley = self.ley
assign_constant_ABD(elem, self.prop)
elem.init_k_KC0 = i*mod.KC0_SPARSE_SIZE
elem.init_k_KCNL = i*mod.KCNL_SPARSE_SIZE
elem.init_k_KG = i*mod.KG_SPARSE_SIZE
elem.init_k_KCfollower = i*mod.KCfollower_SPARSE_SIZE
self.elements.append(elem)
ne = len(self.elements)
self.edofs = np.array([np.concatenate([DOF*nd + np.arange(DOF)
for nd in nodes])
for nodes in zip(n1s, n2s, n3s, n4s)])
INT, DOUBLE = mod.INT, mod.DOUBLE
pts, wts = self.points, self.weights
def arrays(size):
return (np.zeros(size*ne, dtype=INT), np.zeros(size*ne, dtype=INT),
np.zeros(size*ne, dtype=DOUBLE))
r, c, v = arrays(mod.KC0_SPARSE_SIZE)
for elem in self.elements:
mod.update_KC0(elem, pts, wts, r, c, v)
self.KC0 = coo_matrix((v, (r, c)), shape=(N, N)).tocsc()
self.KCNLr, self.KCNLc, self.KCNLv = arrays(mod.KCNL_SPARSE_SIZE)
self.KGr, self.KGc, self.KGv = arrays(mod.KG_SPARSE_SIZE)
# reference pressure, external: negative
self.p = -p_ref(self.t)
self.pp = np.full((nint, nint), self.p)
zero = np.zeros(N)
r, c, v = arrays(mod.KCfollower_SPARSE_SIZE)
self.fp0 = np.zeros(N)
for elem in self.elements:
mod.update_KCfollower(zero, elem, pts, wts, self.pp, r, c, v)
mod.update_fext_follower(zero, elem, pts, wts, self.pp, self.fp0)
self.KF = coo_matrix((v, (r, c)), shape=(N, N)).tocsc()
# boundary conditions
bk = np.zeros(N, dtype=bool)
x0 = isclose(x, 0)
if bc == 'B1':
dofs = (0, 2, 3, 5, 6, 8)
elif bc == 'B2':
if full:
raise NotImplementedError('B2 only for the sector')
dofs = (3, 5, 6, 8)
bk[0] = True # u of the node x = y = 0
else:
raise ValueError(bc)
for dof in dofs:
bk[dof::DOF] |= x0
if not full:
ye = isclose(y, 0) | isclose(y, b)
for dof in (2, 3, 4, 8, 9):
bk[dof::DOF] |= ye
self.bk = bk
self.bu = ~bk
# rows of the interpolation at the integration points, the same for
# every element of this regular mesh, for calc_kG_N and the strains
elem = self.elements[0]
Bm, Bb, Gx, Gy, wp = [], [], [], [], []
for i, xi in enumerate(pts):
for j, eta in enumerate(pts):
elem.update_Bm(xi, eta)
elem.update_Bb(xi, eta)
elem.update_Sw_x(xi, eta)
elem.update_Sw_y(xi, eta)
elem.update_Sv(xi, eta)
Bm.append(np.array(elem.Bm))
Bb.append(np.array(elem.Bb))
Gx.append(np.array(elem.Sw_x))
gy = np.array(elem.Sw_y)
if kinematics == 'sanders':
gy = gy - np.array(elem.Sv)/R
Gy.append(gy)
wp.append(wts[i]*wts[j]*self.lex*self.ley/4)
self.Bm, self.Bb = np.array(Bm), np.array(Bb)
self.Gx, self.Gy = np.array(Gx), np.array(Gy)
self.wp = np.array(wp)
self._rows = np.repeat(self.edofs[:, :, None], 40, axis=2).ravel()
self._cols = np.repeat(self.edofs[:, None, :], 40, axis=1).ravel()
p = self.prop
self.A = np.array([[p.A11, p.A12, p.A16], [p.A12, p.A22, p.A26],
[p.A16, p.A26, p.A66]])
self.B = np.array([[p.B11, p.B12, p.B16], [p.B12, p.B22, p.B26],
[p.B16, p.B26, p.B66]])
self.D = np.array([[p.D11, p.D12, p.D16], [p.D12, p.D22, p.D26],
[p.D16, p.D26, p.D66]])
def expand(self, c):
"""The displacements of all dofs from those of the unknown dofs"""
u = np.zeros(self.N)
u[self.bu] = c
return u
def restrict(self, K):
"""``K`` on the unknown dofs"""
return csc_matrix(K)[self.bu, :][:, self.bu]
def calc_fint(self, u):
"""Internal forces of the displacements ``u``, all dofs"""
fint = np.zeros(self.N)
for elem in self.elements:
self.mod.update_fint(u, elem, self.points, self.weights, fint)
return fint
def calc_fp(self, u, reference=1):
"""Follower force vector of the reference pressure at ``u``, all
dofs; ``reference=0`` leaves out the dead load ``fp0``"""
f = np.zeros(self.N)
for elem in self.elements:
self.mod.update_fext_follower(u, elem, self.points, self.weights,
self.pp, f, 0, reference)
return f
def calc_KCNL(self, u):
"""``KCNL`` of :func:`.update_KCNL` at ``u``, all dofs"""
self.KCNLv *= 0
for elem in self.elements:
self.mod.update_KCNL(u, elem, self.points, self.weights,
self.KCNLr, self.KCNLc, self.KCNLv)
return coo_matrix((self.KCNLv, (self.KCNLr, self.KCNLc)),
shape=(self.N, self.N)).tocsc()
def calc_KG(self, u):
"""``KG`` of :func:`.update_KG` at ``u``, all dofs: the membrane
resultants of the linear strains only"""
self.KGv *= 0
for elem in self.elements:
self.mod.update_KG(u, elem, self.points, self.weights, self.KGr,
self.KGc, self.KGv)
return coo_matrix((self.KGv, (self.KGr, self.KGc)),
shape=(self.N, self.N)).tocsc()
def generalized_strains(self, u):
"""Linear membrane strains, non-linear membrane strains and
curvatures of ``u`` at the integration points, each of shape
``(num_elements, num_points, 3)``"""
ue = u[self.edofs]
epsL = np.einsum('pkj,ej->epk', self.Bm, ue)
kappa = np.einsum('pkj,ej->epk', self.Bb, ue)
gx = np.einsum('pj,ej->ep', self.Gx, ue)
gy = np.einsum('pj,ej->ep', self.Gy, ue)
epsNL = np.stack((gx**2/2, gy**2/2, gx*gy), axis=-1)
return epsL, epsNL, kappa
def membrane_resultants(self, u, nonlinear=True):
"""``(Nxx, Nyy, Nxy)`` of ``u`` at the integration points, with or
without the non-linear strains"""
epsL, epsNL, kappa = self.generalized_strains(u)
eps = epsL + epsNL if nonlinear else epsL
return eps @ self.A.T + kappa @ self.B.T
def calc_kG_N(self, N):
r"""Geometric stiffness matrix of the membrane resultants ``N`` given
at the integration points, shape ``(num_elements, num_points, 3)``:
`\int_A \boldsymbol{G}^\top [N_{xx}\ N_{xy};\ N_{xy}\ N_{yy}]
\boldsymbol{G} dA`, all dofs, with ``G`` of the kinematics"""
Nw = N*self.wp[None, :, None]
Gx, Gy = self.Gx, self.Gy
ke = (np.einsum('ep,pi,pj->eij', Nw[..., 0], Gx, Gx)
+ np.einsum('ep,pi,pj->eij', Nw[..., 1], Gy, Gy)
+ np.einsum('ep,pi,pj->eij', Nw[..., 2], Gx, Gy)
+ np.einsum('ep,pi,pj->eij', Nw[..., 2], Gy, Gx))
rows = np.repeat(self.edofs[:, :, None], 40, axis=2)
cols = np.repeat(self.edofs[:, None, :], 40, axis=1)
return coo_matrix((ke.ravel(), (rows.ravel(), cols.ravel())),
shape=(self.N, self.N)).tocsc()
def calc_KT(self, u):
r"""Tangent stiffness matrix at ``u``, all dofs, ``KC0 + KCNL(u) +
KG(u)``, vectorized over the elements with NumPy: the Hessian of
:meth:`strain_energy` with ``u0 = 0``, the integral of ``B.T A B +
B.T B Bb + Bb.T B B + Bb.T D Bb + G.T N G``, ``B = Bm + BNL(u)``, with
the quadrature of the elements. It gives the matrix of the kernels
:func:`.update_KC0`, :func:`.update_KCNL` and :func:`.update_KG` to
the round-off, see ``tests/test_schweizerhof1984_cylinder.py``, in
a fraction of the time of ``update_KCNL`` of the Sanders-type
element"""
ue = u[self.edofs]
gx = np.einsum('pj,ej->ep', self.Gx, ue)
gy = np.einsum('pj,ej->ep', self.Gy, ue)
Bm = self.Bm[None] + np.stack((
gx[..., None]*self.Gx[None],
gy[..., None]*self.Gy[None],
gx[..., None]*self.Gy[None] + gy[..., None]*self.Gx[None]),
axis=2)
Bb = self.Bb
Nr = self.membrane_resultants(u)*self.wp[None, :, None]
w = self.wp
ke = np.einsum('p,epki,kl,eplj->eij', w, Bm, self.A, Bm,
optimize=True)
ke += np.einsum('p,pki,kl,plj->ij', w, Bb, self.D, Bb,
optimize=True)[None]
if np.any(self.B):
kb = np.einsum('p,epki,kl,plj->eij', w, Bm, self.B, Bb,
optimize=True)
ke += kb + kb.transpose(0, 2, 1)
Gx, Gy = self.Gx, self.Gy
ke += (np.einsum('ep,pi,pj->eij', Nr[..., 0], Gx, Gx)
+ np.einsum('ep,pi,pj->eij', Nr[..., 1], Gy, Gy)
+ np.einsum('ep,pi,pj->eij', Nr[..., 2], Gx, Gy)
+ np.einsum('ep,pi,pj->eij', Nr[..., 2], Gy, Gx))
return coo_matrix((ke.ravel(), (self._rows, self._cols)),
shape=(self.N, self.N)).tocsc()
def strain_energy(self, u, u0=None):
"""Strain energy of the strains ``eps(u + u0) - eps(u0)`` and the
curvatures ``kappa(u)``, independently of the kernels of the
elements"""
if u0 is None:
u0 = np.zeros_like(u)
eL1, eN1, k1 = self.generalized_strains(u + u0)
eL0, eN0, k0 = self.generalized_strains(u0)
e = (eL1 + eN1) - (eL0 + eN0)
k = k1 - k0
dens = (np.einsum('epi,ij,epj->ep', e, self.A, e)
+ 2*np.einsum('epi,ij,epj->ep', e, self.B, k)
+ np.einsum('epi,ij,epj->ep', k, self.D, k))
return 0.5*np.sum(dens*self.wp[None, :])
def w_max(self, u):
"""Largest nodal ``|w|``"""
return np.abs(u[6::DOF]).max()
def symmetrized_KF(model, variant):
r"""Load stiffness matrix symmetrized by eq. (63), (64) or (65) of the
paper, ``variant`` ``63``, ``64`` or ``65``, all dofs
With the skew part ``S = (KF - KF.T)/2``, which couples the in-plane
dofs (``u, ..., v_y``, block ``a``) with those of `w` (block ``b``),
``S = [[0, Kh], [-Kh.T, 0]]``, eq. (62a), the substitute ``[[0, Kh],
[Kh.T, 0]]``, eq. (62b), is added (64) or subtracted (65) to the
symmetric part ``(KF + KF.T)/2`` (63).
"""
KF = csc_matrix(model.KF)
Ksym = (KF + KF.T)/2
if variant == 63:
return Ksym.tocsc()
S = ((KF - KF.T)/2).tocoo()
inplane = (np.arange(model.N) % DOF) < 6
ab = inplane[S.row] & ~inplane[S.col]
Kh = coo_matrix((S.data[ab], (S.row[ab], S.col[ab])),
shape=S.shape).tocsc()
Ksub = Kh + Kh.T
sign = {64: 1., 65: -1.}[variant]
return (Ksym + sign*Ksub).tocsc()
def linear_buckling(model, load='follower', num_eigvalues=4, KF=None):
r"""Buckling factors ``k`` and modes, all dofs, of the linear buckling
problem ``(KC0 + k (KG(u_lin) + KF)) phi = 0``
``u_lin`` is the linear static solution of the dead load of ``k = 1``;
``load='dead'`` leaves out ``KF``, ``'follower'`` adds the unsymmetric
load stiffness matrix of the model or the matrix ``KF`` given, e.g. of
:func:`symmetrized_KF`. Returns the real parts of the load multipliers
and of the modes, and ``u_lin``.
"""
KC0uu = model.restrict(model.KC0)
clin = solve(KC0uu, model.fp0[model.bu], silent=True)
ulin = model.expand(clin)
KG = model.calc_KG(ulin)
symmetric = None
if load == 'follower':
KF = model.KF if KF is None else KF
KG = KG + KF
symmetric = None if KF is not model.KF else False
eigvals, eigvecs = lb(KC0uu, model.restrict(KG),
num_eigvalues=num_eigvalues, silent=True,
symmetric=symmetric)
modes = np.zeros((model.N, eigvecs.shape[1]))
modes[model.bu] = np.real(eigvecs)
return np.real(eigvals), modes, ulin
def num_waves(model, mode):
"""Number of circumferential waves of ``mode``, the dominant harmonic
of `w` along the free edge: ``cos(k y/R)`` and ``sin(k y/R)`` of any
``k`` for the closed cylinder, ``cos(j n y/R)`` for the sector"""
top = np.flatnonzero(isclose(model.x, model.L))
top = top[np.argsort(model.y[top])]
w = mode[DOF*top + 6]
y = model.y[top]
if model.full:
w, y = np.append(w, w[0]), np.append(y, 2*pi*R)
ns = np.arange(0, model.ny//2)
amp = [np.hypot(np.trapezoid(w*np.cos(k*y/R), y),
np.trapezoid(w*np.sin(k*y/R), y)) for k in ns]
else:
ns = model.n*np.arange(0, 4)
amp = [abs(np.trapezoid(w*np.cos(k*y/R), y)) for k in ns]
return int(ns[np.argmax(amp)])
def imperfection(model, mode, xi):
"""The dofs of `w` of ``mode``, scaled to a largest nodal ``|w|`` of
``xi*t``, outwards at the free edge"""
u0 = np.zeros_like(mode)
for dof in (6, 7, 8, 9):
u0[dof::DOF] = mode[dof::DOF]
w = u0[6::DOF]
return u0*xi*model.t/w[np.argmax(np.abs(w))]
class FollowerCylinder(object):
r"""Callables of structsolve.Analysis for the cantilever with the
initial imperfection ``u0`` under the pressure of the load factor ``k``
The callables work on the unknown dofs ``c = u[bu]`` and receive the load
factor ``inc`` from structsolve:
- ``calc_fext(inc) = inc*fp(u0)``, the force vector of the pressure on the
imperfect shell, and ``calc_fext(inc, c) = inc*fp(u + u0)``, the
load vector of the current configuration, ``-dR/dk``, which
structsolve uses in the predictors and in the constraint of the
arc-length methods when ``calc_fext`` declares ``c``;
- ``calc_fint(c, inc) = fint(u + u0) - fint(u0) - KG(N0) u - inc*(fp(u +
u0) - fp(u0))``, such that the residual of structsolve is that of the
follower pressure, ``fint_imp(u) - inc*fp(u + u0)``;
- ``calc_kC(c, NLgeom=True, inc) = KC0 + KCNL(u + u0) - KGNL0 + inc*KF``,
and ``calc_kC()`` the linear stiffness matrix of the imperfect shell,
``KC0 + KCNL(u0) - KGNL0 + KG(u0)``, without ``KF``;
- ``calc_kG(c) = KG(u)``,
see :class:`ImperfectCylinder` of ``simo1990_cylinder.py`` for the terms
of the imperfection. ``load='dead'`` leaves out the terms of ``fp(u)``
and ``KF``. With ``vectorized=True``, the default, ``KC0 + KCNL(u +
u0)`` is computed as ``KT(u + u0) - KG(u + u0)`` with
:meth:`CantileverModel.calc_KT`, the same matrix in a fraction of the
time of ``update_KCNL``; ``False`` calls the kernel.
"""
def __init__(self, model, u0=None, load='follower', vectorized=True):
self.model = m = model
self.vectorized = vectorized
self.load = load
self.u0 = np.zeros(m.N) if u0 is None else np.asarray(u0, float)
self.fint0 = m.calc_fint(self.u0)
N0NL = (m.membrane_resultants(self.u0, nonlinear=True)
- m.membrane_resultants(self.u0, nonlinear=False))
self.KGNL0 = m.calc_kG_N(N0NL)
self.KG0 = m.calc_KG(self.u0)
self.follower = load == 'follower'
f = m.calc_fp(self.u0) if self.follower else m.fp0
self.fextuu = f[m.bu]
self.KFuu = m.restrict(m.KF)
self.K0uu = m.restrict(m.KC0 + m.calc_KCNL(self.u0) - self.KGNL0
+ self.KG0)
def calc_fext(self, inc=1., c=None, silent=True):
if c is None or not self.follower:
return inc*self.fextuu
# the load vector of the current configuration, -dR/dk, for the
# predictors and the constraint of the arc-length methods
m = self.model
return inc*m.calc_fp(m.expand(c) + self.u0)[m.bu]
def calc_fint(self, c, inc=1., silent=True):
m = self.model
u = m.expand(c)
f = (m.calc_fint(u + self.u0) - self.fint0 - self.KG0 @ u
- self.KGNL0 @ u)
if self.follower:
# fp(u + u0) - fp(u0), affine in u for quadratic=0
f -= inc*m.calc_fp(u, reference=0)
return f[m.bu]
def calc_kC(self, c=None, NLgeom=False, inc=1., silent=True):
if c is None or not NLgeom:
# the linear stiffness matrix of the imperfect shell, without
# load stiffness: structsolve uses it for the linear reference
# solution, the convergence weights and the scaling of the
# count of negative eigenvalues
return self.K0uu
else:
m = self.model
u = m.expand(c)
if self.vectorized:
# KG(u + u0) = KG(u) + KG0, KG linear in u
K = m.restrict(m.calc_KT(u + self.u0) - m.calc_KG(u)
- self.KG0 - self.KGNL0)
else:
K = m.restrict(m.KC0 + m.calc_KCNL(u + self.u0)
- self.KGNL0)
if self.follower:
K = K + inc*self.KFuu
return K
def calc_kG(self, c=None, NLgeom=False, silent=True):
m = self.model
return m.restrict(m.calc_KG(m.expand(c)))
def run_riks(shell=1, bc='B1', nx=16, ny=4, xi=0.01, load='follower',
kinematics='sanders', wmax=3., inertia_num_eigvals=6,
silent=True, **kwargs):
r"""Riks analysis of the imperfect cantilever until the largest ``|w|``
reaches ``wmax*t``
The imperfection has the amplitude ``xi*t`` in the shape of the lowest
linear buckling mode of the same ``load``. ``kwargs`` are attributes of
:class:`structsolve.Analysis`. Returns a dict with the load factors
``ks``, the largest ``|w|/t`` ``wt``, from the unloaded state, the
current stiffness parameter ``cs_param`` of each increment, ``dk/|du|``
over its value of the first increment, with the Euclidean norm of the
increment of the unknown dofs (after Bergan, ref. [7] of the paper),
the counts of negative eigenvalues ``neg`` and the determinant signs
``det``, of the converged increments, the linear buckling factor
``k_lin``, the run time ``elapsed``, the analysis ``an``, the model and
the imperfection ``u0``.
"""
model = CantileverModel(shell, bc, nx, ny, kinematics)
eigvals, modes, _ = linear_buckling(model, load)
u0 = imperfection(model, modes[:, 0], xi)
cb = FollowerCylinder(model, u0, load)
an = Analysis(cb.calc_fext, cb.calc_fint, cb.calc_kC, cb.calc_kG)
an.NL_method = 'arc_length_riks'
an.final_load_factor = None
# first increment 0.5 of k; the arc length then grows with the number
# of iterations, the limit of maxInc being lifted
an.initialInc = 0.5
an.maxInc = 100.
an.minInc = 1e-4
an.maxArcLength = 1e9
an.inertia = True
an.inertia_num_eigvals = inertia_num_eigvals
an.linear_solver = 'splu'
an.stop_criterion = lambda an: (
model.w_max(model.expand(an.cs[-1])) >= wmax*model.t)
for k, v in kwargs.items():
setattr(an, k, v)
t0 = time.perf_counter()
an.static(NLgeom=True, silent=silent)
elapsed = time.perf_counter() - t0
wt = np.array([0.] + [model.w_max(model.expand(c)) for c in an.cs])
ks = np.concatenate(([0.], an.increments))
# current stiffness parameter: the stiffness dk/|du| of each increment,
# with the Euclidean norm of the increment of the unknown dofs, over
# that of the first one
cs = np.vstack((np.zeros(model.bu.sum()), np.array(an.cs)))
stiff = np.diff(ks)/np.linalg.norm(np.diff(cs, axis=0), axis=1)
return dict(ks=ks, wt=wt/model.t, cs_param=stiff/stiff[0],
neg=list(an.num_negative_eigvals), det=list(an.det_signs),
k_lin=eigvals[0], elapsed=elapsed, an=an, model=model, u0=u0)