Investigation of the flow field from a serrated nozzle

#### According to measurements from NASA ##### Dj = 0.0494 ------------ Jet Diameter ##### Uid = 300.0 ------------ Reference Velocity ##### Mjet = 0.9 ------------ Jet Mach Number ##### Setpoint = 007 ------------ According to Tanna (1987) ##### NozID = SMC006 ------------ According to NASA -------- ##### Pamb = 97000.0 ------------ Ambient pressure ##### Tamb = 280.2 ------------ Ambient temperature ##### NozzlePressureRatio =1.83 ##### NozzleTemperatureRatio = 1.022 ## References ##### Hao Xia, Paul G. Tucker and Simon Eastwood (2009). Large-eddy simulations of chevron jet flows with noise predictions. International Journal of Heat and Fluid Flow 30 (2009) 1067-1079. ###### Hao Xia and Paul G. Tucker (2011). Numerical Simulation of a Single-Stream Jets from a Serrated Nozzle . Flow Turbulence Combust 2011. ##### More to follow

Define case name

This is the solver case to be analysed

In [1]:
case_name = 'smc006'

Define Data Location

For remote data the interaction will use ssh to securely interact with the data
This uses the reverse connection capability in paraview so that the paraview server can be submitted to a job scheduler
Note: The default paraview server connection will use port 11111

In [2]:
remote_data = True

data_dir='/gpfs/cfms/workarea/projects/hyper_flux/serrated_nozzle/smc_zCFD/smc_jamil'
data_host='[email protected]'

remote_server_auto = True

paraview_cmd='mpiexec /gpfs/cfms/apps/bin/pvserver'
if not remote_server_auto:
    paraview_cmd=None

if not remote_data:
    data_host='localhost'
    paraview_cmd=None

Initialise Environment

In [3]:
%pylab inline
from paraview.simple import *
paraview.simple._DisableFirstRenderCameraReset()
import pylab as pl
import math
import numpy as np
Populating the interactive namespace from numpy and matplotlib
paraview version 4.2.0-75-gff3f889

Data Connection

This starts paraview server on remote host and connects

In [4]:
from zutil.post import pvserver_connect
if remote_data:
    pvserver_connect(data_host=data_host,data_dir=data_dir,paraview_cmd=paraview_cmd)
[[email protected]] Executing task 'port_test'
Selected Port: 12000
[[email protected]] Executing task 'pvserver'
[[email protected]] run: /bin/bash -l -c "cd /gpfs/cfms/workarea/projects/hyper_flux/serrated_nozzle/smc_zCFD/smc_jamil && sleep 2;mpiexec /gpfs/cfms/apps/bin/pvserver -rc --client-host=localhost -sp=12000"
[[email protected]] out: 
[[email protected]] out: 		   _____ ______ __  __  _____ 
[[email protected]] out: 		  / ____|  ____|  \/  |/ ____|
[[email protected]] out: 		 | |    | |__  | \  / | (___  
[[email protected]] out: 		 | |    |  __| | |\/| |\___ \ 
[[email protected]] out: 		 | |____| |    | |  | |____) |
[[email protected]] out: 		  \_____|_|    |_|  |_|_____/ 
[[email protected]] out: 
[[email protected]] out: ++++++++++++++++++++++++++++: System Data :++++++++++++++++++++++++++++
[[email protected]] out: + Hostname = vis03
[[email protected]] out: + Kernel = 2.6.32-358.el6.x86_64
[[email protected]] out: + RHEL Release = Red Hat Enterprise Linux Server release 6.4 (Santiago)
[[email protected]] out: + Uptime = 13:44:26 up 80 days, 3:29, 34 users,
[[email protected]] out: ++++++++++++++++++++++++++++: User Data :++++++++++++++++++++++++++++++
[[email protected]] out: + Username = acimpoeru
[[email protected]] out: +++++++++++++++++++++++: Contact Information :+++++++++++++++++++++++++
[[email protected]] out: + in case of any problems, contact: [email protected]
[[email protected]] out: + for feedback, contact: [email protected] 
[[email protected]] out: +++++++++++++++++++++: Maintenance Information :+++++++++++++++++++++++
[[email protected]] out: + There is no planned maintenance taking place this week
[[email protected]] out: +++++: Group Home Directory Quota Usage (updated every 10 mins)   +++++
[[email protected]] out: + 
[[email protected]] out: +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
[[email protected]] out: Reported: 1 (out of 1) daemons - 1 (out of 1) procs
[[email protected]] out: [[email protected]] rtunnel: opened reverse tunnel: (u'127.0.0.1', 53355) -> ('172.20.1.25', 22) -> ('localhost', 11111)
Connecting to client (reverse connection requested)...
[[email protected]] out: Connection URL: csrc://localhost:12000
[[email protected]] out: Client connected.
[[email protected]] out: 

