Tangent stiffness matrix and its product by a vector

This page documents how the internal force vector and its exact Jacobian, the tangent stiffness matrix, are implemented:

• Tangent stiffness of the CLPT models with Donnell kinematics derives \(\{F_{int}\}\) and \([K_T]\) of the models plate_clpt_donnell and cylshell_clpt_donnell and relates every term to the function of the implementation that evaluates it, including the two integration paths, the grouping of the terms between Shell.calc_kC() and Shell.calc_kG(), and the multi-domain assembly;

• Product of the tangent stiffness matrix by a vector derives the product \([K_T(c)] \{v\}\) for all the geometrically non-linear models, Shell.calc_kTv(), evaluated without assembling \([K_T]\).

Tangent stiffness of the CLPT models with Donnell kinematics

This section documents the internal force vector and the tangent stiffness matrix as they are implemented in the models plate_clpt_donnell and cylshell_clpt_donnell, both based on the classical laminated plate theory (CLPT) with the Donnell-Mushtari-Vlasov (von Kármán) kinematics. Each equation is related to the function of the implementation that evaluates it, such that the reader can move between the equations and the code in both directions. The objective is to show that each term of the tangent stiffness matrix follows from the exact derivative of the internal force vector,

(1)\[\boldsymbol{K}_T(\bar{\boldsymbol{c}}) = \frac{\partial \boldsymbol{F}_{int}}{\partial \bar{\boldsymbol{c}}}\]

which is what delivers the quadratic convergence of the Newton-Raphson iterations \(\boldsymbol{K}_T(\bar{\boldsymbol{c}}_k) \Delta \bar{\boldsymbol{c}} = \boldsymbol{F}_{ext} - \boldsymbol{F}_{int}(\bar{\boldsymbol{c}}_k)\), see Non-linear static analysis with the Newton-Raphson method. In the notation of Product of the tangent stiffness matrix by a vector, \([B(\beta)] = \boldsymbol{B}_L + \boldsymbol{B}_{NL}\), \([G] = \boldsymbol{G}\) and \([F] = \boldsymbol{C}\).

The tangent stiffness matrix is classically decomposed as

(2)\[\boldsymbol{K}_T = \boldsymbol{K}_0 + \boldsymbol{K}_{0L} + \boldsymbol{K}_{L0} + \boldsymbol{K}_{LL} + \boldsymbol{K}_G(\boldsymbol{N})\]

with the constitutive matrices \(\boldsymbol{K}_0\) to \(\boldsymbol{K}_{LL}\) and the geometric stiffness \(\boldsymbol{K}_G\) of the membrane stress \(\boldsymbol{N}\), see (16) and (14). The implementation is organised around a three-way split:

(3)\[\boldsymbol{K}_T(\bar{\boldsymbol{c}}) = \boldsymbol{K}_{C0} + \boldsymbol{K}_{CNL}(\bar{\boldsymbol{c}}) + \boldsymbol{K}_G(\bar{\boldsymbol{c}})\]

whose relation to the grouping of the implementation, stated in a comment of fkC_num in plate_clpt_donnell_num.pyx and in its cylindrical twin, is established in Which stress enters the geometric stiffness, and the degree of homogeneity and tabulated in Symbols and implementation variables. The correspondence is not a renaming only: the boundary between the constitutive and the geometric contributions of (3) sits at a different place from the one of (2), and Where the naming and the code do not line up says exactly where.

The two implementation paths

Each model contributes two modules, registered in panels/modelDB.py under the keys matrices and matrices_num:

• The closed-form path, matrices, integrates the full domain analytically. The integrals of products of the approximation functions are evaluated in closed form, and the resulting matrices are constant with respect to \(\bar{\boldsymbol{c}}\). It provides fk0 and fkG0 of plate_clpt_donnell.pyx and cylshell_clpt_donnell.pyx.

• The numerically integrated path, matrices_num, applies a Gauss-Legendre rule over \(n_x \times n_y\) points, and is the only path that can honour a displacement state \(\bar{\boldsymbol{c}}\), a laminate stiffness varying over the domain, or a partial integration domain. It provides fkC_num, fkG_num and calc_fint of plate_clpt_donnell_num.pyx and cylshell_clpt_donnell_num.pyx.

The selection is made in Shell.calc_kC() and Shell.calc_kG() by two flags, analytical_kC and analytical_kG, initialised to True and cleared by the conditions of the following table.

Conditions that force the numerically integrated path

Condition	Affects
Shell.is_partial_domain() is True	kC and kG
c is given	kC and kG
the argument ABDnxny is given, or the attribute Shell.ABDnxny is set	kC
NLgeom=True	kC and kG
c_cte is given	the kG block inside kC

Shell.is_partial_domain() reports whether the integration limits x1, x2, y1, y2 differ from the edges of the shell, \(0\), \(a\) and \(b\), after Shell.integration_limits() has replaced a limit left as None by the corresponding edge and has validated \(0 \le x_1 < x_2 \le a\) and \(0 \le y_1 < y_2 \le b\). Two consequences are worth stating explicitly:

• Shell.calc_kT() calls both Shell.calc_kC() and Shell.calc_kG() with NLgeom=True. The tangent stiffness matrix therefore always takes the numerically integrated path, and the closed-form path never participates in a non-linear analysis.

• Conversely, the closed-form path is reached only for \(\boldsymbol{K}_{C0}\) and for \(\boldsymbol{K}_G\) evaluated at a uniform prescribed membrane stress state, i.e. for a linear static analysis and for the classical linear buckling problem \((\boldsymbol{K}_0 + \lambda \boldsymbol{K}_G) \boldsymbol{q} = \boldsymbol{0}\), over a full domain with a constant laminate.

Approximation functions

The three mid-surface displacements are approximated using the hierarchic polynomials \(P_i\) of Bardell, see Bardell’s Hierarchical Functions, evaluated by f, fp and fpp of panels/core/include/bardell_functions.hpp, implemented in panels/core/src/bardell_functions.cpp, with the Python wrappers calc_f, calc_fp and calc_fpp of panels/bardell.pyx. Independent boundary flags are used for each displacement component and for each direction:

(4)\[\begin{split}\begin{aligned} u(x,y) &= \sum_{i=1}^{m} \sum_{j=1}^{n} c_{ij}^u \, P_i^u(\xi) P_j^u(\eta) \\ v(x,y) &= \sum_{i=1}^{m} \sum_{j=1}^{n} c_{ij}^v \, P_i^v(\xi) P_j^v(\eta) \\ w(x,y) &= \sum_{i=1}^{m} \sum_{j=1}^{n} c_{ij}^w \, P_i^w(\xi) P_j^w(\eta) \end{aligned}\end{split}\]

where the superscript in \(P_i^u(\xi)\) indicates that the flags \(\delta_{t1}\), \(\delta_{r1}\), \(\delta_{t2}\), \(\delta_{r2}\) of the first four terms refer to the boundary conditions of \(u\) at \(\xi = -1\) and \(\xi = +1\), and similarly for \(P_j^u(\eta)\), \(P_i^v(\xi)\), and so forth. Only the first four terms of vec_f in bardell_functions.cpp carry these flags; the remaining terms vanish together with their first derivative at both ends, which is what makes the set hierarchic. The implementation provides 30 terms. The natural coordinates are:

\[\xi = \frac{2x}{a} - 1 \qquad \eta = \frac{2y}{b} - 1\]

where, for the cylindrical shell, \(y = r\theta\) is the circumferential arc length and \(b\) the corresponding arc length of the panel. Derivatives with respect to the physical coordinates follow from the chain rule:

\[\frac{\partial}{\partial x} = \frac{2}{a}\frac{\partial}{\partial \xi} \qquad \frac{\partial}{\partial y} = \frac{2}{b}\frac{\partial}{\partial \eta}\]

which is the origin of every factor (2/a) and (2/b) appearing in the kernels, for instance in the strains of calc_fint. In vector notation:

\[u = \boldsymbol{S}^u \bar{\boldsymbol{c}} \qquad v = \boldsymbol{S}^v \bar{\boldsymbol{c}} \qquad w = \boldsymbol{S}^w \bar{\boldsymbol{c}}\]

with the alternated ordering of the degrees of freedom:

\[\bar{\boldsymbol{c}} = \left\{ c_{11}^u \quad c_{11}^v \quad c_{11}^w \quad c_{21}^u \quad c_{21}^v \quad c_{21}^w \quad \cdots \quad c_{mn}^u \quad c_{mn}^v \quad c_{mn}^w \right\}^\top\]

