计算结构力学-有限元

1、计算结构力学计算结构力学课课 件（二）件（二）第第2章章 有限单元法有限单元法2.1 概述概述2.2 弹性力学平面问题的矩阵描述弹性力学平面问题的矩阵描述2.3 三结点三角形单元分析三结点三角形单元分析2.4 高精度三角形单元、矩形单元高精度三角形单元、矩形单元2.5 C0连续型单元形函数的构造连续型单元形函数的构造2.6 平面等参数单元分析平面等参数单元分析2.7 有限元程序实现有限元程序实现2.8 平面杆件结构有限元平面杆件结构有限元2.9 板弯矩的有限元板弯矩的有限元2.1 概述概述2.1.1 发展概况发展概况 Courant1943年应用三角形分片插值函数年应用三角形分片插值函数; T

2、urner, Clough等等1956年推广直接刚度法年推广直接刚度法; Clough1960年提出年提出“有限单元法有限单元法”名称；名称； Zienkiewicz等等编写第一本有限元方面专著；编写第一本有限元方面专著； Melosh 证明有限元位移法是里兹法另一形式；证明有限元位移法是里兹法另一形式； 冯康冯康独立证明了有限元法；独立证明了有限元法； Wilson第一个编写通用有限元软件：第一个编写通用有限元软件：SAP； 1960-70年代理论基础研究年代理论基础研究 1960至今：实际工程应用、复杂问题理论研究至今：实际工程应用、复杂问题理论研究 2.1 概述概述 1960年代前：起步

3、阶段年代前：起步阶段 1960年代年代-70初：理论体系建立，快速发展期；初：理论体系建立，快速发展期； 1970初初-80中：巩固期；中：巩固期； 1980中之后：推广、综合应用；中之后：推广、综合应用；2.1.1 发展概况（续）发展概况（续）通用有限元软件：通用有限元软件：SAP, ADINA, NASTRAN, ANSYS, ABAQUS, MIDAS等等2.1.2 有限单元法概念有限单元法概念1) 离散化离散化 划分为有限数目的单元；划分为有限数目的单元； 单元间在指定点连接单元间在指定点连接结点。结点。 单元形状、连接方式可不同。单元形状、连接方式可不同。rg连续体 平面问题的常用单

4、元：平面问题的常用单元：三结点三角形单元三结点三角形单元六结点三角形单元六结点三角形单元矩形单元矩形单元任意四边形单元任意四边形单元8结点曲边结点曲边四边形单元四边形单元2) 单元分析单元分析 假设位移模式，分析单元力学特性：假设位移模式，分析单元力学特性： Fe=ke 体力、面力体力、面力 等效结点荷载等效结点荷载3) 整体分析整体分析 建立建立F=K K= R4) 再次单元分析再次单元分析 求出各单元的应变和应力。求出各单元的应变和应力。静力等效静力等效2.2.1 两类平面问题两类平面问题1) 平面应力问题平面应力问题 z=0, zx=0, zy=0 x= x(x, y)y= y(x, y

5、)xy= xy(x, y)例如：深梁、剪力墙例如：深梁、剪力墙(受竖向力受竖向力)2.2 弹性力学平面问题的矩阵描述弹性力学平面问题的矩阵描述xt/2zOyt/2y体积力体积力表面力表面力平面应力问题平面应力问题形状、受力沿形状、受力沿z向不变向不变 + z向很薄向很薄平面应力问题平面应力问题z向很长向很长平面应变问题平面应变问题2) 平面应变问题平面应变问题 w=0, z=0, zx=0, zy=0 x= x(x, y) y= y(x, y) xy= xy(x, y) x= x(x, y) y= y(x, y) xy= xy(x, y)例如：挡土墙、重力坝例如：挡土墙、重力坝xOyz平面应变

6、问题平面应变问题 1) 基本量基本量2.2.2 基本量及基本方程的矩阵表示基本量及基本方程的矩阵表示TxyyxTyxqqqTxyyxTvufTyxppp应力应力:应变应变:位移位移:体积力体积力:表面力表面力: 几何方程：几何方程： 物理方程：物理方程： (平面应力问题平面应力问题) D: 弹性矩阵弹性矩阵; E: 弹性模量弹性模量; : 泊松比。泊松比。TTxyyxyuxvyvxu2) 基本方程基本方程210001011,2EDD平面应力平面应力 平面应变平面应变 112EE2100010112ED 平衡方程弱形式平衡方程弱形式 能量原理能量原理 VTSTVTVSqfVpfdddSTVTVT

7、SqfVpfVddd21ii) 变分方程变分方程 势能驻值原理势能驻值原理 =(U+UR)=0i) 虚功原理：外力虚功虚功原理：外力虚功=内力虚功内力虚功AqBxy在集中力在集中力F作用下作用下,F=U1 V1 U2 V2 Un VnT =u1 v1 u2 v2 un vnT yxtFTTddd21FVTVT总势能简化：总势能简化：虚功方程简化：虚功方程简化： 单元局部编号：单元局部编号：i, j, m 单元结点位移向量单元结点位移向量: e=ui vi uj vj um vmT 单元结点力向量单元结点力向量: Fe=Ui Vi Uj Vj Um VmT2.3 三结点三角形单元分析三结点三角形

