英语原文共 75 页，剩余内容已隐藏，支付完成后下载完整资料
第五章 三角形常应变单元与二维问题求解
5.1 引言
本章在推导二维有限元求解列式时，按照求解一维问题的步骤进行，位移、面力和分布体力均为位置坐标（x,y）的函数，位移向量u表示如下：
（5.1）
其中u、v分别为u在x和y方向上的分量。应力和应变表示如下：
（5.2）
（5.3）
对于图5.1中描述的二维问题，体力、面力和微元体积表示如下：
（5.4）
其中t是z方向上的厚度。体力f为单位体积上所受的力，而面力T为单位面积上所受的力。应变-位移关系表示如下：
（5.5）
应力和应变之间的关系为（参见式（1.18）和式（1.19））
（5.6）
对求解域进行离散的基本思想是将域内任意一点的位移用离散点的位移值来表示。本章首先介绍三角形单元，然后分别应用能量法和Galerkin方法推导出相应的刚度和载荷。
5.2 有限元模型
将二维求解域离散为若干个直边三角形，图5.2所示的是一种典型的三角形单元划分方式，各三角形的交点称为节点，由三个节点和三条边构成的三角形称为单元（element）。除了边缘部分的微小区域外，单元将覆盖整个区域，由于选用直边三角形，造成曲线边界不能完全逼近，存在一些小的未覆盖区域，如果减小单元尺寸或选择具有曲线边界的单元，这部分未被离散的区域将减小；有限元方法的基本思想是近似地求解连续问题，将连续的求解域离散为单元时所产生的的离散误差是造成近似解误差的一个重要原因。对于图5.2中所示的三角形单元划分，节点序号标在个三角形的角上，单元序号用带圆圈的数字表示。
在本章所讨论的二维问题中，每个节点由x和y两个方向上的位移分量，即每个节点有两个自由度。从分析桁架时采用的编码方案可以看出：节点j在x方向上的位移分量表示为，在y方向上的位移分量表示为，我们将整个位移列阵表示为
（5.7）
其中N是自由度的个数。
在进行计算时，通常用节点坐标及单元的节点连接信息来表示单元的离散信息。节点坐标信息存储在一个二维列阵里，其中行数为最大节点数，每一列对应每个节点的两个坐标。如图5.3中所示，单独分析一个典型单元，可以获得其单元的节点信息。单元三个节点的局部编号如图5.2所示；单元的节点信息通常由一个矩阵表示，矩阵中包含了单元尺寸、数目和各单元三个节点的编号；表5.1是一种典型的单元节点信息表示方法，为了避免三角形面积出现负值，大多数有限元程序采用常规的绕单元逆时针顺序排列节点的方法；然而在本章所给的程序中，对节点的次序并不做特殊的要求。
表5.1建立了局部和整体的节点序号及自由度之间的对应关系。图5.3中的局部节点j在x方向和y方向上的位移分量分别表示为和，因此，单元位移列阵可表示为
（5.8）
注意由表5.1中的单元节点信息，我们能够从整体Q列阵中提取出q列阵，在有限元程序中常有这种运算。同样，可由表5.1看出：对应节点的坐标（）,()和（）是具有整体坐标性质的。节点坐标和自由度的局部表示法为单元表征提供了一种简便而清晰的方法。
5.3 常应变三角形单元（CST）
在进行有限元计算时，单元内部各位置的位移需要用单元节点位移来表示。前面已经讨论过，在单元内进行插值运算时，需要引出形状函数的概念；对于常应变三角形单元，其形状函数是线性的；如图5.4所示，三个形状函数分别对应于节点1、2和3。形状函数在节点1处的值为1，而在节点2和3出呈线性地减小为0；图5.4（a）中用阴影表示的平面就是形状函数的值，分别表示了其形状函数在节点2和3出的值为1而在对边降为0的类似平面；这些形状函数的任何线性组合都代表一种平面。特别地，代表了一个在节点1、2和3处的值均为1的平面，它与三角形123平行。因此，对于任何的，都存在下列关系：
（5.9）
因此并非是线性独立的，它们中仅有两个是独立的，形状函数通常用变量zeta;和eta;表示如下：
（5.10）
其中zeta;、eta;是自然坐标（natural coordinates），参见图5.4.请注意它与一维单元（见第三章）的相似处：在一维问题中，将x坐标映射到zeta;坐标，形状函数定义为zeta;的函数；在二维问题中，将x、y坐标分别映射到zeta;、eta;坐标上，因此形状函数定义为zeta;和eta;的函数。
实际上形状函数可以用面积坐标（area coordinates）表示，如图5.5所示，三角形中的一点（x,y）将三角形划分为三个区域，形状函数可以精确地表示为
（5.11）
其中A是单元面积。显然对于三角形内所有点，关系都成立。
等参表示法
现在，单元内任意一点的位移都可以用形状函数和未知位移场的节点值来表示。由此可得
（5.12a）
或者，由式（5.10）可得
（5.12b）
定义形状函数矩阵后，是（5.12a）可以表示为如下的矩阵形式
（5.13）
和
（5.14）
对于三角形单元，采用相同的形状函数，可以将x,y坐标表示为节点坐标的形式。这就是等参表示法（isoparametric representation），应用等参表示法可以简化推导过程，这使得在拓展到复杂单元时仍可以采用相同的格式。由此可得知坐标的插值关系为
（5.15a）
或者
（5.15b）
利用如下表示法：，式（5.15b）可表示为
（5.15c）
这个公式将x，y坐标和zeta;、eta;坐标联系起来。式（5.12）表明u和v是zeta;和eta;的函数。
例题 5.1
对于例题5.1图中的三角形单元，分别求单元内一点P的形状函数的值。
解答：由等参表示法（式5.15），可得
以上梁式可化简为
求解得xi;=0.3，eta;=0.2，由此可得。在计算单元应变时，要分别计算u和v对于整体坐标x和y的偏导数。从式（5.12）和式（5.15）可以看出，u、v和x、y都是xi;、eta;的函数，即有u=u（x（xi;，eta;），y（xi;，eta;））和v=v（x（xi;，eta;），y（xi;，eta;）），应用求偏导数的链式法则，u的偏导数可以表示如下：
上式可写成矩阵形式，即
（5.16）
将上式中的（2*2）方阵记为Jacobian变换矩阵
（5.17）
关于Jacobian的一些其他特性将在附录中列出，计算出x和y的导数，可得
（5.18）
同样，由式（5.16）可得
（5.19）
其中是J的逆矩阵，有下式给出
（5.20）
（5.21）
由三角形面积的定义可知，detJ的大小是三角形面积的两倍；加入节点1、2和3是按逆时针方向排序的，detJ值为正，则有
（5.22）
其中| |代表数量大小；大多数有限元程序都采用逆时针方向排序，并用detJ来求取面积值。
例题 5.2
确定例题5.1图中的三角形单元的Jacobian矩阵J。
解答：由题意可知
可得知detJ=23.75，它时三角形面积的两倍；假如节点1、2、3按顺时针方向进行排序，detJ将为负值。
由式（5.19）和式（5.20），可得
（5.23a）
将u替换为v，可得到类似的表达式
（5.23b）
由应变-位移关系式（5.5）、（5.12b）和（5.23），可得
（5.24a）
由的定义，可得等；因此，上一个式子可写成如下形式：
（5.24b）
写成矩阵形式为
（5.25）
其中B是（3*6）的单元应变-位移矩阵，表示三个应变和六个节点位移之间的关系，由下式给出
（5.26）
由此可看出B矩阵中的所有元素都是常数，其中每一个数都是通过节点坐标表示的。
例题 5.3
基于例题5.3图所示各三角形单元的局部节点编号，求单元的应变-节点位移关系矩阵Brsquo;。
解答：由题意可知
其中detJ的值为：，利用各三角形单元给出的节点编号，可以求出
势能方法
系统的势能Ⅱ由下式给出
（5.27）
在式（5.27）的最后一项中，i表示受集中载荷作用的点，其中；对i遍历求和即可得到所有集中载荷引起的势能。
由图5.2中所示的离散模型，可得总势能表达式如下：
（5.28a）
或
（5.28b）
其中是单元的应变能。
单元刚度
将单元的应变-位移关系式（5.25）代入单元应变能表达式（5.28b）中，可得
（5.29a）
将单元的厚度取为常数，从前面推导可知：D矩阵和B矩阵中的各项也都是常数，因此有
（5.29b）
这里，其中是单元面积，由此可得
（5.29c）
或
（5.29d）
其中是单元刚度矩阵，由下式给出
（5.30）
对于平面应力或平面应变问题，根据第一章定义的弹性矩阵D，针对相应的材料选取合适的值，通过式（5.30）的乘法运算即可得到单元刚度矩阵，注意由于D具有对称性，故也具有对称性。利用表5.1中所建立的单元节点信息，将单元刚度矩阵组装加入到对应的整体刚度矩阵K中，可得
（5.31）
整体刚度矩阵K具有对称、带状、稀疏的特性，当自由度i和j未通过一个单元相互连接时，刚度值是所有相关单元值得叠加。对于图5.2中的节点排序情况，刚度矩阵的带宽与单元中节点编号的最大差值有关，假设是一个单元e的节点编号，单元节点编号的最大差值可由下式算出
（5.32a）
半带宽可由下式得到
（5.32b）
其中NE是单元个数，2是每个节点的自由度数。
整体刚度矩阵K是在当所有的自由度Q都是自由的情况下得出的，若要处理边界条件，则需要对刚度矩阵进行相应的修改。
载荷项
首先考虑总势能表达式（5.28b）中关于体积力的表达式，将其写成分量形式
应用式（5.12a）中给出的插值关系式，可得到
（5.33）
如图5.4所示，由三角形单元形状函数的定义可知，表示一个底面积为、一个角点上高度为1（无量纲）的四面体的体积。这个四面体的体积由（）进行计算（参见图5.6），即
（5.34）
同样，可得，因此，式（5.33）可写为
（5.35）
其中frsquo;是单元的体力列阵，即
（5.36）
这些单元节点力都将集成到整体载荷列阵F中，在将frsquo;装配到整体载荷F中时，需要再次用到表5.1中的单元节点信息，单元的体力列阵frsquo;是（6*1）维的，而整体载荷列阵F是（N*1）维的，这种装配过程已经在第三章和第四章中进行过探讨，可用以下表达式示意为
（5.37）
面力定义为作用在物体表面上的分布载荷，通常施加在连接边界节点的单元边界上。作用在单元边界上的面力也将对整体载荷列阵F作出贡献，所作贡献的大小由面力表达式来计算；以边为例，假定面力的x、y方向上的分量分别为，如图5.7（a）所示，则有
（5.38）
利用下述插值关系式
（5.39）
以及
（5.40）
可推导出
（5.41）
其中Trsquo;由下式给出
（5.42）
假如是施加在直线12上的法向压力，如图5.7（b）所示，可得
其中
在式（5.42）中，可以同时考虑法向和切向分布载荷。同样在计算整体的载荷列阵F时需要考虑分布载荷的贡献。
本书所给出的程序要求载荷以集中载荷分量的形式给出。对于分布载荷，需要事先确定出等效集中载荷的分量，如以下例题所示。
例题 5.4
如例题5.4图中的二维平面，边7-8-9受线性变化的压力载荷作用，确定节点7、8、9的等效集中载荷。
解答：先分别考虑两个边界7-8和8-9，然后在进行合并。
对于边界7-8，有
这些载荷分别等效到中。
对于边界8-9，有
这些载荷分别等效到中。最后可得到等效集中载荷为
选取受到集中载荷的点作为节点，集中载荷的表达式就很容易得到。假设i为受力作用的节点，则有
（5.43）
因此，在x和y方向上的分量可以直接加到整体的载荷列阵F的第（2i-1）项和第（2i）项中去。
综上所述，体力、面力和集中力对整体的载荷列阵（或称节点力列阵）F的贡献可以通过下式来表达
结合应变能和力的表达式，总势能可表达如下
（5.44）
在处理边界条件时，则需要对刚度矩阵和节点力列阵进行修改。应用第三章和第四章的方法，可得
（5.45）
其中K和F分别是修正后的刚度矩阵和节点力列阵。为了求得位移列阵Q，可以用高斯消元法或其他方法求解方程组。
例题 5.5
全文共11540字，剩余内容已隐藏，支付完成后下载完整资料
资料编号：[144202]，资料为PDF文档或Word文档，PDF文档可免费转换为Word
以上是毕业论文外文翻译，课题毕业论文、任务书、文献综述、开题报告、程序设计、图纸设计等资料可联系客服协助查找。
