Please Help

%----------------------------------------------------------------------- % Impact of particle on a surface (with Fluid Drag) % % This program computes the motion of a particle subjected to % gravity forces and set in motion with initial velocity from % an initial position. The particle impacts a surface that is % defined in the function curve. The impact causes a loss of % energy through the coefficient of restitution. % % Numerical integration is by the generalized trapezoidal rule. % % Author: Keith D. Hjelmstad % Date: January 23, 2019 %----------------------------------------------------------------------- clear; clc; %. Input physical problem data gravity = 9.81; % acceleration of gravity e2 = [0; 1]; % unit vector in direction of gravity mass = 2.0; % mass of particle drag = 0.0; % Coefficient for drag force e = 0.9; % coefficient of restitution type = 1; % type of surface %. Initial conditions and analysis start and end times xo = [ 0.0; 85.0]; % initial position of particle vo = [ 10.0; 10.0]; % initial velocity of particle to = 0.0; % initial time (usually zero) tf = 25.0; % final time (end of analysis) %. Create a function for the surface such that when g > 0 then the % particle position is above the surface. Each row of the array % 'param' has the values of the parameters for the function % The variable 'bottom' will be a lower bound on vertical motion. switch type % param = [a, b, ...] case 1 %g = y - (a*x + b);% plane: Fcn='g(x,y)=y-(a*x+b)'; param = [ .25, 0]; bottom = -10; case 2 %g(x,y)=y-a*x^2-b*x-c Fcn='g(x,y)=y-(a*x^2+b*x+c)'; param=[.05,0,2]; bottom= -param(3); end %. Numerical analysis parameters. Put all parameters that are % associated with the numerical methods here. For implicit methods % you can put the parameters that describe and control Newton's method. dt = 0.01; % analysis time step beta = 0.5; % numerical integration parameter eta = (1-beta)*dt; % time integration parameter nsteps = ceil((tf-to)/dt); % number of time steps in analysis tf = nsteps*dt; % adjust final time to discrete times tol = 1.e-8; % Newton iteration tolerance itmax = 20; % maximum allowable Newton iterations %. Select figures to plot (1=yes, 0=no) dofig1 = 1; % 2D trajectory (static) dofig2 = 0; % Position components vs. time dofig4 = 1; % Energy vs. time dofig5 = 1; % Animation dots = 0; % =1 plot dots, =0 plot lines =2 both %. Set up a history array to save results ever once in a while % This part is simply to avoid plotting every point that we compute % It does not take that many points to give a nice representative % response curve. We need to compute more points to get accuracy in % the solution. nitems = 8; % # of items in 'history' array nreport = 2000; % # of steps to output to 'history' nsteps = ceil((tf-to)/dt); % # of analysis time steps nreport = min(nreport,nsteps); % don't report more than you compute noutput = ceil(nsteps/nreport); % set frequency of output history = zeros(nreport,nitems); % initialize storage for history %. Set time and output counters to get ready to start analysis t = to; out = 0; iout = 0; its = 0; %. Initialize position and velocity (put in 'xold' and 'vold' arrays) xold = xo; vold = vo; %. Compute initial acceleration from the equation of motion aold = -gravity*e2 - drag*norm(vold)*vold/mass; %. Compute motion by numerical integration--loop over time steps for i=1:nsteps %.... Write the values out to history, if appropriate % This segment writes a row to 'history' every 'noutput' steps % Note that for quantities that are not needed to advance the % numerical solution (e.g., speed and xcomp) we only need to % compute them if it is an output step. if (out==0) iout = iout + 1; %...... Compute energy PotEn = mass*gravity*dot(xold,e2); KinEn = 0.5*mass*dot(vold,vold); TotEn = PotEn + KinEn; %...... Store state in 'history' array history(iout,:) = [t,xold',vold',PotEn,KinEn,TotEn]; out = noutput; end out = out - 1; %.... Compute the values for the next time step t = t + dt; err = 1; its = 0; cn = vold + beta*dt*aold; anew = aold; %.... Newton iteration to compute the acceleration while (err > tol) && (its < itmax) its = its + 1; vnew = cn + eta*anew; speed = norm(vnew); mvec = vnew/speed; coeff = drag*speed*eta; g = mass*(anew + gravity*e2) + drag*speed*vnew; A = (mass + coeff)*eye(2) + coeff*(mvec*mvec'); anew = anew - A\g; err = norm(g); end % Newton 'while' loop %.... Compute new velocity and position by trapezoidal rule vnew = vold + dt*(beta*aold + (1-beta)*anew); xnew = xold + dt*(beta*vold + (1-beta)*vnew); %.... Check current position for impact. If contact is detected % then compute new velocity and note time of impact [gold,~] = CP2_plane(xold,param,type); [gnew,~] = CP2_plane(xnew,param,type); if (gnew < 0) && (gold >= 0) [gold,dgold] = CP2_plane(xold,param,type); dx = xnew - xold; dv = vnew - vold; da = anew - aold; frac = 0; err = 1; it = 0; while (err > tol) && (it < itmax) it = it + 1; x = xold + frac*dx; [g,dg] = CP2_plane(x,param,type); frac = frac - g/dot(dg,dx); err = abs(g); end %...... Update state to get to the exact point of contact t = t - dt + frac*dt; % time of contact xnew = xold + frac*dx; % position at contact vnew = vold + frac*dv; % velocity at contact anew = aold + frac*da; % acceleration at contact %...... Resolve contact by computing new velocity vector [g,dg] = CP2_plane(xnew,param,type); nvec = dg/norm(dg); % inward normal to surface vn = dot(vnew,nvec); % normal velocity v = -e*vn; % rebound velocity Eloss = 0.5*mass*(1-e^2)*vn^2; % energy loss from impact P = eye(2) - (nvec*nvec'); % projection matrix vnew = v*nvec + P*vnew; % velocity after contact %...... Write hit statistics and record the state in 'history' fprintf('\n t = %6.3f , Eloss = %9.3e ',t,Eloss) fprintf(' , frac = %6.3f , iter = %2i',frac,it) iout = iout + 1; history(iout,:) = [t,xnew',vnew',PotEn,KinEn,TotEn]; end if xnew(2) < bottom; break; end %.... Put current state (new) into old slot to get ready for new step % Note that this approach allows us to only store two states as % we carry out the calculation: 'old' and 'new' aold = anew; vold = vnew; xold = xnew; end % loop over time steps %. Echo input values to screen fprintf('\r\r%s\n', ' PHYSICAL PROPERTIES ') fprintf('%s%8.3f\n', ' Acceleration of gravity ',gravity) fprintf('%s%8.3f\n', ' Mass ',mass) fprintf('%s%8.3f\n', ' Drag coefficient ',drag) fprintf('%s%8.3f\n', ' Coefficient of restitution ',e) fprintf('%s%8.3f\n', ' Initial time (to) ',to) fprintf('%s%8.3f\n', ' Final time (tf) ',tf) fprintf('\r\r%s\n', ' INTEGRATION OF EQUATIONS ') fprintf('%s%8.3f\n', ' Time increment (dt) ',dt) fprintf('%s%8.3f\n', ' Numerical integration (beta) ',beta) fprintf('%s%8i\n', ' Number of time steps ',nsteps) fprintf('%s%8i\n', ' Number of steps to report ',nreport) fprintf('%s%8i\n', ' Output frequency (no. of steps)',noutput) fprintf('\r\r%s\n', ' INITIAL CONDITIONS ') fprintf('%s%8.3f\n', ' Initial position (x) ',xo(1)) fprintf('%s%8.3f\n', ' Initial position (y) ',xo(2)) fprintf('%s%8.3f\n', ' Initial velocity (x) ',vo(1)) fprintf('%s%8.3f\n', ' Initial velocity (y) ',vo(2)) %. GRAPHICAL OUTPUT SECTION % In this section we create plots to display the results that were % stored in the array history. Each row of history contains the values % of the response at a given time. The value "iout" refers to the last % row number in history. Thus, history(1:iout,1) is time, etc. ii = 1:iout; nCurvePts = 2000; minx = min(history(ii,2)); maxx = max(history(ii,2)); bufx = 0.1*(maxx-minx); maxx = maxx + bufx; minx = minx-bufx; minz = min(history(ii,3)); maxz = max(history(ii,3)); bufz = 0.1*(maxz-minz); maxz = maxz + bufz; minz = minz-bufz; xx = minx:(maxx-minx)/nCurvePts :maxx; xx = xx'; switch type case 1 yy = param(1).*xx + param(2); case 2 yy= param(1).*xx.^2 + param(2).*xx + param(3); end %. Plot 2D trajectory if (dofig1) figure(1); clf; grid on; axis equal; hold on; axis ([minx maxx minz maxz]); xlabel('x'); ylabel('y'); zlabel('z');title('2D Trajectory'); plot(xx,yy,'Color','k','LineWidth',2); if (dots==1) || (dots==2) p = plot(history(ii,2),history(ii,3),'o'); set(p,'MarkerFaceColor','b','Linewidth',1,... 'MarkerEdgeColor','b','MarkerSize',4); end if (dots==0) || (dots==2) p = plot(history(ii,2),history(ii,3),'-'); set(p,'Color','b','LineWidth',1); end plot(xo(1),xo(2),'o',... 'MarkerFaceColor','k','Linewidth',1,... 'MarkerEdgeColor','k','MarkerSize',8); xa = history(iout,2); ya = history(iout,3); plot(xa,ya,'o','MarkerFaceColor','r',... 'Linewidth',1,'MarkerEdgeColor','r','MarkerSize',6); end % dofig1 %. Plot x,y,z components of position vs. time if (dofig2) figure(2); clf; grid on; axis square; hold on; xlabel('t'); ylabel('x (blue), y (red)'); title('Position components vs. Time'); plot(history(ii,1),history(ii,2),'Color','blue','LineWidth',2); plot(history(ii,1),history(ii,3),'Color','red','LineWidth',2); end % dofig2 %. Plot energy if (dofig4==1) figure(4); clf; grid on; axis square; hold on; xlabel('t'); ylabel('U(blue), T(red) E(black)'); title('Energy vs. Time'); plot(history(ii,1),history(ii,6),'Color','b','LineWidth',2); plot(history(ii,1),history(ii,7),'Color','r','LineWidth',2); plot(history(ii,1),history(ii,8),'Color','k','LineWidth',2); end % dofig4 %. Plot 2D simulation as a movie if (dofig5==1) fig5=figure(5); kframes = 3; kforget = 10; numframes = floor(iout/kframes); M(numframes) = getframe(fig5); for k=1:numframes grid on; axis equal; axis ([minx maxx minz maxz]); xlabel('x'); ylabel('y'); zlabel('z');title('2D Trajectory'); hold on; %..... Plot the traces of the path of the two bodies kk = k - kforget; kk = (kk + abs(kk))/2; istart = 1 + kframes*kk; iend = kframes*k; p = plot(history(istart:iend,2),history(istart:iend,3)); set(p,'Color','blue','LineWidth',1); %..... Plot the particle at the end of the trace period xa = history(iend:iend,2); ya = history(iend:iend,3); p = plot(xa,ya,'o','MarkerFaceColor','k',... 'Linewidth',1,'MarkerEdgeColor','k','MarkerSize',8); plot(xx,yy,'Color','k','LineWidth',2); %..... Store frame for movie (set view point) M(k) = getframe(fig5); clf; end % loop over frames %... Play (and record) movie with play speed factor faster that actual clf; N = 1; PlaySpeed = 3; FPS = ceil(PlaySpeed*numframes/tf); movie(fig5,M,N,FPS) end % dofig5 block fprintf('\r%s\r',' Normal Termination of Program') %. End of program Contact_Plane %------------------------------------------------------------------- function [ g,grad_g ] = CP2_plane(z,params,type) % For a surface in 2D we input the z = [x;y] coordinates % and return the values of the function g and grad g where the % surface is defined as g(x,y)=0. The array params = [a,b,c,...] % holds the parameters needed to compute the function. The variable % 'type' defines the specific function. if type > 1; type = 1; end % default back to type 1 switch type case 1 % plane: g(x,y) = y - (a*x + b) a = params(1); b = params(2); x = z(1); y = z(2); g = y - (a*x + b); grad_g = [ -a; 1]; case 2 %plane: g(x,y)=y-a*x^2-b*x-c a = params(1); b = params(2); c= params(3); x = z(1); y = z(2); g=y-(a*x^2+b*x+c); grad_g = [ -2*a*x-b; 1]; end end