Get control dictionary

In [5]:
from zutil.post import get_case_parameters,print_html_parameters
parameters=get_case_parameters(case_name,data_host=data_host,data_dir=data_dir)

Get status file

In [6]:
from zutil.post import get_status_dict
status=get_status_dict(case_name,data_host=data_host,data_dir=data_dir)
num_procs = str(status['num processor'])

Define Jet Conditions

In [7]:
# AMBIENT CONDITIONS
P_amb          =     97000 # [Pa]
T_amb          =     280.2 # [K]
gas_constant   =     287.0
Rho_amb        =     P_amb/(gas_constant * T_amb) # [kg/m^3] - check
mu             =     1.79e-5 # viscosity [Pa.s]
gamma          =     1.4
speed_of_sound =     math.sqrt(gamma*P_amb/Rho_amb) # [m/s]

# JET CONDITIONS 
NPR            =     1.83     # nozzle pressure ratio    P_jet/P_amb
NTR            =     1.022    # nozzle temperature ratio T_jet/T_amb
Mjet           =     0.9      # jet mach number 
Djet           =     49.4     # chevron tip to tip (exit diameter) [mm] 
Ujet           =     300.0    # reference velocity 
Re             =     1.03e6   # based  on  jet exit velocity , jet diameter and ambient conditions for mu and rho
 
from IPython.display import HTML
HTML(print_html_parameters(parameters))
Out[7]:
pressure97000.0
temperature280.2
Reynolds Noundefined
Ref length0.001
Speed0.0
Mach No0.01

Read experimental data

In [8]:
import zutil
import zutil.post as post
from zutil.post import get_csv_data
import os
import glob
from collections import OrderedDict


# TOP FRAME ------ CUT PLANES AT THE SERRATIONS TIPS 
v1 = zutil.post.get_csv_data("/gpfs/cfms/workarea/projects/hyper_flux/serrated_nozzle/smc_zCFD/exp_data/top_frame/data_x_0_25.csv",
                             header=True,remote=True,delim=',')
v2 = zutil.post.get_csv_data("/gpfs/cfms/workarea/projects/hyper_flux/serrated_nozzle/smc_zCFD/exp_data/top_frame/data_x_0_5.csv",
                             header=True,remote=True,delim=',')
v3 = zutil.post.get_csv_data("/gpfs/cfms/workarea/projects/hyper_flux/serrated_nozzle/smc_zCFD/exp_data/top_frame/data_x_1_0.csv",
                             header=True,remote=True,delim=',')
v4 = zutil.post.get_csv_data("/gpfs/cfms/workarea/projects/hyper_flux/serrated_nozzle/smc_zCFD/exp_data/top_frame/data_x_2_5.csv",
                             header=True,remote=True,delim=',')
v5 = zutil.post.get_csv_data("/gpfs/cfms/workarea/projects/hyper_flux/serrated_nozzle/smc_zCFD/exp_data/top_frame/data_x_5_0.csv",
                             header=True,remote=True,delim=',')
v6 = zutil.post.get_csv_data("/gpfs/cfms/workarea/projects/hyper_flux/serrated_nozzle/smc_zCFD/exp_data/top_frame/data_x_7_5.csv",
                             header=True,remote=True,delim=',')
# notch plane

n1 = zutil.post.get_csv_data("/gpfs/cfms/workarea/projects/hyper_flux/serrated_nozzle/smc_zCFD/exp_data/notch_frame/data_x_D=0_25.csv",
                             header=True,remote=True,delim=',')
n2 = zutil.post.get_csv_data("/gpfs/cfms/workarea/projects/hyper_flux/serrated_nozzle/smc_zCFD/exp_data/notch_frame/data_x_D=0_5.csv",
                             header=True,remote=True,delim=',')
n3 = zutil.post.get_csv_data("/gpfs/cfms/workarea/projects/hyper_flux/serrated_nozzle/smc_zCFD/exp_data/notch_frame/data_x_D=1_0.csv",
                             header=True,remote=True,delim=',')
n4 = zutil.post.get_csv_data("/gpfs/cfms/workarea/projects/hyper_flux/serrated_nozzle/smc_zCFD/exp_data/notch_frame/data_x_D=2_5.csv",
                             header=True,remote=True,delim=',')
