%混凝土圆柱对心受压疲劳
%CSPD
%二维模型
%最新修改时间2022/5/18/15:28
%静力加载
%结果完全可用
clc;clear
%% 创建计算结果保存路径
directpath='F:\fatigue_Gc_013\';
if (~exist(directpath,'dir') )
    mkdir(directpath)        % 若不存在，在当前目录中产生一个子目录
end
%% 
len=101.6d-3; %原模型横向长度
thick=0.8*len; %原模型竖向长度
hei=203.2d-3; %原模型高度
xdv=50;  %x向划分
ydv=100;  %y向划分
dl=len/xdv;  %节点之间水平或者竖直方向间距
allnum=xdv*ydv;  %平面节点总数 p表示plane
dot=zeros(allnum,2);  %将初始节点坐标（默认为零）存入dotp矩阵
n=0; %节点编号初始化
% 材料参数
nt=20000;
dt=1.0d0;
dens = 2200.0;%密度
e1 = 27.12d9;   %弹性模量
pratio = 0.2d0;   %泊松比
delta=3.015*dl;    %近场范围尺寸
admg=zeros(nt,1);
nc=zeros(nt,1);
tao=1.0d0;
r=0.1d0;
einf=0;
fc=28.9624d6;
ft=fc/10.1267;
%Cosserat参数
I = 0.8d0*dens;
G=e1/(2*(1+pratio));%剪切模量
Gc = 0.03d0.*G;
lc= 1.0d-4;
kn=6*e1/(thick*pi*(delta.^3)*(1-pratio));
kt=2*e1*(1-3*pratio)/(thick*pi*delta*(1-pratio^2));
ld=sqrt(Gc*(1-pratio^2)/(2*e1*(1-3*pratio)))*delta;
kr=4*Gc/(thick*pi*delta^2)-16*G*(lc^2)/(thick*pi*delta^4);

%输入键的拉伸和压缩时的参数
GF=0.0422*1.6289; %能量释放率
scr1t=sqrt(((1-2*pratio)*pi*GF)/(e1*delta));%临界的变形,超过这个值就会破坏
fmaxt=scr1t*kn;%%%
fmaxn=fmaxt*fc/ft;
%输入混凝土强度包络线的三个系数
mf=0.011094;
nf=0.09974;
pf=0.10907;

a1=0.1492;
m1=21.7391;
a2=2.7;
m2=9.27;
% m=9;
% ita=1.0;
% weival=unifrnd(0,1,[allnum,1]);
% weibull=ita.*(-log(1-weival)).^(1/m);
load('F:\fatigue_2d_01\weibull01')
% testweibull=sort(weibull);
% plot(testweibull)
% scatter_color(dotall(1:allnum,:),weibull)
%% 建模
%主体
for i=1:xdv
    for j=1:ydv
        n=n+1;
        dot(n,1)=(-1.0d0*len./2.0d0)+(i-1).*dl;
        dot(n,2)=(-1.0d0*hei./2.0d0)+(j-1).*dl;
        %从左下角往右上角
    end
end
allnum = n;
n=0;
%加入垫块和加载块
dlens=len;
dwids=3.0d0*dl;
xbdvs=round(dlens/dl);
ybdvs=round(dwids/dl);
alldots=xbdvs*ybdvs; %每个加载行的节点数
dotup0=zeros(alldots,2);
dotdown0=zeros(alldots,2);
dotup=zeros(alldots,2);
dotdown=zeros(alldots,2);
%将节点横纵坐标按编号存入dotup0矩阵
m=0;
for i=1:xbdvs
    for j=1:ybdvs
        m=m+1;
        dotup0(m,1)=(-1.0d0*len./2.0d0)+(i-1).*dl;
        dotup0(m,2)=max(dot(1:allnum,2))+dl+(j-1).*dl;
    end
end
dotupn0=allnum+m;
m=0;
for i=1:xbdvs
    for j=1:ybdvs
        m=m+1;
        dotdown0(m,1)=(-1.0d0*len./2.0d0)+(i-1).*dl;
        dotdown0(m,2)=min(dot(1:allnum,2))-3*dl+(j-1).*dl;
    end
