function NLFEA(PARAM,NOUT,MID,PROP,EXTFORCE,SDISPT,XYZ,LE,CNT,DYN) %****************************************************************************** % MAIN PROGRAM FOR NONLINEAR FINITE ELEMENT ANALYSIS % 2D PLANE-STRAIN QUADRILATERAL ELEMENTS & 3D HEXAHEDRAL ELEMENTS % LINEAR ELASTIC, HYPERELASTIC, ADDITIVE & MULTIPLICATIVE PLASTIC MATERIALS % 2D & 3D FRICTIONAL CONTACT/IMPACT ANALYSIS % IMPLICIT NONLINEAR DYNAMICS 2025 N.H.KIM %****************************************************************************** %% global LTAN UPDATE DISPDD DISPTD ACC FORCE GKF % Global variables [NUMNP,NDOF]=size(XYZ); NE=size(LE,1); NEQ=NDOF*NUMNP; % Problem sizes LCNT=false; LDYN=false; DYNF=1; % Default: no contact & no dynamics if ~isempty(CNT) % Presence of contact nodes LCNT=true; % Logical variable for contact XYZ=[XYZ;CNT.XYZC]; % Augment rigid-body coordinates if CNT.PARM(3)==0, CNT.PARM(2)=0; end % Frictionless contact setup end % End of contact setup if ~isempty(DYN) % Dynamic analysis LDYN=true; % Logical variable for dynamic analysis ACC=zeros(NEQ,1); VEL=ACC; % Initialize the velocity & acceleration if ~isempty(DYN.IV), VEL=repmat(DYN.IV,NUMNP,1); end % Initial velocity end % End of dynamic analysis setup DISPTD=zeros(NEQ,1); DISPDD=zeros(NEQ,1); % Nodal displacement & increment if MID>=0, [PROP,ETAN]=PLSET(PROP,MID,NE,NDOF); end % Initial material prop. CHECKELEM(XYZ,LE,NOUT); % Check element connectivity % Set up load increments------------------------------------------------------- NLOAD=size(PARAM.TIMS,2); ILOAD=1; % # of load cases TIMEF=PARAM.TIMS(1,ILOAD); TIMEI=PARAM.TIMS(2,ILOAD); % Start & end time CUR1=PARAM.TIMS(4,ILOAD); CUR2=PARAM.TIMS(5,ILOAD); % Start & end load factor DELTA=PARAM.TIMS(3,ILOAD); DELTA0=DELTA; % Time increment & save TIME=TIMEF; TDELTA=TIMEI-TIMEF; % Current time & time interval for load step ITOL=1; TARY=zeros(PARAM.NTOL,1); % Bisection level & time stamps % Load increment loop---------------------------------------------------------- ISTEP=-1; FLAG10=1; % Initial step and convergence flag while FLAG10 % Solution has been converged FLAG10=0; FLAG11=1; FLAG20=1; % Initialize flags from convergence if LDYN, VEL=VEL+(DYN.G*DELTA)*ACC; end % Update velocity after converged if(ITOL==1) % No bisection level DELTA=DELTA0; % Original time increment TARY(ITOL)=TIME+DELTA; % Time for current bisection level else % Recover previous bisection ITOL=ITOL-1; ISTEP=ISTEP-1; % Decrease the bisection level & load increment DELTA=TARY(ITOL)-TARY(ITOL+1); % Time increment b/w consecutive bisections TARY(ITOL+1)=0; % Empty converged bisection level end % End of bisection branch TIME0=TIME; % Save the current time % Stress update and print----------------------------------------------------- if ISTEP>=0 % Update stresses except for the 1st step UPDATE=true; LTAN=false; % Update stresses and history variables and store if NDOF==2 % 2D model if MID==0, ELAST2D(ETAN,NE,NDOF,XYZ,LE,DYNF); % Linear elastic material elseif MID>0, PLAST2D(MID,PROP,ETAN,NE,NDOF,XYZ,LE,DYNF); % Elastoplastic elseif MID<0, HYPER2D(PROP,NE,NDOF,XYZ,LE,DYNF); % Hyperleastic material end % End of material model branch else % 3D model if MID==0, ELAST3D(ETAN,NE,NDOF,XYZ,LE,DYNF); % Linear elastic material elseif MID>0, PLAST3D(MID,PROP,ETAN,NE,NDOF,XYZ,LE,DYNF); % Elastoplastic elseif MID<0, HYPER3D(PROP,NE,NDOF,XYZ,LE,DYNF); % Hyperleastic material end % End of material model branch end % End of DOF branch if LCNT, CNTELM(CNT,XYZ,CDISP,NUMNP,DYNF); end % Contact force & stiffness PROUT(NOUT,TIME,NUMNP,NE,NDOF,CNT); % Print results end % End of update stress branch CDISP=DISPTD; % Store converged displacement if LDYN, CVEL=VEL; CACC=ACC; end % Store converged velocity & acceleration TIME=TIME+DELTA; ISTEP=ISTEP+1; % Increase time and load increment % Check time and control bisection-------------------------------------------- while FLAG11 % Bisection loop FLAG11=0; % Turn off bisection flag if ((TIME-TIMEI)>1E-10) % Time passed the end time if ((TIMEI+DELTA-TIME)>1E-10) % One more increment at the end time DELTA=TIMEI+DELTA-TIME; % Time increment to the end DELTA0=DELTA; % Saved time increment TIME=TIMEI; % Current time is the end else % Move to next load step ILOAD=ILOAD+1; % Progress to next load step if(ILOAD>NLOAD), FLAG10=0; break; % Finished final load step, stop program else % Next load step TIME=TIME-DELTA; % Subtract old time increment TIMEF=PARAM.TIMS(1,ILOAD); TIMEI=PARAM.TIMS(2,ILOAD); % Start & end time CUR1=PARAM.TIMS(4,ILOAD); CUR2=PARAM.TIMS(5,ILOAD);% Start & end load fac. DELTA=PARAM.TIMS(3,ILOAD); DELTA0=DELTA; % Time increment & save TIME=TIME+DELTA; TDELTA=TIMEI-TIMEF; % Current time and time interval end % End of if(ILOAD>NLOAD) end % End of if ((TIMEI+DELTA-TIME)>1E-10) end % End of final time branch DISPDD=zeros(NEQ,1); % Zero incremental displacements if LDYN % Dynamic predictor DYNF=DYN.B*(DELTA)^2; % Stiffness multiplier DISPDD=DISPDD+DELTA*VEL+(0.5-DYN.B)*(DELTA^2)*ACC; % Inc. displ. predictor DISPTD=DISPTD+DISPDD; % Displacement predictor VEL=VEL+(1-DYN.G)*DELTA*ACC; % Velocity predictor ACC=zeros(NEQ,1); % Initialize acceleration end % End of dynamic problem branch FACTOR=CUR1+(TIME-TIMEF)/TDELTA*(CUR2-CUR1); % Current load factor SDISP=DELTA*SDISPT(:,3)/TDELTA*(CUR2-CUR1); % Increase fixed displacements if LCNT && ~isempty(CNT.CDISP) % Update contact rigid node coordinates CNT.XYZP=XYZ(NUMNP+1:end,:); % Coord. at previous time XYZ(CNT.CDISP(:,1),CNT.CDISP(:,2))=... % move rigid nodes CNT.XYZC(CNT.CDISP(:,1)-NUMNP,CNT.CDISP(:,2))+FACTOR*CNT.CDISP(:,3); % end % End of moving target rigid nodes % Start convergence iteration------------------------------------------------ ITER=0; % Iteration counter while FLAG20 % N-R Convergence interation loop FLAG20=0; ITER=ITER+1; % Reset iteration flag & increase iteration counter GKF=sparse(NEQ,NEQ); % Initialize global stiffness K FORCE=zeros(NEQ,1); % Initialize residual vector F UPDATE=false; LTAN=true; % Assemble K and F if NDOF==2 % 2D model if MID==0, ELAST2D(ETAN,NE,NDOF,XYZ,LE,DYNF); % Linear elastic material elseif MID>0, PLAST2D(MID,PROP,ETAN,NE,NDOF,XYZ,LE,DYNF); % Elastoplastic elseif MID<0, HYPER2D(PROP,NE,NDOF,XYZ,LE,DYNF); % Hyperleastic material end % End of material model branch if LDYN, MASS2D(NE,NDOF,XYZ,LE,DYN); end % Mass matrice & inertial term else % 3D model if MID==0, ELAST3D(ETAN,NE,NDOF,XYZ,LE,DYNF); % Linear elastic material elseif MID>0, PLAST3D(MID,PROP,ETAN,NE,NDOF,XYZ,LE,DYNF); % Elastoplasatic elseif MID<0, HYPER3D(PROP,NE,NDOF,XYZ,LE,DYNF); % Hyperleastic material end % End of material model branch if LDYN, MASS3D(NE,NDOF,XYZ,LE,DYN); end % Mass matrice & inertial term end % End of DOF branch if LCNT, CNTELM(CNT,XYZ,CDISP,NUMNP,DYNF); end % Contact force & stiffness % Applied loads and prescribed displacements BC----------------------------- if size(EXTFORCE,1)>0 % Increase external force using load factor LOC=NDOF*(EXTFORCE(:,1)-1)+EXTFORCE(:,2); % DOF of force BC FORCE(LOC)=FORCE(LOC) + FACTOR*EXTFORCE(:,3); % Add forces to residual end % End of external force branch NDISP=size(SDISPT,1); % # of prescribed displacements if NDISP~=0 % Fixed displacement BC FIXEDDOF=NDOF*(SDISPT(:,1)-1)+SDISPT(:,2); % DOF of displacement BC GKF(FIXEDDOF,:)=zeros(NDISP,NEQ); % Empty the disp. BC rows GKF(FIXEDDOF,FIXEDDOF)=sum(PROP(1:2))*eye(NDISP); % Replace diagonal term FORCE(FIXEDDOF)=0; % Replace residual term if ITER==1 % only at the 1st iteration FORCE(FIXEDDOF)=sum(PROP(1:2))*SDISP(:)/DYNF; % Replace displacement BC if LDYN % Acceleration BC FORCE(FIXEDDOF)=FORCE(FIXEDDOF)-sum(PROP(1:2))*DISPDD(FIXEDDOF)/DYNF; % end % End of dynamic analysis branch end % End of 1st iteration branch end % End of displacement BC branch % Check convergence and bisection------------------------------------------- if(ITER>1) % Check convergence after first iteration FIXEDDOF=NDOF*(SDISPT(:,1)-1)+SDISPT(:,2);% No residual check for fixed DOF FREEDOF=setdiff(1:NEQ,FIXEDDOF); % Free DOF RESN=max(abs(FORCE(FREEDOF))); % Max norm for residual if isempty(RESN), RESN = 0.0; end % When all DOFs are fixed OUTPUT(ITER-1,RESN,TIME,DELTA) % Print out iteration information if(RESNPARAM.ATOL)||(ITER>=PARAM.ITRA)) % Start bisection ITOL=ITOL+1; % Increase bisection level if(ITOL<=PARAM.NTOL) % Further bisection is allowed DELTA=0.5*DELTA; % Reduce time increment by half TIME=TIME0+DELTA; % Reduce the current time TARY(ITOL)=TIME; % Current bisection time stamp DISPTD=CDISP; % Restore previouse converged displacement if LDYN, VEL=CVEL; ACC=CACC; end % Restore velocity & acceleration if LCNT % % Compute current coordinates of slave nodes CNTCUR end % fprintf(1,'Not converged. Bisecting load increment %3d\n',ITOL); % else % Stop with too many bisections fprintf(1,'\t\t *** Max No. of bisection ***\n'); return; % end % End of no. of bisections branch FLAG11=1; FLAG20=1; break; % Stop iteration and start bisection end % End of bisection branch end % End of check convergence branch % Calculate displacement increment------------------------------------------ if(FLAG11==0) % Solve the system equation SOLN=GKF\FORCE; % Solve matrix equation DISPDD=DISPDD+DYNF*SOLN; % Update incremental displacement DISPTD=DISPTD+DYNF*SOLN; % Update total displacement if LDYN, ACC=ACC+SOLN; end % Dynamic analysis solves acc incr. FLAG20=1; % Continue convergence iteration else % Activate bisection FLAG20=0; % Start bisection. Out of convergence iteration end % End of solving the system equation branch if FLAG10, break; end % Converged. Out of convergence iteration end %FLAG20 Convergence iteration end %FLAG11 Bisection end %FLAG10 Load increment % Successful end of program---------------------------------------------------- fprintf(1,'\t\t *** Successful end of program ***\n'); % fprintf(NOUT,'\t\t *** Successful end of program ***\n'); % end % % function OUTPUT(ITER,RESN,TIME,DELTA) %****************************************************************************** % Print convergence iteration history %****************************************************************************** %% if ITER>1 % Print iteration history fprintf(1,'%27d %14.5e \n',ITER,full(RESN)); % Print residual else % First iteartion to print head fprintf(1,'\n \t Time Time step Iter \t Residual \n'); % Print head fprintf(1,'%10.5f %10.3e %5d %14.5e \n',TIME,DELTA,ITER,full(RESN)); % end % End of print history branch end % % function PROUT(NOUT,TIME,NUMNP,NE,NDOF,CNT) %****************************************************************************** % Print converged displacements and stresses % %****************************************************************************** %% % global SIGMA DISPTD FORCE % Global variables NITG=2^NDOF; % Integraion points (4 for 2D, 8 for 3D) fprintf(NOUT,'\r\n\r\nTIME = %11.3e\r\n\r\nNodal Displacements\r\n',TIME); % fprintf(NOUT,'\r\n Node U1 U2 U3'); % for I=1:NUMNP % Print nodal displacements II=NDOF*(I-1); % DOF location of Node I fprintf(NOUT,'\r\n%5d %11.3e %11.3e %11.3e',I,DISPTD(II+1:II+NDOF)); % end % End of nodal displacement loop fprintf(NOUT,'\r\n\r\nElement Stresses\r\n'); % Print element stress head fprintf(NOUT,'\r\n S11 S22 S33'); % fprintf(NOUT, ' S12 S23 S13'); % for I=1:NE % Print element stress fprintf(NOUT,'\r\nElement %5d',I); % Print element ID II=(I-1)*NITG; % Integration point counter fprintf(NOUT,'\r\n%11.3e %11.3e %11.3e %11.3e %11.3e %11.3e',... % SIGMA(1:6,II+1:II+NITG)); % Print element stress end % End of element stress loop if ~isempty(CNT) % If contact is present fprintf(NOUT,'\r\n\r\nContact Forces\r\n'); % Print contact head fprintf(NOUT,'\r\n Node F1 F2 F3'); % for I=1:size(CNT.CNODE,2) % Loop for all contact nodes II=NDOF*(CNT.CNODE(I)-1); % DOF location of contact node I fprintf(NOUT,'\r\n%5d %11.3e %11.3e %11.3e',CNT.CNODE(I),FORCE(II+1:II+NDOF)); end % End of contac nodes loop end % End of contact branch fprintf(NOUT,'\r\n\r\n'); % end % % function [PROP,ETAN]=PLSET(PROP,MID,NE,NDOF) %****************************************************************************** % Initialize history variables and elastic stiffness matrix % % XQ : 1-6 = Back stress alpha, 7 = Effective plastic strain % % SIGMA : Stress for rate-form plasticity % % : Left Cauchy-Green tensor XB for multiplicative plasticity % % ETAN : Elastic stiffness matrix % %****************************************************************************** %% % global SIGMA XQ % Global variables LAM=PROP(1)*PROP(2)/((1+PROP(2))*(1-2*PROP(2))); % Lame's constant Lambda MU=PROP(1)/(2*(1+PROP(2))); % Lame's constant Mu (shear modulus) PROP(1)=LAM; PROP(2)=MU; % Convert PROP(1:2) to Lame's constants N=(2^NDOF)*NE; % # of total integration points if MID>20 % Storage for multiplicative plasticity SIGMA=zeros(12,N); % stress & left deformation tensor XQ=zeros(4,N); % Principal back stress and eff. plastic strain SIGMA(7:9,:)=1; % Initial left deformation tensor is identity ETAN=[LAM+2*MU LAM LAM ; % Principal elasticity tensor LAM LAM+2*MU LAM ; % LAM LAM LAM+2*MU]; % else % Storage for classical plasticity SIGMA=zeros(6,N); % Stress XQ=zeros(7,N); % Back stress and eff. plastic strain ETAN=[LAM+2*MU LAM LAM 0 0 0; % Elasticity tensor LAM LAM+2*MU LAM 0 0 0; % LAM LAM LAM+2*MU 0 0 0; % 0 0 0 MU 0 0; % 0 0 0 0 MU 0; % 0 0 0 0 0 MU]; % end % End of material model branch end % % function VOLUME=CHECKELEM(XYZ,LE,NOUT) %****************************************************************************** % Check element connectivity and calculate volume % %****************************************************************************** %% % VOLUME=0; NDOF=size(XYZ,2); [NE,NENODE]=size(LE); % Problem parameters for I=1:NE % Loop for all elements ELXY=XYZ(LE(I,:),:); % Element nodal coordinates [~, ~, DET]=SHAPEL(zeros(1,NDOF),ELXY); % Determinant of Jacobian DVOL = (2^NDOF)*DET; % Element volume if DVOL