CN108804783A - 一种毛坯初始残余应力的测评方法 - Google Patents

Abstract

本发明公开了一种毛坯初始残余应力的测评方法,根据零件加工变形要求，确定变形观测位置及其数量；根据精度要求设置零件划分层数，并逐层进行影响系数计算；确定调整后毛坯残余应力系数矩阵；基于旋转算法进行调整后的残余应力求解。本发明所涉及的计算模型及过程通用性好、逻辑性强，易于实现自动化，为残余应力引起加工变形的计算提供理论基础，能够基于工件变形要求对毛坯生产提出初始残余应力要求。

Description

一种毛坯初始残余应力的测评方法

技术领域

本发明涉及整体结构件制造与厚板生产领域，具体是指整体结构件制造与厚板生产领域基于单位力法与旋转算法的毛坯初始残余应力评测方法。

背景技术

航空整体结构件有着重量轻、装配时间短、生产效率高和生产成本低的优势，在航空领域被广泛使用。但由于其尺寸大且在加工过程中材料去除率高达90％以上造成刚度较低，毛坯中初始残余应力的存在容易引发加工变形。虽然通过预拉伸工艺能够大幅度地降低铝厚板中的残余应力水平，但是残留部分的残余应力仍能引发变形问题。为了研究航空整体结构件变形问题，可以通过有限元方法和解析式方法进行分析，其中有限元方法为主要的分析手段。

在有限元分析中，单元生死技术是最常用的材料去除仿真方法。Liu等人在使用层削法获得7085-T7452铝合金锻件初始残余应力的基础上，通过ANSYS有限元软件进行了挡风玻璃框架结构件变形分析，将仿真结果与实验结果比较得出有限元方法能够有效且可靠的预测初始残余应力引起的加工变形。Zhang等人在使用轮廓法获得7050-T7451铝合金厚板初始残余应力的基础上，通过有限元方法研究了T形件在不同位置处的加工变形，总结出适当的厚度方向位置可以有效地降低零件变形。通过单元生死技术进行材料去除，杀死的单元在仿真前已经设定，但是被装夹的工件在加工过程中仍存在局部变形，这将导致无法对工件的加工尺寸精度进行预测。为解决该问题，Cerutti等人以层剥法测量得到毛坯中的初始残余应力，通过在FORGE 有限元软件中加入布朗运算，同时实现了对加工整体变形和加工尺寸精度的预测，并对比实验证明了预测方法的正确性。

在解析式分析中，Nervi等人将初始残余应力拟合为空间坐标系中的函数，通过Navier-Lamè方程对7050-T7451铝合金厚板的加工变形进行了预测，并结合实验对比证明了模型的有效性。但是这种解析式分析方法同样需要结合有限元手段，难以在宏观上揭示残余应力引起加工变形的内在规律。目前的纯解析预测方法局限性较大，仅能分析材料去除方式为简单厚板层削和开单槽的弯曲变形情况。

发明内容

本发明的目的是为了进一步拓展解析式的使用范围，实现宏观上的残余应力引起变形分析，对航空整体结构件的变形进行解析式推导，而提出的结合有限元的单位力计算方法。

本发明采用以下技术方案实现上述目的。一种毛坯初始残余应力的测评方法,其步骤如下：

1)根据零件加工变形要求，确定变形观测位置；

观测位置以实际图纸要求出发，将基准位置和变形测量位置设为观测位置，或在分析中通过将基准位置设为不动点，对变形测量位置进行观测；

2)根据精度要求设置零件划分层数，并逐层进行影响系数计算；

使用有限元分析软件，将零件沿厚度方向分为n层，单元类型使用二阶单元；依次对第i 层的轧制方向和横向方向施加大小为1MPa的残余应力，计算得到影响因子D_Zi(x)和D_Zi(y)；则该点的变形值其中Σ_Xi(z)和Σ_Yi(z)(1≤k≤n)分别为轧制方向和横向方向上的等效应力；

3)确定调整后毛坯残余应力系数矩阵；

