固体力学学报

序号	应力	塑性应变
1	200	0.002
2	246	0.0235
3	294	0.0474
4	374	0.0935
5	437	0.1377
6	480	0.18

序号	应力	塑性应变
1	200	0.002
2	246	0.0235
3	294	0.0474
4	374	0.0935
5	437	0.1377
6	480	0.18

序号	位置坐标（x，y）	σ_x（MPa）	σ_y（MPa）	τ_xy（MPa）
1	（64，84）	37.48	36.55	2.54%	52.83	52.26	1.09%	-38.85	-37.19	4.46%
2	（48，40）	108.73	109.20	0.43%	112.08	111.51	0.51%	9.99	9.15	9.18%
3	（36，58）	142.32	142.35	0.02%	132.25	134.07	1.36%	8.80	8.49	3.65%
4	（10，18）	30.19	34.42	12.29%	29.96	34.50	13.16%	-61.82	-60.69	1.86%
5	（36，86）	71.44	77.41	7.71%	84.20	88.04	4.36%	19.87	21.78	8.77%
6	（60，16）	57.57	57.18	0.68%	50.28	52.15	3.59%	20.93	21.04	0.52%
7	（64，28）	102.91	102.22	0.68%	107.76	105.78	1.87%	26.19	26.78	2.20%
8	（70，22）	87.70	90.86	3.48%	84.99	88.68	4.16%	33.61	35.75	5.99%
9	（80，84）	17.20	16.87	1.96%	27.40	27.30	0.37%	-51.14	-48.52	5.40%
10	（78，58）	47.24	47.56	0.67%	56.02	57.63	2.79%	-19.19	-18.96	1.21%

序号	位置坐标（x，y）	σ_x（MPa）	σ_y（MPa）	τ_xy（MPa）
1	（64，84）	37.48	36.55	2.54%	52.83	52.26	1.09%	-38.85	-37.19	4.46%
2	（48，40）	108.73	109.20	0.43%	112.08	111.51	0.51%	9.99	9.15	9.18%
3	（36，58）	142.32	142.35	0.02%	132.25	134.07	1.36%	8.80	8.49	3.65%
4	（10，18）	30.19	34.42	12.29%	29.96	34.50	13.16%	-61.82	-60.69	1.86%
5	（36，86）	71.44	77.41	7.71%	84.20	88.04	4.36%	19.87	21.78	8.77%
6	（60，16）	57.57	57.18	0.68%	50.28	52.15	3.59%	20.93	21.04	0.52%
7	（64，28）	102.91	102.22	0.68%	107.76	105.78	1.87%	26.19	26.78	2.20%
8	（70，22）	87.70	90.86	3.48%	84.99	88.68	4.16%	33.61	35.75	5.99%
9	（80，84）	17.20	16.87	1.96%	27.40	27.30	0.37%	-51.14	-48.52	5.40%
10	（78，58）	47.24	47.56	0.67%	56.02	57.63	2.79%	-19.19	-18.96	1.21%

序号	位置坐标（x，y）	w（mm）
1	（64，84）	0.1799	0.1656	0.0143	8.6353%
2	（48，40）	0.3682	0.3728	0.0046	1.2339%
3	（36，58）	0.3723	0.3807	0.0084	2.2065%
4	（10，18）	0.0701	0.0737	0.0036	4.8847%
5	（36，86）	0.1878	0.1838	0.0040	2.1763%
6	（60，16）	0.1990	0.1742	0.0248	14.2365%
7	（64，28）	0.2888	0.2722	0.0166	6.0985%
8	（70，22）	0.2317	0.2049	0.0268	13.0796%
9	（80，84）	0.1229	0.0976	0.0253	25.9221%
10	（78，58）	0.2372	0.2175	0.0197	9.0575%

序号	位置坐标（x，y）	w（mm）
1	（64，84）	0.1799	0.1656	0.0143	8.6353%
2	（48，40）	0.3682	0.3728	0.0046	1.2339%
3	（36，58）	0.3723	0.3807	0.0084	2.2065%
4	（10，18）	0.0701	0.0737	0.0036	4.8847%
5	（36，86）	0.1878	0.1838	0.0040	2.1763%
6	（60，16）	0.1990	0.1742	0.0248	14.2365%
7	（64，28）	0.2888	0.2722	0.0166	6.0985%
8	（70，22）	0.2317	0.2049	0.0268	13.0796%
9	（80，84）	0.1229	0.0976	0.0253	25.9221%
10	（78，58）	0.2372	0.2175	0.0197	9.0575%

序号	位置坐标（x，y）	σ_x（MPa）	σ_y（MPa）	τ_xy（MPa）
1	（20，54）	-93.11	-86.19	8.03%	-52.13	-46.95	11.03%	-25.78	-25.30	1.90%
2	（24，90）	31.04	25.43	22.06%	-12.42	-14.94	16.87%	-62.71	-54.93	14.16%
3	（40，60）	-109.67	-95.19	15.21%	-105.19	-91.81	14.57%	-15.38	-13.55	13.51%
4	（44，66）	-106.57	-92.14	15.66%	-105.85	-96.98	9.15%	-15.33	-14.71	4.21%
5	（52，62）	-87.58	-78.42	11.68%	-81.65	-71.33	14.47%	-11.00	-10.79	1.95%
6	（60，28）	-37.84	-35.41	6.86%	-56.91	-52.86	7.66%	-3.56	-3.31	7.55%
7	（62，52）	-66.44	-60.35	10.09%	-57.73	-49.48	16.67%	-2.17	-2.93	25.94%
8	（70，44）	-64.30	-59.24	8.54%	-49.18	-42.74	15.07%	-3.07	-2.78	10.43%
9	（78，60）	-26.63	-23.61	12.79%	38.54	42.07	8.39%	-20.15	-20.59	2.14%
10	（88，38）	-22.56	-21.09	6.97%	50.90	51.57	1.30%	25.08	24.45	2.58%

