一种基于有限测试点的构件残余应力场的预测方法

摘要

本发明公开了一种基于有限测试点的构件残余应力场的预测方法，包括：1、结构离散化；2、单元节点进行分类，建立位移形函数矩阵、单元应力函数矩阵，然后求解单元矩阵；3、根据单元矩阵，建立应力形质点残余应力向量列阵、整体节点载荷列阵、求解整体矩阵；4、根据整体矩阵的秩确定实验测试点个数、分布及测试；5、根据所测残余应力实验结果，作为残余应力边界条件，求得所有应力形质点残余应力；6、根据应力形质点残余应力和单元应力函数矩阵，求解到非应力形质点残余应力，从而可得构件的残余应力场。本发明可精确预测构件的残余应力，而且不受构件材料结构形状的限制，尤其在未知构件制造工艺以及复杂构件残余应力预测具有明显优势。

说明书

技术领域

本发明属于无损检测技术领域，更具体地，涉及一种基于有限测试点的构件残余应力场的预测方法。

背景技术

物体在不受外力时，内部表现出的保持自相平衡的应力系统，通常被称为残余应力。各种机械加工工艺，例如切削、铸造、热处理、焊接、装配等都会不同程度的产生残余应力，残余应力对材料的静态力学性能、应力腐蚀性能、抗疲劳性能、尺寸稳定性以及使用寿命均有着显著的影响，如何精确预测构件内部残余应力场分布对构件的材料及结构设计具有重要意义，已成为业内人士近些年来普遍关注的焦点问题。

目前，构件残余应力获取的方法通常有实验测试法和有限元仿真法。其中，实验测试法按对构件的破坏性有无损检测法和有损检测法两类，X射线法等无损检测法虽然可以检测构件残余应力，但只能检测构件表面或近表面位置残余应力分布；打孔法、剥层法等有损检测法可以检测构件内部一定深度位置残余应力分布，但由于对构件的破坏性，使其在工程实际应用中都受限。有限元仿真法是另一种获知构件残余应力场分布的方法，但目前的研究大都是从构件实际制造过程例如铸造、切削、焊接等方面去模拟得到其分布，由于毛坯铸件自然放置后，残余应力相应发生了变化，仅仅以铸造过程进行模拟无法对残余应力场进行表征，同时，由于缺乏对毛坯初始应力场的正确的表征方法，致使预测结果与实际结果误差较大，同时对于构件制造过程复杂或者不知制造工艺的构件其残余应力预测难于仿真实现。

由此可见，现在迫切需要研发一种适合毛坯件、制造过程复杂、制造工艺未知的各种构件的残余应力场的通用预测方法，为构件的结构设计和制造、使用提供重要支持，以使构件更好地满足工程实际需要。

发明内容

针对现有技术的以上缺陷或改进需求，本发明提供一种基于有限测试点的构件残余应力场的预测方法，以给予构件残余应力场预测，为构件的结构设计和使用提供重要支持。

本发明提供一种基于有限测试点的构件残余应力场的预测方法，包括：

步骤1对构件建模并离散化，将所述构件划分为多个八节点六面体等参单元，其中所述构件整体节点个数为偶数，提取每一节点的坐标值，记录各单元的节点组成信息；

步骤2对各单元的节点进行分类，选取每个单元的四个节点为该单元的应力形质点，未被选取的节点为该单元的非应力形质点，求解单元位移形函数矩阵N、单元应力函数矩阵M，然后根据求得的所述单元位移形函数矩阵N、所述单元应力函数矩阵M，求解单元矩阵k^e；

步骤3建立应力形质点残余应力向量列阵σ、整体节点载荷列阵R，并根据所述应力形质点残余应力向量列阵σ、所述整体节点载荷列阵R和所述单元矩阵k^e，求解整体矩阵K；

步骤4确定实验测试点个数、分布，对所述构件表面残余应力进行无损检测并记录每个实验测试点的6个应力分量；

步骤5根据所述步骤3所得的应力形质点残余应力向量列阵σ、等效载荷列向量R、整体矩阵K，将所述步骤4所得的残余应力测试结果作为残余应力边界条件，代入关系式Kσ＝R求得所有应力形质点残余应力；

步骤6根据所述步骤5所得的应力形质点残余应力和所述步骤2所得的单元应力函数矩阵M，求解非应力形质点残余应力，从而可得所述构件的残余应力场。

总体而言，通过本发明所构思的以上技术方案与现有技术相比，具有以下有益效果：

