scieee Science in your language
[es] (orig)

Desarrollo e implementación en Matlab de un código numérico para la resolución de problemas de Aerodinámica potencial linealizada no estacionaria

Abstract

En el presente trabajo se desarrollado e implementado un método numérico para la resolución de problemas de aerodinámica potencial linealizada en régimen estacionario y no estacionario para alas de cualquier geometría. Dicho método consiste en una reformulación no encontrada antes en la literatura de la ecuación integral que se obtiene tras la imposición de la condición de contorno de impenetrabilidad. Dicho método ha sido implementado por el autor de este trabajo en Matlab basándose en la formulación, ecuación integral y en el programa en C originalmente desarrollados por el Profesor José Manuel Gordillo para la resolución de alas rectangulares en régimen estacionario. Posteriormente y ante la validación de los resultados obtenidos con aquellos encontrados en la literatura para otros métodos, se ha realizado la extensión a alas con geometría en flecha y también para el análisis de flujo no estacionario. El resultado fundamental de este trabajo es que se ha puesto a punto un código numérico basado en una formulación integral robusta y novedosa que permite resolver de manera muy eficiente las ecuaciones de la aerodinámica potencial subsónica en régimen no estacionario para alas de geometría genérica.

Read accessible full text

Desarrollo e implementación en Matlab de un código numérico para la resolución de problemas de Aerodinámica potencial linealizada no estacionaria

Author: Moreno Pino, Fernando
Year: 2019
Source: https://idus.us.es/bitstreams/b0b0329e-2916-4fce-94d9-28505a156d5f/download
i
Dep. Ingenie ía Ae oespacial y Mecánica de Fluidos
Escuela Técnica Supe io de Ingenie ía
Uni e sidad de Se illa
Se illa, 2019
Desa ollo e implemen ación en Ma lab de
un código numé ico pa a la esolución de
p oblemas de Ae odinámica po encial
linealizada no es aciona ia
Au o : Fe nando Mo eno Pino
Tu o : José Manuel Go dillo A ias de Saa ed a
T abajo Fin de Mas e
Mas e Uni e si a io en Ingenie ía Ae onáu ica
ii
iii
T abajo Fin de Mas e
Mas e Uni e si a io en Ingenie ía Ae onáu ica
Desa ollo e implemen ación en Ma lab
de un código numé ico pa a la esolución
de p oblemas de Ae odinámica po encial
linealizada no es aciona ia
Au o :
Fe nando Mo eno Pino
Tu o :
José Manuel Go dillo A ias de Saa ed a
Dep. de Ingenie ía Ae oespacial y Mecánica de Fluidos
Escuela Técnica Supe io de Ingenie ía
Uni e sidad de Se illa
Se illa, 2019
i
T abajo Fin de Mas e : Desa ollo e implemen ación en Ma lab de un código numé ico pa a la
esolución de p oblemas de Ae odinámica po encial linealizada no es aciona ia
Au o :
Fe nando Mo eno Pino
Tu o :
José Manuel Go dillo A ias de Saa ed a
El ibunal nomb ado pa a juzga el P oyec o a iba indicado, compues o po los siguien es miemb os:
P esiden e:
Vocales:
Sec e a io:
Acue dan o o ga le la cali icación de:
Se illa, 2019
El Sec e a io del T ibunal

i
ii
A mi amilia
A mis maes os
A mis amigos
iii
ix
Ag adecimien os
En p ime luga quie o da las g acias a mi pad e, a mi mad e y a mi he mano, pues han sido la
base en la que apoya me desde mucho an es que es udia á el G ado en Ingenie ía Ae oespacial
y po supues o du an e el Mas e en Ingenie ía Ae onáu ica.
Da g acias ambién a mi compañe a de iaje, Bea, que desde que la conozco ha es ado ahí pa a
apoya me, calma me, y saca me una son isa cada ez que lo necesi aba y sin es ue zo apa en e.
Muchas g acias po odo
Quie o ag adece ambién a odas las pe sonas que me han demos ado que me quie en y sé
que puedo con a con ellas pa a cualquie cosa: Guille, Viole a, Dani, Elena, Ál a o, F an, Alex,
Pablo… la lis a pod ía con a con una g an can idad de nomb es de pe sonas cuya o ma de se
consigue que siemp e piense que odo a a sali bien. G acias po hace me eí y po habe
compa ido conmigo an os momen os. Es un hono o ma pa e del g upo “Ae oma ia” con
odos us edes.
Da las g acias a las pe sonas que me espe aban en el pueblo cada ez que he uel o: mis abuelas
Gonzala y An onia, mi abuelo “Juani”, mi abuelo Paco, que en paz descanse, y al es o de mi
amilia.
A mis amigos, los “emp esa ios” po que aunque es én locos y sean unos payasos, siemp e se án
mis payasos.
Finalmen e, ag adece a mi u o y acompañan e en es e p oyec o, José Manuel Go dillo, pues
con él empezó mi andadu a en el mundo de la ae odinámica y con él cie o mi e apa uni e si a ia,
con un p oyec o que a a dicho campo. G acias po habe es ado ahí odos y cada uno de los
días que he enido dudas y me he ace cado a consul a le y po el a o que ha enido conmigo.
Ha sido un place .
x i
10. Mon aje............................................................................................................................... 110

x ii
ÍNDICE DE FIGURAS
Figu a 1: Geome ía Típica de un Ala en Flecha ..................................................................................... 23
Figu a 2: Supe icies de In eg ación....................................................................................................... 37
Figu a 3: No males en el en o no de la Singula idad ............................................................................. 38
Figu a 4: Disc e ización de la geome ía de Ala y Es ela ......................................................................... 45
Figu a 5: Mallado iangula simé ico en Ala y Es ela ........................................................................... 45
Figu a 6: Cambio de a iable pa a la in eg ación. G á ico explica i o I ................................................... 47
Figu a 7: Ángulos empleados en la In eg ación...................................................................................... 47
Figu a 8: Ángulos empleados en la disc e ización. G á ico explica i o II................................................. 51
Figu a 9: De inición de 𝑟(𝜃) .................................................................................................................. 51
Figu a 10: Mallado en el en o no de la singula idad .............................................................................. 52
Figu a 11: Po encial en el bo de de salida y en la es ela en égimen es aciona io pa a un alo de y dado
............................................................................................................................................................. 56
Figu a 12: Mallado iangula en égimen es aciona io .......................................................................... 56
Figu a 13: Condiciones de con o no pa a el Po encial en égimen es aciona io ..................................... 57
Figu a 14: No ación Ma icial Usada...................................................................................................... 58
Figu a 15: Condición de Ku a pa a el caso no es aciona io. .................................................................. 61
Figu a 16: P opagación de los po enciales en la Es ela no es aciona ia .................................................. 63
Figu a 17: Disc e ización en alas ec angula es ..................................................................................... 65
Figu a 18: Po encial sob e el ex adós de alas ec angula es ................................................................. 66
Figu a 19: Velocidad pe u bada en el eje x sob e el ex adós de alas ec angula es. ............................ 66
Figu a 20: Velocidad pe u bada en el eje x sob e el ex adós de alas ec angula es. Re inamien o ...... 67
Figu a 21: Disc e ización en alas con lecha .......................................................................................... 68
Figu a 22: Po encial sob e el ex adós de alas con lecha ...................................................................... 69
Figu a 23: Velocidad pe u bada en el eje x sob e el ex adós de alas con Flecha. ................................. 69
Figu a 24: E olución de la Pendien e de la Cu a de Sus en ación con el iempo. Ala Rec angula ......... 71
igu a 25: E olución de la pendien e de la cu a de sus en ación con el iempo. ala con lecha ............. 72
Figu a 26: Es ela no Es aciona ia en caso de Ala en Flecha .................................................................... 73
Figu a 27; E olución de la Pendien e de la Cu a de Sus en ación con el Ala gamien o. Ala Rec angula 77
Figu a 28: E olución de la Pendien e de la Cu a de Sus en ación con el Ala gamien o. Ψ=30º ......... 77
Figu a 29: E olución de la Pendien e de la Cu a de Sus en ación con el Ala gamien o.Ψ=45º .......... 78
Figu a 30: Mo imien o del Fluido sob e el ex adós de ae ona es con lecha ........................................ 79
Figu a 31: Velocidad pe u bada en el eje x sob e el ex adós de alas con Flecha Nega i a ................... 79
Figu a 32: Respues a an e un escalón ................................................................................................... 81
Figu a 33: E olución del Coe icien e de Sus en ación en un mo imien o oscila o io .............................. 82
Figu a 34: E olución del coe icien e de sus en ación en un mo imien o oscila o io. (KyP) [2] ............... 82
x iii
xix
ÍNDICE DE TABLAS
Tabla 1; Resul ados pa a Ala Rec angula en Régimen Es aciona io ......................................... 67
Tabla 2: Resul ados pa a Ala Rec angula en Régimen Es aciona io. Mallado Re inado ........... 67
Tabla 3: Resul ados pa a ala con Flecha 30º en égimen es aciona io. .................................... 69
Tabla 4: Resul ados pa a Ala Rec angula en Régimen No Es aciona io ................................... 71
Tabla 5: Resul ados pa a Ala con Flecha 45º en Régimen Es aciona io .................................... 72
Tabla 6: Resul ados ob enidos po dis in os mé odos pa a alas ec angula es con dis in o ala gamien o
............................................................................................................................................... 75
Tabla 7: Resul ados ob enidos po dis in os mé odos pa a alas con Ψ= 30º con dis in o ala gamien o
............................................................................................................................................... 75
Tabla 8: Resul ados ob enidos po dis in os mé odos pa a alas con Ψ= 45º con dis in o ala gamien o
............................................................................................................................................... 75
Tabla 9: E o Come ido en Alas Rec angula es en Régimen Es aciona io ................................ 80
Tabla 10: E o Come ido en Alas con Ψ=30º en Régimen Es aciona io ................................ 80
Tabla 10: E o Come ido en Alas con Ψ=45º en Régimen Es aciona io ................................ 80
xx
21
1. INTRODUCCIÓN
Desde la in ención de la a iación y has a nues os días, la ecnología empleada en es e campo
ha e olucionado conside ablemen e. Si se compa an los a iones de hoy en día con sus ances os
de la segunda gue a mundial o incluso con los an e io es biplanos, se obse a que el a ance es
ab umado . Los campos en los que más se no a es a di e encia a simple is a son la a iónica que
lle an a bo do las ae ona es ac uales, la capacidad que son capaces de anspo a , la
p opulsión, e c…
Sin emba go y aunque a simple is a no se obse e, a odas las ae ona es se les pide 2 cosas
undamen almen e: que uelen y que aguan en las ca gas a las que se en some idas.
La disciplina que se enca ga de la segunda pe ición que se les hace a las ae ona es es la
elas icidad y la esis encia de ma e iales, la cual es udia el compo amien o de la es uc u a de
la ae ona e en e a la acción de las dis in as ue zas a las que se e some ida du an e el uelo
y es ablece las ca gas lími e que es a es capaz de aguan a sin pone en iesgo a la segu idad de
los pasaje os, ipulación y es uc u a.
Po o o lado, la disciplina que se enca ga de que una ae ona e uele es la mecánica de luidos
y más especí icamen e la ae odinámica. Es a disciplina se enca ga de es udia cómo es el
mo imien o del luido al ededo de la ae ona e. Es o pe mi e modela cómo es el campo de
p esiones sob e la ae ona e y po an o conoce las ue zas que iene que sopo a la es uc u a.
Es as dos disciplinas son la base pa a lo que undamen almen e se les pide a las ae ona es. Sin
emba go ac ualmen e se ha conseguido e de o ma an na u al el hecho de que una ae ona e
despega y a e iza, que la p eocupación de las compañías aé eas ecae en inco po a y
desa olla sis emas que consigan una mayo e iciencia del uelo.
En consecuencia, la ae oelas icidad, que es la disciplina que es udia la in e acción an o de la
pa e elás ica como de la ae odinámica, es una disciplina undamen al en el diseño de cualquie
ipo de ae ona es.
Aquí en an en juego dos concep os: e icacia y e iciencia.
La e icacia se de ine cómo la capacidad de log a el e ec o que se desea, mien as que la
e iciencia implica hace lo con el mínimo desap o echamien o de los ecu sos disponibles. Una
ez que la a iación ya se consiguió que uese e icaz, aho a se desa ollan concep os y se in es iga
pa a mejo a su e iciencia. Si se analiza la e iciencia de un uelo de un pun o a o o, se puede
obse a que el uelo más e icien e es el que se desa olla en égimen de c uce o y a se posible
a un alo de coe icien e de sus en ación cons an e (c uise climb).

