计算结构力学基础课件

第二章变形体虚位移原理弹性力学基本概念—预备知识变形体虚位移原理和势能原理虚位移原理和势能原理的应用**预备知识（回顾）

线弹性平面问题的平衡方程

小变形平面问题的几何方程

线应变：角应变：**

线弹性平面问题物理方程平面应力：平面应变：平面应力：平面应变：预备知识（回顾）**

平面问题应力边界条件

在应力边界上：

平面问题物理量的矩阵表示应力矩阵应变矩阵体积力矩阵表面力矩阵位移矩阵已知位移矩阵弹性矩阵预备知识（回顾）**

平面问题物理量的矩阵表示

取决于材料性质各相同性、线性弹性时

引入两个算子矩阵微分算子矩阵方向余弦矩阵平面应力平面应变：预备知识（回顾）**基本方程矩阵表示

平衡方程

几何方程

物理方程

边界条件杆系问题的基本方程（作业）

平衡方程如何建立？

几何方程如何建立？

内力和变形间关系如何？

由微段的平衡条件建立

由微段的变形条件建立以上内容必须通过自己动手达到熟练掌握**变形体虚位移原理和势能原理一、变形体虚位移时外力功计算二、变形体虚位移原理表述和证明三、一些名词含义的解释四、势能驻值原理和最小势能原理**变形体虚位移原理和势能原理变形体虚位移时外力功计算内部微元体的受力分析10其余类推**变形体虚位移原理和势能原理变形体虚位移时外力功计算内部微元体的位移分析虚位移算子符号能写出各点的位移吗？提示：连续函数台劳级数展开**变形体虚位移原理和势能原理变形体虚位移时外力功计算内部微元体的外力功计算8y方向的力所做的功等于多少？请大家自行写出内部微元体x向外力的总虚功**变形体虚位移原理和势能原理变形体虚位移时外力功计算内部微元体的外力功计算原形刚性位移变形位移**变形体虚位移原理和势能原理变形体虚位移时外力功计算边界微元体的外力功计算设A点虚位移为**变形体虚位移原理和势能原理变形体虚位移时外力功计算边界微元体的外力功计算不管是否平衡均一样**变形体虚位移原理和势能原理变形体虚位移时外力功计算变形体的外力总虚功计算**变形体虚位移原理和势能原理变形体虚位移时外力功计算矩阵表示变形体的外力总虚功181748**变形体虚位移原理和势能原理虚位移原理的表述与证明虚位移原理的表述

受给定外力作用，变形连续体处于平衡状态的充分必要条件为：对任意虚位移（具有任意、独立性），外力所做的总虚功恒等于变形体所接受的总虚变形功，也即恒满足如下虚功方程26能说出虚位移原理和虚功原理的表述有何区别吗？**变形体虚位移原理和势能原理虚位移原理的表述与证明虚位移原理的必要性证明必要性需证明变形体平衡，虚功方程成立。15**变形体虚位移原理和势能原理虚位移原理的表述与证明虚位移原理的必要性证明必要性需证明变形体平衡，虚功方程成立。15=0格林公式**变形体虚位移原理和势能原理虚位移原理的表述与证明虚位移原理的充分性证明充分性需证明虚功方程恒成立，变形体必平衡。设变形体不平衡，每瞬时均考虑惯性力则动平衡。也即此时，“体积力”为：因为“平衡”，由必要性可得：**变形体虚位移原理和势能原理虚位移原理的表述与证明虚位移原理的充分性证明充分性需证明虚功方程恒成立，变形体必平衡。因为虚功方程恒成立，由此可得：因为虚位移的任意性和独立性，由此可得：这表明在虚功方程恒成立时变形体必无加速度**变形体虚位移原理和势能原理虚位移原理的表述与证明虚位移原理的几点说明1）适用于一切可变形物体（可变性固体、流体等）。2）虚功原理和虚位移原理是不同的。前者只是必要性命题，而后者则是充分必要的命题。3）王光远院士与我曾经证明，当不是取微元体进行研究时，不能证明变形体平衡。4）我们还曾经证明，当虚位移不具有完全任意和独立性时，也不能证明变形体平衡。**变形体虚位移原理和势能原理虚位移原理的表述与证明虚位移原理的几点说明5）只需将面积分改成体积分，线积分改成面积分，即可得到三维问题的虚功方程。6）利用虚位移原理做近似分析时，是应用原理的充分性，认为是用必要性时错误的。7）像虚功原理证明中一样，外力总虚功可分解成荷载与切割面内力的总虚功的和。此时格林公式实质是切割面内力总虚功为零。8）格林公式也可理解成是变形体虚功原理的变形。请大家自行考虑如何从虚功方程出发，用平衡和边界条件推得格林公式。——作业**变形体虚位移原理和势能原理一些名词含义的解释1）任何满足几何方程和位移边界条件的位移，称作可能位移，记作[d]k。2）由可能位移通过几何方程求得的应变，称作可能应变，记作[

]k。3）由可能应变通过物理方程求得的应力，称作可能应力，记作[

]k。4）可能应力在可能应变时所作的功，也即所储存的应变能，称作可能应变能，记作Uk。**变形体虚位移原理和势能原理一些名词含义的解释5）从可能位移退回到初始（也称自然）状态时，外力所作的功，称作外力势能，记作Pf。6）可能应变能和外力势能的总和，称作对应可能位移[d]k的总势能，简称总势能，记作

k。7）可能位移和真实位移的偏差，称作位移的变分，记作

[d]

。由此可得应变、应力的变分。8）可能位移总势能和真实总势能的偏差，其中与位移变分成线性关系的部分，称作势能的一阶变分，记作

。位移变分二次式部分称作势能的二阶变分，记作

2

。余类推。**变形体虚位移原理和势能原理势能驻值原理和最小势能原理势能驻值原理的表述：某一变形可能状态为真实位移状态的充分必要条件是，相应于此位移状态的变形体势能取驻值。也即势能对位移的一阶变分恒等于零。为了证明上述原理，先证明如下的格林公式：式中满足平衡条件，和间满足几何方程，还满足位移边界条件。**变形体虚位移原理和势能原理势能驻值原理和最小势能原理格林公式的证明：格林公式证毕**变形体虚位移原理和势能原理势能驻值原理和最小势能原理势能驻值原理的证明：变形体的总势能可表为：式中势能的一阶变分为：

