当前课程知识点：Finite Element Method (FEM) Analysis and Applications > 7、Finite element analysis of bar and beam structures > 7.1 Construction and MATLAB programming of bar element in local coordinate system > Video 7.1
返回《Finite Element Method (FEM) Analysis and Applications》慕课在线视频课程列表
返回《Finite Element Method (FEM) Analysis and Applications》慕课在线视频列表
同学们好
首先回顾一下上一讲的主要内容
上一讲我们学习了
基于试函数的经典方法和有限元方法
有限元方法中的自然离散和逼近离散
有限元方法的基本步骤
以及试函数经典方法及有限元方法的比较
这一讲我们将介绍
杆单元构建及MATLAB编程
局部坐标系中间的平面纯弯梁构建
及MATLAB编程
局部坐标系中的一般梁单元的构建
梁单元的坐标变换
分布力的处理
门型框架结构的实例分析及MATLAB的编程
最后我们将介绍门型框架结构的ANSYS分析实例
下面构造杆单元
以及进行杆单元的MATLAB编程
单元研究主要要点就是节点描述
单元上的场描述
然后应用最小势能原理或者虚功原理
来建立单元的方程
对于杆单元来说
两个节点
节点1、节点2
它的长度为l
它的横截面面积为A
弹性模量为E
那么节点1的节点位移为u1，节点力为P1
节点2的位移为u2，节点力为P2
那么我们节点的基本变量
首先是几何坐标
在我们局部坐标系里面
也就是说我们把零点设在节点1
然后沿着这个单元的轴向设为x轴
那么我们节点1坐标就是x1=0
节点2就是x2=l
节点位移呢，节点1是u1，节点2是u2
我们把它写成一个列向量qe
可以看出来它是具有两个自由度的节点位移
节点力，同样我们把P1，P2写成一个列向量
我们记为Pe
注意一下，为了表达是单元的参数
所以我们在长度、横截面面积、弹性模量上
我们加了一个element的上标
下面进行单元上的场描述
场的描述也涉及到三大类变量
第一个是位移场
那么位移场要基于节点的位移来进行描述
所以说我们要取位移模式
这个位移模式里的待定系数
要和我们节点的位移值的个数要形成一定的关系
那么我们确定位移模式的两个原则
就是唯一确定性原则
还有一个是从低阶到高阶取多项式
那么我们两个节点，所以说我们取a0和a1
两个系数
那么节点条件u(x)在x=0的时候
它要等于u1
在x=l的时候要等于u2
我们把节点条件代到位移模式里去
就可以把两个待定系数，就是a0和a1
分别求出来
那求出来的这个待定系数再代回位移模式的函数中
我们最后获得的位移场的函数的描述就是
那么我们把u1和u2分别整理一下
我们把u1整理到一块就是1减x除上l
u2呢，前面的系数就是x除上l
那么我们把u1和u2写成节点的列向量，就是qe
那前面的这两项同样也写成两项
就是1减x除上l
另外还有一项是x除上l
这两项就组成形状函数矩阵
当我们得到了位移场的描述以后
我们由几何方程来求应变场
那么我们几何方程对于一维问题来说
它就是对位移求一阶导数
针对x求一阶导数以后
我们就得到了这么一个关系
同样我们把它分别写成两个向量相乘
一个是节点的向量，就是u1,u2
另外还有前面的系数
前面这个系数我们称为B(x)矩阵
那么这个B(x)矩阵很显然
它是针对我们前面位移场函数中的N(x)
求一阶导数所得到的
由几何方程得到了应变场函数以后
我们再来由物理方程求应力场函数
那对于一维问题来说
就是将应变再乘上一个弹性模量就可以得到应力
那么我们得到的表达就是在原来的基础上
前面乘一个弹性模量
同样我们也把它表达成基于节点的位移描述
前面乘上一个相应的矩阵
这个矩阵我们做出来就叫作应力矩阵
具体的表达就是两项
一个是负E除上l
另外还有一项是E除上l
在得到了单元的场的表达的基础上
我们来应用最小势能原理求取单元的描述
也就是单元的刚度方程
那么我们要算一下单元的应变能
应变能对于一维问题来说很简单
它就是应力乘上应变，在整个单元的域里积分
当然前面还有一个1/2
那么我们分别把应力和应变基于节点的描述
也就是说基于q的描述，分别代进去
就得到这么一个表达
注意，我们一维问题里面
应力和应变实际上是单独的两项
所以它就是直接相乘
那对于具有很多分项的应力和应变
两个相乘的话，前面一个一定要转置
才能和后面一个相乘，和对应的项数相乘
那么我们把前面的表达代进去以后
我们再整理一下就可以得到
基于节点的位移的描述
前面是节点位移的转置，后面是节点位移
那么中间这个是一个2X2的矩阵
我们把中间这个矩阵叫作Ke
也叫作单元刚度矩阵
它具体算出来就是
这是一个2X2的矩阵
外力功，对于杆单元来说
两个节点力，也就是P1、P2
分别乘上对应的节点位移
也就是u1和u2
这就是外力功
那么我们把它写成一个矩阵形式
同样也是P矩阵转置乘上一个q矩阵
这个q矩阵就是我们的节点位移
前面这个P矩阵就是P1和P2，就是节点力
在得到了应变能和外力功的表达
我们利用最小势能原理
最小势能原理它要计算这个单元的势能
单元的势能就是应变能减外力功
把这个应变能和外力功分别代进去
就得到这个表达
那么这个单元的势能函数取极小值
由于我们的基本变量是节点位移q
所以说取极值就是针对q来求一阶导数
让它等于0
这样我们就得到了
当然这是单元的
那么这个K就是2X2的
q是节点位移，是2X1的
P是节点力，也是2X1的
这就是单元的刚度方程
我们用虚功原理来构建杆单元的刚度方程
那么虚功原理涉及虚位移
对于节点位移来说，我们是qe
那么虚节点位移就是δqe
这个δ表明是对节点位移有一个微小的变化
那么下面我们来计算单元的虚应变能
单元的虚应变能就是δU
它是等于应力乘上一个虚应变
然后对整个单元这个域进行积分
应力我们前面已经得到了这个表达
它也可以表达成节点位移的这么一个关系
虚应变我们在应变前面加上一个δ
让它作一个微小的变化
同样这个δ也作用在虚的节点位移上
也就是说δq
把这两项代进去以后
我们同样可以得到这么一个表达
前面是节点位移，后面是虚的节点位移
中间这个和前面一样叫Ke矩阵
那么这个Ke矩阵算出来也是这样的
也是一个2X2的
单元的虚外力功同样也是这样的
就是由节点力乘上对应节点的虚位移
那这样我们同样可以表达成
节点位移的列向量转置乘上虚的节点位移δq
那么由虚功原理
虚应变能等于虚外力功
那么我们把前面已经得到的这两项分别代进去
我们就可以得到这么一个方程
这个方程两边都有一个虚节点位移δq
那么我们把它移项
最后写成这么一个表达，它要等于0
由于虚位移
当然这个地方的虚位移是节点的虚位移
它也是具有任意性，不恒为零
要使得的这个关系恒为零
就让前面这一项为0
当然我们把它转置一下
就得到了我们常规的单元的刚度方程
那可以看出
虚功原理和最小势能原理得到的结果
也就是说得到的单元的刚度方程是一样的
对于杆单元
我们前面构建了局部坐标系的单元的刚度方程
那么下面我们讨论一下杆单元的坐标变换
首先对于2D问题
我们局部坐标系ox上面的节点位移u1,u2
整体坐标系我们是xoy，上面带一横
那它的分别的节点位移分量是u1,u2还有v1,v2
分别上面带一横
那么我们把整体坐标系的节点位移分量
向局部坐标系上节点位移进行投影
我们就有这么一个投影的关系
我们把投影关系写成矩阵形式
也就是把u1,u2写成一个列向量
这样我们就得到这么一个转换矩阵乘上一个
节点位移列阵的这么一个关系
其中这个Te就是2D问题单元的坐标转换矩阵
它是cosα,sinα组成的2X4的矩阵
对于3D问题我们也可以得到类似的坐标转换矩阵Te
这个Te是一个2X6的
它是两两坐标轴之间的方向余弦
组成的这么一个2X6的矩阵
有了局部坐标和整体坐标的节点位移的坐标转换关系
我们来计算一下单元的势能
单元的势能大家都知道
势能是一个标量，它和坐标系选择是无关的
所以我们先基于局部坐标系写出单元的势能
也就是说单元的应变能和外力功
然后再把单元的节点位移的刚才的坐标变换
把它从局部坐标变换到整体坐标
变换以后我们就得到基于整体坐标q一横的
应变能减外力功的描述
这个时候中间应变能的矩阵K一横
就是整体坐标系下的刚度矩阵
它就等于T矩阵转置
乘上原来局部坐标系下的K矩阵
再乘上一个转换矩阵T
这样我们就得到这么一个4X4的整体坐标系下的刚度矩阵
同样在整体坐标系的节点力
也是在局部坐标系的节点力的基础上
前乘一个坐标变换阵的转置
对于一维的杆单元
我们用MATLAB来实现它的整个计算的流程
主要有4个MATLAB函数
一个是Stiffness，计算刚度矩阵的
第二个是进行刚度矩阵的装配，叫Assembly
第三个是单元应力的计算，叫Stress
最后一个是单元节点力的计算的函数，叫Force
在MATLAB程序里面
我们这四个实现1D杆单元分析的MATLAB函数Stiffness
它是计算单元的刚度矩阵
输入弹性模量E，横截面面积A和长度L
输出2X2的刚度矩阵K
对于单元组装的函数Assembly
它主要是输入单元的刚度阵K
单元的节点编号i,j
输出整体刚度矩阵KK
对于单元应力计算的函数
它输入单元刚度矩阵K
单元的位移列阵u
以及横截面面积A
输出单元的应力Stress
对于Force这个函数
它是输入单元的刚度矩阵K和单元的位移列阵u
最后可以输出单元节点力的分量Force
对于2D杆单元
它的MATLAB程序也是有四个MATLAB函数
这个和前面不同的地方是有一个坐标变换
同样也是Stiffness,Assembly,Stress,Force这四个函数
对于Stiffness函数，它是计算单元的刚度矩阵
输入弹性模量，横截面面积
输入第一个节点的坐标x1,y1
第二个节点的坐标x2,y2
同样还有一个单元的倾角α
当然这个α其实是可以通过这两个节点来进行计算的
那输出是单元的刚度矩阵
是一个4X4的
对于Assembly这个函数
它是进行单元刚度矩阵的组装
需要输入单元的刚度矩阵K
单元的节点编号i,j
最后输出的是整体的刚度矩阵KK
对于Stress函数，它是计算单元的应力
需要输入弹性模量E
第一个节点坐标还有第二个节点坐标
同样也要输入单元的角度α以及单元节点的位移矢量u
它是计算单元的应力Stress
对于Force这个函数
它是计算单元的应力
输入弹性模量E，横截面面积A
输入第一个节点坐标和第二个节点的坐标
以及这个单元的轴线和整体坐标的角度α
还需要输入单元节点的位移矢量u
它计算单元的节点力Force
我们采用MATLAB编程来进行一个四杆桁架结构的分析
我们这个例子是选用知识点2.4中间的四杆结构
首先还是对这个结构进行离散化和编号
离散化的过程和编号和知识点2.4的相同
我们需要输出弹性模量、横截面面积、四个点的坐标
以及各个单元的α角度
然后分别针对单元1,2,3,4调用四次Stiffness
就是单元的刚度阵计算
这样就可以得到单元的刚度矩阵
那么这是MATLAB窗口的输入、输出情况
那么我们可以看见
这是弹性模量、横截面面积
以及节点1、2、3、4的坐标
以及几个单元的α角
那么这是第一次调用Stiffness计算单元1的刚度阵
这是第二次调用，这是计算单元2的刚度阵
这是第三次、第四次
在得到单元的刚度矩阵以后
我们就要建立整体的刚度方程
由于这个结构有四个节点
那么我们总的刚度矩阵
因为每个节点是两个自由度
所以它是8X8的
首先对整体刚度矩阵进行清零
然后四次调用函数Assembly
也就是说每个单元调用一次
每一次调用过后就进行单元的组装
那么我们看看KK=zeros(8,8)这是清零
然后分别调用四次，针对k1,k2,k3,k4
进行整体的组装
那么这是组装以后的整体刚度阵的输出情况
我们下面进行边界条件的处理
我们看下对于节点1和节点4
它x方向和y方向完全固定
对节点2，y方向是约束的
所以说是u1,v1还有v2,以及u4,v4分别等于0
那对于我们总的节点位移q
从u1,v1一直到u4,v4
那么就让对应的位置为0
就是处理它的刚度方程
实际上就是把所对应的行和列划掉
总的节点力包含节点的支反力
以及所施加的外力
我们叫总的节点力
那对于节点1来说它实际上是x方向和y方向的支反力
对于节点2来说x方向它是外加的力
y方向是支反力
对于节点3来说,x方向外加的力为0
y方向施加了一个-25000N的力
对于节点4来说这两个方向的力都是支反力
那么我们调用MATLAB的总体刚度矩阵
我们把处理边界条件以后
剩下的节点2的x方向还有节点3的x方向和y方向
我们把它对应的刚度矩阵里面对应的数
把它这个矩阵，实际上就是3X3的，给取出来
然后把这个矩阵求逆，再乘上对应的节点力
这样我们就获得了节点2的x方向和
节点3的x方向和y方向的节点位移
再把原来零边界条件的几个节点位移填进去
这样我们就得到整体的
所有节点位移的节点的分量的值
支反力的计算
首先我们把总的节点位移列出来
将总的节点位移乘上总的刚度阵
这样就得到所有的节点力的结果
这个所有的节点力的结果也包括
外面施加的节点力还有我们的支反力
当然你要看相应的位置
那么按照对应的关系
我们就得到相应的支反力
也就是说是1号节点的x方向和y方向的支反力
还有4号节点x方向和y方向的支反力
以及2号节点y方向的支反力
最后我们需要计算各个单元的应力
那么各个单元的应力只需要从整体的节点位移列阵q
提出你所需要求的单元的位移
把它的列阵取出来
然后再调用单元的应力函数Stress
就可以计算各个单元的应力
当然也可以调用Force来计算单元的集中的力
这就是分别针对单元1所计算的单元的应力
这是单元2的单元应力
这是单元3、单元4
-Finite element, infinite capabilities
--Video
-1.1 Classification of mechanics：particle、rigid body、deformed body mechanics
--1.1 Test
-1.2 Main points for deformed body mechanics
--1.2 Test
-1.3 Methods to solve differential equation solving method
--1.3 Test
-1.4 Function approximation
--1.4 Test
-1.5 Function approximation defined on complex domains
--1.5 Test
-1.6 The core of finite element: subdomain function approximation for complex domains
--1.6 Test
-1.7 History and software of FEM development
--1.7 Test
-Discussion
-Homework
-2.1 Principles of mechanic analysis of springs
--2.1 Test
-2.2 Comparison between spring element and bar element
--2.2 Test
-2.3 Coordinate transformation of bar element
--2.3 Test
-2.4 An example of a four-bar structure
--2.4 Test
-2.5 ANSYS case analysis of four-bar structure
--ANSYS
-Discussion
-3.1 Mechanical description and basic assumptions for deformed body
--3.1 Test
-3.2 Index notation
--3.2 Test
-3.3 Thoughts on three major variables and three major equations
--3.3 Test
-3.4 Test
-3.4 Construction of equilibrium Equation of Plane Problem
-3.5 Test
-3.5 Construction of strain-displacement relations for plane problems
-3.6 Test
-3.6 Construction of constitutive relations for plane problems
-3.7 Test
-3.7 Two kinds of boundary conditions
- Discussion
-- Discussion
-4.1 Test
-4.1 Discussion of several special cases
-4.2 Test
-4.2 A complete solution of a simple bar under uniaxial tension based on elastic mechanics
-4.3 Test
-4.3 The description and solution of plane beam under pure bending
-4.4 Test
-4.4 Complete description of 3D elastic problem
-4.5 Test
-4.5 Description and understanding of tensor
-Discussion
-5.1 Test
-5.1Main method classification and trial function method for solving deformed body mechanics equation
-5.2 Test
-5.2 Trial function method for solving pure bending beam: residual value method
-5.3 Test
-5.3How to reduce the order of the derivative of trial function
-5.4 Test
-5.4 The principle of virtual work for solving plane bending beam
-5.5 Test
-5.5 The variational basis of the principle of minimum potential energy for solving the plane bending
-5.6 Test
-5.6 The general energy principle of elastic problem
-Discussion
-6.1Test
-6.1 Classic method and finite element method based on trial function
-6.2 Test
-6.2 Natural discretization and approximated discretization in finite element method
-6.3 Test
-6.3 Basic steps in the finite element method
-6.4 Test
-6.4 Comparison of classic method and finite element method
-Discussion
-7.1 Test
-7.1 Construction and MATLAB programming of bar element in local coordinate system
-7.2 Test
-7.2 Construction and MATLAB programming of plane pure bending beam element in local coordinate syste
-7.3 Construction of three-dimensional beam element in local coordinate system
-7.4 Test
-7.4 Beam element coordinate transformation
-7.5 Test
-7.5 Treatment of distributed force
-7.6 Case Analysis and MATLAB programming of portal frame structure
-7.7 ANSYS case analysis of portal frame structure
-8.1 Test
-8.1 Two-dimensional 3-node triangular element and MATLAB programming
-8.2 Test
-8.2 Two-dimensional 4-node rectangular element and MATLAB programming
-8.3 Test
-8.3 Axisymmetric element
-8.4 Test
-8.4 Treatment of distributed force
-8.5 MATLAB programming of 2D plane rectangular thin plate
-8.6 Finite element GUI operation and command flow of a plane rectangular thin plate on ANSYS softwar
-Discussion
-9.1 Three-dimensional 4-node tetrahedral element and MATLAB programming
-9.2 Three-dimensional 8-node hexahedral element and MATLAB programming
-9.3 Principle of the isoparametric element
-9.4Test
-9.4Numerical integration
-9.5 MATLAB programming for typical 2D problems
-9.6 ANSYS analysis case of typical 3Dl problem
-Discussion
-10.1Test
-10.1Node number and storage bandwidth
-10.2Test
-10.2 Properties of shape function matrix and stiffness matrix
-10.3Test
-10.3 Treatment of boundary conditions and calculation of reaction forces
-10.4Test
-10.4 Requirements for construction and convergence of displacement function
-10.5Test
-10.5C0 element and C1 element
-10.6 Test
-10.6 Patch test of element
-10.7 Test
-10.7 Accuracy and property of numerical solutions of finite element analysis
-10.8Test
-10.8 Error and average processing of element stress calculation result
-10.9 Test
-10.9 Error control and the accuracy improving method of h method and p method
-Discussion
-11.1 Test
-11.1 1D high-order element
-11.2 Test
-11.2 2D high-order element
-11.3 Test
-11.3 3D high-order element
-11.4 Test
-11.4 Bending plate element based on thin plate theory
-11.5 Test
-11.5 Sub-structure and super-element
-12.1Test
-12.1 Finite element analysis for structural vibration: basic principle
-12.2 Test
-12.2 Case of finite element analysis for structural vibration
-12.3 Test
-12.3 Finite element analysis for elastic-plastic problems: basic principle
-12.4 Test
-12.4 Finite element analysis for elastic-plastic problems: solving non-linear equations
-Discussion
-13.1 Test
-13.1 Finite element analysis for heat transfer: basic principle
-13.2 Test
-13.2 Case of finite element analysis for heat transfer
-13.3 Test
-13.3 Finite element analysis for thermal stress problems: basic principle
-13.4 Test
-13.4 Finite element analysis for thermal stress problems: solving non-linear equation
-Discussion
-2D problem: finite element analysis of a 2D perforated plate
-3D problem: meshing control of a flower-shaped chuck
-Modal analysis of vibration: Modal analysis of a cable-stayed bridge
-Elastic-plastic analysis: elastic-plastic analysis of a thick-walled cylinder under internal pressur
-Heat transfer analysis: transient problem of temperature field during steel cylinder cooling process
-Thermal stress analysis: temperature and assembly stress analysis of truss structure
-Probability of structure: Probabilistic design analysis of large hydraulic press frame
-Modeling and application of methods: Modeling and analysis of p-type elements for plane problem
