清华大学-弹性力学有限元大作业
弹性力学及有限元习题参考答案(赵均海、汪梦甫)汇编
(1,1)
4
(0,0)
5
(1,1)
(0,0)
(1,0)
号
③
2
(1,1)
5
(1,0)
6
(2,0)
④
2
(1,1)
3
(2,1)
6
(2,0)
元
②
1
2
1
2
1
2
cm
1
-1
-1
0
0
-1
1
0
-1
K11
K K 21
K31
1
2
K12
K 22
K32
0
1
-1
1
-1
0
1
0
-1
0
-1
1
K13
K 23 (i 1, j 2, m 3)
（3）主方向：
l（ − ) + + = 0
+ ( − ) + = 0
+ + ( − ) = 0
2 + 2 + 2 = 1
第一主方向：将1 = −46 MPa 及个分量代入上式，有：
101l + 40 = 0
0
0
0
0
5 2
7 17
0
0
0
0
5 12
0
0
17
0 12 2
0
0
0
0
0
17
5
5
0
0
0
0 12 5
34
0 12 5
0
0
弹性力学与有限元分析试题及其答案
一、填空题1、弹性力学研究弹性体由于受外力作用、边界约束或温度改变等原因而发生的应力、形变和位移。
2、在弹性力学中规定，线应变以伸长时为正，缩短时为负，与正应力的正负号规定相适应。
3、在弹性力学中规定，切应变以直角变小时为正，变大时为负，与切应力的正负号规定相适应。
4、物体受外力以后，其内部将发生内力，它的集度称为应力。
与物体的形变和材料强度直接有关的，是应力在其作用截面的法线方向和切线方向的分量，也就是正应力和切应力。
应力及其分量的量纲是L -1MT -2。
5、弹性力学的基本假定为连续性、完全弹性、均匀性、各向同性。
6、平面问题分为平面应力问题和平面应变问题。
7、已知一点处的应力分量100=x σMPa ，50=y σMPa ，5010=xy τ MPa ，则主应力=1σ150MPa ，=2σ0MPa ，=1α6135' 。
8、已知一点处的应力分量， 200=x σMPa ，0=y σMPa ，400-=xy τ MPa ，则主应力=1σ512 MPa ，=2σ-312 MPa ，=1α-37°57′。
9、已知一点处的应力分量，2000-=x σMPa ，1000=y σMPa ，400-=xy τ MPa ，则主应力=1σ1052MPa ，=2σ-2052 MPa ，=1α-82°32′。
10、在弹性力学里分析问题，要考虑静力学、几何学和物理学三方面条件，分别建立三套方程。
11、表示应力分量与体力分量之间关系的方程为平衡微分方程。
12、边界条件表示边界上位移与约束，或应力与面力之间的关系式。
分为位移边界条件、应力边界条件和混合边界条件。
13、按应力求解平面问题时常采用逆解法和半逆解法。
14、有限单元法首先将连续体变换成为离散化结构，然后再用结构力学位移法进行求解。
其具体步骤分为单元分析和整体分析两部分。
15、每个单元的位移一般总是包含着两部分：一部分是由本单元的形变引起的，另一部分是由于其他单元发生了形变而连带引起的。
清华大学弹性力学有限元大作业
弹性力学有限元大作业一、模型信息：已知：材料为铝合金。
E=71GPa ，v=0.3.矩形平板的几何参数：板长为480mm ，宽为360mm ，厚度为2mm ；图形如下图；加肋平板：二、matlab 编程实现1、程序相关说明：计算使用的软件为：matlab2010a 主函数：main.m 主要计算部分子函数：Grids.m 生成网格，节点数为：+1*+1I J （）（）、单元数： 2**I J AssembleK.m 将单元刚度矩阵组装成总刚度矩阵（叠加方法）GenerateB.m 生成单元格e B 矩阵 GenerateS.m 生成单元格e S 矩阵 GenerateK.m 生成单元刚度矩阵2、网格划分：利用Grid.m 子函数，取2020I J ==、，即可以得到网格如下： 节点数为：441个，单元格数：800个3、计算过程及结果 （1）、网格划分：通过Grid.m ，生成节点数为：441个、单元格数：800个的网格 （2）、生成总刚度矩阵K ：通过GenerateK.m 、AssembleK.m 生成总刚度矩阵 采用常应变三角单元，e e u N a =,易得=e e B LN由平面应力问题，可以确定2101011002E D νννν⎡⎤⎢⎥⎢⎥=⎢⎥-⎢⎥-⎢⎥⎣⎦即e e S DB =单元刚度矩阵为：e eT e K AtB DB = 总刚度矩阵为：eTe e eK GK G =∑（3）、求解过程：系统平衡方程为：Ka P = 将方程进一步划分为：E EF E E E T F F EFF K K d f r d f K K +⎡⎤⎡⎤⎡⎤=⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦ 通过已知边界条件（位移、载荷），确定E E F d f f 、、 ，从而将K 矩阵划分为四个模块：E EF TEF F K K K K ⎡⎤⎢⎥⎣⎦1()E E E EF F E TF F F EF E r K d K d f d K f K d -=+-=-支反力：部分位移：即整体位移向量为：E F d a d ⎡⎤=⎢⎥⎣⎦整体力边界条件为：E E F f r P f +⎡⎤=⎢⎥⎣⎦（4）后处理：（应力、应变、抹平） a 、单元应力、应变：e e e ee eS a B aσε==b 、抹平得到节点应力、应变：将每个节点参与组成的单元应力、应变叠加，然后除以叠加的单元数，得到抹平后的节点应力、应变。
弹性力学与有限元法习题集
' yx dx
0
' yx
dM dx
S
* 2
I
Q(x) I
n 2 y
y1
dy1
Q(x) I
(n2 4
y2 ) 1 Q(x) (n2 2 2I 4
y2)
2019/7/29
slide14
返回
由剪应力互等定理，
yx
' yx
Q(x) (n2 2I 4
答案 返回
1. 有限单元法的含义？ 答：用有限个单元将连续体离散化，通过对有限个单元作分片插
值求解各种力学、物理问题的一种数值方法。连续体的单元是 各种形状（ 如三角形、四边形、六面体等 ）的单元体。
2.有限元法的解题思路？ 答（1）网格划分； （2）单元分析；（3）整体分析。
3.有限元法的优点？ 答（1）物理概念清晰，便于入门；
13. 已知某单元，其结点编号为i，j，m，其坐标依次为（2， 2）、（6，3）、（5，6），试写出其形函数Ni，Nj，Nm 及单 元的应变矩阵。
2019/7/29
slide22
答案 返回
14. 图示平面应力状态的直角三角形单元及其结点编码，设
1 6
试求：
（1）形态矩阵[N]； （2）几何矩阵[B]及应力转移矩阵[S]; （3）单元刚度矩阵[k]e
6. 应用几何方程推导应变分量应满足下列变形协调方程。
2 x 2 y 2 xy
y2 x2 xy
2019/7/29
slide6
答案 返回
7. 悬臂梁在三角形分布载荷作用下，可以看成平面应力问题，
应力分量表达式为， x
弹性力学有限元习题答案
弹性力学有限元习题答案弹性力学有限元习题答案弹性力学是研究物体在受力作用下产生的变形和应力分布的学科。
有限元方法是一种数值计算方法，用于求解复杂的力学问题。
在弹性力学有限元习题中，我们需要运用弹性力学理论和有限元方法来解答问题。
下面，将给出一些常见的弹性力学有限元习题的解答。
1. 问题描述：一根长为L的均匀梁，两端固定支承，受到均匀分布载荷q作用。
求梁的挠度分布和最大挠度。
解答：首先，我们可以根据弹性力学理论得到梁的挠度方程。
然后，将梁分割为若干个小段，利用有限元方法近似求解挠度分布。
最后，通过计算得到的挠度分布，可以找到最大挠度的位置和数值。
2. 问题描述：一个矩形薄板，边长为a和b，厚度为t。
板的一侧边固定支承，另一侧边受到均匀分布载荷q作用。
求板的应力分布和最大应力。
解答：根据弹性力学理论，可以得到薄板的应力分布方程。
然后，将薄板分割为若干个小单元，利用有限元方法近似求解应力分布。
最后，通过计算得到的应力分布，可以找到最大应力的位置和数值。
3. 问题描述：一个长方体结构，由若干个杆件和节点组成。
杆件的长度、截面积和杨氏模量已知。
节点上的载荷和位移边界条件已知。
求结构的应力分布和变形。
解答：首先，我们可以根据弹性力学理论得到结构的应力分布方程和变形方程。
然后，将结构分割为若干个小单元，利用有限元方法近似求解应力分布和变形。
最后，通过计算得到的应力分布和变形，可以分析结构的受力情况和变形情况。
以上是一些常见的弹性力学有限元习题的解答方法。
在实际应用中，弹性力学有限元方法可以用于求解各种复杂的力学问题，如梁、板、壳体、结构等。
通过运用弹性力学理论和有限元方法，可以得到准确的应力分布和变形情况，为工程设计和分析提供有力的支持。
总结起来，弹性力学有限元习题的解答需要运用弹性力学理论和有限元方法，通过建立适当的数学模型和边界条件，求解应力分布和变形情况。
这些解答方法在实际工程中具有广泛的应用价值，可以帮助工程师和科研人员分析和解决各种力学问题。
清华大学弹性力学-大作业matlab代码
f1=int(N',m,-1,1); f=1/(4*a*b)*2*int(f1,n,-1,1); f=subs(f); %重积分过程求单元结点荷载 f %初始给定 K 为零阵
K=zeros((a+1)*(b+1),(a+1)*(b+1)); for x1=1:(1+a) for y1=1:(1+b) for x2=1:(1+a) for y2=1:(1+b) s1=(y1-1)*(a+1)+x1; s2=(y2-1)*(a+1)+x2; if x1==x2 & y1==y2 if x1==1
T=round(T); K=K(T,T); F=F(T); U=K\F; U(1,:) %A 点位移
end if abs(x1-x2)==1 & abs(y1-y2)==1 if (x1-x2)*(y1-y2)==1 K(s1,s2)=k(2,4); else K(s1,s2)=k(1,3); end end end end end end K; %以上循环生成 K 矩阵 %如果需要可以通过此处查看 K
1/3]; A\b 3.a×b 个矩形单元 %a×b 个矩形单元 clc,clear all a=input('单元格长数'); a=round(a); b=input('单元格宽数'); b=round(b); x=0:1/a:1; y=0:1/b:1; A=[1 -1 1 -1 1111 1 1 -1 -1 1 -1 -1 1]; P=inv(A); syms m n; %m,n 相当于ξ和η N=[1,m,n,m*n]*P; B=[2*a*(P(2,:)+n*P(4,:)) 2*b*(P(3,:)+m*P(4,:))]; D=eye(2); k1=int(B'*D*B,m,-1,1); k=1/(4*a*b)*int(k1,n,-1,1); k=subs(k); %重积分过程求单元刚度矩阵 k
清华大学弹性力学往届考题答案
弹性力学与有限元期末复习资料（清华大学段云岭教授班）1. 弹性力学的任务与研究对象：弹性力学的基本任务就是从宏观上研究弹性体对外因（主要是外力，温度等）作用的反应，确切的说，就是确定此时的应力和应变以及与应变有直接关系的位移。
2. 弹性力学的的基本假设：a. 物体构造的连续性假设；b. 物体的完全弹性假设；c. 物体的均匀性假设；d. 物体的各向同性假设；e. 小变形假设。
Ps.前四条是满足弹性体的假设，第五条是在我们利用精典弹性力学时用到的假设。
3. 简述两种弹性力学的平面问题并举例说明：a ． 平面应力问题。
几何形状特征：物体在一个坐标方向的几何尺寸远远小于其他两个坐标方向的几何尺寸，例如薄板；荷载特征：在薄板的两个侧表面上无表面荷载，作用于边缘的表面力平行于板面，且沿厚度不发生变化，或虽沿厚度变化但对称于板的中间平面，体积力亦平行于板面且沿厚度不变。
简化分析：其中某些量为0，与z 有关的应力都为0.只剩下三个未知量，,,x y xy σστ，并且可以假设这三个应沿板厚度方向是均匀分布的。
实际例子：工程中的梁桥，链条的平面链环，被圆孔或者圆槽削弱的薄板等，薄板等等b ． 平面应变问题。
几何形状特征：荷载特征：简化分析：实际例子：水坝， 辊轴支座的柱形辊轴以及隧道等4. 弹性力学研究问题，需要从哪几方面考虑，并说明将得到什么结果：a. 在物体内部任取一微元体都必须是平衡的，这将导出平衡微分方程；b. 在物体边界上任取一微元体也都必须是平衡的，这将导出力的边界条件。
5. 以位移作为边界未知量，论述弹性力学求解思路和过程：基本方程和边界条件如下：在研究对象的内部（Ω）满足[]{}{}T B F σ=- （1）{}[]{}B f ε= （2）{}[]{}D σε= （3）在研究对象边界处满足[]{}{}L F ∑= 在边界S σ上 （4）其中[]∑=x yx xy y σττσ⎛⎫ ⎪⎝⎭{}{}f f = 在边界u S 上（5）{}{}f f +-= 在内部的孔洞处（6）将（2）（3）代入（1），可以得到方程[][][]{}{}T B D B f F =- （7）如果已知以位移作为边界条件的话，那么关键在于处理方程（5）上。
弹性力学有限元考试卷与答案(AB卷)
2009-2010学年第一学期《弹性力学有限元》课内考试A卷授课班号年级专业学号姓名一、判断正误（×）1. 节点的位置依赖于形态，而并不依赖于载荷的位置（√）2. 对于高压电线的铁塔那样的框架结构的模型化处理使用梁单元（×）3. 不能把梁单元、壳单元和实体单元混合在一起作成模型（√）4. 四边形的平面单元尽可能作成接近正方形形状的单元（×）5. 平面应变单元也好，平面应力单元也好，如果以单位厚来作模型化处理的话会得到一样的答案（×）6. 用有限元法不可以对运动的物体的结构进行静力分析（√）7. 一般应力变化大的地方单元尺寸要划的小才好（×）8. 所谓全约束只要将位移自由度约束住，而不必约束转动自由度（×）9. 线性应力分析也可以得到极大的变形（√）10. 同一载荷作用下的结构，所给材料的弹性模量越大则变形值越小二、填空1．平面应力问题与薄板弯曲问题的弹性体几何形状都是薄板，但前者受力特点是：平行于板面且沿厚度均布载荷作用，变形发生在板面内；后者受力特点是：垂直于板面的力的作用，板将变成有弯有扭的曲面。
（3分）2．平面应力问题与平面应变问题都具有三个独立的应力分量：σx，σy，τxy，三个独立的应变分量：εx，εy，γxy，但对应的弹性体几何形状前者为薄板，后者为长柱体。
（3分）3．位移模式需反映刚体位移 ，反映 常变形 ，满足 单元边界上位移连续 。
（3分）4．单元刚度矩阵的特点有：对称性 ， 奇异性 ，还可按节点分块。
（2分） 5．薄板弯曲问题每个节点有个3自由度，分别是：w 、θx 、θy ，但其中只有 一个是独立的，其余两个可以用它表示为：,x y w wy xθθ∂∂==-∂∂。
（3分） 6．用有限元程序计算分析一结构的强度须提供（4分） ① 几何信息：节点坐标，单元节点组成，板厚度，梁截面等 ② 材料信息：弹性模量，泊松比，密度等 ③ 约束信息：固定约束，对称约束等④ 载荷信息：集中力，集中力矩，分布面力，分布体力等7．轴对称问题单元形状为：三角形或四边形截面的空间环形单元 ，由于轴对称的特性，任意一点变形只发生在子午面上，因此可以作为 二 维问题处理。
清华大学弹性力学-大作业
(1b)
其中 2 为 Laplace 算子， (1,1) (1,1) 为求解区域。 y 1
y 1
u0
-1
O -1
1
x
u 0 n
1
u0
O
u 0 n
图 2 第一象限
x 1
图 1 薄膜小挠度弯曲模型
其势能泛函为：
I (u ) 1 u u [( ) 2 ( ) 2 ]dA f u dA 2 x y x 1 ， y 1 u 0，
1
y 2
y 1
2
1
1
3x 4 (a)
5 (b)
x 3
4 (c)
3
x
图 3 有限元单元种类
(2)
s.t.
(1) 根据最小势能原理，试由式(2)导出式(1)。 (2) 当 f ( x, y ) 2 ,即薄膜受均布荷载作用时，其挠度 u 关于 x 轴和 y 轴均对称， 取第一象限 1 (0,1) (0,1) 为求解区域，如图 2 所示（图中 n 为外法线方向） ， 请分别使用如下三种有限元单元计算原点 O 处的挠度 u (0,0) ， 并比较分析有限 元解与精确解的误差： (a) 2 个常应变三角形单元（如图 3a） ； (b) 4 个常应变三角形单元（如图 3b） ； (c) 1 个四结点四边形单元（如图 3c） 。 （注：精确解为 u (0,0) 0.5893 。 ） y 1 2 1 4
有限元法大作业——薄膜小挠度弯曲问题的有限元求解
如图 1 所示，边长为 2、四边固支的正方形薄膜，受到横向分布荷载 f ( x, y ) 的作用，用挠度 u u ( x, y ) 表示的平衡微分方程为 (1a) 2u f ( x, y ) ， ( x, y )
弹性力学大作业报告
弹性力学文献阅读报告水坝的应力计算问题及其讨论文章题目：Stress Analysis of the Rectangular Retaining Wall Under Water Pressure Based on Elasticity Mechanics作者：Yongwei Wang Yanqin Guo Huanghe Science and Technology College Institute of Technology, Henan Province, china 一、对文章的理解及讨论在工程实践中，我们经常需要了解坝体的应力分布情况，来确定坝体的危险截面，给我们的工程设计和应用带来参考。
这篇文章中将水坝简化为矩形坝体，并用弹性力学的方法分析坝体的应力和位移。
水坝在工程实际中主要受到两种力，即自身体力和液体的压力。
这里取水坝的一个横截面进行分析，将其视为平面应变问题，示意图如图1所示：水坝力学模型1.水坝力学模型图1.这里假设坝体的密度为ρ1，左侧的液体密度为ρ2。
文章中采用的是半逆解法，由应力的表达式假设出应力函数求解，过程如以下公式所示：所以可以假设出应力函数的表达式为：这里文章中并没有解释σx假设成这样的原因，我认为这里之所以将σx假设成这样是因为坝体内的应力由两部分影响组成，即自重ρ1gy和水压ρ2gy,这两部分均为关于y的正比例函数，故σx假设成如上的形式。
将应力函数代入相容方程得：因为上式要对无数y值成立，所以各次项系数都要为0，即：积分可得到下式：将其代入应力函数的表达式,,可得：将其代入应力函数的表达式利用应力函数求出应力分量为：分析主要边界的应力边界条件可得：运用圣维南原理列出次要边界的应力边界条件：这里原文中运用了圣维南原理对上表面的σy 积分为0，其实是对应力边界的一种放松要求，不过由验算得知，这里放松之后的结果和严格遵循边界条件的计算结果相同。
由以上的边界可以解得各应力分量为：再由物理方程求水平位移分量：文章的内容到这里就结束了，在介绍中提到求解应力分量是为了确定危险截面，但是文章正文内容并没有提及。
弹性力学及有限元试题
弹性力学及有限元试题(一) 问答题（20分）1、什么是圣维南原理？举例说明怎样把它应用于工程问题的简化中。
2、什么叫做一点的应力状态？如何表示一点的应力状态（要求具体说明或表达）。
3、何谓逆解法和半逆解法？它们的理论依据是什么？4、什么是平面应力问题？什么是平面应变问题？分别写出弹性力学平面应力问题和平面应变问题的物理方程。
5、要保证有限元方法解答的收敛性，位移模式必须满足那些条件？(二) (10分)1.利用坐标变换从直角坐标的平衡方程推导极坐标下平衡方程（无体力）。
2.利用坐标变换从直角坐标下几何方程推导极坐标下几何方程。
（三）已知，其他应力分量为零，求位移场。
（10分）（四）设有矩形截面的悬臂粱，在自由端受有集中荷载F；体力可以不计。
试根据材料力学公式，写出弯应力σx和切应力τxy的表达式，并取挤压应力σy=0，然后证明，这些表达式满足平衡微分方程和相容方程，再说明，这些表达式是否就表示正确的解答（10分）。
（五）设半平面体在直边界上受有集中力偶，单位宽度上力偶矩为M，试求应力分量（10分）。
提示：单位厚度上的力偶矩M的量纲是LMT-2，应力只能是M/ρ2的形式，所以可假设应力函数由：Φ=Φ(φ).(六) 铅直平面内的正方形薄板，边长为2a，四边固定，图5—18，只受重力的作用。
设μ=0，试取位移分量的表达式为用瑞利—里茨法求解（15分）。
（七）试按图示网格求解结点位移，取t =1m，μ= 0（15分）。
（八）用刚度集成法求下图所示结构的整体刚度矩阵K。
（10分）要求：单元刚度矩阵元素用ek形式表示；单元刚度矩阵用e K形式表ij示，其中e为单元号。
清华大学弹性有限元2006年期末试题
2La弹性力学基础及有限元分析期末试题（A）姓名：________________ 学号：_________________时间： 2007年1月7日 1.问答题（50分）1)根据具体问题边界条件类型的不同，弹性力学问题的求解时常遇到的边值问题有几类？它们分别是什么？2)写出如图1－2所示问题的边界条件。
3)利用有限元位移法求解弹性问题的基本原理是什么？写出相应的平衡方程和位移边界条件的等效积分弱形式的具体表达式。
4)已知如图1－4所示的坝体。
在利用有限元方法求解过程中，为减少计算量而在过渡段采用了变结点单元：i.写出7点和8点在如图所示自然坐标下的插值函数表达式；ii.写出求解时应给出的1点和3点的边界条件。
5)薄板的小变形分析是基于什么假设才成立的？2.计算及证明题（50分）：1)(15分)如图3-1和等值拉压载荷q的共同作用。
距离边缘较远处有一小孔，孔的半径为a。
求：i.板内应力；ii.孔边的最大和最小正应力，并指出最大正应力点的位置。
2)(20分)已知：矩形梁在如图3-2所示的载荷作用下。
已知mL5.0=，2610mNEI⋅=,KNP6=，12kN/mq=，26kN/mq=。
i.利用虚位移方法求2点的转角。
ii.如果利用有限元方法，给出最终求解的总体方程并计算2点在y方向的位移。
3)矩形薄板OABC如图所示，其OA边和OC边为简支边，AB边及BC边为自由边，在B点受横向集中载荷P作用。
试证：w axy=能满足一切条件，并求出挠度、内力及反力。
（15）xOq2Lbii.3.如图3-3所示的桁架结构。
试用最小余能原理求杆12试用有限元方法求点的位移。
AOC3-2第 2 页/共2页。
弹性力学与有限元分析试题及参考答案
弹性力学与有限元分析试题及参考答案四、分析计算题1、试写出无体力情况下平面问题地应力分量存在地必要条件，并考虑下列平面问题地应力分量是否可能在弹性体中存在.资料个人收集整理，勿做商业用途（1）By Ax x，Dy Cx y，Fy Ex xy；（2）)(22y xA x，)(22y xB y，Cxy xy；其中，A ，B ，C ，D ，E ，F 为常数.解：应力分量存在地必要条件是必须满足下列条件：（1）在区域内地平衡微分方程xyyxxyyyxx；（2）在区域内地相容方程02222yx yx；（3）在边界上地应力边界条件sflms f ml ysxy yxs yx x；（4）对于多连体地位移单值条件.资料个人收集整理，勿做商业用途（1）此组应力分量满足相容方程.为了满足平衡微分方程，必须A=-F ，D=-E.此外还应满足应力边界条件.资料个人收集整理，勿做商业用途（2）为了满足相容方程，其系数必须满足A+B=0；为了满足平衡微分方程，其系数必须满足A=B=-C/2.上两式是矛盾地，因此，此组应力分量不可能存在.资料个人收集整理，勿做商业用途2、已知应力分量312x C Qxyx，2223xy C y，y x C yC xy2332，体力不计，Q 为常数.试利用平衡微分方程求系数C1，C2，C3.解：将所给应力分量代入平衡微分方程0xyyxxyyyxx得23033322322212xy C xy C xC yC xC Qy即230333222231xy C C yC Q xC C 由x ，y 地任意性，得23030332231C C C Q C C 由此解得，61Q C ，32Q C ，23Q C 3、已知应力分量q x，q y，0xy，判断该应力分量是否满足平衡微分方程和相容方程.解：将已知应力分量q x，q y，0xy，代入平衡微分方程0Y xyX yxxyyyxx可知，已知应力分量q x，q y，0xy一般不满足平衡微分方程，只有体力忽略不计时才满足.按应力求解平面应力问题地相容方程：yx xyxyxy yx 22222)1(2)()(将已知应力分量q x，q y，0xy代入上式，可知满足相容方程.按应力求解平面应变问题地相容方程：yx xyxyxyyx2222212)1()1(将已知应力分量q x，q y，0xy代入上式，可知满足相容方程.4、试写出平面问题地应变分量存在地必要条件，并考虑下列平面问题地应变分量是否可能存在.（1）Axy x，3By y，2Dy C xy；（2）2Ay x ，y Bx y2，Cxy xy；（3）0x，0y ，Cxy xy ；其中，A ，B ，C ，D 为常数.解：应变分量存在地必要条件是满足形变协调条件，即yx x yxyyx 22222将以上应变分量代入上面地形变协调方程，可知：（1）相容.（2）C By A 22（1分）；这组应力分量若存在，则须满足：B=0，2A=C.（3）0=C ；这组应力分量若存在，则须满足：C=0，则0x，0y，0xy（1分）.5、证明应力函数2by 能满足相容方程，并考察在如图所示地矩形板和坐标系中能解决什么问题（体力不计，0b ）.解：将应力函数2by 代入相容方程24422444yyxx可知，所给应力函数2by 能满足相容方程.由于不计体力，对应地应力分量为b yx222，022xy，2yx xy对于图示地矩形板和坐标系，当板内发生上述应力时，根据边界条件，上下左右四个边上地面力分别为：上边，2h y，0l ，1m，0)(2h yxyxf ，0)(2h yyyf ；下边，2h y，0l ，1m ，0)(2h yxyx f ，0)(2h yyy f ；左边，2l x，1l ，0m ，b f l xxx2)(2，0)(2l xxyy f ；右边，2l x，1l ，0m ，b f l xxx 2)(2，0)(2l xxyy f .l/2l/2h/2h/2yxOOx b 可见，上下两边没有面力，而左右两边分别受有向左和向右地均布面力2b.因此，应力函数2by能解决矩形板在x 方向受均布拉力（b>0）和均布压力（b<0）地问题.资料个人收集整理，勿做商业用途6、证明应力函数axy 能满足相容方程，并考察在如图所示地矩形板和坐标系中能解决什么问题（体力不计，0a ）.解：将应力函数axy 代入相容方程24422444yyxx可知，所给应力函数axy 能满足相容方程.由于不计体力，对应地应力分量为022yx，022xy，ayx xy2对于图示地矩形板和坐标系，当板内发生上述应力时，根据边界条件，上下左右四个边上地面力分别为：上边，2h y，0l ，1m ，a f h yxyx2)(，0)(2h yyyf ；下边，2h y ，0l ，1m ，a f h yxyx 2)(，0)(2h yyy f ；左边，2l x ，1l ，0m ，0)(2l xxxf ，a f l xxyy 2)(；右边，2l x，1l ，0m ，0)(2l xxx f ，a f l xxyy 2)(.可见，在左右两边分别受有向下和向上地均布面力a ，而在上下两边分别受有向右和向左地均布面力 a.因此，应力函数axy 能解决矩形板受均布剪力地问题.资料个人收集整理，勿做商业用途7、如图所示地矩形截面地长坚柱，密度为，在一边侧面上受均布剪力，试求应力分量.解：根据结构地特点和受力情况，可以假定纵向纤维互不挤压，即设0x.由此可知l/2l/2h/2h/2yxO22yx将上式对y 积分两次，可得如下应力函数表达式)()(,21x f y x f yx 将上式代入应力函数所应满足地相容方程则可得)()(424414dxx f d dxx f d y2.3 直角三角形固定在刚性基础上，受齐顶地水压力和自重作用，如图 2.14所示.若按一个单元计算，水地容重g ，三角形平面构件容重g ,取泊松比v =1/6，试求顶点位移和固定面上地反力.资料个人收集整理，勿做商业用途解：按逆时针编码，局部编码与整体编码相同：1-2-3建立坐标)0,0(3)3,0(20,21:a a xoy （1）求形函数矩阵：aa a a 60321a b b a b 303321ac a c c 220321图（2.14）形函数:)(21y c x b a AN i i i i233221aa a A所以：ay ax Na y N a xN 32132321形函数地矩阵为：ay a xay ax ay a x a y a xN NN Nmji321302003210302（2）刚度矩阵333231232221131211KKKK K K K K KKesr sr s r sr s r s r s r sr rsb bc c c b b c b c c b c c b b AEtK21212121142125213531416122aE A Et t 可得：40035353415093532211EKEK251035343127273323531233E KEK215251935313EK41253535323EK431274252151273321352594140012535035250215250254150191009353E Ke（3）位移列向量和右端项由边界条件可确定：Teu a00022水压力和构件厚分别为：10tgh p TTet l q h q h q R 032031020306000001自重为W 与支座反力：Ty x y x e W R R W W R R R 330333112所以：Ty x y x eW R h q R W h q W R R R33363303011由eeeRa K得到下列矩阵方程组：3336300030301122W R h q R W h q W R R u y x y x 化简得：431274252151273321352594140012535035250215250254150191009353E Ke364035353022W h q u E可得：EW E h q u 363567022将22u 代入下式：333425135025103533031122W R h q R W R R u E y x y x 固定面上地反力：ahga gh q 330从而可得支座反力为：43221234120303011h q W Rh q W R W h q R WR y x y x 这是y 地线性方程，但相容方程要求它有无数多地解（全柱内地y 值都应该满足它），可见它地系数和自由项都应该等于零，即资料个人收集整理，勿做商业用途0)(414dxx f d ，)(424dxx f d 这两个方程要求ICx Bx Axx f 231)(，KJx Ex Dxx f 232)(代入应力函数表达式，并略去对应力分量无影响地一次项和常数项后，便得2323)(ExDxCx BxAxy 对应应力分量为22yxgyEDx B Ax y xy26)26(22CBx Axyx xy2322以上常数可以根据边界条件确定.左边，0x ，1l，0m ，沿y 方向无面力，所以有)(C xxy右边，b x ，1l ，0m ，沿y 方向地面力为q ，所以有qBb Ab bxxy23)(2上边，0y ，0l ，1m ，没有水平面力，这就要求xy 在这部分边界上合成地主矢量和主矩均为零，即)(00dx y b xy将xy地表达式代入，并考虑到C=0，则有)23(23232BbAbBxAxdx Bx Ax b b 而00)(dx ybxy自然满足.又由于在这部分边界上没有垂直面力，这就要求y在这部分边界上合成地主矢量和主矩均为零，即0)(00dx y b y，)(00x d x y b y将y地表达式代入，则有02323)26(22Eb DbEx Dx dx E Dx b b22)26(2323EbDbExDxxdx E Dx b b 由此可得2bq A，bq B，0C ，0D ，0E 应力分量为0x,gy bx by q y312,23bx bx q xy虽然上述结果并不严格满足上端面处（y=0）地边界条件，但按照圣维南原理，在稍远离y=0处这一结果应是适用地.资料个人收集整理，勿做商业用途8、证明：如果体力分量虽然不是常量，但却是有势地力，即体力分量可以表示为xV f x，yV f y，其中V 是势函数，则应力分量亦可用应力函数表示为，V yx22，V xy22，yx xy2，试导出相应地相容方程.资料个人收集整理，勿做商业用途证明：在体力为有势力地情况下，按应力求解应力边界问题时，应力分量x，y，xy应当满足平衡微分方程yV xyx V yxxyyyxx（1分）还应满足相容方程y f x f yxy x yx 12222（对于平面应力问题）yf xf yxy x yx 112222（对于平面应变问题）并在边界上满足应力边界条件（1分）.对于多连体，有时还必须考虑位移单值条件.首先考察平衡微分方程.将其改写为0xVyyV x xyyyxx这是一个齐次微分方程组.为了求得通解，将其中第一个方程改写为yxxyVx根据微分方程理论，一定存在某一函数A （x ，y ），使得yA Vx，xAyx同样，将第二个方程改写为yxyxVy（1分）可见也一定存在某一函数B （x ，y ），使得xB Vy，yByx由此得yB xA 因而又一定存在某一函数y x,，使得y A，xB代入以上各式，得应力分量V yx22，V xy22，yx xy2为了使上述应力分量能同量满足相容方程，应力函数y x,必须满足一定地方程，将上述应力分量代入平面应力问题地相容方程，得资料个人收集整理，勿做商业用途Vyx Vx Vyyx2222222222221VyxV yxxy yx222222222222222212简写为V24)1(将上述应力分量代入平面应变问题地相容方程，得Vyx Vx Vyyx22222222222211VyxVyxxy yx2222222222222222112简写为V241219、如图所示三角形悬臂梁只受重力作用，而梁地密度为，试用纯三次地应力函数求解.O解：纯三次地应力函数为3223dycxyy bx ax相应地应力分量表达式为dy cx xf yx x6222,gy by ax yf xy y2622,cybx yx xy222这些应力分量是满足平衡微分方程和相容方程地.现在来考察，如果适当选择各个系数，是否能满足应力边界条件.资料个人收集整理，勿做商业用途上边，0y ，0l ，1m，没有水平面力，所以有2)(bx yxy对上端面地任意x 值都应成立，可见b 同时，该边界上没有竖直面力，所以有6)(ax yy对上端面地任意x 值都应成立，可见a 因此，应力分量可以简化为dy cx x62，gy y，cyxy2斜面，tanx y ，sin 2cosl ，cos cos m ，没有面力，所以有0tantan x y xyyx y yx x lmml 由第一个方程，得sin tan 6sin4costan 2sintan 62dx cx cx dx cx 对斜面地任意x 值都应成立，这就要求tan64d c 由第二个方程，得sin sin tan 2cos tan sintan 2gx cx gx cx 对斜面地任意x 值都应成立，这就要求0tan 2g c （1分）由此解得cot 21g c（1分），2cot31g d从而应力分量为2cot2cot gy gx x,gy y,cotgy xy设三角形悬臂梁地长为l ，高为h ，则lh ta n.根据力地平衡，固定端对梁地约束反力沿x 方向地分量为0，沿y 方向地分量为glh 21.因此，所求x 在这部分边界上合成地主矢应为零，xy应当合成为反力glh 21.资料个人收集整理，勿做商业用途cotcotcot2cot2202gh glh dygy gl dyh lxh xglhgh dygy dyh h lx xy21cot21cot2可见，所求应力分量满足梁固定端地边界条件.10、设有楔形体如图所示，左面铅直，右面与铅直面成角，下端作为无限长，承受重力及液体压力，楔形体地密度为1，液体地密度为2，试求应力分量.资料个人收集整理，勿做商业用途解：采用半逆解法.首先应用量纲分析方法来假设应力分量地函数形式.取坐标轴如图所示.在楔形体地任意一点，每一个应力分量都将由两部分组成：一部分由重力引起，应当与g 1成正比（g 是重力加速度）；另一部分由液体压力引起，应当与g2成正比.此外，每一部分还与，x ，y 有关.由于应力地量纲是L-1MT-2，g 1和g 2地量纲是L-2MT-2，是量纲一地资料个人收集整理，勿做商业用途量，而x 和y 地量纲是L ，因此，如果应力分量具有多项式地解答，那么它们地表达式只可能是gx A 1，gy B 1，gx C2，gy D2四项地组合，而其中地A ，B ，C ，D 是量纲一地量，只与有关.这就是说，各应力分量地表达式只可能是x 和y 地纯一次式.资料个人收集整理，勿做商业用途其次，由应力函数与应力分量地关系式可知，应力函数比应力分量地长度量纲高二次，应该是x 和y 纯三次式，因此，假设资料个人收集整理，勿做商业用途3223dycxyy bx ax相应地应力分量表达式为dy cx xf yx x6222,gy byax yf xy y12226,cybx yx xy222这些应力分量是满足平衡微分方程和相容方程地.现在来考察，如果适当选择各个系数，是否能满足2g1gyxO应力边界条件.资料个人收集整理，勿做商业用途左面，0x ，1l，0m ，作用有水平面力gy 2，所以有gydy xx26)(对左面地任意y 值都应成立，可见62gd同时，该边界上没有竖直面力，所以有2)(cy xxy对左面地任意y 值都应成立，可见c 因此，应力分量可以简化为gy x2，gy byax y126，bxxy2斜面，tan y x ，cos l ，sin2cosm ，没有面力，所以有0tantan y x xy yy x yx x lmm l 由第一个方程，得sin tan 2cos 2by gy 对斜面地任意y 值都应成立，这就要求sin tan 2cos 2b g 由第二个方程，得sin sin4sin tan 6cos tan 2sin 2tan611y g b a by gy byay 对斜面地任意x 值都应成立，这就要求4tan61g ba 由此解得321cot31cot61g g a，22cot21g b 从而应力分量为gy x 2,y g g xg g y 122321cotcot2cot ,22cotgx xy 位移边界条件对称、固定边和简支边上支点地已知位移条件如下：对称轴: 法线转角=0固定边: 挠度=0 (或已知值)边线转角=0 (或已知值)法线转角=0 (或已知值)简支边: 挠度=0 (或已知值)边线转角=0 (或已知值)计算图示四边固定方板方板地边长为l ，厚度为t ，弹性模型量为E ，波松比μ=0.3，全板承受均布法向荷载ｑ，求薄板中地挠度和内力. 资料个人收集整理，勿做商业用途单元划分：为了说明解题方法，采用最简单地网络2×2，即把方板分成四个矩形单元.由于对称性，只需计算一个单元，例如，计算图中有阴影地单元，单元地节点编号为１，２，３，４.此时，单元地a, b 是4l ba 计算节点荷载:由前面地均布荷载计算公式得：Tl l l l l l l l qlR ]21121212[192}{2边界条件：边界23和34为固定边，因此节点2, 3, 4地挠度、边线和法线转角均为零.边界12和14为对称轴，因此θx1 =0、θy1 =0.于是，在4个节点和12个位移分量中,只有一个待求地未知量1w .资料个人收集整理，勿做商业用途结构地代数方程组：这是一个单元地计算题目，单元刚度矩阵在此处即为总刚度矩阵.引入支承条件后，在总刚度矩阵中只取第一行、列元素，在方程组右端项中只保留第一个元素.于是结构地代数方程为：16)681(15815821201120qlw l D w k lD 资料个人收集整理，勿做商业用途同此解出04100148.0D ql w .其中32309158.0)1(12EtEt D 内力:利用式(4-2-6)可求得方板中点力矩为：由表看出，网格越密，计算结果越接近于精确答案.还可看出，位移地精度一般比内力地精度高，这是因为在位移法中，位移是由基本方程直接求出地，而内力则是根据位移间接求出地.资料个人收集整理，勿做商业用途第三章平面问题有限单元法习题答案3-2图示等腰直角三角形单元，设=1/4，记杨氏弹性模量E ，厚度为t ，求形函数矩阵[N]、应变矩阵[B]、应力矩阵[S]与单元刚度矩阵[K]e.资料个人收集整理，勿做商业用途【解】：ijmj imi j ji mm i j i m j m i i m j j m i m j i j m m j ix x c y y b y x y x a x x c y y b y x y x a x x c y y b y x y x a ,,,,,,aj(0,a)aac a a b a aa a a ac b a a a c a a b a a mmmj j j i i i 0,0,0*0*0,00,00**0000,0,0*00*02mjim j i N N N N N N N000),,()(21m j i y c x b a AN i i ii 221001010121aa a Ayxayxy x ay x a Na yx a ay ax aaN a y ay x a N a x y ax a N m j i 000001)(1)00(1)00(122221000310131031001310311103)411(2412100141141411411)4121)(411()411()1(2210011011)21)(1()1(EE E E D321B B B Baa c a ab a aa a a a cb a a ac a a b a a mmmj jj i i i 0,0,0*0*0,000,00**0000,0,0*00*0211011010100001000111110011011000110000110000100212aBaB aB a a a aB b c c b AB mjii i i i i1003101310E D1101101010000100011aB11011313001*********11110100001000110003101310aE a E BD S 1003101310E D11111010000100011aB42311124111331300111011011011013100320211101101010000100011000310131101101010000100011022Et at aE tAB D BKTTe3-3正方形薄板，受力与约束如图所示，划分为两个三角形单元，=1/4，板厚为t ，求各节点位移与应力.【解】：yP 34242311124111331300111011011011013100320Et tAB D BKTe0000000000000000003001310001101100011011001003130031114200111324201Et K4211310024131100111001001303100031013000111001000000000000000000202Et K4211310241311001140023113042011310240111120041300311142001113242021Et KKK 载荷向量：000000P R1001414004040042000000004211310241311001140023113042011310240111120041300311142001113242013344332211P v u Et P v u v u v u v u Et 101414041PEt PEt v u 05010015330050000044332211Et P v u v u v u v u 10003101310E D111101010000100011a B1101110001000010111aB12BB31201010003101325000000110111000100001011000310131033221111atPatPEtP aE v u v u v u BD 1002101031013200500011111000100001110310131022334422atP atP Et P aE v u v u v u BD 3-4三角形单元i,j,m 地j ，m 边作用有如图所示线形分布面载荷，求结点载荷向量.【解】：面力移置公式：tdsp NRTe其中：mjim j i N N N N N N N000),,()(21m j i y c xb a AN i i ii 426,132,62*63*2352,426,26*22*5165,363,213*56*6mmmj j j i i i c b a c b a c b a 213431402212165136122121Aj(6,3)i(2,2)m(5,6)1q 2q yxo)46(131)342(131)321(131y x N y x N y x N mj i 所以:yx yx yx y x yx y x N460342033004603420321131载荷分布函数：0)6(3)(121y q q q p积分函数：])6,5[(213x x y dyy q q q yxyx y x yx yxyxttdsy q q q yx y x yx yxyx y xRe3100)6(3)(460463420034233003211310)6(3)(464634200342330032113112163121dyy q q q y y y y t dyy q q q y yyy yy yy tRe)6(3)(133130013313263130026313000013*3100)6(3)(473160473163283420032834200013*3101216312163dyy q q q q y y q q q q y tdyy q q q q y y q q q q y tRe63121212126312121212))(36(*30))(36(*60027100)3)(2(*133130)3)(2(*26313013*310126323122126312631212632312212631263129292331)(321)36(3)(3)36(299331)(621)36(6)(6)36(q q y yq q yyq q dyy y q q dy y q q q q yyq q yyq q dy y y q q dy y q q 所以:210210031002182902990027100)3)(2(*133130)3)(2(*26313013*310121212126312121212q q q q tq q q q t dyy q q q q y y q q q q ytRe3-5图示悬臂深梁，右端作用均布剪力，合力为P ，取=1/3，厚度为t ，如图示划分四个三角形单元，求整体刚度方程.资料个人收集整理，勿做商业用途【解】：13524612341000420248410012102112)311(23121001311310311311)3121)(311()311()1(22100011011)21)(1()1(EE E E D10420248E D1111101000010001B53411235211442400211011011011024200416211101101010000100011020410241101101010000100018Et t E tAB D B KTTe534112352114424002110110110110242004164321Et KKKK0000000000000000000000000000000000000000000000400042020000010011100000000000000000000000000000410053120000210035140000010011100000200024041K534112352114424002110110110110242004164321Et KKKK00000000000000000000000000000000000000000000001011000100000424002000001253004100001435002100000000000000000000000000002420040000010110001162Et K00000000000000000000000000000000000000000000000080008404000002002220000000000000000000000000000082001062400004200610280000020022200000400048081612Et KK80844000020022200000000000000000000000000000820010624000042006102800000200222000004000480800000000000000000000000000000000000000000000000000001643Et KK808440200222000000000000000000000000000008200186248404420061228222002002220000040004808000000008200106240000420061028000002002220000040004808164321Et K K K K K算例2：正方形薄板平面应力问题地求解已知图示正方形薄板，沿其对角线承受压力作用，载荷沿厚度为均匀分布，P=20kN/m.设泊松比u=0，板厚t=1m ，求此薄板应力.资料个人收集整理，勿做商业用途课本第42页3.7节计算结果如下：变形：76.176.172.388.052.1252.32653321u u v u v v 应力：)/(40.40.2088.021m kN xyy x ；)/(052.1276.122m kN xyy x；)/(08.372.388.023m kN xyy x ；)/(32.172.3024m kN xyy x 1、如图1所示等腰直角三角形单元，其厚度为t ，弹性模量为E ，泊松比0；单元地边长及结点编号见图中所示.求(1)形函数矩阵N(2)应变矩阵B 和应力矩阵S (3)单元刚度矩阵eK1、解：设图1所示地各点坐标为点1（a ，0），点2（a ，a ），点3（0，0）于是，可得单元地面积为12A2a ，及(1)形函数矩阵N 为(7分)123aa12122121(0a a )a 1(00a )a 1(aa 0)aN x y N x y N x y ；123123N N N NI I I N N N (2)应变矩阵B 和应力矩阵S 分别为(7分)12a 010-a a-aaB ，220010a aa 0B ，32-a 0100a-aB ；123B B B B 12a00-a a11-a a 22E S ，22000a a1a 02E S ，32-a 000a10-a 2E S ；123123SD B B B S S S (3)单元刚度矩阵eK(6分)111213T21222331323331102113120111100140202002000201111eEt tAK K K KB DB K K K K K K 2、图2（a ）所示为正方形薄板，其板厚度为t ，四边受到均匀荷载地作用，荷载集度为21/N m ，同时在y 方向相应地两顶点处分别承受大小为2/N m 且沿板厚度方向均匀分布地荷载作用.设薄板材料地弹性模量为E ，泊松比0.试求资料个人收集整理，勿做商业用途(1)利用对称性，取图（b ）所示1/4结构作为研究对象，并将其划分为4个面积大小相等、形状相同地直角三角形单元.给出可供有限元分析地计算模型（即根据对称性条件，在图（b ）中添加适当地约束和荷载，并进行单元编号和结点编号）.资料个人收集整理，勿做商业用途(2)设单元结点地局部编号分别为i 、j 、m ，为使每个单元刚度矩阵eK 相同，试在图（b ）中正确标出每个单元地合理局部编号；并求单元刚度矩阵eK .资料个人收集整理，勿做商业用途(3)计算等效结点荷载.(4)应用适当地位移约束之后，给出可供求解地整体平衡方程（不需要求解）.图13①②③④2、解：(1)对称性及计算模型正确(5分) (2)正确标出每个单元地合理局部编号(3分) (3)求单元刚度矩阵eK(4分) (4)计算等效结点荷载(3分)(5)应用适当地位移约束之后，给出可供求解地整体平衡方程（不需要求解）.(5分)如图3.11所示地平面三角形单元，厚度t=1cm ，弹性模量E=2.0*105mpa ，泊松比γ=0.3，试求插值函数矩阵N ，应变矩阵B ，应力矩阵S ，单元刚度矩阵Ke.资料个人收集整理，勿做商业用途图2j m m mmi ii ij j j 1N /m21N /m 12456对称1011012020031214301201eEt K对称123356322000026121006120146101620212v v u Et tv u u解：此三角形单元可得：2△=（10-2）*4=32，故有a1=1/32*（8u1-5u2-16u3）a2=1/32*（4u1-4u2）a3=1/32*（-8u1+8u3）a4=1/32*（56v1-8v2-16v3）a5=1/32*（-4v1+4v2）a6=1/32*（-8v1+8v3）而b1=y2-y3=-4 b1=x2-x3=-8b1=y3-y1=4 b1=x3-x1=0b1=y1-y2=0 b1=x1-x2=8b1 0 b2 0 b3 0 -4 0 4 0 0[B]=1/2△* 0 c1 0 c2 0 c3 =1/32* 0 -8 0 0 8资料个人收集整理，勿做商业用途c1 b1 c2 b2 c3 b3 -8 4 0 8 01 γ 0 1 0.3 0[D]=[E/(1-γ2)]* γ 1 0 =[E/0.91]* 0.3 1 0资料个人收集整理，勿做商业用途0 0 (1-γ)/2 0 0 0.351 0.3 0 -0.125 0 0.125 0 0[S]=[D]*[B]={E/0.91}* 0.3 1 0 * 0 -0.25 0 0 0.25资料个人收集整理，勿做商业用途0 0 0.35 -0.25 0.125 0 0.25 01.4 0 -1.4 -0.7 0 0.70 4 -0.6 -4 0 0[K]①=BT*D*B①*t*△={E/36.4}* -1.4 -0.6 2.4 1.3 0.6 0.7资料个人收集整理，勿做商业用途-0.7 -4 1.3 -0.6 -1 0.350 0 0.6 -1 -0.6 00.7 0 0.7 -0.35 0 01 0 0 0.6 -1 -0.60 0.35 0.7 0 -0.7 -0.350 0.7 1.4 0 -1.4 -0.7[K]②=BT*D*B②*t*△={E/36.4}* 0.6 0 0 4 -0.6 -4资料个人收集整理，勿做商业用途1 -0.7 -1.4 -0.6 2.4 1.30.6 -0.35 -1.4 -4 1.3 3.53.12 求下图中所示地三角形地单元插值函数矩阵及应变矩阵，u1=2.0mm，v1=1.2mm，u2=2.4mm，v2=1.2mm，u3=2.1mm，v3=1.4mm，求单元内地应变和应力，求出主应力及方向.若在单元jm边作用有线性分布面载荷（x轴），求结点地地载荷分量.资料个人收集整理，勿做商业用途解：如图2△=64/3，解得以下参数：a1=19 a2=-2 a3=6；b1=-3 b2=4 b3=-1；c1=-1 c2=-3 c3=4；资料个人收集整理，勿做商业用途N1={64/3}*(19-3x-y) N2={64/3}*(-2-3x-3y)N3={64/3}*(6-x+4y)故N=Ni 0 Nj 0 Nm 00 Ni0 Nj0 Nm1 0 1 0 1 0 =0 1 0 1 0 1bi 0 bj 0 bm 0[B]={1/2△}* 0ci 0 cj 0 cm ci bi cj bj cm bm -3 0 4 0 -1 0={64/3}*0 -1 0 -3 0 4 -1 -3 -3 4 4 -11 γ 0[D]={E/(1-γ2)}*γ 1 00 (1-γ)/21 γ 0 -3 0 4 0 -1 0 单元应力矩阵[S]=[D]*[B]={E/13(1-γ2)}* γ 1 0* 0 -1 0 -3 0 4资料个人收集整理，勿做商业用途0 (1-γ)/2 -1 -3 -3 4 4 -12 1.1-3 -u 4 3u -1 4u2.4单元应力[δ]=[S]*[q]= {E/13(1-γ2)}* -3u -1 4u -3 -u 4* 1.2资料个人收集整理，勿做商业用途(u-1)/2 (3u-3)/2(3u-3)/2 2-2u 2-2u (u-1)/22.4资料个人收集整理，勿做商业用途1.43.13解：二维单元在x,y 坐标平面内平移到不同位置，单元刚度矩阵相同，在平面矩阵180°时变化，单元作上述变化时，应力矩阵不变化.3.14解：令1t，1p ，而E2.0e 011，1/3，210101102E D（0,1）（2,1）x y①②（0,0）（2,0）12312311223300000b b b Nc c c c b c b c b 2NBA单元①2.250.7500.75 2.25000.75D①②0.500.50000100010.500.51B①-1.125-0.75 1.125000.751.0+011*-0.375-2.250.37502.25-0.75-0.37500.3750.75Se ①S DB1.31250.75-0.5625-0.375-0.75-0.3750.752.4375-0.375-0.1875-0.375-2.25-0.5625-0.3750.562500.375*1.0011-0.375-0.187500.18750.3750-0.75-0.37500.3750.750-0.375-2.250.3752.25kee ①单元②：00.500.50B101001010.50.5②00.75 1.1250.75 1.125002.250.375 2.250.3750*1.00110.7500.750.37500.375Se ②0.7500.750.37500.37502.250.3752.250.3750.750.3751.31250.750.56250.3750.375 2.250.75 2.43750.3750.187500.3750.56250.37510.562500.37500.3750.18750.1875ke②由ke①和ke ②扩充KZ （总刚度阵）1.31250.750.56250.3750.750.375000.752.43750.3750.18750.375 2.25000.56250.3751.312500.75000.3750.3750.18750 2.437502.250.37501.01011*0.750.3750.750 2.06250.750.56250.3750.375 2.250kz e 2.250.75 4.68750.3750.18750000.3750.56250.3750.56250000.37500.3750.187500.1875而Re .kz qe ，其中112211Re22Rx Ry Rx Ry ，112200qex y x y ，化简得：112201.312500.7500.11310 2.43750 2.250.596820.750 2.06250.7500.194702.250.754.687510.42432x y x y 则，11220.56250.3750.750.3750.11130.148100.18750.375 2.250.59680.95170.750.3750.56250.3750.19470.17420.3750.3750.18750.42430.0482Rx Ry Rx Ry 3.15如图所示有限元网格，cm a4，单元厚度mm t 1，弹性模量MPa E5100.2，泊松比3.0.回答下述问题：（1）结点如何编号才能使结构刚度矩阵带宽最小？（2）如何设置位移边界条件才能约束结构地刚体移动？（3）形成单元刚度矩阵并集成结构刚度矩阵.（4）如果施加一定载荷，拟定求解步骤.(1) (2) （3）资料个人收集整理，勿做商业用途解：1、节点编号如图(2)所示；2、如图(3)设置位移边界条件才能约束结构地刚体移动；3、如图(2)所示各节点地坐标为(以m 为单位)：1(0,0)，2(0.08,0)，3(0,0.04)，4(0.08,0.04 )，5(0,0.08)，6(0.08,0.08)，7(0,0.12)，8(0.08,0.12)资料个人收集整理，勿做商业用途解：单元号 1 2 3 4 56相邻结点1 3 4 5 5 72 2 5 4 6 63436 78对于单元号1：04.0321y y b ；04.0132y y b ；0213y y b ；08.0231x x c ；0312x x c ；08.0123x x c ；对于单元号2：04.0423y y b ；0342y y b ；04.0234y y b ；0243x x c ；08.0432x x c ；08.0324x x c ；对于单元号3：04.0354y y b ；0435y y b ；04.0543y y b ；0534x x c ；08.0345x x c ；08.0453x x c ；对于单元号4：04.0645y y b ；0564y y b ；04.0456y y b ；0465x x c ；08.0654x x c ；08.0546x x c ；对于单元号5：04.0765y y b ；04.0576y y b ；0657y y b ；08.0675x x c ；0756x x c ；08.0567x x c ；对于单元号6：04.0867y y b ；0786y y b ；04.0678y y b ；0687x x c ；08.0876x x c ；08.0768x x c ；平面三角形单元地面积均为1112321x x x 23210032.0m y y y 弹性矩阵均为0112E D12/)1(0003.0191.0100.21113.035.000应变矩阵11)5()3()1(021c b BBB110b c 220c b 220b c 330c b 33b c 2505.125.1225005.125.12002500025033)6()4()2(021c b BBB330b c 220c b 220b c 440c b 44b c 005.125.1200250002502505.125.12250应力矩阵)1()5()3()1(BD SSS2308.192418.84725.27100.1116154.99451.544835.1602418.84725.276154.9002308.190009451.544835.16)2()6()4()2(BD SSS2418.84725.27100.1116154.9002308.190009451.544835.162308.192418.84725.276154.99451.544835.16单元刚度矩阵tA SBKKKT)1()1()5()3()1(3297.07692.03846.05495.07143.03187.1100.181978.23846.01923.03297.03901.27143.03297.0005495.03297.05495.003846.01923.001923.03846.007692.03846.003846.07692.01978.2003297.01978.23297.0t A SBKKKT)2()2()6()4()2(3297.05495.03297.0005495.0100.181923.03846.003846.01923.003846.07692.007692.03846.001978.23297.01978.2003297.07143.03187.13297.07692.03846.05495.03901.27143.01978.23846.01923.03297.0结构刚度矩阵为：。
弹性力学大作业
Computational Elasticity and Plasticity Project-2012 FallProject DescriptionMeeting the objectives of the Graduate course, Computational Elasticity and Plasticity in Civil Engineering Materials, the course project as described herein enables the enhancement and practicing of theory of Elasticity and Plasticity learned in the class. As the project is carried out in parallel with the class the students will have to digest the essentials of elastic and plastic constitutive laws that are the core issue of this course.Your tasks:(a) Assuming an isotropic material, derive the elasto-plastic incremental stiffness matrix using V on Mises yield surface, associated flow rule and considering the isotropic strain-hardening rule based on the effective plastic strain. Implement this model with MATLAB. The Von Mises yield surface is as follows:20F J k =-= (1) The uniaxial stress-strain relationship is as follows: 0000p EE E σεσσσσσεσσ=≤-=+> (2) The elasto-plastic stress-strain relationship is as follows:[][][][]x x y y z z ep xy xy yz yz zx zx ep p d d d d d d C d d d d d d C C C σεσεσετγτγτγ⎧⎫⎧⎫⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪=⎨⎬⎨⎬⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎩⎭⎩⎭=+ (3) The elastic moduli matrix is as follows:11.112000[]2(1)(12)120000212000002v v v sym v v v v E C v v v v -⎡⎤⎢⎥-⎢⎥⎢⎥-⎢⎥-⎢⎥=⎢⎥+-⎢⎥-⎢⎥⎢⎥⎢⎥-⎢⎥⎣⎦(4) The plastic moduli matrix is as follows:22222222.1[](1).x x y y x z y z z p x xy y xy z xy xy x yz y yz z yz xy yz yz x zx y zx z zx xy zx yz zx zx xxy xz ij y yz z s s s s sym s s s s s E C H v s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s sym s ⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥=-⎢⎥+⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎡⎤⎢⎥=⎢⎥⎢⎥⎣⎦(5)222()331e e p d E H d vσσε=++ (6) Where, s ij is the stress deviator tensor.Please prove eq.(5) and (6).(b) Apply this newly implemented model to predict the stress development and the stress-strain relationship of one point material behavior subjected to given strain-controlled paths. Analyze the features of your results .path 1: strain-controlled uniaxial monotonic loading::00.05x ε→path 2: strain-controlled uniaxial cyclic loading::00.00100.00100.0100.0100.05x ε→→→-→→→→-→→path 3: strain-controlled pure shear monotonic loading::00.05xy ε→path 4: strain-controlled pure shear cyclic loading::00.00100.00100.0100.0100.05xy ε→→→-→→→→-→→path 5: strain-controlled shear and tension monotonic loading::00.05xy x and εε→path 6: strain-controlled shear and tension cyclic loading::00.00100.00100.0100.0100.05xy x and εε→→→-→→→→-→→Assume that E =200GPa, v =0.3, σ0=300MPa and E p =2000MPa.Project EvaluationWe hope that every student in this course will complete the project independently. Upon the completion every student presents the project in a concise report. The project should be completed before January 25, 2013 and every student presents the project in both printed and electronicversion with your MATLAB coding.Appendix 1: Unified Variable NamelistMATLAB is a high-level language and interactive environment that enables you to perform computationally intensive tasks faster than with traditional programming languages such as C, C++, and FORTRAN. Compared to UMAT (User-defined Material Mechanical Behavior) in ABAQUS, variables can be used in MATLAB without previous declaration and you can use variable names according to your own habits. However, in order to make your program code more readable, we hope that all of you use the following important variable names.DDSDDE(NTENS, NTENS): a NTENS*NTENS matrix representing the elastic-plastic stiffness matrix ∂Δσ/∂Δε，Δσ is incremental stress，Δε is incremental strain. DDSDDE(i, j) defines the change in the I th stress component at the end of the time step caused by an infinitesimal perturbation of the J th component of the strain increment array. Generally speaking it is symmetric matrix.DDEDDS(NTENS, NTENS): a NTENS*NTENS matrix representing the elastic-plastic compliance matrix ∂Δε/∂Δσ，Δσ is incremental stress，Δε is incremental strain. DDEDDS(i, j) defines the change in the I th strain component at the end of the time step caused by an infinitesimal perturbation of the J th component of the stress increment array. Generally speaking it is symmetric matrix.STRESS(NTENS): Stress matrixSTRAIN(NTENS): Strain matrixDSTRESS(NTENS): Stress incremental matrixDSTRAIN(NTENS): Strain incremental matrixPROPS(NPROPS): Material parameter matrixAppendix 2:Formulation of incremental elastic-plastic constitutive relationships, consistency condition check. Please refer to the reference book, Elasticity and Plasticity, 2005, pp460-472, from WF Chen et al.Von Mise yielding surface, equivalent plastic strain and isotropic hardening rules. Please refer to the reference book, Elasticity and Plasticity, 2005, Chapter 7, from WF Chen et al..。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
弹性力学有限元大作业
一、模型信息：
已知：材料为铝合金。
E=71GPa ，v=0.3.
矩形平板的几何参数：板长为480mm ，宽为360mm ，厚度为2mm ；图形如下图；
加肋平板：
二、matlab 编程实现
1、程序相关说明：
计算使用的软件为：matlab2010a 主函数：main.m 主要计算部分
子函数：Grids.m 生成网格，节点数为：+1*+1I J （）（）
、单元数： 2**I J AssembleK.m 将单元刚度矩阵组装成总刚度矩阵（叠加方法）
GenerateB.m 生成单元格e B 矩阵 GenerateS.m 生成单元格e S 矩阵 GenerateK.m 生成单元刚度矩阵
2、网格划分：
利用Grid.m 子函数，取2020I J ==、，即可以得到网格如下： 节点数为：441个，单元格数：800个
3、计算过程及结果 （1）、网格划分：通过Grid.m ，生成节点数为：441个、单元格数：800个的网格 （2）、生成总刚度矩阵K ：通过GenerateK.m 、AssembleK.m 生成总刚度矩阵 采用常应变三角单元，e e u N a =,易得=e e B LN
由平面应力问题，可以确定2101011002E D νννν⎡⎤
⎢⎥
⎢
⎥=⎢⎥-⎢⎥
-⎢⎥⎣⎦
即e e S DB =
单元刚度矩阵为：e eT e K AtB DB = 总刚度矩阵为：eT
e e e
K G
K G =
∑
（3）、求解过程：
系统平衡方程为：Ka P = 将方程进一步划分为：E EF E E E T F F EF
F K K d f r d f K K +⎡⎤⎡⎤⎡⎤
=⎢
⎥⎢⎥⎢⎥
⎣⎦⎣⎦
⎣⎦ 通过已知边界条件（位移、载荷），确定E E F d f f 、、 ，从而将K 矩阵划分为四个模
块：E EF T
EF F K K K K ⎡⎤⎢⎥⎣⎦
1
()
E E E E
F F E T
F F F EF E r K d K d f d K f K d -=+-=-支反力：部分位移：
即整体位移向量为：E F d a d ⎡⎤
=⎢⎥⎣⎦
整体力边界条件为：E E F f r P f +⎡⎤
=⎢
⎥⎣⎦
（4）后处理：（应力、应变、抹平）
a、单元应力、应变：
e e e e e e
S a
B a σ
ε
=
=
b、抹平得到节点应力、应变：将每个节点参与组成的单元应力、应变叠加，然后除以叠加的单元数，得到抹平后的节点应力、应变。
（5）计算结果：
由于K矩阵行、列数过多，故以附件形式放在excel表格中
具体数据见附件：计算结果.xlsx
题目1：
（1）位移场
节点位置图：蓝色为初始位置，红色为最终位置
x、y的位移场：
（2）、应力场
（3）、应变场
题目2
（1）位移场：
节点位置：蓝色为初始位置，红色为最终位置
x、y的位移场：
（2）、应力场：（分别为x y xy σσσ、、 ）
（3）、应变场（分别为x y xy εεε、、 ）
三、Abaqus 分析
1、平板：（文件为：plate1.cae ；plate2.cae ）
题目1：
（1）位移场
（2）应力场（分别为x y xy σσσ、、 ）
（3）应变场（分别为x y xy εεε、、 ）
题目2：
（1）位移场
X 方向： y 方向：
（2）应力场（分别为x y xy σσσ、、 ）
（3）应变场（分别为x y xy εεε、、 ）
2、加肋后的平板（文件为：plate3.cae ；plate4.cae ）
题目1：
（1）位移场
（2）应力场（分别为x y xy σσσ、、 ）
（3）应变场（分别为x y xy εεε、、 ）
题目2：
（1）、位移场：
（2）、应力场：（分别为x y xy σσσ、、 ）
（3）、应变场（分别为x y xy εεε、、 ）
四、误差分析
题目1：
1、位移场
Matlab计算结果：
Abaqus计算结果：
误差分析：由于采用的网格划分都是三角单元，并且节点都比较多，所以x、y方向的位移两者计算结果、变化趋势大致相似。
matlab计算：x方向最大值为0.0017mm，y方向最小值为-0.0100mm，abaqus计算：x方向最大值为0.001758mm，y方向最小值为-0.01mm。
x 方向的变形大致呈鼓形，y方向的变形大致呈线性（层状），可能由于边缘效应，两边会出现一定的角度。
2、应力场
Matlab计算结果：
Abaqus计算结果：
误差分析：两者变化趋势大致相同，但是计算结果误差较大，尤其是角点位置。
以x方向为例，matlab：x方向最大应力值5.4140e+003，最小应力值-8.7261e+005，abaqus：x方向最大应力值为5.98e+003，最小应力值-12.42e+005，这可能跟角点边缘效应、抹平方法在角点效果不大以及abaqus固支定义不同，所以误差比较大。
3、应变场
Matlab计算结果：
Abaqus计算结果：
误差分析：变化趋势大致相同，但是计算结果误差较大。
Matlab：x方向最大应变为：0.0083，y方向最大应变为-0.0252，abaqus：x方向最大应变为0.012，y方向最大应变为：-0.02463。
x方向误差明显比y方向的误差大。
题目2
（1）位移场
Matlab计算结果：
Abaqus 计算结果：
误差分析：变化趋势、计算结果较一致。
（2）应力场（分别为x y xy σσσ、、 ）
Matlab 计算结果：
Abaqus 计算结果：
误差分析：变化趋势一致，但是数值误差较大，abaqus 计算结果的角点影响更明显。
剪应变误差更大。
抹平的方法可能有所不同。
（3）应变场（分别为x y xy εεε、、 ）
Matlab 计算结果：
Abaqus 计算结果：
误差分析：变化趋势一致，但是数值误差较大，尤其是剪应变。