本发明基于构件中的残余应力场、位移连续以及在研究的时刻残余应力场是平衡的假设，利用有限的表面点的残余应力测试值对构件残余应力场进行预测。本发明对构件进行结构离散化为具有偶数个节点的多个8节点六面体等参单元，并对每个单元节点进行分类划分为4个应力形质点和4个非应力形质点，建立了基于八节点的位移形函数矩阵和基于应力形质点的应力函数矩阵，求解得到单元矩阵，确保了应力形质点应力连续性，提高了预测精度；根据节点受力平衡得出整体矩阵，根据整体矩阵的秩得出测试点数目和分布，大大降低了测试点个数并避免了实验测试的盲目性；通过有限的表面测试点，根据得出的整体方程，先预测构件应力形质点残余应力，再得到构件非应力形质点残余应力，从而得到预测构件的残余应力场，避免了对构件进行破坏实验便实现了得知构件内部残余应力的分布，而且不受构件材料结构形状的限制，尤其在未知构件制造工艺以及复杂构件残余应力预测具有明显优势。

附图说明

图1为本发明基于有限测试点的构件残余应力场的预测方法的流程图；

图2为本发明实施例的结构离散化的示意图；

图3为本发明实施例的单元节点坐标系中位置及排序示意图；

图4为本发明实施例的测点划分示意图；

图5为本发明实施例的测点示意图。

具体实施方式

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

本发明基于构件中的残余应力场、位移连续以及在研究的时刻残余应力场是平衡的假设，鉴于构件表面残余应力是可测的，采用有限元理论，基于实验测试值，进行求解构件残余应力场。

图1所示为本发明基于有限测试点的构件残余应力场的预测方法的流程图，具体包括以下步骤：

步骤1结构离散化：

对构件建模并离散，提取整体节点的坐标值，记录各单元的节点组成信息，具体包括以下子步骤：

(1-1)对构件进行建模，编程或借助有限元软件进行网格划分，将构件划分为多个八节点六面体等参单元，整体节点个数为偶数；

(1-2)提取整体节点坐标值、各单元的节点组成信息，将提取的整体坐标系xyz下的坐标值、各单元的节点信息分别以文本格式保存。在本发明实施例中，保存节点坐标值的节点文件每行由节点号、x坐标、y坐标、z坐标组成，四者之间用逗号隔开，保存单元的节点信息的单元文件每行由单元号、第1节点号～第8节点号组成，之间用空格隔开。

步骤2对节点进行分类、求解单元矩阵k^e：

根据步骤1得到的节点文件和单元文件，根据可观可测、应力连续、位移连续假设、有限元理论，对各单元的节点进行分类，选取每个单元的四个节点为该单元的应力形质点，未被选取的节点为该单元的非应力形质点，求解单元位移形函数矩阵N、单元应力函数矩阵M，然后利用求得的单元位移形函数矩阵N、单元应力函数矩阵M，求解单元矩阵k^e，具体包括以下子步骤：

(2-1)对各单元的节点进行分类：

记构件整体节点个数为n，相邻单元每个单元选取4个节点、整体取0.5n个节点，对单元取点，记单元被选取的节点为单元应力形质点，未被选取的节点为单元非应力形质点；

(2-2)求解单元位移形函数矩阵N：

基于每个单元内部的位移具有连续性假设，记节点i的三个位移分量为a_i＝[u_i v_i w_i]^T(i＝1,2，…，n)，单元位移为a＝[u v w]^T、单元内第j个节点的位移分量为$>a_j^e={\left(\begin{array}{ccc}{u}_j^e& v_j^e& w_j^e\end{array}\right)}^T(j=\mathrm{1,2},...,8),$>单元所有节点位移列向量为$>a^e={\left(\begin{array}{cccccccc}{a}_1^e& a_2^e& a_3^e& a_4^e& a_5^e& a_6^e& a_7^e& a_8^e\end{array}\right)}^T,$>将构件离散为八节点六面体等参单元，则在局部坐标系下等参单元顶点坐标分别为(1，-1，1)、(1，1，1)，(1，1，-1)、(1，-1，-1)、(-1，-1，1)、(-1，1，1)、(-1，1，-1)、(-1，-1，-1)。

在局部坐标系(ξ,η,ζ)下，假设每个单元位移具有连续性，即令：

$>\left(\begin{array}{c}u=b_1+b_2\xi +b_3\eta +b_4\zeta +b_5\mathrm{\xi \eta }+b_6\mathrm{\eta \zeta }+b_7\mathrm{\xi \zeta }+b_8\mathrm{\xi \eta \zeta }\\ v=b_9+b_{10}\xi +b_{11}\eta +b_{12}\zeta +b_{13}\mathrm{\xi \eta }+b_{14}\mathrm{\eta \zeta }+b_{15}\mathrm{\xi \zeta }+b_{16}\mathrm{\xi \eta \zeta }\\ w=b_{17}+b_{18}\xi +b_{19}\eta +b_{20}\zeta +b_{21}\mathrm{\xi \eta }+b_{22}\mathrm{\eta \zeta }+b_{23}\mathrm{\xi \zeta }+b_{24}\mathrm{\xi \eta \zeta }\end{array}\right)---\left(1\right)$>

