计算力学大作业报告

计算力学基础大作业报告
航21班苏浩尹思凡
一、求解流程说明
本次作业的求解过程主要分为三个步骤：
1. 前处理
网格划分——采用ABAQUS软件进行画网格，并输出前处理文件。

形成前处理文件——利用matlab编程，形成inmesh文件。

2. 有限元计算
采用上机课上提供的FEATP程序计算，程序略有修改，使其能够计算单元和节点数目更多的算例。

本次大作业所有的图中str1~3,分别代表x（r）、y(theta)方向的正应力，剪力。

3. 后处理
采用Tecplot软件画图，得到位移场和应力场的云图。

二、所编程序及功能简介
二、具体算例
1.悬臂梁问题
1.1 问题描述
图1所示悬臂梁，一端固支，承受集中力 P=1000N，梁的长度为 100mm，高度为10mm，厚度为1mm，材料弹性模量为2×1011Pa，泊松比为0.32，假定在变形过程中界面始终保持完整，利用有限单元法计算位移场分布，并利用求得的位移场计算结点应力，讨论与理论解的误差。

分别使用不同类型不同阶次单元进行计算，比较结果精确度，单元收敛性等。

1.2 解析计算
采用材料力学方法进行计算，建立主轴坐标系，根据材料力学的推导，有如下公式：
()0()z x z
y y z xy z M x y I Q S y I t
σστ=-==-
（1.1）
其中()z M x 是z 轴正方向的弯矩，y Q 是方向为y 轴正方向的剪力，z I 是梁对z 方向的惯性矩，()z S y 是横截面对z 轴的静面矩，t 为厚度。

由题目中的数据，可以得到：
334
222
23
2
21000()()1000(100)()10125012123
111()()()(25)2242y z z y
y z h h Q P N
M x P L x x N mm h t I mm h S y ytdy ty t y y mm --=-=-=--=-⋅-⋅⨯======-=-⎰ （1.2）
将（1.2）中的公式带入（1.1）中，可以得到:
()212(100)0
6(25)()x y xy
x y MPa y MPa σστ⎧=-⎪⎪
=⎨⎪=-⎪⎩ (1.3) 材料力学中的挠度公式和转角公式：
22
()z z z M x d EI dx d dx
ν
νθ=
=
(1.4) 而x 方向的位移：
(,)()z u x y x y θ=- (1.5)
将（1.2）中的计算结果和材料参数带入，并对式（1.4）进行积分，得到：
521523121
610(100)2
1
610(50)
6
d x x c dx x x c x c νν--=-⨯⨯-+=-⨯⨯-++ （1.6）
位移边界条件：0,0,0dv
x dx
ν=== 得到位移场：
525231610(100)2
1
610(50)
6u x x y x x ν--=-⨯⨯-
=-⨯⨯- （1.7） 1.3 问题讨论
（1）位移边界条件：固定位移：左边界——中点 u=0 v=0 其余点 u=0 （2）载荷条件：右上角点
1.4 用不同类型的单元求解位移场和应力场 1.4.1 三角形三节点单元 （1）网格划分
采用三种不同疏密度的网格计算：
（2）计算结果
这里仅选取网格密度最大的第三种情况的结果 位移场,u v
σστ
应力场,,
x y xy
1.4.2 四边形四节点单元（1）网格划分
采用四种不同疏密度的网格计算：
（2）计算结果
位移场,u v
σστ
应力场,,
x y xy
1.4.3 四边形八节点单元
（1）网格划分
采用三种疏密度不同的网格计算
（2）计算结果位移场,u v
σστ应力场,,
x y xy
1.5 单元精度与单元收敛性
1）考虑平面应力问题，理论上v 的挠度位移：
3635PL PL v EI GA
=+ （1.8）
代入数据：P=1000N ，L=100mm E=2×105MPa G=E/(2（1+ν）)=75757.6MPa 得到理论挠度：
3353610001006100010020.1584135575757.6110
321011012
PL PL v mm EI GA ⨯⨯⨯=+=+=⨯⨯⨯⨯⨯⨯⨯⨯
以（
2）再对应力场进行检验
以（50,10）处的正应力值x σ作为收敛对象，理论()121005053000x MPa σ=⨯-⨯=
上表及上图中，三角形单元T3的节点数分别为63、205、1111和4221，四边形单元Q4和Q8的单元数相同，节点数分别为Q4：63、205、1111、4221;Q8:165,569,3221。

