matlab生成粗糙表面到abaqus

function [Nodes,Elements]=Rand_surf(Lx,Ly,Lz,N,sn_x,sn_y,WB_a,WB_b)
%function [Nodes,Elements]=Rand_surf(Lx,Ly,Lz,N,sn_x,sn_y,WB_a,WB_b)
%[Nodes,Elements]=Rand_surf(1,1,0.1,5,30,30,0.05,2.0)  
%复制这行改参数，Lx,Ly,Lz三轴长度
%WB_a,WB_b 威布尔分布参数
%N两个样本点之间参考点数
% m code by Youngbin LIM
% To inquire any queries, Contact: lyb0684@naver.com
%************************%
%Description on the input%
%************************%
% Choose the Length of x and y direction (Lx, Ly)
% Lz is the thickness of deformable solid. It is the distance between 
% the bottom surface and the top (averaged) surface
% Input Lz=0 corresponds to 3D surface
% N is the number of division for distance between sampling points 
% It will decide the mesh size
% The element length in z-direction is approximately the legnth in x,y
% Choose the number of sampling points for x (sn_x) and y (sn_y) direction
% WB_a and WB_b is parameters for Weibull distribution
% where, Cummulative Distribution Function: F=1-exp[-(x/a)^b]
%********************************************%
% Random sampling of Points on Rough surface %
%********************************************%
% Inverse transformation sampling with Weibull distribution is used %
% Rough surface is assumed to be periodic in x & y direction %
% Side walls (4 faces) are generated in rectangular shape %
% It is for periodcity and successful geometry conversion in Abaqus/CAE %
R_sample=-WB_a*gamma(1+1/WB_b)+WB_a*(-log(1-rand(sn_y,sn_x))).^(1/WB_b);
R_sample=[R_sample(:,1)/2 R_sample]; R_sample=[R_sample(1,:)/2; R_sample];
R_sample=[R_sample R_sample(:,size(R_sample,2))/2]; R_sample=[R_sample; R_sample(size(R_sample,1),:)/2];
R_sample=[R_sample R_sample(:,1)-R_sample(:,1)]; R_sample=[R_sample; R_sample(1,:)-R_sample(1,:)];
R_sample=[R_sample(:,1)-R_sample(:,1) R_sample ]; R_sample=[R_sample(1,:)-R_sample(1,:); R_sample];
maxTh=max(R_sample,[],'all')+Lz;
minTh=min(R_sample,[],'all')+Lz;
%************************************%
%Rigid shell generation for Lz=0 case%
%************************************%
if Lz==0 
dx=Lx/(sn_x+3); dy=Ly/(sn_y+3); 
x=0:dx:Lx; y=0:dy:Ly;
x_s=0:dx/N:Lx; y_s=0:dy/N:Ly; x_s=x_s'; y_s=y_s';
%*********************%
%Spline in y-direction%
%*********************%
Ry=[];
for i=1:size(R_sample,2)
    Ry(:,i)=spline(y,R_sample(:,i),y_s);
    Ry=[Ry Ry(:,i)];
end
Ry(:,size(R_sample,2)+1)=[]; 
%***********************%
%Spline for whole points%
%***********************%
R=[];
for i=1:size(Ry,1)
    R(i,:)=spline(x,Ry(i,:),x_s);
    R=[R; R(i,:)];
end
R(size(R,2),:)=[];
%*****************%
%Rearrange R, x, y%
%*****************%
R=R'; z=[]; x=[]; y=[];
for i=1:size(R,2)
    z=[z; R(:,i)];
end
for j=1:size(R,2)
    x=[x;x_s];
end
for j=1:size(R,2)
for i=1:size(x_s)
    y=[y; y_s(j)];
end
end
%*******************************%
%Nodes list for 3D rough surface%
%*******************************%
N_num=1:1:size(z); N_num=N_num';
Nodes=[N_num x y z];
%***************%
%Define Elements%
%***************%
y_elnum=size(y_s,1)-1;x_elnum=size(x_s,1)-1;
El_num=1:1:x_elnum*y_elnum; El_num=El_num';
El_N=zeros(size(El_num,1),4);
%*********************%
%Roughness calculation%
%*********************%
z_avg=mean(Nodes(:,4),'All');
Ra=(1/(Ly*size(y_s,1)))*sum(abs(Nodes(:,4)-z_avg)*(dx/N),'All');
Rq=sqrt((1/(Ly*size(y_s,1)))*sum(((Nodes(:,4)-z_avg).^2)*(dx/N),'All'));
Sa=(1/Lx/Ly)*sum(abs(Nodes(:,4)-z_avg)*(dx/N)*(dy/N),'All');
Sq=sqrt((1/Lx/Ly)*sum(((Nodes(:,4)-z_avg).^2)*(dx/N)*(dy/N),'All'));
fprintf('\nRoughness parameters:')
fprintf('\nRa=%f, Rq(RMS)=%f',Ra,Rq)
fprintf('\nSa=%f, Sq(RMS)=%f',Sa,Sq)
fprintf('\n\nNumber of nodes: %d',size(Nodes,1))
fprintf('\nNumber of elements (R3D4): %d\n\n',size(El_num,1))
if size(El_num,1)>1000000
    fprintf('****************************************************')
    fprintf('\n***WARNING: Number of elements are over 1 million***')
    fprintf('\n****************************************************')
    fprintf('\n\nPress "ENTER" to continue')
    fprintf('\nPress "CTRL+C" to abort\n\n')
    pause;
end
%***************%
%Mesh generation%
%***************%
for k=1:y_elnum
    for j=1:x_elnum
    El_N(j+x_elnum*(k-1),:)=[j+(x_elnum+1)*(k-1) j+(x_elnum+1)*(k-1)+1 j+(x_elnum+1)*(k-1)+(x_elnum+2) j+(x_elnum+1)*(k-1)+(x_elnum+1)];
    end
