中文第14章划分网格的一般原则及计算精度分析

2. 应力磨平： <1> 单元应力磨平
(最简单的一种方法，3节点三角形)
σe const
可视单元内应力平均值，或形心处应力。 所以，平均应力= 1 (σe1 σe2 )
2
得到的平均应力是矩形大单元形心的值。
<2>节点应力磨平
σ 1
m
σe
m e1
<3> 加权磨平
m
σ weσe
应变、应力的精度比位移低一阶。
应力解的近似性体现在： •在单元内部一般不满足平衡方程； •单元与单元之间一般不连续； •在力的边界条件上一般不满足力的边界条件。
§14.5.2 应力精度的改进 1.应力计算公式
• 应力在单元内分布，对于p 完备的位移 模式，应力为p-1次分布
• Gauss积分点应力精度最好。
kd r
Substitute
d T d and r T r
to get
T T kT d r
k
•例：三维实体块单元与板单元之间的连接
板单元节点上有转角自由度，而实体单元没有。 把板单元的转角与实体单元的节点位移协调起 来。
例如：
例1：如图一悬臂梁不同单元的计算结果
例2：计算精度与单元网格划分的关系 如图一悬臂梁，梁的长L,高B, L/B=10，厚度 为1 ；另一端受沿Y方向的均布力P。E=10GPa， P=10KN
理论上，端部v方向的位移可计算得到：
当端部只有两个节点时，等效结点载荷为：P/2，
当端部是3个结点时(T6)，两段的等效节点载荷
是P/6，中间节点上为2P/3。
三角形单元的网格划分情况： M1至M5
Q4，Q8单元的网格划分情况： M1至M5
自由端点挠度的计算结果：
(NDLT：总自由度数)
网格的安排与布置
• 网格的大小和密度对于问题是否合适 (以便减 少计算时间) •选择的单元类型对于问题是否合适 (以便提高计 算精度) •单元的形状是否使单刚近似奇异 •选择的单元和划分的网格是否能够便于描述分 布力载荷
由于
A(σ,σ) 1 M
T
(σ σ) C(σ σ)dV
2 Ve e1
σ 为由节点位移计算得到的应力场。
M
变分得到 A(σ)
T
(σ σ) C σdV 0
Ve
e1
ne
注意到 σ Niσi
i 1
于是有：
A
σi
M

e1
T
(σ σ)
长细比要适当
b/h<=2~4. b 和 h 分别为单元的最长和最短
边长度。
•
•
•
•
Good
• •
Bad
单元形状要保持凸性 单元边的夹角不能接近 0°或 180°
•
•
•
•
•
•
•Hale Waihona Puke Baidu
•
Good
•
•
•
•
•
•
•
•
Bad
不连续处的网格自然划分
在几何、载荷、边界条件突然变化(间断)处布置 节点。
(a) 集中载荷
•对称线或对称面必须是单元的边或面，即对称线或 对称面不能穿越单元。
•在二维问题中，对称线上节点垂直于对称线方向位 移为0，沿对称线的转动分量为0。
•在三维问题中，对称面上节点只能在对称面内有位 移，对称面内的转动分量为0。
允
约
许
束
的
的
位
位
移
移
反对称(Antisymmetry)
反对称问题网格划分要点
建立有限元模型涉及：
•单元类型和形状的选择
•网格的安排与布置 单元类型和形状的选择
(1) 问题的维数：一、二、三维 (2)单元类型：结构单元(梁、板、壳单元)， 实体单元
(3) 单元形状 三角形单元比较适合不规则形状 四边形比较适合规则性状
(4)单元阶次 与求解域内应力变化特点有关
•单元类型要适合于问题的载荷与应力。 •单元要能反映由于载荷、几何和边界条件所导 致的重要的应力。
Ae (σ, σ)

1 2
T
(σ σ) (σ σ)dV
Ve
<6> 利用Gauss点应力计算节点处应力 对于等参元，Gauss积分点的应力值求得后，
将Gauss积分点(例如2×2)应力代入
n
σ(ξI ) Ni (ξI )σi , I 1, 2, 3, 4
i 1
即可解出角节点处的应力值 σi
a

u1 (u1 - u2 )/a - = 0
u2
平衡和协调
•节点上的力和力矩平衡条件满足，因为 KU=F •在节点上变形是协调的，这是由刚度阵的 组装方式决定的。
•在单元之间的边界上，力不平衡，但是随 着网格加密，不平衡程度逐渐减小。
•单元中心处的应力精度比单元边上的应力精 度要高。
•在某些高阶单元以及不同类型单元的连接处 可能存在单元边界上不协调情况，随着网格 加密这种情况的严重程度逐渐减小。
小结：对称性和反对称性
利用对称性以最小代价获得最高精度。 沿对称线或对称面把物体切割，然后在新
的边界上施加对称或反对称边界条件。 要求几何和材料对称。 要求载荷(源项)对称或者反对称。 要求位移边界条件对称或反对称。
§14.3 网格加密的h,p,r方法
网格加密
•在应力梯度比较高的地方使用细网格(通常位于不 连续处)。 •网格密度要平稳过渡，防止相邻网格尺寸相差悬殊。 •网格密度要在计算精度和计算量之间折衷。
例如，二维4节点等参元，利用四个Gauss积分点 计算四个角节点的应力值。
对于4个积分点
4
σ j Niσi,
i 1
Ni

1 4
(1 i )(1 i)
j I , II , III , IV
σI N1(ξI )

σ II



N1 (ξ II
)
σIII N1(ξIII )
Interpolate u5, v5 from u2, v2,u3, v3 :
u5

a L
u2

b L
u3
ab v5 L v2 L v3
5 rotation of edge 2-3
Find:
a

L
0
b L

0
a L
0
0 0 0 0