则单元上的应力场用节点位移表示有：

$>a=\left(\begin{array}{c}u\\ v\\ w\end{array}\right)=\left(\begin{array}{cccccccccc}{N}_1& 0& 0& N_2& 0& 0& ...& N_8& 0& 0\\ 0& N_1& 0& 0& N_2& 0& ...& 0& N_8& 0\\ 0& 0& N_1& 0& 0& N_2& ...& 0& 0& N_8\end{array}\right)\left(\begin{array}{c}{a}_1^e\\ a_2^e\\ a_3^e\\ a_4^e\\ a_5^e\\ a_6^e\\ a_7^e\\ a_8^e\end{array}\right)=N{a}^e---\left(2a\right)$>

$>N=\left(\begin{array}{cccccccccc}{N}_1& 0& 0& N_2& 0& 0& ...& N_8& 0& 0\\ 0& N_1& 0& 0& N_2& 0& ...& 0& N_8& 0\\ 0& 0& N_1& 0& 0& N_2& ...& 0& 0& N_8\end{array}\right)---\left(2b\right)$>

$>N_j=\frac{1}{8}(1+{\xi }_j\xi )(1+{\eta }_j\eta )(1+{\zeta }_j\zeta ),(j=\mathrm{1,2,3,4,5,6,7,8})---\left(2c\right)$>

其中，N为单元位移形函数矩阵，N₁～N₈为单元1～8节点形函数；ξ、η、ζ分别表示在标准单元上的局部坐标，ξ_j、η_j、ζ_j代表单元的第j个节点在局部坐标系ξηζ中的坐标；

同时，记在整体坐标系中x₁～x₈为单元的第1～8节点的x坐标；y₁～y₈为单元的第1～8节点的y坐标；z₁～z₈为单元的第1～8节点的z坐标；单元第j个节点坐标X_j＝[x_j y_j z_j]^T,整体坐标与局部坐标映射关系为：

$>\left(\begin{array}{c}x\\ y\\ z\end{array}\right)=\left(\begin{array}{cccccccccc}{N}_1& 0& 0& N_2& 0& 0& ...& N_8& 0& 0\\ 0& N_1& 0& 0& N_2& 0& ...& 0& N_8& 0\\ 0& 0& N_1& 0& 0& N_2& ...& 0& 0& N_8\end{array}\right)\left(\begin{array}{c}{X}_1\\ X_2\\ X_3\\ X_4\\ X_5\\ X_6\\ X_7\\ X_8\end{array}\right)---\left(3\right)$>

(2-3)求解单元应力函数矩阵M：

假定构件残余应力场的应力是连续性，记节点i的六个应力分量为σ_i＝[σ_ix σ_iy σ_iz τ_xy τ_yz τ_zx]^T,单元应力分量为[σ_x σ_y σ_z τ_xy τ_yz τ_zx]^T，每个单元的四个应力形质点分别记为单元的第h,k,l,m个节点，单元应力形质点残余应力列向量$>{\sigma }^e={\left(\begin{array}{cccc}{\sigma }_h^e& {\sigma }_k^e& {\sigma }_l^e& {\sigma }_m^e\end{array}\right)}^T,$>将单元应力假设为局部坐标系下线性函数，有：

$>\left(\begin{array}{c}{\sigma }_x=c_{11}+c_{12}\xi +c_{13}\eta +c_{14}\zeta \\ {\sigma }_y=c_{21}+c_{22}\xi +c_{13}\eta +c_{24}\zeta \\ {\sigma }_z=c_{31}+c_{32}\xi +c_{33}\eta +c_{34}\zeta \\ {\tau }_{\mathrm{xy}}=c_{41}+c_{42}\xi +c_{43}\eta +c_{44}\zeta \\ {\tau }_{\mathrm{yz}}=c_{51}+c_{52}\xi +c_{53}\eta +c_{54}\zeta \\ {\tau }_{\mathrm{zx}}=c_{61}+c_{62}\xi +c_{63}\eta +c_{64}\zeta \end{array}\right)\left(4\right)$>

其中，σ_x、σ_y、σ_z分别代表单元x、y、z方向应力；τ_xy、τ_yz、τ_zx分别代表单元节点的xy,yz,zx方向切应力；c₁₁～c₆₄表示假定线性系数；

单元上的应力场用4个应力形质点的应力表示，则有：

