面向对象有限元编程|整体结构类(下)

结构对象的重要任务是组装结构刚度矩阵和整体节点力向量。而结构刚度矩阵和整体节点力向量又依赖一个重要的参数--结构总自由度或者总自由度(gdof)。

计算结构总自由度时，可以先将约束(边界)排除，这种方法叫先处理法。也可以不排除边界条件计算结构自由度，这种方法叫后处理法。

▲先处理法结构标识，有约束的地方自由度为0

如何计算结构总自由度？

以先处理法为例，在创建节点对象时，定义了成员变量RX和RY，分别表示节点X和Y方向的约束，若有约束则其值为0，若没有约束则其值为1.如图所示 

初始时对结构总自由度变量gdof赋值为0，对所有的节点遍历，若RX或者RY的值为1，则gdof自加1，与此同时将gdof赋值给RX或者RY。若RX或者RY的值为0，则gdof的值不变。如图所示 

#计算gdof
defgetGdof(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数组
defelemToStruct(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
defpackageStiffnessMatrix(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 ] 

一个完整的桁架结构对象类

classFEModel:
def__init__(self, listNode, listElement):
        self.listNode = listNode
        self.listElement = listElement
        self.FF = None
defgetGdof(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
#组装节点力向量
defpackageForce(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  
defelemToStruct(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
defpackageStiffnessMatrix(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 ]