CN107038270B - 一种表面加工残余应力场引起的加工变形的计算方法 - Google Patents

Abstract

本发明公开了一种表面加工残余应力场引起的加工变形的计算方法。在本方法中，通过对零件表面加工残余应力场的建模和分析，将表面残余应力场对变形的作用等价于一组表面力和边缘力。通过在有限元模型中施加相应表面力和边缘力，完成表面加工残余应力场引起的变形计算。相对于传统方法，本方法需要的单元数目少，对网格划分没有特殊要求，能处理更大尺度的问题，使得在汽轮机和航空发动机的大型叶片、船用大尺度螺旋桨上的应用成为了可能。

Description

一种表面加工残余应力场引起的加工变形的计算方法

技术领域

本发明涉及机械加工的物理仿真领域，尤其是复杂曲面零件的表面加工变形的计算领域，更具体地，涉及一种表面加工残余应力场引起的加工变形的计算方法。

背景技术

加工变形是机械零件加工过程中的普遍现象。零件加工变形的是由于切削力、切削热、基底残余应力的释放、表面加工残余应力场的产生和机床夹具等的共同作用。通常这种加工变形大小在微米量级，可以忽略不计。但对于大尺度构件和薄壁件，这种加工变形常常在毫米级，接近零件的设计公差，会大大增加废品率和影响零件的加工质量。

由表面加工残余应力场引起的加工变形是加工变形的一种主要形式，有时是加工变形的主要部分。预测由表面加工残余应力场引起的加工变形对控制此变形有重要意义。传统的预测方法包括直接将表面残余应力场施加到零件有限元模型中的“直接映射法”和E.Brinksmeier提出的“源应力法”。其中直接映射法处理曲面问题耗费大量计算资源以至于在较大尺度零件上应用有很大困难；源应力法只适用于平面问题。另外如图1b所示，直接映射法还对网格划分有特殊要求，其表面影响层(加工表面下几十至几百微米内)需划分多层网格，以保证网格内部有足够数量的单元来映射残余应力，通常需要划分数百万个计算单元才能得到足够精确的结果。没有经验的使用者使用直接映射法时，不符合要求的网格(表面影响层内没有足够数量的单元来映射残余应力)可能得到不精确的甚至得到错误的结果。

发明内容

针对现有表面加工残余应力场引起的加工变形的计算方法存在的缺点，本发明旨在提供一种能够克服现有两种方法中分别存在的不适用于曲面问题和耗费大量计算资源、无法处理大尺度曲面零件、对网格划分和使用者经验有较高要求等缺点的新的表面加工残余应力场引起的加工变形的计算方法。

为了达到上述目的，本发明提供了一种由表面加工残余应力场引起的加工变形的计算方法，该方法包括以下步骤：

(1)划分单元网格；

(2)读取有限元数据，获取有限元模型信息，包括加工表面上的有限元节点和加工表面边缘的有限元节点；

(3)对于加工表面上的每个有限元节点，计算等效体积力矢量f^e；

(4)对于加工表面上的每个有限元节点，将等效体积力矢量f^e等效成等效面力矢量T：

其中，H为表面变质层厚度；

(5)对于加工表面上的每个有限元节点，计算简化后的等效面力矢量T的等效节点力F_surface：

其中，T_i是有限元节点i对应的等效边力矢量；

N_i是有限元节点i对应的单元形函数；

N是待求取节点的单元形函数；

(6)对于加工表面边缘的每个有限元节点，计算等效面力矢量t^e；

(7)对于加工表面边缘的每个有限元节点，将等效面力矢量t^e等效成等效线载荷P：

其中，H为表面变质层厚度；

(8)对于加工表面边缘的每个有限元节点，计算等效线载荷P的等效节点力F_edge：

其中，P_i是有限元节点i对应的等效线载荷；

N_i是有限元节点i对应的单元形函数；

N是待求取节点的单元形函数；

L是加工表面边缘选取的边线长度；

(9)计算总的节点力F_total＝F_surface+F_edge；