$>\left(\begin{array}{c}\sigma =\left(\begin{array}{cccccccccccccccccccccccc}{M}_h& 0& 0& 0& 0& 0& M_k& 0& 0& 0& 0& 0& M_l& 0& 0& 0& 0& 0& M_m& 0& 0& 0& 0& 0\\ 0& M_h& 0& 0& 0& 0& 0& M_k& 0& 0& 0& 0& 0& M_l& 0& 0& 0& 0& 0& M_m& 0& 0& 0& 0\\ 0& 0& M_h& 0& 0& 0& 0& 0& M_k& 0& 0& 0& 0& 0& M_l& 0& 0& 0& 0& 0& M_m& 0& 0& 0\\ 0& 0& 0& M_h& 0& 0& 0& 0& 0& M_k& 0& 0& 0& 0& 0& M_l& 0& 0& 0& 0& 0& M_m& 0& 0\\ 0& 0& 0& 0& M_h& 0& 0& 0& 0& 0& M_k& 0& 0& 0& 0& 0& M_l& 0& 0& 0& 0& 0& M_m& 0\\ 0& 0& 0& 0& 0& M_h& 0& 0& 0& 0& 0& M_k& 0& 0& 0& 0& 0& M_l& 0& 0& 0& 0& 0& M_m\end{array}\right)\\ =M{\sigma }^e\end{array}\right)\left(\begin{array}{c}{\sigma }_h^e\\ {\sigma }_k^e\\ {\sigma }_l^e\\ {\sigma }_m^e\end{array}\right)---\left(5a\right)$>

$>\underset{(6\times 24)}{M}=\left(\begin{array}{cccccccccccccccccccccccc}{M}_h& 0& 0& 0& 0& 0& M_k& 0& 0& 0& 0& 0& M_l& 0& 0& 0& 0& 0& M_m& 0& 0& 0& 0& 0\\ 0& M_h& 0& 0& 0& 0& 0& M_k& 0& 0& 0& 0& 0& M_l& 0& 0& 0& 0& 0& M_m& 0& 0& 0& 0\\ 0& 0& M_h& 0& 0& 0& 0& 0& M_k& 0& 0& 0& 0& 0& M_l& 0& 0& 0& 0& 0& M_m& 0& 0& 0\\ 0& 0& 0& M_h& 0& 0& 0& 0& 0& M_k& 0& 0& 0& 0& 0& M_l& 0& 0& 0& 0& 0& M_m& 0& 0\\ 0& 0& 0& 0& M_h& 0& 0& 0& 0& 0& M_k& 0& 0& 0& 0& 0& M_l& 0& 0& 0& 0& 0& M_m& 0\\ 0& 0& 0& 0& 0& M_h& 0& 0& 0& 0& 0& M_k& 0& 0& 0& 0& 0& M_l& 0& 0& 0& 0& 0& M_m\end{array}\right)---\left(5b\right)$>

$>M_j=\frac{1}{4}(1+{\xi }_j+{\eta }_j+{\zeta }_j)---\left(5c\right)$>

其中，h,k,l,m表示单元应力形质点单元中序号；M表示单元应力函数矩阵；M_j代表单元的第j个节点的应力函数；

(2-4)求解单元矩阵k^e：

根据步骤1得到的单元文件、节点文件，节点i所受节点力F_i＝[F_ix F_iy F_iz]^T,单元所受节点力列向量为$>F^e={\left(\begin{array}{ccccccc}{F}_{1x}^e& F_{1y}^e& F_{1z}^e& ...& F_{8x}^e& F_{8y}^e& F_{8z}^e]\end{array}\right)}^T,$>根据虚功原理推理得从而有：

$>k^e={\int }_{V^e}\left|J\right|B^T\mathrm{Md\xi d\eta d\zeta }---\left(6\right)$>

其中，k^e表示24×24矩阵；J表示雅可比矩阵；B表示应变矩阵；M表示单元应力函数矩阵；V^e表示单元体积域；用矩阵形式分别表示如下：

$>\underset{(24\times 24)}{k^e}\left(\begin{array}{ccc}{k}_{\mathrm{1,1}}& ...& k_{\mathrm{1,24}}\\ .& & .\\ .& ...& .\\ .& & .\\ k_{\mathrm{24,1}}& ...& k_{\mathrm{24,24}}\end{array}\right)---\left(7a\right)$>

