Python script for ISA atmospheric calculation and velocity conversion

Here is small python script which calculates atmospheric properties and does velocity

conversion: Let me know if you need any explanation.

#!/usr/bin/python2.5

#####################################################################

## Script: Python script to calculate ISA properties and does velocity conversion ##

## ##

## Author: A. Khare ##

## Date: 18.04.2011 ##

## References: ESDU 77022 ##

## ##

#####################################################################

#Import modules

import sys, shutil, os, math

#Define global coefficients

global Gamma, R

Gamma = 1.4

R = 287.05287

#Define atmospheric coefficients

def Coeff(Altitude):

global a, b, c, D, E, F, G, i, j, L, m, n

if Altitude <= 36089.2:

a = 288.15

b = -1.9812 * 10**-3

c = 4.2927085

D = -29.514885 * 10**-6

E = 5.2558797

i = 0.24179285

j = -1.6624675 * 10**-6

L = 4.2558797

elif Altitude >= 36089.2 and Altitude <= 65616.8:

a = 216.65

b = 0

F = 2678.442

G = -48.063462 * 10**-6

m = 4.0012122 * 10**-3

n = -48.063462 * 10**-6

elif Altitude >= 65616.8 and Altitude <= 104987:

a = 196.65

b = 0.3048 * 10**-3

c = 0.79011202

D = 1.2246435 * 10**-6

E = -34.163218

i = 1.1616564

j = 1.8005232 * 10**-6

L = 35.163218

elif Altitude >= 104987 and Altitude <= 154199:

a = 139.05

b = 0.85344 * 10**-3

c = 0.47958369

D = 2.943516 * 10**-6

E = -12.201149

i = 1.3544167

j = 8.3129335 * 10**-6

L = -13.201149

elif Altitude >= 154199 and Altitude <= 164042:

a = 270.65

b = 0

F = 873.60472

G = -38.473855 * 10**-6

m = 1.04466 * 10**-3

n = -38.473855 * 10**-6

#

#Calculate temperature ratio

def Temperature_Ratio(Altitude):

if Altitude > 164042:

print "Altitude is out of valid range, max valid altitude is 164,042 ft"

sys.exit(0)

Coeff(Altitude)

temp_ratio = (a + b*Altitude)/288.15

return temp_ratio

#Calculate pressure ratio

def Pressure_Ratio(Altitude):

if Altitude > 164042:

print "Altitude is out of valid range, max valid altitude is 164,042 ft"

sys.exit(0)

Coeff(Altitude)

if b < 0 or b > 0:

Pressure_Ratio = (c + D*Altitude)**E

else:

Pressure_Ratio = F*math.exp(G*Altitude)

Pressure_Ratio = (Pressure_Ratio*4.4482)/(0.3048**2)/101325

return Pressure_Ratio

#Calculate density ratio

def Density_Ratio(Altitude):

if Altitude > 164042:

print "Altitude is out of valid range, max valid altitude is 164,042 ft"

sys.exit(0)

Density_Ratio = Pressure_Ratio(Altitude)/Temperature_Ratio(Altitude)

return Density_Ratio

#Calculate speed of sound

def Sonic_Velocity(Static_Temp_K):

if Static_Temp_K < 0:

print "Static temperature in Kelvin is out of valid range, Its below zero"

sys.exit(0)

Sonic_Velocity = (Gamma*R*Static_Temp_K)**0.5

return Sonic_Velocity

#Calculate dynamic viscosity ratio

def Dynamic_Viscosity_Ratio(Altitude):

Ts = Temperature_Ratio(Altitude)*288.15

Mu = (1.458 * (10**-6)) * (Ts**1.5) / (Ts + 110.4)

Dynamic_Viscosity_Ratio = Mu / (17.894 * (10**-6))

return Dynamic_Viscosity_Ratio

# Convert KCAS to Mach number

def CAS_to_Mach(Altitude, KCAS):

Ps = Pressure_Ratio(Altitude)

CAS_to_Mach = ((((1 + 0.2 * (KCAS / 661.48)**2)**3.5 - 1) * (1 / Ps) + 1)**(1 / 3.5) - 1)**0.5 * 2.2361

return CAS_to_Mach

# Convert KCAS to KTAS

def CAS_to_TAS(Altitude, KCAS):

Ps = Pressure_Ratio(Altitude)

Ts = Temperature_Ratio(Altitude)*288.15

CAS_to_TAS = ((((1+0.2*(KCAS/661.48)**2)**3.5-1)*(1/Ps)+1)**(1/3.5)-1)**0.5*2.2361*(Gamma*R*Ts)**0.5/0.51444

