第三篇：3d杆单元

序：我要写一期python和Abaqus与有限元的文章，从弹簧单元、杆单元一直到实体单元，通过简单的实例用python编程，Abaqus验证结果。

3d杆同其他维度杆单元类似，手工计算部分考虑到是比较麻烦的重复性工作，此处不再对手工计算结果进行对比。作为杆单元的最后一篇，在这里再提一下自由度，1d杆由2个节点组成，节点自由度为1，2d杆由2节点组成，自由度为2,3d杆由2个节点组成，自由度为3。杆单元只能承受轴向力作用，适于平面或空间桁架、网架仿真。

例：图示桁架系统，已知桁架E=210GPa，A=0.0005m²，其中，节点2,3,4全约束，节点1约束z向自由度。具体载荷以及坐标如下图所示，求：a）节点位移；b）单元应力。

一、有限元法求解

步骤1：离散化

单元	节点i	节点j
1	1	2
2	1	3
3	1	4

步骤2：写单刚

三维空间任意方向杆单元的矩阵公式为：

单元1

代入公式，求得各单元刚度矩阵即可。

步骤3：写总刚

同其他维度杆。

步骤4：边界条件

同其他维度杆。

步骤5：求方程，解u1和v1

同其他维度杆。

步骤6：后处理，求单元应力

三维杆单元应力求解公式：

代入公式求出应力即可。

二、python求解

import numpy as np

import math

# 离散化

dimen = 3  # 3d

x1 = 72

x2 = 0

x3 = 0

x4 = 0

y1 = 0

y2 = 0

y3 = 72

y4 = -48

z1=0

z2=-36

z3=-36

z4=0

ele_node = np.array([[1, 2], [1, 3], [1, 4]])

node_coord = np.array([[x1, y1,z1], [x2, y2,z2], [x3, y3,z3], [x4, y4,z4]])

# 计算单刚

E = 210 * 10 ** 9

A = 0.0005

L1 = math.sqrt((x2 - x1) ** 2 + (y2 - y1) ** 2+(z2 - z1) ** 2)

L2 = math.sqrt((x3 - x1) ** 2 + (y3 - y1) ** 2+(z3 - z1) ** 2)

L3 = math.sqrt((x4 - x1) ** 2 + (y4 - y1) ** 2+(z4 - z1) ** 2)

def fun(Cx, Cy,Cz):

    Ke = 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]])

    return Ke

Cx1 = (x2-x1)/L1

Cy1=(y2-y1)/L1

Cz1=(z2-z1)/L1

T1=fun(Cx1,Cy1,Cz1)

k1 = E * A / L1 * T1

Cx2 = (x3-x1)/L2

Cy2=(y3-y1)/L2

Cz2=(z3-z1)/L2

T2=fun(Cx2,Cy2,Cz2)

k2 = E * A / L2 * T2

Cx3 = (x4-x1)/L3

Cy3=(y4-y1)/L3

Cz3=(z4-z1)/L3

T3=fun(Cx3,Cy3,Cz3)

k3 = E * A / L3 * T3

# print(k1)

# print(k2)

# print(k3)

# 总自由度数

ndof = dimen * len(node_coord)

# 初始化总刚

K = np.zeros((ndof, ndof))

# # 计算总刚

def fun2(K, k, i, j):

    K[2 * i - 2, 2 * i - 2] = K[2 * i - 2, 2 * i - 2] + k[0, 0]

    K[2 * i - 2, 2 * i-1] = K[2 * i - 2, 2 * i-1] + k[0, 1]

    K[2 * i - 2, 2 * j - 2] = K[2 * i - 2, 2 * j - 2] + k[0, 2]

    K[2 * i - 2, 2 * j-1] = K[2 * i - 2, 2 * j-1] + k[0, 3]

    K[2 * i-1, 2 * i - 2] = K[2 * i-1, 2 * i - 2] + k[1, 0]

    K[2 * i-1, 2 * i-1] = K[2 * i-1, 2 * i-1] + k[1, 1]

    K[2 * i-1, 2 * j - 2] = K[2 * i-1, 2 * j - 2] + k[1, 2]

    K[2 * i-1, 2 * j-1] = K[2 * i-1, 2 * j-1] + k[1, 3]

    K[2 * j - 2, 2 * i - 2] = K[2 * j - 2, 2 * i - 2] + k[2, 0]

    K[2 * j - 2, 2 * i-1] = K[2 * j - 2, 2 * i-1] + k[2, 1]

    K[2 * j - 2, 2 * j - 2] = K[2 * j - 2, 2 * j - 2] + k[2, 2]

    K[2 * j - 2, 2 * j-1] = K[2 * j - 2, 2 * j-1] + k[2, 3]

    K[2 * j-1, 2 * i - 2] = K[2 * j-1, 2 * i - 2] + k[3, 0]

    K[2 * j-1, 2 * i-1] = K[2 * j-1, 2 * i-1] + k[3, 1]

    K[2 * j-1, 2 * j - 2] = K[2 * j-1, 2 * j - 2] + k[3, 2]

    K[2 * j-1, 2 * j-1] = K[2 * j-1, 2 * j-1] + k[3, 3]

    return (K)

K = fun2(K, k1, 1, 2)

K= fun2(K, k2, 1, 3)

K= fun2(K, k3, 1, 4)

# print(K)

# # 求节点位移

k=K[0:2, 0:2]

k=np.mat(k)

F=np.mat([0,-1000])

u=k.I*F.T

print(u)

# # 求单元应力

C1=np.mat([-Cx1,-Cy1,-Cz1,Cx1,Cy1,Cz1])

C2=np.mat([-Cx2,-Cy2,-Cz2,Cx2,Cy2,Cz2])

C3=np.mat([-Cx3,-Cy3,-Cz3,Cx3,Cy3,Cz3])

u = np.row_stack((u, [0],[0],[0],[0]))

stress1=E/L1*C1*u

stress2=E/L2*C2*u

stress3=E/L3*C3*u

print(stress1)

print(stress2)

print(stress3)

计算结果如下图

三、Abaqus求解

步骤1：离散化

步骤2：施加边界条件

步骤3：求解，后处理

结果：

1.节点位移u1=6.92×10^-5m，v1=-0.00125m，如下：

2.单元的应力：σ1=161466Pa，σ2=1.71×10⁶Pa，σ3=-1.55×10⁶Pa。，如下图：

结论：

不难发现，python、Abaqus计算结果完全一致。

至此，杆单元结束，但在杆单元有限元分析中，仍有些没提到的问题值得重视，例如：杆单元的方向问题、杆单元的个数对求解位移与应力的影响问题等。

有什么问题欢迎与我交流，微信联系方式：