Follower pressure#
The tests of the follower pressure use a mesh class with the assembly of
the force vector of update_fext_follower(), of the load stiffness
matrix of update_KCfollower(), and of the internal forces and the
tangent stiffness matrix. The method pressure gives the array p of
each element at its integration points:
class Mesh(object):
"""Rectangular mesh of ``nx`` by ``ny`` nodes over ``a`` by ``b``"""
def __init__(self, nx=4, ny=3, a=A, b=B, lam=LAMINATE, nint=NINT):
self.a, self.b = a, b
self.points, self.weights = get_points_weights(nint=nint)
xlin = np.linspace(0, a, nx)
ylin = np.linspace(0, b, ny)
xmesh, ymesh = np.meshgrid(xlin, ylin)
ncoords = np.vstack((xmesh.T.flatten(), ymesh.T.flatten())).T
self.x = ncoords[:, 0]
self.y = ncoords[:, 1]
nids_mesh = np.arange(nx*ny).reshape(nx, ny)
n1s = nids_mesh[:-1, :-1].flatten()
n2s = nids_mesh[1:, :-1].flatten()
n3s = nids_mesh[1:, 1:].flatten()
n4s = nids_mesh[:-1, 1:].flatten()
self.N = DOF*nx*ny
self.elements = []
for k, (n1, n2, n3, n4) in enumerate(zip(n1s, n2s, n3s, n4s)):
plate = BFSCPlate2D()
plate.n1, plate.n2, plate.n3, plate.n4 = n1, n2, n3, n4
plate.c1 = DOF*n1
plate.c2 = DOF*n2
plate.c3 = DOF*n3
plate.c4 = DOF*n4
plate.ABD = lam.ABD
plate.lex = a/(nx - 1)
plate.ley = b/(ny - 1)
plate.init_k_KC0 = k*KC0_SPARSE_SIZE
plate.init_k_KCNL = k*KCNL_SPARSE_SIZE
plate.init_k_KG = k*KG_SPARSE_SIZE
plate.init_k_KCfollower = k*KCfollower_SPARSE_SIZE
self.elements.append(plate)
self.h = lam.h
def origin(self, plate):
return self.x[plate.n1], self.y[plate.n1]
def pressure(self, plate, func):
r"""``p[i, j]`` at ``xi = points[i]``, ``eta = points[j]``"""
x1, y1 = self.origin(plate)
x = x1 + plate.lex*(self.points + 1)/2
y = y1 + plate.ley*(self.points + 1)/2
p = func(x[:, None], y[None, :])*np.ones((x.size, y.size))
return np.ascontiguousarray(p, dtype=DOUBLE)
def fext_follower(self, u, func, quadratic, reference=1):
f = np.zeros(self.N, dtype=DOUBLE)
for plate in self.elements:
update_fext_follower(u, plate, self.points, self.weights,
self.pressure(plate, func), f, quadratic,
reference)
return f
def _matrix(self, size, call):
n = size*len(self.elements)
r = np.zeros(n, dtype=INT)
c = np.zeros(n, dtype=INT)
v = np.zeros(n, dtype=DOUBLE)
for plate in self.elements:
call(plate, r, c, v)
return coo_matrix((v, (r, c)), shape=(self.N, self.N)).tocsc()
def KCfollower(self, u, func, quadratic):
return self._matrix(KCfollower_SPARSE_SIZE,
lambda plate, r, c, v: update_KCfollower(u, plate, self.points,
self.weights, self.pressure(plate, func), r, c, v,
quadratic))
def KC0(self):
return self._matrix(KC0_SPARSE_SIZE,
lambda plate, r, c, v: update_KC0(plate, self.points,
self.weights, r, c, v))
def KT(self, u):
KCNL = self._matrix(KCNL_SPARSE_SIZE,
lambda plate, r, c, v: update_KCNL(u, plate, self.points,
self.weights, r, c, v))
KG = self._matrix(KG_SPARSE_SIZE,
lambda plate, r, c, v: update_KG(u, plate, self.points,
self.weights, r, c, v))
return self.KC0() + KCNL + KG
def fint(self, u):
f = np.zeros(self.N, dtype=DOUBLE)
for plate in self.elements:
update_fint(u, plate, self.points, self.weights, f)
return f
def random_state(self, rng, scale=1.):
"""In-plane displacements of the order of the thickness, deflections
well above it, such that the nonlinear terms carry weight"""
u = scale*self.h*rng.standard_normal(self.N)
u[6::DOF] *= 20
lmin = min(self.a, self.b)
for i in (1, 2, 4, 5, 7, 8):
u[i::DOF] *= 1/lmin
u[7::DOF] *= 20
u[8::DOF] *= 20
u[9::DOF] *= 20/lmin**2
return u
def positions(self, plate):
return np.concatenate([np.arange(c, c + DOF) for c in
(plate.c1, plate.c2, plate.c3, plate.c4)])
Large deflection of a cantilever plate#
The Newton-Raphson method with the residual \(\boldsymbol{F}_{int}(\boldsymbol{u}) - \boldsymbol{F}_p(\boldsymbol{u})\)
of the pressure of the load factor lbd and the tangent stiffness matrix
\(\boldsymbol{K}_{C0} + \boldsymbol{K}_{CNL} + \boldsymbol{K}_{GNL} + \boldsymbol{K}_G + \boldsymbol{K}_F\), or with the dead pressure
\(\boldsymbol{F}_p(\boldsymbol{0})\) for comparison:
def newton_raphson(m, bu, func, load, increments=5, tol=1.e-9,
max_iterations=30):
r"""Solves `\boldsymbol{F}_{int}(\boldsymbol{u}) = \lambda
\boldsymbol{F}_p(\boldsymbol{u})` (``load`` 0 or 1, the truncation of the
follower pressure) or `\boldsymbol{F}_{int}(\boldsymbol{u}) = \lambda
\boldsymbol{F}_p(\boldsymbol{0})` (``load =
'dead'``) for `\lambda` up to 1, returning the solution and the relative
residuals of each increment"""
u = np.zeros(m.N)
histories = []
for lbd in np.linspace(0, 1, increments + 1)[1:]:
pf = lambda x, y: lbd*func(x, y)
history = []
for count in range(max_iterations):
if load == 'dead':
F = m.fext_follower(0*u, pf, 0)
else:
F = m.fext_follower(u, pf, load)
R = m.fint(u) - F
history.append(np.linalg.norm(R[bu])/np.linalg.norm(F[bu]))
if history[-1] < tol:
break
KT = m.KT(u)
if load != 'dead':
KT = KT + m.KCfollower(u, pf, load)
u[bu] -= spsolve(KT[bu, :][:, bu].tocsc(), R[bu])
else:
raise RuntimeError('not converged: %s' % history)
histories.append(history)
return u, histories
A square isotropic plate clamped at \(x = 0\), the other edges free, under a uniform pressure. The load stiffness matrix is unsymmetric because of the free edges, the follower pressure gives a tip deflection measurably larger than the dead pressure, and the two truncations of the area vector give the same result:
def boundary_flags(m, edges):
r"""Constrained degrees of freedom: for each edge in ``edges``, a dict
``{'x1': 'w', ...}`` of the displacements set to zero along it, with
their derivatives along the edge"""
bk = np.zeros(m.N, dtype=bool)
on = {'x1': isclose(m.x, 0), 'x2': isclose(m.x, m.a),
'y1': isclose(m.y, 0), 'y2': isclose(m.y, m.b)}
# value and derivative along the edge: y for the edges x = const
along = {'x': {'u': (0, 2), 'v': (3, 5), 'w': (6, 8)},
'y': {'u': (0, 1), 'v': (3, 4), 'w': (6, 7)}}
for edge, fields in edges.items():
for f in fields:
for i in along[edge[0]][f]:
bk[i::DOF] |= on[edge]
return bk
def test_cantilever_large_deflection_against_panels():
r"""Square isotropic cantilever, clamped at `x = 0`, the other edges
free, under a uniform pressure pushing along `+z`, with the tip
deflection about 20 times the thickness, `w/a = 0.13`. The follower
pressure tilts with the plate and its in-plane component `-p w_{,x}`
pulls the tip towards the clamp: the tip deflection is 0.92% larger and
the tip displacement `u` 1.9% larger in magnitude than under the dead
pressure, against panels within 5e-4 for `w`, i.e. the follower effect
itself within 1%. The load stiffness is unsymmetric (free edges) and
Newton-Raphson converges quadratically with it"""
E, nu, h = 70.e9, 0.3, 2.e-3
a = b = 0.3
nel = 10
m = Mesh(nx=nel + 1, ny=nel + 1, a=a, b=b,
lam=isotropic_plate(thickness=h, E=E, nu=nu))
bk = boundary_flags(m, {'x1': 'uvw'})
bk[7::DOF] |= isclose(m.x, 0)
bk[9::DOF] |= isclose(m.x, 0)
bu = ~bk
func = lambda x, y: 2.e3 + 0*x*y
K = m.KCfollower(np.zeros(m.N), func, 0)
assert asym(K[bu, :][:, bu]) > 1e-3
tip = np.where(isclose(m.x, a) & isclose(m.y, b/2))[0][0]
out = {}
for load in ('dead', 0, 1):
u, histories = newton_raphson(m, bu, func, load)
out[load] = dict(w=u[DOF*tip + 6]/h, u=u[DOF*tip])
# quadratic, e_k+1 <= 10 e_k**1.8, from 1e-3 down to the round-off of
# the residual, about 1e-11 here
checked = 0
for history in histories:
for e0, e1 in zip(history[:-1], history[1:]):
if 1.e-10 < e0 < 1.e-3:
assert e1 <= max(10*e0**1.8, 1.e-10), history
checked += 1
assert checked >= 3, histories
for key in ('w', 'u'):
assert isclose(out[load][key], PANELS[load][key], rtol=1.e-3), (
load, key, out[load][key], PANELS[load][key])
effect = out[0]['w']/out['dead']['w'] - 1
effect_ref = PANELS[0]['w']/PANELS['dead']['w'] - 1
assert effect > 9.e-3
assert isclose(effect, effect_ref, rtol=1.e-2), (effect, effect_ref)
assert out[0]['u']/out['dead']['u'] - 1 > 1.5e-2
# the truncation of the area vector is immaterial, as follower_pressure.md
# of panels predicts for the moderate rotations of the von Karman strains
assert isclose(out[1]['w'], out[0]['w'], rtol=1.e-6)
Tangent stiffness matrix#
\(\boldsymbol{K}_{C0} + \boldsymbol{K}_{CNL} + \boldsymbol{K}_{GNL} + \boldsymbol{K}_G + \boldsymbol{K}_F\) against the finite differences of the residual, for both truncations:
@pytest.mark.parametrize('quadratic', TRUNCATIONS)
def test_tangent_is_jacobian_of_residual(quadratic):
r"""`\boldsymbol{K}_T = \boldsymbol{K}_{C0} + \boldsymbol{K}_{CNL} + \boldsymbol{K}_{GNL} +
\boldsymbol{K}_G + \boldsymbol{K}_F` against the central differences of
`\boldsymbol{R} = \boldsymbol{F}_{int} - \boldsymbol{F}_p`, every entry,
and a Taylor test"""
m = Mesh(nx=3, ny=2)
func = lambda x, y: 100*p_function(x, y)
rng = np.random.default_rng(7)
u = m.random_state(rng)
def R(u):
return m.fint(u) - m.fext_follower(u, func, quadratic)
KT = (m.KT(u) + m.KCfollower(u, func, quadratic)).toarray()
# steps relative to each degree of freedom, large enough to stay clear
# of the round-off of fint, whose stiffness is ~500 times that of the
# load here
step = 1.e-4
J = np.empty((m.N, m.N))
scale = np.abs(u).copy()
scale[scale == 0] = 1.
for j in range(m.N):
e = np.zeros(m.N)
e[j] = step*scale[j]
J[:, j] = (R(u + e) - R(u - e))/(2*e[j])
err = np.linalg.norm(KT - J)/np.linalg.norm(J)
assert err < 1.e-8, err
# also small compared to the load stiffness, which is ~0.2% of KT
Kf = m.KCfollower(u, func, quadratic).toarray()
assert np.linalg.norm(KT - J)/np.linalg.norm(Kf) < 1.e-5
d = m.random_state(rng)
R0 = R(u)
KTd = KT @ d
steps = [1.e-2, 1.e-3, 1.e-4, 1.e-5]
errors = [np.linalg.norm(R(u + h*d) - R0 - h*KTd)/np.linalg.norm(h*KTd)
for h in steps]
for h, prev, e in zip(steps[1:], errors[:-1], errors[1:]):
assert e < 0.2*prev, errors
Symmetry of the load stiffness matrix#
The conditions of Follower pressure under which the load stiffness matrix is symmetric:
def test_symmetry():
r"""Symmetric when `p` is uniform over the whole plate and, on each edge,
`w = 0` or the displacement normal to the edge is zero, unsymmetric
otherwise"""
m = Mesh(nx=5, ny=4)
u0 = np.zeros(m.N)
uniform = lambda x, y: -1.e4 + 0*x*y
ss = {'x1': 'w', 'x2': 'w', 'y1': 'w', 'y2': 'w'}
normal = {'x1': 'u', 'x2': 'u', 'y1': 'v', 'y2': 'v'}
mixed = {'x1': 'w', 'x2': 'u', 'y1': 'v', 'y2': 'w'}
K = m.KCfollower(u0, uniform, 0)
# the element matrix alone is unsymmetric, its edges are interior edges
single = Mesh(nx=2, ny=2)
assert asym(single.KCfollower(np.zeros(single.N), uniform, 0)) > 1e-3
for edges in (ss, normal, mixed):
bu = ~boundary_flags(m, edges)
assert asym(K[bu, :][:, bu]) < 1e-12, edges
# a free edge
free = dict(ss)
free.pop('x2')
bu = ~boundary_flags(m, free)
assert asym(K[bu, :][:, bu]) > 1e-3
# a varying pressure
bu = ~boundary_flags(m, ss)
Kf = m.KCfollower(u0, p_function, 0)
assert asym(Kf[bu, :][:, bu]) > 1e-3
# a pressure on part of the plate
patch = lambda x, y: np.where(x < m.a/2, -1.e4, 0.) + 0*y
Kp = m.KCfollower(u0, patch, 0)
assert asym(Kp[bu, :][:, bu]) > 1e-3
# the complete area vector away from u = 0: with w = 0 on the whole
# boundary and a uniform p the load has the potential -p V, V = int w
# ((1 + u,x) (1 + v,y) - u,y v,x) dA, hence a symmetric matrix also at
# u != 0, unsymmetric with a free edge
rng = np.random.default_rng(2)
u = m.random_state(rng)
u[boundary_flags(m, ss)] = 0
Kq = m.KCfollower(u, uniform, 1)
K0 = m.KCfollower(0*u, uniform, 1)
assert np.linalg.norm((Kq - K0).toarray()) > 1e-2*np.linalg.norm(K0.toarray())
assert asym(Kq[bu, :][:, bu]) < 1e-12
bu = ~boundary_flags(m, free)
assert asym(Kq[bu, :][:, bu]) > 1e-3