序：我要寫一期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計算結果完全一致。

至此，桿單元結束，但在桿單元有限元分析中，仍有些沒提到的問題值得重視，例如：桿單元的方向問題、桿單元的個數對求解位移與應力的影響問題等。

有什么問題歡迎與我交流，微信聯系方式：