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 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(), which needs lex, ley and ABD of each element, and the force vector \(\boldsymbol{F}_{ext}\) of Assembling a model. The tests solve it with scipy.sparse.linalg.spsolve, 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 \(\boldsymbol{F}_{ext}\), such that the membrane forces are those of the plate with its boundary conditions. Since \(\boldsymbol{K}_G\) is linear in \(\boldsymbol{u}_0\), \(\lambda\) multiplies \(\boldsymbol{F}_{ext}\);
update_KG_cte_N(), from membrane forces Nxx, Nyy and Nxy given for each element, positive in tension, without a static analysis: \(\lambda\) multiplies these membrane forces, e.g. Nxx = -1 gives the critical compressive force per unit length.

The tests compute the lowest load multipliers with scipy.sparse.linalg.eigsh, writing the problem as \(\boldsymbol{K}_G\boldsymbol{\phi} = \mu \boldsymbol{K}_{C0}\boldsymbol{\phi}\), with \(\boldsymbol{K}_{C0}\) positive definite once the boundary conditions are applied, and \(\lambda = -1/\mu\):

eigvals, eigvecsu = eigsh(A=KGuu, k=num_eigenvalues, which='SM',
                          M=Kuu, tol=1e-6, sigma=1., mode='cayley')
eigvals = -1./eigvals

see Linear buckling analysis.

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 no integration points, and \(\omega_n\) the natural angular frequencies, in radians per unit time. The test of Natural frequencies solves the problem with the dense scipy.linalg.eigh and takes \(\omega_n\) as the square roots of its eigenvalues.

Nonlinear static analysis#

For large deflections, with the nonlinear strains of von Karman, 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 with Taylor tests and finite differences. The internal forces always include the nonlinear terms, and for small displacements \(\boldsymbol{F}_{int} \approx \boldsymbol{K}_{C0} \boldsymbol{u}\).

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, 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 bfscplate2d, see Installing bfscplate2d. The linear buckling and natural frequency problems above are solved by structsolve.lb(Kuu, KGuu) and structsolve.freq(Kuu, Muu), which return the load multipliers \(\lambda\) and \(-\omega_n^2\), as cross-checked in tests/test_linear_buckling_iso_SSSS.py and tests/test_natural_frequency.py.

The nonlinear static analysis is solved by structsolve.Analysis, with four callables that loop over the elements, call the functions of the element 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 bfscplate2d require, see Installing bfscplate2d. 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 mesh ``m`` under the
    follower pressure ``func(x, y)`` 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.
    """
    def __init__(self, m, bu, func, quadratic):
        self.m = m
        self.bu = bu
        self.quadratic = quadratic
        self.pressures = [m.pressure(plate, func) for plate in m.elements]
        self.KC0 = m.KC0()
        N = m.N
        num_elements = len(m.elements)
        self.KCNLr = np.zeros(KCNL_SPARSE_SIZE*num_elements, dtype=INT)
        self.KCNLc = np.zeros(KCNL_SPARSE_SIZE*num_elements, dtype=INT)
        self.KCNLv = np.zeros(KCNL_SPARSE_SIZE*num_elements, dtype=DOUBLE)
        self.KGr = np.zeros(KG_SPARSE_SIZE*num_elements, dtype=INT)
        self.KGc = np.zeros(KG_SPARSE_SIZE*num_elements, dtype=INT)
        self.KGv = np.zeros(KG_SPARSE_SIZE*num_elements, dtype=DOUBLE)
        self.KCfr = np.zeros(KCfollower_SPARSE_SIZE*num_elements, dtype=INT)
        self.KCfc = np.zeros(KCfollower_SPARSE_SIZE*num_elements, dtype=INT)
        self.KCfv = np.zeros(KCfollower_SPARSE_SIZE*num_elements,
                             dtype=DOUBLE)
        # the dead part of the pressure, the force vector at u = 0
        self.fp0 = np.zeros(N, dtype=DOUBLE)
        for plate, p in zip(m.elements, self.pressures):
            update_fext_follower(np.zeros(N), plate, m.points, m.weights, p,
                                 self.fp0, quadratic)

    def expand(self, c):
        """The displacements of the model from those of the unknown dofs."""
        u = np.zeros(self.m.N, dtype=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):
        m = self.m
        u = self.expand(c)
        fint = np.zeros(m.N, dtype=DOUBLE)
        for plate, p in zip(m.elements, self.pressures):
            update_fint(u, plate, m.points, m.weights, fint)
            # minus the part of the follower pressure that depends on u
            update_fext_follower(u, plate, m.points, m.weights, -inc*p, fint,
                                 self.quadratic, 0)
        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]
        m = self.m
        u = self.expand(c)
        self.KCNLv *= 0
        self.KCfv *= 0
        for plate, p in zip(m.elements, self.pressures):
            update_KCNL(u, plate, m.points, m.weights, self.KCNLr,
                        self.KCNLc, self.KCNLv)
            # the load stiffness matrix of the pressure inc*p
            update_KCfollower(u, plate, m.points, m.weights, inc*p,
                              self.KCfr, self.KCfc, self.KCfv,
                              self.quadratic)
        KCNL = coo_matrix((self.KCNLv, (self.KCNLr, self.KCNLc)),
                          shape=(m.N, m.N))
        KCf = coo_matrix((self.KCfv, (self.KCfr, self.KCfc)),
                         shape=(m.N, m.N))
        return (self.KC0 + KCNL + KCf).tocsc()[bu, :][:, bu]

    def calc_kG(self, c=None, NLgeom=False, silent=False):
        m = self.m
        u = self.expand(c)
        self.KGv *= 0
        for plate in m.elements:
            update_KG(u, plate, m.points, m.weights, self.KGr, self.KGc,
                      self.KGv)
        KG = coo_matrix((self.KGv, (self.KGr, self.KGc)), shape=(m.N, m.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.