Non-linear static analysis with the Newton-Raphson method#
The geometrically non-linear models, based on Donnell’s equations, provide the
internal force vector \(\{F_{int}\}\), calculated with Shell.calc_fint(),
and the tangent stiffness matrix
calculated with Shell.calc_kT(), or with Shell.calc_kC() and
Shell.calc_kG() using NLgeom=True. The tangent stiffness matrix is
the exact derivative of the internal force vector, such that the
Newton-Raphson iterations converge quadratically.
The example below applies a compression of about 37 times the first linear buckling load in a single step, deep in the postbuckling regime, starting from the linear solution, and verifies the order of convergence of the iterations:
def test_nonlinear():
m = 6
n = 6
for model in [
'plate_clpt_donnell',
'cylshell_clpt_donnell',
'cylshell_clpt_sanders',
'plate_fsdt_donnell',
'plate_tsdt_donnell',
]:
print('Testing model: %s' % model)
s = Shell()
s.model = model
s.x1u = 0
s.x1ur = 1
s.x2u = 1
s.x2ur = 1
s.x1v = 1
s.x1vr = 1
s.x2v = 1
s.x2vr = 1
s.x1w = 0
s.x1wr = 1
s.x2w = 0
s.x2wr = 1
s.y1u = 1
s.y1ur = 1
s.y2u = 1
s.y2ur = 1
s.y1v = 0
s.y1vr = 1
s.y2v = 1
s.y2vr = 1
s.y1w = 0
s.y1wr = 1
s.y2w = 0
s.y2wr = 1
s.a = 4.
s.b = 1.
s.r = 1.e15
s.stack = [90, 0, 90, 0]
s.plyt = 1e-3*0.125
E11 = 142.5e9
E22 = E11/20
G12 = G13 = G23 = 0.5*E22
s.laminaprop = (E11, E22, 0.25, G12, G12, G12)
s.m = m
s.n = n
load = 700
Nxx = load/s.b
# distributed axial load
s.add_distr_load_fixed_x(s.a, funcx=lambda y: -Nxx, funcy=None, funcz=None, cte=False)
# perturbation load
s.add_point_load(s.a/2., s.b/2., 0, 0, 0.001, cte=True)
#initial
fext = s.calc_fext()
c = solve(s.calc_kC(), fext, silent=True)
plot_mesh, fields = s.uvw(c=c)
print(' linear wmax', fields['w'].max())
assert np.isclose(fields['w'].max(), 0.0026619, rtol=0.01)
# solving using the full Newton-Raphson method, with the tangent
# stiffness matrix updated at every iteration
D = s.calc_kC().diagonal() # at beginning of load increment
epsilon = 1.e-10
errors = []
while True:
fint = s.calc_fint(c=c)
Ri = fint - fext
crisfield_test = scaling(Ri, D)/max(scaling(fext, D), scaling(fint, D))
errors.append(crisfield_test)
print(' iteration %d, crisfield_test %1.3e' % (len(errors) - 1, crisfield_test))
if crisfield_test < epsilon:
break
if len(errors) > 30:
raise RuntimeError('Not converged!')
KT = s.calc_kT(c=c)
c = c + solve(KT, -Ri, silent=True)
plot_mesh, fields = s.uvw(c=c)
print(' nonlinear wmax', fields['w'].max())
assert np.isclose(fields['w'].max(), 0.004841, rtol=0.01)
# quadratic convergence: once in the asymptotic range, the order
# log(e_k+1)/log(e_k) approaches 2, whereas an inconsistent tangent
# gives an order of 1 and needs hundreds of iterations
assert len(errors) <= 15, 'too many iterations: %d' % len(errors)
orders = [np.log(e1)/np.log(e0) for e0, e1 in zip(errors[:-1], errors[1:])
if 1.e-15 < e1 and e0 < 1.e-2]
print(' convergence orders', orders)
if 'fsdt' in model or 'tsdt' in model:
#NOTE for this very thin plate, a/h = 8000, the residual of the
# shear deformation theories reaches its round-off floor, of
# the order of eps*(G/E)*(a/h)**2 ~ 1e-9, within the asymptotic
# range, where the order cannot be measured; their quadratic
# convergence is pinned by test_tangent_consistency.py
continue
assert len(orders) >= 2
assert min(orders) > 1.6, 'convergence is not quadratic: %s' % orders
The consistency between the tangent stiffness matrix and the internal force vector is verified by a directional Taylor test, where the error of the linear approximation of \(\{F_{int}\}\) must decrease quadratically with the step size:
@pytest.mark.parametrize('model', MODELS)
@pytest.mark.parametrize('seed', [0, 1, 2])
def test_tangent_is_derivative_of_fint(model, seed):
"""Taylor test: the error must fall by about ten for each decade of h."""
s, fint, KT = make_callables(model)
rng = np.random.default_rng(seed)
c = random_state(s, rng)
d = random_state(s, rng)
f0 = fint(c)
KTd = KT(c) @ d
scale = np.linalg.norm(KTd)
assert scale > 0
steps = [1.e-2, 1.e-3, 1.e-4, 1.e-5, 1.e-6]
errors = [np.linalg.norm(fint(c + h*d) - f0 - h*KTd)/(h*scale)
for h in steps]
# a consistent tangent leaves a second order remainder, so the error
# falls by ten for every decade of h; an inconsistent one leaves a first
# order remainder, so the error plateaus instead
for h, prev, err in zip(steps[1:], errors[:-1], errors[1:]):
assert err < 0.2*prev, (
'error did not fall by ten going to h=%.0e: %.3e -> %.3e; the '
'tangent is not the derivative of fint' % (h, prev, err))
# error/h is the constant that multiplies the second derivative, so it
# must not drift with h
coefficients = [err/h for h, err in zip(steps, errors)]
spread = max(coefficients)/min(coefficients)
assert spread < 2., (
'error/h drifted by a factor %.1f over %.0e..%.0e, so the remainder '
'is not second order' % (spread, steps[0], steps[-1]))
assert errors[-1] < 1.e-4, (
'residual %.3e at h=%.0e is too large for a consistent tangent'
% (errors[-1], steps[-1]))
The methods calc_fext, calc_fint, calc_kC and calc_kG of
Shell and MultiDomain have the signatures required by
structsolve.Analysis, which implements the Newton-Raphson method with
load control and the arc-length methods, see Non-linear static analysis with the arc-length methods.