根据实际的加工要求，对零件额外施加一个调整应力Σ_OPTXi，记∑′_Xi(z)＝∑_Xi(z)+∑_OPTXi(z)为调整后的毛坯残余应力，控制零件变形在2Δ_Z范围内，使得

由于毛坯初始残余应力具有自平衡特性，因此调整应力函数也应具有自平衡特性，即存在关系通过引入一个单位变量z_n+1＝1，结合夹逼原则可得AZ≤0；

式中：系数矩阵N＝2n，待定的毛坯应力向量 系数d_ij、δ_Z分别为D_ij＝D_Zj(x_i)、Δ_Z的归一化值，应力z_i＝Σ′_i＝Σ′_Xi(z)，1≤i≤n；

4)基于旋转算法进行调整后的残余应力求解；

采用旋转算法求解变量Z，假定不等式解的一般形式为Z＝Vν+Wω，其中ν为任意的系数向量，ω为任意的非负系数向量，V、W为基向量矩阵；记q、p分别为A＝[a₁,a₂,...,a_p]的行数和列数，则由A既可生成凸锥A_ω＝{x|x＝ω₁a₁+ω₂a₂+...+ω_pa_p}，也可生成线性空间A_v＝{y|y＝ν₁a₁+ν₂a₂+...+ν_pa_p}(ω_i≥0，1≤i≤p)；再记矩阵B＝{a_i|a_i∈A_π，且-a_i∈A_π}，那么求解基向量矩阵V和W的旋转算法如下：

S01：初始化基向量矩阵：

h＝1；V＝I_q×q，W＝Φ，U＝[V,W](I_q×q为单位矩阵，Φ为空集)；I_j＝Φ(1≤j≤q)；

S02：搜索旋转枢轴；

从t＝U^Ta_h中，找出t_k≠0且k≤c_V；其中r_V、r_W、r_U和c_V、c_W、c_U分别为V、W、U的行数和列数，然后执行步骤S03；否则执行步骤S04；

S03：执行旋转运算：

令枢轴号为pivot＝k，将枢轴列标准化为除枢轴(即第pivot列)外，U中其它所有元素进行如下的枢轴变换，即u_k＝u_k+t_ku_pivot,k≠pivot；如果a_h∈B，将u_pivot从V中剔除，即V＝V-u_pivot如果将u_pivot从V中剔除，添加至W中，即将h 添加到集合I_j(j≠pivot,1≤j≤c_V)中，即I_j＝I_j Y{h}后，执行步骤S05；

S04：执行非旋转运算：

将h添加至所有满足t_k＝0的集合I_j中，即I_j＝I_j Y{h}(1≤j≤c_U)后，令矩阵Z＝{w_k(i,j)＝w_i+w_j|t_i<0,t_j>0,w_i∈W,w_j∈W}，则I_k(i,j)＝(I_iΙI_j)Y{h}；从所有I_k(i,j)中找出满足且的集合I_r(a,b)；其中：w_f是满足t_f＝0的列向量，则得Z的子集Z*＝{w_i|i∈I_r(a,b)}；

然后从W中剔除满足t_j>0的列向量w_j；若a_h∈B，再从W中剔除所有满足的列向量w，并将Z*中所有向量加至W，即W＝{W,Z*}；

S05：判断终止条件：

如果h<p，则h＝h+1；计算I_j(1≤j≤r_U)后转至步骤2；否则结束并返回V和W；

最后，按照旋转算法求出基向量矩阵V和W后，将z_n+1＝1(一维)或z_N+1＝1(二维)代入式Z＝Vν+Wω后，经整理可得调整后的等效残余应力。

本发明所涉及的计算模型及过程通用性好、逻辑性强，易于实现自动化，为残余应力引起加工变形的计算提供理论基础，能够基于工件变形要求对毛坯生产提出初始残余应力要求。

具体实施方式

下面结合实施例对本发明进行详细描述：

本发明的工件的材料选用经过拉伸率为δ＝2.4％的7075-T73铝合金板材，其弹性模量为 E＝70.5Gpa，泊松比为ν＝0.33。毛坯尺寸为长L＝200mm，宽W＝20mm，高H＝6mm，而铣削成形的零件尺寸为长l'＝200mm，宽w'＝20mm，高h'＝4.5mm，要求变形不超过0.1 mm。

