Blast and Impact Dynamics

Department of Civil & Structural Engineering

MatLab script for evaluating the cleared blast load acting on a finite target

Description

This MatLab code allows the user to generate a series of keyword files for use with LS-DYNA to predict the blast load and evaluate the effects of blast wave clearing on a finite-sized, deformable target. User documentation is provided.

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

%Hudson_Clearing.m
%==============================================================%
% Purpose written function to establish pressure-time functions on each
% node of a specified plate. Blast wave parameters are etablished using K&B
% empirical relationships and clearing functions given by Hudson method.
% User has option to write DYNA keyword files
%==============================================================%
% Please reference S. E. Rigby, A. Tyas, T. Bennett, J. A. Warren, and 
% S. Fay. Clearing effects on plates subjected to blast loads. Engineering 
% and Computational Mechanics, 166(3):140–148, 2013 when publishing
% academic work that uses this software.
 
% Note: Boundary conditions are set as clamped as default - this can be
% changed at the bottom of 'mesh.k'
%==============================================================%
clear
clc
% ALL UNITS IN N, s, m.
%================================================%
% 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 < 2)
    clc
    disp(['WARNING - Z = ' num2str(Z) ' m/kg^1/3'])
    disp('        - Scaled distance is too small for Hudson method to apply')
    disp('        - Scaled distance must be greater than 2 m/kg^1/3')
    break
elseif (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
 
% User defined mesh.
gridspace = input('Enter number of elements along length of plate (NB must be even): ');
    % Number of elements along lengths of plate - MUST BE EVEN!
    % Assumed to be the same for x and y.
%================================================%
% FIND BLAST WAVE PARAMETERS %
%================================================%
 
load('ConWepValues')
load('HudsonValues')
[psomax,prmax,ta,td,iso,ir] = 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 - ir) / (ir);
            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 - iso) / (iso);
            bso = bso + errso * bso;
        end
 
        clearvars err errso counter tol
%================================================%
% CLEARING CORRECTIONS %
%================================================%
 
%     ______________ x2,y2      ______________ sx2,sy2
%    |                |        |                |
%    |                |        |   ________tx2,ty2
%    |                |        |  |          |  |
%    |     TARGET     |        |  |  TARGET  |  |
%    |                |        |  |__________|  |
%    |                |        | tx1,ty1        |
%    |________________|        |________________|
%  x1,y1                   sx1,sy1
%            1.                       2.
 
targtype = 0;
disp('Select Configuration...')
disp('     ________________          ________________          ________________ ')
disp('    |                |        |                |        |                |')
disp('    |                |        |   __________   |        |                |')
disp('    |                |        |  |          |  |        |  __________    |')
disp('    |     TARGET     |        |  |  TARGET  |  |        | |          |   |')
disp('    |                |        |  |__________|  |        | |  TARGET  |   |')
disp('    |                |        |                |        | |__________|   |')
disp('    |________________|        |________________|        |________________|')
disp('            1.                        2.                        3.')
 
disp('  1. The target forms the entire reflecting surface')
disp('  2. The target is in the centre of a larger surface')
disp('  3. The target is an arbitrary distance from the edge of a larger surface')
while targtype~=1 && targtype ~=2 && targtype ~=3
  targtype = input('     and press Enter: ');
  if targtype ~=1 && targtype ~=2 && targtype ~=3
      disp('Not a valid entry. Please select either 1, 2 or 3')
  end
end
clc
 
% User defined target dimension
%================================================%
if targtype == 1;
    height = input('Enter target height (m): ');
    breadth = input('Enter target breadth (m): ');
 
    surfx1 = -breadth/2;
    surfx2 = breadth/2;
    if bursttype == 1; % surface burst
        surfy1 = 0;
        surfy2 = height;
    elseif bursttype == 2; % free air burst
        surfy1 = -height/2;
        surfy2 = height/2;
    end
 
    targx1 = surfx1;
    targx2 = surfx2;
    targy1 = surfy1;
    targy2 = surfy2;
 