将位移的一阶变分理解为虚位移，则由变形体虚位移原理的虚功方程可证势能一阶变分为零，能保证平衡。因此，势能原理结论正确。16**变形体虚位移原理和势能原理势能驻值原理和最小势能原理最小势能原理：线性、弹性变形体的总势能可表为：

由此可证明，对于一切位移变分，势能的二阶变分恒大于等于零（仅在位移变分为零时才等于零）。因此势能取最小值。

从势能原理证明可见，它和虚位移原理等价。都等价于平衡条件.**虚位移原理和势能原理的应用里兹法里兹法基本思路:

选取满足位移边界条件的函数作为“基函数”，将一个无限自由度的位移设为若干基函数的线性组合，

从而把无限自由度化为有限个自由度问题。

以所设位移作为可能位移，

令体系的总势能一阶变分恒等于零使系统近似平衡，

从而求得组合系数。

代回所设位移场，可进一步确定任意点的位移，利用几何、物理方程，还可求得应变和应力等。

上述近似方法即为里兹法。**虚位移原理和势能原理的应用里兹法里兹法的解题步骤:1）选取满足位移边界条件的函数作为“基函数”2）设近似位移场为基函数的线性组合3）将所设位移场代入势能表达式，从而将势能表为组合系数的函数4）令势能一阶变分（对组合系数偏导）为零，建立求组合系数的线性代数方程组，并求所需量**虚位移原理和势能原理的应用里兹法例一：试求图示悬臂梁挠曲线并计算B点挠度1）选取满足位移边界条件（A点挠度转角为零）的函数作为“基函数”2）对于所选的挠曲线，其虚位移、虚曲率，可能位移对应的弯矩等如下所示：**虚位移原理和势能原理的应用里兹法例一：试求图示悬臂梁挠曲线并计算B点挠度3）外力的总虚功为：4）总虚变形功为：**虚位移原理和势能原理的应用里兹法例一：试求图示悬臂梁挠曲线并计算B点挠度5）当v

仅取一项时，由虚功方程可得：当v

仅取二项时，由虚功方程可得：当v

取三项时，由虚功方程可得：**虚位移原理和势能原理的应用里兹法例一：试求图示悬臂梁挠曲线并计算B点挠度5）当v

仅取一项时，弯矩为：当v

仅取二项时，弯矩为：当v

取三项时，弯矩为：里兹法位移精度高于内力精度当试函数组合包含真解时，结果为精确解，否则为近似解**虚位移原理和势能原理的应用里兹法例二：试求图示桁架结点位移并计算各杆内力(各杆EA相同)设D点水平位移为u，竖向位移为v。在此位移下，体系的应变能U为：体系的外力势能为：体系的总势能为：**虚位移原理和势能原理的应用

根据势能原理，真实位移应使总势能最小，因此由势能对位移的偏导数为零可得由此可解得**虚位移原理和势能原理的应用

由位移可求得各杆变形（伸长）如下：由此可解得各杆的轴力为：从这个例子你能得到什麽结论？（可参考龙书14章）

**虚位移原理和势能原理的应用

例三：用原理推导等直杆AB杆端位移和杆端力关系

取试函数如下：**虚位移原理和势能原理的应用

例三：用原理推导等直杆AB杆端位移和杆端力关系

杆内任意一点的挠度v可用结点位移作参数，用试函数的组合来得到，从试函数示意图可见，此时位移边界条件自动满足

由此位移引起的杆件应变能U为：**虚位移原理和势能原理的应用

由此位移引起的杆件外力势能可用下式求得

由势能原理，总势能对位移的偏导数等于零（这都是数学推演，请大家自己看龙先生的教程）可得

将应变能和外力势能相加，可得杆件总势能为：

式中[k]、{R}见下一页。**虚位移原理和势能原理的应用

由势能原理所推得的“刚度方程”式中刚度矩阵为**虚位移原理和势能原理的应用

由势能原理所推得的“刚度方程”式中等效荷载矩阵为希望大家自行推导上述结果不难验证，此结果和由叠加得到的转角位移方程一样。杆端力矩阵等效结点力矩阵**势能原理在平面问题中的应用

设可能位移为

试用势能原理求图示平面应力问题的位移。显然此位移自动满足位移边界条件。当只取一项时

由此可得应变能为**势能原理在平面问题中的应用

由此可得应变能为

应力边界的位移为由此可得外力势能为**势能原理在平面问题中的应用

由此可得总势能为由最小势能原理可得**势能原理在平面问题中的应用

由此可求得代回可能位移表达式，可得位移为

进一步可以证明，当级数取多项时，其他系数全都等于零，由“可能位移”所求得的解满足平衡和应力边界条件，所以上述结果就是精确解。**18**高斯公式**

平面问题物理量的矩阵表示

引入两个算子矩阵微分算子矩阵方向余弦矩阵15**《结构力学基础》结束语《结构力学基础》到此结束。建议自学《教程》十四、十五、十六章，参考其他教材等作出笔记，从而进一步培养自己获取知识的能力。下学期的内容一定意义上比这学期难，难就难在它要求大家对本学期知识掌握的扎实，对线性代数、常微分方程等数学知识和算法语言有较好的功底。为了学好下学期内容，希望大家切莫学“熊瞎子掰苞米”，不要学一个扔一个。如果已经遗忘一些内容，希望自行复习补上！感谢大家这学期对改革试点工作的支持，请大家在考试前将你对试点工作的希望、意见、建议和体会等交一份书面材料，谢谢大家！预祝大家考试取得好成绩！**第三章杆系结构单元分析*第三章杆系结构单元分析引言等直杆单元的单元分析杆系结构单元分析的实质杆系结构单元分析子程序*3.1引言结点：杆件交汇点、刚度变化点、支承点。有时也取荷载作用点。图中1、2、3、4点均为结点。单元：两结点间的等直杆段。图中1-3、2-4、3-4为单元。编码：黑的结点编号称整体码。红的1、2局限于单元，称局部码。坐标：兰的坐标称整体坐标。红的x、y局限于单元，称局部坐标1342yxxy121122右手系*3.2等直杆单元的单元分析目的：像位移法一样，通过“一拆、一合”来解决结构分析。为此，必须首先掌握单元的特性。杆系最简单，由它介绍思想和方法容易掌握，可为以后学习奠定基础，因此必须深刻理解。3.2.1等直拉压杆结构中拆出的单元如图所示。1）广义坐标法设任意点位移为