步骤1：确定观测位置及数量：

由于零件具有对称性，选择长度方向55mm和100mm位置作为观测点进行分析。

步骤2：零件单元划分：

沿厚度方向Z对工件进行分层,层数n＝20，每层厚度为Δz＝0.3mm。若选取每层中心的残余应力作为该层的应力值，通过静力平衡条件求解出每层的应力值如表1所示。

表1毛坯内每层残余应力值

步骤3：影响系数计算：

有限元分析中单元类型采用CPS8R，网格密度设定为0.15，总共有53320个单元。选择“Static，General”分析步类型后，使用Model change进行“单元杀死”，模拟材料的去除。零件沿厚度方向Z从z＝0至z＝4.5mm共分为15层，考虑距离端面分别为x＝55mm和x＝100mm位置的A、B两点，通过逐层施加单位力，利用有限元方法计算的单位变形(即影响系数)如表2所示。

表2影响系数

步骤4：确定调整后残余应力系数矩阵：

由于薄板件的毛坯初始残余应力分布和结构均呈对称形式，55mm和100mm处两点的变形线性相关，因此需控制100mm处变形不超出Δ_Z＝0.05即可。结合表2的数据，根据式自平衡特性，可知应力调整模型中z₁₆＝1，且系数矩阵A如表3所示。

表3归一化的系数矩阵

步骤5：调整后残余应力求解：

初始化生成元矩阵h＝1。V＝I_16×16，W＝Φ，U＝[V,W]。I_j＝Φ(1≤j≤16)。

由于t＝[-1,-0.8561,-0.7131,-0.5704,-0.4278,-0.2848,-0.1423,0.0003,0.1429, 0.2855,0.4281,0.5707,0.7133,0.8559,0.9985,-0.0008]^T，可见t₁≠0且u₁在V中，故pivot＝1。

执行第一次旋转变换得矩阵如表4所示。然后将h＝1添加到I_j(j≠1)中，即I_j＝{1}(2≤j≤16)。此时将枢轴列u₁从V中剔除并加至W中。

表4第一次旋转运算

由于h＝1，p＝4，h<p，故h＝2。

又，t＝[1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-0.0016]^T，由于t₁₆≠0且 u₁₆列在V中，则枢轴povit＝16。执行第二次旋转变换得矩阵如表5所示。此时将枢轴列 u₁₆从V中剔除并加至W中。

表5第二次旋转运算

由于依然存在h<p，故h＝3。

t＝[0.9992,0.1440,0.2869,0.4296,0.5722,0.7152,0.8578,1.0003,1.1429,1.2855,1.4281,1.5707,1.7133,1.8559,1.9985,-0.5]^T，由于t₂≠0且u₂在V中，则povit ＝2。执行第三次旋转变换得矩阵如表6所示。由于a₃∈B，应将枢轴列u₂从V中剔除。

表6第三次旋转运算

更新生成元矩阵，此时h＝3，p＝4，h<p的关系可进一步令h＝4。

计算得t＝[-1.1430,0,0,0,0,0,0,0,0,0,0,0,0,0,0]^T，无枢轴。由于 t₁<0，t₁₆＝0，故Z＝Φ。因此Z*＝Φ。由于a₄∈B，从W中剔除u₁后得W＝u₁₆。又，因此最终有W＝u₁₆，如表7所示。

表7第四次非旋转运算

至此，h＝4，又p＝4，故基向量矩阵最终为V＝[u₃,u₄,...,u₁₅]和W＝u₁₆。

最后，由于z₁₆＝1，经整理后可得毛坯内的应力X应调整如下后才能保证变形在0.1mm的范围之内，即：

其中v_i(1≤i≤13)为任意数，w≥0。

Claims (1)

1.一种毛坯初始残余应力的测评方法,其步骤如下：

1)根据零件加工变形要求，确定变形观测位置；

观测位置以实际图纸要求出发，将基准位置和变形测量位置设为观测位置，或在分析中通过将基准位置设为不动点，对变形测量位置进行观测；

