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

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