基于临界平面法的多轴疲劳分析

%程序用于读取多个单元中心在单个加载循环的应力应变历程并计算SWT参量clear,clc%*****************************计算参数设置**********************************Elem_num=2400;                                       %待提取节点数量Step_inc_num=51;                                    %分析步增量数量Theta_num=181;                                       %临界平面角度theta分段数Phi_num=181;                                         %临界平面角度phi分段数%**************************************************************************
%*******************************变量初始化**********************************Theta=linspace(0,pi,Theta_num);Phi=linspace(0,pi,Phi_num);SWT=zeros(Elem_num,1);                              %存放每个待评估点的SWTTheta_Critical=zeros(Elem_num,1);                   %存放每个待评估点的临界角度Phi_Critical=zeros(Elem_num,1);                     %存放每个待评估点的临界角度Elem_eval_centroid_coord_X=zeros(Elem_num,1);       %待评估单元中心X坐标Elem_eval_centroid_coord_Y=zeros(Elem_num,1);       %待评估单元中心Y坐标Elem_eval_centroid_coord_Z=zeros(Elem_num,1);       %待评估单元中心Z坐标%**************************************************************************
%*********************读取单元中心应力应变数据*****************************File_input_id=fopen("abaqus.rpt");                    %读取Abaqus数据文件File_output_id=fopen("Stress_Strain_Hist.txt","w");%创建文件存放应力应变数据while ~feof(File_input_id)                    %判断指针是否已指向文件末尾行  tline=fgetl(File_input_id);                 %读取当前指针指向行的文件内容  if ~isempty(tline)                          %跳过空行    if numel(tline)>=11                       %跳过字符数少于10个的文件行      if double(tline(16))>=48&&double(tline(16))<=57        %根据ASCII码判断文件行的第12个字符是否为数字 满足条件则为数据行        fprintf(File_output_id,'%s\n\n',tline); %写入数据到应力应变数据文件      end    end  endendfclose(File_output_id);                         %关闭文件fclose(File_input_id);SE_Data=importdata("Stress_Strain_Hist.txt");   %导入应力应变数据%**************************************************************************
%***************************整理应力应变数据*******************************%初始化元胞数组存放数据Stress_Data=cell(Elem_num,1);                   %存放所有单元中心的应力数据Strain_Data=cell(Elem_num,1);                   %存放所有单元中心的应变数据%格式化存储数据(元胞数组每个元素为一个单元中心在所有增量步的数据)for i=1:Elem_num  for j=1:Step_inc_num    Strain_Data{i}(j,:)=SE_Data(Elem_num*(j-1)+i,2:7);        %读取应变数据    Stress_Data{i}(j,:)=SE_Data(Elem_num*(j-1)+i,8:13);       %读取应力数据  endend%读取待评估单元的单元号Elem_eval_id=SE_Data(1:Elem_num,1);%**************************************************************************
%*****************************计算SWT参数**********************************for i=1:Elem_num%读取当前单元中心的应力和应变分量在一个完整循环的变化历程S11=Stress_Data{i}(:,1);                                    %应力分量S22=Stress_Data{i}(:,2);S33=Stress_Data{i}(:,3);S12=Stress_Data{i}(:,4);S13=Stress_Data{i}(:,5);S23=Stress_Data{i}(:,6);E11=Strain_Data{i}(:,1);                                    %应变分量E22=Strain_Data{i}(:,2);E33=Strain_Data{i}(:,3);E12=Strain_Data{i}(:,4);E13=Strain_Data{i}(:,5);E23=Strain_Data{i}(:,6);%计算不同平面对应的SWT值SWT_Plane=zeros(Theta_num,Phi_num);                  %存放每个平面的SWT取值for j=1:Theta_num                                            %遍历角度theta  for k=1:Phi_num                                            %遍历角度phi    %计算旋转平面之后的应力和应变分量在一个完整循环的变化历程    S11_trans=cos(Theta(j))^2*sin(Phi(k))^2*S11+...      2*sin(Theta(j))*cos(Theta(j))*sin(Phi(k))^2*S12+...      2*cos(Theta(j))*sin(Phi(k))*cos(Phi(k))*S13+...      sin(Theta(j))^2*sin(Phi(k))^2*S22+...      2*sin(Theta(j))*sin(Phi(k))*cos(Phi(k))*S23+...      cos(Phi(k))^2*S33;                                %变换之后的应力分量    E11_trans=cos(Theta(j))^2*sin(Phi(k))^2*E11+...      2*sin(Theta(j))*cos(Theta(j))*sin(Phi(k))^2*E12+...      2*cos(Theta(j))*sin(Phi(k))*cos(Phi(k))*E13+...      sin(Theta(j))^2*sin(Phi(k))^2*E22+...      2*sin(Theta(j))*sin(Phi(k))*cos(Phi(k))*E23+...      cos(Phi(k))^2*E33;                                %变换之后的应变分量    %计算当前平面下的法向应变幅值    Delta_Eps_normal=(max(E11_trans)-min(E11_trans))/2;    %计算当前平面下的最**向应力    Sigma_normal_max=max(abs(S11_trans));    %计算当前平面下的SWT参量    SWT_Plane(j,k)=Sigma_normal_max*Delta_Eps_normal;     endend%计算当前待评估点的SWT及对应的临界角度[SWT_Plane_Column,Theta_max_index]=max(SWT_Plane);[SWT_Max,Phi_max_index]=max(SWT_Plane_Column);SWT(i)=SWT_Max;Theta_Critical(i)=Theta(Theta_max_index(Phi_max_index));Phi_Critical(i)=Phi(Phi_max_index);end%**************************************************************************
%**************************重新构造待评估单元定义信息**********************%读取节点信息数据Node_Data=importdata("Node_id.txt");%读取单元定义信息Elem_Def_Data=importdata("Elem_id.txt");Elem_eval_connect=zeros(Elem_num,5);              %用于存放新的单元定义信息%读取待评估单元在单元定义信息中的编号for i=1:Elem_num  %存放待评估单元号至新的单元定义信息  Elem_eval_connect(i,1)=Elem_eval_id(i);  Elem_eval_index=find(Elem_Def_Data(:,1)==Elem_eval_id(i));  k=2;  for j=1:8    %读取待评估单元的节点编号    Elem_eval_node_id=Elem_Def_Data(Elem_eval_index,j+1);    %根据节点编号查询该节点在节点定义信息中的索引    Elem_eval_node_index=find(Node_Data(:,1)==Elem_eval_node_id);     %读取待评估单元节点坐标    Elem_eval_node_coord_X=Node_Data(Elem_eval_node_index,2);    Elem_eval_node_coord_Y=Node_Data(Elem_eval_node_index,3);    Elem_eval_node_coord_Z=Node_Data(Elem_eval_node_index,4);    %根据Z向坐标判断是否需要剔除(!!!!!!!!!根据不同程序可能需要更改)    if Elem_eval_node_coord_Y==0      Elem_eval_connect(i,k)=Node_Data(Elem_eval_node_index,1);      k=k+1;    else      continue;    end  endend%**************************************************************************
%*************************重新整理待评估单元节点信息***********************%读取待评估单元的所有节点号Node_eval_all=zeros(4*Elem_num,1);      %所有待评估单元节点号(含重复节点号)for i=1:Elem_num  for j=1:4    Node_eval_all(4*(i-1)+j)=Elem_eval_connect(i,j+1);  endend%剔除重复节点并对节点编号重新排序Node_eval_all_unique=unique(Node_eval_all);Node_eval_num=size(Node_eval_all_unique,1);           %待评估单元的节点总数%排序后的节点编号重组形成新的节点信息表Node_eval_list=zeros(Node_eval_num,4);Node_eval_list(:,1)=Node_eval_all_unique;%重新读取节点坐标信息至新的节点信息表for i=1:Node_eval_num  %查询节点在原节点信息表中的索引  Node_index=find(Node_Data(:,1)==Node_eval_list(i,1));  %根据索引读取节点坐标至新的节点信息表  Node_eval_list(i,2)=Node_Data(Node_index,2);                   %读取X坐标  Node_eval_list(i,3)=Node_Data(Node_index,3);                   %读取Y坐标  Node_eval_list(i,4)=Node_Data(Node_index,4);                   %读取Z坐标end%**************************************************************************
%**************重新基于新节点定义信息的待评估单元定义信息******************Elem_eval_connect_new=zeros(Elem_num,5);Elem_eval_connect_new(:,1)=Elem_eval_connect(:,1);for i=1:Elem_num  for j=2:5    %计算待评估单元节点在新节点定义信息中对应的索引    Node_index=find(Node_eval_list(:,1)==Elem_eval_connect(i,j));    Elem_eval_connect_new(i,j)=Node_index;  endend%**************************************************************************
%****************************绘制SWT分布***********************************%采用patch函数绘制 需将数据整理为patch要求的格式Patch_X=zeros(4,Elem_num);Patch_Y=zeros(4,Elem_num);Patch_Z=zeros(4,Elem_num);for i=1:Elem_num     %针对单元的循环  for j=1:4             %针对单元内部节点的循环    %读取坐标值    Patch_X(j,i)=Node_eval_list(Elem_eval_connect_new(i,j+1),2);    Patch_Y(j,i)=Node_eval_list(Elem_eval_connect_new(i,j+1),3);    Patch_Z(j,i)=Node_eval_list(Elem_eval_connect_new(i,j+1),4);  endend%绘制每个临界平面的法向量Plane_Normal_Vector=zeros(Elem_num,3);                %存放每个平面的法向量Elem_Center_Coord_X=zeros(Elem_num,1);                %存放单元中心坐标Elem_Center_Coord_Y=zeros(Elem_num,1);Elem_Center_Coord_Z=zeros(Elem_num,1);for i=1:Elem_num  %计算单元中心临界平面法向量  Plane_Normal_Vector(i,1)=sin(Phi_Critical(i))*cos(Theta_Critical(i));  Plane_Normal_Vector(i,2)=sin(Phi_Critical(i))*sin(Theta_Critical(i));  Plane_Normal_Vector(i,3)=cos(Phi_Critical(i));  %计算单元中心坐标  Elem_Center_Coord_X(i)=mean(Patch_X(:,i));  Elem_Center_Coord_Y(i)=mean(Patch_Y(:,i));  Elem_Center_Coord_Z(i)=mean(Patch_Z(:,i));endpatch(Patch_X,Patch_Z,SWT);axis equal;hold on;quiver(Elem_Center_Coord_X,Elem_Center_Coord_Z,...  Plane_Normal_Vector(:,1),Plane_Normal_Vector(:,3),1,'b');axis equal;xlim([-6 6]);set(gca,'color','none');%**************************************************************************

