In [1]:
%pylab inline
from paraview.simple import *
paraview.simple._DisableFirstRenderCameraReset()
import pylab as pl
import math
import numpy as np
from scipy.optimize import curve_fit
remote_data = True
In [6]:
#from paraview.simple import *
#paraview.simple._DisableFirstRenderCameraReset()
import math
import numpy as np
from scipy.optimize import curve_fit
datum = -150.0 #nose at -150.0m
alpha = 0.0 # degrees
beta = 0.0 # degrees
reference_area = 1.0e6 # mm^2
reference_length = 1.0e3 # mm
reference_point = [0.0,0.0,0.0]
pressure = 87510.53 # Pa
density = pressure/(287.0*288.15) # kg/m^3
gamma = 1.4
speed_of_sound = math.sqrt(gamma*pressure/density) # m/s
#bl_position = 0.97
#face_area = 0.075823
def rotate_vector(vec,alpha_degree,beta_degree):
"""
Rotate vector by alpha and beta based on ESDU definition
"""
alpha = math.radians(alpha_degree)
beta = math.radians(beta_degree)
rot = [0.0,0.0,0.0]
rot[0] = math.cos(alpha)*math.cos(beta)*vec[0] + math.sin(beta)*vec[1] + math.sin(alpha)*math.cos(beta)*vec[2]
rot[1] = -math.cos(alpha)*math.sin(beta)*vec[0] + math.cos(beta)*vec[1] - math.sin(alpha)*math.sin(beta)*vec[2]
rot[2] = -math.sin(alpha)* vec[0] + math.cos(alpha)* vec[2]
return rot
def sum_and_zone_filter_array(input,array_name,ignore_zone,filter=None):
sum = [0.0,0.0,0.0]
p = input.GetCellData().GetArray(array_name)
z = input.GetCellData().GetArray("zone")
numCells = input.GetNumberOfCells()
for x in range(numCells):
zone = z.GetValue(x)
if zone not in ignore_zone:
v = p.GetTuple(x)
if filter == None or filter.test(input,x):
#print 'Zone: %i'%(zone)
for i in range(0,3):
sum[i] += v[i]
return sum
def sum_and_zone_filter(input,array_name,ignore_zone,filter=None):
sum = [0.0,0.0,0.0]
if input.IsA("vtkMultiBlockDataSet"):
iter = input.NewIterator()
iter.UnRegister(None)
iter.InitTraversal()
while not iter.IsDoneWithTraversal():
cur_input = iter.GetCurrentDataObject()
v = sum_and_zone_filter_array(cur_input,array_name,ignore_zone,filter)
for i in range(0,3):
sum[i] += v[i]
iter.GoToNextItem();
else:
sum = sum_and_zone_filter_array(input,array_name,ignore_zone,filter)
return sum
class GeomFilterLT:
def __init__(self,val,idx):
#
self.val = val
self.idx = idx
def test(self,input,x):
centre = input.GetCellData().GetArray("centre").GetTuple(x)
if centre[self.idx] < self.val:
return True
else:
return False
class GeomFilterGT:
def __init__(self,val,idx):
#
self.val = val
self.idx = idx
def test(self,input,x):
centre = input.GetCellData().GetArray("centre").GetTuple(x)
if centre[self.idx] >= self.val:
return True
else:
return False
def calc_force(file_root,ignore_zone,half_model=False,filter=None):
wall = PVDReader( FileName=file_root+'_wall.pvd' )
#CellDatatoPointData1 = CellDatatoPointData(Input=wall)
#CellDatatoPointData1.PassCellData = 1
#wall = CellCenters(Input=wall)
#sum = MinMax(Input=wall)
#sum.Operation = "SUM"
#sum.UpdatePipeline()
wall.UpdatePipeline()
sum_client = servermanager.Fetch(wall)
pforce = sum_and_zone_filter(sum_client,"pressureforce",ignore_zone,filter)
#fforce = sum_client.GetCellData().GetArray("frictionforce").GetTuple(0)
#yplus = wall_slice_client.GetCellData().GetArray("yplus").GetValue(0)
pforce = rotate_vector(pforce,alpha,beta)
#fforce = rotate_vector(fforce,alpha,beta)
if half_model:
for i in range(0,3):
pforce[i] *= 2.0
return pforce
def calc_moment(file_root,ignore_zone,half_model=False):
wall = PVDReader( FileName=file_root+'_wall.pvd' )
#CellDatatoPointData1 = CellDatatoPointData(Input=wall)
#CellDatatoPointData1.PassCellData = 1
#sum = MinMax(Input=wall)
#sum.Operation = "SUM"
#sum.UpdatePipeline()
sum_client = servermanager.Fetch(wall)
pmoment = sum_and_zone_filter(sum_client,"pressuremoment",ignore_zone)
#fforce = sum_client.GetCellData().GetArray("frictionforce").GetTuple(0)
#yplus = wall_slice_client.GetCellData().GetArray("yplus").GetValue(0)
pmoment = rotate_vector(pmoment,alpha,beta)
#fforce = rotate_vector(fforce,alpha,beta)
if half_model:
# This is only valid for X-Z plane reflection
pmoment[0] += -pmoment[0]
pmoment[1] += pmoment[1]
pmoment[2] += -pmoment[2]
return list(pmoment)
def calc_lift_centre_of_action(force,moment,ref_point):
# longitudinal centre xs0 at zs0
# spanwise centre ys0 at zs0
# residual Mz moment (Mx=My=0) mzs0
xs0 = ref_point[0] - moment[1]/force[2]
ys0 = ref_point[1] + moment[0]/force[2]
zs0 = ref_point[2]
mzs0 = moment[2] - force[1]*(xs0-ref_point[0]) + force[0]*(ys0-ref_point[1])
return (xs0,ys0,zs0),mzs0
def get_case(case,num_procs):
return 'C13_2013-04-04-half-'+case+'_P'+num_procs+'_OUTPUT/C13_2013-04-04-half-'+case
#Connect('localhost')
if remote_data:
#print 'Waiting for remote paraview server connection'
ReverseConnect('11111')
num_procs = '24'
case_list = [('m0p5',0.5,0.0,0.0),
('m0p8',0.8,0.0,0.0),
('m1p1',1.1,0.0,0.0),
('m1p3',1.5,0.0,0.0)]
#case_list = [('m0p5',0.5,0.0,0.0)]
force = []
force_fwd = []
force_aft = []
moment = []
ignore_zones = [34,35,101]
#ignore_zones = []
for case in case_list:
case_name = case[0]
mach = case[1]
q = (0.5*gamma*pressure*mach*mach)
alpha = case[2]
beta = case[3]
#f_aft = calc_force(get_case(case_name,num_procs),set(ignore_zones),True,GeomFilterGT((16.0+datum)*1000.0,0))
#for i in range(0,3):
# f_aft[i] /= reference_area
#print f_aft
#break
#calc_test(get_case(case_name,num_procs),set(ignore_zones),True)
f = calc_force(get_case(case_name,num_procs),set(ignore_zones),True)
for i in range(0,3):
f[i] /= reference_area
force.append(f)
m = calc_moment(get_case(case_name,num_procs),set(ignore_zones),True)
for i in range(0,3):
m[i] /= (reference_length*reference_area)
moment.append(m)
p, mz = calc_lift_centre_of_action(f,m,reference_point)
print 'Centre of lift: %f m' % (p[0]-datum)
# Calc load fwd of aero centre
f_fwd = calc_force(get_case(case_name,num_procs),set(ignore_zones),True,GeomFilterLT(p[0]*1000.0,0))
for i in range(0,3):
f_fwd[i] /= reference_area
force_fwd.append(f_fwd)
print f_fwd
# Calc load aft of aero centre
f_aft = calc_force(get_case(case_name,num_procs),set(ignore_zones),True,GeomFilterGT(p[0]*1000.0,0))
for i in range(0,3):
f_aft[i] /= reference_area
force_aft.append(f_aft)
print f_aft
if remote_data:
#print 'Disconnecting from remote paraview server connection'
Disconnect()
x = [0.0,0.5,0.8,1.1,1.3,1.4]
y = []
y.append(force[0][0])
for f in force:
y.append(f[0])
y.append(force[-1][0])
from scipy.interpolate import interp1d
drag_curve = interp1d(x, y, kind='linear')
xnew = np.linspace(0.0, 1.4, 20)
xlabel('Mach')
ylabel('D/q', multialignment='center')
title('Drag Profile')
pl.plot(x,y,xnew,drag_curve(xnew))
In [28]:
from scipy.interpolate import interp1d
# Calculate total thrust for jet and rocket
thrust_factor = 1.0
#jet_thrust = 10000.0 # 10KN thrust for ??
jet_thrust = 4.44822162*20000.0 # 20000lb thrust
jet_end = 75.0 # s
#rocket_thrust = 12000.0 # 12KN for 20 secs
rocket_thrust = 4.44822162*26000.0
rocket_start = 20.0
rocket_end = rocket_start + 20.0
def get_thrust(time_array):
thrust_array = []
for t in time_array:
thrust = jet_thrust*thrust_factor
if t >= rocket_start and t <= rocket_end:
thrust += rocket_thrust
thrust_array.append(thrust)
return thrust_array
# Time sample at 1s intervals
time_sample = 0.1
time_array = np.linspace(0.0, 100.0, 100.0/time_sample+1)
#print time_array
# Get thrust profile
thrust_array = get_thrust(time_array)
thrust_curve = interp1d(time_array, thrust_array, kind='linear')
xlabel('Time (s)')
ylabel('Thrust (N)', multialignment='center')
title('Thrust Profile')
pl.plot(time_array,thrust_array,time_array,thrust_curve(time_array))
Out[28]:
In [21]:
full_weight = 6500 #kg
jet_fuel = 500 #kg
rocket_fuel = 963+181 #kg
thrust_off_weight = -1.0
def get_weight(time_array):
weight_array = []
current_weight = full_weight
for t in time_array:
if t > 0.0 and t <= jet_end:
current_weight -= jet_fuel/jet_end * time_sample
if t >= rocket_start and t <= rocket_end:
current_weight -= rocket_fuel/(rocket_end-rocket_start) * time_sample
weight_array.append(current_weight)
return weight_array
# Get thrust profile
weight_array = get_weight(time_array)
weight_curve = interp1d(time_array, weight_array, kind='linear')
xlabel('Time (s)')
ylabel('Mass (kg)', multialignment='center')
title('Mass Profile')
pl.plot(time_array,weight_array,time_array,weight_curve(time_array))
savefig("mass_profile.png")
from IPython.display import FileLink, display
display(FileLink('mass_profile.png'))
In [8]:
def miles_to_km(miles):
return miles * 1.609344
def km_to_miles(km):
return km / 1.609344
def meters_per_second_to_mph(mps):
return mps*2.23693629
def mph_to_meters_per_second(mph):
return mph/2.23693629
end_of_measured_mile = miles_to_km(5.5) # km
start_of_measured_mile = end_of_measured_mile - miles_to_km(1.0) # km
# Calculate acceleration from a combination of drag and thrust
def get_force(current_time, current_speed, current_distance):
#
global thrust_off_weight
mach = current_speed/speed_of_sound
q = (0.5*gamma*pressure*mach*mach)
# Thrust
thrust = thrust_curve([current_time])
# Switch off jet at end of measured mile
if current_distance > end_of_measured_mile*1000.0:
thrust = 0.0
# Save weight
if thrust_off_weight < 0.0:
thrust_off_weight = weight_curve([current_time])
# Chute 1
if current_speed < mph_to_meters_per_second(579.0):
thrust += -0.9*q
# Chute 2
if current_speed < mph_to_meters_per_second(393.0):
thrust += -1.8*q
# Brakes
if current_speed < mph_to_meters_per_second(100.0):
thrust += -0.2*9.81*thrust_off_weight
# Aero Drag
drag_factor = 1.0
mach = current_speed/speed_of_sound
mach = min(mach,1.4)
drag = drag_curve([mach])[0]
drag *= q*drag_factor
#print thrust, drag
return thrust-drag
def verlet_integrator(r, v, dt, a):
"""Return new position and velocity from current values, time step and acceleration.
Parameters:
r is a numpy array giving the current position vector
v is a numpy array giving the current velocity vector
dt is a float value giving the length of the integration time step
a is a function which takes x as a parameter and returns the acceleration vector as an array
Works with arrays of any dimension as long as they're all the same.
"""
# Deceptively simple (read about Velocity Verlet on wikipedia)
r_new = r + v*dt + a(r)*dt**2/2
v_new = v + (a(r) + a(r_new))/2 * dt
return (r_new, v_new)
def get_performance(time_array):
acceleration_array = []
speed_array = []
distance_array = []
current_speed = 0.0
acceleration = 0.0
current_distance = 0.0
time_at_start = -1.0
time_at_end = -1.0
run_complete = False
# Euler integration
for t in time_array:
# Calculate current speed from previous acceleration
current_speed += acceleration*time_sample
#new_val = verlet_integrator(current_distance,current_speed,time_sample,acceleration)
#current_speed = new_val[1]
if current_speed < 0.0 or run_complete: # Stop the car
current_speed = 0.0
acceleration = 0.0
run_complete = True
else:
#current_speed = max(0.0,current_speed)
current_distance += current_speed*time_sample
#current_distance = new_val[0]
# a = F/m
if thrust_off_weight > 0.0:
weight = thrust_off_weight
else:
weight = weight_curve([t])
acceleration = get_force(t,current_speed, current_distance)/weight
if current_distance >= start_of_measured_mile*1000.0 and time_at_start < 0.0:
time_at_start = t
if current_distance >= end_of_measured_mile*1000.0 and time_at_end < 0.0:
time_at_end = t
acceleration_array.append(acceleration/9.81)
speed_array.append(meters_per_second_to_mph(current_speed))
distance_array.append(km_to_miles(current_distance/1000.0))
return acceleration_array, speed_array, distance_array, time_at_start, time_at_end
acceleration, speed, distance, timer_start, timer_end = get_performance(time_array)
#print timer_start, timer_end
print 'Average speed: %f mph'%(1.0/(timer_end-timer_start)*3600.0)
#pl.plot(time_array,speed)
#pl.plot(distance,speed)
#pl.plot(time_array,acceleration)
timer_start = [timer_start,timer_start]
timer_end = [timer_end, timer_end]
start_timing_line = [km_to_miles(start_of_measured_mile),km_to_miles(start_of_measured_mile)]
end_timing_line = [km_to_miles(end_of_measured_mile),km_to_miles(end_of_measured_mile)]
speed_timing_line = [0.0,1400]
fig = figure(figsize=(25, 10),dpi=100, facecolor='w', edgecolor='k')
ax = fig.add_subplot(1,3,1)
ax.grid(True)
xlabel('time (s)')
ylabel('speed (mph)', multialignment='center')
title('Speed - Time')
ax.plot(time_array, speed, color='r', label='speed (mph)')
ax.plot(timer_start,speed_timing_line,color='b')
ax.plot(timer_end,speed_timing_line,color='b')
ax = fig.add_subplot(1,3,2)
ax.grid(True)
xlabel('distance (miles)')
ylabel('speed (mph)', multialignment='center')
title('Speed - Distance')
ax.plot(distance, speed, color='r', label='speed (mph)')
ax.plot(start_timing_line,speed_timing_line,color='b')
ax.plot(end_timing_line,speed_timing_line,color='b')
ax = fig.add_subplot(1,3,3)
ax.grid(True)
xlabel('time (s)')
ylabel('acceleration (g)', multialignment='center')
title('Acceleration - Time')
ax.plot(time_array, acceleration, color='r', label='acceleration (g)')
Out[8]:
In [ ]:
In [ ]: