Blast and Impact Dynamics

Department of Civil & Structural Engineering

Blast.m - a simple tool for predicting blast pressure parameters

Description

This simple MatLab script will provide pressure, impulse, arrival time and duration predictions for a user defined explosive charge mass and stand-off. The parameter predictions for positive and negative phase are given by digitised data from UFC-3-340-02, Structures to Resist the Effects of Accidental Explosions. The negative phase is approximated with a cubic expression from Granstrom, Loading Characteristics of Air Blasts From Detonating Charges, Technical Report 100, Transactions of the Royal Institute of Technology, Stockholm.

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.

Code

%Blast.m
%==============================================================%
% Purpose written function to establish pressure-time and impulse-time
% curves for a given charge mass and stand-off using UFC-3-340-02
% semi-empirical parameters. Positive and negative phase can be calculated.
% Code written by Samuel E. Rigby
%==============================================================%
clear
clc
%================================================%
% INPUTS %
%================================================%
% User defined charge mass and stand off
W = input('Enter charge mass (kg TNT): ');
R = input('Enter stand-off (m): ');
 
% Calculate scaled distance
Wroot = W^(1/3);
Z = R/Wroot;
 
if (Z > 40)
    clc
    disp(['WARNING - Z = ' num2str(Z) 'm/kg^1/3'])
    disp('        - Scaled distance is too large for blast parameter calculations')
    disp('        - Scaled distance must be less than 40 m/kg^1/3')
    break
end
 
% Select burst type (surface or free air)
bursttype = 0;
disp('Select Configuration...')
disp('  1. Hemispherical Surface Burst')
disp('  2. Spherical Free-Air Burst')
while bursttype~=1 && bursttype ~=2 
  bursttype = input('     and press Enter: ');
  if bursttype ~=1 && bursttype ~=2
      disp('Not a valid entry. Please select either 1 or 2')
  end
end
 
%================================================%
% FIND BLAST WAVE PARAMETERS %
%================================================%
load('ConWepValues')
[psomax,prmax,ta,td,isomax,irmax] = PositiveParameters(R,W,bursttype,air,surface);
[prmin,psomin,negir,negiso] = NegativeParameters(R,W,bursttype);
 
negtd = (16/9) * (negir / prmin); % Integral of negative phase expression
 
        % Calculate Waveform Parameter
        %=========================%
        b = 10; bso = b; err = 100; errso = err; tol = 0.0001;
        counter = 0;
            % Initialise error values and tolerance.
        while (max(abs(err)) > tol) && (counter < 100);
            counter = counter + 1;
            ir_calc = prmax * td * (b - 1 + exp(-b)) / (b^2);
            err = (ir_calc - irmax) / (irmax);
            b = b + err * b;
        end
        counter = 0;
        while (max(abs(errso)) > tol) && (counter < 100);
            counter = counter + 1;
            iso_calc = psomax * td * (bso - 1 + exp(-bso)) / (bso^2);
            errso = (iso_calc - isomax) / (isomax);
            bso = bso + errso * bso;
        end
 
        clearvars err errso counter tol air surface ir_calc iso_calc
 
possteps = 100;
step = td/(possteps-1);
numsteps = round((td+negtd)/step);
ts = linspace(0,step*numsteps,numsteps);
pr = zeros(1,numsteps); pso = zeros(1,numsteps);
pr(ts<=td) = prmax * (1 - (ts(ts<=td) / td)) .* exp(-b*((ts(ts<=td) / td)));
pso(ts<=td) = psomax * (1 - (ts(ts<=td) / td)) .* exp(-bso*((ts(ts<=td) / td)));
 
pr(ts>td) = -prmin * 6.75 * ((ts(ts>td)-(td))/negtd) .* (((1 - (ts(ts>td)-(td))/negtd)).^2);
pso(ts>td) = -psomin * 6.75 * ((ts(ts>td)-(td))/negtd) .* (((1 - (ts(ts>td)-(td))/negtd)).^2);
 
ir = cumtrapz(ts,pr);
iso = cumtrapz(ts,pso);
 
%================================================%
% PLOT CONTROLS %
%================================================%
clc
%Select plot type
plottype = 0;
disp('Select Plot Type...')
disp('  1. Positive phase pressure-time')
disp('  2. Positive and negative phase pressure-time')
disp('  3. Positive phase pressure-time and impulse-time')
disp('  4. Positive and negative phase pressure-time and impulse-time')
disp('  5. No plot')
 