序号	位置坐标（x，y）	σ_x（MPa）	σ_y（MPa）	τ_xy（MPa）
1	（20，54）	-93.11	-86.19	8.03%	-52.13	-46.95	11.03%	-25.78	-25.30	1.90%
2	（24，90）	31.04	25.43	22.06%	-12.42	-14.94	16.87%	-62.71	-54.93	14.16%
3	（40，60）	-109.67	-95.19	15.21%	-105.19	-91.81	14.57%	-15.38	-13.55	13.51%
4	（44，66）	-106.57	-92.14	15.66%	-105.85	-96.98	9.15%	-15.33	-14.71	4.21%
5	（52，62）	-87.58	-78.42	11.68%	-81.65	-71.33	14.47%	-11.00	-10.79	1.95%
6	（60，28）	-37.84	-35.41	6.86%	-56.91	-52.86	7.66%	-3.56	-3.31	7.55%
7	（62，52）	-66.44	-60.35	10.09%	-57.73	-49.48	16.67%	-2.17	-2.93	25.94%
8	（70，44）	-64.30	-59.24	8.54%	-49.18	-42.74	15.07%	-3.07	-2.78	10.43%
9	（78，60）	-26.63	-23.61	12.79%	38.54	42.07	8.39%	-20.15	-20.59	2.14%
10	（88，38）	-22.56	-21.09	6.97%	50.90	51.57	1.30%	25.08	24.45	2.58%

基于弹性理论的残余应力反演和变形计算

PDF下载

顾百骏 , 陶祥泽 , 赵颖涛 ^**

固体力学学报 | 研究论文 2024,45(2): 188-200

收起

固体力学学报 | 研究论文 2024, 45(2): 188-200

基于弹性理论的残余应力反演和变形计算

全屏

顾百骏, 陶祥泽, 赵颖涛^**

作者信息

北京理工大学宇航学院力学系，北京，100081

通讯作者:

^** E-mail：ytzhao623@bit.edu.cn.

Residual Stress Inversion and Deformation Calculation Based on Theory of Elasticity

Baijun Gu, Xiangze Tao, Yingtao Zhao^**

Affiliations

Department of Mechanics, School of Astronautics, Beijing Institute of Technology, Beijing, 100081

出版时间: 2024-04-25 doi: 10.19636/j.cnki.cjsm42-1250/o3.2023.036

文章导航

摘要

收起

毛坯在制造过程中，材料力学性能的非均匀性导致其内部产生残余应力. 残余应力会造成结构破坏，在工件的切削去除过程中，残余应力会逐渐释放并引起变形. 采用有限元的“生死单元”技术模拟材料的切削去除过程，并转化为残余应力的释放，基于板壳理论、小变形理论、弹性理论和叠加原理，将径向基函数插值法与几何方程、物理方程相结合，开发出一种新的解析方法，反演残余应力场并计算变形. 结果表明解析方法的理论解与有限元解高度吻合，能够用于判断残余应力分布并计算残余应力弹性释放引起的变形.

关键词

残余应力 / 残余应力反演 / 残余应力释放 / 材料去除 / 有限元模拟

Abstract

收起

During the manufacturing process of a rough workpiece, the non-uniformity of material mechanical properties can result in the generation of residual stresses within the workpiece, which can lead to structural failure. During the cutting removal process of the workpiece, the residual stresses gradually release, causing deformation. In this study, the "birth-death element" technique of finite element analysis was used to simulate the cutting removal process of the material. This process was then transformed into the release of residual stresses. A novel analytical method combining radial basis function interpolation with geometric and physical equations was proposed based on plate shell theory, small-deformation theory, theory of elasticity, and the superposition principle. The method aimed to invert the residual stress field and calculate the resulting deformation. This study was divided into two parts. In the first part, the influence of stress equilibrium equations was neglected, and the radial basis function interpolation method was used to invert the release of residual stresses in thin plates according to the initial residual stress field and the residual stress field after material removal. Next, the stresses were substituted into the physical equations to calculate the strain. The strain was then substituted into the geometric equations, and the plane displacement was calculated by strain integration from geometric equations. Based on the plate shell theory equations, the buckling deformation was inverted according to the plane displacement. In the second part, the reverse process of the first part was performed. Firstly, the buckling deformation caused by material removal in the thin plate was inverted. Then, the buckling deformation was substituted into the plate shell theory equations to invert the release of residual stress and the reconstructed residual stress field. The results demonstrated the reversibility of these two processes. Furthermore, the analytical solutions showed high agreement with the finite element solutions. This suggested that the analytical method proposed in this paper is applicable to thin-plate structures under elastic conditions and can be used to estimate residual stress distribution and predict deformation in the thin-plate cutting removal process.

Key words

residual stress / residual stress inversion / residual stress release / material removal / finite element simulation

引用本文

顾百骏, 陶祥泽, 赵颖涛. 基于弹性理论的残余应力反演和变形计算. 固体力学学报, 2024 , 45 (2) : 188 -200 . DOI: 10.19636/j.cnki.cjsm42-1250/o3.2023.036

