面向对象有限元编程|整体结构类(下)
结构对象的重要任务是组装结构刚度矩阵和整体节点力向量。而结构刚度矩阵和整体节点力向量又依赖一个重要的参数--结构总自由度或者总自由度(gdof)。
计算结构总自由度时，可以先将约束(边界)排除，这种方法叫先处理法。也可以不排除边界条件计算结构自由度，这种方法叫后处理法。
▲先处理法结构标识，有约束的地方自由度为0
如何计算结构总自由度？
以先处理法为例，在创建节点对象时，定义了成员变量RX和RY，分别表示节点X和Y方向的约束，若有约束则其值为0，若没有约束则其值为1. 如图所示
初始时对结构总自由度变量gdof赋值为0，对所有的节点遍历，若RX或者RY的值为1，则gdof自加1，与此同时将gdof赋值给RX或者RY。若RX或者RY的值为0，则gdof的值不变。如图所示
#计算gdof
def getGdof(self):
gdof = 0
for i in range( len(self.listNode) ):
if self.listNode[i].RX == 1:
gdof = gdof + 1
self.listNode[i].RX = gdof
if self.listNode[i].RY == 1:
gdof = gdof + 1
self.listNode[i].RY = gdof
return gdof
如何组装总刚度矩阵？
以下思路同样适用于组装整体节点力向量。定义一个一维数组elToStr，赋值为每个单元的节点自由度编号，这和单元位移向量v=[v1,θ1,v2,θ2]是一一对应的。比如下图中单元1的elToStr=[0,0,0,1],单元2的elToStr=[0,1,0,2],单元3的elToStr=[0,2,0,3]，它是由单元向结构转化的枢纽。
gdof=3，则结构刚度矩阵就是3行3列。开始组装结构刚度矩阵时，结构刚度矩阵赋值全为0。单元1只有一个节点自由度θ2没有被约束，单元1的elToStr只有第四个元素不为0，借鉴划行划列思路，将单元刚度矩阵中被约束的部分划掉，通过数组elToStr将余下的元素添加到结构刚度矩阵的对应位置 1行1列
此时结构刚度矩阵为
单元2有2个节点自由度θ1,θ2没有被约束，按照上述思路
此时结构刚度矩阵为
单元3有2个节点自由度θ1,θ2没有被约束，按照上述思路 此时结构刚度矩阵为
#elToStr数组
def elemToStruct(self, i):
elToStr = [0, 0, 0, 0]
elToStr[0] = self.listElement[i].node1.RX
elToStr[1] = self.listElement[i].node1.RY
elToStr[2] = self.listElement[i].node2.RX
elToStr[3] = self.listElement[i].node2.RY
return elToStr
#组装结构刚度矩阵KK
def packageStiffnessMatrix(self, elToStr, ek, KK):
for m in range ( len(elToStr) ):
for n in range ( len(elToStr) ):
if (elToStr[m] != 0) and (elToStr[n] != 0) :
KK[ elToStr[m]-1 ][ elToStr[n]-1 ] = \
KK[ elToStr[m]-1 ][ elToStr[n]-1 ] + ek [ m ][ n ]
一个完整的桁架结构对象类
class FEModel:
def __init__ (self, listNode, listElement):
self.listNode = listNode
self.listElement = listElement
self.FF = None
def getGdof(self):
gdof = 0
for i in range( len(self.listNode) ):
if self.listNode[i].RX == 1:
gdof = gdof + 1
self.listNode[i].RX = gdof
if self.listNode[i].RY == 1:
gdof = gdof + 1
self.listNode[i].RY = gdof
return gdof
#组装节点力向量
def packageForce(self):
for i in range( len(self.listNode) ):
if self.listNode[i].RX != 0:
self.FF[self.listNode[i].RX-1 ] = self.listNode[i].fx
if self.listNode[i].RY != 0 :
self.FF[self.listNode[i].RY-1 ] = self.listNode[i].fy
def elemToStruct(self, i):
elToStr = [0, 0, 0, 0]
elToStr[0] = self.listElement[i].node1.RX
elToStr[1] = self.listElement[i].node1.RY
elToStr[2] = self.listElement[i].node2.RX
elToStr[3] = self.listElement[i].node2.RY
return elToStr
def packageStiffnessMatrix(self, i,elToStr, ek, KK):
for m in range ( len(elToStr) ):
for n in range ( len(elToStr) ):
if (elToStr[m] != 0) and (elToStr[n] != 0) :
KK[ elToStr[m]-1 ][ elToStr[n]-1 ] = \
KK[ elToStr[m]-1 ][ elToStr[n]-1 ] + ek [ m ][ n ]