while plottype ~= 1:5
  plottype = input('     and press Enter: ');
    if plottype ~=1:5;
      disp('Not a valid entry. Please select [1] - [5]')
    end
end
clc
 
ts = ts+ta;
ts(3:end+2) = ts(1:end); ts(1) = 0; ts(2) = ta;
pr(3:end+2) = pr(1:end); pr(1:2) = 0;
pso(3:end+2) = pso(1:end); pso(1:2) = 0;
ir(3:end+2) = ir(1:end); ir(1:2) = 0;
iso(3:end+2) = iso(1:end); iso(1:2) = 0;
 
%================================================%
% OUTPUT %
%================================================%
disp('INPUT')
if bursttype ==1
disp('    Hemispherical Surface Burst')
elseif bursttype == 2
disp('    Spherical Free-Air Burst')
end
disp(['    Equivalent charge mass       = ' num2str(W) ' kg'])
disp(['    Stand-off                    = ' num2str(R) ' m'])
disp(['    Scaled distance              = ' num2str(Z) ' m/kg^1/3'])
disp('OUTPUT')
disp(['    Peak incident overpressure   = ' num2str(psomax) ' kPa'])
disp(['    Peak reflected overpressure  = ' num2str(prmax) ' kPa'])
disp(['    Time of arrival              = ' num2str(ta) ' ms'])
disp(['    Positive phase duration      = ' num2str(td) ' ms'])
disp(['    Incident impulse             = ' num2str(isomax) ' kPa.ms'])
disp(['    Reflected impulse            = ' num2str(irmax) ' kPa.ms'])
disp(['    Peak incident underpressure  = ' num2str(psomin) ' kPa'])
disp(['    Peak reflected underpressure = ' num2str(prmin) ' kPa'])
disp(['    Negative phase duration      = ' num2str(negtd) ' ms'])
 
 
if plottype == 1
    plot(ts(1,1:possteps+2),pr(1,1:possteps+2),'k','linewidth',1.5)     
    hold on
    plot(ts(1,1:possteps+2),pso(1,1:possteps+2),'k--','linewidth',1)
    ylabel('Pressure (kPa)','fontsize',14)
    xlabel('Time (ms)','fontsize',14)
    hleg1 = legend('Reflected Pressure','Incident Pressure');
    set(gca,'xlim', [floor(ta) ceil(ta+td)])
 
elseif plottype == 2
    plot(ts,pr,'k','linewidth',1.5)
    hold on
    plot(ts,pso,'k--','linewidth',1)
    ylabel('Pressure (kPa)','fontsize',14)
    xlabel('Time (ms)','fontsize',14)
    hleg1 = legend('Reflected Pressure','Incident Pressure');
    set(gca,'xlim', [floor(ta) ceil(ta+td+negtd)])
 
elseif plottype == 3
    plot(ts(1,1:possteps+2),pr(1,1:possteps+2),'k','linewidth',1.5)     
    hold on
    plot(ts(1,1:possteps+2),pso(1,1:possteps+2),'k--','linewidth',1)
    plot(ts(1,1:possteps+2),ir(1,1:possteps+2),'b','linewidth',1.5)
    plot(ts(1,1:possteps+2),iso(1,1:possteps+2),'b--','linewidth',1)  
    ylabel('Pressure (kPa), Impulse (kPa.ms)','fontsize',14)
    xlabel('Time (ms)','fontsize',14)
    hleg1 = legend('Reflected Pressure','Incident Pressure','Reflected Impulse','Incident Impulse');
    set(gca,'xlim', [floor(ta) ceil(ta+td)])
 
elseif plottype == 4
    plot(ts,pr,'k','linewidth',1.5)     
    hold on
    plot(ts,pso,'k--','linewidth',1)
    plot(ts,ir,'b','linewidth',1.5)
    plot(ts,iso,'b--','linewidth',1)  
    ylabel('Pressure (kPa), Impulse (kPa.ms)','fontsize',14)
    xlabel('Time (ms)','fontsize',14)
    hleg1 = legend('Reflected Pressure','Incident Pressure','Reflected Impulse','Incident Impulse');
    set(gca,'xlim', [floor(ta) ceil(ta+td+negtd)])
end
 
if plottype ~= 5
    set(get(gcf,'CurrentAxes'),'FontSize',12)
    set(legend,'FontSize',12)
end