有限元三角形单元Python实现

有限元三角形单元

1.网格划分

2.三角形单元的刚度矩阵

可以采用下列公式计算：

其中，t是厚度，其它：

平面应力问题：

平面应变问题：

根据传送矩阵得到总刚度

2.荷载向量

F

3.后处理法处理边界

主1副0法

4.求解结果

KD=F

例题求解：（参考的：https://www.bilibili.com/read/cv24507049/）

弹性模量设置为20，泊松比为0.2。以2点为坐标原点，1号单元长直角边为x轴，短直角边为y轴，约束1，2点x向和y向位移。荷载数值为1，作用于5号节点，方向竖直向下。

Python代码实现：

import numpy as np

E, u, t = 20, 0.1, 1  # E，u，t 为E，泊松比，厚度
F_nodes = np.array([[5,2,-1]])     # 节点力，节点号，力类型，力大小
fixed_nodes = [0,1]
ele_num = 3          # 单元数量
node_num = 5         # 节点数量
N = 2*node_num       # 总自由度（未引入约束）
**# 划分单元，网格**
node_xy = np.array([[0,1],
                    [0,0],
                    [2,0],
                    [2,1],
                    [4,1]])    # 节点坐标
ele_node = np.array([[0,1,2],
                    [0,2,3],
                    [2,4,3]])       # 单元的节点号
ele_node = np.array(ele_node,dtype=int)     # 转为整数
# print(ele_node)
**# 求总刚度矩阵**
def total_K(ele_node,node_xy,E,u,t):
    total_k = np.zeros((N,N))              # 总的刚度矩阵
    for ele_index in range(ele_num):      # 遍历每一个单元
        n1,n2,n3 = ele_node[ele_index]    # 单元的节点，并求单元坐标
        ele_xy = np.vstack([node_xy[n1],node_xy[n2],node_xy[n3]])
        ele_k = ele_K(ele_xy, E, u, t)      # 带入后面函数求单元矩阵
        ele_s = ele_S(ele_node[ele_index])  # 带入后面函数求传送矩阵
        total_k = total_k + np.dot(np.dot(ele_s, ele_k), ele_s.T)
    return total_k

#**求单元刚度矩阵**

def ele_K(ele_xy,E,u,t):
    # 三角形单元刚度求解 ele_xy类似于np.array([[1,1],[2,2],[3,3]])
    # E，u，t 为参数弹性模量，泊松比，厚度
    C = np.c_[np.array([1, 1, 1]), ele_xy]  # C矩阵
    A = np.linalg.det(C)/2      # 三角形面积
    C_inv = np.linalg.inv(C)
    B = np.array([[C_inv[1, 0], 0, C_inv[1, 1], 0, C_inv[1, 2], 0],
                  [0, C_inv[2, 0], 0, C_inv[2, 1], 0, C_inv[2, 2]],
                  [C_inv[2, 0], C_inv[1, 0], C_inv[2, 1], C_inv[1, 1], C_inv[2, 2], C_inv[1, 2]]])
    D = E / (1 - u * u) * np.array([[1, u, 0],
                                    [u, 1, 0],
                                    [0, 0, (1 - u) / 2]])
    ele_k =  np.dot(np.dot(B.T,D),B)*t*A
    return ele_k
# 求单元的转化矩阵
def ele_S(ele_3p):
    ele_s = np.zeros((N,6))
    for i in range(3):
        ele_s[2* ele_3p[i], 2*i] = 1
        ele_s[2*ele_3p[i]+1, 2*i+1] = 1
return ele_s

def F_node(F_nodes,N):
    # 根据力列表得到节点力
    F = np.zeros(N)
    for f_index in range(len(F_nodes)):
        node = int(F_nodes[f_index,0])-1
        type = int(F_nodes[f_index,1])-1
        F[2*node+type] = F_nodes[f_index,2]
    return F

F = F_node(F_nodes,N)
total_k  = total_K(ele_node,node_xy,E,u,t)

for item in fixed_nodes:
    # 5个节点
    for i in range(2):
        node_fix = item * 2 + i
        total_k[node_fix] = 0
        total_k[:, node_fix] = 0
        total_k[node_fix, node_fix] = 1
        F[node_fix] = 0

D = np.linalg.solve(total_k,F)
print("displacement:\n",D)

结果为：

符合结果。

Q4-1转为三角形单元计算代码：

import numpy as np

