python求解结构力学桁架系统

已知材料的弹性模量E=210GPa，上，下悬架截面面积为A1=31.2e-2m²，腹杆的截面面积为A2=8.16e-2m²。求解如下图所示的等跨桁架系统（其中红色箭头为1000KN，黑色箭头为1900KN，该图为等跨的桁架系统）

以下为python代码：

import numpy as np                     # 新增，用于 np.sqrtimport matplotlib                      # 先导入顶层 matplotlibmatplotlib.use('TkAgg')                # 设置后端，必须在导入 pyplot 之前from feon.sa import *from feon.tools import pair_wisefrom feon.sa.draw2d import *import matplotlib.pyplot as plt        # 现在导入 pyplotfrom matplotlib.ticker import MultipleLocatorif __name__ == "__main__":    E = 210e9    A1 = 31.2e-2    A2 = 8.16e-2    nds1 = []    nds2 = []for i in range(13):        nds1.append(Node(i,0))for i in range(11):        nds2.append(Node(i+1,-1))els = []for e in pair_wise(nds1):        el.append(Link2D11((e[0],e[1]),E,A1))for e in pair_wise(nds2):        el.append(Link2D11((e[0],e[1]),E,A1))for i in range(6):        el.append(Link2D11((nds1[i],nds2[i]),E,A2))for i in range(6):        el.append(Link2D11((nds2[i+5],nds1[i+7]),E,A2))for i in range(11):        el.append(Link2D11((nds1[i+1],nds2[i]),E,A2))#创建系统

s = System()

#添加节点

 s.add_nodes(nds1,nds2) s.add_elements(el) s.add_node_force(nds1[0].ID,Fy = -1000)  s.add_node_force(nds1[-1].ID,Fy = -1000)for i in range(1,12):        s.add_node_force(nds1[i].ID,Fy = -1900)        s.add_fixed_sup(nds1[0].ID)        s.add_rolled_sup(nds1[-1].ID,"y")        s.solve()#以下代码为绘制轴力和位移图：

    disp = [np.sqrt(nd.disp["Ux"]**2+nd.disp["Uy"]**2) for nd in s.get_nodes()]    eforce = [el.force["N"][0][0] for el in s.get_elements()]    fig = plt.figure()    ax = fig.add_subplot(211)    ax.yaxis.get_major_formatter().set_powerlimits((0,1))     ax2 = fig.add_subplot(212)    ax2.yaxis.get_major_formatter().set_powerlimits((0,1))     ax.set_xlabel(r"$Node ID$")    ax.set_ylabel(r"$Disp/m$")    ax.set_ylim([-0.00005,0.00005])    ax.set_xlim([-1,27])    ax.xaxis.set_minor_locator(MultipleLocator(1))    ax.plot(range(len(disp)),disp,"r*-")    ax2.set_xlabel(r"$Element ID$")    ax2.set_xlim([-1,46])    ax2.set_ylabel(r"$N/kN$")    ax2.set_ylim(-40000,40000)    ax2.xaxis.set_minor_locator(MultipleLocator(1))for i in range(len(eforce)):        ax2.plot([i-0.5,i+0.5],[eforce[i],eforce[i]],"ks-",ms = 3)  plt.show()

在运行完待代码后得出节点的位移和单元轴力如下图所示：

注：也可在交互式解释器上访问节点或单元信息等。