NASA Contractor Report 181790 

Noise Produced by Turbulent Flow 

Into a Rotor: 

Users Manual for Noise Calculation 

R. K. Amiet, C. G. Egolf, and J. C. Sirnonich 

United Technologies Research Center 

East Hartford, CT 

Contract NAS1-17763 

June 1989 

National Aeronautics and 

Space Administration 

Langley Research Center 

Hampton, Virginia 23665-5225 

NA9-29152 

(NASA-CR-181790) NOISt PROOUCEC) RY 

TURBlJLEYT FLOW INTO A ROTOR: USER5 MAHUAL 

fOR NOISE CALCULATION Final Report (United Technologies Research Center) 38 p CSCL ZOA uncl as 

G3/71 0225002 
Table of Contents 

Noise Calculation Programs 

Subroutine NOPMAN ............................................... 1 

Subroutine NOPINP ............................................... 4 

Subroutine NOPPX ................................................ 7 

Subroutine NOPCRS .............................................. 10 

Subroutine NOPKVC .............................................. 11 

Subroutine NOPFNL .............................................. 13 

Subroutine NOPLFT .............................................. 14 

Subroutine NOPNI ............................................... 16 

Appendix A: Computer Program for Noise Produced by Rotor.. Turbulence 

Interaction ............................................ 19 

Appendix B: Subroutine for Calculating the Main Problem Parameters 

and Integrating Over Azimuth ........................... 23 

Subroutine for Summing Over the Turbulence Spectrum for 

Blade-to-Blade Correlation ............................. 26 

Appendix C: 

Appendix D: Subroutine for Calculating Airfoil Response I’unction for 

Leading Edge Noise ..................................... 30 

Appendix E: Subroutine for Calculating the Effect of Contraction on 

the Turbulence Spectrum ................................ 31 

Appendix F: Subroutine for Calculating the Cross-Product of Two 

Vectors ................................................ 33 

Appendix G: Subroutine for Calculating the Fresnel Integrals ....... 34 

Appendix H: Sample Calculation ..................................... 35 

i 
*** 

t 

t 

* 

t 

* 

t 

f 

t 

t 

t 

* 

t 

t 

* 

t 

t 

t 

t 

* 

t 

t 

* 

* 

t 

t 

t 

t 

t 

t 

t 

t 

t 

t 

t 

t 

t 

t 

* 

t 

* 

* 

t 

t 

* 

* 

t 

* 

t 

t 

1 

ROUTINE - NOPHAN 

PURPOSE - The sain program which begins the calculation. 

It also does the summation over radius. 

AUTHOR - Roy K. Aniet 

INPUT 

COMMON BLOCKS 

COMMON IROTCAI 

Pame Type Description 

CL RD ChordJturbulence length scale 

B% RD Number of blades 

COC RD Sound speed 1 chord 

RC ED Blade radiualchord 

R1C RD Far-field distance J chord 

RPS RD Rotational speed revolutions 1 sec 

DXDZ RD Deformation tensor 

ALPHA ED Tip path plane angle of attack 

A0 RD Speed of sound 

T14IN ED Minumum observer polar angle in degrees 

TMXX RD Maximum observer polor angle in degrees 

n RD Diameter of rotor 

DELT ED Desired increment in polar angle between 

successive calculations 

Pnxx RD Minimum observer azimuthal angle in degrees 

(phi=O in direction of mean wind) 

PMAX RD Maximum observer azimuthal angle In degrees 

DELP RD Desired Increment in azimuthal angle between 

successive calculations 

214 RD Axial Mach number 

FN RD Flow Mach numbr in rotor plane 

uu RD Rms turbulence velocity I axial velocity 

GWS ED Geostrophic wind speed (m/s) 

WUINF RD Vertical component of the rms turbulence 

........................................................... 

(theta=O along rotor axis in thrust direction) 

normalized by the stream velocity 

COMMON lNOPCBl 

Hame Type Description 

J I Number of azimuthal integration points 

HF I First harmonic 

HL I Last harmonic 

HUM I Number of points 

IH I Homogenous or nonhomogenous indicator 

0 homogenous case 

1 nonhomogenous 

................................................................... 

NR I Number of radial points 

1 
t 

t 

t 

t 

t 

t 

f 

* 

t 

t 

t 

t 

* 

t 

t 

t 

t 

t 

f 

t 

t 

t 

t 

f 

f 

t 

t 

t 

t 

* 

t 

t 

t 

t 

t 

* 

t 

t 

t 

t 

t 

t 

t 

t 

t 

t 

t 

t 

t 

t 

t 

t 

t 

t 

LOCAL VARIABLES 

Name Type Description 

............................................................. 

DH 

DR 

H 

HOLD 

IARRAY 

I2 

IT 

IP 

ISTXT 

JJ 

K 

LPOUT 

NOPIN2 

NP 

IR 

IT 

NWRDS 

P2 

P22 

PH 

RARBAY 

RCP 

ECP 1. 

SUM 

TH 

OUTPUT 

Name 

RD 

RD 

RD 

RD 

I 

I 

I 

I 

I 

I 

I 

I 

RD 

I 

I 

I 

RD 

RD 

RD 

RD 

RD 

RD 

RD 

RD 

T 

1 

Type 

Fractional harmonic spacing bet-Jeen frequency 

calculat ions 

Increment between radial points 

Frequency of aound in multiples of blade passage 

harmonics 

Last harmonic calculated 

Array containing integer vords from a data 

member record 

Do loop index 

Do loop index 

Number of P values 

Status of calls to data member routines 

Do loop index 

Do loop index 

Logical unit to vrlte 

Array containing data member name of HOP 

input file 

Do loop index 

Number of radial points 

Do loop index 

Humber of vords read in from data mamber record 

Pressure averaged around azimuth 

Previous P2 value 

Angle of flight Mach number 

Array containing the real vords from a data memb 

record 

Radial point 

Previous radial point 

Sum of pressure squared contributions to noise 

Polar angle of observer 

Description 

............................................................... 

H RD Frequency of sound in aultiplea of blade 

passage harmonic 

F1 RD Frequency of sound in Hz 

PSD RD Spectrum level in dB (per unit Hz) 

NP 1 RD Summation counter 

IM 1 RD Summation counter 

N RD Summation counter 

FUNCTIONS 

1. To call subprogram to read in input values 

2. To call subprogram giving spectrum integrated over azimuth 

3. To perform integration of spectrum over radlus 

SUB3ROGBAMS CALLED 

LlNCLOS MIGETR MMOPRD NOPIlP PX XFETCH XSTORE 

2 
* 

* 

t 

t 

t * 

t 

t 

t 

t 

t 

t 

t 

* 

* 

t 

t 

* 

* 

t 

* 

t 

t 

t 

t 

* 

* 

t 

t 

t 

t*t 

CALLING PRQGRA!.IS 

ANOP EXEC 

ERRORS 

NOH-FATAL 

FATAL 

1. Error opening data member for input 

none 

ENTRY 

Print the input quantities for checking 

Write tho values of the deformation tensor 

If beginning harmonic is negative, return to matrix input 

If harmonic spacing is negative, return to previous input 

If number of steps is negative, return to first input 

Calculate frequency step size 

Print headings for the eventual output 

Initialize the frequency variable 

Input of zero for L2 defaults to a single frequency step 

DO 60 K=l,L2 

Increment frequency 

Initialize radius 

Initialize 3um 

DO 22 12=1,19 

Increment blade radlal position 

Calculate tip Mach number at this radius 

Call subroutine giving azimuthally averaged spectrum at this 

radiua 

Halve first value in trapezoidal integral over radius 

Add to integrated spectrum 

END DO 

Calculate PSD 

Calculate haraonic 

Write results 

END DO 

EXIT 

3 
ROUTINE - MOPINP 

~ 

* 

* 

t 

* 

* 

* 

* 

* 

* 

* 

t 

t 

I t 

* 

* * 

* 

t * 

* 

* 

t 

t 

t 

* 

t 

t 

* 

t * * 

t 

* 

t 

* 

* 

t 

I * 

* 

* I * I 

* 

~ * 

PURPOSE - To read in input variables and calculate the inean 

and turbulence properties for an atmospheric boundary 

layer 

AUTHOR - J. C. Simonich 

INPUT 

USER PARAMETERS 

Name 

A0 

D 

C 

BH 

FF 

EPS 

HF 

N UM 

TMIH 

.......... 

nL 

TMAX 

DELT 

PMIN 

PHAX 

DELP 

TY Pe 

........ 

RD 

RD 

RD 

RE) 

RD 

RD 

I 

I 

I 

RD 

RD 

RD 

RD 

RD 

ED 

Description 

Speed of sound (m/s) 

Diameter of rotor (m) 

Blade chord (m) 

Number of blades 

Distance from rotor to observer (m) 

Rotational speed in revlsec 

First harmonic 

Last harmonic 

Number of frequency points 

Minisum observer polar angle in degrees 

(theta=O along rotor axis in thrust direction) 

Maximum observer polar angle in degrees 

Desired increment in polar angle between success 

calculations 

Minisum observer azimuthal angle in degrees 

(phi=O in direction of mean vind) 

