% Danalysis: perform a dimensional analysis.
% Find non-dimensional variables via the basis set of the
% null space of the dimension matrix, D. Several examples,
% starting with finding the drag on a sphere moving at a known
% speed through a viscous fluid, are shown below. The lines
% having comments denoted by the symbol %%% are the
% lines that you will need to change to define a new problem.
% Enter the variables by calling Din; enter one dependent
% variable first, and then as many independent variables
% as appear in the problem. The dimensions are entered as a
% row vector of the exponents on the dimensions Mass, Length,
% time and Temperature, i.e, velocity is [0 1 -1 0]. The
% order of these dimensions is arbitrary, as is the choice
% of the fundamental dimensions themselves. You must be
% consistent within any one problem but otherwise the choice
% is yours.
%
% The raw output is a set of exponents, nsb, that are to be
% associated with each of the input variables. The specific
% order of the exponents, and thus the specific form of the
% non-dimensional variables, is determined by the arbitrary order
% in which you entered the independent variables. Only by
% accident would the form be optimal. Your part of this job is to
% manipulate the non-d variables into a useful form. For example,
% in the first problem of drag on a moving sphere, you might
% guess that inertial drag is the zero order solution, and hence
% you may want to divide Pi(1) by Pi(2) to show this. Note that
% Pi is a symbolic variable, and to do this you could
% enter >>Pi(1) = Pi(1)/Pi(2), and then >>pretty(simplify(Pi))
% (assuming that you have the symbolic toolbox). If you do not
% have the symbolic toolbox (you should get it!), you will still
% get the exponents on each variable returned through nsb.
%
% Written by Jim Price, July, 2001. jprice@whoi.edu, 508-289-2526.
% So far as I know as of late November, 2001, this technique is
% new, but I would be very surprised if it were new to the
% entire world; let me know if you have seen this before.
% Please also let me know if you find a problem on which this
% method fails or if you have suggestions or comments that
% would help improve the script.
%
% Public domain for all personal, educational purposes.
%
function Danalysis
% this has been set up as a function so that all the necessary
% subroutines can be included within one file, Din, Dclear,
% Dnullspace
global D Ds jp nsb Pi
% if you enter 'global Pi nsb' in the command window before
% running Danalysis, the exponent array nsb and the symbolic
% variable array Pi will be available in the command window
% for further use (or at least the last computed values will
% be available)
format compact
format rational
disp(' ')
disp(' *** drag on a sphere moving at a known speed ***') %%%
Dclear % always start with this
Din('Force', [1 1 -2] ) %%% drag (force) = [m*l/t*t]; the dependent variable
Din('rho', [1 -3 0] ) %%% density of the fluid = [m/l*l*l]
Din('radius', [0 1 0] ) %%% radius of the sphere = [l]
Din('mu', [1 -1 -1] ) %%% viscosity = [m/l*t]
Din('speed', [0 1 -1] ) %%% the known speed of the sphere = [l/t]
D; % print out D (or not)
[nsb Pi] = Dnullspace(Ds, D); % now perform the analysis on D
%
% You may prefer other forms of these Pis, as noted above.
%
disp(' hit any key to continue'), disp(' '), disp(' '), pause
% now try with viscosity omitted
disp(' *** drag on a sphere with viscosity omitted ***') %%%
Dclear % always start with this
Din('Force', [1 1 -2] ) %%% drag (force) = [m*l/t*t]; the dependent variable
Din('rho', [1 -3 0] ) %%% density of the fluid = [m/l*l*l]
Din('radius', [0 1 0] ) %%% radius of the sphere = [l]
% Din('mu', [1 -1 -1] ) %%% viscosity = [m/l*t]
Din('speed', [0 1 -1] ) %%% the known speed of the sphere = [l/t]
D; % print out D (or not)
[nsb Pi] = Dnullspace(Ds, D); % now perform the analysis on D
disp(' hit any key to continue'), disp(' '), disp(' '), pause
% now try with density of the fluid omitted
disp(' *** drag on a sphere with density omitted ***') %%%
Dclear % always start with this
Din('Force', [1 1 -2] ) %%% drag (force) = [m*l/t*t]; the dependent variable
% Din('rho', [1 -3 0] ) %%% density of the fluid = [m/l*l*l]
Din('radius', [0 1 0] ) %%% radius of the sphere = [l]
Din('mu', [1 -1 -1] ) %%% viscosity = [m/l*t]
Din('speed', [0 1 -1] ) %%% the known speed of the sphere = [l/t]
D; % print out D (or not)
[nsb Pi] = Dnullspace(Ds, D); % now perform the analysis on D
disp(' hit any key to continue'), disp(' '), disp(' '), pause
% now add g, as if near the sea surface, but w/ viscosity omitted
disp(' *** drag on a sphere, with g but not viscosity ***') %%%
Dclear % always start with this
Din('Force', [1 1 -2] ) %%% drag (force) = [m*l/t*t]; the dependent variable
Din('rho', [1 -3 0] ) %%% density of the fluid = [m/l*l*l]
Din('radius', [0 1 0] ) %%% radius of the sphere = [l]
Din('depth', [0 1 0] ) %%% depth = [l]
Din('g', [0 1 -2]) %%% accel of gravity
% Din('mu', [1 -1 -1] ) %%% viscosity = [m/l*t]
Din('speed', [0 1 -1] ) %%% the known speed of the sphere = [l/t]
D; % print out D (or not)
[nsb Pi] = Dnullspace(Ds, D); % now perform the analysis on D
disp(' hit any key to continue'), disp(' '), disp(' '), pause
% now add everything back in
disp(' *** drag on a moving sphere with all variables ***') %%%
Dclear % always start with this
Din('Force', [1 1 -2] ) %%% drag (force) = [m*l/t*t]; the dependent variable
Din('rho', [1 -3 0] ) %%% density of the fluid = [m/l*l*l]
Din('radius', [0 1 0] ) %%% radius of the sphere = [l]
Din('g', [0 1 -2]) %%% accel of gravity
Din('depth', [0 1 0] ) %%% depth = [l]
Din('mu', [1 -1 -1] ) %%% viscosity = [m/l*t]
Din('speed', [0 1 -1] ) %%% the known speed of the sphere = [l/t]
D; % print out D (or not)
[nsb Pi] = Dnullspace(Ds, D); % now perform the analysis on D
disp(' hit any key to continue'), disp(' '), disp(' '), pause
disp(' *** GI Taylor`s famous analysis of the first A bomb *** ')
% grim, but physically very interesting and historical
Dclear % reinitialize a counter and memory
Din('Energy', [1 2 -2]) %%% energy released (the unknown)
Din('radius', [0 1 0]) %%% observed radius of the blast wave
Din('time', [0 0 1]) %%% time elasped since the blast
Din('rho', [1 -3 0]) %%% density of the ambient fluid
% D; % print out D (or not)
[nsb Pi] = Dnullspace(Ds, D); % perform the analysis on D
%
% some data read off photos in Sedov (1959), Fig 59-61
% R is accurate to no better than about 10-20%
t = 1.e-3*[0.1 0.24 0.38 0.80 1.36 1.93 15.0 127.0]; % time, sec
R = 0.5*[24 36 50 68 80 90 220 340]; % radius of blast wave
E = 8.0e13; % (kinetic) energy released, Joules, found by
% an eyeball fit of the observed and non-d curves. This
% fitting would be plausible if the form (slope) of the
% non-d function were consistent with the data (it seems
% to be in this case). The best fit E corresponds roughly
% to 12 kilotons of TNT (perhaps more than you wanted to know?)
rho = 1.25; % air density
tmod = 1.e-4*[1:10:1000];
Rmod = (E/rho)^0.2*(tmod.^0.4);
% some graphics settings
set(0,'DefaultLineLineWidth',1.0)
set(0,'DefaultTextFontSize',14)
set(0,'DefaultAxesLineWidth',1.1)
set(0,'DefaultAxesFontSize',12)
figure(9)
clf reset
loglog(t, R, '*', tmod, Rmod)
legend('observed radius', 'fit of r(t), E = 8x10^{13} J')
xlabel('elapsed time, sec')
ylabel('radius of the blast wave, m')
title('radius(t) of a very intense blast wave')
disp(' hit any key to continue'), disp(' '), disp(' '), pause
disp (' ')
disp(' *** period of planetary orbits a la Kepler *** ') %%% title
Dclear % reinitialize a counter and memory
Din('T', [0 0 1]) %%% Period, the unknown
Din('R', [0 1 0]) %%% radius of the orbit
Din('m', [1 0 0]) %%% mass of the planet
Din('M', [1 0 0]) %%% mass of the sun
Din('G', [-1 3 -2]) %%% universal gravitational constant
% D; % print out D (or not)
[nsb Pi] = Dnullspace(Ds, D); % perform the nondimensional analysis of D
%
% Kepler's data
% planets are: Mercury, Venus, Earth, Mars, Jupiter, Saturn
Rp = [0.389 0.724 1.0 1.524 5.200 9.510]; % AU
Tp = [87.8 224.7 365.25 687.0 4332.6 10759.2]; % days
mp = [0.054 0.815 1.0 0.108 317.8 95.2]; % mass in units of Earth
Ms = 3.33e5*5.97e24; % mass of the sun
%
Rp = Rp*1.49e11; % convert AU to meters
Tp = Tp*8.64e4; % days to sec
mp = mp*5.97e24; % earth masses to kg
Gp = Rp.^3./(Ms*Tp.^2);
format short;
Gavg = mean(Gp);
G = 6.67e-11; % the universal gravitational constant
figure(33); clf reset
loglog(Rp, Tp, Rp, Tp, 'r*')
xlabel('radius, m'); ylabel('period, sec')
Tmod = 2*3.14159*(Rp.^1.5)/sqrt(G*Ms);
hold on
loglog(Rp, Tmod, 'b')
title('orbital period of the inner planets and Jupiter`s moons')
figure(34); clf reset
Tnd = Tp./((Rp.^1.5)/sqrt(Gavg*Ms));
loglog(mp/5.97e24, Tnd);
axis([0.01 1000 0.1 10.])
xlabel('mass planet/mass Earth'); ylabel('period, non-d')
% the moons of Jupiter; Io, Europa, Ganymede, Callisto
%
Rm = [5.57 8.87 14.15 24.90]*69.9e6; % m
Tm = [1.53 3.07 6.19 14.5]*1.e5; % sec
Mj = 317.8*5.97e24; % mass of Jupiter in kg
G = 6.67e-11;
figure(33); hold on
loglog(Rm, Tm, 'g*')
Tmod = 2*3.14159*(Rm.^1.5)/sqrt(G*Mj);
hold on
loglog(Rm, Tmod, 'b')
hold on
% loglog(Rm, Tm/sqrt(3.33e5/317.8), 'm')
axis([1.e8 1.e14 1.e3 1.e9])
axis('square')
legend('fit of p(r); G = 6.67 x10^{-11}',4)
function [nsb, Pi] = Dnullspace(Ds, D)
%
% Dnullspace: non-dimensional analysis of a dimension
% matrix D. Ds is a cell array of variable names. The
% output is a matrix nsb of exponents on the dimensional
% variables that will give a basis set of the null space
% of D. If symbolic toolbox is available, then the symbolic
% variable Pi is a printable version of the non-d
% variables. The resulting basis set is not unique, and is
% unlikely to conform to your preferred set of non-d
% variables. Jim Price, August, 2001.
disp(' ')
nrank = rank(D);
fprintf(' the rank of the dimension matrix is %g', nrank)
disp(' ')
[nv nd] = size(D);
numndv = nv - nrank;
fprintf(' the number of non-d variables is N = %g', numndv)
disp(' ')
% The next two lines are the key steps. First compute the null space
% basis of the dimensional matrix D by use of null, then put the
% result into reduced row echelon form by rref. The latter insures
% that the first dimensional variable (the dependent variable) will
% appear in one non-d variable only (the first one), and that
% it will have an exponent of 1. Beyond that, it is hard to say
% how nsb (the exponents of the basis set) will be arranged.
nsb = null(D.','r').'; % find a basis set of non-d variables.
% null seems to want the transpose of our D.
nsb = rref(nsb); % sort nsb so that pi(1) has an exponent
% =1 on the dependent variable.
% all that follows is an attempt to make a useful display of
% the result computed just above
[nr nc] = size(nsb);
% list the exponents on the input variables
disp(' ')
disp(' the dimensional variables and the exponents required to make')
disp(' one possible basis set of non-d variables are:')
disp(' ')
for j = 1:nc
fprintf(' %4s ', char(Ds(j)))
end
disp(' '), disp(' ')
for jr=1:nr
exponents = [sprintf('%9.4g ', nsb(jr,:))]
end
disp(' ')
% if the symbolic toolbox is available, then we can
% make a more easily interpreted display of the results
% (if not, you still get the exponents)
Pi = 0; % just in case you can't do the following
if exist('dsolve', 'file') == 2 % do you have symbolic toolbox?
clear Q QQ Pi
if exist('dsolve', 'file') == 2; eval('syms Q QQ Pi'); end
% store the nc variable names
for m=1:nc
QQ(m) = char(Ds(1,m));
end
% assemble the nr non-d variables
for j1=1:nr
Pi(j1) = prod(QQ.^nsb(j1,:));
end
disp(' one possible set of non-d variables (Pi(1) ... Pi(N)) is:')
pretty(Pi)
end % if on symbolic toolbox
disp(' ')
function Dclear
% Dclear: call this once at the start of a non-d
% analysis problem to reset some variables.
global D Ds jp nsb
D = [];
Ds = [];
nsb = [];
jp = 0;
function Din(x, y)
% Din: accept and store input to D and Ds. x is
% a string with the variable name, y is a row matrix
% of the exponents on each of the fundamental dimensions
% that appear in x. The length of y must be the same for
% each call to Din (within a given problem).
global D Ds jp nsb
jp = jp + 1;
if jp == 1
Ds = {x}; % store x in a cell array
else
Ds = [Ds {x}]; % add on to Ds with each call
end
ny = length(y); % the number of fundamental dimensions
D(jp, 1:ny) = y(1:ny); % store the exponents