such that, using zero-based indices \(i\) and \(j\), the coefficients \(c_{ij}^u\), \(c_{ij}^v\) and \(c_{ij}^w\) occupy the positions \(3(jm+i)\), \(3(jm+i)+1\) and \(3(jm+i)+2\) of \(\bar{\boldsymbol{c}}\). This is the mapping col = col0 + DOF*(j*m + i) with DOF = 3 of the kernels, and the size of \(\bar{\boldsymbol{c}}\) is \(3mn\), see Shell.get_size(). The only non-zero components of \(\boldsymbol{S}^u\), \(\boldsymbol{S}^v\) and \(\boldsymbol{S}^w\) are:

\[S^u_{3(jm+i)} = P_i^u(\xi) P_j^u(\eta) \qquad S^v_{3(jm+i)+1} = P_i^v(\xi) P_j^v(\eta) \qquad S^w_{3(jm+i)+2} = P_i^w(\xi) P_j^w(\eta)\]

and, for instance:

\[S^w_{,x} = \frac{2}{a} P_i^{w\prime}(\xi) P_j^w(\eta) \qquad S^w_{,yy} = \frac{4}{b^2} P_i^w(\xi) P_j^{w\prime\prime}(\eta) \qquad S^w_{,xy} = \frac{4}{ab} P_i^{w\prime}(\xi) P_j^{w\prime}(\eta)\]

Throughout the kernels, the pair of indices \((i,j)\) names the row term, denoted \(A\), and \((k,l)\) names the column term, denoted \(B\); the variables fAw, gAweta, fBwxi and so forth follow that convention, read from the tables of the approximation functions evaluated once per integration point, see Approximation functions tabulated at the integration points. Only the upper triangle is assembled, the lower one being recovered by finalize_symmetric_matrix of structsolve, called by Shell.calc_kC() and Shell.calc_kG() with finalize=True, which is why every kernel skips the entries with row > col.

Kinematics as implemented

The Donnell kinematics of a cylindrical shell neglect the terms of the Sanders-Koiter kinematics that bring the in-plane displacements into the rotations and the curvatures, see the model cylshell_clpt_sanders, and the plate kinematics of von Kármán are recovered by making \(1/r = 0\):

(5)\[\begin{split}\boldsymbol{\varepsilon}^{(0)} = \left\{ \begin{matrix} u_{,x} + \frac{1}{2} w_{,x}^2 \\[2pt] v_{,y} + \frac{w}{r} + \frac{1}{2} w_{,y}^2 \\[2pt] u_{,y} + v_{,x} + w_{,x} w_{,y} \end{matrix} \right\} \qquad \boldsymbol{\varepsilon}^{(1)} = \left\{ \begin{matrix} -w_{,xx} \\[2pt] -w_{,yy} \\[2pt] -2 w_{,xy} \end{matrix} \right\}\end{split}\]

(5) is the strain state assembled term by term in calc_fint of cylshell_clpt_donnell_num.pyx, where the hoop term appears as 1./r*cs[col+2]*fAw*gAw and the von Kármán contributions as the products with wxi and weta. The plate, calc_fint of plate_clpt_donnell_num.pyx, differs from it in exactly one place, the absence of that hoop term in eyy; the five remaining strain components are character for character identical. The same single-term difference holds in fkG_num.

On presenting the plate as the limit of the cylinder

The cylinder is the general case and the plate is its \(1/r \to 0\) limit, but this statement must be read as a limit and not as a substitution, for two reasons found in the code.

First, the plate is a separate implementation, not a specialisation evaluated at a particular radius: panels/modelDB.py binds each model to its own pair of modules, and no plate quantity is ever obtained by calling a cylinder kernel.

Second, the limit is \(r \to \infty\) and not \(r = 0\). The cylindrical kernels contain \(1/r\) and \(1/r^2\) in denominators, for instance the term 0.25*A22*a*b*fAwfBw*gAwgBw/(r*r) of fk0 in cylshell_clpt_donnell.pyx and the corresponding terms of fkC_num. Evaluating them at \(r = 0\) produces infinities, which would only surface later, as an assertion inside finalize_symmetric_matrix rejecting the matrix. The radius is therefore validated before any kernel is called. Each model declares in panels/modelDB.py whether it needs one, through the key requires_r, set for the cylinders and cleared for the plates. Shell._check_r() raises a ValueError naming the model when requires_r is set and \(r\) is None or \(r \leq 0\), and it is called by the methods that reach the kernels, e.g. Shell.calc_kC(), Shell.calc_kG(), Shell.calc_kM(), Shell.calc_kA() and Shell.calc_fint(). A shell declared with model='cylshell_clpt_donnell' and r=None is thus rejected with a message pointing at plate_clpt_donnell, rather than being evaluated at a radius that is the opposite of the plate limit.

For the plate model, whose kernels never read r, an unset radius is normalised to zero, the value the remaining code compares against. There \(r = 0\) is not a radius but a marker for the absence of one.

Numerically, the convergence of the cylinder towards the plate is first order in \(1/r\), as expected from (5), where the hoop term enters linearly. Comparing the two implementations on the same \(\bar{\boldsymbol{c}}\), with \(a = 0.3\) m, \(b = 0.2\) m, \(m = n = 4\) and an unsymmetric laminate, gives the relative differences:

Relative difference between the cylindrical and the plate kernels as the radius grows, at a fixed \(\bar{\boldsymbol{c}}\)

\(r\) [m]	\(\|\boldsymbol{K}_C^{cyl} - \boldsymbol{K}_C^{plate}\| / \|\boldsymbol{K}_C^{plate}\|\)	\(\|\boldsymbol{F}_{int}^{cyl} - \boldsymbol{F}_{int}^{plate}\| / \|\boldsymbol{F}_{int}^{plate}\|\)
\(10^{0}\)	\(3.58 \times 10^{-3}\)	\(2.26 \times 10^{-2}\)
\(10^{2}\)	\(3.58 \times 10^{-5}\)	\(2.26 \times 10^{-4}\)
\(10^{4}\)	\(3.58 \times 10^{-7}\)	\(2.26 \times 10^{-6}\)
\(10^{6}\)	\(3.58 \times 10^{-9}\)	\(2.26 \times 10^{-8}\)
\(10^{8}\)	\(3.58 \times 10^{-11}\)	\(2.26 \times 10^{-10}\)

In the equations that follow, the cylinder is written as the general case and the plate is recovered by dropping the terms carrying \(1/r\), which is a statement about the equations and not an instruction about the code.

Constitutive relation

The generalised strain and stress vectors are related by the laminate stiffness:

(6)\[\begin{split}\hat{\boldsymbol{\varepsilon}} = \left\{ \begin{matrix} \boldsymbol{\varepsilon}^{(0)} \\ \boldsymbol{\varepsilon}^{(1)} \end{matrix} \right\} \qquad \hat{\boldsymbol{\sigma}} = \left\{ \begin{matrix} \boldsymbol{N} \\ \boldsymbol{M} \end{matrix} \right\} = \begin{bmatrix} \boldsymbol{A} & \boldsymbol{B} \\ \boldsymbol{B} & \boldsymbol{D} \end{bmatrix} \left\{ \begin{matrix} \boldsymbol{\varepsilon}^{(0)} \\ \boldsymbol{\varepsilon}^{(1)} \end{matrix} \right\} = \boldsymbol{C} \hat{\boldsymbol{\varepsilon}}\end{split}\]

with \(\boldsymbol{N} = \{N_{xx} \quad N_{yy} \quad N_{xy}\}^\top\) and \(\boldsymbol{M} = \{M_{xx} \quad M_{yy} \quad M_{xy}\}^\top\). (6) is applied in calc_fint of both models. The \(6 \times 6\) array \(\boldsymbol{C}\) is the argument named Finput of every kernel, read from the attribute Shell.ABDnxny when it is set, else from the read-only property Shell.ABD, calculated by Shell._get_lam_ABD() from the laminate Shell.lam. When Finput has the shape \((n_x, n_y, 6, 6)\) it is interpreted as one constitutive matrix per integration point, which is how variable stiffness laminates enter, see Shell.integration_points(); a \((6, 6)\) array is used for the whole domain instead. Everything that follows holds pointwise, so a \(\boldsymbol{C}\) varying over the domain changes nothing in the derivation.

