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

1

为了便于计算，桁架弹性模量与截面积乘积 E A = 1 ，所有杆件均为二力杆 \qquad为了便于计算，桁架弹性模量与截面积乘积EA=1，所有杆件均为二力杆为了便于计算，桁架弹性模量与截面积乘积EA=1，所有杆件均为二力杆
只讨论轴向变形。三结构均在左上节点受一水平向右，和竖直向下集中力，均为 10 K N 只讨论轴向变形。三结构均在左上节点受一水平向右，和竖直向下集中力，均为10KN 只讨论轴向变形。三结构均在左上节点受一水平向右，和竖直向下集中力，均为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 轴正向，旋转角不得大于 π , 倒数第二个必须是 [ 2 , 0 ] , 不能是 [ 0 , 2 ] 倒数第一个必须是 [ 1 , 3 ] , 不能是 [ 3 , 1 ] 否则直接影响节点位移计算结果 ! 单元正向逆时针旋转到整体坐标系X轴正向，旋转角不得大于\pi,\\ 倒数第二个必须是[2,0],不能是[0,2]\\ 倒数第一个必须是[1,3],不能是[3,1]\\ 否则直接影响节点位移计算结果! 单元正向逆时针旋转到整体坐标系X轴正向，旋转角不得大于π,倒数第二个必须是[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 指定放松约束，标号在整体结构中顺次递增整体坐标系下单元刚度矩阵的定位向量（单刚中的行列非0元素在总刚中的位置）\\ 定位向量在单刚组装成总刚的时候很重要\\ 每一行存储了每个单元的俩端节点约束信息\\ 分别为u_1，v_1，u_2，v_2四个自由度约束与否的布尔值\\ 定位向量非0元素为结构中尚未约束的自由度编号\\ 定位向量中的最大值，为节点所有未约束自由度的数目\\ 四个自由度，非0指定放松约束，标号在整体结构中顺次递增整体坐标系下单元刚度矩阵的定位向量（单刚中的行列非0元素在总刚中的位置）定位向量在单刚组装成总刚的时候很重要每一行存储了每个单元的俩端节点约束信息分别为u1，v1，u2，v2四个自由度约束与否的布尔值定位向量非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 ] 对应的杆单元的轴力变成了正解相反数上面的定位矩阵第四行[0,0,3,4]改成了[3,4,0,0],\\ 虽然位移计算结果依然正确\\ 但是[0,0,3,4]对应的杆单元的轴力变成了正解相反数上面的定位矩阵第四行[0,0,3,4]改成了[3,4,0,0],虽然位移计算结果依然正确但是[0,0,3,4]对应的杆单元的轴力变成了正解相反数
注意： 1 ）单元正向逆时针旋转到整体坐标系 X 轴正向，旋转角不得大于 π ，这必须满足 2 ）当单元正向与整体坐标系 X 轴正向的夹角小于 90 ° 时，定位向量编号次序与 e l e m N o d e 单元节点次序一致， 3 ）当单元正向与整体坐标系 X 轴正向的夹角大于 90 ° 时，定位向量编号次序与 e l e m N o d e 单元节点次序相反，所以 l a m b d a _ k 第四行单元节点次序为 [ 2 , 0 ] ，但是其单元正向与整体坐标系 X 轴正向的夹角大于 90 ° ，故其定位向量编号次序为 [ 0 , 0 , 3 , 4 ] ，与 [ 2 , 0 ] 对应的 [ 3 , 4 , 0 , 0 ] 相反。注意：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]相反。注意：1）单元正向逆时针旋转到整体坐标系X轴正向，旋转角不得大于π，这必须满足2）当单元正向与整体坐标系X轴正向的夹角小于90°时，定位向量编号次序与elemNode单元节点次序一致，3）当单元正向与整体坐标系X轴正向的夹角大于90°时，定位向量编号次序与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):#坐标转换矩阵Treturn 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]])-1a_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

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,0l,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

对比分析

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

一次超静定一次超静定一次超静定

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

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

Python与Ansys apdl有限元系列一：平面2D桁架竖向集中载荷
1. 二维静定桁架,结构各个参数如下弹性模量200000MPa, 三节点竖向集中力均为100000N, 杆件横截面积4532mm2 桁架高度为3000mm 桁架横杆(包括底部受拉杆,顶部受拉杆)长度 ...
Python与Ansys apdl有限元系列三：单层单跨梁单元受竖向均布力，水平集中力
结构工况 ################################ # Author: GuaDiKaoLa # Email: 582392629@qq.com ############### ...
结构力学计算机矩阵位移法,结构力学-矩阵位移法.ppt
<结构力学-矩阵位移法.ppt>由会员分享,可在线阅读,更多相关<结构力学-矩阵位移法.ppt(42页珍藏版)>请在装配图网上搜索. 1.结构力学,2 / 42,第七节平面刚 ...
有限元matlab_“ANSYS APDL有限元高级分析技术与二次开发”研修班
课程背景 APDL参数化设计语言,作为ANSYS Mechanical高级分析技术之一,是ANSYS高级用户不可或缺的应用技术之一.为提高广大学员利用ANSYS软件解决实际工程问题的能力,宏新 ...
矩阵位移法matlab编程,矩阵位移法_MATLAB_GUI.doc
Matrix_Displacement_Method--by MATLAB GUI PAGE58 / NUMPAGES64 yanfeng39@zju.edu.cn <结构力学>课程设计之 ...
matlab 矩阵位移法编程结构力学,matlab 矩阵位移法编程结构力学
矩阵位移法编程大作业 (091210211) 一.编制原理本程序的原理是基于结构力学矩阵位移法原理,以结构结点位移作基本未知量,将要分析的结构拆成已知节点力-结点力位移关系的单跨梁集合,通过强令结构 ...
矩阵位移法是用于求解杆系结构的计算机方法,结构力学的教学思路
结构力学的教学思路 2019-10-14 版权声明举报文章一.经典位移法(以下简称位移法)和矩阵位移法都是求解杆系结构的基本方法,是结构力学课程中两个十分重要的内容,两种方法都是结构力学课程中讲授 ...
matlab 矩阵位移法编程结构力学,matlab 矩阵位移法编程结构力学.doc
matlab 矩阵位移法编程结构力学.doc 矩阵位移法编程大作业(091210211)一.编制原理本程序的原理是基于结构力学矩阵位移法原理,以结构结点位移作基本未知量,将要分析的结构拆成已知节点力 ...
Python小案例（六）通过熵权法计算指标权重
Python小案例(六)通过熵权法计算指标权重在日常业务中,产品运营需要综合多个指标进行判断,如果没有目标变量进行监督训练的话,很难人为地判断哪个指标更好,综合起来哪个类别更优秀. 这里介绍一种基于 ...

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

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

1

2

2.1 apdl命令流

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

最新文章

热门文章