Maximum observer azimuthal angle in degrees 

Desired increment in azimuthal angle betveen 

successive calculations 

............................................ 

DATA MEMBER ROTNOP(ABLOT1) - This is the output from the FM ABL 

GUS, LWX, WUIHF 

Name Type Description 

GWS RD Geostrophic vind speed (m/s) 

LWX RD Correlation length scale 

WUINF RD Vertical component of the rms turbulence 

................................................................... 

normalized by the velocity 

DATA MEMBER ROTNOP(ROTOT1) - This i3 the output from the FM ROT 

TITLE 

118, J, IH 

R2, ALPHA, R 

STRENG 

If IH=O (homogenous case) read in following: 

T, iT, V, W, X, Y, 2 

XSLL, YSLL, ZSLL, XSLZ, YSLZ, ZSLZ 

USAV, VSAV, USAV 
t 

t 

t 

t 

* 

t 

f 

t 

* 

t 

t 

t 

* 

t 

t 

t 

* 

t 

* 

t 

t 

* 

f 

t 

t 

f 

* 

* 

t 

t 

* 

t 

t 

t 

* 

t 

t 

t 

t 

t 

t 

* 

t 

t 

f 

f 

f 

t 

t 

* 

f 

* 

* 

XA11, xx21, PA31 

XA12, XA22, XA32 

XA13, XA23, XA33 

Name , Type Description 

TITLE ’ A Titls of case run 

NR I lumber of radial points 

J I lumber of azimuthal points 

In I Indicator of homogenous or nonhomogenous case 

R2 RD Maximum radius 

ALPHA RD Rotor tip path plane angle of attack 

R RD Radius of rotor 

STRENG RD Vortex circulation strength 

For homogenous case (IH=O) the folloving input is read, othervise 

the rest of the data member is read in the subroutine NOPPX 

................................. b....~.~.~~~..~....~.~.~..~~9 

T 

U 

V 

W 

X 

Y 

2 

XSLl 

YSLl 

ZSLi 

XSL2 

YSLZ 

ZSL2 

USAV 

vsxv 

ZSAV 

xx12,. . . 

RD 

ED 

RD 

ED 

RD 

RD 

RD 

RD 

RD 

RD 

RD 

RD 

RD 

AD 

ED 

RD 

RD 

Streamline drift time 

U velocity component in radiilsec 

V velocity component in radiilsec 

U velocity component in radiilsec 

Dovnstream X coordinate of streamline in 

rotor plane 

Dovnstream Y coordinate of streamline in 

rotor plane 

Dovnstrean 2 coordinate of streamline in 

rotor plane 

Dovnstream X coordinate of streamline in 

standard coordinate system 

Dovnstream Y coordinate of streamline in 

standard coordinate system 

Dovnstream 2 coordinate of streamline in 

standard coordinate system 

Upstream X coordinate of streamline in 

standard coordinate system 

Upstream Y coordinate of streamline in 

standard coordinate system 

Upstream 2 coordinate of streamline in 

standard coordinate system 

Uptsream U velocity component in radiilsec 

Upstream V velocity component in radiilsec 

Upstream 2 velocity component in radiilsec 

Deformation tensor 

LOCAL VARIABLES 

Name Type Description 

IARRXY I Array to store a data member record 

IHDE I Array containing the length of the largest 

I1 P Do loop index 

I STAT I Status of calls to data member routines 

JJ I Do loop index 

LPOUT I Logical unit to vrite 

.......................................................... 

record and the number of records written 

T ITYPE A Indicator of user parameter type 

5 
* 

t 

f 

t 

t 

* 

t 

* 

t 

t 

* 

t 

t 

t 

* 

* 

t 

* 

t 

* 

t 

t 

t 

* 

t * 

t 

t 

* 

* 

t 

t 

t 

* 

* 

t 

* 

t 

t 

t 

t 

t 

* 

t 

* 

t 

t 

t 

t 

**f 

T MEX I S T A Indicator of existance of data mesber 

ElEL I Indicator of number of elements in usar 

parameter 

NOPIal RD Array containing data sember name of ABL output 

file 

NOPIN2 RD Array containing data sember name of ROT output 

file 

NWRDS I Number of words read in from data member record 

PI R9 3.14159. .. RAREIY RD Storage area for real variables being read in fr 

data member record 

OUTPUT 

COMMON BLOCKS 

NOPCA - sae IOPMAN 

NOPCB - see NOPMAN 

FUN CT I 0 IS 

1. To input variables for NOP. 

2. To input variables from FM ROT. 

3. For homogenous case read in defornation tensor and calculate axial 

mach number, flov mach nuaber and rms tarbulence velocity. 

SUBPROGRAMS CALLED 

HMCLOS, NNGETR, Ml4OPRD, MMVUM, HEUTRAL, STABLE, UNSTAB, XFETCH, 

XSTORE, XGETP , XASKP 

CALLING SUBPROGRAM 

NOPMAN 

ERRORS 

NO# - FATAL 

1. Error finding user parameter in table 

2. Error in opening data member 

3. Error reading in record from data member 

FATAL 

none 

ENTRY 

Input values from data member ROTNOP(NOPIN1) 

If L ) 0 call STABLE 

