Het transformeren van een willekeurige gesloten contour in een cirkel m.b.v. de transformatie van theodorsen

3b. 'i. cheepsbo0

TecisC

HoqeschoOl

Deft

HET TRANSFQRMEREN VAN EEN-.WILLEKEURIGE GESLOTEN CONTOUR IN EEN CIRKEL M.]3.V. DE TRANSFQRMATIE VAN THEODORSEN.

ARCHIEF

juli 1985 JM Koning

:Adm. Ee1rìchsinge1 66

Inhoud. Inleiding De Joukowski Transformatie De Theodorsen Transforrnatie. 3.1 De Transformatie. 3.2 Het iteratieproces

3.3 Bepaling van de integraal. Het Programma.

Resultaten.

1. Inleiding

In dit versiag wordt bekeken hoe een willekeurige gesloten contour getransformeerd wordt in een cirkel.Daarbij wordt gebruik gemaakt van de Theodorsen transformatie.Deze transformatie werkt echter niet als we te maken hebben met slanke contouren en er moet dan

eerst een andere transformatie toegepast worden:De Joukowski

transformatie. Deze transformeert de contour bijna in een cirkel en dat wordt glad gestreken als daarna de Theodorsen transformatie wordt toegepast.

In de eerste paragraaf wordt de Joukowski transformatie besproken en in de tweede paragraaf de Theodorsen's.Daarna wordt het

programma beschreven en tenslotte een drietal voorbeelden bekeken,n.1. een vierkant,een straalbuis en een spanvorm die willekeurig gekozen is.

-1-2.De Joukowski Transformatie.

Deze transformatie ziet er als voigt uit:

(1)

waarbij =x+iy het physische viak voorstelt z'=a*exp(ie) het getransformeerde viak a een constante nog nader te bepalen

Als we nu (i) uitschrijven:

x+iy=a*exp(+i0)+a*exp(_(c+i0))=2acosh(+j0)

=2acosh( cose +2a is inh( s in®

Bekijk het reele en het imaginaire deel:

x=2acosh(,cos0 (2)

y=2asinhcsin0 (3)

Ais nu en ® bekend zijn dan hebben we hiermee meteen x en

y.Echter orn en 9 te bepalen moeten de twee steisels opgelost worden.Uit (2) en (3) voigt: X

cosh=

2 acose M.b.v. cosh-sintd=1 voigt: 2 2 2 2 -

2.2

4a cos 9 4a sin O Dit verder uitwerken:

2.2

2

.2

2.2

2

x sin 0-y (i-sin )=4a sin 0(1-sin

_e)

sin20_((X/a)2+(y/a)2_4)sin20_1/4(y/a)2=Q y

en _{sinh_2asjn0}

-2-1 2 2 .

.2

Stel p={

(x/a)

+(y/a) -4} . Als we dit nu oplossen voor sin

c en er

rekening mee houden dat sin29 groter dan nul is:

2 2 2

2sin e=-p1(p i-(y/a) )=b

We krijgen nu dat e=±arcsinl½b. Er moet wel rekening mee gehouden worden dat,vanwege de arcsin,e ugt tussen -it/2 en it/2. Als dit punt niet in dit gebied ugt moet er t bij opgeteld worden orn de

goede S te krijgen.

Op dezelfde manier kan bepaald worden.

2 2

X _y

+

2 2

2.

2

4a cosh 4) 4a sinh

Dit verder uitwerken geeft:

2.

2 2 . 2

2.

2

. 2

x smb +y (1+smnh

)=4a smrth (1+smnh )

sinh4( (x/a)2+(y/a)2-4)sinh2J/4(y/a)2=O

Dit nu oplossen voor sinh24:

2sinh2=p+1(p2+(y/a)2)=c

si nh =±_1½c

Uit de definitie van sinh voigt:

expc4)-exp(-q))=±2/½c

Vermenigvuldigen met exp geeft:

exp(2) 21(½c)exp-i=O

exp=±1c+1 (½c+1)

-1