2)根据精度要求设置零件划分层数，并逐层进行影响系数计算；

使用有限元分析软件，将零件沿厚度方向分为n层，单元类型使用二阶单元；依次对第i层的轧制方向和横向方向施加大小为1MPa的残余应力，计算得到影响因子D_Zi(x)和D_Zi(y)；则该点的变形值其中Σ_Xi(z)和Σ_Yi(z)(1≤k≤n)分别为轧制方向和横向方向上的等效应力；

3)确定调整后毛坯残余应力系数矩阵；

根据实际的加工要求，对零件额外施加一个调整应力Σ_OPTXi，记∑′_Xi(z)＝∑_Xi(z)+∑_OPTXi(z)为调整后的毛坯残余应力，控制零件变形在2Δ_Z范围内，使得

由于毛坯初始残余应力具有自平衡特性，因此调整应力函数也应具有自平衡特性，即存在关系通过引入一个单位变量z_n+1＝1，结合夹逼原则可得AZ≤0；

式中：系数矩阵N＝2n，待定的毛坯应力向量系数d_ij、δ_Z分别为D_ij＝D_Zj(x_i)、Δ_Z的归一化值，应力z_i＝Σ′_i＝Σ′_Xi(z)，1≤i≤n；

4)基于旋转算法进行调整后的残余应力求解；

采用旋转算法求解变量Z，假定不等式解的一般形式为Z＝Vv+Wω，其中ν为任意的系数向量，ω为任意的非负系数向量，V、W为基向量矩阵；记q、p分别为A＝[a₁,a₂,...,a_p]的行数和列数，则由A既可生成凸锥A_ω＝{x|x＝ω₁a₁+ω₂a₂+...+ω_pa_p}，也可生成线性空间Aν＝{y|y＝ν₁a₁+ν₂a₂+...+ν_pa_p}(ωi≥0，1≤i≤p)；再记矩阵B＝{a_i|a_i∈A_π，且-a_i∈A_π}，那么求解基向量矩阵V和W的旋转算法如下：

S01：初始化基向量矩阵：

h＝1；V＝I_q×q，W＝Φ，U＝[V,W](I_q×q为单位矩阵，Φ为空集)；I_j＝Φ(1≤j≤q)；

S02：搜索旋转枢轴；

从t＝U^Ta_h中，找出t_k≠0且k≤c_V；其中r_V、r_W、r_U和c_V、c_W、c_U分别为V、W、U的行数和列数，然后执行步骤S03；否则执行步骤S04；

S03：执行旋转运算：

令枢轴号为pivot＝k，将枢轴列标准化为除枢轴(即第pivot列)外，U中其它所有元素进行如下的枢轴变换，即u_k＝u_k+t_ku_pivot,k≠pivot；如果a_h∈B，将u_pivot从V中剔除，即V＝V-u_pivot如果将u_pivot从V中剔除，添加至W中，即将h添加到集合I_j(j≠pivot,1≤j≤c_V)中，即I_j＝I_jY{h}后，执行步骤S05；

S04：执行非旋转运算：

将h添加至所有满足t_k＝0的集合I_j中，即I_j＝I_jY{h}(1≤j≤c_U)后，令矩阵Z＝{w_k(i,j)＝w_i+w_j|t_i<0,t_j>0,w_i∈W,w_j∈W}，则I_k(i,j)＝(I_iII_j)Y{h}；从所有I_k(i,j)中找出满足且的集合I_r(a,b)；其中：w_f是满足t_f＝0的列向量，则得Z的子集Z*＝{w_i|i∈I_r(a,b)}；

然后从W中剔除满足t_j>0的列向量w_j；若a_h∈B，再从W中剔除所有满足的列向量w，并将Z*中所有向量加至W，即W＝{W,Z*}；

S05：判断终止条件：

如果h<p，则h＝h+1；计算I_j(1≤j≤r_U)后转至步骤2；否则结束并返回V和W；

最后，按照旋转算法求出基向量矩阵V和W后，将z_n+1＝1(一维)或z_N+1＝1(二维)代入式Z＝Vν+Wω后，经整理可得调整后的等效残余应力。