(10)将F_total施加于有限元模型，计算表面加工残余应力场引起的变形。

进一步地，步骤(3)中计算等效体积力矢量f^e包括以下步骤：

(i)计算曲面第一基本张量a_αβ的导数：

a_αβ为曲面论第一基本张量的协变分量；

r_α,r_β是协变基矢量；

ξ^λ为曲面的参数坐标；

α,β,λ＝1,2，等于1表示曲面的第一个参数坐标方向，等于2表示曲面的第二个参数坐标方向，以此表示物理量的各个方向上的分量；

(ii)计算克里斯托菲尔符号

其中，a^γλ为曲面第一基本张量的逆变分量；

γ＝1,2，等于1表示曲面的第一个参数坐标方向，等于2表示曲面的第二个参数坐标方向；

(iii)计算曲面第二基本张量b_αβ：

其中，n是加工表面外法线方向的单位矢量；

r是位置矢量，r＝x·i+y·j+z·k；

x、y、z为全局坐标系下的坐标；

i、j、k为全局坐标系坐标轴方向的单位矢量；

(Ⅳ)计算

其中，

是残余应力张量的逆变分量；

是残余应力张量的物理分量；

i,j＝1,2,3，等于1表示曲面的第一个参数坐标方向，等于2表示曲面的第二个参数坐标方向，等于3表示曲面的第三个参数坐标方向；

(Ⅴ)计算等效体积力矢量f^e：

其中，r_β是表面的协变基矢量；

n是表面外法线方向的单位矢量；

σ⁰是初始残余应力场张量的实体形式；

是拉普拉斯算子。

进一步地，步骤(6)中计算等效面力矢量t^e具体包括以下步骤：

(i)计算外法线方向单位矢量n：

n＝v_edge×n_out

其中，n_out是被加工表面的单位外法线方向矢量；

v_edge是加工表面边的单位矢量，v_edge的方向选择使得v_edge、n_out和n符合右手定则；

(ii)计算等效面积力矢量t^e：

t^e＝σ⁰·n。

与现有技术相比，本发明提出的计算方法将加工变形计算问题变成等效体积力矢量f^e和等效面力矢量t^e作用下的静态弹性变形问题来进行求解，能快速有效的对由表面加工残余应力场引起的加工变形进行预计。并且，该方法用普通有限元计算的单元网格即可，约需划分数万个计算单元即可实现精确计算，大大降低了对使用者划分单元网格的能力水平要求，从而能够克服现有两种方法中分别存在的不适用于曲面问题和耗费大量计算资源、无法处理大尺度曲面零件、对网格划分和使用者经验有较高要求等缺点。其计算结果对于相关领域工程技术人员总结变形规律，提出优化加工策略，具有重要意义。

附图说明

图1a是本发明所采用的计算网格示意图；

图1b是现有技术中直接映射法采用的特殊网格示意图；

图2是本发明中的等效力计算流程示意图；

图3是等效体积力矢量f^e和等效面力矢量t^e分别简化为等效面力矢量T和等效线载荷P的原理示意图；

图4是依据本发明提出的算法，对加工表面上的每个有限元节点求取的等效面力矢量T的示意图；

图5是将等效面力矢量T转化为节点力F_surface的原理示意图；

图6是依据本发明提出的算法，对加工表面边上的每个有限元节点求取的等效线载荷P的示意图；

图7是将等效线载荷P转化为节点力F_edge的原理示意图；

图8是铣削过程通常选取的物理坐标系示意图；

图9是矢量v_edge、n_out和n的关系示意图；

图10是本发明的一个应用实例的螺旋桨零件三维造型示意图；

图11是图10的螺旋桨零件的表面喷丸残余应力层深曲线；

图12是图10的螺旋桨零件有限元模型；

图13是图10的螺旋桨零件等效面力矢量T和等效线载荷P的分布示意图；

图14是图10的螺旋桨零件变形计算结果图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。此外，下面所描述的本发明各个实施方式中所涉及到的技术特征只要彼此之间未构成冲突就可以相互组合。