b L
0 0 0
T



cos
基本思想：使单元的变形在连接处协调。
4•
3•
y 1•
x
•2
Bilinear plane element (dofs: ui, vi)
Truss element (dofs: ui, vi)
•5
Compatible dofs: no problem
4• y 1•
3•
a
5•
bL
2•
Bilinear plane element
IV

a 1 3 / 2, b 1 / 2, c 1 3 / 2
(b) 分布载荷的突变
(c) 板厚度的突变
(d) 材料性质的突变
疏密网格的过渡 •不同密度划分网格过渡 (a) 采用形状不规则的单元过渡
缺点：可能因单元形状不好而影响局部精度。
(b) 采用三角形单元过渡 缺点：可能因引入不同形式的单元而带来不便。
(c) 采用多点约束方程过渡 需引入约束方程：
e1
面积加权
we
se
m
si
i 1
单元面积占单元片面积的权重
<4>整体磨平
由于位移元得到的应力场不连续，用整体应力 磨平改进计算结果，以使应力场全域连续。
(a)构造一个改进后的应力解，
ne
σ Niσi
i 1
ne是单元节点数， Ni是第i点的插值函数 σi 是改进后的节点应力， σ 是改进后的应力场函数。

sin
cos
sin
0 0 0
6,7 L
L
L
L

0
0
0

0 1 0 0
0
0
0
0 0 1 0
0
0
0
0 0 0 1
is the angle between line 2-3 and y axis
Then follow the usual scheme:
(dofs: ui, vi)
Frame element
(dofs: ui, vi, i)
•6
问题: • 节点 5 不与节点 2 和 3同坐标 • 自由度不同: - 1和4节点仅有平动自由度 - 5和6节点有平动和转动自由度
方法:
把刚架节点5上的三个自由度 ( u5, v5 and 5) 转换 为节点 2和 3上的自由度 (u2, v2, u3, v3).
•反对称面线或反对称面必须是单元的边或面，不能 穿越单元。 •在二维问题中，反对称线上的节点沿反对称线的平 动位移为0，垂直于反对称线的转动位移分量为0。 •在三维问题中，反对称面上的节点在反对称面内的 平动位移为0，垂直于反对称面的转动位移分量为0。
允
约
许
束
的
的
位
位
移
移
其它对称性 •周期对称(cyclic symmetry) •轴对称(axisymmetry)。
特性无关。 位移与材料特性相关。
支反力(力矩)的合力通常与载荷平衡。 对于分布力载荷，计算其合力，检查是否
与支反力平衡。
如果软件提供了误差功能，检查计算误差。
§14.5 应力计算结果的性质和处理
§14.5.1 应力近似解的性质 位移元的有限单元法 得到各节点的位移值 得到e=Bu, s=De=DBu B=LN
§14.2 对称与反对称边界条件
对称性(Symmetry)
•减小有限元问题规模的一个最有效的方法是利 用对称性。 •如果问题几何、载荷和约束都关于一条线或者 一个平面对称，则问题具有对称性。 •利用对称性时，模型需要修改，对称线或对称 面作为边界并施加约束。
对称性的例子
对称问题的网格划分要点
σ
IV


N1(ξ
IV
)
N2(ξI ) N2 (ξII ) N2 (ξIII ) N2(ξI )
N3(ξI ) N3(ξII ) N3(ξIII ) N3(ξI )
N4(ξI ) σ1
N
4
(ξ
II
)

σ2

N4 (ξIII ) σ3
N
4
(ξ
I
)

σ4
Ve
CNidV
0,
i 1, 2,..., n
<5> 利用应力泛函的单元内磨平：
应力在单元内磨平，减少工作量。
原全域应力场A(σ,σ) 1 M
T
(σ σ) C(σ σ)dV
2 Ve e1
现单元应力场
Ae (σ,σ)

1 2
T
(σ σ) C(σ σ)dV
Ve
单元足够小时，令C=I，有

求逆
ξI [ 3 / 3, 3 / 3], ξII [ 3 / 3, 3 / 3], ...
σ1 a b c b σI

σ2


b
a
b
c


σ II

σ3 c b a b σIII
σ4

b
c
b
a

σ
•在单元内部，不满足平衡方程。这是由于 有限元方程 KU = F 不等价于微元体的平 衡方程。
•在单元内部满足协调性条件 ，这是由于 有限元假定的单元位移函数是连续、单值 的。
应力分析结果验证
进行变形、支反力、应力的初步计算。 查看变形形状和主应力向量。 查看等值线的变化。 应力通常依赖于形状和载荷，而与材料
三种加密网格的方法
1. h-refinement (改变单元尺寸大小) 2. p-refinement (使用插值次数较高的高阶单元) 3. r-refinement (移动节点位置)
h-adaptivity
p-adaptivity
•• • • • •
•
•
••
•
• •
•
• •
•
•
•• ••
• • •
计算力学 (力学系本科生)
Chapter 14
划分网格的一般原则及计算精度 分析
General Meshing Guidelines and Accuracy
§14.1 划分有限元网格时需要注意的问题 有限元方法的两大核心
分析过程的有效性
计算结果的可靠性
有限元模型 恰当的分析方案
计算方法的选择
Let d [u5 v5 5 u6 v6 6 ]T
The objective is to find the transformation matrix T such that d T d,
with d is the new dofs of frame element
d=[ u2 v2 u3 v3 u6 v6 6 ]T
• •
• ••
•
•
•
•
• • • • Eight-node elt.
§14.4 不同类型单元之间的连接
•一般来说，具有不同DOF的不同类型单元不能 连接而共享同一自由度。 (例如不能把三维梁 单元与平面应力单元连接)
•有时候允许某些不同类型的单元连接在同一节 点。
•例：平面Q4与杆单元之间的连接