end
allnum0=dotupn0+m;
%将节点横纵坐标按编号存入dotup矩阵
m=0;
for i=1:xbdvs
    for j=1:ybdvs
        m=m+1;
        dotup(m,1)=(-1.0d0*len./2.0d0)+(i-1).*dl;
        dotup(m,2)=max(dot(1:allnum,2))+dl+3*dl+(j-1).*dl;
    end
end
dotupn=allnum0+m;
m=0;
for i=1:xbdvs
    for j=1:ybdvs
        m=m+1;
        dotdown(m,1)=(-1.0d0*len./2.0d0)+(i-1).*dl;
        dotdown(m,2)=min(dot(1:allnum,2))-6*dl+(j-1).*dl;
    end
end
dotdownn=dotupn+m;
dotall=[dot(:,:);dotup0;dotdown0;dotup;dotdown];  %将所有所要的节点存入一个矩阵中
dotn=1:1:allnum; %内点编号
% weibull=ones(dotdownn,1);
%% 观察整体构型 不用时可以注释掉
% k1=ones(dotdownn,1);
% k1(allnum+1:dotdownn,1)=0.5;
% scatter_color(dotall,k1)
%% 其他参数初始化
area = dl * dl;%面积
vol = area * thick;%体积
idist = 0.0d0;%初始键长
fac = 0.0d0;%体积修正系数
radij = dl / 2.0d0;%点的半径

cnode = 0;%当前的点;
nlength  = 0.0d0;%两点间长度在变形后的长度
dforce1 = 0.0d0;%两点间的近场力分量x向
dforce2 = 0.0d0;%两点间的近场力分量y向
dforce3 = 0.0d0;

nodefam=knn(dotall,delta,dl);%总的近场范围点的位置矩阵构造,用的是knn函数
%第一列是统计一共有多少个近场点的，所以循环都应该从第二个点开始
nhorizon=max(nodefam(:,1));%读取nodefam中的列数
fncst=ones(dotdownn,nhorizon);%表面修正系数
fail=zeros(dotdownn,nhorizon);%损伤  1代表没破坏，0代表破坏了
lamda=zeros(dotdownn,nhorizon);%键的剩余寿命 1代表完好，0代表破坏
flag=zeros(dotdownn,nhorizon); %判断每一根键的阶段 1代表1阶段，2代表2阶段
eps=zeros(dotdownn,nhorizon); %每根键的应变
sn=zeros(dotdownn,nhorizon); %每根键的法向作用
st=zeros(dotdownn,nhorizon);%每根键的切向作用
switcht=zeros(dotdownn,nhorizon); %判定过程开关,0表示疲劳模拟过程，1表示裂纹扩展过程
aaa=zeros(dotdownn,nhorizon);
%损伤矩阵及键的寿命矩阵初始化
numbond=0;
for i=1:size(fail,1)
    fail(i,1:nodefam(i,1))=1;
    lamda(i,1:nodefam(i,1))=1;
    flag(i,1:nodefam(i,1))=1;
    numbond=numbond+1;
end
basepoint=length(find(fail==1));
%下述分别为损伤矩阵，力，能量，位移，速度，加速度的初始化
dmg=zeros(dotdownn,1);%损伤
pforce=zeros(dotdownn,3);%PD力
bforce=pforce;%外力
disp=pforce;%位移
vel=pforce;%速度
acc=vel;%加速度
aaamax=0;
%% 初始缺陷
% anum=1:1:allnum;
% P=randperm(size(anum,2));
% Pnum=P(1:100);
% for i=1:100
%     pm=Pnum(i);
%     sbond=round(rand*nodefam(pm,1));
%     for j = 1:sbond%讨论每一点的进场范围内其他点的作用
%         fail(pm,j)=0;
%         fail(j,pm)=0;
%     end
% end
%% ADR算法变量初始化
velhalf=vel;%中间变量
velhalfold=vel;
pforceold=pforce;
massvec=zeros(dotdownn,3); %虚质量
for i = 1:dotdownn%虚拟质量   
    massvec(i,1) = 0.25d0 * dt * dt * (pi * (delta)^2 * thick) * kn / (dl) ;
    massvec(i,2) = 0.25d0 * dt * dt * (pi * (delta)^2 * thick) * kn / (dl) ; 
    massvec(i,3) = (1.0d3)*0.25d0 * dt * dt * (pi * (delta)^2 * thick) * kn / (dl) ;