Baijun Gu, Xiangze Tao, Yingtao Zhao. Residual Stress Inversion and Deformation Calculation Based on Theory of Elasticity[J]. Chinese Journal of Solid Mechanics, 2024 , 45 (2) : 188 -200 . DOI: 10.19636/j.cnki.cjsm42-1250/o3.2023.036

正文

收起

0　引言

收起

航空整体构件是一种材料去除率高达90%-95%的大型薄壁类零件，在材料去除过程中，铣削应力不断引入，初始残余应力逐渐释放，最终导致工件变形^[1]. 但是，化学铣削实验表明在高速加工条件下毛坯初始残余应力释放对变形的影响占到90%以上，是需要考虑的首要因素^[2]. 当加工后留存的残余应力与外载引起的内应力相互叠加，超过材料的屈服极限时，工件将产生二次变形，结构的承载能力和稳定性会降低，部分脆性材料甚至会发生断裂^[3]. 因此，反演残余应力释放引起的变形，推断末状态的残余应力分布，并由此选择性地消除残余应力，成为工程上的一个重要问题.

1　残余应力反演方法

收起

目前残余应力反演主要分为数值和实验两类. 数值方法的特点是根据部分残余应力测量结果，运用插值或拟合方法推断残余应力分布. 其中径向基函数插值法^[4]函数阶数较低，节点配置灵活且计算工作量小，能够拟合空间中大部分曲面且在边界上不会产生龙格现象，插值精度较高. 反距离加权插值法（idw插值法）将距离的倒数做为选择权重的参量，适用于样本点较多且数值连续性较好的情况^[5]. 克里金插值法是idw插值法的改进版，考虑了空间相关性的问题，插值效果更好^[6].

实验上主要通过测量物理参数变化，如位移、应变等反演残余应力. 声表面波技术反演残余应力^[7]的原理是当损伤层中存在残余压应力时，声表面波相速度会降低，声波频率固定时，残余压应力越高，波速下降幅度越大，得到波速变化量与残余应力大小的关系，反演残余应力.

推广Prime提出的方法^[8]是将钢板进行一次，多次或正交多次切割，使钢板变形，释放残余应力，由于残余应力远小于屈服极限，所以认为残余应力的释放和再分布是弹性过程. 该方法测量残余应力释放产生的位移并将其做为外载施加在无应力无变形的同种钢板上，使钢板变形，得到切割位移引起的残余应力释放量，利用叠加原理将该应力场与末状态残余应力场叠加，得到初始残余应力场.

广义残余应力反演方法^[9]给定一组从数字光电技术（DIC）中测量的位移，将残余应力反演问题表述为一个偏微分方程约束的最优化问题，该方法的优点是能够应用于几乎任意几何形状构件的残余应力反演.

基于声弹性理论的残余应力反演方法^[10]通过推导附加无限小弹性变形且含有预应力材料的一般增量本构关系，根据残余应力的均匀性假设，使用超声方法评估材料的残余应力. 其优点是不需要对残余应力的来源或材料存在的状态进行任何假设.

实验方法还包括从表面位移反演残余应力，通过表面位移数据评估材料子域中的塑性损伤区域积累引起的残余应力，并将其描述为一个将表面位移与未知塑性应变联系起来的积分方程组，除此之外，通过各向同性材料中本征应变的分布也可以模拟结构中存在的裂纹^[11].

2　残余应力反演与变形预测

收起

本文的解析方法是数值方法与实验方法的结合，在材料去除这一特定的加工背景下产生. 首先，部分数值方法反演残余应力的精度较高，但是没有引入弹性力学方程，计算得到的应变不满足应变协调条件，无法预测变形. 而解析方法采用插值精度较高的径向基函数反演残余应力，将材料去除过程等效为解除切削部分对剩余部分的作用力，考虑到毛坯的残余应力远小于屈服极限，所以视残余应力释放过程为毛坯发生弹性变形的过程，因此将板壳理论和弹性力学方程引入径向基函数插值中，使应变协调，克服单纯数值反演不能计算变形的局限性^[12].

现阶段，残余应力的有损检测方法耗时且会对材料造成破坏，无损检测费用高昂，所以残余应力反演是减少测量并获取残余应力分布的一种方案^[13]. 解析方法借鉴了实验方法的部分优点，先测量并反演残余应力释放引起的位移，后代入几何方程计算应变，再根据物理方程求出残余应力释放量，由初始残余应力场减去残余应力释放量，得到重构后的残余应力场，可以免去测量过程，节省成本.

2.1　生成残余应力场

本文的研究对象为薄板结构. 首先利用有限元软件建立一个薄板模型，长宽高分别设置为100 mm、100 mm、5 mm，弹性模量和泊松比设置为E=210 GPa，ν=0.3，屈服极限为200 MPa，按照Drucher公设，将总应变分解为弹性应变和塑性应变，在软件中输入屈服应力和对应的塑性应变，如表1所示. 采用Mises屈服准则，并通过周期加载使薄板产生不均匀的塑性变形，生成残余应力场；即：在薄板上表面选中部分区域沿厚度方向施加正弦周期载荷，如图1所示，连续加载2个周期，在加载过程中薄板部分区域的内应力会大于屈服极限，产生塑性变形，当外载撤去后，发生塑性变形的区域无法恢复，就会对剩余部分产生拉伸或压缩作用，满足自平衡条件的初始残余应力场由此生成^[14]. 由于加载形式随机，初始残余应力场也因此变得复杂，这样可以避免偶然性，更好地验证本文方法的反演效果.

有限元软件的Interaction模块内含有“Model change”功能，即“生死单元”技术，该功能可以使选中的单元不再参与后续的计算过程且对其他单元不再产生任何影响，因此能用于模拟材料的切削去除. 本文采用“Model change”模拟薄板上下两层的切割，并将其转化为残余应力的释放过程.

