function PLAST3D(MID,PROP,ETAN,NE,NDOF,XYZ,LE,DYNF) %****************************************************************************** % MAIN PROGRAM COMPUTING GLOBAL STIFFNESS MATRIX RESIDUAL FORCE FOR % % PLASTIC MATERIAL MODELS % %****************************************************************************** %% % global LTAN UPDATE DISPDD DISPTD FORCE GKF SIGMA XQ % Global variables XG=[-1/sqrt(3), 1/sqrt(3)]; WGT=[1.0, 1.0]; % Integration points & weights INTN=0; % Index for history variables (each integration pt) for IE=1:NE % Loop over elements. Assemble K and F ELXY=XYZ(LE(IE,:),:); % Element nodal coordinates IDOF=zeros(1,NDOF*size(LE,2)); % Local to global DOF mapping index for I=1:size(LE,2) % Loop for element nodes IDOF(I*NDOF+[-2 -1 0])=LE(IE,I)*NDOF+[-2 -1 0]; % DOF mapping index end % End of element nodes loop DSP=reshape(DISPTD(IDOF),NDOF,8); % NDOFxNENODE element total displacement DSPD=reshape(DISPDD(IDOF),NDOF,8); % element incremental displacement for LX=1:2, for LY=1:2, for LZ=1:2 %LOOP OVER INTEGRATION POINTS E1=XG(LX); E2=XG(LY); E3=XG(LZ); % Integration point INTN = INTN + 1; % Integration point sequential index [~,SHPD,DET]=SHAPEL([E1 E2 E3], ELXY); % Det & shape function derivatives FAC=WGT(LX)*WGT(LY)*WGT(LZ)*DET; % Integraton weights if MID>20, NALPHA=3; STRESSN=SIGMA(7:12,INTN); % Restore converged stress else, NALPHA=6; STRESSN=SIGMA(1:6,INTN); end % from previous time ALPHAN=XQ(1:NALPHA,INTN); % Restore back stress EPN=XQ(NALPHA+1,INTN); % Restore effective plastic strain if MID>10, F=DSP*SHPD'+eye(3); SHPD=inv(F)'*SHPD; end % Updated Lagrangian DEPS=DSPD*SHPD'; % Incremental strain DDEPS=[DEPS(1,1) DEPS(2,2) DEPS(3,3) ... % Voigt notation DEPS(1,2)+DEPS(2,1) DEPS(2,3)+DEPS(3,2) DEPS(1,3)+DEPS(3,1)]'; % if MID==1 % Infinitesimal plasticity [STRESS,ALPHA,EP,DTAN]=COMBHARD(PROP,ETAN,DDEPS,STRESSN,ALPHAN,EPN,LTAN); % elseif MID==11 % Finite rotation plasticity FAC=FAC*det(F); % updated Lagrangian integration [STRESSN,ALPHAN]=ROTATESTRESS(DEPS,STRESSN,ALPHAN); % Rotate stress [STRESS,ALPHA,EP,DTAN]=COMBHARD(PROP,ETAN,DDEPS,STRESSN,ALPHAN,EPN,LTAN); % elseif MID==21 % Multiplicative plasticity [STRESS,B,ALPHA,EP,DTAN]=MULPLAST(PROP,ETAN,DEPS,STRESSN,ALPHAN,EPN,LTAN);% end % End of plasticity model branch if UPDATE % Update plastic variables SIGMA(1:6,INTN)=STRESS; % Store stress XQ(:,INTN)= [ALPHA; EP]; % Store back stress & effective plastic strain if MID>20, SIGMA(7:12,INTN)=B; end % Store right Cauchy-Green tensor continue; % Do not calculate F and K end % End of update branch BM=zeros(6,24); BG=zeros(9,24); % displacement-strain matrix for I=1:8 % Element nodes loop COL=(I-1)*NDOF+1:(I-1)*NDOF+NDOF; % Node DOFs BM(:,COL)=[SHPD(1,I) 0 0; 0 SHPD(2,I) 0; 0 0 SHPD(3,I); % Displacement- SHPD(2,I) SHPD(1,I) 0; 0 SHPD(3,I) SHPD(2,I); % strain matrix SHPD(3,I) 0 SHPD(1,I)]; % BG(:,COL)=[SHPD(1,I) 0 0; SHPD(2,I) 0 0; SHPD(3,I) 0 0; % Nonlinear 0 SHPD(1,I) 0; 0 SHPD(2,I) 0; 0 SHPD(3,I) 0; % Displacement- 0 0 SHPD(1,I); 0 0 SHPD(2,I); 0 0 SHPD(3,I)]; % strain matrix end % End of element node loop FORCE(IDOF) = FORCE(IDOF) - FAC*BM'*STRESS; % Residual internal forces if LTAN % Calculate tangent stiffness SIG=[STRESS(1) STRESS(4) STRESS(6); % Stress tensor STRESS(4) STRESS(2) STRESS(5); % STRESS(6) STRESS(5) STRESS(3)]; % if MID==1 % Infinitesimal plasticity EKF = BM'*DTAN*BM; % Calculate tangent stiffness elseif MID==11 % Finite rotation plasticity CTAN=[-STRESS(1) STRESS(1) STRESS(1) -STRESS(4) 0 -STRESS(6); % Tangent STRESS(2) -STRESS(2) STRESS(2) -STRESS(4) -STRESS(5) 0; % stiffness STRESS(3) STRESS(3) -STRESS(3) 0 -STRESS(5) -STRESS(6); % from -STRESS(4) -STRESS(4) 0 -0.5*(STRESS(1)+STRESS(2)) ... % objective -0.5*STRESS(6) -0.5*STRESS(5); % stress 0 -STRESS(5) -STRESS(5) -0.5*STRESS(6) ... % rate -0.5*(STRESS(2)+STRESS(3)) -0.5*STRESS(4); % -STRESS(6) 0 -STRESS(6) -0.5*STRESS(5) ... % -0.5*STRESS(4) -0.5*(STRESS(1)+STRESS(3))]; % SHEAD=kron(eye(3),SIG); % Initial stiffness matrix EKF = BM'*(DTAN+CTAN)*BM + BG'*SHEAD*BG; % Calculate tangent stiffness elseif MID==21 % Multiplicative plasticity SHEAD=kron(eye(3),SIG); % Initial stiffness matrix EKF = BM'*DTAN*BM + BG'*SHEAD*BG; % Calculate tangent stiffness end % End of plasticity model branch GKF(IDOF,IDOF)=GKF(IDOF,IDOF)+DYNF*FAC*EKF; % Assemble to global K end % End of LTAN branch end, end, end % End of integration points loop end % End of element loop end %