u=a1+a2x广义坐标，边界条件只两个幂级数简单右手系ijxy12u1,F1u2,F2pEA，l*3.2等直杆单元的单元分析利用边界条件可得

a1＝u1；a2＝(u2-u1)／l将广义坐标代回

u=a1+a2x，整理后可得

u=(1-x／l)u1+u2x／l右手系ijxy12u1,F1u2,F2pEA，l2）形函数及性质形函数自然坐标本点处为1它点处为0

处总和为1任意点的位移可用形函数表为u=(1-x／l)u1+u2x／l=Ｎ1u1＋Ｎ２u２*轴力右手系ijxy12u1,F1u2,F2pEA，l３）用虚位移原理列式3-1）虚位移设结点虚位移为

ui

（i=1，2），则

u=Ｎ1

u1＋Ｎ２

u２3-2）外力虚功3-3）虚变形功3.2等直杆单元的单元分析*右手系ijxy12u1,F1u2,F2pEA，l3-4）用矩阵表示3.2等直杆单元的单元分析*ijxy12u1,F1u2,F2pEA，l3-5）单元刚度方程由虚位移原理可得引入如下矩阵：单元刚度矩阵单元等效荷载则单元刚度方程改写为3.2等直杆单元的单元分析*4）小结4-1）单元位移场可用“广义坐标法”建立。4-2）形函数“本点1，它点0，任意点总和1”。4-3）虚位移原理列式结果单元刚度方程为满跨均布轴力时3.2等直杆单元的单元分析*3.2.2等直杆扭转结构中拆出的单元如图所示。1）试凑法设任意点自然坐标为

，为满足“本1，它0”可设N1=1-

，N2=

。

=N1

1+

N2

2。由性质试凑得到右手系ijxy12

1,M1

2,M2mGJ，l2）势能原理列式2-1）外力势能2-2）应变能3.2等直杆单元的单元分析*右手系ijxy12

1,M1

2,M2mGJ，l2-3）总势能

=U+Pf2-4）势能原理列式结果为3）小结3-1）形函数可根据其性质通过试凑来建立。3-2）将总势能表为结点位移的函数，可由变分为零列式（偏导数为零）得到单元刚度方程。满跨均布扭矩时3.2等直杆单元的单元分析*3.2.3等直弯曲杆单元1）杆端位移（不计轴向变形）2）杆中任意点位移2-1）设挠度为2-2）由此可见，

N1是图示v1=1的挠曲线，因此由积分二次，利用边界条件可得3.2等直杆单元的单元分析*同理，可得3）势能原理单元列式3-1）应变能为则3-2）外力势能为3.2等直杆单元的单元分析*

在图示荷载作用下（坐标正向为正）外力势为3-3）总势能为3.2等直杆单元的单元分析*3-4）单元列式由势能原理可得单元刚度方程单元刚度矩阵经数学推导可得单元刚度矩阵元素如下等效结点荷载满跨均布时3.2等直杆单元的单元分析*3.2等直杆单元的单元分析单元刚度矩阵经数学推导可得单元刚度矩阵元素如下和转角位移方程结果相同29*3-5)小结3-5-1)杆件单元形函数也可由挠曲线微分方程求得。3-5-2)利用虚位移原理或势能原理列式所的结果和用叠加原理建立的杆端位移-力关系一样。3-5-3)所得单元刚度矩阵对称、奇异（1、3行1、3

列是相关的）。3.2.4考虑轴向变形直杆弯曲单元1)单元上任一点的位移1-1)单元上任一点的位移包含轴向和横向位移分量由纯弯单元建立由拉压单元建立3.2等直杆单元的单元分析*1-2)单元形函数由3.2.1和3.2.3组合得到2)应变能1-3)任意点位移式中2-1)应变由两部分引起：拉压和弯曲。3.2等直杆单元的单元分析*拉压弯曲4)单元列式

由应变能和外力势能可见，单元刚度方程可由拉压和弯曲上述已获得的结果组合得到。3)外力势能式中p、u向右为正；q、v向上为正；m和转角逆时针为正。2-2)应变能也包含两部分，拉压和弯曲互不藕联。3.2等直杆单元的单元分析*3.2等直杆单元的单元分析*同理可以组合得到单元刚度方程5)小结5-1)这种单元称作“梁柱自由式单元”，它也具有对称性、奇异性。5-2)单元刚度矩阵和等效结点荷载矩阵可由拉压和弯曲单元组合得到。下面多处在“变形不藕联”条件下，利用简单单元直接构造组合单元。3.2等直杆单元的单元分析*3.2.5有约束单元1)由“梁柱自由式单元”建立有约束单元刚度方程的原则为1-1)刚度矩阵中划去被约束位移对应的行和列。1-2)等效结点荷载矩阵中划去被约束位移对应的行。1-3)单元刚度方程形式不变。单元刚度方程2)连续梁单元

划去1、2、4、5行、列。3.2等直杆单元的单元分析*

单元刚度矩阵和单元等效荷载矩阵为

试自行写出单元刚度和单元等效荷载矩阵。3)简支单元

划去1、2、5行、列。满跨均布时4)小结

如约束能限制刚体位移，单元刚度矩阵非奇异。3.2等直杆单元的单元分析

式中MG

为固端弯矩，顺时针为正。*3.2.6计剪切变形和带有刚域的单元1)计剪切变形的单元（短粗杆）1-1)形函数用广义坐标、试凑和解挠曲线微分方程都比较困难，这里直接给出。3.2等直杆单元的单元分析参考《有限单元法基础》P.41-44*图和式中的

为对细长杆

趋于零，形函数同前。对形函数建立有兴趣的同学可用力法+位移计算来验证3.2等直杆单元的单元分析*1-2)用力法解超静定，求

ij时考虑剪切变形影响。由所建立的形常数、载常数，可利用叠加原理得到杆端力-位移的关系，它即单元刚度方程。1-3)用杆件模型计算高层框架-剪力墙结构时，代表剪力墙的杆件要考虑剪切变形。1-4)单元刚度、等效荷载矩阵的各元素见《教材》。2)带刚域的单元

高层框架-剪力墙结构用杆件模型分析时，与剪力墙连接的单元，要考虑刚域的影响。

这种单元的刚度特性，可由弹性段单元刚度方程，再考虑此段杆端位移、杆端力和单元杆端位移、杆端力间的变换关系来建立。结果见《教材》。3.2等直杆单元的单元分析*3.2.7平面单元分析总结3)虚位移原理、势能原理的结果完全相同。4)单元刚度矩阵对称。自由单元奇异，无刚体位移单元可逆。1)单元分析的关键是：建立单元位移场。2)单元位移场一般可用广义坐标法、试凑法建立。5)无藕联情况下，组合单元的单元特性可由简单单元组装得到。6)单元刚度方程和由形、载常数用叠加原理所得单元杆端位移-杆端力关系完全相同。7)单元上外力是平衡的（静力问题），但由位移求内力并不平衡，杆中内力以后解决。3.2等直杆单元的单元分析*3.2.8空间单元分析总原则

空间单元的单元特性由平面简单单元组装得到。1)空间拉压（桁架）单元

空间桁架单元与平面桁架单元相同。2)交叉梁（格栅）单元

杆端位移、杆端力如图所示

可由无轴向变形弯曲单元和扭转单元组合得到。

由于杆端位移和内力的编号顺序如图所示

因此必须注意，1、3（4、6）力-位移正向规定和弯曲单元的区别，故k13、k16、k46的符号与弯曲单元相反。同理，等效荷载元素的符号也要改变。3.2等直杆单元的单元分析123546*由此可得交叉梁单元的单元刚度矩阵为3.2等直杆单元的单元分析对称16*3)空间梁单元

杆端位移为三个坐标方向线位移和三个绕轴转角，杆端力为三个坐标向力和力偶矩矢。坐标正向为正。空间梁单元正向3.2等直杆单元的单元分析213456781011121212142536x、z平面梁单元12x、y平面梁单元214312扭转单元12与交叉梁单元一样，只要注意正向与图示平面单元的区别，即可根据图示编号由平面单元组合得到。作业试写出图示位移编号下的单元刚度矩阵*3.3杆系结构单元分析的实质3.3.1单元刚度矩阵的性质

空间单元的单元特性由平面简单单元组装得到。1)空间拉压（桁架）单元

空间桁架单元与平面桁架单元相同。2)交叉梁单元

杆端位移、杆端力如图所示

可由无轴向变形弯曲单元和扭转单元组合得到。3)空间梁单元

杆端位移为三个坐标方向线位移和三个绕轴转角，杆端力为三个坐标向力和力偶矩矢。坐标正向为正。

可由拉压单元、两个无轴向变形弯曲单元和扭转单元组合得到。**SubroutineElem_Stiff(······)

说明

Stiff=0.0!单元刚度清零

SelectCase(Type)Case(1)

平面杆系结构单元

Case(2)

空间杆系结构单元

CaseDefault

出错信息

EndSelectEndSubroutineElem_Stiff3.4杆系结构单元分析子程序3.4.1单元刚度总体设计3.4.2说明部分设计Integer,Intent(in)::···

入口整型参数

Real(8),Intent(in)::···

入口实型参数

Real(8),Intent(out)::···

出口实型参数

Real(8)::Work1,···Integer::i,j,k,···

实型和整型工作变量3.4杆系结构单元分析子程序3.4.3平面杆系结构设计SelectCase(Plane)Case(1)

平面桁架元素赋值

Case(2)

平面梁柱元素赋值

Case(3)······CaseDefault

出错信息

EndSelect3.4.4空间杆系结构设计SelectCase(Space)Case(1)

空间桁架元素赋值

Case(2)

空间梁柱元素赋值

Case(3)

交叉梁元素赋值

CaseDefault

出错信息

EndSelect3.4杆系结构单元分析子程序3.4.5有关单元等效结点荷载设计和进一步的考虑1)单元等效结点荷载设计同仿单元刚度。2)从各类单元刚度元素的计算，可看到要用到长度、单元弹性特性、单元截面特性等数据。因此，要确定存放它们的数据结构。要将它们作为出口。3)为计算单元等效结点荷载元素，首先要建立各种荷载情况等效荷载表达式，它们可由积分或载常数表得到。然后要解决荷载信息的存放结构，也要将它们作为出口量。4)单元刚度矩阵、等效结点荷载矩阵都应先清零。4.1杆系结构整体分析

首先就全刚结点平面刚架进行讨论，然后推广。4.1.1总的思路

在单元特性搞清后，将单元拼装回去。在结点处位移自动协调基础上，如果全部结点平衡，则求得的结点位移将是实际结构的解。因此，整体分析就是设法建立结点平衡方程。4.1.2坐标转换

组成结构的杆件可以各个方向，单元分析对局部坐标，因此，必须将物理量转为统一坐标——整体坐标。1)力的转换关系4.1杆系结构整体分析2)位移转换关系3)转换矩阵

转换矩阵是正交矩阵。4.1杆系结构整体分析4)杆端力转换5)杆端位移转换6)刚度方程的转换

如果记称为整体单元刚度矩阵

则

这就是整体坐标下的单元刚度方程。本节以后的讨论认为都是对整体坐标的4.1杆系结构整体分析4.1.3结点平衡方程的建立1)一简单例子（如图）

图中有两套编号，红的是单元杆端编号，黑的是结构整体编号。1-1)结点示意121221①②③

图中蓝色的表示结点荷载（已知），红色的表示杆端力（未知的），、分别1、2单元杆端力子矩阵。对1、4结点“荷载”含有未知反力。21-2)结点平衡

由示意图可见，结构结点的平衡方程为4.1杆系结构整体分析

从例图可见，其全部结点平衡方程为121221①②③若记24.1杆系结构整体分析式中[I]、[0]分别为单位和零矩阵。若引入矩阵记号，则结点平衡方程可改写作