2.2　反演残余应力释放量并计算变形

图2为薄板结构示意图，黄色部分为区域P，代表加工过程中从薄板上去除的多余材料，上下两层各为1.25 mm，绿色部分为区域Q，代表薄板经材料去除后的成形零件，厚度为2.5 mm，由于是薄板结构，所以不考虑厚度方向（Z方向）的作用力. 材料去除前，区域P对区域Q的作用力可表示为F_Px、F_Py、M_Px、M_Py；区域Q对区域P的作用力表示为F_Qx、F_Qy、M_Qx、M_Qy，它们是作用力与反作用力的关系，大小相等方向相反，如图3所示. F_Px、F_Py对区域Q产生面内拉压作用，M_Px、M_Py产生面外翘曲作用，区域Q保持平衡状态，材料P去除后，F_Px、F_Py、M_Px、M_Py瞬时为0，此时区域Q不平衡，残余应力会弹性释放，产生弹性条件下的应力重构，使Q重新达到平衡状态. 定义区域P去除前区域Q的上表面应力为σ_Q₁（up），下表面应力为σ_Q₁（down），区域P去除后，产生应力重构，定义重构后区域Q的上下表面残余应力分别为σ_Q₂（up），σ_Q₂（down）. 定义区域Q上表面和下表面的应力释放量分别为Δσ_Q（up）和Δσ_Q（down），则：

(1)

应力重构使薄板产生拉压和翘曲变形，可以将应力Δσ_Q（up）和Δσ_Q（down）分解为两部分：

(2)

Δσ_（_F_）由F_Px、F_Py引起，使薄板产生面内拉压变形，Δσ_（_M_）由M_Px、M_Py引起，使薄板产生面外翘曲变形. 由于薄板的翘曲变形量通常远大于拉压变形量，所以认为F_Px、F_Py近似为0，即Δσ_（_F_）=0，本文主要计算由Δσ_（_M_）引起的翘曲变形^[15].

在区域Q的中性面内均匀选取9行9列一共81个样本点，记行数m=9，列数n=9，总共样本点个数N=mn=81，坐标分别记为（x_i，y_i），i=1，2，…，81. 读取材料去除前后，样本点在区域Q上表面对应位置处的3个残余应力分量σ_x、σ_y、τ_xy的应力释放量，分别记为Δσ_x_（up）（x_i，y_i）、Δσ_y_（up）（x_i，y_i）、Δτ_xy_（up）（x_i，y_i），以及下表面对应位置处的应力释放量Δσ_x_（down）（x_i，y_i）、Δσ_y_（down）（x_i，y_i）、Δτ_xy_（down）（x_i，y_i），i=1，2，…，81. 由式（2）可得引起翘曲变形的应力表达式为：

(3)

定义引起薄板翘曲变形的残余应力分量分别为Δσ_x_（_M_）（x_i，y_i），Δσ_y_（_M_）（x_i，y_i），Δτ_xy_（_M_）（x_i，y_i），根据（3）计算可得：

(4)

薄板平面内的残余应力释放满足弹性理论，平面内含有应力σ_x、σ_y、τ_xy，应变ε_x、ε_y、γ_xy和位移u、v一共8个物理量. 这些物理量需满足3个物理方程，3个几何方程和2个应力平衡方程一共8个方程，物理量个数与方程个数相等，方程封闭且有唯一解. 但是在实际测量和反演残余应力的过程中，通常会存在误差，所以在每个样本点位置，平衡方程无法精确满足，反演过程可以不考虑应力平衡方程的约束作用，因此方程个数减为6个，存在2个自由的物理量. 将σ_x和σ_y当做自由物理量，用径向基函数插值法进行反演，而后代入几何方程和物理方程中，此时方程包含剩余的6个物理量，6个方程，方程封闭并求得唯一解，再利用计算的平面位移u和v，结合板壳理论反演翘曲变形，整个过程满足应变协调条件，具体过程如下.

将坐标和应力作为数据，用推广的二维Multiquadric径向基函数^[16]反演残余应力释放的场函数. 以σ_x为例，先说明Multiquadric径向基函数的插值原理. 首先给定选取的N（N>2）组样本点数据：（x_i，y_i，Δσ_x_（_M_）（x_i，y_i）），将每组数据分别代入式（5）中：

(5)

其中λ₁，λ₂，…，λ_N为待定系数，c为参数，会影响插值精度，本文选取c=1. 得到由N个方程构成的方程组，如下所示：

(6)

定义：

(7)

则方程（6）转化为矩阵形式如下：

(8)

矩阵满秩，方程组有唯一解，计算待定系数λ₁，λ₂，…，λ_N，获取σ_x的径向基函数表达式^[17]，即为σ_x的反演结果，记为δσ_x_（_M_）（x，y），同理也可以反演得到δσ_y_（_M_）（x，y），它们是薄板表面正应力释放的场函数.

材料去除引起的残余应力释放为弹性过程，根据薄板弹性理论，δσ_x_（_M_）（x，y）和δσ_y_（_M_）（x，y）沿薄板厚度方向呈线性分布且关于中性面对称^[18]. 所以薄板任意厚度的正应力释放场函数可表示为：

(9)

其中h为区域Q的厚度，z为区域Q沿厚度方向的坐标.

接下来反演τ_xy应力场并计算翘曲变形，将δσ_x_（_M_）（x，y）和δσ_y_（_M_）（x，y）代入物理方程中，得到薄板残余应力释放产生的平面内正应变，记为δε_x_（_M_）（x，y）和δε_y_（_M_）（x，y）：