%***********************输出Tecplot可读取的数据文件************************File_Tecplot_id=fopen("SWT_Visual.dat","w");%输出基本信息fprintf(File_Tecplot_id,'title=\"SWT Visualization in 2D FEM\"\n');fprintf(File_Tecplot_id,'variables=\"x\",\"z\",\"SWT\",\"Theta\",\"Phi\"\n');fprintf(File_Tecplot_id,"zone N=%d,",Node_eval_num);fprintf(File_Tecplot_id,"E=%d,",Elem_num);fprintf(File_Tecplot_id,"datapacking=block\n");fprintf(File_Tecplot_id,"varlocation=([3,4,5]=cellcentered)\n");fprintf(File_Tecplot_id,"zonetype=fequadrilateral\n");%输出所有节点的坐标for i=1:Node_eval_num  fprintf(File_Tecplot_id,"%f ",Node_eval_list(i,2));endfprintf(File_Tecplot_id,"\n");for i=1:Node_eval_num  fprintf(File_Tecplot_id,"%f ",Node_eval_list(i,4));endfprintf(File_Tecplot_id,"\n");%输出位于单元中心的SWT值for i=1:Elem_num  if i==Elem_num    fprintf(File_Tecplot_id,"%f\n",SWT(i));  else    fprintf(File_Tecplot_id,"%f ",SWT(i));  endend%输出位于单元中的临界Theta角for i=1:Elem_num  if i==Elem_num    fprintf(File_Tecplot_id,"%f\n",Theta_Critical(i));  else    fprintf(File_Tecplot_id,"%f ",Theta_Critical(i));  endend%输出位于单元中的临界Phi角for i=1:Elem_num  if i==Elem_num    fprintf(File_Tecplot_id,"%f\n",Phi_Critical(i));  else    fprintf(File_Tecplot_id,"%f ",Phi_Critical(i));  endend%输出单元节点构成信息for i=1:Elem_num  for j=2:5    if j==5      fprintf(File_Tecplot_id,"%d\n",Elem_eval_connect_new(i,j));    else      fprintf(File_Tecplot_id,"%d ",Elem_eval_connect_new(i,j));    end  endendfclose(File_Tecplot_id);%**************************************************************************