这一结论虽然是由一个例子得到的，但是显然对一切结构都是成立的。问题在于不同结构，[A]矩阵是不同的。4.1杆系结构整体分析4.1.4杆端位移用结点位移来表示121221①②③仍以简单例子来说明若记

由结点、杆端位移的协调条件，可得[

]、[

]的对应关系为

式中[A]T是前面力关系[A]的转置，因此[A]T称为位移转换矩阵。4.1杆系结构整体分析4.1.5整体刚度方程——结点平衡121221①②③若记

引入位移转换关系，则

这就是整体刚度方程，它的物理实质是结点平衡。[K]称作结构刚度矩阵（或整体刚度矩阵），[R]称作综合等效结点荷载矩阵，它由两部分组成。单元个数4.1杆系结构整体分析4.1.6整体刚度矩阵的建立121221①②③

若将[A]按单元分成图示三个子矩阵

则

由此可见，整体刚度矩阵可由各单元整体刚度矩阵装配累加得到。为说明如何装配，先将单元刚度矩阵进行分割整体结点码

则由矩阵乘法可证明，[A]i[k]i[A]iT的结果是，将刚度矩阵子矩阵按整体结点码r、s

送整体刚度矩阵相应位置。这一装配规则称为“对号入座”。4.1杆系结构整体分析1)任意结构情况

上面结论是通过具体例子（全刚结点平面刚架）得到的，由虚位移原理或势能原理进行整体分析（见讲义），可得任意结构其结论同此例。2)结点位移编号

如果按结点顺序，对结点非零位移进行依次编号，这一序号称作结点位移码。为便于计算机处理并减少结构刚度矩阵的阶次，将零位移的号码变为零。

对图示三铰刚架，当仅用一种单元（梁柱自由是单元）时结点位移编号如图所示。3)单元定位向量

按单元局部结点码顺序，将结点位移码排成的向量，称作单元的定位向量。①②③④①②③④4.1杆系结构整体分析

对图示刚架各单元的定位向量为①(0,0,1,3,4,5)②(0,0,2,10,11,12)③(3,4,5,6,7,8)④(6,7,9,10,11,12)①(0,0,1,2,3)②(0,0,,6,7,8)③(1,2,3,4,5)④(4,5,6,7,8)4)按单元定位向量集装刚度矩阵和综合荷载

前面说明的是分块子矩阵集装,下面说明如何按定为向量来集装.

如果如图所是采用各种不同的单元（一端有铰），则定位向量为①②③④①②③④如何获得带铰的单元刚度矩阵和等效荷载矩阵定位向量①②③④①②③④4.1杆系结构整体分析4-1)刚度集装

以3单元为例来说明定位向量单元局部位移码

根据单元局部位移码和定位向量的对应关系用定位向量位移码送元素。

根据单元局部位移码和定位向量的对应关系用定位向量位移码送元素，定位向量元素为零时不送。①②③④①②③④4.1杆系结构整体分析4-2)荷载集装

以4单元为例来说明定位向量局部位移码此结论同样适用于刚度集装4.1杆系结构整体分析4.1.7整体分析总结1)对局部坐标和整体坐标不一致的单元，要对刚度、荷载进行坐标转换。2)需对“结构”进行结点、位移的局部和整体编号。4)集装所得整体刚度矩阵是对称、带状稀疏矩阵，当支撑条件能限制刚体位移时，矩阵非奇异。3)根据单元局部位移码和定位向量的对应关系用定位向量位移码送元素，定位向量元素为零时不送。据此可集装、累加得到整体刚度矩阵。5)综合荷载由两部分组成，因此首先要将直接作用结点的荷载按结点位移码送入，如果还有单元等效荷载，再按定位向量集装、累加。4.1杆系结构整体分析8)如果有某位移码方向弹性支撑，需进行将弹簧刚度送入位移码对应的对角线元素位置累加。9)如果有某位移码方向已知支撑位移，需进行将“边界条件处理”。具体做法以后介绍。7)整体刚度方程实质是全部结点的平衡条件。6)刚度矩阵带状稀疏，其带宽取决于结点、位移编码。最大半带宽=定位向量中最大元素差+1。4.1杆系结构整体分析2.5.8边界条件的处理10)当用虚位移或势能原理作整体分析时（势能为例），应变能为单元应变能之和，外力势能为单元外力势能之和+结点外力势能。全部杆端力的外力势能彼此抵消。1)乘大数法

设第i个位移为已知值a，N=108或更大的数。乘大数法是将刚度矩阵Kii改为N

Kii，将Ri改为

N

a。请考虑为什麽这样做能使边界条件得到满足？2)置换法（划零置1）

设第i个位移为已知值a。4.1杆系结构整体分析

上述置换工作量大一些，显然可看出边界条件得到精确满足。4.1杆系结构整体分析★3)关于斜边界的处理

如图示意的斜支座情况，有多种处理方案。3-1)通过单元的坐标转换来处理xy

图示有斜支座单元，r结点处以倾角

-

来进行坐标转换，也即在r结点处整体坐标为图示xy。r3-2)通过增加一个单元来处理

图示有斜支座单元，r结点处沿y方向增加一个刚结的单元，此单元有“无穷大”的抗拉刚度、但没有抗弯刚度。单元长度可任意。3-3)对整体刚度矩阵进行处理（参见匡文起教材）最大半带宽2.5杆系结构整体分析2.5.10总刚度矩阵的存储与对应解法

因为总刚度矩阵对称、带状稀疏，利用这一特点可减少存储、提高解算效率。零元素零元素非零元素最大半带宽主对角线元素等半带存储零元素非零元素变带宽一维存储带宽是变的到P.552.5杆系结构整体分析

目前一般都用变带宽存储，下面结合程序说明存储和解法。首先介绍一些F90的语法。

定义导出类型导出类型——结点type::typ_Jointreal::x,y

！坐标

integer::GDOF(3)！整体位移码endtypetyp_Joint1)有关F90语法导出类型

新特性

用结点导出类型作为成员导出单元类型:type::typ_Element

integer::JointNo(2)

！结点编号

type(typ_Joint)::Node(2)

！结点

integer::GlbDOF(6)!定位向量