(10)

定义残余应力释放引起的薄板表面沿坐标系x和y方向的位移分别为δu_（_M_）（x，y）和δv_（_M_）（x，y），应变δε_x_（_M_）（x，y）和δε_y_（_M_）（x，y）与位移δu_（_M_）（x，y）和δv_（_M_）（x，y）的关系为：

(11)

积分得到位移δu_（_M_）（x，y）和δv_（_M_）（x，y），表达式如下：

(12)

其中A（y）和B（x）为待定函数. 定义残余应力释放产生的表面剪应变为δγ_xy_（_M_）（x，y），δγ_xy_（_M_）（x，y）与δu_（_M_）（x，y）和δv_（_M_）（x，y）的关系根据几何方程表示如下：

(13)

将式（12）带入式（13）中，得到剪应变的函数表达式：

(14)

下一步需要确定A′（y）和B′（x）.

将A′（y）和B′（x）设定为待定系数的一维径向基函数，表达式如下：

(15)

其中Y_j为m行样本点每一行对应的y坐标值，X_j为n列样本点每一列对应的x坐标值，m=n=9. 定义薄板表面剪应力释放的场函数为δτ_xy_（_M_）（x，y），与剪应变的关系如（16）所示：

(16)

G为剪切模量，将式（14）和（15）代入式（16）中，得到：

(17)

方程（17）中一共有m+n=18个待定系数，这是经上述方法反演得到的含有待定系数的剪应力释放场函数，将每个样本点坐标（x_i，y_i），i=1，2，…，81代入式（17）中，得到每个样本点含有待定系数的剪应力反演值δτ_xy_（_M_）（x_i，y_i）. 采用让反演得到的每个样本点的δτ_xy_（_M_）（x_i，y_i）和式（4）计算的每个样本点真实的剪应力Δτ_xy_（_M_）（x_i，y_i）一一对应的方式确定待定系数，即：

(18)

将式（17）代入式（18）中得：

(19)

其中i=1，2，…，81，令：

(20)

将式（20）代入式（19）并转化为矩阵形式：

(21)

展开得：

(22)

其中：

(23)

式（22）显示一共有m+n=18个待定系数，N=mn=81个方程，因此（22）是一个超定方程组，根据超定方程组的最小二乘解可以得到最佳系数组合^[19]，计算方法如（24）所示：

(24)

式（24）得到的唯一解即为最小二乘结果，将系数代入式（15）中确定A′（y）和B′（x）的函数表达式，再代入式（17）即可得到δτ_xy_（_M_）（x，y），它是反演得到的薄板表面剪应力释放的场函数，主要引起扭转变形.

将δσ_x_（_M_）（x，y）、δσ_y_（_M_）（x，y）、δτ_xy_（_M_）（x，y）绘制成应力云图并与有限元软件直接计算的结果比较，如图4所示.

将A′（y）和B′（x）积分后代入式（12）中，根据边界条件确定位移的平动项和转动项，得到薄板表面位移函数δu_（_M_）（x，y）和δv_（_M_）（x，y）的表达式. 定义薄板残余应力释放产生的平面外翘曲位移为δw（x，y），与平面内薄板位移δu_（_M_）（x，y）和δv_（_M_）（x，y）的关系：

(25)

取z=h/2，积分得到δw（x，y）：

(26)

将式（26）代入式（25）中，得到：

(27)

同样将c′（y）设定为一维径向基函数的形式：

(28)

将样本点的坐标代入式（27）和（28）：

(29)

式（29）同样利用超定方程组的最小二乘解确定c′（y）系数，代入式（28）得到c′（y）表达式，而后积分并代入式（26），根据位移边界条件得到反演的翘曲函数δw（x，y）. 将δw（x，y）绘制成云图并与真实位移比较，如图5所示：

为了精确地分析反演误差，按照图1所示坐标系，在薄板区域Q的中性面上，靠近四周边界位置选取部分测点，在薄板内部再选取一些测点，同时要避开插值时选取的样本点，避免结果的偶然性，一共选取了10个测点，序号及位置信息如图6（a）所示，读取坐标并计算误差，如表2、表3所示.

σ_x和σ_y的反演误差是由径向基函数插值产生的，而图4和表2显示τ_xy的反演误差稍大于σ_x和σ_y，因为τ_xy是将σ_x和σ_y代入物理方程、几何方程中，经过求解超定方程组最小二乘解得到的，在这个过程中，误差应会逐级放大，且最小二乘拟合不能保证所有的样本点都在拟合曲面上，相比于插值能够保证曲面完全经过样本点，自然会放大误差.

图5显示w的整体反演效果较好，但分析表3发现，远离薄板边界的位置反演误差更小，而靠近边界，误差会逐渐放大. 一方面原因是插值在边界位置的误差会大于内部，另一方面，靠近边界位置薄板的位移绝对值较小，相同的绝对误差，会造成相对误差放大化，例如表3中序号6和序号9测点的对比. 因此，对于薄板位移这种绝对值相对较小的物理量反演，要结合绝对误差和相对误差综合分析反演效果.

对于本文的反演方法，可以通过增加插值所需样本点的密度以及合理选择样本点的分布范围减小误差，保证插值样本点的排列具有一定规律，密度合理，且分布在薄板表面的各个位置，尤其是边界处的样本点数量要充足，才能尽可能减小边界位置的反演误差. 但是，样本点个数过多，径向基函数插值会变得复杂，计算速度减慢，所以要根据实际情况，选取最佳的方案.

