Calcinit.pl

#!/usr/bin/perl -w
# $Id: thesis.tex,v 1.21 1999/11/07 01:40:18 sjm Exp $
# calcinit.pl
#    Calculates initial conditions for use with convect.
use Getopt::Std;
# parameters
$visc0m = 4e-3;# set to reduce numerical fluctuations.
$press1m = $temp1m = 1;
$MAXX = 500; $MAXY = 100;

# given conditions
$leng = 75; $wid = 25;
$rconst = 101325*28.966/1.225014/288.16;
$mach = 10.4;
$gamma = 1.4;
$gammaHeAr = 2.5/1.5;
$MWair = 28.966;
$MWHeAr8515 = 4.628;

# condition 0. test gas
$press0 = 15000;
$temp0 = 288.16;
$rho0 = $press0*$MWair/$rconst/$temp0;

# Calculations from J.K.Wright "Shock Tubes"
$mu = ($gamma-1)/($gamma+1);
$y = $mach*$mach*(1+$mu)-$mu;
$uona0 = (1-$mu)*($mach-1/$mach);
$eta = ($mu+$y)/(1+$mu*$y);
$T1onT0 = $y*(1+$mu*$y)/($mu+$y);

# condition 1. post shock test gas
$press1 = $y*$press0;
$rho1 = $eta*$rho0;
$temp1 = $T1onT0*$temp0;
$a0 = sqrt($gamma*$press0/$rho0);
$u = $uona0*$a0;

# Air calculations
# from Bird, Steward & Lightfoot "Transport Phenomena"
# this is about 2.97
$ktoneair = $temp0/97;

$funcair = 1.048+($ktoneair-2.9)/0.1*(1.039-1.048);
$visc0air = 0.1*0.000026693
                  *sqrt($MWair*$temp0)/3.617/3.617/$funcair;
$reynoldsair = $rho0*$leng*$u/$visc0air;

# Condition 3. drive gas
# reservoir pressure and temperature from Skinner
# reservoir pressure and temperature from Skinner
$temp3 = 4598;
$a3 = sqrt($gammaHeAr * $temp3 * $rconst / $MWHeAr8515);
$press3 = $press0*$y*(1-($gammaHeAr-1)*($a0/$a3)
                 *($y-1)/sqrt(2*$gamma)
                 /sqrt(2*$gamma+($gamma+1)*($y-1)))
                 **(-2*$gammaHeAr/($gammaHeAr-1));

# Condition 2. post contact surface, driver gas
$press2 = $press1;# same as press1
# from Liepmann & Roshko
$temp2 = $temp3*($press2/$press3)**(($gammaHeAr-1)/$gammaHeAr);
$rho2 = $press2*$MWHeAr8515/$rconst/$temp2;

# HeAr8515 calculations
$ktoneHeAr = $temp2/(10.2*0.85+124*0.15);
$funcHeAr = 0.6335+($ktoneHeAr-60)/10*(0.6194-0.6335);
$visc0HeAr = 0.1*0.000026693*sqrt($MWHeAr8515*$temp2)
                      /2.576/2.576/$funcHeAr;
$reynoldsHeAr = $rho2*$leng*$u/$visc0HeAr;

# check command line options
getopts('rsn:');
if (defined $opt_r) {
    undef $opt_r;
    # just use the real conditions for output
    $DX = $wid;
    $DY = $leng;
    $rho1m = $rho1;
    $press1m = $press1;
    $temp1m = $temp1;
    $rho0m = $rho0;
    $um = $u;
    $temp0m = $temp0;
    $rho2m = $rho2;
    $temp2m = $temp2;
    $press2m = $press2;
} else {
    # calculate model conditions
    $rconstm = $MWair;
    $rho1m = $press1m*$MWair/$rconstm/$temp1m;
    $press0m = $press1m/$y;
    $temp0m = $temp1m/$T1onT0;
    $rho0m = $press0m*$MWair/$rconstm/$temp0m;
    $a0m = sqrt($gamma*$press0m/$rho0m);
    $um = $uona0*$a0m;
    $cs = $mach*$a0m;# shock speed

    # Condition 3. drive gas
    # reservoir pressure and temperature from Skinner
    $temp3m = 1;  # arbitrarily chosen
    $a3m = sqrt($gammaHeAr * $temp3m * $rconstm / $MWHeAr8515);
    $press3m = $press0m*$y
*(1-($gammaHeAr-1)*($a0m/$a3m)*($y-1)/
  sqrt(2*$gamma)/sqrt(2*$gamma+($gamma+1)*($y-1)))
    **(-2*$gammaHeAr/($gammaHeAr-1));
    
    # Condition 2. post contact surface, driver gas
    $press2m = $press1m;# same as press1
    # from Liepmann & Roshko
    $temp2m = $temp3m*($press2m/$press3m)
                       **(($gammaHeAr-1)/$gammaHeAr);
    $rho2m = $press2m*$MWHeAr8515/$rconstm/$temp2m;

    $DX = $reynoldsair * $visc0m / $rho0m / $um;
    $DY = $DX * $wid / $leng;
}

sub usage
{
    print <<EOF;
usage: calcinit.pl -s | [-r] -n #
   -s => print shock speed
   -n # => output conditions file for # species (# is 1 or 2).
   -r => print real conditions
EOF
}

# figure out what to output based on command line options
if (defined $opt_s)
{
    undef $opt_s;
    print "shock speed: $cs\n";
} elsif (defined $opt_n)
{
    if ($opt_n == 1)
    {
        # output conditions file for one species
print <<"EOF";
#inputs file generated by calcinit.pl
reynold1 $reynoldsHeAr
reynold2 $reynoldsair
MAXX $MAXX
MAXY $MAXY
DX $DX
DY $DY
RHO0 $rho1m
T0 $temp1m
PRESS0 $press1m
U0 $um
RHOINT $rho0m
TINT $temp0m
EOF
    } elsif ($opt_n ==2)
    {
        # output conditions file
print <<"EOF";
#inputs file generated by calcinit.pl
reynold1 $reynoldsHeAr
reynold2 $reynoldsair
MAXX $MAXX
MAXY $MAXY
DX $DX
DY $DY
RHO0 $rho2m
T0 $temp2m
U0 $um
PRESS0 $press2m
RHOINT $rho0m
TINT $temp0m
EOF
}
else
{
    &usage;
    exit(1);
}
} else
{
    &usage;
    exit(1);
}
EOF
;