$>\underset{(3\times 3)}{J}=\left(\begin{array}{cccccccc}\frac{\partial N_1}{\partial \xi }& \frac{\partial N_2}{\partial \xi }& \frac{\partial N_3}{\partial \xi }& \frac{\partial N_4}{\partial \xi }& \frac{\partial N_5}{\partial \xi }& \frac{\partial N_6}{\partial \xi }& \frac{\partial N_7}{\partial \xi }& \frac{\partial N_8}{\partial \xi }\\ \frac{\partial N_1}{\partial \eta }& \frac{\partial N_2}{\partial \eta }& \frac{\partial N_3}{\partial \eta }& \frac{\partial N_4}{\partial \eta }& \frac{\partial N_5}{\partial \eta }& \frac{\partial N_6}{\partial \eta }& \frac{\partial N_7}{\partial \eta }& \frac{\partial N_8}{\partial \eta }\\ \frac{\partial N_1}{\partial \zeta }& \frac{\partial N_2}{\partial \zeta }& \frac{\partial N_3}{\partial \zeta }& \frac{\partial N_4}{\partial \zeta }& \frac{\partial N_5}{\partial \zeta }& \frac{\partial N_6}{\partial \zeta }& \frac{\partial N_7}{\partial \zeta }& \frac{\partial N_8}{\partial \zeta }\end{array}\right)\left(\begin{array}{ccc}{x}_1& y_1& z_1\\ x_2& y_2& z_2\\ x_3& y_3& z_3\\ x_4& y_4& z_4\\ x_5& y_5& z_5\\ x_6& y_6& z_6\\ x_7& y_7& z_7\\ x_8& y_8& z_8\end{array}\right)---\left(7b\right)$>

$>\underset{(6\times 24)}{B}=\left(\begin{array}{cccccccccc}\frac{\partial N_1}{\partial x}& 0& 0& \frac{\partial N_2}{\partial x}& 0& 0& ......& \frac{\partial N_8}{\partial x}& 0& 0\\ 0& \frac{\partial N_1}{\partial y}& 0& 0& \frac{\partial N_2}{\partial y}& 0& ......& 0& \frac{\partial N_8}{\partial y}& 0\\ 0& 0& \frac{\partial N_1}{\partial z}& 0& 0& \frac{\partial N_2}{\partial z}& ......& 0& 0& \frac{\partial N_8}{\partial z}\\ \frac{\partial N_1}{\partial y}& \frac{\partial N_1}{\partial x}& 0& \frac{\partial N_2}{\partial y}& \frac{\partial N_2}{\partial x}& 0& ......& \frac{\partial N_8}{\partial y}& \frac{\partial N_8}{\partial x}& 0\\ 0& \frac{\partial N_1}{\partial z}& \frac{\partial N_1}{\partial y}& 0& \frac{\partial N_2}{\partial z}& \frac{\partial N_2}{\partial y}& ......& 0& \frac{\partial N_8}{\partial z}& \frac{\partial N_8}{\partial y}\\ \frac{\partial N_1}{\partial z}& 0& \frac{\partial N_1}{\partial x}& \frac{\partial N_2}{\partial z}& 0& \frac{\partial N_2}{\partial x}& ......& \frac{\partial N_8}{\partial z}& 0& \frac{\partial N_8}{\partial x}\end{array}\right)---\left(7c\right)$>

$>\left(\begin{array}{c}\frac{\partial N_i}{\partial x}\\ \frac{\partial N_i}{\partial y}\\ \frac{\partial N_i}{\partial z}\end{array}\right)=J^{-1}\left(\begin{array}{c}\frac{\partial N_i}{\partial \xi }\\ \frac{\partial N_i}{\partial \eta }\\ \frac{\partial N_i}{\partial \zeta }\end{array}\right),(i=\mathrm{1,2}...8)---\left(7d\right)$>

步骤3建立应力形质点残余应力向量列阵σ、整体节点载荷列阵R，求解整体矩阵K：

根据步骤2所得的应力形质点节点号，建立应力形质点残余应力向量列阵σ、整体节点载荷列阵R，并利用建立的应力形质点残余应力向量列阵σ、整体节点载荷列阵R和根据步骤2所得的单元矩阵k^e，求解整体矩阵K，具体包括以下子步骤：

(3-1)建立应力形质点残余应力向量列阵σ：

按照节点号从小到大或从大到小的顺序将应力形质点残余应力向量排列，形成行数为3n的应力向量列阵σ；

(3-2)建立整体节点载荷列阵R：

记所有分布体力和面力移植到节点i所得的等效载荷记为节点等效载荷R_i＝[R_ix R_iy R_iz]^T，根据单元等效载荷与原载荷在虚位移上的虚功相等,有体力f、集中力F'、面力作用下的单元节点载荷列R^e为：由于构件残余应力是自然放置卸载状态，通常仅仅受重力，受支撑面的面力，当构件残余应力求解离散为较多单元时，平均到每个节点载荷较小，按照节点号从小到大或从大到小的顺序将所有节点等效载荷排列，形成行数为3n的所有元素为零的整体节点载荷列阵R：