2.3　由翘曲变形反演重构后的残余应力场

2.2节描述了基于薄板理论由初始残余应力场和重构后的残余应力场，反演薄板残余应力释放量和翘曲变形的方法. 在2.3节中，将计算2.2节的逆过程，以证明本文方法的合理性，即：已知薄板的翘曲位移w，结合相关理论反演残余应力释放量并计算重构后的残余应力场.

首先测量薄板材料去除后，区域Q样本点产生的翘曲位移Δw（x_i，y_i），用径向基函数插值法对位移Δw（x_i，y_i）进行插值，得到翘曲函数δw（x，y），计算翘曲产生的平面内位移：

(30)

将δu_（_M_）（x，y）和δv_（_M_）（x，y）代入几何方程中，得到残余应力释放引起的平面内应变，记为δε_x_（_M_）（x，y），δε_y_（_M_）（x，y）和δγ_xy_（_M_）（x，y）.

(31)

再将式（31）代入物理方程中得到残余应力释放的场函数：

(32)

用材料去除前的初始残余应力场减去式（32）计算的残余应力释放量，反演重构后的残余应力场，绘制成应力云图的形式并与有限元软件的计算结果比较，如图7所示.

选取10个测点进行误差分析，序号及位置分布如图6（b）所示，计算结果如表4所示.

表4表明材料去除后，根据位移反演末状态残余应力的误差基本在20%以内. σ_x、σ_y和τ_xy都是由径向基函数插值的翘曲位移代入几何方程和物理方程计算得到，所以它们的误差基本处于同一数量级，且与2.2节结论相似，绝对值较小的应力，相对误差会放大化.

本节方法的意义是测量材料去除前的残余应力场和残余应力释放产生的翘曲位移，根据板壳理论反演材料去除后薄板的残余应力，免去了重构后残余应力的测量过程，节省成本.

3　结论

收起

（1）本研究利用有限元仿真代替实验，通过周期加载使薄板产生不均匀的塑性变形并生成初始残余应力场，采用“生死单元”技术模拟材料的切削去除并将其转化为残余应力的释放，整个过程均基于薄板理论.

（2）根据薄板材料去除前后样本点的残余应力，将径向基函数插值与板壳理论、弹性力学方程相结合，反演材料去除引起的的残余应力释放量，同时能够在满足应变协调条件的基础上，计算薄板产生的位移.

（3）参考实验方法，测量材料去除后剩余部分产生的变形量，用径向基函数插值得到位移函数，将位移函数代入几何方程、物理方程，得到残余应力释放的场函数，用初始残余应力场减去残余应力释放量，反演薄板材料去除，应力重构后的残余应力场，且误差在合理的范围内.

参考文献

收起

文献

收起

参考文献引证文献

排序方式：

[1]

Chantzis

, Van-der-Veen

, Zettler

. An industrial workflow to minimise part distortion for machining of large monolithic components in aerospace industry[J]. Procedia CIRP, 2013, 8: 281-286.

[2]

Huang

X M

, Sun

, Li

J F

. Finite element simulation and experimental investigation on the residual stress-related monolithic component deformation[J]. International Journal of Advanced Manufacturing Technology, 2015, 77(5-8): 1035.

[3]

高玉魁. 残余应力基础理论及应用[M]. 上海: 上海科学技术出版社, 2019.

(Gao

Y K

. Theory and Application of Residual Stesses[M]. Shanghai: Shanghai Scientific & Technical Publishers, 2019. (in Chinese))

[4]

马利敏. 径向基函数逼近中的若干理论、方法及其应用[D]. 复旦大学, 2009.

[5]

靳国栋, 刘衍聪, 牛文杰. 距离加权反比插值法和克里金插值法的比较[J]. 长春工业大学学报(自然科学版), 2003, 03: 53-57.

(Jin

G D

, Liu

Y C

, Niu

W J

. Comparison between inverse distance weighting method and Kriging[J]. Journal of Changchun University of Technonlgy, 2003, 03: 53-57. (in Chinese))

[6]

李海涛, 邵泽东. 空间插值分析算法综述[J]. 计算机系统应用, 2019, 28(07): 1-8.

(Li

H T

, Shao

Z D

. Review of spatial interpolation analysis algorithm[J]. Computer Systems & Applications, 2019, 28(7): 1-8. (in Chinese))

[7]

Liu

Z W

, Lin

, Liang

X H

, Du

A Y

. Inversion of surface damage and residual stress in ground silicon wafers by laser surface acoustic wave technology[J]. Ultrasonics, 2021, 113: 1-9.

[8]

李建军, 陈少华. 钢板轧后控冷模拟及残余应力反演研究[J]. 机械强度, 2011, 33(04): 601-606.

(Li

J J

, Chen

S H

. Cooling simulation and inversion of residual stresses of steel plate[J]. Journal Mechanical Strength, 2011, 33(04): 601-606. (in Chinese))

[9]

Chen

M J

, Aquino

, Walsh

T F

. A generalized stress inversion approach with application to residual stress estimation[J]. Journal of Applied Mechanics, 2020: 1-39.

[10]

Kuo

M K

, Lee

H T

. Inversion of residual stress[J]. Journal of Mechanics, 2001, 17(2): 10108.

[11]

Gao

Z J

, Mura

. On the inversion of residual stresses from surface displacements[J]. Journal of Applied Mechanics, 1989, 56(3): 1-6.

[12]

秦国华, 郭翊翔, 王华敏, 等. 飞机整体结构件的“加工变形-疲劳寿命”多目标结构优化方法[J]. 工程力学, 2021, 38(8): 222-236.

