matlab 三弯矩方程_程序思维 | 结构矩阵分析与MATLAB

 

Building可视库

一起做些更酷的事吧！

关注

你好，这里是Building可视库，建筑数据可以看得见的小库。城市建筑作为主体，数据信息作为底层，搭载着可视化的小摩托，做着一些酷酷的事情。

回归到小库本专业，这篇推文主要从程序设计角度出发，以有限单元法的原理为例，结合结构矩阵分析与MATLAB的知识，做一期整理与分享。小库本科阶段，已亲测MATLAB程序通过多个样题验证正确。

问题描述

基于结构力学概念、矩阵理论和MATLAB程序，学会快速求解杆系结构受力情况。

结构矩阵分析方法，又称杆系结构的有限元法，是以结构力学中的位移法为基础，并与矩阵理论相结合所形成的一种以计算机作为数值求解工具的方法。

问题背景

传统的手算方法有力法和位移法，以及两者结合的混合法。对于简单杆系结构来说，尚可利用手算进行解决，但是当单元杆件数、约束条件比较复杂的时候，手算工作量将异常巨大，无法有效的得出结果。

因此考虑利用电算进行解决，衍生出矩阵力法和矩阵位移法，由于后者更容易实现计算过程的程式化而广而应用。矩阵位移法通常求解结构力学超静定静力问题。

几个名词

杆单元：只能承受拉压轴力的杆件

梁单元：弯矩、轴力、剪力都可以承受

单元刚度：引起杆件单位位移所需的力

超静定结构：几何不变但存在多余约束的体系

细分流程

1.

结构离散

将结构的每一根具有单一刚度的直杆在两个节点附近“锯”下来，最终分解成为单元杆件、铰结点、刚结点和基础这四类“零件”。当结构中存在曲杆，可用“微分逼近”的方法将其近似为刚架结构。

2.

单元分析

第一步：确定编号

单元杆件：①、②、③、④

节点编号：1、2、3、4

节点位移：(1,2,3)、(0,0,3)

注：位移编号方式有先处理法和后处理法。

先处理法：将节点的已知位移统一编为0号，未知位移编以互不相同的正整数。刚度形成之前就引入约束条件的方法，对应于后处理法。

后处理法：所有的位移均编以互不相同的正整数。

第二步：建立平面直角坐标系XOY

第三步：对各杆件进行刚度分析

注：由于各杆件系从原结构上“锯”下来的，所以各杆件两端的约束均为固端约束，存在固端轴力、剪力和弯矩。另外，凡是与单元有关的物理量，其表达符号上加一短横。

计算各单元的截面积A，惯性矩I、EA/l、EI/l、EI/(l*l)、EI/(l*l*l)。

计算单元刚度矩阵：

简写为：

矩阵Ke的奇异性——如果A为奇异矩阵，则AX=0有无穷解，AX=b有无穷解或者无解。即如果给定杆端力，不能求出杆端位移，因此需要给定杆端位移，求出杆端力。这也给出了我们之后求解未知量的方向。

推导来源：

在一般单元两端加上人为控制的附加约束，使基本杆单元的两端产生单位的轴向位移u、横向位移v，转角位移φ，然后依据位移法或者其他方法(材料力学知识)来推导出两端的轴力N，剪力Q，弯矩M，进而按照规定编号写成矩阵形式的表达式,即局部坐标系中的单元刚度方程。

3.

坐标转换

在单元分析中，杆端位移和杆端力是在单元坐标系中定义的，所以开展整体分析之前，必须将不同单元的坐标系的杆端力和杆端位移等量值换算到统一的坐标系中，即将局部坐标系中的单元刚度矩阵和荷载列阵通过坐标转换矩阵T变换为总体坐标系中的单刚和荷载列阵。

α是整体坐标系和局部坐标系的夹角。

注：T为上面的六阶方阵，为正交矩阵

荷载列阵转换简写为：

位移列阵转换简写为：

整体坐标系下的单元刚度矩阵

经过矩阵理论等效变化：

将结构中的每根单元单刚和坐标转换矩阵T依次代入上式，计算得到整体坐标系下的单元刚度矩阵。

注：用先处理法得到整体坐标系下的的单元刚度矩阵，经约束处理后形成的单元刚度方程，若存在4个未知位移，则 Ke为4阶方阵，即未知位移数目与Ke矩阵阶数同步。

约束处理详情见后面表达。

4.

整体分析

原理：通过虚位移原理或最小势能原理，可以推导出集成总体刚度矩阵后的刚度方程，如下所示。现在主要求解以下各量，思路为先求得K和R，代入上述方程，求得δ。

K为单元集成后的总体刚度方程

δ为结构的节点位移列向量

R为结构的综合节点荷载列向量

注意：这儿的K经过约束处理后，不再是奇异矩阵，与前面奇异矩阵的性质要区分开。

以3个单元，4个节点的刚架为例

为了便于编程，将整体坐标系下的每个单元刚度矩阵用4个子块表示。注意子块表示的下角标对应单元的节点编号i、j。