$>\underset{(3n\times 1)}{R}={\left(\begin{array}{cccccccccc}0& 0& 0& 0& 0& 0& ...& 0& 0& 0\end{array}\right)}^T---\left(8\right)$>

(3-3)求解整体矩阵K：

因为在每个节点受力平衡，因为k^eσ^e＝F^e有其中，为对环绕节点i的所有单元求和，根据首先预置全部元素为零的整体矩阵K，单元矩阵中每个单元的k^e(p,q)元素叠加到K(R|_3n×1中R_i(p)所在行，σ|_3n×1中R|_3n×1σ_i(q)所在行数)位置，求得整体矩阵K：

$>\underset{(3n\times 3n)}{K}\left(\begin{array}{ccc}{K}_{11}& ...& K_{1,3n}\\ .& & .\\ .& ...& .\\ .& & .\\ K_{3n,1}& ...& K_{3n,3n}\end{array}\right)---\left(9\right)$>

步骤4实验测试点个数、分布确定及测试：

根据步骤1所得节点文件和步骤3所得的整体矩阵K，求解整体矩阵K的秩r(K)，确定实验测试点个数必大于等于(3n-r(K))/6。在本发明实施例中，实际测试选取时以所选测试点残余应力必须容易测试，不应为构件边上的点，且以不能共线为原则，对构件表面残余应力进行无损检测，每个测点记录6个应力分量并保存。

步骤5应力形质点残余应力求解：

根据步骤3所得的应力形质点残余应力向量列阵σ、等效载荷列向量R、整体矩阵K，将步骤4所得的残余应力测试结果作为残余应力边界条件，代入关系式Kσ＝R求得所有应力形质点残余应力。

步骤6非应力形质点残余应力求解：

根据步骤5所得的应力形质点残余应力和步骤2所得的单元应力函数矩阵M，求解各单元的非应力形质点残余应力,当节点属于多个单元时，对各个单元分别求解后取平均值，可得到非应力形质点残余应力，从而可得构件的残余应力场。

下面结合实施例规格为30mm×10mm×5mm的试样块残余应力预测方法对上述步骤进行具体说明：

步骤1结构离散化，具体包括以下子步骤：

(1-1)对结构进行编程或借助有限元软件进行网格划分，划分时确保节点个数为偶数。为使节点数为偶数，x、y、z方向单元个数应至少有一个方向为奇数，本发明实施例用abaqus软件作为前处理，将单元整体坐标系xyz建立在构件底部中心处，划分成八节点六面体单元，为了便于描述这里仅划分了9个单元作叙述，而非实际划分网格数。划分结果示意图如图2所示，图上节点旁边数字为节点号，单元面上数字为单元号，在在本发明实施例中，此构件共划分为64个节点、27个单元；

(1-2)利用abaqus软件输出inp文件，提取各节点坐标值，各单元节点，分析各单元组成。构件整体坐标系xyz的原点建立在构件底部中心，x方向为长度方向，y方向为宽度方向，z方向为高度方向，每个单元八节点在整体坐标系xyz和局部坐标系ξηζ中的位置及排序如图3所示，其中节点1(1，-1，1)中括号外的数字代表节点序号，括号内为局部坐标系ξηζ下坐标，每个单元第1节点～第8节点排序特点为：右上前点、右上后点、右下后点、右下前点、左上前点、左上后点、左下后点、左下前点。为编程调用方便，将提取的节点信息(即节点在整体坐标系xyz的坐标值)、各单元的节点信息分别以文本格式保存。在本发明实施例中，节点文件每行由节点号、x坐标、y坐标、z坐标组成，四者之间用逗号隔开；单元文件每行由单元号、单元的第1节点号～第8节点号组成，之间用空格隔开，分别保存为E: ode.txt，E:element.txt，单元文件和节点文件格式示例如下表1、2所示。

表1 单元文件(element.txt)每行格式及示例

表2 节点文件(node.txt)每行格式及示例

步骤2单元节点进行分类、求解单元矩阵k^e，具体包括以下子步骤：

(2-1)对单元节点进行分类：

根据步骤1得到的单元文件、节点文件，根据可观可测、应力连续、位移连续假设、有限元理论，为满足相邻单元每个单元选取4个节点作为应力形质点，整体取0.5n＝0.5×64＝32节点为应力形质点的需求，每个单元选取4个节点且每个单元面上选取2个点，且相邻单元节点选择不同，即选取单元第1、3、6、8节点或2、4、5、7节点作为构建单元应力列向量用点，作为单元应力形质点，未选用的点作为非应力形质点。即令1号单元取第2、4、5、7点即整体节点编号的第18、21、1、6号节点作为应力形质点，则与1号单元相邻的2、4、10单元取第1、3、6、8点即整体节点编号的第17、22、2、5号节点作为应力形质点，最后整体应力形质点选取示意图如图4所示，其中标示出的节点代表应力形质点，整个构件节点号为1、3、6、8、9、11、14、16、18、20、21、23、26、28、29、31、33、35、38、40、41、43、46、48、50、52、53、55、58、60、61、63的32个节点为应力形质点，其余为非应力形质点；

