Follower pressure#
Residual and tangent stiffness matrix#
The residual \(\boldsymbol{F}_{int}(\boldsymbol{u}) - \boldsymbol{F}_p(\boldsymbol{u})\) and its Jacobian
\(\boldsymbol{K}_{C0} + \boldsymbol{K}_{CNL} + \boldsymbol{K}_{GNL} + \boldsymbol{K}_G + \boldsymbol{K}_F\) of one element, for both
kinematics, the modules of Donnell-type and Sanders-type kinematics being selected by the
dictionary KINEMATICS, with a pressure that varies over the element:
def make_element(kinematics, nint=5):
r"""One element of the cylinder of ``test_tangent_consistency.py``,
with a pressure that varies over the element"""
mod, cls = KINEMATICS[kinematics]
R = 0.2032
L = 0.3556
prop = laminated_plate(stack=[45, -45, 0, 90, 90, 0, -45, 45],
laminaprop=(127.629e9, 11.3074e9, 0.300235,
6.00257e9, 6.00257e9, 6.00257e9),
plyt=0.00101539/8)
points, weights = get_points_weights(nint=nint)
elem = getattr(mod, cls)(nint)
elem.c1, elem.c2, elem.c3, elem.c4 = 0, DOF, 2*DOF, 3*DOF
elem.R = R
elem.lex = L/16
elem.ley = 2*pi*R/40
assign_constant_ABD(elem, prop)
xi, eta = np.meshgrid(points, points, indexing='ij')
# p[i, j] at (points[i], points[j]), external pressure
p = np.ascontiguousarray(-2.e5*(1 + 0.3*xi - 0.2*eta**2 + 0.1*xi*eta))
return mod, elem, points, weights, p, prop
def residual_and_tangent(kinematics, quadratic):
mod, elem, points, weights, p, prop = make_element(kinematics)
# a pressure large enough for the load stiffness to be a sizable part
# of the tangent, which a realistic pressure is not compared with the
# membrane stiffness
p = 1e3*p
elem.init_k_KC0 = elem.init_k_KCNL = elem.init_k_KG = 0
def assemble(update, size, u=None):
r = np.zeros(size, dtype=mod.INT)
c = np.zeros(size, dtype=mod.INT)
v = np.zeros(size, dtype=mod.DOUBLE)
if u is None:
update(elem, points, weights, r, c, v)
else:
update(u, elem, points, weights, r, c, v)
return coo_matrix((v, (r, c)), shape=(N, N)).toarray()
KC0 = assemble(mod.update_KC0, mod.KC0_SPARSE_SIZE)
def residual(u):
f = np.zeros(N, dtype=mod.DOUBLE)
mod.update_fint(u, elem, points, weights, f)
return f - fext_follower(mod, elem, points, weights, p, u, quadratic)
def tangent(u):
return (KC0 + assemble(mod.update_KCNL, mod.KCNL_SPARSE_SIZE, u)
+ assemble(mod.update_KG, mod.KG_SPARSE_SIZE, u)
+ KCfollower(mod, elem, points, weights, p, u, quadratic))
def follower(u):
return KCfollower(mod, elem, points, weights, p, u, quadratic)
return residual, tangent, follower, prop
The tangent stiffness matrix against a Taylor test of the residual:
@pytest.mark.parametrize('kinematics', sorted(KINEMATICS))
@pytest.mark.parametrize('quadratic', TRUNCATIONS)
@pytest.mark.parametrize('seed', [0, 1])
def test_tangent_is_jacobian_of_residual(kinematics, quadratic, seed):
residual, tangent, follower, prop = residual_and_tangent(kinematics,
quadratic)
u = random_state(prop, seed)
d = random_state(prop, seed + 10)
KTd = tangent(u) @ d
# the load stiffness carries weight in the tangent
assert np.linalg.norm(follower(u) @ d) > 1e-3*np.linalg.norm(KTd)
r0 = residual(u)
steps = [1e-2, 1e-3, 1e-4, 1e-5]
errors = [np.linalg.norm(residual(u + h*d) - r0 - h*KTd)
/np.linalg.norm(h*KTd) for h in steps]
for prev, err in zip(errors[:-1], errors[1:]):
assert err < 0.2*prev, errors
coefficients = [err/h for h, err in zip(steps, errors)]
assert max(coefficients)/min(coefficients) < 2., errors
Linear buckling of a ring#
A long cylinder under a uniform external pressure buckles as a ring. The
ring is modelled with one element along \(x\), the mesh closed around the
circumference, and symmetry planes that keep the modes with an even number
of circumferential waves; see the docstring of
tests/test_follower_pressure_validation.py. The matrices of the model,
with \(\boldsymbol{K}_G\) of the linear static solution of the dead pressure load and
the load stiffness matrix \(\boldsymbol{K}_{F_{ref}}\) of the same pressure at
\(\boldsymbol{u} = \boldsymbol{0}\):
def cylinder_model(kinematics, L, R, prop, nx, ny, nint=4):
r"""Cylinder of length ``L`` with ``nx`` nodes along `x` and ``ny``
elements around the closed circumference"""
mod, cls = KINEMATICS[kinematics]
points, weights = get_points_weights(nint=nint)
ids = np.arange(nx*ny).reshape(nx, ny)
x = np.repeat(np.linspace(0, L, nx), ny)
y = np.tile(np.arange(ny)*2*pi*R/ny, nx)
elements = []
for ix in range(nx - 1):
for iy in range(ny):
elem = getattr(mod, cls)(nint)
n1, n2 = ids[ix, iy], ids[ix + 1, iy]
n3, n4 = ids[ix + 1, (iy + 1) % ny], ids[ix, (iy + 1) % ny]
elem.n1, elem.n2, elem.n3, elem.n4 = n1, n2, n3, n4
elem.c1, elem.c2, elem.c3, elem.c4 = (DOF*n1, DOF*n2, DOF*n3,
DOF*n4)
elem.R = R
elem.lex = L/(nx - 1)
elem.ley = 2*pi*R/ny
assign_constant_ABD(elem, prop)
elements.append(elem)
return mod, elements, points, weights, x, y, DOF*nx*ny
def assemble(mod, elements, update, tag, size, *args):
n = len(elements)
spsize = getattr(mod, tag + '_SPARSE_SIZE')
r = np.zeros(spsize*n, dtype=mod.INT)
c = np.zeros(spsize*n, dtype=mod.INT)
v = np.zeros(spsize*n, dtype=mod.DOUBLE)
for i, elem in enumerate(elements):
setattr(elem, 'init_k_' + tag, i*spsize)
update(*args[:-1], elem, *args[-1], r, c, v)
return coo_matrix((v, (r, c)), shape=(size, size)).tocsc()
def matrices(mod, elements, points, weights, size, p, u0=None,
Nxx=None, Nyy=None):
r"""``KC0``, ``KG`` (of ``u0`` or of the constant ``Nxx``, ``Nyy``),
``KCfollower`` and the follower force vector at zero displacements, the
dead load"""
pp = np.full((len(points), len(points)), p)
KC0 = assemble(mod, elements, mod.update_KC0, 'KC0', size,
(points, weights))
zero = np.zeros(size)
KCf = assemble(mod, elements, mod.update_KCfollower, 'KCfollower', size,
zero, (points, weights, pp))
f = np.zeros(size)
for elem in elements:
mod.update_fext_follower(zero, elem, points, weights, pp, f)
if u0 is not None:
KG = assemble(mod, elements, mod.update_KG, 'KG', size, u0,
(points, weights))
else:
n = len(elements)
spsize = mod.KG_SPARSE_SIZE
r = np.zeros(spsize*n, dtype=mod.INT)
c = np.zeros(spsize*n, dtype=mod.INT)
v = np.zeros(spsize*n, dtype=mod.DOUBLE)
for i, elem in enumerate(elements):
elem.init_k_KG = i*spsize
mod.update_KG_constant_stress(elem, points, weights, r, c, v,
Nxx, Nyy, 0.)
KG = coo_matrix((v, (r, c)), shape=(size, size)).tocsc()
return KC0, KG, KCf, f
def ring(kinematics, R, h, ny):
prop = isotropic_plate(thickness=h, E=E, nu=NU)
out = cylinder_model(kinematics, 0.1*R, R, prop, 2, ny)
mod, elements, points, weights, x, y, size = out
bk = np.zeros(size, dtype=bool)
for d in (0, 1, 2, 4, 7, 9):
bk[d::DOF] = True
sym = np.zeros(size//DOF, dtype=bool)
for q in range(4):
sym |= np.isclose(y, q*pi*R/2)
bk[3::DOF] |= sym
bk[8::DOF] |= sym
return mod, elements, points, weights, size, ~bk, prop
The linear buckling problem \((\boldsymbol{K}_{C0} + \lambda (\boldsymbol{K}_G +
\boldsymbol{K}_{F_{ref}}))\boldsymbol{\phi} = \boldsymbol{0}\) of the follower pressure, or without
\(\boldsymbol{K}_{F_{ref}}\) for the dead pressure, solved with structsolve.lb,
with the unsymmetric solvers for the follower pressure:
def buckling(KC0, KG, KCf, bu, load):
K = KG + KCf if load == 'follower' else KG
Kuu = KC0[bu, :][:, bu]
KGuu = K[bu, :][:, bu]
eigvals, eigvecs = lb(Kuu, KGuu, silent=True, num_eigvalues=4,
symmetric=False if load == 'follower' else None)
return np.real(eigvals[0]), eigvecs[:, 0]
def ring_buckling_pressure(kinematics, load, R=1., h=0.01, ny=192):
r"""Load multiplier of a unit external pressure on the ring"""
mod, elements, points, weights, size, bu, prop = ring(kinematics, R, h,
ny)
p = -1.
KC0, _, KCf, f = matrices(mod, elements, points, weights, size, p,
u0=np.zeros(size))
u0 = np.zeros(size)
u0[bu] = spsolve(KC0[bu, :][:, bu], f[bu])
_, KG, _, _ = matrices(mod, elements, points, weights, size, p, u0=u0)
lam, mode = buckling(KC0, KG, KCf, bu, load)
return lam, prop
The buckling pressure against the closed-form solution of the ring with the same strains and load terms, and against its inextensional limit, for both kinematics and both pressures:
@pytest.mark.parametrize('kinematics', sorted(KINEMATICS))
@pytest.mark.parametrize('load', ['follower', 'dead'])
def test_isotropic_ring(kinematics, load):
R = 1.
lam, prop = ring_buckling_pressure(kinematics, load, R=R)
A22 = prop.A[1, 1]
D22 = prop.D[1, 1]
ref = -ring_closed_form(kinematics, load, 2, A22, D22, R)
# discretization error for n = 2, relative: 7e-3, 1.5e-4, 2.9e-6 and
# 7.9e-8 (Sanders), 3.9e-3, 8.3e-5, 1.7e-6 and 4.5e-8 (Donnell) with 24,
# 48, 96 and 192 elements around the circumference, see the module
# __main__; the convergence is that of the inextensional mode, which
# the cubic v and w only approximate
assert np.isclose(lam, ref, rtol=2e-7), (lam, ref)
inextensional = {('sanders', 'follower'): 3., ('sanders', 'dead'): 4.,
('donnell', 'follower'): 16/5, ('donnell', 'dead'): 4.}
# R/h = 100: the extensibility is of relative order (h/R)**2
assert np.isclose(lam*R**3/D22, inextensional[(kinematics, load)],
rtol=1e-3)