8、单元分析yui , (Ui)ixvi , (Vi)ujjvjummvm2.3.1 结点力和结点位移结点力和结点位移 整体结点位移向量整体结点位移向量: =u1 v1 u2 v2 un vnT 整体结点力向量整体结点力向量: F =U1 V1 U2 V2 Un VnT, n 结点总数结点总数.2.3.2 单元位移模式单元位移模式 1) 什么是位移模式（位移函数）什么是位移模式（位移函数） 利用单元的结点位移将整个单元的位移分量利用单元的结点位移将整个单元的位移分量表示为坐标的函数。表示为坐标的函数。 yxP(x, y)ijm2) 三结点三角形单元的位移模式三结点三角形单元的位移模式 设：设： u

9、= 1+ 2 x+ 3 y v= 4+ 5 x+ 6 y系数系数 1 6由结点位移由结点位移 ui , vi , uj , vj , um , vm确定确定.将位移模式写成结点位移的显式：将位移模式写成结点位移的显式： u= Niui+ Njuj +Nmum v= Nivi+ Njvj +NmvmNi、Nj、Nm：形函数形函数 (插值函数插值函数)yxPijm),(,111111轮换mjiyxyxyxyxyxyxNmmjjiimmjjiNi(x, y)ijm1ai、bi、ci：分母行列式第：分母行列式第1行元素代数余子式：行元素代数余子式： ai=xjym xmyj bi=yj ym （i,

10、j, m轮换）轮换） ci= xj+xm ),()(21111111轮换mjiycxbaAyxyxyxyxyxyxNiiimmjjiimmjjiNi(x, y)ijm1 形函数的性质形函数的性质: (i) (Ni )i=1，(Ni )j=0，(Ni )m=0 (ii) 单元内任一点：单元内任一点：Ni+Nj+Nm=1emmjjiimmjjiiNvNvNvNuNuNuNvuf3) 位移模式的矩阵表示位移模式的矩阵表示形函数矩阵形函数矩阵mjimjiNNNNNNN000000yxPijm2.3.3 单元应变矩阵、应力矩阵单元应变矩阵、应力矩阵应变矩阵应变矩阵：B= Bi Bj Bm ,)3 , 2

11、 , 1(0021ibccbABiiiiif =N e得到得到代入代入几何方程几何方程 =B ebi=yj ymci= xj+xmB元素均为常数元素均为常数 常应变单元常应变单元 =B e得到得到代入代入物理方程物理方程 =S e应力矩阵应力矩阵：S=DB)3 , 2 , 1(2121)1(22ibccbcbAESiiiiiiibi , ci为常数为常数S元素均为常数元素均为常数 常应力单元常应力单元S= Si Sj Sm ，对平面应力问题：对平面应力问题：2.3.4 有限元方程的建立有限元方程的建立得到得到代入代入 (U+UR)=0有限元方程有限元方程位移模式位移模式f =N e结构离散结构

12、离散eSTTeeATTeeAeTTeSTVTVTeeestqNyxtpNyxtBDBsqfVpfVd)(dd)(dd)(21ddd21令：e=G， G 62n位置矩阵 eATeyxtBDBkdd 单元刚度矩阵单元刚度矩阵eeSTATeqepestqNyxtpNRRRddd 单元等效结点荷载向量单元等效结点荷载向量eeTTeeTTRGGkG)(21eSTTeeATTeeAeTTeeeestqNyxtpNyxtBDBd)(dd)(dd)(21令： 整体刚度矩阵整体刚度矩阵 整体等效结点荷载向量整体等效结点荷载向量eeTTeeTTRGGkG)(21eeTGkGKeeTRGR21RKTT00RK有限元

13、方程有限元方程2.3.5 单元刚度矩阵单元刚度矩阵单刚列式单刚列式mmmjmijmjjjiimijiiTATkkkkkkkkktABDByxtBDBkedd),( ,),(,21212121)1 (42轮换。mjixxcyybmjisrbbcccbbcbccbccbbAEtkmjimjisrsrsrsrsrsrsrsrrs k元素的物理意义元素的物理意义 kpq：第第q个结点位移分量为单位位移（其它结个结点位移分量为单位位移（其它结点位移点位移=0），所引起的第），所引起的第p个结点力分量。个结点力分量。 如如: k26。 k的性质：的性质： (1) 对称性：对称性： kpq= kqp (2)

14、 奇异性；奇异性； (3) 每行（列）元素之和为零。每行（列）元素之和为零。 (4) k取决于单元的形状、方位和弹性常数，与取决于单元的形状、方位和弹性常数，与所在位置（即平移或所在位置（即平移或np p 转动）无关。转动）无关。2) 单刚性质单刚性质ymxk26ij13) 单刚列式推导的另一方法单刚列式推导的另一方法yui , (Ui)ixvi , (Vi)ujjvjummvm由单元平衡条件导出由单元平衡条件导出物理方程物理方程几何方程几何方程位移模式位移模式虚功方程虚功方程f =N e =B e =S e ，S=DBFe=k e，k=BT DBtA结点位移结点位移 位移位移 应变应变 应力