(2-2)求解单元位移形函数矩阵N：

假设每个单元位移具有连续性，根据公式(2c),得：

$>\left(\begin{array}{c}{N}_1=\frac{1}{8}(1+\xi )(1-\eta )(1+\zeta )\\ N_2=\frac{1}{8}(1+\xi )(1+\eta )(1+\zeta )\\ N_3=\frac{1}{8}(1+\xi )(1+\eta )(1-\zeta )\\ N_4=\frac{1}{8}(1+\xi )(1-\eta )(1-\zeta )\\ N_5=\frac{1}{8}(1-\xi )(1-\eta )(1+\zeta )\\ N_6=\frac{1}{8}(1-\xi )(1+\eta )(1+\zeta )\\ N_7=\frac{1}{8}(1-\xi )(1+\eta )(1-\zeta )\\ N_8=\frac{1}{8}(1-\xi )(1-\eta )(1-\zeta )\end{array}\right)---\left(10\right)$>

(2-3)求解单元应力函数矩阵M：

局部坐标系ξηζ中，单元第1到8节点坐标分别为(1，-1，1)、(1，1，1)，(1，1，-1)、(1，-1，-1)、(-1，-1，1)、(-1，1，1)、(-1，1，-1)、(-1，-1，-1)，根据公式(5b)和(5c)得单元应力函数矩阵M为：

(2-3)求解每个单元矩阵k^e：

根据步骤1得到的单元文件、节点文件，记每个单元的应力形质点残余应力矩阵为σ^e＝[σ_h σ_k σ_l σ_m]^T，(h,k,l,m＝1、3、6、8或2、4、5、7)、单元所受节点力列向量分别$>F^e={\left(\begin{array}{ccccccc}{F}_{1x}^e& F_{1y}^e& F_{1z}^e& ...& F_{8x}^e& F_{8y}^e& F_{8z}^e]\end{array}\right)}^T,$>三者满足k^eσ^e＝F^e，根据公式(6)求解每个单元矩阵k^e。在本发明实施例中，雅可比矩阵求解式(7b)、(7c)、(7d)中单元位移形函数N_j(j＝1～8)对局部坐标系导数按照表3所示数值计算。

表3 节点坐标、形函数及形函数导数值

步骤3建立应力形质点残余应力向量列阵σ、整体节点载荷列阵R，求解整体刚度矩阵K，具体包括以下子步骤：

(3-1)建立应力形质点残余应力向量列阵σ：

按照应力形质点节点号从小到大的顺序建立应力向量列阵σ,对本发明实施例，有：

$>\left(\begin{array}{c}\underset{(192\times 1)}{\sigma }=\left(\begin{array}{cccccccccccccccc}{\sigma }_1& {\sigma }_3& {\sigma }_6& {\sigma }_8& {\sigma }_9& {\sigma }_{11}& {\sigma }_{14}& {\sigma }_{16}& {\sigma }_{18}& {\sigma }_{20}& {\sigma }_{21}& {\sigma }_{23}& {\sigma }_{26}& {\sigma }_{28}& {\sigma }_{29}& {\sigma }_{31}\end{array}\right)\\ {\left(\begin{array}{cccccccccccccccc}{\sigma }_{33}& {\sigma }_{35}& {\sigma }_{38}& {\sigma }_{40}& {\sigma }_{41}& {\sigma }_{43}& {\sigma }_{46}& {\sigma }_{48}& {\sigma }_{50}& {\sigma }_{52}& {\sigma }_{53}& {\sigma }_{55}& {\sigma }_{58}& {\sigma }_{60}& {\sigma }_{61}& {\sigma }_{63}\end{array}\right)}^T\end{array}\right)---\left(12\right)$>

(3-2)建立整体节点载荷列阵R：

静止自然放置的构件，通常除包含底面节点的单元受外载大小等于重力的面力和体力外，其他的单元都只受重力，对本发明实施例，实际网格划分比较密时，每个节点体力载荷非常小，计算式忽略不计，从而按照节点号进行从小到大的顺序排序，根据公式(8)得整体节点载荷列阵R：

$>\left(\begin{array}{c}\underset{(192\times 1)}{R}={\left(\begin{array}{ccccccc}{R}_{1x}& R_{1y}& R_{1z}& ......& R_{64x}& R_{64y}& R_{64z}\end{array}\right)}^T\\ ={\left(\begin{array}{ccccccc}0& 0& 0& ......& 0& 0& 0\end{array}\right)}^T\end{array}\right)---\left(13\right)$>

