MB CNS Example Book. Mechanical Engineering Report 2004/10 P. A. Jacobs Centre for Hypersonics The University of Queensland. September 2004, Revised March 2005 Preface MB CNS is a collection of programs for the simulation of transient two-dimensional (or axisymmetric) flows. It is part of the larger collection of compressible flow simulation codes found at http://www.mech.uq.edu.au/cfcfd/. This manual is a collection of example simulations: scripts, results and commentary. It may be convenient for new users of the code to identify an example close to the situation that they wish to model and then adapt the scripts for that example. This report will be updated occasionally with new examples and commentary. As the simulation codes are improved, we will try to maintain compatibility so that older examples are not “broken”, however, this backward compatibility will not be guaranteed.

1

Contents 1 Introduction

3

2 Mach 1.5 flow over a 20o 2.1 .sit file . . . . . . . . . 2.2 Shell scripts . . . . . . 2.3 Notes . . . . . . . . . . 3 Flow of equilibrium-air 3.1 .sit file . . . . . . . . 3.2 Shell scripts . . . . . 3.3 Notes . . . . . . . . . 4 Hypersonic flow of 4.1 .sit file . . . . . 4.2 Shell scripts . . 4.3 Notes . . . . . . 5 Mach 3 flow over 5.1 .sit file . . . . 5.2 Shell scripts . 5.3 Notes . . . . . 6 Flow through a 6.1 .sit file . . . 6.2 Shell scripts 6.3 Notes . . . .

cone 5 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10

over . . . . . . . . .

a sphere . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

ideal air over . . . . . . . . . . . . . . . . . . . . . . . . . . .

a . . .

blunt . . . . . . . . . . . .

12 12 15 16

wedge 18 . . . . . . . . . . . . . . . . . . . 21 . . . . . . . . . . . . . . . . . . . 23 . . . . . . . . . . . . . . . . . . . 24

a sharp-nosed two-dimensional . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

body 26 . . . . . . . . . . . . . . 29 . . . . . . . . . . . . . . 30 . . . . . . . . . . . . . . 31

conical nozzle 32 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38

7 A section of an ideal compressible-flow 7.1 .sit file . . . . . . . . . . . . . . . . . . 7.2 Shell scripts . . . . . . . . . . . . . . . 7.3 Notes . . . . . . . . . . . . . . . . . . . 8 Pressure on a flat-faced 8.1 .sit file . . . . . . . . 8.2 Shell scripts . . . . . 8.3 Notes . . . . . . . . .

vortex . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

40 41 43 44

cylinder . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

46 48 49 51

2

1

Introduction

Setting up a simulation is mostly an exercise in writing a textual description of your flow and its bounding geometry. This is presented to the scriptit program as a text file, sometimes referred to as a “.sit” file or flow-specification file. Once you have prepared your flow specification file, the simulation data is generated in a number of stages: 1a Create the geometry definition with the command. $ scriptit.tcl -f job.sit -do-mpost Input:

Program: -

job.sit

Output:

scriptit.tcl

-

job.p

-

job.bez

-

job.mpost

1b Check the geometry definition (visually) by using Metapost to make a viewable postscript file containing labelled nodes, block boundary curves and blocks. Metapost is distributed as part of the TEX document preparation system. It is most likely already installed on your UNIX/Linux system and there is a stand-alone binary for Win32 systems. $ mpost job.mpost -

job.mpost

-

mpost

job.1

2 Generate a grid and initial flow solution (at t = 0). $ mb prep.exe -f job job.p

6

mb prep.exe

job.bez

-

job.g

-

job.s0

3 Run the simulation code to produce flow data at subsequent times. $ mb cns.exe -f job

3

job.p

job.g

6

mb cns.exe

6

-

job.s

-

job.h

job.s0 4 Extract subsets of the flow solution data for postprocessing. The specific command for this stage depends very much on what you want to do. The flow solution data is cell-averaged data associated with cell centres. You may extract the flow data for all cells at a particular time using mb post.exe and reformat it for a particular plotting program or you may extract data along single grid lines (using mb prof.exe) in a form ready for display with GNU-Plot or for further calculation. See the shell scripts in the examples for ideas on what can be done. Since the output of this stage is always a text file, you may look at the head of the file for hints as to what data is present. Originally, scriptit was a C program which read your script, initialized some data structures and wrote the parameter (job.p) and geometry definition (job.bez) files. Currently, scriptit is implemented as a TCL program (see e.g. Reference [1]) that has a few extra procedures and your specification script is really a TCL script. It is worth the effort to learn just enough Tcl to be dangerous. The Web site http://www.tcl.tk has a good (short) starting guide. After doing some initialization, scriptit(.tcl) sources your script file and assembles the geometry and flow specification data into a form that can be given to the main simulation codes (mb prep.exe and mb cns.exe). The advantage of this approach is that you have the full capability of the Tcl interpreter available to you from within your script. You can perform calculations so that you can parameterize your geometry, for example, or you can use Tcl control structures to make repetitive definitions much more concise. Additionally, you may use Tcl comments to add documentation to the script file. The following examples should be studied in together with the HTML reference for mb cns programs, and scriptit in particular. These hypertext manuals can be found in the documentation section at the URL: http://www.mech.uq.edu.au/cfcfd/.

4

2

Mach 1.5 flow over a 20o cone

This is a small (in both memory and run time) example that is useful for checking that the simulation and plotting programs have been built or installed correctly. Assuming that you have the program executable files built and accessible on your system’s search PATH, try the following commands: $ cd ∼/cfcfd/code/mb cns/examples/cone20 $ ./cone20 run.sh $ ./cone20 plot.sh And, within a minite or so, you should end up with a number of files with various solution data plotted. The grid and initial solution are created and the time-evolution of the flow field is computed for 5 ms (with 1105 time steps being required). The commands invoke the shell scripts displayed in subsection 2.2. NORTH0

1

NORTH1

F

E

D

EAST1

0.8

0.6 y

BLOCK-0

EAST0

WEST0

BLOCK-1

0.4

C 0.2

0

1

TH

SOU

SOUTH0 A 0

B 0.2

0.4

0.6

0.8

1

x

Figure 1: Schematic diagram of the geometry for a cone with 20 degree half-angle. This encapsulated postscript figure was generated from the Metapost file. The free-stream conditions (p∞ = 95.84 kPa, T∞ = 1103 K and u∞ = 1000 m/s) are related to the shock-over-ramp test problem in the original ICASE Report [2] and are set to give a Mach number of 1.5. From Chart 5 in Ref. [3], the expected steady-state shock 5

wave angle is 49o and, from Chart 6, the pressure coefficent is pcone−surf ace − p∞ ≈ 0.387 q∞ and the dynamic pressure for the specified free stream is q∞ = 21 ρ∞ u2∞ ≈ 151.38 kPa. Figure 4 shows the pressure coefficient estimated as Cp =

fx − p ∞ A q∞ A

from the simulated axial force, fx , written into the simulation log file and frontal area of the cone, A.

Figure 2: Pressure data for flow over a cone with 20 degree half-angle. The shock profile is not yet straight and the pressure field near the cone surface is not conically symmetric, although it would become more if we continued the simulation.

6

Figure 3: Shock-sensor data for flow over a cone with 20 degree half-angle. For the adaptive flux calculator, this sensor indicates the regions of the flow where the more dissipative scheme should be used. 20 degree cone in Mach 1.5 flow 0.5

average Cp

0.4

Value from NACA 1135 Chart 6

0.3

0.2

0.1 10x40+30x40 20x80+60x80

0 0

1

2

3 time, ms

4

5

6

Figure 4: Evolution of the axial (drag) force for flow over a cone with 20 degree half-angle for two mesh resolutions. 7