为了便于编程，将整体坐标系下的每个单元刚度矩阵用4个子块表示。实际题目中节点编号为1、2、3、4。即：

注：其中每个子块是一个与单位1的轴向位移、横向位移、转角位移对结构产生的轴力、剪力、弯矩对应的3阶方阵。右上角的角标表示单元编号、右下角的两个角标表示节点编号。

在平面杆系中，每个节点有两个线位移和一个转角位移，该刚架有3个节点，共9个节点位移分量，整理成为结构的节点位移列向量。δ1表示节点1的位移列向量。

同理结构的综合节点荷载列阵为：

集成总体刚度矩阵

实质上是把结构中所有的单元刚度矩阵根据结构的节点位置在相应的总体刚度矩阵中投放。这一点对任何有限元方法均成立。

对号入座法——把每个单元单刚对应的四个子块按照下标的号码投放到总体刚度矩阵的相应位置上。

5.

代入荷载

根据“力=单元刚度×位移”原则，即R=Kδ，等号左边应为节点荷载与节间荷载产生的杆端力的之和表示的列向量，目前已经解决了总体刚度矩阵K的问题，先需要导出R，求出δ,目的是利用δ，求出局部坐标系下的单元刚度方程中的杆端力。

结构的综合节点荷载列阵R

除了节点荷载，还可能包括非节点荷载，即杆间集中力P或均布力q。需要将非节点荷载变换为等效节点荷载。

将作用有非节点荷载的单元，两端变成固定端，求出固定端的反力约束，将原节点节点荷载保持不变，再加上于固定端约束反力反方向的端点力作为附加的节点荷载，称为等效节点荷载，将等效节点荷载加入到节点荷载列阵，建立总体刚度方程。

注意：需要进行坐标变换。

可查阅等截面单元固端力计算公式

6.

解方程组

约束处理的两种方法

前处理法：

在程序设计中，在对结构总体广义位移分量进行编号时，对已经广义位移分量不进行编号，这样自然会没有与已知广义位移分量相关的行和列。

充0置1法：

在总体刚度矩阵中，与节点位移分量为所对应的主对角元素置1，将该主元素所在的行和列的副元素全部充0，同时综合节点荷载列阵中与零节点位移分量所对应的元素也充0，等于是划去了位移为0的方程。这种方法变成简单、计算可靠、广泛应用。

求解总体刚度方程

解上式可以得到整体坐标系下的节点位移分量δ ，经过坐标变换后得到局部坐标系中的单元杆端位移列阵，然后根据局部坐标系下的的单元刚度方程计算出杆端力。

一个人至少拥有一个梦想，有一个理由去坚强。心若没有栖息的地方，到哪里都是在流浪。

——三毛

 

平面桁架的程序设计

开始

↓

输入结构的原始数据

↓

计算单元刚度矩阵

↓

计算总体刚度矩阵

↓

生成结构刚度方程组的右端项

↓

引入约束条件

↓

求解结构刚度方程组得到节点位移

↓

按单元循环求解单元内力和应力

↓

结束

01

平面桁架的主控变量

主控变量：用以控制结构中各种几何和材料数据的正确输入，合理解决程序在运行中主循环的开始和终止等问题。

根据平面桁架的几何和受力特点。程序的主控变量定义为5个，分别为：

N——结点总数(包括固结点)

LR——已知结点位移分量总数

MB——桁架单元总数

LTB——单元的材料特征类型值，即各杆的横截面积和弹性模量相同，值为1

PQ——外荷载总数

以此题为例：

所示平面桁架，水平杆的横截面积A=110mm2，弹性模量E=180GPa，竖直杆横截面积A=100mm2，弹性模量E=190GPa，斜杆的横截面积A=90mm2，弹性模量E=200GPa，外力P=1000N，图中单位为mm，计算该平面桁架的结点位移和单元内力与应力。

N=5，LR=4，MB=7，PQ=4，LTB=3

02

平面桁架的数据描述

(1) 结点坐标数组

在计算单元刚度矩阵、单元应力、结点荷载时需用到结点坐标。它是一个N行3列的数组，数组名为XY，其中N为结点总数。

XY=[1,200e-3,0;2,100e-3,100e-3;

3,100e-3,0;4,0,100e-3;5,0,0]

(2) 单元特性数组

用MB行、4列的二维数组ELB表示，第一列代表单元编号、第二三列表示单元节点的编号(小号节点为第一个节点)，第四列代表单元的材料类型号。

ELB= [1,1,3,1;2,1,2,3;3,2,3,2;

4,3,5,1;5,3,4,3;6,2,4,1;7,4,5,2]

(3) 材料特征数组

在计算单元刚度矩阵和单元杆端力时，需要提供单元材料的特征。

EA=[180e9,110e-6;190e9,100e-6;

200e9,90e-6]

(4) 结点荷载数组