If L (,O call UNSTAB 

If L = 0 call NEUTRL 

Calculate vertical component of the rms turbulence normallzed by the 

stream velocity 

Input data member ROTNOP(ROTOT1) 

If IH=O read in deformation tensor 

EXIT 

6 
t 

t 

t 

* 

* 

t * 

* 

t 

* 

t 

* 

* 

* 

* 

t 

* 

* 

t 

f 

t 

t 

t 

* 

f 

* 

t 

t 

t 

t 

f 

t * 

t 

t 

t 

t 

t 

t 

t 

t 

f 

t 

* 

t 

t 

t 

* 

t 

t 

t 

t 

t 

t 

ROUTIIE - MOPPX 

PURPOSE - To take the cross product of column vectors N2 and N3 

of E and place the result in column N1. 

AUTHOR - Roy iZ. Amiet 

INPUT 

ARGUMENTS 

Name Type Description 

.................................................................. 

NOPINZ RD Array containing data member name 

F1 RD Observer frequency 

RCP RD Radial point 

TH RD Polar angle of observer 

PH RD Angle of flight Mach number 

COMMON BLOCK 

NOPCX 

Name 

CL 

BN 

coc 

R?C 

RPS 

DXDZ 

211 

FH 

uu 

............ Type Description 

RO C3ordjturbulance intagral length scale 

RD Blade number 

RD Sound speedlchord 

RD Far-field diatance/chord 

RD Revolutions/second 

RD Deformation tensor 

RD Axial Mach number 

RD 

AD RHS turbulence velocitylaxial velocity 

.................................................... 

Flight Mach number in rotor plane 

COMMON BLOCK 

UOPCB 

Name Type Description 

............................................................... .J I Number of azimuthal integration points 

IH I Specifies homogenous or 

inhomogenous calculation 

OUTPUT 

ARGUMENT 

Name Type Description 

P2 RD Pressure averaged around azimuth 

tlP 1 I Countsr for upvard summation in NTITRB 

NM 1 I Counter for dovcvard summation in HTITRB 

............................................................... 

subroutine 

7 
t 

t 

* 

t 

t 

t 

t 

* 

* 

t 

t 

t 

* 

t 

t 

* 

t 

t 

t 

t 

f 

t 

* 

t 

t 

* 

t 

t 

t 

t 

t 

t 

t 

t 

* 

* 

t 

t 

* 

* 

t 

* 

* 

t 

t 

* 

* 

t 

* 

* 

* 

f 

t 

, 

aubroutine 

subroutine 

N I Countsr for rezcaling of CVT in YTITRB 

LOCAL VARIABLES 

Name Type Description 

................................................................ 

T M 

AJ 

DEL 

C 

S 

C4 

Sf 

OM 

BFZ2 

SQ 1 

RER 

SUM 

G 

I 

C2 

S2 

cs 

s5 

ALP 

Cl 

s1 

EM 

c3 

X 

Y 

Z 

SG 

FP 

XK 

YK 

T 

T1 

xx 

YY 

22 

T2 

CVT 

ZKO 

AD 

RD 

RD 

RD 

RD 

RD 

RD 

RD 

RD 

RD 

ED 

RD 

RD 

RD 

I 

RD 

ED 

RD 

RD 

RD 

RD 

RD 

RD 

RD 

ED 

RD 

RD 

RD 

ED 

RD 

RD 

RD 

RD 

RD 

ED 

RD 

RD 

RD 

RD 

RD 

SUB P R 0 G A RAM S CALL ED 

NOPNI XOPLFT 

CALLING PROGRAMS 

NOPMAX 

Local tip Mach number 

Float J 

Azimuthal step size 

Cosine theta 

Sine theta 

Cosine psi 

Sine psi 

Expression in analysis 

Expression in analysis 

Expression in analysis 

Retarded radius/actual radius 

Sum of azimuthal spectral contributions 

~zimuthal angle gamma 

DO loop counter 

Cosine gamma 

Sine gamma 

Cosine (gamma + psi) 

Sine (gamma t psi) 

ingle alpha 

Cosine alpha 

Sine alpha 

Mach number component along chord 

Cosine of angle phi in analysis 

Observer coordinate 

Observer coordinate 

Observer coordinate 

Modified radius 

Frequency measured on blade 

Value of x vavenumber 

Value of y vavenumber 

Blade passage time 

The betveen eddy chops 

Specific value of x in analysis 

Specific value of 'J in analysis 

Specific value of z in analysis 

Time T1 plus propagation time difference 

Step size for summation over vavenumber 

Initial value of vavenumber in summation 

Contribution to sum of particular aziinuthal 

station 

8 
t 

t 

* 

t 

* 

t 

* 

t * 

* 

t 

t 

t 

t 

t 

t 

t 

t 

t 

t 

* 

t 

* 

* 

* 

* 

t 

t 

t 

t 

t 

t 

t 

t 

* 

t 

t 

t 

t 

t 

NON -FATAL 

1. Error reading in data mamber record 

FATAL 

none 

ENTRY 

Calculate TN, the local tip Hach number 

Float J 

Calculate DEL, the azimuthal step size 

C, S and C4, SA, are the cosine and sine of theta and psi 

OM Is a factor appearing in the analysis 

BFZ2 is a ?randtl-Glauert factor 

SQl is a factor appearing in the analysis 

RER is the observer retarded distance/actual distance 

Initialize XP1,. ... G 

DO 50 I=l,J 

Increment gamma 

Find cosine and sine of gamma and gamma + psi respectively 

Calculate ALP, the angle of the rotor blade vrt the rotor plane 

Calculate C1 and S1, the cosine and sine of alpha 

Calculate RM, the Mach number of rotor segment vrt to fluid 

Calculate C3, the cosine of an angle in the analysis 

Calculate observer coordinates in rotor fixed coordinates 

Calculate sigma, a sodified observer distance 

Calculate XK and YK, the x and y vavenumbers 

Calculate TI the time between blade passes 

Calculate T1, the time betveen eddy intersection3 

Calculate XX, YY, ZZ which are X, Y, Z values 

Calculate CVT which is 2*Pi/eddy passing distance 2 

Calculate ZfO, the Initial radian frequency in the summation 

Call turbulence summation subroutine 

Call airfoil response subroutine 

Incrsment counters 

Contribution to spectrum from azimuthal integration 

Summation over azimuthal spectrum 

respectlvely 

Integration over azimuthal angle gamma 

. Calculate FP, the frequancy on the blade 

End DO 

Multiply by remaining factors 

EXIT 

9 
t** 

* 

t 

t 

t 

t 

t 

t 

t 

t 

t 

t 

t 

t 

t 

f 

t 

t 

t 

* 

* 

t 

t 

t 

t 

t 

t 

f 

t 

t 

t 

t 

t 

t 

t 

f 

t 

t 

t 

t 

t 

t 

t 

t 

t 

t 

t 

t 

t 

t 

*** 

ROUTINE - MOPCBS 

PURPOSE - To take the cross product of column vectors N2 and N3 of E 

and place the result in column Hi. 

AUTHOR - Roy R. Xmiet 

IBPU? 

ARGUMENT 

Name Type Description 

E RD 3x3 rotation matrix. 

H 1 I Column into which the cross product is to be 

.............................................................. 

placad. 

OUTPUT 

ARG U !f E N T 

lame Type Description 

E ilD Unit vectors aiong coordinates 

............................................................. 

LOCAL VARIXSLES 

Name Type Description 

N2, N3 I Column vectors for which the cross product is found. 

E3 ?M ED lornalizea the output to a magnitude of 1. 

................................................................. 

SUBPROGARAMS CALLED 

Hone 

CALLING PROGRAMS 

NOPKVC 

ERRORS 

None 

ENTRY 

Calculate N2 and N3 so that N1, N2, N3 are in cyclic permutation. 

Calculate croaa product of column vectors l2 and N3; place result 

column N1. 

Find the magnitude of vector N1. 

Horsalize vector N1. 

EXIT 

10 
h 

f 

f 

f 

f 

f 

t 

f * 

t 

* 

t 

t 

t 

t 

* 

* 

t 

* 

t 

f 

t 

t 

R 

* 

t 

t 

t 

t 

t 

t 

t 

t 

t 

t 

t 

t 

t 

t 

t 

t 

i 

* 

i 

t 

* 

t 

t 

t 

ROUTINE - HOPRVC 

PURPOSE - To calculate effect of rapid distortion on turbulence 

spectrum . 

AUTHOR - Roy K. Amiet 

INPUT 

ARGUMEHT 

llame Type Description 

.............................................................. 

DK RD Wavevector 

DXDZ RD Deformation tensor 

VH RD Unit vector normal to blade 

OUTPUT 

XBGUElEHT 

Name Type Description 

.............................................................. 

ED RD Three coordinats vectors downstream of 

E3 RD Three coordinats vectors upstream of 

DKN RD Magnitude of downstream vavevector. 

DKUK ED Ratio of downstream to upstream vavevectors. 

UQDQ RD Ratio of upstream to dovnstream turbulent 

contract ion. 

contraction. 

velocities. 

LOCAL VARIABLES 

Name Type Description 

F(I) RD For temporary storage of EU(1,J). 

AF 1 BD Magnitude of wavevector, to be used for 

normalization. 

................................................................ 

FUNCTIOHS 

1. To calculate the effect of rapid distortion on the wavevector 

and Fourier component amplitude 

SUBPROGARAMS CALLED 

NOPCRS 

CXLLIIG PROGRAMS 

YOPII 

ERRORS 

NONE 

n. m ril 3Y 

Find magnitude of vavevector. 
f 

f 

t 

t 

* 

* 

f 

t 

t 

* 

t 

t 

t 

f 

* 

* 

* 

* 

f 

* 

f 

*it 

DO I = 1,3 

Sat thlrd column of ED equal to vavevector direction. 

Set second column of EL) equal to blade nornai. 

E!lD DO 

Find cross product of first two vectors 

Find cross product of previocs result vlth third vector of ED. 

DO I=1,3 

DO J=1,3 

Initialize EU 

DO K=1,3 

EX3 DO 

Take product of ED with deformation natrix to find E3. 

EID DO 

END DO 

Store third column of EU in F. 

Find the magnitude of first column of EL'. 

loraalize tha first column of EIJ 

Find the cross-product of columns 1 and 2 of EU. 

Find cross-product of previious result with third vector of EU. 

Calculate ratron of dovnstream to upstream wavevector magnitude. 

Calculate ratio of upstream to downstream turbulent velocities. 

BXIT 

12 
* 

* 

* 

* 

* 

I 

t 

A 

t 

d 

4 

f 

* 

* 

B 

* 

Q 

a 

n 

rt 

It 

t 

d 

* 

9 

1' 

* 

* 

* 

* * 

t 

t 

t 

t 

It 

t 

* 

t 

t 

* 

* 

t 

t 

* 

t 

t** 

HOUTIlIE - HOPFNL 

PURPOSE - To calculata the Fresnel integrals C and S. Ref. Abromovitz 

and Stegun, p 302. 

AUTHOR - Roy K. Amiet 

INTUT 

XRGUNENT 

Name Type Description 

............................................................. 

X RD Fresnel integral argusent 

OU'=?UT 

ARGUMENT 

Name Type Description 

............................................................... 

E C Fresnel integrals in form C - IS 

LOCAL VARIABLES 

Nase Type Description 

(2 c Intar3ed:ate variabies from Abromovitz and Stegun. ................................................................ 

- 

H C To rotate G in complsx plane. 

FU B CT I 0 IS 

1. To calculate a value for the Fresnel integrals 

SUBPROGARAMS CALLED 

None 

CALLING PROGRAMS 

'IOPLFT 

EBTRY 

Calculate G from Abronovitz and Stegun equations p 302. 

Calculate H. 

Calculata E = C - IS following Abromovitz and Stegun. 

EXIT 

13 
t 

t 

I * 

f 

* 

t 

t 

ROUTIBE - tJOPLFT 

PURPOSE - To calculate the airfoil gust response function. 

AUTSCR - Roy P.. Amiet 

ItIPUT 

XRGUMEMT 

Name Type Decsription 

x RD Chordvise distance 

SG ED Far-field distance modified by Prandtl-Glauert 

R!4 RD Chordwise Mach number component 

XX BD Chordvise vavenum5er 

YR RD Spanvise wavenumber 

................................................................ 

factor 

OUTPUT 

ARGUMENT 

Name Type Description 

GL2 RD Square of effective lift, including 

.................................................................. 

noncompactness 

LOCAL VARIABLES 

Name Type Description 

B2 RD Column vectors for which the cross product is 

UM RD Normalizea the output to a magnitude of 1. 

RMI ED Similarity Mach number for skeved gusts 

UMI RD Specifies use of lov or high frequency solution. 

T2 RD Intermediate dummy variable. 

c:rl RD Magnitude of E. 

E RD Complex representation of Fresnel integrals; 

................................................................... 

found. 

c - is 

FUNCTZOIIS 

1. To calculata the airfoil response function for given values of 

the wavenumbers 

SUBPRCGARAMS CALLED 

NOPFNL 

CALLING PROGRAMS 

HOPPX 

ENTRY 

Calcnlata Prandtl-Glauert factor 32. 

UM is proportional to airfoil chord/acaust:c vavo:eqtL. 

RMI 13 Graham'~ siailaritjr Hach number for skeved gust;. 

14 
t 

* 

* 

t 

* 

t 

t * 

t 

.p 

tjtt 

UMI used in fo::uving link for branching tsat. 

ea38 af (UMI -.75) NES, ZEERC), 305 

