Matrix Analysis Problem

[feelthismoment18]

Hello all,

I'm currently workingon some code and am trying to figure something out.

I'm using the stiffness method in an iterative solver to simulate a displacement controlled test on a structure. I am rasing the 3rd node in the z direction.

I have the nodes and elements in the following format:

'''
alpha = np.radians(30)

beta = np.radians(60)

a = 1000

E = 210000 #N/mm2

A = 100 #mm2

nodes = np.array([

[0, 0, 0],

[a*np.sin(0.5*np.pi-beta)/np.sin(0.5*np.pi), -a*np.sin(beta)/np.sin(0.5*np.pi), 0],

[a*np.sin(np.pi-beta-alpha)/np.sin(alpha), 0, 0],

[a*np.sin(np.pi-2*beta-alpha)/np.sin(beta+alpha), 0, 0]])

elements = np.array([

[0, 1],

[1, 2],

[0, 2],

[1, 3],

[2, 3]])

'''

My problem is that because all the nodes start on the same plane, the matrix is singular and cannot be used to solve an F = KU relationship in a 3D problem because essentially its a 2D problem at the start.

I cannot just turn it into a 2D problem because I'm assessing the vertical dispalcement/force reltionship.

I tried to start one of the nodes at a tiny fraction higher than the rest of them just to not get a singular matrix but it creates a very strange stiffness matrix that produces force and displacement results way off of what I would be expecting.

Has anyone got any advice for how to deal with this. Also I've attached the small amount of code below I've done for this so far if anyone wants to see it.

import numpy as np
import sympy as sp
import matplotlib.pyplot as plt

alpha = np.radians(30)
beta = np.radians(60)
a = 1000

E = 210000 #N/mm2
A = 100 #mm2

nodes = np.array([
[0, 0, 1],
[a*np.sin(0.5*np.pi-beta)/np.sin(0.5*np.pi), -a*np.sin(beta)/np.sin(0.5*np.pi), 1],
[a*np.sin(np.pi-beta-alpha)/np.sin(alpha), 0, 1+1**(-100)],
[a*np.sin(np.pi-2*beta-alpha)/np.sin(beta+alpha), 0, 1]])

elements = np.array([
[0, 1],
[1, 2],
[0, 2],
[1, 3],
[2, 3]])

U = np.zeros((3*nodes.shape[0], 1))
nodes_0 = nodes.copy()
lengths = np.linalg.norm(nodes[elements[:,1]-1] - nodes[elements[:,0]-1], axis=1)
lengths_0, cos_x, cos_y, cos_z = structure(nodes_0, elements, U)

dof = np.array([0, 1, 2, 5, 7, 10, 11])
restrained = np.array([3, 4, 6, 9])

U_max = 1000
ninc = 100
inc = U_max/ninc

F = outer_force_vector(nodes, restrained)
U_inc = outer_disp_vector(nodes, dof)
U_inc[8] = inc
F_unit = outer_force_vector(nodes, restrained)
F_unit[8] = 1
U = outer_disp_vector(nodes, dof)
N = np.zeros((elements.shape[0], 1))

K_global_list = []

for i, element in enumerate(elements):

K_global = K_global_element(cos_x, cos_y, cos_z, N, lengths, lengths_0, E, A)

K_global_list.append(K_global)

K = assemble_K(K_global_list, nodes, elements)

sp.pprint(K)

equation = K @ U - F_unit
print(U)
print(F_unit)
unknowns = (U.free_symbols).union(F_unit.free_symbols)
solution = sp.solve(equation,*unknowns)
print(solution)

load_ratio = inc/solution["Uz3"]

equation = K @ U_inc - F
unknowns = (U_inc.free_symbols).union(F.free_symbols)
solution = sp.solve(equation,*unknowns)
print(solution)

# Definitions in different cell block


# Define Original Geometric Properties

def structure(nodes_0, elements, U):

U = U.reshape(nodes_0.shape[0], 3)

nodes = nodes_0 + U

lengths = np.linalg.norm(nodes[elements[:,1]-1] - nodes[elements[:,0]-1], axis=1)

cos_x = []
cos_y = []
cos_z = []

for i, element in enumerate(elements):

node1, node2 = nodes_0[elements[i,0]-1], nodes_0[elements[i,1]-1]
cx = (np.array(node2) - np.array(node1))[0]/lengths
cy = (np.array(node2) - np.array(node1))[1]/lengths
cz = (np.array(node2) - np.array(node1))[2]/lengths

cos_x.append(cx)
cos_y.append(cy)
cos_z.append(cz)

lengths = np.array(lengths).reshape(elements.shape[0], 1)
cos_x = np.array(cos_x).reshape(elements.shape[0], 1)
cos_y = np.array(cos_y).reshape(elements.shape[0], 1)
cos_z = np.array(cos_z).reshape(elements.shape[0], 1)

return lengths, cos_x, cos_y, cos_z


# Displacement Vector (outer-loop)

def outer_disp_vector(nodes, dof):

U_symbols = []

for i in range(1, nodes.shape[0]+1):

U_symbols.append(sp.Symbol(f'Ux{i}'))
U_symbols.append(sp.Symbol(f'Uy{i}'))
U_symbols.append(sp.Symbol(f'Uz{i}'))

U = sp.Matrix(U_symbols)

for i in dof:

U = 0

return U


# Displacement Vector (outer-loop)

def outer_force_vector(nodes, restrained):

F_symbols = []

for i in range(1, nodes.shape[0]+1):
F_symbols.append(sp.Symbol(f'Fx{i}'))
F_symbols.append(sp.Symbol(f'Fy{i}'))
F_symbols.append(sp.Symbol(f'Fz{i}'))

F = sp.Matrix(F_symbols)

for i in restrained:

F = 0

return F


# Calculate Stiffness Matrix for Each Element

def K_global_element(cx, cy, cz, N, L, L0, E, A):

K_M = (A * E / L0) * np.array([
[cx*cx, cx*cy, cx*cz, -cx*cx, -cx*cy, -cx*cz],
[cx*cy, cy*cy, cy*cz, -cx*cy, -cy*cy, -cy*cz],
[cx*cz, cy*cz, cz*cz, -cx*cz, -cy*cz, -cz*cz],
[-cx*cx, -cx*cy, -cx*cz, cx*cx, cx*cy, cx*cz],
[-cx*cy, -cy*cy, -cy*cz, cx*cy, cy*cy, cy*cz],
[-cx*cz, -cy*cz, -cz*cz, cx*cz, cy*cz, cz*cz]])

K = K_M

return K


# Assemble Global Stiffness for Entire Structure

def assemble_K(K_global_list, nodes, elements):

K = np.zeros((3*nodes.shape[0], 3*nodes.shape[0]))

for element_idx, element in enumerate(elements):

node1, node2 = element[0]-1, element[1]-1

dof = np.array([3*node1, 3*node1+1, 3*node1+2, 3*node2, 3*node2+1, 3*node2+2])

c = K_global_list[element_idx]

for i in range(6):
for j in range(6):
K[dof, dof[j]] += c[i,j].item()

return K

 

[rb1957]

your problem is that the nodes are co-planar ... so move one node slightly off the plane defined by he other three.

 

[feelthismoment18]

your problem is that the nodes are co-planar ... so move one node slightly off the plane defined by he other three.

I tried to do that as mentioned in the original post but it gave some crazy results

 

[ANE91]

Global stiffness matrix is singular until you remove the rows and columns associated with redundancies/reactions. The kinematic submatrix should be invertible.