一种螺栓连接结构有限元模型的精细建模方法

        在利用有限元方法对螺栓连接结构进行强度和应力分布等方面的数值模拟时，由于螺纹几何的复杂性，建立螺纹连接结构的精细有限元模型往往较为困难。尽管商业有限元前处理软件可以生成四面体自由网格，能够一定程度上解决螺纹网格划分困难的问题。但这样生成的网格往往具有一定随机性，并且容易产生畸变的网格单元。为此，本文提供了一种基于六面体单元的螺纹网格建模方法，该方法建模过程简单便捷，所得模型规律性强，无需针对螺纹网格建模进行复杂的前处理工作。

1

螺栓外螺纹建模方法

1.1  单螺距螺栓建模

        由于整个螺栓的外螺纹网格可以看作是由若干个单螺距螺纹网格堆叠而成，因此本文首先讨论单螺距螺栓的建模方法。此外，标准螺纹的牙型较为复杂，不便于直观地采用数学表达式来描述螺纹的几何形状变化，因此首先以最简单的直角牙型（牙型角α=90º）为例，介绍单螺距螺栓网格建模的基本原理，然后基于该原理将其推广至标准规格牙型，甚至是任意几何形状的牙型。

1.1.1 直角牙型

        假定在一个螺距P的高度内，θ从0线性增长至2π，根据图1.1中的几何关系，很容易得到牙型轮廓上任意一点到螺栓中心轴线的径向距离r与θ的关系可表示为：

式中：d为螺栓的公称直径。

        如果将上式看作是极坐标方程，利用

        将其转换为平面直角坐标系下的坐标，可以得到如图1.2所示的轮廓。实际上，该轮廓即为螺栓某一横截面上的外螺纹线。

图1.2 外螺纹轮廓线示意图(直角牙型)

        假定螺栓横截面上的外螺纹线位于XOY平面，而螺栓长度方向定义为Z向，分别沿Z=0、Z=P/2和Z=P的位置处截取螺栓，可以得到三层截面上的外螺纹轮廓线如图1.3所示。

图1.3 多层外螺纹轮廓线示意图

        在螺纹的形成过程中，第一层截面上的任意一点都将沿着螺栓的长度方向螺旋上升，如图1.3中与X轴成π/3夹角的红色标记点，该点在螺旋上升的过程中，当到达Z=P/2的高度时（即第二层截面），此时与X轴的夹角将恰好为(π/3+π)，在继续上升至Z=P的高度时（即第三层截面），此时红色标志点回到原来的位置，即此时与X轴的夹角为π/3(+2π)。对于外螺纹线上的任意一点，这样的规律都是成立的，因此第二层截面的外螺纹轮廓线可以看作是将第一层截面的外螺纹轮廓线逆时针旋转180º得到的，而第三层截面的外螺纹轮廓线可以看作是将第一层截面的外螺纹轮廓线线逆时针旋转360º得到的。

        实际上，对于距原点高度为z的任意一层截面，其外螺纹轮廓线可以看作是在第一层截面的外螺纹轮廓线的基础上逆时针旋转2π×(z/P)弧度得到的。由于第一层外螺纹轮廓线的数学表达式和几何形状是已知的（见图1.2），因此在第一层外螺纹线的基础上很容易得到任意位置处的外螺纹轮廓线。

        对于采用结构化网格划分的螺栓外螺纹有限元模型，该模型可以看作是由若干个六面体单元离散而成，因此只要根据上述规律确定了各个离散节点之间的位置关系，即可完成螺栓外螺纹的建模工作。本文首先讨论外螺纹最表层网格节点之间的位置关系，假定在未形成外螺纹网格之前的初始网格节点分布如图1.4所示，这些离散节点构成了一个直径为螺栓公称直径d、高度为螺距P的圆带（可以看做是螺杆的外表面），其中螺栓的中心轴线位于全局坐标系的原点。

图1.4 初始网格节点构成的圆带

        假定初始网格中任意一个离散节点的坐标为(x,y,z)，则在该节点所在的横截面内，节点与全局坐标系的X轴的夹角θ_ini可表示为：

