空间梁刚度矩阵计算——有限元基础

1.总体流程

（1）计算单元刚度矩阵；

（2）对单元刚度矩阵进行坐标变换，获得全局坐标系下的单元刚度矩阵：

（3）组装全局刚度矩阵；

2.关键代码

2.1 单元刚度矩阵计算

def beam_element_stiffness(E,G, A,J,Iy,Iz,L):
    """梁单元的刚度矩阵"""
    k_local=np.zeros((12, 12))

    k_local[0, 0] = A*E/L

    k_local[1, 1]=12*E*Iz/L**3

    k_local[2,2]=12*E*Iy/L**3

    k_local[3,3]=G*J/L

    k_local[4,2]=-6*E*Iy/L**2
    k_local[2,4]=k_local[4,2]
    k_local[4,4]=4*E*Iy/L

    k_local[5,1]=6*E*Iz/L**2
    k_local[1,5]=k_local[5,1]
    k_local[5,5]=4*E*Iz/L

    k_local[6,0]=-A*E/L
    k_local[0,6]=k_local[6,0]
    k_local[6,6]=A*E/L

    k_local[7,1]=-12*E*Iz/L**3
    k_local[1,7]=k_local[7,1]
    k_local[7,5]=-6*E*Iz/L**2
    k_local[5,7]=k_local[7,5]
    k_local[7,7]=12*E*Iz/L**3

    k_local[8,2]=-12*E*Iy/L**3
    k_local[2,8]=k_local[8,2]
    k_local[8,4]=6*E*Iy/L**2
    k_local[4,8]=k_local[8,4]
    k_local[8,8]=12*E*Iy/L**3

    k_local[9,3]=-G*J/L
    k_local[3,9]=k_local[9,3]
    k_local[9,9]=G*J/L

    k_local[10,2]=-6*E*Iy/L**2
    k_local[2,10]=k_local[10,2]
    k_local[10,4]=2*E*Iy/L
    k_local[4,10]=k_local[10,4]
    k_local[10,8]=6*E*Iy/L**2
    k_local[8,10]=k_local[10,8]
    k_local[10,10]=4*E*Iy/L

    k_local[11,1]=6*E*Iz/L**2
    k_local[1,11]=k_local[11,1]
    k_local[11,5]=2*E*Iz/L
    k_local[5,11]=k_local[11,5]
    k_local[11,7]=-6*E*Iz/L**2
    k_local[7,11]=k_local[11,7]
    k_local[11,11]=4*E*Iz/L

    return k_local

2.2 坐标变换矩阵计算

首先计算lambda矩阵

def lambda_matrix(node1Coord, node2Coord):
    L = np.linalg.norm(node2Coord - node1Coord)
    x1, y1, z1 = node1Coord
    x2, y2, z2 = node2Coord
    l = (x2 - x1) / L
    m = (y2 - y1) / L
    n = (z2 - z1) / L

    lambdaMatrix = np.zeros((3, 3))
    if l == 0 and m == 0 and n > 0:
        lambdaMatrix[0, 2] = 1
        lambdaMatrix[1, 1] = 1
        lambdaMatrix[2, 0] = -1
    elif l == 0 and m == 0 and n < 0:
        lambdaMatrix[0, 2] = -1
        lambdaMatrix[1, 1] = 1
        lambdaMatrix[2, 0] = 1
    else:
        D = np.sqrt(l ** 2 + m ** 2)
        lambdaMatrix[0, 0] = l
        lambdaMatrix[0, 1] = m
        lambdaMatrix[0, 2] = n
        lambdaMatrix[1, 0] = -m / D
        lambdaMatrix[1, 1] = l / D
        lambdaMatrix[1, 2] = 0
        lambdaMatrix[2, 0] = -l * n / D
        lambdaMatrix[2, 1] = -m * n / D
        lambdaMatrix[2, 2] = D

    return lambdaMatrix

使用lambda矩阵组装坐标变换矩阵：

def beam_element_transformation(node1Coord, node2Coord):
    lambdaMaxtix = lambda_matrix(node1Coord, node2Coord)
    transformMatrix=np.zeros((12,12))
    transformMatrix[0:3,0:3]=lambdaMaxtix
    transformMatrix[3:6,3:6]=lambdaMaxtix
    transformMatrix[6:9,6:9]=lambdaMaxtix
    transformMatrix[9:12,9:12]=lambdaMaxtix
    return transformMatrix

2.3 全局刚度矩阵计算

def build_global_stiffness_matrix(node_positions, element_nodes, E, G, A, I, J):
    num_nodes = len(node_positions)
    K_global = np.zeros((6*num_nodes, 6*num_nodes))
    for row in element_nodes:
        id1, id2 = row
        # 取梁单元的两个节点
        node1_pos = node_positions[id1]
        node2_pos = node_positions[id2]
        element_length = np.linalg.norm(node2_pos - node1_pos)

        # 计算局部刚度矩阵
        k_local = beam_element_stiffness(E, G, A, I, I, J, element_length)

        # 计算空间梁变换矩阵
        transformMatrix = beam_element_transformation(node1_pos, node2_pos)

        # 计算全局刚度矩阵
        k_global = np.dot(np.dot(transformMatrix.T, k_local), transformMatrix)

        # 将局部刚度矩阵加到全局刚度矩阵中
        K_global[6*id1:6*id1+6, 6*id1:6*id1+6] += k_global[0:6, 0:6]
        K_global[6*id1:6*id1+6, 6*id2:6*id2+6] += k_global[0:6, 6:12]
        K_global[6*id2:6*id2+6, 6*id1:6*id1+6] += k_global[6:12, 0:6]
        K_global[6*id2:6*id2+6, 6*id2:6*id2+6] += k_global[6:12, 6:12]
    
    return K_global

参考《有限元方法基础教程》第五版。