Linear and non-linear differentiation operators

The linear part of (5) is written as \(\hat{\boldsymbol{\varepsilon}}_L = \boldsymbol{B}_L \bar{\boldsymbol{c}}\), with:

(7)\[\begin{split}\boldsymbol{B}_L = \begin{bmatrix} \boldsymbol{S}^u_{,x} \\ \boldsymbol{S}^v_{,y} + \frac{1}{r}\boldsymbol{S}^w \\ \boldsymbol{S}^u_{,y} + \boldsymbol{S}^v_{,x} \\ -\boldsymbol{S}^w_{,xx} \\ -\boldsymbol{S}^w_{,yy} \\ -2\boldsymbol{S}^w_{,xy} \end{bmatrix}\end{split}\]

(7) is exactly the set of six accumulators of the strains in fkG_num, which is the one place in the code where the linear operator appears on its own.

For the non-linear part, it is convenient to define the operator \(\boldsymbol{G}\) that extracts the slopes of the deflection:

(8)\[\begin{split}\boldsymbol{G} = \begin{bmatrix} \boldsymbol{G}_1 \\ \boldsymbol{G}_2 \end{bmatrix} = \begin{bmatrix} \boldsymbol{S}^w_{,x} \\ \boldsymbol{S}^w_{,y} \end{bmatrix} \qquad w_{,x} = \boldsymbol{G}_1 \bar{\boldsymbol{c}} \qquad w_{,y} = \boldsymbol{G}_2 \bar{\boldsymbol{c}}\end{split}\]

(8) is the Donnell-Mushtari-Vlasov specialisation of the rotation operators: the Sanders-Koiter kinematics add the terms \(\boldsymbol{S}^u/R_1\) and \(\boldsymbol{S}^v/R_2\) of the in-plane displacements, with the principal radii of curvature \(R_1\) and \(R_2\). The code confirms this: neither fkG_num nor fkG0 contains \(\boldsymbol{S}^v\), and both allocate \(1 \cdot m^2 n^2\) entries instead of the \(9 \cdot m^2 n^2\) of fkC_num (fdim), because only the \(w\)–\(w\) block is populated. This holds for the cylinder as well as for the plate: the \(1/r\) term present in \(\boldsymbol{B}_L\) has no counterpart in \(\boldsymbol{G}\).

In the kernels, the slopes are accumulated in the natural coordinates as wxi and weta, in fkC_num and in calc_fint, and converted as:

\[w_{,x} = \frac{2}{a}\,\texttt{wxi} \qquad w_{,y} = \frac{2}{b}\,\texttt{weta}\]

With (8), the non-linear differentiation operator becomes:

(9)\[\begin{split}\boldsymbol{B}_{NL}(\bar{\boldsymbol{c}}) = \begin{bmatrix} w_{,x} \boldsymbol{G}_1 \\ w_{,y} \boldsymbol{G}_2 \\ w_{,y} \boldsymbol{G}_1 + w_{,x} \boldsymbol{G}_2 \\ \boldsymbol{0} \\ \boldsymbol{0} \\ \boldsymbol{0} \end{bmatrix}\end{split}\]

which gives the strain and its variation:

(10)\[\begin{split}\hat{\boldsymbol{\varepsilon}} = \left( \boldsymbol{B}_L + \frac{1}{2}\boldsymbol{B}_{NL} \right) \bar{\boldsymbol{c}} \qquad \hat{\boldsymbol{\varepsilon}}_{NL} = \frac{1}{2}\boldsymbol{B}_{NL}\bar{\boldsymbol{c}} = \left\{ \begin{matrix} \frac{1}{2} w_{,x}^2 \\[2pt] \frac{1}{2} w_{,y}^2 \\[2pt] w_{,x} w_{,y} \\[2pt] 0 \\ 0 \\ 0 \end{matrix} \right\} \qquad \delta \hat{\boldsymbol{\varepsilon}} = (\boldsymbol{B}_L + \boldsymbol{B}_{NL}) \delta\bar{\boldsymbol{c}}\end{split}\]

The factor \(\frac{1}{2}\) of (10) applies to the shear row as well, and this is visible in the code as an asymmetry that is easy to mistake for an error: in the strains of calc_fint the two normal strains carry an explicit 0.5 whereas the shear strain does not. The reason is that the shear row of \(\boldsymbol{B}_{NL}\) already contains two terms, so \(\frac{1}{2}\boldsymbol{B}_{NL}\bar{\boldsymbol{c}}\) evaluates to \(\frac{1}{2}(2 w_{,x} w_{,y}) = w_{,x}w_{,y}\), which is what the single unscaled product in the code produces. Evaluating the non-linear shear strain as \(\boldsymbol{B}_{NL}\bar{\boldsymbol{c}}\) instead would give \(2 w_{,x} w_{,y}\). That error does not show in \(\boldsymbol{B}_{NL}\) itself, but it corrupts the stress state, and therefore the internal force vector and the geometric stiffness matrix, leaving a tangent stiffness matrix that is no longer the derivative of the internal force vector.

Internal force vector

The internal force vector is the derivative of the strain energy, \(U = \frac{1}{2} \int_{\mathcal{A}} \hat{\boldsymbol{\varepsilon}}^\top \boldsymbol{C} \hat{\boldsymbol{\varepsilon}} \, d\mathcal{A}\), integrated over the mid-surface, with respect to \(\bar{\boldsymbol{c}}\):

(11)\[\boldsymbol{F}_{int}(\bar{\boldsymbol{c}}) = \int_{\mathcal{A}} \left( \boldsymbol{B}_L + \boldsymbol{B}_{NL}(\bar{\boldsymbol{c}}) \right)^\top \hat{\boldsymbol{\sigma}}(\bar{\boldsymbol{c}}) \, d\mathcal{A}\]

which gives, for the three degrees of freedom associated with the term \((i,j)\):

(12)\[\begin{split}\begin{aligned} F_{int}^{u} &= \int_{\mathcal{A}} \left( N_{xx} S^u_{,x} + N_{xy} S^u_{,y} \right) d\mathcal{A} \\ F_{int}^{v} &= \int_{\mathcal{A}} \left( N_{yy} S^v_{,y} + N_{xy} S^v_{,x} \right) d\mathcal{A} \\ F_{int}^{w} &= \int_{\mathcal{A}} \Big( N_{xx} w_{,x} S^w_{,x} + N_{yy} \left( \tfrac{1}{r} S^w + w_{,y} S^w_{,y} \right) + N_{xy} \left( w_{,y} S^w_{,x} + w_{,x} S^w_{,y} \right) \\ &\qquad - M_{xx} S^w_{,xx} - M_{yy} S^w_{,yy} - 2 M_{xy} S^w_{,xy} \Big) d\mathcal{A} \end{aligned}\end{split}\]

(12) is the update of fint in calc_fint of cylshell_clpt_donnell_num.pyx. The hoop contribution \(N_{yy} S^w / r\) appears there as Nyy*(fAw*gAw/r + ...), and the plate form, in plate_clpt_donnell_num.pyx, is the same expression without it. Note that \(\boldsymbol{B}_L\) contributes the hoop term and the three curvature terms to the \(w\) row, whereas \(\boldsymbol{B}_{NL}\) contributes the three products of membrane forces with slopes; the \(u\) and \(v\) rows receive no contribution from \(\boldsymbol{B}_{NL}\), which is why they contain only \(\boldsymbol{N}\) and no slope.

Tangent stiffness matrix

Differentiating (11) with respect to \(\bar{\boldsymbol{c}}\), (1):

(13)\[\boldsymbol{K}_T = \frac{\partial \boldsymbol{F}_{int}}{\partial \bar{\boldsymbol{c}}} = \underbrace{\int_{\mathcal{A}} (\boldsymbol{B}_L + \boldsymbol{B}_{NL})^\top \boldsymbol{C} (\boldsymbol{B}_L + \boldsymbol{B}_{NL}) \, d\mathcal{A}}_{\text{from } \partial \hat{\boldsymbol{\sigma}} / \partial \bar{\boldsymbol{c}}} + \underbrace{\int_{\mathcal{A}} \frac{\partial \boldsymbol{B}_{NL}^\top}{\partial \bar{\boldsymbol{c}}} \hat{\boldsymbol{\sigma}} \, d\mathcal{A}}_{\text{from } \partial \boldsymbol{B}_{NL} / \partial \bar{\boldsymbol{c}}}\]

