Computational Mechanics and Design Group

Department of Civil & Structural Engineering

MatLab code for 1D dynamic wave propagation with gradient elastic dispersion

Description

The code below is a MatLab script showing an animation of the effects of dispersion in a heterogeneous bar (red line); the wavefront in an homogenous bar (blue line) is shown for comparison.

This code is distributed freely for academic use, with the caveat that publications (including Undergraduate, MSc & PhD dissertations, conference and journal papers) using this work fully acknowledge its use.

Research area(s)

Code

% dyn1d.m
% 1D dynamic wave propagation with gradient elastic dispersion
% 2 paramter model:
% Askes, Bennett & Aifantis
% A new formulation and C^{0}-implementation of dynamically consistent gradient elasticity
% IJNME, 2007, 72:111-126
%
% Terry Bennett
% Dept. of Civil and Structural Engineering
% University of Sheffield
 
clear all
clc
format long e
 
disp('1D Gradient elastic dispersion in laminates')
disp(' ' )
disp('Int. J. Fract. (2007) 148:185-193' )
disp('Terry Bennett, Inna M. Gitman & Harm Askes')
disp(' ')
disp('Dept. of Civil and Structural Engineering')
disp('University of Sheffield')
disp(' ')
disp(' ')
 
% input data
bigl = 300.0;               % bar length
nbars = 300;                % no of bars
nnodes = nbars + 1;         % no of nodes
 
l = bigl / nbars;          % length of 1 bar
E = 1.0;                 % Young's modulus
A = 1;                     % cross sectional area
rho = 1.0;               % mass density
 
% laminate material
% Int J Frac, 2007, 148:185-193
E1=1e6; rho1=1.8;
alpha = 0.5;
E2=(E*(1-alpha)*E1)/(E1-E*alpha);
rho2=2*rho-rho1;
 
disp(['volume fraction: alpha ',num2str(alpha)])
disp(['material 1: E = ',num2str(E1),' rho = ',num2str(rho1)])
disp(['material 2: E = ',num2str(E2),' rho = ',num2str(rho2)])
 
top=alpha*(1-alpha)*(E1*rho1-E2*rho2);
bot=(1-alpha)*rho*E1 + alpha*rho*E2;
gamma = (1/12)*(top/bot)^2;
 
unitcell = 1;
 
beta = 1/12 * unitcell;
l_m = (gamma + beta)^0.5 * 1.0;
l_s = beta^0.5 * 1.0;
 
% time steps
ce = sqrt(E/rho);
deltaT = 1.0;              % time step size
nstep  = 280/deltaT;       % no of time steps
 
% external force (nnodes-1 as nth node is fixed)
R = zeros(2*(nnodes-1),1);
R(1,1) = 1.0;  
 
Rfac = zeros(1,nstep);
Rfac(1,1) = 1.0/deltaT;                 % load factor unit jump
 
% stiffness matrix
K = assemK(A,E,l,nnodes,nbars);
 
% mass matrix
M = assemM(A,rho,l,nnodes,nbars, l_s,l_m);
 
% initial conditions
dis_old = zeros(2*(nnodes-1),1);     % displacement
vel_old = zeros(2*(nnodes-1),1);     % velocity
acc_old = zeros(2*(nnodes-1),1);     % acceleration
 
% Newmark integration constants (Bathe notation)
NMdelta = 0.5;
NMalpha = 0.25;
a0 = 1 / (NMalpha * deltaT * deltaT);
a1 = NMdelta / (NMalpha * deltaT);
a2 = 1 / (NMalpha * deltaT);
a3 = ( 1/(2 * NMalpha) ) - 1;
a4 = ( NMdelta / NMalpha ) - 1;
a5 = (deltaT / 2) * ( (NMdelta/NMalpha) - 2 );
a6 = deltaT * (1 - NMdelta);
a7 = NMdelta * deltaT;
 
% effective stiffness matrix
Kstar = K + a0*M;
Kinv = inv(Kstar);
 
% ------------------------------------------------------------------
% ---- time step loop ----------------------------------------------
'time step loop'
 
dis = zeros(2*nnodes,nstep+1); dis(1:2*(nnodes-1),1) = dis_old;
vel = zeros(2*nnodes,nstep+1); vel(1:2*(nnodes-1),1) = vel_old;
acc = zeros(2*nnodes,nstep+1); acc(1:2*(nnodes-1),1) = acc_old;
 
for istep = 1 : nstep
 
  disp(['time step: ',num2str(istep),' time: ',num2str(istep*deltaT)])
 
  % efective loads at t + deltaT    
  Rstar = R*Rfac(1,istep) ...
        + M * (a0*dis_old + a2*vel_old + a3*acc_old);
        %+ C * (a1*dis_old + a4*vel_old + a5*acc_old)
 
  % -----------------------
  % ---- solve ------------
  %dis_new = Kstar \ Rstar;
  dis_new = Kinv*Rstar;
  % ---- solve ------------
  % -----------------------
 
  % new accelerations
  acc_new = a0*(dis_new - dis_old) - a2*vel_old - a3*acc_old;
  % new velocities
  vel_new = vel_old + a6*acc_old + a7*acc_new;
 
  %update
  dis_old = dis_new; vel_old = vel_new; acc_old = acc_new;
 
  % store in array
  dis(1:2*(nnodes-1),istep+1) = dis_new;
  vel(1:2*(nnodes-1),istep+1) = vel_new;
  acc(1:2*(nnodes-1),istep+1) = acc_new;
 
end
 
% ---- time step loop ----------------------------------------------
% ------------------------------------------------------------------
% time snap shot plot
anim8

Publication(s)