2D Truss单元计算实例编程[Python]

#桁架单元:桁架结构作用集中力荷载

import numpy as np
import matplotlib.pyplot as plt

#一、输入信息

#模尺材料参数
E=200000    #弹性模量
A=4532      #截面面积
EA=E*A      #抗拉刚度

#节点坐标及单元节点序号
NodeCoord=np.array([[-4500,0],[-1500,0],[1500,0],[4500,0],[-3000,3000],[0,3000],[3000,3000]])
EleNode=np.array([[1,2],[2,3],[3,4],[5,6],[6,7],[1,5],[2,6],[3,7],[2,5],[3,6],[4,7]])

#获得单元数量和单元节点数量
numEle = EleNode.shape[0]   # 11
numNode = NodeCoord.shape[0] # 7
#节点坐标值(x，y)列向量
xx= NodeCoord[:,0]
yy=NodeCoord[:,1]
#被约束自由度
ConstrainedDof = np.array([1,2,8])
#自由度总数量
numDOF = numNode*2  # 14
#整体位移向量定义
displacement = np.zeros((1, numDOF)) # 1*14 矩阵
#整体荷载向量定义
force = np.zeros((1, numDOF))
force[0][(10-1)]=-100000  # -1转换为索引
force[0][(12-1)]=-100000
force[0][(14-1)]=-100000


#二、计算分析

#定义总体刚度矩阵
#创建一个大小为 mumDOFx mumDOF的全零数组，用于存储整体刚度矩阵
stiffness_sum = np.zeros((numDOF,numDOF))

#遍历单元，求解单元刚度矩阵
for i in range(numEle):
    #求得单元刚度矩阵
    noindex = EleNode[i] #获取第i个单元的节点索引
    deltax = xx[noindex[1] - 1] - xx[noindex[0] - 1] #计算第i个单元在横向上的坐标差
    deltay = yy[noindex[1] - 1] - yy[noindex[0] - 1] #计算第i个单元在纵向上的坐标差
    L= np.sqrt(deltax *deltax + deltay * deltay) #计算第个单元的长度L
    C=deltax/L #计算角度的余弦值C
    S= deltay/L #计算角度的正弦值S
    eleK=(E*A/L)*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]
    ]) #计算单元刚度矩阵 eleK
    #获得节点对应自由度
    eleDof= np.array([noindex[0]*2-2,noindex[0]*2-1, noindex[1]*2-2,noindex[1]*2-1])

    # 将单元刚度矩阵组装到整体刚度矩阵
    for i in range(4):
        for j in range(4):
            stiffness_sum[int(eleDof[i]),int(eleDof[j])] = stiffness_sum[int(eleDof[i]),int(eleDof[j])] + eleK[i][j]

#根据给定的约束条件，删除指定被约束自由度对应的行和列，将整体刚度矩阵进行简化处理
stiffness_sum_sim = np.delete(stiffness_sum, ConstrainedDof-1,0)
stiffness_sum_sim = np.delete(stiffness_sum_sim, ConstrainedDof-1,1)


#解方程,求位移
activeDof = np.delete((np.arange(numDOF)+1), ConstrainedDof-1)
displacement[0][activeDof-1] = np.array(np.asmatrix(stiffness_sum_sim).I * np.asmatrix(force[0][activeDof-1]).T).T

#计算单元轴力
axialForce =np.zeros((numEle,1))
for i in range(numEle):
    noindex = EleNode[i]
    deltax = xx[noindex[1]-1]-xx[noindex[0]-1]
    deltay = yy[noindex[1]-1]-yy[noindex[0]-1]
    L = np.sqrt(deltax * deltax + deltay *deltay)
    C = deltax/L
    S = deltay/L
    eleDof = np.array([noindex[0]*2-1, noindex[0]*2,noindex[1]*2-1,noindex[1]*2])
    axialForce[i]=((EA/L)*np.asmatrix([-C,-S,C,S]))*np.asmatrix(displacement[0][eleDof - 1]).T

#计算支座反力
reation = np.zeros((len(ConstrainedDof),1))
for i in range (len(ConstrainedDof)):
    reation[i]= np.asmatrix(stiffness_sum[(ConstrainedDof[i]-1)]) * np.asmatrix(displacement).T

#三、计算用于绘图相关参数
#计算各单元最大长度maxlen
maxlen =-1
for i in range(numEle):
    noindex= EleNode[i]
    deltax = xx[noindex[1]-1]-xx[noindex[0]-1]
    deltay = yy[noindex[1]-1]-yy[noindex[0]-1]
    L= np.sqrt(deltax*deltax + deltay * deltay)
    if L>maxlen:
        maxlen =L

