% Example 6.3 Central Difference Method A=1e-4;E=1e11;rho=4000;L=1/3;Fmax=0.001*A*E;tmax=3e-3; %Problem parameters Dt=21e-6; Nmax=round(tmax/Dt)+1; % Time step and no. of steps t = linspace(0, tmax, Nmax); % Discrete times K=(A*E/L)*[2 -1 0;-1 2 -1;0 -1 1]; % Stiffness matrix M=(rho*A*L/6)*[4 1 0;1 4 1;0 1 2]; % Consistent mass matrix F=[zeros(2,Nmax); % Force vector linspace(0,Fmax,round(Nmax/2)) Fmax*ones(1,Nmax-round(Nmax/2))]; Q=zeros(3,Nmax); % Nodal displacement vector Minv=(Dt^2)*inv(M); % Inverse of mass matrix Q(:,2)=Q(:,1); % Due to zero initial conditions for n=2:Nmax-1 % Loop for time steps Q(:,n+1) = 2*Q(:,n) - Q(:,n-1) + Minv*(F(:,n)-K*Q(:,n)); % Eq. (6.35) end % End of loop Qs=K\F; %Quasi-static solution plot (t,Q(3,:),t,Qs(3,:)); % Plotting tip displacements % % Example 6.4 Newmark method A=1e-4;E=1e11;rho=4000;L=1/3;Fmax=0.001*A*E;tmax=3e-3; %Problem parameters gamma=0.5;beta=0.25;Dt=42e-6;Nmax=round(tmax/Dt)+1; % Newmark parameters & time t = linspace(0, tmax, Nmax); % Discrete times K=(A*E/L)*[2 -1 0;-1 2 -1;0 -1 1]; % Stiffness matrix M=(rho*A*L/6)*[4 1 0;1 4 1;0 1 2]; % Consistent mass matrix F=[zeros(2,Nmax); % Force vector linspace(0,Fmax,round(Nmax/2)) Fmax*ones(1,Nmax-round(Nmax/2))]; Minv=inv(M+(beta*Dt^2)*K); % Inverse of effective mass matrix Q=zeros(3,Nmax);V=zeros(3,Nmax);A=zeros(3,Nmax); % Displ., vel., acc. vectors for n=1:Nmax-1 % Loop for time steps Vpr = V(:,n)+(1- gamma)*Dt*A(:,n); % Velocity predictor Qpr = Q(:,n)+Dt*V(:,n)+(0.5-beta)*Dt^2*A(:,n); % Displacement predictor A(:,n+1)=Minv*(F(:,n+1)-K*Qpr); % Calculate acceleration V(:,n+1)=Vpr + gamma*Dt*A(:,n+1); % Velocity corrector Q(:,n+1)=Qpr + (beta*Dt^2)*A(:,n+1); % Displacement corrector end % End of time step loop Qs=K\F; % Static solution plot (t,Q(3,:),t,Qs(3,:)); % Plotting tip displacements % % Example 6.5 3D uniaxil bar XYZ=[0 0 0; 0 0.01 0; 0 0.01 0.01; 0 0 0.01; % Nodal coordinates 1/3 0 0;1/3 0.01 0;1/3 0.01 0.01;1/3 0 0.01; % 2/3 0 0;2/3 0.01 0;2/3 0.01 0.01;2/3 0 0.01; % 1 0 0; 1 0.01 0; 1 0.01 0.01; 1 0 0.01]; % LE=[1 2 3 4 5 6 7 8; 5 6 7 8 9 10 11 12;9 10 11 12 13 14 15 16]; % Elements EXTFORCE=[13 1 2500;14 1 2500;15 1 2500;16 1 2500]; % External forces SDISPT=[1 1 0;1 2 0;1 3 0;2 1 0;2 3 0;3 1 0;4 1 0;4 2 0]; % Fixed displacements tmax=3E-3; Dt=42E-6; % Maximum time and time increment TIMS=[0 tmax/2 Dt 0 1; tmax/2 tmax Dt 1 1]'; % Load increments [t0 T Dt f0 fL] MID=0; PROP=[1E11 0.0]; % Elastic material [E NU] PARAM=struct('ITRA',50,'ATOL',1.0E6,'NTOL',6,'TOL',1E-5,'TIMS',TIMS); DYN=struct('B',0.25,'G',0.5,'RHO',4000,'IV',[]); % Dynamic parameters clc; NOUT=fopen('output.txt','w'); % Open output file NLFEA(PARAM,NOUT,MID,PROP,EXTFORCE,SDISPT,XYZ,LE,[],DYN); % Calling NLFEA fclose(NOUT); % Close output file DATA=EXTDIS('output.txt',15,1); % Extract x-displacement of Node 15 plot(DATA(1,:),DATA(2,:)); % Plotting displacements plot(DATA(1,:),DATA(2,:),t,Q(3,:),t,Qs(3,:)); % Plotting displacements % % Problem 6.10 Central Difference Method A=1e-4;E=1e11;rho=4000;L=1/3;Fmax=0.001*A*E;tmax=30e-3; %Problem parameters Dt=21e-6; Nmax=round(tmax/Dt)+1; % Time step and no. of steps t = linspace(0, tmax, Nmax); % Discrete times K=(A*E/L)*[2 -1 0;-1 2 -1;0 -1 1]; % Stiffness matrix M=(rho*A*L/6)*[4 1 0;1 4 1;0 1 2]; % Consistent mass matrix F=[zeros(2,Nmax); % Force vector linspace(0,Fmax,round(Nmax/2)) Fmax*ones(1,Nmax-round(Nmax/2))]; Q=zeros(3,Nmax); % Nodal displacement vector Minv=(Dt^2)*inv(M); % Inverse of mass matrix Q(:,2)=Q(:,1); % Due to zero initial conditions for n=2:Nmax-1 % Loop for time steps Q(:,n+1) = 2*Q(:,n) - Q(:,n-1) + Minv*(F(:,n)-K*Q(:,n)); % Eq. (6.35) end % End of loop Qs=K\F; %Quasi-static solution plot (t,Q(3,:),t,Qs(3,:)); % Plotting tip displacements % % Problem 6.10 Central Difference Method A=1e-4;E=1e11;rho=4000;L=1/3;Fmax=0.001*A*E;tmax=30e-3; %Problem parameters Dt=21e-6; Nmax=round(tmax/Dt)+1; % Time step and no. of steps t = linspace(0, tmax, Nmax); % Discrete times K=(A*E/L)*[2 -1 0;-1 2 -1;0 -1 1]; % Stiffness matrix M=(rho*A*L/6)*[4 1 0;1 4 1;0 1 2]; % Consistent mass matrix F=[zeros(2,Nmax); % Force vector linspace(0,Fmax,round(Nmax/2)) Fmax*ones(1,Nmax-round(Nmax/2))]; Q=zeros(3,Nmax); % Nodal displacement vector Minv=(Dt^2)*inv(M); % Inverse of mass matrix Q(:,2)=Q(:,1); % Due to zero initial conditions for n=2:Nmax-1 % Loop for time steps Q(:,n+1) = 2*Q(:,n) - Q(:,n-1) + Minv*(F(:,n)-K*Q(:,n)); % Eq. (6.35) end % End of loop Qs=K\F; %Quasi-static solution plot (t,Q(3,:),t,Qs(3,:)); % Plotting tip displacements % % Problem 6.14 Bar impact problem (Newmark method) Area=1;E=100;rho=0.01;L=1;tmax=1;gamma=0.5; beta=0; Dt=sqrt(rho/E); Nmax=round(tmax/Dt)+1; t = linspace(0, tmax, Nmax); %Global matrices K=2*eye(10);M=2*eye(10);K(1,1)=1;M(1,1)=1; for n=1:9, K(n,n+1)=-1; K(n+1,n)=-1; end K=(Area*E/L)*K; M=(rho*Area*L/2)*M; Minv=inv(M + (beta*Dt^2)*K); Q=zeros(10,Nmax);V=zeros(10,Nmax);A=zeros(10,Nmax);sigma=zeros(3,Nmax); %Initial conditions V(:,1)=ones(10,1); A(:,1)=-inv(M)*K*Q(:,1);sigma(:,1)=0; %Time integration for n=1:Nmax-1 %Predictor Vpr = V(:,n)+(1-gamma)*Dt*A(:,n); Qpr = Q(:,n)+Dt*V(:,n)+(0.5-beta)*Dt^2*A(:,n); A(:,n+1)=-Minv*K*Qpr; %Corrector V(:,n+1)=Vpr + gamma*Dt*A(:,n+1); Q(:,n+1)=Qpr + (beta*Dt^2)*A(:,n+1); sigma(1,n+1)=E*(Q(2,n+1)-Q(1,n+1))/L; sigma(2,n+1)=E*(Q(6,n+1)-Q(5,n+1))/L; sigma(3,n+1)=E*(-Q(10,n+1))/L; end figure; plot(t,Q(1,:),t,Q(6,:)); axis([0 1 -0.1 0.1]); figure; plot(t,sigma(1,:),t,sigma(2,:),t,sigma(3,:));