end
Elements=[El_num El_N]; 
%Element sets for surface definition%
elset_S1=[1 x_elnum*y_elnum];
%****************%
%Write input file%
%****************%
fileID=fopen('Rough_Surf_Rigid_Shell.inp','w');
fprintf(fileID, '*Part, Name=Rough_surf');
fprintf(fileID, '\n*Node');
fprintf(fileID,'\n%d, %f, %f, %f',Nodes');
fprintf(fileID, '\n*Element, type=R3D4');
fprintf(fileID, '\n%d, %d, %d, %d, %d',Elements');
fprintf(fileID, '\n*Elset, elset=elset_S1, generate');
fprintf(fileID,'\n%d, %d, 1',elset_S1);
fprintf(fileID, '\n*Surface, type=element, name=S1\nelset_S1, S1');
fprintf(fileID, '\n*End Part'); 
fclose(fileID);

elseif Lx<0 || Ly<0 || Lz<0
    fprintf('\nThe length value (Lx,Ly,Lz) should be positive number\n\n')
    Nodes=[]; Elements=[];

%**********************************************%
%Check for intersection of top and bottom layer%
%**********************************************%

elseif minTh<0 && Lz>0
    fprintf('\n**********************************************')
    fprintf('\n************ Mesh generation faild ***********')
    fprintf('\n**********************************************')
    fprintf('\nThe bottom layer intersects with the top layer')
    fprintf('\nPlease increase the thickness\n\n')
    Nodes=[];
    Elements=[];

elseif 0<minTh && maxTh/minTh>10 && Lz>0
    fprintf('\n**********************************************')
    fprintf('\n************ Mesh quality warning ************')
    fprintf('\n**********************************************')
    fprintf('\nThe maximum thickness is 10 times greater than the minimum thickness')
    fprintf('\nThere could be mesh quality issue')
    fprintf('\nPlease increase the thickness\n\n')
    Nodes=[];
    Elements=[];

dx=Lx/(sn_x+3); dy=Ly/(sn_y+3);
x_s=0:dx/N:Lx; y_s=0:dy/N:Ly; x_s=x_s'; y_s=y_s';
avg_mesh_size=(dx/N+dy/N)/2;
Bot_layer_num=floor(Lz/avg_mesh_size);
y_elnum=size(y_s,1)-1;x_elnum=size(x_s,1)-1;
El_num=1:1:Bot_layer_num*x_elnum*y_elnum;
else
dx=Lx/(sn_x+3); dy=Ly/(sn_y+3); 
x=0:dx:Lx; y=0:dy:Ly;
x_s=0:dx/N:Lx; y_s=0:dy/N:Ly; x_s=x_s'; y_s=y_s';
%*********************%
%Spline in y-direction%
%*********************%
Ry=[];
for i=1:size(R_sample,2)
    Ry(:,i)=spline(y,R_sample(:,i),y_s);
    Ry=[Ry Ry(:,i)];
end
Ry(:,size(R_sample,2)+1)=[]; 
%***********************%
%Spline for whole points%
%***********************%
R=[];
for i=1:size(Ry,1)
    R(i,:)=spline(x,Ry(i,:),x_s);
    R=[R; R(i,:)];
end
R(size(R,2),:)=[];
%*****************%
%Rearrange R, x, y%
%*****************%
R=R'; z=[]; x=[]; y=[];
for i=1:size(R,2)
    z=[z; R(:,i)];
end
for j=1:size(R,2)
    x=[x;x_s];
end
for j=1:size(R,2)
for i=1:size(x_s)
    y=[y; y_s(j)];
end
end
%********************************%
%Nodes list for top rough surface%
%********************************%
N_num_top=1:1:size(z); N_num_top=N_num_top';
Nodes_top=[N_num_top x y z];
%****************************%
%Nodes list for bottom layers%
%****************************%
avg_mesh_size=(dx/N+dy/N)/2;
Bot_layer_num=floor(Lz/avg_mesh_size); dz=(Lz+z)/Bot_layer_num;
x_bot=[]; y_bot=[]; z_bot=[]; 
for i=1:Bot_layer_num
    z_bot=[z_bot;z-i*dz];
end
for i=1:Bot_layer_num
    x_bot=[x_bot;x];
    y_bot=[y_bot;y];
end
N_num_bot=(size(z)+1):1:(size(z,1)+size(z_bot,1)); N_num_bot=N_num_bot';
Nodes_bot=[N_num_bot x_bot y_bot z_bot];
%******************%
%Combine Nodes list%
%******************%
Nodes=[Nodes_top;Nodes_bot];
%***************%
%Define Elements%
%***************%
y_elnum=size(y_s,1)-1;x_elnum=size(x_s,1)-1;
El_num=1:1:Bot_layer_num*x_elnum*y_elnum; El_num=El_num';
El_N=zeros(size(El_num,1),8);
N_shift=size(y_s,1)*size(x_s,1);% Node number shifting for bottom layer
%*********************%
%Roughness calculation%
%*********************%
z_avg=mean(Nodes_top(:,4),'All');
Ra=(1/(Ly*size(y_s,1)))*sum(abs(Nodes_top(:,4)-z_avg)*(dx/N),'All');
Rq=sqrt((1/(Ly*size(y_s,1)))*sum(((Nodes_top(:,4)-z_avg).^2)*(dx/N),'All'));
Sa=(1/Lx/Ly)*sum(abs(Nodes_top(:,4)-z_avg)*(dx/N)*(dy/N),'All');
Sq=sqrt((1/Lx/Ly)*sum(((Nodes_top(:,4)-z_avg).^2)*(dx/N)*(dy/N),'All'));
fprintf('\nRoughness parameters:')
fprintf('\nRa=%f, Rq(RMS)=%f',Ra,Rq)
fprintf('\nSa=%f, Sq(RMS)=%f',Sa,Sq)
fprintf('\n\nNumber of nodes: %d',size(Nodes,1))
fprintf('\nNumber of elements (C3D8R): %d\n\n',size(El_num,1))
if size(El_num,1)>1000000
    fprintf('****************************************************')
    fprintf('\n***WARNING: Number of elements are over 1 million***')
    fprintf('\n****************************************************')
    fprintf('\n\nPress "ENTER" to continue')
    fprintf('\nPress "CTRL+C" to abort\n\n')
    pause;
end
%****************%
%Mesh generationr%
%****************%
for i=1:Bot_layer_num
    for k=1:y_elnum
        for j=1:x_elnum
         El_N(j+x_elnum*(k-1)+(i-1)*x_elnum*y_elnum,:)=[j+(x_elnum+1)*(k-1)+i*N_shift j+(x_elnum+1)*(k-1)+1+i*N_shift j+(x_elnum+1)*(k-1)+(x_elnum+2)+i*N_shift j+(x_elnum+1)*(k-1)+(x_elnum+1)+i*N_shift j+(x_elnum+1)*(k-1)+(i-1)*N_shift j+(x_elnum+1)*(k-1)+1+(i-1)*N_shift j+(x_elnum+1)*(k-1)+(x_elnum+2)+(i-1)*N_shift j+(x_elnum+1)*(k-1)+(x_elnum+1)+(i-1)*N_shift];
        end
    end
end
Elements=[El_num El_N];
%Element sets for surface definition%
elset_S1=[x_elnum*y_elnum*(Bot_layer_num-1)+1 size(El_num,1)];
elset_S2=[1 x_elnum*y_elnum];
elset_S3=zeros(Bot_layer_num,2);
for i=1:Bot_layer_num
    elset_S3(i,:)=[1+(x_elnum*y_elnum)*(i-1) x_elnum+(x_elnum*y_elnum)*(i-1)];
end
elset_S4=zeros(Bot_layer_num,3);
for i=1:Bot_layer_num
    elset_S4(i,:)=[x_elnum+(x_elnum*y_elnum)*(i-1) x_elnum*y_elnum+(x_elnum*y_elnum)*(i-1) x_elnum];
end
elset_S5=zeros(Bot_layer_num,2);
for i=1:Bot_layer_num
    elset_S5(i,:)=[x_elnum*(y_elnum-1)+1+(x_elnum*y_elnum)*(i-1) x_elnum*y_elnum+(x_elnum*y_elnum)*(i-1)];
end
elset_S6=zeros(Bot_layer_num,3);
for i=1:Bot_layer_num
    elset_S6(i,:)=[1+(x_elnum*y_elnum)*(i-1) x_elnum*(y_elnum-1)+1+(x_elnum*y_elnum)*(i-1) x_elnum];
end
%****************%
%Write input file%
%****************%
fileID=fopen('Rough_Surf_3D_Solid.inp','w');
fprintf(fileID, '*Part, Name=Rough_surf');
fprintf(fileID, '\n*Node');
fprintf(fileID,'\n%d, %f, %f, %f',Nodes');
fprintf(fileID, '\n*Element, type=C3D8R');
fprintf(fileID, '\n%d, %d, %d, %d, %d, %d, %d, %d, %d',Elements');
fprintf(fileID, '\n*Elset, elset=elset_S1, generate');
fprintf(fileID,'\n%d, %d, 1',elset_S1);
fprintf(fileID, '\n*Elset, elset=elset_S2, generate');
fprintf(fileID,'\n%d, %d, 1',elset_S2);
fprintf(fileID, '\n*Elset, elset=elset_S3, generate');
fprintf(fileID,'\n%d, %d, 1',elset_S3');
fprintf(fileID, '\n*Elset, elset=elset_S4, generate');
fprintf(fileID,'\n%d, %d, %d',elset_S4');
fprintf(fileID, '\n*Elset, elset=elset_S5, generate');
fprintf(fileID,'\n%d, %d, 1',elset_S5');
fprintf(fileID, '\n*Elset, elset=elset_S6, generate');
fprintf(fileID,'\n%d, %d, %d',elset_S6');
fprintf(fileID, '\n*Surface, type=element, name=S1\nelset_S1, S1');
fprintf(fileID, '\n*Surface, type=element, name=S2\nelset_S2, S2');
fprintf(fileID, '\n*Surface, type=element, name=S3\nelset_S3, S3');
fprintf(fileID, '\n*Surface, type=element, name=S4\nelset_S4, S4');
fprintf(fileID, '\n*Surface, type=element, name=S5\nelset_S5, S5');
fprintf(fileID, '\n*Surface, type=element, name=S6\nelset_S6, S6');
fprintf(fileID, '\n*End Part'); 
fclose(fileID);
end