22
Es e abajo iene como obje i o desa olla un mé odo de cálculo que pe mi a ob ene alo es
eales del coe icien e de sus en ación de alas pa a dis in a geome ía y compa a los con los
esul ados expe imen ales de la li e a u a de o ma que se puedan ex ae conclusiones ace ca
de las en ajas e incon enien es de las di e en es geome ías. Además, en es e abajo no sólo
se a a es udia el mo imien o en égimen es aciona io, sino ambién el mo imien o de la
ae ona e en égimen no es aciona io.
Pa a ello en los dis in os pun os del abajo se an a a a los siguien es emas:
- En p ime luga se a a lle a a cabo la deducción de las ecuaciones del p oblema
ae odinámico y las simpli icaciones e hipó esis necesa ias pa a plan ea la ecuación a
esol e con el mé odo que se a a plan ea .
- A con inuación se p esen a á la disc e ización del p oblema y esolución numé ica. Se
explica án los aspec os más impo an es a ene en cuen a pa a esol e los p oblemas
numé icamen e.
- T as es o, se exponen los esul ados ob enidos pa a alas de dis in a geome ía an o
pa a el caso es aciona io como pa a el caso no es aciona io y se compa a án con algunos
ob enidos de la li e a u a así como o os ob enidos haciendo uso de o o so wa e
come cial.
- Finalmen e se expond án las conclusiones ob enidas del abajo y se p opond án
p oyec os de in es igación u u os que pe mi an mejo a lo y complemen a lo.
23
2. GEOMETRÍA
En es e pun o se an a analiza las ca ac e ís icas de la geome ía sob e la que se a a ealiza
el es udio. Dicha geome ía consis e p incipalmen e en un ala con o ma en plan a ec angula ,
pe o pos e io men e se ealiza la ex ensión a alas con geome ía apezoidal gené ica.
Se conside a á que el ala iene su o ma en plan a con enida en el plano Z = 0.
Las ca ac e ís icas p incipales que de inen a es e ipo de alas son:
- La en e gadu a (b): De inida como la dis ancia en e los bo des ma ginales del ala.
- La cue da en la aíz (c ): De inida como la dis ancia en e el bo de de a aque y el bo de
de salida del ala medida en la aíz de la misma. En nues o caso es la cue da del ala en
el plano Y=0.
- Supe icie ala (S): Supe icie de la o ma en plan a del ala, de inida como la p oyección
del ala sob e el plano Z=0.
- c : es la cue da en las pun as.
- 𝛹: es el ángulo de lecha del ala medido desde el bo de de a aque.
- α: es el ángulo de a aque.
Se ecue da que pa a alas con lecha, la geome ía iene dada po :
𝑏=√AR·𝑆 ; 𝑐𝑟=2 √𝑆 /AR
1+𝐸 ; 𝑐𝑡=𝐸𝑐𝑟
FIGURA 1: GEOMETRÍA TÍPICA DE UN ALA EN FLECHA
24
𝑥𝑎(𝑦)=|𝑦|∙ an(𝛹) 𝑥𝑠(𝑦)=𝑐𝑟+𝑏2 an(Ψ)+𝑐𝑡−𝑐𝑟
𝑏2|𝑦|
Donde se de inen los pa áme os ala gamien o y es echamien o como sigue:
AR=𝑏2
𝑆 𝐸=𝑐𝑟
𝑐𝑡
Y se ha u ilizado la siguien e nomencla u a:
- xa: de ine los pun os x del bo de de a aque del ala.
- xs: de ine los pun os x del bo de de salida del ala.
La Fo ma en plan a del ala se puede desc ibi de o ma gene al con la siguien e unción:
𝐹(𝑥,𝑦,𝑧,𝑡)≡𝑍w(𝑥,𝑦,𝑡)−𝑧=0
El ala es a á compues a po pe iles ae odinámicos esbel os, lo cual se á necesa io pa a la
pos e io simpli icación y linealización de las ecuaciones gene ales de la mecánica de luidos,
con el in de ob ene unas ecuaciones más sencillas y esolubles. La condición de esbel ez de
los pe iles ae odinámicos implica que las a iaciones de la geome ía a lo la go de las
di ecciones longi udinal y ans e sal del ala son pequeñas. Ma emá icamen e es o es:
𝜕𝑍𝑤
𝜕𝑥~∆𝑍𝑐
∆𝑥𝑐~𝑧0
𝑐𝑟≪1
𝜕𝑍𝑤
𝜕𝑦~∆𝑍𝑐
∆𝑦𝑐~𝑧0
𝑏≪1
Siendo z0 un alo de espeso ca ac e ís ico del ala.
De o ma gene al, la ecuación que de ine la geome ía del ala se puede exp esa como la suma
de 2 con ibuciones. 𝑍𝑤(𝑥,𝑦,𝑡)=𝑍𝑐(𝑥,𝑦,𝑡)±𝑍𝑒(𝑥,𝑦)
Donde Zc(x,y, ) ep esen a los pun os de las líneas de cu a u a de cada pe il en cada ins an e
de iempo y Ze(x,y) es la unción que ep esen a el espeso de cada pe il. Que el alo de la
unción de espeso de cada pe il no dependa del iempo, es indica i o de que el pe il no su i á
ninguna de o mación a lo la go del p oblema.
25
3. ECUACIONES Y CONDICIONES DE CONTORNO
3.1. Ecuaciones Gene ales
En es e pun o se plan ea án las ecuaciones que gobie nan el p oblema a pa i de las ecuaciones
gene ales de la Mecánica de Fluidos y se impond án las condiciones de con o no del p oblema.
El mo imien o del ai e al ededo del ala iene egido po las ecuaciones de Na ie -S okes, las
cuales se pueden esc ibi de o ma gene al como se mues a a con inuación.
 Ecuación de Con inuidad:
𝜕𝜌
𝜕𝑡+∇·(𝜌𝑣)=0
(3.1.1)
 Ecuación de Can idad de Mo imien o:
𝜌𝜕𝑣
𝜕𝑡+𝜌𝑣·∇𝑣=−∇𝑝+∇·𝜏′+𝜌𝑓𝑚
󰇍
󰇍
󰇍
󰇍

(3.1.2)
 Ecuación de la Ene gía:
𝜌𝑐𝑣𝜕𝑇
𝜕𝑡+𝜌𝑐𝑣𝑣·∇𝑇=−𝑝∇·𝑣+𝜏′:∇𝑣+𝑄𝑟+𝑄𝑞+∇·(𝑘∇𝑇)
(3.1.3)
Es e sis ema con o ma un conjun o de 5 ecuaciones con 6 incógni as:
- 𝜌: Densidad del luido
- 𝑣: Velocidad del luido
- P: P esión del luido
- T: Tempe a u a del luido
Pa a ce a el p oblema se debe comple a el sis ema con la ecuación de es ado del luido bajo
es udio, el cual se a a conside a que es un gas pe ec o.
 Ecuación de Es ado:
𝑝
𝜌=𝑅𝑔𝑇
(3.1.4)
Donde 𝑅𝑔=𝑅
𝑀𝑚=𝑐𝑝−𝑐𝑣=𝑐𝑣(𝛾−1). Siendo 𝑅=8,314 𝐽
𝑚𝑜𝑙·𝐾 la cons an e uni e sal de los
gases y 𝑀𝑚 la masa mola del gas.
Con es o, se iene un sis ema no lineal de ecuaciones en de i adas pa ciales con 6 ecuaciones y
6 incógni as al que hab á que p opo ciona le las condiciones de con o no e iniciales adecuadas
pa a cada p oblema.
32
ℎ=𝐶𝑝𝑇=𝐶𝑝𝑝
𝑅𝑔𝜌=(𝐶𝑝
𝐶𝑣)𝑝
(𝑅𝑔
𝐶𝑣)𝜌=𝛾
𝛾−1𝑝
𝜌=𝑎2
𝛾−1
Se iene que:
𝜕𝜙
𝜕𝑡+|∇𝜙|2
2+𝑎2
𝛾−1=𝐶(𝑡)
En es e p oyec o se conside a á únicamen e el lími e 𝑀∞≪1 lo que pe mi e simpli ica las
ecuaciones aún más. En es e lími e 1
𝜌Δ𝜌<< 1. Lo que implica que el luido puede se conside ado
como incomp esible. En es e caso la ecuación de la con inuidad se simpli ica a:
∇·𝑣′
󰇍
󰇍
󰇍

=0
Los é minos siguien e pueden se desp eciados. Véase Go dillo y Riboux, 2012 [1]
∂ρ
∂ ; 𝑣′
󰇍
󰇍
󰇍

·∇𝜌≪𝜌∇·𝑣′
󰇍
󰇍
󰇍

Pa a el caso 𝜌≈𝑐𝑡𝑒, la ecuación escala pa a el po encial queda:
𝜕𝜙
𝜕𝑡+|∇𝜙|2
2+𝑝
𝜌=𝐶(𝑡)
(3.4.2)
Se obse a que la ecuación de can idad de mo imien o ha pasado a se una ecuación escala
pa a el po encial de elocidades. 𝐶(𝑡) es una cons an e que es únicamen e una unción del
iempo. Conocido el po encial de elocidades, se puede u iliza la ecuación an e io pa a
de e mina la dis ibución de p esiones.
En el caso de 𝜌≈𝑐𝑡𝑒, la sus i ución de la elocidad en unción del po encial en la ecuación de
con inuidad p opo ciona que el po encial de elocidades sa is ace la ecuación de Laplace.
∇2𝜙=0
Con la hipó esis de lujo incomp esible se es a ían come iendo e o es de 𝑂(𝑀2), donde
𝑀=𝑣/𝑎 deno a el núme o de Mach. Dicha demos ación se puede encon a en [1]. Es a
hipó esis es álida pa a alas olando a un Mach<0.3.
Debido a que la ecuación de Laplace es lineal, se puede aplica el mé odo de supe posición de
soluciones elemen ales pa a halla la solución del p oblema gene al.
El p oblema gene al queda de inido po :
∇2𝜙=0
|𝑥|→∞, 𝜙→𝑈∞𝑥
(3.4.3)

33
|𝑥|∈Σ𝑠, 𝑛
󰇍

·(∇𝜙−𝑣𝑝
󰇍
󰇍
󰇍
󰇍

)=0
𝐶𝑜𝑛𝑑𝑖𝑐𝑖𝑜𝑛 𝑑𝑒 𝐾𝑢𝑡𝑡𝑎−𝐽𝑜𝑢𝑘𝑜𝑤𝑠𝑘𝑖
Donde 𝑣𝑝
󰇍
󰇍
󰇍
󰇍

es la elocidad e ical del pe il en caso de encon a nos en égimen no es aciona io
y la condición de Ku a-Joukowski es la que pe mi e de e mina la solución eal a nues o
p oblema de en e las in ini as soluciones exis en es que iene el p oblema del lujo po encial
al ededo de un obje o.
3.5. Ecuaciones Linealizadas
A con inuación se a a p ocede a linealiza las ecuaciones a esol e . Se a a conside a que el
ala es á o mada po pe iles muy esbel os y se mue e a ángulos de a aque pequeños. Es o
pe mi e supone que la capa lími e no se a a desp ende en ningún pun o del ala y en
consecuencia, dicho ala pe u ba poco la co ien e inciden e. Es o quie e deci , que las
pe u baciones en las a iables del luido, o iginadas po la p esencia del ala, son de un o den
de magni ud in e io a su alo aguas a iba.
Las hipó esis a conside a son:
𝛼≪1
ℎ0≪𝑐
𝜕𝑧𝑝
𝜕𝑥≪1, 𝜕𝑧𝑝
𝜕𝑦≪1
(3.5.1)
Donde ℎ0 es el espeso ca ac e ís ico del pe il. Es as hipó esis jun o con la de núme o de
Reynolds muy ele ado, nos pe mi en asegu a que la capa lími e no se desp ende y que la
pe u bación en las a iables luidas es pequeña.
Con es as hipó esis los campos de elocidades, p esión, densidad y empe a u a de las a iables
luidas se puede eesc ibi como la suma de su alo aguas a iba más una pe u bación
o iginada po la p esencia del ala (3.5.2):
𝑣(𝑥)=𝑈∞
󰇍
󰇍
󰇍
󰇍
󰇍

+𝑣′(𝑥), |𝑣′(𝑥)|≪𝑈∞
𝑝(𝑥)=𝑝∞+𝑝′(𝑥), 𝑝′(𝑥)≪𝑝∞
(3.5.2)
Donde 𝑣′(𝑥), 𝑝′(𝑥), 𝜌′(𝑥), 𝑇′(𝑥) son los campos pe u bados de elocidad, p esión, densidad y
empe a u a.
El po encial de elocidades se puede descompone igualmen e como la suma del po encial de
elocidades en el in ini o más el po encial de elocidades pe u badas (3.5.3), que se á la
incógni a de nues o p oblema.
𝜙(𝑥)=𝜙∞+𝜙′(𝑥)
(3.5.3)
34
Debido a la linealidad de la ecuación de Laplace, el po encial de elocidades pe u badas debe
de cumpli (3.5.4):
∇2𝜙′=0
(3.5.4)
La ecuación de can idad de mo imien o en é minos del po encial de elocidades pe u badas
debe de cumpli (3.5.5)
𝜌𝜕𝜙′
𝜕𝑡+𝜌𝑈∞𝜕𝜙′
𝜕𝑥+𝑝′=0
(3.5.5)
La condición de impene abilidad puede simpli ica se g acias a la linealización de las ecuaciones:
𝑛
󰇍

