序：我要寫一期python和Abaqus與有限元的文章，從彈簧單元、桿單元一直到實體單元，通過簡單的實例用python編程，Abaqus驗證結果。

1d桿單元類似于上篇提過的彈簧單元，這里不再寫了，直接上手2d桿單元。桿單元剛度矩陣推導過程中需要遵守以下假定：

1）桿不能承受剪力或彎矩

2）忽略橫向位移影響

3）遵守胡克定律

4）桿中間無外載

例：圖示桁架與彈簧系統，已知桁架E=210GPa，A=0.0005m²，具體載荷如下圖所示，求：a）總剛；b）節點位移；c）單元應力。

一、有限元法求解

步驟1：離散化

單元	節點i	節點j
1	1	2
2	1	3
3	1	4

步驟2：寫單剛

單元1

θ=135°，cosθ=-0.707，sinθ=0.707

求解上式，得

同理，

步驟3：寫總剛

步驟4：邊界條件

本例中，u2=v2=u3=v3=u4=v4=0，F1x=0，F1y=-1000N，代入上述方程

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

利用上述方程不難解出u1=-0.095mm，v1=-0.19mm，具體不再贅述。

步驟6：后處理，求單元1應力與單元2應力

取出相應的方程可求得σ1=2.8Mpa，σ2=-2Mpa。

二、python求解

import numpy as np
import math
# 離散化
dimen = 2  # 2d
x1 = 0
x2 = -5 * math.sqrt(2) / 2
x3 = -10
x4 = 0
y4 = -10
y1 = 0
y2 = 5 * math.sqrt(2) / 2
y3 = 0
k4 = 1000
ele_node = np.array([[1, 2], [1, 3], [1, 4]])
node_coord = np.array([[x1, y1], [x2, y2], [x3, y3], [x4, y4]])
# 計算單剛
E = 210 * 10 ** 9
A = 0.0005
L1 = math.sqrt((x2 - x1) ** 2 + (y2 - y1) ** 2)
L2 = math.sqrt((x3 - x1) ** 2 + (y3 - y1) ** 2)
def fun(C, S):
    Ke = np.array([[C * C, C * S, -C * C, -C * S], [C * S, S * S, -C * S, -S * S],
                   [-C * C, -C * S, C * C, C * S],
                   [-C * S, -S * S, C * S, S * S]])
    return Ke
def fun1(theta):
    x = theta * math.pi / 180
    return x
theta1 = 135
x1 = fun1(theta1)
C1 = round(math.cos(x1),15)
S1 = round(math.sin(x1),15)
theta2 = 180
x2 = fun1(theta2)
C2 = round(math.cos(x2),2)
S2 = round(math.sin(x2),2)
theta3 = 270
x3 = fun1(theta3)
C3 = round(math.cos(x3),2)
S3 = math.sin(x3)
T1 = fun(C1, S1)
T2 = fun(C2, S2)
T3 = fun(C3, S3)
k1 = E * A / L1 * T1
k2 = E * A / L2 * T2
k3 = k4 * 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([-C1,-S1,C1,S1])
C2=np.mat([-C2,-S2,C2,S2])
u = np.row_stack((u, [0],[0]))
stress1=E/L1*C1*u
stress2=E/L2*C2*u
print(stress1)
print(stress2)

計算結果如下圖

三、Abaqus求解

步驟1：離散化

步驟2：創建單元、施加邊界條件

步驟3：求解，后處理

結果：

1.整體剛度矩陣如下圖：

2.節點位移u1=-0.0952mm，v1=-0.1906mm，如下：

3.單元1的σ1=2.8Mpa，σ2=-2Mpa。，如下圖：

結論：

不難發現，手工計算、python計算、Abaqus計算結果完全一致。

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