由对比可以看出，三角形三节点单元的解的误差要远大于四边形单元的解的误差，这是由于悬臂梁问题中，四边形单元的划分方式对位移的描述更符合实际的位移分布；二次单元的解的精度比一次单元解的精度高，且四边形八节点单元在节点数很少时收敛精度就已经很高，由材料力学解可以看出，悬臂梁问题中的位移场u是二次的，二次单元得到的与实际位移场阶数相同，因此二阶单元的精度很高是合理的。

2. 带孔方板问题
2.1 问题描述
图2所示中间含椭圆孔的矩形薄板，在上下两边受均布载荷 p作用，计算沿椭圆轴向的应力分布和应力集中系数。

p 取 50MPa，薄板长宽均为 20mm，圆直径为 5mm。

弹性模量
E=210GPa，泊松比v=0.32。

探讨下列变量（不局限于下列变量）对结果的影响：
（1）椭圆孔不同的长短轴之比 a/b
（2）板长与椭圆轴之比 s/a或s/b
（3）网格密度
2.2 解析计算
对于带圆孔的方板问题，解析解的应力图xσ如图所示（图示解析解的拉力方向是向右拉伸）
对于无限域椭圆孔的单向均匀拉伸问题（Kirsch问题），椭圆孔的应力分布存在解析解，椭圆
的长轴出应该出现应力集中，应力集中系数：2
1a
k
b
=+ 2.3 问题讨论
由于椭圆孔方板的对称性，本题只取板的1/4进行计算，且单元类型均为三角形T3单元。

且位移边界均根据对称性给出，左边界限制u方向位移，下边界限制v方向位移。

同时，对应力集中区域进行了网格加密。

2.3.1 椭圆孔不同长短轴之比a/b（s=20）
(1)a=5,b=1
网格划分：
计算结果：
位移场,u v
σστ应力场,,
x y xy
(2)a=5,b=2
网格划分：计算结果：位移场,u v
σστ应力场,,
x y xy
(3)a=5,b=3
网格划分：
计算结果：
位移场,u v
σστ应力场,,
x y xy
(4)a=5,b=5
网格划分：
位移场,u v
σστ应力场,,
x y xy
（5）分析与讨论
i)由上面四组位移和应力云图可以看出，在椭圆长轴附近出现明显的应力集中，因此椭圆长轴处为危险点，在工程应用中应该格外注意。

另外，从上面的云图可以看出，椭圆方板的短轴边上的x σ上方为拉应力，靠近短轴处存在压应力，这是因为上方均布的拉应力提供了一个相对于短轴点的弯矩，弯矩会在短轴处产生较大的压应力，与理论分析相符合。

ii)应力集中系数和轴向应力分布与椭圆长短轴之比的关系：
从下面的计算结果可以看出，在方板大小不变的情况下，随着椭圆的离心率的减小，椭圆长轴端点的应力集中系数逐渐减小，且趋近kirsch 解得应力集中系数3.0. 显然，由于平板不是无限大，应力集中系数均有一定偏离.
由应力y σ沿轴向分布图可以看出，随着所取得的x 轴上的点偏离长轴，y σ的大小呈指数
的大小已经接近平缓，之后的变化十分缓慢。

趋势迅速衰减，距离长轴点1mm处时，
y
应力沿轴向分布
2.3.2板长与椭圆长轴之比s/a (1)s=10, a=5, b=1
网格划分：
位移场,u v
σστ
应力场,,
x y xy
(2)s=15,a=15,b=1
网格划分：
位移场,u v
σστ应力场,,
x y xy
(3)s=20 a=5 b=1
网格划分：
位移场,u v
σστ
应力场,,
x y xy
(4)s=30,a=5,b=1
网格划分：
位移场,u v
σστ应力场,,
x y xy
(5)分析与讨论
沿x轴向应力分布
应力集中系数比较接近实际值，且板子越小，应力集中系数越接近理论值，这里的原因是所画的网格密度是相同的，板子越大，相对应力集中系数处的网格越稀疏，不利于得到较好的应力集中系数，上表所示板子较小时，由于网格很密，得到的应力集中系数与理论值11符合得很好。