IEG: Use sodified Saars function 

ZERO: Use nodified Sears function 

POS: Use high reduced frequency solution 

Call Fresnel Integrals 

Find absolute value 

Calculate square of effective lift 

End case. 

EXIT 

15 
t 

I * 

* 

* * 

* 

t 

t 

t 

* 

t 

* 

* 

* 

I 

I * 

* 

t 

* 

* 

* 

* 

f 

t 

* 

* 

* 

* 

* 

* 

* 

t 

* 

* 

* 

* 

t 

* 

* 

* 

t 

t 

I 

I 

I 

~ 

ROUTINE - :lC?NI 

* 

* 

PURPOSE - To take tha cro33 product of column vectors 12 and !I3 sf 

E and place the result in column ill 

AUTHOR - Roy R. Amiet 

INPUT 

A R G U M E ?I T 

Name Type Description 

.............................................................. 

CL RD Chord/turbulence integral length scale 

CVT RD Step size for sulnmation Over vavenumber 

XK RD X vavenumber component 

YK RD Y vavenumber component 

ZK RD Z vavenumber component 

ALP RD Angle alpha 

G RD Gamma 

DXDZ RD Deformation tensor; 9 components 

OUTPUT 

A R G WEN T 

Name Type Description 

SYV RD Pressure averaged around azimuth 

NP I Counter for upward summation in NTITRB 

............................................................. 

subroutine 

subroutine 

subroutine 

v IM A Counter for dovnvard summation in NTITRB 

I I Counter for rescaling of CVT in HTITRB 

LOCAL VARIABLES 

Name Type Description 

EK RD Constant in Karaan spectrum 

DK RD Wavevector; 3 components 

C1 RD Cosine alpha 

SI. RD Sine alpha 

C2 RD Cosine gamaa 

S2 RD Sine gamma 

YH RD Normal to blade 

BOT AD Expression in analysis 

DXDZRl RD Intermediate deformation matrix 

DXDZA ED Final deformation matrix after rotation 

CVT 1 AD Dummy for CVT 

ED RD Downstream unit vectors parallel to axes 

EU RD Upstream unit vectors parallel to axes 

DKM R3 Wavevector amplitude downstream at rotor face 

DKUX 23 Ratio of dovnatreaiii to xJstreaa *Jravevector 

............................................................ 

amplitudes 

16 
* 

* 

t 

* 

P 

i' 

f 

.e 

f 

* 

4 

n 

D 

3 

1' 

* 

t 

t 

% 

R 

f 

t 

f 

f 

* 

t 

* 

t 

t 

t 

t 

* 

t 

t 

* 

t 

* 

* 

* 

* 

t 

* 

* 

UQDQ BD 

XKUO RD 

YKUO ED 

XY 0 RD 

DK!4 1 ED 

DKUKl . RD 

UQDQl RD 

XKU RD 

YK'J RD 

ZKU RD 

XY BD 

Dfl BD 

DEL1 ED 

Al RD 

ZKl ED 

ZK2 RD 

Ratlo of upntream to dovnstream velocity 

amplltudss 

Upstrean x value of vavevector; initial value 

Upstream y value of vavevector; initial value 

Expression in Karman spectrum; initial value 

Wavevector amplitude dovnstream at rotor face 

Ratio of dovnstream to upstream vavevector 

amplitudes 

Ratio of upstream to dovnstream velocity 

amplitudes 

Upstream x value of vavevector 

Upstream y value of vavevector 

Upstream z value of vavevector 

Expression in Karman spectrum 

Denominator in Karman spectrum 

Temporary variable for testing end of iteration 

Float N 

Value of x vavenumber 

Value of y vavenumber 

SUBPROGARAMS CALLED 

NOPKVC 

CALL I H G PRO G 2Al-i 

NOPPX 

ERRORS 

IOr(E 

E!ITRY 

Calculate ice vhich equals 0.74683321L equals 0.373417(c/L)/b; 

Define vavevector in -airfoil fixed coordinate system; x along 

Sine and cosine of ang:cs alpha and gamma.' 

Calculate normal to airfoil in 3ame coordinate aystern a3 for DK. 

Calculate rotation matrix. 

chord direction, the second three give the spanvise 

direction and the third three give the normal direction. 

chord, y along span, and z along normal. 

The first three components give the 

DO 50 I=1,3 

DO 30 3=1,3 

Initialize DXDZRl 

DO 30 K=l,3 

EID DO 

Calculate DXDZRL 

EID DO 

DO 70 I=1,3 

EIID no 

DO 60 3=1,3 

Initialize DXDZR 

DO 91 K=1,3 

END DO 

Calculate DXDZR 

END 30 

E:ID DO 

Introduce dummy for CVT 
* 

t 

* 

t 

* 

* 

t 

* 

t 

* 

* 

f 

t 

t 

* 

t 

* 

t 

* 

* 

t 

* 

t 

* 

k 

* 

t 

* 

* 

*** 

~ 

Call deformatron program YOFKVC 

Calculata denominator for xpstream turbulence with kz = 0. 

This vi11 be used for comparisonlater in program. 

l4ake CVT as large as possible without losing accuracy. 

by factors of ~VG and test size of DEL:. 

Initialize sum variable, SVV. 

Begin summation with minimum but still positive value of ZK after 

Initialize counters BP and NU for number of steps up and dovn in kz. 

Sum over increasing values of ZK. 

Increase 

subtracting integer number of CVT values. 

Increment upvard summation counter 

Set third vavevector component 

Call deformation program NOPKVC 

Calculate upstream values of three wavevector components 

Calculate XY, an intermediate variable 

Calculate DN, the denominator of the Karman spectrum 

Add to summation, SVV 

Check for excessive summations 

Increment z vavevector component 

Chck to 3ee hov large the denominator has become. 

End upvard sumnation 

Sua over decreasing values of ZK. 

Increment dovnvard summation counter 

Decrement z wavevector component 

Set third wavevector component 

Call defamation prograa BOPKVC 

Calculate upstream values of three vavevector components 

Calculate XY, an interinediate variable 

Calculate 2t1, the denominator of tha Karman spectrum 

Add to summation, SVV 

Check for excessive summations 

Check to 338 hov large the denominator has Secome. 

End downvard summation 

Hultiply SVV by remaining factors 

EXIT 

18 
AppendiH A: Computer Program for Noise Produced by 

Rot or-Turbulence I n t erac ti on 

A computer program implementing the algorithm discussed in the Theory Manual for Noise 

Calculation is given in the Computer Listings to follow. Computer Listing 1 is the main executive 

portion of the program; it prompts for inputs and prints the output. It also performs the integration 

over the rotor span. 

The inputs to the program, including the symbols used in the program, are: 

Ct 

BN 

C8C 

RC 

R IC 

RPS 

ZM 

FM 

leu 

IM 

SUMS 

PSI 

DH 

AL 1 

L2 

c/L = chord/turbulence length scale 

Blade number 

co/c = sound speedlchord 

Blade radiudchord 

r/c = far-field distancelchord 

Revol u tions/sec 

M, = axial Mach number 

Flow Mach number in rotor plane 

6 = observer angle wrt z axis; observer in x,y plane 

Number of terms in azimuthal summation; N in equation (83) 

$ = angle of Mf wrt z axis; Mf is positive inward as in figure 1 

Frequency step size in fractions of a harmonic 

Initial harmonic number 

Number of frequency steps 

The numerical constant 203.95 occurring in the expression for PSD is obtained by introducing 

equation (82) into equation (63) and multiplying by the blade number. This gives 

where 

Atmospheric pressure Po = 10"0.987 dynes/cm2 and reference pressure P, = 0.0002 dynes/cm*. 

Spp must be multiplied by 417 to convert from a two sided - co < o < co to a one sided 0 < o < Q) 

spectrum and to convert to a unit Hertz rather than a unit radian frequency of measurement. Using the 

relations 