2.1

.sit file

# cone20 . s i t # Mach 1 . 5 f l o w o v e r a 2 0 d e g r e e cone . # S e t up two q u a d r i l a t e r a l s i n t h e ( x , y)− p l a n e be f i r s t d e f i n i n g # t h e c o r n e r nodes , then t h e l i n e s between t h o s c o r n e r s and then # t h e boundary e l e m e n t s f o r t h e b l o c k s . BEGIN GEOMETRY NODE a 0 . 0 0 . 0 NODE b 0 . 2 0 . 0 NODE c 1 . 0 0 . 2 9 1 1 8 NODE d 1 . 0 1 . 0 NODE e 0 . 2 1 . 0 NODE f 0 . 0 1 . 0 LINE LINE LINE LINE LINE LINE LINE

ab bc af be cd fe ed

a b a b c f e

# Define POLYLINE POLYLINE POLYLINE POLYLINE POLYLINE POLYLINE POLYLINE END GEOMETRY

b c f e d e d the boundaries . north0 1 + f e east0 1 + be s o u t h 0 1 + ab west0 1 + af s o u t h 1 1 + bc east1 1 + cd n o r t h 1 1 + ed

BEGIN FLOW # Gas and f l o w p r o p e r t i e s GAS TYPE p e r f a i r 1 4 GAS STATE i n i t i a l 5 9 5 5 . 0 0.0 0.0 304.0 9 5 . 8 4 e3 1 0 0 0 . 0 0 . 0 1 1 0 3 . 0 GAS STATE i n f l o w

1.0 1.0

# S e t t h e boundary d i s c r e t i s a t i o n b e f o r e b u i l d i n g t h e b l o c k s . DISCRETISE n o r t h 0 1 0 0 0 0 . 0 DISCRETISE e a s t 0 40 0 0 0.0 DISCRETISE s o u t h 0 1 0 0 0 0 . 0 DISCRETISE west0 40 0 0 0.0 DISCRETISE n o r t h 1 3 0 0 0 0 . 0 DISCRETISE s o u t h 1 3 0 0 0 0 . 0 DISCRETISE e a s t 1 40 0 0 0.0 # I n f l o w and o u t f l o w b o u n d a r i e s . BOUNDARY SPEC west0 SUP IN i n f l o w BOUNDARY SPEC e a s t 1 EXTRAPOLATE OUT # D e f i n e two b l o c k s with a common boundary . # Although we have used i n t e g e r s i n t h e b l o c k names , # any u n i q u e names would be f i n e . # Same g o e s f o r t h e o t h e r e n t i t i e s such a s nodes , and b o u n d a r i e s . BLOCK b l o c k −0 + n o r t h 0 + e a s t 0 + s o u t h 0 + west0 BLOCK b l o c k −1 + n o r t h 1 + e a s t 1 + s o u t h 1 + e a s t 0 CONNECT BLOCKS b l o c k −0 e a s t b l o c k −1 west # Have s p e c i f i e d an a r e a −o r t h o g o n a l i t y g r i d f o r b l o c k −1 , # j u s t f o r fun . GRID TYPE b l o c k −1 AO # Assign the i n i t i a l gas s t a t e s FILL BLOCK b l o c k −0 i n i t i a l FILL BLOCK b l o c k −1 i n i t i a l END FLOW

8

BEGIN CONTROL TITLE Mach 1 . 5 f l o w o v e r a 2 0 d e g r e e cone . CASE ID 5 AXISYMMETRIC VISCOUS FLUX CALC a d a p t i v e MAX TIME 5 . 0 e−3 MAX STEP 3000 TIME STEP 1 . 0 e−6 DT PLOT 1 . 5 e−3 DT HISTORY 1 0 . 0 e−5 HISTORY CELL b l o c k − 1 1 0 1 END CONTROL # Name t h e o u t p u t f i l e s and b u i l d them . MPOST FILE c o n e 2 0 . mpost MPOST SCALES 0 . 1 2 0 . 1 2 MPOST XAXIS 0 . 0 1 . 0 0 . 2 − 0 . 0 5 MPOST YAXIS 0 . 0 1 . 0 0 . 2 − 0 . 0 4 BEZIER FILE c o n e 2 0 . bez PARAM FILE c o n e 2 0 . p BUILD QUIT

2.2

Shell scripts

#! / b i n / sh # c o n e 2 0 r u n . sh # e x e r c i s e t h e Navier −S t o k e s s o l v e r f o r t h e Cone20 t e s t c a s e . # I t i s assumed t h a t t h e path i s s e t c o r r e c t l y . # G e n e r a t e t h e B e z i e r , I n p u t p a r a m e t e r and MetaPost f i l e s from t h e S c r i p t F i l e . # The MetaPost f i l e p r o v i d e s us with a g r a p h i c a l c h e c k on t h e geometry # s p e c i f i c a t i o n and i t i s worth c r e a t i n g i f metapost ( mpost ) i s a v a i l a b l e . i f [ [ − f / usr / bin / t c l s h ] ] then echo ” Use t h e new s c r i p t i t . ” s c r i p t i t . t c l − f c o n e 2 0 . s i t −do−mpost > c o n e 2 0 . s c r i p t i t −l o g i f [ [ − f c o n e 2 0 . mpost ] ] then mpost c o n e 2 0 . mpost fi else echo ” Use t h e o l d s c r i p t i t . ” s c r i p t i t . exe < cone20 . s i t > cone20 . log fi # G e n e r a t e t h e Grid and I n i t i a l S o l u t i o n F i l e s . mb prep . e x e − f c o n e 2 0 # I n t e g r a t e t h e s o l u t i o n i n time , # r e c o r d i n g t h e a x i a l f o r c e on t h e cone s u r f a c e . # # The f o l l o w i n g e n v i r o n m e n t v a r i a b l e s a l l o w t h e s h a r e d−memory v e r s i o n # o f t h e code t o u s e one t h r e a d f o r each o f t h e two b l o c k s . # The s t a c k s i z e r e q u i r e m e n t s may i n c r e a s e a s t h e code d e v e l o p s and # more e l e m e n t s a r e added t o t h e i n t e r n a l data s t r u c t u r e s . e x p o r t OMP NUM THREADS=2 e x p o r t KMP STACKSIZE=8m time mb cns . e x e − f c o n e 2 0 − f o r c e 1 2 # E x t r a c t t h e s o l u t i o n data and r e f o r m a t . # I f no time i s s p e c i f i e d , t h e f i r s t s o l u t i o n found i s o u t p u t . mb post . e x e − f p c o n e 2 0 . p − f g c o n e 2 0 . g − f s c o n e 2 0 . s − f o c o n e 2 0 − g e n e r i c

9

# E x t r a c t t h e a v e r a g e c o e f f i c i e n t o f p r e s s u r e from t h e a x i a l f o r c e # r e c o r d s t h a t were w r i t t e n t o t h e s i m u l a t i o n l o g f i l e . awk − f cp . awk mb cns . l o g > c o n e 2 0 c p . dat echo ” At t h i s p o i n t , we s h o u l d have a new s o l u t i o n ” echo ” Run c o n e 2 0 p l o t . sh n e x t ”

# c o n e 2 0 p l o t . sh # P i c k up t h e r e f o r m a t t e d data and make p l o t s o f : # 1 . Shock i n d i c a t o r mb cont . e x e − f i c o n e 2 0 . gen − f o c o n e 2 0 S . ps − v a r 9 − ps − c o l o u r \ −x r a n g e 0 . 0 1 . 0 0 . 5 − y r a n g e − 0 . 0 1 . 0 0 . 5 − mesh # 2. Pressure . mb cont . e x e − f i c o n e 2 0 . gen − f o c o n e 2 0 p . ps − v a r 6 − ps −m i r r o r − x r a n g e 0 . 0 1 . 0 0 . 5 − y r a n g e − 1 . 0 mb cont . e x e − f i c o n e 2 0 . gen − f o c o n e 2 0 . g i f − v a r 6 − g i f −m i r r o r − x r a n g e 0 . 0 1 . 0 0 . 5 − y r a n g e − 1 . 0