·(∇𝜙−𝑣𝑝
󰇍
󰇍
󰇍
󰇍

)=∇𝐹𝑒,𝑖
|∇𝐹𝑒,𝑖|·(∇𝜙−𝑣𝑝
󰇍
󰇍
󰇍
󰇍

)=0→∇𝐹𝑒,𝑖·(∇𝜙−𝑣𝑝
󰇍
󰇍
󰇍
󰇍

)
Donde 𝐹𝑒,𝑖 es la unción de los pun os del ex adós y del in adós del ala. Se puede eesc ibi
pues la ecuación de la siguien e o ma:
(𝑒3
󰇍
󰇍
󰇍

−𝜕𝑧𝑝
𝜕𝑥𝑒1
󰇍
󰇍
󰇍

−𝜕𝑧𝑝
𝜕𝑦𝑒2
󰇍
󰇍
󰇍

)·((𝜕𝜙′
𝜕𝑧−𝜕𝑧𝑝
𝜕𝑡)𝑒3
󰇍
󰇍
󰇍

+(𝑈∞+𝜕𝜙′
𝜕𝑥)𝑒1
󰇍
󰇍
󰇍

+𝜕𝜙′
𝜕𝑦𝑒2
󰇍
󰇍
󰇍

)=0
Desp eciando los é minos de segundo o den, se llega a (3.5.6):
𝜕𝜙′
𝜕𝑧(𝑥,𝑦,𝑧=0)=𝑤′(𝑥,𝑦,𝑧=0)=𝑈∞𝜕𝑧𝑝
𝜕𝑥(𝑥,𝑦,𝑡)+𝜕𝑧𝑝
𝜕𝑡(𝑥,𝑦,𝑡)
(3.5.6)
Es a condición se debe ía de impone es ic amen e sob e el ala, pe o debido a que los pe iles
son esbel os y los ángulos que amos a maneja son pequeños, dicha condición se impone sob e
la o ma en plan a del ala (FP). G acias a es as hipó esis, es a ap oximación no supone un e o
impo an e. Véase Go dillo y Riboux, 2012 [1].
El p oblema linealizado iene pues desc i o po las siguien es ecuaciones:
∇2𝜙′=0
|𝑥|→∞, 𝜙′→0
|𝑥|∈𝐹𝑃,𝜕𝜙′
𝜕𝑧(𝑥,𝑦,𝑡)=𝑈∞𝜕𝑧𝑝
𝜕𝑥(𝑥,𝑦,𝑡)+𝜕𝑧𝑝
𝜕𝑡(𝑥,𝑦,𝑡)
𝐶𝑜𝑛𝑑𝑖𝑐𝑖𝑜𝑛 𝑑𝑒 𝐾𝑢𝑡𝑡𝑎−𝐽𝑜𝑢𝑘𝑜𝑤𝑠𝑘𝑖
(3.5.7)
3.6. Condición de Ku a-Joukowky
En es a apa ado se explica la condición de Ku a-Joukowsky y sus implicaciones en el p oblema
ae odinámico.
35
Dicha condición es ablece que de las in ini as soluciones que iene el mo imien o de un luido
al ededo del ala, es aquella pa a la cual el lujo no ebo dea el bo de de salida. Es a condición
ija el alo de la ci culación (Γ), pe mi iendo elegi de en e odas las soluciones aquella que se
da en la ealidad.
La condición de con o no que se debe de impone en el bo de de salida se ob iene de es a
condición, jun o con la condición de que la di e encia de p esiones en e es ados e in adós
debe de se nula a lo la go del bo de de salida po que no exis e sólido que sopo ase dicha
di e encia en el caso de que la hubiese.
En el ex adós y en el in adós del bo de de salida se debe de cumpli que:
𝜌𝜕𝜙+
𝜕𝑡 +𝜌𝑈∞𝜕𝜙+
𝜕𝑥 +𝑝+=0
𝜌𝜕𝜙−
𝜕𝑡 +𝜌𝑈∞𝜕𝜙−
𝜕𝑥 +𝑝−=0
Si se es an ambas ecuaciones y se ienen cuen a que la p esión en el ex adós y en el in adós
son iguales: 𝜕(𝜙+−𝜙−)
𝜕𝑡 +𝑈∞𝜕(𝜙+−𝜙−)
𝜕𝑥 =0
Teniendo en cuen a la an isime ía del campo de elocidades en el p oblema sus en ado ,
Go dilllo y Riboux, 2012, [1], 𝜙+=−𝜙−, se puede eesc ibi la condición an e io pa a
cualquie a de las supe icies ex adós o in adós. Sup imiendo el supe índice + en el po encial,
se iene inalmen e.
𝜕𝜙
𝜕𝑡+𝑈∞𝜕𝜙
𝜕𝑥=𝐷𝜙
𝐷𝑡=0
(3.6.1)
Es a ecuación es la que se a a impone en el bo de de salida. Es o quie e deci que la a iación
po unidad de iempo del po encial de elocidades siguiendo a la pa ícula luida es nula.
En el caso es aciona io, es a condición se aduce en que en el bo de de salida hay un pun o de
emanso (en el caso de bo de de salida anguloso), o bien que la elocidad en el bo de de salida
iene la misma di ección y sen ido an o en el ex adós como en el in adós (caso de bo de de
salida en e oceso).
En el caso es aciona io la condición an e io queda:
𝜕𝜙
𝜕𝑥=0
(3.6.2)
36
37
4. SOLUCIÓN GENERAL DEL PROBLEMA
Se de ine el po encial 𝜓0 de una solución elemen al de la ecuación de Laplace pa a el caso
idimensional co espondien e a una uen e localizada en el pun o (𝑥0,𝑦0,𝑧0) como:
𝜓0=1
√(𝑥−𝑥0)2+(𝑦−𝑦0)2+(𝑧−𝑧0)2
(4.1)
Se puede comp oba que es e ipo de solución cumple la ecuación de Laplace ∇2𝜓0=0.
El po encial de elocidades pe u badas, debe de cumpli ambién la ecuación de Laplace
∇2𝜙′=0.
Reco demos que se iene que cumpli que:
∇2𝜙′=0, 𝑣′=∇𝜙′, → ∇2𝑣′=0
Donde se ha hecho uso de la in e cambiabilidad de las de i adas.
Debido a es o se debe de cumpli la siguien e ecuación:
𝜓0∇2𝑣′−𝑣′ ∇2𝜓0=0
(4.2)
A con inuación se p ocede a in eg al es a ecuación en un olumen de con ol Ω𝑐 limi ado po la
supe icie Σ𝑐=Σ∞∪Σ𝜖∪Σ𝑎𝑙𝑎∪Σ𝑒𝑠𝑡𝑒𝑙𝑎+∪Σ𝑒𝑠𝑡𝑒𝑙𝑎−
. Dichas supe icies se obse an en 2D en
la siguien e igu a:
∫(𝜓0∇2𝑣−𝑣 ∇2𝜓0)𝑑𝜔=
Ω𝑐∫∇·(𝜓0∇2𝑣−𝑣 ∇2𝜓0)𝑑𝜔=
Ω𝑐
FIGURA 2: SUPERFICIES DE INTEGRACIÓN

38
=∫𝜓0𝑛
󰇍

·∇𝑣′−𝑣′𝑛
󰇍

·∇𝜓0 𝑑𝜎=0
Σ𝑐
La in eg al a lo la go de Σ∞ es nula, y las in eg ales a lo la go de Σ𝜖 y Σ𝑎𝑙𝑎 pueden esc ibi se de
la siguien e o ma:
∫𝜓0𝑛
󰇍

·∇𝑣′−𝑣′𝑛
󰇍

·∇𝜓0𝑑𝜎=
Σ𝜖∪Σ𝑎𝑙𝑎
=∫𝜓0𝑛
󰇍

·∇𝑣′−𝑣′𝑛
󰇍

·∇𝜓0 𝑑𝜎+
Σ𝑎𝑙𝑎 ∫𝜓0𝑛
󰇍

𝑒−·∇𝑣′−𝑣′𝑛
󰇍

𝑒−·∇𝜓0 𝑑𝜎+
Σ𝜖−
+∫𝜓0𝑛
󰇍

·∇𝑣′−𝑣′𝑛
󰇍

·∇𝜓0 𝑑𝜎+∫𝜓0(−𝑛
󰇍

𝑒−)·∇𝑣′−𝑣′(−𝑛
󰇍

𝑒−)·∇𝜓0𝑑𝜎=
Σ𝜖−
Σ𝜖
=𝐼𝑎𝑙𝑎+𝐼𝜖
Donde Σ𝜖 es la semies e a supe io que odea la singula idad y cuya no mal posi i a se de ine
hacia el in e io de la es e a, y Σ𝜖− es la semies e a in e io que odea la singula idad y cuya
no mal posi i a se de ine hacia el ex e io de dicha es e a. La in eg al 𝐼𝜖 se puede calcula de la
siguien e mane a.
𝐼𝜖=lim
ϵ→0∫(1ϵ(−𝑒𝑟)·∇𝑣′(𝑥)−𝑣′(𝑥)(𝑒𝑟)·(−𝑒𝑟)
𝜖2)𝜖2𝑑Ω=−4𝜋𝑣′(𝑥0
󰇍
󰇍
󰇍
󰇍

)
4𝜋
0
(4.3)
Donde en es e caso 𝑑Ω=𝑑S/ 2 es el di e encial de ángulo sólido.
T as es e cálculo la ecuación de G een puede esc ibi se de la siguien e o ma:
4𝜋𝑣′(𝑥0)=∫𝜓0𝑛
󰇍

·∇𝑣′−𝑣′𝑛
󰇍

·∇𝜓0𝑑𝜎
Σ𝑎𝑙𝑎∪Σ𝑒𝑠𝑡𝑒𝑙𝑎+∪Σ𝑒𝑠𝑡𝑒𝑙𝑎−
(4.4)
P oyec ando sob e el ec o uni a io 𝑘
󰇍

.
4𝜋𝑤
󰇍
󰇍

′(𝑥0
󰇍
󰇍
󰇍
󰇍

)=∫𝜓0𝑛
󰇍

·∇𝑤
󰇍
󰇍

′−𝑤
󰇍
󰇍

′𝑛
󰇍

·∇𝜓0𝑑𝜎
Σ𝑎𝑙𝑎∪Σ𝑒𝑠𝑡𝑒𝑙𝑎+∪Σ𝑒𝑠𝑡𝑒𝑙𝑎−
(4.5)
Debido a que la an isime ía 𝑤′(𝑥,𝑦,𝑧=0+)=𝑤′(𝑥,𝑦,𝑧=0−) y a que las no males al ala y a
la es ela son ap oximadamen e 𝑘
󰇍

y -𝑘
󰇍

pa a el ex adós y pa a el in adós, se iene pues que:
∫𝑤′(𝑥,𝑦,𝑧=0 )𝑛
󰇍

·∇𝜓0𝑑𝜎=0
Σ𝑒𝑠𝑡𝑒𝑙𝑎+∪Σ𝑒𝑠𝑡𝑒𝑙𝑎−
FIGURA 3: NORMALES EN EL ENTORNO DE LA SINGULARIDAD
39
∫𝑤′(𝑥,𝑦,𝑧=0 )𝑛
󰇍

·∇𝜓0𝑑𝜎=0
Σ𝑎𝑙𝑎+∪Σ𝑎𝑙𝑎−
Luego se llega a la siguien e ecuación
4𝜋𝑤
󰇍
󰇍

′(𝑥0)=∫𝜓0𝑛
󰇍

·∇𝑤
󰇍
󰇍

′𝑑𝜎=
Σ𝑎𝑙𝑎∪Σ𝑒𝑠𝑡𝑒𝑙𝑎+∪Σ𝑒𝑠𝑡𝑒𝑙𝑎−
=−∫𝜓0𝜕𝑤′+
𝜕𝑧 𝑑𝜎++∫𝜓0𝜕𝑤′−
𝜕𝑧 𝑑𝜎−
Σ𝑎𝑙𝑎−∪Σ𝑒𝑠𝑡𝑒𝑙𝑎−
Σ𝑎𝑙𝑎+∪Σ𝑒𝑠𝑡𝑒𝑙𝑎+
Aplicando la ecuación de con inuidad an o en el ex adós como en el in adós se iene:
∇·𝑣′+
󰇍
󰇍
󰇍
󰇍
󰇍
󰇍

=𝜕𝑢′+
𝜕𝑥 +𝜕𝑣′+
𝜕𝑦 +𝜕𝑤′+
𝜕𝑧 =0 → 𝜕𝑤′+
𝜕𝑧 =−𝜕𝑢′+
𝜕𝑥 −𝜕𝑣′+
𝜕𝑦
∇·𝑣′−
󰇍
󰇍
󰇍
󰇍
󰇍
󰇍