本发明的基本原理如下：

含有初始残余应力场张量

的弹性体，如果不考虑其他因素引起的变形，其变形的基本方程组为：

其中，i,j,k和l是张量下标，i,j,k,l＝1,2,3时分别表示x,y,z轴；

σ_ij-应力张量；

ε_kl-应变张量；

-本征应变；

f_i-i轴方向的体积力；

E_ijkl-杨氏模量；

-初始残余应力场张量；

u_i-i方向上的位移；

u_j-j方向上的位移；

n_j-加工表面外法线方向单位向量；

V-弹性体的整个求解域；

S^t-应力边界；

S^u-位移边界。

如果初始残余应力场张量

是可导函数，上述方程组产生的变形效果相当于在弹性体上施加等效体积力和等效表面力

产生的变形效果。

的矢量形式为：

其中，

-拉普拉斯算子；

的实体形式。

的矢量形式为：

t^e＝σ⁰·n

其中，

的实体形式；n-n_j的矢量形式。

而表面加工残余应力场在宏观上可以用残余应力层深曲线表示，是可导函数。于是，通过上述转换，加工变形计算问题变成体积力f_i ^e和表面力

作用下的静态弹性变形问题，可以使用有限元软件简单的进行计算。

图1a是本发明方法所采用的计算网格示意图，包含2.6×10⁴个二阶四面体单元，图1b是现有技术中直接映射法采用的特殊网格示意图，包含1.7×10⁶个同类二阶四面体单元，其中A为施加残余应力而专门细化后的表面网格。

按照本发明的原理，结合有限元方法，本发明提出了一种新的表面加工残余应力场引起的变形的计算方法，包括如图2所示的以下步骤：

(1)划分单元网格；

(2)读取有限元数据，获取有限元模型信息，包括加工表面上的有限元节点和加工表面边缘有限元节点；

(3)对于加工表面上的每个有限元节点，计算等效体积力矢量f^e，包括如下步骤：

(i)计算曲面第一基本张量a_αβ的导数：

(ii)计算克里斯托菲尔符号

(iii)计算曲面第二基本张量b_αβ：

(Ⅳ)计算

其中，

是残余应力张量的逆变分量，

是残余应力张量的物理分量；

(Ⅴ)计算等效体积力矢量f^e，根据表面残余应力场特点、张量理论和微分几何相关知识，有如下算法：

其中，r_β是表面的协变基矢量，n是表面外法线方向的单位矢量。

(4)对于加工表面上的每个有限元节点，将等效体积力矢量f^e等效成等效面力矢量

其中H为表面变质层厚度；

(5)对于加工表面上的每个有限元节点，计算简化后的等效面力矢量T的等效节点力

其中T_i是有限元节点i对应的等效边力矢量，N_i是有限元节点i对应的单元形函数，N是待求取节点的单元形函数；

(6)对于加工表面边缘的每个有限元节点，计算等效面力矢量t^e，包括如下步骤；

(i)计算外法线方向单位矢量n，n＝v_edge×n_out，其中n_out是被加工表面的单位外法线方向矢量，v_edge是加工表面边的单位矢量，v_edge的方向选择使得v_edge、n_out和n符合右手定则；

(ii)计算等效面积力矢量t^e＝σ⁰·n。

(7)对于加工表面边缘的每个有限元节点，将等效面力矢量t^e等效成等效线载荷

其中H为表面变质层厚度；

(8)对于加工表面边缘的每个有限元节点，计算等效线载荷P的等效节点力

其中P_i是有限元节点i对应的等效线载荷，其中N_i是有限元节点i对应的单元形函数，N是待求取节点的单元形函数；

(9)计算总的节点力F_total＝F_surface+F_edge；

(10)将F_total施加于有限元模型，计算表面加工残余应力场引起的变形。