return CAS_to_TAS

# Convert KTAS to Mach number

def TAS_to_Mach(Altitude, KTAS):

Ts = Temperature_Ratio(Altitude)*288.15

TAS_to_Mach = KTAS / (((Gamma * R * Ts) ** 0.5) / 0.51444)

return TAS_to_Mach

# Convert KTAS to KCAS

def TAS_to_CAS(Altitude, KTAS):

Ps = Pressure_Ratio(Altitude)

Ts = Temperature_Ratio(Altitude)*288.15

Mach = KTAS / ((Gamma * R * Ts) ** 0.5 / 0.51444)

TAS_to_CAS = (((((Mach ** 2 * 0.2 + 1) ** 3.5 - 1) * (Ps) + 1) ** 0.2857 - 1) * 5) ** 0.5 * 661.48

return TAS_to_CAS

# Convert Mach number to KCAS

def Mach_to_CAS(Altitude, Mach):

Ps = Pressure_Ratio(Altitude)

Mach_to_CAS = (((((((1 + 0.2 * Mach ** 2) ** 3.5) - 1) * Ps + 1) ** (1 / 3.5)) - 1) * 5 * 661.48 ** 2) ** 0.5

return Mach_to_CAS

# Convert Mach number to KTAS

def Mach_to_TAS(Altitude, Mach):

Ps = Pressure_Ratio(Altitude)

Ts = Temperature_Ratio(Altitude)*288.15

Mach_to_TAS = Mach * (Gamma * R * Ts) ** 0.5 / 0.51444

return Mach_to_TAS

# Convert Mach number to KEAS

def Mach_to_KEAS(Altitude, Mach):

Ps = Pressure_Ratio(Altitude)

Ts = Temperature_Ratio(Altitude)*288.15

tas = Mach * (Gamma * R * Ts) ** 0.5 / 0.51444

Mach_to_KEAS = tas * Density_Ratio(Altitude) ** 0.5

return Mach_to_KEAS

if len(sys.argv) == 3:

func = sys.argv[1]

para1 = float(sys.argv[2])

if func == "Temperature_Ratio":

Out = Temperature_Ratio(para1)

elif func == "Pressure_Ratio":

Out = Pressure_Ratio(para1)

elif func == "Density_Ratio":

Out = Density_Ratio(para1)

elif func == "Sonic_Velocity":

Out = Sonic_Velocity(para1)

elif func == "Dynamic_Viscosity_Ratio":

Out = Dynamic_Viscosity_Ratio(para1)

else:

print "\033[1;31mType of argument for atmospheric property calculation is not valid. Valid arguments are:\033[1;m"

print "\033[1;31mTemperature_Ratio, Pressure_Ratio, Density_Ratio, Sonic_Velocity, Dynamic_Viscosity_Ratio\033[1;m"

sys.exit(0)

if func == "Sonic_Velocity":

print "%s at %g Kelvin is: %.6f" % (func, para1, Out)

else:

print "%s at %g feet is: %.6f" % (func, para1, Out)

elif len(sys.argv) == 4:

func = sys.argv[1]

para1 = float(sys.argv[2])

para2 = float(sys.argv[3])

if func == "CAS_to_Mach":

Out1 = CAS_to_Mach(para1, para2)

Out2 = "KCAS"

Out3 = "Mach"

elif func == "CAS_to_TAS":

Out1 = CAS_to_TAS(para1, para2)

Out2 = "KCAS"

Out3 = "KTAS"

elif func == "TAS_to_Mach":

Out1 = TAS_to_Mach(para1, para2)

Out2 = "KTAS"

Out3 = "Mach"

elif func == "TAS_to_CAS":

Out1 = TAS_to_CAS(para1, para2)

Out2 = "KTAS"

Out3 = "KCAS"

elif func == "Mach_to_CAS":

Out1 = Mach_to_CAS(para1, para2)

Out2 = "Mach"

Out3 = "KCAS"

elif func == "Mach_to_TAS":

Out1 = Mach_to_TAS(para1, para2)

Out2 = "Mach"

Out3 = "KTAS"

elif func == "Mach_to_EAS":

Out1 = Mach_to_KEAS(para1, para2)

Out2 = "Mach"

Out3 = "KEAS"

else:

print "\033[1;31mType of argument for aerodynamic velocity calculation is not valid. Valid arguments are:\033[1;m"

