Exact connections with the null-space method#

The domains of a MultiDomain assembly were originally connected with penalty stiffnesses, see MultiDomain.get_kC_conn() and [castro2017Multidomain], still available with conn_method='penalty'. The penalty method adds, for each quantity \(q_k\) that must be continuous across a connection, the energy:

(1)#\[U_k = \frac{1}{2} k_k \int \left( q_k^{(1)} - q_k^{(2)} \right)^2\]

where \(k_k\) is the translation or rotation penalty constant, 'kt' or 'kr', and the integral is along the connected edges, or over the area of the connection 'SB'. The continuity is then only approximate. Low penalty constants make the assembly too flexible, and high ones make the stiffness matrix ill-conditioned. The default constants of calc_kt_kr() are a compromise, and they can still lead to errors of a few percent, see Comparison with the penalty method.

With conn_method='null-space', the default, the connections are imposed exactly, as with penalty constants going to infinity, without any penalty constant:

md = MultiDomain(panels, conn)  # conn_method='null-space'

and the penalty method is selected with:

md = MultiDomain(panels, conn, conn_method='penalty')

The connections that can be imposed with the null-space method are given in nullspace.NULL_SPACE_FUNCS: 'SSxcte', 'SSycte', 'BFxcte', 'BFycte' and 'SB', for all the models, see Quantities imposed. The damaged connection 'SB_TSL' represents the physical stiffness of a cohesive interface, see Cohesive zone model of the multi-domain approach, and it keeps its penalty stiffness also with conn_method='null-space'.

Theory#

The quantities \(q_k\) of Eq. (1) are linear combinations of the Ritz constants \(\{c\}\) of the two domains. Since Bardell’s functions are polynomials, see Bardell’s Hierarchical Functions, the difference \(q_k^{(1)} - q_k^{(2)}\) is a polynomial along the connection, of the natural coordinate \(\xi\) or \(\eta\) shared by the two domains. It vanishes everywhere when it vanishes at as many points as the number of its coefficients, and the connections become the homogeneous linear constraints:

(2)#\[[B] \{c\} = \{0\}\]

where each row of \([B]\) is \(q_k^{(1)} - q_k^{(2)}\) at one Gauss-Legendre point along the connection, or over the area of 'SB', see nullspace.constraint_matrix(). The first four functions of Bardell are cubic and the function \(i \ge 4\) has the order \(i\), so that \(\max(m_1, m_2, 4)\) points are used along \(x\) and \(\max(n_1, n_2, 4)\) along \(y\). The zero energy of Eq. (1) and the constraints of Eq. (2) are equivalent: the vectors with zero penalty energy are exactly those that satisfy the constraints, which is verified for every connection and model in tests/multidomain/test_null_space.py.

Why not Lagrange multipliers#

The constraints could be imposed with Lagrange multipliers \(\{\lambda\}\), which add one unknown per constraint and give the system:

\[\begin{split}\begin{bmatrix} [K] & [B]^T \\ [B] & [0] \end{bmatrix} \begin{Bmatrix} \{c\} \\ \{\lambda\} \end{Bmatrix} = \begin{Bmatrix} \{f\} \\ \{0\} \end{Bmatrix}\end{split}\]

This matrix is indefinite. For the linear buckling and frequency analyses the geometric stiffness matrix and the mass matrix have zero rows and columns for the multipliers, giving infinite eigenvalues that the shift-and-invert eigenvalue solvers of structsolve do not handle. Instead of adding unknowns, the null-space method removes them, giving the same solution.

Elimination of the constraints#

All the vectors that satisfy Eq. (2) are written with a basis \([T]\) of the null space of \([B]\), \([B] [T] = [0]\):

(3)#\[\{c\} = [T] \{c_r\}\]

where \(\{c_r\}\) are the independent Ritz constants. Substituting Eq. (3) into the total potential energy gives the reduced system:

(4)#\[[K_r] = [T]^T [K] [T], \qquad [K_{G r}] = [T]^T [K_G] [T], \qquad [M_r] = [T]^T [M] [T], \qquad \{f_r\} = [T]^T \{f\}\]

where \([K]\), \([K_G]\), \([M]\) and \(\{f\}\) are the matrices and vectors of the domains, which do not contain the connections imposed exactly. The reduced stiffness matrix is symmetric and positive definite whenever the connected assembly has no rigid-body motion, like the stiffness matrix of the penalty method, and the static, linear buckling, frequency and non-linear solvers of structsolve are used unchanged. The constraints are linear in \(\{c\}\), therefore \([T]\) does not change during a non-linear analysis:

\[\{f_{int, r}\} = [T]^T \{f_{int}\}([T] \{c_r\}), \qquad [K_{T r}] = [T]^T [K_T]([T] \{c_r\}) [T]\]

\([T]\) is built by nullspace.null_space_basis(), one connection at a time:

the rows of \([B]\) of the connection are normalized to a unit norm and written in the independent Ritz constants of the connections already processed, \([B] [T]\);
the rank of the constraints is revealed with a QR decomposition with column pivoting, \([B] [T] [P] = [Q] [R]\), and each pivot larger than the tolerance tol gives one independent constraint;
the Ritz constants of the pivot columns are eliminated, \(\{c_e\} = -[R_{11}]^{-1} [R_{12}] \{c_k\}\), where \(\{c_k\}\) are the remaining constants, and \([T]\) is updated.

Redundant constraints are found in step 2 and dropped. They happen for instance at a corner shared by four domains, around a closed cylinder, or where a boundary condition flag already removes the term that a constraint would remove. Bardell’s functions of order 5 and higher vanish together with their slope at the ends of the domain, \(\xi = \pm 1\), so the constraints of an edge involve only the first four functions across it, and \([T]\) has very few non-zero terms: the elimination is similar to the assembly of the hierarchical elements of the p-version of the finite element method. A connection along an interior line, e.g. 'BFycte' at the middle of the base, involves all the functions across the line.

The terms removed by the boundary condition flags are zero rows and columns of the matrices of the domains, as with the penalty method, and remain in \(\{c_r\}\) with zero rows and columns in the reduced matrices.

\([T]\) is stored in MultiDomain.T and reused while the connections, the number of terms and the boundary condition flags do not change, see MultiDomain.get_T().

Quantities imposed#

The quantities imposed by each connection are those of the corresponding penalty kernels, with p1 and p2 the first and second panels of the connection. For 'SB', the two panels share the axes \(x, y, z\), \(z\) being normal to their mid-surfaces, p1 is the panel on the positive side of the interface along \(z\) and p2 the one on the negative side, e.g. a stiffener base on the side of positive \(z\) of a skin is p1 and the skin p2, see MultiDomain.get_kC_conn():

Connection	Classical laminated plate theory (3 DOFs)	Shear deformation theories (5 DOFs)
`'SSxcte'`	\(u, v, w, w_{,x}\), see `kCSSxcte.pyx`	\(u, v, w, \phi_x, \Phi_y\), and \(w_{,x}\) for the TSDT, see `fkCSSxcte_sdt()`
`'SSycte'`	\(u, v, w, w_{,y}\), see `kCSSycte.pyx`	\(u, v, w, \phi_x, \Phi_y\), and \(w_{,y}\) for the TSDT, see `fkCSSycte_sdt()`
`'BFycte'`	\(u_1 = u_2\), \(v_1 = w_2\), \(w_1 = -v_2\), \(w_{1,y} - v_1/r_1 = w_{2,y} - v_2/r_2\), see `panels.multidomain.connections.kCBFycte`	\(u_1 = u_2\), \(v_1 = w_2\), \(w_1 = -v_2\), \(\Phi_{y1} = \Phi_{y2}\)
`'BFxcte'`	\(u_1 = w_2\), \(v_1 = v_2\), \(w_1 = -u_2\), \(w_{1,x} = w_{2,x}\), see `panels.multidomain.connections.kCBFxcte`	\(u_1 = w_2\), \(v_1 = v_2\), \(w_1 = -u_2\), \(\phi_{x1} = \phi_{x2}\)
`'SB'`	\(u_1 + d_{sb} w_{1,x} = u_2\), \(v_1 + d_{sb} w_{1,y} = v_2\), \(w_1 = w_2\), with \(d_{sb} = h_1/2 + h_2/2\), see `kCSB.pyx`	\(u, v, w\) at the interface, \(z = -h_1/2\) of `p1` and \(z = h_2/2\) of `p2`, see `fkCSB_sdt()`, and \(\phi_x, \Phi_y\) only when the key `'kr'` of the connection is given and not zero

where \(\Phi_y = \phi_y + v/r\) with the Sanders-Koiter kinematics and \(\Phi_y = \phi_y\) otherwise, and the terms \(v/r\) of the classical laminated plate theory are present only with the Sanders-Koiter kinematics. The keys 'kt' and 'kr' of the connections are not used, apart from 'kr' of 'SB' with the shear deformation theories, which tells whether the rotations are connected.

Usage#

The matrices and vectors computed by MultiDomain have the size of all the Ritz constants and do not contain the connections imposed exactly: solving them directly would leave the domains disconnected. They are reduced with MultiDomain.reduce() before solving, and the solutions are expanded to all the Ritz constants with MultiDomain.expand(), which are used by the post-processing methods, e.g. MultiDomain.uvw(), MultiDomain.stress() and MultiDomain.plot(). With conn_method='penalty', \([T]\) is the identity matrix, so that the same code runs with both methods.

Linear static analysis:

from structsolve import solve