%===================%    
elseif targtype == 2;
    surfheight = input('Enter surface height (m): ');
    surfbreadth = input('Enter surface breadth (m): ');
    height = 100;
    while height > surfheight 
        height = input('Enter target height (m): ');
        if height > surfheight;
            disp('Not a valid entry. Target must be smaller than surface')
        end
    end
    breadth = 100;
    while breadth > surfheight 
        breadth = input('Enter target breadth (m): ');
        if breadth > surfheight;
            disp('Not a valid entry. Target must be smaller than surface')
        end
    end
 
    surfx1 = -surfbreadth/2;
    surfx2 = surfbreadth/2;
    if bursttype == 1; % surface burst
        surfy1 = 0;
        surfy2 = surfheight;
    elseif bursttype == 2; % free air burst
        surfy1 = -surfheight/2;
        surfy2 = surfheight/2;
    end
 
    xoffset = (surfbreadth - breadth) /2;
    yoffset = (surfheight - height) /2;
 
    targx1 = surfx1 + xoffset;
    targx2 = targx1 + breadth;
    targy1 = surfy1 + yoffset;
    targy2 = targy1 + height;
 
%===================%    
elseif targtype == 3;
    surfheight = input('Enter surface height (m): ');
    surfbreadth = input('Enter surface breadth (m): ');
    height = 100;
    while height > surfheight 
        height = input('Enter target height (m): ');
        if height > surfheight;
            disp('Not a valid entry. Target must be smaller than surface')
        end
    end
    breadth = 100;
    while breadth > surfheight 
        breadth = input('Enter target breadth (m): ');
        if breadth > surfheight;
            disp('Not a valid entry. Target must be smaller than surface')
        end
    end
 
    disp('   ________________ ')
    disp('  |                |')
    disp('  |   __________   |')
    disp('  |  |          |  |')
    disp('  |  |  TARGET  |  |')
    disp('  |  |__________|  |  _ _')
    disp('  |                |   |')
    disp('  |________________|  _|_ y dist')
    disp('')
    disp('  |__|')
    disp('  |  | x dist')
 
    xoffset = input('Enter x distance to target edge (m): ');
    yoffset = input('Enter y distance to target bottom (m): ');
 
    if xoffset + breadth > surfbreadth;
        disp('WARNING - Target dimensions are too large')
        disp('        - Target exceeds reflecting surface')
        break
    elseif yoffset + height > surfheight;
        disp('WARNING - Target dimensions are too large')
        disp('        - Target exceeds reflecting surface')
        break
    end
 
    surfx1 = -surfbreadth/2;
    surfx2 = surfbreadth/2;
    if bursttype == 1; % surface burst
        surfy1 = 0;
        surfy2 = surfheight;
    elseif bursttype == 2; % free air burst
        surfy1 = -surfheight/2;
        surfy2 = surfheight/2;
    end
 
    targx1 = surfx1 + xoffset;
    targx2 = targx1 + breadth;
    targy1 = surfy1 + yoffset;
    targy2 = targy1 + height;
end
clc
 
% Clearing Lengths
%================================================%
% Calculate clearing lenghts using distance to free edge and positive phase
% duration.
 
if (R / max(height,breadth) < 3) % hypotenuse over 5% larger than charge range
    clc
    disp('WARNING - Target dimensions are too large')
    disp('        - Blast wave is not arriving planar')
    disp('        - Clearing predictions may not be valid')
end
 
xele = breadth/gridspace;
yele = height/gridspace;
 
xnodes = (targx1:xele:targx2);
ynodes = (targy2:-yele:targy1);
 
[Xall,Yall] = meshgrid(xnodes,ynodes);
[X,Y] = meshgrid(xnodes(2:gridspace),ynodes(2:gridspace));
 
xeta1 = (surfx2 - abs(X)) / (td * 340/1000);
yeta1 = (surfy2 - abs(Y)) / (td * 340/1000);
 
%================================================%
% PRESSURE-TIME CALCULATIONS %
%================================================%
% Inputs
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);
 