The first integral uses \(\partial \hat{\boldsymbol{\sigma}} / \partial \bar{\boldsymbol{c}} = \boldsymbol{C}\,(\boldsymbol{B}_L + \boldsymbol{B}_{NL})\), which follows from (10) together with \(\partial (\boldsymbol{B}_{NL}\bar{\boldsymbol{c}}) / \partial \bar{\boldsymbol{c}} = 2\boldsymbol{B}_{NL}\), itself a consequence of \(\boldsymbol{B}_{NL}\) being homogeneous of degree one in \(\bar{\boldsymbol{c}}\). The second integral is evaluated by writing \(\boldsymbol{B}_{NL}^\top \hat{\boldsymbol{\sigma}}\) explicitly with (9):

\[\boldsymbol{B}_{NL}^\top \hat{\boldsymbol{\sigma}} = \boldsymbol{G}_1^\top \left( N_{xx} w_{,x} + N_{xy} w_{,y} \right) + \boldsymbol{G}_2^\top \left( N_{xy} w_{,x} + N_{yy} w_{,y} \right)\]

and differentiating with respect to \(\bar{\boldsymbol{c}}\) at a fixed stress state, using (8):

(14)\[\begin{split}\boldsymbol{K}_G(\boldsymbol{N}) = \int_{\mathcal{A}} \boldsymbol{G}^\top \tilde{\boldsymbol{N}} \boldsymbol{G} \, d\mathcal{A} \qquad \tilde{\boldsymbol{N}} = \begin{bmatrix} N_{xx} & N_{xy} \\ N_{xy} & N_{yy} \end{bmatrix}\end{split}\]

(14) is the geometric stiffness of the Donnell-Mushtari-Vlasov kinematics, and it only has non-zero entries in the \(w\)–\(w\) blocks. For the term pairs \(A = (i,j)\) and \(B = (k,l)\):

(15)\[K_{G\,[3(jm+i)+2,\;3(lm+k)+2]} = \int_{\mathcal{A}} \left( N_{xx} S^w_{A,x} S^w_{B,x} + N_{xy} \left( S^w_{A,x} S^w_{B,y} + S^w_{A,y} S^w_{B,x} \right) + N_{yy} S^w_{A,y} S^w_{B,y} \right) d\mathcal{A}\]

(15) is transcribed literally in fkG_num of cylshell_clpt_donnell_num.pyx, where the four products of shape functions appear in the same order. The plate writes the same expression in an algebraically expanded form; the two are identical after collecting the powers of \(a\) and \(b\).

Which stress enters the geometric stiffness, and the degree of homogeneity

The first integral of (13) expands into the four classical groups of (2):

(16)\[\boldsymbol{K}_0 = \int_{\mathcal{A}} \boldsymbol{B}_L^\top \boldsymbol{C} \boldsymbol{B}_L \, d\mathcal{A} \qquad \boldsymbol{K}_{0L} = \int_{\mathcal{A}} \boldsymbol{B}_L^\top \boldsymbol{C} \boldsymbol{B}_{NL} \, d\mathcal{A} \qquad \boldsymbol{K}_{L0} = \boldsymbol{K}_{0L}^\top \qquad \boldsymbol{K}_{LL} = \int_{\mathcal{A}} \boldsymbol{B}_{NL}^\top \boldsymbol{C} \boldsymbol{B}_{NL} \, d\mathcal{A}\]

The membrane stress of (14) is also split, according to the strain that produces it:

(17)\[\begin{split}\boldsymbol{N} = \boldsymbol{N}^0 + \boldsymbol{N}_L + \boldsymbol{N}_{NL} \qquad \boldsymbol{N}_L = \boldsymbol{A}\boldsymbol{\varepsilon}^{(0)}_L + \boldsymbol{B}\boldsymbol{\varepsilon}^{(1)} \qquad \boldsymbol{N}_{NL} = \boldsymbol{A} \left\{ \begin{matrix} \frac{1}{2} w_{,x}^2 \\[2pt] \frac{1}{2} w_{,y}^2 \\[2pt] w_{,x} w_{,y} \end{matrix} \right\}\end{split}\]

where \(\boldsymbol{N}^0\) is an optional prescribed pre-stress, used for instance to superpose combined load cases, see c_cte and Nxx_cte, Nyy_cte, Nxy_cte in Shell.calc_kC(), and \(\boldsymbol{B}\) does not contribute to \(\boldsymbol{N}_{NL}\) because \(\boldsymbol{\varepsilon}^{(1)}\) is linear for the CLPT. In the code, \(\boldsymbol{N}^0 + \boldsymbol{N}_L\) is computed in fkG_num, from the arguments Nxx0, Nyy0, Nxy0 and from the linear strains of (7); and \(\boldsymbol{N}_{NL}\) is computed separately, in fkC_num, as NxxNL, NyyNL, NxyNL, where the 0.5 factors of (17) are explicit. Since (14) is linear in \(\boldsymbol{N}\):

(18)\[\boldsymbol{K}_G(\boldsymbol{N}) = \boldsymbol{K}_G(\boldsymbol{N}^0 + \boldsymbol{N}_L) + \boldsymbol{K}_{GNL} \qquad \text{with} \qquad \boldsymbol{K}_{GNL} = \boldsymbol{K}_G(\boldsymbol{N}_{NL})\]

With respect to \(\bar{\boldsymbol{c}}\): \(\boldsymbol{K}_0\) is constant; \(\boldsymbol{K}_{0L}\), \(\boldsymbol{K}_{L0}\) and \(\boldsymbol{K}_G(\boldsymbol{N}_L)\) are homogeneous of degree one; and \(\boldsymbol{K}_{LL}\) and \(\boldsymbol{K}_{GNL}\) are homogeneous of degree two. The implementation groups them as:

(19)\[\begin{split}\begin{aligned} \texttt{fkC_num} \;\rightarrow\; \boldsymbol{K}_C(\bar{\boldsymbol{c}}) &= \boldsymbol{K}_0 + \boldsymbol{K}_{0L} + \boldsymbol{K}_{L0} + \boldsymbol{K}_{LL} + \boldsymbol{K}_{GNL} \\ \texttt{fkG_num} \;\rightarrow\; \boldsymbol{K}_G(\bar{\boldsymbol{c}}) &= \boldsymbol{K}_G(\boldsymbol{N}^0 + \boldsymbol{N}_L) \\ \texttt{calc_kT} \;\rightarrow\; \boldsymbol{K}_T(\bar{\boldsymbol{c}}) &= \boldsymbol{K}_C(\bar{\boldsymbol{c}}) + \boldsymbol{K}_G(\bar{\boldsymbol{c}}) \end{aligned}\end{split}\]

(19) is the statement written as a comment in fkC_num of both models, and the summation \(\boldsymbol{K}_T = \boldsymbol{K}_C + \boldsymbol{K}_G\) happens in Shell.calc_kT(). The four constitutive groups of (16) are not assembled separately: the kernel evaluates the whole product \((\boldsymbol{B}_L + \boldsymbol{B}_{NL})^\top \boldsymbol{C} (\boldsymbol{B}_L + \boldsymbol{B}_{NL})\) in expanded symbolic form, in nine blocks, one per pair of displacements, and \(\boldsymbol{K}_{GNL}\) is then added on top of the \(w\)–\(w\) block only, under the comment # KGNL.

Why \(\boldsymbol{K}_{GNL}\) sits where it does. Although \(\boldsymbol{K}_{GNL}\) is geometric in nature, being a \(\boldsymbol{G}^\top \tilde{\boldsymbol{N}} \boldsymbol{G}\) product, it is deliberately collected together with the constitutive terms so that \(\boldsymbol{K}_G(\bar{\boldsymbol{c}})\) stays homogeneous of degree one:

\[\boldsymbol{K}_G(\lambda\bar{\boldsymbol{c}}) = \lambda\,\boldsymbol{K}_G(\bar{\boldsymbol{c}})\]