−colour \ 1.0 0.5 −colour \ 1.0 0.5

# 3 . The mesh a l o n e . mb cont . e x e − f i c o n e 2 0 . gen − f o cone20 mesh . ps − v a r 6 − ps − c o l o u r \ −x r a n g e 0 . 0 1 . 0 0 . 5 − y r a n g e − 0 . 0 1 . 0 0 . 5 − mesh − n o c o n t o u r s # 4 . The a v e r a g e c o e f f i c i e n t o f p r e s s u r e on t h e cone s u r f a c e . # We assume t h a t t h e high−r e s o l u t i o n data f i l e i s a l s o a v a i l a b l e . g n u p l o t s s 3 . r e s u l t

3.3

Notes

• This simulation reaches a final time of 67.95 µs in 5192 steps and, on a Celeron 2.4 Ghz system, this takes 7 min, 53 s of CPU time. This timing will be a bit sensitive to the state of the code because the large data structures appear to be causing a lot of cache misses. If we cut down on the amount of storage for each cell and reduce the size of the temporary arrays, we can achieve significant reductions in the CPU time. • Awk script for extracting the shock location from the stagnation-line flow data. # l o c a t e s h o c k . awk

BEGIN { p old = 0.0; x old = −2.0; # dummy p o s i t i o n y old = −2.0; p t r i g g e r = 2 . 0 e6 ; # s o m e t h i n g midway between f r e e s t r e a m and s t a g n a t i o n shock found = 0; }

16

$1 ! = ” # ” { # f o r any non−comment l i n e , do s o m e t h i n g p new = $7 ; x new = $1 ; y new = $2 ; # p r i n t ” p new =” , p new , ” x new ” , x new , ” y new ” , y new i f ( p new > p t r i g g e r && s h o c k f o u n d = = 0 ) { shock found = 1; f r a c = ( p new − p t r i g g e r ) / ( p new − p o l d ) ; x = x o l d + f r a c ∗ ( x new − x o l d ) ; y = y o l d + f r a c ∗ ( y new − y o l d ) ; p r i n t ” shock−l o c a t i o n = ” , x , y } p o l d = p new ; x o l d = x new ; y o l d = y new ; } END { i f ( shock found == 0 ) { p r i n t ” s h o c k not l o c a t e d ” ; } p r i n t ” done . ” }

17

4

Hypersonic flow of ideal air over a blunt wedge

This example is a partial solution to the CFD exercise for the MECH4470 class in 2004. Because the original specification was given in nondimensional form, an arbitrary 10 mm nose radius has been selected for the inviscid simulation. This is also a reasonable size for a possible wind tunnel experiment. The free-stream condition was specified as having a Mach number of 5 and the gas was specified as ideal air. Choosing particular values of p∞ = 100 kPa, T∞ = 100 K, lead to a free-stream velocity of u∞ = 1002 m/s and a dynamic pressure of q∞ = 1.75 MPa. G

E1

N1

BLK-1

F N0W1

D

N0W1

W0

F1

S1

BLK-0 E1 E0

C S0

E

B

A

Figure 8: Schematic diagram of the geometry for the blunted 10 degree wedge. The simulation is started with low-pressure conditions throughout the flow domain and free-stream conditions applied to the inflow boundary (the west boundary of blk-1 and the north boundary of blk-1). The flow data is allowed to evolve until tf inal = 399 µs, which corresponds to a particle of the free-stream travelling 40 nose radii. The axial force (shown in Fig.10) is seen to settle to a value of 28260 N in that time. This corresponds to a drag coefficient of 0.666. The surface pressure (shown normalised in Fig. 11) has been extracted from the solution file by mb prof by selecting the east-most line of cells of the first block and the south-most line of cells of the second block. The selected data is filtered by an Awk script to produce the normalised data (and the Newtonian reference data) as plotted.

18

Figure 9: Mesh for the blunt wedge exercise.

Blunted wedge: x-force history 35000

total cylinder wedge

30000

x-force, N

25000 20000 15000 10000 5000 0 0

50

100

150

200

250

300

350

400

t, microseconds

Figure 10: History of the axial forces for the blunt-wedge exercise.

19

Blunted wedge: surface pressure coefficient. Pressure Coefficient, (p - p_inf)/q_inf

2

CFD at t=399us Modified Newtonian

1.5

1

0.5

0 0

1

2

3

4

5

6

7

8

9

10

s/Rn

Figure 11: Surface pressure coefficient data for the blunt-wedge exercise.

Figure 12: Mach number data for the blunt-wedge exercise.

20

4.1

.sit file

# bw . s i t # MECH4470/CFD E x e r c i s e : H y p e r s o n i c f l o w o v e r a b l u n t wedge . # PJ , 06 − Oct−04 # Geometry s e t Rn 10.0 e −3; s e t xEnd [ e x p r 8 . 0 ∗ $Rn ] ; s e t alpha [ t o r a d i a n s 1 0 . 0 ] ; s e t d e l t a 1 0 . 0 e −3;

# # # #

radius of c y l i n d r i c a l nose downstream e x t e n t o f wedge a n g l e o f wedge wrt f r e e s t r e a m o f f s e t f o r i n f l o w boundary

# Free stream # I d e a l Air set g gas 1 . 4 ; s e t R gas 2 8 7 . 0 ; s e t M inf 5 . 0 ; # S p e c i f i e d Mach number # Select a st at ic pressure s e t p i n f 1 0 0 . 0 e3 ; set T inf 1 0 0 . 0 ; # and a t e m p e r a t u r e # We need t o d e t e r m i n e v e l o c i t y . s e t a i n f [ e x p r s q r t ( $ T i n f ∗ $R gas ∗ $ g g a s ) ] ; # sound s p e e d # velocity s e t u i n f [ expr $M inf ∗ $ a i n f ] ; # Also , handy t o know dynamic p r e s s u r e f o r n o n d i m e n s i o n a l i z a t i o n # o f t h e p r e s s u r e s and drag f o r c e s . s e t q i n f [ expr 0 . 5 ∗ $g gas ∗ $ p i n f ∗ $M inf ∗ $M inf ] p u t s ” Free−s t r e a m v e l o c i t y , u i n f = $ u i n f ” puts ” static pressure , p inf = $p inf ” puts ” dynamic p r e s s u r e , q i n f = $ q i n f ” # For t r a n s i e n t s i m u l a t i o n , we s t a r t with a low p r e s s u r e . set p init 1000.0; set T init 100.0; BEGIN GEOMETRY # F i r s t , s p e c i f y s u r f a c e o f c y l i n d e r and wedge node a 0 . 0 0 . 0 ; # Centre of curvature f o r nose s e t xb [ e x p r − 1 . 0 ∗ $Rn ] node b $xb 0 . 0 set xc [ e x p r − 1 . 0 ∗ $Rn ∗ s i n ( $ a l p h a ) ] set yc [ e x p r $Rn ∗ c o s ( $ a l p h a ) ] node c $xc $yc a r c bc b c a # Down−s t r e a m end o f wedge s e t yd [ e x p r $yc + ( $xEnd − $xc ) ∗ tan ( $ a l p h a ) ] node d $xEnd $yd p u t s ” h e i g h t a t end o f p l a t e yd= $yd ” l i n e cd c d # Outer−edge o f f l o w domain has t o c o n t a i n t h e s h o c k l a y e r # Allow s u f f i c i e n t f o r s h o c k stand −o f f a t t h e s t a g n a t i o n l i n e . s e t R2 [ e x p r $Rn + $ d e l t a ] set xe [ e x p r − 1 . 0 ∗ $R2 ] node e $xe 0 . 0 # The s h o c k a n g l e f o r a 1 0 d e g r e e ramp with s h a r p l e a d i n g edge # i s 2 0 d e g r e e s ( r e a d from NACA 1 1 3 5 , c h a r t 2 ) , # however , t h e b l u n t n o s e d i s p l a c e s t h e s h o c k a l o n g way out # s o we a l l o w some more s p a c e . # We need t o s e t t h e boundary h i g h enough t o a v o i d t h e s h o c k s e t R3 [ e x p r $Rn + 2 . 0 ∗ $ d e l t a ] set x f [ e x p r − 1 . 0 ∗ $R3 ∗ s i n ( $ a l p h a ) ] set y f [ e x p r $R3 ∗ c o s ( $ a l p h a ) ] node f $xf $yf # Now , put i n i n t e r m e d i a t e c o n t r o l p o i n t s s o t h a t we can u s e # c u b i c B e z i e r c u r v e f o r t h e i n f l o w boundary around t h e n o s e # and a s t r a i g h t l i n e downstream o f p o i n t f . node e1 $xe $ d e l t a set alpha2 [ t o r a d i a n s 4 0 . 0 ] set xf1 [ expr $xf − $delta ∗ cos ( $alpha2 ) ] set yf1 [ expr $yf − $delta ∗ s i n ( $alpha2 ) ] node f 1 $ x f 1 $ y f 1

21

bezier set yg node g li n e fg

e f e e1 f 1 f [ e x p r $ y f + ( $xEnd − $ x f ) ∗ tan ( $ a l p h a 2 ) ] $xEnd $yg f g

# D e f i n e s t r a i g h t −l i n e s e g m e n t s between s u r f a c e # and o u t e r boundary . l i n e eb e b line fc f c l i n e dg d g # Assemble t h e c u r v e s e g m e n t s i n t o p o l y l i n e s # t h a t w i l l become t h e b l o c k b o u n d a r i e s . p o l y l i n e n0w1 1 + f c p o l y l i n e s0 1 + eb p o l y l i n e e0 1 + bc p o l y l i n e w0 1 + ef p o l y l i n e n1 1 + fg p o l y l i n e e1 1 + dg p o l y l i n e s1 1 + cd END GEOMETRY BEGIN FLOW g a s t y p e PERF AIR 14 $p inf gas state inflow gas state in it ia l $p init s e t nnx0 4 0 s e t beta0 1 . 2 s e t nny0 4 0 d i s c r e t i s e n0w1 d i s c r e t i s e s0 d i s c r e t i s e e0 d i s c r e t i s e w0

$nnx0 $nnx0 $nny0 $nny0

0 0 0 0

1 1 0 0

$u inf 0 . 0 $T inf 0.0 0.0 $T init

1.0 1.0

$beta0 $beta0 0.0 0.0

s e t nnx1 1 0 0 s e t nny1 $nnx0 ; # c o n n e c t i n g a n o r t h edge t o a west edge s e t beta1 1 . 2 d i s c r e t i s e n1 $nnx1 1 0 $ b e t a 1 d i s c r e t i s e s1 $nnx1 1 0 $ b e t a 1 d i s c r e t i s e e1 $nny1 0 0 0 . 0 b o u n d a r y s p e c w0 s u p i n i n f l o w b o u n d a r y s p e c n1 s u p i n i n f l o w b o u n d a r y s p e c e1 e x t r a p o l a t e o u t b l o c k blk −0 + n0w1 + e0 b l o c k blk −1 + n1 + e1

+ s0 + s1

+ w0 − n0w1

c o n n e c t b l o c k s blk −0 n o r t h blk −1 west f i l l b l o c k blk −0 i n i t i a l f i l l b l o c k blk −1 i n i t i a l END FLOW BEGIN CONTROL s e t t i t l e s t r i n g ” Blunt Wedge Rn=$Rn , ” append t i t l e s t r i n g ” q i n f =[ f o r m a t ” % 1 2 . 3 e ” $ q i n f ] , ” append t i t l e s t r i n g ” yd=[ f o r m a t ” % 1 0 . 5 f ” $yd ] ” title $title string case id 0 f l u x c a l c adaptive s e t t f i n a l [ e x p r 4 0 . 0 ∗ $Rn / $ u i n f ] p u t s ” F i n a l time = $ t f i n a l ” max time $t final max step 500000 t i m e s t e p 1 . 0 e−8 dt plot $t final h i s t o r y c e l l blk −0 $nnx0 1 d t h i s t o r y [ expr $ t f i n a l / 1 0 0 . 0 ] END CONTROL

22

b e z i e r f i l e bw . bez p a r a m f i l e bw . p m p o s t f i l e bw . mpost mpost scales 1 . 5 1 . 5 build quit

4.2

Shell scripts

# bw prep . sh # s c r i p t i t . t c l − f bw . s i t −do−mpost > bw . s c r i p t i t . l o g mpost bw . mpost mb prep . e x e − f bw mb post . e x e − f p bw . p − f g bw . g − f s bw . s 0 − f o bw 0 − g e n e r i c XMIN=−20.0 e−3 XMAX=100.0 e−3 YMIN=0.0 YMAX=100.0 e−3 TIC=20.0 e−3 mb cont . e x e − f i bw 0 . gen − f o bw 0 mesh . ps − ps − v a r 6 − mesh \ −x r a n g e $XMIN $XMAX $TIC − y r a n g e $YMIN $YMAX $TIC

# bw run . sh # time mb cns . e x e − f bw − f o r c e 0 1 − f o r c e 1 2 mv mb cns . l o g bw . mb cns . l o g echo ” Done”

# bw post . sh TIMES=”399” XMIN=−20.0 e−3 XMAX=100.0 e−3 YMIN=0.0 YMAX=100.0 e−3 TIC=20.0 e−3 f o r TME i n $TIMES do mb post . e x e − f p bw . p − f g bw . g − f s bw . s − f o bw $TME − t $TME. 0 e −6 − g e n e r i c mb cont . e x e −x r a n g e mb cont . e x e −l e v e l s −x r a n g e mb cont . e x e −l e v e l s −x r a n g e

− f i bw $TME . gen − f o bw ”$TME” p . ps − ps − v a r 6 − c o l o u r \ $XMIN $XMAX $TIC − y r a n g e $YMIN $YMAX $TIC − f i bw $TME . gen − f o bw ”$TME” M . ps − ps − v a r 7 − c o l o u r \ 0.0 5.0.0 0.2 \ $XMIN $XMAX $TIC − y r a n g e $YMIN $YMAX $TIC − f i bw $TME . gen − f o bw ”$TME” s o n i c . ps − ps − v a r 7 − c o l o u r \ 0.0 1.0 0.2 \ $XMIN $XMAX $TIC − y r a n g e $YMIN $YMAX $TIC

done

23

# b w f o r c e . sh # P l o t t h e s u r f a c e p r e s s u r e on t h e wedge TME=399 NX=40 mb prof . e x e − f p bw . p − f g bw . g − f s bw . s − f o b w s u r f a c e . dat \ −t $TME. 0 e −6 − x l i n e 0 $NX − y l i n e 1 1 awk − f s u r f a c e p r e s s u r e . awk b w s u r f a c e . dat > b w s u r f a c e p c o e f f . dat g n u p l o t