【应用实例】示例问题为一大尺寸船用螺旋桨的喷丸变形计算，如图10所示，其中，其中字母I表示喷丸表面，字母J表示零件的固定基准面，在仿真计算时在此面上施加固定约束，模拟零件装夹固定，其底面固定，喷丸面为表面的一部分。喷丸残余应力层深曲线如图11所示，忽略喷丸区域不同位置之间的残余应力变化。按照本发明提出的方法对该桨叶喷丸变形计算具体实施方式如下：

(1)对该桨叶划分网格，如图1a、1b和图12所示(划分网格的原则)。图1a是本发明方法所采用的计算网格，包含2.6×10⁴个二阶四面体单元，图1b是现有技术中直接映射法采用的特殊网格，包含1.7×10⁶个同类二阶四面体单元，其中A为施加残余应力而专门细化后的表面网格；

(2)利用本发明提出的算法，通过步骤1、2、3、4、6、7，为喷丸面上和面边缘上的有限元网格节点分别计算的T和p计算结果如图13所示，其中，字母K表示喷丸表面上的每个有限元节点求取的等效面力矢量T，字母L表示加工表面边上的每个有限元节点求取的等效线载荷P。

在有限元模型中通过步骤5、8、9和10施加该T和p对应的节点力以后，该变形问题变成T和p作用下的线弹性静态变形问题，通过有限元软件进行计算，变形计算结果如图14所示。

按照本发明的原理，结合有限元方法，本发明提出了一种新的表面加工残余应力场引起的变形的计算方法，包括如图2所示的以下步骤：

(1)划分单元网格(用普通有限元计算的单元网格即可，约需划分数万个计算单元即可实现精确计算)。

(2)读取有限元数据，获取有限元模型信息。

(3)对于加工表面上的每个有限元节点，计算等效体积力矢量f^e，包括如下步骤：

(i)计算

其中，a_αβ为曲面论第一基本张量的协变分量；

r_α,r_β是协变基矢量；

ξ^λ为曲面的参数坐标；

在以上几个变量中，α、β和λ为上/下标，α,β,λ＝1,2，等于1表示曲面的第一个参数坐标方向，等于2表示曲面的第二个参数坐标方向，以此表示物理量的各个方向上的分量(后面的α、β和λ的定义和本表达式中相同，也为上/下标)。

(ii)利用

计算克里斯托菲尔符号

其中，a^γλ为曲面第一基本张量的逆变分量；

γ为上标，γ＝1,2，等于1表示曲面的第一个参数坐标方向，等于2表示曲面的第二个参数坐标方向。

(iii)计算曲面第二基本张量b_αβ：

其中，n是加工表面外法线方向的单位矢量；

r是位置矢量，r＝x·i+y·j+z·k；

x、y、z为全局坐标系下的坐标，i、j和k为全局坐标系坐标轴方向的单位矢量；全局坐标系即有限元建模时描述节点位置时使用的三维直角坐标系。

(Ⅳ)计算

其中，

是残余应力张量的逆变分量；

是残余应力张量的物理分量；

i和j为上/下标，i,j＝1,2,3，等于1表示曲面的第一个参数坐标方向，等于2表示曲面的第二个参数坐标方向，等于3表示曲面的第三个参数坐标方向。

需要说明的是，此处残余应力张量的物理分量依照选取的物理坐标系的不同而不同，在有的物理坐标系选择下的表达比较简化，在有的物理坐标系选择下

的表达比较复杂，这会给计算带来难度，但不影响最终结果。

物理坐标系的选取因具体加工方法而有所不同。如对于铣削加工，一般选择如图8所示的物理坐标系，选择进给方向、进给方向的垂向和曲面的外法线方向为物理坐标系的三个坐标轴方向，则残余应力在该坐标系中表示较为简化(非对角项为零)。图8中，其中字母B表示加工表面，字母G表示进给方向，字母H表示选取的物理坐标系。

(Ⅴ)计算等效体积力矢量f^e，根据表面残余应力场特点、张量理论和微分几何相关知识，有如下算法：

(4)如图3所示，对于加工表面节点列表中的每个有限元节点，将等效体积力f^e等效成简化后的等效面力矢量T：

