1、有限元分析及应用有限元分析及应用第二讲第二讲 平面问题有限元分析平面问题有限元分析 基本内容基本内容1 1 弹性力学基本方程弹性力学基本方程2 2 虚位移原理虚位移原理3 3 平面问题定义平面问题定义4 4 有限元分析有限元分析5 5 算例算例1 1 弹性力学基本方程弹性力学基本方程 弹性力学中的物理量弹性力学中的物理量载荷载荷：外界作用在弹性体上的力。：外界作用在弹性体上的力。体力体力：分布在整个弹性体体积内的外力。如重力：分布在整个弹性体体积内的外力。如重力和惯性力。和惯性力。TZYXF 面力面力：作用于弹性体表面上的外力，如流体压力和：作用于弹性体表面上的外力，如流体压力和 接触压力。接
2、触压力。集中力集中力：外力作用在某一点上。：外力作用在某一点上。Tzyx 应力：弹性体内某一点作用于某个截面单位面积上的应力：弹性体内某一点作用于某个截面单位面积上的内力，它反映了内力在截面上的分布密度。内力，它反映了内力在截面上的分布密度。应变：外力作用下弹性体产生的变形，包括长度和夹应变：外力作用下弹性体产生的变形，包括长度和夹角的变化。角的变化。位移：弹性体内质点位置的变化。位移：弹性体内质点位置的变化。TzxyzxyzyxTzxyzxyzyxTwvud 一、平衡方程000 xzxyxyxyyzzyzxzXxyzYxyzZxyz弹性力学基本方程zuxwzwywzvyvxvyuxuzxzy
3、zyxyx,二、几何方程二、几何方程三、物理方程三、物理方程zxzxyxzzyzyzxzyyxyxyzyxxEEEEEE)1(2),(1)1(2),(1)1(2),(1弹性力学基本方程矩阵表示二、几何方程三、本构关系四、协调方程五、应力边界条件一、平衡方程0 FAdAT0L0CD0 dd六、位移边界条件dSwZvYuXdxdydzwZvYuXWS)()(1 弹性体处于平衡状态的必要与充分条件：对于任意的、满足相容条件的虚位移，外力所做的功等于弹性体所接受的总虚变形功。总虚变形功：对于平面问题：虚位移原理虚位移原理dxdydzUxyxyzxzxyzyzzzyyxx)(dxdydsvYuXdxdy
4、vYuXxyxyyyxxS)()()(1总外力虚功：UW最小势能原理dSwZvYuXdxdydzZwYvXuVS)()(1 在几何可能的一切容许位移和形变中,真正的位移和形变使总势能取最小值；反之,使总势能取最小值者也必是真正的位移和形变。总 势 能：即：形变势能的变分表达式与虚变形功的表达式完全相同。VU 形变势能：外力势能：dxdydzUxyxyzxzxyzyzzzyyxx)(21dxdydzUxyxyzxzxyzyzzzyyxx)(形变势能变分:dSwZvYuXdxdydzwZvYuXVS)()(1外力势能变分:即：外力势能的变分表达式与外力虚功负值的表达式完全相同。平面问题定义平面问题
5、定义 平面应力问题平面应力问题几何条件：厚度尺寸远远小于截面尺寸，即结构形状呈薄几何条件：厚度尺寸远远小于截面尺寸，即结构形状呈薄板型板型载荷条件：载荷平行于板平面且沿厚度方向均匀分布，而载荷条件：载荷平行于板平面且沿厚度方向均匀分布，而板平面不受任何外力作用。板平面不受任何外力作用。00zyzxz 平面应变问题平面应变问题 几何条件：沿厚度方向的截面形状和大小相同且厚度尺寸几何条件：沿厚度方向的截面形状和大小相同且厚度尺寸远远大于截面尺寸，即结构呈等截面的细长形。远远大于截面尺寸，即结构呈等截面的细长形。载荷条件：载荷垂直于厚度方向（平行横截面）且沿厚度均载荷条件：载荷垂直于厚度方向（平行横
6、截面）且沿厚度均匀分布，两个端面不受力。匀分布，两个端面不受力。yxzzyzxzxyzz00一、结构离散一维二维三维有限元分析 以及，一系列结构单元模型，如：杆单元梁单元管单元膜单元平面单元应当注意实际物体与离散后的单元模型的联系与差异，如：以及，一系列特殊单元模型，如：裂纹单元无界单元波纹板单元选择后，按照单元类型划分实际物体平面问题的四边形单元和三角形单元划分通常把三维实体划分成4面体或6面体单元的网格，平面问题划分成三角形或四边形单元的网格。空间问题的四面体形单元划分空间问题的六面体形单元划分二、位移函数二、位移函数 定义：定义：假设一种函数来近似表示单元内部的实假设一种函数来近似表示单
7、元内部的实际位移分布，该函数称为位移函数。（际位移分布，该函数称为位移函数。（分片插分片插值思想值思想）可以考虑采用多项式作为近似位移函数，多项式函数运算简便，且随其项数的增多，理论上可以逼近任何精确函数。收敛准则 多项式位移模式阶次的选择一、收敛准则1、位移模式必须包含单元的刚体位移满足条件1、2的单元为完备单元二、多项式位移模式阶次的选择按照帕斯卡三角形选2、位移模式必须能包含单元的常应变3、位移模式在单元内要连续、并使相邻单元间的位移必须协调满足条件3的单元为协调单元几何各向同性：位移模式应与局部坐标系的方位无关帕斯卡三角形多项式应有偏惠的坐标方向，多项式项数等于单元边界结点的自由度总数
8、。1yx22yxyx3223yxyyxx432234yxyyxyxx54322345yxyyxyxyxx常数项线性项二次项三次项四次项五次项对称轴 首先，建立首先，建立x xy y坐标系，以三角形单坐标系，以三角形单元为研究对象，其结点编码为元为研究对象，其结点编码为i i，j j，m m，并规定以逆时针方向编码为正向，每个并规定以逆时针方向编码为正向，每个结点可以有结点可以有2 2个位移分量（个位移分量（u ui i，v vi i；u uj j，v vj j；u uk k，v vk k），也称为具有），也称为具有2 2个结点自由度，个结点自由度，如图所示。如图所示。3结点三角形单元结点三角形
9、单元kkkkkyaxaavyaxaau654321三节点三角形单元位移函数三节点位移可表述为kkkkkkjjjjjjiiiiiiyaxaavyaxaauyaxaavyaxaauyaxaavyaxaau654321654321654321kkkjjjiiiyxuyxuyxua211kkjjiiyuyuyua111212kkjjiiuxuxuxa111213kkjjiiyxyxyx1112 kkkkjjjjiiiiuycxbauycxbauycxbau21kkkkjjjjiiiivycxbavycxbavycxbav21kjikjikjikjikjiyyyxxxcccbbbaaa111kjikji
10、kkjjixxcyybyxyxa11,11,式中：),(21kjiycxbaNiiii形函数I 二阶单位阵二阶单位阵,N 形函数矩阵形函数矩阵 eekjiNINININvudiikkjjiiiikkjjiivNvNvNvNvuNuNuNuNuijyxNyxNjjiiii0),1),.1（形函数jjiijjiivNvNvuNuNuij上位移：边界形函数的性质11.2kjiNNN函数和为在单元任一点上三个形0.3kNijijk上的一边在三角形单元ijiixxxxyxN1),(ijijxxxxyxN),(0),(yxNkyx0knji ekkjjiikjikjixyyxbcbcbccccbbbxvy
11、uyvxu00000021),(0021kjibccbBiiiii三、应变 eB kjiBBBB 应变矩阵应变矩阵为常量，单元内应变是常数应变矩阵为常量，单元内应变是常数四、应力 2100010112ED阵：平面应力问题，弹性矩 eeSBDD kjikjiSSSBBBDBDS应力矩阵 应变矩阵为常量，单元内应力也是常数，相邻单元的应变矩阵为常量，单元内应力也是常数，相邻单元的应变与应力将产生突变，但位移确是连续的。应变与应力将产生突变，但位移确是连续的。)1(221010110111)21)(1()1(称对ED平面应变问题五、单元刚度矩阵 Tkkjjiievuvuvu单元结点虚位移eNd单元内
12、虚位移 eB单元内虚应变 eTTeTtdxdyBDBtdxdyU内力虚功 eTeRW外力虚功kjiskjir,;,其中 eeeTekRtdxdyBDBR tBDBktdxdyBDBkTT kkkjkijkjjjiikijiikkkkkkkkkksrsrsrsrsrsrsrsrrsbbcccbbcbccbccbbEtk21212121)1(42UW虚功原理对于平面应力问题：为单元厚度t为单元面积结构整体刚度矩阵 通过对弹性体离散化总位能中的单元刚度矩阵和结点等效载荷列阵实行尺寸扩充变换，增加其阶数和行数，就可形成一个简洁非离散形式的总位能表达式，而且总位能p值保持不变。平面三角形三结点单元的刚度
13、矩阵为6阶矩阵，结点等效载荷列阵为6行列阵，这与单元三结点共有23=6个位移相对应，也就是说，单元刚度矩阵的阶数和结点等效载荷列阵的行数要与单元结点位移数相等。同样，具有n个离散结点的平面弹性体的总刚度矩阵是一个2n阶的矩阵。因此，其单元刚度矩阵和结点等效载荷列阵分别要扩充为2n阶和2n行。以如图所示的有5个结点和3个单元的弹性平板受集中力p作用情况为例，说明尺寸扩充的具体做法。53412123p5结点结点3单元平板的有限元模型单元平板的有限元模型该弹性平板的该弹性平板的3个单元分别由个单元分别由3组结点构成：组结点构成：第一单元（结点第一单元（结点1、2、3）；）；第二单元（结点第二单元（结
14、点1、3、4）；）；第三单元（结点第三单元（结点3、5、4）。）。该平板的总位能表达式可写成 33T22T11T333T222T111T31eeeT31eeeeT31eeeT31eeeT31eeeeT31ee212121PaPaPaaKaaKaaKaPaaKaPaPaaKasfpp 第一、二、三单元结点位移列阵为TTT445533344331123322111vuvuvuavuvuvuavuvuvua 第一、二、三单元刚度矩阵为（行列序号为单元结点序号）44345343354355353334335333334424324123423323121421321122331321311231221
15、2111311211111KKKKKKKKKKKKKKKKKKKKKKKKKKKKKK 第一、二、三单元结点载荷列阵为 从各单元矩阵可以看出：矩阵元素都以所属结点序号为下标，以所属单元序号为上标；TTT434353533333342423232121223131212111111yxyxyxyxyxyxyxyxyxppppppPppppppPppppppP 单元矩阵和列阵尺寸扩充的做法：按结点数的2倍确定所有矩阵和列阵的阶数和行数，本弹性体共有5个结点，因此阶数和行数为10；单元结点位移列阵改写为本弹性体全部结点的10阶位移列阵T5544332211vuvuvuvuvua 单元结点载荷列阵扩充
16、改写为10阶的载荷列阵 321T3T2T1535343433333424232321212313121211111000000000000PPPPPPPyxyxyxyxyxyxyxyxyxpppppppppppppppppp对扩充为10阶的单元刚度矩阵元素实行“对号入座”式的填写，即只填写与所属结点序号对应的行列位置上的元素，其余位置上的元素写零，具体做法（按22子块方式）5535435334534434333533433333442432412342332312142132112233132131123122121113112111110000000000000000000000000000
17、0000000000000000000KKKKKKKKKKKKKKKKKKKKKKKKKKKKKK改写总位能式为PaKaaTTT321TTTT3T2T1T33T22T11T333T222T111T31eeeT31eeeeT31eeeT31eeeT31eeeeT31ee212121212121321 PaaKKKaPaPaPaaKaaKaaKaPaPaPaaKaaKaaKaPaaKaPaPaaKasfpp 上式中，称为结构整体刚度矩阵，即 虽然总体单元数和结点数很多，结构刚度矩阵的阶数很高，但刚度系数中非零系数却很少，这就是刚度矩阵的大型和稀疏性。553543533453443243324123
18、5334323332132131123122121114213211211121321000000KKKKKKKKKKKKKKKKKKKKKKK六、等效结点力、载荷列阵 TkkjjiieYXYXYXR单元结点等效力 tdxdypdtdsqdGdRTTTeTe eeeePQFR功。等于单元上外力所做虚等效结点力所做的虚功eNd单元内虚位移 tdxdypNtdsqNGNRTTTTeeTe tdxdypNPtdsqNQGNFTeTeTe体积力等效结点力表面力等效结点力集中力等效结点力 （1）均质等厚单元的自重 单元的单位体积重量为g，坐标方向如图所示。三角形单元作用体积力 几种常见载荷表达式如下 应
19、有 其中，每个结点的两个等效载荷分量的列矩阵是 （10-29）tdxdyggemjimjief00NNNPPPPfi,j,mgtAyxgtNyxtgNNPPiiiiyixi310dd0dd000eeP 式中 则积分i,j,mycxbaANiiii21eeeeeedddddd2dd2dd21ddiiiiiiiiiiyxycyxxbyxaAgtyxycxbaAgtyxgtycxbaAyxgtN 通过选择三角形单元重心为坐标原点，则有 并且由于ee0ddddiiyxycyxxbAaaamji32 则有积分 单元自重的等效结点载荷是 332222ddegtAAAAgtdxdyaAgtdxdyaAgty
20、xgtNeeiiiT10101031gtAfP（2）均布侧压 侧压q作用在边ij上，q以压为正，如图所示。单元边上作用均布侧压设ij边长为1，与x轴的夹角为，侧压q在x和y方向的分量qx和qy为作用在单元边界上的面积力为 在单元边界上可取局部坐标s，沿ij边插值函数可写作ijyjixxxlqqqyylqqqcossinijjiyxxxyylqqqT01mjiNlsNlsN10.3 结点等效载荷 侧压作用下的单元等效结点载荷为因此 022221mymxijjyjilxlxjjxijiyjilxlxiixpxxqtyyqttdsqlstdsqNxxqtyyqttdsqls-tdsqNPPPPPTs
21、0021ijjiijjixxyyxxyyqtP （3）x向均布力 均布力q作用在ij边，如图所示。单元边上作用的x方向均布力 这时边界上面积力为 单元等效结点载荷平均分配给ij结点为 0qTT00010121eqltP（4）x方向三角形分布载荷载荷作用在ij边，如图所示。单元ij边上作用的x方向三角形分布力 这时边界上面积力写作局部坐标s的函数，即 则单元等效结点载荷为 01qls-TTe0003103221qltP引入位移边界条件 有限元法中，几何边界条件的形式是给定若干个边界结点上的位移值，即 可以是零值或非零值。这样一来有限元方程中出现了恒等式，可以采取简单地去掉这些恒等式的方法，求解只
22、含有未知结点位移的有限元方程，对于刚度矩阵，只需划掉与已知结点位移对应的行列，对于结点位移和结点载荷列阵只需划掉对应行即可。但那样做却不利于编程，因此通常是采取一种不用改变原方程尺寸的对角元素乘大数法引入边界条件。jjaaja 具体做法是，当有结点位移为给定值aj=时，第j个方程作如下修改：对角元素Kjj乘以大数（可取1010左右量级），并将Pj用Kjj取代，即 njjjnjnnnnjnjjjjnnpaaKppaaaaKKKKaKKKKKKKKK212121212222111211 经过修改后的第j个方程为 由于KjjKij（ij），方程左端的Kiiaj项较其他项要大得多，因此近似得到 则有j
23、jjnjnjjjjjaaKaKaaKaKaK2211jjjjjjaaKaaKjjaa对于多个结点（j=c1，c2，cl）的给定位移时，则按序将每个给定位移都作上述修正，得到全部进行修正后的K和P，然后解方程，即可得到包括给定位移在内的全部结点位移值。这个方法使用简单，对任何给定位移（零值或非零值）都适用。采用这种方法引入强制边界条件时方程阶数不变，结点位移顺序不变，编制程序十分方便，因此在有限元法中经常采用。有限元求解方程在引入位移边界条件、消除了K矩阵的奇异性后，就可以从它解得结构的结点位移a。并进一步回到单元中，用已知的位移，求得各个单元的应变和应力。【例】高深悬臂梁平面问题的有限元分析图
24、所示为一高深悬臂梁，在右端部受集中力F作用，材料弹性模量为E、泊松比=1/3，悬臂梁的厚度(板厚)为t，不计体力，试按平面应力问题计算各个节点位移及支座反力。(a)问题描述 (b)有限元分析模型 右端部受集中力作用的高深梁【解】对该问题进行有限元分析的过程如下。（1）结构的离散化与编号 对该结构进行离散，单元编号及结点编号如图10-8（b）所示，即有两个3结点三角形单元（2、3、4和3、2、1）。载荷F按静力等效原则向结点1和结点2移置。结点位移列阵a 结点外载荷列阵FT44332211vuvuvuvua 约束的支反力列阵R 总的结点载荷列阵P 其中 R x3，R y3和R x4，R y4 分
25、别为结点3和结点4的两个方向的支反力。T00002020FFFT34330000yxyxRRRRRT44332020yxyxRRRRFFP （2）各个单元刚度矩阵 若二单元各取图示局部编码（都以最小角顶点为结点i，另一锐角顶点为结点j）时，其单元刚度矩阵将完全相同，即3133443231323437323432143240032323403432031320323103213200132911212213221222223231232233224414314213413313212412312211EtKKKKKKKKKKKKKKKKKKKK代入有限元方程Ka=P中，即有44334433221
26、12020313344323132003437323432100432313003431323234037340321313203431304323213400373234003132432313340032132343437329yxyxRRRRFFvuvuvuvuEt（4）引入边界条件求解有限元方程 该实例的位移边界条件为 将其代入有限元方程中，划去上述已知结点位移所对应的第5至第8小行小列，则有00004433vuvu由此式可求出结点位移如下：20203130432037323443231334323434373292211FFvuvuEt4285019988812211.EtFvuvu
27、（5）支反力的计算将所求得的结点位移式代入总刚度有限元方程中，可求得支反力如下：上述支反力与外载荷构成一个平衡力系。FvuEtRFvuEtRFuvuEtRFvvuEtRyxyx07.1313232923232907.03431323292343232922422421132113 【例】中心具有圆孔的方板，上下两边受y方向均匀拉伸载荷，如图所示，方板边长L=8cm，圆孔半径r=1cm，板厚t=0.1cm，载荷集度q=10Mpa，材料常数E=2.0105MP，=0.3，现采用有限元法进行应力分析，目的是得到孔边的应力集中系数。【解】由于问题的对称性，取板的1/4建立有限元模型。边界条件是：x=0
28、边界上，u=0；y=0边界上，v=0；y=4cm边界上，有沿y方向的均布载荷作用。图所示的是结点数为1110的有限元网格图。中心圆孔方板受均匀拉伸 图（a），（b）所示是用1110结点数的网格计算得到的y（沿y=0边）和x（沿x=0边）的分布图。图中还给出了具有中心小圆孔无穷大板，及无孔板的理论解，以便比较。应力沿坐标轴的分布【例】计算受集中力作用的三角形平板（见图）。已知集中力P=10，弹性模量E=1，泊松比=0.25，板厚t=1，不计体力。（a）三角形平板 （b）有限元模型受集中力作用的三角形平板P=10P=10 xy12345678910123456789【输入数据】结点数目：10 单元
29、数目：9 平面应力问题 PTYPE=1（平面问题类型标志）弹性模量E=1 泊松比=0.25 比重G=0 平板厚度t=1表1 结点坐标结点号码 x y1 +0.0 +0.02 +2.0 +0.03 +4.0 +0.04 +6.0 +0.05 +1.0 +2.06 +10.0 +2.07 +5.0 +2.08 +2.0 +4.09 +4.0 +4.010 +10.0 +6.0矩形双线性单元Tevuvuvuvu44332211TeYXYXYXYXR44332211矩形单元矩形单元结点位移、结点力列阵一、位移模式与形函数87654321aaaavaaaauyx0a1abb234o11234111111
30、1byax,正方形规则单元正方形单元与矩形单元的关系41443322114144332211iiiiiivNvNvNvNvNvuNuNuNuNuNu形函数的性质：本点处值为1，它点处值为00),(0),(0),(1),(441331221111NNNN)1)(1(),()1)(1(),()1)(1(),()1)(1(),(44332211NNNN414321根据形函数的性质有：4,3,2,1,4)1)(1(0000iNiii其中：统一公式：eNvud exyyxBBBBvbuavaubabvaubvbuaxvyuyvxu432111111二、应变 eB4,3,2,1)1()1()1(00)1(
31、410010000ibaababNbNaNaNbabBiiiiiiiii三、应力 eeSBDD 43214321SSSSBBBBDBDS应力矩阵4,3,2,1)1(21)1(21)1()1()1()1()1(40000002ibaabababESiiiiiii平面应力问题四、单元刚度矩阵 Tevuvuvuvu44332211单元结点虚位移：eNd单元内虚位移：eB单元内虚应变：eTTeTdabtdBDBtdxdyU内力虚功：eTeRW外力虚功：4,3,2,1,ji其中：eeeTekRdtabdBDBR ddBDBtabkT 44434241343332312423222114131211kkk
32、kkkkkkkkkkkkkk dtabdBDBkjTiijUW虚功原理为单元厚度t子矩阵为：第4节 六结点三角形单元一、位移模式与形函数2161353244,4,4)3,2,1)12(LLNLLNLLNiLLNiii边中点：（角结点：取三角形顶点和边中点作结点，位移模式为：2121121098726524321yaxyaxayaxaavyaxyaxayaxaau六结点三角形单元yx0321645用面积坐标表示的形函数为：iiiiiivNvuNu6161,exyyxBBBBBBxvyuyvxu654321二、应变)3,2,1()14()14()14(00)14(21iLbLcLcLbBiiiii
33、iiii eNvud其它略自己推导,)6,5,4(iBi第4节 十结点三角形三次单元确定位移模式和形函数3211011393138332723262215121427)13(29),13(29)13(29),13(29)13(29),13(29)3,2,1)23)(13(21LLLNLLLNLLLNLLLNLLLNLLLNLLLNiLLLNiiii（取三角形各边三分点和面积坐标相等的内点作为结点十结点三角形单元。十结点三角形单元yx032164578910第5节 四结点四边形等参单元一、母单元的形函数母单元iiiiiiyNyxNx4141,三、位移模式 eNvudo12341111yx03214四边形单元4,3,2,1,4)1)(1(0000iNiii其中：二、坐标变换iiiiiivNvuNu4141,由此可知：单元的位移场和单元形状用相同的形函数，故称等参数单元（等参元）四、导数的坐标变换yxJyxyxyxyyxxyyxx xxyyJyxdet1iiiiiiiiiiiiyNyxNxxNyxNx,其中：