=𝜕𝑢′−
𝜕𝑥 +𝜕𝑣′−
𝜕𝑦 +𝜕𝑤′−
𝜕𝑧 =0 → 𝜕𝑤′−
𝜕𝑧 =−𝜕𝑢′−
𝜕𝑥 −𝜕𝑣′−
𝜕𝑦
Dado que la placa pa e del eposo, la ci culación Γ es idén icamen e nula en cualquie supe icie
ce ada que englobe al ala como a la es ela según el Teo ema de Bje knes-Kel in, el cual dice:
𝐷Γ
𝐷𝑡=0
(4.6)
Tomando como supe icie, una compues a po el ala y la es ela, se iene:
Γ=∫ ∫ 𝛾(𝑥,𝑦,𝑡) 𝑑𝑥 𝑑𝑦=∫ ∫ (𝑢′+−𝑢′−) 𝑑𝑥 𝑑𝑦=∫ ∫ 𝜕(𝜙+−𝜙−)
𝜕𝑥 𝑑𝑥 𝑑𝑦
𝑥𝑒𝑠𝑡
0
𝑏
2
−𝑏
2
𝑥𝑒𝑠𝑡
0
𝑏
2
−𝑏
2
𝑥𝑒𝑠𝑡
0
𝑏
2
−𝑏
2
Γ=∫(𝜙+(𝑥𝑒𝑠𝑡,𝑦,𝑡)−𝜙−
𝑏
2
−𝑏
2(𝑥𝑒𝑠𝑡,𝑦,𝑡))𝑑𝑦=0
Donde el subíndice “es ” hace e e encia a la es ela y 𝑥𝑒𝑠𝑡 es el inal de la es ela y la a iable
𝛾(𝑥,𝑦,𝑡), es la densidad de ci culación. Se ha u ilizado que el po encial es 0 en el bo de de
a aque, condición que se explica á pos e io men e, cuando se analicen las condiciones de
con o no del p oblema.
Como es a in eg al es nula pa a odo ins an e de iempo, de ello se deduce que el in eg ando
no puede depende del iempo. Es o implica que la a iación con el iempo del po encial en la
coo denada donde acaba la es ela es nula, y po an o las elocidades ho izon ales son iguales.
U ilizando las siguien es igualdades ob enidas de la de i ación de un p oduc o:
𝜓0𝜕𝑢′
𝜕𝑥=𝜕(𝜓0𝑢′)
𝜕𝑥 −𝑢′𝜕𝜓0
𝜕𝑥
40
𝜓0𝜕𝑣′
𝜕𝑦=𝜕(𝜓0𝑣′)
𝜕𝑦 −𝑣′𝜕𝜓0
𝜕𝑦
Y sus i uyendo en la ecuación de G een, se iene que:
4𝜋𝑤
󰇍
󰇍

′(𝑥0)=∫(𝜕(𝜓0𝑢′+)
𝜕𝑥 −𝑢′+𝜕𝜓0
𝜕𝑥 +𝜕(𝜓0𝑣′+)
𝜕𝑦 −𝑣′+𝜕𝜓0
𝜕𝑦)𝑑𝜎+−
Σ𝑎𝑙𝑎+∪Σ𝑒𝑠𝑡𝑒𝑙𝑎+
−∫(𝜕(𝜓0𝑢′−)
𝜕𝑥 −𝑢′−𝜕𝜓0
𝜕𝑥+𝜕(𝜓0𝑣′−)
𝜕𝑦 −𝑣′−𝜕𝜓0
𝜕𝑦)𝑑𝜎−
Σ𝑎𝑙𝑎−∪Σ𝑒𝑠𝑡𝑒𝑙𝑎−
Pa a es e azonamien o se a a ene en cuen a que las in eg ales de supe icie se pueden hace
indis in amen e p ime o en la di ección x y luego en la di ección y o ice e sa. Tal como se ha
esc i o la exp esión an e io se obse a que los é minos son iguales en ambas in eg ales,
simplemen e cambiando la no ación de ex adós a in adós. En onces:
- La suma de las in eg ales del p ime é mino es 0, debido a que al hace la in eg al en
di ección x se iene que en el bo de de a aque dichas elocidades son 0, y a que al inal
de la es ela se cumple 𝑢′+=𝑢′− .
- La suma de las in eg ales del e ce é mino es 0 debido a que al hace la in eg al en
di ección y se debe de ene en cuen a que las elocidades 𝑣′ son an isimé icas
espec o al plano z=0, siendo en onces dichas elocidades idén icas pe o de signo
con a io en los bo des ma ginales.
Eliminando es as in eg ales de la exp esión an e io y haciendo uso de nue o de la an isime ía
del campo de elocidades se llega a:
2𝜋𝑤
󰇍
󰇍

′(𝑥0)=−∫(𝑢′+𝜕𝜓0
𝜕𝑥+𝑣′+𝜕𝜓0
𝜕𝑦)𝑑𝜎+
Σ𝑎𝑙𝑎+∪Σ𝑒𝑠𝑡𝑒𝑙𝑎+
Ecuación que inalmen e queda:
𝑤
󰇍
󰇍

′(𝑥0)=𝑈∞𝜕𝑧𝑝
𝜕𝑥(𝑥0,𝑦0,𝑡)+𝜕𝑧𝑝
𝜕𝑡(𝑥0,𝑦0,𝑡)=−1
2𝜋∫𝑣′(𝑥)·∇𝜓0 𝑑𝜎+
Σ𝑎𝑙𝑎+∪Σ𝑒𝑠𝑡𝑒𝑙𝑎+
(4.7)
Es a es la solución gene al del p oblema ae odinámico idimensional. Se obse a que dicha
exp esión se co esponde con la ecuación de Bio y Sa a . La a iable 𝑣′(𝑥) indica las
componen es de la elocidad en el plano Z=0. El lec o in e esado puede encon a más
in o mación en Go dillo y Riboux, 2012, [1].
De es a o ma se obse a que la componen e e ical del campo de elocidades pe u bado en
cualquie pun o del dominio luido puede se exp esada como la supe posición del campo de
elocidades c eado po una dis ibución con inua de o bellinos si uados sob e la supe icie del
ala y de la es ela, y o ien ados según la en e gadu a, con in ensidad po unidad de longi ud
41
2𝑢′+(𝑥,𝑦), y una dis ibución con inua de o bellinos si uados sob e ala y es ela, de in ensidad
po unidad de longi ud 2𝑣′+(𝑥,𝑦) o ien ados en la di ección del eje x.
T abajos an e io es a an sob e la esolución de es a ecuación in eg al. Sin emba go, dichos
es udios demues an que el sis ema que se ob iene es ines able desde el pun o de is a
numé ico. Véase Flo es Caballe o, 2016 [4]. El mé odo que se a a implemen a co ige es a
ines abilidad median e la ans o mación de dicha ecuación in eg al, ya que en el in eg ando
apa ece á como incógni a el alo del po encial pe u bado en luga de su g adien e.
48
De es a o ma el denominado de la ecuación se puede desa olla como sigue:
|𝑅
󰇍
󰇍
󰇍

−𝑅0
󰇍
󰇍
󰇍
󰇍

| =|(𝑅
󰇍
󰇍
󰇍

−𝑅1
󰇍
󰇍
󰇍
󰇍

)+(𝑅1
󰇍
󰇍
󰇍
󰇍

−𝑅0
󰇍
󰇍
󰇍
󰇍

)| =|𝑟𝑒𝑟
󰇍
󰇍
󰇍

+(𝑅1
󰇍
󰇍
󰇍
󰇍

−𝑅0
󰇍
󰇍
󰇍
󰇍

)| =
=|𝑟𝑒𝑟
󰇍
󰇍
󰇍

+(𝑅1
󰇍
󰇍
󰇍
󰇍

−𝑅0
󰇍
󰇍
󰇍
󰇍

)| =(𝑟2+|𝑅1
󰇍
󰇍
󰇍
󰇍

−𝑅0
󰇍
󰇍
󰇍
󰇍

|2+2𝑟𝑒𝑟
󰇍
󰇍
󰇍

·(𝑅1
󰇍
󰇍
󰇍
󰇍

−𝑅0
󰇍
󰇍
󰇍
󰇍

))1/2=
=((𝑟+𝑒𝑟
󰇍
󰇍
󰇍

·(𝑅1
󰇍
󰇍
󰇍
󰇍

−𝑅0
󰇍
󰇍
󰇍
󰇍

))2+|𝑅1
󰇍
󰇍
󰇍
󰇍

−𝑅0
󰇍
󰇍
󰇍
󰇍

|2−(𝑒𝑟
󰇍
󰇍
󰇍

·(𝑅1
󰇍
󰇍
󰇍
󰇍

−𝑅0
󰇍
󰇍
󰇍
󰇍

))2)1/2=
=(𝑈2+𝑎2)1/2
Donde x y a ienen dadas po las siguien es exp esiones (5.8) y (5.9):
𝑈=𝑟+𝑒𝑟
󰇍
󰇍
󰇍

·(𝑅1
󰇍
󰇍
󰇍
󰇍

−𝑅0
󰇍
󰇍
󰇍
󰇍

)
(5.8)
𝑎2=|𝑅1
󰇍
󰇍
󰇍
󰇍

−𝑅0
󰇍
󰇍
󰇍
󰇍

|2−(𝑒𝑟
󰇍
󰇍
󰇍

·(𝑅1
󰇍
󰇍
󰇍
󰇍

−𝑅0
󰇍
󰇍
󰇍
󰇍

))2
(5.9)
Es a nomencla u a se á u ilizada más adelan e.
U ilizando las nue as a iables 𝜃 y , se llega a la siguien e conclusión.
𝑥−𝑥1=𝑟𝑐𝑜𝑠(𝜃)
𝑦−𝑦1=𝑟𝑠𝑖𝑛(𝜃)
Sus i uyendo en la exp esión an e io , desa ollando el denominado de la in eg al y u ilizando
la de inición de las a iables U y 𝑎2 se llega a lo siguien e.
𝑤𝑖=∫𝜙(𝑥𝑖)
|𝑅𝑖
󰇍
󰇍
󰇍

−𝑅0
󰇍
󰇍
󰇍
󰇍

|3𝑑𝜎𝑖
Σi= 𝜙1∫[ 1
(𝑈2+𝑎2)3/2+𝐶1𝑟𝑐𝑜𝑠(𝜃)
(𝑈2+𝑎2)3/2+𝐷1𝑟𝑠𝑖𝑛(𝜃)
(𝑈2+𝑎2)3/2]𝑑𝜎𝑖
Σi+
+ 𝜙2∫[ 𝐶2𝑟𝑐𝑜𝑠(𝜃)
(𝑈2+𝑎2)3/2+𝐷2𝑟𝑠𝑖𝑛(𝜃)
(𝑈2+𝑎2)3/2]𝑑𝜎𝑖
Σi+
+ 𝜙3∫[ 𝐶3𝑟𝑐𝑜𝑠(𝜃)
(𝑈2+𝑎2)3/2+𝐷3𝑟𝑠𝑖𝑛(𝜃)
(𝑈2+𝑎2)3/2]𝑑𝜎𝑖
Σi
Donde U=U( ), y los lími es de la in eg al an desde =0, has a = (𝜃). Desa ollando la in eg al
doble y sacando los é minos que no dependen de ue a de la in eg al se llega a lo siguien e.
𝜙1∫ [∫ 𝑟𝑑𝑟
(𝑈2+𝑎2)3/2
𝑟(𝜃)
0+(𝐶1𝑐𝑜𝑠(𝜃)+𝐷1𝑠𝑖𝑛(𝜃))∫𝑟2𝑑𝑟
(𝑈2+𝑎2)3/2
𝑟(𝜃)
0]𝑑𝜃+
𝜃3
θ2

49
+ 𝜙2∫ [(𝐶2𝑐𝑜𝑠(𝜃)+𝐷2sin (𝜃))∫𝑟2𝑑𝑟
(𝑈2+𝑎2)3
2
𝑟(𝜃)
0]𝑑𝜃+
𝜃3
θ2
+ 𝜙3∫ [(𝐶3𝑐𝑜𝑠(𝜃)+𝐷3𝑠𝑖𝑛(𝜃))∫𝑟2𝑑𝑟
(𝑈2+𝑎2)3/2
𝑟(𝜃)
0]𝑑𝜃
𝜃3
θ2
Como se puede obse a , las in eg ales a ealiza son simila es en e ellas. Se a a pasa a
analiza las una a una indi idualmen e.
La p ime a se mues a a con inuación:
𝑆1=∫𝑟𝑑𝑟
(𝑈2+𝑎2)3/2
𝑟(𝜃)
0=∫𝑈−𝑒𝑟
󰇍
󰇍
󰇍

·(𝑅0
󰇍
󰇍
󰇍
󰇍

−𝑅1
󰇍
󰇍
󰇍
󰇍

)
(𝑈2+𝑎2)3/2 𝑑𝑈
𝑈(𝜃)
𝑈0
Donde se ha hecho el cambio de a iable mos ado an e io men e. Ope ando con es a exp esión
se ob iene lo siguien e:
∫𝑈−𝑒𝑟
󰇍
󰇍
󰇍

·(𝑅0
󰇍
󰇍
󰇍
󰇍

−𝑅1
󰇍
󰇍
󰇍
󰇍

)
(𝑈2+𝑎2)3/2 𝑑𝑈
𝑈(𝜃)
𝑈0=∫𝑈
(𝑈2+𝑎2)3/2𝑑𝑈
𝑈(𝜃)
𝑈0−𝑒𝑟
󰇍
󰇍
󰇍

·(𝑅0
󰇍
󰇍
󰇍
󰇍

−𝑅1
󰇍
󰇍
󰇍
󰇍

)∫ 𝑑𝑈
(𝑈2+𝑎2)3/2
𝑈(𝜃)
𝑈0
𝑆1=−𝑆0−𝑒𝑟
󰇍
󰇍
󰇍

·(𝑅0
󰇍
󰇍
󰇍
󰇍

−𝑅1
󰇍
󰇍
󰇍
󰇍

)𝑆3
Donde:
𝑆0=−∫𝑈
(𝑈2+𝑎2)3
2𝑑𝑈
𝑈(𝜃)
𝑈0=1
((𝑈(𝜃))2+𝑎2)1
2−1
(𝑈02+𝑎2)1
2
𝑆3=∫𝑑𝑈
(𝑈2+𝑎2)3/2
𝑈(𝜃)
𝑈0=1
𝑎2[𝑈(𝜃)
(𝑈(𝜃)2+𝑎2)1
2−𝑈0
(𝑈02+𝑎2)1
2]
A con inuación se a a analiza el o o ipo de exp esión que apa ece en la in eg al.
𝑆2=∫𝑟2𝑑𝑟
(𝑈2+𝑎2)3/2
𝑟(𝜃)
0=∫(𝑈−𝑒𝑟
󰇍
󰇍
󰇍

·(𝑅0
󰇍
󰇍
󰇍
󰇍

−𝑅1
󰇍
󰇍
󰇍
󰇍

))2
(𝑈2+𝑎2)3/2 𝑑𝑈=
𝑈(𝜃)
𝑈0
=∫𝑈2−2𝑈𝑒𝑟
󰇍
󰇍
󰇍

·(𝑅0
󰇍
󰇍
󰇍
󰇍

−𝑅1
󰇍
󰇍
󰇍
󰇍

)+(𝑒𝑟
󰇍
󰇍
󰇍

·(𝑅0
󰇍
󰇍
󰇍
󰇍

−𝑅1
󰇍
󰇍
󰇍
󰇍

))2
(𝑈2+𝑎2)3/2 𝑑𝑈=
𝑈(𝜃)
𝑈0
=∫𝑈2−2𝑈𝑒𝑟
󰇍
󰇍
󰇍

·(𝑅0
󰇍
󰇍
󰇍
󰇍

−𝑅1
󰇍
󰇍
󰇍
󰇍

)+(𝑒𝑟
󰇍
󰇍
󰇍

·(𝑅0
󰇍
󰇍
󰇍
󰇍

−𝑅1
󰇍
󰇍
󰇍
󰇍

))2±𝑎2
(𝑈2+𝑎2)3/2 𝑑𝑈=
𝑈(𝜃)
𝑈0
Donde en el úl imo paso se ha sumado y es ado el e mino de 𝑎2. A con inuación, desa ollando
la in eg al y u ilizando las in eg ales que se han calculado an e io men e, se ob iene.
50
𝑆2=∫𝑑𝑈
(𝑈2+𝑎2)1/2+2𝑈𝑒𝑟
󰇍
󰇍
󰇍

·(𝑅0
󰇍
󰇍
󰇍
󰇍

−𝑅1
󰇍
󰇍
󰇍
󰇍

)𝑆0+
𝑈(𝜃)
𝑈0[2(𝑒𝑟
󰇍
󰇍
󰇍

·(𝑅0
󰇍
󰇍
󰇍
󰇍

−𝑅1
󰇍
󰇍
󰇍
󰇍

))2−|𝑅0
󰇍
󰇍
󰇍
󰇍

−𝑅1
󰇍
󰇍
󰇍
󰇍

|2]𝑆3
Po comodidad y al igual que se ha hecho an e io men e, se a a u iliza la siguien e
nomencla u a.
𝑆4=∫𝑑𝑈
(𝑈2+𝑎2)1/2=𝐿𝑛(√𝑎2+𝑈(𝜃)2+𝑈(𝜃)
√𝑎2+𝑈02+𝑈0)
𝑈(𝜃)
𝑈0
𝑆2=𝑆4+2𝑈𝑒𝑟
󰇍
󰇍
󰇍

·(𝑅0
󰇍
󰇍
󰇍
󰇍

−𝑅1
󰇍
󰇍
󰇍
󰇍

)𝑆0+[2(𝑒𝑟
󰇍
󰇍
󰇍

·(𝑅0
󰇍
󰇍
󰇍
󰇍

−𝑅1
󰇍
󰇍
󰇍
󰇍

))2−|𝑅0
󰇍
󰇍
󰇍
󰇍

−𝑅1
󰇍
󰇍
󰇍
󰇍

|2]𝑆3
Sus i uyendo las de iniciones de 𝑆1 y 𝑆2 en la ecuación in eg al pa a el iángulo bajo es udio se
end ía que:
𝑤𝑖= 𝜙1∫[𝑆1+(𝐶1𝑐𝑜𝑠(𝜃)+𝐷1sin (𝜃))𝑆2]𝑑𝜃+
𝜃3
θ2
+ 𝜙2∫[(𝐶2𝑐𝑜𝑠(𝜃)+𝐷2sin (𝜃))𝑆2]𝑑𝜃+
𝜃3
θ2
+ 𝜙3∫[(𝐶3𝑐𝑜𝑠(𝜃)+𝐷3sin (𝜃))𝑆2]𝑑𝜃
𝜃3
θ2
A con inuación es impo an e eco da que 𝑈=𝑈(𝑟(𝜃)), donde 𝑈(𝑟) es conocida, pe o aún
hay que de e mina 𝑟(𝜃). Pa a ello hay que u iliza el siguien e iángulo.
51
Donde se puede e que 𝑟(𝜃) es la dis ancia desde un pun o del lado 3-2 al é ice 1 del
iángulo. Aplicando elaciones en e los dis in os ángulos, se pueden ob ene los ángulos
in e io es de la egión indicada.
Aplicando el Teo ema del seno y despejando se ob iene la exp esión 𝑟(𝜃) buscada.
𝑟(𝜃)
sin (𝜋+𝜃2−𝜃4)=|𝑅2
󰇍
󰇍
󰇍
󰇍

−𝑅1
󰇍
󰇍
󰇍
󰇍

|
sin(𝜃4−𝜃)
𝑟(𝜃)=sin (𝜃4−𝜃2)|𝑅2
󰇍
󰇍
󰇍
󰇍

−𝑅1
󰇍
󰇍
󰇍
󰇍

|
sin(𝜃4−𝜃)
FIGURA 8: ÁNGULOS EMPLEADOS EN LA DISCRETIZACIÓN. GRÁFICO
EXPLICATIVO II
FIGURA 9: DEFINICIÓN DE 𝒓(𝜽)
52
Una ez se haya hecho la in eg al en la a iable 𝑟, se p ocede a la in eg ación en la a iable 𝜃, la
cual se hace de o ma numé ica median e un mé odo de colocación.
Una ez hecha la in eg ación pa a la p ime a ecuación, ya se ienen los alo es de los
coe icien es de in luencia que acompañan a los po enciales de elocidades. Es os po enciales
an colocados en los é ices de los iángulos en los que se ha di idido el ala y la es ela.
En el caso en que se enga que in eg a en el en o no de la singula idad, se ealiza á el siguien e
mallado:
Y se supond á que el po encial es de la siguien e o ma:
𝜙′(𝑥)=𝜙1′=𝜙0′
Es deci , en el en o no de la singula idad (8 iángulos) el po encial es cons an e e igual al
po encial de la singula idad. Es o e i a la eno me ines abilidad numé ica que end ían las
ecuaciones si se dejase e oluciona el po encial como en el es o de los casos. La in eg al no es
di e gen e en el en o no de la singula idad sólo si el po encial es cons an e.
La condición de con o no de impene abilidad disc e izada queda de la siguien e o ma:
𝑤(𝑥𝑖
󰇍
󰇍
󰇍

)=[𝐴1,1,𝐴1,2,…𝐴1,𝑖,…𝐴1,(𝑁+2)(𝑀+2)]·
[
𝜙1,1
𝜙1,2
…
𝜙1,𝑖
…
𝜙1,(𝑁+2)(𝑀+2)
]
+
+[𝐸1,1,𝐸1,2,…𝐸1,𝑖,…𝐸1,(𝑁+2)(𝑀 𝑒𝑠𝑡𝑒𝑙𝑎)]·
[
𝜙𝐸1,1
𝜙𝐸1,2
…
𝜙𝐸1,𝑖
…
𝜙𝐸1,(𝑁+2)(𝑀 𝑒𝑠𝑡𝑒𝑙𝑎)
]
FIGURA 10: MALLADO EN EL ENTORNO DE LA
SINGULARIDAD
53
Donde 𝜙𝑖,𝑗 indica el po encial en el nodo i,j y 𝐴𝑖,𝑗 y 𝐸𝑖,𝑗 son los co espondien es coe icien es de
in lencia, calculados empleando las exp esiones analí icas de las in eg ales halladas
an e io men e.

54
55
6. PROBLEMA ESTACIONARIO. RESOLUCIÓN Y
CONDICIONES DE CONTORNO
Una ez ealizada la o mulación de la condición de impene abilidad en su o ma gene al, en
es e apa ado se a a pasa a expone las condiciones de con o no pa a el p oblema es aciona io
y su esolución.
Pa a ce a el p oblema, se deben de expone un o al de ecuaciones igual al núme o de
incógni as, con el obje i o de que la ma iz de coe icien es sea cuad ada e in e ible.
Cómo se ha is o an e io men e que en cualquie pun o del dominio luido donde no haya
sólido, se debe cumpli que el po encial de elocidades pe u badas debe de se 0 al se lo en el
in ini o aguas a iba, en pa icula , sob e el bo de de a aque y sob e los bo des ma ginales se
debe de cumpli que, en i ud de la ecuación de Eule -Be nouilli linealizada, que exp esa que la
de i ada sus ancial del po encial pe u bado debe de se 0, se concluye que dicho po encial es
0 en los bo des de a aque y ma ginales.
 Bo de de a aque y bo des ma ginales:
𝜙′(𝑥𝑏𝑎,𝑦)=𝜙′(𝑥,𝑏2)=𝜙′(𝑥,−𝑏2)=0
Po o a pa e, en el bo de de salida, la condición de Ku a-Joukowski exp esa en el caso
es aciona io que:
 Bo de de salida:
𝜕𝜙′
𝜕𝑥=0→𝜙(𝑥𝑏𝑠−1,𝑦)=𝜙(𝑥𝑏𝑠,𝑦)=𝜙(𝑥𝑏𝑠+1=𝑥𝐸1,𝑦)
Es a condición implica que pa a una línea de coo denada y, el alo del po encial de elocidades
en el bo de de salida en esa línea y es el mismo en el bo de de salida, en la disc e ización en x
an e io al bo de de salida y a lo la go de oda la es ela. G á icamen e se puede e en la igu a
siguien e:
56
Se obse a que en como el po encial es el mismo an o en el bo de de salida, como en el pun o
an e io y como en oda la es ela a lo la go de cada línea y, no iene sen ido ealiza un mallado
a lo la go de la es ela pa a el caso es aciona io. Luego en es e caso, el mallado es mucho más
simple:
FIGURA 11: POTENCIAL EN EL BORDE DE SALIDA Y EN LA ESTELA EN RÉGIMEN ESTACIONARIO PARA UN VALOR
DE Y DADO
FIGURA 12: MALLADO TRIANGULAR EN RÉGIMEN ESTACIONARIO
57
A la ho a de mon a las ecuaciones es e ecin o se puede subdi idi de la siguien e o ma
Únicamen e se debe de impone la condición de con o no de impene abilidad en la zona
in e io al ecuad o ama illo, es deci , un o al de N*M ecuaciones. El po encial en el bo de de
salida y en la es ela ( ecuad o azul) debe de se el mismo que en los pun os an e io es al bo de
de salida, y el po encial en los bo des de a aque y ma ginales ( ecuad os e des) es 0.
Así pues se iene un o al de N*(M+1) incógni as sob e el ala, y N*(M+1) ecuaciones di ididas
en N*M condiciones de impene abilidad, más N condiciones de con o no de Ku a sob e el
bo de de salida. Es o pe mi e o ma una ma iz cuad ada e in e ible con la que esol e el
p oblema.
La ecuación queda ía de la siguien e o ma gene al pa a el p oblema es aciona io (se omi en los
é minos que acompañan a po enciales nulos po la condición de con o no).
𝑤′(𝑥𝑖
󰇍
󰇍
󰇍

)=𝑤𝑖′=[𝐴𝑖,2,𝐴𝑖,3,…𝐴𝑖,𝑀+2,…𝐴𝑖,(𝑁+1)(𝑀+2)]·
[
𝜙2,2
𝜙2,3
…
𝜙2,𝑀+2
…
𝜙2,(𝑁+1)(𝑀+2)
]
+
+[𝐸𝑖,2,𝐸𝑖,3,…𝐸𝑖,𝑖,…𝐸𝑖,(𝑁+1)·]·
[
𝜙𝐸2
𝜙𝐸2
…
𝜙𝐸𝑖
…
𝜙𝐸(𝑁+1)
]
FIGURA 13: CONDICIONES DE CONTORNO PARA EL POTENCIAL EN RÉGIMEN ESTACIONARIO
64
ins an e de iempo simplemen e median e la “con ección” aguas debajo de los po enciales
ob enidos en el bo de de salida.
En el anexo se p opo ciona el p og ama de Ma lab que esuel e el p oblema es aciona io y no
es aciona io pa a geome ías gené icas.