E, = ob/Ub co2 = YP,/po (A2 1 

this can be written (where y = 2nn/N) 

19 
I Converting to dB, 

2 

101Ogl@ [+ ( YnPo pr ) ] =203.95 

The outputs of the program are: 

H 

F1 

P SD 

NP 1,NM l,N 

Frequency of sound expressed in multiple of blade passage harmonic 

Frequency of sound in Hz 

Spectrum level in dB (per unit Hz) 

Counters for certain of the summations for use in diagnosis of output 

The programs in the following appendices are standard Fortran code with added line numbers for 

reference in the documentation. This allows more detailed comments to be used without deleteriously 

affecting the ability to easily read the code. 

The program listing here is for an homogenous turbulence in the rotor plane. This can be easily 

extended to the nonisotropic case by inputting the deformation tensor at each point of the rotor plane 

rather than just once. The turbulence spectrum at each point of the rotor plane is then calculated from 

a rapid contraction of an isotropic turbulence using this deformation tensor. The changes to the main 

program and one subroutine for this generalization to the nonhomogeneous case are listed below. 

Changes to main program: 

12 1 'RADIUS/C, FAR-FIELD/C, RPS, I"/) 

24 IF (IH .EO. 8) READ(l,+) ((DXDZ(I,J), I=1,3), J=1,3) 

25 IF (IH .EQ..B) WRITE(*,2BB) ((DXDZ(I,J), J=1,3), I=1,3) 

26 IF (IH .EQ. B) WRITE(2,2BB) ((DXDZ(I,J), J=1,3), I=1,3) 

28 

59 1 

IF (IH .EQ. 8) CLOSE(1) 

P2,NPl,NMl,M,IH) 

In the subroutine SPX, change the first line of the subroutine to accept the new variable IH and add the 

following line between lines 46 and 47: 