15、应力 结点力结点力 e f Fe虚功原理虚功原理eeAeTTeeTeATeTeyxtBDBFyxtFdd)()(dd)(eATeyxtBDBkdd令：eeekF 应力与结点力关系式：应力与结点力关系式：yui , (Ui)ixvi , (Vi)ujjvjummvm假设单元发生虚位移假设单元发生虚位移 结点力虚功结点力虚功=单元虚变形能单元虚变形能列式列式2.3.6 单元等效结点荷载单元等效结点荷载将非结点荷载等效移置到结点上：将非结点荷载等效移置到结点上：dddPNRstqNRyxtpNRRRRTePSTeqATepeqepeeeyXiixYiXjjYjXmmYmPxPyP可用虚功原理导得可用

16、虚功原理导得2.3.7 整体刚度矩阵整体刚度矩阵建立整体刚度矩阵的方法建立整体刚度矩阵的方法 1) 由最小势能原理建立：由最小势能原理建立：K= (Ge)Tk e Ge 2) 由结点平衡条件建立：由结点平衡条件建立：eieiRFeijijkKK的集成方法：的集成方法：对号入座对号入座 K= (Ge)Tke Ge 。13245(1)(3)(2)Rx3Ry3U3V33RK 总刚总刚K的集成的集成 例例 1: E, t, =常数常数1) 离散离散 4个单元，个单元，6个结点个结点 单元编号：单元编号：(1)(4)； 结点编号：结点编号：1 62) 建立局部建立局部-整体编码关系整体编码关系(1) (

17、2) (4)(3)i j mijmyx1kN/m1m1m1m1m123456(1)(2)(3)(4)Local No.Global No.Element(1)(2)(3)(4)ijm3125242536353) 计算各单刚计算各单刚k)4()2() 1 () 1 () 1 () 1 () 1 () 1 () 1 () 1 () 1 () 1 ( ,kkkkkkkkkkkkmmmjmijmjjjiimijii(1) (2) (4)(3)i j mijmyx1kN/m1m1m1m1m123456(1)(2)(3)(4)4) 换码，对号入座，形成总刚换码，对号入座，形成总刚(1) (2) (4)(3

18、)i j mijmyx1kN/m1m1m1m1m123456(1)(2)(3)(4) 1()3() 1() 1()3() 1()3()2() 1() 1() 1() 1() 1(iimiimijimmiiijjmmmjjijmjjkkkkkkkkkkkkkK)1()1()1()1()1()1()1()1()1()1(mmmjmijmjjjiimijiikkkkkkkkkk对一总刚元素有贡献的单刚：对一总刚元素有贡献的单刚：(1) (2) (4)(3)i j mijmyx1kN/m1m1m1m1m123456(1)(2)(3)(4)4()3()2(55)4()3(53)(,mmjjiimjjme

19、eijpqkkkKkkKkK例例 2: E, t=常数，常数， =0；l12=l13=l35=1m，形成，形成 K。(1)(3)i j m(2)ijm13245(1)(3)(2)232121121232112121002110Symm.1021)1 (22)3()2()1(Etkkkk75. 025. 025. 005 . 025. 075. 025. 05 . 0025. 025. 00025. 05 . 000.5 . 0025. 0SymmEtk(1)(3)i j m(2)ijm13245(1)(3)(2)ijmi j m5 . 125. 025. 05 . 1025. 025. 002

20、5. 025. 005 . 0)3()2()1(33)2()1(32)1(31mmiijjijjijmkkkKkkKkK( =0)5 . 00025. 025. 025. 000005 . 00005 . 0000025. 125. 00 . 125. 025. 000025. 00 . 125. 05 . 025. 05 . 0005 . 125. 0025. 025. 025. 025. 05 . 125. 0005 . 075. 005 . 0075. 025. 025. 0Symm.75. 025. 075. 0EtK(1)(3)i j m(2)ijm13245(1)(3)(2) 整体

21、刚度元素整体刚度元素Kpq的物理意义：的物理意义： 结构第结构第q个结点位移为单位位移（其它结点位个结点位移为单位位移（其它结点位移移=0）时，所引起的第）时，所引起的第p个结点力。个结点力。 K的性质：的性质： (1) 对称性对称性 (Kpq= Kqp), 主对角元素必为正主对角元素必为正； (2) 稀疏性，且一般为带状分布；稀疏性，且一般为带状分布； 平面问题最大半带宽平面问题最大半带宽= 2 (单元结点号之差最大值单元结点号之差最大值+1) (3) 引入约束条件后为正定矩阵。引入约束条件后为正定矩阵。 利用利用对称性对称性、带状稀疏性带状稀疏性， K可用半带宽存储。可用半带宽存储。 12

22、3456789101112对称69686658575655474544363533252423221412110000KKKKKKKKKKKKKKKKKKKKK02.3.8 位移边界条件的引入位移边界条件的引入1.直接引入法（矩阵缩小法）直接引入法（矩阵缩小法）K =R Kaa a=Ra Kab b bababbbaabaaRRKKKK未知结点位移向量未知结点位移向量已知已知结点位移向量结点位移向量解出未知量解出未知量破坏原矩阵，计算机分析中很少采用。破坏原矩阵，计算机分析中很少采用。5 . 00025. 025. 025. 000005 . 00005 . 0000025. 125. 00

23、. 125. 025. 000025. 00 . 125. 05 . 025. 05 . 0005 . 125. 0025. 025. 025. 025. 05 . 125. 0005 . 075. 005 . 0075. 025. 025. 0Symm.75. 025. 075. 0EtK(1)(3)i j m(2)ijm13245(1)(3)(2)5 . 00025. 025. 025. 05 . 00005 . 025. 125. 00 . 125. 00 . 125. 05 . 0Symm5 . 125. 05 . 1EtK(1)(3)i j m(2)ijm13245(1)(3)(2)

24、划行（划行（1-4行）划列（行）划列（1-4列）：列）：2. 对角元素改对角元素改1法（零位移边界）法（零位移边界） i=0 Kii=1 , Kij=0, Kji=0（j i）, Ri=0 NNiNNNNNNRRRKKKKKKKKK0001000021212122221112113. 乘大数法乘大数法 i= i Kii= (大数大数)，Ri= Kii iNiiiNiNNNiNNiNiiiiNiNiRKRRKKKKKKKKKKKKKKKK21212121222221111211第第i个方程：个方程： 等价于：等价于：iiNiNiiiiiKKKKK2211ii2.3.9 位移模式与解答的收敛准则位

25、移模式与解答的收敛准则完备性条件完备性条件协调性条件协调性条件完备性条件完备性条件 = 解答收敛的必要条件解答收敛的必要条件完备完备 + 协调协调 = 解答收敛的充分条件解答收敛的充分条件(1) 位移模式必须能反映单元的刚体位移位移模式必须能反映单元的刚体位移(2) 位移模式必须能反映单元的常量应变位移模式必须能反映单元的常量应变(3) 位移模式应尽可能反映位移的连续性位移模式应尽可能反映位移的连续性1. 收敛准则收敛准则 三结点三角形单元的完备性和协调性三结点三角形单元的完备性和协调性(1) 反映刚体位移：反映刚体位移：平面问题的刚体位移表达式：平面问题的刚体位移表达式：u=u0 w wy，

26、v=v0+w wx(2) 反映常量应变：反映常量应变： x= 2， y= 6， xy= 2+ 3 满足完备性条件满足完备性条件xxyvyyxu2222353564353521(3) 位移连续性：位移连续性： 单元内：单值连续；单元内：单值连续； 相邻单元间：相邻单元间： uij(1)=uij(2)？vij(1)=vij(2) ？ij边的方程边的方程：y=ax+b，则则 uij= 1+ 2 x+ 3(ax+b)= cx+d uij(1)、uij(2)均为坐标的线性函数，故可由均为坐标的线性函数，故可由i、j两两点的结点位移唯一确定。点的结点位移唯一确定。 满足连续性条件满足连续性条件yx(2)i

27、jmp(1)2. 多项式位移模式的选择多项式位移模式的选择 (1) 一般规则：与一般规则：与局部坐标选取无关局部坐标选取无关 几何各向同性几何各向同性 x、y各阶项对称：各阶项对称： xmyn xnym (2) 选择方法选择方法从从Pascal三角形三角形中选项：中选项： 1 x y x2 xy y2 x3 x2y x y2 y3 x4 x3y x2y2 xy3 y4 3. 位移元解答位移元解答位移元：位移元：最小势能原理建立的有限元最小势能原理建立的有限元 基本未知量：结点位移基本未知量：结点位移 位移元近似解位移元近似解=精确解下限精确解下限 2.4.1 高精度三角形单元高精度三角形单元1

28、. 六结点三角形单元（二次单元六结点三角形单元（二次单元T6） 位移模式位移模式（完整二次式完整二次式）：u= 1+ 2x+ 3 y + 4 x2 + 5 xy+ 6 y2v= 1+ 2 x + 3 y + 4 x2 + 5 xy+ 6 y2 满足：满足： 协调单元。协调单元。2.4 高精度的三角形单元、矩形单元高精度的三角形单元、矩形单元yxijm312完备性条件完备性条件连续性条件连续性条件2.4.1 高精度三角形单元高精度三角形单元 (续续)2. 十结点三角形单元（三次单元十结点三角形单元（三次单元T10） 位移模式位移模式（完整完整三三次式次式）：u= 1+ 2x+ 3y+ 4x2+

29、5xy+ 6y2+ 7x3 8x2y+ 9xy2+ 10y3v= 1+ 2x+ 3y+ 4x2+ 5xy+ 6y2+ 7x3 8x2y+ 9xy2+ 10y3 满足：满足： 协调单元。协调单元。完备性条件完备性条件连续性条件连续性条件yxijm3. 面积坐标表达的形函数面积坐标表达的形函数 定义：定义： 三角形单元内任一点三角形单元内任一点P的无量纲面积坐标：的无量纲面积坐标： Li= Ai / A （i, j, m） 3结点三角形单元结点三角形单元 (T3单元单元)： Ni=Li （i, j, m） 6结点三角形单元结点三角形单元 (T6单元单元)： Ni=(2Li 1) Li （i, j,

30、 m） N1=4Lj Lm （1, 2, 3；i, j, m）i(1,0,0)Pj(0,1,0)m(0,0,1)AiAjAm 面积坐标的性质面积坐标的性质 (1) 三角点的面积坐标：三角点的面积坐标： i (1,0,0)、j (0,1,0)、m (0,0,1) (2) 三条边的方程：三条边的方程： jm边：边：Li=0， mi边：边：Lj=0， i j边：边：Lm=0。 (3) 三个面积坐标不独立：三个面积坐标不独立： Li+ Lj + Lm =1i(1,0,0)Pj(0,1,0)m(0,0,1)AiAjAm 面积坐标与直角坐标之间的关系面积坐标与直角坐标之间的关系: i(1,0,0)Pj(0

31、,1,0)m(0,0,1)AiAjAmyxyxcbacbacbaALLLmmmjjjiiimji121mjimjimjiLLLyyyxxxyx1111或(1) 位移位移模式模式 双线性位移模式：双线性位移模式：u= 1+ 2 x+ 3 y+ 4 xyv= 1+ 2 x + 3 y+ 4 xy 满足：满足： 协调单元。协调单元。 结点位移表示的结点位移表示的位移模式：位移模式：2.4.2 四结点矩形单元（四结点矩形单元（R4单元）单元）完备性条件完备性条件连续性条件连续性条件1yx234o4141,iiiiiivNvuNu 形函数的构造形函数的构造： 设：设：N1=A A1(x-x2 ) (y-

32、y3 ) 令令：(N1)1 =1 可得：可得：)1)(1 (414141003211byyaxxbyyaxxNabA4141,iiiiiivNvuNu1yx234 o2a2b)1)(1 (41)1)(1 (41)1)(1 (41004003002byyaxxNbyyaxxNbyyaxxN，(2) 局部坐标下的形函数局部坐标下的形函数 建立局部坐标系建立局部坐标系 (自然坐标自然坐标) o ， 坐标变换式：坐标变换式：byyaxx00，1yx234 o2a2b1234 o = -1 =1 = -1 =1 四角点坐标：四角点坐标：1(-1, -1)， 2(1, -1) 3(1, 1)， 4(-1,

33、 1) 位移模式：位移模式：4141,iiiiiivNvuNu(2) 局部坐标下的形函数局部坐标下的形函数)1)(1 (41),1)(1 (41)1)(1 (41),1)(1 (414321NNNN)4 , 3 , 2 , 1()1)(1 (41iNiii1234 o = -1 =1 = -1 =1 位移模式：位移模式： 形函数：形函数： 形函数统一式：形函数统一式：4141,iiiiiivNvuNu(3) 单元应变矩阵和应力矩阵单元应变矩阵和应力矩阵 单元应变：单元应变： =B e， B=B1 B2 B3 B41234 o = -1 =1 = -1 =1 单元应力：单元应力： =S e=DB

34、 e S= S1 S2 S3 S4，Si=D Bi )4, 3, 2, 1(00100iNbNaNaNbabxNyNyNxNBiiiiiiiii2.5 C0连续型单元形函数的构造连续型单元形函数的构造矩形单元类型矩形单元类型Lagrange矩形单元矩形单元 Serendipity矩形单元矩形单元Lagrange单元单元各网格交点均布置结点各网格交点均布置结点Serendipity单元单元仅单元边界布置结点仅单元边界布置结点2.5 C0连续型单元形函数的构造连续型单元形函数的构造（续）续）2.5.1 Lagrange矩形单元矩形单元1. 一维一维Lagrange插值多项式插值多项式 ), 2 ,

35、 1()(, 1)1(nixxxxxlnijjjijni表示多项式次数表示多项式次数xx1x2x3xnxil1(x)曲线：过曲线：过n个结点个结点例：例：n=2，则，则若令若令x1=0，x2= l，则，则 l1(1)(x)=(l-x)/l， l2(1)(x)=x/l 。121)1(2212)1(1)(,)(xxxxxlxxxxxl可作形函数可作形函数, 即即Ni (x)= li(n-1)(x)2. Lagrange矩形单元的形函数矩形单元的形函数水平向：水平向：r +1个结点，个结点，竖直向：竖直向：p +1个结点。个结点。第第I列列J行结点行结点i的形函数：的形函数：Ni= NIJ =lI(

36、r)( ) lJ(p)( ) (0, 0)(r, p)1lI ( ) 1lJ ( )i (I, J)(,)()(,)(0, 0)(0, 0)(byylaxxlpIjjjJjpJrIjjjIjrINi 能保证边界位移的唯一性和协调性。能保证边界位移的唯一性和协调性。Lagrange矩形单元的优缺点：矩形单元的优缺点：优点：形函数容易构造；为完备协调单元。优点：形函数容易构造；为完备协调单元。缺点：单元自由度较多；内部结点不能有效提高精度。缺点：单元自由度较多；内部结点不能有效提高精度。3 3结点结点Lagrange单元形函数包含的项数单元形函数包含的项数：1x yx2 xy y2 x3 x2y

37、xy2 y3x3y x2y2 xy3x3y2 x2y3 x3y32.5.2 Serendipity矩形单元矩形单元1. 单元特点单元特点只在边界上布置结点，只在边界上布置结点，不同边界上结点数目可不同。不同边界上结点数目可不同。2. 形函数构造方法形函数构造方法R4单元：单元：R5单元：单元：R4单元单元 o1243)4 , 3 , 2 , 1()1)(1 (41iNiii o15243R5单元单元)1)(1 (2125N进行修正，对4321NNNN2.5.2 Serendipity矩形单元矩形单元2. 形函数构造方法（续）形函数构造方法（续）R5单元：单元：进行修正，对4321NNNN o1

38、5243R5单元单元)1)(1 (2125N51121NNN1N51/21011N443352221NNNNNNN，；对对R8单元：单元：R8单元单元 o15247386)1)(1 (21),1)(1 (21)1)(1 (21),1)(1 (2128272625NNNN87447633652285112121,21212121,2121NNNNNNNNNNNNNNNN)4 , 3 , 2 , 1(),(41444iNNNNjjjjiii或2.6 平面等参平面等参数单元数单元规则形状单元规则形状单元2.6.1 四结点四边形等参单元四结点四边形等参单元 在单元内建立局部坐标在单元内建立局部坐标 使

39、局部坐标边界取值为常量使局部坐标边界取值为常量采用局部坐标的双线性位移模式采用局部坐标的双线性位移模式复杂形状单元复杂形状单元坐标变换坐标变换等参数变换：等参数变换：变换函数变换函数=单元形函数单元形函数等参数单元：通过等参变换得到的单元等参数单元：通过等参变换得到的单元1yx234 1234 o = -1 =1 = -1 =1x, y的双线性位移模式的双线性位移模式不满足协调性不满足协调性2. 局部局部-整体坐标映射：整体坐标映射： 映射的有效性验证：映射的有效性验证：1) 四个结点四个结点给出整体坐标；给出整体坐标； 2) 四条边四条边给出整体直线方程。给出整体直线方程。 12边：边：1y

40、x234 1234 o = -1 =1 = -1 =1基本单元（母单元）基本单元（母单元） 实际单元（子单元）实际单元（子单元） 4141,iiiiiiyNyxNx212143212121,21210),1 (21),1 (211yyyxxxNNNN，2.6 平面等参平面等参数单元（续）数单元（续）利用等参变换利用等参变换，可构造更高次四边形曲边单元：，可构造更高次四边形曲边单元： 8结点四边形曲边等参单元结点四边形曲边等参单元； 12结点四边形曲边等参单元结点四边形曲边等参单元； 20结点四边形曲边等参单元结点四边形曲边等参单元8结点四边形等参曲边单元结点四边形等参曲边单元 o1524738

41、6yx152473862.6.2 局部与整体坐标的微分和积分变换式局部与整体坐标的微分和积分变换式根据根据复合函数求导法则：复合函数求导法则：yyxxyyxxyxyxyxyxJ矩阵形式矩阵形式 miiimiiimiiimiiiyNxNyNxNyxyx1111J雅可比雅可比(Jacobi)矩阵矩阵m平面等参单元平面等参单元的结点数目的结点数目2.6.2 局部与整体坐标的微分和积分变换式局部与整体坐标的微分和积分变换式反变换：反变换：面积微分的变换：面积微分的变换：dA=dxdy=|J |d d xxyyyx|11JJ*|1222112111JJJJJJJij*J元素元素Jij 的代数余子式的代数

42、余子式|J|雅可比雅可比(Jacobi)行列式行列式2.6.3 单元刚度矩阵、单元等效结点荷载向量单元刚度矩阵、单元等效结点荷载向量单元应变向量：单元应变向量： =B e单元应变矩阵：单元应变矩阵： B=B1 B2 Bm 单元应力向量：单元应力向量： =S e =DB e单元刚度矩阵：单元刚度矩阵：), 2, 1(00mixNyNyNxNBiiiii 1111dd | ddJtBDByxtBDBkTATe单元等效结点荷载列阵单元等效结点荷载列阵：单元边界线微分：单元边界线微分： eeSTeqTATepeqepestqNRJtpNyxtpNRRRRddd|dd1111面力：体力：)(dd)()(

43、d)(dd)()(d2/1222/122常数常数syxssyxs积分式常用积分式常用Newton-Cotes或或Gauss数值积分计算数值积分计算积分在母单元中进行，可用标准程序实现。积分在母单元中进行，可用标准程序实现。2.6.4 等参变换的应用条件等参变换的应用条件母单元母单元 实际单元实际单元：一一对应关系一一对应关系 局部局部整体坐标变换式：整体坐标变换式：J非奇异或非奇异或|J | 0为确保坐标变换一一对应，应避免：为确保坐标变换一一对应，应避免： 1) 任意两个结点退化为一个，而使任意两个结点退化为一个，而使|d |=0或或|d |=0; 2) 单元过于歪斜，而导致单元过于歪斜，而