-4-Dus: (k=ln{±/½c+/(½c+l)} (5)

We hebben nu twee waarden voor en e.Vanwege het feit dat

(/½c+/(½c+1))(-½c+/(½c+l))=1,correspondeert het punt (x,y) in het ene geval met het punt (a*expc,9) en in bet andere geval met

(a*exp(_,_0)

De transformatie (1) zorgt ervoor dat een lijn dt loopt van -2a naar 2a overgaat in een cirkel met straal a.Het gebied buiten de lijn 4a wordt eenduidig getransformeerd in het gebied buiten

(binnen) de cirkel.Zo zal jeder punt uitwendig een gesloten contour getransformeerd worden in een punt buiten de bijna cirkel met

straal a*exp of inwendig een andere bijna cirkel met straal a*exp(_).Vanaf nu zal alleen het eerste geval bekeken _worden. Nu moet de constante a bepaald worden.Daarvoor wordt gebruik

gemaakt van de kromtestraal p .Als we ervan uitgaan dat de contouren bijna ellipsen zijn,kunnen we de kromtestraal gebruiken _{van een} ellips.Deze is namelijk vrij eenvoudig te bepalen.Stel we hebben de vergelijking voor een ellips: (x/xO)2+(y/yO)2=1 en als

parametervoorstelling wordt genomen x=(xOcosc3,yOsine) .De

kromtestraa.1 in bet punt (xO,0), wat overeenkomt met 9=0 wordt _dan:

1(xI)2+(y1)213/2

p =

x'y-xy'

_xO

We nemen nu xO=2acosh en yO=2asinh dan voigt p2ack2. Verder geldt in het punt (xO,0) dat xO=2acosh=2a(l+½2)=

2a+½p.Hieruit voigt dan de constante a:

1

11

1

3.De Theodorsen Transformatie 3.1 De transformatie.

0m een gesioten contour te transformeren in een cirkel wordt gebruik gemaakt van de voigende transformatie:

z'=zexp( c ¡z

o n

n'

(6)

waarbij z' de contour voorstelt in het physische viak z de cirkel voorstelt

cn=An+iBn de constanten voorstelt die bepaald moeten worden

Als z= dan voigt dat z'=zexpc0.Omdat we willen dat de twee gebieden in het oneindige gelij]c zijn voigt dat c0=O.

Stel dat z'=aexp(+ie) en z=aexp(0+i).We moeten er wei aandenken

dat de a, en _{d.ezeifde waarden hebben als in paragraaf 2,ais}

eerst de Joukowski transformatie wordt toegepast.Wordt die niet toegepast dan nernen we a geiijk aan 1.

We transformeren dus een contour in een cirkel met straal

_aexp0

en met de oorsprong als middeipunt.Ais we nu hiervan gebruik maken en invullen in (6) voigt:

"t,

Neem de iogaritme hiervan:

+i®=0+i+(A+iB)

Schrijf nu z=R(cosc+isind?) met R=aexp0:

(-0)i(e-)=(A +iB )R_n(cosn_isinn)

n n

1

Bekijk het reeie en imaginaire deei:

i

2it A /Rt

_f

cosnd

(9) n it o 2it B /Rr=

_f

sinnd

(10) n it o 2it

f

(11)

Uit (il) voigt dat de lengte van de contour niet verandert. Substitueer (9) en (10) in (8) dan voigt:

waarbij (0-e) een functie van is en we ' gebruiken orn een

verschil aan te duiden met is constant als we integreren.We kunnen dus cosn' en sinn' onder de integratieteken zetten.

i 2it

=-f

(sinn(-)d

o

Er geidt nu de voigende relatie:

n (2n+1)

cosf

sin-' )=½cot(1')

_2sin[

_1)/2]

1

Orndat we te rnaken hebben met een oneindige reeks,krijgen we:

\ /2it

_cos{

(2n+1))

₂

(0_)1im{!

_n

_{f

)cot(')d-f

₍₎

2rc

De eerste integraal is onafhankeiijk van n en de tweede integraal

JE is nui.Stel -0Ece(). Dan geldt er:

-6-e-=

(B/R"cosn-A/Rsinn)

(8)

i

We kunnen nu de reeksen opvatten als Fourrierreeksen.Uit (7) kunnen de constanten _OIAn/R en B/R bepaald worden.

-7-2ic

c(')=f

()cot' ₂ )d (12)

o

Tot nu toe hadden we A/R en B/R bepaald uit (7).We kunnen

echter oak uit gaan van (8).A/Rrl en B/Rr worden nu met -O=c():

2it A /R'= _f

E()sinnd

n it o 2-rc

B /R=-- f

()cosnd

n ic o

Dit invullen in (7) geeft:

2it

i

I

f

c(c){sinntcosn -cosnsinn }d+c0

lt10

_ 2ic

f

10

Op dezelfde manier als het voorgaande voigt:

-

_27ti f

2ic

()cot(' )d+0

(13)

We moeten c(') en ,(4) bepalen uit (12) en (13) en dan is de transformatie rond.Voor het bepalen van (') zullen we

weglaten.Dit kan omdat toch constant is.Voor het bestaan van de integraal in (12) moet

()

continue en differentieerbaar zijn.

Eerst moet er een notatieafspraak worden gemaakt.Als of

bekeken worden als functie van schrijven we c() en E(). Als ze bekeken worden als functie van e schrijven we i(e) en (e)

Per definitie geldt:

(e)( (e))

Omdat 4,-G hebben we verder:

0(4) )4,-c (4)

4, (0 ) =e+

(e)

(15)

Dit wil zeggen dat orn (4,) en E(4,) te bepalen uit i(e) en i(e) alle

punten horizontaal worden verschoven met een afstand (e).Omgekeerd geldt precies hetzelfde.Nu worden alle punten verschoven met een

afstand -c(4,).

Eekijk (12) en ga partieel integreren.

2ic

2 d4,

Bekijk log(sin(4,4, )).Deze is reeel in het interval (4,',4,'+2)

Echter er wordt het interval (O,2it) beschouwd en er wordt nu alleen het reele deel meegenomen.Als we nu weten dat 4,=0+(0)

kunnen we het volgende schrijven:

.( (®+1 )-(9+c1) ') 2

1og(sin=1Qg(5jn(00'))+g(

sin'

(0+)-(0+E)' sin 2 2 log

(e)-(e

+ +10g

(&+)-(0+i)'

C)

Sin

_sin

2 2

Vanwege _(0):(4,(9)) geldt

)4,_(e)9

Vermenigvuldig nu

het linkerlid van (16) met d(4,) het rechterlid met

I

d(0)9

_{en integreer van}

_{O naar 2:}

7t dO

-8-)+

2t

-c( ')(3'

)=

_lt0

f iog(sinO® )d(e)dO+ 2-rc (0+

i

_f

sin --rc log

d(e)de+

+ (e+ ) - (e+ E) 27t sin dO 1 f log 2 di(e) -

(e+

)-(e+

)'

dO n n sin Als we nu definieren:

(e+Ï1)_(e1)'

)d(e)de _{,k=1,..N (17)} (®')=!

_f

_log(sin 2 k it o krijgen we met

(e'

+«-

_n

Uit (17) zien we dat verkregen wordt uit .We definieren nu c1(0)=O+Ek(0) k=1,..,N .Dan geldt:

2ic

kk

_d(k)

* 1

Ek+1(o)=ck+1(k)=- f iog(sin

2 2

dk

dk

o

Door nu (18) te bepalen voor iedere k=1,..N krijgen we uiteindelijk

(e'). Uit (14) en (15) voigt dan c() en In de volgende

paragraaf wordt het hele proces stap voor stap besproken. Eerst bekijken we echter de condities waaraan voldaan moet worden.Omdat we uitgaan van een conforme afbeelding en er een

een-eenduidige verband moet bestaan tussen een punt in het z'-vlak en een punt in bet z-vlak,moet O(e) een monotone toenemende functie zijn,dus =1 dc(4)>Q Dit correspondeert met

Zo ook moet (e) een monotone toenemende functie zijn,dus

dt(e)

of dO 2 -9-(18)

w'

Orndat we willen dat bij jeder stap in bet proes er een een-eenduidige verband bestaat tussen 0 en die we in iedere stap opnieuw bepalen moet er dus gelden:

We moeten er wel aandenken dat j(e) een enkelvoudige functie moet zijnen dat daardoor ook (pCe) enkeivoudig is.Als we nu

uiteindelijk (p') en c,4)((p') bepaald hebben kunnen we mbv. (7) en

(8) leder punt in het z'-vlak afbeelden naar een punt in bet z-vlak en omgekeerd.Stel we hebben een (4 en O gegeven dan krijgen we

een 1) en (p met:

3.2 .Het Iteratieproces.

Het hele proces is gebaseerd op (18).Hierbij gaan we als voigt te werk:We veronderstellien dat i(o) een bekende functie is met

periode 2'it.M.b.v. (17) kunnen we dan 1(t3') bepalen,omdat we

initialiseren op nui.Hoe de integraal berekent wordt komt in de volgende paragraaf ten sprake.Als we c1(O') kennen voigt uit

(15) (p:(p=O+(0) en zo krijgen we 4)(4). Als we het grafisch zouden willen tekenen,betekent dit dat jeder punt van e) met een afstand ii(0) horizontaal verplaatst wordt.We zullen later het grafische proces voor een vierkant wat nader bekijken.M.b.v. (18)

*

kunnen we 2(O')Ec2((pj) bepalen en hierdoor hebben we

(p2:(p2e2(0)U1t (p2 kunnen we dan c2((p2) en _{(2) bepalen,}

(19)

io--4) ( _)+t1)0 (20)

d.w.z. jeder punt van i(o) wordt nu met een afstand c2 horizontaal verplaatst.

Wederom gebruiken we (18) orn 3()E3(e) te bepalen en zo gaan we

door totdat niet meer verandert of wat hetzelfde is totdat het verschil tussen _{k+l(e) en} k(e) naar nul toe gaat.

Als we te maken hebben met slanke contouren en we dus eerst een Joukowski transformatie doen,dan hoeven we slechts 2 à 3 stappen te doen orn goede resultaten te krijgen.

De keuze van o(e)O is belangrijk;Dit betekent namelijk

dat en e bijna gelijk zijn en dat betekent dat a*exp(+ie) al bijna een cirkel representeert.

3.3 Bepaling van de integraal.

In het proces moeten we ledere keer een integraal uitrekenen,die in een algemene vorm te schrijven is als:

2,t

i

r

XX

J log(sin _dx 2iî x-x

-4-f(x)cot

2 dx

Voor x moeten we lezen e of

0m (23) te bepalen moeten we rekening houden met het singuliere punt x' en de integraal in een interval rond x' apart bekijken.We

splitsen de integraal op in 3 stukken:

x'-s1 , xI+s2 _2-t x-x x-x

c(x )-

(x)cot 2 dx+ f (x)cotX_C'dx f (x)cot 2 dx} 2it o 2 (24)

waarbij s1 en s2 klein zijn.

3.3.i.Integraal buiten het interval met het singuliere punt. We gaan de intervallen (O,x'-s1) en (x'+s2,2) opsplitsen in deelintervallen en in jeder deelinterval wordt (x) lineair benader.Bekijk nu een deelinterval [x11x2]. dan geldt

voor

(23)

-il-( X2-X1 X2

_XX'

x-X

Dus: f ((X)cot 2

dxf

cj( cot( 2 )dx+ xi xl (x2)-(x1) X2 _{x_X} x2-X1 f (x-x1)cot 2 dx xi x2 ,(x2)-c(x1) X2

2(x')[1nHs1n

ii] {f

(x-x'jcot

₂ dx+

XX

2 x1 x2-x1 xl X2 X-x f (x'-x1)cot ₂ dx} xi k 2k2 1+2k

Er geldt dat fycotydy=(-1) _(12k)(2k)! y voor Iytt (25)

o

waabij _2k de Bernouilli coefficienten zijn. X-X'

Omdat x ugt tussen O en 2lt ,ligt

2 tussen

-it

en ,t en magen we (25) gebruiken. x2 -x x-_z sin f

c(x)cot

dx=2q(x)1nII 2 _*

-x

x2-X1 X

_sin

--1 -i) 2k

I

B22

1+2k xi-x 1+2k {4 (1+2k)(2k)! 2 2 x2-X' sin 2 +2(x'-x1)lniI Xi-X' i}

sin

₂ -12-2k (26)

3.3.2 Integraal binnen bet interval met bet singuliere punt.

We ontwikkelen cj.I(x) in een Taylorreeks rond x=x':

x,s2

_x-x'

I

x-x

_)-(x')

XX

f

((x)cot(

2 )dx=f {(x')cot 2 dx+(x-x' dx cot 2

X'_S]

XI-si

M.b.v. (25) voigt: s2 xI+s2 I

sin-x-x

f

(x)cot 2 dx=2(x')lnII 2

-Si

XI-s1

sin-.----2k k 2k2 s2 1+2k S1 1+2k 4-(x' )(-1_dx (12k)(2k)!{ (----) } (27)

We moeten

er

wel rekening mee houden dat fycotydy onnauwkeuriger wordt naarmate Ilyll naar 7E toegaat.Daarom integreren we niet van O

naar 27c,maar zodanig dat bet singuliere punt in het raidden van het interval ligt.Verdelen we nu [O,2] in M-1 deelintervallen dan

rangschikken we zo dat er [-] deelintervallen links en rechts van het singuliere punt liggen.Het interval mag verschoven worden omdat we te maken hebben met een periodieke functie onder de integraal. Als we nu (26) en (27) samen nemen krijgen we:

2it

X X'

c(x')

_Ito

_f

(x)cot

₂ dx -1 [(M-1)/2] M-1

=-{

+ }{2(x')lnu i _{[ (M21/2]} i_+1 i B 22k i+1 i 2k i X -x )1+2k x -x 1+2k

(x

)-(x

_{)*[4 (_1)k(l+2k)(2k),l(} 2 2 i+1 i X -x O sin sin 2 2 +2(x'-xt)lnji T]}+2(x')1n L(M_2)/2j_,

i

I

x-x

X

sin

₂

sin

₂ 2k k B2k2 [(M+2)/2J

+4(x' )(-i

2 -x')i+2k

[(M-2)/2],

X 2 )i+2k} -14-i+1 X -x sin

il

X -X

2 11+ (28)

4.Het programma.

Nu de theorie achter de rug is kan een programma geschreven worden die een contour transformeert in een cìrkel.Vanwege de moeilijkheid wordt er eerst een structuurdiagram gemaakt.

Het hele iteratieproces hebben we beschreven in paragraaf 3.2,en door dit te volgen kunnen we een structuurdiagram maken.Voor de Joukowski transformatie gebruiken we (4) en (5) en voor het bepalen van de integraal maken we gebruik van (28).Allereerst een tweetal opmerkingen.

Men moet de y- en z-coordinaten van de contour invoeren _{en wel} zodanig dat ze linksom gerangschikt zijn.Het eerste en het laatste punt moeten met elkaar overeenkomen.Ook moet men invoeren of men met een slank contour te maken heb of niet.

In de subroutines maken we gebruik van twee standaard routines, n.1. icspin en ICSEVU.Deze bepalen de cubic spline coefficienten met periodieke randcondities en de evaluatie van de cubic spline.

Op de volgende bladzijden volgen de structuurdiagrarnmen,die steeds verder verfijnd worden.

-15-

-16-Doe de joukowski (1)

transformatie

Bepaal O en _{(o) rechtstreeks uit} de y-,z-coordinaten. (2)

Bepaal m.b.v. de trapeziumregel. (3)

Begin met 0(0)=O.Dus 1=e1

Bepaal de intervaispunten O die nodig zijn voor de berekening

van de integraal in de punten (4)

Bepaal de integraal _{(G') (5)}

Voor i=l,M

Bepaal de nieuwe intervaispunten en bereken in die punten (4)

Voor k=l,N

- i

*

i

Bepaal de integraal k1(O (5)

Voor i=l,M

k+1=0k+1

(ei)

Bepaal de nieuwe intervaispunten en bepaal in die

punten ₍₄₎

B epa lin ìavan de f o

Ga loop uit

Bepaal cubic spline coefficienten met de punten

+1n+1n+1

Bepaal cubic spline coefficienten met de punten

(1) De Joukowski Transformatie

-17-(la) Bepaling van de hoek e1 Voor i=l,M-1

wort=

soml= _[(y/a)2+(z/a)2_4]+wort

Bepaling van e1 en ervoor zorgen dat de hoek,afhankelijk van de y-,z-coordinaat in de juiste kwadrant komt te

liggen _(la)

som2= _{(y./a)2+(z/a)2_4 /4+wort}

(e1 )ln[ / (½som2 )+/ (½som2+l)]

01=arcsin/ (½soml)

Y>O, zO

F + 2 1 F

T/

yio, z

(2) Bepaal en i(e) rechtstreeks uit de coordinaten

Berekening van de constante

De intervaispunten en -18-Voor i=1,M 81=arctan(z./y.) Voor i=1,M (®1)=1n[/[ (y.)2+(z1)2]] psiO=O Voor i=2,M i i-1 i i-1

psiO=psiO+-(s +0

_{)((e )(e}

₎₎

psiO=psio/211

Voor i=1,M-1

Bepaal de cubic spline coefficienten met de _punten

(5) Berekening van de integraal

-19-Voor i=1,m-1

Bepaal de Bernou11e coefficienten _B2k Voor k=O,23

coef(k)=(_1)k22kB2k/((12k)(2k)L)

Rangschik de intervaispunten en _{(G) zodanig dat}

midden in het hele interval ugt Voor j=1,M-1

pg(j)=lnhIsin[[e1_]/2]/sin[(e_

_)Jii

del(j)=[(S1 )(0J)J/(0J+l0J)

sorn(

_j

)=O Voor k=O,23 som(j)=som(j)+coef(k)([ (eJ+l_1)/2]l+2k_

j

i 1+2k hulpl=O hulp2=O Voor j=l,(M-3)/2

hu1p1=hulp1+2 (e)*pg(j)+del( j)*(4som(j

)+

2(_0)*pg(j))

Voor j=(M+1)/2,M-1

hu1p2=hu1p2+2( e) *pg ( _{j )}+del(

_j

_{)+ ( 4som(}

_j

₎₊

2(_O)*pg(j))

hu1p3=2()*pg[(M_1)/2]+de1[(M_1)/2]*som[(_u)/]

In Appendix i hebben we _{een listing van bet programma en hieronder} voigt een lijst met symbolen die _{we gebruiken in de subroutines}

van het programma.

(A) De subroutine TTRANSF die een gesioten contour _{transformeert} in een cirkel m.b.v. _{Theodorsen transformatie,al dan}

niet voor-afgegaan door de Joukowski _{transformatie als we te hebben}

met een slanke contour.

Aanroep:Call TTRANSF(Y,Z,L,A,TETA,PHI_{,PSI IPSIO,EPSTET,CPSITET,} CPSIPHI , CEPSTET, CEPSPHI)

Invoer:

Y(M+2),Z(M+2): y-,z-coordinaten _{van het physische viak,} waarbij M het aantal punten voorstelt op de contour.

L: L=O _{we hebben te maken met een slank}

contour en er wordt eerst een Joukowski _{transformatie gedaan.}

L=l _{we doen meteen de Theodorsen}

transformatie. A: constante die de brandpuntafstand is van een ellips

Uitvoer:

TETA(M): hoek e in het physische _vlak.

PHI(M): hoek in bet getransformeerde _viak.

PSI(M): Als L=1 dan de logaritme _{van de straal R van}

_{een punt}

in het physische viak.

Ais LzO dan de logaritme _{van R/A van} een punt na Joukowski transformatie. PSIO: de constante

EPSFTET(M):de integraal _{£k+1 als functie van O} CPSITET(M-i,3): cubic spline _{coefficienten van PSI}

als functie van g

CPSIPHI(M-1,3): cubic spline _{coefficienten van PSI als} functie van

-20-CEPSTET(M-1,3): cubic spline coefficienten _{van EPSFTET als} functie van e

CEPSPHI(M-1,3): cubic spline coefficienten _{van EPSFPHI als} functie van

De andere variabelen die in de subroutine _voorkomen:

PI: it

TETA1(M): de intervaispunten _{waarover geintegreerd worden.} Ze bevinden zich te midden van twee singuliere punten.

EPSTET1(M): de integraal _{Ek als functie van e} EPSFPHI(M): de integraal

k+1 als functie van

FOUT(M): Ekl(o)_ck(e)

PSIFTET(M): _{als functie van e}

(B) De subroutine INTGR die de _{integraal bepaald van} in de punten '.

Aanroep: Call INTGR(PHI,TETA1,.PSIFTET,PSI,EPSTET2)

Invoer:

PHI(M): het singuliere punt

TETA1(M): de intervaispunten waarover geintegreerd worden. PSIFTET(M): in de intervaispunten.

PSI(M): 4 in het singuliere punt '.

Uitvoer:

EPSTET2(M): de integraal _{Ek als functie van}

De andere variabelen die in _{de subroutine voorkomen:}

B(O:23): de Bernouille _{coefficienten.}

PI: it

-21-FAC(O:23): kZ, k=O,23

COEF(O:23): (_l)k(22kB(k))/((l+2k)ki) k0,23

TETA2(M): de intervalspunten,zodanig gerangschikt _{dat p'} midden in het hele interval ugt.

HPSIFT(M): de bijbehorende in de punten TETA2(M).

PG(M-l): integraal van cotan[(-')/2] _{geintegreerd over} een deelinterval.

DEL(M-l):

HSOM(M-l): berekening van de integraal y*cotany,geintegreerd over een deelinterval voor een Bernouilli

coefficient.

SOM(M-l): berekening van de integraal y*cotany _geintegreerd over een deelinterval voor alle Bernouilli

coefficienten.

HULP1: integraal van _{tot aan interval met singuliere punt} HULP2: integraal van _{na het interval met singuliere punt.} HULP3: integraal van _{van het interval met singuliere punt}

(C) De subroutine JTRANSF die de Joukowski transformatie toepast als we te maken hebben met slanke _contouren.

Aanroep: Call JTRAJSF(Y,Z,A,TETA,PSI)

Invoer:

Y(M+2),z(M+2): y-,z-coordinaten _{van het physische vlak.} A: constante die

de

_{brandpuntafstand is van een ellips.}

Uitvoer:

TETA(M): de hoek 9 na de Joukowski _{transformatie.}

PSI(M): de logaritme _{van R/A van een punt na de Joukowski} transformatie.

-22-De andere variabelen die in de subroutine voorkomen: PI: it WORT: SOM1:-SOM2: 1/4((y./a)2(z./a)2_4)+wort

-23-5. Resultaten.

Orn het programma te testen hebben we een drietal contouren genornen, die hieronder stuk voor stuk behandeld zullen worden.

een vierkant.

Voor dit voorbeeld hoeft niet eerst een Joukowski transformatie toegepast worden.In Figuur i zien we het resultaat.De punten die genornen zijn orn te controleren of het op een cirkel wordt afgebeeld liggen precies tussen twee punten van het

vierkant.Dit voorbeeld is genomen orn het proces grafisch te bekijken en dan te vergelijken rnet dat van Theodorsen (*) die ook dit voorbeeld heeft genomen.In figuur 2 hebben we 4 als

functie van en als functie van 0 en .Vanwege de

symmetrie wordt alleen het interval [O,it/4] bekeken.Vergelijken met figuur 3,dat overgenomen is van Theodorsen,zien we dat de

overeenkomsten zeer goed zijn. De straalbuis.

Voor dit voorbeeld moet eerst de Joukowski transformatie

toegepast worden orndat we te rnaken hebben met een slank contour (zie figuur 4a).Voor de constante a hebben we de waarde 2.1 genornen.Gernakshalve geven we in tabel i de y- en z-coordinaten die we hebben gebruikt orn de de straalbuis te testen. _Het

resultaat is te zien in figuur 4b.

Tenslotte hebben we een willekeurige spanvorm bekeken. Hierbij treedt een ander probleem op.De Joukowski transforrnatie geldt alleen voor contouren die slank zijn in de lengte richting. Voor deze spanvorm geldt dit niet en daarom rnoet de spanvorrn eerst 90 graden gedraaid worden alvorens we de Joukowski

transformatie kunnen toepassen.In Figuur 5 is dat alvast gedaan en daarrnee wordt er verder gerekend en de resultaten zijn goed. Al met al kunnen we zeggen dat het _{programma goed werkt.}

(*) Referentie:T. _{Theodorsen en I.E. Garrick:general potential} theory of arbitrary wing sections.N.A.C.A. Report _{no. 452.}

-24-n Cn Ci Ci -r Ci Ci T T n

-

: = F

-aci -ci.a a.on U .40

i I

0.i3 1.20

fig.1. Transformatie van een vierkant in een cirkel.

-25--U.4U 0.0Cl û.40

Ci

11

8

rr

8

FIG. 2 FROCES TOEGEPAST 0M EEN VIERKANT TE TRANSFORMEREN IN EEN CIRKEL

-26-.3 .2 Jo .05 0 .1 »2 .3 .4 .5 .6 7 .8 Araument(8,p1,ç2,.3,p in rcdians) 4

Fig.3. Proces toegepast orn een vierkant te transformeren in een cirkel overgenomen uit .A.C.A. Report no. 452.

-27-(a) --'J - () -. (b)

J.(eJ

¡A

,/ AI '»1 j,/ .

"AI' '\

J,; '\ x = J9()J O - jb./e(3)J

/

xi /

//

/

//t&/

ìyI

Tr 3-r c-r .2 .10 -.10 -.20 40 30 .20 .10 O

't t, 't i_1LQ Lß) L.ß4 M L2 LacL Th cl i2 JL

STRAALBLJIS No 1155 [ILL 05 5K MOOEtSCIIAAL OPOI? No P50158 1'llJOELSUfN ¿LIII 111mm

Tabel 1. De y- en z-coordinaten van de straalbuis. -29-y-coordinaten. _{z-coordinaten.} 4.7 -1.05 -2.4 .0 .81795 -.924 4.8 -1.3 -2.2 .12 .8032 -.946 4.9 -1.6 -2. .22 .7855 -.964 5. -1.9 -1.6 .38 .7678 -.976 5.1 -2.15 -1.1 _{.53 .75305 -.988} 5.175 -2.5 -.75 .705 .7324 -1. 5.2 -2.8 -.4 .833 .7147 -1. 5.175 -3.25 -.1 .961 .687 -1. 5.102 -3.7 .2 _{1.069 .6616 -1.} 4.994 -4.1 .5 1.142 .638 -1. 4.866 -4.283 .7 1.167 .6181 -1. 4.65 -4.466 .9 1.1540 .616 -1. 4.5 -4.594 1.2 1.1454 .59 -1. 4.2 -4.702 1.6 1.1277 .518 -1. 3.85 -4.775 2. _{1.1071 .41 -.986} 3.55 -4.8 2.4 1.1089 .282 -.97 3.2 -4.775 2.6 1.0687 .154 -.94 2.9 -4.7 2.8 1.051 .0 _-.919 2.55 -4.6 3. _{1.0304 -.12 -.88} 2.3 -4.5 3.2 1.0156 -.208 -.84 2. -4.4 3.4 .9979 -.28 -.79 1.7 -4.3 3.6 .9802 -.356 -.73 1.45 -4.2 3.8 .9655 -.42 -.66 1.2 -4.1 4. .95 -.476 -.58 .85 -3.9 4.2 .93 -.58 -.48 .7 -3.8 4.3 .9212 -.62 -.389 .5 -3.6 4.4 .9094 -.7 -.3 .2 -3.4 4.5 _{.8912 -.764 -.23} -.1 -3.2 4.6 .874 -.816 -.16 - .3 -3.. 4.7 .8615 -.86 .08 -.6 -2.8 _{.8438 -.896 .0} -.8 -2.6 .832

'N

-

-30--[ i: :1.

Y2fl3

Q Q U-Q Q

-Q C C r rl!

-í

Q

-Q -r - .1j -Q.;fJ - -31-T I _¡- i Ì -u.t] -0.10 0.?Ci

y' rr

i .iCi

fig.5. Transformatie van een willekeurige spanvorm gedraaid over 90 graden in een cirkel.

I I

-1.4Ü -1.Ut3 -Ü.iü -Q. D. 0.Q 1.CQ

6.Appendix 1:Listing van het programma.

Hieronder voigt het programma dat een gesloten contour

transformeert in een cirkei.het hoofdprogramnma is het _{programma dat} gebruikt is voor een spanvorm.

PROORAM AF]3EELD( INPUT1 OUTPUT, TAPE5=INPUT1 _{TA?E=OUTP'JT)} PARAMETER (M=89) DIMENSION Y(M±21 i FSIi(M),TET(M I PSIFPHI(M),CPSI I CERSTET(M-l13) I EPSFTET(M),PHI( PI=4*ATAN(1. 0) C

C DE Y- EN Z-COORDINATEN LOREN LIN.kSOM C

C

READ*, (Y( I) 1=1, 23) (Z( I), 1=1, 23) DO 2 I=2445 I)=-Y(_4e-I) I)=Z(4-I) C ON T i NU E DO 3 I=4&,ó7 'R I )=-Y( I-44> Z (I )=-Z CI-44) 3 CONTINUE DO 4 I=o8, 89

Y( I)=Y(90-I)

Z (I )=-Z (90-I) 4 _CONTINUE

FRINT*, / _{Y-COORCINAAT: "}

(Y(I)

_{1=1, M)}

PRINT*, ' _{Z-COORDINAAT: ', (Z(IY 11,M)}

A=. 7O791

CALL SECOND(TO)

CALL TTRANSF(Y,Z,OsA,TETA,PHI:PSI:PSIO,EPSFTET,

CPSITET1 CPSIPHI, CEPSTETI CEPSPHI)

OND(Ti)

PR1NT*, / Ti-TO:'1 Ti-TO

C C r

ONTROLEREN OF _{VIERk4NT (ETRANSFORMEERD WORDT IN EEN CIRREL.}

C _{PSI IS LO0ARIT'iE VAN DE STRAAL EN TET DE HOER IN GETRANSFORMEERDE}

VLAR

C PARi IS _{LOGARITME VAN DE STRAAL EN PAR2 DE HOER NA JOUROW5RI TRANSFORM.}

C

DO 20 I=1 M-1

Y3( I )=(Y' Ir1

)Ys:

I)) /2. Z3( I )( Z (1+1 +Z C I)) /2.

20 CONTINUE

CALL .JTRANSF'Y2, Z3, A,M1, PAR2 FARi) 2(M±2), PAR1(M)1 Z2(M2), Y2(M+2) TET(M-1, 3) CR81 CEFSPHI(M-1 3), M), P81(M), TETA( PAR2(M) ERSPAR (M-i), PHI (M-1,

3),

'M'S ''t' (M T .,.., S. _Zu L.... St I.' I M)

C

o PSIFPHI:PSI ALS FUNCTIE VAN PHI.

C EPSFTET: INTEQRAAL VAN O NAAR 2P1 DIE BEPAALD MOET WORDEN

C

EN DIE EEN FUNCTIE VAN TETA IS.

C M. 2. V. ICSEVU BEPALEN WE EPSFTET EN PSI IN DE PUNTEN PAR2.

CALL IC8EVU(TETA, EPSFTET, M CEPSTET, M-17 PAR2, EFSPAR, M-1 1ER)

DO 15 I=1M-1

TET( I)=PARZ( I

)-4-EP9PR(I)

15 CONTINUE

CALL ICSEVU(PHI.. PSI.. M CFSIPHL M1: TEL PSIFPHI.. M1: 1ER)

DO 23 1=1.. M-1

FSI1(I)=PARI(I)-PSIFPHI(I)+pSIO

23 CONTINUE

P511 (M)=PSI1 (1)

DO 18 I=1M

Y2( I )=A*EXP

_{(P811(1) )*COS(TET( I))}

Z2( I )=A*EXP (P811(I) )*SIN(TET (I))

18 CONTINUE

CALL SECOND(T2)

PRINT*, / TYD AF2EELDING: ',T2-T1

C C C

C _{TE.ENEÑ VAN DE COORDINATEN IN HET PHYSISCHE VLAV}

C C

C

CALL SCALE(Y, 18. M7 1)

CALL bCAL(Z, 11.

i>

CALL PLOTS

CALL AXIS(O. O.

'Y-AS'.-4, 18.10. Y(N1+l)7Y(M+2))

CALL AXIS(O. O.

'ZAS',4, 11.

90., Z(M+1) Z(M+2))

CALL _{AXIS(O.,A25(Z(M1)/Z(N2))?' '7-1,-lS...Y(M+1)Y(M2))}

CALL AXIS(

_{S(Y(M±i)/Y(M+2)),O.,}

_{'1-17-11.790.7Z(M1)7Z(M2))}

CALL LINE(Y, Z; M; 1; 1; 1)

r

SUSROUTINE TTRANSF(Y1 L L1 A, TETAI PHI, PSI, PSIO, EPSFTETI

i CPSITET, CPSIPHI, CEPSTET, CEPSPHI)

PARAMETER (M=89)

C

C

C

C ROUTINE DIE EEN GE8LOTEN CONTOUR TRANSFORMERT IN EEN CIRKEL M. . V. DE

C METHODE VAN THEODORSEN. DE HOEREN IN HEI PHYSISCHE VLAK MOETEN WEL C GELIURMATIG VERSPREID ZIUN.

C

C INPUT-FRAMETERS: C

C Y,Z: Y- EN Z-000RDINAAT VAN HET PHYSISCHE VLAK.

C M: AANTAL PUNTEN IN HEI PHYSISCHE VLAR OP DE CONTOUR.

C L: L0 WE HEEDEN TE MA1EN NET EEN SLANK CONTOUR EN WE DOEN EERST C EEN UOUWROWSRI TRANSFORMATIE.

C L1 WE DOEN METEEN DE THEODORSEN'B TRANSFORMATIE.

C A:CONSTANTE DIE VLEINER IS DAN EEN KWART VAN DE RARARTERISTIEKE LENGTE.

C C

C OUTPUT-PARAMETERS:

C

C TETA: HOER IN HET PHYSISCHE VLAR.

C PHI: HOER IN HET GETRANSFORMEERDE VLAK.

C PSI:DE LOQARIT1E VAN DE STRAAL VAN EEN PUNT IN HET PHYSISCHE VLAR.

C P510: CONSTANTE WAARVAN DE EXPONENT DE STRAAL VAN GETRANSFORMEERDE C VLAK VOORSTELT.

C _{EPSFTET: INTEGRAAL F(X)=FSI(X)*COTAN((X-PHI)/2) ALS FUNCTIE VAN TETA.}

(' ...,r-j.

.'

r r -r i ri T ri T i r r'ri ' T T E' t- I r- i i A bi r T A r' t-i i b i i-' -r T i I. i r r-i

. ..,.jjL

LJ.rrJ.L.i '1iZ.i Vtii'

rbi

ri'

i

vi

C _{ÇPSIPHI:CUBIC SFLINE COEFFICIENTEN VAN PSI ALS FUNCTIE VAN PHI.}

C _{CEPSTET:CUBIC SPLINE COEFFICIENTEN VAN EPSFTET ALS FUNChE VAN TETA.} C _{CEPSPHI:CUBIC SPLINE COEFFICIENTEN VAN EPSFTET ALS FUNCTIE VAN PHI.} C C C DIMENSION Y(M+2),Z(M2),TETAi(M),WK(*M),EPSFPHI(M), i EPSTET1(M),FOUT(M),CPSITET(M-i,3),CPSIPHI(M-1,3), i CEPSTET(M-113)CEPSPHI(M-113),PSIFTET(M),PSIFPHI(M), i EPSFTET(M),FHI(M),PSI<N)TETA(M) LOGICAL H

PI=4.*ATAN(1.

IF (L.EG.0) THEN C C DE UOUROWSRI TRANSFORMATIE C

CALL JTRANSF(Y, Z, A, M, TETAS PSI)

ELSE

r

C _{GEEN J0UROWSI} TRANSFORMATIE. C DO 10

_I=1..tl-1

TETA( _{I)=ATAN2(Z( I), Y( I))}

IF (TETA( I) . LT. _{0.0) TETA(I)=TETAI)2.} *PI lo CONTINUE TETA(M)=2. *pi 1, C

C BEPAL1NG VAN PSI IN DE PUNTEN TETA.

C C DO 30 i=11M PSI( I )=ALOQ(SGRT(Y(I)**2+Z(I)**2)) 30 CONTINUE ENflE C C

C _{DE HOEEN TETA EN PHI LOPEN LINKSOM VAN O NAAR 2*PI.} C _{WE BEGINNEN MET PHI QELIK AAN TETA.}

C C C DO 15

I=LM

PHI (I )=TETA( I) 15 CONTINUE L. C C

r' T"rA r r rr _T

,.r cn:zci i

iA.or r

_{r cc rz}

Z'I rc i TcurM

L. I I I-1 I 4 I Y ' _I __ - __ y g \' y j #. Il ' I. I . t Y I1 4

C ZICH TE MIDDEN VAN TWEE _{SINGULIERE PUNTEN.}

C C C C PsIo=o. DO 20 1=2, M PSIO=PSIO±(PSI(I)PSI(I-1))*(TETA(I)-TETA(I-1))/2. 20 _CONTINUE PSIO=PSIO/(2. *PI) DO 11 1=1, M-1

TETA1 (I )=(PHI (1±1 )+FHI(I) )/2.

11 CONTI N U E

r

C

C _{BEPALING VAN P810} C

C C

C EPSTET1=INTEQRAAL ALS FUNCTIE VAN TETA.

C EPSFPHI:NIEUWE INTEGRAAL ALS FUNChE VAN DE NIEULJE PHI. C HET VERECHIL TUSSEN EPSTET1 EN EPSFTET MOET NAAR NUL TOE. C

P' Tp'i-' r' i i i. r' r-r' A! rI-. I A . I P' I riTp' P' r' i T Is Ir- r'r' r-rr T t' Tr-Isrr r-pI - r--- r' r- r' r- '-i T

L i L . r Li 'I. 1 C. r r' L. i.I '4 V r' I '4 L '.1 A i L r L i t 'IC. L Li C. r r i L J. C. '1 C.t'I I' L C. r C.r i U L) C. .

i.

C PSIFTET: PSI ALS FUNCTIE VAN TETA. C

C M. 13. V. ICSEVU BEPALEN i-JE PSI IN DE PUNTEN TElAI.

C C

L,

CALL ICSPLN(TETA, PSI, M, CPSITET, M-1g WR, 1ER)

CALL ICSEVU(TETA,PSLM,CPSITET,M-1,TETA1,PSIFTET,M-1, 1ER)

C C

C 3EPALINC VAN INTEGRAAL VOOR DE EERSTE KEER IN DE PUNTEN TETA.

C C

CALL SECOND(T5)

CALL INTOR(?HI TETA1, PSIFTET, PSI EPSFTET) CALL SECOND(T6)

PRINT*, ' TYD INTQR: ',T-T5

C C

C DEPALING VAN NIEUWE PHI. L

C

DO 32 I=l,M

PHI (I )=TETA( I )+EPSFTET( I)

32 CONTINUE

PRINT*, / PHI: '. (PHI (I), 1=1, M)

PRINT*, / EPSFTET: ', (EPSFTET(I), 1=1 M)

C C

P' ri r- r' AI r i I A k I . i rr- i Ir- i i, r-r-i ,.; p' r'i lIs IT- r-L LC. r'LJ.IM Vr'I'4 IM J. C..JWC. jIM I C.Vr'L.r 'JIM I C.L

C C

DO 35 1=1, M-1

TETA1 (I )=(PHI (1+1 )+PHI (I) )/2.

35 CONTINUE

C C C

C SE?ALING VAN PSI IN DE PUNTEN TElAi

C PSIFPHI:PSI ALS FUNCTIE VAN PHI.

C C C

CALL ICSPLN(PHI,PSI,M,CPSIPHI,M-1,WK,IER)

C

C EPSTET2 NOG MAXIMAAL L+2 MEER 3EPALEN.

C C C DO 55 N=1,L+2 DO 48 I=l,M EPSTET1 (I )=EPSFTET( I) 48 CONTINUE C C

C _{SEPALING VAN INTECP.AAL IN DE PUNTEN PHI.}

L

C

CALL INTCR(PHI, TETAIJ PSIFPHI PSI, EPSFTET)

C

5EPALIN VAN NIEUWE PHI

C L DO 102 I=1M rHI I)=TETA(I)+EPSFTET(I) 102 CONTINUE C

C BEPALING VAN DE FOUT IN DE INTEGRAAL EN INTEORAAL ALS _FUNCTIE

C VAN NIEUWE PHI.

C C

DO 135 I=11M

EPSFPHI (I )=EPSFTET( I>

FOUT( I )=EFSFTET( I)-EPSTET1( I)

IF (ABS(FOUT(I)) .QT. lE-OS) H=FALSE.

135 CONTINUE

C

C EPALINC VAN NIEUWE INTERVALBPUNTEN. C

C

DO 125 I=1M-1

TElAi (I) (PHI (1+1 )+PHI CI)) /2.

125 CONTINUE

C C C

C BEPALING VAN PSI IN DE PUNTEN TETA1 C

C C

CALL

ICSPLN(PHI,PSIM,CFSIPHIM-1,W,

1ER)

CALL ICSEVU(PHI,P9I1,CPSIPHI,M-1,TETA1,PSIFPHijM-1,IER) IF (H) GO TO 145

55 CONTINUE

C C

C CUBIC SPLINE INTERPOLATIE MET PERIODE 2P1 C

C

145 CALL ICSPLN(PHI, EPSFPHI1 M CEPSPHI, M-1 WR, 1ER) CALL ICSPLN(TETA, EPSFTET M CE?STETJ M-11 WV 1ER)

RETURN

SUBROUTINE INTQR(PHI TETA1, PSIFTET, PSI. EPSTET2)

PARAMETER (M89)

C C

C

o _{SUBROUTINE DIE INTEQRAAL VAN F(X)=PSI(X)*COTAN((XPHI)/2) BEFAALD} C IN DE PUNTEN PHI. ER WORDT GEINTEGREERD OVER _{INTERVALLEN MET}

o ALS HOEkPUNTEN TETA2 ZODANIG DAT PHI MIDDEN IN HET GEHELE INTERVAL LIG C

C

L

DIMENSION TETA1(M)PHI(M)EPSTET2(M),PSIFTET(M),B(O:23)

1 COEF(O:23),FAC(O:3O),TETA2(M)HPSIFT(M)

i P0(11i ) SOM(M-1 ), HSOMCM-1 ), DEL(M-1 ).. PSI (il)

P1=4. *ATAN( 1.

C

C C

C SEPALINO VAN FACULTEIT 2K

C C C DO i k=O,23 FAC (K)=1 DO 2 I=i;2*K FAC (k)=I*FAC (k) 2 CONTINUE i CONTINUE C C C

C _{SEPALING VAN DE COEFFICIENTEN NODIG VOOR DE 3EESENING VAN DE} C INTEGRAAL: F(X)=X*COTAN(X) C C C DO 3 RO,23 COEF(K)=(-1)**K*(2**(2*K)*B(k))/Ni2*kHFAC(K)) 3 CONTINUE C C

1

3. 333333333EO2 7. 57575757.EO2..

2. 531 135531E01:

1

1. 167E+OO, 7. O921583E0O, 5. 4971 17794EO1,

1

5. 291242424E+02 . 1921231SSE+O3 S. 65SO25311E±04

1

1. 425517167E+06,-2. 7298231O7E+O76. O158O3739E+O8

1

1.

1

4. 883323190E+14, 1. 92957934E+16, 8. 416930476E+17;

DO 10 I=11M-1

C C

C DEPALINO VAN INTERVALSPUNTEN ZODANTO DAT PHI N1IDDEN IN HEI INTERVAL

C LIOT.

C

IF (I-(M+1)/2) 11, 12.. 13

C C

C PHI LICT IN DE EERSTE (M-1)/2 INTERV..LLEN. C C 11 DO 21

J1

(M+1)/2-I TETA2(U)=TETA1(J(M-3)/2I)-2.*PI rn 'r r - n C' t r r

rr

i r- ¡ - t

nr

S

£ rl

(J

i -r

i r .i J-r ( ¿'- -r . 21 CONTINUE DO 22 J=(M+1)/2-I+1..rl TETA2(J)=TETA1(J-(M1)/2I) HPSIFT(J)=PSIFTET(J-(N-i-1)/2I) 22 CONTINUE 00 TO 14 C C

C PHI LIQT IN (M+1)/2 INTERVAL. C C 12 DO 23 J=1J1-1 TETA2(J)=TETA1 (J) HPSIFT(-J)=PSIFTET(J) 23 CONTINUE TETA2(M>=TETM. (1 )+2. *PI HPSIFTUI)=?SIFTET 1) 00 TO 14 C C C

C PHI LIOT IN DE LAATSTE (M-1)/2 INTERVALLEN.

r

C C 13 DO 25 J=1.. (3*M-1)/2-I TETA2(J)=TETA1(J+I-(M1)/2) HPSIFT(J)=PSIFTET(J+I-(M+i)/2) 25 CONTINUE DO 26 J=(3*M-1)/2-I+1,M TETA2(J)=TETA1(J-(3-M-1)/2+I)+2.*Pt HPSIFT(J)=PSIFTET(-J-(3*M-1)/2I) 26 CONTINUE 14 HULP1=0. HULP2=0.

C C

C PG(I): INTEGRAAL VAN G(X)=COTAN( (X-PHI)/2) GEINTEGREERD OVER

C INTERVAL (TETA2(I)TETA2(I+1).

C SOM(I): INTEGRAAL VAN H(X)(X-PHI)*COTAN((X-PHI)/2) GEINTEGREERD C OVER INTERVAL (TETA2(I)1TETA2(I+1).

C C DQ 20 J=1,M-1 P0(J)=ALOG(AJ38(SIN( (TETA2(J+1)-PHI(I) )/2. )/ i SIN((TETA2(J)-PHI(I))/2. ))) DEL(U)=(H.PSIFTJ+1)-HPSIFT(.J))/(TETA2(.J±1)-TETA2(.J)) SOM (3) =0. DO 30 K=023 HSOM(J)=CDEF()*(((TETA2(J+1)-?HI(I))/2.)**2k-.1)_ i ((TETA2(3)-PHI(I))12.)**(2*K+1)) 8Ot1(3)=SOM(J)HSON(U) 30 CONTINUE 20 CONTINUE C C

C INTEORAAL TOT AAN INTERVAL MET SINGULIERE PUNT

C C C DO 40 3=1, (M-3)/2 HULP1=HULP12.*HPSIFT(J)*PG(J)DEL(J)*(2.*(PHI(I)-TETA2(j)) i *PG(J)4.*SOM(J)) 40 CONTINUE C L

C INTEGRAAL NA INTERVAL MET SINGULIERE PUNT. C C DO 50 .J=(M1)/21M-1 HULP2=HULP2+2.*HPSIFT(W*PG(.J)DEL(J)*(2.*(pHI(I)_TETA2(.j)) i .50 CONTINUE C L

C INTEGRAAL INTERVAL MET SINGULIERE PUNT.

C

HULP3=2.PSI(I)?G((M-1)/2)+4.*DEL((M-1),'2)*5oM((M-1),'2)

C C C

C _{BEPALING VAN EPSTET2 EN PHI IN HET I-DE PUNT.} C

C C

EP5TET2( I)=-(HULP1+HULP2+HULP3)/(2. *PI)

10 CONTINUE

EPSTET2(M)=EPSTET2( 1) RETURN

SUBROUTINE JTRANSF(X Y A M1 TETA, PSI)

C

C

C

C SUBROUTINE DIE EEN JOUKOWSRI TRANSFORMATIE TOEPAST ALS WE TE C MA.EN HEBBEN MET SLANKE CONTOUREN.

C DE TRANSFORMATIE IS W=Z(A**2)/Z C WAARBIJ W HET PHYSISCHE VL.AK IS EN

C Z QETRANSFORMEERDE

VLA.

r'

L

C

DIMENSION X(M), Y(M)1 TETA(M),PSI(M)

PI=4*ATANl.

DO lo I=l,M

C C

C BEPALINO VAN DE HOEK. C C WORT=SORT((((XI)!A)**2+(YU)/AH*2-4)**2)/1.(y(I)/A)**2) SOMl=((X(I)/A)**2+(Y(I)/A)**2_4)/4.WORT TETA(I)=ASIN(SQRT(SOM1/2. )) IF (X(I)) 2O3O,3O

30 IF (Y(I>. LT. 0. ) TETA(I)=-TETA(I)+2. *PI

GO TO 0 20 IF (Y(I)) 40,50,50 40 TETA(I)=TETA(I)±PI GO TO 50 TETA(I)=-TETA(I)PI C C

C BEPALIN VAN DE LOGARITME VAN DE STRAAL.

C C SOM2=((X(I)/A)**2(Y(I)/A)**2-4)/4. +WORT PSI(I)=ALOG(SGRT(SOM2!2. )+SGRT(50M2/2.l)) 10 CONTINUE RETURN END