Donnell-type and Sanders-type kinematics#
The tests of the Sanders-type element are those of the Donnell-type element
with the imports of bfsccylinder.sanders, see Donnell-type and Sanders-type kinematics. The
linear buckling analysis of a laminated cylinder under a constant axial
compression of Linear buckling analysis, with the Sanders-type element:
from bfsccylinder.sanders import (BFSCCylinderSanders, update_KC0,
update_KG_constant_stress, DOF, DOUBLE, INT,
KC0_SPARSE_SIZE, KG_SPARSE_SIZE)
from bfsccylinder.quadrature import get_points_weights
from bfsccylinder.utils import assign_constant_ABD
def test_linear_buckling(plot=False):
# geometry Z33 Castro 2014
L = 0.510 # m
R = 0.250 # m
b = 2*pi*R # m
# number of nodes
ny = 50 # circumferential
nx = int(ny*L/b)
# material properties Geier 1997
E11 = 145.5e9
E22 = 8.7e9
nu12 = 0.28
G12 = 5.1e9
plyt = 0.125e-3
laminaprop = (E11, E22, nu12, G12, G12, G12)
stack = [0, 0, 19, -19, 37, -37, 45, -45, 51, -51]
prop = laminated_plate(stack=stack, plyt=plyt, laminaprop=laminaprop)
nids = 1 + np.arange(nx*(ny+1))
nids_mesh = nids.reshape(nx, ny+1)
# closing the cylinder by reassigning last row of node-ids
nids_mesh[:, -1] = nids_mesh[:, 0]
nids = np.unique(nids_mesh)
nid_pos = dict(zip(nids, np.arange(len(nids))))
xlin = np.linspace(0, L, nx)
ytmp = np.linspace(0, b, ny+1)
ylin = np.linspace(0, b-(ytmp[-1] - ytmp[-2]), ny)
xmesh, ymesh = np.meshgrid(xlin, ylin)
xmesh = xmesh.T
ymesh = ymesh.T
# getting nodes
ncoords = np.vstack((xmesh.flatten(), ymesh.flatten())).T
x = ncoords[:, 0]
y = ncoords[:, 1]
n1s = nids_mesh[:-1, :-1].flatten()
n2s = nids_mesh[1:, :-1].flatten()
n3s = nids_mesh[1:, 1:].flatten()
n4s = nids_mesh[:-1, 1:].flatten()
nint = 4
points, weights = get_points_weights(nint=nint)
num_elements = len(n1s)
print('num_elements', num_elements)
elements = []
N = DOF*nx*ny
print('N', N)
init_k_KC0 = 0
init_k_KG = 0
for n1, n2, n3, n4 in zip(n1s, n2s, n3s, n4s):
shell = BFSCCylinderSanders(nint)
shell.n1 = n1
shell.n2 = n2
shell.n3 = n3
shell.n4 = n4
shell.c1 = DOF*nid_pos[n1]
shell.c2 = DOF*nid_pos[n2]
shell.c3 = DOF*nid_pos[n3]
shell.c4 = DOF*nid_pos[n4]
shell.R = R
shell.lex = L/(nx-1)
shell.ley = b/ny
assign_constant_ABD(shell, prop)
shell.init_k_KC0 = init_k_KC0
shell.init_k_KG = init_k_KG
init_k_KC0 += KC0_SPARSE_SIZE
init_k_KG += KG_SPARSE_SIZE
elements.append(shell)
Kr = np.zeros(KC0_SPARSE_SIZE*num_elements, dtype=INT)
Kc = np.zeros(KC0_SPARSE_SIZE*num_elements, dtype=INT)
Kv = np.zeros(KC0_SPARSE_SIZE*num_elements, dtype=DOUBLE)
for shell in elements:
update_KC0(shell, points, weights, Kr, Kc, Kv)
KC0 = coo_matrix((Kv, (Kr, Kc)), shape=(N, N)).tocsc()
Nxx = -1
Nyy = 0
Nxy = 0
KGr = np.zeros(KG_SPARSE_SIZE*num_elements, dtype=INT)
KGc = np.zeros(KG_SPARSE_SIZE*num_elements, dtype=INT)
KGv = np.zeros(KG_SPARSE_SIZE*num_elements, dtype=DOUBLE)
for shell in elements:
update_KG_constant_stress(shell, points, weights, KGr, KGc, KGv, Nxx, Nyy, Nxy)
KG = coo_matrix((KGv, (KGr, KGc)), shape=(N, N)).tocsc()
# applying boundary conditions
bk = np.zeros(N, dtype=bool)
checkBC = isclose(x, 0) | isclose(x, L)
bk[0::DOF] = checkBC
#bk[1::DOF] = checkBC
#bk[2::DOF] = checkBC
bk[3::DOF] = checkBC
#bk[5::DOF] = checkBC
bk[6::DOF] = checkBC
bk[7::DOF] = checkBC
#bk[8::DOF] = checkBC
#bk[9::DOF] = checkBC
bu = ~bk # same as np.logical_not, defining unknown DOFs
Kuu = KC0[bu, :][:, bu]
KGuu = KG[bu, :][:, bu]
print('structural matrices OK')
# A * x[i] = lambda[i] * M * x[i]
num_eigvals = 3
Nu = N - bk.sum()
if True:
#NOTE this works and seems to be the fastest option
PREC = 1/Kuu.diagonal().max()
PREC2 = spilu(PREC*Kuu, diag_pivot_thresh=0, drop_tol=1e-8,
fill_factor=50)
print('spilu OK')
def matvec(x):
return PREC2.solve(x)
Ainv = LinearOperator(matvec=matvec, shape=(Nu, Nu))
maxiter = 1000
X = np.random.default_rng(0).random((Nu, num_eigvals)) - 0.5
X /= np.linalg.norm(X, axis=0)
#NOTE default tolerance is too large
tol = 1e-6
eigvals, eigvecsu, hist = lobpcg(A=PREC*Kuu, B=-PREC*KGuu, X=X, M=Ainv, largest=False,
maxiter=maxiter, retResidualNormsHistory=True, tol=tol)
assert len(hist) <= maxiter, 'did not converge'
load_mult = eigvals
else:
if True:
#NOTE works, but slower than lobpcg
eigvals, eigvecsu = eigsh(A=Kuu, k=num_eigvals, which='SM', M=KGuu,
tol=1e-9, sigma=1., mode='buckling')
load_mult = -eigvals
else:
#NOTE this is giving close but varying results for each run
PREC2 = spilu(PREC*Kuu, diag_pivot_thresh=0)
print('spilu OK')
def matvec(x):
return PREC2.solve(x)
Minv = LinearOperator(matvec=matvec, shape=(Nu, Nu))
eigvals, eigvecsu = eigsh(A=-PREC*KGuu, k=num_eigvals, which='LM',
M=PREC*Kuu, Minv=Minv)
load_mult = 1./eigvals
print('linear buckling analysis OK')
Pcr = load_mult[0]*Nxx*2*np.pi*R
print('load_mult[0]', load_mult[0])
print('Pcr top=', Pcr)
assert np.isclose(Pcr, -215444.50718870337, rtol=0.01)
mode = 0
mode_shape = np.zeros(N, dtype=float)
mode_shape[bu] = eigvecsu[:, mode]
if plot:
w = mode_shape[6::DOF].reshape(nx, ny)
import matplotlib
matplotlib.use('TkAgg')
import matplotlib.pyplot as plt
from matplotlib import cm
plt.gca().set_aspect('equal')
plt.contourf(xmesh, ymesh, w, levels=200, cmap=cm.jet)
plt.colorbar()
plt.show()
The reference critical load of the Sanders-type element, -215444.5 N, is
about 0.7 % smaller in magnitude than the -217000 N of the Donnell-type
element in tests/test_linear_buckling_Z33_constant_stress.py.
The ring of Follower pressure is solved with both kinematics by
the same functions, which take the module and the class from the dictionary
KINEMATICS of tests/test_follower_pressure_validation.py:
KINEMATICS = {
'donnell': (bfsccylinder, 'BFSCCylinder'),
'sanders': (bfsccylinder.sanders, 'BFSCCylinderSanders'),
}