real::EA,EIendtypetyp_Elementtype(typ_Element)::Elem(5)

！定义5个单元类型…!对单元i的端点j的x,y,GDOF(1:3)的赋值Elem(i)%Node(j)=Joint(Elem(i)%JointNo(j))…由导出类型定义新类型由导出类型定义变量real::A(5),B(5,10),C(5)B=0.0

!对B清零A=1.0

!对A赋1:A(i)=1.0,i=1,5C=A+2

!数组与标量运算:A(1:5)+(/2,2,2,2,2/)A=C+A

!数组与数组运算（同形）C=sqrt(A)!数组的函数运算:C(i)=sqrt(A(i),i=1,5数组内部函数：dot_product(vector_a,vector_b)

!点积如：dot_product((/1,2,3/),(/2,3,4/))的值为20（待续）有关F90语法数组运算与赋值：matmul(matrix_a,matrix_b)

矩阵相乘如：locEDisp=matmul(T,glbEDisp)transpose(matrix)

矩阵转置如：glbEDisp=matmul(transpose(T),locEDisp)size(array,dim)

求数组第dim维的长度dim为可选变元:size(a,dim=2)若array为一维时，可不用dim。sum(array,dim,mask)

数组元素求和dim,mask为可选变元；mask=条件表达式sum(a(1:10))

对a的1到10元素求和sum(a(1:10),mask=a>0)

对a(1:10)中大于0的元素求和（续）有关F90语法where结构

新特性例where(C>0)C=0A=B*Dendwherewhere(C>0)A=Bendwhere定义where(数组关系表达式)

数组赋值语句

…elsewhere

数组赋值语句

...endwhere规则:

1)同形数组;2)不许嵌套;3)最多二分叉有关F90语法cycle和exit语句

新特性

用在do循环中

cycle

——

作下一个循环步

exit

——

跳出循环，执行enddo后一条语句等效例

do...if(.not.cond1)then...if(cond2)goto5...endifenddo5...用法do

...if(cond1)cycle ...

if(cond2)exit...enddo...有关F90语法数组构造函数spread语法spread(数组名,dim,ncopies)

将数组沿dim维方向复制ncopies形成新数组

dim,ncopies

—

整型、位置变元、关键字变元（若按位置引用，可略关键字）例:

（仅限一维数组）1）

spread(one,dim=1,ncopies=3)spread(one,1,3)spread(one,ncopies=3,dim=1)

[1,1,1]

或[1,1,1]T2）

ELocVec(1:6)=(/1,0,3,4,5,0/)spread(ELocVec,dim=1,ncopies=3)3）

spread(A(2,2:),dim=1,ncopies=2)如果dim=2呢?有关F90语法指针pointerpointer是变量的属性，可以指向相同类型的变量;被指向的目标必须具有target属性或pointer属性可以将指针变量理解为别名、称号real,target::a,b,EDisp(6)

!可被指针所指

real,pointer::p1,p2!称号：班长、课代表！p1,p2是指针，可以指向实型数据real,pointer::G(:)!先进集体！G是指针，可以指向一维实型数组指针是一种“称号”，上述声明语句建立了“称号”，但并未“授予”某个变量这个称号，因此是指向“空”，并未占用内存a=3.0p1=>a!p1指向a!称号p1授予a，a的数据有两个名：固定名a和流动名p1;既可用p1也可用a（p1—

班长，a—

张三）a=4.0!a的值变为4.0,p1也变为4.0p1=>b!班长换人了G=>EDisp!先进集体有了得主！EDisp(:)的长度就是G(:)的长度，用G和用EDisp同样效果又如：real,target::a,breal,pointer::p,qa=3.14b=2.0p=>a!p=a=3.14q=>b!q=b=2.0q=p!(q指向的数据b)=(p指向的数据a)！即：所有=3.14指针可以指向有名的数据区，也可以指向无名的数据区real,pointer::b(:,:)integer::nread(*,*)nallocate(b(n,n))

!指针指向了一个刚开辟的数组!以下可以当作常规数组用b(1,1)=1.1b(1,2)=1.2...deallocate(b)有关F90语法用指针建立动态数组指针与allocatable数组的区别具备allocatable数组的所有功能还可以用在导出类型中，例如整体刚度矩阵的变带宽存储：type::typ_Kcol

!整体刚度矩阵K的列

real(8),pointer::row(:)

!该列的行元素endtype...type(typ_Kcol),allocatable::Kcol(:)…allocate(Kcol(NGlbDOF))!分配了NGlbDOF列...allocate(Kcol(5)%row(3:5))!第5列只用3至5行(1)NElem,NJoint,NGlbDOF,NJLoad,NELoad

单元数结点数总自由度数结点荷载数单元荷载数[Joint-结点]…NJoint行(2)Joint%X,Joint%Y,GDOF(1:3)X坐标Y坐标结点位移码[Elem-单元]…NElem行(3)JointNo1,JointNo2,EA,EI

起点号终点号刚度[JLoad-结点荷载]…NJLoad行(4)JointNo,LodDOF,LodVal

作用点号局部位移码荷载值[ELoad-单元荷载]…NELoad行(5)ElemNo,Indx,Pos,LodVal

单元号类型号位置荷载值Indx

类型

pos1均布长度2集中位置3...2)某程序输入数据说明3,5,7,1,10,0,0,0,00,4,1,2,34,4,4,5,64,4,4,5,74,0,0,0,01,2,1.0E9,162,3,1.0E9,245,4,1.0E9,122,1,10.0E3

!结点2，自由度1，值为10E32,1,4,-4.0E3

！单元2，均布，长4米，值为-4E32-1)数据文件例子：(2)(1)(3)24135i=6i=4i=310kN4kN/m4m4mEA=109N(1)(2)坐标位移码(3)(4)(5)结点号EA,EIread(5,*)NElem,NJoint,NGlbDOF,NJLoad,NELoadallocate(Joint(NJoint))allocate(Elem(NElem))allocate(JLoad(NJLoad))allocate(ELoad(NELoad))...read(5,*)(Joint(i),i=1,NJoint)read(5,*)(Elem(ie)%JointNo,Elem(ie)%EA,&ELem(ie)%EI,ie=1,NElem)if(NJLoad>0)read(5,*)(JLoad(i),i=1,NJLoad)if(NELoad>0)read(5,*)(ELoad(i),i=1,NELoad)2-2)程序读入计算所需数据：2-3)求始行码和分配带宽子程序!==================================subroutineSetMatBand(Kcol,Elem)!接口简单!==================================