用PQ行、2列的数组表示，PQ=4。第一列代表荷载总体自由度方向(与结点编号i，j有关，数值为2i或2i-1，2i对应于Y轴方向力作用下，2i-1对应于X轴方向力作用下)，第二列代表荷载数值，方向与坐标轴同向为正，反向为负。

NPQ=[1,1000;2,-1000;3,1000;

6,-1000]

(5) 位移荷载数组

用LR行、2列的数组表示，第一列表示结点位移编号所对应的方程编号数，例如已知结点i的横向位移分量Ui，对应的方程编号数为2i-1，竖向位移分量Vi，对应的方程编号数为2i。第二列表示位移值。

SU=[7,0;8,0;9,0;10,0]

03

平面桁架的功能函数

说明：

平面桁架不承受弯矩的作用，在局部坐标系下的单元刚度矩阵为和坐标变换矩阵为：

(1) 单元刚度矩阵函数

(2) 总体刚度矩阵函数

(3) 荷载列阵函数

(4) 约束处理函数

(5) 结点位移函数

(6) 单元内力函数

04

平面桁架的主程序

(1) 前处理程序

(2) 主程序

(3) 后处理程序

05

程序分析的结果展示

(1) 前处理工作区

(2) 结点位移

(3) 单元内力及应力

 

平面刚架的程序设计

01

平面刚架的主控变量

主控变量：用以控制结构中各种几何和材料数据的正确输入，合理解决程序在运行中主循环的开始和终止等问题。

根据平面刚架的几何和受力特点。程序的主控变量定义为5个，分别为：

N——结点总数(包括固结点)

LR——已知结点位移分量总数

MB——刚架单元总数

LTB——单元刚度特征类型数

LPQ——外荷载总数

以此题为例：

如图所示的平面刚架，竖直杆的横截面为直径等于26cm的圆，水平杆为边长为23cm的正方形，弹性模量为150GPa，计算平面刚架的杆端位移和内力。

N=5，LR=2，MB=4，PQ=3，LTB=2

02

平面刚架的数据描述

(1) 结点坐标数组XY

在计算单元刚度矩阵、单元应力、结点荷载时需用到结点坐标。它是一个N行3列的数组，数组名为XY，其中N为结点总数。第一列代表节点的编号，第二列代表节点在整体坐标系下的横坐标值，第三列代表节点在整体坐标系下的纵坐标值。

XY=[1,0,0;2,2,0;3,0,3;

4,2,3;5,2,2];

(2) 单元信息数组ELB

用MB行、4列的二维数组ELB表示，第一列代表单元编号，第二列表示单元左端节点编号，第三列表示单元右端节点编号，第四列代表单元刚度类型号。

ELB=[1,1,3,1;2,3,4,2;

3,2,5,1;4,4,5,1];

(3) 刚度特征数组EAIZ

在计算单元刚度矩阵和单元杆端力时，需要提供单元刚度的特征。刚度特征数组是一个LTB行3列的二维数组，第一列代表弹性模量，第二列代表单元截面面积，第三列代表此截面下的惯性矩。

d1=26e-2;

A1=pi*d1*d1/4;

IZ1=pi*d1^4/64;

a2=23e-2;

A2=a2*a2;IZ2=a2^4/12;

EAIZ=[150e9,A1,IZ1;

150e9,A2,IZ2];

(4) 结点荷载数组ELPQ

用LPQ行、4列的二维数组表示，第一列代表荷载所在的单元编号，第二列代表非节点荷载位置与单元左端的距离C(具体见下表)，第三列代表荷载的数值，第四列代表荷载的类型(均布、集中、轴向力、线性分布)。

ELPQ=[2,2,-40e3,1;

2,1,-25e3,2;3,2,12e3,4];

(5) 位移列阵数组SU

用LR行、2列的数组表示，第一列表示节点位移编号所对应的方程编号数，例如已知结点i的横向位移分量Ui，对应的方程编号数为3i-2，竖向位移分量Vi，对应的方程编号数为3i-1，转角位移φi，对应的方程编号数为3i；第二列表示位移值。

SU=[1,0;2,0;3,0;4,0;5,0;6,0];

03

平面刚架的功能函数

(1) 单元刚度矩阵函数

dxy——2行2列的单元节点数组，第一行存放单元左侧端点的横纵坐标值，第二行存放单元右侧端点的横纵坐标值。

E——杆单元的弹性模量

A——杆单元的横截面面积

IZ——杆单元的惯性矩

KE——6行6列的单元刚度矩阵

T——6行6列的坐标转换矩阵

(2) 总体刚度矩阵函数

(3) 荷载简化函数

与平面桁架不同，平面刚架的外荷载并非都在节点上，因此需要将作用在平面刚架上的外荷载向节点等效简化，具体见上述表格即可。

(4) 荷载列阵函数

(5) 结点位移函数

(6) 单元内力函数

04

平面刚架的主程序

(1) 前处理程序

(2) 主程序

(3) 后处理程序

05

程序分析的结果展示

(1) 结点位移

(2) 单元内力

文字 | 小库Q

图片 | 小库Q    Jason

编辑 | Jason