其中H为表面变质层厚度，dh是H的微分。

如图4所示，其中字母B表示加工表面，字母C表示加工表面上的有限元节点，字母D表示加工表面上的每个有限元节点处按照本发明提出的算法所求取的等效面力矢量T。经上述计算后，加工表面上的每个有限元节点都得到了一个简化后的等效面力矢量T。

(5)如图5所示，对于加工表面节点列表中的每个有限元节点，计算简化后的等效面力矢量T的等效节点力F_surface。图5放大显示了图4所示加工表面上的一个三角形单元，每个有限元节点都计算了一个简化后的等效面力矢量T，节点1、2和3的T分别标记为T₁、T₂和T₃，节点1、2和3的形函数分别标记为N₁、N₂和N₃。在整个三角形表面，等效面力矢量T按照形函数进行分布，在表面上的任意一点，简化后的等效面力矢量T为：

其中，ks是该单元在加工表面上的节点总数，对于图5所示情况，ks＝3；

N_i是有限元节点i对应的单元形函数。

依据虚功原理，对于分布面力在某一面上节点上产生的等效节点载荷可以写作：

其中，N是待求取节点的单元形函数，s是节点1、2和3围成的三角形面积，ds是S的微分。

(6)对于加工表面边上的每个有限元节点，计算等效面力矢量t^e，包括如下步骤：

(i)计算外法线方向单位矢量n：

n＝v_edge×n_out

其中，n_out是被加工表面的单位外法线方向矢量；

v_edge是加工表面边切线方向的单位矢量；

v_edge的方向选择使得v_edge、n_out和n符合右手定则，三个矢量的相互关系如图9所示，其中，其中字母B表示加工表面。

(ii)计算等效面积力矢量t^e：

t^e＝σ⁰·n。

(7)如图3所示，对于加工表面边上的每个有限元节点，将等效面力矢量t^e等效成加工表面边上的等效线载荷P：

其中，H为表面变质层厚度，dh是H的微分。

如图6所示，其中字母B表示加工表面，字母E表示加工表面边上的有限元节点，字母F表示加工表面边上的每个有限元节点处按照本发明所提出的方法求取的等效线载荷P。经上述计算后，加工表面边上的每个有限元节点都得到了一个等效线载荷P。

(8)如图7所示，对于加工表面边上的每个有限元节点，计算P的等效节点力。图7放大显示了图6所示加工表面边上的一个三角形单元，其节点1和2在加工表面边上，节点1和2都计算了一个等效线载荷P，分别标记为P₁和P₂，节点1和2的形函数分别标记为N₁和N₂。

在节点1到节点2的边上，等效线载荷P按照形函数进行分布，在边上的任意一点，等效线载荷P可以写为如下形式：

其中，N_i是有限元节点i对应的单元形函数；

P_i是节点i处的等效线载荷P；

kl-本单元加工表面所求边上的节点数目，对于图7所示情况，kl＝2。

依据虚功原理，等效线载荷P在边缘某一节点上产生的等效节点载荷F_edge可以写作如下形式：

其中，N是待求取节点的单元形函数；

L是单元的边长，在图7所示情况中，即节点1到节点2的距离。

(9)计算总的节点力F_total＝F_surface+F_edge。

(10)将F_total施加于有限元模型，并利用有限元软件计算表面加工残余应力场引起的变形。

通过本发明提出的计算方法能快速有效的对由表面加工残余应力场引起的加工变形进行预计，其计算结果对于相关领域工程技术人员总结变形规律，提出优化加工策略，具有重要意义。

本领域的技术人员容易理解，以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (3)

1.一种由表面加工残余应力场引起的加工变形的计算方法，该方法包括以下步骤：

(1)对加工表面划分单元网格；

(2)读取有限元数据，获取有限元模型信息，包括加工表面上的有限元节点和加工表面边缘的有限元节点；

(3)对于加工表面上的每个有限元节点，计算等效体积力矢量f^e：