44、导致d 与与d 共线。共线。 123,4 = 1 =1 = 1 =1o o 123 = 1 =1 = 1 =1o o42.6.5 例题分析例题分析悬臂平面结构，长悬臂平面结构，长2m，高，高1m。E=常数，常数， =0.3。方法方法1：利用三结点三角形单元进行分析利用三结点三角形单元进行分析12346(1)(2) (3)5(4)yxP jmi(1)(3) jmi(2)(4)1) 集成整体刚度矩阵集成整体刚度矩阵;2) 引进位移边界条件。引进位移边界条件。2.6.5 例题分析（续）例题分析（续）分析步骤：分析步骤：划分单元，结点编号；划分单元，结点编号；单元编号，各单元结点号；单元编号，各单元结

45、点号；计算单元等效结点荷载；计算单元等效结点荷载；计算单刚参数计算单刚参数bi , ci等，形成单刚；等，形成单刚；集成整体刚度矩阵集成整体刚度矩阵;引进位移边界条件；引进位移边界条件；解线性代数方程；解线性代数方程；计算单元应力、应变；计算单元应力、应变；计算结果整理、插值、输出等。计算结果整理、插值、输出等。12346(1)(2) (3)5(4)yxP三角形单元划分应注意的问题三角形单元划分应注意的问题:1. 单元大小及疏密单元大小及疏密根据精度、计算能力综合根据精度、计算能力综合. . 主要部位、应力（位移）变化大的部位划细主要部位、应力（位移）变化大的部位划细 应力误差与单元尺寸成正比