The linear buckling problem \((\boldsymbol{K}_0 + \lambda \boldsymbol{K}_G(\bar{\boldsymbol{c}}_0)) \boldsymbol{q} = \boldsymbol{0}\) needs this property, because there the load multiplier \(\lambda\) scales a stress state obtained from a linear fundamental state \(\bar{\boldsymbol{c}}_0\), and a \(\boldsymbol{K}_G\) containing a degree-two term would make \(\lambda\) scale that state non-linearly, such that the problem would no longer be a linear eigenvalue problem. With the grouping of (19), both requirements are met at once: \(\boldsymbol{K}_T\) is the exact Jacobian of \(\boldsymbol{F}_{int}\), so the Newton-Raphson iterations converge quadratically, and \(\boldsymbol{K}_G\) remains usable as the geometric stiffness of a linear buckling analysis. The consequence is that fkG_num has no switch for \(\boldsymbol{N}_{NL}\) at all: its signature takes no NLgeom argument, and its docstring states why. Shell.calc_kG() still accepts NLgeom, but only to force the numerically integrated path of The two implementation paths; it does not change the value of \(\boldsymbol{K}_G\).

Mapping onto the three-way split

With (19) established, the three-way split of (3) is:

(20)\[\begin{split}\begin{aligned} \boldsymbol{K}_{C0} &= \boldsymbol{K}_0 \\ \boldsymbol{K}_{CNL}(\bar{\boldsymbol{c}}) &= \boldsymbol{K}_{0L} + \boldsymbol{K}_{L0} + \boldsymbol{K}_{LL} + \boldsymbol{K}_{GNL} \\ \boldsymbol{K}_G(\bar{\boldsymbol{c}}) &= \boldsymbol{K}_G(\boldsymbol{N}^0 + \boldsymbol{N}_L) \end{aligned}\end{split}\]

i.e. \(\boldsymbol{K}_{C0}\) is everything that fkC_num produces at \(\bar{\boldsymbol{c}} = \boldsymbol{0}\), \(\boldsymbol{K}_{CNL}\) is everything else that fkC_num produces, and \(\boldsymbol{K}_G\) is what fkG_num produces. The partition is clean, but only under that reading of \(\boldsymbol{K}_{CNL}\), which is not the obvious one; Where the naming and the code do not line up states the caveat.

There is no separate code path for \(\boldsymbol{K}_{C0}\) and \(\boldsymbol{K}_{CNL}\): fkC_num returns their sum. \(\boldsymbol{K}_{C0}\) alone is obtained either from the closed-form fk0 or by calling fkC_num with \(\bar{\boldsymbol{c}} = \boldsymbol{0}\) or with NLgeom=0, both of which make wxi and weta vanish, since their accumulation is guarded by if NLgeom == 1. With \(w_{,x} = w_{,y} = 0\) one has \(\boldsymbol{B}_{NL} = \boldsymbol{0}\) and \(\boldsymbol{N}_{NL} = \boldsymbol{0}\), so \(\boldsymbol{K}_{CNL}\) vanishes identically, as a comment in fkC_num notes.

The closed-form path

The closed-form kernels evaluate the same \(\boldsymbol{K}_0\) and \(\boldsymbol{K}_G\), over the full domain and for a constant laminate, by separating the two directions. Every entry of (16) restricted to \(\boldsymbol{B}_L\) is a sum of products of the form

(21)\[\int_{\mathcal{A}} \frac{\partial^{p} S^{\alpha}_A}{\partial x^{p}} \frac{\partial^{q} S^{\beta}_B}{\partial y^{q}} \, d\mathcal{A} = \frac{ab}{4} \left( \frac{2}{a} \right)^{p} \left( \frac{2}{b} \right)^{q} \int_{-1}^{+1} P_i^{\alpha (p)} P_k^{\beta} \, d\xi \int_{-1}^{+1} P_j^{\alpha} P_l^{\beta (q)} \, d\eta\]

and the two one-dimensional integrals are evaluated analytically by the six functions integral_ff, integral_ffp, integral_ffpp, integral_fpfp, integral_fpfpp and integral_fppfpp of panels/core/include/bardell.hpp, declared at the top of plate_clpt_donnell.pyx and wrapped for Python as calc_integral_ff and so forth in panels/bardell.pyx. They are hoisted out of the inner loops, the \(\xi\) integrals inside the \(i\), \(k\) loops and the \(\eta\) integrals inside the \(j\), \(l\) loops, which is what makes this path cheap: no loop over integration points exists at all.

As an example, the first entry of \(\boldsymbol{K}_0\) in fk0 of plate_clpt_donnell.pyx reads

k0v[c] += A11*b*fAuxifBuxi*gAugBu/a + A16*(fAufBuxi*gAuetagBu
        + fAuxifBu*gAugBueta) + A66*a*fAufBu*gAuetagBueta/b

which is the \(u\)–\(u\) entry of \(\int \boldsymbol{B}_L^\top \boldsymbol{C} \boldsymbol{B}_L \, d\mathcal{A}\), the factor \(b/a\) following from (21) with \(p = q = 1\). The corresponding cylindrical entry, in fk0 of cylshell_clpt_donnell.pyx, is identical, the two implementations of \(\boldsymbol{K}_0\) differing only by the five entries carrying \(1/r\) and \(1/r^2\) that (7) introduces through \(\boldsymbol{S}^w/r\).

Two properties of this path bound its use:

• fk0 takes no \(\bar{\boldsymbol{c}}\) and no constitutive matrix per point; its signature is fk0(shell, size, row0, col0). It can only ever produce \(\boldsymbol{K}_{C0}\).

• fkG0 takes the three membrane stress resultants as scalars, fkG0(Nxx, Nyy, Nxy, shell, size, row0, col0), so it evaluates \(\boldsymbol{K}_G(\boldsymbol{N}^0)\) for a membrane state uniform over the domain. The bodies of fkG0 are identical for the plate and for the cylinder, consistent with \(\boldsymbol{G}\) of (8) carrying no \(1/r\); the only entry assembled is character for character the same in both modules.

The geometric stiffness must therefore be integrated with fkG_num, and not fkG0, whenever the pre-buckling resultants are not uniform, e.g. near the edges of a cylinder whose fundamental state is not a pure membrane state: fkG0 cannot represent it.

Numerical integration

All integrals of the numerically integrated path are computed in natural coordinates using the Gauss-Legendre quadrature, with \(n_x\) and \(n_y\) points obtained from scipy.special.roots_legendre:

(22)\[\int_{\mathcal{A}} f \, d\mathcal{A} = \frac{ab}{4} \int_{-1}^{+1} \int_{-1}^{+1} f(\xi, \eta) \, d\xi \, d\eta \approx \frac{ab}{4} \sum_{p=1}^{n_x} \sum_{q=1}^{n_y} W_p W_q \, f(\xi_p, \eta_q)\]

For integration over a sub-domain \(x_1 \le x \le x_2\) and \(y_1 \le y \le y_2\), the factor \(ab/4\) becomes \((x_2-x_1)(y_2-y_1)/4\), computed as intx*inty, and the Gauss points are mapped into the corresponding interval of \(\xi\) and \(\eta\), the interval itself being given by the physical limits returned by Shell.integration_limits() and mapped unconditionally, which gives exactly \(\xi, \eta = \pm 1\) for the full domain; Shell.integration_points() returns the same points and weights in physical coordinates. The factor \(1/4\) is not always visible in the kernels: it is sometimes folded into the symbolic expansion, as in the entries of fkC_num of plate_clpt_donnell_num.pyx, and sometimes written out, as weight*(intx*inty/4) in fkG_num of cylshell_clpt_donnell_num.pyx. Both forms are the same number.

\(\boldsymbol{F}_{int}\) and \(\boldsymbol{K}_T\) must be integrated with the same quadrature rule, which is why Shell.calc_fint() and Shell.calc_kT() both default \(n_x\) and \(n_y\) to the attributes nx and ny of the shell, or to the shape of ABDnxny when it has one matrix per point. The consistency between them is then exact at the discrete level, whether or not the quadrature integrates the polynomials exactly: (13) is the derivative of (11) pointwise, hence also after applying the same weighted sum to both.

Multi-domain assembly

For an assembly of panels, each Shell writes its own contribution into a global matrix of size \(\sum_p 3 m_p n_p\), offset by the attributes row_start and col_start and passed through as the arguments row0 and col0 of every kernel. Denoting by \(\boldsymbol{L}_p\) the Boolean operator that selects the degrees of freedom of panel \(p\) from the global vector, the assembly is:

(23)\[\begin{split}\begin{aligned} \boldsymbol{F}_{int} &= \sum_{p} \boldsymbol{L}_p^\top \boldsymbol{F}_{int}^{(p)}(\boldsymbol{L}_p \bar{\boldsymbol{c}}) + \boldsymbol{K}_{C}^{conn} \bar{\boldsymbol{c}} \\ \boldsymbol{K}_T &= \sum_{p} \boldsymbol{L}_p^\top \left( \boldsymbol{K}_{C}^{(p)} + \boldsymbol{K}_{G}^{(p)} \right) \boldsymbol{L}_p + \boldsymbol{K}_{C}^{conn} \end{aligned}\end{split}\]