3. 电子器件
3.1 问题描述
“微型电机”是微电子机械 MEMs的重要元件（图 3.1），它通常由多晶硅材料经过光刻得到，实际的“微型电机”比较复杂。

为简单起见，将电机转子等效为二维问题，转子厚度为13微米（具体几何尺寸见图 3.2）材料参数为：弹性模量 169GPa，泊松比 0.262，密度 2300 kg.m-3，假定电机在转动过程中仅受到离心力作用，试计算电机转速为20 转/分时转子内的应力分布。

3.2 问题讨论
微型电机的形状比较复杂，不能进行解析求解。

在r很小时，可以近似认为是厚壁圆筒问题，与轴对称问题的解接近，另外判断出角点应该存在应力集中。

考虑到微型电机的轴对称形
状和载荷（离心力）的轴对称分布，仅对1/4部分的微型电机进行计算。

3.2.1 网格划分
本算例采用三角形T3单元网格计算，采用疏密度不同的四种网格，单元数分别为：1103、2496、3983、5188，其中单元数为5188的网格划分如下图所示
3.2.2边界条件
1. 位移边界条件：由微型电机的对称性结构和离心力的对称条件，可以得到位移边界条件：
==
u y v x
(0,)0,(,0)0
2. 力边界条件
微型电机只受到离心力作用，相当于一个向外的体力作用。

接下来推导微型电机的受力情况：
对于三节点三角形单元有：
12
123123112233
2
2122
1121121,0,,0,,00,,0,,0,0,00,ve
A
e T T T b A
A A P N fdv N ftdA N ft J dL dL N N N N N N N N L N L N L r x f y N x x P t J dL dL L t J dL dL N y y ρωρωρωρω=
==⎡⎤=⎢⎥
⎣⎦===⎡⎤⎡⎤==⎢⎥⎢⎥⎣⎦⎣⎦⎡⎤⎡⎤⎡⎤==⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦
⎣⎦⎰⎰⎰⎰⎰⎰⎰⎰（1）
记常数
2Const t J ρω= （2）
其中
131311232322,,,,,,x y x y x x y y L L J x y x y x x y y L L ξξηη∂∂∂∂⎛⎫
⎛⎫ ⎪ ⎪--⎛⎫
∂∂∂∂ ⎪ ⎪
⎪=== ⎪ ⎪∂∂∂∂ ⎪--⎝⎭ ⎪ ⎪∂∂∂∂⎝⎭⎝⎭
（3） 可以得到：
1112
111223312
1112212312
2
1
1312231312
123()[(1)]()()]111
()
122424
A
A
A
A
P Const N xdL dL Const N N x N x N x dL dL Const L L x L x L L x dL dL Const L x x L L x x L x dL dL Const x x x ==++=++--=-+-+=++⎰⎰⎰⎰ （4）
根据（4）的结果，可以编程实现求解每个节点处的等效载荷。

3. 单位推导
电机为微尺度问题，需要进行单位的统一，选用如下的单位制：
质量[M]：mg ； 长度[L]：μm ；时间[t]：s
利用上述的单位制，可以推导出以下的导出单位： 应力：[][]
[]
1
2
1211M L T mg m s Pa μ----=⋅⋅=
密度： [][]3
3123110M L mg m kg m μ---=⋅=⋅ 所以，得到无量纲参数：
弹性模量916910E =⨯，泊松比0.262ν=，密度12
230010ρ-=⨯，转速23
n π=
因此，由离心体积力公式2F
r ρω=，可以得出离心力的值非常小，为了降低浮点数误差，
将转速提高1010倍，得到的体积力放大了2010，所得到的位移和应力乘以系数2010-即为计算得到的值。

