第二篇：2d杆单元

小仔_0671 2021年4月13日浏览：1717 收藏：5

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

1d杆单元类似于上篇提过的弹簧单元，这里不再写了，直接上手2d杆单元。杆单元刚度矩阵推导过程中需要遵守以下假定：

1）杆不能承受剪力或弯矩

2）忽略横向位移影响

3）遵守胡克定律

4）杆中间无外载

例：图示桁架与弹簧系统，已知桁架E=210GPa，A=0.0005m²，具体载荷如下图所示，求：a）总刚；b）节点位移；c）单元应力。

第二篇：2d杆单元的图1

一、有限元法求解

步骤1：离散化

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

步骤2：写单刚

单元1

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

第二篇：2d杆单元的图2

求解上式，得
第二篇：2d杆单元的图3

同理，
第二篇：2d杆单元的图4

第二篇：2d杆单元的图5

步骤3：写总刚

第二篇：2d杆单元的图6

步骤4：边界条件

第二篇：2d杆单元的图7

本例中，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)