md = MultiDomain(panels, conn)
kC = md.calc_kC()
fext = md.calc_fext()
c = md.expand(solve(md.reduce(kC), md.reduce(fext)))

Linear buckling analysis with the pre-buckling state of a static analysis:

from structsolve import lb, static

kC = md.reduce(md.calc_kC())
incs, cs = static(kC, md.reduce(md.calc_fext()))
c = md.expand(cs[0])
kG = md.reduce(md.calc_kG(c=c))
eigvals, eigvecs = lb(kC, kG)
eigvecs = md.expand(eigvecs)

Frequency analysis:

from structsolve import freq

eigvals, eigvecs = freq(md.reduce(md.calc_kC()), md.reduce(md.calc_kM()))
eigvecs = md.expand(eigvecs)

Geometrically non-linear analysis, with the functions of the reduced problem given by MultiDomain.get_reduced_functions():

from structsolve import Analysis

an = Analysis(*md.get_reduced_functions())
an.static(NLgeom=True)
cs = [md.expand(c_r) for c_r in an.cs]

The functions create_cylinder(), create_cylinder_blade_stiffened(), the analyses of panels.multidomain.cylinder and panels.multidomain.cylinder_blade_stiffened, and tstiff2d_1stiff_freq(), tstiff2d_1stiff_compression() and tstiff2d_1stiff_flutter() take the argument conn_method.

Comparison with the penalty method#

The results below are verified in tests/multidomain/test_null_space.py.

The penalty method is more flexible than the exact connections, and its error decreases with the inverse of the penalty constants. For a thin plate divided in 2x2 domains, with the classical laminated plate theory, the error of the first frequencies is 1.3% with the default constants of calc_kt_kr(), and it becomes 0.13%, 0.013% and 0.0013% with 10, 100 and 1000 times higher constants. The condition number of the stiffness matrix grows with the penalty constants: for the 4-domain plate of tests/multidomain/test_conn_4panels_kt_kr.py, the reduced stiffness matrix of the null-space method has a condition number about 600 times lower than the stiffness matrix of the penalty method with the default constants, and \(6 \times 10^4\) times lower than with 100 times higher constants.

Case	Penalty, default constants	Penalty, higher constants	Null-space
Frequencies of a thick plate divided in 2x2 domains, FSDT, error with respect to the single domain (Noor 1973)	\(7.8 \times 10^{-2}\)	\(9.5 \times 10^{-5}\) (1000 times)	\(3 \times 10^{-12}\)
Laminates \([0/90]\) over \([0/90]\) with `'SB'`, error of the frequencies with respect to \([0/90]_2\), CLPT		\(8 \times 10^{-6}\) (\(10^5\) times)	\(6 \times 10^{-13}\)
Same, large deflection with the von Karman kinematics			\(10^{-15}\)
Blade-stiffened cylinder, buckling load with the pre-buckling state of a static analysis, \(N_{xx}\) (N/m)	-39993.8	-41694.9 (1000 times)	-41696.7
Closed cylinder, buckling load \(N_{xx}\) (N/m)	-47056.0		-47113.0
T-stiffened panel, first frequency (rad/s)	48.2733	48.3858 (100 times)	48.3875

The default penalty constants underestimate the buckling load of the blade-stiffened cylinder by 4%, whereas the higher constants converge to the result of the null-space method.

Limitations#

A connection imposed exactly constrains the domains exactly as the penalty method does with infinite penalty constants, which may constrain the approximation more than the penalty method with the default constants does:

When the two domains have a different number of terms along the connection, the terms of the trace of the richer domain that the other domain cannot represent are forced to vanish.
The boundary condition flags at the ends of an edge of one domain are imposed on the other domain.
With 'SB' and the classical laminated plate theory, \(u_1 = u_2 - d_{sb} w_{,x}\). Restraining the normal in-plane displacement of an edge in both domains, e.g. \(u = 0\) at \(x = 0\), also restrains the rotation \(w_{,x}\) at that edge, which is clamped. The in-plane displacements of the bonded laminates must be restrained consistently with the single laminate that they represent, e.g. restraining only the tangential displacement of the edges, as in test_skin_base_nonlinear_equals_single_laminate.

Only the connections of the domains are imposed exactly: the prescribed displacements of Shell.add_point_pd(), Shell.add_distr_pd_fixed_x() and Shell.add_distr_pd_fixed_y(), and the traction-separation law of 'SB_TSL', remain penalty stiffnesses.

The rank of the constraints is found with the tolerance tol of MultiDomain.get_T(), 1e-9 by default, which applies to the constraints normalized to a unit norm. The independent constraints have pivots many orders of magnitude above the round-off errors of the redundant ones, and a warning is issued when this gap is smaller than \(10^4\).

API#

The methods of MultiDomain for the null-space method are MultiDomain.get_T(), MultiDomain.reduce(), MultiDomain.expand() and MultiDomain.get_reduced_functions().