(23) is MultiDomain.calc_fint() and MultiDomain.calc_kT(). Three points are relevant here:

• The panel-level blocks enter unchanged. MultiDomain.calc_kT() calls p.calc_kC(..., NLgeom=True) and p.calc_kG(..., NLgeom=True) for each panel, i.e. exactly the two calls of Shell.calc_kT(), so (19) holds panel by panel.

• Both are assembled with finalize=False and the symmetrisation is applied once to the sum, because each panel only fills its own upper triangle. The load stiffness of the follower pressure loads, which is unsymmetric, is added after the symmetrisation, see Follower pressure loads: theory.

• The connection matrix \(\boldsymbol{K}_{C}^{conn}\) is added after the loop, to \(\boldsymbol{K}_T\) and, as \(\boldsymbol{K}_{C}^{conn} \bar{\boldsymbol{c}}\), to \(\boldsymbol{F}_{int}\). Its own derivation is out of the scope of this section, see Exact connections with the null-space method; what matters here is that it is added consistently to both, so that the assembled \(\boldsymbol{K}_T\) remains the derivative of the assembled \(\boldsymbol{F}_{int}\) whenever \(\boldsymbol{K}_{C}^{conn}\) is independent of \(\bar{\boldsymbol{c}}\).

Symbols and implementation variables

The following tables relate the symbols used herein to the variables of the implementation, those of the plate kernels unless stated, the cylindrical twins using the same names, and map the matrix terms.

Symbols and implementation variables

Symbol	Variable	Where
\(\bar{\boldsymbol{c}}\)	cs, c	argument of the kernels
\(m\), \(n\)	m, n	attributes of Shell, Shell.get_size()
\(a\), \(b\), \(r\)	a, b, r	attributes of Shell
\(P_i^u(\xi)\), \(P_i^{u\prime}(\xi)\)	fAu, fAuxi	ftab of the kernels
\(P_j^w(\eta)\), \(P_j^{w\prime}\), \(P_j^{w\prime\prime}\)	gAw, gAweta, gAwetaeta	gtab of the kernels
\(P_k^w(\xi)\), \(P_l^w(\eta)\) (term \(B\))	fBw, gBw	ftab, gtab of the kernels
\(\delta_{t1}\), \(\delta_{r1}\), \(\delta_{t2}\), \(\delta_{r2}\)	x1w, x1wr, x2w, x2wr	attributes of Shell
\(\partial w/\partial\xi\), \(\partial w/\partial\eta\)	wxi, weta	fkC_num, calc_fint
\(w_{,x}\), \(w_{,y}\)	wx, wy	fkC_num
\(\boldsymbol{\varepsilon}^{(0)}\), \(\boldsymbol{\varepsilon}^{(1)}\)	exx, eyy, gxy, kxx, kyy, kxy	calc_fint, fkG_num
\(\boldsymbol{C}\), i.e. \(\boldsymbol{A}\), \(\boldsymbol{B}\), \(\boldsymbol{D}\)	Finput, F, A11 … D66	argument of the kernels, Shell.ABD, Shell.ABDnxny
\(\boldsymbol{N}\), \(\boldsymbol{M}\)	Nxx, Nyy, Nxy, Mxx, Myy, Mxy	calc_fint
\(\boldsymbol{N}_{NL}\)	NxxNL, NyyNL, NxyNL	fkC_num
\(\boldsymbol{N}^0\)	Nxx0, Nyy0, Nxy0; Shell.Nxx …	arguments of fkG_num
\(\boldsymbol{N}^0 + \boldsymbol{N}_L\)	Nxx, Nyy, Nxy of fkG_num	fkG_num
\(x_2-x_1\), \(y_2-y_1\)	intx, inty	numerical kernels
\(W_p W_q\)	weight	numerical kernels
\(n_x\), \(n_y\)	nx, ny	attributes of Shell, arguments of the calc methods
row term \(A = (i,j)\), column term \(B = (k,l)\)	i, j, k, l	kernels
\(3(jm+i)\)	DOF*(j*m + i)	row, col of the kernels
offsets of \(\boldsymbol{L}_p\)	row0, col0	row_start, col_start in MultiDomain.calc_kT()

The two last rows of the next table hold the terms whose placement differs between the two namings.

Matrix terms: notation of this section, name used by the code and producing function

Herein	Code name	Function
\(\boldsymbol{K}_{C0}\)	K0	fk0
\(\boldsymbol{K}_{C0}\)	K0	fkC_num with NLgeom=0
\(\boldsymbol{K}_{CNL}\)	K0L + KL0 + KLL	fkC_num with NLgeom=1
\(\boldsymbol{K}_G\)	KG(N0 + N_L)	fkG_num
\(\boldsymbol{K}_G(\boldsymbol{N}^0)\)	KG0	fkG0
\(\boldsymbol{F}_{int}\)	fint	calc_fint
\(\boldsymbol{K}_T\)	kT	Shell.calc_kT()
\(\boldsymbol{K}_{CNL}\), part of	KGNL	fkC_num, under # KGNL
not in \(\boldsymbol{K}_T\)	KG0, KG(c_cte)	Shell.calc_kC(), added to kC

Verification

Three test modules pin the claims of this section. The first two are parametrised over plate_clpt_donnell and cylshell_clpt_donnell, both using an unsymmetric, unbalanced laminate so that \(\boldsymbol{B} \neq \boldsymbol{0}\) and \(A_{16}, A_{26} \neq 0\).

tests/tests_shell/test_tangent_consistency.py verifies the Jacobian claim of (13) directly:

• test_tangent_is_derivative_of_fint is the directional Taylor test. Given an arbitrary state \(\bar{\boldsymbol{c}}\) and direction \(\boldsymbol{d}\), \(\boldsymbol{F}_{int}(\bar{\boldsymbol{c}} + h\boldsymbol{d}) = \boldsymbol{F}_{int}(\bar{\boldsymbol{c}}) + h \boldsymbol{K}_T(\bar{\boldsymbol{c}}) \boldsymbol{d} + O(h^2)\), so the relative error

  \[e(h) = \frac{\left\| \boldsymbol{F}_{int}(\bar{\boldsymbol{c}} + h\boldsymbol{d}) - \boldsymbol{F}_{int}(\bar{\boldsymbol{c}}) - h \boldsymbol{K}_T \boldsymbol{d} \right\|}{\left\| h \boldsymbol{K}_T \boldsymbol{d} \right\|}\]

  must fall proportionally to \(h\). The test asserts a factor of at least five per decade over \(h \in [10^{-6}, 10^{-2}]\), and that \(e(h)/h\) does not drift by more than a factor of two. An inconsistent tangent leaves a remainder of order \(h\) and \(e(h)\) plateaus instead.

• test_tangent_matches_finite_difference_jacobian compares every entry of \(\boldsymbol{K}_T\) against a central finite-difference Jacobian of \(\boldsymbol{F}_{int}\).

• test_kG_is_homogeneous_of_degree_one asserts \(\boldsymbol{K}_G(2\bar{\boldsymbol{c}}) = 2\boldsymbol{K}_G(\bar{\boldsymbol{c}})\), which is the property that the placement of \(\boldsymbol{K}_{GNL}\) in (19) exists to preserve.

• test_fint_and_tangent_in_the_undeformed_state asserts \(\boldsymbol{F}_{int}(\boldsymbol{0}) = \boldsymbol{0}\) and \(\boldsymbol{K}_T(\boldsymbol{0}) = \boldsymbol{K}_{C0}\), which is (20) at \(\bar{\boldsymbol{c}} = \boldsymbol{0}\).

• test_newton_raphson_converges_quadratically manufactures an equilibrium state and checks that the error is squared at each iteration.

tests/tests_shell/test_nonlinear.py exercises the same tangent in a full Newton-Raphson solution, deep in the postbuckling regime at about 37 times the first linear buckling load, and asserts both the converged deflection and a convergence order above \(1.6\).

tests/tests_shell/test_argument_validation.py pins the argument handling around these matrices, and with it three claims of this section:

• test_KGNL_is_in_calc_kC_and_not_in_calc_kG asserts that \(\boldsymbol{K}_G\) is homogeneous of degree one, that NLgeom does not change its value, and that \(\boldsymbol{K}_T\) is exactly the sum of the two matrices of (19).

• test_NLgeom_and_not_c_selects_the_large_displacement_matrix asserts that fkC_num returns \(\boldsymbol{K}_{C0}\) whenever NLgeom=0, whatever \(\bar{\boldsymbol{c}}\), and that the result then matches the closed-form fk0, as stated when mapping the code onto (20).

• test_cylshell_approaches_plate_as_radius_grows asserts the monotonic convergence of On presenting the plate as the limit of the cylinder, and test_cylshell_requires_positive_radius that \(r = 0\), a negative \(r\) and an unset \(r\) are all rejected for the cylinder.

Together the first two modules support every claim made here about \(\boldsymbol{K}_T\) being the exact Jacobian, about \(\boldsymbol{K}_G\) being homogeneous of degree one, and about \(\boldsymbol{K}_T(\boldsymbol{0}) = \boldsymbol{K}_{C0}\). They do not separately test the decomposition of (19) term by term, nor the closed-form path against the numerically integrated one; those were verified independently while preparing this section, by rebuilding \(\boldsymbol{B}_L\), \(\boldsymbol{B}_{NL}\) and \(\boldsymbol{G}\) from (7), (9) and (8) at the Gauss points and comparing against the compiled kernels, with the following results, for both models, \(m = n = 4\), \(n_x = n_y = 10\), an unsymmetric laminate and a random \(\bar{\boldsymbol{c}}\) with deflections of twenty thicknesses.

Relative difference between the compiled kernels and the equations of this section, rebuilt independently

Kernel	Equation	Plate	Cylinder
fkC_num	\(\int (\boldsymbol{B}_L + \boldsymbol{B}_{NL})^\top \boldsymbol{C} (\boldsymbol{B}_L + \boldsymbol{B}_{NL}) + \boldsymbol{K}_{GNL}\)	\(4.0 \times 10^{-16}\)	\(2.3 \times 10^{-16}\)
fkG_num	\(\boldsymbol{K}_G(\boldsymbol{N}^0 + \boldsymbol{N}_L)\), (14)	\(5.7 \times 10^{-16}\)	\(4.9 \times 10^{-16}\)
calc_fint	\(\int (\boldsymbol{B}_L + \boldsymbol{B}_{NL})^\top \hat{\boldsymbol{\sigma}}\), (11)	\(2.2 \times 10^{-16}\)	\(2.2 \times 10^{-16}\)
fk0	fkC_num at \(\bar{\boldsymbol{c}} = \boldsymbol{0}\), \(n_x = n_y = 30\)	\(4.7 \times 10^{-15}\)	\(4.6 \times 10^{-15}\)
fkG0	fkG_num at \(\bar{\boldsymbol{c}} = \boldsymbol{0}\), \(n_x = n_y = 30\)	\(8.6 \times 10^{-15}\)	\(8.8 \times 10^{-15}\)

The last two rows also establish that the closed-form and the numerically integrated paths agree to round-off once the quadrature is fine enough to integrate the polynomials exactly, so the choice between them of The two implementation paths is a matter of cost and of capability, not of accuracy.

Where the naming and the code do not line up

Two points deserve to be stated rather than smoothed over.

The \(\boldsymbol{K}_C\)–\(\boldsymbol{K}_G\) boundary moves. The split of (3) and the split of (2) are not the same split, although both sum to the same \(\boldsymbol{K}_T\). Writing \(\boldsymbol{K}_C\) and \(\boldsymbol{K}_G\) for the classical constitutive and geometric stiffness matrices, defined respectively as the first and the second integral of (13):

(24)\[\boldsymbol{K}_{C0} + \boldsymbol{K}_{CNL} = \boldsymbol{K}_C + \boldsymbol{K}_{GNL} \qquad \boldsymbol{K}_G^{\text{herein}} = \boldsymbol{K}_G - \boldsymbol{K}_{GNL}\]

The term \(\boldsymbol{K}_{GNL}\) is geometric by construction, being \(\int \boldsymbol{G}^\top \tilde{\boldsymbol{N}}_{NL} \boldsymbol{G} \, d\mathcal{A}\), yet under the naming of (3) it belongs to \(\boldsymbol{K}_{CNL}\). The three-way split therefore partitions the five terms of the code cleanly, but only if \(\boldsymbol{K}_{CNL}\) is defined operationally, as everything fkC_num adds beyond \(\bar{\boldsymbol{c}} = \boldsymbol{0}\), and not semantically, as the non-linear constitutive terms. Under the semantic reading, \(\boldsymbol{K}_{GNL}\) has no home in any of the three symbols and the partition fails. This section adopts the operational reading, and (20) states it.

Shell.calc_kC() does not always return \(\boldsymbol{K}_{C0} + \boldsymbol{K}_{CNL}\). If c_cte is given or any of Nxx_cte, Nyy_cte, Nxy_cte is non-zero, a geometric block is added into the matrix returned as kC, either \(\boldsymbol{K}_G(\boldsymbol{N}^{cte})\) through fkG0 or \(\boldsymbol{K}_G\) evaluated at c_cte through fkG_num. This is the superposition device of combined load cases, which keeps a part of the stress state out of reach of the load multiplier \(\lambda\). With NLgeom=True and follower pressure loads, their load stiffness Shell.calc_kCfollower() is added as well. The symbol \(\boldsymbol{K}_{C0} + \boldsymbol{K}_{CNL}\) therefore describes the return value of Shell.calc_kC() only when no constant pre-stress and no follower load is declared. The two modules of Verification that verify the Jacobian set the constant stresses to zero and do not exercise this branch; the third module does, but only to check that it is reached safely, not the value of the block.

Note on symbols

Three groups of symbols are introduced here because the implementation distinguishes quantities that the classical treatment does not:

• \(\boldsymbol{K}_{C0}\) and \(\boldsymbol{K}_{CNL}\), for the constant and the state-dependent parts of the matrix returned by fkC_num. The classical symbol for the first of them is \(\boldsymbol{K}_0\), also used here, in (16); \(\boldsymbol{K}_{C0}\) is kept in (3) because the three-way split it belongs to is the one the implementation is organised around.

• \(\boldsymbol{N}^0\), \(\boldsymbol{N}_L\) and \(\boldsymbol{N}_{NL}\), for the three contributions to the membrane stress entering (14). The decomposition is needed here because the implementation evaluates the third one separately, and places it on the other side of the \(\boldsymbol{K}_C\)–\(\boldsymbol{K}_G\) boundary.

• \(\boldsymbol{L}_p\), for the Boolean selection operator of the multi-domain assembly of (23), which has no counterpart in the single-domain treatment.

Product of the tangent stiffness matrix by a vector

Shell.calc_kTv() and MultiDomain.calc_kTv() return the product \([K_T(c)] \{v\}\) of the tangent stiffness matrix of Shell.calc_kT() by a vector, without assembling \([K_T]\), for all the geometrically non-linear models. Its cost is that of one evaluation of the internal force vector, instead of the assembly of \([K_T]\), whose cost grows with the square of the number of terms of the approximation. It is the exact product, to round-off, and the exact directional derivative of Shell.calc_fint().

The product is what the prescribed displacements imposed exactly need at each configuration, \(\partial R_r/\partial \lambda = [T]^T (\{F_{ext}\} - [K_T(c)] \{c_{p1}\})\), see Prescribed displacements imposed exactly, and what matrix-free iterative solvers and directional checks of the tangent need.

Derivation

All the geometrically non-linear models, the plates and the Donnell and Sanders cylinders of the CLPT, FSDT and TSDT, have at an integration point the generalized strains

\[\{\varepsilon(c)\} = [B(\beta/2)] \{c\}, \qquad \{\beta\} = [G] \{c\}\]

where \(\{\beta\}\) are the two rotations of the von Kármán terms, \(w_{,x}\) and \(w_{,y}\) for the Donnell kinematics, \(w_{,y} - v/r\) replacing \(w_{,y}\) for Sanders, and \([B(\beta)]\) is the strain operator, linear in \(\{\beta\}\), which enters only its membrane rows:

\[\begin{split}\frac{\partial [B]}{\partial \beta_x} = \begin{bmatrix} [G_x] \\ 0 \\ [G_y] \\ 0 \end{bmatrix}, \qquad \frac{\partial [B]}{\partial \beta_y} = \begin{bmatrix} 0 \\ [G_y] \\ [G_x] \\ 0 \end{bmatrix}\end{split}\]

