Python与Ansys apdl有限元系列二：矩阵位移法计算桁架结构

矩阵位移法计算三种桁架（静定，一次，二次超静定）

1，Python部分

$\qquad为了便于计算，桁架弹性模量与截面积乘积EA=1，所有杆件均为二力杆$
$只讨论轴向变形。三结构均在左上节点受一水平向右，和竖直向下集中力，均为10KN$

################################
# Author: GuaDiKaoLa
# Email: 582392629@qq.com
################################
import numpy as np
from math import sqrt,sin,cos,acos
np.set_printoptions(precision=2)
np.set_printoptions(suppress=True)
EA = 1

$单元正向逆时针旋转到整体坐标系X轴正向，旋转角不得大于\pi,\\ 倒数第二个必须是[2,0],不能是[0,2]\\ 倒数第一个必须是[1,3],不能是[3,1]\\ 否则直接影响节点位移计算结果!$

nodeCoord = np.array([[0,-2],[0,0],[2,0],[3.5,-2]])#节点坐标
# 单元节点编号的顺序很重要,决定了单元方向
elemNode = np.array([[1,0],[1,2],[2,3],[2,0],[1,3]])#组成单元的俩节点编号

$整体坐标系下单元刚度矩阵的定位向量 （单刚中的行列非0元素 在总刚中的 位置）\\ 定位向量 在 单刚组装成总刚的时候很重要\\ 每一行存储了每个单元的俩端节点约束信息\\ 分别为u_1，v_1，u_2，v_2四个自由度约束与否的布尔值\\ 定位向量非0元素为 结构中尚未约束的自由度编号\\ 定位向量中的最大值，为节点所有未约束自由度的数目\\ 四个自由度，非0指定放松约束，标号在整体结构中顺次递增$

lambda_k =np.array([[0,0,1,2],
                    [1,2,3,4],
                    [3,4,5,0],
                    [0,0,3,4],
                    [1,2,5,0]])
# 施加载荷
Forces = np.array([10,10,0,0,0])

$上面的定位矩阵第四行[0,0,3,4]改成了[3,4,0,0],\\ 虽然位移计算结果依然正确\\ 但是[0,0,3,4]对应的杆单元 的轴力变成了正解相反数$
$注意：1）单元正向逆时针旋转到整体坐标系X轴正向，旋转角不得大于\pi，这必须满足\\ \qquad\quad2）当单元正向与整体坐标系X轴正向的夹角小于90°时，\\ \qquad\quad定位向量编号次序与elemNode单元节点次序一致，\\ \qquad\quad3）当单元正向与整体坐标系X轴正向的夹角大于90°时，\\ \qquad\quad定位向量编号次序与elemNode单元节点次序相反，\\ 所以lambda\_k第四行单元节点次序为[2,0]，\\ 但是其单元正向与整体坐标系X轴正向的夹角大于90°，\\ 故其定位向量编号次序为[0,0,3,4]，与[2,0]对应的[3,4,0,0]相反。$

nodeCoord = np.array([[0,-2],[0,0],[2,0],[3.5,-2]])#节点坐标
#
# 单元节点编号的顺序很重要,决定了单元方向
elemNode = np.array([[1,0],[1,2],[2,3],[2,0],[1,3]])#组成单元的俩节点编号
lambda_k =np.array([[0,0,1,2],
                    [1,2,3,4],
                    [3,4,0,0],
                    [0,0,3,4],
                    [1,2,0,0]])
 # 施加载荷
Forces = np.array([10,10,0,0])

nodeCoord = np.array([[0,-2],[0,0],[2,0],[3.5,-2],[1,-2]])#节点坐标
#
# 单元节点编号的顺序很重要,决定了单元方向
elemNode = np.array([[1,0],[1,2],[2,3],[2,0],[1,3],[2,4]])#组成单元的俩节点编号
lambda_k =np.array([[0,0,1,2],
                    [1,2,3,4],
                    [3,4,0,0],
                    [0,0,3,4],
                    [1,2,0,0],
                    [0,0,3,4]])
#施加载荷
Forces = np.array([10,10,0,0])

nodeX = nodeCoord[:,0]#节点x坐标
nodeY = nodeCoord[:,1]#节点y坐标

# 局部坐标系下单元刚度矩阵k_e 与 桁架单元长度L 的乘积（因L为变值，故先计算恒定的乘积）
k_eL = EA*np.array([
    [1,0,-1,0],
    [0,0,0,0],
    [-1,0,1,0],
    [0,0,0,0]
])



def TransMatrix(angle):#坐标转换矩阵T
    return np.array([
        [ cos(angle), sin(angle) ,      0   ,      0     ],
        [-sin(angle),  cos(angle),      0   ,      0     ],
        [   0       ,      0     , cos(angle), sin(angle)],
        [   0       ,      0     ,-sin(angle), cos(angle)]
])

def K_e(i):# 整体坐标系下单元刚度矩阵,i为单元编号
    nodeIndex = elemNode[i,:]# 第i个单元的节点编号列表
    node_X1_loc = nodeX[nodeIndex[0]]
    node_Y1_loc = nodeY[nodeIndex[0]]
    node_X2_loc = nodeX[nodeIndex[1]]
    node_Y2_loc = nodeY[nodeIndex[1]]
    L = sqrt((node_X2_loc - node_X1_loc) ** 2 + (node_Y2_loc - node_Y1_loc) ** 2)
    angle = acos((node_X2_loc-node_X1_loc)/L)
    k_e = k_eL / L# 局部坐标系下单元刚度矩阵k_e
    # 整体坐标系下单元刚度矩阵 K_e = [T]' * k_e * [T]
    return np.matmul(np.matmul((TransMatrix(angle).T) , k_e),TransMatrix(angle))

for i in range(len(elemNode)):
    print('')
    print('K_e(%d) = '%i ,'\n', K_e(i))


# 结构整体刚度矩阵 KK，单刚K_e ==组装==>> 总刚KK

# 总自由度 = 节点的个数 * 2
numberOfDof = 2 * len(nodeCoord)
# 总刚的大小为 n*n
#           其中n为结构整体节点的所有未约束的自由度数目
#           也即: 定位向量中元素的最大值
numberOfNodeFreeDof = np.max(lambda_k)
KK = np.zeros((numberOfNodeFreeDof,numberOfNodeFreeDof))
##############################################单刚K_e ==组装==>> 总刚KK
# 将单刚矩阵对应的定位向量 置于单刚的行首列首
# 将单刚矩阵中 定位向量元素为0的所在行列全部划掉
# 剩下的元素 依次添加到总刚中，在总刚中的位置就是定位向量
#

for i in range(len(elemNode)):
    for j in lambda_k[i,:]:
        if j != 0:
            loc_kej = np.argwhere(lambda_k[i, :] == j)
            for k in lambda_k[i,:]:
                if k != 0:
                    loc_kek = np.argwhere(lambda_k[i, :] == k)
                    KK[j-1,k-1] += K_e(i)[loc_kej,loc_kek]
print('='*20)
print('KK = ','\n',KK)
#############################单刚K_e ==组装==>> 总刚KK

# 求解位移
displacement = np.linalg.inv(KK).dot(Forces)
print('\n','displacement:','\n',displacement)

# 求解矩阵a_e，各个单元的结点位移向量合并而成
# 目标为 将单元定位向量矩阵的非0自由度编号替换为对应的节点位移
a_e = np.zeros((lambda_k.shape),dtype=float)
a_e += lambda_k
print(a_e)

# 上面的a_e为各个单元的结点位移矩阵合并而成，
# 节点位移矩阵使用 单元定位向量矩阵为原本，将非0的自由度编号替换为对应的节点位移即可
# 1，单元定位向量非0元素的索引为NoZeroIndex
# 2，a_e每一行根据非0元素索引 拿到非0元素
# 3，这些非0元素减一 又是 活动节点位移向量displacement 的索引
# 目标为 将单元定位向量矩阵的非0自由度编号替换为对应的节点位移
# 结点位移==>> displacement索引
#        ==>> a_e非0元素值减一
#        ==>> a_e非0元素值索引
#        ==>> np.argwhere(a_e[i, :] != 0)
#


for i in range(len(elemNode)):
    # 整体坐标系下单元杆端位移
    NoZeroIndex_raw = np.argwhere(lambda_k[i, :] != 0)
    NoZeroNum = NoZeroIndex_raw.shape[0] * NoZeroIndex_raw.shape[1]
    # 因为NoZeroIndex_raw是一个多行一列的数组，但是必须将其展平为一行多列
    # 才能使用其元素作为索引寻址，故拿到它的元素个数后使用reshape展平
    NoZeroIndex = NoZeroIndex_raw.reshape(1, NoZeroNum)
    for j in range(NoZeroNum):
        indexDof = int(a_e[i,:][NoZeroIndex[0,j]])-1
        a_e[i, :][NoZeroIndex[0, j]] = displacement[indexDof]
    print(a_e[i,:].shape)#(4,)

print(a_e)

for i in range(len(elemNode)):
    AxialForceXY[i, :] = np.matmul(K_e(i), a_e[i, :])
print('AxialForceXY','\n',AxialForceXY)

2，Ansys APDL 部分

2.1 apdl命令流

finish
/clear
/prep7
k,1,0,-2,0
k,2,0,0,0
k,3,2,0,0
k,4,3.5,-2,0
k,5,1,-2,0

l,1,2
l,2,3
l,3,4
l,1,3
l,2,4
l,3,5!二次超静定这行保留，

ET,1,LINK1
R,1,4532
MP,EX,1,1/4532
MP,PRXY,1,0.3
LATT,1,1,1
LESIZE,ALL,,,1
LMESH,ALL
DK,1,UX,,,,UY
DK,4,UX,,,,UY
!DK,4,UY
DK,5,UX,,,,UY!约束可自行设置
FORCES = 10
FK,2,FX,FORCES
FK,2,FY,-FORCES
/SOLU
ANTYPE,0
SOLVE
FINISH

$静定结构（很惊讶滑动铰支座水平位移如此之大！）$

$一次超静定$

$二次超静定$
$哈哈，数据基本上是吻合的……$

转载：https://blog.csdn.net/weixin_43717845/article/details/105515372