几何非线性杆单元，弧长法求解程序
本帖最后由 xiaoqiushui318 于 2010-9-28 09:01 编辑%
% program for geometric nonlinear truss element
%
clear
colordef white
tic;
% referenc:
% 1 programing the finite element method with matlab
% by Jack Chessa
% Northwestern University
% 2 finite element proceduresby K.J. Bathe
% the example considered P494-495 Example 6.3
% the formulation can be found on page 542-547
% 3 nonlinear finite element analysis of solids and structures
% by M.A. Criesfied
% the algorithm for arc-length method can be found in Chapter 9.3
%
%***************************************************
% INPUT
%***************************************************
disp('*********************************************')
disp(' Running Start')
disp('*********************************************')
disp()
% material property
E0=100000; %
nu0=0.3; %
% **************************************************
% generate finite element mesh
% **************************************************
% data structure
% node - is a matrix of the node coordinates, i.e. node(I,j) -> x_Ij
% element - is a matrix of element connectivities, i.e. the connectivity
% of element e is given by > element(e,:) ->
%
disp()
ncoord = 3;
l = 10; %length of the element
a = 0.5; %section aera
%
% nodal coordinates
%
node=[0.0 0.0 0.0;
l*cos(pi/12) l*sin(pi/12) 0.0;
2.0*l*cos(pi/12) 0.0 0.0 ];
%
% element connectivity
%
element=[12;
23 ];
numnode = size(node,1); % number of nodes
numelem = size(element,1); % number of elements
enode = 2;
edof = ncoord*enode; % number of dofs per element
% generate mesh
clf
for lmn = 1:numelem
sctr = element(lmn,:);
xcor = node(sctr,1);
ycor = node(sctr,2);
zcor = node(sctr,3);
axis equal;
axis();
plot(xcor,ycor,'-o','LineWidth',2,'MarkerEdgeColor','m','MarkerFaceColor','c','MarkerSize',10);
hold on
end
% set equilibrium iteration controlling data
nstep = 30;
tol = 0.0001;
maxit = 30; % maximum number of iterations
% nflg convergence indicator
% = 0 not converged
% = 1 converged
nflg = 0;
numdof = ncoord*numnode; % total number of dofs
ft = zeros(numdof,1); % total applied external load
ft(5) = -2000.;
kg = zeros(numdof); % global stiffness matrix
ke = zeros(ncoord*enode); % element stiffness matrix
fu = zeros(numdof,1); % unbalanced load vector
%react = zeros(numdof,1); % reaction force
u = zeros(numdof,1); % total displacement vector
delu = zeros(numdof,1); % incremental displacement vector
facto = 1.0/nstep; % incremental loading factor for load control
tfact = 0.0;
SPK2 = zeros(numelem,1); % local element P-K stress
fe = zeros(ncoord*enode,1);% element internal force
fixedNode = '; % a vector of the node numbers which are fixed
fixedNodeX = fixedNode;
fixedNodeY = fixedNode + numnode;
fixedNodeZ = ';
udofs = fixedNodeX; % global indecies of the fixed x displacements
vdofs = fixedNodeY; % global indecies of the fixed y displacements
wdofs = fixedNodeZ + numnode*2; % all dofs in z direction are fixed
xlgic = ones(numdof,1); % xlgic is a logic vector indicating the active dofs
% xlgic(i) = 1 --- i'th dof is active
xlgic(udofs) = 0;
xlgic(vdofs) = 0;
xlgic(wdofs) = 0;
xlgic = logical(xlgic);
uthy = ;
fthy = ;
uout = ;
fstep = ;
itnum = [];
slcol = [];
% ********************************************************************
% set arclength parameter
% ********************************************************************
lamda = 0.0;% load factor
dlamda = 0; % increment of lamda
slmax = 0.5;% maximum arc length allowed
slmin = 0.1;% minimum arc length aloowed
sl = 0.2; % arc length used
for kstep = 1:nstep
fu = lamda*ft; % currently applied external load
iter = 0;
ucorr = zeros(numdof,1);
nflg = 0; % convergence indicator
% evaluate the tangent stiffness matrix
for iele = 1:numelem
sctr = element(iele,:); % element scatter vector
sctrB=; % vector that scatters a B matrix
nn=length(sctr);
dHdr = [ -1/2 1/2 0 0 0 0;
0 0 -1/2 1/2 0 0;
0 0 0 0 -1/2 1/2];
x0 = node(sctr,:);
xv = x0(2,:) - x0(1,:);
L0 = xv(1,:)*xv(1,:)';
L0 = sqrt(L0);
drds = 2/L0;
x0_hat = ';
ut_hat = ';
BL0t = drds*drds*(x0_hat'*dHdr'*dHdr + ut_hat'*dHdr'*dHdr); % linear
BNL0t = drds *dHdr; % nonlinear
ke = BL0t'*E0*BL0t*L0*a + BNL0t'*SPK2(iele)*BNL0t*L0*a;
fe = BL0t'*SPK2(iele)*L0*a;
kg(sctrB,sctrB) = kg(sctrB,sctrB) + ke;
fu(sctrB) = fu(sctrB) - fe;
end
while((nflg==0) && (iter<maxit))
iter = iter + 1;
% modify the golobal stiffness and force vector to accommodate boudnary condtitions
kg(udofs,:)=0; % zero out the rows and columns of the K matrix
kg(vdofs,:)=0;
kg(wdofs,:)=0;
kg(:,udofs)=0;
kg(:,vdofs)=0;
kg(:,wdofs)=0;
kg(udofs,udofs) = speye(length(udofs));
kg(vdofs,vdofs) = speye(length(vdofs));
kg(wdofs,wdofs) = speye(length(wdofs));
fu(udofs) = 0;
fu(vdofs) = 0;
fu(wdofs) = 0;
% evaluate the displacement increment due to unbalanced force
delu = kg\fu;
% evaluate the displacement due to total external force
u_bar = kg\ft;
% update the displacement vector
% u = u + delu;
if iter==1
iii = 0;
if det(kg)>0
dlamda = sl/sqrt(u_bar'*u_bar);
else
dlamda = -sl/sqrt(u_bar'*u_bar);
end
% dl = u_bar'*u_bar;
end
% determine the incremental lamda
%
if iter~=1
c1 = u_bar(xlgic)'*u_bar(xlgic);
c2 = 2*(ucorr(xlgic)+delu(xlgic))'*u_bar(xlgic);
c3 = (ucorr(xlgic)+delu(xlgic))'*(ucorr(xlgic)+delu(xlgic))-sl^2;
c4 = c2^2-4*c1*c3;
if c1 == 0 && c2 ~= 0
dlamda = -c3/c2;
disp('only one root find');
elseif c1 ~=0 && c4>=0
c4 = sqrt(c4);
dlamda1 = (-c2 + c4)/(2*c1);
dlamda2 = (-c2 - c4)/(2*c1);
d1 = ucorr(xlgic)'*(ucorr(xlgic)+delu(xlgic)+dlamda1*u_bar(xlgic));
d2 = ucorr(xlgic)'*(ucorr(xlgic)+delu(xlgic)+dlamda2*u_bar(xlgic));
dlamda = dlamda1;
if d1<d2
dlamda =dlamda2;
end
elseif c1~=0 && c4<0
disp('no real root find');
exit;
end
end
lamda = lamda + dlamda;
delu = delu + dlamda*u_bar;
u = u + delu;
ucorr = ucorr+delu;
fu = lamda*ft;
% loop over all elements to evaluate the element internal force
for iele = 1:numelem
sctr = element(iele,:); % element scatter vector
sctrB=; % vector that scatters a B matrix
nn=length(sctr);
dHdr = [ -1/2 1/2 0 0 0 0;
0 0 -1/2 1/2 0 0;
0 0 0 0 -1/2 1/2];
x0 = node(sctr,:);
xv0 = x0(2,:) - x0(1,:);
L0 = dot(xv0,xv0);
L0 = sqrt(L0);
drds = 2/L0;
x0_hat = ';
ut_hat = ';
xv = x0(2,:) + ut_hat(2:2:edof)' - (x0(1,:) + ut_hat(1:2:edof)');
Ln = dot(xv,xv);
Ln = sqrt(Ln);
delL = Ln - L0;
Epsilon0 = delL/L0;
Sigama0 = E0*Epsilon0;
SPK2(iele) = Sigama0;
BL0t = drds*drds*(x0_hat'*dHdr'*dHdr + ut_hat'*dHdr'*dHdr); % linear
fe = BL0t'*Sigama0*L0*a;
fu(sctrB) = fu(sctrB) - fe;
end
%judge whethere convergence is attained
residu = sum(abs(fu(xlgic)))/sum(abs(ft(xlgic)));
if residu < tol
nflg = 1;
end
%evaluate the geometric stiffness matrix
kg = zeros(numdof);
%loop over all elements to evaluate the element geometric stiffness matrix
for iele = 1:numelem
sctr = element(iele,:); % element scatter vector
sctrB=; % vector that scatters a B matrix
x0 = node(sctr,:);
xv = x0(2,:) - x0(1,:);
L0 = dot(xv,xv);
L0 = sqrt(L0);
drds = 2/L0;
dHdr = [ -1/2 1/2 0 0 0 0;
0 0 -1/2 1/2 0 0;
0 0 0 0 -1/2 1/2];
x0_hat = ';
ut_hat = ';
Sigama0 = SPK2(iele);
BL0t = drds*drds*(x0_hat'*dHdr'*dHdr + ut_hat'*dHdr'*dHdr); % linear
BNL0t = drds *dHdr; % nonlinear
ke = BL0t'*E0*BL0t*L0*a + BNL0t'*Sigama0*BNL0t*L0*a;
kg(sctrB,sctrB) = kg(sctrB,sctrB) + ke;
end
%evaluation geometric stiffness matrix finished
end
%convergence attained at current step
disp([' step ', num2str(kstep), ' passed '])
uout = ;
fstep = ;
% fstep = ;
itnum = ;
itero = iter;
uo = u;
lamdao = lamda;
if sl>slmax
sl = slmax;
elseif sl<slmin
sl = slmin;
end
slcol = ;
% %evaluate analytical solution according to finite element procedures K.J.
% %Bathe P 495
% deltau = abs(u(5));
% xxt = (-1. + 1/sqrt(1-2*deltau/Ln*sin(pi/12)+(deltau/Ln)^2))*(sin(pi/12)-deltau/Ln);
% fthy = ;
% % xxt = (-1. + 1/sqrt(1-2*delL/Ln*sin(pi/12)+(delL/Ln)^2))*(sin(pi/12)-delL/Ln);
% xxt = u(5)/Ln;
% uthy = ;
end
uout = -uout';
fstep = fstep';
itnum = itnum';
fthy = fthy';
uthy = -uthy';
slcol = slcol';
figure
hold on
plot(uout,fstep,'-ob')
xlabel('displacement')
ylabel('load')
% figure
% hold on
% plot(uthy,fthy,'-or')
% xlabel('displacement')
% ylabel('load')
% 笑脸为冒号加单括号 好的学习资料，谢谢楼主分享！！！ 谢谢啊！很实用！ 虽然现在还有点不懂，但是我认为本论坛的强人还是不少的。感谢分享程序！
页:
[1]