1 SUBROUTINE PX (F 1, ZM, CL ,BN, CBC , UU , RC , R 1C , RPS , FM, TH, J, PSI , DXOZ , 

46a IF (IH .EQ. 1) READ(1,') ((DXDZ(I,J), J=1,3), I=1,31 

In the above code the parameter IH is an input designating whether a homogeneous or inhomogeneous 

calculation should be performed. If IH = 0 is input, homogeneous turbulence is assumed. If IH = 1 is 

input, inhomogeneous turbulence is assumed, and the above read and write statements are deferred to 

the subroutine SPX. 

20 
1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

'/sa 

2 1 

-3 

24 

25 

26 

27 

128 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 

41 

42 

43 

44 

45 

46 

47 

48 

49 

50 

3 3 .L L 

+ 

500 

1 

30 

100 

1 

300 

1 

1 

1 

10 

200 

20 

40 0 

900 

70 

700 

1 

1 

1 

1 

800 

PXFINI 12/3/84 R. K. Amiet Integrates over spon. 

DIMENSION DXD2(3,3) 

CHARACTER*10 FNAME 

CHARACTER"10 FNAMEl 

WRITE(+,500) 

FORMAT(' Turbulence ingestion program with integration over' 

WRITE(','(A\)') ' OUTPUT FILE NAME -' 

READ(*,'(A)') FNAMEl 

OPEN(2,FILE=FNAMEl,STATUS='NEW') 

WRITE(*,100) 

FORMAT(//' INPUT C/L, BLADE NO., C0/C, ' 

READ(*,+) CL, BN, CBC, RC, RlC, RPS 

TM1 = 2.+3.141593*RPS'RC/C0C 

WRITE(*,300) TM1, CL, BN, CBC, RC, RlC, RPS 

WRITE(2,300) TM1, CL, BN, CBC, RC, RlC, RPS 

FORMAT(//' M TIP = ',F4.3,' CHORD/TURB SCALE =',F7.4/ 

' BLADE NO. =',F4.0,' SOUND SPEED/CHORD =',F7.1/ 

' RADIUS/CHORD =',F5.1,' FAR-FIELD/CHORD =', 

F6.1/' REV/SEC =',FS.l/) 

* span'/' ond nonisotropic turbulence. R. K. Amiet 12/84'/) 

'RADIUS/C, FAR-FIELD/C, RPS'/') 

WRITE(*,'(A\)') ' INPUT FILE NAME -' 

READ(*,'(A)') FNAME 

OPEN(l,FILE=FNAME,STATUS='OLD') 

READ(l,*) ((DXDZ(I,J), 1=1,3), J=1,3) 

WRITE(*,200) ((DXDZ(I,J), J=1,3), I=1,3) 

WRITE(2,200) ((DXDZ(I,J), J=1,3), I=1,3) 

FORMAT(lX,3F5.2) 

CLOSE ( 1 1 

WRITE(+,400) 

FORMAT(/' INPUT Mi!, FM, UU, THETA, SUMS, PHI'/) 

READ(*,+) ZM, FM, UU, TH, J, PH 

WRITE(",900) ZM, UU, TH, FM, J, PH 

WRITE(2,900) ZM, UU, TH, FM, J, PH 

FORMAT(' MZ =',F4.3,' UU =',F6.3,' THETA =',F5.0/ 

WRITE(*,700) 

FORMAT(/' INPUT'/ 

' MF =',F4.3,' SUMS =',I4,' PSI =',F5.1/1 

BEGINNING HARMONIC: If 0, go to matrix input.'/ 

FRACTIONAL HARMONIC SPACING: If < 0, go to previous input.'/ 

NUMBER STEPS: If < 0, go to first input.'/) 

READ(*,+) AL1, OH, L2 

IF (AL1 .LT. 0.1 GO TO 10 

IF (OH .LT. 0.) GO TO 20 

IF (L2 .LT. 0) GO TO 30 

OF = RPS'BN'DH 

WRITE(*,800) 

WRITE(2,8001 

FORMAT(/2X,'HRMNC'.SX,'FREQ',6X,'PSD',4X,'SM P' ,lX,'SM M'/) 

F1 = (ALl/DH-l.)*DF 

IF (L2 .EQ. 0) L2=1 

Computer Listing 1: Executive program for turbulence ingestion. 

21 
51 

52 

53 

54 

55 

56 

57 

58 

59 

60 

61 

62 

63 

64 

65 

66 

67 

68 

69 

70 

w 

7-9 

10-13 

14 

15-20 

21-23 

24-28 

29-35 

36-44 

45 

46-48 

49 

50 

51-69 

52 

53-54 

55-63 

56 

57 

58-59 

60 

61 

63 

64 

65-67 

1 

22 

600 

60 

DO 60 K=l,L2 

F1 = F1 + OF 

RCP = RC'1.05 

SUM = B. 

DO 22 I2=0,19 

RCP = RCP - .05*RC 

TM = TMl*RCP/RC 

CALL PX(Fl,ZM,CL,BN,CBC,UU,RCP,RlC,RPS,FM,TH,J,PH,DXD~, 

P2,NPl,NMl,N) 

IF (12 .EO. a) ~2 = ~212. 

SUM = SUM + P2 

CONTINUE 

PSD = 10.*ALOG10(SUM*RC*.05) 

H = Fl/(RPS*BN) 

WRITE(*,60@) H, F1, PSD, NP1, NM1, N 

WRITE(2,608) H, F1, PSD, NP1, NM1, N 

FORMAT(lX,F6.1,1X,F8.l,lX,F9.2,lX,I5,216) 

CONT I NUE 

END 

GO TO 70 

Comment 

Prompt for a filename into which to write the calculated results 

Prompts for inputs 

Calculation of tip Mach number 

Prints the input quantities for checking 

Prompt for input containing the deformation tensor values 

Read and write the values for the deformation tensor 

Prompt for more input parameters; writes out values to file 

Prompt for final input values; if any are negative the program returns to one of the 

previous input prompts 

DF = frequency step size 

Prints headings for the eventual printout 

Initialize the frequency variable 

Input of zero for L2 gives one iteration with this line 

Stepping through frequency range 

Increment frequency variable 

Initialize radius and summation variables 

Integration over rad i us 

RCP = local radius 

TM = Mt local tip Mach number 

Call routine for summation over Y 

Divide the first term in the summation by 2; the last term also should halved, but since it 

approaches zero, it needn't be explicitly set to zero. 

Perform radial summation 

Introduce remaining factors in equation (69) 

Calculate harmonic number 

Write outputs 

Computer Listing 1: Executive program for turbulence ingestion. 

, 22 
Appendix B: Subroutine Calculating the Main Problem Parameters 

and Integrating ouer Azimuth 

This subroutine is called by the main turbulence ingestion program in Appendix A. It calculates 

geometry and integrates over azimuthal angle. The inputs to the subroutine are: 

F1 

ZM 

BN 

coc 

uu 

FZ: 

R1 C 

RPS 

FM 

M 

J 

PSI 

DXDZ 

a 

Observer frequency 

Axial Mach number; see figure 1 

Chordhurbulence integral length scale 

Blade number . 

Sound speed/chord 

RMS turbulence intensity/axial velocity 

Local blade radiuskhord 

Far-fie Id distancekho rd 

Revolutions per second 

Flight Mach number in rotor plane; see figure 1 

Polar angle of observer; see figure 1 

Number of azimuthal integration points 

Angle of flight Mach number with the y axis; see figure 1 

Deformation tensor 

The outputs are: 

P2 

NP1 ,NM1 ,N 

Sound level from a radial blade segment 

Counters used in the summations: used for diagnostic purposes 

The program calls the subroutines TRBNl and LEFF, the turbulence and the airfoil response 

calculations. The subroutine performs the azimuthal integration over 0 in Eq. (63). 

23 
1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 

21 

22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 

41 

42 

43 

44 

4s 

46 

47 

Subroutine SPX 

Program written by R. K. Amiet. 

SUBROUTINE PX (F 1, ZM, CL ,BN ,C0C ,UU ,RC ,R 1C ,RPS ,FM, TH, J, PSI ,DXDZ, 

DIMENSION DXDZ(3,3) 

PI = 3.14159 

TM = 2."PI*RPS"RC/CBC 

AJ = J 

DEL = 360./AJ 

C = COS(TH/57.2958) 

S = SIN(TH/57.2958) 

C4 = COS(PSI/57.2958) 

S4 = SIN(PSI/57.2958) 

QM = FM*S+S4+ZM+C 

SQ1 = SQRT(QM**2+BFZ21 

RER = (SQl+QM)/BFZ2 

NP1 = 0 

NM1 = B 

SUM = 0. 

G = 90. - DEL 

DO 50 I=l,J 

G = G + DEL 

C2 = COS(G/57.2958) 

S2 = SIN(G/57.2958) 

C5 = COS((G+PSI)/57.2958) 

S5 = SIN((G+PSI)/57.2958) 

ALP = 57.2958+ATAN2(ZM,TM+FM+C5) 

C1 = COS(ALP/57.2958) 

S1 = SIN(ALP/57.2958) 

RM = SQRT(ZM""2 + (TM + FM+C5)**2) 

c3 = C'S1 - S'Cl'S2 

X = RlC"(TM"RER"C1 - C3) 

Y = RlC"(S'C2 + FM"RER"S5) 

Z = RlC"(C'C1 + S"Sl"S2 + TM'RER'Sl) 

SG = SQRT(X**2,+(1.-RM*'2)*~Y**2+Z++2)) 

FP = Fl*(l.+TM*(S2"S-FM"RER*C5)/SQl~ 

XK = PI*FP/(CBC*RM) 

YK = PI+FP+Y/(SG*C0C) 

T = 2.*PI*RC/(TM+C0C'BN) 

T1 = T*TM*Sl*Cl/ZM 

1 P2,NPl,NMl,N) 

BFZ2 = l.-FM*+2-ZM*"2 

Next line added to make peaks occur at bpf when ZM = 0. 

T1 = T/(l.+ MF*C5/TM) 

XX = TM*(T-Tl) 

YY = Tl'FM'S5 

ZZ = T*TM"CBC*Sl 

CVT = PI/ZZ 

ZK0 = PI*T2*Fl/ZZ 

CALL TRBNI(CL,CVT,XK,YK,ZK0,ALP,G,DXDZ,SVV,NP,NM,N) 

T2 = T1 + XX"(RM-X/SG)/(l.- RM'"2) + YY+Y/SG 

Computer Listing 2: Subroutine for calculating principle parameters in problem 
48 

49 

50 

51 

52 

53 

54 

55 

56 

w 

1-2 

5 

6 

7 

8-11 

12 

13 

14 

15 

14- 19 

20-53 

21 

22-25 

26 

27-28 

29 

30 

31-33 

34 

35 

36-37 

38 

39-40 

41-43 

44 

45 

46 

47 

48 

49-50 

51 

52 

54 

CALL LFT (X, SG ,RM ,XK , YK ,GL2 1 

NP1 = NP1 + NP 

NM1 = NM1 + NM 

AD = GL2*SVV*RM+(RM+Z+FP*XK/(Fl*SG**2))**2 

' SUM = SUM + AD 

50 CONTINUE 

P2 = 2.48E20*SUM*BN'(UU*ZM)++2/(C0C*A~~ 

RETURN 

END 

Comment 

Inputs 

TM = Mt = local tip Mach number 

Float J 

DEL = azimuthal step size 

C, S and C4, S4 = cosine and sine of 8 and # respectively 

OM = M, cos 0; Eq. (43 

BF22 = 1 - MS2; Eq. (43) 

SO1 = ( 1 - MS2 sin2 @)'I2 

RER = re/r; Eq. (42 

In i t i a1 i ze variables 

Integration over azimuthal angle 7 

Increment 7 

C2, S2 and C5, SS = cosine and sine of 7 and y + Ip respectively 

ALP = a; Eq. (34) 

C1, S1 = cosine and sine of a 

RM = Mb; Eq. (36) 

C3 = cos*; Eq. (49) 

X, Y, Z = ~3, y3, 23; Eq= (48 

SG = Q 

FP = frequency on blade; Eq. (55 

XK, YK = specific wavenumbers I?, Ry; Eqs. (91, ( 14) 

T = time between blade passes 

T1 = time between eddy intersections; Eq. (62) for line 39. Equation (32) is not 

appropriate for eddies stretched in axial direction; for this case the diagonal distance in 

figure 3 is not appropriate; given here is the result for an eddy of infinite axial length 

corresponding to the limit M,+ 0; Eq. (62). 

XX, YY, Zt = X, Y, Z in Eqs. (621, (64) 

CVT = 2n/Z = step size in Eq. (70 

ZKO = ooT2/2 = initial radian frequency in Eq. (70) 

Call turbulence summation subroutine 

Call airfoil response subroutine 

Increment counters 

Contribution to spectrum from azimuthal integration; Eqs, (50) and (69) 

Summation over azimuthal spectrum 

Multiplication by remaining factors in Eq. (69) 

T2 = T2; Eq. (65) 

Computer Listfng 2: Subroutine for calculating principle parameters in problem 

25 
AppendiH C: Subroutine for Summing ouer the Turbulence Spectrum for 

Blade to Blade Correlation 

This subroutine calculates the summation in Eq. (69) over the turbulence spectrum. The output 

SVV of this program is 

svv 

The inputs to the 

CL = c/L 

XK = K, 

program and the equivalent symbols used in the report are 

CVT = 2nb/Z 

YK = Ky 

(El) 

ZKO = aOTZ/t DXDZ(1,J) = i3xi/i3(, (xi = precontraction coordinate; tj = post 

contraction coordinate) 

Other parameters in the program are: EK = k,b where k, is given in Eq. (32). CVT is the step size in 

the summation. Since the spectrum is written in a form with the wavenumbers k, and ky normalized 

by k,, CVT must also be so normalized. Note that this changes CVT in the main program, but this is of 

no concern since CVT is used only in this subroutine. 

The summation over n is over the range *a. In numerically carrying out this summation it is 

easiest to begin where the individual terms are maximum and sum to either side of this point until the 

terms become negligible. Thus, the integer N is the optimum value of n at which to begin the 

summation. The optimum value would be n = 0 except for the offset ooT2/Z occurring in KZ(") (see 

Eq. (69 1). 

If the step size CVT is too small, very many terms in the summation would be required. For this 

case, the summation can be replaced by an integral. For the case of an isotropic turbulence the 

integral can be carried out in closed form as in Eq. (3 1). 

The parameters ZK 1 and ZK2 represent K,(')/ k,. ZK 1 is used for increasing values of n and 

1 ZK2 for decreasing values. 

The integers NP and NM count the number of summation steps in the upward and downward 

directions respectively. If either goes over 100 the summation is terminated, although this case has 

never been encountered. If the output specifies that NP or NM reached the value of 100, further 

checking should be done to ensure that a sufficient number of steps were taken. Generally the 

summation is terminated when the term DN (denoting the denominator in Eq. (32)) becomes greater 

than 200 times the initial value of the denominator. 

Finally, the coefficient 0.1 1561 = 55 r(5/6)/[9$.rrI'( 1/31 4nl. The factor 4n comes from 

Eq. (30) and the remaining factors come from the expression for I. 

26 
1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

1 '4 

15 

16 

17 

18 

1. 9 

?a 

21 

22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 

41 

42 

43 

44 

45 

46 

47 

48 

* EO0YNI.FTN; Summation over turbulence spectrum. 

Written by R. K. Amiet; final version 1985 

SUBROUTINE TRBNI(CL,CVT,XK,YK,ZK,ALP,G,DXOZ,SVV,NP,NM,N) 

DIMENSION EU(3,3), ED(3,3), VN(31, DK(31, DXDZ(3,3),' 

DXDZR(3,3), ROT(3,3), DXDZR1(3,3) 

EK = .373417*CL 

DK(1) = XK 

DK(2) = YK 

DK(3) = 0. 

C1 = COS(ALP/57.2958) 

S1 = SIN(ALP/57.2958) 

C2 = COS(G/57.2958) 

S2 = SIN(G/57.2958) 

VN(1) = 0. 

VN(2) = 0. 

VN(3) = 1. 

ROT(1.1) = Cl*S2 

ROT(2,l) = C2 

ROT(3,l) = Sl*S2 

ROT(1,2) = -Cl*C2 

ROT(2.2) = S2 

ROT(3,2) = -Sl*C2 

ROT(1,3) = -S1 

ROT(2,3) = 0. 

ROT(3,3) = C1 

DO 50 I=1,3 

DO 30 J=1,3 

DXDZRl(1,J) = 0. 

DO 98 K=1,3 

OXOZRl(1,J) = DXDZ(I,K)*ROT(J,K) + DXDZRl(1.J) 

90 CONTINUE 

30 CONTINUE 

50 CONTINUE 

DO 70 I=1,3 

DO 60 J=1,3 

OXOZR(1.J) = Q. 

DO 91 K=1,3 

DXOZR(1,J) = DXDZRl(K,J)*ROT(I,K) + DXOZR(1,J) 

91 CONTINUE 

60 CONTINUE 

70 CONTINUE 

CVTl = CVT 

CALL KVEC(DK,OXOZR,VN,ED,EU,OKM,DKUK,UQDQ) 

XKU0 = EU(1,3)*DKM/DKUK 

YKU0 = EU(2,3)*DKM/OKUK 

XY0 = EK**2 + XKUO**2 + YKU0**2 

00 = XY0**2.83333 

N=O 

10 N=N+l 

CVT = CVTl 

Computer Listing 3: Subroutine for calculating summation over turbulence spectrum. 

27 
49 

50 

51 

52 

53 

54 

55 

56 

57 

58 

59 

60 

61 

62 

63 

64 

65 

66 

67 

68 

69 

70 

71 

72 

73 

74 

75 

76 

77 

79 

I ae 

78 

ai 

a2 

as 

aa 

I 83 

84 

86 

87 

89 

90 

91 

92 

IF (N .GT. 1) CVTl = 2.*CVTl 

DK(3) = CVTl 

CALL KVEC(DK,DXDZR,VN,ED,EU,DKMl,DKUKl,UQDQl) 

XKU = EU(1,3)'DKMl/DKUKl 

YKU = EU(2,3)'DKMl/DKUKl 

ZKU = EU(3,3)*DKMl/DKUKl 

XY = EK*+2 + XKU"2 + YKU+*2 

DN = (XY + ZKU**2)"2.83333 

DEL1 = ABS(1.- DN/DB) 

IF (DEL1 .LT. .el) GO TO 18 

N = ZK/CVT 

AN = N 

ZK1 = ZK - AN+CVT 

ZK2 = ZK1 

NM = 0 

DK(3) = ZK1 

CALL KVEC(DK,DXDZR,VN,ED,EU,DKM,DKUK,UQDQ) 

XKU = EU(l,S)*DKM/DKUK 

YKU = EU(2,3)*DKM/DKUK 

ZKU = EU(3,3)*DKM/DKUK 

XY = EK"2 + XKU"2 + YKU"2 

DN = (XY + ZKU**2)**2.83333 

SVV = SVV + (XY - EK+'2)/(UQDQ"2+DN) 

IF (NP .GT. 100) GO TO 40 

ZK1 = ZK1 + CVT 

IF (ON .LT. 00*200.) GO TO 20 

svv = e. 

NP = a 

2a NP = NP + 1 

40 NM = NM + 1 

ZK2 = ZK2 - CVT 

DK(3) = ZK2 

CALL KVEC(DK,DXDZR,VN,ED,EU,DKM,DKUK,UQDQ) 

XKU = EU(1,3)*DKM/DKUK 

YKU = EU(2,3)*DKM/DKUK 

ZKU = EU(3,3)tDKM/DKUK 

XY = EK*+2 + XKU*+2 + YKU'"2 

ON = (XY + ZKU**2)**2.83333 

SVV = SVV + (XY - EK**2)/(UQDQ"2'DN) 

IF (ON .LT. 00*280.) GO TO 40 

80 SVV = .11561*SVV+CVT*EK**.666667 

RETURN 

END 

IF (NM .GT. lee) GO TO a0 

~ 

Computer Llstlng 3: Subroutine for calculating summation over turbulence spectrum. 

28 

.' 
llnel 

1 

4 

5-7 

8-11 

12-14 

15-23 

24-39 

40-45 

46-58 

59 

68-63 

64-65 

66-77 

78-89 

90 

fanmm 

Inputs 

k, = 0.74683421L = 0.373417(c/L)/b; 0.7468342 = I'(5/6)fi/I'( 1/31. 

Define wavevector k in coordinate system fixed to airfoil; x along chord, y along span, 

and z along normal. k, will take different values. 

Sine and cosine of angles a and Y. 

Normal to airfoil in same coordinate system as in lines 5-7. 

Rotation matrix; dzi/dx3 with the definitions in Eq.( 47). The first three components 

give the chord direction, the second three give the spanwise direction and the third three 

give the normal direction. 

2 1 VN( 1) written in the unrotatedgl system. 

22 VN( 2) written in the unrotated gl system. 

23 VN( 3) written in the unrotatedgl system. 

Calculate deformation tensor in coordinate system defined in lines 5-7. First 

premultiply, then postmultiply, DXDZ by ROT. 

Denominator for upstream turbulence with k, = 0. Used for comparison in lines 77 and 

89. 

If CVT is too small, will take too many summations. Make CVT larger until start getting 

significant variation in either the ratio of wavevectors between downstream and 

upstream. 

Initialize sum variable. 

Begin with minimum but still positive value after subtracting integer number of CVT 

values. 

Initialize counters for number of steps up and down in k, 

Sum over increasing values of k,. Lines 69-7 1 calculate upstream values of 

wavevector. These expressions are placed into isotropic turbulence in 72-83. Line 75: 

check on number of iterations. Line 76: increment k,. Line 77: check to see how large 

the denominator has become. 

Same comments as lines 53-64 except summation is in downward direction. 

Factors multiply sum. Note that CVT is included in the product so that it is o'k to 

increase CVT in lines 46-58 making the summation step size bigger. 

0.1 1561 = ssr(sm/[36~~3/* r( 11311 

Computer Lfttfng 3: Subroutine for cslculatlng summation over turbulence spectrum. 

29 
AppendiH 0: Program for Calculating Airfoil Response Function, L, 

for leading Edge Noise 

Computer Listing 4 gives the subroutine for calculating the effective lift ILl for noise produced 

by airfoil- turbulence interaction. The inputs to the subroutine and the equivalent notation used in 

this report are: 

Proaramm 

X x3 

SG U 

RM Mbiequation (36) 

XK kx; x wavenumber component 

YK k ; y wavenumber component Y 

The output GL2 is the effective lift lLI2 and is calculated in a low frequency and a high frequency 

regime. The low frequency result is given by equation (27) and the high frequency result by Eqs. 

(28). 

* Lift function subroutine 

* Written by R. K. Amiet; final version 1978 

1 SUBROUTINE LFT(X,SG,RM,XK,YK,GL2) 

2 COMPLEX E 

4 UM = RM*XK/B2 

5 RMI = SQRT(RM**2-B2*(YK/XK)+'2) 

6 UMI = UM*RMI/RM 

7 IF (UMI-.75) 10, 10, 20 

8 10 GL2 = l./(B2*(1./(1.+2.4+XK/B2)+6.28319*XK/B2)) 

9 GO TO 30 

3 82 = l.-RM"2 

le 20 T1 = ABS(UM+(X/SG-RMI/RM)) 

11 CALL FRNL(SQRT(1.27324+Tl),E) 

12 CM = CABS(E1 

13 GL2 = .28264*CM**2/(XK*Tl*(l.+RMI)) 

14 30 RETURN 

15 END 

** 

ule 

1 

3 

4 

5 

6 

7 

8 

10-13 

10 

13 

ciunnmh 

Inputs and output 

UM = u 

RMI = M,; equation (24) 

UMI = u-; equation (24) 

Check whether to use low or high frequency solution. 

Low frequency solution; equation (27 

High frequency solution; equation (28a) 

T 1 = Oi; equation (29 

Equation (28a); .20264 = 2/n2 

82 = t3b2 

Computer Listing 4: Subroutine for calculating lift response for leading edge noise 
tlppendix E: Program for Calculating the Effect of Contraction on 

the Turbulence Spectrum 

Computer Listing 5 gives the subroutine for calculating the effect of a contraction of a 

turbulence spectrum. This subroutine has as input a 3x3 deformation matrix representing a flow 

contraction. The program then accepts a flow direction &, and a wavevector k both for the downstream 

flow, and calculates the.upstream Wavevector and flow vector. The inputs to the subroutine and the 

equivalent notation used in this report are: 

Proaramm 

DK 

DXDZ 

VN 

ED, EU 

- kd = downstream wavevector components 

dxi/dtj = deformation tensor; equation (71 1 

Upwash velocity component normal to the airfoil 

Tensors representing the three vector components of the upstream and downstream 

coordinate systems. Each such vector is represented by a column of the matrix. 

1 

2 

3 

4 

5 

6 

7 

9 

11 

12 

13 

14 

1s 

16 

17 

18 

19 

20 

21 

22 

23 

24 

2s 

a 

ie 

The outputs are: 

* 

* 

ie 

38 

20 

68 

50 

48 

Wavevector amplitude downstream at the rotor face 

Ratio of downstream to upstream wavevector amplitudes 

Ratio of upstream to downstream velocity amplitudes 

Turbulence contraction calculation. 

Written by R. K. Amiet 9/7/83 

SUBROUTINE KVEC(DK,DXDZ,VN,ED,EU,DKM,DKUK,UODO) 

DIMENSION E0(3,3). EU(3,3), VN(31, F3(3), OK(31, OXOZ(3.3) 

DKM = SORT(OK(1)**2 + OK(2)**2 + OK(3)*+2) 

ED(I,3) = DK(I)/DKM 

EO(I,2) = VN(1) 

CONT I NUE 

CALL CROSS(ED,l) 

CALL CROSS(ED,2) 

DO le 1=1,3 

DO 28 1=i,3 

00 38 ~=i,3 

EU(1,J) = e.** 

CONT I NUE 

CONT I NUE 

DO 50 J=1,3 

DO 68 K=1,3 

EU(1,J) = EU(1,J) + DXDZ(K,I)*ED(K,J) 

CONT I NUE 

CONT I NUE 

CONT I NU€ 

F3(1) = EU(1,3) 

F3(2) = EU(2,3) 

F3(31 = EU(3,3) 

AF1 = SQRT(EU(1,1)**2 + EU(2,1)**2 + EU(3,1)**2) 

DO 48 1=1,3 

Computer Listing 5: Subroutine for calculating effect of contractlon on turbulence 

spectr um 

31 
26 

27 

28 

29 

30 

31 

32 

33 

34 

3 

4- 7 

8 

9 

10-14 

15-21 

22-24 

25 

26-28 

29 

30 

31 

EU(1,l) = EU(l,l)/AFl 

EU(2,l) = EU(2,1)/AFl 

EU(3,l) = EU(3,1)/AF1 

CALL CROSS(EU,3) 

CALL CROSS(EU,2) 

OKUK = EU(1,3)*F3(1) + EU(2,3)*F3(2) + EU(3,3)*F3(3) 

UQDO = AFI'DKUK 

RETURN 

END 

Camments 

Calculates the magnitude of the downstream wavevector kd. 

Reads yn into the 2'nd column of ed and a unit vector along kd into the 3'rd column. 

Calculates the direction of the vorticity determined by the specified wavevector k and 

velocity Vn directions. 

Calculates the direction of the velocity produced by the wavevector. !In is a component of 

this velocity. 

Initialize matrix eu. 

Calculates the values of eu from the downstream values ed by using equation( 7 1 ) on each of 

the 3 column vectors making up ed. This is not the final form for eu. The vector GIu 

representing the vorticity direction will keep its direction but be normalized to a unit 

vector. GSu will be found by taking the cross product of vectors 1 and 2, and GZu is 

then determined by taking the cross product of GSu and GIu. 

- F represents a vector found by taking a unit vector alongkd and transforming according to 

equation(6 1 lb). This is used to determine kd/ku using equation(72). 

The magnitude of the transformed value of the vorticity, needed in equation(73). 

Unit vector in the direction of the upstream vorticity. 

Determines the unit vector isu from the cross product of ilu and&". 

Determines the unit vector 

Determines kd/ku using equation( 72). 

from the cross product of and diu. 

Computer Listing 5: Subroutine for calculating effect of contraction on turbulence 

spectrum ~ 

I v' 

32 
AppendiH F: Subroutine for Calculating the Cross Product of Two Uectors 

Computer Listing 6 gives the subroutine for calculating the cross product of two vectors. The 

subroutine takes a 3x3 matrix E(i,j) composed of 3 column vectors ii and replaces the Nl'th vector 

by the normalized cross product of the other two in proper cyclic order.The inputs to the subroutine 

and the equivalent notation used in this report are: 

Proaramm 

EM) 

N1 

Tensor with columns composed of vectors: also represents the output 

Column not involved in the cross product b 

* 

+ 

1 

2 

3 

5 

6 

'I 

t3 

3 

10 

11 

12 

13 

14 

15 

4 

Subroutine for calculating cross product of two vectors 

Written by R. K. Amiet 9/1/83 

SUBROUTINE CROSS(E,Nll 

DIMENSION E(3,3) 

N2 = N1 + 1 

N3 = N1 + 2 

IF (N2 .GT. 3) N2 = N2 - 3 

IF (N3 .GT. 3) N3 = N3 - 3 

E(1,Nl) = E(2,N2)+E(3,N3) - E(3,N2l*E(2,N3) 

E(2,Nl) = E(3,N2l*E(l,N3) - E(l,N2)*E(3,N3) 

E(3,Nl) = E(l,N2)*E(2,N3) - E(2,N2)*E(l,N3) 

ENlM = SQRT(E(l,N1)+*2 + E(2,N1)**2 + E(3,N1)+*2) 

E(1,Nl) = E(l.Nl)/ENlM 

E(2,Nl) = E(2,Nl)/ENlM 

E(3,Nl) = E(3,Nl)/ENlM 

RETURN 

END 

WfhuMu 

2-6 

7-9 

10- 13 

Defines the cyclic sequence 1,2,3 or 2,3,1 or 3,1,2 where N 1 is the first number of 

the sequence. 

Replaces the N l'th vector by the cross product of the other two. 

Normalizes the Nl'!h vector to have a magnitude of 1. 

Computer Listing 6: Subroutine for calculating cross product of two vectors 

33 
Appendix G: Subroutine for Calculating the Fresnel Integral 

Computer Listing 7 gives the subroutine for calculating the Fresnel integrals. The program uses 

the algorithm given by equations 7.3.9, 7.3.32, 7.3.33 of Abromowitz and Stegun, Handbook of 

Mathematical Functions, Dover Publications, New York, 1968. 

1 

2 

3 

4 1 

5 

6 

7 

8 

FRNL.FOR Subroutine for calculation of Fresnel integrols 

Written by R. K. Amiet 1976 

SUBROUTINE FRNL(X,E) 

COMPLEX E, G, H 

G = CMPLX((1.+.926*X)/(2.+ 1.792"X + 3.1@4*X**2), 

H = CMPLX(SIN(1.5708*X**2~,COS(l.57~8*X**2)) 

E = G*H + CMPLX(.5,-.5) 

RETURN 

END 

1./(2.+ 4.142*X + 3.492*X**2 + 6.67*X**3)) 

, 

Computer Listing 7: Subroutine for calculating the Fresnel integrals. 

34 
. 

Appendix H: Sample Calculation 

These two test cases use the same inputs, listed below. The only difference between the cases is 

the deformation tensor. For the first case the tensor represents no deformation. The second 

represents a deformation by a factor of 4 in t h e a x i a l direction. 

Rev/sec = 10. 

Blade# = 2. 

Chord/Turb scale = .0 1 

Sound speed/Chord = 1000. 

Mach no. Mf = 0. 

Far field/chord = 100. 

Radius/chord = 10. rms turb vel/U = .010 

Axial Mach no. = M, = 0.1 

Case 1 

Deformation tensor = 

1.00 .OO .OO 

.OO 1.00 .OO 

.OO .OO 1.00 

HRMNC FREQ 

.5 10.0 

1.0 20.0 

1.5 30.0 

2.a 40.0 

2.5 50.0 

3.0 60.0 

3.5 70.0 

4.0 80.0 

4.5 90.0 

5.0 100.0 

PSD 

31.94 

62.82 

44.56 

56.06 

47.12 

51.93 

47.24 

49.11 

46.55 

47.09 

HRMNC 

5.5 

6.0 

6.5 

7.0 

7.5 

8.0 

8.5 

9.0 

9.5 

10.0 

Theta8 = 0. 

FREQ 

110.0 

120.0 

130.0 

140.0 

150.0 

160.0 

170.0 

180.0 

190.0 

200.0 

PSD 

45.62 

45.58 

44.66 

44.39 

43.73 

43.36 

42.84 

42.45 

42.00 

41.62 

HRMNC FREQ PSD 

.5 10.0 12.34 

1.0 20.0 62.08 

1.5 30.0 28.14 

2.0 40.0 55.32 

2.5 50.0 33.18 

3.0 60.0 51.06 

3.5 70.0 35.89 

4.0 80.0 47.95 

4.5 90.0 37.82 

5.0 100.0 45.51 

HRMNC FREQ 

5.5 110.0 

6.0 120.0 

6.5 130.0 

7.0 140.0 

7.5 150.0 

8.0 160.0 

8.5 170.0 

9.0 180.0 

9.5 190.0 

10.0 200.0 

35 

PSD 

39.34 

43.53 

40.56 

41.88 

41.56 

40.45 

42.34 

39.26 

42.91 

38.29 
Report Documentation Page 

1 Report No 2 Government Accession No 

NASA CR-181790 - ' 4 Title and Subtitle 1 Noise Produced by Turbulent Flow Into a Rotor: 

3 Recipient's Catalog No 

5 Report Date 

June 1989 

1 User Manual for Noise Calculation 

9. Performing Organization Name and Address 

United Technologies Research Center 

East Hartford, CT 06108 

12. Sponsoring Agency Name and Address 

National Aeronautics and Space Administration 

Langley Research Center 

Hampton, VA 23665-5225 

8. Performing Organization Report No. 

R. K. Amiet, C.G. Egolf and J.C. Simonich 

____.- 

10 Work Unit No 

505-63-5 1 

__- 

11. Contract or Grant No 

NAS1-17763 

13. Type of Report and Period Covered 

Cohtractor Report 

14 Sponsoring bgency Code 

__. 

6. Performing Organization Code t17. Key Words (Suggested by AuthorW) noise prediction 

helicopter noise, turbulence ingestion 

noise, rotor noise 

18 Distribution Statement 

Unclassified - Unlimited 

Subj ect Category 7 I 

19. Security Classif (of this report) 20 Security Classif (of this page) 21 No of pages 

Unclassified Unclassified 37 

22 Price 

A0 3