type(typ_Kcol),intent(inout)::Kcol(:)

！总刚列

type(typ_Element),intent(in))::Elem(:)！单元

integer(ikind)::minDOF,ELocVec(6)integer(ikind)::Row1(size(Kcol,dim=1))

！Row1为自动数组，子程序结束后自动释放。！这样做可使接口简单，减少了数组的控制变量。

integer(ikind)::ie,j,NGlbDOF,NElemNGlbDOF=size(Kcol,dim=1)！使接口简单NElem=size(Elem,dim=1)Row1=NGlbDOF

!先设所有始行码同终行码

!确定各列始行码向量

doie=1,Nelem!对单元循环

!确定定位向量

ELocVec(:)=Elem(ie)%GlbDOF(:)!寻找定位向量中大于零的最小值

minDOF=minval(ELocVec,mask=ELocVec>0)

!屏蔽定位向量中小于等于零的

where(ELocVec>0)

!寻找Row1(ELocVec)和minDOF中的最小值

Row1(ELocVec)=min(Row1(ELocVec),&minDOF)endwhereenddo!为各列的带宽分配空间

doj=1,NGlbDOF

!对总自由度数循环

allocate(Kcol(j)%row(Row1(j):j)

!给Kcol(j)%row从Row1(j)到j个空间

Kcol(j)%row=zero

!总刚元素清零

enddoreturnendsubroutineSetMatBand3)整体刚度矩阵的集成

doie=1,NElem…!计算单刚

EK(6,6),ELocVec(6)doj=1,6!对单元逐列集成

JGDOF=ELocVec(j)!取出位移码

if(JGDOF=0)cycle!作下一循环步

where(ELocVec>0.And.ELocVec<=JGDOF)

!位移码非零同时小于第j个局部码对应的位移码

Kcol(JGDOF)%row(ELocVec)=&Kcol(JGDOF)%row(ELocVec)+EK(:,j)!集成

endwhereenddoenddo局部码位移码4)变带宽矩阵的分解求解4-1)LDLT分解法求解[A]{x}

={b}若对称正定，则可分解为=单位上三角阵对角阵原方程变为求解步骤：1.分解：2.解y：3.解x：LDLT分解法Gauss消去前消去处理右端项回代(不同的b只做一次分解)存放：主对角—

上三角—

时不必求和(上三角:i<j)不动存在处存在处(第j列系数)上三角:i<j

4-2)分解一般情况：逐列分解对角线:i=

j>1(第j列系数)第

i列中第1个非零元素的行码为：4-3)变列宽存贮的分解修正：

第

j

列中第1个非零元素的行码为：

则分解顺序4-4)F90实现!三角分解diag(1:ncol)=(/(Kcol(j)%row(j),j=1,ncol)/)doj=2,ncolrow1=lbound(Kcol(j)%row,1)!i_1jdoi=row1,j-1row_1=max(row1,lbound(Kcol(i)%row,1))!i_1s=sum(diag(row_1:i-1)*Kcol(i)%row(row_1:i-1)*&Kcol(j)%row(row_1:i-1))

!求和部分

Kcol(j)%row(i)=(Kcol(j)%row(i)-s)/diag(i)enddo

！第j列系数分解完毕

s=sum(diag(row1:j-1)*Kcol(j)%row(row1:j-1)**2)diag(j)=diag(j)-s

!第j列的主对角元素enddo4-5)

一般情况解

y公式可不动可存在处4-6)变列宽存储解

y的修正上三角的第

i列从第行元素开始!回代步骤1:自上而下doi=2,ncolrow1=lbound(Kcol(i)%row,dim=1)GP(i)=GP(i)-sum(Kcol(i)%row(row1:i-1)&*GP(row1:i-1))enddo

求出后，不再用，可将存在处4-7)一般解

x的公式自下而上：这样做的缺欠：自下而上逐行计算，行中遇到0元素需跳过，不方便！4-8)改为自右向左逐列计算！记解出对y向量预处理再解出第n列乘上后移到右边去，修正y向量!回代步骤2:自右向左GP(:)=GP(:)/diag(:)doj=ncol,2,-1row1=lbound(Kcol(j)%row,dim=1)GP(row1:j-1)=GP(row1:j-1)-GP(j)*&Kcol(j)%row(row1:j-1)enddo!现在GP就是解4-9)小结：

无须一维存贮，无须行列码转换定位公式与程序直接对应翻译，直截了当求和采用内部函数进行数组运算动态内存，用多大、开多大数据封装性好，接口简单:

subroutineSetMatBand(Kcol,Elem)subroutineVarBandSolv(Disp,Kcol,GP)

通用性强2.6杆系结构内力计算2.6.1杆端内力计算公式

解方程的结果可以得到结点位移，有了位移就可以进一步求各单元的内力。

解算步骤：

从[

]根据定位向量取出单元杆端位移。

由单元倾角确定是算还是。

减去等效结点荷载或加上固端力矩阵。2.6杆系结构内力计算2.6.2杆中任意截面内力计算公式

注意事项：

由图示隔离体图，列杆段的平衡方程即可得到任意截面的内力计算公式。请自行写出。

按上述隔离体图所求得的内力是精确的。

求得单元结点位移后，代入单元位移场、求应变、求内力，这样所得的内力一般不满足平衡条件。只是近似解。2.7杆系结构静力分析程序2.7杆系结构静力分析程序程序演示本程序可作以下结构计算:平面和空间桁架计算(网架视作空间桁架)平面和空间刚架计算各种组合结构计算高和不等高三铰拱计算多跨梁(静定、超静定）计算2.8程序调试中关键变量的

速算方法2.8.1总刚度矩阵元素的确定1)总刚度矩阵元素的物理意义

整体刚度方程为如果，则可见总刚度矩阵元素的物理意义为：当且仅当时，在处所需施加对应于的广义力。或理解为：当且仅当时，在限制位移的约束上所产生的约束反力。2.8程序调试中关键变量的

