%弧长法. 两结点单元. 载荷步参数4.
clc;clear;
tol=0.001; %收敛误差.
factor=0.5; %起始因子.
endfactor=4; %结束因子.
q=4/endfactor; %体力
n=10;%单元数
nDOF=n+1;%自由度数
le=1/n;%单元长度
index=zeros(n,2);%所有单元的自由度索引向量
for i=1:n
index(i,1)=i;
index(i,2)=i+1;
end
FreeDOF=2:1:n+1; %没有约束的自由度(只有自由度1有位移约束)
d=[1;1.5*ones(n,1)];%指定初始位移值
RP=[zeros(n,1);2/endfactor]; %指定载荷向量
fd=zeros(n,1); %
iter=1; %总迭代数
count=1;
fconv=1.0;
oldfactor=0;
while ((abs(factor-endfactor) > tol)||(abs(factor-oldfactor) > tol)||(fconv > tol))
incD=zeros(nDOF,1); %位移增量
KT=zeros(nDOF,nDOF);%总体切线刚度矩阵
RI=zeros(nDOF,1);%总体内力向量
RQ=zeros(nDOF,1);%体力q的总体载荷向量.
R=zeros(nDOF,1);%总体不平衡力向量.
for i=1:n
indexi=index(i,:);%提取单元自由度索引
di=d(indexi); %提取单元位移向量
%迭代i时的切线刚度矩阵kt.
kti=[di(1)^2/le -di(2)^2/le;
-di(1)^2/le di(2)^2/le];
%迭代i时的内力向量ri.
rii=[(di(1)^3-di(2)^3)/(3*le);(-di(1)^3+di(2)^3)/(3*le)];
%迭代i时的载荷向量rq.
rqi=[q*le/2; q*le/2];
%集成单元kt,ri,rq成总体KT,RI,RQ.
KT(indexi,indexi)= KT(indexi,indexi) + kti;
RI(indexi)=RI(indexi)+rii;
RQ(indexi)=RQ(indexi)+rqi;
end
%计算不平衡力向量R.
RE=RQ+RP;
fRE=factor*RE; %载荷因子乘以外载得到当前载荷步载荷向量
R=fRE-RI;
%施加位移边界条件.第一个自由度约束,d(1)为零.
KKT=KT(FreeDOF,FreeDOF);
RR=R(FreeDOF);
RRE=fRE(FreeDOF);
fconv=norm(RR,2)^2/(1+norm(RRE,2)^2);
%如果此载荷步收敛，则进行下一下载荷步.
if(fconv <0.01*tol)
s=s*5/count;
factor=2*factor;
if(factor > endfactor) factor=endfactor;end
fd=zeros(n,1);
if(abs(factor-oldfactor) > tol) fconv=1;end
count=1;
else
%求解a 及 b 向量.
a=inv(KKT)*RR;
b=inv(KKT)*RE(FreeDOF);
%确定dfactor.
%第一次迭代时, 设置s 及 dfactor值.
if(iter == 1)
s=factor*sqrt(b'*b);
dfactor=0.0;
else %否则，用弧长法确定dfactor.
a1=(fd+a)'*(fd+a)-s^2;
a2=2*(fd+a)'*b;
a3=b'*b;
a4=a2^2-4*a1*a3;
if (a4>0)
df1=(-a2+sqrt(a4))/(2*a3);
df2=(-a2-sqrt(a4))/(2*a3);
dd1=a+b*df1;
dd2=a+b*df2;
c1=(fd+dd1)'*fd;
c2=(fd+dd2)'*fd;
if(c1>0 && c2<0)
dfactor=df1;
elseif(c1<0 && c2>0)
dfactor=df2;
elseif(c1<0 && c2<0)
dfactor=0;
else
dflin=-a1/a2;
if(abs(df1-dflin)<abs(df2-dflin))
dfactor=df1;
else
dfactor=df2;
end
end
else
alph=fconv/tol;
if(alph<0.1) alph = 0.1;end
if(alph>0.5) alph = 0.5;end
dfactor=-alph*factor;
end
if(factor+dfactor>endfactor)dfactor=0;end
end
%计算位移增量
dd=a+b*dfactor;
incD(FreeDOF)=incD(FreeDOF)+dd;
%更新位移.
fd=fd+dd;
oldfactor=factor;
factor=factor+dfactor;
d=d+incD;
count=count+1;
iter=iter+1;
end
if(iter>100)break;end
end
%绘近似解与准确解图形.
x=0:0.02:1;
u=(-6*x.^2+18*x+1).^(1/3);
plot(x,u);
hold on;
n=(0:1/n:1)';
plot(n,d,'r--');
legend('exact solution','FEM solution');
xlabel('x');
ylabel('u');