E, u, t, a, b = 1, 0.1,0.01, 2, 3  # E，u，t 为E，泊松比，厚度,a,b
nx,ny = 10,30                        # 网格划分nx,ny 单元数目
F_nodes = np.array([[200,2,50]])    # 节点力，节点号，力类型，力大小
fixed_nodes = [0, 1, 2, 3, 4]
ele_num = nx * ny * 2           # 单元数量
node_num = (nx + 1) * (ny + 1)  # 节点数量
N = 2*node_num                  # 总自由度（未引入  约束）

# 划分单元，网格
node_xy = np.zeros((node_num,2))    # 节点坐标
ele_node = np.zeros((ele_num,3))    # 单元的节点号
w = a / nx      # 计算x向划分的长度
h = b / ny      # 计算y向划分的长度
# 开始定义节点
for j in range(ny + 1):         # 分行
    for i in range(nx + 1):     # 分列
        node_xy[j*(nx+1) + i] = np.array([w*i,h*j]) # 节点坐标
couter = 0                  # 单元计数，划分单元，并记录单元的节点
for j in range(ny):         # 分层
    for i in range(nx):     # 分列
        ele_node[couter] = np.array([j*(nx+1)+i,j*(nx+1)+i+1,(j+1)*(nx+1)+i])
        ele_node[couter+1] = np.array([j*(nx+1)+i+1,(j+1)*(nx+1)+i+1,(j+1)*(nx+1)+i])
        couter += 2

ele_node = np.array(ele_node,dtype=int)     # 转为整数

# 求总刚度矩阵
def total_K(ele_node,node_xy,E,u,t):
    total_k = np.zeros((N,N))              # 总的刚度矩阵
    for ele_index in range(ele_num):      # 遍历每一个单元
        n1,n2,n3 = ele_node[ele_index]    # 单元的节点，并求单元坐标
        ele_xy = np.vstack([node_xy[n1],node_xy[n2],node_xy[n3]])
        ele_k = ele_K(ele_xy, E, u, t)      # 带入后面函数求单元矩阵
        ele_s = ele_S(ele_node[ele_index])  # 带入后面函数求传送矩阵
        total_k = total_k + np.dot(np.dot(ele_s, ele_k), ele_s.T)
    return total_k

# 求单元刚度矩阵
def ele_K(ele_xy,E,u,t):
    # 三角形单元刚度求解 ele_xy类似于np.array([[1,1],[2,2],[3,3]])
    # E，u，t 为参数弹性模量，泊松比，厚度
    C = np.c_[np.array([1, 1, 1]), ele_xy]  # C矩阵
    A = abs(np.linalg.det(C)/2)      # 三角形面积
    C_inv = np.linalg.inv(C)
    B = np.array([[C_inv[1, 0], 0, C_inv[1, 1], 0, C_inv[1, 2], 0],
                  [0, C_inv[2, 0], 0, C_inv[2, 1], 0, C_inv[2, 2]],
                  [C_inv[2, 0], C_inv[1, 0], C_inv[2, 1], C_inv[1, 1], C_inv[2, 2], C_inv[1, 2]]])
    D = E / (1 - u * u) * np.array([[1, u, 0],
                                    [u, 1, 0],
                                    [0, 0, (1 - u) / 2]])
    ele_k =  np.dot(np.dot(B.T,D),B)*t*A
    return ele_k

# 求单元的转化矩阵
def ele_S(ele_3p):
    ele_s = np.zeros((N,6))
    for i in range(3):
        ele_s[2* ele_3p[i], 2*i] = 1
        ele_s[2*ele_3p[i]+1, 2*i+1] = 1
    return ele_s

def F_node(F_nodes,zyd):
    # 根据力列表得到节点力
    F = np.zeros(zyd)
    for f_index in range(len(F_nodes)):
        node = int(F_nodes[f_index,0])-1
        type = int(F_nodes[f_index,1])-1
        F[2*node+type] = F_nodes[f_index,2]
    return F

F = F_node(F_nodes,N)
total_k = total_K(ele_node,node_xy,E,u,t)
# print(total_k)
# 后处理约束，主1副0法
for item in fixed_nodes:
    for i in range(2):
        node_fix = item * 2 + i
        total_k[node_fix] = 0
        total_k[:, node_fix] = 0
        total_k[node_fix, node_fix] = 1
        F[node_fix] = 0

D = np.linalg.solve(total_k,F)
print("displacement:\n",D)

结果如下，比较接近，认为合理（考虑误差可能来自单元不同导致的误差，或者节点荷载可能对应节点号加减了一位，误差也可能来自单元精度不够）。
自己的三角形：

老师的四边形：