下载文档原格式
1.题目概况
矩形板尺寸如下图1，板厚为5mm。
材料弹性模量为
松比μ= 0.27 。
施加约束和载荷并讨论：
图
1 计算简图
1.1基本数据
E = 2⨯105N/mm2，泊
序号载荷约束备注42 向下集中载荷F=800N, 作用于cd 边3/4 处（近d） c d 点简支
1.2分析任务/分析工况
讨论板上开孔、切槽等对于应力分布的影响。
（载荷约束组合不变）。
提示：各种圆孔，椭圆孔随大小、形状、数量，分布位置变化引起的应力分布变化；各种形状，大小的切槽及不同位置引起应力分布的变化等，选择二至三种情况讨论，并思考其与机械零部件的构型的相对应关系。
2.模型建立
2.1单元选择及其分析
由于平板长宽分别为300x100，故可取网格单元大小为1。
如图：
2.2模型建立及网格划分
模型按单元为1 划分后的网格大小如图所示：
2.3载荷处理
向下集中载荷F=800N, 作用于cd 边3/4 处（近d） c d 点简支
3.计算分析
3.1位移分布及其分析
（1）位移分布如图：。
250250试题 5：图示为带方孔（边长为 80mm ）的悬臂梁，其上受部分均布载荷（P=10KN/m ）作用，试采用一种平面单元，对图示两种结构进行有限元分析，并就方孔的布置（即方位）进行分析比较，如将方孔设计为圆孔，结果有何变化？（板厚为 1mm ，材料为钢）。
3001KN9003001KN图6-1一、几何建模与分析由图6-1及问题描述可知，板的长宽尺寸远远大于厚度，研究结构为一很薄的等厚度薄板，满足平面应力的几何条件；作用于薄板上的载荷平行于板平面且沿厚度方向均匀分布，而在两板面上无外力作用，满足平面应力的载荷条件。
故该问题属于平面应力问题，薄板所受的载荷为面载荷，分布情况及方向如图6-1所示，建立几何模型，进行求解。
薄板的材料为钢，则其材料参数：弹性模量E=2.1e11，泊松比σ=0.3。
二、有限元分析及其计算结果选取PLANE182作为分析的单元，来分析薄板的位移和应力，由于此问题是平面应力问题，并在K3选择str w/thk ，设置THK 为1。
1）方孔竖直制，划分方式采用自由方式，划分后网格的模型如图6-2所示。
计算得到的位移和应力分布如图6-3所示。
图6-2 方孔竖直的网格划分图6-3 位移及应力分布云图2）方孔正直制，划分方式采用自由方式，划分后网格的模型如图6-4所示。
计算得到的位移和应力分布如图6-5所示。
图6-4 方孔正直的网格划分图6-5 位移及应力分布云图3）圆孔按图6-1所示模型进行建模。
并用PLANE182单元进行划分网格，网格大小采用全局网格控制，划分方式采用自由方式，划分后网格的模型如图6-6所示。
计算得到的位移和应力分布如图6-7所示。
图6-4 方孔正直的网格划分图6-5 位移及应力分布云图根据以上的模型分析的位移和应力图，可以得出方孔竖直、方孔正直、圆孔的最大最小位移应力的分布如表6-1所示。
三、比较与分析1）方孔竖直与方孔正直的比较，发现方孔正直的位移变形较小，应力相差不大2）圆孔与方孔比较，发现圆孔的位移变性最小，应力也最小，故可以得出圆孔的布置结构对整体布置的效果最好。
超静定梁的有限元分析本文分别通过材料力学解法和有限元解法，求出了超静定梁的支反力、最大位移及最大位移出现位置，并对两者进行了比较和误差分析。
一、超静定梁的材料力学解法梁的约束反力数目超过了有效平衡方程数，单纯使用静力平衡不能确定全部未知力的梁称为超静定梁。
超静定梁比静定梁有许多优点，如可用较少材料获得较大的刚度和强度，个别约束破坏后仍可工作等。
因而超静定梁在工程中得到较多的应用。
超静定梁的解法有很多种，本文采用力法的一种——变形比较法求解未知量。
图1图2选取C 点的支座为多余约束，Rc 为多余支座反力，则相应的基本静定梁为一外伸梁，如图2所示，其上受集中载荷P 、均布载荷q 和多余支座反力Rc 的作用。
相应的变形条件为：c cP cq cRc f f f f =++=其中316cP B Pl f l EI θ=⨯= 4724cq ql f EI =-323c cRc R l f EI =则316Pl EI 4724ql EI -+323c R l EI=0 将已知数据带入可求得 6.25c R =- 负号表示c R 的方向与假设的方向相反。
再列出平衡方程：0X =∑AX R =0A M =∑ 232022B C ql Pl R l R l ---=0C M =∑ 232022AY B ql PllR R l +--=带入已知条件求得：AX R = 393.75AY R = 812.5B R =用叠加法求最大位移：最大的向下位移在A 与B 两点中间：334410.7910481632C R l Pl ql f EI EI EI -=-++=-⨯最大的向上位移在B 与C 两点中间：3344213490.22525103248512C R l Pl ql f EI EI -=--=⨯二、超静定梁的有限元解法在ANSYS 平台上，求解超静定梁。
建模、单元划分、加载后结果如图3所示。
图3求解后可以通过图形和列表两种方式查看结果。
1．推导有限元计算格式，理解有限元原理：建立图示受拉直杆在自重（设单位长度重度为q ，截面积为A ）和外力P 作用下的拉伸问题的微分方程，并分别利用不同的原理（变分求极值（最小势能或虚功原理）、加权残值法）推导有限元计算格式（取两个单元）。
手工求出端点的位移（自己给定参数值）。
设杆长为L ，截面面积为A(x)，弹性模数为E,单位长重量q ，受拉杆x 处的位移为u(x)。
取微元dx 的力平衡，建立受拉杆位移所满足的微分方程()du x dx ε=，()du x E E dxσε== dx 上下截面内力与微元自重相等得()*()()*()A x dx x dx A x x dx qdx σσ++-+=-(()())dA x x q dxσ∴=- (())d duEA x q dx dx=- 0x L << ()0u x = 0x =()duEA x p dx= x L = 得解析解：2()2q x P u Lx x EA EA=-+将其分为两个单元，节点为1,2,3，得22382qL PL u EA EA=+232qL PL u EA EA=+有限元法：1)位移函数01u α= 2111u u l α-=得1211(1)x x u u u l l =-+ 令11(1)x N l =-21x N l = 11122122u u N u N u N N u⎧⎫⎪⎪⎡⎤=+=⎨⎬⎣⎦⎪⎪⎩⎭{}1u N d ⎡⎤=⎣⎦ 2)应变、应力表达{}{}111211du dN d d dx dx l l ε⎡⎤⎡⎤===-⎢⎥⎣⎦⎣⎦{}1B d ε⎡⎤=⎣⎦ {}1E E B d σε⎡⎤==⎣⎦ {}1S d σ⎡⎤=⎣⎦3)势能表示{}{}(){}{}(){}{}{}{}{}1111''112211''121112210111111111111111121221222T V ll T T T T T U W D dV F u F u qdx u u d B E d Adx F u F u ql EA EA ql l l d d d F d EA EA ql l l εε⎡⎤=-=-+-⎣⎦+⎡⎤=-+-⎣⎦⎡⎤⎡⎤-⎢⎥⎢⎥⎢⎥⎢⎥=--⎢⎥⎢⎥⎢⎥-⎢⎥⎣⎦⎢⎥⎣⎦∏⎰⎰⎰4)单元平衡方程 a)最小势能原理110u ∂=∂∏120u ∂=∂∏111111212112112ql F u AE l u ql F ⎧⎫-⎪⎪⎧⎫⎡⎤-⎪⎪⎪⎪=⎨⎬⎨⎬⎢⎥-⎣⎦⎪⎪⎪⎪⎩⎭+⎪⎪⎩⎭b)虚位移原理{}(){}(){}TeTdd F qdx d δδδεσΩ+=Ω⎰⎰{}{}1B d σεδ⎡⎤=⎣⎦ {}1E E B d σεδ⎡⎤==⎣⎦{}(){}{}(){}111111TTT l d F d B E B d Adxδδ⎡⎤=⎣⎦⎰ 由虚位移任意性得，{}{}1111T lF B E B Adxd ⎡⎤=⎣⎦⎰ 积分得111111212112112ql F u AE l u ql F ⎧⎫-⎪⎪⎧⎫⎡⎤-⎪⎪⎪⎪=⎨⎬⎨⎬⎢⎥-⎣⎦⎪⎪⎪⎪⎩⎭+⎪⎪⎩⎭ 记为{}{}111k d F ⎡⎤=⎣⎦ 同理222212323112112ql F u AE l u ql F ⎧⎫-⎪⎪⎧⎫⎡⎤-⎪⎪⎪⎪=⎨⎬⎨⎬⎢⎥-⎣⎦⎪⎪⎪⎪⎩⎭+⎪⎪⎩⎭{}{}222k d F ⎡⎤=⎣⎦ {}{}ei i eF R =∑ 12220F F += 23F P =11111112211223222022202EAEAql F l l u ql ql EA EA EA EA u l l l l u ql EAEA P l l ⎡⎤⎧⎫-⎢⎥+⎪⎪⎢⎥⎧⎫⎪⎪⎢⎥⎪⎪⎪⎪⎢-+-⎥=+⎨⎬⎨⎬⎢⎥⎪⎪⎪⎪⎢⎥⎩⎭⎪⎪+⎢⎥⎪⎪--⎢⎥⎩⎭⎣⎦可得：22382qL PLu EA EA=+232qL PL u EA EA=+与解析解结果一致。
有限元分析大作业报告一、引言有限元分析是工程领域中常用的数值模拟方法，通过将连续的物理问题离散为有限个子区域，然后利用数学方法求解，最终得到数值解。
有限元分析的快速发展和广泛应用，为工程领域提供了一种强大的工具。
本报告将介绍在大作业中所进行的有限元分析工作及结果。
二、有限元模型建立本次大作业的研究对象是工程结构的应力分析。
首先，通过对结构进行几何建模，确定了结构的尺寸和形状。
然后，将结构离散为有限个单元，每个单元又可以看作一个小的子区域。
接下来，为了求解结构的应力分布，需要为每个单元确定适当的单元类型和单元属性。
最后，根据结构的边界条件，建立整个有限元模型。
三、材料属性和加载条件在建立有限元模型的过程中，需要为材料和加载条件确定适当的参数。
本次大作业中，通过实验获得了结构材料的弹性模量、泊松比等参数，并将其输入到有限元模型中。
对于加载条件，我们选取了其中一种常见的加载方式，并将其施加到有限元模型中。
四、数值计算和结果分析为了求解结构的应力分布，需要进行数值计算。
在本次大作业中，我们选用了一种常见的有限元求解器进行计算。
通过输入模型的几何形状、材料属性和加载条件，求解器可以根据有限元方法进行计算，并得到结构的应力分布。
最后，我们通过对计算结果进行分析，得出了结论。
五、结果讨论和改进方法根据计算结果，我们可以对结构的应力分布进行分析和讨论。
根据分析结果，我们可以得出结论是否满足设计要求以及结构的强度情况。
同时，根据分析结果，我们还可以提出改进方法，针对结构的特点和问题进行相应的优化设计。
六、结论通过对工程结构进行有限元分析，我们得到了结构的应力分布，并根据分析结果进行了讨论和改进方法的提出。
有限元分析为工程领域提供了一种有效的数值模拟方法，可以帮助工程师进行结构设计和分析工作，提高设计效率和设计质量。
【1】XXX，XXXX。
【2】XXX，XXXX。
以上是本次大作业的有限元分析报告，总结了在建立有限元模型、确定材料属性和加载条件、数值计算和结果分析等方面的工作，并对计算结果进行讨论和改进方法的提出。
有限元大作业一题目要求：图1所示为一悬臂梁，在端部承受载荷，材料弹性模量为E，泊松比为1/3，悬臂梁的厚度（板厚）为t，若该粱被划分为两个单元，单元和节点编号如图所示，试按平面应力问题计算各个节点位移计支反力。
一、单元划分1.计算简图及单元划分如下所示：2.进行节点及单元编号节点i j m单元① 2 3 4② 3 2 13.节点坐标值节点号1 2 3 4坐标值X 2 2 0 0Y 1 0 1 0二、计算单元刚度矩阵1、计算每个单元面积△以及i b ，i c （m j i i ,,=） ①②单元的面积相等，即12121=⨯⨯=∆ 单元①的i b ，i c⎩⎨⎧=--==-=0)(1m j i m j i y x c y y b ⎩⎨⎧=--==-=2)(0i m ji m j x x c y y b ⎩⎨⎧-=--=-=-=2)(1j i mj i m y x c y y b 对平面应力问题，其表达式为[]⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡-+-+-+-+∆-=s r s r sr s r s r s r s r s r b b uc c cb u b uc b c u c ub c c u b b u Et Krs 21212121)1(42 然后对单元①求解单元刚度子矩阵2==i r 2==i s []⎥⎦⎤⎢⎣⎡=3/1001329)1(22Et K 2==i r 3==j s []⎥⎦⎤⎢⎣⎡=03/23/20329)1(23Et K2==i r 4==m s []⎥⎦⎤⎢⎣⎡----=3/13/23/21329)1(24Et K 3==j r 3==j s []⎥⎦⎤⎢⎣⎡=4003/4329)1(33Et K 3==j r 2==i s []⎥⎦⎤⎢⎣⎡=03/23/20329)1(32Et K 3==j r 4==m s []⎥⎦⎤⎢⎣⎡----=43/23/23/4329)1(34Et K 4==m r 4==m s []⎥⎦⎤⎢⎣⎡=3/133/43/43/7329)1(44Et K 4==m r 2==i s []⎥⎦⎤⎢⎣⎡----=3/13/23/21329)1(42Et K 4==m r 3==j s []⎥⎦⎤⎢⎣⎡----=43/23/23/4329)1(43Et K由子矩阵[]e rs K 合成单元刚度矩阵[]⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡----------------=3/133/443/23/13/23/43/73/23/43/2143/24003/23/23/403/43/203/13/203/23/103/213/2001329)1(Et K将单元①的单元刚度矩阵补零升阶变为单元刚度矩阵，其在总体刚度矩阵中的位置为：节点号→单元②的i b ，i c⎩⎨⎧=--=-=-=0)(1m j im j i y x c y y b ⎩⎨⎧-=--==-=2)(0i m ji m j x x c y y b ⎩⎨⎧=--==-=2)(1j i mj i m y x c y y b 然后对单元 求解单元刚度子矩阵：3==i r 3==i s []⎥⎦⎤⎢⎣⎡=3/1001329)2(33Et K 3==i r 2==j s []⎥⎦⎤⎢⎣⎡=03/23/20329)2(32Et K 3==i r 1==m s []⎥⎦⎤⎢⎣⎡----=3/13/23/21329)2(31Et K 1 2 3 412[])1(22K[])1(23K[])1(24K3[])1(32K[])1(33K[])1(34K4[])1(42K[])1(43K[])1(44K2==j r 2==j s []⎥⎦⎤⎢⎣⎡=4003/4329)2(22Et K 2==j r 3==i s []⎥⎦⎤⎢⎣⎡=03/23/20329)2(23Et K 2==j r 1==m s []⎥⎦⎤⎢⎣⎡----=43/23/23/4329)2(21Et K 1==m r 1==m s []⎥⎦⎤⎢⎣⎡=3/133/43/43/7329)2(11Et K 1==m r 3==i s []⎥⎦⎤⎢⎣⎡----=3/13/23/21329)2(13Et K 1==m r 2==j s []⎥⎦⎤⎢⎣⎡----=43/23/23/4329)2(12Et K 由子矩阵[]e rs K 合成单元刚度矩阵[]⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡----------------=3/133/443/23/13/23/43/73/23/43/2143/24003/23/23/403/43/203/13/203/23/103/213/2001329)2(Et K将单元②的单元刚度矩阵补零升阶变为单元贡献矩阵，其在总体刚度矩阵中的位置为：节点号→1 2 3 41 [])2(11K[])2(12K[])2(13K2 [])2(21K[])2(22K[])2(23K3 [])2(31K [])2(32K [])2(33K 4三、计算总体刚度矩阵总体刚度矩阵是由各单元的贡献矩阵迭加而成)2()1(][][][][K K K K e +==∑四、进行节点约束处理根据节点约束情况，在总刚矩阵中可采用划行划列处理约束的方法，由题目易知，节点3和4的已知水平位移和垂直位移都为零，划去其相对应的行和列，则总刚矩阵由8阶变为4阶，矩阵如下：⎪⎪⎭⎪⎪⎬⎫⎪⎪⎩⎪⎪⎨⎧--=⎪⎪⎭⎪⎪⎬⎫⎪⎪⎩⎪⎪⎨⎧⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡--------2/02/03/13043/203/73/23/443/23/133/43/23/43/43/73292211p p v u v u Et329][Et K =1 2 3 413/133/43/43/743/23/23/4----3/13/23/21----000243/23/23/4----3/13003/73/43/403/13/23/21----33/13/23/21----3/43/403/13003/743/23/23/4----40003/13/23/21----43/23/23/4----3/133/43/43/7化简⎪⎪⎭⎪⎪⎬⎫⎪⎪⎩⎪⎪⎨⎧--=⎪⎪⎭⎪⎪⎬⎫⎪⎪⎩⎪⎪⎨⎧⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡--------Et p Et p v u v u 3/1603/160130122072412213424472211 五、求解线性方程组方法：采用LU 分解法 1.求解矩阵[]U 各元素⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡-----⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡------⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡--------75/10775/640075/6475/353007/767/27/7502447~7/877/87/7607/87/337/207/767/27/7502447~13012207241221342447⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡----353/44900075/6475/353007/767/27/7502447~ 得到的[]U 矩阵如下:[]⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡----=353/44900075/6475/353007/767/27/7502447U 2.求解矩阵[]L 各元素⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡--------=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡----⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡----13012207241221342447353/44900075/6475/353007/767/27/75024471353/6475/767/20175/27/40017/40001 得到的[]L 矩阵如下：[]⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡--------=13012207241221342447L3.进行求解⎪⎪⎭⎪⎪⎬⎫⎪⎪⎩⎪⎪⎨⎧---=⇒⎪⎪⎭⎪⎪⎬⎫⎪⎪⎩⎪⎪⎨⎧--=Et p Et p Et p y Et p Et p Ly 79425/850800225/323/1603/1603/160⎪⎪⎭⎪⎪⎬⎫⎪⎪⎩⎪⎪⎨⎧---=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡----⇒=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡Et p Et p Et p v u v u y v u v u U 79425/850800225/323/160353/44900075/6475/353007/7675/27/750244722112211 解得Et p v /422.82-= Et p u /497.12-= Et p v /028.91-= Et p u /897.11=于是求得各节点的位移为：⎩⎨⎧-==Etp v Etp u /028.9/897.111 ⎩⎨⎧-=-=Etp v Etp u /422.8/497.122 ⎩⎨⎧==033v u ⎩⎨⎧==044v u 六、求解相应的支反力（运用静力学的平衡方程进行求解）3号节点和4号节点的支反力如下图所示：。
有限元计算分析报告***********院&……*设计专业2008年六月目录试题一 (1)试题二 (4)试题三 (6)试题四 (8)试题一问题描述：图示为一带圆孔的单位厚度的正方形平板，在x 方向作用均布压力0.25Mpa ，试用三节点常应变单元和六节点三角形单元对平板进行有限元分析。
图1.（1）分析：（1）该平板的受力属于平面应力问题，在对称外载荷的作用下，模型的受力也是对称的。
故而对该问题也可以简化为4/1平板倒圆角的模型。
在对称均布压力作用下，平板的内部应力图也应该是对称的，平板受到沿x 方向的正应力x σ和y 方向的应力y σ，板面上没有力的作用，即：0=zy τ，0=zx τ。
垂直于板面没有正应力作用，z σ=0。
即该数学模型与z 无关，仅仅是x 、y 的函数。
（2）简化为41平板计算时，可以将倒圆左边界视作x 约束而y 方向无约束，下边界视作y 方向有约束而x 方向没有约束。
分别用三节点和六节点画出单元网格。
（3） 由于平板的中间孔存在集中应力，所以在孔的附近的有限元网格需要细化，而远离孔的网格就可以不画那么细了。
有限元网格划分结果如下图1.（2）所示；数学建模：按照1/4计算，取点（0，0，0），由矢量（0.024，0.024，0）通过surface 创建平面，再建立2DArcangle 画出90度圆弧，将平面打断删去多余部分便得到了几何模型。
六节点应力图六节点应变图三节点应力图三节点应变图三节点不同数目网格应力图三节点不同数目网格数目应变图试题二问题描述：图示 2.（1）为带方孔（边长为80mm）的悬臂梁，其上受部分均布载荷（p=10Kn/m）作用，试采用一种平面单元，对图示两种结构进行有限元分析，并就方孔的布置进行分析比较，如将方孔设计为圆孔，结果有何变化？（板厚为1mm，材料为钢）分析：（1）该题同第一个问题一样是属于平面问题，是平面的应变问题。
（2）平面及面内无z方向应力分量，且限制了z向位移。
有限元计算结构力学fortran程序计算结构力学程序计算结构力学编程大作业时间,2007年6月!!!***************************************************************** ***********!!! 关于程序的说明!!!***************************************************************** ***********!一、功能:! 1、可计算包括节点力，一般非节点力，支座沉降、温度荷载作用、制造误差的平! 面桁架、梁、刚架及其组合结构的节点位移与杆端力;! 2、可同时计算多种工况下的节点位移与杆端力。
!******************************************************************* **********!******************************************************************* ***********!! 二、变量说明:! NE——单元数;! N——结构中自由度数;! NJ——节点数;! NS——特殊节点数，包括支座节点、主从节点(1节点不做主节点)、连接桁架的铰节点(没有转角);! NAI——结构的单元截面类型数;! MT——单元截面类型号;! NL——荷载工况数;! H——截面高度;! E——弹性模量;! JC——单元定位向量数组;! X(NJ),Y(NJ)——节点的X，Y坐标值;! JE(NE,2)——单元两端节点码数组;! AI(NAI,2)——按单元类型顺序存放A与I，AI(I，1)—第I类单元的截面积，AI(I，2)—第I类单元的! 惯性矩;! MT(NE)——单元所属单元类型号;! JS(NS,4)——特殊节点信息，JS(I,1)—结点码;JS(I,2)，JS(I,3)，JS(I,4)—U，V，CETA约束信息，! 有约束为1，没有约束为0;从节点某位移同主节点时位移时，该位移约束信息填主节点码;!! PJ(NP,3)——节点荷载信息数组;PJ(I，1)—节点力所在节点号;PJ(I，2)—节点力作用坐标方向:! 坐标方向U，V，M分别为1，2，3; PJ(I，3)—节点力的大小(含正负号);U，V方向集中力时，! 与坐标轴正向同向为正，M按右手法则为正;本程序推导过程取y轴向下为正。
有限元大作业魏博宇力学111 3111631031一、划分单元，确定半带宽。
x13 24 6 85 7 9 11 13 15 10 12 14 1618 20 22 2417 19 21 23 25 27 29 31 33 3526 28 30 32 34 36y单元数：36 ； 节点数 28。
1.单元局部节点编号： y x ijm ijm yx12 3 4 56 7 8 9 10 11121314 1516 17 18 19 20 212223 24 25 26 27 28、单元号 1 2 3 4 5 6 7 8 9 10 11 12i 2 4 5 5 7 8 8 9 9 11 12 12j 3 5 3 6 8 5 9 6 10 12 8 13 m 1 2 2 3 4 4 5 5 6 7 7 8 单元号13 14 15 16 17 18 19 20 21 22 23 24i 13 13 14 14 16 17 17 18 18 19 19 20j 9 14 10 15 17 12 18 13 19 14 20 15 m 8 9 9 10 11 11 12 12 13 13 14 14 单元号25 26 27 28 29 30 31 32 33 34 35 36i 20 22 23 23 24 24 25 25 26 26 27 27j 21 23 17 24 18 25 19 26 20 27 21 28 m 15 16 16 17 17 18 18 19 19 20 20 212.节点坐标结点号 1 2 3 4 5 6 7 8 9 10 11 12 13 14 x 0 0 1 0 1 2 0 1 2 3 0 1 2 3 y 0 1.5 1.5 3 3 3 4.5 4.5 4.5 4.5 6 6 6 6 结点号15 16 17 18 19 20 21 22 23 24 25 26 27 28 x 4 0 1 2 3 4 5 0 1 2 3 4 5 6 y 6 7.5 7.5 7.5 7.5 7.5 7.5 9 9 9 9 9 9 93.带状性1 5 10 15 20 25 28可求半带宽D=（7+1）× 2 = 16二、载荷向节点移置。
SHANGHAI JIAO TONG UNIVERSITY 平面应力问题解的Matlab实现姓名: heiya168学号: 帆哥班级：指导老师：目录1绪论 (1)2平面问题的四节点四边形单元 (2)2.1单元的构造 (2)2.2等参变换 (5)2.3边界条件的处理——置“1”法 (8)3有限元分析流程 (10)3.1程序原理和流程 (10)3.2使用的函数 (11)3.3文件管理 (11)3.4数据文件格式 (11)4算例——开方孔的矩形板拉伸分析 (13)4.1问题的具体参数与载荷 (13)4.2Matlab程序计算 (13)4.3ANSYS建模计算 (15)4.4误差分析 (17)5总结 (18)参考文献 (18)附录 (19)1绪论有限元方法(finite element method)，是求取复杂微分方程近似解的一种非常有效的工具，是现代数字化科技的一种重要基础性原理。
将它用于在科学研究中，可成为探究物质客观规律的先进手段。
将它应用于工程技术中，可成为工程设计和分析的可靠工具。
弹性体在载荷作用下，其基本方程可写成以下的三类方程和两种边界条件。
平衡方程——应力与外载荷的关系；几何方程——应变位移关系；物理方程——应力应变关系；力的边界条件；几何边界条件。
应用最小位能原理，并利用上述关系，最终建立由刚度方程，节点位移和等效节点载荷所构成的求解方程。
带入边界条件求解方程，就可以得出弹性力学问题的一般性解答。
本次大作业基于有限元方法的基本原理，使用Matlab这一平台，针对平面应力问题，采用四节点四边形单元编写了求解单元节点位移的程序。
主要内容包括：1）介绍有限元的基本原理；2）编程基本思路及流程介绍；3）程序原理及说明； 4）具体算例这四个部分。
2 平面问题的四节点四边形单元2.1 单元的构造（1） 单元的几何和节点描述平面4节点矩形单元如图4-6所示，单元的节点位移共有8个自由度。
节点的编号为1、2、3、4，各自的位置坐标为(x i, y i), i=1,2,3,4，各个节点的位移(分别沿x 方向和y 方向)为(u i,v i), i=1,2,3,4。
图2.1 平面4节点矩形单元若采用无量纲坐标, =x ya bξη=（2.1） 则单元4个节点的几何位置为112233441 =11 =11 =11 =-1ξηξηξηξη==-=--=（2.2）将所有节点上的位移组成一个列阵，记作e q ；同样，将所有节点上的各个力也组成一个列阵，记作e P ，那么11223344(81)11223344(81)[ v v u v u v ][ ]e TeTx y x y x y x y q u u P P P P P P P P P ⨯⨯== （2.3）若该单元承受分布外载，可以将其等效到节点上，也可以表示为如式(2.3)所示的节点力。
利用函数插值、几何方程、物理方程以及势能计算公式，可以将单元的所有力学参量用节点位移列阵及相关的插值函数e q 来表示；下面进行具体的推导。
（2） 单元位移场的表达从图2-1可以看出，节点条件共有8个，即x 方向4个1234(,,,)u u u u ，y 方向4个1234(,,,)v v v v ，因此，x 和y 方向的位移场可以各有4个待定系数，即取以下多项式作为单元的位移场模式01230123(,)(,)u x y a a x a y a xy v x y b b x b y b xy=+++=+++ （2.4）它们是具有完全一次项的非完全二次项，以上两式中右端的第四项是考虑到x 方向和y 方向的对称性而取的，除此外xy 项还有个重要特点，就是“双线性”，当x 或y 不变时，沿y 或x 方向位移函数呈线性变化，这与前面的线性项最为相容，而2x 或2y 项是二次曲线变化的。
因此，未选2x 或2y 项。
由节点条件，在x =x i ，y =y i 处，有(,)(,)i i i i i iu x y u v x y v == 1,2,3,4i = （2.5）将式(2.5)代入式(2.4)中，可以求解出待定系数a 0，…，a 3和b 0，…，b 3，然后代回式(4-52)中，经整理后有1122334411223344(,)(,)(,)(,)(,)(,)(,)(,)(,)(,)u x y N x y u N x y u N x y u N x y u v x y N x y v N x y v N x y v N x y v =+++=+++ （2.6）其中12341(,)(1)(1)41(,)(1)(1)41(,)(1)(1)41(,)(1)(1)4x yN x y a b x yN x y a bx yN x y a b x yN x y a b=++=-+=--=+- （2.7）如以无量纲坐标系(2.1)来表达，式(2.7)可以写成1(1)(1)4i i i N ξξηη=++1,2,3,4i = （2.8） 将式(2.6)写成矩阵形式，有11212342(21)(28)(81)12343344(,)u (,)=(,)eu v u N N N N v u x y x y N qN N N N u v x y v u v ⨯⨯⨯⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎡⎤⎡⎤⎢⎥==⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦（2.9） 其中，(,)N x y 为该单元的形状函数矩阵。
（3） 单元应变场的表达由弹性力学平面问题的几何方程(矩阵形式)，有单元应变的表达[][](38)(31)(21)(28)(81)(81)(32)(32)(,)u xx e e yy xy x y N q B q εεεγ⨯⨯⨯⨯⨯⨯⨯⨯⎡⎤⎢⎥==∂=∂⋅=⋅⎢⎥⎢⎥⎣⎦（2.10）其中几何矩阵B(x,y)为[]1234(38)(28)(32)1234(,)x N N N N B x y N N N N N y yx ⨯⨯⨯⎡⎤∂⎢⎥∂⎢⎥⎢⎥⎡⎤∂=∂=⎢⎥⎢⎥∂⎣⎦⎢⎥⎢⎥∂∂⎢⎥∂∂⎣⎦=1234(32)(32)(32)(32)B B B B ⨯⨯⨯⨯⎡⎤⎢⎥⎣⎦（2.11） 式(4-59)中的子矩阵i B 为(32)1,2,3,4i i i i i N x N B i y N N yx ⨯⎡⎤∂⎢⎥∂⎢⎥⎢⎥∂==⎢⎥∂⎢⎥⎢⎥∂∂⎢⎥∂∂⎣⎦（2.12）（4） 单元应力场的表达由弹性力学中平面问题的物理方程，可得到单元的应力表达式(33)(33)(38)(31)(31)(38)(81)(81)e e D D B q S q σε⨯⨯⨯⨯⨯⨯⨯⨯=== （2.13）（5） 单元应力场的表达以上已将单元的三大基本变量(,,)u εσ用基于节点位移列阵来进行表达，见式(2.9)、式(2.10)及式(2.13)；将其代入单元的势能表达式中，有12e eT e eeT e q K q P q ∏=- （2.14） 其中e K 是4节点矩形单元的刚度矩阵。
将单元的势能对节点位移e q 取一阶极值，可得到单元的刚度方程(88)(81)(81)e e e K q P ⨯⨯⨯⋅= （2.15）2.2 等参变换由前面的单元构造过程可以看出，一个单元的关键就是计算它的刚度矩阵，而由刚度矩阵的构成可知要实现两个坐标系中单元刚度矩阵的变换，必须计算两个坐标系之间的三种映射关系：坐标映射(,)(,)x y ξη⇒ （2.16）偏导数映射(,)(,)x y ξη∂∂∂∂⇒∂∂∂∂ （2.17） 面积映射ee A A dxdy d d ξη⇒⎰⎰ （2.18）图2.2矩形单元映射为任意四边形单元（1） 两个坐标系之间的函数映射设如图4-17所示的两个坐标系的坐标映射关系为(,)(,)x x y y ξηξη==（2.19）x 和y 方向上可以分别写出各包含有4个待定系数的多项式，即01230123(,)(,)x a a a a y b b b b ξηξηξηξηξηξη=+++=+++ （2.20）其中待定系数a 0，…，a 3和b 0，…，b 3可由节点映射条件(2.19)来唯一确定。
对照前面4节点矩形单元的单元位移函数式(2.5)，映射函数式(2.20)具有完全相同的形式，同样，将求出的待定系数再代回式(2.20)中，重写该式为1122334411223344(,)(,)(,)(,)(,)(,)(,)(,)(,)(,)x N x N x N x N x y N y N y N y N y ξηξηξηξηξηξηξηξηξηξη=+++=+++ （2.21）其中1(1)(1) i=1,2,3,44i i i N ξξηη=++ （2.22） 11212342(21)(28)(81)12343344(,)x (,)=(,)(,)e x y x N N N N y x N q N N N N x y y x y ξηξηξηξη⨯⨯⨯⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎡⎤⎡⎤⎢⎥==⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦（2.23） 这就可以实现两个坐标系间的映射。
（2）两个坐标系之间的偏导数映射对物理坐标系(x,y)中的任意一个函数Φ(x,y)，求它的偏导数，有x yx y x yx y ξξξηηη∂Φ∂Φ∂∂Φ∂=+∂∂∂∂∂∂Φ∂Φ∂∂Φ∂=+∂∂∂∂∂ （2.24）则偏导数的变换关系为x y x y x y x y ξξξηηη∂∂∂∂∂=+∂∂∂∂∂∂∂∂∂∂=+∂∂∂∂∂ （2.25）写成矩阵形式，有x J y ξη∂⎡⎤∂⎡⎤⎢⎥⎢⎥∂∂⎢⎥⎢⎥=∂∂⎢⎥⎢⎥⎢⎥⎢⎥∂∂⎣⎦⎣⎦（2.26） 其中x y J x y ξξηη∂∂⎡⎤⎢⎥∂∂⎢⎥=∂∂⎢⎥⎢⎥∂∂⎣⎦（2.27） 两个坐标系的偏导数映射关系11y y x J x x y J ηξξηηξξη⎛⎫∂∂∂∂∂=- ⎪∂∂∂∂∂⎝⎭⎛⎫∂∂∂∂∂=-+ ⎪∂∂∂∂∂⎝⎭（2.28）（3）两个坐标系之间的面积元映射如图2-2所示，在物理坐标系(x , y )中，由d ξ和d η所围成的微小平行四边形，其面积为dA d d ξη=⨯ （2.29）由于d ξ和d η在物理坐标系(x , y )中的分量为x yd d i d j x y d d i d jξξξξξηηηηη∂∂=⋅+⋅∂∂∂∂=⋅+⋅∂∂ （2.30）其中i 和j 分别为物理坐标系(x , y )中的x 方向和y 方向的单位向量。
由(2.29)式，则有面积微元的变换计算xy d d dA J d d x y d d ξξξξξηηηηη∂∂∂∂==∂∂∂∂ （2.31）2.3 边界条件的处理——置“1”法设边界条件为第r 个自由度的位移为零，即0r q =；可置整体刚度矩阵中所对应对角元素位置的1rr k =，而该行和该列的其它元素为零，即0()rs sr k k r s ==≠，同时也置对应的载荷元素 则0r p =，则（2.32）进行以上设置后，这时方程(2.32)应等价于原方程加上了边界条件0r q =，下面考察这种等价性，就(2.32)中的第r 行，有rr r r k q p ⋅= （2.33）由于置1,0rr r k p ==，则有0r q =即为所要求的位移边界条件。