end
flag2point=0;
%% 开始积分求解
% K=zeros(nt,4);
% F1=zeros(nt,4);
% F=zeros(nt,4);
% P1=zeros(nt,4);
% P2=zeros(nt,4);
% sigmaxx=zeros(nt,1);
% sigmaxy=zeros(nt,1);
% sigmayy=zeros(nt,1);
% alldmg=zeros(nt,1);
appbf = (0.198*fc*thick*len/3/(alldots*vol)); %面力转体力
sbond=0.0d0; %只要有键的临界伸长率超过scr0，sbond就不为0
testbond=0.0d0;
timebond=zeros(nt,1); %记录每一次循环后的完好键数
kao=0;
tell=0;
crack=zeros(nt,1);
beta=0.3;

nforce=0;
forceback=zeros(nt,1);
for tt =  1:nt    
    nforce=0;
    fprintf('tt=%d\n',tt);  
    nc(tt,1)=tt/tao;
% 施加荷载与约束
    if tt<=1000
       bforce(allnum0+1:dotupn,2)= -appbf*tt/1000;
       bforce(dotupn+1:dotdownn,2)= appbf*tt/1000;
    else
       bforce(allnum0+1:dotupn,2)= -appbf;
       bforce(dotupn+1:dotdownn,2)= appbf;
    end
    
    for i = 1:allnum0
        if (dotall(i,2)>=max(dotall(1:allnum0,2)))
            disp(i,1)=0;
        end
    end
    for i = 1:allnum0
        if (dotall(i,2)<=min(dotall(1:allnum0,2)))
            disp(i,1)=0;
        end
    end      
    for i = 1:dotdownn  %对内部的每一个点都做循环
        dmgpar1 = 0.0d0;%损伤变量初始化
        dmgpar2 = 0.0d0;
        pforce(i,1) = 0.0d0;%近场力初始化
        pforce(i,2) = 0.0d0;
        pforce(i,3) = 0.0d0;
        for j = 1:nodefam(i,1)%讨论每一点的进场范围内其他点的作用
            cnode = nodefam(i,j+1);
            idist = sqrt((dotall(cnode,1) - dotall(i,1))^2 + (dotall(cnode,2) - dotall(i,2))^2);
            nlength = sqrt((dotall(cnode,1) + disp(cnode,1) - dotall(i,1) - disp(i,1))^2 + (dotall(cnode,2) + disp(cnode,2) - dotall(i,2) - disp(i,2))^2);
            eps(i,j)=(1-r)*((nlength - idist) / idist);
            Abond=pi*delta^2*thick/(464*idist);
            %体积修正
            if (idist<=delta-radij)
                fac = 1.0d0;
            elseif (idist<=delta+radij)
                fac = (delta+radij-idist)/(2.0d0*radij);
            else
                fac = 0.0d0;
            end
            if (abs(dotall(cnode,2) - dotall(i,2))<=1.0d-10)
                theta = 0.0d0;
            elseif (abs(dotall(cnode,1) - dotall(i,1))<=1.0d-10)
                theta = 90.0d0 * pi / 180.0d0;
            else
                theta = atan(abs(dotall(cnode,2) - dotall(i,2)) / abs(dotall(cnode,1) - dotall(i,1)));
            end
            
            fncst(i,j)=2*nhorizon/(nodefam(i,1)+nodefam(cnode,1));
            scr=fncst(i,j);
            
            theat1=(dotall(cnode,1)+disp(cnode,1)-dotall(i,1)-disp(i,1))/nlength;%键的角度值，cos
            theat2=(dotall(cnode,2)+disp(cnode,2)-dotall(i,2)-disp(i,2))/nlength;%键的角度值，sin
              k11=((6*e1)./(thick*pi*idist*(1-pratio)*(delta.^3)));
              k22=(((2*e1*(1-3*pratio))./(thick*pi*delta*(idist.^3)*(1-pratio.^2))));
              k23=(-(sqrt((2*Gc*e1*(1-3*pratio))./(1-pratio.^2))./(thick*pi*(idist.^3))));
              k24=(-(sqrt((2*Gc*e1*(1-3*pratio))./(1-pratio.^2))./(thick*pi*(idist.^3))));
              k32=(-(sqrt((2*Gc*e1*(1-3*pratio))./(1-pratio.^2))./(thick*pi*(idist.^3))));
              k33=(((Gc*delta)./(thick*pi*(idist.^3))-(4*Gc)./(thick*pi*(delta.^2))+(16*G*lc.^2)./(thick*pi*delta.^4)));
              k34=(((Gc*delta)./(thick*pi*(idist.^3))-(4*Gc)./(thick*pi*(delta.^2))-(16*G*lc.^2)./(thick*pi*delta.^4)));           
              v1=(disp(cnode,1)-disp(i,1));
              v2=(disp(cnode,2)-disp(i,2));
              v3=disp(i,3);
              v4=disp(cnode,3);
              F1=(k11*theat1.^2+k22*theat2.^2)*v1+(k11*theat1*theat2-k22*theat1*theat2)*v2-k23*theat2*v3-k24*theat2*v4;
              F2=(k11*theat1*theat2-k22*theat1*theat2)*v1+(k11*theat2.^2+k22*theat1.^2)*v2+k23*theat1*v3+k24*theat1*v4;
              F3=k32*(-theat2)*v1+k32*theat1*v2+k33*v3+k34*v4;   
              F1c=k11*v1*theat1+k11*v2*theat2;
              F2c=-k22*v1*theat2+k22*theat1*v2+k23*v3+k24*v4;
              F3c=-k32*v1*theat2+k32*theat1*v2+k33*v3+k34*v4;
            if (fail(i,j)==1)%如果这两点间的键没有坏，则有近场力作用
               dforce1 = F1 * vol * fac * scr;
               dforce2 = F2 * vol * fac * scr;
               dforce3 = F3 * vol * fac;
               dc1 = F1c * fac * scr;
               dc2 = F2c * fac * scr;
               dc3 = F3c * fac;                      
            else%否则没有
                dforce1 = 0.0d0;
                dforce2 = 0.0d0;
                dforce3 = 0.0d0;
                dc1 = 0.0d0;
                dc2 = 0.0d0;
                dc3 = 0.0d0;             
            end

            % 预计算迭代因子aaa
            if i<=allnum&&cnode<=allnum
               sn(i,j)=sqrt(dc1^2)/fmaxn;
               st(i,j)=beta*sqrt(dc2^2+(dc3/(idist))^2)/fmaxn; 
            end

            if sn(i,j)==0&&st(i,j)==0%表示处于坐标原点
               aaa(i,j)=0;
            elseif sn(i,j)==0&&st(i,j)~=0% sigma n =0,按纯剪切考虑
               aaa(i,j)=st(i,j)/sqrt(mf);
            elseif eps(i,j)>0&&st(i,j)==0% sigma t =0, 按纯拉考虑
               aaa(i,j)=-sn(i,j)/((nf/2/pf)-sqrt((nf/2/pf)^2+mf/pf));
            elseif eps(i,j)<=0&&st(i,j)==0% sigma t =0, 按纯压考虑
               aaa(i,j)=sn(i,j)/((nf/2/pf)+sqrt((nf/2/pf)^2+mf/pf));
            elseif eps(i,j)>0&&st(i,j)~=0&&sn(i,j)~=0%受拉为负
               k=-st(i,j)/sn(i,j);
               stmax=0.5*k*nf/(k^2+pf)+sqrt((0.5*k*nf/(k^2+pf)).^2+mf*k^2/(k^2+pf));
               snmax=stmax/k;
               aaa(i,j)=-sn(i,j)/snmax;
            elseif eps(i,j)<=0&&st(i,j)~=0&&sn(i,j)~=0%受压为正
               k=st(i,j)/sn(i,j);
               stmax=0.5*k*nf/(k^2+pf)+sqrt((0.5*k*nf/(k^2+pf)).^2+mf*k^2/(k^2+pf));
               snmax=stmax/k;
               aaa(i,j)=sn(i,j)/snmax;                                 
            end

            if aaamax<aaa(i,j)
               aaamax=aaa(i,j);
               aaamaxi=i;
               aaamaxj=j;
            end 
            
            
            %求力矩阵，将每个点的作用力叠加起来
            pforce(i,1) = pforce(i,1) + dforce1;
            pforce(i,2) = pforce(i,2) + dforce2;
            pforce(i,3) = pforce(i,3) + dforce3;
            
           %进行键的伸长判断，键伸长达到一定值发生破坏,不去判断那些已经损伤的键     
%            if tt>100&&(fail(i,j)==1)&&i<=allnum&&eps(i,j)>=scr1t&&nof(i)>0            
%                  fail(i,j)=0;%0代表键已经坏了
% %                  flag(i,j)=0;
%                  switcht(i,j)=1.0d0;   %只要有键被拉坏，则跳过所有疲劳过程
%            end
%            进行键的寿命判断，键的寿命在疲劳过程中进行迭代
        if tt>=1000&&(fail(i,j)==1)&&(i<=allnum)&&(cnode<=allnum)&&switcht(i,j)<=0.5   %判定是否还进行疲劳迭代过程
           if dmg(i)>0.398||dmg(cnode)>0.398  %有问题
              for k=1:nodefam(i,1)
                 if (flag(i,k)>0)&&(flag(i,k)<1.5)%1阶段的键进入循环
                     lamda(i,k)=1.0;
                     flag(i,k)=2.0;%这些键进入2阶段
                 end   
              end   
           end
        end
       end
    end
    for i = 1:dotdownn  %对内部的每一个点都做循环
        dmgpar1 = 0.0d0;%损伤变量初始化
        dmgpar2 = 0.0d0;
      for j = 1:nodefam(i,1)%讨论每一点的进场范围内其他点的作用
          cnode = nodefam(i,j+1);
          if (cnode<=allnum)&&(i<=allnum)
              avalue=0.5*(weibull(i)+weibull(cnode));
          else
              avalue=1.0;
          end     
          
        if tt>=1000&&(flag(i,j)>0)&&(flag(i,j)<1.5)&&i<=allnum
            lamda(i,j)=lamda(i,j)-(nc(tt,1)-nc(tt-1,1))*a1*((aaa(i,j)/avalue)^m1);           
        elseif flag(i,j)>1.5 %二阶段的键进行二阶段迭代
            lamda(i,j)=lamda(i,j)-(nc(tt,1)-nc(tt-1,1))*a2*((aaa(i,j)/avalue)^m2);                
        end
              
        if (lamda(i,j)<=0)&&(fail(i,j)==1)
            fail(i,j)=0;
            flag(i,j)=0;
        end
         dmgpar1 = dmgpar1 + fail(i,j) *fac* vol ;%综合起来即为一点处的损伤值，通过损伤值的云图分布来反映具体的裂纹现象
         dmgpar2 = dmgpar2 + vol*fac  ;     
       end
        %计算损伤因子
        dmg(i,1) = 1.0d0 - dmgpar1 / dmgpar2;%损伤因子越大，损伤越大
        if i>=dotupn
           nforce=nforce+vol*pforce(i,2);
        end
    end    