3.2.3 计算结果（只取网格最密的情况） 1. 节点数2769，单元数5188 2. 笛卡尔坐标系下的计算结果 （1）位移,x y u v
σστ（2）应力,,
x y xy
3.极坐标转换
由于微型电机是在结构上是轴对称的（不是轴对称问题），因此得到极坐标下的位移和应力对于分析来说更有意义，利用坐标转换公式进行坐标转换：
cos sin sin cos r u u v v θθθθ
θ⎧⎫⎡⎤⎧⎫
=⎨⎬⎨⎬⎢
⎥-⎣⎦⎩⎭
⎩⎭ cos sin cos sin sin cos sin cos x
xy r r xy y r θθ
θστστθ
θθ
θτστσθθθθ-⎡⎤⎡⎤⎡⎤⎡⎤
=⎢⎥⎢⎥⎢⎥⎢⎥-⎣⎦⎣⎦
⎣⎦⎣⎦
得到：
cos sin sin cos r u u v v u v θθθ
θθ
=+⎧⎨
=-+⎩ ()()2222
22cos sin 2sin cos sin cos 2sin cos sin cos cos sin r x y xy x y xy r x y xy θθσσθσθτθθσσθσθτθθ
τσσθθτθθ⎧=++⎪⎪=+-⎨⎪=--+-⎪⎩
得到极坐标下的结果：
（1）位移,r u v θ
（2）应力,,r r θθσστ
3.2.4 分析与讨论 选取重要截面（10,22.5,45θ=），沿着这些路径确定r σ的值并进行分析。

从上图中可以看出：
1.在r很小的起始点，三个截面的径向应力rσ几乎相等，在r很小时，可以认为是轴对称，各截面的rσ应该相等。

2.在转子10°方向，即角点处有很大的应力集中，容易造成破坏，是实际应用中的危险点。

4. 厚壁圆筒
4.1 问题描述
承受均布内压的厚壁圆筒，尺寸如图 4 所示（单位 mm），弹性模量为3×105MPa，泊松比为0.35，内压为 20MPa，假定材料始终处于线弹性情况，试用FEATP程序的平面应变和轴对称问题进行求解，给出所用边界条件及其描述，并将计算结果和解析解进行比较，同时比较不同单元类型和网格划分对计算精度的影响。

4.2 解析计算 厚壁圆筒的lame 解：
2222222222222222222211110
i r i r o o p b a b a r p b a p a b b a r p a b a r r a b b θθσστ⎛⎫=- ⎪-⎝⎭⎛⎫=⎛⎫-- ⎪
-⎝⎭⎛⎫-+ ⎪-⎝⎝⎭
+ ⎪-⎭=
取0,15,30o p a b ===有：
22222
222222222222201530203011(1())3015320301(1())
30
i r i r p a b b a r r r p a b b a r r
θθσστ⎛⎫⎛⎫⨯=-=-=- ⎪ ⎪
--⎝⎭⎝⎭⎛⎫=+=+ ⎪-⎝⎭=
Matlab 编程实现：（参见程序：lilunjie4.m ）
4.3 平面应变解法
4.3.1 边界条件
位移边界条件：圆筒外部固定。

力边界条件：圆筒内壁受到侧压。

4.3.2 三角形单元T3
（1）网格划分
采用三种疏密度不同的网格进行计算，网格划分如下图所示：
三种网格的单元数分别为1600、2414和3454。

（2）计算结果
这里仅选取网格密度最大的第三种情况的结果笛卡尔坐标系下的计算结果
u v
（1）位移,x y
σστ
（2）应力,,
x y xy
与第三题做法相同，进行极坐标转换，得到极坐标下的结果：（1）位移,r u v
（2）应力,,r r θθσστ
4.3.3 四边形T4单元T8单元
（1）网格划分
采用三种疏密度不同的网格进行计算，网格划分如下图所示：
三种网格的单元数分别为729、1100、1560
（2）计算结果
四边形T4和T8单元得到的结果与三角形T3单元几乎相同，为节省篇幅在此不再列出。

4.4 轴对称解法
4.4.1 边界条件
位移边界条件：圆筒外部固定，直立圆筒高度取为60，下端固定。

力边界条件：圆筒内壁受到侧压。

4.3.2 三角形单元T3
（1）网格划分
采用三种疏密度不同的网格进行计算，网格划分如下图所示：
三种网格的单元数分别为1600、2414和3454。

