Bifurcación de Hopf y comportamiento
no lineal de conductos flexibles bajo la
influencia de la gravedad
Javier Rivero Rodríguez
Ingeniería Aerospacial y Mecánica de Fluidos
Universidad de Sevilla
Trabajo Fin de Máster
Máster Diseno Avanzado de Sistemas Mecánicos
5 de Diciembre de 2012
A mis padres...
Resumen
En este trabajo se estudia la pérdida de estabilidad de una manguera empotrada-
libre bajo un campo gravitatorio con un caudal impuesto a través de una bifurcación
tipo Hopf, y el efecto que tiene sobre ésta la orientación respecto a la gravedad. El
conducto se supondrá muy esbelto e inextensible, de manera que la deformación
axial sea despreciable frente a la de flexión. Se admite una ley de comportamiento
a flexión lineal, en la cual curvatura y flexión sean proporcionales. Se modela el
sistema con la teoría de Kirchhof para medios esbeltos, evitando así el uso de coorde-
nadas cartesianas. Se explica el efecto de la gravedad según esta modifica la solución
de equilibrio, flectando el conducto e imponiendo una tracción o compresión según
sea la orientación. En la realización de este trabajo, se han usado conocimientos
adquiridos durante el Máster tales como estabilidad de sistemas dinámicos, métodos
de continuación, parámetros de Euler y cálculo no lineal de estructuras esbeltas.
Índice
Índice iii
Lista de Figuras v
1 Introducción 1
2 Formulación del problema 4
2.1 Ecuaciones de movimiento . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
2.1.1 Cinemática . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
2.1.2 Dinámica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
2.2 Linealización de las ecuaciones . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
2.3 Transitorio . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.4 Órbitas periódicas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
3 Resultados 14
3.1 Estabilidad . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
3.1.1 Caso de gravedad despreciable . . . . . . . . . . . . . . . . . . . . . . . . 14
3.1.2 Conducto en posición vertical . . . . . . . . . . . . . . . . . . . . . . . . 15
3.1.2.1 Descarga hacia abajo . . . . . . . . . . . . . . . . . . . . . . . . 16
3.1.2.2 Descarga hacia arriba . . . . . . . . . . . . . . . . . . . . . . . 16
3.1.3 Conducto en posición horizontal . . . . . . . . . . . . . . . . . . . . . . . 16
3.1.4 Conducto en posición oblicua . . . . . . . . . . . . . . . . . . . . . . . . 17
3.2 Comportamiento inestable no lineal . . . . . . . . . . . . . . . . . . . . . . . . . 19
3.2.1 Transitorio . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
3.2.2 Órbitas periódicas no lineales . . . . . . . . . . . . . . . . . . . . . . . . 19
4 Conclusiones y trabajos futuros 21
4.1 Chorros electrificados . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
4.1.1 Ecuaciones . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
4.1.2 Resultados experimentales y numéricos . . . . . . . . . . . . . . . . . . . 23
4.2 Conclusiones . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
A Problema de Valores de Contorno 25
A.1 Método de shooting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
iii
ÍNDICE
A.2 Diferencias finitas centradas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
Bibliografía 28
iv
Lista de Figuras
2.1 Triedro intrínseco . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
2.2 Balances de cantidad de movimiento . . . . . . . . . . . . . . . . . . . . . . . . 6
2.3 Autovalores en bifuración de Hopf . . . . . . . . . . . . . . . . . . . . . . . . . . 9
3.1 Transición a inestable para gravedad despreciable . . . . . . . . . . . . . . . . . 15
3.2 Transición a inestable para conductos verticales hacia abajo . . . . . . . . . . . 16
3.3 Transición a inestable para conductos verticales hacia arriba . . . . . . . . . . . 17
3.4 Equilibrio conducto horizontal . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
3.5 Transición a inestable para conductos horizontales . . . . . . . . . . . . . . . . . 18
3.6 Transición a inestable para conductos oblicuos entre horizontal y hacia arriba . . 18
3.7 Transición a inestable para conductos oblicuos entre horizontal y hacia abajo . . 19
3.8 Transitorio estable . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
3.9 Transitorio inestable . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
3.10 Órbita periódica G = 0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
4.1 Inestabilidades en chorros cargados eléctricamente . . . . . . . . . . . . . . . . . 21
4.2 Electrospinning numérico y experimental . . . . . . . . . . . . . . . . . . . . . . 23
v
Capítulo 1
Introducción
La interacción fluido-estructura es un problema muy interesante para la ingeniera ya que explica
muchos fenómenos de inestabilidades como flameo de alas de avin, colapso de estructuras como
en el conocido caso del puente de Tacoma, o las vibraciones violenta que experimenta una
manguera de jardn cuando descarga caudal suficientemente elevado de agua. Estos problemas,
bien conocidos por todos, tienen una caracterstica comn y es que desarrollan inestabilidades
con carcter oscilatorio En este trabajo nos vamos a centrar en el problema de las inestabilidades
que puede presentar un conducto flexible que puede moverse en tres dimensiones cuando está
sometido sometido a la acción de la gravedad y por el que circula un cierto caudal de agua.
El problema de inestabilidades en conductos ha sido ampliamente estudiado, aparentemente
desde alrededor de 1963. Benjamin se puede considerar el primer investigador que resolvió prob-
lema análogo que exhibe el mismo comportamiento Benjamin [1961a,b], aunque no considerara
un sistema continuo sino dos tuberas articuladas, . Tres anos más tarde, Roth [1964] estableció
las ecuaciones para una multitud de casos, desde tuberas de sección variable hasta conductos
apoyados en un lecho el ástico, pero estas ecuaciones fueron planteadas en desplazamientos
transversales, por lo que están restringidas sólo a casos débilmente no lineales. Posteriormente,
otros autores establecieron de nuevo las ecuaciones, entre ellos, G. X. Li [1994] que realizó una
comparativa entre las formulaciones realizadas anteriormente y demostró la incoherencia de
algunas de ellas. En este último trabajo se discuten las formas de obtener las ecuaciones, tanto
mediante el método de Newton o por el de Lagrange, así como las coordenadas apropiadas y la
condiciíon de inextensibilidad que permite despreciar la deformación axial. Trabajos anteriores
y referenciados por G. X. Li [1994] haban perseguido diversos objetivos entre los que cabe citar
los trabajo P.R. Sethna and Bajaj. [1979, 1980], en los que se determina las condiciones para
la aparición de la inestabilidad con un caudal de líquido impuesto y constante en el tiempo,
así como el trabajo G. Herrman [1980] donde se considera qu ocurre bajo una diferencia de
presiones impuesta entre la toma y la sección de salida del conducto. En estos trabajos las
condiciones de contorno impuestas a la tubería son las que se considerarn también en este
trabajo son las de un extremo empotrado y el otro libre, siendo este último el correspondiente
a la seccin de salida del líquido. Otro trabajo interesante en la materia se debe a J.S.Jensen
[1997] en la que se estudia cmo forzando la vibracin en un conducto empotrado-libre con toma
1
1. Introducción
a un depósito, se tiene una bomba sin partes móviles que pueden resultar de pueden tener gran
interés práctico en aplicaciones médicas. Una referencia clásica para este tipo de problemas así
como otros concernientes a placas y otros casos concernientes a elementos bidimensionales se
debe a Paidoussis [a,b], quien posee una gran experiencia en este campo.
En este trabajo se estudia la pérdida de estabilidad de una manguera empotrada-libre bajo
un campo gravitatorio con un caudal impuesto a través de una bifurcación tipo Hopf, y el
efecto que tiene sobre ésta la orientación respecto a la gravedad. El conducto se supondrá
muy esbelto e inextensible, de manera que la deformación axial sea despreciable frente a la de
flexión. Se admite una ley de comportamiento a flexión lineal, en la cual curvatura y flexión
sean proporcionales.
Se han usado conocimientos adquiridos en el Máster tales como aquellos de estabilidad
de sistemas dinámicos, métodos de continuación, parámetros de Euler y cálculo no lineal de
estructuras esbeltas.
En el capítulo 2 se obtienen las ecuaciones generales que gobiernan la cinemática y la
dinámica del problema. Las ecuaciones de movimiento se obtienen aplicando un balance de
fuerzas y momentos a una rebanada del conducto y se definen los parámetros y variables
adimensionales. La descripciń cinemática del problema se ha tratado de manera diferente para
incluir el efecto tridimensional. En lugar de usar los ángulos girado por la linea media de
la sección (2D), se han usado modelos de Kirchhof’s rod como ya se han aplicado al mismo
campo M.T.Mathew [2010]. En este capítulo se obtienen también las ecuaciones linealizadas
del problema, el esquema numérico que se utiliza para resolver el transitorio y las ecuaciones
para las órbitas periódicas que se alcanzan cuando el sistema es inestable. Finalmente, se aplica
un balance de energía en un volumen de control que contiene al sistema completo obteniéndose
así la ecuación de la energía que permite analizar la transferencia de energía en el sistema.
En el capítulo 3 se exponen los resultados de estabilidad para un caso de referencia, sea
ste en ausencia de gravedad. Se resuelven igualmente los casos de conducto vertical, en ambas
orientaciones. Estos estn incluidos en el estudio bidimensional de Paidousss. Se comparan
los resultados bidimensionales y tridimensionales. Se estudia el caso de conducto horizontal,
as como configuraciones intermedias barriendo as cualquier orientacin. Se simulan algunos
transitorios, donde se observan dos comportamiento principales: estable y flameo. Se obtienen
tambin las rbitas peridicas para casos con forma de equilibrio rectas.
En el capítulo 4 se escriben las ecuaciones para un chorro cargado eléctricamente y se
obtiene la bifurcación de Hopf, observándose comportamientos cualitativamente similares a los
observados experimentalmente en el laboratorio, G.Riboux [2011].
En el capítulo 4 se expone una breve introducción al problema del electrospinning, de
interés en el departamento, y se fórmula el problema. Los métodos numéricos desarrollados
para obtener las bifurcaciones de Hopf en conductos tienen aplicación en dicho problema que
presenta geometrías muy parecidas a la del conducto. En el mismo capítulo se expone las
conclusiones a las que se llega del presente trabajo.
En el apéndice A, se dan detalles sobre la resolución numérica de las ecuaciones de valores
2
1. Introducción
de contorno por el método de shooting y de diferencias finitas.
3
Capítulo 2
Formulación del problema
En este capítulo se deducen las ecuaciones que gobiernan el problema, aplicando las leyes de
Newton a cada rebanada infinitesimal del conducto y del fluido y usando el modelo de Kirchhof
(J.Valverde et al. [2006]), basado en el triedro de Cosserat, para caracterizar la geometría del
problema. Posteriormente se adimensionaliza n las ecuaciones y se linealizan para estudiar la
pérdida de estabilidad por bifurcación de Hopf (ver Strogatz). Por último se analizarán también
la evolución transitoria hacia las órbitas periódicas que se produce después de dicha bifurcación.
2.1 Ecuaciones de movimiento
En esta sección se hará primero una descripción cinmática del problema en términos del vector
posición respecto de la sección de entrada del conducto r(s, t), donde s la longitud del tramo
de conducto hasta el punto considerado, la velocidad del conducto v(s, t), su velocidad an-
gular y su curvatura, ω(s, t) y κ(s, t) respectivamente, y la velocidad impuesta del fluido, u.
A continuación, en la formulación de la dinámica del problema, se considerarán también los
esfuerzos integrados en cada sección, caracterizados por la fuerza f (s, t) y el momento m(s, t),
así como las condiciones de contorno y la ley de comportamiento lineal usada. Por último,
una vez adimensionalizadas las ecuaciones, se linealizan éstas y se escriben en términos de los
parámetros adimensionales del problema: la relación masa de fluido a masa total del sistema,
β, la velocidad adimensional impuesta del fluido, U , y gravedad adimensional G.
2.1.1 Cinemática
Se usará una base vectorial intrínseca di(s, t) (figure 2.1) anclada al conducto y tal que
d3(s, t) = ∂sr(s, t). (2.1)
Las variaciones espaciales y temporales de los vectores de esta base están dados por la curvatura
y la velocidad angular del conducto, κ y ω, mediante las relaciones:
∂sdi = κ × di, (2.2a)
4
2. Formulación del problema
∂tdi = ω × di. (2.2b)
Se ha comprobado que esta base es más conveniente que el sistema de coordenadas carte-
sianas xi (base ai) tradicionalmente usado para este tipo de problemas en la literatura, ya que
éstas complican mucho el tratamiento numérico y oscurecen la formulación del problema. La
velocidad de cada sección del conducto del conducto es
v(s, t) = ∂tr(s, t), (2.3)
y la velocidad absoluta del fluido es v + ud3.
Figure 2.1: Triedro intrínseco
Las derivadas derivadas espaciales y temporales de un vector genérico expresado en la base
intrínseca, b =∑
i bidi, son:
∂sb =∑
i
∂s(bidi) =∑
i
(∂sbi)di +∑
i
bi∂sdi =∑
i
(∂sbi)di +∑
i
bi∂sdi = b′ + κ × b (2.4)
y, anaĺogamente,
∂tb = b + ω × b, (2.5)
donde, para simplificar la notación, se han definido los vectores, b′ ≡ ∑
i(∂sbi)di y b ≡∑
i(∂tbi)di.
Obsérvese que las variables v, κ y ω deben cumplir las siguientes ecuaciones, o condiciones
de compatibilidad, deducidas de la igualdad de derivadas cruzadas del vector de posición r(s, t)
y del triedro di(s, t):
∂s∂tr = ∂t∂sr : v′ + κ × v = ω × d3, (2.6a)
∂s∂tdi = ∂t∂sdi : ω′ = κ + ω × κ. (2.6b)
Finalmente, para el estudio de la evolución no lineal del problema ha resultado muy útil
el uso de los parámetros de Euler q0, q1, q2 y q3, en términos de los que la matriz de cosenos
directores de la base cartesiana respecto al triedro intrínseco, Dij(s, t) = di(s, t) ·aj , se escribe
5
2. Formulación del problema
como
D =
q20 + q2
1 − q22 − q2
3 2(q0q3 + q1q2) 2(−q0q2 + q1q3)
2(−q0q3 + q1q2) q20 − q2
1 + q22 − q2
3 2(q0q1 + q2q3)
2(q0q2 + q1q3) 2(−q0q1 + q2q3) q20 − q2
1 − q22 + q2
3
, (2.7)
donde los parámetros de Euler deben satisfacer q20 + q2
1 + q22 + q2
3 = 1.
A partir de las ecuaciones que relacionan las derivadas de los vectores del triedro intrínseco
con la curvatura y la velocidad angular 2.2, se obtienen las ecuaciones:
∂s
q0
q1
q2
q3
=12
−q1 −q2 −q3
q0 −q3 q2
q3 q0 −q1
−q2 q1 q0
κ1
κ2
κ3
, ∂t
q0
q1
q2
q3
=12
−q1 −q2 −q3
q0 −q3 q2
q3 q0 −q1
−q2 q1 q0
ω1
ω2
ω3
. (2.8)
Nótese que se conserva la ortonormalidad del triedro tal que se satisface la condición de nor-
malidad de los parámetros de Euler, (ver W.Arne [2010])
En cualquier caso, existe una alternativa al uso de parámetros de Euler. En efecto, debido
a que un vector cartesiano no varía con s,
∂sai = a′
i + κ × ai = 0, (2.9)
se puede obtener dicho vector en la base intrínseca integrando la ecuación diferencial
a′
i = −κ × ai (2.10)
sujeta a las condiciones iniciales que se obtienen de hacer coincidir en s = 0 del triedro di y la
base ai.
2.1.2 Dinámica
Las ecuaciones del movimiento se obtienen mediante un balance de cantidad de movimiento
linear y angular para una sección genérica del conducto como la mostrada en la figura 2.2. Sea
Ap, Ip, ρp son, respectivamente, el área, momento de inercia y densidad del conducto y Af , If ,
ρf las correpondientes magnitudes del fluido,
Figure 2.2: Balances de cantidad de movimiento
6
2. Formulación del problema
las ecuaciones de balance de cantidad de movimiento para el conducto se escriben como
∂sfp + ρp Ap g + ffp = ρp Ap ∂tv, (2.11a)
∂smp + d3 × fp + mfp = ρpIp · ∂tω, (2.11b)
siendo ffp y mfp las fuerzas y momentos del fluido sobre el conducto, mientras que las corre-
spondientes ecuaciones para el fluido son:
∂sff + ρf Af g + fpf = ρf Af dt(v + ud3), (2.12a)
∂smf + d3 × ff + mpf = ρf If · dt(ω + uκ), (2.12b)
con fpf y mpf las fuerzas y momentos del conducto sobre el fluido. Nótese que la aceleración
del fluido es la derivada sustancial (derivada siguiendo el movimiento relativo al conducto de
una partícula de fluido) de la velocidad, estando el operador derivada sustancial definido como:
dt = ∂t + u∂s, (2.13)
Sumando ahora las ecuaciones de balance del conducto 2.11 y las de fluido 2.12, reagrupando
los esfuerzos como f = ff + fp y m = mf + mp y considerando la ley de acción y reacción:
ffp + fpf = 0 y mfp + mpf = 0, las ecuaciones de la dinámica para el sistema resultan:
∂sf + g = ρf Af dt(v + ud3) + ρp Ap ∂tv (2.14a)
∂sm + d3 × f = ρf If · dt(ω + uκ) + ρpIp · ∂tω. (2.14b)
Para cerrar el problema es necesario hacer uso de una ley de comportamiento que relacione
los momentos flectores con la curvatura del conducto y la torsión del conducto. Como es usual,
se escogerá ley de comportamiento lineal:
m = EIκ, (2.15)
en el tensor I las componentes I1 = I2 son las correspondientes a flexión y I3 = 2I1 = 2I2 la
correspondiente a torsión. Obsérvese que se ha despreciado el momento debido a los esfuerzos de
viscosidad del fluido, puesto que, de acuerdo la ley constitutiva mf = 3µIfdiag([1, 1, 2/3]) ·∂sω
normalmente usada en W.Arne [2010], se tiene
mf
mp∼ ff
fp∼ µIf L/ω0
EIp/L<< 1, (2.16)
7
2. Formulación del problema
donde se han considerado tiempos característicos del orden del periodo natural de flexión, ω−10 ,
siendo ω0 ∼√
EI/(ρAL4). Por tanto, m ≈ mp.
Las ecuaciones !!!! deben resolverse sujetas a las condiciones de contorno correspondientes
a un conducto empotrado en s = 0 y libre en s = L :
v(0, t) = 0 ω(0, t) = 0, (2.17a)
f (L, t) = 0 κ(0, t) = 0. (2.17b)
Para el análisis y la resolución numérica de las ecuaciones anteriores es conveniente escribir-
las en términos adimensionales mediante las sustituciones
s → Ls, f → ρA(ω0L)2f , v → ω0Lv, ω → ω0ω, κ → κ/L (2.18)
donde L es la longitud del conducto, ρA...., I el momento de inercia de la sección del con-
ducto, y ω0 ∼√
EI/(ρAL4) la frecuencia característica natural de flexion, que proporcionan
las sigueintes ecuaciones adimensionales:
f
v
ω
κ
′
=
0 0 I 0
0 0 0 0
0 I 0 0
0 0 0 0
˙
f
v
ω
κ
+
−κ × f + ω × v + 2βUω × d3 + βU2κ × d3 − G
−κ × v + ω × d3
ω × κ
−d3 × f
. (2.19)
donde G = g/ω20L es la gravedad adimensional, β = ρfAf /(ρfAf +ρpAp) es el ratio de masas y
U = u/(ω0L) es la velocidad del fluido adimensional. Obsérvese que los términos de inercia en
la ecuación de momentos se pueden despreciar por ser del orden I/A L2 ≪ 1 para conductos y
UI/A L2 << 1 para el fluido, donde como veremos más adelante los valores típicos de U están
en el orden de 10.
Obsérvese que definiendo y = [f , v, ω, κ]T el sistema de ecuaciones anterior se puede
escribir en forma compacta como
y′ = M y + F (y). (2.20)
Para un conducto típico de 1m de sección circular, de diámetro exterior φi = 10mm,
diámetro interior φe = 20mm, de silicona tiene los siguientes valores adimensionales,
ω0 ≈ 0.36Hz, v0 ≈ 0.36m/s, f0 ≈ 0.049N, m0 ≈ 0.049Nm, v0 ≈ 0.36m/s, β ≈ 0.23, G ≈ 72.
(2.21)
8
2. Formulación del problema
2.2 Linealización de las ecuaciones
Como es usual, para la estudiar la estabilidad lineal de una solución de equilibrio f0(s),
v0, ω0, κ0 del sistema 2.19 se perturba dicha solución en la forma
f = f0 + ǫf1eiΩt + O(ǫ2) (2.22)
v = v0 + ǫv1eiΩt + O(ǫ2) (2.23)
ω = ω0 + ǫω1eiΩt + O(ǫ2) (2.24)
κ = κ0 + ǫκ1eiΩt + O(ǫ2), (2.25)
ya que en este trabajo estamos interesados en la pérdida de estabilidad que se desarrolla cuando
los autovalores de la función cortan al eje imaginario desde la parte real negativa hasta hacia
la positiva (±iΩ), siendo Ω la frecuencia de la oscilación. Este tipo de inestabilidad se conoce
como bifurcación de Hopf (Strogatz).
ℑ(λ)
ℜ(λ)
+i Ω
−i Ω
Figure 2.3: Autovalores en bifuración de Hopf
Para linearlizar las ecuaciones 2.19 es necesario perturbar los vectores del triedro intrínseco
y un vector expresado en dicho triedro. En primer lugar, la derivada de un vector del triedro
se puede escribir
di = d(0)i + ǫd
(1)i + O(ǫ2), (2.26)
y haciendo unso de la condición de ortonormalidad
di · dj = d(0)i · d(0)
j + ǫ (d(1)i · d(0)
j + d(1)j · d(0)
i ) + O(ǫ2) = δij + O(ǫ2), (2.27)
se deduce
d(1)i = α × d
(0)i . (2.28)
Para obtener una expresión que defina α, se deriva espacialmente el triedro, obteniéndose
d′
i = (κ0 + ǫκ1 + O(ǫ)) × (d(0)i + ǫα × d
(0)i + O(ǫ)) (2.29)
= κ0 × d(0)i + ǫ[κ0 × (α × d
(0)i ) + κ1 × d
(0)i ] (2.30)
9
2. Formulación del problema
y
d′
i = d(0)′
i + ǫ(α × d(0)i )′ = (κ0 + ǫκ1) × d
(0)i + ǫ(α′ × d
(0)i + α × (κ0 × d
(0)i )) (2.31)
= κ0 × d(0)i + ǫ[α′ × di + α × (κ0 × d
(0)i )] (2.32)
por tanto, desarrollando los productos triples y teniendo en cuenta que se cumple para i = 1, 2, 3,
se obtiene
α′ = κ1 + α × κ0. (2.33)
En segundo lugar, un vector b expresado en el triedro intrínseco se perturba como
b = b(0)i d0
i + ǫb(1)i d0
i = b0 + ǫb1. (2.34)
Una vez conocido como perturbar un vector y el triedro, se introduce la perturbación 2.25
en las ecuaciones 2.19, de forma que se obtiene las ecuaciones de equilibrio, orden ǫ0:
f ′
0 + κ0 × f0 + g = βU2κ × d3, (2.35a)
κ′
0 + d3 × f0 = 0, (2.35b)
donde se omiten v y ω por ser idénticamente 0. Y las ecuaciones para las perturbaciones,
orden ǫ1, son:
f ′
1 + κ1 × f0 + κ0 × f1 = iΩv1 + 2βUω1 × d3 + βU2[κ1 × d3 + κ0 × (α × d3)], (2.36a)
κ′
1 + d3 × f1 + (α × d3) × f0 = 0, (2.36b)
v′
1 + κ0 × v1 + κ1 × v0 = ω1 × d3 (2.36c)
ω′
1 = iΩκ1 + ω1 × κ0 (2.36d)
Comparando 2.33 y 2.36d, se obtiene que ω = iΩα, por tanto el problema se puede escribir
como
f ′
1 + κ1 × f0 + κ0 × f1 = iΩv1 + 2iΩβUα × d3 + βU2[κ1 × d3 + κ0 × (α × d3)], (2.37a)
κ′
1 + d3 × f1 + (α × d3) × f0 = 0, (2.37b)
v′
1 + κ0 × v1 + κ1 × v0 = iΩα × d3 (2.37c)
10
2. Formulación del problema
α′ = κ1 + α × κ0 (2.37d)
Las ecuaciones 2.37 se pueden escribir en forma compacta como
z′
1 = M(µ)z1 + K(f0,κ0;µ)z1, (2.38)
donde se han reagrupado las variables como z1 = [f1,v1,α,κ1]T y los parámetros en µ. Sep-
arando partes real e imaginaria del vector de variables, z1 = zr1 + izi
1 el problema se escribe
como
zr1
zi1
′
=
K iΩM
−iΩ K
yr1
yi1
. (2.39)
junto con las condiciones de contorno
g[z1(0), z1(1);µ] = 0, (2.40)
separados en parte real e imaginaria. Hay que considerar una condición de normalización (dos
realmente, una para la parte real y otra para la parte imaginaria) ya que tenemos dos parámetros
libres Ω y una función escalar de µ. Es importante que esta condición de normalización no sea
incompatible con alguna condición de contorno.
Para la resolución de estas ecuaciones se ha implementado en MATLAB un programa que
están expresadas en forma general en el capítulo 6 del libro Seydel.
Nota numérica: En nuestro caso se fija G y β para calcular U y Ω. Posteriormente, para
una G fijada, se obtiene una curva en el plano β y U , partiendo de la solución anteriormente
calculada y continuando por el método de pseudo-arclength. Este método consiste en calcular
la curva µ = µ(ξ) sujeto a unas condiciones, en este caso, la curva en el espacio paramétrico
donde ocurre la bifurcación de Hopf.
2.3 Transitorio
En esta sección se expone el método numérico con el que se integra 2.19 en el tiempo para
obtener el comportamiento transitorio del problema. Debido a que la matriz M no es de rango
completo, es necesario integrar en el tiempo de forma implícita. Se discretizan las ecuaciones
2.20 con diferencias finitas centradas cortas tanto en el tiempo como en el espacio. Las derivadas
espaciales temporales y el vector de variables queda discretizado como
y′(sn+1/2, tm+1/2) =y(sn+1, tm+1) + y(sn+1, tm) − y(sn, tm+1) − y(sn, tm)
2(sn+1 − sn)(2.41a)
y(sn+1/2, tm+1/2) =y(sn+1, tm+1) − y(sn+1, tm) + y(sn, tm+1) − y(sn, tm)
2(tm+1 − tm)(2.41b)
11
2. Formulación del problema
y(sn+1/2, tm+1/2) =y(sn+1, tm+1) + y(sn+1, tm) + y(sn, tm+1) + y(sn, tm)
4(2.41c)
Que sustituyendo en el sistema 2.19 se obtiene una función de las variables en tm y tm+1
Ψ[y(sn+1/2, tm+1/2)] ≡ y′(sn+1/2, tm+1/2) − M y(sn+1/2, tm+1/2) − f [y(sn+1/2, tm+1/2)] = 0.
(2.42)
Entonces, conocido y(sn, tm) para n = 1, .., N , y reagrupando para tm+1
Y (tm+1) =
y(s1, tm+1)...
y(sn, tm+1)...
y(sN , tm+1)
(2.43)
se obtiene una expresión para calcular las variables en el paso de tiempo siguiente, tm+1,
Ψ[yν(sn+1/2, tm+1/2)] = Ψ[Y ν+1(tm+1)] = 0, que se puede resolver con el método de Newton-
Raphson con la siguiente expresión
Ψ[Y ν+1(tm+1)] ≈ Ψ[Y ν(tm+1)] +∂Ψ[Y ν(tm+1)]
∂Y ν(tm+1)∆Y ν(tm+1) = 0, (2.44)
y actualizando Y ν+1(tm+1) = Y ν(tm+1) + ∆Y ν(tm+1).
2.4 Órbitas periódicas
Después de que el sistema pierde la estabilidad, se observará que en el transitorio se desarrolla
hacia una órbitas periódicas. Entonces, se buscan éstas aprovechando la simetría del sistema,
es decir, tales que el conducto giran alrededor de la dirección a3 con una velocidad Ω. Se toma
un sistema de referencia giratorio que vea la manguera de forma estacionaria,
vΩ = v − Ωa3 × r (2.45a)
ωΩ = ω − Ωa3 (2.45b)
Imponiendo la estacionareidad, vΩ = 0 y ωΩ = 0, las ecuaciones de compatibilidad 2.6
se cumplen automáticamente. Los vectores de posición y de velocidad se expresan en la base
rotatoria y se usan los parámetros de Euler para relacionar ambas bases.
La ecuación 2.1 expresando el vector d3 en función de los parámetros de Euler
r = d3(qi), (2.46)
12
2. Formulación del problema
y la expresión 2.8 definen la geometría del sistema. Las ecuaciones de conservación de cantidad
de movimiento (2.14), escritas en las nuevas variables
f ′ = −κ × f + D [Ω × (Ω × r)] + 2βU(D Ω) × d3 + βU2κ × d3 − Dg (2.47a)
κ′ + d3 × f = 0, (2.47b)
y las condiciones de contorno cierran el problema.
Nota numérica: para resolver este problema de valores de contorno se usan métodos itera-
tivos en los que se requiere un ’initial guess’, los cuales se pueden tomar las funciones calculadas
en el problema linealizado que convergen con más facilidad e ir adaptando la solución paso a
paso hasta el conjunto de valores de β y U que se quiera calcular.
13
Capítulo 3
Resultados
3.1 Estabilidad
Se resuelve la bifurcación de Hopf obteniéndose como resultado el diagrama de estabilidad, los
autovalores y autofunciones. Se ha considerado conveniente, expresar los diagramas de estabil-
idad en el plano βU2-βU , debido a que en las ecuaciones de movimiento β y U solo aparecen
formando esos dos grupos. Recuérdese que βU2 está relacionado con las fuerza centrífuga y βU
con la fuerza de Coriolis.
Se resuelven los casos de gravedad despreciable (conductos muy rígidos, G ≈ 0) y con
gravedad importante (G > 0) para conductos orientados paralelos a la gravedad, es decir,
en posición vertical considerando las dos posibles orientaciones posibles. Estos casos están
resueltos en la literatura por Paidoussis y se utilizan para testar el programa, así como de
referencia para las variaciones que se estudian. Posteriormente, se resuelven los casos en los
que el conducto no está alineado con la gravedad, de manera que el conducto se encuentra
deformado en la posición de equilibrio. Esto tiene consecuencias en la estabilidad del sistema
y en los modos oscilatorios inestables.
También se ha resuelto el transitorio del problema para dos condiciones, una estable y otra
inestable. En la inestable se observa como el sistema alcanza una órbita periódica, la cual se
estudia para casos en que la forma de equilibrio es recta.
3.1.1 Caso de gravedad despreciable
El caso sin gravedad se tomar como un caso de referencia. Ya que está reportado en la literatura
se utilizar para testar el programa y se amplia el caso al estudio del sistema con libertad de
movimiento en 3D.
En 3.1(a) se observa que existe una región donde β > 1, lo cual no tiene sentido físico.
Observamos que la responsable de la inestabilidad es la fuerza centrífuga, βU2, mientras la
fuerzas de Coriolis estabilizan el sistema, βU . Se observa igualmente que cuando aumenta el
flujo, ya sea βU2 (como aparece en 3.1(b)) o βU , las oscilaciones del sistema se vuelven más
lentas al igual que los modos tienen longitudes de onda más cortas (3.1(c)-3.1(d)). En estas
14
3. Resultados
0 75 150 225 300 375 4500
4
8
12
16
20
βU2
βU
Sta
ble Flutter
β > 1
A
B
(a) Diagrama de estabilidad
0 75 150 225 300 375 4500
30
60
90
βU2
Ω
A
B
(b) Autovalor Ω
−1−0.8
−0.6−0.4
−0.20
−1
−0.5
0
0.5
1−0.5
0
0.5
xy
z
(c) Modo inestable para valores en punto A
−1−0.8
−0.6−0.4
−0.20
−1
−0.5
0
0.5
1−0.5
0
0.5
xy
z
(d) Modo inestable para valores en punto B
Figure 3.1: Transición a inestable para gravedad despreciable
figuras se representas las vista tridimensional del modo y sus proyecciones en los tres planos.
Se representa solo la parte real del modo, ya que al ser existir simentría de resolución, la parte
imaginaria es igual que la real desfasada 90o. Nótese que las perturbaciones son por definición
pequenas y se ha utilizado las ecuaciones linealizadas, pero se han representados modos con
amplitudes mayores para apreciarlos mejor.
3.1.2 Conducto en posición vertical
Se estudian dos casos de orientación vertical, descargando el chorro en sentidos opuestos, es
decir, hacia abajo y hacia arriba. En estos casos, la forma de equilibrio es recta y el conducto
está sometido a axiles de tracción y de compresión respectivamente. Como se verá, éste rigidiza
o flexibiliza el sistema de forma que se estabiliza o desestabiliza, respectivamente.
15
3. Resultados
3.1.2.1 Descarga hacia abajo
En el caso de que el conducto descarga hacia abajo, observamos que la gravedad estabiliza
(figura 3.2(a)) y disminuye la frecuencia de oscilación ligeramente (figura 3.2(b)). También se
endereza el modo en la zona cercana al empotramiento, donde la tracción debida a la gravedad
son más elevados que en el resto del conducto. Este efecto es ínfimo y se omite la representación
de los modos ya que esta diferencia no es muy importante.
0 75 150 2250
3
6
9
12
15
βU2
βU
Sta
ble
Flutter
β > 1
(a) Diagrama de estabilidad
0 75 150 2250
10
20
30
40
50
βU2
Ω
(b) Autovalor Ω
Figure 3.2: Transición a inestable para conductos verticales descargando hacia abajo, G = 0( ), G = 15 ( ), G = 30 ( )
3.1.2.2 Descarga hacia arriba
En el caso de gravedad vertical hacia arriba, el efecto de la gravedad es el contrario al anterior
(figura 3.3). En este trabajo no se ha estudiado el pandeo debido a que ya está reportado en
la literatura, y resolvemos este caso para tenerlo como referencia en orientaciones oblicuas.
3.1.3 Conducto en posición horizontal
Cuando el conducto se encuentra en horizontal, hay una flexión del conducto en el caso de
equilibrio que afectará a la estabilidad del sistema. De las ecuaciones 2.35, se deduce que dicho
equilibrio sólo depende de G y de βU2. En las figuras 3.4 se representa dicho efecto. Se observa
que la gravedad tiende a flectar el conducto, al igual que lo tracciona al alinearse parcialmente
el conducto y la gravedad. El fluido tiende a enderezar el conducto, compensando el efecto de
la gravedad.
En 3.5(a) podemos observar como el efecto de la gravedad es estabilizante, debido al efecto
de la curvatura del conducto y a la tracción que como habíamos visto anteriormente, es esta-
bilizante. Este efecto, es importante para valores bajos de βU2 ya que para valores mayores
16
3. Resultados
0 75 150 2250
3
6
9
12
15
βU2
βU
Sta
ble
Flutter
β > 1
(a) Diagrama de estabilidad
0 75 150 2250
10
20
30
40
50
βU2
Ω
(b) Autovalor Ω
Figure 3.3: Transición a inestable para conductos verticales descargando hacia arriba, G = 0( ), G = 15 ( ), G = 30 ( )
G = 0
G = 30
(a) Efecto de la gravedad, βU2 = 0
βU2 = 150
βU2 = 0
(b) Efecto del flujo, G = 30
Figure 3.4: Equilibrio de conductos horizontales
el efecto de la gravedad en el equilibrio es inferior según éste aumenta. Los autovalores Ω
disminuyen con G, figura 3.5(b).
3.1.4 Conducto en posición oblicua
Una vez estudiado el caso horizontal, se estudian orientaciones intermedias entre el caso hori-
zontal y cada uno de los verticales. En primer lugar, se expone el comportamiento del caso con
orientación entre la horizontal y la vertical hacia arriba. En la figura 3.6 se representa como una
orientación intermedia se comporta similar al conducto vertical hacia arriba para flujos βU2
altos. Esto se debe a que para flujos altos, el conducto está en una posición muy enderezada
por lo que la componente de la gravedad transversal al conducto juga un papel irrelevante al
ser compensado por el flujo. Sin embargo, el comportamiento es similar al conducto horizontal
para flujos pequenos, debido a que el papel que juega la componente transversal de la gravedad
17
3. Resultados
0 75 150 2250
3
6
9
12
15
βU2
βU
Sta
ble
Flutter
β > 1
(a) Diagrama de estabilidad
0 75 150 2250
10
20
30
40
50
βU2
Ω
(b) Autovalor Ω
Figure 3.5: Transición a inestable para conductos horizontales, G = 0 ( ), G = 15 ( ),G = 30 ( )
es más importante respecto al flujo y la deformada se aleja de la recta. La compresión a la que
está sometida el conducto, ayuda a que el sistema se aleje de la posición recta.
0 75 150 2250
3
6
9
12
15
βU2
βU
Sta
ble
Flutter
β > 1
(a) Diagrama de estabilidad
0 75 150 2250
10
20
30
40
50
βU2
Ω
(b) Autovalor Ω
Figure 3.6: Transición a inestable para conductos oblicuos entre horizontal y hacia arriba,G = 0 ( ), G = 30 horizontal ( ), G = 30 hacia arriba ( ), G = 30 oblicuo ( )
En segundo lugar, se expone el comportamiento del caso con orientación entre la horizontal
y la vertical hacia abajo. En la figura 3.7 se representa como una orientación intermedia
se comporta similar al comportamiento para un conducto vertical hacia abajo. Para flujos
pequenos, el comportamiento en horizontal y vertical hacia abajo es similar, pero no lo es así
para βU2 altos debido, otra vez al efecto enderezador del flujo. Al girar el conducto desde
una posición horizontal a una vertical hacia abajo, se ayuda a enderezar el conducto y por
18
3. Resultados
tanto el comportamiento de un conducto en posición intermedia tiende rápidamente a aquel
del conducto vertical hacia abajo.
Es importante remarcar que el efecto que tiene la gravedad se debe exclusivamente a la
modificación de la forma de equilibrio, tanto a la parte de momentos y cortantes que determinan
la forma de equilibrio como a la parte de los axiles que rigidizan o flexibilizan el conducto de
igual manera que la tensión en las cuerdas de guitarra.
0 75 150 2250
3
6
9
12
15
βU2
βU
Sta
ble
Flutter
β > 1
(a) Diagrama de estabilidad
0 75 150 2250
10
20
30
40
50
βU2
Ω
(b) Autovalor Ω
Figure 3.7: Transición a inestable para conductos oblicuos entre horizontal y hacia abajo, G = 0( ), G = 30 horizontal ( ), G = 30 hacia arriba ( ), G = 30 oblicuo ( )
3.2 Comportamiento inestable no lineal
En esta sección se estudia el transitorio en dos casos, estable e inestable, y las órbitas periódicas
que se alcanzan en casos con forma de equilibrio recta.
3.2.1 Transitorio
En la figura 3.8 se representa un transitorio en condiciones estables en ausencia de gravedad,
observándose como una perturbación es amortigua.
En la figura 3.9 se representa un transitorio en condiciones inestables en ausencia de
gravedad, observándose como una perturbación provoca que el sistema vaya a una órbita per-
iódica.
3.2.2 Órbitas periódicas no lineales
Las órbitas periódicas no lineales se estudian para casos en que el conducto en equilibrio es
recto, en caso contrario no existiría simetría que se ha utilizado para obtener las ecuaciones
19
3. Resultados
Figure 3.8: Transitorio estable
Figure 3.9: Transitorio inestable
(2.46-2.8-2.47). En la figura 3.10 se representan los valores para los que se calcula la órbita
periódica, y la forma del conducto en varios inestantes de tiempo.
0 75 150 225 300 375 4500
4
8
12
16
20
βU2
βU
Sta
ble Flutter
β > 1
Figure 3.10: Órbita periódica G = 0
20
Capítulo 4
Conclusiones y trabajos futuros
4.1 Chorros electrificados
Un método novedoso y barato para producir fibras con diámetros tan pequenos como 100nm
está relacionado con un proceso electrohidromecánico llamado electrospinning. Una aguja de
diámetro interior ∼ 100µm conectada a un potencial del orden de 1kV/cm por la que se bombea
un líquido (< 1ml/h), un cristal líquido, un polímero, emulsiones, etc, hacia un contraelectrodo
exhibe una diversidad de comportamientos aun sin explorar completamente hoy día y está bajo
una intensa investigación. Existen diversos modos según el tipo de inestabilidad que predomine:
dripping (figura 4.1(a)), jetting(figura 4.1(b)) y electrospray (figura 4.1(c)) y electrospinning
(figura 4.2).
(a) Dripping (b) Jetting
(c) Electrospray (d) Electrospinning
Figure 4.1: Inestabilidades en chorros cargados eléctricamente
21
4. Conclusiones y trabajos futuros
En el caso del electrospinning, la inestabilidad predominante es no-axilsimétrica debida a
la repulsión lateral de las cargas y es conocida como whipping. Esta peculiaridad hace que el
problema no pueda tratarse como el resto de casos en los que existe simetría de revolución.
Como podemos observar en la figura 4.2, la forma que exhibe el whipping presenta bastante
similitud con ’el problema de la manguera’. Por ello, se ha empezado con un problema más
simple y más conocido que permite desarrollar herramientas numéricas para capturar la forma
del chorro o manguera, teniendo en cuenta que ambos tienen una longitud predominante.
El primer estudio del problema se centro en el cono de Taylor que fue resuelto en Taylor
[1964]. V.M.Entov [1984] es el primero en modelar la dinámica de un chorro no recto, conc-
retamente inestabilidades no axilsimétricas de origen aerodinámico. Posteriormente, el mismo
autor ha modelado el comportamiento del electrospinning para polimeros viscoelásticos en co-
ordenadas cartesianas, A.L.Yarin [2000, 2001]. Otros estudios de dinámica de chorros no rectos
son N.M.Ribe [2004, 2006] que consideran el problema de ’coiling’ en líquidos viscosos como
la miel, o el proceso de fabricación conocido como ’rotational spinning’ en Marheineke and
Wegener [2009]; W.Arne [2010]; W.Arne. Un trabajo analítico que estudia inestabilidades en
chorros infinitos cargados eléctricamente, que también consideran las inestabilidades laterales
que ocurren en el electrospinning, es M.M.Hohman [2001a] completado con un estudio experi-
mental M.M.Hohman [2001b]. F.J.Higuera resuelve el problema eléctrico y la forma del chorro
recto, lo que para nosotros es la solución de equilibrio. Actualmente estamos aplicando las
herramientas desarrolladas para el problema del conducto al fenómeno del electrospinning.
4.1.1 Ecuaciones
Las ecuaciones del sistema son las mismas que para el conducto pero eliminando los términos
de inercia del conducto y permitiendo variaciones del radio de la sección,
v′ + κ × v = ω × d3, (4.1a)
ω′ = κ + ω × κ, (4.1b)
∂sff + ρf Af g + τe = ρ A dt(v + ud3), (4.1c)
∂smf + d3 × ff = ρ I · dt(ω + uκ), (4.1d)
con conservación de masa uA = Q. Las leyes de comportamiento son leyes viscosas,
f · d3 = 3µA∂su, (4.2a)
m = 3µIdiag([1, 1, 2/3]) · ∂s(ω + uκ), (4.2b)
22
4. Conclusiones y trabajos futuros
(a) Numérico (b) Numérico
Figure 4.2: Electrospinning
con una fuerza aplicada al chorro de origen eléctrico τe = 2√
πAEσ. La carga se acumula en
la superficie y la convección es dominante, la conservación de la carga impone√
Auσ = cte.
4.1.2 Resultados experimentales y numéricos
En la figura 4.2, se representan dos imagenes del electrospinning, la de la izquierda es una
imagen de laboratorio en la que se ha resltado el fondo y el chorro G.Riboux [2011]. La de
la derecha es una imagen de una simulación numérica en la que se resuelven las ecuaciones
linealizadas de 4.1-4.2.
Una vez comprendido el origen de la inestabilidad lateral y desarrollado una herramienta
para calcularla numéricamente, se pretende resolver el problema eléctrico en vez de usar una
solución impuesta a priori. Una vez modelada la dinámica transversal y longitudinal del chorro,
así como el problema eléctrico, se pretende reproducir los resultados obtenidos en el laboratorio
del departamento y reportados en G.Riboux [2011].
4.2 Conclusiones
Hemos resarrollado una heramienta numérica para capturar la bifurcación de Hopf para un
conducto en voladizo. Posteriormente, se ha implementado un código para resolver la dinámica
tridimensional nolineal, considerando la geometría de forma exacta. Se ha implementado un
código para integrar en el tiempo las ecuaciones de movimiento del conducto, resolviendo el
problema de forma implícita debido a que no es posible hacerlo de manera explícita, y se han
calculado las órbitas periódicas.
Se han estudiado varios casos ya resuelto, comprobando así el código, y posteriormente se ha
resuelto casos con orientaciones de gravedad no consideradas anteriormente en la literatura. Se
observa como las fuerzas centrífugas desestabilizan y las de Coriolis estabilizan. Posteriormente
se estudia el efecto que tiene la gravedad y el flujo en la solución de equilibrio. Seguidamente,
se comprueba la importancia que tiene en la estabilidad del sistema la forma de equilibrio así
23
4. Conclusiones y trabajos futuros
como los esfuerzos axiles a los que queda sometida en dicha configuranción, comprobándose
que tanto los axiles de tracción como la curvatura estabilizan. Con ambos estudios, se puede
llegar a las siguientes conclusiones. El efecto de la gravedad depende de su orientación y tiene
dos efectos principales. La componente paralela al conducto, contribuye a los axiles de manera
que es estabilizante o desestabilizante según provoque tracción o compresión, respectivamente.
Mientras la componente transversal contribuye claramente, a la flexión del conducto que se ve
enderezada al aumentar el flujo (fuerza centrífuga). En este sentido la gravedad estabiliza así
como lo hace la fuerza centrífuga que desestabiliza el sistema en términos generales pero que
estabiliza según endereza el conducto. Se concluye que la orientación de la gravedad es muy
importante y condiciona la estabilidad del sistema.
Respecto a la solución numérica de las ecuaciones, se ha encontrado la difultad que tienen
los métodos usados para converger y que es necesario buscar una solución conocida o fácil de
calcular, y utilizar métodos de continuación.
24
Apéndice A
Problema de Valores de Contorno
Sea el sistema a resolver
y′ = f(x, y; µ), (A.1a)
g(y(a), y(b); µ) = 0, (A.1b)
con a ≤ x ≤ b.
Se han implementado dos métodos para resolver el problema de valores de contorno con
parmetros libres, el método de shooting y con diferencias finitas centradas.
A.1 Método de shooting
Este método convierte el BVP en un problema de valores iniciales que se resuelve haciendo uso
de un integrador para ODEs, entonces se integra f con condiciones iniciales y(a) = ya hasta
s = b donde y(b) = yb y se obtiene
yb = ya +∫ b
af(x, y; µ) dx (A.2)
derivando respecto de las condiciones iniciales
dyi(b)dyj(a)
= δij +∫ b
a
df
dy(x)dy(x)dy(a)
dx, (A.3a)
dyi(b)dµ
=∫ b
a(
df
dy(x)dy(x)
dµ+
df
dµ) dx. (A.3b)
Imponiendo las condiciones de contorno
g(ya, yb(ya, µ); µ) = 0, (A.4)
25
y resolviendo por el método de Newton-Raphson
[
dgdya
+ dgdyb
dyb
dya
dgdµ
+ dgdyb
dyb
dµ
]
∆ya
∆µ
= g(ya, yb; µ) (A.5)
y actualizando
yν+1a = yν
a + ∆ya
µν+1 = µν + ∆µ
A.2 Diferencias finitas centradas
En este método se discretiza la derivada y se resulve el problema algebraico resultante. Se usa
la siguiente abreviacón, yi = y(xi)
ǫi = [yi+1 − yi] − (xi+1 − xi)f(xi+1/2, yi+1/2; µ), (A.6a)
g(y1, yN ; µ) = 0, (A.6b)
resolviendo por Newton-Raphson
ǫν+1i = ǫν
i +∂ǫi
∂yj∆yj +
∂ǫi
∂p∆p = 0, (A.7a)
gν+1 = gν +∂g
∂y1∆y1 +
∂g
∂yN∆yN +
∂g
∂p∆p = 0 (A.7b)
M
∆y1
∆yi
∆yN
∆µ
=
y2 − y1 − h f1+1/2
yi+1 − yi − h fi+1/2
yN − yN−1 − h fN−1+1/2
g(y1, yN ; µ)
(A.8)
con
M =
−I − h2
∂f1+1/2
∂y+I − h
2
∂f1+1/2
∂y
. . . −∂f1+1/2
∂µ. . . −I − h
2
∂fi+1/2
∂y+I − h
2
∂fi+1/2
∂y
. . . −∂fi+1/2
∂µ. . . −I − h
2
∂fN−1+1/2
∂y+I − h
2
∂fN−1+1/2
∂y−∂fN−1+1/2
∂µ∂gy1
∂gyN
∂g∂µ
(A.9)
26
y actualizando
yν+1i = yν
a + ∆yi
µν+1 = µν + ∆µ
27
Bibliografía
D.H.Reneker A.L.Yarin, S.Koombhongse. Bending instability of electrically charged liquid jets
of polymer solution in electrospinning. Journal of Applied Physics, 87(9):4531–4547, 2000.
22
D.H.Reneker A.L.Yarin, S.Koombhongse. Bending instability in electrospinning of nanofibers.
Journal of Applied Physics, 89(5):3018–3026, 2001. 22
T.B. Benjamin. Dynamics of a system of articulated pipes conveying fluid. i: theory. Proceedings
of the Royal Society of London. Series A, Mathetmatical and Physical Science., 261 (1307):
457–486, 1961a. 1
T.B. Benjamin. Dynamics of a system of articulated pipes conveying fluid. ii: experiments.
Proceedings of the Royal Society of London. Series A, Mathetmatical and Physical Science.,
261 (1307):457–486, 1961b. 1
F.J.Higuera. Electrodispersion of a liquid of finite electrical conductivity in an immiscible
dielectric liquid. Journal of Fluids, 22. 22
J. Rousselet. G. Herrman. Dynamic behavior of continuous cantilevered pipes conveying fluid
near critical velocities. Journal of Applied Mechanics, 48:213–230, 1980. 1
M.P.Paidoussis G. X. Li, C. Semler. The nonlinear equations of motion of pipes conveying
fluid. Journal of Sound and Vibration, 169 (5):577–599, 1994. 1
I.G.Loscertales A.Barrero G.Riboux, A.G.Marín. Whipping instability characterization of an
electrified visco-capillary jet. Journal of Fluid Mechanics, 671:226–253, 2011. 2, 23
J.S.Jensen. Fluid transport due to nonlinear fluid-structure interaction. Journal of Fluids and
Structures, 11:327–344, 1997. 1
J.Valverde, J.L.Escalona, J.Dominguez, and A.R.Champneys. Stability and bifurcation analysis
of a spinning space tether. Journal of Nonlinear Science, 16(5):507–542, 2006. 4
N. Marheineke and R. Wegener. Asymptotic model for the dynamics of curved viscous fibers
with surface tension. Journal of Fluid Mechanics, 622:345–369, 2009. 22
G.Rutledge M:P.Brenner M.M.Hohman, M.Shin. Electrospinning and electrically forced jets. i.
stability theory. Physics of fluids, 13(8):2201–2220, 2001a. 22
28
BIBLIOGRAFÍA
G.Rutledge M:P.Brenner M.M.Hohman, M.Shin. Electrospinning and electrically forced jets.
ii. applications. Physics of fluids, 13(8):2221–2236, 2001b. 22
A. Goriely M.T.Mathew, A. beauregard. The nonlinear dynamics of elastic tubes conveying
fluid. International Journal of Solids and Structures, 47:161–168, 2010. 2
N.M.Ribe. Coiling of viscous jets. Proceedings of the Royal Society of London, A., 2005(5):
3223–3239, 2004. 22
D.Bonn N.M.Ribe, M.Habibi. Stability of liquid rope coiling. Physics of fluids, 18:084–102,
2006. 22
M.P. Paidoussis. Fluid-Structure Interactions: Slender Structures and Axial Flow., volume 1.
London: Academic Press, a. 2
M.P. Paidoussis. Fluid-Structure Interactions: Slender Structures and Axial Flow., volume 2.
London: Academic Press, b. 2
T.S.Lundgren P.R. Sethna and A.K. Bajaj. Stability boundaries for flow induced motions of
tubes with and inclined nozzle. Journal of Sound and Vibration, 64:553–571, 1979. 1
T.S.Lundgren P.R. Sethna and A.K. Bajaj. Hopf bifurcation phenomena in tubes carrying a
fluid. SIAM Journal of Applied Mathematics, 39:213–230, 1980. 1
W. Roth. Instabilitaet durchstroemte rohre. Ingenieur-Archiv, 33:236–263, 1964. 1
Rüdiger Seydel. Practical Bifurcation and Stability Analysis. 11
Steven H. Strogatz. Nonlinear dynamics and chaos : with applications to physics, biology,
chemistry, and engineering. 4, 9
Geoffrey Taylor. Disintegration of water droplets in an electric field. Proc. Roy. Soc. London.
Ser. A, 280(1382):383, 1964. 22
A.L.Yarin. V.M.Entov. The dynamics of thin liquid jets in air. Journal of Fluid Mecanics, 140:
91–111, 1984. 22
A.Meister R.Wegener W.Arne, N.Marheineke. Numerical analysis of cosserat rod and string
models for viscous jets in rotational spinning processes. Mathematical Models Methods Applied
Science, 20(10):1941–1965, 2010. 6, 7, 22
J.Schnebele R.Wegener W.Arne, N.Marheineke. Fluid-fiber-interactions in rotational spinning
process of glass wool production. Journal of Mathematics in Industry. 22
29
Top Related