Nonlinear static analysis with structsolve#
The plate of Nonlinear static analysis, solved with structsolve.Analysis
of structsolve, see
Solving with structsolve. The four callables loop over
the elements with the functions of Nonlinear static analysis and work on
the unknown degrees of freedom bu, the partition of the free degrees of
freedom; calc_fint also keeps the internal forces of each iteration:
class Callables(object):
"""The callables of structsolve.Analysis for the model ``m`` under the
external forces ``fext``.
They loop over the elements with the update functions and assemble, and
work on the unknown degrees of freedom only, the partition of the free
dofs of structsolve: the solution vector of the solvers is ``c =
u[bu]``, and the force vectors and matrices are restricted to ``bu``.
``calc_fint`` also keeps the load factor and the internal forces of each
iteration, for the residual history.
"""
def __init__(self, m, fext):
self.m = m
self.fext = fext
self.iterations = []
def expand(self, c):
"""The displacements of the model from those of the unknown dofs."""
u = np.zeros(self.m.N)
u[self.m.bu] = c
return u
def calc_fext(self, inc=1., silent=False):
return inc*self.fext[self.m.bu]
def calc_fint(self, c, inc=1., silent=False):
fint = calc_fint(self.m, self.expand(c))[self.m.bu]
self.iterations.append((inc, fint))
return fint
def calc_kC(self, c=None, NLgeom=False, silent=False):
if not NLgeom:
return self.m.KC0uu
bu = self.m.bu
KCNL = calc_KCNL(self.m, self.expand(c))
return (self.m.KC0 + KCNL)[bu, :][:, bu]
def calc_kG(self, c=None, NLgeom=False, silent=False):
bu = self.m.bu
return calc_KG(self.m, self.expand(c))[bu, :][:, bu]
The analysis applies the load in equal increments, with
Analysis.initialInc and Analysis.maxInc, checks that none was cut,
and computes Crisfield’s relative residual, the convergence criterion of
structsolve, at each iteration of the last increment:
def newton_raphson(m, fext, increments=1):
"""Full Newton-Raphson of structsolve.Analysis, with ``increments``
equal load increments. Returns the solution and the residual history of
the last increment.
The tangent is rebuilt every iteration, so the residual history is the
direct evidence of whether KT is the derivative of fint: a consistent
tangent squares the residual each step, an inconsistent one multiplies
it by a constant factor instead. The residual is measured with
Crisfield's criterion, the convergence criterion of structsolve, with
its default tolerance of 1e-9.
"""
callables = Callables(m, fext)
an = Analysis(calc_fext=callables.calc_fext,
calc_fint=callables.calc_fint,
calc_kC=callables.calc_kC,
calc_kG=callables.calc_kG)
an.initialInc = 1./increments
an.maxInc = 1./increments
increments_done, cs = an.static(NLgeom=True)
# every increment converged, none was cut
assert np.allclose(increments_done,
np.linspace(0, 1, increments + 1)[1:]), (
'Not converged! load factors %s' % increments_done)
u = callables.expand(cs[-1])
w = crisfield_weights(m.KC0uu)
f = fext[m.bu]
history = [crisfield_zeta(inc*f - fint, inc*f, fint, w)
for inc, fint in callables.iterations if inc == 1.]
return u, history
A consistent tangent stiffness matrix squares the residual at the last iteration:
def assert_quadratic_convergence(history, epsilon=1.e-9):
"""The last Newton step must square the residual, not merely shrink it.
Quadratic convergence is an asymptotic statement, so it is the final
transition that carries it: the step that takes the iterate from above
the tolerance to below it. Earlier steps are not a fair test here,
because the residual is normalised by the whole external load, which the
axial part may satisfy at the first step while the out-of-plane part is
still far from converged.
With an inconsistent tangent the residual falls by a roughly constant
factor throughout, so no step ever squares it and this fails.
"""
assert history[-1] < epsilon, 'did not reach the tolerance'
assert len(history) <= 6, (
'took %d iterations to converge; a consistent tangent needs a '
'handful' % len(history))
assert len(history) >= 2, 'no Newton step was taken'
prev, curr = history[-2], history[-1]
# the floor keeps the test meaningful once prev**2 drops below the
# round-off level of the residual itself, about 1e-14 here: the last
# increment of test_large_deflection goes 2.7e-8 -> 1.5e-14
assert curr < max(10.*prev**2, 1.e-13), (
'the last step took the residual %.3e -> %.3e, which is linear '
'rather than quadratic; the tangent is not the derivative of fint'
% (prev, curr))
With a small transverse force the solution in one increment reproduces the linearised one, \((\boldsymbol{K}_{C0} + \boldsymbol{K}_G)^{-1}\boldsymbol{F}_{ext}\), with \(\boldsymbol{K}_G\) of the axial load:
def test_static(plot=False):
"""Small transverse load: the linearised buckling amplification."""
load = 300. # N
transverse = 0.01 # N
m = build_model()
fext, ftrans = forces(m, load, transverse)
load_cr, u_axial = linear_buckling_load(m, load)
print('load_cr', load_cr)
assert isclose(load_cr, 501.59022126865693, rtol=1e-3)
assert load < load_cr
# the exact linearised answer: the transverse load resisted by the
# bending stiffness reduced by the geometric stiffness of the axial load
KG = calc_KG(m, u_axial)
u_lin = np.zeros(m.N)
u_lin[m.bu] = spsolve((m.KC0 + KG)[m.bu, :][:, m.bu], ftrans[m.bu])
w_lin = u_lin[6::DOF].max()
# and the transverse load on its own, for the amplification factor
u_bend = np.zeros(m.N)
u_bend[m.bu] = spsolve(m.KC0uu, ftrans[m.bu])
w_bend = u_bend[6::DOF].max()
u, history = newton_raphson(m, fext + ftrans)
assert_quadratic_convergence(history)
w = u[6::DOF].reshape(nx, ny).T
print('w min max', w.min(), w.max())
print('w linearised ', w_lin)
print('amplification ', w.max()/w_bend)
print('amplification ceiling ', 1./(1. - load/load_cr))
if plot:
import matplotlib
matplotlib.use('TkAgg')
import matplotlib.pyplot as plt
plt.gca().set_aspect('equal')
levels = np.linspace(w.min(), w.max(), 300)
plt.contourf(m.xmesh, m.ymesh, w, levels=levels)
plt.colorbar()
plt.show()
# w/h is about 0.011 here, so the membrane stretching is negligible and
# the nonlinear answer must reproduce the linearised one
assert isclose(w.max(), w_lin, rtol=1e-3)
# a point load excites many modes, so the amplification stays below the
# single-mode ceiling 1/(1 - P/Pcr); it must still exceed one
amplification = w.max()/w_bend
assert 1. < amplification < 1./(1. - load/load_cr)
assert isclose(w.max(), 1.1104539457685102e-05, rtol=1e-3)
Follower pressure#
The callables of a follower pressure, FollowerCallables, which receive
the load factor inc, are given in Solving with structsolve.
The cantilever of Follower pressure, with free edges and an
unsymmetric load stiffness matrix, solved with these callables in five
load increments, against the reference values of panels:
@pytest.mark.parametrize('quadratic', TRUNCATIONS)
def test_cantilever_large_deflection_against_panels(quadratic):
r"""The square isotropic cantilever of ``test_follower_pressure.py``,
clamped at `x = 0`, the other edges free, under a uniform follower
pressure pushing along `+z`, in five equal load increments, against
panels within 1e-3 for the tip deflection `w` and the tip displacement
`u`"""
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
tip = np.where(isclose(m.x, a) & isclose(m.y, b/2))[0][0]
callables = FollowerCallables(m, bu, func, quadratic)
an = Analysis(calc_fext=callables.calc_fext,
calc_fint=callables.calc_fint,
calc_kC=callables.calc_kC,
calc_kG=callables.calc_kG)
an.initialInc = 0.2
an.maxInc = 0.2
increments, cs = an.static(NLgeom=True)
# every increment converged, none was cut
assert np.allclose(increments, [0.2, 0.4, 0.6, 0.8, 1.]), increments
u = callables.expand(cs[-1])
w_tip = u[DOF*tip + 6]/h
u_tip = u[DOF*tip]
print('w/h', w_tip, 'u', u_tip)
assert isclose(w_tip, PANELS[quadratic]['w'], rtol=1.e-3)
assert isclose(u_tip, PANELS[quadratic]['u'], rtol=1.e-3)