（2）计算结果
这里仅选取网格密度最大的第三种情况的结果。

笛卡尔坐标系下的计算结果
u v
（1）位移,x y
σστ
（2）应力,,
x y xy
4.4.3 四边形T4单元T8单元（1）网格划分
采用三种疏密度不同的网格进行计算，网格划分如下图所示：
三种网格的单元数分别为36、100、900
（2）计算结果
四边形T4和T8单元得到的结果与三角形T3单元几乎相同，为节省篇幅在此不再列出。

4.5 分析与讨论
4.5.1 网格密度对计算结果的影响。

以三角形平面应变情况为例，在圆筒上取一个截面，其径向位移和径向正应力随径向坐标的关系如下图所示
由上两图可以看出，不同网格密度的单元对径向位移几乎没有影响，对于径向正应力有一定的影响，但比较小。

说明对于轴对称问题，单元收敛性是比较好的。

4.5.2 单元类型对计算结果的影响
分别取三角形T3单元和四边形Q4单元，采用网格密度最大的单元取相同截面，其径向位移和径向正应力随径向坐标的关系如下图所示
由上面两图可以看出，不同的单元类型对轴向位移和轴向正应力的影响都很小，说明对于厚壁圆筒问题，三角形T3单元和四边形T4单元所得到的结果都是比较好的。

4.5.3 平面问题和轴对称问题的比较
由上面两图可以看出，平面应变问题和轴对称问题算出的结果几乎是相同的，因此平面应变问题与轴对称问题在厚壁圆筒上的适用性相同。

5. 温度场计算
5.1 问题描述
5.2 解析计算
(1)根据问题的条件，二维稳态热传导无内热元问题的控制方程简化为：
2
2
22()()0(
)()0k k x x y y k k k x x y y x y
φφφφφφ
∂∂∂∂+=∂∂∂∂∂∂∂∂∂∂+++=∂∂∂∂∂∂
将10sin()sin()80.4k x y a
b
ππ
=+带入得
1210cos()sin()10sin()cos()k x y k x a a b
k x y k x b a b πππ
πππ
∂==∂∂==∂
即
221222()0k k k x y x y φφφφ∂∂∂∂+++=∂∂∂∂
(2)边界条件 左边界：0,100x φ== 右边界：0,200x φ== 上边界：,
0y b y
φ
∂==∂ 下边界：100,
80.4
y y φ∂==-∂ (3)基于有限差分法的节点离散方程
1. 内部节点离散方程，取x y ∆=∆
221222()0k k k x y x y φφφφ∂∂∂∂+++=∂∂∂∂
1,1,,1,1
1,,1,,1,,1
1
2
2
2
22(
)022()()m n m n
m n m n m n m n m n
m n m n m n t t t t t t t t t t k k k x
y
x y +-+--+-+---+-++++
=∆∆∆∆
2. 边界节点离散方程（用平衡法建立） 上边界：
1,,,1,1,,1,,11,,0
22
2(2,3, (1)
4m N m n m N m N m N m N m N m N m N
m n t t t t t t y y
k
k x k x y x t t t t m M --+--+---∆∆+∆+=∆∆∆++==-
下边界：
1,1,,2,11,,1,,11,,0
22
22/(2,3, (1)
4m m n m m m N m N m N m N m N m n t t t t t t y y
k
k x k q x x y x t t t q k t m M -+--+---∆∆+∆++∆=∆∆∆+++∆==-
左右边界：
(4)代数方程组的求解方法
⎧++++=⎪
++++=⎪⎨
⎪⎪++++=⎩1111221331121122223322
11
2233.............
...n n n n n n n nn n n a x a x a x a x b a x a x a x a x b a x a x a x a x b
将代数方程组改写为：1
112213311122222233221n 223311
(...)1(...). (1)
(...)n n n n n n n nn n n x b a x a x a x a x b a x a x a x a x b a x a x a x a ⎧=----⎪⎪
⎪=----⎪⎨
⎪
⎪⎪=----⎪⎩
1.设定初值 xi(0) i=1,2,3,…,n
2.将初值代入方程，计算新值 xi(1) i=1,2,3,…,n。