65
8. SOLUCIÓN EN RÉGIMEN ESTACIONARIO.
8.1. Alas Rec angula es
A con inuación se a a expone la solución pa a alas ec angula es en égimen es aciona io, y
pa a dis in os ala gamien os. Pa a ob ene una solución asequible es necesa io alcanza un
equilib io en e el núme o de di isiones del mallado y el iempo compu acional, de mane a que
se in en e alcanza la mejo solución (o una solución con p ecisión su icien e) en un iempo de
cálculo azonable.
Pa a alcanza es e obje i o es p eciso eco da que la dis ibución de p esiones en el ex adós
del ala su e unos g adien es muy acen uados en las p oximidades del bo de de a aque y de los
bo des ma ginales. Es e enómeno lle a a pensa que es necesa io e ina el mallado en es as
zonas, con obje i o de ob ene unos alo es azonables del coe icien e de sus en ación y de la
pendien e de la cu a de sus en ación (Se supond á en odo momen o a iación lineal del
coe icien e de sus en ación con el ángulo de a aque). El mallado en el eje x e y se ob iene con
las siguien es exp esiones ( eniendo en cuen a el núme o de di isiones):
𝑚𝑦=0.5𝑏·cos((𝑗−1)𝜋
𝑁+1) 𝑐𝑜𝑛 𝑗∈[1,𝑁+2]
𝑚𝑥=(1−cos((𝑖−1)𝜋
2(𝑀−1))) 𝑐𝑜𝑛 𝑖∈[1,𝑀]
El mallado pa a el caso de alas en lecha es simila .
Se mues a un ejemplo de cómo se ían las di isiones del mallado. Se ha p esen ado el ejemplo
pa a que se obse e el ala gamien o del ala y se pueda in ui la longi ud de la es ela (la es ela
FIGURA 17: DISCRETIZACIÓN EN ALAS RECTANGULARES
66
llega ía más lejos, pe o la idea con es a imagen es mos a al lec o las dimensiones del
p oblema).
En la igu a an e io se obse a cómo el mallado es más ino en los pun os del bo de de a aque
y ma ginales pa a cap a bien los g adien es.
Resol iendo el p oblema pa a un ángulo de a aque de 5º, un ala gamien o de 5 (b=5, c=1) y
es echamien o nulo, se ob ienen las siguien es dis ibuciones del po encial sob e la supe icie
del ala.
Si se ealiza la de i ación del po encial espec o a la coo denada x, se puede ob ene el campo
de elocidades sob e el ex adós.
Es as dos igu as mues an el eno me g adien e de elocidades que exis e an o en los bo des
ma ginales, como en el bo de de a aque y po an o jus i ican el mallado ealizado.
La pendien e de la cu a de sus en ación ob enida pa a es a con igu ación es de:
FIGURA 18: POTENCIAL SOBRE EL EXTRADÓS DE ALAS RECTANGULARES
FIGURA 19: VELOCIDAD PERTURBADA EN EL EJE X SOBRE EL EXTRADÓS DE ALAS RECTANGULARES.
67
𝐶𝐿𝛼
Nx
Ny
AR
3,8741
15
21
5
TABLA 1: RESULTADOS PARA ALA RECTANGULAR EN RÉGIMEN ESTACIONARIO
A con inuación se ha ealizado el mismo cálculo pe o aumen ando el mallado has a un alo de
Nx=25 y Ny=41 pa a obse a las di e encias.
Se mues a a con inuación una igu a del campo de elocidades pa a es e caso, de o ma que se
ap ecie mejo el g adien e que se p oduce en el bo de de a aque y en los bo des ma ginales.
Los esul ados pa a es e caso son:
𝐶𝐿𝛼
Nx
Ny
AR
3,9096
25
41
5
TABLA 2: RESULTADOS PARA ALA RECTANGULAR EN RÉGIMEN ESTACIONARIO. MALLADO REFINADO
El iempo compu acional ha sido mayo pa a es e caso, y la mejo a se no a en la e ce a ci a
signi ica i a.
Se puede obse a que la pendien e de la cu a de sus en ación iende a aumen a con o me el
mallado se hace más ino. Es o se debe a que se cap a mejo las e oluciones de las a iables
luidas en el en o no de los bo des ma ginales y del bo de de a aque, pe mi iendo que la
in eg ación de luga a un esul ado más p eciso.
De aquí en adelan e se ealiza á el mallado con los siguien es alo es de Nx=15 y Ny=8*b+1. La
azón de es os alo es es po que cuando se esuel a el p oblema no es aciona io, el iempo
FIGURA 20: VELOCIDAD PERTURBADA EN EL EJE X SOBRE EL EXTRADÓS DE ALAS RECTANGULARES.
REFINAMIENTO
68
compu acional aumen a eno memen e pa a mallados muy g andes si se quie en ob ene
alo es que mejo an poco en p ecisión con espec o a es os. De o ma que con es e mé odo
aco amos el e o a la e ce a ci a signi ica i a mien as man enemos un iempo azonable.
La azón de la ecuación Ny=8*b+1, es iba en que en el análisis de alas con g an ala gamien o,
se debe án hace más di isiones en la coo denada y debido a que al aumen a la en e gadu a,
ambién lo hacen los paneles en los que se di ide el ala en esa di ección y po an o se pe de ía
p ecisión si no se aumen a el mallado.
8.2. Alas con Flecha
En es e apa ado se a a la solución del p oblema es aciona io pa a el caso de alas con lecha.
De es a o ma, dando alo es de S, AR, 𝛹 y E, queda de inida el ala.
Pa a el caso de alas en lecha, se p obó a u iliza un mallado simila al caso de alas ec angula es,
sin emba go la p ecisión adop ada en compa ación con los alo es ípicos de la li e a u a no e a
su icien e.
La azón es iba en que las alas en lecha poseen un g adien e de elocidades en las
p oximidades de la cue da aíz del ala, que esul aba se la zona donde el mallado an e io es
menos p eciso (elemen os más g andes, cap aban mal la e olución).
Pa a la comp obación de que es e código e a e ec i o, puede comp oba el lec o in e esado
que si se u iliza el código pa a el ala en lecha con las pa icula idades de E=1 y ángulo de lecha
nulo, se uel e a ob ene los alo es que se ob u ie on pa a el ala con o ma en plan a
ec angula .
Pa a sol en a es e p oblema se ecu ió a un e inado ambién po dicha zona simila a lo que
se había hecho pa a los bo des de a aque y ma ginales. Se mues a el mallado pa icula
ealizado pa a las alas en lecha:
FIGURA 21: DISCRETIZACIÓN EN ALAS CON FLECHA
69
El po encial de elocidades sob e es e ipo de alas iene la siguien e o ma:
De i ando el po encial, se puede ob ene el campo de elocidades pe u badas sob e la
supe icie de ex adós, el cual se mues a en la siguien e igu a:
En es a igu a se puede obse a el g adien e de elocidades que se p oduce en las p oximidades
de la cue da aíz del ala.
𝐶𝐿𝛼
Nx
Ny
AR
3,5368
15
41
5
TABLA 3: RESULTADOS PARA ALA CON FLECHA 30º EN RÉGIMEN ESTACIONARIO.
Se obse a que el alo de la pendien e de la cu a de sus en ación pa a un ala en lecha es
meno que pa a el ala ec angula , como e a de espe a .
FIGURA 22: POTENCIAL SOBRE EL EXTRADÓS DE ALAS CON FLECHA
FIGURA 23: VELOCIDAD PERTURBADA EN EL EJE X SOBRE EL EXTRADÓS DE ALAS CON FLECHA.

70
71
9. SOLUCIÓN EN RÉGIMEN NO ESTACIONARIO.
9.1. Alas Rec angula es
A con inuación se a a mos a los esul ados ob enidos pa a alas ec angula es en égimen no
es aciona io. Pa a ello, se a a ealiza el mismo análisis, con el mismo mallado en el ala, pe o
pa a égimen no es aciona io.
El inc emen o de iempo Δ𝑡 debe de se in e io a 𝑐/𝑈∞, pa a que el e ec o no es aciona io se
obse e co ec amen e.
El iempo de cálculo se ha es ingido a 200 𝑐/𝑈∞, con un Δ𝑡 de 0.25𝑐/𝑈∞, en e cada
i e ación. Las ecuaciones han sido in eg adas en el iempo u ilizando un mé odo de Eule de
p ime o den.
El alo ob enido se mues a a con inuación jun o con la cu a que indica la e olución que ha
seguido el alo de la pendien e de la cu a de sus en ación desde el ins an e inicial has a el
inal pa a el caso de un ala ec angula con elación de aspec o AR=5 y suponiendo que la
elocidad adimensional en el in ini o es una unción escalón en el iempo.
𝐶𝐿𝛼
Nx
Ny
AR
3,8788
15
41
5
TABLA 4: RESULTADOS PARA ALA RECTANGULAR EN RÉGIMEN NO ESTACIONARIO
Cabe esal a el hecho de que la pendien e de la cu a de sus en ación alcanza el alo ob enido
en égimen es aciona io an e io men e.
FIGURA 24: EVOLUCIÓN DE LA PENDIENTE DE LA CURVA DE SUSTENTACIÓN CON EL TIEMPO. ALA
RECTANGULAR
72
9.2. Alas con Flecha
Pa a las alas en lecha se ha seguido un p ocedimien o simila al de las alas ec angula es, con
la pa icula idad de hace un mallado más ino en la zona ce cana a la cue da aíz del ala.
Una ez se iene implemen ado el código pa a la esolución no es aciona ia de alas
ec angula es, simplemen e se modi ica la geome ía pa a ealiza el es udio del
compo amien o no es aciona io pa a alas en lecha.
Se ha ealizado de nue o una es icción en el iempo de cálculo de 200 𝑐/𝑈∞, con un Δ𝑡 de
0.25𝑐/𝑈∞. Es o da paneles de longi ud 0.25 en la es ela y da alo es lo su icien emen e p ecisos
aunque con un iempo compu acional ele ado.
𝐶𝐿𝛼
Nx
Ny
AR
3,5303
15
41
5
TABLA 5: RESULTADOS PARA ALA CON FLECHA 45º EN RÉGIMEN ESTACIONARIO
Pa a que el lec o enga una idea de las dimensiones que iene la es ela no es aciona ia analizada
pa a el case de un ala en lecha se adjun a la igu a siguien e:
FIGURA 25: EVOLUCIÓN DE LA PENDIENTE DE LA CURVA DE SUSTENTACIÓN CON EL TIEMPO. ALA CON FLECHA
73
Cada di isión en el eje x co esponde a un ins an e de iempo. En o al se ienen 200 ins an es
de iempo en di isiones de 0.25 s.
FIGURA 26: ESTELA NO ESTACIONARIA EN CASO DE ALA EN FLECHA
80
El e o ob enido en e el cálculo es aciona io y el ob enido del g á ico se mues an a
con inuación:
ALA RECTANGULAR
AR
𝐶𝐿𝛼𝐸𝑠𝑡𝑎𝑐𝑖𝑜𝑛𝑎𝑟𝑖𝑜
KyP
E o (%)
4
3,5656
3,6
0,95
5
3,8741
3,9
0,66
6
4,1321
4,2
1,61
7
4,3236
4,3
0,54
TABLA 9: ERROR COMETIDO EN ALAS RECTANGULARES EN RÉGIMEN ESTACIONARIO
ALA FLECHA 30º
AR
𝐶𝐿𝛼𝐸𝑠𝑡𝑎𝑐𝑖𝑜𝑛𝑎𝑟𝑖𝑜
KyP
E o (%)
4
3,2911
3,35
1,75
5
3,5368
3,7
4,41
6
3,7222
3,85
3,31
7
3,8712
4
3,22
TABLA 10: ERROR COMETIDO EN ALAS CON 𝚿=𝟑𝟎º EN RÉGIMEN ESTACIONARIO
ALA FLECHA 45º
AR
𝐶𝐿𝛼𝐸𝑠𝑡𝑎𝑐𝑖𝑜𝑛𝑎𝑟𝑖𝑜
KyP
E o (%)
4
2,922
3
2,6
5
3,0942
3,2
3,30625
6
3,2197
3,35
3,8895
7
3,3231
3,5
5,054
TABLA 11: ERROR COMETIDO EN ALAS CON 𝚿=𝟒𝟓º EN RÉGIMEN ESTACIONARIO
Po lo que se puede obse a que el mé odo desa ollado iene un e o máximo del 5% pa a el
caso de alas con un al o alo del ángulo de lecha. Es o es debido a que el mallado pa a es e
ipo de alas debe ía de es a mejo adap ado y e inado. Pe o hace es os e inamien os
aumen a el iempo compu acional.