#计算节点最大绝对位移
maxabsDisp =0
for i in range(numNode):
    indexU=2*(i+1)-1
    indexV= 2*(i+1)
    dispU= displacement[0][indexU-1]
    dispV= displacement[0][indexV-1]
    Disp =np.sqrt(dispU*dispU + dispV* dispV)
    if Disp> maxabsDisp:
        maxabsDisp = Disp

#获取变形图变形放大系数
scalefactor =0.2
factor=1.1
if maxabsDisp> 1e-30:
    factor = scalefactor * maxlen/maxabsDisp

#四、绘制图形
#遍历单元，绘制变形前结构形状
for i in range(numEle):
    noindex =EleNode[i]
    x=xx[noindex[:]-1]
    y=yy[noindex[:]-1]
    plt.plot(x, y,color='0.6', marker='o', mfc='w', linestyle='--')

#遍历单元节点，标注单元节点编号
for i in range(numNode):
    plt.text(xx[i]+200,yy[i],i+1,color='r', ha='center', va='center')

#定义变形后节点坐标值(x，y)列向量
xx_Deformed= np.empty_like(xx)
yy_Deformed =np.empty_like(yy)

# 日遍单元，绘制变形后结构形状
for i in range(numEle):
    noindex= EleNode[i]
    indexU= 2*noindex[:]-1
    indexV=2*noindex[:]
    dispU= displacement[0][indexU-1]
    dispV= displacement[0][indexV-1]
    xxDeformed = xx[noindex[:]-1]+ dispU*factor
    yyDeformed = yy[noindex[:]-1] + dispV*factor
    xx_Deformed[noindex-1]= xxDeformed
    yy_Deformed[noindex-1]= yyDeformed
    plt.plot(xxDeformed, yyDeformed, color='black',marker='o', mfc='w',)

#设置坐标格式，使得图形以适合的比例进行显示
max_xx = max(max(xx),max(xx_Deformed))
max_yy = max(max(yy),max(yy_Deformed))
min_xx = min(min(xx),min(xx_Deformed))
min_yy = min(min(yy), min(yy_Deformed))
minx =-0.1*(max_xx-min_xx)+min_xx
maxx =0.1*(max_xx-min_xx)+max_xx
miny = -0.1*(max_yy-min_yy)+min_yy
maxy =0.1*(max_yy-min_yy)+max_yy
plt.xlim(minx,maxx)
plt.ylim(miny,maxy)
plt.show()

#绘制内力图
#遍历单元，绘制变形前结构形状
for i in range(numEle):
    noindex = EleNode[i]
    x=xx[noindex[:]-1]
    y=yy[noindex[:]-1]
    plt.plot(x, y, color='0.6', marker='o', mfc='w', linestyle='--')

#定义变形后节点坐标值(x，y)列向量
xx_Deformed= np.empty_like(xx)
yy_Deformed = np.empty_like(yy)
#将内力值归一化，仅绘图用
absminForce = 1e10
axialForceStd = np.zeros((numEle, 1))
for i in range(numEle):
    if abs(axialForce[i]) <absminForce:
        absminForce = axialForce[i]

for i in range(numEle):
    axialForceStd[i] = axialForce[i]/absminForce

#遍历单元，绘制变形后结构形状
for i in range(numEle):
    noindex = EleNode[i]
    indexU= 2*noindex[:]-1
    indexV= 2*noindex[:]
    dispU= displacement[0][indexU-1]
    dispV= displacement[0][indexV-1]
    xxDeformed = xx[noindex[:]-1] + dispU*factor
    yyDeformed = yy[noindex[:]-1] + dispV*factor
    xx_Deformed[noindex-1]= xxDeformed
    yy_Deformed[noindex-1]= yyDeformed

    if axialForce[i]>0:
        plt.plot(xxDeformed, yyDeformed, color='green', marker='o', mfc='w', linewidth = 2 * axialForceStd[i][0])
    else:
        plt.plot(xxDeformed, yyDeformed, color='red', marker='o', mfc='w', linewidth = -2 * axialForceStd[i][0])

#遍历单元，标注单元编号
for i in range(numEle):
    noindex = EleNode[i]
    x=(xx[noindex[0]-1]+ xx[noindex[1]-1])/2
    y=(yy[noindex[0]-1]+yy[noindex[1]-1])/2
    plt.text(x, y, i+1, color='b', ha='center', va='center')

#设置坐标格式，使得图形以适合的比例进行显示
max_xx = max(max(xx), max(xx_Deformed))
max_yy = max(max(yy), max(yy_Deformed))
min_xx = min(min(xx), min(xx_Deformed))
min_yy = min(min(yy), min(yy_Deformed))
minx= -0.1*(max_xx-min_xx)+min_xx
maxx = 0.1*(max_xx-min_xx)+max_xx
miny= -0.1*(max_yy-min_yy)+min_yy
maxy = 0.1*(max_yy-min_yy)+max_yy
plt.xlim(minx,maxx)
plt.ylim(miny, maxy)
plt.show()