(Qin

G H

, Guo

Y X

, Wang

H M

. Multiple objective structural optimization on “machining deformation-fatigue life” of aeronautical monolithic components[J]. Engineering Mechanics, 2021, 38(8): 222-236. (in Chinese))

[13]

Abvabi

, Rolfe

, Hodgson

P D

. An inverse routine to predict residual stress in sheet material[J]. Materials Science and Engineering: A, 2016, 652: 99-104.

[14]

Withers

P J

. Residual stress and its role in failure[J]. Reports on Progress in Physics, 2007, 70(12): 2211.

[15]

叶海潮, 秦国华, 林锋, 等. 面向7075铝厚板加工变形的初始残余应力释放与再分布模型[J]. 稀有金属材料与工程, 2019, 48(01): 123-132.

(Ye

H C

, Qin

G H

, Lin

. Model of release and redistribution of initial residual stresses for machining deformations of 7075 aluminum alloy thick plate[J]. Rare Metal Materials and Engineering, 2019, 48(01): 123-132. (in Chinese))

[16]

Sekii

. Two-dimensional inversion of rotational splitting data[J]. Progress of Seismology of the Sun and Stars, 1990, 367: 1-4.

[17]

蒋伟. 航空毛坯件残余应力场预测[D]. 南京: 南京航空航天大学, 2008: 1-61.

(Jiang

. Measurement and Calculation of Residual Stress Field in Aeronautic Blank[D]. Najing: Nanjing Uinversity of Aeronautics and Astronautics, 2008: 1-61. (in Chinese))

[18]

徐芝纶. 弹性力学[M]. 北京: 高等教育出版社, 2016: 1-10.

(Xu

Z L

. Elasticity[M]. Beijing: Higher Education Press, 2016: 1-10. (in Chinese))

[19]

丁丽娟, 程杞元. 数值计算方法[M]. 北京: 高等教育出版社, 2011: 182-184.

(Ding

L J

, Cheng

Q Y

. Numerical Calculation Method[M]. Beijing: Higher Education Press, 2011: 182-184. (in Chinese))

2024年第45卷第2期

PDF下载

引用本文

BibTeX

文章信息

doi: 10.19636/j.cnki.cjsm42-1250/o3.2023.036

接收时间：2023-07-26
首发时间：2026-04-01
出版时间：2024-04-25

补充材料

相关文章

文章信息

作者

出版历史

收稿日期：2023-07-26

基金

作者信息

北京理工大学宇航学院力学系，北京，100081

通讯作者:

^** E-mail：ytzhao623@bit.edu.cn.

参考文献

分享链接

https://castjournals.cast.org.cn/joweb/gtlxxb/CN/10.19636/j.cnki.cjsm42-1250/o3.2023.036

分享至

全文二维码

扫描看全文

引用本文

BibTeX

本文的引用情况

2种不同金属材料的力学参数

科 Family	属数 Number of genus	种数 Number of species	占总种数比例 Percentage of total species (%)	属 Genus	种数 Number of species	占总种数比例 Percentage of total species (%)
鹅膏菌科Amanitaceae	2	11	5.26	鹅膏菌属 Amanita	10	4.78
小菇科 Mycenaceae	2	12	5.74	丝盖伞属 Inocybe	5	2.39
多孔菌科 Polyporaceae	8	14	6.70	蜡蘑属 Laccaria	5	2.39
红菇科 Russulaceae	3	23	11.00	小皮伞属 Marasmius	6	2.87
				小菇属 Mycena	11	5.26
				光柄菇属 Pluteus	5	2.39
				红菇属 Russula	17	8.13
				栓菌属 Trametes	5	2.39

关闭全屏

BibTeX
EndNote
RefWorks
TxT

序号	应力	塑性应变
1	200	0.002
2	246	0.0235
3	294	0.0474
4	374	0.0935
5	437	0.1377
6	480	0.18

序号

应力

塑性应变

200

0.002

246

0.0235

294

0.0474

374

0.0935

437

0.1377

480

0.18

序号	位置坐标（x，y）	σ_x（MPa）	σ_y（MPa）	τ_xy（MPa）
1	（64，84）	37.48	36.55	2.54%	52.83	52.26	1.09%	-38.85	-37.19	4.46%
2	（48，40）	108.73	109.20	0.43%	112.08	111.51	0.51%	9.99	9.15	9.18%
3	（36，58）	142.32	142.35	0.02%	132.25	134.07	1.36%	8.80	8.49	3.65%
4	（10，18）	30.19	34.42	12.29%	29.96	34.50	13.16%	-61.82	-60.69	1.86%
5	（36，86）	71.44	77.41	7.71%	84.20	88.04	4.36%	19.87	21.78	8.77%
6	（60，16）	57.57	57.18	0.68%	50.28	52.15	3.59%	20.93	21.04	0.52%
7	（64，28）	102.91	102.22	0.68%	107.76	105.78	1.87%	26.19	26.78	2.20%
8	（70，22）	87.70	90.86	3.48%	84.99	88.68	4.16%	33.61	35.75	5.99%
9	（80，84）	17.20	16.87	1.96%	27.40	27.30	0.37%	-51.14	-48.52	5.40%
10	（78，58）	47.24	47.56	0.67%	56.02	57.63	2.79%	-19.19	-18.96	1.21%

序号

位置坐标（x，y）

σ_x（MPa）

σ_y（MPa）

τ_xy（MPa）

反演值

真实值

误差

反演值

真实值

误差

反演值

真实值

误差

（64，84）

37.48

36.55

2.54%

52.83

52.26

1.09%

-38.85

-37.19

4.46%

（48，40）

108.73

109.20

0.43%

112.08

111.51

0.51%

9.99

9.15

9.18%

（36，58）

142.32

142.35

0.02%

132.25

134.07

1.36%

8.80

8.49

3.65%

（10，18）

30.19

34.42

12.29%

29.96

34.50

13.16%

-61.82

-60.69

1.86%

（36，86）

71.44

77.41

7.71%

84.20

88.04

4.36%

19.87

21.78

8.77%

（60，16）

57.57

57.18

0.68%

50.28

52.15

3.59%

20.93

21.04

0.52%

（64，28）

102.91

102.22

0.68%

107.76

105.78

1.87%

26.19

26.78

2.20%

（70，22）

87.70

90.86

3.48%

84.99

88.68

4.16%

33.61

35.75

5.99%

（80，84）

17.20

16.87

1.96%

27.40

27.30

0.37%

-51.14

-48.52

5.40%

（78，58）

47.24

47.56

0.67%

56.02

57.63

2.79%

-19.19

-18.96

1.21%

序号	位置坐标（x，y）	w（mm）
1	（64，84）	0.1799	0.1656	0.0143	8.6353%
2	（48，40）	0.3682	0.3728	0.0046	1.2339%
3	（36，58）	0.3723	0.3807	0.0084	2.2065%
4	（10，18）	0.0701	0.0737	0.0036	4.8847%
5	（36，86）	0.1878	0.1838	0.0040	2.1763%
6	（60，16）	0.1990	0.1742	0.0248	14.2365%
7	（64，28）	0.2888	0.2722	0.0166	6.0985%
8	（70，22）	0.2317	0.2049	0.0268	13.0796%
9	（80，84）	0.1229	0.0976	0.0253	25.9221%
10	（78，58）	0.2372	0.2175	0.0197	9.0575%

序号

位置坐标（x，y）

w（mm）

反演值

真实值

绝对误差

相对误差

（64，84）

0.1799

0.1656

0.0143

8.6353%

（48，40）

0.3682

0.3728

0.0046

1.2339%

（36，58）

0.3723

0.3807

0.0084

2.2065%

（10，18）

0.0701

0.0737

0.0036

4.8847%

（36，86）

0.1878

0.1838

0.0040

2.1763%

（60，16）

0.1990

0.1742

0.0248

14.2365%

（64，28）

0.2888

0.2722

0.0166

6.0985%

（70，22）

0.2317

0.2049

0.0268

13.0796%

（80，84）

0.1229

0.0976

0.0253

25.9221%

（78，58）

0.2372

0.2175

0.0197

9.0575%

序号	位置坐标（x，y）	σ_x（MPa）	σ_y（MPa）	τ_xy（MPa）
1	（20，54）	-93.11	-86.19	8.03%	-52.13	-46.95	11.03%	-25.78	-25.30	1.90%
2	（24，90）	31.04	25.43	22.06%	-12.42	-14.94	16.87%	-62.71	-54.93	14.16%
3	（40，60）	-109.67	-95.19	15.21%	-105.19	-91.81	14.57%	-15.38	-13.55	13.51%
4	（44，66）	-106.57	-92.14	15.66%	-105.85	-96.98	9.15%	-15.33	-14.71	4.21%
5	（52，62）	-87.58	-78.42	11.68%	-81.65	-71.33	14.47%	-11.00	-10.79	1.95%
6	（60，28）	-37.84	-35.41	6.86%	-56.91	-52.86	7.66%	-3.56	-3.31	7.55%
7	（62，52）	-66.44	-60.35	10.09%	-57.73	-49.48	16.67%	-2.17	-2.93	25.94%
8	（70，44）	-64.30	-59.24	8.54%	-49.18	-42.74	15.07%	-3.07	-2.78	10.43%
9	（78，60）	-26.63	-23.61	12.79%	38.54	42.07	8.39%	-20.15	-20.59	2.14%
10	（88，38）	-22.56	-21.09	6.97%	50.90	51.57	1.30%	25.08	24.45	2.58%

序号

位置坐标（x，y）

σ_x（MPa）

σ_y（MPa）

τ_xy（MPa）

反演值

真实值

误差

反演值

真实值

误差

反演值

真实值

误差

（20，54）

-93.11

-86.19

8.03%

-52.13

-46.95

11.03%

-25.78

-25.30

1.90%

（24，90）

31.04

25.43

22.06%

-12.42

-14.94

16.87%

-62.71

-54.93

14.16%

（40，60）

-109.67

-95.19

15.21%

-105.19

-91.81

14.57%

-15.38

-13.55

13.51%

（44，66）

-106.57

-92.14

15.66%

-105.85

-96.98

9.15%

-15.33

-14.71

4.21%

（52，62）

-87.58

-78.42

11.68%

-81.65

-71.33

14.47%

-11.00

-10.79

1.95%

（60，28）

-37.84

-35.41

6.86%

-56.91

-52.86

7.66%

-3.56

-3.31

7.55%

（62，52）

-66.44

-60.35

10.09%

-57.73

-49.48

16.67%

-2.17

-2.93

25.94%

（70，44）

-64.30

-59.24

8.54%

-49.18

-42.74

15.07%

-3.07

-2.78

10.43%

（78，60）

-26.63

-23.61

12.79%

38.54

42.07

8.39%

-20.15

-20.59

2.14%

（88，38）

-22.56

-21.09

6.97%

50.90

51.57

1.30%

25.08

24.45

2.58%