式中：θ_ini的取值范围为[0,2π]，并且定义为逆时针旋转为正。

        由图1.3可知，在形成外螺纹网格之后，点(x,y,0)到螺栓中心轴线的径向距离r_ini可表示为：

        在螺纹的螺旋上升过程中，由于高度上升引起的角度增量θ_add为：

        因此，在形成外螺纹网格之后，点(x,y,z)到螺栓中心轴线的径向距离r应在r_ini的基础上包含外螺纹轮廓线旋转上升引起的变化量，令

则有

        点(x,y,z)在修改之后的坐标为(x',y',z')，可表示为：

        需要注意的是，由于在修改节点的过程中节点并未发生旋转运动，只是假想外螺纹轮廓线在上升的过程中会发生旋转，因此上式中的角坐标仍为θ_ini而非θ。

        基于上述公式，遍历图1.4中的所有离散节点以修改节点坐标，最终得到形成外螺纹之后的离散节点分布如图1.5所示。

图1.5 形成外螺纹之后的离散节点分布（直角牙型）

        图1.5中的三条红色曲线为公称直径线，对比离散节点可以看出，位于第二条红色曲线所在截面的那一层离散节点是在最底层离散节点的基础上旋转180º得到的。

        在得到形成外螺纹的最表层离散节点后，内部离散节点的位置很容易通过最表层离散节点的坐标线性插值获得，线性插值的计算过程如图1.6所示。

（a）初始网格

（b）形成外螺纹后的网格

图1.6 内部节点线性插值示意图

        图1.6(a)中的离散节点排布是在图1.4中圆带的基础上向内部填充节点获得的，此时初始网格为一个横截面为矩形的圆环。假定圆环内径为r_inter，外径为r_exter，圆环内部任意一个离散节点到螺栓中心轴线的径向距离为r_dist。为了便于形成外螺纹网格之后螺栓芯部网格的填充，假定变形之后的网格的内径仍为r_inter，由于形成外螺纹后最边缘节点到螺栓中心轴线的径向距离r是已知的，因此内部节点在变形之后的径向距离r_true满足如下关系：

        在获得内部节点到螺栓中心轴线的径向距离后，很容易就能得到内部节点在形成外螺纹之后的节点坐标。

        基于上述原理，本文首先在有限元前处理软件HyperMesh中建立高度为单个螺距P的圆环网格，通过输出Abaqus求解文件(inp文件)获得网格节点的定义信息，基于节点坐标对这些节点进行处理，获得形成外螺纹网格之后的节点坐标，最后将修改之后的节点坐标覆盖于原来的求解文件中，最后导入到HyperMesh即可获得直角牙型的螺栓外螺纹网格模型，如图1.7所示。

图1.7 螺栓外螺纹网格（直角牙型）

1.1.2 标准规格牙型

        从前面的分析可以看到，建立螺栓外螺纹网格的关键在于根据沿螺栓中心轴线的螺纹截面的几何关系，确定牙型轮廓任意一点到螺栓中心轴线的径向距离r与θ的数学表达式。对于ISO标准给出的标准规格的螺纹牙型，沿螺栓中心轴线的螺纹截面如图1.8所示。为了防止螺纹齿根处出现过度的应力集中现象，要求齿根半径ρ不得小于0.125P。

图1.8 沿螺栓中心轴线的螺纹截面（标准牙型轮廓）

        根据图1.8中的几何关系，可以推导得到采用标准规格牙型的外螺纹截面中径向坐标r与周向坐标θ的表达式为：

其中

式中：H为螺纹工作高度。

        利用上述公式可以计算得到外螺纹的轮廓线如图1.9所示。

图1.9 外螺纹轮廓线示意图（标准牙型）

        同样基于前面的处理方法，可以得到标准牙型下的外螺纹网格如图1.10所示。

图1.10 螺栓外螺纹网格（标准牙型）

1.2  螺栓退刀槽建模

        前面主要介绍了螺栓外螺纹网格的建模，对于螺栓网格模型而言，一般外螺纹网格还需要与无螺纹的螺杆网格连接在一起，为了满足网格的连续性要求，需要保证接合面两侧节点重合。而从图1.11中可以看出，螺栓外螺纹网格并不能直接与螺杆网格连接，而是需要建立一段过渡网格将外螺纹网格与螺杆网格接合在一起，即螺栓的退刀槽网格。这一点从图1.9中也可以看出，因为外螺纹轮廓线和公称直径仅在右侧的BF段是重合的。

图1.11 外螺纹和螺杆网格不匹配示意图

        为了解决外螺纹和螺杆网格的匹配问题，将外螺纹网格和螺杆网格通过一段退刀槽网格进行连接，假定退刀槽的螺距为P_groove，在该段距离内，退刀槽表层的节点到螺栓中心轴线的径向距离应缓慢从r线性过渡至d/2(即r_exter)，如图1.12所示。

（a）初始网格

（b）形成退刀槽后的网格

图1.12 退刀槽表层节点线性插值示意图

        显然，对于退刀槽最表层任意一个到坐标原点高度为z的离散节点，该节点到螺栓中心轴线的径向距离r_true应满足如下所示的关系：

        为了便于区分外螺纹网格和退刀槽网格，可以在生成初始网格时将外螺纹网格放置于全局坐标系Z轴的正半轴，退刀槽网格放置于Z轴的负半轴。在计算得到退刀槽最表层的节点后，同样利用1.1.1节中的方法计算退刀槽内部节点的坐标，即可得到包含退刀槽的螺栓外螺纹网格，如图1.13所示。由于该网格最底层的离散节点构成了一个公称直径大小的圆，因此可以与无螺纹的螺杆网格直接相连。

图1.13 包含退刀槽的螺栓外螺纹网格

1.3  任意长度螺栓建模

        以上内容主要讨论单螺距螺栓中外螺纹和退刀槽网格的建模方法，尽管多螺距的螺栓网格可以通过若干个单螺距的螺栓网格堆叠得到，但为了便于建模，本文仍希望通过上述方法建立任意长度的初始网格即可获得对应的螺栓外螺纹网格。前面已经提到，螺纹在螺旋上升过程中径向长度r的变化是通过高度上升引起的角度增量θ_add来体现的。由于高度每上升一个螺距P，角度增量θ_add即增加2π，这将使得外螺纹的轮廓线重新回到原来的位置，显然对于距坐标原点为任意高度z的离散节点，由于高度上升引起的角度增量θ_add可表示为：

式中：mod表示取余运算。

        对于任意长度的螺栓建模，除以上计算公式有所区别以外，其余计算流程与单螺距螺栓的建模方法完全相同。

1.4  螺栓芯部网格填充方法

        在采用本文所述方法建立螺栓的外螺纹网格之后，还需要填充螺栓芯部的空白区域。一种填充方法为通过拉伸或旋转的方式生成结构化的六面体网格进行逐步过渡，另一种填充方法为直接在空白区域内采用二维网格进行包覆，形成封闭几何体，然后利用HyperMesh的Tera Mesh功能生成四面体网格。两种填充方案得到的螺栓网格如图1.14所示。

（a）六面体填充

（b）四面体填充

图1.14 螺栓芯部网格填充方案

2

螺母内螺纹建模方法

        螺母内螺纹的建模方法与螺栓外螺纹的建模方法基本相同。对于ISO标准给出的标准规格的螺纹牙型，沿螺母中心轴线的内螺纹截面如图2.1所示。

图2.1 沿螺母中心轴线的内螺纹截面

        根据图2.1中的几何关系，同样可以推导得到采用标准规格牙型的内螺纹截面中径向坐标r与周向坐标θ的表达式为：

式中：

        与螺栓外螺纹的建模方式有所不同的是，由于螺母的螺纹为内螺纹，因此初始网格的内径应和螺栓的公称直径一致，并且需要首先计算内表面离散节点到螺母中心轴线的径向距离r，然后通过插值的方式计算外部节点到螺母中心轴线的径向距离r_true，此时满足的关系式变为：

        基于上述方法建立螺母内螺纹网格模型，与前面建立的螺栓外螺纹网格装配后即可得到螺栓连接结构模型，如图2.2所示。

图2.2 螺栓连接结构网格模型

        这里对本文提出的螺栓连接结构网格的建模方法做几点说明：首先，由于螺栓的螺纹是螺旋上升的，因此整个网格不存在几何对称性；其次，本文讨论的计算方法只涉及对离散节点坐标的计算，并未严格限制初始网格的划分方式，仅要求初始网格为圆柱体，例如，图2.3给出了初始网格完全为四面体网格时获得的模型，可以看到该建模方法具有较高的通用性；最后，虽然本文针对标准牙型的外螺纹和内螺纹均给出了相应的计算表达式，但显然这并非是应用该方法的必要条件，只要建立了径向坐标r和角坐标θ之间的关系（可以事先建立一组描述该关系的离散数据，然后通过插值实时获取相应的取值），就可以实现任意牙型的螺纹网格建模。

图2.3 螺栓连接结构网格模型（四面体网格）

螺栓连接结构建模演示

附录

%************************************************************************
%程序名:Bolt_modeling.m
%程序用于建立标准螺栓外螺纹的精细模型(支持任意长度的螺纹建模)
%螺栓含有退刀槽(退刀槽位于Z轴负半轴,长度为一个退刀槽螺距)
%************************************************************************
clear,clc
%********************************建模参数输入******************************
d=12;               %螺栓的公称直径
r_exter=6;          %有限元模型外径(和螺栓的公称直径对应)
r_inter=4;          %有限元模型内径
P=1.75;             %螺距
P_groove=1;         %退刀槽螺距
Rho=0.13*P;         %牙型半径
H=sqrt(3)/2*P;
Theta1=pi/8;
Theta2=(1-sqrt(3)/P*Rho)*pi;
Theta3=(1+sqrt(3)/P*Rho)*pi;
Theta4=15/8*pi;
%**************************************************************************
%*************************绘制外螺纹轮廓线*********************************
theta_plot=linspace(0,2*pi,1001)';    %极坐标的角度样点
%绘制公称直径线
x_plot=d/2*cos(theta_plot);
y_plot=d/2*sin(theta_plot);
figure;plot(x_plot,y_plot,'--','LineWidth',1);axis equal;hold on;
%绘制外螺纹轮廓线
r_plot=zeros(1001,1);                  %极坐标的半径样点
fori=1:1001
if (theta_plot(i)>=0)&&(theta_plot(i)<=Theta1)
        r_plot(i)=d/2;
elseif (theta_plot(i)>=Theta1)&&(theta_plot(i)<=Theta2)
        r_plot(i)=d/2+H/8-H/pi*theta_plot(i);
elseif (theta_plot(i)>=Theta2)&&(theta_plot(i)<=pi)
        r_plot(i)=d/2-7/8*H+2*Rho-...
sqrt(Rho^2-(pi-theta_plot(i))^2/(4*pi^2)*P^2);
elseif (theta_plot(i)>=pi)&&(theta_plot(i)<=Theta3)
        r_plot(i)=d/2-7/8*H+2*Rho-...
sqrt(Rho^2-(theta_plot(i)-pi)^2/(4*pi^2)*P^2);
elseif (theta_plot(i)>=Theta3)&&(theta_plot(i)<=Theta4)
        r_plot(i)=d/2-15/8*H+theta_plot(i)/pi*H;
else
        r_plot(i)=d/2;
end
end
x_plot=r_plot.*cos(theta_plot);
y_plot=r_plot.*sin(theta_plot);
plot(x_plot,y_plot,'r','LineWidth',2);hold off;
legend('公称直径','外螺纹轮廓线');title('外螺纹轮廓线示意图');
figure
plot(r_plot,theta_plot/(2*pi)*P,'r','LineWidth',2);axis equal;
title('沿螺栓中心轴线的螺纹截面');
%**************************************************************************
%****************************处理有限元网格节点****************************
%读取初始网格节点
Mesh_Data=importdata('MESH_DATA.txt');
Node_num=size(Mesh_Data,1);             %网格总节点数
%绘制初始网格节点分布
figure;
scatter3(Mesh_Data(:,2),Mesh_Data(:,3),Mesh_Data(:,4));
%逐节点循环修改节点坐标
Mesh_Data_Modified=zeros(Node_num,4);   %用于存放修改之后的节点坐标
fori=1:Node_num
%读取初始网格节点坐标
    Coord_X=Mesh_Data(i,2);             %X坐标
    Coord_Y=Mesh_Data(i,3);             %Y坐标
    Coord_Z=Mesh_Data(i,4);             %Z坐标
%计算节点在对应螺纹截面的角度(假定螺纹截面位于XOY平面)
if Coord_X>=0&&Coord_Y>=0%节点位于截面的第一象限
        theta_ini=atan(Coord_Y/Coord_X);
elseif Coord_X<0&&Coord_Y>=0%节点位于截面的第二象限
        theta_ini=pi+atan(Coord_Y/Coord_X);
elseif Coord_X<0&&Coord_Y<0%节点位于截面的第三象限
        theta_ini=pi+atan(Coord_Y/Coord_X);                  
else%节点位于截面的第四象限
        theta_ini=2*pi+atan(Coord_Y/Coord_X);
end
%计算由于螺纹截面高度变化引起的节点角度增量
    theta_add=2*pi*(Coord_Z/P);
%计算节点真实的角度
    theta=theta_ini+theta_add;
    theta=mod(theta,2*pi);       %角度超过2pi后需要进行取余运算
%计算最外层节点修改后的径向坐标
if Coord_Z>=0%节点位于正常螺纹部分
if (theta>=0)&&(theta<=Theta1)
            r=d/2;
elseif (theta>=Theta1)&&(theta<=Theta2)
            r=d/2+H/8-H/pi*theta;
elseif (theta>=Theta2)&&(theta<=pi)
            r=d/2-7/8*H+2*Rho-...
sqrt(Rho^2-(pi-theta)^2/(4*pi^2)*P^2);
elseif (theta>=pi)&&(theta<=Theta3)
            r=d/2-7/8*H+2*Rho-...
sqrt(Rho^2-(theta-pi)^2/(4*pi^2)*P^2);
elseif (theta>=Theta3)&&(theta<=Theta4)
            r=d/2-15/8*H+theta/pi*H;
else
            r=d/2;
end
else%节点位于退刀槽部分
%计算Z=0位置处的螺纹线
if (theta_ini>=0)&&(theta_ini<=Theta1)
            r_ini=d/2;
elseif (theta_ini>=Theta1)&&(theta_ini<=Theta2)
            r_ini=d/2+H/8-H/pi*theta_ini;
elseif (theta_ini>=Theta2)&&(theta_ini<=pi)
            r_ini=d/2-7/8*H+2*Rho-...
sqrt(Rho^2-(pi-theta_ini)^2/(4*pi^2)*P^2);
elseif (theta_ini>=pi)&&(theta_ini<=Theta3)
            r_ini=d/2-7/8*H+2*Rho-...
sqrt(Rho^2-(theta_ini-pi)^2/(4*pi^2)*P^2);
elseif (theta_ini>=Theta3)&&(theta_ini<=Theta4)
            r_ini=d/2-15/8*H+theta_ini/pi*H;
else
            r_ini=d/2;
end
%根据Z=0处的螺纹线线性插值得到退刀槽对应的螺纹线
ifabs(Coord_Z)<=P_groove    %仅在一个退刀槽长度内更改网格节点
            r=(d/2-r_ini)/P_groove*abs(Coord_Z)+r_ini;
else
            r=d/2;                   %超出退刀槽长度的网格保持不变
end
end
%计算未修改节点到原点的径向距离
    r_dist=sqrt(Coord_X^2+Coord_Y^2); 
%计算内层网格到原点的真实径向距离(线性插值)
    r=(r_dist-r_inter)/(r_exter-r_inter)*(r-r_inter)+r_inter;     
%根据节点的径向坐标和角坐标计算笛卡尔坐标
    Coord_X_Modified=r*cos(theta_ini);
    Coord_Y_Modified=r*sin(theta_ini);
%存放修改后的节点信息到新的矩阵
    Mesh_Data_Modified(i,1)=Mesh_Data(i,1);
    Mesh_Data_Modified(i,2)=Coord_X_Modified;
    Mesh_Data_Modified(i,3)=Coord_Y_Modified;
    Mesh_Data_Modified(i,4)=Coord_Z;    
end
%绘制修改之后的网格节点分布
figure;
scatter3(Mesh_Data_Modified(:,2),Mesh_Data_Modified(:,3),...
    Mesh_Data_Modified(:,4));
axis equal;hold on;
scatter3(d/2*cos(theta_plot),d/2*sin(theta_plot),zeros(1001,1));
%**************************************************************************
%************************修改后的网格节点信息输出**************************
fileID=fopen('MESH_DATA_Modified.txt','w');
fori=1:Node_num
    fprintf(fileID,'%8d,',Mesh_Data_Modified(i,1));          %节点编号
    fprintf(fileID,'%20.15g,',Mesh_Data_Modified(i,2));      %X坐标
    fprintf(fileID,'%20.15g,',Mesh_Data_Modified(i,3));      %Y坐标
    fprintf(fileID,'%20.15g\n',Mesh_Data_Modified(i,4));     %Z坐标
end
fclose(fileID);
%**************************************************************************

%************************************************************************
%程序名:Nut_modeling.m
%程序用于建立标准螺母内螺纹的精细模型(支持任意长度的螺纹建模)
clear,clc
%********************************建模参数输入******************************
d=12;                           %螺栓的公称直径
r_exter=8;                      %有限元模型外径
r_inter=6;                      %有限元模型内径(和螺栓的公称直径对应)
P=1.75;                         %螺距
Rho=0.12;                       %牙型半径
H=sqrt(3)/2*P;
d1=d-2*5/8*H;                   %螺纹小径
Theta1=sqrt(3)*Rho/P*pi;
Theta2=3*pi/4;
Theta3=5*pi/4;
Theta4=(2-sqrt(3)*Rho/P)*pi;
%**************************************************************************
%*************************绘制内螺纹轮廓线*********************************
theta_plot=linspace(0,2*pi,1001)';    %极坐标的角度样点
%绘制公称直径线
x_plot=d/2*cos(theta_plot);
y_plot=d/2*sin(theta_plot);
figure;plot(x_plot,y_plot,'--','LineWidth',1);axis equal;hold on;
%绘制螺纹小径轮廓线
x_plot=d1/2*cos(theta_plot);
y_plot=d1/2*sin(theta_plot);
plot(x_plot,y_plot,'--','LineWidth',1)
%绘制内螺纹轮廓线
r_plot=zeros(1001,1);                  %极坐标的半径样点
fori=1:1001
if (theta_plot(i)>=0)&&(theta_plot(i)<=Theta1)
        r_plot(i)=d/2+H/8-2*Rho+...
sqrt(Rho^2-theta_plot(i)^2/(4*pi^2)*P^2);
elseif (theta_plot(i)>=Theta1)&&(theta_plot(i)<=Theta2)
        r_plot(i)=d/2+H/8-H/pi*theta_plot(i);
elseif (theta_plot(i)>=Theta2)&&(theta_plot(i)<=Theta3)
        r_plot(i)=d1/2;
elseif (theta_plot(i)>=Theta3)&&(theta_plot(i)<=Theta4)
        r_plot(i)=d/2-15*H/8+H/pi*theta_plot(i);
else
        r_plot(i)=d/2+H/8-2*Rho+...
sqrt(Rho^2-(2*pi-theta_plot(i))^2/(4*pi^2)*P^2);
end
end
x_plot=r_plot.*cos(theta_plot);
y_plot=r_plot.*sin(theta_plot);
plot(x_plot,y_plot,'r','LineWidth',2);hold off;
legend('公称直径','螺纹小径','内螺纹轮廓线');title('内螺纹轮廓线示意图');
figure
plot(r_plot,theta_plot/(2*pi)*P,'r','LineWidth',2);axis equal;
title('沿螺母中心轴线的螺纹截面');
%**************************************************************************
%****************************处理有限元网格节点*******************************
%读取初始网格节点
Mesh_Data=importdata('MESH_DATA.txt');
Node_num=size(Mesh_Data,1);             %网格总节点数
%绘制初始网格节点分布
figure;
scatter3(Mesh_Data(:,2),Mesh_Data(:,3),Mesh_Data(:,4));
%逐节点循环修改节点坐标
Mesh_Data_Modified=zeros(Node_num,4);   %用于存放修改之后的节点坐标
fori=1:Node_num
%读取初始网格节点坐标
    Coord_X=Mesh_Data(i,2);             %X坐标
    Coord_Y=Mesh_Data(i,3);             %Y坐标
    Coord_Z=Mesh_Data(i,4);             %Z坐标
%计算节点在对应螺纹截面的角度(假定螺纹截面位于XOY平面)
if Coord_X>=0&&Coord_Y>=0%节点位于截面的第一象限
        theta_ini=atan(Coord_Y/Coord_X);
elseif Coord_X<0&&Coord_Y>=0%节点位于截面的第二象限
        theta_ini=pi+atan(Coord_Y/Coord_X);
elseif Coord_X<0&&Coord_Y<0%节点位于截面的第三象限
        theta_ini=pi+atan(Coord_Y/Coord_X);                  
else%节点位于截面的第四象限
        theta_ini=2*pi+atan(Coord_Y/Coord_X);
end
%计算由于螺纹截面高度变化引起的节点角度增量
    theta_add=2*pi*(Coord_Z/P);
%计算节点真实的角度
    theta=theta_ini+theta_add;
    theta=mod(theta,2*pi);      %角度超过2pi后需要进行取余运算
%计算最内层节点修改后的坐标(不考虑退刀槽)
if (theta>=0)&&(theta<=Theta1)
        r=d/2+H/8-2*Rho+sqrt(Rho^2-theta^2/(4*pi^2)*P^2);
elseif (theta>=Theta1)&&(theta<=Theta2)
        r=d/2+H/8-H/pi*theta;
elseif (theta>=Theta2)&&(theta<=Theta3)
        r=d1/2;
elseif (theta>=Theta3)&&(theta<=Theta4)
        r=d/2-15*H/8+H/pi*theta;
else
        r=d/2+H/8-2*Rho+sqrt(Rho^2-(2*pi-theta)^2/(4*pi^2)*P^2);
end
%计算未修改节点到原点的径向距离
    r_dist=sqrt(Coord_X^2+Coord_Y^2); 
%计算外层网格到原点的真实径向距离(线性插值)
    r=r_exter-(r_exter-r_dist)/(r_exter-r_inter)*(r_exter-r);    
%根据节点的径向坐标和角坐标计算笛卡尔坐标
    Coord_X_Modified=r*cos(theta_ini);
    Coord_Y_Modified=r*sin(theta_ini);
%存放修改后的节点信息到新的矩阵
    Mesh_Data_Modified(i,1)=Mesh_Data(i,1);
    Mesh_Data_Modified(i,2)=Coord_X_Modified;
    Mesh_Data_Modified(i,3)=Coord_Y_Modified;
    Mesh_Data_Modified(i,4)=Coord_Z;
end
%绘制修改之后的网格节点分布
figure;
scatter3(Mesh_Data_Modified(:,2),Mesh_Data_Modified(:,3),...
    Mesh_Data_Modified(:,4));
axis equal;hold on;
scatter3(d/2*cos(theta_plot),d/2*sin(theta_plot),zeros(1001,1));
%**************************************************************************
%************************修改后的网格节点信息输出******************************
fileID=fopen('MESH_DATA_Modified.txt','w');
fori=1:Node_num
    fprintf(fileID,'%8d,',Mesh_Data_Modified(i,1));          %节点编号
    fprintf(fileID,'%20.15g,',Mesh_Data_Modified(i,2));      %X坐标
    fprintf(fileID,'%20.15g,',Mesh_Data_Modified(i,3));      %Y坐标
    fprintf(fileID,'%20.15g\n',Mesh_Data_Modified(i,4));     %Z坐标
end
fclose(fileID);
%**************************************************************************