print "\033[1;31mCAS_to_Mach, CAS_to_TAS, TAS_to_Mach, TAS_to_CAS, Mach_to_CAS, Mach_to_TAS, Mach_to_EAS\033[1;m"

sys.exit(0)

print "%.6f %s at %g feet is: %.6f %s" % (para2, Out2, para1, Out1, Out3)

else:

print "\033[1;31mWrong number of arguments to the script\033[1;m"

print "Use for Atmospheric Properties:"

print "cfd_atm Property_Type Altitude[feets]"

print "Example:"

print "cfd_atm Pressure_Ratio 30000\n"

print "Use for Aerodynamic Velocities:"

print "cfd_atm Converson_Type Altitude[feets] Input[Velocity/Mach]"

print "Example:"

print "cfd_atm CAS_to_TAS 30000 250"

sys.exit(0)

#End of script

#####################################################################

Script can be downloaded from here:

atm_velo_calculator 

Python script to process Cl values in list of directories

Following is the small python script which post process all the directories inside a given directories and store values of Alpha vs CL in a file “Alpha_CL.dat”.

This is helpful when your are running some polar calculations in CFD++.

Use: script common_name_of_directories

Example: If you have AoA polar calculations with directories name as AoA00, AoA02, AoA04 and so on run the following script as:

script AoA

#!/usr/bin/python2.5

#####################################################################

## Script: Python script to post process CL values for all the ##

## cases in given directory ##

## ##

## Author: A. Khare ##

## Date: 03.03.2011 ##

## ##

#####################################################################

#Import modules

import sys, os, glob, commands

#Check for argument. If argument is not given print the use and exit

if len(sys.argv) != 2:

print "Use: get_cl prefix_for_cases"

exit(1)

wdir = os.getcwd()

#Open file for writing CL values

f = open("Alpha_CL.dat", 'w')

#Reading the prefix and finding out all the directories

list = str(sys.argv[1] + "*")

dir_list = glob.glob(list)

#Loop over all the directories and process CL

for i in dir_list:

j = wdir + '/' + i

os.chdir(j)

alpha = float(commands.getoutput("cat infout1f.inp | grep alpha | head -1 | awk \'{print $2}\'"))

cl = float(commands.getoutput("infout1f 1 | tail -21 | head -1 | awk \'{print $9}\'"))

f.write('%.2f\t%.3f\n' % (alpha, cl))

os.chdir(wdir)

f.close()

#End of script

########################################################################

Script to generate cut sections for a CFD++ solution

This script will generate cut sections in CFD++/Tecplot format from a CFD++ solution file.

Techplot cut sections can be visualized in ParaView.

How to use this script:

Save the below script with a name say cfd_cuts and execute it in following way:

./cfd_cuts Min Max Number_of_cuts Axis

For example to make 10 cuts in X direction starting from 1.2 to 3 issue following command:

./cfd_cuts 1.2 3 10 x

#!/usr/bin/tclsh8.4

#User inputs

#Minimum value of axis

 set Min [format %0.3f [lindex $argv 0]]

#Maximum value of axis

 set Max [format %0.3f [lindex $argv 1]]

#Number of cuts in between Min and max value

 set Number [format %0.3f [lindex $argv 2]]

#Axis in which we need these cuts

 set Direction [lindex $argv 3]

#Calculating the intervals for cuts

 set Interval [expr (($Max-$Min)/($Number))]

#Opening files for writing CFD++ commands for cuts

 set fp [open "section_cuts.sh" w] 

 set fpt [open "section_cuts_tec.sh" w] 

#Loops for number of cuts

 set j $Min

 for {set i 0 } {$i <= $Number} {incr i} {

 puts $fp "npfcutpl $Direction= [expr $j] cellsin.bin  pltosout.bin mpf $Direction=[expr $j]m"

 puts $fpt "mpf3dtecp cut_$Direction=[expr $j]m.mpf3d 

cut_$Direction=[expr $j]m_tec"

 set j [expr $j+$Interval]

 }

#Closing files

 close $fp

 close $fpt

#Executing cuts commands for generating cuts in CFD++ format

 exec /usr/bin/sh section_cuts.sh

 puts "\033\[01;32mIf you want to convert CFD++ sections files into Tecplot files press y otherwise press n\n\033\[0m"

 set input [gets stdin]

 scan $input "%s" answer

 if {$answer == "y"} {

#Executing cuts commands for generating cuts in Tecplot format, these files can be visualize in ParaView

 exec /usr/bin/sh section_cuts_tec.sh

  } elseif {$answer == "n"} {

    exit 0

    } else {

    exit 0

    }

Let me know if you need any explanation.