n5 = zutil.post.get_csv_data("/gpfs/cfms/workarea/projects/hyper_flux/serrated_nozzle/smc_zCFD/exp_data/notch_frame/data_x_D=5_0.csv",
                             header=True,remote=True,delim=',')
n6 = zutil.post.get_csv_data("/gpfs/cfms/workarea/projects/hyper_flux/serrated_nozzle/smc_zCFD/exp_data/notch_frame/data_x_D=7_5.csv",
                             header=True,remote=True,delim=',')


station_values = OrderedDict([("S01" , 0.25* Djet), 
                              ("S02" , 0.5 * Djet), 
                              ("S03" , 1.0 * Djet), 
                              ("S04" , 2.5 * Djet),
                              ("S05",  5.0 * Djet),
                              ("S06",  7.5 * Djet)])
exp_values = OrderedDict([("E01" , (0.25,v1[0],v1[1])), 
                          ("E02" , (0.5,v2[0],v2[1])), 
                          ("E03" , (1.0,v3[0],v3[1])), 
                          ("E04" , (2.5,v4[0],v4[1])),
                          ("E05",  (5.0,v5[0],v5[1])),
                          ("E06",  (7.5,v6[0],v6[1]))])

exp_values_n = OrderedDict([("N01" , (0.25,n1[0],n1[1])),
                            ("N02" , (0.5,n2[0],n2[1])), 
                            ("N03" , (1.0,n3[0],n3[1])), 
                            ("N04" , (2.5,n4[0],n4[1])),
                            ("N05",  (5.0,n5[0],n5[1])),
                            ("N06",  (7.5,n6[0],n6[1]))])


#isinstance(station_values['S01'], float)
#print n1.keys()
//Applications/paraview.app/Contents/Python/paraview/vtk/numpy_interface/dataset_adapter.py:126: FutureWarning: comparison to `None` will result in an elementwise object comparison in the future.
  if array == None:
In [9]:
from zutil.post import for_each

def plot_velocity_profile(ax,file_root):
    
    smc006 = PVDReader( FileName=file_root+'_0'+'.pvd' )
    smc006.UpdatePipeline()
    
    slice1 = Slice( SliceType="Plane", Input=smc006 )
    slice1.SliceType.Normal = [0.0, 0.0, 1.0]
    slice1.SliceType.Origin = [0.0, 0.0, 0.0]    
    slice1.UpdatePipeline()
    merge = MergeBlocks(Input=slice1)
    merge.UpdatePipeline()
    
    point_data = CellDatatoPointData(Input=merge)
    point_data.PassCellData = 0
    point_data.UpdatePipeline()
    count = 0    
    
    for station in station_values:
        slice2 = Slice( SliceType="Plane", Input=point_data )
        slice2.SliceType.Normal = [1.0, 0.0, 0.0]
        slice2.SliceType.Origin = [station_values[station], 0.0, 0.0]    
        slice2.UpdatePipeline()
         
        #print point_data
        slice_client = servermanager.Fetch(slice2)
        for_each(slice_client,func=plot_array,axis=ax,station=count)
        count +=1
        print 'I am working on   ', station 
    
def plot_array(data_array,pts_array,**kwargs):
    ax = kwargs['axis']
    #loc = kwargs['loc']
    station=kwargs['station']
    data = []
    vel = 0
    axis_array = pts_array.GetPoints()[:,1]/49.4
    vel_array = data_array.GetPointData()['V'][:,0]/300.0 + station
    ax.plot(vel_array, axis_array,'b.',label='zCFD'+ '-station' + str(station+1),markersize=20)
    #ax.legend(loc = 'upper left',numpoints=1,prop = prop)
    
    
In [10]:
def plot_velocity_profile_notch(ax,file_root):
    
    smc006_n = PVDReader( FileName=file_root+'_0'+'.pvd' )
    smc006_n.UpdatePipeline()
    
    slice3 = Slice( SliceType="Plane", Input=smc006_n )
    slice3.SliceType.Normal = [0.0, 1.0, 0.0]
    slice3.SliceType.Origin = [0.0, 0.0, 0.0]    
    slice3.UpdatePipeline()
    merge2 = MergeBlocks(Input=slice3)
    merge2.UpdatePipeline()
    
    point_data2 = CellDatatoPointData(Input=merge2)
    point_data2.PassCellData = 0
    point_data2.UpdatePipeline()
    count2 = 0    
    
    for station in station_values:
        slice4 = Slice( SliceType="Plane", Input=point_data2 )
        slice4.SliceType.Normal = [1.0, 0.0, 0.0]
        slice4.SliceType.Origin = [station_values[station], 0.0, 0.0]    
        slice4.UpdatePipeline()
         
        #print point_data
        slice_client2 = servermanager.Fetch(slice4)
        for_each(slice_client2,func=plot_array2,axis=ax,station=count2)
        count2 +=1
        print 'I am working on   ', station 
    