(3-3)求解整体矩阵K：

预置一个整体矩阵并全部元素预置为零，每个单元矩阵中每个单元的k^e(p,q)元素叠加到K(中(p)所在行，中(q)所在行数)位置，求得整体矩阵K：

$>\underset{(192\times 192)}{K}\left(\begin{array}{ccc}{K}_{11}& ...& K_{1,192}\\ .& & .\\ .& ...& .\\ .& & .\\ K_{192,1}& ...& K_{192,192}\end{array}\right)---\left(14\right)$>

步骤4确定实验测试点个数、分布并测试，具体包括：

根据步骤1所得节点文件和所步骤3所得的整体矩阵K，求解整体矩阵K的秩r(K)，在本发明实施例中，n＝64,r(K)＝183，则选取实验测试点个数为2×(3×64-193)/3＝6个，可选取18、35、6、11、21、41六个节点进行测试，测试点示意图如图5所示，其中标示出的节点代表测试点位置，旁边数字代表编号；根据步骤1所得的节点坐标文件查找节点号对应坐标，对实验测试点进行精确定位，对构件进行无损检测，将得到的每个点六个应力分量保存，保存格式为txt文件，每行为：节点号(空格)1～6数字(1代表σ_x，2代表σ_y，3代表σ_z，4代表τ_xy，5代表τ_yz，6代表τ_zx)(空格)残余应力数值。示例如表4所示。在本发明实施例中，利用X射线仪器进行实验数据测试，保存为constrain.txt文件，记录格式及示例如表4所示。

表4 测试数据文件(constrain.txt)每行格式及示例

步骤5求解应力形质点残余应力，具体包括：

根据步骤3所得的整体矩阵K、等效载荷列向量R和步骤4所测残余应力实验结果，代入关系式Kσ＝R，可求得应力向量中的未知量。在本发明实施例中从而32个应力形质点1、3、6、8、9、11、14、16、18、20、21、23、26、28、29、31、33、35、38、40、41、43、46、48、50、52、53、55、58、60、61、63的应力都可求得。

步骤6求解非应力形质点残余应力，具体包括：

根据步骤5所得的应力形质点残余应力σ和形函数M，本发明实施例单元号为奇数的所有单元都选取了2、4、5、7点，则单元1、3、6、8点残余应力按公式(15a)求解，选取1、3、6、8点的所有偶数单元上，2、4、5、7点残余应力按公式(15b)求解：

$>\left(\begin{array}{c}{\sigma }_1=\frac{1}{2}({\sigma }_2+{\sigma }_4+{\sigma }_5-{\sigma }_7),{\sigma }_3=\frac{1}{2}({\sigma }_2+{\sigma }_4-{\sigma }_5+{\sigma }_7),\\ {\sigma }_6=\frac{1}{2}({\sigma }_2-{\sigma }_4+{\sigma }_5+{\sigma }_7),{\sigma }_8=\frac{1}{2}(-{\sigma }_2+{\sigma }_4+{\sigma }_5+{\sigma }_7)\end{array}\right)---\left(15a\right)$>

$>\left(\begin{array}{c}{\sigma }_2=\frac{1}{2}({\sigma }_1+{\sigma }_3+{\sigma }_6-{\sigma }_8),{\sigma }_4=\frac{1}{2}({\sigma }_1+{\sigma }_3-{\sigma }_6+{\sigma }_8),\\ {\sigma }_5=\frac{1}{2}({\sigma }_1-{\sigma }_3+{\sigma }_6+{\sigma }_8),{\sigma }_7=\frac{1}{2}(-{\sigma }_1+{\sigma }_3+{\sigma }_6+{\sigma }_8)\end{array}\right)---\left(15b\right)$>

例如，节点3属于单元1、2，则在单元1、2中分别求解，然后取平均值，依次类推，从而可得到全部非应力形质点：2、4、5、7、10、12、13、15、17、19、22、24、25、27、30、32、34、33、37、39、42、44、45、47、49、51、54、56、57、59、62、64的残余应力，从而可得到所有节点的残余应力，即得到构件的残余应力场。

由于采用了上述技术方案，本发明的方法利用位移连续、应力连续假设，求解单元4节点应力函数和8节点位移形函数，求解得到单元矩阵和整体矩阵，根据整体矩阵的秩进一步判定了测试点数目和分布，避免了测试的盲目性，不仅能够预测构件的残余应力场，而且不受构件材料结构形状的限制，尤其在未知构件制造工艺以及复杂构件残余应力预测具有明显优势。

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