Prescribed displacements imposed exactly#
A prescribed displacement can be imposed exactly as a linear constraint on the degrees of freedom, \(\boldsymbol{c} = \boldsymbol{Z} \boldsymbol{c}_r + \boldsymbol{c}_{p0} + \lambda \boldsymbol{c}_{p1}\), where \(\boldsymbol{Z}\) is a basis of the null space of the constraints and \(\boldsymbol{c}_{p0}\), \(\boldsymbol{c}_{p1}\) are particular solutions for the constant and the incremented prescribed values, as done by the null-space method of panels. The reduced residual is
and minus its derivative with respect to the load factor,
depends on the configuration and on the load factor. calc_fext receives
both when it declares the keyword arguments c and lbd, see
callbacks. Under a prescribed displacement alone, \(F_{ext} = 0\), the
reference vector calc_fext(inc=1.) and the reduced internal force vector
both vanish at equilibrium, and the solvers use \(-\partial \boldsymbol{R}/\partial
\lambda\) for the displacement scale of the arc-length methods, the predictor
of the first Newton-Raphson increment and the convergence
criterion.
Three hardening springs in series, with the displacement of the last node prescribed:
class SpringsProblem(object):
def __init__(self, k=(100., 150., 200.), a=(20., 40., 60.), d0=0.,
d1=1., F=0., with_lbd=True):
self.k = np.asarray(k)
self.a = np.asarray(a)
self.cp0 = np.array([0., 0., d0])
self.cp1 = np.array([0., 0., d1])
self.Fext = np.array([F, 0., 0.])
self.with_lbd = with_lbd
self.an = None
# (lbd, c_r) of each call to calc_fext with lbd, and (inc, ||R_r||,
# number of converged increments) of each call to calc_fint
self.fext_calls = []
self.fint_calls = []
def full(self, c_r, lbd):
return Z @ c_r + self.cp0 + lbd*self.cp1
def elongations(self, c):
return np.array([c[0], c[1] - c[0], c[2] - c[1]])
def fint_full(self, c):
k, a = self.k, self.a
f = k*self.elongations(c) + a*self.elongations(c)**3
return np.array([f[0] - f[1], f[1] - f[2], f[2]])
def kT_full(self, c):
ke = self.k + 3*self.a*self.elongations(c)**2
return np.array([[ke[0] + ke[1], -ke[1], 0.],
[-ke[1], ke[1] + ke[2], -ke[2]],
[0., -ke[2], ke[2]]])
def residual(self, c_r, lbd):
return Z.T @ (self.fint_full(self.full(c_r, lbd)) - lbd*self.Fext)
def load_vector(self, c_r, lbd):
"""-dR_r/dlbd"""
c = self.full(c_r, lbd)
return Z.T @ (self.Fext - self.kT_full(c) @ self.cp1)
def reaction(self, c_r, lbd):
r"""Force of the support on node 3, `\boldsymbol{c}_{p1}^\top
\boldsymbol{r} / d_1`, with `\boldsymbol{r} =
\boldsymbol{F}_{int}(\boldsymbol{c}) - \lambda \boldsymbol{F}_{ext}`"""
c = self.full(c_r, lbd)
r = self.fint_full(c) - lbd*self.Fext
return self.cp1 @ r/self.cp1[2]
# callables of structsolve.Analysis
def calc_fext(self, inc=1., silent=True, c=None, lbd=None):
if c is None:
return inc*(Z.T @ self.Fext)
self.fext_calls.append((lbd, c.copy()))
return inc*self.load_vector(c, lbd)
def calc_fint(self, c, inc=1., silent=True):
step = len(self.an.increments) if self.an is not None else 0
self.fint_calls.append((inc, np.linalg.norm(self.residual(c, inc)),
step))
return Z.T @ self.fint_full(self.full(c, inc))
def calc_kC(self, c=None, NLgeom=False, inc=1., silent=True):
if c is None or not NLgeom:
return csr_matrix(Z.T @ self.kT_full(np.zeros(3)) @ Z)
return csr_matrix(Z.T @ self.kT_full(self.full(c, inc)) @ Z)
def calc_kG(self, c=None, NLgeom=False, silent=True):
return csr_matrix((2, 2))
def analysis(self, **kwargs):
if self.with_lbd:
calc_fext = self.calc_fext
else:
calc_fext = lambda inc=1., silent=True: self.calc_fext(inc)
an = Analysis(calc_fext, self.calc_fint, self.calc_kC, self.calc_kG)
for key, value in kwargs.items():
setattr(an, key, value)
self.an = an
return an
def exact(self, lbd):
r"""Solution without the load on node 1: the three springs carry the
same force `P`, found from `\sum_i e_i(P) = x_3`"""
assert self.Fext[0] == 0.
x3 = self.cp0[2] + lbd*self.cp1[2]
def e(P, i):
return brentq(lambda ei: self.k[i]*ei + self.a[i]*ei**3 - P,
-1e4, 1e4, xtol=1e-15, rtol=1e-15)
def gap(P):
return sum(e(P, i) for i in range(3)) - x3
P = 0. if x3 == 0 else brentq(gap, -1e4, 1e4, xtol=1e-14,
rtol=1e-15)
e1, e2 = e(P, 0), e(P, 1)
return np.array([e1, e1 + e2]), P
def solve_reference(self, lbd, tol=1e-14):
"""Dense Newton-Raphson with continuation on the load factor"""
c_r = np.zeros(2)
for l in np.linspace(0, lbd, 21)[1:]:
for _ in range(100):
R = self.residual(c_r, l)
if np.linalg.norm(R) <= tol*max(1., np.linalg.norm(
self.load_vector(c_r, l))):
break
kT = Z.T @ self.kT_full(self.full(c_r, l)) @ Z
c_r = c_r + np.linalg.solve(kT, -R)
return c_r
@pytest.mark.parametrize('method', ['NR', 'arc_length_riks',
'arc_length_crisfield'])
@pytest.mark.parametrize('convergence', ['crisfield', 'relative'])
def test_displacement_control_relative_criterion(method, convergence):
r"""The prescribed displacement alone, with the tolerance of the relative
force imbalance only, ``absTOL=None``: every converged increment is the
exact solution, and ``calc_fext`` receives the load factor of the
configuration"""
pb = SpringsProblem(d1=0.8)
an = pb.analysis(NL_method=method, initialInc=0.2, relTOL=1e-10,
absTOL=None, convergence=convergence)
an.static(NLgeom=True, silent=True)
assert np.isclose(an.increments[-1], 1.)
assert len(an.increments) > 2
for lbd, c_r in zip(an.increments, an.cs):
ref, P = pb.exact(lbd)
assert np.allclose(c_r, ref, rtol=1e-8, atol=1e-12), (lbd, c_r, ref)
assert np.isclose(pb.reaction(c_r, lbd), P, rtol=1e-8)
assert (np.linalg.norm(pb.residual(c_r, lbd))
<= 1e-9*np.linalg.norm(lbd*pb.load_vector(c_r, lbd)))
# the arc-length solvers evaluate -dR/dlbd at the configurations and load
# factors of their iterations
if method != 'NR':
assert len(pb.fext_calls) > len(an.increments)
lbds = [lbd for lbd, _ in pb.fext_calls]
assert lbds[0] == 0.
assert all(lbd is not None for lbd in lbds)
assert len(set(lbds)) > 2
A constant prescribed displacement together with an incremented one and a load:
@pytest.mark.parametrize('method', ['NR', 'arc_length_riks',
'arc_length_crisfield'])
def test_constant_and_incremented_prescriptions(method):
r"""A constant prescribed displacement `d_0`, held during the analysis,
with a load `\lambda F` on node 1 and an incremented displacement
`\lambda d_1`"""
pb = SpringsProblem(d0=-0.3, d1=0.6, F=40.)
an = pb.analysis(NL_method=method, initialInc=0.2, relTOL=1e-11)
an.static(NLgeom=True, silent=True)
assert np.isclose(an.increments[-1], 1.)
for lbd, c_r in zip(an.increments, an.cs):
ref = pb.solve_reference(lbd)
assert np.allclose(c_r, ref, rtol=1e-8, atol=1e-12), (lbd, c_r, ref)
# at lbd = 0 the constant displacement alone is not in equilibrium with
# c_r = 0: the first increment starts from a non-null residual
assert np.linalg.norm(pb.residual(np.zeros(2), 0.)) > 1.
Both together: the null-space method and constraint_values#
A displacement prescribed by the null-space method, through calc_fext
with c and lbd, combines with the non-homogeneous constraints
\(\boldsymbol{C}^\top \boldsymbol{c} = \lambda \boldsymbol{d}\) of
Analysis.constraints and Analysis.constraint_values on the reduced
problem: the solutions of \(-\partial \boldsymbol{R}/\partial \lambda\) have
the border \(\boldsymbol{d}\), those of the residual a null border. The
springs above, with a load on the node 2 and either the node 1 or the
relative displacement of the nodes 1 and 2 prescribed as well, against
the problem with both prescriptions eliminated by hand, a single unknown
found by scipy.optimize.brentq(), for all the non-linear solvers:
def reference(p, case, lbd):
"""The configuration of the reduced problem at ``lbd``, see the module
docstring"""
x3 = lbd*D1
if case == 'node1':
x1 = lbd*E
def g(x2):
return (spring_force(p, 1, x2 - x1) - spring_force(p, 2, x3 - x2)
- lbd*F)
x2 = brentq(g, -5., 5., xtol=1e-15, rtol=1e-15)
return np.array([x1, x2])
def g(x1):
x2 = x1 + lbd*E
return spring_force(p, 0, x1) - spring_force(p, 2, x3 - x2) - lbd*F
x1 = brentq(g, -5., 5., xtol=1e-15, rtol=1e-15)
return np.array([x1, x1 + lbd*E])
@pytest.mark.parametrize('case', list(CASES))
@pytest.mark.parametrize('method', METHODS)
def test_path(method, case):
p, an = make(case, method)
an.static(NLgeom=True, silent=True)
assert an.increments[-1] == 1.
assert len(an.increments) > 3
# calc_fext received the load factor
assert p.fext_calls and any(lbd for lbd, _ in p.fext_calls)
C = CASES[case]
for lbd, c in zip(an.increments, an.cs):
assert abs(C @ c - lbd*E) <= 1e-12
np.testing.assert_allclose(c, reference(p, case, lbd), rtol=1e-8,
atol=1e-12)
# the residual of the free direction vanishes
Z = np.array([1., 1.]) if case == 'relative' else np.array([0., 1.])
R = p.residual(c, lbd)
assert abs(Z @ R) <= 1e-7*np.abs(p.load_vector(c, lbd)).max()
# far from the linear springs
lin = reference(SpringsProblem(a=(0., 0., 0.)), case, 1.)
assert np.abs(an.cs[-1] - lin).max() > 0.05