deltatime = zeros(size(xeta1)) - step/td;
po = zeros(size(xeta1)); ptemp = po;
p = zeros(size(xeta1,1),size(xeta1,2),numsteps);
    % p is an 3 dimensional array of an nxn matrix for each timestep. Numbers
    % indicate row,column and timestep respectively. i.e. p(10,10,:) is the
    % pressure-time history of node (10,10).
xclear = p; yclear = p;
i = 0;
 
presplot = 0;
disp('Animate Pressure-Time History?')
disp('  1. Yes')
disp('  2. No')
while presplot~=1 && presplot ~=2 
  presplot = input('     Select 1 or 2 and press Enter: ');
  if presplot ~=1 && presplot ~=2
      disp('Not a valid entry. Please select either 1 or 2')
  end
end
 
if presplot == 1
    az = -40;
    el = 40;
end
 
for t = ts;
    i = i+1;
    deltatime = deltatime + step/td;
 
    xdelta1 = deltatime - xeta1;
    ydelta1 = deltatime - yeta1;
 
    clearfunctx1 = interp2(delta,eta,z,xdelta1,xeta1);
    clearfuncty1 = interp2(delta,eta,z,ydelta1,yeta1);
 
    p(:,:,i) = po + pr(i) + ((clearfunctx1 + clearfuncty1) * psomax);
    xclear(:,:,i) = (clearfunctx1) * psomax;
    yclear(:,:,i) = (clearfuncty1) * psomax; 
 
    ptemp = p(:,:,i);
    if i > possteps
        p(:,:,i) = (ptemp<0) .* p(:,:,i);
    end
 
    if presplot == 1;           
        surf(X,Y,p(:,:,i))
        view(az, el);
        set(gca,'zlim',[-prmax,prmax])
        set(gca,'xlim',[targx1,targx2])
        set(gca,'ylim',[targy1,targy2])
        shading interp
        colormap(jet(256))
        xlabel('X')
        ylabel('Y')
        zlabel('Pressure (kPa)')
        pause(0.01)
    end
end
 
clearvars b bso delta deltatime eta i po xnodes ynodes z
 
tpzero = find(sum(sum(p))==0,1);
p(:,:,tpzero:end) = [];
ts(tpzero:end) = [];
xclear(:,:,tpzero:end) = [];
yclear(:,:,tpzero:end) = [];
clearvars tpzero
 
% LS-DYNA KEYWORD FILES
%================================================%
clc
keyword = 0;
disp('Write DYNA Keyword Files?')
disp('  1. Yes')
disp('  2. No')
while keyword~=1 && keyword ~=2 
  keyword = input('     Select 1 or 2 and press Enter: ');
  if keyword ~=1 && keyword ~=2
      disp('Not a valid entry. Please select either 1 or 2')
  end
end
 
if keyword == 1;
    % Mesh Generator - based on target geometry
        [elelist] = mesh(Xall,Yall);
        %==========================%
 
    xclear = xclear(1,:,:);
    yclear = yclear(:,1,:); 
    % Load Curve Generator - values taken from Hudson clearing calculations and
    % split into xdistance and ydistance loads. Pressure at a node is given as
    % superposition of inf + xdist + ydist. Uniform load stored in yload.k
        nodeloads(xclear,yclear,gridspace,ts,pr);
        %========================================%
 
    % Node Set Generator - All nodes along same vertical line have same x
    % clearing load and are in the same node set.
        nodesets(gridspace,xele,yele);
        %============================%
 
    % Segment Generator - Generates plate surface as a segment and applies
    % LOAD BLAST on segment
        segment(elelist,R,W,bursttype,ta);
        %================================%
 
    % Node load and node set ids named in the following convention:
    % 1. Boundary Nodes
    % 2. All Nodes (for uniform load)
    % 1001. First Node set for x distance
    % 2001. First Node set for y distance
end
clearvars Xall Yall keyword

Publication(s)

(2013). Clearing effects on plates subjected to blast loads. Proceedings of the Institution of Civil Engineers-Engineering and Computational Mechanics, 166 (3), (Full Text).