速算方法2)指定总刚度矩阵元素的速算方法

根据总刚元素的物理意义，令仅仅产生,利用位移法中的形常数作弯矩图，象位移法一样即可求得指定总刚元素值，为校核总刚集成的正确性提供测试数据。

注意：（1）要牢记总刚元素的物理意义。（2）仅仅产生。（3）实质上这里纯粹是用位移法来求解。2.8程序调试中关键变量的

速算方法2.8.2综合等效荷载元素的确定1)综合等效荷载元素的组成

综合等效荷载为也即，它由直接结点荷载和等效结点荷载组成。2)综合等效荷载元素的确定

直接结点荷载只需将外荷载坐标方向投影即可，因此关键是确定等效结点荷载。

由2.2.7之6)已知，单元刚度方程和象位移法用叠加所得力-位移关系（转角位移方程）一样，因此2.8程序调试中关键变量的

速算方法只要熟记载常数，即可将单元荷载转化为单元等效结点荷载，再经过往坐标方向的投影，即可获得作用在结点的元素。由此和直接结点荷载相加即得到需求的综合等效荷载元素。

注意：（1）坐标正向的“荷载”为正。（2）建议先按载常数确定固端力的实际方向和数值,然后反方向得到等效荷载实际方向(局部坐标方向)。（3）将所有单元荷载的等效荷载作用到结点，同时考虑直接结点荷载（斜杆需投影）即可得需求值。2.8程序调试中关键变量的

速算方法2.8.3单元杆端内力元素的确定1)单元杆端位移的确定

整体刚度方程求解结果，所得到的是整体坐标下的结点位移，为求单元杆端内力，需作两件事：

从整体位移矩阵中根据定位向量驱除单元位移；将整体坐标的位移往单元局部坐标方向投影。这样即可获得单元局部坐标下的单元杆端位移。2)由单元刚度方程来求“固端力”2.8程序调试中关键变量的

速算方法

如果只需求某指定内力，实际并不需要作整个矩阵乘。

按杆端力方程求要作矩阵运算，为避免它，可在获得局部坐标位移后，利用形、载常数通过叠加来得到某指定内力值。3)由形、载常数叠加来求

注意：（1）如果要求整体坐标下的内力该怎麽办？（2）“固端力”符号规定和位移法有区别。（3）建议用叠加法求。2.9几点重要说明（1）本章方法、思路具有普遍性。特别是整体分析，其方法、结论完全适用于其他有限元分析。（2）为用有限元分析实际结构，首先要做离散化：建立两套坐标、确定结点、单元、位移编号。此时要注意尽可能使半带宽最小。（3）有限元分析的关键问题是：建立合适的位移模式。对一般问题可用广义坐标法或试凑法。对杆系问题也可由挠曲线微分方程积分得到形函数。（4）可用虚位移原理或势能原理进行单元分析。（5）可用虚位移原理或势能原理进行整体分析.结论是：整体刚度矩阵、综合等效荷载可按定位向量由单元集装得到。“综合=直接+等效”。实质是结点平衡。

三、分支稳定和增量变刚度极限分析及程序

三、分支稳定和增量变刚度极限分析及程序结构稳定性分析基本概念分支稳定分析程序结构极限分析基本概念增量变刚度极限分析程序3.1结构稳定性分析基本概念3.1.1一些基本概念薄壁、高强、受压结构，设计不当容易产生部件或整个结构丧失稳定。因此，结构设计除关心强度、刚度外，对易失稳的结构还要进行稳定验算。结构稳定分静力和动力稳定两大类，本章只讨论静力稳定。一些问题可以抽象为受压杆都是“理想中心受压直杆”计算模型，这称为完善体系。如果结构杆件不满足上述“理想中心受压直杆”假定(不直或有偏心)，此系统称非完善体系。完善体系从稳定到不稳定，其受力、变形状态将变化，也即随荷载变大有分叉点，称分支点稳定。3.1结构稳定性分析基本概念

非完善体系，一般受力、变形性质不发生改变。但随着荷载增大存在一极值荷载（此后变形增大荷载反而减少），这类稳定现象称极值点稳定。一些扁平拱式结构还可能产生从受压跳转到受拉的“急跳”现象，当然实际结构不允许出现这情况。本章以讨论分支点稳定问题临界荷载为主，也介绍一点其他内容。由于实际结构刚度都很大，变形和杆件尺寸相比十分微小，因此作受力分析列平衡方程时都忽略变形影响。因此线弹性材料力-位移成正比，叠加原理适用。

在作稳定分析时，必须考虑变形的影响，这时叠加原理不再适用。SMCAI3.1结构稳定性分析基本概念3.1.2稳定问题分析基本方法一：静力法通过考虑失稳状态下的平衡关系，利用两类稳定问题的特征，确定临界荷载的方法——静力法。1)分支点稳定问题1-1)分析步骤设定约束所允许的可能失稳状态建立平衡方程用分支点稳定的平衡两重性（可在两状态平衡）建立特征方程求特征方程的非零解，从而得到临界荷载。1-2)按大、小挠度分析及结论待切换到SMCAI。SMCAI3.1结构稳定性分析基本概念3.1.3稳定问题分析基本方法二：能量法通过考虑失稳状态下的应变能、外力时能，利用稳定问题的能量特征（总势能取驻值），确定临界荷载的方法——能量法（只讨论分支点问题）。1)分析步骤设定约束所允许的可能失稳状态通过求应变能、外力时能确定总势能用分支点稳定的能量准则（总势能取驻值）建立特征方程求特征方程的非零解，从而得到临界荷载。2)举例待切换到SMCAI。SMCAI《结构力学教程》和SMCAI上还有许多内容限于授课学时希望学有余力的同学能自学。3.2分支稳定分析程序

用有限单元法可以解决两类稳定问题，还可以考虑弹性后的屈曲问题。但这些有待大家今后继续学习，本节只介绍弹性分支点稳定问题。这等于假定结构直到临界状态材料是处于弹性阶段，荷载达临界值以前结构处于无弯矩状态、可忽略轴向变形影响。3.2.1基本原理1)刚度方程的建立在第二章建立单元刚度方程时，不计杆内压力对