使用python和直接刚度矩阵法计算节点位移

Question

我使用python和直接刚度矩阵方法分析了以下框架。

我使用以下坐标系定义：

我得到的位移值是正确的，但是节点位移的方向在x方向上是负的。这是不正确的。我已经通过其他 FEM 程序检查了结果。 x 方向的位移应为正。

你能看看我下面的代码吗？也许你会看错

import numpy as np

e = 3
x_1 = 0
y_1 = 180

x_2 = 90
y_2 = 180

x_3 = 180
y_3 = 180

x_4 = 90+e
y_4 = 90

x_5 = 90
y_5 = 0
# ======================================
E_beam = 210000 # Young's modulus (MPa)
A_beam = 120  # cross-sectional area (mm^2)
I_beam = 4000  # moment of inertia (mm^4)

E_column = 210000 # Young's modulus (MPa)
A_column = 120  # cross-sectional area (mm^2)
I_column = 360  # moment of inertia (mm^4)
# ======================================

# Initialize the global stiffness matrix
N = 5  # total number of nodes
dofs_per_node = 3  # 3 degrees of freedom per node (u, v, and theta)
K_glob = np.zeros((N * dofs_per_node, N * dofs_per_node))

# ================================== element stiffness matrix 1
L_1 = np.sqrt((x_2 - x_1) ** 2 + (y_2 - y_1) ** 2)  # calculate the length of the element
cos_theta_1 = (x_2 - x_1) / L_1
sin_theta_1 = (y_2 - y_1) / L_1
# Arguments in element stiffness matrix
a_1 = E_beam * A_beam / L_1
b_1 = 12 * E_beam * I_beam / L_1 ** 3
c_1 = 6 * E_beam * I_beam / L_1 ** 2
d_1 = 4 * E_beam * I_beam / L_1
e_1 = 2 * E_beam * I_beam / L_1

k_e_beam_1 = np.array([[a_1, 0, 0, -a_1, 0, 0],
                 [0, b_1, c_1, 0, -b_1, c_1],
                 [0, c_1, d_1, 0, -c_1, e_1],
                 [-a_1, 0, 0, a_1, 0, 0],
                 [0, -b_1, -c_1, 0, b_1, -c_1],
                 [0, c_1, e_1, 0, -c_1, d_1]])
# Compute the transformation matrix
T_beam_1 = np.array([[cos_theta_1, sin_theta_1, 0, 0, 0, 0],
                 [-sin_theta_1, cos_theta_1, 0, 0, 0, 0],
                 [0, 0, 1, 0, 0, 0],
                 [0, 0, 0, cos_theta_1, sin_theta_1, 0],
                 [0, 0, 0, -sin_theta_1, cos_theta_1, 0],
                 [0, 0, 0, 0, 0, 1]])

k_glob_beam_1 = np.matmul(np.matmul(T_beam_1, k_e_beam_1), np.transpose(T_beam_1))

# ================================== element stiffness matrix 2
L_2 = np.sqrt((x_3 - x_2) ** 2 + (y_3 - y_2) ** 2)  # calculate the length of the element
cos_theta_2 = (x_3 - x_2) / L_2
sin_theta_2 = (y_3 - y_2) / L_2

# Arguments in element stiffness matrix
a_2 = E_beam * A_beam / L_2
b_2 = 12 * E_beam * I_beam / L_2 ** 3
c_2 = 6 * E_beam * I_beam / L_2 ** 2
d_2 = 4 * E_beam * I_beam / L_2
e_2 = 2 * E_beam * I_beam / L_2

k_e_beam_2 = np.array([[a_2, 0, 0, -a_2, 0, 0],
                 [0, b_2, c_2, 0, -b_2, c_2],
                 [0, c_2, d_2, 0, -c_2, e_2],
                 [-a_2, 0, 0, a_2, 0, 0],
                 [0, -b_2, -c_2, 0, b_2, -c_2],
                 [0, c_2, e_2, 0, -c_2, d_2]])

# Compute the transformation matrix
T_beam_2 = np.array([[cos_theta_2, sin_theta_2, 0, 0, 0, 0],
                 [-sin_theta_2, cos_theta_2, 0, 0, 0, 0],
                 [0, 0, 1, 0, 0, 0],
                 [0, 0, 0, cos_theta_2, sin_theta_2, 0],
                 [0, 0, 0, -sin_theta_2, cos_theta_2, 0],
                 [0, 0, 0, 0, 0, 1]])

k_glob_beam_2 = np.matmul(np.matmul(T_beam_2, k_e_beam_2), np.transpose(T_beam_2))

# ================================== element stiffness matrix 3
L_3 = np.sqrt((x_4 - x_2) ** 2 + (y_4 - y_2) ** 2)  # calculate the length of the element
cos_theta_3 = (x_4 - x_2) / L_3
sin_theta_3 = (y_4 - y_2) / L_3

# Arguments in element stiffness matrix
a_3 = E_column * A_column / L_3
b_3 = 12 * E_column * I_column / L_3 ** 3
c_3 = 6 * E_column * I_column / L_3 ** 2
d_3 = 4 * E_column * I_column / L_3
e_3 = 2 * E_column * I_column / L_3