def plot_array2(data_array,pts_array,**kwargs):
    ax = kwargs['axis']
    #loc = kwargs['loc']
    station=kwargs['station']
    #data = []
    #vel = 0
    axis_array2 = pts_array.GetPoints()[:,2]/49.4
    vel_array2 = data_array.GetPointData()['V'][:,0]/300.0 + station
    ax.plot(vel_array2, axis_array2,'b.',label='zCFD'+ '-station' + str(station+1),markersize=20)
    #ax.legend(loc = 'upper left',numpoints=1,prop = prop)

Velocity profiles (tip-to-tip) respectively notch planes

In [11]:
from zutil.post import get_case_root
import matplotlib.font_manager as fm
from matplotlib import colors
prop=fm.FontProperties(size=20)

#             TIP TO TIP PLANE
fig = pl.figure(figsize=(30, 10),dpi=300, facecolor='w', edgecolor='k')
fig.suptitle('$SMC006 $ $Velocity $ $Profiles$ $y-z plane $', fontsize=35, fontweight='bold', style='italic')
ax = fig.add_subplot(1,1,1)
ax1=fig.add_subplot(1,1,1)
ax.grid(True)
ax.set_xlabel('$U$ / $U_j$ $[ ]$',fontsize=35,fontweight='bold',style='italic')
ax.set_ylabel('$R$ / $D_j$ $[ ]$',fontsize=35,fontweight='bold',style='italic')
ax.axis([0.0,7.0,0.0,1.0])  
plot_velocity_profile(ax,get_case_root(case_name,num_procs))
#ax.legend(loc='upper right', shadow=True,prop=prop)
colors  =[["r"],["g"],["m"],["c"],["b"]]
for key in exp_values:
    ax1.plot(exp_values[key][1],exp_values[key][2],'b.',color='r',label='Exp x/D = ' + str(exp_values[key][0]), markersize=20)
ax.legend(loc = 'upper right',numpoints=1,shadow=True, prop=prop)


#             NOTCH PLANE
fig2= pl.figure(figsize=(30, 10),dpi=300, facecolor='w', edgecolor='k')
fig2.suptitle('$SMC006 $ $Velocity $ $Profiles$ $x-z plane $', fontsize=35, fontweight='bold', style='italic')
ax2 = fig2.add_subplot(1,1,1)
ax3=fig2.add_subplot(1,1,1)
ax2.grid(True)
ax2.set_xlabel('$U$ / $U_j$ $[ ]$',fontsize=35,fontweight='bold',style='italic')
ax2.set_ylabel('$R$ / $D_j$ $[ ]$',fontsize=35,fontweight='bold',style='italic')
ax2.axis([0.0,7.0,0.0,2.0])
plot_velocity_profile_notch(ax2,get_case_root(case_name,num_procs))
colors  =[["r"],["g"],["m"],["c"],["b"]]
for key in exp_values_n:
    ax2.plot(exp_values_n[key][1],exp_values_n[key][2],'b.',color='r',label='Exp x/D = ' + str(exp_values_n[key][0]), markersize=20)
ax2.legend(loc = 'upper right',numpoints=1,shadow=True, prop=prop)

# Save the figures 

from matplotlib.backends.backend_pdf import PdfPages
pp1 = PdfPages('SMC006_TIP_TO_TIP.pdf')
pp1.savefig()
pp1.close()
fig.savefig("SMC006_TIP_TO_TIP.png")
#show()
pp2 = PdfPages('SMC006_NOTCH.pdf')
pp2.savefig()
pp2.close()
fig2.savefig("SMC006_NOTCH.png")
show()

from IPython.display import FileLink, display 
display(FileLink('SMC006_TIP_TO_TIP.png'))
display(FileLink('SMC006_NOTCH.png'))


    
    
    
    
    
I am working on    S01
I am working on    S02
I am working on    S03
I am working on    S04
I am working on    S05
I am working on    S06
I am working on    S01
I am working on    S02
I am working on    S03
I am working on    S04
I am working on    S05
I am working on    S06

Residuals L2 norms

In [22]:
from zutil.post import residual_plot, get_case_report
residual_plot(get_case_report(case_name))
show()

Cleaning up

In [12]:
if remote_data:
    #print 'Disconnecting from remote paraview server connection'
    Disconnect()
    pass
Exiting...
[[email protected]] out: