% This script is a short driver to validate the direct time % integration functions. "direct" and "forcer" by comparing % the results to an analytic solution, by Biggs, for a % multi-degree of freedom system. % M*d2u(t)/dt2 + C*du(t)/dt + K*u(t) = F(t) % ------------------------------------------------------------ % The Bigg's example involves 3 masses moving wrt a fixed mass % at joint 1. The mass there is not important, but is included % to illustrate how a term is sometimes included but then % removed by a boundary condition. % % u1=0 -->u2 -->u3 -->u4 % M1---K1---M2---K2---M3---K3---M4 % Fixed ==>F2 ==>F3 <==F4 % % Forces decrease to zero linearly in 0.1 sec % System initially at rest at t=0 % % Typical data here are the spring assembly data, % Spring Joint1 Joint2 Stiffness (lb/in) % 1 1 2 6000 % 2 2 3 4000 % 3 3 4 2000 % % The mass assembly data, % Mass Joint Value (lb sec^2/in) % 1 1 4 (any) % 2 2 2 % 3 3 1 % 4 4 1 % % The force assembly data, % Force Joint Value (lb) % 2 2 3000 % 3 3 4000 % 4 4 -2000 % % The joints being restraned (to zero) % Restraint Joint % 1 1 % % Before restraint the matrices are: % Stiffness Mass % 6000 -6000 0 0 4 0 0 0 % -6000 10000 -4000 0 0 2 0 0 % 0 -4000 6000 -2000 0 0 1 0 % 0 0 -2000 2000 0 0 0 1 % % After the one restraint they have a new first row, column % Stiffness Mass % 1 0 0 0 1 0 0 0 % 0 10000 -4000 0 0 2 0 0 % 0 -4000 6000 -2000 0 0 1 0 % 0 0 -2000 2000 0 0 0 1 % % Alternatively, one could have run the reduced problem, with % no restraints given during integration using the matrices % Stiffness Mass % 10000 -4000 0 2 0 0 % -4000 6000 -2000 0 1 0 % 0 -2000 2000 0 0 1 % % The analytic time history solution was given by % Briggs, J. M., Introduction to Structural Dynamics", % McGraw Hill, 1964. % ----------------------------------------------------------------- function direct_forcer_validate % test of direct integration programs, full matrices % am * d2u(t)/dt2 + ac * du(t)/dt + ak * u(t) = f(t) % By WIlson algorithm (unconditionally stable Newmark method) % transient test problem Bigg's p. 123 n_d = 4 ; % Number of degrees of freedom n_c = 1 ; % Number of constraints k_print = 1 ; % Number of time steps between printouts n_steps = 81 ; % Total number of time steps dt = 0.0025 ; % Time step size omega = 1.25 ; % Wilson algorithm constant ak = zeros (n_d, n_d) ; % zero stiffess matrix am = zeros (n_d, n_d) ; % zero mass matrix ac = zeros (n_d, n_d) ; % zero damping matrix u = zeros (n_d, 1) ; % zero displacements f = zeros (n_d, 1) ; % zero applied forces vel = zeros (n_d, 1) ; % zero velocity acc = zeros (n_d, 1) ; % zero acceleration up = zeros (n_d, 1) ; % zero a working vector ibc = zeros (n_d, 1) ; % zero boundary condition flag % fix the top mass (binary flag) ibc (1) = 1 ; % set non-zero stiffness terms ak = [6000 -6000 0 0 ; ... -6000 10000 -4000 0 ; ... 0 -4000 6000 -2000 ; ... 0 0 -2000 2000 ] ; % set non-zero mass terms am = [4 0 0 0 ; ... 0 2 0 0 ; ... 0 0 1 0 ; ... 0 0 0 1] ; % % direct time integration, Wilson method direct (n_d, am, ac, ak, up, u, vel, acc, ... f, omega, dt, n_steps, k_print, n_c, ibc) end % direct_forcer_validate function direct (n_d, a, b, c, drp, r, dr, d2r, p, ... omega, dt, n_steps, i_print, n_c, ibc) % * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * % step by step integration of matrix equations: % a*d2r(t)/dt2 + b*dr(t)/dt + c*r(t) = p(t) % from initial conditions: r(0), dr(0)/dt, to linear system % s*r(t) = f(t) at each time. (by the Wilson method) % * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * % allocate a (n_d, n_d), b (n_d, n_d), ibc (n_c), ... % r (n_d), dr (n_d), d2r (n_d), p (n_d), ... % drp (n_d), c (n_d, n_d) zero = 0. ; % initial values of r and dr are passed thru arguments % a = mass (or capacitor) matrix % b = damping (or resistor) matrix % c = stiffness (or inductor) matrix % drp = value of dr at time = t + delt = tau % dt = time step size % f = linear system net forcing, stored over p % ibc = array containing the n_c dof's with zero bc % i_print = no. of integration steps between printing % n_c = no. dof with specified values of zero (if any) % n_d = total number of dof in system % n_steps = no. of integration steps, t_final = dt*n_steps % omega = 1.25 is suggested Wilson extrapolation amount % p = forcing function % r,dr,d2r = 0,1,2 order deriv. of r wrt t at time=t % s = linear system square matrix, stored over b % tau = extrapolated time step, omega*dt % ** initial calculations ** disp ('direct step by step integration') if i_print < 1 i_print = 1 ; % data check end n_steps = (n_steps/i_print)*i_print ; if n_steps <= 0 n_steps = i_print ; % data check end % if validation delt = (omega - 1)*dt ; % extrapolation increment tau = omega*dt ; % extrapolation amount if n_c < 1 disp ('warning: no constraints in direct \n') end % if warning % Set one node to display time history node_h = 3 ; % node for valudation h_id = fopen ('node_history.txt', 'w') ; % initialize arrays for loads< velocity< acceleration p = zeros(n_d,1); drp = zeros(n_d,1); d2r = zeros(n_d,1); % overwrite "b" with "s" matrix to save storage b = b + 2*a/tau + tau*c/3 ; % name change % ** apply boundary conditions ( to s ) ** if n_c > 0 % then apply essential boundary conditions for i = 1:n_c % BC value is zero here [b, p] = modify_full (n_d, ibc(i), zero, b, p) ; end % for on essential boundary conditions end % if boundary conditions needed % ** invert b (once) ** b = inv (b); % actually, factor if s is large % (approximate the initial value of d2r) % ** calculate solution at time t ** i_count = i_print - 2 ; for i_step = 0: n_steps % step zero is initial acceleration ===>> ===>> i_count = i_count + 1 ; t = dt*(i_step-1) ; t_plus = t + delt ; % extrapolate for the load if i_step <= 1 t_plus = zero ; % unless initial condition end % if if i_count == i_print fprintf(1,'\n step number = %5i time = %e \n', i_step, t); fprintf(1,'\n i r(i) dr/dt d2r/dt2 \n') end % if % forcer is to define the forcing function p(t) p = forcer (t_plus, n_d) ; % new force value % form modified forcing function, f, at t + delt for i = 1:n_d p(i) = p(i) + sum ( a(i,:)*( 2*dr(:)/tau + d2r(:)) ) ... - sum ( c(i,:)*(r(:) + 2*tau*dr(:)/3 + tau*tau*d2r(:)/6) ); end % for net forcing vector % ** apply boundary conditions ( to f ) ** if n_c > 0 p( ibc(1:n_c) ) = 0 ; % for zero value EBC end % if bc needed % solve for drp at time t+delt, s_inverse*f drp = b*p' ; % forward and back substitute if s is large % using data at t-dt and t+delt calculate values at t for i = 1:n_d drdt = (1-1/omega)*dr(i) + drp(i)/omega ; d2rdt2 = d2r(i)*(1-2/omega)+2*(drdt-dr(i))/omega/dt ; % approximate the initial value of d2r, if i_step <= 1 if i_step > 1 r(i) = r(i) + dt*(2*dr(i)+drdt)/3 + dt*dt*d2r(i)/6 ; dr(i) = drdt ; end % if it is a starter step d2r(i) = d2rdt2 ; end % for each degree of freedom % output results for time t if ( i_step > 1 ) fprintf(h_id,'%e %e %e %e \n', t, r(node_h), dr(node_h), d2r(node_h)); end % if skip initial acceleration at -del_t if i_count == i_print % then print this step for k = 1:n_d fprintf (1,'%5i %e %e %e \n', k, r(k), dr(k), d2r(k) ) end % for dof result listing at this time i_count = 0 ; end % if print increment end % for i_step over all times <<=== <<=== <<=== <<=== <<=== <<=== fclose (h_id) ; fprintf ('Saved Wilson results in file node_history.txt \n') end % function direct function [p] = forcer (t, n_d) % p = active force vector at time t % t = current time % n_d = number of terms in p % * * * * * * * * * * * * * * * * * * * * * * * * * * * * % define forcing function for step by step integration % NOTE: this requires application dependent functions % * * * * * * * * * * * * * * * * * * * * * * * * * * * * % transient test problem bigg's p. 123, % linear decreasing time ramp for forces f_of_t = 0.0 ; % default forcing value at t t_one = 0.1 ; % length of time ramp if ( t < t_one ) ; % force changing with time f_of_t = 1-t/t_one ; % current time multipler end % if within time ramp f2=3000; f3=4000; f4=-2000 ; % initial force values p(1) = 0 ; % new force at time t p(2) = f2*f_of_t ; % new force at time t p(3) = f3*f_of_t ; % new force at time t p(4) = f4*f_of_t ; % new force at time t end % function forcer function [s, c] = modify_full (n_d, n, value, s, c) %* * * * * * * * * * * * * * * * * * * * * * * * * * % apply an essential b.c. to full symmetric eqs % s*d = c, d(n) = value %* * * * * * * * * * * * * * * * * * * * * * * * * * % c = full column matrix % n = dof number of constrained parameter % n_d = total number of equations % s = full square matrix % value = given value of dof number n % note: reaction data are lost % s_nn = s(n, n) ; % get current diagonal term c = c - value*s (:, n) ; % move given column to right side s (:, n) = zeros (n_d, 1) ; % zero square matrix rows s (n, :) = zeros (1, n_d) ; % zero square matrix columns s (n, n) = s_nn ; % set identity in row n c (n) = value*s_nn ; % set identity in row n end % function modify_full