面向对象有限元编程|优化框架_python

面向对象有限元编程|优化框架

上一次提出的有限元程序框架面临一个问题：把单元都放在一个模块文件modelements，导致代码过多，不方便查找。并且也不能体现出这个框架的物理含义，即单元的多样性。同样，后处理模块是一种单元对应一种后处理方法，所以也会定义很多方法。现在把单元拆分，一种单元放在一个模块中。再把所有的单元定义为一个子包Elements，里面除了定义的单元外，还有一个文件__init__.py，这是包的标志。

Python中的包和模块有两种导入方式：精确导入和模糊导入：

精确导入方式：

from Elements import MdTruss

就是从子包Elements导入模块MdTruss,MdTruss.pre 就是调用里面的方法了。

模糊导入方式：

from Elements import *

import numpy as np
from Materials import MdMaterial
from Elements import MdBeam2D
from FiniteModel import MdFEModel
from Numerics import MdLinerSolver
from Posts import MdPost

mat1 = MdMaterial.Material( 1 )
#节点对象
nd1 = MdBeam2D.NodeBeam2D( 1,  0,         0,        1,   0,   0)
nd2 = MdBeam2D.NodeBeam2D( 2,  4,         0,        1,   0,   0)
nd3 = MdBeam2D.NodeBeam2D( 3,  9,         0,        1,   0,   0)
#单元对象
elem1 = MdBeam2D.ElemBeam2D( 1, 1, mat1, nd1, nd2, 1,   -4.5,    -4.5 )
elem2 = MdBeam2D.ElemBeam2D( 2, 1, mat1, nd2, nd3, 1,   -3.2,   -3.2 )

listNode = [nd1, nd2, nd3]
listElem = [elem1, elem2]

fm = MdFEModel.FEModel(listNode, listElem) #有限元模型实例
ff = fm.getForceVector()
KK = fm.getStructStiffnessMatrix()
cls =  MdLinerSolver.LinerSolver(KK, ff, 1e-4, 500) #创建一个求解器的实例cls
delta = cls.CGsolver() #调用共轭梯度法求解

post = MdPost.Post(fm, delta) #后处理实例
post.outputResult()

post.Visualize()

来源：数值分析与有限元编程

Python画弯矩图剪力图(一)

Python画弯矩图，都以分段直线近似代替曲线，分段越多，越接近真实的弯矩曲线，就像这样import matplotlib.pyplot as pltimport numpy as npfig, axs = plt.subplots(1, 1, constrained_layout=True)axs.plot([0,1,5],[0,2,1], linewidth=2, color="dimgrey")plt.show() 除了大致轮廓，还可以把细节完善一下。比如两端封闭起来，填充颜色等等。以受均布荷载的简支梁为例import matplotlib.pyplot as pltimport numpy as npplt.rcParams['font.sans-serif'] = ['SimHei'] # 正常显示中文标签plt.rcParams['axes.unicode_minus'] = False # 正常显示负号l = 5q = 20x = np.linspace(0,l,20) #将梁20等分M = q/2*(l*x-x**2) #计算等分点弯矩MV = q*(l/2-x) #计算等分点剪力Vfig, axs = plt.subplots(2, 1, constrained_layout=True,figsize=(8,4) )axs[0].invert_yaxis() axs[0].plot(x,M, label="M1", linewidth=3) #画弯矩图axs[0].fill_between(x, 0, M, alpha=0.3, where=(M<=0) ) axs[0].fill_between(x, 0, M, alpha=0.3, where=(M>=0) )axs[0].plot([0,5],[0,0], linewidth = 2, color = "dimgrey") #画梁axs[0].set_ylabel('弯矩(kN.m)', fontsize = 16)fig.suptitle('内力图', fontsize = 18)x1 = np.array( [0] )x2 = np.array( [l] )x3 = np.hstack( (x1, x, x2) )v1 = np.array( [ 0 ] )v3 = np.hstack( (v1, V, v1) )axs[1].plot(x3, v3, label="M1", linewidth=3)#画剪力图axs[1].set_ylabel('剪力(kN)', fontsize = 16)axs[1].fill_between(x, 0, V, alpha=0.3, where=(V<=0) ) #填充颜色axs[1].fill_between(x, 0, V, alpha=0.3, where=(V>=0) )axs[1].plot([0,5],[0,0], linewidth = 2, color = "dimgrey") #画梁fig.savefig('./f2.png', dpi = 500) #保存图片 plt.show() 再看一个多跨梁弯矩图注意事项：正弯矩画在杆件的下方，负弯矩画在杆件的上方。使杆件下部受拉的弯矩为正,上部受拉的弯矩为负。弯矩图画在杆件纤维受拉的一侧。剪力图注意事项：正剪力画在杆件的上方；负剪力画在杆件的下方；使杆件截面顺时针方向转动的剪力为正剪力；使杆件截面逆时针方向转动的剪力为负剪力；一般情况下，剪力与杆件所受外力的方向相反。弯矩图是一条表示杆件不同截面弯矩的曲线。这里所说的曲线是广义的，它包括直线、折线和一般意义的曲线。弯矩图是对构件弯矩的图形表示，弯矩图画在受拉侧，无须标正负号。来源：数值分析与有限元编程