Analyses#

The matrices and vectors of Assembling a model, restricted to the unknown degrees of freedom bu, give the systems of the analyses below, solved with SciPy as in the tests, or with structsolve, see Solving with structsolve. The functions are those of both kinematics, see Donnell-type and Sanders-type kinematics, and the theory of each matrix is in the Notes of its function.

Linear static analysis#

\[\boldsymbol{K}_{C0}\boldsymbol{u} = \boldsymbol{F}_{ext}\]

with \(\boldsymbol{K}_{C0}\) of update_KC0() and the force vector or the prescribed displacements of Assembling a model. The tests solve it with scipy.sparse.linalg.spsolve or with the conjugate gradient method of scipy.sparse.linalg.cg, with the matrix and the right-hand side scaled by a scalar, see Linear static analysis.

Linear buckling analysis#

\[(\boldsymbol{K}_{C0} + \lambda \boldsymbol{K}_G)\boldsymbol{\phi} = \boldsymbol{0}\]

The critical loads are \(\lambda\) times the load that gives \(\boldsymbol{K}_G\), and \(\boldsymbol{\phi}\) are the buckling modes. \(\boldsymbol{K}_G\) is computed in one of two ways:

update_KG(), from the displacements \(\boldsymbol{u}_0\) of the linear static solution of a reference load or prescribed displacement, such that the membrane stress resultants are those of the shell with its boundary conditions. Since \(\boldsymbol{K}_G\) is linear in \(\boldsymbol{u}_0\), \(\lambda\) multiplies the reference load, e.g. the sum of the reaction forces of a prescribed shortening;
update_KG_constant_stress(), from the membrane stress resultants Nxx, Nyy and Nxy given for each element, the last arguments of the function, without a static analysis: \(\lambda\) multiplies these resultants, e.g. Nxx = -1 for an axial compression, of total force \(-2 \pi R\) along the circumference.

The tests compute the lowest load multipliers with scipy.sparse.linalg.lobpcg, for the problem \(\boldsymbol{K}_{C0}\boldsymbol{\phi} = \lambda (-\boldsymbol{K}_G)\boldsymbol{\phi}\), with an incomplete LU factorization of \(\boldsymbol{K}_{C0}\), spilu, as preconditioner, see Linear buckling analysis. The matrices can also be given to the linear buckling solver of structsolve, structsolve.lb, as in the validation tests of the follower pressure, see Follower pressure.

Natural frequencies#

\[(\boldsymbol{K}_{C0} - \omega_n^2 \boldsymbol{M})\boldsymbol{\phi} = \boldsymbol{0}\]

with the mass matrix \(\boldsymbol{M}\) of update_M(), which needs lex, ley, h and rho of each element, and R for the Sanders-type kinematics, and no integration points, and \(\omega_n\) the natural angular frequencies, in radians per unit time. The lowest frequencies of a cylinder are clustered; the tests use the function lowest_frequencies of tests/natural_frequency_solver.py, which solves the shift-inverted problem with ARPACK from a fixed starting vector and verifies the residual of each eigenpair and that no eigenvalue is missing, see Natural frequencies.

Nonlinear static analysis#

For large displacements, with the nonlinear strains of the kinematics, the equilibrium of the internal forces of update_fint() with the external forces is solved for the displacements, e.g. with the Newton-Raphson method:

\[\boldsymbol{R}(\boldsymbol{u}) = \boldsymbol{F}_{int}(\boldsymbol{u}) - \lambda \boldsymbol{F}_{ext}, \qquad \boldsymbol{K}_T(\boldsymbol{u}) = \boldsymbol{K}_{C0} + \boldsymbol{K}_{CNL}(\boldsymbol{u}) + \boldsymbol{K}_{GNL}(\boldsymbol{u}) + \boldsymbol{K}_G(\boldsymbol{u}), \qquad \boldsymbol{K}_T \Delta \boldsymbol{u} = -\boldsymbol{R}\]

with the load factor \(\lambda\) increased in steps and the tangent stiffness matrix \(\boldsymbol{K}_T\) of update_KC0(), update_KCNL() and update_KG(), all computed with the current displacements and the same integration points as \(\boldsymbol{F}_{int}\). \(\boldsymbol{K}_T\) is the exact Jacobian of \(\boldsymbol{F}_{int}\), such that the Newton-Raphson iterations converge quadratically; tests/test_tangent_consistency.py checks it for both kinematics with Taylor tests and finite differences. \(\boldsymbol{K}_G\) keeps only the part linear in the displacements, the one of the linear buckling analysis, and the rest of the geometric terms is \(\boldsymbol{K}_{GNL}\), computed with \(\boldsymbol{K}_{CNL}\) by update_KCNL(), see its Notes.