81
11. CONCLUSIONES
A la is a de los esul ados se puede conclui que el mé odo p esen ado en es e abajo supe a
ampliamen e las expec a i as iniciales. Se ha desa ollado un mé odo que consigue una
p ecisión de cálculo al a, que a oja alo es simila es a los ob enidos expe imen almen e y a los
que se pueden ob ene con o os p og amas.
El mé odo desa ollado p esen a una g an obus ez, pe mi iendo el análisis de alas de cualquie
geome ía y en ambos egímenes: es aciona io y no es aciona io.
También se han ealizado p uebas cuando el ala se some e a dis in os mo imien os en pleno
uelo. Po ejemplo las siguien es si uaciones ísicas.
- Ala a anca y pasa de es a en eposo a un ángulo de a aque 𝛼0, a mo e se con una
elocidad 𝑈∞. En un ins an e 0, el ángulo de a aque pasa a se ins an áneamen e el
doble. (Compo amien o an e un escalón)
- Ala uela a una elocidad 𝑈∞y iene además un mo imien o pe iódico en la di ección
e ical Z( ). Angulo de a aque pe manece nulo.
Pa a el caso de espues a an e un escalón en un ins an e de e minado, en el cual el ángulo de
a aque pasa de ale α=5º a un alo de α=10º, se ha ob enido la siguien e e olución de la
pendien e de la cu a de sus en ación.
FIGURA 32: RESPUESTA ANTE UN ESCALÓN
Se obse a que as la pe u bación se p oduce un descenso acusado del alo de la pendien e
de la cu a de sus en ación. Es o es lógico debido a que el ángulo de a aque pasa a se el doble
de o ma ins an ánea. La azón de que la pendien e de la cu a de sus en ación no sea
exac amen e la mi ad del alo pa a el ala en lecha, es po que el coe icien e de sus en ación
aumen a, pe o lo hace e iden emen e de o ma menos ab up a que el ángulo de a aque.
AR=4
E=1
Ψ=0
82
T as la pe u bación se puede e como el alo de la pendien e de la cu a de sus en ación
uel e a ende al alo que end ía en égimen es aciona io como e a de espe a (se puede
comp oba que el alo lími e coincide con el de la abla, ya que se ecue da que es e alo es
cons an e en égimen es aciona io según la eo ía po encial linealizada.
En cuan o a la e olución pa a el caso de mo imien o oscila o io a ángulo de a aque dado (5º)
se ha ob enido la siguien e g á ica
FIGURA 33: EVOLUCIÓN DEL COEFICIENTE DE SUSTENTACIÓN EN UN MOVIMIENTO OSCILATORIO
Donde: Ω= 𝑤𝑐
2𝑈∞
Es la ecuencia adimensional y 𝑤 la ecuencia dimensional.
Es as p uebas ya no se han podido es udia a ondo po al a de iempo y po que quedaban
ue a del obje i o p incipal del abajo. Se deja pa a es udios pos e io es la alidación de es os
esul ados.
Las condiciones mos adas an e io men e son las que se mues an en la igu a siguien e, la
cual se puede encon a en Ka z y Plo kin, 1991 [2]
FIGURA 34: EVOLUCIÓN DEL COEFICIENTE DE SUSTENTACIÓN EN UN MOVIMIENTO OSCILATORIO. (KYP) [2]
AR=4
E=1
Ψ=0
H0=0.1
Ω=0.5
α =-5º
83
12. DESARROLLOS FUTUROS
T as es e p oyec o se p oponen los siguien es desa ollos:
- Validación del mé odo con esul ados ob enidos de la li e a u a pa a el caso de
mo imien os e icales Z( ) y de cabeceo 𝛼0( ) y es udio del caso de de lexión de
supe icies hipe sus en ado as simples como laps o sla s
- Implemen ación de las ecuaciones de la elas icidad pa a esol e el p oblema
ae oelás ico pa a alas con dis in a geome ía y calcula elocidades de Flu e y de
di e gencia. Es udio del e ec o de la lecha en dichas elocidades.
- Op imización de la posición de los ale ones pa a ob ene un momen o de balance
deseado den o de los lími es de la no ma i a pa a alas de dis in a geome ía.
- Implemen ación de las ecuaciones de la elas icidad pa a esol e el p oblema
ae oelás ico pa a alas con dis in a geome ía y calcula elocidades de in e sión de
mando.
84
85
REFERENCIAS
[1]
José Manuel Go dillo A ias de Saa ed a y Guillaume Riboux Ache , "In oducción a la
Ae odinámica Po encial" Pa anin o,2012.
[2]
Joseph Ka z y Allen Plo kin, “Low-Speed Ae odynamics: F om Wing Theo y o Panel
Me hods”, McG aw-Hill, 1991.
[3]
Dowell, E. H., E. F. C awley, H. C. Cu iss, J ., D. A. Pe e s, R. H. Scanlan, and F. Sis o,
“A Mode n Cou se in Ae oelas ici y”, 3 d ed., Kluwe Academic Publishe s, 1995.
[4]
Flo es Caballe o, Míguel Ángel, “Aplicación de la Teo ía Po encial Linealizada pa a el
es udio de las ue zas eje cidas po una co ien e inciden e sob e placas ec angula es
oscilan es”, Se illa 2016

86
87
ANEXOS
1. CÓDIGO PRINCIPAL DE CÁLCULO PARA EL CASO NO ESTACIONARIO DE ALA EN
FLECHA.
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% CÁLCULO DE COEFICIENTE DE SUSTENTACIÓN PARA ALAS EN FLECHA %%%%%%%%%%%
%% %%%%%%%%% RÉGIMEN NO ESTACIONARIO %%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% AUTOR: Fe nando Mo eno Pino %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% TUTOR : José Manuel Go dillo %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clea all; close all; clc
%% %%% INTRODUCCIÓN DE LOS PARÁMETROS
alpha0 = 5*pi/180;
ecx=[-0.238619186083197; -0.661209386466265; -0.932469514203152;0.238619186083197;
0.661209386466265; 0.932469514203152];
ecw=[0.467913934572691; 0.360761573048139; 0.171324492379170; 0.467913934572691;
0.360761573048139; 0.171324492379170];
dospi=1/2/pi;
o i=1:6
ec he ak(i)=0.5*pi*( ecx(i)+1);
end
%% GEOMETRÍA DEFINICION DE NUMERO DE PANELES Y DE INCOGNITAS
AR=5; % Ala gamien o
S=AR; % Supe icie
E=1; % Es echamien o
b=sq (S*AR); % en e gadu a
psi = 0*pi/180; % Angulo de lecha
c = 2*sq (S/AR)/(1+E); % Cue da
iala = 15; % Nume o de líneas a lo la go del eje x en el ala
Nyincog = 4*b+1; % Nume o de incógni as a lo la go del eje y
Nyincog2 = (Nyincog+1)/2;
c =c*E; % Cue da en los ex emos
i in = iala+1; % P ime a di ision de la es ela
d he ax = pi/(2*(iala-1)); % Di isión de angulo que ba e las coo denadas X
d he ay = pi/(Nyincog+1); % Di isión de ángulo que ba e las coo denadas Y
d he ay2= pi/(Nyincog2);
del ay = b/(Nyincog+1); % Sepa ación en el mallado
del ax = c/(iala-1);
Nxec = iala-2;
Nyec = Nyincog;
88
NTec = Nxec*Nyec;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% GEOMETRÍA. DEFINICION DE LA FORMA EN PLANTA Y MALLADO %%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% MALLADO DOBLE EN Y
y(1:Nyincog2+1)=(0.25*b+0.25*b*cos(((1:Nyincog2+1)-1)*d he ay2));
y(Nyincog2+1:Nyincog+2)=(-0.25*b+0.25*b*cos(((1:Nyincog2+1)-1)*d he ay2));
xa = abs(y)'* an(psi); % pun os del bo de de a aque
xs = c+(b/2* an(psi)+c -c)/(b/2)*abs(y)'; % bo de de salida
i in = iala+1; % Final de la es ela
Nxec = iala-2;
Nyec = Nyincog;
NTec = Nxec*Nyec; % El nume o o al de ecuaciones donde impone la elocidad
e ical
mx = ze os (Nyincog+2, iala+1); % Ma iz que almacena á odas las componen es x de los pun os del
mallado
my = ze os (Nyincog+2, iala+1); % Ma iz que almacena á odas las componen es y de los pun os del
mallado
UNSTEADY=1;
i UNSTEADY ==0
ies ela=1;
inT=0;
D =20*b/ies ela;
else
D =0.25;
ies ela=200; % longi ud de la es ela
inT=ies ela;
end
conT=0;
o j=1:iala
my(:,j)=y';
end
o j= 1:Nyincog+2
o i=1:iala
mx(j,i) = xa(j)+(xs(j)-xa(j))*(1-cos((i-1)*d he ax));
end
end
o j=1:Nyincog+2
o i=1:ies ela-1
i UNSTEADY ==0;
ecy(j) = my(Nyincog+3-j, iala);
89
mx (j,iala+1)=mx(j,iala)+100*b;
my(j,iala+1) = (0.5*b*cos((j-1)*d he ay));
else
mx (j,iala+i)=mx(j,iala)+i*D ;
my(j,iala+i) = my(j,iala);
end
end
end
mxS=mx(:,iala:iala+1); % P ime a anja de es ela en x
myS=my(:,iala:iala+1); % P ime a anja de es ela en y
mxE=mx(:,iala+1:end); %ma ices de la es ela en x
myE=my(:,iala+1:end); % ma ices de la es ela en y
Dxbs=mx(1,iala)-mx(1,iala-1);
%% PINTAR ALA %%
igu e (10)
o i = 1:iala % Es e bucle se enca ga de pin a las lineas en e ical de los paneles
plo ( mx(:,i), my(:, i),' ')
hold on
end
o j = 1:Nyincog+2 % Es e bucle pin a las lineas ho izon ales de los paneles
plo (mx(j,1:iala), my(j,1:iala))
hold on
end
%% PINTAR ESTELA %%
igu e(11)
o i = iala: ies ela+iala-1 % Es e bucle se enca ga de pin a las lineas en e ical de los paneles
plo ( mx(:,i), my(:, i),' ')
hold on
end
o j = 1:Nyincog+2 % Es e bucle pin a las lineas ho izon ales de los paneles
plo (mx(j,iala:ies ela+iala-1 ), my(j,iala:ies ela+iala-1))
hold on
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% CREACION DE MATRICES %%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%
sid = ((iala-2)*(Nyincog)+Nyincog); %Dimensiones de la ma iz de coe icien es en ala
%(en "sid" pun os se impone la condicion de con o no de
%impene abilidad)
sidE=((iala-2)*(Nyincog)); % Dimensión de la ma iz de coe icien es en es ela
wEIm1=ze os(sid,1); % Velocidad e ical c eada po es ela
% debido a o bellinos inyec ados en la es ela has a el ins an e -1
o j=1:Nyincog
96
ila(1) = j;
columna(1)=i;
ila(2)=j-1;
columna(2)=i;
ila(3)=j-1;
columna(3)=i-1;
ila(4)=j;
columna(4)=i-1;
end
end
o iangulo=1:2
x1=mx( ila(1), columna(1));
y1=my( ila(1), columna(1));
x2=mx( ila( iangulo+1), columna( iangulo+1));
y2=my( ila( iangulo+1), columna( iangulo+1));
x3=mx( ila( iangulo+2), columna( iangulo+2));
y3=my( ila( iangulo+2), columna( iangulo+2));
i i<=iala %No es oy en la es ela
columnaphi1= columna(1)-1+( ila(1)-2)*(iala-1);
columnaphi2= columna( iangulo+1)-1+( ila( iangulo+1)-2)*(iala-1);
columnaphi3= columna( iangulo+2)-1+( ila( iangulo+2)-2)*(iala-1);
end
D=(x2-x1)*(y3-y1)-(y2-y1)*(x3-x1);
c1=(y2-y3)/D;
d1=(x3-x2)/D;
c2=(y3-y1)/D;
d2=(x1-x3)/D;
c3 =(y1-y2)/D;
d3 =(x2-x1)/D;
s
acphi1=dospi*( I(1)+c1* I(3)+d1* I(2));
acphi2=dospi*(c2* I(3)+d2* I(2));
acphi3=dospi*(c3* I(3)+d3* I(2));
i ila(1)>1 && ila(1)<Nyincog+2 && columna(1)>1
Ma ( 0,columnaphi1)=Ma ( 0,columnaphi1)+ acphi1;
End
i ila( iangulo+1)>1 && ila( iangulo+1)<Nyincog+2 && columna( iangulo+1)>1
Ma ( 0,columnaphi2)=Ma ( 0,columnaphi2)+ acphi2;
End
i ila( iangulo+2)>1 && ila( iangulo+2)<Nyincog+2 && columna( iangulo+2)>1
Ma ( 0,columnaphi3)=Ma ( 0,columnaphi3)+ acphi3;
end
end

97
3. Mon aje_SOLAPE
%% FUNCION el mon aje en el bo de de salida y comienzo de la es ela
x0=mx(l, k); y0=my(l,k); % Coo denadas de la singula idad
0 =k-1+(l-2)*(iala-2);
Nyincog2=(Nyincog-1)/2;
ila=ze os(1,4);
columna=ze os(1,4);
i j<Nyincog2+3
i mod(i,2)==0
ipo ec angulo=1;
ila(1) = j;
columna(1)=i;
ila(2)=j-1;
columna(2)=i;
ila(3)=j-1;
columna(3)=i-1;
ila(4)=j;
columna(4)=i-1;
else
ipo ec angulo=2;
ila(4) = j;
columna(4)=i;
ila(1)=j-1;
columna(1)=i;
ila(2)=j-1;
columna(2)=i-1;
ila(3)=j;
columna(3)=i-1;
end
else
i mod(i,2)==0
ipo ec angulo=2;
ila(4) = j;
columna(4)=i;
ila(1)=j-1;
columna(1)=i;
ila(2)=j-1;
columna(2)=i-1;
ila(3)=j;
columna(3)=i-1;
else
ipo ec angulo=1;
ila(1) = j;
columna(1)=i;
98
ila(2)=j-1;
columna(2)=i;
ila(3)=j-1;
columna(3)=i-1;
ila(4)=j;
columna(4)=i-1;
end
end
o iangulo=1:2
x1=mxS( ila(1), columna(1));
y1=myS( ila(1), columna(1));
x2=mxS( ila( iangulo+1), columna( iangulo+1));
y2=myS( ila( iangulo+1), columna( iangulo+1));
x3=mxS( ila( iangulo+2), columna( iangulo+2));
y3=myS( ila( iangulo+2), columna( iangulo+2));
i i<=2 %No es oy en la es ela
columnaphi1= columna(1)+( ila(1)-2)*(2);
columnaphi2= columna( iangulo+1)+( ila( iangulo+1)-2)*(2);
columnaphi3= columna( iangulo+2)+( ila( iangulo+2)-2)*(2);
end
D=(x2-x1)*(y3-y1)-(y2-y1)*(x3-x1);
c1=(y2-y3)/D;
d1=(x3-x2)/D;
c2=(y3-y1)/D;
d2=(x1-x3)/D;
c3 =(y1-y2)/D;
d3 =(x2-x1)/D;
s
acphi1=dospi*( I(1)+c1* I(3)+d1* I(2));
acphi2=dospi*(c2* I(3)+d2* I(2));
acphi3=dospi*(c3* I(3)+d3* I(2));
i ila(1)>1 && ila(1)<Nyincog+2 && columna(1)>1
Ma SOLAPE( 0,columnaphi1)=Ma SOLAPE( 0,columnaphi1)+ acphi1;
End
i ila( iangulo+1)>1 && ila( iangulo+1)<Nyincog+2 && columna( iangulo+1)>0
Ma SOLAPE( 0,columnaphi2)=Ma SOLAPE( 0,columnaphi2)+ acphi2;
End
i ila( iangulo+2)>1 && ila( iangulo+2)<Nyincog+2 && columna( iangulo+2)>0
Ma SOLAPE( 0,columnaphi3)=Ma SOLAPE( 0,columnaphi3)+ acphi3;
end
end
99
4. Mon aje_Es ela_2
%% FUNCION PARA COLOCAR EN LA MATRIZ DE LA ESTELA
x0=mx(l, k); y0=my(l,k); % Coo denadas de la singula idad
0 =k-1+(l-2)*(iala-2);
Nyincog2=(Nyincog-1)/2;
ila=ze os(1,4);
columna=ze os(1,4);
i j<Nyincog2+3
i mod(i,2)==0
ipo ec angulo=1;
ila(1) = j;
columna(1)=i;
ila(2)=j-1;
columna(2)=i;
ila(3)=j-1;
columna(3)=i-1;
ila(4)=j;
columna(4)=i-1;
else
ipo ec angulo=2;
ila(4) = j;
columna(4)=i;
ila(1)=j-1;
columna(1)=i;
ila(2)=j-1;
columna(2)=i-1;
ila(3)=j;
columna(3)=i-1;
end
else
i mod(i,2)==0
ipo ec angulo=2;
ila(4) = j;
columna(4)=i;
ila(1)=j-1;
columna(1)=i;
ila(2)=j-1;
columna(2)=i-1;
ila(3)=j;
columna(3)=i-1;
else
ipo ec angulo=1;
ila(1) = j;
columna(1)=i;
100
ila(2)=j-1;
columna(2)=i;
ila(3)=j-1;
columna(3)=i-1;
ila(4)=j;
columna(4)=i-1;
end
end
o iangulo=1:2
x1=mxE( ila(1), columna(1));
y1=myE( ila(1), columna(1));
x2=mxE( ila( iangulo+1), columna( iangulo+1));
y2=myE( ila( iangulo+1), columna( iangulo+1));
x3=mxE( ila( iangulo+2), columna( iangulo+2));
y3=myE( ila( iangulo+2), columna( iangulo+2));
i i<=ies ela-1 %No es oy en la es ela
columnaphi1= columna(1)+( ila(1)-2)*(ies ela-1);
columnaphi2= columna( iangulo+1)+( ila( iangulo+1)-2)*(ies ela-1);
columnaphi3= columna( iangulo+2)+( ila( iangulo+2)-2)*(ies ela-1);
end
D=(x2-x1)*(y3-y1)-(y2-y1)*(x3-x1);
c1=(y2-y3)/D;
d1=(x3-x2)/D;
c2=(y3-y1)/D;
d2=(x1-x3)/D;
c3 =(y1-y2)/D;
d3 =(x2-x1)/D;
s
acphi1=dospi*( I(1)+c1* I(3)+d1* I(2));
acphi2=dospi*(c2* I(3)+d2* I(2));
acphi3=dospi*(c3* I(3)+d3* I(2));
i ila(1)>1 && ila(1)<Nyincog+2 && columna(1)>1
Ma E( 0,columnaphi1)=Ma E( 0,columnaphi1)+ acphi1;
End
i ila( iangulo+1)>1 && ila( iangulo+1)<Nyincog+2 && columna( iangulo+1)>0
Ma E( 0,columnaphi2)=Ma E( 0,columnaphi2)+ acphi2;
end
i ila( iangulo+2)>1 && ila( iangulo+2)<Nyincog+2 && columna( iangulo+2)>0
Ma E( 0,columnaphi3)=Ma E( 0,columnaphi3)+ acphi3;
end
end
101
5. Mon aje_0
%% FUNCION QUE COLOCA LOS VALORES DE LA MATRIZ CUANDO ANALIZAMOS LA SINGULARIDAD
x0=mx(l,k);
y0=my(l,k);
x1=x0;
y1=y0;
ila(1)=l;
columna(1)=k;
ila(2)=l;
columna(2)=k+1;
ila(3)=l-1;
columna(3)=k+1;
ila(4)=l-1;
columna(4)=k;
ila(5)=l-1;
columna(5)=k-1;
ila(6)=l;
columna(6)=k-1;
ila(7)=l+1;
columna(7)=k-1;
ila(8)=l+1;
columna(8)=k;
ila(9)=l+1;
columna(9)=k+1;
ila(10)=l;
columna(10)=k+1;
ila0=k-1+(l-2)*(iala-2);
columna0=k-1+(l-2)*(iala-1);
o iangulo=1:8
x2=mx( ila( iangulo+1), columna( iangulo+1));
y2=my( ila( iangulo+1), columna( iangulo+1));
x3=mx( ila( iangulo+2), columna( iangulo+2));
y3=my( ila( iangulo+2), columna( iangulo+2));
columnaphi2= columna( iangulo+1)-1+( ila( iangulo+1)-2)*(iala-1);
columnaphi3= columna( iangulo+2)-1+( ila( iangulo+2)-2)*(iala-1);
D=(x2-x1)*(y3-y1)-(y2-y1)*(x3-x1);
c1=(y2-y3)/D;
d1=(x3-x2)/D;
c2=(y3-y1)/D;
d2=(x1-x3)/D;
c3=(y1-y2)/D;
d3=(x2-x1)/D;
s_0
acphi1=dospi*( I(1));
Ma ( ila0,columna0)=Ma ( ila0,columna0)+ acphi1;
end

102
6. FS
%% PROGRAMA QUE REALIZA LA INTEGRACIÓN EN CASO DE ESTAR EN UN PUNTO QUE NO ES LA
%%SINGULARIDAD %%%
1 0x=x1-x0;
1 0y=y1-y0;
2 1x=x2-x1;
2 1y=y2-y1;
3 1x=x3-x1;
3 1y=y3-y1;
3 2x=x3-x2;
3 2y=y3-y2;
L12=sq ( 2 1x* 2 1x+ 2 1y* 2 1y);
L10=sq ( 1 0x* 1 0x+ 1 0y* 1 0y);
he a2=angulo( 2 1x, 2 1y);
he a3=angulo( 3 1x, 3 1y);
i (abs( he a3)<1e-12 && he a3< he a2)
he a3=2*pi;
end
he a4=angulo( 3 2x, 3 2y);
i he a2<0|| he a3<0|| he a4<0
disp('pa a')
pause
end
I1k=0;
I2k=0;
I3k=0;
o p=1:6
he ak=0.5*( he a3- he a2)*( ecx(p)+1)+ he a2;
k=sin( he a4- he a2)*L12/(sin( he a4)*cos( he ak)-cos( he a4)*sin( he ak));
x0k=cos( he ak)* 1 0x+sin( he ak)* 1 0y;
xk= k+ x0k;
a2k=L10*L10- x0k* x0k;
S0k=1/sq ( xk* xk+ a2k)-1/L10;
i (abs( a2k)>1e-8)
S3k=(1/ a2k)*( xk/sq ( a2k+ xk* xk)- x0k/sq ( a2k+ x0k* x0k));
S4k=log((sq ( a2k+ xk* xk)+ xk)/(sq ( a2k+ x0k* x0k)+ x0k));
S1k=-S0k- x0k*S3k;
S2k=S4k-2* x0k*(-S0k- x0k*S3k)-L10*L10*S3k;
else
S3k=-0.5*(1/(( k+L10)*( k+L10))-1/(L10*L10));
S4k=log((L10+ k)/L10);
S1k=-S0k- x0k*S3k;
S2k=S4k-2* x0k*(-S0k- x0k*S3k)-L10*L10*S3k;
103
end
I1k= I1k+S1k* ecw(p);
I2k= I2k+S2k*sin( he ak)* ecw(p);
I3k= I3k+S2k*cos( he ak)* ecw(p);
end
I(1)=0.5*( he a3- he a2)* I1k;
I(2)=0.5*( he a3- he a2)* I2k;
I(3)=0.5*( he a3- he a2)* I3k;
104
7. FS_0
%% INTEGRACION EN CASO DE SINGULARIDAD
2 1x= x2-x1;
2 1y= y2-y1;
3 1x= x3-x1;
3 1y=y3-y1;
3 2x=x3-x2;
3 2y=y3-y2;
L12 = sq ( 2 1x* 2 1x+ 2 1y* 2 1y);
he a2=angulo( 2 1x, 2 1y);
he a3=angulo( 3 1x, 3 1y);
i (abs( he a3)<1e-12 && he a3< he a2)
he a3=2*pi;
end
he a4=angulo( 3 2x, 3 2y);
I1k=0;
o p=1:6
he ak=0.5*( he a3- he a2)*( ecx(p)+1)+ he a2;
k=sin( he a4- he a2)*L12/(sin( he a4)*cos( he ak)-cos( he a4)*sin( he ak));
I1k= I1k+(-1/ k)* ecw(p);
end
I(1)=0.5*( he a3- he a2)* I1k;
105
8. Ángulo
%% FUNCION QUE COGE EL ÁNGULO CORRECTO PARA LA INTEGRACION %%
unc ion [ he a]=angulo( x, y)
he a0=a an(abs( y)/abs( x));
i (abs( x)>1e-15)
he a0=a an(abs( y)/abs( x));
else
i ( y>0)
he a=pi/2.;
else
he a=3*pi/2.;
end
end
i (abs( y)<1e-15)
i ( x>0)
he a=0.;
else
he a=pi;
end
else
i ( y>1e-15 && x>1e-15)
he a= he a0;
end
i ( y>1e-15 && x<-1e-15)
he a=pi- he a0;
end
i ( y<1e-15 && x<-1e-15)
he a=pi+ he a0;
end
i ( y<1e-15 && x>1e-15)
he a=2*pi- he a0;
end
end
end
112
acphi2=dospi*(c2* I(3)+d2* I(2));
acphi3=dospi*(c3* I(3)+d3* I(2));
i ila(1)>1 && ila(1)<Nyincog+2 && columna(1)>1
Ma ( 0,columnaphi1)=Ma ( 0,columnaphi1)+ acphi1;
End
i ila( iangulo+1)>1 && ila( iangulo+1)<Nyincog+2 && columna( iangulo+1)>1
Ma ( 0,columnaphi2)=Ma ( 0,columnaphi2)+ acphi2;
End
i ila( iangulo+2)>1 && ila( iangulo+2)<Nyincog+2 && columna( iangulo+2)>1
Ma ( 0,columnaphi3)=Ma ( 0,columnaphi3)+ acphi3;
end
end