2
$\begingroup$

Solver.py:

import numpy as np
import math 

E=3000000
A=2.0
truss = [
    np.array([[0, 0], [0, 8]]),
    np.array([[0, 0], [8, 8]]),
    np.array([[0, 8], [8, 8]])
]
lengths=[]
for i, t in enumerate(truss):
    lengths.append(math.hypot(t[0][0] - t[1][0], t[0][1] - t[1][1]))

n_joints = len(truss) + 1  # assuming the first point is also a joint
n_reactions = 0

k_w=[]
k_w.append((E * A / lengths[0]) * np.array([
    [1, 0, -1, 0],
    [0, 0, 0, 0],
    [-1, 0, 1, 0],
    [0, 0, 0, 0]
]))

k_w.append((E * A / lengths[1]) * np.array([
    [1, 0, -1, 0],
    [0, 0, 0, 0],
    [-1, 0, 1, 0],
    [0, 0, 0, 0]
]))

k_w.append((E * A / lengths[2]) * np.array([
    [1, 0, -1, 0],
    [0, 0, 0, 0],
    [-1, 0, 1, 0],
    [0, 0, 0, 0]
]))

# Initialize the global stiffness matrix
K = np.zeros((4*len(truss), 4*len(truss)))
theta = []

for i in range(len(truss)):
    dy = truss[i][0][0] - truss[i][1][0]
    dx = truss[i][0][1] - truss[i][1][1]
    theta.append(math.atan2(dy, dx))

# Assemble the global stiffness matrix
for i, t in enumerate(truss):
    L = np.array([
        [1, 0, 0, 0],
        [0, 1, 0, 0],
        [0, 0, math.cos(theta[i]), math.sin(theta[i])],
        [0, 0, -math.sin(theta[i]), math.cos(theta[i])]
    ])

    # Apply the transformation to the local stiffness matrix
    k_local = np.dot(np.dot(L.T, k_w[i]), L)

    # Add the transformed local stiffness matrix to the global stiffness matrix
    K[i * 4:(i + 1) * 4, i * 4:(i + 1) * 4] += k_local

    # Check if the current truss member has a support reaction
    if i == 0 or i == 4:
        n_reactions += 2

P = np.zeros(4*len(truss))
P[0]=1
P[1]=1
# Set boundary conditions
K[2, :] = 0
K[3, :] = 0
K[:, 2] = 0
K[:, 3] = 0
K[2, 2] = 1
K[3, 3] = 1
P[2] = 0
P[3] = 0

# Add rows and columns for reaction forces
#K = np.vstack((np.hstack((K, np.zeros((n_reactions, 4*len(truss))))), np.zeros((4*len(truss), n_reactions))))
#P = np.hstack((P, np.zeros(n_reactions)))

# Impose equilibrium
#for i in range(n_reactions):
#    row = 4*len(truss) + i
#    K[row, row] = 1
#    if i < 2:
#        K[row, i*4+1] = -1
#    else:
#        K[row, i*4+3] = -1

# Solve for the displacements
#displacements = np.linalg.solve(K, P)

# Add rows and columns for reaction forces
#n_reactions = 2  # assuming there are two reaction forces
#K = np.vstack((np.hstack((K, np.zeros((n_reactions, 4*len(truss))))), np.zeros((4*len(truss), n_reactions))))
#P = np.hstack((P, np.zeros(n_reactions)))

# Impose equilibrium
#for i in range(n_reactions):
#    row = 4*len(truss) + i
#    K[row, row] = 1
#    if i < 2:
#        K[row, i*4+1] = -1
#    else:
#        K[row, i*4+3] = -1

# Solve for the displacements
displacements = np.linalg.solve(K, P)

if(len(truss)+n_reactions<(n_joints*2)): print('statically determinant...')

So I am attempting so far to assign Dirichlet boundary conditions to nodes at 0,0 and 0,8 as fixed or zero reaction. As far as assembling K I think it's generally on the right track but I don't feel sure about it yet. Mostly at this point what I am trying to solve is the displacements whatever they start out as so far. So far the use of solve leads to numpy.linalg.LinAlgError: Singular matrix. and I am really not quite sure yet why that is....As far as the commented out bit about reactions and equilibrium I could maybe use a bit of help understading that too. Mostly I don'y understand yet why the singular matrix?

$\endgroup$

1 Answer 1

4
$\begingroup$

It looks as if you are making entire rows and columns of the stiffness matrix zero. This makes the matrix singular. If the solution's value is fixed at some entry, then remove that row/column from the system and substitute its value wherever it appears in the other equations and add it to the righthand side.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.