f^e＝-▽·σ⁰

其中，▽是拉普拉斯算子；σ⁰是初始残余应力场张量

的实体形式；i,j是张量下标，i,j＝1,2,3，等于1时表示x轴，等于2时表示y轴，等于3时表示z轴；

(4)对于加工表面上的每个有限元节点，将等效体积力矢量f^e等效成等效面力矢量T：

其中，H为表面变质层厚度；

(5)对于加工表面上的每个有限元节点，计算简化后的等效面力矢量T的等效节点力F_surface：

其中，T_i是有限元节点i对应的等效面力矢量；

N_i是有限元节点i对应的单元形函数；

N是待求取节点的单元形函数；

ks是有限元节点i所属单元在加工表面上的节点总数；

S是节点1、节点2和节点3围成的三角形面积，ds是S的微分；

(6)对于加工表面边缘的每个有限元节点，计算等效面力矢量t^e：

t^e＝σ⁰·n

其中，n是加工表面外法线方向单位向量n_j的矢量形式，σ⁰是初始残余应力场张量

的实体形式；i,j是张量下标，i,j＝1,2,3，等于1时表示x轴，等于2时表示y轴，等于3时表示z轴；

(7)对于加工表面边缘的每个有限元节点，将等效面力矢量t^e等效成等效线载荷P：

其中，H为表面变质层厚度；

(8)对于加工表面边缘的每个有限元节点，计算等效线载荷P的等效节点力F_edge：

其中，P_i是有限元节点i对应的等效线载荷；

N_i是有限元节点i对应的单元形函数；

N是待求取节点的单元形函数；

L是加工表面边缘选取的边线长度；

kl是当前单元加工表面所求边上的节点数目；

(9)计算总的节点力F_total＝F_surface+F_edge；

(10)将F_total施加于有限元模型，计算表面加工残余应力场引起的变形。

2.如权利要求1所述的一种由表面加工残余应力场引起的加工变形的计算方法，其特征在于，步骤(3)中计算等效体积力矢量f^e包括以下步骤：

(i)计算曲面第一基本张量a_αβ的导数：

a_αβ为曲面第一基本张量的协变分量；

r_α,r_β是协变基矢量；

ξ^λ为曲面的参数坐标；

α,β,λ＝1,2，等于1表示曲面的第一个参数坐标方向，等于2表示曲面的第二个参数坐标方向，以此表示物理量的各个方向上的分量；

(ii)计算克里斯托菲尔符号

其中，a^γλ为曲面第一基本张量的逆变分量；

γ＝1,2，等于1表示曲面的第一个参数坐标方向，等于2表示曲面的第二个参数坐标方向；

(iii)计算曲面第二基本张量b_αβ：

其中，n是加工表面外法线方向的单位矢量；

r是位置矢量，r＝x·i+y·j+z·k；

x、y、z为全局坐标系下的坐标；

i、j、k为全局坐标系坐标轴方向的单位矢量；

(Ⅳ)计算

其中，

是残余应力张量的逆变分量；

是残余应力张量的物理分量；

i,j＝1,2,3，等于1表示曲面的第一个参数坐标方向，等于2表示曲面的第二个参数坐标方向，等于3表示曲面的第三个参数坐标方向；

(Ⅴ)计算等效体积力矢量f^e：

其中，r_β是表面的协变基矢量；

n是表面外法线方向的单位矢量；

σ⁰是初始残余应力场张量的实体形式；

▽是拉普拉斯算子。

3.如权利要求1或2所述的一种由表面加工残余应力场引起的加工变形的计算方法，其特征在于，步骤(6)中计算等效面力矢量t^e具体包括以下步骤：

(i)计算外法线方向单位矢量n：

n＝v_edge×n_out

其中，n_out是被加工表面的单位外法线方向矢量；

v_edge是加工表面边的单位矢量，v_edge的方向选择使得v_edge、n_out和n符合右手定则；

(ii)计算等效面积力矢量t^e：

t^e＝σ⁰·n。

Images