材料成形的数值分析方法

1 应力平衡方程 作用在等厚度 t 的平面单元体上的应力和单位体积上的体力如图 2.1 所示：
图 2.1 作用在平面单元体上的力
这里逗号表示偏微分，例如： σ x , x = ∂σ x / ∂x 。下面推导单元体的平衡方程：
ΣX = 0 : −σ xtdy − τ xy tdx + (σ x + σ x , x dx)tdy + (τ xy + τ xy , y dy )tdx + f x tdxdy = 0
图 1.1 解析法
图 1.2 数值法
科学家和工程师寻求另一种求解问题的方法－数值解法。 近几十年来随着计算机技术的 飞速发展，数值分析方法已经成为解决工程实际问题的主要工具。
2 数值分析方法分类 ⑴ 有限差分法。求解步骤是首先将求解域划分为网格，然后在网格的节点上用差分方程近 似微分方程。特别适合求解建立于空间坐标系的流体流动问题，在流体力学领域，至今仍占 据支配地位。弹用于求解域几何形状复杂时，精度降低或求解困难。 ⑵ 第二大类数值分析方法。首先建立和原问题微分方程及边界条件相等效的积分提法，然 后建立近似解法。例如：配点法、最小二乘法、Gelekin 法、里兹法等。 这些方法在不同领域、 不同类型的问题中得到成功应用， 但只限于几何形状规则的问题。 其基本原因是它们都是在整个求解域上假设近似函数。 对于几何形状复杂的问题不能建立合 乎要求的近似函数。
8
4 边界条件 弹性体的全部边界用 S 表示。一部分边界已知外力 T x 、 T y 、 T z ，称为力边界条件， 用 Sσ 表示；一部分边界已知位移 u 、 v 、 w ，称为位移边界条件，用 Su 表示，如图 2.3 所 示。
图 2.3 弹性体上的力边界及位移边界
有： S = Sσ + Su ⑴ 力边界条件。由力学平衡方程，有力边界条件：
对于三维的情况，应力单元体的平衡方程为：
5
⎧σ x , x + τ xy , y + τ xz , z + f x = 0 ⎪ ⎪ ⎨τ xy , x + σ y , y + τ yz , z + f y = 0 ⎪ ⎪ ⎩τ zx , x + τ zy , y + σ z , z + f z = 0
HYPERMESH: 大型有限元前处理软件，有多种 CAD 格式接口并可以和多种 CAE 软件实现数据 共享，适用于零件的网格划分，网格划分后零件的组装。其软件包 Hyperworks 可以集成多 种求解器如 Optistruct、Radioss 等，可进行结构的静、动力、非线性及优化分析。
2 专用有限元分析软件 Ls-Dyna 非线性有限元分析求解器，其中显式分析最为有名。可用于碰撞分析、冲压分析、 固液耦合及爆炸计算等。 Dynaform: 板料成形专用有限元分析前、后处理软件，求解器采用 Ls-Dyna。该软件有美
⎤ 0⎥ ⎥ ∂ ⎥ ⎧u ⎫ ⎨ ⎬ ∂y ⎥ ⎥ ⎩v ⎭ ∂⎥ ⎥ ∂x ⎥ ⎦
⎤ 0⎥ ⎥ 0⎥ ⎥ ⎥ ∂⎥ ⎧u ⎫ ∂z ⎥ ⎪ ⎪ ⎥ ⎨v ⎬ ⎥ ⎪ w⎪ 0 ⎩ ⎥ ⎭ ⎥ ∂⎥ ∂y ⎥ ⎥ ∂⎥ ∂x ⎥ ⎦
T
可推广得到空间问题的应变－位移关系：
0 ∂ ∂y 0 ∂ ∂x ∂ ∂z 0
6
同理：
ε y = v, y
按工程上的定义，工程剪应变为直角的改变量。在小变形条件下， β1 ≈ tg β1 ， β 2 ≈ tg β 2 ，
γ xy = β1 + β 2 =
(u + u
,y
dy ) − u
dy
+
(v + v
,x
dx ) − v
dx
= u, y + v, x
故得平面问题的应变－位移关系：
⎧Tx = T x ⎪ ⎪ ⎨Ty = T y ⎪ ⎪ ⎩Tz = T z
（弹性体的内力和外力平衡）
设边界外法线方向余弦为 nx 、 n y 、 nz ，则边界上弹性体内力可表示为：
⎧Tx = nxσ x + n yτ yx + nzτ zx ⎪ ⎨Ty = nxτ xy + nyσ y + nzτ zy ⎪T = n τ + n τ + n σ x xz y yz z z ⎩ z
国 ETA 软件技术公司开发，可用于板料成形的成形及回弹分析，回弹分析功能强。 Autoform: 由瑞士软件公司开发，板料成形有限元分析专用软件，隐式算法，可用于板料成 形起皱、拉裂预测，回弹精度差。该软件的特色是模面设计功能强大。 Pamstamp: 板料成形专用有限元分析软件，法国 ESI 软件技术公司开发，显式分析，可用于 板料成形起皱、拉裂及回弹分析。回弹分析精度不及 Dynaform。 Deform: 由 SFTC 公司开发，是一套基于工艺模拟系统的有限元分析软件，专门用于分析各 种金属成形过程中的三维流动， 提供集邮价值的工艺分析数据， 及有关成形过程中的材料和 温度流动。
⎡∂ ⎢ ⎧ ε x ⎫ ⎢ ∂x ⎪ ⎪ ⎢ ⎨ε y ⎬ = ⎢ 0 ⎪γ ⎪ ⎢ ⎩ xy ⎭ ⎢ ∂ ⎢ ⎢ ⎣ ∂y
⎡∂ ⎢ ∂x ⎢ ⎢0 ⎧εx ⎫ ⎢ ⎪ε ⎪ ⎢ ⎪ y⎪ ⎢0 ⎪εz ⎪ ⎪ ⎢ ⎪ ⎨ ⎬=⎢ ∂ ⎪γ xy ⎪ ⎢ ⎪γ yz ⎪ ⎢ ∂y ⎪ ⎪ ⎢ ⎪γ zx ⎭ ⎪ ⎢0 ⎩ ⎢ ⎢ ⎢∂ ⎢ ⎣ ∂z
3 发展趋势 工程分析的作用已从单纯工程分析扩展到优化设计，并和计算机辅助设计技术相结合， 发展了一门全新的学科 C3P，即 CAD/CAE/CAM/PLM，出现了专门的 CAE 工程师。
4
第二章 弹性力学基本方程与变分原理
§2.1 弹性力学基本方程
弹性力学的基本方程可分为三大类，即应力平衡方程、几何变形方程和材料的物理方 程。我们先从弹性力学的二维问题入手推导三类基本方程，然后将其推广到三维问题。
图 2.2 平面问题中的变形表达
上图给出了一个平面问题的一般应变场，由位形 O12 转变为位形 O 12 。位移 u、v 是 坐标的函数，增量 u, x dx 与 u 或 v 相比是无穷小量。由定义：
' ' '
εx =
dx + ( u + u, x dx ) − u ⎤ Lo '2' − Lo 2 ⎡ ⎦ − dx = u =⎣ ,x Lo 2 dx
弹性体上的力边界条件可用矩阵形式表示为:
T=T
(在 Sσ 上)
其中：
⎡ nx ⎢ T = nσ ， n = ⎢ 0 ⎢0 ⎣
0 ny 0
0 0 nz
ny nx 0
0 nz ny
9
nz ⎤ ⎥ 0⎥ nx ⎥ ⎦
⑵ 位移边界条件。弹性体上的位移边界条件可表示为：
⎧u = u ⎪ ⎪ ⎨v=v ⎪ ⎪ ⎩w = w
CAE 理论基础及软件应用
讲义
材料工程学院 王金彦
第一篇
线性有限元分析基础
1
第一章 概论
§1.1 有限元的发展历史
1 数值分析方法来源于工程实际 在科学技术及工程实际问题中， 有许许多多的实际问题可以归结为微分方程 （常微分方 程或偏微分方程）及边值问题。 解析方法：方程性质简单、求解域几何形状规则的问题。如 图 1.1 所示。 微分方程及边界条件 数值解法：方程性质及求解域几何形状复杂的问题，不能得 到解析解。如图 1.2 所示。
简化得： 同理由 ΣY = 0 得：
σ x , x + τ xy , y + f x = 0 τ xy , x + σ y , y + f y = 0
⎧ ⎪σ x , x + τ xy , y + f x = 0 ⎨ ⎪ ⎩τ xy , x + σ y , y + f y = 0
故平面应力单元体的平衡方程为：
其中 E 为弹性模量， G 为剪切弹性模量， μ 为泊松比，且有以下关系：
G=
E 2 (1 + μ )
将以上二维平面问题的物理方程写成矩阵形式，有
ε = cσ 或 σ = Dε ，这里： c = D-1 c 称为材料的柔度矩阵，D 称为材料的刚度矩阵。可将二维平面问题的物理方程推广到三维
情况，有：
σx μ ⎧ ⎪ε x = E − E (σ y + σ z ) ⎪ ⎪ε = σ y − μ σ + σ x) ⎪ y E E( z ⎪ ⎪ ε = σ z − μ (σ + σ ) y ⎪ z E E x ⎨ τ ⎪ γ xy = xy ⎪ G ⎪ τ ⎪ γ yz = yz ⎪ G ⎪ τ ⎪ γ zx = zx ⎪ G ⎩
其逆形式为： σ = Dε ，定义拉梅常数：
G=
E (1 − μ ) E Eμ ，λ = ，则： λ + 2G = 2 (1 + μ ) (1 + μ )(1 − 2μ ) (1 + μ )(1 − 2μ )
λ λ 0 0 0⎤ ⎡λ + 2G ⎢ λ + 2G λ 0 0 0⎥ ⎢ ⎥ ⎢ λ + 2G 0 0 0 ⎥ D=⎢ ⎥ G 0 0⎥ ⎢ ⎢ G 0⎥ ⎢ ⎥ G⎦ ⎢ ⎥ ⎣
3
§1.2 常用有限元软件介绍
1 大型通用有限元分析软件 NASTRAN: 适用于大型通用结构静、动力及非线性分析，可以计算结构位移、应变、应力、 固有频率及模态分析、响应分析、谱分析及随机振动分析。NASTRAN 属于 Code(即求解器)， 动力学分析是它的特色，非线性功能相对较弱。 ANSYS: 适用于大型结构静、动力及非线性分析，特别适合多场耦合分析。 ABAQUS: 适用于大型结构静、动力及非线性分析，非线性分析功能强大。有用户程序接口， 可二次开发。 MARC: 适用于大型结构静、动力、非线性及多场耦合分析，非线性功能强大。
插值函数 问题的有限元分析 （单元场函数） 满足收敛要求 单元数目增加 （节点增加） 单元自由度增加（插值函数精度提高） 解的近似程度 收敛于精确解。 节点自由度 单元内任一点值 求解域的值 ；
3 有限元法发展史 ⑴ 有限元法的基本思想由特朗（Courant）在 1943 年提出，他第一次尝试把定义在三角形 区域上的分片连续函数和最小位能原理相结合来求解圣维南扭转问题－变分问题的里兹方 法。 ⑵ 工程上真正利用计算机采用有限元法解决工程实际问题的是特纳（Turner）克拉夫 (Clough)所做的工作。1956 年他们在分析飞机结构时，第一次给出了用三角形单元求解平 面应力问题的正确解答。 三角形单元的单元特性由弹性理论方程确定， 采用直接刚度法推导 单元刚度阵。1960 年，Clough 进一步处理了平面弹性问题，并第一次提出了“有限单元法” 的名称。 ⑶ 有限元法在工程上的应用范围： 弹性力学平面问题 空间问题 板壳问题 几何非线性（大变形） 静力平衡问题 动力学问题 稳定问题 非线性问题 材料非线性 弹塑性、塑性 粘弹性、 粘塑性 粘弹塑性 流体力学 传热学 固体力学 流固耦合 热力耦合 电磁场、多场耦合 复合材料
2
⑶ 有限元法。有限元法的出现是数值分析方法研究领域内重大突破性进展。 a: 可以近似几何形状复杂的求解域。基本思想是将连续的求解域离散为一组有限个且按 一定方式相互联结在一起的单元的组合体。 由于单元能按不同联结方式组合且单元本身又可 以有不同形状，故可以近似几何形状复杂的区域。 b: 在每一个单元内假设近似函数分片地表示全求解域上的未知场函数。单元内的近似函 数由未知场函数或其导数在单元的各个节点数值（自由度）和其插值函数表达。
平衡方程的矩阵形式：
Aσ + f = 0
A 是微分算子，
⎡∂ ⎢ ∂x ⎢ ⎢ A=⎢ 0 ⎢ ⎢ ⎢0 ⎢ ⎣
0 ∂ ∂y 0
0 0 ∂ ∂z
yLeabharlann Baidu
∂ ∂y ∂ ∂x 0
0 ∂ ∂z ∂ ∂y
∂⎤ ∂z ⎥ ⎥ ⎥ 0⎥ ⎥ ∂⎥ ⎥ ∂x ⎥ ⎦
f 是体积力向量，
f =⎡ ⎣fx
T
f
f z⎤ ⎦
2 几何变形方程（应变－位移关系）
用矩阵形式表示为： u = u （在 Su 上） 。
写成矩阵形式： ε = Lu ， L = A
3 物理方程（应力－应变关系） 由广义胡克定律，有二维平面应力情况下的物理方程：
1 ⎧ = ε (σ x − μσ y ) x ⎪ E ⎪ 1 ⎪ ⎨ε y = (σ y − μσ x ) E ⎪ 1 ⎪ γ xy = τ xy ⎪ G ⎩
其逆形式为：
7
E ⎧ ⎪σ x = 1 − μ 2 ( ε x − με y ) ⎪ E ⎪ (ε y − με x ) ⎨σ y = 1− μ2 ⎪ ⎪ τ xy = Gγ xy ⎪ ⎩