%% 找裂纹长度的代码
    snumber=0;
    jud_cl=zeros(10000,1);
    for ii=1:allnum
        if dotall(ii,1) <= max(dotall(1:allnum,1))-6*dl
           if (dmg(ii,1)>0.6)
              snumber=snumber+1;
              jud_cl(snumber,1)=ii;
           end
        end
    end
    jud_cl=jud_cl(jud_cl~=0);
     if (tt==1)
         crack(tt,1)=0;
     end
     if (tt~=1)
         if(max(jud_cl)~=0)
             crack(tt,1)=max(sqrt((dotall(jud_cl,1)-max(dotall(1:allnum,1))).^2+(dotall(jud_cl,2)-max(dotall(1:allnum,2))).^2));
         else
             crack(tt,1)=0.0d0;
         end
     end        
    %% ADR计算部分
    cn1 = 0.0d0;
    cn2 = 0.0d0;
    
    for i = 1:allnum
        if (velhalfold(i,1)~=0.0d0)
            cn1 = cn1 - disp(i,1) * disp(i,1) * (pforce(i,1) / massvec(i,1) - pforceold(i,1) / massvec(i,1)) / (dt * velhalfold(i,1));
        end
        if (velhalfold(i,2)~=0.0d0)
            cn1 = cn1 - disp(i,2) * disp(i,2) * (pforce(i,2) / massvec(i,2) - pforceold(i,2) / massvec(i,2)) / (dt * velhalfold(i,2));
        end
        if (velhalfold(i,3)~=0.0d0)
            cn1 = cn1 - disp(i,3) * disp(i,3) * (pforce(i,3) / massvec(i,3) - pforceold(i,3) / massvec(i,3)) / (dt * velhalfold(i,3));
        end      
        cn2 = cn2 + disp(i,1) * disp(i,1);
        cn2 = cn2 + disp(i,2) * disp(i,2);
        cn2 = cn2 + disp(i,3) * disp(i,3);
    end
    
    if (cn2~=0.0d0)
        if ((cn1 / cn2) > 0.0d0)
            cn = 2.0d0 * sqrt(cn1 / cn2);
        else
            cn = 0.0d0;
        end
    else
        cn = 0.0d0;
    end
    
    if (cn > 2.0d0)
        cn = 1.9d0;
    end
    
    for i = 1:dotdownn
        if (tt==1)
            velhalf(i,1) = 1.0d0 * dt / massvec(i,1) * (pforce(i,1) + bforce(i,1)) / 2.0d0	;
            velhalf(i,2) = 1.0d0 * dt / massvec(i,2) * (pforce(i,2) + bforce(i,2)) / 2.0d0;           
            velhalf(i,3) = 1.0d0 * dt / massvec(i,3) * (pforce(i,3) + bforce(i,3)) / 2.0d0;
        else
            velhalf(i,1) = ((2.0d0 - cn * dt) * velhalfold(i,1) + 2.0d0 * dt / massvec(i,1) * (pforce(i,1) + bforce(i,1))) / (2.0d0 + cn * dt);
            velhalf(i,2) = ((2.0d0 - cn * dt) * velhalfold(i,2) + 2.0d0 * dt / massvec(i,2) * (pforce(i,2) + bforce(i,2))) / (2.0d0 + cn * dt);            
            velhalf(i,3) = ((2.0d0 - cn * dt) * velhalfold(i,3) + 2.0d0 * dt / massvec(i,3) * (pforce(i,3) + bforce(i,3))) / (2.0d0 + cn * dt);
        end        
        vel(i,1) = 0.5d0 * (velhalfold(i,1) + velhalf(i,1))	;
        vel(i,2) = 0.5d0 * (velhalfold(i,2) + velhalf(i,2))	; 
        vel(i,3) = 0.5d0 * (velhalfold(i,3) + velhalf(i,3))	; 
        disp(i,1) = disp(i,1) + velhalf(i,1) * dt	;
        disp(i,2) = disp(i,2) + velhalf(i,2) * dt	;  
        disp(i,3) = disp(i,3) + velhalf(i,3) * dt	;
        velhalfold(i,1) = velhalf(i,1);
        velhalfold(i,2) = velhalf(i,2);   
        velhalfold(i,3) = velhalf(i,3);
        pforceold(i,1) = pforce(i,1);
        pforceold(i,2) = pforce(i,2);  
        pforceold(i,3) = pforce(i,3); 
    end
    % 基本数据输出，损伤分布以及位移分布
    if(mod(tt,10)==0)%只要时间步长是10的倍数，就进行输出，输出的文件里面有各点的位移和损伤值，若想输出其他的可自己定义。         
         str1=[directpath,'dispdmg',num2str(tt),'.mat'];
         dispdmg=[disp(dotn,:),dmg(dotn,1)];
         save(str1,'dispdmg')     
%          subplot(2,1,1)
%          scatter_color(dotall(1:allnum,:),dispdmg(:,2))     
%          subplot(2,1,2)
%          scatter_color(dotall(1:allnum,:),dispdmg(:,4))         
    end
   timebond(tt,1)=numbond;
   forceback(tt,1)=forceback(tt,1)+nforce;
end
save([directpath,'forceback013',],'forceback');