几何非线性杆单元，弧长法求解程序

xiaoqiushui318

%
% 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 procedures  by 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([num2str(toc),' START'])
% 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, -> [n1 n2 n3 ...]
%
disp([num2str(toc),' GENERATING MESH'])
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=[  1  2;
           2  3 ];
      
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([0 20 -1 4]);
    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 = [1 3]';          % a vector of the node numbers which are fixed
  fixedNodeX = fixedNode;
  fixedNodeY = fixedNode + numnode;
  fixedNodeZ = [1 2 3]';
  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 = [0];
  fthy = [0];
  uout = [0];
  fstep = [0];
  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=[sctr sctr+numnode sctr+numnode+numnode]; % 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 = [x0(1,1) x0(2,1) x0(1,2) x0(2,2) x0(1,3) x0(2,3)]';
        ut_hat = [u(sctr)' u(sctr+numnode)' u(sctr+numnode+numnode)']';
        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=[sctr sctr+numnode sctr+numnode+numnode]; % 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 = [x0(1,1) x0(2,1) x0(1,2) x0(2,2) x0(1,3) x0(2,3)]';
           ut_hat = [u(sctr)' u(sctr+numnode)' u(sctr+numnode+numnode)']';         
           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=[sctr sctr+numnode sctr+numnode+numnode]; % 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 = [x0(1,1) x0(2,1) x0(1,2) x0(2,2) x0(1,3) x0(2,3)]';
           ut_hat = [u(sctr)' u(sctr+numnode)' u(sctr+numnode+numnode)']';         
           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 = [uout u(5)/Ln];
   fstep = [fstep -lamda*ft(5)/(2*E0*a)];
%    fstep = [fstep lamda];
   itnum = [itnum iter];
   itero = iter;   
   uo = u;
   lamdao = lamda;
   
   if sl>slmax
       sl = slmax;
   elseif sl<slmin
       sl = slmin;
   end
   slcol = [slcol sl];
% %  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 = [fthy,xxt];
% %    xxt = (-1. + 1/sqrt(1-2*delL/Ln*sin(pi/12)+(delL/Ln)^2))*(sin(pi/12)-delL/Ln);
%   xxt = u(5)/Ln;
%   uthy = [uthy xxt];
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')
%

xiaoqiushui318

笑脸为冒号加单括号

jsl024110

好的学习资料，谢谢楼主分享！！！

renyz518

谢谢啊！很实用！

zh123168

虽然现在还有点不懂，但是我认为本论坛的强人还是不少的。感谢分享程序！