46、应力误差与单元尺寸成正比 位移误差与位移误差与(单元尺寸单元尺寸)2成正比成正比2. 三边尽量接近三边尽量接近应力、位移误差反比最小内角之正弦应力、位移误差反比最小内角之正弦. .3. 尺寸或材料突变处划作单元边界，附近单元划小；尺寸或材料突变处划作单元边界，附近单元划小；4. 荷载突变或集中荷载处布置结点，附近单元划小。荷载突变或集中荷载处布置结点，附近单元划小。可应用于其他类型单元，如四边形单元可应用于其他类型单元，如四边形单元.2.6.5 例题分析（续）例题分析（续）悬臂平面结构，长悬臂平面结构，长2m，高，高1m。E, =常数常数。方法方法2：利用四结点矩形单元进行分析利用四结点矩形单

47、元进行分析12346(1)(2)5yxP1) 集成整体刚度矩阵集成整体刚度矩阵;2) 引进位移边界条件。引进位移边界条件。悬臂平面结构，长悬臂平面结构，长2m，高，高1m。 E, =常数常数。方法方法3：利用：利用四结点四边形等参单元四结点四边形等参单元计算步骤：计算步骤： 离散化离散化;求求 ，J，|J|，J-1；求求B在各积分点的数值在各积分点的数值Big；利用高斯积分计算并形成利用高斯积分计算并形成k；集成集成K、P；引进位移边界条件。引进位移边界条件。12346(1)(2)5yxP0.75m 1.25m1.25m 0.75m)4 , 3 , 2 , 1(,iNNii悬臂平面结构，长悬臂

48、平面结构，长2m，高，高1m。 E, =常数常数。四结点四边形等参单元四结点四边形等参单元单元分析书面作业单元分析书面作业（单独完成）（单独完成） :求求(1)、(2)单元内单元内( =0, =0)点对应的点对应的 (x, y)坐标坐标;已知各结点已知各结点ui , vi (i=3, 4, 5, 6)，求，求(x=0.625, y=0)处的位移。处的位移。12346(1)(2)5yxP0.75m 1.25m1.25m 0.75m悬臂平面结构，长悬臂平面结构，长2m，高，高1m。 E, =常数常数。方法方法4：利用八结点等参单元进行分析利用八结点等参单元进行分析136813(1)(2)11yxP

49、245712910上机作业上机作业:（单独完成，与编程作业任选一）（单独完成，与编程作业任选一） 重力坝问题重力坝问题三层框架三层框架-楼板体系楼板体系 (1) 对计算结果进行分析、比较；对计算结果进行分析、比较； (2) 分别加密网格，分析、比较结果的改进情况。分别加密网格，分析、比较结果的改进情况。2.7 有限元程序实现有限元程序实现有限元程序有限元程序 前处理程序前处理程序生成网格及数据文件生成网格及数据文件主体分析程序主体分析程序核心计算分析核心计算分析后处理程序后处理程序结果处理，生成结果文件结果处理，生成结果文件2.7.1 程序设计一般步骤程序设计一般步骤算法描述和列式推导算法描述

50、和列式推导框图设计框图设计代码编写代码编写上机调试、考核上机调试、考核编写应用说明编写应用说明修改、补充、完善修改、补充、完善 2.7 有限元程序实现有限元程序实现（续）（续）程序设计一般要求程序设计一般要求 功能较齐全功能较齐全通用性通用性、可移植性较强可移植性较强较好可扩充性较好可扩充性良好可读性良好可读性足够可靠性足够可靠性良好的适应性良好的适应性2.7.2 输入数据及分类输入数据及分类1. 控制数据：控制数据：结点总数、单元总数、约束总数、结点总数、单元总数、约束总数、 荷载总数、问题类型数荷载总数、问题类型数等等;2. 几何数据：几何数据：结点坐标、单元信息结点坐标、单元信息 (各单

51、元的结点各单元的结点 编号编号)、约束条件、单元类型数约束条件、单元类型数(E, , 厚度不同为一类厚度不同为一类);3. 物性数据：物性数据：E、 、t（厚度）（厚度）4. 荷载数据：荷载数据：荷载类型荷载类型（集中、分布）、（集中、分布）、位置、位置、 方向、大小方向、大小等等.2.7.3 程序设计方法程序设计方法1. 模块化方法模块化方法1个主程序个主程序+若干子程序；若干子程序；2. 动态数组存储技术：动态数组存储技术：数组共享、覆盖技术数组共享、覆盖技术 编程作业编程作业:（1-3人分组完成，与上机作业任选一）人分组完成，与上机作业任选一） 编制平面问题的四结点或八结点等参单元分析程

52、序，编制平面问题的四结点或八结点等参单元分析程序，提交源程序和提交源程序和2个算例。个算例。2.8 平面平面杆件杆件结构的有限元结构的有限元2.8.1 等截面直梁单元（忽略剪切变形）等截面直梁单元（忽略剪切变形） 1. 基本方程基本方程几何关系：几何关系：内力内力-位移关系：位移关系：22dd,ddxvxux22dd00dd,xxLfLx矩阵形式矩阵形式 22dd,ddxvEIEIMxuEAEANxEIEADDMN00,矩阵形式矩阵形式 xyp(x)q(x)2.8.1 等截面直梁单元（续等截面直梁单元（续1） 1. 基本方程基本方程平衡关系：平衡关系：边界条件：边界条件：)(dddd,dddd

53、),(dd223322xqxMxQxvEIxMQxpxuEA00:;,)(00:;,)(00:;,)(MQMMQQMvMMvvvdxdvvvssssss，自由如，简支如，固支如2. 离散化离散化将一平面杆件结构离散为将一平面杆件结构离散为ne个单元，个单元，n个结点个结点2.8.1 等截面直梁单元（续等截面直梁单元（续2） 2. 离散化离散化离散为离散为ne个单元，个单元，n个结点：个结点： 结点位移向量：结点位移向量： =u1 v1 1 u2 v2 2 un vn nT 结点力向量：结点力向量： F=U1 V1 M1 U2 V2 M2 Un Vn MnT 3. 位移模式位移模式局部坐标下局部

54、坐标下 单元结点位移向量：单元结点位移向量： 单元结点力向量：单元结点力向量：TjjjiiieTjjjiiieMQNMQNFvuvuxyxyviuivjuj i j 2.8.1 等截面直梁单元（续等截面直梁单元（续3）3. 位移模式位移模式结点荷载下：结点荷载下：梁轴向位移沿梁轴呈线性分布，梁轴向位移沿梁轴呈线性分布， 横向位移呈三次曲线分布。横向位移呈三次曲线分布。位移模式：位移模式： u= 1+ 2x v= 3+ 4 x+ 5 x2+ 6 x3 或或 u=N1ui+N4 uj v=N2vi+ N3 i+ N5vj+ N6 j23263322542323332221,23,2,231,1lx

55、lxNlxlxNlxNlxlxxNlxlxNlxN2.8.1 等截面直梁单元（续等截面直梁单元（续4）3. 位移模式位移模式 矩阵表达式：矩阵表达式： u、v独立插值，独立插值， 不独立插值不独立插值要求要求C1连续。连续。653241000000NNNNNNNNvufe2.8.1 等截面直梁单元（续等截面直梁单元（续5）4. 单元刚度矩阵单元刚度矩阵单元应变矩阵、应力矩阵单元应变矩阵、应力矩阵f =Ne代入代入几何方程几何方程=Be代入代入物理方程物理方程=DBeexBfL22dd00dd,xxLNLBEIEAD002.8.1 等截面直梁单元（续等截面直梁单元（续6）4. 单元刚度矩阵单元刚

56、度矩阵局部坐标下局部坐标下:lTxBDBkdlEIlEIlEIlEIlEIlEIlEIlEIlEAlEAlEIlEIlEIlEIlEIlEIlEIlEIlEAlEAk4602606120612000002604606120612000002223232223232.8.1 等截面直梁单元（续等截面直梁单元（续7）5. 单元刚度矩阵的坐标变换单元刚度矩阵的坐标变换结点位移向量：结点位移向量： e=TT e结点力向量：结点力向量： Fe =TTFe刚度矩阵：刚度矩阵： k =TTkT1000cossin0sincos,0066tttT坐标变换矩阵：坐标变换矩阵： 反变换：反变换：eeeeFTFT

57、, xyxyviuivjuj i j 2.8.2 考虑剪切变形的梁单元考虑剪切变形的梁单元基本假定：平截面，不垂直于变形后轴线。基本假定：平截面，不垂直于变形后轴线。挠度计算：挠度计算：v=vb+vs； 杆轴线斜率：杆轴线斜率：dv/dx= + 。lxNlxN87,1dxdx截面转角截面转角 剪应变剪应变 xyq(x)dx1. 在经典梁单元基础上引入剪切变形影响在经典梁单元基础上引入剪切变形影响vb用三次式插值；用三次式插值；vs用线性插值：用线性插值： vs= N7vsi + N8vsj，1. 在经典梁单元基础上引入剪切变形影响在经典梁单元基础上引入剪切变形影响单元刚度矩阵：单元刚度矩阵：d

58、x)1 ()4()1 (60)1 ()2()1 (60)1 (120)1 (6)1 (12000)1 ()4)1 (60)1 (1202232323lEIlEIlEIlEIlEIlEIlEIlEAlEAlEIlEIlEIlEAk（对称212GAlEIk2. Timoshenko梁单元梁单元 特点特点对对u、v、 均均独立插值：独立插值： u=N1u1+N2 u2位移模式：位移模式： v= N1v1+ N2v2 =N1 1+ N2 2 1c2l/2l/22,)(2),1 (21),1 (212121xxxlxxNNcc320620102100000320100100100000000010000

59、222lllllllklGAIAIAIAlEIkkksb对称对称2. Timoshenko梁单元梁单元优势：位移模式简单；优势：位移模式简单；缺点：缺点：l /h 时时, 只能得到零解只能得到零解剪切自锁剪切自锁(Shear locking) 原因：原因：v、 同阶插值，同阶插值，dv/dx与与 不同阶不同阶 =dv/dx- =0不能满足，除非不能满足，除非 =常数。常数。解决方法：解决方法：减缩积分减缩积分数值积分时采用比精确积分要求少的积分数值积分时采用比精确积分要求少的积分点数，例如对两结点梁单元采用一点积分。点数，例如对两结点梁单元采用一点积分。假设剪切应变假设剪切应变对剪应变对剪应变

60、g另行假定插值形式；另行假定插值形式；替代插值函数替代插值函数对对 采用低一阶的插值函数，如两结点采用低一阶的插值函数，如两结点梁单元其插值函数为常数梁单元其插值函数为常数1/2，即，即:2211NN2. Timoshenko梁单元梁单元 两结点两结点Timoshenko梁单元：梁单元： 1) 包含横向刚体位移包含横向刚体位移 (v=c)、刚体转动、刚体转动 ( =dv/dx =c); 2) 包含常剪切应变包含常剪切应变 ( =0，dv/dx = =c); 3) 不包含常弯曲应变的位移状态不包含常弯曲应变的位移状态 ( =cx，v=0.5cx2) 不能分析纯弯问题（伴随剪切应变）。不能分析纯弯