Half the rotations in \([B]\) give the quadratic terms \(\beta_x^2/2\), \(\beta_y^2/2\) and \(\beta_x \beta_y\). The internal force vector is \(\{F_{int}\} = \int [B(\beta)]^T [F] \{\varepsilon(c)\} dA\), and its directional derivative along \(\{v\}\) is

(25)\[[K_T] \{v\} = \frac{d}{dt} \{F_{int}(c + t v)\} \Big|_{t=0} = \int \left( [B(\beta)]^T [F] [B(\beta)] \{v\} + [G]^T [N] [G] \{v\} \right) dA\]

because the strains are a quadratic form whose bilinear form is \([B(\beta)]\), such that \(\frac{d}{dt} \{\varepsilon(c + t v)\} = [B(\beta)] \{v\}\), and \(\frac{d}{dt} [B]^T \{\sigma\} = [G]^T [N] [G] \{v\}\) from the derivatives of \([B]\) above, with \([N]\) the membrane stresses. The two terms are those assembled by fkC_num and fkG_num of each model. The constant membrane stresses Nxx_cte, Nyy_cte and Nxy_cte, whose geometric stiffness Shell.calc_kC() adds to the tangent stiffness matrix, are added to \([N]\).

At each integration point, with the approximation functions tabulated once, the kernels calc_kTv of panels/models/<model>_num.pyx compute

1. \(\{\beta\} = [G] \{c\}\) and \([G] \{v\}\);

2. \(\{\varepsilon\} = [B(\beta/2)] \{c\}\) and \(\{\delta \varepsilon\} = [B(\beta)] \{v\}\);

3. the stresses \([F] \{\varepsilon\}\), of which \([N]\), and \([F] \{\delta \varepsilon\}\);

4. the sum of (25) over the terms.

Each step is a loop over the \(m n\) terms, against a loop over the \(m^2 n^2\) pairs of terms for the assembly of \([K_T]\).

The expressions are generated by theory/shells/tangent_product, which also verifies with SymPy, for every model, that (25) is the derivative of the internal force vector of one term, with a symbolic constitutive matrix; \([B]\) and \([G]\) being linear in the degrees of freedom, the identity extends to any number of terms.

Follower pressure loads

With follower pressure loads, Shell.calc_kTv() adds the product of their load stiffness Shell.calc_kCfollower() as the central difference

\[-\frac{\partial \{F_p\}}{\partial \{c\}} \{v\} = -\frac{\{F_p(c + t v)\} - \{F_p(c - t v)\}}{2 t}\]

which is exact for any \(t\), the force vector \(\{F_p(c)\}\) being affine in \(\{c\}\) for follower='linear' and quadratic for follower='quadratic', see Follower pressure loads: theory. \(t = \max(\|c\|/\|v\|, 1)\) avoids the cancellation of large and small terms.

Verification

The checks of the tangent stiffness matrix of tests/tests_shell/test_tangent_consistency.py are repeated for the product, for all 9 models, in tests/tests_shell/test_tangent_product.py:

• the product equals calc_kT(c) @ v with a relative difference below \(10^{-12}\), also with a constant membrane stress state, one constitutive matrix per integration point, a partial integration domain, an offset in a larger vector, and dead and follower pressure loads at a load factor inc;

• the directional Taylor test of Shell.calc_fint(), whose remainder falls as \(O(h^2)\);

• the matrix built column by column from the product against the central difference of Shell.calc_fint();

• in the undeformed state the product is \([K_0] \{v\}\);

• Newton-Raphson with the tangent built from the product converges quadratically:

@pytest.mark.parametrize('model', MODELS)
def test_newton_raphson_converges_quadratically(model):
    """Newton-Raphson around a known equilibrium, with the tangent built from
    the product, must square the error; see the test of the same name in
    ``test_tangent_consistency.py``"""
    s, fint, kTv, KT = make_callables(model)

    def KT_from_product(c):
        N = c.shape[0]
        K = np.empty((N, N))
        for j in range(N):
            e = np.zeros(N)
            e[j] = 1.
            K[:, j] = kTv(c, e)
        return K

    for seed in range(13, 33):
        rng = np.random.default_rng(seed)
        c_eq = random_state(s, rng)
        used = ~np.isclose(KT(c_eq).diagonal(), 0)
        c_eq[~used] = 0
        K0 = s.calc_kC(c=np.zeros_like(c_eq)).toarray()[np.ix_(used, used)]
        lam = eigh(KT(c_eq)[np.ix_(used, used)], K0, eigvals_only=True)
        if np.abs(lam).min() > 0.05:
            break
    else:
        raise AssertionError('no regular equilibrium state found')
    fext = fint(c_eq)

    c = c_eq + 1.e-2*random_state(s, rng)
    c[~used] = 0
    errors = [np.linalg.norm(c - c_eq)/np.linalg.norm(c_eq)]
    while errors[-1] > 1.e-8 and len(errors) < 12:
        K = KT_from_product(c)[np.ix_(used, used)]
        R = fint(c) - fext
        c[used] -= np.linalg.solve(K, R[used])
        errors.append(np.linalg.norm(c - c_eq)/np.linalg.norm(c_eq))

    assert len(errors) - 1 <= 6, 'too many iterations: %s' % errors
    floor = 1.e-9 if modelDB.db[model]['dofs'] == 5 else 1.e-12
    orders = [np.log(e2/e1)/np.log(e1/e0)
              for e0, e1, e2 in zip(errors[:-2], errors[1:-1], errors[2:])
              if e2 > floor]
    assert len(orders) >= 1, 'no step above the round-off floor: %s' % errors
    assert min(orders) > 1.5, 'convergence is not quadratic: %s' % errors

tests/multidomain/test_tangent_product_multidomain.py checks MultiDomain.calc_kTv() with penalty connections and follower loads, and \(\partial R_r/\partial \lambda\) of the prescribed displacements against the assembled tangent stiffness matrix.

Performance

Single shell, laminate [0, 45, -45, 90], nx = 2m, ny = 2n integration points, one core; Shell.calc_kT() with the tabulated kernels, see Approximation functions tabulated at the integration points:

model	\(m = n\)	size	calc_kT	calc_fint	calc_kTv	ratio	difference
plate_clpt_donnell	8 / 12 / 16	192 / 432 / 768	0.035 / 0.35 / 1.9 s	1.5 / 8.1 / 28 ms	1.8 / 9.6 / 32 ms	20 / 37 / 59	\(2 \times 10^{-15}\)
cylshell_fsdt_sanders	8 / 12 / 16	320 / 720 / 1280	0.16 / 1.7 / 9.1 s	2.6 / 14 / 47 ms	2.9 / 15 / 51 ms	54 / 109 / 180	\(2 \times 10^{-15}\)
plate_tsdt_donnell	8 / 12 / 16	320 / 720 / 1280	0.28 / 3.0 / 17 s	2.8 / 15 / 50 ms	3.1 / 16 / 55 ms	91 / 184 / 304	\(2 \times 10^{-15}\)

The product costs about one Shell.calc_fint(), and the ratio grows with the number of terms. With the previous kernels, whose approximation functions were evaluated in the innermost loops, calc_kT took 7.0, 10.9 and 11.9 s for \(m = n = 12\), i.e. 700 times the product.

Displacement-controlled post-buckling of the plate of tests/multidomain/test_null_space_pd.py with \(m = n = 12\) (407 reduced degrees of freedom), prescribed displacements imposed exactly, the four combinations of panels before (0.12.0, with the previous kernels and the assembly of \([K_T]\) for \(\partial R_r/\partial \lambda\)) and after these developments with structsolve 0.7.0 and 0.8.0. The reactions at \(\lambda = 1\) agree to all the digits printed:

panels	structsolve	Newton-Raphson	modified Newton-Raphson	Riks
before	0.7.0	285 s	705 s	363 s
before	0.8.0	229 s	158 s	350 s
after	0.7.0	11.9 s	9.1 s	18.7 s
after	0.8.0	11.8 s	7.8 s	18.4 s

Before, almost all the time was in the assembly of \([K_T]\) inside calc_fext, for \(\partial R_r/\partial \lambda\); after, it is in the assembly of \([K_T]\) for the Newton iterations, 32, 20 and 50 assemblies, the products costing 0.2 to 0.5 s in total.