k_e_column_3 = np.array([[a_3, 0, 0, -a_3, 0, 0],
                 [0, b_3, c_3, 0, -b_3, c_3],
                 [0, c_3, d_3, 0, -c_3, e_3],
                 [-a_3, 0, 0, a_3, 0, 0],
                 [0, -b_3, -c_3, 0, b_3, -c_3],
                 [0, c_3, e_3, 0, -c_3, d_3]])

# Compute the transformation matrix
T_column_3 = np.array([[cos_theta_3, sin_theta_3, 0, 0, 0, 0],
                 [-sin_theta_3, cos_theta_3, 0, 0, 0, 0],
                 [0, 0, 1, 0, 0, 0],
                 [0, 0, 0, cos_theta_3, sin_theta_3, 0],
                 [0, 0, 0, -sin_theta_3, cos_theta_3, 0],
                 [0, 0, 0, 0, 0, 1]])

k_glob_column_3 = np.matmul(np.matmul(T_column_3, k_e_column_3), np.transpose(T_column_3))

# ================================== element stiffness matrix 4
L_4 = np.sqrt((x_5 - x_4) ** 2 + (y_5 - y_4) ** 2)  # calculate the length of the element
cos_theta_4 = (x_5 - x_4) / L_4
sin_theta_4 = (y_5 - y_4) / L_4
print(cos_theta_4, sin_theta_4)

# Arguments in element stiffness matrix
a_4 = E_column * A_column / L_4
b_4 = 12 * E_column * I_column / L_4 ** 3
c_4 = 6 * E_column * I_column / L_4 ** 2
d_4 = 4 * E_column * I_column / L_4
e_4 = 2 * E_column * I_column / L_4

k_e_column_4 = np.array([[a_4, 0, 0, -a_4, 0, 0],
                 [0, b_4, c_4, 0, -b_4, c_4],
                 [0, c_4, d_4, 0, -c_4, e_4],
                 [-a_4, 0, 0, a_4, 0, 0],
                 [0, -b_4, -c_4, 0, b_4, -c_4],
                 [0, c_4, e_4, 0, -c_4, d_4]])

# Compute the transformation matrix
T_column_4 = np.array([[cos_theta_4, sin_theta_4, 0, 0, 0, 0],
                 [-sin_theta_4, cos_theta_4, 0, 0, 0, 0],
                 [0, 0, 1, 0, 0, 0],
                 [0, 0, 0, cos_theta_4, sin_theta_4, 0],
                 [0, 0, 0, -sin_theta_4, cos_theta_4, 0],
                 [0, 0, 0, 0, 0, 1]])

k_glob_column_4 = np.matmul(np.matmul(T_column_4, k_e_column_4), np.transpose(T_column_4))

K_glob[0:6, 0:6] += k_glob_beam_1
K_glob[3:9, 3:9] += k_glob_beam_2
# ===================
K_glob[3:6, 3:6] += k_glob_column_3[0:3, 0:3]
K_glob[9:12, 3:6] += k_glob_column_3[3:6, 0:3]
K_glob[3:6, 9:12] += k_glob_column_3[0:3, 3:6]
K_glob[9:12, 9:12] += k_glob_column_3[3:6, 3:6]
K_glob[9:15, 9:15] += k_glob_column_4

# External Loads
P_0 = -10000  # (N)

# Create the zero force vector F_global
F_global = np.zeros((N * dofs_per_node, 1))

# Specify the position to add P_0 value to F_global vector
i = 4

# Add the value to x at the specified position
F_global[i, 0] += P_0

# ====================================== Applying of BC
# BC for Global stiffness matrix
K_glob_BC = K_glob
K_glob_BC = np.delete(K_glob_BC, [0, 1, 7, -2, -3], 0)  # Delete rows
K_glob_BC = np.delete(K_glob_BC, [0, 1, 7, -2, -3], 1)  # Delete columns

# BC for Global Force Vector
F_vec_glob_BC = F_global
F_vec_glob_BC = np.delete(F_vec_glob_BC, [0, 1, 7, -2, -3], 0)  # Delete rows

# ====================================== Calculation of displacement (Linear solution)
U_glob_BC = np.matmul(np.linalg.inv(K_glob_BC), F_vec_glob_BC)

# add zero values on the position of BC
U_glob = np.insert(np.asarray(U_glob_BC), (0, 0, -1, -1, -5), values=[0], axis=0)
U_glob_resh = np.reshape(U_glob, (5, 3))
print(U_glob_resh)

结果：

[[ 0.00000000e+00  0.00000000e+00 -1.79379767e-03]
[-4.07387190e-04 -9.66284063e-02  3.66648471e-04]
[-4.07387190e-04  0.00000000e+00  1.42714920e-03]
[-4.49222495e-01 -4.83345801e-02  2.03579898e-03]
[ 0.00000000e+00  0.00000000e+00 -8.52345411e-03]]