Any comments are welcome. 

This script can be downloaded from here:

cfd_cuts 

 

Script to create Mach number and Alpha polar for CFD++

I have written a TCL script which creates Mach number and Alpha polar for CFD++.
To use this script you have to specify your Mach number range, your AoA range and number of processor on which you want to run your jobs.
You also need to set up your first case using GUI for starting mach number and starting AoA .
In your polar directory you should have following files:


cellsin.bin   
exbcsin.bin
nodesin.bin 
mcfd.bc
mcfd.inp
infout1f.inp
cfd_sub (Your job launching file based on your scheduler SGE , PBS etc.)
polar (Polar script copy and paste from below)
inputs.dat (Input for polar script)


inputs.dat file should be in following format:


MACH_LIST 0.12 0.5 1 1.17

ALPHA_LIST 0.2 2.3 3.9 6

NPROCS 128


polar file is as follows:

#!/usr/bin/tclsh8.4

  foreach fl {mcfd.inp inputs.dat infout1f.inp cellsin.bin exbcsin.bin nodesin.bin mcfd.bc cfd_sub} {

  if {[file exists $fl] == 0} {

     puts "$fl file does not exist. Can not continue. Aborting ...."

     exit 0

        }

  }

 set fui [open "inputs.dat" r]

 set file_data [read $fui]

 set data [split $file_data "\n"]

 set mkey [lindex $data 0]

 set ml [llength $mkey]

 set akey [lindex $data 1]

 set al [llength $akey]

 set nprocs [lindex $data 2]

 set NPROCS [lindex $nprocs 1]

 close $fui

 set Old_mach_dir [glob -nocomplain Mach_*]

 foreach filename $Old_mach_dir {

 if {[file exists $filename]} {

    puts "Previous $filename directory found. Deleting it.\n"

    file delete -force $filename

    puts "Deleted previous $filename directory\n"

    } else {

       puts "No Previous $filename file found\n"

           }

 }

  set Pres [exec cat mcfd.inp | grep -A2 primitive | cut -d " " -f2 | head -2 | tail -1]

  set Temp [exec cat mcfd.inp | grep -A2 primitive | cut -d " " -f3 | head -2 | tail -1]

  set Uin [exec cat mcfd.inp | grep -A2 primitive | cut -d " " -f4 | head -2 | tail -1]

  set Vin [exec cat mcfd.inp | grep -A2 primitive | cut -d " " -f5 | head -2 | tail -1]

  set Win [exec cat mcfd.inp | grep -A2 primitive | cut -d " " -f6 | head -2 | tail -1]

  set Win [exec cat mcfd.inp | grep -A2 primitive | cut -d " " -f6 | head -2 | tail -1]

  set U [expr {sqrt($Uin*$Uin+$Vin*$Vin+$Win*$Win)}]

  set Alti [exec cat mcfd.inp | grep aero_altid | awk {{print $2}}]

  exec sed -e {s/^.*aero_pres.*/aero_pres PRESS/} -e {s/^.*aero_temp.*/aero_temp TEMP/} -e {s/^.*aero_u.*/aero_u VELOX/} -e {s/^.*aero_v.*/aero_v VELOY/} -e {s/^.*aero_w.*/aero_w VELOZ/} -e {s/^.*aero_ma.*/aero_ma MACH/} -e {s/^.*aero_alpha.*/aero_alpha ALPHA/} -e s/$Pres/PRESS/g -e s/$Temp/TEMP/g -e s/$Uin\ $Vin\ $Win/VELOX\ VELOY\ VELOZ/g < mcfd.inp > mcfd.data

  exec sed -e {s/^.*alpha.*/alpha ALPHA/} -e {s/^.*uref.*/uref UREF/} < infout1f.inp > infout1f.data

  set Sound [exec /opt/METACOMP/CFD++10.1/mlib/mcfd.10.1/exec/atmosprp 0.0 kilometers | tr -s ' ' | grep Sound | tr -s ' ' | cut -d " " -f11]

  set WDIR [exec pwd]

   set jf [open "$WDIR/jobs.sh" w]

 puts "\n"

 puts "=========================================================="

 puts "You have selected following inputs for your simulations =="

 puts "=========================================================="

 puts "\n"

 puts [exec cat inputs.dat]

 puts "\n"

 puts "\033\[01;32mType y/n and press enter to continue/quit\n\033\[0m"

 set input [gets stdin]

 scan $input "%s" answer

 if {$answer == "y"} {

    puts "\nContinuing ...\n"

  } elseif {$answer == "n"} {

    puts "\nQuitting ...\n"

    exit 0

    } else {

    puts "\nQuitting ...\n"

    exit 0

    }

 if {[file exists mcpusin.bin.$NPROCS]} {

 puts "\033\[01;32mmcpusin.bin.$NPROCS already exist. If you would like to use it press y otherwise press n\n\033\[0m"

 set finput [gets stdin]

 scan $finput "%s" fanswer

 if {$fanswer == "y"} {

 puts "\nPrevious mcpusin.bin.$NPROCS will be used\n"

 } elseif {$fanswer == "n"} {

 puts "\nCrating $NPROCS partitions\n"

 exec /opt/METACOMP/CFD++10.1/mlib/mcfd.10.1/exec/tometis pmetis $NPROCS

 }

 } else {

 puts "\nHICrating $NPROCS partitions\n"

 exec /opt/METACOMP/CFD++10.1/mlib/mcfd.10.1/exec/tometis pmetis $NPROCS

 }

 puts "Creating directory structure\n"

 for {set j 1} {$j<$ml} {incr j} {

  set Mach [lindex $mkey $j]

  set DIR Mach_$Mach

  file mkdir $DIR

  cd $DIR

 for {set i 1} {$i<$al} {incr i} {

   set Alpha [lindex $akey $i]

   set Starting_Alpha [lindex $akey 1]

   set DIR Alpha_$Alpha

   file mkdir $DIR

   cd $DIR

   set CDIR [exec pwd] 

   set U [expr {$Mach*$Sound}]

   set u [expr {$U*cos(($Alpha*3.14159)/180)}]

   set v 0.0

   set w [expr {$U*sin(($Alpha*3.14159)/180)}]

   file copy ../../mcfd.data .

   file copy ../../infout1f.data .

   file copy ../../cfd_sub .

   file link -symbolic cellsin.bin ../../cellsin.bin

   file link -symbolic nodesin.bin ../../nodesin.bin

   file link -symbolic exbcsin.bin ../../exbcsin.bin

   file link -symbolic mcfd.bc ../../mcfd.bc

   file link -symbolic mcpusin.bin.$NPROCS ../../mcpusin.bin.$NPROCS

   puts "Creating mcfd.inp file for Mach=$Mach and Alpha=$Alpha\n"

   exec  sed -e s/PRESS/$Pres/ -e s/TEMP/$Temp/ -e s/VELOX/$u/ -e s/VELOY/$v/ -e s/VELOZ/$w/ -e s/MACH/$Mach/ -e s/ALPHA/$Alpha/ < mcfd.data > mcfd.inp

   exec  sed -e s/ALPHA/$Alpha/ -e s/UREF/$U/ < infout1f.data > infout1f.inp

   file delete mcfd.data

   file delete infout1f.data

   if {$Alpha != $Starting_Alpha} {

   exec sed -e s/mc_filecopy\ cdepsout.bin\ cdepsin.bin/mc_filecopy\ \.\.\\/Alpha_[lindex $akey [expr $i-1]]\\/cdepsout.bin\ cdepsin.bin/ -e s/istart\ 0/istart\ 1/ < mcfd.inp > mcfd.inp.new

   file copy -force mcfd.inp.new mcfd.inp

   file delete mcfd.inp.new

   }

   puts $jf "cd $CDIR"

   puts $jf "./cfd_sub"

   puts $jf "sleep 30s"

   cd ..

  }

  cd ../

 }

 close $jf

 #Deleting unused files

 file delete mcfd_metis.graph

 file delete mcfd.data

 file delete infout1f.data

 puts "Directory structure has been created and all required files are placed accordingly\n"

 puts "A job file called jobs.sh has been created in your working directory\n"

 puts "\033\[01;32mIf you want to fire your jobs press y otherwise press n\n\033\[0m"

  set input_j [gets stdin]

 scan $input_j "%s" answer

 if {$answer == "y"} {

 puts "\nContinuing ...\n"

 exec sh jobs.sh

  } elseif {$answer == "n"} {

    exit 0

    } else {

    exit 0

    }

Let me know if you need any explanation.

This is a preliminary script hence there is lot of scope to improve it.
So please post your comments.

To check the script with out installing CFD++ do following things:
Create following files in a directory using touch command:
touch cellsin.bin nodesin.bin exbcsin.bin mcfd.bc mcpusin.bin.128 cfd_sub
Download following files:

polar
inputs.dat
infout1f.inp
mcfd.inp 

Now to generate the Mach and Alpha polar run polar script by ./polar