Nonlinear static analysis#
A rectangular isotropic plate, with \(w\) constrained at the nodes of the
four edges, compressed along \(x\) below its critical load, with a transverse force at the centre large enough to reach a
deflection of about 2.6 times the thickness. The model keeps the value
arrays of \(\boldsymbol{K}_{CNL} + \boldsymbol{K}_{GNL}\) and \(\boldsymbol{K}_G\), which are zeroed and recomputed at each
iteration with update_KCNL() and update_KG(), and the
internal forces of update_fint():
def build_model():
points, weights = get_points_weights(nint=4)
lam = isotropic_plate(thickness=h, E=E, nu=nu)
xlin = np.linspace(0, a, nx)
ylin = np.linspace(0, b, ny)
xmesh, ymesh = np.meshgrid(xlin, ylin)
# getting nodes
ncoords = np.vstack((xmesh.T.flatten(), ymesh.T.flatten())).T
x = ncoords[:, 0]
y = ncoords[:, 1]
nids = 1 + np.arange(ncoords.shape[0])
nid_pos = dict(zip(nids, np.arange(len(nids))))
nids_mesh = nids.reshape(nx, ny)
n1s = nids_mesh[:-1, :-1].flatten()
n2s = nids_mesh[1:, :-1].flatten()
n3s = nids_mesh[1:, 1:].flatten()
n4s = nids_mesh[:-1, 1:].flatten()
num_elements = len(n1s)
print('num_elements', num_elements)
N = DOF*nx*ny
KC0r = np.zeros(KC0_SPARSE_SIZE*num_elements, dtype=INT)
KC0c = np.zeros(KC0_SPARSE_SIZE*num_elements, dtype=INT)
KC0v = np.zeros(KC0_SPARSE_SIZE*num_elements, dtype=DOUBLE)
KCNLr = np.zeros(KCNL_SPARSE_SIZE*num_elements, dtype=INT)
KCNLc = np.zeros(KCNL_SPARSE_SIZE*num_elements, dtype=INT)
KCNLv = np.zeros(KCNL_SPARSE_SIZE*num_elements, dtype=DOUBLE)
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)
init_k_KC0 = 0
init_k_KCNL = 0
init_k_KG = 0
elements = []
for n1, n2, n3, n4 in zip(n1s, n2s, n3s, n4s):
plate = BFSCPlate2D()
plate.n1 = n1
plate.n2 = n2
plate.n3 = n3
plate.n4 = n4
plate.c1 = DOF*nid_pos[n1]
plate.c2 = DOF*nid_pos[n2]
plate.c3 = DOF*nid_pos[n3]
plate.c4 = DOF*nid_pos[n4]
plate.ABD = lam.ABD
plate.lex = a/(nx - 1)
plate.ley = b/(ny - 1)
plate.init_k_KC0 = init_k_KC0
plate.init_k_KCNL = init_k_KCNL
plate.init_k_KG = init_k_KG
update_KC0(plate, points, weights, KC0r, KC0c, KC0v)
init_k_KC0 += KC0_SPARSE_SIZE
init_k_KCNL += KCNL_SPARSE_SIZE
init_k_KG += KG_SPARSE_SIZE
elements.append(plate)
KC0 = coo_matrix((KC0v, (KC0r, KC0c)), shape=(N, N)).tocsc()
# applying boundary conditions
# simply supported
bk = np.zeros(N, dtype=bool)
checkSS = isclose(x, 0) | isclose(x, a)
bk[3::DOF] = checkSS
bk[6::DOF] = checkSS
checkSS = isclose(y, 0) | isclose(y, b)
bk[6::DOF] += checkSS
check = isclose(x, a/2.)
bk[0::DOF] = check
bu = ~bk
print('boundary conditions bu', bu.sum())
m = Model()
m.points, m.weights, m.lam = points, weights, lam
m.ncoords, m.x, m.y = ncoords, x, y
m.xmesh, m.ymesh = xmesh, ymesh
m.nid_pos = nid_pos
m.elements = elements
m.N = N
m.KC0 = KC0
m.KC0uu = KC0[bu, :][:, bu]
m.bk, m.bu = bk, bu
m.KCNLr, m.KCNLc, m.KCNLv = KCNLr, KCNLc, KCNLv
m.KGr, m.KGc, m.KGv = KGr, KGc, KGv
return m
def calc_KCNL(m, u):
"""The nonlinear part KCNL(u) of the tangent, K_CNL + K_GNL."""
m.KCNLv *= 0
for plate in m.elements:
update_KCNL(u, plate, m.points, m.weights, m.KCNLr, m.KCNLc, m.KCNLv)
return coo_matrix((m.KCNLv, (m.KCNLr, m.KCNLc)),
shape=(m.N, m.N)).tocsc()
def calc_KG(m, u):
"""The geometric stiffness KG(u), of the tangent and of the linearised
comparisons."""
m.KGv *= 0
for plate in m.elements:
update_KG(u, plate, m.points, m.weights, m.KGr, m.KGc, m.KGv)
return coo_matrix((m.KGv, (m.KGr, m.KGc)), shape=(m.N, m.N)).tocsc()
def calc_fint(m, u):
fint = np.zeros(m.N)
for plate in m.elements:
update_fint(u, plate, m.points, m.weights, fint)
return fint
The full Newton-Raphson method of structsolve, with the tangent stiffness
matrix \(\boldsymbol{K}_{C0} + \boldsymbol{K}_{CNL} + \boldsymbol{K}_{GNL} + \boldsymbol{K}_G\) rebuilt at every iteration and the
external forces applied in four increments, newton_raphson(m, fext + ftrans,
increments=4), see Nonlinear static analysis with structsolve. The test checks the
quadratic convergence and that the membrane stretching makes the plate
stiffer than the linearised prediction \((\boldsymbol{K}_{C0} + \boldsymbol{K}_G)^{-1}\boldsymbol{F}_{ext}\), with
\(\boldsymbol{K}_G\) of the axial load:
def test_large_deflection(plot=False):
"""Transverse load reaching w/h of about 2.6: membrane stiffening."""
load = 300. # N
transverse = 5. # N
m = build_model()
fext, ftrans = forces(m, load, transverse)
load_cr, u_axial = linear_buckling_load(m, load)
# the linearised answer the plate would give if the von Karman terms did
# not stiffen it as it deflects
KG = calc_KG(m, u_axial)
u_lin = np.zeros(m.N)
u_lin[m.bu] = spsolve((m.KC0 + KG)[m.bu, :][:, m.bu], ftrans[m.bu])
w_lin = u_lin[6::DOF].max()
u, history = newton_raphson(m, fext + ftrans, increments=4)
assert_quadratic_convergence(history)
w = u[6::DOF].reshape(nx, ny).T
print('w min max', w.min(), w.max())
print('w/h', w.max()/h)
print('w linearised', w_lin)
if plot:
import matplotlib
matplotlib.use('TkAgg')
import matplotlib.pyplot as plt
plt.gca().set_aspect('equal')
levels = np.linspace(w.min(), w.max(), 300)
plt.contourf(m.xmesh, m.ymesh, w, levels=levels)
plt.colorbar()
plt.show()
# the regime is genuinely nonlinear
assert 2. < w.max()/h < 3.
# and membrane stretching makes the plate stiffer than the linearised
# prediction, by a wide margin at this deflection
assert w.max() < 0.6*w_lin
assert isclose(w.max(), 2.5826494310546447e-03, rtol=1e-3)