The value arrays of \(\boldsymbol{K}_{CNL} + \boldsymbol{K}_{GNL}\) and \(\boldsymbol{K}_G\) are zeroed before each assembly, since the functions add to them, see Assembling a model. The test of Nonlinear static analysis solves the iterations with structsolve, with the convergence criterion of the scaled norms of the residual and of the forces, see Solving with structsolve and Nonlinear static analysis with structsolve. A follower pressure adds its force vector to the residual and its load stiffness matrix to \(\boldsymbol{K}_T\), see Follower pressure.

Solving with structsolve#

The solvers of structsolve receive the matrices and vectors of the model, assembled by the user. It is used by the tests and the examples, and it is not a requirement of bfsccylinder, see Installing bfsccylinder. The linear buckling problem above is solved by structsolve.lb, as in the validation tests of the follower pressure.

The nonlinear static analysis is solved by structsolve.Analysis, with four callables that loop over the elements, call the functions of the element, of either kinematics, and assemble, see Nonlinear static analysis with structsolve:

calc_fext(inc=1., silent=False), the external forces \(\boldsymbol{F}_{ext}\) times the load factor inc;
calc_fint(c, silent=False), the internal forces of update_fint() at the displacements c;
calc_kC(c=None, NLgeom=False, silent=False), \(\boldsymbol{K}_{C0}\) with NLgeom=False, and \(\boldsymbol{K}_{C0} + \boldsymbol{K}_{CNL}(\boldsymbol{c}) + \boldsymbol{K}_{GNL}(\boldsymbol{c})\) with NLgeom=True;
calc_kG(c=None, NLgeom=False, silent=False), \(\boldsymbol{K}_G(\boldsymbol{c})\) of update_KG(), called with NLgeom=True;

such that calc_kC + calc_kG is the tangent stiffness matrix \(\boldsymbol{K}_T\) and the residual of structsolve is the one above, \(\boldsymbol{R} = \boldsymbol{F}_{int}(\boldsymbol{c}) - \lambda \boldsymbol{F}_{ext}\), with the Newton-Raphson corrections \(\boldsymbol{K}_T \Delta \boldsymbol{c} = -\boldsymbol{R}\), the convention of structsolve 0.9.3 or newer, which the examples of bfsccylinder require, see Installing bfsccylinder. The value arrays are zeroed in each call. The load factor goes from 0 to 1 in increments, the first one of size Analysis.initialInc and the next ones growing up to Analysis.maxInc, each starting from a predictor extrapolated from the previous increment, the linear solution of \(\boldsymbol{K}_{C0}\) in the first one; Analysis.static returns the load factors and the solutions of the converged increments.

structsolve imposes no boundary condition itself. The constrained degrees of freedom are removed in one of two ways:

the partition of the free degrees of freedom, as in the tests: the callables receive c = u[bu] and return the vectors restricted to bu and the matrices to bu by bu, and the solutions are expanded with u[bu] = c;
null rows and columns: the callables return the vectors and matrices of the whole model, with the rows and columns of the constrained degrees of freedom zeroed in all the matrices. The solvers remove them and return vectors with the size of the model, zero at these degrees of freedom.

The convergence criterion is, by default, Crisfield’s: the residual and the reference forces, the external forces and the internal forces, are scaled with the diagonal of the linear stiffness matrix calc_kC() over the free degrees of freedom, and the iterations stop when their ratio is below Analysis.relTOL = 1e-9. The scaled norm does not depend on the units of the degrees of freedom, which here mix displacements, their first derivatives and \(w_{,xy}\), and leaves out the reactions at the constrained ones; Analysis.convergence = 'relative' selects the Euclidean norm instead. With a consistent tangent the criterion is reached in a few iterations.

A follower pressure \(\lambda p\), see Follower pressure, enters through the optional keyword argument inc, the load factor, which structsolve passes to the callables that declare it. The external forces are the dead part \(\boldsymbol{F}_p(\boldsymbol{0})\) of p, calc_fint subtracts the part that depends on the displacements, reference=0, and calc_kC adds the load stiffness matrix of \(\lambda p\), as in tests/test_follower_pressure_structsolve.py, see Nonlinear static analysis with structsolve:

class FollowerCallables(object):
    r"""The callables of structsolve.Analysis for the ``elements`` of the
    module ``mod`` under the follower pressures ``pressures``, one array of
    shape ``(nint, nint)`` per element, times the load factor

    They loop over the elements with the update functions and assemble, and
    work on the unknown degrees of freedom ``bu`` 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``.
    ``quadratic`` is the truncation of the area vector. ``calc_fint`` also
    keeps the load factor and the internal forces of each iteration, for the
    residual history.
    """
    def __init__(self, mod, elements, points, weights, N, bu, pressures,
                 quadratic=0):
        self.mod = mod
        self.elements = elements
        self.points = points
        self.weights = weights
        self.N = N
        self.bu = bu
        self.pressures = pressures
        self.quadratic = quadratic
        self.iterations = []
        self.KC0 = assemble(mod, elements, mod.update_KC0, 'KC0', N,
                            (points, weights))
        num_elements = len(elements)
        for i, shell in enumerate(elements):
            shell.init_k_KCNL = i*mod.KCNL_SPARSE_SIZE
            shell.init_k_KG = i*mod.KG_SPARSE_SIZE
            shell.init_k_KCfollower = i*mod.KCfollower_SPARSE_SIZE
        size = mod.KCNL_SPARSE_SIZE*num_elements
        self.KCNLr = np.zeros(size, dtype=mod.INT)
        self.KCNLc = np.zeros(size, dtype=mod.INT)
        self.KCNLv = np.zeros(size, dtype=mod.DOUBLE)
        size = mod.KG_SPARSE_SIZE*num_elements
        self.KGr = np.zeros(size, dtype=mod.INT)
        self.KGc = np.zeros(size, dtype=mod.INT)
        self.KGv = np.zeros(size, dtype=mod.DOUBLE)
        size = mod.KCfollower_SPARSE_SIZE*num_elements
        self.KCfr = np.zeros(size, dtype=mod.INT)
        self.KCfc = np.zeros(size, dtype=mod.INT)
        self.KCfv = np.zeros(size, dtype=mod.DOUBLE)
        # the dead part of the pressure, the force vector at u = 0
        self.fp0 = np.zeros(N, dtype=mod.DOUBLE)
        for shell, p in zip(elements, pressures):
            mod.update_fext_follower(np.zeros(N), shell, points, weights, p,
                                     self.fp0, quadratic)

    def expand(self, c):
        """The displacements of the model from those of the unknown dofs."""
        u = np.zeros(self.N, dtype=self.mod.DOUBLE)
        u[self.bu] = c
        return u

    def calc_fext(self, inc=1., silent=False):
        return inc*self.fp0[self.bu]

    def calc_fint(self, c, inc=1., silent=False):
        mod = self.mod
        u = self.expand(c)
        fint = np.zeros(self.N, dtype=mod.DOUBLE)
        for shell, p in zip(self.elements, self.pressures):
            mod.update_fint(u, shell, self.points, self.weights, fint)
            # minus the part of the follower pressure that depends on u
            mod.update_fext_follower(u, shell, self.points, self.weights,
                                     -inc*p, fint, self.quadratic, 0)
        self.iterations.append((inc, fint[self.bu]))
        return fint[self.bu]

    def calc_kC(self, c=None, NLgeom=False, inc=1., silent=False):
        bu = self.bu
        if not NLgeom:
            return self.KC0[bu, :][:, bu]
        mod = self.mod
        u = self.expand(c)
        self.KCNLv *= 0
        self.KCfv *= 0
        for shell, p in zip(self.elements, self.pressures):
            mod.update_KCNL(u, shell, self.points, self.weights, self.KCNLr,
                            self.KCNLc, self.KCNLv)
            # the load stiffness matrix of the pressure inc*p
            mod.update_KCfollower(u, shell, self.points, self.weights,
                                  inc*p, self.KCfr, self.KCfc, self.KCfv,
                                  self.quadratic)
        KCNL = coo_matrix((self.KCNLv, (self.KCNLr, self.KCNLc)),
                          shape=(self.N, self.N))
        KCf = coo_matrix((self.KCfv, (self.KCfr, self.KCfc)),
                         shape=(self.N, self.N))
        return (self.KC0 + KCNL + KCf).tocsc()[bu, :][:, bu]

    def calc_kG(self, c=None, NLgeom=False, silent=False):
        mod = self.mod
        u = self.expand(c)
        self.KGv *= 0
        for shell in self.elements:
            mod.update_KG(u, shell, self.points, self.weights, self.KGr,
                          self.KGc, self.KGv)
        KG = coo_matrix((self.KGv, (self.KGr, self.KGc)),
                        shape=(self.N, self.N))
        return KG.tocsc()[self.bu, :][:, self.bu]

The tangent stiffness matrix is then unsymmetric in general, which the solvers of structsolve support with LU factorizations. For the arc-length methods, Analysis.NL_method, calc_fext also declares the configuration c and returns the follower force vector \(\boldsymbol{F}_p(\boldsymbol{c})\) of p, minus the derivative of the residual \(\boldsymbol{R} = \boldsymbol{F}_{int} - \lambda \boldsymbol{F}_p(\boldsymbol{c})\) with respect to the load factor.